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