source: sasview/sansmodels/src/c_models/parameters.cpp @ acd0fd10

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 acd0fd10 was 67424cd, checked in by Mathieu Doucet <doucetm@…>, 13 years ago

Reorganizing models in preparation of cpp cleanup

  • Property mode set to 100644
File size: 11.5 KB
RevLine 
[fca6936]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
[83bf68e1]19#if defined(_MSC_VER)
20#include "gamma_win.h"
21#endif
22
[fca6936]23/**
24 * TODO: normalize all dispersion weight lists
25 */
26
27
28/**
[836fe6e]29 * Weight points
[fca6936]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/**
[836fe6e]42 * Dispersion models
[fca6936]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)();
[59b9b675]64        double sig;
[fca6936]65        if (npts<2) {
66                weights.insert(weights.end(), WeightPoint(value, 1.0));
67        } else {
68                for(int i=0; i<npts; i++) {
69
[59b9b675]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);
[fca6936]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/**
[836fe6e]93 * Gaussian dispersion
[fca6936]94 */
95
96GaussianDispersion :: GaussianDispersion() {
[a5e14410]97        npts  = 100;
[fca6936]98        width = 0.0;
[a5e14410]99        nsigmas = 10;
[fca6936]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;
[8dc02d8b]112    expo_value = -vary*vary/(2.0*sigma*sigma);
[fca6936]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;
[a5e14410]129                nsigmas = 10;
[fca6936]130        }
131
132        Parameter* par = (Parameter*)param;
133        double value = (*par)();
[59b9b675]134        double sig;
[fca6936]135        if (npts<2) {
136                weights.insert(weights.end(), WeightPoint(value, 1.0));
137        } else {
138                for(int i=0; i<npts; i++) {
[59b9b675]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                        }
[fcd8a80e]147                        // We cover n(nsigmas) times sigmas on each side of the mean
[59b9b675]148                        double val = value + sig * (2.0*nsigmas*double(i)/double(npts-1) - nsigmas);
[fca6936]149                        if ( ((*par).has_min==false || val>(*par).min)
150                          && ((*par).has_max==false || val<(*par).max)  ) {
[59b9b675]151                                double _w = gaussian_weight(value, sig, val);
[fca6936]152                                weights.insert(weights.end(), WeightPoint(val, _w));
153                        }
154                }
155        }
156}
157
[eba9885]158
159/**
[8dc02d8b]160 * Flat dispersion
161 */
162
163RectangleDispersion :: RectangleDispersion() {
[a5e14410]164        npts  = 100;
[8dc02d8b]165        width = 0.0;
[0ad5703]166        nsigmas = 1.73205;
[8dc02d8b]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) {
[0ad5703]177    double wid = fabs(sigma) * sqrt(3.0);
178    if (x>= (mean-wid) && x<=(mean+wid)){
[8dc02d8b]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;
[0ad5703]198                nsigmas = 1.73205;
[8dc02d8b]199        }
200
201        Parameter* par = (Parameter*)param;
202        double value = (*par)();
[59b9b675]203        double sig;
[8dc02d8b]204        if (npts<2) {
205                weights.insert(weights.end(), WeightPoint(value, 1.0));
206        } else {
207                for(int i=0; i<npts; i++) {
[59b9b675]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                        }
[8dc02d8b]216                        // We cover n(nsigmas) times sigmas on each side of the mean
[59b9b675]217                        double val = value + sig * (2.0*nsigmas*double(i)/double(npts-1) - nsigmas);
[8dc02d8b]218                        if ( ((*par).has_min==false || val>(*par).min)
219                          && ((*par).has_max==false || val<(*par).max)  ) {
[59b9b675]220                                double _w = rectangle_weight(value, sig, val);
[8dc02d8b]221                                weights.insert(weights.end(), WeightPoint(val, _w));
222                        }
223                }
224        }
225}
226
227
228/**
[eba9885]229 * LogNormal dispersion
230 */
231
232LogNormalDispersion :: LogNormalDispersion() {
[a5e14410]233        npts  = 100;
[eba9885]234        width = 0.0;
[a5e14410]235        nsigmas = 10.0;
[eba9885]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) {
[c5607fa]246
[1d78e4b]247        double sigma2 = pow(sigma, 2.0);
[a5e14410]248        return 1.0/(x*sigma) * exp( -pow((log(x) - mean), 2.0) / (2.0*sigma2));
[c5607fa]249
[eba9885]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;
[a5e14410]264                nsigmas = 10.0;
[eba9885]265        }
266
267        Parameter* par = (Parameter*)param;
268        double value = (*par)();
[59b9b675]269        double sig;
[eba9885]270        if (npts<2) {
271                weights.insert(weights.end(), WeightPoint(value, 1.0));
272        } else {
273                for(int i=0; i<npts; i++) {
[0ad5703]274                        // Note that the definition of sigma is different from Gaussian
[59b9b675]275                        if ((*par).has_min==false){
[0ad5703]276                                // sig  for angles
[a5e14410]277                                sig = width;
[59b9b675]278                        }
279                        else{
[0ad5703]280                                // by lognormal definition, PD is same as sigma
[a5e14410]281                                sig = width * value;
[59b9b675]282                        }
[a5e14410]283                       
[eba9885]284                        // We cover n(nsigmas) times sigmas on each side of the mean
[a5e14410]285                        //constant bin in real space
[59b9b675]286                        double val = value + sig * (2.0*nsigmas*double(i)/double(npts-1) - nsigmas);
[a5e14410]287                        // sigma in the lognormal function is in ln(R) space, thus needs converting
288                        sig = fabs(sig/value); 
[eba9885]289                        if ( ((*par).has_min==false || val>(*par).min)
290                          && ((*par).has_max==false || val<(*par).max)  ) {
[a5e14410]291                                double _w = lognormal_weight(log(value), sig, val);
[eba9885]292                                weights.insert(weights.end(), WeightPoint(val, _w));
[a5e14410]293                                //printf("%g \t %g\n",val,_w);
294
[eba9885]295                        }
296                }
297        }
298}
299
300
301
302/**
303 * Schulz dispersion
304 */
305
306SchulzDispersion :: SchulzDispersion() {
[a5e14410]307        npts  = 100;
[eba9885]308        width = 0.0;
[a5e14410]309        nsigmas = 10.0;
[eba9885]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) {
[1d78e4b]320    double z = pow(mean/ sigma, 2.0)-1.0;
[eba9885]321        double R= x/mean;
[1d78e4b]322        double zz= z+1.0;
[c5607fa]323        double expo;
324        expo = zz*log(zz)+z*log(R)-R*zz-log(mean)-lgamma(zz);
325        return  exp(expo);
[eba9885]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;
[a5e14410]340                nsigmas = 10.0;
[eba9885]341        }
342
343        Parameter* par = (Parameter*)param;
344        double value = (*par)();
[59b9b675]345        double sig;
[eba9885]346        if (npts<2) {
347                weights.insert(weights.end(), WeightPoint(value, 1.0));
348        } else {
349                for(int i=0; i<npts; i++) {
[59b9b675]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                        }
[eba9885]358                        // We cover n(nsigmas) times sigmas on each side of the mean
[59b9b675]359                        double val = value + sig * (2.0*nsigmas*double(i)/double(npts-1) - nsigmas);
[eba9885]360
361                        if ( ((*par).has_min==false || val>(*par).min)
362                          && ((*par).has_max==false || val<(*par).max)  ) {
[59b9b675]363                                double _w = schulz_weight(value, sig, val);
[eba9885]364                                weights.insert(weights.end(), WeightPoint(val, _w));
365                        }
366                }
367        }
368}
369
370
371
372
[fca6936]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
[07da749]391        if (npts<2) {
392                weights.insert(weights.end(), WeightPoint(value, 1.0));
393        } else {
[fca6936]394                for(int i=0; i<npts; i++) {
[0fff338f]395                        double val = _values[i]; //+ value;  //ToDo: Talk to Paul and put back the 'value'.
[59b9b675]396
[0fff338f]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]));
[fca6936]400                }
[07da749]401        }
[fca6936]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/**
[836fe6e]414 * Parameters
[fca6936]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;
[34c2649]466        return value;
[fca6936]467}
Note: See TracBrowser for help on using the repository browser.