Changeset 789 for sans/XOP_Dev/MonteCarlo/DebyeSpheres.c

Timestamp:

Feb 10, 2011 1:33:02 PM (12 years ago)

Author:

srkline

Message:

Added XOP-ized functions for doing the distance-binned Debye Spheres calculation. Two versions are added, one that includes the SLDs and one that does not.

Also added a pseudo-random sequence generator, SobolX, that will generate pseudo-random 2D or 3D distributions. This fills space more uniformly than a random generation.

File:

• sans/XOP_Dev/MonteCarlo/DebyeSpheres.c (modified) (2 diffs)

sans/XOP_Dev/MonteCarlo/DebyeSpheres.c

@@ -7 +7 @@
 #include "DebyeSpheres.h"
 
-#pragma XOP_SET_STRUCT_PACKING                  // All structures are 2-byte-aligned.
+//#pragma XOP_SET_STRUCT_PACKING                        // All structures are 2-byte-aligned.
 
 // Prototypes
@@ -158 +158 @@
 
 
-#pragma XOP_RESET_STRUCT_PACKING                        // All structures are 2-byte-aligned.
-// All structures are 2-byte-aligned.
+/*      
+ 
+ given the distances XYZ as a triplet (on a unit grid)
+ return the maximum distance. The calling program must multiply by
+ the grid dimension to get real distance
+ 
+ */
+int
+maxDistanceX(DistParamPtr p)
+{
+        double dmax,dij;                                //output dmax value, dij
+        double *xv,*yv,*zv;             //pointers to input xyz coordinates
+        int i,j;
+    int npt;
+        
+        
+        
+        // check for all of the required waves
+        if (p->zwavH == NIL) {
+                SetNaN64(&p->result);
+                return NON_EXISTENT_WAVE;
+        }
+        if (p->ywavH == NIL) {
+                SetNaN64(&p->result);
+                return NON_EXISTENT_WAVE;
+        }
+        if (p->xwavH == NIL) {
+                SetNaN64(&p->result);
+                return NON_EXISTENT_WAVE;
+        }
+
+    //check to see that all are double
+    if(WaveType(p->zwavH) != NT_FP64 ) {
+        SetNaN64(&p->result);
+        return kExpectedNT_FP64;
+    }
+    if(WaveType(p->ywavH) != NT_FP64 ) {
+        SetNaN64(&p->result);
+        return kExpectedNT_FP64;
+    }
+    if(WaveType(p->xwavH) != NT_FP64 ) {
+        SetNaN64(&p->result);
+        return kExpectedNT_FP64;
+    }
+        
+        //
+    npt = (int) WavePoints(p->xwavH);   //wavePoints returns long, number of XYZ points
+        xv = WaveData(p->xwavH);                //xyz locations
+        yv = WaveData(p->ywavH);
+        zv = WaveData(p->zwavH);
+        
+        dmax = 0;
+        //do the i!=j double loop, keeping the maximum distance
+
+        for(i=0;i<npt;i+=1) {
+                for(j=(i+1);j<npt;j+=1) {
+//                      dij=XYZDistance(xv[i],xv[j],yv[i],yv[j],zv[i],zv[j]);
+                        dij = (xv[i]-xv[j])*(xv[i]-xv[j]) + (yv[i]-yv[j])*(yv[i]-yv[j]) + (zv[i]-zv[j])*(zv[i]-zv[j]);
+                        if(dij > dmax) {
+                                dmax = dij;
+                        }
+                }
+        }
+        
+    p->result = dmax;
+        
+    return 0;
+        
+}
+
+/*      
+ 
+ given the distances XYZ as a triplet (on a unit grid)
+ return the binned histogram of distances
+ 
+ */
+int
+binDistanceX(BinParamPtr p)
+{
+        double *xv,*yv,*zv,*bv;         //pointers to input xyz coordinates
+        int i,j;
+    int npt,numBins,binIndex;
+        double grid,binWidth,val;
+        
+        
+        
+        // check for all of the required waves
+        if (p->bwavH == NIL) {
+                SetNaN64(&p->result);
+                return NON_EXISTENT_WAVE;
+        }
+        if (p->zwavH == NIL) {
+                SetNaN64(&p->result);
+                return NON_EXISTENT_WAVE;
+        }
+        if (p->ywavH == NIL) {
+                SetNaN64(&p->result);
+                return NON_EXISTENT_WAVE;
+        }
+        if (p->xwavH == NIL) {
+                SetNaN64(&p->result);
+                return NON_EXISTENT_WAVE;
+        }
+        
+    //check to see that all are double
+        if(WaveType(p->bwavH) != NT_FP64 ) {
+        SetNaN64(&p->result);
+        return kExpectedNT_FP64;
+    }
+    if(WaveType(p->zwavH) != NT_FP64 ) {
+        SetNaN64(&p->result);
+        return kExpectedNT_FP64;
+    }
+    if(WaveType(p->ywavH) != NT_FP64 ) {
+        SetNaN64(&p->result);
+        return kExpectedNT_FP64;
+    }
+    if(WaveType(p->xwavH) != NT_FP64 ) {
+        SetNaN64(&p->result);
+        return kExpectedNT_FP64;
+    }
+        
+        //
+    npt = (int) WavePoints(p->xwavH);   //wavePoints returns long, number of XYZ points
+    numBins = (int) WavePoints(p->bwavH);       //wavePoints returns long, number of points in bin wave
+        
+        xv = WaveData(p->xwavH);                //xyz locations
+        yv = WaveData(p->ywavH);
+        zv = WaveData(p->zwavH);
+        bv = WaveData(p->bwavH);
+        
+        grid = p->grid;
+        binWidth = p->binWidth;
+        
+        //do the i!=j double loop,      
+        for(i=0;i<npt;i+=1) {
+                for(j=(i+1);j<npt;j+=1) {
+                        val = XYZDistance(xv[i],xv[j],yv[i],yv[j],zv[i],zv[j])*grid;
+                        binIndex = (int)(val/binWidth-0.5);
+                        if(binIndex > numBins -1 ) {
+                                //Print "bad index"
+                        } else {
+                                bv[binIndex] += 1;
+                        }
+                        
+                }
+        }
+        
+    p->result = 0;
+        
+    return 0;
+        
+}
+
+
+/*      
+ 
+ given the distances XYZ as a triplet (on a unit grid) and the SLD at each point,
+ return the binned histogram of distances for each of the parwise interactions
+
+ The returned binning is a matrix, and has to be assigned as such
+ 
+ */
+int
+binSLDDistanceX(BinSLDParamPtr p)
+{
+        double *xv,*yv,*zv;             //pointers to input xyz coordinates
+        double *rho,*SLDLook,*PSFid;    // rho and the SLD lookup vector
+        int i,j;
+    int npt,numBins,binIndex;
+        double grid,binWidth,val,retVal;
+        
+// for accessing the 2D wave data to write the results  
+        waveHndl wavH,PSFwavH;
+//      long numDimensions;
+//      long dimensionSizes[MAX_DIMENSIONS+1];
+        double value[2];                                // Pointers used for double data.
+        long indices[MAX_DIMENSIONS];
+//      
+        long rhoi,rhoj,rii,rji,PSFIndex;
+        
+        
+        // check for all of the required waves
+        if (p->PSFidH == NIL) {
+                SetNaN64(&p->result);
+                return NON_EXISTENT_WAVE;
+        }
+        if (p->SLDLookH == NIL) {
+                SetNaN64(&p->result);
+                return NON_EXISTENT_WAVE;
+        }
+        if (p->bwavH == NIL) {
+                SetNaN64(&p->result);
+                return NON_EXISTENT_WAVE;
+        }
+        if (p->rhowavH == NIL) {
+                SetNaN64(&p->result);
+                return NON_EXISTENT_WAVE;
+        }
+        if (p->zwavH == NIL) {
+                SetNaN64(&p->result);
+                return NON_EXISTENT_WAVE;
+        }
+        if (p->ywavH == NIL) {
+                SetNaN64(&p->result);
+                return NON_EXISTENT_WAVE;
+        }
+        if (p->xwavH == NIL) {
+                SetNaN64(&p->result);
+                return NON_EXISTENT_WAVE;
+        }
+        
+    //check to see that all are double
+        if(WaveType(p->PSFidH) != NT_FP64 ) {
+        SetNaN64(&p->result);
+        return kExpectedNT_FP64;
+    }
+        if(WaveType(p->SLDLookH) != NT_FP64 ) {
+        SetNaN64(&p->result);
+        return kExpectedNT_FP64;
+    }
+        if(WaveType(p->bwavH) != NT_FP64 ) {
+        SetNaN64(&p->result);
+        return kExpectedNT_FP64;
+    }
+        if(WaveType(p->rhowavH) != NT_FP64 ) {
+        SetNaN64(&p->result);
+        return kExpectedNT_FP64;
+    }
+    if(WaveType(p->zwavH) != NT_FP64 ) {
+        SetNaN64(&p->result);
+        return kExpectedNT_FP64;
+    }
+    if(WaveType(p->ywavH) != NT_FP64 ) {
+        SetNaN64(&p->result);
+        return kExpectedNT_FP64;
+    }
+    if(WaveType(p->xwavH) != NT_FP64 ) {
+        SetNaN64(&p->result);
+        return kExpectedNT_FP64;
+    }
+        
+        
+        // access the 2D wave data for writing using the direct method
+        wavH = p->bwavH;
+        if (wavH == NIL)
+                return NOWAV;
+        //
+        PSFwavH = p->PSFidH;
+        
+    npt = (int) WavePoints(p->xwavH);   //wavePoints returns long, number of XYZ points
+    numBins = (int) WavePoints(p->bwavH);       //wavePoints returns long, number of points in bin wave
+        
+        xv = WaveData(p->xwavH);                //xyz locations
+        yv = WaveData(p->ywavH);
+        zv = WaveData(p->zwavH);
+        rho = WaveData(p->rhowavH);
+        SLDLook = WaveData(p->SLDLookH);
+        PSFid = WaveData(p->PSFidH);                    //this one is 2D
+        
+        grid = p->grid;
+        binWidth = p->binWidth;
+        
+        //do the i!=j double loop,      
+        for(i=0;i<npt;i+=1) {
+                for(j=(i+1);j<npt;j+=1) {
+                        val = XYZDistance(xv[i],xv[j],yv[i],yv[j],zv[i],zv[j])*grid;
+                        binIndex = (int)(val/binWidth-0.5);
+                        if(binIndex > numBins -1 ) {
+                                //Print "bad index"
+                        } else {
+                                rhoi = (long) rho[i];                           //get the rho value at i and j
+                                rhoj = (long) rho[j];
+                                rii = (long) SLDLook[rhoi];                     //rho i index
+                                rji = (long) SLDLook[rhoj];                     //rho j index
+                                MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
+                                indices[0] = rii;
+                                indices[1] = rji;
+                                if (retVal = MDGetNumericWavePointValue(PSFwavH, indices, value))
+                                        return retVal;
+                                //PSFIndex = (long) PSFid[rii][rji];            //doesn't work
+                                PSFIndex = (long) value[0];
+                                
+                                //now do the assignment to the 2D
+                                // equivalent to binMatrix[binIndex][PSFIndex]
+                                
+                                MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions.
+                                indices[0] = binIndex;
+                                indices[1] = PSFIndex;
+                                if (retVal = MDGetNumericWavePointValue(wavH, indices, value))
+                                        return retVal;
+                                value[0] += 1; // Real part
+                                if (retVal = MDSetNumericWavePointValue(wavH, indices, value))
+                                        return retVal;
+                                
+                        }
+                        
+                }
+        }
+        
+    p->result = 0;
+        
+    return 0;
+        
+}
+
+
+///// this is directly from Numerical Recipes
+// -- I did change the float to double, since Igor treats all as double
+// and n is an int, not a pointer (seemed unnecessary)
+//
+#define MAXBIT 30
+#define MAXDIM 6
+static int iminarg1,iminarg2;
+#define IMIN(a,b) (iminarg1=(a),iminarg2=(b),(iminarg1) < (iminarg2) ? (iminarg1) : (iminarg2))
+
+int
+SobolX(SobolParamPtr p)
+{
+        int j,k,l;
+        unsigned long i,im,ipp;
+        static double fac;
+        static unsigned long in,ix[MAXDIM+1],*iu[MAXBIT+1];
+        static unsigned long mdeg[MAXDIM+1]={0,1,2,3,3,4,4};
+        static unsigned long ip[MAXDIM+1]={0,0,1,1,2,1,4};
+        static unsigned long iv[MAXDIM*MAXBIT+1]={
+                0,1,1,1,1,1,1,3,1,3,3,1,1,5,7,7,3,3,5,15,11,5,15,13,9};
+        
+        static int initDone=0;
+        char buf[256];
+
+        int n=0;
+        double *x;              //output x vector
+        
+        // check for all of the required waves
+        if (p->bwavH == NIL) {
+                SetNaN64(&p->result);
+                return NON_EXISTENT_WAVE;
+        }
+        
+    //check to see that all are double
+        if(WaveType(p->bwavH) != NT_FP64 ) {
+        SetNaN64(&p->result);
+        return kExpectedNT_FP64;
+    }
+        x = WaveData(p->bwavH);
+        n = (int)(p->nIn);                      // not sure that the negative input will be properly cast to int
+        
+//      sprintf(buf, "input, recast n = %g  %d\r",p->nIn, n);
+//      XOPNotice(buf);
+        
+        if (n < 0) {
+                
+                if(initDone) {
+                        sprintf(buf, "Don't re-initialize\r");
+                        XOPNotice(buf);
+                        return 0;
+                }
+                
+                for (j=1,k=0;j<=MAXBIT;j++,k+=MAXDIM) iu[j] = &iv[k];
+                for (k=1;k<=MAXDIM;k++) {
+                        for (j=1;j<=mdeg[k];j++) iu[j][k] <<= (MAXBIT-j);
+                        for (j=mdeg[k]+1;j<=MAXBIT;j++) {
+                                ipp=ip[k];
+                                i=iu[j-mdeg[k]][k];
+                                i ^= (i >> mdeg[k]);
+                                for (l=mdeg[k]-1;l>=1;l--) {
+                                        if (ipp & 1) i ^= iu[j-l][k];
+                                        ipp >>= 1;
+                                }
+                                iu[j][k]=i;
+                        }
+                }
+                fac=1.0/(1L << MAXBIT);
+                in=0;
+                
+                initDone=1;
+                
+                sprintf(buf, "Initialization loop done\r");
+                XOPNotice(buf);
+                
+        } else {
+                im=in;
+                for (j=1;j<=MAXBIT;j++) {
+                        if (!(im & 1)) break;
+                        im >>= 1;
+                }
+                if (j > MAXBIT) {
+                        sprintf(buf, "MAXBIT too small in sobseq\r");
+                        XOPNotice(buf);
+                }
+                im=(j-1)*MAXDIM;
+                for (k=1;k<=IMIN(n,MAXDIM);k++) {
+                        ix[k] ^= iv[im+k];
+                        x[k-1]=ix[k]*fac;                       /// this is a real array to send back, count this one from zero
+                        //sprintf(buf, "calculate x[%d] = %g\r",k,ix[k]*fac);
+                        //XOPNotice(buf);
+                }
+                in++;
+        }
+        
+//      sprintf(buf, "x[0],x[1] = %g  %g\r",x[0],x[1]);
+//      XOPNotice(buf);
+
+        
+        p->result = 0;
+        
+    return 0;
+}
+
+#undef MAXBIT
+#undef MAXDIM
+
+
+//#pragma XOP_RESET_STRUCT_PACKING                      // All structures are 2-byte-aligned.