ADMB Documentation  11.1x.2706
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2lp8.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2lp8.cpp 2615 2014-11-10 22:35:40Z 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   for (int i=2;i<=num_separable_calls;i++)
00142   {
00143     if (itmp(i) != itmp(i-1))
00144     {
00145       cerr << "At present can only use antithetical rv's when "
00146               "all separable calls are the same size" << endl;
00147       ad_exit(1);
00148     }
00149   }
00150   int n=itmp(1);
00151   int samplesize=num_importance_samples;
00152 
00153   // mesh size
00154   double delta=0.01;
00155   // maximum of distribution is near here
00156   double mid=sqrt(double(n-1));
00157   dmatrix weights(1,2*n,1,2);
00158   double spread=15;
00159   if (mid-spread<=0.001)
00160     spread=mid-0.1;
00161   double ssum=0.0;
00162   double x=0.0;
00163   double tmax=(n-1)*log(mid)-0.5*mid*mid;
00164   for (x=mid-spread;x<=mid+spread;x+=delta)
00165   {
00166     ssum+=exp((n-1)*log(x)-0.5*x*x-tmax);
00167   }
00168   double tsum=0;
00169   dvector dist(1,samplesize+1);
00170   dist.initialize();
00171   int is=0;
00172   int ii;
00173   for (x=mid-spread;x<=mid+spread;x+=delta)
00174   {
00175     tsum+=exp((n-1)*log(x)-0.5*x*x-tmax)/ssum*samplesize;
00176     int ns=int(tsum);
00177     for (ii=1;ii<=ns;ii++)
00178     {
00179       dist(++is)=x;
00180     }
00181     tsum-=ns;
00182   }
00183   if (is==samplesize-1)
00184   {
00185     dist(samplesize)=mid;
00186   }
00187   else if (is<samplesize-1)
00188   {
00189     cerr << "This can't happen" << endl;
00190     exit(1);
00191   }
00192 
00193   // get random numbers
00194 
00195   random_number_generator rng(rseed);
00196   if (antiepsilon)
00197   {
00198     if (allocated(*antiepsilon))
00199     {
00200       delete antiepsilon;
00201       antiepsilon=0;
00202     }
00203   }
00204   antiepsilon=new dmatrix(1,samplesize,1,n);
00205   dmatrix & M=*antiepsilon;
00206   M.fill_randn(rng);
00207 
00208   for (int i=1;i<=samplesize;i++)
00209   {
00210     M(i)=M(i)/norm(M(i));
00211   }
00212   int nvar=(samplesize-1)*n;
00213   independent_variables xx(1,nvar);
00214   ii=0;
00215   for (int i=2;i<=samplesize;i++)
00216   {
00217     for (int j=1;j<=n;j++)
00218     {
00219       xx(++ii)=M(i,j);
00220     }
00221   }
00222 
00223   fmmt1 fmc(nvar,5);
00224   //fmm fmc(nvar,5);
00225   fmc.noprintx=1;
00226   fmc.iprint=10;
00227   fmc.maxfn=2500;
00228   fmc.crit=1.e-6;
00229 
00230   double f;
00231   double fbest=1.e+50;;
00232   dvector g(1,nvar);
00233   dvector gbest(1,nvar);
00234   dvector xbest(1,nvar);
00235 
00236   gbest.fill_seqadd(1.e+50,0.);
00237   {
00238     while (fmc.ireturn>=0)
00239     {
00240       //int badflag=0;
00241       fmc.fmin(f,xx,g);
00242       if (fmc.ihang)
00243       {
00244         //int hang_flag=fmc.ihang;
00245         //double maxg=max(g);
00246         //double ccrit=fmc.crit;
00247         //int current_ifn=fmc.ifn;
00248       }
00249       if (fmc.ireturn>0)
00250       {
00251          f=fcomp1(xx,dist,samplesize,n,g,M);
00252          if (f < fbest)
00253          {
00254            fbest=f;
00255            gbest=g;
00256            xbest=xx;
00257          }
00258        }
00259      }
00260      xx=xbest;
00261   }
00262   ii=0;
00263   for (int i=2;i<=samplesize;i++)
00264   {
00265     for (int j=1;j<=n;j++)
00266     {
00267       M(i,j)=xx(++ii);
00268     }
00269   }
00270   for (int i=1;i<=samplesize;i++)
00271   {
00272     M(i)*=dist(i)/norm(M(i));
00273   }
00274 }
00275 
00280 double fcomp1(dvector x,dvector d,int samplesize,int n,dvector & g,
00281   dmatrix& M)
00282 {
00283   dmatrix VM(1,samplesize,1,n);
00284   dmatrix C(1,samplesize,1,samplesize);
00285   dmatrix VM0(1,samplesize,1,n);
00286   dvector N(1,samplesize);
00287   dmatrix dfVM(1,samplesize,1,n);
00288   dmatrix dfVM0(1,samplesize,1,n);
00289   dvector dfN(1,samplesize);
00290   dfVM.initialize();
00291   dfVM0.initialize();
00292   dfN.initialize();
00293 
00294   double f=0.0;
00295   int ii=0;
00296   VM0(1)=M(1);
00297   for (int i=2;i<=samplesize;i++)
00298   {
00299     for (int j=1;j<=n;j++)
00300     {
00301       VM0(i,j)=x(++ii);
00302     }
00303   }
00304   for (int i=1;i<=samplesize;i++)
00305   {
00306     N(i)=norm(VM0(i));
00307     VM(i)=VM0(i)*(d(i)/N(i));
00308   }
00309   for (int 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 (int 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 (int 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 (int i=2;i<=samplesize;i++)
00346   {
00347     for (int 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         block_diagonal_vch = new dvar3_array(1,num_separable_calls,
00719           1,itmp,1,itmp);
00720 
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 
00729         if (block_diagonal_hessian)
00730         {
00731           delete block_diagonal_hessian;
00732           block_diagonal_hessian=0;
00733         }
00734         block_diagonal_hessian = new d3_array(1,num_separable_calls,
00735           1,itmp,1,itmp);
00736         if (block_diagonal_hessian ==0)
00737         {
00738           cerr << "error_allocating d3_array" << endl;
00739           ad_exit(1);
00740         }
00741 
00742         if (block_diagonal_re_list)
00743         {
00744           delete block_diagonal_re_list;
00745           block_diagonal_re_list = 0;
00746         }
00747         block_diagonal_re_list = new imatrix(1,num_separable_calls,
00748           1,itmp);
00749         if (block_diagonal_re_list == 0)
00750         {
00751           cerr << "error_allocating imatrix" << endl;
00752           ad_exit(1);
00753         }
00754 
00755         if (block_diagonal_fe_list)
00756         {
00757           delete block_diagonal_fe_list;
00758           block_diagonal_fe_list = 0;
00759         }
00760         block_diagonal_fe_list = new imatrix(1,num_separable_calls,
00761           1,itmpf);
00762         if (block_diagonal_fe_list ==0)
00763         {
00764           cerr << "error_allocating imatrix" << endl;
00765           ad_exit(1);
00766         }
00767 
00768         // ****************************************************
00769         if (block_diagonal_Dux)
00770         {
00771           delete block_diagonal_Dux;
00772           block_diagonal_Dux=0;
00773         }
00774         block_diagonal_Dux = new d3_array(1,num_separable_calls,
00775           1,itmp,1,itmpf);
00776         if (block_diagonal_Dux ==0)
00777         {
00778           cerr << "error_allocating d3_array" << endl;
00779           ad_exit(1);
00780         }
00781 
00782         // ****************************************************
00783         // ****************************************************
00784         if (block_diagonal_vhessian)
00785         {
00786           delete block_diagonal_vhessian;
00787           block_diagonal_vhessian=0;
00788         }
00789         block_diagonal_vhessian = new dvar3_array(1,num_separable_calls,
00790           1,itmp,1,itmp);
00791         if (block_diagonal_vhessian ==0)
00792         {
00793           cerr << "error_allocating d3_array" << endl;
00794           ad_exit(1);
00795         }
00796 
00797         if (block_diagonal_vhessianadjoint)
00798         {
00799           delete block_diagonal_vhessianadjoint;
00800           block_diagonal_vhessianadjoint=0;
00801         }
00802         block_diagonal_vhessianadjoint = new d3_array(1,num_separable_calls,
00803           1,itmp,1,itmp);
00804         if (block_diagonal_vhessianadjoint ==0)
00805         {
00806           cerr << "error_allocating d3_array" << endl;
00807           ad_exit(1);
00808         }
00809       }
00810       funnel_init_var::lapprox=0;
00811       block_diagonal_flag=0;
00812       pen.deallocate();
00813     }
00814   }
00815 }
00816 
00821 void laplace_approximation_calculator::allocate_block_diagonal_stuff(void)
00822 {
00823   const ivector & itmp=(*num_local_re_array)(1,num_separable_calls);
00824   const ivector & itmpf=(*num_local_fixed_array)(1,num_separable_calls);
00825 
00826   // ****************************************************
00827   // ****************************************************
00828   if (block_diagonal_vch)
00829   {
00830     delete block_diagonal_vch;
00831     block_diagonal_vch=0;
00832   }
00833   block_diagonal_vch = new dvar3_array(1,num_separable_calls,
00834           1,itmp,1,itmp);
00835   if (block_diagonal_ch)
00836   {
00837     delete block_diagonal_ch;
00838     block_diagonal_ch=0;
00839   }
00840   block_diagonal_ch = new d3_array(1,num_separable_calls,
00841     1,itmp,1,itmp);
00842   if (block_diagonal_hessian)
00843   {
00844     delete block_diagonal_hessian;
00845     block_diagonal_hessian=0;
00846   }
00847   block_diagonal_hessian = new d3_array(1,num_separable_calls,
00848     1,itmp,1,itmp);
00849   if (block_diagonal_hessian ==0)
00850   {
00851     cerr << "error_allocating d3_array" << endl;
00852     ad_exit(1);
00853   }
00854   if (block_diagonal_re_list)
00855   {
00856     delete block_diagonal_re_list;
00857     block_diagonal_re_list=0;
00858   }
00859   block_diagonal_re_list = new imatrix(1,num_separable_calls,
00860     1,itmp);
00861   if (block_diagonal_re_list ==0)
00862   {
00863     cerr << "error_allocating imatrix" << endl;
00864     ad_exit(1);
00865   }
00866   if (block_diagonal_fe_list)
00867   {
00868     delete block_diagonal_fe_list;
00869     block_diagonal_fe_list=0;
00870   }
00871   block_diagonal_fe_list = new imatrix(1,num_separable_calls,
00872     1,itmpf);
00873   if (block_diagonal_fe_list ==0)
00874   {
00875     cerr << "error_allocating imatrix" << endl;
00876     ad_exit(1);
00877   }
00878   // ****************************************************
00879   if (block_diagonal_Dux)
00880   {
00881     delete block_diagonal_Dux;
00882     block_diagonal_Dux=0;
00883   }
00884   block_diagonal_Dux = new d3_array(1,num_separable_calls,
00885     1,itmp,1,itmpf);
00886   if (block_diagonal_Dux ==0)
00887   {
00888     cerr << "error_allocating d3_array" << endl;
00889     ad_exit(1);
00890   }
00891 
00892   // ****************************************************
00893   // ****************************************************
00894   if (block_diagonal_vhessian)
00895   {
00896     delete block_diagonal_vhessian;
00897     block_diagonal_vhessian=0;
00898   }
00899   block_diagonal_vhessian = new dvar3_array(1,num_separable_calls,
00900     1,itmp,1,itmp);
00901   if (block_diagonal_vhessian ==0)
00902   {
00903     cerr << "error_allocating d3_array" << endl;
00904     ad_exit(1);
00905   }
00906 
00907   if (block_diagonal_vhessianadjoint)
00908   {
00909     delete block_diagonal_vhessianadjoint;
00910     block_diagonal_vhessianadjoint=0;
00911   }
00912   block_diagonal_vhessianadjoint = new d3_array(1,num_separable_calls,
00913     1,itmp,1,itmp);
00914   if (block_diagonal_vhessianadjoint ==0)
00915   {
00916     cerr << "error_allocating d3_array" << endl;
00917     ad_exit(1);
00918   }
00919 }
00920 
00924 void save_number_of_local_effects(int num_separable_calls,
00925   ivector ** num_local_re_array, ivector ** num_local_fixed_array,
00926   int num_local_re,int num_fixed_effects)
00927   //ivector& lre_index,ivector& lfe_index)
00928 {
00929   if (*num_local_re_array==0)
00930   {
00931     *num_local_re_array=new ivector(1,1000);
00932     if (*num_local_re_array==0)
00933     {
00934       cerr << "error allocating ivector" << endl;
00935       ad_exit(1);
00936     }
00937   }
00938 
00939   if (num_separable_calls> (*num_local_re_array)->indexmax())
00940   {
00941     if (num_separable_calls != (*num_local_re_array)->indexmax()+1)
00942     {
00943       cerr << "this can't happen" << endl;
00944       ad_exit(1);
00945     }
00946     int old_max=(*num_local_re_array)->indexmax();
00947 #ifdef OPT_LIB
00948     int new_max=old_max+100+(int)(10.0*sqrt(double(old_max)));
00949 #else
00950     double sqrt_oldmax = 10.0 * sqrt(double(old_max));
00951     assert(sqrt_oldmax <= INT_MAX);
00952     int new_max=old_max+100+(int)sqrt_oldmax;
00953 #endif
00954     ivector tmp(1,old_max);
00955     tmp=(**num_local_re_array);
00956     (*num_local_re_array)=new ivector(1,new_max);
00957 
00958     delete *num_local_re_array;
00959     *num_local_re_array=new ivector(1,new_max);
00960     if (*num_local_re_array==0)
00961     {
00962       cerr << "error allocating ivector" << endl;
00963       ad_exit(1);
00964     }
00965     (**num_local_re_array)(1,old_max)=tmp;
00966   }
00967   (**num_local_re_array)(num_separable_calls)=num_local_re;
00968 
00969  //***********************************************************
00970  //***********************************************************
00971  //***********************************************************
00972  //***********************************************************
00973 
00974   if (*num_local_fixed_array==0)
00975   {
00976     *num_local_fixed_array=new ivector(1,1000);
00977     if (*num_local_fixed_array==0)
00978     {
00979       cerr << "error allocating ivector" << endl;
00980       ad_exit(1);
00981     }
00982   }
00983 
00984   if (num_separable_calls> (*num_local_fixed_array)->indexmax())
00985   {
00986     if (num_separable_calls != (*num_local_fixed_array)->indexmax()+1)
00987     {
00988       cerr << "this can't happen" << endl;
00989       ad_exit(1);
00990     }
00991     int old_max=(*num_local_fixed_array)->indexmax();
00992 #ifdef OPT_LIB
00993     int new_max=old_max+100+(int)(10.0*sqrt(double(old_max)));
00994 #else
00995     double sqrt_oldmax = 10.0 * sqrt(double(old_max));
00996     assert(sqrt_oldmax <= INT_MAX);
00997     int new_max=old_max+100+(int)sqrt_oldmax;
00998 #endif
00999     ivector tmp(1,old_max);
01000     tmp=(**num_local_fixed_array);
01001     (*num_local_fixed_array)=new ivector(1,new_max);
01002 
01003     delete *num_local_fixed_array;
01004     *num_local_fixed_array=new ivector(1,new_max);
01005     if (*num_local_fixed_array==0)
01006     {
01007       cerr << "error allocating ivector" << endl;
01008       ad_exit(1);
01009     }
01010     (**num_local_fixed_array)(1,old_max)=tmp;
01011   }
01012   (**num_local_fixed_array)(num_separable_calls)=num_fixed_effects;
01013 }
01014 
01015 
01020 void laplace_approximation_calculator::
01021   do_separable_stuff_hessian_type_information(void)
01022 {
01023   df1b2_gradlist::set_no_derivatives();
01024 
01025   imatrix& list=*funnel_init_var::plist;
01026   int num_local_re=0;
01027   int num_fixed_effects=0;
01028 #ifndef OPT_LIB
01029   assert(funnel_init_var::num_active_parameters <= INT_MAX);
01030 #endif
01031   ivector lre_index(1, (int)funnel_init_var::num_active_parameters);
01032   ivector lfe_index(1, (int)funnel_init_var::num_active_parameters);
01033 
01034   for (int i=1;i<=(int)funnel_init_var::num_active_parameters;i++)
01035   {
01036     if (list(i,1)>xsize)
01037     {
01038       lre_index(++num_local_re)=i;
01039     }
01040     else if (list(i,1)>0)
01041     {
01042       lfe_index(++num_fixed_effects)=i;
01043     }
01044   }
01045 
01046   //if (num_local_re > 0)
01047   {
01048     switch(hesstype)
01049     {
01050     case 3:
01051       num_separable_calls++;
01052       save_number_of_local_effects(num_separable_calls,
01053         &num_local_re_array,&num_local_fixed_array,num_local_re,
01054         num_fixed_effects); //,lre_index,lfe_index);
01055       for (int i=1;i<=num_local_re;i++)
01056       {
01057         int lrei=lre_index(i);
01058         for (int j=1;j<=num_local_re;j++)
01059         {
01060           int lrej=lre_index(j);
01061           int i1=list(lrei,1)-xsize;
01062           int j1=list(lrej,1)-xsize;
01063           if (i1>=j1)
01064           {
01065             //(*bHess)(i1,j1)+=locy(i2).u_bar[j2-1];
01066             int w=i1-j1+1;
01067             if (bw<w) bw=w;
01068           }
01069         }
01070       }
01071 
01072       if (sparse_hessian_flag)
01073       {
01074         if (allocated(Hess))
01075         {
01076           Hess.deallocate();
01077          /*
01078           if (Hess.indexmax() !=usize)
01079           {
01080             Hess.deallocate();
01081             Hess.allocate(1,usize,1,usize);
01082             Hess.initialize();
01083           }
01084           */
01085         }
01086        /*
01087         else
01088         {
01089           Hess.allocate(1,usize,1,usize);
01090           Hess.initialize();
01091         }
01092         */
01093         int dim= num_local_re*num_local_re;
01094         imatrix tmp(1,2,1,dim);
01095 
01096         int ii=0;
01097         for (int i=1;i<=num_local_re;i++)
01098         {
01099           int lrei=lre_index(i);
01100           for (int j=1;j<=num_local_re;j++)
01101           {
01102             int lrej=lre_index(j);
01103             int i1=list(lrei,1)-xsize;
01104             int j1=list(lrej,1)-xsize;
01105             if (i1<=0)
01106             {
01107               cout << "cant happen?" << endl;
01108             }
01109             if (i1<=j1)
01110             {
01111               //Hess(i1,j1)=1;
01112               ii++;
01113               tmp(1,ii)=i1;
01114               tmp(2,ii)=j1;
01115               //(*triplet_information)(1,num_separable_calls)=i1;
01116               //(*triplet_information)(2,num_separable_calls)=j1;
01117             }
01118           }
01119         }
01120 
01121         if (allocated((*triplet_information)(num_separable_calls)))
01122         {
01123           (*triplet_information)(num_separable_calls).deallocate();
01124         }
01125         if (ii>0)
01126         {
01127           (*triplet_information)(num_separable_calls).allocate(1,2,1,ii);
01128           (*triplet_information)(num_separable_calls)(1)=tmp(1)(1,ii);
01129           (*triplet_information)(num_separable_calls)(2)=tmp(2)(1,ii);
01130         }
01131       }
01132     }
01133   }
01134   if (derindex)
01135   {
01136     if (num_separable_calls> derindex->indexmax())
01137     {
01138        cerr << "Need to increase the maximum number of separable calls allowed"
01139             << " to at least " << num_separable_calls
01140             << endl << "Current value is " <<  derindex->indexmax() << endl;
01141        cerr << "Use the -ndi N command line option" << endl;
01142        ad_exit(1);
01143     }
01144 
01145     if (allocated((*derindex)(num_separable_calls)))
01146       (*derindex)(num_separable_calls).deallocate();
01147     (*derindex)(num_separable_calls).allocate(1,num_local_re);
01148     for (int i=1;i<=num_local_re;i++)
01149     {
01150       int lrei=lre_index(i);
01151       int i1=list(lrei,1)-xsize;
01152       //int i1=list(lrei,1);
01153       (*derindex)(num_separable_calls)(i)=i1;
01154     }
01155   }
01156 
01157   f1b2gradlist->reset();
01158   f1b2gradlist->list.initialize();
01159   f1b2gradlist->list2.initialize();
01160   f1b2gradlist->list3.initialize();
01161   f1b2gradlist->nlist.initialize();
01162   f1b2gradlist->nlist2.initialize();
01163   f1b2gradlist->nlist3.initialize();
01164   funnel_init_var::num_vars=0;
01165   funnel_init_var::num_active_parameters=0;
01166   funnel_init_var::num_inactive_vars=0;
01167 }
01168 
01173 imatrix laplace_approximation_calculator::check_sparse_matrix_structure(void)
01174 {
01175   ivector rowsize(1,usize);
01176   rowsize.initialize();
01177 
01178   imatrix tmp(1,usize,1,usize);
01179   tmp.initialize();
01180   for (int ii=1;ii<=num_separable_calls;ii++)
01181   {
01182     if (allocated((*derindex)(ii)))
01183     {
01184       ivector   iv = sort((*derindex)(ii));
01185       int n=iv.indexmax();
01186       if (n>1)                     // so we have off diagonal elements
01187       {
01188         for (int i=1;i<=n;i++)
01189         {
01190           int row=iv(i);
01191           for (int j=1;j<=n;j++)
01192           {
01193             if (i != j)
01194             {
01195               int col=iv(j);
01196               int  foundmatch=0;
01197               for (int k=1;k<=rowsize(row);k++)
01198               {
01199                 if (tmp(row,k)==col)
01200                 {
01201                   foundmatch=1;
01202                   break;
01203                 }
01204               }
01205               if (foundmatch==0)  // add a new element to tmp(row)
01206               {
01207                 rowsize(row)++;
01208                 if (rowsize(row)> tmp(row).indexmax())
01209                 {
01210                   tmp(row).reallocate(2);  // double the size
01211                 }
01212                 tmp(row,rowsize(row))=col;
01213               }
01214             }
01215           }
01216         }
01217       }
01218     }
01219   }
01220   imatrix M(1,usize,1,rowsize);
01221   ofstream ofs("sparseness.info");
01222   ofs << "# Number of parameters" << endl;
01223   ofs << usize << endl;
01224   ofs << "# Number of off diagonal elements in each row" << endl;
01225   ofs << rowsize << endl;
01226   ofs << "# The nonempty rows of M" << endl;
01227   for (int i=1;i<=usize;i++)
01228   {
01229     if (rowsize(i)>0)
01230     {
01231       M(i)=sort(tmp(i)(1,rowsize(i)));
01232       ofs << setw(4) << M(i) << endl;
01233     }
01234   }
01235   exit(1);
01236   return M;
01237 }