Levenbergmarquardt algorithm Eigen Implementation (c + +) __ algorithm

VS2013 PCL1.7.2 uses a eigen library

#include <iostream> #include <Eigen/Dense> #include <unsupported/Eigen/NonLinearOptimization> struct Myfunctor {int operator () (const EIGEN::VECTORXF &x, EIGEN::VECTORXF &fvec) const {Double
        TMP1, TMP2, Tmp3; static const double Y[15] = {1.4e-1, 1.8e-1, 2.2e-1, 2.5e-1, 2.9e-1, 3.2e-1, 3.5e-1, 3.9e-1, 3.7e-1, 5.8e-1,

        7.3e-1, 9.6e-1, 1.34, 2.1, 4.39};
            for (int i = 0; I < values (); i++) {tmp1 = i + 1;
            TMP2 = 16-i-1; Tmp3 = (i >= 8)?
            TMP2:TMP1;
        Fvec[i] = y[i]-(x[0] + tmp1/(x[1) * TMP2 + x[2] * Tmp3));
    return 0; 
        } int df (const EIGEN::VECTORXF &x, EIGEN::MATRIXXF &fjac) const {Double tmp1, tmp2, Tmp3, Tmp4;
            for (int i = 0; I < values (); i++) {tmp1 = i + 1;
            TMP2 = 16-i-1; Tmp3 = (i >= 8)?
            TMP2:TMP1; Tmp4 = (x[1] * TMP2 + X[2] * Tmp3);
            Tmp4 = Tmp4*tmp4;
            Fjac (i, 0) =-1;
            Fjac (i, 1) = Tmp1*tmp2/tmp4;
        Fjac (i, 2) = Tmp1*tmp3/tmp4;
    return 0;
    int inputs () const {return 3;}


int values () const {return 3;}//Number of constraints};
    int main (int argc, char *argv[]) {EIGEN::VECTORXF x (3);
    X (0) = 10;
    X (1) =-2;
    X (2) = 3;

    X (2) = 0;

    Std::cout << "x:" << x << Std::endl;
    Myfunctor functor;

    Eigen::levenbergmarquardt<myfunctor, float> lm (functor);

    Lm.minimize (x);

    Std::cout << "x that minimizes the function:" << x << Std::endl;
    System ("pause");
return 0; }

Results

Sample2: (for reference only)

#include <iostream> #include <Eigen/Dense> #include <unsupported/Eigen/NonLinearOptimization> #
    Include "Mex.h" class Ph_module {//functor private:double H, H2co3, HCO3, CO2, CO3, NH3, NH4, HS, H2S, OH;
    Double Fe, Ca, NO3, SO4, PO4;

Double Kc0, KC1, KC2, KNH, KHS, Kw;
        Public:ph_module (eigen::vectorxf x0, eigen::vectorxf other_conc) {H = x0 (0) *x0 (0);
        HCO3 = x0 (1) *x0 (1);
        CO2 = x0 (2) *x0 (2);
        CO3 = x0 (3) *x0 (3);
        NH3 = x0 (4) *x0 (4);
        NH4 = x0 (5) *x0 (5);
        HS = x0 (6) *x0 (6);
        H2S = x0 (7) *x0 (7);
        OH = x0 (8) *x0 (8);
        H2CO3 = x0 (9) *x0 (9);
        Fe = Other_conc (0) *other_conc (0);
        Ca = Other_conc (1) *other_conc (1);
        NO3 = Other_conc (2) *other_conc (2);
        SO4 = Other_conc (3) *other_conc (3);
        PO4 = Other_conc (4) *other_conc (4);
        KC0 = 1.7 *pow (10.0,-3.0);
        KC1 = 5.01 *pow (10.0, -7.0) *pow (10.0, 6.0); KC2 = 4.78 *pow (10.0, -11.0) *pow (10.0, 6.0);
        KNH = 5.62 *pow (10.0, -10.0) *pow (10.0, 6.0);
        KHS = 1.3 *pow (10.0, -7.0) *pow (10.0, 6.0);
    Kw = POW ( -14) *pow (10, 12);

    };
    int operator () (const EIGEN::VECTORXF &x, EIGEN::VECTORXF &fvec) const;
    int DF (const EIGEN::VECTORXF &x, EIGEN::MATRIXXF &fjac) const;
    int inputs () const;//inputs is the dimension of X. int values () const;

"Values" is the dimension of F};  int Ph_module::operator () (const EIGEN::VECTORXF &x, EIGEN::VECTORXF &fvec) const {//Equations:fvec (0) =
    X (9) *x (9)-kc0*x (2) *x (2);
    Fvec (1) = x (0) *x (0) *x (1) *x (1)-kc1*x (9) *x (9);
    Fvec (2) = x (0) *x (0) *x (3) *x (3)-kc2*x (1) *x (1);
    Fvec (3) = x (0) *x (0) *x (4) *x (4)-knh*x (5) *x (5);
    Fvec (4) = x (0) *x (0) *x (6) *x (6)-khs*x (7) *x (7);
    Fvec (5) = x (0) *x (0) *x (8) *x (8)-Kw;
    Fvec (6) = x (4) *x (4) + x (5) *x (5)-nh4-nh3;
    Fvec (7) = x (6) *x (6) + x (7) *x (7)-HS-H2S; Fvec (8) = x (3) *x (3) + x (1) *x (1) + x (2) *x (2) + x (9) *x (9)-Co3-hco3- Co2-h2co3;  Fvec (9) = x (0) *x (0) + x (5) *x (5) + 2 * Fe + 2 * Ca-(x (1) *x (1) + 2 * x (3) *x (3) + x (6) *x (6) + x (8) *x (8) + NO3 + 2 * SO4 + 3
    * PO4);
return 0; int Ph_module::d f (const EIGEN::VECTORXF &x, EIGEN::MATRIXXF &fjac) const {//Implement Jacobian//RO W #0 Fjac (0, 0) = 0; Fjac (0, 1) = 0; Fjac (0, 2) =-2 * KC0*X (2); Fjac (0, 3) = 0; Fjac (0, 4) = 0; Fjac (0, 5) = 0; Fjac (0, 6) = 0; Fjac (0, 7) = 0; Fjac (0, 8) = 0;

    Fjac (0, 9) = 2 * x (9); Row #1 Fjac (1, 0) = 2 * x (0) *x (1) *x (1); Fjac (1, 1) = 2 * x (0) *x (0) *x (1); Fjac (1, 2) = 0; Fjac (1, 3) = 0; Fjac (1, 4) = 0; Fjac (1, 5) = 0; Fjac (1, 6) = 0; Fjac (1, 7) = 0; Fjac (1, 8) = 0;

    Fjac (1, 9) =-2 * KC1*X (9); Row #2 Fjac (2, 0) = 2 * x (0) *x (3) *x (3); Fjac (2, 1) =-2 * kc2*x (1); Fjac (2, 2) = 0; Fjac (2, 3) = 2 * x (0) *x (0) *x (3); Fjac (2, 4) = 0; Fjac (2, 5) = 0; Fjac (2, 6) = 0; Fjac (2, 7) = 0; Fjac (2, 8) = 0;

    Fjac (2, 9) = 0; Row #3 Fjac (3, 0) = 2 * x (0) *x (4) *x (4); Fjac (3, 1) = 0; Fjac (3, 2) = 0; Fjac (3, 3) = 0; Fjac (3, 4) = 2 * x (0) *x (0) *x (4); Fjac (3, 5) =-2 * KNH*X (5); Fjac (3, 6) = 0; Fjac (3, 7) = 0; Fjac (3, 8) = 0;

    Fjac (3, 9) = 0; Row #4 Fjac (4, 0) = 2 * x (0) *x (6) *x (6); Fjac (4, 1) = 0; Fjac (4, 2) = 0; Fjac (4, 3) = 0; Fjac (4, 4) = 0; Fjac (4, 5) = 0; Fjac (4, 6) = 2 * x (0) *x (0) *x (6); Fjac (4, 7) = 2 * KHS*X (7); Fjac (4, 8) = 0;

    Fjac (4, 9) = 0; Row #5 Fjac (5, 0) = 2 * x (0) *x (8) *x (8); Fjac (5, 1) = 0; Fjac (5, 2) = 0; Fjac (5, 3) = 0; Fjac (5, 4) = 0; Fjac (5, 5) = 0; Fjac (5, 6) = 0; Fjac (5, 7) = 0; Fjac (5, 8) = 2 * x (0) *x (0) *x (8);

    Fjac (5, 9) = 0; Row #6 Fjac (6, 0) = 0; Fjac (6, 1) = 0; Fjac (6, 2) = 0; Fjac (6, 3) = 0; Fjac (6, 4) = 2 * x (4); Fjac (6, 5) = 2 * x (5); Fjac (6, 6) = 0; Fjac (6, 7) = 0; Fjac (6, 8) = 0;

    Fjac (6, 9) = 0; Row #7 Fjac (7, 0) = 0; Fjac (7, 1) = 0; Fjac (7, 2) = 0; Fjac (7, 3) = 0; Fjac (7, 4) = 0; Fjac (7, 5) = 0; Fjac (7, 6) = 2 * x (6); Fjac (7, 7) = 2 * x (7); Fjac (7, 8) = 0;

    Fjac (7, 9) = 0; Row #8 Fjac (8, 0) = 0; Fjac (8, 1) = 2 * x (1); Fjac (8, 2) = 2 * x (2); Fjac (8, 3) = 2 * x (3); Fjac (8, 4) = 0; Fjac (8, 5) = 0; Fjac (8, 6) = 0; Fjac (8, 7) = 0; Fjac (8, 8) = 0;

    Fjac (8, 9) = 2 * x (9); Row #9 Fjac (9, 0) = 2 * x (0); Fjac (9, 1) =-2 * x (1); Fjac (9, 2) = 0; Fjac (9, 3) =-4 * x (3); Fjac (9, 4) = 0; Fjac (9, 5) = 2 * x (5); Fjac (9, 6) =-2 * x (6); Fjac (9, 7) = 0; Fjac (9, 8) =-2 * x (8);


    Fjac (9, 9) = 0;
return 0; the int ph_module::inputs () const {return]//inputs is the dimension of X. int ph_module::values () const {return 10 ;

    }//"Values" is the number of f_i void mexfunction (int nlhs, Mxarray *plhs[], int nrhs, const Mxarray *prhs[]) {

    Double *data = (double *) mxgetdata (prhs[0]);
    EIGEN::VECTORXF x0 (10);

    EIGEN::VECTORXF Other_conc (5);
    Forming VEC of unknowns for (int i = 0; i < ++i) {x0 (i) = data[i];
    }//Vec of ions for (int i = 0; i < 5; ++i) {Other_conc (i) = Data[i + 10]; } ph_module Functor (x0, other_conc);
    Eigen::levenbergmarquardt<ph_module, float> lm (functor);
    Lm.parameters.xtol = 1e-8;
    Lm.parameters.ftol = 1e-8;

    int info = lm.minimize (x0);
    Assign output PLHS[0] = Mxcreatedoublematrix (1, mxreal);
    Double * z = MXGETPR (plhs[0]);
    for (int i = 0; i < ++i) {Z[i] = x0 (i); }
}