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