ADMB Documentation  11.1.2547
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2lp8.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2lp8.cpp 2496 2014-10-23 18:07:23Z 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         block_diagonal_vch = new dvar3_array(1,num_separable_calls,
00723           1,itmp,1,itmp);
00724 
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 
00733         if (block_diagonal_hessian)
00734         {
00735           delete block_diagonal_hessian;
00736           block_diagonal_hessian=0;
00737         }
00738         block_diagonal_hessian = new d3_array(1,num_separable_calls,
00739           1,itmp,1,itmp);
00740         if (block_diagonal_hessian ==0)
00741         {
00742           cerr << "error_allocating d3_array" << endl;
00743           ad_exit(1);
00744         }
00745 
00746         if (block_diagonal_re_list)
00747         {
00748           delete block_diagonal_re_list;
00749           block_diagonal_re_list = 0;
00750         }
00751         block_diagonal_re_list = new imatrix(1,num_separable_calls,
00752           1,itmp);
00753         if (block_diagonal_re_list == 0)
00754         {
00755           cerr << "error_allocating imatrix" << endl;
00756           ad_exit(1);
00757         }
00758 
00759         if (block_diagonal_fe_list)
00760         {
00761           delete block_diagonal_fe_list;
00762           block_diagonal_fe_list = 0;
00763         }
00764         block_diagonal_fe_list = new imatrix(1,num_separable_calls,
00765           1,itmpf);
00766         if (block_diagonal_fe_list ==0)
00767         {
00768           cerr << "error_allocating imatrix" << endl;
00769           ad_exit(1);
00770         }
00771 
00772         // ****************************************************
00773         if (block_diagonal_Dux)
00774         {
00775           delete block_diagonal_Dux;
00776           block_diagonal_Dux=0;
00777         }
00778         block_diagonal_Dux = new d3_array(1,num_separable_calls,
00779           1,itmp,1,itmpf);
00780         if (block_diagonal_Dux ==0)
00781         {
00782           cerr << "error_allocating d3_array" << endl;
00783           ad_exit(1);
00784         }
00785 
00786         // ****************************************************
00787         // ****************************************************
00788         if (block_diagonal_vhessian)
00789         {
00790           delete block_diagonal_vhessian;
00791           block_diagonal_vhessian=0;
00792         }
00793         block_diagonal_vhessian = new dvar3_array(1,num_separable_calls,
00794           1,itmp,1,itmp);
00795         if (block_diagonal_vhessian ==0)
00796         {
00797           cerr << "error_allocating d3_array" << endl;
00798           ad_exit(1);
00799         }
00800 
00801         if (block_diagonal_vhessianadjoint)
00802         {
00803           delete block_diagonal_vhessianadjoint;
00804           block_diagonal_vhessianadjoint=0;
00805         }
00806         block_diagonal_vhessianadjoint = new d3_array(1,num_separable_calls,
00807           1,itmp,1,itmp);
00808         if (block_diagonal_vhessianadjoint ==0)
00809         {
00810           cerr << "error_allocating d3_array" << endl;
00811           ad_exit(1);
00812         }
00813       }
00814       funnel_init_var::lapprox=0;
00815       block_diagonal_flag=0;
00816       pen.deallocate();
00817     }
00818   }
00819 }
00820 
00821 #include <cassert>
00826 void laplace_approximation_calculator::allocate_block_diagonal_stuff(void)
00827 {
00828   const ivector & itmp=(*num_local_re_array)(1,num_separable_calls);
00829   const ivector & itmpf=(*num_local_fixed_array)(1,num_separable_calls);
00830 
00831   // ****************************************************
00832   // ****************************************************
00833   if (block_diagonal_vch)
00834   {
00835     delete block_diagonal_vch;
00836     block_diagonal_vch=0;
00837   }
00838   block_diagonal_vch = new dvar3_array(1,num_separable_calls,
00839           1,itmp,1,itmp);
00840   if (block_diagonal_ch)
00841   {
00842     delete block_diagonal_ch;
00843     block_diagonal_ch=0;
00844   }
00845   block_diagonal_ch = new d3_array(1,num_separable_calls,
00846     1,itmp,1,itmp);
00847   if (block_diagonal_hessian)
00848   {
00849     delete block_diagonal_hessian;
00850     block_diagonal_hessian=0;
00851   }
00852   block_diagonal_hessian = new d3_array(1,num_separable_calls,
00853     1,itmp,1,itmp);
00854   if (block_diagonal_hessian ==0)
00855   {
00856     cerr << "error_allocating d3_array" << endl;
00857     ad_exit(1);
00858   }
00859   if (block_diagonal_re_list)
00860   {
00861     delete block_diagonal_re_list;
00862     block_diagonal_re_list=0;
00863   }
00864   block_diagonal_re_list = new imatrix(1,num_separable_calls,
00865     1,itmp);
00866   if (block_diagonal_re_list ==0)
00867   {
00868     cerr << "error_allocating imatrix" << endl;
00869     ad_exit(1);
00870   }
00871   if (block_diagonal_fe_list)
00872   {
00873     delete block_diagonal_fe_list;
00874     block_diagonal_fe_list=0;
00875   }
00876   block_diagonal_fe_list = new imatrix(1,num_separable_calls,
00877     1,itmpf);
00878   if (block_diagonal_fe_list ==0)
00879   {
00880     cerr << "error_allocating imatrix" << endl;
00881     ad_exit(1);
00882   }
00883   // ****************************************************
00884   if (block_diagonal_Dux)
00885   {
00886     delete block_diagonal_Dux;
00887     block_diagonal_Dux=0;
00888   }
00889   block_diagonal_Dux = new d3_array(1,num_separable_calls,
00890     1,itmp,1,itmpf);
00891   if (block_diagonal_Dux ==0)
00892   {
00893     cerr << "error_allocating d3_array" << endl;
00894     ad_exit(1);
00895   }
00896 
00897   // ****************************************************
00898   // ****************************************************
00899   if (block_diagonal_vhessian)
00900   {
00901     delete block_diagonal_vhessian;
00902     block_diagonal_vhessian=0;
00903   }
00904   block_diagonal_vhessian = new dvar3_array(1,num_separable_calls,
00905     1,itmp,1,itmp);
00906   if (block_diagonal_vhessian ==0)
00907   {
00908     cerr << "error_allocating d3_array" << endl;
00909     ad_exit(1);
00910   }
00911 
00912   if (block_diagonal_vhessianadjoint)
00913   {
00914     delete block_diagonal_vhessianadjoint;
00915     block_diagonal_vhessianadjoint=0;
00916   }
00917   block_diagonal_vhessianadjoint = new d3_array(1,num_separable_calls,
00918     1,itmp,1,itmp);
00919   if (block_diagonal_vhessianadjoint ==0)
00920   {
00921     cerr << "error_allocating d3_array" << endl;
00922     ad_exit(1);
00923   }
00924 }
00925 
00929 void save_number_of_local_effects(int num_separable_calls,
00930   ivector ** num_local_re_array, ivector ** num_local_fixed_array,
00931   int num_local_re,int num_fixed_effects)
00932   //ivector& lre_index,ivector& lfe_index)
00933 {
00934   if (*num_local_re_array==0)
00935   {
00936     *num_local_re_array=new ivector(1,1000);
00937     if (*num_local_re_array==0)
00938     {
00939       cerr << "error allocating ivector" << endl;
00940       ad_exit(1);
00941     }
00942   }
00943 
00944   if (num_separable_calls> (*num_local_re_array)->indexmax())
00945   {
00946     if (num_separable_calls != (*num_local_re_array)->indexmax()+1)
00947     {
00948       cerr << "this can't happen" << endl;
00949       ad_exit(1);
00950     }
00951     int old_max=(*num_local_re_array)->indexmax();
00952 #ifdef OPT_LIB
00953     int new_max=old_max+100+(int)(10.0*sqrt(double(old_max)));
00954 #else
00955     double sqrt_oldmax = 10.0 * sqrt(double(old_max));
00956     assert(sqrt_oldmax <= INT_MAX);
00957     int new_max=old_max+100+(int)sqrt_oldmax;
00958 #endif
00959     ivector tmp(1,old_max);
00960     tmp=(**num_local_re_array);
00961     (*num_local_re_array)=new ivector(1,new_max);
00962 
00963     delete *num_local_re_array;
00964     *num_local_re_array=new ivector(1,new_max);
00965     if (*num_local_re_array==0)
00966     {
00967       cerr << "error allocating ivector" << endl;
00968       ad_exit(1);
00969     }
00970     (**num_local_re_array)(1,old_max)=tmp;
00971   }
00972   (**num_local_re_array)(num_separable_calls)=num_local_re;
00973 
00974  //***********************************************************
00975  //***********************************************************
00976  //***********************************************************
00977  //***********************************************************
00978 
00979   if (*num_local_fixed_array==0)
00980   {
00981     *num_local_fixed_array=new ivector(1,1000);
00982     if (*num_local_fixed_array==0)
00983     {
00984       cerr << "error allocating ivector" << endl;
00985       ad_exit(1);
00986     }
00987   }
00988 
00989   if (num_separable_calls> (*num_local_fixed_array)->indexmax())
00990   {
00991     if (num_separable_calls != (*num_local_fixed_array)->indexmax()+1)
00992     {
00993       cerr << "this can't happen" << endl;
00994       ad_exit(1);
00995     }
00996     int old_max=(*num_local_fixed_array)->indexmax();
00997 #ifdef OPT_LIB
00998     int new_max=old_max+100+(int)(10.0*sqrt(double(old_max)));
00999 #else
01000     double sqrt_oldmax = 10.0 * sqrt(double(old_max));
01001     assert(sqrt_oldmax <= INT_MAX);
01002     int new_max=old_max+100+(int)sqrt_oldmax;
01003 #endif
01004     ivector tmp(1,old_max);
01005     tmp=(**num_local_fixed_array);
01006     (*num_local_fixed_array)=new ivector(1,new_max);
01007 
01008     delete *num_local_fixed_array;
01009     *num_local_fixed_array=new ivector(1,new_max);
01010     if (*num_local_fixed_array==0)
01011     {
01012       cerr << "error allocating ivector" << endl;
01013       ad_exit(1);
01014     }
01015     (**num_local_fixed_array)(1,old_max)=tmp;
01016   }
01017   (**num_local_fixed_array)(num_separable_calls)=num_fixed_effects;
01018 }
01019 
01020 
01025 void laplace_approximation_calculator::
01026   do_separable_stuff_hessian_type_information(void)
01027 {
01028   df1b2_gradlist::set_no_derivatives();
01029 
01030   imatrix& list=*funnel_init_var::plist;
01031   int num_local_re=0;
01032   int num_fixed_effects=0;
01033 
01034   ivector lre_index(1,funnel_init_var::num_active_parameters);
01035   ivector lfe_index(1,funnel_init_var::num_active_parameters);
01036 
01037   for (int i=1;i<=funnel_init_var::num_active_parameters;i++)
01038   {
01039     if (list(i,1)>xsize)
01040     {
01041       lre_index(++num_local_re)=i;
01042     }
01043     else if (list(i,1)>0)
01044     {
01045       lfe_index(++num_fixed_effects)=i;
01046     }
01047   }
01048 
01049   //if (num_local_re > 0)
01050   {
01051     switch(hesstype)
01052     {
01053     case 3:
01054       num_separable_calls++;
01055       save_number_of_local_effects(num_separable_calls,
01056         &num_local_re_array,&num_local_fixed_array,num_local_re,
01057         num_fixed_effects); //,lre_index,lfe_index);
01058       for (int i=1;i<=num_local_re;i++)
01059       {
01060         int lrei=lre_index(i);
01061         for (int j=1;j<=num_local_re;j++)
01062         {
01063           int lrej=lre_index(j);
01064           int i1=list(lrei,1)-xsize;
01065           int j1=list(lrej,1)-xsize;
01066           if (i1>=j1)
01067           {
01068             //(*bHess)(i1,j1)+=locy(i2).u_bar[j2-1];
01069             int w=i1-j1+1;
01070             if (bw<w) bw=w;
01071           }
01072         }
01073       }
01074 
01075       if (sparse_hessian_flag)
01076       {
01077         if (allocated(Hess))
01078         {
01079           Hess.deallocate();
01080          /*
01081           if (Hess.indexmax() !=usize)
01082           {
01083             Hess.deallocate();
01084             Hess.allocate(1,usize,1,usize);
01085             Hess.initialize();
01086           }
01087           */
01088         }
01089        /*
01090         else
01091         {
01092           Hess.allocate(1,usize,1,usize);
01093           Hess.initialize();
01094         }
01095         */
01096         int dim= num_local_re*num_local_re;
01097         imatrix tmp(1,2,1,dim);
01098 
01099         int ii=0;
01100         for (int i=1;i<=num_local_re;i++)
01101         {
01102           int lrei=lre_index(i);
01103           for (int j=1;j<=num_local_re;j++)
01104           {
01105             int lrej=lre_index(j);
01106             int i1=list(lrei,1)-xsize;
01107             int j1=list(lrej,1)-xsize;
01108             if (i1<=0)
01109             {
01110               cout << "cant happen?" << endl;
01111             }
01112             if (i1<=j1)
01113             {
01114               //Hess(i1,j1)=1;
01115               ii++;
01116               tmp(1,ii)=i1;
01117               tmp(2,ii)=j1;
01118               //(*triplet_information)(1,num_separable_calls)=i1;
01119               //(*triplet_information)(2,num_separable_calls)=j1;
01120             }
01121           }
01122         }
01123 
01124         if (allocated((*triplet_information)(num_separable_calls)))
01125         {
01126           (*triplet_information)(num_separable_calls).deallocate();
01127         }
01128         if (ii>0)
01129         {
01130           (*triplet_information)(num_separable_calls).allocate(1,2,1,ii);
01131           (*triplet_information)(num_separable_calls)(1)=tmp(1)(1,ii);
01132           (*triplet_information)(num_separable_calls)(2)=tmp(2)(1,ii);
01133         }
01134       }
01135     }
01136   }
01137   if (derindex)
01138   {
01139     if (num_separable_calls> derindex->indexmax())
01140     {
01141        cerr << "Need to increase the maximum number of separable calls allowed"
01142             << " to at least " << num_separable_calls
01143             << endl << "Current value is " <<  derindex->indexmax() << endl;
01144        cerr << "Use the -ndi N command line option" << endl;
01145        ad_exit(1);
01146     }
01147 
01148     if (allocated((*derindex)(num_separable_calls)))
01149       (*derindex)(num_separable_calls).deallocate();
01150     (*derindex)(num_separable_calls).allocate(1,num_local_re);
01151     for (int i=1;i<=num_local_re;i++)
01152     {
01153       int lrei=lre_index(i);
01154       int i1=list(lrei,1)-xsize;
01155       //int i1=list(lrei,1);
01156       (*derindex)(num_separable_calls)(i)=i1;
01157     }
01158   }
01159 
01160   f1b2gradlist->reset();
01161   f1b2gradlist->list.initialize();
01162   f1b2gradlist->list2.initialize();
01163   f1b2gradlist->list3.initialize();
01164   f1b2gradlist->nlist.initialize();
01165   f1b2gradlist->nlist2.initialize();
01166   f1b2gradlist->nlist3.initialize();
01167   funnel_init_var::num_vars=0;
01168   funnel_init_var::num_active_parameters=0;
01169   funnel_init_var::num_inactive_vars=0;
01170 }
01171 
01176 imatrix laplace_approximation_calculator::check_sparse_matrix_structure(void)
01177 {
01178   ivector rowsize(1,usize);
01179   rowsize.initialize();
01180 
01181   imatrix tmp(1,usize,1,usize);
01182   tmp.initialize();
01183   for (int ii=1;ii<=num_separable_calls;ii++)
01184   {
01185     if (allocated((*derindex)(ii)))
01186     {
01187       ivector   iv = sort((*derindex)(ii));
01188       int n=iv.indexmax();
01189       if (n>1)                     // so we have off diagonal elements
01190       {
01191         for (int i=1;i<=n;i++)
01192         {
01193           int row=iv(i);
01194           for (int j=1;j<=n;j++)
01195           {
01196             if (i != j)
01197             {
01198               int col=iv(j);
01199               int  foundmatch=0;
01200               for (int k=1;k<=rowsize(row);k++)
01201               {
01202                 if (tmp(row,k)==col)
01203                 {
01204                   foundmatch=1;
01205                   break;
01206                 }
01207               }
01208               if (foundmatch==0)  // add a new element to tmp(row)
01209               {
01210                 rowsize(row)++;
01211                 if (rowsize(row)> tmp(row).indexmax())
01212                 {
01213                   tmp(row).reallocate(2);  // double the size
01214                 }
01215                 tmp(row,rowsize(row))=col;
01216               }
01217             }
01218           }
01219         }
01220       }
01221     }
01222   }
01223   imatrix M(1,usize,1,rowsize);
01224   ofstream ofs("sparseness.info");
01225   ofs << "# Number of parameters" << endl;
01226   ofs << usize << endl;
01227   ofs << "# Number of off diagonal elements in each row" << endl;
01228   ofs << rowsize << endl;
01229   ofs << "# The nonempty rows of M" << endl;
01230   for (int i=1;i<=usize;i++)
01231   {
01232     if (rowsize(i)>0)
01233     {
01234       M(i)=sort(tmp(i)(1,rowsize(i)));
01235       ofs << setw(4) << M(i) << endl;
01236     }
01237   }
01238   exit(1);
01239   return M;
01240 }