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