2 #include <Matematica.h> 5 #include <gsl/gsl_randist.h> 7 return gsl_rng_uniform(rng);
13 return Media + gsl_ran_gaussian(rng,Scarto);
24 double Rand = genrand_real1() ;
25 double Length = -log( Rand ) ;
26 double Theta = 2. * M_PI * genrand_real2() ;
27 double kdX = sqrt( 2. * Length ) * cos( Theta ) ;
28 return Scarto * kdX + Media ;
29 return RandomGaussian(Media,Scarto);
32 return genrand_real2();
33 return RandomUniform();
34 int IA=16807,IM=2147483647,IQ=127773,IR=2836,NDIV=1+(IM-1)/NTAB,kcas,jcas;
41 static int seme = 2351653;
44 seme = -seme>1 ? -seme : 1;
45 for(jcas=NTAB+8;jcas>=1;jcas--){
47 seme=IA*(seme-kcas*IQ)-IR*kcas;
48 if(seme<=0) seme += IM;
49 if(jcas<=NTAB) iv[jcas]=seme;
54 seme=IA*(seme-kcas*IQ)-IR*kcas;
59 double casuale=AM*iy < RNMX ? AM*iy : RNMX;
66 double InvNBin = 1./(double)NBin;
68 for(
int i=0;i<NBin-1;i++){
69 if(Ran > Prob[i] && Ran <= Prob[i+1]){
74 double xi = Ran - floor(Ran*NBin)*InvNBin;
75 double m = (Prob[j+1] - Prob[j])*NBin;
76 double y = j*InvNBin + m*xi;
82 fftw_complex *uscita = (fftw_complex *) fftw_malloc(
sizeof(fftw_complex)*(NMax) );
83 fftw_complex *entrata = (fftw_complex *) fftw_malloc(
sizeof(fftw_complex)*(NMax) );
85 for(
int i=0;i<NMax;i++){
86 entrata[0][i] = st[i];
90 p = fftw_plan_dft_1d(NMax,entrata,uscita,FFTW_FORWARD,FFTW_ESTIMATE);
92 for(
int i=0;i<NMax;i++){
93 sw[i] = QUAD(uscita[i][0])+QUAD(uscita[i][1]);
102 fftw_complex *uscita = (fftw_complex *) fftw_malloc(
sizeof(fftw_complex)*(NMax*NMax) );
103 fftw_plan p = fftw_plan_dft_r2c_2d(NMax,NMax,st,uscita,FFTW_MEASURE);
105 for(
int i=0;i<NMax;i++){
107 for(
int j=0;j<NMax;j++){
108 sw[i*NMax+j] = SQR(uscita[i*NMax+j][0])+SQR(uscita[i*NMax+j][0]);
111 fftw_destroy_plan(p);
116 double NMaxInv = 1./(double)NMax;
117 int NHalf = (int)(NMax/2.);
120 fftw_complex *uscita = (fftw_complex *)calloc(NMax*NMax,
sizeof(fftw_complex));
122 p = fftw_plan_dft_r2c_2d(NMax,NMax,st,uscita,FFTW_MEASURE);
124 for(
int i=0;i<NMax;i++){
125 for(
int j=0;j<NMax;j++){
126 sw[i][j] = QUAD(uscita[0][i*NMax+j])+QUAD(uscita[1][i*NMax+j]);
127 sw[i][j] *= QUAD(NMaxInv);
131 fftw_destroy_plan(p);
137 double Re1=0.,Im1=0.,Re2=0.,Im2=0.;
138 double FMass=(double)NMax/DUE_PI;
139 double dNMax = 1./(double)NMax;
140 for(
int j=0;j<NMax;j++){
141 Re1=0.;Re2=0.;Im1=0.;Im2=0.;
142 for(
int i=0;i<NMax;i++){
143 Re1+=st[i]*cos(DUE_PI*dNMax*(i)*j);
144 Im1-=st[i]*sin(DUE_PI*dNMax*(i)*j);
146 sw[j]= sqrt(QUAD((Re1+Re2))*dNMax+QUAD((Im1+Im2))*dNMax);
158 double dNMax = 1./(double)NMax;
159 int NHalf = (int)(NMax/2.);
160 for(
int kx=-NMax/2;kx<NMax/2;kx++){
161 double qx = kx*dNMax;
162 for(
int ky=-NMax/2;ky<NMax/2;ky++){
163 double qy = ky*dNMax;
164 double Re2=0.,Im2=0.;
165 double Re1=0.,Im1=0.;
166 for(
int lx=0;lx<NMax;lx++){
167 for(
int ly=0;ly<NMax;ly++){
168 double Arg = 2.*M_PI*(lx*kx + ly*ky);
169 double cy = cos(Arg);
170 double sy = sin(Arg);
171 Re1 += st[lx*NMax+ly]*cy;
172 Im1 += st[lx*NMax+ly]*sy;
175 int kkx = kx + NMax/2;
176 int kky = ky + NMax/2;
177 sw[kkx*NMax+kky] = SQR(Re1*dNMax) + SQR(Im1*dNMax);
206 double dNMax = 1./(double)NMax;
207 int NHalf = (int)(NMax/2.);
208 for(
int l=0;l<NMax;l++){
209 for(
int k=0;k<NMax;k++){
210 double Re1=0.,Im1=0.,Re2=0.,Im2=0.;
211 for(
int n=0;n<NMax;n++){
212 double Cosn = cos(DUE_PI*l*n*dNMax);
213 double Sinn = sin(DUE_PI*l*n*dNMax);
214 for(
int m=0;m<NMax;m++){
215 double Arg = DUE_PI*((l-NHalf)*n*dNMax + (k-NHalf)*m*dNMax);
218 Re1 += Cosn*cos(DUE_PI*k*m*dNMax)*st[n*NMax + m];
219 Re2 -= Sinn*sin(DUE_PI*k*m*dNMax)*st[n*NMax + m];
220 Im1 += Cosn*sin(DUE_PI*k*m*dNMax)*st[n*NMax + m];
221 Im2 += Sinn*cos(DUE_PI*k*m*dNMax)*st[n*NMax + m];
224 sw[l][k] = QUAD((Re1+Re2)*dNMax*dNMax) + QUAD((Im1+Im2)*dNMax*dNMax);
239 for(
int j=0;j<N;j++){
243 sw[j] = sqrt(-st[j]);
247 for(
int i=0;i<NMax;i++){
249 for(
int j=0,k=0;j<NMax;j++){
253 sAutocor[i] += (double)(st[k]*st[j]);
255 sAutocor[i] /= (double)(NMax);
262 double PartiInv = 1./(double)Parti;
263 int Nw = (int) (NMax/(
double)Parti);
264 for(
int i=0;i<Nw;i++)
266 for(
int i=0,j=0;i<NMax;i++){
278 int Nw = (int) (NMax/(
double)NParti);
279 double InvNParti = 1./(double)NParti;
280 for(
int i=0;i<Nw;i++){
284 for(
int i=0,j=0;i<NMax;i++){
286 sErr[j] += SQR((st[i]));
287 if( ((i+1)%NParti)==0 ){
289 sErr[j] = sqrt( (sErr[j] - SQR(sw[j])*NParti)*InvNParti );
297 int PuntiMass = NMax - Punti;
298 for(
int i=0;i<PuntiMass;i++){
299 sw[i] = (st[i] + st[i+Punti]) / 2.;
304 for(
int j=0;j<Potenze;j++){
307 for(
int i=0;i<NMax;i++){
309 for(
int j=0;j<Potenze;j++){
314 for(
int j = 0;j<Potenze;j++){
315 sw[j] = log10(POS(sw[j]))/log10((
double)NMax);
325 for(
int i = 1;i<NMax;i++){
332 Media /= (double)NMax;
333 for(
int i=0;i<NMax;i++){
334 sw[i] = (st[i]-yMin)/(yMass-yMin);
335 if(st[i] > Media && st[i-1] < Media)
338 Campioni = (int)(NMax/(
double)
Zeri);
347 for(
int i = 1;i<NMax;i++){
354 Media /= (double)NMax;
355 for(
int i=0;i<NMax;i++){
356 st[i] = (st[i]-yMin)/(yMass-yMin);
357 if(st[i] > Media && st[i-1] < Media)
360 Campioni = (int)(NMax/(
double)
Zeri);
365 for(
int i = 0;i<NMax;i++){
368 for(
int i=0;i<NMax;i++){
377 for(
int i=0;i<NMax;i++){
379 if(st[i] < m1.
Min) m1.
Min = st[i];
380 if(st[i] > m1.
Max) m1.
Max = st[i];
382 m1.
Uno /= (double)(NMax);
384 for(
int i=0;i<NMax;i++){
385 m1.
Due+=QUAD((st[i]-m1.
Uno));
386 m1.
Tre+=QUAD((st[i]-m1.
Uno))*(st[i]-m1.
Uno);
389 m1.
Due=sqrt(m1.
Due/(
double)(NMax-1));
390 m1.
Tre = pow(m1.
Tre / (
double)(NMax - 2),.33333);
395 for(
int i=0;i<NBin;i++)
397 double Add = 1./(double)NMax;
398 for(
int i=0;i<NMax;i++){
399 int j = (int)((st[i]-Confine[0])/(Confine[1]-Confine[0])*NBin);
400 if(j >= NBin || j < 0)
continue;
403 for(
int i=0;i<NBin;i++){
412 for(
int i=0;i<NBin;i++) Distr[i]=0.;
413 for(
int i=0;i<NMax;i++) m1.
Uno += st[i];
414 m1.
Uno/=(
double)(NMax);
416 double Add = !IfNorm ? 1. : 1./(double)NMax;
417 for(
int i=0;i<NMax;i++){
418 m1.
Due += QUAD(st[i]-m1.
Uno);
420 int j = (int)((st[i]-Confine[0])/(Confine[1]-Confine[0])*NBin);
421 if(j >= NBin || j < 0)
continue;
427 for(
int v=0;v<NBin;v++){
428 if(m1.
yMin > Distr[v]) m1.
yMin = Distr[v];
429 if(m1.
yMax < Distr[v]) m1.
yMax = Distr[v];
431 m1.
Due=sqrt(m1.
Due/((
double)(NMax-1)));
432 m1.
Tre = pow(m1.
Tre,.33333) / (double)(NMax - 2);
442 for(
int i=0;i<NBin;i++) Distr[i] = 0.;
443 for(
int i=0;i<NMax;i++){
444 if(st[i] < m1.
Min) m1.
Min = st[i];
445 if(st[i] > m1.
Max) m1.
Max = st[i];
448 Sum3 += st[i]*st[i]*st[i];
450 double Add = !IfNorm ? 1. : 1./(double)NMax;
452 for(
int i=0;i<NMax;i++){
453 int j = (int)((st[i]-m1.
Min)/(m1.
Max-m1.
Min)*NBin);
454 if(j >= NBin)
continue;
460 for(
int v=0;v<NBin;v++){
461 if(m1.
yMin > Distr[v]) m1.
yMin = Distr[v];
462 if(m1.
yMax < Distr[v]) m1.
yMax = Distr[v];
464 m1.
Uno/=(double)(NMax);
465 m1.
Due = sqrt((Sum2 - QUAD(m1.
Uno)*NMax)/((
double)NMax-1.));
466 m1.
Tre = Sum3-3.*Sum2*m1.
Uno+3.*NMax*QUAD(m1.
Uno)-NMax*m1.
Uno;
467 m1.
Tre = pow( m1.
Tre , 1./3.)/(double)(NMax-1);
471 void Matematica::DistrSample(
double *Px,
double *Py,
int NMax,
double **Distr,
int NBin,
const int NSample,
int IfNorm,
double *yBorder){
474 double Norm[NSample];
477 for(
int n=0;n<NMax;n++){
478 if(xBorder[1] < Px[n]) xBorder[1] = Px[n];
479 if(xBorder[0] > Px[n]) xBorder[0] = Px[n];
481 double Dx = (xBorder[1] - xBorder[0]) > 0. ? 1./(xBorder[1] - xBorder[0]) : 1.;
483 for(
int s=0;s<NSample;s++){
487 for(
int n=0;n<NMax;n++){
488 int sx = (int)((Px[n]-xBorder[0])*Dx*NSample);
489 if(sx < 0 || sx >= NSample)
continue;
490 if(yBorder[0] > Py[n]) yBorder[0] = Py[n];
491 if(yBorder[1] < Py[n]) yBorder[1] = Py[n];
494 yDelta = yBorder[1] - yBorder[0] > 0. ? 1./(yBorder[1] - yBorder[0]) : 1.;
495 for(
int n=0;n<NMax;n++){
496 int sx = (int)((Px[n]-xBorder[0])*Dx*NSample);
497 if(sx < 0 || sx >= NSample)
continue;
498 int by = (int)((Py[n]-yBorder[0])*yDelta*NBin);
499 if(by < 0 || by >= NBin)
continue;
503 for(
int sx=0;sx<NSample;sx++){
505 for(
int by=0;by<NBin;by++){
506 Norm[sx] += Distr[sx][by];
509 for(
int sx=0;sx<NSample;sx++){
510 Norm[sx] = Norm[sx] > 1. ? 1./Norm[sx] : 1.;
511 for(
int by=0;by<NBin;by++){
512 Distr[sx][by] *= Norm[sx];
521 for(
int v=0;v<NBin;v++){
523 m1.
Chi += QUAD( (Distr[v] - dInt[v])/m1.
Due);
532 for(
int i=0;i<NBin;i++){
534 m1.
Chi += QUAD( Distr[i] - dInt[i]);
544 for(
int i=0;i<NMax;i++){
545 if(sx[i] < m1.
Min) m1.
Min = sx[i];
546 if(sx[i] > m1.
Max) m1.
Max = sx[i];
547 m1.
Uno += sx[i]*sy[i];
549 m1.
Due += sy[i] > 0. ? 1./sy[i] : 0.;
556 void Matematica::WeightHisto(
double **hist,
double *Border,
int NBin,
int NHisto,
double tolerance,
double *OrPos,
double *kSpring){
559 cout<<
"do_wham_c <TEMPERATURE> <MIN> <MAX> <BINS> <TOLERANCE> <OUTPUT> <INPUT1> <INPUT2> ....\n\n";
560 cout<<
"TEMPERATURE\t\tTemperature that simulations were performed at\n";
561 cout<<
"MIN\t\t\tMinimum value used in histogram. All histograms must have been computed with matching value.\n";
562 cout<<
"MAX\t\t\tMaximum value used in histogram. All histograms must have been computed with matching value.\n";
563 cout<<
"NBIN\t\t\tNumber of NBin in histogram. All histograms must have been computed with matching value.\n";
564 cout<<
"TOLERANCE\t\tTolerance value used to determine convergence. Value of 0.00001 seems to work fine.\n";
565 cout<<
"OUTPUT\t\t\tName of output file.\n";
566 cout<<
"INPUT\t\t\tNames of input files. These are histogram files produced from the .pdo files using make_histo.pl.\n";
569 double temperature = 1.;
570 double min = Border[0];
571 double max = Border[1];
572 double RT = temperature * 8.314472e-3;
573 ofstream output(
"WeightHistoAnal.dat");
574 double*
F =
new double[NHisto];
575 double* old_F =
new double[NHisto];
576 double* result =
new double[NBin];
577 double* N =
new double[NHisto];
578 double* energy =
new double[NBin];
579 double * numerator =
new double[NBin];
580 double ** expU =
new double *[NHisto];
581 double ** NexpU =
new double *[NHisto];
583 for(
int h = 0;h < NHisto;++h) {
587 expU[h] =
new double[NBin];
588 NexpU[h] =
new double[NBin];
589 for(
int i = 0; i < NBin; ++i) {
591 double pos = (i + 0.5)/NBin * (max - min) + min;
592 expU[h][i] = exp(-(0.5 * kSpring[h] * (pos - OrPos[h]) * (pos - OrPos[h]))/RT);
594 for(
int i=0; i < NBin; ++i) {
595 NexpU[h][i] = N[h] * expU[h][i];
598 for(
int bin = 0; bin < NBin; ++bin) {
602 for(
int bin = 0; bin < NBin; ++bin) {
603 numerator[bin] = 0.0;
604 for(
int file = 0; file < NHisto; ++file) {
605 numerator[bin] += hist[file][bin];
608 cerr<<
"Done loading histograms. Doing WHAM."<<endl;
613 for(
int bin=0; bin < NBin; ++bin) {
615 for(
int file = 0; file < NHisto; ++file) {
616 denom += NexpU[file][bin] / F[file];
619 if(denom > 1.e9) IfContinue = 0;
620 denom = denom <= 0.0000001 ? 1. : denom;
621 result[bin] = numerator[bin] / denom;
624 for(
int bin=0; bin < NBin; ++bin) norm += result[bin];
625 norm = norm > 0. ? norm : 1.;
626 for(
int bin=0; bin < NBin; ++bin) result[bin] /= norm;
627 for(
int file = 0; file < NHisto; ++file) {
629 for (
int bin = 0; bin < NBin; ++bin){
630 Z += result[bin] * expU[file][bin];
632 if(isnan(Z)){IfContinue = 0;
break;}
633 if(isinf(Z)){IfContinue = 0;
break;}
634 if(Z <= 1.e-9){IfContinue = 0;
break;}
635 double dZ = (max - min)/(
double)NBin;
637 old_F[file] = F[file];
641 if(!IfContinue)
break;
643 for(
int file = 0; file < NHisto; ++file) {
644 temp = RT * fabs(log(F[file]) - log(old_F[file]));
649 cerr<<
"Step: "<<count<<
"\tTolerance: "<<d_max<<endl;
651 if(count == 2)
break;
652 }
while(d_max > tolerance);
653 for(
int bin = 0; bin < NBin; ++bin) {
656 for(
int bin = 0; bin < NBin; ++bin) {
657 double res = result[bin] > 0. ? log(result[bin]) : 0.;
658 energy[bin] = -RT*res;
659 double pos = (double)bin / NBin * (max - min) + min;
660 output << pos <<
" " << result[bin] <<
" " << energy[bin] << endl;
669 for(
int current = 0;current < NHisto;++current) {
670 delete [] expU[current];
671 delete [] NexpU[current];
677 for(
int i=0;i<NMax;i++){
678 for(
int j=i;j>0;j--){
679 if(Sign[j] < Sign[j-1])
687 double Temp = Sign[j];
693 for(
int d=0;d<NDim;d++){
694 Temp[d] = s[si*NDim+d];
695 s[si*NDim+d] = t[ti*NDim+d];
696 t[ti*NDim+d] = Temp[d];
700 for(
int i=0;i<NMax;i++){
701 for(
int j=i;j>0;j--){
702 if(Sign[j] < Sign[j-1])
715 double *st1 = (
double *)calloc(NMax,
sizeof(
double));
716 for(
int n=0;n<NMax;n++){
719 for(
int n=0;n<NMax;n++){
721 if(n+WIndex[0] < 0 || n+WIndex[NWeight-1] >= NMax){
725 for(
int w=0;w<NWeight;w++){
726 st[n] += sw[w]*st1[n+WIndex[w]];
732 int Half = NWeight/2;
734 for(
int w=0;w<NWeight;w++){
736 double x = CutOff*(w-Half)/(
double)NWeight;
738 double Gauss = exp(-r2*.5/SQR(Sigma));
742 for(
int w=0;w<NWeight;w++){
747 for(
int i=0;i<NMax;i++){
752 if(strlen(cmd) == 0)
return x;
762 for(
int i=0;i<NOp;i++){
763 cPart = strpbrk(cPart+1,
"*/^+-");
764 if(cPart == NULL)
break;
765 Comm[i] = (int)*cPart;
766 PosOp[i] = cPart-cmd+1;
767 PosNum[nOp] = PosOp[i];
771 PosNum[nOp++] = strlen(cmd);
773 for(
int i=0;i<nOp-1;i++){
774 if( *(cmd+PosNum[i]) ==
'x'){
777 if( *(cmd+PosNum[i]) ==
'y'){
781 sscanf(cmd+PosNum[i],
"%lf",&Numb);
814 if(strlen(cmd) == 0)
return st[0][n];
824 for(
int i=0;i<NOp;i++){
825 cPart = strpbrk(cPart+1,
"*/^+-");
826 if(cPart == NULL)
break;
827 Comm[i] = (int)*cPart;
828 PosOp[i] = cPart-cmd+1;
829 PosNum[nOp] = PosOp[i];
833 PosNum[nOp++] = strlen(cmd);
836 for(
int i=0;i<nOp-1;i++){
837 if( *(cmd+PosNum[i]) ==
'x'){
838 sscanf(cmd+PosNum[i]+1,
"%d",&NVar);
842 sscanf(cmd+PosNum[i],
"%lf",&Numb);
double RandDiscrProb(double *Prob, int NBin)
Random number following a discrete probability.
void Radice(double *st, double *sw, int N)
Compute the root of the signal.
void FillWeightGauss(double *st, int *WIndex, int NWeight, double CutOff, double Sigma)
Fill the weight array with a gaussian fuction.
int Zeri(double a, double b, double *Radici, int NRadici)
Find the.
MOMENTI DistribuzioneMaxwell(const double *st, int NMass, double *Intervalli, double *dInt, int Valori, int IfNorm)
Look for the Maxwellian distribution.
bool InizializzaGaussiano(double Scarto, int N)
Initialize the Gaussian number generator.
void ConvWeight(double *st, int NMax, double *sw, int *WIndex, int NWeight)
Convolute with a weight.
double yMin
Minimum of the ordinate.
MOMENTI Distribuzione(const double *st, int NMass)
Moments of a signal.
MOMENTI WeightAverage(const double *sx, const double *sy, int NMax)
Calculate the weighted average.
double yMax
Maximum of the ordinate.
Moments of a distribution.
void Sort(double *Sign, int NMass)
Sort.
void WeightHisto(double **hist, double *Border, int NBin, int NHisto, double tolerance, double *OrPos, double *kSpring)
Weighted histogram analysis.
MOMENTI DistrErr(const double *st, int NMass, double *Intervalli, double *Err, int Valori, double *Confine, int IfNorm)
Moments and histogram of a signal between two values.
double Casuale()
Random uniform number.
int NormalizeArea(double *st, int NMass)
Normalize.
void Swap(int i, int j, double *Sign)
Swap to indices.
void DistrSample(double *Px, double *Py, int NMax, double **Distr, int NBin, const int NSample, int IfNorm, double *xBound)
Compare the distribution of a sample of data.
MOMENTI DistribuzioneGauss(const double *st, int NMass, double *Intervalli, double *dInt, int Valori, int IfNorm)
Look for the Gaussian distribution.
int CorrelaDuePunti(double *st, int NMass, double *sw, int Punti)
Two points correlation.
void Spettro2d(double *st, double *sw, int NMass)
Compute the 2d spectrum of.
double Gauss(double Media, double Scarto, double x)
Gaussian.
double Gaussiano(double Media, double Scarto)
Gaussian random number.
double F(double TD, double T)
Boh.
double ExecFormula(double x, double y, char *cmd)
Execute a formula.
double Norm(double *st, int NMass)
Norm of an array.
int Num
Number of points considered.
void Autocor(bool *st, double *sAuto, int N)
Compute the autocorrelation of a boolean signal.
void ExecCommand(double *st, double *st1, int NMass, char *cmd)
Execute a command defined in string.
int Normalizza(double *st, int NMass)
Normalize.
double Delta
Step interval.
void Spettro(double *st, double *sw, int NMass)
Compute the spectrum.
void Autosimilarita(double *st, int NMass, double *sw, int Valori)
Self similarity.
void MediaMobile(double *st, int NMass, double *sw, int Parti)
Running average.