[Developers] Control files and basic MCMC

[dave fournier]

On 14-01-10 12:00 PM, Arni Magnusson wrote:

I  added a bit of code to admb to do more or less what you want. Read 
the txt file.

     Dave

> Happy new year to you, core developers!
>
> As Johnoel mentions, there is a wishlist item #11 in the Reykjavik 
> Workshop Report that says:
>
> "Flag/switch class Simplify the use of control and MCMC files, perhaps 
> MFCL-like object capability"
>
> This wishlist item is really about two things:
>
> ---
>
> (1) Sometimes we want to make ADMB applications general and 
> user-friendly. By user-friendly, I don't mean beginner-friendly, but 
> an efficient interface for running a variety of related models without 
> recompiling. Fixing parameters and setting bounds, phases, and initial 
> values, and various model-specific settings. This is sometimes done 
> inside the main data file and sometimes in a separate control file.
>
> Take http://www.hafro.is/~arnima/pella/pella.tpl for example. The data 
> and parameter sections are rather convoluted, in order to provide a 
> PLUI (phase,lower,upper,init) interface for controlling the parameter 
> estimation.
>
> Of course the exact implementation will vary from model to model. But 
> writing the PLUI code feels unnecessarily clunky and repetitive. The 
> code for switching to the control file is rather ugly, and and so is 
> the code to set the parameter phases, bounds and inits. Perhaps these 
> tasks are so common and so useful that we can provide some new 
> functionality to make them easier for the ADMB programmer.
>
> This example is a tiny model with very few parameters. The importance 
> of this wishlist item becomes greater with dozens of parameters.
>
> ---
>
> (2) MCMC is an important feature of ADMB. I like telling potential and 
> new users that it's a piece of cake to run MCMC for any ADMB model: 
> you basically just type mymodel -mcmc.
>
> Well, sort of.
>
> The http://www.hafro.is/~arnima/pella/pella.tpl example shows a common 
> approach. First we have if(mceval_phase()) write_mcmc(); and later we 
> have FUNCTION write_mcmc, with a long workaround to put a header for 
> each column in the MCMC output.
>
> Wouldn't it be cool if the user could actually type mymodel -mcmc 
> without all that stuff, and get a text file with column headers, ready 
> to plot. For the most advanced users, the PSV file is just that, but 
> for the majority of users that's not a very intuitive or accessible 
> format. Again, a task that is so common and so useful could be made 
> easier to use.
>
> A related idea is an MCMC_SECTION with a stream object that's ready to 
> use, mcmc<<"blah"<<..., much like the report stream object in the 
> REPORT_SECTION.
>
> ---
>
> As you can tell, I'm interested in both of these wishlist items. They 
> may not be top priority, but potential areas where ADMB could be a bit 
> smoother around the edges. We can revisit them at the next workshop, 
> to decide whether and how some of this could be done. We did not look 
> at this at the Reykjavik workshop, but included it in the report as a 
> reminder of a discussion item after the workshop.
>
> All the best,
>
> Arni
>
>
>
> On Fri, 10 Jan 2014, Johnoel Ancheta wrote:
>
>> Hi all,
>>
>> Happy New Years.
>>
>> I looked over the workshop report and noticed a priority task 
>> "Flag/Switch class to simplify the use of control and MCMC files..."
>>
>> It says someone is working on this, but did not say who?
>>
>> Johnoel
>>
> _______________________________________________
> Developers mailing list
> Developers at admb-project.org
> http://lists.admb-project.org/mailman/listinfo/developers
>

-------------- next part --------------
I can show you how to change all the crap for
a variable to one line of code like

  init_bounded_number logr(logr_plui)

instead of

  !! phz = (int) logr_plui(1);
  !! lb  =       logr_plui(2);
  !! ub  =       logr_plui(3);
  init_bounded_number logr(lb,ub,phz)
It involves changin tpl2cpp to permit the new syntax
and adding a new allocate function in ADMB to implement it.

First, in tpl2cpp.lex search on init_bounded_number.

we find

<DEFINE_PARAMETERS>init_bounded_number {
        prior_found=1;
    BEGIN INIT_BOUNDED_NUMBER_DEF;
    fprintf(fdat,"%s","  param_init_bounded_number ");
                     }


the important thing here is 

    BEGIN INIT_BOUNDED_NUMBER_DEF;

That determines how the line containing  init_bounded_vector ...
will be parsed, so search for

   INIT_BOUNDED_NUMBER_DEF 

until we find


<INIT_BOUNDED_NUMBER_DEF>({name}\({float_num_exp},{float_num_exp},{num_exp}\)) { 

    before_part(tmp_string,yytext,'(');  // get x in x(1,4)
    fprintf(fdat,"%s",tmp_string);
    fprintf(fall,"  %s",tmp_string);
    fprintf(fall,"%s",".allocate");
    after_part(tmp_string1,yytext,'(');  // get x in x(1,4)
    before_part(tmp_string2,tmp_string1,')');
    fprintf(fall,"%s,\"%s\")",tmp_string2,tmp_string);
    fprintf(fdat,"%s",";\n");
    fprintf(fall,"%s",";\n");
    if (!params_defined)
    {
      BEGIN DEFINE_DATA;
    }
    else
    {
      if(prior_found) {
        if(prior_counter<MAX_PRIOR_CHECK) sprintf(prior_checker[prior_counter++],"%s",tmp_string);
        prior_found=0;
      }
      BEGIN DEFINE_PARAMETERS;
    }

                            }

Now there is nothing in this code that requires more than one arguement so we can add another case to it to process, one that looks like

   init_bounded_number xxx(yyy)

change the beginning of the above to

<INIT_BOUNDED_NUMBER_DEF>({name}\({float_num_exp},{float_num_exp},{num_exp}\))| 
<INIT_BOUNDED_NUMBER_DEF>({name}\({name}\)) {


This enables us to write

PARAMETER_SECTION
  // Estimated
  init_bounded_number logr(logr_plui)
  init_bounded_number logk(logk_plui)
  init_bounded_number loga(loga_plui)
  init_bounded_number logp(logp_plui)
  init_bounded_number logq(logq_plui)
  init_bounded_number logsigma(logsigma_plui)

Of course now we need to modify the admb source to add an allocate function 
which takes a dvector. We need to find the type of logr. In pella.htp we
find
 
  param_init_bounded_number logr;

So we need the new allocate function for the param_init_bounded_number class.
find the class declaration by a little search through the header files
and we get
// need forward declaration of data_vector
class data_vector;

class param_init_bounded_number: public param_init_number
{
public:
  double get_minb(void);
  void set_minb(double b);
  double get_maxb(void);
  void set_maxb(double b);
protected:
  double minb;
  double maxb;
  void allocate(const data_vector& v,const char * s="UNNAMED");
  void allocate(double _minb,double _maxb,int phase_start=1,
    const char * s="UNNAMED");
  void allocate(double _minb,double _maxb,const char * s="UNNAMED");
  // ...

where I added

  void allocate(const data_vector& v,const char * s="UNNAMED");

and the forward declaration.

  class data_vector;

Finally we need to write the code for the allocate function.
This goes into model.cpp because the other allocate member functions are there.

  void param_init_bounded_number::allocate(const data_vector & v,
    const char * _s)
  {
   int phz=int(v(1));
   double lb=v(2);
   double ub=v(3);
   allocate(lb,ub,phz);
  }


-------------- next part --------------
A non-text attachment was scrubbed...
Name: pella_stuff.zip
Type: application/zip
Size: 21349 bytes
Desc: not available
URL: <http://lists.admb-project.org/pipermail/developers/attachments/20140113/c8d341b6/attachment-0001.zip>