ADMB Documentation  11.2.2853
 All Classes Files Functions Variables Typedefs Friends Defines
quadpri.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: quadpri.cpp 2789 2014-12-11 20:05:11Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #include <df1b2fnl.h>
00012 
00013 #ifndef OPT_LIB
00014   #include <cassert>
00015   #include <climits>
00016 #endif
00017 
00018 int quadratic_prior::in_qp_calculations=0;
00019 
00020 // this should be a resizeable array
00021 quadratic_prior * quadratic_prior::ptr[100];
00022 
00023 int quadratic_prior::num_quadratic_prior=0;
00024 const int quadratic_prior::max_num_quadratic_prior=100;
00025 
00030  void quadratic_prior::get_M_calculations(void)
00031  {
00032    for (int i=0;i<num_quadratic_prior;i++)
00033    {
00034      //if (ptr[i]->get_num_active_parameters()>0)
00035      {
00036        ptr[i]->get_cM();
00037      }
00038    }
00039  }
00040 
00045 dvector evaluate_function_with_quadprior(const dvector& x,int usize,
00046   function_minimizer * pfmin)
00047 {
00048   int xsize = initial_params::nvarcalc();
00049   dvector g(1,xsize);
00050   gradcalc(0,g);
00051   //double f=0.0;
00052   independent_variables u(1,xsize);
00053   u=x;
00054   dvariable vf=0.0;
00055   initial_params::reset(dvar_vector(u));
00056   //vf=0.0;
00057   dvar_matrix Hess_all(1,usize,1,usize);
00058   *objective_function_value::pobjfun=0.0;
00059   // so that dvar_matrix Hessian contributions are calculated
00060   laplace_approximation_calculator::where_are_we_flag=3;
00061   pfmin->AD_uf_inner();
00062   if ( quadratic_prior::get_num_quadratic_prior()>0)
00063   {
00064     quadratic_prior::get_M_calculations();
00065   }
00066   laplace_approximation_calculator::where_are_we_flag=0;
00067 
00068   *objective_function_value::pobjfun=0.0;
00069   Hess_all=pfmin->lapprox->Hess;
00070   for (int i=0;i<quadratic_prior::get_num_quadratic_prior();i++)
00071   {
00072     //Hess_all += quadratic_prior::get_ptr(i)->get_vHessian();
00073     unsigned int nv =
00074       df1b2quadratic_prior::get_ptr(i)->get_num_active_parameters();
00075     if (nv>0)
00076       quadratic_prior::get_ptr(i)->get_vHessian(Hess_all,xsize);
00077     else
00078       quadratic_prior::get_ptr(i)->get_cHessian(Hess_all,xsize);
00079   }
00080   int sgn;
00081   vf=0.5*ln_det(Hess_all,sgn);
00082   gradcalc(xsize,g);
00083   return g;
00084 }
00085 
00090 void quadratic_prior::add_to_list(void)
00091 {
00092   if (num_quadratic_prior>=max_num_quadratic_prior)
00093   {
00094     cerr << "Error[" << __FILE__ << ':' << __LINE__
00095          << "]: Max size exceeded.\n";
00096 
00097     ad_exit(1);
00098   }
00099   else
00100   {
00101     xmyindex=num_quadratic_prior;
00102     ptr[num_quadratic_prior++]=this;
00103   }
00104 }
00105 
00110   dvariable quadratic_prior::get_function(void)
00111   {
00112     return (*pu)*((*pMinv)*(*pu));
00113   }
00114 
00119   dvar_matrix quadratic_prior::get_Hessian(void)
00120   {
00121     return *pMinv;
00122   }
00123 
00128   int quadratic_prior::get_offset(int xs)
00129   {
00130     df1b2_init_vector * fpu=df1b2quadratic_prior::ptr[get_myindex()]->pu;
00131     int mmin=(*fpu)(fpu->indexmin()).get_ind_index();
00132     return mmin-xs-1;
00133   }
00134   //dmatrix quadratic_prior::get_cHessian(void)
00135 
00140   void quadratic_prior::get_cHessian(dmatrix H,int xsize)
00141   {
00142     int offset=get_offset(xsize);
00143     int imin=pMinv->indexmin();
00144     int imax=pMinv->indexmax();
00145     if (offset==0)
00146     {
00147       int i,j;
00148       switch(old_style_flag)
00149       {
00150       case 0:
00151         for (i=imin;i<=imax;i++)
00152           for (j=imin;j<=imax;j++)
00153             H(i,j)+=(*pMinv)(i,j);
00154         break;
00155       case 1:
00156         for (i=imin;i<=imax;i++)
00157           for (j=imin;j<=imax;j++)
00158             H(i,j)+=2.0*(*pMinv)(i,j);
00159         break;
00160       case 2:
00161         for (i=imin;i<=imax;i++)
00162           for (j=imin;j<=imax;j++)
00163             H(i,j)+=2.0*(*pMinv)(i,j);
00164          break;
00165       default:
00166          cerr << "Illegal value in switch statement" << endl;
00167          ad_exit(1);
00168       }
00169     }
00170     else
00171     {
00172       int i,j;
00173       switch(old_style_flag)
00174       {
00175       case 0:
00176         for (i=imin;i<=imax;i++)
00177           for (j=imin;j<=imax;j++)
00178             H(offset+i,offset+j)+=(*pMinv)(i,j);
00179         break;
00180       case 1:
00181         for (i=imin;i<=imax;i++)
00182           for (j=imin;j<=imax;j++)
00183             H(offset+i,offset+j)+=2.0*(*pMinv)(i,j);
00184         break;
00185       case 2:
00186         for (i=imin;i<=imax;i++)
00187           for (j=imin;j<=imax;j++)
00188             H(offset+i,offset+j)+=2.0*(*pMinv)(i,j);
00189          break;
00190       default:
00191          cerr << "Illegal value in switch statement" << endl;
00192          ad_exit(1);
00193       }
00194     }
00195   }
00196 
00201   void quadratic_prior::get_cHessian(dvar_matrix H,int xsize)
00202   {
00203     int offset=get_offset(xsize);
00204     int imin=pMinv->indexmin();
00205     int imax=pMinv->indexmax();
00206     if (offset==0)
00207     {
00208       int i,j;
00209       switch(old_style_flag)
00210       {
00211       case 0:
00212         for (i=imin;i<=imax;i++)
00213           for (j=imin;j<=imax;j++)
00214             H(i,j)+=(*pMinv)(i,j);
00215         break;
00216       case 1:
00217         for (i=imin;i<=imax;i++)
00218           for (j=imin;j<=imax;j++)
00219             H(i,j)+=2.0*(*pMinv)(i,j);
00220         break;
00221       case 2:
00222         for (i=imin;i<=imax;i++)
00223           for (j=imin;j<=imax;j++)
00224             H(i,j)+=2.0*(*pMinv)(i,j);
00225          break;
00226       default:
00227          cerr << "Illegal value in switch statement" << endl;
00228          ad_exit(1);
00229       }
00230     }
00231     else
00232     {
00233       int i,j;
00234       switch(old_style_flag)
00235       {
00236       case 0:
00237         for (i=imin;i<=imax;i++)
00238           for (j=imin;j<=imax;j++)
00239             H(offset+i,offset+j)+=(*pMinv)(i,j);
00240         break;
00241       case 1:
00242         for (i=imin;i<=imax;i++)
00243           for (j=imin;j<=imax;j++)
00244             H(offset+i,offset+j)+=2.0*(*pMinv)(i,j);
00245         break;
00246       case 2:
00247         for (i=imin;i<=imax;i++)
00248           for (j=imin;j<=imax;j++)
00249             H(offset+i,offset+j)+=2.0*(*pMinv)(i,j);
00250          break;
00251       default:
00252          cerr << "Illegal value in switch statement" << endl;
00253          ad_exit(1);
00254       }
00255     }
00256   }
00257 
00262 void quadratic_prior::get_vHessian(dvar_matrix H,int xsize)
00263 {
00264   if (!dfpMinv)
00265   {
00266       cerr << "This can't happen" << endl;
00267       ad_exit(1);
00268   }
00269   else
00270   {
00271     int imin=dfpMinv->indexmin();
00272     int imax=dfpMinv->indexmax();
00273     int offset=get_offset(xsize);
00274     if (offset==0)
00275     {
00276       switch(old_style_flag)
00277       {
00278       case 0:
00279         for (int i=imin;i<=imax;i++)
00280           for (int j=imin;j<=imax;j++)
00281             H(i,j)+=(*dfpMinv)(i,j);
00282         break;
00283       case 1:
00284         for (int i=imin;i<=imax;i++)
00285           for (int j=imin;j<=imax;j++)
00286             H(i,j)+=2.0*(*dfpMinv)(i,j);
00287         break;
00288       case 2:
00289         for (int i=imin;i<=imax;i++)
00290           for (int j=imin;j<=imax;j++)
00291             H(i,j)+=2.0*(*dfpMinv)(i,j);
00292          break;
00293       default:
00294          cerr << "Illegal valueinswitch statement" << endl;
00295          ad_exit(1);
00296       }
00297     }
00298     else
00299     {
00300       switch(old_style_flag)
00301       {
00302       case 0:
00303         for (int i=imin;i<=imax;i++)
00304           for (int j=imin;j<=imax;j++)
00305             H(offset+i,offset+j)+=(*dfpMinv)(i,j);
00306         break;
00307       case 1:
00308         for (int i=imin;i<=imax;i++)
00309           for (int j=imin;j<=imax;j++)
00310             H(offset+i,offset+j)+=2.0*(*dfpMinv)(i,j);
00311         break;
00312       case 2:
00313         for (int i=imin;i<=imax;i++)
00314           for (int j=imin;j<=imax;j++)
00315             H(offset+i,offset+j)+=2.0*(*dfpMinv)(i,j);
00316          break;
00317       default:
00318          cerr << "Illegal valueinswitch statement" << endl;
00319          ad_exit(1);
00320       }
00321     }
00322   }
00323 }
00324 
00325  /*
00326   dvar_matrix quadratic_prior::get_vHessian(void)
00327   {
00328     return *dfpMinv;
00329     //return value(*pMinv);
00330   }
00331  */
00332 
00337    void quadratic_prior::get_cHessian_from_vHessian(dmatrix H,int xs)
00338   {
00339     int offset=get_offset(xs);
00340     int imin=dfpMinv->indexmin();
00341     int imax=dfpMinv->indexmax();
00342     if (offset==0)
00343     {
00344       int i,j;
00345       switch(old_style_flag)
00346       {
00347       case 0:
00348         for (i=imin;i<=imax;i++)
00349           for (j=imin;j<=imax;j++)
00350             H(i,j)+=value((*dfpMinv)(i,j));
00351         break;
00352       case 1:
00353         for (i=imin;i<=imax;i++)
00354           for (j=imin;j<=imax;j++)
00355             H(i,j)+=2.0*value((*dfpMinv)(i,j));
00356         break;
00357       case 2:
00358         for (i=imin;i<=imax;i++)
00359           for (j=imin;j<=imax;j++)
00360             H(i,j)+=2.0*value((*dfpMinv)(i,j));
00361         break;
00362          break;
00363       default:
00364          cerr << "Illegal valueinswitch statement" << endl;
00365          ad_exit(1);
00366       }
00367     }
00368     else
00369     {
00370       int i,j;
00371       switch(old_style_flag)
00372       {
00373       case 0:
00374         for (i=imin;i<=imax;i++)
00375           for (j=imin;j<=imax;j++)
00376             H(offset+i,offset+j)+=value((*dfpMinv)(i,j));
00377         break;
00378       case 1:
00379         for (i=imin;i<=imax;i++)
00380           for (j=imin;j<=imax;j++)
00381             H(offset+i,offset+j)+=2.0*value((*dfpMinv)(i,j));
00382         break;
00383       case 2:
00384         for (i=imin;i<=imax;i++)
00385           for (j=imin;j<=imax;j++)
00386             H(offset+i,offset+j)+=2.0*value((*dfpMinv)(i,j));
00387          break;
00388       default:
00389          cerr << "Illegal valueinswitch statement" << endl;
00390          ad_exit(1);
00391       }
00392     }
00393     //return value(*dfpMinv);
00394   }
00395 
00400   dvar_vector quadratic_prior::get_gradient(void)
00401   {
00402     return ((*pMinv)*(*pu));
00403   }
00404 
00409   void quadratic_prior::get_cgradient(dvector g,int xs)
00410   {
00411     int offset=get_offset(xs);
00412     dvector tg=((*pMinv)*value(*pu));
00413     int imin=pMinv->indexmin();
00414     int imax=pMinv->indexmax();
00415     if (offset==0)
00416     {
00417       int i;
00418       switch(old_style_flag)
00419       {
00420       case 0:
00421         for (i=imin;i<=imax;i++)
00422           g(i)+=tg(i);
00423         break;
00424       case 1:
00425         for (i=imin;i<=imax;i++)
00426           g(i)+=2.0*tg(i);
00427         break;
00428       case 2:
00429         for (i=imin;i<=imax;i++)
00430           g(i)+=2.0*tg(i);
00431         break;
00432       default:
00433         cerr << "Illegal valueinswitch statement" << endl;
00434         ad_exit(1);
00435       }
00436     }
00437     else
00438     {
00439       int i;
00440       switch(old_style_flag)
00441       {
00442       case 0:
00443         for (i=imin;i<=imax;i++)
00444           g(offset+i)+=tg(i);
00445         break;
00446       case 1:
00447         for (i=imin;i<=imax;i++)
00448           g(offset+i)+=2.0*tg(i);
00449         break;
00450       case 2:
00451         for (i=imin;i<=imax;i++)
00452           g(offset+i)+=2.0*tg(i);
00453         break;
00454       default:
00455         cerr << "Illegal valueinswitch statement" << endl;
00456         ad_exit(1);
00457       }
00458     }
00459     //return ((*pMinv)*value(*pu));
00460   }
00461 
00466   quadratic_prior::quadratic_prior(void)
00467   {
00468     pMinv=0;
00469     dfpMinv=0;
00470     pu=0;
00471     add_to_list();
00472   }
00473 
00478   quadratic_prior::~quadratic_prior(void)
00479   {
00480     if (pMinv) delete pMinv;
00481     pMinv=0;
00482     if (pu) delete pu;
00483     pu=0;
00484     if (dfpMinv) delete pMinv;
00485     dfpMinv=0;
00486   }
00487 
00492   void quadratic_prior::allocate( const dvar_vector & _u,const char * s)
00493   {
00494     allocate(_u);
00495   }
00496 
00501   void quadratic_prior::allocate(const dvar_vector & _u)
00502   {
00503     if (!allocated(_u))
00504     {
00505       cerr << "You must put random effects vector before"
00506        " quadtratic prior in the TPL file" << endl;
00507       ad_exit(1);
00508     }
00509     pu = new dvar_vector((dvar_vector&)(_u));
00510   }
00511 
00516   void quadratic_prior::allocate(const dvar_matrix & _M,
00517     const dvar_vector & _u,const char * s)
00518   {
00519     allocate(_M,_u);
00520   }
00521 
00526   void quadratic_prior::allocate(const dvar_matrix & _M,
00527     const dvar_vector & _u)
00528   {
00529     pMinv =new dmatrix(value(inv(_M)));
00530     pu = new dvar_vector((dvar_vector&)(_u));
00531   }
00532 
00537  dvariable quadratic_prior::get_quadratic_priors(void)
00538  {
00539    dvariable f=0.0;
00540    for (int i=0;i<num_quadratic_prior;i++)
00541    {
00542      f+=ptr[i]->get_function();
00543    }
00544    return f;
00545  }
00546 
00551  void quadratic_prior::get_cgradient_contribution(dvector g,int xs)
00552  {
00553    for (int i=0;i<num_quadratic_prior;i++)
00554    {
00555      ptr[i]->get_cgradient(g,xs);
00556     /*
00557      if (old_style_flag)
00558      {
00559        return ptr[i]->get_cgradient();
00560      }
00561      else
00562      {
00563        return ptr[i]->get_cgradient();
00564      }
00565     */
00566    }
00567    //return f;
00568  }
00569 /*
00570  dvar_vector quadratic_prior::get_gradient_contribution(void)
00571  {
00572    for (int i=0;i<num_quadratic_prior;i++)
00573    {
00574      return ptr[i]->get_gradient();
00575    }
00576    //return f;
00577  }
00578 */
00579 
00584 void quadratic_prior::get_cHessian_contribution(dmatrix H,int xsize)
00585 {
00586   for (int i=0;i<num_quadratic_prior;i++)
00587   {
00588     if (!ptr[i])
00589     {
00590        cerr << "ptr["<<i<<"] = 0 in"
00591          " quadratic_prior::get_cHessian_contribution" << endl;
00592        ad_exit(1);
00593     }
00594     else if (!ptr[i]->pMinv)
00595     {
00596        cerr << "ptr["<<i<<"]->pMinv = 0 in"
00597          " quadratic_prior::get_cHessian_contribution" << endl;
00598        ad_exit(1);
00599     }
00600     else if (!allocated(*(ptr[i]->pMinv)))
00601     {
00602        cerr << "*ptr["<<i<<"] is unallocated in"
00603          " quadratic_prior::get_cHessian_contribution" << endl;
00604        ad_exit(1);
00605     }
00606     else
00607     {
00608       ptr[i]->get_cHessian(H,xsize);
00609       /*
00610       if (old_style_flag)
00611       {
00612         return 2.0*ptr[i]->get_cHessian();
00613       }
00614       else
00615       {
00616         return ptr[i]->get_cHessian();
00617       }
00618       */
00619     }
00620   }
00621   //return f;
00622 }
00623 
00628  void quadratic_prior::get_vHessian_contribution(dvar_matrix H,int xs)
00629  {
00630    for (int i=0;i<num_quadratic_prior;i++)
00631    {
00632      ptr[i]->get_vHessian(H,xs);
00633     /*
00634      if (old_style_flag)
00635      {
00636        return 2.0*ptr[i]->get_vHessian();
00637      }
00638      else
00639      {
00640        return ptr[i]->get_vHessian();
00641      }
00642     */
00643    }
00644    //return f;
00645  }
00646  /*
00647  dvar_matrix quadratic_prior::get_Hessian_contribution(void)
00648  {
00649    for (int i=0;i<num_quadratic_prior;i++)
00650    {
00651      return ptr[i]->get_Hessian();
00652    }
00653    //return f;
00654  }
00655  */
00656 
00661 void quadratic_prior::get_cHessian_contribution_from_vHessian(dmatrix Hess,
00662   int xsize)
00663  {
00664    for (int i=0;i<num_quadratic_prior;i++)
00665    {
00666      unsigned int nv=df1b2quadratic_prior::get_ptr(i)->
00667        get_num_active_parameters();
00668      if (nv)
00669        ptr[i]->get_cHessian_from_vHessian(Hess,xsize);
00670      else
00671        ptr[i]->get_cHessian(Hess,xsize);
00672    }
00673    //return f;
00674  }
00675 
00680  void quadratic_prior::operator = (const dvar_matrix & _M)
00681  {
00682    dvariable lndet;
00683    dvariable sgn;
00684 
00685    switch (quadratic_prior::old_style_flag)
00686    {
00687    case 0:
00688      *objective_function_value::pobjfun+=0.5*(*pu)*(solve(_M,*pu,lndet,sgn));
00689      *objective_function_value::pobjfun+=0.5*lndet;
00690      //*objective_function_value::pobjfun+=0.5*(*pu)*(solve(_M,*pu));
00691      //*objective_function_value::pobjfun+=0.5*ln_det(_M);
00692      break;
00693    case 1:
00694      *objective_function_value::pobjfun+=(*pu)*(solve(_M,*pu));
00695      break;
00696    case 2:
00697      *objective_function_value::pobjfun+=(*pu) * ( _M * (*pu) );
00698      break;
00699    default:
00700      cerr << "Illegal value for quadratic_prior::old_style_flag"
00701           << endl;
00702      ad_exit(1);
00703    }
00704    if (pMinv)
00705    {
00706      delete pMinv;
00707      pMinv=0;
00708    }
00709    if (dfpMinv)
00710    {
00711      delete dfpMinv;
00712      dfpMinv=0;
00713    }
00714    switch (quadratic_prior::old_style_flag)
00715    {
00716    case 0:
00717    case 1:
00718      if (laplace_approximation_calculator::where_are_we_flag==2)
00719      {
00720        pMinv = new dmatrix(inv(value(_M)));
00721        if (pMinv==0)
00722        {
00723          cerr << "Error allocating dmatrix" << endl;
00724          ad_exit(1);
00725        }
00726      }
00727      if (laplace_approximation_calculator::where_are_we_flag==3)
00728      {
00729        dfpMinv = new dvar_matrix(inv(_M));
00730        if (dfpMinv==0)
00731        {
00732          cerr << "Error allocating dvar_matrix" << endl;
00733          ad_exit(1);
00734        }
00735      }
00736      break;
00737    case 2:
00738      if (laplace_approximation_calculator::where_are_we_flag==2)
00739      {
00740        pMinv = new dmatrix(value(_M));
00741        if (pMinv==0)
00742        {
00743          cerr << "Error allocating dmatrix" << endl;
00744          ad_exit(1);
00745        }
00746      }
00747      if (laplace_approximation_calculator::where_are_we_flag==3)
00748      {
00749        unsigned int nv =
00750          df1b2quadratic_prior::get_ptr(xmyindex)->get_num_active_parameters();
00751        //if (nv==0)
00752        if (nv!=0)
00753        {
00754          dfpMinv = new dvar_matrix(_M);
00755          if (dfpMinv==0)
00756          {
00757            cerr << "Error allocating dvar_matrix" << endl;
00758            ad_exit(1);
00759          }
00760        }
00761        else
00762        {
00763          pMinv = new dmatrix(value(_M));
00764          if (pMinv==0)
00765          {
00766            cerr << "Error allocating dmatrix" << endl;
00767            ad_exit(1);
00768          }
00769        }
00770      }
00771      break;
00772    default:
00773      cerr << "Illegal value for quadratic_prior::old_style_flag"
00774           << endl;
00775      ad_exit(1);
00776    }
00777  }
00778 
00783  void quadratic_prior::operator = (const dmatrix & _M)
00784  {
00785    dvariable lndet;
00786    dvariable sgn;
00787 
00788    switch (quadratic_prior::old_style_flag)
00789    {
00790    case 0:
00791      cerr << " can't get here " << endl;
00792      ad_exit(1);
00793      break;
00794    case 1:
00795      cerr << " can't get here " << endl;
00796      ad_exit(1);
00797      break;
00798    case 2:
00799      *objective_function_value::pobjfun+=(*pu) * ( _M * (*pu) );
00800      break;
00801    default:
00802      cerr << "Illegal value for quadratic_prior::old_style_flag"
00803           << endl;
00804      ad_exit(1);
00805    }
00806    if (pMinv)
00807    {
00808      delete pMinv;
00809      pMinv=0;
00810    }
00811    if (dfpMinv)
00812    {
00813      delete dfpMinv;
00814      dfpMinv=0;
00815    }
00816    switch (quadratic_prior::old_style_flag)
00817    {
00818    case 0:
00819    case 1:
00820      cerr << " can't get here " << endl;
00821      ad_exit(1);
00822      break;
00823    case 2:
00824      if (laplace_approximation_calculator::where_are_we_flag==2 ||
00825        laplace_approximation_calculator::where_are_we_flag==3)
00826      {
00827        pMinv = new dmatrix(_M);
00828        if (pMinv==0)
00829        {
00830          cerr << "Error allocating dmatrix" << endl;
00831          ad_exit(1);
00832        }
00833      }
00834      break;
00835    default:
00836      cerr << "Illegal value for quadratic_prior::old_style_flag"
00837           << endl;
00838      ad_exit(1);
00839    }
00840  }
00841 
00846 normal_quadratic_prior::normal_quadratic_prior(void)
00847 {
00848   set_old_style_flag();
00849 }
00850 
00855 void normal_quadratic_prior::set_old_style_flag(void)
00856 {
00857   old_style_flag=0;
00858 }
00859 
00864 void normal_quadratic_prior::operator = (const dvar_matrix & M)
00865 {
00866   quadratic_prior::operator = (M);
00867 }
00868 
00873 quadratic_re_penalty::quadratic_re_penalty(void)
00874 {
00875   set_old_style_flag();
00876 }
00877 
00882 void quadratic_re_penalty::set_old_style_flag(void)
00883 {
00884   old_style_flag=2;
00885 }
00886 
00891 void quadratic_re_penalty::operator = (const dvar_matrix & M)
00892 {
00893   quadratic_prior::operator = (M);
00894 }
00895 
00900 void quadratic_re_penalty::operator = (const dmatrix & M)
00901 {
00902   quadratic_prior::operator = (M);
00903 }
00904 
00909 constant_quadratic_re_penalty::constant_quadratic_re_penalty(void)
00910 {
00911   set_old_style_flag();
00912 }
00913 
00918 void constant_quadratic_re_penalty::set_old_style_flag(void)
00919 {
00920   old_style_flag=2;
00921 }
00922 
00927 void constant_quadratic_re_penalty::operator = (const dmatrix & M)
00928 {
00929   quadratic_prior::operator = (M);
00930 }