ADMB Documentation  11.4.2891
 All Classes Files Functions Variables Typedefs Friends Defines
quadpri.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  */
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 
00465 quadratic_prior::quadratic_prior()
00466 {
00467   pMinv = NULL;
00468   dfpMinv = NULL;
00469   pu = NULL;
00470   add_to_list();
00471 }
00475 quadratic_prior::~quadratic_prior()
00476 {
00477   if (pMinv)
00478   {
00479     delete pMinv;
00480     pMinv = NULL;
00481   }
00482   if (pu)
00483   {
00484     delete pu;
00485     pu = NULL;
00486   }
00487   if (dfpMinv)
00488   {
00489     delete dfpMinv;
00490     dfpMinv = NULL;
00491   }
00492 }
00493 
00498   void quadratic_prior::allocate( const dvar_vector & _u,const char * s)
00499   {
00500     allocate(_u);
00501   }
00502 
00507   void quadratic_prior::allocate(const dvar_vector & _u)
00508   {
00509     if (!allocated(_u))
00510     {
00511       cerr << "You must put random effects vector before"
00512        " quadtratic prior in the TPL file" << endl;
00513       ad_exit(1);
00514     }
00515     pu = new dvar_vector((dvar_vector&)(_u));
00516   }
00517 
00522   void quadratic_prior::allocate(const dvar_matrix & _M,
00523     const dvar_vector & _u,const char * s)
00524   {
00525     allocate(_M,_u);
00526   }
00527 
00532   void quadratic_prior::allocate(const dvar_matrix & _M,
00533     const dvar_vector & _u)
00534   {
00535     pMinv =new dmatrix(value(inv(_M)));
00536     pu = new dvar_vector((dvar_vector&)(_u));
00537   }
00538 
00543  dvariable quadratic_prior::get_quadratic_priors(void)
00544  {
00545    dvariable f=0.0;
00546    for (int i=0;i<num_quadratic_prior;i++)
00547    {
00548      f+=ptr[i]->get_function();
00549    }
00550    return f;
00551  }
00552 
00557  void quadratic_prior::get_cgradient_contribution(dvector g,int xs)
00558  {
00559    for (int i=0;i<num_quadratic_prior;i++)
00560    {
00561      ptr[i]->get_cgradient(g,xs);
00562     /*
00563      if (old_style_flag)
00564      {
00565        return ptr[i]->get_cgradient();
00566      }
00567      else
00568      {
00569        return ptr[i]->get_cgradient();
00570      }
00571     */
00572    }
00573    //return f;
00574  }
00575 /*
00576  dvar_vector quadratic_prior::get_gradient_contribution(void)
00577  {
00578    for (int i=0;i<num_quadratic_prior;i++)
00579    {
00580      return ptr[i]->get_gradient();
00581    }
00582    //return f;
00583  }
00584 */
00585 
00590 void quadratic_prior::get_cHessian_contribution(dmatrix H,int xsize)
00591 {
00592   for (int i=0;i<num_quadratic_prior;i++)
00593   {
00594     if (!ptr[i])
00595     {
00596        cerr << "ptr["<<i<<"] = 0 in"
00597          " quadratic_prior::get_cHessian_contribution" << endl;
00598        ad_exit(1);
00599     }
00600     else if (!ptr[i]->pMinv)
00601     {
00602        cerr << "ptr["<<i<<"]->pMinv = 0 in"
00603          " quadratic_prior::get_cHessian_contribution" << endl;
00604        ad_exit(1);
00605     }
00606     else if (!allocated(*(ptr[i]->pMinv)))
00607     {
00608        cerr << "*ptr["<<i<<"] is unallocated in"
00609          " quadratic_prior::get_cHessian_contribution" << endl;
00610        ad_exit(1);
00611     }
00612     else
00613     {
00614       ptr[i]->get_cHessian(H,xsize);
00615       /*
00616       if (old_style_flag)
00617       {
00618         return 2.0*ptr[i]->get_cHessian();
00619       }
00620       else
00621       {
00622         return ptr[i]->get_cHessian();
00623       }
00624       */
00625     }
00626   }
00627   //return f;
00628 }
00629 
00634  void quadratic_prior::get_vHessian_contribution(dvar_matrix H,int xs)
00635  {
00636    for (int i=0;i<num_quadratic_prior;i++)
00637    {
00638      ptr[i]->get_vHessian(H,xs);
00639     /*
00640      if (old_style_flag)
00641      {
00642        return 2.0*ptr[i]->get_vHessian();
00643      }
00644      else
00645      {
00646        return ptr[i]->get_vHessian();
00647      }
00648     */
00649    }
00650    //return f;
00651  }
00652  /*
00653  dvar_matrix quadratic_prior::get_Hessian_contribution(void)
00654  {
00655    for (int i=0;i<num_quadratic_prior;i++)
00656    {
00657      return ptr[i]->get_Hessian();
00658    }
00659    //return f;
00660  }
00661  */
00662 
00667 void quadratic_prior::get_cHessian_contribution_from_vHessian(dmatrix Hess,
00668   int xsize)
00669  {
00670    for (int i=0;i<num_quadratic_prior;i++)
00671    {
00672      unsigned int nv=df1b2quadratic_prior::get_ptr(i)->
00673        get_num_active_parameters();
00674      if (nv)
00675        ptr[i]->get_cHessian_from_vHessian(Hess,xsize);
00676      else
00677        ptr[i]->get_cHessian(Hess,xsize);
00678    }
00679    //return f;
00680  }
00681 
00686  void quadratic_prior::operator = (const dvar_matrix & _M)
00687  {
00688    dvariable lndet;
00689    dvariable sgn;
00690 
00691    switch (quadratic_prior::old_style_flag)
00692    {
00693    case 0:
00694      *objective_function_value::pobjfun+=0.5*(*pu)*(solve(_M,*pu,lndet,sgn));
00695      *objective_function_value::pobjfun+=0.5*lndet;
00696      //*objective_function_value::pobjfun+=0.5*(*pu)*(solve(_M,*pu));
00697      //*objective_function_value::pobjfun+=0.5*ln_det(_M);
00698      break;
00699    case 1:
00700      *objective_function_value::pobjfun+=(*pu)*(solve(_M,*pu));
00701      break;
00702    case 2:
00703      *objective_function_value::pobjfun+=(*pu) * ( _M * (*pu) );
00704      break;
00705    default:
00706      cerr << "Illegal value for quadratic_prior::old_style_flag"
00707           << endl;
00708      ad_exit(1);
00709    }
00710    if (pMinv)
00711    {
00712      delete pMinv;
00713      pMinv=0;
00714    }
00715    if (dfpMinv)
00716    {
00717      delete dfpMinv;
00718      dfpMinv=0;
00719    }
00720    switch (quadratic_prior::old_style_flag)
00721    {
00722    case 0:
00723    case 1:
00724      if (laplace_approximation_calculator::where_are_we_flag==2)
00725      {
00726        pMinv = new dmatrix(inv(value(_M)));
00727        if (pMinv==0)
00728        {
00729          cerr << "Error allocating dmatrix" << endl;
00730          ad_exit(1);
00731        }
00732      }
00733      if (laplace_approximation_calculator::where_are_we_flag==3)
00734      {
00735        dfpMinv = new dvar_matrix(inv(_M));
00736        if (dfpMinv==0)
00737        {
00738          cerr << "Error allocating dvar_matrix" << endl;
00739          ad_exit(1);
00740        }
00741      }
00742      break;
00743    case 2:
00744      if (laplace_approximation_calculator::where_are_we_flag==2)
00745      {
00746        pMinv = new dmatrix(value(_M));
00747        if (pMinv==0)
00748        {
00749          cerr << "Error allocating dmatrix" << endl;
00750          ad_exit(1);
00751        }
00752      }
00753      if (laplace_approximation_calculator::where_are_we_flag==3)
00754      {
00755        unsigned int nv =
00756          df1b2quadratic_prior::get_ptr(xmyindex)->get_num_active_parameters();
00757        //if (nv==0)
00758        if (nv!=0)
00759        {
00760          dfpMinv = new dvar_matrix(_M);
00761          if (dfpMinv==0)
00762          {
00763            cerr << "Error allocating dvar_matrix" << endl;
00764            ad_exit(1);
00765          }
00766        }
00767        else
00768        {
00769          pMinv = new dmatrix(value(_M));
00770          if (pMinv==0)
00771          {
00772            cerr << "Error allocating dmatrix" << endl;
00773            ad_exit(1);
00774          }
00775        }
00776      }
00777      break;
00778    default:
00779      cerr << "Illegal value for quadratic_prior::old_style_flag"
00780           << endl;
00781      ad_exit(1);
00782    }
00783  }
00784 
00789  void quadratic_prior::operator = (const dmatrix & _M)
00790  {
00791    dvariable lndet;
00792    dvariable sgn;
00793 
00794    switch (quadratic_prior::old_style_flag)
00795    {
00796    case 0:
00797      cerr << " can't get here " << endl;
00798      ad_exit(1);
00799      break;
00800    case 1:
00801      cerr << " can't get here " << endl;
00802      ad_exit(1);
00803      break;
00804    case 2:
00805      *objective_function_value::pobjfun+=(*pu) * ( _M * (*pu) );
00806      break;
00807    default:
00808      cerr << "Illegal value for quadratic_prior::old_style_flag"
00809           << endl;
00810      ad_exit(1);
00811    }
00812    if (pMinv)
00813    {
00814      delete pMinv;
00815      pMinv=0;
00816    }
00817    if (dfpMinv)
00818    {
00819      delete dfpMinv;
00820      dfpMinv=0;
00821    }
00822    switch (quadratic_prior::old_style_flag)
00823    {
00824    case 0:
00825    case 1:
00826      cerr << " can't get here " << endl;
00827      ad_exit(1);
00828      break;
00829    case 2:
00830      if (laplace_approximation_calculator::where_are_we_flag==2 ||
00831        laplace_approximation_calculator::where_are_we_flag==3)
00832      {
00833        pMinv = new dmatrix(_M);
00834        if (pMinv==0)
00835        {
00836          cerr << "Error allocating dmatrix" << endl;
00837          ad_exit(1);
00838        }
00839      }
00840      break;
00841    default:
00842      cerr << "Illegal value for quadratic_prior::old_style_flag"
00843           << endl;
00844      ad_exit(1);
00845    }
00846  }
00847 
00852 normal_quadratic_prior::normal_quadratic_prior(void)
00853 {
00854   set_old_style_flag();
00855 }
00856 
00861 void normal_quadratic_prior::set_old_style_flag(void)
00862 {
00863   old_style_flag=0;
00864 }
00865 
00870 void normal_quadratic_prior::operator = (const dvar_matrix & M)
00871 {
00872   quadratic_prior::operator = (M);
00873 }
00874 
00879 quadratic_re_penalty::quadratic_re_penalty(void)
00880 {
00881   set_old_style_flag();
00882 }
00883 
00888 void quadratic_re_penalty::set_old_style_flag(void)
00889 {
00890   old_style_flag=2;
00891 }
00892 
00897 void quadratic_re_penalty::operator = (const dvar_matrix & M)
00898 {
00899   quadratic_prior::operator = (M);
00900 }
00901 
00906 void quadratic_re_penalty::operator = (const dmatrix & M)
00907 {
00908   quadratic_prior::operator = (M);
00909 }
00910 
00915 constant_quadratic_re_penalty::constant_quadratic_re_penalty(void)
00916 {
00917   set_old_style_flag();
00918 }
00919 
00924 void constant_quadratic_re_penalty::set_old_style_flag(void)
00925 {
00926   old_style_flag=2;
00927 }
00928 
00933 void constant_quadratic_re_penalty::operator = (const dmatrix & M)
00934 {
00935   quadratic_prior::operator = (M);
00936 }