ADMB Documentation  11.1x.2717
 All Classes Files Functions Variables Typedefs Friends Defines
quadpri.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: quadpri.cpp 2602 2014-11-10 18:14:01Z 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) ad_exit(1);
00093     xmyindex=num_quadratic_prior;
00094     ptr[num_quadratic_prior++]=this;
00095   }
00096 
00101   dvariable quadratic_prior::get_function(void)
00102   {
00103     return (*pu)*((*pMinv)*(*pu));
00104   }
00105 
00110   dvar_matrix quadratic_prior::get_Hessian(void)
00111   {
00112     return *pMinv;
00113   }
00114 
00119   int quadratic_prior::get_offset(int xs)
00120   {
00121     df1b2_init_vector * fpu=df1b2quadratic_prior::ptr[get_myindex()]->pu;
00122     int mmin=(*fpu)(fpu->indexmin()).get_ind_index();
00123     return mmin-xs-1;
00124   }
00125   //dmatrix quadratic_prior::get_cHessian(void)
00126 
00131   void quadratic_prior::get_cHessian(dmatrix H,int xsize)
00132   {
00133     int offset=get_offset(xsize);
00134     int imin=pMinv->indexmin();
00135     int imax=pMinv->indexmax();
00136     if (offset==0)
00137     {
00138       int i,j;
00139       switch(old_style_flag)
00140       {
00141       case 0:
00142         for (i=imin;i<=imax;i++)
00143           for (j=imin;j<=imax;j++)
00144             H(i,j)+=(*pMinv)(i,j);
00145         break;
00146       case 1:
00147         for (i=imin;i<=imax;i++)
00148           for (j=imin;j<=imax;j++)
00149             H(i,j)+=2.0*(*pMinv)(i,j);
00150         break;
00151       case 2:
00152         for (i=imin;i<=imax;i++)
00153           for (j=imin;j<=imax;j++)
00154             H(i,j)+=2.0*(*pMinv)(i,j);
00155          break;
00156       default:
00157          cerr << "Illegal value in switch statement" << endl;
00158          ad_exit(1);
00159       }
00160     }
00161     else
00162     {
00163       int i,j;
00164       switch(old_style_flag)
00165       {
00166       case 0:
00167         for (i=imin;i<=imax;i++)
00168           for (j=imin;j<=imax;j++)
00169             H(offset+i,offset+j)+=(*pMinv)(i,j);
00170         break;
00171       case 1:
00172         for (i=imin;i<=imax;i++)
00173           for (j=imin;j<=imax;j++)
00174             H(offset+i,offset+j)+=2.0*(*pMinv)(i,j);
00175         break;
00176       case 2:
00177         for (i=imin;i<=imax;i++)
00178           for (j=imin;j<=imax;j++)
00179             H(offset+i,offset+j)+=2.0*(*pMinv)(i,j);
00180          break;
00181       default:
00182          cerr << "Illegal value in switch statement" << endl;
00183          ad_exit(1);
00184       }
00185     }
00186   }
00187 
00192   void quadratic_prior::get_cHessian(dvar_matrix H,int xsize)
00193   {
00194     int offset=get_offset(xsize);
00195     int imin=pMinv->indexmin();
00196     int imax=pMinv->indexmax();
00197     if (offset==0)
00198     {
00199       int i,j;
00200       switch(old_style_flag)
00201       {
00202       case 0:
00203         for (i=imin;i<=imax;i++)
00204           for (j=imin;j<=imax;j++)
00205             H(i,j)+=(*pMinv)(i,j);
00206         break;
00207       case 1:
00208         for (i=imin;i<=imax;i++)
00209           for (j=imin;j<=imax;j++)
00210             H(i,j)+=2.0*(*pMinv)(i,j);
00211         break;
00212       case 2:
00213         for (i=imin;i<=imax;i++)
00214           for (j=imin;j<=imax;j++)
00215             H(i,j)+=2.0*(*pMinv)(i,j);
00216          break;
00217       default:
00218          cerr << "Illegal value in switch statement" << endl;
00219          ad_exit(1);
00220       }
00221     }
00222     else
00223     {
00224       int i,j;
00225       switch(old_style_flag)
00226       {
00227       case 0:
00228         for (i=imin;i<=imax;i++)
00229           for (j=imin;j<=imax;j++)
00230             H(offset+i,offset+j)+=(*pMinv)(i,j);
00231         break;
00232       case 1:
00233         for (i=imin;i<=imax;i++)
00234           for (j=imin;j<=imax;j++)
00235             H(offset+i,offset+j)+=2.0*(*pMinv)(i,j);
00236         break;
00237       case 2:
00238         for (i=imin;i<=imax;i++)
00239           for (j=imin;j<=imax;j++)
00240             H(offset+i,offset+j)+=2.0*(*pMinv)(i,j);
00241          break;
00242       default:
00243          cerr << "Illegal value in switch statement" << endl;
00244          ad_exit(1);
00245       }
00246     }
00247   }
00248 
00253 void quadratic_prior::get_vHessian(dvar_matrix H,int xsize)
00254 {
00255   if (!dfpMinv)
00256   {
00257       cerr << "This can't happen" << endl;
00258       ad_exit(1);
00259   }
00260   else
00261   {
00262     int imin=dfpMinv->indexmin();
00263     int imax=dfpMinv->indexmax();
00264     int offset=get_offset(xsize);
00265     if (offset==0)
00266     {
00267       switch(old_style_flag)
00268       {
00269       case 0:
00270         for (int i=imin;i<=imax;i++)
00271           for (int j=imin;j<=imax;j++)
00272             H(i,j)+=(*dfpMinv)(i,j);
00273         break;
00274       case 1:
00275         for (int i=imin;i<=imax;i++)
00276           for (int j=imin;j<=imax;j++)
00277             H(i,j)+=2.0*(*dfpMinv)(i,j);
00278         break;
00279       case 2:
00280         for (int i=imin;i<=imax;i++)
00281           for (int j=imin;j<=imax;j++)
00282             H(i,j)+=2.0*(*dfpMinv)(i,j);
00283          break;
00284       default:
00285          cerr << "Illegal valueinswitch statement" << endl;
00286          ad_exit(1);
00287       }
00288     }
00289     else
00290     {
00291       switch(old_style_flag)
00292       {
00293       case 0:
00294         for (int i=imin;i<=imax;i++)
00295           for (int j=imin;j<=imax;j++)
00296             H(offset+i,offset+j)+=(*dfpMinv)(i,j);
00297         break;
00298       case 1:
00299         for (int i=imin;i<=imax;i++)
00300           for (int j=imin;j<=imax;j++)
00301             H(offset+i,offset+j)+=2.0*(*dfpMinv)(i,j);
00302         break;
00303       case 2:
00304         for (int i=imin;i<=imax;i++)
00305           for (int j=imin;j<=imax;j++)
00306             H(offset+i,offset+j)+=2.0*(*dfpMinv)(i,j);
00307          break;
00308       default:
00309          cerr << "Illegal valueinswitch statement" << endl;
00310          ad_exit(1);
00311       }
00312     }
00313   }
00314 }
00315 
00316  /*
00317   dvar_matrix quadratic_prior::get_vHessian(void)
00318   {
00319     return *dfpMinv;
00320     //return value(*pMinv);
00321   }
00322  */
00323 
00328    void quadratic_prior::get_cHessian_from_vHessian(dmatrix H,int xs)
00329   {
00330     int offset=get_offset(xs);
00331     int imin=dfpMinv->indexmin();
00332     int imax=dfpMinv->indexmax();
00333     if (offset==0)
00334     {
00335       int i,j;
00336       switch(old_style_flag)
00337       {
00338       case 0:
00339         for (i=imin;i<=imax;i++)
00340           for (j=imin;j<=imax;j++)
00341             H(i,j)+=value((*dfpMinv)(i,j));
00342         break;
00343       case 1:
00344         for (i=imin;i<=imax;i++)
00345           for (j=imin;j<=imax;j++)
00346             H(i,j)+=2.0*value((*dfpMinv)(i,j));
00347         break;
00348       case 2:
00349         for (i=imin;i<=imax;i++)
00350           for (j=imin;j<=imax;j++)
00351             H(i,j)+=2.0*value((*dfpMinv)(i,j));
00352         break;
00353          break;
00354       default:
00355          cerr << "Illegal valueinswitch statement" << endl;
00356          ad_exit(1);
00357       }
00358     }
00359     else
00360     {
00361       int i,j;
00362       switch(old_style_flag)
00363       {
00364       case 0:
00365         for (i=imin;i<=imax;i++)
00366           for (j=imin;j<=imax;j++)
00367             H(offset+i,offset+j)+=value((*dfpMinv)(i,j));
00368         break;
00369       case 1:
00370         for (i=imin;i<=imax;i++)
00371           for (j=imin;j<=imax;j++)
00372             H(offset+i,offset+j)+=2.0*value((*dfpMinv)(i,j));
00373         break;
00374       case 2:
00375         for (i=imin;i<=imax;i++)
00376           for (j=imin;j<=imax;j++)
00377             H(offset+i,offset+j)+=2.0*value((*dfpMinv)(i,j));
00378          break;
00379       default:
00380          cerr << "Illegal valueinswitch statement" << endl;
00381          ad_exit(1);
00382       }
00383     }
00384     //return value(*dfpMinv);
00385   }
00386 
00391   dvar_vector quadratic_prior::get_gradient(void)
00392   {
00393     return ((*pMinv)*(*pu));
00394   }
00395 
00400   void quadratic_prior::get_cgradient(dvector g,int xs)
00401   {
00402     int offset=get_offset(xs);
00403     dvector tg=((*pMinv)*value(*pu));
00404     int imin=pMinv->indexmin();
00405     int imax=pMinv->indexmax();
00406     if (offset==0)
00407     {
00408       int i;
00409       switch(old_style_flag)
00410       {
00411       case 0:
00412         for (i=imin;i<=imax;i++)
00413           g(i)+=tg(i);
00414         break;
00415       case 1:
00416         for (i=imin;i<=imax;i++)
00417           g(i)+=2.0*tg(i);
00418         break;
00419       case 2:
00420         for (i=imin;i<=imax;i++)
00421           g(i)+=2.0*tg(i);
00422         break;
00423       default:
00424         cerr << "Illegal valueinswitch statement" << endl;
00425         ad_exit(1);
00426       }
00427     }
00428     else
00429     {
00430       int i;
00431       switch(old_style_flag)
00432       {
00433       case 0:
00434         for (i=imin;i<=imax;i++)
00435           g(offset+i)+=tg(i);
00436         break;
00437       case 1:
00438         for (i=imin;i<=imax;i++)
00439           g(offset+i)+=2.0*tg(i);
00440         break;
00441       case 2:
00442         for (i=imin;i<=imax;i++)
00443           g(offset+i)+=2.0*tg(i);
00444         break;
00445       default:
00446         cerr << "Illegal valueinswitch statement" << endl;
00447         ad_exit(1);
00448       }
00449     }
00450     //return ((*pMinv)*value(*pu));
00451   }
00452 
00457   quadratic_prior::quadratic_prior(void)
00458   {
00459     pMinv=0;
00460     dfpMinv=0;
00461     pu=0;
00462     add_to_list();
00463   }
00464 
00469   quadratic_prior::~quadratic_prior(void)
00470   {
00471     if (pMinv) delete pMinv;
00472     pMinv=0;
00473     if (pu) delete pu;
00474     pu=0;
00475     if (dfpMinv) delete pMinv;
00476     dfpMinv=0;
00477   }
00478 
00483   void quadratic_prior::allocate( const dvar_vector & _u,const char * s)
00484   {
00485     allocate(_u);
00486   }
00487 
00492   void quadratic_prior::allocate(const dvar_vector & _u)
00493   {
00494     if (!allocated(_u))
00495     {
00496       cerr << "You must put random effects vector before"
00497        " quadtratic prior in the TPL file" << endl;
00498       ad_exit(1);
00499     }
00500     pu = new dvar_vector((dvar_vector&)(_u));
00501   }
00502 
00507   void quadratic_prior::allocate(const dvar_matrix & _M,
00508     const dvar_vector & _u,const char * s)
00509   {
00510     allocate(_M,_u);
00511   }
00512 
00517   void quadratic_prior::allocate(const dvar_matrix & _M,
00518     const dvar_vector & _u)
00519   {
00520     pMinv =new dmatrix(value(inv(_M)));
00521     pu = new dvar_vector((dvar_vector&)(_u));
00522   }
00523 
00528  dvariable quadratic_prior::get_quadratic_priors(void)
00529  {
00530    dvariable f=0.0;
00531    for (int i=0;i<num_quadratic_prior;i++)
00532    {
00533      f+=ptr[i]->get_function();
00534    }
00535    return f;
00536  }
00537 
00542  void quadratic_prior::get_cgradient_contribution(dvector g,int xs)
00543  {
00544    for (int i=0;i<num_quadratic_prior;i++)
00545    {
00546      ptr[i]->get_cgradient(g,xs);
00547     /*
00548      if (old_style_flag)
00549      {
00550        return ptr[i]->get_cgradient();
00551      }
00552      else
00553      {
00554        return ptr[i]->get_cgradient();
00555      }
00556     */
00557    }
00558    //return f;
00559  }
00560 /*
00561  dvar_vector quadratic_prior::get_gradient_contribution(void)
00562  {
00563    for (int i=0;i<num_quadratic_prior;i++)
00564    {
00565      return ptr[i]->get_gradient();
00566    }
00567    //return f;
00568  }
00569 */
00570 
00575  void quadratic_prior::get_cHessian_contribution(dmatrix H,int xsize)
00576  {
00577    for (int i=0;i<num_quadratic_prior;i++)
00578    {
00579      if (!ptr[i])
00580      {
00581        cerr << "ptr["<<i<<"] = 0 in"
00582          " quadratic_prior::get_cHessian_contribution" << endl;
00583        ad_exit(1);
00584      }
00585      if (!ptr[i]->pMinv)
00586      {
00587        cerr << "ptr["<<i<<"]->pMinv = 0 in"
00588          " quadratic_prior::get_cHessian_contribution" << endl;
00589        ad_exit(1);
00590      }
00591      if (!allocated(*(ptr[i]->pMinv)))
00592      {
00593        cerr << "*ptr["<<i<<"] is unallocated in"
00594          " quadratic_prior::get_cHessian_contribution" << endl;
00595        ad_exit(1);
00596      }
00597      ptr[i]->get_cHessian(H,xsize);
00598     /*
00599      if (old_style_flag)
00600      {
00601        return 2.0*ptr[i]->get_cHessian();
00602      }
00603      else
00604      {
00605        return ptr[i]->get_cHessian();
00606      }
00607     */
00608    }
00609    //return f;
00610  }
00611 
00616  void quadratic_prior::get_vHessian_contribution(dvar_matrix H,int xs)
00617  {
00618    for (int i=0;i<num_quadratic_prior;i++)
00619    {
00620      ptr[i]->get_vHessian(H,xs);
00621     /*
00622      if (old_style_flag)
00623      {
00624        return 2.0*ptr[i]->get_vHessian();
00625      }
00626      else
00627      {
00628        return ptr[i]->get_vHessian();
00629      }
00630     */
00631    }
00632    //return f;
00633  }
00634  /*
00635  dvar_matrix quadratic_prior::get_Hessian_contribution(void)
00636  {
00637    for (int i=0;i<num_quadratic_prior;i++)
00638    {
00639      return ptr[i]->get_Hessian();
00640    }
00641    //return f;
00642  }
00643  */
00644 
00649 void quadratic_prior::get_cHessian_contribution_from_vHessian(dmatrix Hess,
00650   int xsize)
00651  {
00652    for (int i=0;i<num_quadratic_prior;i++)
00653    {
00654      unsigned int nv=df1b2quadratic_prior::get_ptr(i)->
00655        get_num_active_parameters();
00656      if (nv)
00657        ptr[i]->get_cHessian_from_vHessian(Hess,xsize);
00658      else
00659        ptr[i]->get_cHessian(Hess,xsize);
00660    }
00661    //return f;
00662  }
00663 
00668  void quadratic_prior::operator = (const dvar_matrix & _M)
00669  {
00670    dvariable lndet;
00671    dvariable sgn;
00672 
00673    switch (quadratic_prior::old_style_flag)
00674    {
00675    case 0:
00676      *objective_function_value::pobjfun+=0.5*(*pu)*(solve(_M,*pu,lndet,sgn));
00677      *objective_function_value::pobjfun+=0.5*lndet;
00678      //*objective_function_value::pobjfun+=0.5*(*pu)*(solve(_M,*pu));
00679      //*objective_function_value::pobjfun+=0.5*ln_det(_M);
00680      break;
00681    case 1:
00682      *objective_function_value::pobjfun+=(*pu)*(solve(_M,*pu));
00683      break;
00684    case 2:
00685      *objective_function_value::pobjfun+=(*pu) * ( _M * (*pu) );
00686      break;
00687    default:
00688      cerr << "Illegal value for quadratic_prior::old_style_flag"
00689           << endl;
00690      ad_exit(1);
00691    }
00692    if (pMinv)
00693    {
00694      delete pMinv;
00695      pMinv=0;
00696    }
00697    if (dfpMinv)
00698    {
00699      delete dfpMinv;
00700      dfpMinv=0;
00701    }
00702    switch (quadratic_prior::old_style_flag)
00703    {
00704    case 0:
00705    case 1:
00706      if (laplace_approximation_calculator::where_are_we_flag==2)
00707      {
00708        pMinv = new dmatrix(inv(value(_M)));
00709        if (pMinv==0)
00710        {
00711          cerr << "Error allocating dmatrix" << endl;
00712          ad_exit(1);
00713        }
00714      }
00715      if (laplace_approximation_calculator::where_are_we_flag==3)
00716      {
00717        dfpMinv = new dvar_matrix(inv(_M));
00718        if (dfpMinv==0)
00719        {
00720          cerr << "Error allocating dvar_matrix" << endl;
00721          ad_exit(1);
00722        }
00723      }
00724      break;
00725    case 2:
00726      if (laplace_approximation_calculator::where_are_we_flag==2)
00727      {
00728        pMinv = new dmatrix(value(_M));
00729        if (pMinv==0)
00730        {
00731          cerr << "Error allocating dmatrix" << endl;
00732          ad_exit(1);
00733        }
00734      }
00735      if (laplace_approximation_calculator::where_are_we_flag==3)
00736      {
00737        unsigned int nv =
00738          df1b2quadratic_prior::get_ptr(xmyindex)->get_num_active_parameters();
00739        //if (nv==0)
00740        if (nv!=0)
00741        {
00742          dfpMinv = new dvar_matrix(_M);
00743          if (dfpMinv==0)
00744          {
00745            cerr << "Error allocating dvar_matrix" << endl;
00746            ad_exit(1);
00747          }
00748        }
00749        else
00750        {
00751          pMinv = new dmatrix(value(_M));
00752          if (pMinv==0)
00753          {
00754            cerr << "Error allocating dmatrix" << endl;
00755            ad_exit(1);
00756          }
00757        }
00758      }
00759      break;
00760    default:
00761      cerr << "Illegal value for quadratic_prior::old_style_flag"
00762           << endl;
00763      ad_exit(1);
00764    }
00765  }
00766 
00771  void quadratic_prior::operator = (const dmatrix & _M)
00772  {
00773    dvariable lndet;
00774    dvariable sgn;
00775 
00776    switch (quadratic_prior::old_style_flag)
00777    {
00778    case 0:
00779      cerr << " can't get here " << endl;
00780      ad_exit(1);
00781      break;
00782    case 1:
00783      cerr << " can't get here " << endl;
00784      ad_exit(1);
00785      break;
00786    case 2:
00787      *objective_function_value::pobjfun+=(*pu) * ( _M * (*pu) );
00788      break;
00789    default:
00790      cerr << "Illegal value for quadratic_prior::old_style_flag"
00791           << endl;
00792      ad_exit(1);
00793    }
00794    if (pMinv)
00795    {
00796      delete pMinv;
00797      pMinv=0;
00798    }
00799    if (dfpMinv)
00800    {
00801      delete dfpMinv;
00802      dfpMinv=0;
00803    }
00804    switch (quadratic_prior::old_style_flag)
00805    {
00806    case 0:
00807    case 1:
00808      cerr << " can't get here " << endl;
00809      ad_exit(1);
00810      break;
00811    case 2:
00812      if (laplace_approximation_calculator::where_are_we_flag==2 ||
00813        laplace_approximation_calculator::where_are_we_flag==3)
00814      {
00815        pMinv = new dmatrix(_M);
00816        if (pMinv==0)
00817        {
00818          cerr << "Error allocating dmatrix" << endl;
00819          ad_exit(1);
00820        }
00821      }
00822      break;
00823    default:
00824      cerr << "Illegal value for quadratic_prior::old_style_flag"
00825           << endl;
00826      ad_exit(1);
00827    }
00828  }
00829 
00834 normal_quadratic_prior::normal_quadratic_prior(void)
00835 {
00836   set_old_style_flag();
00837 }
00838 
00843 void normal_quadratic_prior::set_old_style_flag(void)
00844 {
00845   old_style_flag=0;
00846 }
00847 
00852 void normal_quadratic_prior::operator = (const dvar_matrix & M)
00853 {
00854   quadratic_prior::operator = (M);
00855 }
00856 
00861 quadratic_re_penalty::quadratic_re_penalty(void)
00862 {
00863   set_old_style_flag();
00864 }
00865 
00870 void quadratic_re_penalty::set_old_style_flag(void)
00871 {
00872   old_style_flag=2;
00873 }
00874 
00879 void quadratic_re_penalty::operator = (const dvar_matrix & M)
00880 {
00881   quadratic_prior::operator = (M);
00882 }
00883 
00888 void quadratic_re_penalty::operator = (const dmatrix & M)
00889 {
00890   quadratic_prior::operator = (M);
00891 }
00892 
00897 constant_quadratic_re_penalty::constant_quadratic_re_penalty(void)
00898 {
00899   set_old_style_flag();
00900 }
00901 
00906 void constant_quadratic_re_penalty::set_old_style_flag(void)
00907 {
00908   old_style_flag=2;
00909 }
00910 
00915 void constant_quadratic_re_penalty::operator = (const dmatrix & M)
00916 {
00917   quadratic_prior::operator = (M);
00918 }