source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Packages/GeneticOptimization/NCNR_GenFitUtils.ipf @ 705

Last change on this file since 705 was 705, checked in by srkline, 13 years ago

Fix conditional compile error if GenCurveFit? not installed

New version number 7.02 for Summer School release

File size: 9.6 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion=6.1
3
4// genetic optimization, uses XOP supplied by Andy Nelson, ANSTO
5// http://www.igorexchange.com/project/gencurvefit
6
7
8
9// TO USE GENETIC OPTIMIZATION
10//      Variable/G root:Packages:NIST:gUseGenCurveFit = 0                       //set to 1 to use genetic optimization
11// -- this flag is set by a menu item SANS Models->Packages->Enable Genetic Optimization
12//
13// ** it is currently only available for single fits - not for global fits (even the simple ones)
14// incorporating into WM's global fit would be u-g-l-y
15//
16//
17//  to complete the addition of genetic curve fitting:
18//
19// X add a check to make sure that the XOP is installed before the switch is set (force 0 if not present)
20// X add a global variable as a switch
21// X parse the limits. All must be filled in. If not, fit aborts and asks for user input to fill them in. Limits
22//   on fixed parameters are +/- 1% of the value. Has no effect on fit (duh), but filled in for completeness.
23// X create a mask wave for use when cursors are selected. [a,b] subranges can't be used w/Andy's XOP
24// X fitYw not behaving correctly - probably need to hand-trim
25// X use the switch in the wrapper, building the structure and function call as necessary
26// X be sure that the smeared function name is printed out as it's used, not the generic wrapper
27// X test for speed. the smeared fit is especially SLOW, WAY slower than it should be...AAO vs. point function??
28// X decide which options to add (/N /DUMP /TOL, robust fitting, etc - see Andy's list
29// X add the XOP to the distribution, or instructions (better to NOT bundle...)
30// X odd bug where first "fit" fails on AppendToGraph as GenCurveFit is started. May need to disable /D flag
31//
32// -- need to pass back the chi-squared and number of points. "V_" globals don't appear
33//
34// for the speed test. try writing my own wrapper for an unsmeared calculation, and see if it's still dog-slow.
35// --- NOPE, the wrapper is not a problem, tested this, and they seem to be both the same speed, each only about 5 seconds.
36// -- timing results for fitting apoferritin:
37//
38// L-M method, unsmeared: 0.95 s
39// L-M method, smeared: 5.85 s
40//
41// Genetic method, unsmeared: 7.32 s (2199 function evaluations) = 0.0033 sec/eval
42// Genetic method, smeared: 416 s (2830 function evaluations x20 for smearing) = 0.0074 sec/eval !
43//
44// -- even correcting the Gen-smeared fit for the 20pt quadrature, calculations there are 2x slower than
45// the same unsmeared fit. Why?? Both converge in the same number of iterations (~30 to 40). The number of function
46// evaluations is different, due to the random start and random mutations. Setting the seed is unlikely
47// to change the results.
48
49Static Constant kGenOp_tol=0.001
50
51
52Function Init_GenOp()
53        if(!exists("GenCurveFit"))
54                DoAlert 1,"The genetic optimiztion XOP is not installed. Do you want to open the web page to download the installer?"
55                if(V_flag == 1)
56                        BrowseURL "http://www.igorexchange.com/project/gencurvefit"
57                endif
58        else
59                DoAlert 0,"Genetic Optimization has been enabled."
60                Variable/G root:Packages:NIST:gUseGenCurveFit = 1                       //set to 1 to use genetic optimization
61        endif
62        BuildMenu "SANS Models"
63End
64
65// uncheck the flag, and change the menu
66Function UnSet_GenOp()
67        DoAlert 0,"Genetic Optimization has been disabled"
68        Variable/G root:Packages:NIST:gUseGenCurveFit = 0                       //set to 1 to use genetic optimization
69        BuildMenu "SANS Models"
70End
71
72Function/S GenOpFlagEnable()
73        Variable flag = NumVarOrDefault("root:Packages:NIST:gUseGenCurveFit",0)
74//      NVAR/Z flag = root:Packages:NIST:gUseGenCurveFit
75//      if(!NVAR_Exists(flag))
76//              return("")              //to catch the initial menu build that occurs before globals are created
77//      endif
78       
79        if(flag)
80                return("!"+num2char(18) + " ")
81        else
82                return("")
83        endif
84        BuildMenu "SANS Models"
85end
86// this is the default selection if nvar does not exist
87Function/S GenOpFlagDisable()
88        Variable flag = NumVarOrDefault("root:Packages:NIST:gUseGenCurveFit",0)
89
90//      NVAR/Z flag = root:Packages:NIST:gUseGenCurveFit
91//      if(!NVAR_Exists(flag))
92//              Variable/G root:Packages:NIST:gUseGenCurveFit = 0
93//      endif
94       
95        if(!flag)
96                return("!"+num2char(18) + " ")
97        else
98                return("")
99        endif
100        BuildMenu "SANS Models"
101end
102
103
104
105// the structure must be named fitFuncStruct, or the XOP will report an error and not run
106//
107Structure fitFuncStruct   
108        Wave w
109        wave y
110        wave x[50]
111        int16 numVarMD
112        wave ffsWaves[50]
113        wave ffsTextWaves[10]
114        variable ffsvar[5]
115        string ffsstr[5]
116        nvar ffsnvars[5]
117        svar ffssvars[5]
118        funcref SANSModelAAO_proto ffsfuncrefs[10]              //this is an AAO format
119        uint32 ffsversion    // Structure version.
120EndStructure
121
122Function GeneticFit_SmearedModel(s) : FitFunc
123        Struct fitFuncStruct &s
124
125        FUNCREF SANSModelSTRUCT_proto foo=$s.ffsstr[0]
126        NVAR num = root:num_evals
127       
128        STRUCT ResSmearAAOStruct fs
129        WAVE fs.coefW = s.w     
130        WAVE fs.yW = s.y
131        WAVE fs.xW = s.x[0]
132        WAVE fs.resW = s.ffsWaves[0]
133       
134        Variable err
135        err = foo(fs)           //this is the smeared STRUCT function
136               
137        num += 1
138        return(0)
139End
140
141Function GeneticFit_UnSmearedModel(s) : FitFunc
142        Struct fitFuncStruct &s
143
144        FUNCREF SANSModelAAO_proto foo=$s.ffsstr[0]
145        NVAR num = root:num_evals
146
147        foo(s.w,s.y,s.x[0])                     //this is the unsmeared function
148       
149        num += 1
150       
151        return(0)
152End
153
154// need to pass back the chi-squared and number of points. "V_" globals don't appear
155//
156// function will only compile if the GenCurveFit XOP is installed - else it will return zero
157//
158Function DoGenCurveFit(useRes,useCursors,sw,fitYw,fs,funcStr,holdStr,val,lolim,hilim,pt1,pt2)
159        Variable useRes,useCursors
160        WAVE sw,fitYw
161        STRUCT ResSmearAAOStruct &fs
162        String &funcStr,holdStr
163        Variable &val
164        WAVE/T lolim,hilim
165        Variable pt1,pt2                //already sorted if cursors are needed
166
167#if exists("GenCurveFit")
168
169        // initialise the structure you will use
170        struct fitFuncStruct bar
171
172        // we must set the version of the structure (currently 1000)
173        bar.ffsversion = 1000
174       
175        // numVarMD is the number of dependent variables you are fitting
176        // this must be correct, or Gencurvefit won't run.
177        bar.numVarMD=1         
178
179        // fill in the details and waves that I need
180        bar.ffsstr[0] = funcStr //generate the reference as needed for smeared or unsmeared
181        WAVE bar.w = fs.coefW
182        WAVE bar.y =fs.yW
183        WAVE bar.x[0] = fs.xW
184        WAVE/Z bar.ffsWaves[0] = fs.resW                //will not exist for 3-column data sets
185       
186
187        //need to parse limits, or make up some defaults
188        // limits is (n,2)
189        Variable nPnts = numpnts(fs.coefW),i,multip=1.01,isUSANS=0,tol=0.01
190        Make/O/D/N=(nPnts,2) limits
191        Wave limits=limits
192        for (i=0; i < nPnts; i += 1)
193       
194                if (strlen(lolim[i]) > 0)
195                        limits[i][0] = str2num(lolim[i])
196                else
197                        if(cmpstr(holdStr[i],"0")==0)           //no limit, not held
198                                Abort "You must enter low and high coefficient limits for all free parameters"
199                        else
200                                limits[i][0] = fs.coefW[i]/multip               //fixed parameter, just stick something in
201                        endif
202                endif
203               
204                if (strlen(hilim[i]) > 0)
205                        limits[i][1] = str2num(hilim[i])
206                else
207                        if(cmpstr(holdStr[i],"0")==0)           //no limit, not held
208                                Abort "You must enter low and high coefficient limits for all free parameters"
209                        else   
210                                limits[i][1] = fs.coefW[i] * multip
211                        endif
212                endif
213        endfor
214
215        if (dimsize(fs.resW,1) > 4)
216                isUSANS = 1
217        endif
218       
219//      generate a mask wave if needed (1=use, 0=mask)
220// currently, the mask is not used, since smeared USANS is not handled correctly
221// temporary, trimmed data sets are used instead
222        if(useCursors)
223                npnts = numpnts(fs.xW)
224//              Make/O/D/N=(npnts) GenOpMask
225//              Wave GenOpMask=GenOpMask
226//              GenOpMask = 1
227//              for(i=0;i<pt1;i+=1)
228//                      GenOpMask[i] = 0
229//              endfor
230//              for(i=pt2;i<npnts-1;i+=1)
231//                      GenOpMask[i] = 0
232//              endfor
233
234                Redimension/N=(pt2-pt1+1) fitYw
235               
236                Make/O/D/N=(pt2-pt1+1) trimY,trimX,trimS
237                WAVE trimY=trimY
238                trimY = fs.yW[p+pt1]
239                WAVE trimX=trimX
240                trimX = fs.xW[p+pt1]
241                WAVE trimS=trimS
242                trimS = sw[p+pt1]
243                //trim all of the waves, don't use the mask
244                WAVE bar.y = trimY
245                WAVE bar.x[0] = trimX
246               
247        endif   
248
249        Variable t0 = stopMStimer(-2)
250
251
252       
253        Variable/G root:num_evals               //for my own tracking
254        NVAR num = root:num_evals
255        num=0
256
257       
258        // other useful flags:  /N /DUMP /TOL /D=fitYw
259        //  /OPT=1 seems to make no difference whether it's used or not
260       
261        // append the fit
262        //do this only because GenCurveFit tries to append too quickly?
263        String traces=TraceNameList("", ";", 1 )                //"" as first parameter == look on the target graph
264        if(strsearch(traces,"FitYw",0) == -1)
265                if(useCursors)
266                        AppendtoGraph fitYw vs trimX
267                else
268                        AppendToGraph FitYw vs fs.xw
269                endif
270        else
271                RemoveFromGraph FitYw
272                if(useCursors)
273                        AppendtoGraph fitYw vs trimX
274                else
275                        AppendToGraph FitYw vs fs.xw
276                endif
277        endif
278        ModifyGraph lsize(FitYw)=2,rgb(FitYw)=(0,0,0)
279
280        do
281                       
282                if(useRes && useCursors)
283                        GenCurveFit/MAT /STRC=bar /X=bar.x[0] /I=1 /TOL=(kGenOp_tol) /W=trimS /D=fitYw GeneticFit_SmearedModel,bar.y,bar.w,holdStr,limits
284                        break
285                endif
286                if(useRes)
287                        GenCurveFit/MAT /STRC=bar /X=bar.x[0] /I=1 /TOL=(kGenOp_tol) /W=sw /D=fitYw GeneticFit_SmearedModel,bar.y,bar.w,holdStr,limits
288                        break
289                endif
290               
291                //no resolution
292                if(!useRes && useCursors)
293                        GenCurveFit/MAT /STRC=bar /X=bar.x[0] /I=1 /TOL=(kGenOp_tol) /W=trimS /D=fitYw GeneticFit_UnSmearedModel,bar.y,bar.w,holdStr,limits
294                        break
295                endif
296                if(!useRes)
297                        GenCurveFit/MAT /STRC=bar /X=bar.x[0] /I=1 /TOL=(kGenOp_tol) /W=sw /D=fitYw GeneticFit_UnSmearedModel,bar.y,bar.w,holdStr,limits
298                        break
299                endif
300               
301        while(0)       
302       
303//      NVAR V_chisq = V_chisq
304//      NVAR V_npnts = V_npnts
305//      NVAR V_fitIters = V_fitIters
306        WAVE/Z W_sigma = W_sigma
307       
308        val = V_npnts           //return this as a parameter
309       
310        t0 = (stopMSTimer(-2) - t0)*1e-6
311        Printf  "fit time = %g seconds\r\r",t0
312        Print W_sigma   
313        Print "number of iterations = ",V_fitIters
314        Print "number of function evaluations = ",num
315        Print "Chi-squared = ",V_chisq
316       
317        return(V_chisq)
318#else
319        return(0)
320#endif 
321end
Note: See TracBrowser for help on using the repository browser.