{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "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.Web\n",
    "import Data.Manifold.DifferentialEquation\n",
    "import Math.LinearMap.Category\n",
    "import Data.VectorSpace\n",
    "import Data.Basis (Basis)\n",
    "import Linear(V2(..), ex, ey)\n",
    "import Data.Semigroup\n",
    "import qualified Data.Foldable as Hask\n",
    "import Control.Lens\n",
    ":opt no-lint\n",
    "import Control.Category.Constrained.Prelude\n",
    "import Control.Arrow.Constrained"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From [diagrams](http://projects.haskell.org/diagrams):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import Diagrams.Prelude (p2, circle, (&), (^&), moveTo, opacity, fromVertices, Point(P))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "type X = ℝ\n",
    "type T = ℝ\n",
    "type U = ℝ\n",
    "type Ðx'U = ℝ\n",
    "type x × y = ℝ²\n",
    "et = ey :: Basis ℝ²"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From [dynamic-plot](http://hackage.haskell.org/package/dynamic-plot):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import Graphics.Dynamic.Plot.R2\n",
    "import Data.Colour.Names\n",
    "import Data.Colour.Manifold\n",
    "\n",
    "colourscheme :: Shade' ℝ -> Shade (Colour ℝ)\n",
    "colourscheme (Shade' u du) = interp (Shade u $ dualNorm du :: Shade ℝ)\n",
    " where Just interp = rangeOnGeodesic darkblue orange\n",
    " \n",
    "prettyWebPlot :: PointsWeb ℝ² y -> DynamicPlottable\n",
    "prettyWebPlot w = plot [ diagramPlot . opacity 0.5 $ fromVertices [P r₁, P r₂]\n",
    "                       | ((r₁@(V2 x₁ y₁),_),(r₂@(V2 x₂ y₂),_)) <- edg ]\n",
    " where edg = webEdges w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "colourscheme_heat :: Shade' (U, Ðx'U) -> Shade (Colour ℝ)\n",
    "colourscheme_heat = cm . factoriseShade\n",
    " where cm :: (Shade' U, Shade' Ðx'U) -> Shade (Colour ℝ)\n",
    "       cm (Shade' u du, _) = interp (Shade u $ dualNorm du :: Shade ℝ)\n",
    "       Just interp = rangeOnGeodesic darkblue orange"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "μ :: LocalLinear (X × T) (U, Ðx'U) +> (U, Ðx'U)\n",
    "μ = arr.LinearFunction $\n",
    "        \\(LinearMap (V2 (ðx'u, ðx'ðx'u)\n",
    "                        (ðt'u, ðt'ðx'u))) -> (ðx'ðx'u - ðt'u, ðx'u)\n",
    "\n",
    "heatEqn :: DifferentialEqn (X × T) (U, Ðx'U)\n",
    "      -- Shade (X×T, U,Ðx'U) -> Shade' (X×T +> (U,Ðx'U)) -> ?(Shade' (X×T +> (U,Ðx'U)))\n",
    "heatEqn = constLinearPDE $ μ"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Option {getOption = Just (Shade' {_shade'Ctr = (0.0,0.0), _shade'Narrowness = spanNorm [(0.0,1.0)]},Shade' {_shade'Ctr = Left ().<V2 0.0 0.3333333333333333 ^+^ Right ().<V2 0.6666666666666667 0.0, _shade'Narrowness = spanNorm [ex.⊗(1.3107943462579543,0.0) ^+^ ey.⊗(0.0,0.0),ex.⊗(0.0,1.3107943462579543) ^+^ ey.⊗(-0.6866065623255952,0.0),ex.⊗(0.0,0.0) ^+^ ey.⊗(1.116581052478165,0.0),ex.⊗(0.0,0.0) ^+^ ey.⊗(0.0,0.904534033733291)]})}"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "--(`linMapAt` V2 0 2) <$>\n",
    "heatEqn ((V2 0 0,(0,0)):±[(V2 1 0,(0,0)), (V2 0 1,(0,0)), (V2 0 0,(1,0)), (V2 0 0,(0,1))])\n",
    "        (Shade' (LinearMap $ V2 (0,1) (0,0)) . spanNorm $ Tensor <$>\n",
    "                         [V2 (1,0) (0,0), V2 (0,1) (0,0), V2 (0,0) (1,0), V2 (0,0) (0,1)] )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "GraphWindowSpecR2{lBound=-1.333333333333333, rBound=1.333333333333333, bBound=-1.0328396847046188, tBound=1.0328396847046166, xResolution=640, yResolution=480}"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "initState_heat :: X -> (U, Ðx'U)\n",
    "initState_heat x = ( tanh (s * (1 - x^2))\n",
    "                   , - 2 * s*x / cosh (s * (1 - x^2))^2 )\n",
    " where s = 0.5\n",
    "\n",
    "-- plotWindow $ continFnPlot <$> [fst . initState_heat, snd . initState_heat]\n",
    "\n",
    "tf_heat :: Needle X -> Needle T -> PointsWeb (X × T) (Shade' (U, Ðx'U))\n",
    "tf_heat δx₀ δt₀ = fromWebNodes euclideanMetric\n",
    "      $  [(V2 x 0, initState_heat x|±|[(0.1,0), (0,0.1)]) | x<-[-2, δx₀-2 .. 0] ]\n",
    "      ++ [(V2 x t, zeroV|±|[(1, 0),  (0,1)  ]) | x<-[-2, δx₀-2 .. 0], t<-[δt₀, δt₀*2 .. 1] ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "startSt_heat = tf_heat 0.13 0.13\n",
    "forM_ [ iterateFilterDEqn_static (HighlightInconsistencies $ (-1, 0)|±|[(10, 0),(0, 10)]) heatEqn startSt_heat ]\n",
    "  $ \\tfs ->\n",
    "    plotWindow\n",
    "       [ prettyWebPlot $ head tfs\n",
    "       , plotLatest [ plot (fmap colourscheme_heat tfi)\n",
    "                       & legendName (\"i = \"++show i)\n",
    "                    | (i,tfi) <- zip [0..] tfs ] ]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Initial state from which to start refining solution of the heat PDE](https://raw.githubusercontent.com/leftaroundabout/manifolds/master/manifolds/images/examples/PDE-solution-filter/HeatEqn-InitialState.png)"
   ]
  },
  {
   "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": "7.10.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}