Overview
Examples
Screenshots
Comparisons
Applications
Download
Documentation
Tutorials
Bazaar
Status & Roadmap
FAQ
Authors & License
Forums
Funding Ultimate++
Search on this site











SourceForge.net Logo

Eigen demo

 

Eigen package demo

 

 

eigen_demo.cpp

 

#include <Core/Core.h>

#include <plugin/Eigen/Eigen.h>

 

using namespace Upp;

 

using namespace Eigen;

 

void NonLinearTests();

void FFTTests();

 

struct SerialTest {

    MatrixXd m;

    VectorXd v;

    SerialTest() : m(2, 2), v(3){}

    void Print() {

        Cout() << "\nHere is the matrix:\n" << m;

        Cout() << "\nHere is the vector:\n" << v;

    }

    void Serialize(Stream& stream) {

        ::Serialize(stream, m);

        ::Serialize(stream, v);

    }

    void Jsonize(JsonIO &jio) {

        jio("matrix", m)("vector", v);

    }

    void Xmlize(XmlIO &xml) {

        xml("matrix", m)("vector", v);

    }

};

        

CONSOLE_APP_MAIN

{

    StdLogSetup(LOG_COUT|LOG_FILE);

    

    UppLog() << "Eigen library demo";

    

    // https://eigen.tuxfamily.org/dox/group__TutorialMatrixClass.html

    UppLog() << "\n\nTutorial page 1 - The Matrix class";

    

    UppLog() << "\n\nCoefficient accessors";

    {

        MatrixXd m(2,2);

        m(0,0) = 3;

        m(1,0) = 2.5;

        m(0,1) = -1;

        m(1,1) = m(1,0) + m(0,1);

        UppLog() << "\nHere is the matrix m:\n" << m;

        

        VectorXd v(2);

        v(0) = 4;

        v(1) = v(0) - 1;

        UppLog() << "\nHere is the vector v:\n" << v;

    }

    UppLog() << "\n\nResizing";

    {

        MatrixXd m(2,5);

        m.resize(4,3);

        UppLog() << "\nThe matrix m is of size " << m.rows() << "x" << m.cols();

        UppLog() << "\nIt has " << m.size() << " coefficients";

 

        VectorXd v(2);

        v.resize(5);

        UppLog() << "\nThe vector v is of size " << v.size();

        UppLog() << "\nAs a matrix, v is of size " << v.rows() << "x" << v.cols();

    }

    UppLog() << "\n\nAssignment and resizing";

    {

        double _dat[] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17};// Assignment from C vector

        VectorXd dat = Map<VectorXd>(_dat, sizeof(_dat)/sizeof(double));  

        UppLog() << "\nC array data is " << dat.transpose();

        

        const int dec = 5;

       VectorXd decimated = Map<VectorXd, 0, InnerStride<dec>>(dat.data(), 1+((dat.size()-1)/dec));

       UppLog() << "\nDecimated " << decimated.transpose();

        

        VectorXd even = Map<VectorXd, 0, InnerStride<2>>(dat.data(), (dat.size()+1)/2);

        UppLog() << "\nEven " << even.transpose();

        

        VectorXd odd = Map<VectorXd, 0, InnerStride<2>>(dat.data()+1, dat.size()/2);

        UppLog() << "\nOdd  " << odd.transpose();

    

        MatrixXf a(2,2);

        UppLog() << "\na is of size " << a.rows() << "x" << a.cols();

        MatrixXf b(3,3);

        a = b;

        UppLog() << "\na is now of size " << a.rows() << "x" << a.cols();

    }

    

    // https://eigen.tuxfamily.org/dox-devel/group__TutorialMatrixArithmetic.html

    UppLog() << "\n\nTutorial page 2 - Matrix and vector arithmetic";

    

    UppLog() << "\n\nAddition and subtraction";

    {

        Matrix2d a;

        a << 1, 2,

             3, 4;

        MatrixXd b(2,2);

        b << 2, 3,

             1, 4;

 

        UppLog() << "\na + b =\n" << a + b;

        UppLog() << "\na - b =\n" << a - b;

        UppLog() << "\nDoing a += b;";

        a += b;

        UppLog() << "\nNow a =\n" << a;

 

        Vector3d v(1,2,3);

        Vector3d w(1,0,0);

        UppLog() << "\n-v + w - v =\n" << -v + w - v;

    }

    UppLog() << "\n\nScalar multiplication and division";

    {

        Matrix2d a;

        a << 1, 2,

             3, 4;

        Vector3d v(1,2,3);

        UppLog() << "\na * 2.5 =\n" << a * 2.5;

        UppLog() << "\n0.1 * v =\n" << 0.1 * v;

        UppLog() << "\nDoing v *= 2;";

        v *= 2;

        UppLog() << "\nNow v =\n" << v;

    }

    UppLog() << "\n\nTransposition and conjugation";

    {

        MatrixXcf a = MatrixXcf::Random(2,2);

        UppLog() << "\nHere is the matrix a\n" << a;

        UppLog() << "\nHere is the matrix a^T\n" << a.transpose();

        UppLog() << "\nHere is the conjugate of a\n" << a.conjugate();

        UppLog() << "\nHere is the matrix a^*\n" << a.adjoint();

        

        VectorXd v(5);

        v << 1, 2, 3, 4, 5;

        UppLog() << "\n\nInitial vector  " << v.transpose();

        UppLog() <<   "\nReversed vector " << v.reverse().transpose();

    }

    UppLog() << "\n\nMatrix-matrix and matrix-vector multiplication";

    {

        Matrix2d mat;

        mat << 1, 2,

               3, 4;

        Vector2d u(-1,1), v(2,0);

        UppLog() << "\nHere is mat*mat:\n" << mat*mat;

        UppLog() << "\nHere is mat*u:\n" << mat*u;

        UppLog() << "\nHere is u^T*mat:\n" << u.transpose()*mat;

        UppLog() << "\nHere is u^T*v:\n" << u.transpose()*v;

        UppLog() << "\nHere is u*v^T:\n" << u*v.transpose();

        UppLog() << "\nLet's multiply mat by itself";

        mat *= mat;

        UppLog() << "\nNow mat is mat:\n" << mat;

    }

    UppLog() << "\n\nDot product and cross product";

    {

        Vector3d v(1,2,3);

          Vector3d w(0,1,2);

    

          UppLog() << "\nDot product: " << v.dot(w);

          double dp = v.adjoint()*w; // automatic conversion of the inner product to a scalar

          UppLog() << "\nDot product via a matrix product: " << dp;

          UppLog() << "\nCross product:\n" << v.cross(w);

    }

    UppLog() << "\n\nBasic arithmetic reduction operations";

    {

        Eigen::Matrix2d mat;

        mat << 1, 2,

               3, 4;

        UppLog() << "\nHere is mat.sum():       " << mat.sum();

        UppLog() << "\nHere is mat.prod():      " << mat.prod();

        UppLog() << "\nHere is mat.mean():      " << mat.mean();

        UppLog() << "\nHere is mat.minCoeff():  " << mat.minCoeff();

        UppLog() << "\nHere is mat.maxCoeff():  " << mat.maxCoeff();

        UppLog() << "\nHere is mat.trace():     " << mat.trace();

        

        Matrix3f m = Matrix3f::Random();

        ptrdiff_t i, j;

        float minOfM = m.minCoeff(&i, &j);

        UppLog() << "\nHere is the matrix m:\n" << m;

        UppLog() << "\nIts minimum coefficient (" << minOfM

               << ") is at position (" << i << "," << j << ")\n";

        

        RowVector4i v = RowVector4i::Random();

        ptrdiff_t maxOfV = v.maxCoeff(&i);

        UppLog() << "\nHere is the vector v: " << v;

        UppLog() << "\nIts maximum coefficient (" << maxOfV

               << ") is at position " << i;

    }

    

    // https://eigen.tuxfamily.org/dox/group__TutorialArrayClass.html

    UppLog() << "\n\nTutorial page 3 - The Array class and coefficient-wise operations ";

    

    UppLog() << "\n\nAccessing values inside an Array";

    {

        ArrayXXf  m(2,2);

        

        // assign some values coefficient by coefficient

        m(0,0) = 1.0; m(0,1) = 2.0;

        m(1,0) = 3.0; m(1,1) = m(0,1) + m(1,0);

        

        // print values to standard output

        UppLog() << "\n" << m;

        

        // using the comma-initializer is also allowed

        m << 1.0,2.0,

                3.0,4.0;

        

        // print values to standard output

        UppLog() << "\n" << m;

    }

    UppLog() << "\n\nAddition and subtraction";

    {

        ArrayXXf a(3,3);

        ArrayXXf b(3,3);

        a << 1,2,3,

                4,5,6,

                7,8,9;

        b << 1,2,3,

                1,2,3,

                1,2,3;

           

        // Adding two arrays

        UppLog() << "\na + b = " << "\n" << a + b;

        

        // Subtracting a scalar from an array

        UppLog() << "\na - 2 = " << "\n" << a - 2;

    }

    UppLog() << "\n\nArray multiplication";

    {

        ArrayXXf a(2,2);

        ArrayXXf b(2,2);

        a << 1,2,

                3,4;

        b << 5,6,

                7,8;

        UppLog() << "\na * b = " << "\n" << a * b;

    }

    UppLog() << "\n\nOther coefficient-wise operations";

    {

        ArrayXf a = ArrayXf::Random(5);

        a *= 2;

        UppLog() << "\na =" << "\n" << a;

        UppLog() << "\na.abs() =" << "\n" << a.abs();

        UppLog() << "\na.abs().sqrt() =" << "\n" << a.abs().sqrt();

        UppLog() << "\na.min(a.abs().sqrt()) =" << "\n" << a.min(a.abs().sqrt());

    }

    UppLog() << "\n\nConverting between array and matrix expressions";

    {

        MatrixXf m(2,2);

        MatrixXf n(2,2);

        m << 1,2,

                3,4;

        n << 5,6,

                7,8;

        

        UppLog() << "\n-- Matrix m*n: --" << "\n" << m * n;

        UppLog() << "\n-- Array m*n: --" << "\n" << m.array() * n.array();

        UppLog() << "\n-- With cwiseProduct: --" << "\n" << m.cwiseProduct(n);

        UppLog() << "\n-- Array m + 4: --" << "\n" << m.array() + 4;

    }

    {

        MatrixXf m(2,2);

        MatrixXf n(2,2);

        m << 1,2,

             3,4;

        n << 5,6,

             7,8;

        

        UppLog() << "\n-- Combination 1: --" << "\n" << (m.array() + 4).matrix() * m;

        UppLog() << "\n-- Combination 2: --" << "\n" << (m.array() * n.array()).matrix() * m;

    }

    

    // https://eigen.tuxfamily.org/dox/group__TutorialBlockOperations.html

    UppLog() << "\n\nTutorial page 4 - Block operations";

    

    UppLog() << "\n\nUsing block operations";

    {

        Eigen::MatrixXf m(4,4);

        m << 1, 2, 3, 4,

             5, 6, 7, 8,

             9,10,11,12,

                13,14,15,16;

        UppLog() << "\nBlock in the middle\n";

        UppLog() << m.block<2,2>(1,1);

        for (ptrdiff_t i = 1; i <= 3; ++i) {

            UppLog() << "\nBlock of size " << i << "x" << i << "\n";

            UppLog() << m.block(0, 0, i, i);

        }

    }

    {

        Array22d m;

        m << 1,2,

                3,4;

        Array44d a = Array44d::Constant(0.6);

        UppLog() << "\nHere is the array a:\n" << a;

        a.block<2,2>(1,1) = m;

        UppLog() << "\nHere is now a with m copied into its central 2x2 block:\n" << a;

        a.block(0,0,2,3) = a.block(2,1,2,3);

        UppLog() << "\nHere is now a with bottom-right 2x3 block copied into top-left 2x2 block:\n" << a;

    }

    UppLog() << "\n\nColumns and rows";

    {

        Eigen::MatrixXf m(3,3);

        m << 1,2,3,

             4,5,6,

                7,8,9;

        UppLog() << "\nHere is the matrix m:\n" << m;

        UppLog() << "\n2nd Row: " << m.row(1);

        m.col(2) += 3 * m.col(0);

        UppLog() << "\nAfter adding 3 times the first column into the third column, the matrix m is:\n";

        UppLog() << m;

    }

    UppLog() << "\n\nCorner-related operations";

    {

        Eigen::Matrix4f m;

        m << 1, 2, 3, 4,

                5, 6, 7, 8,

                9, 10,11,12,

                13,14,15,16;

        UppLog() << "\nm.leftCols(2) =\n" << m.leftCols(2);

        UppLog() << "\nm.bottomRows<2>() =\n" << m.bottomRows<2>();

        m.topLeftCorner(1,3) = m.bottomRightCorner(3,1).transpose();

        UppLog() << "\nAfter assignment, m = \n" << m;

    }

    UppLog() << "\n\nBlock operations for vectors";

    {

        Eigen::ArrayXf v(6);

        v << 1, 2, 3, 4, 5, 6;

        UppLog() << "\nv.head(3) =\n" << v.head(3);

        UppLog() << "\nv.tail<3>() = \n" << v.tail<3>();

        v.segment(1,4) *= 2;

        UppLog() << "\nafter 'v.segment(1,4) *= 2', v =\n" << v;

    }

    

    // https://eigen.tuxfamily.org/dox/group__TutorialAdvancedInitialization.html

    UppLog() << "\n\nTutorial page 5 - Advanced initialization";

 

    UppLog() << "\n\nThe comma initializer";

    {

        RowVectorXd vec1(3);

        vec1 << 1, 2, 3;

        UppLog() << "\nvec1 = " << vec1;

        

        RowVectorXd vec2(4);

        vec2 << 1, 4, 9, 16;

        UppLog() << "\nvec2 = " << vec2;

        

        RowVectorXd joined(7);

        joined << vec1, vec2;

        UppLog() << "\njoined = " << joined;

        

        MatrixXf matA(2, 2);

        matA << 1, 2, 3, 4;

        MatrixXf matB(4, 4);

        matB << matA, matA/10, matA/10, matA;

        UppLog() << matB;

        

        Matrix3f m;

        m.row(0) << 1, 2, 3;

        m.block(1,0,2,2) << 4, 5, 7, 8;

        m.col(2).tail(2) << 6, 9;                   

        UppLog() << m;

    }

    UppLog() << "\n\nSpecial matrices and arrays";

    {

        UppLog() << "\nA fixed-size array:\n";

        Array33f a1 = Array33f::Zero();

        UppLog() << a1 << "\n\n";

        

        UppLog() << "\nA one-dimensional dynamic-size array:\n";

        ArrayXf a2 = ArrayXf::Zero(3);

        UppLog() << a2 << "\n\n";

        

        UppLog() << "\nA two-dimensional dynamic-size array:\n";

        ArrayXXf a3 = ArrayXXf::Zero(3, 4);

        UppLog() << a3 << "\n";

        

        UppLog() << "\nA two-dimensional dynamic-size array set to 1.23:\n";

        MatrixXd a4 = MatrixXd::Constant(3, 4, 1.23);

        UppLog() << a4 << "\n";

        

        ArrayXXd table(10, 4);

        table.col(0) = ArrayXd::LinSpaced(10, 0, 90);

        table.col(1) = M_PI / 180 * table.col(0);

        table.col(2) = table.col(1).sin();

        table.col(3) = table.col(1).cos();

        UppLog() << "\n  Degrees   Radians      Sine    Cosine\n";

        UppLog() << table;

        

        const ptrdiff_t size = 6;

        MatrixXd mat1(size, size);

        mat1.topLeftCorner(size/2, size/2)     = MatrixXd::Zero(size/2, size/2);

        mat1.topRightCorner(size/2, size/2)    = MatrixXd::Identity(size/2, size/2);

        mat1.bottomLeftCorner(size/2, size/2)  = MatrixXd::Identity(size/2, size/2);

        mat1.bottomRightCorner(size/2, size/2) = MatrixXd::Zero(size/2, size/2);

        UppLog() << "\n" << mat1;

        

        MatrixXd mat2(size, size);

        mat2.topLeftCorner(size/2, size/2).setZero();

        mat2.topRightCorner(size/2, size/2).setIdentity();

        mat2.bottomLeftCorner(size/2, size/2).setIdentity();

        mat2.bottomRightCorner(size/2, size/2).setZero();

        UppLog() << "\n" << mat2;

        

        MatrixXd mat3(size, size);

        mat3 << MatrixXd::Zero(size/2, size/2), MatrixXd::Identity(size/2, size/2),

                MatrixXd::Identity(size/2, size/2), MatrixXd::Zero(size/2, size/2);

        UppLog() << "\n" << mat3;

    }

    UppLog() << "\n\nUsage as temporary objects";

    {

        MatrixXd m = MatrixXd::Random(3,3);

        m = (m + MatrixXd::Constant(3,3,1.2)) * 50;

        UppLog() << "\nm =\n" << m;

        VectorXd v(3);

        v << 1, 2, 3;

        UppLog() << "\nm * v =\n" << m * v;

    }

    {

        MatrixXf mat = MatrixXf::Random(2, 3);

        UppLog() << mat;

        mat = (MatrixXf(2,2) << 0, 1, 1, 0).finished() * mat;

        UppLog() << mat;

    }

    

    // https://eigen.tuxfamily.org/dox/group__TutorialLinearAlgebra.html

    UppLog() << "\n\nTutorial page 6 - Linear algebra and decompositions";

 

    UppLog() << "\n\nBasic linear solving    Ax = b";

    {

        Matrix3f A;

        Vector3f b;

        A << 1, 2, 3,  

             4, 5, 6,  

             7, 8,10;

        b << 3, 3, 4;

        UppLog() << "\nHere is the matrix A:\n" << A;

        UppLog() << "\nHere is the vector b:\n" << b;

        Vector3f x = A.colPivHouseholderQr().solve(b);

        UppLog() << "\nThe solution is:\n" << x;

    }

    {

        Matrix2f A, b;

        A << 2, -1, -1, 3;

        b << 1, 2, 3, 1;

        UppLog() << "\nHere is the matrix A:\n" << A;

        UppLog() << "\nHere is the right hand side b:\n" << b;

        Matrix2f x = A.ldlt().solve(b);

        UppLog() << "\nThe solution is:\n" << x;

    }

    UppLog() << "\n\nChecking if a solution really exists";

    {

        MatrixXd A = MatrixXd::Random(100,100);

        MatrixXd b = MatrixXd::Random(100,50);

        MatrixXd x = A.fullPivLu().solve(b);

        double relative_error = (A*x - b).norm() / b.norm(); // norm() is L2 norm

        UppLog() << "\nThe relative error is:\n" << relative_error;

    }

    UppLog() << "\n\nComputing eigenvalues and eigenvectors";

    {

        Matrix2f A;

        A << 1, 2, 2, 3;

        UppLog() << "\nHere is the matrix A:\n" << A;

        SelfAdjointEigenSolver<Matrix2f> eigensolver(A);

        UppLog() << "\nThe eigenvalues of A are:\n" << eigensolver.eigenvalues();

        UppLog() << "\nHere's a matrix whose columns are eigenvectors of A "

            << "corresponding to these eigenvalues:\n"

            << eigensolver.eigenvectors();

    }

    UppLog() << "\n\nComputing inverse and determinant";

    {

        Matrix3f A;

        A << 1, 2, 1,

             2, 1, 0,

            -1, 1, 2;

        UppLog() << "\nHere is the matrix A:\n" << A;

        UppLog() << "\nThe determinant of A is " << A.determinant();

        UppLog() << "\nThe inverse of A is:\n" << A.inverse();

    }

    UppLog() << "\n\nLeast squares solving";

    {

        MatrixXf A = MatrixXf::Random(5, 2);

        UppLog() << "\nHere is the matrix A:\n" << A;

        VectorXf b = VectorXf::Random(5);

        UppLog() << "\nHere is the right hand side b:\n" << b;

        UppLog() << "\nThe least-squares solution is:\n"

               << A.jacobiSvd(ComputeThinU | ComputeThinV).solve(b);

    }

    UppLog() << "\n\nSeparating the computation from the construction";

    {

        Matrix2f A, b;

        LLT<Matrix2f> llt;

        A << 2, -1, -1, 3;

        b << 1, 2, 3, 1;

        UppLog() << "\nHere is the matrix A:\n" << A;

        UppLog() << "\nHere is the right hand side b:\n" << b;

        UppLog() << "\nComputing LLT decomposition...";

        llt.compute(A);

        UppLog() << "\nThe solution is:\n" << llt.solve(b);

        A(1,1)++;

        UppLog() << "\nThe matrix A is now:\n" << A;

        UppLog() << "\nComputing LLT decomposition...";

        llt.compute(A);

        UppLog() << "\nThe solution is now:\n" << llt.solve(b);

    }

    UppLog() << "\n\nRank-revealing decompositions";

    {

        Matrix3f A;

        A << 1, 2, 5,

             2, 1, 4,

             3, 0, 3;

        UppLog() << "\nHere is the matrix A:\n" << A;

        FullPivLU<Matrix3f> lu_decomp(A);

        UppLog() << "\nThe rank of A is " << lu_decomp.rank();

        UppLog() << "\nHere is a matrix whose columns form a basis of the null-space of A:\n"

            << lu_decomp.kernel();

        UppLog() << "\nHere is a matrix whose columns form a basis of the column-space of A:\n"

            << lu_decomp.image(A); // yes, have to pass the original A

    }

    {

        Matrix2d A;

        A << 2, 1,

             2, 0.9999999999;

        FullPivLU<Matrix2d> lu(A);

        UppLog() << "\nBy default, the rank of A is found to be " << lu.rank();

        lu.setThreshold(1e-5);

        UppLog() << "\nWith threshold 1e-5, the rank of A is found to be " << lu.rank();

    }

    

    // https://eigen.tuxfamily.org/dox/group__TutorialReductionsVisitorsBroadcasting.html

    UppLog() << "\n\nTutorial page 7 - Reductions, visitors and broadcasting";

        

    UppLog() << "\n\nReductions";

    {

        Eigen::Matrix2d mat;

        mat << 1, 2,

               3, 4;

        UppLog() << "\nHere is mat.sum():       " << mat.sum();

        UppLog() << "\nHere is mat.prod():      " << mat.prod();

        UppLog() << "\nHere is mat.mean():      " << mat.mean();

        UppLog() << "\nHere is mat.minCoeff():  " << mat.minCoeff();

        UppLog() << "\nHere is mat.maxCoeff():  " << mat.maxCoeff();

        UppLog() << "\nHere is mat.trace():     " << mat.trace();

    }    

    UppLog() << "\n\nNorm computations";

    {

        VectorXf v(2);

        MatrixXf m(2,2), n(2,2);

        

        v << -1,

                 2;

        

        m << 1,-2,

            -3, 4;

        

        UppLog() << "\nv.squaredNorm() = " << v.squaredNorm();

        UppLog() << "\nv.norm() = " << v.norm();

        UppLog() << "\nv.lpNorm<1>() = " << v.lpNorm<1>();

        UppLog() << "\nv.lpNorm<Infinity>() = " << v.lpNorm<Infinity>();

        

        UppLog() << "\n";

        UppLog() << "\nm.squaredNorm() = " << m.squaredNorm();

        UppLog() << "\nm.norm() = " << m.norm();

        UppLog() << "\nm.lpNorm<1>() = " << m.lpNorm<1>();

        UppLog() << "\nm.lpNorm<Infinity>() = " << m.lpNorm<Infinity>();

    }

    UppLog() << "\n\nBoolean reductions";

    {

        ArrayXXf a(2,2);

        

        a << 1,2,

                3,4;

        

        UppLog() << "\n(a > 0).all()   = " << (a > 0).all();

        UppLog() << "\n(a > 0).any()   = " << (a > 0).any();

        UppLog() << "\n(a > 0).count() = " << (a > 0).count();

        UppLog() << "\n";

        UppLog() << "\n(a > 2).all()   = " << (a > 2).all();

        UppLog() << "\n(a > 2).any()   = " << (a > 2).any();

        UppLog() << "\n(a > 2).count() = " << (a > 2).count();

    }

    UppLog() << "\n\nVisitors";

    {

        Eigen::MatrixXf m(2,2);

        

        m << 1, 2,

                3, 4;

        

        //get location of maximum

        MatrixXf::Index maxRow, maxCol;

        float max = m.maxCoeff(&maxRow, &maxCol);

        

        //get location of minimum

        MatrixXf::Index minRow, minCol;

        float min = m.minCoeff(&minRow, &minCol);

        

        UppLog() << "\nMax: " << max <<  ", at: " << maxRow << "," << maxCol;

        UppLog() << "\nMin: " << min << ", at: " << minRow << "," << minCol;

    }

    UppLog() << "\n\nPartial reductions";

    {

        Eigen::MatrixXf mat(2,4);

        mat << 1, 2, 6, 9,

               3, 1, 7, 2;

        

        UppLog() << "\nColumn's maximum: \n" << mat.colwise().maxCoeff();

    }

    {

        Eigen::MatrixXf mat(2,4);

        mat << 1, 2, 6, 9,

               3, 1, 7, 2;

        

        UppLog() << "\nRow's maximum: \n" << mat.rowwise().maxCoeff();

    }

    UppLog() << "\n\nCombining partial reductions with other operations";

    {

        MatrixXf mat(2,4);

        mat << 1, 2, 6, 9,

               3, 1, 7, 2;

        

        MatrixXf::Index   maxIndex;

        float maxNorm = mat.colwise().sum().maxCoeff(&maxIndex);

        

        UppLog() << "\nMaximum sum at position " << maxIndex;

        

        UppLog() << "\nThe corresponding vector is: ";

        UppLog() << "\n" << mat.col( maxIndex );

        UppLog() << "\nAnd its sum is is: " << maxNorm;

    }

    UppLog() << "\n\nBroadcasting";

    {

        Eigen::MatrixXf mat(2,4);

        Eigen::VectorXf v(2);

        

        mat << 1, 2, 6, 9,

               3, 1, 7, 2;

             

        v << 0,

             1;

           

        //add v to each column of m

        mat.colwise() += v;

        

        UppLog() << "\nBroadcasting result: ";

        UppLog() << "\n" << mat;

    }

    {

        Eigen::MatrixXf mat(2,4);

        Eigen::VectorXf v(4);

        

        mat << 1, 2, 6, 9,

               3, 1, 7, 2;

             

        v << 0,1,2,3;

           

        //add v to each row of m

        mat.rowwise() += v.transpose();

        

        UppLog() << "\nBroadcasting result: ";

        UppLog() << "\n" << mat;

    }

    UppLog() << "\n\nCombining broadcasting with other operations";

    {

        Eigen::MatrixXf m(2,4);

        Eigen::VectorXf v(2);

        

        m << 1, 23, 6, 9,

                3, 11, 7, 2;

           

        v << 2,

             3;

        

        MatrixXf::Index index;

        // find nearest neighbour

        (m.colwise() - v).colwise().squaredNorm().minCoeff(&index);

        

        UppLog() << "\nNearest neighbour is column " << index << ":";

        UppLog() << "\n" << m.col(index);

    }

        

    UppLog() << "\n\nSerializing tests";

    {

        SerialTest serialTest, serialTest_j, serialTest_x, serialTest_s;

        serialTest.m << 1, 2,

                         4, 8;

        serialTest.v << 1, 2, 4;

 

        StoreAsJsonFile(serialTest, GetExeDirFile("Json.txt"));

        LoadFromJsonFile(serialTest_j, GetExeDirFile("Json.txt"));

        UppLog() << "\nJSON demo";

        serialTest_j.Print();

        

        StoreAsXMLFile(serialTest, "XMLdata", GetExeDirFile("Xml.txt"));

        LoadFromXMLFile(serialTest_x, GetExeDirFile("Xml.txt"));

        UppLog() << "\nXML demo";

        serialTest_x.Print();

        

        StoreToFile(serialTest, GetExeDirFile("Serial.dat"));

        LoadFromFile(serialTest_s, GetExeDirFile("Serial.dat"));

        UppLog() << "\nSerialization demo";

        serialTest_s.Print();

    }

    

    NonLinearTests();

    

    FFTTests();

    

    #ifdef flagDEBUG

    Cout() << "\nPress enter key to end";

    ReadStdIn();

    #endif   

}

 

 

 

 

 

 

fft.cpp

 

#include <Core/Core.h>

#include <plugin/Eigen/Eigen.h>

 

using namespace Upp;

 

#ifdef USE_FFTW

#include <fftw3.h>

#endif

 

using namespace Eigen;

 

 

void FFTTests()

{

    UppLog() << "\nFFT sample\nGets the FFT of equation"

              "\n f(t) = 2*sin(2*PI*t/50 - PI/3) + 5*sin(2*PI*t/30 - PI/2) + 30*sin(2*PI*t/10 - PI/5)"

              "\nsampled with a frequency of 14 samples/second";

    

    int numData = 8000;

    double samplingFrecuency = 14;    

 

   // Filling the data series

   VectorXd timebuf(numData);

   {

        double t = 0;

        for (int i = 0; i < numData; ++i, t = i/samplingFrecuency)

               timebuf[i] = 2*sin(2*M_PI*t/50 - M_PI/3) + 5*sin(2*M_PI*t/30 - M_PI/2) + 30*sin(2*M_PI*t/10 - M_PI/5);

   }

   

   // FFT

   VectorXcd freqbuf;

   FFT<double> fft;

   fft.SetFlag(fft.HalfSpectrum);

   fft.fwd(freqbuf, timebuf);

    

    // Filter the FFT. Frequencies between 1/25 and 1/35 Hz are removed

    // Original FFT is not changed for saving it later

    VectorXcd freqbuf2(freqbuf.size());

    {

        for (int i = 0; i < freqbuf.size(); ++i) {

            double freq = i*samplingFrecuency/numData;

            double T = 1/freq;

            if (T > 25 && T < 35)

                freqbuf2[i] = 0;

            else

                freqbuf2[i] = freqbuf[i];

        }

    }

    

    // Inverse filtered FFT to get filtered series

    VectorXd timebuf2(numData);

    fft.inv(timebuf2, freqbuf2);

    

    String csvSep = ";";

 

    // Saving original and filtered FFT

    {

        String str;

        str << "Frec" << csvSep << "T" << csvSep << "fft" << csvSep << "Filtered fft";

        for (int i = 0; i < freqbuf.size(); ++i) {

            double freq = i*samplingFrecuency/numData;

            double T = 1/freq;

            str << "\n" << freq << csvSep << (freq > 0 ? FormatDouble(T) : "") << csvSep

                        << 2*std::abs(freqbuf[i])/numData << csvSep

                        << 2*std::abs(freqbuf2[i])/numData;

        }

        String fftFileName = GetExeDirFile("fft.csv");

        UppLog() << "\nFFT saved in '" << fftFileName << "'";

        VERIFY(SaveFile(fftFileName, str));

    }

    

    // Saving original and filtered series

    {

        String str;

        str << "Time" << csvSep << "Data" << csvSep << "Filtered data";

        double t = 0;

        for (int i = 0; i < numData; ++i, t = i*1/samplingFrecuency)

               str << "\n" << t << csvSep << timebuf[i] << csvSep << timebuf2[i];;

        String dataFileName = GetExeDirFile("data.csv");

        UppLog() << "\nSource data saved in '" << dataFileName << "'";

        VERIFY(SaveFile(dataFileName, str));

   }

}

 

 

 

 

non-linear.cpp

 

#include <Core/Core.h>

#include <plugin/Eigen/Eigen.h>

 

using namespace Upp;

 

using namespace Eigen;

 

#define VerifyIsApprox(a, b)    (fabs(a-(b)) < 0.0001)

 

struct Eckerle4_functor : NonLinearOptimizationFunctor<double> {

    Eckerle4_functor() : NonLinearOptimizationFunctor(3, 35) {}

    static const double x[35];

    static const double y[35];

    int operator()(const VectorXd &b, VectorXd &fvec) const {

        ASSERT(b.size() == unknowns);

        ASSERT(fvec.size() == datasetLen);

        for(int i = 0; i < 35; i++)

            fvec[i] = b[0]/b[1] * exp(-0.5*(x[i]-b[2])*(x[i]-b[2])/(b[1]*b[1])) - y[i];

        return 0;

    }

};

const double Eckerle4_functor::x[35] = {400.0, 405.0, 410.0, 415.0, 420.0, 425.0, 430.0, 435.0, 436.5, 438.0, 439.5, 441.0, 442.5, 444.0, 445.5, 447.0, 448.5, 450.0, 451.5, 453.0, 454.5, 456.0, 457.5, 459.0, 460.5, 462.0, 463.5, 465.0, 470.0, 475.0, 480.0, 485.0, 490.0, 495.0, 500.0};

const double Eckerle4_functor::y[35] = {0.0001575, 0.0001699, 0.0002350, 0.0003102, 0.0004917, 0.0008710, 0.0017418, 0.0046400, 0.0065895, 0.0097302, 0.0149002, 0.0237310, 0.0401683, 0.0712559, 0.1264458, 0.2073413, 0.2902366, 0.3445623, 0.3698049, 0.3668534, 0.3106727, 0.2078154, 0.1164354, 0.0616764, 0.0337200, 0.0194023, 0.0117831, 0.0074357, 0.0022732, 0.0008800, 0.0004579, 0.0002345, 0.0001586, 0.0001143, 0.0000710 };

 

struct Thurber_functor : NonLinearOptimizationFunctor<double> {

    static VectorXd _x, _y;

    

    Thurber_functor() : NonLinearOptimizationFunctor() {

        double x[] = {-3.067, -2.981, -2.921, -2.912, -2.840, -2.797, -2.702, -2.699, -2.633, -2.481, -2.363, -2.322, -1.501, -1.460, -1.274, -1.212, -1.100, -1.046, -0.915, -0.714, -0.566, -0.545, -0.400, -0.309, -0.109, -0.103, 0.010, 0.119, 0.377, 0.790, 0.963, 1.006, 1.115, 1.572, 1.841, 2.047, 2.200};

        double y[] = {80.574, 84.248, 87.264, 87.195, 89.076, 89.608, 89.868, 90.101, 92.405, 95.854, 100.696, 101.060, 401.672, 390.724, 567.534, 635.316, 733.054, 759.087, 894.206, 990.785, 1090.109, 1080.914, 1122.643, 1178.351, 1260.531, 1273.514, 1288.339, 1327.543, 1353.863, 1414.509, 1425.208, 1421.384, 1442.962, 1464.350, 1468.705, 1447.894, 1457.628};        

        _x = Map<VectorXd>(x, sizeof(x)/sizeof(double));

        _y = Map<VectorXd>(y, sizeof(y)/sizeof(double));

        unknowns = 7;

        datasetLen = _x.size();

    }

    

    int operator()(const VectorXd &b, VectorXd &fvec) const {

        ASSERT(b.size() == unknowns);

        ASSERT(fvec.size() == datasetLen);

        for(int i = 0; i < datasetLen; i++) {

            double x = _x[i], xx = x*x, xxx = xx*x;

            fvec[i] = (b[0] + b[1]*x + b[2]*xx + b[3]*xxx)/(1. + b[4]*x + b[5]*xx + b[6]*xxx) - _y[i];

        }

        return 0;

    }

};

 

VectorXd Thurber_functor::_x, Thurber_functor::_y;

 

void NonLinearOptimization() {

    UppLog() << "\n\nNon linear equations optimization using Levenberg Marquardt based on Minpack\n"

         "(Given a set of non linear equations and a set of data longer than the equations, "

         "the program finds the equation coefficients that better fit with the equations)";

    

    {

        UppLog() << "\n\nEckerle4 equation\nSee: http://www.itl.nist.gov/div898/strd/nls/data/eckerle4.shtml";

        

        VectorXd x(3);

        

        x << 1., 10., 500.;        // Initial values

        

        Eckerle4_functor functor;

        NumericalDiff<Eckerle4_functor> numDiff(functor);

        LevenbergMarquardt<NumericalDiff<Eckerle4_functor> > lm(numDiff);

        int ret = lm.minimize(x);

        if (ret == LevenbergMarquardtSpace::ImproperInputParameters ||

            ret == LevenbergMarquardtSpace::TooManyFunctionEvaluation)

            UppLog() << "\nNo convergence!: " << ret;

        else {

            if (VerifyIsApprox(lm.fvec.squaredNorm(), 1.4635887487E-03))

                UppLog() << "\nNorm^2 is right";

            else

                UppLog() << "\nNorm^2 is NOT right";

            if (VerifyIsApprox(x[0], 1.5543827178) &&

                VerifyIsApprox(x[1], 4.0888321754) &&

                VerifyIsApprox(x[2], 4.5154121844E+02))

                UppLog() << "\nNon-linear function optimization is right!";

            else

                UppLog() << "\nNon-linear function optimization is NOT right!";

        }

    }

    {

        UppLog() << "\n\nThis is a simpler way, using NonLinearOptimization()";

        double x[] = {400.0, 405.0, 410.0, 415.0, 420.0, 425.0, 430.0, 435.0, 436.5, 438.0, 439.5, 441.0, 442.5, 444.0, 445.5, 447.0, 448.5, 450.0, 451.5, 453.0, 454.5, 456.0, 457.5, 459.0, 460.5, 462.0, 463.5, 465.0, 470.0, 475.0, 480.0, 485.0, 490.0, 495.0, 500.0};

        double f[] = {0.0001575, 0.0001699, 0.0002350, 0.0003102, 0.0004917, 0.0008710, 0.0017418, 0.0046400, 0.0065895, 0.0097302, 0.0149002, 0.0237310, 0.0401683, 0.0712559, 0.1264458, 0.2073413, 0.2902366, 0.3445623, 0.3698049, 0.3668534, 0.3106727, 0.2078154, 0.1164354, 0.0616764, 0.0337200, 0.0194023, 0.0117831, 0.0074357, 0.0022732, 0.0008800, 0.0004579, 0.0002345, 0.0001586, 0.0001143, 0.0000710};    

        int num = sizeof(x)/sizeof(double);

        VectorXd y(3);

        y << 1., 10., 500.;

        VERIFY(NonLinearOptimization(y, num, [&](const VectorXd &y, VectorXd &residual)->int {

            for(int i = 0; i < num; i++)

                residual[i] = y[0]/y[1] * exp(-0.5*(x[i]-y[2])*(x[i]-y[2])/(y[1]*y[1])) - f[i];

            return 0;    

        }));

        

        VERIFY(VerifyIsApprox(y[0], 1.5543827178) &&

                VerifyIsApprox(y[1], 4.0888321754) &&

                VerifyIsApprox(y[2], 4.5154121844E+02));

    }

    {

        UppLog() << "\n\nThurber equation\nSee: http://www.itl.nist.gov/div898/strd/nls/data/thurber.shtml\n";

        

        VectorXd x(7);

        x << 1000, 1000, 400, 40, 0.7, 0.3, 0.0;    // Initial values

        

        Thurber_functor functor;    // Vars and equations known at run time

        NumericalDiff<Thurber_functor> numDiff(functor);

        LevenbergMarquardt<NumericalDiff<Thurber_functor> > lm(numDiff);

        lm.parameters.ftol = 1.E4*NumTraits<double>::epsilon();

        lm.parameters.xtol = 1.E4*NumTraits<double>::epsilon();

        int ret = lm.minimize(x);

        

        VERIFY(!(ret == LevenbergMarquardtSpace::ImproperInputParameters ||

            ret == LevenbergMarquardtSpace::TooManyFunctionEvaluation));

        

        VERIFY(VerifyIsApprox(lm.fvec.squaredNorm(), 5.64270823E+03));

        

        VERIFY(VerifyIsApprox(x[0], 1.2881396E+03) &&

                VerifyIsApprox(x[1], 1.4910792E+03) &&

                VerifyIsApprox(x[2], 5.8323836E+02) &&

                VerifyIsApprox(x[3], 7.5416644E+01) &&

                VerifyIsApprox(x[4], 9.6629502E-01) &&

                VerifyIsApprox(x[5], 3.9797285E-01) &&

                VerifyIsApprox(x[6], 4.9727297E-02));

    }

}

 

struct Hybrd_functor : NonLinearOptimizationFunctor<double>

{

    Hybrd_functor() : NonLinearOptimizationFunctor<double>(9, 9) {}

    int operator()(const VectorXd &x, VectorXd &fvec) const {

        double temp, temp1, temp2;

        const ptrdiff_t n = x.size();

        

        ASSERT(fvec.size() == n);

        for (ptrdiff_t k = 0; k < n; k++) {

            temp = (3. - 2.*x[k])*x[k];

            temp1 = 0.;

            if (k)

                temp1 = x[k-1];

            temp2 = 0.;

            if (k != n-1)

                temp2 = x[k+1];

            fvec[k] = temp - temp1 - 2.*temp2 + 1.;

        }

        return 0;

    }

};

 

void NonLinearSolving() {

    UppLog() << "\n\nNon linear equation solving using the Powell hybrid method (\"dogleg\") based on Minpack. "

           << "(Finds a zero of a system of n nonlinear equations in n variables)";

    

    const int n = 9;

    VectorXd x;

    

    x.setConstant(n, -1.);        // Initial values

    

    Hybrd_functor functor;

    HybridNonLinearSolver<Hybrd_functor> solver(functor);

    int ret = solver.solveNumericalDiff(x);

    

    VERIFY(!(ret == HybridNonLinearSolverSpace::ImproperInputParameters ||

        ret == HybridNonLinearSolverSpace::TooManyFunctionEvaluation ||

        ret == HybridNonLinearSolverSpace::NotMakingProgressJacobian ||

        ret == HybridNonLinearSolverSpace::NotMakingProgressIterations));

 

    VERIFY(VerifyIsApprox(solver.fvec.blueNorm(), 1.192636e-08));

            

    VERIFY(VerifyIsApprox(x[0], -0.5706545) &&

            VerifyIsApprox(x[1], -0.6816283) &&     

            VerifyIsApprox(x[2], -0.7017325) &&

            VerifyIsApprox(x[3], -0.7042129) &&  

            VerifyIsApprox(x[4], -0.701369) &&

            VerifyIsApprox(x[5], -0.6918656) &&

            VerifyIsApprox(x[6], -0.665792) &&

            VerifyIsApprox(x[7], -0.5960342) &&  

            VerifyIsApprox(x[8], -0.4164121));

    

    UppLog() << "\n\nThis is a simpler way, using SolveNonLinearEquations()";

    

    x.setConstant(n, -1.);        // Initial values

    

    VERIFY(SolveNonLinearEquations(x, [&](const VectorXd &x, VectorXd &residual)->int {

        const ptrdiff_t n = x.size();

        

        ASSERT(residual.size() == n);

        for (ptrdiff_t k = 0; k < n; k++) {

            double temp1 = 0.;

            if (k)

                temp1 = x[k-1];

            double temp2 = 0.;

            if (k != n-1)

                temp2 = x[k+1];

            residual[k] = (3. - 2.*x[k])*x[k] - temp1 - 2.*temp2 + 1.;

        }

        return 0;

    }));    // No convergence!

    VERIFY(VerifyIsApprox(x[0], -0.5706545) &&

            VerifyIsApprox(x[1], -0.6816283) &&     

            VerifyIsApprox(x[2], -0.7017325) &&

            VerifyIsApprox(x[3], -0.7042129) &&  

            VerifyIsApprox(x[4], -0.701369) &&

            VerifyIsApprox(x[5], -0.6918656) &&

            VerifyIsApprox(x[6], -0.665792) &&

            VerifyIsApprox(x[7], -0.5960342) &&  

            VerifyIsApprox(x[8], -0.4164121));

}

 

void NonLinearTests() {

    NonLinearSolving();

    NonLinearOptimization();

}

 

 

 

 

 

 

 

Do you want to contribute?