Changeset 622
- Timestamp:
- 09/16/09 22:52:57 (15 years ago)
- Files:
-
- 7 modified
-
applications/bdmtoolbox/mex/estimator.cpp (modified) (5 diffs)
-
applications/bdmtoolbox/mex/merger.cpp (modified) (2 diffs)
-
applications/bdmtoolbox/mex/mxstruct2config.cpp (modified) (1 diff)
-
applications/bdmtoolbox/mex/simulator.cpp (modified) (6 diffs)
-
applications/bdmtoolbox/tutorial/arx_mex_test.cfg (modified) (3 diffs)
-
applications/bdmtoolbox/tutorial/arx_test.cfg (modified) (3 diffs)
-
library/bdm/itpp_ext.cpp (modified) (13 diffs)
Legend:
- Unmodified
- Added
- Removed
-
applications/bdmtoolbox/mex/estimator.cpp
r611 r622 9 9 {rank="same"; "Data Source"; "Bayesian Model"} 10 10 "Data Source" -> "Bayesian Model" [label="data"]; 11 "Bayesian Model" -> "Result Logger" [label="estimated\n statistics"]; 11 "Bayesian Model" -> "Result Logger" [label="estimated\n statistics"]; 12 12 "Data Source" -> "Result Logger" [label="Simulated\n data"]; 13 13 } 14 14 \enddot 15 15 16 Here, 16 Here, 17 17 \li Data Source is an object (class DS) providing sequential data, \f$ [d_1, d_2, \ldots d_t] \f$. 18 18 \li Bayesian Model is an object (class BM) performing Bayesian filtering, … … 65 65 void mexFunction ( int n_output, mxArray *output[], int n_input, const mxArray *input[] ) { 66 66 // Check the number of inputs and output arguments 67 if (n_input<2) mexErrMsgTxt("Usage:\n"68 "result=estimator(system, estimators, experiment, logger)\n"69 " system = struct('class','datasource',...); % Estimated system\n"70 " estimators = {struct('class','estimator',...), % Estimators\n"71 " struct('class','estimator',...),...} \n"72 " === optional ==="73 " experiment = struct('ndat',10); % number of data in experiment, full length of finite datasources, 10 otherwise \n"74 " logger = struct('class','mexlogger'); % How to store results, default=mexlog, i.e. matlab structure\n\n"75 "see documentation of classes datasource, BM, and mexlogger and their offsprings in BDM.");76 67 if ( n_input<2 ) mexErrMsgTxt ( "Usage:\n" 68 "result=estimator(system, estimators, experiment, logger)\n" 69 " system = struct('class','datasource',...); % Estimated system\n" 70 " estimators = {struct('class','estimator',...), % Estimators\n" 71 " struct('class','estimator',...),...} \n" 72 " === optional ===" 73 " experiment = struct('ndat',10); % number of data in experiment, full length of finite datasources, 10 otherwise \n" 74 " logger = struct('class','mexlogger'); % How to store results, default=mexlog, i.e. matlab structure\n\n" 75 "see documentation of classes datasource, BM, and mexlogger and their offsprings in BDM." ); 76 77 77 RV::clear_all(); 78 //CONFIG 78 //CONFIG 79 79 UImxArray Cfg; 80 try {81 Cfg.addGroup(input[0],"system");82 Cfg.addList(input[1],"estimators");83 if (n_input>2){84 Cfg.addGroup(input[2],"experiment");85 }86 if (n_input>3){87 Cfg.addGroup(input[3],"logger");88 }/*else{80 try { 81 Cfg.addGroup ( input[0],"system" ); 82 Cfg.addList ( input[1],"estimators" ); 83 if ( n_input>2 ) { 84 Cfg.addGroup ( input[2],"experiment" ); 85 } 86 if ( n_input>3 ) { 87 Cfg.addGroup ( input[3],"logger" ); 88 }/*else{ 89 89 // define logger as mexlog 90 90 Setting &S=Cfg.getRoot(); … … 97 97 S["logger"]["maxlen"]=maxlen; 98 98 }*/ 99 } catch(SettingException e){it_error("error: "+string(e.getPath()));} 99 } catch ( SettingException e ) { 100 it_error ( "error: "+string ( e.getPath() ) ); 101 } 100 102 101 103 //DBG 102 Cfg.writeFile ("estimator.cfg");104 Cfg.writeFile ( "estimator.cfg" ); 103 105 104 106 #else … … 112 114 UIFile Cfg ( fname ); 113 115 #endif 114 115 shared_ptr<DS> Ds = UI::build<DS>( Cfg, "system" ); 116 Array<shared_ptr<BM> > Es; UI::get(Es,Cfg, "estimators" ); 116 117 shared_ptr<DS> Ds = UI::build<DS> ( Cfg, "system" ); 118 Array<shared_ptr<BM> > Es; 119 UI::get ( Es,Cfg, "estimators" ); 117 120 long Ndat=10; 118 if ( Cfg.exists("experiment")) {119 if ( Cfg.lookupValue ( "experiment.ndat",Ndat )) {120 bdm_assert (Ndat<=Ds->max_length(), "Data source has less data then required");121 if ( Cfg.exists ( "experiment" ) ) { 122 if ( Cfg.lookupValue ( "experiment.ndat",Ndat ) ) { 123 bdm_assert ( Ndat<=Ds->max_length(), "Data source has less data then required" ); 121 124 }; 122 125 } else { 123 if ( Ds->max_length() < inf) {126 if ( Ds->max_length() < inf ) { 124 127 Ndat=Ds->max_length(); 125 };// else Ndat=10; 128 } 129 ;// else Ndat=10; 126 130 } 127 shared_ptr<logger> L = UI::build<logger> ( Cfg, "logger",UI::optional);128 if ( !L) {129 131 shared_ptr<logger> L = UI::build<logger> ( Cfg, "logger",UI::optional ); 132 if ( !L ) { 133 #ifdef MEX 130 134 //mex logger has only from_setting constructor - we have to fill it... 131 L=new mexlog (Ndat);132 135 L=new mexlog ( Ndat ); 136 #else 133 137 L=new stdlog(); 134 138 #endif 135 139 } 136 140 137 141 Ds->log_add ( *L ); 138 142 string Ename; 139 143 Setting &S=Cfg; 140 for ( int i=0; i<Es.length(); i++){141 try {142 UI::get (Ename, S["estimators"][i], "name");143 } catch ( UIException e){144 Ename="Est"+num2str (i);144 for ( int i=0; i<Es.length(); i++ ) { 145 try { 146 UI::get ( Ename, S["estimators"][i], "name" ); 147 } catch ( UIException e ) { 148 Ename="Est"+num2str ( i ); 145 149 } 146 147 Es (i)->log_add(*L,Ename); // estimate150 151 Es ( i )->log_add ( *L,Ename ); // estimate 148 152 } 149 153 L->init(); 150 154 151 155 vec dt=zeros ( Ds->_drv()._dsize() ); //data variable 152 Array<datalink*> Dls (Es.length());153 for ( int i=0; i<Es.length(); i++){154 Dls (i)=new datalink( Es(i)->_drv(),Ds->_drv() ); //datalink between a datasource and estimator156 Array<datalink*> Dls ( Es.length() ); 157 for ( int i=0; i<Es.length(); i++ ) { 158 Dls ( i ) =new datalink ( Es ( i )->_drv(),Ds->_drv() ); //datalink between a datasource and estimator 155 159 } 156 160 157 161 for ( int tK=0;tK<Ndat;tK++ ) { 158 162 Ds->getdata ( dt ); // read data 159 163 Ds->logit ( *L ); 160 161 for ( int i=0; i<Es.length(); i++){162 Es (i)->bayes ( Dls(i)->pushdown ( dt ) ); // update estimates163 Es (i)->logit (*L);164 165 for ( int i=0; i<Es.length(); i++ ) { 166 Es ( i )->bayes ( Dls ( i )->pushdown ( dt ) ); // update estimates 167 Es ( i )->logit ( *L ); 164 168 } 165 169 L->step(); … … 171 175 172 176 #ifdef MEX 173 mexlog* mL=dynamic_cast<mexlog*> (L.get());177 mexlog* mL=dynamic_cast<mexlog*> ( L.get() ); 174 178 175 if ( mL) { // user wants output!!179 if ( mL ) { // user wants output!! 176 180 if ( n_output<1 ) mexErrMsgTxt ( "Wrong number of output variables!" ); 177 181 output[0] = mL->toCell(); -
applications/bdmtoolbox/mex/merger.cpp
r569 r622 29 29 30 30 #ifdef MEX 31 void mexFunction (int n_output, mxArray *output[], int n_input, const mxArray *input[]) 32 { 31 void mexFunction ( int n_output, mxArray *output[], int n_input, const mxArray *input[] ) { 33 32 // Check the number of inputs and output arguments 34 if ( n_input != 3) mexErrMsgTxt ("Usage:\n"35 "result=merger(sources, support, merger)\n"36 " sources= { struct('class','epdf'),... }; % cell of pdfs (epdfs or mpdfs) to be merged,\n"37 " support= struct(\n"38 " grid = {[dim1_start,dim1_end], [dim2_start, dim2_end]...} %support boundary \n"39 " nbins = [bins_in_dim1, bins_in_dim2,...] %fixed \n"40 " === OR ==\n"41 " pdf = struct('class','epdf'); % pdf to draw samples from\n"42 " nsamples= 100; % number of samples\n"43 " );\n"44 " If all elements are present, (grid,nbins) is used;\n"45 " merger = struct('class','merger_*'); % object to be used for merging,\n\n"46 "see documentation of classes epdf, mpdf, merger_base and their offsprings in BDM.");33 if ( n_input != 3 ) mexErrMsgTxt ( "Usage:\n" 34 "result=merger(sources, support, merger)\n" 35 " sources= { struct('class','epdf'),... }; % cell of pdfs (epdfs or mpdfs) to be merged,\n" 36 " support= struct(\n" 37 " grid = {[dim1_start,dim1_end], [dim2_start, dim2_end]...} %support boundary \n" 38 " nbins = [bins_in_dim1, bins_in_dim2,...] %fixed \n" 39 " === OR ==\n" 40 " pdf = struct('class','epdf'); % pdf to draw samples from\n" 41 " nsamples= 100; % number of samples\n" 42 " );\n" 43 " If all elements are present, (grid,nbins) is used;\n" 44 " merger = struct('class','merger_*'); % object to be used for merging,\n\n" 45 "see documentation of classes epdf, mpdf, merger_base and their offsprings in BDM." ); 47 46 RV::clear_all(); 48 47 // LOAD CONFIG 49 48 UImxArray Cfg; 50 Cfg.addList ( input[0], "Sources");51 Cfg.addGroup ( input[1], "Support");52 Cfg.addGroup ( input[2], "Merger");49 Cfg.addList ( input[0], "Sources" ); 50 Cfg.addGroup ( input[1], "Support" ); 51 Cfg.addGroup ( input[2], "Merger" ); 53 52 54 53 //DBG 55 Cfg.writeFile ( "merger.cfg");54 Cfg.writeFile ( "merger.cfg" ); 56 55 #else 57 56 int main() { 58 UIFile Cfg ( "merger.cfg");57 UIFile Cfg ( "merger.cfg" ); 59 58 #endif 60 59 // Sources 61 60 Array<shared_ptr<mpdf> > Sources; 62 61 //abuse Mer to store sources 63 Setting& _Sources = Cfg.lookup ( "Sources");62 Setting& _Sources = Cfg.lookup ( "Sources" ); 64 63 int Slen = _Sources.getLength(); 65 Sources.set_size ( Slen);66 for ( int i = 0; i < Slen; i++) {64 Sources.set_size ( Slen ); 65 for ( int i = 0; i < Slen; i++ ) { 67 66 try { 68 shared_ptr<mpdf> mtmp = UI::build<mpdf> ( _Sources, i);69 Sources ( i) = mtmp;70 } catch ( UIException) {67 shared_ptr<mpdf> mtmp = UI::build<mpdf> ( _Sources, i ); 68 Sources ( i ) = mtmp; 69 } catch ( UIException ) { 71 70 // it is not mpdf - see if it is epdf 72 71 try { 73 shared_ptr<epdf> etmp = UI::build<epdf> ( _Sources, i);74 if ( etmp) {75 Sources ( i) = new mepdf (etmp); // hopefully OK72 shared_ptr<epdf> etmp = UI::build<epdf> ( _Sources, i ); 73 if ( etmp ) { 74 Sources ( i ) = new mepdf ( etmp ); // hopefully OK 76 75 } 77 } catch ( UIException &e) {78 it_error ( "No mpdfs or epdfs found! " + string (e.what()));79 } catch ( std::exception e) {80 it_error ( "Error in UI at " + _Sources[i].getPath());76 } catch ( UIException &e ) { 77 it_error ( "No mpdfs or epdfs found! " + string ( e.what() ) ); 78 } catch ( std::exception e ) { 79 it_error ( "Error in UI at " + _Sources[i].getPath() ); 81 80 } 82 } catch ( std::exception &e) {83 it_error ( "Error in UI at " + _Sources[i].getPath());81 } catch ( std::exception &e ) { 82 it_error ( "Error in UI at " + _Sources[i].getPath() ); 84 83 } 85 84 } 86 85 87 shared_ptr<merger_base> Merger = UI::build<merger_base> ( Cfg, "Merger");86 shared_ptr<merger_base> Merger = UI::build<merger_base> ( Cfg, "Merger" ); 88 87 89 88 // Support 90 89 try { 91 shared_ptr<rectangular_support> RecSup = UI::build<rectangular_support> (Cfg, "Support"); 92 Merger->set_support(*RecSup); 93 } 94 catch (UIException &e) { 95 shared_ptr<discrete_support> DisSup = UI::build<discrete_support> (Cfg, "Support"); 96 Merger->set_support (*DisSup); 90 shared_ptr<rectangular_support> RecSup = UI::build<rectangular_support> ( Cfg, "Support" ); 91 Merger->set_support ( *RecSup ); 92 } catch ( UIException &e ) { 93 shared_ptr<discrete_support> DisSup = UI::build<discrete_support> ( Cfg, "Support" ); 94 Merger->set_support ( *DisSup ); 97 95 } 98 96 // COMPUTE RESULTS 99 Merger->set_sources ( Sources);97 Merger->set_sources ( Sources ); 100 98 Merger->merge(); 101 99 102 100 // save results 103 Array<vec> source_vals ( Sources.length());104 for ( int i = 0;i < Sources.length();i++) {101 Array<vec> source_vals ( Sources.length() ); 102 for ( int i = 0;i < Sources.length();i++ ) { 105 103 datalink_m2e dl; 106 dl.set_connection ( Sources (i)->_rv(), Sources (i)->_rvc(), Merger->_rv());104 dl.set_connection ( Sources ( i )->_rv(), Sources ( i )->_rvc(), Merger->_rv() ); 107 105 108 vec ll ( Merger->_Smp()._samples().length());109 for ( int j = 0; j < Merger->_Smp()._samples().length(); j++) {110 vec dt = dl.pushdown ( Merger->_Smp()._samples() (j));111 vec dtc = dl.get_cond ( Merger->_Smp()._samples() (j));112 ll ( j) = Sources (i)->evallogcond (dt, dtc);106 vec ll ( Merger->_Smp()._samples().length() ); 107 for ( int j = 0; j < Merger->_Smp()._samples().length(); j++ ) { 108 vec dt = dl.pushdown ( Merger->_Smp()._samples() ( j ) ); 109 vec dtc = dl.get_cond ( Merger->_Smp()._samples() ( j ) ); 110 ll ( j ) = Sources ( i )->evallogcond ( dt, dtc ); 113 111 } 114 112 115 vec sll = exp ( ll);113 vec sll = exp ( ll ); 116 114 117 source_vals ( i) = sll / sum (sll);115 source_vals ( i ) = sll / sum ( sll ); 118 116 } 119 117 120 merger_mix* MerMix=dynamic_cast<merger_mix*> (Merger.get());118 merger_mix* MerMix=dynamic_cast<merger_mix*> ( Merger.get() ); 121 119 vec mix_val; 122 123 if ( MerMix){124 vec ll ( Merger->_Smp()._samples().length());125 for ( int j = 0; j < Merger->_Smp()._samples().length(); j++) {126 ll ( j) = Merger->evallog (Merger->_Smp()._samples() (j));120 121 if ( MerMix ) { 122 vec ll ( Merger->_Smp()._samples().length() ); 123 for ( int j = 0; j < Merger->_Smp()._samples().length(); j++ ) { 124 ll ( j ) = Merger->evallog ( Merger->_Smp()._samples() ( j ) ); 127 125 } 128 129 vec sll = exp ( ll);130 131 mix_val = sll / sum ( sll);126 127 vec sll = exp ( ll ); 128 129 mix_val = sll / sum ( sll ); 132 130 } 133 131 … … 135 133 mxArray* tmp ; 136 134 // Save results 137 if ( n_output > 0) {138 tmp = mxCreateStructMatrix ( 1, 1, 0, NULL);135 if ( n_output > 0 ) { 136 tmp = mxCreateStructMatrix ( 1, 1, 0, NULL ); 139 137 //support 140 138 Array<vec> &samples = Merger->_Smp()._samples(); 141 if ( samples.size() > 0) {142 mxArray* fld = mxCreateDoubleMatrix ( samples (0).length(), samples.size(), mxREAL);143 Arrayvec2mxArray ( samples, fld);144 mxReplaceFieldNM ( tmp, "support", fld);139 if ( samples.size() > 0 ) { 140 mxArray* fld = mxCreateDoubleMatrix ( samples ( 0 ).length(), samples.size(), mxREAL ); 141 Arrayvec2mxArray ( samples, fld ); 142 mxReplaceFieldNM ( tmp, "support", fld ); 145 143 } 146 144 147 145 //weights 148 146 vec &w = Merger->_Smp()._w(); 149 mxArray* fldw = mxCreateDoubleMatrix ( 1, w.length(), mxREAL);150 vec2mxArray ( w, fldw);151 mxReplaceFieldNM ( tmp, "weights", fldw);147 mxArray* fldw = mxCreateDoubleMatrix ( 1, w.length(), mxREAL ); 148 vec2mxArray ( w, fldw ); 149 mxReplaceFieldNM ( tmp, "weights", fldw ); 152 150 153 151 //mixture values 154 if ( mix_val.length()>0){155 mxArray* fldm = mxCreateDoubleMatrix ( 1, w.length(), mxREAL);156 vec2mxArray ( mix_val, fldm);157 mxReplaceFieldNM ( tmp, "mix", fldm);152 if ( mix_val.length() >0 ) { 153 mxArray* fldm = mxCreateDoubleMatrix ( 1, w.length(), mxREAL ); 154 vec2mxArray ( mix_val, fldm ); 155 mxReplaceFieldNM ( tmp, "mix", fldm ); 158 156 } 159 157 // sources 160 158 char srcstr[20]; 161 for ( int i = 0;i < Sources.length();i++) {162 sprintf ( srcstr, "source%d", i + 1);163 mxArray* fldw = mxCreateDoubleMatrix ( 1, source_vals (i).length(), mxREAL);164 vec2mxArray ( source_vals (i), fldw);165 mxReplaceFieldNM ( tmp, srcstr, fldw);159 for ( int i = 0;i < Sources.length();i++ ) { 160 sprintf ( srcstr, "source%d", i + 1 ); 161 mxArray* fldw = mxCreateDoubleMatrix ( 1, source_vals ( i ).length(), mxREAL ); 162 vec2mxArray ( source_vals ( i ), fldw ); 163 mxReplaceFieldNM ( tmp, srcstr, fldw ); 166 164 } 167 165 -
applications/bdmtoolbox/mex/mxstruct2config.cpp
r391 r622 2 2 3 3 void mexFunction ( int n_output, mxArray *output[], int n_input, const mxArray *input[] ) { 4 UImxArray F ( input[0]);5 string filename = mxArray2string (input[1]);6 F.writeFile (filename.c_str());4 UImxArray F ( input[0] ); 5 string filename = mxArray2string ( input[1] ); 6 F.writeFile ( filename.c_str() ); 7 7 } -
applications/bdmtoolbox/mex/simulator.cpp
r617 r622 12 12 \enddot 13 13 14 Here, 14 Here, 15 15 \li Data Source is an object (class DS) providing sequential data, \f$ [d_1, d_2, \ldots d_t] \f$. 16 16 \li Result Logger is an object (class logger) dedicated to storing important data from the experiment. … … 56 56 void mexFunction ( int n_output, mxArray *output[], int n_input, const mxArray *input[] ) { 57 57 // Check the number of inputs and output arguments 58 if (n_input<1) mexErrMsgTxt("Usage:\n"59 "result=simulator(system, experiment, logger)\n"60 " system = struct('class','datasource',...); % Estimated system\n"61 " === optional ==="62 " experiment = struct('ndat',10); % number of data in experiment, full length of finite datasources, 10 otherwise \n"63 " logger = struct('class','mexlogger'); % How to store results, default=mexlog, i.e. matlab structure\n\n"64 "see documentation of classes datasource and mexlogger and their offsprings in BDM.");65 58 if ( n_input<1 ) mexErrMsgTxt ( "Usage:\n" 59 "result=simulator(system, experiment, logger)\n" 60 " system = struct('class','datasource',...); % Estimated system\n" 61 " === optional ===" 62 " experiment = struct('ndat',10); % number of data in experiment, full length of finite datasources, 10 otherwise \n" 63 " logger = struct('class','mexlogger'); % How to store results, default=mexlog, i.e. matlab structure\n\n" 64 "see documentation of classes datasource and mexlogger and their offsprings in BDM." ); 65 66 66 RV::clear_all(); 67 //CONFIG 67 //CONFIG 68 68 UImxArray Cfg; 69 try {70 Cfg.addGroup(input[0],"system");71 if (n_input>2){72 Cfg.addGroup(input[1],"experiment");73 }74 if (n_input>3){75 Cfg.addGroup(input[2],"logger");76 }/*else{69 try { 70 Cfg.addGroup ( input[0],"system" ); 71 if ( n_input>2 ) { 72 Cfg.addGroup ( input[1],"experiment" ); 73 } 74 if ( n_input>3 ) { 75 Cfg.addGroup ( input[2],"logger" ); 76 }/*else{ 77 77 // define logger as mexlog 78 78 Setting &S=Cfg.getRoot(); … … 85 85 S["logger"]["maxlen"]=maxlen; 86 86 }*/ 87 } catch(SettingException e){it_error("error: "+string(e.getPath()));} 87 } catch ( SettingException e ) { 88 it_error ( "error: "+string ( e.getPath() ) ); 89 } 88 90 89 91 //DBG 90 Cfg.writeFile ("simulator.cfg");92 Cfg.writeFile ( "simulator.cfg" ); 91 93 92 94 #else … … 100 102 UIFile Cfg ( fname ); 101 103 #endif 102 103 shared_ptr<DS> Ds = UI::build<DS> ( Cfg, "system" );104 105 shared_ptr<DS> Ds = UI::build<DS> ( Cfg, "system" ); 104 106 long Ndat=10; 105 if ( Cfg.exists("experiment")) {106 if ( Cfg.lookupValue ( "experiment.ndat",Ndat )) {107 bdm_assert (Ndat<=Ds->max_length(), "Data source has less data then required");107 if ( Cfg.exists ( "experiment" ) ) { 108 if ( Cfg.lookupValue ( "experiment.ndat",Ndat ) ) { 109 bdm_assert ( Ndat<=Ds->max_length(), "Data source has less data then required" ); 108 110 }; 109 111 } else { 110 if ( Ds->max_length() < std::numeric_limits< int >::max()) {112 if ( Ds->max_length() < std::numeric_limits< int >::max() ) { 111 113 Ndat=Ds->max_length(); 112 };// else Ndat=10; 114 } 115 ;// else Ndat=10; 113 116 } 114 shared_ptr<logger> L = UI::build<logger> ( Cfg, "logger",UI::optional);115 if ( !L) {116 117 shared_ptr<logger> L = UI::build<logger> ( Cfg, "logger",UI::optional ); 118 if ( !L ) { 119 #ifdef MEX 117 120 //mex logger has only from_setting constructor - we have to fill it... 118 L=new mexlog (Ndat);119 121 L=new mexlog ( Ndat ); 122 #else 120 123 L=new stdlog(); 121 124 #endif 122 125 } 123 126 124 127 Ds->log_add ( *L ); 125 128 L->init(); … … 130 133 Ds->getdata ( dt ); // read data 131 134 Ds->logit ( *L ); 132 135 133 136 L->step(); 134 137 Ds->step(); // simulator step … … 139 142 140 143 #ifdef MEX 141 mexlog* mL=dynamic_cast<mexlog*> (L.get());144 mexlog* mL=dynamic_cast<mexlog*> ( L.get() ); 142 145 143 if ( mL) { // user wants output!!146 if ( mL ) { // user wants output!! 144 147 if ( n_output<1 ) mexErrMsgTxt ( "Wrong number of output variables!" ); 145 148 output[0] = mL->toCell(); -
applications/bdmtoolbox/tutorial/arx_mex_test.cfg
r568 r622 1 1 //Data generating system 2 system = { type="external"; filename="arx_test.cfg"; path="system";}; 2 system = { type="external"; 3 filename="arx_test.cfg"; 4 path="system"; 5 }; 3 6 4 7 //store results … … 6 9 type= "mexlog"; 7 10 maxlen = 90; // 8 11 //dirname = "exp/ax"; 9 12 }; 10 13 11 14 //estimation 12 15 estimator = { 13 type = "ARX"; 14 y = {type="RV"; names=("y"); }; 16 type = "ARX"; 17 y = {type="RV"; 18 names= ( "y" ); 19 }; 15 20 rgr = {type="RV"; 16 names = ("y","y","y","u");17 18 };21 names = ( "y","y","y","u" ); 22 times = [-1, -2, -3, -1]; 23 }; 19 24 20 25 //optional fields 21 dV0 = matrix (5,1,[1e-3, 1e-5, 1e-5, 1e-5, 1e-5]); //default: 1e-3 for y, 1e-5 for rgr26 dV0 = matrix ( 5,1,[1e-3, 1e-5, 1e-5, 1e-5, 1e-5] ); //default: 1e-3 for y, 1e-5 for rgr 22 27 //nu0 = 8.; //default: rgrlen + 2 23 28 frg = .9991; // forgetting, default frg=1.0 … … 25 30 26 31 //experiment description 27 experiment: {32 experiment: { 28 33 ndat = 90; 29 34 }; -
applications/bdmtoolbox/tutorial/arx_test.cfg
r357 r622 2 2 system = { 3 3 class = "ArxDS"; 4 y = {class="RV"; names=("y", "u");}; 5 u = {class="RV"; names=(); }; 4 y = {class="RV"; 5 names= ( "y", "u" ); 6 }; 7 u = {class="RV"; 8 names=(); 9 }; 6 10 rgr = {class="RV"; 7 names = ("y","y","y","u");8 9 };11 names = ( "y","y","y","u" ); 12 times = [-1, -2, -3, -1]; 13 }; 10 14 //AR parameters 11 theta = ( 2,4,[0.8, -0.3, 0.4, 1.0,12 0.0, 0.0, 0.0, 0.0]);15 theta = ( 2,4,[0.8, -0.3, 0.4, 1.0, 16 0.0, 0.0, 0.0, 0.0] ); 13 17 // offset 14 18 offset = [0.0, 0.0]; 15 19 //variance 16 r = ( 2,2,[0.1, 0.0,17 0.0, 1.0] );20 r = ( 2,2,[0.1, 0.0, 21 0.0, 1.0] ); 18 22 // log also theta 19 23 opt="L_theta"; … … 30 34 estimator = { 31 35 class = "ARX"; 32 y = {class="RV"; names=("y"); }; 36 y = {class="RV"; 37 names= ( "y" ); 38 }; 33 39 rgr = {class="RV"; 34 names = ("y","y","y","u");35 36 };40 names = ( "y","y","y","u" ); 41 times = [-1, -2, -3, -1]; 42 }; 37 43 38 44 //optional fields … … 43 49 44 50 //experiment description 45 experiment: {51 experiment: { 46 52 ndat = 9000; 47 53 }; -
library/bdm/itpp_ext.cpp
r587 r622 18 18 // from algebra/lapack.h 19 19 20 extern "C" { /* QR factorization of a general matrix A */ 21 void dgeqrf_ ( int *m, int *n, double *a, int *lda, double *tau, double *work, 22 int *lwork, int *info ); 20 extern "C" /* QR factorization of a general matrix A */ 21 { 22 void dgeqrf_ (int *m, int *n, double *a, int *lda, double *tau, double *work, 23 int *lwork, int *info); 23 24 }; 24 25 25 namespace itpp { 26 Array<int> to_Arr ( const ivec &indices ) { 27 Array<int> a ( indices.size() ); 28 for ( int i = 0; i < a.size(); i++ ) { 29 a ( i ) = indices ( i ); 26 namespace itpp 27 { 28 Array<int> to_Arr (const ivec &indices) 29 { 30 Array<int> a (indices.size()); 31 for (int i = 0; i < a.size(); i++) { 32 a (i) = indices (i); 30 33 } 31 34 return a; 32 35 } 33 36 34 ivec linspace ( int from, int to ) { 37 ivec linspace (int from, int to) 38 { 35 39 int n = to - from + 1; 36 40 int i; 37 it_assert_debug ( n > 0, "wrong linspace");38 ivec iv ( n);39 for ( i = 0; i < n; i++ ) iv ( i) = from + i;41 it_assert_debug (n > 0, "wrong linspace"); 42 ivec iv (n); 43 for (i = 0; i < n; i++) iv (i) = from + i; 40 44 return iv; 41 45 }; 42 46 43 void set_subvector ( vec &ov, const ivec &iv, const vec &v ) { 44 it_assert_debug ( ( iv.length() <= v.length() ), 47 void set_subvector (vec &ov, const ivec &iv, const vec &v) 48 { 49 it_assert_debug ( (iv.length() <= v.length()), 45 50 "Vec<>::set_subvector(ivec, vec<Num_T>): Indexing out " 46 "of range of v" ); 47 for ( int i = 0; i < iv.length(); i++ ) { 48 it_assert_debug ( iv ( i ) < ov.length(), 49 "Vec<>::set_subvector(ivec, vec<Num_T>): Indexing out " 50 "of range of v" ); 51 ov ( iv ( i ) ) = v ( i ); 52 } 53 } 54 55 vec get_vec ( const vec &v, const ivec &indexlist ) { 51 "of range of v"); 52 for (int i = 0; i < iv.length(); i++) { 53 it_assert_debug (iv (i) < ov.length(), 54 "Vec<>::set_subvector(ivec, vec<Num_T>): Indexing out " 55 "of range of v"); 56 ov (iv (i)) = v (i); 57 } 58 } 59 60 vec get_vec (const vec &v, const ivec &indexlist) 61 { 56 62 int size = indexlist.size(); 57 vec temp ( size);58 for ( int i = 0; i < size; ++i) {59 temp ( i ) = v._data() [indexlist ( i) ];63 vec temp (size); 64 for (int i = 0; i < size; ++i) { 65 temp (i) = v._data() [indexlist (i) ]; 60 66 } 61 67 return temp; … … 68 74 #define R_FINITE std::isfinite 69 75 70 bvec operator> ( const vec &t1, const vec &t2 ) { 71 it_assert_debug ( t1.length() == t2.length(), "Vec<>::operator>(): different size of vectors" ); 72 bvec temp ( t1.length() ); 73 for ( int i = 0; i < t1.length(); i++ ) 74 temp ( i ) = ( t1[i] > t2[i] ); 76 bvec operator> (const vec &t1, const vec &t2) 77 { 78 it_assert_debug (t1.length() == t2.length(), "Vec<>::operator>(): different size of vectors"); 79 bvec temp (t1.length()); 80 for (int i = 0; i < t1.length(); i++) 81 temp (i) = (t1[i] > t2[i]); 75 82 return temp; 76 83 } 77 84 78 bvec operator< ( const vec &t1, const vec &t2 ) { 79 it_assert_debug ( t1.length() == t2.length(), "Vec<>::operator>(): different size of vectors" ); 80 bvec temp ( t1.length() ); 81 for ( int i = 0; i < t1.length(); i++ ) 82 temp ( i ) = ( t1[i] < t2[i] ); 85 bvec operator< (const vec &t1, const vec &t2) 86 { 87 it_assert_debug (t1.length() == t2.length(), "Vec<>::operator>(): different size of vectors"); 88 bvec temp (t1.length()); 89 for (int i = 0; i < t1.length(); i++) 90 temp (i) = (t1[i] < t2[i]); 83 91 return temp; 84 92 } 85 93 86 94 87 bvec operator& ( const bvec &a, const bvec &b ) { 88 it_assert_debug ( b.size() == a.size(), "operator&(): Vectors of different lengths" ); 89 90 bvec temp ( a.size() ); 91 for ( int i = 0; i < a.size(); i++ ) { 92 temp ( i ) = a ( i ) & b ( i ); 95 bvec operator& (const bvec &a, const bvec &b) 96 { 97 it_assert_debug (b.size() == a.size(), "operator&(): Vectors of different lengths"); 98 99 bvec temp (a.size()); 100 for (int i = 0; i < a.size(); i++) { 101 temp (i) = a (i) & b (i); 93 102 } 94 103 return temp; 95 104 } 96 105 97 bvec operator| ( const bvec &a, const bvec &b ) { 98 it_assert_debug ( b.size() != a.size(), "operator&(): Vectors of different lengths" ); 99 100 bvec temp ( a.size() ); 101 for ( int i = 0; i < a.size(); i++ ) { 102 temp ( i ) = a ( i ) | b ( i ); 106 bvec operator| (const bvec &a, const bvec &b) 107 { 108 it_assert_debug (b.size() != a.size(), "operator&(): Vectors of different lengths"); 109 110 bvec temp (a.size()); 111 for (int i = 0; i < a.size(); i++) { 112 temp (i) = a (i) | b (i); 103 113 } 104 114 return temp; … … 106 116 107 117 //#if 0 108 Gamma_RNG::Gamma_RNG ( double a, double b ) { 109 setup ( a, b ); 110 } 111 double Gamma_RNG::sample() { 118 Gamma_RNG::Gamma_RNG (double a, double b) 119 { 120 setup (a, b); 121 } 122 double Gamma_RNG::sample() 123 { 112 124 //A copy of rgamma code from the R package!! 113 125 // … … 147 159 double scale = 1.0 / beta; 148 160 149 if ( !R_FINITE ( a ) || !R_FINITE ( scale ) || a < 0.0 || scale <= 0.0) {150 it_error ( "Gamma_RNG wrong parameters");151 } 152 153 if ( a < 1. ) {/* GS algorithm for parameters a < 1 */154 if ( a == 0)161 if (!R_FINITE (a) || !R_FINITE (scale) || a < 0.0 || scale <= 0.0) { 162 it_error ("Gamma_RNG wrong parameters"); 163 } 164 165 if (a < 1.) { /* GS algorithm for parameters a < 1 */ 166 if (a == 0) 155 167 return 0.; 156 168 e = 1.0 + exp_m1 * a; 157 for ( ;; ) {//VS repeat169 for (;;) { //VS repeat 158 170 p = e * unif_rand(); 159 if ( p >= 1.0) {160 x = -log ( ( e - p ) / a);161 if ( exp_rand() >= ( 1.0 - a ) * log ( x ))171 if (p >= 1.0) { 172 x = -log ( (e - p) / a); 173 if (exp_rand() >= (1.0 - a) * log (x)) 162 174 break; 163 175 } else { 164 x = exp ( log ( p ) / a);165 if ( exp_rand() >= x)176 x = exp (log (p) / a); 177 if (exp_rand() >= x) 166 178 break; 167 179 } … … 173 185 174 186 /* Step 1: Recalculations of s2, s, d if a has changed */ 175 if ( a != aa) {187 if (a != aa) { 176 188 aa = a; 177 189 s2 = a - 0.5; 178 s = sqrt ( s2);190 s = sqrt (s2); 179 191 d = sqrt32 - s * 12.0; 180 192 } … … 186 198 x = s + 0.5 * t; 187 199 ret_val = x * x; 188 if ( t >= 0.0)200 if (t >= 0.0) 189 201 return scale * ret_val; 190 202 191 203 /* Step 3: u = 0,1 - uniform sample. squeeze acceptance (s) */ 192 204 u = unif_rand(); 193 if ( ( d * u ) <= ( t * t * t ))205 if ( (d * u) <= (t * t * t)) 194 206 return scale * ret_val; 195 207 196 208 /* Step 4: recalculations of q0, b, si, c if necessary */ 197 209 198 if ( a != aaa) {210 if (a != aaa) { 199 211 aaa = a; 200 212 r = 1.0 / a; 201 q0 = ( ( ( ( ( ( q7 * r + q6 ) * r + q5 ) * r + q4 ) * r + q3) * r202 + q2 ) * r + q1) * r;213 q0 = ( ( ( ( ( (q7 * r + q6) * r + q5) * r + q4) * r + q3) * r 214 + q2) * r + q1) * r; 203 215 204 216 /* Approximation depending on size of parameter a */ … … 206 218 /* were established by numerical experiments */ 207 219 208 if ( a <= 3.686) {220 if (a <= 3.686) { 209 221 b = 0.463 + s + 0.178 * s2; 210 222 si = 1.235; 211 223 c = 0.195 / s - 0.079 + 0.16 * s; 212 } else if ( a <= 13.022) {224 } else if (a <= 13.022) { 213 225 b = 1.654 + 0.0076 * s2; 214 226 si = 1.68 / s + 0.275; … … 222 234 /* Step 5: no quotient test if x not positive */ 223 235 224 if ( x > 0.0) {236 if (x > 0.0) { 225 237 /* Step 6: calculation of v and quotient q */ 226 v = t / ( s + s);227 if ( fabs ( v ) <= 0.25)228 q = q0 + 0.5 * t * t * ( ( ( ( ( ( a7 * v + a6 ) * v + a5 ) * v + a4) * v229 + a3 ) * v + a2 ) * v + a1) * v;238 v = t / (s + s); 239 if (fabs (v) <= 0.25) 240 q = q0 + 0.5 * t * t * ( ( ( ( ( (a7 * v + a6) * v + a5) * v + a4) * v 241 + a3) * v + a2) * v + a1) * v; 230 242 else 231 q = q0 - s * t + 0.25 * t * t + ( s2 + s2 ) * log ( 1.0 + v);243 q = q0 - s * t + 0.25 * t * t + (s2 + s2) * log (1.0 + v); 232 244 233 245 234 246 /* Step 7: quotient acceptance (q) */ 235 if ( log ( 1.0 - u ) <= q)247 if (log (1.0 - u) <= q) 236 248 return scale * ret_val; 237 249 } 238 250 239 for ( ;; ) {//VS repeat251 for (;;) { //VS repeat 240 252 /* Step 8: e = standard exponential deviate 241 253 * u = 0,1 -uniform deviate … … 244 256 u = unif_rand(); 245 257 u = u + u - 1.0; 246 if ( u < 0.0)258 if (u < 0.0) 247 259 t = b - si * e; 248 260 else 249 261 t = b + si * e; 250 262 /* Step 9: rejection if t < tau(1) = -0.71874483771719 */ 251 if ( t >= -0.71874483771719) {263 if (t >= -0.71874483771719) { 252 264 /* Step 10: calculation of v and quotient q */ 253 v = t / ( s + s);254 if ( fabs ( v ) <= 0.25)265 v = t / (s + s); 266 if (fabs (v) <= 0.25) 255 267 q = q0 + 0.5 * t * t * 256 ( ( ( ( ( ( a7 * v + a6 ) * v + a5 ) * v + a4 ) * v + a3) * v257 + a2 ) * v + a1) * v;268 ( ( ( ( ( (a7 * v + a6) * v + a5) * v + a4) * v + a3) * v 269 + a2) * v + a1) * v; 258 270 else 259 q = q0 - s * t + 0.25 * t * t + ( s2 + s2 ) * log ( 1.0 + v);271 q = q0 - s * t + 0.25 * t * t + (s2 + s2) * log (1.0 + v); 260 272 /* Step 11: hat acceptance (h) */ 261 273 /* (if q not positive go to step 8) */ 262 if ( q > 0.0) {274 if (q > 0.0) { 263 275 // TODO: w = expm1(q); 264 w = exp ( q) - 1;276 w = exp (q) - 1; 265 277 /* ^^^^^ original code had approximation with rel.err < 2e-7 */ 266 278 /* if t is rejected sample again at step 8 */ 267 if ( ( c * fabs ( u ) ) <= ( w * exp ( e - 0.5 * t * t ) ))279 if ( (c * fabs (u)) <= (w * exp (e - 0.5 * t * t))) 268 280 break; 269 281 } … … 275 287 276 288 277 bool qr ( const mat &A, mat &R ) { 289 bool qr (const mat &A, mat &R) 290 { 278 291 int info; 279 292 int m = A.rows(); 280 293 int n = A.cols(); 281 294 int lwork = n; 282 int k = std::min ( m, n);283 vec tau ( k);284 vec work ( lwork);295 int k = std::min (m, n); 296 vec tau (k); 297 vec work (lwork); 285 298 286 299 R = A; … … 288 301 // perform workspace query for optimum lwork value 289 302 int lwork_tmp = -1; 290 dgeqrf_ ( 291 &info);292 if ( info == 0) {293 lwork = static_cast<int> ( work ( 0 ));294 work.set_size ( lwork, false);295 } 296 dgeqrf_ ( &m, &n, R._data(), &m, tau._data(), work._data(), &lwork, &info);303 dgeqrf_ (&m, &n, R._data(), &m, tau._data(), work._data(), &lwork_tmp, 304 &info); 305 if (info == 0) { 306 lwork = static_cast<int> (work (0)); 307 work.set_size (lwork, false); 308 } 309 dgeqrf_ (&m, &n, R._data(), &m, tau._data(), work._data(), &lwork, &info); 297 310 298 311 // construct R 299 for ( int i = 0; i < m; i++)300 for ( int j = 0; j < std::min ( i, n ); j++)301 R ( i, j) = 0;302 303 return ( info == 0);312 for (int i = 0; i < m; i++) 313 for (int j = 0; j < std::min (i, n); j++) 314 R (i, j) = 0; 315 316 return (info == 0); 304 317 } 305 318 306 319 //#endif 307 std::string num2str ( double d ) { 320 std::string num2str (double d) 321 { 308 322 char tmp[20];//that should do 309 sprintf ( tmp, "%f", d);310 return std::string ( tmp);323 sprintf (tmp, "%f", d); 324 return std::string (tmp); 311 325 }; 312 std::string num2str ( int i ) { 326 std::string num2str (int i) 327 { 313 328 char tmp[10];//that should do 314 sprintf ( tmp, "%d", i);315 return std::string ( tmp);329 sprintf (tmp, "%d", i); 330 return std::string (tmp); 316 331 }; 317 332 … … 321 336 #define el 0.5772156649015329 322 337 323 double psi ( double x ) { 338 double psi (double x) 339 { 324 340 double s, ps, xa, x2; 325 341 int n, k; … … 335 351 }; 336 352 337 xa = fabs ( x);353 xa = fabs (x); 338 354 s = 0.0; 339 if ( ( x == ( int ) x ) && ( x <= 0.0 )) {355 if ( (x == (int) x) && (x <= 0.0)) { 340 356 ps = 1e308; 341 357 return ps; 342 358 } 343 if ( xa == ( int ) xa) {359 if (xa == (int) xa) { 344 360 n = xa; 345 for ( k = 1; k < n; k++) {361 for (k = 1; k < n; k++) { 346 362 s += 1.0 / k; 347 363 } 348 364 ps = s - el; 349 } else if ( ( xa + 0.5 ) == ( ( int ) ( xa + 0.5 ) )) {365 } else if ( (xa + 0.5) == ( (int) (xa + 0.5))) { 350 366 n = xa - 0.5; 351 for ( k = 1; k <= n; k++) {352 s += 1.0 / ( 2.0 * k - 1.0);367 for (k = 1; k <= n; k++) { 368 s += 1.0 / (2.0 * k - 1.0); 353 369 } 354 370 ps = 2.0 * s - el - 1.386294361119891; 355 371 } else { 356 if ( xa < 10.0) {357 n = 10 - ( int) xa;358 for ( k = 0; k < n; k++) {359 s += 1.0 / ( xa + k);372 if (xa < 10.0) { 373 n = 10 - (int) xa; 374 for (k = 0; k < n; k++) { 375 s += 1.0 / (xa + k); 360 376 } 361 377 xa += n; 362 378 } 363 x2 = 1.0 / ( xa * xa);364 ps = log ( xa ) - 0.5 / xa + x2 * ( ( ( ( ( ( ( a[7] * x2 + a[6] ) * x2 + a[5]) * x2 +365 a[4] ) * x2 + a[3] ) * x2 + a[2] ) * x2 + a[1] ) * x2 + a[0]);379 x2 = 1.0 / (xa * xa); 380 ps = log (xa) - 0.5 / xa + x2 * ( ( ( ( ( ( (a[7] * x2 + a[6]) * x2 + a[5]) * x2 + 381 a[4]) * x2 + a[3]) * x2 + a[2]) * x2 + a[1]) * x2 + a[0]); 366 382 ps -= s; 367 383 } 368 if ( x < 0.0)369 ps = ps - M_PI * std::cos ( M_PI * x ) / std::sin ( M_PI * x) - 1.0 / x;384 if (x < 0.0) 385 ps = ps - M_PI * std::cos (M_PI * x) / std::sin (M_PI * x) - 1.0 / x; 370 386 return ps; 371 387 } 372 388 373 void triu(mat &A){ 374 for(int i=1;i<A.rows();i++) { // row cycle 375 for (int j=0; j<i; j++) {A(i,j)=0;} 376 } 377 } 378 379 class RandunStorage{ 380 const int A; 381 const int M; 382 static double seed; 383 static int counter; 389 void triu (mat &A) 390 { 391 for (int i = 1;i < A.rows();i++) { // row cycle 392 for (int j = 0; j < i; j++) {A (i, j) = 0;} 393 } 394 } 395 396 class RandunStorage 397 { 398 const int A; 399 const int M; 400 static double seed; 401 static int counter; 384 402 public: 385 RandunStorage() : A(16807), M(2147483647) {};386 void set_seed (double seed0){seed=seed0;}387 double get() {403 RandunStorage() : A (16807), M (2147483647) {}; 404 void set_seed (double seed0) {seed = seed0;} 405 double get() { 388 406 long long tmp = A * seed; 389 407 tmp = tmp % M; 390 408 seed = tmp; 391 counter++; 392 return seed/M;} 409 counter++; 410 return seed / M; 411 } 393 412 }; 394 413 static RandunStorage randun_global_storage; 395 double RandunStorage::seed=1111111; 396 int RandunStorage::counter=0; 397 double randun(){return randun_global_storage.get();}; 398 vec randun(int n){vec res(n); for(int i=0;i<n;i++){res(i)=randun();}; return res;}; 399 mat randun(int n, int m){mat res(n,m); for(int i=0;i<n*m;i++){res(i)=randun();}; return res;}; 414 double RandunStorage::seed = 1111111; 415 int RandunStorage::counter = 0; 416 double randun() {return randun_global_storage.get();}; 417 vec randun (int n) {vec res (n); for (int i = 0;i < n;i++) {res (i) = randun();}; return res;}; 418 mat randun (int n, int m) {mat res (n, m); for (int i = 0;i < n*m;i++) {res (i) = randun();}; return res;}; 419 400 420 ivec unique (const ivec &in) 401 421 {
