Context Navigation

source: sans/XOP_Dev/MonteCarlo/MonteCarlo2.c @ 505

Last change on this file since 505 was 458, checked in by srkline, 14 years ago

Added poor man's threading to the MonteCarlo? calculation.

My guess is that the ran() function from NR is not thread safe (it is non-reentrant). So I simply duplicated Monte_SANSX to Monte_SANSX2, where each incarnation uses a different random number generator, either ran1() or ran3(). This means that currently only two processors are supported. Not a big deal. At least it works.

File size: 15.2 KB

Line
1	/*
2	* MonteCarlo.c
3	* SANSAnalysis
4	*
5	* Created by Steve Kline on 10/16/08.
6	* Copyright 2008 __MyCompanyName__. All rights reserved.
7	*
8	*/
9
10
11	#include "XOPStandardHeaders.h" // Include ANSI headers, Mac headers, IgorXOP.h, XOP.h and XOPSupport.h
12	#include "MonteCarlo.h"
13
14	static int gCallSpinProcess = 1; // Set to 1 to all user abort (cmd dot) and background processing.
15
16	//////////
17	// PROGRAM Monte_SANS
18	// PROGRAM simulates multiple SANS.
19	// revised 2/12/99 JGB
20	// added calculation of random deviate, and 2D 10/2008 SRK
21
22	// N1 = NUMBER OF INCIDENT NEUTRONS.
23	// N2 = NUMBER INTERACTED IN THE SAMPLE.
24	// N3 = NUMBER ABSORBED.
25	// THETA = SCATTERING ANGLE.
26
27	// works fine in the single-threaded case.
28	//
29	/// currently crashes if threaded. apparently something here is either doing an unknown callback, or is accessing
30	// a bad place in memory.
31	//
32	// the operations SpinProcess() and WaveHandleModified() are callbacks and MUST be commented out before threading.
33	// - some locks are non-existent
34	// - supposedly safe wave access routines are used
35	//
36	// random number generators are not thread-safe, and can give less than random results, but is this enough to crash?
37	// -- a possible workaround is to define multiple versions (poor man's threading)
38	//
39	//
40	//
41
42	// version X uses ran3
43	// version X2 uses ran1
44
45	int
46	Monte_SANSX2(MC_ParamsPtr p) {
47	double inputWave; / pointer to double precision wave data */
48	double ran_dev; / pointer to double precision wave data */
49	double nt; / pointer to double precision wave data */
50	double j1; / pointer to double precision wave data */
51	double j2; / pointer to double precision wave data */
52	double nn; / pointer to double precision wave data */
53	// double MC_linear_data; / pointer to double precision wave data */
54	double results; / pointer to double precision wave data */
55	double result; //return value
56
57	long imon;
58	double r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,sig_sas;
59	long ind,index,n_index;
60	double qmax,theta_max,q0,zpow;
61	long n1,n2,n3;
62	double dth,zz,xx,yy,phi;
63	double theta,ran,ll,rr,ttot;
64	long done,find_theta,err; //used as logicals
65	long xPixel,yPixel;
66	double vx,vy,vz,theta_z;
67	double sig_abs,ratio,sig_total;
68	double testQ,testPhi,left,delta,dummy,pi;
69	double sigabs_0,num_bins;
70	long NSingleIncoherent,NSingleCoherent,NScatterEvents,incoherentEvent,coherentEvent;
71	long NDoubleCoherent,NMultipleScatter,isOn,xCtr_long,yCtr_long;
72
73	// for accessing the 2D wave data, direct method (see the WaveAccess example XOP)
74	waveHndl wavH;
75	int waveType,hState;
76	long numDimensions;
77	long dimensionSizes[MAX_DIMENSIONS+1];
78	char* dataStartPtr;
79	long dataOffset;
80	long numRows, numColumns,numRows_ran_dev;
81	double dp0, dp, value[2]; // Pointers used for double data.
82	long seed;
83	long indices[MAX_DIMENSIONS];
84
85	char buf[256];
86
87	/* check that wave handles are all valid */
88	if (p->inputWaveH == NIL) {
89	SetNaN64(&p->result); /* return NaN if wave is not valid */
90	return(NON_EXISTENT_WAVE);
91	}
92	if (p->ran_devH == NIL) {
93	SetNaN64(&p->result); /* return NaN if wave is not valid */
94	return(NON_EXISTENT_WAVE);
95	}
96	if (p->ntH == NIL) {
97	SetNaN64(&p->result); /* return NaN if wave is not valid */
98	return(NON_EXISTENT_WAVE);
99	}
100	if (p->j1H == NIL) {
101	SetNaN64(&p->result); /* return NaN if wave is not valid */
102	return(NON_EXISTENT_WAVE);
103	}
104	if (p->j2H == NIL) {
105	SetNaN64(&p->result); /* return NaN if wave is not valid */
106	return(NON_EXISTENT_WAVE);
107	}
108	if (p->nnH == NIL) {
109	SetNaN64(&p->result); /* return NaN if wave is not valid */
110	return(NON_EXISTENT_WAVE);
111	}
112	if (p->MC_linear_dataH == NIL) {
113	SetNaN64(&p->result); /* return NaN if wave is not valid */
114	return(NON_EXISTENT_WAVE);
115	}
116	if (p->resultsH == NIL) {
117	SetNaN64(&p->result); /* return NaN if wave is not valid */
118	return(NON_EXISTENT_WAVE);
119	}
120
121	p->result = 0;
122
123	// trusting that all inputs are DOUBLE PRECISION WAVES!!!
124	inputWave = WaveData(p->inputWaveH);
125	ran_dev = WaveData(p->ran_devH);
126	nt = WaveData(p->ntH);
127	j1 = WaveData(p->j1H);
128	j2 = WaveData(p->j2H);
129	nn = WaveData(p->nnH);
130	// MC_linear_data = WaveData(p->MC_linear_dataH);
131	results = WaveData(p->resultsH);
132
133	seed = (long)results[0];
134
135	// sprintf(buf, "input seed = %ld\r", seed);
136	// XOPNotice(buf);
137
138	if(seed >= 0) {
139	seed = -1234509876;
140	}
141
142	dummy = ran1(&seed); //initialize the random sequence by passing in a negative value
143	seed = 12348765; //non-negative after that does nothing
144
145	imon = (int)inputWave[0];
146	r1 = inputWave[1];
147	r2 = inputWave[2];
148	xCtr = inputWave[3];
149	yCtr = inputWave[4];
150	sdd = inputWave[5];
151	pixSize = inputWave[6];
152	thick = inputWave[7];
153	wavelength = inputWave[8];
154	sig_incoh = inputWave[9];
155	sig_sas = inputWave[10];
156	xCtr_long = round(xCtr);
157	yCtr_long = round(yCtr);
158
159	dummy = MDGetWaveScaling(p->ran_devH, 0, &delta, &left); //0 is the rows
160	if (result = MDGetWaveDimensions(p->ran_devH, &numDimensions, dimensionSizes))
161	return result;
162	numRows_ran_dev = dimensionSizes[0];
163
164	pi = 4.0*atan(1.0);
165
166	// access the 2D wave data for writing using the direct method
167	wavH = p->MC_linear_dataH;
168	if (wavH == NIL)
169	return NOWAV;
170
171	// waveType = WaveType(wavH);
172	// if (waveType & NT_CMPLX)
173	// return NO_COMPLEX_WAVE;
174	// if (waveType==TEXT_WAVE_TYPE)
175	// return NUMERIC_ACCESS_ON_TEXT_WAVE;
176	// if (result = MDGetWaveDimensions(wavH, &numDimensions, dimensionSizes))
177	// return result;
178	// numRows = dimensionSizes[0];
179	// numColumns = dimensionSizes[1];
180
181	// if (result = MDAccessNumericWaveData(wavH, kMDWaveAccessMode0, &dataOffset))
182	// return result;
183
184	// hState = MoveLockHandle(wavH); // So wave data can't move. Remember to call HSetState when done.
185	// dataStartPtr = (char)(wavH) + dataOffset;
186	// dp0 = (double*)dataStartPtr; // Pointer to the start of the 2D wave data.
187
188	//scattering power and maximum qvalue to bin
189	// zpow = .1 //scattering power, calculated below
190	qmax = 4.0*pi/wavelength; //maximum Q to bin 1D data. (A-1) (not really used)
191	sigabs_0 = 0.0; // ignore absorption cross section/wavelength [1/(cm A)]
192	n_index = 50; // maximum number of scattering events per neutron
193	num_bins = 200; //number of 1-D bins (not really used)
194
195	//c total SAS cross-section
196	//
197	zpow = sig_sas*thick; //since I now calculate the sig_sas from the model
198	sig_abs = sigabs_0 * wavelength;
199	sig_total = sig_abs + sig_sas + sig_incoh;
200	// Print "The TOTAL XSECTION. (CM-1) is ",sig_total
201	// Print "The TOTAL SAS XSECTION. (CM-1) is ",sig_sas
202	// results[0] = sig_total;
203	// results[1] = sig_sas;
204	// RATIO = SIG_ABS / SIG_TOTAL
205	ratio = sig_incoh / sig_total;
206
207	theta_max = wavelengthqmax/(2pi);
208	//C SET Theta-STEP SIZE.
209	dth = theta_max/num_bins;
210	// Print "theta bin size = dth = ",dth
211
212	//C INITIALIZE COUNTERS.
213	n1 = 0;
214	n2 = 0;
215	n3 = 0;
216	NSingleIncoherent = 0;
217	NSingleCoherent = 0;
218	NDoubleCoherent = 0;
219	NMultipleScatter = 0;
220	NScatterEvents = 0;
221	isOn = 0;
222
223	//C MONITOR LOOP - looping over the number of incedent neutrons
224	//note that zz, is the z-position in the sample - NOT the scattering power
225	// NOW, start the loop, throwing neutrons at the sample.
226	do {
227	////SpinProcess() IS A CALLBACK, and not good for Threading!
228	// if ((n1 % 1000 == 0) && gCallSpinProcess && SpinProcess()) { // Spins cursor and allows background processing.
229	// result = -1; // User aborted.
230	// break;
231	// }
232
233	vx = 0.0; // Initialize direction vector.
234	vy = 0.0;
235	vz = 1.0;
236
237	theta = 0.0; // Initialize scattering angle.
238	phi = 0.0; // Intialize azimuthal angle.
239	n1 += 1; // Increment total number neutrons counter.
240	done = 0; // True when neutron is absorbed or when scattered out of the sample.
241	index = 0; // Set counter for number of scattering events.
242	zz = 0.0; // Set entering dimension of sample.
243	incoherentEvent = 0;
244	coherentEvent = 0;
245
246	do { // Makes sure position is within circle.
247	ran = ran1(&seed); //[0,1]
248	xx = 2.0r1(ran-0.5); //X beam position of neutron entering sample.
249	ran = ran1(&seed); //[0,1]
250	yy = 2.0r1(ran-0.5); //Y beam position ...
251	rr = sqrt(xxxx+yyyy); //Radial position of neutron in incident beam.
252	} while(rr>r1);
253
254	do { //Scattering Loop, will exit when "done" == 1
255	// keep scattering multiple times until the neutron exits the sample
256	ran = ran1(&seed); //[0,1] RANDOM NUMBER FOR DETERMINING PATH LENGTH
257	ll = path_len(ran,sig_total);
258	//Determine new scattering direction vector.
259	err = NewDirection(&vx,&vy,&vz,theta,phi); //vx,vy,vz updated, theta, phi unchanged by function
260
261	//X,Y,Z-POSITION OF SCATTERING EVENT.
262	xx += ll*vx;
263	yy += ll*vy;
264	zz += ll*vz;
265	rr = sqrt(xxxx+yyyy); //radial position of scattering event.
266
267	//sprintf(buf, "xx,yy,zz,vx,vy,vz,ll = %g %g %g %g %g %g %g\r",xx,yy,zz,vx,vy,vz,ll);
268	//XOPNotice(buf);
269
270	//Check whether interaction occurred within sample volume.
271	if (((zz > 0.0) && (zz < thick)) && (rr < r2)) {
272	//NEUTRON INTERACTED.
273	//sprintf(buf,"neutron interacted\r");
274	//XOPNotice(buf);
275
276	index += 1; //Increment counter of scattering events.
277	if (index == 1) {
278	n2 += 1; //Increment # of scat. neutrons
279	}
280	ran = ran1(&seed); //[0,1]
281	//Split neutron interactions into scattering and absorption events
282	if (ran > ratio ) { //C NEUTRON SCATTERED coherently
283	//sprintf(buf,"neutron scatters coherently\r");
284	//XOPNotice(buf);
285	coherentEvent = 1;
286	find_theta = 0; //false
287	do {
288	// pick a q-value from the deviate function
289	// pnt2x truncates the point to an integer before returning the x
290	// so get it from the wave scaling instead
291	// q0 =left + binarysearchinterp(ran_dev,ran1(seed))*delta;
292
293	q0 =left + locate_interp(ran_dev,numRows_ran_dev,ran1(&seed))*delta;
294	theta = q0/2/pi*wavelength; //SAS approximation
295
296	find_theta = 1; //always accept
297
298	//sprintf(buf, "after locate_interp call q0 = %g, theta = %g,left = %g,delta = %g\r",q0,theta,left,delta);
299	//XOPNotice(buf);
300
301	} while(!find_theta);
302
303	ran = ran1(&seed); //[0,1]
304	phi = 2.0piran; //Chooses azimuthal scattering angle.
305	} else {
306	//NEUTRON scattered incoherently
307	//sprintf(buf,"neutron scatters incoherent\r");
308	//XOPNotice(buf);
309	incoherentEvent = 1;
310	// phi and theta are random over the entire sphere of scattering
311	// !can't just choose random theta and phi, won't be random over sphere solid angle
312
313	ran = ran1(&seed); //[0,1]
314	theta = acos(2.0*ran-1);
315
316	ran = ran1(&seed); //[0,1]
317	phi = 2.0piran; //Chooses azimuthal scattering angle.
318	} //(ran > ratio)
319	} else {
320	//NEUTRON ESCAPES FROM SAMPLE -- bin it somewhere
321	done = 1; //done = true, will exit from loop
322	//Increment #scattering events array
323	MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
324	indices[0] =index; //this sets access to nn[index]
325	if (index <= n_index) {
326	if (result = MDGetNumericWavePointValue(p->nnH, indices, value))
327	return result;
328	value[0] += 1; // add one to the value
329	if (result = MDSetNumericWavePointValue(p->nnH, indices, value))
330	return result;
331	// nn[index] += 1;
332	}
333
334	if( index != 0) { //neutron was scattered, figure out where it went
335	theta_z = acos(vz); // Angle (= 2theta) WITH respect to z axis.
336	testQ = 2pisin(theta_z)/wavelength;
337
338	// pick a random phi angle, and see if it lands on the detector
339	// since the scattering is isotropic, I can safely pick a new, random value
340	// this would not be true if simulating anisotropic scattering.
341	testPhi = ran1(&seed)2pi;
342
343	// is it on the detector?
344	FindPixel(testQ,testPhi,wavelength,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel);
345
346	if(xPixel != -1 && yPixel != -1) {
347	isOn += 1;
348	MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
349	indices[0] = xPixel;
350	indices[1] = yPixel;
351	if (result = MDGetNumericWavePointValue(wavH, indices, value))
352	return result;
353	value[0] += 1; // Real part
354	if (result = MDSetNumericWavePointValue(wavH, indices, value))
355	return result;
356	//if(index==1) // only the single scattering events
357	//dp = dp0 + xPixel + yPixel*numColumns; //offset the pointer to the exact memory location
358	//*dp += 1; //increment the value there
359	//endif
360	}
361
362
363	/* is this causing me a problem since I'm not locking these? Probably not, since it crashes even if I comment these out... */
364	if(theta_z < theta_max) {
365	//Choose index for scattering angle array.
366	//IND = NINT(THETA_z/DTH + 0.4999999)
367	ind = round(theta_z/dth + 0.4999999); //round is eqivalent to nint()
368	nt[ind] += 1; //Increment bin for angle.
369	//Increment angle array for single scattering events.
370	if (index == 1) {
371	j1[ind] += 1;
372	}
373	//Increment angle array for double scattering events.
374	if (index == 2) {
375	j2[ind] += 1;
376	}
377	}
378	/**/
379
380	// increment all of the counters now since done==1 here and I'm sure to exit and get another neutron
381	NScatterEvents += index; //total number of scattering events
382	if(index == 1 && incoherentEvent == 1) {
383	NSingleIncoherent += 1;
384	}
385	if(index == 1 && coherentEvent == 1) {
386	NSingleCoherent += 1;
387	}
388	if(index == 2 && coherentEvent == 1 && incoherentEvent == 0) {
389	NDoubleCoherent += 1;
390	}
391	if(index > 1) {
392	NMultipleScatter += 1;
393	}
394	} else { // index was zero, neutron must be transmitted, so just increment the proper counters and data
395	isOn += 1;
396	nt[0] += 1;
397	MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
398	indices[0] = xCtr_long;
399	indices[1] = yCtr_long;
400	if (result = MDGetNumericWavePointValue(wavH, indices, value))
401	return result;
402	value[0] += 1; // Real part
403	if (result = MDSetNumericWavePointValue(wavH, indices, value))
404	return result;
405	}
406	}
407	} while (!done);
408	} while(n1 < imon);
409
410	// assign the results to the wave
411
412	MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
413	value[0] = (double)n1;
414	indices[0] = 0;
415	if (result = MDSetNumericWavePointValue(p->resultsH, indices, value))
416	return result;
417	value[0] = (double)n2;
418	indices[0] = 1;
419	if (result = MDSetNumericWavePointValue(p->resultsH, indices, value))
420	return result;
421	value[0] = (double)isOn;
422	indices[0] = 2;
423	if (result = MDSetNumericWavePointValue(p->resultsH, indices, value))
424	return result;
425	value[0] = (double)NScatterEvents;
426	indices[0] = 3;
427	if (result = MDSetNumericWavePointValue(p->resultsH, indices, value))
428	return result;
429	value[0] = (double)NSingleCoherent;
430	indices[0] = 4;
431	if (result = MDSetNumericWavePointValue(p->resultsH, indices, value))
432	return result;
433	value[0] = (double)NDoubleCoherent;
434	indices[0] = 5;
435	if (result = MDSetNumericWavePointValue(p->resultsH, indices, value))
436	return result;
437	value[0] = (double)NMultipleScatter;
438	indices[0] = 6;
439	if (result = MDSetNumericWavePointValue(p->resultsH, indices, value))
440	return result;
441
442
443	// HSetState((Handle)wavH, hState); //release the handle of the 2D data wave
444	// WaveHandleModified(wavH); // Inform Igor that we have changed the wave. (CALLBACK! needed, but not allowed in Threading)
445
446	return(0);
447	}
448	//////// end of main function for calculating multiple scattering
449

Note: See TracBrowser for help on using the repository browser.

Download in other formats: