simulate.hpp

Go to the documentation of this file.
00001 /***************************************************************************/
00020 #include <assert.h>
00021 
00022 namespace aureservoir
00023 {
00024 
00026 
00027 
00028 template <typename T>
00029 SimBase<T>::SimBase(ESN<T> *esn)
00030 {
00031   esn_=esn;
00032   reallocate();
00033 }
00034 
00035 template <typename T>
00036 void SimBase<T>::reallocate()
00037 {
00038   last_out_.resize(esn_->outputs_, 1);
00039   t_.resize(esn_->neurons_);
00040 }
00041 
00042 template <typename T>
00043 void SimBase<T>::setBPCutoffConst(T f1, T f2) throw(AUExcept)
00044 {
00045   std::string str = "SimBase::setBPCutoffConst: ";
00046   str += "this is not implemented in standard ESNs, ";
00047   str += "use e.g. SIM_BP !";
00048 
00049   throw AUExcept( str );
00050 }
00051 
00052 template <typename T>
00053 void SimBase<T>::setBPCutoff(const typename ESN<T>::DEVector &f1,
00054                              const typename ESN<T>::DEVector &f2)
00055   throw(AUExcept)
00056 {
00057   std::string str = "SimBase::setBPCutoff: ";
00058   str += "this is not implemented in standard ESNs, ";
00059   str += "use e.g. SIM_BP !";
00060 
00061   throw AUExcept( str );
00062 }
00063 
00064 template <typename T>
00065 void SimBase<T>::setIIRCoeff(const typename DEMatrix<T>::Type &B,
00066                              const typename DEMatrix<T>::Type &A,
00067                              int series)
00068   throw(AUExcept)
00069 {
00070   std::string str = "SimBase::setIIRCoeff: ";
00071   str += "this is not implemented in standard ESNs, ";
00072   str += "use e.g. SIM_FILTER !";
00073 
00074   throw AUExcept( str );
00075 }
00076 
00077 template <typename T>
00078 void SimBase<T>::initDelayLine(int index,
00079                                const typename DEVector<T>::Type &initbuf)
00080   throw(AUExcept)
00081 {
00082   std::string str = "SimBase::initDelayLines: ";
00083   str += "this is not implemented in standard ESNs, ";
00084   str += "use ESN with delay&sum readout, e.g. SIM_FILTER_DS !";
00085 
00086   throw AUExcept( str );
00087 }
00088 
00089 template <typename T>
00090 typename DEMatrix<T>::Type SimBase<T>::getDelays()
00091   throw(AUExcept)
00092 {
00093   std::string str = "SimBase::getDelays: ";
00094   str += "this is not implemented in standard ESNs, ";
00095   str += "use ESN with delay&sum readout, e.g. SIM_FILTER_DS !";
00096 
00097   throw AUExcept( str );
00098 }
00099 
00100 template <typename T>
00101 typename DEVector<T>::Type &SimBase<T>::getDelayBuffer(int output, int nr)
00102     throw(AUExcept)
00103 {
00104   std::string str = "SimBase::getDelayBuffer: ";
00105   str += "this is not implemented in standard ESNs, ";
00106   str += "use ESN with delay&sum readout, e.g. SIM_FILTER_DS !";
00107 
00108   throw AUExcept( str );
00109 }
00110 
00112 
00113 
00114 
00115 template <typename T>
00116 void SimStd<T>::simulate(const typename ESN<T>::DEMatrix &in,
00117                          typename ESN<T>::DEMatrix &out)
00118 {
00119   assert( in.numRows() == esn_->inputs_ );
00120   assert( out.numRows() == esn_->outputs_ );
00121   assert( in.numCols() == out.numCols() );
00122   assert( last_out_.numRows() == esn_->outputs_ );
00123 
00124   int steps = in.numCols();
00125   typename ESN<T>::DEMatrix::View
00126     Wout1 = esn_->Wout_(_,_(1, esn_->neurons_)),
00127     Wout2 = esn_->Wout_(_,_(esn_->neurons_+1, esn_->neurons_+esn_->inputs_));
00128 
00132 
00137 
00138 
00139   // First run with output from last simulation
00140 
00141   t_ = esn_->x_; // temp object needed for BLAS
00142   esn_->x_ = esn_->Win_*in(_,1) + esn_->W_*t_ + esn_->Wback_*last_out_(_,1);
00143   // add noise
00144   Rand<T>::uniform(t_, -1.*esn_->noise_, esn_->noise_);
00145   esn_->x_ += t_;
00146   esn_->reservoirAct_( esn_->x_.data(), esn_->x_.length() );
00147 
00148   // output = Wout * [x; in]
00149   last_out_(_,1) = Wout1*esn_->x_ + Wout2*in(_,1);
00150 
00151   // output activation
00152   esn_->outputAct_( last_out_.data(),
00153                     last_out_.numRows()*last_out_.numCols() );
00154   out(_,1) = last_out_(_,1);
00155 
00156 
00157   // the rest
00158 
00159   for(int n=2; n<=steps; ++n)
00160   {
00161     t_ = esn_->x_; // temp object needed for BLAS
00162     esn_->x_ = esn_->Win_*in(_,n) + esn_->W_*t_ + esn_->Wback_*out(_,n-1);
00163     // add noise
00164     Rand<T>::uniform(t_, -1.*esn_->noise_, esn_->noise_);
00165     esn_->x_ += t_;
00166     esn_->reservoirAct_( esn_->x_.data(), esn_->x_.length() );
00167 
00168     // output = Wout * [x; in]
00169     last_out_(_,1) = Wout1*esn_->x_ + Wout2*in(_,n);
00170 
00171     // output activation
00172     esn_->outputAct_( last_out_.data(),
00173                       last_out_.numRows()*last_out_.numCols() );
00174     out(_,n) = last_out_(_,1);
00175   }
00176 }
00177 
00179 
00180 
00181 
00182 template <typename T>
00183 void SimLI<T>::simulate(const typename ESN<T>::DEMatrix &in,
00184                         typename ESN<T>::DEMatrix &out)
00185 {
00186   assert( in.numRows() == esn_->inputs_ );
00187   assert( out.numRows() == esn_->outputs_ );
00188   assert( in.numCols() == out.numCols() );
00189   assert( last_out_.numRows() == esn_->outputs_ );
00190 
00191   int steps = in.numCols();
00192   typename ESN<T>::DEMatrix::View
00193     Wout1 = esn_->Wout_(_,_(1, esn_->neurons_)),
00194     Wout2 = esn_->Wout_(_,_(esn_->neurons_+1, esn_->neurons_+esn_->inputs_));
00195 
00197 
00198 
00199   // First run with output from last simulation
00200 
00201   t_ = esn_->x_; // temp object needed for BLAS
00202   esn_->x_ = (1. - esn_->init_params_[LEAKING_RATE])*t_ +
00203              esn_->Win_*in(_,1) + esn_->W_*t_ + esn_->Wback_*last_out_(_,1);
00204   // add noise
00205   Rand<T>::uniform(t_, -1.*esn_->noise_, esn_->noise_);
00206   esn_->x_ += t_;
00207   esn_->reservoirAct_( esn_->x_.data(), esn_->x_.length() );
00208 
00209   // output = Wout * [x; in]
00210   last_out_(_,1) = Wout1*esn_->x_ + Wout2*in(_,1);
00211 
00212   // output activation
00213   esn_->outputAct_( last_out_.data(),
00214                     last_out_.numRows()*last_out_.numCols() );
00215   out(_,1) = last_out_(_,1);
00216 
00217 
00218   // the rest
00219 
00220   for(int n=2; n<=steps; ++n)
00221   {
00222     t_ = esn_->x_; // temp object needed for BLAS
00223     esn_->x_ = (1. - esn_->init_params_[LEAKING_RATE])*t_ +
00224                esn_->Win_*in(_,n) + esn_->W_*t_ + esn_->Wback_*out(_,n-1);
00225     // add noise
00226     Rand<T>::uniform(t_, -1.*esn_->noise_, esn_->noise_);
00227     esn_->x_ += t_;
00228     esn_->reservoirAct_( esn_->x_.data(), esn_->x_.length() );
00229 
00230     // output = Wout * [x; in]
00231     last_out_(_,1) = Wout1*esn_->x_ + Wout2*in(_,n);
00232 
00233     // output activation
00234     esn_->outputAct_( last_out_.data(),
00235                       last_out_.numRows()*last_out_.numCols() );
00236     out(_,n) = last_out_(_,1);
00237   }
00238 }
00239 
00241 
00242 
00243 
00244 template <typename T>
00245 void SimBP<T>::setBPCutoffConst(T f1, T f2) throw(AUExcept)
00246 {
00247   typename ESN<T>::DEVector f1vec(esn_->neurons_);
00248   typename ESN<T>::DEVector f2vec(esn_->neurons_);
00249 
00250   std::fill_n( f1vec.data(), f1vec.length(), f1 );
00251   std::fill_n( f2vec.data(), f2vec.length(), f2 );
00252   filter_.setBPCutoff(f1vec,f2vec);
00253 }
00254 
00255 template <typename T>
00256 void SimBP<T>::setBPCutoff(const typename ESN<T>::DEVector &f1,
00257                            const typename ESN<T>::DEVector &f2)
00258   throw(AUExcept)
00259 {
00260   if( f1.length() != esn_->neurons_ )
00261     throw AUExcept("SimBP: f1 must have same length as reservoir neurons!");
00262   if( f2.length() != esn_->neurons_ )
00263     throw AUExcept("SimBP: f2 must have same length as reservoir neurons!");
00264 
00265   filter_.setBPCutoff(f1,f2);
00266 }
00267 
00268 template <typename T>
00269 void SimBP<T>::simulate(const typename ESN<T>::DEMatrix &in,
00270                         typename ESN<T>::DEMatrix &out)
00271 {
00272   assert( in.numRows() == esn_->inputs_ );
00273   assert( out.numRows() == esn_->outputs_ );
00274   assert( in.numCols() == out.numCols() );
00275   assert( last_out_.numRows() == esn_->outputs_ );
00276 
00277   int steps = in.numCols();
00278   typename ESN<T>::DEMatrix::View
00279     Wout1 = esn_->Wout_(_,_(1, esn_->neurons_)),
00280     Wout2 = esn_->Wout_(_,_(esn_->neurons_+1, esn_->neurons_+esn_->inputs_));
00281 
00283 
00284   // First run with output from last simulation
00285 
00286   // calc neuron activation
00287   t_ = esn_->x_;
00288   esn_->x_ = esn_->Win_*in(_,1) + esn_->W_*t_ + esn_->Wback_*last_out_(_,1);
00289   // add noise
00290   Rand<T>::uniform(t_, -1.*esn_->noise_, esn_->noise_);
00291   esn_->x_ += t_;
00292   esn_->reservoirAct_( esn_->x_.data(), esn_->x_.length() );
00293 
00294   // Bandpass Filtering
00295   filter_.calc(esn_->x_);
00296 
00297   // output = Wout * [x; in]
00298   last_out_(_,1) = Wout1*esn_->x_ + Wout2*in(_,1);
00299 
00300   // output activation
00301   esn_->outputAct_( last_out_.data(),
00302                     last_out_.numRows()*last_out_.numCols() );
00303   out(_,1) = last_out_(_,1);
00304 
00305 
00306   // the rest
00307 
00308   for(int n=2; n<=steps; ++n)
00309   {
00310     t_ = esn_->x_; // temp object needed for BLAS
00311     esn_->x_ = esn_->Win_*in(_,n) + esn_->W_*t_ + esn_->Wback_*out(_,n-1);
00312     // add noise
00313     Rand<T>::uniform(t_, -1.*esn_->noise_, esn_->noise_);
00314     esn_->x_ += t_;
00315     esn_->reservoirAct_( esn_->x_.data(), esn_->x_.length() );
00316 
00317     // Bandpass Filtering
00318     filter_.calc(esn_->x_);
00319 
00320     // output = Wout * [x; in]
00321     last_out_(_,1) = Wout1*esn_->x_ + Wout2*in(_,n);
00322 
00323     // output activation
00324     esn_->outputAct_( last_out_.data(),
00325                       last_out_.numRows()*last_out_.numCols() );
00326     out(_,n) = last_out_(_,1);
00327   }
00328 }
00329 
00331 
00332 
00333 
00334 template <typename T>
00335 void SimFilter<T>::setIIRCoeff(const typename DEMatrix<T>::Type &B,
00336                            const typename DEMatrix<T>::Type &A,
00337                            int series)
00338   throw(AUExcept)
00339 {
00340   if( B.numRows() != esn_->neurons_ )
00341     throw AUExcept("SimFilter: B must have same rows as reservoir neurons!");
00342   if( A.numRows() != esn_->neurons_ )
00343     throw AUExcept("SimFilter: A must have same rows as reservoir neurons!");
00344 
00345   filter_.setIIRCoeff(B,A,series);
00346 }
00347 
00348 template <typename T>
00349 void SimFilter<T>::simulate(const typename ESN<T>::DEMatrix &in,
00350                             typename ESN<T>::DEMatrix &out)
00351 {
00352   assert( in.numRows() == esn_->inputs_ );
00353   assert( out.numRows() == esn_->outputs_ );
00354   assert( in.numCols() == out.numCols() );
00355   assert( last_out_.numRows() == esn_->outputs_ );
00356 
00357   int steps = in.numCols();
00358   typename ESN<T>::DEMatrix::View
00359     Wout1 = esn_->Wout_(_,_(1, esn_->neurons_)),
00360     Wout2 = esn_->Wout_(_,_(esn_->neurons_+1, esn_->neurons_+esn_->inputs_));
00361 
00363 
00364   // First run with output from last simulation
00365 
00366   // calc neuron activation
00367   t_ = esn_->x_;
00368   esn_->x_ = esn_->Win_*in(_,1) + esn_->W_*t_ + esn_->Wback_*last_out_(_,1);
00369   // add noise
00370   Rand<T>::uniform(t_, -1.*esn_->noise_, esn_->noise_);
00371   esn_->x_ += t_;
00372   esn_->reservoirAct_( esn_->x_.data(), esn_->x_.length() );
00373 
00374   // IIR Filtering
00375   filter_.calc(esn_->x_);
00376 
00377   // output = Wout * [x; in]
00378   last_out_(_,1) = Wout1*esn_->x_ + Wout2*in(_,1);
00379 
00380   // output activation
00381   esn_->outputAct_( last_out_.data(),
00382                     last_out_.numRows()*last_out_.numCols() );
00383   out(_,1) = last_out_(_,1);
00384 
00385 
00386   // the rest
00387 
00388   for(int n=2; n<=steps; ++n)
00389   {
00390     t_ = esn_->x_; // temp object needed for BLAS
00391     esn_->x_ = esn_->Win_*in(_,n) + esn_->W_*t_ + esn_->Wback_*out(_,n-1);
00392     // add noise
00393     Rand<T>::uniform(t_, -1.*esn_->noise_, esn_->noise_);
00394     esn_->x_ += t_;
00395     esn_->reservoirAct_( esn_->x_.data(), esn_->x_.length() );
00396 
00397     // IIR Filtering
00398     filter_.calc(esn_->x_);
00399 
00400     // output = Wout * [x; in]
00401     last_out_(_,1) = Wout1*esn_->x_ + Wout2*in(_,n);
00402 
00403     // output activation
00404     esn_->outputAct_( last_out_.data(),
00405                       last_out_.numRows()*last_out_.numCols() );
00406     out(_,n) = last_out_(_,1);
00407   }
00408 }
00409 
00411 
00412 
00413 
00414 template <typename T>
00415 void SimFilter2<T>::simulate(const typename ESN<T>::DEMatrix &in,
00416                             typename ESN<T>::DEMatrix &out)
00417 {
00418   assert( in.numRows() == esn_->inputs_ );
00419   assert( out.numRows() == esn_->outputs_ );
00420   assert( in.numCols() == out.numCols() );
00421   assert( last_out_.numRows() == esn_->outputs_ );
00422 
00423   int steps = in.numCols();
00424   typename ESN<T>::DEMatrix::View
00425     Wout1 = esn_->Wout_(_,_(1, esn_->neurons_)),
00426     Wout2 = esn_->Wout_(_,_(esn_->neurons_+1, esn_->neurons_+esn_->inputs_));
00427 
00429 
00430   // First run with output from last simulation
00431 
00432   t_ = esn_->x_; // temp object needed for BLAS
00433   esn_->x_ = esn_->Win_*in(_,1) + esn_->W_*t_ + esn_->Wback_*last_out_(_,1);
00434 
00435   // IIR Filtering
00436   filter_.calc(esn_->x_);
00437 
00438   // add noise
00439   Rand<T>::uniform(t_, -1.*esn_->noise_, esn_->noise_);
00440   esn_->x_ += t_;
00441   esn_->reservoirAct_( esn_->x_.data(), esn_->x_.length() );
00442 
00443   // output = Wout * [x; in]
00444   last_out_(_,1) = Wout1*esn_->x_ + Wout2*in(_,1);
00445 
00446   // output activation
00447   esn_->outputAct_( last_out_.data(),
00448                     last_out_.numRows()*last_out_.numCols() );
00449   out(_,1) = last_out_(_,1);
00450 
00451 
00452   // the rest
00453 
00454   for(int n=2; n<=steps; ++n)
00455   {
00456     t_ = esn_->x_; // temp object needed for BLAS
00457     esn_->x_ = esn_->Win_*in(_,n) + esn_->W_*t_ + esn_->Wback_*out(_,n-1);
00458 
00459     // IIR Filtering
00460     filter_.calc(esn_->x_);
00461 
00462     // add noise
00463     Rand<T>::uniform(t_, -1.*esn_->noise_, esn_->noise_);
00464     esn_->x_ += t_;
00465     esn_->reservoirAct_( esn_->x_.data(), esn_->x_.length() );
00466 
00467     // output = Wout * [x; in]
00468     last_out_(_,1) = Wout1*esn_->x_ + Wout2*in(_,n);
00469 
00470     // output activation
00471     esn_->outputAct_( last_out_.data(),
00472                       last_out_.numRows()*last_out_.numCols() );
00473     out(_,n) = last_out_(_,1);
00474   }
00475 }
00476 
00478 
00479 
00480 
00481 template <typename T>
00482 void SimFilterDS<T>::reallocate()
00483 {
00484   last_out_.resize(esn_->outputs_, 1);
00485   t_.resize(esn_->neurons_);
00486   dellines_.resize( (esn_->neurons_+esn_->inputs_)*esn_->outputs_ );
00487   intmp_.resize(esn_->inputs_,1);
00488 }
00489 
00490 template <typename T>
00491 void SimFilterDS<T>::initDelayLine(int index,
00492                                const typename DEVector<T>::Type &initbuf)
00493   throw(AUExcept)
00494 {
00495   assert( index >= 0 );
00496   assert( index < esn_->outputs_*(esn_->inputs_+esn_->neurons_) );
00497 
00498   dellines_[index].initBuffer(initbuf);
00499 }
00500 
00501 template <typename T>
00502 typename DEMatrix<T>::Type SimFilterDS<T>::getDelays() throw(AUExcept)
00503 {
00504   typename DEMatrix<T>::Type del(esn_->outputs_,
00505                                  esn_->inputs_+esn_->neurons_);
00506 
00507   for(int i=1; i<=esn_->outputs_;++i) {
00508   for(int j=1; j<=esn_->inputs_+esn_->neurons_; ++j) {
00509     del(i,j) = T( dellines_[(i-1)*(esn_->neurons_+esn_->inputs_)+j-1].delay_ );
00510   } }
00511 
00512   return del;
00513 }
00514 
00515 template <typename T>
00516 typename DEVector<T>::Type &SimFilterDS<T>::getDelayBuffer(int output, int nr)
00517     throw(AUExcept)
00518 {
00519   return dellines_[output*(esn_->neurons_+esn_->inputs_)+nr].buffer_;
00520 }
00521 
00522 template <typename T>
00523 void SimFilterDS<T>::simulate(const typename ESN<T>::DEMatrix &in,
00524                               typename ESN<T>::DEMatrix &out)
00525 {
00526   assert( in.numRows() == esn_->inputs_ );
00527   assert( out.numRows() == esn_->outputs_ );
00528   assert( in.numCols() == out.numCols() );
00529   assert( last_out_.numRows() == esn_->outputs_ );
00530 
00531   int steps = in.numCols();
00532   typename ESN<T>::DEMatrix::View
00533     Wout1 = esn_->Wout_(_,_(1, esn_->neurons_)),
00534     Wout2 = esn_->Wout_(_,_(esn_->neurons_+1, esn_->neurons_+esn_->inputs_));
00535 
00537 
00538   // First run with output from last simulation
00539 
00540   // calc neuron activation
00541   t_ = esn_->x_;
00542   esn_->x_ = esn_->Win_*in(_,1) + esn_->W_*t_ + esn_->Wback_*last_out_(_,1);
00543   // add noise
00544   Rand<T>::uniform(t_, -1.*esn_->noise_, esn_->noise_);
00545   esn_->x_ += t_;
00546   esn_->reservoirAct_( esn_->x_.data(), esn_->x_.length() );
00547 
00548   // IIR Filtering
00549   filter_.calc(esn_->x_);
00550 
00551   // delay states and inputs for all individual outputs
00552   for(int i=1; i<=esn_->outputs_; ++i)
00553   {
00554     // delay x_ vector and store into t_
00555     for(int j=1; j<=esn_->neurons_; ++j)
00556       t_(j) =  dellines_[(i-1)*(esn_->neurons_+esn_->inputs_)+j-1].tic(
00557                           esn_->x_(j) );
00558 
00559     // store correct delayed input vector in intmp_
00560     for(int j=1; j<=esn_->inputs_; ++j)
00561       intmp_(j,1) = dellines_[ (i-1)*(esn_->neurons_+esn_->inputs_)
00562                                 +esn_->neurons_+j-1 ].tic( in(j,1) );
00563 
00564     // calc  Wout * [x; in] for current output with delayed values
00565     last_out_(i,1) = Wout1(i,_)*t_ + Wout2(i,_)*intmp_(_,1);
00566   }
00567 
00568   // output activation
00569   esn_->outputAct_( last_out_.data(),
00570                     last_out_.numRows()*last_out_.numCols() );
00571   out(_,1) = last_out_(_,1);
00572 
00573 
00574   // the rest
00575 
00576   for(int n=2; n<=steps; ++n)
00577   {
00578     t_ = esn_->x_; // temp object needed for BLAS
00579     esn_->x_ = esn_->Win_*in(_,n) + esn_->W_*t_ + esn_->Wback_*out(_,n-1);
00580     // add noise
00581     Rand<T>::uniform(t_, -1.*esn_->noise_, esn_->noise_);
00582     esn_->x_ += t_;
00583     esn_->reservoirAct_( esn_->x_.data(), esn_->x_.length() );
00584 
00585     // IIR Filtering
00586     filter_.calc(esn_->x_);
00587 
00588     // delay states and inputs for all individual outputs
00589     for(int i=1; i<=esn_->outputs_; ++i)
00590     {
00591       // delay x_ vector and store into t_
00592       for(int j=1; j<=esn_->neurons_; ++j)
00593         t_(j) =  dellines_[(i-1)*(esn_->neurons_+esn_->inputs_)+j-1].tic(
00594                             esn_->x_(j) );
00595 
00596       // store correct delayed input vector in intmp_
00597       for(int j=1; j<=esn_->inputs_; ++j)
00598         intmp_(j,1) = dellines_[ (i-1)*(esn_->neurons_+esn_->inputs_)
00599                                   +esn_->neurons_+j-1 ].tic( in(j,n) );
00600 
00601       // calc  Wout * [x; in] for current output with delayed values
00602       last_out_(i,1) = Wout1(i,_)*t_ + Wout2(i,_)*intmp_(_,1);
00603     }
00604 
00605     // output activation
00606     esn_->outputAct_( last_out_.data(),
00607                       last_out_.numRows()*last_out_.numCols() );
00608     out(_,n) = last_out_(_,1);
00609   }
00610 }
00611 
00613 
00614 
00615 
00616 template <typename T>
00617 void SimSquare<T>::reallocate()
00618 {
00619   last_out_.resize(esn_->outputs_, 1);
00620   t_.resize(esn_->neurons_);
00621   t2_.resize(esn_->neurons_);
00622   dellines_.resize( (esn_->neurons_+esn_->inputs_)*esn_->outputs_ );
00623   intmp_.resize(esn_->inputs_,1);
00624   insq_.resize(esn_->inputs_);
00625 }
00626 
00627 template <typename T>
00628 void SimSquare<T>::simulate(const typename ESN<T>::DEMatrix &in,
00629                             typename ESN<T>::DEMatrix &out)
00630 {
00631   assert( in.numRows() == esn_->inputs_ );
00632   assert( out.numRows() == esn_->outputs_ );
00633   assert( in.numCols() == out.numCols() );
00634   assert( last_out_.numRows() == esn_->outputs_ );
00635 
00636   // we have to resize Wout_ to also have connections for
00637   // the squared states
00638   esn_->Wout_.resize(esn_->outputs_, 2*(esn_->neurons_+esn_->inputs_));
00639 
00640   int steps = in.numCols();
00641   typename ESN<T>::DEMatrix::View
00642     Wout1 = esn_->Wout_(_,_(1, esn_->neurons_)),
00643     Wout2 = esn_->Wout_(_,_(esn_->neurons_+1,esn_->neurons_+esn_->inputs_)),
00644     Wout3 = esn_->Wout_(_,_(esn_->neurons_+esn_->inputs_+1,
00645                             2*esn_->neurons_+esn_->inputs_)),
00646     Wout4 = esn_->Wout_(_,_(2*esn_->neurons_+esn_->inputs_+1,
00647                             2*(esn_->neurons_+esn_->inputs_)));
00648 
00649 
00650   // First run with output from last simulation
00651 
00652   t_ = esn_->x_; // temp object needed for BLAS
00653   esn_->x_ = esn_->Win_*in(_,1) + esn_->W_*t_ + esn_->Wback_*last_out_(_,1);
00654   // add noise
00655   Rand<T>::uniform(t_, -1.*esn_->noise_, esn_->noise_);
00656   esn_->x_ += t_;
00657   esn_->reservoirAct_( esn_->x_.data(), esn_->x_.length() );
00658 
00659   // IIR Filtering
00660   filter_.calc(esn_->x_);
00661 
00662   // delay states and inputs for all individual outputs
00663   for(int i=1; i<=esn_->outputs_; ++i)
00664   {
00665     // delay x_ vector and store into t_
00666     for(int j=1; j<=esn_->neurons_; ++j)
00667       t_(j) =  dellines_[(i-1)*(esn_->neurons_+esn_->inputs_)+j-1].tic(
00668                           esn_->x_(j) );
00669 
00670     // store correct delayed input vector in intmp_
00671     for(int j=1; j<=esn_->inputs_; ++j)
00672       intmp_(j,1) = dellines_[ (i-1)*(esn_->neurons_+esn_->inputs_)
00673                                 +esn_->neurons_+j-1 ].tic( in(j,1) );
00674 
00675     // calculate squared state version
00676     for(int j=1; j<=t_.length(); ++j)
00677       t2_(j) = pow( t_(j), 2 );
00678     // calculate squared input version
00679     for(int j=1; j<=insq_.length(); ++j)
00680       insq_(j) = pow( intmp_(j,1), 2 );
00681 
00682     // output = Wout * [x; in; x^2; in^2]
00683     last_out_(i,1) = Wout1(i,_)*t_ + Wout2(i,_)*intmp_(_,1)
00684                      + Wout3(i,_)*t2_ + Wout4(i,_)*insq_;
00685   }
00686 
00687   // output activation
00688   esn_->outputAct_( last_out_.data(),
00689                     last_out_.numRows()*last_out_.numCols() );
00690   out(_,1) = last_out_(_,1);
00691 
00692 
00693   // the rest
00694 
00695   for(int n=2; n<=steps; ++n)
00696   {
00697     t_ = esn_->x_; // temp object needed for BLAS
00698     esn_->x_ = esn_->Win_*in(_,n) + esn_->W_*t_ + esn_->Wback_*out(_,n-1);
00699     // add noise
00700     Rand<T>::uniform(t_, -1.*esn_->noise_, esn_->noise_);
00701     esn_->x_ += t_;
00702     esn_->reservoirAct_( esn_->x_.data(), esn_->x_.length() );
00703 
00704     // IIR Filtering
00705     filter_.calc(esn_->x_);
00706 
00707     // delay states and inputs for all individual outputs
00708     for(int i=1; i<=esn_->outputs_; ++i)
00709     {
00710       // delay x_ vector and store into t_
00711       for(int j=1; j<=esn_->neurons_; ++j)
00712         t_(j) =  dellines_[(i-1)*(esn_->neurons_+esn_->inputs_)+j-1].tic(
00713                             esn_->x_(j) );
00714 
00715       // store correct delayed input vector in intmp_
00716       for(int j=1; j<=esn_->inputs_; ++j)
00717         intmp_(j,1) = dellines_[ (i-1)*(esn_->neurons_+esn_->inputs_)
00718                                 +esn_->neurons_+j-1 ].tic( in(j,n) );
00719 
00720       // calculate squared state version
00721       for(int j=1; j<=t_.length(); ++j)
00722         t2_(j) = pow( t_(j), 2 );
00723       // calculate squared input version
00724       for(int j=1; j<=insq_.length(); ++j)
00725         insq_(j) = pow( intmp_(j,1), 2 );
00726 
00727       // output = Wout * [x; in; x^2; in^2]
00728       last_out_(i,1) = Wout1(i,_)*t_ + Wout2(i,_)*intmp_(_,1)
00729                        + Wout3(i,_)*t2_ + Wout4(i,_)*insq_;
00730     }
00731 
00732     // output activation
00733     esn_->outputAct_( last_out_.data(),
00734                       last_out_.numRows()*last_out_.numCols() );
00735     out(_,n) = last_out_(_,1);
00736   }
00737 }
00738 
00740 
00741 } // end of namespace aureservoir

Generated on Wed Mar 12 21:16:05 2008 for aureservoir by  doxygen 1.5.3