source: sasview/sansmodels/src/c_models/parameters.cpp @ 96277b9

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since 96277b9 was 0528cd6, checked in by Jae Cho <jhjcho@…>, 13 years ago

less number of points as default

  • Property mode set to 100644
File size: 11.5 KB
Line 
1/**
2        This software was developed by the University of Tennessee as part of the
3        Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
4        project funded by the US National Science Foundation.
5
6        If you use DANSE applications to do scientific research that leads to
7        publication, we ask that you acknowledge the use of the software with the
8        following sentence:
9
10        "This work benefited from DANSE software developed under NSF award DMR-0520547."
11
12        copyright 2008, University of Tennessee
13 */
14#include "parameters.hh"
15#include <stdio.h>
16#include <math.h>
17using namespace std;
18
19#if defined(_MSC_VER)
20#include "gamma_win.h"
21#endif
22
23/**
24 * TODO: normalize all dispersion weight lists
25 */
26
27
28/**
29 * Weight points
30 */
31WeightPoint :: WeightPoint() {
32        value = 0.0;
33        weight = 0.0;
34}
35
36WeightPoint :: WeightPoint(double v, double w) {
37        value = v;
38        weight = w;
39}
40
41/**
42 * Dispersion models
43 */
44DispersionModel :: DispersionModel() {
45        npts  = 1;
46        width = 0.0;
47};
48
49void DispersionModel :: accept_as_source(DispersionVisitor* visitor, void* from, void* to) {
50        visitor->dispersion_to_dict(from, to);
51}
52void DispersionModel :: accept_as_destination(DispersionVisitor* visitor, void* from, void* to) {
53        visitor->dispersion_from_dict(from, to);
54}
55void DispersionModel :: operator() (void *param, vector<WeightPoint> &weights){
56        // Check against zero width
57        if (width<=0) {
58                width = 0.0;
59                npts  = 1;
60        }
61
62        Parameter* par = (Parameter*)param;
63        double value = (*par)();
64        double sig;
65        if (npts<2) {
66                weights.insert(weights.end(), WeightPoint(value, 1.0));
67        } else {
68                for(int i=0; i<npts; i++) {
69
70                        if ((*par).has_min==false){
71                                // width = sigma for angles
72                                sig = width;
73                        }
74                        else{
75                                //width = polydispersity (=sigma/value) for length
76                                sig = width * value;
77                        }
78                        double val = value + sig * (1.0*double(i)/double(npts-1) - 0.5);
79                        if ( ((*par).has_min==false || val>(*par).min)
80                          && ((*par).has_max==false || val<(*par).max)  )
81                                weights.insert(weights.end(), WeightPoint(val, 1.0));
82                }
83        }
84}
85
86/**
87 * Method to set the weights
88 * Not implemented for this class
89 */
90void DispersionModel :: set_weights(int npoints, double* values, double* weights){}
91
92/**
93 * Gaussian dispersion
94 */
95
96GaussianDispersion :: GaussianDispersion() {
97        npts  = 35;
98        width = 0.0;
99        nsigmas = 3;
100};
101
102void GaussianDispersion :: accept_as_source(DispersionVisitor* visitor, void* from, void* to) {
103        visitor->gaussian_to_dict(from, to);
104}
105void GaussianDispersion :: accept_as_destination(DispersionVisitor* visitor, void* from, void* to) {
106        visitor->gaussian_from_dict(from, to);
107}
108
109double gaussian_weight(double mean, double sigma, double x) {
110        double vary, expo_value;
111    vary = x-mean;
112    expo_value = -vary*vary/(2.0*sigma*sigma);
113    //return 1.0;
114    return exp(expo_value);
115}
116
117/**
118 * Gaussian dispersion
119 * @param mean: mean value of the Gaussian
120 * @param sigma: standard deviation of the Gaussian
121 * @param x: value at which the Gaussian is evaluated
122 * @return: value of the Gaussian
123 */
124void GaussianDispersion :: operator() (void *param, vector<WeightPoint> &weights){
125        // Check against zero width
126        if (width<=0) {
127                width = 0.0;
128                npts  = 1;
129                nsigmas = 3;
130        }
131
132        Parameter* par = (Parameter*)param;
133        double value = (*par)();
134        double sig;
135        if (npts<2) {
136                weights.insert(weights.end(), WeightPoint(value, 1.0));
137        } else {
138                for(int i=0; i<npts; i++) {
139                        if ((*par).has_min==false){
140                                // width = sigma for angles
141                                sig = width;
142                        }
143                        else{
144                                //width = polydispersity (=sigma/value) for length
145                                sig = width * value;
146                        }
147                        // We cover n(nsigmas) times sigmas on each side of the mean
148                        double val = value + sig * (2.0*nsigmas*double(i)/double(npts-1) - nsigmas);
149                        if ( ((*par).has_min==false || val>(*par).min)
150                          && ((*par).has_max==false || val<(*par).max)  ) {
151                                double _w = gaussian_weight(value, sig, val);
152                                weights.insert(weights.end(), WeightPoint(val, _w));
153                        }
154                }
155        }
156}
157
158
159/**
160 * Flat dispersion
161 */
162
163RectangleDispersion :: RectangleDispersion() {
164        npts  = 35;
165        width = 0.0;
166        nsigmas = 1.73205;
167};
168
169void RectangleDispersion :: accept_as_source(DispersionVisitor* visitor, void* from, void* to) {
170        visitor->rectangle_to_dict(from, to);
171}
172void RectangleDispersion :: accept_as_destination(DispersionVisitor* visitor, void* from, void* to) {
173        visitor->rectangle_from_dict(from, to);
174}
175
176double rectangle_weight(double mean, double sigma, double x) {
177    double wid = fabs(sigma) * sqrt(3.0);
178    if (x>= (mean-wid) && x<=(mean+wid)){
179        return 1.0;
180    }
181    else{
182        return 0.0;
183    }
184}
185
186/**
187 * Flat dispersion
188 * @param mean: mean value
189 * @param sigma: half width of top hat function
190 * @param x: value at which the Flat is evaluated
191 * @return: value of the Flat
192 */
193void RectangleDispersion :: operator() (void *param, vector<WeightPoint> &weights){
194        // Check against zero width
195        if (width<=0) {
196                width = 0.0;
197                npts  = 1;
198                nsigmas = 1.73205;
199        }
200
201        Parameter* par = (Parameter*)param;
202        double value = (*par)();
203        double sig;
204        if (npts<2) {
205                weights.insert(weights.end(), WeightPoint(value, 1.0));
206        } else {
207                for(int i=0; i<npts; i++) {
208                        if ((*par).has_min==false){
209                                // width = sigma for angles
210                                sig = width;
211                        }
212                        else{
213                                //width = polydispersity (=sigma/value) for length
214                                sig = width * value;
215                        }
216                        // We cover n(nsigmas) times sigmas on each side of the mean
217                        double val = value + sig * (2.0*nsigmas*double(i)/double(npts-1) - nsigmas);
218                        if ( ((*par).has_min==false || val>(*par).min)
219                          && ((*par).has_max==false || val<(*par).max)  ) {
220                                double _w = rectangle_weight(value, sig, val);
221                                weights.insert(weights.end(), WeightPoint(val, _w));
222                        }
223                }
224        }
225}
226
227
228/**
229 * LogNormal dispersion
230 */
231
232LogNormalDispersion :: LogNormalDispersion() {
233        npts  = 80;
234        width = 0.0;
235        nsigmas = 8.0;
236};
237
238void LogNormalDispersion :: accept_as_source(DispersionVisitor* visitor, void* from, void* to) {
239        visitor->lognormal_to_dict(from, to);
240}
241void LogNormalDispersion :: accept_as_destination(DispersionVisitor* visitor, void* from, void* to) {
242        visitor->lognormal_from_dict(from, to);
243}
244
245double lognormal_weight(double mean, double sigma, double x) {
246
247        double sigma2 = pow(sigma, 2.0);
248        return 1.0/(x*sigma) * exp( -pow((log(x) - mean), 2.0) / (2.0*sigma2));
249
250}
251
252/**
253 * Lognormal dispersion
254 * @param mean: mean value of the LogNormal
255 * @param sigma: standard deviation of the LogNormal
256 * @param x: value at which the LogNormal is evaluated
257 * @return: value of the LogNormal
258 */
259void LogNormalDispersion :: operator() (void *param, vector<WeightPoint> &weights){
260        // Check against zero width
261        if (width<=0) {
262                width = 0.0;
263                npts  = 1;
264                nsigmas = 8.0;
265        }
266
267        Parameter* par = (Parameter*)param;
268        double value = (*par)();
269        double sig;
270        if (npts<2) {
271                weights.insert(weights.end(), WeightPoint(value, 1.0));
272        } else {
273                for(int i=0; i<npts; i++) {
274                        // Note that the definition of sigma is different from Gaussian
275                        if ((*par).has_min==false){
276                                // sig  for angles
277                                sig = width;
278                        }
279                        else{
280                                // by lognormal definition, PD is same as sigma
281                                sig = width * value;
282                        }
283                       
284                        // We cover n(nsigmas) times sigmas on each side of the mean
285                        //constant bin in real space
286                        double val = value + sig * (2.0*nsigmas*double(i)/double(npts-1) - nsigmas);
287                        // sigma in the lognormal function is in ln(R) space, thus needs converting
288                        sig = fabs(sig/value); 
289                        if ( ((*par).has_min==false || val>(*par).min)
290                          && ((*par).has_max==false || val<(*par).max)  ) {
291                                double _w = lognormal_weight(log(value), sig, val);
292                                weights.insert(weights.end(), WeightPoint(val, _w));
293                                //printf("%g \t %g\n",val,_w);
294
295                        }
296                }
297        }
298}
299
300
301
302/**
303 * Schulz dispersion
304 */
305
306SchulzDispersion :: SchulzDispersion() {
307        npts  = 80;
308        width = 0.0;
309        nsigmas = 8.0;
310};
311
312void SchulzDispersion :: accept_as_source(DispersionVisitor* visitor, void* from, void* to) {
313        visitor->schulz_to_dict(from, to);
314}
315void SchulzDispersion :: accept_as_destination(DispersionVisitor* visitor, void* from, void* to) {
316        visitor->schulz_from_dict(from, to);
317}
318
319double schulz_weight(double mean, double sigma, double x) {
320    double z = pow(mean/ sigma, 2.0)-1.0;
321        double R= x/mean;
322        double zz= z+1.0;
323        double expo;
324        expo = zz*log(zz)+z*log(R)-R*zz-log(mean)-lgamma(zz);
325        return  exp(expo);
326}
327
328/**
329 * Schulz dispersion
330 * @param mean: mean value of the Schulz
331 * @param sigma: standard deviation of the Schulz
332 * @param x: value at which the Schulz is evaluated
333 * @return: value of the Schulz
334 */
335void SchulzDispersion :: operator() (void *param, vector<WeightPoint> &weights){
336        // Check against zero width
337        if (width<=0) {
338                width = 0.0;
339                npts  = 1;
340                nsigmas = 8.0;
341        }
342
343        Parameter* par = (Parameter*)param;
344        double value = (*par)();
345        double sig;
346        if (npts<2) {
347                weights.insert(weights.end(), WeightPoint(value, 1.0));
348        } else {
349                for(int i=0; i<npts; i++) {
350                        if ((*par).has_min==false){
351                                // width = sigma for angles
352                                sig = width;
353                        }
354                        else{
355                                //width = polydispersity (=sigma/value) for length
356                                sig = width * value;
357                        }
358                        // We cover n(nsigmas) times sigmas on each side of the mean
359                        double val = value + sig * (2.0*nsigmas*double(i)/double(npts-1) - nsigmas);
360
361                        if ( ((*par).has_min==false || val>(*par).min)
362                          && ((*par).has_max==false || val<(*par).max)  ) {
363                                double _w = schulz_weight(value, sig, val);
364                                weights.insert(weights.end(), WeightPoint(val, _w));
365                        }
366                }
367        }
368}
369
370
371
372
373/**
374 * Array dispersion based on input arrays
375 */
376
377void ArrayDispersion :: accept_as_source(DispersionVisitor* visitor, void* from, void* to) {
378        visitor->array_to_dict(from, to);
379}
380void ArrayDispersion :: accept_as_destination(DispersionVisitor* visitor, void* from, void* to) {
381        visitor->array_from_dict(from, to);
382}
383
384/**
385 * Method to get the weights
386 */
387void ArrayDispersion :: operator() (void *param, vector<WeightPoint> &weights) {
388        Parameter* par = (Parameter*)param;
389        double value = (*par)();
390
391        if (npts<2) {
392                weights.insert(weights.end(), WeightPoint(value, 1.0));
393        } else {
394                for(int i=0; i<npts; i++) {
395                        double val = _values[i]; //+ value;  //ToDo: Talk to Paul and put back the 'value'.
396
397                        if ( ((*par).has_min==false || val>(*par).min)
398                          && ((*par).has_max==false || val<(*par).max)  )
399                                weights.insert(weights.end(), WeightPoint(val, _weights[i]));
400                }
401        }
402}
403/**
404 * Method to set the weights
405 */
406void ArrayDispersion :: set_weights(int npoints, double* values, double* weights){
407        npts = npoints;
408        _values = values;
409        _weights = weights;
410}
411
412
413/**
414 * Parameters
415 */
416Parameter :: Parameter() {
417        value = 0;
418        min   = 0.0;
419        max   = 0.0;
420        has_min = false;
421        has_max = false;
422        has_dispersion = false;
423        dispersion = new GaussianDispersion();
424}
425
426Parameter :: Parameter(double _value) {
427        value = _value;
428        min   = 0.0;
429        max   = 0.0;
430        has_min = false;
431        has_max = false;
432        has_dispersion = false;
433        dispersion = new GaussianDispersion();
434}
435
436Parameter :: Parameter(double _value, bool disp) {
437        value = _value;
438        min   = 0.0;
439        max   = 0.0;
440        has_min = false;
441        has_max = false;
442        has_dispersion = disp;
443        dispersion = new GaussianDispersion();
444}
445
446void Parameter :: get_weights(vector<WeightPoint> &weights) {
447        (*dispersion)((void*)this, weights);
448}
449
450void Parameter :: set_min(double value) {
451        has_min = true;
452        min = value;
453}
454
455void Parameter :: set_max(double value) {
456        has_max = true;
457        max = value;
458}
459
460double Parameter :: operator()() {
461        return value;
462}
463
464double Parameter :: operator=(double _value){
465        value = _value;
466        return value;
467}
Note: See TracBrowser for help on using the repository browser.