contrib/mul/mbl/mbl_ar_process.txx
Go to the documentation of this file.
00001 // This is mul/mbl/mbl_ar_process.txx
00002 #ifndef mbl_ar_process_txx_
00003 #define mbl_ar_process_txx_
00004 //:
00005 // \file
00006 
00007 #include "mbl_ar_process.h"
00008 
00009 #include <vnl/algo/vnl_symmetric_eigensystem.h>
00010 #include <vnl/algo/vnl_qr.h>
00011 #include <vnl/algo/vnl_matrix_inverse.h>
00012 #include <vnl/vnl_vector_fixed.h>
00013 #include <vnl/vnl_random.h>
00014 #include <vcl_cmath.h>
00015 
00016 //: Constructor
00017 template<class T>
00018 mbl_ar_process<T>::mbl_ar_process()
00019 {
00020 }
00021 
00022 //: Destructor
00023 template<class T>
00024 mbl_ar_process<T>::~mbl_ar_process()
00025 {
00026 }
00027 
00028 //: Version number for I/O
00029 template<class T>
00030 short mbl_ar_process<T>::version_no() const
00031 {
00032     return 1;
00033 }
00034 
00035 //: Print class to os
00036 template<class T>
00037 void mbl_ar_process<T>::print_summary(vcl_ostream& /*os*/) const
00038 {
00039 #if 0
00040     os << data_; // example of data output
00041 #endif
00042     vcl_cerr << "mbl_ar_process<T>::print_summary() is NYI\n";
00043 }
00044 
00045 //: Save class to binary file stream
00046 template  <class T>
00047 void mbl_ar_process<T>::b_write(vsl_b_ostream& /*bfs*/) const
00048 {
00049   vcl_cout<<"mbl_ar_process<T>::b_write - NYI !\n";
00050 }
00051 
00052 //: Load class from binary file stream
00053 template  <class T>
00054 void mbl_ar_process<T>::b_read(vsl_b_istream& /*bfs*/)
00055 {
00056   vcl_cout<<"mbl_ar_process<T>::b_read - NYI !\n";
00057 }
00058 
00059 //: Does the name of the class match the argument?
00060 template<class T>
00061 bool mbl_ar_process<T>::is_class(vcl_string const& s) const
00062 {
00063   return s==mbl_ar_process<T>::is_a();
00064 }
00065 
00066 //: Prediction
00067 // of a vcl_vector given the two previous vectors
00068 template<class T>
00069 vnl_vector<T> mbl_ar_process<T>::predict(vnl_vector<T>& Xm1,
00070                                          vnl_vector<T>& Xm2,
00071                                          vnl_random *rng/*=0*/)
00072 {
00073   vnl_vector<T> Xm0; // uninitialised
00074   if (Xm1.size()!=Xm2.size() || Xm1.size()!=Xm.size()) return Xm0;
00075 
00076   vnl_vector<T> wk(Xm.size());
00077 
00078   if (!rng)
00079   {
00080     static vnl_random mz_random;
00081     rng = &mz_random;
00082   }
00083 
00084   for (unsigned int i=0;i<Xm.size();i++)
00085   {
00086     wk[i]=(T)rng->normal();
00087   }
00088 
00089   return Xm+A_2*(Xm2-Xm)+A_1*(Xm1-Xm)+B_0*wk;
00090 }
00091 
00092 //: Learning using Burg's algorithm
00093 template<class T>
00094 void mbl_ar_process<T>::learn_burg(vcl_vector<vnl_vector<T> >& data)
00095 {
00096   if (data.size()<2) return;
00097 
00098   Xm=data[0];
00099   for (unsigned int i=1;i<data.size();i++)
00100     Xm+=data[i];
00101 
00102   Xm/=(T)data.size();
00103 
00104   for (unsigned int i=0;i<data.size();i++)
00105     data[i]-=Xm;
00106 
00107   unsigned int dim=data.size();
00108 
00109   vnl_vector<T> Ef(dim);
00110   vnl_vector<T> Eb(dim);
00111 
00112   A_1.set_size(data[0].size(),data[0].size());
00113   A_1.fill((T)0.0);
00114   A_2.set_size(data[0].size(),data[0].size());
00115   A_2.fill((T)0.0);
00116   B_0.set_size(data[0].size(),data[0].size());
00117   B_0.fill((T)0.0);
00118 
00119   for (unsigned int j=0;j<data[0].size();j++)
00120   {
00121     T E;
00122 
00123     for (unsigned int i=0;i<dim;i++)
00124     {
00125       Ef[i]=data[i][j];
00126       Eb[i]=data[i][j];
00127     }
00128 
00129     E=dot_product(Ef,Ef)/((T)dim);
00130 
00131     vnl_vector_fixed<T,3> a(T(1),T(0),T(0));
00132     T km;
00133     for (unsigned int i=0;i<2;i++)
00134     {
00135       vnl_vector<T> Efp(dim-i-1);
00136       vnl_vector<T> Ebp(dim-i-1);
00137       for (unsigned int k=0;k<dim-i-1;k++)
00138         {
00139         Efp[k]=Ef[k+1];
00140         Ebp[k]=Eb[k];
00141         }
00142 
00143       km=(-((T)(2.0))*dot_product(Ebp,Efp)/
00144         (dot_product(Ebp,Ebp)+dot_product(Efp,Efp)));
00145       for (unsigned int k=0;k<dim-i-1;k++)
00146         {
00147         Ef[k]=Efp[k]+km*Ebp[k];
00148         Eb[k]=Ebp[k]+km*Efp[k];
00149         }
00150       vnl_vector_fixed<T,3> b=a;
00151       for (unsigned int k=0;k<i+1;k++)
00152       {
00153         b[k+1]+=km*a[i-k];
00154       }
00155       a=b;
00156       E=(1-km*km)*E;
00157     }
00158     A_1[j][j]=-a[1];
00159     A_2[j][j]=-a[2];
00160     B_0[j][j]=E>((T)0.0)?vcl_sqrt(E):((T)0.0);
00161   }
00162 
00163   for (unsigned int i=0;i<data.size();i++)
00164     data[i]+=Xm;
00165 }
00166 
00167 //: Dynamic learning
00168 template<class T>
00169 void mbl_ar_process<T>::learn(vcl_vector<vnl_vector<T> >& data)
00170 {
00171   if (data.size()==0) return;
00172   unsigned int dim=data[0].size();
00173 
00174   vnl_vector<T> R0(dim,T(0)),R1(dim,T(0)),R2(dim,T(0)); // initialise to 0
00175   vnl_matrix<T> R00(dim,dim,T(0)), R01(dim,dim,T(0)), R02(dim,dim,T(0)),
00176                 R10(dim,dim,T(0)), R11(dim,dim,T(0)), R12(dim,dim,T(0)),
00177                 R20(dim,dim,T(0)), R21(dim,dim,T(0)), R22(dim,dim,T(0));
00178 
00179   Xm=data[0];
00180   for (unsigned int i=1;i<data.size();i++)
00181     Xm+=data[i];
00182   Xm/=(T)data.size();
00183 
00184   for (unsigned int i=0;i<data.size();i++)
00185     data[i]-=Xm;
00186 
00187   for (unsigned int i=2;i<data.size();i++)
00188   {
00189     R0+=data[i];
00190     R1+=data[i-1];
00191     R2+=data[i-2];
00192     R00+=outer_product(data[i],data[i]);
00193     R01+=outer_product(data[i],data[i-1]);
00194     R02+=outer_product(data[i],data[i-2]);
00195     R10+=outer_product(data[i-1],data[i]);
00196     R11+=outer_product(data[i-1],data[i-1]);
00197     R12+=outer_product(data[i-1],data[i-2]);
00198     R20+=outer_product(data[i-2],data[i]);
00199     R21+=outer_product(data[i-2],data[i-1]);
00200     R22+=outer_product(data[i-2],data[i-2]);
00201   }
00202 
00203   T coef=(T)1.0/((T)data.size()-(T)2.0);
00204 
00205   vnl_matrix<T> Rp01=R01-coef*outer_product(R0,R1);
00206   vnl_matrix<T> Rp02=R02-coef*outer_product(R0,R2);
00207   vnl_matrix<T> Rp11=R11-coef*outer_product(R1,R1);
00208   vnl_matrix<T> Rp12=R12-coef*outer_product(R1,R2);
00209   vnl_matrix<T> Rp21=R21-coef*outer_product(R2,R1);
00210   vnl_matrix<T> Rp22=R22-coef*outer_product(R2,R2);
00211 
00212   vnl_matrix_inverse<double> ti1A2(Rp11);
00213   vnl_matrix<T> t1A2=ti1A2*Rp12;
00214   vnl_matrix_inverse<double> ti2A2(Rp22-Rp21*t1A2);
00215 
00216   A_2=(Rp02-Rp01*t1A2)*ti2A2;
00217   A_1=(Rp01-A_2*Rp21)*ti1A2;
00218 
00219   vnl_vector<T> D=coef*(R0-A_2*R2-A_1*R1);
00220   vnl_matrix<T> C=coef*(R00-A_2*R20-A_1*R10-outer_product(D,R0));
00221   vnl_symmetric_eigensystem<T> srB(C);
00222   B_0=srB.square_root();
00223 
00224   vnl_qr<double> qr_A_1(A_1);
00225   vnl_qr<double> qr_A_2(A_2);
00226 
00227   for (unsigned int i=0;i<data.size();i++)
00228     data[i]+=Xm;
00229 }
00230 
00231 //: Write to binary stream
00232 template<class T>
00233 void vsl_b_write(vsl_b_ostream& os, const mbl_ar_process<T>* p)
00234 {
00235   if (p==0)
00236   {
00237     vsl_b_write(os, false); // Indicate null pointer stored
00238   }
00239   else
00240   {
00241     vsl_b_write(os,true); // Indicate non-null pointer stored
00242     p->b_write(os);
00243     //vsl_b_write(os,*p);
00244   }
00245 }
00246 
00247 //: Read data from binary stream
00248 template<class T>
00249 void vsl_b_read(vsl_b_istream& is, mbl_ar_process<T>* & v)
00250 {
00251   delete v;
00252   bool not_null_ptr;
00253   vsl_b_read(is, not_null_ptr);
00254   if (not_null_ptr)
00255   {
00256     v = new mbl_ar_process<T>();
00257     v->b_read(is);
00258     //vsl_b_read(is, *v);
00259   }
00260   else
00261     v = 0;
00262 }
00263 
00264 //: Print class to os
00265 template<class T>
00266 void vsl_print_summary(vcl_ostream& os, const mbl_ar_process<T>* p)
00267 {
00268   p->print_summary(os);
00269 }
00270 
00271 #undef MBL_AR_PROCESS_INSTANTIATE
00272 #define MBL_AR_PROCESS_INSTANTIATE(T) \
00273 VCL_DEFINE_SPECIALIZATION vcl_string mbl_ar_process<T >::is_a() const \
00274 { return vcl_string("mbl_ar_process<" #T ">"); } \
00275 template class mbl_ar_process<T >; \
00276 template void vsl_b_write(vsl_b_ostream& s, const mbl_ar_process<T >* arp); \
00277 template void vsl_b_read(vsl_b_istream& s, mbl_ar_process<T >* & arp); \
00278 template void vsl_print_summary(vcl_ostream& s,const mbl_ar_process<T >* arp)
00279 
00280 #endif //mbl_ar_process_txx_