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