Changeset b85e9db in tol


Ignore:
Timestamp:
May 28, 2009, 2:53:03 PM (17 years ago)
Author:
Víctor de Buen <vdebuen@…>
Branches:
1941-docker-image-to-test, devel, master
Children:
7742d71
Parents:
8c67d66
Message:

New method get.evaluator to avoid more modular handling in draw.scalar and custom.partial.draw, that replaces old filter method marginalDraw

git-svn-id: ​https://www.tol-project.org/svn/tolp/trunk@1045 be4846ae-f051-483a-831f-b64e7d08527f

File:
1 edited

Legend:

Unmodified
Added
Removed

  • tol/stdlib/math/stat/models/bayesian/bysMcmc/bsr/gibbs/_nonLinBlock.tol

    r8c67d66 rb85e9db  
    1212"given BSR model.\n"
    1313"Arguments are:\n"
    14 "  NameBlock filterEval: Evaluates the non linear filter of inputs or "
     14"  NameBlock usrFilter: Evaluates the non linear filter of inputs or "
    1515   "output. Must have next mandatory methods:\n"
    1616"    Text get.name(Real unused): Identifies the filter\n"
    … …  
    2929  BlockSampler NonLinBlock(
    3030    //Filter evaluator
    31     NameBlock filterEval,
     31    NameBlock usrFilter,
    3232    //Basic BSR that will be modified by filter
    3333    NameBlock bsr,
    … …  
    4747  NameBlock _.bsr = bsr;
    4848  //Filter evaluator
    49   NameBlock _.filterEval = filterEval;
     49  NameBlock _.filter = usrFilter;
    5050 
    5151  //Defining the block in the Markov Chain
    5252  Real _defined = define
    5353  (
    54     _.filterEval::get.name(0),
    55     _.filterEval::get.colNames(0),
     54    _.filter::get.name(0),
     55    _.filter::get.colNames(0),
    5656    bsr::_.mcmcPartialHandler,
    5757    numBlock,
    … …  
    6060  );
    6161  //Name of affected segment in BSR
    62   Text _.segmentName = _.filterEval::get.segmentName(0);
     62  Text _.segmentName = _.filter::get.segmentName(0);
    6363  //Number of affected segment in BSR
    6464  Real _.segmentNum = FindIndexByName((_.bsr::_.modelDef)->NoiseDistrib,
    … …  
    7878  });
    7979  //For input filters is +1, for output ones -1
    80   Real _.isInputFilter = If(FunctionExist("Set","_.filterEval::get.linBlkNames"),
    81   {
    82     If(Card(_.filterEval::get.linBlkNames(0)),1,0)
     80  Real _.isInputFilter = If(FunctionExist("Set","_.filter::get.linBlkNames"),
     81  {
     82    If(Card(_.filter::get.linBlkNames(0)),1,0)
    8383  },
    8484  {
    8585    0
    8686  });
    87   //Default value for marginal draw
    88   Matrix _defaultMarginalDraw = Rand(Card(_.colNames), 1, ?, ?);
     87  //Default value for custom partial draw
     88  Matrix _default.custom.partial.draw = Rand(Card(_.colNames), 1, ?, ?);
    8989  //Optional filter method to be called before to call draw.scalar
    9090  //This method will be called just before to call draw method of corresponding
    9191  //internal non linear block in order to select which parameters will be drawn.
    92   //It will return a matrix with marginal sampling in selected indexes and ? in
     92  //It will return a matrix with custom partial sampling in selected indexes and ? in
    9393  //other ones that will be sampled by non linear block in the standard way
    94   Code _.marginalDraw = If(FunctionExist("Matrix","_.filterEval::marginalDraw"),
    95   {
    96     WriteLn(_MID+"Using specified marginalDraw method");
    97     _.filterEval::marginalDraw
     94  Code _.custom.partial.draw = If(
     95  FunctionExist("Matrix","_.filter::custom.partial.draw"),
     96  {
     97    WriteLn(_MID+"Using specified custom.partial.draw method");
     98    _.filter::custom.partial.draw
    9899  },
    99100  {
    100     WriteLn(_MID+"Using default marginalDraw method");
    101     Matrix marginalDraw(Real unused) { _defaultMarginalDraw }
     101    WriteLn(_MID+"Using default custom.partial.draw method");
     102    Matrix custom.partial.draw(Real unused) { _default.custom.partial.draw }
    102103  });
    103104  //Optional filter method that returns log of prior density for one parameter
    104105  //Real get.priorLogDens(Real paramIdx)
    105   Code _.get.priorLogDens = If(FunctionExist("Real","_.filterEval::get.priorLogDens"),
    106     _.filterEval::get.priorLogDens,
     106  Code _.get.priorLogDens = If(FunctionExist("Real","_.filter::get.priorLogDens"),
     107    _.filter::get.priorLogDens,
    107108    Real (Real paramIdx) { 0 } );
    108109  //Names of affected parameters in linear block
    … …  
    110111  Set _.linBlkNames = If(!_.isInputFilter, Copy(Empty),
    111112  {
    112     _.filterEval::get.linBlkNames(0)
     113    _.filter::get.linBlkNames(0)
    113114  });
    114115  //Positions of columns of affected parameters in linear block
    … …  
    131132  ////////////////////////////////////////////////////////////////////////////
    132133  {
    133     Real ok = _.filterEval::initialize(0);
    134     Real setStore(_.filterEval::get.parameter(0));
     134    Real ok = _.filter::initialize(0);
     135    Real setStore(_.filter::get.parameter(0));
    135136    ok
    136137  };
    137138 
     139  ///////////////////////////////////////////////////////////////////////////
     140    NameBlock get.evaluator(Real  unused)
     141  ///////////////////////////////////////////////////////////////////////////
     142  {
     143    Set segment = _.segment;
     144    Set aux     = _.segment->ArimaAuxInfo;
     145    Set arima   = _.segment->Arima;
     146    NameBlock sigBlk    = _.bsr::_.sig.blk;
     147    VMatrix si = SubRow(SubCol(sigBlk::_.sigmaInv,
     148                _.segment->EquIdx), _.segment->EquIdx);
     149
     150    NameBlock bsr = _.bsr;
     151    NameBlock filter = _.filter;
     152    Matrix u = VMat2Mat(_.store);
     153    Matrix omega= If(!_.isInputFilter, Col(1),
     154      SubRow(VMat2Mat(_.bsr::_.mln.blk::_.store), _.filterColIdx));
     155
     156    Matrix _.W := VMat2Mat(SubRow(bsr::_.noise,_.segment->EquIdx));
     157    Matrix _.Z := _.W + _.Fi * omega;
     158    Matrix Z = _.Z;
     159
     160    Code _getPriorLogDens = _.get.priorLogDens;
     161
     162    NameBlock [[
     163    Set       _.segment   = segment;
     164    VMatrix   _.fullNoise = bsr::_.noise;
     165    NameBlock _.arimaBlk  = bsr::_.arm.blk;
     166    NameBlock _.sigBlk    = bsr::_.sig.blk;
     167    Code _filterEval     = filter::eval;
     168    Code _get.priorLogDens = _getPriorLogDens;
     169    Matrix _.u = u;
     170    Matrix _.omega= omega;
     171    Matrix _.Z = Z;
     172
     173    VMatrix _.si = si;
     174    Polyn _.dif   = ARIMAGetDIF(arima);
     175    Real _.d = aux::difDegree;
     176    Real _.p = aux::p;
     177    Real _.q = aux::q;
     178    Real _.m = aux::m;
     179    Real _.n = _.m - _.d;
     180
     181    //Parameter index
     182    Real _.i = 1;
     183   
     184    Matrix EvalFilteredNoise(Real x)
     185    {
     186      Real old = MatDat(_.u,_.i,1);
     187      Real Do PutMatDat(_.u,_.i,1,x);
     188      Matrix Fi = _filterEval(_.u);
     189      Real Do PutMatDat(_.u,_.i,1,old);
     190      Matrix fn = VMat2Mat(_.si*Mat2VMat(_.Z-Fi*_.omega));
     191      fn
     192    };
     193    Real LogDens.WhiteNoise(Real x)
     194    {
     195      Matrix ei = EvalFilteredNoise(x);
     196      Real ld = -0.5*MatDat(MtMSqr(ei),1,1);
     197    //WriteLn("TRACE LogDens.WhiteNoise("+_.colNames[_.i]+","<<x+")="<<ld);
     198      ld+_get.priorLogDens(_.i)
     199    };
     200    Real LogDens.ARMA(Real x)
     201    {
     202      Matrix Zi = EvalFilteredNoise(x);
     203      VMatrix Zi_ = Mat2VMat(Zi);
     204      VMatrix ei_ = _.arimaBlk::filterJustBlock(_.segment, Zi_);
     205      Matrix  ei  = VMat2Mat(ei_);
     206      Real ld = -0.5*MatDat(MtMSqr(ei),1,1);
     207    //WriteLn("TRACE LogDens.ARMA("+_.colNames[_.i]+","<<x+")="<<ld);
     208      ld+_get.priorLogDens(_.i)
     209    };
     210    Real LogDens.ARIMA(Real x)
     211    {
     212      Matrix Zi = EvalFilteredNoise(x);
     213    //WriteLn("TRACE Zi("<<Rows(Zi)+","<<Columns(Zi)+")");
     214      Matrix Wi = Rand(_.d,1,0,0) << Sub(DifEq(_.dif/1,Zi),_.d+1,1,_.n,1);
     215    //WriteLn("TRACE Wi("<<Rows(Wi)+","<<Columns(Wi)+")");
     216      VMatrix Wi_ = Mat2VMat(Wi);
     217      VMatrix ei_ = _.arimaBlk::filterJustBlock(_.segment, Wi_);
     218    //WriteLn("TRACE ei_("<<VRows(ei_)+","<<VColumns(ei_)+")");
     219      Matrix  ei  = VMat2Mat(Sub(ei_,_.d+1,1,_.n,1));
     220    //WriteLn("TRACE ei("<<Rows(ei)+","<<Columns(ei)+")");
     221      Real ld = -0.5*MatDat(MtMSqr(ei),1,1);
     222    //WriteLn("TRACE LogDens.ARIMA("+_.colNames[_.i]+","<<x+")="<<ld);
     223      ld+_get.priorLogDens(_.i)
     224    };
     225    Code LogDens = Case(
     226      Or(Not(_.arimaBlk::enabled),And(Not(_.p),Not(_.q))),
     227        LogDens.WhiteNoise,
     228      _.d==0,                     
     229        LogDens.ARMA,
     230      1==1,                     
     231        LogDens.ARIMA )
     232  ]]};
     233
    138234  ////////////////////////////////////////////////////////////////////////////
    139235  Text _.autodoc.member.draw.scalar =
    … …  
    144240  {
    145241  //WriteLn("TRACE"+_MID+"draw.scalar 1 ");
    146     VMatrix   fullNoise = _.bsr::_.noise;
    147     NameBlock arimaBlk  = _.bsr::_.arm.blk;
    148     NameBlock sigBlk    = _.bsr::_.sig.blk;
    149     Matrix u = VMat2Mat(_.store);
    150   //WriteLn("TRACE"+_MID+"draw.scalar 1.1 _.isInputFilter = "<<_.isInputFilter);
    151     Matrix omega= If(!_.isInputFilter, Col(1),
    152       SubRow(VMat2Mat(_.bsr::_.mln.blk::_.store), _.filterColIdx));
    153   //WriteLn("TRACE"+_MID+"draw.scalar 1.2 omega("<<Rows(omega)+"x"<<Columns(omega));
    154   //WriteLn("TRACE"+_MID+"draw.scalar 1.3 _.Fi("<<Rows(_.Fi)+"x"<<Columns(_.Fi));
    155     Matrix _.W := VMat2Mat(SubRow(fullNoise,_.segment->EquIdx));
    156   //WriteLn("TRACE"+_MID+"draw.scalar 1.4 _.W("<<Rows(_.W)+"x"<<Columns(_.W));
    157     Matrix _.Z := _.W + _.Fi * omega;
    158   //WriteLn("TRACE"+_MID+"draw.scalar 1.5 _.Z("<<Rows(_.Z)+"x"<<Columns(_.Z));
    159     Set segment = _.segment;
    160     Set aux     = _.segment->ArimaAuxInfo;
    161     Set arima   = _.segment->Arima;
    162     Polyn dif   = ARIMAGetDIF(arima);
    163     //
    164     VMatrix si = SubRow(SubCol(sigBlk::_.sigmaInv,_.segment->EquIdx),_.segment->EquIdx);
    165     Real d = aux::difDegree;
    166     Real p = aux::p;
    167     Real q = aux::q;
    168     Real m = aux::m;
    169     Real n = m - d;
     242
    170243  //WriteLn("TRACE"+_MID+"draw.scalar 1.6 d="<<d+" p="<<p+" q="<<q+" m="<<m+" n="<<n);
    171244  //WriteLn("TRACE"+_MID+"draw.scalar 2");
    172     Matrix mrgDraw = _.marginalDraw(0);
     245   
     246    NameBlock ev = get.evaluator(0);
     247    Matrix mrgDraw = _.custom.partial.draw(0);
    173248    Set draw.all = For(1, Rows(mrgDraw), Real(Real i)
    174249    {
    175250      Real mrg = MatDat(mrgDraw,i,1);
    176251      If(Not(IsUnknown(mrg)), mrg, {
    177     //WriteLn("TRACE"+_MID+"draw.scalar.one 2.["+IntText(i)+"].1");
    178       Set bounds = _.filterEval::get.bounds(i,u);
    179       Code fltEv = _.filterEval::eval;
    180       Code get.priorLogDens = _.get.priorLogDens;
    181       Matrix Z = _.Z;
     252      Real ev::_.i := i;
     253    //WriteLn("TRACE"+_MID+"draw.scalar.one 2.["+IntText(i)+"].3");
     254      Real x = MatDat(ev::_.u,i,1);
     255    //WriteLn("TRACE"+_MID+"draw.scalar.one 2.["+IntText(i)+"].4");
     256      Set bounds = _.filter::get.bounds(i,ev::_.u);
    182257      Real x0 = bounds[1];
    183258      Real x1 = bounds[2];
    184     //WriteLn("TRACE"+_MID+"draw.scalar.one 2.["+IntText(i)+"].2");
    185     //WriteLn("TRACE"+_MID+"draw.scalar.one arimaBlk::enabled=["<<arimaBlk::enabled+"]");
    186     //WriteLn("TRACE"+_MID+"draw.scalar.one d=["<<d+"]");
    187     //WriteLn("TRACE"+_MID+"draw.scalar.one p=["<<p+"]");
    188     //WriteLn("TRACE"+_MID+"draw.scalar.one q=["<<q+"]");
    189     //WriteLn("TRACE"+_MID+"draw.scalar.one bounds=["<<x0+","<<x1+"]");
    190     //WriteLn("TRACE"+_MID+"draw.scalar.one omega="<<omega);
    191       Matrix EvalFilteredNoise(Real x)
    192       {
    193         Real old = MatDat(u,i,1);
    194         Real Do PutMatDat(u,i,1,x);
    195         Matrix Fi = fltEv(u);
    196         Real Do PutMatDat(u,i,1,old);
    197         Matrix fn = VMat2Mat(si*Mat2VMat(Z-Fi*omega));
    198         fn
    199       };
    200       Real LogDens.WhiteNoise(Real x)
    201       {
    202         Matrix ei = EvalFilteredNoise(x);
    203         Real ld = -0.5*MatDat(MtMSqr(ei),1,1);
    204       //WriteLn("TRACE LogDens.WhiteNoise("+_.colNames[i]+","<<x+")="<<ld);
    205         ld+get.priorLogDens(i)
    206       };
    207       Real LogDens.ARMA(Real x)
    208       {
    209         Matrix Zi = EvalFilteredNoise(x);
    210         VMatrix Zi_ = Mat2VMat(Zi);
    211         VMatrix ei_ = arimaBlk::filterJustBlock(segment, Zi_);
    212         Matrix  ei  = VMat2Mat(ei_);
    213         Real ld = -0.5*MatDat(MtMSqr(ei),1,1);
    214       //WriteLn("TRACE LogDens.ARMA("+_.colNames[i]+","<<x+")="<<ld);
    215         ld+get.priorLogDens(i)
    216       };
    217       Real LogDens.ARIMA(Real x)
    218       {
    219         Matrix Zi = EvalFilteredNoise(x);
    220       //WriteLn("TRACE Zi("<<Rows(Zi)+","<<Columns(Zi)+")");
    221         Matrix Wi = Rand(d,1,0,0) << Sub(DifEq(dif/1,Zi),d+1,1,n,1);
    222       //WriteLn("TRACE Wi("<<Rows(Wi)+","<<Columns(Wi)+")");
    223         VMatrix Wi_ = Mat2VMat(Wi);
    224         VMatrix ei_ = arimaBlk::filterJustBlock(segment, Wi_);
    225       //WriteLn("TRACE ei_("<<VRows(ei_)+","<<VColumns(ei_)+")");
    226         Matrix  ei  = VMat2Mat(Sub(ei_,d+1,1,n,1));
    227       //WriteLn("TRACE ei("<<Rows(ei)+","<<Columns(ei)+")");
    228         Real ld = -0.5*MatDat(MtMSqr(ei),1,1);
    229       //WriteLn("TRACE LogDens.ARIMA("+_.colNames[i]+","<<x+")="<<ld);
    230         ld+get.priorLogDens(i)
    231       };
    232       Code LogDens = Case(
    233         Or(Not(arimaBlk::enabled),And(Not(p),Not(q))),
    234           LogDens.WhiteNoise,
    235         d==0,                     
    236           LogDens.ARMA,
    237         1==1,                     
    238           LogDens.ARIMA );
    239     //WriteLn("TRACE"+_MID+"draw.scalar.one 2.["+IntText(i)+"].3");
    240       Real x = MatDat(u,i,1);
    241     //WriteLn("TRACE"+_MID+"draw.scalar.one 2.["+IntText(i)+"].4");
    242259      Real If(Or(x<x0,x>x1), x := Rand(x0,x1) );
    243260    //WriteLn("TRACE"+_MID+"draw.scalar.one 2.["+IntText(i)+"].5");
    244       Real x:= _scalarSampler(LogDens,x0,x1,x);
     261      Real x:= _scalarSampler(ev::LogDens,x0,x1,x);
    245262    //WriteLn("TRACE draw.scalar x= "<<x);
    246263    //WriteLn("TRACE"+_MID+"draw.scalar.one 2.["+IntText(i)+"].6");
    247       Real PutMatDat(u,i,1,x);
     264      Real PutMatDat(ev::_.u,i,1,x);
    248265    //WriteLn("TRACE"+_MID+"draw.scalar.one x="<<x);
    249266      x
    … …  
    297314    {
    298315    //WriteLn("TRACE "+_MID+"::setStore 3");
    299       Matrix _.Fi := _.filterEval::eval(VMat2Mat(_.store));
     316      Matrix _.Fi := _.filter::eval(VMat2Mat(_.store));
    300317    //WriteLn("TRACE "+_MID+"::setStore 4");
    301318      Real ok = If(!_.isInputFilter, 0,
Note: See TracChangeset for help on using the changeset viewer.