Changeset 4db909c in tol for tol_tests


Ignore:
Timestamp:
Apr 9, 2010, 9:52:03 AM (17 years ago)
Author:
Víctor de Buen <vdebuen@…>
Branches:
1941-docker-image-to-test, devel, master
Children:
c7cc179
Parents:
aa7db24
Message:

división de draw_multiple_try en sub-modulos para ver detalles del profiler
Aceleración de arima2col y col2arima

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

File:
1 edited

Legend:

Unmodified
Added
Removed

  • tol_tests/tol/ARIMA/mcmc/mtrymh.tol

    raa7db24 r4db909c  
    179179    If(_.k==1, draw_single_try(void), {
    180180    Real _.num_draw := _.num_draw + 1;
     181    Real lwYjx_max = ?;
     182    Real lwXjy_max = ?;
     183    Real gYjx_sum  = ?;
     184    Real gXjy_sum = ?;
    181185    VMatrix J = DifEq(1/(1-B),Constant(_.k,1,1));
    182     Set Y = For(1,_.k, VMatrix (Real j) { Q_draw(_.x) });
    183     VMatrix lwYjx = Mat2VMat(SetCol(For(1,_.k, Real(Real j)
    184     {
    185       Real _.num_logWeight := _.num_logWeight+1;
    186       logWeight(Y[j],_.x)
    187     })));
    188   //WriteLn("TRACE draw lwYjx=\n"<<Matrix Tra(VMat2Mat(lwYjx)));
    189     Real lwYjx_max = VMatMax(lwYjx);
    190   //WriteLn("TRACE draw lwYjx_max=\n"<<lwYjx_max);
    191     VMatrix gYjx = Exp(lwYjx-lwYjx_max);
    192   //WriteLn("TRACE draw gYjx=\n"<<Matrix Tra(VMat2Mat(gYjx)));
    193     Real gYjx_sum = VMatSum(gYjx);
    194   //WriteLn("TRACE draw gYjx_sum=\n"<<gYjx_sum);
    195     Matrix fYjx = VMat2Mat(gYjx / gYjx_sum);
    196   //WriteLn("TRACE draw fYjx=\n"<<Matrix Tra(fYjx));
    197     Matrix perm = Sort(VMat2Mat(lwYjx),[[1]]);
    198   //WriteLn("TRACE draw perm=\n"<<Matrix Tra(perm));
    199     Matrix fxYp = PivotByRows(fYjx,perm);
    200   //WriteLn("TRACE draw fxYp=\n"<<Matrix Tra(fxYp));
    201     Matrix FxYp = DifEq(1/(1-B),fxYp);
    202   //WriteLn("TRACE draw FxYp=\n"<<Matrix Tra(FxYp));
    203     Real u = Rand(0,1);
    204   //WriteLn("TRACE draw Constant(_.k,1,u)=\n"<<Matrix Tra(Constant(_.k,1,u)));
    205     Matrix FLTu = LT(FxYp,Constant(_.k,1,u));
    206   //WriteLn("TRACE draw FLTu=\n"<<Matrix Tra(FLTu));
    207     Real sFLTu = MatSum(FLTu);
    208     Real h = Min(_.k,1+sFLTu);
    209     Real j = MatDat(perm,h,1);
    210   //WriteLn("TRACE draw u="<<u+" sFLTu="<<sFLTu+" h="<<h+" j="<<j);
    211     VMatrix y = Y[j];
    212   //WriteLn("TRACE draw pi_matchesDomain(y)="<<pi_matchesDomain(y));
    213     Set X = For(1,_.k-1, VMatrix (Real j) { Q_draw(y) }) << [[_.x]];
    214     VMatrix lwXjy = Mat2VMat(SetCol(For(1,_.k, Real(Real j)
    215     {
    216       Real _.num_logWeight := _.num_logWeight+1;
    217       logWeight(X[j],y)
    218     })));
    219   //WriteLn("TRACE draw lwXjy=\n"<<Matrix Tra(VMat2Mat(lwXjy)));
    220     Real lwXjy_max = VMatMax(lwXjy);
    221   //WriteLn("TRACE draw lwXjy_max=\n"<<lwXjy_max);
    222     VMatrix gXjy = Exp(lwXjy-lwXjy_max);
    223     Real gXjy_sum = VMatSum(gXjy);
    224   //WriteLn("TRACE draw gXjy_sum=\n"<<gXjy_sum);
     186    Set Y = For(1,_.k, VMatrix Y_Q_draw(Real j) { Q_draw(_.x) });
     187   
     188    Real calc_numerator_Yjx(Real void)
     189    {
     190      VMatrix lwYjx = Mat2VMat(SetCol(For(1,_.k, Real(Real j)
     191      {
     192        Real _.num_logWeight := _.num_logWeight+1;
     193        logWeight(Y[j],_.x)
     194      })));
     195    //WriteLn("TRACE draw lwYjx=\n"<<Matrix Tra(VMat2Mat(lwYjx)));
     196      Real lwYjx_max := VMatMax(lwYjx);
     197    //WriteLn("TRACE draw lwYjx_max=\n"<<lwYjx_max);
     198      VMatrix gYjx = Exp(lwYjx-lwYjx_max);
     199    //WriteLn("TRACE draw gYjx=\n"<<Matrix Tra(VMat2Mat(gYjx)));
     200      Real gYjx_sum := VMatSum(gYjx);
     201    //WriteLn("TRACE draw gYjx_sum=\n"<<gYjx_sum);
     202      True
     203    };
     204    Real calc_numerator_Yjx(0);
     205
     206    VMatrix draw_y(Real void)
     207    {
     208      Matrix fYjx = VMat2Mat(gYjx / gYjx_sum);
     209    //WriteLn("TRACE draw fYjx=\n"<<Matrix Tra(fYjx));
     210      Matrix perm = Sort(VMat2Mat(lwYjx),[[1]]);
     211    //WriteLn("TRACE draw perm=\n"<<Matrix Tra(perm));
     212      Matrix fxYp = PivotByRows(fYjx,perm);
     213    //WriteLn("TRACE draw fxYp=\n"<<Matrix Tra(fxYp));
     214      Matrix FxYp = DifEq(1/(1-B),fxYp);
     215    //WriteLn("TRACE draw FxYp=\n"<<Matrix Tra(FxYp));
     216      Real u = Rand(0,1);
     217    //WriteLn("TRACE draw Constant(_.k,1,u)=\n"<<Matrix Tra(Constant(_.k,1,u)));
     218      Matrix FLTu = LT(FxYp,Constant(_.k,1,u));
     219    //WriteLn("TRACE draw FLTu=\n"<<Matrix Tra(FLTu));
     220      Real sFLTu = MatSum(FLTu);
     221      Real h = Min(_.k,1+sFLTu);
     222      Real j = MatDat(perm,h,1);
     223    //WriteLn("TRACE draw u="<<u+" sFLTu="<<sFLTu+" h="<<h+" j="<<j);
     224      Y[j]
     225    };
     226    VMatrix y = draw_y(0);
     227
     228    Set X = For(1,_.k-1, VMatrix X_Q_draw(Real j) { Q_draw(y) }) << [[_.x]];
     229
     230    Real calc_denominator_Xjy(Real void)
     231    {
     232    //WriteLn("TRACE draw pi_matchesDomain(y)="<<pi_matchesDomain(y));
     233      VMatrix lwXjy = Mat2VMat(SetCol(For(1,_.k, Real(Real j)
     234      {
     235        Real _.num_logWeight := _.num_logWeight+1;
     236        logWeight(X[j],y)
     237      })));
     238    //WriteLn("TRACE draw lwXjy=\n"<<Matrix Tra(VMat2Mat(lwXjy)));
     239      Real lwXjy_max := VMatMax(lwXjy);
     240    //WriteLn("TRACE draw lwXjy_max=\n"<<lwXjy_max);
     241      VMatrix gXjy = Exp(lwXjy-lwXjy_max);
     242      Real gXjy_sum := VMatSum(gXjy);
     243      True
     244    };
     245    Real calc_denominator_Xjy(0);
     246
     247    //WriteLn("TRACE draw gXjy_sum=\n"<<gXjy_sum);
    225248    Real _.r := Min(1,Exp(lwYjx_max-lwXjy_max)*(gYjx_sum/gXjy_sum));
    226249  //WriteLn("TRACE draw _.r=\n"<<_.r);
    … …  
    444467  Real _.sigma;
    445468  Set _.arima = Copy(Empty);
     469  Set _.arma = Copy(Empty);
    446470  Polyn  _.ar = 1;
    447471  Polyn  _.ma = 1;
    … …  
    450474  VMatrix _.z0 = Rand(0,0,0,0);
    451475  VMatrix _.a0 = Rand(0,0,0,0);
    452 
    453476 
    454477  ////////////////////////////////////////////////////////////////////////////
    … …  
    456479  ////////////////////////////////////////////////////////////////////////////
    457480  {
    458     Mat2VMat(SetCol(ARIMAToParam(arima)))
     481  //Mat2VMat(SetCol(ARIMAToParam(arima)))
     482    Set arma = EvalSet(arima, Set(ARIMAStruct f)
     483    {
     484      SetOfPolyn(1-f->AR, 1-f->MA)
     485    });
     486    Mat2VMat(Tra(GetNumeric(arma)))
    459487  };
    460488 
    … …  
    463491  ////////////////////////////////////////////////////////////////////////////
    464492  {
    465     ParamToARIMA(_.arima, MatSet(Tra(VMat2Mat(x)))[1])
    466   };
    467  
     493//  ParamToARIMA(_.arima, MatSet(Tra(VMat2Mat(x)))[1])
     494    Set arma = DeepCopy(_.arma, VMat2Mat(x));
     495    Set arima =  DeepCopy(_.arima);
     496    Set For(1,Card(_.arima), Real(Real k)
     497    {
     498      Polyn arima[k]->AR := 1-arma[k][1];
     499      Polyn arima[k]->MA := 1-arma[k][2];
     500      True
     501    });
     502    arima
     503  };
     504
    468505  ////////////////////////////////////////////////////////////////////////////
    469506  Real setArima(Real void)
    … …  
    484521  {
    485522    Set _.arima := DeepCopy(arima);
    486     VMatrix _.x := arima2col(_.arima);
     523    Set _.arma := EvalSet(_.arima, Set(ARIMAStruct f)
     524    {
     525      SetOfPolyn(1-f->AR, 1-f->MA)
     526    });
     527        VMatrix _.x := arima2col(_.arima);
    487528    Real ok = setArima(0);
    488529    Real _.n := VRows(_.x);
    … …  
    657698Real TolOprProfiler.Enabled := True;
    658699While(numSim<2500, append_draw(mtm_arima::draw_multiple_try(0)));
    659 Real TolOprProfiler.Dump("draw_multiple_try.csv");
     700Real TolOprProfiler.Dump("draw_multiple_try.02.csv");
    660701Real mtm.time = (Copy(Time) - mtm.t0)/numSim;
    661702WriteLn("Time by MTM simulation "<<mtm.time);
Note: See TracChangeset for help on using the changeset viewer.