[Developers] fixing init_number_vector without a plan

[dave fournier]

At our workshop Mark pointed out some problem with the 
init_number_vector class for RE models.
I have identified some serious problems that are not so easy to fix so I 
wondered if, as an exercise
it would be easy to use some C++ capabilities to introduce a new class 
which works better,
without having to modify the existing code.  This sort of programming 
typically uses virtual
functions.

One starts with the  df1b2_init_vector class which is what gets invoked 
when you use an init_vector
in an RE model. Now the problem with a df1b2_init_vector is that it only
has an int for its starting phase.  Since we want the different elements
of the vector to have their own starting phase we will add that to the 
class.
Sow we derive a new   class xdf1b2_init_number_vector from
public df1b2_init_vector and add what we need to it.

class xdf1b2_init_number_vector : public df1b2_init_vector
{
     index_type * it;

    // ... other stuff

};

all that is needed is a pointer to index_type.  This will permit the use
of a vector or a number to specify the starting phases. If a number is used
they will all be the same.

Now what you can do is to use the class in a toy program.
Start with a tpl file like this.
********************************************************
********************************************************

DATA_SECTION
   int n
  !! n=2;
  ivector st(1,n)
  !! st(1)=1;
  !! st(2)=-1;


PARAMETER_SECTION
   init_number_vector beta(1,n,st)
   random_effects_vector eps(1,1)

   objective_function_value f

PROCEDURE_SECTION
   f+=0.5*norm2(eps);
   for (int i=1;i<=n;i++)
   {
     f+= square(beta(i)-1.0);
   }

GLOBALS_SECTION


#include <admodel.h>
   #include <df1b2fun.h>
   #include <adrndeff.h>
   class xdf1b2_init_number_vector : public df1b2_init_vector
   {
     index_type * it;
   public:
     xdf1b2_init_number_vector(void);

     void allocate(int imin,int imax,const index_type& ps,char * s)
     {
       df1b2_init_vector::allocate(imin,imax,s);
       it=new index_type(ps);
     }

   };

   xdf1b2_init_number_vector::xdf1b2_init_number_vector(void):
     df1b2_init_vector()
   {
     it=0;
   }

********************************************************
********************************************************
********************************************************
********************************************************

Run tpl2rem on this tpl file and change

   df1b2_init_number_vector  to  xdf1b2_init_number_vector

in the cpp and htp files and compile with debug and safe mode.
When you run the model you will ifnd that it runs out of an array
in
     void df1b2_init_vector set_value(const init_df1b2vector& _x,const 
int& _ii,
       const df1b2variable& pen);

The reason is that it is using the starting phase number from that class.
To change that behaviour we just need override the set_value function
with a new one in the derived class. It looks like this.



   void xdf1b2_init_number_vector::set_value(const init_df1b2vector& x,
     const int& _ii,const df1b2variable& pen)
   {
     ADUNCONST(int,ii)
     int mmin=indexmin();
     int mmax=indexmax();
     for (int i=mmin;i<=mmax;i++)
     {
       if (withinbound(0,ad_integer((*it)(i)),current_phase))
         (*this)(i) = (x(ii++));
     }
   }

The change is

   if (withinbound(0,ad_integer((*it)(i)),current_phase)

which has been added to the function.
The class declaration is extended to

   class xdf1b2_init_number_vector : public df1b2_init_vector
   {
     index_type * it;
   public:
     xdf1b2_init_number_vector(void);

     void allocate(int imin,int imax,const index_type& ps,char * s)
     {
       df1b2_init_vector::allocate(imin,imax,s);
       it=new index_type(ps);
     }
     virtual void set_value(const init_df1b2vector& _x,const int& _ii,
       const df1b2variable& pen);
   };

It seems to work, although there may be some other functions which will need
to be overridden as well.

What's my point?  Well no planning or committees were needed. All I had
to do was to run the code with the debugger and it told me what to do. 
Entire
process took about two hours.

Time for a drink.