Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages

SOptimizer.cpp

Go to the documentation of this file.
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 //#define SOPT_TAU      0.0f
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     // construct the model space with UPPER and LOWER bound
00057     Vector<double> v(modelsize), upBound(ub),lowBound(lb);
00058     //v = 1.0f;
00059     v = minimize(M,y);  // use unconstrained solution as initial guess
00060     v.constrain(lb,ub);
00061     
00062     Vector<double> step((SVector&)gradstep);
00063     //step = 0.1;                               // for numerical gradient
00064     
00065     // constructing the initial model without UPPER and LOWER bound
00066     //Model<double> m(v);
00067     // constructing the initial model with UPPER and LOWER bound
00068     Model<double> m(upBound,lowBound,v);
00069 
00070     // Read the Right hand side
00071     //Vector<double>& rhs = h;
00072     
00073     //Constrcut optimization objects
00074     //LinearForward* lop = new LinearForward(&matA);
00075 
00076     //QuadraticObjFcn* qobj = new QuadraticObjFcn(modelsize, &matA, &h, 0.0f);
00077     SpectralObjFcn* obj = new SpectralObjFcn(&matA, h);
00078     
00079     /****IF USING ART, UNCOMMENT THE FOLLOWING LINE**********/
00080     // QuadraticOptima* opt = new ART(modelsize, lop, rhs, itMax, tol, verbose);
00081     /****IF USING SIRT, UNCOMMENT THE FOLLOWING LINE**********/
00082     // QuadraticOptima* opt = new SIRT(modelsize, lop, rhs, itMax, tol, verbose);
00083 
00084     /****IF USING IRLS, UNCOMMENT THE FOLLOWING LINE**********/
00085     //QuadraticOptima* opt = new IterativeReweightedLS( modelsize, lop, &rhs, 
00086     //                                                1, 2,itMax, tol, tol, verbose);
00087 
00088     /****IF USING CGLS, UNCOMMENT THE FOLLOWING LINE**********/
00089     //QuadraticOptima* opt = new LSConjugateGradient(modelsize, lop, &rhs,
00090     //                                             itMax, tol, verbose);
00091 
00092     /****IF USING CG, UNCOMMENT THE FOLLOWING LINE**********/
00093     CubicLineSearch* ls = new CubicLineSearch(obj, m_ItMaxLine, &step);
00094     Optima* opt = new ConjugateGradient(ls, m_ItMax, m_Tol, m_Verbose);
00095 
00096     //crashes in model??: Simplex* opt = new Simplex(obj, &m, itMax, 1, .5, 2, verbose);
00097 
00098     // DO the optimization, output is in mOpt
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         // Output Results and relative information
00109         ostream &outfp = cout;
00110         //outfp << "The best model" << endl << mOpt <<endl;
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; //<<opt->allResidue() << 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         // zeropad matrix to make enough rows 
00135         // (rank doesn't matter for the algorithm)
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 //#define FOR_SOME_REASON_IT_WORKS
00152 #ifdef FOR_SOME_REASON_IT_WORKS
00153         //this is beautiful, but it doesn't work!!!!!
00154         // Solution of Ax=b:
00155         linalg::Vector x = SVD_inv_mult(svd,b,SOPT_TAU);
00156         v = SVector(x);
00157 /* aparently, the matrix M is pretty ill conditioned sometimes, so that values
00158    on S, slightly above zero still matter. Matlab can handle that, but
00159    our algorithm fails here :-(
00160 */
00161 
00162 #else
00163         //so we do it by hand... (which improves results a little)
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         //create an inverted S and skip the nearly-zero elements
00172         // on the diagonal
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         //this sucks, too...
00183         //linalg::Matrix B = V * (S * (linalg::Matrix&)(transposed(U)));
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 }

Generated on Wed Dec 15 21:20:31 2004 for vuVolume by  doxygen 1.3.9.1