#include <Imagine/Common.h>
#include <Imagine/LinAlg.h>
#include <cmath>
#include <iostream>
#include <vector>
#include <utility>
using namespace std;
template <typename T>
void vectors() {
cout << "Vector functions" << endl;
T t[]={1,2,3};
T *t2=new T[3];
c=b;
c=a+b;
c+=b;
c=a-b;
c-=b;
c=a+1;
c+=1;
c=a-1;
c-=1;
c=-a;
c=a*2;
c*=2;
c=a/2;
c/=2;
c=1+a;
c=1-a;
c=2*a;
T x;
x=a*b;
x=norm2(a);
x=norm(a);
(void)x;
c=normalized(a);
x=maxNorm(a);
x=sum(a);
x=mean(a);
c=convolution(K,a);
c=truncConvolution(K,a,1);
a[1]=2;b[2]=2;
x=correlation(a,b);
x=normalizedCorrelation(a,b);
dist(a,b);
}
template <typename T>
void matrices() {
cout << "Matrix functions" << endl;
for(int i=0; i<3; i++)
for(int j=0; j<4; j++)
B(i,j)=rand()/(T)RAND_MAX;
T tt[]={1,2,3,4,5,6};
T *tt2=new T[6];
for(int i=0; i<6; i++)
tt2[i]=rand()/(T)RAND_MAX;
A=B;
A(1,2)=3;
w=A*v;
A(1,2)=2;
A34=A32*A24;
C=A+B;
C+=B;
C=A-B;
C-=B;
C=-A;
C=A*2;
C*=2;
C=A/2;
C/=2;
C=2*A;
T a;
a=dot(A,C);
a=norm(A);
a=norm(A,'1');
a=norm(A,'I');
a=norm(A,'M');
a=rcond(A);
(void)a;
cout << "SVD check 1: "
cout << "SVD check 2: "
cout << "SVD check 3: "
T coeff[] = {1,2,2,1};
cout << "Eigen values of matrix E=" << E;
cout << "Real parts: " << lambdaR << " (should be 3, -1)" << endl;
cout << "Imaginary parts: " << lambdaI << " (should be 0)" << endl;
for (int i=0;i<5;i++)
for (int j=0;j<5;j++)
for (int j=0;j<5;j++)
cout << "Square system: " << norm(A*x-b) << endl;
cout << "Over determined system: " << norm(C*x-b) << endl;
cout << "Minimum norm under determined system: "
<< norm(C*x-b.getSubVect(0,4)) << " " << norm(x) << endl;
cout << "Pseudo inverse, square: " << norm(A*x-b) << endl;
cout << "Pseudo inverse, over determined: " << norm(C*x-b) << endl;
cout << "Pseudo inverse, under determined: "
<< norm(C*x-b.getSubVect(0,4)) << " " << norm(x) << endl;
cout << "QR " << norm(D-Q*R) << endl;
cout << "QR " << norm(D-Q*R) << endl;
cout <<
"Lower Cholesky " << norm(B-L*
transpose(L)) << endl;
cout <<
"Upper Cholesky " << norm(B-
transpose(U)*U) << endl;
cout << "A=" << A << endl;
cout <<
"Det: " <<
det(A) << endl;
ofstream out("tmp.txt");
out << A << endl;
out.close();
ifstream in("tmp.txt");
in >> A;
in.close();
}
template <typename T>
void symMatrices() {
cout << "SymMatrix functions" << endl;
T t[]={1,2,3,4,5};
T *t2=new T[5];
F=B;
F=B.clone();
F.setSize(5);
F.nrow(),F.ncol();
F.fill(2);
T a=F(2,2);
F(2,3)=a;
F.
fill(1);F(0,0)=2;F(0,1)=3;F(1,2)=2;F(2,3)=4;F(3,4)=-1;
cout <<
"Determinant: " <<
det(F)*
det(invF)-1 << endl;
EVD(F,Q,Lambda);
cout << "Check EVD: "
F=F*F;
invF=posDefInverse(F);
P=F*G;
P=P*F;
P=F*P;
F=F+G;
F+=G;
F=F-G;
F-=G;
F=F*2;
F*=2;
F=2*F;
F=F/2;
F/=2;
F=-F;
cout << norm(F*x-b) << endl;
ofstream out("tmp.bin",ios::binary);
F.write(out);
out.close();
ifstream in("tmp.bin",ios::binary);
F.read(in);
in.close();
out.open("tmp.txt");
out << F << endl;
out.close();
in.open("tmp.txt");
in >> F;
in.close();
}
template <typename T>
void fmatrices() {
cout << "FMatrix functions" << endl;
A(1,1)=A(1,2)=A(2,1)=T(3);
A2(1,2)=2;
A3(1,2)=12;A3(1,1)=0;A3(0,0) = 4;
cout << "A3=" << A3 << endl;
cout << "A3*A4=" << A3*A4 << endl;
cout << "SVD check 1: "
for(int i=0; i < 3; i++)
S2(i,i) = S[i];
cout << "SVD check 2: "
<< norm(A-U*S2*Vt2) << endl;
cout << "SVD check 3: "
for (int i=0;i<5;i++)
for (int j=0;j<5;j++)
for (int j=0;j<5;j++)
cout << "Square system: " << norm(B*x-b) << endl;
cout << "Over determined system: "
<< norm(C*x2-b) << endl;
cout << "Minimum norm under determined system: "
<< norm(D*x3-d) << " " << norm(x3) << endl;
cout << "Pseudo inverse, square: "
<< norm(B*x-b) << endl;
cout << "Pseudo inverse, over determined: "
<< norm(C*x2-b) << endl;
cout << "Pseudo inverse, under determined: "
<< norm(D*x3-d) << " " << norm(x3) << endl;
T tt[]={1,2,3,4,5,6};
cout << "QR " << norm(D2-Q*R) << endl;
cout << "QR " << norm(D2-Q2*R2) << endl;
cout << "Lower Cholesky "
cout << "Upper Cholesky "
cout << "Determinants: "
cout << "Determinants: "
}
int main(int argc, char* argv[]) {
typedef void (*FUNC)();
vector< pair<string,FUNC> > tests;
#define ADD_TEST(func) \
tests.push_back(make_pair(string(#func), func));
ADD_TEST(vectors<float>);
ADD_TEST(vectors<double>);
ADD_TEST(matrices<float>);
ADD_TEST(matrices<double>);
ADD_TEST(symMatrices<float>);
ADD_TEST(symMatrices<double>);
ADD_TEST(fmatrices<float>);
ADD_TEST(fmatrices<double>);
std::vector< pair<string,FUNC> >::const_iterator it;
if(argc<=1)
for(it=tests.begin(); it!=tests.end(); ++it) {
cout << "---" << it->first << "---" << std::endl;
it->second();
}
else {
bool fail=false;
for(int i=1; i<argc; i++) {
for(it=tests.begin(); it!=tests.end(); ++it) {
if(it->first == argv[i]) {
cout << "---" << it->first << "---" << std::endl;
it->second();
break;
}
}
if(it==tests.end()) {
std::cout << "Unknown test " << argv[i] << endl;
fail = true;
}
}
if(fail) {
cout << "--- Available tests:" << endl;
for(it=tests.begin(); it!=tests.end(); ++it)
cout << it->first << ' ';
cout << endl;
return 0;
}
}
return 0;
}
const T * data() const
Data pointer (read).
Definition: FArray.h:188
Matrix of fixed dimension.
Definition: FMatrix.h:17
const T * data() const
Data pointer (read).
Definition: FMatrix.h:182
Matrix of variable size.
Definition: Matrix.h:23
Vector< T > getCol(int j) const
Get column.
Definition: Matrix.h:149
Vector< T > getRow(int i) const
Get row.
Definition: Matrix.h:162
int nrow() const
Number of rows.
Definition: Matrix.h:101
Matrix getSubMat(int i0, int m, int j0, int n) const
Get sub matrix.
Definition: Matrix.h:256
int ncol() const
Number of columns.
Definition: Matrix.h:108
Vector< T > getColRef(int j)
Get column by reference.
Definition: Matrix.h:229
Matrix clone() const
Cloning.
Definition: Matrix.h:89
Matrix & setDiagonal(const Vector< T > &v)
Set diagonal.
Definition: Matrix.h:214
Matrix & fill(T x)
Filling.
Definition: Matrix.h:116
Vector< T > getDiagonal() const
Get diagonal.
Definition: Matrix.h:174
Vector< T > getSubColRef(int j, int i0, int m)
Get sub column by reference.
Definition: Matrix.h:243
Matrix & setCol(int j, const Vector< T > &v)
Set column.
Definition: Matrix.h:189
Matrix getColsRef(int j0, int n)
Get part of columns by reference.
Definition: Matrix.h:272
Matrix & setRow(int i, const Vector< T > &v)
Set row.
Definition: Matrix.h:202
void setSize(const Coords< dim > &sz)
Change sizes.
Definition: MultiArray.h:146
Symmetric Matrix of variable size.
Definition: SymMatrix.h:22
Vector of variable size.
Definition: Vector.h:24
Vector & fill(T x)
Filling.
Definition: Vector.h:93
Vector clone() const
Cloning.
Definition: Vector.h:80
Vector getSubVect(size_t offset, size_t size) const
Sub vector.
Definition: Vector.h:114
Vector getSubVectRef(size_t offset, size_t size)
Sub vector by reference.
Definition: Vector.h:126
Vector & normalize()
Euclidean in-place normalization.
Definition: Vector.h:340
FMatrix< T, M, O > multt(const FMatrix< T, M, N > &A, const FMatrix< T, O, N > &B)
Product.
Definition: FMatrix.h:630
void initRandom(unsigned int s)
Init with seed.
Definition: Random.h:18
FMatrix< T, M, M > Diagonal(const FVector< T, M > &d)
Diagonal.
Definition: FMatrix.h:743
T det(const FMatrix< T, M, M > &A)
Determinant.
Definition: FMatrix.h:717
FMatrix< T, M, M > inverse(const FMatrix< T, M, M > &A)
Inverse.
Definition: FMatrix.h:678
FMatrix< T, N, M > transpose(const FMatrix< T, M, N > &A)
Transpose.
Definition: FMatrix.h:564
double doubleRandom()
Uniform double.
Definition: Random.h:38
void waitKey(const char *message="Press <enter> to continue...")
Pause in program execution until key press.
Definition: IO.h:13
FMatrix< T, N, O > tmult(const FMatrix< T, M, N > &A, const FMatrix< T, M, O > &B)
Product.
Definition: FMatrix.h:615
FMatrix< T, N, O > tmultt(const FMatrix< T, M, N > &A, const FMatrix< T, O, M > &B)
Product.
Definition: FMatrix.h:645
T detFMatrix(const FMatrix< T, N, N > &A)
Determinant.
Definition: Matrix.h:1099
FMatrix< T, N, M > pseudoInverse(const FMatrix< T, M, N > &A, T tolrel=0)
Pseudo inverse.
Definition: Matrix.h:995
void eigenvalues(const Matrix< T > &A, Vector< T > &wr, Vector< T > &wi)
Eigenvalues calculation Given square matrix, only eigenvalues are obtained.
Definition: Matrix.h:907
FVector< T, N > linSolve(const FMatrix< T, M, N > &A, const FVector< T, M > &b)
Linear system.
Definition: Matrix.h:935
FMatrix< T, N, N > cholesky(const FMatrix< T, N, N > &A, bool low=true)
Cholesky decomposition.
Definition: Matrix.h:1067
bool QRAll(const FMatrix< T, M, N > &A, FMatrix< T, M, M > &Q, FMatrix< T, M, N > &R)
QR decomposition.
Definition: Matrix.h:1044
bool QR(const FMatrix< T, M, N > &A, FMatrix< T, M, N > &Q, FMatrix< T, N, N > &R)
QR decomposition.
Definition: Matrix.h:1022
FMatrix< T, N, N > inverseFMatrix(const FMatrix< T, N, N > &A)
Inverse.
Definition: Matrix.h:853
void svd(const FMatrix< T, M, N > &A, FMatrix< T, M, M > &U, FVector< T, M > &S, FMatrix< T, N, N > &Vt, bool all=false)
SVD.
Definition: Matrix.h:886
Matrix< T > outerProduct(const Vector< T > &v1, const Vector< T > &v2)
Outer product matrix.
Definition: Matrix.h:830