goal-core 0.1 → 0.20

raw patch · 17 files changed

+2374/−628 lines, 17 filesdep +asyncdep +bytestringdep +cassavadep −Chartdep −Chart-cairodep −Chart-gtkdep ~basedep ~containersdep ~goal-corePVP ok

version bump matches the API change (PVP)

Dependencies added: async, bytestring, cassava, criterion, deepseq, directory, finite-typelits, ghc-typelits-knownnat, ghc-typelits-natnormalise, hmatrix, hmatrix-gsl, math-functions, mwc-probability, mwc-random, optparse-applicative, primitive, process, vector, vector-sized

Dependencies removed: Chart, Chart-cairo, Chart-gtk, cairo, colour, data-default-class, gtk, lens

Dependency ranges changed: base, containers, goal-core

API changes (from Hackage documentation)

Files

• Goal/Core.hs +57/−145
• Goal/Core/Circuit.hs +215/−0
• Goal/Core/Plot.hs +0/−264
• Goal/Core/Plot/Contour.hs +0/−163
• Goal/Core/Project.hs +185/−0
• Goal/Core/Util.hs +265/−0
• Goal/Core/Vector/Boxed.hs +245/−0
• Goal/Core/Vector/Generic.hs +233/−0
• Goal/Core/Vector/Generic/Internal.hs +11/−0
• Goal/Core/Vector/Generic/Mutable.hs +11/−0
• Goal/Core/Vector/Storable.hs +741/−0
• README.md +12/−0
• benchmarks/convolutions.hs +110/−0
• benchmarks/multiplications.hs +129/−0
• benchmarks/outer-products.hs +71/−0
• goal-core.cabal +89/−29
• scripts/contours.hs +0/−27

@@ -1,39 +1,43 @@+-- | This module re-exports a number of (compatible) modules from across base+-- and other libraries, as well as most of the modules in @goal-core@. It does+-- not re-export the Vector modules, which should be imported with+-- qualification. module Goal.Core     ( -- * Module Exports-      module Goal.Core.Plot+      module Goal.Core.Util+    , module Goal.Core.Project+    , module Goal.Core.Circuit     , module Data.Function+    , module Data.Functor+    , module Data.Foldable+    , module Data.Traversable     , module Data.Ord-    , module Data.Monoid-    , module Data.List     , module Data.Maybe     , module Data.Either-    , module Data.Default.Class+    , module Data.Finite+    , module Data.Csv+    , module Data.Proxy+    , module Data.Kind+    , module Data.Functor.Identity+    , module Data.Type.Equality     , module Control.Applicative     , module Control.Monad+    , module Control.Monad.Primitive     , module Control.Monad.ST     , module Control.Arrow-    , module Control.Lens.Type-    , module Control.Lens.Getter-    , module Control.Lens.Setter-    , module Control.Lens.TH     , module Control.Concurrent+    , module Control.DeepSeq     , module Numeric+    , module Numeric.SpecFunctions+    , module Options.Applicative+    , module GHC.TypeNats+    , module GHC.Generics     , module Debug.Trace-    -- * Lists-    , takeEvery-    , breakEvery-    -- * Low-Level-    , traceGiven-    -- * Numeric-    , roundSD-    , toPi-    -- ** Functions-    , logistic-    , logit-    -- ** Lists-    , mean-    , range-    , discretizeFunction+    , module System.Directory+    -- * (Re)names+    , NatNumber+    , ByteString+    , orderedHeader     ) where  @@ -42,140 +46,48 @@  -- Re-exports -- -import Goal.Core.Plot hiding (empty,over)+import Goal.Core.Util+import Goal.Core.Project+import Goal.Core.Circuit +import Data.Csv hiding (Parser,Field,header)+import qualified Data.Csv as CSV+import Data.Functor+import Data.Foldable+import Data.Traversable import Data.Ord-import Data.Function-import Data.Monoid hiding (Dual)-import Data.List hiding (sum)+import Data.Function hiding ((&)) import Data.Maybe import Data.Either+import Data.Proxy+import Data.Finite+import Data.Kind (Type)+import Data.Functor.Identity+import Data.Type.Equality -import Control.Applicative+import Control.Applicative hiding (empty) import Control.Arrow hiding ((<+>))-import Control.Monad+import Control.Monad hiding (join) import Control.Monad.ST-import Control.Lens.Type-import Control.Lens.Getter-import Control.Lens.Setter hiding (Identity)-import Control.Lens.TH import Control.Concurrent--import Debug.Trace-import Data.Default.Class-import Numeric------ General Functions ------takeEvery :: Int -> [x] -> [x]--- | Takes every nth element, starting with the head of the list.-takeEvery m = map snd . filter (\(x,_) -> mod x m == 0) . zip [0..]--breakEvery :: Int -> [x] -> [[x]]--- | Break the list up into lists of length n.-breakEvery _ [] = []-breakEvery n xs = take n xs : breakEvery n (drop n xs)--traceGiven :: Show a => a -> a--- | Runs traceShow on the given element.-traceGiven a = traceShow a a------ Numeric ------roundSD :: (Floating x, RealFrac x) => Int -> x -> x--- | Roundest the number to the specified significant digit.-roundSD n x = (/10^n) . fromIntegral . round $ 10^n * x--toPi :: (Floating x, RealFrac x) => x -> x--- | Modulo pi thingy.-toPi x =-    let xpi = x / pi-        n = floor xpi-        f = xpi - fromIntegral n-    in if even n then pi * f else -(pi * (1 - f))--logistic :: Floating x => x -> x--- | A standard sigmoid function.-logistic x = 1 / (1 + exp(negate x))--logit :: Floating x => x -> x--- | The inverse of the logistic.-logit x = log $ x / (1 - x)---- Lists ----mean :: Fractional x => [x] -> x--- | Average value of a list of numbers.-mean = uncurry (/) . foldr (\e (s,c) -> (e+s,c+1)) (0,0)--range :: Double -> Double -> Int -> [Double]--- | Returns n  numbers from mn to mx.-range _ _ 0 = []-range mn mx 1 = [(mn + mx) / 2]-range mn mx n =-    [ x * mx + (1 - x) * mn | x <- (/ (fromIntegral n - 1)) . fromIntegral <$> [0 .. n-1] ]+import Control.DeepSeq hiding (force)+import Control.Monad.Primitive hiding (internal) -discretizeFunction :: Double -> Double -> Int -> (Double -> Double) -> [(Double,Double)]--- | Takes range information in the form of a minimum, maximum, and sample count,--- a function to sample, and returns a list of pairs (x,f(x)) over the specified--- range.-discretizeFunction mn mx n f =-    let rng = range mn mx n-    in zip rng $ f <$> rng+import Options.Applicative --- Graveyard --+import GHC.TypeNats hiding (Mod)+import GHC.Generics (Generic) -{--parMapWH :: (a -> b) -> [a] -> [b]--- | ParMap using rseq. WH stands for Weak Head (normal form).-parMapWH = parMap rseq+import Debug.Trace+import System.Directory+import Numeric hiding (log1p,expm1)+import Numeric.SpecFunctions -parMapDS :: NFData b => (a -> b) -> [a] -> [b]-parMapDS = parMap rdeepseq+import Numeric.Natural+import Data.ByteString (ByteString) -gridSearch-    :: (a -> Double) -- ^ The error function on the model-    -> (x -> a) -- ^ A constructor from the parameter being tested to the model-    -> [x] -- ^ The list of parameter values to test-    -> (x,[(x,Double)]) -- ^ The best parameter and model, and accompanying statistics+type NatNumber = Natural --- | A general implementation of a grid search. Returns a triple where the first--- element is the best parameter, the second is the best model, and the third is a list--- of the errors calculated at each parameter value.-gridSearch errorfun constructor xs =-    let as = constructor <$> xs-        errs = parMapDS errorfun as-        x = fst . minimumBy (comparing snd) $ zip xs errs-    in (x,zip xs errs)+orderedHeader :: [ByteString] -> Header+orderedHeader = CSV.header -iterativeOptimization-    :: (a -> a -> Double) -- ^ Difference measure-    -> Double -- ^ Differential error threshold-    -> Int -- ^ Maximum number of iterations (< 1 is interpreted as infinity)-    -> (a -> a) -- ^ Iterator-    -> a -- ^ Initial Value-    -> (a,[(a,Double)]) -- ^ The final value and associated descent--- | Iterates a value, stopping when a stopping criterion is satisfied, or the maximum--- number of iterations is reached. The algorithm returns the resulting value as well as--- the descent preceding it.------ The error accompanying each element of the descent corresponds to the error between--- that element and the following element. The last error in the descent will therefore--- correspond to the measured difference between the last element of the descent, and the--- returned singleton in the pair.------ If the last error in the returned descent is less then the given threshold, then the--- threshold was reached. Otherwise, n will have been reached as the terminal condition.-iterativeOptimization difference thrsh n iterator a =-    let itrs = iterate iterator a-        itrs' = if n > 0 then take (n+1) itrs else itrs-        zps = zip itrs' $ tail itrs'-        zps' = case break ((< thrsh) . snd) . zip zps $ uncurry difference <$> zps of-                   (zps0,[]) -> zps0-                   (zps0,rst) -> zps0 ++ [head rst]-    in (snd . fst . last $ zps', first fst <$> zps')-    -}

@@ -0,0 +1,215 @@+{-# LANGUAGE Arrows,LambdaCase #-}+-- | A set of functions for working with the 'Arrow' known as a Mealy automata,+-- here referred to as 'Circuit's. Circuits are essentialy a way of building+-- composable fold and iterator operations, where some of the values being+-- processed can be hidden.+module Goal.Core.Circuit+    ( -- * Circuits+    Circuit (Circuit, runCircuit)+    , accumulateFunction+    , accumulateCircuit+    , streamCircuit+    , iterateCircuit+    , loopCircuit+    , loopAccumulator+    , arrM+    -- * Chains+    , Chain+    , chain+    , chainCircuit+    , streamChain+    , iterateChain+    , skipChain+    , skipChain0+    -- ** Recursive Computations+    , iterateM+    , iterateM'+    ) where++--- Imports ---+++-- Unqualified --++import Control.Arrow++-- Qualified --++import qualified Control.Category as C++--- Circuits ---++-- | An arrow which takes an input, monadically produces an output, and updates+-- an (inaccessable) internal state.+newtype Circuit m a b = Circuit+    { runCircuit :: a -> m (b, Circuit m a b) }++-- | Takes a function from a value and an accumulator (e.g. just a sum value or+-- an evolving set of parameters for some model) to a value and an accumulator.+-- The accumulator is then looped back into the function, returning a Circuit+-- from a to b, which updates the accumulator every step.+accumulateFunction :: Monad m => acc -> (a -> acc -> m (b,acc)) -> Circuit m a b+{-# INLINE accumulateFunction #-}+accumulateFunction acc f = Circuit $ \a -> do+    (b,acc') <- f a acc+    return (b,accumulateFunction acc' f)++-- | accumulateCircuit takes a 'Circuit' and an inital value and loops it.+accumulateCircuit :: Monad m => acc -> Circuit m (a,acc) (b,acc) -> Circuit m a b+{-# INLINE accumulateCircuit #-}+accumulateCircuit acc0 mly0 = accumulateFunction (acc0,mly0) $ \a (acc,Circuit crcf) -> do+    ((b,acc'),mly') <- crcf (a,acc)+    return (b,(acc',mly'))++-- | Takes a Circuit and an inital value and loops it, but continues+-- to return both the output and the accumulated value.+loopCircuit :: Monad m => acc -> Circuit m (a,acc) (b,acc) -> Circuit m a (b,acc)+{-# INLINE loopCircuit #-}+loopCircuit acc0 mly0 = accumulateFunction (acc0,mly0) $ \a (acc,Circuit crcf) -> do+    ((b,acc'),mly') <- crcf (a,acc)+    return ((b,acc'),(acc',mly'))++-- | Takes a Circuit which only produces an accumulating value, and loops it.+loopAccumulator :: Monad m => acc -> Circuit m (a,acc) acc -> Circuit m a acc+{-# INLINE loopAccumulator #-}+loopAccumulator acc0 mly0 = accumulateFunction (acc0,mly0) $ \a (acc,Circuit crcf) -> do+    (acc',mly') <- crcf (a,acc)+    return (acc',(acc',mly'))++-- | Feeds a list of inputs into a 'Circuit' and returns the (monadic) list of outputs.+streamCircuit :: Monad m => Circuit m a b -> [a] -> m [b]+{-# INLINE streamCircuit #-}+streamCircuit _ [] = return []+streamCircuit (Circuit mf) (a:as) = do+    (b,crc') <- mf a+    (b :) <$> streamCircuit crc' as++-- | Feeds a list of inputs into a Circuit automata and returns the final+-- monadic output. Throws an error on the empty list.+iterateCircuit :: Monad m => Circuit m a b -> [a] -> m b+{-# INLINE iterateCircuit #-}+iterateCircuit _ [] = error "Empty list fed to iterateCircuit"+iterateCircuit (Circuit mf) [a] = fst <$> mf a+iterateCircuit (Circuit mf) (a:as) = do+    (_,crc') <- mf a+    iterateCircuit crc' as++-- | Turn a monadic function into a circuit.+arrM :: Monad m => (a -> m b) -> Circuit m a b+{-# INLINE arrM #-}+arrM mf = Circuit $ \a -> do+    b <- mf a+    return (b, arrM mf)+++--- Chains ---+++-- | A 'Chain' is an iterator built on a 'Circuit'. 'Chain' constructors are+-- designed to ensure that the first value returned is the initial value of the+-- iterator (this is not entirely trivial).+type Chain m x = Circuit m () x++-- | Creates a 'Chain' from an initial state and a transition function. The+-- first step of the chain returns the initial state, and then continues with+-- generated states.+chain+    :: Monad m+    => x -- ^ The initial state+    -> (x -> m x) -- ^ The transition function+    -> Chain m x -- ^ The resulting 'Chain'+{-# INLINE chain #-}+chain x0 mf = accumulateFunction x0 $ \() x -> do+    x' <- mf x+    return (x,x')++-- | Creates a 'Chain' from an initial state and a transition circuit. The+-- first step of the chain returns the initial state, and then continues with+-- generated states.+chainCircuit+    :: Monad m+    => x -- ^ The initial state+    -> Circuit m x x -- ^ The transition circuit+    -> Chain m x -- ^ The resulting 'Chain'+{-# INLINE chainCircuit #-}+chainCircuit x0 crc = accumulateCircuit x0 $ proc ((),x) -> do+    x' <- crc -< x+    returnA -< (x,x')++-- | Returns the specified number of the given 'Chain's output.+streamChain :: Monad m => Int -> Chain m x -> m [x]+{-# INLINE streamChain #-}+streamChain n chn = streamCircuit chn $ replicate (n+1) ()++-- | Returns the given 'Chain's output at the given index.+iterateChain :: Monad m => Int -> Chain m x -> m x+{-# INLINE iterateChain #-}+iterateChain 0 (Circuit mf) = fst <$> mf ()+iterateChain k (Circuit mf) = mf () >>= iterateChain (k-1) . snd++-- | Modify the given 'Chain' so that it returns the initial value, and then+-- skips the specified number of outputs before producing each subsequent output.+skipChain :: Monad m => Int -> Chain m x -> Chain m x+{-# INLINE skipChain #-}+skipChain n (Circuit mf) = Circuit $ \() -> do+    (x',crc') <- mf ()+    return (x', skipChain0 n crc')++-- | Modify the given 'Chain' so that it skips the specified number of outputs+-- before producing each subsequent output (this skips the initial output too).+skipChain0 :: Monad m => Int -> Chain m x -> Chain m x+{-# INLINE skipChain0 #-}+skipChain0 n crc = Circuit $ \() -> do+    (Circuit mf) <- iterateM' n iterator crc+    (x',crc') <- mf ()+    return (x', skipChain0 n crc')+        where iterator (Circuit mf') = snd <$> mf' ()+++-- | Iterate a monadic action the given number of times, returning the complete+-- sequence of values.+iterateM :: Monad m => Int -> (x -> m x) -> x -> m [x]+{-# INLINE iterateM #-}+iterateM n mf x0 = streamChain n $ chain x0 mf++-- | Iterate a monadic action the given number of times, returning the final value.+iterateM' :: Monad m => Int -> (x -> m x) -> x -> m x+{-# INLINE iterateM' #-}+iterateM' n mf x0 = iterateChain n $ chain x0 mf++++--- Instances ---+++instance Monad m => C.Category (Circuit m) where+    --id :: Circuit a a+    {-# INLINE id #-}+    id = Circuit $ \a -> return (a,C.id)+    --(.) :: Circuit b c -> Circuit a b -> Circuit a c+    {-# INLINE (.) #-}+    (.) = dot+        where dot (Circuit crc1) (Circuit crc2) = Circuit $ \a -> do+                  (b, crcA2') <- crc2 a+                  (c, crcA1') <- crc1 b+                  return (c, crcA1' `dot` crcA2')++instance Monad m => Arrow (Circuit m) where+    --arr :: (a -> b) -> Circuit a b+    {-# INLINE arr #-}+    arr f = Circuit $ \a -> return (f a, arr f)+    --first :: Circuit a b -> Circuit (a,c) (b,c)+    {-# INLINE first #-}+    first (Circuit crcf) = Circuit $ \(a,c) -> do+        (b, crcA') <- crcf a+        return ((b,c), first crcA')++instance Monad m => ArrowChoice (Circuit m) where+    --left :: Circuit a b -> Circuit (Either a c) (Either b c)+    {-# INLINE left #-}+    left crcA@(Circuit crcf) = Circuit $+        \case+          Left a -> do+              (b,crcA') <- crcf a+              return (Left b,left crcA')+          Right c -> return (Right c,left crcA)

@@ -1,264 +0,0 @@-module Goal.Core.Plot-    ( -- * Module Exports-      module Graphics.Rendering.Chart-    , module Data.Colour-    , module Data.Colour.Names-    , module Data.Colour.SRGB.Linear-    , module Graphics.Rendering.Chart.Backend.Cairo-    , module Graphics.Rendering.Chart.Grid-    , module Graphics.Rendering.Chart.Gtk-    , module Graphics.Rendering.Chart.State-    , module Goal.Core.Plot.Contour-    -- * Plots-    -- ** PixMap-    , pixMapPlot-    -- ** Histograms-    , histogramPlot-    , histogramPlot0-    , logHistogramPlot-    , logHistogramPlot0-    -- * Layouts-    -- ** PixMap-    , pixMapLayout-    -- ** Histogram-    , histogramLayout-    , logHistogramLayout-    , histogramLayoutLR-    -- * Util-    , rgbaGradient-    -- * Rendering-    , renderableToAspectWindow-    ) where------ Imports -----import Data.List hiding (sum)--import Control.Monad-import Control.Lens.Getter-import Control.Lens.Setter hiding (Identity)--import Data.Default.Class-import Numeric----- Re-exports ----import Graphics.Rendering.Chart hiding (x0,y0,Point)-import Data.Colour-import Data.Colour.Names-import Data.Colour.SRGB.Linear-import Graphics.Rendering.Chart.Backend.Cairo-import Graphics.Rendering.Chart.State-import Graphics.Rendering.Chart.Grid-import Graphics.Rendering.Chart.Gtk---- Scientific ----import Goal.Core.Plot.Contour---- Qualified ----import qualified Graphics.UI.Gtk as G----- Unqualified ----import Graphics.Rendering.Cairo (liftIO)------ Util -----rgbaGradient :: (Double, Double, Double, Double) -> (Double, Double, Double, Double) -> Int-    -> [AlphaColour Double]--- | Returns an ordered list of colours useful for plotting.-rgbaGradient (rmn,gmn,bmn,amn) (rmx,gmx,bmx,amx) n =-    zipWith (flip withOpacity) [amn,amn + astp .. amx]-    $ zipWith3 rgb [rmn,rmn + rstp .. rmx] [gmn,gmn + gstp .. gmx] [bmn,bmn + bstp .. bmx]-    where rstp = (rmx - rmn) / fromIntegral n-          gstp = (gmx - gmn) / fromIntegral n-          bstp = (bmx - bmn) / fromIntegral n-          astp = (amx - amn) / fromIntegral n------ Plots -------- PixMap ----pixMapPlot :: (Double,Double) -> [[AlphaColour Double]] -> Plot Double Double--- | Returns a pixmap representation of a matrix style set of doubles. Based on the--- defaults, the list of colours are assumed to be in (y,x) coordinates, where the--- origin is at the lower left of the image. If matrix style coordinates are desired,--- The containing layout should be given a reversed y axis, so that the origin is at the--- top left of the image.----  The pair of doubles indicates where corner of the given image should be located.-pixMapPlot (x0,y0) pss = foldr1 joinPlot-    $ concat [ [ boxPlot r c . solidFillStyle $ p | (c,p) <- zip [x0..] ps ] | (r,ps) <- zip [y0..] pss ]-    where boxPlot y x stl = toPlot-              $ plot_fillbetween_style .~ stl-              $ plot_fillbetween_values .~ [(x-0.5,(y-0.5,y+0.5)),(x+0.5,(y-0.5,y+0.5))]-              $ def--pixMapLayout :: Int -> Int -> Layout Double Double -> Layout Double Double--- | A nice base layout for a pixMap, with a box around the pixmap with one 'pixel'--- padding, and a reversed y axis for matrix style coordinates.-pixMapLayout rws0 cls0 lyt =-    layout_top_axis_visibility .~ AxisVisibility True True False-    $ layout_left_axis_visibility .~ AxisVisibility True True False-    $ layout_right_axis_visibility .~ AxisVisibility True True False-    $ layout_bottom_axis_visibility .~ AxisVisibility True True False-    $ layout_y_axis . laxis_reverse .~ True-    $ layout_y_axis . laxis_generate .~-        const (makeAxis (const "") ([-1.5,rws+0.5],[-1.5,rws+0.5],[-1.5,rws+0.5]))-    $ layout_x_axis . laxis_generate .~-        const (makeAxis (const "") ([-1.5,cls+0.5],[-1.5,cls+0.5],[-1.5,cls+0.5]))-    $ lyt-    where cls = fromIntegral cls0-          rws = fromIntegral rws0---- Histogram ----histogramPlot0 :: (Num a,BarsPlotValue a) =>Int -> [[Double]] -> PlotBars Double a -> PlotBars Double a--- | Generates a histogram plot where the min and max bin value is taken from the data set.-histogramPlot0 n xss plt =-    let mx = maximum $ maximum <$> xss-        mn = minimum $ minimum <$> xss-    in histogramPlot n mn mx xss plt--logHistogramPlot0 :: Int -> [[Double]] -> PlotBars Double Double -> PlotBars Double Double--- | Generates a histogram plot where the min and max bin value is taken from the data set.-logHistogramPlot0 n xss =-    let mx = maximum $ maximum <$> xss-        mn = minimum $ minimum <$> xss-    in logHistogramPlot n mn mx xss--histogramPlot-    :: (Num a,BarsPlotValue a)-    => Int -- ^ Number of bins-    -> Double -- ^ Min range-    -> Double -- ^ Max range-    -> [[Double]] -- ^ Data set-    -> PlotBars Double a -- ^ Plot-    -> PlotBars Double a -- ^ New Plot--- | Creates a histogram out of a data set. The data set is a list of list of values, where--- each sublist is a collection of data along an axis. Under and overflow is put into the--- first and last bin, respectively. The bars are centered at the mid point between each--- pair of bins.-histogramPlot n mn mx xss plt =-    let bns = range mn mx (n+1)-        vls = transpose $ toHistogram (tail $ take n bns) . sort <$> xss-        stp = (head (tail bns) - head bns) / 2-        bns' = (+ stp) <$> take n bns-     in plot_bars_alignment .~ BarsCentered $ plot_bars_values .~ zip bns' vls $ plt-    where toHistogram _ [] = repeat 0-          toHistogram [] xs = [genericLength xs]-          toHistogram (bn:bns') xs =-              let (hds,xs') = span (< bn) xs-              in genericLength hds : toHistogram bns' xs'--range _ _ 0 = []-range mn mx 1 = [(mn + mx) / 2]-range mn mx n =-    [ x * mx + (1 - x) * mn | x <- (/ (fromIntegral n - 1)) . fromIntegral <$> [0 .. n-1] ]--logHistogramPlot-    :: Int -- ^ Number of bins-    -> Double -- ^ Min range-    -> Double -- ^ Max range-    -> [[Double]] -- ^ Data set-    -> PlotBars Double Double -- ^ Plot-    -> PlotBars Double Double -- ^ New Plot-logHistogramPlot n mn mx xss plt =-    let bplt = histogramPlot n mn mx xss plt-        vls = modVals <$> bplt ^. plot_bars_values-        --cl = fromIntegral . ceiling . maximum . concat $ snd <$> vls-    in plot_bars_values .~ vls $ bplt-    where lbs = 10-          modVals (x,ys) = (x, logBase lbs . (+1) <$> ys)--histogramLayout :: BarsPlotValue a => PlotBars Double a -> Layout Double a -> Layout Double a--- | The base layout for a histogram.-histogramLayout pbrs lyt =-    let bns = fst <$> pbrs ^. plot_bars_values-        stp = (head (tail bns) - head bns) / 2-        bns' = (head bns - stp) : map (+stp) bns-        rng = abs $ maximum bns'-        labelFun x-            | rng >= 1000 = showEFloat (Just 2) x ""-            | rng <= 0.01 = showEFloat (Just 2) x ""-            | otherwise = reverse . dropWhile (== '.') . dropWhile (== '0') . reverse $ showFFloat (Just 2) x ""-    in layout_x_axis . laxis_generate .~ const (makeAxis labelFun (bns',[],bns'))-       $ layout_plots %~ (plotBars pbrs:)-       $ lyt--logHistogramLayout :: PlotBars Double Double -> Layout Double Double -> Layout Double Double--- | Base layout for a log-histogram.-logHistogramLayout pbrs lyt =-    let vls = pbrs ^. plot_bars_values-        bns = fst <$> vls-        stp = (head (tail bns) - head bns) / 2-        bns' = (head bns - stp) : map (+stp) bns-        cl = fromIntegral . ceiling . maximum . concat $ snd <$> vls-        rng = abs $ maximum bns'-        xLabelFun x-            | rng >= 1000 = showEFloat (Just 2) x ""-            | rng <= 0.01 = showEFloat (Just 2) x ""-            | otherwise = reverse . dropWhile (== '.') . dropWhile (== '0') . reverse $ showFFloat (Just 2) x ""-    in layout_plots %~ (plotBars pbrs:)-        $ layout_y_axis . laxis_generate .~ const (makeAxis yLabelFun ([0..cl],[],[0..cl]))-        $ layout_x_axis . laxis_generate .~ const (makeAxis xLabelFun (bns',[],bns'))-        $ lyt-    where lbs = 10-          yLabelFun 0 = show 0-          yLabelFun x = show (round lbs) ++ "e" ++ show (round x) ++ "-1"--histogramLayoutLR :: (BarsPlotValue a,PlotValue b) => PlotBars Double a -> LayoutLR Double a b -> LayoutLR Double a b--- | The base layout for a histogramLR.-histogramLayoutLR pbrs lyt =-    let bns = fst <$> pbrs ^. plot_bars_values-        stp = (head (tail bns) - head bns) / 2-        bns' = (head bns - stp) : map (+stp) bns-        rng = abs $ maximum bns'-        labelFun x-            | rng >= 1000 = showEFloat (Just 2) x ""-            | rng <= 0.01 = showEFloat (Just 2) x ""-            | otherwise = reverse . dropWhile (== '.') . dropWhile (== '0') . reverse $ showFFloat (Just 2) x ""-    in layoutlr_x_axis . laxis_generate .~ const (makeAxis labelFun (bns',[],bns'))-       $ layoutlr_plots %~ (Left (plotBars pbrs):)-       $ lyt----- IO ------renderableToAspectWindow-    :: Bool -- ^ Display Full Screen-    -> Int -- ^ Image width-    -> Int -- ^ Image height-    -> Renderable a -- ^ The Renderable-    -> IO () -- ^ Renders the renderable to the screen---- | Displays a renderable in a GTK aspect window.-renderableToAspectWindow fs wdth hght rnbl = do--    G.initGUI--    win <- G.windowNew-    afrm <- G.aspectFrameNew 0.5 0.5 . Just $ realToFrac (fromIntegral wdth / fromIntegral hght)-    da <- G.drawingAreaNew--    G.set afrm [ G.containerChild G.:= da ]-    G.set win [ G.containerChild G.:= afrm ]-    when fs $ G.windowFullscreen win-    G.onDestroy win G.mainQuit--    (da `G.on` G.exposeEvent) . liftIO $ do--        updateCanvas rnbl da-        return True--    G.widgetShowAll win-    G.mainGUI-

@@ -1,163 +0,0 @@-module Goal.Core.Plot.Contour (contours) where------ Imports -------- General ----import Control.Monad-import Data.List---- Qualified ----import qualified Data.Map as M----- Unqualified ----import Data.Tuple (swap)------ Contour Plot ------contours-    :: (Double,Double,Int) -- ^ The range along the x axis.-    -> (Double,Double,Int) -- ^ The range along the y axis.-    -> Int -- ^ The number of isolevels.-    -> (Double -> Double -> Double) -- ^ The function to contour.-    -> [(Double,[[(Double,Double)]])] -- ^ (isolevel, (x,y))---- | Given various parameters and a function to analyze, contourPlot returns a--- list of PlotLines each of which is an isoline.-contours (xmn,xmx,nx) (ymn,ymx,ny) nlvls f =-    let lsts = functionToLists f (xmn,ymn) (xmx,ymx) stps-        stps = ((xmx - xmn) / fromIntegral nx,(ymx - ymn) / fromIntegral ny)-        (mn,mx) = (minimum $ minimum <$> lsts,maximum $ maximum <$> lsts)-        isostp = (mx - mn) / fromIntegral nlvls-        isolvls = take nlvls [mn + (isostp / 2),mn + (isostp * 1.5)..]-     in contour lsts (xmx,ymx) stps <$> isolvls------ Internal ------contour :: [[Double]] -> (Double,Double) -> (Double,Double) -> Double -> (Double,[[(Double,Double)]])-contour lsts mxs stps isolvl = (isolvl,linksToLines mxs stps <$> listsToLinks lsts isolvl)------ Contour Plot Internal -------- Types ----type Link = (Int,Int)-type ContourPair = (Link,Link)-data ContourBox =-    {- Styled like the indicies of a 3x3 Matrix. The lines are drawn from-       left to right, and if need be, top to bottom -}-    Empty-    | Line ContourPair-    | DLine ContourPair ContourPair-    deriving (Show,Eq)--data PairMap = PM (M.Map Link Link) (M.Map Link Link) deriving Show---- Type Functions ----contourBoxesToPairMap :: [ContourBox] -> PairMap-contourBoxesToPairMap clns =-    let cprs = concatMap toPairs clns-    in PM (M.fromList cprs) (M.fromList $ map swap cprs)-    where toPairs (Line prpr) = [prpr]-          toPairs (DLine lprpr rprpr) = [lprpr,rprpr]--deletePair :: ContourPair -> PairMap -> PairMap-{- Deletes a left oriented line segmented -}-deletePair (lpr,rpr) (PM lmp rmp) =-    let lmp' = if M.lookup lpr lmp == Just rpr then M.delete lpr lmp else lmp-        rmp' = if M.lookup rpr rmp == Just lpr then M.delete rpr rmp else rmp-    in PM lmp' rmp'--popLink :: Link -> PairMap -> (Maybe Link,PairMap)-popLink lnk pmp@(PM lmp rmp) =-    let lft = M.lookup lnk lmp-        rgt = M.lookup lnk rmp-    in case (lft,rgt) of-           (Just rlnk,_) -> (Just rlnk,deletePair (lnk,rlnk) pmp)-           (_,Just llnk) -> (Just llnk,deletePair (llnk,lnk) pmp)-           _ -> (Nothing,pmp)---- Core Algorithm ----functionToLists :: (Double -> Double -> Double) -> (Double,Double) -> (Double,Double)-    -> (Double,Double) -> [[Double]]-functionToLists f (xmn,ymn) (xmx,ymx) (xstp,ystp) =-    map (uncurry f) <$> [ [ (x,y) | x <- [xmx,(xmx - xstp)..xmn] ] | y <- [ymx,(ymx - ystp)..ymn] ]--situationTable :: Double -> (((Bool,Double),(Bool,Double)),((Bool,Double),(Bool,Double))) -> ContourBox-situationTable _ (((True,_),(True,_)),((True,_),(True,_))) = Empty-situationTable _ (((True,_),(True,_)),((False,_),(True,_))) = Line ((1,0),(2,1))-situationTable _ (((True,_),(True,_)),((True,_),(False,_))) = Line ((2,1),(1,2))-situationTable _ (((True,_),(True,_)),((False,_),(False,_))) = Line ((1,0),(1,2))-situationTable _ (((True,_),(False,_)),((True,_),(True,_))) = Line ((0,1),(1,2))-situationTable isolvl (((True,ul),(False,ur)),((False,ll),(True,lr)))-    | sum [ul,ur,ll,lr] / 4 < isolvl = DLine ((1,0),(0,1)) ((2,1),(1,2))-    | otherwise = DLine ((1,0),(2,1)) ((0,1),(1,2))-situationTable _ (((True,_),(False,_)),((True,_),(False,_))) = Line ((0,1),(2,1))-situationTable _ (((True,_),(False,_)),((False,_),(False,_))) = Line ((1,0),(0,1))-situationTable _ (((False,_),(True,_)),((True,_),(True,_))) = Line ((1,0),(0,1))-situationTable _ (((False,_),(True,_)),((False,_),(True,_))) = Line ((0,1),(2,1))-situationTable isolvl (((False,ul),(True,ur)),((True,ll),(False,lr)))-    | sum [ul,ur,ll,lr] / 4 < isolvl = DLine ((1,0),(2,1)) ((0,1),(1,2))-    | otherwise = DLine ((1,0),(0,1)) ((2,1),(1,2))-situationTable _ (((False,_),(True,_)),((False,_),(False,_))) = Line ((0,1),(1,2))-situationTable _ (((False,_),(False,_)),((True,_),(True,_))) = Line ((1,0),(1,2))-situationTable _ (((False,_),(False,_)),((False,_),(True,_))) = Line ((2,1),(1,2))-situationTable _ (((False,_),(False,_)),((True,_),(False,_))) = Line ((1,0),(2,1))-situationTable _ (((False,_),(False,_)),((False,_),(False,_))) = Empty--listsToContourBoxes :: [[Double]] -> Double -> [ContourBox]-listsToContourBoxes lsts isolvl = do-    (stnrw,r) <- zip stnlsts [0..]-    (stn,c) <- zip stnrw [0..]-    let bx = situationTable isolvl stn-    guard (bx /= Empty)-    return $ repositionBox (r,c) bx-    where stnlsts = rowZipper $ elementPairs . map threshold <$> lsts-          threshold x = (x >= isolvl,x)-          elementPairs lst = zip lst $ tail lst-          rowZipper rws = uncurry zip <$> elementPairs rws-          repositionContourPair (r,c) ((lr,lc),(rr,rc)) =-              ((lr + 2 * r, lc + 2 * c),(rr + 2 * r, rc + 2 * c))-          repositionBox rc (Line pr) =-              Line $ repositionContourPair rc pr-          repositionBox rc (DLine pr1 pr2) =-              DLine (repositionContourPair rc pr1) (repositionContourPair rc pr2)--traceLinks :: ContourPair -> PairMap -> ([Link],PairMap)-traceLinks (llnk,rlnk) pmp =-    let (llnks,pmp') = tracer [] (Just llnk,pmp)-        (rlnks,pmp'') = tracer [] (Just rlnk,pmp')-    in (llnks ++ reverse rlnks,pmp'')-    where tracer lnks (Just lnk,pmp') =-              tracer (lnk:lnks) $ popLink lnk pmp'-          tracer lnks (Nothing,pmp') = (lnks,pmp')--popLinks :: PairMap -> Maybe ([Link],PairMap)-popLinks pmp@(PM lmp _)-    | M.null lmp = Nothing-    | otherwise =-        let cpr = M.findMin lmp-        in Just (traceLinks cpr $ deletePair cpr pmp)--listsToLinks :: [[Double]] -> Double -> [[(Int,Int)]]-listsToLinks lsts isolvl =-    unfoldr popLinks .  contourBoxesToPairMap $ listsToContourBoxes lsts isolvl--linksToLines :: (Double,Double) -> (Double,Double) -> [(Int,Int)] -> [(Double,Double)]-linksToLines (xmx,ymx) (xstp,ystp) lnks =-   (\(r,c) -> (xmx - fromIntegral c * xstp / 2,ymx - fromIntegral r * ystp / 2)) <$> lnks--

@@ -0,0 +1,185 @@+{-# LANGUAGE OverloadedStrings #-}+-- | This module provides functions for incorporating Goal into a+-- data-processing project. In particular, this module provides tools for+-- managing CSV files, and connecting them with gnuplot scripts for plotting.+-- CSV management is powered by @cassava@.+module Goal.Core.Project+    (+    -- * CSV+      goalImport+    , goalImportNamed+    , goalExport+    , goalExportLines+    , goalExportNamed+    , goalExportNamedLines+    -- ** CSV Instances+    , goalCSVParser+    , goalCSVNamer+    , goalCSVOrder+    -- * Util+    , runGnuplot+    , runGnuplotWithVariables+    ) where+++--- Imports ---+++-- Unqualified --++import System.Process+import System.Directory+import Data.Csv+import Data.Char+import GHC.Generics++-- Qualified --++import qualified Data.Vector as V+import qualified Data.ByteString.Lazy as BS+import qualified Data.ByteString as BSI+++--- Experiments ---+++--- Import/Export ---+++-- | Runs @gnuplot@ on the given @.gpi@, passing it a @load_path@ variable to+-- help it find Goal-generated csvs.+runGnuplot+    :: FilePath -- ^ Gnuplot loadpath+    -> String -- ^ Gnuplot script+    -> IO ()+runGnuplot ldpth gpipth = do+    let cmd = concat [ "gnuplot ", " -e \"load_path='", ldpth, "'\" ",gpipth,".gpi" ]+    putStrLn $ "Running Command: " ++ cmd+    callCommand cmd++-- | Runs @gnuplot@ on the given @.gpi@, passing it a @load_path@ variable to+-- help it find Goal-generated csvs, and a list of variables.+runGnuplotWithVariables+    :: FilePath -- ^ Gnuplot loadpath+    -> String -- ^ Gnuplot script+    -> [(String,String)] -- ^ Arguments+    -> IO ()+runGnuplotWithVariables ldpth gpipth args = do+    let cmd = concat $ [ "gnuplot ", " -e \"load_path='", ldpth, "'" ]+            ++ (mapArgs <$> args) ++ [ "\" ",gpipth,".gpi" ]+    putStrLn $ "Running Command: " ++ cmd+    callCommand cmd+        where mapArgs (nm,val) = concat ["; ",nm,"='",val,"'"]+++-- | Load the given CSV file. The @.csv@ extension is automatically added.+goalImport+    :: FromRecord r+    => FilePath+    -> IO (Either String [r]) -- ^ CSVs+goalImport flpth = do+    bstrm <- decode NoHeader <$> BS.readFile (flpth ++ ".csv")+    case bstrm of+      Right as -> return . Right $ V.toList as+      Left str -> return $ Left str++-- | Load the given CSV file with headers. The @.csv@ extension is automatically added.+goalImportNamed+    :: FromNamedRecord r+    => FilePath+    -> IO (Either String [r]) -- ^ CSVs+goalImportNamed flpth = do+    bstrm <- decodeByName <$> BS.readFile (flpth ++ ".csv")+    case bstrm of+      Right as -> return . Right . V.toList $ snd as+      Left str -> return $ Left str++filePather :: FilePath -> FilePath -> IO FilePath+filePather ldpth flnm = do+    createDirectoryIfMissing True ldpth+    return $ concat [ldpth,"/",flnm,".csv"]++-- | Export the given CSVs to a file in the given directory. The @.csv@+-- extension is automatically added to the file name.+goalExport+    :: ToRecord r+    => FilePath -- load_path+    -> String -- File Name+    -> [r] -- ^ CSVs+    -> IO ()+goalExport ldpth flnm csvs = do+    flpth <- filePather ldpth flnm+    BS.writeFile flpth $ encode csvs++-- | Export the given list of CSVs to a file in the given directory, seperating+-- each set of CSVs by a single line. This causes gnuplot to the read CSV as a+-- collection of line segments. The @.csv@ extension is automatically added to+-- the file name.+goalExportLines+    :: ToRecord r+    => FilePath+    -> FilePath+    -> [[r]] -- ^ CSVss+    -> IO ()+goalExportLines ldpth flnm csvss = do+    flpth <- filePather ldpth flnm+    BS.writeFile flpth . BS.tail . BS.tail . BS.concat $ BS.append "\r\n" . encode <$> csvss++-- | Export the named CSVs to a file in the given directory, adding a header to+-- the @.csv@ file.+goalExportNamed+    :: (ToNamedRecord r, DefaultOrdered r)+    => FilePath+    -> FilePath+    -> [r] -- ^ CSVs+    -> IO ()+goalExportNamed ldpth flnm csvs = do+    flpth <- filePather ldpth flnm+    BS.writeFile flpth $ encodeDefaultOrderedByName csvs++-- | Export the given list of named CSVs to a file, breaking it into a set of+-- line segments (with headers).+goalExportNamedLines+    :: (ToNamedRecord r, DefaultOrdered r)+    => FilePath+    -> FilePath+    -> [[r]] -- ^ CSVss+    -> IO ()+goalExportNamedLines ldpth flnm csvss = do+    flpth <- filePather ldpth flnm+    BS.writeFile flpth . BS.concat $ BS.append "\r\n" . encodeDefaultOrderedByName <$> csvss+++--- Util ---+++deCamelCaseLoop :: String -> String+deCamelCaseLoop "" = ""+deCamelCaseLoop (c:wrds) =+    let (wrd,wrds') = span isLower wrds+     in (c:wrd) ++ ' ' : deCamelCaseLoop wrds'++deCamelCase :: String -> String+deCamelCase (c:wrds) = init $ deCamelCaseLoop (toUpper c : wrds)+deCamelCase "" = error "How is deCamelCase being run on an empty string?"++deCamelCaseCSV :: Options+deCamelCaseCSV = defaultOptions { fieldLabelModifier = deCamelCase }++-- | A generic @.csv@ parser which reorganizes a header name in camel case into+-- "human readable" text. Useful for instantiating 'FromNamedRecord'.+goalCSVParser :: (Generic a, GFromNamedRecord (Rep a)) => NamedRecord -> Parser a+goalCSVParser = genericParseNamedRecord deCamelCaseCSV++-- | A generic @.csv@ namer which reorganizes a header name in camel case into+-- "human readable" text. Useful for instantiating 'ToNamedRecord'.+goalCSVNamer+    :: (Generic a, GToRecord (Rep a) (BSI.ByteString, BSI.ByteString)) => a -> NamedRecord+goalCSVNamer = genericToNamedRecord deCamelCaseCSV++-- | A generic @.csv@ order which reorganizes a header name in camel case into+-- "human readable" text. Useful for instantiating 'DefaultOrdered'.+goalCSVOrder :: (Generic a, GToNamedRecordHeader (Rep a)) => a -> Header+goalCSVOrder = genericHeaderOrder deCamelCaseCSV++

@@ -0,0 +1,265 @@+-- | A collection of generic numerical and list manipulation functions.+module Goal.Core.Util+    ( -- * List Manipulation+      takeEvery+    , breakEvery+    , kFold+    , kFold'+    -- * Numeric+    , roundSD+    , toPi+    , circularDistance+    , integrate+    , logistic+    , logit+    , square+    , triangularNumber+    -- ** List Numerics+    , average+    , weightedAverage+    , circularAverage+    , weightedCircularAverage+    , range+    , discretizeFunction+    , logSumExp+    , logIntegralExp+    -- * Tracing+    , traceGiven+    -- * TypeNats+    , finiteInt+    , natValInt+    , Triangular+    -- ** Type Rationals+    , Rat+    , type (/)+    , ratVal+    ) where+++--- Imports ---+++-- Unqualified --++import Numeric+import Data.Ratio+import Data.Proxy+import Debug.Trace+import Data.Finite+import GHC.TypeNats++-- Qualified --++import qualified Numeric.GSL.Integration as I+import qualified Data.List as L++--- General Functions ---+++-- | Takes every nth element, starting with the head of the list.+takeEvery :: Int -> [x] -> [x]+{-# INLINE takeEvery #-}+takeEvery m = map snd . filter (\(x,_) -> mod x m == 0) . zip [0..]++-- | Break the list up into lists of length n.+breakEvery :: Int -> [x] -> [[x]]+{-# INLINE breakEvery #-}+breakEvery _ [] = []+breakEvery n xs = take n xs : breakEvery n (drop n xs)++-- | Runs traceShow on the given element.+traceGiven :: Show a => a -> a+{-# INLINE traceGiven #-}+traceGiven a = traceShow a a+++--- Numeric ---++-- | Numerically integrates a 1-d function over an interval.+integrate+    :: Double -- ^ Error Tolerance+    -> (Double -> Double) -- ^ Function+    -> Double -- ^ Interval beginning+    -> Double -- ^ Interval end+    -> (Double,Double) -- ^ Integral+{-# INLINE integrate #-}+integrate errbnd = I.integrateQAGS errbnd 10000++-- | Rounds the number to the specified significant digit.+roundSD :: RealFloat x => Int -> x -> x+{-# INLINE roundSD #-}+roundSD n x =+    let n' :: Int+        n' = round $ 10^n * x+     in fromIntegral n'/10^n++-- | Value of a point on a circle, minus rotations.+toPi :: RealFloat x => x -> x+{-# INLINE toPi #-}+toPi x =+    let xpi = x / (2*pi)+        f = xpi - fromIntegral (floor xpi :: Int)+     in 2 * pi * f++-- | Distance between two points on a circle, removing rotations.+circularDistance :: RealFloat x => x -> x -> x+{-# INLINE circularDistance #-}+circularDistance x y =+    let x' = toPi x+        y' = toPi y+     in min (toPi $ x' - y') (toPi $ y' - x')++-- | A standard sigmoid function.+logistic :: Floating x => x -> x+{-# INLINE logistic #-}+logistic x = 1 / (1 + exp (negate x))++-- | The inverse of the logistic.+logit :: Floating x => x -> x+{-# INLINE logit #-}+logit x = log $ x / (1 - x)++-- | The square of a number (for avoiding endless default values).+square :: Floating x => x -> x+{-# INLINE square #-}+square x = x^(2::Int)++-- | Triangular number.+triangularNumber :: Int -> Int+{-# INLINE triangularNumber #-}+triangularNumber n = flip div 2 $ n * (n+1)+++-- Lists --++-- | Average value of a 'Traversable' of 'Fractional's.+average :: (Foldable f, Fractional x) => f x -> x+{-# INLINE average #-}+average = uncurry (/) . L.foldl' (\(s,c) e -> (e+s,c+1)) (0,0)++-- | Weighted Average given a 'Traversable' of (weight,value) pairs.+weightedAverage :: (Foldable f, Fractional x) => f (x,x) -> x+{-# INLINE weightedAverage #-}+weightedAverage = uncurry (/) . L.foldl' (\(sm,nrm) (w,x) -> (sm + w*x,nrm + w)) (0,0)++-- | Circular average value of a 'Traversable' of radians.+circularAverage :: (Traversable f, RealFloat x) => f x -> x+{-# INLINE circularAverage #-}+circularAverage rds =+    let snmu = average $ sin <$> rds+        csmu = average $ cos <$> rds+     in atan2 snmu csmu++-- | Returns k (training,validation) pairs. k should be greater than or equal to 2.+kFold :: Int -> [x] -> [([x],[x])]+{-# INLINE kFold #-}+kFold k xs =+    let nvls = ceiling . (/(fromIntegral k :: Double)) . fromIntegral $ length xs+     in L.unfoldr unfoldFun ([], breakEvery nvls xs)+    where unfoldFun (_,[]) = Nothing+          unfoldFun (hds,tl:tls) = Just ((concat $ hds ++ tls,tl),(tl:hds,tls))++-- | Returns k (training,test,validation) pairs for early stopping algorithms. k+-- should be greater than or equal to 3.+kFold' :: Int -> [x] -> [([x],[x],[x])]+{-# INLINE kFold' #-}+kFold' k xs =+    let nvls = ceiling . (/(fromIntegral k :: Double)) . fromIntegral $ length xs+        brks = breakEvery nvls xs+     in L.unfoldr unfoldFun ([], brks)+    where unfoldFun (hds,tl:tl':tls) = Just ((concat $ hds ++ tls,tl,tl'),(tl:hds,tl':tls))+          unfoldFun (hds,tl:tls) =+              let (tl0:hds') = reverse hds+               in Just ((concat $ reverse hds' ++ tls,tl,tl0),(tl:hds,tls))+          unfoldFun (_,[]) = Nothing+++-- | Weighted Circular average value of a 'Traversable' of radians.+weightedCircularAverage :: (Traversable f, RealFloat x) => f (x,x) -> x+{-# INLINE weightedCircularAverage #-}+weightedCircularAverage wxs =+    let snmu = weightedAverage $ sinPair <$> wxs+        csmu = weightedAverage $ cosPair <$> wxs+     in atan2 snmu csmu+    where sinPair (w,rd) = (w,sin rd)+          cosPair (w,rd) = (w,cos rd)++-- | Returns n numbers which uniformly partitions the interval [mn,mx].+range+    :: RealFloat x => x -> x -> Int -> [x]+{-# INLINE range #-}+range _ _ 0 = []+range mn mx 1 = [(mn + mx) / 2]+range mn mx n =+    [ x * mx + (1 - x) * mn | x <- (/ (fromIntegral n - 1)) . fromIntegral <$> [0 .. n-1] ]++-- | Takes range information in the form of a minimum, maximum, and sample count,+-- a function to sample, and returns a list of pairs (x,f(x)) over the specified+-- range.+discretizeFunction :: Double -> Double -> Int -> (Double -> Double) -> [(Double,Double)]+{-# INLINE discretizeFunction #-}+discretizeFunction mn mx n f =+    let rng = range mn mx n+    in zip rng $ f <$> rng++-- | Given a set of values, computes the "soft maximum" by way of taking the+-- exponential of every value, summing the results, and then taking the+-- logarithm. Incorporates some tricks to improve numerical stability.+logSumExp :: (Ord x, Floating x, Traversable f) => f x -> x+{-# INLINE logSumExp #-}+logSumExp xs =+    let mx = maximum xs+     in (+ mx) . log1p . subtract 1 . sum $ exp . subtract mx <$> xs++-- | Given a function, computes the "soft maximum" of the function by computing+-- the integral of the exponential of the function, and taking the logarithm of+-- the result. The maximum is first approximated on a given set of samples to+-- improve numerical stability. Pro tip: If you want to compute the normalizer+-- of a an exponential family probability density, provide the unnormalized+-- log-density to this function.+logIntegralExp+    :: Traversable f+    => Double -- ^ Error Tolerance+    -> (Double -> Double) -- ^ Function+    -> Double -- ^ Interval beginning+    -> Double -- ^ Interval end+    -> f Double -- ^ Samples (for approximating the max)+    -> Double -- ^ Log-Integral-Exp+{-# INLINE logIntegralExp #-}+logIntegralExp err f mnbnd mxbnd xsmps =+    let mx = maximum $ f <$> xsmps+        expf x = exp $ f x - mx+     in (+ mx) . log1p . subtract 1 . fst $ integrate err expf mnbnd mxbnd+++--- TypeLits ---+++-- | Type-level triangular number.+type Triangular n = Div (n * (n + 1)) 2++-- | Type level rational numbers. This implementation does not currently permit negative numbers.+data Rat (n :: Nat) (d :: Nat)++-- | Infix 'Rat'.+type (/) n d = Rat n d++-- | Recover a rational value from a 'Proxy'.+ratVal :: (KnownNat n, KnownNat d) => Proxy (n / d) -> Rational+{-# INLINE ratVal #-}+ratVal = ratVal0 Proxy Proxy+++-- | 'natVal and 'fromIntegral'.+natValInt :: KnownNat n => Proxy n -> Int+{-# INLINE natValInt #-}+natValInt = fromIntegral . natVal++-- | 'getFinite' and 'fromIntegral'.+finiteInt :: KnownNat n => Finite n -> Int+{-# INLINE finiteInt #-}+finiteInt = fromIntegral . getFinite++ratVal0 :: (KnownNat n, KnownNat d) => Proxy n -> Proxy d -> Proxy (n / d) -> Rational+{-# INLINE ratVal0 #-}+ratVal0 prxyn prxyd _ = fromIntegral (natVal prxyn) % fromIntegral (natVal prxyd)

@@ -0,0 +1,245 @@+ {-# OPTIONS_GHC -fplugin=GHC.TypeLits.KnownNat.Solver -fplugin=GHC.TypeLits.Normalise -fconstraint-solver-iterations=10 #-}+-- | Vectors and Matrices with statically-typed dimensions based on boxed vectors.++module Goal.Core.Vector.Boxed+    ( -- * Vector+      module Data.Vector.Sized+      -- ** Construction+    , doubleton+    , range+    , breakStream+    , breakEvery+    -- ** Deconstruction+    , toPair+    , concat+    -- * Matrix+    , Matrix+    , nRows+    , nColumns+    -- ** Construction+    , fromRows+    , fromColumns+    , matrixIdentity+    , outerProduct+    , diagonalConcat+    -- ** Deconstruction+    , toRows+    , toColumns+    -- ** Manipulation+    , columnVector+    , rowVector+    -- ** BLAS+    , dotProduct+    , matrixVectorMultiply+    , matrixMatrixMultiply+    , inverse+    , transpose+    ) where+++--- Imports ---+++-- Goal --++import Goal.Core.Util hiding (breakEvery,range)+import qualified Goal.Core.Util (breakEvery)++import qualified Goal.Core.Vector.Generic as G++-- Unqualified --++import Prelude hiding (concat,zipWith,(++),replicate)+import qualified Data.Vector as B+import qualified Data.Vector.Mutable as BM+import qualified Control.Monad.ST as ST+import qualified Data.Vector.Generic.Sized.Internal as I++-- Qualified --++import Data.Vector.Sized+import GHC.TypeNats+import Data.Proxy++-- Qualified Imports --++-- | Flatten a 'Vector' of 'Vector's.+concat :: KnownNat n => Vector m (Vector n x) -> Vector (m*n) x+{-# INLINE concat #-}+concat = G.concat++-- | Create a 'Vector' of length 2.+doubleton :: x -> x -> Vector 2 x+{-# INLINE doubleton #-}+doubleton = G.doubleton++-- | Partition of an interval.+range :: (KnownNat n, Fractional x) => x -> x -> Vector n x+{-# INLINE range #-}+range = G.range++-- | Cycles a list of elements and breaks it up into an infinite list of 'Vector's.+breakStream :: forall n a. KnownNat n => [a] -> [Vector n a]+{-# INLINE breakStream #-}+breakStream as =+    I.Vector . B.fromList <$> Goal.Core.Util.breakEvery (natValInt (Proxy :: Proxy n)) (cycle as)++-- | Converts a length two 'Vector' into a pair of elements.+toPair :: Vector 2 x -> (x,x)+{-# INLINE toPair #-}+toPair = G.toPair++-- | Breaks a 'Vector' into a Vector of Vectors.+breakEvery :: (KnownNat n, KnownNat k) => Vector (n*k) a -> Vector n (Vector k a)+{-# INLINE breakEvery #-}+breakEvery = G.breakEvery+++--- Matrices ---+++-- | Matrices with static dimensions (boxed).+type Matrix = G.Matrix B.Vector++-- | The number of rows in the 'Matrix'.+nRows :: forall m n a . KnownNat m => Matrix m n a -> Int+{-# INLINE nRows #-}+nRows = G.nRows++-- | The number of columns in the 'Matrix'.+nColumns :: forall m n a . KnownNat n => Matrix m n a -> Int+{-# INLINE nColumns #-}+nColumns = G.nColumns++-- | Convert a 'Matrix' into a 'Vector' of 'Vector's of rows.+toRows :: (KnownNat m, KnownNat n) => Matrix m n x -> Vector m (Vector n x)+{-# INLINE toRows #-}+toRows = G.toRows++-- | Convert a 'Matrix' into a 'Vector' of 'Vector's of columns.+toColumns :: (KnownNat m, KnownNat n) => Matrix m n x -> Vector n (Vector m x)+{-# INLINE toColumns #-}+toColumns = G.toColumns+++-- | Turn a 'Vector' into a single column 'Matrix'.+columnVector :: Vector n a -> Matrix n 1 a+{-# INLINE columnVector #-}+columnVector = G.columnVector++-- | Turn a 'Vector' into a single row 'Matrix'.+rowVector :: Vector n a -> Matrix 1 n a+{-# INLINE rowVector #-}+rowVector = G.rowVector++-- | Create a 'Matrix' from a 'Vector' of row 'Vector's.+fromRows :: KnownNat n => Vector m (Vector n x) -> Matrix m n x+{-# INLINE fromRows #-}+fromRows = G.fromRows++-- | Create a 'Matrix' from a 'Vector' of column 'Vector's.+fromColumns :: (KnownNat n, KnownNat m) => Vector n (Vector m x) -> Matrix m n x+{-# INLINE fromColumns #-}+fromColumns = G.fromColumns++-- | Diagonally concatenate two matrices, padding the gaps with zeroes (pure implementation).+diagonalConcat+    :: (KnownNat n, KnownNat m, KnownNat o, KnownNat p, Num a)+    => Matrix n m a -> Matrix o p a -> Matrix (n+o) (m+p) a+{-# INLINE diagonalConcat #-}+diagonalConcat mtx1 mtx2 =+    let rws1 = (++ replicate 0) <$> toRows mtx1+        rws2 = (replicate 0 ++) <$> toRows mtx2+     in fromRows $ rws1 ++ rws2++-- | Pure implementation of the dot product.+dotProduct :: Num x => Vector n x -> Vector n x -> x+{-# INLINE dotProduct #-}+dotProduct = G.dotProduct++-- | Pure implementation of the outer product.+outerProduct+    :: (KnownNat m, KnownNat n, Num x)+    => Vector m x -> Vector n x -> Matrix m n x+{-# INLINE outerProduct #-}+outerProduct = G.outerProduct++-- | Pure implementation of 'Matrix' transposition.+transpose+    :: (KnownNat m, KnownNat n, Num x)+    => Matrix m n x -> Matrix n m x+{-# INLINE transpose #-}+transpose = G.transpose++-- | Pure 'Matrix' x 'Vector' multiplication.+matrixVectorMultiply+    :: (KnownNat m, KnownNat n, Num x)+    => Matrix m n x -> Vector n x -> Vector m x+{-# INLINE matrixVectorMultiply #-}+matrixVectorMultiply mtx = G.toVector . matrixMatrixMultiply mtx . columnVector++-- | The identity 'Matrix'.+matrixIdentity :: (KnownNat n, Num a) => Matrix n n a+{-# INLINE matrixIdentity #-}+matrixIdentity =+    fromRows $ generate (\i -> generate (\j -> if finiteInt i == finiteInt j then 1 else 0))++-- | Pure implementation of matrix inversion.+inverse :: forall a n. (Fractional a, Ord a, KnownNat n) => Matrix n n a -> Maybe (Matrix n n a)+{-# INLINE inverse #-}+inverse mtx =+    let rws = fromSized $ fromSized <$> zipWith (++) (toRows mtx) (toRows matrixIdentity)+        n = natValInt (Proxy :: Proxy n)+        rws' = B.foldM' eliminateRow rws $ B.generate n id+     in G.Matrix . I.Vector . B.concatMap (B.drop n) <$> rws'++-- | Pure 'Matrix' x 'Matrix' multiplication.+matrixMatrixMultiply+    :: forall m n o a. (KnownNat m, KnownNat n, KnownNat o, Num a)+    => Matrix m n a -> Matrix n o a -> Matrix m o a+{-# INLINE matrixMatrixMultiply #-}+matrixMatrixMultiply (G.Matrix (I.Vector v)) wm =+    let n = natValInt (Proxy :: Proxy n)+        o = natValInt (Proxy :: Proxy o)+        (G.Matrix (I.Vector w')) = G.transpose wm+        f k = let (i,j) = divMod (finiteInt k) o+                  slc1 = B.unsafeSlice (i*n) n v+                  slc2 = B.unsafeSlice (j*n) n w'+               in G.weakDotProduct slc1 slc2+     in G.Matrix $ G.generate f+++--- Internal ---+++eliminateRow :: (Ord a, Fractional a) => B.Vector (B.Vector a) -> Int -> Maybe (B.Vector (B.Vector a))+eliminateRow mtx k = do+    mtx' <- pivotRow k mtx+    return . nullifyRows k $ normalizePivot k mtx'++pivotRow :: (Fractional a, Ord a) => Int -> B.Vector (B.Vector a) -> Maybe (B.Vector (B.Vector a))+pivotRow k rws =+    let l = (+k) . B.maxIndex $ abs . flip B.unsafeIndex k . B.take (B.length rws) <$> B.drop k rws+        ak = B.unsafeIndex rws k B.! l+     in if abs ak < 1e-10 then Nothing+                  else ST.runST $ do+                           mrws <- B.thaw rws+                           BM.unsafeSwap mrws k l+                           Just <$> B.freeze mrws++normalizePivot :: Fractional a => Int -> B.Vector (B.Vector a) -> B.Vector (B.Vector a)+normalizePivot k rws = ST.runST $ do+    let ak = recip . flip B.unsafeIndex k $ B.unsafeIndex rws k+    mrws <- B.thaw rws+    BM.modify mrws ((*ak) <$>) k+    B.freeze mrws++nullifyRows :: Fractional a => Int -> B.Vector (B.Vector a) -> B.Vector (B.Vector a)+nullifyRows k rws =+    let rwk = B.unsafeIndex rws k+        ak = B.unsafeIndex rwk k+        generator i = if i == k then 0 else B.unsafeIndex (B.unsafeIndex rws i) k / ak+        as = B.generate (B.length rws) generator+     in B.zipWith (B.zipWith (-)) rws $ (\a -> (*a) <$> rwk) <$> as++

@@ -0,0 +1,233 @@+{-# LANGUAGE StandaloneDeriving,GeneralizedNewtypeDeriving #-}+ {-# OPTIONS_GHC -fplugin=GHC.TypeLits.KnownNat.Solver -fplugin=GHC.TypeLits.Normalise -fconstraint-solver-iterations=10 #-}+-- | Vectors and Matrices with statically typed dimensions.+module Goal.Core.Vector.Generic+    ( -- * Vector+      module Data.Vector.Generic.Sized+    , VectorClass+    -- * Construction+    , doubleton+    , range+    , breakEvery+    -- * Deconstruction+    , concat+    -- * Matrix+    , Matrix (Matrix,toVector)+    , nRows+    , nColumns+    -- ** Construction+    , fromRows+    , fromColumns+    -- ** Deconstruction+    , toPair+    , toRows+    , toColumns+    -- ** Manipulation+    , columnVector+    , rowVector+    -- ** BLAS+    , transpose+    , dotProduct+    , weakDotProduct+    , outerProduct+    , matrixVectorMultiply+    , matrixMatrixMultiply+    ) where+++--- Imports ---+++-- Goal --++import Goal.Core.Util hiding (breakEvery,range)++-- Unqualified --++import GHC.TypeNats+import Data.Proxy+import Control.DeepSeq+import Data.Vector.Generic.Sized+import Data.Vector.Generic.Sized.Internal+import Foreign.Storable+import Prelude hiding (concatMap,concat,map,sum,replicate)++-- Qualified --++import qualified Data.Vector.Generic as G+import qualified Data.Vector.Storable as S++import Numeric.LinearAlgebra (Numeric)++--- Vector ---+++type VectorClass = G.Vector++-- | Create a 'Matrix' from a 'Vector' of 'Vector's which represent the rows.+concat :: (KnownNat n, G.Vector v x, G.Vector v (Vector v n x)) => Vector v m (Vector v n x) -> Vector v (m*n) x+{-# INLINE concat #-}+concat = concatMap id++-- | Collect two values into a length 2 'Vector'.+doubleton :: G.Vector v x => x -> x -> Vector v 2 x+{-# INLINE doubleton #-}+doubleton x1 x2 = cons x1 $ singleton x2++-- | Breaks a 'Vector' into a Vector of Vectors.+breakEvery+    :: forall v n k a . (G.Vector v a, G.Vector v (Vector v k a), KnownNat n, KnownNat k)+    => Vector v (n*k) a -> Vector v n (Vector v k a)+{-# INLINE breakEvery #-}+breakEvery v0 =+    let k = natValInt (Proxy :: Proxy k)+        v = fromSized v0+     in generate (\i -> Vector $ G.unsafeSlice (finiteInt i*k) k v)++-- | Reshapes a length 2 'Vector' into a pair of values.+toPair :: G.Vector v a => Vector v 2 a -> (a,a)+{-# INLINE toPair #-}+toPair v = (unsafeIndex v 0, unsafeIndex v 1)++-- | Uniform partition of an interval into a 'Vector'.+range+    :: forall v n x. (G.Vector v x, KnownNat n, Fractional x)+    => x -> x -> Vector v n x+{-# INLINE range #-}+range mn mx =+    let n = natValInt (Proxy :: Proxy n)+        stp = (mx - mn)/fromIntegral (n-1)+     in enumFromStepN mn stp+++--- Matrix ---+++-- | Matrices with static dimensions.+newtype Matrix v (m :: Nat) (n :: Nat) a = Matrix { toVector :: Vector v (m*n) a }+    deriving (Eq,Show,NFData)++deriving instance (KnownNat m, KnownNat n, Storable x) => Storable (Matrix S.Vector m n x)+deriving instance (KnownNat m, KnownNat n, Numeric x, Num x)+  => Num (Matrix S.Vector m n x)+deriving instance (KnownNat m, KnownNat n, Numeric x, Fractional x)+  => Fractional (Matrix S.Vector m n x)+deriving instance (KnownNat m, KnownNat n, Numeric x, Floating x)+  => Floating (Matrix S.Vector m n x)++-- | Turn a 'Vector' into a single column 'Matrix'.+columnVector :: Vector v n a -> Matrix v n 1 a+{-# INLINE columnVector #-}+columnVector = Matrix++-- | Turn a 'Vector' into a single row 'Matrix'.+rowVector :: Vector v n a -> Matrix v 1 n a+{-# INLINE rowVector #-}+rowVector = Matrix++-- | Create a 'Matrix' from a 'Vector' of 'Vector's which represent the rows.+fromRows :: (G.Vector v x, G.Vector v (Vector v n x), KnownNat n) => Vector v m (Vector v n x) -> Matrix v m n x+{-# INLINE fromRows #-}+fromRows = Matrix . concat++-- | Create a 'Matrix' from a 'Vector' of 'Vector's which represent the columns.+fromColumns+    :: (G.Vector v x, G.Vector v Int, G.Vector v (Vector v n x), G.Vector v (Vector v m x), KnownNat n, KnownNat m)+    => Vector v n (Vector v m x) -> Matrix v m n x+{-# INLINE fromColumns #-}+fromColumns = transpose . fromRows++-- | The number of rows in the 'Matrix'.+nRows :: forall v m n a . KnownNat m => Matrix v m n a -> Int+{-# INLINE nRows #-}+nRows _ = natValInt (Proxy :: Proxy m)++-- | The number of columns in the 'Matrix'.+nColumns :: forall v m n a . KnownNat n => Matrix v m n a -> Int+{-# INLINE nColumns #-}+nColumns _ = natValInt (Proxy :: Proxy n)++-- | Convert a 'Matrix' into a 'Vector' of 'Vector's of rows.+toRows :: (G.Vector v a, G.Vector v (Vector v n a), KnownNat n, KnownNat m)+       => Matrix v m n a -> Vector v m (Vector v n a)+{-# INLINE toRows #-}+toRows (Matrix v) = breakEvery v++-- | Convert a 'Matrix' into a 'Vector' of 'Vector's of columns.+toColumns+    :: (G.Vector v a, G.Vector v (Vector v m a), KnownNat m, KnownNat n, G.Vector v Int)+    => Matrix v m n a -> Vector v n (Vector v m a)+{-# INLINE toColumns #-}+toColumns = toRows . transpose+++--- BLAS ---+++-- | Pure implementation of 'Matrix' transposition.+transpose+    :: forall v m n a . (KnownNat m, KnownNat n, G.Vector v Int, G.Vector v a, G.Vector v (Vector v m a))+    => Matrix v m n a -> Matrix v n m a+{-# INLINE transpose #-}+transpose (Matrix v) =+    let n = natValInt (Proxy :: Proxy n)+     in fromRows $ generate (\j -> generate (\i -> unsafeIndex v $ finiteInt j + finiteInt i*n) :: Vector v m a)++-- | Pure implementation of the dot product.+dotProduct :: (G.Vector v x, Num x) => Vector v n x -> Vector v n x -> x+{-# INLINE dotProduct #-}+dotProduct v1 v2 = weakDotProduct (fromSized v1) (fromSized v2)++-- | Pure implementation of the outer product.+outerProduct+    :: ( KnownNat m, KnownNat n, Num x+       , G.Vector v Int, G.Vector v x, G.Vector v (Vector v n x), G.Vector v (Vector v m x), G.Vector v (Vector v 1 x) )+     => Vector v n x -> Vector v m x -> Matrix v n m x+{-# INLINE outerProduct #-}+outerProduct v1 v2 = matrixMatrixMultiply (columnVector v1) (rowVector v2)++-- | Pure implementation of the dot product on standard vectors.+weakDotProduct :: (G.Vector v x, Num x) => v x -> v x -> x+{-# INLINE weakDotProduct #-}+weakDotProduct v1 v2 = G.foldl foldFun 0 (G.enumFromN 0 (G.length v1) :: S.Vector Int)+    where foldFun d i = d + G.unsafeIndex v1 i * G.unsafeIndex v2 i++-- | Pure 'Matrix' x 'Vector' multiplication.+matrixVectorMultiply+    :: (KnownNat m, KnownNat n, G.Vector v x, G.Vector v (Vector v n x), Num x)+    => Matrix v m n x+    -> Vector v n x+    -> Vector v m x+{-# INLINE matrixVectorMultiply #-}+matrixVectorMultiply mtx v =+    map (dotProduct v) $ toRows mtx++-- | Pure 'Matrix' x 'Matrix' multiplication.+matrixMatrixMultiply+    :: ( KnownNat m, KnownNat n, KnownNat o, Num x+       , G.Vector v Int, G.Vector v x, G.Vector v (Vector v m x), G.Vector v (Vector v n x), G.Vector v (Vector v o x) )+    => Matrix v m n x+    -> Matrix v n o x+    -> Matrix v m o x+{-# INLINE matrixMatrixMultiply #-}+matrixMatrixMultiply mtx1 mtx2 =+    fromColumns . map (matrixVectorMultiply mtx1) $ toColumns mtx2+++--- Numeric Classes ---+++--instance (Storable x, Numeric x, KnownNat n, KnownNat m)+--  => Num (Matrix S.Vector n m x) where+--    {-# INLINE (+) #-}+--    (+) (Matrix (Vector v1)) (Matrix (Vector v2)) = Matrix $ Vector (H.add v1 v2)+--    {-# INLINE (*) #-}+--    (*) (Matrix xs) (Matrix xs') = Matrix $ xs * xs'+--    {-# INLINE negate #-}+--    negate (Matrix (Vector v)) = Matrix $ Vector (H.scale (-1) v)+--    {-# INLINE abs #-}+--    abs (Matrix xs) = Matrix $ abs xs+--    {-# INLINE signum #-}+--    signum (Matrix xs) = Matrix $ signum xs+--    {-# INLINE fromInteger #-}+--    fromInteger x = Matrix . replicate $ fromInteger x

@@ -0,0 +1,11 @@+-- | Forbidden access to the static vector constructor.+module Goal.Core.Vector.Generic.Internal+    ( -- * Vector+      module Data.Vector.Generic.Sized.Internal+    ) where+++--- Imports ---+++import Data.Vector.Generic.Sized.Internal

@@ -0,0 +1,11 @@+-- | Mutable generic static vectors.+module Goal.Core.Vector.Generic.Mutable+    ( -- * Vector+      module Data.Vector.Generic.Mutable.Sized+    ) where+++--- Imports ---+++import Data.Vector.Generic.Mutable.Sized

@@ -0,0 +1,741 @@+ {-# OPTIONS_GHC -fplugin=GHC.TypeLits.KnownNat.Solver -fplugin=GHC.TypeLits.Normalise -fconstraint-solver-iterations=10 #-}+-- | Vectors and Matrices with statically typed dimensions based on storable vectors and using HMatrix where possible.+module Goal.Core.Vector.Storable+    ( -- * Vector+      module Data.Vector.Storable.Sized+      -- ** Construction+    , doubleton+    , range+      -- ** Deconstruction+    , concat+    , breakEvery+    , toPair+    -- ** Computation+    , average+    , zipFold+    -- * Matrix+    , Matrix+    , nRows+    , nColumns+    -- ** Construction+    , fromRows+    , fromColumns+    , matrixIdentity+    , outerProduct+    , sumOuterProduct+    , averageOuterProduct+    , weightedAverageOuterProduct+    , diagonalMatrix+    , fromLowerTriangular+    -- ** Deconstruction+    , toRows+    , toColumns+    , lowerTriangular+    , takeDiagonal+    -- ** Manipulation+    , columnVector+    , rowVector+    , combineTriangles+    , diagonalConcat+    , horizontalConcat+    , verticalConcat+    -- ** Computation+    , trace+    , withMatrix+    -- *** BLAS+    , scale+    , add+    , dotProduct+    , dotMap+    , matrixVectorMultiply+    , matrixMatrixMultiply+    , matrixMap+    , eigens+    , isSemiPositiveDefinite+    , determinant+    , inverseLogDeterminant+    , inverse+    , pseudoInverse+    , matrixRoot+    , transpose+    -- *** Least Squares+    , linearLeastSquares+    , meanSquaredError+    , rSquared+    , l2Norm+    , unsafeCholesky+    -- *** Convolutions+    , crossCorrelate2d+    , convolve2d+    , kernelOuterProduct+    , kernelTranspose+    -- ** Miscellaneous+    , prettyPrintMatrix+    ) where+++--- Imports ---+++-- Goal --++import Goal.Core.Util hiding (average,breakEvery,range)++-- Unqualified --++import Data.Proxy+import Data.Complex+import Foreign.Storable+import Data.Vector.Storable.Sized+import Numeric.LinearAlgebra (Field,Numeric)+import GHC.TypeNats+import Prelude hiding (concat,foldr1,concatMap,replicate,(++),length,map,sum,zip,and)++-- Qualified --++import qualified Prelude+import qualified Data.Vector.Storable as S+import qualified Goal.Core.Vector.Generic as G+import qualified Data.Vector.Generic.Sized.Internal as G+import qualified Numeric.LinearAlgebra as H+import qualified Data.List as L+++--- Generic ---+++-- | Matrices with static dimensions (storable).+type Matrix = G.Matrix S.Vector++-- | A fold over pairs of elements of 'Vector's of equal length.+zipFold :: (KnownNat n, Storable x, Storable y)+        => (z -> x -> y -> z)+        -> z+        -> Vector n x+        -> Vector n y+        -> z+{-# INLINE zipFold #-}+zipFold f z0 xs ys =+    let n = length xs+        foldfun z i = f z (unsafeIndex xs i) (unsafeIndex ys i)+     in L.foldl' foldfun z0 [0..n-1]++-- | Concatenates a vector of vectors into a single vector.+concat :: (KnownNat n, Storable x) => Vector m (Vector n x) -> Vector (m*n) x+{-# INLINE concat #-}+concat = G.concat++-- | Collect two values into a length 2 'Vector'.+doubleton :: Storable x => x -> x -> Vector 2 x+{-# INLINE doubleton #-}+doubleton = G.doubleton++-- | The number of rows in the 'Matrix'.+nRows :: forall m n a . KnownNat m => Matrix m n a -> Int+{-# INLINE nRows #-}+nRows = G.nRows++-- | The columns of columns in the 'Matrix'.+nColumns :: forall m n a . KnownNat n => Matrix m n a -> Int+{-# INLINE nColumns #-}+nColumns = G.nColumns++-- | Convert a 'Matrix' into a 'Vector' of 'Vector's of rows.+toRows :: (KnownNat m, KnownNat n, Storable x) => Matrix m n x -> Vector m (Vector n x)+{-# INLINE toRows #-}+toRows = G.toRows++-- | Turn a 'Vector' into a single column 'Matrix'.+columnVector :: Vector n a -> Matrix n 1 a+{-# INLINE columnVector #-}+columnVector = G.columnVector++-- | Turn a 'Vector' into a single row 'Matrix'.+rowVector :: Vector n a -> Matrix 1 n a+{-# INLINE rowVector #-}+rowVector = G.rowVector++-- | Create a 'Matrix' from a 'Vector' of 'Vector's which represent the rows.+fromRows :: (KnownNat n, Storable x) => Vector m (Vector n x) -> Matrix m n x+{-# INLINE fromRows #-}+fromRows = G.fromRows++-- | Uniform partition of an interval into a 'Vector'.+range :: (KnownNat n, Fractional x, Storable x) => x -> x -> Vector n x+{-# INLINE range #-}+range = G.range++-- | Reshapes a length 2 'Vector' into a pair of values.+toPair :: Storable x => Vector 2 x -> (x,x)+{-# INLINE toPair #-}+toPair = G.toPair+++--- HMatrix ---+++-- | Converts a pure, Storable-based 'Matrix' into an HMatrix matrix.+toHMatrix+    :: forall m n x . (KnownNat n, KnownNat m, H.Element x, Storable x)+    => Matrix m n x+    -> H.Matrix x+{-# INLINE toHMatrix #-}+toHMatrix (G.Matrix mtx) =+    let n = natValInt (Proxy :: Proxy n)+        m = natValInt (Proxy :: Proxy m)+     in if n == 0+           then H.fromRows $ Prelude.replicate m S.empty+           else H.reshape n $ fromSized mtx++-- | Converts an HMatrix matrix into a pure, Storable-based 'Matrix'.+fromHMatrix :: Numeric x => H.Matrix x -> Matrix m n x+{-# INLINE fromHMatrix #-}+fromHMatrix = G.Matrix . G.Vector . H.flatten++-- | Convert a 'Matrix' into a 'Vector' of 'Vector's of columns.+toColumns :: (KnownNat m, KnownNat n, Numeric x) => Matrix m n x -> Vector n (Vector m x)+{-# INLINE toColumns #-}+toColumns = toRows . transpose++-- | Create a 'Matrix' from a 'Vector' of 'Vector's which represent the columns.+fromColumns :: (KnownNat m, KnownNat n, Numeric x) => Vector n (Vector m x) -> Matrix m n x+{-# INLINE fromColumns #-}+fromColumns = transpose . fromRows++-- | Breaks a 'Vector' into a Vector of Vectors.+breakEvery :: forall n k a . (KnownNat n, KnownNat k, Storable a) => Vector (n*k) a -> Vector n (Vector k a)+{-# INLINE breakEvery #-}+breakEvery v0 =+    let k = natValInt (Proxy :: Proxy k)+        v = fromSized v0+     in generate (\i -> G.Vector $ S.unsafeSlice (finiteInt i*k) k v)+++--- BLAS ---+++-- | The sum of two 'Vector's.+add :: Numeric x => Vector n x -> Vector n x -> Vector n x+{-# INLINE add #-}+add (G.Vector v1) (G.Vector v2) = G.Vector (H.add v1 v2)++-- | Scalar multiplication of a 'Vector'.+scale :: Numeric x => x -> Vector n x -> Vector n x+{-# INLINE scale #-}+scale x (G.Vector v) = G.Vector (H.scale x v)++-- | Apply a 'Vector' operation to a 'Matrix'.+withMatrix :: (Vector (n*m) x -> Vector (n*m) x) -> Matrix n m x -> Matrix n m x+{-# INLINE withMatrix #-}+withMatrix f (G.Matrix v) = G.Matrix $ f v++-- | Returns the lower triangular part of a square matrix.+lowerTriangular :: forall n x . (Storable x, H.Element x, KnownNat n) => Matrix n n x -> Vector (Triangular n) x+{-# INLINE lowerTriangular #-}+lowerTriangular mtx =+    let hmtx = toHMatrix mtx+        rws = H.toRows hmtx+        rws' = Prelude.zipWith S.take [1..] rws+     in G.Vector $ S.concat rws'+--    let n = natValInt (Proxy :: Proxy n)+--        idxs = G.Vector . S.fromList+--            $ Prelude.concat [ from2Index n <$> Prelude.zip (repeat k) [0..k] | k <- [0..n-1] ]+--     in backpermute xs idxs++-- | Constructs a `Matrix` from a lower triangular part.+fromLowerTriangular :: forall n x . (Storable x, KnownNat n) => Vector (Triangular n) x -> Matrix n n x+{-# INLINE fromLowerTriangular #-}+fromLowerTriangular xs =+    let n = natValInt (Proxy :: Proxy n)+        idxs = generate (toTriangularIndex . to2Index n . finiteInt)+     in G.Matrix $ backpermute xs idxs++-- | Build a matrix with the given diagonal, lower triangular part given by the+-- first matrix, and upper triangular part given by the second matrix.+combineTriangles+    :: (KnownNat k, Storable x)+    => Vector k x -- ^ Diagonal+    -> Matrix k k x -- ^ Lower triangular source+    -> Matrix k k x -- ^ Upper triangular source+    -> Matrix k k x+{-# INLINE combineTriangles #-}+combineTriangles (G.Vector diag) crs1 crs2 =+    fromRows $ generate (generator (toRows crs1) (toRows crs2))+        where+            generator rws1 rws2 fnt =+                let (G.Vector rw1) = index rws1 fnt+                    (G.Vector rw2) = index rws2 fnt+                    i = fromIntegral fnt+                 in G.Vector $ S.take i rw1 S.++ S.cons (diag S.! i) (S.drop (i+1) rw2)++-- | The average of a 'Vector' of elements.+average :: (Numeric x, Fractional x) => Vector n x -> x+{-# INLINE average #-}+average (G.Vector v) = H.sumElements v / fromIntegral (S.length v)++-- | The dot product of two numerical 'Vector's.+dotProduct :: Numeric x => Vector n x -> Vector n x -> x+{-# INLINE dotProduct #-}+dotProduct v1 v2 = H.dot (fromSized v1) (fromSized v2)++-- | The determinant of a 'Matrix'.+diagonalMatrix :: forall n x . (KnownNat n, Field x) => Vector n x -> Matrix n n x+{-# INLINE diagonalMatrix #-}+diagonalMatrix v =+    let n = natValInt (Proxy :: Proxy n)+     in fromHMatrix $ H.diagRect 0 (fromSized v) n n++-- | The determinant of a 'Matrix'.+takeDiagonal :: (KnownNat n, Field x) => Matrix n n x -> Vector n x+{-# INLINE takeDiagonal #-}+takeDiagonal = G.Vector . H.takeDiag . toHMatrix++-- | The determinant of a 'Matrix'.+trace :: (KnownNat n, Field x) => Matrix n n x -> x+{-# INLINE trace #-}+trace = S.sum . H.takeDiag . toHMatrix++-- | Returns the eigenvalues and eigenvectors 'Matrix'.+eigens :: (KnownNat n, Field x) => Matrix n n x -> (Vector n (Complex Double), Vector n (Vector n (Complex Double)))+{-# INLINE eigens #-}+eigens mtx =+    let (exs,evs) = H.eig $ toHMatrix mtx+     in (G.Vector exs, G.Vector . S.fromList $ G.Vector <$> H.toColumns evs)++-- | Test if the matrix is semi-positive definite.+isSemiPositiveDefinite :: (KnownNat n, Field x) => Matrix n n x -> Bool+{-# INLINE isSemiPositiveDefinite #-}+isSemiPositiveDefinite =+    and . map ((0 <=) . realPart) . fst . eigens++-- | Returns the inverse, the logarithm of the absolute value of the+-- determinant, and the sign of the determinant of a given matrix.+inverseLogDeterminant :: (KnownNat n, Field x) => Matrix n n x -> (Matrix n n x, x, x)+{-# INLINE inverseLogDeterminant #-}+inverseLogDeterminant mtx =+    let (imtx,(ldet,sgn)) = H.invlndet $ toHMatrix mtx+     in (fromHMatrix imtx, ldet, sgn)++-- | The determinant of a 'Matrix'.+determinant :: (KnownNat n, Field x) => Matrix n n x -> x+{-# INLINE determinant #-}+determinant = H.det . toHMatrix++-- | Transpose a 'Matrix'.+transpose+    :: forall m n x . (KnownNat m, KnownNat n, Numeric x)+    => Matrix m n x+    -> Matrix n m x+{-# INLINE transpose #-}+transpose (G.Matrix mtx) =+    G.Matrix $ withVectorUnsafe (H.flatten . H.tr . H.reshape (natValInt (Proxy :: Proxy n))) mtx++-- | Diagonally concatenate two matrices, padding the gaps with zeroes.+diagonalConcat+    :: (KnownNat n, KnownNat m, KnownNat o, KnownNat p, Numeric x)+    => Matrix n m x -> Matrix o p x -> Matrix (n+o) (m+p) x+{-# INLINE diagonalConcat #-}+diagonalConcat mtx10 mtx20 =+    let mtx1 = toHMatrix mtx10+        mtx2 = toHMatrix mtx20+     in fromHMatrix $ H.diagBlock [mtx1,mtx2]++-- | Diagonally concatenate two matrices, padding the gaps with zeroes.+horizontalConcat+    :: (KnownNat n, KnownNat m, KnownNat o, Numeric x)+    => Matrix n m x -> Matrix n o x -> Matrix n (m+o) x+{-# INLINE horizontalConcat #-}+horizontalConcat mtx10 mtx20 =+    let mtx1 = toHMatrix mtx10+        mtx2 = toHMatrix mtx20+     in fromHMatrix $ mtx1 H.||| mtx2++-- | Diagonally concatenate two matrices, padding the gaps with zeroes.+verticalConcat+    :: (KnownNat n, KnownNat m, KnownNat o, Numeric x)+    => Matrix n o x -> Matrix m o x -> Matrix (n+m) o x+{-# INLINE verticalConcat #-}+verticalConcat mtx10 mtx20 =+    let mtx1 = toHMatrix mtx10+        mtx2 = toHMatrix mtx20+     in fromHMatrix $ mtx1 H.=== mtx2++-- | Invert a 'Matrix'.+inverse :: forall n x . (KnownNat n, Field x) => Matrix n n x -> Matrix n n x+{-# INLINE inverse #-}+inverse (G.Matrix mtx) =+    G.Matrix $ withVectorUnsafe (H.flatten . H.inv . H.reshape (natValInt (Proxy :: Proxy n))) mtx++-- | Pseudo-Invert a 'Matrix'.+pseudoInverse :: forall n x . (KnownNat n, Field x) => Matrix n n x -> Matrix n n x+{-# INLINE pseudoInverse #-}+pseudoInverse (G.Matrix mtx) =+    G.Matrix $ withVectorUnsafe (H.flatten . H.pinv . H.reshape (natValInt (Proxy :: Proxy n))) mtx++-- | Square root of a 'Matrix'.+matrixRoot :: forall n x . (KnownNat n, Field x) => Matrix n n x -> Matrix n n x+{-# INLINE matrixRoot #-}+matrixRoot (G.Matrix mtx) =+    G.Matrix $ withVectorUnsafe (H.flatten . H.sqrtm . H.reshape (natValInt (Proxy :: Proxy n))) mtx++-- | The outer product of two 'Vector's.+outerProduct :: (KnownNat m, KnownNat n, Numeric x) => Vector m x -> Vector n x -> Matrix m n x+{-# INLINE outerProduct #-}+outerProduct v1 v2 =+    fromHMatrix $ H.outer (fromSized v1) (fromSized v2)++-- | The summed outer product of two lists of 'Vector's.+sumOuterProduct :: (KnownNat m, KnownNat n, Fractional x, Numeric x) => [(Vector m x,Vector n x)] -> Matrix m n x+{-# INLINE sumOuterProduct #-}+sumOuterProduct v12s =+    let (v1s,v2s) = L.unzip v12s+        mtx1 = H.fromColumns $ fromSized <$> v1s+        mtx2 = H.fromRows $ fromSized <$> v2s+     in fromHMatrix (mtx1 H.<> mtx2)++-- | The average outer product of two lists of 'Vector's.+averageOuterProduct :: (KnownNat m, KnownNat n, Fractional x, Numeric x) => [(Vector m x,Vector n x)] -> Matrix m n x+{-# INLINE averageOuterProduct #-}+averageOuterProduct v12s =+    let (v1s,v2s) = L.unzip v12s+        mtx1 = H.fromColumns $ fromSized <$> v1s+        (_,n) = H.size mtx1+        mtx2 = H.scale (1/fromIntegral n) . H.fromRows $ fromSized <$> v2s+     in fromHMatrix (mtx1 H.<> mtx2)++-- | The average outer product of two lists of 'Vector's.+weightedAverageOuterProduct+    :: ( KnownNat m, KnownNat n, Fractional x, Numeric x )+    => [(x,Vector m x,Vector n x)]+    -> Matrix m n x+{-# INLINE weightedAverageOuterProduct #-}+weightedAverageOuterProduct wv12s =+    let (ws,v1s,v2s) = L.unzip3 wv12s+        v1s' = L.zipWith H.scale ws $ fromSized <$> v1s+        mtx1 = H.fromColumns v1s'+        mtx2 = H.fromRows $ fromSized <$> v2s+     in fromHMatrix (mtx1 H.<> mtx2)++-- | The identity 'Matrix'.+matrixIdentity :: forall n x . (KnownNat n, Numeric x, Num x) => Matrix n n x+{-# INLINE matrixIdentity #-}+matrixIdentity =+    fromHMatrix . H.ident $ natValInt (Proxy :: Proxy n)++-- | The dot products of one vector with a list of vectors.+dotMap :: (KnownNat n, Numeric x) => Vector n x -> [Vector n x] -> [x]+{-# INLINE dotMap #-}+dotMap v vs =+    let mtx' = H.fromRows $ fromSized <$> vs+     in H.toList $ mtx' H.#> fromSized v+--     in if S.null w+--           then replicate 0+--           else fmap G.Vector . H.toColumns $ toHMatrix mtx H.<> mtx'++-- | Map a linear transformation over a list of 'Vector's.+matrixMap :: (KnownNat m, KnownNat n, Numeric x)+                     => Matrix m n x -> [Vector n x] -> [Vector m x]+{-# INLINE matrixMap #-}+matrixMap mtx vs =+    let mtx' = H.fromColumns $ fromSized <$> vs+     in fmap G.Vector . H.toColumns $ toHMatrix mtx H.<> mtx'+--     in if S.null w+--           then replicate 0+--           else fmap G.Vector . H.toColumns $ toHMatrix mtx H.<> mtx'+++-- | Apply a linear transformation to a 'Vector'.+matrixVectorMultiply :: (KnownNat m, KnownNat n, Numeric x)+                     => Matrix m n x -> Vector n x -> Vector m x+{-# INLINE matrixVectorMultiply #-}+matrixVectorMultiply mtx v =+    G.Vector $ toHMatrix mtx H.#> fromSized v+--    let w = toHMatrix mtx H.#> fromSized v+--     in if S.null w+--           then replicate 0+--           else G.Vector w++-- | Multiply a 'Matrix' with a second 'Matrix'.+matrixMatrixMultiply+    :: (KnownNat m, KnownNat n, KnownNat o, Numeric x)+    => Matrix m n x+    -> Matrix n o x+    -> Matrix m o x+{-# INLINE matrixMatrixMultiply #-}+matrixMatrixMultiply mtx1 mtx2 = fromHMatrix $ toHMatrix mtx1 H.<> toHMatrix mtx2++-- | Pretty print the values of a 'Matrix' (for extremely simple values of pretty).+prettyPrintMatrix :: (KnownNat m, KnownNat n, Numeric a, Show a) => Matrix m n a -> IO ()+prettyPrintMatrix = print . toHMatrix++-- | The Mean Squared difference between two vectors.+meanSquaredError+    :: KnownNat k+    => Vector k Double+    -> Vector k Double+    -> Double+{-# INLINE meanSquaredError #-}+meanSquaredError ys yhts = average $ map square (ys - yhts)++-- | L2 length of a vector.+l2Norm+    :: KnownNat k+    => Vector k Double+    -> Double+{-# INLINE l2Norm #-}+l2Norm (G.Vector xs) = H.norm_2 xs++-- | Computes the coefficient of determintation for the given outputs and model+-- predictions.+rSquared+    :: KnownNat k+    => Vector k Double -- ^ Dependent variable observations+    -> Vector k Double -- ^ Predicted Values+    -> Double -- ^ R-squared+{-# INLINE rSquared #-}+rSquared ys yhts =+    let ybr = average ys+        ssres = sum $ map square (ys - yhts)+        sstot = sum $ map (square . subtract ybr) ys+     in 1 - (ssres/sstot)++-- | Solves the linear least squares problem.+linearLeastSquares+    :: KnownNat l+    => [Vector l Double] -- ^ Independent variable observations+    -> [Double] -- ^ Dependent variable observations+    -> Vector l Double -- ^ Parameter estimates+{-# INLINE linearLeastSquares #-}+linearLeastSquares as xs =+    G.Vector $ H.fromRows (fromSized <$> as) H.<\> S.fromList xs+++unsafeCholesky+    :: (KnownNat n, Field x, Storable x)+    => Matrix n n x+    -> Matrix n n x+unsafeCholesky =+    transpose . fromHMatrix . H.chol . H.trustSym . toHMatrix+++--- Convolutions ---+++-- | 2d cross-correlation of a kernel over a matrix of values.+crossCorrelate2d+    :: forall nk rdkr rdkc mr mc md x+    . ( KnownNat rdkr, KnownNat rdkc, KnownNat md, KnownNat mr, KnownNat mc+      , KnownNat nk, Numeric x, Storable x )+      => Proxy rdkr -- ^ Number of Kernel rows+      -> Proxy rdkc -- ^ Number of Kernel columns+      -> Proxy mr -- ^ Number of Matrix/Image rows+      -> Proxy mc -- ^ Number of Kernel/Image columns+      -> Matrix nk (md*(2*rdkr+1)*(2*rdkc+1)) x -- ^ Kernels (nk is their number)+      -> Matrix md (mr*mc) x -- ^ Image (md is the depth)+      -> Matrix nk (mr*mc) x -- ^ Cross-correlated image+{-# INLINE crossCorrelate2d #-}+crossCorrelate2d prdkr prdkc pmr pmc krns (G.Matrix v) =+    let pmd = Proxy :: Proxy md+        mtx = im2col prdkr prdkc pmd pmr pmc v+     in matrixMatrixMultiply krns mtx++-- | The transpose of a convolutional kernel.+kernelTranspose+    :: (KnownNat nk, KnownNat md, KnownNat rdkr, KnownNat rdkc, Numeric x, Storable x)+    => Proxy nk+    -> Proxy md+    -> Proxy rdkr+    -> Proxy rdkc+    -> Matrix nk (md*(2*rdkr+1)*(2*rdkc+1)) x -- ^ Kernels (nk is their number)+    -> Matrix md (nk*(2*rdkr+1)*(2*rdkc+1)) x -- ^ Kernels (nk is their number)+{-# INLINE kernelTranspose #-}+kernelTranspose pnk pmd prdkr prdkc (G.Matrix kv) = G.Matrix . backpermute kv $ kernelTransposeIndices pnk pmd prdkr prdkc++-- | 2d convolution of a kernel over a matrix of values. This is the adjoint of crossCorrelate2d.+convolve2d+    :: forall nk rdkr rdkc md mr mc x+    . ( KnownNat rdkr, KnownNat rdkc, KnownNat mr, KnownNat mc+      , KnownNat md, KnownNat nk, Numeric x, Storable x )+      => Proxy rdkr -- ^ Number of Kernel rows+      -> Proxy rdkc -- ^ Number of Kernel columns+      -> Proxy mr -- ^ Number of Matrix/Image rows+      -> Proxy mc -- ^ Number of Kernel/Image columns+      -> Matrix nk (md*(2*rdkr+1)*(2*rdkc+1)) x -- ^ Kernels (nk is their number)+      -> Matrix nk (mr*mc) x -- ^ Dual image (nk is its depth)+      -> Matrix md (mr*mc) x -- ^ Convolved image+{-# INLINE convolve2d #-}+convolve2d prdkr prdkc pmr pmc krn mtxs =+    let pnk = Proxy :: Proxy nk+        pmd = Proxy :: Proxy md+        krn' = kernelTranspose pnk pmd prdkr prdkc krn+     in crossCorrelate2d prdkr prdkc pmr pmc krn' mtxs++-- | The outer product of an image and a dual image to produce a convolutional kernel.+kernelOuterProduct+    :: forall nk rdkr rdkc md mr mc x+    . ( KnownNat rdkr, KnownNat rdkc, KnownNat mr, KnownNat mc+      , KnownNat md, KnownNat nk, Numeric x, Storable x )+      => Proxy rdkr -- ^ Number of Kernel rows+      -> Proxy rdkc -- ^ Number of Kernel columns+      -> Proxy mr -- ^ Number of Matrix/Image rows+      -> Proxy mc -- ^ Number of Kernel/Image columns+      -> Matrix nk (mr*mc) x -- ^ Dual image (nk is its depth)+      -> Matrix md (mr*mc) x -- ^ Image (md is the depth)+      -> Matrix nk (md*(2*rdkr+1)*(2*rdkc+1)) x -- ^ Kernels+{-# INLINE kernelOuterProduct #-}+kernelOuterProduct prdkr prdkc pmr pmc omtx (G.Matrix v) =+    let pmd = Proxy :: Proxy md+        imtx = im2col prdkr prdkc pmd pmr pmc v+     in matrixMatrixMultiply omtx $ transpose imtx+++--- Internal ---+++toTriangularIndex :: (Int,Int) -> Int+toTriangularIndex (i,j)+    | i >= j = triangularNumber i + j+    | otherwise = toTriangularIndex (j,i)++to2Index :: Int -> Int -> (Int,Int)+to2Index nj ij = divMod ij nj++to3Index :: Int -> Int -> Int -> (Int,Int,Int)+{-# INLINE to3Index #-}+to3Index nj nk ijk =+    let nj' = nj*nk+        (i,jk) = divMod ijk nj'+        (j,k) = divMod jk nk+     in (i,j,k)++from3Index :: Int -> Int -> (Int,Int,Int) -> Int+{-# INLINE from3Index #-}+from3Index nj nk (i,j,k) =+    let nj' = nj*nk+     in i*nj' + j*nk + k++to4Index :: Int -> Int -> Int -> Int -> (Int,Int,Int,Int)+{-# INLINE to4Index #-}+to4Index nj nk nl ijkl =+    let nk' = nl*nk+        nj' = nj*nk'+        (i,jkl) = divMod ijkl nj'+        (j,kl) = divMod jkl nk'+        (k,l) = divMod kl nl+     in (i,j,k,l)++from4Index :: Int -> Int -> Int -> (Int,Int,Int,Int) -> Int+{-# INLINE from4Index #-}+from4Index nj nk nl (i,j,k,l) =+    let nk' = nl*nk+        nj' = nj*nk'+     in i*nj' + j*nk' + k*nl + l++kernelTransposeIndices+    :: (KnownNat nk, KnownNat md, KnownNat rdkr, KnownNat rdkc)+    => Proxy nk+    -> Proxy md+    -> Proxy rdkr+    -> Proxy rdkc+    -> Vector (nk*md*(2*rdkr+1)*(2*rdkc+1)) Int+{-# INLINE kernelTransposeIndices #-}+kernelTransposeIndices pnk pmd prdkr prdkc =+    let nkrn = natValInt pnk+        md = natValInt pmd+        rdkr = natValInt prdkr+        rdkc = natValInt prdkc+        dmkr = 2*rdkr+1+        dmkc = 2*rdkc+1+        nl = dmkc+        nk = dmkr+        nj = nkrn+        nl' = dmkc+        nk' = dmkr+        nj' = md+        reIndex idx =+            let (i,j,k,l) = to4Index nj nk nl idx+             in from4Index nj' nk' nl' (j,i,nk-1-k,nl-1-l)+     in generate (reIndex . fromIntegral)++im2colIndices+    :: forall rdkr rdkc mr mc md+     . (KnownNat rdkr, KnownNat rdkc, KnownNat mr, KnownNat mc, KnownNat md)+    => Proxy rdkr+    -> Proxy rdkc+    -> Proxy md+    -> Proxy mr+    -> Proxy mc+    -> Vector (((2*rdkr+1)*(2*rdkc+1)*md)*(mr*mc)) Int+{-# INLINE im2colIndices #-}+im2colIndices prdkr prdkc _ pmr pmc =+    let rdkr = natValInt prdkr+        rdkc = natValInt prdkc+        nj = (2*rdkr + 1)+        nk = (2*rdkc + 1)+        reWindow idx =+            let (i,j,k) = to3Index nj nk idx+             in windowIndices prdkr prdkc pmr pmc i j k+          in (concatMap reWindow :: Vector ((2*rdkr+1)*(2*rdkc+1)*md) Int -> Vector (((2*rdkr+1)*(2*rdkc+1)*md)*(mr*mc)) Int) $ generate finiteInt++im2col+    :: forall rdkr rdkc md mr mc x+    . (KnownNat rdkr, KnownNat rdkc, KnownNat mc, KnownNat md, KnownNat mr, Num x, Storable x)+    => Proxy rdkr+    -> Proxy rdkc+    -> Proxy md+    -> Proxy mr+    -> Proxy mc+    -> Vector (md*mr*mc) x+    -> Matrix (md*(2*rdkr+1)*(2*rdkc+1)) (mr*mc) x+{-# INLINE im2col #-}+im2col prdkr prdkc pmd pmr pmc mtx =+    let idxs = im2colIndices prdkr prdkc pmd pmr pmc+        mtx' = padMatrix prdkr prdkc pmd pmr pmc mtx+     in G.Matrix $ backpermute mtx' idxs++windowIndices+    :: forall rdkr rdkc mr mc . (KnownNat rdkr, KnownNat rdkc, KnownNat mr, KnownNat mc)+    => Proxy rdkr+    -> Proxy rdkc+    -> Proxy mr+    -> Proxy mc+    -> Int+    -> Int+    -> Int+    -> Vector (mr*mc) Int+{-# INLINE windowIndices #-}+windowIndices prdkr prdkc pmr pmc kd kr kc =+    let rdkr = natValInt prdkr+        rdkc = natValInt prdkc+        mr = natValInt pmr+        mc = natValInt pmc+        mrc = mr*mc+        nj' = mr + 2*rdkr+        nk' = mc + 2*rdkc+        reIndex idx =+            let (j,k) = divMod idx mc+             in from3Index nj' nk' (kd,j+kr,k+kc)+     in G.Vector $ S.generate mrc reIndex++padMatrix+    :: forall rdkr rdkc mr mc md x+    . (KnownNat rdkr, KnownNat rdkc, KnownNat md, KnownNat mr, KnownNat mc, Num x, Storable x)+    => Proxy rdkr+    -> Proxy rdkc+    -> Proxy md+    -> Proxy mr+    -> Proxy mc+    -> Vector (md*mr*mc) x+    -> Vector (md*(mr + 2*rdkr)*(mc + 2*rdkc)) x+{-# INLINE padMatrix #-}+padMatrix _ _ _ _ _ v =+    let mtxs :: Vector md (Matrix mr mc x)+        mtxs = map G.Matrix $ breakEvery v+        pdrs :: Vector rdkr (Vector mc x)+        pdrs = replicate $ replicate 0+        mtxs' = map (\mtx -> fromRows $ pdrs ++ toRows mtx ++ pdrs) mtxs+        pdcs :: Vector rdkc (Vector (mr + 2*rdkr) x)+        pdcs = replicate $ replicate 0+     in concatMap G.toVector $ map (\mtx' -> G.fromColumns $ pdcs ++ G.toColumns mtx' ++ pdcs) mtxs'++

@@ -0,0 +1,12 @@+This is the least interesting package in the Goal libraries, and serves simply+to re-export existing libraries and provide essential utility factors in a+manner that is compatible with Goal. Nevertheless, there are a few modules worth+highlighting.++**Goal.Core.Circuit**: Provides an implementation of monadic Mealy automata to+facilitate simple stream-based processing.++**Goal.Core.Vector.Storable**: Combines the+[vector-sized](https://hackage.haskell.org/package/vector-sized) and+[hmatrix](https://hackage.haskell.org/package/hmatrix) libraries to provide+efficient linear algebra with static sizes.

@@ -0,0 +1,110 @@+{-# LANGUAGE+    TypeOperators,+    NoStarIsType,+    DataKinds+    #-}++import Goal.Core+import qualified Goal.Core.Vector.Generic as G+import qualified Goal.Core.Vector.Storable as S+import qualified Goal.Core.Vector.Boxed as B++import qualified Numeric.LinearAlgebra as H+import qualified Criterion.Main as C+import qualified System.Random.MWC.Probability as P++++--- Globals ---+++-- Sizes --++type KRadius = 2+type KDiameter = (2*KRadius + 1)+type KNumber = 50+type MSize = 50+type MDepth = 50++pkr :: Proxy KRadius+pkr = Proxy++pkd :: Proxy KDiameter+pkd = Proxy++pms :: Proxy MSize+pms = Proxy++kdmt,ms,krd :: Int+krd = natValInt pkr+kdmt = natValInt pkd+ms = natValInt pms++goalCorr+    :: ( S.Matrix KNumber (KDiameter*KDiameter*MDepth) Double+       , S.Matrix MDepth (MSize*MSize) Double )+    -> S.Matrix KNumber (MSize*MSize) Double+goalCorr (krns,mtx) = S.crossCorrelate2d pkr pkr pms pms krns mtx++goalConv+    :: ( S.Matrix KNumber (KDiameter*KDiameter*MDepth) Double+       , S.Matrix KNumber (MSize*MSize) Double )+    -> S.Matrix MDepth (MSize*MSize) Double+goalConv (krns,mtx) = S.convolve2d pkr pkr pms pms krns mtx++hmatrixCorr+    :: (B.Vector KNumber (B.Vector MDepth (H.Matrix Double)), B.Vector MDepth (H.Matrix Double))+    -> B.Vector KNumber (H.Matrix Double)+hmatrixCorr (krnss,mtxs) = fromJust . B.fromList+    $ [ foldr1 H.add [ H.corr2 krn mtx | (krn,mtx) <- B.toList $ B.zip krns mtxs ] | krns <- B.toList krnss ]++repadHMatrix :: H.Matrix Double -> H.Matrix Double+repadHMatrix hmtx =+    let rw' = replicate krd . H.fromList $ replicate ms 0+        hmtx' = H.fromRows . (rw' ++) . (++ rw') $ H.toRows hmtx+        cl' = replicate krd . H.fromList $ replicate (ms + 2*krd) 0+     in H.fromColumns . (cl' ++) . (++ cl') $ H.toColumns hmtx'++--depadHMatrix :: H.Matrix Double -> H.Matrix Double+--depadHMatrix hmtx =+--    let hmtx' = H.fromRows . take mn . drop kdmt $ H.toRows hmtx+--     in H.fromColumns . take mn . drop kdmt $ H.toColumns hmtx'+++-- Benchmark++main :: IO ()+main = do++    let rnd :: P.Prob IO Double+        rnd = P.uniformR (-1,1)++    mtxv <- P.withSystemRandom . P.sample $ S.replicateM rnd+    krnv <- P.withSystemRandom . P.sample $ S.replicateM rnd+    mtxz <- P.withSystemRandom . P.sample $ S.replicateM rnd++    let krn = G.Matrix krnv+        mtx = G.Matrix mtxv+        crr = goalCorr (krn,mtx)++    let hmtxs = H.reshape ms . G.fromSized <$> G.convert (S.toRows mtx)+        hkrnss = fmap (H.reshape kdmt . G.fromSized . G.convert) . B.breakEvery . G.convert <$> G.convert (S.toRows krn)+        hcrr = hmatrixCorr (hkrnss,repadHMatrix <$> hmtxs)++    putStrLn "Squared Error between HMatrix and Goal solutions:"+    print . sum $+        B.zipWith (\mtx1 mtx2 -> H.sumElements . H.cmap (^(2 :: Int)) . H.add mtx1 $ H.scale (-1) mtx2)+         hcrr $ H.reshape ms . G.fromSized <$> G.convert (S.toRows crr)+    putStrLn ""++    putStrLn "Transpose Dot Product Difference:"+    let mtx' = G.Matrix mtxz+        cnv = goalConv (krn,mtx')+    print $ S.dotProduct (G.toVector crr) (G.toVector mtx')+        - S.dotProduct (G.toVector cnv) (G.toVector mtx)+    putStrLn ""++    C.defaultMain+       [ C.bench "goal-corr" $ C.nf goalCorr (krn,mtx)+       , C.bench "goal-conv" $ C.nf goalConv (krn,mtx')+       , C.bench "hmatrix-corr" $ C.nf hmatrixCorr (hkrnss,hmtxs) ]

@@ -0,0 +1,129 @@+{-# LANGUAGE DataKinds #-}++import Goal.Core+import qualified Goal.Core.Vector.Generic as G+import qualified Goal.Core.Vector.Storable as S+import qualified Goal.Core.Vector.Boxed as B++import qualified Numeric.LinearAlgebra as H+import qualified Criterion.Main as C+import qualified System.Random.MWC.Probability as P+++--- Globals ---+++-- Sizes --++type M = 1000+type N = 10++n,m :: Int+m = 1000+n = 10++-- Matrices --++goalMatrix1 :: S.Matrix M M Double+goalMatrix1 = G.Matrix $ S.generate fromIntegral++goalMatrix2 :: S.Matrix M N Double+goalMatrix2 = G.Matrix $ S.generate fromIntegral++bGoalMatrix1 :: B.Matrix M M Double+bGoalMatrix1 = G.Matrix $ B.generate fromIntegral++bGoalMatrix2 :: B.Matrix M N Double+bGoalMatrix2 = G.Matrix $ B.generate fromIntegral++goalVal :: (S.Matrix M M Double,S.Matrix M N Double) -> Double+goalVal (m1,m2) =+    let G.Matrix v = S.matrixMatrixMultiply m1 m2+     in S.sum v++goalVal2 :: (B.Matrix M M Double, B.Matrix M N Double) -> Double+goalVal2 (m1,m2) =+    let G.Matrix v = B.matrixMatrixMultiply m1 m2+     in sum v++hmatrixMatrix1 :: H.Matrix Double+hmatrixMatrix1 = H.fromLists . take m . breakEvery m $ [0..]++hmatrixMatrix2 :: H.Matrix Double+hmatrixMatrix2 = H.fromLists . take m . breakEvery n $ [0..]++hmatrixVal :: (H.Matrix Double,H.Matrix Double) -> Double+hmatrixVal (m1,m2) =+    let m3 = m1 H.<> m2+     in H.sumElements m3++-- Benchmark+main :: IO ()+main = do++    let rnd :: P.Prob IO Double+        rnd = P.uniformR (-1,1)++    v1 <- P.withSystemRandom . P.sample $ S.replicateM rnd+    v2 <- P.withSystemRandom . P.sample $ S.replicateM rnd++    let m1 = G.Matrix v1+        m2 = G.Matrix v2++    let bm1 = G.Matrix $ G.convert v1+        bm2 = G.Matrix $ G.convert v2++    let m1'' = H.fromLists . take m . breakEvery m $!! S.toList v1+        m2'' = H.fromLists . take m . breakEvery n $!! S.toList v2++    C.defaultMain+       [ C.bench "generative-goal" $ C.nf goalVal (goalMatrix1,goalMatrix2)+       , C.bench "generative-goal2" $ C.nf goalVal2 (bGoalMatrix1,bGoalMatrix2)+       , C.bench "generative-hmatrix" $ C.nf hmatrixVal (hmatrixMatrix1,hmatrixMatrix2)+       , C.bench "random-goal" $ C.nf goalVal (m1,m2)+       , C.bench "random-goal2" $ C.nf goalVal2 (bm1,bm2)+       , C.bench "random-hmatrix" $ C.nf hmatrixVal (m1'',m2'') ]++-- Sanity Check+--sanityCheck :: IO ()+--sanityCheck = do+--+--    let rnd :: P.Prob IO Double+--        rnd = P.uniformR (-1,1)+--+--    putStrLn "Goal 1:"+--    print $ S.matrixMatrixMultiply goalMatrix1 goalMatrix2+--    putStrLn "Goal 2:"+--    print $ B.matrixMatrixMultiply bGoalMatrix1 bGoalMatrix2+--    putStrLn "Matrix:"+--    print $ M.multStd2 matrixMatrix1 matrixMatrix2+--    putStrLn "HMatrix:"+--    print $ hmatrixMatrix1 H.<> hmatrixMatrix2+--+--    v1 <- P.withSystemRandom . P.sample $ S.replicateM rnd+--    v2 <- P.withSystemRandom . P.sample $ S.replicateM rnd+--+--    let m1 = Matrix v1+--        m2 = Matrix v2+--+--    let bm1 = Matrix $ G.convert v1+--        bm2 = Matrix $ G.convert v2+--+--    let m1' = M.fromLists . take m . breakEvery m $!! S.toList v1+--        m2' = M.fromLists . take m . breakEvery n $!! S.toList v2+--+--    let m1'' = H.fromLists . take m . breakEvery m $!! S.toList v1+--        m2'' = H.fromLists . take m . breakEvery n $!! S.toList v2+--+--    putStrLn "Goal 1:"+--    print $ goalVal (m1,m2)+--    print $ S.matrixMatrixMultiply m1 m2+--    putStrLn "Goal 2:"+--    print $ goalVal2 (bm1,bm2)+--    print $ B.matrixMatrixMultiply bm1 bm2+--    putStrLn "Matrix:"+--    print $ matrixVal (m1',m2')+--    print $ M.multStd2 m1' m2'+--    putStrLn "HMatrix:"+--    print $ hmatrixVal (m1'',m2'')+--    print $ m1'' H.<> m2''

@@ -0,0 +1,71 @@+{-# LANGUAGE DataKinds,ScopedTypeVariables,FlexibleContexts #-}++import Goal.Core+import qualified Goal.Core.Vector.Generic as G+import qualified Goal.Core.Vector.Storable as S+import qualified Goal.Core.Vector.Boxed as B++import qualified Numeric.LinearAlgebra as H+import qualified Criterion.Main as C+import qualified System.Random.MWC.Probability as P+++--- Globals ---+++type Rows = 1000+type Columns = 1000++n :: Int+n = 20++addMatrices :: (KnownNat n, KnownNat k) => S.Matrix n k Double -> S.Matrix n k Double -> S.Matrix n k Double+addMatrices (G.Matrix mtx1) (G.Matrix mtx2) = G.Matrix $ S.add mtx1 mtx2+++--- Function ---+++averageOuterProduct2 v12s =+    let ln = fromIntegral $ length v12s+     in S.withMatrix (S.scale ln) $ foldr foldfun (G.Matrix $ S.replicate 0) v12s+    where foldfun (v1,v2) mtx = addMatrices mtx $ S.outerProduct v1 v2++averageWeightedOuterProduct2 wv12s =+    foldr foldfun (G.Matrix $ S.replicate 0) wv12s+        where foldfun (w,v1,v2) mtx = addMatrices mtx . S.outerProduct v1 $ S.scale w v2+++--- Main ---+++main :: IO ()+main = do++    let rnd :: P.Prob IO Double+        rnd = P.uniformR (-1,1)++    v1s :: [S.Vector Rows Double]+        <- P.withSystemRandom . P.sample . replicateM n $ S.replicateM rnd+    v2s :: [S.Vector Columns Double]+        <- P.withSystemRandom . P.sample . replicateM n $ S.replicateM rnd++    let ws :: [Double]+        ws0 = [1..fromIntegral n]+        ws = (/sum ws0) <$> ws0++    let v12s = zip v1s v2s+        wv12s = zip3 ws v1s v2s+    let v11s = zip v1s v1s+        wv11s = zip3 ws v1s v1s++    C.defaultMain+       [ C.bench "average-outer-product" $ C.nf averageOuterProduct2 v12s+       , C.bench "bulk-outer-product" $ C.nf S.averageOuterProduct v12s+       , C.bench "average-weighted-outer-product" $ C.nf averageWeightedOuterProduct2 wv12s+       , C.bench "bulk-weighted-outer-product" $ C.nf S.weightedAverageOuterProduct wv12s+       , C.bench "average-outer-product0" $ C.nf averageOuterProduct2 v11s+       , C.bench "bulk-outer-product0" $ C.nf S.averageOuterProduct v11s+       , C.bench "average-weighted-outer-product0" $ C.nf averageWeightedOuterProduct2 wv11s+       , C.bench "bulk-weighted-outer-product0" $ C.nf S.weightedAverageOuterProduct wv11s+       , C.bench "triangular-weighted-outer-product0" $ C.nf (S.lowerTriangular . S.weightedAverageOuterProduct) wv11s ]

@@ -1,43 +1,103 @@+cabal-version: 3.0 name: goal-core-version: 0.1-synopsis: Core imports for Geometric Optimization Libraries.-description: Core provides a bunch of convenience functions, basic exports, as-    well as plotting functionality, which are independent of the rest of the-    library.-license: BSD3+synopsis: Common, non-geometric tools for use with Goal+description: goal-core re-exports a number of other libraries, and provides a set of additional utility functions useful for scientific computing. In particular, implementations of Mealy Automata (Circuits), tools for working with CSV files and gnuplot, and a module which combines vector-sized vectors with hmatrix.+license: BSD-3-Clause license-file: LICENSE+extra-source-files: README.md+version: 0.20 author: Sacha Sokoloski-maintainer: sokolo@mis.mpg.de+maintainer: sacha.sokoloski@mailbox.org+homepage: https://gitlab.com/sacha-sokoloski/goal category: Math build-type: Simple-cabal-version: >=1.10  library     exposed-modules:         Goal.Core,-        Goal.Core.Plot,-        Goal.Core.Plot.Contour+        Goal.Core.Util,+        Goal.Core.Project,+        Goal.Core.Circuit,+        Goal.Core.Vector.Storable,+        Goal.Core.Vector.Generic,+        Goal.Core.Vector.Generic.Internal,+        Goal.Core.Vector.Generic.Mutable,+        Goal.Core.Vector.Boxed     build-depends:-        base==4.*,-        data-default-class==0.0.*,-        lens==4.*,-        containers==0.5.*,-        colour==2.3.*,-        gtk==0.14.*,-        cairo==0.13.*,-        Chart==1.5.*,-        Chart-cairo==1.5.*,-        Chart-gtk==1.5.*-    default-extensions: TemplateHaskell+        base >= 4.13 && < 4.15,+        directory,+        containers,+        vector,+        math-functions,+        hmatrix,+        vector-sized,+        finite-typelits,+        ghc-typelits-knownnat,+        ghc-typelits-natnormalise,+        deepseq,+        process,+        hmatrix-gsl,+        primitive,+        bytestring,+        cassava,+        async,+        criterion,+        optparse-applicative     default-language: Haskell2010-    ghc-options: -O2 -Wall -fno-warn-type-defaults -fno-warn-missing-signatures+    default-extensions:+        ScopedTypeVariables,+        ExplicitNamespaces,+        TypeOperators,+        KindSignatures,+        DataKinds,+        RankNTypes,+        TypeFamilies,+        FlexibleContexts,+        MultiParamTypeClasses,+        ConstraintKinds,+        GADTs,+        NoStarIsType,+        FlexibleInstances+    ghc-options: -Wall -O2 -executable contours-    main-is: contours.hs-    hs-source-dirs: scripts-    ghc-options: -Wall -O2 -threaded -rtsopts -fno-warn-type-defaults-        -fno-warn-missing-signatures -fno-warn-unused-do-bind+benchmark outer-products+    type: exitcode-stdio-1.0+    main-is: outer-products.hs+    hs-source-dirs: benchmarks     build-depends:-        base==4.*,-        goal-core==0.1+        base,+        hmatrix,+        mwc-random,+        mwc-probability,+        criterion,+        goal-core     default-language: Haskell2010+    ghc-options: -Wall -O2++benchmark convolutions+    type: exitcode-stdio-1.0+    main-is: convolutions.hs+    hs-source-dirs: benchmarks+    build-depends:+        base,+        hmatrix,+        mwc-random,+        mwc-probability,+        criterion,+        goal-core+    default-language: Haskell2010+    ghc-options: -Wall -O2++benchmark multiplications+    type: exitcode-stdio-1.0+    main-is: multiplications.hs+    hs-source-dirs: benchmarks+    build-depends:+        base,+        hmatrix,+        mwc-random,+        mwc-probability,+        criterion,+        goal-core+    default-language: Haskell2010+    ghc-options: -Wall -O2

@@ -1,27 +0,0 @@-import Goal.Core----- Test -----main :: IO ()-main = do-    let f x y = 0.2 * x + 0.2 * y + sin x + sin y-        sz = 7.5-        nstp = 1000-        rng = (-sz,sz,nstp)-        niso = 20-        cntrs = contours rng rng niso f-        clrs = rgbaGradient (0,0,1,1) (1,0,0,1) niso-        lyt = toRenderable . execEC $ do--            layout_title .= "Contours of a 2D Sin Curve"--            sequence $ do--                ((_,cntr),clr) <- zip cntrs clrs--                return . plot . liftEC $ do--                    plot_lines_style .= solidLine 3 clr-                    plot_lines_values .= cntr--    renderableToAspectWindow False 800 600 lyt