miutils.c

互信息盲源分离

#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
void make_box1(double *x, int N, double scal, int bs, 
	       int *box, int *lis, int *mxi) {
  int i, ix;
  for (i=0;i<=bs;i++) box[i]=-1;
  for (i=0;i<=bs;i++) mxi[i]=0;
  for(i=0;i<N;i++) {ix=(int)(x[i]*scal);lis[i]=box[ix]; box[ix]=i; mxi[ix]++;}
  for (i=1;i<=bs;i++) mxi[i]+=mxi[i-1];
}
void make_box2(double **x, int dim, int N, int comp1, int comp2, int bs, int inveps, 
	       int **box, int *lis) {
  int d,i,ix,iy,ixy;
  int ib=bs-1;
  double **xx;

  xx=(double **)calloc(dim, sizeof(double*)); 
  for (d=0;d<dim;d++) {
    xx[d]=(double *)calloc(N, sizeof(double)); 
    memcpy(xx[d],x[d],N*sizeof(double));
  }
  for (ix=0;ix<bs;ix++) for (iy=0;iy<bs;iy++) box[ix][iy] = -1;
  
  for (i=0;i<N;i++) {
    ix=(int)(x[comp1][i]*inveps)&ib;iy=(int)(x[comp2][i]*inveps)&ib;
    lis[i]=box[ix][iy];box[ix][iy]=i;
  }
  i=-1;
  for (ix=0;ix<bs;ix++) for (iy=0;iy<bs;iy++) {
    ixy=box[ix][iy];
    while(ixy>=0) {
      i++;
      for (d=0;d<dim;d++) x[d][i]=xx[d][ixy];
      ixy=lis[ixy];
    }
    box[ix][iy]=-1;
  }
  for (i=0;i<N;i++) {
    ix=(int)(x[comp1][i]*inveps)&ib;iy=(int)(x[comp2][i]*inveps)&ib;
    lis[i]=box[ix][iy];box[ix][iy]=i;
  }
  for (d=0;d<dim;d++) free(xx[d]); free(xx);
}
void make_box2ind(double **x, int dim, int N, int comp1, int comp2, int bs, int inveps, 
		  int *ind,  int **box, int *lis) {
  int d,i,ix,iy,ixy;
  int ib=bs-1;
  double **xx;

  xx=(double **)calloc(dim, sizeof(double*)); 
  for (d=0;d<dim;d++) {
    xx[d]=(double *)calloc(N, sizeof(double)); 
    memcpy(xx[d],x[d],N*sizeof(double));
  }
  for (ix=0;ix<bs;ix++) for (iy=0;iy<bs;iy++) box[ix][iy] = -1;
  
  for (i=0;i<N;i++) {
    ix=(int)(x[comp1][i]*inveps)&ib;iy=(int)(x[comp2][i]*inveps)&ib;
    lis[i]=box[ix][iy];box[ix][iy]=i;
  }
  i=-1;
  for (ix=0;ix<bs;ix++) for (iy=0;iy<bs;iy++) {
    ixy=box[ix][iy];
    while(ixy>=0) {
      i++;
      for (d=0;d<dim;d++) x[d][i]=xx[d][ixy];
      ind[ixy]=i;ixy=lis[ixy];
    }
    box[ix][iy]=-1;
  }
  for (i=0;i<N;i++) {
    ix=(int)(x[comp1][i]*inveps)&ib;iy=(int)(x[comp2][i]*inveps)&ib;
    lis[i]=box[ix][iy];box[ix][iy]=i;
  }
  for (d=0;d<dim;d++) free(xx[d]); free(xx);
}
int neiE1(double *x, int i, double scal, int bs, double eps,  int *box, int *lis, int *mxi) {
  double dd;
  int mi,mp,mm,nx=0;
  double xc=x[i];

  mp=int((xc+eps)*scal); if(mp>bs) mp=bs;
  mm=int((xc-eps)*scal); if(mm<0)  mm=0;
  mi=box[mp];
  while(mi>=0) {
    dd=x[mi]-xc;if(fabs(dd)<=eps) nx++;
    mi=lis[mi];
  }
  if(mm>=mp) return nx-1;
  mi=box[mm];
  while(mi>=0) {
    dd=xc-x[mi];if(fabs(dd)<=eps) nx++;
    mi=lis[mi];
  }
  nx+=mxi[mp-1]-mxi[mm];
  return nx-1;
}
int neiE(double **x, int i, int comp1, int comp2, int dim, int bs, double epsgrid, double eps, int **box, int *lis) {
  int ix,iy,ix1,iy1,ix2,jj,step,d;
  int el,nx,ib=bs-1;
  double dd,dy;
  double *xx;

  xx=(double *)calloc(dim,sizeof(double));
  for (d=0;d<dim;d++) xx[d]=x[d][i];

  ix=(int)(xx[comp1]/epsgrid)&ib; iy=(int)(xx[comp2]/epsgrid)&ib;
  jj=0; nx=0;
  while (eps>epsgrid*(jj-1)) {
    step = (jj) ? 2*jj : 1;
    for (ix1=ix-jj;ix1<=ix+jj;ix1++) {
      ix2=ix1&ib;
      for (iy1=iy-jj;iy1<=iy+jj;iy1+=step) {
	el=box[ix2][iy1&ib];
	while (el != -1) { 
	  dd=fabs(xx[0]-x[0][el]); 
	  for (d=1;d<dim;d++) if ( (dy=fabs(xx[d]-x[d][el]))>dd ) if ( (dd=dy) > eps ) break;
	  if (dd<=eps) nx++;
	  el=lis[el];
	}
      }
    }
    for (ix1=ix-jj;ix1<=ix+jj;ix1+=step) {
      ix2=ix1&ib;
      for (iy1=iy-jj+1;iy1<=iy+jj-1;iy1++) {
	el=box[ix2][iy1&ib];
	while (el != -1) { 
	  dd=fabs(xx[0]-x[0][el]); 
	  for (d=1;d<dim;d++) if ( (dy=fabs(xx[d]-x[d][el]))>dd ) if ( (dd=dy) > eps ) break;
	  if (dd<=eps) nx++;
	  el=lis[el];
	}
      }
    }
    jj++;
    if (jj==(bs/2)) break;
  }
  if ( jj==(bs/2) ) { //half of the layer
    for (ix1=ix-jj;ix1<ix+jj;ix1++) {
      ix2=ix1&ib; iy1=iy-jj;
      el=box[ix2][iy1&ib];
      while (el != -1) { 
	dd=fabs(xx[0]-x[0][el]); 
	for (d=1;d<dim;d++) if ( (dy=fabs(xx[d]-x[d][el]))>dd ) if ( (dd=dy) > eps ) break;
	if (dd<=eps) nx++;
	el=lis[el];
      }
    }
    ix1=ix-jj; ix2=ix1&ib;
    for (iy1=iy-jj+1;iy1<=iy+jj-1;iy1++) {
      el=box[ix2][iy1&ib];
      while (el != -1) { 
	dd=fabs(xx[0]-x[0][el]); 
	for (d=1;d<dim;d++) if ( (dy=fabs(xx[d]-x[d][el]))>dd ) if ( (dd=dy) > eps ) break;
	if (dd<=eps) nx++;
	el=lis[el];
      }
    }
  }
  free(xx);
  return nx-1;
}

void neiEK(double **x, int i, int comp1, int comp2, int dim, int K,
	   int bs, double epsgrid, double *eps, int **box, int *lis,
	   int *nx) {
  int ix,iy,ix1,iy1,ix2,jj,step,d,ik;
  int el,ib=bs-1;
  double dd,dy;
  double *xx;

  xx=(double *)calloc(dim,sizeof(double));
  for (d=0;d<dim;d++) xx[d]=x[d][i];

  ix=(int)(xx[comp1]/epsgrid)&ib; iy=(int)(xx[comp2]/epsgrid)&ib;
  jj=0; 
  for (ik=0;ik<K;ik++) nx[ik]=0;
  while (eps[K-1]>epsgrid*(jj-1)) {
    step = (jj) ? 2*jj : 1;
    for (ix1=ix-jj;ix1<=ix+jj;ix1++) {
      ix2=ix1&ib;
      for (iy1=iy-jj;iy1<=iy+jj;iy1+=step) {
	el=box[ix2][iy1&ib];
	while (el != -1) { 
	  dd=fabs(xx[0]-x[0][el]); 
	  for (d=1;d<dim;d++) if ( (dy=fabs(xx[d]-x[d][el]))>dd ) if ( (dd=dy) > eps[K-1] ) break;; 
	  for (ik=0;ik<K;ik++) {
	    if (dd<=eps[ik]) nx[ik]++;
	  }
	  el=lis[el];
	}
      }
    }
    for (ix1=ix-jj;ix1<=ix+jj;ix1+=step) {
      ix2=ix1&ib;
      for (iy1=iy-jj+1;iy1<=iy+jj-1;iy1++) {
	el=box[ix2][iy1&ib];
	while (el != -1) { 
	  dd=fabs(xx[0]-x[0][el]); 
	  for (d=1;d<dim;d++) if ( (dy=fabs(xx[d]-x[d][el]))>dd ) if ( (dd=dy) > eps[K-1] ) break;
	  for (ik=0;ik<K;ik++) {
	    if (dd<=eps[ik]) nx[ik]++;
	  }
	  el=lis[el];
	}
      }
    }
    jj++;
    if (jj==(bs/2)) break;
  }
  if ( jj==(bs/2) ) { //half of the layer
    for (ix1=ix-jj;ix1<ix+jj;ix1++) {
      ix2=ix1&ib; iy1=iy-jj;
      el=box[ix2][iy1&ib];
      while (el != -1) { 
	dd=fabs(xx[0]-x[0][el]); 
	for (d=1;d<dim;d++) if ( (dy=fabs(xx[d]-x[d][el]))>dd ) if ( (dd=dy) > eps[K-1] ) break;
	for (ik=0;ik<K;ik++) {
	  if (dd<=eps[ik]) nx[ik]++;
	}
	el=lis[el];
      }
    }
    ix1=ix-jj; ix2=ix1&ib;
    for (iy1=iy-jj+1;iy1<=iy+jj-1;iy1++) {
      el=box[ix2][iy1&ib];
      while (el != -1) { 
	dd=fabs(xx[0]-x[0][el]); 
	for (d=1;d<dim;d++) if ( (dy=fabs(xx[d]-x[d][el]))>dd ) if ( (dd=dy) > eps[K-1] ) break;
	for (ik=0;ik<K;ik++) {
	  if (dd<=eps[ik]) nx[ik]++;
	}
	el=lis[el];
      }
    }
  }
  free(xx);
  for (ik=0;ik<K;ik++) nx[ik]-=1;
}


void neiK(double **x, int dim, int comp1, int comp2, int i, 
	  int bs, double epsgrid, int K, int **box, int *lis,
	  int *nn) {
  double *dn,*xx;
  int k,ix,iy,ix1,iy1,ix2,jj,step,ib=bs-1;
  int el;
  double dd,dy;
  int d;

  dn=(double*)calloc(K+1,sizeof(double));
  xx=(double*)calloc(dim,sizeof(double));
  for(k=0;k<dim;k++) xx[k]=x[k][i];
  dn[0]=0;
  for(k=1;k<=K;k++) dn[k]=1e30;
  ix=(int)(xx[comp1]/epsgrid)&ib; iy=(int)(xx[comp2]/epsgrid)&ib;
  jj=0; 
  while (dn[K]>epsgrid*(jj-1)) {
    step = (jj) ? 2*jj : 1;
    for (ix1=ix-jj;ix1<=ix+jj;ix1++) {
      ix2=ix1&ib;
      for (iy1=iy-jj;iy1<=iy+jj;iy1+=step) {
	el=box[ix2][iy1&ib];
	while (el != -1) { if (el!=i) {
	  dd=fabs(xx[0]-x[0][el]); 
	  for (d=1;d<dim;d++) if ( (dy=fabs(xx[d]-x[d][el]))>dd ) dd=dy;
	  if(dd<dn[K]) {
	    k=K; while(dd<dn[k]) { if (k<K) { dn[k+1]=dn[k]; nn[k+1]=nn[k]; } k--; }
	    dn[k+1]=dd; nn[k+1]=el;
	  }
	} el=lis[el];
	}
      }
    }
    for (ix1=ix-jj;ix1<=ix+jj;ix1+=step) {
      ix2=ix1&ib;
      for (iy1=iy-jj+1;iy1<=iy+jj-1;iy1++) {
	el=box[ix2][iy1&ib];
	while (el != -1) { if (el!=i) {
	  dd=fabs(xx[0]-x[0][el]); 
	  for (d=1;d<dim;d++) if ( (dy=fabs(xx[d]-x[d][el]))>dd ) dd=dy;
	  if ( (dy=fabs(xx[1]-x[1][el]))>dd ) dd=dy;
	  if(dd<dn[K]) {
	    k=K; while(dd<dn[k]) { if (k<K) { dn[k+1]=dn[k]; nn[k+1]=nn[k]; } k--; }
	    dn[k+1]=dd; nn[k+1]=el;
	  } 
	} el=lis[el];
	}
      }
    }
    jj++;
    if (jj==bs/2) break;
  }

  if ( jj==(bs/2) ) { //half of the layer
    for (ix1=ix-jj;ix1<ix+jj;ix1++) {
      ix2=ix1&ib; iy1=iy-jj;
      el=box[ix2][iy1&ib];
      while (el != -1) { 
	if (el!=i) {
	  dd=fabs(xx[0]-x[0][el]); 
	  for (d=1;d<dim;d++) if ( (dy=fabs(xx[d]-x[d][el]))>dd ) dd=dy;
	  if(dd<dn[K]) {
	    k=K; while(dd<dn[k]) { if (k<K) { dn[k+1]=dn[k]; nn[k+1]=nn[k]; } k--; }
	    dn[k+1]=dd; nn[k+1]=el;
	  }
	}
	el=lis[el];
      }
    }
    ix1=ix-jj; ix2=ix1&ib;
    for (iy1=iy-jj+1;iy1<=iy+jj-1;iy1++) {
      el=box[ix2][iy1&ib];
      while (el != -1) { 
	if (el!=i) {
	  dd=fabs(xx[0]-x[0][el]); 
	  for (d=1;d<dim;d++) if ( (dy=fabs(xx[d]-x[d][el]))>dd ) dd=dy;
	  if(dd<dn[K]) {
	    k=K; while(dd<dn[k]) { if (k<K) { dn[k+1]=dn[k]; nn[k+1]=nn[k]; } k--; }
	    dn[k+1]=dd; nn[k+1]=el;
	  }
	} 
	el=lis[el];
      }
    }
  }
  free(dn);free(xx);
}


/*********/
void mi2(double **x, 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;

  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));
  n1=(int*)calloc(dim,sizeof(int));
  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]);      
  
  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]]; }
      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 %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;
  
  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 mi2c(double **x, 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;

  int dim=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;

  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,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=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 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