packages feed

goal-geometry (empty) → 0.1

raw patch · 14 files changed

+1375/−0 lines, 14 filesdep +basedep +goal-coredep +goal-geometrysetup-changed

Dependencies added: base, goal-core, goal-geometry, hmatrix, vector

Files

+ Goal/Geometry.hs view
@@ -0,0 +1,27 @@+module Goal.Geometry+    (+    -- * Re-Exports+      module Goal.Geometry.Set+    , module Goal.Geometry.Manifold+    , module Goal.Geometry.Linear+    , module Goal.Geometry.Map+    , module Goal.Geometry.Map.Multilinear+    , module Goal.Geometry.Differential+    , module Goal.Geometry.Differential.Convex+    , module Goal.Geometry.Plot+    ) where+++-- Imports --+++-- Re-exports --++import Goal.Geometry.Set+import Goal.Geometry.Manifold+import Goal.Geometry.Linear+import Goal.Geometry.Map+import Goal.Geometry.Map.Multilinear+import Goal.Geometry.Differential+import Goal.Geometry.Differential.Convex+import Goal.Geometry.Plot
+ Goal/Geometry/Differential.hs view
@@ -0,0 +1,236 @@+-- | This module provides tools for working with differential and Riemannian+-- geometry.+module Goal.Geometry.Differential (+    -- * Tangent Spaces+    -- ** Types+      Tangent (Tangent, removeTangent)+    , Bundle (Bundle, removeBundle)+    , Partials (Partials)+    , Differentials (Differentials)+    -- ** Functions+    , gradientStep+    , projectTangent+    , tangentToBundle+    , bundleToTangent+    -- * Riemannian Manifolds+    , Riemannian (metric, flat, sharp)+    -- ** Gradient Pursuit+    , gradientAscent+    , vanillaGradientAscent+    , gradientDescent+    , vanillaGradientDescent+    ) where+++--- Imports ---+++import Prelude hiding (map,minimum,maximum)++-- Package --++import Goal.Core++import Goal.Geometry.Set+import Goal.Geometry.Manifold+import Goal.Geometry.Linear+import Goal.Geometry.Map+import Goal.Geometry.Map.Multilinear++-- Qualified --++import qualified Data.Vector.Storable as C+import qualified Numeric.LinearAlgebra.HMatrix as H++--import Data.Vector.Storable.UnsafeSerialize+++--- Differentiable Manifolds ---+++-- | 'Tangent' spaces on 'Manifold's are the basis for differential geometry.+-- 'Tangent' spaces are defined at each point on a differentiable 'Manifold'.+newtype Tangent c m = Tangent { removeTangent :: c :#: m } deriving (Eq, Read, Show)++-- | A 'Tangent' 'Bundle' is the original 'Manifold' combined with all its+-- 'Tangent' spaces.+newtype Bundle c m = Bundle { removeBundle :: m } deriving (Eq, Read, Show)++-- | The 'Partials' coordinate system is defined as the partial derivatives of+-- the coordinate functions at a particular point.+data Partials = Partials deriving (Eq, Read, Show)++-- | The 'Differentials' coordinate system represents the set of linear+-- functionals on the 'Tangent' space.+data Differentials = Differentials deriving (Eq, Read, Show)++gradientStep :: Manifold m => Double -> Partials :#: Tangent c m -> c :#: m+-- | 'gradientStep' follows takes a gradient in a particular tangent space and+-- transforms the point underlying the given tangent space by shifting it+-- slightly in the direction of the gradient.+gradientStep eps f' =+    let (Tangent p) = manifold f'+        x' = coordinates $ eps .> f'+     in fromCoordinates (manifold p) (coordinates p + x')++projectTangent :: d :#: Tangent c m -> c :#: m+-- | Returns the underlying 'Point' from a 'Tangent' vector.+projectTangent = removeTangent . manifold++bundleToTangent :: Manifold m => c :#: Bundle d m -> c :#: Tangent d m+-- | Converts a 'Point' on a 'Tangent' 'Bundle' into a 'Tangent' vector.+bundleToTangent p =+    let (cs,dcs) = C.splitAt (div (dimension $ manifold p) 2) $ coordinates p+        (Bundle m) = manifold p+     in fromCoordinates (Tangent $ fromCoordinates m cs) dcs++tangentToBundle :: Manifold m => c :#: Tangent d m -> c :#: Bundle d m+-- | Converts  a 'Tangent' vector into a 'Point' on a 'Tangent' 'Bundle'.+tangentToBundle cm =+    let (Tangent dm) = manifold cm+        m = manifold dm+     in fromCoordinates (Bundle m) $ coordinates dm C.++ coordinates cm++replicatedTangents :: Manifold m => d :#: Tangent c (Replicated m) -> [d :#: Tangent c m]+-- | Converts a 'Tangent' vector on a 'Replicated' 'Manifold' into a list of+-- 'Tangent' vectors.+replicatedTangents dp =+    let (Tangent p) = manifold dp+        ts = mapReplicated Tangent p+        cs = listCoordinates dp+     in zipWith fromList ts $ breakEvery (dimension $ head ts) cs+++-- Gradient Pursuit --++gradientAscent :: (Riemannian c m, Manifold m)+    => Double -- ^ Step size+    -> (c :#: m -> Differentials :#: Tangent c m) -- ^ Gradient calculator+    -> (c :#: m) -- ^ The initial point+    -> [c :#: m] -- ^ The gradient ascent+gradientAscent eps f' = iterate (gradientStep eps . sharp . f')++vanillaGradientAscent :: Manifold m+    => Double -- ^ Step size+    -> (c :#: m -> Differentials :#: Tangent c m) -- ^ Gradient calculator+    -> (c :#: m) -- ^ The initial point+    -> [c :#: m] -- ^ The gradient ascent+vanillaGradientAscent eps f' = iterate (gradientStep eps . breakChart . f')++gradientDescent :: (Riemannian c m, Manifold m)+    => Double -- ^ Step size+    -> (c :#: m -> Differentials :#: Tangent c m) -- ^ Gradient calculator+    -> (c :#: m) -- ^ The initial point+    -> [c :#: m] -- ^ The gradient ascent+gradientDescent eps = gradientAscent (-eps)++vanillaGradientDescent :: Manifold m+    => Double -- ^ Step size+    -> (c :#: m -> Differentials :#: Tangent c m) -- ^ Gradient calculator+    -> (c :#: m) -- ^ The initial point+    -> [c :#: m] -- ^ The gradient ascent+vanillaGradientDescent eps = vanillaGradientAscent (-eps)+++--- Riemannian Manifolds ---+++-- | 'Riemannian' 'Manifold's are differentiable 'Manifold's where associated+-- with each point in the 'Manifold' is a 'Tangent' space with a smoothly+-- varying inner product. 'flat' and 'sharp' correspond to lowering and+-- raising the indices via the musical isomorphism determined by the metric+-- tensor.+--+-- A 'Riemannian' 'Manifold' should should satisfy the law+--+-- > flat $ sharp p = p+--+class Manifold m => Riemannian c m where+    metric :: c :#: m -> Function Partials Differentials :#: Tensor (Tangent c m) (Tangent c m)+    flat :: Partials :#: Tangent c m -> Differentials :#: Tangent c m+    flat p = matrixApply (metric $ projectTangent p) p+    sharp :: Differentials :#: Tangent c m -> Partials :#: Tangent c m+    sharp p = matrixApply (matrixInverse . metric $ projectTangent p) p+++--- Instances ---+++-- Replicated --++instance (Manifold m, Riemannian c m) => Riemannian c (Replicated m) where+    metric p =+        let mtxs = mapReplicated (toHMatrix . metric) p+         in fromHMatrix (Tensor (Tangent p) (Tangent p)) $ H.diagBlock mtxs+    flat dp =+        fromCoordinates (manifold dp) . C.concat $ coordinates . flat <$> replicatedTangents dp+    sharp dp =+        fromCoordinates (manifold dp) . C.concat $ coordinates . sharp <$> replicatedTangents dp++-- Euclidean --++instance Riemannian Cartesian Continuum where+    metric p = fromList (Tensor (Tangent p) (Tangent p)) [1]+    flat = breakChart+    sharp = breakChart++instance Riemannian Cartesian Euclidean where+    metric p = fromHMatrix (Tensor (Tangent p) (Tangent p)) . H.ident . dimension $ manifold p+    flat = breakChart+    sharp = breakChart++-- Trivial higher order spaces --++instance (Manifold m, Riemannian c m) => Riemannian Partials (Tangent c m) where+    metric dp =+        fromCoordinates (Tensor (Tangent dp) (Tangent dp)) . coordinates . metric $ projectTangent dp+    sharp ddp = fromCoordinates (manifold ddp) . coordinates+        . sharp . fromCoordinates (manifold $ projectTangent ddp) $ coordinates ddp+    flat pdd = fromCoordinates (manifold pdd) . coordinates+        . flat . fromCoordinates (manifold $ projectTangent pdd) $ coordinates pdd++-- Tangent Spaces --++instance Manifold m => Manifold (Tangent c m) where+    dimension (Tangent p) = dimension $ manifold p++instance Manifold m => Manifold (Bundle c m) where+    dimension (Bundle m) = 2 * dimension m++-- Tanget Space Coordinates --++instance Primal Partials where+    type Dual Partials = Differentials++instance Primal Differentials where+    type Dual Differentials = Partials+++--- Graveyard ---++++{-+--- Functions ---+++pushForward :: (Manifold m, Manifold n)+    => Function c d :#: Tensor n m+    -> c :#: m+    -> Function Partials Partials :#: Tensor (Tangent d n) (Tangent c m)+-- | 'pushForward' takes a 'Map' between 'Manifold's and turns it into a map+-- between the 'Tangent' spaces of the 'Manifold's. Although this ought to be a+-- class, right now it's simply the trivial 'pushForward' as applied to linear+-- maps.+pushForward pq q = fromCoordinates (Tensor (Tangent $ matrixApply pq q) (Tangent q)) $ coordinates pq++pushForward0 :: (Manifold m, Manifold n)+    => Function c d :#: Tensor n m+    -> c :#: m+    -> d :#: n+    -> Function Partials Partials :#: Tensor (Tangent d n) (Tangent c m)+-- | 'pushForward0' takes a 'Map' between 'Manifold's and turns it into a map+-- between the 'Tangent' spaces of the 'Manifold's. In this version we can+-- specify the target space more directly.+pushForward0 pq q p = fromCoordinates (Tensor (Tangent p) (Tangent q)) $ coordinates pq+-}
+ Goal/Geometry/Differential/Convex.hs view
@@ -0,0 +1,52 @@+-- | Tools are also provided for convex analysis, as the dual structures of+-- convex analysis are equivalent to Riemannian manifolds with certain+-- properties.+module Goal.Geometry.Differential.Convex where+++--- Imports ---+++-- Goal --++import Goal.Geometry.Set+import Goal.Geometry.Manifold+import Goal.Geometry.Linear+import Goal.Geometry.Differential++--- Dually Flat Manifolds ---++-- | Although convex analysis is usually developed seperately from differential+-- geometry, it arrises naturally out of the theory of dually flat 'Manifold's.+--+-- A 'Manifold' is 'Legendre' for a particular coordinated system if it is+-- associated with a particular convex function on points of the manifold known+-- as a 'potential'.+class (Primal c, Manifold m) => Legendre c m where+    potential :: (c :#: m) -> Double+    potentialDifferentials :: (c :#: m) -> Differentials :#: Tangent c m++potentialMapping :: Legendre c m => (c :#: m) -> Dual c :#: m+potentialMapping p = fromCoordinates (manifold p) . coordinates $ potentialDifferentials p++-- | Computes the 'divergence' between two points.+divergence :: (Primal c, Legendre c m, Legendre (Dual c) m) => (c :#: m) -> (Dual c :#: m) -> Double+divergence pp dq = potential pp + potential dq - (pp <.> dq)++legendreFlat :: (Legendre c m, Riemannian c m) => c :#: m -> c :#: m -> Dual c :#: m+-- | Applies 'flat' to the second input, based on the tangent space at the first input.+legendreFlat mp err = fromCoordinates (manifold mp) . coordinates . flat . fromCoordinates (Tangent mp) $ coordinates err+++--- Instances ---+++-- Generic --++-- Direct Sums --++instance Legendre c m => Legendre c (Replicated m) where+    potential ps = sum $ mapReplicated potential ps+    potentialDifferentials ps =+        let dps = mapReplicated potentialDifferentials ps+        in fromCoordinates (Tangent ps) . coordinates $ joinReplicated dps
+ Goal/Geometry/Linear.hs view
@@ -0,0 +1,108 @@+-- | The 'Linear' module provides the tools for treating a given 'Manifold' as a+-- linear space.+module Goal.Geometry.Linear (+    -- * Vector Spaces+      (<+>)+    , (.>)+    , (<->)+    , (/>)+    , meanPoint+    -- * Dual Spaces+    , Primal+    , Dual+    , (<.>)+    ) where++--- Imports ---++import Prelude hiding (map,minimum,maximum)++-- Package --++import Goal.Core hiding (dot)+import Goal.Geometry.Manifold++-- Unqualified --++import Numeric.LinearAlgebra.HMatrix hiding (Field,(><),(<>),(<.>))++--import Data.Vector.Storable.UnsafeSerialize+++--- Vector Spaces on Manifolds ---+++infixl 6 <+>+(<+>) :: Manifold m => c :#: m -> c :#: m -> c :#: m+-- | Vector addition of points on a manifold.+(<+>) p p' = fromCoordinates (manifold p) (coordinates p' + coordinates p)+{-+  | m == manifold p' = fromCoordinates m (coordinates p' + coordinates p)+  | otherwise = error "Attempting to add points from distinct manifolds."+    where m = manifold p+          -}++infixl 6 <->+(<->) :: Manifold m => c :#: m -> c :#: m -> c :#: m+-- | Vector subtraction of points on a manifold.+(<->) p p' = fromCoordinates (manifold p) (coordinates p - coordinates p')+{-+  | m == manifold p' = fromCoordinates m (coordinates p - coordinates p')+  | otherwise = error "Attempting to subtract points from distinct manifolds."+    where m = manifold p+          -}+++infix 7 .>+(.>) :: Manifold m => Double -> c :#: m -> c :#: m+-- | Scalar multiplication of points on a manifold.+(.>) a = alterCoordinates (*a)++infix 7 />+(/>) :: Manifold m => Double -> c :#: m -> c :#: m+-- | Scalar division of points on a manifold.+(/>) a v = recip a .> v+++--- Dual Spaces ---+++-- | 'Primal' charts have a 'Dual' coordinate system. The 'Dual' coordinate+-- system is the system which determines the dual basis of the dual vector+-- space via the restriction that the inner product '<.>' be the dot product.+--+-- Since finite dimensional vector spaces are isomorphic to their dual spaces+-- through the dual basis,  vector space duality is handled purely at the level+-- of coordinates in Goal -- that is, 'Primal' and 'Dual' coordinates are+-- considered different ways of describing the same fundamental objects. In+-- practice, encoding this relationship purely at the level of Charts saves a+-- great deal of computational effort.+class (Dual (Dual c)) ~ c => Primal c where+    type Dual c :: *++infix 7 <.>+(<.>) :: c :#: m -> Dual c :#: m -> Double+-- | '<.>' is the inner product between a dual pair of 'Point's. The defining+-- property of 'Dual' coordinate systems is that the inner product can be+-- expressed as a dot product.+(<.>) p q = dot (coordinates p) (coordinates q)++-- Utility --++meanPoint :: Manifold m => [c :#: m] -> c :#: m+-- | Finds the midpoint amongst a set of vectors in a convex set.+meanPoint ps = fromCoordinates (manifold $ head ps) . mean $ coordinates <$> ps+  {-+  | all (== m) (manifold <$> ps) = fromCoordinates m . mean $ coordinates <$> ps+  | otherwise = error "Attempting to add points from distinct manifolds."+    where m = manifold $ head ps+          -}+++--- Instances ---+++-- Cartesian Spaces --++instance Primal Cartesian where+    type Dual Cartesian = Cartesian
+ Goal/Geometry/Manifold.hs view
@@ -0,0 +1,284 @@+-- | This module provides the core mathematical definitions used by the rest of+-- Goal. In Goal, all mathematical structures are 'Manifold's, even when they are+-- not especially complicated ones; 'Manifold's may indicate highly articulated+-- structures, but may also indicate simpler concepts such as (vector) spaces.+--+-- 'Manifold's are sets of points which can be described locally as 'Euclidean'+-- spaces. In geometry, a point is typically a member of the actual 'Manifold'.+-- However, arbitrary types of points will often be difficult to represent+-- directly, and so points in Goal are always represented in terms of their+-- 'Coordinates' in terms of a given chart.+--+-- Charts are in turn represented by phantom types. Mathematically, charts are+-- maps between the 'Manifold' and the relevant 'Cartesian' coordinate system.+-- However, since we do not represent the points of a 'Manifold' explicility,+-- we also cannot represent Charts explicitly. As such, Atlases merely index a+-- point so as to indicate how to interpret its particular 'Coordinates'.+module Goal.Geometry.Manifold+    ( -- * Manifolds+      Manifold (dimension)+    , Transition (transition)+    -- ** Sets+    , Embedded (Embedded, disembed)+    -- ** Points+    , Coordinates+    , (:#:) (coordinates, manifold)+    , coordinate+    , chart+    , breakChart+    , alterChart+    , listCoordinates+    , alterCoordinates+    , toPair+    -- ** Charts+    , Cartesian (Cartesian)+    , Polar (Polar)+    -- ** Constructors+    , fromList+    , fromCoordinates+    , euclideanPoint+    , realNumber+    -- * Direct Sums+    -- ** Replicated+    , mapReplicated+    , joinReplicated+    , concatReplicated+    -- ** DirectSum+    , joinPair+    , splitPair+    , joinPair'+    , splitPair'+    , joinTriple+    , splitTriple+    , joinTriple'+    , splitTriple'+    ) where+++--- Imports ---+++-- Goal --++import Goal.Core++import Goal.Geometry.Set++-- Qualified --++import qualified Data.Vector.Storable as C+++--- Manifolds ---+++-- | A geometric object with a certain 'dimension'. We assume that a 'Manifold'+-- somehow represents all the geometric, coordinate independent structure under+-- consideration. 'Manifold's should satisfy+--+-- > dimension m = length $ coordinates (Point m cs)+--+class Eq m => Manifold m where+    dimension :: m -> Int++-- | A point is an element of a 'Manifold' 'm' in terms of a particular+-- chart 'c'.+data c :#: m = Point+    { coordinates :: !Coordinates+    , manifold :: m } deriving (Eq, Read, Show)++infixr 1 :#:++coordinate :: Int -> c :#: m -> Double+coordinate n (Point cs _) = cs C.! n++data Embedded m c = Embedded { disembed :: m } deriving (Eq, Read, Show)++chart :: Manifold m => c -> c :#: m -> c :#: m+-- | 'chart' allows one to specify the Atlas of a new point. This is often+-- necessary when typeclass methods are used to generate points under a+-- variety of coordinate systems.+chart _ = id++breakChart :: Manifold m => c :#: m -> d :#: m+breakChart p = Point (coordinates p) (manifold p)++alterChart :: Manifold m => d -> c :#: m -> d :#: m+-- | Combines 'breakChart' and 'chart'.+alterChart _ = breakChart++toPair :: c :#: m -> (Double,Double)+toPair p = (coordinate 0 p,coordinate 1 p)++alterCoordinates :: Manifold m => (Double -> Double) -> c :#: m -> c :#: m+-- | 'alterCoordinates' allows one to map a function over the 'coordinates' of a+-- point without changing the chart.+alterCoordinates f (Point cs m) = Point (C.map f cs) m++listCoordinates :: c :#: m -> [Double]+-- | Returns the 'Coordinates' of the point in list form.+listCoordinates (Point cs _) = C.toList cs++-- | A 'transition' involves taking a point represented by the chart 'c',+-- and re-representing in terms of the chart 'd'. This will usually require+-- recomputation of the 'Coordinates'. 'Transition's should satisfy the law+--+-- > transition $ transition p = p+--+class Transition c d m where+    transition :: c :#: m -> d :#: m++fromList :: Manifold m => m -> [Double] -> c :#: m+-- | 'fromList' builds points without the need to work with vectors.+fromList m cs = fromCoordinates m $ C.fromList cs++fromCoordinates :: Manifold m => m -> Coordinates -> c :#: m+fromCoordinates m cs -- = Point cs m+    | dimension m == C.length cs = Point cs m+    | otherwise = error+        $ "Coordinate dimension (" ++ show (C.length cs) ++ ") does not match Manifold dimension (" ++ show (dimension m) ++ ")."++euclideanPoint :: [Double] -> Cartesian :#: Euclidean+-- | A convenience function for building 'Euclidean' vectors.+euclideanPoint xs = fromList (Euclidean $ length xs) xs++realNumber :: Double -> Cartesian :#: Continuum+-- | A convenience function for building elements of a 'Continuum'.+realNumber x = fromList Continuum [x]++--- Construction ---+++-- Euclidean --++-- | The 'Cartesian' coordinate system.+data Cartesian = Cartesian++-- | The 'Polar' coordinate system.+data Polar = Polar++-- | A function to map functions over a point on a 'Replicated' 'Manifold'.+mapReplicated :: Manifold m => (c :#: m -> x) -> c :#: Replicated m -> [x]+mapReplicated pf ps =+    let (Replicated m k) = manifold ps+        cs = coordinates ps+        b = dimension m+     in [ pf . fromCoordinates m $ C.slice (i * b) b cs | i <- [0.. k -1 ] ]++joinReplicated :: Manifold m => [c :#: m] -> c :#: Replicated m+-- | Joins a list of distributions into a 'Replicated' 'Manifold'. Be advised that this function assumes+-- that the families of the individual distributions are equal.+joinReplicated ps =+    Point (foldl1' (C.++) (coordinates <$> ps)) $ Replicated (manifold $ head ps) (length ps)++concatReplicated :: c :#: Replicated m -> c :#: Replicated m -> c :#: Replicated m+-- | Joins two 'Replicated' 'Manifold's.+concatReplicated (Point cs (Replicated m x)) (Point cs' (Replicated _ y)) = Point (cs C.++ cs') $ Replicated m (x + y)++-- Direct Sums --++joinPair :: (Manifold m, Manifold n) => c :#: m -> d :#: n -> (c,d) :#: (m,n)+-- | Joins a pair of Points into a Point on the the direct sum of the underlying Charts and 'Manifold's.+joinPair = unsafeJoinPair++splitPair :: (Manifold m, Manifold n) => (c,d) :#: (m,n) -> (c :#: m, d :#: n)+-- | Splits a direct sum pair.+splitPair = unsafeSplitPair++joinPair' :: (Manifold m, Manifold n) => c :#: m -> c :#: n -> c :#: (m,n)+-- | Alternative version where we assume that the Charts are shared.+joinPair' = unsafeJoinPair++splitPair' :: (Manifold m, Manifold n) => c :#: (m,n) -> (c :#: m, c :#: n)+-- | Alternative version where we assume that the Charts are shared.+splitPair' = unsafeSplitPair++unsafeJoinPair :: (Manifold m, Manifold n) => c :#: m -> d :#: n -> e :#: (m,n)+unsafeJoinPair cm dn =+    fromCoordinates (manifold cm,manifold dn) $ coordinates cm C.++ coordinates dn++unsafeSplitPair :: (Manifold m, Manifold n) => c :#: (m,n) -> (d :#: m, e :#: n)+unsafeSplitPair cmn =+    let (m,n) = manifold cmn+        cs = coordinates cmn+        (mcs,ncs) = C.splitAt (dimension m) cs+     in (fromCoordinates m mcs, fromCoordinates n ncs)++joinTriple :: (Manifold m, Manifold n, Manifold o) => c :#: m -> d :#: n -> e :#: o -> (c,d,e) :#: (m,n,o)+-- | Joins a triple of Points into a Point on the the direct sum of the underlying Charts and 'Manifold's.+joinTriple = unsafeJoinTriple++splitTriple :: (Manifold m, Manifold n, Manifold o) => (c,d,e) :#: (m,n,o) -> (c :#: m, d :#: n, e :#: o)+-- | Splits a direct sum triple.+splitTriple = unsafeSplitTriple++joinTriple' :: (Manifold m, Manifold n, Manifold o) => c :#: m -> c :#: n -> c :#: o -> c :#: (m,n,o)+-- | Alternative version where we assume that the Charts are shared.+joinTriple' = unsafeJoinTriple++splitTriple' :: (Manifold m, Manifold n, Manifold o) => c :#: (m,n,o) -> (c :#: m, c :#: n, c :#: o)+-- | Alternative version where we assume that the Charts are shared.+splitTriple' = unsafeSplitTriple++unsafeJoinTriple :: (Manifold m, Manifold n, Manifold o) => c :#: m -> d :#: n -> e :#: o -> f :#: (m,n,o)+unsafeJoinTriple cm dn eo =+    fromCoordinates (manifold cm, manifold dn, manifold eo) $ coordinates cm C.++ coordinates dn C.++ coordinates eo++unsafeSplitTriple :: (Manifold m, Manifold n, Manifold o) => c :#: (m,n,o) -> (d :#: m, e :#: n, f :#: o)+unsafeSplitTriple cmno =+    let (m,n,o) = manifold cmno+        (mcs,cs') = C.splitAt (dimension m) $ coordinates cmno+        (ncs,ocs) = C.splitAt (dimension n) cs'+     in (fromCoordinates m mcs, fromCoordinates n ncs, fromCoordinates o ocs)+++--- Instances ---+++instance Transition c c m where+    transition = id++-- Embedded --++instance Manifold m => Set (Embedded m c) where+    type Element (Embedded m c) = c :#: m++-- Euclidean --++instance Manifold Euclidean where+    dimension (Euclidean n) = n++instance Manifold Continuum where+    dimension _ = 1++instance Transition Polar Cartesian Euclidean where+    transition p =+        let r:phis = listCoordinates p+            phiss = reverse . tails $ reverse phis+            m = manifold p+            xs = [ r * cos phi * product (sin <$> phis') | (phi,phis') <- zip phis phiss ]+         in fromList m $ xs ++ [r * product (sin <$> phis)]++instance Transition Cartesian Polar Euclidean where+    transition p =+        let (Euclidean n) = manifold p+            xs = listCoordinates p+            xs2 = listCoordinates $ alterCoordinates (^2) p+            r = sqrt $ sum xs2+            (phis,phin0:_) = splitAt (n-2) [ acos $ xi / sqrt (sum xs2i) | (xi,xs2i) <- zip xs (tails xs2) ]+            xn = last xs+            phin = if xn > 0 then phin0 else 2*pi - phin0+         in fromList (Euclidean n) $ r : (phis ++ [phin])++-- DirectSum --++instance (Manifold m, Manifold n) => Manifold (m,n) where+    dimension (m,n) = dimension m + dimension n++instance (Manifold m, Manifold n, Manifold o) => Manifold (m,n,o) where+    dimension (m,n,o) = dimension m + dimension n + dimension o++-- Replicated --++instance Manifold m => Manifold (Replicated m) where+    dimension (Replicated m rn) = dimension m * rn
+ Goal/Geometry/Map.hs view
@@ -0,0 +1,76 @@+-- | The Map module provides tools for developing function space 'Manifold's.+-- A map is a 'Manifold' where the 'Point's of the Manifold represent+-- parametric functions between 'Manifold's. The defining feature of 'Map's is+-- that they have a particular 'Domain' and 'Codomain', which themselves are+-- 'Manifold's.++module Goal.Geometry.Map (+    -- * Maps+      Map (Domain, domain, Codomain, codomain)+    , Apply ((>.>), (>$>))+    -- * Map Charts+    , Function (Function)+    ) where+++--- Imports ---+++-- Goal --++import Goal.Geometry.Manifold++--- Maps between Manifolds ---++-- Charts on Maps --++data Function c d = Function c d+-- | 'Function' Charts help track Charts on the 'Domain' and 'Codomain'. The+-- first Chart corresponds to the 'Domain's chart.++class Manifold m => Map m where+    type Domain m :: *+    domain :: m -> Domain m+    type Codomain m :: *+    codomain :: m -> Codomain m++class Map m => Apply c d m where+    -- | 'Map' application.+    (>.>) :: Function c d :#: m -> c :#: Domain m -> d :#: Codomain m+    (>.>) f x = head $ f >$> [x]+    -- | 'Map' list application. May sometimes have a more efficient implementation+    -- than simply list-mapping (>.>).+    (>$>) :: Function c d :#: m -> [c :#: Domain m] -> [d :#: Codomain m]+    (>$>) f = map (f >.>)++infix 8 >.>+infix 8 >$>++++{-+--- Tables ---+++newtype Table s = Table s deriving (Eq, Read, Show)+++--- Instances ---+++-- Table --++instance Discrete s => Manifold (Table s) where+    dimension (Table s) = length $ elements s++instance Discrete s => Function Cartesian (Table s) where+    type Domain Cartesian (Table s) = s+    domain cm = let (Table s) = manifold cm in s+    type Codomain Cartesian (Table s) = Continuum+    codomain _ = Continuum+    (>.>) cm k =+        let ctgs = listCoordinates cm+            Just (ctg,_) = find ((==k) . snd) . zip ctgs . elements $ domain cm+         in ctg+    (>$>) cm ks = (cm >.>) <$> ks+-}
+ Goal/Geometry/Map/Multilinear.hs view
@@ -0,0 +1,211 @@+-- | The Map module provides tools for developing function space 'Manifold's.+-- A map is a 'Manifold' where the 'Point's of the Manifold represent+-- parametric functions between 'Manifold's. The defining feature of 'Map's is+-- that they have a particular 'Domain' and 'Codomain', which themselves are+-- 'Manifold's.++module Goal.Geometry.Map.Multilinear (+    -- * Tensors+      Tensor (Tensor)+    -- ** Construction+    , (>.<)+    -- ** Matrix Operations+    , (<#>)+    , matrixRank+    , matrixInverse+    , matrixTranspose+    , matrixSquareRoot+    , matrixApply+    , matrixMap+    , matrixDiagonalConcatenate+    -- ** Cartesian+    , coordinateTransform+    , linearProjection+    -- ** HMatrix Conversion+    , toHMatrix+    , fromHMatrix+    -- * Affine Functions+    , Affine (Affine)+    , splitAffine+    , joinAffine+    ) where++--- Imports ---++import Prelude hiding (map,minimum,maximum)++-- Package --++import Goal.Core++import Goal.Geometry.Set+import Goal.Geometry.Manifold+import Goal.Geometry.Linear+import Goal.Geometry.Map++-- Qualified --++import qualified Data.Vector.Storable as C+import qualified Numeric.LinearAlgebra.HMatrix as H++--import Data.Vector.Storable.UnsafeSerialize++++--- Affine Functions ---+++-- | 'Manifold's of 'Affine' functions.+data Affine m n = Affine m n deriving (Eq, Read, Show)++splitAffine :: (Manifold m, Manifold n) => Function c d :#: Affine m n -> (d :#: m, Function c d :#: Tensor m n)+-- | Splits an 'Point' on an 'Affine' space into a matrix and a constant.+splitAffine aff =+    let (Affine m n) = manifold aff+        tns = Tensor m n+        css = coordinates aff+        (mcs,mtxcs) = C.splitAt (dimension m) css+     in (fromCoordinates m mcs, fromCoordinates tns mtxcs)++joinAffine :: (Manifold m, Manifold n) => d :#: m -> Function c d :#: Tensor m n -> Function c d :#: Affine m n+-- | Combines a matrix and a constant into 'Point' on an 'Affine' space.+joinAffine dm mtx =+    let (Tensor m n) = manifold mtx+     in fromCoordinates (Affine m n) $ coordinates dm C.++ coordinates mtx++-- Tensor Products --++-- | 'Manifold' of 'Tensor's given by the tensor product of the underlying pair of 'Manifold's.+data Tensor m n = Tensor m n deriving (Eq, Read, Show)++toHMatrix :: Manifold n => c :#: Tensor m n -> H.Matrix Double+-- | Converts a point on a 'Tensor' product manifold to a matrix for snappy+-- calculation.+toHMatrix pq =+    let (Tensor _ m) = manifold pq+     in H.reshape (dimension m) $ coordinates pq++fromHMatrix :: (Manifold m, Manifold n) => Tensor m n -> H.Matrix Double -> c :#: Tensor m n+fromHMatrix tns = fromCoordinates tns . H.flatten++matrixRank :: (Manifold m, Manifold n) => c :#: Tensor m n -> Int+matrixRank = H.rank . toHMatrix++(>.<) :: (Manifold m, Manifold n) => d :#: m -> c :#: n -> Function (Dual c) d :#: Tensor m n+-- | '>.<' denotes the outer product between two points. It provides a way of+-- constructing matrices of the 'Tensor' product space.+(>.<) p q = fromHMatrix (Tensor (manifold p) $ manifold q) $ coordinates p `H.outer` coordinates q++(<#>) :: (Manifold m, Manifold n, Manifold o)+      => Function d e :#: Tensor m n -> Function c d :#: Tensor n o -> Function c e :#: Tensor m o+-- | Tensor product composition.+(<#>) p q =+    let (Tensor m _) = manifold p+        (Tensor _ o) = manifold q+     in fromHMatrix (Tensor m o) $ toHMatrix p <> toHMatrix q++matrixSquareRoot :: Manifold m => c :#: Tensor m m -> c :#: Tensor m m+-- | The square root of a matrix.+matrixSquareRoot pq = fromHMatrix (manifold pq) . H.sqrtm $ toHMatrix pq++matrixInverse :: (Manifold n, Manifold m) => Function c d :#: Tensor m n -> Function d c :#: Tensor n m+-- | The inverse of a given 'Tensor' point.+matrixInverse pq =+    let Tensor m n = manifold pq+     in fromHMatrix (Tensor n m) . H.inv $ toHMatrix pq++matrixTranspose :: (Manifold m, Manifold n) => Function c d :#: Tensor m n -> Function (Dual d) (Dual c) :#: Tensor n m+-- | The transpose of a given 'Tensor' point.+matrixTranspose pq =+    let Tensor m n = manifold pq+     in fromHMatrix (Tensor n m) . H.tr $ toHMatrix pq++matrixDiagonalConcatenate :: (Manifold m, Manifold n, Manifold o, Manifold p)+    => Function c d :#: Tensor m n+    -> Function e f :#: Tensor o p+    -> Function (c,e) (d,f) :#: Tensor (m,o) (n,p)+-- | Creates a block diagonal matrix.+matrixDiagonalConcatenate cdmn efop =+    let (Tensor m n) = manifold cdmn+        (Tensor o p) = manifold efop+     in fromHMatrix (Tensor (m,o) (n,p)) $ H.diagBlock [toHMatrix cdmn, toHMatrix efop]+++coordinateTransform :: Manifold m => [c :#: m] -> Function Cartesian c :#: Tensor m Euclidean+-- | Returns the coordinate transformation from 'Euclidean' space into the space+-- defined by the given basis vectors. This is a glorified fromColumns function.+coordinateTransform bss =+    fromHMatrix (Tensor (manifold $ head bss) . Euclidean $ length bss) . H.fromColumns $ coordinates <$> bss++linearProjection :: Manifold m => [Cartesian :#: m] -> Function Cartesian Cartesian :#: Tensor m m+-- | Returns the linear projection operator for the given subset of basis vectors.+linearProjection bss =+    let mtx = coordinateTransform bss+        mtxt = matrixTranspose mtx+     in mtx <#> matrixInverse (mtxt <#> mtx) <#> mtxt++matrixApply :: (Manifold m, Manifold n) => (Function c d :#: Tensor n m) -> (c :#: m) -> d :#: n+-- | Matrix vector multiplication.+matrixApply pq p =+    let (Tensor n _) = manifold pq+     in fromCoordinates n $ toHMatrix pq H.#> coordinates p+    {-+    let (Tensor n m) = manifold pq+     in if m == manifold p+          then fromCoordinates n $ toHMatrix pq H.#> coordinates p+          else error "matrix applied to wrong Manifold"+          -}++matrixMap :: (Manifold m, Manifold n) => (Function c d :#: Tensor m n) -> [c :#: n] -> [d :#: m]+-- | Mapped matrix vector multiplication, where we first turn the input vectors into a matrix itself (this can greatly improve computation time).+matrixMap pq ps =+    let (Tensor n _) = manifold pq+        mtx = toHMatrix pq+        xs = H.fromColumns $ coordinates <$> ps+     in map (fromCoordinates n) . H.toColumns $ mtx <> xs+    {-+    let (Tensor n m) = manifold pq+        mtx = toHMatrix pq+        xs = H.fromColumns $ coordinates <$> ps+     in if all (== m) $ manifold <$> ps+           then map (fromCoordinates n) . H.toColumns $ mtx <> xs+           else error "matrix applied to wrong Manifold"+           -}+++--- Instances ---+++-- Tensor Products --++instance (Manifold m, Manifold n) => Manifold (Tensor n m) where+    dimension (Tensor n m) = dimension m * dimension n++instance (Manifold m, Manifold n) => Map (Tensor m n) where+    type Domain (Tensor m n) = n+    domain (Tensor _ n) = n+    type Codomain (Tensor m n) = m+    codomain (Tensor m _) = m++instance (Manifold m, Manifold n) => Apply c d (Tensor m n) where+    (>.>) = matrixApply+    (>$>) = matrixMap++-- Affine Map --++instance (Manifold m, Manifold n) => Manifold (Affine m n) where+    dimension (Affine m n) = dimension m * dimension n + dimension m++instance (Manifold m, Manifold n) => Map (Affine m n) where+    type Domain (Affine m n) = n+    domain (Affine _ n) = n+    type Codomain (Affine m n) = m+    codomain (Affine m _) = m++instance (Manifold m, Manifold n) => Apply c d (Affine m n) where+    (>.>) p x =+        let (b,mtx) = splitAffine p+         in mtx >.> x <+> b+    (>$>) p xs =+        let (b,mtx) = splitAffine p+         in map (<+> b) $ mtx >$> xs
+ Goal/Geometry/Plot.hs view
@@ -0,0 +1,29 @@+module Goal.Geometry.Plot where+++--- Imports ---+++-- Goal --++import Goal.Core+import Goal.Geometry.Set++import qualified Data.Vector.Storable as C++coordinateLogHistogram :: Int -> String -> [String] -> [Coordinates] -> Layout Double Double+coordinateLogHistogram nbns ttl ttls css =+    let bplt = plot_bars_titles .~ ttls $ logHistogramPlot0 nbns (C.toList <$> css) def+     in layout_title .~ ttl+        $ layout_y_axis . laxis_override .~ axisGridHide+        $ layout_x_axis . laxis_override .~ axisGridHide+        $ logHistogramLayout bplt def++coordinateHistogram :: Int -> String -> [String] -> [Coordinates] -> Layout Double Double+coordinateHistogram nbns ttl ttls css =+    let bplt = plot_bars_titles .~ ttls+            $ histogramPlot0 nbns (C.toList <$> css) def+     in layout_title .~ ttl+        $ layout_y_axis . laxis_override .~ axisGridHide+        $ layout_x_axis . laxis_override .~ axisGridHide+        $ histogramLayout bplt def
+ Goal/Geometry/Set.hs view
@@ -0,0 +1,134 @@+-- | A module for describing 'Set's of 'Element's. Necessary in a few cases (such as discrete sets) that 'Manifold's don't handle well.+module Goal.Geometry.Set+    ( -- * Sets+      Set+    , Element+    , Discrete (elements)+    -- * Instances+    -- ** Discrete+    , Boolean (Boolean)+    , NaturalNumbers (NaturalNumbers)+    , Integers (Integers)+    -- ** Continuous+    , Coordinates+    , Euclidean (Euclidean)+    , Continuum (Continuum)+    -- * Combinators+    -- ** Replicated+    , Replicated (Replicated)+    ) where+++--- Imports ---+++-- Goal --++import Goal.Core++-- Qualified --++import qualified Data.Vector.Storable as C+++--- Classes ---+++-- | 'Set's are collections of distinguishable 'Element's.+class (Eq s, Eq (Element s)) => Set s where+    type Element s :: *+++-- | A 'Discrete' 'Set' is one where we can list its elements. The+-- returned list should satisfy the law+--+-- > elements s = nub $ elements s+--+class Set s => Discrete s where+    elements :: s -> [Element s]+++--- Types ---+++-- Discrete --++-- | The set of natural numbers.+data NaturalNumbers = NaturalNumbers deriving (Eq,Read,Show)++-- | The set of integers.+data Integers = Integers deriving (Eq,Read,Show)++-- | 'True' and 'False'.+data Boolean = Boolean deriving (Eq,Read,Show)++-- Continuous  --++-- | 'Euclidean' space.+newtype Euclidean = Euclidean Int deriving (Eq,Read,Show)++-- | One dimensional 'Euclidean' space.+data Continuum = Continuum deriving (Eq,Read,Show)++-- | 'Element's of 'Euclidean' spaces are referred to as 'Coordinates'.+type Coordinates = C.Vector Double++-- Replicated --++-- | A 'Replicated' set is a single set multiplied a specified number of times+-- via the Cartesian product.+data Replicated m = Replicated !m !Int deriving (Eq,Read,Show)+++--- Instances ---+++-- Discrete --++instance Set NaturalNumbers where+    type Element NaturalNumbers = Int++instance Discrete NaturalNumbers where+    elements _ = [0..]++instance Set Integers where+    type Element Integers = Int++instance Discrete Integers where+    elements _ = (0:) $ concat [ [-k,k] | k <- [1..] ]++instance Set Boolean where+    type Element Boolean = Bool++instance Discrete Boolean where+    elements _ = [True,False]++instance Eq k => Set [k] where+    type Element [k] = k++instance Eq k => Discrete [k] where+    elements = id++-- Continuous --++instance Set Continuum where+    type Element Continuum = Double++instance Set Euclidean where+    type Element Euclidean = Coordinates+++-- Replicated --++instance Set s => Set (Replicated s) where+    type Element (Replicated s) = [Element s]+++instance Discrete s => Discrete (Replicated s) where+    elements (Replicated s n) = replicateM n $ elements s++-- Direct Sums --++instance (Set s, Set r) => Set (s,r) where+    type Element (s,r) = (Element s,Element r)+
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright (c) 2014, Sacha Sokoloski++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++    * Redistributions of source code must retain the above copyright+      notice, this list of conditions and the following disclaimer.++    * Redistributions in binary form must reproduce the above+      copyright notice, this list of conditions and the following+      disclaimer in the documentation and/or other materials provided+      with the distribution.++    * Neither the name of Sacha Sokoloski nor the names of other+      contributors may be used to endorse or promote products derived+      from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ goal-geometry.cabal view
@@ -0,0 +1,59 @@+name: goal-geometry+version: 0.1+synopsis: Scientific computing on geometric objects+description: This library provides all the types and classes essential for+    defining manifolds. This includes definitions and algorithms for sets,+    points, linear and multilinear algebra, function spaces, basic differential+    geometry, and convex analysis.+license: BSD3+license-file: LICENSE+author: Sacha Sokoloski+maintainer: sokolo@mis.mpg.de+category: Math+build-type: Simple+cabal-version: >=1.10++library+    exposed-modules:+        Goal.Geometry,+        Goal.Geometry.Set,+        Goal.Geometry.Manifold,+        Goal.Geometry.Linear,+        Goal.Geometry.Map,+        Goal.Geometry.Map.Multilinear,+        Goal.Geometry.Differential,+        Goal.Geometry.Differential.Convex,+        Goal.Geometry.Plot+    build-depends:+        base==4.*,+        goal-core==0.1,+        vector==0.11.*,+        hmatrix==0.17.*+    default-extensions: TypeOperators, TypeFamilies, MultiParamTypeClasses,+        FlexibleInstances, FlexibleContexts+    default-language: Haskell2010+    ghc-options: -O2 -Wall -fno-warn-type-defaults -fno-warn-missing-signatures++executable coordinates+    main-is: coordinates.hs+    hs-source-dirs: scripts+    ghc-options: -Wall -O2 -threaded -rtsopts -fno-warn-type-defaults+        -fno-warn-missing-signatures -fno-warn-unused-do-bind+    build-depends:+        base==4.*,+        goal-core==0.1,+        goal-geometry==0.1+    default-language: Haskell2010++executable gradient-descent+    main-is: gradient-descent.hs+    hs-source-dirs: scripts+    ghc-options: -Wall -O2 -threaded -rtsopts -fno-warn-type-defaults+        -fno-warn-missing-signatures -fno-warn-unused-do-bind+    build-depends:+        base==4.*,+        goal-core==0.1,+        goal-geometry==0.1+    default-language: Haskell2010++
+ scripts/coordinates.hs view
@@ -0,0 +1,52 @@+--- Imports ---+++import Goal.Core+import Goal.Geometry+++--- Program ---+++-- Globals --++mxx = pi+mnx = -pi+mxy = pi+mny = -pi+nstps = 10+npnts = 50++axprms = LinearAxisParams (show . round) 5 5++hlns = [ [ [x,y] | x <- range mnx mxx npnts ] | y <- range mnx mxx nstps ]+vlns = [ [ [x,y] | y <- range mny mxy npnts ] | x <- range mny mxy nstps ]+lns0 = hlns ++ vlns++eclds = map euclideanPoint <$> lns0+plrs = map (chart Cartesian . transition . chart Polar . fromList (Euclidean 2)) <$> lns0++layoutMaker lns = execEC $ do++    layout_x_axis . laxis_override .= axisGridHide+    layout_x_axis . laxis_generate .= scaledAxis axprms (-2,2)+    layout_x_axis . laxis_title .= "x"++    layout_y_axis . laxis_override .= axisGridHide+    layout_y_axis . laxis_generate .= scaledAxis axprms (-2,2)+    layout_y_axis . laxis_title .= "y"++    plot . liftEC $ do++        plot_lines_values .= (map toPair <$> lns)+        plot_lines_style .= solidLine 3 (opaque black)+++-- Main --++main = do+    let lyt1 = layoutMaker eclds+        lyt2 = layoutMaker plrs+        rnbl = toRenderable $ StackedLayouts [StackedLayout lyt1, StackedLayout lyt2] False+    --renderableToAspectWindow False 400 800 rnbl+    void $ renderableToFile (FileOptions (200,400) PDF) "coordinates.pdf" rnbl
+ scripts/gradient-descent.hs view
@@ -0,0 +1,75 @@+--- Imports ---+++import Goal.Core+import Goal.Geometry+++--- Globals ---+++-- Functions --++f p = let (x,y) = toPair p in x^2 + 2*y^2 + (x-y)^2++df p =+    let (x,y) = toPair p+        x' = 2*x + 2*(x-y)+        y' = 4*y - 2*(x-y)+     in fromList (Tangent p) [x',y']++-- Plot --++res = 400+mn = -4+mx = 4+niso = 10+cntrf x y = f $ euclideanPoint [x,y]+rng = (mn,mx,res)+clrs = rgbaGradient (0.9,0,0,1) (0,0,0,1) niso+axprms = LinearAxisParams (show . round) 5 5++-- Gradient Descent --++p0 = euclideanPoint [-4,2]+eps = 0.01+nstps = 200+grds = take nstps $ gradientDescent eps df p0+++--- Main ---+++main = do++    -- Contour plots+    let rnbl = toRenderable . execEC $ do++            let cntrs = contours rng rng niso cntrf++            sequence_ $ do++                ((_,cntr),clr) <- zip cntrs clrs++                return . plot . liftEC $ do++                    plot_lines_style .= solidLine 3 clr+                    plot_lines_values .= cntr++            layout_x_axis . laxis_generate .= scaledAxis axprms (mn,mx)+            layout_x_axis . laxis_override .= axisGridHide+            layout_x_axis . laxis_title .= "x"+            layout_y_axis . laxis_generate .= scaledAxis axprms (mn,mx)+            layout_y_axis . laxis_override .= axisGridHide+            layout_y_axis . laxis_title .= "y"++            plot . liftEC $ do+                plot_points_style .= filledCircles 5 (opaque red)+                plot_points_values .= [(0,0)]++            plot . liftEC $ do+                plot_lines_style .= solidLine 3 (opaque black)+                plot_lines_values .= [toPair <$> grds]++    --void $ renderableToAspectWindow False 800 800 rnbl+    void $ renderableToFile (FileOptions (200,200) PDF) "gradient-descent.pdf" rnbl