#define LastUpd "Last Update 06-01-1996" #define FREE_ARG char* #include #include #include #include #include #include #include "SW.h" #include "free_ragged.h" #include "../utilities/nrutil.h" #include "../utilities/ranmar.h" int main(int argc, char *argv[]) { int N; /* Number of points */ unsigned char *Spin; /* Spin values: 0,1,2, ..., Q-1 */ float **P; /* Deletion Probabilities (for a satisfied */ /* bond; i.e. if S[i] = S[j] otherwise P = 1 */ /* It is no longer simmetric. */ float **J; /* Contens the interactions with its */ /* neigboors */ unsigned short **NK; /* Contens the neighboors of each point. */ /* NK[i][j], i=1,N; j=1,... are the */ /* label of neighboring points for each i */ /* the last element in the array NK[i] is */ /* set MAX_USHRT to indicate its end. */ unsigned short **KN; /* KN[i][k] = m, means that point i us the */ /* m-th neighbor of point j = N[i][k] */ char **Bond; /* Bond[i][k] takes value 1 (0) if bond */ /* between spins i and its k-th neighbor is */ /* frozen (deleted). */ unsigned short *ClusterSize; /* ClusterSize[n] contains the number of */ /* points belonging to cluster n in the */ /* current configuration. The clusters are */ /* ordered fron the biggest (n=1) to the */ /* smallest */ unsigned short *Block; /* Block[i] Is the number of the cluster to */ /* which Spin[i] belongs */ float **X; /* X[d][i] is the d-th coordinate of the */ /* i-th point */ unsigned short **CorrN; /* Is the Correlation function; i.e. the */ /* estimated probability that two spins be */ /* to the same cluster */ /* exists between spins i and j, due the */ /* "whole system" (CorrN). */ float *Size1; /* Sum of all clustersizes for a given T */ float *Size2; /* Sum of all clustersizes^2 for a given T */ int Q; /* Potts Spin order; Si = 0, 1, ..., Q-1 */ float T; /* Temperature */ float Tmin,Tmax,dT; int D; /* Dimension of the point vectors, D=0 means */ /* that the distances are provided (instead */ /* the coordinates of the vectors). */ int nc; /* present number of clusters */ int nb; /* present number of frozen bonds */ float nb1,nb2,nc1,nc2; /* sum of nb, nb^2, nc, nc^2 */ int NS; /* Number of clusters size to be printed out */ /* calculating the magnetization */ int i; int ncy, cyc; float ChD; /* Characteristic Distance */ float Jmed; /* average interaction between neighbors */ int nT; float nn; /* Average number of edges: edges per site */ float Tc; float thN; /* distintos thresholds */ float th_MIN, th_MAX,dth; /* delta-threshold */ int nth; /* numbers the different threshold files */ int wrLabels, wrClusters; int T_updp; /* for choosing the temperatures: */ /* 0 -> multiplicative T_n = Tmin * dt^n */ /* 1 -> additive T_n = Tmin + d * n */ char OutFile[80]; char DataFile[80]; /* file containing data infrmation, */ /* if D==0 -> proximity matrix */ /* D >0 -> vector ccordinates */ char EdgesFile[80]; char thFile[80]; /* threshold file */ float SWfract; /* fraction of the SW trials that are used to */ /* compute the averages. The first */ /* (1-fraction) cycles are considered as */ /* transient states */ float Ener[2]; /* Ener[0] = ; Ener[1] = Var(Energy) */ float m1, m2; /* Actual Magnetization */ float M1, M2; /* Magnetization cumulants */ int nq = 0; /* random sampling for -> magnetization3 */ float M_fraction_points=0.9; /* minimal fraction points to be considered */ float xi[2]; /* xi[0] low T suceptibility, xi[1] high T */ float fk; int ij,kl; /* seed for the random number generator */ start_time(); rmarin2(1802,9373); ReadRunParam(argv[1],OutFile,DataFile,EdgesFile,&N,&D,&Q,&cyc,&Tmin,&Tmax,&dT, &T_updp,&SWfract,&th_MIN,&th_MAX,&dth,&NS,&wrLabels,&wrClusters); ShowParameters(OutFile,DataFile,EdgesFile,N,D,Q,cyc,Tmin,Tmax,dT,T_updp,SWfract, th_MIN,th_MAX,dth,NS,wrLabels,wrClusters); /* --------------------------------------------------------------*/ /* shows parameters (on standard output) */ /*---------------------------------------------------------------*/ Check(N,D,Q,NS,cyc); /* --------------------------------------------------------------*/ /* checks parameters (the function - below) */ /*---------------------------------------------------------------*/ /* ------------------------------------------------------------- */ /* Setup of the physical problem: */ /* points interactions and initial state */ /* ------------------------------------------------------------- */ NK = Edge1(EdgesFile, N, &nn); KN = Edge2(N, NK); if(D == 0) J = Interaction_Dist(DataFile, N, NK, nn, &ChD); else { X = matrix(0,(long)(D-1),0,(long)(N-1)); ReadData(DataFile,D,N,X); J = Interaction(N, D, NK, nn, X, &ChD); free_matrix(X,0,(long)(D-1),0,(long)(N-1)); } Spin = ucvector(0,(long)(N-1)); Tc = CriticalTemp(Q,N,J,NK); Jmed = AverageInteraction(N, J, NK); PrintParameters(OutFile,N,D,Q,cyc,SWfract,ChD,Jmed,nn,th_MIN,th_MAX,Tc); InitialSpinConfig(N,Spin,Q); for(T = Tmin, nT = 0; T < Tmax; nT++){ /* T LOOP */ if(nT != 0){ if (T_updp == 1) T += dT; else if (T_updp == 0) T *= dT; else nrerror("invalid t_updp value\n"); } fprn_time(OutFile); ij = IUnif(0,RANMAR_IJ_MAX); kl = IUnif(0,RANMAR_KL_MAX); rmarin2(ij,kl); /* Memory allocations: */ CorrN = CorrMatrix(N,NK); Bond = BondMatrix(N,NK); P = ProbMatrix(N,NK); ClusterSize = usvector(0,(long)(N-1)); Block = usvector(0,(long)(N-1)); Size1 = vector(0,(long) NS); Size2 = vector(0,(long) NS); DeletionProbabilities(N,NK,T,J,P); /* Transient of Monte Carlo (not included in averages) */ for( ncy = 0; ncy < (1-SWfract) * cyc ; ncy++){ nb = SetBond(N,P,Spin,Bond,NK,KN); nc = Coarsening(N,Bond,Block,NK,ClusterSize); NewSpinConfig(N,Spin,Block,nc,Q); } /* Estimation of Monte Carlo averages */ for(i = 0; i < NS; i++) Size1[i] = Size2[i] = zero; xi[1] = M1 = M2 = zero; nb1 = nb2 = nc1 = nc2 = 0.0000; for(ncy = 0; ncy <= SWfract * cyc ; ncy++){ nb = SetBond(N,P,Spin,Bond,NK,KN); nc = Coarsening(N,Bond,Block,NK,ClusterSize); nb1 += (float)nb; nb2 += (float)(nb*nb); nc1 += (float)nc; nc2 += (float)(nc*nc); NewSpinConfig(N,Spin,Block,nc,Q); OrderingClusters(N,nc,Block,ClusterSize,NS); for(i = 0; i < NS; i++){ Size1[i] += (float)(ClusterSize[i]); Size2[i] += (float)(ClusterSize[i] * ClusterSize[i]); } GlobalCorrelation(N,CorrN,NK,Block); m1 = magnetization3(N, Q, nc, ClusterSize, nq, M_fraction_points, &m2); M1 += m1; M2 += m2; for(i = 0; i < nc; i++) { fk = (float)ClusterSize[i] / (float)N; xi[1] += fk * fk; } } ncy--; nb1 /= ncy; nb2 = nb2/ncy - nb1 * nb1; nc1 /= ncy; nc2 = nc2/ncy - nc1 * nc1; PrintSWaverages(OutFile,".SWa",NS,nT,T,N,nb1,nb2,nc1,nc2); ClusterAverage(ncy,NS,Size1,Size2); PrintAverages(OutFile,".ave",NS,nT,T,Jmed,nb1/N,Size1); /* PrintAverages(OutFile,".var",NS,nT,T,Jmed,nb2/N,Size2); */ /* PrintBoth(OutFile,".bot",T,Jmed,Size1,Size2,NS); */ /* Energy(N, Q, ncy, J, CorrN, NK, Ener); */ /* PrintEnergy(OutFile,nT,T,Jmed,Ener); */ susceptibility(N,Q,ncy,xi,M1,M2); PrintMagnet(OutFile,nT,T,Jmed,M1/ncy,xi); /* PrintCorrN(N,CorrN,NK,ncy); */ free_float_ragged(P, N); free_char_ragged(Bond, N); free_vector(Size1,0,(long) NS); free_vector(Size2,0,(long) NS); /* distintos thresholds son usados para encontrar el "core" */ /* de los clusters. El mejor threshold es aquel que mejor */ /* aproxima el tama#o de los SW-clusters */ nth = 0; for(thN = th_MIN; thN <= th_MAX; thN += dth){ nth ++; SLabel(nth, ".th_", thFile); nc = Thresholding(N,ncy,NS,thN,CorrN,NK,Block,ClusterSize); PrintSizes(OutFile,thFile,NS,nT,T,Jmed,nc,ClusterSize); /* if (wrClusters != 0) */ /* WriteClusters(DataFile,N,D,Block,NS,OutFile,nT,thFile); */ /* MY CHANGE!!! */ /* if (wrLabels != 0) */ /* WriteLables(N, Block, OutFile,thFile, nT); */ /* WriteLables1(N, Block, ClusterSize, NS, OutFile,thFile, nT); */ } /* threshod + directed graph con distintos thresholds, por ahora a mano */ nth = 0; for(thN = th_MIN; thN <= th_MAX; thN += dth){ nth ++; SLabel(nth, ".dg_", thFile); nc = directed_graph(N,ncy,NS,thN,CorrN,NK,KN,Block,ClusterSize); PrintSizes(OutFile,thFile,NS,nT,T,Jmed,nc,ClusterSize); /* if (wrClusters != 0) WriteClusters(DataFile,N,D,Block,NS,OutFile,nT,thFile); E.P. */ /* MY CHANGE!!!! */ if (wrLabels != 0) /* WriteLables(N, Block, OutFile,thFile, nT); */ WriteLables1(N, Block, ClusterSize, NS, OutFile,thFile, nT); } free_us_ragged(CorrN, N); free_usvector(Block, 0, (long)(N-1)); free_usvector(ClusterSize, 0, (long)(N-1)); } /* END OF T LOOP */ free_float_ragged(J, N); free_us_ragged(NK, N); free_us_ragged(KN, N); free_ucvector(Spin, 0, (long)(N-1)); fprn_time(OutFile); return(0); } /* -------------------------------------------------------------------- */ /* Check the input parameters related to array dimensionalization */ /* in order to prevent errors in the memory assignation. */ /* -------------------------------------------------------------------- */ void Check(int N, int D, int Q, int NS, int cy) { int error = 0; if (N < 1) { printf(" N = %d is not alloewd\n", N); error = 1;} if (D < 0) { printf(" D = %d is not allowed\n", D); error = 1;} if (Q > UCHAR_MAX) {printf("Q = %d while Qmax = %d\n",Q,UCHAR_MAX); error = 1;} if (cy > USHRT_MAX-1) { printf("cyc = %d while SHRT_MAX = %d\n", cy,SHRT_MAX); printf("if cyc must be bigger than %d then change: \n"); printf(" unsigned short Corr -> unsignrd int Corr\n", cy); error = 1; } if (N > USHRT_MAX-1) { printf("N = %d while USHRT_MAX = %d\n", N,USHRT_MAX); printf("Many 'short' arrays may lead to wrong results\n"); error = 1; } if(error != 0) nrerror("Check >>> input error detected"); return; } /* ----------------------------------------------------------- */ /* given the two "cumulants" the mean value and the VARIANCE */ /* are returned */ /* ----------------------------------------------------------- */ void ClusterAverage(int ncy, int NS, float *Size1, float *Size2) { int i; for(i = 0; i < NS; i++) Size1[i] /= (float)ncy; for(i = 0; i < NS; i++) Size2[i] /= (float)ncy; for(i = 0; i < NS; i++) Size2[i] -= (Size1[i] * Size1[i]); }