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..4eacd59 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. @@ -211,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 @@ -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, 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 @@ -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. @@ -326,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 @@ -366,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 @@ -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 """ @@ -481,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 @@ -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. @@ -545,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 @@ -584,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 @@ -632,7 +649,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. @@ -645,18 +662,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 @@ -685,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 @@ -734,7 +752,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. @@ -747,18 +765,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 @@ -790,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 @@ -835,7 +853,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 +863,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 +905,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 +948,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. @@ -936,18 +961,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 """ @@ -978,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 @@ -1022,7 +1047,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 +1083,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 +1124,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 +1135,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 +1162,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 +1173,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 +1211,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. @@ -1187,18 +1224,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 @@ -1264,7 +1301,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,27 +1331,50 @@ 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): + """ + 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): @@ -1323,32 +1382,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 @@ -1358,17 +1416,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: """ @@ -1403,7 +1463,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: @@ -1412,10 +1472,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: @@ -1496,9 +1557,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/DepthFunc.py b/depth/model/DepthFunc.py index 1985bf6..04e2607 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. @@ -545,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 ----- @@ -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]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 e73dab4..e100128 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, @@ -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,14 +106,15 @@ 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): 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/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/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..099d894 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,9 @@ 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, + seed=2801): if exact: points_list=data.flatten() @@ -39,8 +37,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)) @@ -53,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/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/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/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/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 new file mode 100644 index 0000000..63b7814 --- /dev/null +++ b/depth/model/multivariate/isInConvexes.py @@ -0,0 +1,79 @@ +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, d = X.shape[0], 1 + n_z = z.shape[0] + + + 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/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 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/ProjectionDepths.cpp b/depth/src/ProjectionDepths.cpp index ce0ad3d..c363823 100644 --- a/depth/src/ProjectionDepths.cpp +++ b/depth/src/ProjectionDepths.cpp @@ -43,7 +43,7 @@ using namespace dyMatrixClass; const int debug = 1; // seed for the RNG -#define _seed 1234 +// #define _seed 1234 // error indicator used when computing the zonoid depth int error; @@ -80,8 +80,8 @@ function 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/ddalpha.cpp b/depth/src/ddalpha.cpp index 7e271f7..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. @@ -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++){ @@ -62,6 +70,7 @@ void IsInConvexes(double *points, int *dimension, int *cardinalities, int *numCl } 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++){ @@ -69,6 +78,8 @@ void IsInConvexes(double *points, int *dimension, int *cardinalities, int *numCl } } + + void ZDepth(double *points, double *objects, int *numPoints, int *numObjects, int *dimension, int *seed, double *depths){ setSeed(*seed); TMatrix x(numPoints[0]); @@ -92,7 +103,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, @@ -262,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; diff --git a/depth/src/ddalpha.dll b/depth/src/ddalpha.dll index c197ec2..bfc6886 100644 Binary files a/depth/src/ddalpha.dll and b/depth/src/ddalpha.dll differ 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); diff --git a/setup.py b/setup.py index 3e9e4f2..1814b93 100644 --- a/setup.py +++ b/setup.py @@ -62,6 +62,7 @@ def build_extensions(self): # 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",],