OpErgodicity_Compare.cpp File Reference

Functions
double	HistogramIntersection (double Dist1[Nbin], Integer_t Dist1Ndata[Nbin], double Dist2[Nbin])
double	KLDivergence (double Dist1[Nbin], Integer_t Dist1Ndata[Nbin], double Dist2[Nbin])
double	JSDivergence (double Dist1[Nbin], Integer_t Dist1Ndata[Nbin], double Dist2[Nbin], Integer_t Dist2Ndata[Nbin])
double	L1Norm (double Dist1[Nbin], Integer_t Dist1Ndata[Nbin], double Dist2[Nbin])
double	L2Norm (double Dist1[Nbin], Integer_t Dist1Ndata[Nbin], double Dist2[Nbin])
int	main (int argc, char **argv)

Variables
Integer_t const	Nbin = 100
Integer_t const	SistemSizeMAX = 20
Integer_t const	Nitem = 4

Function Documentation

HistogramIntersection()

double HistogramIntersection	(	double	Dist1[Nbin],
		Integer_t	Dist1Ndata[Nbin],
		double	Dist2[Nbin]
	)

                                                                                                 {
  Integer_t i;
  double count = 0;
  double res = 0;
  for(i = 0;i < Nbin; ++i) {
    if( Dist1Ndata[i] == 0 ) continue;
    if( Dist1[i] > Dist2[i] ) res += Dist2[i];
    else res += Dist1[i];
    count += Dist1[i];
    // if( res ) printf("KLDivergence: Dist1[%lld]=%f, Dist2[%lld]=%f\n", i, Dist1[i], i, Dist2[i]);
  }
  return res/count;
}

JSDivergence()

double JSDivergence	(	double	Dist1[Nbin],
		Integer_t	Dist1Ndata[Nbin],
		double	Dist2[Nbin],
		Integer_t	Dist2Ndata[Nbin]
	)

                                                                                                                    {
  double res = KLDivergence(Dist1, Dist1Ndata, Dist2);
  res += KLDivergence(Dist2, Dist2Ndata, Dist1);
  return res /= 2.0;
}

KLDivergence()

double KLDivergence	(	double	Dist1[Nbin],
		Integer_t	Dist1Ndata[Nbin],
		double	Dist2[Nbin]
	)

                                                                                        {
  Integer_t i;
  double count = 0;
  double res = 0;
  for(i = 0;i < Nbin; ++i) {
    if( Dist1[i] == 0 || Dist2[i] == 0 ) continue;
    res   += Dist1[i]*log(Dist1[i]/Dist2[i]);
    // if( isnan(res) ) printf("KLDivergence: Dist1[%lld]=%f, Dist2[%lld]=%f\n", i, Dist1[i], i, Dist2[i]);
    count += Dist1[i];
  }
  return res/(double)count;
}

L1Norm()

double L1Norm	(	double	Dist1[Nbin],
		Integer_t	Dist1Ndata[Nbin],
		double	Dist2[Nbin]
	)

                                                                                  {
  Integer_t i;
  double count = 0;
  double res = 0;
  for(i = 0;i < Nbin; ++i) {
    if( Dist1Ndata[i] == 0 ) continue;
    res += fabs(Dist1[i] -Dist2[i]);
    // if( isnan(res) ) printf("L1Norm: Dist1[%lld]=%f, Dist2[%lld]=%f\n", i, Dist1[i], i, Dist2[i]);
    count += 1.0;
  }
  // printf("L1Norm: res=%f, count=%f\n", res, count);
  return res/count;
}

L2Norm()

double L2Norm	(	double	Dist1[Nbin],
		Integer_t	Dist1Ndata[Nbin],
		double	Dist2[Nbin]
	)

                                                                                  {
  Integer_t i;
  double count = 0;
  double res = 0;
  for(i = 0;i < Nbin; ++i) {
    if( Dist1Ndata[i] == 0 ) continue;
    res += (Dist1[i] -Dist2[i])*(Dist1[i] -Dist2[i]);
    // if( isnan(res) ) printf("L2Norm: Dist1[%lld]=%f, Dist2[%lld]=%f\n", i, Dist1[i], i, Dist2[i]);
    count += 1.0;
  }
  // printf("L2Norm: res=%f, count=%f\n", res, count);
  return sqrt(res/count);
}

main()

int main	(	int	argc,
		char **	argv
	)

                               {
  if(argc < 2) {
    printf("Usage: 1.This 2.Directory\n");
    exit(EXIT_FAILURE);
  }
  Integer_t i,j;
  Integer_t Nsample = 0;
  Integer_t n, bin, dimShell, itemp;
  Integer_t NnanfE[SistemSizeMAX];
  Integer_t Ntotal[SistemSizeMAX];
  double ratio, expval, stddev, fE, dos, MCave, dtemp;
 
  struct stat statBuf;
  struct Distance {
    double HI;
    double KLdivergence;
    double JSdivergence;
    double L1norm;
    double L2norm;
  };
  struct Distance distance[SistemSizeMAX][Nitem]; // 分布間の距離の指標
 
  struct EnsembleAverage {
    Integer_t N[Nbin];
    double ExpValue[Nbin];
    double Stddev[Nbin];
    double funcF[Nbin];
    double DOS[Nbin];
  };
  struct EnsembleAverage Average[SistemSizeMAX];
  for(i = 0; i < SistemSizeMAX; ++i) for(j = 0;j < Nbin; ++j) {
    Average[i].N[j] = 0;
    Average[i].ExpValue[j] = 0;
    Average[i].Stddev[j]   = 0;
    Average[i].funcF[j]    = 0;
    Average[i].DOS[j] = 0;
  }
 
  Integer_t const BUFSIZE = 1000;
  char dstr[BUFSIZE], str[BUFSIZE], stemp[BUFSIZE];
  char *filename = str;
  char *ext;
  strcpy(dstr, argv[1]);
  DIR *dir;
  FILE *Input, *Output;
    struct dirent *dp;
 
  //********** アンサンブル平均の読み込み ********** //
  strcpy(stemp, dstr);
  strcat(stemp, "/Output");
  strcpy(str, stemp);
  strcat(str, "/");
  strcat(str, "EnsembleAverage.txt");
  printf("Input file: %s\n", filename);
  Input = fopen(filename, "r");
  if(Input == NULL) {
    printf("Can't open a file %s\n", filename);
    exit(EXIT_FAILURE);
  }
  skip_header(Input);
  while( fscanf(Input, "%lld %lf %le %le %le %le %lld", &n, &ratio, &expval, &stddev, &fE, &dos, &itemp) != EOF ) {
    if(ratio == 1) bin = Nbin-1; else bin = (Integer_t)(ratio*Nbin);
    Average[n].N[bin] = itemp;
    Average[n].ExpValue[bin] = expval;
    Average[n].Stddev[bin]   = stddev;
    Average[n].funcF[bin]    = fE;
    Average[n].DOS[bin] = dos;
  }
  fclose(Input);
  //********** (END)アンサンブル平均の読み込み ********** //
 
  //********** Output fileを作成 ********** //
  strcpy(stemp, dstr);
  strcat(stemp, "/Output");
  mkdir(stemp, 0777);
  strcpy(str, stemp);
  strcat(str, "/");
  strcat(str, "ErgodicityFluctuation_Compare.txt");
  strcpy(stemp, str);
  printf("Output file: %s\n", filename);
  Output = fopen(filename, "w");
  if(Output == NULL) {
    printf("Can't open a file %s\n", filename);
    exit(EXIT_FAILURE);
  }
  fprintf(Output, "# 1.Nsample 2.N 3.funcF L1norm 4.funcF L2norm\n\n");
  //********** (END)Output fileを作成 ********** //
 
  dir = opendir(dstr);
  for(dp = readdir(dir);dp != NULL;dp = readdir(dir)) {
    if( strstr(dp->d_name, "Operators_") == NULL ) continue;
    ext = strrchr(dp->d_name, '.');
    if( strcmp(ext, ".txt") != 0 ) continue;
 
    strcpy(str, dstr);
    strcat(str, "/");
    strcat(str, dp->d_name);
    if( stat(filename, &statBuf) == -1 ){
      printf("Can't get file information: %s.", filename);
      continue;
    }
    printf("Input file: %s, size: %lld\n", filename, statBuf.st_size);
    if( statBuf.st_size==0 ) continue;
    Input = fopen(filename, "r");
    if(Input == NULL) {
      printf("Can't open a file %s\n", filename);
      continue;
    }
    //********** 変数の初期化 **********//
      for(n = 0;n < SistemSizeMAX; ++n) {
        Ntotal[n] = 0;
        NnanfE[n] = 0;
      }
      for(n = 0;n < SistemSizeMAX; ++n) for(j = 0;j < Nitem; ++j) {
        distance[n][j].HI = 0;
        distance[n][j].KLdivergence = 0;
        distance[n][j].JSdivergence = 0;
        distance[n][j].L1norm = 0;
        distance[n][j].L2norm = 0;
      }
    //********** (END)変数の初期化 **********//
 
    Nsample += 1;
    skip_header(Input);
    //********** 関数の間の距離を計算 **********//
    while( fscanf(Input, "%lld %lld %le %le %le %le %le %le %lld", &itemp, &n, &dtemp, &ratio, &expval, &MCave, &fE, &dtemp, &dimShell) != EOF ) {
      if( isnan(fE) ) {
        NnanfE[n] += 1;
        continue;
      }
      Ntotal[n] += 1;
      if(ratio == 1) bin = Nbin-1; else bin = (Integer_t)(ratio*Nbin);
      // funcF
      distance[n][2].L1norm += fabs(fE -Average[n].funcF[bin]);
      distance[n][2].L2norm += (fE -Average[n].funcF[bin])*(fE -Average[n].funcF[bin]);
    }
    //********** (END)関数の間の距離を計算 **********//
    fclose(Input);
 
    //********** ファイルへの書き込み **********//
    for(n = 0; n < SistemSizeMAX; ++n) {
      if( Ntotal[n] == 0) continue;
      fprintf(Output, "%lld %lld ", Nsample, n);
      fprintf(Output, "%le ", distance[n][2].L1norm/(double)( Ntotal[n] ) );
      fprintf(Output, "%le ", sqrt( distance[n][2].L2norm/(double)( Ntotal[n] ) ) );
      fprintf(Output, "\n");
      printf("Ntotal[%lld]=%lld\n", n, Ntotal[n]);
    }
    //********** (END)ファイルへの書き込み **********//
 
    }
 
  printf("Output file: %s\n", stemp);
  fclose(Output);
 
  return 0;
}

Variable Documentation

Nbin

Integer_t const Nbin = 100

Nitem

Integer_t const Nitem = 4

SistemSizeMAX

Integer_t const SistemSizeMAX = 20