Ignore:
Timestamp:
Jun 16, 2010 2:10:47 PM (12 years ago)
Author:
srkline
Message:

SA_Includes_v410 : now include Smear_2D

PeakGauss_2D, Sphere_2D : included threaded resolution smearing calculations for testing

DataSetHandling? : Included a quick and dirty batch converter for XML->6-col. See the top
of the file for the command to run

GaussUtils? : re-define the ResSemear_2D_AAOStruct. Relocated q-value and phi calculations from
RawWindowHook? to this file so they would be available for reduction and analysis

Smear_2D : now has a generic (non-threaded) smearing routine. Threading must be done in
individual functions since FUNCREF can't be passed to threads (plus a few other issues)

PlotUtils_2D : updated loader for new QxQy? columns. Fixes to Wrapper_2D to enable smeared fits

RawWindowHook? : removed Q-value calculation functions and moved these to GaussUtils?

WriteQIS : now writes out 8-columns for QxQy? data, defining the resolution
function in terms of directions parallel and perpendicular to Q. TEMPORARILY in the data
file an error in intensity is generated that is SQRT(I), being careful to
replace any NaN, inf, or zero with an average error value

MultiScatter_MonteCarlo_2D : 4-processor aware

NCNR_Utils : 2D resolution calculation is now in terms of parallel and perpendicular
rather than x and y. Gravity is included in the y-component

Location:
sans/Dev/trunk/NCNR_User_Procedures/Analysis
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/Models_2D/PeakGauss2D_v40.ipf

    r570 r708  
    157157 
    158158//non-threaded version of the function 
    159 Function PeakGauss2D_noThread(cw,zw,xw,yw) 
     159ThreadSafe Function PeakGauss2D_noThread(cw,zw,xw,yw) 
    160160        WAVE cw,zw, xw,yw 
    161161         
    162162#if exists("PeakGauss2DX")                      //to hide the function if XOP not installed 
    163          
    164163        zw= PeakGauss2DX(cw,xw,yw) 
    165164#else 
     
    175174        Struct ResSmear_2D_AAOStruct &s 
    176175         
    177 //      Smear_2DModel_20(PeakGauss2D_noThread,s) 
    178         Smear_2DModel_5(PeakGauss2D_noThread,s) 
     176        //no threads 
     177//      Smear_2DModel_PP(PeakGauss2D_noThread,s,10) 
     178 
     179        //or threaded... 
     180        SmearPeakGauss2D_THR(s,10) 
     181         
    179182        return(0) 
    180183end 
     
    190193        WAVE qy = $(DF+str+"_qy") 
    191194        WAVE qz = $(DF+str+"_qz") 
    192         WAVE sigQx = $(DF+str+"_sigQx") 
    193         WAVE sigQy = $(DF+str+"_sigQy") 
     195        WAVE sQpl = $(DF+str+"_sQpl") 
     196        WAVE sQpp = $(DF+str+"_sQpp") 
    194197        WAVE shad = $(DF+str+"_fs") 
    195198         
     
    197200        WAVE s.coefW = coefW     
    198201        WAVE s.zw = resultW      
    199         WAVE s.qx = qx 
    200         WAVE s.qy = qy 
     202        WAVE s.xw[0] = qx 
     203        WAVE s.xw[1] = qy 
    201204        WAVE s.qz = qz 
    202         WAVE s.sigQx = sigQx 
    203         WAVE s.sigQy = sigQy 
     205        WAVE s.sQpl = sQpl 
     206        WAVE s.sQpp = sQpp 
    204207        WAVE s.fs = shad 
    205208         
     
    220223//      retval = fPeak_Gauss_model(w,qval) 
    221224         
    222         RETURN(retVal) 
    223 END 
     225        return(retVal) 
     226End 
     227 
     228/////////////// 
     229// 
     230// this is the threaded version, that dispatches the calculation out to threads 
     231// 
     232// must be written specific to each 2D function 
     233//  
     234Function SmearPeakGauss2D_THR(s,nord) 
     235        Struct ResSmear_2D_AAOStruct &s 
     236        Variable nord 
     237         
     238        String weightStr,zStr 
     239         
     240// create all of the necessary quadrature waves here - rather than inside a threadsafe function 
     241        switch(nord)     
     242                case 5:          
     243                        weightStr="gauss5wt" 
     244                        zStr="gauss5z" 
     245                        if (WaveExists($weightStr) == 0) 
     246                                Make/O/D/N=(nord) $weightStr,$zStr 
     247                                Make5GaussPoints($weightStr,$zStr)       
     248                        endif 
     249                        break                            
     250                case 10:                 
     251                        weightStr="gauss10wt" 
     252                        zStr="gauss10z" 
     253                        if (WaveExists($weightStr) == 0) 
     254                                Make/O/D/N=(nord) $weightStr,$zStr 
     255                                Make10GaussPoints($weightStr,$zStr)      
     256                        endif 
     257                        break                            
     258                case 20:                 
     259                        weightStr="gauss20wt" 
     260                        zStr="gauss20z" 
     261                        if (WaveExists($weightStr) == 0) 
     262                                Make/O/D/N=(nord) $weightStr,$zStr 
     263                                Make20GaussPoints($weightStr,$zStr)      
     264                        endif 
     265                        break 
     266                default:                                                         
     267                        Abort "Smear_2DModel_PP_Threaded called with invalid nord value"                                         
     268        endswitch 
     269         
     270        Wave/Z wt = $weightStr 
     271        Wave/Z xi = $zStr               // wave references to pass 
     272 
     273        Variable npt=numpnts(s.xw[0]) 
     274        Variable i,nthreads= ThreadProcessorCount 
     275        variable mt= ThreadGroupCreate(nthreads) 
     276 
     277        Variable t1=StopMSTimer(-2) 
     278         
     279        for(i=0;i<nthreads;i+=1) 
     280        //      Print (i*npt/nthreads),((i+1)*npt/nthreads-1) 
     281                ThreadStart mt,i,SmearPeakGauss2D_T(s.coefW,s.xw[0],s.xw[1],s.qz,s.sQpl,s.sQpp,s.fs,s.zw,wt,xi,(i*npt/nthreads),((i+1)*npt/nthreads-1),nord) 
     282        endfor 
     283 
     284        do 
     285                variable tgs= ThreadGroupWait(mt,100) 
     286        while( tgs != 0 ) 
     287 
     288        variable dummy= ThreadGroupRelease(mt) 
     289         
     290// comment out the threading + uncomment this for testing to make sure that the single thread works 
     291//      nThreads=1 
     292//      SmearSphere2D_T(s.coefW,s.xw[0],s.xw[1],s.qz,s.sQpl,s.sQpp,s.fs,s.zw,wt,xi,(i*npt/nthreads),((i+1)*npt/nthreads-1),nord) 
     293 
     294        Print "elapsed time = ",(StopMSTimer(-2) - t1)/1e6 
     295         
     296        return(0) 
     297end 
     298 
     299// 
     300// - worker function for threads of PeakGauss2D 
     301// 
     302ThreadSafe Function SmearPeakGauss2D_T(coef,qxw,qyw,qzw,sxw,syw,fsw,zw,wt,xi,pt1,pt2,nord) 
     303        WAVE coef,qxw,qyw,qzw,sxw,syw,fsw,zw,wt,xi 
     304        Variable pt1,pt2,nord 
     305         
     306// now passed in.... 
     307//      Wave wt = $weightStr 
     308//      Wave xi = $zStr          
     309 
     310        Variable ii,jj,kk,num 
     311        Variable qx,qy,qz,qval,sx,sy,fs 
     312        Variable qy_pt,qx_pt,res_x,res_y,answer,sumIn,sumOut 
     313         
     314        Variable normFactor,phi,theta,maxSig,numStdDev=3 
     315         
     316/// keep these waves local 
     317        Make/O/D/N=(nord) fcnRet,xptW,res_tot,yptW 
     318         
     319        // now just loop over the points as specified 
     320         
     321        answer=0 
     322         
     323        Variable spl,spp,apl,app,bpl,bpp,phi_pt,qpl_pt 
     324 
     325        //loop over q-values 
     326        for(ii=pt1;ii<(pt2+1);ii+=1) 
     327                 
     328                qx = qxw[ii] 
     329                qy = qyw[ii] 
     330                qz = qzw[ii] 
     331                qval = sqrt(qx^2+qy^2+qz^2) 
     332                spl = sxw[ii] 
     333                spp = syw[ii] 
     334                fs = fsw[ii] 
     335                 
     336                normFactor = 2*pi*spl*spp 
     337                 
     338                phi = FindPhi(qx,qy) 
     339                 
     340                apl = -numStdDev*spl + qval             //parallel = q integration limits 
     341                bpl = numStdDev*spl + qval 
     342                app = -numStdDev*spp + phi              //perpendicular = phi integration limits 
     343                bpp = numStdDev*spp + phi 
     344                 
     345                //make sure the limits are reasonable. 
     346                if(apl < 0) 
     347                        apl = 0 
     348                endif 
     349                // do I need to specially handle limits when phi ~ 0? 
     350         
     351                 
     352                sumOut = 0 
     353                for(jj=0;jj<nord;jj+=1)         // call phi the "outer' 
     354                        phi_pt = (xi[jj]*(bpp-app)+app+bpp)/2 
     355 
     356                        sumIn=0 
     357                        for(kk=0;kk<nord;kk+=1)         //at phi, integrate over Qpl 
     358 
     359                                qpl_pt = (xi[kk]*(bpl-apl)+apl+bpl)/2 
     360                                 
     361                                FindQxQy(qpl_pt,phi_pt,qx_pt,qy_pt)             //find the corresponding QxQy to the Q,phi 
     362                                yPtw[kk] = qy_pt                                        //phi is the same in this loop, but qy is not 
     363                                xPtW[kk] = qx_pt                                        //qx is different here too, as we're varying Qpl 
     364                                 
     365                                res_tot[kk] = exp(-0.5*( (qpl_pt-qval)^2/spl/spl + (phi_pt-phi)^2/spp/spp ) ) 
     366                                res_tot[kk] /= normFactor 
     367//                              res_tot[kk] *= fs 
     368 
     369                        endfor 
     370                         
     371                        PeakGauss2D_noThread(coef,fcnRet,xptw,yptw)                     //fcn passed in is an AAO 
     372                         
     373                        //sumIn += wt[jj]*wt[kk]*res_tot*fcnRet[0] 
     374                        fcnRet *= wt[jj]*wt*res_tot 
     375                        // 
     376                        answer += (bpl-apl)/2.0*sum(fcnRet)             // 
     377                endfor 
     378 
     379                answer *= (bpp-app)/2.0 
     380                zw[ii] = answer 
     381        endfor 
     382         
     383        return(0) 
     384end 
     385 
  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/Models_2D/Sphere_2D_v40.ipf

    r570 r708  
    108108End 
    109109 
    110  
    111 //threaded version of the function 
    112 ThreadSafe Function Sphere2D_T(cw,zw,xw,yw,p1,p2) 
    113         WAVE cw,zw,xw,yw 
    114         Variable p1,p2 
    115          
    116 #if exists("Sphere_2DX")                        //to hide the function if XOP not installed 
    117          
    118         zw[p1,p2]= Sphere_2DX(cw,xw,yw) 
    119  
    120 #endif 
    121  
    122         return 0 
    123 End 
    124  
    125110// 
    126111//  Fit function that is actually a wrapper to dispatch the calculation to N threads 
     
    140125 
    141126#endif 
     127 
     128        return(0) 
     129End 
     130 
     131 
     132/// 
     133//// keep this section as an example 
     134// 
    142135 
    143136//      Variable npt=numpnts(yw) 
     
    160153////    Print "elapsed time = ",(StopMSTimer(-2) - t1)/1e6 
    161154//       
    162         return(0) 
    163 End 
    164  
    165  
    166 //non-threaded version of the function 
    167 Function Sphere2D_noThread(cw,zw,xw,yw) 
     155// 
     156////// end example of threading 
     157 
     158 
     159 
     160//threaded version of the function 
     161//ThreadSafe Function Sphere2D_T(cw,zw,xw,yw,p1,p2) 
     162//      WAVE cw,zw,xw,yw 
     163//      Variable p1,p2 
     164//       
     165//#if exists("Sphere_2DX")                      //to hide the function if XOP not installed 
     166//       
     167//      zw[p1,p2]= Sphere_2DX(cw,xw,yw) 
     168// 
     169//#endif 
     170// 
     171//      return 0 
     172//End 
     173 
     174 
     175//non-threaded version of the function, necessary for the smearing calculation 
     176// -- the smearing calculation can only calculate (nord) points at a time. 
     177// 
     178ThreadSafe Function Sphere2D_noThread(cw,zw,xw,yw) 
    168179        WAVE cw,zw, xw,yw 
    169180         
    170181#if exists("Sphere_2DX")                        //to hide the function if XOP not installed 
    171          
    172182        zw= Sphere_2DX(cw,xw,yw) 
    173  
    174183#endif 
    175184 
     
    177186End 
    178187 
    179 // I think I did this because when I do the quadrature loops I'm calling the AAO with 1-pt waves, so threading 
    180 // would just be a slowdown 
     188 
     189//  when I do the quadrature loops I'm calling the AAO with 1-pt waves, so threading 
     190// is just a slowdown (by more than a factor of 2 !!) 
     191// 
     192// - the threading gives a clean speedup of 2 for N=2, even for this simple calculation 
     193// 
     194// nord = 5,10,20 allowed 
     195// 
    181196Function SmearedSphere2D(s) 
    182197        Struct ResSmear_2D_AAOStruct &s 
    183198         
    184         Smear_2DModel_5(Sphere2D_noThread,s) 
    185 //      Smear_2DModel_20(Sphere2D_noThread,s) 
     199//// non-threaded, but generic calculation 
     200//// the last param is nord      
     201//      Smear_2DModel_PP(Sphere2D_noThread,s,10) 
     202 
     203//// the threaded version must be specifically written, since 
     204//// FUNCREF can't be passed into  a threaded calc (or structures) 
     205//// the last param is nord 
     206        SmearSphere2D_THR(s,10)          
     207 
    186208        return(0) 
    187209end 
     
    197219        WAVE qy = $(DF+str+"_qy") 
    198220        WAVE qz = $(DF+str+"_qz") 
    199         WAVE sigQx = $(DF+str+"_sigQx") 
    200         WAVE sigQy = $(DF+str+"_sigQy") 
     221        WAVE sQpl = $(DF+str+"_sQpl") 
     222        WAVE sQpp = $(DF+str+"_sQpp") 
    201223        WAVE shad = $(DF+str+"_fs") 
    202224         
     
    204226        WAVE s.coefW = coefW     
    205227        WAVE s.zw = resultW      
    206         WAVE s.qx = qx 
    207         WAVE s.qy = qy 
     228        WAVE s.xw[0] = qx 
     229        WAVE s.xw[1] = qy 
    208230        WAVE s.qz = qz 
    209         WAVE s.sigQx = sigQx 
    210         WAVE s.sigQy = sigQy 
     231        WAVE s.sQpl = sQpl 
     232        WAVE s.sQpp = sQpp 
    211233        WAVE s.fs = shad 
    212234         
     
    216238        return (0) 
    217239End 
     240 
     241 
     242 
     243 
     244// 
     245// this is the threaded version, that dispatches the calculation out to threads 
     246// 
     247// must be written specific to each 2D function 
     248//  
     249Function SmearSphere2D_THR(s,nord) 
     250        Struct ResSmear_2D_AAOStruct &s 
     251        Variable nord 
     252         
     253        String weightStr,zStr 
     254         
     255// create all of the necessary quadrature waves here - rather than inside a threadsafe function 
     256        switch(nord)     
     257                case 5:          
     258                        weightStr="gauss5wt" 
     259                        zStr="gauss5z" 
     260                        if (WaveExists($weightStr) == 0) 
     261                                Make/O/D/N=(nord) $weightStr,$zStr 
     262                                Make5GaussPoints($weightStr,$zStr)       
     263                        endif 
     264                        break                            
     265                case 10:                 
     266                        weightStr="gauss10wt" 
     267                        zStr="gauss10z" 
     268                        if (WaveExists($weightStr) == 0) 
     269                                Make/O/D/N=(nord) $weightStr,$zStr 
     270                                Make10GaussPoints($weightStr,$zStr)      
     271                        endif 
     272                        break                            
     273                case 20:                 
     274                        weightStr="gauss20wt" 
     275                        zStr="gauss20z" 
     276                        if (WaveExists($weightStr) == 0) 
     277                                Make/O/D/N=(nord) $weightStr,$zStr 
     278                                Make20GaussPoints($weightStr,$zStr)      
     279                        endif 
     280                        break 
     281                default:                                                         
     282                        Abort "Smear_2DModel_PP_Threaded called with invalid nord value"                                         
     283        endswitch 
     284         
     285        Wave/Z wt = $weightStr 
     286        Wave/Z xi = $zStr               // wave references to pass 
     287 
     288        Variable npt=numpnts(s.xw[0]) 
     289        Variable i,nthreads= ThreadProcessorCount 
     290        variable mt= ThreadGroupCreate(nthreads) 
     291 
     292        Variable t1=StopMSTimer(-2) 
     293         
     294        for(i=0;i<nthreads;i+=1) 
     295//              Print (i*npt/nthreads),((i+1)*npt/nthreads-1) 
     296                ThreadStart mt,i,SmearSphere2D_T(s.coefW,s.xw[0],s.xw[1],s.qz,s.sQpl,s.sQpp,s.fs,s.zw,wt,xi,(i*npt/nthreads),((i+1)*npt/nthreads-1),nord) 
     297        endfor 
     298 
     299        do 
     300                variable tgs= ThreadGroupWait(mt,100) 
     301        while( tgs != 0 ) 
     302 
     303        variable dummy= ThreadGroupRelease(mt) 
     304         
     305// comment out the threading + uncomment this for testing to make sure that the single thread works 
     306//      nThreads=1 
     307//      SmearSphere2D_T(s.coefW,s.xw[0],s.xw[1],s.qz,s.sQpl,s.sQpp,s.fs,s.zw,wt,xi,(i*npt/nthreads),((i+1)*npt/nthreads-1),nord) 
     308 
     309        Print "elapsed time = ",(StopMSTimer(-2) - t1)/1e6 
     310         
     311        return(0) 
     312end 
     313 
     314// 
     315// - worker function for threads of Sphere2D 
     316// 
     317ThreadSafe Function SmearSphere2D_T(coef,qxw,qyw,qzw,sxw,syw,fsw,zw,wt,xi,pt1,pt2,nord) 
     318        WAVE coef,qxw,qyw,qzw,sxw,syw,fsw,zw,wt,xi 
     319        Variable pt1,pt2,nord 
     320         
     321// now passed in.... 
     322//      Wave wt = $weightStr 
     323//      Wave xi = $zStr          
     324 
     325        Variable ii,jj,kk,num 
     326        Variable qx,qy,qz,qval,sx,sy,fs 
     327        Variable qy_pt,qx_pt,res_x,res_y,answer,sumIn,sumOut 
     328         
     329        Variable normFactor,phi,theta,maxSig,numStdDev=3 
     330         
     331/// keep these waves local 
     332        Make/O/D/N=(nord) fcnRet,xptW,res_tot,yptW 
     333         
     334        // now just loop over the points as specified 
     335         
     336        answer=0 
     337         
     338        Variable spl,spp,apl,app,bpl,bpp,phi_pt,qpl_pt 
     339 
     340        //loop over q-values 
     341        for(ii=pt1;ii<(pt2+1);ii+=1) 
     342                 
     343                qx = qxw[ii] 
     344                qy = qyw[ii] 
     345                qz = qzw[ii] 
     346                qval = sqrt(qx^2+qy^2+qz^2) 
     347                spl = sxw[ii] 
     348                spp = syw[ii] 
     349                fs = fsw[ii] 
     350                 
     351                normFactor = 2*pi*spl*spp 
     352                 
     353                phi = FindPhi(qx,qy) 
     354                 
     355                apl = -numStdDev*spl + qval             //parallel = q integration limits 
     356                bpl = numStdDev*spl + qval 
     357                app = -numStdDev*spp + phi              //perpendicular = phi integration limits 
     358                bpp = numStdDev*spp + phi 
     359                 
     360                //make sure the limits are reasonable. 
     361                if(apl < 0) 
     362                        apl = 0 
     363                endif 
     364                // do I need to specially handle limits when phi ~ 0? 
     365         
     366                 
     367                sumOut = 0 
     368                for(jj=0;jj<nord;jj+=1)         // call phi the "outer' 
     369                        phi_pt = (xi[jj]*(bpp-app)+app+bpp)/2 
     370 
     371                        sumIn=0 
     372                        for(kk=0;kk<nord;kk+=1)         //at phi, integrate over Qpl 
     373 
     374                                qpl_pt = (xi[kk]*(bpl-apl)+apl+bpl)/2 
     375                                 
     376                                FindQxQy(qpl_pt,phi_pt,qx_pt,qy_pt)             //find the corresponding QxQy to the Q,phi 
     377                                yPtw[kk] = qy_pt                                        //phi is the same in this loop, but qy is not 
     378                                xPtW[kk] = qx_pt                                        //qx is different here too, as we're varying Qpl 
     379                                 
     380                                res_tot[kk] = exp(-0.5*( (qpl_pt-qval)^2/spl/spl + (phi_pt-phi)^2/spp/spp ) ) 
     381                                res_tot[kk] /= normFactor 
     382//                              res_tot[kk] *= fs 
     383 
     384                        endfor 
     385                         
     386                        Sphere2D_noThread(coef,fcnRet,xptw,yptw)                        //fcn passed in is an AAO 
     387                         
     388                        //sumIn += wt[jj]*wt[kk]*res_tot*fcnRet[0] 
     389                        fcnRet *= wt[jj]*wt*res_tot 
     390                        // 
     391                        answer += (bpl-apl)/2.0*sum(fcnRet)             // 
     392                endfor 
     393 
     394                answer *= (bpp-app)/2.0 
     395                zw[ii] = answer 
     396        endfor 
     397         
     398        return(0) 
     399end 
     400 
     401 
  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Packages/Wrapper_v40.ipf

    r704 r708  
    10441044         
    10451045        if(yesReport) 
    1046                 String parStr=GetWavesDataFolder(cw,1)+ WaveList("*param*"+suffix, "", "TEXT:1," )              //this is *hopefully* one wave 
     1046                String parStr=GetWavesDataFolder(cw,1)+ WaveList("*param*"+suffix, "", "TEXT:1," )              // this is *hopefully* one wave 
    10471047                String topGraph= WinName(0,1)   //this is the topmost graph 
    10481048         
  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/SA_includes_v410.ipf

    r693 r708  
    2828 
    2929#include "DataSetHandling"                                      //added Nov 2009 AJJ 
     30#include "Smear_2D"                                             //for 2D resolution smearing, May 2010 
    3031 
    3132Menu "SANS Models" 
  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/WriteModelData_v40.ipf

    r607 r708  
    193193                PathInfo/S catPathName 
    194194                fullPath = DoSaveFileDialog("Save data as",fname=baseStr+".txt") 
     195                Print fullPath 
    195196                If(cmpstr(fullPath,"")==0) 
    196197                        //user cancel, don't write out a file 
     
    200201                //Print "dialog fullpath = ",fullpath 
    201202        Endif 
    202          
     203 
    203204        Open refnum as fullpath 
    204205         
Note: See TracChangeset for help on using the changeset viewer.