00001
00002 #ifndef mbl_ar_process_txx_
00003 #define mbl_ar_process_txx_
00004
00005
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
00017 template<class T>
00018 mbl_ar_process<T>::mbl_ar_process()
00019 {
00020 }
00021
00022
00023 template<class T>
00024 mbl_ar_process<T>::~mbl_ar_process()
00025 {
00026 }
00027
00028
00029 template<class T>
00030 short mbl_ar_process<T>::version_no() const
00031 {
00032 return 1;
00033 }
00034
00035
00036 template<class T>
00037 void mbl_ar_process<T>::print_summary(vcl_ostream& ) const
00038 {
00039 #if 0
00040 os << data_;
00041 #endif
00042 vcl_cerr << "mbl_ar_process<T>::print_summary() is NYI\n";
00043 }
00044
00045
00046 template <class T>
00047 void mbl_ar_process<T>::b_write(vsl_b_ostream& ) const
00048 {
00049 vcl_cout<<"mbl_ar_process<T>::b_write - NYI !\n";
00050 }
00051
00052
00053 template <class T>
00054 void mbl_ar_process<T>::b_read(vsl_b_istream& )
00055 {
00056 vcl_cout<<"mbl_ar_process<T>::b_read - NYI !\n";
00057 }
00058
00059
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
00067
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)
00072 {
00073 vnl_vector<T> Xm0;
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
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
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));
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
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);
00238 }
00239 else
00240 {
00241 vsl_b_write(os,true);
00242 p->b_write(os);
00243
00244 }
00245 }
00246
00247
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
00259 }
00260 else
00261 v = 0;
00262 }
00263
00264
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_