~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

Linux Cross Reference
Tina4/src/math/covar/covar.c

Version: ~
Architecture: ~ [ i386 ] ~ [ alpha ] ~ [ m68k ] ~ [ mips ] ~ [ ppc ] ~ [ sparc ] ~ [ sparc64 ] ~

  1 /* @(#)Jacobian & evaluate the inverse covariance matrix from the
  2  * jacobian */
  3 
  4 #include <math.h>
  5 #include <tina/sys.h>
  6 #include <tina/sysfuncs.h>
  7 #include <tina/math.h>
  8 #include <tina/mathfuncs.h>
  9 
 10 #define LOWACC 0.5
 11 #define TARGET 1.0
 12 #define HIGACC 2.0
 13 #define TOL 0.000001
 14 #define ITMAX 3
 15 
 16 /* A numerical routine to determine the jacobian `jf' of a function
 17  * `func' for `n' data points with respect to `m' parameters `a' using
 18  * finite differencing at a scale determined by `da'. The floating
 19  * point accuracy is determined by TARGET as the fractional change
 20  * required on the evaluated function. This should give accurate
 21  * derivative evaluations for chi-squared evaluations at a minimum. The
 22  * chi-squared function is assumed to be approximately quadratic to
 23  * allow tuning of the step parameters if the fractional function
 24  * change does not fall between the specified parameters. NAT 26/2/91 */
 25 void    jacob(double *a, double *da, Matrix * jf, double (*func) ( /* ??? */ ))
 26 {
 27     int     i, j, iter;
 28     double  acc, chisq, chisq1, aj;
 29     double *f0 = NULL, *f1, *f2;
 30 
 31     if (jf == NULL || a == NULL || da == NULL || func == NULL)
 32         return;
 33 
 34     f1 = (double *) ralloc((unsigned) jf->m * sizeof(double));
 35     f2 = (double *) ralloc((unsigned) jf->m * sizeof(double));
 36 
 37     chisq = (*func) (f0, jf->m, a, jf->n);
 38 
 39     for (j = 0; j < jf->n; j++)
 40     {
 41         iter = 0;
 42         aj = a[j];
 43         if (fabs(da[j] / a[j]) < TOL)
 44             da[j] = a[j] * TOL;
 45         acc = TARGET;
 46         while (iter == 0 || (iter < ITMAX && (acc < LOWACC || acc > HIGACC)))
 47         {
 48             if (acc < TOL)
 49                 da[j] /= TOL;
 50             else
 51                 da[j] *= sqrt(TARGET / acc);
 52             a[j] = aj + da[j];
 53             chisq1 = (*func) (f1, jf->m, a, jf->n);
 54             acc = (fabs(chisq1 - chisq) / chisq);
 55             iter++;
 56         }
 57         a[j] = aj - da[j];
 58         (void) (*func) (f2, jf->m, a, jf->n); 
 59         a[j] = aj;
 60         for (i = 0; i < jf->m; i++)
 61             jf->el.double_v[i][j] = (f1[i] - f2[i]) / (2.0 * da[j]);
 62     }
 63     /* reset the a parameters to their starting values */
 64     chisq = (*func) (f0, jf->m, a, jf->n);
 65     rfree((void *) f1);
 66     rfree((void *) f2);
 67 }
 68 
 69 Covar  *invcov(int m, double *a, double (*func) ( /* ??? */ ), int n)
 70 
 71 
 72 
 73 
 74 /* function to evaluate the inverse covariance matrix from the jacobian
 75  * NAT 27/2/91 */
 76 {
 77     int     i, j, k;
 78     double *f, *da;
 79     Matrix *jf;
 80     Matrix *matrix_alloc();
 81     Covar  *incov;
 82     Covar  *covar_alloc(int n);
 83     Vector *vector_alloc();
 84 
 85     if (a == NULL || func == NULL)
 86         return (NULL);
 87 
 88     jf = matrix_alloc(n, m, matrix_full, double_v);
 89     incov = covar_alloc(m);
 90     incov->mat = matrix_alloc(m, m, matrix_full, double_v);
 91     incov->vec = vector_alloc(m, double_v);
 92     f = (double *) ralloc((unsigned) n * sizeof(double));
 93     da = (double *) ralloc((unsigned) m * sizeof(double));
 94 
 95     for (i = 0; i < m; i++)
 96         da[i] = 0.001;
 97     jacob(a, da, jf, func);
 98 
 99     for (k = 0; k < m; k++)
100     {
101         VECTOR_DOUBLE(incov->vec, k) = a[k];
102         for (j = 0; j < m; j++)
103         {
104             incov->mat->el.double_v[j][k] = 0.0;
105             for (i = 0; i < n; i++)
106             {
107                 incov->mat->el.double_v[j][k] += jf->el.double_v[i][j] * jf->el.double_v[i][k];
108             }
109         }
110     }
111     rfree((void *) f);
112     rfree((void *) da);
113     matrix_free(jf);
114     return (incov);
115 }
116 
117 Covar  *covar(int m, double *a, double (*func) ( /* ??? */ ), int n, double condition)
118 {
119     Matrix *matrix_alloc();
120     Matrix *jf;
121     Matrix *v = NULL;
122     Matrix *u = NULL;
123     double *w = NULL, *da, *f;
124     double  sum;
125     Covar  *cov;
126     Covar  *covar_alloc(int n);
127     int     i, j, k;
128     void    matrix_svd();
129     Vector *vector_alloc();
130 
131     if (a == NULL || func == NULL)
132         return (NULL);
133 
134     da = (double *) ralloc((unsigned) m * sizeof(double));
135     f = (double *) ralloc((unsigned) n * sizeof(double));
136     jf = matrix_alloc(n, m, matrix_full, double_v);
137     cov = covar_alloc(m);
138     cov->mat = matrix_alloc(m, m, matrix_full, double_v);
139     cov->vec = vector_alloc(m, double_v);
140 
141     for (i = 0; i < m; i++)
142         da[i] = 0.001;
143     jacob(a, da, jf, func);
144     matrix_svd(jf, &u, &v, &w, condition);
145 
146     /* Sum contributions to the covariance matrix */
147     for (i = 0; i < m; i++)
148     {
149         VECTOR_DOUBLE(cov->vec, i) = a[i];
150         for (j = 0; j <= i; j++)
151         {
152             for (sum = 0.0, k = 0; k < m; k++)
153                 if (w[k])
154                     sum += v->el.double_v[i][k] * v->el.double_v[j][k] / (w[k] * w[k]);
155             cov->mat->el.double_v[j][i] = cov->mat->el.double_v[i][j] = sum;
156         }
157     }
158 
159     rfree((void *) w);
160     rfree((void *) da);
161     rfree((void *) f);
162     matrix_free(jf);
163     matrix_free(v);
164     matrix_free(u);
165 
166     return (cov);
167 }
168 
169 Covar  *covar_update(Covar * inv_acov, Covar * inv_bcov, double condition)
170 {
171     int     i;
172     Matrix *matrix_prod();
173     Matrix *matrix_invsvd();
174     Matrix *matrix_alloc();
175     Matrix *matrix_add();
176     Matrix *delta;
177     Matrix *new_delta;
178     Matrix *comb;
179     Covar  *tot_invc;
180     Covar  *covar_alloc(int n);
181     Matrix *tot_cov;
182 
183     if (inv_acov == NULL)
184         return (NULL);
185     if (inv_bcov == NULL)
186         return (inv_acov);
187 
188     delta = matrix_alloc(1, inv_acov->n, matrix_full, double_v);
189     for (i = 0; i < inv_acov->n; i++)
190         delta->el.double_v[0][i] = VECTOR_DOUBLE(inv_acov->vec, i)
191             - VECTOR_DOUBLE(inv_bcov->vec, i);
192 
193     tot_invc = covar_alloc(inv_acov->n);
194     tot_invc->n = inv_acov->n;
195     tot_invc->label = inv_acov->label;
196     tot_invc->vec = vector_alloc(inv_acov->n, double_v);
197     tot_invc->mat = matrix_add(inv_acov->mat, inv_bcov->mat);
198 
199     tot_cov = matrix_invsvd(tot_invc->mat, condition);
200     comb = matrix_prod(delta, inv_acov->mat);
201     new_delta = matrix_prod(comb, tot_cov);
202 
203     for (i = 0; i < inv_acov->n; i++)
204     {
205         VECTOR_DOUBLE(tot_invc->vec, i) = VECTOR_DOUBLE(inv_bcov->vec, i)
206             + new_delta->el.double_v[0][i];
207     }
208 
209     matrix_free(delta);
210     matrix_free(comb);
211     matrix_free(tot_cov);
212     matrix_free(new_delta);
213 
214     return (tot_invc);
215 }
216 
217 Covar  *covar_add(Covar * inv_acov, Covar * inv_bcov)
218 {
219     Covar  *tot_invc;
220     Covar  *covar_alloc(int n);
221     Vector *vector_alloc();
222     Matrix *matrix_add();
223     Matrix *matrix_copy();
224     int     i;
225 
226     if (inv_acov == NULL && inv_bcov == NULL)
227         return (NULL);
228     if (inv_acov == NULL && inv_bcov != NULL)
229     {
230         inv_acov = inv_bcov;
231         inv_bcov = NULL;
232     }
233     if (inv_acov != NULL && inv_bcov != NULL && (inv_acov->label != inv_bcov->label))
234         return (NULL);
235 
236     tot_invc = covar_alloc(inv_acov->n);
237     tot_invc->n = inv_acov->n;
238     tot_invc->label = inv_acov->label;
239     tot_invc->vec = vector_alloc(inv_acov->n, double_v);
240     if (inv_bcov != NULL)
241         tot_invc->mat = matrix_add(inv_acov->mat, inv_bcov->mat);
242     else
243         tot_invc->mat = matrix_copy(inv_acov->mat);
244 
245     for (i = 0; i < inv_acov->n; i++)
246     {
247         VECTOR_DOUBLE(tot_invc->vec, i) = VECTOR_DOUBLE(inv_acov->vec, i);
248     }
249 
250     return (tot_invc);
251 }
252 
253 Covar  *covar_alloc(int n)
254 {
255     Covar  *cov = ts_ralloc(Covar);
256 
257     if (cov == NULL)
258         return (NULL);
259     cov->n = n;
260     cov->mat = NULL;
261     cov->vec = NULL;
262     return (cov);
263 }
264 
265 void    covar_free(Covar * cov)
266 {
267     if (cov == NULL)
268         return;
269     if (cov->mat)
270         matrix_free(cov->mat);
271     if (cov->vec)
272         vector_free(cov->vec);
273     rfree((void *) cov);
274 }
275 

~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

This page was automatically generated by the LXR engine.
Visit the LXR main site for more information.