miutils.c

互信息盲源分离

  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*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));
  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);
}

void red(double **x, int dim, int N, int K, 
	 double *psi,
	 double *scal,
	 double *mic, double *mir) {
  int i,k,*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 *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*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));

  make_box2(x,dim,N,0,dim-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);
    //    for (k=1;k<=K;k++) {
    //      fprintf(stdout,"k%d %f %f %f %f\n",k,x[0][i],x[0][nn[k]],x[1][i],x[1][nn[k]]);
    //    }
    
    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]]; }
      else { 
	n1[d]=neiE1(x[d],i,scal[d],BOX1,Eps,box1[d],lis1[d],mxi[d]); 
	dxy1+=psi[n1[d]+1];
      }
      dxy2+=psi[n2[d]];
    }
    //    fprintf(stdout,"%f %f %e %e %e %d %d %d %d %f %d\n",xc[0],xc[1],eps[0],eps[1],(eps[0]-eps[1]),n1[0],n1[1],n2[0],n2[1],Eps,nn[K]);
  }
  dxy1/=N;*mic=(dim-1)*psi[N]+psi[K]-dxy1;
  dxy2/=N;*mir=(dim-1)*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(n1);free(n2);
  free(phi);
}


void redc(double **x, int dim, int N, int K, 
	  double *psi,
	  double *scal,
	  double *mic) {
  
  int i,k,*n1;
  double *xc,dx;
  double *eps,Eps;
  int maxdim;
  double dxy1;
  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;

  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));

  make_box2(x,dim,N,0,dim-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=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++) {
      n1[d]=neiE1(x[d],i,scal[d],BOX1,Eps,box1[d],lis1[d],mxi[d]); 
      if (d==maxdim) { 
	dxy1+=psi[n1[d]];
      } else { 
	dxy1+=psi[n1[d]+1];
      }
    }
    //    fprintf(stdout,"%f %f %d %d %d %f %d\n",xc[0],xc[1],(eps[0]<=eps[1]),n1[0],n1[1],Eps,nn[K]);
  }
  dxy1/=N;*mic=(dim-1)*psi[N]+psi[K]-dxy1;
  
  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);
}

  


void redr(double **x, int dim, 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;

  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*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));
  n2=(int*)calloc(dim,sizeof(int));

  make_box2(x,dim,N,0,dim-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=(dim-1)*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);
}
  

void redr_embed(double **x, int dim, int edim, int tau, int N, int K, 
		double *psi,
		double *mir) {
  int i,k,*n2;
  double *xc,dx;
  double *eps;
  double dxy2;
  int *nn;
  int d,de;
  int N_;

  double **xx;
  double ***xxx;
  
  int BOX;
  int **box,*lis; // two dimensional boxes
  int ***boxxx; // twodimensional boxes for each dim
  int **lisxx; // lists for each dim
  double epsilon;
  int inveps;

  int *ind;
  int **indxx;

  double *phi;

  //  if (edim<2) Just call usual redundancy 

  N_=N-(edim-1)*tau;
  xx=(double **)calloc(dim*edim,sizeof(double*));
  // reservation, and embedding
  for (d=0; d<dim*edim; d++) {
    xx[d]=(double*)calloc(N_,sizeof(double));
    memcpy(xx[d],x[d/edim]+tau*(d%edim),N_*sizeof(double));
  }
  xxx=(double ***)calloc(dim,sizeof(double**));
  // just copy, make bo will destroy an order
  for (d=0; d<dim; d++) {
    xxx[d]=(double**)calloc(edim,sizeof(double*));
    for (de=0; de<edim; de++) {
      xxx[d][de]=(double*)calloc(N_,sizeof(double));
      memcpy(xxx[d][de],xx[d*edim+de],N_*sizeof(double));
    }
  }

  ind=(int*)calloc(N_,sizeof(int));
  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*edim+1,sizeof(double));
  BOX=1; while (0.5*BOX*BOX*K<N) BOX*=2;
  epsilon=4.0/BOX;
  inveps=BOX/4;

  boxxx=(int***)calloc(dim,sizeof(int**)); 
  lisxx=(int**)calloc(dim,sizeof(int*)); 
  for (d=0;d<dim;d++) {
    boxxx[d]=(int**)calloc(BOX,sizeof(int*)); 
    for (i=0;i<BOX;i++) boxxx[d][i]=(int*)calloc(BOX,sizeof(int));
    lisxx[d]=(int*)calloc(N_,sizeof(int)); 
  }
  indxx=(int**)calloc(dim,sizeof(int*));
  for (d=0;d<dim;d++) {
    indxx[d]=(int*)calloc(N_,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));

  //for searching neighbours in product (dim*edim) space
  make_box2ind(xx,dim*edim,N_,0,edim,BOX,inveps,ind,box,lis); 
  //second component is arbitratry, lets take the first component of second vector 
  // it should be the same for neiK procedure

  for (d=0;d<dim;d++) 
    make_box2ind(xxx[d],edim,N_,0,edim-1,BOX,inveps,indxx[d],boxxx[d],lisxx[d]);
  
  dxy2=0.0;
  for (i=0;i<N_;i++) {
    for (d=0;d<dim*edim;d++) xc[d]=xx[d][ind[i]];
    neiK(xx,dim*edim,0,edim,ind[i],BOX,epsilon,K,box,lis,nn);
    
    for (d=0;d<dim;d++) {
      eps[d]=0;
      for (de=0;de<edim;de++) 
	for(k=1;k<=K;k++) 
	  if( (dx=fabs(xc[d*edim+de]-xx[d*edim+de][nn[k]]))>eps[d] ) eps[d]=dx;
    }	
    for (d=0;d<dim;d++) {
      n2[d]=neiE(xxx[d],indxx[d][i],0,edim-1,edim,BOX,epsilon,eps[d],boxxx[d],lisxx[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=(dim-1)*psi[N_]+phi[K]-dxy2;

  for (d=0; d<dim*edim; d++) free(xx[d]);
  free(xx);

  for (d=0; d<dim; d++) {
    for (de=0; de<edim; de++) free(xxx[d][de]);
    for (i=0;i<BOX;i++) free(boxxx[d][i]);
    free(boxxx[d]);free(lisxx[d]);free(xxx[d]);free(indxx[d]);
  }
  free(xxx);free(boxxx);free(lisxx);free(indxx);
  free(ind);free(phi);free(nn);
  free(xc);

  for (i=0;i<BOX;i++) free(box[i]);
  free(box);free(lis);
  free(eps);free(n2);
}
  

void mi_xnyn(double **x, int dimx, int dimy, int N, int K, 
	     double *psi, 
	     double *scal,
	     double *mic, double *mir) {
  int i,k,nx1,ny1,nx2,ny2;
  double *xc,dy,dx;
  double epsx,epsy,Eps;
  double dxy1,dxy2;
  double **xx,**yy;;
  double scalx, scaly;
  int *nn;
  int d;

  int BOX,BOX1;
  int **box,**boxy,**boxx,*lis; // two dimensional boxes
  int *lisy,*lisx; // lists for two dimensions
  int *boxx1, *boxy1; // onedimensional boxes
  int *lisy1,*lisx1; // lists for one dimensions
  int *mxi, *myi; //accumulative lists of points in oned boxes
  int *ind,*indx,*indy; //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(dimx+dimy,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); 

  dxy1=dxy2=0.0;
  for (i=0;i<N;i++) {