Line data Source code
1 : /**
2 : Bojan Nikolic <bojan@bnikolic.co.uk>
3 : Initial version 2005
4 :
5 : This file is part of BNMin1 and is licensed under GNU General
6 : Public License version 2.
7 :
8 : \file lmminutils.cxx
9 : Renamed to lmminutils.cc 2023
10 :
11 : */
12 :
13 : #include "lmminutils.h"
14 :
15 : #include <memory>
16 : #include <cmath>
17 :
18 : #include <gsl/gsl_permutation.h>
19 : #include <gsl/gsl_vector.h>
20 : #include <gsl/gsl_matrix.h>
21 : #include <gsl/gsl_cblas.h>
22 : #include <gsl/gsl_blas.h>
23 :
24 : #include "lmmin.h"
25 :
26 : namespace Minim {
27 :
28 :
29 :
30 0 : gsl_permutation * GetPerm( LMMin & l )
31 : {
32 :
33 0 : gsl_permutation * res( gsl_permutation_alloc (l.ipvt.size()));
34 :
35 : // Permutations are taken from F77 hence probably 1-relative
36 : // rather than 0-relative which are taken by gsl.
37 0 : for (unsigned i (0) ; i < l.ipvt.size(); ++i )
38 0 : res->data[i] = l.ipvt[i]-1;
39 :
40 0 : return res;
41 : }
42 :
43 0 : gsl_matrix * GetR ( LMMin &l )
44 : {
45 :
46 0 : unsigned n = l.NParam();
47 0 : unsigned m = l.NRes();
48 :
49 0 : gsl_matrix * res( gsl_matrix_alloc (m,n ));
50 :
51 :
52 0 : for (unsigned i (0) ; i < m ; ++i )
53 0 : for (unsigned j(0) ; j < n ; ++j )
54 0 : gsl_matrix_set( res , i, j , l.fjac[ i + j * m ] );
55 0 : return res;
56 :
57 : }
58 :
59 0 : std::vector<double> * CVMatrix ( LMMin &l , double epsrel )
60 : {
61 :
62 :
63 :
64 0 : unsigned n = l.NParam();
65 :
66 0 : std::unique_ptr< std::vector<double> > res ( new std::vector<double>( n * n ));
67 :
68 0 : gsl_matrix * r = GetR(l);
69 :
70 0 : gsl_matrix * covar (gsl_matrix_alloc (n,n ) );
71 :
72 0 : gsl_permutation *perm = GetPerm(l);
73 :
74 :
75 :
76 :
77 :
78 :
79 0 : double tolr = epsrel * fabs(gsl_matrix_get(r, 0, 0));
80 :
81 : unsigned k,j,i,kmax;
82 0 : kmax=0;
83 :
84 0 : for (k = 0 ; k < n ; k++)
85 : {
86 0 : double rkk = gsl_matrix_get(r, k, k);
87 :
88 0 : if (fabs(rkk) <= tolr)
89 : {
90 0 : break;
91 : }
92 :
93 0 : gsl_matrix_set(r, k, k, 1.0/rkk);
94 :
95 0 : for (j = 0; j < k ; j++)
96 : {
97 0 : double t = gsl_matrix_get(r, j, k) / rkk;
98 0 : gsl_matrix_set (r, j, k, 0.0);
99 :
100 0 : for (i = 0; i <= j; i++)
101 : {
102 0 : double rik = gsl_matrix_get (r, i, k);
103 0 : double rij = gsl_matrix_get (r, i, j);
104 :
105 0 : gsl_matrix_set (r, i, k, rik - t * rij);
106 : }
107 : }
108 0 : kmax = k;
109 : }
110 :
111 : /* Form the full upper triangle of the inverse of R^T R in the full
112 : upper triangle of R */
113 :
114 0 : for (k = 0; k <= kmax ; k++)
115 : {
116 0 : for (j = 0; j < k; j++)
117 : {
118 0 : double rjk = gsl_matrix_get (r, j, k);
119 :
120 0 : for (i = 0; i <= j ; i++)
121 : {
122 0 : double rij = gsl_matrix_get (r, i, j);
123 0 : double rik = gsl_matrix_get (r, i, k);
124 :
125 0 : gsl_matrix_set (r, i, j, rij + rjk * rik);
126 : }
127 : }
128 :
129 : {
130 0 : double t = gsl_matrix_get (r, k, k);
131 :
132 0 : for (i = 0; i <= k; i++)
133 : {
134 0 : double rik = gsl_matrix_get (r, i, k);
135 :
136 0 : gsl_matrix_set (r, i, k, t * rik);
137 : };
138 : }
139 : }
140 :
141 : /* Form the full lower triangle of the covariance matrix in the
142 : strict lower triangle of R and in w */
143 :
144 0 : for (j = 0 ; j < n ; j++)
145 : {
146 0 : size_t pj = gsl_permutation_get (perm, j);
147 :
148 0 : for (i = 0; i <= j; i++)
149 : {
150 0 : size_t pi = gsl_permutation_get (perm, i);
151 :
152 : double rij;
153 :
154 0 : if (j > kmax)
155 : {
156 0 : gsl_matrix_set (r, i, j, 0.0);
157 0 : rij = 0.0 ;
158 : }
159 : else
160 : {
161 0 : rij = gsl_matrix_get (r, i, j);
162 : }
163 :
164 0 : if (pi > pj)
165 : {
166 0 : gsl_matrix_set (r, pi, pj, rij);
167 : }
168 0 : else if (pi < pj)
169 : {
170 0 : gsl_matrix_set (r, pj, pi, rij);
171 : }
172 :
173 : }
174 :
175 : {
176 0 : double rjj = gsl_matrix_get (r, j, j);
177 0 : gsl_matrix_set (covar, pj, pj, rjj);
178 : }
179 : }
180 :
181 :
182 : /* symmetrize the covariance matrix */
183 :
184 0 : for (j = 0 ; j < n ; j++)
185 : {
186 0 : for (i = 0; i < j ; i++)
187 : {
188 0 : double rji = gsl_matrix_get (r, j, i);
189 :
190 0 : gsl_matrix_set (covar, j, i, rji);
191 0 : gsl_matrix_set (covar, i, j, rji);
192 : }
193 : }
194 :
195 :
196 0 : for (i = 0 ; i < n ; i++)
197 : {
198 0 : for (j = 0; j < n ; j++)
199 : {
200 0 : (*res)[j*n+i]= gsl_matrix_get (covar, i, j);
201 : }
202 : }
203 :
204 0 : gsl_matrix_free(r);
205 0 : gsl_matrix_free(covar);
206 0 : gsl_permutation_free(perm);
207 0 : return res.release();
208 0 : }
209 :
210 :
211 :
212 : }
213 :
214 :
|
