source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Alpha/Tinker/FFT_Debye_Spheres.ipf @ 798

Last change on this file since 798 was 798, checked in by srkline, 11 years ago

I'm not sure this is a great idea, but I'm putting the FFT / Debye sphere work that I have completed into SVN. It's all really rough - the math, I believe is correct, but the interface if really, really rough. But it's not going to develop without help.

File size: 29.9 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2
3
4// with the Debye Sphere method -- as of Jan 2011
5// the discretization seems to work now, and is significantly faster than the old double sum
6//
7// now I need to XOP-it (done)
8
9
10// routines to calculate I(q) and p(r) from the 3d arrangement of spheres
11// using Debye's method
12//
13
14
15
16
17//calculate I(q) given a set of spheres as xyz
18//
19// this is the REALLY slow way
20//
21Function CalcIQRfromXYZ(xv,yv,zv,qval,rho,rval,grid)
22        Wave xv,yv,zv
23        Variable qval,rho,rval,grid
24       
25        Variable num=numpnts(xv),ii,kk,dik
26        Variable Iqr,vol,dum,fQR
27       
28        //NVAR dCall=dCall
29        vol=4*Pi/3*rval*rval*rval
30        Iqr=0
31        fQR=PhiQR(qval,rval)
32        //do i=j sum
33        for(ii=0;ii<num;ii+=1)
34                dum = rho*vol*fQR
35                Iqr += dum*dum
36        endfor
37        //do i != j (double) sum
38        for(ii=0;ii<num;ii+=1)
39                for(kk=(ii+1);kk<num;kk+=1)
40                        dik=AV_Distance(xv[ii],xv[kk],yv[ii],yv[kk],zv[ii],zv[kk]) * grid
41                        //dCall += 1
42                        dum = rho*vol*fQR
43                        Iqr += 2*dum*dum*sinc(dik*qval)
44                endfor
45        endfor
46               
47        return(Iqr)
48End
49
50//calculate I(q) given a set of spheres as xyz
51//
52// this is the REALLY slow way
53//
54Function CalcIQRfromMat(mat,qval,rho,rval,grid)
55        Wave mat
56        Variable qval,rho,rval,grid
57       
58        Variable num,ii,kk,dik
59        Variable Iqr,vol,dum,fQR,dCall=0
60       
61        ParseMatrix3D_rho(mat)
62        WAVE xv=x3d
63        WAVE yv=y3d
64        WAVE zv=z3d
65        WAVE rho3d=rho3d
66        num=numpnts(xv)
67       
68        vol=4*Pi/3*rval*rval*rval
69        Iqr=0
70        fQR=PhiQR(qval,rval)
71        //do i=j sum
72        for(ii=0;ii<num;ii+=1)
73                dum = rho3d[ii]*vol*fQR
74                Iqr += dum*dum
75        endfor
76               
77        //do i != j (double) sum
78        for(ii=0;ii<num;ii+=1)
79                for(kk=(ii+1);kk<num;kk+=1)
80                        dik=AV_Distance(xv[ii],xv[kk],yv[ii],yv[kk],zv[ii],zv[kk]) * grid
81                        //dCall += 1
82                        dum = vol*fQR
83                        Iqr += 2*dum*dum*rho3d[ii]*rho3d[kk]*sinc(dik*qval)
84                endfor
85        endfor
86       
87//      Print "num, dCall = ", num,dCall
88        return(Iqr)
89End
90
91//calculate I(q) given a set of spheres as xyz
92//
93// this is the improved way, doing distance binning
94// -- substantially speeded up with two XOPs that take care of the double loops
95// binning into the histogram is the bottleneck at this point, but I see no
96// easy way to make this multi-processor aware. Supposedly it should be easy
97// to do, but I don't see the fast way to dispatch it to threads without writing four
98// functions like in the MC simulation code. Maybe data folders?
99//
100// some of the work is still done in Igor, but it's the really fast stuff, after checking the timing of each step.
101//
102// can't use the SLD information here, since one of the major assumptions of this simplification
103// is that the SLDs of the spheres are all the same
104Function CalcIQRfromMat_bin(qW,iqr,xv,yv,zv,rho,rval,grid)
105        Wave qW,iqr,xv,yv,zv
106        Variable rho,rval,grid
107       
108        Variable num,ii,kk,dik,iter,t1
109        Variable vol,dum,fQR,qval,F2Q
110        Variable dmax,binWidth,Qmax,numBins,binIndex,val
111        Variable nthreads
112
113        t1 = ticks
114       
115        num=numpnts(xv)
116       
117        //find maximum distance
118       
119//// I don't have enough space for the NxN matrix once there are 10k or so points.
120////    Make/O/D/N=(num,num) distMat
121////    distMat = 0
122//
123//      iter = 0
124//      for(ii=0;ii<num;ii+=1)
125//              for(kk=(ii+1);kk<num;kk+=1)
126////                    dist[ii]=AV_Distance(xv[ii],xv[kk],yv[ii],yv[kk],zv[ii],zv[kk]) * grid
127//                      distMat[ii][kk]=AV_Distance(xv[ii],xv[kk],yv[ii],yv[kk],zv[ii],zv[kk]) * grid
128//              endfor
129//      endfor
130
131//      // so this is a quick way to find the maximum distance. entering a value would be easier
132//      Make/O/D/N=(num-1) dist
133//      for(ii=0;ii<num;ii+=1)
134//              for(kk=(ii+1);kk<num;kk+=1)
135//                      dist[ii]=AV_Distance(xv[ii],xv[kk],yv[ii],yv[kk],zv[ii],zv[kk]) * grid
136//              endfor
137//      endfor
138//     
139//      dmax = WaveMax(dist)
140//      Print "dmax = ",dmax
141//      Printf "Finding dmax Igor = %g seconds\r",(ticks-t1)/60.15
142//      t1 = ticks
143       
144        ////// Find maxiumum distance
145        // write an XOP that takes the xyz and returns the square of the maximum distance, must sqrt then multiply by grid distance
146       
147        // the  num < 10000 points is an empirical value for where the overhead of threading is OK
148        nthreads=ThreadProcessorCount
149        if(nthreads == 1 || num < 10000 )               
150                dmax = maxDistanceX(xv,yv,zv,0,numpnts(xv))
151//              Print "dmax = ",dmax
152//              Printf "Finding dmax XOP = %g seconds\r",(ticks-t1)/60.15
153//              t1 = ticks
154        else
155                dmax = maxDistance_Threaded(xv,yv,zv)
156//              Print "dmax = ",dmax
157        endif
158       
159        dmax = sqrt(dmax)
160        dmax *= grid
161//      Print "dmax = ",dmax
162
163//      Printf "Finding dmax Threaded = %g seconds\r",(ticks-t1)/60.15
164//      t1 = ticks
165
166       
167        // use a bin width Dmax / 10000 as suggested in Otto's book, pg 160
168        //
169        Variable ndiv=100000
170        binWidth = dmax/ndiv
171//      Print "binWidth = ",binWidth
172        numBins = ndiv
173        Make/O/D/N=(numBins) distBins
174       
175        distBins = 0
176       
177//      Printf "Initializing distBins = %g seconds\r",(ticks-t1)/60.15
178//      t1 = ticks
179       
180//      for(ii=0;ii<num;ii+=1)
181//              for(kk=(ii+1);kk<num;kk+=1)
182//                      val = AV_Distance(xv[ii],xv[kk],yv[ii],yv[kk],zv[ii],zv[kk]) * grid
183//                      binIndex = trunc(val/binWidth-0.5)
184//                      if(binIndex > numBins -1 )
185//                              Print "bad index"
186//                      else
187//                              distBins[binIndex] += 1
188//                      endif
189//              endfor
190//      endfor
191
192
193///// threading of the bins
194
195        //write an XOP that takes the xyz and bins, and fills them
196        //distBins is returned
197        nthreads=ThreadProcessorCount
198
199        if(nthreads == 1)
200                binDistanceX(xv, yv, zv, distBins, grid, binWidth,0,numpnts(xv))                //the end point is the numpnts, not n-1
201        else
202                binDistance_Threaded(xv, yv, zv, distBins, grid, binWidth)
203        endif
204
205//      Printf "Binning = %g seconds\r",(ticks-t1)/60.15
206//      t1 = ticks
207       
208//      WaveStats/Q distBins
209//      Print "binned values, numSpheres = ",V_avg*V_npnts,num
210       
211        Duplicate/O distBins dist_at_bin
212        for(ii=0;ii<numBins;ii+=1)
213                dist_at_bin[ii] =  ii*binWidth
214        endfor
215       
216       
217        // remove all of the bins with zero distance, then reset the number of Bins
218        RemoveZerosXY(distBins, dist_at_bin)
219        numBins = numpnts(distBins)
220       
221//      Printf "Removing zeroes = %g seconds\r",(ticks-t1)/60.15
222//      t1 = ticks
223       
224        vol=4*Pi/3*rval*rval*rval
225       
226        Iqr=0
227       
228        Variable nq = numpnts(qW)
229
230//      Make/O/D/N=(nq) Sij_Bin
231       
232        for(kk=0;kk<nq;kk+=1)           // loop over the q-values
233       
234                qval = qW[kk]
235                fQR=PhiQR(qval,rval)
236                F2Q = fQR*fQR*rho*rho*vol*vol
237               
238                for(ii=0;ii<numBins;ii+=1)             
239                        Iqr[kk] +=  distBins[ii]*sinc(qval*dist_at_bin[ii])
240                endfor
241               
242//              Sij_Bin[kk] = num + 2*Iqr[kk]
243               
244                Iqr[kk] = F2Q*(num + 2*Iqr[kk])
245                       
246        endfor
247       
248//      Printf "Loop over q = %g seconds\r",(ticks-t1)/60.15
249       
250        return(0)
251End
252
253// support up to 8 threads at this time
254Function binDistance_Threaded(xv, yv, zv, distBins, grid, binWidth)
255        Wave xv, yv, zv, distBins
256        Variable grid, binWidth
257       
258        Variable nthreads,mt,left,right
259        Variable ii,num
260
261        nthreads=ThreadProcessorCount
262
263        if(nthreads > 8)
264                nthreads = 8
265        endif
266       
267        mt = ThreadGroupCreate(nthreads)
268       
269        ii=0
270        left = 0
271        num = numpnts(xv)
272       
273        for(ii=0;ii<nthreads;ii+=1)
274                Duplicate/O distBins $("distBins"+num2str(ii))
275               
276//  this is an even spreading of the points - not appropriate here for the triangle
277//                      Print (ii*num/nthreads),((ii+1)*num/nthreads)
278//                      ThreadStart mt,i,Cyl_PolyRadius_T(cw,yw,xw,(ii*num/nthreads),((ii+1)*num/nthreads))
279
280// this splits up the triangle into equal area chunks trapezoid left(ii) -> right(ii)
281// be sure that for the last iteration of ii, right = 1 *(num)
282                right = 1 - sqrt(1- (ii+1)/nthreads)
283
284//                      Print left,right
285               
286                if(ii==0)
287                        Wave distBins0
288                        distBins0 = 0
289                        //Print (ii*num/nthreads),((ii+1)*num/nthreads)
290                        //ThreadStart mt,ii,binDistance_WF(xv, yv, zv, distBins0, grid, binWidth,(ii*num/nthreads),((ii+1)*num/nthreads))
291                        ThreadStart mt,ii,binDistance_WF(xv, yv, zv, distBins0, grid, binWidth, left*num, right*num)
292                endif
293                if(ii==1)
294                        Wave distBins1
295                        distBins1 = 0
296                        ThreadStart mt,ii,binDistance_WF(xv, yv, zv, distBins1, grid, binWidth, left*num, right*num)
297                endif
298                if(ii==2)
299                        Wave distBins2
300                        distBins2 = 0
301                        ThreadStart mt,ii,binDistance_WF(xv, yv, zv, distBins2, grid, binWidth, left*num, right*num)
302                endif
303                if(ii==3)
304                        Wave distBins3
305                        distBins3 = 0
306                        ThreadStart mt,ii,binDistance_WF(xv, yv, zv, distBins3, grid, binWidth, left*num, right*num)
307                endif
308                if(ii==4)
309                        Wave distBins4
310                        distBins4 = 0
311                        ThreadStart mt,ii,binDistance_WF(xv, yv, zv, distBins4, grid, binWidth, left*num, right*num)
312                endif
313                if(ii==5)
314                        Wave distBins5
315                        distBins5 = 0
316                        ThreadStart mt,ii,binDistance_WF(xv, yv, zv, distBins5, grid, binWidth, left*num, right*num)
317                endif
318                if(ii==6)
319                        Wave distBins6
320                        distBins6 = 0
321                        ThreadStart mt,ii,binDistance_WF(xv, yv, zv, distBins6, grid, binWidth, left*num, right*num)
322                endif
323                if(ii==7)
324                        Wave distBins7
325                        distBins7 = 0
326                        ThreadStart mt,ii,binDistance_WF(xv, yv, zv, distBins7, grid, binWidth, left*num, right*num)
327                endif
328               
329                left = right
330
331        endfor
332       
333        // wait until done
334        do
335                variable tgs= ThreadGroupWait(mt,100)
336        while( tgs != 0 )
337        variable dummy= ThreadGroupRelease(mt)
338        mt=0
339//      Print "done with all threads"
340
341        //      then add them all back together
342        if(nthreads == 1)
343                distBins = distBins0            // add up each instance
344        endif
345        if(nthreads == 2)
346                distBins = distBins0    + distBins1
347        endif
348        if(nthreads == 3)
349                distBins = distBins0    + distBins1 + distBins2
350        endif
351        if(nthreads == 4)
352                distBins = distBins0    + distBins1 + distBins2 + distBins3
353        endif
354        if(nthreads == 5)
355                distBins = distBins0    + distBins1 + distBins2 + distBins3 + distBins4
356        endif
357        if(nthreads == 6)
358                distBins = distBins0    + distBins1 + distBins2 + distBins3 + distBins4 + distBins5
359        endif
360        if(nthreads == 7)
361                distBins = distBins0    + distBins1 + distBins2 + distBins3 + distBins4 + distBins5 + distBins6
362        endif
363        if(nthreads == 8)
364                distBins = distBins0    + distBins1 + distBins2 + distBins3 + distBins4 + distBins5 + distBins6 + distBins7
365        endif
366       
367        // then clean up
368        KillWaves/Z distBins0,distBins1,distBins2,distBins3
369        KillWaves/Z distBins4,distBins5,distBins6,distBins7
370       
371        return(0)
372End
373
374// this is just a worker function to get the ThreadStart operation to compile
375ThreadSafe Function binDistance_WF(xv, yv, zv, distBins, grid, binWidth,p1,p2)
376        Wave xv, yv, zv, distBins
377        Variable grid, binWidth,p1,p2
378
379        Variable ret
380        ret = binDistanceX(xv, yv, zv, distBins, grid, binWidth,p1,p2)
381       
382        return(0)
383End
384
385
386// support up to 8 threads at this time
387Function maxDistance_Threaded(xv, yv, zv)
388        Wave xv, yv, zv
389       
390        Variable nthreads,mt,left,right
391        Variable ii,num
392
393        nthreads=ThreadProcessorCount
394
395        if(nthreads > 8)
396                nthreads = 8
397        endif
398       
399        mt = ThreadGroupCreate(nthreads)
400       
401        ii=0
402        left = 0
403        num = numpnts(xv)
404       
405        for(ii=0;ii<nthreads;ii+=1)
406               
407//  this is an even spreading of the points - not appropriate here for the triangle
408//                      Print (ii*num/nthreads),((ii+1)*num/nthreads)
409//                      ThreadStart mt,i,Cyl_PolyRadius_T(cw,yw,xw,(ii*num/nthreads),((ii+1)*num/nthreads))
410
411// this splits up the triangle into equal area chunks trapezoid left(ii) -> right(ii)
412// be sure that for the last iteration of ii, right = 1 *(num)
413                right = 1 - sqrt(1- (ii+1)/nthreads)
414
415//                      Print left,right
416               
417                ThreadStart mt,ii,maxDistance_WF(xv, yv, zv, left*num, right*num)
418               
419                left = right
420
421        endfor
422       
423        // wait until done
424        do
425                variable tgs= ThreadGroupWait(mt,100)
426        while( tgs != 0 )
427
428
429//      Get the return values, and find the maximum to return
430        Variable maxValue = 0
431       
432        for(ii=0;ii<nthreads;ii+=1)
433//                      Print ThreadReturnValue(mt,ii)
434                        maxValue = max(maxValue,ThreadReturnValue(mt,ii))
435        endfor
436
437        // now release the threads
438        variable dummy= ThreadGroupRelease(mt)
439        mt=0
440//      Print "done with all threads"
441       
442        return(maxValue)
443End
444
445// this is just a worker function to get the ThreadStart operation to compile
446ThreadSafe Function maxDistance_WF(xv, yv, zv, p1, p2)
447        Wave xv, yv, zv
448        Variable p1,p2
449
450        Variable ret
451        ret = maxDistanceX(xv,yv,zv,p1,p2)
452       
453        return(ret)
454End
455
456// distance separating two xyz points
457// (this is one of the big time-consuming steps)
458Function AV_Distance(x1,x2,y1,y2,z1,z2)
459        Variable x1,x2,y1,y2,z1,z2
460       
461        Variable retval
462        retval=sqrt( (x1-x2)^2 + (y1-y2)^2 + (z1-z2)^2 )
463        return(retval)
464End
465
466//bessel function...
467
468Function PhiQR(qval,rval)
469        Variable qval,rval
470       
471        Variable retval,qr
472       
473        qr=qval*rval
474        retval=(sin(qr)-qr*cos(qr))/qr/qr/qr
475       
476        return(3*retval)
477End
478
479
480
481//called from the FFT method panel
482//
483// calls the full double sum that can take the SLD wave into account, if filled
484// - as a result, it's slower than the binned calculation
485//
486Proc DoSpheresCalcFFTPanel(num,qMin,qMax)
487        Variable num=100,qmin=0.004,qmax=0.5
488       
489        Variable t1
490        String qStr="qval_full",iStr="ival_full"                //default wave names, always overwritten
491        Variable grid
492       
493        grid=root:FFT_T
494
495        Make/O/D/N=(num) $qStr,$iStr
496        $qStr = alog(log(qmin) + x*((log(qmax)-log(qmin))/num))         
497       
498        Variable estTime,nx
499        String str = ""
500
501        nx = NonZeroValues(mat)
502       
503        estTime = EstimatedTime(nx,num,0)               // 0 =  XOP, 1 = no XOP
504        sprintf str, "Estimated time for the calculation is %g seconds. Proceed?",estTime
505        DoAlert 1,str
506        if(V_Flag==1)           //yes, proceed
507                t1=ticks
508                fDoCalc($qStr,$iStr,grid,0,1)
509//              Printf "Elapsed AltiSpheres time = %g seconds\r\r",(ticks-t1)/60.15
510        Endif
511End
512
513
514//called from the FFT method panel
515//
516// in this method, the distances are binned as by Otto Glatter, and has been partially XOPed
517//
518Proc DoBinnedSpheresCalcFFTPanel(num,qMin,qMax)
519        Variable num=100,qmin=0.004,qmax=0.5
520       
521        Variable t1
522        String qStr="qval_XOP",iStr="ival_XOP"          //default wave names, always overwritten
523        Variable grid
524       
525        grid=root:FFT_T
526
527        Make/O/D/N=(num) $qStr,$iStr
528        $qStr = alog(log(qmin) + x*((log(qmax)-log(qmin))/num))         
529       
530        Variable estTime,nx
531        String str = ""
532
533        nx = NonZeroValues(mat)
534       
535        estTime = EstimatedTime(nx,num,2)               // 0 =  XOP, 1 = no XOP, 2 = binned distances
536        sprintf str, "Estimated time for the calculation is %g seconds. Proceed?",estTime
537        DoAlert 1,str
538        if(V_Flag==1)           //yes, proceed
539                t1=ticks
540                fDoCalc($qStr,$iStr,grid,2,1)
541//              Printf "Elapsed AltiSpheres time = %g seconds\r\r",(ticks-t1)/60.15
542        Endif
543End
544
545//called from the FFT method panel
546//
547// in this method, the distances are binned as by Otto Glatter, and has been partially XOPed
548//
549Proc DoBinnedSLDCalcFFTPanel(num,qMin,qMax)
550        Variable num=100,qmin=0.004,qmax=0.5
551       
552        Variable t1
553        String qStr="qval_SLD",iStr="ival_SLD"          //default wave names, always overwritten
554        Variable grid
555       
556        grid=root:FFT_T
557
558        Make/O/D/N=(num) $qStr,$iStr
559        $qStr = alog(log(qmin) + x*((log(qmax)-log(qmin))/num))         
560       
561        Variable estTime,nx
562        String str = ""
563
564        nx = NonZeroValues(mat)
565       
566        estTime = EstimatedTime(nx,num,3)               // 0 =  XOP, 1 = no XOP, 2 = binned distances
567        sprintf str, "Estimated time for the calculation is %g seconds. Proceed?",estTime
568        DoAlert 1,str
569        if(V_Flag==1)           //yes, proceed
570                t1=ticks
571                fDoCalc($qStr,$iStr,grid,3,1)
572//              Printf "Elapsed AltiSpheres time = %g seconds\r\r",(ticks-t1)/60.15
573        Endif
574End
575
576
577
578////////////////////
579
580//waves must exist at this point
581// ival_(ext), qval_(ext) are what is passed at this point
582//
583// switch the type of calculation based on a flag
584//
585// flag = 0 = full calculation, using double sum
586// flag = 2 = binned distances
587//                      = 3 = binned distances with SLD
588//              = 12 = Ken's generation of XYZ, skip my matrix filling and parsing
589//
590//
591Function fDoCalc(qval,ival,grid,flag,verbose)
592        Wave qval,ival
593        Variable grid,flag,verbose
594
595        Variable t1,val
596        WAVE mat = root:mat
597        // convert the matrix to XYZ triplet for calculation
598        t1=ticks
599       
600        if(flag == 12 || flag == 13)                    //skip the parsing
601        // when using Ken's method, there is no matrix to parse, or I can put
602        // his xyz (non-integer) onto the mat, then parse it back off (works fine, but not smart)       
603                Duplicate/O root:xOutW,x3d
604                Duplicate/O root:yOutW,y3d
605                Duplicate/O root:zOutW,z3d
606                Duplicate/O root:sldW,rho3d
607                flag -= 10              //switch to the binned method, either 2 or 3
608        else
609                ParseMatrix3D_rho(mat)
610                if(verbose)
611                        Printf "ParseMatrix3D time = %g seconds\r",(ticks-t1)/60.15
612                endif
613        endif
614       
615        WAVE x3d=x3d
616        WAVE y3d=y3d
617        WAVE z3d=z3d
618        WAVE rho3d=rho3d
619       
620        Variable num,nx
621        num=numpnts(qval)
622        nx=numpnts(x3d)
623       
624        if(verbose)
625                Printf "Estimated  calculation time = %g seconds for %d qvals and %d spheres\r",        EstimatedTime(nx,num,flag),num,nx
626        endif
627        t1=ticks
628       
629        NVAR Tscale=root:FFT_T
630        NVAR solventSLD = root:FFT_SolventSLD
631       
632        Variable vol,rval,frac
633        rval = grid*0.62
634        NVAR delRho = root:FFT_delRho
635        vol = nx*(Tscale)^3
636       
637        frac = nx/DimSize(mat,0)^3
638       
639        switch(flag)    // numeric switch
640                case 0:         // full double sum
641                        Duplicate/O rho3D tmpRho3d
642                        tmpRho3d = rho3d - solventSLD
643                        MultiThread     ival = DebyeSpheresX(qval,x3d,y3d,z3d,tmpRho3d,0.62*grid,grid)
644                        ival *= delRho*delRho
645                        ival /= vol
646                        ival *= 1e8
647       
648                        ival *= frac
649                       
650                        KillWaves/Z tmpRho3d
651                        break                                   
652                case 2:         // binned distances
653                                        //////improved version uses distance binning, and is now AAO
654                        // don't need to duplicate the wave since I'm passing only a single value for the SLD
655                        CalcIQRfromMat_bin(qval,ival,x3d,y3d,z3d,rho3D[0]-solventSLD, rval, grid)
656                        ival *= delRho*delRho
657                        ival /= vol
658                        ival *= 1e8
659                       
660                        ival *= frac
661                       
662                        break
663                case 3:         // binned distances + SLD
664                                        //////improved version uses distance binning, and is now AAO
665                        Duplicate/O rho3D tmpRho3d
666                        tmpRho3d = rho3d - solventSLD
667                       
668                        CalcIQRfromMat_bin_SLD(qval,ival,x3d,y3d,z3d,tmpRho3d, rval, grid)
669                        ival *= delRho*delRho
670                        ival /= vol
671                        ival *= 1e8
672                       
673                        ival *= frac
674                       
675                        KillWaves/Z tmpRho3d
676                        break
677                default:                                                        // optional default expression executed
678                        Print "No case matched in fDoCalc"                              // when no case matches
679        endswitch
680       
681        if(verbose)
682                Printf "Elapsed time for %d qvals and %d spheres = %g seconds\r",num,nx,(ticks-t1)/60.15
683        endif
684//      Duplicate/O ival_XOP ival_Igor,ival_binned
685       
686//////
687//      Printf "Estimated Igor Function time = %g seconds for %d qvals and %d spheres\r",       EstimatedTime(nx,num,1),num,nx
688//      t1=ticks
689//      ival_Igor = CalcIQRfromMat(mat,qval,1, rval, grid)
690//     
691////    vol = nx*4*Pi/3*rval*rval*rval
692//      ival_Igor *= delRho*delRho
693//      ival_Igor /= vol
694//      ival_Igor *= 1e8
695//     
696//      Printf "Igor function time for %d qvals and %d spheres = %g seconds\r",num,nx,(ticks-t1)/60.15
697///////
698        return(0)
699End
700
701
702
703// Based on the WM procedure
704// RemoveNaNsXY(theXWave, theYWave)
705//      Removes all Zeroes in an XY pair if the X wave has the value 0.
706//      Returns the number of points removed.
707Function RemoveZerosXY(theXWave, theYWave)
708        Wave theXWave
709        Wave theYWave
710
711        Variable p, numPoints, numNaNs
712        Variable xval, yval
713       
714        numNaNs = 0
715        p = 0                                                                                   // the loop index
716        numPoints = numpnts(theXWave)                   // number of times to loop
717
718        do
719                xval = theXWave[p]
720                yval = theYWave[p]
721                if (xval == 0)          //is xWave == 0?
722                        numNaNs += 1
723                else                                                                            // if not an outlier
724                        theXWave[p - numNaNs] = xval            // copy to input wave
725                        theYWave[p - numNaNs] = yval            // copy to input wave
726                endif
727                p += 1
728        while (p < numPoints)
729       
730        // Truncate the wave
731        DeletePoints numPoints-numNaNs, numNaNs, theXWave, theYWave
732       
733        return(numNaNs)
734End
735
736
737///////// binned Debye, with multiple SLDs
738// as in E. Pantos et.al. J. Molec. Struct. 383 (1996) 303-308.
739//
740//calculate I(q) given a set of spheres as xyz
741//
742// this is the improved way, doing distance binning
743// - in this method, there can be more than one SLD (different from the solvent)
744//
745// -(done)- substantially speeded up with two XOPs that take care of the double loops
746// binning into the histogram is the bottleneck at this point, but I see no
747// easy way to make this multi-processor aware. Supposedly it should be easy
748// to do, but I don't see the fast way to dispatch it to threads without writing four
749// functions like in the MC simulation code. Maybe data folders?
750//
751// some of the work is still done in Igor, but it's the really fast stuff, after checking the timing of each step.
752//
753//----
754// The bookkeeping here is a bit of a mess to keep track of the Gij binning, so here are the details. Note
755// that as ugly and as many loops there are in the function, the time spent in these loops
756// is practically nothing (much less than a second) compared to the binning.
757//
758// SLD_id : a vector with the values of the SLD that are not solvent. The point value of this wave
759//                              is used to refer to the "material" with the given SLD, such that SLD_id[index] = SLD_Value
760// SLDLookup :  vector that does the reverse correlation to pair up the SLD with its index, such that
761//                                              SLDLookup[SLD_Value] = index
762// psf_id_vec : a vector of the ij pairs, numbered starting from 11, 12, 13, 22, 23, 33, etc.
763// psf_id_mat : a matrix such that psf_id_mat[i][j] gives the INDEX of the (i+1)(j+1) PSF pair. For example,
764//                                              psf_id_mat[0][2] corresponds to PSF 13, with and index (in psf_id_vec) of 2.
765// Nij          :       a vector of the number of Ni type spheres. Index matches the psf_id_vec. Values are zero
766//                                      for i != j, and are the number of spheres for i==j
767//
768// binMatrix : a matrix of /N=(numBins,NumDiffPSF) that contains the Gij binned pairwise based on the SLD
769//                              the columns (numDiffPSF) are indexed to the psf_id_vec
770//
771// The Gijx and Sijx waves are calculated and saved as a function of r or q, where x is the numerical
772// index that corresponds to psf_id_vec. Note that the Sijx are not exactly what you expect for a partial
773// structure factor. Some conversion is necessary and the details are still to be worked out, but this is
774// all within Igor, and not in the XOP.
775//
776Function CalcIQRfromMat_bin_SLD(qW,iqr,xv,yv,zv,rho,rval,grid)
777        Wave qW,iqr,xv,yv,zv,rho
778        Variable rval,grid
779       
780        Variable num,ii,jj,kk,dik,iter,t1,np,index,nq
781        Variable vol,dum,fQR,qval,F2Q
782        Variable dmax,binWidth,Qmax,numBins,binIndex,val
783        Variable rhoi,rhok,newNumBins
784        Variable nthreads
785
786        NVAR    solventSLD = root:FFT_SolventSLD
787
788//      t1 = ticks
789       
790        num=numpnts(xv)
791       
792//      Printf "\r\rParsing = %g seconds\r",(ticks-t1)/60.15
793//      t1 = ticks
794               
795        ////// Find maxiumum distance
796        // write an XOP that takes the xyz and returns the square of the maximum distance, must sqrt then multiply by grid distance
797        nthreads=ThreadProcessorCount
798        if(nthreads == 1)
799                dmax = maxDistanceX(xv,yv,zv,0,numpnts(xv))
800        else
801                dmax = maxDistance_Threaded(xv,yv,zv)
802        endif
803       
804//      dmax = maxDistanceX(xv,yv,zv,0,numpnts(xv))
805        dmax = sqrt(dmax)
806        dmax *= grid
807//      Print "dmax = ",dmax
808
809//      Printf "Finding dmax XOP = %g seconds\r",(ticks-t1)/60.15
810//      t1 = ticks
811
812// Now since there are multiple SLDs, we need:
813// - a count of the number of different SLDs
814// - the number of different PSFs
815// - make the histogram matrix, keeping a lookup table of the psf_id
816
817// how many different SLDs are there?
818        //loop through the rho wave, disregarding the solvent
819        Duplicate/O rho tmpSLD
820        Make/O/D/N=0 SLD_id
821        SLD_id = 0
822        np = 0
823        do
824                val = WaveMin(tmpSLD)
825                if(val != solventSLD && val != 1e10)
826                        np += 1
827                        InsertPoints numpnts(SLD_id),1,SLD_id
828                        SLD_id[np-1] = val
829                endif
830                tmpSLD = (tmpSLD == val) ? 1e10 : tmpSLD                // set the values to 1e10
831        while(val != 1e10)
832        KillWaves/Z tmpSLD     
833
834        Variable numDiffSLD = numpnts(SLD_id)
835        Variable numDiffPSF, PSFIndex
836        numDiffPSF = NumDiffSLD*(numDiffSLD + 1)/2
837
838// make a lookup table to get from rho@xyz to the index
839        Variable maxSLD = WaveMax(SLD_id)
840        Make/O/D/N=(maxSLD+1) SLDlookup
841        SLDLookup = NaN
842        ii=0
843        do
844                SLDLookup[SLD_id[ii]] = ii
845                ii+=1
846        while(ii<numpnts(SLD_id))
847
848//      Printf "Find number of different SLDs = %g seconds\r",(ticks-t1)/60.15
849//      t1 = ticks
850
851        // use a bin width Dmax / 10000 as suggested in Otto's book, pg 160
852        //
853        numBins=100000
854        binWidth = dmax/numBins
855       
856        Make/O/D/N=(numBins,NumDiffPSF) binMatrix
857        binMatrix=0
858       
859        Make/O/D/N=(numDiffSLD,numDiffSLD) psf_id_mat
860        Make/O/D/N=(numDiffPSF) psf_id_vec,Nij
861        psf_id_mat = 0
862        psf_id_vec = 0
863        Nij = 0
864       
865        index = 0
866        for(ii=1;ii<=numDiffSLD;ii+=1)
867                for(kk=ii;kk<=numDiffSLD;kk+=1)         
868                        //Print ii,kk   
869                        psf_id_vec[index] = 10*ii+kk
870                        psf_id_mat[ii-1][kk-1] = index
871                        psf_id_mat[kk-1][ii-1] = index
872                        index += 1     
873                endfor
874        endfor
875       
876//      Printf "Set up SLD id matrix, making waves = %g seconds\r",(ticks-t1)/60.15
877//      t1 = ticks
878
879        for(ii=0;ii<num;ii+=1)
880                rhoi = rho[ii]
881                // find the Nij. Not sure if these are what I really need...
882                PSFIndex = psf_id_mat[ SLDLookup[rhoi] ][ SLDlookup[rhoi] ]
883                Nij[PSFIndex] += 1                      //these are just the Ni counts
884        endfor
885       
886//      Printf "get the Ni counts = %g seconds\r",(ticks-t1)/60.15
887//      t1 = ticks
888       
889// this is the double loop that has been XOPed 
890//      for(ii=0;ii<num;ii+=1)
891//              for(kk=(ii+1);kk<num;kk+=1)
892//                      val = AV_Distance(xv[ii],xv[kk],yv[ii],yv[kk],zv[ii],zv[kk]) * grid
893//                      binIndex = trunc(val/binWidth-0.5)
894//                      if(binIndex > numBins -1 )
895//                              Print "bad index"
896//                      else
897//                              //figure out which PSF the pair distance belongs to
898//                              rhoi = rho[ii]
899//                              rhok = rho[kk]
900//                              PSFIndex = psf_id_mat[ SLDLookup[rhoi] ][ SLDlookup[rhok] ]
901//                              binMatrix[binIndex][PSFIndex] += 1
902//                             
903//                      endif
904//              endfor         
905//      endfor
906
907// binMatrix is returned
908
909        nthreads=ThreadProcessorCount
910        if(nthreads == 1)
911                binSLDDistanceX(xv, yv, zv, rho, binMatrix, SLDLookup, psf_id_mat, grid, binWidth,0,numpnts(xv))
912        else
913                binSLDDistance_SLD_Threaded(xv, yv, zv, rho, binMatrix, SLDLookup, psf_id_mat, grid, binWidth)
914        endif
915       
916//      Printf "binned the distances and SLDs = %g seconds\r",(ticks-t1)/60.15
917//      t1 = ticks
918
919// Now for each Gij
920        // -- loop over q, calculating Sij
921
922        nq = numpnts(qW)
923       
924        for(jj=0;jj<numDiffPSF;jj+=1)
925       
926                Make/O/D/N=(numBins) dist_at_bin, Gij
927
928                for(ii=0;ii<numBins;ii+=1)
929                        dist_at_bin[ii] =  ii*binWidth
930                endfor
931                               
932                Gij = binMatrix[p][jj]
933               
934                // remove all of the bins with zero distance, then reset the number of Bins
935                RemoveZerosXY(Gij, dist_at_bin)
936               
937                Duplicate/O Iqr, Sij
938                Sij = 0
939                newNumBins = numpnts(Gij)
940                               
941                for(kk=0;kk<nq;kk+=1)
942                        qval = qW[kk]
943                        for(ii=0;ii<newNumBins;ii+=1)           
944                                Sij[kk] += Gij[ii]*sinc(qval*dist_at_bin[ii])
945                        endfor
946                        Sij[kk] *= 2
947//                      Sij[kk] += Nij[jj]                      // these may not be the proper Nij
948                endfor
949               
950                // keep copies before looping to the next Gij
951                Duplicate/O Sij,$("Sij"+num2str(jj))
952                Duplicate/O Gij, $("Gij"+num2str(jj))
953                Duplicate/O dist_at_bin,$("dist_at_bin"+num2str(jj))
954       
955        endfor
956
957//      Printf "Calculate Sij's and Gij's = %g seconds\r",(ticks-t1)/60.15
958        t1 = ticks
959       
960        // now add up everything to get the intensity
961        // use the psf_id_mat[][] to get the Gij index
962       
963        vol=4*Pi/3*rval*rval*rval
964        Iqr = 0
965       
966        for(kk=0;kk<nq;kk+=1)           // loop over the q-values
967       
968                qval = qW[kk]
969                fQR = PhiQR(qval,rval)
970
971                for(ii=0;ii<numDiffSLD;ii+=1)
972                        index = psf_id_mat[ii][ii]
973                       
974                        F2Q = fQR*fQR*SLD_id[ii]*SLD_id[ii]*vol*vol
975                       
976                        Iqr[kk] += Nij[index]*F2Q                       //these are just the "ii" numbers
977               
978                        for(jj=ii;jj<numDiffSLD;jj+=1) 
979                                index = psf_id_mat[ii][jj]
980                               
981                                WAVE Sij = $("root:Sij"+num2str(index))         //switch the Sij
982                               
983                                Iqr[kk] += fQR*fQR*SLD_id[ii]*SLD_id[jj]*vol*vol*Sij[kk]                                //factor of 2 removed, not sure why
984                               
985                        endfor
986                endfor
987               
988        endfor
989               
990//      Printf "Last loop to calculate I(q) = %g seconds\r",(ticks-t1)/60.15
991
992        return (0)
993End
994
995
996// support up to 8 threads at this time
997Function binSLDDistance_SLD_Threaded(xv, yv, zv, rho, binMatrix, SLDLookup, psf_id_mat, grid, binWidth)
998        Wave xv, yv, zv, rho, binMatrix, SLDLookup, psf_id_mat
999        Variable grid, binWidth
1000       
1001        Variable nthreads,mt,left,right
1002        Variable ii,num
1003
1004        nthreads=ThreadProcessorCount
1005
1006        if(nthreads > 8)
1007                nthreads = 8
1008        endif
1009       
1010        mt = ThreadGroupCreate(nthreads)
1011       
1012        ii=0
1013        left = 0
1014        num = numpnts(xv)
1015       
1016        for(ii=0;ii<nthreads;ii+=1)
1017                Duplicate/O binMatrix $("binMatrix"+num2str(ii))
1018               
1019//  this is an even spreading of the points - not appropriate here for the triangle
1020//                      Print (ii*num/nthreads),((ii+1)*num/nthreads)
1021//                      ThreadStart mt,i,Cyl_PolyRadius_T(cw,yw,xw,(ii*num/nthreads),((ii+1)*num/nthreads))
1022
1023// this splits up the triangle into equal area chunks trapezoid left(ii) -> right(ii)
1024// be sure that for the last iteration of ii, right = 1 *(num)
1025                right = 1 - sqrt(1- (ii+1)/nthreads)
1026
1027//                      Print left,right
1028               
1029                if(ii==0)
1030                        Wave binMatrix0
1031                        binMatrix0 = 0
1032                        //Print (ii*num/nthreads),((ii+1)*num/nthreads)
1033                        //ThreadStart mt,ii,binDistance_WF(xv, yv, zv, binMatrix0, grid, binWidth,(ii*num/nthreads),((ii+1)*num/nthreads))
1034                        ThreadStart mt,ii,binDistance_SLD_WF(xv, yv, zv, rho, binMatrix0, SLDLookup, psf_id_mat, grid, binWidth, left*num, right*num)
1035                endif
1036                if(ii==1)
1037                        Wave binMatrix1
1038                        binMatrix1 = 0
1039                        ThreadStart mt,ii,binDistance_SLD_WF(xv, yv, zv, rho, binMatrix1, SLDLookup, psf_id_mat, grid, binWidth, left*num, right*num)
1040                endif
1041                if(ii==2)
1042                        Wave binMatrix2
1043                        binMatrix2 = 0
1044                        ThreadStart mt,ii,binDistance_SLD_WF(xv, yv, zv, rho, binMatrix2, SLDLookup, psf_id_mat, grid, binWidth, left*num, right*num)
1045                endif
1046                if(ii==3)
1047                        Wave binMatrix3
1048                        binMatrix3 = 0
1049                        ThreadStart mt,ii,binDistance_SLD_WF(xv, yv, zv, rho, binMatrix3, SLDLookup, psf_id_mat, grid, binWidth, left*num, right*num)
1050                endif
1051                if(ii==4)
1052                        Wave binMatrix4
1053                        binMatrix4 = 0
1054                        ThreadStart mt,ii,binDistance_SLD_WF(xv, yv, zv, rho, binMatrix4, SLDLookup, psf_id_mat, grid, binWidth, left*num, right*num)
1055                endif
1056                if(ii==5)
1057                        Wave binMatrix5
1058                        binMatrix5 = 0
1059                        ThreadStart mt,ii,binDistance_SLD_WF(xv, yv, zv, rho, binMatrix5, SLDLookup, psf_id_mat, grid, binWidth, left*num, right*num)
1060                endif
1061                if(ii==6)
1062                        Wave binMatrix6
1063                        binMatrix6 = 0
1064                        ThreadStart mt,ii,binDistance_SLD_WF(xv, yv, zv, rho, binMatrix6, SLDLookup, psf_id_mat, grid, binWidth, left*num, right*num)
1065                endif
1066                if(ii==7)
1067                        Wave binMatrix7
1068                        binMatrix7 = 0
1069                        ThreadStart mt,ii,binDistance_SLD_WF(xv, yv, zv, rho, binMatrix7, SLDLookup, psf_id_mat, grid, binWidth, left*num, right*num)
1070                endif
1071               
1072                left = right
1073
1074        endfor
1075       
1076        // wait until done
1077        do
1078                variable tgs= ThreadGroupWait(mt,100)
1079        while( tgs != 0 )
1080        variable dummy= ThreadGroupRelease(mt)
1081        mt=0
1082//      Print "done with all threads"
1083
1084        //      then add them all back together
1085        if(nthreads == 1)
1086                binMatrix = binMatrix0          // add up each instance
1087        endif
1088        if(nthreads == 2)
1089                binMatrix = binMatrix0  + binMatrix1
1090        endif
1091        if(nthreads == 3)
1092                binMatrix = binMatrix0  + binMatrix1 + binMatrix2
1093        endif
1094        if(nthreads == 4)
1095                binMatrix = binMatrix0  + binMatrix1 + binMatrix2 + binMatrix3
1096        endif
1097        if(nthreads == 5)
1098                binMatrix = binMatrix0  + binMatrix1 + binMatrix2 + binMatrix3 + binMatrix4
1099        endif
1100        if(nthreads == 6)
1101                binMatrix = binMatrix0  + binMatrix1 + binMatrix2 + binMatrix3 + binMatrix4 + binMatrix5
1102        endif
1103        if(nthreads == 7)
1104                binMatrix = binMatrix0  + binMatrix1 + binMatrix2 + binMatrix3 + binMatrix4 + binMatrix5 + binMatrix6
1105        endif
1106        if(nthreads == 8)
1107                binMatrix = binMatrix0  + binMatrix1 + binMatrix2 + binMatrix3 + binMatrix4 + binMatrix5 + binMatrix6 + binMatrix7
1108        endif
1109       
1110        // then clean up
1111        KillWaves/Z binMatrix0,binMatrix1,binMatrix2,binMatrix3
1112        KillWaves/Z binMatrix4,binMatrix5,binMatrix6,binMatrix7
1113       
1114        return(0)
1115End
1116
1117// this is just a worker function to get the ThreadStart operation to compile
1118ThreadSafe Function binDistance_SLD_WF(xv, yv, zv, rho, binMatrix, SLDLookup, psf_id_mat, grid, binWidth, p1,p2)
1119        Wave xv, yv, zv, rho, binMatrix, SLDLookup, psf_id_mat
1120        Variable grid, binWidth, p1,p2
1121
1122        Variable ret
1123        ret = binSLDDistanceX(xv, yv, zv, rho, binMatrix, SLDLookup, psf_id_mat, grid, binWidth, p1,p2)
1124
1125        return(0)
1126End
Note: See TracBrowser for help on using the repository browser.