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