From 8c403140229764d69da1b20054c72ea6d98fea0c Mon Sep 17 00:00:00 2001 From: Sammy Sharief Date: Mon, 28 Oct 2024 14:05:13 -0400 Subject: [PATCH 1/7] Update to PQMass I have done the following 1) Updated the README to include more detail about PQMass as well as examples (in distribution and out of distribution) as well as advanced options 2) Made PQMass GPU compatible that have the same results as the cpu (default) option 3) Renamed "whiten" to "z_score_norm" as whiten doesn't tell the user what they are working with 4) Updated notebooks to include device parameter to showcase that it could be 'cuda' if the user has access to GPU To Do Later 1) More encompassing test cases --- README.md | 143 +++++++++++++---- media/Voronoi.png | Bin 0 -> 24838 bytes notebooks/mnist.ipynb | 8 +- notebooks/test.ipynb | 14 +- notebooks/time_series.ipynb | 10 +- requirements.txt | 3 +- src/pqm/pqm.py | 298 +++++++++++++++++++++++++++--------- src/pqm/test_gaussian.py | 2 +- 8 files changed, 364 insertions(+), 114 deletions(-) create mode 100644 media/Voronoi.png diff --git a/README.md b/README.md index 20dd82f..06c7503 100644 --- a/README.md +++ b/README.md @@ -6,63 +6,152 @@ ![PyPI - Downloads](https://img.shields.io/pypi/dm/pqm) [![arXiv](https://img.shields.io/badge/arXiv-2402.04355-b31b1b.svg)](https://arxiv.org/abs/2402.04355) -Implementation of the PQMass two sample test from Lemos et al. 2024 [here](https://arxiv.org/abs/2402.04355) + + +[PQMass](https://arxiv.org/abs/2402.04355) is a new sample-based method for evaluating the quality of generative models as well assessing distribution shifts. ## Install -Just do: +To install PQMass, run the following: -``` +```python pip install pqm ``` ## Usage -This is the main use case: +PQMass takes in $x$ and $y$ two datasets and determines if they come from the same underlying distribution. For instance, in the case of generative models, $x$ represents the samples generated by your model, while $y$ corresponds to the real data or test set. + +![Headline plot showing an example tessellation for PQMass](media/Voronoi.png "") +PQMass partitions the space by taking reference points from $x$ and $y$ and creating Voronoi tessellations around the reference points. On the left is an example of one such region, which we note follows a Binomial Distribution; the samples are either inside or outside the region. On the right is the entire space partitioned, allowing us to see that this is a multinomial distribution, a given sample can be in region P or any other region. This is crucial as it allows for two metrics to be defined that can be used to determine if $x$ and $y$ come from the same underlying distribution. The first is the $\chi_{PQM}^2$ +$$\chi_{PQM}^2 \equiv \sum_{i = 1}^{n_R} \left[ \frac{(k({\bf x}, R_i) - \hat{N}_{x, i})^2}{\hat{N}_{x, i}} + \frac{(k({\bf y}, R_i) - \hat{N}_{y, i})^2}{\hat{N}_{y, i}} \right]$$ + +and the second is the $\text{p-value}(\chi_{PQM}^2)$ +$$\text{p-value}(\chi_{PQM}^2) \equiv \int_{-\infty}^{\chi^2_{\rm {PQM}}} \chi^2_{n_R - 1}(z) dz$$ + +For $\chi_{PQM}^2$ metric, given your two sets of samples, if they come from the same +distribution, the histogram of your $\chi_{PQM}^2$ values should follow the $\chi^2$ +distribution. The degrees of freedom (DoF) will equal `DoF = num_refs - 1` The +peak of this distribution will be at `DoF - 2`, the mean will equal `DoF`, and +the standard deviation will be `sqrt(2 * DoF)`. If your $\chi_{PQM}^2$ values are too +high (`chi^2 / DoF > 1`), it suggests that the samples are out of distribution. +Conversely, if the values are too low (`chi^2 / DoF < 1`), it indicates +potential duplication of samples between `x` and `y` (i.e. +memorization for generative models). + +If your two samples are drawn from the same distribution, then the $\text{p-value}(\chi_{PQM}^2)$ +should be drawn from the random $\mathcal{U}(0,1)$ distribution. This means that if +you get a very small value (i.e., 1e-6), then you have failed the null +hypothesis test, and the two samples are not drawn from the same distribution. +If you get values approximately equal to 1 every time then that suggests +potential duplication of samples between `x` and `y`. + +PQMass can work for any two datasets as it measures the distribution shift between the $x$ and $y$, which we show below. + +## Example + +We are using 100 regions. Thus, the DoF is 99, our expected $\chi^2$ peak of the distribution is 97, the median is 99, and the standard deviation should be 14.07. With this in mind, we set up our example. For the p-value, we expect to be between 0 and 1 and a significantly small p-value (e.g., $< 0.05$ or $< 0.01$) would mean we reject the null hypothesis and thus $x$ and $y$ do not come from the same distribution. + +Our expected p-value should be around 0.5 to pass the null hypothesis test; any significant deviation away from this would indicate failure of the null hypothesis test. + +Given two distributions, $x$ and $y$, sampling from a $\mathcal{N}(0, 1)$ in 10 dimensions, the goal is to determine if they come from the same underlying distribution. This is considered the null test as we know they come from the same distribution, but we show how one would use PQMass to determine this. ```python from pqm import pqm_pvalue, pqm_chi2 import numpy as np -x_sample = np.random.normal(size = (500, 10)) -y_sample = np.random.normal(size = (400, 10)) +p = np.random.normal(size = (500, 10)) +q = np.random.normal(size = (400, 10)) + +# To get chi^2 from PQMass +chi2_stat = pqm_chi2(p, q, re_tessellation = 1000) +print(np.mean(chi2_stat), np.std(chi2_stat)) # 98.51, 11.334 # To get pvalues from PQMass -pvalues = pqm_pvalue(x_sample, y_sample, num_refs = 100, re_tessellation = 50) -print(np.mean(pvalues), np.std(pvalues)) +pvalues = pqm_pvalue(p, q, re_tessellation = 1000) +print(np.mean(pvalues), np.std(pvalues)) # 0.50, 0.26 +``` + +We see that both $\chi_{PQM}^2$ and $\text{p-value}(\chi_{PQM}^2)$ follow the expected $\chi^2$ indicatiing that both $x$ and $y$ come from the same underlying distribution. + +Another such example in which we do $\textit{not}$ expect $x$ and $y$ to come from the same distribution is if $x$ is again sampled from a $\mathcal{N}(0, 1)$ in 10 dimensions whereas $y$ is sampled from a $\mathcal{U}(0, 1)$ in 10 dimensions. + +```python +from pqm import pqm_pvalue, pqm_chi2 +import numpy as np + +p = np.random.normal(size = (500, 10)) +q = np.random.uniform(size = (400, 10)) # To get chi^2 from PQMass -chi2_stat = pqm_chi2(x_sample, y_sample, num_refs = 100, re_tessellation = 50) -print(np.mean(chi2_stat), np.std(chi2_stat)) +chi2_stat = pqm_chi2(p, q, re_tessellation = 1000) +print(np.mean(chi2_stat), np.std(chi2_stat)) # 577.29, 25.74 + +# To get pvalues from PQMass +pvalues = pqm_pvalue(p, q, re_tessellation = 1000) +print(np.mean(pvalues), np.std(pvalues)) # 3.53e-56, 8.436e-55 ``` -If your two samples are drawn from the same distribution, then the p-value -should be drawn from the random uniform(0,1) distribution. This means that if -you get a very small value (i.e., 1e-6), then you have failed the null -hypothesis test, and the two samples are not drawn from the same distribution. -If you get values approximately equal to 1 every time then that suggests -potential duplication of samples between `x_samples` and `y_samples`. +Here it is clear that both $\chi_{PQM}^2$ and $\text{p-value}(\chi_{PQM}^2)$ are not close to their expected results, thus showing that $x$ and $y$ do $\textbf{not}$ come from the same underlying distribution. -For the chi^2 metric, given your two sets of samples, if they come from the same -distribution, the histogram of your chi^2 values should follow the chi^2 -distribution. The degrees of freedom (DoF) will equal `DoF = num_refs - 1` The -peak of this distribution will be at `DoF - 2`, the mean will equal `DoF`, and -the standard deviation will be `sqrt(2 * DoF)`. If your chi^2 values are too -high (`chi^2 / DoF > 1`), it suggests that the samples are out of distribution. -Conversely, if the values are too low (`chi^2 / DoF < 1`), it indicates -potential duplication of samples between `x_samples` and `y_samples` (i.e. -memorization for generative models). +Thus, PQMass can be used to identify if any two distributions come from the same underlying distributions if enough samples are given. We encourage users to look through the paper to see the varying experiments and use cases for PQMass! + +## Advanced Usage + +Depending on the data you are working with we show other uses of the parameters for PQMass. + +### Z-Score Normalization +If you determine that you need to normalize $x$ and $y$, there is a z-score normalization function built into PQMass, and one can call it by setting `z_score_norm = True`: + +```python +chi2_stat = pqm_chi2(p, q, re_tessellation = 1000, z_score_norm = True) +pvalues = pqm_pvalue(p, q, re_tessellation = 1000, z_score_norm = True) +``` + +### Modification to how references points are selected + +The default setup for selecting reference points is to take the number of regions and then sample from $x$ and $y$ proportional to each length, respectively. However, if, for your case, you want to only sample the reference points from $x$ by setting `x_frac = 1.0`: + +```python +chi2_stat = pqm_chi2(p, q, re_tessellation = 1000, x_frac = 1.0) +pvalues = pqm_pvalue(p, q, re_tessellation = 1000, x_frac = 1.0) +``` + +Alternatively, you can sample the reference points only from $y$ by setting `x_frac = 0`: + +```python +chi2_stat = pqm_chi2(p, q, re_tessellation = 1000, x_frac = 0) +pvalues = pqm_pvalue(p, q, re_tessellation = 1000, x_frac = 0) +``` + +Similary you can sample reference points equally from both $x$ and $y$ by setting `x_frac = 0.5`: + +```python +chi2_stat = pqm_chi2(p, q, re_tessellation = 1000, x_frac = 0.5) +pvalues = pqm_pvalue(p, q, re_tessellation = 1000, x_frac = 0.5) +``` + +Lastly one could not sample reference points from either $x$ or $y$ but instead sample from a Guassian by using the `guass_frac = 1.0`: + +```python +chi2_stat = pqm_chi2(p, q, re_tessellation = 1000, guass_frac = 1.0) +pvalues = pqm_pvalue(p, q, re_tessellation = 1000, guass_frac = 1.0) +``` + +### GPU Compatibility + +PQMass now works on both CPU and GPU. All that is needed is to pass what device you are on via `device = 'cuda'` or `device = 'cpu'` ## Developing If you're a developer then: -``` +```python git clone git@github.com:Ciela-Institute/PQM.git cd PQM git checkout -b my-new-branch pip install -e . ``` -But make an issue first so we can discuss implementation ideas. \ No newline at end of file +But make an issue first so we can discuss implementation ideas. diff --git a/media/Voronoi.png b/media/Voronoi.png new file mode 100644 index 0000000000000000000000000000000000000000..978679a300e5309b25032aaa2851a7763a1aff3c GIT binary patch literal 24838 zcmbSy^Lu1n&~0pMVjC0Nw(X9MiJeI@v2EM7ZJRT(ZQp*sd%u6+{&3E-`<$+Rp4xkN z?W(m_g)1pYBEsRqfq;M@N=u2UfPjG30mo}FP{4acPD?rP0hN`gsFIDSq^N_fgOjSG zk%^hGse`eZilhiLH#Zju2z9&(0H7j8Pd#RY1^|pr)6u{=xvPYQMX3P%`}>B82m8kI zrt(tK^z{B=qx>5LfndVL^y!3BgM~@WR;M~H3NTQe^3Is(tCI#^ZcW+(;H7kFiSU)PZN;zDP=eQRS zoCK5vB?3a8V`$kIJ`9qgL_I4EGlWpsbUXlXV>r>*cLz;kzA)KWClneG;7=0pZUJhyu8b{iFwPBXH)p%(SG<<>f(WfnyjD&`>K72;c}5 zc;Nys5D>6Du>XAmT9*g@zhjWK|6Vo+_mY5s2!Tk839EU4Uiw0Nplhx?#oGqiF%c}< zEqY(0z;u(Nz)%!HFAztQXHK_{Q~tGgw3@23zkAr>Z40giRc|7Q?~AlC!_|I7iv?}M`*m?TE0{EsQx zf2O?uf2IK*%_w}I4~s5`6PUEREma$`rrhB|Q^a*PkvNSdqbK=;*G|wJXts-s(L*o^LPvvFhs8 z836&Fh=hD;860*M6~@|9q$sj5d@uikDmuL{H`u-I#;Md(S55jA1wWG{V(_KaE7ayq z?jxE{BN==0goB`zFE^Sa30Ci8gP@UBh9fY`oen3w>`RX9HXALF8MNwzT)Qi9u*n&j zYM(AQD{J(yl3<45z~Us9&?azw@tds{Ln<_D6{>YxH8wB9l6l=Ou!o0-ap0W)MRJUi zR_nE^xV9eRrdgMn3ja6L!o=+C$un8pg&#|gSm+K@7ttc2$V4TrJ6@C2)6tf^6yTB6 zA9qJ+peCeB)jI4v-cMB|xPSA6!Fl{XxDF2waSFOS-LFf`r!#O0s^jg8tTBy?@My0c z=ppJn{LYqW>QU^fQoJ4PL`1-L_w6`FpSFGO{!B2jjv{x3`I32m+Kj z1}p-jd?Ka%A6(^YgsQr_I9yJ9`Y+PA2E2o#gaDiaVVvr8CjA=2UV-GF`Wu$=&l?8V ztVSUS=*L?6rqDulaf7Tgv$Lft#WJ45a>!KmNKez8e^2sh3^q}Sh!h7y5bt#)a=L|7horwa}z zQ*CsuFjgsXv@itq*4e{_s8N{`NT9KQsyTbLN^Es`HqPn&xoa0!-A)J)UtCP`{T11l z#qDf45`|+w+s+HzOo9vk-|vR*r3cJZ5K*QwvF3vfin9n81Hte~iHRoEp*E$fB&3}{ zYX~>bMMsi4k6YvX1RQqjq%&P4B7v-xTElBiF_7hCNoC(;ODij9fnNM&WjBpApQ=+< zK|_V!#bUG9%cuP0aye0A&aFlkc==7}d4L)ygW9ieF`L716{P*TP%KOEEEcx`4voa+ zH<_IPLu}18b>i&oESs}*irO8L%;WVyO@sR73YYyelL9C0Cz0$0dxeD}IwRok*)v%z z#J9EI)4B!eaJch7y&^~vHV-lLllemE?(eVny^%4pW4WPfBOI)BQ)nV1;_+(DkZ;cXe9uMuB+bNeu5*-mm=;n$(GE zy#HD8qtc?HFbs+~j|Z_ZoyB?+Nq+Wm6uWqYDU?_m4qVUoOLLu0j)1Rc{~!wkIAbMr zkv%YyiRaJfzl`(<4m(}a-!G$rWy*yTf#W1%YSlW;Ns}-{Z5iw~Jbse|&oJ2XV4{q~ zr^NVaIVU~uSA|hg2>lh4O%F&E~Mf7oYt@4$-q6g#lvSAzwwO(L7qp zr);>}xEP6viT_=r-5Uu2kjj98)kXf}W477qkP7Mp?t!fW1`opK8U72UuL?^*U*?2r zhy_A8x#&-aNRyKmIjIdJObD#aY7OGphk+W3=y#N5_*$Kj=)}atg1i@ImyJXiRr%Bi# zK=X#jw36?6bQ=;BCMmrlCcOG%_i_vQN-yA+mWBzUATS`{7Xhz%p=hY>hsqoC9oV3p z0R?K94#ZmL?`_1e82ozff`XT`D%90}Kp*&z2&s(4!8*&ghf`ASGRCMw&u17hNO17Q ziSc@6r-o>7*Nrt{br$B;Q9F~#j==ppGx8q=72ho2cVz%^-eoAA_~3(#IMd);SCFKbr9r!L}74H8r2!$wCXId3Hl zD8=aNv?@nd9AAEB!Wcj#c zgTT5OS=%FXp)bS9trL=>l6Ie2XFLoSs;K4W<|ei(l(ku$`ngankAJczhJynzKEu{R zD@FZH^Ris2`M7|%e5q6^A=dYKfex9J<{bWWOR9*@bKKNO-8fGLo3xjHcRXV!5UvtR zfeNxPFyItIR-a=qAye@PPm{Vwu0iXYg@W=-1+Me$vemTPPagndUXqM)l%;UERIW-& zxR>eBLvwZidYtQgH?_S-#$va@HhLdq)!2mK+?Q`X8dmaff{mSqD~d$TV5`Z9J_}AG z&@XoM@~@~)zte4grCOI@9-tOL51%X~#8QVAC$<`J7XD+*I-)et2X7gW9LJple&}=| z$X^{yv+!f(1cDiSQ3TFez7$My6KLFMUDeSxeRnun4h+HrNgL~)I2FZimm8=mwOU&0 z6di`~kl;fGP=g9dOG~CzH|ESTG2pe@mo^-xHJ7Q) zQDmoU6kujdt*W{j(HiJHW0VF*@<(+xGzdv>ij@_MIXYbDimM72y+f(Q3Umkz7wmz( zs4I8x7{c;@5lTx)B#?)^?`C$U)Lz@5tFcxe{pIhTwOIa7hn!7J?CtNz5$967vjr+T z@=R5ofZ7RPt)0}OIV7U1u2eNJWmMDTLMX1$Z&cRH3fH&s`75UlJvlaZqlLYoQ4WR3 zm-BoKvC?heu?rbkp`+pmV|kN~Py>2r%?rlqu;Jl5LiynsL)d(_ti(hfOBcNg&Vocx z@ZkH7?~h|EW0)%c9b{y1h1=B_H=9SQ=kW~jO zX&Y5~f9$KfMTn6voESYte#ryHx?oS*=CuD#h}N5;kBzPmL0;uz5Lc2;S8G*yRWvsfZO_Z(9U zhXHkD!$7l?*2fONaZ9T%N56}rEAiGtwout;32o&vAPh{9NVWQ8;0Qs>MbrZFn|~W} zWa9C+(geues(W$}eeg)PJCP3nj1ZD(A*ocR@R!q>D^EnGu!MOqt4xILxJ<`T-mW^Y zFILjJ-16g`X)(#0#j^PHsHFM9ie*#HUG-T=I=sXYDGAT6^}5K_uEWDQ>_hq6gsXL0 zbDVDzTK)1bt*op_yJhZ=rX8|j`U$c5R(i9dLqr@i)5=v?n22IG%-e;zL%F&`qK8SL zg&l|ceh2ID>5{RRs;E;b>7^^gpnY8dwkz%RJ0H}0_;lfM^*G!((L-riu-QDI3;4Kh zQT%uBD!XiWFdgD`I$Y0>8twJ2Yy1Tmqwc^WM4icTu!$99SsTr0o@_-51NaSrM^R`- z0Ku`NJrh1G7TKTKAQSexM$$OZs)T;Q#C(6-k~S8#l5yOK+K1n5+Is0VNn57bFR~T9 z(K71Z@m1}eG@R!+4SI$@)PfQc!+fN@N;o$fJBheFUqO?^#jeZ#NsA^0Y) zgj3{Z7D&*0Qg%r)F$YS3pLsN6)VS6;ZPkUNQ&rfdY&wO#XJp%5-K$Yx$HfrCeNnc`6>?drdwhfp!Mq+WCPoo_hKY|Rvv)gZJ(MQpp zag-q|=C(i!RU2&ENr!-fva_Ljgi4iZ9SyKTnDUY;K50-qakF=%PCu$ zmPQQO)>_2wyvm`K19}dRD3vGlO4=UGrg=r$M{@v2DT=q4A@@A%*qF`3Zm5mgEVY)t z;KTv_r88xq*1HEWD2n-{)rAg24QEQ85~%TgagJ#AC}~-u7$cJMv!?CwG#M6?f(ow| zI}|n~$O~KiQ>`%xH2ECKHJ^$RitYNd`f+xVDi95g)mc-Rev0NsZ7AdR0|cNFa9Yn* zV7Gk$wL%F*)OR!+BlV!Luf6ik=?r$Vs%Gg0q<}hEKv;A}@^Apt>kJb(|1XiwnLh=T znaL(OhZiX7qp5gOlo})%=y*v7=!sIA4S0|s~;fte_m9RB;H)B-|49V?T^1er4n+L17(Ek`4Fbc}SBzjD6T zw@ldw3WXZBL@}+PNcgFt?w+CkWKp?YRcTHi~NzG%tC?^p+K|?EFyQ*)%@F+JPb4xt-Pv0M-+HEs= zAeh`17N<#Rn8=y=QLweB5lKXrCB??lSj-|EmHe-}pSg1O5c{X%%y9jsI-LVodFzlV zLm7NV0A}Yg{q2O{lzU0oC-rF zD+8wC1Ji{-s6k*8HXKn-&d73Rb1)Q{HztUCz(pVlX+Z!FuK;Wgq>q78Wx;<4TPy5l zsYoYtdZ6BW_U>(ER*wk!_Asv-zd3a|xZSp^|5Pho?qgH93roXsSC_VnPnU%Wn)+=^ zxahl3vi=GVlOy~kaRFbu(*p8bREuOF$)MUxzO(aPQ`OwOL@U*Ly`fmDe!n&VQeaDo zX0H|qFBmuoBsy=BocI5MMej;sU67tqMVWweLkRma%x3Kg%^a{FLcKzQauS?;TuyR! zgppPIHdlZQG28pyGg$s_PfgPGs7S+HAXSSh0uw~ltPM6L5B@|~7n zAJ=8liEM9!L|=Wk;M|*LY2_SLb_7fS5cv$Ya09PX3c}5+UOFRtjlR}OeY}KJX>8PF zMj}UcVSQ>GEmR+^IQ@@vMVI>#_6dP(8T6cvY>Q(ouiaoSoi<0J>7>7 z8T0kF+KZN8@y7Bs_kJJa@=gxEPX!WLZkL&=oNMn|Zk_D(Fb1x9a$kGd$QDTvL+DE# z+u{0jRc2eSm0qu%c8YmNl)3cd=pe;kWiFmd;ONiIxw@vBtV%eY+Jm_J* zHOR4I;1S!Cn3*|hLK6D-Q*UB=SI{E32`KUm>+N=oT}KIdN2OwjbtoP5!IX zaTm9E)lRs!>y9(VCT+q`*Rng>`B(X=xaB0^(`Knm`Qb3taG|ocXTn+Yl`82Y7=BpUsXx-bRmvrv*m)`Tb-L6R!1lvO22L7bX$^3PZL>v zSs3-w(zJ9wxM^Vl#}B8C4>R7Z9mzn%gN`68JK)w? z&}jl7h=*ITr{rn7-W0XbL~qZ`UiX(fS&o5U#kXpL^SRmU-wad3zha6`)E zFvH1KXNV^s%@larpC{b97-MHceST5{sDNE>4x0icV&mcxm8Tu|hcM0wuv>2m)M+g{ z6Lcc!D$qY?<7x6|FU{QbVND-b1l(9 zLwNnSLY*#tKwP5Y876)eUAe!{c5EIEgoLY<;u@lOb^Zr9-SoMooo-TLvrP)b9Gz#l z7_krpiUFuU=v)`Kp=z%Ig9m{Ym2z<8+g!Jd47;5+rr!p^#|iU=l8ZgQ@0P$) zXSdx5V|(}pVQNIV5FQEt-^Cby%YfoWgB518eF^e%zRhALM}@d1LTXloSF;pP@2wfbXvxgO0MZ-d4??m3AQh8m}?n5qdns;@nX}*0+?+n5H~is==++ z?^a5CeRt-h05BAULvZCk9lJHyVh{tS$4T&N2elk_Fm(MpngCavOPv;~Pk#d;Yre{) z_a)RCb*$bVkH`Hukv!pJ@$2hx`~V5%8+bs(IJd+7TH&wXd-Z?Bh(zBU(CMJ)?yz-h z{&4!N+EF4J+xrrBdSCT_kQSe5D;+)YR?U@4=_<#ToV?nV6sCvM6Aj2MvZ-;_0!`qG zcAci&zGQ%Q$ABPJ_*hR`K7RvO6UHoF;my|SLIeKJvK~q_NJ*q?S;5&{h1D7ZC0}m{ zznPg9i|!FFxaREPxqQ)e|4XoHgKiyHyuxZ?A18SDb!iQX(*@+I5F|&*q2G5)sN!#X z!8-kaB8uv`FhRv`6S8=_6W`pdj>_G(@vpafIHnC7!f-jP>7C$St<#nx@`|5M2%eI# z-puBD-biPmV}B#=aq&)w$nJ*!c^FOx0~c#+%(DfOY`?5q68?u1BK75Pp^Y2n5u&4Z7@f!|jC|56B7v`${fDPb2i5(N zXm0WF9`wTjNSfSK`=+;`H;nX3boDOx3D!9Z@WwhXDu6}smLb`yGv+a1pj+9o1_;aoI7sW6*n z@hjCA{fT_tWUycmKyzTQ@cCuTV?I!>4^2+udXhTp$f0QH1A=C~3M(~d=HFapPor0f zo;2h*kFW5zBEMda$j`PN7?%}n_8Rq*=`y*s7+8LwP!)v`@bHh zg0G-(vnGSBsBnL;SQ90ntz4tLO>X9*Q7Y0TPNCu^`I<#yaQ^juCICXU2`Ti-4qYX@ zP$4*E2&_VBuo>IKenvBfV{%CFd>W-Jy-6G8zJlR&Z*Nmnzm_1E(E}6?h91 zjcqO_zd`If_4s8KD<(e_J>lM zkK}_(yuS&Es$VFoUOKcCBiSYV-fp2U3au5PXa7!=JBg;|3&-N!&a-l!wZI37K+xcC zup$}g;cvY92L}xC(N{vX_w%AC>J1Q}O8UJ!j8C7J@(jc)Y6kW+5h=N4TM`ZtTDi+H z^gms4GGqgPvKo6hWPC6BpT4pGO{Wa((N>1+;$xh;>mSnF*1Zun>~yQvg7-QrFMaW4 zC|tG+s)4TLlB)1NO^B|R_k5KA?I5gl6wk`e@qgJUnc{!&Tn zvRaVM+oHP7Lk~uUsm;(qcQ9IizL;{q5y%fx`ly4nomL+)5!1dMiwi~Mlu4y%z1zhG zs-p|z(L>Gf_GgP0oaJud>BExA=XN&I$=dL0&r$u*(U+X&*CTjxrXyAQ^-nhO`PbL-@!MndFNQ^|VB?(e)l*Z!sNAQ>ze)o$luySs?gW zagHD<+=O#Jme!?FWvMw1?ff^HPIKXCI+I&-lWqzj$V9j!z7Hf!Yjv-j6)QITJMRl7 z1d!ZlGYV(Dte4qLCjQp2U5CtrcYZZU5v*RJm2##wdptaY*e_J@rZd%ZGd+kvFR(K+ zre=u08_1L(!ZN32YyjqDt~RqjS|3}h6xZZ-h+h$(F+{Y^BfMhTi@%_0gcV5%ck9%@ zJUc)ZjM(1z!@^ff1|}F2sX}FV3T;ijr9oG|(cD+A_kDYEU=Io5ym{3C{pg-)ZH}xz z8YGH-Op*CN^@90qjyyV*!jjm;5+*qOa?obCH`TICvr;>TT^2)jF=O#`1V+h+M(6Z6 zv+A>7>wKYNvr#0sYf3!9hWT8N{v>zTw0d(ZY+DxTxT$GcxdA`DN#C+?nap+@ zNSU8D&VKdN4hY`0ZHxN8xfJw%UOEi8bh%&e3}8dWbTOI5YKYvpnCAPvot6SOHwWla zP~u z)D{|f-iLhE_3)(jVm@EqW+Y6#jPu59>G0w;RNi#gqe;tueVslj$SojgDaObsw)%Ho zz9b-`Kd?^#Pyc|DZbWBQGma8HDc95G7umtH^hKaD}vJq$kg zzTT$bmTYe93iR@%LAgm}{}_>odT~dH>{qrc=$b(evA*ZKq1Sk%!IX zien@ne7Q<{@v-TH0re&vI9oI3fCoQBHlhZ|avhE_W5L{P&fp^-cjV zYfqm*xK@b}HaQc0l|}L~2m(MSmj4MmGrKYoCSREbx^bPF414ZBHh-zt-`|L(*6bO~ zTSupb_;*JgcEs?|i$|vsC{#gbv(OOOR)ZDJ<5fZQLopj>)P$jVvEF%?yh4j3l&%Kj zN5({`Jh0mBv51;@{KaKv(7wyN7Mvh^# ziwOj=(nnbWvZz@H1SJkn31j)Q#CtuPAt@+ys6(kjzTYSO6vbunC391%)?_BJ6pIDG z3F&|BEP_hO0aZIF30QZYK#W}Je5LwyhWGj)T*yo@)sv4f8V$RrKiB6jJdW@E@%*^{ zIv^(u8+@<|xko-gX}QAO?W7qt5>SS+&CIyueKJQPBi$uZuHI}02xwYgWEb@94e(1d zarhh7*}#wK?9E5UpX_$EJ)I)wT`rwiVm~<+G#?FOQwbnp)PwrCq>0j1 z=hd1GOlo#L*I4n(DgdH*%I{B?Ru^maq{bnm@V;4np#q^6C|)eN&-140rSXJN^CDO6 zRBWiaqyv2l*$tIn2?@6N>I02$v*;o)_)4`o{1WM{&FRzWlV!4aQk-v;{y3@54jnFK z^1XE!8ohqtsNYksS+%}CFaKh}M#YA_9-E=|Dg`_Rx?qF?Cng&1#Ad{nD>K(X`o7;w zi_`xcK(4|lqW6D#b1p1TmcEGkhJJIx>|?4|9xM(%=fa%Cu-hyZxm8=wWf2MFdlPvL zD?_g|#^7`B0o~$}i*G<-Hed#yooziv1P@`i_83`EuAf9Is%`i%TI7L9nzS2ubN>2x zb-6o`_Bt(2^$3Jmu6e?6gM#YC_i!WoM?XMFNu{g3_0;R!LmwDR$S-s25@NkPr*WA7 zNd{9hLYq9Pe7VHhWadoidbYJ(JHd1>F+;{wlSpM}Yi}>2^j$^;uPdsZZbL=flLV9{ z;8YosbRpNX0rsV%zKIlkn0&SEXLH}#5l=QV1Z%!xRmu8!U8IP-_NN`nB(G6jc8zNQ zI$E;>9T;g?4V0a_vnrH!L;A3r$v0~?G^YQkXq`Q%?sSF$Vf+b%3Ng7VnU ze^%7bYkC@NBe?Hs6dpbr6_r|Z2JsV8G9KxpvfH>i6C0+Fi0;8qB=o);jG?9o7)_;% zlF%BX-y8hD{W1i6JL0-T?2^!c^POy1tz+Kdff|AxR$u(-jZn}2eZ(xIXj9{R$W!%e zwkJ`oFvWige!EO(nS2kAlon(qtGnF9b?x$ZhN*ZO#S3Izjif2>oc0p9KM&-iw_8j2 zbW?uB^OeG$V)u;M@ocu3uXn3s*+kZ?_@}OO;_las@i5l>c~R))S_y5K z4r;wtcZ%I?yQoMH0==;%1f`ev9C6?p@emvUEpoyVe@UDWy>Q!oJew}o>h6HHz>d)E z@(t%7=daxode#0t{5$i{VZ0DSC1Nv_2DvZII#L`pdoBa@-*5kKAeMM8q2J*Gi!JKU z7{vqCVs?*fyEVsg9|{ zrRMif=6V6t3*I1G=v%FFS&*<$(oZ;s=DLJ>HX!2XsYDT0hBFIwY(lrD`$8^&=_A78(EG2E}@H8 z=;OvRe3RQW+Sz7l=1%5Wqp$Okdf`$>P?=7*Fq@72xEZ7zMCr7i1et_jf5pO>c994e zeDiRQV%-uJAUU5UB&DJKnT(3I8<{FST(^~Duvn zI>u!%4Lc}E2o(nX@uG3| z?a;e|y{5G%L1z_z+gN1`yqg4fIS@Jyb^ zV@e{s)5RSh-i3-YH5-K+Kb}1NcVJ(GVTnp4<#L&f*?wm}vF>opqinf^t9Q~K;-n_^ zHpV6-e5y0V+l*SP-?{PV=V!p9^-S8Hy{EdmI*ej41|4iznJoygqo0P*9Vh8K^d+m_GK z7M~EBq)236_Mip*F7GGm*z9bv;g9Kp<6PKU3T^6s%Y_pjv&+}!A%>Y8Ud754q7uRG z$78hP5ANTV)KbzxWgODdl!UT%CN~cBIMTi_M!UcorZ8VLR3#FN;Uj>S+zgR1fJA$1 zC%dxD>>da=b@3nWQ_U_eXH9(UHQi}U9djdahy@7aA`S^2fNTPh-6rdGgb69d{0yV8 zJQS-o|Fq^bC`3Z}L%JfU<)q2f{NcJ$bn%Q*JkF2W8S1k^5sg)$JL}CBhR2Kb#qXGq ze{h4Wb=kbFs!*XzH2X;e+z?#iT{1&4dcCIzuGhby-`gOc*Wx4*&%uHE8Lsa3=^7q z4WD~CPh_YKj=aT&2@QKG73KiJE&>lZKJUF$4a8fNrC;U=4bzR4-=6?hX{ogDsQ(4|)V#^ijkVo43ao~9hl3B!KtC{A+OhHJ(o1i9l~sFRwImHqv2k6s>Wpt5Z&&vly@{WX zB*`(-*zp(`e+#GTYIP(rX-kzepM2TRb%4sD9IXs+ zq!2_XtT1%emA`fPzki?Hcqlc4^ynOCo^+HbqEMJa;3KdNs)%9J7CvVbO2MJPX_ylV>XL?I@8z?F~gqQ#-vH5g<*1mX~iRLbzXUh{yF% zHWW4xHj3k|!2SB$RaY z_M^btIUmn{rd(mq2{due3tu|>>r6-+zN<`cycsXW_2F{ZVxx$R4Mlr0IK{up+?`(| zD%9nC0|_Ox1Wj7=Y#wq^832&0+3rfxX{E-&l?cU#UM_M?rz$f^4ft?^s;xz?Ai*_6 zNc}P1#RkRx`^@c9OB&KMwM>2TWG37F&XmJ$-7n$__tgc2I`}C;NIguZ63|ooc)SbL z;yv_36U~7m^YtZ8drefkw6c?f^SYd%P*CZWoRwp`K$oQ;^7tIEX-~P9{o1G0(V#m| zfZ;+(6SQCfvTG8W@b@tzR&A_yBegE8Vf9IR68_dN0GzZCrtpe?T8M_VK;^7*r`u&y zjltp@!gybct~uD8 ze|Bg$TUaBli(CP$?i}bzuN+h zS34WlL(Y&{U-`eOMC#?J8lD=EcL>i`C_3bi{a;-ni&t%=qchJ~nW>gS&A1W6dv{Ak zZm(u?1#@GF1QlWhzuG#a5jSs-iTI>ega%?Qmq10C>r^kwQ{DHmLp$e7#iJx<|B%_^ z7h>7HotbNtsP{f7A+xaRdXCUFFgMl+YFK@3NiWq}utq2{RgpBzO1g(5q5`Wy$>*)- zz}bd)+Pu`qS|As@5P?Z|(tey{(u(D}vle`@CPy!5VhoV_wNqgQsXwJ5 z@>tg}t5P>0qB78qpF_GYK{p;gkt3(@}QVp?NEg&ah? z;<**&VqF$(@|3C;N0^a4Ed-y9ydz~53?+?ZyzBo9GT#%_unny^4Jp%%?@1b=RnhWc z_#9psVjbKQXwmv7=yVVkCJtQ!aRzJzAR+rH@t9#RY;ntp=)cdaVB8b7y`Bj=jHb1A1ud_ zH7F~Ka{Y6cmbCd?Gi73ACc)k~j>|R?(Q7~AxXEqGTl@;9dACb?nXZa zxKsEu(W~BN>20G!CV^kcbfM~SLAi)9a#oLL%Scdv@MY7=x^GAD0zFBN@qB&(TS7$W zG;1(M`+g^t!N4%Al*r@tC-ohm(BS(^dG2P^gTUXQ7Gb)YH0Nage83Zy99HLJ9S87zE{Qvy0g7iTpKT&?W&Gb9<-A;oul04cnsz=NG3a{{S3{2o<{pRV9XX3*fO?5oy*?&q* z-!aR|ITn(lE?A|dcR(s@rSEw%{A6I5?*DCFPYy$K3heX27|mK{Mi%^V40ne{MFWdY zP^$64^!g}owYXr;X7T1#p_H>9 z>zTi`#rO8hfRg`v)4e+h?ji}1$W~)o$5YO#*ee7B!(J1o>ps{2iyMfDsmjHCh?C_| zNprY|5b}*Q#mg?YjyWPA=nqacxnF~3DOO+tlhxEG0%PZY3OGqB8YZp*Uigv8K2OCd zP2T|RoT9Xooq7TqXGmUy)-Ii@sQcJSP5ZfBbV`NWJl~&`X>B^5m!nf99a`VanpN1D zD%h1ibME|cd-p+0TvrSG)W?JKE0oh64H!>EIxmalUgwwPOTCKDXRF5?`c?e~y@XV0 zKOD7M>q*t8vE?9BkGz7-zQB-}S6g_CSC32PiSG|agQ-a-Z>xbnupMbd5F5ycJff^jE?^z2 zE4EGoXFJCy#`fSD-0oD_Jif}OhElDXf*GRg{C=|SHf}#UXFH1NURw3jYE-y$CAO0F zFS_nyO-5(urvufjSNLvl$kHe(n`Glied}u`{9-H=;s)vn$)!pEMF=G?@jQG?5-r@! z&n{~I{aNZz@*yQ9jAq?XFY{>H`s36zHWlU_Zn->LRxygd;-l!A%E9HLb2FXWh8eNP z89|sImzqR0P!%12nR+30AB2^KVT)o-+aM{W0>8{6x#8dl!vU8KM4s+v_$LayV&q zeKgEwJC<&IoaKM9(VaSz)vDAhjLc2+WRn|ALMVwXTL;qZN)v_2goB}QBW?oADv_38 z!Urk3dk4-7Z0`cnvf4!w^FG=9*dT#jD{(e}owux<9EIkqT%k-(yCNb#9l1-LltRjA zYb$#k=f}m|KNLtCDjr;!P8t_^ML7%v@|8Lixi`nrZdU>d$(79o8_nSEO{Rl9@tW}3 zMvMkcx1?R;IVFG8L+0i;x5CVCHsi=qZ$gAl7slD%0#eIuV)Gv^SN%(u&nXwdpPx!M znytM$r>>pj+J1K5GOTm?;GQhtsg_NZZT{)dEw8epn>`<;HZdD{{zRs;ngahGOdQs= zmOS)g*F~c4w>jfLIyX)2YiP71wIHZrZ){5VXC4LD?1N$gyHfUxuTx&gJ|6|m$I>XB z$9{{;hz(;vH@N$(EWR%*KFC$GbVi^1d6w0Mi#Pd*eobBqp&S$T^xUvfi^|?YiNgTi zIeEZm8-8K;iJ!B~rgf($!izuF6?6(Keiu7K#KLNJsKX1}q3w%XEyBPig^B`9ta&kG zj_&&%^YP@Kxcm85FN4`x;IO_pR|Dd7QfH`3sq}F}zJuAh*>IUrzbB+Ztx~dj146 z-AfIF{7sg_kUEWDDE)tFwf<}V=)j}SWI#>JEPE*cC;~zIN6YVRh-h4b`Sd4CgUas$L zM;f46V_Oh~fWW8T8oMr{#i;hGAEh1D<=SN{LxIcr#mKhX-aEe)?0pgx+gs>-{L`az z=)67io}*6HGsgsh#$?>-4Q*TR-!?$B++${oqG}jhXxU~(6A7{_7^+Z%K@MX5kQBji zL`fC^JT^(WI%Bb!Nda3W?d zd^q{~yv#|JT5<-zQxphJCA&piZ_|i2x<#PXdCMgUi;eEI-P-9rdRA|xHZ~+_31szL|`g@4-E8x=EQj{QwfOq7^BX zDJG`hwr1#e5nTYbvz%W|&c|*UohKI7E5te3QH$E&)nCHoZdNic@PCj9_i8OgI>tP3 z#~$x+zft>aQ6ulnCensR%OnGpT_puN`fc$=bNvmMfsk@B6f2m zk1p#bCUxX^7!my}T5JFKpk!i>AzF&mi-?L8Qi{x=+fw3v^2{+8c-?1 z7jW*?6-SW6+G<#^>4oIPPvIYm0*c)sbT=_8ndFwRpc$t{W?9L^aD1y)u7c?GH0(dm zkNkz~d%tb{c6Xlcg86hpQX>C`y+!mO5>XOS61@}r#{2WTzJKm@uAP}{&dfRU+|PZ#UcEzE+c9oV z1RBNe(6(->b^*2gkNZGDtCz=+78K5)5f$~fxF=o38i#37F#8LKTmOw}msKfZKSY4T z3@%PSrhD#OGt_*wWKX-T$X|+w$EJP50$5xaF?eLx`; zgNCQ)4=;rRs{^tA*)&}l9F;QPNc~gV~(yszl1PMt_J~@NuTS> z6T{2%F6Ku6KjJF!GjAUK9wAM2vVQE~4u`9qhv(2h|CIm>8_QCt7vJA;EHq}~MfBTI zlh`TD+mK4hV3|_Uz86^o=?%2y zjqlZKq7M;up~}G|3hHd39Vz8{g}kUPzP0@3th(>!(a8x7WHv%;2elXUvL)Sm`x7DUa&k?oa#8scX+Fqle#O*+);nvmJEKhbGqFt`DakRND1KOWJm# zy`-$}p1ha3oVcIf=n*@J9C5R1{#90d{%}RWdt_q`4*WJ|!c(fDzL^dAoW|ft5f59_QQ}aa<*{>vIr!BY>D_fh#oPbaip{vX5 z>7+dpey>XuaEe1rE=O2A-Zy4^n_dTQW4w){VI6a^c6x~Mw>ci zi}RjWi<4<;hU;5?*Q(8aZyrPtA3IAtx!@f%{jSdtsrw43L>DhGW^NL6Pq{6~k6_vK zJVarW;Y8$fyn+eVA!n=Th%l0XC+(K>M*>J;w0iCdO$lo`@#pA?Hk$^QpXEL7<9vTL z+f53yOqE{wS(a@#ADx)%6#(!F0s`k}EE)-+k@-a*%OtXf-!`peMe^;tK0eA#%rGjK z^ZhGscHes;i}d-W8(@cQ%quxWoP2*4$$kr`ofcLZw4co{^<*MVeepzSq4YOkA7CAp zv%5QRV$@e&^KiOo$fH$fYP+(1zXi$80ye4`oQWQ5zfXk;%r77Pu2~p*reM|a^y|+W ztm`X@Q0S}AB=kGuKW}k(aWD`+U+ToqmGiCIJ?hR47D{)t4N(j}O*-1TmwiqVlsEIF zVY2Si@zcqVX|iD=Rvmg8dbt&<1GOPPlKNZNyfrr#Ne!<#hUXY_oM_}#%< z>6(u(>2_r>j}uHFS*(dn`y>89KV~$--n`8(O!q1SoT0=H!%E&_K4Q7NqAHw9YMwy!`!BGHALX$zh@%8VibcvlRif3QgC^QZ30&Y9F@Z#aPbAvGS(2 zzgq;hzWLcgXVe>AS4A5`e~3u$SA0H8+n3!!3P`!ld_U;h#{f{4C8Z+s#MZ;fH-!`D zQizGyHo>|zxn+|JJ?Br0hq-;$#xMBGES)sIKAzi~jqE^Mp_eMYm+8=N&KtY>X80*- z+mb<)-9O=6_U`y4l_o8)mG-?6nGzM$(;u@|4Sx$KE0*0y#i8wJ*1vmzV@h+?_l6SJ zH|P;Xf9~oS6yj4cy+af*QM!;V-jV{D=8Yk)^8jv%2*j}JJ&(z-CaG=@uT`D30;TC8 zi-t=FyS4PKH!qE~K=#C8rKPk1OR0;7m(ls_funfAr@J)aHqMf%2*}*YVyI|Pe}6I? z+kNWCll1&yE(AYov>`I;EY8|M*8vi?M7=&)!^SKq&(Co`Pcc}oW!UF=lR2sT#$;Ke z+pI_c;Qeu82i+ah7l+?18dSaE&nQx(aos>CbOY#laB2!qsJLQe-g$ge_c95&I$M){ z!SyIP!{`QPqpbw7BB3T1zg|)foZemE%>8>G!#)m#f;Yc5ysqz5U7FXVu~+`Rv-l+W zHIpM}P|gIKLO=%3g96<#`2)-$=TUf*qpfvkT(qGya^(?^^D^HLwTJ6Qe&h)Q+u@B> zwMqc#<#A-^Qdfj*!gH%-!qrhz-I8zqwb29NgSK7xw|dR?EhjG-JFGFhisU?)zw~CF zUK(bq*Xt$8>!}*7%ks}Rw&MpK)WDj=bf4=8Bx<{he?^#3bP8@y_4tQ(4ioZUkHZ3s zOddr2@{Uu=rf%8kxIiwMYw|4f2tKN<_m+DiK_XuNbXz#}i;4%QKMUQi%v>q zPVe*`Lf?&izaLsavJVv)lP%24m_5qYo|%&rl`l&O6V$Ph{&n^Ht2GNTs{~gQDm$^k z)B?W5-W~;03u3=apY5e1*!>U$pmF%`-FpSt$F}^|$YpkGq^NImJ-9-y=>1Yw-L>+J zL7=iT2dGN>yNmbMhR`J`1K$VOLUn<>UpZi4I$QW?)&@|76cD~}s!@N7^&JdVa5vg+ zLs0MGzg6Y6WLo0>xUTr=J^S6K!NAjys*dM+7gYawWTX@@-c@;0?q7d3Yo$q-`C|H& zTwpS*k`GqP0CgR3oM7meihUjLy*h;c`KRtJFBRkzPDg+lVdJN$@|!{1UTQOH6jQ>>IWmiyiIvX7=|B#1y}H$C z-cr+?v-fghs>QwYvn9!j0)*`Tgwv~eK90}vk#`aK+DWKh|zZP@H?ButZ2zsc3u4toN062NgnJvtoRG2}|io7fWDH3QtaA zRtB>?Hk6r?l|w(h&FdtKPpViHMFr|bDw16i-01!_ZM*f4$TH{ghs8@hN7?5c(Cph3 zcj4a8Fx)e`r`%Dprt1S8InR9iK1^hP|K*+Ey{12?^NEl~|^M)ug*^}Ng(vGf%FJxf^%s{J0%g7>)WLBTl*^(l?e_KyNEp%^@O!;B5 zAOcFP4)(&?yzVg0ncQsbJ0Gl@P-@utkBgm;g`2Z=J1CVaqXe|0t!paTBUx52(!Y)T zS^!N;Kc=(Zt_K8FojxXo^6a8&f_=GapE`9H7OaccJ$g{6n12o=*zU;rj2li-%+kFl zc5cz&VHN@?TrH|g%Q0-_LGRk~flczmRIBW}kjnrXJBLm*eydBm0W$|9jV0X5HQ<^MxZsGn>vL=lzgZFI}PmAYa_k94r(gLo}4?K-`Uknt?#g zdCQAO-4Ua&|Ki=$n# zOQ~Pp7pP5FyIJH{z%qI94^aYm0JWG+%0#U@^%|xjUih+_u_;!DUup&5w^(#?H z-jUQU;0D^NCU}{${x_*7Kk^jMenFE6_2p&ID4q-p_-K>cEZLDA-wrgvcYuXZRo@a+ zdONk!E(Je1=t(ajRVqJ{UC)waC!jz!zUZB+G4GGftTiqIbtHN*IQ1)DE-5S2#700) zges0a{iD4@T%Dff<;~l zZow0Cb0eF$>mDj;X96=WKT~wizONOY%o~Ks#BufCenZo)ZKu2_drFq}LUsw@w(w>* z0%#`Scl$Kj_rNCPEQ)xSDNi#BLlW6{3+PwF+u-MaFu5X7lB4Z)yI638I{&pGI=*(9kXAoDDQQN#u&dODpSl zuy??DYg1cuO24Wwmx6L<{FbtpPCikG&(-hR1eIxbqa*4&C`d&xMsYidkM96Ge(fNJ z;@O4%IKX=8Yqe0IW1k*av1O;2jO(57;|#a*u%7DAodae=`}K>bOl091vxT!A?7{r}@gQhA zD%;shgPTy$c9B-+`8$ggTeSpIHgZS#XTPd@6ocTm5%B7))gC|0&~LY{`K|oKXZoof zd4}bj^o*r&EnEcNbQn<#6*8kihanKZNI2U8Aht{b zpL|itbH*kKu399NEYu;7YRvv{D$1y#c+Q2=+=rXs_m`67h>7=D*p^gID_Zu>(A;ik z`#e|`A&Z$_jt;*gqEgZ9`4S_TZj{~gG4dIe*EK_&<2Kx8KRb({X}rZLw03w$R6?>P zQM(k}hR6RTyDZ^h&}co+_%GEi91N`cZpw;@Phr8h=`u$~`k9Z}w8{cy zzx?8IYc6w?pslR|Hyde6OZSzId&%7YlB#_Dwy8=tG^<0^7nrj3gS!?%JO}^6;@_Az z&5ffw@`e3!sB$7{P5zdCq5`A{vRk5-EcURc`FrjrqNQ+YoawZ?NLKU-FEJ=ojh$Pi z9V*mkT&d#Srl%hcPna&cf*jG^!Fy+0VD^SKX+?-n=^V{@R7d3X!?RErX1 z+epY0f$6pv&R8U=Sjs0w{Q2@Tj#9_EW2xuD#c6}QN)nSV^CR9DR?Frf-Qj$Ht*WOL z0brU^DWIFZYXt%Ey4m%_IQ7iW^MxhVL`;Yf`96m>Ja$vEWts3jO_;{TfK2qqH&G7C z{bC(h^JaQtZ~39yi~Z9^vjx$Ie!DgwQl!T|4#R?#gOp6C*~XSX_3;iV9p}6A^qWwq;;eV%B7*~r1Ir&qL&x>F-wb(gdB;>U zwy-RdH$3}$&wc*ysBBM0h9BOK|NMiE{6c0>SoVs7C*#Q%MMDL=1rlPClPwj;DmokC z;4XIZu{_)>nxCt?o$UElq-Q`1BEID*-Y7cqp`vcX-oN^2{$eRF1gK=96){5V;ct4% z8IO>k^+HCBw=%o_>5NEa=Z)nvMcW0bozOutf)^Arx}dfQ;P>FMz0 zJ!4M~uU&XFE`kckCcG*nGBUV=NOF{Iet~9XR@rRp#rh1Ce;>Q@@S16CMEfkEM{%aJXbqgS1w+ zrd@0Ws|SjKC;YYE!?>n!{Gb}g^yBp~F zQ!G-8iHE~~vFT+u%Ot%~vDyyNHoYGINGOhrVyxS>NkUjYrbVWLAb62BEhWyT!>rhT zlHSi!GeSO=f;(zYNj3IPSEGNC9U_V3)%b$y{z?`eg@z!%ZpFP*RkBab_J#Isj`}!T zxH@T=jqN|B7kFER=3U==S*SvB$Zn>l^{=L@JA)F*2tBr|i$%%Q2>NNKN`qT9X1H+k zckV;@5EV!NUnMkVHoXBdz!$t@^}li0orf=;Iy-wnO5t#N(;jKIF9C|TEbu~x^0_<`4ilPZ>$$a8_k$_`mIiLDmL$9Zlm_$%YDY+DzSA)lEV6I9Nikdu>-YhvJp{ZO};_#pb>Rk0SKM+KwwfP<=q+&(ob6dSU!BOps=~{bP!I+9G zGUF-9mzqCwnfJOl@d~JTc!neh%tCbr2fZq~GSarb=5xlBiE;J(Ti6NV?N;GmvutLT zb}usj$I96qg8+SI7L|S5m79q&{H%O`2Di;U)BPSl)kkOd^>x=vV{9JzG2Q-o3ak^O z{ZQ5G2eEzq$^&L>lH79)z1Kq>v@k}BTMLB2;FSihiIx66vE^N(ydyUr!xl1FWyY&2 z@_Muz3@$$y{P$P1O9r6w2;yviF&JLW(0p12ShK6I*=!eaMv;h-?T4JEBT9xwLxfBr zmLXp~N#^KY(AclHQLYoRspgywnO}3xsom`Py0{C{>*x0q}dmj za#b7ew=8?F>gvC_7=E~JsW@{C9KVUTRyh>ck<9%Lz=x>JMP!9b#F;&D_s4gD09Ei~ zrJI82cT3Niq{n3OOYB3XK5*aNTt0K#gq??5%$#OECHyHg!-)smv;ybM%hWD zr5S?eFb25)*YmjVe->5yNr1Z-CaC^EwozK%?a zY8h68OC3hqZo|Ua%Ep)7NF|d3ewC);e_g!{2KX8J5*>Kagyotn@iLl#WFv>fb3L7R z?~`=y1;6WciV196WJ=RNIB+y{Pw0nmjJ)z!c!+$j|3b7{En)Cn;1e@0TjWkoBDnO& zH;N7N2qkZyNZER@EYSnUL`@LlxkDrVZOynoajFv#l$J0{^~Yl-|LLB~Uv_cmk6^w^ z=lN#O``1A1J9;_a$5*t*&ITGeKoAswkRzFl!(MZ1lPObk{9GJq=r31(S4N)VdeAS_aCpSSCdr?~n1HxBE98nET z_HN~_DT8Mpw@I`+`WEW}@v0mo!_6`J}z33Ock%xfHkzD7JPH1>XN<6)@li?)>1W>&_k`jihh1eda>=Y$f6NF^v}MW&+e>(DcV zuyA6nSOi}b_|WxBU~XgYGj7}RdUM*UKoR1hTPQPm=0YU17OBvDK}F0}tl1Ozknfq{ z2kBL;`De<641WD&C4tD^L5sG$IO$ZzL*zP`H1sxY703W}~$Xq_2~Y$F;6!=Y zf5h@>7}`j{^YY!am52^eEIOl;0HkNl^WP9H$g}>{@J{G;VDnjbhWOJr3Iu6kF3FYE z27o>jD+UrUuHs6w1{5mZ{s7jm+uFo*fFQZR4M-Jy%1StOOxiX>F*87`UJ)Bo*9Oq& z0U9?F(R515)TgHUT?PK3&hD6|w)7i)qz}w49t>C~hd3noBe(cpO)Q>hw9_b5sR}@Jl$yYpw3h;~Ypdkb>|4Q;La@q0avZzrZ!%bF!y& z+vSOii-VvNDuIbV)Q37TS0Ja0n=vaGIt4H@N`{rzqJQ^SB=vgdA82Z^k+dr)er0WE z$I~!rihK*y`6vs^ml6yF^(rGj|5@xaUR*DAI&`{Cr@0cF6c~<)n$@FCYAq{)D2pBk zwFFqnxAdj{F)$&`^DA%j@M#8@BXOiG_=?ufhP!_@0EhQ%1vhkz&fS&#)6@SW$d8i~ zJiVnW`atEKDsWeC81nyA0JZnyQd-o_mEOS2MFDBMJc}3h?SC`IMGyi5R3aO%q3-{G zlp=ZnBtP{jJ!R{E40;IuN9q)c5z!V$;FSk%FFOrz7C_Qn)f#+$!>GUcn+;D^10k7X^e24p7j`#jIj_ec#T_N 0 + O = O[mask] + E = E[mask] + + # Compute chi-squared statistic + chi2_stat = ((O - E) ** 2 / E).sum() + + # Move dof and chi2_stat to the same device + dof = torch.tensor(dof, dtype=torch.float32, device=chi2_stat.device) + + # Compute p-value using the survival function (1 - CDF) + p_value = torch.special.gammaincc(dof / 2, chi2_stat / 2).item() + + return chi2_stat, p_value, dof, E + + def _pqm_test( - x_samples: np.ndarray, - y_samples: np.ndarray, + x_samples: torch.Tensor, + y_samples: torch.Tensor, num_refs: int, - whiten: bool, - x_frac: Optional[float] = None, - gauss_frac: float = 0.0, + z_score_norm: bool, + x_frac: Optional[float], + gauss_frac: float, + device: str, ): """ Helper function to perform the PQM test and return the results from @@ -38,14 +122,22 @@ def _pqm_test( Parameters ---------- x_samples : np.ndarray - Samples from the first distribution, test samples. + Samples from the first distribution. Must have shape (N, *D) N is the + number of x samples, and D is the dimensionality of the samples. y_samples : np.ndarray - Samples from the second distribution, reference samples. + Samples from the second distribution. Must have shape (M, *D) M is the + number of y samples, and D is the dimensionality of the samples. num_refs : int - Number of reference samples to use. - whiten : bool - If True, whiten the samples by subtracting the mean and dividing by the - standard deviation. + Number of reference samples to use. These samples will be drawn from + x_samples, y_samples, and/or a Gaussian distribution, see the note + below. + re_tessellation : Optional[int] + Number of times pqm_pvalue is called, re-tesselating the space. No + re_tessellation if None (default). + z_score_norm : bool + If True, z_score_norm the samples by subtracting the mean and dividing by the + standard deviation. mean and std are calculated from the combined + x_samples and y_samples. x_frac : float Fraction of x_samples to use as reference samples. ``x_frac = 1`` will use only x_samples as reference samples, ``x_frac = 0`` will use only @@ -58,6 +150,21 @@ def _pqm_test( support of the reference samples if pathological behavior is expected. Default: 0.0 no gaussian samples. + Note + ---- + When using ``x_frac`` and ``gauss_frac``, note that the number of + reference samples from the x_samples, y_samples, and Gaussian + distribution will be determined by a multinomial distribution. This + means that the actual number of reference samples from each distribution + may not be exactly equal to the requested fractions, but will on average + equal those numbers. The mean relative number of reference samples drawn + from x_samples, y_samples, and Gaussian is ``Nx=x_frac*(1-gauss_frac)``, + ``Ny=(1-x_frac)*(1-gauss_frac)``, and ``Ng=gauss_frac`` respectively. + For best results, we suggest using a large number of re-tessellations, + though this is our recommendation in any case. + device : str + Device to use for computation. Default: 'cpu'. If 'cuda' is selected, + Note ---- When using ``x_frac`` and ``gauss_frac``, note that the number of @@ -76,8 +183,8 @@ def _pqm_test( tuple Results from scipy.stats.chi2_contingency function. """ - nx, *D = x_samples.shape - ny, *D = y_samples.shape + nx = x_samples.shape[0] + ny = y_samples.shape[0] if (nx + ny) <= num_refs + 2: raise ValueError( "Number of reference samples (num_ref) must be less than the number of x/y samples. Ideally much less." @@ -86,7 +193,7 @@ def _pqm_test( warnings.warn( "Number of samples is small (less than twice the number of reference samples). Result will have high variance and/or be non-discriminating." ) - if whiten: + if z_score_norm: mean, std = _mean_std(x_samples, y_samples) y_samples = (y_samples - mean) / std x_samples = (x_samples - mean) / std @@ -94,48 +201,73 @@ def _pqm_test( # Determine fraction of x_samples to use as reference samples x_frac = nx / (nx + ny) if x_frac is None else x_frac - # Determine number of samples from each distribution (x_samples, y_samples, gaussian) - Nx, Ny, Ng = np.random.multinomial( - num_refs, - [x_frac * (1.0 - gauss_frac), (1.0 - x_frac) * (1.0 - gauss_frac), gauss_frac], + # Determine number of samples from each distribution + probs = torch.tensor( + [ + x_frac * (1.0 - gauss_frac), + (1.0 - x_frac) * (1.0 - gauss_frac), + gauss_frac, + ], + device=device, + ) + + counts = Multinomial(total_count=num_refs, probs=probs).sample() + counts = counts.round().long() + Nx, Ny, Ng = counts.tolist() + assert (Nx + Ny + Ng) == num_refs, ( + f"Something went wrong. Nx={Nx}, Ny={Ny}, Ng={Ng} should sum to num_refs={num_refs}" ) - assert (Nx + Ny + Ng) == num_refs, f"Something went wrong. Nx={Nx}, Ny={Ny}, Ng={Ng} should sum to num_refs={num_refs}" # fmt: skip # Collect reference samples from x_samples - xrefs = np.random.choice(nx, Nx, replace=False) - xrefs, x_samples = x_samples[xrefs], np.delete(x_samples, xrefs, axis=0) + x_indices = torch.randperm(nx, device=device) + if Nx > nx: + raise ValueError("Cannot sample more references from x_samples than available") + xrefs_indices = x_indices[:Nx] + x_samples_indices = x_indices[Nx:] + + xrefs = x_samples[xrefs_indices] + x_samples = x_samples[x_samples_indices] # Collect reference samples from y_samples - yrefs = np.random.choice(ny, Ny, replace=False) - yrefs, y_samples = y_samples[yrefs], np.delete(y_samples, yrefs, axis=0) + y_indices = torch.randperm(ny, device=device) + if Ny > ny: + raise ValueError("Cannot sample more references from y_samples than available") + yrefs_indices = y_indices[:Ny] + y_samples_indices = y_indices[Ny:] + + yrefs = y_samples[yrefs_indices] + y_samples = y_samples[y_samples_indices] # Join the full set of reference samples - refs = np.concatenate([xrefs, yrefs], axis=0) + refs = torch.cat([xrefs, yrefs], dim=0) - # get gaussian reference points if requested + # Get gaussian reference points if requested if Ng > 0: m, s = _mean_std(x_samples, y_samples) - gauss_refs = np.random.normal( - loc=m, - scale=s, - size=(Ng, *D), + gauss_refs = torch.normal( + mean=m.repeat(Ng, 1), + std=s.repeat(Ng, 1), ) - refs = np.concatenate([refs, gauss_refs], axis=0) + refs = torch.cat([refs, gauss_refs], dim=0) - # Build KDtree to measure distances - tree = KDTree(refs) + num_refs = refs.shape[0] - idx = tree.query(x_samples, k=1, workers=-1)[1] - counts_x = np.bincount(idx, minlength=num_refs) + # Compute nearest reference for x_samples + distances = torch.cdist(x_samples, refs) + idx = distances.argmin(dim=1) + counts_x = torch.bincount(idx, minlength=num_refs) - idx = tree.query(y_samples, k=1, workers=-1)[1] - counts_y = np.bincount(idx, minlength=num_refs) + # Compute nearest reference for y_samples + distances = torch.cdist(y_samples, refs) + idx = distances.argmin(dim=1) + counts_y = torch.bincount(idx, minlength=num_refs) # Remove reference samples with no counts C = (counts_x > 0) | (counts_y > 0) - counts_x, counts_y = counts_x[C], counts_y[C] + counts_x = counts_x[C] + counts_y = counts_y[C] - n_filled_bins = np.sum(C) + n_filled_bins = C.sum().item() if n_filled_bins == 1: raise ValueError( """ @@ -155,17 +287,20 @@ def _pqm_test( """ ) - return chi2_contingency(np.stack([counts_x, counts_y])) + # Perform chi-squared test + counts = torch.stack([counts_x, counts_y]) + return _chi2_contingency(counts, device) def pqm_pvalue( - x_samples: np.ndarray, - y_samples: np.ndarray, + x_samples, + y_samples, num_refs: int = 100, re_tessellation: Optional[int] = None, - whiten: bool = False, + z_score_norm: bool = False, x_frac: Optional[float] = None, gauss_frac: float = 0.0, + device: str = 'cpu', ): """ Perform the PQM test of the null hypothesis that `x_samples` and `y_samples` @@ -187,8 +322,8 @@ def pqm_pvalue( re_tessellation : Optional[int] Number of times pqm_pvalue is called, re-tesselating the space. No re_tessellation if None (default). - whiten : bool - If True, whiten the samples by subtracting the mean and dividing by the + z_score_norm : bool + If True, z_score_norm the samples by subtracting the mean and dividing by the standard deviation. mean and std are calculated from the combined x_samples and y_samples. x_frac : float @@ -215,6 +350,8 @@ def pqm_pvalue( ``Ny=(1-x_frac)*(1-gauss_frac)``, and ``Ng=gauss_frac`` respectively. For best results, we suggest using a large number of re-tessellations, though this is our recommendation in any case. + device : str + Device to use for computation. Default: 'cpu'. If 'cuda' is selected, Returns ------- @@ -222,30 +359,40 @@ def pqm_pvalue( pvalue(s). Null hypothesis that both samples are drawn from the same distribution. """ + # Move samples to torch tensors on the selected device + x_samples = torch.tensor(x_samples, device=device) + y_samples = torch.tensor(y_samples, device=device) + if re_tessellation is not None: return [ pqm_pvalue( x_samples, y_samples, num_refs=num_refs, - whiten=whiten, + z_score_norm=z_score_norm, x_frac=x_frac, gauss_frac=gauss_frac, + device=device, ) for _ in range(re_tessellation) ] - _, pvalue, _, _ = _pqm_test(x_samples, y_samples, num_refs, whiten, x_frac, gauss_frac) - return pvalue + chi2_stat, p_value, dof, _ = _pqm_test( + x_samples, y_samples, num_refs, z_score_norm, x_frac, gauss_frac, device + ) + + # Return p-value as a float + return p_value if isinstance(p_value, float) else float(p_value) def pqm_chi2( - x_samples: np.ndarray, - y_samples: np.ndarray, + x_samples, + y_samples, num_refs: int = 100, re_tessellation: Optional[int] = None, - whiten: bool = False, + z_score_norm: bool = False, x_frac: Optional[float] = None, gauss_frac: float = 0.0, + device: str = 'cpu', ): """ Perform the PQM test of the null hypothesis that `x_samples` and `y_samples` @@ -265,10 +412,10 @@ def pqm_chi2( x_samples, y_samples, and/or a Gaussian distribution, see the note below. re_tessellation : Optional[int] - Number of times pqm_pvalue is called, re-tesselating the space. No + Number of times pqm_chi2 is called, re-tesselating the space. No re_tessellation if None (default). - whiten : bool - If True, whiten the samples by subtracting the mean and dividing by the + z_score_norm : bool + If True, z_score_norm the samples by subtracting the mean and dividing by the standard deviation. mean and std are calculated from the combined x_samples and y_samples. x_frac : float @@ -282,6 +429,8 @@ def pqm_chi2( determined from the combined x_samples/y_samples. This ensures full support of the reference samples if pathological behavior is expected. Default: 0.0 no gaussian samples. + device : str + Device to use for computation. Default: 'cpu'. If 'cuda' is selected, Note ---- @@ -312,26 +461,31 @@ def pqm_chi2( float or list chi2 statistic(s). """ + # Move samples to torch tensors on the selected device + x_samples = torch.tensor(x_samples, device=device) + y_samples = torch.tensor(y_samples, device=device) + if re_tessellation is not None: return [ pqm_chi2( x_samples, y_samples, num_refs=num_refs, - whiten=whiten, + z_score_norm=z_score_norm, x_frac=x_frac, gauss_frac=gauss_frac, + device=device, ) for _ in range(re_tessellation) ] - chi2_stat, _, dof, _ = _pqm_test(x_samples, y_samples, num_refs, whiten, x_frac, gauss_frac) + chi2_stat, _, dof, _ = _pqm_test( + x_samples, y_samples, num_refs, z_score_norm, x_frac, gauss_frac, device + ) + + # Rescale chi2 statistic if necessary if dof != num_refs - 1: - # Rescale chi2 to new value which has the same cumulative probability - if chi2_stat / dof < 10: - cp = chi2.sf(chi2_stat, dof) - chi2_stat = chi2.isf(cp, num_refs - 1) - else: - chi2_stat = chi2_stat * (num_refs - 1) / dof - dof = num_refs - 1 - return chi2_stat + chi2_stat = rescale_chi2(chi2_stat, dof, num_refs - 1, device) + + # Return chi2_stat as a float + return chi2_stat.item() if isinstance(chi2_stat, torch.Tensor) else float(chi2_stat) \ No newline at end of file diff --git a/src/pqm/test_gaussian.py b/src/pqm/test_gaussian.py index b67398b..f99cd6d 100644 --- a/src/pqm/test_gaussian.py +++ b/src/pqm/test_gaussian.py @@ -8,6 +8,6 @@ def test(): y_samples = np.random.normal(size=(500, 50)) x_samples = np.random.normal(size=(250, 50)) - new.append(pqm_pvalue(x_samples, y_samples)) + new.append(pqm_pvalue(x_samples, y_samples, device = 'cpu')) assert np.abs(np.mean(new) - 0.5) < 0.15 From 4cbe01dd841af052a5aa6e1708e0123e49061959 Mon Sep 17 00:00:00 2001 From: Sammy Sharief Date: Mon, 28 Oct 2024 14:35:25 -0400 Subject: [PATCH 2/7] Updated test_gaussian to change from whiten to z_score_norm --- tests/test_gaussian.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/tests/test_gaussian.py b/tests/test_gaussian.py index f4d011c..0006746 100644 --- a/tests/test_gaussian.py +++ b/tests/test_gaussian.py @@ -3,16 +3,16 @@ import pytest -@pytest.mark.parametrize("whiten", [True, False]) +@pytest.mark.parametrize("z_score_norm", [True, False]) @pytest.mark.parametrize("num_refs", [20, 100]) @pytest.mark.parametrize("ndim", [1, 50]) -def test_pass_pvalue(whiten, num_refs, ndim): +def test_pass_pvalue(z_score_norm, num_refs, ndim): new = [] for _ in range(50): y_samples = np.random.normal(size=(500, ndim)) x_samples = np.random.normal(size=(250, ndim)) - new.append(pqm_pvalue(x_samples, y_samples, whiten=whiten, num_refs=num_refs)) + new.append(pqm_pvalue(x_samples, y_samples, z_score_norm=z_score_norm, num_refs=num_refs)) # Check for roughly uniform distribution of p-values assert np.abs(np.mean(new) - 0.5) < 0.15 @@ -45,17 +45,17 @@ def test_fail_pvalue(num_refs, ndim): assert np.mean(new) < 1e-3 -@pytest.mark.parametrize("whiten", [True, False]) +@pytest.mark.parametrize("z_score_norm", [True, False]) @pytest.mark.parametrize("num_refs", [20, 100]) @pytest.mark.parametrize("ndim", [1, 50]) -def test_fail_chi2(whiten, num_refs, ndim): +def test_fail_chi2(z_score_norm, num_refs, ndim): new = [] for _ in range(100): y_samples = np.random.normal(size=(500, ndim)) y_samples[:, 0] += 5 # one dim off by 5sigma x_samples = np.random.normal(size=(250, ndim)) - new.append(pqm_chi2(x_samples, y_samples, whiten=whiten, num_refs=num_refs)) + new.append(pqm_chi2(x_samples, y_samples, z_score_norm=z_score_norm, num_refs=num_refs)) new = np.array(new) assert (np.mean(new) / (num_refs - 1)) > 1.5 From 99a8e9b673720092afd31dc4c6f2276fd0c7701d Mon Sep 17 00:00:00 2001 From: Sammy Sharief Date: Sun, 3 Nov 2024 19:27:19 -0500 Subject: [PATCH 3/7] Updated PQM with better handling between CPU (numpy) and GPU (torch) --- src/pqm/pqm.py | 562 ++++++++++++++++++++++++++++++++++--------------- 1 file changed, 388 insertions(+), 174 deletions(-) diff --git a/src/pqm/pqm.py b/src/pqm/pqm.py index 1460760..8b40432 100644 --- a/src/pqm/pqm.py +++ b/src/pqm/pqm.py @@ -3,12 +3,12 @@ import torch import numpy as np from scipy.stats import chi2_contingency, chi2 -from scipy.spatial import KDTree +from scipy.spatial.distance import cdist from torch.distributions import Multinomial +from typing import Optional, Union, Tuple __all__ = ("pqm_chi2", "pqm_pvalue") - def _mean_std(sample1, sample2, dim=0): """Get the mean and std of two combined samples without actually combining them.""" n1 = sample1.shape[dim] @@ -29,15 +29,17 @@ def _mean_std(sample1, sample2, dim=0): ) return m, s - def rescale_chi2(chi2_stat, orig_dof, target_dof, device): """ Rescale chi2 statistic using appropriate methods depending on the device. """ - - # Move tensors to CPU and convert to NumPy - chi2_stat_cpu = chi2_stat.cpu().item() # Convert to float - orig_dof_cpu = orig_dof.cpu().item() # Convert to float + if device.type == 'cuda': + # Move tensors to CPU and convert to NumPy + chi2_stat_cpu = chi2_stat.cpu().item() # Convert to float + orig_dof_cpu = orig_dof.cpu().item() # Convert to float + else: + chi2_stat_cpu = chi2_stat + orig_dof_cpu = orig_dof if orig_dof_cpu == target_dof: return chi2_stat_cpu @@ -48,83 +50,316 @@ def rescale_chi2(chi2_stat, orig_dof, target_dof, device): return chi2.isf(cp, target_dof) else: # Use simple scaling for large values - return chi2_stat_cpu * target_dof / orig_dof_cpu - + return chi2_stat_cpu * target_dof / orig_dof_cpu +def _chi2_contingency_torch( + counts_x: torch.Tensor, + counts_y: torch.Tensor +) -> Tuple[torch.Tensor, float, torch.Tensor, torch.Tensor]: + """ + Perform chi-squared contingency test using PyTorch tensors. + + Returns: + chi2_stat (torch.Tensor): Chi-squared statistic. + p_value (float): p-value. + dof (torch.Tensor): Degrees of freedom. + expected (torch.Tensor): Expected frequencies. + """ + counts = torch.stack([counts_x, counts_y]) + + # Observed counts + O = counts.float() + + # Row sums and column sums + row_sums = O.sum(dim=1, keepdim=True) # shape (2, 1) + col_sums = O.sum(dim=0, keepdim=True) # shape (1, N) + total = O.sum() + + # Expected counts under the null hypothesis of independence + E = row_sums @ col_sums / total # shape (2, N) + + # Degrees of freedom + dof = (O.size(0) - 1) * (O.size(1) - 1) + + # Avoid division by zero + mask = E > 0 + O_masked = O[mask] + E_masked = E[mask] + + # Compute chi-squared statistic + chi2_stat = ((O_masked - E_masked) ** 2 / E_masked).sum() + + # Move dof and chi2_stat to the same device + dof = torch.tensor(dof, dtype=torch.float32, device=chi2_stat.device) + + # Compute p-value using the survival function (1 - CDF) + p_value = torch.special.gammaincc(dof / 2, chi2_stat / 2).item() + + return chi2_stat, p_value, dof, E -def _chi2_contingency(counts, device): +def _sample_reference_indices_numpy(Nx, nx, Ny, ny, Ng, x_samples, y_samples): """ - Computes the chi-squared statistic and p-value for a contingency table. + Helper function to sample references for CPU-based NumPy computations. Parameters ---------- - counts: torch.Tensor - 2xN tensor of counts for each category. - device : str - Device to use for computation. Default: 'cpu'. If 'cuda' is selected, + Nx : int + Number of references to sample from x_samples. + nx : int + Number of samples in x_samples. + Ny : int + Number of references to sample from y_samples. + ny : int + Number of samples in y_samples. + Ng : int + Number of references to sample from a Gaussian distribution. Returns ------- - tuple - chi2_stat, p_value, dof, expected + np.ndarray + References samples. """ - if device == 'cpu': - counts_np = counts.cpu().numpy() - chi2_stat, p_value, dof, expected = chi2_contingency(counts_np) - chi2_stat = torch.tensor(chi2_stat, device=device) - dof = torch.tensor(dof, device=device) - return chi2_stat, p_value, dof, expected - else: - # Observed counts - O = counts.float() - # Row sums and column sums - row_sums = O.sum(dim=1, keepdim=True) # shape (2, 1) - col_sums = O.sum(dim=0, keepdim=True) # shape (1, N) - total = O.sum() + if Nx > nx: + raise ValueError("Cannot sample more references from x_samples than available") + if Ny > ny: + raise ValueError("Cannot sample more references from y_samples than available") + + # Reference samples from x_samples + xrefs_indices = np.random.choice(nx, Nx, replace=False) + xrefs = x_samples[xrefs_indices] + x_samples = np.delete(x_samples, xrefs_indices, axis=0) + + # Reference samples from y_samples + yrefs_indices = np.random.choice(ny, Ny, replace=False) + yrefs = y_samples[yrefs_indices] + y_samples = np.delete(y_samples, yrefs_indices, axis=0) + + # Combine references + refs = np.concatenate([xrefs, yrefs], axis=0) + + # Gaussian references + if Ng > 0: + m, s = _mean_std(x_samples, y_samples) + gauss_refs = np.random.normal( + loc=m, + scale=s, + size=(Ng, ) + tuple(x_samples.shape[1:]) + ) + refs = np.concatenate([refs, gauss_refs], axis=0) + + return refs - # Expected counts under the null hypothesis of independence - E = row_sums @ col_sums / total # shape (2, N) +def _compute_distances_numpy(x_samples, y_samples, refs, current_num_refs, num_refs): + """ + Helper function to calculate distances for CPU-based NumPy computations. - # Degrees of freedom - dof = (O.size(0) - 1) * (O.size(1) - 1) + Parameters + ---------- + x_samples : np.ndarray + Samples from the first distribution. Must have shape (N, *D) N is the + number of x samples, and D is the dimensionality of the samples. + y_samples : np.ndarray + Samples from the second distribution. Must have shape (M, *D) M is the + number of y samples, and D is the dimensionality of the samples. + refs : np.ndarray + Reference samples. Must have shape (num_refs, *D) where D is the + dimensionality of the samples. + current_num_refs : int + Number of reference samples used in the test. + num_refs : int + Number of reference samples to use. - # Avoid division by zero - mask = E > 0 - O = O[mask] - E = E[mask] + Returns + ------- + tuple + Results from scipy.stats.chi2_contingency. + """ + + # Compute distances + distances_x = cdist(x_samples, refs, metric='euclidean') + distances_y = cdist(y_samples, refs, metric='euclidean') + + # Nearest references + idx_x = np.argmin(distances_x, axis=1) + idx_y = np.argmin(distances_y, axis=1) + + # Counts + counts_x = np.bincount(idx_x, minlength=current_num_refs) + counts_y = np.bincount(idx_y, minlength=current_num_refs) + + # Remove references with no counts + C = (counts_x > 0) | (counts_y > 0) + counts_x = counts_x[C] + counts_y = counts_y[C] + + n_filled_bins = np.sum(C) + if n_filled_bins == 1: + raise ValueError( + """ + Only one Voronoi cell has samples, so chi^2 cannot + be computed. This is likely due to a small number + of samples or a pathological distribution. If possible, + increase the number of x_samples and y_samples. + """ + ) + if n_filled_bins < (num_refs // 2): + warnings.warn( + """ + Less than half of the Voronoi cells have any samples in them. + Possibly due to a small number of samples or a pathological + distribution. Result may be unreliable. If possible, increase the + number of x_samples and y_samples. + """ + ) + + # Perform chi-squared test using SciPy + contingency_table = np.stack([counts_x, counts_y]) + return chi2_contingency(contingency_table) - # Compute chi-squared statistic - chi2_stat = ((O - E) ** 2 / E).sum() +def _sample_reference_indices_torch(Nx, nx, Ny, ny, Ng, x_samples, y_samples, device): + """ + Helper function to sample references for GPU-based Torch computations. - # Move dof and chi2_stat to the same device - dof = torch.tensor(dof, dtype=torch.float32, device=chi2_stat.device) + Parameters + ---------- + Nx : int + Number of references to sample from x_samples. + nx : int + Number of samples in x_samples. + Ny : int + Number of references to sample from y_samples. + ny : int + Number of samples in y_samples. + Ng : int + Number of references to sample from a Gaussian distribution. - # Compute p-value using the survival function (1 - CDF) - p_value = torch.special.gammaincc(dof / 2, chi2_stat / 2).item() + Returns + ------- + np.ndarray + References samples. + """ + + if Nx > nx: + raise ValueError("Cannot sample more references from x_samples than available") + if Ny > ny: + raise ValueError("Cannot sample more references from y_samples than available") + + # Reference samples from x_samples + x_indices = torch.randperm(nx, device=device) + xrefs_indices = x_indices[:Nx] + x_samples_indices = x_indices[Nx:] + xrefs = x_samples[xrefs_indices] + x_samples = x_samples[x_samples_indices] + + # Reference samples from y_samples + y_indices = torch.randperm(ny, device=device) + yrefs_indices = y_indices[:Ny] + y_samples_indices = y_indices[Ny:] + yrefs = y_samples[yrefs_indices] + y_samples = y_samples[y_samples_indices] + + # Combine references + refs = torch.cat([xrefs, yrefs], dim=0) + + # Gaussian references + if Ng > 0: + m, s = _mean_std(x_samples, y_samples) + # Ensure m and s have the correct shape + if m.dim() == 1: + m = m.unsqueeze(0) + if s.dim() == 1: + s = s.unsqueeze(0) + gauss_refs = torch.normal( + mean=m.repeat(Ng, 1), + std=s.repeat(Ng, 1), + ) + refs = torch.cat([refs, gauss_refs], dim=0) + return refs + +def _compute_distances_torch(x_samples, y_samples, refs, current_num_refs, num_refs): + """ + Helper function to calculate distances for GPU-based Torch computations. - return chi2_stat, p_value, dof, E + Parameters + ---------- + x_samples : torch.Tensor + Samples from the first distribution. Must have shape (N, *D) N is the + number of x samples, and D is the dimensionality of the samples. + y_samples : torch.Tensor + Samples from the second distribution. Must have shape (M, *D) M is the + number of y samples, and D is the dimensionality of the samples. + refs : torch.Tensor + Reference samples. Must have shape (num_refs, *D) where D is the + dimensionality of the samples. + current_num_refs : int + Number of reference samples used in the test. + num_refs : int + Number of reference samples to use. + Returns + ------- + tuple + Results from the PyTorch implementation of chi2_contingency. + """ + + # Compute distances and find nearest references + distances_x = torch.cdist(x_samples, refs) + idx_x = distances_x.argmin(dim=1) + counts_x = torch.bincount(idx_x, minlength=current_num_refs) + + distances_y = torch.cdist(y_samples, refs) + idx_y = distances_y.argmin(dim=1) + counts_y = torch.bincount(idx_y, minlength=current_num_refs) + + # Remove references with no counts + C = (counts_x > 0) | (counts_y > 0) + counts_x = counts_x[C] + counts_y = counts_y[C] + + n_filled_bins = C.sum().item() + if n_filled_bins == 1: + raise ValueError( + """ + Only one Voronoi cell has samples, so chi^2 cannot + be computed. This is likely due to a small number + of samples or a pathological distribution. If possible, + increase the number of x_samples and y_samples. + """ + ) + if n_filled_bins < (num_refs // 2): + warnings.warn( + """ + Less than half of the Voronoi cells have any samples in them. + Possibly due to a small number of samples or a pathological + distribution. Result may be unreliable. If possible, increase the + number of x_samples and y_samples. + """ + ) + + # Perform chi-squared test using the PyTorch implementation + chi2_stat, p_value, dof, expected = _chi2_contingency_torch(counts_x, counts_y) + return chi2_stat, p_value, dof, expected def _pqm_test( - x_samples: torch.Tensor, - y_samples: torch.Tensor, + x_samples: Union[np.ndarray, torch.Tensor], + y_samples: Union[np.ndarray, torch.Tensor], num_refs: int, z_score_norm: bool, x_frac: Optional[float], gauss_frac: float, - device: str, -): + device: str = 'cpu', +) -> Tuple: """ Helper function to perform the PQM test and return the results from - chi2_contingency. + chi2_contingency (using SciPy or a PyTorch implementation). Parameters ---------- - x_samples : np.ndarray + y_samples : np.ndarray or torch.Tensor Samples from the first distribution. Must have shape (N, *D) N is the number of x samples, and D is the dimensionality of the samples. - y_samples : np.ndarray + y_samples : np.ndarray or torch.Tensor Samples from the second distribution. Must have shape (M, *D) M is the number of y samples, and D is the dimensionality of the samples. num_refs : int @@ -149,19 +384,6 @@ def _pqm_test( determined from the combined x_samples/y_samples. This ensures full support of the reference samples if pathological behavior is expected. Default: 0.0 no gaussian samples. - - Note - ---- - When using ``x_frac`` and ``gauss_frac``, note that the number of - reference samples from the x_samples, y_samples, and Gaussian - distribution will be determined by a multinomial distribution. This - means that the actual number of reference samples from each distribution - may not be exactly equal to the requested fractions, but will on average - equal those numbers. The mean relative number of reference samples drawn - from x_samples, y_samples, and Gaussian is ``Nx=x_frac*(1-gauss_frac)``, - ``Ny=(1-x_frac)*(1-gauss_frac)``, and ``Ng=gauss_frac`` respectively. - For best results, we suggest using a large number of re-tessellations, - though this is our recommendation in any case. device : str Device to use for computation. Default: 'cpu'. If 'cuda' is selected, @@ -181,116 +403,84 @@ def _pqm_test( Returns ------- tuple - Results from scipy.stats.chi2_contingency function. + Results from scipy.stats.chi2_contingency or the PyTorch implementation. """ + + # Determine if we're working with NumPy or PyTorch + is_numpy = isinstance(x_samples, np.ndarray) and isinstance(y_samples, np.ndarray) + is_torch = isinstance(x_samples, torch.Tensor) and isinstance(y_samples, torch.Tensor) + + if not (is_numpy or is_torch): + raise TypeError("x_samples and y_samples must both be either NumPy arrays or PyTorch tensors.") + + # Validate sample sizes nx = x_samples.shape[0] ny = y_samples.shape[0] if (nx + ny) <= num_refs + 2: raise ValueError( - "Number of reference samples (num_ref) must be less than the number of x/y samples. Ideally much less." + "Number of reference samples (num_refs) must be less than the number of x/y samples. Ideally much less." ) elif (nx + ny) < 2 * num_refs: warnings.warn( - "Number of samples is small (less than twice the number of reference samples). Result will have high variance and/or be non-discriminating." + "Number of samples is small (less than twice the number of reference samples). " + "Result will have high variance and/or be non-discriminating." ) + + # Z-score normalization if z_score_norm: mean, std = _mean_std(x_samples, y_samples) - y_samples = (y_samples - mean) / std - x_samples = (x_samples - mean) / std - + if is_numpy: + x_samples = (x_samples - mean) / std + y_samples = (y_samples - mean) / std + elif is_torch: + x_samples = (x_samples - mean) / std + y_samples = (y_samples - mean) / std + # Determine fraction of x_samples to use as reference samples - x_frac = nx / (nx + ny) if x_frac is None else x_frac - + if x_frac is None: + x_frac = nx / (nx + ny) + # Determine number of samples from each distribution - probs = torch.tensor( - [ - x_frac * (1.0 - gauss_frac), - (1.0 - x_frac) * (1.0 - gauss_frac), - gauss_frac, - ], - device=device, - ) - - counts = Multinomial(total_count=num_refs, probs=probs).sample() - counts = counts.round().long() - Nx, Ny, Ng = counts.tolist() - assert (Nx + Ny + Ng) == num_refs, ( - f"Something went wrong. Nx={Nx}, Ny={Ny}, Ng={Ng} should sum to num_refs={num_refs}" - ) - - # Collect reference samples from x_samples - x_indices = torch.randperm(nx, device=device) - if Nx > nx: - raise ValueError("Cannot sample more references from x_samples than available") - xrefs_indices = x_indices[:Nx] - x_samples_indices = x_indices[Nx:] - - xrefs = x_samples[xrefs_indices] - x_samples = x_samples[x_samples_indices] - - # Collect reference samples from y_samples - y_indices = torch.randperm(ny, device=device) - if Ny > ny: - raise ValueError("Cannot sample more references from y_samples than available") - yrefs_indices = y_indices[:Ny] - y_samples_indices = y_indices[Ny:] - - yrefs = y_samples[yrefs_indices] - y_samples = y_samples[y_samples_indices] - - # Join the full set of reference samples - refs = torch.cat([xrefs, yrefs], dim=0) - - # Get gaussian reference points if requested - if Ng > 0: - m, s = _mean_std(x_samples, y_samples) - gauss_refs = torch.normal( - mean=m.repeat(Ng, 1), - std=s.repeat(Ng, 1), + if is_numpy: + counts = np.random.multinomial( + num_refs, + [x_frac * (1.0 - gauss_frac), (1.0 - x_frac) * (1.0 - gauss_frac), gauss_frac], ) - refs = torch.cat([refs, gauss_refs], dim=0) - - num_refs = refs.shape[0] - - # Compute nearest reference for x_samples - distances = torch.cdist(x_samples, refs) - idx = distances.argmin(dim=1) - counts_x = torch.bincount(idx, minlength=num_refs) - - # Compute nearest reference for y_samples - distances = torch.cdist(y_samples, refs) - idx = distances.argmin(dim=1) - counts_y = torch.bincount(idx, minlength=num_refs) - - # Remove reference samples with no counts - C = (counts_x > 0) | (counts_y > 0) - counts_x = counts_x[C] - counts_y = counts_y[C] - - n_filled_bins = C.sum().item() - if n_filled_bins == 1: - raise ValueError( - """ - Only one Voronoi cell has samples, so chi^2 cannot - be computed. This is likely due to a small number - of samples or a pathological distribution. If possible, - increase the number of x_samples and y_samples. - """ + Nx, Ny, Ng = counts + elif is_torch: + probs = torch.tensor( + [ + x_frac * (1.0 - gauss_frac), + (1.0 - x_frac) * (1.0 - gauss_frac), + gauss_frac, + ], + device=device, ) - if n_filled_bins < (num_refs // 2): - warnings.warn( - """ - Less than half of the Voronoi cells have any samples in them. - Possibly due to a small number of samples or a pathological - distribution. Result may be unreliable. If possible, increase the - number of x_samples and y_samples. - """ + counts_tensor = Multinomial(total_count=num_refs, probs=probs).sample() + counts = counts_tensor.round().long().cpu().numpy() + Nx, Ny, Ng = counts.tolist() + + # Validate counts + if Nx + Ny + Ng != num_refs: + raise ValueError( + f"Something went wrong. Nx={Nx}, Ny={Ny}, Ng={Ng} should sum to num_refs={num_refs}" ) - - # Perform chi-squared test - counts = torch.stack([counts_x, counts_y]) - return _chi2_contingency(counts, device) - + + # Sampling reference indices + if is_numpy: + refs = _sample_reference_indices_numpy(Nx, nx, Ny, ny, Ng, x_samples, y_samples) + elif is_torch: + refs = _sample_reference_indices_torch(Nx, nx, Ny, ny, Ng, x_samples, y_samples, device) + + # Update num_refs in case Gaussian samples were added + current_num_refs = refs.shape[0] + + # Compute nearest references and counts + if is_numpy: + return _compute_distances_numpy(x_samples, y_samples, refs, current_num_refs, num_refs) + + elif is_torch: + return _compute_distances_torch(x_samples, y_samples, refs, current_num_refs, num_refs) def pqm_pvalue( x_samples, @@ -320,7 +510,7 @@ def pqm_pvalue( x_samples, y_samples, and/or a Gaussian distribution, see the note below. re_tessellation : Optional[int] - Number of times pqm_pvalue is called, re-tesselating the space. No + Number of times _pqm_test is called, re-tesselating the space. No re_tessellation if None (default). z_score_norm : bool If True, z_score_norm the samples by subtracting the mean and dividing by the @@ -337,6 +527,8 @@ def pqm_pvalue( determined from the combined x_samples/y_samples. This ensures full support of the reference samples if pathological behavior is expected. Default: 0.0 no gaussian samples. + device : str + Device to use for computation. Default: 'cpu'. Note ---- @@ -350,8 +542,7 @@ def pqm_pvalue( ``Ny=(1-x_frac)*(1-gauss_frac)``, and ``Ng=gauss_frac`` respectively. For best results, we suggest using a large number of re-tessellations, though this is our recommendation in any case. - device : str - Device to use for computation. Default: 'cpu'. If 'cuda' is selected, + Returns ------- @@ -359,9 +550,20 @@ def pqm_pvalue( pvalue(s). Null hypothesis that both samples are drawn from the same distribution. """ - # Move samples to torch tensors on the selected device - x_samples = torch.tensor(x_samples, device=device) - y_samples = torch.tensor(y_samples, device=device) + # Check the device and convert to the respective type (Numpy or Torch) and call their respective _pqm_test function + + if device.type == 'cpu': + # Check if x_samples and y_samples are not already NumPy arrays + if not isinstance(x_samples, np.ndarray): + x_samples = x_samples.cpu().numpy() + if not isinstance(y_samples, np.ndarray): + y_samples = y_samples.cpu().numpy() + elif device.type == 'cuda': + # Check if x_samples and y_samples are not already torch tensors + if not torch.is_tensor(x_samples): + x_samples = torch.tensor(x_samples, device=device) + if not torch.is_tensor(y_samples): + y_samples = torch.tensor(y_samples, device=device) if re_tessellation is not None: return [ @@ -376,14 +578,14 @@ def pqm_pvalue( ) for _ in range(re_tessellation) ] - chi2_stat, p_value, dof, _ = _pqm_test( + + _, p_value, _, _ = _pqm_test( x_samples, y_samples, num_refs, z_score_norm, x_frac, gauss_frac, device ) # Return p-value as a float return p_value if isinstance(p_value, float) else float(p_value) - def pqm_chi2( x_samples, y_samples, @@ -401,10 +603,10 @@ def pqm_chi2( Parameters ---------- - x_samples : np.ndarray + x_samples : np.ndarray or torch.Tensor Samples from the first distribution. Must have shape (N, *D) N is the number of x samples, and D is the dimensionality of the samples. - y_samples : np.ndarray + y_samples : np.ndarray or torch.Tensor Samples from the second distribution. Must have shape (M, *D) M is the number of y samples, and D is the dimensionality of the samples. num_refs : int @@ -412,7 +614,7 @@ def pqm_chi2( x_samples, y_samples, and/or a Gaussian distribution, see the note below. re_tessellation : Optional[int] - Number of times pqm_chi2 is called, re-tesselating the space. No + Number of times _pqm_test is called, re-tesselating the space. No re_tessellation if None (default). z_score_norm : bool If True, z_score_norm the samples by subtracting the mean and dividing by the @@ -430,7 +632,7 @@ def pqm_chi2( support of the reference samples if pathological behavior is expected. Default: 0.0 no gaussian samples. device : str - Device to use for computation. Default: 'cpu'. If 'cuda' is selected, + Device to use for computation. Default: 'cpu'. Note ---- @@ -461,9 +663,20 @@ def pqm_chi2( float or list chi2 statistic(s). """ - # Move samples to torch tensors on the selected device - x_samples = torch.tensor(x_samples, device=device) - y_samples = torch.tensor(y_samples, device=device) + + # Check the device and convert to the respective type (Numpy or Torch) and call their respective _pqm_test function + if device.type == 'cpu': + # Check if x_samples and y_samples are not already NumPy arrays + if not isinstance(x_samples, np.ndarray): + x_samples = x_samples.cpu().numpy() + if not isinstance(y_samples, np.ndarray): + y_samples = y_samples.cpu().numpy() + elif device.type == 'cuda': + # Check if x_samples and y_samples are not already torch tensors + if not torch.is_tensor(x_samples): + x_samples = torch.tensor(x_samples, device=device) + if not torch.is_tensor(y_samples): + y_samples = torch.tensor(y_samples, device=device) if re_tessellation is not None: return [ @@ -479,6 +692,7 @@ def pqm_chi2( for _ in range(re_tessellation) ] + chi2_stat, _, dof, _ = _pqm_test( x_samples, y_samples, num_refs, z_score_norm, x_frac, gauss_frac, device ) From 2642cc63f654e52a0c1c32e4b9ee2b2752cb46d5 Mon Sep 17 00:00:00 2001 From: Sammy Sharief Date: Sun, 3 Nov 2024 19:30:45 -0500 Subject: [PATCH 4/7] Updated test_guassian --- README.md | 16 +++++++++++----- src/pqm/test_gaussian.py | 4 ++-- tests/test_gaussian.py | 2 +- 3 files changed, 14 insertions(+), 8 deletions(-) diff --git a/README.md b/README.md index 06c7503..44b7427 100644 --- a/README.md +++ b/README.md @@ -6,15 +6,13 @@ ![PyPI - Downloads](https://img.shields.io/pypi/dm/pqm) [![arXiv](https://img.shields.io/badge/arXiv-2402.04355-b31b1b.svg)](https://arxiv.org/abs/2402.04355) - - -[PQMass](https://arxiv.org/abs/2402.04355) is a new sample-based method for evaluating the quality of generative models as well assessing distribution shifts. +[PQMass](https://arxiv.org/abs/2402.04355) is a new sample-based method for evaluating the quality of generative models as well as assessing distribution shifts to determine if two datasets come from the same underlying distribution. ## Install To install PQMass, run the following: -```python +```bash pip install pqm ``` @@ -50,7 +48,7 @@ PQMass can work for any two datasets as it measures the distribution shift betwe ## Example -We are using 100 regions. Thus, the DoF is 99, our expected $\chi^2$ peak of the distribution is 97, the median is 99, and the standard deviation should be 14.07. With this in mind, we set up our example. For the p-value, we expect to be between 0 and 1 and a significantly small p-value (e.g., $< 0.05$ or $< 0.01$) would mean we reject the null hypothesis and thus $x$ and $y$ do not come from the same distribution. +We are using 100 regions. Thus, the DoF is 99, our expected $\chi^2$ peak of the distribution is 97, the mean is 99, and the standard deviation should be 14.07. With this in mind, we set up our example. For the p-value, we expect to be between 0 and 1 and a significantly small p-value (e.g., $< 0.05$ or $< 0.01$) would mean we reject the null hypothesis and thus $x$ and $y$ do not come from the same distribution. Our expected p-value should be around 0.5 to pass the null hypothesis test; any significant deviation away from this would indicate failure of the null hypothesis test. @@ -96,6 +94,14 @@ Here it is clear that both $\chi_{PQM}^2$ and $\text{p-value}(\chi_{PQM}^2)$ are Thus, PQMass can be used to identify if any two distributions come from the same underlying distributions if enough samples are given. We encourage users to look through the paper to see the varying experiments and use cases for PQMass! +## How to Intrept Result + +We have shown what to expect for PQMass when working with $\chi_{PQM}^2$ or $\text{p-value}(\chi_{PQM}^2)$ however when working with $\chi_{PQM}^2$, there is the case in which it will return 0's. There are a couple reasons in why this could happen + +- For Generative Models; 0's indicate memorization. Samples are duplicates of the data it has been trained on. +- For non generative model scenario, it is typically due to lack of samples espically in high dimensions. Increasing samples should alleviate the issue. +- Another scenario in which one could get 0's in a non generative model case is that it can also be an inidcator of duplicate samples in $x$ and $y$. + ## Advanced Usage Depending on the data you are working with we show other uses of the parameters for PQMass. diff --git a/src/pqm/test_gaussian.py b/src/pqm/test_gaussian.py index f99cd6d..db9f7de 100644 --- a/src/pqm/test_gaussian.py +++ b/src/pqm/test_gaussian.py @@ -8,6 +8,6 @@ def test(): y_samples = np.random.normal(size=(500, 50)) x_samples = np.random.normal(size=(250, 50)) - new.append(pqm_pvalue(x_samples, y_samples, device = 'cpu')) + new.append(pqm_pvalue(x_samples, y_samples)) - assert np.abs(np.mean(new) - 0.5) < 0.15 + assert np.abs(np.mean(new) - 0.5) < 0.15 \ No newline at end of file diff --git a/tests/test_gaussian.py b/tests/test_gaussian.py index 0006746..170900f 100644 --- a/tests/test_gaussian.py +++ b/tests/test_gaussian.py @@ -68,4 +68,4 @@ def test_fracs(x_frac, gauss_frac): x_samples[:, 0] += 5 # one dim off by 5sigma pval = pqm_pvalue(x_samples, y_samples, x_frac=x_frac, gauss_frac=gauss_frac) - assert pval < 1e-3 + assert pval < 1e-3 \ No newline at end of file From cd8344d41b0ce692384b4009968d2f608a55607e Mon Sep 17 00:00:00 2001 From: Sammy Sharief Date: Sun, 3 Nov 2024 19:35:51 -0500 Subject: [PATCH 5/7] Updated the default device from 'cpu' to torch.device('cpu') --- src/pqm/pqm.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/pqm/pqm.py b/src/pqm/pqm.py index 8b40432..3c02a3c 100644 --- a/src/pqm/pqm.py +++ b/src/pqm/pqm.py @@ -348,7 +348,7 @@ def _pqm_test( z_score_norm: bool, x_frac: Optional[float], gauss_frac: float, - device: str = 'cpu', + device: str = torch.device("cpu"), ) -> Tuple: """ Helper function to perform the PQM test and return the results from @@ -490,7 +490,7 @@ def pqm_pvalue( z_score_norm: bool = False, x_frac: Optional[float] = None, gauss_frac: float = 0.0, - device: str = 'cpu', + device: str = torch.device("cpu"), ): """ Perform the PQM test of the null hypothesis that `x_samples` and `y_samples` @@ -594,7 +594,7 @@ def pqm_chi2( z_score_norm: bool = False, x_frac: Optional[float] = None, gauss_frac: float = 0.0, - device: str = 'cpu', + device: str = torch.device("cpu"), ): """ Perform the PQM test of the null hypothesis that `x_samples` and `y_samples` From 34879274f487b0a2099bdf77f332f0ca41295fc4 Mon Sep 17 00:00:00 2001 From: Sammy Sharief Date: Sun, 3 Nov 2024 19:46:45 -0500 Subject: [PATCH 6/7] Updated _mean_std function and passed all test locally --- src/pqm/pqm.py | 56 +++++++++++++++++++++++++++++++++++--------------- 1 file changed, 40 insertions(+), 16 deletions(-) diff --git a/src/pqm/pqm.py b/src/pqm/pqm.py index 3c02a3c..224a3c8 100644 --- a/src/pqm/pqm.py +++ b/src/pqm/pqm.py @@ -11,23 +11,47 @@ def _mean_std(sample1, sample2, dim=0): """Get the mean and std of two combined samples without actually combining them.""" - n1 = sample1.shape[dim] - n2 = sample2.shape[dim] - # Get mean/std of combined sample - mx = torch.mean(sample1, dim=dim) - sx = torch.std(sample1, dim=dim, unbiased=True) - my = torch.mean(sample2, dim=dim) - sy = torch.std(sample2, dim=dim, unbiased=True) - m = (n1 * mx + n2 * my) / (n1 + n2) - s = torch.sqrt( - ( - (n1 - 1) * (sx ** 2) - + (n2 - 1) * (sy ** 2) - + n1 * n2 * (mx - my) ** 2 / (n1 + n2) + # Check if both samples are PyTorch tensors + if isinstance(sample1, torch.Tensor) and isinstance(sample2, torch.Tensor): + n1 = sample1.shape[dim] + n2 = sample2.shape[dim] + + mx = torch.mean(sample1, dim=dim) + sx = torch.std(sample1, dim=dim, unbiased=True) + my = torch.mean(sample2, dim=dim) + sy = torch.std(sample2, dim=dim, unbiased=True) + + m = (n1 * mx + n2 * my) / (n1 + n2) + s = torch.sqrt( + ( + (n1 - 1) * (sx ** 2) + + (n2 - 1) * (sy ** 2) + + n1 * n2 * (mx - my) ** 2 / (n1 + n2) + ) + / (n1 + n2 - 1) ) - / (n1 + n2 - 1) - ) - return m, s + return m, s + + # Check if both samples are NumPy arrays + elif isinstance(sample1, np.ndarray) and isinstance(sample2, np.ndarray): + n1 = sample1.shape[dim] + n2 = sample2.shape[dim] + + mx = np.mean(sample1, axis=dim) + sx = np.std(sample1, axis=dim, ddof=1) + my = np.mean(sample2, axis=dim) + sy = np.std(sample2, axis=dim, ddof=1) + + m = (n1 * mx + n2 * my) / (n1 + n2) + s = np.sqrt( + ( + (n1 - 1) * (sx ** 2) + + (n2 - 1) * (sy ** 2) + + n1 * n2 * (mx - my) ** 2 / (n1 + n2) + ) + / (n1 + n2 - 1) + ) + return m, s def rescale_chi2(chi2_stat, orig_dof, target_dof, device): """ From 34e2ae64bfabe2e9ff71fa75ad69cdeb26f44c38 Mon Sep 17 00:00:00 2001 From: Sammy Sharief Date: Mon, 4 Nov 2024 14:39:50 -0500 Subject: [PATCH 7/7] Updated repo post comments --- README.md | 3 +-- src/pqm/pqm.py | 30 ++++++++++++++---------------- 2 files changed, 15 insertions(+), 18 deletions(-) diff --git a/README.md b/README.md index 44b7427..3e26801 100644 --- a/README.md +++ b/README.md @@ -34,8 +34,7 @@ peak of this distribution will be at `DoF - 2`, the mean will equal `DoF`, and the standard deviation will be `sqrt(2 * DoF)`. If your $\chi_{PQM}^2$ values are too high (`chi^2 / DoF > 1`), it suggests that the samples are out of distribution. Conversely, if the values are too low (`chi^2 / DoF < 1`), it indicates -potential duplication of samples between `x` and `y` (i.e. -memorization for generative models). +potential duplication of samples between `x` and `y`. If your two samples are drawn from the same distribution, then the $\text{p-value}(\chi_{PQM}^2)$ should be drawn from the random $\mathcal{U}(0,1)$ distribution. This means that if diff --git a/src/pqm/pqm.py b/src/pqm/pqm.py index 224a3c8..5d441a7 100644 --- a/src/pqm/pqm.py +++ b/src/pqm/pqm.py @@ -57,7 +57,7 @@ def rescale_chi2(chi2_stat, orig_dof, target_dof, device): """ Rescale chi2 statistic using appropriate methods depending on the device. """ - if device.type == 'cuda': + if device.type == 'cuda' or device == 'cuda': # Move tensors to CPU and convert to NumPy chi2_stat_cpu = chi2_stat.cpu().item() # Convert to float orig_dof_cpu = orig_dof.cpu().item() # Convert to float @@ -329,11 +329,12 @@ def _compute_distances_torch(x_samples, y_samples, refs, current_num_refs, num_r # Compute distances and find nearest references distances_x = torch.cdist(x_samples, refs) - idx_x = distances_x.argmin(dim=1) - counts_x = torch.bincount(idx_x, minlength=current_num_refs) - distances_y = torch.cdist(y_samples, refs) + + idx_x = distances_x.argmin(dim=1) idx_y = distances_y.argmin(dim=1) + + counts_x = torch.bincount(idx_x, minlength=current_num_refs) counts_y = torch.bincount(idx_y, minlength=current_num_refs) # Remove references with no counts @@ -391,7 +392,7 @@ def _pqm_test( x_samples, y_samples, and/or a Gaussian distribution, see the note below. re_tessellation : Optional[int] - Number of times pqm_pvalue is called, re-tesselating the space. No + Number of times _pqm_test is called, re-tesselating the space. No re_tessellation if None (default). z_score_norm : bool If True, z_score_norm the samples by subtracting the mean and dividing by the @@ -453,12 +454,9 @@ def _pqm_test( # Z-score normalization if z_score_norm: mean, std = _mean_std(x_samples, y_samples) - if is_numpy: - x_samples = (x_samples - mean) / std - y_samples = (y_samples - mean) / std - elif is_torch: - x_samples = (x_samples - mean) / std - y_samples = (y_samples - mean) / std + + x_samples = (x_samples - mean) / std + y_samples = (y_samples - mean) / std # Determine fraction of x_samples to use as reference samples if x_frac is None: @@ -497,7 +495,7 @@ def _pqm_test( refs = _sample_reference_indices_torch(Nx, nx, Ny, ny, Ng, x_samples, y_samples, device) # Update num_refs in case Gaussian samples were added - current_num_refs = refs.shape[0] + current_num_refs = Nx + Ny + Ng # Compute nearest references and counts if is_numpy: @@ -576,13 +574,13 @@ def pqm_pvalue( """ # Check the device and convert to the respective type (Numpy or Torch) and call their respective _pqm_test function - if device.type == 'cpu': + if device.type == 'cpu' or device == 'cpu': # Check if x_samples and y_samples are not already NumPy arrays if not isinstance(x_samples, np.ndarray): x_samples = x_samples.cpu().numpy() if not isinstance(y_samples, np.ndarray): y_samples = y_samples.cpu().numpy() - elif device.type == 'cuda': + elif device.type == 'cuda' or device == 'cuda': # Check if x_samples and y_samples are not already torch tensors if not torch.is_tensor(x_samples): x_samples = torch.tensor(x_samples, device=device) @@ -689,13 +687,13 @@ def pqm_chi2( """ # Check the device and convert to the respective type (Numpy or Torch) and call their respective _pqm_test function - if device.type == 'cpu': + if device.type == 'cpu' or device == 'cpu': # Check if x_samples and y_samples are not already NumPy arrays if not isinstance(x_samples, np.ndarray): x_samples = x_samples.cpu().numpy() if not isinstance(y_samples, np.ndarray): y_samples = y_samples.cpu().numpy() - elif device.type == 'cuda': + elif device.type == 'cuda' or device == 'cuda': # Check if x_samples and y_samples are not already torch tensors if not torch.is_tensor(x_samples): x_samples = torch.tensor(x_samples, device=device)