source: sasview/sansmodels/igor_wrapper/src/c_disperser.c @ c21d6b9

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 c21d6b9 was 25a60dc1, checked in by Jae Cho <jhjcho@…>, 13 years ago

moving a folder

  • Property mode set to 100644
File size: 5.9 KB
Line 
1/**
2 * Straight C disperser
3 *
4 * This code was written as part of the DANSE project
5 * http://danse.us/trac/sans/
6 *
7 * Copyright 2007: University of Tennessee, for the DANSE project
8 */
9#include "math.h"
10#include <stdio.h>
11#include <stdlib.h>
12
13/**
14 * Weight distribution to give to each point in the dispersion
15 * The distribution is a Gaussian with
16 *
17 * @param mean: mean value of the Gaussian
18 * @param sigma: sigma of the Gaussian
19 * @param x: point to evaluate at
20 * @return: weight value
21 *
22 */
23double c_disperser_weight(double mean, double sigma, double x) {
24        double vary, expo_value;
25    vary = x-mean;
26    expo_value = -vary*vary/(2*sigma*sigma);
27    return exp(expo_value);
28}
29
30/**
31 * Function to apply dispersion to a list of parameters.
32 *
33 * This function is re-entrant. It should be called with iPar=0.
34 * It will then call itself with increasing values for iPar until
35 * all parameters to be dispersed have been dealt with.
36 *
37 *
38 * @param eval: pointer to the function used to evaluate the model at
39 *                              a particular point.
40 * @param dp: complete array of parameter values for the model.
41 * @param n_pars: number of parameters to apply dispersion to.
42 * @param idList: list of parameter indices for the parameters to apply
43 *                              dispersion to. For a given parameter, its index is the
44 *                              index of its position in the parameter vector of the model
45 *                              function.
46 * @param sigmaList: list of sigma values for the parameters to apply
47 *                              dispersion to.
48 * @param centers: list of mean values for the parameters to apply
49 *                              dispersion to.
50 * @param n_pts: number of points to use when applying dispersion.
51 * @param q: q-value to evaluate the model at.
52 * @param phi: angle of the q-vector with the q_x axis.
53 * @param iPar: index of the parameter to apply dispersion to (should
54 *                              be 0 when called by the user).
55 *
56 */
57double c_disperseParam( double (*eval)(double[], double, double), double dp[], int n_pars, 
58                                        int *idList, double *sigmaList, double *centers, 
59                                        int n_pts, double q, double phi, int iPar ) {
60        double min_value, max_value;
61        double step;
62        double prev_value;
63        double value_sum;
64        double gauss_sum;
65        double gauss_value;
66        double func_value;
67        double error_sys;
68        double value;
69        int n_sigma;
70        int i;
71        // Number of std variations to average over     
72        n_sigma = 2;   
73    if( iPar < n_pars ) {
74           
75                // Average over Gaussian distribution (2 sigmas)
76        value_sum = 0.0;
77        gauss_sum = 0.0;
78           
79        // Average over 4 sigmas wide
80        min_value = centers[iPar] - n_sigma*sigmaList[iPar];
81        max_value = centers[iPar] + n_sigma*sigmaList[iPar];
82           
83        // Calculate step size
84        step = (max_value - min_value)/(n_pts-1);
85           
86        // If we are not changing the parameter, just return the
87        // value of the function
88        if (step == 0.0) {
89            return c_disperseParam(eval, dp, n_pars, idList, sigmaList, 
90                                                        centers, n_pts, q, phi, iPar+1);
91        }
92           
93        // Compute average
94        prev_value = 0.0;
95        error_sys  = 0.0;
96        for( i=0; i<n_pts; i++ ) {
97            // Set the parameter value           
98            value = min_value + (double)i*step;
99            dp[idList[iPar]] = value;
100               
101                        gauss_value = c_disperser_weight(centers[iPar], sigmaList[iPar], value);
102            func_value = c_disperseParam(eval, dp, n_pars, idList, sigmaList, 
103                                                        centers, n_pts, q, phi, iPar+1);
104
105            value_sum += gauss_value * func_value;
106                gauss_sum += gauss_value;       
107        }   
108        return value_sum/gauss_sum;
109       
110    } else {
111         return (*eval)(dp, q, phi);
112    }
113       
114}
115
116/**
117 * Function to add dispersion to a model.
118 * The dispersion is Gaussian around the value of given parameters.
119 *
120 * @param eval: pointer to the function used to evaluate the model at
121 *                              a particular point.
122 * @param n_pars: number of parameters to apply dispersion to.
123 * @param idList: list of parameter indices for the parameters to apply
124 *                              dispersion to. For a given parameter, its index is the
125 *                              index of its position in the parameter vector of the model
126 *                              function.
127 * @param sigmaList: list of sigma values for the parameters to apply
128 *                              dispersion to.
129 * @param n_pts: number of points to use when applying dispersion.
130 * @param q: q-value to evaluate the model at.
131 * @param phi: angle of the q-vector with the q_x axis.
132 *
133 */
134double c_disperser( double (*eval)(double[], double, double), double dp[], int n_pars, 
135                                        int *idList, double *sigmaList, int n_pts, double q, double phi ) {
136        double *centers;
137        double value;
138        int i;
139       
140        // Allocate centers array
141        if( n_pars > 0 ) {
142            centers = (double *)malloc(n_pars * sizeof(double));
143                if(centers==NULL) {
144                    printf("c_disperser could not allocate memory\n");
145                        return 0.0;
146                }
147        }
148       
149        // Store current values in centers array
150        for(i=0; i<n_pars; i++) {
151                centers[i] = dp[idList[i]];     
152        }
153       
154       
155        if( n_pars > 0 ) {
156                value = c_disperseParam(eval, dp, n_pars, idList, sigmaList, centers, n_pts, q, phi, 0);
157    } else {
158        value = (*eval)(dp, q, phi);
159        }
160        free(centers);
161        return value;
162}
163
164/**
165 *
166 * Angles are in radian.
167 *
168 *
169 */
170double weight_dispersion( double (*eval)(double[], double, double), 
171                                                  double *par_values, double *weight_values, 
172                                                  int npts, int i_par,
173                                                  double dp[], double q, double phi ) {
174        int i;
175        double value;
176        double norma;
177                               
178        value = 0.0;
179        norma = 0.0;           
180
181        // If we have an empty array of points, just
182        // evaluate the function
183        if(npts == 0) {
184                return (*eval)(dp, q, phi);
185        } else {
186                for(i=0; i<npts; i++) {
187                        dp[i_par] = par_values[i];
188                        value += weight_values[i] * (*eval)(dp, q, phi);
189                        //dp[i_par] = -par_values[i];
190                        //value += weight_values[i] * (*eval)(dp, q, phi);
191                       
192                        norma += weight_values[i];
193                }
194        }
195        return value/norma/2.0;
196                                                       
197}
Note: See TracBrowser for help on using the repository browser.