ADMB Documentation  11.1x.2711
 All Classes Files Functions Variables Typedefs Friends Defines
shared.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: shared.cpp 2587 2014-11-07 19:38:15Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00007 #include <admodel.h>
00008 
00009 #  if defined(USE_SHARE_FLAGS)
00010   static int integer(const index_type& it)
00011   {
00012     return it.integer();
00013   }
00014 
00015   int initial_params::shared_size_count(void)
00016   {
00017     cerr << " Error -- shared_size_count not defined for this class"
00018          << endl;
00019     ad_exit(1);
00020     return 0;
00021   }
00022 
00023   void initial_params::shared_set_value(const dvar_vector& x,
00024     const int& ii,const dvariable& pen)
00025   {
00026     cerr << " Error -- shared_scaled_set_value_inv not defined for this class"
00027          << endl;
00028     ad_exit(1);
00029   }
00030 
00031   void initial_params::shared_set_value_inv(const dvector& x,const int& ii)
00032   {
00033     cerr << " Error -- shared_scaled_set_value_inv not defined for this class"
00034          << endl;
00035     ad_exit(1);
00036   }
00037   int param_init_matrix::shared_size_count(void)
00038   {
00039     if (share_flags->get_current_phase() != current_phase)
00040     {
00041       share_flags->get_inv_matrix_shared(current_phase);
00042     }
00043     return share_flags->get_maxshare();
00044   }
00045 
00046   int param_init_vector::shared_size_count(void)
00047   {
00048     if (share_flags->get_current_phase() != current_phase)
00049     {
00050       share_flags->get_inv_vector_shared(current_phase);
00051     }
00052     return share_flags->get_maxshare();
00053   }
00054 
00055   void param_init_matrix::shared_set_value_inv(const dvector& _x,const int& _ii)
00056   {
00057     ADUNCONST(int,ii)
00058     ADUNCONST(dvector,x)
00059     index_type& it=*(share_flags->get_invflags());
00060     //int mmin=share_flags->invflags->indexmin();
00061     //int mmax=share_flags->invflags->indexmax();
00062     int mmin=it.indexmin();
00063     int mmax=it.indexmax();
00064     for (int i=mmin;i<=mmax;i++)
00065     {
00066       x(ii++)=value((*this)(it(i)(1).integer(),it(i)(2).integer()));
00067     }
00068   }
00069 
00070   void param_init_vector::shared_set_value_inv(const dvector& _x,const int& _ii)
00071   {
00072     ADUNCONST(int,ii)
00073     ADUNCONST(dvector,x)
00074     index_type& it=*(share_flags->get_invflags());
00075     int mmin=it.indexmin();
00076     int mmax=it.indexmax();
00077     for (int i=mmin;i<=mmax;i++)
00078     {
00079       x(ii++)=value((*this)(it(i).integer()));
00080     }
00081   }
00082 
00083   void param_init_bounded_matrix::shared_set_value_inv
00084     (const dvector& _x,const int& _ii)
00085   {
00086     ADUNCONST(int,ii)
00087     ADUNCONST(dvector,x)
00088     index_type& it=*(share_flags->get_invflags());
00089     int mmin=it.indexmin();
00090     int mmax=it.indexmax();
00091     for (int i=mmin;i<=mmax;i++)
00092     {
00093       x(ii++)=
00094         boundpin((*this)(it(i)(1).integer(),it(i)(2).integer()),
00095         minb,maxb);
00096     }
00097   }
00098 
00099   void param_init_matrix::shared_set_value(const dvar_vector& _x,
00100     const int& _ii,const dvariable& pen)
00101   {
00102     ADUNCONST(int,ii)
00103     ADUNCONST(dvar_vector,x)
00104     int mmin=indexmin();
00105     int mmax=indexmax();
00106     index_type& sf=*(share_flags->get_shareflags());
00107     index_type& af=*(share_flags->get_activeflags());
00108     for (int i=mmin;i<=mmax;i++)
00109     {
00110       int jmin=(*this)(i).indexmin();
00111       int jmax=(*this)(i).indexmax();
00112       for (int j=jmin;j<=jmax;j++)
00113       {
00114         int afvalue=integer(af(i)(j));
00115         if (afvalue && afvalue<=current_phase)
00116         {
00117           int is=sf(i)(j).integer();
00118           (*this)(i,j)=x(ii+is-1);
00119         }
00120       }
00121     }
00122     ii+=share_flags->get_maxshare();
00123   }
00124 
00125   void param_init_vector::shared_set_value(const dvar_vector& _x,
00126     const int& _ii,const dvariable& pen)
00127   {
00128     ADUNCONST(int,ii)
00129     ADUNCONST(dvar_vector,x)
00130     int mmin=indexmin();
00131     int mmax=indexmax();
00132     index_type& sf=*(share_flags->get_shareflags());
00133     index_type& af=*(share_flags->get_activeflags());
00134     for (int i=mmin;i<=mmax;i++)
00135     {
00136       int afvalue=integer(af(i));
00137       if (afvalue && afvalue<=current_phase)
00138       {
00139         int is=sf(i).integer();
00140         (*this)(i)=x(ii+is-1);
00141       }
00142     }
00143     ii+=share_flags->get_maxshare();
00144   }
00145 
00146   void param_init_bounded_matrix::shared_set_value(const dvar_vector& _x,
00147     const int& _ii,const dvariable& pen)
00148   {
00149     ADUNCONST(int,ii)
00150     ADUNCONST(dvar_vector,x)
00151     int mmin=indexmin();
00152     int mmax=indexmax();
00153     index_type& sf=*(share_flags->get_shareflags());
00154     index_type& af=*(share_flags->get_activeflags());
00155     for (int i=mmin;i<=mmax;i++)
00156     {
00157       int jmin=(*this)(i).indexmin();
00158       int jmax=(*this)(i).indexmax();
00159       for (int j=jmin;j<=jmax;j++)
00160       {
00161         int afvalue=integer(af(i)(j));
00162         if (afvalue && afvalue<=current_phase)
00163         {
00164           int is=sf(i)(j).integer();
00165           //(*this)(i,j)=x(ii+is-1);
00166           (*this)(i,j)=boundp(x(ii+is-1),minb,maxb,pen);
00167         }
00168       }
00169     }
00170     ii+=share_flags->get_maxshare();
00171   }
00172 
00173   void initial_params::setshare(const index_type& sf,
00174     const index_type& af)
00175   {
00176     cerr << " setshare not yet defined for this class " << endl;
00177     ad_exit(1);
00178   }
00179 
00180   shareinfo::~shareinfo(void)
00181   {
00182     delete sflags;
00183     delete original_sflags;
00184     delete aflags;
00185     delete invflags;
00186     delete bmap;
00187     sflags=0;
00188     aflags=0;
00189     original_sflags=0;
00190     invflags=0;
00191     current_phase=-1;
00192     maxshare=0;
00193     bmap=0;
00194   }
00195 
00196   int &  shareinfo::get_current_phase(void)
00197   {
00198     return current_phase;
00199   }
00200   index_type * shareinfo::get_original_shareflags(void)
00201   {
00202     return original_sflags;
00203   }
00204   index_type * shareinfo::get_shareflags(void)
00205   {
00206     return sflags;
00207   }
00208   i3_array & shareinfo::get_bmap(void)
00209   {
00210     return *bmap;
00211   }
00212   index_type * shareinfo::get_invflags(void)
00213   {
00214     return invflags;
00215   }
00216   index_type * shareinfo::get_activeflags(void)
00217   {
00218     return aflags;
00219   }
00220   void shareinfo::set_invflags(const index_type& sf)
00221   {
00222     invflags=new index_type(sf);
00223   }
00224   void shareinfo::set_bmap(const i3_array & _bmap)
00225   {
00226     bmap=new i3_array(_bmap);
00227   }
00228   void shareinfo::reset_bmap(const i3_array & _bmap)
00229   {
00230     if (bmap)
00231     {
00232       delete bmap;
00233       bmap=0;
00234     }
00235     bmap=new i3_array(_bmap);
00236   }
00237   void shareinfo::reset_shareflags(const index_type& sf)
00238   {
00239     if (sflags)
00240     {
00241       delete sflags;
00242       sflags=0;
00243     }
00244     sflags=new index_type(sf);
00245   }
00246   void shareinfo::set_original_shareflags(const index_type& sf)
00247   {
00248     original_sflags=new index_type(sf);
00249   }
00250   void shareinfo::set_shareflags(const index_type& sf)
00251   {
00252     sflags=new index_type(sf);
00253   }
00254   void shareinfo::set_activeflags(const index_type& af)
00255   {
00256     aflags=new index_type(af);
00257   }
00258   shareinfo::shareinfo(const index_type& sf,const index_type& af)
00259   {
00260     set_shareflags(sf);
00261     set_original_shareflags(sf);
00262     set_activeflags(af);
00263     invflags=0;
00264     current_phase=-1;
00265   }
00266 
00267   void shareinfo::get_inv_matrix_shared(int _cf)
00268   {
00269     if (current_phase != _cf)
00270     {
00271       current_phase = _cf;
00272       const index_type& sf=*(get_original_shareflags());
00273       const index_type& af=*(get_activeflags());
00274       if (sf.dimension() !=2)
00275       {
00276         cerr << "error " << endl;
00277         ad_exit(1);
00278       }
00279       int imin= sf.indexmin();
00280       int imax= sf.indexmax();
00281       int fmin = 0, fmax = 0;
00282       int ibreak=0;
00283       // get intial values for min and max active flag values
00284       for (int i=imin;i<=imax;i++)
00285       {
00286         int jmin= sf(i).indexmin();
00287         int jmax= sf(i).indexmax();
00288         for (int j=jmin;j<=jmax;j++)
00289         {
00290           int afvalue=integer(af(i)(j));
00291           if (afvalue && afvalue<=current_phase)
00292           {
00293             fmax=integer(sf(i)(j));
00294             fmin=integer(sf(i)(j));
00295             ibreak=1;
00296             break;
00297           }
00298         }
00299         if (ibreak)
00300           break;
00301       }
00302       // get initial values for minimum and maximum shared
00303       // flags -- not just active ones
00304       int fmax1=integer(sf(imin)(sf(imin).indexmin()));
00305       int fmin1=integer(sf(imin)(sf(imin).indexmin()));
00306       for (int i=imin;i<=imax;i++)
00307       {
00308         int jmin= sf(i).indexmin();
00309         int jmax= sf(i).indexmax();
00310         for (int j=jmin;j<=jmax;j++)
00311         {
00312           fmax1=max(fmax1,integer(sf(i)(j)));
00313           fmin1=min(fmin1,integer(sf(i)(j)));
00314           int afvalue=integer(af(i)(j));
00315           if (afvalue && afvalue<=current_phase)
00316           {
00317             fmax=max(fmax,integer(sf(i)(j)));
00318             fmin=min(fmin,integer(sf(i)(j)));
00319           }
00320         }
00321       }
00322       // set up info for sanity test on shared and active flags
00323       ivector icount2(fmin1,fmax1);
00324       icount2.initialize();
00325       for (int i=imin;i<=imax;i++)
00326       {
00327         int jmin= sf(i).indexmin();
00328         int jmax= sf(i).indexmax();
00329         for (int j=jmin;j<=jmax;j++)
00330         {
00331           int sfvalue=integer(sf(i)(j));
00332           icount2(sfvalue)++;
00333         }
00334       }
00335       i3_array bmap2(fmin1,fmax1,1,icount2,1,2);
00336       icount2.initialize();
00337       for (int i=imin;i<=imax;i++)
00338       {
00339         int jmin= sf(i).indexmin();
00340         int jmax= sf(i).indexmax();
00341         for (int j=jmin;j<=jmax;j++)
00342         {
00343           int sfvalue=integer(sf(i)(j));
00344           int ind=++icount2(sfvalue);
00345           bmap2(sfvalue,ind,1)=i;
00346           bmap2(sfvalue,ind,2)=j;
00347         }
00348       }
00349       for (int k=fmin1;k<=fmax1;k++)
00350       {
00351         int lmin=bmap2(k).indexmin();
00352         int lmax=bmap2(k).indexmax();
00353 
00354         if (lmax>0)
00355         {
00356           int i1=bmap2(k,lmin,1);
00357           int j1=bmap2(k,lmin,2);
00358           int a1=integer(af(i1)(j1));
00359 
00360           for (int l=lmin+1;l<=lmax;l++)
00361           {
00362             int i2=bmap2(k,l,1);
00363             int j2=bmap2(k,l,2);
00364             int a2=integer(af(i2)(j2));
00365             if (a1 !=a2)
00366             {
00367               cerr << "Sanity error in grouping/active flags"
00368                << endl << "for indices (" << i2 << "," << j2 << ")"
00369                << endl;
00370               ad_exit(1);
00371             }
00372           }
00373         }
00374       }
00375 
00376       // indirect will cotain pointers for the number
00377       // of active parameters it starts out with the
00378       // number of shared flags and removes the inactive ones
00379       ivector indirect(fmin1,fmax1);
00380       ivector processed_flag(fmin1,fmax1);
00381       processed_flag.initialize();
00382       ivector mindx(imin,imax);
00383       ivector maxx(imin,imax);
00384       indirect.fill_seqadd(1,1);
00385       for (int i=imin;i<=imax;i++)
00386       {
00387         int jmin= sf(i).indexmin();
00388         int jmax= sf(i).indexmax();
00389         mindx(i)=jmin;
00390         maxx(i)=jmax;
00391         for (int j=jmin;j<=jmax;j++)
00392         {
00393           int afvalue=integer(af(i)(j));
00394           if (afvalue==0 || afvalue>current_phase)
00395           {
00396             int in=integer(sf(i)(j));
00397             if (processed_flag(in)==0)
00398             {
00399               processed_flag(in)=1;
00400               for (int k=in;k<=fmax1;k++)
00401               {
00402                 indirect(k)-=1;
00403               }
00404             }
00405           }
00406         }
00407       }
00408       imatrix tmp1(imin,imax,mindx,maxx);
00409       for (int i=imin;i<=imax;i++)
00410       {
00411         int jmin= sf(i).indexmin();
00412         int jmax= sf(i).indexmax();
00413         for (int j=jmin;j<=jmax;j++)
00414         {
00415           int afvalue=integer(af(i)(j));
00416           if (afvalue && afvalue<=current_phase)
00417           {
00418              tmp1(i,j)=indirect(integer(sf(i)(j)));
00419           }
00420           else
00421           {
00422             tmp1(i,j)=0;
00423           }
00424         }
00425       }
00426       int itmp=indirect(fmax1);
00427       imatrix tmp(1,itmp,1,2);
00428       ivector counter(1,itmp);
00429       counter.initialize();
00430       for (int i=imin;i<=imax;i++)
00431       {
00432         int jmin= sf(i).indexmin();
00433         int jmax= sf(i).indexmax();
00434         for (int j=jmin;j<=jmax;j++)
00435         {
00436           int afvalue=integer(af(i)(j));
00437           if (afvalue && afvalue<=current_phase)
00438           {
00439             if (++counter(tmp1(i,j))==1)
00440             {
00441               tmp(tmp1(i,j),1)=i;
00442               tmp(tmp1(i,j),2)=j;
00443             }
00444           }
00445         }
00446       }
00447       i3_array _bmap(1,itmp,1,counter,1,2);
00448 
00449       counter.initialize();
00450       _bmap.initialize();
00451       for (int i=imin;i<=imax;i++)
00452       {
00453         int jmin= sf(i).indexmin();
00454         int jmax= sf(i).indexmax();
00455         for (int j=jmin;j<=jmax;j++)
00456         {
00457           int afvalue=integer(af(i)(j));
00458           if (afvalue && afvalue<=current_phase)
00459           {
00460             int ind=++counter(tmp1(i,j));
00461             _bmap(tmp1(i,j),ind,1)=i;
00462             _bmap(tmp1(i,j),ind,2)=j;
00463           }
00464         }
00465       }
00466       {
00467         cout <<  endl;
00468         for (int i=1;i<=itmp;i++)
00469           cout << _bmap(i) << endl << endl;
00470       }
00471       set_bmap(_bmap);
00472 
00473       get_maxshare()=itmp;
00474       reset_shareflags(tmp1);
00475       set_invflags(tmp);
00476       cout << tmp1 << endl;
00477       //return tmp1;
00478     }
00479   }
00480 
00481   void shareinfo::get_inv_vector_shared(int _cf)
00482   {
00483     if (current_phase != _cf)
00484     {
00485       current_phase = _cf;
00486       const index_type& sf=*(get_original_shareflags());
00487       const index_type& af=*(get_activeflags());
00488       if (sf.dimension() !=1)
00489       {
00490         cerr << "error " << endl;
00491         ad_exit(1);
00492       }
00493       int imin= sf.indexmin();
00494       int imax= sf.indexmax();
00495       //int ibreak=0;
00496       // get intial values for min and max active flag values
00497       // get initial values for minimum and maximum shared
00498       // flags -- not just active ones
00499       int fmax1=integer(sf(imin));
00500       int fmin1=integer(sf(imin));
00501       for (int i=imin;i<=imax;i++)
00502       {
00503         fmax1=max(fmax1,integer(sf(i)));
00504         fmin1=min(fmin1,integer(sf(i)));
00505       }
00506 
00507       ivector indirect(fmin1,fmax1);
00508       ivector processed_flag(fmin1,fmax1);
00509       processed_flag.initialize();
00510       indirect.fill_seqadd(1,1);
00511       for (int i=imin;i<=imax;i++)
00512       {
00513         {
00514           int afvalue=integer(af(i));
00515           if (afvalue==0 || afvalue>current_phase)
00516           {
00517             int in=integer(sf(i));
00518             if (processed_flag(in)==0)
00519             {
00520               processed_flag(in)=1;
00521               for (int k=in;k<=fmax1;k++)
00522               {
00523                 indirect(k)-=1;
00524               }
00525             }
00526           }
00527         }
00528       }
00529       ivector tmp1(imin,imax);
00530       for (int i=imin;i<=imax;i++)
00531       {
00532         {
00533           int afvalue=integer(af(i));
00534           if (afvalue && afvalue<=current_phase)
00535           {
00536              tmp1(i)=indirect(integer(sf(i)));
00537           }
00538           else
00539           {
00540             tmp1(i)=0;
00541           }
00542         }
00543       }
00544       int itmp=indirect(fmax1);
00545       ivector tmp(1,itmp);
00546       ivector counter(1,itmp);
00547       counter.initialize();
00548       for (int i=imin;i<=imax;i++)
00549       {
00550         int afvalue=integer(af(i));
00551         if (afvalue && afvalue<=current_phase)
00552         {
00553           if (++counter(tmp1(i))==1)
00554           {
00555             tmp(tmp1(i))=i;
00556           }
00557         }
00558       }
00559       get_maxshare()=itmp;
00560       reset_shareflags(tmp1);
00561       set_invflags(tmp);
00562     }
00563   }
00564 
00565   void param_init_vector::setshare(const index_type& sf,const index_type& af)
00566   {
00567     share_flags = new shareinfo(sf,af);
00568     int idim1= share_flags->get_shareflags()->dimension();
00569     share_flags->get_dimension()=idim1;
00570     //int idim2= share_flags->get_activeflags()->dimension();
00571     if (idim1==2)
00572     {
00573       cout << " Error dim too high" << endl;
00574       ad_exit(1);
00575     }
00576        // check rationality
00577     /*
00578       int mmin1= share_flags->get_shareflags()->indexmin();
00579       int mmax1= share_flags->get_shareflags()->indexmax();
00580       int mmin2= share_flags->get_activeflags()->indexmin();
00581       int mmax2= share_flags->get_activeflags()->indexmax();
00582       int mmin=indexmin();
00583       int mmax=indexmax();
00584       if (mmin1 != mmin || mmax1 != mmax ||
00585         mmin2 != mmin || mmax2 != mmax)
00586       {
00587         cerr << "sanity error" << endl;
00588         ad_exit(1);
00589       }
00590     */
00591     share_flags->get_inv_vector_shared(current_phase);
00592   }
00593 
00594   void param_init_matrix::setshare(const index_type& sf,const index_type& af)
00595   {
00596     share_flags = new shareinfo(sf,af);
00597     int idim1= share_flags->get_shareflags()->dimension();
00598     share_flags->get_dimension()=idim1;
00599     int idim2= share_flags->get_activeflags()->dimension();
00600     switch (idim1)
00601     {
00602     case 1:
00603       share_flags->get_inv_vector_shared(current_phase);
00604       break;
00605     case 2:
00606       {
00607         cout << idim1 << " " << idim2 << endl;
00608          // check rationality
00609         int mmin1= share_flags->get_shareflags()->indexmin();
00610         int mmax1= share_flags->get_shareflags()->indexmax();
00611         int mmin2= share_flags->get_activeflags()->indexmin();
00612         int mmax2= share_flags->get_activeflags()->indexmax();
00613         int mmin=indexmin();
00614         int mmax=indexmax();
00615         if (mmin1 != mmin || mmax1 != mmax ||
00616           mmin2 != mmin || mmax2 != mmax)
00617         {
00618           cerr << "sanity error" << endl;
00619           ad_exit(1);
00620         }
00621         share_flags->get_inv_matrix_shared(current_phase);
00622       }
00623       break;
00624     default:
00625       cerr << "Error" << endl;
00626       ad_exit(1);
00627     }
00628   }
00629 #  endif