ADMB Documentation  11.2.2828
 All Classes Files Functions Variables Typedefs Friends Defines
fquadpri.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: fquadpri.cpp 2760 2014-12-05 20:21:18Z johnoel $
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 }
00089 
00090   df1b2quadratic_prior::df1b2quadratic_prior(void)
00091   {
00092     add_to_list();
00093     num_active_parameters=0;
00094     M=0;
00095     pu=0;
00096     index=0;
00097     Lxu=0;
00098   }
00099   df1b2quadratic_prior::~df1b2quadratic_prior(void)
00100   {
00101     if (index) delete index;
00102     index=0;
00103     if (Lxu) delete Lxu;
00104     Lxu=0;
00105     if (M) delete M;
00106     M=0;
00107     if (pu) delete pu;
00108     pu=0;
00109   }
00110 
00111 void df1b2quadratic_prior::allocate(const df1b2_init_vector& _u,
00112   const char* s)
00113   {
00114     allocate(_u);
00115   }
00116   void df1b2quadratic_prior::allocate(const df1b2matrix & _M,
00117     const df1b2_init_vector & _u, const char * s)
00118   {
00119     allocate(_M,_u);
00120   }
00121   void df1b2quadratic_prior::allocate(const df1b2_init_vector & _u)
00122   {
00123     pu = new df1b2_init_vector(_u);
00124   }
00125   void df1b2quadratic_prior::allocate(const df1b2matrix & _M,
00126     const df1b2_init_vector & _u)
00127   {
00128     M =new df1b2matrix(_M);
00129     pu = new df1b2_init_vector(_u);
00130   }
00131   void df1b2quadratic_prior::operator = (const df1b2matrix & M)
00132   {
00133     quadratic_prior::in_qp_calculations=0;
00134     num_active_parameters=funnel_init_var::num_vars;
00135     df1b2_gradlist::no_derivatives=1;
00136     dvector cu=value(*pu);
00137 
00138     if (laplace_approximation_calculator::where_are_we_flag==3)
00139     {
00140       df1b2variable::noallocate=1;
00141       df1b2vector v(M.indexmin(),M.indexmax());
00142       df1b2variable::noallocate=0;
00143       switch (old_style_flag)
00144       {
00145       case 0:
00146       case 1:
00147         v = solve(M,cu);
00148         break;
00149       case 2:
00150         v = M*cu;
00151         break;
00152       default:
00153         cerr << "Illegal value for quadratic_prior::old_style_flag"
00154              << endl;
00155         ad_exit(1);
00156       }
00157       int mmin=v.indexmin();
00158       int mmax=v.indexmax();
00159 
00160 #ifndef OPT_LIB
00161       assert(num_active_parameters <= INT_MAX);
00162 #endif
00163       if (index)
00164       {
00165         if (index->indexmax() != (int)num_active_parameters)
00166         {
00167           delete index;
00168           index=0;
00169         }
00170       }
00171 
00172       if (num_active_parameters>0)
00173       {
00174         if (!index)
00175         {
00176           index=new ivector(column(*funnel_init_var::plist,1));
00177         }
00178         if (Lxu)
00179         {
00180           int tmin = Lxu->indexmin();
00181           if ((Lxu->indexmin() != mmin)
00182               || (Lxu->indexmax() != mmax)
00183               || ((*Lxu)(tmin).indexmin() != 1)
00184               || ((*Lxu)(tmin).indexmax() != (int)num_active_parameters))
00185           {
00186             delete Lxu;
00187             Lxu=0;
00188           }
00189         }
00190         if (!Lxu)
00191         {
00192           Lxu=new dmatrix(1, (int)num_active_parameters, mmin-1, mmax);
00193         }
00194         for (int i=1;i<=(int)num_active_parameters;i++)
00195         {
00196           (*Lxu)(i,mmin-1)=(*funnel_init_var::plist)(i,1);
00197         }
00198         for (int j=mmin;j<=mmax;j++)
00199         {
00200           for (int i=1;i<=(int)num_active_parameters;i++)
00201           {
00202             switch (old_style_flag)
00203             {
00204             case 0:
00205               (*Lxu)(i,j)=v(j).get_u_dot()[i-1];
00206               break;
00207             case 1:
00208             case 2:
00209               (*Lxu)(i,j)=2.0*v(j).get_u_dot()[i-1];
00210               break;
00211             default:
00212               cerr << "Illegal value for quadratic_prior::old_style_flag"
00213                    << endl;
00214               ad_exit(1);
00215             }
00216           }
00217         }
00218       }
00219       else
00220       {
00221         if (Lxu)
00222         {
00223           delete Lxu;
00224           Lxu=0;
00225         }
00226       }
00227     }
00228     df1b2_gradlist::no_derivatives=0;
00229   }
00230 
00231  void df1b2quadratic_prior::get_Lxu_contribution(dmatrix& M)
00232  {
00233    for (int i=0;i<num_quadratic_prior;i++)
00234    {
00235      //cout << ptr[i]->get_num_active_parameters() << endl;
00236      //if (ptr[i]->get_num_active_parameters()>0)
00237      {
00238        ptr[i]->get_Lxu(M);
00239      }
00240    }
00241  }
00242 
00243 normal_df1b2quadratic_prior::normal_df1b2quadratic_prior(void)
00244 {
00245   set_old_style_flag();
00246 }
00247 
00248 void normal_df1b2quadratic_prior::set_old_style_flag(void)
00249 {
00250   old_style_flag=0;
00251 }
00252 void normal_df1b2quadratic_prior::operator = (const df1b2matrix & M)
00253 {
00254   df1b2quadratic_prior::operator = (M);
00255 }
00256 
00257 void df1b2quadratic_re_penalty::set_old_style_flag(void)
00258 {
00259   old_style_flag=2;
00260 }
00261 void df1b2quadratic_re_penalty::operator = (const df1b2matrix & M)
00262 {
00263   df1b2quadratic_prior::operator = (M);
00264 }
00265 void df1b2quadratic_re_penalty::operator = (const dmatrix & M)
00266 {
00267   df1b2quadratic_prior::operator = (M);
00268 }
00269 
00270 
00271 df1b2quadratic_re_penalty::df1b2quadratic_re_penalty(void)
00272 {
00273   set_old_style_flag();
00274 }
00275 // *******************************************************
00276 // *******************************************************
00277 // *******************************************************
00278 // *******************************************************
00279 constant_df1b2quadratic_re_penalty::constant_df1b2quadratic_re_penalty(void)
00280 {
00281   set_old_style_flag();
00282 }
00283 
00284 void constant_df1b2quadratic_re_penalty::set_old_style_flag(void)
00285 {
00286   old_style_flag=2;
00287 }
00288 void constant_df1b2quadratic_re_penalty::operator = (const dmatrix & M)
00289 {
00290   //df1b2quadratic_prior::operator = (M);
00291 }
00292 void df1b2quadratic_prior::operator = (const dmatrix & M)
00293 {
00294   quadratic_prior::in_qp_calculations=0;
00295   num_active_parameters=funnel_init_var::num_vars;
00296   df1b2_gradlist::no_derivatives=1;
00297   dvector cu=value(*pu);
00298 
00299   if (laplace_approximation_calculator::where_are_we_flag==3)
00300   {
00301     df1b2variable::noallocate=1;
00302     //df1b2vector v(M.indexmin(),M.indexmax());
00303     df1b2variable::noallocate=0;
00304     switch (old_style_flag)
00305     {
00306     case 0:
00307     case 1:
00308       cout << "this can't happen" << endl;
00309       ad_exit(1);
00310       break;
00311     case 2:
00312       //v = M*cu;
00313       break;
00314     default:
00315       cerr << "Illegal value for quadratic_prior::old_style_flag"
00316            << endl;
00317       ad_exit(1);
00318     }
00319     //int mmin=v.indexmin();
00320     //int mmax=v.indexmax();
00321 
00322     if (index)
00323     {
00324       //if (index->indexmax() != nvar)
00325       cout << "this can't happen" << endl;
00326       ad_exit(1);
00327       delete index;
00328       index=0;
00329     }
00330 
00331     if (Lxu)
00332     {
00333       cout << "this can't happen" << endl;
00334       ad_exit(1);
00335       delete Lxu;
00336       Lxu=0;
00337     }
00338   }
00339   df1b2_gradlist::no_derivatives=0;
00340 }
00341 
00342