{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "{-# LANGUAGE ScopedTypeVariables, TypeOperators #-}\n",
    "import Prelude ()\n",
    "import Control.Category.Constrained.Prelude\n",
    "import Control.Arrow.Constrained\n",
    "\n",
    "import Data.Manifold\n",
    "import Data.Manifold.TreeCover\n",
    "import Data.Manifold.Web\n",
    "import Data.Manifold.Shade\n",
    "import Data.Manifold.Types\n",
    "import Data.Function.Affine\n",
    "import Math.LinearMap.Category\n",
    "import Math.LinearMap.Category.Derivatives\n",
    "import Data.VectorSpace\n",
    "import Linear (V2(..), _x, _y)\n",
    "import Control.Lens\n",
    ":opt no-lint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "test x = x`seq`putStrLn \"Ok\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Ok"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Ok"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Ok"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import Graphics.Dynamic.Plot.R2\n",
    "import Data.Colour.Names\n",
    "import Data.Colour.Manifold\n",
    "import Diagrams.Prelude ((^&))\n",
    "\n",
    "colourscheme :: Shade' ℝ -> Shade (Colour ℝ)\n",
    "colourscheme = f `seq` \\(Shade' y e) -> f (Shade y $ dualNorm e)\n",
    " where Just (f :: Shade ℝ -> Shade (Colour ℝ))\n",
    "          = rangeWithinVertices (1, maxc)\n",
    "                              [ (0, black) ]\n",
    "       Just maxc = toInterior (grey :: Colour ℝ)\n",
    "test colourscheme\n",
    "\n",
    "colourscheme_nabla :: Shade' (V2 ℝ+>ℝ) -> Shade (Colour ℝ)\n",
    "colourscheme_nabla = f `seq` \\(Shade' y e) -> f (Shade y $ dualNorm e)\n",
    " where Just (f :: Shade (V2 ℝ+>ℝ) -> Shade (Colour ℝ))\n",
    "          = rangeWithinVertices (zeroV, neutral)\n",
    "                              [ (V2 2 0-+|>1, red)\n",
    "                              , (V2 0 2-+|>1, green) ]\n",
    "       Just neutral = toInterior (grey :: Colour ℝ)\n",
    "test colourscheme_nabla\n",
    "\n",
    "colourscheme_laplace :: Shade' (V2 ℝ⊗〃+>ℝ) -> Shade (Colour ℝ)\n",
    "colourscheme_laplace = f `seq` \\(Shade' y e) -> f (Shade y $ dualNorm e)\n",
    " where Just (f :: Shade (V2 ℝ⊗〃+>ℝ) -> Shade (Colour ℝ))\n",
    "          = rangeWithinVertices (zeroV, neutral)\n",
    "                              [ ((squareV (V2 1 0)^+^squareV (V2 0 1)) -+|> 1, white)\n",
    "                              , ((squareV (V2 10 10)^-^squareV (V2 10 (-10))) -+|> 1, red)\n",
    "                              , ((squareV (V2 10 0)^-^squareV (V2 0 10)) -+|> 1, blue) ]\n",
    "       Just neutral = toInterior (grey :: Colour ℝ)\n",
    "test colourscheme_laplace"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[3.100000000000001,2.900000000000001,2.7000000000000006,2.5000000000000004,2.3000000000000003,2.1,1.9,1.7000000000000002,0.9000000000000001,0.7000000000000002,1.1,1.3000000000000003,1.5,0.30000000000000004,0.5,0.1]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import Data.Manifold.Function.LocalModel\n",
    "\n",
    "osc :: PointsWeb ℝ (Shade' ℝ)\n",
    "osc = fromWebNodes euclideanMetric\n",
    "           $ [(x, sin x|±|[0.001]) | x<-[0,0.2 .. pi]]\n",
    "fst <$> (localModels_CGrid osc :: [(ℝ, QuadraticModel ℝ ℝ)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "GraphWindowSpecR2{lBound=-0.5333333333333335, rBound=3.7333333333333343, bBound=-1.3363228594363021, tBound=1.3332730976811917, xResolution=640, yResolution=480}"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plotWindow [ plot osc  & legendName \"sin\"\n",
    "           , plot (linIsoTransformShade (lfun ($ 1)) <$> differentiateUncertainWebFunction osc)\n",
    "                       & legendName \"sin' (nodes)\"\n",
    "           , plot (linIsoTransformShade (lfun ($ squareV 1)) <$> differentiate²UncertainWebFunction osc)\n",
    "                       & legendName \"sin'' (nodes)\"\n",
    "           , plot (fromWebNodes euclideanMetric\n",
    "               $ second (linIsoTransformShade (lfun ($ 1)) . fst . snd . quadraticModel_derivatives)\n",
    "                              <$> localModels_CGrid osc )\n",
    "                       & legendName \"sin' (links)\"\n",
    "           , plot (fromWebNodes euclideanMetric\n",
    "               $ second (linIsoTransformShade (lfun ($ squareV 1)) . snd . snd . quadraticModel_derivatives)\n",
    "                              <$> localModels_CGrid osc )\n",
    "                       & legendName \"sin'' (links)\"\n",
    "           , continFnPlot (negate . sin)\n",
    "                       & legendName \"-sin\"\n",
    "           ]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Numerically differentiating a function sampled on a 1D  web](https://raw.githubusercontent.com/leftaroundabout/manifolds/master/manifolds/images/examples/WebSampling/1D-differentiation.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "derivatives1Plot, derivatives2Plot :: PointsWeb ℝ (Shade' ℝ²) -> [DynamicPlottable]\n",
    "derivatives1Plot signal\n",
    "         = [ plot [ projectShade(lensEmbedding _x) . fst <$>signalLM\n",
    "                  , projectShade(lensEmbedding _y) . fst <$>signalLM :: PointsWeb ℝ (Shade' ℝ) ]\n",
    "           , plot [ projectShade(lensEmbedding (1*∂_x/∂id)) . fst . snd <$>signalLM\n",
    "                  , projectShade(lensEmbedding (1*∂_y/∂id)) . fst . snd <$>signalLM :: PointsWeb ℝ (Shade' ℝ) ]\n",
    "                  & legendName \"∂\"\n",
    "           ]\n",
    " where signalLM = fmap quadraticModel_derivatives . fromWebNodes euclideanMetric $ localModels_CGrid signal\n",
    "derivatives2Plot signal\n",
    "         = [ plot [ projectShade(lensEmbedding _x)<$>signal\n",
    "                  , projectShade(lensEmbedding _y)<$>signal :: PointsWeb ℝ (Shade' ℝ) ]\n",
    "           , plot [ projectShade(lensEmbedding (1*∂_x/∂id)) . fst . snd <$>signalLM\n",
    "                  , projectShade(lensEmbedding (1*∂_y/∂id)) . fst . snd <$>signalLM :: PointsWeb ℝ (Shade' ℝ) ]\n",
    "                  & legendName \"∂\"\n",
    "           , plot [ projectShade(lensEmbedding (1*∂_x/∂id.∂id)) . snd . snd <$>signalLM\n",
    "                  , projectShade(lensEmbedding (1*∂_y/∂id.∂id)) . snd . snd <$>signalLM :: PointsWeb ℝ (Shade' ℝ) ]\n",
    "                  & legendName \"∂²\"\n",
    "           ]\n",
    " where signalLM = fmap quadraticModel_derivatives . fromWebNodes euclideanMetric $ localModels_CGrid signal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "GraphWindowSpecR2{lBound=-1.033333333333334, rBound=7.233333333333339, bBound=-8.195384039338357, tBound=8.304633783682512, xResolution=640, yResolution=480}"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "twinSignal :: PointsWeb ℝ (Shade' ℝ²)\n",
    "twinSignal = fromWebNodes euclideanMetric\n",
    "               $ [(x, sin (x+sin x) ^& cos (x-cos x)\n",
    "                      |±|[2e-3^&0, 0^&(5e-3*(x+0.1))]) | x<-[0,0.2 .. 2*pi]]\n",
    "\n",
    "plotWindow $ derivatives2Plot twinSignal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Numerically differentiating a function sampled on a 1D  web](https://raw.githubusercontent.com/leftaroundabout/manifolds/master/manifolds/images/examples/WebSampling/1D-differentiation_uncertain.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "GraphWindowSpecR2{lBound=-0.2999999999999998, rBound=3.2999999999999985, bBound=-1.8175505050505059, tBound=3.722853535353538, xResolution=640, yResolution=480}"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "jumpFunction :: PointsWeb ℝ (Shade' ℝ²)\n",
    "jumpFunction = fromWebNodes euclideanMetric\n",
    "               $ [(x, signum (x-1) ^& signum (x-2)\n",
    "                      |±|[1e-2^&0, 0^&1e-2]) | x<-[0,0.3 .. 3]]\n",
    "plotWindow $ derivatives1Plot jumpFunction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "gaußianPeak :: [ℝ²] -> ℝ -> PointsWeb ℝ² (Shade' ℝ)\n",
    "gaußianPeak ps δy = fromWebNodes euclideanMetric\n",
    "       [ (xy, exp (-2*magnitudeSq xy) |±|[δy * if x > 0 then 10 else 1]) | xy@(V2 x _) <- ps ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "δ = 0.01\n",
    "grid = [V2 x y | x<-[-2,-1.8..2], y<-[-1,-0.8..2]]\n",
    "hexGrid = [V2 x y | x<-[-2,-1.7..2], y<-[-1,-0.6..2]]\n",
    "        ++ [V2 x y | x<-[-1.85,-1.55..2], y<-[-0.8,-0.4..2]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "forM_ [grid{-, hexGrid-}] $ \\g -> do\n",
    " let gaußianPeak_loc = fmap quadraticModel_derivatives\n",
    "                 . fromWebNodes euclideanMetric\n",
    "                 . localModels_CGrid $ gaußianPeak g δ\n",
    "\n",
    " plotWindow [plot . fmap (colourscheme . fst) $ gaußianPeak_loc, dynamicAxes]\n",
    " plotWindow [plot . fmap (colourscheme_nabla . fst . snd) $ gaußianPeak_loc, dynamicAxes]\n",
    " plotWindow [plot . fmap (colourscheme_laplace . snd . snd) $ gaußianPeak_loc, dynamicAxes]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Evaluating a sampled function on the links of a 2D web](https://raw.githubusercontent.com/leftaroundabout/manifolds/master/manifolds/images/examples/WebSampling/noisyGaussianPeak_links.png)\n",
    "![Numerically differentiating (nabla) a sampled function on the links of a 2D  web](https://raw.githubusercontent.com/leftaroundabout/manifolds/master/manifolds/images/examples/WebSampling/noisyGaussianPeak-nabla_links.png)\n",
    "![Numerically differentiating (Laplace) a sampled function on the links of a 2D  web](https://raw.githubusercontent.com/leftaroundabout/manifolds/master/manifolds/images/examples/WebSampling/noisyGaussianPeak-Laplace_links.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "GraphWindowSpecR2{lBound=-2.666666666666667, rBound=2.666666666666666, bBound=-1.4999999999999996, tBound=2.499999999999999, xResolution=640, yResolution=480}"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "GraphWindowSpecR2{lBound=-2.666666666666666, rBound=2.6666666666666656, bBound=-1.4999999999999998, tBound=2.4999999999999987, xResolution=640, yResolution=480}"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "GraphWindowSpecR2{lBound=-2.666666666666666, rBound=2.6666666666666656, bBound=-1.5, tBound=2.4999999999999987, xResolution=640, yResolution=480}"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "GraphWindowSpecR2{lBound=-2.666666666666666, rBound=2.666666666666666, bBound=-1.5, tBound=2.499999999999999, xResolution=640, yResolution=480}"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plotWindow [plot . fmap colourscheme $ gaußianPeak grid δ ]\n",
    "plotWindow [plot . fmap colourscheme_nabla . differentiateUncertainWebFunction\n",
    "              $ gaußianPeak grid δ ]\n",
    "plotWindow [plot . fmap (colourscheme_nabla . projectShade (lensEmbedding (1*∂id/∂_x)))\n",
    "   . differentiateUncertainWebFunction . differentiateUncertainWebFunction\n",
    "       $ gaußianPeak grid δ ]\n",
    "plotWindow [plot . fmap (colourscheme_nabla . projectShade (lensEmbedding ((1*∂id/∂_x)/∂id)))\n",
    "   . differentiateUncertainWebFunction . differentiateUncertainWebFunction\n",
    "       $ gaußianPeak grid δ ]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Evaluating a function sampled on a 2D web](https://raw.githubusercontent.com/leftaroundabout/manifolds/master/manifolds/images/examples/WebSampling/noisyGaussianPeak_nodes.png)\n",
    "![Numerically differentiating a function on the nodes of a 2D web: first nabla, then directionally](https://raw.githubusercontent.com/leftaroundabout/manifolds/master/manifolds/images/examples/WebSampling/noisyGaussianPeak-∂x.∇_nodes.png)\n",
    "![Numerically differentiating a function on the nodes of a 2D web: first directionally, then nabla](https://raw.githubusercontent.com/leftaroundabout/manifolds/master/manifolds/images/examples/WebSampling/noisyGaussianPeak-∇.∂x_nodes.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": "8.2.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}