1 /* Linear algebra operations in double-precision matrices.
2 *
3 * Implements ESL_DMATRIX (double-precision matrix) and
4 * ESL_PERMUTATION (permutation matrix) objects.
5 *
6 * Table of contents:
7 * 1. The ESL_DMATRIX object
8 * 2. Debugging/validation code for ESL_DMATRIX
9 * 3. Visualization tools
10 * 4. The ESL_PERMUTATION object
11 * 5. Debugging/validation code for ESL_PERMUTATION
12 * 6. The rest of the dmatrix API
13 * 7. Optional: Interoperability with GSL
14 * 8. Optional: Interfaces to LAPACK
15 * 9. Unit tests
16 * 10. Test driver
17 * 11. Examples
18 *
19 * To do:
20 * - eventually probably want additional matrix types
21 * - unit tests poor
22 */
23 #include "esl_config.h"
24
25 #include <stdio.h>
26 #include <stdlib.h>
27 #include <string.h>
28 #include <math.h>
29
30 #include "easel.h"
31 #include "esl_vectorops.h"
32 #include "esl_dmatrix.h"
33
34
35 /*****************************************************************
36 * 1. The ESL_DMATRIX object.
37 *****************************************************************/
38
39 /* Function: esl_dmatrix_Create()
40 *
41 * Purpose: Creates a general <n> x <m> matrix (<n> rows, <m>
42 * columns).
43 *
44 * Args: <n> - number of rows; $>= 1$
45 * <m> - number of columns; $>= 1$
46 *
47 * Returns: a pointer to a new <ESL_DMATRIX> object. Caller frees
48 * with <esl_dmatrix_Destroy()>.
49 *
50 * Throws: <NULL> if an allocation failed.
51 */
52 ESL_DMATRIX *
esl_dmatrix_Create(int n,int m)53 esl_dmatrix_Create(int n, int m)
54 {
55 ESL_DMATRIX *A = NULL;
56 int r;
57 int status;
58
59 ESL_ALLOC(A, sizeof(ESL_DMATRIX));
60 A->mx = NULL;
61 A->n = n;
62 A->m = m;
63
64 ESL_ALLOC(A->mx, sizeof(double *) * n);
65 A->mx[0] = NULL;
66
67 ESL_ALLOC(A->mx[0], sizeof(double) * n * m);
68 for (r = 1; r < n; r++)
69 A->mx[r] = A->mx[0] + r*m;
70
71 A->type = eslGENERAL;
72 A->ncells = n * m;
73 return A;
74
75 ERROR:
76 esl_dmatrix_Destroy(A);
77 return NULL;
78 }
79
80
81 /* Function: esl_dmatrix_CreateUpper()
82 * Incept: SRE, Wed Feb 28 08:45:45 2007 [Janelia]
83 *
84 * Purpose: Creates a packed upper triangular matrix of <n> rows and
85 * <n> columns. Caller may only access cells $i \leq j$.
86 * Cells $i > j$ are not stored and are implicitly 0.
87 *
88 * Not all matrix operations in Easel can work on packed
89 * upper triangular matrices.
90 *
91 * Returns: a pointer to a new <ESL_DMATRIX> object of type
92 * <eslUPPER>. Caller frees with <esl_dmatrix_Destroy()>.
93 *
94 * Throws: <NULL> if allocation fails.
95 *
96 * Xref: J1/10
97 */
98 ESL_DMATRIX *
esl_dmatrix_CreateUpper(int n)99 esl_dmatrix_CreateUpper(int n)
100 {
101 int status;
102 ESL_DMATRIX *A = NULL;
103 int r; /* counter over rows */
104 int nc; /* cell counter */
105
106 /* matrix structure allocation */
107 ESL_ALLOC(A, sizeof(ESL_DMATRIX));
108 A->mx = NULL;
109 A->n = n;
110 A->m = n;
111
112 /* n row ptrs */
113 ESL_ALLOC(A->mx, sizeof(double *) * n);
114 A->mx[0] = NULL;
115
116 /* cell storage */
117 ESL_ALLOC(A->mx[0], sizeof(double) * n * (n+1) / 2);
118
119 /* row pointers set in a tricksy overlapping way, so
120 * mx[i][j] access works normally but only i<=j are valid.
121 * xref J1/10.
122 */
123 nc = n; /* nc is the number of valid cells assigned to rows so far */
124 for (r = 1; r < n; r++) {
125 A->mx[r] = A->mx[0] + nc - r; /* -r overlaps this row w/ previous row */
126 nc += n-r;
127 }
128 A->type = eslUPPER;
129 A->ncells = n * (n+1) / 2;
130 return A;
131
132 ERROR:
133 esl_dmatrix_Destroy(A);
134 return NULL;
135 }
136
137
138 /* Function: esl_dmatrix_Destroy()
139 *
140 * Purpose: Frees an <ESL_DMATRIX> object <A>.
141 */
142 int
esl_dmatrix_Destroy(ESL_DMATRIX * A)143 esl_dmatrix_Destroy(ESL_DMATRIX *A)
144 {
145 if (A != NULL && A->mx != NULL && A->mx[0] != NULL) free(A->mx[0]);
146 if (A != NULL && A->mx != NULL) free(A->mx);
147 if (A != NULL) free(A);
148 return eslOK;
149 }
150
151
152 /* Function: esl_dmatrix_Copy()
153 *
154 * Purpose: Copies <src> matrix into <dest> matrix. <dest> must
155 * be allocated already by the caller.
156 *
157 * You may copy to a matrix of a different type, so long as
158 * the copy makes sense. If <dest> matrix is a packed type
159 * and <src> is not, the values that should be zeros must
160 * be zero in <src>, else the routine throws
161 * <eslEINCOMPAT>. If the <src> matrix is a packed type and
162 * <dest> is not, the values that are implicitly zeros are
163 * set to zeros in the <dest> matrix.
164 *
165 * Returns: <eslOK> on success.
166 *
167 * Throws: <eslEINCOMPAT> if <src>, <dest> are different sizes,
168 * or if their types differ and <dest> cannot represent
169 * <src>.
170 */
171 int
esl_dmatrix_Copy(const ESL_DMATRIX * src,ESL_DMATRIX * dest)172 esl_dmatrix_Copy(const ESL_DMATRIX *src, ESL_DMATRIX *dest)
173 {
174 int i,j;
175
176 if (dest->n != src->n || dest->m != src->m)
177 ESL_EXCEPTION(eslEINCOMPAT, "matrices of different size");
178
179 if (src->type == dest->type) /* simple case. */
180 memcpy(dest->mx[0], src->mx[0], src->ncells * sizeof(double));
181
182 else if (src->type == eslGENERAL && dest->type == eslUPPER)
183 {
184 for (i = 1; i < src->n; i++)
185 for (j = 0; j < i; j++)
186 if (src->mx[i][j] != 0.)
187 ESL_EXCEPTION(eslEINCOMPAT, "general matrix isn't upper triangular, can't be copied/packed");
188 for (i = 0; i < src->n; i++)
189 for (j = i; j < src->m; j++)
190 dest->mx[i][j] = src->mx[i][j];
191 }
192
193 else if (src->type == eslUPPER && dest->type == eslGENERAL)
194 {
195 for (i = 1; i < src->n; i++)
196 for (j = 0; j < i; j++)
197 dest->mx[i][j] = 0.;
198 for (i = 0; i < src->n; i++)
199 for (j = i; j < src->m; j++)
200 dest->mx[i][j] = src->mx[i][j];
201 }
202
203 return eslOK;
204 }
205
206
207 /* Function: esl_dmatrix_Clone()
208 * Incept: SRE, Tue May 2 14:38:45 2006 [St. Louis]
209 *
210 * Purpose: Duplicates matrix <A>, making a copy in newly
211 * allocated space.
212 *
213 * Returns: a pointer to the copy. Caller frees with
214 * <esl_dmatrix_Destroy()>.
215 *
216 * Throws: <NULL> on allocation failure.
217 */
218 ESL_DMATRIX *
esl_dmatrix_Clone(const ESL_DMATRIX * A)219 esl_dmatrix_Clone(const ESL_DMATRIX *A)
220 {
221 ESL_DMATRIX *new;
222
223 switch (A->type) {
224 case eslUPPER: if ( (new = esl_dmatrix_CreateUpper(A->n)) == NULL) return NULL; break;
225 default: case eslGENERAL: if ( (new = esl_dmatrix_Create(A->n, A->m)) == NULL) return NULL; break;
226 }
227 esl_dmatrix_Copy(A, new);
228 return new;
229 }
230
231
232 /* Function: esl_dmatrix_Compare()
233 *
234 * Purpose: Compares matrix <A> to matrix <B> element by element,
235 * using <esl_DCompare()> on each cognate element pair,
236 * with relative equality defined by a fractional tolerance
237 * <tol>. If all elements are equal, return <eslOK>; if
238 * any elements differ, return <eslFAIL>.
239 *
240 * <A> and <B> may be of different types; for example,
241 * a packed upper triangular matrix A is compared to
242 * a general matrix B by assuming <A->mx[i][j] = 0.> for
243 * all $i>j$.
244 */
245 int
esl_dmatrix_Compare(const ESL_DMATRIX * A,const ESL_DMATRIX * B,double tol)246 esl_dmatrix_Compare(const ESL_DMATRIX *A, const ESL_DMATRIX *B, double tol)
247 {
248 int i,j,c;
249 double x1,x2;
250
251 if (A->n != B->n) return eslFAIL;
252 if (A->m != B->m) return eslFAIL;
253
254 if (A->type == B->type)
255 { /* simple case. */
256 for (c = 0; c < A->ncells; c++) /* can deal w/ packed or unpacked storage */
257 if (esl_DCompare(A->mx[0][c], B->mx[0][c], tol) == eslFAIL) return eslFAIL;
258 }
259 else
260 { /* comparing matrices of different types */
261 for (i = 0; i < A->n; i++)
262 for (j = 0; j < A->m; j++)
263 {
264 if (A->type == eslUPPER && i > j) x1 = 0.;
265 else x1 = A->mx[i][j];
266
267 if (B->type == eslUPPER && i > j) x2 = 0.;
268 else x2 = B->mx[i][j];
269
270 if (esl_DCompare(x1, x2, tol) == eslFAIL) return eslFAIL;
271 }
272 }
273 return eslOK;
274 }
275
276
277 /* Function: esl_dmatrix_CompareAbs()
278 *
279 * Purpose: Compares matrix <A> to matrix <B> element by element,
280 * using <esl_DCompareAbs()> on each cognate element pair,
281 * with absolute equality defined by a absolute difference tolerance
282 * <tol>. If all elements are equal, return <eslOK>; if
283 * any elements differ, return <eslFAIL>.
284 *
285 * <A> and <B> may be of different types; for example,
286 * a packed upper triangular matrix A is compared to
287 * a general matrix B by assuming <A->mx[i][j] = 0.> for
288 * all $i>j$.
289 */
290 int
esl_dmatrix_CompareAbs(const ESL_DMATRIX * A,const ESL_DMATRIX * B,double tol)291 esl_dmatrix_CompareAbs(const ESL_DMATRIX *A, const ESL_DMATRIX *B, double tol)
292 {
293 int i,j,c;
294 double x1,x2;
295
296 if (A->n != B->n) return eslFAIL;
297 if (A->m != B->m) return eslFAIL;
298
299 if (A->type == B->type)
300 { /* simple case. */
301 for (c = 0; c < A->ncells; c++) /* can deal w/ packed or unpacked storage */
302 if (esl_DCompareAbs(A->mx[0][c], B->mx[0][c], tol) == eslFAIL) return eslFAIL;
303 }
304 else
305 { /* comparing matrices of different types */
306 for (i = 0; i < A->n; i++)
307 for (j = 0; j < A->m; j++)
308 {
309 if (A->type == eslUPPER && i > j) x1 = 0.;
310 else x1 = A->mx[i][j];
311
312 if (B->type == eslUPPER && i > j) x2 = 0.;
313 else x2 = B->mx[i][j];
314
315 if (esl_DCompareAbs(x1, x2, tol) == eslFAIL) return eslFAIL;
316 }
317 }
318 return eslOK;
319 }
320
321
322 /* Function: esl_dmatrix_Set()
323 *
324 * Purpose: Set all elements $a_{ij}$ in matrix <A> to <x>,
325 * and returns <eslOK>.
326 */
327 int
esl_dmatrix_Set(ESL_DMATRIX * A,double x)328 esl_dmatrix_Set(ESL_DMATRIX *A, double x)
329 {
330 int i;
331 for (i = 0; i < A->ncells; i++) A->mx[0][i] = x;
332 return eslOK;
333 }
334
335
336 /* Function: esl_dmatrix_SetZero()
337 *
338 * Purpose: Sets all elements $a_{ij}$ in matrix <A> to 0,
339 * and returns <eslOK>.
340 */
341 int
esl_dmatrix_SetZero(ESL_DMATRIX * A)342 esl_dmatrix_SetZero(ESL_DMATRIX *A)
343 {
344 int i;
345 for (i = 0; i < A->ncells; i++) A->mx[0][i] = 0.;
346 return eslOK;
347 }
348
349
350 /* Function: esl_dmatrix_SetIdentity()
351 *
352 * Purpose: Given a square matrix <A>, sets all diagonal elements
353 * $a_{ii}$ to 1, and all off-diagonal elements $a_{ij},
354 * j \ne i$ to 0. Returns <eslOK> on success.
355 *
356 * Throws: <eslEINVAL> if the matrix isn't square.
357 */
358 int
esl_dmatrix_SetIdentity(ESL_DMATRIX * A)359 esl_dmatrix_SetIdentity(ESL_DMATRIX *A)
360 {
361 int i;
362
363 if (A->n != A->m) ESL_EXCEPTION(eslEINVAL, "matrix isn't square");
364 esl_dmatrix_SetZero(A);
365 for (i = 0; i < A->n; i++) A->mx[i][i] = 1.;
366 return eslOK;
367 }
368
369 /*****************************************************************
370 * 2. Debugging, validation code
371 *****************************************************************/
372
373 /* Function: esl_dmatrix_Dump()
374 * Incept: SRE, Mon Nov 29 19:21:20 2004 [St. Louis]
375 *
376 * Purpose: Given a matrix <A>, dump it to output stream <ofp> in human-readable
377 * format.
378 *
379 * If <rowlabel> or <collabel> are non-NULL, they specify a
380 * string of single-character labels to put on the rows and
381 * columns, respectively. (For example, these might be a
382 * sequence alphabet for a 4x4 or 20x20 rate matrix or
383 * substitution matrix.) Numbers <1..ncols> or <1..nrows> are
384 * used if <collabel> or <rowlabel> are passed as <NULL>.
385 *
386 * Args: ofp - output file pointer; stdout, for example.
387 * A - matrix to dump.
388 * rowlabel - optional: NULL, or character labels for rows
389 * collabel - optional: NULL, or character labels for cols
390 *
391 * Returns: <eslOK> on success.
392 */
393 int
esl_dmatrix_Dump(FILE * ofp,const ESL_DMATRIX * A,const char * rowlabel,const char * collabel)394 esl_dmatrix_Dump(FILE *ofp, const ESL_DMATRIX *A, const char *rowlabel, const char *collabel)
395 {
396 int a,b;
397
398 fprintf(ofp, " ");
399 if (collabel != NULL)
400 for (b = 0; b < A->m; b++) fprintf(ofp, " %c ", collabel[b]);
401 else
402 for (b = 0; b < A->m; b++) fprintf(ofp, "%8d ", b+1);
403 fprintf(ofp, "\n");
404
405 for (a = 0; a < A->n; a++) {
406 if (rowlabel != NULL) fprintf(ofp, " %c ", rowlabel[a]);
407 else fprintf(ofp, "%5d ", a+1);
408
409 for (b = 0; b < A->m; b++) {
410 switch (A->type) {
411 case eslUPPER:
412 if (a > b) fprintf(ofp, "%8s ", "");
413 else fprintf(ofp, "%8.4f ", A->mx[a][b]);
414 break;
415
416 default: case eslGENERAL:
417 fprintf(ofp, "%8.4f ", A->mx[a][b]);
418 break;
419 }
420 }
421 fprintf(ofp, "\n");
422 }
423 return eslOK;
424 }
425
426
427 /*****************************************************************
428 * 3. Visualization tools
429 *****************************************************************/
430
431 /* Function: esl_dmatrix_PlotHeatMap()
432 * Synopsis: Export a heat map visualization, in PostScript
433 *
434 * Purpose: Export a heat map visualization of the matrix in <D>
435 * to open stream <fp>, in PostScript format.
436 *
437 * All values between <min> and <max> in <D> are rescaled
438 * linearly and assigned to shades. Values below <min>
439 * are assigned to the lowest shade; values above <max>, to
440 * the highest shade.
441 *
442 * The plot is hardcoded to be a full US 8x11.5" page,
443 * with at least a 20pt margin.
444 *
445 * Several color schemes are enumerated in the code
446 * but all but one is commented out. The currently enabled
447 * scheme is a 10-class scheme consisting of the 9-class
448 * Reds from colorbrewer2.org plus a blue background class.
449 *
450 * Note: Binning rules basically follow same convention as
451 * esl_histogram. nb = xmax-xmin/w, so w = xmax-xmin/nb;
452 * picking bin is (int) ceil((x - xmin)/w) - 1. (xref
453 * esl_histogram_Score2Bin()). This makes bin b contain
454 * values bw+min < x <= (b+1)w+min. (Which means that
455 * min itself falls in bin -1, whoops - but we catch
456 * all bin<0 and bin>=nshades and put them in the extremes.
457 *
458 * Returns: <eslOK> on success.
459 *
460 * Throws: (no abnormal error conditions)
461 */
462 int
esl_dmatrix_PlotHeatMap(FILE * fp,ESL_DMATRIX * D,double min,double max)463 esl_dmatrix_PlotHeatMap(FILE *fp, ESL_DMATRIX *D, double min, double max)
464 {
465 #if 0
466 /*
467 * This color scheme roughly follows Tufte, Envisioning Information,
468 * p.91, where he shows a beautiful bathymetric chart. The CMYK
469 * values conjoin two recommendations from ColorBrewer (Cindy Brewer
470 * and Mark Harrower, colorbrewer2.org), specifically the 9-class
471 * sequential2 Blues and 9-class sequential YlOrBr.
472 */
473 int nshades = 18;
474 double cyan[] = { 1.00, 1.00, 0.90, 0.75, 0.57, 0.38, 0.24, 0.13, 0.03,
475 0.00, 0.00, 0.00, 0.00, 0.00, 0.07, 0.20, 0.40, 0.60};
476 double magenta[] = { 0.55, 0.45, 0.34, 0.22, 0.14, 0.08, 0.06, 0.03, 0.01,
477 0.00, 0.03, 0.11, 0.23, 0.40, 0.55, 0.67, 0.75, 0.80};
478 double yellow[] = { 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00,
479 0.10, 0.25, 0.40, 0.65, 0.80, 0.90, 1.00, 1.00, 1.00};
480 double black[] = { 0.30, 0.07, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00,
481 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00};
482 #endif
483 #if 0
484 /* colorbrewer2.org 5-class YlOrBr scheme: sequential, multihue, 5-class, CMYK */
485 int nshades = 5;
486 double cyan[] = { 0.00, 0.00, 0.00, 0.15, 0.40 };
487 double magenta[] = { 0.00, 0.15, 0.40, 0.60, 0.75 };
488 double yellow[] = { 0.17, 0.40, 0.80, 0.95, 1.00 };
489 double black[] = { 0, 0, 0, 0, 0 };
490 #endif
491 #if 0
492 /* colorbrewer2.org 9-class YlOrBr scheme, +zero class */
493 int nshades = 10;
494 double cyan[] = { 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.07, 0.20, 0.40, 0.60 };
495 double magenta[] = { 0.00, 0.00, 0.03, 0.11, 0.23, 0.40, 0.55, 0.67, 0.75, 0.80 };
496 double yellow[] = { 0.00, 0.10, 0.25, 0.40, 0.65, 0.80, 0.90, 1.00, 1.00, 1.00 };
497 double black[] = { 0.05, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 };
498 #endif
499 /* colorbrewer2.org 9-class Reds + zero class as dim blue */
500 int nshades = 10;
501 double cyan[] = { 0.30, 0.00, 0.00, 0.00, 0.00, 0.00, 0.05, 0.20, 0.35, 0.60 };
502 double magenta[] = { 0.03, 0.04, 0.12, 0.27, 0.43, 0.59, 0.77, 0.90, 0.95, 1.00 };
503 double yellow[] = { 0.00, 0.04, 0.12, 0.27, 0.43, 0.59, 0.72, 0.80, 0.85, 0.90 };
504 double black[] = { 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 };
505
506 int pageheight = 792;
507 int pagewidth = 612;
508 double w;
509 int i,j;
510 int bin;
511 float boxsize; /* box size in points */
512 float xcoord, ycoord; /* postscript coords in points */
513 float leftmargin;
514 float bottommargin;
515
516 /* Set some defaults that might become arguments later.
517 */
518 leftmargin = 20.;
519 bottommargin = 20.;
520
521 /* Determine some working parameters
522 */
523 w = (max-min) / (double) nshades; /* w = bin size for assigning values->colors*/
524 boxsize = ESL_MIN( (pageheight - (bottommargin * 2.)) / (float) D->n,
525 (pagewidth - (leftmargin * 2.)) / (float) D->m);
526
527 /* or start from j=i, to do diagonals */
528 for (i = 0; i < D->n; i++)
529 for (j = 0; j < D->m; j++)
530 {
531 xcoord = (float) j * boxsize + leftmargin;
532 ycoord = (float) (D->n-i+1) * boxsize + bottommargin;
533
534 if (D->mx[i][j] == -eslINFINITY) bin = 0;
535 else if (D->mx[i][j] == eslINFINITY) bin = nshades-1;
536 else {
537 bin = (int) ceil((D->mx[i][j] - min) / w) - 1;
538 if (bin < 0) bin = 0;
539 if (bin >= nshades) bin = nshades-1;
540 }
541
542 fprintf(fp, "newpath\n");
543 fprintf(fp, " %.2f %.2f moveto\n", xcoord, ycoord);
544 fprintf(fp, " 0 %.2f rlineto\n", boxsize);
545 fprintf(fp, " %.2f 0 rlineto\n", boxsize);
546 fprintf(fp, " 0 -%.2f rlineto\n", boxsize);
547 fprintf(fp, " closepath\n");
548 fprintf(fp, " %.2f %.2f %.2f %.2f setcmykcolor\n",
549 cyan[bin], magenta[bin], yellow[bin], black[bin]);
550 fprintf(fp, " fill\n");
551 }
552 fprintf(fp, "showpage\n");
553 return eslOK;
554 }
555
556
557
558 /*****************************************************************
559 * 4. The ESL_PERMUTATION object.
560 *****************************************************************/
561
562 /* Function: esl_permutation_Create()
563 *
564 * Purpose: Creates a new permutation "matrix" of size <n> for
565 * permuting <n> x <n> square matrices; returns a
566 * pointer to it.
567 *
568 * A permutation matrix consists of 1's and 0's such that
569 * any given row or column contains only one 1. We store it
570 * more efficiently as a vector; each value $p_i$
571 * represents the column $j$ that has the 1. Thus, on
572 * initialization, $p_i = i$ for all $i = 0..n-1$.
573 *
574 * Returns: a pointer to a new <ESL_PERMUTATION> object. Free with
575 * <esl_permutation_Destroy()>.
576 *
577 * Throws: <NULL> if allocation fails.
578 */
579 ESL_PERMUTATION *
esl_permutation_Create(int n)580 esl_permutation_Create(int n)
581 {
582 int status;
583 ESL_PERMUTATION *P = NULL;
584
585 ESL_DASSERT1(( n > 0 ));
586
587 ESL_ALLOC(P, sizeof(ESL_PERMUTATION));
588 P->pi = NULL;
589 P->n = n;
590 ESL_ALLOC(P->pi, sizeof(int) * n);
591
592 esl_permutation_Reuse(P); /* initialize it */
593 return P;
594
595 ERROR:
596 esl_permutation_Destroy(P);
597 return NULL;
598 }
599
600 /* Function: esl_permutation_Destroy()
601 *
602 * Purpose: Frees an <ESL_PERMUTATION> object <P>.
603 */
604 int
esl_permutation_Destroy(ESL_PERMUTATION * P)605 esl_permutation_Destroy(ESL_PERMUTATION *P)
606 {
607 if (P != NULL && P->pi != NULL) free(P->pi);
608 if (P != NULL) free(P);
609 return eslOK;
610 }
611
612 /* Function: esl_permutation_Reuse()
613 *
614 * Purpose: Resets a permutation matrix <P> to
615 * $p_i = i$ for all $i = 0..n-1$.
616 *
617 * Returns: <eslOK> on success.
618 */
619 int
esl_permutation_Reuse(ESL_PERMUTATION * P)620 esl_permutation_Reuse(ESL_PERMUTATION *P)
621 {
622 int i;
623 for (i = 0; i < P->n; i++)
624 P->pi[i] = i;
625 return eslOK;
626 }
627
628
629 /*****************************************************************
630 * 5. Debugging/validation for ESL_PERMUTATION.
631 *****************************************************************/
632
633 /* Function: esl_permutation_Dump()
634 *
635 * Purpose: Given a permutation matrix <P>, dump it to output stream <ofp>
636 * in human-readable format.
637 *
638 * If <rowlabel> or <collabel> are non-NULL, they represent
639 * single-character labels to put on the rows and columns,
640 * respectively. (For example, these might be a sequence
641 * alphabet for a 4x4 or 20x20 rate matrix or substitution
642 * matrix.) Numbers 1..ncols or 1..nrows are used if
643 * <collabel> or <rowlabel> are NULL.
644 *
645 * Args: ofp - output file pointer; stdout, for example
646 * P - permutation matrix to dump
647 * rowlabel - optional: NULL, or character labels for rows
648 * collabel - optional: NULL, or character labels for cols
649 *
650 * Returns: <eslOK> on success.
651 */
652 int
esl_permutation_Dump(FILE * ofp,const ESL_PERMUTATION * P,const char * rowlabel,const char * collabel)653 esl_permutation_Dump(FILE *ofp, const ESL_PERMUTATION *P, const char *rowlabel, const char *collabel)
654 {
655 int i,j;
656
657 fprintf(ofp, " ");
658 if (collabel != NULL)
659 for (j = 0; j < P->n; j++) fprintf(ofp, " %c ", collabel[j]);
660 else
661 for (j = 0; j < P->n; j++) fprintf(ofp, "%3d ", j+1);
662 fprintf(ofp, "\n");
663
664 for (i = 0; i < P->n; i++) {
665 if (rowlabel != NULL) fprintf(ofp, " %c ", rowlabel[i]);
666 else fprintf(ofp, "%3d ", i+1);
667
668 for (j = 0; j < P->n; j++)
669 fprintf(ofp, "%3d ", (j == P->pi[i]) ? 1 : 0);
670 fprintf(ofp, "\n");
671 }
672 return eslOK;
673 }
674
675 /*****************************************************************
676 * 6. The rest of the dmatrix API.
677 *****************************************************************/
678
679
680
681 /* Function: esl_dmx_Max()
682 * Incept: SRE, Thu Mar 1 14:46:48 2007 [Janelia]
683 *
684 * Purpose: Returns the maximum value of all the elements $a_{ij}$ in matrix <A>.
685 */
686 double
esl_dmx_Max(const ESL_DMATRIX * A)687 esl_dmx_Max(const ESL_DMATRIX *A)
688 {
689 int i;
690 double best;
691
692 best = A->mx[0][0];
693 for (i = 0; i < A->ncells; i++)
694 if (A->mx[0][i] > best) best = A->mx[0][i];
695 return best;
696 }
697
698 /* Function: esl_dmx_Min()
699 * Incept: SRE, Thu Mar 1 14:49:29 2007 [Janelia]
700 *
701 * Purpose: Returns the minimum value of all the elements $a_{ij}$ in matrix <A>.
702 */
703 double
esl_dmx_Min(const ESL_DMATRIX * A)704 esl_dmx_Min(const ESL_DMATRIX *A)
705 {
706 int i;
707 double best;
708
709 best = A->mx[0][0];
710 for (i = 0; i < A->ncells; i++)
711 if (A->mx[0][i] < best) best = A->mx[0][i];
712 return best;
713 }
714
715
716 /* Function: esl_dmx_MinMax()
717 * Incept: SRE, Wed Mar 14 16:58:03 2007 [Janelia]
718 *
719 * Purpose: Finds the maximum and minimum values of the
720 * elements $a_{ij}$ in matrix <A>, and returns
721 * them in <ret_min> and <ret_max>.
722 *
723 * Returns: <eslOK> on success.
724 *
725 */
726 int
esl_dmx_MinMax(const ESL_DMATRIX * A,double * ret_min,double * ret_max)727 esl_dmx_MinMax(const ESL_DMATRIX *A, double *ret_min, double *ret_max)
728 {
729 double min, max;
730 int i;
731
732 min = max = A->mx[0][0];
733 for (i = 0; i < A->ncells; i++) {
734 if (A->mx[0][i] < min) min = A->mx[0][i];
735 if (A->mx[0][i] > max) max = A->mx[0][i];
736 }
737 *ret_min = min;
738 *ret_max = max;
739 return eslOK;
740 }
741
742
743
744 /* Function: esl_dmx_Sum()
745 * Incept: SRE, Thu Mar 1 16:45:16 2007
746 *
747 * Purpose: Returns the scalar sum of all the elements $a_{ij}$ in matrix <A>,
748 * $\sum_{ij} a_{ij}$.
749 */
750 double
esl_dmx_Sum(const ESL_DMATRIX * A)751 esl_dmx_Sum(const ESL_DMATRIX *A)
752 {
753 int i;
754 double sum = 0.;
755
756 for (i = 0; i < A->ncells; i++)
757 sum += A->mx[0][i];
758 return sum;
759 }
760
761
762 /* Function: esl_dmx_FrobeniusNorm()
763 * Incept: SRE, Thu Mar 15 17:59:35 2007 [Janelia]
764 *
765 * Purpose: Calculates the Frobenius norm of a matrix, which
766 * is the element-wise equivalant of a
767 * Euclidean vector norm:
768 * $ = \sqrt(\sum a_{ij}^2)$
769 *
770 * Args: A - matrix
771 * ret_fnorm - Frobenius norm.
772 *
773 * Returns: <eslOK> on success, and the Frobenius norm
774 * is in <ret_fnorm>.
775 */
776 int
esl_dmx_FrobeniusNorm(const ESL_DMATRIX * A,double * ret_fnorm)777 esl_dmx_FrobeniusNorm(const ESL_DMATRIX *A, double *ret_fnorm)
778 {
779 double F = 0.;
780 int i;
781
782 for (i = 0; i < A->ncells; i++)
783 F += A->mx[0][i] * A->mx[0][i];
784 *ret_fnorm = sqrt(F);
785 return eslOK;
786 }
787
788
789 /* Function: esl_dmx_Multiply()
790 *
791 * Purpose: Matrix multiplication: calculate <AB>, store result in <C>.
792 * <A> is $n times m$; <B> is $m \times p$; <C> is $n \times p$.
793 * Matrix <C> must be allocated appropriately by the caller.
794 *
795 * Not supported for anything but general (<eslGENERAL>)
796 * matrix type, at present.
797 *
798 * Throws: <eslEINVAL> if matrices don't have compatible dimensions,
799 * or if any of them isn't a general (<eslGENERAL>) matrix.
800 */
801 int
esl_dmx_Multiply(const ESL_DMATRIX * A,const ESL_DMATRIX * B,ESL_DMATRIX * C)802 esl_dmx_Multiply(const ESL_DMATRIX *A, const ESL_DMATRIX *B, ESL_DMATRIX *C)
803 {
804 int i, j, k;
805
806 if (A->m != B->n) ESL_EXCEPTION(eslEINVAL, "can't multiply A,B");
807 if (A->n != C->n) ESL_EXCEPTION(eslEINVAL, "A,C # of rows not equal");
808 if (B->m != C->m) ESL_EXCEPTION(eslEINVAL, "B,C # of cols not equal");
809 if (A->type != eslGENERAL) ESL_EXCEPTION(eslEINVAL, "A isn't of type eslGENERAL");
810 if (B->type != eslGENERAL) ESL_EXCEPTION(eslEINVAL, "B isn't of type eslGENERAL");
811 if (C->type != eslGENERAL) ESL_EXCEPTION(eslEINVAL, "B isn't of type eslGENERAL");
812
813 /* i,k,j order should optimize stride, relative to a more textbook
814 * order for the indices
815 */
816 esl_dmatrix_SetZero(C);
817 for (i = 0; i < A->n; i++)
818 for (k = 0; k < A->m; k++)
819 for (j = 0; j < B->m; j++)
820 C->mx[i][j] += A->mx[i][k] * B->mx[k][j];
821
822 return eslOK;
823 }
824
825
826 /*::cexcerpt::function_comment_example::begin::*/
827 /* Function: esl_dmx_Exp()
828 * Synopsis: Calculates matrix exponential $\mathbf{P} = e^{t\mathbf{Q}}$.
829 * Incept: SRE, Thu Mar 8 18:41:38 2007 [Janelia]
830 *
831 * Purpose: Calculates the matrix exponential $\mathbf{P} = e^{t\mathbf{Q}}$,
832 * using a scaling and squaring algorithm with
833 * the Taylor series approximation \citep{MolerVanLoan03}.
834 *
835 * <Q> must be a square matrix of type <eslGENERAL>.
836 * Caller provides an allocated <P> matrix of the same size and type as <Q>.
837 *
838 * A typical use of this function is to calculate a
839 * conditional substitution probability matrix $\mathbf{P}$
840 * (whose elements $P_{xy}$ are conditional substitution
841 * probabilities $\mathrm{Prob}(y \mid x, t)$ from time $t$
842 * and instantaneous rate matrix $\mathbf{Q}$.
843 *
844 * Args: Q - matrix to exponentiate (an instantaneous rate matrix)
845 * t - time units
846 * P - RESULT: $e^{tQ}$.
847 *
848 * Returns: <eslOK> on success.
849 *
850 * Throws: <eslEMEM> on allocation error.
851 *
852 * Xref: J1/19.
853 */
854 int
esl_dmx_Exp(const ESL_DMATRIX * Q,double t,ESL_DMATRIX * P)855 esl_dmx_Exp(const ESL_DMATRIX *Q, double t, ESL_DMATRIX *P)
856 {
857 /*::cexcerpt::function_comment_example::end::*/
858 ESL_DMATRIX *Qz = NULL; /* Q/2^z rescaled matrix*/
859 ESL_DMATRIX *Qpow = NULL; /* keeps running product Q^k */
860 ESL_DMATRIX *C = NULL; /* tmp storage for matrix multiply result */
861 double factor = 1.0;
862 double fnorm;
863 int z;
864 double zfac;
865 int k;
866 int status;
867
868 /* Contract checks */
869 if (Q->type != eslGENERAL) ESL_EXCEPTION(eslEINVAL, "Q isn't general");
870 if (Q->n != Q->m) ESL_EXCEPTION(eslEINVAL, "Q isn't square");
871 if (P->type != Q->type) ESL_EXCEPTION(eslEINVAL, "P isn't of same type as Q");
872 if (P->n != P->m) ESL_EXCEPTION(eslEINVAL, "P isn't square");
873 if (P->n != Q->n) ESL_EXCEPTION(eslEINVAL, "P isn't same size as Q");
874
875 /* Allocation of working space */
876 if ((Qz = esl_dmatrix_Create(Q->n, Q->n)) == NULL) { status = eslEMEM; goto ERROR; }
877 if ((Qpow = esl_dmatrix_Create(Q->n, Q->n)) == NULL) { status = eslEMEM; goto ERROR; }
878 if ((C = esl_dmatrix_Create(Q->n, Q->n)) == NULL) { status = eslEMEM; goto ERROR; }
879
880 /* Figure out how much to scale the matrix down by. This is not
881 * magical; we're just knocking its magnitude down in an ad hoc way.
882 */
883 esl_dmx_FrobeniusNorm(Q, &fnorm);
884 zfac = 1.;
885 z = 0;
886 while (t*fnorm*zfac > 0.1) { zfac /= 2.; z++; }
887
888 /* Make a scaled-down copy of Q in Qz.
889 */
890 esl_dmatrix_Copy(Q, Qz);
891 esl_dmx_Scale(Qz, zfac);
892
893 /* Calculate e^{t Q_z} by the Taylor, to complete convergence. */
894 esl_dmatrix_SetIdentity(P);
895 esl_dmatrix_Copy(Qz, Qpow); /* Qpow is now Qz^1 */
896 for (k = 1; k < 100; k++)
897 {
898 factor *= t/k;
899 esl_dmatrix_Copy(P, C); /* C now holds the previous P */
900 esl_dmx_AddScale(P, factor, Qpow); /* P += factor*Qpow */
901 if (esl_dmatrix_Compare(C, P, 0.) == eslOK) break;
902
903 esl_dmx_Multiply(Qpow, Qz, C); /* C = Q^{k+1} */
904 esl_dmatrix_Copy(C, Qpow); /* Qpow = C = Q^{k+1} */
905 }
906
907 /* Now square it back up: e^{tQ} = [e^{tQ_z}]^{2^z} */
908 while (z--) {
909 esl_dmx_Multiply(P, P, C);
910 esl_dmatrix_Copy(C, P);
911 }
912
913 esl_dmatrix_Destroy(Qz);
914 esl_dmatrix_Destroy(Qpow);
915 esl_dmatrix_Destroy(C);
916 return eslOK;
917
918 ERROR:
919 if (Qz != NULL) esl_dmatrix_Destroy(Qz);
920 if (Qpow != NULL) esl_dmatrix_Destroy(Qpow);
921 if (C != NULL) esl_dmatrix_Destroy(C);
922 return status;
923 }
924
925
926 /* Function: esl_dmx_Transpose()
927 *
928 * Purpose: Transpose a square matrix <A> in place.
929 *
930 * <A> must be a general (<eslGENERAL>) matrix type.
931 *
932 * Throws: <eslEINVAL> if <A> isn't square, or if it isn't
933 * of type <eslGENERAL>.
934 */
935 int
esl_dmx_Transpose(ESL_DMATRIX * A)936 esl_dmx_Transpose(ESL_DMATRIX *A)
937 {
938 int i,j;
939 double swap;
940
941 if (A->n != A->m) ESL_EXCEPTION(eslEINVAL, "matrix isn't square");
942 if (A->type != eslGENERAL) ESL_EXCEPTION(eslEINVAL, "A isn't of type eslGENERAL");
943
944 for (i = 0; i < A->n; i++)
945 for (j = i+1; j < A->m; j++)
946 { swap = A->mx[i][j]; A->mx[i][j] = A->mx[j][i]; A->mx[j][i] = swap; }
947 return eslOK;
948 }
949
950
951 /* Function: esl_dmx_Add()
952 *
953 * Purpose: <A = A+B>; adds matrix <B> to matrix <A> and leaves result
954 * in matrix <A>.
955 *
956 * <A> and <B> may be of any type. However, if <A> is a
957 * packed upper triangular matrix (type
958 * <eslUPPER>), all values $i>j$ in <B> must be
959 * zero (i.e. <B> must also be upper triangular, though
960 * not necessarily packed upper triangular).
961 *
962 * Throws: <eslEINVAL> if matrices aren't the same dimensions, or
963 * if <A> is <eslUPPER> and any cell $i>j$ in
964 * <B> is nonzero.
965 */
966 int
esl_dmx_Add(ESL_DMATRIX * A,const ESL_DMATRIX * B)967 esl_dmx_Add(ESL_DMATRIX *A, const ESL_DMATRIX *B)
968 {
969 int i,j;
970
971 if (A->n != B->n) ESL_EXCEPTION(eslEINVAL, "matrices of different size");
972 if (A->m != B->m) ESL_EXCEPTION(eslEINVAL, "matrices of different size");
973
974 if (A->type == B->type) /* in this case, can just add cell by cell */
975 {
976 for (i = 0; i < A->ncells; i++)
977 A->mx[0][i] += B->mx[0][i];
978 }
979 else if (A->type == eslUPPER || B->type == eslUPPER)
980 {
981 /* Logic is: if either matrix is upper triangular, then the operation is
982 * to add upper triangles only. If we try to add a general matrix <B>
983 * to packed UT <A>, make sure all lower triangle entries in <B> are zero.
984 */
985 if (B->type != eslUPPER) {
986 for (i = 1; i < A->n; i++)
987 for (j = 0; j < i; j++)
988 if (B->mx[i][j] != 0.) ESL_EXCEPTION(eslEINVAL, "<B> has nonzero cells in lower triangle");
989 }
990 for (i = 0; i < A->n; i++)
991 for (j = i; j < A->m; j++)
992 A->mx[i][j] += B->mx[i][j];
993 }
994 return eslOK;
995 }
996
997 /* Function: esl_dmx_Scale()
998 *
999 * Purpose: Calculates <A = kA>: multiply matrix <A> by scalar
1000 * <k> and leave answer in <A>.
1001 */
1002 int
esl_dmx_Scale(ESL_DMATRIX * A,double k)1003 esl_dmx_Scale(ESL_DMATRIX *A, double k)
1004 {
1005 int i;
1006
1007 for (i = 0; i < A->ncells; i++) A->mx[0][i] *= k;
1008 return eslOK;
1009 }
1010
1011
1012 /* Function: esl_dmx_AddScale()
1013 *
1014 * Purpose: Calculates <A + kB>, leaves answer in <A>.
1015 *
1016 * Only defined for matrices of the same type (<eslGENERAL>
1017 * or <eslUPPER>).
1018 *
1019 * Throws: <eslEINVAL> if matrices aren't the same dimensions, or
1020 * of different types.
1021 */
1022 int
esl_dmx_AddScale(ESL_DMATRIX * A,double k,const ESL_DMATRIX * B)1023 esl_dmx_AddScale(ESL_DMATRIX *A, double k, const ESL_DMATRIX *B)
1024 {
1025 int i;
1026
1027 if (A->n != B->n) ESL_EXCEPTION(eslEINVAL, "matrices of different size");
1028 if (A->m != B->m) ESL_EXCEPTION(eslEINVAL, "matrices of different size");
1029 if (A->type != B->type) ESL_EXCEPTION(eslEINVAL, "matrices of different type");
1030
1031 for (i = 0; i < A->ncells; i++) A->mx[0][i] += k * B->mx[0][i];
1032 return eslOK;
1033 }
1034
1035
1036 /* Function: esl_dmx_Permute_PA()
1037 *
1038 * Purpose: Computes <B = PA>: do a row-wise permutation of a square
1039 * matrix <A>, using the permutation matrix <P>, and put
1040 * the result in a square matrix <B> that the caller has
1041 * allocated.
1042 *
1043 * Throws: <eslEINVAL> if <A>, <B>, <P> do not have compatible dimensions,
1044 * or if <A> or <B> is not of type <eslGENERAL>.
1045 */
1046 int
esl_dmx_Permute_PA(const ESL_PERMUTATION * P,const ESL_DMATRIX * A,ESL_DMATRIX * B)1047 esl_dmx_Permute_PA(const ESL_PERMUTATION *P, const ESL_DMATRIX *A, ESL_DMATRIX *B)
1048 {
1049 int i,ip,j;
1050
1051 if (A->n != P->n) ESL_EXCEPTION(eslEINVAL, "matrix dimensions not compatible");
1052 if (A->n != B->n) ESL_EXCEPTION(eslEINVAL, "matrix dimensions not compatible");
1053 if (A->n != A->m) ESL_EXCEPTION(eslEINVAL, "matrix dimensions not compatible");
1054 if (B->n != B->m) ESL_EXCEPTION(eslEINVAL, "matrix dimensions not compatible");
1055 if (A->type != eslGENERAL) ESL_EXCEPTION(eslEINVAL, "matrix A not of type eslGENERAL");
1056 if (B->type != eslGENERAL) ESL_EXCEPTION(eslEINVAL, "matrix B not of type eslGENERAL");
1057
1058 for (i = 0; i < A->n; i++)
1059 {
1060 ip = P->pi[i];
1061 for (j = 0; j < A->m; j++)
1062 B->mx[i][j] = A->mx[ip][j];
1063 }
1064 return eslOK;
1065 }
1066
1067 /* Function: esl_dmx_LUP_decompose()
1068 *
1069 * Purpose: Calculates a permuted LU decomposition of square matrix
1070 * <A>; upon return, <A> is replaced by this decomposition,
1071 * where <U> is in the lower triangle (inclusive of the
1072 * diagonal) and <L> is the upper triangle (exclusive of
1073 * diagonal, which is 1's by definition), and <P> is the
1074 * permutation matrix. Caller provides an allocated
1075 * permutation matrix <P> compatible with the square matrix
1076 * <A>.
1077 *
1078 * Implements Gaussian elimination with pivoting
1079 * \citep[p.~759]{Cormen99}.
1080 *
1081 * Throws: <eslEINVAL> if <A> isn't square, or if <P> isn't the right
1082 * size for <A>, or if <A> isn't of general type.
1083 */
1084 int
esl_dmx_LUP_decompose(ESL_DMATRIX * A,ESL_PERMUTATION * P)1085 esl_dmx_LUP_decompose(ESL_DMATRIX *A, ESL_PERMUTATION *P)
1086 {
1087 int i,j,k;
1088 int kpiv = 0; // initialization serves to quiet overzealous static analyzers
1089 double max;
1090 double swap;
1091
1092 if (A->n != A->m) ESL_EXCEPTION(eslEINVAL, "matrix isn't square");
1093 if (P->n != A->n) ESL_EXCEPTION(eslEINVAL, "permutation isn't the right size");
1094 if (A->type != eslGENERAL) ESL_EXCEPTION(eslEINVAL, "matrix isn't of general type");
1095 esl_permutation_Reuse(P);
1096
1097 for (k = 0; k < A->n-1; k++)
1098 {
1099 /* Identify our pivot; find row with maximum value in col[k].
1100 * This is guaranteed to succeed and set <kpiv>
1101 * (no matter what a static analyzer tells you)
1102 */
1103 max = 0.;
1104 for (i = k; i < A->n; i++)
1105 if (fabs(A->mx[i][k]) > max) {
1106 max = fabs(A->mx[i][k]);
1107 kpiv = i;
1108 }
1109 if (max == 0.) ESL_EXCEPTION(eslEDIVZERO, "matrix is singular");
1110
1111 /* Swap those rows (k and kpiv);
1112 * and keep track of that permutation in P. (misuse j for swapping integers)
1113 */
1114 j = P->pi[k]; P->pi[k] = P->pi[kpiv]; P->pi[kpiv] = j;
1115 for (j = 0; j < A->m; j++)
1116 { swap = A->mx[k][j]; A->mx[k][j] = A->mx[kpiv][j]; A->mx[kpiv][j] = swap; }
1117
1118 /* Gaussian elimination for all rows k+1..n.
1119 */
1120 for (i = k+1; i < A->n; i++)
1121 {
1122 A->mx[i][k] /= A->mx[k][k];
1123 for (j = k+1; j < A->m; j++)
1124 A->mx[i][j] -= A->mx[i][k] * A->mx[k][j];
1125 }
1126 }
1127 return eslOK;
1128 }
1129
1130
1131 /* Function: esl_dmx_LU_separate()
1132 *
1133 * Purpose: Separate a square <LU> decomposition matrix into its two
1134 * triangular matrices <L> and <U>. Caller provides two
1135 * allocated <L> and <U> matrices of same size as <LU> for
1136 * storing the results.
1137 *
1138 * <U> may be an upper triangular matrix in either unpacked
1139 * (<eslGENERAL>) or packed (<eslUPPER>) form.
1140 * <LU> and <L> must be of <eslGENERAL> type.
1141 *
1142 * Throws: <eslEINVAL> if <LU>, <L>, <U> are not of compatible dimensions,
1143 * or if <LU> or <L> aren't of general type.
1144 */
1145 int
esl_dmx_LU_separate(const ESL_DMATRIX * LU,ESL_DMATRIX * L,ESL_DMATRIX * U)1146 esl_dmx_LU_separate(const ESL_DMATRIX *LU, ESL_DMATRIX *L, ESL_DMATRIX *U)
1147 {
1148 int i,j;
1149
1150 if (LU->n != LU->m) ESL_EXCEPTION(eslEINVAL, "LU isn't square");
1151 if (L->n != L->m) ESL_EXCEPTION(eslEINVAL, "L isn't square");
1152 if (U->n != U->m) ESL_EXCEPTION(eslEINVAL, "U isn't square");
1153 if (LU->n != L->n) ESL_EXCEPTION(eslEINVAL, "LU, L have incompatible dimensions");
1154 if (LU->n != U->n) ESL_EXCEPTION(eslEINVAL, "LU, U have incompatible dimensions");
1155 if (LU->type != eslGENERAL) ESL_EXCEPTION(eslEINVAL, "matrix isn't of general type");
1156 if (L->type != eslGENERAL) ESL_EXCEPTION(eslEINVAL, "matrix isn't of general type");
1157
1158 esl_dmatrix_SetZero(L);
1159 esl_dmatrix_SetZero(U);
1160
1161 for (i = 0; i < LU->n; i++)
1162 for (j = i; j < LU->m; j++)
1163 U->mx[i][j] = LU->mx[i][j];
1164
1165 for (i = 0; i < LU->n; i++)
1166 {
1167 L->mx[i][i] = 1.;
1168 for (j = 0; j < i; j++)
1169 L->mx[i][j] = LU->mx[i][j];
1170 }
1171 return eslOK;
1172 }
1173
1174 /* Function: esl_dmx_Invert()
1175 *
1176 * Purpose: Calculates the inverse of square matrix <A>, and stores the
1177 * result in matrix <Ai>. Caller provides an allocated
1178 * matrix <Ai> of same dimensions as <A>. Both must be
1179 * of type <eslGENERAL>.
1180 *
1181 * Peforms the inversion by LUP decomposition followed by
1182 * forward/back-substitution \citep[p.~753]{Cormen99}.
1183 *
1184 * Throws: <eslEINVAL> if <A>, <Ai> do not have same dimensions,
1185 * if <A> isn't square, or if either isn't of
1186 * type <eslGENERAL>.
1187 * <eslEMEM> if internal allocations (for LU, and some other
1188 * bookkeeping) fail.
1189 */
1190 int
esl_dmx_Invert(const ESL_DMATRIX * A,ESL_DMATRIX * Ai)1191 esl_dmx_Invert(const ESL_DMATRIX *A, ESL_DMATRIX *Ai)
1192 {
1193 ESL_DMATRIX *LU = NULL;
1194 ESL_PERMUTATION *P = NULL;
1195 double *y = NULL; /* column vector, intermediate calculation */
1196 double *b = NULL; /* column vector of permuted identity matrix */
1197 int i,j,k;
1198 int status;
1199
1200 if (A->n != A->m) ESL_EXCEPTION(eslEINVAL, "matrix isn't square");
1201 if (A->n != Ai->n || A->m != Ai->m) ESL_EXCEPTION(eslEINVAL, "matrices are different size");
1202 if (A->type != eslGENERAL) ESL_EXCEPTION(eslEINVAL, "matrix A not of general type");
1203 if (Ai->type != eslGENERAL) ESL_EXCEPTION(eslEINVAL, "matrix B not of general type");
1204
1205 /* Copy A to LU, and do an LU decomposition.
1206 */
1207 if ((LU = esl_dmatrix_Create(A->n, A->m)) == NULL) { status = eslEMEM; goto ERROR; }
1208 if ((P = esl_permutation_Create(A->n)) == NULL) { status = eslEMEM; goto ERROR; }
1209 if (( status = esl_dmatrix_Copy(A, LU)) != eslOK) goto ERROR;
1210 if (( status = esl_dmx_LUP_decompose(LU, P)) != eslOK) goto ERROR;
1211
1212 /* Now we have:
1213 * PA = LU
1214 *
1215 * to invert a matrix A, we want A A^-1 = I;
1216 * that's PAx = Pb, for columns x of A^-1 and b of the identity matrix;
1217 * and that's n equations LUx = Pb;
1218 *
1219 * so, solve Ly = Pb for y by forward substitution;
1220 * then Ux = y by back substitution;
1221 * x is then a column of A^-1.
1222 *
1223 * Do that for all columns.
1224 */
1225 ESL_ALLOC(b, sizeof(double) * A->n);
1226 ESL_ALLOC(y, sizeof(double) * A->n);
1227
1228 for (k = 0; k < A->m; k++) /* for each column... */
1229 {
1230 /* build Pb for column j of the identity matrix */
1231 for (i = 0; i < A->n; i++)
1232 if (P->pi[i] == k) b[i] = 1.; else b[i] = 0.;
1233
1234 /* forward substitution
1235 */
1236 for (i = 0; i < A->n; i++)
1237 {
1238 y[i] = b[i];
1239 for (j = 0; j < i; j++) y[i] -= LU->mx[i][j] * y[j];
1240 }
1241
1242 /* back substitution
1243 */
1244 for (i = A->n-1; i >= 0; i--)
1245 {
1246 Ai->mx[i][k] = y[i];
1247 for (j = i+1; j < A->n; j++) Ai->mx[i][k] -= LU->mx[i][j] * Ai->mx[j][k];
1248 Ai->mx[i][k] /= LU->mx[i][i];
1249 }
1250 }
1251
1252 free(b);
1253 free(y);
1254 esl_dmatrix_Destroy(LU);
1255 esl_permutation_Destroy(P);
1256 return eslOK;
1257
1258 ERROR:
1259 if (y != NULL) free(y);
1260 if (b != NULL) free(b);
1261 if (LU != NULL) esl_dmatrix_Destroy(LU);
1262 if (P != NULL) esl_permutation_Destroy(P);
1263 return status;
1264 }
1265
1266
1267 /*****************************************************************
1268 * 7. Optional: interoperability with GSL
1269 *****************************************************************/
1270 #ifdef HAVE_LIBGSL
1271
1272 #include <gsl/gsl_matrix.h>
1273
1274 int
esl_dmx_MorphGSL(const ESL_DMATRIX * E,gsl_matrix ** ret_G)1275 esl_dmx_MorphGSL(const ESL_DMATRIX *E, gsl_matrix **ret_G)
1276 {
1277 gsl_matrix *G = NULL;
1278 int i,j;
1279
1280 if (E->type != eslGENERAL) ESL_EXCEPTION(eslEINVAL, "can only morph general matrices to GSL right now");
1281
1282 G = gsl_matrix_alloc(E->m, E->n);
1283 for (i = 0; i < E->m; i++)
1284 for (j = 0; j < E->n; j++)
1285 gsl_matrix_set(G, i, j, E->mx[i][j]);
1286 *ret_G = G;
1287 return eslOK;
1288 }
1289
1290 int
esl_dmx_UnmorphGSL(const gsl_matrix * G,ESL_DMATRIX ** ret_E)1291 esl_dmx_UnmorphGSL(const gsl_matrix *G, ESL_DMATRIX **ret_E)
1292 {
1293 ESL_DMATRIX *E = NULL;
1294 int i,j;
1295
1296 if ((E = esl_dmatrix_Create(G->size1, G->size2)) == NULL) return eslEMEM;
1297 for (i = 0; i < G->size1; i++)
1298 for (j = 0; j < G->size2; j++)
1299 E->mx[i][j] = gsl_matrix_get(G, i, j);
1300 *ret_E = E;
1301 return eslOK;
1302 }
1303
1304 #endif /*HAVE_LIBGSL*/
1305
1306 /*****************************************************************
1307 * 8. Optional: Interfaces to LAPACK
1308 *****************************************************************/
1309 #ifdef HAVE_LIBLAPACK
1310
1311 /* To include LAPACK code, you need to:
1312 * 1. declare the C interface to the Fortran routine,
1313 * appending _ to the Fortran routine's name (dgeev becomes dgeev_)
1314 *
1315 * 2. Remember to transpose matrices into column-major
1316 * Fortran form
1317 *
1318 * 3. everything must be passed by reference, not by value
1319 *
1320 * 4. you don't need any include files, just lapack.a
1321 *
1322 * 5. Add -llapack to the compile line.
1323 * (It doesn't appear that blas or g2c are needed?)
1324 */
1325
1326 /* Declare the C interface to the Fortran77 dgeev routine
1327 * provided by the LAPACK library:
1328 */
1329 extern void dgeev_(char *jobvl, char *jobvr, int *n, double *a,
1330 int *lda, double *wr, double *wi, double *vl,
1331 int *ldvl, double *vr, int *ldvr,
1332 double *work, int *lwork, int *info);
1333
1334
1335 /* Function: esl_dmx_Diagonalize()
1336 * Incept: SRE, Thu Mar 15 09:28:03 2007 [Janelia]
1337 *
1338 * Purpose: Given a square real matrix <A>, diagonalize it:
1339 * solve for $U^{-1} A U = diag(\lambda_1... \lambda_n)$.
1340 *
1341 * Upon return, <ret_Er> and <ret_Ei> are vectors
1342 * containing the real and complex parts of the eigenvalues
1343 * $\lambda_i$; <ret_UL> is the $U^{-1}$ matrix containing
1344 * the left eigenvectors; and <ret_UR> is the $U$ matrix
1345 * containing the right eigenvectors.
1346 *
1347 * <ret_UL> and <ret_UR> are optional; pass <NULL> for
1348 * either if you don't want that set of eigenvectors.
1349 *
1350 * This is a C interface to the <dgeev()> routine in the
1351 * LAPACK linear algebra library.
1352 *
1353 * Args: A - square nxn matrix to diagonalize
1354 * ret_Er - RETURN: real part of eigenvalues (0..n-1)
1355 * ret_Ei - RETURN: complex part of eigenvalues (0..n-1)
1356 * ret_UL - optRETURN: nxn matrix of left eigenvectors
1357 * ret_UR - optRETURN:
1358 *
1359 * Returns: <eslOK> on success.
1360 * <ret_Er> and <ret_Ei> (and <ret_UL>,<ret_UR> when they are
1361 * requested) are allocated here, and must be free'd by the caller.
1362 *
1363 * Throws: <eslEMEM> on allocation failure.
1364 * In this case, the four return pointers are returned <NULL>.
1365 *
1366 * Xref: J1/19.
1367 */
1368 int
esl_dmx_Diagonalize(const ESL_DMATRIX * A,double ** ret_Er,double ** ret_Ei,ESL_DMATRIX ** ret_UL,ESL_DMATRIX ** ret_UR)1369 esl_dmx_Diagonalize(const ESL_DMATRIX *A, double **ret_Er, double **ret_Ei,
1370 ESL_DMATRIX **ret_UL, ESL_DMATRIX **ret_UR)
1371 {
1372 int status;
1373 double *Er = NULL;
1374 double *Ei = NULL;
1375 ESL_DMATRIX *At = NULL;
1376 ESL_DMATRIX *UL = NULL;
1377 ESL_DMATRIX *UR = NULL;
1378 double *work = NULL;
1379 char jobul, jobur;
1380 int lda;
1381 int ldul, ldur;
1382 int lwork;
1383 int info;
1384
1385 if (A->n != A->m) ESL_EXCEPTION(eslEINVAL, "matrix isn't square");
1386
1387 if ((At = esl_dmatrix_Clone(A)) == NULL) { status = eslEMEM; goto ERROR; }
1388 if ((UL = esl_dmatrix_Create(A->n,A->n)) == NULL) { status = eslEMEM; goto ERROR; }
1389 if ((UR = esl_dmatrix_Create(A->n,A->n)) == NULL) { status = eslEMEM; goto ERROR; }
1390 ESL_ALLOC(Er, sizeof(double) * A->n);
1391 ESL_ALLOC(Ei, sizeof(double) * A->n);
1392 ESL_ALLOC(work, sizeof(double) * 8 * A->n);
1393
1394 jobul = (ret_UL == NULL) ? 'N' : 'V'; /* do we want left eigenvectors? */
1395 jobur = (ret_UR == NULL) ? 'N' : 'V'; /* do we want right eigenvectors? */
1396 lda = A->n;
1397 ldul = A->n;
1398 ldur = A->n;
1399 lwork = 8*A->n;
1400
1401 /* Fortran convention is colxrow, not rowxcol; so transpose
1402 * a copy of A before passing it to a Fortran routine.
1403 */
1404 esl_dmx_Transpose(At);
1405
1406 /* The Fortran77 interface call to LAPACK's dgeev().
1407 * All args must be passed by reference.
1408 * Fortran 2D arrays are 1D: so pass the A[0] part of a DSMX.
1409 */
1410 dgeev_(&jobul, &jobur, &(At->n), At->mx[0], &lda, Er, Ei,
1411 UL->mx[0], &ldul, UR->mx[0], &ldur, work, &lwork, &info);
1412 if (info < 0) ESL_XEXCEPTION(eslEINVAL, "argument %d to LAPACK dgeev is invalid", -info);
1413 if (info > 0) ESL_XEXCEPTION(eslEINVAL,
1414 "diagonalization failed; only eigenvalues %d..%d were computed",
1415 info+1, At->n);
1416
1417 /* Now, UL, UR are transposed (col x row), so transpose them back to
1418 * C language convention.
1419 */
1420 esl_dmx_Transpose(UL);
1421 esl_dmx_Transpose(UR);
1422
1423 esl_dmatrix_Destroy(At);
1424 if (ret_UL != NULL) *ret_UL = UL; else esl_dmatrix_Destroy(UL);
1425 if (ret_UR != NULL) *ret_UR = UR; else esl_dmatrix_Destroy(UR);
1426 if (ret_Er != NULL) *ret_Er = Er; else free(Er);
1427 if (ret_Ei != NULL) *ret_Ei = Ei; else free(Ei);
1428 free(work);
1429 return eslOK;
1430
1431 ERROR:
1432 if (ret_UL != NULL) *ret_UL = NULL;
1433 if (ret_UR != NULL) *ret_UR = NULL;
1434 if (ret_Er != NULL) *ret_Er = NULL;
1435 if (ret_Ei != NULL) *ret_Ei = NULL;
1436 if (At != NULL) esl_dmatrix_Destroy(At);
1437 if (UL != NULL) esl_dmatrix_Destroy(UL);
1438 if (UR != NULL) esl_dmatrix_Destroy(UR);
1439 if (Er != NULL) free(Er);
1440 if (Ei != NULL) free(Ei);
1441 if (work != NULL) free(work);
1442 return status;
1443 }
1444
1445
1446 #endif /*HAVE_LIBLAPACK*/
1447
1448 /*****************************************************************
1449 * 9. Unit tests
1450 *****************************************************************/
1451 #ifdef eslDMATRIX_TESTDRIVE
1452
1453 static void
utest_misc_ops(void)1454 utest_misc_ops(void)
1455 {
1456 char *msg = "miscellaneous unit test failed";
1457 ESL_DMATRIX *A, *B, *C;
1458 int n = 20;
1459
1460 if ((A = esl_dmatrix_Create(n,n)) == NULL) esl_fatal(msg);
1461 if ((B = esl_dmatrix_Create(n,n)) == NULL) esl_fatal(msg);
1462 if ((C = esl_dmatrix_Create(n,n)) == NULL) esl_fatal(msg);
1463
1464 if (esl_dmatrix_SetIdentity(A) != eslOK) esl_fatal(msg); /* A=I */
1465 if (esl_dmx_Invert(A, B) != eslOK) esl_fatal(msg); /* B=I^-1=I */
1466 if (esl_dmx_Multiply(A,B,C) != eslOK) esl_fatal(msg); /* C=I */
1467 if (esl_dmx_Transpose(A) != eslOK) esl_fatal(msg); /* A=I still */
1468
1469 if (esl_dmx_Scale(A, 0.5) != eslOK) esl_fatal(msg); /* A= 0.5I */
1470 if (esl_dmx_AddScale(B, -0.5, C) != eslOK) esl_fatal(msg); /* B= 0.5I */
1471
1472 if (esl_dmx_Add(A,B) != eslOK) esl_fatal(msg); /* A=I */
1473 if (esl_dmx_Scale(B, 2.0) != eslOK) esl_fatal(msg); /* B=I */
1474
1475 if (esl_dmatrix_Compare(A, B, 1e-12) != eslOK) esl_fatal(msg);
1476 if (esl_dmatrix_Compare(A, C, 1e-12) != eslOK) esl_fatal(msg);
1477 if (esl_dmatrix_Copy(B, C) != eslOK) esl_fatal(msg);
1478 if (esl_dmatrix_Compare(A, C, 1e-12) != eslOK) esl_fatal(msg);
1479
1480 esl_dmatrix_Destroy(A);
1481 esl_dmatrix_Destroy(B);
1482 esl_dmatrix_Destroy(C);
1483 return;
1484 }
1485
1486
1487 static void
utest_Invert(ESL_DMATRIX * A)1488 utest_Invert(ESL_DMATRIX *A)
1489 {
1490 char *msg = "Failure in matrix inversion unit test";
1491 ESL_DMATRIX *Ai = NULL;
1492 ESL_DMATRIX *B = NULL;
1493 ESL_DMATRIX *I = NULL;
1494
1495 if ((Ai = esl_dmatrix_Create(A->n, A->m)) == NULL) esl_fatal(msg);
1496 if ((B = esl_dmatrix_Create(A->n, A->m)) == NULL) esl_fatal(msg);
1497 if ((I = esl_dmatrix_Create(A->n, A->m)) == NULL) esl_fatal(msg);
1498 if (esl_dmx_Invert(A, Ai) != eslOK) esl_fatal("matrix inversion failed");
1499 if (esl_dmx_Multiply(A,Ai,B) != eslOK) esl_fatal(msg);
1500 if (esl_dmatrix_SetIdentity(I) != eslOK) esl_fatal(msg);
1501 if (esl_dmatrix_Compare(B,I, 1e-12) != eslOK) esl_fatal("inverted matrix isn't right");
1502
1503 esl_dmatrix_Destroy(Ai);
1504 esl_dmatrix_Destroy(B);
1505 esl_dmatrix_Destroy(I);
1506 return;
1507 }
1508
1509
1510 #endif /*eslDMATRIX_TESTDRIVE*/
1511
1512
1513
1514 /*****************************************************************
1515 * 10. Test driver
1516 *****************************************************************/
1517
1518 /* gcc -g -Wall -o test -I. -L. -DeslDMATRIX_TESTDRIVE esl_dmatrix.c -leasel -lm
1519 */
1520 #ifdef eslDMATRIX_TESTDRIVE
1521 #include "easel.h"
1522 #include "esl_dmatrix.h"
1523 #include "esl_random.h"
1524
main(void)1525 int main(void)
1526 {
1527 ESL_RANDOMNESS *r;
1528 ESL_DMATRIX *A;
1529 int n = 30;
1530 int seed = 42;
1531 int i,j;
1532 double range = 100.;
1533
1534 /* Create a square matrix with random values from -range..range */
1535 if ((r = esl_randomness_Create(seed)) == NULL) esl_fatal("failed to create random source");
1536 if ((A = esl_dmatrix_Create(n, n)) == NULL) esl_fatal("failed to create matrix");
1537 for (i = 0; i < n; i++)
1538 for (j = 0; j < n; j++)
1539 A->mx[i][j] = esl_random(r) * range * 2.0 - range;
1540
1541 utest_misc_ops();
1542 utest_Invert(A);
1543
1544 esl_randomness_Destroy(r);
1545 esl_dmatrix_Destroy(A);
1546 return 0;
1547 }
1548 #endif /*eslDMATRIX_TESTDRIVE*/
1549
1550
1551 /*****************************************************************
1552 * 11. Examples
1553 *****************************************************************/
1554
1555 /* gcc -g -Wall -o example -I. -DeslDMATRIX_EXAMPLE esl_dmatrix.c easel.c -lm
1556 */
1557 #ifdef eslDMATRIX_EXAMPLE
1558 /*::cexcerpt::dmatrix_example::begin::*/
1559 #include "easel.h"
1560 #include "esl_dmatrix.h"
1561
main(void)1562 int main(void)
1563 {
1564 ESL_DMATRIX *A, *B, *C;
1565
1566 A = esl_dmatrix_Create(4,4);
1567 B = esl_dmatrix_Create(4,4);
1568 C = esl_dmatrix_Create(4,4);
1569
1570 esl_dmatrix_SetIdentity(A);
1571 esl_dmatrix_Copy(A, B);
1572
1573 esl_dmx_Multiply(A,B,C);
1574
1575 esl_dmatrix_Dump(stdout, C, NULL, NULL);
1576
1577 esl_dmatrix_Destroy(A);
1578 esl_dmatrix_Destroy(B);
1579 esl_dmatrix_Destroy(C);
1580 return 0;
1581 }
1582 /*::cexcerpt::dmatrix_example::end::*/
1583 #endif /*eslDMATRIX_EXAMPLE*/
1584
1585
1586
1587
1588