miutils.c

互信息盲源分离

  Eps=(double*)calloc(K,sizeof(double));

  BOX=1; while (0.5*BOX*BOX*K<N) BOX*=2;
  epsilon=4.0/BOX;
  inveps=BOX/4;
  BOX1=N-5;

  if (dimx>1) {
    xx=(double**)calloc(dimx,sizeof(double*));
    for (d=0;d<dimx;d++) xx[d]=(double*)calloc(N,sizeof(double));
    boxx=(int**)calloc(BOX,sizeof(int*)); 
    for (i=0;i<BOX;i++) boxx[i]=(int*)calloc(BOX,sizeof(int));
    lisx=(int*)calloc(N,sizeof(int));
  } else { boxx1=(int*)calloc(BOX1+1,sizeof(int)); 
  lisx1=(int*)calloc(N,sizeof(int)); mxi=(int*)calloc(BOX1+1,sizeof(int)); }
  if (dimy>1) {
    yy=(double**)calloc(dimy,sizeof(double*));
    for (d=0;d<dimy;d++) yy[d]=(double*)calloc(N,sizeof(double));
    boxy=(int**)calloc(BOX,sizeof(int*)); 
    for (i=0;i<BOX;i++) boxy[i]=(int*)calloc(BOX,sizeof(int));
    lisy=(int*)calloc(N,sizeof(int));
  } else { boxy1=(int*)calloc(BOX1+1,sizeof(int)); 
  lisy1=(int*)calloc(N,sizeof(int)); myi=(int*)calloc(BOX1+1,sizeof(int)); }
 
  box=(int**)calloc(BOX,sizeof(int*));
  for (i=0;i<BOX;i++) box[i]=(int*)calloc(BOX,sizeof(int));
  lis=(int*)calloc(N,sizeof(int));

  ind=(int*)calloc(N,sizeof(int));
  indx=(int*)calloc(N,sizeof(int));
  indy=(int*)calloc(N,sizeof(int));

  //save x if it would be reordered
  if (dimx>1) for (d=0;d<dimx;d++) memcpy(xx[d],x[d],N*sizeof(double)); 
  if (dimy>1) for (d=0;d<dimy;d++) memcpy(yy[d],x[d+dimx],N*sizeof(double));

  make_box2ind(x,dimx+dimy,N,0,dimx,BOX,inveps,ind,box,lis); 
  //for searching neighbours in product space

  if (dimx==1) {scalx=scal[0]; make_box1(x[0],N,scalx,BOX1,boxx1,lisx1,mxi);}
  else make_box2ind(xx,dimx,N,0,dimx-1,BOX,inveps,indx,boxx,lisx); 
  if (dimy==1) {scaly=scal[dimx]; make_box1(x[dimx],N,scaly,BOX1,boxy1,lisy1,myi); }
  else make_box2ind(yy,dimy,N,0,dimy-1,BOX,inveps,indy,boxy,lisy); 

  for (ik=0;ik<K;ik++) {dxy1[ik]=dxy2[ik]=0.0;}
  for (i=0;i<N;i++) {
    for (d=0;d<dimx+dimy;d++) xc[d]=x[d][ind[i]];

    neiK(x,dimx+dimy,0,dimx,ind[i],BOX,epsilon,K,box,lis,nn);
    for (ik=0;ik<K;ik++) {
      epsx[ik]=0; for (d=0;d<dimx;d++) 
	for(k=1;k<=ik+1;k++) if( (dx=fabs(xc[d]-x[d][nn[k]]))>epsx[ik] ) epsx[ik]=dx;
      epsy[ik]=0; for (d=dimx;d<dimx+dimy;d++) 
	for(k=1;k<=ik+1;k++) if( (dy=fabs(xc[d]-x[d][nn[k]]))>epsy[ik] ) epsy[ik]=dy;
      if (dimx>1) nx2[ik]=neiE(xx,indx[i],0,dimx-1,dimx,BOX,epsilon,epsx[ik],boxx,lisx);
      else  nx2[ik]=neiE1(x[0],ind[i],scalx,BOX1,epsx[ik],boxx1,lisx1,mxi);
      if (dimy>1) ny2[ik]=neiE(yy,indy[i],0,dimy-1,dimy,BOX,epsilon,epsy[ik],boxy,lisy);
      else ny2[ik]=neiE1(x[dimx],ind[i],scaly,BOX1,epsy[ik],boxy1,lisy1,myi);
    
      if (epsx[ik]>epsy[ik]) {
	Eps[ik]=epsx[ik];nx1[ik]=nx2[ik];
	if (dimy>1) ny1[ik]=neiE(yy,indy[i],0,dimy-1,dimy,BOX,epsilon,Eps[ik],boxy,lisy);
	else ny1[ik]=neiE1(x[dimx],ind[i],scaly,BOX1,Eps[ik],boxy1,lisy1,myi);
	dxy1[ik]+=psi[nx1[ik]]+psi[ny1[ik]+1];
      } else {
	Eps[ik]=epsy[ik];ny1[ik]=ny2[ik];
	if (dimx>1) nx1[ik]=neiE(xx,indx[i],0,dimx-1,dimx,BOX,epsilon,Eps[ik],boxx,lisx);
	else nx1[ik]=neiE1(x[0],ind[i],scalx,BOX1,Eps[ik],boxx1,lisx1,mxi);
	dxy1[ik]+=psi[nx1[ik]+1]+psi[ny1[ik]];
      }
      dxy2[ik]+=psi[nx2[ik]]+psi[ny2[ik]];
    }
  }
  for (ik=0;ik<K;ik++) {
    dxy1[ik]/=N;mi_cr[2*ik]=psi[N]+psi[ik+1]-dxy1[ik];
    dxy2[ik]/=N;mi_cr[2*ik+1]=psi[N]+phi[ik+1]-dxy2[ik];
  }

  free(xc);free(nn);free(lis);
  for (i=0;i<BOX;i++) free(box[i]); free(box);
  free(ind);free(indx);free(indy);
  if (dimx==1) {free(mxi);free(boxx1);free(lisx1);} 
  else { for (i=0;i<BOX;i++) free(boxx[i]); free(boxx); free(lisx); for (d=0;d<dimx;d++) free(xx[d]); free(xx); }
  if (dimy==1) {free(myi);free(boxy1);free(lisy1);} 
  else { for (i=0;i<BOX;i++) free(boxy[i]); free(boxy); free(lisy); for (d=0;d<dimy;d++) free(yy[d]); free(yy); }
  free(phi);
  free(epsx);free(epsy);free(Eps);
  free(nx1);free(nx2);free(ny1);free(ny2);
  free(dxy1);free(dxy2);
}

void mi_xnynKembed(double **x, int dim, int N, int K, 
		   double **xx, double **yy, 
		   int **boxx, int *lisx, int *indx,
		   int **boxy, int *lisy, int *indy,
		   double *psi, 
		   double *mi_cr) {

  int i,ik,k,*nx1,*ny1,*nx2,*ny2;
  double *xc,dy,dx;
  double *epsx,*epsy,*Eps;
  char *maxdim;
  double *dxy1,*dxy2;
  int *nn;
  int d;

  int BOX;
  int **box,*lis; // two dimensional boxes
  int *ind; //indeces of original data (the data resorted during box creating)
  double epsilon;
  int inveps;

  double *phi;

  phi=(double*)calloc(K+1,sizeof(double));
  for (i=1;i<=K;i++) phi[i]=psi[i]-1/(double(i));   // 

  nn=(int*)calloc(K+1,sizeof(int));

  xc=(double*)calloc(2*dim,sizeof(double));
  dxy1=(double*)calloc(K,sizeof(double));
  dxy2=(double*)calloc(K,sizeof(double));
  nx1=(int*)calloc(K,sizeof(int));
  nx2=(int*)calloc(K,sizeof(int));
  ny1=(int*)calloc(K,sizeof(int));
  ny2=(int*)calloc(K,sizeof(int));
  epsx=(double*)calloc(K,sizeof(double));
  epsy=(double*)calloc(K,sizeof(double));
  Eps=(double*)calloc(K,sizeof(double));
  maxdim=(char*)calloc(K,sizeof(char));

  BOX=1; while (0.5*BOX*BOX*K<N) BOX*=2;
  epsilon=4.0/BOX;
  inveps=BOX/4;

  box=(int**)calloc(BOX,sizeof(int*));
  for (i=0;i<BOX;i++) box[i]=(int*)calloc(BOX,sizeof(int));
  lis=(int*)calloc(N,sizeof(int));

  ind=(int*)calloc(N,sizeof(int));


  make_box2ind(x,2*dim,N,0,dim,BOX,inveps,ind,box,lis); 

  for (ik=0;ik<K;ik++) {dxy1[ik]=dxy2[ik]=0.0;}
  for (i=0;i<N;i++) {
    for (d=0;d<2*dim;d++) xc[d]=x[d][ind[i]];

    neiK(x,2*dim,0,dim,ind[i],BOX,epsilon,K,box,lis,nn);
    for (ik=0;ik<K;ik++) {
      epsx[ik]=0; for (d=0;d<dim;d++) 
	for(k=1;k<=ik+1;k++) if( (dx=fabs(xc[d]-x[d][nn[k]]))>epsx[ik] ) epsx[ik]=dx;
      epsy[ik]=0; for (d=dim;d<2*dim;d++) 
	for(k=1;k<=ik+1;k++) if( (dy=fabs(xc[d]-x[d][nn[k]]))>epsy[ik] ) epsy[ik]=dy;
    }
    neiEK(xx,indx[i],0,dim-1,dim,K,BOX,epsilon,epsx,boxx,lisx,nx2);
    neiEK(yy,indy[i],0,dim-1,dim,K,BOX,epsilon,epsy,boxy,lisy,ny2);
    for (ik=0;ik<K;ik++) {
      maxdim[ik]=0;
      if (epsx[ik]>epsy[ik]) { Eps[ik]=epsx[ik]; } else { Eps[ik]=epsy[ik]; maxdim[ik]=1; }
    }
    neiEK(yy,indy[i],0,dim-1,dim,K,BOX,epsilon,Eps,boxy,lisy,ny1);
    neiEK(xx,indx[i],0,dim-1,dim,K,BOX,epsilon,Eps,boxx,lisx,nx1);

    for (ik=0;ik<K;ik++) {
      if (maxdim[ik]) dxy1[ik]+=psi[nx1[ik]+1]+psi[ny1[ik]]; else dxy1[ik]+=psi[nx1[ik]]+psi[ny1[ik]+1];
      dxy2[ik]+=psi[nx2[ik]]+psi[ny2[ik]];
    }
   
  }
  for (ik=0;ik<K;ik++) {
    dxy1[ik]/=N;mi_cr[2*ik]=psi[N]+psi[ik+1]-dxy1[ik];
    dxy2[ik]/=N;mi_cr[2*ik+1]=psi[N]+phi[ik+1]-dxy2[ik];
  }

  free(xc);free(nn);free(lis);
  for (i=0;i<BOX;i++) free(box[i]); free(box);
  free(ind);
  free(phi);
  free(epsx);free(epsy);free(Eps);free(maxdim);
  free(nx1);free(nx2);free(ny1);free(ny2);
  free(dxy1);free(dxy2);
}



void mi2h(double **x, int N, int K,
	  double *psi, 
	  double *scal,
	  double *mic, double *mir, double *hc, double *hr) {
  
  int i,k,*n1,*n2;
  double *xc,dx;
  double *eps,Eps;
  double *h1,*h2,hj1,hj2;
  int maxdim;
  double dxy1,dxy2;
  int *nn;
  int d;
  
  int BOX,BOX1;
  int **box,*lis; // two dimensional boxes
  int **box1; // onedimensional boxes
  int **lis1; // lists for one dimensions
  int **mxi; //accumulative lists of points in oned boxes
  double epsilon;
  int inveps;

  int dim=2;

  double *phi;
  double t_d1,t_d2;

  phi=(double*)calloc(N+1,sizeof(double));
  for (i=1;i<=N;i++) phi[i]=psi[i]-(dim-1)/(double(i));   // 


  nn=(int*)calloc(K+1,sizeof(int));
  xc=(double*)calloc(dim+1,sizeof(double));
  BOX=1; while (0.5*BOX*BOX*K<N) BOX*=2;
  epsilon=4.0/BOX;
  inveps=BOX/4;

  BOX1=N-5;

  box1=(int**)calloc(dim,sizeof(int*)); 
  lis1=(int**)calloc(dim,sizeof(int*)); 
  mxi=(int**)calloc(dim,sizeof(int*)); 
  for (d=0;d<dim;d++) {
    box1[d]=(int*)calloc(BOX1+1,sizeof(int)); 
    lis1[d]=(int*)calloc(N,sizeof(int)); 
    mxi[d]=(int*)calloc(BOX1+1,sizeof(int));
  }
 
  box=(int**)calloc(BOX,sizeof(int*));
  for (i=0;i<BOX;i++) box[i]=(int*)calloc(BOX,sizeof(int));
  lis=(int*)calloc(N,sizeof(int));

  eps=(double*)calloc(dim,sizeof(double));
  n1=(int*)calloc(dim,sizeof(int));
  n2=(int*)calloc(dim,sizeof(int));

  h1=(double*)calloc(dim,sizeof(double));
  h2=(double*)calloc(dim,sizeof(double));
  for (d=0;d<dim;d++) h1[d]=h2[d]=0;
  hj1=hj2=0;

  make_box2(x,dim,N,0,1,BOX,inveps,box,lis); //for searching neighbours in prodict space
  for (d=0;d<dim;d++) make_box1(x[d],N,scal[d],BOX1,box1[d],lis1[d],mxi[d]);      
  
  dxy1=dxy2=0.0;
  for (i=0;i<N;i++) {
    for (d=0;d<dim;d++) xc[d]=x[d][i];
    neiK(x,dim,0,dim-1,i,BOX,epsilon,K,box,lis,nn);
    
    Eps=0;maxdim=-1;
    for (d=0;d<dim;d++) {
      eps[d]=0;
      for(k=1;k<=K;k++) {if( (dx=fabs(xc[d]-x[d][nn[k]]))>eps[d] ) eps[d]=dx; }
      if (eps[d]>Eps) {Eps=eps[d];maxdim=d;}
    }	
    for (d=0;d<dim;d++) {
      n2[d]=neiE1(x[d],i,scal[d],BOX1,eps[d],box1[d],lis1[d],mxi[d]);
      if (d==maxdim) { 
	n1[d]=n2[d]; dxy1+=psi[n1[d]];
	h1[d]+=log(Eps)-psi[n1[d]];
      }
      else { 
	n1[d]=neiE1(x[d],i,scal[d],BOX1,Eps,box1[d],lis1[d],mxi[d]); 
	dxy1+=psi[n1[d]+1];
	h1[d]+=log(Eps)-psi[n1[d]+1];
      }
      dxy2+=psi[n2[d]];

      h2[d]+=log(eps[d])-phi[n2[d]];
      hj1+=log(Eps);
      hj2+=log(eps[d]);
    }
    //    fprintf(stdout,"%f %f %d %d %d %d %d %f %d\n",xc[0],xc[1],(eps[0]<=eps[1]),n1[0],n1[1],n2[0],n2[1],Eps,nn[K]);
    
  }
  dxy1/=N;*mic=psi[N]+psi[K]-dxy1;
  dxy2/=N;*mir=psi[N]+phi[K]-dxy2;
  hj1/=N; *hc=psi[N]-psi[K]+hj1;//dim-1
  hj2/=N; *hr=psi[N]-phi[K]+hj2;//dim-1
  t_d1=+1e30;
  t_d2=+1e30;
  fprintf(stdout,  "hj_c\t%1.3f\t\t\t\thj_r\t%1.3f\n",*hc,*hr);
  for (d=0;d<dim;d++) {
    h1[d]/=N;h1[d]+=psi[N];
    h2[d]/=N;h2[d]+=psi[N];
    if (h1[d]<t_d1) t_d1=h1[d];
    if (h2[d]<t_d2) t_d2=h2[d];
    if ((*hc-h1[d])<t_d1) t_d1=*hc-h1[d];
    if ((*hr-h2[d])<t_d2) t_d2=*hr-h2[d];
    fprintf(stdout,"h[%d]_c\t%1.3f\thcond[%d]_c\t%1.3f\th[%d]_r\t%1.3f\thcond[%d]_r\t%1.3f\n",
	    d,h1[d],d,*hc-h1[d],d,h2[d],d,*hr-h2[d]);
  }
  *hc-=2*t_d1;
  *hr-=2*t_d2;
  
  free(xc);free(nn);
  for (i=0;i<BOX;i++) free(box[i]); free(box);
  free(lis);
  for (d=0;d<dim;d++) {
    free(box1[d]);free(lis1[d]);free(mxi[d]);
  }
  free(box1);free(lis1);free(mxi);
  free(eps);free(n1);free(n2);
  free(phi);
  free(h1);free(h2);
}










void mi_d(double **x, int dim, int N, int K, float *mi_cr,
	  double *psi, double *phi, double minx, double maxx, double miny, double maxy) {
  
  int i,k,ik,**n1,**n2;
  double *xc,dx;
  double **eps,*Eps;
  int *maxdim;
  double *dxy1,*dxy2;
  int *nn;
  int d;
  
  int BOX,BOX1;
  int **box,*lis; // two dimensional boxes
  int **box1; // onedimensional boxes
  int **lis1; // lists for one dimensions
  int **mxi; //accumulative lists of points in oned boxes
  double epsilon;
  int inveps;

  double scal[2];

  nn=(int*)calloc(K+1,sizeof(int));
  xc=(double*)calloc(dim+1,sizeof(double));
  BOX=1; while (0.5*BOX*BOX*K<N) BOX*=2;
  epsilon=4.0/BOX;
  inveps=BOX/4;

  BOX1=N-5;
  scal[0]=BOX1/(maxx-minx);
  scal[1]=BOX1/(maxy-miny);

  box1=(int**)calloc(dim,sizeof(int*)); 
  lis1=(int**)calloc(dim,sizeof(int*)); 
  mxi=(int**)calloc(dim,sizeof(int*)); 
  for (d=0;d<dim;d++) {
    box1[d]=(int*)calloc(BOX1+1,sizeof(int)); 
    lis1[d]=(int*)calloc(N,sizeof(int)); 
    mxi[d]=(int*)calloc(BOX1+1,sizeof(int));
  }
 
  box=(int**)calloc(BOX,sizeof(int*));
  for (i=0;i<BOX;i++) box[i]=(int*)calloc(BOX,sizeof(int));
  lis=(int*)calloc(N,sizeof(int));

  eps=(double**)calloc(K,sizeof(double*));
  Eps=(double*)calloc(K,sizeof(double));
  n1=(int**)calloc(K,sizeof(int*));
  n2=(int**)calloc(K,sizeof(int*));
  maxdim=(int*)calloc(K,sizeof(int));

  for (ik=0;ik<K;ik++) {
    eps[ik]=(double*)calloc(dim,sizeof(double));
    n1[ik]=(int*)calloc(dim,sizeof(int));
    n2[ik]=(int*)calloc(dim,sizeof(int));
    mi_cr[ik*2]=0;
    mi_cr[ik*2+1]=0;
  }

  dxy1=(double*)calloc(K,sizeof(double));
  dxy2=(double*)calloc(K,sizeof(double));

  make_box2(x,dim,N,0,1,BOX,inveps,box,lis); //for searching neighbours in prodict space
  for (d=0;d<dim;d++) {
    make_box1(x[d],N,scal[d&1],BOX1,box1[d],lis1[d],mxi[d]);      
  }
  
  for (ik=0;ik<K;ik++) {dxy1[ik]=dxy2[ik]=0.0;}
  for (i=0;i<N;i++) {
    for (d=0;d<dim;d++) xc[d]=x[d][i];
    neiK(x,dim,0,1,i,BOX,epsilon,K,box,lis,nn);
    
    for (ik=0;ik<K;ik++) {
      Eps[ik]=0;maxdim[ik]=-1;
      for (d=0;d<dim;d++) {
	eps[ik][d]=0;
	for(k=1;k<=ik+1;k++) {if( (dx=fabs(xc[d]-x[d][nn[k]]))>eps[ik][d] ) eps[ik][d]=dx; }
	if (eps[ik][d]>Eps[ik]) {Eps[ik]=eps[ik][d];maxdim[ik]=d;}
      }	
    }
    for (ik=0;ik<K;ik++) {
      for (d=0;d<dim;d++) {
	n2[ik][d]=neiE1(x[d],i,scal[d&1],BOX1,eps[ik][d],box1[d],lis1[d],mxi[d]);
	if (d==maxdim[ik]) { n1[ik][d]=n2[ik][d]; dxy1[ik]+=psi[n1[ik][d]]; }
	else { 
	  n1[ik][d]=neiE1(x[d],i,scal[d&1],BOX1,Eps[ik],box1[d],lis1[d],mxi[d]); 
	  dxy1[ik]+=psi[n1[ik][d]+1];
	}
	dxy2[ik]+=psi[n2[ik][d]];
      }
    }
  }
  for (ik=0;ik<K;ik++) {
    dxy1[ik]/=N;mi_cr[ik*2]=float((dim-1)*psi[N]+psi[ik+1]-dxy1[ik]);
    dxy2[ik]/=N;mi_cr[ik*2+1]=float((dim-1)*psi[N]+phi[ik+1]-dxy2[ik]);
  }

  free(xc);free(nn);
  for (i=0;i<BOX;i++) free(box[i]); free(box);
  free(lis);
  for (d=0;d<dim;d++) {
    free(box1[d]);free(lis1[d]);free(mxi[d]);
  }
  free(box1);free(lis1);free(mxi);
  
  for (ik=0;ik<K;ik++) {
    free(eps[ik]);
    free(n1[ik]);
    free(n2[ik]);
  }
  free(eps);free(n1);free(n2);
  free(Eps);free(maxdim);

  free(dxy1); free(dxy2);
}




void mi2r_(double **x, int N, int K,
	  double *psi, 
	  double *scal,
	  double *mir) {
  
  int i,k,*n2;
  double *xc,dx;
  double *eps;
  double dxy2;
  int *nn;
  int d;
  
  int BOX,BOX1;
  int **box,*lis; // two dimensional boxes
  int **box1; // onedimensional boxes
  int **lis1; // lists for one dimensions
  int **mxi; //accumulative lists of points in oned boxes
  double epsilon;
  int inveps;

  int dim=2;

  double *phi;

  phi=(double*)calloc(K+1,sizeof(double));
  for (i=1;i<=K;i++) phi[i]=psi[i]-(dim-1)/(double(i));   // 


  nn=(int*)calloc(K+1,sizeof(int));
  xc=(double*)calloc(dim+1,sizeof(double));
  BOX=1; while (0.5*BOX*BOX<N) BOX*=2;
  epsilon=4.0/BOX;
  inveps=BOX/4;

  BOX1=N-5;

  box1=(int**)calloc(dim,sizeof(int*)); 
  lis1=(int**)calloc(dim,sizeof(int*)); 
  mxi=(int**)calloc(dim,sizeof(int*)); 
  for (d=0;d<dim;d++) {
    box1[d]=(int*)calloc(BOX1+1,sizeof(int)); 
    lis1[d]=(int*)calloc(N,sizeof(int)); 
    mxi[d]=(int*)calloc(BOX1+1,sizeof(int));
  }
 
  box=(int**)calloc(BOX,sizeof(int*));
  for (i=0;i<BOX;i++) box[i]=(int*)calloc(BOX,sizeof(int));
  lis=(int*)calloc(N,sizeof(int));

  eps=(double*)calloc(dim,sizeof(double));
  n2=(int*)calloc(dim,sizeof(int));

  make_box2(x,dim,N,0,1,BOX,inveps,box,lis); //for searching neighbours in prodict space
  for (d=0;d<dim;d++) make_box1(x[d],N,scal[d],BOX1,box1[d],lis1[d],mxi[d]);      
  
  dxy2=0.0;
  for (i=0;i<N;i++) {
    for (d=0;d<dim;d++) xc[d]=x[d][i];
    neiK(x,dim,0,dim-1,i,BOX,epsilon,K,box,lis,nn);
    
    for (d=0;d<dim;d++) {
      eps[d]=0;
      for(k=1;k<=K;k++) {if( (dx=fabs(xc[d]-x[d][nn[k]]))>eps[d] ) eps[d]=dx; }
    }	
    for (d=0;d<dim;d++) {
      n2[d]=neiE1(x[d],i,scal[d],BOX1,eps[d],box1[d],lis1[d],mxi[d]);
      dxy2+=psi[n2[d]];
    }
    //    fprintf(stdout,"%f %f %d %d %d %d\n",xc[0],xc[1],(eps[0]<=eps[1]),n2[0],n2[1],nn[K]);
  }
  dxy2/=N;*mir=psi[N]+phi[K]-dxy2;
  
  free(xc);free(nn);
  for (i=0;i<BOX;i++) free(box[i]); free(box);
  free(lis);
  for (d=0;d<dim;d++) {
    free(box1[d]);free(lis1[d]);free(mxi[d]);
  }
  free(box1);free(lis1);free(mxi);
  free(eps);free(n2);
  free(phi);
}