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

astyle applied all over the library

Files:
1 modified

Legend:

Unmodified
Added
Removed

  • library/bdm/stat/exp_family.cpp

    r1063 r1064  
    1515 
    1616void BMEF::bayes ( const vec &yt, const vec &cond ) { 
    17         this->bayes_weighted ( yt, cond, 1.0 ); 
     17    this->bayes_weighted ( yt, cond, 1.0 ); 
    1818}; 
    1919 
    2020void egiw::set_parameters ( int dimx0, ldmat V0, double nu0 ) { 
    21         dimx = dimx0; 
    22         nPsi = V0.rows()-dimx; 
    23         V = V0; 
    24         if ( nu0 < 0 ) { 
    25                 nu = 0.1 + nPsi + 2 * dimx + 2; // +2 assures finite expected value of R 
    26                 // terms before that are sufficient for finite normalization 
    27         } else { 
    28                 nu = nu0; 
    29         } 
     21    dimx = dimx0; 
     22    nPsi = V0.rows()-dimx; 
     23    V = V0; 
     24    if ( nu0 < 0 ) { 
     25        nu = 0.1 + nPsi + 2 * dimx + 2; // +2 assures finite expected value of R 
     26        // terms before that are sufficient for finite normalization 
     27    } else { 
     28        nu = nu0; 
     29    } 
    3030} 
    3131 
    3232vec egiw::sample() const { 
    33         mat M; 
    34         chmat R; 
    35         sample_mat ( M, R ); 
    36  
    37         return concat ( cvectorize ( M ), cvectorize ( R.to_mat() ) ); 
     33    mat M; 
     34    chmat R; 
     35    sample_mat ( M, R ); 
     36 
     37    return concat ( cvectorize ( M ), cvectorize ( R.to_mat() ) ); 
    3838} 
    3939 
    4040mat egiw::sample_mat ( int n ) const { 
    41         // TODO - correct approach - convert to product of norm * Wishart 
    42         mat M; 
    43         ldmat Vz; 
    44         ldmat Lam; 
    45         factorize ( M, Vz, Lam ); 
    46  
    47         chmat ChLam ( Lam.to_mat() ); 
    48         chmat iChLam; 
    49         ChLam.inv ( iChLam ); 
    50  
    51         eWishartCh Omega; //inverse Wishart, result is R, 
    52         Omega.set_parameters ( iChLam, nu - 2*nPsi - dimx ); // 2*nPsi is there to match numercial simulations - check if analytically correct 
    53         Omega.validate();        
    54  
    55         mat OmChi; 
    56         mat Z ( M.rows(), M.cols() ); 
    57  
    58         mat Mi; 
    59         mat RChiT; 
    60         mat tmp ( dimension(), n ); 
    61         M=M.T();// ugly hack == decide what to do with M. 
    62         for ( int i = 0; i < n; i++ ) { 
    63                 OmChi = Omega.sample_mat(); 
    64                 RChiT = inv ( OmChi ); 
    65                 Z = randn ( M.rows(), M.cols() ); 
    66                 Mi = M + RChiT * Z * inv ( Vz._L().T() * diag ( sqrt ( Vz._D() ) ) ); 
    67  
    68                 tmp.set_col ( i, concat ( cvectorize ( Mi ), cvectorize ( RChiT*RChiT.T() ) ) ); 
    69         } 
    70         return tmp; 
     41    // TODO - correct approach - convert to product of norm * Wishart 
     42    mat M; 
     43    ldmat Vz; 
     44    ldmat Lam; 
     45    factorize ( M, Vz, Lam ); 
     46 
     47    chmat ChLam ( Lam.to_mat() ); 
     48    chmat iChLam; 
     49    ChLam.inv ( iChLam ); 
     50 
     51    eWishartCh Omega; //inverse Wishart, result is R, 
     52    Omega.set_parameters ( iChLam, nu - 2*nPsi - dimx ); // 2*nPsi is there to match numercial simulations - check if analytically correct 
     53    Omega.validate(); 
     54 
     55    mat OmChi; 
     56    mat Z ( M.rows(), M.cols() ); 
     57 
     58    mat Mi; 
     59    mat RChiT; 
     60    mat tmp ( dimension(), n ); 
     61    M=M.T();// ugly hack == decide what to do with M. 
     62    for ( int i = 0; i < n; i++ ) { 
     63        OmChi = Omega.sample_mat(); 
     64        RChiT = inv ( OmChi ); 
     65        Z = randn ( M.rows(), M.cols() ); 
     66        Mi = M + RChiT * Z * inv ( Vz._L().T() * diag ( sqrt ( Vz._D() ) ) ); 
     67 
     68        tmp.set_col ( i, concat ( cvectorize ( Mi ), cvectorize ( RChiT*RChiT.T() ) ) ); 
     69    } 
     70    return tmp; 
    7171} 
    7272 
    7373void egiw::sample_mat ( mat &Mi, chmat &Ri ) const { 
    7474 
    75         // TODO - correct approach - convert to product of norm * Wishart 
    76         mat M; 
    77         ldmat Vz; 
    78         ldmat Lam; 
    79         factorize ( M, Vz, Lam ); 
    80  
    81         chmat Ch; 
    82         Ch.setCh ( Lam._L() *diag ( sqrt ( Lam._D() ) ) ); 
    83         chmat iCh; 
    84         Ch.inv ( iCh ); 
    85  
    86         eWishartCh Omega; //inverse Wishart, result is R, 
    87         Omega.set_parameters ( iCh, nu - 2*nPsi - dimx ); // 2*nPsi is there to match numercial simulations - check if analytically correct 
    88         Omega.validate(); 
    89  
    90         chmat Omi; 
    91         Omi.setCh ( Omega.sample_mat() ); 
    92  
    93         if (M._datasize()>0){ 
    94                 mat Z = randn ( M.rows(), M.cols() ); 
    95                 Mi = M + Omi._Ch() * Z * inv ( Vz._L() * diag ( sqrt ( Vz._D() ) ) ); 
    96         } 
    97         Omi.inv ( Ri ); 
     75    // TODO - correct approach - convert to product of norm * Wishart 
     76    mat M; 
     77    ldmat Vz; 
     78    ldmat Lam; 
     79    factorize ( M, Vz, Lam ); 
     80 
     81    chmat Ch; 
     82    Ch.setCh ( Lam._L() *diag ( sqrt ( Lam._D() ) ) ); 
     83    chmat iCh; 
     84    Ch.inv ( iCh ); 
     85 
     86    eWishartCh Omega; //inverse Wishart, result is R, 
     87    Omega.set_parameters ( iCh, nu - 2*nPsi - dimx ); // 2*nPsi is there to match numercial simulations - check if analytically correct 
     88    Omega.validate(); 
     89 
     90    chmat Omi; 
     91    Omi.setCh ( Omega.sample_mat() ); 
     92 
     93    if (M._datasize()>0) { 
     94        mat Z = randn ( M.rows(), M.cols() ); 
     95        Mi = M + Omi._Ch() * Z * inv ( Vz._L() * diag ( sqrt ( Vz._D() ) ) ); 
     96    } 
     97    Omi.inv ( Ri ); 
    9898} 
    9999 
    100100double egiw::evallog_nn ( const vec &val ) const { 
    101         bdm_assert_debug(val.length()==dimx*(nPsi+dimx),"Incorrect cond in egiw::evallog_nn" ); 
    102          
    103         int vend = val.length() - 1; 
    104  
    105         if ( dimx == 1 ) { //same as the following, just quicker. 
    106                 double r = val ( vend ); //last entry! 
    107                 if ( r < 0 ) return -inf; 
    108                 vec Psi ( nPsi + dimx ); 
    109                 Psi ( 0 ) = -1.0; 
    110                 Psi.set_subvector ( 1, val ( 0, vend - 1 ) ); // fill the rest 
    111  
    112                 double Vq = V.qform ( Psi ); 
    113                 return -0.5* ( nu*log ( r ) + Vq / r ); 
    114         } else { 
    115                 mat Tmp; 
    116                 if (nPsi>0){ 
    117                         mat Th = reshape ( val ( 0, nPsi * dimx - 1 ), nPsi, dimx ); 
    118                         Tmp = concat_vertical ( -eye ( dimx ), Th ); 
    119                 } else { 
    120                         Tmp = -eye(dimx); 
    121                 } 
    122                 fsqmat R ( reshape ( val ( nPsi*dimx, vend ), dimx, dimx ) ); 
    123                 double ldetR = R.logdet(); 
    124                 if ( !std::isfinite(ldetR) ) return -inf; 
    125                 fsqmat iR ( dimx ); 
    126                 R.inv ( iR ); 
    127  
    128                 return -0.5* ( nu*ldetR + trace ( iR.to_mat() *Tmp.T() *V.to_mat() *Tmp ) ); 
    129         } 
     101    bdm_assert_debug(val.length()==dimx*(nPsi+dimx),"Incorrect cond in egiw::evallog_nn" ); 
     102 
     103    int vend = val.length() - 1; 
     104 
     105    if ( dimx == 1 ) { //same as the following, just quicker. 
     106        double r = val ( vend ); //last entry! 
     107        if ( r < 0 ) return -inf; 
     108        vec Psi ( nPsi + dimx ); 
     109        Psi ( 0 ) = -1.0; 
     110        Psi.set_subvector ( 1, val ( 0, vend - 1 ) ); // fill the rest 
     111 
     112        double Vq = V.qform ( Psi ); 
     113        return -0.5* ( nu*log ( r ) + Vq / r ); 
     114    } else { 
     115        mat Tmp; 
     116        if (nPsi>0) { 
     117            mat Th = reshape ( val ( 0, nPsi * dimx - 1 ), nPsi, dimx ); 
     118            Tmp = concat_vertical ( -eye ( dimx ), Th ); 
     119        } else { 
     120            Tmp = -eye(dimx); 
     121        } 
     122        fsqmat R ( reshape ( val ( nPsi*dimx, vend ), dimx, dimx ) ); 
     123        double ldetR = R.logdet(); 
     124        if ( !std::isfinite(ldetR) ) return -inf; 
     125        fsqmat iR ( dimx ); 
     126        R.inv ( iR ); 
     127 
     128        return -0.5* ( nu*ldetR + trace ( iR.to_mat() *Tmp.T() *V.to_mat() *Tmp ) ); 
     129    } 
    130130} 
    131131 
    132132double egiw::lognc() const { 
    133         const vec& D = V._D(); 
    134  
    135         double m = nu - nPsi - dimx - 1; 
     133    const vec& D = V._D(); 
     134 
     135    double m = nu - nPsi - dimx - 1; 
    136136#define log2  0.693147180559945286226763983 
    137137#define logpi 1.144729885849400163877476189 
    … …  
    139139#define Inf std::numeric_limits<double>::infinity() 
    140140 
    141         double nkG = 0.5 * dimx * ( -nPsi * log2pi + sum ( log ( D.mid ( dimx, nPsi ) ) ) ); 
    142         // temporary for lgamma in Wishart 
    143         double lg = 0; 
    144         for ( int i = 0; i < dimx; i++ ) { 
    145                 lg += lgamma ( 0.5 * ( m - i ) ); 
    146         } 
    147  
    148         double nkW = 0.5 * ( m * sum ( log ( D ( 0, dimx - 1 ) ) ) ) \ 
    149                      - 0.5 * dimx * ( m * log2 + 0.5 * ( dimx - 1 ) * log2pi )  - lg; 
     141    double nkG = 0.5 * dimx * ( -nPsi * log2pi + sum ( log ( D.mid ( dimx, nPsi ) ) ) ); 
     142    // temporary for lgamma in Wishart 
     143    double lg = 0; 
     144    for ( int i = 0; i < dimx; i++ ) { 
     145        lg += lgamma ( 0.5 * ( m - i ) ); 
     146    } 
     147 
     148    double nkW = 0.5 * ( m * sum ( log ( D ( 0, dimx - 1 ) ) ) ) \ 
     149                 - 0.5 * dimx * ( m * log2 + 0.5 * ( dimx - 1 ) * log2pi )  - lg; 
    150150 
    151151//      bdm_assert_debug ( ( ( -nkG - nkW ) > -Inf ) && ( ( -nkG - nkW ) < Inf ), "ARX improper" ); 
    152         if ( -nkG - nkW == Inf ) { 
    153                 cout << "??" << endl; 
    154         } 
    155         return -nkG - nkW; 
     152    if ( -nkG - nkW == Inf ) { 
     153        cout << "??" << endl; 
     154    } 
     155    return -nkG - nkW; 
    156156} 
    157157 
    158158vec egiw::est_theta() const { 
    159         if ( dimx == 1 ) { 
    160                 const mat &L = V._L(); 
    161                 int end = L.rows() - 1; 
    162  
    163                 mat iLsub = ltuinv ( L ( dimx, end, dimx, end ) ); 
    164  
    165                 vec L0 = L.get_col ( 0 ); 
    166  
    167                 return iLsub * L0 ( 1, end ); 
    168         } else { 
    169                 bdm_error ( "ERROR: est_theta() not implemented for dimx>1" ); 
    170                 return vec(); 
    171         } 
     159    if ( dimx == 1 ) { 
     160        const mat &L = V._L(); 
     161        int end = L.rows() - 1; 
     162 
     163        mat iLsub = ltuinv ( L ( dimx, end, dimx, end ) ); 
     164 
     165        vec L0 = L.get_col ( 0 ); 
     166 
     167        return iLsub * L0 ( 1, end ); 
     168    } else { 
     169        bdm_error ( "ERROR: est_theta() not implemented for dimx>1" ); 
     170        return vec(); 
     171    } 
    172172} 
    173173 
    174174void egiw::factorize ( mat &M, ldmat &Vz, ldmat &Lam ) const { 
    175         const mat &L = V._L(); 
    176         const vec &D = V._D(); 
    177         int end = L.rows() - 1; 
    178  
    179         Lam = ldmat ( L ( 0, dimx - 1, 0, dimx - 1 ), D ( 0, dimx - 1 ) );  //exp val of R 
    180  
    181         if (dimx<=end){ 
    182                 Vz = ldmat ( L ( dimx, end, dimx, end ), D ( dimx, end ) ); 
    183                 mat iLsub = ltuinv ( Vz._L() ); 
    184                 // set mean value 
    185                 mat Lpsi = L ( dimx, end, 0, dimx - 1 ); 
    186                 M = iLsub * Lpsi; 
    187         } 
    188 /*      if ( 0 ) { // test with Peterka 
    189                 mat VF = V.to_mat(); 
    190                 mat Vf = VF ( 0, dimx - 1, 0, dimx - 1 ); 
    191                 mat Vzf = VF ( dimx, end, 0, dimx - 1 ); 
    192                 mat VZ = VF ( dimx, end, dimx, end ); 
    193  
    194                 mat Lam2 = Vf - Vzf.T() * inv ( VZ ) * Vzf; 
    195         }*/ 
     175    const mat &L = V._L(); 
     176    const vec &D = V._D(); 
     177    int end = L.rows() - 1; 
     178 
     179    Lam = ldmat ( L ( 0, dimx - 1, 0, dimx - 1 ), D ( 0, dimx - 1 ) );  //exp val of R 
     180 
     181    if (dimx<=end) { 
     182        Vz = ldmat ( L ( dimx, end, dimx, end ), D ( dimx, end ) ); 
     183        mat iLsub = ltuinv ( Vz._L() ); 
     184        // set mean value 
     185        mat Lpsi = L ( dimx, end, 0, dimx - 1 ); 
     186        M = iLsub * Lpsi; 
     187    } 
     188    /*  if ( 0 ) { // test with Peterka 
     189                mat VF = V.to_mat(); 
     190                mat Vf = VF ( 0, dimx - 1, 0, dimx - 1 ); 
     191                mat Vzf = VF ( dimx, end, 0, dimx - 1 ); 
     192                mat VZ = VF ( dimx, end, dimx, end ); 
     193 
     194                mat Lam2 = Vf - Vzf.T() * inv ( VZ ) * Vzf; 
     195        }*/ 
    196196} 
    197197 
    198198ldmat egiw::est_theta_cov() const { 
    199         if ( dimx == 1 ) { 
    200                 const mat &L = V._L(); 
    201                 const vec &D = V._D(); 
    202                 int end = D.length() - 1; 
    203  
    204                 mat Lsub = L ( 1, end, 1, end ); 
     199    if ( dimx == 1 ) { 
     200        const mat &L = V._L(); 
     201        const vec &D = V._D(); 
     202        int end = D.length() - 1; 
     203 
     204        mat Lsub = L ( 1, end, 1, end ); 
    205205//              mat Dsub = diag ( D ( 1, end ) ); 
    206206 
    207                 ldmat LD ( inv ( Lsub ).T(), 1.0 / D ( 1, end ) ); 
    208                 return LD; 
    209  
    210         } else { 
    211                 bdm_error ( "ERROR: est_theta_cov() not implemented for dimx>1" ); 
    212                 return ldmat(); 
    213         } 
     207        ldmat LD ( inv ( Lsub ).T(), 1.0 / D ( 1, end ) ); 
     208        return LD; 
     209 
     210    } else { 
     211        bdm_error ( "ERROR: est_theta_cov() not implemented for dimx>1" ); 
     212        return ldmat(); 
     213    } 
    214214 
    215215} 
    … …  
    217217vec egiw::mean() const { 
    218218 
    219         if ( dimx == 1 ) { 
    220                 const vec &D = V._D(); 
    221                 int end = D.length() - 1; 
    222  
    223                 vec m ( dim ); 
    224                 m.set_subvector ( 0, est_theta() ); 
    225                 m ( end ) = D ( 0 ) / ( nu - nPsi - 2 * dimx - 2 ); 
    226                 return m; 
    227         } else { 
    228                 mat M; 
    229                 mat R; 
    230                 mean_mat ( M, R ); 
    231                 return concat ( cvectorize ( M ), cvectorize ( R ) ); 
    232         } 
     219    if ( dimx == 1 ) { 
     220        const vec &D = V._D(); 
     221        int end = D.length() - 1; 
     222 
     223        vec m ( dim ); 
     224        m.set_subvector ( 0, est_theta() ); 
     225        m ( end ) = D ( 0 ) / ( nu - nPsi - 2 * dimx - 2 ); 
     226        return m; 
     227    } else { 
     228        mat M; 
     229        mat R; 
     230        mean_mat ( M, R ); 
     231        return concat ( cvectorize ( M ), cvectorize ( R ) ); 
     232    } 
    233233 
    234234} 
    235235 
    236236vec egiw::variance() const { 
    237         int l = V.rows(); 
    238         // cut out rest of lower-right part of V 
    239         // invert it 
    240         ldmat itmp; 
    241         if (dimx<l){ 
    242                 const ldmat tmp ( V, linspace ( dimx, l - 1 ) ); 
    243                 tmp.inv ( itmp ); 
    244         }  
    245         // following Wikipedia notation 
    246         // m=nu-nPsi-dimx-1, p=dimx 
    247         double mp1p = nu - nPsi - 2 * dimx; // m-p+1 
    248         double mp1m = mp1p - 2;     // m-p-1 
    249          
    250         if ( dimx == 1 ) { 
    251                 double cove = V._D() ( 0 ) / mp1m ; 
    252                  
    253                 vec var ( l ); 
    254                 var.set_subvector ( 0, diag ( itmp.to_mat() ) *cove ); 
    255                 var ( l - 1 ) = cove * cove / ( mp1m - 2 ); 
    256                 return var; 
    257         } else { 
    258                 ldmat Vll ( V, linspace ( 0, dimx - 1 ) ); // top-left part of V 
    259                 mat Y = Vll.to_mat(); 
    260                 mat varY ( Y.rows(), Y.cols() ); 
    261                  
    262                 double denom = ( mp1p - 1 ) * mp1m * mp1m * ( mp1m - 2 );         // (m-p)(m-p-1)^2(m-p-3) 
    263                  
    264                 int i, j; 
    265                 for ( i = 0; i < Y.rows(); i++ ) { 
    266                         for ( j = 0; j < Y.cols(); j++ ) { 
    267                                 varY ( i, j ) = ( mp1p * Y ( i, j ) * Y ( i, j ) + mp1m * Y ( i, i ) * Y ( j, j ) ) / denom; 
    268                         } 
    269                 } 
    270                 vec mean_dR = diag ( Y ) / mp1m; // corresponds to cove 
    271                 vec var_th = diag ( itmp.to_mat() ); 
    272                 vec var_Th ( mean_dR.length() *var_th.length() ); 
    273                 // diagonal of diag(mean_dR) \kron diag(var_th) 
    274                 for ( int i = 0; i < mean_dR.length(); i++ ) { 
    275                         var_Th.set_subvector ( i*var_th.length(), var_th*mean_dR ( i ) ); 
    276                 } 
    277                  
    278                 return concat ( var_Th, cvectorize ( varY ) ); 
    279         } 
     237    int l = V.rows(); 
     238    // cut out rest of lower-right part of V 
     239    // invert it 
     240    ldmat itmp; 
     241    if (dimx<l) { 
     242        const ldmat tmp ( V, linspace ( dimx, l - 1 ) ); 
     243        tmp.inv ( itmp ); 
     244    } 
     245    // following Wikipedia notation 
     246    // m=nu-nPsi-dimx-1, p=dimx 
     247    double mp1p = nu - nPsi - 2 * dimx; // m-p+1 
     248    double mp1m = mp1p - 2;     // m-p-1 
     249 
     250    if ( dimx == 1 ) { 
     251        double cove = V._D() ( 0 ) / mp1m ; 
     252 
     253        vec var ( l ); 
     254        var.set_subvector ( 0, diag ( itmp.to_mat() ) *cove ); 
     255        var ( l - 1 ) = cove * cove / ( mp1m - 2 ); 
     256        return var; 
     257    } else { 
     258        ldmat Vll ( V, linspace ( 0, dimx - 1 ) ); // top-left part of V 
     259        mat Y = Vll.to_mat(); 
     260        mat varY ( Y.rows(), Y.cols() ); 
     261 
     262        double denom = ( mp1p - 1 ) * mp1m * mp1m * ( mp1m - 2 );         // (m-p)(m-p-1)^2(m-p-3) 
     263 
     264        int i, j; 
     265        for ( i = 0; i < Y.rows(); i++ ) { 
     266            for ( j = 0; j < Y.cols(); j++ ) { 
     267                varY ( i, j ) = ( mp1p * Y ( i, j ) * Y ( i, j ) + mp1m * Y ( i, i ) * Y ( j, j ) ) / denom; 
     268            } 
     269        } 
     270        vec mean_dR = diag ( Y ) / mp1m; // corresponds to cove 
     271        vec var_th = diag ( itmp.to_mat() ); 
     272        vec var_Th ( mean_dR.length() *var_th.length() ); 
     273        // diagonal of diag(mean_dR) \kron diag(var_th) 
     274        for ( int i = 0; i < mean_dR.length(); i++ ) { 
     275            var_Th.set_subvector ( i*var_th.length(), var_th*mean_dR ( i ) ); 
     276        } 
     277 
     278        return concat ( var_Th, cvectorize ( varY ) ); 
     279    } 
    280280} 
    281281 
    282282void egiw::mean_mat ( mat &M, mat&R ) const { 
    283         const mat &L = V._L(); 
    284         const vec &D = V._D(); 
    285         int end = L.rows() - 1; 
    286  
    287         ldmat ldR ( L ( 0, dimx - 1, 0, dimx - 1 ), D ( 0, dimx - 1 ) / ( nu - nPsi - 2*dimx - 2 ) ); //exp val of R 
    288  
    289         // set mean value 
    290         if (dimx<=end){ 
    291                 mat iLsub = ltuinv ( L ( dimx, end, dimx, end ) ); 
    292                 mat Lpsi = L ( dimx, end, 0, dimx - 1 ); 
    293                 M = iLsub * Lpsi; 
    294         } 
    295         R = ldR.to_mat()  ; 
     283    const mat &L = V._L(); 
     284    const vec &D = V._D(); 
     285    int end = L.rows() - 1; 
     286 
     287    ldmat ldR ( L ( 0, dimx - 1, 0, dimx - 1 ), D ( 0, dimx - 1 ) / ( nu - nPsi - 2*dimx - 2 ) ); //exp val of R 
     288 
     289    // set mean value 
     290    if (dimx<=end) { 
     291        mat iLsub = ltuinv ( L ( dimx, end, dimx, end ) ); 
     292        mat Lpsi = L ( dimx, end, 0, dimx - 1 ); 
     293        M = iLsub * Lpsi; 
     294    } 
     295    R = ldR.to_mat()  ; 
    296296} 
    297297 
    298298void egiw::log_register ( bdm::logger& L, const string& prefix ) { 
    299         epdf::log_register ( L, prefix ); 
    300         if ( log_level[logvartheta] ) {  
    301                 int th_dim = dim - dimx*dimx; // dimension - dimension of cov 
    302                 L.add_vector( log_level, logvartheta, RV ( th_dim ), prefix );  
    303         } 
     299    epdf::log_register ( L, prefix ); 
     300    if ( log_level[logvartheta] ) { 
     301        int th_dim = dim - dimx*dimx; // dimension - dimension of cov 
     302        L.add_vector( log_level, logvartheta, RV ( th_dim ), prefix ); 
     303    } 
    304304} 
    305305 
    306306void egiw::log_write() const { 
    307         epdf::log_write(); 
    308         if ( log_level[logvartheta] ) {  
    309                 mat M; 
    310                 ldmat Lam; 
    311                 ldmat Vz; 
    312                 factorize ( M, Vz, Lam ); 
    313                 if( log_level[logvartheta] ) 
    314                         log_level.store( logvartheta, cvectorize ( est_theta_cov().to_mat() ) ); 
    315         } 
     307    epdf::log_write(); 
     308    if ( log_level[logvartheta] ) { 
     309        mat M; 
     310        ldmat Lam; 
     311        ldmat Vz; 
     312        factorize ( M, Vz, Lam ); 
     313        if( log_level[logvartheta] ) 
     314            log_level.store( logvartheta, cvectorize ( est_theta_cov().to_mat() ) ); 
     315    } 
    316316} 
    317317 
    318318void egiw::from_setting ( const Setting &set ) { 
    319                 epdf::from_setting ( set ); 
    320                 UI::get ( dimx, set, "dimx", UI::compulsory ); 
    321                 if ( !UI::get ( nu, set, "nu", UI::optional ) ) { 
    322                         nu = -1; 
    323                 } 
    324                 mat Vful; 
    325                 if (!UI::get(V, set, "V", UI::optional)){ 
    326                         if ( !UI::get ( Vful, set, "V", UI::optional ) ) { 
    327                                 vec dV; 
    328                                 UI::get ( dV, set, "dV", UI::compulsory ); 
    329                                 set_parameters ( dimx, ldmat ( dV ), nu ); 
    330  
    331                         } else { 
    332                                 set_parameters ( dimx, Vful, nu ); 
    333                         } 
    334                 } 
    335         } 
     319    epdf::from_setting ( set ); 
     320    UI::get ( dimx, set, "dimx", UI::compulsory ); 
     321    if ( !UI::get ( nu, set, "nu", UI::optional ) ) { 
     322        nu = -1; 
     323    } 
     324    mat Vful; 
     325    if (!UI::get(V, set, "V", UI::optional)) { 
     326        if ( !UI::get ( Vful, set, "V", UI::optional ) ) { 
     327            vec dV; 
     328            UI::get ( dV, set, "dV", UI::compulsory ); 
     329            set_parameters ( dimx, ldmat ( dV ), nu ); 
     330 
     331        } else { 
     332            set_parameters ( dimx, Vful, nu ); 
     333        } 
     334    } 
     335} 
    336336 
    337337void egiw::to_setting ( Setting& set ) const { 
    338                 epdf::to_setting ( set ); 
    339                 UI::save ( dimx, set, "dimx" ); 
    340                 UI::save ( V, set, "V" ); 
    341                 UI::save ( nu, set, "nu" ); 
    342         }; 
     338    epdf::to_setting ( set ); 
     339    UI::save ( dimx, set, "dimx" ); 
     340    UI::save ( V, set, "V" ); 
     341    UI::save ( nu, set, "nu" ); 
     342}; 
    343343 
    344344void egiw::validate() { 
    345         eEF::validate(); 
    346         nPsi = V.rows() - dimx; 
    347         dim = dimx * ( dimx + nPsi );    
    348          
    349         if ( nu < 0 ) { 
    350                 nu = 0.1 + nPsi + 2 * dimx + 2; // +2 assures finite expected value of R 
    351                 // terms before that are sufficient for finite normalization 
    352         } 
    353          
    354             // check sizes, rvs etc. 
    355                 // also check if RV are meaningful!!! 
    356                 // meaningful =  rv for theta  and rv for r are split! 
    357         } 
     345    eEF::validate(); 
     346    nPsi = V.rows() - dimx; 
     347    dim = dimx * ( dimx + nPsi ); 
     348 
     349    if ( nu < 0 ) { 
     350        nu = 0.1 + nPsi + 2 * dimx + 2; // +2 assures finite expected value of R 
     351        // terms before that are sufficient for finite normalization 
     352    } 
     353 
     354    // check sizes, rvs etc. 
     355    // also check if RV are meaningful!!! 
     356    // meaningful =  rv for theta  and rv for r are split! 
     357} 
    358358void multiBM::bayes ( const vec &yt, const vec &cond ) { 
    359         if ( frg < 1.0 ) { 
    360                 beta *= frg; 
    361                 last_lognc = est.lognc(); 
    362         } 
    363         beta += yt; 
    364         if ( evalll ) { 
    365                 ll = est.lognc() - last_lognc; 
    366         } 
     359    if ( frg < 1.0 ) { 
     360        beta *= frg; 
     361        last_lognc = est.lognc(); 
     362    } 
     363    beta += yt; 
     364    if ( evalll ) { 
     365        ll = est.lognc() - last_lognc; 
     366    } 
    367367} 
    368368 
    369369double multiBM::logpred ( const vec &yt ) const { 
    370         eDirich pred ( est ); 
    371         vec &beta = pred._beta(); 
    372  
    373         double lll; 
    374         if ( frg < 1.0 ) { 
    375                 beta *= frg; 
    376                 lll = pred.lognc(); 
    377         } else if ( evalll ) { 
    378                 lll = last_lognc; 
    379         } else { 
    380                 lll = pred.lognc(); 
    381         } 
    382  
    383         beta += yt; 
    384         return pred.lognc() - lll; 
     370    eDirich pred ( est ); 
     371    vec &beta = pred._beta(); 
     372 
     373    double lll; 
     374    if ( frg < 1.0 ) { 
     375        beta *= frg; 
     376        lll = pred.lognc(); 
     377    } else if ( evalll ) { 
     378        lll = last_lognc; 
     379    } else { 
     380        lll = pred.lognc(); 
     381    } 
     382 
     383    beta += yt; 
     384    return pred.lognc() - lll; 
    385385} 
    386386void multiBM::flatten ( const BMEF* B, double weight=1.0 ) { 
    387         const multiBM* E = dynamic_cast<const multiBM*> ( B ); 
    388         // sum(beta) should be equal to sum(B.beta) 
    389         const vec &Eb = E->beta;//const_cast<multiBM*> ( E )->_beta(); 
    390         beta *= ( sum ( Eb ) / sum ( beta ) * weight); 
    391         if ( evalll ) { 
    392                 last_lognc = est.lognc(); 
    393         } 
     387    const multiBM* E = dynamic_cast<const multiBM*> ( B ); 
     388    // sum(beta) should be equal to sum(B.beta) 
     389    const vec &Eb = E->beta;//const_cast<multiBM*> ( E )->_beta(); 
     390    beta *= ( sum ( Eb ) / sum ( beta ) * weight); 
     391    if ( evalll ) { 
     392        last_lognc = est.lognc(); 
     393    } 
    394394} 
    395395 
    396396vec egamma::sample() const { 
    397         vec smp ( dim ); 
    398         int i; 
    399  
    400         for ( i = 0; i < dim; i++ ) { 
    401                 if ( beta ( i ) > std::numeric_limits<double>::epsilon() ) { 
    402                         GamRNG.setup ( alpha ( i ), beta ( i ) ); 
    403                 } else { 
    404                         GamRNG.setup ( alpha ( i ), std::numeric_limits<double>::epsilon() ); 
    405                 } 
     397    vec smp ( dim ); 
     398    int i; 
     399 
     400    for ( i = 0; i < dim; i++ ) { 
     401        if ( beta ( i ) > std::numeric_limits<double>::epsilon() ) { 
     402            GamRNG.setup ( alpha ( i ), beta ( i ) ); 
     403        } else { 
     404            GamRNG.setup ( alpha ( i ), std::numeric_limits<double>::epsilon() ); 
     405        } 
    406406#pragma omp critical 
    407                 smp ( i ) = GamRNG(); 
    408         } 
    409  
    410         return smp; 
     407        smp ( i ) = GamRNG(); 
     408    } 
     409 
     410    return smp; 
    411411} 
    412412 
    413413 
    414414double egamma::evallog ( const vec &val ) const { 
    415         double res = 0.0; //the rest will be added 
    416         int i; 
    417  
    418         if ( any ( val <= 0. ) ) return -inf; 
    419         if ( any ( beta <= 0. ) ) return -inf; 
    420         for ( i = 0; i < dim; i++ ) { 
    421                 res += ( alpha ( i ) - 1 ) * std::log ( val ( i ) ) - beta ( i ) * val ( i ); 
    422         } 
    423         double tmp = res - lognc();; 
    424         bdm_assert_debug ( std::isfinite ( tmp ), "Infinite value" ); 
    425         return tmp; 
     415    double res = 0.0; //the rest will be added 
     416    int i; 
     417 
     418    if ( any ( val <= 0. ) ) return -inf; 
     419    if ( any ( beta <= 0. ) ) return -inf; 
     420    for ( i = 0; i < dim; i++ ) { 
     421        res += ( alpha ( i ) - 1 ) * std::log ( val ( i ) ) - beta ( i ) * val ( i ); 
     422    } 
     423    double tmp = res - lognc();; 
     424    bdm_assert_debug ( std::isfinite ( tmp ), "Infinite value" ); 
     425    return tmp; 
    426426} 
    427427 
    428428double egamma::lognc() const { 
    429         double res = 0.0; //will be added 
    430         int i; 
    431  
    432         for ( i = 0; i < dim; i++ ) { 
    433                 res += lgamma ( alpha ( i ) ) - alpha ( i ) * std::log ( beta ( i ) ) ; 
    434         } 
    435  
    436         return res; 
     429    double res = 0.0; //will be added 
     430    int i; 
     431 
     432    for ( i = 0; i < dim; i++ ) { 
     433        res += lgamma ( alpha ( i ) ) - alpha ( i ) * std::log ( beta ( i ) ) ; 
     434    } 
     435 
     436    return res; 
    437437} 
    438438 
    439439void egamma::from_setting ( const Setting &set ) { 
    440         epdf::from_setting ( set ); // reads rv 
    441         UI::get ( alpha, set, "alpha", UI::compulsory ); 
    442         UI::get ( beta, set, "beta", UI::compulsory ); 
    443 } 
    444          
     440    epdf::from_setting ( set ); // reads rv 
     441    UI::get ( alpha, set, "alpha", UI::compulsory ); 
     442    UI::get ( beta, set, "beta", UI::compulsory ); 
     443} 
     444 
    445445void egamma::to_setting ( Setting &set ) const 
    446         {                        
    447                 epdf::to_setting( set ); 
    448                 UI::save( alpha, set, "alpha" ); 
    449                 UI::save( beta, set, "beta" ); 
    450         }  
    451          
    452          
     446{ 
     447    epdf::to_setting( set ); 
     448    UI::save( alpha, set, "alpha" ); 
     449    UI::save( beta, set, "beta" ); 
     450} 
     451 
     452 
    453453void egamma::validate() { 
    454                 eEF::validate(); 
    455                 bdm_assert ( alpha.length() == beta.length(), "parameters do not match" ); 
    456                 dim = alpha.length(); 
    457         } 
     454    eEF::validate(); 
     455    bdm_assert ( alpha.length() == beta.length(), "parameters do not match" ); 
     456    dim = alpha.length(); 
     457} 
    458458 
    459459void mgamma::set_parameters ( double k0, const vec &beta0 ) { 
    460         k = k0; 
    461         iepdf.set_parameters ( k * ones ( beta0.length() ), beta0 ); 
     460    k = k0; 
     461    iepdf.set_parameters ( k * ones ( beta0.length() ), beta0 ); 
    462462} 
    463463 
    464464 
    465465void mgamma::from_setting ( const Setting &set ) { 
    466                 pdf::from_setting ( set ); // reads rv and rvc 
    467                 vec betatmp; // ugly but necessary 
    468                 UI::get ( betatmp, set, "beta", UI::compulsory ); 
    469                 UI::get ( k, set, "k", UI::compulsory ); 
    470                 set_parameters ( k, betatmp ); 
    471         } 
     466    pdf::from_setting ( set ); // reads rv and rvc 
     467    vec betatmp; // ugly but necessary 
     468    UI::get ( betatmp, set, "beta", UI::compulsory ); 
     469    UI::get ( k, set, "k", UI::compulsory ); 
     470    set_parameters ( k, betatmp ); 
     471} 
    472472 
    473473void mgamma::to_setting  (Setting  &set) const { 
    474                 pdf::to_setting(set); 
    475                 UI::save( _beta, set, "beta"); 
    476                 UI::save( k, set, "k"); 
    477  
    478         } 
     474    pdf::to_setting(set); 
     475    UI::save( _beta, set, "beta"); 
     476    UI::save( k, set, "k"); 
     477 
     478} 
    479479 
    480480void mgamma::validate() { 
    481                 pdf_internal<egamma>::validate(); 
    482  
    483                 dim = _beta.length(); 
    484                 dimc = _beta.length(); 
    485         } 
     481    pdf_internal<egamma>::validate(); 
     482 
     483    dim = _beta.length(); 
     484    dimc = _beta.length(); 
     485} 
    486486 
    487487void eEmp::resample ( RESAMPLING_METHOD method ) { 
    488         ivec ind = zeros_i ( n ); 
    489         bdm::resample(w,ind,method); 
    490         // copy the internals according to ind 
    491         for (int i = 0; i < n; i++ ) { 
    492                 if ( ind ( i ) != i ) { 
    493                         samples ( i ) = samples ( ind ( i ) ); 
    494                 } 
    495                 w ( i ) = 1.0 / n; 
    496         } 
     488    ivec ind = zeros_i ( n ); 
     489    bdm::resample(w,ind,method); 
     490    // copy the internals according to ind 
     491    for (int i = 0; i < n; i++ ) { 
     492        if ( ind ( i ) != i ) { 
     493            samples ( i ) = samples ( ind ( i ) ); 
     494        } 
     495        w ( i ) = 1.0 / n; 
     496    } 
    497497} 
    498498 
    499499void resample ( const vec &w, ivec &ind, RESAMPLING_METHOD method ) { 
    500         int n = w.length(); 
    501         ind = zeros_i ( n ); 
    502         ivec N_babies = zeros_i ( n ); 
    503         vec cumDist = cumsum ( w ); 
    504         vec u ( n ); 
    505         int i, j, parent; 
    506         double u0; 
    507          
    508         switch ( method ) { 
    509                 case MULTINOMIAL: 
    510                         u ( n - 1 ) = pow ( UniRNG.sample(), 1.0 / n ); 
    511                          
    512                         for ( i = n - 2; i >= 0; i-- ) { 
    513                                 u ( i ) = u ( i + 1 ) * pow ( UniRNG.sample(), 1.0 / ( i + 1 ) ); 
    514                         } 
    515                          
    516                         break; 
    517                          
    518                 case STRATIFIED: 
    519                          
    520                         for ( i = 0; i < n; i++ ) { 
    521                                 u ( i ) = ( i + UniRNG.sample() ) / n; 
    522                         } 
    523                          
    524                         break; 
    525                          
    526                 case SYSTEMATIC: 
    527                         u0 = UniRNG.sample(); 
    528                          
    529                         for ( i = 0; i < n; i++ ) { 
    530                                 u ( i ) = ( i + u0 ) / n; 
    531                         } 
    532                          
    533                         break; 
    534                          
    535                 default: 
    536                         bdm_error ( "PF::resample(): Unknown resampling method" ); 
    537         } 
    538          
    539         // U is now full 
    540         j = 0; 
    541          
    542         for ( i = 0; i < n; i++ ) { 
    543                 while ( u ( i ) > cumDist ( j ) ) j++; 
    544                  
    545                 N_babies ( j ) ++; 
    546         } 
    547         // We have assigned new babies for each Particle 
    548         // Now, we fill the resulting index such that: 
    549         // * particles with at least one baby should not move * 
    550         // This assures that reassignment can be done inplace; 
    551          
    552         // find the first parent; 
    553         parent = 0; 
    554         while ( N_babies ( parent ) == 0 ) parent++; 
    555          
    556         // Build index 
    557         for ( i = 0; i < n; i++ ) { 
    558                 if ( N_babies ( i ) > 0 ) { 
    559                         ind ( i ) = i; 
    560                         N_babies ( i ) --; //this index was now replicated; 
    561                 } else { 
    562                         // test if the parent has been fully replicated 
    563                         // if yes, find the next one 
    564                         while ( ( N_babies ( parent ) == 0 ) || ( N_babies ( parent ) == 1 && parent > i ) ) parent++; 
    565                          
    566                         // Replicate parent 
    567                         ind ( i ) = parent; 
    568                          
    569                         N_babies ( parent ) --; //this index was now replicated; 
    570                 } 
    571                  
    572         } 
     500    int n = w.length(); 
     501    ind = zeros_i ( n ); 
     502    ivec N_babies = zeros_i ( n ); 
     503    vec cumDist = cumsum ( w ); 
     504    vec u ( n ); 
     505    int i, j, parent; 
     506    double u0; 
     507 
     508    switch ( method ) { 
     509    case MULTINOMIAL: 
     510        u ( n - 1 ) = pow ( UniRNG.sample(), 1.0 / n ); 
     511 
     512        for ( i = n - 2; i >= 0; i-- ) { 
     513            u ( i ) = u ( i + 1 ) * pow ( UniRNG.sample(), 1.0 / ( i + 1 ) ); 
     514        } 
     515 
     516        break; 
     517 
     518    case STRATIFIED: 
     519 
     520        for ( i = 0; i < n; i++ ) { 
     521            u ( i ) = ( i + UniRNG.sample() ) / n; 
     522        } 
     523 
     524        break; 
     525 
     526    case SYSTEMATIC: 
     527        u0 = UniRNG.sample(); 
     528 
     529        for ( i = 0; i < n; i++ ) { 
     530            u ( i ) = ( i + u0 ) / n; 
     531        } 
     532 
     533        break; 
     534 
     535    default: 
     536        bdm_error ( "PF::resample(): Unknown resampling method" ); 
     537    } 
     538 
     539    // U is now full 
     540    j = 0; 
     541 
     542    for ( i = 0; i < n; i++ ) { 
     543        while ( u ( i ) > cumDist ( j ) ) j++; 
     544 
     545        N_babies ( j ) ++; 
     546    } 
     547    // We have assigned new babies for each Particle 
     548    // Now, we fill the resulting index such that: 
     549    // * particles with at least one baby should not move * 
     550    // This assures that reassignment can be done inplace; 
     551 
     552    // find the first parent; 
     553    parent = 0; 
     554    while ( N_babies ( parent ) == 0 ) parent++; 
     555 
     556    // Build index 
     557    for ( i = 0; i < n; i++ ) { 
     558        if ( N_babies ( i ) > 0 ) { 
     559            ind ( i ) = i; 
     560            N_babies ( i ) --; //this index was now replicated; 
     561        } else { 
     562            // test if the parent has been fully replicated 
     563            // if yes, find the next one 
     564            while ( ( N_babies ( parent ) == 0 ) || ( N_babies ( parent ) == 1 && parent > i ) ) parent++; 
     565 
     566            // Replicate parent 
     567            ind ( i ) = parent; 
     568 
     569            N_babies ( parent ) --; //this index was now replicated; 
     570        } 
     571 
     572    } 
    573573} 
    574574 
    575575void eEmp::set_statistics ( const vec &w0, const epdf &epdf0 ) { 
    576         dim = epdf0.dimension(); 
    577         w = w0; 
    578         w /= sum ( w0 );//renormalize 
    579         n = w.length(); 
    580         samples.set_size ( n ); 
    581  
    582         for ( int i = 0; i < n; i++ ) { 
    583                 samples ( i ) = epdf0.sample(); 
    584         } 
     576    dim = epdf0.dimension(); 
     577    w = w0; 
     578    w /= sum ( w0 );//renormalize 
     579    n = w.length(); 
     580    samples.set_size ( n ); 
     581 
     582    for ( int i = 0; i < n; i++ ) { 
     583        samples ( i ) = epdf0.sample(); 
     584    } 
    585585} 
    586586 
    587587void eEmp::set_samples ( const epdf* epdf0 ) { 
    588         w = 1; 
    589         w /= sum ( w );//renormalize 
    590  
    591         for ( int i = 0; i < n; i++ ) { 
    592                 samples ( i ) = epdf0->sample(); 
    593         } 
     588    w = 1; 
     589    w /= sum ( w );//renormalize 
     590 
     591    for ( int i = 0; i < n; i++ ) { 
     592        samples ( i ) = epdf0->sample(); 
     593    } 
    594594} 
    595595 
    596596void migamma_ref::from_setting ( const Setting &set ) { 
    597         migamma::from_setting(set); 
    598         vec ref; 
    599         double k,l; 
    600  
    601         UI::get ( ref, set, "ref" , UI::compulsory ); 
    602         UI::get( k, set, "k", UI::compulsory ); 
    603         UI::get( l, set, "l", UI::compulsory ); 
    604         set_parameters ( k, ref, l ); 
     597    migamma::from_setting(set); 
     598    vec ref; 
     599    double k,l; 
     600 
     601    UI::get ( ref, set, "ref" , UI::compulsory ); 
     602    UI::get( k, set, "k", UI::compulsory ); 
     603    UI::get( l, set, "l", UI::compulsory ); 
     604    set_parameters ( k, ref, l ); 
    605605} 
    606606 
    607607void  migamma_ref::to_setting  (Setting  &set) const { 
    608         migamma::to_setting(set); 
    609         UI::save ( pow ( refl, 1/(1.0 - l) ), set, "ref");       
    610         UI::save(l,set,"l"); 
    611         UI::save(k,set,"k"); 
     608    migamma::to_setting(set); 
     609    UI::save ( pow ( refl, 1/(1.0 - l) ), set, "ref"); 
     610    UI::save(l,set,"l"); 
     611    UI::save(k,set,"k"); 
    612612} 
    613613 
    614614void mlognorm::from_setting ( const Setting &set ) { 
    615         pdf_internal<elognorm>::from_setting(set); 
    616         vec mu0; 
    617         double k; 
    618         UI::get ( mu0, set, "mu0", UI::compulsory ); 
    619         UI::get( k, set, "k", UI::compulsory ); 
    620         set_parameters ( mu0.length(), k ); 
    621         condition ( mu0 ); 
     615    pdf_internal<elognorm>::from_setting(set); 
     616    vec mu0; 
     617    double k; 
     618    UI::get ( mu0, set, "mu0", UI::compulsory ); 
     619    UI::get( k, set, "k", UI::compulsory ); 
     620    set_parameters ( mu0.length(), k ); 
     621    condition ( mu0 ); 
    622622} 
    623623 
    624624void mlognorm::to_setting  (Setting  &set) const { 
    625         pdf_internal<elognorm>::to_setting(set); 
    626         UI::save ( exp(mu + sig2), set, "mu0"); 
    627  
    628         // inversion of sig2 = 0.5 * log ( k * k + 1 ); 
    629         double k = sqrt( exp( 2 * sig2 ) - 1  ); 
    630         UI::save(k,set,"k"); 
     625    pdf_internal<elognorm>::to_setting(set); 
     626    UI::save ( exp(mu + sig2), set, "mu0"); 
     627 
     628    // inversion of sig2 = 0.5 * log ( k * k + 1 ); 
     629    double k = sqrt( exp( 2 * sig2 ) - 1  ); 
     630    UI::save(k,set,"k"); 
    631631} 
    632632 
    633633 
    634634void mlstudent::condition ( const vec &cond ) { 
    635         if ( cond.length() > 0 ) { 
    636                 iepdf._mu() = A * cond + mu_const; 
    637         } else { 
    638                 iepdf._mu() =  mu_const; 
    639         } 
    640         double zeta; 
    641         //ugly hack! 
    642         if ( ( cond.length() + 1 ) == Lambda.rows() ) { 
    643                 zeta = Lambda.invqform ( concat ( cond, vec_1 ( 1.0 ) ) ); 
    644         } else { 
    645                 zeta = Lambda.invqform ( cond ); 
    646         } 
    647         _R = Re; 
    648         _R *= ( 1 + zeta );// / ( nu ); << nu is in Re!!!!!! 
     635    if ( cond.length() > 0 ) { 
     636        iepdf._mu() = A * cond + mu_const; 
     637    } else { 
     638        iepdf._mu() =  mu_const; 
     639    } 
     640    double zeta; 
     641    //ugly hack! 
     642    if ( ( cond.length() + 1 ) == Lambda.rows() ) { 
     643        zeta = Lambda.invqform ( concat ( cond, vec_1 ( 1.0 ) ) ); 
     644    } else { 
     645        zeta = Lambda.invqform ( cond ); 
     646    } 
     647    _R = Re; 
     648    _R *= ( 1 + zeta );// / ( nu ); << nu is in Re!!!!!! 
    649649} 
    650650 
    651651void eEmp::qbounds ( vec &lb, vec &ub, double perc ) const { 
    652         // lb in inf so than it will be pushed below; 
    653         lb.set_size ( dim ); 
    654         ub.set_size ( dim ); 
    655         lb = std::numeric_limits<double>::infinity(); 
    656         ub = -std::numeric_limits<double>::infinity(); 
    657         int j; 
    658         for ( int i = 0; i < n; i++ ) { 
    659                 for ( j = 0; j < dim; j++ ) { 
    660                         if ( samples ( i ) ( j ) < lb ( j ) ) { 
    661                                 lb ( j ) = samples ( i ) ( j ); 
    662                         } 
    663                         if ( samples ( i ) ( j ) > ub ( j ) ) { 
    664                                 ub ( j ) = samples ( i ) ( j ); 
    665                         } 
    666                 } 
    667         } 
     652    // lb in inf so than it will be pushed below; 
     653    lb.set_size ( dim ); 
     654    ub.set_size ( dim ); 
     655    lb = std::numeric_limits<double>::infinity(); 
     656    ub = -std::numeric_limits<double>::infinity(); 
     657    int j; 
     658    for ( int i = 0; i < n; i++ ) { 
     659        for ( j = 0; j < dim; j++ ) { 
     660            if ( samples ( i ) ( j ) < lb ( j ) ) { 
     661                lb ( j ) = samples ( i ) ( j ); 
     662            } 
     663            if ( samples ( i ) ( j ) > ub ( j ) ) { 
     664                ub ( j ) = samples ( i ) ( j ); 
     665            } 
     666        } 
     667    } 
    668668} 
    669669 
    670670void eEmp::to_setting ( Setting &set ) const { 
    671                 epdf::to_setting( set ); 
    672                 UI::save ( samples, set, "samples" ); 
    673                 UI::save ( w, set, "w" ); 
    674         } 
     671    epdf::to_setting( set ); 
     672    UI::save ( samples, set, "samples" ); 
     673    UI::save ( w, set, "w" ); 
     674} 
    675675 
    676676void eEmp::from_setting ( const Setting &set ) { 
    677         epdf::from_setting( set ); 
    678                  
    679         UI::get( samples, set, "samples", UI::compulsory ); 
    680         UI::get ( w, set, "w", UI::compulsory ); 
    681 } 
    682  
    683 void eEmp::validate (){ 
    684         epdf::validate(); 
    685         bdm_assert (samples.length()==w.length(),"samples and weigths are of different lengths"); 
    686         n = w.length(); 
    687         if (n>0) 
    688         pdf::dim = samples ( 0 ).length(); 
    689 } 
    690          
     677    epdf::from_setting( set ); 
     678 
     679    UI::get( samples, set, "samples", UI::compulsory ); 
     680    UI::get ( w, set, "w", UI::compulsory ); 
     681} 
     682 
     683void eEmp::validate () { 
     684    epdf::validate(); 
     685    bdm_assert (samples.length()==w.length(),"samples and weigths are of different lengths"); 
     686    n = w.length(); 
     687    if (n>0) 
     688        pdf::dim = samples ( 0 ).length(); 
     689} 
     690 
    691691void eDirich::from_setting ( const Setting &set ) { 
    692         epdf::from_setting ( set ); 
    693         UI::get ( beta, set, "beta", UI::compulsory ); 
     692    epdf::from_setting ( set ); 
     693    UI::get ( beta, set, "beta", UI::compulsory ); 
    694694} 
    695695 
    696696void eDirich::validate() { 
    697                 //check rv 
    698                 eEF::validate(); 
    699                 dim = beta.length(); 
    700         } 
     697    //check rv 
     698    eEF::validate(); 
     699    dim = beta.length(); 
     700} 
    701701 
    702702void eDirich::to_setting ( Setting &set ) const 
    703         {                        
    704                 eEF::to_setting( set ); 
    705                 UI::save( beta, set, "beta" ); 
    706         }  
     703{ 
     704    eEF::to_setting( set ); 
     705    UI::save( beta, set, "beta" ); 
     706} 
    707707 
    708708void euni::from_setting ( const Setting &set ) { 
    709                 epdf::from_setting ( set ); // reads rv and rvc 
    710  
    711                 UI::get ( high, set, "high", UI::compulsory ); 
    712                 UI::get ( low, set, "low", UI::compulsory ); 
    713                 set_parameters ( low, high ); 
    714                  
    715         } 
    716          
     709    epdf::from_setting ( set ); // reads rv and rvc 
     710 
     711    UI::get ( high, set, "high", UI::compulsory ); 
     712    UI::get ( low, set, "low", UI::compulsory ); 
     713    set_parameters ( low, high ); 
     714 
     715} 
     716 
    717717void    euni::to_setting  (Setting  &set) const { 
    718                 epdf::to_setting ( set ); 
    719                 UI::save ( high, set, "high" ); 
    720                 UI::save ( low, set, "low" ); 
    721         } 
    722          
     718    epdf::to_setting ( set ); 
     719    UI::save ( high, set, "high" ); 
     720    UI::save ( low, set, "low" ); 
     721} 
     722 
    723723void euni::validate() { 
    724                 epdf::validate(); 
    725                 bdm_assert ( high.length() == low.length(), "Incompatible high and low vectors" ); 
    726                 dim = high.length(); 
    727                 bdm_assert ( min ( distance ) > 0.0, "bad support" ); 
    728         } 
    729  
    730 void mgdirac::from_setting(const Setting& set){ 
    731                         pdf::from_setting(set); 
    732                         g=UI::build<fnc>(set,"g",UI::compulsory); 
    733                         validate(); 
    734                 } 
    735 void mgdirac::to_setting(Setting &set) const{ 
    736                         pdf::to_setting(set); 
    737                         UI::save(g.get(), set, "g"); 
    738                 } 
     724    epdf::validate(); 
     725    bdm_assert ( high.length() == low.length(), "Incompatible high and low vectors" ); 
     726    dim = high.length(); 
     727    bdm_assert ( min ( distance ) > 0.0, "bad support" ); 
     728} 
     729 
     730void mgdirac::from_setting(const Setting& set) { 
     731    pdf::from_setting(set); 
     732    g=UI::build<fnc>(set,"g",UI::compulsory); 
     733    validate(); 
     734} 
     735void mgdirac::to_setting(Setting &set) const { 
     736    pdf::to_setting(set); 
     737    UI::save(g.get(), set, "g"); 
     738} 
    739739void mgdirac::validate() { 
    740                         pdf::validate(); 
    741                         dim = g->dimension(); 
    742                         dimc = g->dimensionc(); 
    743                 } 
     740    pdf::validate(); 
     741    dim = g->dimension(); 
     742    dimc = g->dimensionc(); 
     743} 
    744744 
    745745void mDirich::from_setting ( const Setting &set ) { 
    746                 pdf::from_setting ( set ); // reads rv and rvc 
    747                 if ( _rv()._dsize() > 0 ) { 
    748                         rvc = _rv().copy_t ( -1 ); 
    749                 } 
    750                 vec beta0; 
    751                 if ( !UI::get ( beta0, set, "beta0", UI::optional ) ) { 
    752                         beta0 = ones ( _rv()._dsize() ); 
    753                 } 
    754                 if ( !UI::get ( betac, set, "betac", UI::optional ) ) { 
    755                         betac = 0.1 * ones ( _rv()._dsize() ); 
    756                 } 
    757                 _beta = beta0; 
    758  
    759                 UI::get ( k, set, "k", UI::compulsory ); 
    760         } 
     746    pdf::from_setting ( set ); // reads rv and rvc 
     747    if ( _rv()._dsize() > 0 ) { 
     748        rvc = _rv().copy_t ( -1 ); 
     749    } 
     750    vec beta0; 
     751    if ( !UI::get ( beta0, set, "beta0", UI::optional ) ) { 
     752        beta0 = ones ( _rv()._dsize() ); 
     753    } 
     754    if ( !UI::get ( betac, set, "betac", UI::optional ) ) { 
     755        betac = 0.1 * ones ( _rv()._dsize() ); 
     756    } 
     757    _beta = beta0; 
     758 
     759    UI::get ( k, set, "k", UI::compulsory ); 
     760} 
    761761 
    762762void mDirich::to_setting  (Setting  &set) const { 
    763                 pdf::to_setting(set); 
    764                 UI::save( _beta, set, "beta0"); 
    765                 UI::save( betac, set, "betac"); 
    766                 UI::save ( k, set, "k" ); 
    767         } 
     763    pdf::to_setting(set); 
     764    UI::save( _beta, set, "beta0"); 
     765    UI::save( betac, set, "betac"); 
     766    UI::save ( k, set, "k" ); 
     767} 
    768768 
    769769 
    770770void mDirich::validate() { 
    771                 pdf_internal<eDirich>::validate(); 
    772                 bdm_assert ( _beta.length() == betac.length(), "beta0 and betac are not compatible" ); 
    773                 if ( _rv()._dsize() > 0 ) { 
    774                         bdm_assert ( ( _rv()._dsize() == dimension() ) , "Size of rv does not match with beta" ); 
    775                 } 
    776                 dimc = _beta.length(); 
    777         } 
     771    pdf_internal<eDirich>::validate(); 
     772    bdm_assert ( _beta.length() == betac.length(), "beta0 and betac are not compatible" ); 
     773    if ( _rv()._dsize() > 0 ) { 
     774        bdm_assert ( ( _rv()._dsize() == dimension() ) , "Size of rv does not match with beta" ); 
     775    } 
     776    dimc = _beta.length(); 
     777} 
    778778 
    779779 
    780780void mBeta::from_setting ( const Setting &set ) { 
    781         pdf::from_setting ( set ); // reads rv and rvc 
    782         if ( _rv()._dsize() > 0 ) { 
    783                 rvc = _rv().copy_t ( -1 ); 
    784         } 
    785         if ( !UI::get ( iepdf.beta, set, "beta", UI::optional ) ) { 
    786                 iepdf.beta = ones ( _rv()._dsize() ); 
    787         } 
    788         if ( !UI::get ( iepdf.alpha, set, "alpha", UI::optional ) ) { 
    789                 iepdf.alpha = ones ( _rv()._dsize() ); 
    790         } 
    791         if ( !UI::get ( betac, set, "betac", UI::optional ) ) { 
    792                 betac = 0.1 * ones ( _rv()._dsize() ); 
    793         } 
    794          
    795         UI::get ( k, set, "k", UI::compulsory ); 
     781    pdf::from_setting ( set ); // reads rv and rvc 
     782    if ( _rv()._dsize() > 0 ) { 
     783        rvc = _rv().copy_t ( -1 ); 
     784    } 
     785    if ( !UI::get ( iepdf.beta, set, "beta", UI::optional ) ) { 
     786        iepdf.beta = ones ( _rv()._dsize() ); 
     787    } 
     788    if ( !UI::get ( iepdf.alpha, set, "alpha", UI::optional ) ) { 
     789        iepdf.alpha = ones ( _rv()._dsize() ); 
     790    } 
     791    if ( !UI::get ( betac, set, "betac", UI::optional ) ) { 
     792        betac = 0.1 * ones ( _rv()._dsize() ); 
     793    } 
     794 
     795    UI::get ( k, set, "k", UI::compulsory ); 
    796796} 
    797797 
    798798void mBeta::to_setting  (Setting  &set) const { 
    799         pdf::to_setting(set); 
    800         UI::save( iepdf.beta, set, "beta"); 
    801         UI::save( betac, set, "betac"); 
    802         UI::save ( k, set, "k" ); 
    803 } 
    804  
    805 } 
     799    pdf::to_setting(set); 
     800    UI::save( iepdf.beta, set, "beta"); 
     801    UI::save( betac, set, "betac"); 
     802    UI::save ( k, set, "k" ); 
     803} 
     804 
     805}