htswanalysis/MACS/lib/gsl/gsl-1.11/multifit/covar.c

   1 /* multifit/covar.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004, 2007 Brian Gough
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 #include <config.h>
  21 #include <gsl/gsl_math.h>
  22 #include <gsl/gsl_errno.h>
  23 #include <gsl/gsl_permutation.h>
  24 #include <gsl/gsl_linalg.h>
  25 #include <gsl/gsl_multifit_nlin.h>
  26
  27 /* Compute the covariance matrix
  28
  29    cov = inv (J^T J)
  30
  31    by QRP^T decomposition of J
  32 */
  33
  34 int
  35 gsl_multifit_covar (const gsl_matrix * J, double epsrel, gsl_matrix * covar)
  36 {
  37   double tolr;
  38
  39   size_t i, j, k;
  40   size_t kmax = 0;
  41
  42   gsl_matrix * r;
  43   gsl_vector * tau;
  44   gsl_vector * norm;
  45   gsl_permutation * perm;
  46
  47   size_t m = J->size1, n = J->size2 ;
  48
  49   if (m < n)
  50     {
  51       GSL_ERROR ("Jacobian be rectangular M x N with M >= N", GSL_EBADLEN);
  52     }
  53
  54   if (covar->size1 != covar->size2 || covar->size1 != n)
  55     {
  56       GSL_ERROR ("covariance matrix must be square and match second dimension of jacobian", GSL_EBADLEN);
  57     }
  58
  59   r = gsl_matrix_alloc (m, n);
  60   tau = gsl_vector_alloc (n);
  61   perm = gsl_permutation_alloc (n) ;
  62   norm = gsl_vector_alloc (n) ;
  63
  64   {
  65     int signum = 0;
  66     gsl_matrix_memcpy (r, J);
  67     gsl_linalg_QRPT_decomp (r, tau, perm, &signum, norm);
  68   }
  69
  70
  71   /* Form the inverse of R in the full upper triangle of R */
  72
  73   tolr = epsrel * fabs(gsl_matrix_get(r, 0, 0));
  74
  75   for (k = 0 ; k < n ; k++)
  76     {
  77       double rkk = gsl_matrix_get(r, k, k);
  78
  79       if (fabs(rkk) <= tolr)
  80         {
  81           break;
  82         }
  83
  84       gsl_matrix_set(r, k, k, 1.0/rkk);
  85
  86       for (j = 0; j < k ; j++)
  87         {
  88           double t = gsl_matrix_get(r, j, k) / rkk;
  89           gsl_matrix_set (r, j, k, 0.0);
  90
  91           for (i = 0; i <= j; i++)
  92             {
  93               double rik = gsl_matrix_get (r, i, k);
  94               double rij = gsl_matrix_get (r, i, j);
  95
  96               gsl_matrix_set (r, i, k, rik - t * rij);
  97             }
  98         }
  99       kmax = k;
 100     }
 101
 102   /* Form the full upper triangle of the inverse of R^T R in the full
 103      upper triangle of R */
 104
 105   for (k = 0; k <= kmax ; k++)
 106     {
 107       for (j = 0; j < k; j++)
 108         {
 109           double rjk = gsl_matrix_get (r, j, k);
 110
 111           for (i = 0; i <= j ; i++)
 112             {
 113               double rij = gsl_matrix_get (r, i, j);
 114               double rik = gsl_matrix_get (r, i, k);
 115
 116               gsl_matrix_set (r, i, j, rij + rjk * rik);
 117             }
 118         }
 119
 120       {
 121         double t = gsl_matrix_get (r, k, k);
 122
 123         for (i = 0; i <= k; i++)
 124           {
 125             double rik = gsl_matrix_get (r, i, k);
 126
 127             gsl_matrix_set (r, i, k, t * rik);
 128           };
 129       }
 130     }
 131
 132   /* Form the full lower triangle of the covariance matrix in the
 133      strict lower triangle of R and in w */
 134
 135   for (j = 0 ; j < n ; j++)
 136     {
 137       size_t pj = gsl_permutation_get (perm, j);
 138
 139       for (i = 0; i <= j; i++)
 140         {
 141           size_t pi = gsl_permutation_get (perm, i);
 142
 143           double rij;
 144
 145           if (j > kmax)
 146             {
 147               gsl_matrix_set (r, i, j, 0.0);
 148               rij = 0.0 ;
 149             }
 150           else
 151             {
 152               rij = gsl_matrix_get (r, i, j);
 153             }
 154
 155           if (pi > pj)
 156             {
 157               gsl_matrix_set (r, pi, pj, rij);
 158             }
 159           else if (pi < pj)
 160             {
 161               gsl_matrix_set (r, pj, pi, rij);
 162             }
 163
 164         }
 165
 166       {
 167         double rjj = gsl_matrix_get (r, j, j);
 168         gsl_matrix_set (covar, pj, pj, rjj);
 169       }
 170     }
 171
 172
 173   /* symmetrize the covariance matrix */
 174
 175   for (j = 0 ; j < n ; j++)
 176     {
 177       for (i = 0; i < j ; i++)
 178         {
 179           double rji = gsl_matrix_get (r, j, i);
 180
 181           gsl_matrix_set (covar, j, i, rji);
 182           gsl_matrix_set (covar, i, j, rji);
 183         }
 184     }
 185
 186   gsl_matrix_free (r);
 187   gsl_permutation_free (perm);
 188   gsl_vector_free (tau);
 189   gsl_vector_free (norm);
 190
 191   return GSL_SUCCESS;
 192 }