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