test.cpp

最近点对问题

#include <iostream>
#include <fstream>
#include <string>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <cmath>
#include <cassert>
#include "RTree.h"
using namespace std;
#define MAXINT 32767
#define sqr(a) ((a)*(a))
#define DIM  10
int NUMBER=0;
int Statistics_Count;
double* Statistics_SumX;
double* Statistics_SumXsquared;
double* Statistics_SumProduct;
void InitStatistics(int Dimensions)
//   ==============
// initializes variables for the statistics
{
  Statistics_SumX = new double[Dimensions];
  Statistics_SumXsquared = new double[Dimensions];
  Statistics_SumProduct = new double[Dimensions*Dimensions];
  
  Statistics_Count = 0;
  for (int d=0; d<Dimensions; d++) {
    Statistics_SumX[d]=0.0;
    Statistics_SumXsquared[d]=0.0;
    for (int dd=0; dd<Dimensions; dd++) Statistics_SumProduct[d*Dimensions+dd] = 0.0;
  }
}


void EnterStatistics(int Dimensions,double* x)
//   ===============
// counts the vector "x" in the statistics
{
  Statistics_Count++;
  for (int d=0; d<Dimensions; d++) {
    Statistics_SumX[d] += x[d];
    Statistics_SumXsquared[d] += x[d]*x[d];
    for (int dd=0; dd<Dimensions; dd++) Statistics_SumProduct[d*Dimensions+dd] += x[d]*x[dd];
  }
}  


void OutputStatistics(int Dimensions)
//   ================
// outputs the statistics
{
  for (int d=0; d<Dimensions; d++) {
    double E = Statistics_SumX[d] / Statistics_Count;
    double V = Statistics_SumXsquared[d]/Statistics_Count - E*E;
    double s = sqrt(V);
    
  }
  
  for ( d=0; d<Dimensions; d++) {
    for (int dd=0; dd<Dimensions; dd++) {
      double Kov = (Statistics_SumProduct[d*Dimensions+dd]/Statistics_Count) -
                   (Statistics_SumX[d]/Statistics_Count) * (Statistics_SumX[dd]/Statistics_Count);
      double Cor = Kov /
                   sqrt(Statistics_SumXsquared[d]/Statistics_Count - sqr(Statistics_SumX[d] / Statistics_Count)) /
                   sqrt(Statistics_SumXsquared[dd]/Statistics_Count - sqr(Statistics_SumX[dd] / Statistics_Count));
      
    }
   
  }
  
}


double RandomEqual(double min,double max)
//     ===========
// generates a random number equally distributed in the interval [min,max[
{
  // double x = (double)rand()/MAXINT;   <-- patch
  double x = (double)rand()/RAND_MAX;
  return x*(max-min)+min;
}


double RandomPeak(double min,double max,int dim)
//     ==========
// generates a random number in the interval [min,max[
// as a sum of "dim" equally distributes random numbers
{
  double sum = 0.0;
  for (int d=0; d<dim; d++) sum += RandomEqual(0,1);
  sum /= dim;
  return sum*(max-min)+min;
}


double RandomNormal(double med,double var)
//     ============
// generates a normally distributed random number with mean "med"
// in the interval ]med-var,med+var[
{
  return RandomPeak(med-var,med+var,12);
}


void GenerateDataEqually(FILE* f,int Count,int Dimensions)
//   ===================
// generates "Count" independently equally distributed vectors
// in the file "f"
{
  InitStatistics(Dimensions);
  for (int i=0; i<Count; i++) {
    double x[DIM];
    for (int d=0; d<Dimensions; d++) {
      x[d] = RandomEqual(0,1);
      fprintf(f,"%8.6f ",x[d]);
    }
    EnterStatistics(Dimensions,x);
    fprintf(f,"\n");
  }
  OutputStatistics(Dimensions);
}


void GenerateDataCorrelated(FILE* f,int Count,int Dimensions)
//   ======================
// generates "Count" correlated vectors in the file "f"
{
  InitStatistics(Dimensions);
  double x[DIM];
  for (int i=0; i<Count; i++) {
again:
    double v = RandomPeak(0,1,Dimensions);
    for (int d=0; d<Dimensions; d++) x[d] = v;
    double l = v<=0.5 ? v:1.0-v;
    for ( d=0; d<Dimensions; d++) {
      double h = RandomNormal(0,l);
      x[d] += h;
      x[(d+1)%Dimensions] -= h;
    }
    for ( d=0; d<Dimensions; d++) if (x[d]<0 || x[d]>=1) goto again;
    for ( d=0; d<Dimensions; d++) fprintf(f,"%8.6f ",x[d]);
    EnterStatistics(Dimensions,x);
    fprintf(f,"\n");
  }
  OutputStatistics(Dimensions);
}


void GenerateDataAnticorrelated(FILE* f,int Count,int Dimensions)
//   ==========================
// generates "Count" anti-correlated vectors in the file "f"
{
  InitStatistics(Dimensions);
  double x[DIM];
  for (int i=0; i<Count; i++) {
again:
    double v = RandomNormal(0.5,0.25);
    for (int d=0; d<Dimensions; d++) x[d] = v;
    double l = v<=0.5 ? v:1.0-v;
    for ( d=0; d<Dimensions; d++) {
      double h = RandomEqual(-l,l);
      x[d] += h;
      x[(d+1)%Dimensions] -= h;
    }
    for (d=0; d<Dimensions; d++) if (x[d]<0 || x[d]>=1) goto again;
    for ( d=0; d<Dimensions; d++) fprintf(f,"%8.6f ",x[d]);
    EnterStatistics(Dimensions,x);
    fprintf(f,"\n");
  }
  OutputStatistics(Dimensions);
}


void GenerateData(int Dimensions,char Distribution,int Count,char* FileName)
//   ============
// generates a file with random data
{
  if (Count <= 0) {
    printf("Ungtige Anzahl von Punkten.\n");
    return;
  }
  if (Dimensions < 2) {
    printf("Ungtige Anzahl von Dimensionen.\n");
    return;
  }
  switch (Distribution) {
    case 'E':
    case 'e': Distribution = 'E'; break;
    case 'C':
    case 'c': Distribution = 'C'; break;
    case 'A':
    case 'a': Distribution = 'A'; break;
    default: printf("Invalid distribution.\n"); return;
  }
  
  FILE* f = fopen(FileName,"wt");
  if (f == NULL) {
    printf("Cannot create file \"%s\".\n",FileName);
    return;
  }
  fprintf(f,"%d %d\n",Count,Dimensions);
  switch (Distribution) {
    case 'E': GenerateDataEqually(f,Count,Dimensions); break;
    case 'C': GenerateDataCorrelated(f,Count,Dimensions); break;
    case 'A': GenerateDataAnticorrelated(f,Count,Dimensions); break;
  }
  fclose(f);
  printf("%d vectors generated, filename \"%s\".\n",Count,FileName);
}

class NNTree:public RTree<long,int,2>{
public:
	          int MINDIST(Branch *branch,int[]);		
	          int MINMAXDIST(Branch *branch,int[]);	
	          void Push(Node* node,int[]);	
	          void NNSEARCH(int[]);	
     private:
	          struct Heap		 
			  {
		           Branch* branch;
		           int level;
			  };
	          int E[10000][2];	
	          Heap heap[1000];		
};
struct info
{
	int MinDist;
	int no;
	
};
struct info nn[100];
int NNTree::MINDIST(Branch* branch,int lookup[])
{
	int r,distant=0;
	for(int i=0;i<2;i++)
	{
		
		if(lookup[i]<branch->m_rect.m_min[i])
		    r=branch->m_rect.m_min[i];
		else if(lookup[i]>branch->m_rect.m_max[i])
			r=branch->m_rect.m_max[i];
		else
			r=lookup[i];
        distant=distant+(lookup[i]-r)*(lookup[i]-r);
	}
    return distant;
}
int NNTree::MINMAXDIST(Branch* branch,int lookup[])
{
	int r,distant=0,dis[2];
	for(int i=0;i<2;i++)
	{
		distant=0;		
		for(int j=0;j<2;j++)
		{

			if(j==i)
			{
				if(lookup[j]<=(branch->m_rect.m_min[j]+branch->m_rect.m_max[j])/2)
					r=branch->m_rect.m_min[j];
				else r=branch->m_rect.m_max[j];
			}
			else if(j!=i)
			{
				if(lookup[j]>=(branch->m_rect.m_min[j]+branch->m_rect.m_max[j])/2)
					r=branch->m_rect.m_min[j];
				else r=branch->m_rect.m_max[j];
			}
			distant=distant+(lookup[j]-r)*(lookup[j]-r);		
		}
        dis[i]=distant;	
	}
    if(dis[0]<=dis[1])
	distant=dis[0];
	else distant=dis[1];
	return distant;
}
void NNTree::Push(Node* node,int lookup[])
{
	
	int i,j=0,min_dist,node_no,k,num,change;
    for(i=0;i<MAXNODES;i++)
	{
		if((node->m_branch[i].m_rect.m_min[0]<0)&&(node->m_branch[i].m_rect.m_max[0]<0)) 
		continue;
		nn[j].MinDist=MINDIST(node->m_branch+i,lookup);
		nn[j].no=i;//结点标号
		j++;
	}
	num=j;//数组nn中结点数
	for(i=num-1,change=1;i>=1&&change;--i)
	{
		change=0;
		for(j=0;j<i;++j)
		{
			if(nn[j].MinDist>nn[j+1].MinDist)
			{
               min_dist=nn[j].MinDist;
			   node_no=nn[j].no;
			   nn[j].MinDist=nn[j+1].MinDist;
			   nn[j].no=nn[j+1].no;
			   nn[j+1].MinDist=min_dist;
			   nn[j+1].no=node_no;
			   change=1;
			}
		}
	}
	for(k=0;k<=num;k++)
	{
		i=nn[k].no;
		heap[NUMBER].branch=node->m_branch+i;
		heap[NUMBER].level=node->m_level;
		NUMBER++;
		
	}
}


void NNTree::NNSEARCH(int lookup[])
{
	int i=0,j=0,k=0,MinMaxDist,dist,Best=999999999,Edistant[10000];
	Push(m_root,lookup);//根结点压栈
    while(NUMBER>0)
	{
		NUMBER--;
		if((MINDIST(heap[NUMBER].branch,lookup))>Best)
		continue;
	    if((MinMaxDist=MINMAXDIST(heap[NUMBER].branch,lookup))<Best)
		Best=MinMaxDist;
        if(heap[NUMBER].level==0)  //是叶子结点
		{
			for(j=0;j<2;j++)
			{
				E[i][j]=heap[NUMBER].branch->m_rect.m_min[j];
			}
			i++;
		}		
		else{//是叶子结点
		Push(heap[NUMBER].branch->m_child,lookup);
		}
        for(k=0;k<i;k++)
		{
		     Edistant[k]=0;
             for(j=0;j<2;j++)
			 {
			 Edistant[k]+=(E[k][j]-lookup[j])*(E[k][j]-lookup[j]);
			 }
        Edistant[k]=sqrt(Edistant[k]);
		}
	}
    dist=Edistant[0];
	for(k=0;k<i;k++)
	{
		if(dist>Edistant[k])
		{
			dist=Edistant[k];
			j=k;
		}	
	}
	cout<<"离所要查询的点的最近的点是:";
	for(i=0;i<2;i++)
	{
		cout<<E[j][i]<<"  ";
	}
    cout<<endl;
    
}
int main()
//  ====
// main program
{  int a,lookup[2],dot[2],j=0;
   float number;
  	ifstream infile;
	NNTree nntree;
   cout<<"请输入需要的数据量:\n";
   cin>>a;
	GenerateData(2,'A',a,"a.txt");
    infile.open("a.txt");
	infile>>a;
	while(!infile.eof())
	{
		for(int i=0;i<2;i++)
		{
			infile>>number;
			dot[i]=number*1000;
		}
		nntree.Insert(dot,dot,j);
		j++;
	}
	cout<<"请输入要查找的点的坐标：\n";
   for(a=0;a<2;a++)
   {
	   cin>>lookup[a];
   }
    nntree.NNSEARCH(lookup);
  return 0;
}