ADMB Documentation  11.1.2497
 All Classes Files Functions Variables Typedefs Friends Defines
quadpri.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: quadpri.cpp 1919 2014-04-22 22:02:01Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #include <df1b2fnl.h>
00012 //#include <df1b2fun.h>
00013 
00014 int quadratic_prior::in_qp_calculations=0;
00015 
00016 // this should be a resizeable array
00017 quadratic_prior * quadratic_prior::ptr[100];
00018 
00019 int quadratic_prior::num_quadratic_prior=0;
00020 const int quadratic_prior::max_num_quadratic_prior=100;
00021 
00026  void quadratic_prior::get_M_calculations(void)
00027  {
00028    for (int i=0;i<num_quadratic_prior;i++)
00029    {
00030      //if (ptr[i]->get_num_active_parameters()>0)
00031      {
00032        ptr[i]->get_cM();
00033      }
00034    }
00035  }
00036 
00041 dvector evaluate_function_with_quadprior(const dvector& x,int usize,
00042   function_minimizer * pfmin)
00043 {
00044   int xsize=initial_params::nvarcalc();
00045   dvector g(1,xsize);
00046   gradcalc(0,g);
00047   //double f=0.0;
00048   independent_variables u(1,xsize);
00049   u=x;
00050   dvariable vf=0.0;
00051   initial_params::reset(dvar_vector(u));
00052   //vf=0.0;
00053   dvar_matrix Hess_all(1,usize,1,usize);
00054   *objective_function_value::pobjfun=0.0;
00055   // so that dvar_matrix Hessian contributions are calculated
00056   laplace_approximation_calculator::where_are_we_flag=3;
00057   pfmin->AD_uf_inner();
00058   if ( quadratic_prior::get_num_quadratic_prior()>0)
00059   {
00060     quadratic_prior::get_M_calculations();
00061   }
00062   laplace_approximation_calculator::where_are_we_flag=0;
00063 
00064   *objective_function_value::pobjfun=0.0;
00065   Hess_all=pfmin->lapprox->Hess;
00066   for (int i=0;i<quadratic_prior::get_num_quadratic_prior();i++)
00067   {
00068     //Hess_all += quadratic_prior::get_ptr(i)->get_vHessian();
00069     int nv=df1b2quadratic_prior::get_ptr(i)->get_num_active_parameters();
00070     if (nv>0)
00071       quadratic_prior::get_ptr(i)->get_vHessian(Hess_all,xsize);
00072     else
00073       quadratic_prior::get_ptr(i)->get_cHessian(Hess_all,xsize);
00074   }
00075   int sgn;
00076   vf=0.5*ln_det(Hess_all,sgn);
00077   gradcalc(xsize,g);
00078   return g;
00079 }
00080 
00085   void quadratic_prior::add_to_list(void)
00086   {
00087     if (num_quadratic_prior>=max_num_quadratic_prior) ad_exit(1);
00088     xmyindex=num_quadratic_prior;
00089     ptr[num_quadratic_prior++]=this;
00090   }
00091 
00096   dvariable quadratic_prior::get_function(void)
00097   {
00098     return (*pu)*((*pMinv)*(*pu));
00099   }
00100 
00105   dvar_matrix quadratic_prior::get_Hessian(void)
00106   {
00107     return *pMinv;
00108   }
00109 
00114   int quadratic_prior::get_offset(int xs)
00115   {
00116     df1b2_init_vector * fpu=df1b2quadratic_prior::ptr[get_myindex()]->pu;
00117     int mmin=(*fpu)(fpu->indexmin()).get_ind_index();
00118     return mmin-xs-1;
00119   }
00120   //dmatrix quadratic_prior::get_cHessian(void)
00121 
00126   void quadratic_prior::get_cHessian(dmatrix H,int xsize)
00127   {
00128     int offset=get_offset(xsize);
00129     int imin=pMinv->indexmin();
00130     int imax=pMinv->indexmax();
00131     if (offset==0)
00132     {
00133       int i,j;
00134       switch(old_style_flag)
00135       {
00136       case 0:
00137         for (i=imin;i<=imax;i++)
00138           for (j=imin;j<=imax;j++)
00139             H(i,j)+=(*pMinv)(i,j);
00140         break;
00141       case 1:
00142         for (i=imin;i<=imax;i++)
00143           for (j=imin;j<=imax;j++)
00144             H(i,j)+=2.0*(*pMinv)(i,j);
00145         break;
00146       case 2:
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       default:
00152          cerr << "Illegal value in switch statement" << endl;
00153          ad_exit(1);
00154       }
00155     }
00156     else
00157     {
00158       int i,j;
00159       switch(old_style_flag)
00160       {
00161       case 0:
00162         for (i=imin;i<=imax;i++)
00163           for (j=imin;j<=imax;j++)
00164             H(offset+i,offset+j)+=(*pMinv)(i,j);
00165         break;
00166       case 1:
00167         for (i=imin;i<=imax;i++)
00168           for (j=imin;j<=imax;j++)
00169             H(offset+i,offset+j)+=2.0*(*pMinv)(i,j);
00170         break;
00171       case 2:
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       default:
00177          cerr << "Illegal value in switch statement" << endl;
00178          ad_exit(1);
00179       }
00180     }
00181   }
00182 
00187   void quadratic_prior::get_cHessian(dvar_matrix H,int xsize)
00188   {
00189     int offset=get_offset(xsize);
00190     int imin=pMinv->indexmin();
00191     int imax=pMinv->indexmax();
00192     if (offset==0)
00193     {
00194       int i,j;
00195       switch(old_style_flag)
00196       {
00197       case 0:
00198         for (i=imin;i<=imax;i++)
00199           for (j=imin;j<=imax;j++)
00200             H(i,j)+=(*pMinv)(i,j);
00201         break;
00202       case 1:
00203         for (i=imin;i<=imax;i++)
00204           for (j=imin;j<=imax;j++)
00205             H(i,j)+=2.0*(*pMinv)(i,j);
00206         break;
00207       case 2:
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       default:
00213          cerr << "Illegal value in switch statement" << endl;
00214          ad_exit(1);
00215       }
00216     }
00217     else
00218     {
00219       int i,j;
00220       switch(old_style_flag)
00221       {
00222       case 0:
00223         for (i=imin;i<=imax;i++)
00224           for (j=imin;j<=imax;j++)
00225             H(offset+i,offset+j)+=(*pMinv)(i,j);
00226         break;
00227       case 1:
00228         for (i=imin;i<=imax;i++)
00229           for (j=imin;j<=imax;j++)
00230             H(offset+i,offset+j)+=2.0*(*pMinv)(i,j);
00231         break;
00232       case 2:
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       default:
00238          cerr << "Illegal value in switch statement" << endl;
00239          ad_exit(1);
00240       }
00241     }
00242   }
00243 
00248   void quadratic_prior::get_vHessian(dvar_matrix H,int xsize)
00249   {
00250     int offset=get_offset(xsize);
00251     if (dfpMinv==0)
00252     {
00253       cerr << "This can't happen" << endl;
00254       ad_exit(1);
00255     }
00256     int imin=dfpMinv->indexmin();
00257     int imax=dfpMinv->indexmax();
00258     if (offset==0)
00259     {
00260       int i,j;
00261       switch(old_style_flag)
00262       {
00263       case 0:
00264         for (i=imin;i<=imax;i++)
00265           for (j=imin;j<=imax;j++)
00266             H(i,j)+=(*dfpMinv)(i,j);
00267         break;
00268       case 1:
00269         for (i=imin;i<=imax;i++)
00270           for (j=imin;j<=imax;j++)
00271             H(i,j)+=2.0*(*dfpMinv)(i,j);
00272         break;
00273       case 2:
00274         for (i=imin;i<=imax;i++)
00275           for (j=imin;j<=imax;j++)
00276             H(i,j)+=2.0*(*dfpMinv)(i,j);
00277          break;
00278       default:
00279          cerr << "Illegal valueinswitch statement" << endl;
00280          ad_exit(1);
00281       }
00282     }
00283     else
00284     {
00285       int i,j;
00286       switch(old_style_flag)
00287       {
00288       case 0:
00289         for (i=imin;i<=imax;i++)
00290           for (j=imin;j<=imax;j++)
00291             H(offset+i,offset+j)+=(*dfpMinv)(i,j);
00292         break;
00293       case 1:
00294         for (i=imin;i<=imax;i++)
00295           for (j=imin;j<=imax;j++)
00296             H(offset+i,offset+j)+=2.0*(*dfpMinv)(i,j);
00297         break;
00298       case 2:
00299         for (i=imin;i<=imax;i++)
00300           for (j=imin;j<=imax;j++)
00301             H(offset+i,offset+j)+=2.0*(*dfpMinv)(i,j);
00302          break;
00303       default:
00304          cerr << "Illegal valueinswitch statement" << endl;
00305          ad_exit(1);
00306       }
00307     }
00308   }
00309  /*
00310   dvar_matrix quadratic_prior::get_vHessian(void)
00311   {
00312     return *dfpMinv;
00313     //return value(*pMinv);
00314   }
00315  */
00316 
00321    void quadratic_prior::get_cHessian_from_vHessian(dmatrix H,int xs)
00322   {
00323     int offset=get_offset(xs);
00324     int imin=dfpMinv->indexmin();
00325     int imax=dfpMinv->indexmax();
00326     if (offset==0)
00327     {
00328       int i,j;
00329       switch(old_style_flag)
00330       {
00331       case 0:
00332         for (i=imin;i<=imax;i++)
00333           for (j=imin;j<=imax;j++)
00334             H(i,j)+=value((*dfpMinv)(i,j));
00335         break;
00336       case 1:
00337         for (i=imin;i<=imax;i++)
00338           for (j=imin;j<=imax;j++)
00339             H(i,j)+=2.0*value((*dfpMinv)(i,j));
00340         break;
00341       case 2:
00342         for (i=imin;i<=imax;i++)
00343           for (j=imin;j<=imax;j++)
00344             H(i,j)+=2.0*value((*dfpMinv)(i,j));
00345         break;
00346          break;
00347       default:
00348          cerr << "Illegal valueinswitch statement" << endl;
00349          ad_exit(1);
00350       }
00351     }
00352     else
00353     {
00354       int i,j;
00355       switch(old_style_flag)
00356       {
00357       case 0:
00358         for (i=imin;i<=imax;i++)
00359           for (j=imin;j<=imax;j++)
00360             H(offset+i,offset+j)+=value((*dfpMinv)(i,j));
00361         break;
00362       case 1:
00363         for (i=imin;i<=imax;i++)
00364           for (j=imin;j<=imax;j++)
00365             H(offset+i,offset+j)+=2.0*value((*dfpMinv)(i,j));
00366         break;
00367       case 2:
00368         for (i=imin;i<=imax;i++)
00369           for (j=imin;j<=imax;j++)
00370             H(offset+i,offset+j)+=2.0*value((*dfpMinv)(i,j));
00371          break;
00372       default:
00373          cerr << "Illegal valueinswitch statement" << endl;
00374          ad_exit(1);
00375       }
00376     }
00377     //return value(*dfpMinv);
00378   }
00379 
00384   dvar_vector quadratic_prior::get_gradient(void)
00385   {
00386     return ((*pMinv)*(*pu));
00387   }
00388 
00393   void quadratic_prior::get_cgradient(dvector g,int xs)
00394   {
00395     int offset=get_offset(xs);
00396     dvector tg=((*pMinv)*value(*pu));
00397     int imin=pMinv->indexmin();
00398     int imax=pMinv->indexmax();
00399     if (offset==0)
00400     {
00401       int i;
00402       switch(old_style_flag)
00403       {
00404       case 0:
00405         for (i=imin;i<=imax;i++)
00406           g(i)+=tg(i);
00407         break;
00408       case 1:
00409         for (i=imin;i<=imax;i++)
00410           g(i)+=2.0*tg(i);
00411         break;
00412       case 2:
00413         for (i=imin;i<=imax;i++)
00414           g(i)+=2.0*tg(i);
00415         break;
00416       default:
00417         cerr << "Illegal valueinswitch statement" << endl;
00418         ad_exit(1);
00419       }
00420     }
00421     else
00422     {
00423       int i;
00424       switch(old_style_flag)
00425       {
00426       case 0:
00427         for (i=imin;i<=imax;i++)
00428           g(offset+i)+=tg(i);
00429         break;
00430       case 1:
00431         for (i=imin;i<=imax;i++)
00432           g(offset+i)+=2.0*tg(i);
00433         break;
00434       case 2:
00435         for (i=imin;i<=imax;i++)
00436           g(offset+i)+=2.0*tg(i);
00437         break;
00438       default:
00439         cerr << "Illegal valueinswitch statement" << endl;
00440         ad_exit(1);
00441       }
00442     }
00443     //return ((*pMinv)*value(*pu));
00444   }
00445 
00450   quadratic_prior::quadratic_prior(void)
00451   {
00452     pMinv=0;
00453     dfpMinv=0;
00454     pu=0;
00455     add_to_list();
00456   }
00457 
00462   quadratic_prior::~quadratic_prior(void)
00463   {
00464     if (pMinv) delete pMinv;
00465     pMinv=0;
00466     if (pu) delete pu;
00467     pu=0;
00468     if (dfpMinv) delete pMinv;
00469     dfpMinv=0;
00470   }
00471 
00476   void quadratic_prior::allocate( const dvar_vector & _u,const char * s)
00477   {
00478     allocate(_u);
00479   }
00480 
00485   void quadratic_prior::allocate(const dvar_vector & _u)
00486   {
00487     if (!allocated(_u))
00488     {
00489       cerr << "You must put random effects vector before"
00490        " quadtratic prior in the TPL file" << endl;
00491       ad_exit(1);
00492     }
00493     pu = new dvar_vector((dvar_vector&)(_u));
00494   }
00495 
00500   void quadratic_prior::allocate(const dvar_matrix & _M,
00501     const dvar_vector & _u,const char * s)
00502   {
00503     allocate(_M,_u);
00504   }
00505 
00510   void quadratic_prior::allocate(const dvar_matrix & _M,
00511     const dvar_vector & _u)
00512   {
00513     pMinv =new dmatrix(value(inv(_M)));
00514     pu = new dvar_vector((dvar_vector&)(_u));
00515   }
00516 
00521  dvariable quadratic_prior::get_quadratic_priors(void)
00522  {
00523    dvariable f=0.0;
00524    for (int i=0;i<num_quadratic_prior;i++)
00525    {
00526      f+=ptr[i]->get_function();
00527    }
00528    return f;
00529  }
00530 
00535  void quadratic_prior::get_cgradient_contribution(dvector g,int xs)
00536  {
00537    for (int i=0;i<num_quadratic_prior;i++)
00538    {
00539      ptr[i]->get_cgradient(g,xs);
00540     /*
00541      if (old_style_flag)
00542      {
00543        return ptr[i]->get_cgradient();
00544      }
00545      else
00546      {
00547        return ptr[i]->get_cgradient();
00548      }
00549     */
00550    }
00551    //return f;
00552  }
00553 /*
00554  dvar_vector quadratic_prior::get_gradient_contribution(void)
00555  {
00556    for (int i=0;i<num_quadratic_prior;i++)
00557    {
00558      return ptr[i]->get_gradient();
00559    }
00560    //return f;
00561  }
00562 */
00563 
00568  void quadratic_prior::get_cHessian_contribution(dmatrix H,int xsize)
00569  {
00570    for (int i=0;i<num_quadratic_prior;i++)
00571    {
00572      if (!ptr[i])
00573      {
00574        cerr << "ptr["<<i<<"] = 0 in"
00575          " quadratic_prior::get_cHessian_contribution" << endl;
00576        ad_exit(1);
00577      }
00578      if (!ptr[i]->pMinv)
00579      {
00580        cerr << "ptr["<<i<<"]->pMinv = 0 in"
00581          " quadratic_prior::get_cHessian_contribution" << endl;
00582        ad_exit(1);
00583      }
00584      if (!allocated(*(ptr[i]->pMinv)))
00585      {
00586        cerr << "*ptr["<<i<<"] is unallocated in"
00587          " quadratic_prior::get_cHessian_contribution" << endl;
00588        ad_exit(1);
00589      }
00590      ptr[i]->get_cHessian(H,xsize);
00591     /*
00592      if (old_style_flag)
00593      {
00594        return 2.0*ptr[i]->get_cHessian();
00595      }
00596      else
00597      {
00598        return ptr[i]->get_cHessian();
00599      }
00600     */
00601    }
00602    //return f;
00603  }
00604 
00609  void quadratic_prior::get_vHessian_contribution(dvar_matrix H,int xs)
00610  {
00611    for (int i=0;i<num_quadratic_prior;i++)
00612    {
00613      ptr[i]->get_vHessian(H,xs);
00614     /*
00615      if (old_style_flag)
00616      {
00617        return 2.0*ptr[i]->get_vHessian();
00618      }
00619      else
00620      {
00621        return ptr[i]->get_vHessian();
00622      }
00623     */
00624    }
00625    //return f;
00626  }
00627  /*
00628  dvar_matrix quadratic_prior::get_Hessian_contribution(void)
00629  {
00630    for (int i=0;i<num_quadratic_prior;i++)
00631    {
00632      return ptr[i]->get_Hessian();
00633    }
00634    //return f;
00635  }
00636  */
00637 
00642 void quadratic_prior::get_cHessian_contribution_from_vHessian(dmatrix Hess,
00643   int xsize)
00644  {
00645    for (int i=0;i<num_quadratic_prior;i++)
00646    {
00647      int nv=df1b2quadratic_prior::get_ptr(i)->
00648        get_num_active_parameters();
00649      if (nv)
00650        ptr[i]->get_cHessian_from_vHessian(Hess,xsize);
00651      else
00652        ptr[i]->get_cHessian(Hess,xsize);
00653    }
00654    //return f;
00655  }
00656 
00661  void quadratic_prior::operator = (const dvar_matrix & _M)
00662  {
00663    dvariable lndet;
00664    dvariable sgn;
00665 
00666    switch (quadratic_prior::old_style_flag)
00667    {
00668    case 0:
00669      *objective_function_value::pobjfun+=0.5*(*pu)*(solve(_M,*pu,lndet,sgn));
00670      *objective_function_value::pobjfun+=0.5*lndet;
00671      //*objective_function_value::pobjfun+=0.5*(*pu)*(solve(_M,*pu));
00672      //*objective_function_value::pobjfun+=0.5*ln_det(_M);
00673      break;
00674    case 1:
00675      *objective_function_value::pobjfun+=(*pu)*(solve(_M,*pu));
00676      break;
00677    case 2:
00678      *objective_function_value::pobjfun+=(*pu) * ( _M * (*pu) );
00679      break;
00680    default:
00681      cerr << "Illegal value for quadratic_prior::old_style_flag"
00682           << endl;
00683      ad_exit(1);
00684    }
00685    if (pMinv)
00686    {
00687      delete pMinv;
00688      pMinv=0;
00689    }
00690    if (dfpMinv)
00691    {
00692      delete dfpMinv;
00693      dfpMinv=0;
00694    }
00695    switch (quadratic_prior::old_style_flag)
00696    {
00697    case 0:
00698    case 1:
00699      if (laplace_approximation_calculator::where_are_we_flag==2)
00700      {
00701        pMinv = new dmatrix(inv(value(_M)));
00702        if (pMinv==0)
00703        {
00704          cerr << "Error allocating dmatrix" << endl;
00705          ad_exit(1);
00706        }
00707      }
00708      if (laplace_approximation_calculator::where_are_we_flag==3)
00709      {
00710        dfpMinv = new dvar_matrix(inv(_M));
00711        if (dfpMinv==0)
00712        {
00713          cerr << "Error allocating dvar_matrix" << endl;
00714          ad_exit(1);
00715        }
00716      }
00717      break;
00718    case 2:
00719      if (laplace_approximation_calculator::where_are_we_flag==2)
00720      {
00721        pMinv = new dmatrix(value(_M));
00722        if (pMinv==0)
00723        {
00724          cerr << "Error allocating dmatrix" << endl;
00725          ad_exit(1);
00726        }
00727      }
00728      if (laplace_approximation_calculator::where_are_we_flag==3)
00729      {
00730        int nv =
00731          df1b2quadratic_prior::get_ptr(xmyindex)->get_num_active_parameters();
00732        //if (nv==0)
00733        if (nv!=0)
00734        {
00735          dfpMinv = new dvar_matrix(_M);
00736          if (dfpMinv==0)
00737          {
00738            cerr << "Error allocating dvar_matrix" << endl;
00739            ad_exit(1);
00740          }
00741        }
00742        else
00743        {
00744          pMinv = new dmatrix(value(_M));
00745          if (pMinv==0)
00746          {
00747            cerr << "Error allocating dmatrix" << endl;
00748            ad_exit(1);
00749          }
00750        }
00751      }
00752      break;
00753    default:
00754      cerr << "Illegal value for quadratic_prior::old_style_flag"
00755           << endl;
00756      ad_exit(1);
00757    }
00758  }
00759 
00764  void quadratic_prior::operator = (const dmatrix & _M)
00765  {
00766    dvariable lndet;
00767    dvariable sgn;
00768 
00769    switch (quadratic_prior::old_style_flag)
00770    {
00771    case 0:
00772      cerr << " can't get here " << endl;
00773      ad_exit(1);
00774      break;
00775    case 1:
00776      cerr << " can't get here " << endl;
00777      ad_exit(1);
00778      break;
00779    case 2:
00780      *objective_function_value::pobjfun+=(*pu) * ( _M * (*pu) );
00781      break;
00782    default:
00783      cerr << "Illegal value for quadratic_prior::old_style_flag"
00784           << endl;
00785      ad_exit(1);
00786    }
00787    if (pMinv)
00788    {
00789      delete pMinv;
00790      pMinv=0;
00791    }
00792    if (dfpMinv)
00793    {
00794      delete dfpMinv;
00795      dfpMinv=0;
00796    }
00797    switch (quadratic_prior::old_style_flag)
00798    {
00799    case 0:
00800    case 1:
00801      cerr << " can't get here " << endl;
00802      ad_exit(1);
00803      break;
00804    case 2:
00805      if (laplace_approximation_calculator::where_are_we_flag==2 ||
00806        laplace_approximation_calculator::where_are_we_flag==3)
00807      {
00808        pMinv = new dmatrix(_M);
00809        if (pMinv==0)
00810        {
00811          cerr << "Error allocating dmatrix" << endl;
00812          ad_exit(1);
00813        }
00814      }
00815      break;
00816    default:
00817      cerr << "Illegal value for quadratic_prior::old_style_flag"
00818           << endl;
00819      ad_exit(1);
00820    }
00821  }
00822 
00827 normal_quadratic_prior::normal_quadratic_prior(void)
00828 {
00829   set_old_style_flag();
00830 }
00831 
00836 void normal_quadratic_prior::set_old_style_flag(void)
00837 {
00838   old_style_flag=0;
00839 }
00840 
00845 void normal_quadratic_prior::operator = (const dvar_matrix & M)
00846 {
00847   quadratic_prior::operator = (M);
00848 }
00849 
00854 quadratic_re_penalty::quadratic_re_penalty(void)
00855 {
00856   set_old_style_flag();
00857 }
00858 
00863 void quadratic_re_penalty::set_old_style_flag(void)
00864 {
00865   old_style_flag=2;
00866 }
00867 
00872 void quadratic_re_penalty::operator = (const dvar_matrix & M)
00873 {
00874   quadratic_prior::operator = (M);
00875 }
00876 
00881 void quadratic_re_penalty::operator = (const dmatrix & M)
00882 {
00883   quadratic_prior::operator = (M);
00884 }
00885 
00890 constant_quadratic_re_penalty::constant_quadratic_re_penalty(void)
00891 {
00892   set_old_style_flag();
00893 }
00894 
00899 void constant_quadratic_re_penalty::set_old_style_flag(void)
00900 {
00901   old_style_flag=2;
00902 }
00903 
00908 void constant_quadratic_re_penalty::operator = (const dmatrix & M)
00909 {
00910   quadratic_prior::operator = (M);
00911 }