www.pudn.com > MutualInformationICA.zip > miutils.C
#include
#include
#include
#include
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=0) {
i++;
for (d=0;d=0) {
i++;
for (d=0;dbs) 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;depsgrid*(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;ddd ) 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;ddd ) 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;ix1dd ) 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;ddd ) 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;depsgrid*(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;ddd ) if ( (dd=dy) > eps[K-1] ) break;;
for (ik=0;ikdd ) if ( (dd=dy) > eps[K-1] ) break;
for (ik=0;ikdd ) if ( (dd=dy) > eps[K-1] ) break;
for (ik=0;ikdd ) if ( (dd=dy) > eps[K-1] ) break;
for (ik=0;ikepsgrid*(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;ddd ) dd=dy;
if(dddd ) dd=dy;
if ( (dy=fabs(xx[1]-x[1][el]))>dd ) dd=dy;
if(dddd ) dd=dy;
if(dddd ) dd=dy;
if(ddeps[d] ) eps[d]=dx; }
if (eps[d]>Eps) {Eps=eps[d];maxdim=d;}
}
for (d=0;deps[d] ) eps[d]=dx; }
if (eps[d]>Eps) {Eps=eps[d];maxdim=d;}
}
for (d=0;deps[d] ) eps[d]=dx; }
}
for (d=0;deps[d] ) eps[d]=dx; }
if (eps[d]>Eps) {Eps=eps[d];maxdim=d;}
}
for (d=0;deps[d] ) eps[d]=dx; }
if (eps[d]>Eps) {Eps=eps[d];maxdim=d;}
}
for (d=0;deps[d] ) eps[d]=dx; }
}
for (d=0;deps[d] ) eps[d]=dx;
}
for (d=0;d1) {
xx=(double**)calloc(dimx,sizeof(double*));
for (d=0;d1) {
yy=(double**)calloc(dimy,sizeof(double*));
for (d=0;d1) for (d=0;d1) for (d=0;depsx ) epsx=dx;
epsy=0; for (d=dimx;depsy ) epsy=dy;
if (dimx>1) nx2=neiE(xx,indx[i],0,dimx-1,dimx,BOX,epsilon,epsx,boxx,lisx);
else nx2=neiE1(x[0],ind[i],scalx,BOX1,epsx,boxx1,lisx1,mxi);
if (dimy>1) ny2=neiE(yy,indy[i],0,dimy-1,dimy,BOX,epsilon,epsy,boxy,lisy);
else ny2=neiE1(x[dimx],ind[i],scaly,BOX1,epsy,boxy1,lisy1,myi);
if (epsx>epsy) {
Eps=epsx;nx1=nx2;
if (dimy>1) ny1=neiE(yy,indy[i],0,dimy-1,dimy,BOX,epsilon,Eps,boxy,lisy);
else ny1=neiE1(x[dimx],ind[i],scaly,BOX1,Eps,boxy1,lisy1,myi);
dxy1+=psi[nx1]+psi[ny1+1];
} else {
Eps=epsy;ny1=ny2;
if (dimx>1) nx1=neiE(xx,indx[i],0,dimx-1,dimx,BOX,epsilon,Eps,boxx,lisx);
else nx1=neiE1(x[0],ind[i],scalx,BOX1,Eps,boxx1,lisx1,mxi);
dxy1+=psi[nx1+1]+psi[ny1];
}
dxy2+=psi[nx2]+psi[ny2];
}
dxy1/=N;*mic=psi[N]+psi[K]-dxy1;
dxy2/=N;*mir=psi[N]+phi[K]-dxy2;
free(xc);free(nn);free(lis);
for (i=0;i1) {
xx=(double**)calloc(dimx,sizeof(double*));
for (d=0;d1) {
yy=(double**)calloc(dimy,sizeof(double*));
for (d=0;d1) for (d=0;d1) for (d=0;depsx ) epsx=dx;
epsy=0; for (d=dimx;depsy ) epsy=dy;
if (epsx>epsy) { Eps=epsx;} else {Eps=epsy;}
if (dimy>1) ny1=neiE(yy,indy[i],0,dimy-1,dimy,BOX,epsilon,Eps,boxy,lisy);
else ny1=neiE1(x[dimx],ind[i],scaly,BOX1,Eps,boxy1,lisy1,myi);
if (dimx>1) nx1=neiE(xx,indx[i],0,dimx-1,dimx,BOX,epsilon,Eps,boxx,lisx);
else nx1=neiE1(x[0],ind[i],scalx,BOX1,Eps,boxx1,lisx1,mxi);
if (epsx>epsy) {
dxy1+=psi[nx1]+psi[ny1+1];
} else {
dxy1+=psi[nx1+1]+psi[ny1];
}
}
dxy1/=N;*mic=psi[N]+psi[K]-dxy1;
free(xc);free(nn);free(lis);
for (i=0;i1) {
xx=(double**)calloc(dimx,sizeof(double*));
for (d=0;d1) {
yy=(double**)calloc(dimy,sizeof(double*));
for (d=0;d1) for (d=0;d1) for (d=0;depsx ) epsx=dx;
epsy=0; for (d=dimx;depsy ) epsy=dy;
if (dimx>1) nx2=neiE(xx,indx[i],0,dimx-1,dimx,BOX,epsilon,epsx,boxx,lisx);
else nx2=neiE1(x[0],ind[i],scalx,BOX1,epsx,boxx1,lisx1,mxi);
if (dimy>1) ny2=neiE(yy,indy[i],0,dimy-1,dimy,BOX,epsilon,epsy,boxy,lisy);
else ny2=neiE1(x[dimx],ind[i],scaly,BOX1,epsy,boxy1,lisy1,myi);
dxy2+=psi[nx2]+psi[ny2];
}
dxy2/=N;*mir=psi[N]+phi[K]-dxy2;
free(xc);free(nn);free(lis);
for (i=0;ieps[ik][d] ) eps[ik][d]=dx; }
if (eps[ik][d]>Eps[ik]) {Eps[ik]=eps[ik][d];maxdim[ik]=d;}
}
}
for (ik=0;ikeps[ik][d] ) eps[ik][d]=dx; }
if (eps[ik][d]>Eps[ik]) {Eps[ik]=eps[ik][d];maxdim[ik]=d;}
}
}
for (ik=0;ik1) {
xx=(double**)calloc(dimx,sizeof(double*));
for (d=0;d1) {
yy=(double**)calloc(dimy,sizeof(double*));
for (d=0;d1) for (d=0;d1) for (d=0;depsx[ik] ) epsx[ik]=dx;
epsy[ik]=0; for (d=dimx;depsy[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;ikepsx[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;ikepsy[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;ikeps[d] ) eps[d]=dx; }
if (eps[d]>Eps) {Eps=eps[d];maxdim=d;}
}
for (d=0;deps[ik][d] ) eps[ik][d]=dx; }
if (eps[ik][d]>Eps[ik]) {Eps[ik]=eps[ik][d];maxdim[ik]=d;}
}
}
for (ik=0;ikeps[d] ) eps[d]=dx; }
}
for (d=0;d