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

Last change on this file since 458 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
14static 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
45int
46Monte_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 = wavelength*qmax/(2*pi);
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.0*r1*(ran-0.5);          //X beam position of neutron entering sample.
249                        ran = ran1(&seed);              //[0,1]
250                        yy = 2.0*r1*(ran-0.5);          //Y beam position ...
251                        rr = sqrt(xx*xx+yy*yy);         //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(xx*xx+yy*yy);         //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.0*pi*ran;                       //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.0*pi*ran;                       //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 = 2*pi*sin(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)*2*pi;
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.