Show
Ignore:
Timestamp:
06/09/10 14:00:40 (15 years ago)
Author:
mido
Message:

astyle applied all over the library

Location:
library/tests/stresssuite
Files:
19 modified

Legend:

Unmodified
Added
Removed

  • library/tests/stresssuite/additive_generator.cpp

    r717 r1064  
    66 
    77void additive_generator::from_setting ( const Setting &set ) { 
    8         int sz; 
    9         UI::get ( sz, set, "size", UI::compulsory ); 
     8    int sz; 
     9    UI::get ( sz, set, "size", UI::compulsory ); 
    1010 
    11         mat a0 = randu ( sz, sz ); 
    12         a = a0 * a0.T(); 
    13         vec v = randu ( sz ); 
    14         v2 = outer_product ( v, v ); 
     11    mat a0 = randu ( sz, sz ); 
     12    a = a0 * a0.T(); 
     13    vec v = randu ( sz ); 
     14    v2 = outer_product ( v, v ); 
    1515 
    16         UI::get ( lambda, set, "lambda", UI::optional ); 
     16    UI::get ( lambda, set, "lambda", UI::optional ); 
    1717} 
    1818 
    1919mat additive_generator::next() { 
    20         mat b = a; 
    21         a = ( 1 - lambda ) * a + lambda * v2; 
    22         return b; 
     20    mat b = a; 
     21    a = ( 1 - lambda ) * a + lambda * v2; 
     22    return b; 
    2323} 

  • library/tests/stresssuite/additive_generator.h

    r717 r1064  
    2121class additive_generator : public generator { 
    2222private: 
    23         itpp::mat a; 
    24         itpp::mat v2; 
    25         double lambda; 
     23    itpp::mat a; 
     24    itpp::mat v2; 
     25    double lambda; 
    2626 
    2727public: 
    28         additive_generator() : lambda ( 0.5 ) { } 
     28    additive_generator() : lambda ( 0.5 ) { } 
    2929 
    30         itpp::mat next(); 
     30    itpp::mat next(); 
    3131 
    32         //! Load from structure with elements: 
    33         //!  \code 
    34         //! { size = 7; // size (rows == cols) of the generated matrix 
    35         //!   lambda = 0.5; // weight of the added vector 
    36         //! } 
    37         //! \endcode 
    38         //! size is mandatory, lambda optional (with default as shown). 
    39         void from_setting ( const Setting &set ); 
     32    //! Load from structure with elements: 
     33    //!  \code 
     34    //! { size = 7; // size (rows == cols) of the generated matrix 
     35    //!   lambda = 0.5; // weight of the added vector 
     36    //! } 
     37    //! \endcode 
     38    //! size is mandatory, lambda optional (with default as shown). 
     39    void from_setting ( const Setting &set ); 
    4040}; 
    4141 

  • library/tests/stresssuite/arx_elem_stress.cpp

    r1009 r1064  
    55 
    66TEST ( arx_elem_stress ) { 
    7         // Setup model : ARX for 1D Gaussian 
    8         //Test constructor 
    9         mat V0 = 0.00001 * eye ( 2 ); 
    10         V0 ( 0, 0 ) = 0.1; // 
    11         ARX Ar; 
    12         Ar.set_statistics ( 1, V0, -1.0 ); 
    13         Ar.set_constant ( true ); 
    14         Ar.validate(); 
     7    // Setup model : ARX for 1D Gaussian 
     8    //Test constructor 
     9    mat V0 = 0.00001 * eye ( 2 ); 
     10    V0 ( 0, 0 ) = 0.1; // 
     11    ARX Ar; 
     12    Ar.set_statistics ( 1, V0, -1.0 ); 
     13    Ar.set_constant ( true ); 
     14    Ar.validate(); 
    1515 
    16         mat mu ( 1, 1 ); 
    17         mat R ( 1, 1 ); 
    18         Ar.posterior().mean_mat ( mu, R ); 
    19         cout << "Prior moments: mu=" << mu << ", R=" << R << endl; 
     16    mat mu ( 1, 1 ); 
     17    mat R ( 1, 1 ); 
     18    Ar.posterior().mean_mat ( mu, R ); 
     19    cout << "Prior moments: mu=" << mu << ", R=" << R << endl; 
    2020 
    21         int ndat = 200; 
    22         vec smp = randn ( ndat ); 
    23         // 
    24         mat Smp = ones ( 2, ndat ); 
    25         Smp.set_row ( 0, smp ); 
    26         // 
    27         Ar.bayes_batch ( Smp ); 
    28         // Ar is now filled with estimates of N(0,1); 
    29         cout << "Empirical moments: mu=" << sum ( smp ) / ndat << ", R=" << sum_sqr ( smp ) / ndat - pow ( sum ( smp ) / ndat, 2 ) << endl; 
    30         Ar.posterior().mean_mat ( mu, R ); 
    31         cout << "Posterior moments: mu=" << mu << ", R=" << R << endl; 
     21    int ndat = 200; 
     22    vec smp = randn ( ndat ); 
     23    // 
     24    mat Smp = ones ( 2, ndat ); 
     25    Smp.set_row ( 0, smp ); 
     26    // 
     27    Ar.bayes_batch ( Smp ); 
     28    // Ar is now filled with estimates of N(0,1); 
     29    cout << "Empirical moments: mu=" << sum ( smp ) / ndat << ", R=" << sum_sqr ( smp ) / ndat - pow ( sum ( smp ) / ndat, 2 ) << endl; 
     30    Ar.posterior().mean_mat ( mu, R ); 
     31    cout << "Posterior moments: mu=" << mu << ", R=" << R << endl; 
    3232 
    33         //////// TEST prediction 
    34         vec x = linspace ( -3.0, 3.0, 100 ); 
    35         double xstep = 6.0 / 100.0; 
    36         mat X ( 1, 100 ); 
    37         mat X2 ( 2, 100 ); 
    38         X.set_row ( 0, x ); 
    39         X2.set_row ( 0, x ); 
     33    //////// TEST prediction 
     34    vec x = linspace ( -3.0, 3.0, 100 ); 
     35    double xstep = 6.0 / 100.0; 
     36    mat X ( 1, 100 ); 
     37    mat X2 ( 2, 100 ); 
     38    X.set_row ( 0, x ); 
     39    X2.set_row ( 0, x ); 
    4040 
    41         mlstudent* Ap = Ar.predictor_student(); 
    42         vec Ap_x = Ap->evallogcond_mat ( X, empty_vec ); 
    43         vec ll_x = Ar.logpred_mat ( X2 , empty_vec); 
     41    mlstudent* Ap = Ar.predictor_student(); 
     42    vec Ap_x = Ap->evallogcond_mat ( X, empty_vec ); 
     43    vec ll_x = Ar.logpred_mat ( X2 , empty_vec); 
    4444 
    45         cout << "normalize : " << xstep*sum ( exp ( Ap_x ) ) << endl; 
    46         cout << "normalize : " << xstep*sum ( exp ( ll_x ) ) << endl; 
     45    cout << "normalize : " << xstep*sum ( exp ( Ap_x ) ) << endl; 
     46    cout << "normalize : " << xstep*sum ( exp ( ll_x ) ) << endl; 
    4747 
    48         it_file it ( "arx_elem_test.it" ); 
    49         it << Name ( "Ap_x" ) << Ap_x; 
    50         it << Name ( "ll_x" ) << ll_x; 
     48    it_file it ( "arx_elem_test.it" ); 
     49    it << Name ( "Ap_x" ) << Ap_x; 
     50    it << Name ( "ll_x" ) << ll_x; 
    5151} 

  • library/tests/stresssuite/arx_stress.cpp

    r766 r1064  
    1717 
    1818TEST ( arx_stress ) { 
    19         // Setup model 
    20         vec th ( "0.8 -0.3 0.4 0.01" ); 
    21         int ord = th.length(); //auxiliary variable 
    22         double sqr = 0.1; 
     19    // Setup model 
     20    vec th ( "0.8 -0.3 0.4 0.01" ); 
     21    int ord = th.length(); //auxiliary variable 
     22    double sqr = 0.1; 
    2323 
    24         //Test constructor 
    25         mat V0 = 0.00001 * eye ( ord + 1 ); 
    26         V0 ( 0 ) = 1; // 
    27         double nu0 = ord + 5.0; 
     24    //Test constructor 
     25    mat V0 = 0.00001 * eye ( ord + 1 ); 
     26    V0 ( 0 ) = 1; // 
     27    double nu0 = ord + 5.0; 
    2828 
    29         ARX Ar; 
    30         Ar.set_statistics ( 1, V0, nu0 );               // Estimator 
    31         Ar.set_constant ( false ); 
    32         Ar.validate(); 
    33         const epdf& f_thr = Ar.posterior();          // refrence to posterior of the estimator 
     29    ARX Ar; 
     30    Ar.set_statistics ( 1, V0, nu0 );               // Estimator 
     31    Ar.set_constant ( false ); 
     32    Ar.validate(); 
     33    const epdf& f_thr = Ar.posterior();          // refrence to posterior of the estimator 
    3434 
    35         //Test estimation 
    36         int ndat = 100;                 // number of data records 
    37         vec Yt ( ndat );                // Store generated data 
    38         Yt.set_subvector ( 0, randn ( ord ) ); //initial values 
    39         vec rgr ( ord );                // regressor 
     35    //Test estimation 
     36    int ndat = 100;                 // number of data records 
     37    vec Yt ( ndat );                // Store generated data 
     38    Yt.set_subvector ( 0, randn ( ord ) ); //initial values 
     39    vec rgr ( ord );                // regressor 
    4040 
    41         //print moments of the prior distribution 
    42         cout << "prior mean: " << f_thr.mean() << endl; 
    43         cout << "prior variance: " << f_thr.variance() << endl; 
     41    //print moments of the prior distribution 
     42    cout << "prior mean: " << f_thr.mean() << endl; 
     43    cout << "prior variance: " << f_thr.variance() << endl; 
    4444 
    45         // cycle over time: 
    46         for ( int t = ord; t < ndat; t++ ) { 
    47                 //Generate regressor 
    48                 for ( int j = 0; j < ( ord ); j++ ) { 
    49                         rgr ( j ) = Yt ( t - j - 1 ); 
    50                 } 
    51                 //model 
    52                 Yt ( t ) = th * rgr + sqr * NorRNG(); 
     45    // cycle over time: 
     46    for ( int t = ord; t < ndat; t++ ) { 
     47        //Generate regressor 
     48        for ( int j = 0; j < ( ord ); j++ ) { 
     49            rgr ( j ) = Yt ( t - j - 1 ); 
     50        } 
     51        //model 
     52        Yt ( t ) = th * rgr + sqr * NorRNG(); 
    5353 
    54                 Ar.bayes ( vec_1 ( Yt ( t ) ), rgr );       // Bayes rule 
     54        Ar.bayes ( vec_1 ( Yt ( t ) ), rgr );       // Bayes rule 
    5555 
    56                 // Build predictor 
    57                 mlstudent*      Pr = Ar.predictor_student ( ); 
    58                 // Test similarity of likelihoods from the Bayes rule and the predictor 
    59                 cout << "BR log-lik: " << Ar._ll(); 
    60                 cout << " , predictor ll: " <<  Pr->evallogcond ( vec_1 ( Yt ( t ) ), rgr )  << endl; 
    61                 delete Pr; 
    62         } 
    63         //print posterior moments 
    64         cout << "posterior mean: " << f_thr.mean() << endl; 
    65         cout << "posterior variance: " << f_thr.variance() << endl; 
     56        // Build predictor 
     57        mlstudent*      Pr = Ar.predictor_student ( ); 
     58        // Test similarity of likelihoods from the Bayes rule and the predictor 
     59        cout << "BR log-lik: " << Ar._ll(); 
     60        cout << " , predictor ll: " <<  Pr->evallogcond ( vec_1 ( Yt ( t ) ), rgr )  << endl; 
     61        delete Pr; 
     62    } 
     63    //print posterior moments 
     64    cout << "posterior mean: " << f_thr.mean() << endl; 
     65    cout << "posterior variance: " << f_thr.variance() << endl; 
    6666 
    67         // Test brute-froce structure estimation 
     67    // Test brute-froce structure estimation 
    6868 
    69         cout << "Structure estimation: " << endl; 
    70         cout << Ar.structure_est ( egiw ( 1, V0, nu0 ) ) << endl; 
     69    cout << "Structure estimation: " << endl; 
     70    cout << Ar.structure_est ( egiw ( 1, V0, nu0 ) ) << endl; 
    7171} 

  • library/tests/stresssuite/blas_stress.cpp

    r870 r1064  
    1010 
    1111mat matmul ( mat &A, mat &B ) { 
    12         mat C ( A.rows(), B.cols() ); 
    13         int i, j, k; 
    14         double sum; 
     12    mat C ( A.rows(), B.cols() ); 
     13    int i, j, k; 
     14    double sum; 
    1515 
    16         for ( i = 0; i < A.rows(); i++ ) { 
    17                 for ( j = 0; j < A.cols(); j++ ) { 
    18                         sum = 0.0; 
    19                         for ( k = 0; k < A.cols(); k++ ) { 
    20                                 sum += A._elem ( i, k ) * B._elem ( k, j ); 
    21                         } 
    22                         C ( i, j ) = sum; 
    23                 } 
    24         } 
    25         return C; 
     16    for ( i = 0; i < A.rows(); i++ ) { 
     17        for ( j = 0; j < A.cols(); j++ ) { 
     18            sum = 0.0; 
     19            for ( k = 0; k < A.cols(); k++ ) { 
     20                sum += A._elem ( i, k ) * B._elem ( k, j ); 
     21            } 
     22            C ( i, j ) = sum; 
     23        } 
     24    } 
     25    return C; 
    2626} 
    2727 
    2828void matmul2 ( int n,  double *A, double *B, double *C ) { 
    29         int i, j, k; 
    30         double sum; 
     29    int i, j, k; 
     30    double sum; 
    3131 
    32         for ( i = 0; i < n; i++ ) { 
    33                 for ( j = 0; j < n; j++ ) { 
    34                         sum = 0.0; 
    35                         for ( k = 0; k < n; k++ ) { 
    36                                 sum += A [ i*n+k ] * B [ k*n+j ]; 
    37                         } 
    38                         C[ i*n+j] = sum; 
    39                 } 
    40         } 
     32    for ( i = 0; i < n; i++ ) { 
     33        for ( j = 0; j < n; j++ ) { 
     34            sum = 0.0; 
     35            for ( k = 0; k < n; k++ ) { 
     36                sum += A [ i*n+k ] * B [ k*n+j ]; 
     37            } 
     38            C[ i*n+j] = sum; 
     39        } 
     40    } 
    4141//      return C; 
    4242} 
    4343 
    4444TEST ( blas_stress ) { 
    45         Real_Timer tt; 
    46         vec exec_times ( 4 ); 
    47         vec exec_times_b ( 4 ); 
    48         vec exec_times_c ( 4 ); 
     45    Real_Timer tt; 
     46    vec exec_times ( 4 ); 
     47    vec exec_times_b ( 4 ); 
     48    vec exec_times_c ( 4 ); 
    4949 
    50         mat A; 
    51         mat B; 
    52         mat C; 
     50    mat A; 
     51    mat B; 
     52    mat C; 
    5353 
    54         vec vn = "5 50 200 500"; 
    55         int n; 
     54    vec vn = "5 50 200 500"; 
     55    int n; 
    5656 
    57         for ( int i = 0; i < vn.length(); i++ ) { 
    58                 n = (int) vn ( i ); 
    59                 A = randu ( n, n ); 
    60                 B = randu ( n, n ); 
     57    for ( int i = 0; i < vn.length(); i++ ) { 
     58        n = (int) vn ( i ); 
     59        A = randu ( n, n ); 
     60        B = randu ( n, n ); 
    6161 
    62                 tt.tic(); 
    63                 for ( int ii = 0; ii < 10; ii++ ) { 
    64                         C = matmul ( A, B ); 
    65                 } 
    66                 exec_times ( i ) = tt.toc(); 
     62        tt.tic(); 
     63        for ( int ii = 0; ii < 10; ii++ ) { 
     64            C = matmul ( A, B ); 
     65        } 
     66        exec_times ( i ) = tt.toc(); 
    6767 
    68                 tt.tic(); 
    69                 for ( int ii = 0; ii < 10; ii++ ) { 
    70                         C = A * B; 
    71                 } 
    72                 exec_times_b ( i ) = tt.toc(); 
     68        tt.tic(); 
     69        for ( int ii = 0; ii < 10; ii++ ) { 
     70            C = A * B; 
     71        } 
     72        exec_times_b ( i ) = tt.toc(); 
    7373 
    74                 C = zeros ( n, n ); 
    75                 tt.tic(); 
    76                 for ( int ii = 0; ii < 10; ii++ ) { 
    77                         matmul2 ( n, A._data(), B._data(), C._data() ); 
    78                 } 
    79                 exec_times_c ( i ) = tt.toc(); 
    80         } 
    81         cout << exec_times << endl; 
    82         cout << exec_times_b << endl; 
    83         cout << exec_times_c << endl; 
     74        C = zeros ( n, n ); 
     75        tt.tic(); 
     76        for ( int ii = 0; ii < 10; ii++ ) { 
     77            matmul2 ( n, A._data(), B._data(), C._data() ); 
     78        } 
     79        exec_times_c ( i ) = tt.toc(); 
     80    } 
     81    cout << exec_times << endl; 
     82    cout << exec_times_b << endl; 
     83    cout << exec_times_c << endl; 
    8484 
    85         it_file itf ( "blas_test.it" ); 
    86         itf << Name ( "exec_times" ) << exec_times; 
    87         itf << Name ( "exec_times_b" ) << exec_times_b; 
    88         itf << Name ( "exec_times_c" ) << exec_times_c; 
     85    it_file itf ( "blas_test.it" ); 
     86    itf << Name ( "exec_times" ) << exec_times; 
     87    itf << Name ( "exec_times_b" ) << exec_times_b; 
     88    itf << Name ( "exec_times_c" ) << exec_times_c; 
    8989} 

  • library/tests/stresssuite/generator.h

    r717 r1064  
    2121class generator : public bdm::root { 
    2222public: 
    23         generator(); 
     23    generator(); 
    2424 
    25         //! Generates a matrix. Returned matrices should become 
    26         //! progressively more complicated. 
    27         virtual itpp::mat next() = 0; 
     25    //! Generates a matrix. Returned matrices should become 
     26    //! progressively more complicated. 
     27    virtual itpp::mat next() = 0; 
    2828}; 
    2929 

  • library/tests/stresssuite/kalman_QR_stress.cpp

    r722 r1064  
    1111 
    1212TEST ( kalman_QR_stress ) { 
    13         // Klaman filter 
    14         mat A, B, C, D, R, Q, P0; 
    15         vec mu0; 
    16         mat Mu0;// read from matlab 
    17         // input from Matlab 
    18         it_file fin ( "kalman_stress.it" ); 
     13    // Klaman filter 
     14    mat A, B, C, D, R, Q, P0; 
     15    vec mu0; 
     16    mat Mu0;// read from matlab 
     17    // input from Matlab 
     18    it_file fin ( "kalman_stress.it" ); 
    1919 
    20         mat Dt, XQRt, eR, eQ; 
    21         int Ndat; 
     20    mat Dt, XQRt, eR, eQ; 
     21    int Ndat; 
    2222 
    23         bool xxx = fin.seek ( "d" ); 
    24         if ( !xxx ) { 
    25                 bdm_error ( "kalman_stress.it not found" ); 
    26         } 
    27         fin >> Dt; 
    28         fin.seek ( "A" ); 
    29         fin >> A; 
    30         fin.seek ( "B" ); 
    31         fin >> B; 
    32         fin.seek ( "C" ); 
    33         fin >> C; 
    34         fin.seek ( "D" ); 
    35         fin >> D; 
    36         fin.seek ( "R" ); 
    37         fin >> R; 
    38         fin.seek ( "Q" ); 
    39         fin >> Q; 
    40         fin.seek ( "P0" ); 
    41         fin >> P0; 
    42         fin.seek ( "mu0" ); 
    43         fin >> Mu0; 
    44         mu0 = Mu0.get_col ( 0 ); 
     23    bool xxx = fin.seek ( "d" ); 
     24    if ( !xxx ) { 
     25        bdm_error ( "kalman_stress.it not found" ); 
     26    } 
     27    fin >> Dt; 
     28    fin.seek ( "A" ); 
     29    fin >> A; 
     30    fin.seek ( "B" ); 
     31    fin >> B; 
     32    fin.seek ( "C" ); 
     33    fin >> C; 
     34    fin.seek ( "D" ); 
     35    fin >> D; 
     36    fin.seek ( "R" ); 
     37    fin >> R; 
     38    fin.seek ( "Q" ); 
     39    fin >> Q; 
     40    fin.seek ( "P0" ); 
     41    fin >> P0; 
     42    fin.seek ( "mu0" ); 
     43    fin >> Mu0; 
     44    mu0 = Mu0.get_col ( 0 ); 
    4545 
    46         Ndat = Dt.cols(); 
    47         XQRt = zeros ( 5, Ndat ); 
    48         mat Xt = zeros ( 2, Ndat ); 
     46    Ndat = Dt.cols(); 
     47    XQRt = zeros ( 5, Ndat ); 
     48    mat Xt = zeros ( 2, Ndat ); 
    4949 
    5050//      cout << KF; 
    51         RV rx ( "{x }", "2" ); 
    52         RV ru ( "{u }", "1" ); 
    53         RV ry ( "{y }", "1" ); 
    54         RV rQR ( "{Q,R }", "3" ); 
    55         // 
    56         KFcondQR KF; 
    57         // KF with R unknown 
    58         KF.set_parameters ( A, B, C, D, ldmat ( R ), ldmat ( Q ) ); 
    59         KF.set_est ( mu0, ldmat ( P0 ) ); 
    60         // 
    61         Kalman<ldmat> KFtr; 
    62         // KF with R unknown 
    63         KFtr.set_parameters ( A, B, C, D, ldmat ( R ), ldmat ( Q ) ); 
    64         KFtr.set_est ( mu0, ldmat ( P0 ) ); 
     51    RV rx ( "{x }", "2" ); 
     52    RV ru ( "{u }", "1" ); 
     53    RV ry ( "{y }", "1" ); 
     54    RV rQR ( "{Q,R }", "3" ); 
     55    // 
     56    KFcondQR KF; 
     57    // KF with R unknown 
     58    KF.set_parameters ( A, B, C, D, ldmat ( R ), ldmat ( Q ) ); 
     59    KF.set_est ( mu0, ldmat ( P0 ) ); 
     60    // 
     61    Kalman<ldmat> KFtr; 
     62    // KF with R unknown 
     63    KFtr.set_parameters ( A, B, C, D, ldmat ( R ), ldmat ( Q ) ); 
     64    KFtr.set_est ( mu0, ldmat ( P0 ) ); 
    6565 
    66         mgamma evolQR; 
    67         evolQR.set_parameters ( 10.0, "1 1 1" ); //sigma = 1/10 mu 
     66    mgamma evolQR; 
     67    evolQR.set_parameters ( 10.0, "1 1 1" ); //sigma = 1/10 mu 
    6868 
    69         MPF<KFcondQR> KF_QR; 
    70         KF_QR.set_parameters ( &evolQR, &evolQR, 100 ); 
    71         evolQR.condition ( "1 1 1" ); 
    72         egamma kfinit = egamma ( "10 10 10", "1 1 1" ); 
    73         KF_QR.set_statistics ( egamma ( "10 10 10", "1 1 1" ), &KF ); 
    74         const epdf& mpost = KF_QR.posterior(); 
    75         const epdf& mposttr = KFtr.posterior(); 
     69    MPF<KFcondQR> KF_QR; 
     70    KF_QR.set_parameters ( &evolQR, &evolQR, 100 ); 
     71    evolQR.condition ( "1 1 1" ); 
     72    egamma kfinit = egamma ( "10 10 10", "1 1 1" ); 
     73    KF_QR.set_statistics ( egamma ( "10 10 10", "1 1 1" ), &KF ); 
     74    const epdf& mpost = KF_QR.posterior(); 
     75    const epdf& mposttr = KFtr.posterior(); 
    7676 
    77         XQRt.set_col ( 0, mpost.mean() ); 
    78         Xt.set_col ( 0, mposttr.mean() ); 
    79         for ( int t = 1; t < Ndat; t++ ) { 
    80                 KF_QR.bayes ( Dt.get_col ( t ) ); 
    81                 KFtr.bayes ( Dt.get_col ( t ) ); 
     77    XQRt.set_col ( 0, mpost.mean() ); 
     78    Xt.set_col ( 0, mposttr.mean() ); 
     79    for ( int t = 1; t < Ndat; t++ ) { 
     80        KF_QR.bayes ( Dt.get_col ( t ) ); 
     81        KFtr.bayes ( Dt.get_col ( t ) ); 
    8282 
    83                 XQRt.set_col ( t, mpost.mean() ); 
    84                 Xt.set_col ( t, mposttr.mean() ); 
    85         } 
     83        XQRt.set_col ( t, mpost.mean() ); 
     84        Xt.set_col ( t, mposttr.mean() ); 
     85    } 
    8686 
    87         it_file fou ( "kalman_stress_QR_res.it" ); 
    88         fou << Name ( "xqrth" ) << XQRt; 
    89         fou << Name ( "xth" ) << Xt; 
     87    it_file fou ( "kalman_stress_QR_res.it" ); 
     88    fou << Name ( "xqrth" ) << XQRt; 
     89    fou << Name ( "xth" ) << Xt; 
    9090} 

  • library/tests/stresssuite/kalman_QRexh_stress.cpp

    r721 r1064  
    1111 
    1212TEST ( kalman_QRexh_stress ) { 
    13         // Klaman filter 
    14         mat A, B, C, D, R, Q, P0; 
    15         vec mu0; 
    16         mat Mu0;// read from matlab 
    17         // input from Matlab 
    18         it_file fin ( "kalman_stress.it" ); 
     13    // Klaman filter 
     14    mat A, B, C, D, R, Q, P0; 
     15    vec mu0; 
     16    mat Mu0;// read from matlab 
     17    // input from Matlab 
     18    it_file fin ( "kalman_stress.it" ); 
    1919 
    20         mat Dt, XQRt, eR, eQ; 
    21         int Ndat; 
     20    mat Dt, XQRt, eR, eQ; 
     21    int Ndat; 
    2222 
    23         bool xxx = fin.seek ( "d" ); 
     23    bool xxx = fin.seek ( "d" ); 
    2424 
    25         if ( !xxx ) { 
    26                 bdm_error ( "kalman_stress.it not found" ); 
    27         } 
     25    if ( !xxx ) { 
     26        bdm_error ( "kalman_stress.it not found" ); 
     27    } 
    2828 
    29         fin >> Dt; 
     29    fin >> Dt; 
    3030 
    31         fin.seek ( "A" ); 
    32         fin >> A; 
    33         fin.seek ( "B" ); 
    34         fin >> B; 
    35         fin.seek ( "C" ); 
    36         fin >> C; 
    37         fin.seek ( "D" ); 
    38         fin >> D; 
    39         fin.seek ( "R" ); 
    40         fin >> R; 
    41         fin.seek ( "Q" ); 
    42         fin >> Q; 
    43         fin.seek ( "P0" ); 
    44         fin >> P0; 
    45         fin.seek ( "mu0" ); 
    46         fin >> Mu0; 
    47         mu0 = Mu0.get_col ( 0 ); 
     31    fin.seek ( "A" ); 
     32    fin >> A; 
     33    fin.seek ( "B" ); 
     34    fin >> B; 
     35    fin.seek ( "C" ); 
     36    fin >> C; 
     37    fin.seek ( "D" ); 
     38    fin >> D; 
     39    fin.seek ( "R" ); 
     40    fin >> R; 
     41    fin.seek ( "Q" ); 
     42    fin >> Q; 
     43    fin.seek ( "P0" ); 
     44    fin >> P0; 
     45    fin.seek ( "mu0" ); 
     46    fin >> Mu0; 
     47    mu0 = Mu0.get_col ( 0 ); 
    4848 
    49         vec vQ1 = "0.01:0.04:1"; 
    50         vec vQ2 = "0.01:0.04:1"; 
     49    vec vQ1 = "0.01:0.04:1"; 
     50    vec vQ2 = "0.01:0.04:1"; 
    5151 
    52         mat LL = zeros ( vQ1.length(),  vQ2.length() ); 
     52    mat LL = zeros ( vQ1.length(),  vQ2.length() ); 
    5353 
    54         Ndat = Dt.cols(); 
     54    Ndat = Dt.cols(); 
    5555 
    56         RV rx ( "{x}", "2" ); 
    57         RV ru ( "{u}", "1" ); 
    58         RV ry ( "{y}", "2" ); 
    59         // 
    60         // 
    61         Kalman<ldmat> KFtr; 
     56    RV rx ( "{x}", "2" ); 
     57    RV ru ( "{u}", "1" ); 
     58    RV ry ( "{y}", "2" ); 
     59    // 
     60    // 
     61    Kalman<ldmat> KFtr; 
    6262 
    63         for ( int i = 0; i < vQ1.length(); i++ ) { 
     63    for ( int i = 0; i < vQ1.length(); i++ ) { 
    6464 
    65                 for ( int j = 0; j < vQ2.length(); j++ ) { 
    66                         // KF with R unknown 
    67                         mat Qj = Q; 
    68                         Qj ( 0, 0 ) = vQ1 ( i ); 
    69                         Qj ( 1, 1 ) = vQ2 ( j ); 
    70                         KFtr.set_parameters ( A, B, C, D, ldmat ( R ), ldmat ( Qj ) ); 
    71                         KFtr.set_est ( mu0, ldmat ( P0 ) ); 
     65        for ( int j = 0; j < vQ2.length(); j++ ) { 
     66            // KF with R unknown 
     67            mat Qj = Q; 
     68            Qj ( 0, 0 ) = vQ1 ( i ); 
     69            Qj ( 1, 1 ) = vQ2 ( j ); 
     70            KFtr.set_parameters ( A, B, C, D, ldmat ( R ), ldmat ( Qj ) ); 
     71            KFtr.set_est ( mu0, ldmat ( P0 ) ); 
    7272 
    73                         for ( int t = 1; t < Ndat; t++ ) { 
    74                                 KFtr.bayes ( Dt.get_col ( t ) ); 
    75                                 LL ( i, j ) += KFtr._ll(); 
    76                         } 
    77                 } 
    78         } 
     73            for ( int t = 1; t < Ndat; t++ ) { 
     74                KFtr.bayes ( Dt.get_col ( t ) ); 
     75                LL ( i, j ) += KFtr._ll(); 
     76            } 
     77        } 
     78    } 
    7979 
    80         it_file fou ( "kalman_stress_QR_exh.it" ); 
     80    it_file fou ( "kalman_stress_QR_exh.it" ); 
    8181 
    82         fou << Name ( "LL" ) << LL; 
    83         fou << Name ( "Q1" ) << vQ1; 
    84         fou << Name ( "Q2" ) << vQ2; 
     82    fou << Name ( "LL" ) << LL; 
     83    fou << Name ( "Q1" ) << vQ1; 
     84    fou << Name ( "Q2" ) << vQ2; 
    8585} 

  • library/tests/stresssuite/kalman_stress.cpp

    r722 r1064  
    1010 
    1111TEST ( kalman_stress ) { 
    12         // Kalman filter 
    13         mat A, B, C, D, R, Q, P0; 
    14         vec mu0; 
    15         mat Mu0;; 
    16         // input from Matlab 
    17         it_file fin ( "kalman_stress.it" ); 
     12    // Kalman filter 
     13    mat A, B, C, D, R, Q, P0; 
     14    vec mu0; 
     15    mat Mu0;; 
     16    // input from Matlab 
     17    it_file fin ( "kalman_stress.it" ); 
    1818 
    19         mat Dt; 
    20         int Ndat; 
     19    mat Dt; 
     20    int Ndat; 
    2121 
    22         bool xxx = fin.seek ( "d" ); 
    23         if ( !xxx ) { 
    24                 bdm_error ( "kalman_stress.it not found" ); 
    25         } 
    26         fin >> Dt; 
    27         fin.seek ( "A" ); 
    28         fin >> A; 
    29         fin.seek ( "B" ); 
    30         fin >> B; 
    31         fin.seek ( "C" ); 
    32         fin >> C; 
    33         fin.seek ( "D" ); 
    34         fin >> D; 
    35         fin.seek ( "R" ); 
    36         fin >> R; 
    37         fin.seek ( "Q" ); 
    38         fin >> Q; 
    39         fin.seek ( "P0" ); 
    40         fin >> P0; 
    41         fin.seek ( "mu0" ); 
    42         fin >> Mu0; 
    43         mu0 = Mu0.get_col ( 0 ); 
     22    bool xxx = fin.seek ( "d" ); 
     23    if ( !xxx ) { 
     24        bdm_error ( "kalman_stress.it not found" ); 
     25    } 
     26    fin >> Dt; 
     27    fin.seek ( "A" ); 
     28    fin >> A; 
     29    fin.seek ( "B" ); 
     30    fin >> B; 
     31    fin.seek ( "C" ); 
     32    fin >> C; 
     33    fin.seek ( "D" ); 
     34    fin >> D; 
     35    fin.seek ( "R" ); 
     36    fin >> R; 
     37    fin.seek ( "Q" ); 
     38    fin >> Q; 
     39    fin.seek ( "P0" ); 
     40    fin >> P0; 
     41    fin.seek ( "mu0" ); 
     42    fin >> Mu0; 
     43    mu0 = Mu0.get_col ( 0 ); 
    4444 
    45         Ndat = Dt.cols(); 
    46         int dimx = A.rows(); 
     45    Ndat = Dt.cols(); 
     46    int dimx = A.rows(); 
    4747 
    48         // Prepare for Kalman filters in BDM: 
    49         RV rx ( "{x }", vec_1 ( A.cols() ) ); 
    50         RV ru ( "{u }", vec_1 ( B.cols() ) ); 
    51         RV ry ( "{y }", vec_1 ( C.rows() ) ); 
     48    // Prepare for Kalman filters in BDM: 
     49    RV rx ( "{x }", vec_1 ( A.cols() ) ); 
     50    RV ru ( "{u }", vec_1 ( B.cols() ) ); 
     51    RV ry ( "{y }", vec_1 ( C.rows() ) ); 
    5252 
    5353//      // LDMAT 
    … …  
    5959//      Xt.set_col( 0,KFep.mean() ); 
    6060 
    61         //Chol 
    62         KalmanCh KF; 
    63         KF.set_parameters ( A, B, C, D, chmat ( Q ), chmat ( R ) ); 
    64         KF.set_statistics ( mu0, chmat ( P0 ) ); //prediction! 
    65         KF.set_evalll ( false ); 
    66         KF.validate(); 
    67         const epdf& KFep = KF.posterior(); 
    68         mat Xt ( dimx, Ndat ); 
    69         Xt.set_col ( 0, KFep.mean() ); 
     61    //Chol 
     62    KalmanCh KF; 
     63    KF.set_parameters ( A, B, C, D, chmat ( Q ), chmat ( R ) ); 
     64    KF.set_statistics ( mu0, chmat ( P0 ) ); //prediction! 
     65    KF.set_evalll ( false ); 
     66    KF.validate(); 
     67    const epdf& KFep = KF.posterior(); 
     68    mat Xt ( dimx, Ndat ); 
     69    Xt.set_col ( 0, KFep.mean() ); 
    7070 
    71         // FULL 
    72         KalmanFull KF2; 
    73         KF2.set_parameters ( A, B, C, D,  Q, R ); 
    74         KF2.set_statistics ( mu0, P0 ); 
    75         KF2.set_evalll ( false ); 
    76         KF2.validate(); 
    77         mat Xt2 ( dimx, Ndat ); 
    78         Xt2.set_col ( 0, mu0 ); 
     71    // FULL 
     72    KalmanFull KF2; 
     73    KF2.set_parameters ( A, B, C, D,  Q, R ); 
     74    KF2.set_statistics ( mu0, P0 ); 
     75    KF2.set_evalll ( false ); 
     76    KF2.validate(); 
     77    mat Xt2 ( dimx, Ndat ); 
     78    Xt2.set_col ( 0, mu0 ); 
    7979 
    8080 
    81         // EKF 
    82         shared_ptr<bilinfn> fxu = new bilinfn ( A, B ); 
    83         shared_ptr<bilinfn> hxu = new bilinfn ( C, D ); 
    84         EKFCh KFE; 
    85         KFE.set_parameters ( fxu, hxu, Q, R ); 
    86         KFE.set_statistics ( mu0, chmat ( P0 ) ); 
    87         KFE.set_evalll ( false ); 
    88         KFE.validate(); 
    89         const epdf& KFEep = KFE.posterior(); 
    90         mat XtE ( dimx, Ndat ); 
    91         XtE.set_col ( 0, KFEep.mean() ); 
     81    // EKF 
     82    shared_ptr<bilinfn> fxu = new bilinfn ( A, B ); 
     83    shared_ptr<bilinfn> hxu = new bilinfn ( C, D ); 
     84    EKFCh KFE; 
     85    KFE.set_parameters ( fxu, hxu, Q, R ); 
     86    KFE.set_statistics ( mu0, chmat ( P0 ) ); 
     87    KFE.set_evalll ( false ); 
     88    KFE.validate(); 
     89    const epdf& KFEep = KFE.posterior(); 
     90    mat XtE ( dimx, Ndat ); 
     91    XtE.set_col ( 0, KFEep.mean() ); 
    9292 
    93         //test performance of each filter 
    94         Real_Timer tt; 
    95         vec exec_times ( 3 ); // KF, KF2, KFE 
     93    //test performance of each filter 
     94    Real_Timer tt; 
     95    vec exec_times ( 3 ); // KF, KF2, KFE 
    9696 
    97         vec dt; 
    98         tt.tic(); 
    99         for ( int t = 1; t < Ndat; t++ ) { 
    100                 dt = Dt.get_col ( t ); 
    101                 KF.bayes ( dt.get ( 0, C.rows() - 1 ), dt.get ( C.rows(), dt.length() - 1 ) ); 
    102                 Xt.set_col ( t, KFep.mean() ); 
    103         } 
    104         exec_times ( 0 ) = tt.toc(); 
     97    vec dt; 
     98    tt.tic(); 
     99    for ( int t = 1; t < Ndat; t++ ) { 
     100        dt = Dt.get_col ( t ); 
     101        KF.bayes ( dt.get ( 0, C.rows() - 1 ), dt.get ( C.rows(), dt.length() - 1 ) ); 
     102        Xt.set_col ( t, KFep.mean() ); 
     103    } 
     104    exec_times ( 0 ) = tt.toc(); 
    105105 
    106         tt.tic(); 
    107         for ( int t = 1; t < Ndat; t++ ) { 
    108                 dt = Dt.get_col ( t ); 
    109                 KF2.bayes ( dt.get ( 0, C.rows() - 1 ), dt.get ( C.rows(), dt.length() - 1 ) ); 
    110                 Xt2.set_col ( t, KF2.posterior().mean() ); 
    111         } 
    112         exec_times ( 1 ) = tt.toc(); 
     106    tt.tic(); 
     107    for ( int t = 1; t < Ndat; t++ ) { 
     108        dt = Dt.get_col ( t ); 
     109        KF2.bayes ( dt.get ( 0, C.rows() - 1 ), dt.get ( C.rows(), dt.length() - 1 ) ); 
     110        Xt2.set_col ( t, KF2.posterior().mean() ); 
     111    } 
     112    exec_times ( 1 ) = tt.toc(); 
    113113 
    114         tt.tic(); 
    115         for ( int t = 1; t < Ndat; t++ ) { 
    116                 dt = Dt.get_col ( t ); 
    117                 KFE.bayes ( dt.get ( 0, C.rows() - 1 ), dt.get ( C.rows(), dt.length() - 1 ) ); 
    118                 XtE.set_col ( t, KFEep.mean() ); 
    119         } 
    120         exec_times ( 2 ) = tt.toc(); 
     114    tt.tic(); 
     115    for ( int t = 1; t < Ndat; t++ ) { 
     116        dt = Dt.get_col ( t ); 
     117        KFE.bayes ( dt.get ( 0, C.rows() - 1 ), dt.get ( C.rows(), dt.length() - 1 ) ); 
     118        XtE.set_col ( t, KFEep.mean() ); 
     119    } 
     120    exec_times ( 2 ) = tt.toc(); 
    121121 
    122122 
    123         it_file fou ( "kalman_stress_res.it" ); 
    124         fou << Name ( "xth" ) << Xt; 
    125         fou << Name ( "xth2" ) << Xt2; 
    126         fou << Name ( "xthE" ) << XtE; 
    127         fou << Name ( "exec_times" ) << exec_times; 
     123    it_file fou ( "kalman_stress_res.it" ); 
     124    fou << Name ( "xth" ) << Xt; 
     125    fou << Name ( "xth2" ) << Xt2; 
     126    fou << Name ( "xthE" ) << XtE; 
     127    fou << Name ( "exec_times" ) << exec_times; 
    128128} 

  • library/tests/stresssuite/merger_2d_stress.cpp

    r886 r1064  
    1212TEST ( merger_2d_stress ) { 
    1313 
    14         RNG_randomize(); 
     14    RNG_randomize(); 
    1515 
    16         RV x ( "{x }", "1" ); 
    17         RV y ( "{y }", "1" ); 
     16    RV x ( "{x }", "1" ); 
     17    RV y ( "{y }", "1" ); 
    1818 
    19         RV xy = x; 
    20         xy.add ( y ); 
     19    RV xy = x; 
     20    xy.add ( y ); 
    2121 
    22         enorm_fsqmat_ptr f1; 
    23         f1->set_rv ( xy ); 
    24         enorm_fsqmat_ptr f2; 
    25         f2->set_rv ( xy ); 
     22    enorm_fsqmat_ptr f1; 
     23    f1->set_rv ( xy ); 
     24    enorm_fsqmat_ptr f2; 
     25    f2->set_rv ( xy ); 
    2626 
    27         mat R1 ( "0.5 0.48; 0.48 0.5" ); 
    28         mat R2 ( "0.5 0; 0 0.1" ); 
    29         f1->set_parameters ( "1 1", R1 ); 
    30         f2->set_parameters ( "1 1", mat ( "0.5 0; 0 0.1" ) ); 
     27    mat R1 ( "0.5 0.48; 0.48 0.5" ); 
     28    mat R2 ( "0.5 0; 0 0.1" ); 
     29    f1->set_parameters ( "1 1", R1 ); 
     30    f2->set_parameters ( "1 1", mat ( "0.5 0; 0 0.1" ) ); 
    3131 
    32         pdf_array A ( 2 ); 
    33         A ( 0 ) = f1; 
    34         A ( 1 ) = f2; 
     32    pdf_array A ( 2 ); 
     33    A ( 0 ) = f1; 
     34    A ( 1 ) = f2; 
    3535 
    36         int Npoints = 100; 
    37         mat x_grid ( 1, Npoints ); 
    38         x_grid.set_row ( 0, linspace ( -2.0, 4.0, Npoints ) ); 
    39         mat y_grid ( 1, Npoints ); 
    40         y_grid.set_row ( 0, linspace ( -2.0, 4.0, Npoints ) ); 
     36    int Npoints = 100; 
     37    mat x_grid ( 1, Npoints ); 
     38    x_grid.set_row ( 0, linspace ( -2.0, 4.0, Npoints ) ); 
     39    mat y_grid ( 1, Npoints ); 
     40    y_grid.set_row ( 0, linspace ( -2.0, 4.0, Npoints ) ); 
    4141 
    42         mat Grid ( 2, Npoints*Npoints ); 
    43         Grid.set_submatrix ( 0, 0, kron ( x_grid, ones ( 1, Npoints ) ) ); 
    44         Grid.set_submatrix ( 1, 0, kron ( ones ( 1, Npoints ), y_grid ) ); 
     42    mat Grid ( 2, Npoints*Npoints ); 
     43    Grid.set_submatrix ( 0, 0, kron ( x_grid, ones ( 1, Npoints ) ) ); 
     44    Grid.set_submatrix ( 1, 0, kron ( ones ( 1, Npoints ), y_grid ) ); 
    4545 
    46         merger_mix M ( A ); 
    47         enorm<fsqmat> g0; 
    48         g0.set_rv ( xy ); 
    49         g0.set_parameters ( vec ( "1 1" ), mat ( "10 0; 0 10" ) ); 
     46    merger_mix M ( A ); 
     47    enorm<fsqmat> g0; 
     48    g0.set_rv ( xy ); 
     49    g0.set_parameters ( vec ( "1 1" ), mat ( "10 0; 0 10" ) ); 
    5050 
    51         M.set_parameters ( 1 ); 
    52         M.set_method ( LOGNORMAL, 1e8 ); 
    53         M.set_support ( g0, 400 ); 
    54         M.merge(); 
     51    M.set_parameters ( 1 ); 
     52    M.set_method ( LOGNORMAL, 1e8 ); 
     53    M.set_support ( g0, 400 ); 
     54    M.merge(); 
    5555 
    56         MixEF &MM = M._Mix(); 
    57         emix* MP = MM.epredictor(); //xy 
     56    MixEF &MM = M._Mix(); 
     57    emix* MP = MM.epredictor(); //xy 
    5858 
    59         vec Res1 = M.evallog_mat ( Grid ); 
    60         mat Res2 = ( MP )->evallog_coms ( Grid ); 
     59    vec Res1 = M.evallog_mat ( Grid ); 
     60    mat Res2 = ( MP )->evallog_coms ( Grid ); 
    6161 
    62         it_file it ( "merger_2d_test.it" ); 
    63         it << Name ( "Npoints" ) << Npoints; 
    64         it << Name ( "Grid" ) << Grid; 
    65         it << Name ( "Res1" ) << Res1; 
    66         it << Name ( "Res2" ) << Res2; 
    67         it << Name ( "S1" ) << f1->evallog_mat ( Grid ); 
    68         it << Name ( "S2" ) << f2->evallog_mat ( Grid ); 
     62    it_file it ( "merger_2d_test.it" ); 
     63    it << Name ( "Npoints" ) << Npoints; 
     64    it << Name ( "Grid" ) << Grid; 
     65    it << Name ( "Res1" ) << Res1; 
     66    it << Name ( "Res2" ) << Res2; 
     67    it << Name ( "S1" ) << f1->evallog_mat ( Grid ); 
     68    it << Name ( "S2" ) << f2->evallog_mat ( Grid ); 
    6969//      cout << ( ( enorm<ldmat>* ) ( MP->_Coms ( 0 ).get() ) )->_R().to_mat() << endl; 
    7070} 

  • library/tests/stresssuite/merger_iter_stress.cpp

    r741 r1064  
    1111TEST ( merger_iter_stress ) { 
    1212 
    13         RNG_randomize(); 
     13    RNG_randomize(); 
    1414 
    15         RV x ( "{x }", "1" ); 
    16         RV y ( "{y }", "1" ); 
     15    RV x ( "{x }", "1" ); 
     16    RV y ( "{y }", "1" ); 
    1717 
    18         RV xy = x; 
    19         xy.add ( y ); 
     18    RV xy = x; 
     19    xy.add ( y ); 
    2020 
    21         enorm_fsqmat_ptr f1; 
    22         f1->set_rv ( xy ); 
    23         enorm_fsqmat_ptr f2; 
    24         f2->set_rv ( xy ); 
    25         enorm_fsqmat_ptr f3; 
    26         f3->set_rv ( y ); 
     21    enorm_fsqmat_ptr f1; 
     22    f1->set_rv ( xy ); 
     23    enorm_fsqmat_ptr f2; 
     24    f2->set_rv ( xy ); 
     25    enorm_fsqmat_ptr f3; 
     26    f3->set_rv ( y ); 
    2727 
    28         f1->set_parameters ( "4 3", mat ( "0.4 0.3; 0.3 0.4" ) ); 
    29         f1->validate(); 
    30         f2->set_parameters ( "1 3", mat ( "0.3 -0.2; -0.2 0.3" ) ); 
    31         f2->validate(); 
    32         f3->set_parameters ( "2", mat ( "0.4" ) ); 
    33         f3->validate(); 
    34          
    35         pdf_array A ( 3 ); 
    36         A ( 0 ) = f1; 
    37         A ( 1 ) = f2; 
    38         A ( 2 ) = f3; 
     28    f1->set_parameters ( "4 3", mat ( "0.4 0.3; 0.3 0.4" ) ); 
     29    f1->validate(); 
     30    f2->set_parameters ( "1 3", mat ( "0.3 -0.2; -0.2 0.3" ) ); 
     31    f2->validate(); 
     32    f3->set_parameters ( "2", mat ( "0.4" ) ); 
     33    f3->validate(); 
    3934 
    40         int Npoints = 100; 
    41         mat x_grid ( 1, Npoints ); 
    42         x_grid.set_row ( 0, linspace ( -2.0, 4.0, Npoints ) ); 
    43         mat y_grid ( 1, Npoints ); 
    44         y_grid.set_row ( 0, linspace ( -2.0, 4.0, Npoints ) ); 
     35    pdf_array A ( 3 ); 
     36    A ( 0 ) = f1; 
     37    A ( 1 ) = f2; 
     38    A ( 2 ) = f3; 
    4539 
    46         mat Grid ( 2, Npoints*Npoints ); 
    47         Grid.set_submatrix ( 0, 0, kron ( x_grid, ones ( 1, Npoints ) ) ); 
    48         Grid.set_submatrix ( 1, 0, kron ( ones ( 1, Npoints ), y_grid ) ); 
     40    int Npoints = 100; 
     41    mat x_grid ( 1, Npoints ); 
     42    x_grid.set_row ( 0, linspace ( -2.0, 4.0, Npoints ) ); 
     43    mat y_grid ( 1, Npoints ); 
     44    y_grid.set_row ( 0, linspace ( -2.0, 4.0, Npoints ) ); 
    4945 
    50         merger_mix M ( A ); 
    51         enorm<fsqmat> g0; 
    52         g0.set_rv ( xy ); 
    53         g0.set_parameters ( vec ( "4 4" ), mat ( "1 0; 0 1" ) ); 
    54         g0.validate(); 
    55          
    56         M.set_parameters ( 5 ); 
    57         M.set_method ( LOGNORMAL, 1.2 ); 
    58         M.set_support ( g0, 400 ); 
    59         M.validate(); 
    60         M.merge(); 
    61          
     46    mat Grid ( 2, Npoints*Npoints ); 
     47    Grid.set_submatrix ( 0, 0, kron ( x_grid, ones ( 1, Npoints ) ) ); 
     48    Grid.set_submatrix ( 1, 0, kron ( ones ( 1, Npoints ), y_grid ) ); 
    6249 
    63         MixEF &MM = M._Mix(); 
    64         epdf* MP = MM.epredictor();//xy 
     50    merger_mix M ( A ); 
     51    enorm<fsqmat> g0; 
     52    g0.set_rv ( xy ); 
     53    g0.set_parameters ( vec ( "4 4" ), mat ( "1 0; 0 1" ) ); 
     54    g0.validate(); 
    6555 
    66         vec Res1 = M.evallog_mat ( Grid ); 
    67         mat Res2 = ( ( emix* ) MP )->evallog_coms ( Grid ); 
     56    M.set_parameters ( 5 ); 
     57    M.set_method ( LOGNORMAL, 1.2 ); 
     58    M.set_support ( g0, 400 ); 
     59    M.validate(); 
     60    M.merge(); 
    6861 
    69         it_file it ( "merger_iter_test.it" ); 
    70         it << Name ( "Npoints" ) << Npoints; 
    71         it << Name ( "Grid" ) << Grid; 
    72         it << Name ( "Res1" ) << Res1; 
    73         it << Name ( "Res2" ) << Res2; 
     62 
     63    MixEF &MM = M._Mix(); 
     64    epdf* MP = MM.epredictor();//xy 
     65 
     66    vec Res1 = M.evallog_mat ( Grid ); 
     67    mat Res2 = ( ( emix* ) MP )->evallog_coms ( Grid ); 
     68 
     69    it_file it ( "merger_iter_test.it" ); 
     70    it << Name ( "Npoints" ) << Npoints; 
     71    it << Name ( "Grid" ) << Grid; 
     72    it << Name ( "Res1" ) << Res1; 
     73    it << Name ( "Res2" ) << Res2; 
    7474} 

  • library/tests/stresssuite/mixtures_stress.cpp

    r1009 r1064  
    1111 
    1212void disp_mix2D ( const MixEF &Mi ) { 
    13         const eprod &ep = ( ( const eprod& ) ( Mi.posterior() ) ); 
    14         int i; 
    15         mat Var ( 2, 2 ), M ( 1, 2 ); 
    16         for ( i = 0; i < ep._rv().length() - 1; i++ ) { // -1 => last one is the weight 
    17                 ( ( egiw* ) ep.factor ( i ) )->mean_mat ( M, Var ); 
    18                 cout << "Mean: " << M << endl << "Variance: " << endl << Var << endl; 
    19         } 
    20         cout << "Weights: " << ep.factor ( i )->mean() << endl; 
     13    const eprod &ep = ( ( const eprod& ) ( Mi.posterior() ) ); 
     14    int i; 
     15    mat Var ( 2, 2 ), M ( 1, 2 ); 
     16    for ( i = 0; i < ep._rv().length() - 1; i++ ) { // -1 => last one is the weight 
     17        ( ( egiw* ) ep.factor ( i ) )->mean_mat ( M, Var ); 
     18        cout << "Mean: " << M << endl << "Variance: " << endl << Var << endl; 
     19    } 
     20    cout << "Weights: " << ep.factor ( i )->mean() << endl; 
    2121} 
    2222 
    2323mat grid2D ( const MixEF &M, const vec &xb, const vec &yb, int Ngr = 100 ) { 
    24         mat PPdf ( Ngr + 1, Ngr + 1 ); 
    25         vec rgr ( 3 ); 
     24    mat PPdf ( Ngr + 1, Ngr + 1 ); 
     25    vec rgr ( 3 ); 
    2626 
    27         rgr ( 2 ) = 1; 
    28         int i = 0, j = 0; 
    29         double xstep = ( xb ( 1 ) - xb ( 0 ) ) / Ngr; 
    30         double ystep = ( yb ( 1 ) - yb ( 0 ) ) / Ngr; 
     27    rgr ( 2 ) = 1; 
     28    int i = 0, j = 0; 
     29    double xstep = ( xb ( 1 ) - xb ( 0 ) ) / Ngr; 
     30    double ystep = ( yb ( 1 ) - yb ( 0 ) ) / Ngr; 
    3131 
    32         for ( double x = xb ( 0 ); x <= xb ( 1 ); x += xstep, i++ ) { 
    33                 rgr ( 0 ) = x; 
    34                 j = 0; 
    35                 for ( double y = yb ( 0 ); y <= yb ( 1 ); y += ystep, j++ ) { 
    36                         rgr ( 1 ) = y; 
    37                         PPdf ( i, j ) = exp ( M.logpred ( rgr, empty_vec ) ); 
    38                 } 
    39         } 
    40         return PPdf; 
     32    for ( double x = xb ( 0 ); x <= xb ( 1 ); x += xstep, i++ ) { 
     33        rgr ( 0 ) = x; 
     34        j = 0; 
     35        for ( double y = yb ( 0 ); y <= yb ( 1 ); y += ystep, j++ ) { 
     36            rgr ( 1 ) = y; 
     37            PPdf ( i, j ) = exp ( M.logpred ( rgr, empty_vec ) ); 
     38        } 
     39    } 
     40    return PPdf; 
    4141} 
    4242 
    4343TEST ( mixtures_stress ) { 
    44         RV x ( "{x }", "2" ); 
     44    RV x ( "{x }", "2" ); 
    4545 
    4646 
    47         cout << "Test for estimation of a 2D Gaussian mixture" << endl; 
    48         /////////////////////////////////// 
    49         // Setup model 
     47    cout << "Test for estimation of a 2D Gaussian mixture" << endl; 
     48    /////////////////////////////////// 
     49    // Setup model 
    5050 
    51         vec m1 ( "-2 -2" ); 
    52         vec m2 ( "2 2" ); 
     51    vec m1 ( "-2 -2" ); 
     52    vec m2 ( "2 2" ); 
    5353 
    54         fsqmat V1 ( mat ( "1 0.5; 0.5 1" ) ); 
    55         fsqmat V2 ( mat ( "2 -0.1; -0.1 2" ) ); 
     54    fsqmat V1 ( mat ( "1 0.5; 0.5 1" ) ); 
     55    fsqmat V2 ( mat ( "2 -0.1; -0.1 2" ) ); 
    5656 
    57         enorm_fsqmat_ptr C1; 
    58         C1->set_rv ( x ); 
    59         C1->set_parameters ( m1, V1 ); 
    60         enorm_fsqmat_ptr C2; 
    61         C2->set_rv ( x ); 
    62         C2->set_parameters ( m2, V2 ); 
     57    enorm_fsqmat_ptr C1; 
     58    C1->set_rv ( x ); 
     59    C1->set_parameters ( m1, V1 ); 
     60    enorm_fsqmat_ptr C2; 
     61    C2->set_rv ( x ); 
     62    C2->set_parameters ( m2, V2 ); 
    6363 
    64         epdf_array Sim ( 2 ); 
    65         Sim ( 0 ) = C1; 
    66         Sim ( 1 ) = C2; 
    67         emix Simul; 
    68         Simul.set_rv ( x ); 
    69         Simul._w() = "0.5 0.6"; 
    70         Simul._Coms() = Sim; 
    71         Simul.validate(); 
     64    epdf_array Sim ( 2 ); 
     65    Sim ( 0 ) = C1; 
     66    Sim ( 1 ) = C2; 
     67    emix Simul; 
     68    Simul.set_rv ( x ); 
     69    Simul._w() = "0.5 0.6"; 
     70    Simul._Coms() = Sim; 
     71    Simul.validate(); 
    7272 
    73         // Sample parameters 
    74         int ndat = 100; 
    75         mat Smp = Simul.sample_mat ( ndat ); 
     73    // Sample parameters 
     74    int ndat = 100; 
     75    mat Smp = Simul.sample_mat ( ndat ); 
    7676 
    77         cout << "Simulated means: " << m1 << " and " << m2 << endl; 
    78         cout << "Simulated covariances: " << endl << V1 << " and " << endl << V2 << endl; 
     77    cout << "Simulated means: " << m1 << " and " << m2 << endl; 
     78    cout << "Simulated covariances: " << endl << V1 << " and " << endl << V2 << endl; 
    7979 
    80         ////////////////////////////// 
    81         // Prepare estimator manually 
     80    ////////////////////////////// 
     81    // Prepare estimator manually 
    8282 
    83         // Initial values of components 
    84         mat V0 = 1e-4 * eye ( 3 ); 
    85         V0 ( 0, 0 ) *= 100; 
    86         V0 ( 1, 1 ) *= 100; 
    87         mat Vp = "0 0 1; 0 0 1; 1 1 0"; 
    88         Vp *= 5 * 1e-5; 
     83    // Initial values of components 
     84    mat V0 = 1e-4 * eye ( 3 ); 
     85    V0 ( 0, 0 ) *= 100; 
     86    V0 ( 1, 1 ) *= 100; 
     87    mat Vp = "0 0 1; 0 0 1; 1 1 0"; 
     88    Vp *= 5 * 1e-5; 
    8989 
    90         ARX M1; 
    91         M1.set_statistics ( 2, V0 - Vp, 8 ); 
    92         M1.validate(); 
    93         ARX M2; 
    94         M2.set_statistics ( 2, V0 + Vp, 8 ); 
    95         M2.validate(); 
     90    ARX M1; 
     91    M1.set_statistics ( 2, V0 - Vp, 8 ); 
     92    M1.validate(); 
     93    ARX M2; 
     94    M2.set_statistics ( 2, V0 + Vp, 8 ); 
     95    M2.validate(); 
    9696 
    97         // Build mixture model 
    98         Array<BMEF*> A ( 2 ); 
    99         A ( 0 ) = &M1; 
    100         A ( 1 ) = &M2; 
    101         MixEF Post ( A, "0.5 0.5" ); 
    102         cout << "Initial mixture:" << endl; 
    103         disp_mix2D ( Post ); 
     97    // Build mixture model 
     98    Array<BMEF*> A ( 2 ); 
     99    A ( 0 ) = &M1; 
     100    A ( 1 ) = &M2; 
     101    MixEF Post ( A, "0.5 0.5" ); 
     102    cout << "Initial mixture:" << endl; 
     103    disp_mix2D ( Post ); 
    104104 
    105         // Add ones for constant coefficients 
    106         mat Data = concat_vertical ( Smp, ones ( 1, Smp.cols() ) ); 
    107         Post.bayes ( Data , empty_vec ); 
     105    // Add ones for constant coefficients 
     106    mat Data = concat_vertical ( Smp, ones ( 1, Smp.cols() ) ); 
     107    Post.bayes ( Data , empty_vec ); 
    108108 
    109         cout << "Posterior mixture:" << endl; 
    110         disp_mix2D ( Post ); 
     109    cout << "Posterior mixture:" << endl; 
     110    disp_mix2D ( Post ); 
    111111 
    112         //TEST random initialization 
    113         RNG_randomize(); 
    114         ARX Aflat; 
    115         Aflat.set_statistics ( 2, V0, 7 ); 
    116         Aflat.validate(); 
    117         MixEF RND; 
    118         RND.init ( &Aflat, Data, 10 ); // already initialized! 
    119         cout << endl << "== Randomly initialized mixture ==" << endl; 
    120         cout << endl << "== INIT ==" << endl; 
    121         disp_mix2D ( RND ); 
    122         mat PPdf_I = grid2D ( RND, vec_2 ( -5.0, 5.0 ), vec_2 ( -5.0, 5.0 ) );; 
     112    //TEST random initialization 
     113    RNG_randomize(); 
     114    ARX Aflat; 
     115    Aflat.set_statistics ( 2, V0, 7 ); 
     116    Aflat.validate(); 
     117    MixEF RND; 
     118    RND.init ( &Aflat, Data, 10 ); // already initialized! 
     119    cout << endl << "== Randomly initialized mixture ==" << endl; 
     120    cout << endl << "== INIT ==" << endl; 
     121    disp_mix2D ( RND ); 
     122    mat PPdf_I = grid2D ( RND, vec_2 ( -5.0, 5.0 ), vec_2 ( -5.0, 5.0 ) );; 
    123123 
    124         RND.bayes ( Data, empty_vec ); 
    125         cout << endl << "== BAYES ==" << endl; 
    126         disp_mix2D ( RND ); 
     124    RND.bayes ( Data, empty_vec ); 
     125    cout << endl << "== BAYES ==" << endl; 
     126    disp_mix2D ( RND ); 
    127127 
    128128 
    129         it_file of ( "mixef_test.it", true ); 
    130         of << Name ( "Smp" ) << Smp; 
    131         of << Name ( "PPdf" ) << grid2D ( Post, vec_2 ( -5.0, 5.0 ), vec_2 ( -5.0, 5.0 ) ); 
    132         of << Name ( "PPdf_I" ) << PPdf_I; 
    133         of << Name ( "PPdf_RND" ) << grid2D ( RND, vec_2 ( -5.0, 5.0 ), vec_2 ( -5.0, 5.0 ) ); 
     129    it_file of ( "mixef_test.it", true ); 
     130    of << Name ( "Smp" ) << Smp; 
     131    of << Name ( "PPdf" ) << grid2D ( Post, vec_2 ( -5.0, 5.0 ), vec_2 ( -5.0, 5.0 ) ); 
     132    of << Name ( "PPdf_I" ) << PPdf_I; 
     133    of << Name ( "PPdf_RND" ) << grid2D ( RND, vec_2 ( -5.0, 5.0 ), vec_2 ( -5.0, 5.0 ) ); 
    134134 
    135         cout << endl << "Variables written to file mixef_test.it:" << endl; 
    136         cout << "Smp  : data sampled from the simulator " << endl; 
    137         cout << "PPdf : posterior density on a grid <-5,5><-5,5>" << endl; 
    138         cout << "PPdf_RND : posterior density of randomly initilaized mixture on a grid <-5,5><-5,5>" << endl; 
     135    cout << endl << "Variables written to file mixef_test.it:" << endl; 
     136    cout << "Smp  : data sampled from the simulator " << endl; 
     137    cout << "PPdf : posterior density on a grid <-5,5><-5,5>" << endl; 
     138    cout << "PPdf_RND : posterior density of randomly initilaized mixture on a grid <-5,5><-5,5>" << endl; 
    139139} 

  • library/tests/stresssuite/particle_stress.cpp

    r887 r1064  
    1111 
    1212TEST ( particle_stress ) { 
    13         RV x ( "1" ); 
    14         RV xm = x; 
    15         xm.t_plus ( -1 ); 
    16         const 
    17         RV y ( "2" ); 
     13    RV x ( "1" ); 
     14    RV xm = x; 
     15    xm.t_plus ( -1 ); 
     16    const 
     17    RV y ( "2" ); 
    1818 
    19         mat A = "1"; 
    20         vec vR = "1"; 
    21         ldmat R ( vR ); 
     19    mat A = "1"; 
     20    vec vR = "1"; 
     21    ldmat R ( vR ); 
    2222 
    23         eEmp emp; 
    24         euni eun; 
    25         eun.set_parameters ( "0", "1" ); 
    26         emp.set_statistics ( ones ( 10 ), eun ); 
    27         vec &v = emp._w(); 
    28         Array<vec> &S = emp._samples(); 
     23    eEmp emp; 
     24    euni eun; 
     25    eun.set_parameters ( "0", "1" ); 
     26    emp.set_statistics ( ones ( 10 ), eun ); 
     27    vec &v = emp._w(); 
     28    Array<vec> &S = emp._samples(); 
    2929 
    30         for ( int i = 0; i < 10; i++ ) { 
    31                 v ( i ) = exp ( -0.5 * sum ( pow ( S ( i ) - 1, 2.0 ) ) * 10 ); 
    32         } 
    33         v /= sum ( v ); 
     30    for ( int i = 0; i < 10; i++ ) { 
     31        v ( i ) = exp ( -0.5 * sum ( pow ( S ( i ) - 1, 2.0 ) ) * 10 ); 
     32    } 
     33    v /= sum ( v ); 
    3434 
    35         cout << "p:" << S << endl; 
    36         cout << "w:" << v << endl; 
     35    cout << "p:" << S << endl; 
     36    cout << "w:" << v << endl; 
    3737 
    38         ivec ind; 
    39         resample ( v, ind ); 
     38    ivec ind; 
     39    resample ( v, ind ); 
    4040 
    41         cout << ind << endl; 
     41    cout << ind << endl; 
    4242} 

  • library/tests/stresssuite/resample_stress.cpp

    r887 r1064  
    1212TEST ( resample_stress ) { 
    1313 
    14         RV x ( "1" ); 
    15         RV xm = x; 
    16         xm.t_plus ( -1 ); 
    17         const 
    18         RV y ( "2" ); 
     14    RV x ( "1" ); 
     15    RV xm = x; 
     16    xm.t_plus ( -1 ); 
     17    const 
     18    RV y ( "2" ); 
    1919 
    20         mat A = "1"; 
    21         vec vR = "1"; 
    22         ldmat R ( vR ); 
     20    mat A = "1"; 
     21    vec vR = "1"; 
     22    ldmat R ( vR ); 
    2323 
    24         eEmp emp; 
    25         euni eun; 
    26         eun.set_parameters ( "0", "1" ); 
    27         emp.set_statistics ( ones ( 10 ), eun ); 
    28         vec &v = emp._w(); 
    29         Array<vec> &S = emp._samples(); 
     24    eEmp emp; 
     25    euni eun; 
     26    eun.set_parameters ( "0", "1" ); 
     27    emp.set_statistics ( ones ( 10 ), eun ); 
     28    vec &v = emp._w(); 
     29    Array<vec> &S = emp._samples(); 
    3030 
    31         for ( int i = 0; i < 10; i++ ) { 
    32                 v ( i ) = exp ( -0.5 * sum ( pow ( S ( i ) - 1, 2.0 ) ) * 10 ); 
    33         } 
    34         v /= sum ( v ); 
     31    for ( int i = 0; i < 10; i++ ) { 
     32        v ( i ) = exp ( -0.5 * sum ( pow ( S ( i ) - 1, 2.0 ) ) * 10 ); 
     33    } 
     34    v /= sum ( v ); 
    3535 
    36         cout << "p:" << S << endl; 
    37         cout << "w:" << v << endl; 
     36    cout << "p:" << S << endl; 
     37    cout << "w:" << v << endl; 
    3838 
    39         ivec ind; 
    40         resample ( v,ind ); 
     39    ivec ind; 
     40    resample ( v,ind ); 
    4141 
    42         cout << ind << endl; 
     42    cout << ind << endl; 
    4343} 

  • library/tests/stresssuite/size_generator.cpp

    r722 r1064  
    66 
    77void size_generator::from_setting ( const Setting &set ) { 
    8         UI::get ( sz, set, "size", UI::optional ); 
    9         UI::get ( step, set, "step", UI::optional ); 
     8    UI::get ( sz, set, "size", UI::optional ); 
     9    UI::get ( step, set, "step", UI::optional ); 
    1010} 
    1111 
    1212mat size_generator::next() { 
    13         mat A0 = randu ( sz, sz ); 
    14         mat A = A0 * A0.T(); 
    15         sz *= step; 
    16         return A; 
     13    mat A0 = randu ( sz, sz ); 
     14    mat A = A0 * A0.T(); 
     15    sz *= step; 
     16    return A; 
    1717} 

  • library/tests/stresssuite/size_generator.h

    r717 r1064  
    2121class size_generator : public generator { 
    2222private: 
    23         int step; 
    24         int sz; 
     23    int step; 
     24    int sz; 
    2525 
    2626public: 
    27         size_generator() : step ( 7 ), sz ( step ) { } 
     27    size_generator() : step ( 7 ), sz ( step ) { } 
    2828 
    29         itpp::mat next(); 
     29    itpp::mat next(); 
    3030 
    31         //! Load from structure with elements: 
    32         //!  \code 
    33         //! { size = 7; // size (rows == cols) of the first generated matrix 
    34         //!   step = 7; // how many times to increase the generated matrix in every iteration 
    35         //! } 
    36         //! \endcode 
    37         //! Both elements are optional, with defaults as shown. 
    38         void from_setting ( const Setting &set ); 
     31    //! Load from structure with elements: 
     32    //!  \code 
     33    //! { size = 7; // size (rows == cols) of the first generated matrix 
     34    //!   step = 7; // how many times to increase the generated matrix in every iteration 
     35    //! } 
     36    //! \endcode 
     37    //! Both elements are optional, with defaults as shown. 
     38    void from_setting ( const Setting &set ); 
    3939}; 
    4040 

  • library/tests/stresssuite/square_mat_prep.cpp

    r717 r1064  
    2222 
    2323int main ( int argc, char const *argv[] ) { 
    24         RNG_randomize(); 
     24    RNG_randomize(); 
    2525 
    26         bool unknown = false; 
    27         int update_next = 0; // 1 generator file, 2 agenda file, 3 agenda length 
    28         const char **param = argv + 1; 
    29         while ( *param && !unknown ) { 
    30                 if ( update_next ) { 
    31                         if ( update_next == 1 ) { 
    32                                 generator_file_name = *param; 
    33                         } else if ( update_next == 2 ) { 
    34                                 agenda_file_name = *param; 
    35                         } else { 
    36                                 int length = atoi ( *param ); 
    37                                 if ( length > 0 ) { 
    38                                         agenda_length = length; 
    39                                 } else { 
    40                                         cerr << "invalid agenda length value ignored" << endl; 
    41                                 } 
    42                         } 
     26    bool unknown = false; 
     27    int update_next = 0; // 1 generator file, 2 agenda file, 3 agenda length 
     28    const char **param = argv + 1; 
     29    while ( *param && !unknown ) { 
     30        if ( update_next ) { 
     31            if ( update_next == 1 ) { 
     32                generator_file_name = *param; 
     33            } else if ( update_next == 2 ) { 
     34                agenda_file_name = *param; 
     35            } else { 
     36                int length = atoi ( *param ); 
     37                if ( length > 0 ) { 
     38                    agenda_length = length; 
     39                } else { 
     40                    cerr << "invalid agenda length value ignored" << endl; 
     41                } 
     42            } 
    4343 
    44                         update_next = 0; 
    45                 } else { 
    46                         if ( !strcmp ( *param, "-a" ) ) { 
    47                                 update_next = 2; 
    48                         } else if ( !strcmp ( *param, "-g" ) ) { 
    49                                 update_next = 1; 
    50                         } else if ( !strcmp ( *param, "-l" ) ) { 
    51                                 update_next = 3; 
    52                         } else { 
    53                                 unknown = true; 
    54                         } 
    55                 } 
     44            update_next = 0; 
     45        } else { 
     46            if ( !strcmp ( *param, "-a" ) ) { 
     47                update_next = 2; 
     48            } else if ( !strcmp ( *param, "-g" ) ) { 
     49                update_next = 1; 
     50            } else if ( !strcmp ( *param, "-l" ) ) { 
     51                update_next = 3; 
     52            } else { 
     53                unknown = true; 
     54            } 
     55        } 
    5656 
    57                 ++param; 
    58         } 
     57        ++param; 
     58    } 
    5959 
    60         if ( unknown || update_next ) { 
    61                 cerr << "usage: " << argv[0] << " [ -g generator.cfg ] [ -a agenda_output.cfg ] [ -l agenda_length ]" << endl; 
    62         } else { 
    63                 Array<shared_ptr<square_mat_point> > mag ( agenda_length ); 
     60    if ( unknown || update_next ) { 
     61        cerr << "usage: " << argv[0] << " [ -g generator.cfg ] [ -a agenda_output.cfg ] [ -l agenda_length ]" << endl; 
     62    } else { 
     63        Array<shared_ptr<square_mat_point> > mag ( agenda_length ); 
    6464 
    65                 UIFile gspec ( generator_file_name ); 
    66                 shared_ptr<generator> gen ( UI::build<generator> ( gspec, "generator", UI::compulsory ) ); 
    67                 for ( int i = 0; i < agenda_length; ++i ) { 
    68                         mat m = gen->next(); 
    69                         square_mat_point *p = new square_mat_point(); 
    70                         p->set_parameters ( m, randu ( m.rows() ), randu() ); 
    71                         mag ( i ) = p; 
    72                 } 
     65        UIFile gspec ( generator_file_name ); 
     66        shared_ptr<generator> gen ( UI::build<generator> ( gspec, "generator", UI::compulsory ) ); 
     67        for ( int i = 0; i < agenda_length; ++i ) { 
     68            mat m = gen->next(); 
     69            square_mat_point *p = new square_mat_point(); 
     70            p->set_parameters ( m, randu ( m.rows() ), randu() ); 
     71            mag ( i ) = p; 
     72        } 
    7373 
    74                 UIFile fag; 
    75                 UI::save ( mag, fag, "agenda" ); 
    76                 fag.save ( agenda_file_name ); 
    77         } 
     74        UIFile fag; 
     75        UI::save ( mag, fag, "agenda" ); 
     76        fag.save ( agenda_file_name ); 
     77    } 
    7878} 

  • library/tests/stresssuite/square_mat_stress.cpp

    r766 r1064  
    3737                  const itpp::mat &actual, double tolerance, 
    3838                  TestDetails const& details ) { 
    39         if ( !AreClose ( expected, actual, tolerance ) ) { 
    40                 MemoryOutStream stream; 
    41                 stream << "failed at " << expected.rows() 
    42                 << " x " << expected.cols(); 
    43  
    44                 results.OnTestFailure ( details, stream.GetText() ); 
    45         } 
     39    if ( !AreClose ( expected, actual, tolerance ) ) { 
     40        MemoryOutStream stream; 
     41        stream << "failed at " << expected.rows() 
     42               << " x " << expected.cols(); 
     43 
     44        results.OnTestFailure ( details, stream.GetText() ); 
     45    } 
    4646} 
    4747 
    … …  
    5252template<typename TMatrix> 
    5353void test_matrix ( int index, square_mat_point *point ) { 
    54         Real_Timer tt; 
    55  
    56         cout << "agenda[" << index << "]:" << endl; 
    57         mat A = point->get_matrix(); 
    58         int sz = A.rows(); 
    59         CHECK_EQUAL ( A.cols(), sz ); 
    60  
    61         tt.tic(); 
    62         TMatrix sq_mat ( A ); 
    63         double elapsed = tt.toc(); 
    64         cout << "ctor(" << sz << " x " << sz << "): " << elapsed << " s" << endl; 
    65  
    66         tt.tic(); 
    67         mat res = sq_mat.to_mat(); 
    68         elapsed = tt.toc(); 
    69  
    70         if ( !fast ) { 
    71                 CHECK_CLOSE ( A, res, epsilon ); 
    72         } 
    73  
    74         cout << "to_mat: " << elapsed << " s" << endl; 
    75  
    76         vec v = point->get_vector(); 
    77         double w = point->get_scalar(); 
    78         TMatrix sq_mat2 = sq_mat; 
    79  
    80         tt.tic(); 
    81         sq_mat2.opupdt ( v, w ); 
    82         elapsed = tt.toc(); 
    83  
    84         if ( !fast ) { 
    85                 mat expA = A + w * outer_product ( v, v ); 
    86                 CHECK_CLOSE ( expA, sq_mat2.to_mat(), epsilon ); 
    87         } 
    88  
    89         cout << "opupdt: " << elapsed << " s" << endl; 
    90  
    91         TMatrix invmat ( sz ); 
    92  
    93         tt.tic(); 
    94         sq_mat.inv ( invmat ); 
    95         elapsed = tt.toc(); 
    96  
    97         mat invA; 
    98         if ( !fast ) { 
    99                 invA = inv ( A ); 
    100                 CHECK_CLOSE ( invA, invmat.to_mat(), epsilon ); 
    101         } 
    102  
    103         cout << "inv: " << elapsed << " s" << endl; 
    104  
    105         tt.tic(); 
    106         double ld = sq_mat.logdet(); 
    107         elapsed = tt.toc(); 
    108  
    109         if ( !fast ) { 
    110                 double d = det ( A ); 
    111                 CHECK_CLOSE ( log ( d ), ld, epsilon ); 
    112         } 
    113  
    114         cout << "logdet: " << elapsed << " s" << endl; 
    115  
    116         tt.tic(); 
    117         double q = sq_mat.qform ( ones ( sz ) ); 
    118         elapsed = tt.toc(); 
    119  
    120         if ( !fast ) { 
    121                 CHECK_CLOSE ( sumsum ( A ), q, epsilon ); 
    122         } 
    123  
    124         cout << "qform(1): " << elapsed << " s" << endl; 
    125  
    126         tt.tic(); 
    127         q = sq_mat.qform ( v ); 
    128         elapsed = tt.toc(); 
    129  
    130         if ( !fast ) { 
    131                 double r = ( A * v ) * v; 
    132                 CHECK_CLOSE ( r, q, epsilon ); 
    133         } 
    134  
    135         cout << "qform(v): " << elapsed << " s" << endl; 
    136  
    137         tt.tic(); 
    138         q = sq_mat.invqform ( v ); 
    139         elapsed = tt.toc(); 
    140  
    141         if ( !fast ) { 
    142                 double r = ( invA * v ) * v; 
    143                 CHECK_CLOSE ( r, q, epsilon ); 
    144         } 
    145  
    146         cout << "invqform: " << elapsed << " s" << endl; 
    147  
    148         TMatrix twice = sq_mat; 
    149  
    150         tt.tic(); 
    151         twice += sq_mat; 
    152         elapsed = tt.toc(); 
    153  
    154         if ( !fast ) { 
    155                 res = 2 * A; 
    156                 CHECK_CLOSE ( res, twice.to_mat(), epsilon ); 
    157         } 
    158  
    159         cout << "+=: " << elapsed << " s" << endl; 
    160  
    161         sq_mat2 = sq_mat; 
    162  
    163         tt.tic(); 
    164         sq_mat2.mult_sym ( A ); 
    165         elapsed = tt.toc(); 
    166  
    167         if ( !fast ) { 
    168                 res = ( A * A ) * A.T(); 
    169                 CHECK_CLOSE ( res, sq_mat2.to_mat(), epsilon ); 
    170         } 
    171  
    172         cout << "^2: " << elapsed << " s" << endl; 
     54    Real_Timer tt; 
     55 
     56    cout << "agenda[" << index << "]:" << endl; 
     57    mat A = point->get_matrix(); 
     58    int sz = A.rows(); 
     59    CHECK_EQUAL ( A.cols(), sz ); 
     60 
     61    tt.tic(); 
     62    TMatrix sq_mat ( A ); 
     63    double elapsed = tt.toc(); 
     64    cout << "ctor(" << sz << " x " << sz << "): " << elapsed << " s" << endl; 
     65 
     66    tt.tic(); 
     67    mat res = sq_mat.to_mat(); 
     68    elapsed = tt.toc(); 
     69 
     70    if ( !fast ) { 
     71        CHECK_CLOSE ( A, res, epsilon ); 
     72    } 
     73 
     74    cout << "to_mat: " << elapsed << " s" << endl; 
     75 
     76    vec v = point->get_vector(); 
     77    double w = point->get_scalar(); 
     78    TMatrix sq_mat2 = sq_mat; 
     79 
     80    tt.tic(); 
     81    sq_mat2.opupdt ( v, w ); 
     82    elapsed = tt.toc(); 
     83 
     84    if ( !fast ) { 
     85        mat expA = A + w * outer_product ( v, v ); 
     86        CHECK_CLOSE ( expA, sq_mat2.to_mat(), epsilon ); 
     87    } 
     88 
     89    cout << "opupdt: " << elapsed << " s" << endl; 
     90 
     91    TMatrix invmat ( sz ); 
     92 
     93    tt.tic(); 
     94    sq_mat.inv ( invmat ); 
     95    elapsed = tt.toc(); 
     96 
     97    mat invA; 
     98    if ( !fast ) { 
     99        invA = inv ( A ); 
     100        CHECK_CLOSE ( invA, invmat.to_mat(), epsilon ); 
     101    } 
     102 
     103    cout << "inv: " << elapsed << " s" << endl; 
     104 
     105    tt.tic(); 
     106    double ld = sq_mat.logdet(); 
     107    elapsed = tt.toc(); 
     108 
     109    if ( !fast ) { 
     110        double d = det ( A ); 
     111        CHECK_CLOSE ( log ( d ), ld, epsilon ); 
     112    } 
     113 
     114    cout << "logdet: " << elapsed << " s" << endl; 
     115 
     116    tt.tic(); 
     117    double q = sq_mat.qform ( ones ( sz ) ); 
     118    elapsed = tt.toc(); 
     119 
     120    if ( !fast ) { 
     121        CHECK_CLOSE ( sumsum ( A ), q, epsilon ); 
     122    } 
     123 
     124    cout << "qform(1): " << elapsed << " s" << endl; 
     125 
     126    tt.tic(); 
     127    q = sq_mat.qform ( v ); 
     128    elapsed = tt.toc(); 
     129 
     130    if ( !fast ) { 
     131        double r = ( A * v ) * v; 
     132        CHECK_CLOSE ( r, q, epsilon ); 
     133    } 
     134 
     135    cout << "qform(v): " << elapsed << " s" << endl; 
     136 
     137    tt.tic(); 
     138    q = sq_mat.invqform ( v ); 
     139    elapsed = tt.toc(); 
     140 
     141    if ( !fast ) { 
     142        double r = ( invA * v ) * v; 
     143        CHECK_CLOSE ( r, q, epsilon ); 
     144    } 
     145 
     146    cout << "invqform: " << elapsed << " s" << endl; 
     147 
     148    TMatrix twice = sq_mat; 
     149 
     150    tt.tic(); 
     151    twice += sq_mat; 
     152    elapsed = tt.toc(); 
     153 
     154    if ( !fast ) { 
     155        res = 2 * A; 
     156        CHECK_CLOSE ( res, twice.to_mat(), epsilon ); 
     157    } 
     158 
     159    cout << "+=: " << elapsed << " s" << endl; 
     160 
     161    sq_mat2 = sq_mat; 
     162 
     163    tt.tic(); 
     164    sq_mat2.mult_sym ( A ); 
     165    elapsed = tt.toc(); 
     166 
     167    if ( !fast ) { 
     168        res = ( A * A ) * A.T(); 
     169        CHECK_CLOSE ( res, sq_mat2.to_mat(), epsilon ); 
     170    } 
     171 
     172    cout << "^2: " << elapsed << " s" << endl; 
    173173} 
    174174 
    175175void test_agenda ( FTestMatrix test ) { 
    176         UIFile fag ( agenda_file_name ); 
    177         Array<shared_ptr<square_mat_point> > mag; 
    178         UI::get ( mag, fag, "agenda", UI::compulsory ); 
    179         int sz = mag.size(); 
    180         CHECK ( sz > 0 ); 
    181         for ( int i = 0; i < sz; ++i ) { 
    182                 test ( i, mag ( i ).get() ); 
    183         } 
     176    UIFile fag ( agenda_file_name ); 
     177    Array<shared_ptr<square_mat_point> > mag; 
     178    UI::get ( mag, fag, "agenda", UI::compulsory ); 
     179    int sz = mag.size(); 
     180    CHECK ( sz > 0 ); 
     181    for ( int i = 0; i < sz; ++i ) { 
     182        test ( i, mag ( i ).get() ); 
     183    } 
    184184} 
    185185 
    186186SUITE ( ldmat ) { 
    187         TEST ( agenda ) { 
    188                 test_agenda ( test_matrix<ldmat> ); 
    189         } 
     187    TEST ( agenda ) { 
     188        test_agenda ( test_matrix<ldmat> ); 
     189    } 
    190190} 
    191191 
    192192SUITE ( fsqmat ) { 
    193         TEST ( agenda ) { 
    194                 test_agenda ( test_matrix<fsqmat> ); 
    195         } 
     193    TEST ( agenda ) { 
     194        test_agenda ( test_matrix<fsqmat> ); 
     195    } 
    196196} 
    197197 
    198198SUITE ( chmat ) { 
    199         TEST ( agenda ) { 
    200                 test_agenda ( test_matrix<chmat> ); 
    201         } 
     199    TEST ( agenda ) { 
     200        test_agenda ( test_matrix<chmat> ); 
     201    } 
    202202} 
    203203 
    204204int main ( int argc, char const *argv[] ) { 
    205         bool unknown = false; 
    206         int update_next = 0; // 1 suite, 2 epsilon, 3 agenda file 
    207         const char *suite = "ldmat"; 
    208         const char **param = argv + 1; 
    209         while ( *param && !unknown ) { 
    210                 if ( update_next ) { 
    211                         if ( update_next == 1 ) { 
    212                                 suite = *param; 
    213                         } else if ( update_next == 2 ) { 
    214                                 double eps = atof ( *param ); 
    215                                 if ( eps > 0 ) { 
    216                                         epsilon = eps; 
    217                                 } else { 
    218                                         cerr << "invalid epsilon value ignored" << endl; 
    219                                 } 
    220                         } else { 
    221                                 agenda_file_name = *param; 
    222                         } 
    223  
    224                         update_next = 0; 
    225                 } else { 
    226                         if ( !strcmp ( *param, "-a" ) ) { 
    227                                 update_next = 3; 
    228                         } else if ( !strcmp ( *param, "-c" ) ) { 
    229                                 update_next = 1; 
    230                         } else if ( !strcmp ( *param, "-e" ) ) { 
    231                                 update_next = 2; 
    232                         } else if ( !strcmp ( *param, "-f" ) ) { 
    233                                 fast = true; 
    234                         } else { 
    235                                 unknown = true; 
    236                         } 
    237                 } 
    238  
    239                 ++param; 
    240         } 
    241  
    242         if ( unknown || update_next ) { 
    243                 cerr << "usage: " << argv[0] << " [ -f ] [ -e epsilon ] [ -a agenda_input.cfg ] [ -c class ]" << endl; 
    244         } else { 
    245                 UnitTest::TestReporterStdout reporter; 
    246                 UnitTest::TestRunner runner ( reporter ); 
    247                 return runner.RunTestsIf ( UnitTest::Test::GetTestList(), 
    248                                            suite, 
    249                                            UnitTest::True(), 
    250                                            0 ); 
    251         } 
    252 } 
     205    bool unknown = false; 
     206    int update_next = 0; // 1 suite, 2 epsilon, 3 agenda file 
     207    const char *suite = "ldmat"; 
     208    const char **param = argv + 1; 
     209    while ( *param && !unknown ) { 
     210        if ( update_next ) { 
     211            if ( update_next == 1 ) { 
     212                suite = *param; 
     213            } else if ( update_next == 2 ) { 
     214                double eps = atof ( *param ); 
     215                if ( eps > 0 ) { 
     216                    epsilon = eps; 
     217                } else { 
     218                    cerr << "invalid epsilon value ignored" << endl; 
     219                } 
     220            } else { 
     221                agenda_file_name = *param; 
     222            } 
     223 
     224            update_next = 0; 
     225        } else { 
     226            if ( !strcmp ( *param, "-a" ) ) { 
     227                update_next = 3; 
     228            } else if ( !strcmp ( *param, "-c" ) ) { 
     229                update_next = 1; 
     230            } else if ( !strcmp ( *param, "-e" ) ) { 
     231                update_next = 2; 
     232            } else if ( !strcmp ( *param, "-f" ) ) { 
     233                fast = true; 
     234            } else { 
     235                unknown = true; 
     236            } 
     237        } 
     238 
     239        ++param; 
     240    } 
     241 
     242    if ( unknown || update_next ) { 
     243        cerr << "usage: " << argv[0] << " [ -f ] [ -e epsilon ] [ -a agenda_input.cfg ] [ -c class ]" << endl; 
     244    } else { 
     245        UnitTest::TestReporterStdout reporter; 
     246        UnitTest::TestRunner runner ( reporter ); 
     247        return runner.RunTestsIf ( UnitTest::Test::GetTestList(), 
     248                                   suite, 
     249                                   UnitTest::True(), 
     250                                   0 ); 
     251    } 
     252} 

  • library/tests/stresssuite/stresssuite.cpp

    r722 r1064  
    1010 
    1111int main ( int argc, char const *argv[] ) { 
    12         if ( argc > 1 ) { 
    13                 if ( !strcmp ( argv[1], "print" ) ) 
    14                         print_test_list(); 
    15                 else 
    16                         pick_selected_tests ( argc, argv ); 
    17         } else { 
    18                 cout << endl << "STRESSSUITE - a program covering all BDM stress tests." << endl << endl 
    19                      << argv[0] << " ....................................... run all stress tests" << endl 
    20                      << argv[0] << " particular_test_1 particular_test_2 ... run selected stress tests" << endl 
    21                      << argv[0] << " print ................................. print all the implemented stress tests" << endl; 
    22         } 
     12    if ( argc > 1 ) { 
     13        if ( !strcmp ( argv[1], "print" ) ) 
     14            print_test_list(); 
     15        else 
     16            pick_selected_tests ( argc, argv ); 
     17    } else { 
     18        cout << endl << "STRESSSUITE - a program covering all BDM stress tests." << endl << endl 
     19             << argv[0] << " ....................................... run all stress tests" << endl 
     20             << argv[0] << " particular_test_1 particular_test_2 ... run selected stress tests" << endl 
     21             << argv[0] << " print ................................. print all the implemented stress tests" << endl; 
     22    } 
    2323 
    24         return run_selected_tests(); 
     24    return run_selected_tests(); 
    2525}