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