ADMB Documentation  11.1.1927
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2lp8.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2lp8.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 #if defined(USE_LAPLACE)
00012 #  include <admodel.h>
00013 #  include <df1b2fun.h>
00014 #  include <adrndeff.h>
00015 double fcomp1(dvector x,dvector d,int samplesize,int n,dvector & g,
00016   dmatrix& M);
00017 
00022 void laplace_approximation_calculator::make_sparse_triplet(void)
00023 {
00024   //int i;
00025   /*
00026   int mmax=Hess.indexmax();
00027   int nz=sum(Hess);
00028   if (sparse_triplet)
00029   {
00030     delete sparse_triplet;
00031     sparse_triplet=0;
00032   }
00033   */
00034   int nz2=0;
00035   if (compressed_triplet_information)
00036   {
00037     imatrix & cti = *compressed_triplet_information;
00038     int mmin=cti(cti.indexmin()).indexmin();
00039     int mmax=cti(cti.indexmin()).indexmax();
00040     nz2=1;
00041     int lmin=cti(2,1);
00042     for (int i=mmin+1;i<=mmax;i++)
00043     {
00044       if (cti(1,i)>cti(1,i-1))
00045       {
00046         nz2++;
00047         lmin=cti(2,i);
00048       }
00049       else if (cti(2,i)>lmin)
00050       {
00051         lmin=cti(2,i);
00052         nz2++;
00053       }
00054     }
00055   }
00056   //cout << nz2-nz << endl;
00057 
00058   if (sparse_triplet)
00059   {
00060     delete sparse_triplet;
00061     sparse_triplet=0;
00062   }
00063   //sparse_triplet = new dcompressed_triplet(1,nz2,usize,usize);
00064   if (sparse_triplet2)
00065   {
00066     delete sparse_triplet2;
00067     sparse_triplet2=0;
00068   }
00069   sparse_triplet2 = new dcompressed_triplet(1,nz2,usize,usize);
00070 
00071   if (compressed_triplet_information)
00072   {
00073     imatrix & cti = *compressed_triplet_information;
00074     int mmin=cti(cti.indexmin()).indexmin();
00075     int mmax=cti(cti.indexmin()).indexmax();
00076     if (sparse_iterator)
00077     {
00078       delete sparse_iterator;
00079       sparse_iterator=0;
00080     }
00081     sparse_iterator=new ivector(mmin,mmax);
00082     int ii=1;
00083     int lmin=cti(2,1);
00084     (*sparse_triplet2)(1,ii)=cti(1,1);
00085     (*sparse_triplet2)(2,ii)=cti(2,1);
00086     (*sparse_iterator)(cti(3,1))=ii;
00087     for (int i=mmin+1;i<=mmax;i++)
00088     {
00089       if (cti(1,i)>cti(1,i-1))
00090       {
00091         ii++;
00092         (*sparse_triplet2)(1,ii)=cti(1,i);
00093         (*sparse_triplet2)(2,ii)=cti(2,i);
00094         lmin=cti(2,i);
00095       }
00096       else if (cti(2,i)>lmin)
00097       {
00098         ii++;
00099         (*sparse_triplet2)(1,ii)=cti(1,i);
00100         (*sparse_triplet2)(2,ii)=cti(2,i);
00101         lmin=cti(2,i);
00102       }
00103       (*sparse_iterator)(cti(3,i))=ii;
00104     }
00105   }
00106   //cout << setw(8) << setprecision(2) << setscientific() << rowsum(Hess)
00107   //     << endl;
00108   //cout << setw(8) << setprecision(2) << setscientific() << Hess << endl;
00109  /*
00110   int ii=0;
00111   for (i=1;i<=mmax;i++)
00112   {
00113     for (int j=i;j<=mmax;j++)
00114     {
00115       if (Hess(i,j) != 0.0)
00116       {
00117         ++ii;
00118         (*sparse_triplet)(1,ii)=i;
00119         (*sparse_triplet)(2,ii)=j;
00120         (*sparse_triplet)(ii)=0.0;
00121       }
00122     }
00123   }
00124   */
00125   //sparse_symbolic = new hs_symbolic(*sparse_triplet,1);
00126   sparse_symbolic2 = new hs_symbolic(*sparse_triplet2,1);
00127 }
00128 
00133 void laplace_approximation_calculator::generate_antithetical_rvs()
00134 {
00135   // number of random vectors
00136   const ivector & itmp=(*num_local_re_array)(1,num_separable_calls);
00137   //const ivector & itmpf=(*num_local_fixed_array)(1,num_separable_calls);
00138   int i;
00139   for (i=2;i<=num_separable_calls;i++)
00140   {
00141     if (itmp(i) != itmp(i-1))
00142     {
00143       cerr << "At present can only use antithetical rv's when "
00144               "all separable calls are the same size" << endl;
00145       ad_exit(1);
00146     }
00147   }
00148   int n=itmp(1);
00149   double samplesize=num_importance_samples;
00150 
00151   // mesh size
00152   double delta=0.01;
00153   // maximum of distribution is near here
00154   double mid=sqrt(double(n-1));
00155   dmatrix weights(1,2*n,1,2);
00156   double spread=15;
00157   if (mid-spread<=0.001)
00158     spread=mid-0.1;
00159   double ssum=0.0;
00160   double x=0.0;
00161   double tmax=(n-1)*log(mid)-0.5*mid*mid;
00162   for (x=mid-spread;x<=mid+spread;x+=delta)
00163   {
00164     ssum+=exp((n-1)*log(x)-0.5*x*x-tmax);
00165   }
00166   double tsum=0;
00167   dvector dist(1,samplesize+1);
00168   dist.initialize();
00169   int is=0;
00170   int ii;
00171   for (x=mid-spread;x<=mid+spread;x+=delta)
00172   {
00173     tsum+=exp((n-1)*log(x)-0.5*x*x-tmax)/ssum*samplesize;
00174     int ns=int(tsum);
00175     for (ii=1;ii<=ns;ii++)
00176     {
00177       dist(++is)=x;
00178     }
00179     tsum-=ns;
00180   }
00181   if (is==samplesize-1)
00182   {
00183     dist(samplesize)=mid;
00184   }
00185   else if (is<samplesize-1)
00186   {
00187     cerr << "This can't happen" << endl;
00188     exit(1);
00189   }
00190 
00191   // get random numbers
00192 
00193   random_number_generator rng(rseed);
00194   if (antiepsilon)
00195   {
00196     if (allocated(*antiepsilon))
00197     {
00198       delete antiepsilon;
00199       antiepsilon=0;
00200     }
00201   }
00202   antiepsilon=new dmatrix(1,samplesize,1,n);
00203   dmatrix & M=*antiepsilon;
00204   M.fill_randn(rng);
00205 
00206   for (i=1;i<=samplesize;i++)
00207   {
00208     M(i)=M(i)/norm(M(i));
00209   }
00210   int nvar=(samplesize-1)*n;
00211   independent_variables xx(1,nvar);
00212   ii=0;
00213   for (i=2;i<=samplesize;i++)
00214   {
00215     int j;
00216     for (j=1;j<=n;j++)
00217     {
00218       xx(++ii)=M(i,j);
00219     }
00220   }
00221 
00222   fmmt1 fmc(nvar,5);
00223   //fmm fmc(nvar,5);
00224   fmc.noprintx=1;
00225   fmc.iprint=10;
00226   fmc.maxfn=2500;
00227   fmc.crit=1.e-6;
00228 
00229   double f;
00230   double fbest=1.e+50;;
00231   dvector g(1,nvar);
00232   dvector gbest(1,nvar);
00233   dvector xbest(1,nvar);
00234 
00235   gbest.fill_seqadd(1.e+50,0.);
00236   {
00237     while (fmc.ireturn>=0)
00238     {
00239       //int badflag=0;
00240       fmc.fmin(f,xx,g);
00241       if (fmc.ihang)
00242       {
00243         //int hang_flag=fmc.ihang;
00244         //double maxg=max(g);
00245         //double ccrit=fmc.crit;
00246         //int current_ifn=fmc.ifn;
00247       }
00248       if (fmc.ireturn>0)
00249       {
00250          f=fcomp1(xx,dist,samplesize,n,g,M);
00251          if (f < fbest)
00252          {
00253            fbest=f;
00254            gbest=g;
00255            xbest=xx;
00256          }
00257        }
00258      }
00259      xx=xbest;
00260   }
00261   ii=0;
00262   for (i=2;i<=samplesize;i++)
00263   {
00264     int j;
00265     for (j=1;j<=n;j++)
00266     {
00267       M(i,j)=xx(++ii);
00268     }
00269   }
00270   for (i=1;i<=samplesize;i++)
00271   {
00272     M(i)*=dist(i)/norm(M(i));
00273   }
00274 }
00275 
00280 double fcomp1(dvector x,dvector d,int samplesize,int n,dvector & g,
00281   dmatrix& M)
00282 {
00283   dmatrix VM(1,samplesize,1,n);
00284   dmatrix C(1,samplesize,1,samplesize);
00285   dmatrix VM0(1,samplesize,1,n);
00286   dvector N(1,samplesize);
00287   dmatrix dfVM(1,samplesize,1,n);
00288   dmatrix dfVM0(1,samplesize,1,n);
00289   dvector dfN(1,samplesize);
00290   dfVM.initialize();
00291   dfVM0.initialize();
00292   dfN.initialize();
00293 
00294   double f=0.0;
00295   int ii=0;
00296   int i,j;
00297   VM0(1)=M(1);
00298   for (i=2;i<=samplesize;i++)
00299   {
00300     for (j=1;j<=n;j++)
00301     {
00302       VM0(i,j)=x(++ii);
00303     }
00304   }
00305   for (i=1;i<=samplesize;i++)
00306   {
00307     N(i)=norm(VM0(i));
00308     VM(i)=VM0(i)*(d(i)/N(i));
00309   }
00310   for (i=1;i<=samplesize;i++)
00311   {
00312     for (ii=i+1;ii<=samplesize;ii++)
00313     {
00314       //C(i,ii)=1.0/(0.01+norm2(VM(i)-VM(ii)));
00315       //f+=C(i,ii);
00316       C(i,ii)=norm2(VM(i)-VM(ii));
00317       f-=C(i,ii);
00318       //f+=1.0/(0.01+norm2(VM(i)-VM(ii)));
00319     }
00320     f+=100.0*square(log(N(i)));
00321   }
00322   for (i=1;i<=samplesize;i++)
00323   {
00324     //f+=100.0*square(log(N(i)));
00325     dfN(i)+=200*log(N(i))/N(i);
00326     for (ii=i+1;ii<=samplesize;ii++)
00327     {
00328       //f+=1.0/(0.01+norm2(VM(i)-VM(ii)));
00329       //double tmp=-1.0/square(0.01+norm2(VM(i)-VM(ii)));
00330       //double tmp=-square(C(i,ii));
00331       dvector vtmp=-2.0*(VM(i)-VM(ii));
00332       dfVM(i)+=vtmp;
00333       dfVM(ii)-=vtmp;
00334     }
00335   }
00336   for (i=1;i<=samplesize;i++)
00337   {
00338     //VM(i)=VM0(i)*(d(i)/N(i));
00339     dfVM0(i)=dfVM(i)*d(i)/N(i);
00340     dfN(i)-=(dfVM(i)*VM0(i))*d(i)/square(N(i));
00341 
00342     //N(i)=norm(VM0(i));
00343     dfVM0(i)+=dfN(i)/N(i)*VM0(i);
00344   }
00345   ii=0;
00346   for (i=2;i<=samplesize;i++)
00347   {
00348     for (j=1;j<=n;j++)
00349     {
00350       //VM0(i,j)=vx(++ii);
00351       g(++ii)=dfVM0(i,j);
00352     }
00353   }
00354   return f;
00355 }
00356 
00361 void laplace_approximation_calculator::check_hessian_type(const dvector& _x,
00362   function_minimizer * pfmin)
00363 {
00364   pfmin->pre_user_function();
00365 }
00366 
00371 void function_minimizer::pre_user_function(void)
00372 {
00373   if (lapprox)
00374   {
00375     if (lapprox->hesstype==2)
00376     {
00377       lapprox->separable_calls_counter=0;
00378     }
00379   }
00380   user_function();
00381  /*
00382   if (lapprox)
00383   {
00384     if (lapprox->hesstype==2)
00385     {
00386       lapprox->nested_shape.trim();
00387       cout << lapprox->nested_shape;
00388       lapprox->nested_indices.allocate(lapprox->nested_shape);
00389       lapprox->separable_calls_counter=0;
00390       user_function();
00391     }
00392   }
00393  */
00394 }
00395 
00400 void laplace_approximation_calculator::
00401   check_hessian_type(function_minimizer * pfmin)
00402 {
00403   int ip = 0;
00404   if (quadratic_prior::get_num_quadratic_prior()>0)
00405   {
00406     hesstype=4;
00407     if (allocated(Hess))
00408     {
00409       if (Hess.indexmax()!=usize)
00410       {
00411         Hess.deallocate();
00412         Hess.allocate(1,usize,1,usize);
00413       }
00414     }
00415     else
00416     {
00417        Hess.allocate(1,usize,1,usize);
00418     }
00419     if (allocated(Hessadjoint))
00420     {
00421       if (Hessadjoint.indexmax()!=usize)
00422       {
00423         Hessadjoint.deallocate();
00424         Hessadjoint.allocate(1,usize,1,usize);
00425       }
00426     }
00427     else
00428     {
00429        Hessadjoint.allocate(1,usize,1,usize);
00430     }
00431     return;
00432   }
00433   else
00434   {
00435     int nv=initial_df1b2params::set_index();
00436     if (allocated(used_flags))
00437     {
00438       if (used_flags.indexmax() != nv)
00439       {
00440         used_flags.safe_deallocate();
00441       }
00442     }
00443     if (!allocated(used_flags))
00444     {
00445       used_flags.safe_allocate(1,nv);
00446     }
00447 
00448     //for (ip=1;ip<=num_der_blocks;ip++)
00449     {
00450       used_flags.initialize();
00451       // do we need to reallocate memory for df1b2variables?
00452       check_for_need_to_reallocate(ip);
00453       df1b2_gradlist::set_no_derivatives();
00454       //cout << re_objective_function_value::pobjfun << endl;
00455       //cout << re_objective_function_value::pobjfun->ptr << endl;
00456       (*re_objective_function_value::pobjfun)=0;
00457       df1b2variable pen=0.0;
00458       df1b2variable zz=0.0;
00459 
00460       initial_df1b2params::reset(y,pen);
00461       // call function to do block diagonal newton-raphson
00462       // the step vector from the newton-raphson is in the vector step
00463       df1b2_gradlist::set_no_derivatives();
00464 
00465       funnel_init_var::lapprox=this;
00466       block_diagonal_flag=5;
00467 
00468       quadratic_prior::in_qp_calculations=1;
00469 
00470       if (sparse_hessian_flag)
00471       {
00472         // just to get the number of separable calls
00473         separable_calls_counter=0;
00474         pfmin->AD_uf_inner();
00475         // allocate space for uncompressed sparse hessian information
00476 
00477         //num_separable_calls=separable_calls_counter;
00478         if (triplet_information==0)
00479         {
00480           triplet_information =new i3_array(1,separable_calls_counter);
00481         }
00482         else if ( triplet_information->indexmax() != separable_calls_counter)
00483         {
00484           delete triplet_information;
00485           triplet_information =new i3_array(1,separable_calls_counter);
00486         }
00487         triplet_information->initialize();
00488         separable_calls_counter=0;
00489       }
00490 
00491       pfmin->pre_user_function();
00492 
00493 
00494       if (sparse_hessian_flag)
00495       {
00496         // turn triplet_informaiton into  compressed_triplet_information
00497         int mmin= triplet_information->indexmin();
00498         int mmax= triplet_information->indexmax();
00499         int ndim=0;
00500         for (int i=mmin;i<=mmax;i++)
00501         {
00502           if (allocated((*triplet_information)(i)))
00503           {
00504             ndim+=(*triplet_information)(i,1).indexmax();
00505           }
00506         }
00507         if (compressed_triplet_information)
00508         {
00509           delete compressed_triplet_information;
00510           compressed_triplet_information=0;
00511         }
00512         compressed_triplet_information=new imatrix(1,ndim,1,3);
00513         (*compressed_triplet_information)(3).fill_seqadd(1,1);
00514         int ii=0;
00515         for (int i=mmin;i<=mmax;i++)
00516         {
00517           if (allocated((*triplet_information)(i)))
00518           {
00519             int jmin=(*triplet_information)(i,1).indexmin();
00520             int jmax=(*triplet_information)(i,1).indexmax();
00521             for (int j=jmin;j<=jmax;j++)
00522             {
00523               ii++;
00524               (*compressed_triplet_information)(ii,1)=
00525                 (*triplet_information)(i,1,j);
00526               (*compressed_triplet_information)(ii,2)=
00527                 (*triplet_information)(i,2,j);
00528               (*compressed_triplet_information)(ii,3)=ii;
00529             }
00530           }
00531         }
00532         imatrix & cti= *compressed_triplet_information;
00533         cti=sort(cti,1);
00534         int lmin=1;
00535         int lmax=0;
00536         for (int i=2;i<=ndim;i++)
00537         {
00538           if (cti(i,1)>cti(i-1,1))
00539           {
00540             lmax=i-1;
00541             cti.sub(lmin,lmax)=sort(cti.sub(lmin,lmax),2);
00542             lmin=i;
00543           }
00544         }
00545         cti.sub(lmin,ndim)=sort(cti.sub(lmin,ndim),2);
00546         imatrix tmp=trans(cti);
00547         delete compressed_triplet_information;
00548         compressed_triplet_information=new imatrix(tmp);
00549       }
00550 
00551       quadratic_prior::in_qp_calculations=0;
00552 
00553       int non_block_diagonal=0;
00554       for (int i=xsize+1;i<=xsize+usize;i++)
00555       {
00556         if (used_flags(i)>1)
00557         {
00558           non_block_diagonal=1;
00559           break;
00560         }
00561       }
00562       if (non_block_diagonal)
00563       {
00564         if (bw< usize/2 && sparse_hessian_flag==0)
00565         {
00566           hesstype=3;  //banded
00567           if (bHess)
00568           {
00569             if (bHess->bandwidth() !=bw)
00570             {
00571               delete bHess;
00572               bHess = new banded_symmetric_dmatrix(1,usize,bw);
00573               if (bHess==0)
00574               {
00575                 cerr << "Error allocating banded_symmetric_dmatrix" << endl;
00576                 ad_exit(1);
00577               }
00578             }
00579           }
00580           else
00581           {
00582             bHess = new banded_symmetric_dmatrix(1,usize,bw);
00583             if (bHess==0)
00584             {
00585               cerr << "Error allocating banded_symmetric_dmatrix" << endl;
00586               ad_exit(1);
00587             }
00588           }
00589           if (bHessadjoint)
00590           {
00591             if (bHessadjoint->bandwidth() !=bw)
00592             {
00593               delete bHessadjoint;
00594               bHessadjoint = new banded_symmetric_dmatrix(1,usize,bw);
00595               if (bHessadjoint==0)
00596               {
00597                 cerr << "Error allocating banded_symmetric_dmatrix" << endl;
00598                 ad_exit(1);
00599               }
00600             }
00601           }
00602           else
00603           {
00604             bHessadjoint = new banded_symmetric_dmatrix(1,usize,bw);
00605             if (bHessadjoint==0)
00606             {
00607               cerr << "Error allocating banded_symmetric_dmatrix" << endl;
00608               ad_exit(1);
00609             }
00610           }
00611         }
00612         else
00613         {
00614           //check_sparse_matrix_structure();
00615           hesstype=4;  // band is so wide so use full matrix
00616           if (bHess)
00617           {
00618             delete bHess;
00619             bHess=0;
00620           }
00621 
00622           if (bHessadjoint)
00623           {
00624             delete bHessadjoint;
00625             bHessadjoint=0;
00626           }
00627 
00628           if (allocated(Hess))
00629           {
00630             if (sparse_hessian_flag)
00631             {
00632               Hess.deallocate();
00633             }
00634             else
00635             {
00636               if (Hess.indexmax() != usize)
00637               {
00638                 Hess.deallocate();
00639                 Hess.allocate(1,usize,1,usize);
00640               }
00641             }
00642           }
00643           else
00644           {
00645             if (sparse_hessian_flag==0)
00646               Hess.allocate(1,usize,1,usize);
00647           }
00648           if (sparse_hessian_flag)
00649           {
00650             make_sparse_triplet();
00651           }
00652 
00653           if (allocated(Hessadjoint))
00654           {
00655             if (sparse_hessian_flag)
00656             {
00657               Hess.deallocate();
00658             }
00659             else
00660             {
00661               if (Hessadjoint.indexmax() != usize)
00662               {
00663                 Hessadjoint.deallocate();
00664                 Hessadjoint.allocate(1,usize,1,usize);
00665               }
00666             }
00667           }
00668           else
00669           {
00670             if (sparse_hessian_flag==0)
00671               Hessadjoint.allocate(1,usize,1,usize);
00672           }
00673         }
00674       }
00675       else
00676       {
00677         hesstype=2;
00678       }
00679       if (hesstype==2 && num_importance_samples>0)
00680       {
00681         if (importance_sampling_components)
00682         {
00683           delete importance_sampling_components;
00684           importance_sampling_components=0;
00685         }
00686         importance_sampling_components=
00687           new dvar_matrix(1,pmin->lapprox->num_separable_calls,
00688             1,num_importance_samples);
00689       }
00690 
00691       if (hesstype==2 && (num_importance_samples>0 || use_gauss_hermite>0))
00692       {
00693         const ivector & itmp=(*num_local_re_array)(1,num_separable_calls);
00694         const ivector & itmpf=(*num_local_fixed_array)(1,num_separable_calls);
00695 
00696         // ****************************************************
00697         // ****************************************************
00698         if (antiflag>0)
00699         {
00700           // generate antithetical rv's
00701           generate_antithetical_rvs();
00702         }
00703         if (use_gauss_hermite>0)
00704         {
00705           if (gh)
00706           {
00707             delete gh;
00708             gh=0;
00709           }
00710           gh=new gauss_hermite_stuff(this,use_gauss_hermite,
00711             num_separable_calls,itmp);
00712         }
00713 
00714         if (block_diagonal_vch)
00715         {
00716           delete block_diagonal_vch;
00717           block_diagonal_vch=0;
00718         }
00719 
00720         block_diagonal_vch = new dvar3_array(1,num_separable_calls,
00721           1,itmp,1,itmp);
00722         if (block_diagonal_ch)
00723         {
00724           delete block_diagonal_ch;
00725           block_diagonal_ch=0;
00726         }
00727         block_diagonal_ch = new d3_array(1,num_separable_calls,
00728           1,itmp,1,itmp);
00729         if (block_diagonal_hessian)
00730         {
00731           delete block_diagonal_hessian;
00732           block_diagonal_hessian=0;
00733         }
00734         block_diagonal_hessian = new d3_array(1,num_separable_calls,
00735           1,itmp,1,itmp);
00736         if (block_diagonal_hessian ==0)
00737         {
00738           cerr << "error_allocating d3_array" << endl;
00739           ad_exit(1);
00740         }
00741         block_diagonal_re_list = new imatrix(1,num_separable_calls,
00742           1,itmp);
00743         if (block_diagonal_re_list ==0)
00744         {
00745           cerr << "error_allocating imatrix" << endl;
00746           ad_exit(1);
00747         }
00748         block_diagonal_fe_list = new imatrix(1,num_separable_calls,
00749           1,itmpf);
00750         if (block_diagonal_fe_list ==0)
00751         {
00752           cerr << "error_allocating imatrix" << endl;
00753           ad_exit(1);
00754         }
00755         // ****************************************************
00756         if (block_diagonal_Dux)
00757         {
00758           delete block_diagonal_Dux;
00759           block_diagonal_Dux=0;
00760         }
00761         block_diagonal_Dux = new d3_array(1,num_separable_calls,
00762           1,itmp,1,itmpf);
00763         if (block_diagonal_Dux ==0)
00764         {
00765           cerr << "error_allocating d3_array" << endl;
00766           ad_exit(1);
00767         }
00768 
00769         // ****************************************************
00770         // ****************************************************
00771         if (block_diagonal_vhessian)
00772         {
00773           delete block_diagonal_vhessian;
00774           block_diagonal_vhessian=0;
00775         }
00776         block_diagonal_vhessian = new dvar3_array(1,num_separable_calls,
00777           1,itmp,1,itmp);
00778         if (block_diagonal_vhessian ==0)
00779         {
00780           cerr << "error_allocating d3_array" << endl;
00781           ad_exit(1);
00782         }
00783 
00784         if (block_diagonal_vhessianadjoint)
00785         {
00786           delete block_diagonal_vhessianadjoint;
00787           block_diagonal_vhessianadjoint=0;
00788         }
00789         block_diagonal_vhessianadjoint = new d3_array(1,num_separable_calls,
00790           1,itmp,1,itmp);
00791         if (block_diagonal_vhessianadjoint ==0)
00792         {
00793           cerr << "error_allocating d3_array" << endl;
00794           ad_exit(1);
00795         }
00796       }
00797       funnel_init_var::lapprox=0;
00798       block_diagonal_flag=0;
00799       pen.deallocate();
00800     }
00801   }
00802 }
00803 
00808 void laplace_approximation_calculator::allocate_block_diagonal_stuff(void)
00809 {
00810   const ivector & itmp=(*num_local_re_array)(1,num_separable_calls);
00811   const ivector & itmpf=(*num_local_fixed_array)(1,num_separable_calls);
00812 
00813   // ****************************************************
00814   // ****************************************************
00815   if (block_diagonal_vch)
00816   {
00817     delete block_diagonal_vch;
00818     block_diagonal_vch=0;
00819   }
00820   block_diagonal_vch = new dvar3_array(1,num_separable_calls,
00821           1,itmp,1,itmp);
00822   if (block_diagonal_ch)
00823   {
00824     delete block_diagonal_ch;
00825     block_diagonal_ch=0;
00826   }
00827   block_diagonal_ch = new d3_array(1,num_separable_calls,
00828     1,itmp,1,itmp);
00829   if (block_diagonal_hessian)
00830   {
00831     delete block_diagonal_hessian;
00832     block_diagonal_hessian=0;
00833   }
00834   block_diagonal_hessian = new d3_array(1,num_separable_calls,
00835     1,itmp,1,itmp);
00836   if (block_diagonal_hessian ==0)
00837   {
00838     cerr << "error_allocating d3_array" << endl;
00839     ad_exit(1);
00840   }
00841   block_diagonal_re_list = new imatrix(1,num_separable_calls,
00842     1,itmp);
00843   if (block_diagonal_re_list ==0)
00844   {
00845     cerr << "error_allocating imatrix" << endl;
00846     ad_exit(1);
00847   }
00848   block_diagonal_fe_list = new imatrix(1,num_separable_calls,
00849     1,itmpf);
00850   if (block_diagonal_fe_list ==0)
00851   {
00852     cerr << "error_allocating imatrix" << endl;
00853     ad_exit(1);
00854   }
00855   // ****************************************************
00856   if (block_diagonal_Dux)
00857   {
00858     delete block_diagonal_Dux;
00859     block_diagonal_Dux=0;
00860   }
00861   block_diagonal_Dux = new d3_array(1,num_separable_calls,
00862     1,itmp,1,itmpf);
00863   if (block_diagonal_Dux ==0)
00864   {
00865     cerr << "error_allocating d3_array" << endl;
00866     ad_exit(1);
00867   }
00868 
00869   // ****************************************************
00870   // ****************************************************
00871   if (block_diagonal_vhessian)
00872   {
00873     delete block_diagonal_vhessian;
00874     block_diagonal_vhessian=0;
00875   }
00876   block_diagonal_vhessian = new dvar3_array(1,num_separable_calls,
00877     1,itmp,1,itmp);
00878   if (block_diagonal_vhessian ==0)
00879   {
00880     cerr << "error_allocating d3_array" << endl;
00881     ad_exit(1);
00882   }
00883 
00884   if (block_diagonal_vhessianadjoint)
00885   {
00886     delete block_diagonal_vhessianadjoint;
00887     block_diagonal_vhessianadjoint=0;
00888   }
00889   block_diagonal_vhessianadjoint = new d3_array(1,num_separable_calls,
00890     1,itmp,1,itmp);
00891   if (block_diagonal_vhessianadjoint ==0)
00892   {
00893     cerr << "error_allocating d3_array" << endl;
00894     ad_exit(1);
00895   }
00896 }
00897 
00902 void save_number_of_local_effects(int num_separable_calls,
00903   ivector ** num_local_re_array, ivector ** num_local_fixed_array,
00904   int num_local_re,int num_fixed_effects)
00905   //ivector& lre_index,ivector& lfe_index)
00906 {
00907   if (*num_local_re_array==0)
00908   {
00909     *num_local_re_array=new ivector(1,1000);
00910     if (*num_local_re_array==0)
00911     {
00912       cerr << "error allocating ivector" << endl;
00913       ad_exit(1);
00914     }
00915   }
00916 
00917   if (num_separable_calls> (*num_local_re_array)->indexmax())
00918   {
00919     if (num_separable_calls != (*num_local_re_array)->indexmax()+1)
00920     {
00921       cerr << "this can't happen" << endl;
00922       ad_exit(1);
00923     }
00924     int old_max=(*num_local_re_array)->indexmax();
00925     int new_max=old_max+100+10*sqrt(double(old_max));
00926     ivector tmp(1,old_max);
00927     tmp=(**num_local_re_array);
00928     (*num_local_re_array)=new ivector(1,new_max);
00929 
00930     delete *num_local_re_array;
00931     *num_local_re_array=new ivector(1,new_max);
00932     if (*num_local_re_array==0)
00933     {
00934       cerr << "error allocating ivector" << endl;
00935       ad_exit(1);
00936     }
00937     (**num_local_re_array)(1,old_max)=tmp;
00938   }
00939   (**num_local_re_array)(num_separable_calls)=num_local_re;
00940 
00941  //***********************************************************
00942  //***********************************************************
00943  //***********************************************************
00944  //***********************************************************
00945 
00946   if (*num_local_fixed_array==0)
00947   {
00948     *num_local_fixed_array=new ivector(1,1000);
00949     if (*num_local_fixed_array==0)
00950     {
00951       cerr << "error allocating ivector" << endl;
00952       ad_exit(1);
00953     }
00954   }
00955 
00956   if (num_separable_calls> (*num_local_fixed_array)->indexmax())
00957   {
00958     if (num_separable_calls != (*num_local_fixed_array)->indexmax()+1)
00959     {
00960       cerr << "this can't happen" << endl;
00961       ad_exit(1);
00962     }
00963     int old_max=(*num_local_fixed_array)->indexmax();
00964     int new_max=old_max+100+10*sqrt(double(old_max));
00965     ivector tmp(1,old_max);
00966     tmp=(**num_local_fixed_array);
00967     (*num_local_fixed_array)=new ivector(1,new_max);
00968 
00969     delete *num_local_fixed_array;
00970     *num_local_fixed_array=new ivector(1,new_max);
00971     if (*num_local_fixed_array==0)
00972     {
00973       cerr << "error allocating ivector" << endl;
00974       ad_exit(1);
00975     }
00976     (**num_local_fixed_array)(1,old_max)=tmp;
00977   }
00978   (**num_local_fixed_array)(num_separable_calls)=num_fixed_effects;
00979 }
00980 
00981 
00986 void laplace_approximation_calculator::
00987   do_separable_stuff_hessian_type_information(void)
00988 {
00989   df1b2_gradlist::set_no_derivatives();
00990 
00991   imatrix& list=*funnel_init_var::plist;
00992   int num_local_re=0;
00993   int num_fixed_effects=0;
00994 
00995   ivector lre_index(1,funnel_init_var::num_active_parameters);
00996   ivector lfe_index(1,funnel_init_var::num_active_parameters);
00997 
00998   for (int i=1;i<=funnel_init_var::num_active_parameters;i++)
00999   {
01000     if (list(i,1)>xsize)
01001     {
01002       lre_index(++num_local_re)=i;
01003     }
01004     else if (list(i,1)>0)
01005     {
01006       lfe_index(++num_fixed_effects)=i;
01007     }
01008   }
01009 
01010   //if (num_local_re > 0)
01011   {
01012     switch(hesstype)
01013     {
01014     case 3:
01015       num_separable_calls++;
01016       save_number_of_local_effects(num_separable_calls,
01017         &num_local_re_array,&num_local_fixed_array,num_local_re,
01018         num_fixed_effects); //,lre_index,lfe_index);
01019       for (int i=1;i<=num_local_re;i++)
01020       {
01021         int lrei=lre_index(i);
01022         for (int j=1;j<=num_local_re;j++)
01023         {
01024           int lrej=lre_index(j);
01025           int i1=list(lrei,1)-xsize;
01026           int j1=list(lrej,1)-xsize;
01027           if (i1>=j1)
01028           {
01029             //(*bHess)(i1,j1)+=locy(i2).u_bar[j2-1];
01030             int w=i1-j1+1;
01031             if (bw<w) bw=w;
01032           }
01033         }
01034       }
01035 
01036       if (sparse_hessian_flag)
01037       {
01038         if (allocated(Hess))
01039         {
01040           Hess.deallocate();
01041          /*
01042           if (Hess.indexmax() !=usize)
01043           {
01044             Hess.deallocate();
01045             Hess.allocate(1,usize,1,usize);
01046             Hess.initialize();
01047           }
01048           */
01049         }
01050        /*
01051         else
01052         {
01053           Hess.allocate(1,usize,1,usize);
01054           Hess.initialize();
01055         }
01056         */
01057         int dim= num_local_re*num_local_re;
01058         imatrix tmp(1,2,1,dim);
01059 
01060         int ii=0;
01061         for (int i=1;i<=num_local_re;i++)
01062         {
01063           int lrei=lre_index(i);
01064           for (int j=1;j<=num_local_re;j++)
01065           {
01066             int lrej=lre_index(j);
01067             int i1=list(lrei,1)-xsize;
01068             int j1=list(lrej,1)-xsize;
01069             if (i1<=0)
01070             {
01071               cout << "cant happen?" << endl;
01072             }
01073             if (i1<=j1)
01074             {
01075               //Hess(i1,j1)=1;
01076               ii++;
01077               tmp(1,ii)=i1;
01078               tmp(2,ii)=j1;
01079               //(*triplet_information)(1,num_separable_calls)=i1;
01080               //(*triplet_information)(2,num_separable_calls)=j1;
01081             }
01082           }
01083         }
01084 
01085         if (allocated((*triplet_information)(num_separable_calls)))
01086         {
01087           (*triplet_information)(num_separable_calls).deallocate();
01088         }
01089         if (ii>0)
01090         {
01091           (*triplet_information)(num_separable_calls).allocate(1,2,1,ii);
01092           (*triplet_information)(num_separable_calls)(1)=tmp(1)(1,ii);
01093           (*triplet_information)(num_separable_calls)(2)=tmp(2)(1,ii);
01094         }
01095       }
01096     }
01097   }
01098   if (derindex)
01099   {
01100     if (num_separable_calls> derindex->indexmax())
01101     {
01102        cerr << "Need to increase the maximum number of separable calls allowed"
01103             << " to at least " << num_separable_calls
01104             << endl << "Current value is " <<  derindex->indexmax() << endl;
01105        cerr << "Use the -ndi N command line option" << endl;
01106        ad_exit(1);
01107     }
01108 
01109     if (allocated((*derindex)(num_separable_calls)))
01110       (*derindex)(num_separable_calls).deallocate();
01111     (*derindex)(num_separable_calls).allocate(1,num_local_re);
01112     for (int i=1;i<=num_local_re;i++)
01113     {
01114       int lrei=lre_index(i);
01115       int i1=list(lrei,1)-xsize;
01116       //int i1=list(lrei,1);
01117       (*derindex)(num_separable_calls)(i)=i1;
01118     }
01119   }
01120 
01121   f1b2gradlist->reset();
01122   f1b2gradlist->list.initialize();
01123   f1b2gradlist->list2.initialize();
01124   f1b2gradlist->list3.initialize();
01125   f1b2gradlist->nlist.initialize();
01126   f1b2gradlist->nlist2.initialize();
01127   f1b2gradlist->nlist3.initialize();
01128   funnel_init_var::num_vars=0;
01129   funnel_init_var::num_active_parameters=0;
01130   funnel_init_var::num_inactive_vars=0;
01131 }
01132 
01137 imatrix laplace_approximation_calculator::check_sparse_matrix_structure(void)
01138 {
01139   ivector rowsize(1,usize);
01140   rowsize.initialize();
01141 
01142   imatrix tmp(1,usize,1,usize);
01143   tmp.initialize();
01144   for (int ii=1;ii<=num_separable_calls;ii++)
01145   {
01146     if (allocated((*derindex)(ii)))
01147     {
01148       ivector   iv = sort((*derindex)(ii));
01149       int n=iv.indexmax();
01150       if (n>1)                     // so we have off diagonal elements
01151       {
01152         for (int i=1;i<=n;i++)
01153         {
01154           int row=iv(i);
01155           for (int j=1;j<=n;j++)
01156           {
01157             if (i != j)
01158             {
01159               int col=iv(j);
01160               int  foundmatch=0;
01161               for (int k=1;k<=rowsize(row);k++)
01162               {
01163                 if (tmp(row,k)==col)
01164                 {
01165                   foundmatch=1;
01166                   break;
01167                 }
01168               }
01169               if (foundmatch==0)  // add a new element to tmp(row)
01170               {
01171                 rowsize(row)++;
01172                 if (rowsize(row)> tmp(row).indexmax())
01173                 {
01174                   tmp(row).reallocate(2);  // double the size
01175                 }
01176                 tmp(row,rowsize(row))=col;
01177               }
01178             }
01179           }
01180         }
01181       }
01182     }
01183   }
01184   imatrix M(1,usize,1,rowsize);
01185   ofstream ofs("sparseness.info");
01186   ofs << "# Number of parameters" << endl;
01187   ofs << usize << endl;
01188   ofs << "# Number of off diagonal elements in each row" << endl;
01189   ofs << rowsize << endl;
01190   ofs << "# The nonempty rows of M" << endl;
01191   for (int i=1;i<=usize;i++)
01192   {
01193     if (rowsize(i)>0)
01194     {
01195       M(i)=sort(tmp(i)(1,rowsize(i)));
01196       ofs << setw(4) << M(i) << endl;
01197     }
01198   }
01199   exit(1);
01200   return M;
01201 }
01202 #endif // if defined(USE_LAPLACE)