mirror of
https://github.com/jkriege2/JKQtPlotter.git
synced 2024-11-16 02:25:50 +08:00
869 lines
34 KiB
C++
869 lines
34 KiB
C++
/*
|
|
Copyright (c) 2008-2022 Jan W. Krieger (<jan@jkrieger.de>)
|
|
|
|
last modification: $LastChangedDate$ (revision $Rev$)
|
|
|
|
This software is free software: you can redistribute it and/or modify
|
|
it under the terms of the GNU Lesser General Public License (LGPL) as published by
|
|
the Free Software Foundation, either version 2.1 of the License, or
|
|
(at your option) any later version.
|
|
|
|
This program is distributed in the hope that it will be useful,
|
|
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
GNU Lesser General Public License (LGPL) for more details.
|
|
|
|
You should have received a copy of the GNU Lesser General Public License (LGPL)
|
|
along with this program. If not, see <http://www.gnu.org/licenses/>.
|
|
*/
|
|
|
|
|
|
#ifndef JKQTPLINALGTOOLS_H_INCLUDED
|
|
#define JKQTPLINALGTOOLS_H_INCLUDED
|
|
|
|
#include <stdint.h>
|
|
#include <cmath>
|
|
#include <stdlib.h>
|
|
#include <string.h>
|
|
#include <iostream>
|
|
#include <stdio.h>
|
|
#include <stdlib.h>
|
|
#include <limits>
|
|
#include <vector>
|
|
#include <utility>
|
|
#include <cfloat>
|
|
#include <ostream>
|
|
#include <iomanip>
|
|
#include <sstream>
|
|
#include "jkqtcommon_statistics_and_math/jkqtcommon_statistics_and_math_imexport.h"
|
|
#include "jkqtcommon_statistics_and_math/jkqtparraytools.h"
|
|
#include "jkqtcommon/jkqtpmathtools.h"
|
|
#include "jkqtcommon/jkqtpstringtools.h"
|
|
|
|
#ifdef _OPENMP
|
|
# include <omp.h>
|
|
#endif
|
|
|
|
|
|
#ifndef JKQTP_ALIGNMENT_BYTES
|
|
#define JKQTP_ALIGNMENT_BYTES 32
|
|
#endif
|
|
|
|
|
|
#ifdef JKQTP_STATISTICS_TOOLS_MAY_USE_EIGEN3
|
|
# include <Eigen/Core>
|
|
# include <Eigen/SVD>
|
|
# include <Eigen/Jacobi>
|
|
# include <Eigen/LU>
|
|
# include <Eigen/QR>
|
|
#endif
|
|
|
|
|
|
|
|
|
|
|
|
/** \brief calculate the index of the entry in line \a l and column \a c
|
|
* in a row-major matrix with \a C columns
|
|
* \ingroup jkqtptools_math_linalg
|
|
*
|
|
* You can use this to access a matrix with L rows and C columns:
|
|
* \code
|
|
* for (long l=0; l<L; l++) {
|
|
* for (long c=0; c<C; c++) {
|
|
* matrix[jkqtplinalgMatIndex(l,c,C)=0;
|
|
* if (l==c) matrix[jkqtplinalgMatIndex(l,c,C)=1;
|
|
* }
|
|
* }
|
|
* \endcode
|
|
*/
|
|
#define jkqtplinalgMatIndex(l,c,C) ((l)*(C)+(c))
|
|
|
|
|
|
/** \brief print the given LxC matrix to std::cout
|
|
* \ingroup jkqtptools_math_linalg
|
|
*
|
|
* \tparam type of the matrix cells (typically double or float)
|
|
* \param matrix the matrix to print out
|
|
* \param L number of lines/rows in the matrix
|
|
* \param C number of columns in the matrix
|
|
* \param width width (in characters) of each cell in the output (used for formatting)
|
|
* \param precision precision (in digits) for string-conversions in the output (used for formatting)
|
|
* \param mode if \c =='f' the mode \c std::fixed is used for output, otherwise \c std::scientific is used
|
|
*/
|
|
template <class T>
|
|
inline void jkqtplinalgPrintMatrix(const T* matrix, long L, long C, int width=9, int precision=3, char mode='f') {
|
|
for (long l=0; l<L; l++) {
|
|
for (long c=0; c<C; c++) {
|
|
if (c>0) std::cout<<", ";
|
|
std::cout.precision(precision);
|
|
std::cout.width(width);
|
|
if (mode=='f') std::cout<<std::fixed<<std::right<<matrix[jkqtplinalgMatIndex(l,c,C)];
|
|
else std::cout<<std::scientific<<std::right<<matrix[jkqtplinalgMatIndex(l,c,C)];
|
|
}
|
|
std::cout<<std::endl;
|
|
}
|
|
}
|
|
|
|
|
|
|
|
/** \brief convert the given LxC matrix to std::string
|
|
* \ingroup jkqtptools_math_linalg
|
|
*
|
|
* \tparam type of the matrix cells (typically double or float)
|
|
* \param matrix the matrix to convert
|
|
* \param L number of lines/rows in the matrix
|
|
* \param C number of columns in the matrix
|
|
* \param width width (in characters) of each cell in the output (used for formatting)
|
|
* \param precision precision (in digits) for string-conversions in the output (used for formatting)
|
|
* \param mode the (printf()) string conversion mode for output of the cell values
|
|
*/
|
|
template <class T>
|
|
inline std::string jkqtplinalgMatrixToString(T* matrix, long L, long C, int width=9, int precision=3, const std::string& mode=std::string("g")) {
|
|
std::string format="%"+jkqtp_inttostr(width)+std::string(".")+jkqtp_inttostr(precision)+std::string("l")+mode;
|
|
std::ostringstream ost;
|
|
for (long l=0; l<L; l++) {
|
|
for (long c=0; c<C; c++) {
|
|
if (c>0) ost<<", ";
|
|
char buf[500];
|
|
sprintf(buf, format.c_str(), jkqtp_todouble(matrix[jkqtplinalgMatIndex(l,c,C)]));
|
|
ost<<buf;
|
|
/*ost.precision(precision);
|
|
ost.width(width);
|
|
if (mode=='f') ost<<std::fixed<<std::right<<matrix[jkqtplinalgMatIndex(l,c,C)];
|
|
else ost<<std::scientific<<std::right<<matrix[jkqtplinalgMatIndex(l,c,C)];*/
|
|
|
|
}
|
|
ost<<std::endl;
|
|
}
|
|
return ost.str();
|
|
}
|
|
|
|
/** \brief maps the number range -1 ... +1 to a color-scale lightblue - white - lightred (used for coloring matrices!)
|
|
* \ingroup jkqtptools_math_linalg
|
|
*
|
|
* \param val the value to convert
|
|
* \param[out] r returns the red value (0..255)
|
|
* \param[out] g returns the green value (0..255)
|
|
* \param[out] b returns the blue value (0..255)
|
|
*/
|
|
inline void jkqtplinalgPM1ToRWBColors(double val, uint8_t& r, uint8_t& g, uint8_t& b){
|
|
r=255;
|
|
g=255;
|
|
b=255;
|
|
const double fval=fabs(val);
|
|
if (val<0 && val>=-1) {
|
|
r=jkqtp_boundedRoundTo<uint8_t>(0,255.0-fval*127.0,255);
|
|
g=jkqtp_boundedRoundTo<uint8_t>(0,255.0-fval*127.0,255);
|
|
} else if (val>0 && val<=1) {
|
|
b=jkqtp_boundedRoundTo<uint8_t>(0,255.0-fval*127.0,255);
|
|
g=jkqtp_boundedRoundTo<uint8_t>(0,255.0-fval*127.0,255);
|
|
} else if (val<-1) {
|
|
r=127;
|
|
g=127;
|
|
b=255;
|
|
} else if (val>1) {
|
|
r=255;
|
|
g=127;
|
|
b=127;
|
|
}
|
|
}
|
|
|
|
/** \brief maps the number range -1 ... +1 to a non-linear color-scale lightblue - white - lightred (used for coloring matrices!)
|
|
* \ingroup jkqtptools_math_linalg
|
|
*
|
|
* \param val the value to convert
|
|
* \param[out] r returns the red value (0..255)
|
|
* \param[out] g returns the green value (0..255)
|
|
* \param[out] b returns the blue value (0..255)
|
|
* \param gamma a gamma-value for non-linear color scaling
|
|
*/
|
|
inline void jkqtplinalgPM1ToNonlinRWBColors(double val, uint8_t& r, uint8_t& g, uint8_t& b, double gamma=0.5){
|
|
if (val<0) {
|
|
jkqtplinalgPM1ToRWBColors(-1.0*pow(-val,gamma),r,g,b);
|
|
} else {
|
|
jkqtplinalgPM1ToRWBColors(pow(val,gamma),r,g,b);
|
|
}
|
|
}
|
|
/** \brief convert the given LxC matrix to std::string, encoded as HTML table
|
|
* \ingroup jkqtptools_math_linalg
|
|
*
|
|
*
|
|
* \tparam type of the matrix cells (typically double or float)
|
|
* \param matrix the matrix to convert
|
|
* \param L number of lines/rows in the matrix
|
|
* \param C number of columns in the matrix
|
|
* \param width width (in characters) of each cell in the output (used for formatting)
|
|
* \param precision precision (in digits) for string-conversions in the output (used for formatting)
|
|
* \param mode the (printf()) string conversion mode for output of the cell values
|
|
* \param tableformat this is inserted into the \c <table...> tag (may contain HTML property definitions)
|
|
* \param prenumber this is inserted before each number (may contain HTML markup)
|
|
* \param postnumber this is inserted after each number (may contain HTML markup)
|
|
* \param colorcoding if \c true, teh cell backgrounds are color-coded
|
|
* \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)
|
|
* \param[out] colorlabel outputs a label explaining the auto-generated color-coding
|
|
* \param nonlinColors if \c true, a non-linear color-coding is used
|
|
* \param nonlinColorsGamma gamma-value for a non-linear color-coding
|
|
* \param colortableformat lie \a tableformat, but for the legend table output in \a colorLabel
|
|
*
|
|
* \see jkqtplinalgPM1ToNonlinRWBColors() and jkqtplinalgPM1ToRWBColors()
|
|
*/
|
|
template <class T>
|
|
inline 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()) {
|
|
std::ostringstream ost;
|
|
ost<<"<table "<<tableformat<<">\n";
|
|
std::string format="%"+jkqtp_inttostr(width)+std::string(".")+jkqtp_inttostr(precision)+std::string("l")+mode;
|
|
double minv=0, maxv=0;
|
|
if (colorcoding) {
|
|
jkqtpstatMinMax(matrix, L*C, minv, maxv);
|
|
}
|
|
for (long l=0; l<L; l++) {
|
|
ost<<" <tr>";
|
|
for (long c=0; c<C; c++) {
|
|
const double val=matrix[jkqtplinalgMatIndex(l,c,C)];
|
|
std::string cols="";
|
|
if (colorcoding) {
|
|
double valrel=0;
|
|
uint8_t r=255,g=255,b=255;
|
|
if (zeroIsWhite){
|
|
if (val<0) valrel=-1.0*fabs(val/minv);
|
|
if (val>0) valrel=fabs(val/maxv);
|
|
} else {
|
|
valrel=((val-minv)/(maxv-minv)-0.5)*2.0;
|
|
}
|
|
if (nonlinColors) {
|
|
jkqtplinalgPM1ToNonlinRWBColors(valrel, r,g,b, nonlinColorsGamma);
|
|
} else {
|
|
jkqtplinalgPM1ToRWBColors(valrel, r,g,b);
|
|
}
|
|
char buf[500];
|
|
sprintf(buf, " bgcolor=\"#%02X%02X%02X\"", int(r), int(g), int(b));
|
|
cols=buf;
|
|
}
|
|
ost<<"<td align=\"center\" valign=\"middle\" width=\""<<(100.0/double(C))<<"%\" "<<cols<<"><nobr>";
|
|
ost.precision(precision);
|
|
ost.width(width);
|
|
char buf[500];
|
|
sprintf(buf, format.c_str(), val);
|
|
ost<<prenumber<<buf<<postnumber;
|
|
ost<<"</nobr></td>";
|
|
}
|
|
ost<<"</tr>\n";
|
|
}
|
|
ost<<"</table>";
|
|
|
|
if (colorcoding && colorlabel) {
|
|
char buf[8192];
|
|
uint8_t rm=255,gm=255,bm=255;
|
|
uint8_t rmc=255,gmc=255,bmc=255;
|
|
uint8_t rc=255,gc=255,bc=255;
|
|
uint8_t rcp=255,gcp=255,bcp=255;
|
|
uint8_t rp=255,gp=255,bp=255;
|
|
double vm=minv;
|
|
double vc=0;
|
|
double vp=maxv;
|
|
if (!zeroIsWhite) {
|
|
vc=(maxv+minv)/2.0;
|
|
}
|
|
if (nonlinColors) {
|
|
jkqtplinalgPM1ToNonlinRWBColors(-1, rm, gm, bm, nonlinColorsGamma);
|
|
jkqtplinalgPM1ToNonlinRWBColors(-0.5, rmc, gmc, bmc, nonlinColorsGamma);
|
|
jkqtplinalgPM1ToNonlinRWBColors(0, rc, gc, bc, nonlinColorsGamma);
|
|
jkqtplinalgPM1ToNonlinRWBColors(0.5, rcp, gcp, bcp, nonlinColorsGamma);
|
|
jkqtplinalgPM1ToNonlinRWBColors(1, rp, gp, bp, nonlinColorsGamma);
|
|
} else {
|
|
jkqtplinalgPM1ToRWBColors(-1, rm, gm, bm);
|
|
jkqtplinalgPM1ToRWBColors(-0.5, rmc, gmc, bmc);
|
|
jkqtplinalgPM1ToRWBColors(0, rc, gc, bc);
|
|
jkqtplinalgPM1ToRWBColors(0.5, rcp, gcp, bcp);
|
|
jkqtplinalgPM1ToRWBColors(1, rp, gp, bp);
|
|
}
|
|
sprintf(buf, "<table %s cellpadding=\"2\" cellspacing=\"0\" border=\"1\"><tr><td><table width=\"100%%\" cellpadding=\"3\" cellspacing=\"0\" border=\"0\"><tr>"
|
|
"<td bgcolor=\"#%02X%02X%02X\" width=\"20%%\"><nobr> %9.3lg </nobr></td>"
|
|
"<td bgcolor=\"#%02X%02X%02X\" width=\"20%%\"><nobr> — </nobr></td>"
|
|
"<td bgcolor=\"#%02X%02X%02X\" width=\"20%%\"><nobr> %9.3lg </nobr></td>"
|
|
"<td bgcolor=\"#%02X%02X%02X\" width=\"20%%\"><nobr> — </nobr></td>"
|
|
"<td bgcolor=\"#%02X%02X%02X\" width=\"20%%\"><nobr> %9.3lg </nobr></td>"
|
|
"</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);
|
|
(*colorlabel)=std::string(buf);
|
|
}
|
|
return ost.str();
|
|
}
|
|
|
|
/** \brief dot-product between two vectors \a vec1 and \a vec2, each with a length of \a N entries
|
|
* \ingroup jkqtptools_math_linalg
|
|
*
|
|
* \tparam T of the vector cells (typically double or float)
|
|
* \param vec1 first vector
|
|
* \param vec2 second vector
|
|
* \param N number of entries in \a vec1 and \a vec2
|
|
*/
|
|
template <class T>
|
|
inline T jkqtplinalgDotProduct(const T* vec1, const T* vec2, long N) {
|
|
T res=0;
|
|
for (long l=0; l<N; l++) {
|
|
res=res+vec1[l]*vec2[l];
|
|
}
|
|
return res;
|
|
}
|
|
|
|
|
|
/** \brief transpose the given NxN matrix
|
|
* \ingroup jkqtptools_math_linalg
|
|
*
|
|
* \tparam T of the matrix cells (typically double or float)
|
|
* \param matrix the matrix to transpose
|
|
* \param N number of rows and columns in the matrix
|
|
*
|
|
*/
|
|
template <class T>
|
|
inline void jkqtplinalgTransposeMatrix(T* matrix, long N) {
|
|
for (long l=0; l<N; l++) {
|
|
for (long c=l+1; c<N; c++) {
|
|
jkqtpArraySwap(matrix, jkqtplinalgMatIndex(l,c,N), jkqtplinalgMatIndex(c,l,N));
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
/** \brief transpose the given LxC matrix
|
|
* \ingroup jkqtptools_math_linalg
|
|
*
|
|
* \tparam T of the matrix cells (typically double or float)
|
|
* \param matrix the matrix to transpose
|
|
* \param L number of rows in the matrix
|
|
* \param C number of columns in the matrix
|
|
*
|
|
* \note The output is interpreted as CxL matrix!!!
|
|
*/
|
|
template <class T>
|
|
inline void jkqtplinalgTransposeMatrix(T* matrix, long L, long C) {
|
|
JKQTPArrayScopedPointer<T> t(jkqtpArrayDuplicate<T>(matrix, L*C));
|
|
for (long l=0; l<L; l++) {
|
|
for (long c=0; c<C; c++) {
|
|
matrix[jkqtplinalgMatIndex(c,l,L)]=t[jkqtplinalgMatIndex(l,c,C)];
|
|
}
|
|
}
|
|
}
|
|
|
|
/** \brief swap two lines in a matrix
|
|
* \ingroup jkqtptools_math_linalg
|
|
*
|
|
* \tparam T of the matrix cells (typically double or float)
|
|
* \param m the matrix to work on
|
|
* \param l1 the row to swap with row \a l2
|
|
* \param l2 the row to swap with row \a l1
|
|
* \param C number of columns in the matrix
|
|
*/
|
|
template <class T>
|
|
inline void jkqtplinalgMatrixSwapLines(T* m, long l1, long l2, long C) {
|
|
for (long c=0; c<C; c++) {
|
|
jkqtpArraySwap(m, jkqtplinalgMatIndex(l1, c, C), jkqtplinalgMatIndex(l2, c, C));
|
|
}
|
|
}
|
|
|
|
|
|
/*! \brief matrix-matrix product
|
|
\ingroup jkqtptools_math_linalg
|
|
|
|
\tparam T of the matrix cells (typically double or float)
|
|
\param M1 matrix 1, size: L1xC1
|
|
\param L1 number of rows in the matrix \a M1
|
|
\param C1 number of columns in the matrix \a M1
|
|
\param M2 matrix 1, size: L2xC2
|
|
\param L2 number of rows in the matrix \a M2
|
|
\param C2 number of columns in the matrix \a M2
|
|
\param[out] M output matrix M=M1*M2, size: L1xC2
|
|
*/
|
|
template <class T>
|
|
inline void jkqtplinalgMatrixProduct(const T* M1, long L1, long C1, const T* M2, long L2, long C2, T* M) {
|
|
if (M1!=M &&M2!=M) {
|
|
|
|
for (long l=0; l<L1; l++) {
|
|
for (long c=0; c<C2; c++) {
|
|
double s=T(0);
|
|
for (long i=0; i<C1; i++) {
|
|
s = s + M1[jkqtplinalgMatIndex(l, i, C1)]*M2[jkqtplinalgMatIndex(i,c, C2)];
|
|
}
|
|
M[jkqtplinalgMatIndex(l, c, C2)]=s;
|
|
}
|
|
}
|
|
} else if (M1==M && M2!=M) {
|
|
JKQTPArrayScopedPointer<T> MM(jkqtpArrayDuplicate(M1, L1*C1));
|
|
jkqtplinalgMatrixProduct(MM.data(),L1,C1,M2,L2,C2,M);
|
|
} else if (M1!=M && M2==M) {
|
|
JKQTPArrayScopedPointer<T> MM(jkqtpArrayDuplicate(M1, L1*C1));
|
|
jkqtplinalgMatrixProduct(M1,L1,C1,MM.data(),L2,C2,M);
|
|
} else if (M1==M && M2==M) {
|
|
JKQTPArrayScopedPointer<T> MM(jkqtpArrayDuplicate(M1, L1*C1));
|
|
jkqtplinalgMatrixProduct(MM.data(),L1,C1,MM.data(),L2,C2,M);
|
|
}
|
|
}
|
|
|
|
|
|
/*! \brief matrix-matrix product of two NxN matrices
|
|
\ingroup jkqtptools_math_linalg
|
|
|
|
\param M1 matrix 1, size: NxN
|
|
\param M2 matrix 1, size: NxN
|
|
\param N number os rows/columns in the matrices \a M1, \a M2 and \a M
|
|
\param[out] M output matrix M=M1*M2, size: NxN
|
|
*/
|
|
template <class T>
|
|
inline void jkqtplinalgMatrixProduct(const T* M1, const T* M2, long N, T* M) {
|
|
jkqtplinalgMatrixProduct(M1,N,N,M2,N,N,M);
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*! \brief performs a Gauss-Jordan eliminaion on a LxC matrix m
|
|
\ingroup jkqtptools_math_linalg
|
|
|
|
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
|
|
\f[ \left[A | I_L\right] \f] where \f$ I_L \f$ is the unit matrix with LxL entries.
|
|
|
|
\tparam T of the matrix cells (typically double or float)
|
|
\param m the matrix
|
|
\param L number of rows in the matrix
|
|
\param C number of columns in the matrix
|
|
|
|
\see http://www.virtual-maxim.de/matrix-invertieren-in-c-plus-plus/
|
|
*/
|
|
template <class T>
|
|
inline bool jkqtplinalgGaussJordan(T* m, long L, long C) {
|
|
|
|
const long N=L;
|
|
|
|
//std::cout<<"\nstep 0:\n";
|
|
//linalgPrintMatrix(m, N, C);
|
|
|
|
// first we perform a Gauss-elimination, which transforms the matrix in the left half of m into upper triangular form
|
|
for (long k=0; k<N-1; k++) {
|
|
//std::cout<<"\nstep G"<<k<<": pivot="<<m[jkqtpstatisticsMatIndex(k,k,C)]<<"\n";
|
|
if (m[jkqtplinalgMatIndex(k,k,C)]==0.0) {
|
|
// if pivot(m[k,k])==0, then swap this line with a line, which does not have a 0 in the k-th column
|
|
for (long ks=k+1; ks<N; ks++) {
|
|
if (m[jkqtplinalgMatIndex(ks,k,C)]!=0.0) {
|
|
jkqtplinalgMatrixSwapLines(m, ks, k, C);
|
|
break;
|
|
} else if (ks==N-1) {
|
|
// if no such element is found, the matrix may not be inverted!
|
|
return false;
|
|
}
|
|
}
|
|
}
|
|
|
|
// now we can eliminate all entries i below the pivot line p, by subtracting
|
|
// the pivot line, scaled by s, from every line, where
|
|
// s=m[i,p]/m[p,p]
|
|
for (long i=k+1; i<N; i++) {
|
|
const T s=m[jkqtplinalgMatIndex(i,k,C)]/m[jkqtplinalgMatIndex(k,k,C)];
|
|
for (long c=k; c<C; c++) {
|
|
m[jkqtplinalgMatIndex(i,c,C)] -= m[jkqtplinalgMatIndex(k,c,C)]*s;
|
|
}
|
|
}
|
|
|
|
//linalgPrintMatrix(m, N, C);
|
|
}
|
|
|
|
// now we can caluate the determinant of the left half-matrix, which can be used to determine, whether matrix
|
|
// is invertible at all: if det(T)==0.0 -> matrix is not invertible
|
|
// the determinant of an upper triangular marix equals the product of the diagonal elements
|
|
T det=1.0;
|
|
for (long k=0; k<N; k++) {
|
|
det = det * m[jkqtplinalgMatIndex(k,k,C)];
|
|
}
|
|
//linalgPrintMatrix(m, N, C);
|
|
//std::cout<<"\nstep 2: det(M)="<<det<<"\n";
|
|
if (fabs(det)<DBL_MIN*10.0) return false;
|
|
|
|
|
|
// if the matrix may be inverted, we can go on with the JOrdan part of the algorithm:
|
|
// we work the Nx(2N) matrix from bottom to top and transform the left side into a unit matrix
|
|
// - the last row is left unchanged
|
|
// - the last row is subtracted from every row i above, scaled by m[i,N]/m[N,N]
|
|
// then we repeat this for the (N-1)*(N-1) left upper matrix, which has again full triangular form
|
|
for (long k=N-1; k>0; k--) {
|
|
//std::cout<<"\nstep J"<<k<<":\n";
|
|
for (long i=k-1; i>=0; i--) {
|
|
const T s=m[jkqtplinalgMatIndex(i,k,C)]/m[jkqtplinalgMatIndex(k,k,C)];
|
|
for (long c=k; c<C; c++) {
|
|
m[jkqtplinalgMatIndex(i,c,C)] -= m[jkqtplinalgMatIndex(k,c,C)]*s;
|
|
}
|
|
}
|
|
//linalgPrintMatrix(m, N, C);
|
|
}
|
|
// finally each line is normalized to 1 by dividing by the diagonal entry in the left NxN matrix
|
|
// and copy the result to matrix_out
|
|
for (long k=0; k<N; k++) {
|
|
|
|
const T s=m[jkqtplinalgMatIndex(k,k,C)];
|
|
for (long c=N; c<C; c++) {
|
|
m[jkqtplinalgMatIndex(k,c,C)] = m[jkqtplinalgMatIndex(k,c,C)] / s;
|
|
}
|
|
m[jkqtplinalgMatIndex(k,k,C)]=T(1);
|
|
}
|
|
|
|
|
|
return true;
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
/*! \brief invert the given NxN matrix using Gauss-Jordan elimination
|
|
\ingroup jkqtptools_math_linalg
|
|
|
|
\tparam T of the matrix cells (typically double or float)
|
|
\param matrix the matrix to invert
|
|
\param[out] matrix_out target for the inverted matrix
|
|
\param N number of rows and columns in the matrix
|
|
|
|
\return \c true on success and the inverted matrix in matrix_out.
|
|
|
|
\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!
|
|
\note matrix and matrix_out have to be of size N*N. Matrices are interpreted as row-major!
|
|
|
|
\see http://www.virtual-maxim.de/matrix-invertieren-in-c-plus-plus/
|
|
*/
|
|
template <class T>
|
|
inline bool jkqtplinalgMatrixInversion(const T* matrix, T* matrix_out, long N) {
|
|
#ifdef JKQTP_STATISTICS_TOOLS_MAY_USE_EIGEN3
|
|
if (N==2) {
|
|
Eigen::Map<const Eigen::Matrix<T,2,2,Eigen::RowMajor> > eA(matrix);
|
|
Eigen::Map<Eigen::Matrix<T,2,2,Eigen::RowMajor> > eO(matrix_out);
|
|
eO=eA.inverse();
|
|
//std::cout<<"\n--------------------------------------\n2x2 input matrix\n"<<eA<<"\n--------------------------------------\n";
|
|
return eO.allFinite();
|
|
} else if (N==3) {
|
|
Eigen::Map<const Eigen::Matrix<T,3,3,Eigen::RowMajor> > eA(matrix);
|
|
Eigen::Map<Eigen::Matrix<T,3,3,Eigen::RowMajor> > eO(matrix_out);
|
|
//std::cout<<"\n--------------------------------------\n3x3 input matrix\n"<<eA<<"\n--------------------------------------\n";
|
|
eO=eA.inverse();
|
|
return eO.allFinite();
|
|
} else if (N==4) {
|
|
Eigen::Map<const Eigen::Matrix<T,4,4,Eigen::RowMajor> > eA(matrix);
|
|
Eigen::Map<Eigen::Matrix<T,4,4,Eigen::RowMajor> > eO(matrix_out);
|
|
//std::cout<<"\n--------------------------------------\n4x4 input matrix\n"<<eA<<"\n--------------------------------------\n";
|
|
eO=eA.inverse();
|
|
return eO.allFinite();
|
|
} else {
|
|
Eigen::Map<const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic,Eigen::RowMajor> > eA(matrix,N,N);
|
|
Eigen::Map<Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic,Eigen::RowMajor> > eO(matrix_out,N,N);
|
|
eO=eA.inverse();
|
|
//std::cout<<"\n--------------------------------------\nNxN input matrix\n"<<eA<<"\n--------------------------------------\n";
|
|
return eO.allFinite();
|
|
}
|
|
return true;
|
|
#else
|
|
// first build a N*(2N) matrix of the form
|
|
//
|
|
// m11 m12 ... | 1 0 0
|
|
// m21 m22 ... | 0 1 0
|
|
// ... ... ... | .....
|
|
//
|
|
const long msize=N*N*2;
|
|
std::vector<T> m;
|
|
m.resize(msize);
|
|
for (long i=0; i<msize; i++) m[i]=T(0); // init with 0
|
|
for (long l=0; l<N; l++) {
|
|
for (long c=0; c<N; c++) { // init left half with matrix
|
|
m[jkqtplinalgMatIndex(l,c,2*N)]=matrix[jkqtplinalgMatIndex(l,c,N)];
|
|
}
|
|
// init right half with unit matrix
|
|
m[jkqtplinalgMatIndex(l,N+l,2*N)]=T(1);
|
|
}
|
|
|
|
|
|
bool ok=jkqtplinalgGaussJordan(m.data(), N, 2*N);
|
|
|
|
if (ok) {
|
|
// finally we copy the result to matrix_out
|
|
for (long k=0; k<N; k++) {
|
|
for (long c=N; c<2*N; c++) {
|
|
matrix_out[jkqtplinalgMatIndex(k,c-N,N)] = m[jkqtplinalgMatIndex(k,c,2*N)];
|
|
}
|
|
}
|
|
}
|
|
|
|
return ok;
|
|
#endif
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
/*! \brief invert the given NxN matrix using Gauss-Jordan elimination
|
|
\ingroup jkqtptools_math_linalg
|
|
|
|
\tparam T of the matrix cells (typically double or float)
|
|
\param[in,out] matrix the matrix to invert (at the same time the target)
|
|
\param N number of rows and columns in the matrix
|
|
*/
|
|
template <class T>
|
|
inline bool jkqtplinalgMatrixInversion(T* matrix, long N) {
|
|
return jkqtplinalgMatrixInversion(matrix, matrix, N);
|
|
}
|
|
|
|
|
|
|
|
/*! \brief solve a system of N linear equations \f$ A\cdot\vec{x}=B \f$ simultaneously for C columns in B
|
|
\ingroup jkqtptools_math_linalg
|
|
|
|
\param A an NxN matrix of coefficients
|
|
\param B an NxC marix
|
|
\param N number of equations
|
|
\param C number of columns in B
|
|
\param result_out a NxC matrix with the results after the inversion of the system of equations
|
|
\return \c true on success
|
|
|
|
\note This function uses the Gauss-Jordan algorithm
|
|
\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!
|
|
*/
|
|
template <class T>
|
|
inline bool jkqtplinalgLinSolve(const T* A, const T* B, long N, long C, T* result_out) {
|
|
#if defined(JKQTP_STATISTICS_TOOLS_MAY_USE_EIGEN3) && (!defined(STATISTICS_TOOLS_linalgLinSolve_EIGENMETHOD_noeigen))
|
|
if (N==2 && C==1) {
|
|
Eigen::Map<const Eigen::Matrix<T,2,2,Eigen::RowMajor> > eA(A);
|
|
Eigen::Matrix<T,2,1> eB;
|
|
Eigen::Map<Eigen::Matrix<T,2,1> > x(result_out);
|
|
|
|
eB=Eigen::Map<const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic,Eigen::RowMajor> >(B,2,1);
|
|
# ifdef STATISTICS_TOOLS_linalgLinSolve_EIGENMETHOD_fullPivLu
|
|
x=eA.fullPivLu().solve(eB);
|
|
# elif defined(STATISTICS_TOOLS_linalgLinSolve_EIGENMETHOD_householderQr)
|
|
x=eA.householderQr().solve(eB);
|
|
# elif defined(STATISTICS_TOOLS_linalgLinSolve_EIGENMETHOD_fullPivHouseholderQr)
|
|
x=eA.fullPivHouseholderQr().solve(eB);
|
|
# elif defined(STATISTICS_TOOLS_linalgLinSolve_EIGENMETHOD_jacobiSvd)
|
|
x=eA.jacobiSVD().solve(eB);
|
|
# elif defined(STATISTICS_TOOLS_linalgLinSolve_EIGENMETHOD_colPivHouseholderQr)
|
|
x=eA.colPivHouseholderQr().solve(eB);
|
|
# else
|
|
x=eA.fullPivLu().solve(eB);
|
|
# endif
|
|
} else if (N==3 && C==1) {
|
|
Eigen::Map<const Eigen::Matrix<T,3,3,Eigen::RowMajor> > eA(A);
|
|
Eigen::Matrix<T,3,1> eB;
|
|
Eigen::Map<Eigen::Matrix<T,3,1> > x(result_out);
|
|
|
|
eB=Eigen::Map<const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic,Eigen::RowMajor> >(B,3,1);
|
|
|
|
# ifdef STATISTICS_TOOLS_linalgLinSolve_EIGENMETHOD_fullPivLu
|
|
x=eA.fullPivLu().solve(eB);
|
|
# elif defined(STATISTICS_TOOLS_linalgLinSolve_EIGENMETHOD_householderQr)
|
|
x=eA.householderQr().solve(eB);
|
|
# elif defined(STATISTICS_TOOLS_linalgLinSolve_EIGENMETHOD_fullPivHouseholderQr)
|
|
x=eA.fullPivHouseholderQr().solve(eB);
|
|
# elif defined(STATISTICS_TOOLS_linalgLinSolve_EIGENMETHOD_jacobiSvd)
|
|
x=eA.jacobiSVD().solve(eB);
|
|
# elif defined(STATISTICS_TOOLS_linalgLinSolve_EIGENMETHOD_colPivHouseholderQr)
|
|
x=eA.colPivHouseholderQr().solve(eB);
|
|
# else
|
|
x=eA.fullPivLu().solve(eB);
|
|
# endif
|
|
} else {
|
|
Eigen::Map<const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic,Eigen::RowMajor> > eA(A,N,N);
|
|
Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic,Eigen::RowMajor> eB(N,C);
|
|
Eigen::Map<Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic,Eigen::RowMajor> > x(result_out,N,C);
|
|
|
|
eB=Eigen::Map<const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic,Eigen::RowMajor> >(B,N,C);
|
|
|
|
# ifdef STATISTICS_TOOLS_linalgLinSolve_EIGENMETHOD_fullPivLu
|
|
x=eA.fullPivLu().solve(eB);
|
|
# elif defined(STATISTICS_TOOLS_linalgLinSolve_EIGENMETHOD_householderQr)
|
|
x=eA.householderQr().solve(eB);
|
|
# elif defined(STATISTICS_TOOLS_linalgLinSolve_EIGENMETHOD_fullPivHouseholderQr)
|
|
x=eA.fullPivHouseholderQr().solve(eB);
|
|
# elif defined(STATISTICS_TOOLS_linalgLinSolve_EIGENMETHOD_jacobiSvd)
|
|
x=eA.jacobiSVD().solve(eB);
|
|
# elif defined(STATISTICS_TOOLS_linalgLinSolve_EIGENMETHOD_colPivHouseholderQr)
|
|
x=eA.colPivHouseholderQr().solve(eB);
|
|
# else
|
|
x=eA.fullPivLu().solve(eB);
|
|
# endif
|
|
}
|
|
return true;
|
|
#else
|
|
// first build a N*(N+C) matrix of the form
|
|
//
|
|
// <---- N ----> <---- C ---->
|
|
// ^ A11 A12 ... | B11 B12 ...
|
|
// | A21 A22 ... | B21 B22 ...
|
|
// N ... ... ... | .....
|
|
// | ... ... ... | .....
|
|
// v ... ... ... | .....
|
|
//
|
|
const long msize=N*(N+C);
|
|
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
|
|
for (long l=0; l<N; l++) {
|
|
for (long c=0; c<N; c++) { // init left half with matrix A
|
|
m[jkqtplinalgMatIndex(l,c,N+C)]=A[jkqtplinalgMatIndex(l,c,N)];
|
|
}
|
|
// init right half with B
|
|
for (long c=0; c<C; c++) { // init left half with matrix B
|
|
m[jkqtplinalgMatIndex(l,N+c,N+C)]=B[jkqtplinalgMatIndex(l,c,C)];
|
|
}
|
|
}
|
|
|
|
|
|
bool ok=jkqtplinalgGaussJordan(m.data(), N, N+C);
|
|
|
|
if (ok) {
|
|
for (long k=0; k<N; k++) {
|
|
for (long c=N; c<(N+C); c++) {
|
|
if (!JKQTPIsOKFloat(m[jkqtplinalgMatIndex(k,c,N+C)])) {
|
|
ok=false;
|
|
break;
|
|
}
|
|
}
|
|
if (!ok) break;
|
|
}
|
|
if (ok) {
|
|
// finally we copy the result to matrix_out
|
|
for (long k=0; k<N; k++) {
|
|
for (long c=N; c<(N+C); c++) {
|
|
result_out[jkqtplinalgMatIndex(k,c-N,C)] = m[jkqtplinalgMatIndex(k,c,N+C)];
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
return ok;
|
|
#endif
|
|
}
|
|
|
|
/*! \brief solve a system of N linear equations \f$ A\cdot\vec{x}=B \f$ simultaneously for C columns in B
|
|
\ingroup jkqtptools_math_linalg
|
|
|
|
\param A an NxN matrix of coefficients
|
|
\param[in,out] B an NxC marix (also receives the result)
|
|
\param N number of equations
|
|
\param C number of columns in B
|
|
\return \c true on success
|
|
|
|
\note This function uses the Gauss-Jordan algorithm
|
|
\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!
|
|
*/
|
|
template <class T>
|
|
inline bool jkqtplinalgLinSolve(const T* A, T* B, long N, long C) {
|
|
return jkqtplinalgLinSolve(A,B,N,C,B);
|
|
}
|
|
|
|
|
|
/*! \brief solve a system of N linear equations \f$ A\cdot\vec{x}=\vec{b} \f$
|
|
\ingroup jkqtptools_math_linalg
|
|
|
|
\param A an NxN matrix of coefficients
|
|
\param b an N-entry vector
|
|
\param N number of rows and columns in \a A
|
|
\param[out] result_out a N-entry vector with the result
|
|
\return \c true on success
|
|
|
|
\note This function uses the Gauss-Jordan algorithm
|
|
\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!
|
|
*/
|
|
template <class T>
|
|
inline bool jkqtplinalgLinSolve(const T* A, const T* b, long N, T* result_out) {
|
|
return jkqtplinalgLinSolve(A, b, N, 1, result_out);
|
|
}
|
|
|
|
/*! \brief solve a system of N linear equations \f$ A\cdot\vec{x}=\vec{b} \f$
|
|
\ingroup jkqtptools_math_linalg
|
|
|
|
\param A an NxN matrix of coefficients
|
|
\param[in,out] b an N-entry vector (also receives the result)
|
|
\param N number of rows and columns in \a A
|
|
\return \c true on success
|
|
|
|
\note This function uses the Gauss-Jordan algorithm
|
|
\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!
|
|
*/
|
|
template <class T>
|
|
inline bool jkqtplinalgLinSolve(const T* A, T* b, long N) {
|
|
return jkqtplinalgLinSolve(A,b,N,1,b);
|
|
}
|
|
|
|
|
|
|
|
|
|
/*! \brief determinant the given NxN matrix
|
|
\ingroup jkqtptools_math_linalg
|
|
|
|
\tparam T of the matrix cells (typically double or float)
|
|
\param a the matrix for which to calculate the determinant
|
|
\param N number of rows and columns in the matrix
|
|
\return the determinant of \a a
|
|
|
|
\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!
|
|
*/
|
|
template <class T>
|
|
inline T jkqtplinalgMatrixDeterminant(const T* a, long N) {
|
|
#ifdef JKQTP_STATISTICS_TOOLS_MAY_USE_EIGEN3
|
|
if (N < 1) { /* Error */
|
|
return NAN;
|
|
} else if (N == 1) { /* Shouldn't get used */
|
|
return a[jkqtplinalgMatIndex(0,0,N)];
|
|
} else if (N == 2) {
|
|
return a[jkqtplinalgMatIndex(0,0,N)] * a[jkqtplinalgMatIndex(1,1,N)] - a[jkqtplinalgMatIndex(1,0,N)] * a[jkqtplinalgMatIndex(0,1,N)];
|
|
} else if (N==3){
|
|
Eigen::Map<const Eigen::Matrix<T,3,3,Eigen::RowMajor> > eA(a);
|
|
return eA.determinant();
|
|
} else if (N==4){
|
|
Eigen::Map<const Eigen::Matrix<T,4,4,Eigen::RowMajor> > eA(a);
|
|
return eA.determinant();
|
|
} else {
|
|
Eigen::Map<const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic,Eigen::RowMajor> > eA(a,N,N);
|
|
return eA.determinant();
|
|
}
|
|
#else
|
|
long i,j,j1,j2;
|
|
T det = 0;
|
|
|
|
if (N < 1) { /* Error */
|
|
det=NAN;
|
|
} else if (N == 1) { /* Shouldn't get used */
|
|
det = a[jkqtplinalgMatIndex(0,0,N)];
|
|
} else if (N == 2) {
|
|
det = a[jkqtplinalgMatIndex(0,0,N)] * a[jkqtplinalgMatIndex(1,1,N)] - a[jkqtplinalgMatIndex(1,0,N)] * a[jkqtplinalgMatIndex(0,1,N)];
|
|
} else {
|
|
det = 0;
|
|
for (j1=0;j1<N;j1++) {
|
|
JKQTPArrayScopedPointer<T> m(static_cast<T*>(calloc((N-1)*(N-1),sizeof(T *))));
|
|
|
|
for (i=1;i<N;i++) {
|
|
j2 = 0;
|
|
for (j=0;j<N;j++) {
|
|
if (j != j1){
|
|
m[jkqtplinalgMatIndex(i-1,j2,N-1)] = a[jkqtplinalgMatIndex(i,j,N)];
|
|
j2++;
|
|
}
|
|
}
|
|
}
|
|
//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));
|
|
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);
|
|
|
|
}
|
|
}
|
|
|
|
//printf("\n\n jkqtplinalgMatrixDeterminant(%d):\n",N);
|
|
//linalgPrintMatrix(a,N,N);
|
|
//printf(" jkqtplinalgMatrixDeterminant(%d) = %lf\n", N, det);
|
|
return(det);
|
|
#endif
|
|
}
|
|
|
|
|
|
|
|
|
|
#endif // JKQTPLINALGTOOLS_H_INCLUDED
|
|
|
|
|