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