JKQTPlotter trunk/v5.0.0
an extensive Qt5+Qt6 Plotter framework (including a feature-richt plotter widget, a speed-optimized, but limited variant and a LaTeX equation renderer!), written fully in C/C++ and without external dependencies
Loading...
Searching...
No Matches
jkqtplinalgtools.h
1/*
2 Copyright (c) 2008-2024 Jan W. Krieger (<jan@jkrieger.de>)
3
4 last modification: $LastChangedDate$ (revision $Rev$)
5
6 This software is free software: you can redistribute it and/or modify
7 it under the terms of the GNU Lesser General Public License (LGPL) as published by
8 the Free Software Foundation, either version 2.1 of the License, or
9 (at your option) any later version.
10
11 This program is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU Lesser General Public License (LGPL) for more details.
15
16 You should have received a copy of the GNU Lesser General Public License (LGPL)
17 along with this program. If not, see <http://www.gnu.org/licenses/>.
18*/
19
20
21#ifndef JKQTPLINALGTOOLS_H_INCLUDED
22#define JKQTPLINALGTOOLS_H_INCLUDED
23
24#include <stdint.h>
25#include <cmath>
26#include <stdlib.h>
27#include <string.h>
28#include <iostream>
29#include <stdio.h>
30#include <stdlib.h>
31#include <limits>
32#include <vector>
33#include <utility>
34#include <cfloat>
35#include <ostream>
36#include <iomanip>
37#include <sstream>
38#include "jkqtmath/jkqtmath_imexport.h"
39#include "jkqtmath/jkqtparraytools.h"
40#include "jkqtcommon/jkqtpmathtools.h"
41#include "jkqtcommon/jkqtpstringtools.h"
42
43#ifdef _OPENMP
44# include <omp.h>
45#endif
46
47
48#ifndef JKQTP_ALIGNMENT_BYTES
49#define JKQTP_ALIGNMENT_BYTES 32
50#endif
51
52
53#ifdef JKQTP_STATISTICS_TOOLS_MAY_USE_EIGEN3
54# include <Eigen/Core>
55# include <Eigen/SVD>
56# include <Eigen/Jacobi>
57# include <Eigen/LU>
58# include <Eigen/QR>
59#endif
60
61
62
63
64
65/** \brief calculate the index of the entry in line \a l and column \a c
66 * in a row-major matrix with \a C columns
67 * \ingroup jkqtptools_math_linalg
68 *
69 * You can use this to access a matrix with L rows and C columns:
70 * \code
71 * for (long l=0; l<L; l++) {
72 * for (long c=0; c<C; c++) {
73 * matrix[jkqtplinalgMatIndex(l,c,C)=0;
74 * if (l==c) matrix[jkqtplinalgMatIndex(l,c,C)=1;
75 * }
76 * }
77 * \endcode
78 */
79#define jkqtplinalgMatIndex(l,c,C) ((l)*(C)+(c))
80
81
82/** \brief print the given LxC matrix to std::cout
83 * \ingroup jkqtptools_math_linalg
84 *
85 * \tparam type of the matrix cells (typically double or float)
86 * \param matrix the matrix to print out
87 * \param L number of lines/rows in the matrix
88 * \param C number of columns in the matrix
89 * \param width width (in characters) of each cell in the output (used for formatting)
90 * \param precision precision (in digits) for string-conversions in the output (used for formatting)
91 * \param mode if \c =='f' the mode \c std::fixed is used for output, otherwise \c std::scientific is used
92*/
93template <class T>
94inline void jkqtplinalgPrintMatrix(const T* matrix, long L, long C, int width=9, int precision=3, char mode='f') {
95 for (long l=0; l<L; l++) {
96 for (long c=0; c<C; c++) {
97 if (c>0) std::cout<<", ";
98 std::cout.precision(precision);
99 std::cout.width(width);
100 if (mode=='f') std::cout<<std::fixed<<std::right<<matrix[jkqtplinalgMatIndex(l,c,C)];
101 else std::cout<<std::scientific<<std::right<<matrix[jkqtplinalgMatIndex(l,c,C)];
102 }
103 std::cout<<std::endl;
104 }
105}
106
107
108
109/** \brief convert the given LxC matrix to std::string
110 * \ingroup jkqtptools_math_linalg
111 *
112 * \tparam type of the matrix cells (typically double or float)
113 * \param matrix the matrix to convert
114 * \param L number of lines/rows in the matrix
115 * \param C number of columns in the matrix
116 * \param width width (in characters) of each cell in the output (used for formatting)
117 * \param precision precision (in digits) for string-conversions in the output (used for formatting)
118 * \param mode the (printf()) string conversion mode for output of the cell values
119*/
120template <class T>
121inline std::string jkqtplinalgMatrixToString(T* matrix, long L, long C, int width=9, int precision=3, const std::string& mode=std::string("g")) {
122 std::string format="%"+jkqtp_inttostr(width)+std::string(".")+jkqtp_inttostr(precision)+std::string("l")+mode;
123 std::ostringstream ost;
124 for (long l=0; l<L; l++) {
125 for (long c=0; c<C; c++) {
126 if (c>0) ost<<", ";
127 char buf[500];
128 snprintf(buf, 500, format.c_str(), jkqtp_todouble(matrix[jkqtplinalgMatIndex(l,c,C)]));
129 ost<<buf;
130 /*ost.precision(precision);
131 ost.width(width);
132 if (mode=='f') ost<<std::fixed<<std::right<<matrix[jkqtplinalgMatIndex(l,c,C)];
133 else ost<<std::scientific<<std::right<<matrix[jkqtplinalgMatIndex(l,c,C)];*/
134
135 }
136 ost<<std::endl;
137 }
138 return ost.str();
139}
140
141/** \brief maps the number range -1 ... +1 to a color-scale lightblue - white - lightred (used for coloring matrices!)
142 * \ingroup jkqtptools_math_linalg
143 *
144 * \param val the value to convert
145 * \param[out] r returns the red value (0..255)
146 * \param[out] g returns the green value (0..255)
147 * \param[out] b returns the blue value (0..255)
148 */
149inline void jkqtplinalgPM1ToRWBColors(double val, uint8_t& r, uint8_t& g, uint8_t& b){
150 r=255;
151 g=255;
152 b=255;
153 const double fval=fabs(val);
154 if (val<0 && val>=-1) {
155 r=jkqtp_boundedRoundTo<uint8_t>(0,255.0-fval*127.0,255);
156 g=jkqtp_boundedRoundTo<uint8_t>(0,255.0-fval*127.0,255);
157 } else if (val>0 && val<=1) {
158 b=jkqtp_boundedRoundTo<uint8_t>(0,255.0-fval*127.0,255);
159 g=jkqtp_boundedRoundTo<uint8_t>(0,255.0-fval*127.0,255);
160 } else if (val<-1) {
161 r=127;
162 g=127;
163 b=255;
164 } else if (val>1) {
165 r=255;
166 g=127;
167 b=127;
168 }
169}
170
171/** \brief maps the number range -1 ... +1 to a non-linear color-scale lightblue - white - lightred (used for coloring matrices!)
172 * \ingroup jkqtptools_math_linalg
173 *
174 * \param val the value to convert
175 * \param[out] r returns the red value (0..255)
176 * \param[out] g returns the green value (0..255)
177 * \param[out] b returns the blue value (0..255)
178 * \param gamma a gamma-value for non-linear color scaling
179*/
180inline void jkqtplinalgPM1ToNonlinRWBColors(double val, uint8_t& r, uint8_t& g, uint8_t& b, double gamma=0.5){
181 if (val<0) {
182 jkqtplinalgPM1ToRWBColors(-1.0*pow(-val,gamma),r,g,b);
183 } else {
184 jkqtplinalgPM1ToRWBColors(pow(val,gamma),r,g,b);
185 }
186}
187/** \brief convert the given LxC matrix to std::string, encoded as HTML table
188 * \ingroup jkqtptools_math_linalg
189 *
190 *
191 * \tparam type of the matrix cells (typically double or float)
192 * \param matrix the matrix to convert
193 * \param L number of lines/rows in the matrix
194 * \param C number of columns in the matrix
195 * \param width width (in characters) of each cell in the output (used for formatting)
196 * \param precision precision (in digits) for string-conversions in the output (used for formatting)
197 * \param mode the (printf()) string conversion mode for output of the cell values
198 * \param tableformat this is inserted into the \c <table...> tag (may contain HTML property definitions)
199 * \param prenumber this is inserted before each number (may contain HTML markup)
200 * \param postnumber this is inserted after each number (may contain HTML markup)
201 * \param colorcoding if \c true, teh cell backgrounds are color-coded
202 * \param zeroIsWhite if \c the color-coding is forced to white for 0 and then encodes in positive/negative direction with colors (red and blue)
203 * \param[out] colorlabel outputs a label explaining the auto-generated color-coding
204 * \param nonlinColors if \c true, a non-linear color-coding is used
205 * \param nonlinColorsGamma gamma-value for a non-linear color-coding
206 * \param colortableformat lie \a tableformat, but for the legend table output in \a colorLabel
207 *
208 * \see jkqtplinalgPM1ToNonlinRWBColors() and jkqtplinalgPM1ToRWBColors()
209*/
210template <class T>
211inline std::string jkqtplinalgMatrixToHTMLString(T* matrix, long L, long C, int width=9, int precision=3, const std::string& mode=std::string("g"), const std::string& tableformat=std::string(), const std::string& prenumber=std::string(), const std::string& postnumber=std::string(), bool colorcoding=false, bool zeroIsWhite=true, std::string* colorlabel=nullptr, bool nonlinColors=false, double nonlinColorsGamma=0.25, const std::string& colortableformat=std::string()) {
212 std::ostringstream ost;
213 ost<<"<table "<<tableformat<<">\n";
214 std::string format="%"+jkqtp_inttostr(width)+std::string(".")+jkqtp_inttostr(precision)+std::string("l")+mode;
215 double minv=0, maxv=0;
216 if (colorcoding) {
217 jkqtpstatMinMax(matrix, L*C, minv, maxv);
218 }
219 for (long l=0; l<L; l++) {
220 ost<<" <tr>";
221 for (long c=0; c<C; c++) {
222 const double val=matrix[jkqtplinalgMatIndex(l,c,C)];
223 std::string cols="";
224 if (colorcoding) {
225 double valrel=0;
226 uint8_t r=255,g=255,b=255;
227 if (zeroIsWhite){
228 if (val<0) valrel=-1.0*fabs(val/minv);
229 if (val>0) valrel=fabs(val/maxv);
230 } else {
231 valrel=((val-minv)/(maxv-minv)-0.5)*2.0;
232 }
233 if (nonlinColors) {
234 jkqtplinalgPM1ToNonlinRWBColors(valrel, r,g,b, nonlinColorsGamma);
235 } else {
236 jkqtplinalgPM1ToRWBColors(valrel, r,g,b);
237 }
238 char buf[500];
239 snprintf(buf, 500, " bgcolor=\"#%02X%02X%02X\"", int(r), int(g), int(b));
240 cols=buf;
241 }
242 ost<<"<td align=\"center\" valign=\"middle\" width=\""<<(100.0/double(C))<<"%\" "<<cols<<"><nobr>";
243 ost.precision(precision);
244 ost.width(width);
245 char buf[500];
246 snprintf(buf, 500, format.c_str(), val);
247 ost<<prenumber<<buf<<postnumber;
248 ost<<"</nobr></td>";
249 }
250 ost<<"</tr>\n";
251 }
252 ost<<"</table>";
253
254 if (colorcoding && colorlabel) {
255 char buf[8192];
256 uint8_t rm=255,gm=255,bm=255;
257 uint8_t rmc=255,gmc=255,bmc=255;
258 uint8_t rc=255,gc=255,bc=255;
259 uint8_t rcp=255,gcp=255,bcp=255;
260 uint8_t rp=255,gp=255,bp=255;
261 double vm=minv;
262 double vc=0;
263 double vp=maxv;
264 if (!zeroIsWhite) {
265 vc=(maxv+minv)/2.0;
266 }
267 if (nonlinColors) {
268 jkqtplinalgPM1ToNonlinRWBColors(-1, rm, gm, bm, nonlinColorsGamma);
269 jkqtplinalgPM1ToNonlinRWBColors(-0.5, rmc, gmc, bmc, nonlinColorsGamma);
270 jkqtplinalgPM1ToNonlinRWBColors(0, rc, gc, bc, nonlinColorsGamma);
271 jkqtplinalgPM1ToNonlinRWBColors(0.5, rcp, gcp, bcp, nonlinColorsGamma);
272 jkqtplinalgPM1ToNonlinRWBColors(1, rp, gp, bp, nonlinColorsGamma);
273 } else {
274 jkqtplinalgPM1ToRWBColors(-1, rm, gm, bm);
275 jkqtplinalgPM1ToRWBColors(-0.5, rmc, gmc, bmc);
276 jkqtplinalgPM1ToRWBColors(0, rc, gc, bc);
277 jkqtplinalgPM1ToRWBColors(0.5, rcp, gcp, bcp);
278 jkqtplinalgPM1ToRWBColors(1, rp, gp, bp);
279 }
280 snprintf(buf, 8192, "<table %s cellpadding=\"2\" cellspacing=\"0\" border=\"1\"><tr><td><table width=\"100%%\" cellpadding=\"3\" cellspacing=\"0\" border=\"0\"><tr>"
281 "<td bgcolor=\"#%02X%02X%02X\" width=\"20%%\"><nobr>&nbsp;%9.3lg&nbsp;</nobr></td>"
282 "<td bgcolor=\"#%02X%02X%02X\" width=\"20%%\"><nobr>&nbsp;&nbsp;&nbsp;&mdash;&nbsp;&nbsp;&nbsp;</nobr></td>"
283 "<td bgcolor=\"#%02X%02X%02X\" width=\"20%%\"><nobr>&nbsp;%9.3lg&nbsp;</nobr></td>"
284 "<td bgcolor=\"#%02X%02X%02X\" width=\"20%%\"><nobr>&nbsp;&nbsp;&nbsp;&mdash;&nbsp;&nbsp;&nbsp;</nobr></td>"
285 "<td bgcolor=\"#%02X%02X%02X\" width=\"20%%\"><nobr>&nbsp;%9.3lg&nbsp;</nobr></td>"
286 "</tr></table></td></tr></table>", colortableformat.c_str(), int(rm), int(gm), int(bm), vm, int(rmc), int(gmc), int(bmc), int(rc), int(gc), int(bc), vc, int(rcp), int(gcp), int(bcp), int(rp), int(gp), int(bp), vp);
287 (*colorlabel)=std::string(buf);
288 }
289 return ost.str();
290}
291
292/** \brief dot-product between two vectors \a vec1 and \a vec2, each with a length of \a N entries
293 * \ingroup jkqtptools_math_linalg
294 *
295 * \tparam T of the vector cells (typically double or float)
296 * \param vec1 first vector
297 * \param vec2 second vector
298 * \param N number of entries in \a vec1 and \a vec2
299 */
300template <class T>
301inline T jkqtplinalgDotProduct(const T* vec1, const T* vec2, long N) {
302 T res=0;
303 for (long l=0; l<N; l++) {
304 res=res+vec1[l]*vec2[l];
305 }
306 return res;
307}
308
309
310/** \brief transpose the given NxN matrix
311 * \ingroup jkqtptools_math_linalg
312 *
313 * \tparam T of the matrix cells (typically double or float)
314 * \param matrix the matrix to transpose
315 * \param N number of rows and columns in the matrix
316 *
317 */
318template <class T>
319inline void jkqtplinalgTransposeMatrix(T* matrix, long N) {
320 for (long l=0; l<N; l++) {
321 for (long c=l+1; c<N; c++) {
323 }
324 }
325}
326
327
328/** \brief transpose the given LxC matrix
329 * \ingroup jkqtptools_math_linalg
330 *
331 * \tparam T of the matrix cells (typically double or float)
332 * \param matrix the matrix to transpose
333 * \param L number of rows in the matrix
334 * \param C number of columns in the matrix
335 *
336 * \note The output is interpreted as CxL matrix!!!
337 */
338template <class T>
339inline void jkqtplinalgTransposeMatrix(T* matrix, long L, long C) {
340 JKQTPArrayScopedPointer<T> t(jkqtpArrayDuplicate<T>(matrix, L*C));
341 for (long l=0; l<L; l++) {
342 for (long c=0; c<C; c++) {
343 matrix[jkqtplinalgMatIndex(c,l,L)]=t[jkqtplinalgMatIndex(l,c,C)];
344 }
345 }
346}
347
348/** \brief swap two lines in a matrix
349 * \ingroup jkqtptools_math_linalg
350 *
351 * \tparam T of the matrix cells (typically double or float)
352 * \param m the matrix to work on
353 * \param l1 the row to swap with row \a l2
354 * \param l2 the row to swap with row \a l1
355 * \param C number of columns in the matrix
356*/
357template <class T>
358inline void jkqtplinalgMatrixSwapLines(T* m, long l1, long l2, long C) {
359 for (long c=0; c<C; c++) {
361 }
362}
363
364
365/*! \brief matrix-matrix product
366 \ingroup jkqtptools_math_linalg
367
368 \tparam T of the matrix cells (typically double or float)
369 \param M1 matrix 1, size: L1xC1
370 \param L1 number of rows in the matrix \a M1
371 \param C1 number of columns in the matrix \a M1
372 \param M2 matrix 1, size: L2xC2
373 \param L2 number of rows in the matrix \a M2
374 \param C2 number of columns in the matrix \a M2
375 \param[out] M output matrix M=M1*M2, size: L1xC2
376*/
377template <class T>
378inline void jkqtplinalgMatrixProduct(const T* M1, long L1, long C1, const T* M2, long L2, long C2, T* M) {
379 if (M1!=M &&M2!=M) {
380
381 for (long l=0; l<L1; l++) {
382 for (long c=0; c<C2; c++) {
383 double s=T(0);
384 for (long i=0; i<C1; i++) {
385 s = s + M1[jkqtplinalgMatIndex(l, i, C1)]*M2[jkqtplinalgMatIndex(i,c, C2)];
386 }
387 M[jkqtplinalgMatIndex(l, c, C2)]=s;
388 }
389 }
390 } else if (M1==M && M2!=M) {
392 jkqtplinalgMatrixProduct(MM.data(),L1,C1,M2,L2,C2,M);
393 } else if (M1!=M && M2==M) {
395 jkqtplinalgMatrixProduct(M1,L1,C1,MM.data(),L2,C2,M);
396 } else if (M1==M && M2==M) {
398 jkqtplinalgMatrixProduct(MM.data(),L1,C1,MM.data(),L2,C2,M);
399 }
400}
401
402
403/*! \brief matrix-matrix product of two NxN matrices
404 \ingroup jkqtptools_math_linalg
405
406 \param M1 matrix 1, size: NxN
407 \param M2 matrix 1, size: NxN
408 \param N number os rows/columns in the matrices \a M1, \a M2 and \a M
409 \param[out] M output matrix M=M1*M2, size: NxN
410*/
411template <class T>
412inline void jkqtplinalgMatrixProduct(const T* M1, const T* M2, long N, T* M) {
413 jkqtplinalgMatrixProduct(M1,N,N,M2,N,N,M);
414}
415
416
417
418
419
420
421/*! \brief performs a Gauss-Jordan eliminaion on a LxC matrix m
422 \ingroup jkqtptools_math_linalg
423
424 For a matrix equation \f[ A\cdot\vec{x}=\vec{b} \f] the input matrix is \f[ \left[A | \vec{b}\right] \f] for matrix inversion it is
425 \f[ \left[A | I_L\right] \f] where \f$ I_L \f$ is the unit matrix with LxL entries.
426
427 \tparam T of the matrix cells (typically double or float)
428 \param m the matrix
429 \param L number of rows in the matrix
430 \param C number of columns in the matrix
431
432 \see http://www.virtual-maxim.de/matrix-invertieren-in-c-plus-plus/
433*/
434template <class T>
435inline bool jkqtplinalgGaussJordan(T* m, long L, long C) {
436
437 const long N=L;
438
439 //std::cout<<"\nstep 0:\n";
440 //linalgPrintMatrix(m, N, C);
441
442 // first we perform a Gauss-elimination, which transforms the matrix in the left half of m into upper triangular form
443 for (long k=0; k<N-1; k++) {
444 //std::cout<<"\nstep G"<<k<<": pivot="<<m[jkqtpstatisticsMatIndex(k,k,C)]<<"\n";
445 if (m[jkqtplinalgMatIndex(k,k,C)]==0.0) {
446 // if pivot(m[k,k])==0, then swap this line with a line, which does not have a 0 in the k-th column
447 for (long ks=k+1; ks<N; ks++) {
448 if (m[jkqtplinalgMatIndex(ks,k,C)]!=0.0) {
449 jkqtplinalgMatrixSwapLines(m, ks, k, C);
450 break;
451 } else if (ks==N-1) {
452 // if no such element is found, the matrix may not be inverted!
453 return false;
454 }
455 }
456 }
457
458 // now we can eliminate all entries i below the pivot line p, by subtracting
459 // the pivot line, scaled by s, from every line, where
460 // s=m[i,p]/m[p,p]
461 for (long i=k+1; i<N; i++) {
462 const T s=m[jkqtplinalgMatIndex(i,k,C)]/m[jkqtplinalgMatIndex(k,k,C)];
463 for (long c=k; c<C; c++) {
464 m[jkqtplinalgMatIndex(i,c,C)] -= m[jkqtplinalgMatIndex(k,c,C)]*s;
465 }
466 }
467
468 //linalgPrintMatrix(m, N, C);
469 }
470
471 // now we can caluate the determinant of the left half-matrix, which can be used to determine, whether matrix
472 // is invertible at all: if det(T)==0.0 -> matrix is not invertible
473 // the determinant of an upper triangular marix equals the product of the diagonal elements
474 T det=1.0;
475 for (long k=0; k<N; k++) {
476 det = det * m[jkqtplinalgMatIndex(k,k,C)];
477 }
478 //linalgPrintMatrix(m, N, C);
479 //std::cout<<"\nstep 2: det(M)="<<det<<"\n";
480 if (fabs(det)<DBL_MIN*10.0) return false;
481
482
483 // if the matrix may be inverted, we can go on with the JOrdan part of the algorithm:
484 // we work the Nx(2N) matrix from bottom to top and transform the left side into a unit matrix
485 // - the last row is left unchanged
486 // - the last row is subtracted from every row i above, scaled by m[i,N]/m[N,N]
487 // then we repeat this for the (N-1)*(N-1) left upper matrix, which has again full triangular form
488 for (long k=N-1; k>0; k--) {
489 //std::cout<<"\nstep J"<<k<<":\n";
490 for (long i=k-1; i>=0; i--) {
491 const T s=m[jkqtplinalgMatIndex(i,k,C)]/m[jkqtplinalgMatIndex(k,k,C)];
492 for (long c=k; c<C; c++) {
493 m[jkqtplinalgMatIndex(i,c,C)] -= m[jkqtplinalgMatIndex(k,c,C)]*s;
494 }
495 }
496 //linalgPrintMatrix(m, N, C);
497 }
498 // finally each line is normalized to 1 by dividing by the diagonal entry in the left NxN matrix
499 // and copy the result to matrix_out
500 for (long k=0; k<N; k++) {
501
502 const T s=m[jkqtplinalgMatIndex(k,k,C)];
503 for (long c=N; c<C; c++) {
504 m[jkqtplinalgMatIndex(k,c,C)] = m[jkqtplinalgMatIndex(k,c,C)] / s;
505 }
506 m[jkqtplinalgMatIndex(k,k,C)]=T(1);
507 }
508
509
510 return true;
511}
512
513
514
515
516
517
518
519
520/*! \brief invert the given NxN matrix using Gauss-Jordan elimination
521 \ingroup jkqtptools_math_linalg
522
523 \tparam T of the matrix cells (typically double or float)
524 \param matrix the matrix to invert
525 \param[out] matrix_out target for the inverted matrix
526 \param N number of rows and columns in the matrix
527
528 \return \c true on success and the inverted matrix in matrix_out.
529
530 \note It is save to call \c jkqtpstatMatrixInversion(A,A,N) with the same argument for in and out matrix. Then the input will be overwritten with the new matrix!
531 \note matrix and matrix_out have to be of size N*N. Matrices are interpreted as row-major!
532
533 \see http://www.virtual-maxim.de/matrix-invertieren-in-c-plus-plus/
534*/
535template <class T>
536inline bool jkqtplinalgMatrixInversion(const T* matrix, T* matrix_out, long N) {
537#ifdef JKQTP_STATISTICS_TOOLS_MAY_USE_EIGEN3
538 if (N==2) {
539 Eigen::Map<const Eigen::Matrix<T,2,2,Eigen::RowMajor> > eA(matrix);
540 Eigen::Map<Eigen::Matrix<T,2,2,Eigen::RowMajor> > eO(matrix_out);
541 eO=eA.inverse();
542 //std::cout<<"\n--------------------------------------\n2x2 input matrix\n"<<eA<<"\n--------------------------------------\n";
543 return eO.allFinite();
544 } else if (N==3) {
545 Eigen::Map<const Eigen::Matrix<T,3,3,Eigen::RowMajor> > eA(matrix);
546 Eigen::Map<Eigen::Matrix<T,3,3,Eigen::RowMajor> > eO(matrix_out);
547 //std::cout<<"\n--------------------------------------\n3x3 input matrix\n"<<eA<<"\n--------------------------------------\n";
548 eO=eA.inverse();
549 return eO.allFinite();
550 } else if (N==4) {
551 Eigen::Map<const Eigen::Matrix<T,4,4,Eigen::RowMajor> > eA(matrix);
552 Eigen::Map<Eigen::Matrix<T,4,4,Eigen::RowMajor> > eO(matrix_out);
553 //std::cout<<"\n--------------------------------------\n4x4 input matrix\n"<<eA<<"\n--------------------------------------\n";
554 eO=eA.inverse();
555 return eO.allFinite();
556 } else {
557 Eigen::Map<const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic,Eigen::RowMajor> > eA(matrix,N,N);
558 Eigen::Map<Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic,Eigen::RowMajor> > eO(matrix_out,N,N);
559 eO=eA.inverse();
560 //std::cout<<"\n--------------------------------------\nNxN input matrix\n"<<eA<<"\n--------------------------------------\n";
561 return eO.allFinite();
562 }
563 return true;
564#else
565 // first build a N*(2N) matrix of the form
566 //
567 // m11 m12 ... | 1 0 0
568 // m21 m22 ... | 0 1 0
569 // ... ... ... | .....
570 //
571 const long msize=N*N*2;
572 std::vector<T> m;
573 m.resize(msize);
574 for (long i=0; i<msize; i++) m[i]=T(0); // init with 0
575 for (long l=0; l<N; l++) {
576 for (long c=0; c<N; c++) { // init left half with matrix
577 m[jkqtplinalgMatIndex(l,c,2*N)]=matrix[jkqtplinalgMatIndex(l,c,N)];
578 }
579 // init right half with unit matrix
580 m[jkqtplinalgMatIndex(l,N+l,2*N)]=T(1);
581 }
582
583
584 bool ok=jkqtplinalgGaussJordan(m.data(), N, 2*N);
585
586 if (ok) {
587 // finally we copy the result to matrix_out
588 for (long k=0; k<N; k++) {
589 for (long c=N; c<2*N; c++) {
590 matrix_out[jkqtplinalgMatIndex(k,c-N,N)] = m[jkqtplinalgMatIndex(k,c,2*N)];
591 }
592 }
593 }
594
595 return ok;
596#endif
597}
598
599
600
601
602
603/*! \brief invert the given NxN matrix using Gauss-Jordan elimination
604 \ingroup jkqtptools_math_linalg
605
606 \tparam T of the matrix cells (typically double or float)
607 \param[in,out] matrix the matrix to invert (at the same time the target)
608 \param N number of rows and columns in the matrix
609*/
610template <class T>
611inline bool jkqtplinalgMatrixInversion(T* matrix, long N) {
612 return jkqtplinalgMatrixInversion(matrix, matrix, N);
613}
614
615
616
617/*! \brief solve a system of N linear equations \f$ A\cdot\vec{x}=B \f$ simultaneously for C columns in B
618 \ingroup jkqtptools_math_linalg
619
620 \param A an NxN matrix of coefficients
621 \param B an NxC marix
622 \param N number of equations
623 \param C number of columns in B
624 \param result_out a NxC matrix with the results after the inversion of the system of equations
625 \return \c true on success
626
627 \note This function uses the Gauss-Jordan algorithm
628 \note It is save to call \c jkqtpstatLinSolve(A,B,N,C,B) with the same argument for B and result_out. Then the input will be overwritten with the new matrix!
629*/
630template <class T>
631inline bool jkqtplinalgLinSolve(const T* A, const T* B, long N, long C, T* result_out) {
632#if defined(JKQTP_STATISTICS_TOOLS_MAY_USE_EIGEN3) && (!defined(STATISTICS_TOOLS_linalgLinSolve_EIGENMETHOD_noeigen))
633 if (N==2 && C==1) {
634 Eigen::Map<const Eigen::Matrix<T,2,2,Eigen::RowMajor> > eA(A);
635 Eigen::Matrix<T,2,1> eB;
636 Eigen::Map<Eigen::Matrix<T,2,1> > x(result_out);
637
638 eB=Eigen::Map<const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic,Eigen::RowMajor> >(B,2,1);
639# ifdef STATISTICS_TOOLS_linalgLinSolve_EIGENMETHOD_fullPivLu
640 x=eA.fullPivLu().solve(eB);
641# elif defined(STATISTICS_TOOLS_linalgLinSolve_EIGENMETHOD_householderQr)
642 x=eA.householderQr().solve(eB);
643# elif defined(STATISTICS_TOOLS_linalgLinSolve_EIGENMETHOD_fullPivHouseholderQr)
644 x=eA.fullPivHouseholderQr().solve(eB);
645# elif defined(STATISTICS_TOOLS_linalgLinSolve_EIGENMETHOD_jacobiSvd)
646 x=eA.jacobiSVD().solve(eB);
647# elif defined(STATISTICS_TOOLS_linalgLinSolve_EIGENMETHOD_colPivHouseholderQr)
648 x=eA.colPivHouseholderQr().solve(eB);
649# else
650 x=eA.fullPivLu().solve(eB);
651# endif
652 } else if (N==3 && C==1) {
653 Eigen::Map<const Eigen::Matrix<T,3,3,Eigen::RowMajor> > eA(A);
654 Eigen::Matrix<T,3,1> eB;
655 Eigen::Map<Eigen::Matrix<T,3,1> > x(result_out);
656
657 eB=Eigen::Map<const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic,Eigen::RowMajor> >(B,3,1);
658
659# ifdef STATISTICS_TOOLS_linalgLinSolve_EIGENMETHOD_fullPivLu
660 x=eA.fullPivLu().solve(eB);
661# elif defined(STATISTICS_TOOLS_linalgLinSolve_EIGENMETHOD_householderQr)
662 x=eA.householderQr().solve(eB);
663# elif defined(STATISTICS_TOOLS_linalgLinSolve_EIGENMETHOD_fullPivHouseholderQr)
664 x=eA.fullPivHouseholderQr().solve(eB);
665# elif defined(STATISTICS_TOOLS_linalgLinSolve_EIGENMETHOD_jacobiSvd)
666 x=eA.jacobiSVD().solve(eB);
667# elif defined(STATISTICS_TOOLS_linalgLinSolve_EIGENMETHOD_colPivHouseholderQr)
668 x=eA.colPivHouseholderQr().solve(eB);
669# else
670 x=eA.fullPivLu().solve(eB);
671# endif
672 } else {
673 Eigen::Map<const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic,Eigen::RowMajor> > eA(A,N,N);
674 Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic,Eigen::RowMajor> eB(N,C);
675 Eigen::Map<Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic,Eigen::RowMajor> > x(result_out,N,C);
676
677 eB=Eigen::Map<const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic,Eigen::RowMajor> >(B,N,C);
678
679# ifdef STATISTICS_TOOLS_linalgLinSolve_EIGENMETHOD_fullPivLu
680 x=eA.fullPivLu().solve(eB);
681# elif defined(STATISTICS_TOOLS_linalgLinSolve_EIGENMETHOD_householderQr)
682 x=eA.householderQr().solve(eB);
683# elif defined(STATISTICS_TOOLS_linalgLinSolve_EIGENMETHOD_fullPivHouseholderQr)
684 x=eA.fullPivHouseholderQr().solve(eB);
685# elif defined(STATISTICS_TOOLS_linalgLinSolve_EIGENMETHOD_jacobiSvd)
686 x=eA.jacobiSVD().solve(eB);
687# elif defined(STATISTICS_TOOLS_linalgLinSolve_EIGENMETHOD_colPivHouseholderQr)
688 x=eA.colPivHouseholderQr().solve(eB);
689# else
690 x=eA.fullPivLu().solve(eB);
691# endif
692 }
693 return true;
694#else
695 // first build a N*(N+C) matrix of the form
696 //
697 // <---- N ----> <---- C ---->
698 // ^ A11 A12 ... | B11 B12 ...
699 // | A21 A22 ... | B21 B22 ...
700 // N ... ... ... | .....
701 // | ... ... ... | .....
702 // v ... ... ... | .....
703 //
704 const long msize=N*(N+C);
705 JKQTPArrayScopedPointer<T> m(static_cast<T*>(malloc(msize*sizeof(T)))); // use scoped pointer to ensure, that m is free'd, when the function is ending
706 for (long l=0; l<N; l++) {
707 for (long c=0; c<N; c++) { // init left half with matrix A
708 m[jkqtplinalgMatIndex(l,c,N+C)]=A[jkqtplinalgMatIndex(l,c,N)];
709 }
710 // init right half with B
711 for (long c=0; c<C; c++) { // init left half with matrix B
712 m[jkqtplinalgMatIndex(l,N+c,N+C)]=B[jkqtplinalgMatIndex(l,c,C)];
713 }
714 }
715
716
717 bool ok=jkqtplinalgGaussJordan(m.data(), N, N+C);
718
719 if (ok) {
720 for (long k=0; k<N; k++) {
721 for (long c=N; c<(N+C); c++) {
722 if (!JKQTPIsOKFloat(m[jkqtplinalgMatIndex(k,c,N+C)])) {
723 ok=false;
724 break;
725 }
726 }
727 if (!ok) break;
728 }
729 if (ok) {
730 // finally we copy the result to matrix_out
731 for (long k=0; k<N; k++) {
732 for (long c=N; c<(N+C); c++) {
733 result_out[jkqtplinalgMatIndex(k,c-N,C)] = m[jkqtplinalgMatIndex(k,c,N+C)];
734 }
735 }
736 }
737 }
738
739 return ok;
740#endif
741}
742
743/*! \brief solve a system of N linear equations \f$ A\cdot\vec{x}=B \f$ simultaneously for C columns in B
744 \ingroup jkqtptools_math_linalg
745
746 \param A an NxN matrix of coefficients
747 \param[in,out] B an NxC marix (also receives the result)
748 \param N number of equations
749 \param C number of columns in B
750 \return \c true on success
751
752 \note This function uses the Gauss-Jordan algorithm
753 \note It is save to call \c jkqtpstatLinSolve(A,B,N,C,B) with the same argument for B and result_out. Then the input will be overwritten with the new matrix!
754*/
755template <class T>
756inline bool jkqtplinalgLinSolve(const T* A, T* B, long N, long C) {
757 return jkqtplinalgLinSolve(A,B,N,C,B);
758}
759
760
761/*! \brief solve a system of N linear equations \f$ A\cdot\vec{x}=\vec{b} \f$
762 \ingroup jkqtptools_math_linalg
763
764 \param A an NxN matrix of coefficients
765 \param b an N-entry vector
766 \param N number of rows and columns in \a A
767 \param[out] result_out a N-entry vector with the result
768 \return \c true on success
769
770 \note This function uses the Gauss-Jordan algorithm
771 \note It is save to call \c jkqtpstatLinSolve(A,B,N,C,B) with the same argument for B and result_out. Then the input will be overwritten with the new matrix!
772*/
773template <class T>
774inline bool jkqtplinalgLinSolve(const T* A, const T* b, long N, T* result_out) {
775 return jkqtplinalgLinSolve(A, b, N, 1, result_out);
776}
777
778/*! \brief solve a system of N linear equations \f$ A\cdot\vec{x}=\vec{b} \f$
779 \ingroup jkqtptools_math_linalg
780
781 \param A an NxN matrix of coefficients
782 \param[in,out] b an N-entry vector (also receives the result)
783 \param N number of rows and columns in \a A
784 \return \c true on success
785
786 \note This function uses the Gauss-Jordan algorithm
787 \note It is save to call \c jkqtpstatLinSolve(A,B,N,C,B) with the same argument for B and result_out. Then the input will be overwritten with the new matrix!
788*/
789template <class T>
790inline bool jkqtplinalgLinSolve(const T* A, T* b, long N) {
791 return jkqtplinalgLinSolve(A,b,N,1,b);
792}
793
794
795
796
797/*! \brief determinant the given NxN matrix
798 \ingroup jkqtptools_math_linalg
799
800 \tparam T of the matrix cells (typically double or float)
801 \param a the matrix for which to calculate the determinant
802 \param N number of rows and columns in the matrix
803 \return the determinant of \a a
804
805 \note for large matrices this algorithm is very slow, think about defining the macro JKQTP_STATISTICS_TOOLS_MAY_USE_EIGEN3 to use faster methods from the EIGEN library!
806*/
807template <class T>
808inline T jkqtplinalgMatrixDeterminant(const T* a, long N) {
809#ifdef JKQTP_STATISTICS_TOOLS_MAY_USE_EIGEN3
810 if (N < 1) { /* Error */
811 return NAN;
812 } else if (N == 1) { /* Shouldn't get used */
813 return a[jkqtplinalgMatIndex(0,0,N)];
814 } else if (N == 2) {
815 return a[jkqtplinalgMatIndex(0,0,N)] * a[jkqtplinalgMatIndex(1,1,N)] - a[jkqtplinalgMatIndex(1,0,N)] * a[jkqtplinalgMatIndex(0,1,N)];
816 } else if (N==3){
817 Eigen::Map<const Eigen::Matrix<T,3,3,Eigen::RowMajor> > eA(a);
818 return eA.determinant();
819 } else if (N==4){
820 Eigen::Map<const Eigen::Matrix<T,4,4,Eigen::RowMajor> > eA(a);
821 return eA.determinant();
822 } else {
823 Eigen::Map<const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic,Eigen::RowMajor> > eA(a,N,N);
824 return eA.determinant();
825 }
826#else
827 long i,j,j1,j2;
828 T det = 0;
829
830 if (N < 1) { /* Error */
831 det=NAN;
832 } else if (N == 1) { /* Shouldn't get used */
833 det = a[jkqtplinalgMatIndex(0,0,N)];
834 } else if (N == 2) {
835 det = a[jkqtplinalgMatIndex(0,0,N)] * a[jkqtplinalgMatIndex(1,1,N)] - a[jkqtplinalgMatIndex(1,0,N)] * a[jkqtplinalgMatIndex(0,1,N)];
836 } else {
837 det = 0;
838 for (j1=0;j1<N;j1++) {
839 JKQTPArrayScopedPointer<T> m(static_cast<T*>(calloc((N-1)*(N-1),sizeof(T *))));
840
841 for (i=1;i<N;i++) {
842 j2 = 0;
843 for (j=0;j<N;j++) {
844 if (j != j1){
845 m[jkqtplinalgMatIndex(i-1,j2,N-1)] = a[jkqtplinalgMatIndex(i,j,N)];
846 j2++;
847 }
848 }
849 }
850 //printf("%d: %lf (%lf %lf)\n",j1,pow(-1.0,1.0+double(j1)+1.0),a[jkqtplinalgMatIndex(0,j1,N)], jkqtplinalgMatrixDeterminant(m,N-1));
851 det = det + double(((1+j1+1)%2==0)?1.0:-1.0)/* pow(-1.0,1.0+double(j1)+1.0)*/ * a[jkqtplinalgMatIndex(0,j1,N)] * jkqtplinalgMatrixDeterminant(m,N-1);
852
853 }
854 }
855
856 //printf("\n\n jkqtplinalgMatrixDeterminant(%d):\n",N);
857 //linalgPrintMatrix(a,N,N);
858 //printf(" jkqtplinalgMatrixDeterminant(%d) = %lf\n", N, det);
859 return(det);
860#endif
861}
862
863
864
865
866#endif // JKQTPLINALGTOOLS_H_INCLUDED
867
868
this class ensures that the given pointer is freed when the class is destroyed.
Definition jkqtparraytools.h:108
T * data() const
Definition jkqtparraytools.h:138
void jkqtpArraySwap(T *a, long long l, long long r)
swap two elements l and r in an array a
Definition jkqtparraytools.h:67
T * jkqtpArrayDuplicate(const T *dataIn, long long N)
duplicate an array of data
Definition jkqtparraytools.h:92
constexpr double jkqtp_todouble(const T &d)
converts a boolean to a double, is used to convert boolean to double by JKQTPDatastore
Definition jkqtpmathtools.h:113
bool JKQTPIsOKFloat(T v)
check whether the dlotaing point number is OK (i.e. non-inf, non-NAN)
Definition jkqtpmathtools.h:496
T jkqtplinalgMatrixDeterminant(const T *a, long N)
determinant the given NxN matrix
Definition jkqtplinalgtools.h:808
#define jkqtplinalgMatIndex(l, c, C)
calculate the index of the entry in line l and column c in a row-major matrix with C columns
Definition jkqtplinalgtools.h:79
void jkqtplinalgTransposeMatrix(T *matrix, long N)
transpose the given NxN matrix
Definition jkqtplinalgtools.h:319
bool jkqtplinalgGaussJordan(T *m, long L, long C)
performs a Gauss-Jordan eliminaion on a LxC matrix m
Definition jkqtplinalgtools.h:435
std::string jkqtplinalgMatrixToString(T *matrix, long L, long C, int width=9, int precision=3, const std::string &mode=std::string("g"))
convert the given LxC matrix to std::string
Definition jkqtplinalgtools.h:121
std::string jkqtplinalgMatrixToHTMLString(T *matrix, long L, long C, int width=9, int precision=3, const std::string &mode=std::string("g"), const std::string &tableformat=std::string(), const std::string &prenumber=std::string(), const std::string &postnumber=std::string(), bool colorcoding=false, bool zeroIsWhite=true, std::string *colorlabel=nullptr, bool nonlinColors=false, double nonlinColorsGamma=0.25, const std::string &colortableformat=std::string())
convert the given LxC matrix to std::string, encoded as HTML table
Definition jkqtplinalgtools.h:211
bool jkqtplinalgMatrixInversion(const T *matrix, T *matrix_out, long N)
invert the given NxN matrix using Gauss-Jordan elimination
Definition jkqtplinalgtools.h:536
void jkqtplinalgPrintMatrix(const T *matrix, long L, long C, int width=9, int precision=3, char mode='f')
print the given LxC matrix to std::cout
Definition jkqtplinalgtools.h:94
bool jkqtplinalgLinSolve(const T *A, const T *B, long N, long C, T *result_out)
solve a system of N linear equations simultaneously for C columns in B
Definition jkqtplinalgtools.h:631
void jkqtplinalgPM1ToRWBColors(double val, uint8_t &r, uint8_t &g, uint8_t &b)
maps the number range -1 ... +1 to a color-scale lightblue - white - lightred (used for coloring matr...
Definition jkqtplinalgtools.h:149
void jkqtplinalgMatrixProduct(const T *M1, long L1, long C1, const T *M2, long L2, long C2, T *M)
matrix-matrix product
Definition jkqtplinalgtools.h:378
void jkqtplinalgMatrixSwapLines(T *m, long l1, long l2, long C)
swap two lines in a matrix
Definition jkqtplinalgtools.h:358
T jkqtplinalgDotProduct(const T *vec1, const T *vec2, long N)
dot-product between two vectors vec1 and vec2, each with a length of N entries
Definition jkqtplinalgtools.h:301
void jkqtplinalgPM1ToNonlinRWBColors(double val, uint8_t &r, uint8_t &g, uint8_t &b, double gamma=0.5)
maps the number range -1 ... +1 to a non-linear color-scale lightblue - white - lightred (used for co...
Definition jkqtplinalgtools.h:180
void jkqtpstatMinMax(InputIt first, InputIt last, double &min, double &max, InputIt *minPos=nullptr, InputIt *maxPos=nullptr, size_t *Noutput=nullptr)
calculates the minimum and maximum values in the given data range first ... last
Definition jkqtpstatbasics.h:168
JKQTCOMMON_LIB_EXPORT std::string jkqtp_inttostr(long data)
convert an integer to a string