ADMB Documentation  11.5.3254
 All Classes Files Functions Variables Typedefs Friends Defines
fquadpri.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id$
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00007 #include <df1b2fun.h>
00008 #ifndef OPT_LIB
00009   #include <cassert>
00010   #include <climits>
00011 #endif
00012 
00013 // this should be a resizeable array
00014 df1b2quadratic_prior* df1b2quadratic_prior::ptr[100];
00015 
00016  int df1b2quadratic_prior::num_quadratic_prior=0;
00017  const int df1b2quadratic_prior::max_num_quadratic_prior=100;
00018 
00019 
00020 void df1b2quadratic_prior::add_to_list(void)
00021 {
00022   if (num_quadratic_prior >= max_num_quadratic_prior)
00023   {
00024     cerr << "Error[" << __FILE__ << ':' << __LINE__
00025          << "]: Max size exceeded.\n";
00026     ad_exit(1);
00027   }
00028   else
00029   {
00030     xmyindex=num_quadratic_prior;
00031     ptr[num_quadratic_prior++]=this;
00032   }
00033 }
00034 void df1b2quadratic_prior::get_Lxu(dmatrix& M)
00035 {
00036   bool isallocated = Lxu && allocated(*Lxu) && index;
00037 #ifndef OPT_LIB
00038   assert(isallocated);
00039 #endif
00040   if (isallocated)
00041   {
00042     int mmin=(*pu)(pu->indexmin()).get_ind_index();
00043     int size=pu->size();
00044     int offset=mmin-M(M.indexmin()).indexmax()-1;  // subtract x offset
00045     int nvar = index->indexmax();
00046     {
00047       switch(old_style_flag)
00048       {
00049       case 0:
00050         for (int i=1;i<=nvar;i++)
00051         {
00052           int jcol=(*index)(i);
00053           for (int ii=1;ii<=size;ii++)
00054           {
00055             M(ii+offset,jcol)+=(*Lxu)(i,ii);
00056           }
00057         }
00058         break;
00059       case 1:
00060         break;
00072       case 2:
00073         for (int i=1;i<=nvar;i++)
00074         {
00075           int jcol=(*index)(i);
00076           for (int ii=1;ii<=size;ii++)
00077           {
00078             M(ii+offset,jcol)+=(*Lxu)(i,ii);
00079           }
00080         }
00081         break;
00082       default:
00083          cerr << "Illegal valueinswitch statement" << endl;
00084          ad_exit(1);
00085       }
00086     }
00087   }
00088 }
00092 df1b2quadratic_prior::df1b2quadratic_prior():
00093   CM(NULL)
00094 {
00095   add_to_list();
00096   num_active_parameters=0;
00097   M = NULL;
00098   pu = NULL;
00099   index = NULL;
00100   Lxu = NULL;
00101 }
00105 df1b2quadratic_prior::~df1b2quadratic_prior(void)
00106 {
00107   if (index)
00108   {
00109     delete index;
00110     index = NULL;
00111   }
00112   if (Lxu)
00113   {
00114     delete Lxu;
00115     Lxu = NULL;
00116   }
00117   if (M)
00118   {
00119     delete M;
00120     M = NULL;
00121   }
00122   if (pu)
00123   {
00124     delete pu;
00125     pu = NULL;
00126   }
00127 }
00128 
00129 void df1b2quadratic_prior::allocate(const df1b2_init_vector& _u,
00130   const char* s)
00131   {
00132     allocate(_u);
00133   }
00134   void df1b2quadratic_prior::allocate(const df1b2matrix & _M,
00135     const df1b2_init_vector & _u, const char * s)
00136   {
00137     allocate(_M,_u);
00138   }
00139   void df1b2quadratic_prior::allocate(const df1b2_init_vector & _u)
00140   {
00141     pu = new df1b2_init_vector(_u);
00142   }
00143   void df1b2quadratic_prior::allocate(const df1b2matrix & _M,
00144     const df1b2_init_vector & _u)
00145   {
00146     M =new df1b2matrix(_M);
00147     pu = new df1b2_init_vector(_u);
00148   }
00149   void df1b2quadratic_prior::operator = (const df1b2matrix & M)
00150   {
00151     quadratic_prior::in_qp_calculations=0;
00152     num_active_parameters=funnel_init_var::num_vars;
00153     df1b2_gradlist::no_derivatives=1;
00154     dvector cu=value(*pu);
00155 
00156     if (laplace_approximation_calculator::where_are_we_flag==3)
00157     {
00158       df1b2variable::noallocate=1;
00159       df1b2vector v(M.indexmin(),M.indexmax());
00160       df1b2variable::noallocate=0;
00161       switch (old_style_flag)
00162       {
00163       case 0:
00164       case 1:
00165         v = solve(M,cu);
00166         break;
00167       case 2:
00168         v = M*cu;
00169         break;
00170       default:
00171         cerr << "Illegal value for quadratic_prior::old_style_flag"
00172              << endl;
00173         ad_exit(1);
00174       }
00175       int mmin=v.indexmin();
00176       int mmax=v.indexmax();
00177 
00178 #ifndef OPT_LIB
00179       assert(num_active_parameters <= INT_MAX);
00180 #endif
00181       if (index)
00182       {
00183         if (index->indexmax() != (int)num_active_parameters)
00184         {
00185           delete index;
00186           index=0;
00187         }
00188       }
00189 
00190       if (num_active_parameters>0)
00191       {
00192         if (!index)
00193         {
00194           index=new ivector(column(*funnel_init_var::plist,1));
00195         }
00196         if (Lxu)
00197         {
00198           int tmin = Lxu->indexmin();
00199           if ((Lxu->indexmin() != mmin)
00200               || (Lxu->indexmax() != mmax)
00201               || ((*Lxu)(tmin).indexmin() != 1)
00202               || ((*Lxu)(tmin).indexmax() != (int)num_active_parameters))
00203           {
00204             delete Lxu;
00205             Lxu=0;
00206           }
00207         }
00208         if (!Lxu)
00209         {
00210           Lxu=new dmatrix(1, (int)num_active_parameters, mmin-1, mmax);
00211         }
00212         for (int i=1;i<=(int)num_active_parameters;i++)
00213         {
00214           (*Lxu)(i,mmin-1)=(*funnel_init_var::plist)(i,1);
00215         }
00216         for (int j=mmin;j<=mmax;j++)
00217         {
00218           for (int i=1;i<=(int)num_active_parameters;i++)
00219           {
00220             switch (old_style_flag)
00221             {
00222             case 0:
00223               (*Lxu)(i,j)=v(j).get_u_dot()[i-1];
00224               break;
00225             case 1:
00226             case 2:
00227               (*Lxu)(i,j)=2.0*v(j).get_u_dot()[i-1];
00228               break;
00229             default:
00230               cerr << "Illegal value for quadratic_prior::old_style_flag"
00231                    << endl;
00232               ad_exit(1);
00233             }
00234           }
00235         }
00236       }
00237       else
00238       {
00239         if (Lxu)
00240         {
00241           delete Lxu;
00242           Lxu=0;
00243         }
00244       }
00245     }
00246     df1b2_gradlist::no_derivatives=0;
00247   }
00248 
00249  void df1b2quadratic_prior::get_Lxu_contribution(dmatrix& M)
00250  {
00251    for (int i=0;i<num_quadratic_prior;i++)
00252    {
00253      //cout << ptr[i]->get_num_active_parameters() << endl;
00254      //if (ptr[i]->get_num_active_parameters()>0)
00255      {
00256        ptr[i]->get_Lxu(M);
00257      }
00258    }
00259  }
00260 
00261 normal_df1b2quadratic_prior::normal_df1b2quadratic_prior(void)
00262 {
00263   set_old_style_flag();
00264 }
00265 
00266 void normal_df1b2quadratic_prior::set_old_style_flag(void)
00267 {
00268   old_style_flag=0;
00269 }
00270 void normal_df1b2quadratic_prior::operator = (const df1b2matrix & M)
00271 {
00272   df1b2quadratic_prior::operator = (M);
00273 }
00274 
00275 void df1b2quadratic_re_penalty::set_old_style_flag(void)
00276 {
00277   old_style_flag=2;
00278 }
00279 void df1b2quadratic_re_penalty::operator = (const df1b2matrix & M)
00280 {
00281   df1b2quadratic_prior::operator = (M);
00282 }
00283 void df1b2quadratic_re_penalty::operator = (const dmatrix & M)
00284 {
00285   df1b2quadratic_prior::operator = (M);
00286 }
00287 
00288 
00289 df1b2quadratic_re_penalty::df1b2quadratic_re_penalty(void)
00290 {
00291   set_old_style_flag();
00292 }
00293 // *******************************************************
00294 // *******************************************************
00295 // *******************************************************
00296 // *******************************************************
00297 constant_df1b2quadratic_re_penalty::constant_df1b2quadratic_re_penalty(void)
00298 {
00299   set_old_style_flag();
00300 }
00301 
00302 void constant_df1b2quadratic_re_penalty::set_old_style_flag(void)
00303 {
00304   old_style_flag=2;
00305 }
00306 void constant_df1b2quadratic_re_penalty::operator = (const dmatrix & M)
00307 {
00308   //df1b2quadratic_prior::operator = (M);
00309 }
00310 void df1b2quadratic_prior::operator = (const dmatrix & M)
00311 {
00312   quadratic_prior::in_qp_calculations=0;
00313   num_active_parameters=funnel_init_var::num_vars;
00314   df1b2_gradlist::no_derivatives=1;
00315   dvector cu=value(*pu);
00316 
00317   if (laplace_approximation_calculator::where_are_we_flag==3)
00318   {
00319     //df1b2variable::noallocate=1;
00320     //df1b2vector v(M.indexmin(),M.indexmax());
00321     df1b2variable::noallocate=0;
00322     switch (old_style_flag)
00323     {
00324     case 0:
00325     case 1:
00326       cout << "this can't happen" << endl;
00327       ad_exit(1);
00328       break;
00329     case 2:
00330       //v = M*cu;
00331       break;
00332     default:
00333       cerr << "Illegal value for quadratic_prior::old_style_flag"
00334            << endl;
00335       ad_exit(1);
00336     }
00337     //int mmin=v.indexmin();
00338     //int mmax=v.indexmax();
00339 
00340     if (index)
00341     {
00342       //if (index->indexmax() != nvar)
00343       cout << "this can't happen" << endl;
00344       ad_exit(1);
00345       delete index;
00346       index=0;
00347     }
00348 
00349     if (Lxu)
00350     {
00351       cout << "this can't happen" << endl;
00352       ad_exit(1);
00353       delete Lxu;
00354       Lxu=0;
00355     }
00356   }
00357   df1b2_gradlist::no_derivatives=0;
00358 }
00359 
00360