Changeset 789


Ignore:
Timestamp:
Feb 10, 2011 1:33:02 PM (11 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.

Location:
sans/XOP_Dev/MonteCarlo
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • sans/XOP_Dev/MonteCarlo/DebyeSpheres.c

    r758 r789  
    77#include "DebyeSpheres.h" 
    88 
    9 #pragma XOP_SET_STRUCT_PACKING                  // All structures are 2-byte-aligned. 
     9//#pragma XOP_SET_STRUCT_PACKING                        // All structures are 2-byte-aligned. 
    1010 
    1111// Prototypes 
     
    158158 
    159159 
    160 #pragma XOP_RESET_STRUCT_PACKING                        // All structures are 2-byte-aligned. 
    161 // All structures are 2-byte-aligned. 
     160/*       
     161  
     162 given the distances XYZ as a triplet (on a unit grid) 
     163 return the maximum distance. The calling program must multiply by 
     164 the grid dimension to get real distance 
     165  
     166 */ 
     167int 
     168maxDistanceX(DistParamPtr p) 
     169{ 
     170        double dmax,dij;                                //output dmax value, dij 
     171        double *xv,*yv,*zv;             //pointers to input xyz coordinates 
     172        int i,j; 
     173    int npt; 
     174         
     175         
     176         
     177        // check for all of the required waves 
     178        if (p->zwavH == NIL) { 
     179                SetNaN64(&p->result); 
     180                return NON_EXISTENT_WAVE; 
     181        } 
     182        if (p->ywavH == NIL) { 
     183                SetNaN64(&p->result); 
     184                return NON_EXISTENT_WAVE; 
     185        } 
     186        if (p->xwavH == NIL) { 
     187                SetNaN64(&p->result); 
     188                return NON_EXISTENT_WAVE; 
     189        } 
     190 
     191    //check to see that all are double 
     192    if(WaveType(p->zwavH) != NT_FP64 ) { 
     193        SetNaN64(&p->result); 
     194        return kExpectedNT_FP64; 
     195    } 
     196    if(WaveType(p->ywavH) != NT_FP64 ) { 
     197        SetNaN64(&p->result); 
     198        return kExpectedNT_FP64; 
     199    } 
     200    if(WaveType(p->xwavH) != NT_FP64 ) { 
     201        SetNaN64(&p->result); 
     202        return kExpectedNT_FP64; 
     203    } 
     204         
     205        // 
     206    npt = (int) WavePoints(p->xwavH);   //wavePoints returns long, number of XYZ points 
     207        xv = WaveData(p->xwavH);                //xyz locations 
     208        yv = WaveData(p->ywavH); 
     209        zv = WaveData(p->zwavH); 
     210         
     211        dmax = 0; 
     212        //do the i!=j double loop, keeping the maximum distance 
     213 
     214        for(i=0;i<npt;i+=1) { 
     215                for(j=(i+1);j<npt;j+=1) { 
     216//                      dij=XYZDistance(xv[i],xv[j],yv[i],yv[j],zv[i],zv[j]); 
     217                        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]); 
     218                        if(dij > dmax) { 
     219                                dmax = dij; 
     220                        } 
     221                } 
     222        } 
     223         
     224    p->result = dmax; 
     225         
     226    return 0; 
     227         
     228} 
     229 
     230/*       
     231  
     232 given the distances XYZ as a triplet (on a unit grid) 
     233 return the binned histogram of distances 
     234  
     235 */ 
     236int 
     237binDistanceX(BinParamPtr p) 
     238{ 
     239        double *xv,*yv,*zv,*bv;         //pointers to input xyz coordinates 
     240        int i,j; 
     241    int npt,numBins,binIndex; 
     242        double grid,binWidth,val; 
     243         
     244         
     245         
     246        // check for all of the required waves 
     247        if (p->bwavH == NIL) { 
     248                SetNaN64(&p->result); 
     249                return NON_EXISTENT_WAVE; 
     250        } 
     251        if (p->zwavH == NIL) { 
     252                SetNaN64(&p->result); 
     253                return NON_EXISTENT_WAVE; 
     254        } 
     255        if (p->ywavH == NIL) { 
     256                SetNaN64(&p->result); 
     257                return NON_EXISTENT_WAVE; 
     258        } 
     259        if (p->xwavH == NIL) { 
     260                SetNaN64(&p->result); 
     261                return NON_EXISTENT_WAVE; 
     262        } 
     263         
     264    //check to see that all are double 
     265        if(WaveType(p->bwavH) != NT_FP64 ) { 
     266        SetNaN64(&p->result); 
     267        return kExpectedNT_FP64; 
     268    } 
     269    if(WaveType(p->zwavH) != NT_FP64 ) { 
     270        SetNaN64(&p->result); 
     271        return kExpectedNT_FP64; 
     272    } 
     273    if(WaveType(p->ywavH) != NT_FP64 ) { 
     274        SetNaN64(&p->result); 
     275        return kExpectedNT_FP64; 
     276    } 
     277    if(WaveType(p->xwavH) != NT_FP64 ) { 
     278        SetNaN64(&p->result); 
     279        return kExpectedNT_FP64; 
     280    } 
     281         
     282        // 
     283    npt = (int) WavePoints(p->xwavH);   //wavePoints returns long, number of XYZ points 
     284    numBins = (int) WavePoints(p->bwavH);       //wavePoints returns long, number of points in bin wave 
     285         
     286        xv = WaveData(p->xwavH);                //xyz locations 
     287        yv = WaveData(p->ywavH); 
     288        zv = WaveData(p->zwavH); 
     289        bv = WaveData(p->bwavH); 
     290         
     291        grid = p->grid; 
     292        binWidth = p->binWidth; 
     293         
     294        //do the i!=j double loop,       
     295        for(i=0;i<npt;i+=1) { 
     296                for(j=(i+1);j<npt;j+=1) { 
     297                        val = XYZDistance(xv[i],xv[j],yv[i],yv[j],zv[i],zv[j])*grid; 
     298                        binIndex = (int)(val/binWidth-0.5); 
     299                        if(binIndex > numBins -1 ) { 
     300                                //Print "bad index" 
     301                        } else { 
     302                                bv[binIndex] += 1; 
     303                        } 
     304                         
     305                } 
     306        } 
     307         
     308    p->result = 0; 
     309         
     310    return 0; 
     311         
     312} 
     313 
     314 
     315/*       
     316  
     317 given the distances XYZ as a triplet (on a unit grid) and the SLD at each point, 
     318 return the binned histogram of distances for each of the parwise interactions 
     319 
     320 The returned binning is a matrix, and has to be assigned as such 
     321  
     322 */ 
     323int 
     324binSLDDistanceX(BinSLDParamPtr p) 
     325{ 
     326        double *xv,*yv,*zv;             //pointers to input xyz coordinates 
     327        double *rho,*SLDLook,*PSFid;    // rho and the SLD lookup vector 
     328        int i,j; 
     329    int npt,numBins,binIndex; 
     330        double grid,binWidth,val,retVal; 
     331         
     332// for accessing the 2D wave data to write the results   
     333        waveHndl wavH,PSFwavH; 
     334//      long numDimensions; 
     335//      long dimensionSizes[MAX_DIMENSIONS+1]; 
     336        double value[2];                                // Pointers used for double data. 
     337        long indices[MAX_DIMENSIONS]; 
     338//       
     339        long rhoi,rhoj,rii,rji,PSFIndex; 
     340         
     341         
     342        // check for all of the required waves 
     343        if (p->PSFidH == NIL) { 
     344                SetNaN64(&p->result); 
     345                return NON_EXISTENT_WAVE; 
     346        } 
     347        if (p->SLDLookH == NIL) { 
     348                SetNaN64(&p->result); 
     349                return NON_EXISTENT_WAVE; 
     350        } 
     351        if (p->bwavH == NIL) { 
     352                SetNaN64(&p->result); 
     353                return NON_EXISTENT_WAVE; 
     354        } 
     355        if (p->rhowavH == NIL) { 
     356                SetNaN64(&p->result); 
     357                return NON_EXISTENT_WAVE; 
     358        } 
     359        if (p->zwavH == NIL) { 
     360                SetNaN64(&p->result); 
     361                return NON_EXISTENT_WAVE; 
     362        } 
     363        if (p->ywavH == NIL) { 
     364                SetNaN64(&p->result); 
     365                return NON_EXISTENT_WAVE; 
     366        } 
     367        if (p->xwavH == NIL) { 
     368                SetNaN64(&p->result); 
     369                return NON_EXISTENT_WAVE; 
     370        } 
     371         
     372    //check to see that all are double 
     373        if(WaveType(p->PSFidH) != NT_FP64 ) { 
     374        SetNaN64(&p->result); 
     375        return kExpectedNT_FP64; 
     376    } 
     377        if(WaveType(p->SLDLookH) != NT_FP64 ) { 
     378        SetNaN64(&p->result); 
     379        return kExpectedNT_FP64; 
     380    } 
     381        if(WaveType(p->bwavH) != NT_FP64 ) { 
     382        SetNaN64(&p->result); 
     383        return kExpectedNT_FP64; 
     384    } 
     385        if(WaveType(p->rhowavH) != NT_FP64 ) { 
     386        SetNaN64(&p->result); 
     387        return kExpectedNT_FP64; 
     388    } 
     389    if(WaveType(p->zwavH) != NT_FP64 ) { 
     390        SetNaN64(&p->result); 
     391        return kExpectedNT_FP64; 
     392    } 
     393    if(WaveType(p->ywavH) != NT_FP64 ) { 
     394        SetNaN64(&p->result); 
     395        return kExpectedNT_FP64; 
     396    } 
     397    if(WaveType(p->xwavH) != NT_FP64 ) { 
     398        SetNaN64(&p->result); 
     399        return kExpectedNT_FP64; 
     400    } 
     401         
     402         
     403        // access the 2D wave data for writing using the direct method 
     404        wavH = p->bwavH; 
     405        if (wavH == NIL) 
     406                return NOWAV; 
     407        // 
     408        PSFwavH = p->PSFidH; 
     409         
     410    npt = (int) WavePoints(p->xwavH);   //wavePoints returns long, number of XYZ points 
     411    numBins = (int) WavePoints(p->bwavH);       //wavePoints returns long, number of points in bin wave 
     412         
     413        xv = WaveData(p->xwavH);                //xyz locations 
     414        yv = WaveData(p->ywavH); 
     415        zv = WaveData(p->zwavH); 
     416        rho = WaveData(p->rhowavH); 
     417        SLDLook = WaveData(p->SLDLookH); 
     418        PSFid = WaveData(p->PSFidH);                    //this one is 2D 
     419         
     420        grid = p->grid; 
     421        binWidth = p->binWidth; 
     422         
     423        //do the i!=j double loop,       
     424        for(i=0;i<npt;i+=1) { 
     425                for(j=(i+1);j<npt;j+=1) { 
     426                        val = XYZDistance(xv[i],xv[j],yv[i],yv[j],zv[i],zv[j])*grid; 
     427                        binIndex = (int)(val/binWidth-0.5); 
     428                        if(binIndex > numBins -1 ) { 
     429                                //Print "bad index" 
     430                        } else { 
     431                                rhoi = (long) rho[i];                           //get the rho value at i and j 
     432                                rhoj = (long) rho[j]; 
     433                                rii = (long) SLDLook[rhoi];                     //rho i index 
     434                                rji = (long) SLDLook[rhoj];                     //rho j index 
     435                                MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions. 
     436                                indices[0] = rii; 
     437                                indices[1] = rji; 
     438                                if (retVal = MDGetNumericWavePointValue(PSFwavH, indices, value)) 
     439                                        return retVal; 
     440                                //PSFIndex = (long) PSFid[rii][rji];            //doesn't work 
     441                                PSFIndex = (long) value[0]; 
     442                                 
     443                                //now do the assignment to the 2D 
     444                                // equivalent to binMatrix[binIndex][PSFIndex] 
     445                                 
     446                                MemClear(indices, sizeof(indices)); // Must be 0 for unused dimensions. 
     447                                indices[0] = binIndex; 
     448                                indices[1] = PSFIndex; 
     449                                if (retVal = MDGetNumericWavePointValue(wavH, indices, value)) 
     450                                        return retVal; 
     451                                value[0] += 1; // Real part 
     452                                if (retVal = MDSetNumericWavePointValue(wavH, indices, value)) 
     453                                        return retVal; 
     454                                 
     455                        } 
     456                         
     457                } 
     458        } 
     459         
     460    p->result = 0; 
     461         
     462    return 0; 
     463         
     464} 
     465 
     466 
     467///// this is directly from Numerical Recipes 
     468// -- I did change the float to double, since Igor treats all as double 
     469// and n is an int, not a pointer (seemed unnecessary) 
     470// 
     471#define MAXBIT 30 
     472#define MAXDIM 6 
     473static int iminarg1,iminarg2; 
     474#define IMIN(a,b) (iminarg1=(a),iminarg2=(b),(iminarg1) < (iminarg2) ? (iminarg1) : (iminarg2)) 
     475 
     476int 
     477SobolX(SobolParamPtr p) 
     478{ 
     479        int j,k,l; 
     480        unsigned long i,im,ipp; 
     481        static double fac; 
     482        static unsigned long in,ix[MAXDIM+1],*iu[MAXBIT+1]; 
     483        static unsigned long mdeg[MAXDIM+1]={0,1,2,3,3,4,4}; 
     484        static unsigned long ip[MAXDIM+1]={0,0,1,1,2,1,4}; 
     485        static unsigned long iv[MAXDIM*MAXBIT+1]={ 
     486                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}; 
     487         
     488        static int initDone=0; 
     489        char buf[256]; 
     490 
     491        int n=0; 
     492        double *x;              //output x vector 
     493         
     494        // check for all of the required waves 
     495        if (p->bwavH == NIL) { 
     496                SetNaN64(&p->result); 
     497                return NON_EXISTENT_WAVE; 
     498        } 
     499         
     500    //check to see that all are double 
     501        if(WaveType(p->bwavH) != NT_FP64 ) { 
     502        SetNaN64(&p->result); 
     503        return kExpectedNT_FP64; 
     504    } 
     505        x = WaveData(p->bwavH); 
     506        n = (int)(p->nIn);                      // not sure that the negative input will be properly cast to int 
     507         
     508//      sprintf(buf, "input, recast n = %g  %d\r",p->nIn, n); 
     509//      XOPNotice(buf); 
     510         
     511        if (n < 0) { 
     512                 
     513                if(initDone) { 
     514                        sprintf(buf, "Don't re-initialize\r"); 
     515                        XOPNotice(buf); 
     516                        return 0; 
     517                } 
     518                 
     519                for (j=1,k=0;j<=MAXBIT;j++,k+=MAXDIM) iu[j] = &iv[k]; 
     520                for (k=1;k<=MAXDIM;k++) { 
     521                        for (j=1;j<=mdeg[k];j++) iu[j][k] <<= (MAXBIT-j); 
     522                        for (j=mdeg[k]+1;j<=MAXBIT;j++) { 
     523                                ipp=ip[k]; 
     524                                i=iu[j-mdeg[k]][k]; 
     525                                i ^= (i >> mdeg[k]); 
     526                                for (l=mdeg[k]-1;l>=1;l--) { 
     527                                        if (ipp & 1) i ^= iu[j-l][k]; 
     528                                        ipp >>= 1; 
     529                                } 
     530                                iu[j][k]=i; 
     531                        } 
     532                } 
     533                fac=1.0/(1L << MAXBIT); 
     534                in=0; 
     535                 
     536                initDone=1; 
     537                 
     538                sprintf(buf, "Initialization loop done\r"); 
     539                XOPNotice(buf); 
     540                 
     541        } else { 
     542                im=in; 
     543                for (j=1;j<=MAXBIT;j++) { 
     544                        if (!(im & 1)) break; 
     545                        im >>= 1; 
     546                } 
     547                if (j > MAXBIT) { 
     548                        sprintf(buf, "MAXBIT too small in sobseq\r"); 
     549                        XOPNotice(buf); 
     550                } 
     551                im=(j-1)*MAXDIM; 
     552                for (k=1;k<=IMIN(n,MAXDIM);k++) { 
     553                        ix[k] ^= iv[im+k]; 
     554                        x[k-1]=ix[k]*fac;                       /// this is a real array to send back, count this one from zero 
     555                        //sprintf(buf, "calculate x[%d] = %g\r",k,ix[k]*fac); 
     556                        //XOPNotice(buf); 
     557                } 
     558                in++; 
     559        } 
     560         
     561//      sprintf(buf, "x[0],x[1] = %g  %g\r",x[0],x[1]); 
     562//      XOPNotice(buf); 
     563 
     564         
     565        p->result = 0; 
     566         
     567    return 0; 
     568} 
     569 
     570#undef MAXBIT 
     571#undef MAXDIM 
     572 
     573 
     574//#pragma XOP_RESET_STRUCT_PACKING                      // All structures are 2-byte-aligned. 
  • sans/XOP_Dev/MonteCarlo/DebyeSpheres.h

    r758 r789  
    99#include <stdio.h> 
    1010 
     11// All structures passed to Igor are two-byte aligned. 
     12#include "XOPStructureAlignmentTwoByte.h" 
     13 
    1114// Parameters passed from Igor to the function, in reverse order of the parameter list 
    1215// result is the last parameter, always. 
     
    1518    double   Rprimary;  //primary sphere radius 
    1619        waveHndl rhowavH;       // rho at xyz!!! 
    17         waveHndl zwavH; // z coordinate. ALL are expected to be SP waves 
     20        waveHndl zwavH; // z coordinate. ALL are expected to be DP waves 
    1821        waveHndl ywavH; // y coordinate. 
    1922        waveHndl xwavH; // x coordinate. 
     
    2225        double result; 
    2326}AltiParams, *AltiParamsPtr;    
    24   
    25   
     27 
     28 
     29// Parameters passed from Igor to the function, in reverse order of the parameter list 
     30// result is the last parameter, always. 
     31typedef struct DistParam { 
     32        waveHndl zwavH; // z coordinate. ALL are expected to be DP waves 
     33        waveHndl ywavH; // y coordinate. 
     34        waveHndl xwavH; // x coordinate. 
     35//      void* tp;                       //unused void for threadsafe functions 
     36        double result; 
     37}DistParam, *DistParamPtr;    
     38 
     39// Parameters passed from Igor to the function, in reverse order of the parameter list 
     40// result is the last parameter, always. 
     41typedef struct BinParam { 
     42        double binWidth; 
     43        double grid; 
     44        waveHndl bwavH; // bin wave 
     45        waveHndl zwavH; // z coordinate. ALL are expected to be DP waves 
     46        waveHndl ywavH; // y coordinate. 
     47        waveHndl xwavH; // x coordinate. 
     48        //      void* tp;                       //unused void for threadsafe functions 
     49        double result; 
     50}BinParam, *BinParamPtr;    
     51 
     52// Parameters passed from Igor to the function, in reverse order of the parameter list 
     53// result is the last parameter, always. 
     54typedef struct BinSLDParam { 
     55        double binWidth; 
     56        double grid; 
     57        waveHndl PSFidH;        // psf id matrix, this is 2D 
     58        waveHndl SLDLookH;      // SLD lookup wave, this is 1D 
     59        waveHndl bwavH; // bin wave, this is 2D 
     60        waveHndl rhowavH;               //rho @ XYZ 
     61        waveHndl zwavH; // z coordinate. ALL are expected to be DP waves 
     62        waveHndl ywavH; // y coordinate. 
     63        waveHndl xwavH; // x coordinate. 
     64        //      void* tp;                       //unused void for threadsafe functions 
     65        double result; 
     66}BinSLDParam, *BinSLDParamPtr;    
     67 
     68 
     69// Parameters passed from Igor to the function, in reverse order of the parameter list 
     70// result is the last parameter, always. 
     71typedef struct SobolParam { 
     72        waveHndl bwavH; // result wave 
     73        double nIn;             // initialization, or dimension - this is really an int 
     74        //      void* tp;                       //unused void for threadsafe functions 
     75        double result; 
     76}SobolParam, *SobolParamPtr;    
     77 
     78 
     79 
     80#include "XOPStructureAlignmentReset.h" 
     81 
     82 
    2683double PhiQR(double qval, double rval); 
    2784double XYZDistance(double x1, double x2,double y1, double y2,double z1, double z2); 
    2885int DebyeSpheresX(AltiParamsPtr p); 
     86int maxDistanceX(DistParamPtr p); 
     87int binDistanceX(BinParamPtr p); 
     88int binSLDDistanceX(BinSLDParamPtr p); 
     89int SobolX(SobolParamPtr p); 
     90 
  • sans/XOP_Dev/MonteCarlo/MonteCarlo.c

    r711 r789  
    114114        seed = (long)results[0]; 
    115115         
    116 //      sprintf(buf, "input seed = %ld\r", seed); 
    117 //      XOPNotice(buf); 
     116        //      sprintf(buf, "input seed = %ld\r", seed); 
     117        //      XOPNotice(buf); 
    118118         
    119119        if(seed >= 0) { 
  • sans/XOP_Dev/MonteCarlo/MonteCarlo.r

    r758 r789  
    104104                }, 
    105105                 
     106                "maxDistanceX",         /* function name */ 
     107                F_UTIL | F_EXTERNAL,                                    /* function category == UTILITY */ 
     108                NT_FP64,                                                        /* return value type */                  
     109                { 
     110                        NT_FP64 + WAVE_TYPE, 
     111                        NT_FP64 + WAVE_TYPE, 
     112                        NT_FP64 + WAVE_TYPE, 
     113                }, 
     114                 
     115                "binDistanceX",         /* function name */ 
     116                F_UTIL | F_EXTERNAL,                                    /* function category == UTILITY */ 
     117                NT_FP64,                                                        /* return value type */                  
     118                { 
     119                        NT_FP64 + WAVE_TYPE,                    // XYZ waves 
     120                        NT_FP64 + WAVE_TYPE, 
     121                        NT_FP64 + WAVE_TYPE, 
     122                        NT_FP64 + WAVE_TYPE,                    // distance Bin wave (this is returned) 
     123                        NT_FP64,                                                // grid 
     124                        NT_FP64,                                                // bin width 
     125                }, 
     126                 
     127                "SobolX",               /* function name */ 
     128                F_UTIL | F_EXTERNAL,                                    /* function category == UTILITY */ 
     129                NT_FP64,                                                        /* return value type */                  
     130                { 
     131                        NT_FP64,                                                // requested dimension (this is really an integer) 
     132                        NT_FP64 + WAVE_TYPE,                    // result wave 
     133                }, 
     134                 
     135                "binSLDDistanceX",              /* function name */ 
     136                F_UTIL | F_EXTERNAL,                                    /* function category == UTILITY */ 
     137                NT_FP64,                                                        /* return value type */                  
     138                { 
     139                        NT_FP64 + WAVE_TYPE,                    // XYZ waves 
     140                        NT_FP64 + WAVE_TYPE, 
     141                        NT_FP64 + WAVE_TYPE, 
     142                        NT_FP64 + WAVE_TYPE,                    // rho wave 
     143                        NT_FP64 + WAVE_TYPE,                    // distance Bin wave, 2D (this is returned) 
     144                        NT_FP64 + WAVE_TYPE,                    // SLD lookup 
     145                        NT_FP64 + WAVE_TYPE,                    // PSF id matrix 
     146                        NT_FP64,                                                // grid 
     147                        NT_FP64,                                                // bin width 
     148                }, 
     149                 
    106150        } 
    107151}; 
  • sans/XOP_Dev/MonteCarlo/MonteCarlo_Main.c

    r758 r789  
    371371                        return((long)Monte_SANSX4); 
    372372                        break; 
     373                case 5:                                         //  
     374                        return((long)maxDistanceX); 
     375                        break; 
     376                case 6:                                         //  
     377                        return((long)binDistanceX); 
     378                        break; 
     379                case 7:                                         //  
     380                        return((long)SobolX); 
     381                        break; 
     382                case 8:                                         //  
     383                        return((long)binSLDDistanceX); 
     384                        break; 
    373385        } 
    374386        return(NIL); 
Note: See TracChangeset for help on using the changeset viewer.