From 40e183aa05a3fe41100e835fcdd603475927321e Mon Sep 17 00:00:00 2001 From: Leonardo Leone Date: Tue, 26 May 2026 15:42:56 +0200 Subject: [PATCH 1/8] Put IsInConvex function to DepthEucl --- .gitignore | 1 + depth/model/DepthEucl.py | 20 +++++++ .../model/multivariate/Depth_approximation.py | 2 +- depth/model/multivariate/__init__.py | 4 +- depth/model/multivariate/isInConvexes.py | 55 ++++++++++++++++++ depth/src/ddalpha.cpp | 52 +++++++++++++++-- depth/src/ddalpha.dll | Bin 255139 -> 255383 bytes 7 files changed, 128 insertions(+), 6 deletions(-) create mode 100644 depth/model/multivariate/isInConvexes.py diff --git a/.gitignore b/.gitignore index edc0c83..3814d42 100644 --- a/.gitignore +++ b/.gitignore @@ -4,6 +4,7 @@ __pycache__/ # C extensions *.so +*.dll auto_examples diff --git a/depth/model/DepthEucl.py b/depth/model/DepthEucl.py index 139016d..3aacb39 100644 --- a/depth/model/DepthEucl.py +++ b/depth/model/DepthEucl.py @@ -1316,6 +1316,26 @@ def ACA(self,dim:int=2, line_solver=line_solver, bound_gc=bound_gc) return ACA_tab + + def IsInConvexes(self,z): + """ + Checks the belonging to at least one of class convex hulls of the training dataset + + Parameters + ---------- + z: array_like, + points to check the belongingness + + + Results + -------- + belongingness : array_like + The return respresents directions that best represents anomalies in the dataset. + + """ + return mtv.IsInConvexes(self.data,z,self.distribution,self.seed) + + ## Det and MCD def _calcDet(self,mat:np.ndarray): diff --git a/depth/model/multivariate/Depth_approximation.py b/depth/model/multivariate/Depth_approximation.py index e73dab4..41f2c0a 100644 --- a/depth/model/multivariate/Depth_approximation.py +++ b/depth/model/multivariate/Depth_approximation.py @@ -3,7 +3,7 @@ from math import ceil import sys, os, glob import platform -from .import_CDLL import libApprox +from .import_CDLL import libApprox, libExact def depth_approximation(z, X, diff --git a/depth/model/multivariate/__init__.py b/depth/model/multivariate/__init__.py index 264287d..18105ad 100644 --- a/depth/model/multivariate/__init__.py +++ b/depth/model/multivariate/__init__.py @@ -19,6 +19,8 @@ try:from .CUDA_approximation import cudaApprox except:cudaApprox=None from .ACA_wrapper import ACA +from .import_CDLL import libExact,libApprox +from .isInConvexes import IsInConvexes __all__ = ["depth_approximation", "betaSkeleton", "cexpchull", "cexpchullstar", "geometrical", "halfspace", "L2", "mahalanobis", "potential", "projection", "aprojection", "qhpeeling", "simplicial", "simplicialVolume", "spatial", "zonoid", "depth_mesh", "depth_plot2d", "calcDet", - "MCD", "cudaApprox","ACA"] + "MCD", "cudaApprox","ACA", "IsInConvexes"] diff --git a/depth/model/multivariate/isInConvexes.py b/depth/model/multivariate/isInConvexes.py new file mode 100644 index 0000000..0b1d50b --- /dev/null +++ b/depth/model/multivariate/isInConvexes.py @@ -0,0 +1,55 @@ +import numpy as np +from ctypes import * +from math import ceil +import sys, os, glob +import platform +from .import_CDLL import libExact + +def IsInConvexes(X,z,distributions,seed): + """ + Check if points are inside the convex hull + """ + try: + n, d = X.shape + except ValueError: + n = X.shape[0] + d = 1 + n_z = z.shape[0] + points_list=X.flatten() + points=(c_double*len(points_list))(*points_list) + objects_list=z.flatten() + objects=(c_double*len(objects_list))(*objects_list) + distrSeq_list=distributions.flatten() + distrSeq=(c_int*len(distrSeq_list))(*distrSeq_list) + points=pointer(points) + objects=pointer(objects) + distrSeq=pointer(distrSeq) + + distr=np.unique(distributions,return_counts=True)[1] + distribution_list=distr.flatten() + distribution=(c_int*len(distribution_list))(*distribution_list) + distribution=pointer(distribution) + + CSum=np.zeros(distr.shape,dtype=int) + CSum[1:]=distr.cumsum(dtype=int)[:-1] + cumSum_list=CSum.flatten() + cumSum=(c_int*len(cumSum_list))(*cumSum_list) + cumSum=pointer(cumSum) + numPoints=pointer(c_int(n)) + numObjects=pointer(c_int(n_z)) + dimension=pointer(c_int(d)) + seed=pointer((c_int(seed))) + numClasses=pointer(c_int(distr.shape[0])) + belongs=pointer((c_int*len(z))(*np.zeros(distr.shape[0],dtype=int))) + libExact.IsInConvexes(points,dimension,distribution,numClasses, + objects,numObjects,seed,belongs,cumSum, distrSeq,) + res=np.zeros((distr.shape[0],len(z))) + for i in range(distr.shape[0]): + for j in range(len(z)): + res[i][j]=belongs[i][j] + + + return res + + + \ No newline at end of file diff --git a/depth/src/ddalpha.cpp b/depth/src/ddalpha.cpp index 7e271f7..878c723 100644 --- a/depth/src/ddalpha.cpp +++ b/depth/src/ddalpha.cpp @@ -38,17 +38,25 @@ void setSeed(int random_seed){ } } -void IsInConvexes(double *points, int *dimension, int *cardinalities, int *numClasses, double *objects, int *numObjects, int *seed, int *isInConvexes){ +void IsInConvexes(double *points, int *dimension, int *cardinalities, + int *numClasses, double *objects, int *numObjects, int *seed, int *isInConvexes, + int *cumSum, int *distrSeq){ setSeed(*seed); - int numPoints = 0;for (int i = 0; i < numClasses[0]; i++){numPoints += cardinalities[i];} + int numPoints = 0; + for (int i = 0; i < numClasses[0]; i++){numPoints += cardinalities[i];} TMatrix x(numPoints); + TPoint disCount(numClasses[0]); + int posX; for (int i = 0; i < numPoints; i++){x[i] = TPoint(dimension[0]);} for (int i = 0; i < numPoints; i++){ + posX=cumSum[distrSeq[i]] + disCount[distrSeq[i]]; for (int j = 0; j < dimension[0]; j++){ - x[i][j] = points[i * dimension[0] + j]; + x[posX][j] = points[i * dimension[0] + j]; } + disCount[distrSeq[i]] += 1; } + TMatrix o(numObjects[0]); for (int i = 0; i < numObjects[0]; i++){o[i] = TPoint(dimension[0]);} for (int i = 0; i < numObjects[0]; i++){ @@ -69,6 +77,42 @@ void IsInConvexes(double *points, int *dimension, int *cardinalities, int *numCl } } +// void IsInConvexes(double *points, int *dimension, int *cardinalities, +// int *numClasses, double *objects, int *numObjects, int *seed, int *isInConvexes, +// int *cumSum, int *distrSeq){ +// setSeed(*seed); +// for(int i = 0; i < 2; i++){ +// cumSum[i]; +// } +// int numPoints = 0;for (int i = 0; i < numClasses[0]; i++){numPoints += cardinalities[i];} +// TMatrix x(numPoints); +// for (int i = 0; i < numPoints; i++){x[i] = TPoint(dimension[0]);} +// for (int i = 0; i < numPoints; i++){ +// for (int j = 0; j < dimension[0]; j++){ +// x[i][j] = points[i * dimension[0] + j]; +// } +// } +// TMatrix o(numObjects[0]); +// for (int i = 0; i < numObjects[0]; i++){o[i] = TPoint(dimension[0]);} +// for (int i = 0; i < numObjects[0]; i++){ +// for (int j = 0; j < dimension[0]; j++){ +// o[i][j] = objects[i * dimension[0] + j]; +// } +// } +// TVariables cars(numClasses[0]); +// for (int i = 0; i < numClasses[0]; i++){ +// cars[i] = cardinalities[i]; +// } +// TIntMatrix answers(o.size()); +// int error = 0; +// InConvexes(x, cars, o, error, &answers); +// for (int i = 0; i < numObjects[0]; i++) +// for (int j = 0; j < numClasses[0]; j++){ +// isInConvexes[numClasses[0]*i+j] = answers[i][j]; +// } +// } + + void ZDepth(double *points, double *objects, int *numPoints, int *numObjects, int *dimension, int *seed, double *depths){ setSeed(*seed); TMatrix x(numPoints[0]); @@ -92,7 +136,7 @@ void ZDepth(double *points, double *objects, int *numPoints, int *numObjects, in depths[i] = ZonoidDepth(x, z[i], error); } } - + void HDepthSpaceEx(double *points, double *objects, int *cardinalities, int *numClasses, int *numObjects, diff --git a/depth/src/ddalpha.dll b/depth/src/ddalpha.dll index c197ec284390e563e10236fd3cbd090b2ce89396..c759f5286019c7eb93b9970a1a07b772de5d871e 100644 GIT binary patch delta 36632 zcmce86*UE~t z%&j6RM@qQiQhNTp3<*ymsUVo1%W?P$CqQ7fcvn5E_TyySoQ{H%{ zC}z9Y;xECr4hsoW2aXl_uASZdobo@5)2?At@4id&dBvu@%Zj3ux+0IF#>FFx-Rt^k zrl=Q2GbOKh;S96Q*F($b^&7IN9}d*WLqI-`HrR}Z zg_~P&gCwRTTJKxGHe2tghySUninEr0XA`~MVpx(`>ejT8QzGgtNsNRzI8k2#)gQGeI)BB66?0Z_p^0V*cfrP?n(Am(cATUF*ZsM3C z*@*6JZCYi(wY7e31ri!IH~VsCPaW<{LIC;k+A z$=f)qsD0ay7;7M|whQxZz?7)s8yDkps;V*%R#-^?&we5xGMV{_jL1&--W-{P@82Uk zD6yt}0}TD<_T4NH$)7s4ttlXvB$T)j zal_+=4R_A>K(npWfW=HnupY`gYg+Rf_~x{;iD@e69Y8l*zqHs}pD~yGXt3Bj&Fd`M zbO>%R85-ku(|R@Aeczr>QYC6tJu$6Alt*juDA}aG*w`UF_|JM;5;YDSGk(IjxUq5L z;>K582GtY*xl3w!J#o83kk5XSvl;e!XiB}7cM^de{dkLdBD&)-R#)&&@vO7x*D0Fa z64{-S*@mJsotAOdtY}bla~9;q6lZO64o!btO5B8lG`*tLPrSrSF+G@)MO-c6DLQQ+(BJEh`mMyH96ti7&f<+b??tQ!4(# zx+ay8SVVW#Ja>>n9_l1@5t2~XX$~ZF{$Zv(cVgLF-$IL&3Kdz^dhX)49tV24QI{OF zSXFb0+TE1*G5Cqpzt#$Mtj2?xDUB*N1;j$GeyKU9m|58Mm?K$tpufjp@q5oc_5PTK z(&9}iM`FIs+H1D{sXnMB`uB=rwZzh11Mz*L*GlFphQ{`(4Z8}y*?L%A?J73L_QG>})Xe}n8SGg#Q4{1}zDG$xQ$T1b4 zA)yT+E#S|ILrYm3Ka1VCEJOq)G!Jsu6dBAmvub9SY=Z-Hjx)7Bge&IJ^7xC1$qDVE z_Dn(U@G^7BeZ!58Xku-&|4(zyL8jIsy>{RJsWjins?*Hc$i$(9E)mCAb!}*7RqRA^ zYChF=H!&VhrDuahXx|WCkBI?&d$=E(jJ0fU;}tB{^-b`dxR z8&%GV+p0tdbA-n^QB=MC%dt(}&1q!>qW|+cqrBMg8%i zC^7qWw?88P+S}x933F3VfbvK|%f|}R9!Zh{Z}p9bqF!PcU-wY-NNnems&QSb7hF)F zSeDq7ANWUXNL(UJa)v-iYN%Rv?{Rwqqy{Zmmc8b25Jl9spa zGK%(#Va8ulu%|{~_FAlyLz1zB+^?wmyQ(T7uQ;o&Mxv>^QEVC*%Z`e(!BKq!e9Rf z3gb*;JGPw~X9_!t7If+dI<2U3bIxHmG~kH3@E;K~s$)2bqNhW`K&&M3Lu%{Cw0Kh( z?WU!w|3Bi@QFB;qk#Tf8##}_FF6lL(d~)NIs3cFXYz}j z{V6s~U72Oe4GDx%Pq1e&*v4p>{&?$QGpJ#pCZ*jbSI(S3)K5ZI(R_te>?1Q4<)SG4 z?X@*?H_#q(>prdUUEgV)*Xc+v_K4Zv-thzE*0P6LSULQVNK(}&MO_7P%nC|Hs?Q-2 zZ~Zl%(g)T3$EvE5(!RsWNQ8A=)X|)F7|A_k6fE{wo@5`&EY^$SHqy2V;aVzEiq z+ZLoVQ2PYygG6ho#X2@*F*cE$$swsq-mR?06_Lne3K*45McO=aRA$QUL~3%byCm3M zT;gL&GF5vr%PC>z_14-sF&CD$G^ZWdL-hoaTp9NK$g4lnyP>R6J^ z#UzH!=wvhxu(@;*Yi7h}kvo#Cnf)ZHqqlkthJlDFG^Fup7tIn22$2>TbeikxuH#iz zxy1|P;$wm{r;%*5);OD+Rn5C)TKki>oq~!wTHW$)E!+s#O8~nPB~ys)T;LX|bUP zldOLvT2CkC9W$*eLUl0=#lI$+B1)1lcUEOPDaVL%6s|Sq1KRACm~+tW*JSvCeePRg zHiVFrDeqO?I|x(9d_%Q@pP&vV40}={V4|&#Iqz7etJ)oje3SPlA^DWgdtXQaAnd%D zIII42Z!j&^0~TuuYI4kc=6cTe%^~oQR7XycOG@A!YW21q0Yv6D_^uuw|KyqDCWG&raV==klvmvm&LC1*G1}V z3-|w8d^kI}jyLJDdo_uH&1~?s_zs-cE{QgCnsJ}IMM-mfIsf&yqRAOBc$dV>^DKPz zdU17Lt6D#Ig#6VDqW=6I$OA>W;edoXEBHKW9#3SHz%&U$ricf-A$cY}H=j$uscwVnL$% ze?feO`h4LTF?dnD*GO<-m9ctN%Nen8QBS_%mbkiTRHI{M#JEJB{S)T;Bsp`JM&)`SD!_G92)*>{!>N6mr(u_V%*Yqn~!;o=ElQVZPl9A(1E~g zpBHPfW6_$ezp7*27QxGwv&%&rm+fSY*B=Gj=7*8{R_nI3@8Az5!jr7JCgx@};8Tx^ zRaw#OrYOmp#t)T>i0tgJ(KtC%cMSo=B+sO^Ux8+spJ-OcQgJr>n%~`Ts;VlcQ2#u* zf^Iub*k6Z>PCowuX9tUxV4` zc>4?Tqp%{#p-XF~Ibv$pFT~ANgS_^9PR0C6OZmgkh1uG&ZXPJwf!bodO6`4fm3ZEI zV$7EPXvc6`HP$0HnqcWTnqVpH1XZ(hkOV{B{lJ+xpEV`EOAI>b9I9c;v=U zP?1S0=I$5oznH+4{o?kE;emM{Q3-#-&L+UT9}27P+=}!25JXwKPjs>E;O_gxkG2*3 z`On1c)tmV5dqrJ)u>Y5Pq1bLz+K~SEkid$0)alRe6-jn0_uDHj+Q)In9?@mZ>{h9J zw7NVscxIy8+dSI}O(lb^-rFs{T$9c}-Yt5p9nFXD7O$-x-f~SDRq>8mXE${n{P`!( zW)*iGsOiTJ;hS&c7KhlFKZXCeOSrEa$ye_Zlh!S)`{{nIjcVvFtr}|6U7})L<5~)m z>Sv!8dA>B9w>?>oI?x$R8=ODQ;Rn<`M7j-IjLj5g6eni%(%6tdHS(lb^YS3mGrh_A z@|_M~Yr3MhxcTx}9@bm*evq~r(rzW%>cWs1hj@!mWT?+j@TrYA zFzZpqs$Sx{Y}~6KsOtQS*z=pFI`wn{s9?~bpYp>jHe(4QKxoGS082Rbjr(TnIdhtC z+g>7meTzmOTE@mv$FQlK#Sw5m#C|g^_nK)Ex$ArK2$MLpKCR8E9ys*nw$Y}yu^5Sp zAmDTi3_pO0x)~Jp9rV8?&9_<4qG7MaF#gdsvG(FL{Pi))Ry}qAzzRwIvrA2P0ehx8;zcvnHKZyRDdWAms zF7j(ZxZj;YxM+1cNMuP0a})^E`w=Ll>-%YZP=R>+{Z77N3c^cKcQLW5a6bPtDea+ZEGVjY zpR;-m-==CDWN|T9vQ&zHxVhb5 zU2OeeNc%glp%9u573xnw`&(UZCt?LB#Qc(VgNSE{Z?BCVfaDIVC9jD#AKv59&xjTu zO$}bWo~lSa@Ld+Rv^(lFA|TrY4gRT)TrWQQD4w5@;*XC)xM7}XT)3$F7^&&hP9cf* z%s-HWPE7G2V@P^fw8M2N;@Y&qLo^N1rauvjDpiA|I8~U>-@hv+d>kC!^&>DbdePcW zrY?9@QEq%p-ZI8CH{>X?Qpo#H^{m+b@xZ2th-j0>9P=Bx6Xd^Mab}NHD_;?QJDam> zB4%eG-!xy0+1X^w6yn-_-BC-58sQB087Np!JDiNXzld>3ZG+?yadzh@Uh=YNQ`CmP^s*RR)QB&7Su8GU!N>nzR8Z8M@dzOfeUi<$ zyd?U5`h)M~b;t>y-5hg6p9t#Fb)x4kKVG;_jMx=i&-O?_Z`zV|V)d@Dp2?u!*lD)G zjfumcR-A@SQHRdNEId-tji%)>)t9o;Hj}e2q!_JWSzV%z!JMfN^2M*a;`o_-(aF)e z(PuNXOqz>hQZ@W@S}_Ci#dD6qv4J=XXwCdn^#zJq^iyKaA)vIeo?cM(T=vvWY_d{-;FIw!G+GX(?t^S*+iy z1Bikl?XyN4+!N3D+lBkyK=zx6*!xWP0e0~8GjEOA))8cyZQXBC*8e)&eC*=ey#;*b zO`2^HjmT&{o^3;M=>D&>ZRZqm^0TC-H*64xje0cQz9s@qPPNQN(~TdREV}Lt;5pO9 zbNiZ%Ng%e}x8gd0bH2R+3g+8tCu2$_F)pcINIJW;o|94dIL8c79?iF#FH(j(s^vBD z>%N}sK~dEHSjJBjiv^$0WWN+$`FtefaeyBKqzGxyvs z#DOS&aho`DAea|z6E_Y7@mL}192~^&<%$Uho78?5>l`5tI60-cV&%c-_>1ogmqSN+ z@n&)Qkbi?kIht8}{!no}MR8iTgX+4?;-5qD%^QP9z4Dtj0S}t(t3xW?u@*1C;G|#9 z7iovr2j6-AaWLkH+9AEF3XzzEyeiY$D4gAB8vHDp9SPQ*O<0*8yg6AzCz;0?YL4UTnXH^h)*{o58Sbc(NnV*_Pg#enK(m|)oT z?^mQfKyfw(hH~SQ`e2zjek_)S3D4urvLaom;gscSQ$8{=PV5g6Bmy+ zZ5FHESPB)|JpP_Abp9nk+TYGB74=R;41EvDc7VEoOOeN8J(rmB) zhqJ;jhhzh>u;bM3#R*Yo<8mG%G%^+bmBqRgfe(suh*>8i+$oHLQX1Y9@11OgW80aN zZ@8}3a*MgAqOwvfw#gLGdpQJq4;lu?!D8K@MJW*7AuBNr8;r%al+JWl=~zc`6x#i} zCt1n^2S$|u)8hU$ z@eXJamf5iU6lg$z+HAcSb=Qo!lxne_(?VU^IWHj!0U8{85j>=OjJRp2N6greu!Y|w zT7NU6v8iV3L$mc48?;nVkVT#PY!D92H3D!`^Pjk=OSKs*#mA>3+=Ad_+S+{KCVn{G z-1m~AMB&QU+qW68q_kTWd+200(df%!-gU9~=F5g1E`pp_+}m1T$!v`SJ}l7VBT=Mk zMPMDd$6`RKi%F5|MG;mOUjN)PG}2m#8ViCL$NnRxltuEj&j?u-#4^R6vaT&pe1&Cd zfhg^ev1SY(`sARt(*9CYb~z{ausqS|%_Xl)E$En1(+df;kPAyF8OM)};;*w!`9GPW;kgho zx0`1z-OUlb&W+}O=ZiPaJ;#^l3;*-Y`H+0k{e08V?F*=e+I<#IJ#a;isI=cE;T-`_ z930^rsp}Vr<>&kGFS$5)K7uz`AbvR?(IAKRSy~2Iqh@Pv$ZjkP^{e?J@}3kmSqGAT-wJv)#2A=>?lnz;(uJ4$}$Ha;uwD%+2dD8ZHOOxJI5q{j6VO zMU|RX48?&%9p4j{GVQCVVzY?;&kVjGLu~p_OoQkQ$mQM5^rlWgWA2#4e^rAs#Bcv; z)Bg9l;N={|S@g=nAIy2TGQCI$!&sqS(pZQQxKJkqnq%Ta7Mt?Mh#?o-anBw_ITt%J zo-{}7xfIH)?uqhCzF7ftNZAo6y9Z@(Aq!9`vvs;##c8ybj?ZMB52)Y4z*KxjB)gaA zDb(AQGNRmAc*t63-%#B@aFV8zTk?Poj(X+QU$g}I0S%ys?B_4f5j-d|zwTRW$w}y2 zI9rr_J%T5^B^q9C>vt+$kJXUX4K&-)J~ir1G466R{{1Sk>~d>f^rqN;Iht?XAbz|Y zZgPdZv&k5pxFw?LH+y-fS>n<+xIMfA7BdDcv^_S>`uFJd30>P5_?^6#k40kQ5y z<+}R3)=Uv_rMG*lpUK*&OJ|6gS6<+OGsNvHBYEj`(eJwz%>%Ae_B-mJ>6Gp15J{Km z;+OAm;kZ)NDR1ab@sJzus=KC%$npS$-KdUw4`WVI;L2j#5mF3@c(%oQ!fd@@p>s*D zZozXywvY_c*_e{rMn?*ck|HcLL_Kpu3X!n6sgZx7Vrlj0jeu6PUemC*RMH*CbZk3o zqVAdkvSJlOBy75iIi_skALiY|ogbq)?KH`yYD@eY<;AK^#j0ITD;U!#74ec-oZ%0n z4%BFEi&E98&d7Uzd;MtVsz)^owQ9bsdDRqSD2?@|OeodgBv18Z?pd^bS1(L+nn7FJ z@u*};1@6^HQt1f6HRa8I8?t&#nc5a>rRT=Em@ja88epXIwL;4q9CY>u|^pi*U2Dcrb~i;3MDipJmj%F52n7@?I% zdnMNMCm|OhWb;n4{zOCOiJc!|kBT=zgvk+KT=ma_Vlwg3=28%8hLoD^X19U%CfWeS zvdqB59Mz3Z*<gyVS<2bO5J817rbaq2jcXmbfk#1cd)-9r*F#%@WF^+pe8dQ$iY zOXdD#XT!9!Qew^_fAmK@f@N#&Wx+`_@5VkyACoz4IQD=t(+ly`6R2`tKadlmE+<6c zM0Mj;bE$JIsHWypSRD<)+D$g6X%l)CF))wtARhfC^;dMo3J@_(aJ~4WCXa=cLKUJG z!|}uoDPG`-q-H1PJYbm-7JJ|O&IZS79j+eCcD6c}Ha@f7^kZ7nYbDx7qT!}Q>mW3H z!V|4_K_PHO{4Fs{agfeUc%rG6v?+=uraap8!VV+tIM+;hKb!DM3DpS zt>2?&lN%YEhGvr?&DC}vOz$wQi_G>U7b0RiQA!g?8Wn$ld0q|$N%n*XQD<>0roM*W zvf-r-?EYrljoyb?UU9$W{$@CI*+}3}lLvqzD%9~>&;DU)l^$X?a4zUO&^sc*j!#>Ww@Ee#b~_MWu?s!QTM1`uAZNW%3udgq#ZN~J7^t?^)J-b5`JCXO+1j`X5f~fLft}XTYMI+TDI&)L`vcievyy3|-A}uiGOLrFN9qQibTY1FGs!GMKg?rhv;2Qo zIAkIQ&AI361&~|i&Zug`V4$^y_WU8*VQPlKDsH9H0*j~FU|KVcTx(h)*)TYl zasORHm$?7XaSdKT9Mk6P&5CZ2CEv%AmH@M22}#K^W1OSmIAhD9wM_nHK`};{jw+fP znPfdgYaz+r{dshGVj4@b*FOk1-`fMXfXjQ+tU@XxJ`Z_QaPQf=*u0)?b-tcqGMqxo|y|9bIwyPqgjn4me15 ziOey5PQ@@8LRE+bKkgzqUUsjo`?0Lw9k1bU626{ryfUX- zi)sULV=IK0lE5?70#3_uV;dN0#2D|?GLrBiPJFzK-DTMXeZK&a%$4v z<8;0t95-^e+7{O2(d(0HSFu1XuDhoT+V|L@QQgOhgo^fAq0ppt2YO6BI+hxMoy?q; ztG!M`){A=AdHh2Q|Pijcb*sJ`+;kY@&+wTt4@GKqV*y? z6?+qZ$|~+wc(boWt?M&byqI^rsp-fVDr2?gAgVrN#hcd~3m@+w1FT@tzU$v}K6j0H z`&KK{3v1BY9gd9&R`Ww)m=yeZYW(vR`vXT`uf-|{n`6t(#^(14d+2mX=Q zxIQMXt0m^g`TOH+!$Yioan^(o?Kb6pp}76Wc&|zu9iDKUy`zRM5u@%d;W2sQ!rjH} zH8Jqsb4|}JrSdy1MU+2ocP;8>yRYYm0FciTKizBN{_#Lcm8kECCilDhjYaCw)o&mj zJixjq=G^b>9;gd?ix2Mis68ncCx5(;a7j%$B<|k-oUhL*+Vy8sj+eR5|2>+|c}_$= zNWp8j)(1AWTloHyz~&c?`DX^>VR_=n!@pQ$QAyPeBa4t97+EV8;kaaEFEYM0Q%w z*SfLEd{mm_4>#74u@;U94>q5(V8@%aSzpFyjCcI##Tqfz{8ba{!{3oI!NmG|eK7w~ zUtUd>8%%62kDD#MyjeJ(I#Wh_v*)5l%|a>gP#=wdJkJmBzIj2kog**1n za9s0dZ@Sezk>+erzTP0uiBdFVaXfmG{Hh`AShx4T_%Ec&CXHCnx;rN#A2tgbViIPFeM*=3^&HSMR-# z;V!X$tq%O*VO5fCA?DGrkaRL;`)@v^_n0o3Z9R4t!d>cAhHq>$orXP+v9mR>_r>zW zf~j#@p9v0*>PJi~-@uf#wb_Uq(l3;SutM1_lwIY6vSmOEHr0(Ls?F=zhtl4H1!W;+ z%H0GLj=Sa8k+=~^=V-gv3KGLvdJP;o-0PAYLIoKF(!9|JJ~#pO=|Jo!z%oy z;{u`t@TZgUq-fECPi$@cv;%*Kc$=@=2XY%;`wvN724CZrY0VNGXi&mDxPA5{hS(Y% zcD2U5wM zF?5HthO^OZySx_8a(ogt(p@7~|AucLR@vJ0+%6YKuxadDSsuZB**1AMf;IC}AJR^G z483U#eNG0oVkvGQQL?t<`Buz>vCDFGYu1kaB=@yu19{hHW$iXBp8et&+=k(ETMli@ z?CiGV`?kz%G|t1cYbo1yWW$UZKucgsvOw4MG z{n_t;B^#AlW}hIM0aFynOSnJ9nn&ib_{)E5zbWUE8CJMcYj=s%3-Yar0h zSAe9t8$|T-opd?72b*r32GVo#>mKY`zV&s-@Sbcf=hFwtpJFkx0gkXZcFtf-h0O3* z6If>rYb3iRFtg7D5Kn3Py^auGB1;lj{f3#evD&@fBNCMd{58T~CU%^h*_S0}9Y25= zXns}bP91L4;Y&JvL5K5oI7x>?bl6*m?RD5(hjn$xb@*$Evz)6sJnMx1p$9bN(Yt?q4-cORQ{a7&XkS1I8V?FuMVmYfH>yqfPhCJj< z*Mpj`Txz^)$E6?|Tmhq2m7nvMA=5MOY-X*4POqZE24heo?5q6>iAwD-N#5wk!rZsc zpi(<+oFSW-@hdH(Ofj?A5h4CW^*yC|F(vAc4UvS~AEvz!tGW```gI}&ZAw({M(E}u zNaRn=cmC=F=jcv;sFvA!0&$Er_~i<68{O>kwwd+x&7g^Z00G5Lexi}G89!Yu;}V&F ziyTZI+sN8+yCby`tc{4Gaph`=z%B+N@Z;}a<=g9XnU%mh;whQGf$ zkfu7WrR`pm$H|EPtj*+CAcg3rkJ&!O4F{Gu?KdWI6UH4ZH7){+`d)npy{w%OvT;f{ zY@*udlvk?PX<^uU{D6swmtJlX2tQP6?1OCdG21DqV>NcO$DC#D=qGGAoG%@Xe!TK@w}!XYXjV;sc7y*}eJ_{toc`>i2s< z9K1?Vs(d%?(ON>MD)@?-9q{YC-5ROKx;@fsAZry!oE+yHW9F*RsbD%MB7Xc9Z8M%; zDN_cr#H^L%xVc%Zv2InqKkp`;zSBXYs8bjUTS_?0H~MBU&??`{kh$@q7U8*(hu>jf z<=Bip!NNd$ARLVnu^tqJjSzwH>Y?RXtLgRfCS5pCw}l50iR01@uc`Sma1d)*YXxae z-7N9>qu8>y-v8W~>nM;iZ0Q!^~nT}QYra5G*WEN8k{zjE= zH;0^+%s%HsM@jrxY!yG2B0n9>mh;3E*=Yz&w!6$6!m_w^r2GR(o;Xr=8p_u2JtO4d zp=>|@Yo{Z37~9AAyF+B}5o{VCFhp)2!5Z?wA@axw)`@>VSpGAD4dUwu%Yh@AAI}~v zr;cO?_$SFSB84q(^FlI>TDyz$>)_kC6e&|q3}P>I2b@7`G1Y>e#?LT@`OZz2U#754 zO=7Hq3_@ZLdIGbPiV1VoG#aBVve`J+gw>VZ#<6j}IBs1!aPioc ztCtR4{Nc*6DINm&s{3&Mc65X)v1I`^s%;ER$bJa0E_fbs66l zFFQ_QE%^L+Id%$;KbClTUGprh78_9}Dzvn>)k<)%ff)O9l9Kgx_}Sdpu}4e-{E@WpHubN#Uu5^p&i z&$4XB-gHDi$JTMSP99&zy0d&ot>w(_%6E5@ujaDxJf)kwo6Gw1+TCRDJT^3|!2#M` zyb42T)%#uxfTnkg;k!CbZH)~{Gse*-pG z5caQhiX6O_1q7v_teQ>r{#QCsu3O6nWv%)bIsaej_5Vs2{3~tNOV8isEG}Dz)$OU~ z7wGiZC;1*qtgfKgF(aP|6M)^T{l+NU_* zS()T@vj3~BZ*3Y0>Wi}K5uxFYWDd5T4KQSFSg0u9oMTD>VDD!X~fwd z*?=!+DM}OI1%L%Fpk{;v7GU#ECLEAL3Dbd-_t6_ZQ6}(YKoeZQT%rSR!czde5O5;A z*iX*0EF}X-GMXJzyTlP83tUzI~GguqyWzb%zsr;(t#7K zyB>uB4+KnkO;I)gPX!$K1_}e749MO>S-`ykQ{O?0u946?it@%rMfn|g0bpwwP?xWm z5(#*3v!c`ouDq`(A8%2V?!b!yeYc}v;DLZc3Q-Z@sepNS3W4VXb}mA}z+(Y#?oyQV zz$*cx;U?ZBdKQo=cwGL&lng+>eTw1+JP`2nV)*F5iwPc7l*z!80T<#~3OpOI^ASbK z1|AC-cT`ahkqqFv<7g%D0>HQv7-FIW_B@G#Ngizty5T|~D4Zu$?* z3-Cg~u3uwZfyV;=@q?nQ2d-ShQ2eMU_kpJZ+J8pgOK36R&9H^99CFD2uJ z`)@$UPg0u?VM@vs3=`l+JgLA70l&nP4!jKTM?4w8D* za=|DE>^&SU0Zyi=JcTJ6hz{6)6jMGSI^Y#Nhlmb11AkWIB=8Kt`D0LN;MssVc)kIi z54aCcIq+gYZ~S46YrqL^&cc_9L?K}0IHueOPH-BY2f#A`Z{ShBLkj>~j7M((j|AL< z#|?Nf;GzjksSUgg@GCstz{>#xQ<+j9IKi!W0)ZC-x=+Ip0MGIU^2>ClbVQ;Ouyh7f zx&to*{2Na!aAhV_LT6!^fJXwJ#*<7kfWP1w2D}o`3x5S*EO2kY|KLdlUJlrAHiiRu z>>Q>n$FmrCHsIqd&cJj2 zTl5g%F+AS@F9STc2zr2*1Dc;f_XAG_?D#Ai20RvU3!chv&~wjbG3Bl2(A7v30J<$j zMSuqaK3s;0h78JbOiDc0K+gd5%tlv(?hW|<^QaK;Lcq)ypb&UAV4EED1n@|}&vKd4 z8hA0_@)b;3N-}^`Rxu?c3$4ip($I=70wWNxh+IKE!~POckmXh8e&rd=CZ zsF4kJw0Ix!Ha1v}*uu^ML~Lcx0Bqcf1NC58cN?1puzVZqmN=LxRaLC>PQC|SwcGVI zwyBO5 z8=KWqsU7K=>XG4|>{jer=2Ffp6=w;9XJLmbvqbT4x z9>$&Zlx$Dqin_{AURPZd<;%LnyvGQX&3THguJM(5$_w?3AI?$!Q{VV-wzALH zcz3q)y|3|Xx^mdhxH?@~*TDENL%H9;XrHUtn!wMQudE8V?uOhY;VAuFcskJ4t&J<@ zDLHMjITCNTqr|ax#`iLnt5I3TFXt+y?TyFhDqA`j56n?s>S$cGK&k3vTt8p=rnB*n z`O1+l#=Q%a^U=mHGL)?`#=RNJ^{$0=KFv_ROEgx@RW=SV-ddn!TZ~`KQwlA{oePxx zNtX=TIAM^Ngrp z{yY%&YN)!v_=WWRfX%@KSn>f&1Gw^mHY)~x$R^-QSaIo%UDn;8xKHPb%Ott&Lv~KC z`-qLjth)UXd&peqPbNvXoh&emHfP$lh;z%5fwTI<-eigII@zfM6Odd z3drsTItnAVrvMf3jsV02CZ`}4{Me~zSRd3p6%6tw$kU)r*&Briwgg@T2q%D?0`jlC zspgO+4-7AOc#fP%slGw6^lOM&c4BeO>4jI0=0FS1%>t;kA| z`$}$WHlBPu1$c@HhZQ0x!V&TbQ+;~}HiFlSzmv7to1-DlPYP)HpdGD2D+a9`v}~Wp za^}WdSzk~816#(d^>B!5es<{TLpuCMhc|WT=IfLX)M1zoJL}M_!)l!=dU~o3m*{Yn z6J{wd>4XhB+@ixm9q!fPAswF7;a55=*P*JzN*(^CLl-}1h3n|hUx&>#hzZqNC$!gL zbM&P_1v=cP!=pO9sl)gN&Vr}waFq_X=&+I?%2y)&or2jq+@QmKIy|aF zmxfL`ZyiSKaIy|FbhuQ9t8|#J!}U7+M29BJQQo;= zXi+cnH&5T3u{nG5s?Di8rtc`+v3rNZ=K?EWuYS#9-OKej#xLuF8pzbIS)@B1 z$EmY65vA~Jmh6!|)JaJ-t_$*$ftQ(&TZ)tNl4jrJz{{+I%NUU4_v4&@Pq5%Ji*h3h z6#XPwyunxAzRc>mOdz3;bDctg-#`;lAQVfbC4=Ce>LgUSQeoL7GR;ZRoxwsPO?Q(1 zL%GU{FvCe`MFj6}L73?zggQOQSR(0Rrng8YqXr>>2n96=Bb}Mbh&0Ek#{!vGuRKmi zI~&#Z3ajTf&nYsU)ZnD&<2K((I7I}U*nHd;I0@+v1i(hJUF9skAuk{7==&$@VQg5a zcX$B0fj(@*IXE@-z8oSqJjBINq554n-}iSg|WF z8{()i^42a5x9H&@K4HS*^a)Fph*mez$L@w=+17(M@{ir=td1AhlP646dW6FxbvGQ9 zsUF;qmC0v3cr*)?g&w>?Lkf=3_&8NfR^HcD?Uv;p{1}Rq+dcUJ{^JmN-;)RMjG?k& zZT<|e93o$;&HH)=>QhGv94ddP%~SZSp^hG2d^Gc~4q=hv8K|R2B(`36LxKFL4iELO z4uBCi4cywZHEzD+W*y#u@$8`vqcYOg*)j$()H>@hs8_iNkS_y59?6XGb z#91R`P#}->tR}dCkRFJ_7U>dHSPP%Wg{>T+73PT`veJ}Pm5-3$2cnjMK|#p)OO1?| zYh+YLQpQ32(~%NC1#D07ZV+$F_l}e;nxnT*j+FhH^R}Lg^ak)nDKfh`s#Z;K0pUb* zl<<6we!xlRQ9r!o?1#QlsG9p1HTvPOQ^?(rqW6X3Szd#UesnhkR)1hZ|!i$7s31`#>ct! z>huU`@vJ@EDfDQB?(0(+wT>Moo3-G6<7Vp_)+2-aIz2HLiIbYv!G?8(SsVlE7OA+m zuc7aT#zu3Bt95$zuEEw?S|Hc8;LS1Z~6@e*6MzIs0?b!efh!wwB<2fbmO_* z6(O)TdLj)81ND6D_s;%Ml%Q2^I5v)UCfNt0lJ|f_y(m7c`zXE4T8a{Y za!@Z8WnhDuS*i1}&MSxFo=#qkI(uChB{xO!5WJyvG>Qj9-_0l<)LPHe1$kT>>WSt^ zv}~m(Py)|HgX;W|B4gU~W^9xrr9BT~Y>Zsefd{fNjyF2+2}nBr?#N&9WaH(!SZ)Uh zh~p6e{o{Cy$9TQge~p*+I6fETf8uZ{H(rkKtqECs^BANK_U0qmIN7)l4+V(p!&`zl zqYod5)b2if9zaAqZ|OcxABT=|U2k4r&WYz;z}g%S13FHAS5u;)QzAHlcV}Z|Y66b~ z*qDGq$I4R)D0Hl>)t8S2nB12y0JzwfM|h0Y^_$1a+WmNRLp^=OQ&AASQ2fv(3(-aE z#>#R1_$l;6D>Ls2Fw+eEW8@Aq|C~*5%u3|L87^B&`t$urW(~k!Nkmdu_+^&rn483Z zWh_-r8pQW`PSX2FSvN_BC!=<$a#%7?U=!u*PEelATX;^amdcnYeFpQ^kFym|l4cNC znp`-Tw_uayhQYj*=j3Yn?8)-%U>*i_e-FmcO_WWC@CGe)b!KcdyrEtLI$o{9R>x>w zA>f1kPL3G@&6DJ+Av~ysF4hK<`+z1kn6JZ+ZW^XJi_&!UPLp2_;r@6tO{JV^(sw9t z4iG&QWlxr)onR>uC(E~oVrQ8w9ZupU0#oGQ#G4|U4Z~JA#i^M}9EeGeI#G%Ivng`m zFimP6K(887J044!98+u5>nA-yZFGUUnr%c)gM4bgo6A6_&;2T;x%)KfG8`)}RrVgv z$9U=<62iu*jxED6)z~!0_apc=K)E`FkMqzSCc~O(^0yS;yt$rMP@Gqpo;ZpP^#%x9 z*Cxv@qj+dN-K!!wrA`ibvEz6iV>9K|@jMxz!vqvPQ!bvs zTh!6rxv!zE;g*?l#{?9eFPo;K=<5?ubPax^M!KW75yCLvEQ-_g$#+27hp}0Xb`znV z&63%Z_#}WIChSR*A>oNK5ba{LxZvve^&*T==vEwZMBAes*Viy01 zu{mT!oINivUnjM3{ZCw@8L1uSq}SRz8tlPcdvG!Uz6k0 zcw2d55mxg8$Nfe8>o9VE&%ezx82O@o-{Gyu_gwG}x4{?PwGsX-Ig9?AV6x#XI&AN9 z1a?y#?``3)a~R0j?R*C|j9c4z4%V)Bgjt}=@_Fzu=C|G%J4)EarVM|Irs!Nmo}xmF)th<`z zLDi=;_+8q^f^Xaf;Wo*c=Q(Gd4@!9l_YB=HG04HCJirH2O{s)v^|5Ze{P+~9KgvDH zJL`HBo+dHUB_|OWEnhr_i9JRhKgQcd=?+C} zXPuP#fETeCigh|de%dEqr+bf;5yyEyH+@&jmLOiJdB6I)Qei%DpGBCZ@^WCou(P$c-nllh2f2 zpG5Vl&y8C5c~JL7!m^jml0#1MUJZ3_Q@u1AybQ3t)8)2P(2*{$o`NqiM|M5U0{}*! zM%i=Zi>Fcce8<_-{A5#h)baLBelQ68REeJ@@C?-tNrP+=UU&98*k{u|b7tbn~ z0yB?Huu=GIfCrpBttMqRDA}`WQVxPrfavz)TwfDqeoe}6M1h0zSf;KCUHeQ;N-QYJ zS&M6u#*h%E%j0aTKq-8J@;N9OOKZyf1xhRqE01%v#fh&nSxnzO(R1y|7{vxlY1HzWe#zJ>*3<^vFau`Q^wA4W=quYn%oRfif724O`vqU zhS%hlgOV{vMm9DChd!aSEBw!7+#)XBPW2A1hs)f@j6tvv+4E&4F>q=7m?4O< zSZ)KuytpQ#CCq5SOEpV<7L>@BYqoABnOa$s@;)f3@=|j{KvpTx;`KF|7eEPoqo&H6 zBwkQc<^d@APf)yY7>(UjQzi(M%FQ+P#E{I^nlkvjLRsCnz6Lu1s0CnISk4moSg}@;5MIch$`7fg^A7?wXW9Pzv_eq(p;KHd!u* zx+arrs>}38SY#?u$nTjh*VQrjMtC%=CRQ7k*Fj1C%&E~ilcqO9J+kl{AwzH$F3lfT zuu&5*a8dP`5(-M;=QSyPLGeCVlQNTJN^2Id43t=0PCeFh6qG>RKsoDLAAm4Sj=}(g zPd$xMUEgwS3x)IK`2a)kNM z@t~w$s!5pzinrF2)NiA-ev=u224lS^y7v;qGA7GafmjkxF!~|tn7vGX4o2N&)fJ?| zzXzjG4h}Q~bZ-_2G_59cI4H3{)T~E3D8)@{%B%q;dyyO%gd!K!WPB2Y%u8f0F_zS1 zoNtbP4Ul_-44pc+tjSmqj0)h!^zjgU0ZRE3l=GmZYWtmxZEkSKKN6`Segn+nY$wH_ z(8m?qc(I(_9D853TnE7OU1i2t?4Bo^8)~yyTGqj=K>kA@Rt5zd+8TQ>CEP`(b~DtM zsa*^nGP|xJM7|nqNH^MQ<76dE5d4-+LJXlET2VA05n4oCriK_?*lu|y#L!?}TQ`{- zVrbvQ4JI8{@L!*ZP_p^bBizu;H5&M#+OjQhV?W@La&S2OfDu41lCA*P=Yg7CW$|*< zyD;3~X50b7H7D;sK!?@j@fPG(P(FARJ}uxlj0fd`yb}sduLDi5Dg6})wdM0IAYDto zL(Hc&XQ-fL5J4hxtiOe30`teD?>B+a=0Ok70dnMhJO6~ zV}&v#!Vt#4sX<5w!Sfp8ltxSi1B-Sp;ZI1cNu5WnkPnCuQbz_(!&rB3Wr)PAoZSk_ zQ)|oo-k=q%#Jr6p=dOwDHO^40g%U>-9y*n_GbT^WoEXw2-V-nJ)ueP<u4JA?OA8M=E)@79Lq_>G*QZi|yoWxyW} z_*Vs^J>^Ww@Dnn~czo3ld`hOHr#wV5v2CDq{ZmTugSCGtjr5ew+91R0$ROV$rIAm` zwDy#XN#=J_y6P#V(NARv^OR>P!@#!4u=J@6u}@_P@szFGBEuHSknvOo^HUiDJ>?3@ zz}q3iKXOST^@~I3NPi1%@KZ-I0%hv57Us>zlZC?D8{C(-(9`Y+7Uw_L} z^4W3nVJj4$3AQXJnN?58tcsHr+(+y$111JNS^}BH$z%GW7o%NysM$HJlfgNs9t1mr}N-EhE3`SZ_4Rwcz3}+ z&MWLK8+9=>=jCy-XBWdvxvv|%gGH?jbsWWDQ*0epxv@;`X>gM_qjBoal0Gpgp-g_< z(-7=Pi7^B+L^HK>xkJ9$72Zf$Z|y`b%exv}yepq5eQQVgxu2m~l74%u9EI~iPK`m) zC`;gBKcaxF6jN>cRM)?)Zr;l@owVa@f{am4tPAI*Ss!>b>${vvUrBC5wT5VhtCw13aGLOJf2@qtkj4sxZxGjV4uID?Da z(aZ1&t916?&oU3;ji>ZpL`{7L7#hD?2Ls|KZ^b!fQWLOMNA@#o-=d#42X{Vr3 z)+E%u_>p#M%{)-NE!E`_Wd$gOPf#|2f`5~z*NJ3410~}L%4tw4pP*c!Tmx(7`nv&0 z#g9nUjmsY(Ur#d3G?qUa97o_lgN4CsLNkJ8{vboB%OEI{_oic9_res&39AhMUs+cJ z994CNU-pH#1qw?xK#JBZ5Ny;M0_9JBR}n~o%2H6*1hVUcwW-wQ#WUSsqZ06&-m0y7$l}k&WTZL;gU514xc)b z6r)FKMwVE+QT>KUz>y;P1PCti$R3meN)&NvEu5W@_FVWrPF@r0IztJn$r`bG9lDT% zeW?**$HRJliF4}^>wYPAeLeW!6nUiDN9i?5ls$o9`wSNu{_yxj`x9!t!!uUAy8z?O z`6MJ#Jr6yJymBZbMPdcOgHHimDt#4mN8&YlE2SotCgfy zag#`|66pmSDI`hTMaO24Jexs!*+y!Wq?q`MNJoft+(tsxr9S?PmJD4DINXDSyW7U= zlDvJy%WA~oCE}g3@#2yy5}@e-2PNxT^Ga&T~; zwegNiUOw?|5$^>$B%0N1nmMS8f21|tb^yYKgL}}XLh?|v{@o0JJx`=z8>v?+n?TB1 zLpY4a!F}DvOG@55v}nCTydM!S(^`({xzi_kmBfosqB(JJkFoLkCC^E`FNwFCc)2#- zfaJY71C$UYWF`)7`l<;}vy|kmC*Bvt+ey4{+jxVLcfAOX>=OS8V_R#cMbg)lzY>!p zsQccn>ZF_$R*}4^;f zn-pJxBdw^XaQ2JfHpFK>#u(-lEwn=v>*{SdT*e`1jx8EWndq-oXNe{vULs<+1;hdy z(Ww)MDRs-i5yM3E*oZEj7$)KcA~rn(qSHpqEfdKGOspGlS_~0q$TkOibY(3x0M8Mz zz7>aC_{~wEpSyfb&Ykjk5phG&6{R604dV$Hck=8cBzdc^|LB&lxj<*_*SHH$#d#~ zj4DdVDPjywD`hnP4kDk$WBBNF?N9hRqUt1c99NnX~7d3c~SZ7|w&vtM-mF$DFC{WNZl$ok*H>)r_32-s zRpNTn5XP$z2FndN<6gzmv`$G*B@nUlhiZ0I?lu|BpHpY43PF!#b_4P2c1q+75U+mT zp#(~*8+7?43z zHJy7Eh)a)&?AeDvXs!eeQiw?dzIm-2H~UIb{{+XQjekS)mq3y#9v`v#;lQJ1M(#<~ z6M&TCdgJzOL#z;pg1f;MvVsu2ELU#?;$yw9)i|CsshvPzvw}ZdZ>T;3h+Ewa#3+S& z&#LjnrD)rW8>44%GGum}EVu+DY5JIuVIcGjF3#2 zHO5&K1$3A-b^u5hL*4>HFBwk+y(-j(bYiYK1x|u>BSlgy=AQ$RyL%qJ8#5Ln1(Zzj z(mq@!Ketd(`dGsUfG8~7uMmeReIG~>OQbJ=kf7y>%tiQ!7Zx3Q!`VPQOv)Dp)X$pV z0>sape2^5ev^fbR!Cd+tkU`ei-+*|T(f=SRSR->CWaZ|fD5B0BIAJlaHJGC5K(4Vo zSPW!+FlX%tQpB1)3dG0Ee**{))#)gpe4yyL zo(AWE#7xyy75I}WRTvMfk3-;iA!5bq3J~sMr&m$RS%k*}F*Zu5Jw=)9jlxDoDg?=@ zR5hV+T~nk+NwDTugX3pXB|wtqRI;ui7z{zZlX)An7A(f3P`5)ifg z(F-JokJixn6A(TI&I9RY0lEw%M@g3&^9XLc>hr7YeLj$A5(P|-L6moTrbS7xOf3RQ z-d`|G>;&S0YRkk!KzQXl2E=Oy5oF~QkQB?00U%vbmLUHnMXax`M^P|Y&n5vGfIkIC)Fa4(o;&NEbch=|D5|gKPy9GfCl~O#MD)^>%QMtLa9^fpCLfB|WV9-vXiU zG+JH#6i5h2)Pg<(6f;SI+O|+Ls!J_v{xUdZtyL}VdJJ=%DMj5_3*pX2|P6T&7F4MRpWo)>xs2 zPJuMY#QqF~D!1jq&w%iHdjm+Eb;Gq>R%lD=ULbzfjr-{@{`wIeD_drRlw{4X01|@% zvae)H6_B`g73R@U{D~wW^H?J61ag=mzW~z5G{3tXBQRl-w8&Ip@S4u2#`+Tde@ORwlseC01KV&$b0RWO>QV#%$dH_LECp`*C-T*dwRstlKd8i&p zh`BKW#KV&C$3R+5sX(NzF*s+I64RuNwvG~q{`-mL{5Jr3VX|e=CqQaVMbwY~0io%Q zdL?@{u9)!j5q7@~nGYnFwN(rxl*uY^6`;6&BSNYU0qKWo%ZymxXKv;xVNH!Em?<_XEmV$aa9^ zgkRGt@KHcs=IWC`LM*O>K*$+ZfUW}RH945<8iFVo{j4$9O03S9R34BdA@D!?gR9mG zTtFyArXqY$ib%{X7lbSa!Y8|0AjQn+79eRga$A#j+1kYn+C?1AAha)I_6i^uv+*Pl zUS!?@qL?PO=*zP~ruwhocv*g21yXLBSQ8EgsVm2&)B}3+>MAt%P*?qrirV+$G+pLKd7UGCy5=N=m~X1&-C}yK4;I=J~e-2R@!s$BdzD) jdzSr04P}lF)HQ|z>5s~fvM7oex~NX=`S`4Q=Ij3le;v=B delta 35891 zcmc%ydtB5-^goWjW|l=jU;zd3&Sg<9ir^K@%c7u*E(%`pj-iBiR5U9S5>eN66`gX- zEKR$pSy@=CW(DFkO-Zw~i=o|Kgr!*-mdW>dUb_g^`}6y~|M-3$-+erG&zUo4&di)S zbLPx>y-<1Du=A9mWVPpSU(e1}e09pP4I2jYzb0<@qBrvpf$Sj`A~IMD_ay;b@#7mp zTgklTIuD*zC^oWT?0~q+hO^Ik8C#{1)FGr6m0IrWrpG}6-7DWjy8%KKM)*t zpP$Y$MZGwVDaB=rW|(b(KWPnpE&!-+q-o?4AeTDHTZ!Cp0)9rA7$k8>Dz3|a(rnYm zi-}Ff2DSu)O^2Z^3We~8NO;hAak5Ffu<=YO$uYD~R+O=0Mw*MR$Dz%Roe%;Mn1fCK zvheT-Z((6dsSl=n9*N8%mU^^p<>ZJzPaI>xk!-UVJJqaH zTU$h_M-WBtT8&0k|1OKpJ{ z-I%#_YO9Nk4HCWDe9Yp6zU{D}YAaK$m&(=l?9wN-e_py}h(4*NB#0Sp#~3D|>Xsq% z{4K-eBdV#NINP=z8vWchhPjEB?WVNu%M@*zIfffw)>P#fA}ekg{Ap=snS^NfWOL~K zmch@b?P{i=vJ$vHR2bX$Wc|g2_MKS=VQp{5@A>wf@LShD9b-)oY3*6Q(%J5^$PZc1 zB7|FLNBnjQjbN?Ch|ouw5Fdr6ur%Qpmd?UUr-r4n=3!XjqT?Q^iGMEo$qX$;s;k4q z2jR|lzoObg@_!5z?c>r}pva5sf!|l+((wCxTu=PQ z#1CN2#oYMBHq9}v?2paX3UlJw1%u4??ph~gN3)uXium4OyB8lBKy0>TW270TS9CO? z+PAq#=-My187$M*+-R}8{h?VbYONJ&NSs*LwFPS?Ug`Rf1uP{wIkATJD~Ll$8k00S zY1C+^y@!}>Jz|zHCB=HQ_`G4=v#2*`pHIzJL2m}S*?P`m@AQee{Cl0n-eZ0b(YafA zi;0k!bb$7&*&euU0dbY9e+Gz0y2X1LQAfq*&BTswQQ^M^Xhn=TWWwYrlaeMTO-h4epV^z7~Z& z(%IV5PkKDbnQ!UvgfJHB%@pToijL+s^pZ5Cg3K#H-R~_nCiZ3CqB3zPe%ti!?(^Mj zY@}pkADDStKQXCyH+DvB=-r>>_I6{;OMLdwdgc&0eRA14ajwtTgY#xErTQ=IYZ4iU zO>|o=_5jIv&qY!fqX=<5CWEu+A7&`_q$+!-pD`n)T18W}nTPnT?~(pZNF@~(t7sRN4@JVhV+kT?qX^GVfa1We-+b5Z2*5d{TFZ!O9TSzw#?B`hwIS<^n^0Yfg0_UT zf#0PX+REDcS?u~J!$oLHSZMWhqU+2yvub9SZ5bg&A2HP%%+>R0d)!jQw3IIK&k|Qu zrMbLTSJNF+^v3+}n2RcyTBo4K9ylb6?2WuSS=JUs98KvJT~2k{Xvngv_n??lLZjVB z6%D5JSeS?yWaLd1F?3L0&zGiQFWWo&gozDlE?JAXC$IREhqU zp`19ln=TV@-r%{t4ay!EiW84kA^)EEX>fdS2ukn2w#WHV#XYB> zqI;s7IofNAt2q3gc*GnR?So>{0i9FMW0SN=rCN`xum6KqHdw{418k!r-PEuD5#{Ef zKE?ORw+0yeERh~+1t|9gbZihX;XWq`@KePG2_qvuDIfO9g8w)}~E&uT_*w`ZvU%Y0F|e;*8EQ zjM9ITi8D1C)@!j&Gp6GNsjaU2yRI&!xU3*RBa!KD6g!87vp>a&VPAURCY}76YCWnJ z-WG2Re}vr@LFo%a;KX3oFuIDnq9!w}({RGK+=aj-!^G~f zIZ1}dlbAt|UZB&CsxlWH^S}g-tBdZ6#BtrDh!vVPMuK>g*pI568feLeNIFeVsKIx| zrg3vwT&aG17shm=$AsSa&705-zb{RQ!0$&BV(@!=LS!JP0|l;T9vk*{1;QUc|TUmdz?@E7}8phdQrT){_a8@g}=Bz3( z+g9MTRg`4gVIwS9NwH@!GrVLR4uNFrv9#>pDVR&QSZ^jz^YRrLEVBvvh`ZB^(viMU!5wbt4l4q$C7c# znP3-&m1&S`&PKl{57OSZ*fV*WJ%d@SCoR_Vsn%2GsE^D=<;-k&^GN*p)IWc?v^R-r7v9D~xL~s?_WgX3RM=)$Vr0V(;UHp_x%t zh$@Ud{uCA^au>_@J4nan3EJk0Ik7yyxz*cdu_rN1fDf{|6@iwz=(=sKCX+^K!f?xy ztX@ghlgWuOh86Kp0jWRNg5L6Iq)W_2(~Mb4@y~{JBOn~k56#l8YE4J+qHoaXl4Kh% zQkh}J9BNtgliL8hcVtrH$^0Mf=`0Cnig|M~I!?T})M(B=q75(V6P=~#Zhf3!Z_cdE zi|F)8Joe%H*Ya&Mh|`FP@7wKv2I3kL8`5p=JwO9H!aqQ#yIs9Tmxy{Y{5^0 z#d1_l;4Nb{5n$h}wba+}XnXOC7{C=}o zzcE{n*92%GY;t03{v={9Aa1P{%_rGPEm!Q&B(>1Sm3pGAmJXV3^}|YuKZQYQgVPQ; zZA1|m(TcMC-R7cXV-|z-MLTH)vNZqnx7m7ySaWH^Lhq-kza?fH#+6opkjeTdR@7qs zB&p~cqa`zcQf=}4MsxA0e0TNlFBGN35a74pSVArR4jZ=v zVrn8}R`Y^$sC8g}3ftO)KpT@nndo+RxSm3deaYZv`oKkXPp@m$cIQ*W3t zO*ah_UF``?(hO-9M`&KwR-jm{r_!u&pf^ro=1E0W%;MR(oTiV7W}d2MV8I5i{sulS z)%v%26_g3pRnMe&NWF|%h^La9iFaX>IUjMW^yoZz;Y9I! zo`qkRV!-^ECWE_!zxSsie}3QeFE5aON5iKy3Y$-;$u~)fRZh2d5u0_0y%YQ>tw6Q< zD%y~%colNY*1uqnh>RmsioY|gI4RhI)`mP1ij5gp^f$KMLJ_kd#E=SRi*>d!9WC4D zcNgOq^z`ZkT5{sIhIL28mIdvFcgL_6nq8;bW?%H6Q$5*7J5$ee6XzF{`k$c1vduPT zl0n3&H$M@Z7Pj{~_X&0Ot9s%S@!`T8-m9zVnEwd#73=dq>-1+FA}ppgs}EMaD`5sE zkZ1-ss9r%qednCmuqfF_o+Dv*QmrS{wdcgOMg90!H6m>BxK^H8#bvs1JaCHqdbaeL z#qON%J|muAl;9)JG_*hCY_E#2$IxD@_#fJzE`8=PYJW;R|HMmS&uu1qyl2ujtzjMA zQ_S}HCW}28v$oz+pL;=MFI~o3mtI@?Hfw$9II3+ySi7%jZBn!08_H7=;P(Gkyt%9e zf91G1xhw%kDJyt{yB!xZ3koB*9HSBbs%``glAPKIozPBf@$gE%kotbTK^@k}4@GDU&!qg(&6A$B)-s%kQTbk+71 z8Pcz6&4<)&1$d~M_%0JeraJR0@pN%F{``mHRB>0%z7qcw&*T?N#Eg|~Sy1V!m3W%@&K+%lJV!cN5 z{d$!sUwwLl`X1&nn)VI6Xd7&vKU-rn`~)ibHWN5eq1OK8EJNp))f!hJooxRK6(CMM zMeumd4-m1RM7;W*_|cleP49`=HBlkQ-k=`tz^$d=$T=EWuiWb9=t(_yNIbe`FHb!r zB5W&p`9blzZ71)3P>fg`9&9)W!FIjUndB!ML#pSK(hnXG_O(_%=70#cPvU<%#H02( zG20#5SpKtktBG#!{8$X8`UAD@?WjtPvdC==ryPNNv;=9b*RTYjl!i|D0+IfFsc_H;CUi zMfeNU(;Qz63efu*7^Vtp*auxwpGe~eHA0jO_J)^cZJtvMRoSwC2 z>jiUm;M{&lGBLLxI!nH0lTp^@lzf$W@9%vjZphEq!`(BI7r%iF40($9=|a%lr<1?!maN6$qdW zE20rCk&g^VQX&3!O8j^6b+{Gy`DA;G=e9(Lyn2+5@-8pD@n2KE=4Si(Clz62;rAv(YdC_mz9RBNK@%pZw?GG-ltFs=f9)h0Z&r|iu zKQuxj2x^~q1mE3-H`^+@?M~y@pA*Y>FXWFtC;r_1B(L2f7Vr7CQ_L3UEX|4coU@Ee zW;JRcsMT+~6n+!_7TRg+A&Osnt<@-`haf6Nzt!%lBuP*{@1)k@i*3Vh4CbsQYQY68X#)+4Qo3T)nqj0#NVJ!kb4tBi3W^LA%`exum>b}RnwMseWn4t(8jrB!c-F+OvH zaNA$V)e^CC|F?my5xhgo7S$_8$4t>-V;N?_3VOw#f;&_DtQTG1P2&FR z#iDmRwfZ|(YovuZ235m9pdDjdC*FHE!<4m-q`->qsH1_xihf8fItnZLv7pAOoPHVy zLdTuyw-X2<9H^M2_PgUp-&s9LVW)0Jr7`O^XOvE|dDa3p1+%7;ZDF zjVS26Xsys4Oy>XCMC!p1JUN+p@XDQX?CaF@f3&s?n`nDzJ9pHO zwMCC7r`2F>V$J`-+HOq~-tVQgjY0hb=GKY`nz{L{A#>wy)5PO%wdWO&i1*%WGojdu zYI|U-p8%Zpb^#RF+fipl@eQiDtfrypROv`(#ntNjEnWoWzP)v#YS`OaDY_o+$2yhH zKWt*W>4&29{h6#wX}b@`GXD8E@#u$xSQqiZhi3lnVc}7>m@hjZ3d@EB@4)73xb|Sf zAb(M=dcP>Hl|}OAb5tWB&@(V}cB%bobeCDnzqTk6DeDD=9_GE7sBQ~BK5?ht;WIuqw2Eu%90;+$*z;NoSk(*^U()bcddl&8_MbrdTjxddNjNyI!C z9ZI!+X+a(#-Hb$1RvMB9)-XW8<74tb{@Z_WZCc`whT=}(ve=v^P(2@JO&sZ1o%{V>1xL*) zb@roleQy8YGsG;LUctjJq!sb#(2P{kigBXy z_u2|=SN^5u?x*c{S)EoSqAR!XyKdrmWvo|>2dIfB4a=*sSd|~)By2VaAwL>J0(Nos>0E|ezohQ;iZN(SZ#q3WrS*|$tX)?lVuZw;8 zKl4QD#jZSKt$6a{#DeT{tg{_{l^VH!fc0%nwny22oVs{k@+aP6u>OU^A*w5JKCO`{ z2z(YX#)`5<7j61ddPZgoTup$iJfNtIn>S+-+>Vf&H6iO{wtfLCK4Dln65C2$))ziA z`?L5mv&j5x1{ZTh#b=2vGUtM?_;$V@DFKuD)f{z89XMAs`@D14=2TZyfm`X7Mc zf6n(ICM_#*XAV_i!@FswZ%#z2(oj50Ecm<&kLgpo^YiYE&zUVQUy9&{+Krb33;NHd z2FD@mPsl@?@{4wytrE|F(TQJpPMrQCfuGtU^j}69I)L9fV=T_pC8F<_hxq6j z;?9?o`0ZSg`PFE?EmyqxRd$O#I1p)gLvpGAU)2GC@J1qTLv!9SSM>jSfM?SmvGD8F zt&fN`UoYpTM}+T{vHZsz@yL~xVRhBi{8x2e4mJC4u4HVEFnrUN-&`qTziH|D%=?h^ zlKRgJ>MnJZ#)9=3EUgpxk&3oiFLm zvj7$3qf@O{Y8Gog+TIC1g=FGs677nLX(*f^5@;tZc-msNo8fDh8WUK4u*IH&!PCQ= z1U3FCQVY2(Jj9~M<0Jjd)`!g2p&l0NP(S$X{HEptq%N?!1tb*KpPKkheysTbVte`+ zQnQ;IwYj#U8Dg%*b3w+b8IRb%X0N8BBwHIJC~#NnpG7#53cF8Ym(yVS*R^`Gpg40j z28*+a!oUVaCe|7XmS56S*Gr9=@MOPHJVR%g+1|!vF8@}iZGwqNWS+}L4TC3z7t#U| zC~}{sg003fz(&WW`Gm<26bLN~Hu#%Vb7BQ)SX*{Hgfd&VY7Wfimujsm4#2uMg(bxQ zT#XP=%V~-auLc*OJLlsYr=}6^k(P(qy|uBLpy|{i)uUe1j)KOTZeq5L^|M%enXSX2 z>ESL-I~A>`dC*vJpVE?G*$VP$>w*wQHdcG8OYfDW5Yqb?r{0lCfa0Gg!?0wZpl2v~ zspVC=Mzb$T%`Q*18B*=dk7)09h;82n6`a75+(}KmvFHu?OwnM2`YEgvC2)v|;w&1|h8rk-1oJlqvCYNb*PqSV{;;`YxCs2E)+5F?r-=~Q^ zzCHydhBwPoY;GG_B_7M+sc4Edxfr@itStBePnT1R%7c>K(|PoFv(3^g}SFjap&#dR(rl3YhGH4-Ntkh__ZY<|CpM(tLEh%3+)J7JK z2Zndl!?@~^5);(VRy+6BSwM)x19h_=6n@4Y3u=oMEB;=kpLnF?Is!!2hWDp|-Qw< z#J8M-UT`qoU;ktp1riRlWpmZuh_PHIQ%!FtK0eFwoj-ftBfvJ@IiV6-0*w0f!&0Gbdks+&x0&)#Ow z%De5^GW|+W&N>Q0SQ=xi<;S7S*z}1AxEVbqYsWe6hO#K8Uje3YM|3zl#Mlz~y^%e` z-j+{9ur>vFuAn}(kVd)IsA80=F%Q9d^PCdIx3kqBcVoHBtzW3azPneKW?KY18fBbK zPhH;s`5xW$E}LzA-^P2YwjP!E#r3X7@^Vbo*1|prdlMU`(VcHRYH;B@ZdmmkOw!hQ z2U^IW4$R2jk`HxY*LeS@Wc!XR$Ky?s;B)FtY46BF3s5qw*hztpRlfy=nh?5<*?m?L z8?I#QQixqoe7JDt5yp_|DE4rNs(W?#(}jx&3;fco>G)=9m(KmR&Ox;A zY~7M=fgW$jJy9&on7)jp8`eFJa0@-mN4RHCV>n(9z=3I#S+`8yjAE1epaz*oyqrJr z0#y9TVsFeRaDbDF+||T#o1SYxQfLEie`iR%1HJXiQd*yf0tjNE|OZ4mt zj$!MSdv&(X{r1WwG3*g`MOMYIK=zuv9mCrBRNr%^j_|ZFaH$OK#4U3KG6&yh2G<0lfU<6gS!&ygA6g6`4Y~6AM}9!u%2r`u-Nx%4o5uzBDDO;966^S%hl(C zccJ{EAA5{jw>n1mXX`nS%#c4MVP!)ckptKToqj5K(l!lZJvA&qJ~W7#{j)$krpfm? zPWVAtK8Q7M2`@uKf0;9dUA=cflhX#GklCVjJ}sE~^K#HvTP!lY-Jc%r=vU z*&66ETn3q0I1isDA2PFkJiA;zVP?HjTdk)HL%w^3c5ytRcSjkPrBqMFY)gv%(&aZT zKA+#D!WlkP)5gCaV=IdGhc(n|c((l8%pyIXrTm4x$HtkmQ!4)GC{X66GSisu!9)#w zA3-M_S%0(yqy)A0@~0X^wbn0Etv9JMq!mq?@yQbXZ3`r;BeJ*0sR! z&|~mAc`udq3oN8X#}yZW?1~>Sa8w<$%Zwo`xPyS9*v5J%9f%{k-5WV5y|fdGgezX{ z3=Xaf9Q>hNKZJ$yDR%kx5Y~A@530v?J(bdKBmL?Eze0G0{r+xU;2vOD|8Eg$!E2Er z)15H>w)RQ_ME`|*j%}<5SXP4tSJ~iTbr_j7a@TgBQIlowp{(u1y7TSL<`tjnN*+WNkB zm;&jBi=Rz@b&c$i2E#e9MvhNoUHG??r7ey1WxvZ)X>1-_ArptO(JWT3AI94AO%vq4 zVXW22J&4oP*D1|o_vuOai@-}7?rTnKs()7mOc_>1sSAAWkTx?T3J4POyWtPe_GzTP z)>H@zxNHPO+L8IFl zGB-X;35V;3t}=ku1zrG8%^EEQQ&WsvAa;^X-w#z-crPS0(IfVUf(%a_n)lV(tkUHA z_>nAJ8DkM}#^EAU<3?^c+^}gx-!4Ew^?S zppOLV13PZfFE;%mjWho{m0mTf`5LFW{$Xe1;^k;!(~m-_sFp1{0>d+{l2_7Md>fF= zC-gc1`kc@+7giTI;ef=SJ|#9mtWg)(?SPz>!9L*K#> z*VQtABrD)EGUXpA@~BMNV-#D*pB*cYjbiWfTcwT_quF7`SB#Vc#vCyM(XZPt2penMJuYywp7^aCnA1H;(OW`#5f6)fcs!XENUHpf_*Za*`73GcR4PAJ5wI zyma}-cy^hO9xkmD*u3c95%yLO$5qXB$b&cV5@?Tm;heo>SkWC-Ifk);XNJkJi7c9L z9wvuP#Hfc4liMb;juE%hph6`SiiC#mprE}+P4}M=b$k=;@NRoe)=XrR4S%E7u*wcr z3Qzu!nUh!>elkrip2WJhT%*Ye(<=SZ7&p!@;R|P_$%;v=M~CkBtIjYG`IntvOpvcnW;QRQMN^92 zXpsY^u&$o!P-j_{PflSc`QV{)a28AI@ZFF`bNdm6MswSY{wWZs93tPyVnMCSP+jx9 zbC=O=8K=Q|s^2_BUddu%eD)B@r!q%CO{$_)g-~tzX*ab&u@P3A-=sI#a`sw*t zMI5YAGJP7ZGWoK28k}gG!E((!%w2RtVWi`)g{%p8Um1_WAV-tMY!h?;sS65wJG6b*n{&ofG zp+tCggFjM=FuvRHzGt4jhYXIs93Fu_2LB#=fU4gw2ifkk47DIV3KE!>x|q}x&f4^~ zT5Y1mS>3O`dN-=kila;~5Qjj}ZJXB3-@~2VkbACP4z`k#YMTa1>OYRDm)JHY>sB#; zmL>gHvp6S9(SIxmad{bTB@TcqQ z2IAWb-xz#f+E`a7BW+C8S7o4;oVtp&kx6S=3w_C9S2^5KxR!0_8KzNms5E@ywmZA6 zrz}Z-u-+QpE$hX$B;RU(>V$1tidC9-Zh5sri zxO&ff*C{^B1?zjRuP=1bLmsU6QsP_!wmYmPOz4M&H0MgzyG$`bLL~5zM!3m^2ZL^U z$I0UdJkEt91lB%f>z#v1FOf4hu|YoC8bDqdtp%;%H0EHN*-%|U=we0r>H<^Z0Ixl& zD0hM11Z=TXQQSUdN(kVnLPcodX{3%nALBk0Ko?gw}d z3C9&g2h2v?w;gyMU?@(X{h#6wo#KGh;(G=a3WDZ0*DK0*!1n>}!J7+i7l8wQfbS^a zWq`e(QIt&JCct#5D6@eRJii&80j~lK-l8a5frkM4ZbfIn{Q%2fKwrQs0Yi4+nap>@ zv;*npor>~1@GQXVFQMblpb0?7tBTSbcp2dQ-HOr&cp>1$*U&NWD!_oZ;E{oc0H)x( z4|qD@ZKT8Y16TGd%3g<}dqIX=Kp26q+vivUz|RgViU;s2z-8qK%7GUG zdLLDkX~6ved*HhixC!v~aYZQvPX9;%JEbT`i3c#_G-e4r3y__`5)&Qp9=_L!|6`o9 z=g|LM6!L(4d;!*W32Fg+{t}!K@a=%nUt(Q>#{pjaPEj@kuLA7;17-o-1UN@U+slvx z*zG4pF#$ILj=zOv0-go9zY|mT11|$?9gD$y32nt<{ZI5H;piCfWxQgY4tyWr6?_i^ zuL5*4Gvz37KfuewnbQ3$@TW5+0DoF(x&lFf8}T&*-wya8zN3Jb0bay66L=Ni1pIN@ zMBrI~+4$yODZr?KoW*xGDk=fHj$+CZ;N+Gn$1-IFH3l4#$&{_Y(*dvGyO-#IGsau`d58y$3Zvrm^Ts)a6zXPuX{0!e(;8lPjQbkmf@QPyby5h9ISsX3fob*l!qk%o<5%`Z5QA#1&9u~ z6W_>6j2h4rf1~Gy+4up@!1wUi&=BA$e2)UJ1iY{a>kYgL(7c!_(|~6Ic3*;az)gU= z@ICq!G`FOHDK9*ZhA3mqSy);{e|)f(+ngfXj+8LE-_NzKSW{1(;1Ckd~`4GE{^BcDKT4fSUjp zuf^<%2e9ip7!J_^N3VyrfM)^D#dkXfL2xO)F99dG2H#!4O8^hxyA=FofWz=s!wTT( zfbOI5v-BHTo}z0Z^DeQS@9tm`dY10!@G6pNEM1P-#V!Cu?`Dqz?AVPPbh-@KgAgEH zF5AN%N=;|@8#vbUZGI5b+~xi(6U~%=n&|%sP=52!9}iG&d+N&rl+{fmmIo*+{fU0v ze>6OOGd&m;r6u0Nce@58UGmr1L^=lFZj9&214a&0~`hN_{ z%bxn5e3eb!B|1epWYDAfxIvHF^BP`DxIFzD>s{cKvzcFGuVQfP-Soxtl|yc&P5;SL zYA}>Mes#^;ycVLHa9s%B2?ipXVwyt=F=` z1$g`BCl9oKF-or}rBSniUXRiLkf(gxsgR>k7fXfRUG&BI%8~d2{hRZYSG(%>%u`l( z)9;w8-0G%R7bxHK(AUma4)@ewT%f$tOTT5IvOht;El;r~>bK-6ANAhn`%IqlVXFSr zJZ1S%{l^QHs-gOA^OX{de&a&r`LxTrm*y)MCh0dXP>xL2KfgfP{_rmhT%#~cjv`MEU5nmj$mcQX%!E%OOIxlq4NhQH3{Vj(xa&awe+zpkxr&Kqos zy!-~MY({G{00)StPVvm;ird3T=9v%v;5R-RNKNksfcekODeu>3pWM zn9gQ81n4X#ryo)U%f(?q{+td9I@8IYQ;=5H^6izGKm?vdOgxkGY; zGAs1YHZAev7@YIt-HI!#K|9%qmJiw`&f8?1IN# zu+jxDxnPwG{^WwcyWl++^ayeeJV1lth`C*aFc+*Jd7P^}N_0f+e_JILvrNtf&rYwfG-2-y+2BW_2+4!@T}UX~p$v5=t1SK?kV zy<&bP{gtd&a$hm+HSaCiyLoSYn&UD{@~m{F5N>P;Z6UM1U~!(>O&BY;lPLSXVCi0Y zBU~4F{f5vW8S*9b_Z;h_Y=l1wlDL?}yNw4$_8+AS+z{t3w}awI6v+8OM?5#t%55_7 zln-(8gnY&PJ%|Dx6ZMr2g69+`q1v4q6%tF1lj4fR_7Uk3C+Tx)g@DuFBiBiYAp$~K ze~%eXLWDDXL;UFP;Y#7WKstQ08ir4q zvSjv@$Cc=q>*&B!w_V11@m9ehZ#z5mLG`pLQGoy&*T$EIQshV@hrGLbwq+^V2oaIICnjDh4M{b9uZuhc%s^D)OzC*gxWmEbzk0s z@w^cZy&oUx%T1%?_;%dF-ybD+wc~l*Z>((9p7$d-zCG{RvbD=fUkXH?!VlA6BRG6g z(@rukl(&{g+w)HG6)qm17FxZhF1OUFAI>s~S&8U`))Pj%h!aN3&=78FT2F8TVRi`b z*(%4yL49@bZ|LjM(Q5zafe=s(7(r~W?F%{cxf zs)7neIX(>M8yT-0AxB2Q4EBwbuSf9Mf!~d6;GPoXY_SjxesR&GA*JbESKoCh*LbnMmakrc{N%!rn9|XX49^X zsrK7(}8zS*y^(BX}*ZPoi@X|DxR;o$SEM(G;)TgE~{Ru=b;fYv?CAX zCr8MHj=T*&J3@}>$e(0m9iMjOr`*{%`9dt84De4ZPweG#oUgenice=(AsYo9Y2I|z zIesJxS9>VR>*Jh7_WC$Erweb>%vGM^1=WWUN7*>JwF__A#wAFraPR9XjiEl~jF+Eu z;mcU2OpoK47{%^5-p&B;`==i%kf2Jv$BsZA29iV*@j|Lc$#1p+Hxkmlh zBxz6L^FaP5iDv;!9-uK<2k=CcDhBYexNf%|$Rhxf2J(&|&KSstp>$v%pAQh7%sYBc zbge^oxnTfrF6SomUZ{E{8SZnU{H8HSOD9Ko3h%=v$gC8e1h69oolcM+r=ZgbvdJJm z5n$RNz7XKjARg^C!6n~3L3$772`yda<4qOilK|`zsALgTv|)mrG?;%3O~jb-?=y^- zGtH1cUhXyX57=bKtW-Xl;VDM>5dJ=j1w(l=Dhdn#l1*{UOXIf~?z9gN=ZBk4b?Hah zFjYpSV{}vGsC1r!Cm~y$pemhrXqr{em4_!J8N5?Nv$Cnu3<7&tF3R8?Shn1n!DE_c z*Yg)<%kvpL65{^Oz|v*Owj+3pjxKR#98tWbYXo!?T8)E_kzFx{om3_#jDX~+a`gxv z+QG%v8OD7?W6R(h@b?cbA9i-7$->)v=SJ{gRwh+y`LGNe$-@8=MxyU*Io=7D5;0qD z8;NrzTRNP?%LJy$zo~ATY&QxA;WVdY>Twv19%G^&`Fqpkuu&Sx0nJc6sKE<~P zAwlyTId2?KW{=4I<9MuZV>!W($eZJ!+c-IOJdcpE> zfv;k@vT_2y2=K~8{tlbrm^g_SGd4qBo6OSzx=lgXGvtyfyo0YRmiq$ITJD-5_fA3A zd9rO5y8dYjx^7h8s+B9s+W}^{Zx+wVeDWON_GfISqsvr?XESBt!~9`@Zy)9Z0Xk;$ zE&yY*c^Pb*9wkz){b16y@+2Ih!RTX23BaF$oz38w_uFjXYn;`j^o2w{7uH@$a!<%wB|@ThsOdO zpTi#mFwNzK;3=OA?W35Q$9Di+&x4}p%J}*4c5~&}`P||)*EtP8jU;*!ZzkGR~+`WKD`MUUN+1f6fC#x3l<_ItUSit|m>ixEm zuVnKbkLSa`u=%oV5zhb!Sj_u+x#DN`;e0u6G4E3!K0hlbXY*KjdNFqM0!Qs){zW7O zzn^a7d5j{_K`-(ciaZy-$Zd#3-`#;|mO{ngop9L*6&<#hI8v%(9WU?VTR0q~={3HW zjd%R~8tdJEJkj)b@QuX^mQ=mT6b*zA#mP+qPac!ejeF7;+a8%Vr4ru|SG>PKW}L?1()k2$&9`UBAt!ijX#E9^;x77F@Xdslt6XiK zaJG5<1n=gV<%&ymGUEhq?+;T`Zi2P`MmJf$^)ZP*$(vGi*83y^Q$*E|pX85vxI(^} zGX5jpPPRA&d(V`;Pr<--KV$u*~7*yH6#r+AlmS3uE;Mubmud@o@$l)320&S{@y z7u|1yjQ)rZE^wV{h0fah&Ej0AYp;LN^IY^N{zbp(qWAn4JJFRQt`6dp?GYS@PP)2r=f!-e-AxfbnP1_Z+$AEc%}BIDeL( zY0Gvwwq55Hp>nxF7osnNfJ#{&q-!r<@zn(}KY7(x*Iw_5g;WR)XE(K4I2)Pes|(W) zMNNsb!E6H>Y}au8OCawfw6)`iLH7$|`y6li>kuAPNx;^de03hoBx9QCj5sLMn?dfw z(v+ra=h)axH-|+aABs$;frftX-9 z4T4sKl80<|1LXry;&5$gp!@}jvaB&B)(?8b-O$N1a?u=Za{E@)Q`@!Q{+XDcYC2PkiXQikUx z41Ma@2TWGA(lsCU>3wSb=Kc*7#nD)nR|`nT`y36r4*@0Xy~dOTP%5)!NjqJ5o9xD` zksplVPp|JbS8njt1x9n0sKTkRFRy-DN9Fo9y1*H&+SL=O{le+_Bd1TFrc`}eFLC(HshYPieO6DQQ4K?@ zyi1KKlR+_k(U>v|6hBR`w0Pq*Ey}zQoxa%v%DzmZvgPU!Y?cQq2DgXFER`RiqRrCA z^SFwNeGk{q!zs8=yO4Tfy>`l69@|h-QOQ)~2zjb)Xb3lEXrE;cf9uu0dD(LouDl z>ig)a(=0D|a9BOHzLQTu327$}h3b0r?9f=i!f=cSPevO=eh5m{1C&oeF=>aQG==Fr z@y`G>Vy{7E*#oDXI}ZEpa!wcy$1J%4K%a+1g-%ZGplc@2gz3Clh`da-h4LQ)Au=>v z7pw1scSqdhzDc^~@=!0Gmn;m>8Re#M-E93rZ`{nJe}XPhwlV4=ytLkEv7)tXy38`_ z+;D*RjnK8&5bGgV7AnYOm5cl(@GK|UyCa;Ti`*CZ zY^U70KovQu0|v1T=n@$psS9xf=I11OMdHmfCu=g$5GS<|Xo-`07-*c6dJX7SC)Fd0 zmPHdg2Iwx0>i!bYbRSt3pbK(N@FNiRJL_%&&GeD|A~2hc(d4E{9q!#@@P1h1CZ|C- z%@kyjJhgrDYxo_Oz$Mx(^$PhXM`EZ8KQcL&ck;@Ltx zx&Pq_G0C~a(+5K3@5E#I4^N~?o+O@4#4`*(&6EGb6K9fPJ-~C5cyJBnn?|g{x2FF^^$j~!BJ{Je+&2juZJ=(IVK(r z+QsVvTI~Iw8_1}IxIvVGLoOk;$oZn3sE_;89yz)a{u>S zHfwC}%>~;}#I}Rj4*Z90yT;a+*phpI?JlwH{SVtNjqSrZV0(es#^9$#$$!}PX>3Wv zR!eOA=%-CjT|lz~RP%t&R?DQG(3ux%1N8=F}h%Rcc9K&c1Y0Og;Rqg3zvfuaV33R&Ps&zvLr`IqHX{~ zVuhTXc4>WE{@Gir(yps5wxE})2m0L&iXZhGp0Eq&txnU4_5k#8LHga1ZT@ySweU(W24P_u=0)yFE`6K3t>8!`}WNd$(8MxMbU{bF9us_<3FC-X_KNknR zH4Yq^K3V4`-A%gvtg=B3*W@>*2PDNxwU;j7KL6Y#-F~ik$nE`ftu(#)%fSPnx4Z%J zh0)qU=KdvSUpYXQ4A6Dwret}1fG!N9q11zZWDkUCPP7nc76mq1+c5)mGg+1_ABetZ z%PGTj;SP^vT{7cwT32B*F9pi#C7&TMxV|fSDMc5>%#shn?GJCeH5f_9o(M)$>nZgs zP&HT%8Kj%3ryPP_-pGWVv`^4AaaFV)m0la`WPyc2n>-hlrVYrA;epwl3?3a zlXM=kodv}`PwDX50_m{Q44srT(n+q|ovJk^8*Y=6w_rNraiws`zsFNJgi<6y@|U5y z-twd6x>oF{oHs=GH!pJyW!+GSSt?5jxQ0@hfv67yxWsQZ>hYj#omSstc)|US$Qb8C z6@JL(G$>~!DEl6u>;xsw(wOHxQs4uWv!L93fO3Uer8RE#Hz;Lk0rjMMy_XD;zo%h= zPdc!&k1zzZ@uOj5zF787?{nB6q`E1om9*_HSQVGEB})*Ue(X zoXe!?W?woMr$YvmnJ$wuU{4q2bm$B`nn%$Bn#y7FjSLK^med_C`;5?CLd>fTDq$pe zf0L64G-}ecs(zqwC|8G~(O=Va{rFv93KoN$?xy=FxSM?>W-tqukBx?e;bf}e&4&M5 zBO$}(r=u|d57~bVxXT=?#^?r85=|S!FJsZPpKO*1FvdCT5t+KqzzZ|cC{EU7VqR?= zhH=m)_w%CN7$&!j2Ya^skib~Yixu?hg_TBLmo)6I#4ebgKV{y6N6PL~$UaiQlXssI zo4e={^21rswbRALK0bZHJWz5QDAVRalA=s)Xf->YL@5&*D8$nuzpJwSK8YAI(Akkp zxp{&P53Xg?V{=e@V@?eXmlS@M6sbMCq4 ze!jaO$%{YpEqOi>tC!3fdE&5N79^+p#B0sEMvPaQ#~TSfCz)|VFHDm| zbDLyF$$5Sfn{etYct7M}#G~rW zjVY*;LY)p#qaMYAE>z=IsIws|?@^n~4O^kcg^C6nLhnv-%%#ljtwj0?jJdrV@WK!{ z?tue)V7>uk-V|_d2%PZ1H^t-20@mKak5sokGwSU^A2XRd>}9jZIEd33oWwS5x0|mP_mY7; z%`$ssC{I%!Rbq|^^-Ws2d{Kz9%%xlGWeoJow*eg&s3HVK9Vqfisd;!CGvx-lV0l@H zDs!o+QnRoRYE-DH(Arz!Qm0GJFN9hrW8(xrw=XKctuPN4 zCdDe=CA6Jcmv^;%n?@Mh&VC>+Gw%X3^*xL?nlJ7s`0$pJ+-;_L2PQ|&t2@{Vw>#_5 zf}x8uvK>gcw%II3UFIQqO1R}a`_5BYUQ$b#+MV`^Y4@3{cX7GFU!_X<|G+M8H`3nH z7QKt5jT;RooYeP$kD4FzbeVSoJHu*k?Pz6Lw<-grYx6bjC~@~iJZHOm|87ccG@ZNc z*Lh++nQvNT$SvAuUtqdg$=8GYTK=7&TY0HIVg9&>QnSwP-Cam?I7l;J=62e#S#oAp zopm=W_O$8dxz+ZZ>CPdU>=5LZC35YR=GnXLMQr(lgQ^j!(IqC9u%oxinU}k3A!zHu zmF!q~tCLou>Q+sd*&!tsQkIzAd)ae2ThV3i{p&I_xYxd5&WgW$cEh?f6)|l~=&rxJ z#-4el_Xboe<@uz8^0<<6E2;_|NR%6$XwAWWybu&4WSMwewvQ`Lcl4EmC{^BX5v#s- zQ+su9zqL~Ndl*I@;=}dx6iSUBMSX4OwtlWvDy!#Vq8U|u9%)c1|3az6Gzf9#wJ1JQ z-`>@~cmbYN?%pFxqo#~M`QWfXc$#&aZQZzCnxKY#jZ!R$j6 zi&sbkW!6`Pbg#PsU4~MvY9^n)B7h>R4JgE1q%@;Q(3iUP5Q9FE@ae|r)NB)<6E^St z!2VL^AQ!9sY771dihE)c^U^nrBA;!g3ZC);N}o#2>nIa`A>wBOMW=5{{3vF6<~ok4 ztBRX=BT7~&+fd};wkM6=l`A8ubczkpnU0}*KGw!kH zWp4QzH=e5CjVP9C?xZlPx%Z>QR7d(PiWml6;<)_EgKdOY@Ea)3S(YEl#~GxRDy8#Q z!K+TP9OaOz)6FOo%5@XUN)=)!N|v&N9&tBHPGu{%iUBjO96b#*rCNOyWmttchBBy> z_faZUJ}$l6*-{VIDuJ8MZI5KtKts}W6s=5&wh^o z$|D$KCr~O?d;SyUkYB9O=ia~%m8^?VWP=l=Z!t=hihqSzsiwLmgR~G;dYaqYk|+t^ zbyGt-Qyj&sI`+cEz=OiFD7ub6La6}=I@r%p3MU(GS&!b3;UkI)8U`vXWf6-LrPj`L zvz1fDpMa@WRwq#k8&AHEsXduYQK6{HVx-D_lT{GVqADi9taIE^7;DRba#nFZ)}Zte zZeaFZ6g>u#D1)kjcA(6*imeWzOnG}_xAwnA$&4^vxg(A#Pd?pYWmQ*s1}LFiM^}@4 zt2o3Mik|1c9|ZjOB@2NZJzC910N zGRj_AfFndM|19f$q&^?&PP$!xR7E}GM)IKo`T~j$XbDO}6<>vtS5})*T2L}U?bAqo zK2+b*pXy?EX;H=RfRSi}+2THw3ExVp@i&wlOyD^@`+p2mCRNQ?F|YRLtR-sYXbwt+ zKZID(?e9X3K!sNaT}8eEN*-=K*!nfd z*00KMAh~t~rpiyX%P^xT=4lu`fQx4K^C;fgYcTL~=BEdlk^dw<)!6wErN-~h>4r`z z+8S4hD7lF~hbpl~>-qzf1J2D< z%q#3Hiu}RypswewgHg&?P)3yJ_3IdshkVFsnI_Cs8f!6Q6)V~~Z$=qajs{SqvO)3t zQBESx@B-&EnPZx5M!?dm$8n?@)t^qHl&K7zMoFnOo<)f%tBbx)v-(!`>DG?q z@}*YRu_~0c7$&7}Nmd}mW)!cq%_ud#BjU7ni8K6-+`^tQBR_CjU|2^`l67 zYX3b-m5TZ$6e%RA>w75j2$VA@6~0w{vb(#JgZ_kSvY9tyLSBnzv vF<#<~kJz6dnmlGd9yNCzvzHDX8MWu0J9OqXyR*c+Gj7ix`t$4d@z4Dq%?J+a From f920b8f6c9d03b073ed4aa68313cf6c5df558d90 Mon Sep 17 00:00:00 2001 From: Leonardo Leone Date: Wed, 27 May 2026 11:03:11 +0200 Subject: [PATCH 2/8] Fix kernel kill in isInConvexes --- depth/model/multivariate/import_CDLL.py | 2 +- depth/model/multivariate/isInConvexes.py | 92 ++++++++++++++--------- depth/src/ddalpha.dll | Bin 255383 -> 255383 bytes 3 files changed, 59 insertions(+), 35 deletions(-) diff --git a/depth/model/multivariate/import_CDLL.py b/depth/model/multivariate/import_CDLL.py index b052498..25c403b 100644 --- a/depth/model/multivariate/import_CDLL.py +++ b/depth/model/multivariate/import_CDLL.py @@ -32,7 +32,7 @@ def import_CDLL(): libACA=ct.CDLL(ACA_approx[0]) if sys.platform=='win32': - site_packages = [p for p in sys.path if ('site-packages' in p) or ("dist-packages" in p)] #Add search dist-packages + site_packages = [p for p in sys.path if ('site-packages' in p) or ("dist-packages" in p)] #Add search dist-packages for i in site_packages: dll_path = os.path.join(i, 'depth', 'src') if os.path.isdir(dll_path): #check it is a real dir diff --git a/depth/model/multivariate/isInConvexes.py b/depth/model/multivariate/isInConvexes.py index 0b1d50b..63b7814 100644 --- a/depth/model/multivariate/isInConvexes.py +++ b/depth/model/multivariate/isInConvexes.py @@ -12,44 +12,68 @@ def IsInConvexes(X,z,distributions,seed): try: n, d = X.shape except ValueError: - n = X.shape[0] - d = 1 + n, d = X.shape[0], 1 n_z = z.shape[0] - points_list=X.flatten() - points=(c_double*len(points_list))(*points_list) - objects_list=z.flatten() - objects=(c_double*len(objects_list))(*objects_list) - distrSeq_list=distributions.flatten() - distrSeq=(c_int*len(distrSeq_list))(*distrSeq_list) - points=pointer(points) - objects=pointer(objects) - distrSeq=pointer(distrSeq) - - distr=np.unique(distributions,return_counts=True)[1] - distribution_list=distr.flatten() - distribution=(c_int*len(distribution_list))(*distribution_list) - distribution=pointer(distribution) - - CSum=np.zeros(distr.shape,dtype=int) - CSum[1:]=distr.cumsum(dtype=int)[:-1] - cumSum_list=CSum.flatten() - cumSum=(c_int*len(cumSum_list))(*cumSum_list) - cumSum=pointer(cumSum) - numPoints=pointer(c_int(n)) - numObjects=pointer(c_int(n_z)) - dimension=pointer(c_int(d)) - seed=pointer((c_int(seed))) - numClasses=pointer(c_int(distr.shape[0])) - belongs=pointer((c_int*len(z))(*np.zeros(distr.shape[0],dtype=int))) - libExact.IsInConvexes(points,dimension,distribution,numClasses, - objects,numObjects,seed,belongs,cumSum, distrSeq,) - res=np.zeros((distr.shape[0],len(z))) - for i in range(distr.shape[0]): - for j in range(len(z)): - res[i][j]=belongs[i][j] + + + distr_uniques, counts= np.unique(distributions, return_counts=True) + numClasses=int(counts.shape[0]) + + cumsum_arr= np.zeros(numClasses, dtype=np.int32) + cumsum_arr[1:]=counts.cumsum()[:-1] + + c_points=(c_double*X.size)(*X.flatten().astype(np.float64)) + c_objects=(c_double*z.size)(*z.flatten().astype(np.float64)) + c_distrSeq=(c_int*len(distributions))(*distributions.flatten().astype(np.int32)) + c_cardin=(c_int*numClasses)(*counts.astype(np.int32)) + c_cumsum=(c_int*numClasses)(*cumsum_arr) + + output_size= numClasses*n_z + c_belongs= (c_int*output_size)(*([0]*output_size)) + + c_dim=c_int(d) + c_numClasses=c_int(numClasses) + c_numObjects=c_int(n_z) + c_seed=c_int(int(seed)) + + libExact.IsInConvexes.restype=None + libExact.IsInConvexes.argtypes=[ + POINTER(c_double), # points + POINTER(c_int), # dimension + POINTER(c_int), # cardinalities + POINTER(c_int), # number distributions + POINTER(c_double), # objects + POINTER(c_int), # number objects + POINTER(c_int), # seed + POINTER(c_int), # output + POINTER(c_int), # cumulative sum + POINTER(c_int), # distribution seq + ] + + libExact.IsInConvexes( + c_points, + byref(c_dim), + c_cardin, + byref(c_numClasses), + c_objects, + byref(c_numObjects), + byref(c_seed), + c_belongs, + c_cumsum, + c_distrSeq, + ) + + + res=np.zeros((n_z,numClasses), dtype=np.int32) + + for i in range(n_z): + for j in range(numClasses): + res[i,j]=c_belongs[numClasses*i+j] + return res + \ No newline at end of file diff --git a/depth/src/ddalpha.dll b/depth/src/ddalpha.dll index c759f5286019c7eb93b9970a1a07b772de5d871e..bfc6886cfb2a499c4d473c8565d7b3a2fd747c01 100644 GIT binary patch delta 43 wcmbP!nt%Ff{s|q-RvX18cKb3t323~v^%kSnMvy=|`zA&pX4=laiFtxC0Fz}AApigX delta 43 wcmbP!nt%Ff{s|q--?Ky~cKb3#A8EX`^%kSnMvy=|`zA&pX4=laiFtxC0H~W0_5c6? From 71df6334837394b16dba161881dde923d7bc1f17 Mon Sep 17 00:00:00 2001 From: Leonardo Leone Date: Fri, 5 Jun 2026 17:07:00 +0200 Subject: [PATCH 3/8] Fix errors in MCD and implement it into notions --- depth/model/DepthEucl.py | 173 +++++++++++------- depth/model/multivariate/ACA_wrapper.py | 15 +- depth/model/multivariate/BetaSkeleton.py | 17 +- .../model/multivariate/Depth_approximation.py | 2 +- depth/model/multivariate/L2.py | 16 +- depth/model/multivariate/MCD.py | 105 +++++++++-- depth/model/multivariate/Mahalanobis.py | 13 +- depth/model/multivariate/SimplicialVolume.py | 15 +- depth/model/multivariate/Spatial.py | 14 +- depth/src/MCD.cpp | 143 ++++++++------- depth/src/MCD.h | 22 ++- depth/src/ddalpha.cpp | 48 +---- 12 files changed, 335 insertions(+), 248 deletions(-) diff --git a/depth/model/DepthEucl.py b/depth/model/DepthEucl.py index 3aacb39..336dd0c 100644 --- a/depth/model/DepthEucl.py +++ b/depth/model/DepthEucl.py @@ -129,23 +129,24 @@ def load_dataset(self,data:np.ndarray=None,distribution:np.ndarray=None, CUDA:bo Not used, present for API consistency by convention. Returns - --------- + ------- self : DepthEucl model object. Returns the instance itself. """ if type(data)==None: raise Exception("You must load a dataset") assert(type(data)==np.ndarray), "The dataset must be a numpy array" - self._nSamples=data.shape[0] # define dataset size - n - self._spaceDim=data.shape[1] # define space dimension - d + # define dataset size - n if type(distribution)!=type(None): if distribution.shape[0]!=data.shape[0]: raise Exception(f"distribution and dataset must have same length, {distribution.shape[0]}!={data.shape[0]}") self.distribution=distribution # define distributions - self.distRef=np.unique(distribution) # define unique dist else: self.distribution=np.repeat(0,data.shape[0]) - self.distRef=np.array([0]) # define unique dist + + self.distRef,self._nSamples=np.unique(self.distribution,return_counts=True) # define unique dist + self._spaceDim=data.shape[1] # define space dimension - d + if type(y)!=type(None): if y.shape[0]!=data.shape[0]: @@ -182,7 +183,7 @@ def mahalanobis(self, x: np.ndarray = None, exact: bool = True, mah_estimate: Li cap_size=1, start="mean", space= "sphere", line_solver="goldensection", bound_gc= True, output_option:Literal["lowest_depth","final_depth_dir", - "all_depth","all_depth_directions"]="lowest_depth", evaluate_dataset:bool=False): + "all_depth","all_depth_directions"]="lowest_depth", evaluate_dataset:bool=False,): """ Mahalanobis depth @@ -198,7 +199,7 @@ def mahalanobis(self, x: np.ndarray = None, exact: bool = True, mah_estimate: Li - ``"all_depth_directions"`` : tuple of numpy arrays Returns - ---------- + ------- array_like or tuple of array_like The first return is the lowest comuted depth regarding all explored directions in space. The second return is the direction that best represents the analyzed point, the direction corresponfing to the lowest depth. @@ -243,6 +244,12 @@ def mahalanobis(self, x: np.ndarray = None, exact: bool = True, mah_estimate: Li alpha_Dirichlet=alpha_Dirichlet, cooling_factor=cooling_factor, cap_size=cap_size, output_option=output_option, ) # check if parameters are valid + if mah_estimate.lower()=="mcd": + if type(self.MCD)==type(None): + self.computeMCD(h=mah_parMcd) + + + option=self._determine_option(x,NRandom,output_option,exact=exact) # determine option number if option>=2: if evaluate_dataset:self.mahalanobisDirDS=np.empty((self.distRef.shape[0],x.shape[0],x.shape[1])) @@ -253,6 +260,8 @@ def mahalanobis(self, x: np.ndarray = None, exact: bool = True, mah_estimate: Li self.allDirections=np.empty((self.distRef.shape[0],x.shape[0],NRandom,x.shape[1])) for ind, d in enumerate(self.distRef): + try:covMCD=self.MCD[ind] + except:covMCD=None DM=mtv.mahalanobis( x,self.data[self.distribution==d],exact,mah_estimate.lower(),mah_parMcd, solver=solver, NRandom=NRandom, @@ -260,6 +269,7 @@ def mahalanobis(self, x: np.ndarray = None, exact: bool = True, mah_estimate: Li alpha_Dirichlet=alpha_Dirichlet, cooling_factor=cooling_factor, cap_size=cap_size, start=start, space=space, line_solver=line_solver, bound_gc=bound_gc,option=option, + covMCD=covMCD ) #compute depth value if evaluate_dataset==False: if exact or option==1:self.mahalanobisDepth[ind]=DM # assign value - exact or option 1 @@ -313,7 +323,7 @@ def aprojection(self,x:np.ndarray=None,solver: str = "neldermead", NRandom: int - ``"all_depth_directions"`` : tuple of numpy arrays Returns - ---------- + ------- array_like or tuple of array_like The first return is the lowest comuted depth regarding all explored directions in space. The second return is the direction that best represents the analyzed point, the direction corresponfing to the lowest depth. @@ -410,7 +420,7 @@ def betaSkeleton(self,x:np.ndarray=None, beta:int=2,distance: str = "Lp", Samples matrix to compute depth Results - ---------- + ------- Beta-skeleton depth : array_like """ @@ -421,10 +431,16 @@ def betaSkeleton(self,x:np.ndarray=None, beta:int=2,distance: str = "Lp", else: self.betaSkeletonDepth=np.empty((self.distRef.shape[0],x.shape[0])) self._check_variables(x=x,mah_estimate=mah_estimate,mah_parMcd=mah_parMcd) #check validity - + + if mah_estimate.lower()=="mcd": + if type(self.MCD)==type(None): + self.computeMCD(h=mah_parMcd) for ind,d in enumerate(self.distRef): + try:covMCD=self.MCD[ind] + except:covMCD=None DB=mtv.betaSkeleton(x=x,data=self.data[self.distribution==d],beta=beta,distance=distance, Lp_p=Lp_p, - mah_estimate=mah_estimate,mah_parMcd=mah_parMcd) # compute depth + mah_estimate=mah_estimate,mah_parMcd=mah_parMcd, + covMCD=covMCD) # compute depth if evaluate_dataset==False: self.betaSkeletonDepth[ind]=DB if evaluate_dataset==True: self.betaSkeletonDepthDS[ind]=DB if self.distRef.shape[0]==1: @@ -449,7 +465,7 @@ def cexpchull(self,x: np.ndarray=None,solver:str= "neldermead",NRandom:int = 100 Samples matrix to compute depth Results - ---------- + ------- Continuous explected convex hull depth : array_like """ @@ -532,7 +548,7 @@ def cexpchullstar(self,x: np.ndarray=None, solver: str = "neldermead", NRandom: - ``"all_depth_directions"`` : tuple of numpy arrays Returns - ---------- + ------- array_like or tuple of array_like The first return is the lowest comuted depth regarding all explored directions in space. The second return is the direction that best represents the analyzed point, the direction corresponfing to the lowest depth. @@ -632,7 +648,7 @@ def geometrical(self,x:np.ndarray=None,solver: str = "neldermead", NRandom: int - ``"all_depth_directions"`` : tuple of numpy arrays Returns - ---------- + ------- array_like or tuple of array_like The first return is the lowest comuted depth regarding all explored directions in space. The second return is the direction that best represents the analyzed point, the direction corresponfing to the lowest depth. @@ -734,7 +750,7 @@ def halfspace(self, x:np.ndarray=None,exact: bool = True,method: str = "recursiv - ``"all_depth_directions"`` : tuple of numpy arrays Returns - ---------- + ------- array_like or tuple of array_like The first return is the lowest comuted depth regarding all explored directions in space. The second return is the direction that best represents the analyzed point, the direction corresponfing to the lowest depth. @@ -835,7 +851,7 @@ def L2(self,x: np.ndarray=None, mah_estimate: str = 'moment', mah_parMcd: float Samples matrix to compute depth Results - ---------- + ------- L2 depth : array_like """ if evaluate_dataset==True: # Dataset evaluation @@ -845,9 +861,16 @@ def L2(self,x: np.ndarray=None, mah_estimate: str = 'moment', mah_parMcd: float else: # create self self.L2Depth=np.zeros((self.distRef.shape[0], x.shape[0])) self._check_variables(x=x,mah_estimate=mah_estimate, mah_parMcd=mah_parMcd) # check if parameters are valid + + if mah_estimate.lower()=="mcd": + if type(self.MCD)==type(None): + self.computeMCD(h=mah_parMcd) for ind,d in enumerate(self.distRef): # run distributions + try:covMCD=self.MCD[ind] + except:covMCD=None DL2=mtv.L2(x=x,data=self.data[self.distribution==d], - mah_estimate=mah_estimate, mah_parMcd=mah_parMcd) + mah_estimate=mah_estimate, mah_parMcd=mah_parMcd, + covMCD=covMCD) if evaluate_dataset:self.L2DepthDS[ind]=DL2 else:self.L2Depth[ind]=DL2 if self.distRef.shape[0]==1: # Fix size @@ -880,7 +903,7 @@ def potential(self,x:np.ndarray=None, pretransform: str = "1Mom", kernel: str = the single bandwidth parameter of the kernel. If ``0`` - the Scott`s rule of thumb is used. Results - ---------- + ------- Potential depth : array_like """ if evaluate_dataset==True: # Dataset evaluation @@ -923,7 +946,7 @@ def projection(self,x:np.ndarray=None,solver: str = "neldermead",NRandom: int = - ``"all_depth_directions"`` : tuple of numpy arrays Returns - ---------- + ------- array_like or tuple of array_like The first return is the lowest comuted depth regarding all explored directions in space. The second return is the direction that best represents the analyzed point, the direction corresponfing to the lowest depth. @@ -1022,7 +1045,7 @@ def qhpeeling(self,x:np.ndarray=None, evaluate_dataset:bool=False): Samples matrix to compute depth Results - ---------- + ------- Convex hull peeling depth : array_like """ @@ -1058,7 +1081,7 @@ def simplicial(self,x:np.ndarray=None,exact:bool=True,k:float=0.05,evaluate_data but the calculation precision stays approximately the same. Results - ---------- + ------- Simplicial depth : array_like """ @@ -1099,7 +1122,7 @@ def simplicialVolume(self,x:np.ndarray=None,exact: bool = True, k: float = 0.05, but the calculation precision stays approximately the same. Results - ---------- + ------- Simplicial volume depth : array_like """ @@ -1110,9 +1133,15 @@ def simplicialVolume(self,x:np.ndarray=None,exact: bool = True, k: float = 0.05, else: # create self self.simplicialVolumeDepth=np.zeros((self.distRef.shape[0], x.shape[0])) self._check_variables(x=x)# check if parameters are valid + if mah_estimate.lower()=="mcd": + if type(self.MCD)==type(None): + self.computeMCD(h=mah_parMCD) for ind,d in enumerate(self.distRef): + try:covMCD=self.MCD[ind] + except:covMCD=None DS=mtv.simplicialVolume(x=x,data=self.data[self.distribution==d], - exact=exact,k=k,mah_estimate=mah_estimate,mah_parMCD=mah_parMCD,seed=self.seed) + exact=exact,k=k,mah_estimate=mah_estimate,mah_parMCD=mah_parMCD,seed=self.seed, + covMCD=covMCD) if evaluate_dataset==True:self.simplicialVolumeDepthDS[ind]=DS elif evaluate_dataset==False:self.simplicialVolumeDepth[ind]=DS if self.distRef.shape[0]==1: # Fix size @@ -1131,7 +1160,7 @@ def spatial(self,x:np.ndarray=None,mah_estimate:str='moment',mah_parMcd:float=0. Samples matrix to compute depth Results - ---------- + ------- Spatial depth : array_like """ @@ -1142,8 +1171,14 @@ def spatial(self,x:np.ndarray=None,mah_estimate:str='moment',mah_parMcd:float=0. else: # create self self.spatialDepth=np.zeros((self.distRef.shape[0], x.shape[0])) self._check_variables(x=x,mah_estimate=mah_estimate,mah_parMcd=mah_parMcd) # check if parameters are valid + if mah_estimate.lower()=="mcd": + if type(self.MCD)==type(None): + self.computeMCD(h=mah_parMcd) for ind,d in enumerate(self.distRef): - DS=mtv.spatial(x,self.data[self.distribution==d],mah_estimate=mah_estimate,mah_parMcd=mah_parMcd) + try:covMCD=self.MCD[ind] + except:covMCD=None + DS=mtv.spatial(x,self.data[self.distribution==d],mah_estimate=mah_estimate,mah_parMcd=mah_parMcd, + covMCD=covMCD) if evaluate_dataset==False:self.spatialDepth[ind]=DS if evaluate_dataset==True:self.spatialDepthDS[ind]=DS if self.distRef.shape[0]==1: # Fix size @@ -1174,7 +1209,7 @@ def zonoid(self,x:np.ndarray=None, exact:bool=True, - ``"all_depth_directions"`` : tuple of numpy arrays Returns - ---------- + ------- array_like or tuple of array_like The first return is the lowest comuted depth regarding all explored directions in space. The second return is the direction that best represents the analyzed point, the direction corresponfing to the lowest depth. @@ -1264,7 +1299,6 @@ def zonoid(self,x:np.ndarray=None, exact:bool=True, def ACA(self,dim:int=2, sample_size: None = None, - sample: None = None, notion: str = "projection", solver: str = "neldermead", NRandom: int = 100, @@ -1295,26 +1329,29 @@ def ACA(self,dim:int=2, Chosen notion for depth computation Results - -------- + ------- ACA directions for dimensional reduction : array_like The return respresents directions that best represents anomalies in the dataset. """ - ACA_tab=mtv.ACA(X=self.data,dim=dim, - sample_size=sample_size, - sample=sample, - notion=notion, - solver=solver, - NRandom=NRandom, - n_refinements=n_refinements, - sphcap_shrink=sphcap_shrink, - alpha_Dirichlet=alpha_Dirichlet, - cooling_factor=cooling_factor, - cap_size=cap_size, - start=start, - space=space, - line_solver=line_solver, - bound_gc=bound_gc) + + ACA_tab=np.zeros((self.distRef.shape[0],self._spaceDim,dim)) + for ind, i in enumerate(self.distRef): + ACA_tab[ind]=mtv.ACA(X=self.data[self.distribution==i],dim=dim, + sample_size=sample_size, + notion=notion, + solver=solver, + NRandom=NRandom, + n_refinements=n_refinements, + sphcap_shrink=sphcap_shrink, + alpha_Dirichlet=alpha_Dirichlet, + cooling_factor=cooling_factor, + cap_size=cap_size, + start=start, + space=space, + line_solver=line_solver, + bound_gc=bound_gc) + ACA_tab = ACA_tab[0] if self.distRef.shape[0]==1 else ACA_tab return ACA_tab def IsInConvexes(self,z): @@ -1328,7 +1365,7 @@ def IsInConvexes(self,z): Results - -------- + ------- belongingness : array_like The return respresents directions that best represents anomalies in the dataset. @@ -1343,32 +1380,31 @@ def _calcDet(self,mat:np.ndarray): Computes the determinant of a matrix Parametres - ----------- + ---------- mat: {array-like} Matrix to compute the determinant Results - ----------- + ------- Det: float determinant of the matrix """ # self._check_variables return mtv.calcDet(mat) - def computeMCD(self,mat:np.ndarray=None, h:float=1, mfull: int = 10, nstep: int = 7, hiRegimeCompleteLastComp: bool = True)->None: + def computeMCD(self, h:float=1., mfull: int = 10, nstep: int = 7, hiRegimeCompleteLastComp: bool = True)->None: + """ Compute Minimum Covariance Determinant (MCD) Parametres - ----------- + ---------- mat: {array-like} or None, default=None Matrix to compute MCD. If set to None, compute the MCD of the loaded dataset - h: int or float, default=1 + h: int or float, default=1. Represents the amount of data of the dataset used to compute the MCD. - If the value is in the interval [0,1], it is treated as the percentage of dataset, - if the value is in the interval [n/2,n], it is treated as the amount of sample points. - It in the interval ]1,n/2[, the amount is rounded to n/2. + The value is in the interval [0,1], and it is treated as the percentage of dataset mfull: int, default=10 @@ -1378,17 +1414,19 @@ def computeMCD(self,mat:np.ndarray=None, h:float=1, mfull: int = 10, nstep: int hiRegimeCompleteLastComp: bool, default=True Results - ----------- + ------- Minimum Covariance Determinant (MCD): {array-like} """ - self._check_variables(h) # check if h is in the acceptable range - if h>0 and h<=1: # transform h in the good value for MCD function - h=int(h*self._nSamples) - elif hNone: """ @@ -1423,7 +1461,7 @@ def change_dataset(self,newDataset:np.ndarray,newY:np.ndarray=None, newDistribut except:pass try: self.distribution=np.concatenate((self.distribution,newDistribution)) # try for distribution - self.distRef=np.unique(self.distribution) + # self.distRef=np.unique(self.distribution) except: self.distribution=np.concatenate((self.distribution,np.repeat(0,newDataset.shape[0]))) # try for distribution else: @@ -1432,10 +1470,11 @@ def change_dataset(self,newDataset:np.ndarray,newY:np.ndarray=None, newDistribut except:pass try: self.distribution=newDistribution # try for distribution - self.distRef=np.unique(self.distribution) + # self.distRef=np.unique(self.distribution) except: self.distribution=np.repeat(0,newDataset.shape[0]) - self.distRef=np.unique(self.distribution) + # self.distRef=np.unique(self.distribution) + self.distRef,self._nSamples=np.unique(self.distribution,return_counts=True) return self #### auxiliar functions #### def set_seed(self,seed:int=2801)->None: @@ -1516,9 +1555,9 @@ def _check_variables(self,**kwargs)->None: if value not in self.approxOption: raise ValueError(f"Only output_option possibilities are {self.approxOption}, got {value}.") if key=="h": - assert type(value)==int or type(value)==float, f"h must be a float or int, got {type(value)}" - if value<=0 or value>self._nSamples: - raise ValueError(f"h must be in the range from 0 to {self._nSamples}, got {value}.") + assert type(value)==float, f"h must be a float, got {type(value)}" + # if value<=0 or value>np.min(self._nSamples): + # raise ValueError(f"h must be in the range from 0 to {self._nSamples}, got {value}.") def _check_CUDA(self,CUDA,solver): if solver not in ["simplerandom", "refinedrandom"] and CUDA==True: diff --git a/depth/model/multivariate/ACA_wrapper.py b/depth/model/multivariate/ACA_wrapper.py index 398e76a..af8bed3 100644 --- a/depth/model/multivariate/ACA_wrapper.py +++ b/depth/model/multivariate/ACA_wrapper.py @@ -9,22 +9,18 @@ import sys, os, glob from .import_CDLL import libACA -def ACA(X, dim = 1, sample_size = None, sample = None, notion = "projection", # Can't use halfspace with NelderMead +def ACA(X, dim = 1, sample_size = None, notion = "projection", # Can't use halfspace with NelderMead solver = "neldermead", NRandom = 100, n_refinements = 10, sphcap_shrink = 0.5, alpha_Dirichlet = 1.25, cooling_factor = 0.95, cap_size = 1, start = "mean", space = "sphere", line_solver = "goldensection", bound_gc = True): z=X.copy() - if(sample_size != None and sample == None): # Run method on a (specified) sample + if(sample_size != None): # Run method on a (specified) sample ind = np.random.default_rng().choice(X.shape[0], size=sample_size, replace=False) X = X[ind] - elif(sample_size == None and sample is not None): - ind = sample - X = X[sample] - elif(sample_size != None and sample is not None): - print("Can't give size of uniform sampling and your own index for sampling") - return(None) + else: + pass # Check arguments depth_indice = check_depth(notion) @@ -37,8 +33,7 @@ def ACA(X, dim = 1, sample_size = None, sample = None, notion = "projection", # try: n, d = X.shape except ValueError: - n = X.shape[0] - d = 1 + n,d = X.shape[0],1 basis = np.eye(d, dtype=np.double) d_aca = d iter_aca = dim diff --git a/depth/model/multivariate/BetaSkeleton.py b/depth/model/multivariate/BetaSkeleton.py index eb17e4b..c682d9f 100644 --- a/depth/model/multivariate/BetaSkeleton.py +++ b/depth/model/multivariate/BetaSkeleton.py @@ -2,17 +2,18 @@ from ctypes import * from multiprocessing import * import math -import sklearn.covariance as sk +# import sklearn.covariance as sk import sys, os, glob import platform from .import_CDLL import libExact -def MCD_fun(data,alpha,NeedLoc=False): - cov = sk.MinCovDet(support_fraction=alpha).fit(data) - if NeedLoc:return([cov.covariance_,cov.location_]) - else:return(cov.covariance_) +# def MCD_fun(data,alpha,NeedLoc=False): +# cov = sk.MinCovDet(support_fraction=alpha).fit(data) +# if NeedLoc:return([cov.covariance_,cov.location_]) +# else:return(cov.covariance_) -def betaSkeleton(x, data, beta = 2, distance = "Lp", Lp_p = 2, mah_estimate = "moment", mah_parMcd = 0.75): +def betaSkeleton(x, data, beta = 2, distance = "Lp", Lp_p = 2, mah_estimate = "moment", mah_parMcd = 0.75, + covMCD=None): points_list=data.flatten() objects_list=x.flatten() if (distance == "Mahalanobis"): @@ -22,8 +23,8 @@ def betaSkeleton(x, data, beta = 2, distance = "Lp", Lp_p = 2, mah_estimate = "m else: if(mah_estimate == "moment"): tmpCov = np.cov(np.transpose(data)) - elif (mah_estimate == "MCD"): - tmpCov = MCD_fun(data, mah_parMcd) + elif (mah_estimate.lower() == "mcd"): + tmpCov = covMCD else: print("Wrong argument \"mah_estimate\", should be one of \"moment\", \"MCD\", \"none\"") diff --git a/depth/model/multivariate/Depth_approximation.py b/depth/model/multivariate/Depth_approximation.py index 41f2c0a..02e6f11 100644 --- a/depth/model/multivariate/Depth_approximation.py +++ b/depth/model/multivariate/Depth_approximation.py @@ -111,7 +111,7 @@ def depth_approximation(z, for i in range(n_z): if(np.sum(z[i]*best_directions[i]) < np.sum(z[i]*(-best_directions[i]))): best_directions[i] = -best_directions[i] - # print(option) + depths=np.nan_to_num(depths, nan=0) if(option == 1): return depths diff --git a/depth/model/multivariate/L2.py b/depth/model/multivariate/L2.py index d446a3d..293d941 100644 --- a/depth/model/multivariate/L2.py +++ b/depth/model/multivariate/L2.py @@ -1,16 +1,16 @@ import numpy as np from ctypes import * from multiprocessing import * -import sklearn.covariance as sk +# import sklearn.covariance as sk import sys, os, glob import platform -def MCD_fun(data,alpha,NeedLoc=False): - cov = sk.MinCovDet(support_fraction=alpha).fit(data) - if NeedLoc:return([cov.covariance_,cov.location_]) - else:return(cov.covariance_) +# def MCD_fun(data,alpha,NeedLoc=False): +# cov = sk.MinCovDet(support_fraction=alpha).fit(data) +# if NeedLoc:return([cov.covariance_,cov.location_]) +# else:return(cov.covariance_) -def L2(x, data,mah_estimate='moment',mah_parMcd=0.75): +def L2(x, data,mah_estimate='moment',mah_parMcd=0.75,covMCD=None): points_list=data.flatten() objects_list=x.flatten() @@ -19,8 +19,8 @@ def L2(x, data,mah_estimate='moment',mah_parMcd=0.75): else: if mah_estimate=='moment': cov=np.cov(np.transpose(data)) - elif mah_estimate=='MCD': - cov=MCD_fun(data, mah_parMcd) + elif mah_estimate.lower()=='mcd': + cov=covMCD else : print("Wrong argument \"mah.estimate\", should be one of \"moment\", \"MCD\", \"none\"") print("moment is used") diff --git a/depth/model/multivariate/MCD.py b/depth/model/multivariate/MCD.py index 47d6cfc..2def537 100644 --- a/depth/model/multivariate/MCD.py +++ b/depth/model/multivariate/MCD.py @@ -3,34 +3,28 @@ from scipy.stats import chi2 from .import_CDLL import libExact -def MCD(data, h, seed=None, mfull = 10, nstep = 7, hiRegimeCompleteLastComp = True): +def MCD(data, h, seed=1, mfull = 10, nstep = 7, hiRegimeCompleteLastComp = True): try: n, d = data.shape except ValueError: - n = data.shape[0] - d = 1 + n, d = data.shape[0], 1 hParam = pointer(c_int(h)) numPoints = pointer(c_int(n)) dimension = pointer(c_int(d)) - points_list=data.flatten() + points_list=np.ascontiguousarray(data, dtype=np.float64).flatten() points=(c_double*len(points_list))(*points_list) - points=pointer(points) + c_points=cast(points, POINTER(c_double)) + # points=pointer(points) - if seed==None: - seeded = False - seed=0 - else: - seeded = True - - seed=pointer((c_int(seed))) - c_seeded = c_bool(seeded) + c_seed=(c_int(seed)) cov_size = d*d # print("cov_size",cov_size) - mat_MCD=pointer((c_double*(cov_size))(*np.zeros((cov_size)))) + mat_MCD=(c_double*cov_size)(*([0.0]*cov_size)) + c_mat_MCD=cast(mat_MCD, POINTER(c_double)) chisqr05 = chi2(d).isf(0.5) chisqr0975 = chi2(d).isf(0.025) @@ -42,15 +36,94 @@ def MCD(data, h, seed=None, mfull = 10, nstep = 7, hiRegimeCompleteLastComp = Tr c_nstep = c_int(nstep) c_hiRegimeCompleteLastComp = c_bool(hiRegimeCompleteLastComp) - libExact.MinimumCovarianceDeterminantEstim(points, numPoints, dimension, hParam, seed, mat_MCD,c_chisqr05,c_chisqr0975,c_mfull,c_nstep,c_hiRegimeCompleteLastComp,c_seeded) + libExact.MinimumCovarianceDeterminantEstim.restype=None + libExact.MinimumCovarianceDeterminantEstim.argtypes = [ + POINTER(c_double), # points (flattened row-major matrix) + POINTER(c_int), # numPoints + POINTER(c_int), # dimension + POINTER(c_int), # hParam + POINTER(c_int), # seed + POINTER(c_double), # mat_MCD (output, d*d) + c_double, # chisqr05 + c_double, # chisqr0975 + c_int, # mfull + c_int, # nstep + c_bool, # hiRegimeCompleteLastComp + ] + libExact.MinimumCovarianceDeterminantEstim( + c_points, # POINTER(c_double) + byref(c_int(n)), # POINTER(c_int) + byref(c_int(d)), # POINTER(c_int) + byref(c_int(h)), # POINTER(c_int) + byref(c_seed), # POINTER(c_int) + c_mat_MCD, # POINTER(c_double) + c_chisqr05, # plain c_double + c_chisqr0975,# plain c_double + c_mfull, # plain c_int + c_nstep, # plain c_int + c_hiRegimeCompleteLastComp, + ) + res = np.zeros((d,d)) for i in range(d): for j in range(d): - res[i,j]=mat_MCD[0][i*d+j] + res[i,j]=c_mat_MCD[i*d+j] return res +# def MCD(data, h, seed=2801, mfull = 10, nstep = 7, hiRegimeCompleteLastComp = True): + +# try: +# n, d = data.shape +# except ValueError: +# n,d = data.shape[0],1 + +# c_hParam = c_int(h) +# c_numPoints = c_int(n) +# c_dimension = c_int(d) + +# points=(c_double*data.size)(*data.flatten().astype(np.float64)) +# # points=pointer(points) + +# c_seed=c_int(int(seed)) + +# cov_size = d*d +# c_mat_MCD=(c_double*(cov_size))(*([0]*cov_size)) +# chisqr05 = chi2(d).isf(0.5) +# chisqr0975 = chi2(d).isf(0.025) +# c_chisqr05 = c_double(chisqr05) +# c_chisqr0975 = c_double(chisqr0975) +# c_mfull = c_int(mfull) +# c_nstep = c_int(nstep) +# c_hiRegimeCompleteLastComp = c_bool(hiRegimeCompleteLastComp) + +# #MinimumCovarianceDeterminantEstim(double *points, int *numPoints, int *dimension, int *hParam, int *seed, double *mat_MCD, double chisqr05, double chisqr0975, int mfull, +# # int nstep, bool hiRegimeCompleteLastComp) + +# libExact.MinimumCovarianceDeterminantEstim( +# points, +# byref(c_numPoints), +# byref(c_dimension), +# byref(c_hParam), +# byref(c_seed), +# c_mat_MCD, +# byref(c_chisqr05), +# byref(c_chisqr0975), +# byref(c_mfull), +# byref(c_nstep), +# byref(c_hiRegimeCompleteLastComp), +# ) + +# res = np.zeros((d,d)) +# print(c_mat_MCD[8]) +# # for i in range(d): +# # for j in range(d): +# # print(c_mat_MCD[i]) +# # res[i,j]=c_mat_MCD[0][i*d+j] + +# return res + MCD.__doc__= """ Description diff --git a/depth/model/multivariate/Mahalanobis.py b/depth/model/multivariate/Mahalanobis.py index dd83192..98cd118 100644 --- a/depth/model/multivariate/Mahalanobis.py +++ b/depth/model/multivariate/Mahalanobis.py @@ -1,15 +1,11 @@ import numpy as np from ctypes import * -import sklearn.covariance as sk +# import sklearn.covariance as sk from .Depth_approximation import depth_approximation import sys, os, glob import platform from .import_CDLL import libExact,libApprox -def MCD_fun(data,alpha,NeedLoc=False): - cov = sk.MinCovDet(support_fraction=alpha).fit(data) - if NeedLoc:return([cov.covariance_,cov.location_]) - else:return(cov.covariance_) def mahalanobis(x, data, exact=True, mah_estimate="moment", mah_parMcd = 0.75, solver = "neldermead", @@ -23,7 +19,8 @@ def mahalanobis(x, data, exact=True, mah_estimate="moment", mah_parMcd = 0.75, start = "mean", space = "sphere", line_solver = "goldensection", - bound_gc = True): + bound_gc = True, + covMCD=None): if exact: points_list=data.flatten() @@ -39,8 +36,8 @@ def mahalanobis(x, data, exact=True, mah_estimate="moment", mah_parMcd = 0.75, dimension=pointer(c_int(len(data[0]))) if mah_estimate=='moment': # compute cov based on user choice PY_MatMCD=np.cov(np.transpose(data)) - else: # compute cov based on user choice - PY_MatMCD=MCD_fun(data,mah_parMcd) + elif mah_estimate=="mcd": # compute cov based on user choice + PY_MatMCD=covMCD PY_MatMCD=PY_MatMCD.flatten(order='C') mat_MCD=pointer((c_double*len(PY_MatMCD))(*PY_MatMCD)) diff --git a/depth/model/multivariate/SimplicialVolume.py b/depth/model/multivariate/SimplicialVolume.py index 445a839..ae70144 100644 --- a/depth/model/multivariate/SimplicialVolume.py +++ b/depth/model/multivariate/SimplicialVolume.py @@ -1,7 +1,7 @@ import numpy as np from ctypes import * from multiprocessing import * -import sklearn.covariance as sk +# import sklearn.covariance as sk import scipy.special as scspecial import sys, os, glob import platform @@ -13,10 +13,10 @@ def longtoint(k): k2 = int(k - k1*limit) return np.array([k1,k2]) -def MCD_fun(data,alpha,NeedLoc=False): - cov = sk.MinCovDet(support_fraction=alpha).fit(data) - if NeedLoc:return([cov.covariance_,cov.location_]) - else:return(cov.covariance_) +# def MCD_fun(data,alpha,NeedLoc=False): +# cov = sk.MinCovDet(support_fraction=alpha).fit(data) +# if NeedLoc:return([cov.covariance_,cov.location_]) +# else:return(cov.covariance_) def calcDet(A): dim_pointer = pointer(c_int(A.shape[0])) @@ -27,7 +27,8 @@ def calcDet(A): return res_pointer[0] def simplicialVolume(x, data, exact = True, k = 0.05, - mah_estimate = "moment", mah_parMCD = 0.75, seed = 0): + mah_estimate = "moment", mah_parMCD = 0.75, seed = 0, + covMCD=None): points_list=data.flatten() objects_list=x.flatten() if (mah_estimate == "none"): @@ -38,7 +39,7 @@ def simplicialVolume(x, data, exact = True, k = 0.05, covEst=np.cov(np.transpose(data)) elif (mah_estimate == "MCD") : useCov = 2 - covEst = MCD_fun(data, mah_parMCD) + covEst = covMCD else: print("Wrong argument \"mah.estimate\", should be one of \"moment\", \"MCD\", \"none\"") print("moment is use") diff --git a/depth/model/multivariate/Spatial.py b/depth/model/multivariate/Spatial.py index 07ffc3d..f0386d6 100644 --- a/depth/model/multivariate/Spatial.py +++ b/depth/model/multivariate/Spatial.py @@ -5,12 +5,12 @@ import platform import sklearn.covariance as sk -def MCD_fun(data,alpha,NeedLoc=False): - cov = sk.MinCovDet(support_fraction=alpha).fit(data) - if NeedLoc:return([cov.covariance_,cov.location_]) - else:return(cov.covariance_) +# def MCD_fun(data,alpha,NeedLoc=False): +# cov = sk.MinCovDet(support_fraction=alpha).fit(data) +# if NeedLoc:return([cov.covariance_,cov.location_]) +# else:return(cov.covariance_) -def spatial(x, data,mah_estimate='moment',mah_parMcd=0.75): +def spatial(x, data,mah_estimate='moment',mah_parMcd=0.75,covMCD=None): depths_tab=[] if mah_estimate=='none': @@ -19,8 +19,8 @@ def spatial(x, data,mah_estimate='moment',mah_parMcd=0.75): cov[:]=np.nan elif mah_estimate=='moment': cov=np.cov(np.transpose(data)) - elif mah_estimate=='MCD': - cov=MCD_fun(data,mah_parMcd) + elif mah_estimate.lower()=='mcd': + cov=covMCD if np.sum(np.isnan(cov))==0: w,v=np.linalg.eig(cov) lambda1=np.linalg.inv(np.matmul(v,np.diag(np.sqrt(w)))) diff --git a/depth/src/MCD.cpp b/depth/src/MCD.cpp index c25b984..e34ac5c 100644 --- a/depth/src/MCD.cpp +++ b/depth/src/MCD.cpp @@ -9,21 +9,18 @@ vector split(int n){//split n into (at most 5) bins and return the vector v val = 300; } else{ - k = (int) n/300; - if (n % 300 == 0) k = k-1; + k = std::max(1, (int)(n / 300)); + if (n % 300 == 0 && k > 1) k = k-1; val = (int) n/k; } vector v(k,val);//will contain number of elements in each bin // Adding extra numbers to add up exactly to n + remainder = n - k*val; if (n<1500){ - remainder = n - k*val; - for (int i=0; i < v.size(); i++){ - if (remainder>0){ - v[i] += 1; - remainder -= 1; - } + for (int i=0; i < (int)v.size() && remainder > 0; i++, remainder--){ + v[i] ++; } } return v; @@ -73,6 +70,7 @@ double det(TDMatrix M, int d){ } if (amax < eps_pivot) { delete[] colp; + deleteM(A); return 0; } // Column swap @@ -313,14 +311,16 @@ void ExactUnivariateMcd(TDMatrix X, int n, int h, double* T, TDMatrix M){ for(int i=1; i index(n, 0); std::random_device rd; std::mt19937 g(rd()); - if (seeded){ - g.seed(*seed); - } - int bestIndex; + g.seed(*seed); + int bestIndex=0; double finalDet = DBL_MAX; double tempDet; if (h==n){ @@ -391,9 +389,9 @@ void Mcd(TDMatrix X, int n, int d, int h, double* mat_MCD, double chisqr05, doub best_indices[i] = best(10, all_det[i] , rep); } - // Merging + // // Merging - // First count the occurences + // // First count the occurences int counters[1500] = {0}; for (int i = 0;i 0 ) nMerged +=1; } - - // Now build Xmerged + + // // Now build Xmerged TDMatrix Xmerged = newM(nMerged, d); - // Reinitialise counters, to then check first occurence or not + // // Reinitialise counters, to then check first occurence or not for (int i=0;i<1500;i++) counters[i] = 0; - // Fill in Xmerged + // // Fill in Xmerged int movingInd = 0;// index for Xmerged filling step by step int tempoInd;// current index when checking through the best solutions to avoid recomputation of the index at each call for (int i = 0;i mergedBest = best(mfull, merged_all_det , nSelect); - // Full dataset computation + // // Full dataset computation - // for simplicity recopy the previous best results (T,S) in a fresh compilation + // // for simplicity recopy the previous best results (T,S) in a fresh compilation vector full_all_cov(mfull); vector full_all_T(mfull); int convertedM, convertedI; // i=2,n<=600 - // initialise index vector + // // initialise index vector for (int i = 0 ; i != index.size() ; i++) { // assume index.size() is n ! index[i] = i; } - // big loop 500 + // // big loop 500 TDMatrix* all_cov = new TDMatrix[500]; vector* all_index = new vector[500]; double* all_det = new double[500]; @@ -529,8 +534,9 @@ void Mcd(TDMatrix X, int n, int d, int h, double* mat_MCD, double chisqr05, doub // update Cov & T according to most recent update of index MeanCovUp(all_index[i], T, all_cov[i], X, Xh, n, d, h); all_det[i] = det(all_cov[i],d); - } - // take the 10 best results + } + + // // take the 10 best results vector index500(500, 0); //index to help find 10 best iota(index500.begin(), index500.end(), 0); std::nth_element(index500.begin(),index500.begin()+9,index500.end(), @@ -539,6 +545,7 @@ void Mcd(TDMatrix X, int n, int d, int h, double* mat_MCD, double chisqr05, doub } ); //ten first elements now are index of ten smallest det int ind; + // run until convergence for each of the 10 best for(int i=0;i<10;i++){ ind = index500[i]; @@ -556,43 +563,43 @@ void Mcd(TDMatrix X, int n, int d, int h, double* mat_MCD, double chisqr05, doub } - cout << " M cov " << endl; - for (int k=0; k < d; k++){ - for (int p=0; p < d; p++){ - std::cout << M[k][p] << " " ; - } - std::cout << std::endl ; - } + // cout << " M cov " << endl; + // for (int k=0; k < d; k++){ + // for (int p=0; p < d; p++){ + // std::cout << M[k][p] << " " ; + // } + // std::cout << std::endl ; + // } - // Ultimate reweighting - // std::cout << "det M before reweighting " << det(M,d) << std::endl; + // // Ultimate reweighting + // // std::cout << "det M before reweighting " << det(M,d) << std::endl; DistanceUp(X, n, d, distTab, T, M); double medi = DataDepth::med(distTab,n); double medi2 = medi*medi; - // if (medi2==0){ - // for(int i=0;i split(int n); vector best(int p, vector& all_det ,int rep); -void biased_cov(TDMatrix X, int n, int d, TDMatrix S); -void IndexUp(vector& index, double* distTab); -void MeanCovUp(vector& index, double* T,TDMatrix S, TDMatrix X, TDMatrix Xh,int n, int d, int h); -void DistanceUp(TDMatrix X, int n, int d, double* distTab, double* T,TDMatrix S); -void cstep(vector& index, double* distTab, double* T,TDMatrix S, TDMatrix X, TDMatrix Xh,int n, int d, int h); -void cstep_TSstart(vector& index, double* distTab, double* T,TDMatrix S, TDMatrix X, TDMatrix Xh,int n, int d, int h); -void mcd_routine(vector& index, double* distTab, double* T,TDMatrix S, TDMatrix X, TDMatrix Xh,int n, int d, int h); -void ExactUnivariateMcd(TDMatrix X, int n, int h, double* T, TDMatrix M); -void Mcd(TDMatrix X, int n, int d, int h, double* mat_MCD, double chisqr05, double chisqr0975, int mfull, int nstep, bool hiRegimeCompleteLastComp,int *seed, bool seeded); +// void biased_cov(TDMatrix X, int n, int d, TDMatrix S); +// void IndexUp(vector& index, double* distTab); +// void MeanCovUp(vector& index, double* T,TDMatrix S, TDMatrix X, TDMatrix Xh,int n, int d, int h); +// void DistanceUp(TDMatrix X, int n, int d, double* distTab, double* T,TDMatrix S); +// void cstep(vector& index, double* distTab, double* T,TDMatrix S, TDMatrix X, TDMatrix Xh,int n, int d, int h); +// void cstep_TSstart(vector& index, double* distTab, double* T,TDMatrix S, TDMatrix X, TDMatrix Xh,int n, int d, int h); +// void mcd_routine(vector& index, double* distTab, double* T,TDMatrix S, TDMatrix X, TDMatrix Xh,int n, int d, int h); +// void ExactUnivariateMcd(TDMatrix X, int n, int h, double* T, TDMatrix M); +void Mcd(TDMatrix X, int n, int d, int h, double* mat_MCD, double chisqr05, + double chisqr0975, int mfull, int nstep, bool hiRegimeCompleteLastComp,int *seed); +// void Mcd(TDMatrix X, int* n, int* d, int* h, double* mat_MCD, double* chisqr05, double* chisqr0975, int* mfull, int* nstep, bool* hiRegimeCompleteLastComp,int *seed) \ No newline at end of file diff --git a/depth/src/ddalpha.cpp b/depth/src/ddalpha.cpp index 878c723..9f0f5ca 100644 --- a/depth/src/ddalpha.cpp +++ b/depth/src/ddalpha.cpp @@ -2,7 +2,7 @@ File: ddalpha.cpp Created by: Pavlo Mozharovskyi, Oleksii Pokotylo, Arturo Castellanos First published: 28.02.2013 - Last revised: 24.10.2024 + Last revised: 28.05.2026 Defines the exported functions for the former 'ddalpha'-package, now 'data-depth'-library. @@ -70,6 +70,7 @@ void IsInConvexes(double *points, int *dimension, int *cardinalities, } TIntMatrix answers(o.size()); int error = 0; + // int InConvexes(TMatrix &points, TVariables &cardinalities, TMatrix &objects, int &Error, TIntMatrix *areInConvexes) InConvexes(x, cars, o, error, &answers); for (int i = 0; i < numObjects[0]; i++) for (int j = 0; j < numClasses[0]; j++){ @@ -77,40 +78,6 @@ void IsInConvexes(double *points, int *dimension, int *cardinalities, } } -// void IsInConvexes(double *points, int *dimension, int *cardinalities, -// int *numClasses, double *objects, int *numObjects, int *seed, int *isInConvexes, -// int *cumSum, int *distrSeq){ -// setSeed(*seed); -// for(int i = 0; i < 2; i++){ -// cumSum[i]; -// } -// int numPoints = 0;for (int i = 0; i < numClasses[0]; i++){numPoints += cardinalities[i];} -// TMatrix x(numPoints); -// for (int i = 0; i < numPoints; i++){x[i] = TPoint(dimension[0]);} -// for (int i = 0; i < numPoints; i++){ -// for (int j = 0; j < dimension[0]; j++){ -// x[i][j] = points[i * dimension[0] + j]; -// } -// } -// TMatrix o(numObjects[0]); -// for (int i = 0; i < numObjects[0]; i++){o[i] = TPoint(dimension[0]);} -// for (int i = 0; i < numObjects[0]; i++){ -// for (int j = 0; j < dimension[0]; j++){ -// o[i][j] = objects[i * dimension[0] + j]; -// } -// } -// TVariables cars(numClasses[0]); -// for (int i = 0; i < numClasses[0]; i++){ -// cars[i] = cardinalities[i]; -// } -// TIntMatrix answers(o.size()); -// int error = 0; -// InConvexes(x, cars, o, error, &answers); -// for (int i = 0; i < numObjects[0]; i++) -// for (int j = 0; j < numClasses[0]; j++){ -// isInConvexes[numClasses[0]*i+j] = answers[i][j]; -// } -// } void ZDepth(double *points, double *objects, int *numPoints, int *numObjects, int *dimension, int *seed, double *depths){ @@ -306,12 +273,17 @@ void BetaSkeletonDepth(double *points, double *objects, int *numPoints, int *num delete[] s; } -void MinimumCovarianceDeterminantEstim(double *points, int *numPoints, int *dimension, int *hParam, int *seed, double *mat_MCD, double chisqr05, double chisqr0975, int mfull, int nstep, bool hiRegimeCompleteLastComp, bool seeded){ +void MinimumCovarianceDeterminantEstim(double *points, int *numPoints, + int *dimension, int *hParam, int *seed, double *mat_MCD, + double chisqr05, double chisqr0975, int mfull, int nstep, bool hiRegimeCompleteLastComp){ TDMatrix X = asMatrix(points, *numPoints, *dimension); - Mcd(X, *numPoints,*dimension, *hParam, mat_MCD, chisqr05, chisqr0975, mfull, nstep, hiRegimeCompleteLastComp, seed, seeded); - delete[] X; + + Mcd(X, *numPoints,*dimension, *hParam, mat_MCD, chisqr05, chisqr0975, mfull, nstep, + hiRegimeCompleteLastComp, seed); + // delete[] X; } + int main() { std::cout << "Hello Ddalpha!"; return 0; From 5405122f158e6a698a3f883b4ae3a4faaee84d3f Mon Sep 17 00:00:00 2001 From: Leonardo Leone Date: Fri, 12 Jun 2026 16:57:48 +0200 Subject: [PATCH 4/8] fix notion typo --- depth/model/DepthFunc.py | 52 ++++++++++++++++++++++++++++++++++++---- 1 file changed, 47 insertions(+), 5 deletions(-) diff --git a/depth/model/DepthFunc.py b/depth/model/DepthFunc.py index 1985bf6..a0c9fcb 100644 --- a/depth/model/DepthFunc.py +++ b/depth/model/DepthFunc.py @@ -493,10 +493,6 @@ def projection_based_func_depth(self, query, Number of random projections or optimization restarts used in computing projection-based depth. - notion - {"mahalanobis", "halfspace", "zonoid", "projection", "aprojection", "cexpchullstar", "cexpchull", "geometrical"}, - Which depth will be computed. - n_refinements For ``solver`` = ``refinedrandom`` or ``refinedgrid``, set the maximum of iteration for computing the depth of one point. @@ -592,8 +588,54 @@ def projection_based_func_depth(self, query, if option==1:return depth_array - else:return depth_array,direction_array + else:return depth_array,direction_array + + def general_func_depth(self,query,notion='halfspace', **kwargs): + """ + Compute non projection-based functional depth for query functional data with respect to a reference dataset. + + This function computes depth values of functional observations (in `query`) relative to a + reference dataset (`df`) using projection-based methods such as halfspace depth. + Each function (trajectory) is represented by a sequence of multivariate values over time. + + Parameters + ---------- + query : pandas.DataFrame + Query dataset containing functional observations whose depth will be computed + relative to `df`. Must have the same column structure as `df`. + + notion + {"mahalanobis", "halfspace", "zonoid", "projection", "aprojection", "cexpchullstar", "cexpchull", "geometrical"}, + Which depth will be computed. + + + Returns + ------- + depth_array : np.ndarray of shape (n_query,) + Array of depth values, where `n_query` is the number of functional observations + (unique `case_id`s) in the `query` dataset. + The return is the lowest comuted depth regarding all explored directions in space. + Notes + ----- + - If `timestamp` is of type `datetime64`, it is converted internally to seconds + relative to the global minimum timestamp (`t_min`). + - Duplicate timestamps within each `case_id` group are automatically dropped. + - Interpolation uses linear extrapolation outside the observed time range. + """ + self._check_depth(notion) + if type(query)==np.ndarray: + query=self._3Dnp_tp_pd(query,self.TSnp, self.CInp) + if query[self.timestamp_col].max()>self.t_max: + print(f"Values with {self.timestamp_col} greater the base set domain are excluded") + query.drop(query[query[self.timestamp_col]>self.t_max].index, inplace=True) + if query[self.timestamp_col].min()>self.t_min: + print(f"Values with {self.timestamp_col} smaller the base set domain are excluded") + query.drop(query[query[self.timestamp_col] Date: Fri, 3 Jul 2026 09:37:34 +0200 Subject: [PATCH 5/8] modify package dir search --- .gitignore | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.gitignore b/.gitignore index 3814d42..96ba493 100644 --- a/.gitignore +++ b/.gitignore @@ -4,7 +4,7 @@ __pycache__/ # C extensions *.so -*.dll + auto_examples From f0a0e091b378d51cdcb07c39c5228ad35eab17d2 Mon Sep 17 00:00:00 2001 From: Leonardo Leone Date: Fri, 3 Jul 2026 09:38:04 +0200 Subject: [PATCH 6/8] modify setup --- .gitignore | 2 +- setup.py | 7 +++---- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/.gitignore b/.gitignore index 96ba493..3814d42 100644 --- a/.gitignore +++ b/.gitignore @@ -4,7 +4,7 @@ __pycache__/ # C extensions *.so - +*.dll auto_examples diff --git a/setup.py b/setup.py index 3e9e4f2..26434bf 100644 --- a/setup.py +++ b/setup.py @@ -61,10 +61,9 @@ def build_extensions(self): # description="The package provides many procedures for calculating the depth of points in an empirical distribution for many notions of data depth", # long_description="The package provides many procedures for calculating the depth of points in an empirical distribution for many notions of data depth", # long_description_content_type="text/markdown", - packages=find_packages(), - # install_requires=['numpy','scipy','scikit-learn','matplotlib', - # "torch", - # "torchvision",], + #packages=find_packages(), + package_dir = {"": "depth"}, + packages=find_namespace_packages(where='depth'), include_package_data=True, data_files=[('depth/src', glob.glob("depth/docs/*"))], zip_safe=False From 7ca28a29b9db1d226b03ef98cb7faef95687e970 Mon Sep 17 00:00:00 2001 From: Leonardo Leone Date: Mon, 3 Aug 2026 16:44:30 +0200 Subject: [PATCH 7/8] Fix direction word typo --- depth/model/DepthEucl.py | 42 ++++++++++++++++++++-------------------- depth/model/DepthFunc.py | 2 +- setup.py | 8 +++++--- 3 files changed, 27 insertions(+), 25 deletions(-) diff --git a/depth/model/DepthEucl.py b/depth/model/DepthEucl.py index 336dd0c..3552dbb 100644 --- a/depth/model/DepthEucl.py +++ b/depth/model/DepthEucl.py @@ -212,18 +212,18 @@ def mahalanobis(self, x: np.ndarray = None, exact: bool = True, mah_estimate: Li If ``output_option=="final_depth_dir"`` returns: Tuple of array_like - Lowest Mahalanobis Detph - - Lowest depth respective sirection + - Lowest depth respective direction If ``output_option=="all_depth"`` returns: array_like - Lowest Mahalanobis Detph - - Lowest depth respective sirection + - Lowest depth respective direction - All computed depths If ``output_option=="all_depth_directions"`` returns: array_like - Lowest Mahalanobis Detph - - Lowest depth respective sirection + - Lowest depth respective direction - All computed depths - All respective directions @@ -336,18 +336,18 @@ def aprojection(self,x:np.ndarray=None,solver: str = "neldermead", NRandom: int If ``output_option=="final_depth_dir"`` returns: Tuple of array_like - Lowest Asymmetrical Projection Detph - - Lowest depth respective sirection + - Lowest depth respective direction If ``output_option=="all_depth"`` returns: array_like - Lowest Asymmetrical Projection Detph - - Lowest depth respective sirection + - Lowest depth respective direction - All computed depths If ``output_option=="all_depth_directions"`` returns: array_like - Lowest Asymmetrical Projection Detph - - Lowest depth respective sirection + - Lowest depth respective direction - All computed depths - All respective directions @@ -561,18 +561,18 @@ def cexpchullstar(self,x: np.ndarray=None, solver: str = "neldermead", NRandom: If ``output_option=="final_depth_dir"`` returns: Tuple of array_like - Lowest Continuous Modified Explected Convex Hull Detph - - Lowest depth respective sirection + - Lowest depth respective direction If ``output_option=="all_depth"`` returns: array_like - Lowest Continuous Modified Explected Convex Hull Detph - - Lowest depth respective sirection + - Lowest depth respective direction - All computed depths If ``output_option=="all_depth_directions"`` returns: array_like - Lowest Continuous Modified Explected Convex Hull Detph - - Lowest depth respective sirection + - Lowest depth respective direction - All computed depths - All respective directions @@ -661,18 +661,18 @@ def geometrical(self,x:np.ndarray=None,solver: str = "neldermead", NRandom: int If ``output_option=="final_depth_dir"`` returns: Tuple of array_like - Lowest Geometrical Detph - - Lowest depth respective sirection + - Lowest depth respective direction If ``output_option=="all_depth"`` returns: array_like - Lowest Geometrical Detph - - Lowest depth respective sirection + - Lowest depth respective direction - All computed depths If ``output_option=="all_depth_directions"`` returns: array_like - Lowest Geometrical Detph - - Lowest depth respective sirection + - Lowest depth respective direction - All computed depths - All respective directions @@ -763,18 +763,18 @@ def halfspace(self, x:np.ndarray=None,exact: bool = True,method: str = "recursiv If ``output_option=="final_depth_dir"`` returns: Tuple of array_like - Lowest Halfspace (Tukey) Detph - - Lowest depth respective sirection + - Lowest depth respective direction If ``output_option=="all_depth"`` returns: array_like - Lowest Halfspace (Tukey) Detph - - Lowest depth respective sirection + - Lowest depth respective direction - All computed depths If ``output_option=="all_depth_directions"`` returns: array_like - Lowest Halfspace (Tukey) Detph - - Lowest depth respective sirection + - Lowest depth respective direction - All computed depths - All respective directions @@ -959,18 +959,18 @@ def projection(self,x:np.ndarray=None,solver: str = "neldermead",NRandom: int = If ``output_option=="final_depth_dir"`` returns: Tuple of array_like - Lowest Projection Detph - - Lowest depth respective sirection + - Lowest depth respective direction If ``output_option=="all_depth"`` returns: array_like - Lowest Projection Detph - - Lowest depth respective sirection + - Lowest depth respective direction - All computed depths If ``output_option=="all_depth_directions"`` returns: array_like - Lowest Projection Detph - - Lowest depth respective sirection + - Lowest depth respective direction - All computed depths - All respective directions """ @@ -1222,18 +1222,18 @@ def zonoid(self,x:np.ndarray=None, exact:bool=True, If ``output_option=="final_depth_dir"`` returns: Tuple of array_like - Lowest Zonoid Detph - - Lowest depth respective sirection + - Lowest depth respective direction If ``output_option=="all_depth"`` returns: array_like - Lowest Zonoid Detph - - Lowest depth respective sirection + - Lowest depth respective direction - All computed depths If ``output_option=="all_depth_directions"`` returns: array_like - Lowest Zonoid Detph - - Lowest depth respective sirection + - Lowest depth respective direction - All computed depths - All respective directions diff --git a/depth/model/DepthFunc.py b/depth/model/DepthFunc.py index a0c9fcb..04e2607 100644 --- a/depth/model/DepthFunc.py +++ b/depth/model/DepthFunc.py @@ -541,7 +541,7 @@ def projection_based_func_depth(self, query, If ``output_option=="final_depth_dir"`` returns: Tuple of array_like - Lowest Asymmetrical Projection Detph - - Lowest depth respective sirection + - Lowest depth respective direction Notes ----- diff --git a/setup.py b/setup.py index 26434bf..1814b93 100644 --- a/setup.py +++ b/setup.py @@ -61,9 +61,11 @@ def build_extensions(self): # description="The package provides many procedures for calculating the depth of points in an empirical distribution for many notions of data depth", # long_description="The package provides many procedures for calculating the depth of points in an empirical distribution for many notions of data depth", # long_description_content_type="text/markdown", - #packages=find_packages(), - package_dir = {"": "depth"}, - packages=find_namespace_packages(where='depth'), + packages=find_packages(), + + # install_requires=['numpy','scipy','scikit-learn','matplotlib', + # "torch", + # "torchvision",], include_package_data=True, data_files=[('depth/src', glob.glob("depth/docs/*"))], zip_safe=False From 5fda464a554e4aa9454d5ac398a8185ab36d7547 Mon Sep 17 00:00:00 2001 From: Leonardo Leone Date: Thu, 13 Aug 2026 12:08:17 +0200 Subject: [PATCH 8/8] Put seed in all approximative functions py + cpp --- depth/model/DepthEucl.py | 22 ++++++++++--------- depth/model/multivariate/Aprojection.py | 9 +++++--- .../model/multivariate/CUDA_approximation.py | 4 ++-- depth/model/multivariate/Cexpchull.py | 5 +++-- depth/model/multivariate/Cexpchullstar.py | 5 +++-- .../model/multivariate/Depth_approximation.py | 9 +++++--- depth/model/multivariate/Geometrical.py | 5 +++-- depth/model/multivariate/Halfspace.py | 7 +++--- depth/model/multivariate/Mahalanobis.py | 5 +++-- depth/model/multivariate/Projection.py | 7 +++--- depth/model/multivariate/Zonoid.py | 4 ++-- depth/src/ACA_wrapper.cpp | 15 +++++++++++-- depth/src/ProjectionDepths.cpp | 6 ++--- depth/src/ProjectionDepths.h | 2 +- depth/src/depth_wrapper.cpp | 17 +++++++++++--- 15 files changed, 79 insertions(+), 43 deletions(-) diff --git a/depth/model/DepthEucl.py b/depth/model/DepthEucl.py index 3552dbb..4eacd59 100644 --- a/depth/model/DepthEucl.py +++ b/depth/model/DepthEucl.py @@ -269,7 +269,7 @@ def mahalanobis(self, x: np.ndarray = None, exact: bool = True, mah_estimate: Li alpha_Dirichlet=alpha_Dirichlet, cooling_factor=cooling_factor, cap_size=cap_size, start=start, space=space, line_solver=line_solver, bound_gc=bound_gc,option=option, - covMCD=covMCD + covMCD=covMCD, seed=self.seed ) #compute depth value if evaluate_dataset==False: if exact or option==1:self.mahalanobisDepth[ind]=DM # assign value - exact or option 1 @@ -376,10 +376,10 @@ def aprojection(self,x:np.ndarray=None,solver: str = "neldermead", NRandom: int for ind,d in enumerate(self.distRef): if CUDA and self.CUDA:DAP=mtv.aprojection(x=x,data=self.dataCuda[:,self.distribution==d],solver=solver,NRandom=NRandom,option=option, n_refinements=n_refinements, sphcap_shrink=sphcap_shrink, alpha_Dirichlet=alpha_Dirichlet, cooling_factor=cooling_factor, - cap_size=cap_size,start=start,space=space,line_solver=line_solver,bound_gc=bound_gc,CUDA=CUDA,device=self.device) #compute depth value + cap_size=cap_size,start=start,space=space,line_solver=line_solver,bound_gc=bound_gc,CUDA=CUDA,device=self.device, seed=self.seed) #compute depth value else:DAP=mtv.aprojection(x=x,data=self.data[self.distribution==d],solver=solver,NRandom=NRandom,option=option, n_refinements=n_refinements, sphcap_shrink=sphcap_shrink, alpha_Dirichlet=alpha_Dirichlet, cooling_factor=cooling_factor, - cap_size=cap_size,start=start,space=space,line_solver=line_solver,bound_gc=bound_gc,CUDA=CUDA,device=self.device) #compute depth value + cap_size=cap_size,start=start,space=space,line_solver=line_solver,bound_gc=bound_gc,CUDA=CUDA,device=self.device, seed=self.seed) #compute depth value if evaluate_dataset==False: if option==1:self.aprojectionDepth[ind]=DAP # assign val option 1 elif option==2:self.aprojectionDepth[ind],self.aprojectionDir[ind]=DAP # assign value option 2 @@ -497,7 +497,7 @@ def cexpchull(self,x: np.ndarray=None,solver:str= "neldermead",NRandom:int = 100 x=x, data=self.data[self.distribution==d],solver=solver,NRandom=NRandom,option=option,n_refinements=n_refinements, sphcap_shrink=sphcap_shrink,alpha_Dirichlet =alpha_Dirichlet,cooling_factor=cooling_factor, cap_size =cap_size,start =start,space =space,line_solver =line_solver,bound_gc =bound_gc, - ) # compute depth + seed=self.seed) # compute depth if evaluate_dataset==False: if option==1:self.cexpchullDepth[ind]=DC # assign value elif option==2:self.cexpchullDepth[ind],self.cexpchullDir[ind]=DC # assign value @@ -600,7 +600,8 @@ def cexpchullstar(self,x: np.ndarray=None, solver: str = "neldermead", NRandom: for ind,d in enumerate(self.distRef): DC=mtv.cexpchullstar(x=x,data=self.data[self.distribution==d], solver=solver, NRandom=NRandom, option=option, n_refinements=n_refinements, sphcap_shrink=sphcap_shrink, alpha_Dirichlet=alpha_Dirichlet, cooling_factor=cooling_factor, - cap_size=cap_size,start=start, space=space, line_solver=line_solver, bound_gc=bound_gc) + cap_size=cap_size,start=start, space=space, line_solver=line_solver, bound_gc=bound_gc, + seed=self.seed) if evaluate_dataset==False: if option==1:self.cexpchullstarDepth[ind]=DC # assign value elif option==2:self.cexpchullstarDepth[ind],self.cexpchullstarDir[ind]=DC # assign value @@ -701,7 +702,8 @@ def geometrical(self,x:np.ndarray=None,solver: str = "neldermead", NRandom: int for ind,d in enumerate(self.distRef): DG=mtv.geometrical(x=x,data=self.data[self.distribution==d], solver=solver, NRandom=NRandom, option=option, n_refinements=n_refinements, sphcap_shrink=sphcap_shrink, alpha_Dirichlet=alpha_Dirichlet, cooling_factor=cooling_factor, - cap_size=cap_size,start=start, space=space, line_solver=line_solver, bound_gc=bound_gc) + cap_size=cap_size,start=start, space=space, line_solver=line_solver, bound_gc=bound_gc, + seed=self.seed) if evaluate_dataset==False: if option==1:self.geometricalDepth[ind]=DG # assign value elif option==2:self.geometricalDepth[ind],self.geometricalDir[ind]=DG # assign value @@ -806,12 +808,12 @@ def halfspace(self, x:np.ndarray=None,exact: bool = True,method: str = "recursiv if CUDA==True and self.CUDA==True:DH=mtv.halfspace(x=x,data=self.dataCuda[:,self.distribution==d],exact=exact,method=method, solver=solver,NRandom=NRandom,option=option,n_refinements=n_refinements,sphcap_shrink=sphcap_shrink, alpha_Dirichlet=alpha_Dirichlet,cooling_factor=cooling_factor,cap_size=cap_size,start=start, - space=space,line_solver=line_solver,bound_gc=bound_gc,CUDA=CUDA, device=self.device, + space=space,line_solver=line_solver,bound_gc=bound_gc,CUDA=CUDA, device=self.device,seed=self.seed ) elif CUDA==False:DH=mtv.halfspace(x=x,data=self.data[self.distribution==d],exact=exact,method=method, solver=solver,NRandom=NRandom,option=option,n_refinements=n_refinements,sphcap_shrink=sphcap_shrink, alpha_Dirichlet=alpha_Dirichlet,cooling_factor=cooling_factor,cap_size=cap_size,start=start, - space=space,line_solver=line_solver,bound_gc=bound_gc,CUDA=CUDA, + space=space,line_solver=line_solver,bound_gc=bound_gc,CUDA=CUDA,seed=self.seed ) if evaluate_dataset==False: if option==1 or exact==True:self.halfspaceDepth[ind]=DH # assign value @@ -1001,12 +1003,12 @@ def projection(self,x:np.ndarray=None,solver: str = "neldermead",NRandom: int = if CUDA and self.CUDA:DP=mtv.projection(x=x,data=self.dataCuda[:,self.distribution==d],solver=solver,NRandom=NRandom,option=option, n_refinements=n_refinements,sphcap_shrink=sphcap_shrink, alpha_Dirichlet=alpha_Dirichlet,cooling_factor=cooling_factor,cap_size=cap_size,start=start, - space=space,line_solver=line_solver,bound_gc=bound_gc,CUDA=CUDA,device=self.device, + space=space,line_solver=line_solver,bound_gc=bound_gc,CUDA=CUDA,device=self.device,seed=self.seed ) else:DP=mtv.projection(x=x,data=self.data[self.distribution==d],solver=solver,NRandom=NRandom,option=option, n_refinements=n_refinements,sphcap_shrink=sphcap_shrink, alpha_Dirichlet=alpha_Dirichlet,cooling_factor=cooling_factor,cap_size=cap_size,start=start, - space=space,line_solver=line_solver,bound_gc=bound_gc,CUDA=CUDA + space=space,line_solver=line_solver,bound_gc=bound_gc,CUDA=CUDA,seed=self.seed ) if evaluate_dataset==False: if option==1:self.projectionDepth[ind]=DP # assign value diff --git a/depth/model/multivariate/Aprojection.py b/depth/model/multivariate/Aprojection.py index 9af3894..3d725f5 100644 --- a/depth/model/multivariate/Aprojection.py +++ b/depth/model/multivariate/Aprojection.py @@ -19,13 +19,16 @@ def aprojection(x, data, line_solver = "goldensection", bound_gc = True, CUDA=False, - device=None): + device=None, + seed=2801): if CUDA==False: return depth_approximation(x, data, "aprojection", solver, NRandom, option, n_refinements, - sphcap_shrink, alpha_Dirichlet, cooling_factor, cap_size, start, space, line_solver, bound_gc,) + sphcap_shrink, alpha_Dirichlet, cooling_factor, cap_size, start, space, line_solver, + bound_gc,seed=seed) if CUDA==True: - return cudaApprox(data,x, "aprojection", solver, option,NRandom, n_refinements, sphcap_shrink,device=device) + return cudaApprox(data,x, "aprojection", solver, option,NRandom, n_refinements, sphcap_shrink, + device=device, seed=seed) aprojection.__doc__=""" diff --git a/depth/model/multivariate/CUDA_approximation.py b/depth/model/multivariate/CUDA_approximation.py index 300bcb7..c82ece5 100644 --- a/depth/model/multivariate/CUDA_approximation.py +++ b/depth/model/multivariate/CUDA_approximation.py @@ -13,10 +13,10 @@ # device = torch.device("cpu") def cudaApprox(data:torch.Tensor,x:torch.Tensor,notion:str, solver:str,option:int,NRandom:int,n_refinements:int,sphcap_shrink:float, - step:int=10000,device="cpu")->torch.Tensor: + step:int=10000,device="cpu", seed=2801)->torch.Tensor: """Main function to compute approximated depth based on chosen notion """ - torch.manual_seed(2801) + torch.manual_seed(seed) # IMPORTANT TO REMEMBER: data is a transposed matrix, spaceDim x nSamples if len(x.shape)==1:x=x.reshape(1,-1) # xCUDA=torch.tensor(x,dtype=torch.float32,device=device) # transfert x to cuda diff --git a/depth/model/multivariate/Cexpchull.py b/depth/model/multivariate/Cexpchull.py index 2cba237..49be065 100644 --- a/depth/model/multivariate/Cexpchull.py +++ b/depth/model/multivariate/Cexpchull.py @@ -16,10 +16,11 @@ def cexpchull(x, data, start = "mean", space = "sphere", line_solver = "goldensection", - bound_gc = True): + bound_gc = True, + seed=2801): return depth_approximation(x, data, "cexpchull", solver, NRandom, option, n_refinements, - sphcap_shrink, alpha_Dirichlet, cooling_factor, cap_size, start, space, line_solver, bound_gc) + sphcap_shrink, alpha_Dirichlet, cooling_factor, cap_size, start, space, line_solver, bound_gc,seed) cexpchull.__doc__=""" diff --git a/depth/model/multivariate/Cexpchullstar.py b/depth/model/multivariate/Cexpchullstar.py index 0495266..3e165e8 100644 --- a/depth/model/multivariate/Cexpchullstar.py +++ b/depth/model/multivariate/Cexpchullstar.py @@ -16,10 +16,11 @@ def cexpchullstar(x, data, start = "mean", space = "sphere", line_solver = "goldensection", - bound_gc = True): + bound_gc = True, + seed=2801): return depth_approximation(x, data, "cexpchullstar", solver, NRandom, option, n_refinements, - sphcap_shrink, alpha_Dirichlet, cooling_factor, cap_size, start, space, line_solver, bound_gc) + sphcap_shrink, alpha_Dirichlet, cooling_factor, cap_size, start, space, line_solver, bound_gc,seed) cexpchullstar.__doc__=""" diff --git a/depth/model/multivariate/Depth_approximation.py b/depth/model/multivariate/Depth_approximation.py index 02e6f11..e100128 100644 --- a/depth/model/multivariate/Depth_approximation.py +++ b/depth/model/multivariate/Depth_approximation.py @@ -19,7 +19,8 @@ def depth_approximation(z, start = "mean", space = "sphere", line_solver = "goldensection", - bound_gc = True): + bound_gc = True, + seed = 2801): depth_indice = check_depth(notion) check_space(space) @@ -27,6 +28,7 @@ def depth_approximation(z, start_indice = check_start(start) line_solver_indice = check_line_solver(line_solver) check_bound(bound_gc) + try: n, d = X.shape @@ -80,7 +82,7 @@ def depth_approximation(z, objects=(c_double*len(objects_list))(*objects_list) points=pointer(points) objects=pointer(objects) - + seed = pointer((c_int(seed))) libApprox.depth_approximation( objects, @@ -104,7 +106,8 @@ def depth_approximation(z, c_void_p(depths_iter.ctypes.data), c_void_p(directions.ctypes.data), c_void_p(directions_card.ctypes.data), - c_void_p(best_directions.ctypes.data) + c_void_p(best_directions.ctypes.data), + seed ) if(option == 2 or option == 3 or option == 4): diff --git a/depth/model/multivariate/Geometrical.py b/depth/model/multivariate/Geometrical.py index 5642dc9..9ae429e 100644 --- a/depth/model/multivariate/Geometrical.py +++ b/depth/model/multivariate/Geometrical.py @@ -16,10 +16,11 @@ def geometrical(x, data, start = "mean", space = "sphere", line_solver = "goldensection", - bound_gc = True): + bound_gc = True, + seed=2801): return depth_approximation(x, data, "geometrical", solver, NRandom, option, n_refinements, - sphcap_shrink, alpha_Dirichlet, cooling_factor, cap_size, start, space, line_solver, bound_gc) + sphcap_shrink, alpha_Dirichlet, cooling_factor, cap_size, start, space, line_solver, bound_gc,seed) geometrical.__doc__=""" diff --git a/depth/model/multivariate/Halfspace.py b/depth/model/multivariate/Halfspace.py index 55946ac..a8b05a0 100644 --- a/depth/model/multivariate/Halfspace.py +++ b/depth/model/multivariate/Halfspace.py @@ -20,7 +20,8 @@ def halfspace(x, data, exact=True, method="recursive", line_solver = "goldensection", bound_gc = True, CUDA=False, - device=None): + device=None, + seed=2801): if exact: if (method =="recursive" or method==1): method=1 @@ -56,10 +57,10 @@ def halfspace(x, data, exact=True, method="recursive", return res else: if CUDA==False:return depth_approximation(x, data, "halfspace", solver, NRandom ,option, n_refinements, - sphcap_shrink, alpha_Dirichlet, cooling_factor, cap_size, start, space, line_solver, bound_gc) + sphcap_shrink, alpha_Dirichlet, cooling_factor, cap_size, start, space, line_solver, bound_gc,seed) if CUDA==True: return cudaApprox(data,x, "halfspace", solver=solver, option=option,NRandom=NRandom, n_refinements=n_refinements, - sphcap_shrink=sphcap_shrink,device=device) + sphcap_shrink=sphcap_shrink,device=device,seed=seed) halfspace.__doc__=""" diff --git a/depth/model/multivariate/Mahalanobis.py b/depth/model/multivariate/Mahalanobis.py index 98cd118..099d894 100644 --- a/depth/model/multivariate/Mahalanobis.py +++ b/depth/model/multivariate/Mahalanobis.py @@ -20,7 +20,8 @@ def mahalanobis(x, data, exact=True, mah_estimate="moment", mah_parMcd = 0.75, space = "sphere", line_solver = "goldensection", bound_gc = True, - covMCD=None): + covMCD=None, + seed=2801): if exact: points_list=data.flatten() @@ -50,7 +51,7 @@ def mahalanobis(x, data, exact=True, mah_estimate="moment", mah_parMcd = 0.75, return res else: return depth_approximation(x, data, "mahalanobis", solver, NRandom, option, n_refinements, - sphcap_shrink, alpha_Dirichlet, cooling_factor, cap_size, start, space, line_solver, bound_gc) + sphcap_shrink, alpha_Dirichlet, cooling_factor, cap_size, start, space, line_solver, bound_gc,seed) mahalanobis.__doc__= """ diff --git a/depth/model/multivariate/Projection.py b/depth/model/multivariate/Projection.py index 831c8e2..50e67ba 100644 --- a/depth/model/multivariate/Projection.py +++ b/depth/model/multivariate/Projection.py @@ -19,13 +19,14 @@ def projection(x, data, line_solver = "goldensection", bound_gc = True, CUDA=False, - device=None): + device=None, + seed=2801): if CUDA==False: #check cuda return depth_approximation(x, data, "projection", solver, NRandom, option, n_refinements, - sphcap_shrink, alpha_Dirichlet, cooling_factor, cap_size, start, space, line_solver, bound_gc) + sphcap_shrink, alpha_Dirichlet, cooling_factor, cap_size, start, space, line_solver, bound_gc,seed) else: return cudaApprox(data,x, "projection", solver, option,NRandom, n_refinements, - sphcap_shrink,device=device,) # return for depth cuda + sphcap_shrink,device=device,seed=seed) # return for depth cuda projection.__doc__=""" diff --git a/depth/model/multivariate/Zonoid.py b/depth/model/multivariate/Zonoid.py index e581302..2e49bc1 100644 --- a/depth/model/multivariate/Zonoid.py +++ b/depth/model/multivariate/Zonoid.py @@ -16,7 +16,7 @@ def zonoid(x, data, seed=0, exact=True, solver="neldermead", start="mean", space="sphere", line_solver="goldensection", - bound_gc=True): + bound_gc=True,): if exact: points_list=data.flatten() objects_list=x.flatten() @@ -39,7 +39,7 @@ def zonoid(x, data, seed=0, exact=True, solver="neldermead", return res else: return depth_approximation(x, data, "zonoid", solver, NRandom, option, n_refinements, - sphcap_shrink, alpha_Dirichlet, cooling_factor, cap_size, start, space, line_solver, bound_gc) + sphcap_shrink, alpha_Dirichlet, cooling_factor, cap_size, start, space, line_solver, bound_gc,seed) zonoid.__doc__= """ diff --git a/depth/src/ACA_wrapper.cpp b/depth/src/ACA_wrapper.cpp index d82645b..3cab377 100644 --- a/depth/src/ACA_wrapper.cpp +++ b/depth/src/ACA_wrapper.cpp @@ -58,11 +58,22 @@ int SetDepthPars(cProjection& depthObj, int n_refinements, } +void setSeed(int random_seed){ + if (random_seed != 0) { + std::seed_seq seq{random_seed}; + } + else { + std::seed_seq seq{time(NULL)}; + } +} + void ACA(double *z, double *x, int notion, int solver, int NRandom, int n_refinements, double sphcap_shrink, double alpha_Dirichlet, double cooling_factor, double cap_size, int start, int line_solver, int bound_gc, int n, int d, - int n_z, double *depths, double *best_directions, double *basis_py, int d_aca, int option){ + int n_z, double *depths, double *best_directions, double *basis_py, int d_aca, int option,int *seed){ + setSeed(*seed); + unsigned int seed_un = static_cast(*seed); dyMatrixClass::cMatrix X(n ,d); // Fill data matrice with numpy array for(int i=0; i uniDepths[] // constructor -cProjection::cProjection(const cMatrix& x, int n, int d, int NRandom) - : x{ x }, n{ n }, d{ d }, _nProjections{ 0 }, _nRandom{ NRandom }, gen{ _seed }, rnd{ 0.0, 1.0 }, +cProjection::cProjection(const cMatrix& x, int n, int d, int NRandom, unsigned int seed) + : x{ x }, n{ n }, d{ d }, _nProjections{ 0 }, _nRandom{ NRandom }, gen{ seed }, rnd{ 0.0, 1.0 }, xp{ new double[n] {} } { // initialization of array 'Method' that contains the different approximations methods diff --git a/depth/src/ProjectionDepths.h b/depth/src/ProjectionDepths.h index 8c55305..40ca3c9 100644 --- a/depth/src/ProjectionDepths.h +++ b/depth/src/ProjectionDepths.h @@ -148,7 +148,7 @@ class cProjection { std::function MultiDepth; public: // constructor - cProjection(const dyMatrixClass::cMatrix& x, int n = 0, int d = 0, int NRandom = 1000); + cProjection(const dyMatrixClass::cMatrix& x, int n = 0, int d = 0, int NRandom = 1000, unsigned int seed = 2801); // destructor ~cProjection(); // setter methods for most of the private members declared above diff --git a/depth/src/depth_wrapper.cpp b/depth/src/depth_wrapper.cpp index a5ab6d5..5598275 100644 --- a/depth/src/depth_wrapper.cpp +++ b/depth/src/depth_wrapper.cpp @@ -19,6 +19,16 @@ extern "C" { #endif +void setSeed(int random_seed){ + if (random_seed != 0) { + std::seed_seq seq{random_seed}; + } + else { + std::seed_seq seq{time(NULL)}; + } +} + + int SetDepthPars(cProjection& depthObj, int n_refinements, double sphcap_shrink, double alpha_Dirichlet, double cooling_factor, double cap_size, @@ -66,16 +76,17 @@ int depth_approximation(double *z, double *x, int notion, int solver, int line_solver, int bound_gc, int n, int d, int n_z, double *depths, double *depths_iter, double *directions, int *directions_card, - double *best_directions){ + double *best_directions,int *seed){ + setSeed(*seed); dyMatrixClass::cMatrix X(n ,d); // Fill data matrice with numpy array for(int i=0; i(*seed); + cProjection Proj(X, n, d, NRandom, seed_un); Proj.SetDepthNotion((eDepth)notion); Proj.SetMethod((eProjMeth)solver); Proj.SetD_ACA(d);