1/*___________________________________________________________________________________________________________*/ 2 3function PrintTableToFile (dataMatrix, titleMatrix, promptOrNot) 4{ 5 SetDialogPrompt ("Export tab separated data to:"); 6 7 if (promptOrNot) 8 { 9 fprintf (baseResPath, CLEAR_FILE, KEEP_OPEN); 10 } 11 12 fprintf (baseResPath, "Sample\t", titleMatrix[0][0]); 13 14 for (counter1=1; counter1<Columns(titleMatrix); counter1 = counter1+1) 15 { 16 fprintf (baseResPath, "\t", titleMatrix[0][counter1]); 17 } 18 19 fprintf (baseResPath, "\n"); 20 21 for (counter1=0; counter1<Rows(dataMatrix); counter1 = counter1 + 1) 22 { 23 fprintf (baseResPath,counter1+1); 24 for (counter2 = 0; counter2 < Columns (titleMatrix); counter2 = counter2+1) 25 { 26 fprintf (baseResPath,"\t",dataMatrix[counter1][counter2]); 27 } 28 fprintf (baseResPath,"\n"); 29 } 30 if (promptOrNot) 31 { 32 fprintf (baseResPath, CLOSE_FILE); 33 } 34 35 return 1; 36} 37/*___________________________________________________________________________________________________________*/ 38 39 40 41fprintf (stdout, "\n\nComputing likelihood weights (this may take a while)...\n"); 42 43scores = {SAMPLE_N,2}; 44bestScore = -1e25; 45extras = {}; 46 47 48timer = Time (1); 49 50extraPopulator = 0; 51 52 53if (MPI_NODE_ID == 0 && Abs(_ADD_TO_SRS_) == 2 && Type (_ADD_TO_SRS_) == "AssociativeList") 54{ 55 extraPopulator = 1; 56} 57 58if (MPI_NODE_COUNT > 1 && MPI_NODE_ID == 0 && SAMPLE_N >= MPI_NODE_COUNT) 59{ 60 PRESERVE_SLAVE_NODE_STATE = 1; 61 samplesPerNode = SAMPLE_N$MPI_NODE_COUNT; 62 nodeRanges = {MPI_NODE_COUNT,2}; 63 done = 0; 64 for (itCount = 1; itCount < MPI_NODE_COUNT; itCount = itCount + 1) 65 { 66 if (extraPopulator) 67 { 68 LF_NEXUS_EXPORT_EXTRA = _ADD_TO_SRS_["INIT"]; 69 LF_EXTRA = "(MPI_NEXUS_FILE_RETURN[1])[itCount]=" + _ADD_TO_SRS_["EVAL"] + ";\n"; 70 } 71 else 72 { 73 LF_NEXUS_EXPORT_EXTRA = ""; 74 LF_EXTRA = ""; 75 } 76 77 LF_NEXUS_EXPORT_EXTRA += "LFCompute("+ 78 LF_NAME+ 79 ",LF_START_COMPUTE);SAMPLE_N="+ 80 samplesPerNode+ 81 ";assignmentString=\""+ 82 assignmentString+"\";"; 83 84 from = done; 85 to = done+samplesPerNode-1; 86 mpiGenSamples = generatedSamples[{{from,0}}][{{to,Columns(generatedSamples)-1}}]; 87 LF_NEXUS_EXPORT_EXTRA = LF_NEXUS_EXPORT_EXTRA + "generatedSamples="+mpiGenSamples+";MPI_NEXUS_FILE_RETURN={};MPI_NEXUS_FILE_RETURN[0]={SAMPLE_N,1};MPI_NEXUS_FILE_RETURN[1]={};\nfor (itCount = 0; itCount < SAMPLE_N; itCount = itCount + 1)\n\t" 88 +"{ExecuteCommands (assignmentString);LFCompute ("+LF_NAME+",lfWeight);(MPI_NEXUS_FILE_RETURN[0])[itCount]=lfWeight;`LF_EXTRA`}\nLFCompute ("+LF_NAME+",LF_DONE_COMPUTE);\n"; 89 90 ExecuteCommands ("Export(lfExport,"+LF_NAME+");"); 91 92 MPISend (itCount, lfExport); 93 94 done = done+samplesPerNode; 95 nodeRanges [itCount][0] = from; 96 nodeRanges [itCount][1] = to; 97 fprintf (stdout, "[SENT RANGE ", from+1, " - ", to + 1, " TO MPI NODE ", itCount, "]\n"); 98 } 99 ExecuteCommands ("LFCompute ("+LF_NAME+",LF_START_COMPUTE);"); 100 toDo = SAMPLE_N-done+1; 101 if (extraPopulator) 102 { 103 ExecuteCommands (_ADD_TO_SRS_["INIT"]); 104 } 105 for (itCount = done; itCount < SAMPLE_N; itCount = itCount + 1) 106 { 107 ExecuteCommands (assignmentString+"LFCompute ("+LF_NAME+",lfWeight);"); 108 scores[itCount][0] = lfWeight; 109 if (extraPopulator) 110 { 111 extras[itCount] = Eval (_ADD_TO_SRS_["EVAL"]); 112 } 113 if ((1+itCount-done) % 100 == 0) 114 { 115 fprintf (stdout, itCount+1, "/", SAMPLE_N, " evaluations done. Estimated remaining time: ",Format (((SAMPLE_N-itCount-1)/(itCount+1-done))*(Time(1)-timer),5,2)," seconds \n"); 116 } 117 } 118 for (itCount = 1; itCount < MPI_NODE_COUNT; itCount = itCount + 1) 119 { 120 MPIReceive (-1,fromNode,res); 121 ExecuteCommands ("mpiRes="+res); 122 for (sampleCount = 0; sampleCount < Rows(mpiRes[0]); sampleCount = sampleCount + 1) 123 { 124 scores[sampleCount+nodeRanges[fromNode][0]][0] = (mpiRes[0])[sampleCount]; 125 extras[sampleCount+nodeRanges[fromNode][0]][0] = (mpiRes[1])[sampleCount]; 126 } 127 fprintf (stdout, "[GOT RANGE ", nodeRanges[fromNode][0]+1, " - ", nodeRanges[fromNode][0]+Rows(mpiRes[0]), " FROM MPI NODE ", fromNode, "]\n"); 128 } 129 for (itCount = 0; itCount < SAMPLE_N; itCount = itCount + 1) 130 { 131 bestScore = Max(bestScore,scores[itCount][0]); 132 } 133 fprintf (stdout, "[BEST SCORE = ",bestScore,"]\n"); 134} 135else 136{ 137 ExecuteCommands ("LFCompute ("+LF_NAME+",LF_START_COMPUTE);"); 138 if (extraPopulator) 139 { 140 ExecuteCommands (_ADD_TO_SRS_["INIT"]); 141 } 142 for (itCount = 0; itCount < SAMPLE_N; itCount = itCount + 1) 143 { 144 ExecuteCommands (assignmentString+"LFCompute ("+LF_NAME+",lfWeight);"); 145 scores[itCount][0] = lfWeight; 146 if (extraPopulator) 147 { 148 extras[itCount] = Eval (_ADD_TO_SRS_["EVAL"]); 149 } 150 151 if (sMarginals>1) 152 { 153 ExecuteCommands ("ConstructCategoryMatrix(catmat,"+LF_NAME+",COMPLETE);"); 154 for (sit=0; sit < Columns(catVarList); sit = sit+1) 155 { 156 ExecuteCommands ("GetInformation(catVarInfo," + catVarList[sit] + ");"); 157 fprintf (marginalOutFileAll,"\n",catVarInfo); 158 } 159 fprintf (marginalOutFileAll,"\n",catmat); 160 fprintf (marginalOutFileLF, "\n",lfWeight); 161 } 162 163 if (lfWeight > bestScore) 164 { 165 bestScore = lfWeight; 166 } 167 168 if ((1+itCount) % 100 == 0) 169 { 170 fprintf (stdout, itCount+1, "/", SAMPLE_N, " evaluations done. Estimated remaining time: ",Format (((SAMPLE_N-itCount-1)/(itCount+1))*(Time(1)-timer),5,2)," seconds \n"); 171 } 172 } 173} 174 175/* restore to original values */ 176 177for (k=0; k<varCount; k=k+1) 178{ 179 aKey = usedVars[k]; 180 ExecuteCommands (aKey + "=" + stashedValues[k][0] + ";"); 181} 182 183_resamplerReturn = {}; 184_resamplerReturn["RAW_SAMPLES"] = generatedSamples; 185_resamplerReturn["SCORES"] = scores; 186_resamplerReturn["EXTRAS"] = 0; 187 188if (extraPopulator) 189{ 190 _resamplerReturn["EXTRAS"] = {}; 191} 192 193ExecuteCommands ("LFCompute ("+LF_NAME+",LF_DONE_COMPUTE);"); 194 195fprintf (stdout, "\n\nResampling...\n"); 196coord = 0; /* this is the current sum of all weighting factors */ 197for (itCount = 0; itCount < SAMPLE_N; itCount = itCount + 1) 198{ 199 scores[itCount][0] = Exp(scores[itCount][0]-bestScore); 200 coord = coord + scores[itCount][0]; 201} 202 203N_eff = 0; 204 205for (itCount = 0; itCount < SAMPLE_N; itCount = itCount + 1) 206{ 207 N_eff = N_eff + (scores[itCount][0]/coord)^2; 208} 209 210fprintf (stdout, "Estimated effective sample size: ", 1/N_eff, "\n\n"); 211 212j = 0; 213sampledPoints = {SAMPLE_M,varCount+1}; 214 215timer = Time (1); 216 217if (sMarginals) 218{ 219 ExecuteCommands ("LFCompute ("+LF_NAME+",LF_START_COMPUTE);"); 220} 221 222/* resampling loop */ 223 224for (itCount = 0; itCount < SAMPLE_M; itCount = itCount + 1) 225{ 226 sum = Random(0,1); 227 runningSum = 0; 228 229 local_scores = scores; 230 local_scale = 0; 231 232 for (k=0; k<SAMPLE_N;k=k+1) 233 { 234 if (scores[k][1] == 0) 235 { 236 local_scores[k][0] = scores[k][0]/(coord-(itCount+1)*scores[k][0]); 237 if (local_scores[k][0] < 0) 238 { 239 local_scores = scores; 240 local_scale = coord; 241 break; 242 } 243 local_scale = local_scale + local_scores[k][0]; 244 } 245 } 246 247 lastValidIndex = 0; 248 249 for (k=0; k<SAMPLE_N && runningSum < sum;k=k+1) 250 { 251 if (local_scores[k][1] == 0) 252 /* can still sample this point */ 253 { 254 lastValidIndex = k; 255 runningSum = runningSum + local_scores[k][0] / local_scale; 256 } 257 } 258 259 sampledPoints[j][0] = Log(scores[lastValidIndex][0])+bestScore; 260 261 for (k = 0; k < varCount; k=k+1) 262 { 263 sampledPoints[j][k+1] = generatedSamples[lastValidIndex][k]; 264 } 265 266 if (extraPopulator) 267 { 268 (_resamplerReturn["EXTRAS"])[j] = extras[lastValidIndex]; 269 } 270 271 coord = coord - scores[lastValidIndex][0]; 272 scores[lastValidIndex][1] = 1; 273 j = j+1; 274 275 if (sMarginals) 276 { 277 sit = itCount; 278 itCount = lastValidIndex; 279 ExecuteCommands (assignmentString+"LFCompute ("+LF_NAME+",lfWeight);ConstructCategoryMatrix(catmat,"+LF_NAME+",COMPLETE);"); 280 itCount = sit; 281 for (sit=0; sit < Columns(catVarList); sit = sit+1) 282 { 283 ExecuteCommands ("GetInformation(catVarInfo," + catVarList[sit] + ");"); 284 fprintf (marginalOutFile,"\n",catVarInfo); 285 } 286 fprintf (marginalOutFile,"\n",catmat); 287 } 288 289 if ((1+itCount) % 100 == 0) 290 { 291 fprintf (stdout, itCount+1, "/", SAMPLE_M, " samples drawn. Estimated remaining time: ",Format (((SAMPLE_M-itCount-1)/(itCount+1))*(Time(1)-timer),5,2)," seconds \n"); 292 } 293} 294 295 296if (sMarginals) 297{ 298 for (k=0; k<varCount; k=k+1) 299 { 300 aKey = usedVars[k]; 301 ExecuteCommands (aKey + "=" + stashedValues[k][0] + ";"); 302 } 303 ExecuteCommands ("LFCompute ("+LF_NAME+",LF_DONE_COMPUTE);"); 304} 305 306 307labelMatrix = {1,varCount+1}; 308 309labelMatrix[0] = "-Log(L)"; 310 311_resamplerReturn ["LABELS"] = {}; 312(_resamplerReturn ["LABELS"])[0] = labelMatrix[0]; 313 314for (k = 0; k < varCount; k=k+1) 315{ 316 labelMatrix[k+1] = usedVars[k]; 317 (_resamplerReturn ["LABELS"])[k+1] = usedVars[k]; 318} 319 320/* open a chart window */ 321 322PrintTableToFile (sampledPoints, labelMatrix, 1); 323 324OpenWindow (CHARTWINDOW,{{"Sampled parameter values"} 325 {"labelMatrix"}, 326 {"sampledPoints"}, 327 {""}, 328 {"Index"}, 329 {""}, 330 {""}, 331 {""}, 332 {""}, 333 {"3"}}, 334 "(SCREEN_WIDTH-200);(SCREEN_HEIGHT-200);80;100"); 335 336 337_resamplerReturn ["VALUES"] = sampledPoints; 338 339return _resamplerReturn; 340