00001 #include <iostream>
00002 #include <stdio.h>
00003
00004 #include "SOptimizer.h"
00005
00006 #include "Coool/C++.hh"
00007 #include "Coool/ObjFcns.hh"
00008 #include "Coool/Optim.hh"
00009 #include "SpecObj.hh"
00010
00011 #include "LinAlg/svd.h"
00012
00013
00014 #define SOPT_TAU 0.0000001f
00015
00016 SOptimizer::SOptimizer()
00017 {
00018 m_Verbose = 0;
00019 m_ItMax = 100;
00020 m_ItMaxLine = 15;
00021 m_Tol = 0.001;
00022 m_NIter = 0;
00023 m_Residue = 0;
00024 m_GradStep = 0.1;
00025 }
00026
00027 SOptimizer::~SOptimizer()
00028 {
00029 }
00030
00038 SVector SOptimizer::minimize(const SMatrix& M, const SVector& y,
00039 const SVector& lb, const SVector& ub)
00040 {
00041 SVector step(M.getNCols(), m_GradStep);
00042 return minimize(M, y, lb, ub, step);
00043 }
00044
00045 SVector SOptimizer::minimize(const SMatrix& M, const SVector& y,
00046 const SVector& lb, const SVector& ub,
00047 const SVector& gradstep)
00048 {
00049 using namespace coool;
00050
00051 int modelsize = M.getNCols();
00052
00053 DensMatrix<double> matA(M);
00054 Vector<double> h(y);
00055
00056
00057 Vector<double> v(modelsize), upBound(ub),lowBound(lb);
00058
00059 v = minimize(M,y);
00060 v.constrain(lb,ub);
00061
00062 Vector<double> step((SVector&)gradstep);
00063
00064
00065
00066
00067
00068 Model<double> m(upBound,lowBound,v);
00069
00070
00071
00072
00073
00074
00075
00076
00077 SpectralObjFcn* obj = new SpectralObjFcn(&matA, h);
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093 CubicLineSearch* ls = new CubicLineSearch(obj, m_ItMaxLine, &step);
00094 Optima* opt = new ConjugateGradient(ls, m_ItMax, m_Tol, m_Verbose);
00095
00096
00097
00098
00099 Model<double> mOpt(m);
00100 if(m_Verbose) cout << "starting optimization ..." << endl;
00101
00102 mOpt = opt->optimizer(m);
00103
00104 if(m_Verbose)
00105 {
00106 cout <<"The number of iterations: "<< opt->numIterations() << endl;
00107
00108
00109 ostream &outfp = cout;
00110
00111 outfp <<"The number of iterations: "<< opt->numIterations() << endl;
00112 outfp <<"The first residue: " << opt->firstResidue() << endl;
00113 outfp <<"The final residue: " << opt->finalResidue() << endl;
00114 outfp <<"The residue history: <available>" << endl;
00115 }
00116
00117 m_NIter = opt->numIterations();
00118 m_Residue = opt->finalResidue();
00119
00120 SVector ret(mOpt.modParam());
00121 return ret;
00122 }
00123
00128 SVector SOptimizer::minimize(const SMatrix& M, const SVector& y)
00129 {
00130 using namespace linalg;
00131 SVector v;
00132 try {
00133 if(!m_Verbose) freopen("/dev/null","a+",stderr);
00134
00135
00136 int d = M.getNCols() - M.getNRows();
00137 if(d<0) d=0;
00138
00139 linalg::Matrix mat(SMatrix(M.getNRows()+d,
00140 M.getNCols(),0.0f).insert(M));
00141 linalg::Vector b((SVector)SVector(M.getNRows()+d).insert(v));
00142
00143 SVD svd(mat);
00144
00145 if(m_Verbose) cout << "condition number of matrix A " <<
00146 svd.q_cond_number() << endl;
00147
00148 if(m_Verbose) cout << svd.q_U().q_nrows() << " != ? " <<
00149 b.q_nrows() << endl;
00150
00151
00152 #ifdef FOR_SOME_REASON_IT_WORKS
00153
00154
00155 linalg::Vector x = SVD_inv_mult(svd,b,SOPT_TAU);
00156 v = SVector(x);
00157
00158
00159
00160
00161
00162 #else
00163
00164 linalg::Vector sig(svd.q_sig());
00165 linalg::Matrix U(svd.q_U());
00166 linalg::Matrix V(svd.q_V());
00167
00168 linalg::Matrix S(sig.q_lwb(),sig.q_upb(),sig.q_lwb(),sig.q_upb());
00169 linalg::MatrixDA eth(S);
00170
00171
00172
00173 for(register int i=eth.q_row_lwb(); i<=eth.q_row_upb(); i++)
00174 for(register int j=eth.q_col_lwb(); j<=eth.q_col_upb(); j++)
00175 {
00176 if(i == j) {
00177 if(sig(i)>SOPT_TAU) eth(i,j) = 1/sig(i);
00178 else eth(i,j) = 0;
00179 } else eth(i,j) = 0;
00180 }
00181
00182
00183
00184 linalg::Matrix Ut = transposed(U);
00185 SMatrix mUt(Ut);
00186 SMatrix mV(V);
00187 SMatrix mS1(S);
00188 SMatrix mS = SMatrix(mV.getNCols(), mUt.getNRows(),0.0f).insert(mS1);
00189
00190 SMatrix B = mV * mS * mUt;
00191 v = B*SVector(B.getNCols(),0.0f).insert(y);
00192
00193 #define SHOW_SUV
00194 #ifdef SHOW_SUV
00195 if(m_Verbose) {
00196 cout << " M: " << M << endl;
00197 cout << mV.getNRows() << " " << mV.getNCols() << endl;
00198 cout << mS.getNRows() << " " << mS.getNCols() << endl;
00199 cout << mUt << endl;
00200 cout << mUt.getNRows() << " " << mUt.getNCols() << endl;
00201
00202 cout << "S^-1: " << mS << endl;
00203 cout << "Ut: " << mUt << endl;
00204 cout << "V: " << mV << endl;
00205 cout << "V(S^-1)U': " << B << endl;
00206
00207 cout << "colour: " << endl << M*v << endl;
00208 cout << "y : " << y << endl;
00209 cout << "v : " << v << " fine." << endl;
00210 }
00211 #endif //of SHOW_SUV
00212
00213 #endif //of FOR_SOME_REASON_IT_WORKS
00214
00215 m_NIter = 0;
00216 m_Residue = ((M*v) - y).norm();
00217
00218 } catch(void*)
00219 {
00220 if(!m_Verbose) freopen("/dev/stdout","a+",stderr);
00221 }
00222 if(!m_Verbose) freopen("/dev/stdout","a+",stderr);
00223
00224 return v;
00225 }