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

Linux Cross Reference
Tina5/tina-libs/tina/math/mathStat_covar.c

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

  1 /**********
  2  *
  3  * Copyright (c) 2003, Division of Imaging Science and Biomedical Engineering,
  4  * University of Manchester, UK.  All rights reserved.
  5  * 
  6  * Redistribution and use in source and binary forms, with or without modification, 
  7  * are permitted provided that the following conditions are met:
  8  * 
  9  *   . Redistributions of source code must retain the above copyright notice, 
 10  *     this list of conditions and the following disclaimer.
 11  *    
 12  *   . Redistributions in binary form must reproduce the above copyright notice,
 13  *     this list of conditions and the following disclaimer in the documentation 
 14  *     and/or other materials provided with the distribution.
 15  * 
 16  *   . Neither the name of the University of Manchester nor the names of its
 17  *     contributors may be used to endorse or promote products derived from this 
 18  *     software without specific prior written permission.
 19  * 
 20  * 
 21  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
 22  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 
 23  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 
 24  * ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE 
 25  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
 26  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 
 27  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 28  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 
 29  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 
 30  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 
 31  * POSSIBILITY OF SUCH DAMAGE.
 32  *
 33  **********
 34  *
 35  * Program :    TINA
 36  * File    :  $Source: /home/tina/cvs/tina-libs/tina/math/mathStat_covar.c,v $
 37  * Date    :  $Date: 2009/03/19 18:36:26 $
 38  * Version :  $Revision: 1.5 $
 39  * CVS Id  :  $Id: mathStat_covar.c,v 1.5 2009/03/19 18:36:26 paul Exp $
 40  *
 41  * Author  : Legacy TINA
 42  *
 43  * Notes :Jacobian & evaluate the inverse covariance matrix from the
 44  * jacobian
 45  *
 46  *********
 47 */
 48 /** 
 49  *  @file
 50  *  @brief  Functions for evaluating the Jacobian and inverse covariance matrices.      
 51  *
 52  *  A numerical routine to determine the jacobian `jf' of a function `func'
 53  *  for `n' data points with respect to `m' parameters `a' using finite
 54  *  differencing at a scale determined by `da'.
 55  * 
 56 */
 57 
 58 #include "mathStat_covar.h"
 59 
 60 
 61 #if HAVE_CONFIG_H
 62 #include <config.h>
 63 #endif
 64 
 65 #include <math.h>
 66 #include <tina/math/math_VecDef.h>
 67 #include <tina/math/math_VecPro.h>
 68 #include <tina/math/math_StatDef.h>
 69 #include <tina/sys/sysPro.h>
 70 #include <tina/math/math_MatrPro.h>
 71 
 72 
 73 /*
 74  * A numerical routine to determine the jacobian `jf' of a function `func'
 75  * for `n' data points with respect to `m' parameters `a' using finite
 76  * differencing at a scale determined by `da'. The floating point accuracy is
 77  * determined by TARGET as the fractional change required on the evaluated
 78  * function. This should give accurate derivative evaluations for chi-squared
 79  * evaluations at a minimum. The chi-squared function is assumed to be
 80  * approximately quadratic to allow tuning of the step parameters if the
 81  * fractional function change does not fall between the specified parameters.
 82  * NAT 26/2/91
 83  */
 84 
 85 
 86 Covar          *covar_alloc(int n)
 87 {
 88         Covar          *cov = ts_ralloc(Covar);
 89 
 90         if (cov == NULL)
 91                 return (NULL);
 92         cov->n = n;
 93         cov->mat = NULL;
 94         cov->vec = NULL;
 95         return (cov);
 96 }
 97 
 98 
 99 void            covar_free(Covar * cov)
100 {
101         if (cov == NULL)
102                 return;
103         if (cov->mat)
104                 matrix_free(cov->mat);
105         if (cov->vec)
106                 vector_free(cov->vec);
107         rfree((void *) cov);
108 }
109 
110 
111 void            jacob(double *a, double *da, Matrix * jf, double (*func) ( /* ??? */ ))
112 {
113         int             i, j, iter;
114         double          acc, chisq, chisq1, aj;
115         double         *f0 = NULL, *f1, *f2;
116 
117         if (jf == NULL || a == NULL || da == NULL || func == NULL)
118                 return;
119 
120         f1 = (double *) ralloc((unsigned) jf->m * sizeof(double));
121         f2 = (double *) ralloc((unsigned) jf->m * sizeof(double));
122 
123         chisq = (*func) (f0, jf->m, a, jf->n);
124 
125         for (j = 0; j < jf->n; j++)
126         {
127                 iter = 0;
128                 aj = a[j];
129                 if (fabs(da[j] / a[j]) < TOL)
130                         da[j] = a[j] * TOL;
131                 acc = TARGET;
132                 while (iter == 0 || (iter < ITMAX && (acc < LOWACC || acc > HIGACC)))
133                 {
134                         if (acc < TOL)
135                                 da[j] /= TOL;
136                         else
137                                 da[j] *= sqrt(TARGET / acc);
138                         a[j] = aj + da[j];
139                         chisq1 = (*func) (f1, jf->m, a, jf->n);
140                         acc = (fabs(chisq1 - chisq) / chisq);
141                         iter++;
142                 }
143                 a[j] = aj - da[j];
144                 (void) (*func) (f2, jf->m, a, jf->n);
145                 a[j] = aj;
146                 for (i = 0; i < jf->m; i++)
147                         jf->el.double_v[i][j] = (f1[i] - f2[i]) / (2.0 * da[j]);
148         }
149         /* reset the a parameters to their starting values */
150         chisq = (*func) (f0, jf->m, a, jf->n);
151         rfree((void *) f1);
152         rfree((void *) f2);
153 }
154 
155 Covar          *invcov(int m, double *a, double (*func) ( /* ??? */ ), int n)
156 {
157         /* Function to evaluate the inverse covariance matrix from the jacobian NAT 27/2/91 */
158 
159         int             i, j, k;
160         double         *f, *da;
161         Matrix         *jf;
162         Covar          *incov;
163 
164         if (a == NULL || func == NULL)
165                 return (NULL);
166 
167         jf = matrix_alloc(n, m, matrix_full, double_v);
168         incov = covar_alloc(m);
169         incov->mat = matrix_alloc(m, m, matrix_full, double_v);
170         incov->vec = vector_alloc(m, double_v);
171         f = (double *) ralloc((unsigned) n * sizeof(double));
172         da = (double *) ralloc((unsigned) m * sizeof(double));
173 
174         for (i = 0; i < m; i++)
175                 da[i] = 0.001;
176         jacob(a, da, jf, func);
177 
178         for (k = 0; k < m; k++)
179         {
180                 VECTOR_DOUBLE(incov->vec, k) = a[k];
181                 for (j = 0; j < m; j++)
182                 {
183                         incov->mat->el.double_v[j][k] = 0.0;
184                         for (i = 0; i < n; i++)
185                         {
186                                 incov->mat->el.double_v[j][k] += jf->el.double_v[i][j] * jf->el.double_v[i][k];
187                         }
188                 }
189         }
190         rfree((void *) f);
191         rfree((void *) da);
192         matrix_free(jf);
193         return (incov);
194 }
195 
196 Covar          *covar(int m, double *a, double (*func) ( /* ??? */ ), int n, double condition)
197 {
198         Matrix         *jf;
199         Matrix         *v = NULL;
200         Matrix         *u = NULL;
201         double         *w = NULL, *da, *f;
202         double          sum;
203         Covar          *cov;
204         int             i, j, k;
205 
206         if (a == NULL || func == NULL)
207                 return (NULL);
208 
209         da = (double *) ralloc((unsigned) m * sizeof(double));
210         f = (double *) ralloc((unsigned) n * sizeof(double));
211         jf = matrix_alloc(n, m, matrix_full, double_v);
212         cov = covar_alloc(m);
213         cov->mat = matrix_alloc(m, m, matrix_full, double_v);
214         cov->vec = vector_alloc(m, double_v);
215 
216         for (i = 0; i < m; i++)
217                 da[i] = 0.001;
218         jacob(a, da, jf, func);
219         matrix_svd(jf, &u, &v, &w, condition);
220 
221         /* Sum contributions to the covariance matrix */
222         for (i = 0; i < m; i++)
223         {
224                 VECTOR_DOUBLE(cov->vec, i) = a[i];
225                 for (j = 0; j <= i; j++)
226                 {
227                         for (sum = 0.0, k = 0; k < m; k++)
228                                 if (w[k])
229                                         sum += v->el.double_v[i][k] * v->el.double_v[j][k] / (w[k] * w[k]);
230                         cov->mat->el.double_v[j][i] = cov->mat->el.double_v[i][j] = sum;
231                 }
232         }
233 
234         rfree((void *) w);
235         rfree((void *) da);
236         rfree((void *) f);
237         matrix_free(jf);
238         matrix_free(v);
239         matrix_free(u);
240 
241         return (cov);
242 }
243 
244 Covar          *covar_update(Covar * inv_acov, Covar * inv_bcov, double condition)
245 {
246         int             i;
247         Matrix         *delta;
248         Matrix         *new_delta;
249         Matrix         *comb;
250         Covar          *tot_invc;
251         Matrix         *tot_cov;
252 
253         if (inv_acov == NULL)
254                 return (NULL);
255         if (inv_bcov == NULL)
256                 return (inv_acov);
257 
258         delta = matrix_alloc(1, inv_acov->n, matrix_full, double_v);
259         for (i = 0; i < inv_acov->n; i++)
260                 delta->el.double_v[0][i] = VECTOR_DOUBLE(inv_acov->vec, i)
261                         - VECTOR_DOUBLE(inv_bcov->vec, i);
262 
263         tot_invc = covar_alloc(inv_acov->n);
264         tot_invc->n = inv_acov->n;
265         tot_invc->label = inv_acov->label;
266         tot_invc->vec = vector_alloc(inv_acov->n, double_v);
267         tot_invc->mat = matrix_add(inv_acov->mat, inv_bcov->mat);
268 
269         tot_cov = matrix_invsvd(tot_invc->mat, condition);
270         comb = matrix_prod(delta, inv_acov->mat);
271         new_delta = matrix_prod(comb, tot_cov);
272 
273         for (i = 0; i < inv_acov->n; i++)
274         {
275                 VECTOR_DOUBLE(tot_invc->vec, i) = VECTOR_DOUBLE(inv_bcov->vec, i)
276                         + new_delta->el.double_v[0][i];
277         }
278 
279         matrix_free(delta);
280         matrix_free(comb);
281         matrix_free(tot_cov);
282         matrix_free(new_delta);
283 
284         return (tot_invc);
285 }
286 
287 Covar          *covar_add(Covar * inv_acov, Covar * inv_bcov)
288 {
289         Covar          *tot_invc;
290         int             i;
291 
292         if (inv_acov == NULL && inv_bcov == NULL)
293                 return (NULL);
294         if (inv_acov == NULL && inv_bcov != NULL)
295         {
296                 inv_acov = inv_bcov;
297                 inv_bcov = NULL;
298         }
299         if (inv_acov != NULL && inv_bcov != NULL && (inv_acov->label != inv_bcov->label))
300                 return (NULL);
301 
302         tot_invc = covar_alloc(inv_acov->n);
303         tot_invc->n = inv_acov->n;
304         tot_invc->label = inv_acov->label;
305         tot_invc->vec = vector_alloc(inv_acov->n, double_v);
306         if (inv_bcov != NULL)
307                 tot_invc->mat = matrix_add(inv_acov->mat, inv_bcov->mat);
308         else
309                 tot_invc->mat = matrix_copy(inv_acov->mat);
310 
311         for (i = 0; i < inv_acov->n; i++)
312         {
313                 VECTOR_DOUBLE(tot_invc->vec, i) = VECTOR_DOUBLE(inv_acov->vec, i);
314         }
315 
316         return (tot_invc);
317 }
318 

~ [ 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.