From d3fd2df1f58546454fa6da30ebb4851107db2891 Mon Sep 17 00:00:00 2001 From: josemanuel22 Date: Sat, 9 Sep 2023 11:37:22 +0200 Subject: [PATCH 01/12] Initial time_series_blocklearning --- examples/autoregressive-process/Project.toml | 8 ++ examples/autoregressive-process/README.md | 18 +++ examples/autoregressive-process/loss.png | Bin 0 -> 22347 bytes examples/autoregressive-process/model.jl | 132 +++++++++++++++++++ examples/autoregressive-process/utils.jl | 41 ++++++ src/CustomLossFunction.jl | 2 +- 6 files changed, 200 insertions(+), 1 deletion(-) create mode 100644 examples/autoregressive-process/Project.toml create mode 100644 examples/autoregressive-process/README.md create mode 100644 examples/autoregressive-process/loss.png create mode 100644 examples/autoregressive-process/model.jl create mode 100644 examples/autoregressive-process/utils.jl diff --git a/examples/autoregressive-process/Project.toml b/examples/autoregressive-process/Project.toml new file mode 100644 index 0000000..53b9832 --- /dev/null +++ b/examples/autoregressive-process/Project.toml @@ -0,0 +1,8 @@ +[deps] +Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" + +[compat] +Flux = "0.13.9" +julia = "1.6" \ No newline at end of file diff --git a/examples/autoregressive-process/README.md b/examples/autoregressive-process/README.md new file mode 100644 index 0000000..91bbf1b --- /dev/null +++ b/examples/autoregressive-process/README.md @@ -0,0 +1,18 @@ +# Autoregressive Model + +An [autoregressive (AR) process](https://en.wikipedia.org/wiki/Autoregressive_model) is a stochastic process with an autoregressive structure, i.e., past realizations influence its future realizations. + +This model-zoo example illustrates how to use Flux's recurrent layers to model an AR process. + +The example contains the following files: ++ [utils.jl](utils.jl): + + `generate_process`: generates an AR process + + `batch_timeseries`: transforms a vector into the proper format for recurrent layers in Flux and allows to batch the time series as required. + ++ [model.jl](model.jl): creates and trains the recurrent model to predict the generated AR process. + +## Example loss + +Running the model with the hyperparameters currently given in the example, we obtain the following train and test losses. We see that the model begins to overfit after around 30 epochs. + +![loss](loss.png) \ No newline at end of file diff --git a/examples/autoregressive-process/loss.png b/examples/autoregressive-process/loss.png new file mode 100644 index 0000000000000000000000000000000000000000..a2c97a5d226db75117daeb6d7dd4d92d7c39f678 GIT binary patch literal 22347 zcmb4rbyQW|7wu6&E+rsHr*tFT2uLd>jescK-6g4nAl*oJcb9ZYH%NDPzs>g>@BjC% zcL3wYKKtyw_KG>@@`t>v1S%3f5(Gi0Qj+f!A?T?Z1U(6ZJq54ass0!Ozn*`Qk$49^ zJpTLLm>mH@q>$7*VI`-;-Fe5)m}?IRha0u^{aLPZ8Z)fY)2>m<=N7_&-HLyh18KNh z7z1?>eliELMpBY~e?MI=U$9Fie{O*{&nFc@&@8l#29JalzdEvdcbzX)r+UDcPyUg3!^t(ttZqw?795S7f;7fr->;mV--eUEozgBdry_TvG(}KfI&nL zh2tw(X>?!kZeMaaJ{%bYd7Y8Mo znz)FF_ps5J5xWh}Yl8CW#E1ynl}=RC#JOSbF!bdad!zO2p%^Xx+`v3$M@9(3CWW;z zH{Xtw=!(`N|N3BHkagzaCP6?>;VbVq(zqDwt-OedaU=W;Y9Z(TjQ&!DJ!GbQe}U4| zBDx+deu?8xJYQe*izm!r9bZRd_h`W&Xo@2MRYdGO_Nhar*?`*aZu5wa%EQ9pk2MVjr#ASyn)(cj;~ z^78V98kN1x+O;-o?alc)Ih*wdP30f_LOY(|j4k>o)H{W^@fq&l^oB|9)<1{j?LF_# zRIhOCcP5f}?YBS53D|E?`n0yN1WRYqJ!quLTg@#BV-ga&?a$SErbzs^+l4Bo`EB&; zmy9AJ-&qTyi$ot9C`W(L^lBkPUe9LA^nsaa`M)5-$HTiCSJffnw)!ob#Cx_giH3+M zB`r;Nz%kI>{r&rQ8lfiZmCn(u5BVRnW!2QwLWu>4rwdJ#e=fde*uYMSbYCi7PNM%* zPt^kkeycoB6E%o6wW(2fm)n=oQA8wnvV^Ow3BNed zA_t|I-jM@>lx<#zy1Dm?*Gt}yCEE`VOAuq@;qh4gdn*EKXL<7ENp5bg=k1Y%hey+} z5x?8z;aykr>4=PUW;pst?#JwM!y%1^-Co;xQVL)Di}UyKvGt$h4=$Ox7eCkbJ%OBT z1VoL865-7KAsei`el zuVcz9wC=5&hx zBvj8=VHA;6QW|3$;+9=!F+Vq_?Xs15{OtZDX}QUMxjA$m+WatZ4K&G~VY1J@#@c%{fQYc!TKoF*1z^>tyk+VX_ok~4t= zKLlQGKbW!2X8$f_!|kc9$1Qrb$y3P6V35GqSFt!-M^&uTz=LPO?!#bQJ#fwB8~bBl zaI1-jXhhR(Jma_Vh2Bl-ddNtLHivW3!K?HB?Fa8 z_VKLS5_7a$>7jzl&i5mkGK5|3M;$oDP1oCb%J~NUv2R>XH)%g+p^Mnq+9Hx6V>7@| zPztPe6!6Gk2kj+H<}Stw_DlEHmc8|LoTv|s(YtE^P)_Y zVo~B(zVxa^i7dKp`}6gPWXZ|NJv}|XVfY*-kJBb4{cSjue1Ct>Npsdi3ZoZ~GO4x@ z)v8Bf;&W8;WXLoceJy5{ktM9DAHo~hVwG2NQFsEil;)-pfM<_ZhfrC0xm+^8^Ug%! zalqHEpRdHMdn0Md$jBZK(N<>6%a<>)88s*@f1!Dky5aWPY{z1QXZW-pzw4jas^-%c z6+yP7lNU;G5_I0CHpw0Z0bDyveKpPs*v;=DE5$c&-hjO&;jvj27z!l>zq)^2U0o?C zC{W8L!qCug#1ybPsHa}^yBscP(JG~mjae4vF2i|D90lm5L==+M6iKnxe^A1USPQ14XgTnF=pI@=VWo`tb@D7Agucy$UNo&l z0|O##sj*)9^5qMR?C9tS98Dc?^@T=PA8+qo$7ze{vNeBH0#3^XLST*I;mccFQn2Yt z)5)1iQ)ZneGuAl0qWcFj7^AP^#C(GNMH3e(BAKrw6?uy5Y@&5fq=&XALxDot}DStQ2U975A?p%TW&@cAA0Vj=y} zrPDuytgs!}zm3N!MB5sV(*H|x?hVT+z>7HcOd`Wt@72PBygt4~pTWpes_#TYbaEwL zS&lwib86S|?AF47yjDcfXNbh}Z4cOHTzALXP?3+LWn+3#Y*gW(f87PF;zq1-3$6!^ zBO_dbgX+=%1ntAfA?i%`;JDp#+iYk{U z&%Ju$l~`rD=y8ART$>*@aOdb)Mw$;tHgtm09!$W+L`%CxCPPIjh0o_uP+#u~GA)kY zoc$b;tzRsI=<_U=OYB+!aFa|vYeg_<^$0BcPy*M`!oH#D*CBz6#Hc~DDAF;L$`p@l zS~6e7>#%l@%cTI9vuy=w=j3W}yr06r#jqAQ42F^%)>iW)G`_=C*$;$TChCB*Y;<+3 z&9_o33&GF^xh6zSAMRquaDABgN0O!Y%NX|Ar($>Gfg*|d-xbz=odqC#!b}EF^G(h8~h5@ zO4CY8SSTpKQ6^VcSNA4O6m+k1+nBgt4EvVVjV)#=5oDu%hF*u=UJ-D1-%_tr-uROrGB5fb-Fc56_Kmk?oZgo zs8K<~#pO7XA@#qg#iu zO<7%^Zuy3p&sKqSV{L5>tjWBso{>=wDNKc~UBMW?iu_$~+}+Qing_b4Cf!;uptl?b zg`{&8#-n#~vWeVB_tc9{E2#A2|EAbcXE7CekdQhBd!Bw@8iV**F39RWJS!FrJ=4jUjDdI@VY8Ojt>)w=QEDPP1mqqLcP z3d!G*!ic_^ZV2b`ez4v65Y?}?gbW=zq5HmQOMf^i$hKSc ze|C6;ccOb=%+}9U-2Iqqa)1%u28PiLwKA{WpQZxG_*gTXzYyMx1@~*3U_Tgl@#{~* zS`j=5>0~gXXw4yvf&jVuK+DIp|G$@TGPOmi<19ac$fvOns>od&@4vqkGG`3?5h&|b zqLYM}9*1=3ed2RiCA^Qci@3A3!SO|Ue<m$eC>*ZG-bgLP-Jj zQ$5OH(!(2v*YHqHN;= z{DX}J=DCu;c6SYr3FH(;7YZ1whW;*{*^ym!yW-bkA`=-ZB3Ytk$wnArEpqmcOtHsh z{>)6s3Ca^~y70Lm*HosHbKTS>6OSXY;oVdB%=l~bijZaC-l#9lAm}a3dKV1|7xQW2c8>b*hTvuQpPFTuvh3ma6ZJ(+ zXusM>!qmC`$QqQ9i)Oz2;j4ggFm>UM@OZ(T^33ygmrytQot+qE413_aSBsuQZw&@t z`mW)-j&?41cFliiw4;#V{UyPfA)}Ls&Wt5R##HeOXOHk*{qs>tE;`Fxl@6CnC#ny< zcUQI#lZLEJ7B(@@AX4fLq&83Fk&)jlGhDj|0S%H&*ir#D8AMzep`kzMqW>U_?zdK` z(xhAe(Fy0xs;#pAv@o(IIOO>*hdf&lYH8=W#MNE1M#)_f4v=^g#9P9Zyv}K7xa0hH z+OXrF{`xR&+h9WY3S)cAXRtrQ=z@~L2=0{8(+d$YGD_@+X_mxTA~8pngy@4fyz<1X z1sC6YD0}kXs}A*QAwgbd=)R%oIh2@R$6OaG$mrO(;v2_f$r z%cG67Ql)(j2wMF{T2N3&d^)*!T z=^U<%Qo|gvO=WW5wR#?`G#EJhv?ef9kuoO_Y5jIKj^eidGa6G^8{IH%x92u1HD>L> zY+1df)KcF7a~4wealX!tJO=vug6+HgHaV*y5y*-YCepWwE2(OyD@RH9w4!1Io6O=_ z^KQ5rRkzl}ik0!KCWd^v)rv?h@)M|p7A}<%-lXRyB#vIR=0|D&?a7xZo09|tKjd-4 z4!^E{F{sq{EXh0!kEfx^hqQwKC)u<1xIKZ|u_i*wjJxN89+=I`HXLf!tw#3RL+=N5 z&MW03|3eTsy-6$5Xt=JOkJn4Va@FJ~9n%O_yQQ1AC z5c|T-c9tW*!*xH)UjNCuwM4s}8pmAyOOvzMqgN?-A#7xd|1&tRgK%@`ZYEA3k&=n0 z{MR&TR_$*8EgBPCYJZ#EX%O~q6or^5^NcgOkX{T?Vb=vN|AN)A+)AA2Barb~3?@?l zFOM{WfHlUSCjMrmDB)6=Z}0UBybo6#sX30#%2_yAj#hO0)dh7=-FI9?2+^ z&!G}3INh_ih!VLG4g&5oI=$c}tOSzP!YChJ#gIAUU~mb z_1)TxPM(XMc^?pR8~kx%J}QJRvPiC3p<&DLBt&MS;l+P5BSNLrGb37mD`YXiQF`Zc zulVcLXxJH|4G{~c0mkR)KDlb**`m+t~LqMzL=)OQZuhN3;h$*ABUHTP^Ga0P`3$RQah`&#?99I zPrmfV74dSl=1~@kdGaoyF7;M-Pah9y9ywg<%E!{x^7Y-&y=mPZ?X*iP2db>>ij@?q{6% zxr89QbjFDz*izJY4FDBFn00;S^dtOA^7_IH=?$Txz*1RiOH8lq;v$CDVWnDz|iFvUB?BzZWZ8KG=NJVB1Jx-zv`2 zJx}!wQ!sQcMN7u!-SLtFF1o9aW8B*f7fbr_z|B1qhWVnuT$28HW33vPc|UvgvRr>K z7o~V-;bJ;5MGz~*r7sZ!YC%HxJ>){8oE}R8(e!e_t;jDc+(gq(S!`lRhc`)5rCI1T z)Dj_#-X_HD-)}y24olseJi)`BaqlEQ5-Y8!G3v4jN<@WPmeGB!aHP8fK1|_kvse@O zBKu)?i}c!Pr!UZcnAdyHCx2EB43p3so$gSY;Y1fRYm{oqsR+I-%>Lw*_)KE>Cc_& zB%w(2Zwz+E^B=Y9tryxJUMl2VPeTv`yf6%Hlr55JyX9=YlceGuwhM~WHx{$jx|BQ* zItY>hH${F5OGW0hx~o!8)9Pl?CeMmY;?sz)8-AGyCi}KS*lBAd6QD%Vo{X@zM}QFE z2?tF#vgm(v*MePq|H4~;S~S01O3X7LOP_EN5~?)X87V{VVpJr6-o8>M^mV^EzuT)` z>}^gKbiaOthz=G6j7KsiwOt4NE58wK;Yk}->fczb$*n4LAoCnR#R`@ijQ5E{wZ^|@ zASf@*o0J_F#c(hY%mSJ0U3dnlD0D#?1YibDQ0ap*MOx;~&X5!;ilIS;fvfo|CLC1` zWMz)8^DZbuS4YoEE4p7-K7sg4k-S-igi_F#4-VvDQZh0!a&iD6p&-l9Y4!lH_V##w zRU!Ir_eYMn3w$;$F2w4*2G<;!1nu~#2_D;pi1*ly94t`F*E~{KV?zT?CV=aF!^ZPf zP1)lo3N=!n>{b~=MRK_S5(*m&$BTO_ZdcvU34>_m*ZdM|rqH>k;_|Vu`S-F2A!(ci6&eU?r4Q^L908g&A zm$9Dur6oeuX}Jf3a`jYsHvmHdB>vTN zT`8#_vf1^R^)=eMyEh^E*Ed{?b@6OFjG?W}rrXt}lIBre z8w+$&+L=Oi63X?ZMOr$_dUp~YL3*5ca(m1I743Fs6(G!LhAtH(WTQ&oalxeqqx-5T zJ6vI|zuK9xBVrNj-oT+(^JwllxUG2tlFc79a>6H~Brbkr%+lX|$n zLnNE8cN78#V=VV0z>kxE{jxq@m6E2P#XJwgULwL_`ed}wa23s{sZp7CxmT=OtW6a` z%g^6%s|!e;=PwAIHk2!1`g$z5SV$e2ibXpP2D+NLwx5+JagMo|uo+{s2u{n&WIqNe zU80L<8}#JGy&S3+qF^M~#vXb2bL= zVXoWddh4b3avL!x(H|$LY>!#RYJ~%9VT`rKKEuUTHl|1mjny~%A1A24^Lqrz=PC)Z zlG1f=+}_@Xhlc|;3ta>W895>%0vH94$ym+{BHpnAbvANxa(4EL@K^jJ0|{KNXWM7P z5+njN@Q7p#ed6NcoF=CiUEluxm58)w=1RL@S(A@*_u<_ z!JNtql7hwh!IBHd;s-Q-5tyzzC(Wl3l;weQ2}#=|Loq^Vu$c<5_byn@qp zCMf0CaM+(S=1i!!-w}r)zIZ|9t7l-KK?h@FV_R8SdE|?9+k68Evy|xS9S`n*zY*~h zNtLGoQ`_xe{M9XTm8cW>vsprL#YvCDi&jbY2UF^KLK``0Nj4~(^jyH|RuXdJP zA9Atwu3=jnGAa*%@U?L+Q9qB-5k_}$b#--dDFVR+Kw!u?Oa;J!IajZS5OBc&dmwio zrxTv>4B-W!Mc?L+8yFgPb#4-{}Y z;5L+MJbJy)|NLJqK=S<*@qyu@s{dr;2}Sy_vmi(6?s2ZSa=8tTTYpYPa1cSdJtYHq z+<*sWlbfV!VAt;>bvCN#jDnBVTjp0bnSt*UI785vcXni8SNAsueJQtF*(q~#bNvRk zev3At0J-bv>XML`{sp2PibE}%A4 z3N_AttAYp&FP!sU>a(>=$jJ}>#EEN<(6t&FU8YK-7=8OiVF(*jo5p&{O;%$Zf^FNwXvjw zh&^KEHOvN^)E)K`kLDCw<*&O(<83V#+_j?I)zrUXkVbiV- z0TuJhmG3-V3d6~P7>D>byP56_#IFm+dfz-TN-Jm_Tf2g1I!rJ1@?4IcJHm4r9;y%e zxtN?oJz`E!i(GW0>P@a9$8`8Hb;j6#mS1kYQ5P4FWjI04+S#yK_xN|Q&6Ow@wfkL@il4r@cekO^XK z#y@y-C^Vso{;9%{5jE7hP`78Go0!{+EN~ext3OGE03$|>i;GJ@Fqz7!XxdcA6&;PXGfU|T~$2BKCnk)7o}Lf=ll1IMRxSJ z`73&7*1&&~&67=iI_DqGgQY;Wyq;Gj`>dr%N^v24Yew8!Wgc5}wki}r>bpGj9&$ZH zLkEE8F)^VddzUI7sZ;53cTK=$QSomE(Ml<#acj3~w`H&G+r3#>Ri3&eY`Ei0Rg8<_ z_AWbIz;ym5c!b^PqMHNFzxkdfo2jOic%U)JR*@z`)%zAVnU)Lv_aB&PYd=;OPzHN0^-iPO@0p>*%ucCcgT9k*m;2T=K7&ProoSyRuRB2?qt3J%$J)5l*w>{ill%Sd0C)5&V@3D^YE%>)vpMwJh1 zr`Oijf<*3~C((SSg3NcW!FjY(7&eYYsV7W0yX| z=4$_VR$rM11zHY%MfLwmdzJ5Th<_emXK~`HJyNyV3qRJa#WLWw60Uin>$YDQK5u#| z;h6tb^Q!4GH_c6_!X#%I{){e)j)Q}PlkI^?}Veci&Pi;@MI-Bld+W9jNrzq--nIuYp#3%7bbc?OETlLKntr5Zf$TxP3b zZ^VtNi3PqH@Ju2c zwdtxC@N{J#5q5ga$rLFpOW<*wB zb{|v9%j0sOjd1BMojWg`I_+!-r^w8Be5`J4)YsNlnvCaZR$G*h%@d`xEYJltN!-HUB5Dlj0{tYSrmz=*j4Y%3eSGoSqIo|`Xq%7Huh#X!Jvd#) z1ywP#5y~792_+40?|Ul0yiA-vpleS9m91(V@73HaPx-?7t?0OYLar%S>HNC4ia%)q z1Y&mRD_`F*TR+9Q=N`6>MN(0YiPlQ`F|ybN$sI;^#E2`zf(UX8q>;TS{+PL$0S>X2 zpC;JNvdt&dUnvZovwfzFC^hN}AMxRvy9u^2o}R~s|1L7v5c)%LZ^V_|lTO+5bAlez zXw_!t6K01z1t6QRKMFyGTDL2e9?$XEjs zumRHqs2;1j$v2Zy&wgsK^(zxMuihV0;x00) zbD?rvkl06j80OqRycEi4EVn=hS{|E~&XCYsf1o~6J3Nm1Mp~gM2blOU)eGd9zzlmW zp63H?2$*M-Vwy3VXa_dN}R&2^fz>^-alFi0esLZYRWv(J>%y4Au6Gk{auXR7_9k?X|iy+|X6 zBjCxH!({0|9JC*!Z^WUgmdTi zH+HB%*!xakInw6^y6jAKKM}!&7;5U}hbJ=qF+H!FZsF4%oXe8*DavV2`i}M+{Y~-j z1)eit>ERq`fOy$uxGT|H($$pY4JJ!GrQ)zu!?Thmub6}K69kM9!xt_>H{Vzr23L{` z>)xWPK*BWX)PP6s{?It~C5c7o=N@*j%2GO_Uu_4U$2iG+q4CsUrZbMd6H4lshA>Uw zdI8F4WqIdlKra*FiktV?8^*CS4H{8K=ShN{aEg_!XC3vEhw6BFInR4bW`Obri|kkE zu-||BBl$2!tTkl6+Ar3`Demo&9>b1pF+lp1+lFnLm2YUf;khd0Tm^Wth0->_LdB@e z+Y@2-%4xnzYxu58# z1sea#En|xkIkosy9rhSfqKaTSt(*L~A;`+6c$LWM*0rWKArxRu-1Zp-xjuH+4;Xz? zlk=Py_c=S%{UiXAL;EBHM^3n$UhQYtb8xq_enyu_ii8G;*hinD48*A2 zlhhqL%t5s7+{zqk4SEvqw{T7a9JfPv4L*H}A{!12lhl{Yc+I^{W|^DT^p1fMEAYhc z)-07qj;0#}sqAV_Q)G~_tt3&76sV_QKb8E&X+_d-L!LEp$QVRib?TQ^$g^*Yea+Py z3%U+>{i!?Oz?}k3@BxC%?b2j04=z=Yb1Rg}l%RuWvMD$EYx;-kCuUz2-3>gW*{z#D zFU@}B=G8ZlBo}*KwkBvFk#t~l_^kA z`E7J4RQ}EKe2>xxFGgcB_ImAtGK^wap6-aM=Mz=jKzMDiqBI9naOc^A;9IDMz#Dxv z<73#k)^rW3o2cc?6oKW0@>^Q*1Sq**i=CM%WVPS!5lOq_dkY{x%zu26|zy*CgjQIM0b%+=bip|!WQ0l*79 z4hYe)*SlS7y^~_cqs^L!a#(;dM<+@;=!^bE!++#Yxi+J6A5r>`YMp z?vs|f!yg2@3Xud_bH7|os69(zn1)bG9dT(#JShx_x&>@jzsqLJNK5DGG&h;|KHbEB zzi9mDe3oi=r>jJcE|H4tITJ6b&>!X0fd2| zq^A;%`xmKQ)Ro8-7VdCFP-n|fC;6m>cYhmStK8rNK+CrtER`O0+5qSfxH}zHE&3$& z;R7lfS_-#y8ljl?$#oTn>B}Cz0Qg73I^awWJuuB_sm%|-9P`cwXWQdUOiYiGcr-N2 z(X0W(6n6(8c)GoHrY2V zZB0$2kl9EEA{m8%zeb#-*z;9;(I{nM*o#K;|6H2M; z^#n^8S)=PklWhxNCTKHLQ&U+Z5oL=jNs>sqDEI#jVUZ#~X$jwt=eR!VA_08D?Be3$ zzyKo`SGDD0)7pIMS~RnccC&}OqoZR24&qxFnS0x24CC~Ef?uQiG(K-XPzv?c(8B-_ z7@e5N4{%;@Z*PE2%f+*A4a#I@W{Qc40g$ND*$;Za#qPk0m7E-CTM$fa7T-64TFB<= ze>VY8ZEvQs{(&0@Ck7O!(-{lBxD<2_Fd*kO%B}VF@55j%XkVUk7;S8^Crn=Dn3g4LCohLl6mzjU*78mHRy}S z)*s5EGi?>B_Nv~1>T1{qFdm=7&r$>Rh6=+HaV&9NcQ%e-Rd_r6lhX%NN znLg?dlZO(ve~q6&q1nso>EFJ6TN`M-ySd;po!AB76^Wo*O$hz79DX5Gb)uoBMc=K; zs1zGXmULqn6m?#*h6xOLy21Gb90xez?>~Q5segdv46NX6jUYz`36tvNhZW9ve@LG8 z)m-YrXBM-m5?#>Az~{J+f{5sR4In_z`>XD5$ypqP!x|B-N1i?HWO);-RZeultmPE- zEAD$*crS!{dKj=2KFfs$<$M(kk(Zd500UoHUf!R}#`&%B+ZV`1g1;N6=)3}wg(iDR zxAZ3wN)f|Khq-9bzkO0dK{qEA6%`SYw%@`&cDzZ`W(z=|D6oi5LSor|3?-wB_WXt! z019*`>vUXdY2008gEY$%b%r|lNTE-Q5Sof`6`~+O@y~p01F1U|Ui8ZJm2b3rk(T6B zbIV#rtM$q8LwwERTuYDxvhMOzBU_z3gS8?wXw^X#{gFw#-LHpUbbb``7kks5CK=bV zUq`!FN84mA|GqB3TU zof4rY9eMa|%LVthh`EEpIj(%Yyu-=~jWo7*KoIp|(>iYIZ#7kP27RZLc zQ{7hCZA*@%_^QDYCbf0oY_MD|onr8I(59OjwQ*E-~u;W}nw0KlREJMt{oB z2`}0?p*x78a%n2*qr8rn{f?N+<7rM6xTIXSs8m+9-otK}`5<~i;o3w0pB6BnFDYD3 zI1N`_Cc=d^+xnKC?jH*>0#~f3N?b4z2zisBGf{H2VLaoW?qsg~^@JLvL=Il>v-%ij zsZl!XViUw1I52s1Qwov%j(>~(R88?2bNonvYbJ3-iI0RF?b}K(Qn(Tp>F)1F>!+lC z0fB#tSUuep*^z7};UF}B@E#Hx0|zJrg}nO{FxZlYaZODnwE@ET|O6AI>*hpROy7}Ey$prM$Bx>u{kg7~ zq`HFKm0dOOD&P52o>d(*APQaSj+#RWCR-Ak*QM2{yNPUHZm9&v{9c6UCFe=`R7^D! zyd^fmc-l?+R#ecaEdx#Yy=iB1cf@w_GN0wDColE-=URBLfmAUVDM562QEP=MizZb` z!brzjOy(q7Y=YD1x21kBEBtM5&evXypet@Q`%6!YmBDB3H%PnZ1k3UHBNAP{FjZcb zL_B$)iLSi`}~=sKEf5eMIa9MW8PLq!}Y^%sn}SF-dk=WH6P<835*{QHZaQ) z2?oj}AEOOY2XRvyb_wt?!gAGa#)xbGH+Lp1THaybA#5M;G1JJfEWwIgv>WS>HnnOR zjvubKeOVaPw{creUr@T6guN}CqkO7l3ls4&NzbBB!j+-!WOdeGz39_Vuw=C);tbd_ z8j_!(k2GGBP_MvzLQLX~{`>k!@mDHGPYF$%Rp6caNS58vnvkC*ZXJF4x)4MWR13JhYG{{WI!Q_i6a@M#)xt#4=e~_shF%i zJ$6aqUF6i)+L36Ktw!>;xTm=?=%e(J=Scc2@i2L^IG>W58Z|&|E73XkK~@wfj43*z zApcm!?TJRQ{#tjsE^I3x$4lcag_7O6$qj#rV`?AYaC_Vb_@MLSHf{wiM=}${fNZe0{wN^df)DFPbzUY1;N57w8Gix`121kEg(6bSMf>-8rv}Fo|52+!2o!SY*x1;h1#zaCx^&B~4x|YPM4uEDkCy^SJOI(Yh9(h$e>P)UJnyVr zJbA)D#1c!!^vtLAw>%B^3)@hAc@>pTmVQxQExMtguph&#Q8eL+Jq&e!n(V)yzrdD| z2xN5sYpC(&v9oK$r_;S9&ncNJVZfVI>DZ}I<)+&!Y=g0Xaxw;8QB7@gd07NTOG9&b za^m{8)xyG}++aZZKRJTG?wAROw>&mO-)q=fnmi4|D__S`N6Nc&$*NzsZo-6~YI-o| zrWb6l8(*!w z?vE3m+jZC-CZs)9Jt@&6yNagh(ZcnLQrjuNU1|>iOMU8cceK(4)FeUJjDU($Df(0k zn0^w^+h6Fc7IQUdf|sOpblYuy$RK~^bGs~wqI=>+Ju2u`h4%TRZz?n9&yO45m&yE^ zL3qkR&yC+?KPL>V^tS2=rv!o#&_47(amo;omw?tfCAz5meABBZpsj_DiHVq1?|VQ% zSDAiaT6+4UohkWc2L^}gc6PA@*k(9y{JGBet-1#S7~=zZ@oyN2RxqWqn-`&f{`>)b z0s=n4+A_;oocb5~Qp5E9 zcU6eZ+SsTKb5T$Dsq9Riv_yObopBxhYj3J%A2oxDSStX$%3&hS>%cS zJCa=iUt7~nu-sk1WxzPB<#30>tyf; zkS+hf=>Kk%6(S^-LOX>E`n2_I4U4uy6peetjCZwLDZ-bA`@hrR?+3nxBM_dXfx4Ne z{A7XR6LmS--oQ3kG05JnE9BPusR^?zw-|%jx{ffdn}Kb@n|+>^+C%&tG?I;oMhYVdkXdX(R_99TigD=hJ$X_ z(y1`F-$*>SBB(D1oat0cwbsy(06li>Y7q({49qp}+qL#-U$hwzHEW_y?{Q^(Ut3S@ zX^TdYdK*Klg!lLf)Sr)4hNae2pD5Pe-PTyT$R_KFDS~hFaah3ND#$M6@0{|wgWmKv zdCcuEBu0}zt3pOqJ<&UtOyMu`NdnVoY5By(eYcLnFFK~B@v~Tp{oJTyC?4H!UUm6D2`N#c<8yM>JXuwT5 zbM~4fh~Pj9zWVx{)h_^K242qs?b1i%OD7|+fFwQ>r{k~w`i~pG3O9E2)Th7>tkeia zpZl{VxLO5V+w#J1m5Ldi@Ch9dI#~DsEaFMakzf5tXhNB4M2#asdBR905XXN#NAjGD z{wo6K%{T$IIG048B4sv;Ouv6tf0=gKM&d0+L7pCmG0Z-1Np1(a%r_tZUx4=kNwbMK z8eqD1w;VtV{eFgYKSp-f+-J;=k)-*@!`jRBF>YG0CfvpAO?CzZoL_U6{p$LMaQ=Y# zs-7cKm&%0zA>>o9@{gr@R&eE+2#!wBMk+;PHBzC+2Izq}1BdS02QXvkUPV+Qn^V7- z!uo*Az_XdjkBLog>Y_yXszGu>3)9tdE#-r7wh4!S4=V&Y*0N<&HAcp{ssu0onY76k z5RHYYzTB5?w;%JNoRsMq%J#k5$cVRCy_)F{WS`h93Pqd$k@pBryP10vDY|-izW>-UXhzE+L*l@6<>&T*8ojDPatCP)-Ycy zDze)Y!{am1CRBufKFFqEaSEiV8Cq>iPZDKyZ}3y4|ArwzPW()S3JVttJnZ4dByYQR zJ3i~4fusg12?iEph@9cKVU3()@?0w&NHmL-&g7@%i1@+=5)~>z zT)qk&Fl8txcS9<>T5pBxWM6u_4I|0BhE7tW>GYHNdOH(6sz41Wg%;1{6@aC#D4tiVcTnO+v+O*Ief{}uznJNwA69@IE&54%eX zvRo?Ql#JJ9#&^@3G0Y%FFe-rko|AlMWn9$RQ+iwVi6Wx-bsN ze@SG6C=F_nZMZk_;h~)8K*;@tf7900uO|Rb`1L(Leh!%j{`D7>QAXUA0A+8mzd6Kd z^TyI7KheMZlgbjV(P^18PnON$vW8%D($OfuwjQH!n z*d7(NXy*Ge?RL0=ZtUQzUgY1cGV5J!uk=^tf7-FgTwkTJ-^-;|A>1{gY7SvqU$Q(d zO4m~E>v?bED}KcvB|*&JLFtuk!U|kX5mIeXxhgI>9Ax!AG7MuFoOM41QhU{7A)OYF zg;$^;fw)$cVO+S@ z0Yx?R4~yHMpQwkBp!XUcFB*hMegr(H7OcNpgezy9v$I#m;gqjEa%NUmW+5aI z*;I*`#yeu{C@YBNB5EU`~7*nU-S8VhI*Fx z){e%{GQ)x(Eyf??;%Uz@{I(LBw(xxidSjoBa)cPsU;rCT94TR+?x~-OkW$DH^V)HA z03`_CFcJ-z9AcgAY0>mc_tSaSu%WO=VlmRETW{Hl%fnUxVSWVh*v0d{8m4n~y?<@F zWAbB0qg4Oqfa*+c9>ugakPtVTl!CPteGZhXcU3h087g@F-+D4XzKLC5>ucujLwyld zd+wl=Gbi8|2lf>4wlE5n>hhj()0l?zJ0cd!= zUvH8C22Cni+d!67$)E!Vgit#c0|nkzoy_V3tpo(a)tZ>YnK-kzz+?jIE*0WT=8jj7 zzR8z4>ISLK;_d5yC!P;qbC;m(6e}C-3InXMhs4BL}V36!sm3Zq2yaDeTx$zW~nSekh5vTDOa-@VM6_<46?7hkr1{qF$Bs^89kCouvCY@}LEd{v{Fi*cx!v-BCa;q>?C zfc~s$!2_nrp;!k$JDvvp4a>O?E&SVS2d9PDCHm68?cN?!UwpPNBC_umGLR>H&TZ{hq^tSdlH3SMxd&4=(}He}U9rW9?$ok^N6|eYdj7 z&B)4I9(&u2lrP(oX0Bej-%2Zm>U@|tQuW-RK?%0|Kd$-Ml^s*m5e3+G0 z2mI-`Z{OC}*H<%a|CkQnS{h7CV-5-V1EmXyWF!c|vYhMhZ@e?DeDiS5)^GB(inm}@?VaO?RsoL7RPX)o6$Z*LNjNHi1ILHa52Xlk~WmzOV&?dq>9~&Ek9-jsJ82lGS?JkcwIZiKFos=nmj|8CxlS&mOC2`h-#W6sh#zaMJ z?1p6;)?OWY3T2XW6G{EnHdbAo47H_zHB2o4q1a=eY4~JlnG?12TRXC#K#)DL)9v++ zjyxw9*CiY$mOXL(jPc;ol6cCTJ?TFDOwcPbCk`uV47LU$|=Bs0CI4m<@P>(ox6`ui~s!F`MxX+5daA05>k(ZYj zDz`lbm_;R&j>H+VOhnM$c zUDdsN_ndwiSw?hKz&jz4Oj<>hlsG${ z!7wd2Bm^UhWtYO?(o&X^GOUSY@}KE>{~KmzX$Q{^=6E-yBqhOCY$|Mj!fffBO9<_R z*Ne(Z8F_gtEPLIH7p221PIq=`Xd_ELZLf`2U1*qHT>LyVWa;GO{IUtg9RM#? z4pZAp%m+&D%RSKmVLpc2`kZSugZ?`rS2e zZW)wW`8QlB3~O?|aPGkOEYx))!MJEt>r#?WE7&K zR|ylVk|A^-2p+(N7P$|6gxstweWwdzf*=Ql#jZAhs>0&euYT(j$~@hfhN7w6uL2es zc?Tb)fmjHPve13sKUWRBdsnQ;R|IZbMO8K4?ixI2$A?DCZD)rc9^96aB2G;?@^s5y zZ2Y#m3KSH_pwyU~sb<7s$Hu_;?8*3EH$7clHa50F-g|^cn&BEZ=KFhl&xY5pEi5eb zWyv0bOon{4I}32H?BQc$Q!?lgmrlpT2R$d8kPwyNn|Xg&H^VfQq57S_AYv=A^z@0<+G1C0$I{bTer5A>ah(ktOH53prlyWBo6*+R2486Z3Oc}&G44|{AqwMdL+EhCR{%S)|r?{tg6ySDyyh$|Ngxk zs8U!C3Md~j7b3$0_qW&JeF3UCL-d=bt^o5<8k+W#`4+Hm@c6(vh{N-ApUk31{e!Nt zt*D#2`5iQ|-gw&CapgptnLM$G@0HN!9Ws>aAafzCbfgS1j0C|teBW2 zL#rB~+jyYz+<*X8U-}OTjP*oJ#l5_FA&pvd&^o{4i@vnAw(X> z$91}ZYM6k8#MQ<3xnM^DK|#^^k*r*Nd>Xg5D-Nhxii(OL7XT5ZOV7XPR3-OrjnT>e zI~QtbcwCNwoTq+EiRV3ZO<@WdF9K>U;8^q!=jrHTq32pQyaE>b*6kEeHyA1p z22*snC6z516chwvCbw_j#uwnk#Pkkl2OvL@l#V0x~j;d9HbQAj^ql$)SW_Ei`XrYpSV1B%hp`TL1Awvw?p$ zCYO>GMuhj8{XkPeEU={de<2}~l8v?%YeA)-ytZ~& z_6lxh425`cR|nY^I1W?~O-7T0ou95e-D|6>UDgna#WM|MvijVgcaN+D0wK{67$4xy zs%twsIwUg_y62Kf+}zwZ_qNKKo7GY2bjA#S8EKrKip$6>{0tQp6@_F4La-X2tcXZ) zQBmYI91i#FavwWfcHo*gIE;f2ZE4Y{f!xPRhyW?w-80VK-fOe7g_@T@-EFGcl`gBV zwyrKDI2ZsJV#K$7eL&n);z>Cp;QagxXm{bN1aT5@9@!Fl`uebA3k!?T(0XGrm~iJ9 zi_tSN*${_yB1cC@p{p*ZuHIZw@DROS`RhkE=k}$K>sL~iQn|>S?Il6!{0t2 z@5jc(^dtkXw79jkwW6Y;sVNMQzu-zEBO{Roa3R28K(S|{F^E`JwhWw8*yPaR*pbH= zMntv@^xj=tBVlPe(sF(0{QXaT$rIm&l*yN;s(4@M7EDpIJp+ixnX5y=CRWW z3`ScfvD1Y;F$iwb(vpXZH4qvht;6?r0##b05tmq zqv8@Pp`L)W5Lm26CMIl&ow!$6oEGI6eFy}C>@rrE8QBA>Mmu2|glayGusCM_z(5f% zup3fEMcjB4sL8#?8MkgL?~#4t9&7^14RFIUbpR{|y*fa5?15+)lnA0MXJP2ydgCsF zlf7Oa+4k<&(x9+ESNBhNq`LdXMvDLAva z3*=EnpNXA1S!zTFR9_Jfp1gQ*?*h2Ey^W1$B_${E9-RqnNI~*c-kV9S@bK`c420^_-WQWed_*k{7BnO+Y@udQwMW%!ht$I>I2A!d0)F0H$PNz4q~Og zo}NbNhJdK(WpHW`nY58KbaW81w_s1Erlv$|Qh8-1gu6d5km+=hGe5z?B1AA4y%y}| zCK7WAAo`p!(SQz7%Ldm13HY0?uA8BuNhv9LN`6#=-+*}4D@Et*UfA%Q5%Kg}bh92B z515LM!*$`C3-7oH20Hg(xCVLguK~!y3=ELFHef>#-`v(fX%YkFjFu->T6!Q)J@F7| z0;AQgZbpTG(}KzgijwqK+DN~K2^d5fPpeY3p%AXgRJhwfyIT%)?UO@kBQMtb3aY3y zR#y*BPLg30NQ#bSbx~DNU?50WY2^Qyojnz7NhXs)hBEfalli~9o6=aU7^m%(D<$w| zfpHuVAljA+AQ5-SU{q9{-&kfmbw^LbIyHk#902Xull}4S~?a&goOQaP8*+u+U z;sG0?>xfDpc(wn0R<=Y~SXfO>%~n2ivy7?1A}Fq{o$LRm!2j#-$Q AR(3)) + proclen::Int = 2000 # Process length + # Recurrent net parameters + dev = cpu # Device: cpu or gpu + opt = ADAM # Optimizer + η::Float64 = 1e-3 # Learning rate + hidden_nodes::Int = 10 # Number of hidden nodes + hidden_layers::Int = 2 # Number of hidden layers + layer = RNN # Type of layer, should be one of LSTM, GRU, RNN + epochs::Int = 100 # Number of epochs + seqlen::Int = 49 # Sequence length to use as input + seqshift::Int = 1 # Shift between sequences (see utils.jl) + train_ratio::Float64 = 0.5 # Percentage of data in the train set + verbose::Bool = true # Whether we log the results during training or not + noise_model = Normal(0.0f0, 1.0f0) # Noise to add to the data + K = 5 # Number of simulted observations +end + +# Creates a model according to the pre-defined hyperparameters `args` +function build_model(args) + Chain( + args.layer(1, args.hidden_nodes), + [args.layer(args.hidden_nodes, args.hidden_nodes) for _ ∈ 1:args.hidden_layers-1]..., + Dense(args.hidden_nodes, 1, identity) + ) |> args.dev +end + +nn_model = Chain( + RNN(1 => 32, relu), + Dense(32 => 1, identity) +) + +# Creates training and testing samples according to hyperparameters `args` +function generate_train_test_data(args) + # Generate full AR process + data = generate_process(args.ϕ, args.proclen) + # Create input X and output y (series shifted by 1) + X, y = data[1:end-1], data[2:end] + # Split data into training and testing sets + idx = round(Int, args.train_ratio * length(X)) + Xtrain, Xtest = X[1:idx], X[idx+1:end] + ytrain, ytest = y[1:idx], y[idx+1:end] + + return (Xtrain, Xtest, ytrain, ytest) + # Transform data to time series batches and return + #map(x -> batch_timeseries(x, args.seqlen, args.seqshift) |> args.dev, + # (Xtrain, Xtest, ytrain, ytest)) +end + +function ts_adaptative_block_learning(nn_model, x, y, hparams) + # Warm up recurrent model on first observation + #nn_model(x[1]) + + #@assert length(data) == hparams.samples + + losses = [] + optim = Flux.setup(Flux.Adam(hparams.η), nn_model) + @showprogress for epoch in 1:(hparams.epochs) + Flux.reset!(nn_model) + loss, grads = Flux.withgradient(nn_model) do nn + aₖ = zeros(hparams.K + 1) + for yₜ in y + xₖ = rand(hparams.noise_model, hparams.K) + yₖ = nn(xₖ') + aₖ += generate_aₖ(yₖ, yₜ) + end + scalar_diff(aₖ ./ sum(aₖ)) + end + Flux.update!(optim, nn_model, grads[1]) + push!(losses, loss) + end + return losses +end + +# Trains and outputs the model according to the chosen hyperparameters `args` +function train_model(hparams) + Random.seed!(hparams.seed) + # Create recurrent model + nn_model = build_model(hparams) + # Get data + Xtrain, Xtest, ytrain, ytest = generate_train_test_data(hparams) + + X = [[x] for x ∈ Xtrain] + + loss = ts_adaptative_block_learning(nn_model, X, ytrain, hparams) + return model +end + +cd(@__DIR__) + +hparams = HyperParamsTS() # Set up hyperparameters +m = train_model(args) # Train and output model + +# Generate AR process + +data = CSV.read("/Users/jmfrutos/Desktop/HistoricalQuotes.csv", DataFrame, stripwhitespace=true) +price = tryparse.(Float32, replace.(data.Close, "\$" => "")) +X = price[1:end-1] +Y = price[2:end] + +X = [[x] for x ∈ X] + +epochs = 100 +opt = ADAM() +θ = Flux.params(nn_model) # Keep track of the model parameters +@showprogress for epoch ∈ 1:epochs # Training loop + Flux.reset!(nn_model) # Reset the hidden state of the RNN + # Compute the gradient of the mean squared error loss + ∇ = gradient(θ) do + nn_model(X[1]) # Warm-up the model + sum(Flux.Losses.mse.([nn_model(x)[1] for x ∈ X[2:end]], Y[2:end])) + end + Flux.update!(opt, θ, ∇) # Update the parameters +end + +loss = ts_adaptative_block_learning(nn_model, X, Y, hparams) diff --git a/examples/autoregressive-process/utils.jl b/examples/autoregressive-process/utils.jl new file mode 100644 index 0000000..93bf8d9 --- /dev/null +++ b/examples/autoregressive-process/utils.jl @@ -0,0 +1,41 @@ +# Generates an AR(p) process with coefficients `ϕ`. +# `ϕ` should be provided as a vector and it represents the coefficients of the AR model. +# Hence the order of the generated process is equal to the length of `ϕ`. +# `s` indicates the total length of the series to be generated. +function generate_process(ϕ::AbstractVector{Float32}, s::Int) + s > 0 || error("s must be positive") + # Generate white noise + ϵ = randn(Float32, s) + # Initialize time series + X = zeros(Float32, s) + p = length(ϕ) + X[1] = 10.0f0 + # Reverse the order of the coefficients for multiplication later on + ϕ = reverse(ϕ) + # Fill first p observations + for t ∈ 1:p-1 + X[t+1] = X[1:t]'ϕ[1:t] + ϵ[t+1] + end + # Compute values iteratively + for t ∈ p+1:s + X[t] = X[t-p:t-1]'ϕ + ϵ[t] + end + X +end + +# Create batches of a time series `X` by splitting the series into +# sequences of length `s`. Each new sequence is shifted by `r` steps. +# When s == r, the series is split into non-overlapping batches. +function batch_timeseries(X, s::Int, r::Int) + r > 0 || error("r must be positive") + # If X is passed in format T×1, reshape it + if isa(X, AbstractVector) + X = permutedims(X) + end + T = size(X, 2) + s ≤ T || error("s cannot be longer than the total series") + # Ensure uniform sequence lengths by dropping the first observations until + # the total sequence length matches a multiple of the batchsize + X = X[:, ((T - s) % r)+1:end] + [X[:, t:r:end-s+t] for t ∈ 1:s] # Output +end diff --git a/src/CustomLossFunction.jl b/src/CustomLossFunction.jl index 6e5dcc0..fc94965 100644 --- a/src/CustomLossFunction.jl +++ b/src/CustomLossFunction.jl @@ -4,7 +4,7 @@ Sigmoid function centered at `y`. """ function _sigmoid(ŷ::Matrix{T}, y::T) where {T<:AbstractFloat} - return sigmoid_fast.((y .- ŷ) .* 20) + return sigmoid_fast.((y .- ŷ) .* 10) end; function _leaky_relu(ŷ::Matrix{T}, y::T) where {T<:AbstractFloat} From 2e024cf35ad092d59a942dc299ef4a41fea8dab9 Mon Sep 17 00:00:00 2001 From: josemanuel22 Date: Wed, 13 Sep 2023 10:16:33 +0200 Subject: [PATCH 02/12] Initial time_series_blocklearning --- examples/autoregressive-process/model.jl | 231 ++++++++++++++--------- examples/autoregressive-process/utils.jl | 17 +- 2 files changed, 161 insertions(+), 87 deletions(-) diff --git a/examples/autoregressive-process/model.jl b/examples/autoregressive-process/model.jl index d0495a4..a60c4c2 100644 --- a/examples/autoregressive-process/model.jl +++ b/examples/autoregressive-process/model.jl @@ -7,79 +7,130 @@ using Distributions using DataFrames using CSV using Plots - -include("utils.jl") - -# Hyperparameters and configuration of AR process -@Base.kwdef mutable struct HyperParamsTS - seed::Int = 72 # Random seed - # AR process parameters - ϕ::Vector{Float32} = [.7f0, .2f0, .1f0] # AR coefficients (=> AR(3)) - proclen::Int = 2000 # Process length - # Recurrent net parameters - dev = cpu # Device: cpu or gpu - opt = ADAM # Optimizer - η::Float64 = 1e-3 # Learning rate - hidden_nodes::Int = 10 # Number of hidden nodes - hidden_layers::Int = 2 # Number of hidden layers - layer = RNN # Type of layer, should be one of LSTM, GRU, RNN - epochs::Int = 100 # Number of epochs - seqlen::Int = 49 # Sequence length to use as input - seqshift::Int = 1 # Shift between sequences (see utils.jl) - train_ratio::Float64 = 0.5 # Percentage of data in the train set - verbose::Bool = true # Whether we log the results during training or not - noise_model = Normal(0.0f0, 1.0f0) # Noise to add to the data - K = 5 # Number of simulted observations -end - -# Creates a model according to the pre-defined hyperparameters `args` -function build_model(args) - Chain( - args.layer(1, args.hidden_nodes), - [args.layer(args.hidden_nodes, args.hidden_nodes) for _ ∈ 1:args.hidden_layers-1]..., - Dense(args.hidden_nodes, 1, identity) - ) |> args.dev +include("../../benchmarks/benchmark_utils.jl") + +# AR process parameters +Base.@kwdef mutable struct ARParams + ϕ::Vector{Float32} = [0.4f0, 0.3f0, 0.2f0] # AR coefficients (=> AR(3)) + proclen::Int = 10000 # Process length + x₁::Float32 = 0.0f0 # Initial value + noise = Normal(0.0f0, 1.0f0) # Noise to add to the data + seqshift::Int = 1 # Shift between sequences (see utils.jl) + train_ratio::Float64 = 0.8 # Percentage of data in the train set end -nn_model = Chain( - RNN(1 => 32, relu), - Dense(32 => 1, identity) +# Generates an AR(p) process with coefficients `ϕ`. +# `ϕ` should be provided as a vector and it represents the coefficients of the AR model. +# Hence the order of the generated process is equal to the length of `ϕ`. +# `s` indicates the total length of the series to be generated. +function generate_process( + ϕ::AbstractVector{Float32}, s::Int, x₁::Float32=0.0f0, noise=Normal(0.0f0, 1.0f0) ) + s > 0 || error("s must be positive") + # Generate white noise + ϵ = Float32.(rand(noise, s)) + # Initialize time series + X = zeros(Float32, s) + p = length(ϕ) + X[1] = x₁ + # Reverse the order of the coefficients for multiplication later on + ϕ = reverse(ϕ) + # Fill first p observations + for t in 1:(p - 1) + X[t + 1] = X[1:t]'ϕ[1:t] + ϵ[t + 1] + end + # Compute values iteratively + for t in (p + 1):s + X[t] = X[(t - p):(t - 1)]'ϕ + ϵ[t] + end + return X +end # Creates training and testing samples according to hyperparameters `args` -function generate_train_test_data(args) +function generate_train_test_data(ARParams) # Generate full AR process - data = generate_process(args.ϕ, args.proclen) + data = generate_process(ARParams.ϕ, ARParams.proclen, ARParams.x₁, ARParams.noise) # Create input X and output y (series shifted by 1) - X, y = data[1:end-1], data[2:end] + X, y = data[1:(end - 1)], data[2:end] # Split data into training and testing sets - idx = round(Int, args.train_ratio * length(X)) - Xtrain, Xtest = X[1:idx], X[idx+1:end] - ytrain, ytest = y[1:idx], y[idx+1:end] + idx = round(Int, ARParams.train_ratio * length(X)) + Xtrain, Xtest = X[1:idx], X[(idx + 1):end] + ytrain, ytest = y[1:idx], y[(idx + 1):end] return (Xtrain, Xtest, ytrain, ytest) - # Transform data to time series batches and return - #map(x -> batch_timeseries(x, args.seqlen, args.seqshift) |> args.dev, - # (Xtrain, Xtest, ytrain, ytest)) end -function ts_adaptative_block_learning(nn_model, x, y, hparams) - # Warm up recurrent model on first observation - #nn_model(x[1]) +function generate_batch_train_test_data(hparams, arparams) + Random.seed!(hparams.seed) - #@assert length(data) == hparams.samples + # Get data + Xtrain = [] + Ytrain = [] + Xtest = [] + Ytest = [] + for _ in 1:(hparams.epochs) + xtrain, xtest, ytrain, ytest = generate_train_test_data(arparams) + append!(Xtrain, xtrain) + append!(Ytrain, ytrain) + append!(Xtest, xtest) + append!(Ytest, ytest) + end + loaderXtrain = Flux.DataLoader( + Xtrain; + batchsize=round(Int, arparams.train_ratio * arparams.proclen), + shuffle=false, + partial=false, + ) + loaderYtrain = Flux.DataLoader( + Ytrain; + batchsize=round(Int, arparams.train_ratio * (arparams.proclen - 1)), + shuffle=false, + partial=false, + ) + loaderXtest = Flux.DataLoader( + Xtest; + batchsize=round(Int, (1 - arparams.train_ratio) * (arparams.proclen)), + shuffle=false, + partial=false, + ) + loaderYtest = Flux.DataLoader( + Ytest; + batchsize=round(Int, (1 - arparams.train_ratio) * (arparams.proclen - 1)), + shuffle=false, + partial=false, + ) + + return (loaderXtrain, loaderYtrain, loaderXtest, loaderYtest) +end + +# Hyperparameters for the method `ts_adaptative_block_learning` +Base.@kwdef mutable struct HyperParamsTS + seed::Int = 72 # Random seed + dev = cpu # Device: cpu or gpu + η::Float64 = 1e-3 # Learning rate + epochs::Int = 100 # Number of epochs + noise_model = Normal(0.0f0, 1.0f0) # Noise to add to the data + window_size = 100 # Window size for the histogram + K = 10 # Number of simulted observations +end + +# Train and output the model according to the chosen hyperparameters `args` +function ts_adaptative_block_learning(nn_model, Xₜ, Xₜ₊₁, hparams) losses = [] optim = Flux.setup(Flux.Adam(hparams.η), nn_model) - @showprogress for epoch in 1:(hparams.epochs) + @showprogress for batch in Xₜ₊₁ + j = 0 Flux.reset!(nn_model) + nn_model([Xₜ.data[1]]) # Warm up recurrent model on first observation loss, grads = Flux.withgradient(nn_model) do nn aₖ = zeros(hparams.K + 1) - for yₜ in y + for i in (1:(hparams.window_size)) xₖ = rand(hparams.noise_model, hparams.K) yₖ = nn(xₖ') - aₖ += generate_aₖ(yₖ, yₜ) + aₖ += generate_aₖ(yₖ, batch[j + i]) end + j += hparams.window_size scalar_diff(aₖ ./ sum(aₖ)) end Flux.update!(optim, nn_model, grads[1]) @@ -88,45 +139,57 @@ function ts_adaptative_block_learning(nn_model, x, y, hparams) return losses end -# Trains and outputs the model according to the chosen hyperparameters `args` -function train_model(hparams) - Random.seed!(hparams.seed) - # Create recurrent model - nn_model = build_model(hparams) - # Get data - Xtrain, Xtest, ytrain, ytest = generate_train_test_data(hparams) - - X = [[x] for x ∈ Xtrain] +function ts_mse_learning(nn_model, Xₜ, Xₜ₊₁, hparams) + losses = [] + optim = Flux.setup(Flux.Adam(hparams.η), nn_model) + @showprogress for epoch in (1:(hparams.epochs)) + Flux.reset!(nn_model) # Reset the hidden state of the RNN + loss, grads = Flux.withgradient(nn_model) do nn + nn(Xₜ[1]) # Warm-up the model + sum(Flux.Losses.mse.([nn(x)[1] for x in Xₜ[1:(end - 1)]], Xₜ₊₁[1:end])) + end + Flux.update!(optim, nn_model, grads[1]) + push!(losses, loss) + end + return losses +end - loss = ts_adaptative_block_learning(nn_model, X, ytrain, hparams) - return model +function get_stats(nn_model, data_xₜ, data_xₜ₊₁, hparams) + losses = [] + @showprogress for batch in data_xₜ₊₁ + Flux.reset!(nn_model) + nn_model([data_xₜ.data[1]]) + aₖ = zeros(hparams.K + 1) + for i in 1:(hparams.window_size) + xₖ = rand(hparams.noise_model, hparams.K) + yₖ = nn_model(xₖ') + aₖ += generate_aₖ(yₖ, batch[i]) + end + push!(losses, scalar_diff(aₖ ./ sum(aₖ))) + end + return losses end -cd(@__DIR__) +function plot_ts(nn_model, xₜ, xₜ₊₁, hparams) + Flux.reset!(nn_model) + nn_model([xₜ.data[1]]) + plot(xₜ.data[1:1000]; seriestype=:scatter) + return plot!(vec(-nn_model.([xₜ.data[1:1000]]')...) .+3 ; seriestype=:scatter) +end -hparams = HyperParamsTS() # Set up hyperparameters -m = train_model(args) # Train and output model +@test_experiments "testing" begin + ar_hparams = ARParams(; ϕ = [0.5f0, 0.3f0, 0.2f0], x₁ = rand(Normal(0.0f0, 0.5f0)), proclen=1000, noise=Normal(0.0f0, 0.5f0)) + hparams = HyperParamsTS(; η = 1e-3, epochs = 500, window_size = 10, K = 5) -# Generate AR process + nn_model = Chain(RNN(1 => 32, relu), RNN(32 => 32, relu), Dense(32 => 1, identity)) -data = CSV.read("/Users/jmfrutos/Desktop/HistoricalQuotes.csv", DataFrame, stripwhitespace=true) -price = tryparse.(Float32, replace.(data.Close, "\$" => "")) -X = price[1:end-1] -Y = price[2:end] + loaderXtrain, loaderYtrain, loaderXtest, loaderYtest = generate_batch_train_test_data( + hparams, ar_hparams + ) -X = [[x] for x ∈ X] + loss = ts_adaptative_block_learning( + nn_model, loaderXtrain, loaderYtrain, hparams + ) -epochs = 100 -opt = ADAM() -θ = Flux.params(nn_model) # Keep track of the model parameters -@showprogress for epoch ∈ 1:epochs # Training loop - Flux.reset!(nn_model) # Reset the hidden state of the RNN - # Compute the gradient of the mean squared error loss - ∇ = gradient(θ) do - nn_model(X[1]) # Warm-up the model - sum(Flux.Losses.mse.([nn_model(x)[1] for x ∈ X[2:end]], Y[2:end])) - end - Flux.update!(opt, θ, ∇) # Update the parameters + plot_ts(nn_model, loaderXtrain, loaderYtrain, hparams) end - -loss = ts_adaptative_block_learning(nn_model, X, Y, hparams) diff --git a/examples/autoregressive-process/utils.jl b/examples/autoregressive-process/utils.jl index 93bf8d9..3c4bbdd 100644 --- a/examples/autoregressive-process/utils.jl +++ b/examples/autoregressive-process/utils.jl @@ -1,15 +1,26 @@ +# AR process parameters +Base.@kwdef mutable struct ARParams + ϕ::Vector{Float32} = [0.4f0, 0.3f0, 0.2f0] # AR coefficients (=> AR(3)) + proclen::Int = 10000 # Process length + x₁::Float32 = 0.0f0 # Initial value + noise = Normal(0.0f0, 1.0f0) # Noise to add to the data + seqshift::Int = 1 # Shift between sequences (see utils.jl) + train_ratio::Float64 = 0.8 # Percentage of data in the train set +end + + # Generates an AR(p) process with coefficients `ϕ`. # `ϕ` should be provided as a vector and it represents the coefficients of the AR model. # Hence the order of the generated process is equal to the length of `ϕ`. # `s` indicates the total length of the series to be generated. -function generate_process(ϕ::AbstractVector{Float32}, s::Int) +function generate_process(ϕ::AbstractVector{Float32}, s::Int, x₁::Float32=0.0f0, noise=Normal(0.0f0, 1.0f0)) s > 0 || error("s must be positive") # Generate white noise - ϵ = randn(Float32, s) + ϵ = Float32.(rand(noise, s)) # Initialize time series X = zeros(Float32, s) p = length(ϕ) - X[1] = 10.0f0 + X[1] = x₁ # Reverse the order of the coefficients for multiplication later on ϕ = reverse(ϕ) # Fill first p observations From afc901b6aae1d53a52c575ed54821b35a7de3ee2 Mon Sep 17 00:00:00 2001 From: josemanuel22 Date: Wed, 13 Sep 2023 12:33:33 +0200 Subject: [PATCH 03/12] Initial time_series_blocklearning --- examples/autoregressive-process/model.jl | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/examples/autoregressive-process/model.jl b/examples/autoregressive-process/model.jl index a60c4c2..f27c2ae 100644 --- a/examples/autoregressive-process/model.jl +++ b/examples/autoregressive-process/model.jl @@ -127,7 +127,8 @@ function ts_adaptative_block_learning(nn_model, Xₜ, Xₜ₊₁, hparams) aₖ = zeros(hparams.K + 1) for i in (1:(hparams.window_size)) xₖ = rand(hparams.noise_model, hparams.K) - yₖ = nn(xₖ') + nn_cp = deepcopy(nn) + yₖ = nn_cp(xₖ') aₖ += generate_aₖ(yₖ, batch[j + i]) end j += hparams.window_size @@ -162,7 +163,8 @@ function get_stats(nn_model, data_xₜ, data_xₜ₊₁, hparams) aₖ = zeros(hparams.K + 1) for i in 1:(hparams.window_size) xₖ = rand(hparams.noise_model, hparams.K) - yₖ = nn_model(xₖ') + nn_cp = deepcopy(nn) + yₖ = nn_cp(xₖ') aₖ += generate_aₖ(yₖ, batch[i]) end push!(losses, scalar_diff(aₖ ./ sum(aₖ))) @@ -174,12 +176,18 @@ function plot_ts(nn_model, xₜ, xₜ₊₁, hparams) Flux.reset!(nn_model) nn_model([xₜ.data[1]]) plot(xₜ.data[1:1000]; seriestype=:scatter) - return plot!(vec(-nn_model.([xₜ.data[1:1000]]')...) .+3 ; seriestype=:scatter) + return plot!(-vec(nn_model.([xₜ.data[1:1000]]')...) ; seriestype=:scatter) +end + +function get_density(nn, t) + + + end @test_experiments "testing" begin ar_hparams = ARParams(; ϕ = [0.5f0, 0.3f0, 0.2f0], x₁ = rand(Normal(0.0f0, 0.5f0)), proclen=1000, noise=Normal(0.0f0, 0.5f0)) - hparams = HyperParamsTS(; η = 1e-3, epochs = 500, window_size = 10, K = 5) + hparams = HyperParamsTS(; η = 1e-3, epochs = 500, window_size = 100, K = 5) nn_model = Chain(RNN(1 => 32, relu), RNN(32 => 32, relu), Dense(32 => 1, identity)) From a240fca8cfbe3acd3c570a482cebca712e8493fc Mon Sep 17 00:00:00 2001 From: josemanuel22 Date: Thu, 14 Sep 2023 10:58:17 +0200 Subject: [PATCH 04/12] Initial time_series_blocklearning --- examples/autoregressive-process/model.jl | 55 +++++++++++++++++++----- 1 file changed, 44 insertions(+), 11 deletions(-) diff --git a/examples/autoregressive-process/model.jl b/examples/autoregressive-process/model.jl index f27c2ae..3bff6ac 100644 --- a/examples/autoregressive-process/model.jl +++ b/examples/autoregressive-process/model.jl @@ -131,10 +131,13 @@ function ts_adaptative_block_learning(nn_model, Xₜ, Xₜ₊₁, hparams) yₖ = nn_cp(xₖ') aₖ += generate_aₖ(yₖ, batch[j + i]) end - j += hparams.window_size scalar_diff(aₖ ./ sum(aₖ)) end Flux.update!(optim, nn_model, grads[1]) + for i in (1:(hparams.window_size)) + nn_model([batch[j + i]]) + end + j += hparams.window_size push!(losses, loss) end return losses @@ -146,8 +149,8 @@ function ts_mse_learning(nn_model, Xₜ, Xₜ₊₁, hparams) @showprogress for epoch in (1:(hparams.epochs)) Flux.reset!(nn_model) # Reset the hidden state of the RNN loss, grads = Flux.withgradient(nn_model) do nn - nn(Xₜ[1]) # Warm-up the model - sum(Flux.Losses.mse.([nn(x)[1] for x in Xₜ[1:(end - 1)]], Xₜ₊₁[1:end])) + nn([Xₜ[1]]') # Warm-up the model + sum(Flux.Losses.mse.([nn([x]')[1] for x in Xₜ[1:(end - 1)]], Xₜ₊₁[1:end])) end Flux.update!(optim, nn_model, grads[1]) push!(losses, loss) @@ -172,22 +175,43 @@ function get_stats(nn_model, data_xₜ, data_xₜ₊₁, hparams) return losses end +function get_window_of_Aₖ(transform, model, data, K::Int64) + Flux.reset!(model) + res = [] + for d in data + xₖ = rand(transform, K) + model_cp = deepcopy(model) + yₖ = model_cp(xₖ') + push!(res, yₖ .< d) + end + return [count(x -> x == i, count.(res)) for i in 0:K] +end; + function plot_ts(nn_model, xₜ, xₜ₊₁, hparams) Flux.reset!(nn_model) nn_model([xₜ.data[1]]) - plot(xₜ.data[1:1000]; seriestype=:scatter) - return plot!(-vec(nn_model.([xₜ.data[1:1000]]')...) ; seriestype=:scatter) + plot(xₜ.data[1:800]; seriestype=:scatter) + return plot!(vec(nn_model.([xₜ.data[1:800]]')...) ; seriestype=:scatter) end -function get_density(nn, t) - - - +function get_density(nn, data, t, m) + Flux.reset!(nn) + res = [] + for d in data[1:t-1] + nn([d]) + end + for _ in 1:m + nn_cp = deepcopy(nn) + xₜ = rand(Normal(0.0f0, 1.0f0)) + yₜ = nn_cp([xₜ]) + append!(res, yₜ) + end + return res end @test_experiments "testing" begin - ar_hparams = ARParams(; ϕ = [0.5f0, 0.3f0, 0.2f0], x₁ = rand(Normal(0.0f0, 0.5f0)), proclen=1000, noise=Normal(0.0f0, 0.5f0)) - hparams = HyperParamsTS(; η = 1e-3, epochs = 500, window_size = 100, K = 5) + ar_hparams = ARParams(; ϕ = [0.7f0, 0.1f0, 0.2f0], x₁ = rand(Normal(0.0f0, 0.5f0)), proclen=1000, noise=Normal(0.0f0, 0.5f0)) + hparams = HyperParamsTS(; seed = 50, η = 1e-3, epochs = 800, window_size = 100, K = 5) nn_model = Chain(RNN(1 => 32, relu), RNN(32 => 32, relu), Dense(32 => 1, identity)) @@ -195,9 +219,18 @@ end hparams, ar_hparams ) + loss = ts_mse_learning(nn_model, collect(loaderXtrain)[1], collect(loaderYtrain)[1], hparams) + loss = ts_adaptative_block_learning( nn_model, loaderXtrain, loaderYtrain, hparams ) plot_ts(nn_model, loaderXtrain, loaderYtrain, hparams) + + @gif for i in 1:100 + histogram(get_density(nn_model, collect(loaderXtrain)[1], i, 1000); bins=(-25:0.2:20), normalize=:pdf, label="t=$i") + println("$i") + end every 2 + + bar(get_window_of_Aₖ(Normal(0.0f0, 1.0f0), nn_model, collect(loaderXtrain)[1], 5)) end From b126ccc8049f05f7e055a9a271fd65a658df9e6b Mon Sep 17 00:00:00 2001 From: josemanuel22 Date: Thu, 14 Sep 2023 15:41:20 +0200 Subject: [PATCH 05/12] Initial time_series_blocklearning --- examples/autoregressive-process/model.jl | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/examples/autoregressive-process/model.jl b/examples/autoregressive-process/model.jl index 3bff6ac..905e6bb 100644 --- a/examples/autoregressive-process/model.jl +++ b/examples/autoregressive-process/model.jl @@ -84,19 +84,19 @@ function generate_batch_train_test_data(hparams, arparams) ) loaderYtrain = Flux.DataLoader( Ytrain; - batchsize=round(Int, arparams.train_ratio * (arparams.proclen - 1)), + batchsize=round(Int, arparams.train_ratio * arparams.proclen), shuffle=false, partial=false, ) loaderXtest = Flux.DataLoader( Xtest; - batchsize=round(Int, (1 - arparams.train_ratio) * (arparams.proclen)), + batchsize=round(Int, arparams.train_ratio * arparams.proclen), shuffle=false, partial=false, ) loaderYtest = Flux.DataLoader( Ytest; - batchsize=round(Int, (1 - arparams.train_ratio) * (arparams.proclen - 1)), + batchsize=round(Int, arparams.train_ratio * arparams.proclen), shuffle=false, partial=false, ) @@ -119,7 +119,7 @@ end function ts_adaptative_block_learning(nn_model, Xₜ, Xₜ₊₁, hparams) losses = [] optim = Flux.setup(Flux.Adam(hparams.η), nn_model) - @showprogress for batch in Xₜ₊₁ + @showprogress for (batch_Xₜ, batch_Xₜ₊₁) in zip(Xₜ, Xₜ₊₁) j = 0 Flux.reset!(nn_model) nn_model([Xₜ.data[1]]) # Warm up recurrent model on first observation @@ -129,14 +129,15 @@ function ts_adaptative_block_learning(nn_model, Xₜ, Xₜ₊₁, hparams) xₖ = rand(hparams.noise_model, hparams.K) nn_cp = deepcopy(nn) yₖ = nn_cp(xₖ') - aₖ += generate_aₖ(yₖ, batch[j + i]) + aₖ += generate_aₖ(yₖ, batch_Xₜ₊₁[j + i]) + nn([batch_Xₜ[j + i]]) end scalar_diff(aₖ ./ sum(aₖ)) end Flux.update!(optim, nn_model, grads[1]) - for i in (1:(hparams.window_size)) - nn_model([batch[j + i]]) - end + #for i in (1:(hparams.window_size)) + # nn_model([batch[j + i]]) + #end j += hparams.window_size push!(losses, loss) end @@ -210,8 +211,8 @@ function get_density(nn, data, t, m) end @test_experiments "testing" begin - ar_hparams = ARParams(; ϕ = [0.7f0, 0.1f0, 0.2f0], x₁ = rand(Normal(0.0f0, 0.5f0)), proclen=1000, noise=Normal(0.0f0, 0.5f0)) - hparams = HyperParamsTS(; seed = 50, η = 1e-3, epochs = 800, window_size = 100, K = 5) + ar_hparams = ARParams(; ϕ = [0.5f0, 0.3f0, 0.2f0], x₁ = rand(Normal(0.0f0, 0.5f0)), proclen=2000, noise=Normal(0.0f0, 0.5f0)) + hparams = HyperParamsTS(; seed = 50, η = 1e-3, epochs = 1000, window_size = 100, K = 5) nn_model = Chain(RNN(1 => 32, relu), RNN(32 => 32, relu), Dense(32 => 1, identity)) From d8a70040fb0f8277bd9719c292d0cd80cacd327f Mon Sep 17 00:00:00 2001 From: josemanuel22 Date: Sat, 16 Sep 2023 00:47:30 +0200 Subject: [PATCH 06/12] Initial time_series_blocklearning --- examples/autoregressive-process/model.jl | 81 ++++++++++----- examples/autoregressive-process/utils.jl | 122 +++++++++++++++-------- 2 files changed, 136 insertions(+), 67 deletions(-) diff --git a/examples/autoregressive-process/model.jl b/examples/autoregressive-process/model.jl index 905e6bb..7908e83 100644 --- a/examples/autoregressive-process/model.jl +++ b/examples/autoregressive-process/model.jl @@ -84,19 +84,19 @@ function generate_batch_train_test_data(hparams, arparams) ) loaderYtrain = Flux.DataLoader( Ytrain; - batchsize=round(Int, arparams.train_ratio * arparams.proclen), + batchsize=round(Int, arparams.train_ratio * arparams.proclen - 1), shuffle=false, partial=false, ) loaderXtest = Flux.DataLoader( Xtest; - batchsize=round(Int, arparams.train_ratio * arparams.proclen), + batchsize=round(Int, (1 - arparams.train_ratio) * arparams.proclen), shuffle=false, partial=false, ) loaderYtest = Flux.DataLoader( Ytest; - batchsize=round(Int, arparams.train_ratio * arparams.proclen), + batchsize=round(Int, (1 - arparams.train_ratio) * arparams.proclen - 1), shuffle=false, partial=false, ) @@ -161,16 +161,19 @@ end function get_stats(nn_model, data_xₜ, data_xₜ₊₁, hparams) losses = [] - @showprogress for batch in data_xₜ₊₁ + @showprogress for (batch_xₜ, batch_xₜ₊₁) in zip(data_xₜ, data_xₜ₊₁) + j = 0 Flux.reset!(nn_model) - nn_model([data_xₜ.data[1]]) + nn_model([batch_xₜ[1]]) aₖ = zeros(hparams.K + 1) for i in 1:(hparams.window_size) xₖ = rand(hparams.noise_model, hparams.K) - nn_cp = deepcopy(nn) + nn_cp = deepcopy(nn_model) yₖ = nn_cp(xₖ') - aₖ += generate_aₖ(yₖ, batch[i]) + aₖ += generate_aₖ(yₖ, batch_xₜ₊₁[j + i]) + nn_model([batch_xₜ[j + i]]) end + j += hparams.window_size push!(losses, scalar_diff(aₖ ./ sum(aₖ))) end return losses @@ -188,18 +191,11 @@ function get_window_of_Aₖ(transform, model, data, K::Int64) return [count(x -> x == i, count.(res)) for i in 0:K] end; -function plot_ts(nn_model, xₜ, xₜ₊₁, hparams) - Flux.reset!(nn_model) - nn_model([xₜ.data[1]]) - plot(xₜ.data[1:800]; seriestype=:scatter) - return plot!(vec(nn_model.([xₜ.data[1:800]]')...) ; seriestype=:scatter) -end - function get_density(nn, data, t, m) Flux.reset!(nn) res = [] - for d in data[1:t-1] - nn([d]) + for d in data[1:(t - 1)] + nn([d]) end for _ in 1:m nn_cp = deepcopy(nn) @@ -211,8 +207,13 @@ function get_density(nn, data, t, m) end @test_experiments "testing" begin - ar_hparams = ARParams(; ϕ = [0.5f0, 0.3f0, 0.2f0], x₁ = rand(Normal(0.0f0, 0.5f0)), proclen=2000, noise=Normal(0.0f0, 0.5f0)) - hparams = HyperParamsTS(; seed = 50, η = 1e-3, epochs = 1000, window_size = 100, K = 5) + ar_hparams = ARParams(; + ϕ=[0.5f0, 0.3f0, 0.2f0], + x₁=rand(Normal(0.0f0, 0.5f0)), + proclen=2000, + noise=Normal(0.0f0, 0.5f0), + ) + hparams = HyperParamsTS(; seed=1234, η=1e-3, epochs=1000, window_size=100, K=5) nn_model = Chain(RNN(1 => 32, relu), RNN(32 => 32, relu), Dense(32 => 1, identity)) @@ -220,18 +221,50 @@ end hparams, ar_hparams ) - loss = ts_mse_learning(nn_model, collect(loaderXtrain)[1], collect(loaderYtrain)[1], hparams) - - loss = ts_adaptative_block_learning( - nn_model, loaderXtrain, loaderYtrain, hparams + loss = ts_mse_learning( + nn_model, collect(loaderXtrain)[1], collect(loaderYtrain)[1], hparams ) + loss = ts_adaptative_block_learning(nn_model, loaderXtrain, loaderYtrain, hparams) + plot_ts(nn_model, loaderXtrain, loaderYtrain, hparams) @gif for i in 1:100 - histogram(get_density(nn_model, collect(loaderXtrain)[1], i, 1000); bins=(-25:0.2:20), normalize=:pdf, label="t=$i") + histogram( + get_density(nn_model, collect(loaderXtrain)[1], i, 1000); + bins=(-25:0.2:20), + normalize=:pdf, + label="t=$i", + ) println("$i") end every 2 - bar(get_window_of_Aₖ(Normal(0.0f0, 1.0f0), nn_model, collect(loaderXtrain)[1], 5)) + loss = get_stats(nn_model, loaderXtrain, loaderYtrain, hparams) + + bar(get_window_of_Aₖ(Normal(0.0f0, 1.0f0), nn_model, collect(loaderXtrain)[1], 2)) + + Flux.reset!(nn_model) + for data in collect(loaderXtrain)[2] + nn_model.([[data]]) + end + + prediction = Vector{Float32}() + for data in collect(loaderXtest)[2] + y = nn_model.([[data]]) + append!(prediction, y[1]) + end + + plot(prediction; seriestype=:scatter) + plot!(Float32.(collect(loaderXtest)[1]); seriestype=:scatter) + + ND(Float32.(collect(loaderXtest)[1])[1:200], prediction[1:200]) + + RMSE(Float32.(collect(loaderXtest)[1])[1:200], prediction[1:200]) + + y = collect(loaderYtest)[1] + Flux.reset!(nn_model) + nn_model.([collect(loaderXtest)[1]']) + collect(loaderYtrain)[1] + + get_watson_durbin_test(y, ŷ) end diff --git a/examples/autoregressive-process/utils.jl b/examples/autoregressive-process/utils.jl index 3c4bbdd..91da745 100644 --- a/examples/autoregressive-process/utils.jl +++ b/examples/autoregressive-process/utils.jl @@ -1,52 +1,88 @@ -# AR process parameters -Base.@kwdef mutable struct ARParams - ϕ::Vector{Float32} = [0.4f0, 0.3f0, 0.2f0] # AR coefficients (=> AR(3)) - proclen::Int = 10000 # Process length - x₁::Float32 = 0.0f0 # Initial value - noise = Normal(0.0f0, 1.0f0) # Noise to add to the data - seqshift::Int = 1 # Shift between sequences (see utils.jl) - train_ratio::Float64 = 0.8 # Percentage of data in the train set +using LinearAlgebra +using ToeplitzMatrices + +ND(xₜ, x̂ₜ) = sum(abs.(xₜ .- x̂ₜ)) / sum(abs.(xₜ)) + +function RMSE(xₜ, x̂ₜ) + return sqrt((1 / length(xₜ)) * sum((xₜ .- x̂ₜ) .^ 2)) / ((1 / length(xₜ)) * sum(xₜ)) end +function QLρ(xₜ, x̂ₜ; ρ=0.5) + return 2 * + (sum(abs.(xₜ))^-1) * + sum(ρ .* (xₜ .- x̂ₜ) .* (xₜ .> x̂ₜ) .+ (1 - ρ) .* (x̂ₜ .- xₜ) .* (xₜ .<= x̂ₜ)) +end -# Generates an AR(p) process with coefficients `ϕ`. -# `ϕ` should be provided as a vector and it represents the coefficients of the AR model. -# Hence the order of the generated process is equal to the length of `ϕ`. -# `s` indicates the total length of the series to be generated. -function generate_process(ϕ::AbstractVector{Float32}, s::Int, x₁::Float32=0.0f0, noise=Normal(0.0f0, 1.0f0)) - s > 0 || error("s must be positive") - # Generate white noise - ϵ = Float32.(rand(noise, s)) - # Initialize time series - X = zeros(Float32, s) - p = length(ϕ) - X[1] = x₁ - # Reverse the order of the coefficients for multiplication later on - ϕ = reverse(ϕ) - # Fill first p observations - for t ∈ 1:p-1 - X[t+1] = X[1:t]'ϕ[1:t] + ϵ[t+1] +function get_watson_durbin_test(y, ŷ) + e = [] + for (yₜ, ŷₜ) in zip(y, ŷ) + append!(e, yₜ - ŷₜ) end - # Compute values iteratively - for t ∈ p+1:s - X[t] = X[t-p:t-1]'ϕ + ϵ[t] + sum = 0 + for i in 2:2:length(e) + sum += (e[i] - e[i - 1])^2 end - X + return sum / sum(e .^ 2) end -# Create batches of a time series `X` by splitting the series into -# sequences of length `s`. Each new sequence is shifted by `r` steps. -# When s == r, the series is split into non-overlapping batches. -function batch_timeseries(X, s::Int, r::Int) - r > 0 || error("r must be positive") - # If X is passed in format T×1, reshape it - if isa(X, AbstractVector) - X = permutedims(X) +function yule_(x::Vector{Float64}; + order::Int64=1, + method="adjusted", + df::Union{Nothing,Int64}=nothing, + inv=false, + demean=true +) + method in ("adjusted", "mle") || + throw(ArgumentError("ACF estimation method must be 'adjusted' or 'MLE'")) + + x = copy(x) + if demean + x .-= mean(x) + end + n = isnothing(df) ? length(x) : df + + adj_needed = method == "adjusted" + + if ndims(x) > 1 || size(x, 2) != 1 + throw(ArgumentError("Expecting a vector to estimate AR parameters")) + end + + r = zeros(Float64, order + 1) + r[1] = sum(x .^ 2) / n + for k in 2:order+1 + x[1:(end - k)] + x[k : (end)] + r[k] + println(x[1:(end - k)]) + r[k] = sum(x[1:(end - k)] .* x[k : (end)] ) / (n - k * adj_needed) + end + R = Toeplitz(r[1:end-1], conj(r[1:end-1])) + + rho = 0 + try + rho = R \ r[2:end] + catch err + if occursin("Singular matrix", string(err)) + @warn "Matrix is singular. Using pinv." + rho = pinv(R) * r[2:end] + else + throw(err) + end + end + + sigmasq = r[1] - dot(r[2:end], rho) + sigma = isnan(sigmasq) || sigmasq <= 0 ? NaN : sqrt(sigmasq) + + if inv + return rho, sigma, inv(R) + else + return rho, sigma end - T = size(X, 2) - s ≤ T || error("s cannot be longer than the total series") - # Ensure uniform sequence lengths by dropping the first observations until - # the total sequence length matches a multiple of the batchsize - X = X[:, ((T - s) % r)+1:end] - [X[:, t:r:end-s+t] for t ∈ 1:s] # Output +end + +function plot_ts(nn_model, xₜ, xₜ₊₁, hparams) + Flux.reset!(nn_model) + nn_model([xₜ.data[1]]) + plot(xₜ.data[1:15000]; seriestype=:scatter) + return plot!(vec(nn_model.([xₜ.data[1:15000]]')...); seriestype=:scatter) end From 6b8ccf1f99634fa0b05bdb227fbf87cdbc40c884 Mon Sep 17 00:00:00 2001 From: josemanuel22 Date: Sat, 16 Sep 2023 12:11:11 +0200 Subject: [PATCH 07/12] Initial time_series_blocklearning --- examples/autoregressive-process/test_utils.jl | 112 ++++++++++++++++++ examples/autoregressive-process/utils.jl | 15 +-- 2 files changed, 118 insertions(+), 9 deletions(-) create mode 100644 examples/autoregressive-process/test_utils.jl diff --git a/examples/autoregressive-process/test_utils.jl b/examples/autoregressive-process/test_utils.jl new file mode 100644 index 0000000..49d8038 --- /dev/null +++ b/examples/autoregressive-process/test_utils.jl @@ -0,0 +1,112 @@ +using Test + +tol = 1e-5 + +include("./utils.jl") + +@testset "ND" begin + @test ND([1, 2, 3], [1, 2, 3]) == 0.0 + @test ND([1, 2, 4], [1, 2, 3]) == 1 / 7 + @test ND([1, 2, 3], [1, 2, 4]) == 1 / 6 + @test ND([2, 3, 7], [4, 2, 5]) == 5 / 12 +end; + +@testset "RMSE" begin + @test RMSE([1, 2, 3], [1, 2, 3]) == 0.0 +end; + +@testset "QLρ" begin + @test QLρ([1, 2, 3], [1, 2, 3]) == 0.0 + @test QLρ([1, 2, 3], [1, 2, 3]; ρ=0.9) == 0.0 + @test QLρ([2, 3, 7], [4, 2, 5]; ρ=0.5) == 5 / 12 + x = rand(5) + x̂ = rand(5) + @test QLρ(x, x̂; ρ=0.5) == ND(x, x̂) + + @test QLρ([2, 3, 7], [4, 2, 5]; ρ=0.9) ≈ 29 / 60 +end; + +@testset "yule_walker" begin + rho, sigma = yule_walker([1.0, 2, 3]; order=1) + @test rho == [0.0] + + rho, sigma = yule_walker([1.0, 2, 3]; order=2) + @test rho == [0.0, -1.5] + + x = [0.9901178, -0.74795127, 0.44612542, 1.1362954, -0.04040932] + rho, sigma = yule_walker(x; order=3, method="mle") + @test rho ≈ [-0.9418963, -0.90335955, -0.33267884] + @test sigma ≈ 0.44006365345695164 + + rho, sigma = yule_walker(x; order=3) + @test isapprox(rho, [0.10959317, 0.05242324, 1.06587676], atol=tol) + @test isapprox(sigma, 0.15860522671108127, atol=tol) + + rho, sigma = yule_walker(x; order=5, method="mle") + @test isapprox( + rho, [-1.24209771, -1.56893346, -1.16951484, -0.79844781, -0.27598787], atol=tol + ) + @test isapprox(sigma, 0.3679474002175471, atol=tol) + + x = [ + 0.9901178, + -0.74795127, + 0.44612542, + 1.1362954, + -0.04040932, + 0.28625813, + 0.88901716, + -0.1079814, + -0.33231995, + 0.4607741, + ] + + rho, sigma = yule_walker(x; order=3, method="mle") + @test isapprox( + rho, [-0.4896151627237206, -0.5724647370433921, 0.09083516892540627], atol=tol + ) + @test isapprox(sigma, 0.4249693094713215, atol=tol) + + x = [ + 0.9901178, + -0.74795127, + 0.44612542, + 1.1362954, + -0.04040932, + 0.28625813, + 0.88901716, + -0.1079814, + -0.33231995, + 0.4607741, + 0.7729643, + -1.0998684, + 1.098167, + 1.0105597, + -1.3370227, + 1.239718, + -0.01393661, + -0.4790918, + 1.5009186, + -1.1647809, + ] + + rho, sigma = yule_walker(x; order=3, method="mle") + @test isapprox(rho, [-0.82245705, -0.57029742, 0.12166898], atol=tol) + @test isapprox(sigma, 0.5203501608988023, atol=tol) + + rho, sigma = yule_walker(x; order=3) + @test isapprox(rho, [-0.93458149, -0.68653741, 0.10161722], atol=tol) + @test isapprox(sigma, 0.4269012058667671, atol=tol) + + rho, sigma = yule_walker(x; order=5, method="mle") + @test isapprox( + rho, [-0.83107755, -0.56407764, 0.20950143, 0.1232321, 0.10249279], atol=tol + ) + @test isapprox(sigma, 0.5172269743102993, atol=tol) + + rho, sigma = yule_walker(x; order=5) + @test isapprox( + rho, [-0.96481241, -0.65359486, 0.31587079, 0.28403115, 0.1913565], atol=tol + ) + @test isapprox(sigma, 0.41677565377507053, atol=tol) +end; diff --git a/examples/autoregressive-process/utils.jl b/examples/autoregressive-process/utils.jl index 91da745..b6bad4a 100644 --- a/examples/autoregressive-process/utils.jl +++ b/examples/autoregressive-process/utils.jl @@ -25,12 +25,13 @@ function get_watson_durbin_test(y, ŷ) return sum / sum(e .^ 2) end -function yule_(x::Vector{Float64}; +function yule_walker( + x::Vector{Float64}; order::Int64=1, method="adjusted", df::Union{Nothing,Int64}=nothing, inv=false, - demean=true + demean=true, ) method in ("adjusted", "mle") || throw(ArgumentError("ACF estimation method must be 'adjusted' or 'MLE'")) @@ -49,14 +50,10 @@ function yule_(x::Vector{Float64}; r = zeros(Float64, order + 1) r[1] = sum(x .^ 2) / n - for k in 2:order+1 - x[1:(end - k)] - x[k : (end)] - r[k] - println(x[1:(end - k)]) - r[k] = sum(x[1:(end - k)] .* x[k : (end)] ) / (n - k * adj_needed) + for k = 1:order + r[k+1] = sum(x[1:end-k] .* x[k+1:end]) / (n - k * adj_needed) end - R = Toeplitz(r[1:end-1], conj(r[1:end-1])) + R = Toeplitz(r[1:(end - 1)], conj(r[1:(end - 1)])) rho = 0 try From a9df17a65be7be7aa0bf83b5c0721c617ad38ef6 Mon Sep 17 00:00:00 2001 From: josemanuel22 Date: Tue, 19 Sep 2023 23:28:30 +0200 Subject: [PATCH 08/12] Initial time_series_blocklearning --- examples/autoregressive-process/model.jl | 243 ++++++++++++++++-- examples/autoregressive-process/test_utils.jl | 8 + examples/autoregressive-process/utils.jl | 7 +- src/CustomLossFunction.jl | 2 +- 4 files changed, 240 insertions(+), 20 deletions(-) diff --git a/examples/autoregressive-process/model.jl b/examples/autoregressive-process/model.jl index 7908e83..d8c87fb 100644 --- a/examples/autoregressive-process/model.jl +++ b/examples/autoregressive-process/model.jl @@ -8,6 +8,7 @@ using DataFrames using CSV using Plots include("../../benchmarks/benchmark_utils.jl") +include("utils.jl") # AR process parameters Base.@kwdef mutable struct ARParams @@ -122,7 +123,7 @@ function ts_adaptative_block_learning(nn_model, Xₜ, Xₜ₊₁, hparams) @showprogress for (batch_Xₜ, batch_Xₜ₊₁) in zip(Xₜ, Xₜ₊₁) j = 0 Flux.reset!(nn_model) - nn_model([Xₜ.data[1]]) # Warm up recurrent model on first observation + nn_model([batch_Xₜ[1]]) # Warm up recurrent model on first observation loss, grads = Flux.withgradient(nn_model) do nn aₖ = zeros(hparams.K + 1) for i in (1:(hparams.window_size)) @@ -144,6 +145,78 @@ function ts_adaptative_block_learning(nn_model, Xₜ, Xₜ₊₁, hparams) return losses end + +function ts_covariates_adaptative_block_learning(nn_model, Xₜ, Xₜ₊₁, hparams) + losses = [] + optim = Flux.setup(Flux.Adam(hparams.η), nn_model) + @showprogress for _ in 1:(hparams.epochs) + j = 0 + Flux.reset!(nn_model) + nn_model(Xₜ[1]) + loss, grads = Flux.withgradient(nn_model) do nn + aₖ = zeros(hparams.K + 1) + for i in (1:(hparams.window_size)) + xₖ = rand(hparams.noise_model, hparams.K) + nn_cp = deepcopy(nn) + yₖ = nn_cp.([vcat(x, Xₜ[j+i][2:end]) for x in xₖ]) + aₖ += generate_aₖ(cat(yₖ..., dims=2), Xₜ₊₁[j + i][1]) + nn(Xₜ[j+i]) + end + scalar_diff(aₖ ./ sum(aₖ)) + end + Flux.update!(optim, nn_model, grads[1]) + #for i in (1:(hparams.window_size)) + # nn_model([batch[j + i]]) + #end + j += hparams.window_size + push!(losses, loss) + end + return losses +end + +function ts_covariates_mse_learning(nn_model, Xₜ, Xₜ₊₁, hparams) + losses = [] + optim = Flux.setup(Flux.Adam(hparams.η), nn_model) + @showprogress for epoch in (1:(hparams.epochs)) + Flux.reset!(nn_model) # Reset the hidden state of the RNN + loss, grads = Flux.withgradient(nn_model) do nn + nn(Xₜ[1]) # Warm-up the model + sum(Flux.Losses.mse.([nn(x)[1] for x in Xₜ[1:end]], Xₜ₊₁[1:end])) + end + Flux.update!(optim, nn_model, grads[1]) + push!(losses, loss) + end + return losses +end + +function ts_adaptative_block_learning(nn_model, , tsₖ, s, hparams) + losses = [] + optim = Flux.setup(Flux.Adam(hparams.η), nn_model) + @showprogress for _ in 1:(hparams.epochs) + j = 0 + Flux.reset!(nn_model) + nn_model(insert!([x[1] for x in tsₖ], s, ts₁[1])) + loss, grads = Flux.withgradient(nn_model) do nn + aₖ = zeros(hparams.K + 1) + for i in (1:(hparams.window_size)) + xₖ = rand(hparams.noise_model, hparams.K) + nn_cp = deepcopy(nn) + yₖ = nn_cp(insert!([[x[j + i] for x in tsₖ]], s, xₖ)) + aₖ += generate_aₖ(yₖ, [ts₁[j + i + 1], [x[j + i + 1] for x in tsₖ]]) + nn([ts₁[j + i], [x[j + i] for x in tsₖ]]) + end + scalar_diff(aₖ ./ sum(aₖ)) + end + Flux.update!(optim, nn_model, grads[1]) + #for i in (1:(hparams.window_size)) + # nn_model([batch[j + i]]) + #end + j += hparams.window_size + push!(losses, loss) + end + return losses +end + function ts_mse_learning(nn_model, Xₜ, Xₜ₊₁, hparams) losses = [] optim = Flux.setup(Flux.Adam(hparams.η), nn_model) @@ -208,12 +281,12 @@ end @test_experiments "testing" begin ar_hparams = ARParams(; - ϕ=[0.5f0, 0.3f0, 0.2f0], + ϕ=[0.7f0, 0.2f0, 0.1f0], x₁=rand(Normal(0.0f0, 0.5f0)), - proclen=2000, + proclen=10000, noise=Normal(0.0f0, 0.5f0), ) - hparams = HyperParamsTS(; seed=1234, η=1e-3, epochs=1000, window_size=100, K=5) + hparams = HyperParamsTS(; seed=1234, η=1e-3, epochs=1000, window_size=200, K=5) nn_model = Chain(RNN(1 => 32, relu), RNN(32 => 32, relu), Dense(32 => 1, identity)) @@ -221,14 +294,14 @@ end hparams, ar_hparams ) - loss = ts_mse_learning( - nn_model, collect(loaderXtrain)[1], collect(loaderYtrain)[1], hparams - ) - loss = ts_adaptative_block_learning(nn_model, loaderXtrain, loaderYtrain, hparams) plot_ts(nn_model, loaderXtrain, loaderYtrain, hparams) + loss = ts_mse_learning( + nn_model, collect(loaderXtrain)[1], collect(loaderYtrain)[1], hparams + ) + @gif for i in 1:100 histogram( get_density(nn_model, collect(loaderXtrain)[1], i, 1000); @@ -243,15 +316,20 @@ end bar(get_window_of_Aₖ(Normal(0.0f0, 1.0f0), nn_model, collect(loaderXtrain)[1], 2)) - Flux.reset!(nn_model) - for data in collect(loaderXtrain)[2] - nn_model.([[data]]) - end + res = [] + for i in (1:(hparams.epochs - 1)) + Flux.reset!(nn_model) + for data in collect(loaderXtrain)[i] + nn_model.([[data]]) + end - prediction = Vector{Float32}() - for data in collect(loaderXtest)[2] - y = nn_model.([[data]]) - append!(prediction, y[1]) + prediction = Vector{Float32}() + for data in collect(loaderXtest)[i] + y = nn_model.([[data]]) + append!(prediction, y[1]) + end + r, _ = yule_walker(Float64.(prediction); order=3) + push!(res, r) end plot(prediction; seriestype=:scatter) @@ -261,6 +339,9 @@ end RMSE(Float32.(collect(loaderXtest)[1])[1:200], prediction[1:200]) + yule_walker(Float64.(collect(loaderYtest)[2]); order=3) + yule_walker(Float64.(prediction); order=3) + y = collect(loaderYtest)[1] Flux.reset!(nn_model) nn_model.([collect(loaderXtest)[1]']) @@ -268,3 +349,133 @@ end get_watson_durbin_test(y, ŷ) end + +@test_experiments "testing" begin + csv_file_path = "/Users/jmfrutos/Desktop/LD2011_2014.txt" + + df = CSV.File( + csv_file_path; + delim=';', + header=true, + decimal=',', + types=Dict( + "MT_330" => Float32, + "MT_331" => Float32, + "MT_332" => Float32, + "MT_333" => Float32, + "MT_334" => Float32, + "MT_335" => Float32, + ), + ) + + hparams = HyperParamsTS(; seed=1234, η=1e-3, epochs=1, window_size=1000, K=5) + + nn_model = Chain(RNN(1 => 32, relu), RNN(32 => 32, relu), Dense(32 => 1, identity)) + + loaderXtrain = Flux.DataLoader( + df.MT_334[1:40000]; batchsize=round(Int, 40000), shuffle=false, partial=false + ) + + loaderYtrain = Flux.DataLoader( + df.MT_334[2:40001]; batchsize=round(Int, 40000), shuffle=false, partial=false + ) + + res = Vector{Float32}() + @showprogress for i in 1:200 + loss = ts_adaptative_block_learning(nn_model, loaderXtrain, loaderYtrain, hparams) + append!(res, loss) + end + + Flux.reset!(nn_model) + for data in collect(loaderXtrain)[1] + nn_model.([[data]]) + end + + prediction = Vector{Float32}() + for data in df.MT_334[39900:40100] + y = nn_model.([[data]]) + append!(prediction, y[1]) + end +end + +@test_experiments "testing" begin + csv_file_path = "/Users/jmfrutos/Desktop/LD2011_2014.txt" + + df = CSV.File( + csv_file_path; + delim=';', + header=true, + decimal=',', + types=Dict( + "MT_331" => Float32, + "MT_332" => Float32, + "MT_333" => Float32, + "MT_334" => Float32, + "MT_335" => Float32, + "MT_336" => Float32, + "MT_338" => Float32, + ), + ) + + hparams = HyperParamsTS(; seed=1234, η=1e-3, epochs=100, window_size=100, K=5) + + nn_model = Chain(RNN(5 => 32, relu), RNN(32 => 32, relu), Dense(32 => 1, identity)) + + loaderXtrain = Flux.DataLoader( + [ + [df.MT_333[i], df.MT_334[i], df.MT_335[i], df.MT_336[i], df.MT_338[i]] for + i in 1:length(df.MT_335[1:40000]) + ]; + batchsize=round(Int, 40000), + shuffle=false, + partial=false, + ) + + loaderYtrain = Flux.DataLoader( + [ + [df.MT_333[i], df.MT_334[i], df.MT_335[i], df.MT_336[i], df.MT_338[i]] for + i in 1:length(df.MT_335[1:40000]) + ]; + batchsize=round(Int, 200), + shuffle=false, + partial=false, + ) + + loaderXtest = Flux.DataLoader( + [ + [df.MT_333[40000+i], df.MT_334[40000+i], df.MT_335[40000+i], df.MT_336[40000+i], df.MT_338[40000+i]] for + i in 1:length(df.MT_335[40000:40200]) + ]; + batchsize=round(Int, 40000), + shuffle=false, + partial=false, + ) + + #loss = ts_covariates_mse_learning(nn_model, collect(loaderXtrain)[1], collect(loaderYtrain)[1], hparams) + + loss = ts_covariates_adaptative_block_learning(nn_model, collect(loaderXtrain)[1], collect(loaderYtrain)[1], hparams) + + res = [] + Flux.reset!(nn_model) + for data in collect(loaderXtrain)[1] + y = nn_model.([data]) + append!(prediction, y[1]) + end + + prediction = Vector{Float32}() + for data in collect(loaderXtest)[1] + y = nn_model.([data]) + append!(prediction, y[1]) + end + + prediction = Vector{Float32}() + data = collect(loaderXtest)[1] + y = nn_model(data[1]) + append!(prediction, y[1]) + for i in 1:200 + y = nn_model(vcat(y, data[i][2:end])) + append!(prediction, y[1]) + end + r, _ = yule_walker(Float64.(prediction); order=3) + push!(res, r) +end diff --git a/examples/autoregressive-process/test_utils.jl b/examples/autoregressive-process/test_utils.jl index 49d8038..7a35c38 100644 --- a/examples/autoregressive-process/test_utils.jl +++ b/examples/autoregressive-process/test_utils.jl @@ -9,10 +9,14 @@ include("./utils.jl") @test ND([1, 2, 4], [1, 2, 3]) == 1 / 7 @test ND([1, 2, 3], [1, 2, 4]) == 1 / 6 @test ND([2, 3, 7], [4, 2, 5]) == 5 / 12 + @test ND([2, -3, 7], [4, -2, 5]) == 5 / 12 end; @testset "RMSE" begin @test RMSE([1, 2, 3], [1, 2, 3]) == 0.0 + @test RMSE([1, 2, 3], [2, 3, 4]) == sqrt(3)/2 + @test RMSE([1, 2, 3, 5], [2, 3, 4, 0]) == 8*sqrt(7)/11 + @test RMSE([1, 2, 3, -5], [2, 3, 4, 0]) == 8*sqrt(7)/11 end; @testset "QLρ" begin @@ -23,6 +27,10 @@ end; x̂ = rand(5) @test QLρ(x, x̂; ρ=0.5) == ND(x, x̂) + x = .-rand(5) + x̂ = .-rand(5) + @test QLρ(x, x̂; ρ=0.5) == ND(x, x̂) + @test QLρ([2, 3, 7], [4, 2, 5]; ρ=0.9) ≈ 29 / 60 end; diff --git a/examples/autoregressive-process/utils.jl b/examples/autoregressive-process/utils.jl index b6bad4a..c81d560 100644 --- a/examples/autoregressive-process/utils.jl +++ b/examples/autoregressive-process/utils.jl @@ -4,7 +4,8 @@ using ToeplitzMatrices ND(xₜ, x̂ₜ) = sum(abs.(xₜ .- x̂ₜ)) / sum(abs.(xₜ)) function RMSE(xₜ, x̂ₜ) - return sqrt((1 / length(xₜ)) * sum((xₜ .- x̂ₜ) .^ 2)) / ((1 / length(xₜ)) * sum(xₜ)) + return sqrt((1 / length(xₜ)^2) * sum((xₜ .- x̂ₜ) .^ 2)) / + ((1 / length(xₜ)^2) * sum(abs.(xₜ))) end function QLρ(xₜ, x̂ₜ; ρ=0.5) @@ -50,8 +51,8 @@ function yule_walker( r = zeros(Float64, order + 1) r[1] = sum(x .^ 2) / n - for k = 1:order - r[k+1] = sum(x[1:end-k] .* x[k+1:end]) / (n - k * adj_needed) + for k in 1:order + r[k + 1] = sum(x[1:(end - k)] .* x[(k + 1):end]) / (n - k * adj_needed) end R = Toeplitz(r[1:(end - 1)], conj(r[1:(end - 1)])) diff --git a/src/CustomLossFunction.jl b/src/CustomLossFunction.jl index fc94965..d5cf42c 100644 --- a/src/CustomLossFunction.jl +++ b/src/CustomLossFunction.jl @@ -4,7 +4,7 @@ Sigmoid function centered at `y`. """ function _sigmoid(ŷ::Matrix{T}, y::T) where {T<:AbstractFloat} - return sigmoid_fast.((y .- ŷ) .* 10) + return sigmoid_fast.((y .- ŷ) .* 10.0f0) end; function _leaky_relu(ŷ::Matrix{T}, y::T) where {T<:AbstractFloat} From 51aad9c77c2ce57ab8beef71ae1d45edcb0664ae Mon Sep 17 00:00:00 2001 From: josemanuel22 Date: Tue, 19 Sep 2023 23:52:32 +0200 Subject: [PATCH 09/12] refactor time_series_blocklearning --- GAN/MMD_GAN/mmd_gan_1d.jl | 9 +- benchmarks/benchmark_multimodal.jl | 9 +- benchmarks/benchmark_unimodal.jl | 54 +++++--- benchmarks/benchmark_utils.jl | 4 +- examples/autoregressive-process/model.jl | 97 +------------- examples/autoregressive-process/test_utils.jl | 120 ------------------ examples/autoregressive-process/utils.jl | 86 ------------- 7 files changed, 53 insertions(+), 326 deletions(-) delete mode 100644 examples/autoregressive-process/test_utils.jl delete mode 100644 examples/autoregressive-process/utils.jl diff --git a/GAN/MMD_GAN/mmd_gan_1d.jl b/GAN/MMD_GAN/mmd_gan_1d.jl index 2fbcf75..ee1a3ef 100644 --- a/GAN/MMD_GAN/mmd_gan_1d.jl +++ b/GAN/MMD_GAN/mmd_gan_1d.jl @@ -38,9 +38,9 @@ end epochs::Int = 1000 num_gen::Int = 1 num_enc_dec::Int = 1 - lr_enc::Float64 = 1.0e-10 - lr_dec::Float64 = 1.0e-10 - lr_gen::Float64 = 1.0e-10 + lr_enc::Float64 = 1.0e-3 + lr_dec::Float64 = 1.0e-3 + lr_gen::Float64 = 1.0e-3 lambda_AE::Float64 = 8.0 @@ -73,6 +73,8 @@ function train_mmd_gan_1d(enc, dec, gen, hparams::HyperParamsMMD1D) @showprogress for epoch in 1:(hparams.epochs) for _ in 1:(hparams.num_enc_dec) loss, grads = Flux.withgradient(enc, dec) do enc, dec + Flux.reset!(enc) + Flux.reset!(dec) target = Float32.(rand(hparams.target_model, hparams.batch_size)) noise = Float32.(rand(hparams.noise_model, hparams.batch_size)) encoded_target = enc(target') @@ -93,6 +95,7 @@ function train_mmd_gan_1d(enc, dec, gen, hparams::HyperParamsMMD1D) end for _ in 1:(hparams.num_gen) loss, grads = Flux.withgradient(gen) do gen + Flux.reset!(gen) target = Float32.(rand(hparams.target_model, hparams.batch_size)) noise = Float32.(rand(hparams.noise_model, hparams.batch_size)) encoded_target = enc(target') diff --git a/benchmarks/benchmark_multimodal.jl b/benchmarks/benchmark_multimodal.jl index eb2ef6f..252a604 100644 --- a/benchmarks/benchmark_multimodal.jl +++ b/benchmarks/benchmark_multimodal.jl @@ -13,7 +13,7 @@ include("benchmark_utils.jl") dscr = Chain( Dense(1, 11), elu, Dense(11, 29), elu, Dense(29, 11), elu, Dense(11, 1, σ) ) - target_model = MixtureModel([Normal(4.0f0, 2.0f0), Normal(-2.0f0, 1.0f0)]) + target_model = hparams = HyperParamsVanillaGan(; data_size=100, batch_size=1, @@ -29,7 +29,7 @@ include("benchmark_utils.jl") train_vanilla_gan(dscr, gen, hparams) hparams = HyperParams(; - samples=1000, K=100, epochs=100, η=1e-2, transform=noise_model + samples=1000, K=10, epochs=1000, η=1e-2, transform=noise_model ) #hparams = AutoAdaptativeHyperParams(; # max_k=20, samples=1200, epochs=10000, η=1e-3, transform=noise_model @@ -39,6 +39,7 @@ include("benchmark_utils.jl") #save_gan_model(gen, dscr, hparams) + adaptative_block_learning_1(gen, loader, hparams) ksd = KSD(noise_model, target_model, n_samples, 18:0.1:28) @@ -53,13 +54,13 @@ include("benchmark_utils.jl") #save_gan_model(gen, dscr, hparams) plot_global( - x -> -quantile.(target_model, cdf(noise_model, x)), + x -> quantile.(target_model, cdf(noise_model, x)), noise_model, target_model, gen, n_samples, (-3:0.1:3), - (5:0.2:15), + (:0.2:10), ) #@test js_divergence(hist1.weights, hist2.weights)/hparams.samples < 0.03 diff --git a/benchmarks/benchmark_unimodal.jl b/benchmarks/benchmark_unimodal.jl index 5261085..9b3c6cb 100644 --- a/benchmarks/benchmark_unimodal.jl +++ b/benchmarks/benchmark_unimodal.jl @@ -92,7 +92,7 @@ include("benchmark_utils.jl") dscr = Chain( Dense(1, 11), elu, Dense(11, 29), elu, Dense(29, 11), elu, Dense(11, 1, σ) ) - target_model = Cauchy(23.0f0, 1.0f0) + target_model = MixtureModel([Normal(-10.0, 1.0), Uniform(-5.0,5.0), Pareto(3.0, 10.0)]) hparams = HyperParamsVanillaGan(; data_size=100, batch_size=1, @@ -108,7 +108,7 @@ include("benchmark_utils.jl") train_vanilla_gan(dscr, gen, hparams) hparams = AutoAdaptativeHyperParams(; - max_k=10, samples=1000, epochs=400, η=1e-2, transform=noise_model + max_k=10, samples=1000, epochs=1000, η=1e-2, transform=noise_model ) train_set = Float32.(rand(target_model, hparams.samples)) loader = Flux.DataLoader(train_set; batchsize=-1, shuffle=true, partial=false) @@ -130,7 +130,7 @@ include("benchmark_utils.jl") gen, n_samples, (-3:0.1:3), - (18:0.1:55), + (-20:0.2:30), ) end @@ -313,9 +313,7 @@ end dscr = Chain( Dense(1, 11), elu, Dense(11, 29), elu, Dense(29, 11), elu, Dense(11, 1, σ) ) - target_model = MixtureModel([ - Normal(5.0f0, 2.0f0), Normal(-1.0f0, 1.0f0), Normal(-7.0f0, 0.4f0) - ]) + target_model = Pareto(1.0f0, 2.0f0) hparams = HyperParamsWGAN(; noise_model=noise_model, @@ -323,7 +321,7 @@ end data_size=100, batch_size=1, epochs=1e3, - n_critic=4, + n_critic=2, lr_dscr=1e-2, #lr_gen = 1.4e-2, lr_gen=1e-2, @@ -331,11 +329,11 @@ end loss = train_wgan(dscr, gen, hparams) - hparams = HyperParams(; samples=100, K=10, epochs=2000, η=1e-3, noise_model) + hparams = HyperParams(; samples=100, K=10, epochs=1000, η=1e-3, noise_model) train_set = rand(target_model, hparams.samples) loader = Flux.DataLoader(train_set; batchsize=-1, shuffle=true, partial=false) - adaptative_block_learning(gen, loader, hparams) + auto_adaptative_block_learning(gen, loader, hparams) ksd = KSD(noise_model, target_model, n_samples, 20:0.1:25) mae = min( @@ -347,6 +345,9 @@ end MSE(noise_model, x -> .-x .+ 23, n_sample), ) + save_gan_model(gen, dscr, hparams) + + #@test js_divergence(hist1.weights, hist2.weights)/hparams.samples < 0.03 end @@ -567,6 +568,8 @@ end mse = MSE( noise_model, x -> quantile.(target_model, cdf(noise_model, x)), n_sample ) + + save_adaptative_model(gen, hparams) end @test_experiments "Uniform(-1,1) to Pareto(1,23)" begin @@ -618,23 +621,34 @@ end dec = Chain(Dense(29, 11), elu, Dense(11, 1)) gen = Chain(Dense(1, 7), elu, Dense(7, 13), elu, Dense(13, 7), elu, Dense(7, 1)) - target_model = Normal(23.0f0, 1.0f0) + target_model = Normal(4.0f0, 2.0f0) hparams = HyperParamsMMD1D(; noise_model=noise_model, target_model=target_model, - data_size=1, + data_size=100, batch_size=1, num_gen=1, num_enc_dec=5, - epochs=1e5, - lr_dec=1.0e-2, - lr_enc=1.0e-2, - lr_gen=1.0e-2, + epochs=1000000, + lr_dec=1.0e-3, + lr_enc=1.0e-3, + lr_gen=1.0e-3, ) train_mmd_gan_1d(enc, dec, gen, hparams) + plot_global( + x -> quantile.(target_model, cdf(noise_model, x)), + noise_model, + target_model, + gen, + n_samples, + (-3:0.1:3), + (-5:0.2:10), + ) + + hparams = HyperParams(; samples=100, K=10, epochs=2000, η=1e-3, noise_model) train_set = rand(target_model, hparams.samples) loader = Flux.DataLoader(train_set; batchsize=-1, shuffle=true, partial=false) @@ -651,6 +665,16 @@ end MSE(noise_model, x -> .-x .+ 23, n_sample), ) + plot_global( + x -> quantile.(target_model, cdf(noise_model, x)), + noise_model, + target_model, + gen, + n_samples, + (-3:0.1:3), + (-5:0.2:10), + ) + #@test js_divergence(hist1.weights, hist2.weights)/hparams.samples < 0.03 end diff --git a/benchmarks/benchmark_utils.jl b/benchmarks/benchmark_utils.jl index b1030b1..26cdd0f 100644 --- a/benchmarks/benchmark_utils.jl +++ b/benchmarks/benchmark_utils.jl @@ -88,7 +88,7 @@ function plot_transformation(real_transform, gen, range) linecolor=:redsblues, ) y = gen(range') - return plot!(range, vec(y); legend=:bottomright, label="neural network", linecolor=get(ColorSchemes.rainbow, 0.2), ylims=(-20,20)) + return plot!(range, vec(y); legend=:bottomright, label="neural network", linecolor=get(ColorSchemes.rainbow, 0.2), ylims=(-10,10)) end function plot_global( @@ -149,7 +149,7 @@ function save_gan_model(gen, dscr, hparams) function getName(hparams) gan = gans[typeof(hparams)] lr_gen = hparams.lr_gen - dscr_steps = hparams.dscr_steps + dscr_steps = hparams.n_critic noise_model = replace(strip(string(hparams.noise_model)), "\n" => "", r"(K = .*)" => "", r"components\[.*\] " => "", r"prior = " => "", "μ=" => "", "σ=" => "", r"\{Float.*\}" => "") target_model = replace(strip(string(hparams.target_model)), "\n" => "", r"(K = .*)" => "", r"components\[.*\] " => "", r"prior = " => "", "μ=" => "", "σ=" => "", r"\{Float.*\}" => "") basename = "$gan-$noise_model-$target_model-lr_gen=$lr_gen-dscr_steps=$dscr_steps" diff --git a/examples/autoregressive-process/model.jl b/examples/autoregressive-process/model.jl index d8c87fb..d0fae86 100644 --- a/examples/autoregressive-process/model.jl +++ b/examples/autoregressive-process/model.jl @@ -8,102 +8,7 @@ using DataFrames using CSV using Plots include("../../benchmarks/benchmark_utils.jl") -include("utils.jl") - -# AR process parameters -Base.@kwdef mutable struct ARParams - ϕ::Vector{Float32} = [0.4f0, 0.3f0, 0.2f0] # AR coefficients (=> AR(3)) - proclen::Int = 10000 # Process length - x₁::Float32 = 0.0f0 # Initial value - noise = Normal(0.0f0, 1.0f0) # Noise to add to the data - seqshift::Int = 1 # Shift between sequences (see utils.jl) - train_ratio::Float64 = 0.8 # Percentage of data in the train set -end - -# Generates an AR(p) process with coefficients `ϕ`. -# `ϕ` should be provided as a vector and it represents the coefficients of the AR model. -# Hence the order of the generated process is equal to the length of `ϕ`. -# `s` indicates the total length of the series to be generated. -function generate_process( - ϕ::AbstractVector{Float32}, s::Int, x₁::Float32=0.0f0, noise=Normal(0.0f0, 1.0f0) -) - s > 0 || error("s must be positive") - # Generate white noise - ϵ = Float32.(rand(noise, s)) - # Initialize time series - X = zeros(Float32, s) - p = length(ϕ) - X[1] = x₁ - # Reverse the order of the coefficients for multiplication later on - ϕ = reverse(ϕ) - # Fill first p observations - for t in 1:(p - 1) - X[t + 1] = X[1:t]'ϕ[1:t] + ϵ[t + 1] - end - # Compute values iteratively - for t in (p + 1):s - X[t] = X[(t - p):(t - 1)]'ϕ + ϵ[t] - end - return X -end - -# Creates training and testing samples according to hyperparameters `args` -function generate_train_test_data(ARParams) - # Generate full AR process - data = generate_process(ARParams.ϕ, ARParams.proclen, ARParams.x₁, ARParams.noise) - # Create input X and output y (series shifted by 1) - X, y = data[1:(end - 1)], data[2:end] - # Split data into training and testing sets - idx = round(Int, ARParams.train_ratio * length(X)) - Xtrain, Xtest = X[1:idx], X[(idx + 1):end] - ytrain, ytest = y[1:idx], y[(idx + 1):end] - - return (Xtrain, Xtest, ytrain, ytest) -end - -function generate_batch_train_test_data(hparams, arparams) - Random.seed!(hparams.seed) - - # Get data - Xtrain = [] - Ytrain = [] - Xtest = [] - Ytest = [] - for _ in 1:(hparams.epochs) - xtrain, xtest, ytrain, ytest = generate_train_test_data(arparams) - append!(Xtrain, xtrain) - append!(Ytrain, ytrain) - append!(Xtest, xtest) - append!(Ytest, ytest) - end - - loaderXtrain = Flux.DataLoader( - Xtrain; - batchsize=round(Int, arparams.train_ratio * arparams.proclen), - shuffle=false, - partial=false, - ) - loaderYtrain = Flux.DataLoader( - Ytrain; - batchsize=round(Int, arparams.train_ratio * arparams.proclen - 1), - shuffle=false, - partial=false, - ) - loaderXtest = Flux.DataLoader( - Xtest; - batchsize=round(Int, (1 - arparams.train_ratio) * arparams.proclen), - shuffle=false, - partial=false, - ) - loaderYtest = Flux.DataLoader( - Ytest; - batchsize=round(Int, (1 - arparams.train_ratio) * arparams.proclen - 1), - shuffle=false, - partial=false, - ) - - return (loaderXtrain, loaderYtrain, loaderXtest, loaderYtest) -end +include("ts_utils.jl") # Hyperparameters for the method `ts_adaptative_block_learning` Base.@kwdef mutable struct HyperParamsTS diff --git a/examples/autoregressive-process/test_utils.jl b/examples/autoregressive-process/test_utils.jl deleted file mode 100644 index 7a35c38..0000000 --- a/examples/autoregressive-process/test_utils.jl +++ /dev/null @@ -1,120 +0,0 @@ -using Test - -tol = 1e-5 - -include("./utils.jl") - -@testset "ND" begin - @test ND([1, 2, 3], [1, 2, 3]) == 0.0 - @test ND([1, 2, 4], [1, 2, 3]) == 1 / 7 - @test ND([1, 2, 3], [1, 2, 4]) == 1 / 6 - @test ND([2, 3, 7], [4, 2, 5]) == 5 / 12 - @test ND([2, -3, 7], [4, -2, 5]) == 5 / 12 -end; - -@testset "RMSE" begin - @test RMSE([1, 2, 3], [1, 2, 3]) == 0.0 - @test RMSE([1, 2, 3], [2, 3, 4]) == sqrt(3)/2 - @test RMSE([1, 2, 3, 5], [2, 3, 4, 0]) == 8*sqrt(7)/11 - @test RMSE([1, 2, 3, -5], [2, 3, 4, 0]) == 8*sqrt(7)/11 -end; - -@testset "QLρ" begin - @test QLρ([1, 2, 3], [1, 2, 3]) == 0.0 - @test QLρ([1, 2, 3], [1, 2, 3]; ρ=0.9) == 0.0 - @test QLρ([2, 3, 7], [4, 2, 5]; ρ=0.5) == 5 / 12 - x = rand(5) - x̂ = rand(5) - @test QLρ(x, x̂; ρ=0.5) == ND(x, x̂) - - x = .-rand(5) - x̂ = .-rand(5) - @test QLρ(x, x̂; ρ=0.5) == ND(x, x̂) - - @test QLρ([2, 3, 7], [4, 2, 5]; ρ=0.9) ≈ 29 / 60 -end; - -@testset "yule_walker" begin - rho, sigma = yule_walker([1.0, 2, 3]; order=1) - @test rho == [0.0] - - rho, sigma = yule_walker([1.0, 2, 3]; order=2) - @test rho == [0.0, -1.5] - - x = [0.9901178, -0.74795127, 0.44612542, 1.1362954, -0.04040932] - rho, sigma = yule_walker(x; order=3, method="mle") - @test rho ≈ [-0.9418963, -0.90335955, -0.33267884] - @test sigma ≈ 0.44006365345695164 - - rho, sigma = yule_walker(x; order=3) - @test isapprox(rho, [0.10959317, 0.05242324, 1.06587676], atol=tol) - @test isapprox(sigma, 0.15860522671108127, atol=tol) - - rho, sigma = yule_walker(x; order=5, method="mle") - @test isapprox( - rho, [-1.24209771, -1.56893346, -1.16951484, -0.79844781, -0.27598787], atol=tol - ) - @test isapprox(sigma, 0.3679474002175471, atol=tol) - - x = [ - 0.9901178, - -0.74795127, - 0.44612542, - 1.1362954, - -0.04040932, - 0.28625813, - 0.88901716, - -0.1079814, - -0.33231995, - 0.4607741, - ] - - rho, sigma = yule_walker(x; order=3, method="mle") - @test isapprox( - rho, [-0.4896151627237206, -0.5724647370433921, 0.09083516892540627], atol=tol - ) - @test isapprox(sigma, 0.4249693094713215, atol=tol) - - x = [ - 0.9901178, - -0.74795127, - 0.44612542, - 1.1362954, - -0.04040932, - 0.28625813, - 0.88901716, - -0.1079814, - -0.33231995, - 0.4607741, - 0.7729643, - -1.0998684, - 1.098167, - 1.0105597, - -1.3370227, - 1.239718, - -0.01393661, - -0.4790918, - 1.5009186, - -1.1647809, - ] - - rho, sigma = yule_walker(x; order=3, method="mle") - @test isapprox(rho, [-0.82245705, -0.57029742, 0.12166898], atol=tol) - @test isapprox(sigma, 0.5203501608988023, atol=tol) - - rho, sigma = yule_walker(x; order=3) - @test isapprox(rho, [-0.93458149, -0.68653741, 0.10161722], atol=tol) - @test isapprox(sigma, 0.4269012058667671, atol=tol) - - rho, sigma = yule_walker(x; order=5, method="mle") - @test isapprox( - rho, [-0.83107755, -0.56407764, 0.20950143, 0.1232321, 0.10249279], atol=tol - ) - @test isapprox(sigma, 0.5172269743102993, atol=tol) - - rho, sigma = yule_walker(x; order=5) - @test isapprox( - rho, [-0.96481241, -0.65359486, 0.31587079, 0.28403115, 0.1913565], atol=tol - ) - @test isapprox(sigma, 0.41677565377507053, atol=tol) -end; diff --git a/examples/autoregressive-process/utils.jl b/examples/autoregressive-process/utils.jl deleted file mode 100644 index c81d560..0000000 --- a/examples/autoregressive-process/utils.jl +++ /dev/null @@ -1,86 +0,0 @@ -using LinearAlgebra -using ToeplitzMatrices - -ND(xₜ, x̂ₜ) = sum(abs.(xₜ .- x̂ₜ)) / sum(abs.(xₜ)) - -function RMSE(xₜ, x̂ₜ) - return sqrt((1 / length(xₜ)^2) * sum((xₜ .- x̂ₜ) .^ 2)) / - ((1 / length(xₜ)^2) * sum(abs.(xₜ))) -end - -function QLρ(xₜ, x̂ₜ; ρ=0.5) - return 2 * - (sum(abs.(xₜ))^-1) * - sum(ρ .* (xₜ .- x̂ₜ) .* (xₜ .> x̂ₜ) .+ (1 - ρ) .* (x̂ₜ .- xₜ) .* (xₜ .<= x̂ₜ)) -end - -function get_watson_durbin_test(y, ŷ) - e = [] - for (yₜ, ŷₜ) in zip(y, ŷ) - append!(e, yₜ - ŷₜ) - end - sum = 0 - for i in 2:2:length(e) - sum += (e[i] - e[i - 1])^2 - end - return sum / sum(e .^ 2) -end - -function yule_walker( - x::Vector{Float64}; - order::Int64=1, - method="adjusted", - df::Union{Nothing,Int64}=nothing, - inv=false, - demean=true, -) - method in ("adjusted", "mle") || - throw(ArgumentError("ACF estimation method must be 'adjusted' or 'MLE'")) - - x = copy(x) - if demean - x .-= mean(x) - end - n = isnothing(df) ? length(x) : df - - adj_needed = method == "adjusted" - - if ndims(x) > 1 || size(x, 2) != 1 - throw(ArgumentError("Expecting a vector to estimate AR parameters")) - end - - r = zeros(Float64, order + 1) - r[1] = sum(x .^ 2) / n - for k in 1:order - r[k + 1] = sum(x[1:(end - k)] .* x[(k + 1):end]) / (n - k * adj_needed) - end - R = Toeplitz(r[1:(end - 1)], conj(r[1:(end - 1)])) - - rho = 0 - try - rho = R \ r[2:end] - catch err - if occursin("Singular matrix", string(err)) - @warn "Matrix is singular. Using pinv." - rho = pinv(R) * r[2:end] - else - throw(err) - end - end - - sigmasq = r[1] - dot(r[2:end], rho) - sigma = isnan(sigmasq) || sigmasq <= 0 ? NaN : sqrt(sigmasq) - - if inv - return rho, sigma, inv(R) - else - return rho, sigma - end -end - -function plot_ts(nn_model, xₜ, xₜ₊₁, hparams) - Flux.reset!(nn_model) - nn_model([xₜ.data[1]]) - plot(xₜ.data[1:15000]; seriestype=:scatter) - return plot!(vec(nn_model.([xₜ.data[1:15000]]')...); seriestype=:scatter) -end From 8d58f0a3bb9a6dfa8fdfc6d6597679638706cad4 Mon Sep 17 00:00:00 2001 From: josemanuel22 Date: Wed, 20 Sep 2023 00:38:55 +0200 Subject: [PATCH 10/12] refactor time_series_blocklearning --- examples/autoregressive-process/Project.toml | 8 - examples/autoregressive-process/README.md | 18 - examples/autoregressive-process/loss.png | Bin 22347 -> 0 bytes examples/autoregressive-process/model.jl | 386 ------------------- src/CustomLossFunction.jl | 158 +++++++- 5 files changed, 148 insertions(+), 422 deletions(-) delete mode 100644 examples/autoregressive-process/Project.toml delete mode 100644 examples/autoregressive-process/README.md delete mode 100644 examples/autoregressive-process/loss.png delete mode 100644 examples/autoregressive-process/model.jl diff --git a/examples/autoregressive-process/Project.toml b/examples/autoregressive-process/Project.toml deleted file mode 100644 index 53b9832..0000000 --- a/examples/autoregressive-process/Project.toml +++ /dev/null @@ -1,8 +0,0 @@ -[deps] -Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c" -Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" -Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" - -[compat] -Flux = "0.13.9" -julia = "1.6" \ No newline at end of file diff --git a/examples/autoregressive-process/README.md b/examples/autoregressive-process/README.md deleted file mode 100644 index 91bbf1b..0000000 --- a/examples/autoregressive-process/README.md +++ /dev/null @@ -1,18 +0,0 @@ -# Autoregressive Model - -An [autoregressive (AR) process](https://en.wikipedia.org/wiki/Autoregressive_model) is a stochastic process with an autoregressive structure, i.e., past realizations influence its future realizations. - -This model-zoo example illustrates how to use Flux's recurrent layers to model an AR process. - -The example contains the following files: -+ [utils.jl](utils.jl): - + `generate_process`: generates an AR process - + `batch_timeseries`: transforms a vector into the proper format for recurrent layers in Flux and allows to batch the time series as required. - -+ [model.jl](model.jl): creates and trains the recurrent model to predict the generated AR process. - -## Example loss - -Running the model with the hyperparameters currently given in the example, we obtain the following train and test losses. We see that the model begins to overfit after around 30 epochs. - -![loss](loss.png) \ No newline at end of file diff --git a/examples/autoregressive-process/loss.png b/examples/autoregressive-process/loss.png deleted file mode 100644 index a2c97a5d226db75117daeb6d7dd4d92d7c39f678..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 22347 zcmb4rbyQW|7wu6&E+rsHr*tFT2uLd>jescK-6g4nAl*oJcb9ZYH%NDPzs>g>@BjC% zcL3wYKKtyw_KG>@@`t>v1S%3f5(Gi0Qj+f!A?T?Z1U(6ZJq54ass0!Ozn*`Qk$49^ zJpTLLm>mH@q>$7*VI`-;-Fe5)m}?IRha0u^{aLPZ8Z)fY)2>m<=N7_&-HLyh18KNh z7z1?>eliELMpBY~e?MI=U$9Fie{O*{&nFc@&@8l#29JalzdEvdcbzX)r+UDcPyUg3!^t(ttZqw?795S7f;7fr->;mV--eUEozgBdry_TvG(}KfI&nL zh2tw(X>?!kZeMaaJ{%bYd7Y8Mo znz)FF_ps5J5xWh}Yl8CW#E1ynl}=RC#JOSbF!bdad!zO2p%^Xx+`v3$M@9(3CWW;z zH{Xtw=!(`N|N3BHkagzaCP6?>;VbVq(zqDwt-OedaU=W;Y9Z(TjQ&!DJ!GbQe}U4| zBDx+deu?8xJYQe*izm!r9bZRd_h`W&Xo@2MRYdGO_Nhar*?`*aZu5wa%EQ9pk2MVjr#ASyn)(cj;~ z^78V98kN1x+O;-o?alc)Ih*wdP30f_LOY(|j4k>o)H{W^@fq&l^oB|9)<1{j?LF_# zRIhOCcP5f}?YBS53D|E?`n0yN1WRYqJ!quLTg@#BV-ga&?a$SErbzs^+l4Bo`EB&; zmy9AJ-&qTyi$ot9C`W(L^lBkPUe9LA^nsaa`M)5-$HTiCSJffnw)!ob#Cx_giH3+M zB`r;Nz%kI>{r&rQ8lfiZmCn(u5BVRnW!2QwLWu>4rwdJ#e=fde*uYMSbYCi7PNM%* zPt^kkeycoB6E%o6wW(2fm)n=oQA8wnvV^Ow3BNed zA_t|I-jM@>lx<#zy1Dm?*Gt}yCEE`VOAuq@;qh4gdn*EKXL<7ENp5bg=k1Y%hey+} z5x?8z;aykr>4=PUW;pst?#JwM!y%1^-Co;xQVL)Di}UyKvGt$h4=$Ox7eCkbJ%OBT z1VoL865-7KAsei`el zuVcz9wC=5&hx zBvj8=VHA;6QW|3$;+9=!F+Vq_?Xs15{OtZDX}QUMxjA$m+WatZ4K&G~VY1J@#@c%{fQYc!TKoF*1z^>tyk+VX_ok~4t= zKLlQGKbW!2X8$f_!|kc9$1Qrb$y3P6V35GqSFt!-M^&uTz=LPO?!#bQJ#fwB8~bBl zaI1-jXhhR(Jma_Vh2Bl-ddNtLHivW3!K?HB?Fa8 z_VKLS5_7a$>7jzl&i5mkGK5|3M;$oDP1oCb%J~NUv2R>XH)%g+p^Mnq+9Hx6V>7@| zPztPe6!6Gk2kj+H<}Stw_DlEHmc8|LoTv|s(YtE^P)_Y zVo~B(zVxa^i7dKp`}6gPWXZ|NJv}|XVfY*-kJBb4{cSjue1Ct>Npsdi3ZoZ~GO4x@ z)v8Bf;&W8;WXLoceJy5{ktM9DAHo~hVwG2NQFsEil;)-pfM<_ZhfrC0xm+^8^Ug%! zalqHEpRdHMdn0Md$jBZK(N<>6%a<>)88s*@f1!Dky5aWPY{z1QXZW-pzw4jas^-%c z6+yP7lNU;G5_I0CHpw0Z0bDyveKpPs*v;=DE5$c&-hjO&;jvj27z!l>zq)^2U0o?C zC{W8L!qCug#1ybPsHa}^yBscP(JG~mjae4vF2i|D90lm5L==+M6iKnxe^A1USPQ14XgTnF=pI@=VWo`tb@D7Agucy$UNo&l z0|O##sj*)9^5qMR?C9tS98Dc?^@T=PA8+qo$7ze{vNeBH0#3^XLST*I;mccFQn2Yt z)5)1iQ)ZneGuAl0qWcFj7^AP^#C(GNMH3e(BAKrw6?uy5Y@&5fq=&XALxDot}DStQ2U975A?p%TW&@cAA0Vj=y} zrPDuytgs!}zm3N!MB5sV(*H|x?hVT+z>7HcOd`Wt@72PBygt4~pTWpes_#TYbaEwL zS&lwib86S|?AF47yjDcfXNbh}Z4cOHTzALXP?3+LWn+3#Y*gW(f87PF;zq1-3$6!^ zBO_dbgX+=%1ntAfA?i%`;JDp#+iYk{U z&%Ju$l~`rD=y8ART$>*@aOdb)Mw$;tHgtm09!$W+L`%CxCPPIjh0o_uP+#u~GA)kY zoc$b;tzRsI=<_U=OYB+!aFa|vYeg_<^$0BcPy*M`!oH#D*CBz6#Hc~DDAF;L$`p@l zS~6e7>#%l@%cTI9vuy=w=j3W}yr06r#jqAQ42F^%)>iW)G`_=C*$;$TChCB*Y;<+3 z&9_o33&GF^xh6zSAMRquaDABgN0O!Y%NX|Ar($>Gfg*|d-xbz=odqC#!b}EF^G(h8~h5@ zO4CY8SSTpKQ6^VcSNA4O6m+k1+nBgt4EvVVjV)#=5oDu%hF*u=UJ-D1-%_tr-uROrGB5fb-Fc56_Kmk?oZgo zs8K<~#pO7XA@#qg#iu zO<7%^Zuy3p&sKqSV{L5>tjWBso{>=wDNKc~UBMW?iu_$~+}+Qing_b4Cf!;uptl?b zg`{&8#-n#~vWeVB_tc9{E2#A2|EAbcXE7CekdQhBd!Bw@8iV**F39RWJS!FrJ=4jUjDdI@VY8Ojt>)w=QEDPP1mqqLcP z3d!G*!ic_^ZV2b`ez4v65Y?}?gbW=zq5HmQOMf^i$hKSc ze|C6;ccOb=%+}9U-2Iqqa)1%u28PiLwKA{WpQZxG_*gTXzYyMx1@~*3U_Tgl@#{~* zS`j=5>0~gXXw4yvf&jVuK+DIp|G$@TGPOmi<19ac$fvOns>od&@4vqkGG`3?5h&|b zqLYM}9*1=3ed2RiCA^Qci@3A3!SO|Ue<m$eC>*ZG-bgLP-Jj zQ$5OH(!(2v*YHqHN;= z{DX}J=DCu;c6SYr3FH(;7YZ1whW;*{*^ym!yW-bkA`=-ZB3Ytk$wnArEpqmcOtHsh z{>)6s3Ca^~y70Lm*HosHbKTS>6OSXY;oVdB%=l~bijZaC-l#9lAm}a3dKV1|7xQW2c8>b*hTvuQpPFTuvh3ma6ZJ(+ zXusM>!qmC`$QqQ9i)Oz2;j4ggFm>UM@OZ(T^33ygmrytQot+qE413_aSBsuQZw&@t z`mW)-j&?41cFliiw4;#V{UyPfA)}Ls&Wt5R##HeOXOHk*{qs>tE;`Fxl@6CnC#ny< zcUQI#lZLEJ7B(@@AX4fLq&83Fk&)jlGhDj|0S%H&*ir#D8AMzep`kzMqW>U_?zdK` z(xhAe(Fy0xs;#pAv@o(IIOO>*hdf&lYH8=W#MNE1M#)_f4v=^g#9P9Zyv}K7xa0hH z+OXrF{`xR&+h9WY3S)cAXRtrQ=z@~L2=0{8(+d$YGD_@+X_mxTA~8pngy@4fyz<1X z1sC6YD0}kXs}A*QAwgbd=)R%oIh2@R$6OaG$mrO(;v2_f$r z%cG67Ql)(j2wMF{T2N3&d^)*!T z=^U<%Qo|gvO=WW5wR#?`G#EJhv?ef9kuoO_Y5jIKj^eidGa6G^8{IH%x92u1HD>L> zY+1df)KcF7a~4wealX!tJO=vug6+HgHaV*y5y*-YCepWwE2(OyD@RH9w4!1Io6O=_ z^KQ5rRkzl}ik0!KCWd^v)rv?h@)M|p7A}<%-lXRyB#vIR=0|D&?a7xZo09|tKjd-4 z4!^E{F{sq{EXh0!kEfx^hqQwKC)u<1xIKZ|u_i*wjJxN89+=I`HXLf!tw#3RL+=N5 z&MW03|3eTsy-6$5Xt=JOkJn4Va@FJ~9n%O_yQQ1AC z5c|T-c9tW*!*xH)UjNCuwM4s}8pmAyOOvzMqgN?-A#7xd|1&tRgK%@`ZYEA3k&=n0 z{MR&TR_$*8EgBPCYJZ#EX%O~q6or^5^NcgOkX{T?Vb=vN|AN)A+)AA2Barb~3?@?l zFOM{WfHlUSCjMrmDB)6=Z}0UBybo6#sX30#%2_yAj#hO0)dh7=-FI9?2+^ z&!G}3INh_ih!VLG4g&5oI=$c}tOSzP!YChJ#gIAUU~mb z_1)TxPM(XMc^?pR8~kx%J}QJRvPiC3p<&DLBt&MS;l+P5BSNLrGb37mD`YXiQF`Zc zulVcLXxJH|4G{~c0mkR)KDlb**`m+t~LqMzL=)OQZuhN3;h$*ABUHTP^Ga0P`3$RQah`&#?99I zPrmfV74dSl=1~@kdGaoyF7;M-Pah9y9ywg<%E!{x^7Y-&y=mPZ?X*iP2db>>ij@?q{6% zxr89QbjFDz*izJY4FDBFn00;S^dtOA^7_IH=?$Txz*1RiOH8lq;v$CDVWnDz|iFvUB?BzZWZ8KG=NJVB1Jx-zv`2 zJx}!wQ!sQcMN7u!-SLtFF1o9aW8B*f7fbr_z|B1qhWVnuT$28HW33vPc|UvgvRr>K z7o~V-;bJ;5MGz~*r7sZ!YC%HxJ>){8oE}R8(e!e_t;jDc+(gq(S!`lRhc`)5rCI1T z)Dj_#-X_HD-)}y24olseJi)`BaqlEQ5-Y8!G3v4jN<@WPmeGB!aHP8fK1|_kvse@O zBKu)?i}c!Pr!UZcnAdyHCx2EB43p3so$gSY;Y1fRYm{oqsR+I-%>Lw*_)KE>Cc_& zB%w(2Zwz+E^B=Y9tryxJUMl2VPeTv`yf6%Hlr55JyX9=YlceGuwhM~WHx{$jx|BQ* zItY>hH${F5OGW0hx~o!8)9Pl?CeMmY;?sz)8-AGyCi}KS*lBAd6QD%Vo{X@zM}QFE z2?tF#vgm(v*MePq|H4~;S~S01O3X7LOP_EN5~?)X87V{VVpJr6-o8>M^mV^EzuT)` z>}^gKbiaOthz=G6j7KsiwOt4NE58wK;Yk}->fczb$*n4LAoCnR#R`@ijQ5E{wZ^|@ zASf@*o0J_F#c(hY%mSJ0U3dnlD0D#?1YibDQ0ap*MOx;~&X5!;ilIS;fvfo|CLC1` zWMz)8^DZbuS4YoEE4p7-K7sg4k-S-igi_F#4-VvDQZh0!a&iD6p&-l9Y4!lH_V##w zRU!Ir_eYMn3w$;$F2w4*2G<;!1nu~#2_D;pi1*ly94t`F*E~{KV?zT?CV=aF!^ZPf zP1)lo3N=!n>{b~=MRK_S5(*m&$BTO_ZdcvU34>_m*ZdM|rqH>k;_|Vu`S-F2A!(ci6&eU?r4Q^L908g&A zm$9Dur6oeuX}Jf3a`jYsHvmHdB>vTN zT`8#_vf1^R^)=eMyEh^E*Ed{?b@6OFjG?W}rrXt}lIBre z8w+$&+L=Oi63X?ZMOr$_dUp~YL3*5ca(m1I743Fs6(G!LhAtH(WTQ&oalxeqqx-5T zJ6vI|zuK9xBVrNj-oT+(^JwllxUG2tlFc79a>6H~Brbkr%+lX|$n zLnNE8cN78#V=VV0z>kxE{jxq@m6E2P#XJwgULwL_`ed}wa23s{sZp7CxmT=OtW6a` z%g^6%s|!e;=PwAIHk2!1`g$z5SV$e2ibXpP2D+NLwx5+JagMo|uo+{s2u{n&WIqNe zU80L<8}#JGy&S3+qF^M~#vXb2bL= zVXoWddh4b3avL!x(H|$LY>!#RYJ~%9VT`rKKEuUTHl|1mjny~%A1A24^Lqrz=PC)Z zlG1f=+}_@Xhlc|;3ta>W895>%0vH94$ym+{BHpnAbvANxa(4EL@K^jJ0|{KNXWM7P z5+njN@Q7p#ed6NcoF=CiUEluxm58)w=1RL@S(A@*_u<_ z!JNtql7hwh!IBHd;s-Q-5tyzzC(Wl3l;weQ2}#=|Loq^Vu$c<5_byn@qp zCMf0CaM+(S=1i!!-w}r)zIZ|9t7l-KK?h@FV_R8SdE|?9+k68Evy|xS9S`n*zY*~h zNtLGoQ`_xe{M9XTm8cW>vsprL#YvCDi&jbY2UF^KLK``0Nj4~(^jyH|RuXdJP zA9Atwu3=jnGAa*%@U?L+Q9qB-5k_}$b#--dDFVR+Kw!u?Oa;J!IajZS5OBc&dmwio zrxTv>4B-W!Mc?L+8yFgPb#4-{}Y z;5L+MJbJy)|NLJqK=S<*@qyu@s{dr;2}Sy_vmi(6?s2ZSa=8tTTYpYPa1cSdJtYHq z+<*sWlbfV!VAt;>bvCN#jDnBVTjp0bnSt*UI785vcXni8SNAsueJQtF*(q~#bNvRk zev3At0J-bv>XML`{sp2PibE}%A4 z3N_AttAYp&FP!sU>a(>=$jJ}>#EEN<(6t&FU8YK-7=8OiVF(*jo5p&{O;%$Zf^FNwXvjw zh&^KEHOvN^)E)K`kLDCw<*&O(<83V#+_j?I)zrUXkVbiV- z0TuJhmG3-V3d6~P7>D>byP56_#IFm+dfz-TN-Jm_Tf2g1I!rJ1@?4IcJHm4r9;y%e zxtN?oJz`E!i(GW0>P@a9$8`8Hb;j6#mS1kYQ5P4FWjI04+S#yK_xN|Q&6Ow@wfkL@il4r@cekO^XK z#y@y-C^Vso{;9%{5jE7hP`78Go0!{+EN~ext3OGE03$|>i;GJ@Fqz7!XxdcA6&;PXGfU|T~$2BKCnk)7o}Lf=ll1IMRxSJ z`73&7*1&&~&67=iI_DqGgQY;Wyq;Gj`>dr%N^v24Yew8!Wgc5}wki}r>bpGj9&$ZH zLkEE8F)^VddzUI7sZ;53cTK=$QSomE(Ml<#acj3~w`H&G+r3#>Ri3&eY`Ei0Rg8<_ z_AWbIz;ym5c!b^PqMHNFzxkdfo2jOic%U)JR*@z`)%zAVnU)Lv_aB&PYd=;OPzHN0^-iPO@0p>*%ucCcgT9k*m;2T=K7&ProoSyRuRB2?qt3J%$J)5l*w>{ill%Sd0C)5&V@3D^YE%>)vpMwJh1 zr`Oijf<*3~C((SSg3NcW!FjY(7&eYYsV7W0yX| z=4$_VR$rM11zHY%MfLwmdzJ5Th<_emXK~`HJyNyV3qRJa#WLWw60Uin>$YDQK5u#| z;h6tb^Q!4GH_c6_!X#%I{){e)j)Q}PlkI^?}Veci&Pi;@MI-Bld+W9jNrzq--nIuYp#3%7bbc?OETlLKntr5Zf$TxP3b zZ^VtNi3PqH@Ju2c zwdtxC@N{J#5q5ga$rLFpOW<*wB zb{|v9%j0sOjd1BMojWg`I_+!-r^w8Be5`J4)YsNlnvCaZR$G*h%@d`xEYJltN!-HUB5Dlj0{tYSrmz=*j4Y%3eSGoSqIo|`Xq%7Huh#X!Jvd#) z1ywP#5y~792_+40?|Ul0yiA-vpleS9m91(V@73HaPx-?7t?0OYLar%S>HNC4ia%)q z1Y&mRD_`F*TR+9Q=N`6>MN(0YiPlQ`F|ybN$sI;^#E2`zf(UX8q>;TS{+PL$0S>X2 zpC;JNvdt&dUnvZovwfzFC^hN}AMxRvy9u^2o}R~s|1L7v5c)%LZ^V_|lTO+5bAlez zXw_!t6K01z1t6QRKMFyGTDL2e9?$XEjs zumRHqs2;1j$v2Zy&wgsK^(zxMuihV0;x00) zbD?rvkl06j80OqRycEi4EVn=hS{|E~&XCYsf1o~6J3Nm1Mp~gM2blOU)eGd9zzlmW zp63H?2$*M-Vwy3VXa_dN}R&2^fz>^-alFi0esLZYRWv(J>%y4Au6Gk{auXR7_9k?X|iy+|X6 zBjCxH!({0|9JC*!Z^WUgmdTi zH+HB%*!xakInw6^y6jAKKM}!&7;5U}hbJ=qF+H!FZsF4%oXe8*DavV2`i}M+{Y~-j z1)eit>ERq`fOy$uxGT|H($$pY4JJ!GrQ)zu!?Thmub6}K69kM9!xt_>H{Vzr23L{` z>)xWPK*BWX)PP6s{?It~C5c7o=N@*j%2GO_Uu_4U$2iG+q4CsUrZbMd6H4lshA>Uw zdI8F4WqIdlKra*FiktV?8^*CS4H{8K=ShN{aEg_!XC3vEhw6BFInR4bW`Obri|kkE zu-||BBl$2!tTkl6+Ar3`Demo&9>b1pF+lp1+lFnLm2YUf;khd0Tm^Wth0->_LdB@e z+Y@2-%4xnzYxu58# z1sea#En|xkIkosy9rhSfqKaTSt(*L~A;`+6c$LWM*0rWKArxRu-1Zp-xjuH+4;Xz? zlk=Py_c=S%{UiXAL;EBHM^3n$UhQYtb8xq_enyu_ii8G;*hinD48*A2 zlhhqL%t5s7+{zqk4SEvqw{T7a9JfPv4L*H}A{!12lhl{Yc+I^{W|^DT^p1fMEAYhc z)-07qj;0#}sqAV_Q)G~_tt3&76sV_QKb8E&X+_d-L!LEp$QVRib?TQ^$g^*Yea+Py z3%U+>{i!?Oz?}k3@BxC%?b2j04=z=Yb1Rg}l%RuWvMD$EYx;-kCuUz2-3>gW*{z#D zFU@}B=G8ZlBo}*KwkBvFk#t~l_^kA z`E7J4RQ}EKe2>xxFGgcB_ImAtGK^wap6-aM=Mz=jKzMDiqBI9naOc^A;9IDMz#Dxv z<73#k)^rW3o2cc?6oKW0@>^Q*1Sq**i=CM%WVPS!5lOq_dkY{x%zu26|zy*CgjQIM0b%+=bip|!WQ0l*79 z4hYe)*SlS7y^~_cqs^L!a#(;dM<+@;=!^bE!++#Yxi+J6A5r>`YMp z?vs|f!yg2@3Xud_bH7|os69(zn1)bG9dT(#JShx_x&>@jzsqLJNK5DGG&h;|KHbEB zzi9mDe3oi=r>jJcE|H4tITJ6b&>!X0fd2| zq^A;%`xmKQ)Ro8-7VdCFP-n|fC;6m>cYhmStK8rNK+CrtER`O0+5qSfxH}zHE&3$& z;R7lfS_-#y8ljl?$#oTn>B}Cz0Qg73I^awWJuuB_sm%|-9P`cwXWQdUOiYiGcr-N2 z(X0W(6n6(8c)GoHrY2V zZB0$2kl9EEA{m8%zeb#-*z;9;(I{nM*o#K;|6H2M; z^#n^8S)=PklWhxNCTKHLQ&U+Z5oL=jNs>sqDEI#jVUZ#~X$jwt=eR!VA_08D?Be3$ zzyKo`SGDD0)7pIMS~RnccC&}OqoZR24&qxFnS0x24CC~Ef?uQiG(K-XPzv?c(8B-_ z7@e5N4{%;@Z*PE2%f+*A4a#I@W{Qc40g$ND*$;Za#qPk0m7E-CTM$fa7T-64TFB<= ze>VY8ZEvQs{(&0@Ck7O!(-{lBxD<2_Fd*kO%B}VF@55j%XkVUk7;S8^Crn=Dn3g4LCohLl6mzjU*78mHRy}S z)*s5EGi?>B_Nv~1>T1{qFdm=7&r$>Rh6=+HaV&9NcQ%e-Rd_r6lhX%NN znLg?dlZO(ve~q6&q1nso>EFJ6TN`M-ySd;po!AB76^Wo*O$hz79DX5Gb)uoBMc=K; zs1zGXmULqn6m?#*h6xOLy21Gb90xez?>~Q5segdv46NX6jUYz`36tvNhZW9ve@LG8 z)m-YrXBM-m5?#>Az~{J+f{5sR4In_z`>XD5$ypqP!x|B-N1i?HWO);-RZeultmPE- zEAD$*crS!{dKj=2KFfs$<$M(kk(Zd500UoHUf!R}#`&%B+ZV`1g1;N6=)3}wg(iDR zxAZ3wN)f|Khq-9bzkO0dK{qEA6%`SYw%@`&cDzZ`W(z=|D6oi5LSor|3?-wB_WXt! z019*`>vUXdY2008gEY$%b%r|lNTE-Q5Sof`6`~+O@y~p01F1U|Ui8ZJm2b3rk(T6B zbIV#rtM$q8LwwERTuYDxvhMOzBU_z3gS8?wXw^X#{gFw#-LHpUbbb``7kks5CK=bV zUq`!FN84mA|GqB3TU zof4rY9eMa|%LVthh`EEpIj(%Yyu-=~jWo7*KoIp|(>iYIZ#7kP27RZLc zQ{7hCZA*@%_^QDYCbf0oY_MD|onr8I(59OjwQ*E-~u;W}nw0KlREJMt{oB z2`}0?p*x78a%n2*qr8rn{f?N+<7rM6xTIXSs8m+9-otK}`5<~i;o3w0pB6BnFDYD3 zI1N`_Cc=d^+xnKC?jH*>0#~f3N?b4z2zisBGf{H2VLaoW?qsg~^@JLvL=Il>v-%ij zsZl!XViUw1I52s1Qwov%j(>~(R88?2bNonvYbJ3-iI0RF?b}K(Qn(Tp>F)1F>!+lC z0fB#tSUuep*^z7};UF}B@E#Hx0|zJrg}nO{FxZlYaZODnwE@ET|O6AI>*hpROy7}Ey$prM$Bx>u{kg7~ zq`HFKm0dOOD&P52o>d(*APQaSj+#RWCR-Ak*QM2{yNPUHZm9&v{9c6UCFe=`R7^D! zyd^fmc-l?+R#ecaEdx#Yy=iB1cf@w_GN0wDColE-=URBLfmAUVDM562QEP=MizZb` z!brzjOy(q7Y=YD1x21kBEBtM5&evXypet@Q`%6!YmBDB3H%PnZ1k3UHBNAP{FjZcb zL_B$)iLSi`}~=sKEf5eMIa9MW8PLq!}Y^%sn}SF-dk=WH6P<835*{QHZaQ) z2?oj}AEOOY2XRvyb_wt?!gAGa#)xbGH+Lp1THaybA#5M;G1JJfEWwIgv>WS>HnnOR zjvubKeOVaPw{creUr@T6guN}CqkO7l3ls4&NzbBB!j+-!WOdeGz39_Vuw=C);tbd_ z8j_!(k2GGBP_MvzLQLX~{`>k!@mDHGPYF$%Rp6caNS58vnvkC*ZXJF4x)4MWR13JhYG{{WI!Q_i6a@M#)xt#4=e~_shF%i zJ$6aqUF6i)+L36Ktw!>;xTm=?=%e(J=Scc2@i2L^IG>W58Z|&|E73XkK~@wfj43*z zApcm!?TJRQ{#tjsE^I3x$4lcag_7O6$qj#rV`?AYaC_Vb_@MLSHf{wiM=}${fNZe0{wN^df)DFPbzUY1;N57w8Gix`121kEg(6bSMf>-8rv}Fo|52+!2o!SY*x1;h1#zaCx^&B~4x|YPM4uEDkCy^SJOI(Yh9(h$e>P)UJnyVr zJbA)D#1c!!^vtLAw>%B^3)@hAc@>pTmVQxQExMtguph&#Q8eL+Jq&e!n(V)yzrdD| z2xN5sYpC(&v9oK$r_;S9&ncNJVZfVI>DZ}I<)+&!Y=g0Xaxw;8QB7@gd07NTOG9&b za^m{8)xyG}++aZZKRJTG?wAROw>&mO-)q=fnmi4|D__S`N6Nc&$*NzsZo-6~YI-o| zrWb6l8(*!w z?vE3m+jZC-CZs)9Jt@&6yNagh(ZcnLQrjuNU1|>iOMU8cceK(4)FeUJjDU($Df(0k zn0^w^+h6Fc7IQUdf|sOpblYuy$RK~^bGs~wqI=>+Ju2u`h4%TRZz?n9&yO45m&yE^ zL3qkR&yC+?KPL>V^tS2=rv!o#&_47(amo;omw?tfCAz5meABBZpsj_DiHVq1?|VQ% zSDAiaT6+4UohkWc2L^}gc6PA@*k(9y{JGBet-1#S7~=zZ@oyN2RxqWqn-`&f{`>)b z0s=n4+A_;oocb5~Qp5E9 zcU6eZ+SsTKb5T$Dsq9Riv_yObopBxhYj3J%A2oxDSStX$%3&hS>%cS zJCa=iUt7~nu-sk1WxzPB<#30>tyf; zkS+hf=>Kk%6(S^-LOX>E`n2_I4U4uy6peetjCZwLDZ-bA`@hrR?+3nxBM_dXfx4Ne z{A7XR6LmS--oQ3kG05JnE9BPusR^?zw-|%jx{ffdn}Kb@n|+>^+C%&tG?I;oMhYVdkXdX(R_99TigD=hJ$X_ z(y1`F-$*>SBB(D1oat0cwbsy(06li>Y7q({49qp}+qL#-U$hwzHEW_y?{Q^(Ut3S@ zX^TdYdK*Klg!lLf)Sr)4hNae2pD5Pe-PTyT$R_KFDS~hFaah3ND#$M6@0{|wgWmKv zdCcuEBu0}zt3pOqJ<&UtOyMu`NdnVoY5By(eYcLnFFK~B@v~Tp{oJTyC?4H!UUm6D2`N#c<8yM>JXuwT5 zbM~4fh~Pj9zWVx{)h_^K242qs?b1i%OD7|+fFwQ>r{k~w`i~pG3O9E2)Th7>tkeia zpZl{VxLO5V+w#J1m5Ldi@Ch9dI#~DsEaFMakzf5tXhNB4M2#asdBR905XXN#NAjGD z{wo6K%{T$IIG048B4sv;Ouv6tf0=gKM&d0+L7pCmG0Z-1Np1(a%r_tZUx4=kNwbMK z8eqD1w;VtV{eFgYKSp-f+-J;=k)-*@!`jRBF>YG0CfvpAO?CzZoL_U6{p$LMaQ=Y# zs-7cKm&%0zA>>o9@{gr@R&eE+2#!wBMk+;PHBzC+2Izq}1BdS02QXvkUPV+Qn^V7- z!uo*Az_XdjkBLog>Y_yXszGu>3)9tdE#-r7wh4!S4=V&Y*0N<&HAcp{ssu0onY76k z5RHYYzTB5?w;%JNoRsMq%J#k5$cVRCy_)F{WS`h93Pqd$k@pBryP10vDY|-izW>-UXhzE+L*l@6<>&T*8ojDPatCP)-Ycy zDze)Y!{am1CRBufKFFqEaSEiV8Cq>iPZDKyZ}3y4|ArwzPW()S3JVttJnZ4dByYQR zJ3i~4fusg12?iEph@9cKVU3()@?0w&NHmL-&g7@%i1@+=5)~>z zT)qk&Fl8txcS9<>T5pBxWM6u_4I|0BhE7tW>GYHNdOH(6sz41Wg%;1{6@aC#D4tiVcTnO+v+O*Ief{}uznJNwA69@IE&54%eX zvRo?Ql#JJ9#&^@3G0Y%FFe-rko|AlMWn9$RQ+iwVi6Wx-bsN ze@SG6C=F_nZMZk_;h~)8K*;@tf7900uO|Rb`1L(Leh!%j{`D7>QAXUA0A+8mzd6Kd z^TyI7KheMZlgbjV(P^18PnON$vW8%D($OfuwjQH!n z*d7(NXy*Ge?RL0=ZtUQzUgY1cGV5J!uk=^tf7-FgTwkTJ-^-;|A>1{gY7SvqU$Q(d zO4m~E>v?bED}KcvB|*&JLFtuk!U|kX5mIeXxhgI>9Ax!AG7MuFoOM41QhU{7A)OYF zg;$^;fw)$cVO+S@ z0Yx?R4~yHMpQwkBp!XUcFB*hMegr(H7OcNpgezy9v$I#m;gqjEa%NUmW+5aI z*;I*`#yeu{C@YBNB5EU`~7*nU-S8VhI*Fx z){e%{GQ)x(Eyf??;%Uz@{I(LBw(xxidSjoBa)cPsU;rCT94TR+?x~-OkW$DH^V)HA z03`_CFcJ-z9AcgAY0>mc_tSaSu%WO=VlmRETW{Hl%fnUxVSWVh*v0d{8m4n~y?<@F zWAbB0qg4Oqfa*+c9>ugakPtVTl!CPteGZhXcU3h087g@F-+D4XzKLC5>ucujLwyld zd+wl=Gbi8|2lf>4wlE5n>hhj()0l?zJ0cd!= zUvH8C22Cni+d!67$)E!Vgit#c0|nkzoy_V3tpo(a)tZ>YnK-kzz+?jIE*0WT=8jj7 zzR8z4>ISLK;_d5yC!P;qbC;m(6e}C-3InXMhs4BL}V36!sm3Zq2yaDeTx$zW~nSekh5vTDOa-@VM6_<46?7hkr1{qF$Bs^89kCouvCY@}LEd{v{Fi*cx!v-BCa;q>?C zfc~s$!2_nrp;!k$JDvvp4a>O?E&SVS2d9PDCHm68?cN?!UwpPNBC_umGLR>H&TZ{hq^tSdlH3SMxd&4=(}He}U9rW9?$ok^N6|eYdj7 z&B)4I9(&u2lrP(oX0Bej-%2Zm>U@|tQuW-RK?%0|Kd$-Ml^s*m5e3+G0 z2mI-`Z{OC}*H<%a|CkQnS{h7CV-5-V1EmXyWF!c|vYhMhZ@e?DeDiS5)^GB(inm}@?VaO?RsoL7RPX)o6$Z*LNjNHi1ILHa52Xlk~WmzOV&?dq>9~&Ek9-jsJ82lGS?JkcwIZiKFos=nmj|8CxlS&mOC2`h-#W6sh#zaMJ z?1p6;)?OWY3T2XW6G{EnHdbAo47H_zHB2o4q1a=eY4~JlnG?12TRXC#K#)DL)9v++ zjyxw9*CiY$mOXL(jPc;ol6cCTJ?TFDOwcPbCk`uV47LU$|=Bs0CI4m<@P>(ox6`ui~s!F`MxX+5daA05>k(ZYj zDz`lbm_;R&j>H+VOhnM$c zUDdsN_ndwiSw?hKz&jz4Oj<>hlsG${ z!7wd2Bm^UhWtYO?(o&X^GOUSY@}KE>{~KmzX$Q{^=6E-yBqhOCY$|Mj!fffBO9<_R z*Ne(Z8F_gtEPLIH7p221PIq=`Xd_ELZLf`2U1*qHT>LyVWa;GO{IUtg9RM#? z4pZAp%m+&D%RSKmVLpc2`kZSugZ?`rS2e zZW)wW`8QlB3~O?|aPGkOEYx))!MJEt>r#?WE7&K zR|ylVk|A^-2p+(N7P$|6gxstweWwdzf*=Ql#jZAhs>0&euYT(j$~@hfhN7w6uL2es zc?Tb)fmjHPve13sKUWRBdsnQ;R|IZbMO8K4?ixI2$A?DCZD)rc9^96aB2G;?@^s5y zZ2Y#m3KSH_pwyU~sb<7s$Hu_;?8*3EH$7clHa50F-g|^cn&BEZ=KFhl&xY5pEi5eb zWyv0bOon{4I}32H?BQc$Q!?lgmrlpT2R$d8kPwyNn|Xg&H^VfQq57S_AYv=A^z@0<+G1C0$I{bTer5A>ah(ktOH53prlyWBo6*+R2486Z3Oc}&G44|{AqwMdL+EhCR{%S)|r?{tg6ySDyyh$|Ngxk zs8U!C3Md~j7b3$0_qW&JeF3UCL-d=bt^o5<8k+W#`4+Hm@c6(vh{N-ApUk31{e!Nt zt*D#2`5iQ|-gw&CapgptnLM$G@0HN!9Ws>aAafzCbfgS1j0C|teBW2 zL#rB~+jyYz+<*X8U-}OTjP*oJ#l5_FA&pvd&^o{4i@vnAw(X> z$91}ZYM6k8#MQ<3xnM^DK|#^^k*r*Nd>Xg5D-Nhxii(OL7XT5ZOV7XPR3-OrjnT>e zI~QtbcwCNwoTq+EiRV3ZO<@WdF9K>U;8^q!=jrHTq32pQyaE>b*6kEeHyA1p z22*snC6z516chwvCbw_j#uwnk#Pkkl2OvL@l#V0x~j;d9HbQAj^ql$)SW_Ei`XrYpSV1B%hp`TL1Awvw?p$ zCYO>GMuhj8{XkPeEU={de<2}~l8v?%YeA)-ytZ~& z_6lxh425`cR|nY^I1W?~O-7T0ou95e-D|6>UDgna#WM|MvijVgcaN+D0wK{67$4xy zs%twsIwUg_y62Kf+}zwZ_qNKKo7GY2bjA#S8EKrKip$6>{0tQp6@_F4La-X2tcXZ) zQBmYI91i#FavwWfcHo*gIE;f2ZE4Y{f!xPRhyW?w-80VK-fOe7g_@T@-EFGcl`gBV zwyrKDI2ZsJV#K$7eL&n);z>Cp;QagxXm{bN1aT5@9@!Fl`uebA3k!?T(0XGrm~iJ9 zi_tSN*${_yB1cC@p{p*ZuHIZw@DROS`RhkE=k}$K>sL~iQn|>S?Il6!{0t2 z@5jc(^dtkXw79jkwW6Y;sVNMQzu-zEBO{Roa3R28K(S|{F^E`JwhWw8*yPaR*pbH= zMntv@^xj=tBVlPe(sF(0{QXaT$rIm&l*yN;s(4@M7EDpIJp+ixnX5y=CRWW z3`ScfvD1Y;F$iwb(vpXZH4qvht;6?r0##b05tmq zqv8@Pp`L)W5Lm26CMIl&ow!$6oEGI6eFy}C>@rrE8QBA>Mmu2|glayGusCM_z(5f% zup3fEMcjB4sL8#?8MkgL?~#4t9&7^14RFIUbpR{|y*fa5?15+)lnA0MXJP2ydgCsF zlf7Oa+4k<&(x9+ESNBhNq`LdXMvDLAva z3*=EnpNXA1S!zTFR9_Jfp1gQ*?*h2Ey^W1$B_${E9-RqnNI~*c-kV9S@bK`c420^_-WQWed_*k{7BnO+Y@udQwMW%!ht$I>I2A!d0)F0H$PNz4q~Og zo}NbNhJdK(WpHW`nY58KbaW81w_s1Erlv$|Qh8-1gu6d5km+=hGe5z?B1AA4y%y}| zCK7WAAo`p!(SQz7%Ldm13HY0?uA8BuNhv9LN`6#=-+*}4D@Et*UfA%Q5%Kg}bh92B z515LM!*$`C3-7oH20Hg(xCVLguK~!y3=ELFHef>#-`v(fX%YkFjFu->T6!Q)J@F7| z0;AQgZbpTG(}KzgijwqK+DN~K2^d5fPpeY3p%AXgRJhwfyIT%)?UO@kBQMtb3aY3y zR#y*BPLg30NQ#bSbx~DNU?50WY2^Qyojnz7NhXs)hBEfalli~9o6=aU7^m%(D<$w| zfpHuVAljA+AQ5-SU{q9{-&kfmbw^LbIyHk#902Xull}4S~?a&goOQaP8*+u+U z;sG0?>xfDpc(wn0R<=Y~SXfO>%~n2ivy7?1A}Fq{o$LRm!2j#-$Q x == i, count.(res)) for i in 0:K] -end; - -function get_density(nn, data, t, m) - Flux.reset!(nn) - res = [] - for d in data[1:(t - 1)] - nn([d]) - end - for _ in 1:m - nn_cp = deepcopy(nn) - xₜ = rand(Normal(0.0f0, 1.0f0)) - yₜ = nn_cp([xₜ]) - append!(res, yₜ) - end - return res -end - -@test_experiments "testing" begin - ar_hparams = ARParams(; - ϕ=[0.7f0, 0.2f0, 0.1f0], - x₁=rand(Normal(0.0f0, 0.5f0)), - proclen=10000, - noise=Normal(0.0f0, 0.5f0), - ) - hparams = HyperParamsTS(; seed=1234, η=1e-3, epochs=1000, window_size=200, K=5) - - nn_model = Chain(RNN(1 => 32, relu), RNN(32 => 32, relu), Dense(32 => 1, identity)) - - loaderXtrain, loaderYtrain, loaderXtest, loaderYtest = generate_batch_train_test_data( - hparams, ar_hparams - ) - - loss = ts_adaptative_block_learning(nn_model, loaderXtrain, loaderYtrain, hparams) - - plot_ts(nn_model, loaderXtrain, loaderYtrain, hparams) - - loss = ts_mse_learning( - nn_model, collect(loaderXtrain)[1], collect(loaderYtrain)[1], hparams - ) - - @gif for i in 1:100 - histogram( - get_density(nn_model, collect(loaderXtrain)[1], i, 1000); - bins=(-25:0.2:20), - normalize=:pdf, - label="t=$i", - ) - println("$i") - end every 2 - - loss = get_stats(nn_model, loaderXtrain, loaderYtrain, hparams) - - bar(get_window_of_Aₖ(Normal(0.0f0, 1.0f0), nn_model, collect(loaderXtrain)[1], 2)) - - res = [] - for i in (1:(hparams.epochs - 1)) - Flux.reset!(nn_model) - for data in collect(loaderXtrain)[i] - nn_model.([[data]]) - end - - prediction = Vector{Float32}() - for data in collect(loaderXtest)[i] - y = nn_model.([[data]]) - append!(prediction, y[1]) - end - r, _ = yule_walker(Float64.(prediction); order=3) - push!(res, r) - end - - plot(prediction; seriestype=:scatter) - plot!(Float32.(collect(loaderXtest)[1]); seriestype=:scatter) - - ND(Float32.(collect(loaderXtest)[1])[1:200], prediction[1:200]) - - RMSE(Float32.(collect(loaderXtest)[1])[1:200], prediction[1:200]) - - yule_walker(Float64.(collect(loaderYtest)[2]); order=3) - yule_walker(Float64.(prediction); order=3) - - y = collect(loaderYtest)[1] - Flux.reset!(nn_model) - nn_model.([collect(loaderXtest)[1]']) - collect(loaderYtrain)[1] - - get_watson_durbin_test(y, ŷ) -end - -@test_experiments "testing" begin - csv_file_path = "/Users/jmfrutos/Desktop/LD2011_2014.txt" - - df = CSV.File( - csv_file_path; - delim=';', - header=true, - decimal=',', - types=Dict( - "MT_330" => Float32, - "MT_331" => Float32, - "MT_332" => Float32, - "MT_333" => Float32, - "MT_334" => Float32, - "MT_335" => Float32, - ), - ) - - hparams = HyperParamsTS(; seed=1234, η=1e-3, epochs=1, window_size=1000, K=5) - - nn_model = Chain(RNN(1 => 32, relu), RNN(32 => 32, relu), Dense(32 => 1, identity)) - - loaderXtrain = Flux.DataLoader( - df.MT_334[1:40000]; batchsize=round(Int, 40000), shuffle=false, partial=false - ) - - loaderYtrain = Flux.DataLoader( - df.MT_334[2:40001]; batchsize=round(Int, 40000), shuffle=false, partial=false - ) - - res = Vector{Float32}() - @showprogress for i in 1:200 - loss = ts_adaptative_block_learning(nn_model, loaderXtrain, loaderYtrain, hparams) - append!(res, loss) - end - - Flux.reset!(nn_model) - for data in collect(loaderXtrain)[1] - nn_model.([[data]]) - end - - prediction = Vector{Float32}() - for data in df.MT_334[39900:40100] - y = nn_model.([[data]]) - append!(prediction, y[1]) - end -end - -@test_experiments "testing" begin - csv_file_path = "/Users/jmfrutos/Desktop/LD2011_2014.txt" - - df = CSV.File( - csv_file_path; - delim=';', - header=true, - decimal=',', - types=Dict( - "MT_331" => Float32, - "MT_332" => Float32, - "MT_333" => Float32, - "MT_334" => Float32, - "MT_335" => Float32, - "MT_336" => Float32, - "MT_338" => Float32, - ), - ) - - hparams = HyperParamsTS(; seed=1234, η=1e-3, epochs=100, window_size=100, K=5) - - nn_model = Chain(RNN(5 => 32, relu), RNN(32 => 32, relu), Dense(32 => 1, identity)) - - loaderXtrain = Flux.DataLoader( - [ - [df.MT_333[i], df.MT_334[i], df.MT_335[i], df.MT_336[i], df.MT_338[i]] for - i in 1:length(df.MT_335[1:40000]) - ]; - batchsize=round(Int, 40000), - shuffle=false, - partial=false, - ) - - loaderYtrain = Flux.DataLoader( - [ - [df.MT_333[i], df.MT_334[i], df.MT_335[i], df.MT_336[i], df.MT_338[i]] for - i in 1:length(df.MT_335[1:40000]) - ]; - batchsize=round(Int, 200), - shuffle=false, - partial=false, - ) - - loaderXtest = Flux.DataLoader( - [ - [df.MT_333[40000+i], df.MT_334[40000+i], df.MT_335[40000+i], df.MT_336[40000+i], df.MT_338[40000+i]] for - i in 1:length(df.MT_335[40000:40200]) - ]; - batchsize=round(Int, 40000), - shuffle=false, - partial=false, - ) - - #loss = ts_covariates_mse_learning(nn_model, collect(loaderXtrain)[1], collect(loaderYtrain)[1], hparams) - - loss = ts_covariates_adaptative_block_learning(nn_model, collect(loaderXtrain)[1], collect(loaderYtrain)[1], hparams) - - res = [] - Flux.reset!(nn_model) - for data in collect(loaderXtrain)[1] - y = nn_model.([data]) - append!(prediction, y[1]) - end - - prediction = Vector{Float32}() - for data in collect(loaderXtest)[1] - y = nn_model.([data]) - append!(prediction, y[1]) - end - - prediction = Vector{Float32}() - data = collect(loaderXtest)[1] - y = nn_model(data[1]) - append!(prediction, y[1]) - for i in 1:200 - y = nn_model(vcat(y, data[i][2:end])) - append!(prediction, y[1]) - end - r, _ = yule_walker(Float64.(prediction); order=3) - push!(res, r) -end diff --git a/src/CustomLossFunction.jl b/src/CustomLossFunction.jl index d5cf42c..12a175c 100644 --- a/src/CustomLossFunction.jl +++ b/src/CustomLossFunction.jl @@ -8,10 +8,9 @@ function _sigmoid(ŷ::Matrix{T}, y::T) where {T<:AbstractFloat} end; function _leaky_relu(ŷ::Matrix{T}, y::T) where {T<:AbstractFloat} - return min.(0.001 .* (y .- ŷ) .+ 1., leakyrelu.((y .- ŷ) .* 10, 0.001)) + return min.(0.001 .* (y .- ŷ) .+ 1.0, leakyrelu.((y .- ŷ) .* 10, 0.001)) end; - """ ψₘ(y, m) @@ -140,7 +139,7 @@ function adaptative_block_learning(nn_model, data, hparams) @showprogress for epoch in 1:(hparams.epochs) loss, grads = Flux.withgradient(nn_model) do nn aₖ = zeros(hparams.K + 1) - for i in 1:hparams.samples + for i in 1:(hparams.samples) x = rand(hparams.transform, hparams.K) yₖ = nn(x') aₖ += generate_aₖ(yₖ, data.data[i]) @@ -161,7 +160,7 @@ function adaptative_block_learning_1(nn_model, loader, hparams) @showprogress for data in loader loss, grads = Flux.withgradient(nn_model) do nn aₖ = zeros(hparams.K + 1) - for i in 1:hparams.samples + for i in 1:(hparams.samples) x = rand(hparams.transform, hparams.K) yₖ = nn(x') aₖ += generate_aₖ(yₖ, data[i]) @@ -203,9 +202,7 @@ end; Generate a window of the rv's `Aₖ` for a given model and target function. """ function get_window_of_Aₖ(transform, model, data, K::Int64) - window = count.([ - model(rand(transform, K)') .< d for d in data - ]) + window = count.([model(rand(transform, K)') .< d for d in data]) return [count(x -> x == i, window) for i in 0:K] end; @@ -221,7 +218,7 @@ end; function get_better_K(nn_model, data, min_K, hparams) K = hparams.max_k - for k in min_K:hparams.max_k + for k in min_K:(hparams.max_k) if !convergence_to_uniform(get_window_of_Aₖ(hparams.transform, nn_model, data, k)) K = k break @@ -261,7 +258,7 @@ function auto_adaptative_block_learning(nn_model, data, hparams) end loss, grads = Flux.withgradient(nn_model) do nn aₖ = zeros(K + 1) - for i in 1:hparams.samples + for i in 1:(hparams.samples) x = rand(hparams.transform, K) yₖ = nn(x') aₖ += generate_aₖ(yₖ, data.data[i]) @@ -290,7 +287,7 @@ function auto_adaptative_block_learning_1(nn_model, loader, hparams) end loss, grads = Flux.withgradient(nn_model) do nn aₖ = zeros(K + 1) - for i in 1:hparams.samples + for i in 1:(hparams.samples) x = rand(hparams.transform, K) yₖ = nn(x') aₖ += generate_aₖ(yₖ, data[i]) @@ -302,3 +299,144 @@ function auto_adaptative_block_learning_1(nn_model, loader, hparams) end return losses end; + +# Hyperparameters for the method `ts_adaptative_block_learning` + +""" + HyperParamsTS + +Hyperparameters for the method `ts_adaptative_block_learning` +""" +Base.@kwdef mutable struct HyperParamsTS + seed::Int = 72 # Random seed + dev = cpu # Device: cpu or gpu + η::Float64 = 1e-3 # Learning rate + epochs::Int = 100 # Number of epochs + noise_model = Normal(0.0f0, 1.0f0) # Noise to add to the data + window_size = 100 # Window size for the histogram + K = 10 # Number of simulted observations +end + +# Train and output the model according to the chosen hyperparameters `hparams` + +""" + ts_adaptative_block_learning(nn_model, Xₜ, Xₜ₊₁, hparams) + +Custom loss function for the model `nn_model` on time series data `Xₜ` and `Xₜ₊₁`. +""" +function ts_adaptative_block_learning(nn_model, Xₜ, Xₜ₊₁, hparams) + losses = [] + optim = Flux.setup(Flux.Adam(hparams.η), nn_model) + @showprogress for (batch_Xₜ, batch_Xₜ₊₁) in zip(Xₜ, Xₜ₊₁) + j = 0 + Flux.reset!(nn_model) + nn_model([batch_Xₜ[1]]) # Warm up recurrent model on first observation + loss, grads = Flux.withgradient(nn_model) do nn + aₖ = zeros(hparams.K + 1) + for i in (1:(hparams.window_size)) + xₖ = rand(hparams.noise_model, hparams.K) + nn_cp = deepcopy(nn) + yₖ = nn_cp(xₖ') + aₖ += generate_aₖ(yₖ, batch_Xₜ₊₁[j + i]) + nn([batch_Xₜ[j + i]]) + end + scalar_diff(aₖ ./ sum(aₖ)) + end + Flux.update!(optim, nn_model, grads[1]) + j += hparams.window_size + push!(losses, loss) + end + return losses +end + +""" + ts_covariates_adaptative_block_learning(nn_model, Xₜ, Xₜ₊₁, hparams) + +Custom loss function for the model `nn_model` on time series data `Xₜ` and `Xₜ₊₁`. +Where `Xₜ` is a vector of vectors where each vector [x₁, x₂, ..., xₙ] is a time series where +we want to predict the value at position x₁ using the values [x₂, ..., xₙ] as covariate and +`Xₜ₊₁` is a vector. +""" +function ts_covariates_adaptative_block_learning(nn_model, Xₜ, Xₜ₊₁, hparams) + losses = [] + optim = Flux.setup(Flux.Adam(hparams.η), nn_model) + @showprogress for _ in 1:(hparams.epochs) + j = 0 + Flux.reset!(nn_model) + nn_model(Xₜ[1]) + loss, grads = Flux.withgradient(nn_model) do nn + aₖ = zeros(hparams.K + 1) + for i in (1:(hparams.window_size)) + xₖ = rand(hparams.noise_model, hparams.K) + nn_cp = deepcopy(nn) + yₖ = nn_cp.([vcat(x, Xₜ[j + i][2:end]) for x in xₖ]) + aₖ += generate_aₖ(cat(yₖ...; dims=2), Xₜ₊₁[j + i][1]) + nn(Xₜ[j + i]) + end + scalar_diff(aₖ ./ sum(aₖ)) + end + Flux.update!(optim, nn_model, grads[1]) + j += hparams.window_size + push!(losses, loss) + end + return losses +end + +""" + get_proxy_histogram_loss_ts(nn_model, data_xₜ, data_xₜ₊₁, hparams) + +Get the loss of the model `nn_model` on the data `data_xₜ` and `data_xₜ₊₁` using the +proxy histogram method. +""" +function get_proxy_histogram_loss_ts(nn_model, data_xₜ, data_xₜ₊₁, hparams) + losses = [] + @showprogress for (batch_xₜ, batch_xₜ₊₁) in zip(data_xₜ, data_xₜ₊₁) + j = 0 + Flux.reset!(nn_model) + nn_model([batch_xₜ[1]]) + aₖ = zeros(hparams.K + 1) + for i in 1:(hparams.window_size) + xₖ = rand(hparams.noise_model, hparams.K) + nn_cp = deepcopy(nn_model) + yₖ = nn_cp(xₖ') + aₖ += generate_aₖ(yₖ, batch_xₜ₊₁[j + i]) + nn_model([batch_xₜ[j + i]]) + end + j += hparams.window_size + push!(losses, scalar_diff(aₖ ./ sum(aₖ))) + end + return losses +end + +function get_window_of_Aₖ_ts(transform, model, data, K::Int64) + Flux.reset!(model) + res = [] + for d in data + xₖ = rand(transform, K) + model_cp = deepcopy(model) + yₖ = model_cp(xₖ') + push!(res, yₖ .< d) + end + return [count(x -> x == i, count.(res)) for i in 0:K] +end; + +""" + get_density(nn, data, t, m) + +Generate `m` samples from the model `nn` given the data `data` up to time `t`. +Can be used to generate the histogram at time `t` of the model. +""" +function get_density(nn, data, t, m) + Flux.reset!(nn) + res = [] + for d in data[1:(t - 1)] + nn([d]) + end + for _ in 1:m + nn_cp = deepcopy(nn) + xₜ = rand(Normal(0.0f0, 1.0f0)) + yₜ = nn_cp([xₜ]) + append!(res, yₜ) + end + return res +end From 76c52e48ce63012764014461989878abb8eb687a Mon Sep 17 00:00:00 2001 From: josemanuel22 Date: Wed, 20 Sep 2023 12:29:58 +0200 Subject: [PATCH 11/12] refactor time_series_blocklearning --- src/AdaptativeBlockLearning.jl | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/AdaptativeBlockLearning.jl b/src/AdaptativeBlockLearning.jl index 0d8fac8..5f9f150 100644 --- a/src/AdaptativeBlockLearning.jl +++ b/src/AdaptativeBlockLearning.jl @@ -29,6 +29,12 @@ export _sigmoid, AutoAdaptativeHyperParams, auto_adaptative_block_learning, adaptative_block_learning_1, - auto_adaptative_block_learning_1 + auto_adaptative_block_learning_1, + HyperParamsTS, + ts_adaptative_block_learning, + ts_covariates_adaptative_block_learning, + get_proxy_histogram_loss_ts, + get_window_of_Aₖ_ts, + get_density end From 3b8863c5c9f672156fde2975520d07d4f8c1b30f Mon Sep 17 00:00:00 2001 From: josemanuel22 Date: Wed, 20 Sep 2023 12:30:28 +0200 Subject: [PATCH 12/12] refactor time_series_blocklearning --- examples/time_series_predictions/model.jl | 374 ++++++++++++++++++++++ 1 file changed, 374 insertions(+) create mode 100644 examples/time_series_predictions/model.jl diff --git a/examples/time_series_predictions/model.jl b/examples/time_series_predictions/model.jl new file mode 100644 index 0000000..4115903 --- /dev/null +++ b/examples/time_series_predictions/model.jl @@ -0,0 +1,374 @@ +using Flux +using Random +using Statistics + +using AdaptativeBlockLearning +using Distributions +using DataFrames +using CSV +using Plots + +include("../../benchmarks/benchmark_utils.jl") +include("ts_utils.jl") + +function ts_covariates_mse_learning(nn_model, Xₜ, Xₜ₊₁, hparams) + losses = [] + optim = Flux.setup(Flux.Adam(hparams.η), nn_model) + @showprogress for epoch in (1:(hparams.epochs)) + Flux.reset!(nn_model) + loss, grads = Flux.withgradient(nn_model) do nn + nn(Xₜ[1]) + sum(Flux.Losses.mse.([nn(x)[1] for x in Xₜ[1:end]], Xₜ₊₁[1:end])) + end + Flux.update!(optim, nn_model, grads[1]) + push!(losses, loss) + end + return losses +end + +function ts_mse_learning(nn_model, Xₜ, Xₜ₊₁, hparams) + losses = [] + optim = Flux.setup(Flux.Adam(hparams.η), nn_model) + @showprogress for epoch in (1:(hparams.epochs)) + Flux.reset!(nn_model) + loss, grads = Flux.withgradient(nn_model) do nn + nn([Xₜ[1]]') + sum(Flux.Losses.mse.([nn([x]')[1] for x in Xₜ[1:(end - 1)]], Xₜ₊₁[1:end])) + end + Flux.update!(optim, nn_model, grads[1]) + push!(losses, loss) + end + return losses +end + +@test_experiments "testing AutoRegressive Model" begin + ar_hparams = ARParams(; + ϕ=[0.7f0, 0.2f0, 0.1f0], + x₁=rand(Normal(0.0f0, 0.5f0)), + proclen=10000, + noise=Normal(0.0f0, 0.5f0), + ) + hparams = HyperParamsTS(; seed=1234, η=1e-3, epochs=1000, window_size=200, K=5) + + nn_model = Chain(RNN(1 => 32, relu), RNN(32 => 32, relu), Dense(32 => 1, identity)) + + loaderXtrain, loaderYtrain, loaderXtest, loaderYtest = generate_batch_train_test_data( + hparams, ar_hparams + ) + + loss = ts_adaptative_block_learning(nn_model, loaderXtrain, loaderYtrain, hparams) + + plot_ts(nn_model, loaderXtrain, loaderYtrain, hparams) + + loss = ts_mse_learning( + nn_model, collect(loaderXtrain)[1], collect(loaderYtrain)[1], hparams + ) + + @gif for i in 1:100 + histogram( + get_density(nn_model, collect(loaderXtrain)[1], i, 1000); + bins=(-25:0.2:20), + normalize=:pdf, + label="t=$i", + ) + println("$i") + end every 2 + + loss = get_stats(nn_model, loaderXtrain, loaderYtrain, hparams) + + bar(get_window_of_Aₖ(Normal(0.0f0, 1.0f0), nn_model, collect(loaderXtrain)[1], 2)) + + res = [] + for i in (1:(hparams.epochs - 1)) + Flux.reset!(nn_model) + for data in collect(loaderXtrain)[i] + nn_model.([[data]]) + end + + prediction = Vector{Float32}() + for data in collect(loaderXtest)[i] + y = nn_model.([[data]]) + append!(prediction, y[1]) + end + r, _ = yule_walker(Float64.(prediction); order=3) + push!(res, r) + end + + plot(prediction; seriestype=:scatter) + plot!(Float32.(collect(loaderXtest)[1]); seriestype=:scatter) + + ND(Float32.(collect(loaderXtest)[1])[1:200], prediction[1:200]) + + RMSE(Float32.(collect(loaderXtest)[1])[1:200], prediction[1:200]) + + yule_walker(Float64.(collect(loaderYtest)[2]); order=3) + yule_walker(Float64.(prediction); order=3) + + y = collect(loaderYtest)[1] + Flux.reset!(nn_model) + nn_model.([collect(loaderXtest)[1]']) + collect(loaderYtrain)[1] + + get_watson_durbin_test(y, ŷ) +end + +@test_experiments "testing al-pv-2006" begin + #= + https://www.nrel.gov/grid/solar-power-data.html + =# + + csv1 = "examples/time_series_predictions/data/al-pv-2006/Actual_30.45_-88.25_2006_UPV_70MW_5_Min.csv" + csv2 = "examples/time_series_predictions/data/al-pv-2006/Actual_30.55_-87.55_2006_UPV_80MW_5_Min.csv" + csv3 = "examples/time_series_predictions/data/al-pv-2006/Actual_30.55_-87.75_2006_DPV_36MW_5_Min.csv" + csv4 = "examples/time_series_predictions/data/al-pv-2006/Actual_30.55_-88.15_2006_DPV_38MW_5_Min.csv" + csv5 = "examples/time_series_predictions/data/al-pv-2006/Actual_30.55_-88.25_2006_DPV_38MW_5_Min.csv" + + df1 = CSV.File( + csv1; delim=',', header=true, decimal='.', types=Dict("Power(MW)" => Float32) + ) + + df2 = CSV.File( + csv2; delim=',', header=true, decimal='.', types=Dict("Power(MW)" => Float32) + ) + + df3 = CSV.File( + csv3; delim=',', header=true, decimal='.', types=Dict("Power(MW)" => Float32) + ) + + df4 = CSV.File( + csv4; delim=',', header=true, decimal='.', types=Dict("Power(MW)" => Float32) + ) + + df5 = CSV.File( + csv5; delim=',', header=true, decimal='.', types=Dict("Power(MW)" => Float32) + ) + + hparams = HyperParamsTS(; seed=1234, η=1e-3, epochs=100, window_size=1000, K=5) + + nn_model = Chain(RNN(5 => 32, relu), RNN(32 => 32, relu), Dense(32 => 1, identity)) + + num_training_data = 10000 + + loaderXtrain = Flux.DataLoader( + [ + [ + df1["Power(MW)"][i], + df1["Power(MW)"][i], + df1["Power(MW)"][i], + df1["Power(MW)"][i], + df1["Power(MW)"][i], + ] for i in 1:length(df1["Power(MW)"][1:num_training_data]) + ]; + batchsize=round(Int, num_training_data), + shuffle=false, + partial=false, + ) + + loaderYtrain = Flux.DataLoader( + [ + [ + df1["Power(MW)"][i], + df1["Power(MW)"][i], + df1["Power(MW)"][i], + df1["Power(MW)"][i], + df1["Power(MW)"][i], + ] for i in 1:length(df1["Power(MW)"][2:(num_training_data + 1)]) + ]; + batchsize=round(Int, num_training_data), + shuffle=false, + partial=false, + ) + + num_test_data = 200 + loaderXtest = Flux.DataLoader( + [ + [ + df1["Power(MW)"][i], + df1["Power(MW)"][i], + df1["Power(MW)"][i], + df1["Power(MW)"][i], + df1["Power(MW)"][i], + ] for i in + 1:length( + df1["Power(MW)"][num_training_data:(num_training_data + num_test_data)] + ) + ]; + batchsize=round(Int, num_training_data), + shuffle=false, + partial=false, + ) + + loss = ts_covariates_adaptative_block_learning( + nn_model, collect(loaderXtrain)[1], collect(loaderYtrain)[1], hparams + ) + + Flux.reset!(nn_model) + for data in collect(loaderXtrain)[1] + nn_model.([data]) + end + + prediction = Vector{Float32}() + for data in collect(loaderXtest)[1] + y = nn_model.([data]) + append!(prediction, y[1]) + end +end + +@test_experiments "testing 30 years european wind generation" begin + #= + https://www.nrel.gov/grid/solar-power-data.html + =# + + csv1 = "examples/time_series_predictions/data/30-years-of-european-wind-generation/TS.CF.N2.30yr.csv" + + df1 = CSV.File( + csv1; + delim=',', + header=true, + decimal='.', + stripwhitespace=true, + types=Dict("AT11" => Float32), + ) + + hparams = HyperParamsTS(; seed=1234, η=1e-3, epochs=100, window_size=1000, K=5) + + nn_model = Chain(RNN(5 => 32, relu), RNN(32 => 32, relu), Dense(32 => 1, identity)) + + num_training_data = 10000 + + loaderXtrain = Flux.DataLoader( + map( + x -> Float32.(x), + cat( + [ + [df1.FR43[i], df1.FR82[i], df1.FR71[i], df1.FR72[i], df1.FR26[i]] for + i in 1:length(df1[1:num_training_data]) + ]; + dims=1, + ), + ); + batchsize=round(Int, num_training_data), + shuffle=false, + partial=false, + ) + + loaderYtrain = Flux.DataLoader( + map( + x -> Float32.(x), + cat( + [ + [df1.FR43[i], df1.FR82[i], df1.FR71[i], df1.FR72[i], df1.FR26[i]] for + i in 1:length(df1[2:(num_training_data + 1)]) + ]; + dims=1, + ), + ); + batchsize=round(Int, num_training_data), + shuffle=false, + partial=false, + ) + + loss = ts_covariates_adaptative_block_learning( + nn_model, collect(loaderXtrain)[1], collect(loaderYtrain)[1], hparams + ) + + Flux.reset!(nn_model) + for data in collect(loaderXtrain)[1] + nn_model.([data]) + end + + prediction = Vector{Float32}() + for data in collect(loaderXtest)[1] + y = nn_model.([data]) + append!(prediction, y[1]) + end + + real_seriel = [x[1] for x in collect(loaderXtest)[1]] +end + +@test_experiments "testing LD2011_2014" begin + csv_file_path = "/AdaptativeBlockLearning/examples/time_series_predictions" + + df = CSV.File( + csv_file_path; + delim=';', + header=true, + decimal=',', + types=Dict( + "MT_331" => Float32, + "MT_332" => Float32, + "MT_333" => Float32, + "MT_334" => Float32, + "MT_335" => Float32, + "MT_336" => Float32, + "MT_338" => Float32, + ), + ) + + hparams = HyperParamsTS(; seed=1234, η=1e-3, epochs=100, window_size=100, K=5) + + nn_model = Chain(RNN(5 => 32, relu), RNN(32 => 32, relu), Dense(32 => 1, identity)) + + loaderXtrain = Flux.DataLoader( + [ + [df.MT_333[i], df.MT_334[i], df.MT_335[i], df.MT_336[i], df.MT_338[i]] for + i in 1:length(df.MT_335[1:40000]) + ]; + batchsize=round(Int, 40000), + shuffle=false, + partial=false, + ) + + loaderYtrain = Flux.DataLoader( + [ + [df.MT_333[i], df.MT_334[i], df.MT_335[i], df.MT_336[i], df.MT_338[i]] for + i in 1:length(df.MT_335[1:40000]) + ]; + batchsize=round(Int, 200), + shuffle=false, + partial=false, + ) + + loaderXtest = Flux.DataLoader( + [ + [ + df.MT_333[40000 + i], + df.MT_334[40000 + i], + df.MT_335[40000 + i], + df.MT_336[40000 + i], + df.MT_338[40000 + i], + ] for i in 1:length(df.MT_335[40000:40200]) + ]; + batchsize=round(Int, 40000), + shuffle=false, + partial=false, + ) + + loss = ts_covariates_adaptative_block_learning( + nn_model, collect(loaderXtrain)[1], collect(loaderYtrain)[1], hparams + ) + + res = [] + Flux.reset!(nn_model) + for data in collect(loaderXtrain)[1] + y = nn_model.([data]) + append!(prediction, y[1]) + end + + prediction = Vector{Float32}() + for data in collect(loaderXtest)[1] + y = nn_model.([data]) + append!(prediction, y[1]) + end + + prediction = Vector{Float32}() + data = collect(loaderXtest)[1] + y = nn_model(data[1]) + append!(prediction, y[1]) + for i in 1:200 + y = nn_model(vcat(y, data[i][2:end])) + append!(prediction, y[1]) + end + r, _ = yule_walker(Float64.(prediction); order=3) + push!(res, r) +end