ADMB Documentation  11.1.2547
 All Classes Files Functions Variables Typedefs Friends Defines
quadpri.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: quadpri.cpp 2535 2014-10-31 12:28:33Z 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     int nv=df1b2quadratic_prior::get_ptr(i)->get_num_active_parameters();
00074     if (nv>0)
00075       quadratic_prior::get_ptr(i)->get_vHessian(Hess_all,xsize);
00076     else
00077       quadratic_prior::get_ptr(i)->get_cHessian(Hess_all,xsize);
00078   }
00079   int sgn;
00080   vf=0.5*ln_det(Hess_all,sgn);
00081   gradcalc(xsize,g);
00082   return g;
00083 }
00084 
00089   void quadratic_prior::add_to_list(void)
00090   {
00091     if (num_quadratic_prior>=max_num_quadratic_prior) ad_exit(1);
00092     xmyindex=num_quadratic_prior;
00093     ptr[num_quadratic_prior++]=this;
00094   }
00095 
00100   dvariable quadratic_prior::get_function(void)
00101   {
00102     return (*pu)*((*pMinv)*(*pu));
00103   }
00104 
00109   dvar_matrix quadratic_prior::get_Hessian(void)
00110   {
00111     return *pMinv;
00112   }
00113 
00118   int quadratic_prior::get_offset(int xs)
00119   {
00120     df1b2_init_vector * fpu=df1b2quadratic_prior::ptr[get_myindex()]->pu;
00121     int mmin=(*fpu)(fpu->indexmin()).get_ind_index();
00122     return mmin-xs-1;
00123   }
00124   //dmatrix quadratic_prior::get_cHessian(void)
00125 
00130   void quadratic_prior::get_cHessian(dmatrix H,int xsize)
00131   {
00132     int offset=get_offset(xsize);
00133     int imin=pMinv->indexmin();
00134     int imax=pMinv->indexmax();
00135     if (offset==0)
00136     {
00137       int i,j;
00138       switch(old_style_flag)
00139       {
00140       case 0:
00141         for (i=imin;i<=imax;i++)
00142           for (j=imin;j<=imax;j++)
00143             H(i,j)+=(*pMinv)(i,j);
00144         break;
00145       case 1:
00146         for (i=imin;i<=imax;i++)
00147           for (j=imin;j<=imax;j++)
00148             H(i,j)+=2.0*(*pMinv)(i,j);
00149         break;
00150       case 2:
00151         for (i=imin;i<=imax;i++)
00152           for (j=imin;j<=imax;j++)
00153             H(i,j)+=2.0*(*pMinv)(i,j);
00154          break;
00155       default:
00156          cerr << "Illegal value in switch statement" << endl;
00157          ad_exit(1);
00158       }
00159     }
00160     else
00161     {
00162       int i,j;
00163       switch(old_style_flag)
00164       {
00165       case 0:
00166         for (i=imin;i<=imax;i++)
00167           for (j=imin;j<=imax;j++)
00168             H(offset+i,offset+j)+=(*pMinv)(i,j);
00169         break;
00170       case 1:
00171         for (i=imin;i<=imax;i++)
00172           for (j=imin;j<=imax;j++)
00173             H(offset+i,offset+j)+=2.0*(*pMinv)(i,j);
00174         break;
00175       case 2:
00176         for (i=imin;i<=imax;i++)
00177           for (j=imin;j<=imax;j++)
00178             H(offset+i,offset+j)+=2.0*(*pMinv)(i,j);
00179          break;
00180       default:
00181          cerr << "Illegal value in switch statement" << endl;
00182          ad_exit(1);
00183       }
00184     }
00185   }
00186 
00191   void quadratic_prior::get_cHessian(dvar_matrix H,int xsize)
00192   {
00193     int offset=get_offset(xsize);
00194     int imin=pMinv->indexmin();
00195     int imax=pMinv->indexmax();
00196     if (offset==0)
00197     {
00198       int i,j;
00199       switch(old_style_flag)
00200       {
00201       case 0:
00202         for (i=imin;i<=imax;i++)
00203           for (j=imin;j<=imax;j++)
00204             H(i,j)+=(*pMinv)(i,j);
00205         break;
00206       case 1:
00207         for (i=imin;i<=imax;i++)
00208           for (j=imin;j<=imax;j++)
00209             H(i,j)+=2.0*(*pMinv)(i,j);
00210         break;
00211       case 2:
00212         for (i=imin;i<=imax;i++)
00213           for (j=imin;j<=imax;j++)
00214             H(i,j)+=2.0*(*pMinv)(i,j);
00215          break;
00216       default:
00217          cerr << "Illegal value in switch statement" << endl;
00218          ad_exit(1);
00219       }
00220     }
00221     else
00222     {
00223       int i,j;
00224       switch(old_style_flag)
00225       {
00226       case 0:
00227         for (i=imin;i<=imax;i++)
00228           for (j=imin;j<=imax;j++)
00229             H(offset+i,offset+j)+=(*pMinv)(i,j);
00230         break;
00231       case 1:
00232         for (i=imin;i<=imax;i++)
00233           for (j=imin;j<=imax;j++)
00234             H(offset+i,offset+j)+=2.0*(*pMinv)(i,j);
00235         break;
00236       case 2:
00237         for (i=imin;i<=imax;i++)
00238           for (j=imin;j<=imax;j++)
00239             H(offset+i,offset+j)+=2.0*(*pMinv)(i,j);
00240          break;
00241       default:
00242          cerr << "Illegal value in switch statement" << endl;
00243          ad_exit(1);
00244       }
00245     }
00246   }
00247 
00252   void quadratic_prior::get_vHessian(dvar_matrix H,int xsize)
00253   {
00254     int offset=get_offset(xsize);
00255     if (dfpMinv==0)
00256     {
00257       cerr << "This can't happen" << endl;
00258       ad_exit(1);
00259     }
00260     int imin=dfpMinv->indexmin();
00261     int imax=dfpMinv->indexmax();
00262     if (offset==0)
00263     {
00264       int i,j;
00265       switch(old_style_flag)
00266       {
00267       case 0:
00268         for (i=imin;i<=imax;i++)
00269           for (j=imin;j<=imax;j++)
00270             H(i,j)+=(*dfpMinv)(i,j);
00271         break;
00272       case 1:
00273         for (i=imin;i<=imax;i++)
00274           for (j=imin;j<=imax;j++)
00275             H(i,j)+=2.0*(*dfpMinv)(i,j);
00276         break;
00277       case 2:
00278         for (i=imin;i<=imax;i++)
00279           for (j=imin;j<=imax;j++)
00280             H(i,j)+=2.0*(*dfpMinv)(i,j);
00281          break;
00282       default:
00283          cerr << "Illegal valueinswitch statement" << endl;
00284          ad_exit(1);
00285       }
00286     }
00287     else
00288     {
00289       int i,j;
00290       switch(old_style_flag)
00291       {
00292       case 0:
00293         for (i=imin;i<=imax;i++)
00294           for (j=imin;j<=imax;j++)
00295             H(offset+i,offset+j)+=(*dfpMinv)(i,j);
00296         break;
00297       case 1:
00298         for (i=imin;i<=imax;i++)
00299           for (j=imin;j<=imax;j++)
00300             H(offset+i,offset+j)+=2.0*(*dfpMinv)(i,j);
00301         break;
00302       case 2:
00303         for (i=imin;i<=imax;i++)
00304           for (j=imin;j<=imax;j++)
00305             H(offset+i,offset+j)+=2.0*(*dfpMinv)(i,j);
00306          break;
00307       default:
00308          cerr << "Illegal valueinswitch statement" << endl;
00309          ad_exit(1);
00310       }
00311     }
00312   }
00313  /*
00314   dvar_matrix quadratic_prior::get_vHessian(void)
00315   {
00316     return *dfpMinv;
00317     //return value(*pMinv);
00318   }
00319  */
00320 
00325    void quadratic_prior::get_cHessian_from_vHessian(dmatrix H,int xs)
00326   {
00327     int offset=get_offset(xs);
00328     int imin=dfpMinv->indexmin();
00329     int imax=dfpMinv->indexmax();
00330     if (offset==0)
00331     {
00332       int i,j;
00333       switch(old_style_flag)
00334       {
00335       case 0:
00336         for (i=imin;i<=imax;i++)
00337           for (j=imin;j<=imax;j++)
00338             H(i,j)+=value((*dfpMinv)(i,j));
00339         break;
00340       case 1:
00341         for (i=imin;i<=imax;i++)
00342           for (j=imin;j<=imax;j++)
00343             H(i,j)+=2.0*value((*dfpMinv)(i,j));
00344         break;
00345       case 2:
00346         for (i=imin;i<=imax;i++)
00347           for (j=imin;j<=imax;j++)
00348             H(i,j)+=2.0*value((*dfpMinv)(i,j));
00349         break;
00350          break;
00351       default:
00352          cerr << "Illegal valueinswitch statement" << endl;
00353          ad_exit(1);
00354       }
00355     }
00356     else
00357     {
00358       int i,j;
00359       switch(old_style_flag)
00360       {
00361       case 0:
00362         for (i=imin;i<=imax;i++)
00363           for (j=imin;j<=imax;j++)
00364             H(offset+i,offset+j)+=value((*dfpMinv)(i,j));
00365         break;
00366       case 1:
00367         for (i=imin;i<=imax;i++)
00368           for (j=imin;j<=imax;j++)
00369             H(offset+i,offset+j)+=2.0*value((*dfpMinv)(i,j));
00370         break;
00371       case 2:
00372         for (i=imin;i<=imax;i++)
00373           for (j=imin;j<=imax;j++)
00374             H(offset+i,offset+j)+=2.0*value((*dfpMinv)(i,j));
00375          break;
00376       default:
00377          cerr << "Illegal valueinswitch statement" << endl;
00378          ad_exit(1);
00379       }
00380     }
00381     //return value(*dfpMinv);
00382   }
00383 
00388   dvar_vector quadratic_prior::get_gradient(void)
00389   {
00390     return ((*pMinv)*(*pu));
00391   }
00392 
00397   void quadratic_prior::get_cgradient(dvector g,int xs)
00398   {
00399     int offset=get_offset(xs);
00400     dvector tg=((*pMinv)*value(*pu));
00401     int imin=pMinv->indexmin();
00402     int imax=pMinv->indexmax();
00403     if (offset==0)
00404     {
00405       int i;
00406       switch(old_style_flag)
00407       {
00408       case 0:
00409         for (i=imin;i<=imax;i++)
00410           g(i)+=tg(i);
00411         break;
00412       case 1:
00413         for (i=imin;i<=imax;i++)
00414           g(i)+=2.0*tg(i);
00415         break;
00416       case 2:
00417         for (i=imin;i<=imax;i++)
00418           g(i)+=2.0*tg(i);
00419         break;
00420       default:
00421         cerr << "Illegal valueinswitch statement" << endl;
00422         ad_exit(1);
00423       }
00424     }
00425     else
00426     {
00427       int i;
00428       switch(old_style_flag)
00429       {
00430       case 0:
00431         for (i=imin;i<=imax;i++)
00432           g(offset+i)+=tg(i);
00433         break;
00434       case 1:
00435         for (i=imin;i<=imax;i++)
00436           g(offset+i)+=2.0*tg(i);
00437         break;
00438       case 2:
00439         for (i=imin;i<=imax;i++)
00440           g(offset+i)+=2.0*tg(i);
00441         break;
00442       default:
00443         cerr << "Illegal valueinswitch statement" << endl;
00444         ad_exit(1);
00445       }
00446     }
00447     //return ((*pMinv)*value(*pu));
00448   }
00449 
00454   quadratic_prior::quadratic_prior(void)
00455   {
00456     pMinv=0;
00457     dfpMinv=0;
00458     pu=0;
00459     add_to_list();
00460   }
00461 
00466   quadratic_prior::~quadratic_prior(void)
00467   {
00468     if (pMinv) delete pMinv;
00469     pMinv=0;
00470     if (pu) delete pu;
00471     pu=0;
00472     if (dfpMinv) delete pMinv;
00473     dfpMinv=0;
00474   }
00475 
00480   void quadratic_prior::allocate( const dvar_vector & _u,const char * s)
00481   {
00482     allocate(_u);
00483   }
00484 
00489   void quadratic_prior::allocate(const dvar_vector & _u)
00490   {
00491     if (!allocated(_u))
00492     {
00493       cerr << "You must put random effects vector before"
00494        " quadtratic prior in the TPL file" << endl;
00495       ad_exit(1);
00496     }
00497     pu = new dvar_vector((dvar_vector&)(_u));
00498   }
00499 
00504   void quadratic_prior::allocate(const dvar_matrix & _M,
00505     const dvar_vector & _u,const char * s)
00506   {
00507     allocate(_M,_u);
00508   }
00509 
00514   void quadratic_prior::allocate(const dvar_matrix & _M,
00515     const dvar_vector & _u)
00516   {
00517     pMinv =new dmatrix(value(inv(_M)));
00518     pu = new dvar_vector((dvar_vector&)(_u));
00519   }
00520 
00525  dvariable quadratic_prior::get_quadratic_priors(void)
00526  {
00527    dvariable f=0.0;
00528    for (int i=0;i<num_quadratic_prior;i++)
00529    {
00530      f+=ptr[i]->get_function();
00531    }
00532    return f;
00533  }
00534 
00539  void quadratic_prior::get_cgradient_contribution(dvector g,int xs)
00540  {
00541    for (int i=0;i<num_quadratic_prior;i++)
00542    {
00543      ptr[i]->get_cgradient(g,xs);
00544     /*
00545      if (old_style_flag)
00546      {
00547        return ptr[i]->get_cgradient();
00548      }
00549      else
00550      {
00551        return ptr[i]->get_cgradient();
00552      }
00553     */
00554    }
00555    //return f;
00556  }
00557 /*
00558  dvar_vector quadratic_prior::get_gradient_contribution(void)
00559  {
00560    for (int i=0;i<num_quadratic_prior;i++)
00561    {
00562      return ptr[i]->get_gradient();
00563    }
00564    //return f;
00565  }
00566 */
00567 
00572  void quadratic_prior::get_cHessian_contribution(dmatrix H,int xsize)
00573  {
00574    for (int i=0;i<num_quadratic_prior;i++)
00575    {
00576      if (!ptr[i])
00577      {
00578        cerr << "ptr["<<i<<"] = 0 in"
00579          " quadratic_prior::get_cHessian_contribution" << endl;
00580        ad_exit(1);
00581      }
00582      if (!ptr[i]->pMinv)
00583      {
00584        cerr << "ptr["<<i<<"]->pMinv = 0 in"
00585          " quadratic_prior::get_cHessian_contribution" << endl;
00586        ad_exit(1);
00587      }
00588      if (!allocated(*(ptr[i]->pMinv)))
00589      {
00590        cerr << "*ptr["<<i<<"] is unallocated in"
00591          " quadratic_prior::get_cHessian_contribution" << endl;
00592        ad_exit(1);
00593      }
00594      ptr[i]->get_cHessian(H,xsize);
00595     /*
00596      if (old_style_flag)
00597      {
00598        return 2.0*ptr[i]->get_cHessian();
00599      }
00600      else
00601      {
00602        return ptr[i]->get_cHessian();
00603      }
00604     */
00605    }
00606    //return f;
00607  }
00608 
00613  void quadratic_prior::get_vHessian_contribution(dvar_matrix H,int xs)
00614  {
00615    for (int i=0;i<num_quadratic_prior;i++)
00616    {
00617      ptr[i]->get_vHessian(H,xs);
00618     /*
00619      if (old_style_flag)
00620      {
00621        return 2.0*ptr[i]->get_vHessian();
00622      }
00623      else
00624      {
00625        return ptr[i]->get_vHessian();
00626      }
00627     */
00628    }
00629    //return f;
00630  }
00631  /*
00632  dvar_matrix quadratic_prior::get_Hessian_contribution(void)
00633  {
00634    for (int i=0;i<num_quadratic_prior;i++)
00635    {
00636      return ptr[i]->get_Hessian();
00637    }
00638    //return f;
00639  }
00640  */
00641 
00646 void quadratic_prior::get_cHessian_contribution_from_vHessian(dmatrix Hess,
00647   int xsize)
00648  {
00649    for (int i=0;i<num_quadratic_prior;i++)
00650    {
00651      int nv=df1b2quadratic_prior::get_ptr(i)->
00652        get_num_active_parameters();
00653      if (nv)
00654        ptr[i]->get_cHessian_from_vHessian(Hess,xsize);
00655      else
00656        ptr[i]->get_cHessian(Hess,xsize);
00657    }
00658    //return f;
00659  }
00660 
00665  void quadratic_prior::operator = (const dvar_matrix & _M)
00666  {
00667    dvariable lndet;
00668    dvariable sgn;
00669 
00670    switch (quadratic_prior::old_style_flag)
00671    {
00672    case 0:
00673      *objective_function_value::pobjfun+=0.5*(*pu)*(solve(_M,*pu,lndet,sgn));
00674      *objective_function_value::pobjfun+=0.5*lndet;
00675      //*objective_function_value::pobjfun+=0.5*(*pu)*(solve(_M,*pu));
00676      //*objective_function_value::pobjfun+=0.5*ln_det(_M);
00677      break;
00678    case 1:
00679      *objective_function_value::pobjfun+=(*pu)*(solve(_M,*pu));
00680      break;
00681    case 2:
00682      *objective_function_value::pobjfun+=(*pu) * ( _M * (*pu) );
00683      break;
00684    default:
00685      cerr << "Illegal value for quadratic_prior::old_style_flag"
00686           << endl;
00687      ad_exit(1);
00688    }
00689    if (pMinv)
00690    {
00691      delete pMinv;
00692      pMinv=0;
00693    }
00694    if (dfpMinv)
00695    {
00696      delete dfpMinv;
00697      dfpMinv=0;
00698    }
00699    switch (quadratic_prior::old_style_flag)
00700    {
00701    case 0:
00702    case 1:
00703      if (laplace_approximation_calculator::where_are_we_flag==2)
00704      {
00705        pMinv = new dmatrix(inv(value(_M)));
00706        if (pMinv==0)
00707        {
00708          cerr << "Error allocating dmatrix" << endl;
00709          ad_exit(1);
00710        }
00711      }
00712      if (laplace_approximation_calculator::where_are_we_flag==3)
00713      {
00714        dfpMinv = new dvar_matrix(inv(_M));
00715        if (dfpMinv==0)
00716        {
00717          cerr << "Error allocating dvar_matrix" << endl;
00718          ad_exit(1);
00719        }
00720      }
00721      break;
00722    case 2:
00723      if (laplace_approximation_calculator::where_are_we_flag==2)
00724      {
00725        pMinv = new dmatrix(value(_M));
00726        if (pMinv==0)
00727        {
00728          cerr << "Error allocating dmatrix" << endl;
00729          ad_exit(1);
00730        }
00731      }
00732      if (laplace_approximation_calculator::where_are_we_flag==3)
00733      {
00734        int nv =
00735          df1b2quadratic_prior::get_ptr(xmyindex)->get_num_active_parameters();
00736        //if (nv==0)
00737        if (nv!=0)
00738        {
00739          dfpMinv = new dvar_matrix(_M);
00740          if (dfpMinv==0)
00741          {
00742            cerr << "Error allocating dvar_matrix" << endl;
00743            ad_exit(1);
00744          }
00745        }
00746        else
00747        {
00748          pMinv = new dmatrix(value(_M));
00749          if (pMinv==0)
00750          {
00751            cerr << "Error allocating dmatrix" << endl;
00752            ad_exit(1);
00753          }
00754        }
00755      }
00756      break;
00757    default:
00758      cerr << "Illegal value for quadratic_prior::old_style_flag"
00759           << endl;
00760      ad_exit(1);
00761    }
00762  }
00763 
00768  void quadratic_prior::operator = (const dmatrix & _M)
00769  {
00770    dvariable lndet;
00771    dvariable sgn;
00772 
00773    switch (quadratic_prior::old_style_flag)
00774    {
00775    case 0:
00776      cerr << " can't get here " << endl;
00777      ad_exit(1);
00778      break;
00779    case 1:
00780      cerr << " can't get here " << endl;
00781      ad_exit(1);
00782      break;
00783    case 2:
00784      *objective_function_value::pobjfun+=(*pu) * ( _M * (*pu) );
00785      break;
00786    default:
00787      cerr << "Illegal value for quadratic_prior::old_style_flag"
00788           << endl;
00789      ad_exit(1);
00790    }
00791    if (pMinv)
00792    {
00793      delete pMinv;
00794      pMinv=0;
00795    }
00796    if (dfpMinv)
00797    {
00798      delete dfpMinv;
00799      dfpMinv=0;
00800    }
00801    switch (quadratic_prior::old_style_flag)
00802    {
00803    case 0:
00804    case 1:
00805      cerr << " can't get here " << endl;
00806      ad_exit(1);
00807      break;
00808    case 2:
00809      if (laplace_approximation_calculator::where_are_we_flag==2 ||
00810        laplace_approximation_calculator::where_are_we_flag==3)
00811      {
00812        pMinv = new dmatrix(_M);
00813        if (pMinv==0)
00814        {
00815          cerr << "Error allocating dmatrix" << endl;
00816          ad_exit(1);
00817        }
00818      }
00819      break;
00820    default:
00821      cerr << "Illegal value for quadratic_prior::old_style_flag"
00822           << endl;
00823      ad_exit(1);
00824    }
00825  }
00826 
00831 normal_quadratic_prior::normal_quadratic_prior(void)
00832 {
00833   set_old_style_flag();
00834 }
00835 
00840 void normal_quadratic_prior::set_old_style_flag(void)
00841 {
00842   old_style_flag=0;
00843 }
00844 
00849 void normal_quadratic_prior::operator = (const dvar_matrix & M)
00850 {
00851   quadratic_prior::operator = (M);
00852 }
00853 
00858 quadratic_re_penalty::quadratic_re_penalty(void)
00859 {
00860   set_old_style_flag();
00861 }
00862 
00867 void quadratic_re_penalty::set_old_style_flag(void)
00868 {
00869   old_style_flag=2;
00870 }
00871 
00876 void quadratic_re_penalty::operator = (const dvar_matrix & M)
00877 {
00878   quadratic_prior::operator = (M);
00879 }
00880 
00885 void quadratic_re_penalty::operator = (const dmatrix & M)
00886 {
00887   quadratic_prior::operator = (M);
00888 }
00889 
00894 constant_quadratic_re_penalty::constant_quadratic_re_penalty(void)
00895 {
00896   set_old_style_flag();
00897 }
00898 
00903 void constant_quadratic_re_penalty::set_old_style_flag(void)
00904 {
00905   old_style_flag=2;
00906 }
00907 
00912 void constant_quadratic_re_penalty::operator = (const dmatrix & M)
00913 {
00914   quadratic_prior::operator = (M);
00915 }