{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "{-# LANGUAGE TypeOperators, FlexibleContexts, TypeFamilies #-}\n",
    "import Prelude ()\n",
    "import Data.Manifold.TreeCover\n",
    "import Data.Random\n",
    "import Data.Random.Manifold\n",
    "import Data.Manifold\n",
    "import Data.Manifold.Shade\n",
    "import Data.Manifold.Web\n",
    "import Data.Manifold.DifferentialEquation\n",
    "import Math.LinearMap.Category\n",
    "import Data.Function.Affine\n",
    "import Data.VectorSpace\n",
    "import Linear(V2(..), _x, _y)\n",
    "import Data.Semigroup\n",
    "import qualified Data.Foldable as Hask\n",
    "import qualified Control.Monad as Hask\n",
    "import Control.Lens\n",
    "\n",
    "import Control.Category.Constrained.Prelude\n",
    "import Control.Arrow.Constrained\n",
    ":opt no-lint  -- lint gives bogus warnings with constrained-categories"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From [diagrams](http://projects.haskell.org/diagrams):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import Diagrams.Prelude (p2, circle, (&), (^&), moveTo, opacity, fromVertices)\n",
    "import qualified Diagrams.Prelude as Dia"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From [dynamic-plot](http://hackage.haskell.org/package/dynamic-plot):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import Graphics.Dynamic.Plot.R2\n",
    "\n",
    "type T = ℝ\n",
    "type X = ℝ\n",
    "\n",
    "viewRange = plot [forceXRange (-2,4), forceYRange (-1,3)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "μ :: LocalLinear T X +> X\n",
    "μ = arr.LinearFunction $ \\d -> d$1\n",
    "\n",
    "deq :: ODE T X\n",
    "deq = constLinearODE μ"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "tf :: Needle T -> PointsWeb T (Shade' X)\n",
    "tf δt₀ = fromWebNodes euclideanMetric $ \n",
    "        [] --  [(t, exp t|±|[0.001]) | t<-take 2 [-δt₀, -2*δt₀]]\n",
    "       ++ (0, 1|±|[0.0001])\n",
    "        : [(t, 0|±|[3]) | t<-[δt₀, 2*δt₀ .. 1.2] ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "forM_ [ Hask.toList $ iterateFilterDEqn_static (inconsistencyAware intersectShade's) id deq (tf 0.1)\n",
    "      , Hask.toList $ iterateFilterDEqn_static_selective\n",
    "                          (inconsistencyAware intersectShade's) id (euclideanVolGoal 0.01) deq (tf 0.05)\n",
    "      , iterateFilterDEqn_adaptive euclideanMetric AbortOnInconsistency\n",
    "                   deq (euclideanVolGoal 0.01) (tf 0.2)]\n",
    "  $ \\tfs -> do\n",
    "    plotWindow\n",
    "       [ plot ((1^&exp 1) :± [0.1^&0, 0^&0.1] :: Shade ℝ²) -- Euler's number as reference for x=1\n",
    "       , plotLatest [plot st & legendName (show i) | (i,st) <- zip [0..] tfs]\n",
    "         & plotDelay 0.5\n",
    "       , plot ((0.4^&exp 0.4) :± [0.1^&0, 0^&0.1] :: Shade ℝ²)\n",
    "       , xAxisLabel \"𝑡\", yAxisLabel \"exp 𝑡\"\n",
    "       ]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Static resolution:\n",
    "![filtering the exponential function as a solution to an ODE](https://raw.githubusercontent.com/leftaroundabout/manifolds/master/manifolds/images/examples/ODE-solution-filter/simple-exponential.gif)\n",
    "Adaptive resolution:\n",
    "![filtering the exponential function as a solution to an ODE (adaptive resolution)](https://raw.githubusercontent.com/leftaroundabout/manifolds/master/manifolds/images/examples/ODE-solution-filter/simple-exponential-adaptive.gif)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import qualified Control.Comonad.Cofree as Cofree\n",
    "import Data.Foldable(Foldable)\n",
    "import Data.List.NonEmpty(NonEmpty(..))\n",
    "import Control.Monad.Trans.Except\n",
    "\n",
    "tf_bad :: Needle T -> PointsWeb T (Shade' X)\n",
    "tf_bad δt₀ = fromWebNodes euclideanMetric $ \n",
    "        (0, 0.5|±|[0.003]) : [(t, (1-t)|±|[1]) | t<-[δt₀, 2*δt₀ .. 1.2] ]\n",
    "\n",
    "tfs_inconsistent = iterateFilterDEqn_pathwise (indicateInconsistencies intersectShade's)\n",
    "                     id deq (tf_bad 0.05)\n",
    "plotWindow [ plotLatest [plot st & legendName (show i) | (i,st) <- zip [0..] $ Hask.toList tfs_inconsistent]\n",
    "              -- & plotDelay 2\n",
    "           ]\n",
    "\n",
    "findErr :: (Hask.Monad m, Hask.Foldable m)\n",
    "          => Cofree.Cofree m a -> m ()\n",
    "findErr (a Cofree.:< q) = case Hask.toList q of\n",
    "        [] -> const () <$> q\n",
    "        l  -> foldr1 (>>) $ findErr<$>l\n",
    "\n",
    "case runExcept $ findErr tfs_inconsistent of Left e -> print e"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "type Y = ℝ\n",
    "type XY = ℝ²\n",
    "\n",
    "μ₂ :: LocalLinear T XY +> XY\n",
    "μ₂ = arr.LinearFunction $ ($1) >>> \\(V2 dx dy) -> V2 dy (-dx)\n",
    "\n",
    "deq₂ :: ODE T XY\n",
    "deq₂ = constLinearODE μ₂"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import Data.Manifold.Function.LocalModel\n",
    "\n",
    "tf₂ :: Needle T -> PointsWeb T (Shade' XY)\n",
    "tf₂ δt₀ = fromWebNodes euclideanMetric $\n",
    "        [ (t, (1^&0)|±|[0^&δ, δ^&0])\n",
    "        | t<-(**γ)<$>[0, δt₀**(1/γ) .. tEnd**(1/γ)]\n",
    "        , let δ = 1e-4 + 4*tanh t ]\n",
    " where γ = 1.5\n",
    "       tEnd = 15"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "forM_ [Hask.toList $ iterateFilterDEqn_pathwise (inconsistencyAware intersectShade's) id deq₂ (tf₂ 0.03)]\n",
    "   $ \\tfs₂ -> do\n",
    " plotWindow [ continFnPlot sin, continFnPlot cos\n",
    "  , plotLatest [ plotMultiple [ plot [fmap fst tffacts, fmap snd tffacts] & legendName (show i) & tint Dia.green\n",
    "                              , plot (take 0 [fmap fst tffactsCG, fmap snd tffactsCG]) & tint Dia.purple ]\n",
    "                   & plotDelay 1\n",
    "               | (i,tf') <- zip [0..] tfs₂\n",
    "               , let tffacts = fmap factoriseShade (coerceShade <$> tf' :: PointsWeb T (Shade' (X,Y)))\n",
    "                     tffactsCG = fromWebNodes euclideanMetric\n",
    "                               $ second (factoriseShade . coerceShade . fst . quadraticModel_derivatives )\n",
    "                               <$> localModels_CGrid tf'\n",
    "                                 :: PointsWeb T (Shade' X, Shade' Y)\n",
    "               ]\n",
    "            , xAxisLabel \"𝑡\", yAxisLabel \"cos 𝑡, sin 𝑡\"\n",
    "            , forceXRange (-0.01,6) ]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "![filtering the sin and cosine function as a solution to an ODE](https://raw.githubusercontent.com/leftaroundabout/manifolds/master/manifolds/images/examples/ODE-solution-filter/simple-trigonometric.gif)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "--case runExcept $ findErr (iterateFilterDEqn_static (indicateInconsistencies intersectShade's) id deq₂ (tf₂ 0.004)) of\n",
    "  -- Left e -> print e\n",
    "\n",
    "--[(x, cos x) | x <- (^2)<$>[0, sqrt 0.004 .. sqrt 0.1]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import Data.Manifold.Function.LocalModel\n",
    "import Math.LinearMap.Category.Derivatives\n",
    "\n",
    "deq₃Init :: T -> Shade' XY\n",
    "deq₃Init t = (f t ^& 2)|±|[1e-1^&0, 0^&δ]\n",
    " where δ | t > -4      = 0.8\n",
    "         | otherwise  = 1e-1\n",
    "       f t = sin (t + 3*sin (t/2))\n",
    "\n",
    "-- Simple “differentials equalisation” system: make it so the functions x(t) and y(t) move in parallel.\n",
    "deq₃ :: DifferentialEqn AffineModel T XY\n",
    "deq₃ (Shade (t,_) _) = LocalDifferentialEqn {\n",
    "   _rescanDifferentialEqn = \\(AffineModel d⁰ d¹)\n",
    "      -> let x' = projectShade (lensEmbedding (1*∂_x/∂id)) (dualShade d¹) :: Shade' ℝ\n",
    "             y' = projectShade (lensEmbedding (1*∂_y/∂id)) (dualShade d¹) :: Shade' ℝ\n",
    "         in ( (if t>0 then intersectShade's . (:|[dualShade d⁰]) else return)\n",
    "                $ deq₃Init t\n",
    "            , return $\n",
    "                  embedShade (lensEmbedding (1*∂_y/∂id)) x'\n",
    "              {-:|[ embedShade (lensEmbedding (1*∂_x/∂id)) y' ]-})\n",
    "   }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import Control.Applicative ((<|>))\n",
    "\n",
    "tf₃ :: Needle T -> PointsWeb T (Shade' XY)\n",
    "tf₃ δt₀ = fromWebNodes euclideanMetric\n",
    "        [ (t, deq₃Init t)\n",
    "        | t <- subtract (δt₀/2) . (^2)<$>[0, sqrt δt₀ .. 3] ]\n",
    "\n",
    "tfs₃ = iterateFilterDEqn_static (indicateInconsistencies $ \\l -> intersectShade's l)\n",
    "                     id deq₃ (tf₃ 0.01)\n",
    "tfs₃_s = iterateFilterDEqn_static_selective (indicateInconsistencies $ \\l -> intersectShade's l)\n",
    "                     id (euclideanVolGoal 0.001) deq₃ (tf₃ 0.01)\n",
    "\n",
    "forM_ [ Hask.toList tfs₃, Hask.toList tfs₃_s ]\n",
    "   $ \\tfs -> do\n",
    " plotWindow [\n",
    "    plotLatest [ plot [fmap fst tffacts, fmap snd tffacts] & legendName (show i)\n",
    "               | (i,tf') <- zip [0..] tfs\n",
    "               , let tffacts = fmap factoriseShade (coerceShade <$> tf' :: PointsWeb T (Shade' (X,Y))) ] ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "-- case runExcept $ findErr tfs₃ of Left e -> print e"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Haskell",
   "language": "haskell",
   "name": "haskell"
  },
  "language_info": {
   "codemirror_mode": "ihaskell",
   "file_extension": ".hs",
   "name": "haskell",
   "version": "8.2.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}