validation_capr.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Retrieve research results (validation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To inspect the reported research results you must install the following julia packages:\n",
    "**JLD**, **PGFPlotsX**, **Plots**, **LaTeXStrings**, **CSV**, **DataFrames**\n",
    "\n",
    "Notational conventions: prefixes **m** and **v** denote mean and variance of the random variables respectively.\n",
    "\n",
    "For a sake of brevity, throughout the code we use $x$ and $z$ instead of $x^{(1)}$ and $x^{(2)}$ respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "┌ Warning: backend `PGFPlotsX` is not installed.\n",
      "└ @ Plots C:\\Users\\mkpap\\.julia\\packages\\Plots\\sxUvK\\src\\backends.jl:43\n"
     ]
    }
   ],
   "source": [
    "using JLD, Plots, PGFPlotsX, SparseArrays, CSV, DataFrames, LaTeXStrings\n",
    "pgfplotsx()\n",
    "push!(PGFPlotsX.CUSTOM_PREAMBLE, raw\"\\usepgfplotslibrary{fillbetween}\")\n",
    "\n",
    "# df = CSV.File(\"data/AAPL.csv\") |> DataFrame\n",
    "df = CSV.File(\"data/capr_prl_one_subject.csv\") |> DataFrame\n",
    "capr_one_subject = df[!, :reward_tally];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "shgf3l_results = JLD.load(\"dump/results_3shgf_stocks_mixture.jld\")\n",
    "shgf2l_results = JLD.load(\"dump/results_2shgf_stocks_mixture.jld\")\n",
    "shgf2l3c_results = JLD.load(\"dump/results_2l3cshgf_stocks_mixture.jld\")\n",
    "hgf2l_results  = JLD.load(\"dump/results_hgf_stocks_mixture.jld\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fe_shgf3l = shgf3l_results[\"results\"][\"fe\"][2:end]\n",
    "fe_shgf2l = shgf2l_results[\"results\"][\"fe\"][2:end]\n",
    "fe_shgf2l3c = shgf2l3c_results[\"results\"][\"fe\"][2:end]\n",
    "fe_hgf2l  = hgf2l_results[\"results\"][\"fe\"][2:end]\n",
    "\n",
    "fe_max_iter = 40\n",
    "free_energy_axis = @pgf Axis({xlabel=\"iteration\",\n",
    "           ylabel=\"free-enery [nats]\",\n",
    "           legend_pos = \"north east\",\n",
    "           grid = \"major\",\n",
    "           mark_options = {scale=0.3}, xmin=-5.0, xmax=45,\n",
    "    },\n",
    "    Plot({mark = \"square*\", \"blue\"},Coordinates(collect(1:fe_max_iter),fe_shgf3l[1:fe_max_iter])),\n",
    "    LegendEntry(\"3L-SHGF-2\"),\n",
    "    Plot({mark = \"triangle*\", \"red\"},Coordinates(collect(1:fe_max_iter),fe_shgf2l[1:fe_max_iter])),\n",
    "    LegendEntry(\"2L-SHGF-2\"),\n",
    "    Plot({mark = \"square\", \"green\"},Coordinates(collect(1:fe_max_iter),fe_shgf2l3c[1:fe_max_iter])),\n",
    "    LegendEntry(\"2L-SHGF-3\"),\n",
    "    Plot({mark = \"o\", \"black\"},Coordinates(collect(1:fe_max_iter),fe_hgf2l[1:fe_max_iter])),\n",
    "    LegendEntry(\"2L-HGF\"),\n",
    ")\n",
    "pgfsave(\"figures/validation_results_free_energy.tikz\", free_energy_axis)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mz_shgf3l = shgf3l_results[\"results\"][\"mz1\"]\n",
    "mz_shgf2l = shgf2l_results[\"results\"][\"mz\"]\n",
    "mz_shgf2l3c = shgf2l3c_results[\"results\"][\"mz\"]\n",
    "mz_hgf2l  = hgf2l_results[\"results\"][\"mz\"]\n",
    "vz_shgf3l = shgf3l_results[\"results\"][\"vz1\"]\n",
    "vz_shgf2l = shgf2l_results[\"results\"][\"vz\"]\n",
    "vz_shgf2l3c = shgf2l3c_results[\"results\"][\"vz\"]\n",
    "vz_hgf2l  = hgf2l_results[\"results\"][\"vz\"]\n",
    "n_samples = length(AAPL_stocks)\n",
    "axis3 = @pgf Axis({xlabel=L\"t\",\n",
    "           ylabel=L\"x_t^{(2)}\",\n",
    "        legend_pos = \"north west\",\n",
    "        mark_options = {scale=0.3},\n",
    "        grid=\"major\",\n",
    "        style = {thick}\n",
    "    },\n",
    "    \n",
    "    Plot({no_marks,color=\"blue\"},Coordinates(collect(1:n_samples), mz_shgf3l)),\n",
    "    LegendEntry(\"3L-SHGF-2\"),\n",
    "    Plot({no_marks,color=\"red\"}, Coordinates(collect(1:n_samples), mz_shgf2l)),\n",
    "    LegendEntry(\"2L-SHGF-2\"),\n",
    "    Plot({no_marks,color=\"green\"},Coordinates(collect(1:n_samples), mz_shgf2l3c)),\n",
    "    LegendEntry(\"2L-SHGF-3\"),\n",
    "    Plot({no_marks,color=\"black\"},Coordinates(collect(1:n_samples), mz_hgf2l)),\n",
    "    LegendEntry(\"2L-HGF\"),\n",
    "    \n",
    "    Plot({ \"name path=f1\", no_marks,color=\"blue\",opacity=0.2 }, Coordinates(collect(1:n_samples), mz_shgf3l .+  sqrt.(vz_shgf3l))),\n",
    "    Plot({ \"name path=g1\", no_marks, color=\"blue\",opacity=0.2}, Coordinates(collect(1:n_samples), mz_shgf3l .-  sqrt.(vz_shgf3l))),\n",
    "    Plot({ thick, color = \"blue\", fill = \"blue\", opacity = 0.2 },\n",
    "              raw\"fill between [of=f1 and g1]\"),\n",
    "    \n",
    "    Plot({ \"name path=f2\", no_marks,color=\"red\",opacity=0.2 }, Coordinates(collect(1:n_samples), mz_shgf2l .+  sqrt.(vz_shgf2l))),\n",
    "    Plot({ \"name path=g2\", no_marks, color=\"red\",opacity=0.2}, Coordinates(collect(1:n_samples), mz_shgf2l .-  sqrt.(vz_shgf2l))),\n",
    "    Plot({ thick, color = \"blue\", fill = \"red\", opacity = 0.2 },\n",
    "               raw\"fill between [of=f2 and g2]\"),\n",
    "    \n",
    "    \n",
    "    Plot({ \"name path=f3\", no_marks,color=\"green\",opacity=0.2 }, Coordinates(collect(1:n_samples), mz_shgf2l3c .+  sqrt.(vz_shgf2l3c))),\n",
    "    Plot({ \"name path=g3\", no_marks, color=\"green\",opacity=0.2}, Coordinates(collect(1:n_samples), mz_shgf2l3c .-  sqrt.(vz_shgf2l3c))),\n",
    "    Plot({ thick, color = \"green\", fill = \"green\", opacity = 0.2 },\n",
    "               raw\"fill between [of=f3 and g3]\"),\n",
    "    \n",
    "   \n",
    "    Plot({ \"name path=f4\", no_marks,color=\"black\",opacity=0.2 }, Coordinates(collect(1:n_samples), mz_hgf2l .+  sqrt.(vz_hgf2l))),\n",
    "    Plot({ \"name path=g4\", no_marks, color=\"black\",opacity=0.2}, Coordinates(collect(1:n_samples), mz_hgf2l .-  sqrt.(vz_hgf2l))),\n",
    "    Plot({ thick, color = \"blue\", fill = \"black\", opacity = 0.2 },\n",
    "            raw\"fill between [of=f4 and g4]\"),\n",
    ")\n",
    "\n",
    "pgfsave(\"figures/validation_results_vol.tikz\", axis3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mx_shgf3l = [\"mx1\"]\n",
    "mx_shgf2l = shgf2l_results[\"results\"][\"mx\"]\n",
    "mx_hgf2l  = hgf2l_results[\"results\"][\"mx\"]\n",
    "vx_shgf3l = shgf3l_results[\"results\"][\"vx\"]\n",
    "vx_shgf2l = shgf2l_results[\"results\"][\"vx\"]\n",
    "vx_hgf2l  = hgf2l_results[\"results\"][\"vx\"];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ms = shgf2l_results[\"results\"][\"ms\"]\n",
    "categories = [x[2] for x in findmax.(ms)]\n",
    "maxup = maximum(AAPL_stocks) + 3.0\n",
    "mindown = minimum(AAPL_stocks) - 3.0\n",
    "axis4 = @pgf Axis({xlabel=L\"t\",\n",
    "           ylabel=L\"x_t^{(1)}\",\n",
    "        legend_pos = \"north west\",\n",
    "        mark_options = {scale=0.3},\n",
    "        grid=\"major\",\n",
    "        style = {thick}\n",
    "    },\n",
    "    Plot({no_marks,color=\"magenta\"},Coordinates(collect(1:n_samples), mx_shgf2l)),\n",
    "    LegendEntry(\"2L-SHGF\"),\n",
    "    Plot({ \"name path=f\", no_marks,color=\"magenta\",opacity=0.2 }, Coordinates(collect(1:n_samples), mx_shgf2l .+  3*sqrt.(vx_shgf2l))),\n",
    "    Plot({ \"name path=g\", no_marks, color=\"magenta\",opacity=0.2}, Coordinates(collect(1:n_samples), mx_shgf2l .-  3*sqrt.(vx_shgf2l))),\n",
    "    Plot({ thick, color = \"magenta\", fill = \"magenta\", opacity = 0.2 },\n",
    "               raw\"fill between [of=f and g]\"),\n",
    "    Plot({only_marks,color=\"black\",opacity=0.3}, Coordinates(collect(1:n_samples),AAPL_stocks)),\n",
    "    Plot(\n",
    "       {only_marks, scatter, scatter_src = \"explicit\"},\n",
    "       Table(\n",
    "           {x = \"x\", y = \"y\", meta = \"col\"},\n",
    "            x = collect(1:n_samples), y = mindown*ones(n_samples), col = categories\n",
    "       ),\n",
    "    ),\n",
    ")\n",
    "\n",
    "#pgfsave(\"figures/validation_prices2.tikz\", axis4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mω = shgf3l_results[\"results\"][\"mω2\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Please restart kernel to use this snippet\n",
    "import Plots\n",
    "using JLD\n",
    "\n",
    "shgf3l_results = JLD.load(\"dump/results_3shgf_stocks_mixture.jld\")\n",
    "shgf2l_results = JLD.load(\"dump/results_2shgf_stocks_mixture.jld\")\n",
    "shgf2l3c_results = JLD.load(\"dump/results_2l3cshgf_stocks_mixture.jld\")\n",
    "hgf2l_results  = JLD.load(\"dump/results_hgf_stocks_mixture.jld\");\n",
    "\n",
    "fe_shgf3l = shgf3l_results[\"results\"][\"fe\"][2:end]\n",
    "fe_shgf2l = shgf2l_results[\"results\"][\"fe\"][2:end]\n",
    "fe_shgf2l3c = shgf2l3c_results[\"results\"][\"fe\"][2:end]\n",
    "fe_hgf2l  = hgf2l_results[\"results\"][\"fe\"][2:end]\n",
    "\n",
    "Plots.plot(fe_shgf3l, label=\"3L-SHGF\")\n",
    "Plots.plot!(fe_shgf2l, label=\"2L-SHGF\")\n",
    "Plots.plot!(fe_shgf2l3c, label=\"2L-SHGF, 3 cats\")\n",
    "Plots.plot!(fe_hgf2l, label=\"2L-HGF\", background_color=:grey95, ylabel=\"Free energy [nats]\", xlabel=\"iteration\")\n",
    "Plots.savefig(\"fe_comparison.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Implemetation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To run the inference, the following packages must be installed:\n",
    "**ForneyLab**, **ProgressMeter**, **CSV**, **DataFrames**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: replacing module SGCV.\n",
      "WARNING: using SGCV.prod! in module Main conflicts with an existing identifier.\n",
      "WARNING: replacing module GCV.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "prod! (generic function with 82 methods)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "using ForneyLab\n",
    "include(\"sgcv/SGCV.jl\")\n",
    "using Main.SGCV\n",
    "include(\"gcv/GCV.jl\")\n",
    "using Main.GCV\n",
    "using Plots\n",
    "using SparseArrays\n",
    "using ProgressMeter\n",
    "include(\"compatibility.jl\") # resolves compatibility issues with ForneyLab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: replacing module SGCV.\n",
      "WARNING: using SGCV.prod! in module Main conflicts with an existing identifier.\n",
      "WARNING: replacing module GCV.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "prod! (generic function with 82 methods)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "include(\"sgcv/SGCV.jl\")\n",
    "using Main.SGCV\n",
    "include(\"gcv/GCV.jl\")\n",
    "using Main.GCV\n",
    "using Plots\n",
    "using SparseArrays\n",
    "using ProgressMeter\n",
    "include(\"compatibility.jl\") # resolves compatibility issues with ForneyLab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "pad (generic function with 1 method)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Auxialry function for graph initializtion\n",
    "pad(sym::Symbol, t::Int) = sym*:_*Symbol(lpad(t,3,'0')) # Left-pads a number with zeros, converts it to symbol and appends to sym"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "using CSV, DataFrames\n",
    "# df = CSV.File(\"data/AAPL.csv\") |> DataFrame\n",
    "df = CSV.File(\"data/capr_prl_one_subject.csv\") |> DataFrame\n",
    "# series = df[!, :Open];\n",
    "series = df[!, :reward_tally];\n",
    "# scatter(series, markershape=:xcross, markersize=1.5, markeralpha=0.5, background_color=:grey97, label=\"APPL\")\n",
    "scatter(series, markershape=:xcross, markersize=1.5, markeralpha=0.5, background_color=:grey97, label=\"capr_prl_one_subject\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "ename": "PGFPlotsX.MissingExternalProgramError",
     "evalue": "No LaTeX installation found, figures will not be generated. Make sure either pdflatex, xelatex or lualatex are installed and that the PATH variable is correctly set.",
     "output_type": "error",
     "traceback": [
      "No LaTeX installation found, figures will not be generated. Make sure either pdflatex, xelatex or lualatex are installed and that the PATH variable is correctly set.\n",
      "\n",
      "Stacktrace:\n",
      "  [1] latexengine()\n",
      "    @ PGFPlotsX C:\\Users\\mkpap\\.julia\\packages\\PGFPlotsX\\wtVR1\\src\\build.jl:19\n",
      "  [2] show(io::IOStream, #unused#::MIME{Symbol(\"image/png\")}, p::TikzDocument)\n",
      "    @ PGFPlotsX C:\\Users\\mkpap\\.julia\\packages\\PGFPlotsX\\wtVR1\\src\\tikzdocument.jl:339\n",
      "  [3] _show(io::IOStream, mime::MIME{Symbol(\"image/png\")}, plt::Plots.Plot{Plots.PGFPlotsXBackend})\n",
      "    @ Plots C:\\Users\\mkpap\\.julia\\packages\\Plots\\sxUvK\\src\\backends\\pgfplotsx.jl:1371\n",
      "  [4] #invokelatest#2\n",
      "    @ .\\essentials.jl:819 [inlined]\n",
      "  [5] invokelatest\n",
      "    @ .\\essentials.jl:816 [inlined]\n",
      "  [6] show\n",
      "    @ C:\\Users\\mkpap\\.julia\\packages\\Plots\\sxUvK\\src\\output.jl:232 [inlined]\n",
      "  [7] #344\n",
      "    @ C:\\Users\\mkpap\\.julia\\packages\\Plots\\sxUvK\\src\\output.jl:6 [inlined]\n",
      "  [8] open(::Plots.var\"#344#345\"{Plots.Plot{Plots.PGFPlotsXBackend}}, ::String, ::Vararg{String}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})\n",
      "    @ Base .\\io.jl:395\n",
      "  [9] open\n",
      "    @ .\\io.jl:392 [inlined]\n",
      " [10] png\n",
      "    @ C:\\Users\\mkpap\\.julia\\packages\\Plots\\sxUvK\\src\\output.jl:6 [inlined]\n",
      " [11] frame(anim::Animation, plt::Plots.Plot{Plots.PGFPlotsXBackend})\n",
      "    @ Plots C:\\Users\\mkpap\\.julia\\packages\\Plots\\sxUvK\\src\\animation.jl:26\n",
      " [12] frame(anim::Animation)\n",
      "    @ Plots C:\\Users\\mkpap\\.julia\\packages\\Plots\\sxUvK\\src\\animation.jl:24\n",
      " [13] macro expansion\n",
      "    @ C:\\Users\\mkpap\\.julia\\packages\\Plots\\sxUvK\\src\\animation.jl:232 [inlined]\n",
      " [14] macro expansion\n",
      "    @ c:\\xampp\\htdocs\\SGCV\\validation_capr.ipynb:4 [inlined]\n",
      " [15] top-level scope\n",
      "    @ C:\\Users\\mkpap\\.julia\\packages\\Plots\\sxUvK\\src\\animation.jl:251"
     ]
    }
   ],
   "source": [
    "shift = 0.5\n",
    "@gif for i in 1:length(series)\n",
    "    scatter(series[1:i], markershape=:xcross, markersize=1.5, markeralpha=0.5, xlim = (0, length(series)), ylim = (minimum(series)-shift, maximum(series)+shift), label=\"capr_prl_one_subject\", xlabel=:trial, ylabel=\"reward_tally\", background_color=:grey95)\n",
    "end every 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n",
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"600\" height=\"400\" viewBox=\"0 0 2400 1600\">\n",
       "<defs>\n",
       "  <clipPath id=\"clip300\">\n",
       "    <rect x=\"47\" y=\"47\" width=\"2307\" height=\"1507\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<defs>\n",
       "  <clipPath id=\"clip301\">\n",
       "    <rect x=\"0\" y=\"0\" width=\"2400\" height=\"1600\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polyline clip-path=\"url(#clip301)\" style=\"stroke:#000000; stroke-linecap:round; stroke-linejoin:round; stroke-width:4; stroke-opacity:1; fill:none\" points=\"47.2441,1510.15 2352.76,1510.15 \"/>\n",
       "<polyline clip-path=\"url(#clip301)\" style=\"stroke:#000000; stroke-linecap:round; stroke-linejoin:round; stroke-width:4; stroke-opacity:1; fill:none\" points=\"1200,1552.76 1200,47.2441 \"/>\n",
       "<polyline clip-path=\"url(#clip300)\" style=\"stroke:#add8e6; stroke-linecap:round; stroke-linejoin:round; stroke-width:20; stroke-opacity:1; fill:none\" points=\"112.494,1508.74 119.575,1508.61 126.655,1508.46 192.122,1506.41 257.589,1502.26 296.932,1498.09 336.275,1492.04 369.539,1484.97 402.803,1475.59 420.871,1469.33 438.938,1462.12 457.005,1453.84 475.072,1444.4 491.728,1434.55 508.384,1423.51 525.04,1411.17 541.696,1397.45 559.554,1381.08 577.411,1362.89 595.268,1342.75 613.125,1320.57 633.815,1292.2 654.504,1260.82 675.193,1226.35 695.883,1188.74 714.269,1152.64 732.655,1114.07 751.041,1073.06 769.428,1029.72 785.551,989.875 801.674,948.441 817.797,905.548 833.92,861.355 869.077,761.401 904.234,658.448 939.6,554.635 974.966,453.734 984.741,426.809 994.517,400.439 1004.29,374.697 1014.07,349.655 1023.84,325.386 1033.62,301.96 1043.39,279.447 1053.17,257.914 1062.47,238.397 1071.77,219.882 1081.07,202.421 1090.37,186.064 1099.67,170.858 1108.98,156.848 1118.28,144.073 1127.58,132.572 1137.25,122.002 1146.92,112.879 1156.59,105.234 1166.26,99.0911 1175.93,94.4708 1185.6,91.3881 1195.27,89.8529 1204.94,89.8704 1213.09,91.0889 1221.23,93.4048 1229.37,96.8127 1237.51,101.305 1245.65,106.871 1253.79,113.497 1261.93,121.17 1270.07,129.87 1280.01,141.857 1289.95,155.306 1299.89,170.172 1309.83,186.405 1319.77,203.951 1329.71,222.752 1339.64,242.748 1349.58,263.874 1367.79,305.287 1386,349.832 1395.11,373.138 1404.21,397.055 1413.32,421.525 1422.42,446.489 1456.06,542.023 1489.69,640.559 1523.77,740.603 1557.84,838.321 1576.93,891.232 1596.02,942.402 1615.12,991.556 1634.21,1038.47 1651.51,1078.9 1668.81,1117.24 1686.11,1153.42 1703.41,1187.4 1723.22,1223.58 1743.03,1256.88 1762.84,1287.33 1782.65,1315.01 1800.99,1338.28 1819.33,1359.36 1837.67,1378.38 1856.01,1395.45 1874.09,1410.5 1892.18,1423.9 1910.26,1435.79 1928.34,1446.28 1945.9,1455.25 1963.45,1463.13 1981.01,1470.02 1998.56,1476.03 2034.56,1485.95 2070.55,1493.25 2109.99,1498.94 2149.43,1502.85 2206.71,1506.36 2264,1508.25 2275.75,1508.51 2287.51,1508.74 \"/>\n",
       "</svg>\n"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using StatsPlots\n",
    "import Distributions\n",
    "normal = Distributions.Normal(0.0, 1.0)\n",
    "plot(normal, background_color=:transparent, foreground_color=:black, grid=false, ticks = false, legend=false, linewidth=5.0, color=:lightblue, framestyle = :origin)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2-L SHGF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "generate_mp_2l (generic function with 1 method)"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Function for generating 2-L SHGF and message-passing (MP) algorithm\n",
    "function generate_mp_2l(ndim, n_samples)\n",
    "    fg = FactorGraph()\n",
    "    z = Vector{Variable}(undef, n_samples)\n",
    "    x = Vector{Variable}(undef, n_samples)\n",
    "    y = Vector{Variable}(undef, n_samples)\n",
    "    s = Vector{Variable}(undef, n_samples)\n",
    "    @RV A ~ Dirichlet(ones(ndim, ndim))\n",
    "    @RV ω ~ GaussianMeanPrecision(placeholder(:mω, dims=(ndim, )), placeholder(:wω, dims=(ndim, ndim)))\n",
    "    @RV [id=pad(:z,1)] z[1] ~ GaussianMeanPrecision(placeholder(:mz_prior1), placeholder(:wz_prior1))\n",
    "    @RV [id=pad(:x,1)] x[1] ~ GaussianMeanPrecision(placeholder(:mx_prior1), placeholder(:wx_prior1))\n",
    "    @RV [id=pad(:y,1)] y[1] ~ GaussianMeanPrecision(x[1], placeholder(:wy_prior1))\n",
    "    @RV [id=pad(:s,1)] s[1] ~ ForneyLab.Categorical(ones(ndim) ./ ndim)\n",
    "    placeholder(y[1], :y, index = 1)\n",
    "    for t in 2:n_samples\n",
    "        @RV [id=pad(:s,t)] s[t] ~ Transition(s[t-1], A)\n",
    "        @RV [id=pad(:z,t)] z[t] ~ GaussianMeanPrecision(z[t - 1], placeholder(pad(:wz_transition, t)))\n",
    "        @RV [id=pad(:x,t)] x[t] ~ SwitchingGaussianControlledVariance(x[t - 1], z[t],ones(ndim), ω,s[t])\n",
    "        @RV [id=pad(:y,t)] y[t] ~ GaussianMeanPrecision(x[t], placeholder(pad(:wy_transition, t)))\n",
    "        placeholder(y[t], :y, index = t)\n",
    "    end\n",
    "    q = PosteriorFactorization(x, z, ω, s, A, ids=[:X :Z :Ω :S :A])\n",
    "    algo = messagePassingAlgorithm(free_energy=true)\n",
    "    src_code = algorithmSourceCode(algo, free_energy=true);\n",
    "    return src_code\n",
    "end\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "mp_2l (generic function with 1 method)"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "# Function for running MP-based inference for 2-L SHGF\n",
    "function mp_2l(obs;\n",
    "    ndims,\n",
    "    n_its = 100,\n",
    "    wy_prior1 = 1.0,\n",
    "    κ_m_prior = ones(ndims),\n",
    "    ω_m_prior = omegas,\n",
    "    κ_w_prior =  huge .* diageye(ndims),\n",
    "    ω_w_prior = 1.0 * diageye(ndims),\n",
    "    z_m_prior = 0.0,\n",
    "    z_w_prior = 1.0,\n",
    "    x_m_prior = 0.0,\n",
    "    x_w_prior = 1.0,\n",
    "    x_x_m_prior = zeros(ndims),\n",
    "    x_x_w_prior = 1.0*diageye(ndims),\n",
    "    z_z_m_prior = zeros(ndims),\n",
    "    z_z_w_prior = 1.0*diageye(ndims),\n",
    "    z_w_transition_prior = 100.0,\n",
    "    y_w_transition_prior =  1/mnv,\n",
    ")\n",
    "    n_samples = length(obs)\n",
    "    marginals = Dict()\n",
    "    marginals[:A] = ProbabilityDistribution(ForneyLab.MatrixVariate, Dirichlet, a=ones(ndims, ndims))\n",
    "    marginals[:κ] = ProbabilityDistribution(ForneyLab.Multivariate, GaussianMeanPrecision, m = κ_m_prior, w = κ_w_prior)\n",
    "    marginals[:ω] = ProbabilityDistribution(ForneyLab.Multivariate, GaussianMeanPrecision, m = ω_m_prior, w = ω_w_prior)\n",
    "    marginals[pad(:z,1)] = vague(GaussianMeanPrecision)\n",
    "    marginals[pad(:x,1)] = vague(GaussianMeanPrecision)\n",
    "    marginals[pad(:s,1)] = vague(Categorical, ndims)\n",
    "    for t = 2:n_samples\n",
    "        marginals[pad(:z,t)] = ProbabilityDistribution(ForneyLab.Univariate, GaussianMeanPrecision, m = z_m_prior, w = z_w_prior)\n",
    "        marginals[pad(:x,t)] = ProbabilityDistribution(ForneyLab.Univariate, GaussianMeanPrecision, m = x_m_prior, w = x_w_prior)\n",
    "        marginals[pad(:s,t)] = ProbabilityDistribution(Categorical, p = ones(ndims) ./ ndims)\n",
    "        marginals[pad(:s,t)*:_*pad(:s,t-1)] = ProbabilityDistribution(Contingency,p=ones(ndims, ndims) ./ ndims)\n",
    "        marginals[pad(:z,t)*:_*pad(:z,t-1)] = ProbabilityDistribution(ForneyLab.Multivariate,GaussianMeanPrecision, m = z_z_m_prior, w = z_z_w_prior)\n",
    "        marginals[pad(:x,t)*:_*pad(:x,t-1)] = ProbabilityDistribution(ForneyLab.Multivariate,GaussianMeanPrecision, m = x_x_m_prior, w = x_x_w_prior)\n",
    "    end\n",
    "    data = Dict()\n",
    "    data[:y] = obs\n",
    "    data[:mz_prior1] = z_m_prior\n",
    "    data[:wz_prior1] = z_w_prior\n",
    "    data[:mx_prior1] = x_m_prior\n",
    "    data[:wx_prior1] = x_w_prior\n",
    "    data[:wy_prior1] = wy_prior1\n",
    "    data[:mω] = ω_m_prior\n",
    "    data[:wω] = ω_w_prior\n",
    "    for t = 1:n_samples\n",
    "        data[pad(:wz_transition, t)] = z_w_transition_prior\n",
    "        data[pad(:wy_transition, t)] = y_w_transition_prior\n",
    "    end\n",
    "\n",
    "\n",
    "    fe = Vector{Float64}(undef, n_its)\n",
    "\n",
    "    @showprogress \"Iterations\" for i = 1:n_its\n",
    "\n",
    "        stepX!(data, marginals)\n",
    "        stepS!(data, marginals)\n",
    "        stepA!(data, marginals)\n",
    "        stepΩ!(data, marginals)\n",
    "        stepZ!(data, marginals)\n",
    "\n",
    "        fe[i] = freeEnergy(data, marginals)\n",
    "    end\n",
    "\n",
    "    mz = [ForneyLab.unsafeMean(marginals[pad(:z,t)]) for t=1:n_samples]\n",
    "    vz = [ForneyLab.unsafeVar(marginals[pad(:z,t)]) for t=1:n_samples]\n",
    "    mω = ForneyLab.unsafeMean(marginals[:ω])\n",
    "    vω = ForneyLab.unsafeCov(marginals[:ω])\n",
    "    mx = [ForneyLab.unsafeMean(marginals[pad(:x,t)]) for t=1:n_samples]\n",
    "    vx = [ForneyLab.unsafeVar(marginals[pad(:x,t)]) for t=1:n_samples]\n",
    "    ms = [ForneyLab.unsafeMean(marginals[pad(:s,t)]) for t=1:n_samples]\n",
    "    return mz,vz,mω, vω, mx,vx,ms,fe\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Vector{Float64}:\n",
       " -8.0\n",
       " 10.0"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Empirically selected prior\n",
    "omegas = [-8.0, 10.0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "No applicable StructuredVariationalRule{SwitchingGaussianControlledVariance} for interface :z of switchinggaussiancontrolledvariance_1 with inbound types:\n[Distribution, Nothing, Distribution, Distribution, Distribution]",
     "output_type": "error",
     "traceback": [
      "No applicable StructuredVariationalRule{SwitchingGaussianControlledVariance} for interface :z of switchinggaussiancontrolledvariance_1 with inbound types:\n[Distribution, Nothing, Distribution, Distribution, Distribution]",
      "",
      "Stacktrace:",
      " [1] error(s::String)",
      "   @ Base .\\error.jl:33",
      " [2] selectApplicableRule(rule_type::Type, entry::ScheduleEntry, inbound_types::Vector{Type}, applicable_rules::Vector{Type})",
      "   @ ForneyLab C:\\Users\\mkpap\\.julia\\packages\\ForneyLab\\EnMB2\\src\\algorithms\\inference_algorithm.jl:137",
      " [3] inferUpdateRule!(entry::ScheduleEntry, rule_type::Type{StructuredVariationalRule{SwitchingGaussianControlledVariance}}, inferred_outbound_types::Dict{Union{Nothing, Interface}, Type})",
      "   @ ForneyLab C:\\Users\\mkpap\\.julia\\packages\\ForneyLab\\EnMB2\\src\\algorithms\\structured_variational_bayes.jl:29",
      " [4] inferUpdateRules!(schedule::Vector{ScheduleEntry}; inferred_outbound_types::Dict{Union{Nothing, Interface}, Type})",
      "   @ ForneyLab C:\\Users\\mkpap\\.julia\\packages\\ForneyLab\\EnMB2\\src\\message_passing.jl:211",
      " [5] messagePassingSchedule(pf::PosteriorFactor)",
      "   @ ForneyLab C:\\Users\\mkpap\\.julia\\packages\\ForneyLab\\EnMB2\\src\\algorithms\\posterior_factor.jl:75",
      " [6] messagePassingAlgorithm(target_variables::Vector{Variable}, pfz::PosteriorFactorization; ep_sites::Vector{Tuple}, id::Symbol, free_energy::Bool)",
      "   @ ForneyLab C:\\Users\\mkpap\\.julia\\packages\\ForneyLab\\EnMB2\\src\\algorithms\\inference_algorithm.jl:82",
      " [7] #messagePassingAlgorithm#248",
      "   @ C:\\Users\\mkpap\\.julia\\packages\\ForneyLab\\EnMB2\\src\\algorithms\\inference_algorithm.jl:108 [inlined]",
      " [8] generate_mp_2l(ndim::Int64, n_samples::Int64)",
      "   @ Main .\\In[24]:23",
      " [9] top-level scope",
      "   @ In[26]:2"
     ]
    }
   ],
   "source": [
    "n_cats = length(omegas)\n",
    "code = generate_mp_2l(n_cats, length(series))\n",
    "eval(Meta.parse(code))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "ename": "UndefVarError",
     "evalue": "UndefVarError: `stepX!` not defined",
     "output_type": "error",
     "traceback": [
      "UndefVarError: `stepX!` not defined\n",
      "\n",
      "Stacktrace:\n",
      " [1] macro expansion\n",
      "   @ c:\\xampp\\htdocs\\SGCV\\validation_capr.ipynb:57 [inlined]\n",
      " [2] macro expansion\n",
      "   @ C:\\Users\\mkpap\\.julia\\packages\\ProgressMeter\\sN2xr\\src\\ProgressMeter.jl:938 [inlined]\n",
      " [3] mp_2l(obs::Vector{Float64}; ndims::Int64, n_its::Int64, wy_prior1::Float64, κ_m_prior::Vector{Float64}, ω_m_prior::Vector{Float64}, κ_w_prior::LinearAlgebra.Diagonal{Float64, Vector{Float64}}, ω_w_prior::LinearAlgebra.Diagonal{Float64, Vector{Float64}}, z_m_prior::Float64, z_w_prior::Float64, x_m_prior::Float64, x_w_prior::Float64, x_x_m_prior::Vector{Float64}, x_x_w_prior::LinearAlgebra.Diagonal{Float64, Vector{Float64}}, z_z_m_prior::Vector{Float64}, z_z_w_prior::LinearAlgebra.Diagonal{Float64, Vector{Float64}}, z_w_transition_prior::Float64, y_w_transition_prior::Float64)\n",
      "   @ Main c:\\xampp\\htdocs\\SGCV\\validation_capr.ipynb:55\n",
      " [4] top-level scope\n",
      "   @ c:\\xampp\\htdocs\\SGCV\\validation_capr.ipynb:1"
     ]
    }
   ],
   "source": [
    "mz,vz,mω, vω, mx,vx,ms,fe = mp_2l(series, x_m_prior=series[1], ndims=n_cats, z_w_prior=100.0, z_z_w_prior=100.0*diageye(2), z_w_transition_prior=100.0, ω_m_prior=omegas, ω_w_prior=diageye(2), y_w_transition_prior=1e6)\n",
    "\n",
    "results =   Dict(\"mz\" => mz, \"vz\" => vz,\n",
    "                  \"mx\" => mx, \"vx\" => vx,\n",
    "                  \"ms\" => ms, \"fe\" => fe,\n",
    "                  \"mω\" => mω, \"vω\" => vω,\n",
    "                  \"ωprior\" => omegas);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "ename": "UndefVarError",
     "evalue": "UndefVarError: `mω` not defined",
     "output_type": "error",
     "traceback": [
      "UndefVarError: `mω` not defined\n"
     ]
    }
   ],
   "source": [
    "mω"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "ename": "UndefVarError",
     "evalue": "UndefVarError: `scatter` not defined",
     "output_type": "error",
     "traceback": [
      "UndefVarError: `scatter` not defined\n",
      "\n",
      "Stacktrace:\n",
      " [1] top-level scope\n",
      "   @ c:\\xampp\\htdocs\\SGCV\\validation_capr.ipynb:4"
     ]
    }
   ],
   "source": [
    "xlims_ = 1, 350\n",
    "switch_line = 52\n",
    "shift = 0.05\n",
    "scatter(series, background_color=:grey95, xlabel=:Day, ylabel=\"Price \\$\", label=:APPL, markershape=:xcross, markersize=2, markeralpha=1.0, color=:black)\n",
    "plot!(mx, ribbon=sqrt.(vx), label=false)\n",
    "scatter!(switch_line*ones(length(ms)), color=argmax.(ms), markersize=5.0, markeralpha=0.5, label=false)\n",
    "xlims!(xlims_)\n",
    "ylims!(mx[xlims_[1]]-shift, mx[xlims_[2]]+shift)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "ename": "UndefVarError",
     "evalue": "UndefVarError: `series` not defined",
     "output_type": "error",
     "traceback": [
      "UndefVarError: `series` not defined\n",
      "\n",
      "Stacktrace:\n",
      " [1] top-level scope\n",
      "   @ c:\\xampp\\htdocs\\SGCV\\validation_capr.ipynb:2"
     ]
    }
   ],
   "source": [
    "shift = 1.0\n",
    "plot(mz, ribbon=sqrt.(vz), label=false, xlim = (0, length(series)), ylim = (minimum(mz)-shift, maximum(mz)+shift) ,background_color=:grey95, xlabel=:Day)\n",
    "xlims!(xlims_)\n",
    "ylims!(mz[xlims_[1]]-shift, mz[xlims_[2]]+shift)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "ename": "UndefVarError",
     "evalue": "UndefVarError: `ms` not defined",
     "output_type": "error",
     "traceback": [
      "UndefVarError: `ms` not defined\n",
      "\n",
      "Stacktrace:\n",
      " [1] top-level scope\n",
      "   @ c:\\xampp\\htdocs\\SGCV\\validation_capr.ipynb:1"
     ]
    }
   ],
   "source": [
    "scatter(argmax.(ms), color=argmax.(ms), markersize=3.0, markeralpha=0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "LoadError: UndefVarError: `@gif` not defined\nin expression starting at c:\\xampp\\htdocs\\SGCV\\validation_capr.ipynb:2",
     "output_type": "error",
     "traceback": [
      "LoadError: UndefVarError: `@gif` not defined\n",
      "in expression starting at c:\\xampp\\htdocs\\SGCV\\validation_capr.ipynb:2\n"
     ]
    }
   ],
   "source": [
    "shift = 10\n",
    "@gif for i in 1:length(series)\n",
    "    scatter(series[1:i], markershape=:xcross, markersize=1.5, markeralpha=0.5, xlim = (0, length(series)), ylim = (minimum(series)-shift, maximum(series)+shift), label=:APPL, xlabel=:Day, ylabel=\"Price \\$\", background_color=:grey95)\n",
    "    scatter!(50*ones(length(ms[1:i])), color=argmax.(ms[1:i]), markersize=5.0, markeralpha=0.5, label=false)    \n",
    "    plot!(mx[1:i], ribbon=sqrt.(vx[1:i]), label=false)\n",
    "end every 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "LoadError: UndefVarError: `@gif` not defined\nin expression starting at c:\\xampp\\htdocs\\SGCV\\validation_capr.ipynb:2",
     "output_type": "error",
     "traceback": [
      "LoadError: UndefVarError: `@gif` not defined\n",
      "in expression starting at c:\\xampp\\htdocs\\SGCV\\validation_capr.ipynb:2\n"
     ]
    }
   ],
   "source": [
    "shift = 1\n",
    "@gif for i in 1:length(series)\n",
    "    plot(mz[1:i], ribbon=sqrt.(vz[1:i]), label=false, xlim = (0, length(series)), ylim = (minimum(mz)-shift, maximum(mz)+shift) ,background_color=:grey95, xlabel=:Day)\n",
    "end every 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "ename": "UndefVarError",
     "evalue": "UndefVarError: `results` not defined",
     "output_type": "error",
     "traceback": [
      "UndefVarError: `results` not defined\n",
      "\n",
      "Stacktrace:\n",
      " [1] top-level scope\n",
      "   @ c:\\xampp\\htdocs\\SGCV\\validation_capr.ipynb:2"
     ]
    }
   ],
   "source": [
    "using JLD\n",
    "JLD.save(\"dump/results_2shgf_stocks_mixture.jld\",\"results\",results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3-L SHGF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "LoadError: UndefVarError: `@RV` not defined\nin expression starting at c:\\xampp\\htdocs\\SGCV\\validation_capr.ipynb:10",
     "output_type": "error",
     "traceback": [
      "LoadError: UndefVarError: `@RV` not defined\n",
      "in expression starting at c:\\xampp\\htdocs\\SGCV\\validation_capr.ipynb:10\n"
     ]
    }
   ],
   "source": [
    "# Function for generating 3-L SHGF and message-passing (MP) algorithm\n",
    "function generate_mp_3l(n_cats1, n_cats2, n_samples)\n",
    "    fg = FactorGraph()\n",
    "    z2 = Vector{Variable}(undef, n_samples)\n",
    "    s2 = Vector{Variable}(undef, n_samples)\n",
    "    z1 = Vector{Variable}(undef, n_samples)\n",
    "    s1 = Vector{Variable}(undef, n_samples)\n",
    "    x = Vector{Variable}(undef, n_samples)\n",
    "    y = Vector{Variable}(undef, n_samples)\n",
    "    @RV A2 ~ Dirichlet(ones(n_cats2, n_cats2))\n",
    "    @RV A1 ~ Dirichlet(ones(n_cats1, n_cats1))\n",
    "    @RV ω2 ~ GaussianMeanPrecision(placeholder(:mω2, dims=(n_cats2, )), placeholder(:wω2, dims=(n_cats2, n_cats2)))\n",
    "    @RV ω1 ~ GaussianMeanPrecision(placeholder(:mω1, dims=(n_cats1, )), placeholder(:wω1, dims=(n_cats1, n_cats1)))\n",
    "    @RV [id=pad(:z1,1)] z1[1] ~ GaussianMeanPrecision(placeholder(:mz1_prior1), placeholder(:wz1_prior1))\n",
    "    @RV [id=pad(:z2,1)] z2[1] ~ GaussianMeanPrecision(placeholder(:mz2_prior1), placeholder(:wz2_prior1))\n",
    "    @RV [id=pad(:x,1)] x[1] ~ GaussianMeanPrecision(placeholder(:mx_prior1), placeholder(:wx_prior1))\n",
    "    @RV [id=pad(:y,1)] y[1] ~ GaussianMeanPrecision(x[1], placeholder(:wy_prior1))\n",
    "    @RV [id=pad(:s1,1)] s1[1] ~ ForneyLab.Categorical(ones(n_cats1) ./ n_cats1)\n",
    "    @RV [id=pad(:s2,1)] s2[1] ~ ForneyLab.Categorical(ones(n_cats2) ./ n_cats2)\n",
    "    placeholder(y[1], :y, index = 1)\n",
    "    for t in 2:n_samples\n",
    "        @RV [id=pad(:s2, t)] s2[t] ~ Transition(s2[t-1], A2)\n",
    "        @RV [id=pad(:s1, t)] s1[t] ~ Transition(s1[t-1], A1)\n",
    "        @RV [id=pad(:z2,t)] z2[t] ~ GaussianMeanPrecision(z2[t - 1], placeholder(pad(:wz2_transition, t)))\n",
    "        @RV [id=pad(:z1,t)] z1[t] ~ SwitchingGaussianControlledVariance(z1[t - 1], z2[t], ones(n_cats2), ω2,s2[t])\n",
    "        @RV [id=pad(:x,t)] x[t] ~ SwitchingGaussianControlledVariance(x[t - 1], z1[t], ones(n_cats1), ω1,s1[t])\n",
    "        @RV [id=pad(:y,t)] y[t] ~ GaussianMeanPrecision(x[t], placeholder(pad(:wy_transition, t)))\n",
    "        placeholder(y[t], :y, index = t)\n",
    "    end\n",
    "    q = PosteriorFactorization(x, z1, s1, z2, s2, A1, A2, ω1, ω2, ids=[:X :Z1 :S1 :Z2 :S2 :A1 :A2 :Ω1 :Ω2])\n",
    "    algo = messagePassingAlgorithm(free_energy=true)\n",
    "    src_code = algorithmSourceCode(algo, free_energy=true);\n",
    "    return src_code\n",
    "end\n",
    "\n",
    "# Function for running MP-based inference for 3-L SHGF\n",
    "function mp_3l(obs;\n",
    "    n_cats1, n_cats2,\n",
    "    n_its = 20,\n",
    "    wy_prior1 = 1.0,\n",
    "    κ1_m_prior = ones(n_cats1),\n",
    "    ω1_m_prior = omegas1,\n",
    "    κ1_w_prior =  huge .* diageye(n_cats1),\n",
    "    ω1_w_prior = 1.0 * diageye(n_cats1),\n",
    "    z1_m_prior = 0.0,\n",
    "    z1_w_prior = 1.0,\n",
    "    κ2_m_prior = ones(n_cats2),\n",
    "    ω2_m_prior = omegas2,\n",
    "    κ2_w_prior =  huge .* diageye(n_cats2),\n",
    "    ω2_w_prior = 1.0 * diageye(n_cats2),\n",
    "    z2_m_prior = 0.0,\n",
    "    z2_w_prior = 1.0,\n",
    "    x_m_prior = 0.0,\n",
    "    x_w_prior = 1.0,\n",
    "    x_x_m_prior = zeros(n_cats1),\n",
    "    x_x_w_prior = 1.0*diageye(n_cats1),\n",
    "    z1_z1_m_prior = zeros(n_cats1),\n",
    "    z1_z1_w_prior = 1.0*diageye(n_cats1),\n",
    "    z1_w_transition_prior = 1000.0,\n",
    "    z2_z2_m_prior = zeros(n_cats2),\n",
    "    z2_z2_w_prior = 1.0*diageye(n_cats2),\n",
    "    z2_w_transition_prior = 10.0,\n",
    "    y_w_transition_prior =  1/mnv,\n",
    ")\n",
    "    n_samples = length(obs)\n",
    "    marginals = Dict()\n",
    "\n",
    "    # second  layer\n",
    "    marginals[:A1] = ProbabilityDistribution(ForneyLab.MatrixVariate, Dirichlet, a=ones(n_cats1, n_cats1))\n",
    "    marginals[:κ1] = ProbabilityDistribution(ForneyLab.Multivariate, GaussianMeanPrecision, m = κ1_m_prior, w = κ1_w_prior)\n",
    "    marginals[:ω1] = ProbabilityDistribution(ForneyLab.Multivariate, GaussianMeanPrecision, m = ω1_m_prior, w = ω1_w_prior)\n",
    "    marginals[pad(:z1,1)] = vague(GaussianMeanPrecision)\n",
    "\n",
    "    # third layer\n",
    "    marginals[:A2] = ProbabilityDistribution(ForneyLab.MatrixVariate, Dirichlet, a=ones(n_cats2, n_cats2))\n",
    "    marginals[:κ2] = ProbabilityDistribution(ForneyLab.Multivariate, GaussianMeanPrecision, m = κ2_m_prior, w = κ2_w_prior)\n",
    "    marginals[:ω2] = ProbabilityDistribution(ForneyLab.Multivariate, GaussianMeanPrecision, m = ω2_m_prior, w = ω2_w_prior)\n",
    "    marginals[pad(:z2,1)] = vague(GaussianMeanPrecision)\n",
    "\n",
    "    marginals[pad(:x,1)] = vague(GaussianMeanPrecision)\n",
    "    marginals[pad(:s1,1)] = vague(Categorical, n_cats1)\n",
    "    marginals[pad(:s2,1)] = vague(Categorical, n_cats2)\n",
    "    for t = 2:n_samples\n",
    "        marginals[pad(:z1,t)] = ProbabilityDistribution(ForneyLab.Univariate, GaussianMeanPrecision, m = z1_m_prior, w = z1_w_prior)\n",
    "        marginals[pad(:z2,t)] = ProbabilityDistribution(ForneyLab.Univariate, GaussianMeanPrecision, m = z2_m_prior, w = z2_w_prior)\n",
    "        marginals[pad(:x,t)] = ProbabilityDistribution(ForneyLab.Univariate, GaussianMeanPrecision, m = x_m_prior, w = x_w_prior)\n",
    "        marginals[pad(:s1,t)] = ProbabilityDistribution(Categorical, p = ones(n_cats1) ./ n_cats1)\n",
    "        marginals[pad(:s2,t)] = ProbabilityDistribution(Categorical, p = ones(n_cats2) ./ n_cats2)\n",
    "        marginals[pad(:s1,t)*:_*pad(:s1,t-1)] = ProbabilityDistribution(Contingency,p=ones(n_cats1, n_cats1) ./ n_cats1)\n",
    "        marginals[pad(:s2,t)*:_*pad(:s2,t-1)] = ProbabilityDistribution(Contingency,p=ones(n_cats2, n_cats2) ./ n_cats2)\n",
    "        marginals[pad(:z1,t)*:_*pad(:z1,t-1)] = ProbabilityDistribution(ForneyLab.Multivariate,GaussianMeanPrecision, m = z1_z1_m_prior, w = z2_z2_m_prior)\n",
    "        marginals[pad(:z2,t)*:_*pad(:z2,t-1)] = ProbabilityDistribution(ForneyLab.Multivariate,GaussianMeanPrecision, m = z2_z2_m_prior, w = z2_z2_w_prior)\n",
    "        marginals[pad(:x,t)*:_*pad(:x,t-1)] = ProbabilityDistribution(ForneyLab.Multivariate,GaussianMeanPrecision, m = x_x_m_prior, w = x_x_w_prior)\n",
    "    end\n",
    "    data = Dict()\n",
    "    data[:y] = obs\n",
    "    data[:mz1_prior1] = z1_m_prior\n",
    "    data[:wz1_prior1] = z1_w_prior\n",
    "    data[:mz2_prior1] = z2_m_prior\n",
    "    data[:wz2_prior1] = z2_w_prior\n",
    "    data[:mx_prior1] = x_m_prior\n",
    "    data[:wx_prior1] = x_w_prior\n",
    "    data[:wy_prior1] = wy_prior1\n",
    "    data[:mω1] = ω1_m_prior\n",
    "    data[:mω2] = ω2_m_prior\n",
    "    data[:wω1] = ω1_w_prior\n",
    "    data[:wω2] = ω2_w_prior\n",
    "    for t = 1:n_samples\n",
    "        data[pad(:wz1_transition, t)] = z1_w_transition_prior\n",
    "        data[pad(:wz2_transition, t)] = z2_w_transition_prior\n",
    "        data[pad(:wy_transition, t)] = y_w_transition_prior\n",
    "    end\n",
    "\n",
    "\n",
    "    fe = Vector{Float64}(undef, n_its)\n",
    "\n",
    "    @showprogress \"Iterations\" for i = 1:n_its\n",
    "\n",
    "        stepX!(data, marginals)\n",
    "        stepS1!(data, marginals)\n",
    "        stepA1!(data, marginals)\n",
    "        stepΩ1!(data, marginals)\n",
    "        stepZ1!(data, marginals)\n",
    "        stepS2!(data, marginals)\n",
    "        stepA2!(data, marginals)\n",
    "        stepΩ2!(data, marginals)\n",
    "        stepZ2!(data, marginals)\n",
    "\n",
    "        fe[i] = freeEnergy(data, marginals)\n",
    "    end\n",
    "\n",
    "    mz1 = [ForneyLab.unsafeMean(marginals[pad(:z1,t)]) for t=1:n_samples]\n",
    "    vz1 = [ForneyLab.unsafeVar(marginals[pad(:z1,t)]) for t=1:n_samples]\n",
    "    mω1 = ForneyLab.unsafeMean(marginals[:ω1])\n",
    "    vω1 = ForneyLab.unsafeCov(marginals[:ω1])\n",
    "    mz2 = [ForneyLab.unsafeMean(marginals[pad(:z2,t)]) for t=1:n_samples]\n",
    "    vz2 = [ForneyLab.unsafeVar(marginals[pad(:z2,t)]) for t=1:n_samples]\n",
    "    mω2 = ForneyLab.unsafeMean(marginals[:ω2])\n",
    "    vω2 = ForneyLab.unsafeCov(marginals[:ω2])\n",
    "    mx = [ForneyLab.unsafeMean(marginals[pad(:x,t)]) for t=1:n_samples]\n",
    "    vx = [ForneyLab.unsafeVar(marginals[pad(:x,t)]) for t=1:n_samples]\n",
    "    ms1 = [ForneyLab.unsafeMean(marginals[pad(:s1,t)]) for t=1:n_samples]\n",
    "    ms2 = [ForneyLab.unsafeMean(marginals[pad(:s2,t)]) for t=1:n_samples]\n",
    "    return mz1, vz1, mω1, vω1, mz2, vz2, mω2, vω2, mx, vx, ms1, ms2, fe\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Vector{Float64}:\n",
       " -3.0\n",
       "  1.0"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "omegas1 = [-1.0, 4.0]\n",
    "omegas2 = [-3.0, 1.0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "ename": "UndefVarError",
     "evalue": "UndefVarError: `series` not defined",
     "output_type": "error",
     "traceback": [
      "UndefVarError: `series` not defined\n",
      "\n",
      "Stacktrace:\n",
      " [1] top-level scope\n",
      "   @ c:\\xampp\\htdocs\\SGCV\\validation_capr.ipynb:1"
     ]
    }
   ],
   "source": [
    "code = generate_mp_3l(2, 2, length(series))\n",
    "eval(Meta.parse(code))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "ename": "UndefVarError",
     "evalue": "UndefVarError: `series` not defined",
     "output_type": "error",
     "traceback": [
      "UndefVarError: `series` not defined\n",
      "\n",
      "Stacktrace:\n",
      " [1] top-level scope\n",
      "   @ c:\\xampp\\htdocs\\SGCV\\validation_capr.ipynb:1"
     ]
    }
   ],
   "source": [
    "mz1, vz1, mω1, vω1, mz2, vz2, mω2, vω2, mx1, vx1, ms1, ms2, fe2 = mp_3l(series, n_its = 100, n_cats1=2, n_cats2=2,\n",
    "                                                                        x_m_prior=series[1],\n",
    "                                                                        ω1_m_prior=omegas1, ω2_m_prior=omegas2,\n",
    "                                                                        ω1_w_prior = 1.0 * diageye(2),\n",
    "                                                                        ω2_w_prior = 10.0 * diageye(2),\n",
    "                                                                        z1_w_transition_prior = 100.0,\n",
    "                                                                        z1_w_prior = 100.0,\n",
    "                                                                        z2_w_prior = 100.0,\n",
    "                                                                        y_w_transition_prior=1.0)\n",
    "\n",
    "results =   Dict(\"mz1\" => mz1, \"vz1\" => vz1,\n",
    "                 \"mz2\" => mz2, \"vz2\" => vz2,\n",
    "                 \"mx1\" => mx1, \"vx\" => vx1,\n",
    "                 \"ms1\" => ms1, \"ms2\" => ms2,\n",
    "                 \"fe\" => fe2,\n",
    "                 \"mω1\" => mω1, \"vω1\" => vω1,\n",
    "                 \"mω2\" => mω2, \"vω2\" => vω2,\n",
    "                 \"ωprior1\" => omegas1,\n",
    "                 \"ωprior2\" => omegas2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "ename": "UndefVarError",
     "evalue": "UndefVarError: `scatter` not defined",
     "output_type": "error",
     "traceback": [
      "UndefVarError: `scatter` not defined\n",
      "\n",
      "Stacktrace:\n",
      " [1] top-level scope\n",
      "   @ c:\\xampp\\htdocs\\SGCV\\validation_capr.ipynb:4"
     ]
    }
   ],
   "source": [
    "xlims_ = 1, 350\n",
    "switch_line = 52\n",
    "shift = 0.1\n",
    "scatter(series, background_color=:grey95, xlabel=:Day, ylabel=\"Price \\$\", label=:APPL, markershape=:xcross, markersize=2, markeralpha=1.0, color=:black)\n",
    "plot!(mx, ribbon=sqrt.(vx), label=false)\n",
    "scatter!(switch_line*ones(length(ms)), color=argmax.(ms), markersize=5.0, markeralpha=0.5, label=false)\n",
    "xlims!(xlims_)\n",
    "ylims!(mx[xlims_[1]]-shift, mx[xlims_[2]]+shift)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "ename": "UndefVarError",
     "evalue": "UndefVarError: `results` not defined",
     "output_type": "error",
     "traceback": [
      "UndefVarError: `results` not defined\n",
      "\n",
      "Stacktrace:\n",
      " [1] top-level scope\n",
      "   @ c:\\xampp\\htdocs\\SGCV\\validation_capr.ipynb:2"
     ]
    }
   ],
   "source": [
    "using JLD\n",
    "JLD.save(\"dump/results_3shgf_stocks_mixture.jld\",\"results\",results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "ename": "UndefVarError",
     "evalue": "UndefVarError: `scatter` not defined",
     "output_type": "error",
     "traceback": [
      "UndefVarError: `scatter` not defined\n",
      "\n",
      "Stacktrace:\n",
      " [1] top-level scope\n",
      "   @ c:\\xampp\\htdocs\\SGCV\\validation_capr.ipynb:4"
     ]
    }
   ],
   "source": [
    "xlims_ = 1, 350\n",
    "switch_line = 52\n",
    "shift = 30.0\n",
    "scatter(series, background_color=:grey95, xlabel=:Day, ylabel=\"Price \\$\", label=:APPL, markershape=:xcross, markersize=2, markeralpha=1.0, color=:black)\n",
    "plot!(mx1, ribbon=sqrt.(vx1), label=false)\n",
    "scatter!(switch_line*ones(length(ms1)), color=argmax.(ms1), markersize=5.0, markeralpha=0.5, label=false)\n",
    "xlims!(xlims_)\n",
    "ylims!(mx1[xlims_[1]]-shift, mx1[xlims_[2]]+shift)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2-L HGF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "prod! (generic function with 97 methods)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# NOTE: To run the experiments for the regular HGF you must have to kill Julia first (this happens due to compatibility issues)\n",
    "# Then start running the code from this block\n",
    "using ForneyLab\n",
    "include(\"gcv/GCV.jl\")\n",
    "using Main.GCV\n",
    "using CSV\n",
    "using DataFrames\n",
    "using Plots\n",
    "using ProgressMeter\n",
    "include(\"compatibility.jl\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "pad (generic function with 1 method)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Auxialry function for graph initializtion\n",
    "pad(sym::Symbol, t::Int) = sym*:_*Symbol(lpad(t,3,'0')) # Left-pads a number with zeros, converts it to symbol and appends to sym"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "using CSV, DataFrames\n",
    "df = CSV.File(\"data/AAPL.csv\") |> DataFrame\n",
    "series = df[!, :Open];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "mp (generic function with 1 method)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Function for generating 2-L HGF and message-passing (MP) algorithm\n",
    "function generate_mp(n_samples)\n",
    "    fg = FactorGraph()\n",
    "    z = Vector{Variable}(undef, n_samples)\n",
    "    x = Vector{Variable}(undef, n_samples)\n",
    "    y = Vector{Variable}(undef, n_samples)\n",
    "\n",
    "    @RV ω ~ GaussianMeanPrecision(placeholder(:mω), placeholder(:wω))\n",
    "    @RV [id=pad(:z,1)] z[1] ~ GaussianMeanPrecision(placeholder(:mz_prior1), placeholder(:wz_prior1))\n",
    "    @RV [id=pad(:x,1)] x[1] ~ GaussianMeanPrecision(placeholder(:mx_prior1), placeholder(:wx_prior1))\n",
    "    @RV [id=pad(:y,1)] y[1] ~ GaussianMeanPrecision(x[1], placeholder(:wy_prior1))\n",
    "    placeholder(y[1], :y, index = 1)\n",
    "    for t in 2:n_samples\n",
    "        @RV [id=pad(:z,t)] z[t] ~ GaussianMeanPrecision(z[t - 1], placeholder(pad(:wz_transition, t)))\n",
    "        @RV [id=pad(:x,t)] x[t] ~ GaussianControlledVariance(x[t - 1], z[t], 1.0, ω)\n",
    "        @RV [id=pad(:y,t)] y[t] ~ GaussianMeanPrecision(x[t], placeholder(pad(:wy_transition, t)))\n",
    "        placeholder(y[t], :y, index = t)\n",
    "    end\n",
    "    q = PosteriorFactorization(x, z, ω, ids=[:X :Z :Ω])\n",
    "    algo = messagePassingAlgorithm(free_energy=true)\n",
    "    src_code = algorithmSourceCode(algo, free_energy=true);\n",
    "    return src_code\n",
    "end\n",
    "\n",
    "# Function for running MP-based inference for 2-L HGF\n",
    "function mp(obs;\n",
    "    n_its = 100,\n",
    "    wy_prior1 = 1.0,\n",
    "    ω_m_prior = omega,\n",
    "    ω_w_prior = 1.0,\n",
    "    z_m_prior = 0.0,\n",
    "    z_w_prior = 100.0,\n",
    "    x_m_prior = 0.0,\n",
    "    x_w_prior = 1.0,\n",
    "    x_x_m_prior = zeros(2),\n",
    "    x_x_w_prior = 1.0*diageye(2),\n",
    "    z_z_m_prior = zeros(2),\n",
    "    z_z_w_prior = 100.0*diageye(2),\n",
    "    z_w_transition_prior = 1000.0,\n",
    "    y_w_transition_prior =  1/mnv,\n",
    ")\n",
    "    n_samples = length(obs)\n",
    "    marginals = Dict()\n",
    "    marginals[:ω] = ProbabilityDistribution(ForneyLab.Univariate, GaussianMeanPrecision, m = ω_m_prior, w = ω_w_prior)\n",
    "    marginals[pad(:z,1)] = vague(GaussianMeanPrecision)\n",
    "    marginals[pad(:x,1)] = vague(GaussianMeanPrecision)\n",
    "    for t = 2:n_samples\n",
    "        marginals[pad(:z,t)] = ProbabilityDistribution(ForneyLab.Univariate, GaussianMeanPrecision, m = z_m_prior, w = z_w_prior)\n",
    "        marginals[pad(:x,t)] = ProbabilityDistribution(ForneyLab.Univariate, GaussianMeanPrecision, m = x_m_prior, w = x_w_prior)\n",
    "        marginals[pad(:z,t)*:_*pad(:z,t-1)] = ProbabilityDistribution(ForneyLab.Multivariate,GaussianMeanPrecision, m = z_z_m_prior, w = z_z_w_prior)\n",
    "        marginals[pad(:x,t)*:_*pad(:x,t-1)] = ProbabilityDistribution(ForneyLab.Multivariate,GaussianMeanPrecision, m = x_x_m_prior, w = x_x_w_prior)\n",
    "    end\n",
    "    data = Dict()\n",
    "    data[:y] = obs\n",
    "    data[:mz_prior1] = z_m_prior\n",
    "    data[:wz_prior1] = z_w_prior\n",
    "    data[:mx_prior1] = x_m_prior\n",
    "    data[:wx_prior1] = x_w_prior\n",
    "    data[:wy_prior1] = wy_prior1\n",
    "    data[:mω] = ω_m_prior\n",
    "    data[:wω] = ω_w_prior\n",
    "    for t = 1:n_samples\n",
    "        data[pad(:wz_transition, t)] = z_w_transition_prior\n",
    "        data[pad(:wy_transition, t)] = y_w_transition_prior\n",
    "    end\n",
    "\n",
    "\n",
    "    fe = Vector{Float64}(undef, n_its)\n",
    "\n",
    "    @showprogress \"Iterations\" for i = 1:n_its\n",
    "\n",
    "        stepX!(data, marginals)\n",
    "        stepZ!(data, marginals)\n",
    "        stepΩ!(data, marginals)\n",
    "\n",
    "        fe[i] = freeEnergy(data, marginals)\n",
    "    end\n",
    "\n",
    "    mz = [ForneyLab.unsafeMean(marginals[pad(:z,t)]) for t=1:n_samples]\n",
    "    vz = [ForneyLab.unsafeVar(marginals[pad(:z,t)]) for t=1:n_samples]\n",
    "    mω = ForneyLab.unsafeMean(marginals[:ω])\n",
    "    vω = ForneyLab.unsafeCov(marginals[:ω])\n",
    "    mx = [ForneyLab.unsafeMean(marginals[pad(:x,t)]) for t=1:n_samples]\n",
    "    vx = [ForneyLab.unsafeVar(marginals[pad(:x,t)]) for t=1:n_samples]\n",
    "    return mz,vz,mω, vω, mx,vx,fe\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "ename": "UndefVarError",
     "evalue": "UndefVarError: `results` not defined",
     "output_type": "error",
     "traceback": [
      "UndefVarError: `results` not defined\n",
      "\n",
      "Stacktrace:\n",
      " [1] top-level scope\n",
      "   @ c:\\xampp\\htdocs\\SGCV\\validation_capr.ipynb:2"
     ]
    }
   ],
   "source": [
    "using JLD\n",
    "JLD.save(\"dump/results_hgf_stocks_mixture.jld\", \"results\", results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "ename": "ErrorException",
     "evalue": "Multiple applicable SumProductRule{Equality} for interface :3 of equ_z_002_1 with inbound types:\n[Message{GaussianMeanVariance}, Message{Function}, Nothing]\n[SPEqualityFnG, Main.GCV.SPEqualityGaussianGCV]",
     "output_type": "error",
     "traceback": [
      "Multiple applicable SumProductRule{Equality} for interface :3 of equ_z_002_1 with inbound types:\n",
      "[Message{GaussianMeanVariance}, Message{Function}, Nothing]\n",
      "[SPEqualityFnG, Main.GCV.SPEqualityGaussianGCV]\n",
      "\n",
      "Stacktrace:\n",
      "  [1] error(s::String)\n",
      "    @ Base .\\error.jl:35\n",
      "  [2] selectApplicableRule(rule_type::Type, entry::ScheduleEntry, inbound_types::Vector{Type}, applicable_rules::Vector{Type})\n",
      "    @ ForneyLab C:\\Users\\mkpap\\.julia\\packages\\ForneyLab\\EnMB2\\src\\algorithms\\inference_algorithm.jl:139\n",
      "  [3] inferUpdateRule!(entry::ScheduleEntry, rule_type::Type{SumProductRule{Equality}}, inferred_outbound_types::Dict{Union{Nothing, Interface}, Type})\n",
      "    @ ForneyLab C:\\Users\\mkpap\\.julia\\packages\\ForneyLab\\EnMB2\\src\\algorithms\\sum_product.jl:73\n",
      "  [4] inferUpdateRules!(schedule::Vector{ScheduleEntry}; inferred_outbound_types::Dict{Union{Nothing, Interface}, Type})\n",
      "    @ ForneyLab C:\\Users\\mkpap\\.julia\\packages\\ForneyLab\\EnMB2\\src\\message_passing.jl:211\n",
      "  [5] inferUpdateRules!\n",
      "    @ C:\\Users\\mkpap\\.julia\\packages\\ForneyLab\\EnMB2\\src\\message_passing.jl:205 [inlined]\n",
      "  [6] messagePassingSchedule(pf::PosteriorFactor)\n",
      "    @ ForneyLab C:\\Users\\mkpap\\.julia\\packages\\ForneyLab\\EnMB2\\src\\algorithms\\posterior_factor.jl:75\n",
      "  [7] messagePassingAlgorithm(target_variables::Vector{Variable}, pfz::PosteriorFactorization; ep_sites::Vector{Tuple}, id::Symbol, free_energy::Bool)\n",
      "    @ ForneyLab C:\\Users\\mkpap\\.julia\\packages\\ForneyLab\\EnMB2\\src\\algorithms\\inference_algorithm.jl:82\n",
      "  [8] messagePassingAlgorithm\n",
      "    @ C:\\Users\\mkpap\\.julia\\packages\\ForneyLab\\EnMB2\\src\\algorithms\\inference_algorithm.jl:52 [inlined]\n",
      "  [9] #messagePassingAlgorithm#248\n",
      "    @ C:\\Users\\mkpap\\.julia\\packages\\ForneyLab\\EnMB2\\src\\algorithms\\inference_algorithm.jl:108 [inlined]\n",
      " [10] messagePassingAlgorithm\n",
      "    @ C:\\Users\\mkpap\\.julia\\packages\\ForneyLab\\EnMB2\\src\\algorithms\\inference_algorithm.jl:108 [inlined]\n",
      " [11] generate_mp(n_samples::Int64)\n",
      "    @ Main c:\\xampp\\htdocs\\SGCV\\validation_capr.ipynb:20\n",
      " [12] top-level scope\n",
      "    @ c:\\xampp\\htdocs\\SGCV\\validation_capr.ipynb:3"
     ]
    }
   ],
   "source": [
    "omega  = 1.0\n",
    "kappa = 1.0\n",
    "src_code = generate_mp(length(series))\n",
    "eval(Meta.parse(src_code));\n",
    "mz0, vz0, mω0, vω0, mx0,vx0,fe0 = mp(series, x_m_prior=series[1], z_z_w_prior=100.0 * diageye(2), z_w_transition_prior=100.0, ω_m_prior=omega, ω_w_prior=1.0, y_w_transition_prior=1.0);\n",
    "\n",
    "results =   Dict(\"mz\" => mz0, \"vz\" => vz0,\n",
    "                 \"mx\" => mx0, \"vx\" => vx0,\n",
    "                 \"fe\" => fe0,\n",
    "                 \"mω\" => mω0, \"vω\" => vω0,\n",
    "                 \"ωprior\" => omega)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "ename": "ArgumentError",
     "evalue": "ArgumentError: unsupported symbol APPL passed to `label`",
     "output_type": "error",
     "traceback": [
      "ArgumentError: unsupported symbol APPL passed to `label`\n",
      "\n",
      "Stacktrace:\n",
      "  [1] label_to_string\n",
      "    @ C:\\Users\\mkpap\\.julia\\packages\\Plots\\3BCH5\\src\\args.jl:2027 [inlined]\n",
      "  [2] _broadcast_getindex_evalf\n",
      "    @ .\\broadcast.jl:683 [inlined]\n",
      "  [3] _broadcast_getindex\n",
      "    @ .\\broadcast.jl:656 [inlined]\n",
      "  [4] getindex\n",
      "    @ .\\broadcast.jl:610 [inlined]\n",
      "  [5] copy\n",
      "    @ .\\broadcast.jl:888 [inlined]\n",
      "  [6] materialize(bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{0}, Nothing, typeof(Plots.label_to_string), Tuple{Base.RefValue{Symbol}, Int64}})\n",
      "    @ Base.Broadcast .\\broadcast.jl:873\n",
      "  [7] _update_series_attributes!(plotattributes::RecipesPipeline.DefaultsDict, plt::Plots.Plot{Plots.GRBackend}, sp::Plots.Subplot{Plots.GRBackend})\n",
      "    @ Plots C:\\Users\\mkpap\\.julia\\packages\\Plots\\3BCH5\\src\\args.jl:2111\n",
      "  [8] add_series!(plt::Plots.Plot{Plots.GRBackend}, plotattributes::RecipesPipeline.DefaultsDict)\n",
      "    @ Plots C:\\Users\\mkpap\\.julia\\packages\\Plots\\3BCH5\\src\\pipeline.jl:377\n",
      "  [9] _process_seriesrecipe(plt::Any, plotattributes::Any)\n",
      "    @ RecipesPipeline C:\\Users\\mkpap\\.julia\\packages\\RecipesPipeline\\BGM3l\\src\\series_recipe.jl:46\n",
      " [10] _process_seriesrecipes!(plt::Any, kw_list::Any)\n",
      "    @ RecipesPipeline C:\\Users\\mkpap\\.julia\\packages\\RecipesPipeline\\BGM3l\\src\\series_recipe.jl:27\n",
      " [11] recipe_pipeline!(plt::Any, plotattributes::Any, args::Any)\n",
      "    @ RecipesPipeline C:\\Users\\mkpap\\.julia\\packages\\RecipesPipeline\\BGM3l\\src\\RecipesPipeline.jl:99\n",
      " [12] _plot!(plt::Plots.Plot, plotattributes::Any, args::Any)\n",
      "    @ Plots C:\\Users\\mkpap\\.julia\\packages\\Plots\\3BCH5\\src\\plot.jl:223\n",
      " [13] plot(args::Any; kw::Base.Pairs{Symbol, V, Tuple{Vararg{Symbol, N}}, NamedTuple{names, T}} where {V, N, names, T<:Tuple{Vararg{Any, N}}})\n",
      "    @ Plots C:\\Users\\mkpap\\.julia\\packages\\Plots\\3BCH5\\src\\plot.jl:102\n",
      " [14] top-level scope\n",
      "    @ c:\\xampp\\htdocs\\SGCV\\validation_capr.ipynb:1"
     ]
    }
   ],
   "source": [
    "plot(series, background_color=:grey95, xlabel=:Day, ylabel=\"Price \\$\", label=:APPL)\n",
    "plot!(mx0, ribbon=sqrt.(vx0), label=false)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "ename": "ArgumentError",
     "evalue": "ArgumentError: unsupported symbol APPL passed to `label`",
     "output_type": "error",
     "traceback": [
      "ArgumentError: unsupported symbol APPL passed to `label`\n",
      "\n",
      "Stacktrace:\n",
      "  [1] label_to_string\n",
      "    @ C:\\Users\\mkpap\\.julia\\packages\\Plots\\3BCH5\\src\\args.jl:2027 [inlined]\n",
      "  [2] _broadcast_getindex_evalf\n",
      "    @ .\\broadcast.jl:683 [inlined]\n",
      "  [3] _broadcast_getindex\n",
      "    @ .\\broadcast.jl:656 [inlined]\n",
      "  [4] getindex\n",
      "    @ .\\broadcast.jl:610 [inlined]\n",
      "  [5] copy\n",
      "    @ .\\broadcast.jl:888 [inlined]\n",
      "  [6] materialize(bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{0}, Nothing, typeof(Plots.label_to_string), Tuple{Base.RefValue{Symbol}, Int64}})\n",
      "    @ Base.Broadcast .\\broadcast.jl:873\n",
      "  [7] _update_series_attributes!(plotattributes::RecipesPipeline.DefaultsDict, plt::Plots.Plot{Plots.GRBackend}, sp::Plots.Subplot{Plots.GRBackend})\n",
      "    @ Plots C:\\Users\\mkpap\\.julia\\packages\\Plots\\3BCH5\\src\\args.jl:2111\n",
      "  [8] add_series!(plt::Plots.Plot{Plots.GRBackend}, plotattributes::RecipesPipeline.DefaultsDict)\n",
      "    @ Plots C:\\Users\\mkpap\\.julia\\packages\\Plots\\3BCH5\\src\\pipeline.jl:377\n",
      "  [9] _process_seriesrecipe(plt::Any, plotattributes::Any)\n",
      "    @ RecipesPipeline C:\\Users\\mkpap\\.julia\\packages\\RecipesPipeline\\BGM3l\\src\\series_recipe.jl:46\n",
      " [10] _process_seriesrecipes!(plt::Any, kw_list::Any)\n",
      "    @ RecipesPipeline C:\\Users\\mkpap\\.julia\\packages\\RecipesPipeline\\BGM3l\\src\\series_recipe.jl:27\n",
      " [11] recipe_pipeline!(plt::Any, plotattributes::Any, args::Any)\n",
      "    @ RecipesPipeline C:\\Users\\mkpap\\.julia\\packages\\RecipesPipeline\\BGM3l\\src\\RecipesPipeline.jl:99\n",
      " [12] _plot!(plt::Plots.Plot, plotattributes::Any, args::Any)\n",
      "    @ Plots C:\\Users\\mkpap\\.julia\\packages\\Plots\\3BCH5\\src\\plot.jl:223\n",
      " [13] plot(args::Any; kw::Base.Pairs{Symbol, V, Tuple{Vararg{Symbol, N}}, NamedTuple{names, T}} where {V, N, names, T<:Tuple{Vararg{Any, N}}})\n",
      "    @ Plots C:\\Users\\mkpap\\.julia\\packages\\Plots\\3BCH5\\src\\plot.jl:102\n",
      " [14] scatter(args::Any; kw::Base.Pairs{Symbol, V, Tuple{Vararg{Symbol, N}}, NamedTuple{names, T}} where {V, N, names, T<:Tuple{Vararg{Any, N}}})\n",
      "    @ Plots C:\\Users\\mkpap\\.julia\\packages\\RecipesBase\\BRe07\\src\\RecipesBase.jl:427\n",
      " [15] macro expansion\n",
      "    @ c:\\xampp\\htdocs\\SGCV\\validation_capr.ipynb:3 [inlined]\n",
      " [16] top-level scope\n",
      "    @ C:\\Users\\mkpap\\.julia\\packages\\Plots\\3BCH5\\src\\animation.jl:251"
     ]
    }
   ],
   "source": [
    "shift = 10\n",
    "@gif for i in 1:length(series)\n",
    "    scatter(series[1:i], markershape=:xcross, markersize=1.5, markeralpha=0.5, xlim = (0, length(series)), ylim = (minimum(series)-shift, maximum(series)+shift), label=:APPL, xlabel=:Day, ylabel=\"Price \\$\", background_color=:grey95)\n",
    "#     scatter(series[1:i], markershape=:xcross, markersize=1.5, markeralpha=0.5, xlim = (0, length(series)), ylim = (minimum(series)-shift, maximum(series)+shift), label=:APPL, xlabel=:Day, ylabel=\"Price \\$\", background_color=:grey95)\n",
    "    plot!(mx0[1:i], ribbon=sqrt.(vx0[1:i]), label=false)\n",
    "end every 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "ename": "UndefVarError",
     "evalue": "UndefVarError: `mz0` not defined",
     "output_type": "error",
     "traceback": [
      "UndefVarError: `mz0` not defined\n",
      "\n",
      "Stacktrace:\n",
      " [1] macro expansion\n",
      "   @ c:\\xampp\\htdocs\\SGCV\\validation_capr.ipynb:3 [inlined]\n",
      " [2] top-level scope\n",
      "   @ C:\\Users\\mkpap\\.julia\\packages\\Plots\\3BCH5\\src\\animation.jl:251"
     ]
    }
   ],
   "source": [
    "shift = 1\n",
    "@gif for i in 1:length(series)\n",
    "    plot(mz0[1:i], ribbon=sqrt.(vz0[1:i]), label=false, xlim = (0, length(series)), ylim = (minimum(mz0)-shift, maximum(mz0)+shift) ,background_color=:grey95)\n",
    "end every 1"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.9.3",
   "language": "julia",
   "name": "julia-1.9"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.9.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}