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