text:
{-# LANGUAGE BangPatterns, LambdaCase #-}

module Main where

import           Help
import           Input
import           Fractals
import           Colors
import           Data.Complex
import qualified Graphics.Image                as I
import qualified Graphics.Image.Interface      as I
import           System.IO
import           Data.Word
import           System.Environment
import           System.Directory
import           Data.Time

-- Entry point for program
main :: IO ()
main = do
   args <- getArgs
   if not (length args == 0) &&
      (  head args == "h"
      || head args == "-h"
      || head args == "help"
      || head args == "-help"
      || head args == "?"
      || head args == "-?"
      )
     then do
        putStrLn $ help args
        hFlush stdout
     else do
        let options = processArgs args
        putStrLn $ pPrintOptions options
        hFlush stdout
        fractal args options

fractal :: [String] -> Options -> IO ()
fractal args options@(Options fractalType (numRows, numColumns) (centerReal :+ centerImag)
                      range iterations power aaIn normalization color animation
                      cValue setColor framesIn startingFrame)
   = do
      currentDir <- getCurrentDirectory
      {- Path is generated by finding the current date/time, replacing colons (an invalid character in
       Windows paths) with an alternate Unicode colon, and removing the last 12 digits (which contain
       the fractional part of the seconds and the "UTC" code) -}
      path  <- fmap
         ( ((currentDir ++ "\\Generated Fractals\\") ++)
         . reverse . drop 12 . reverse
         . map (\case ':' -> '\xA789'; x -> x)
         . show
         )
         getCurrentTime
      createDirectoryIfMissing True path
      setCurrentDirectory path
      writeFile "data.txt" $ (pPrintOptions options) ++ '\n' : "Arguments:       " ++ (unwords args)
      putStr $ (show $ startingFrame - 1) ++ '/' : (show $ floor frames)
      hFlush stdout
      mapM_
         (\frame -> do
            I.writeImageExact (I.PNG) []
                  (path ++ '\\' : (take ((countDigits $ floor frames) - (countDigits $ frame)) $ repeat '0')
                   ++ show (frame) ++ ".png") $
                  I.toManifest $ I.toWord8I $ I.makeImageR I.RPS
                                       (numRows, numColumns)
                                       (\point -> colorFunc $ fractalFunction frame point)
            putStr $ '\r' : (show $ frame + 1) ++ '/' : (show $ floor frames)
            hFlush stdout
         )
         [startingFrame - 1 .. (truncate frames) - 1]
      putStrLn ""
 where
   frames = case animation of
      NoAnimation -> 1
      otherwise   -> fromIntegral framesIn
   colorFunc =
      case color of
            Greyscale                  -> colorGrey setColor
            Hue                        -> colorHue setColor
            Gradient isCircular colors -> colorGrad setColor isCircular colors
         . fmap
            (case normalization of
               Linear               -> normLinear iterations
               Sigmoid center power -> normSigmoid center power
               Periodic period      -> normPeriodic period
               Sine     period      -> normSine period
            )
   aa = case aaIn of
      AAEnabled  -> True
      AADisabled -> False
   animVals = map
      (case animation of
         Zoom final _ -> ((exp $ (log (final / range)) / (frames - 1)) **)
         Grid rIni rFnl _ _ rNum _ ->
            interpolate rIni rFnl . (/ (fromIntegral rNum - 1)) . fromIntegral . (`mod` rNum) . floor
         otherwise ->
            case animation of
                NoAnimation      -> \x -> 1
                Power final      -> interpolate power final
                Iterations final -> interpolate (fromIntegral iterations) (fromIntegral final)
                Theta final _    -> interpolate (phase cValue) (final * pi / 180)
                LinearC final    -> interpolate (realPart cValue) (realPart final)
            . if frames /= 1 then (/ (frames - 1)) else id
      )
      [fromIntegral startingFrame - 1 .. frames - 1]
   altAnimVals = map
      (case animation of
         Grid _ _ iIni iFnl rNum iNum ->
            interpolate iIni iFnl . (/ (fromIntegral iNum - 1)) . fromIntegral . (`div` rNum) . floor
         otherwise ->
            case animation of
                Zoom _ (Just finalIter) -> interpolate (fromIntegral iterations) (fromIntegral finalIter)
                Zoom _ Nothing          -> \x -> fromIntegral iterations
                Theta _ (Just final)    -> interpolate (magnitude cValue) final
                Theta _ Nothing         -> \x -> magnitude cValue
                LinearC final           -> interpolate (imagPart cValue) (imagPart final)
            . if frames /= 1 then (/ (frames - 1)) else id
      )
      [fromIntegral startingFrame - 1 .. frames - 1]
   fractalFunction :: Int -> (Int, Int) -> Maybe Double
   fractalFunction frame (r, c) = case (fractalType, animation) of
      -- Mandelbrot animation functions
      (Mandelbrot, NoAnimation) ->
         pMandelbrot aa pixelSize iterations power $ pairToComplex (r, c)
      (Mandelbrot, Power _) ->
         pMandelbrot aa pixelSize iterations value $ pairToComplex (r, c)
      (Mandelbrot, Zoom _ _) ->
         pMandelbrot aa (pixelSize * value) (round $ altValue) power $ pairToComplexZ (r, c) value
      (Mandelbrot, Iterations _) ->
         pMandelbrot aa pixelSize (round value) power $ pairToComplex (r, c)
      (Mandelbrot, Theta _ _) ->
         error "Theta animations can only be generated for Julia fractals"
      (Mandelbrot, LinearC _) ->
         error "Linear C-Value animations can only be generated for Julia fractals"
      (Mandelbrot, Grid _ _ _ _ _ _) ->
         error "Grided c-values can only be generated for Julia fractals"
      -- Julia animation functions
      (Julia, NoAnimation) ->
         pJulia aa pixelSize iterations power cValue $ pairToComplex (r, c)
      (Julia, Power _) ->
         pJulia aa pixelSize iterations value cValue $ pairToComplex (r, c)
      (Julia, Zoom _ _) ->
         pJulia aa (pixelSize * value) (round $ altValue) power cValue $ pairToComplexZ (r, c) value
      (Julia, Iterations _) ->
         pJulia aa pixelSize (round value) power cValue $ pairToComplex (r, c)
      (Julia, Theta _ _) ->
         pJulia aa pixelSize iterations power (mkPolar altValue value) $ pairToComplex (r, c)
      (Julia, LinearC _) ->
         pJulia aa pixelSize iterations power (value :+ altValue) $ pairToComplex (r, c)
      (Julia, Grid _ _ _ _ _ _) ->
         pJulia aa pixelSize iterations power (value :+ altValue) $ pairToComplex (r, c)
      where
         -- The value being animated
         value     = animVals !! frame
         altValue  = altAnimVals !! frame
         -- Other values
         pixelSize = (2 * range) / (fromIntegral $ numColumns - 1)
         halfV     = 0.5 * fromIntegral numRows
         halfH     = 0.5 * fromIntegral numColumns
         -- Function to get a complex point from a pixel coordinate pair
         pairToComplex (r, c) =
            (centerReal - range + fromIntegral c * pixelSize)
               :+ (centerImag + (pixelSize * (fromIntegral (numRows - 1) / 2)) - fromIntegral r * pixelSize)
         -- Same as above but for zoom animations
         pairToComplexZ (r, c) zoomFactor =
            ((fromIntegral c - halfH) * pixelSize * zoomFactor + centerReal)
               :+ ((fromIntegral r - halfV) * pixelSize * zoomFactor - centerImag)

-- Helper function for interpolating from the start to the end of an animation's value range
interpolate :: RealFloat a => a -> a -> a -> a
interpolate i f = (i +) . (*) (f - i)

-- Helper funtion to count the number of digits in an integer when expressed in base 10
countDigits :: Integral a => a -> a
countDigits = (+ 1) . floor . logBase 10 . fromIntegral

max_stars_repo_name: MaygeKyatt/Fractals

max_issues_repo_name: MaygeKyatt/Fractals

max_forks_repo_name: MaygeKyatt/Fractals


text:
{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}
{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
{-# OPTIONS_GHC -fno-warn-name-shadowing #-}

{-# LANGUAGE ConstraintKinds     #-}
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE GADTs               #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TemplateHaskell     #-}
{-# LANGUAGE TypeOperators       #-}
{-# LANGUAGE TypeApplications    #-}

module Test.Grenade.Layers.Transpose where

import           Grenade.Core
import           Grenade.Layers.Transpose
import           Data.Serialize
import           Data.Either
import           Test.Utils.Rnf

import           Numeric.LinearAlgebra        hiding (R, konst, uniformSample,
                                               (===))
import           Numeric.LinearAlgebra.Static as H hiding ((===), create)
import           System.Random.MWC             (create)
import           Hedgehog

prop_transpose_layer_calls_transpose_correctly = withTests 1 $ property $ do
  -- TEST 1
  let input :: L 12 4 = H.fromList [ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47 ]
      perms :: R 4    = H.fromList [ 2, 0, 1, 3 ]
      expected_output = (12 >< 4)  [ 0, 1, 2, 3, 12, 13, 14, 15, 24, 25, 26, 27, 36, 37, 38, 39, 4, 5, 6, 7, 16, 17, 18, 19, 28, 29, 30, 31, 40, 41, 42, 43, 8, 9, 10, 11, 20, 21, 22, 23, 32, 33, 34, 35, 44, 45, 46, 47 ]

      input' = S4D input :: S ('D4 2 2 3 4)
      layer = Transpose perms :: Transpose 4 ('D4 2 2 3 4) ('D4 3 2 2 4)

      S4D output = snd $ runForwards layer input' :: S ('D4 3 2 2 4)

  H.extract output === expected_output

  -- TEST 2
  let input :: L 12 4 = H.fromList [ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47 ]
      perms :: R 4    = H.fromList [ 3, 2, 1, 0 ]
      expected_output = (24 >< 2)  [ 0, 24, 12, 36, 4, 28, 16, 40, 8, 32, 20, 44, 1, 25, 13, 37, 5, 29, 17, 41, 9, 33, 21, 45, 2, 26, 14, 38, 6, 30, 18, 42, 10, 34, 22, 46, 3, 27, 15, 39, 7, 31, 19, 43, 11, 35, 23, 47 ]

      input' = S4D input :: S ('D4 2 2 3 4)
      layer = Transpose perms :: Transpose 4 ('D4 2 2 3 4) ('D4 4 3 2 2)

      S4D output = snd $ runForwards layer input' :: S ('D4 4 3 2 2)

  H.extract output === expected_output

prop_transpose_can_update_and_batch = withTests 1 $ property $ do
  gen <- evalIO create
  layer :: Transpose 4 ('D4 1 2 3 4) ('D4 4 3 2 1) <- evalIO $ createRandomWith UniformInit gen
  runUpdate defSGD layer () `seq` success
  runUpdate defAdam layer () `seq` success
  reduceGradient @(Transpose 4 ('D4 1 2 3 4) ('D4 4 3 2 1)) [()] `seq` success

prop_can_show_transpose = withTests 1 $ property $ do
  gen <- evalIO create
  layer :: Transpose 4 ('D4 1 2 3 4) ('D4 4 3 2 1) <- evalIO $ createRandomWith UniformInit gen
  show layer `seq` success

prop_transpose_rnf = withTests 1 $ property $ do
  let layer  = Transpose (H.fromList [3, 2, 1, 0]) :: Transpose 4 ('D4 1 2 3 4) ('D4 4 3 2 1)
      layer' = Transpose (error expectedErrStr) :: Transpose 4 ('D4 1 2 3 4) ('D4 4 3 2 1)
  err0 <- evalIO $ tryEvalRnf layer'
  assert $ rnfRaisedErr err0

  noErrStr <- evalIO $ tryEvalRnf layer
  assert $ rnfNoError noErrStr

prop_transpose_is_serializable = withTests 1 $ property $ do
  gen <- evalIO create
  layer :: Transpose 4 ('D4 1 2 3 4) ('D4 4 3 2 1) <- evalIO $ createRandomWith UniformInit gen
  let Transpose perms = layer
      bs = encode layer
      dec = decode bs :: Either String (Transpose 4 ('D4 1 2 3 4) ('D4 4 3 2 1))
  assert $ isRight dec
  let Right (Transpose perms') = dec
  (H.extract perms') === (H.extract perms)

tests :: IO Bool
tests = checkParallel $$(discover)

ext: hs

max_stars_repo_name: th-char/grenade

max_issues_repo_name: th-char/grenade

max_forks_repo_name: th-char/grenade


text:
module Main where

-- import           Lib

import qualified Control.Exception.Safe as Exception
import Data.Complex (Complex (..))
import qualified Data.Complex as Complex
import Data.Foldable (foldlM)
import Data.Functor (($>))
import qualified Data.Generics.Sum as GSum
import qualified Data.List as List
import Data.Map.Strict ((!?))
import qualified Data.Map.Strict as Map
import qualified Data.Maybe as Maybe
import Data.Ratio ((%))
import qualified Data.Ratio as Ratio
import Data.Text (Text)
import qualified Data.Text as Text
import qualified Data.Text.IO as Text
import Data.Void (Void)
import GHC.Generics (Generic)
import qualified Optics as Op
import System.Environment (getArgs)

-- import Text.Megaparsec
-- import qualified Text.Megaparsec.Char as MChar
-- import qualified Text.Megaparsec.Char.Lexer as L

data LispVal
  = Atom Text
  | List [LispVal]
  | DottedList [LispVal] LispVal
  | Number LispNumber
  | Float Float
  | String Text
  | Bool Bool
  deriving (Generic, Eq, Show)

data LispNumber
  = Complex LispComplex
  | Real Double
  | Rational Rational
  | Integer Integer
  deriving (Generic, Eq, Show)

type LispComplex = Complex Double

-- lispNumComplex = Op.iso toComplex toLispNumber

lispNumToComplex :: LispNumber -> LispComplex
lispNumToComplex (Complex x) = x
lispNumToComplex (Real x) = x :+ 0
lispNumToComplex (Rational x) = fromRational x :+ 0
lispNumToComplex (Integer x) = fromInteger x :+ 0

lispNumToDouble :: Exception.MonadThrow m => LispNumber -> m Double
lispNumToDouble x@(Complex _) = Exception.throw [TypeMismatch "Real" $ Number x]
lispNumToDouble (Real x) = pure x
lispNumToDouble (Rational x) = pure $ fromRational x
lispNumToDouble (Integer x) = pure $ fromInteger x

lispNumToRational :: Exception.MonadThrow m => LispNumber -> m Rational
lispNumToRational x@(Complex _) = Exception.throw [TypeMismatch "Rational" $ Number x]
lispNumToRational x@(Real _) = Exception.throw [TypeMismatch "Rational" $ Number x]
lispNumToRational (Rational x) = pure x
lispNumToRational (Integer x) = pure $ fromInteger x

lispNumToInteger :: Exception.MonadThrow m => LispNumber -> m Integer
lispNumToInteger x@(Complex _) = Exception.throw [TypeMismatch "Integer" $ Number x]
lispNumToInteger x@(Real _) = Exception.throw [TypeMismatch "Integer" $ Number x]
lispNumToInteger x@(Rational _) = Exception.throw [TypeMismatch "Integer" $ Number x]
lispNumToInteger (Integer x) = pure x

-- toLispNumber :: LispComplex -> LispNumber
-- toLispNumber (r :+ 0)
--   | Ratio.denominator r == 0 = Integer $ Ratio.numerator r
--   | otherwise = Rational r
-- toLispNumber c = Complex c

-- apNumOpToLispNumber ::
--   (Complex Rational -> Complex Rational -> Complex Rational) ->
--   (LispNumber -> LispNumber -> LispNumber)
-- apNumOpToLispNumber op x y = toLispNumber $ toComplex x `op` toComplex y

-- instance Num LispNumber where
--   (+) = apNumOpToLispNumber (+)
--   (-) = apNumOpToLispNumber (-)
--   (*) = apNumOpToLispNumber (*)
--   abs = Op.over lispNumComplex abs
--   signum = Op.over lispNumComplex signum
--   fromInteger = toLispNumber . fromInteger

-- instance Fractional LispNumber where
--   fromRational = toLispNumber . fromRational
--   recip = Op.over lispNumComplex recip

data LispError
  = NumArgs Integer [LispVal]
  | TypeMismatch Text LispVal
  | --   | Parser ParseError
    BadSpecialForm Text LispVal
  | -- | NotFunction T.Text T.Text
    UnboundVar Text
  | Default Text
  deriving (Generic, Eq, Show)

instance Exception.Exception [LispError]

type LispThrowErr = Either [LispError]

-- type Parser = Parsec Void T.Text

-- symbolChar :: Parser Char
-- symbolChar = oneOf "!#$%&|*+-/:<=>?@^_~"

-- spaceConsumer :: Parser ()
-- spaceConsumer = L.space MChar.space1 empty empty

-- lexeme :: Parser a -> Parser a
-- lexeme = L.lexeme spaceConsumer

-- symbol :: T.Text -> Parser T.Text
-- symbol = L.symbol' spaceConsumer

-- -- readExpr :: T.Text -> T.Text
-- -- readExpr input =
-- --   case parse (space >> symbolChar) "lisp" input of
-- --     Left err -> "No match: " <> T.pack (errorBundlePretty err)
-- --     Right val -> "Found value: " <> val

-- parseExpr :: Parser LispVal
-- parseExpr = lexeme $ Atom . T.singleton <$> symbolChar

-- parseString :: Parser LispVal
-- parseString = do
--   x <- between (symbol "\"") (symbol "\"") $ takeWhileP Nothing (/= '\'')
--   return $ String x

-- parseAtom :: Parser LispVal
-- parseAtom = do
--   first <- MChar.letterChar <|> symbolChar
--   rest <- many (letter <|> digit <|> symbol)
--   let atom = first : rest
--   return $ case atom of
--     "#t" -> Bool True
--     "#f" -> Bool False
--     _ -> Atom atom

-- parseNumber :: Parser LispVal
-- parseNumber = Number . read <$> many1 digit

foldNumOp :: Exception.MonadThrow m => (LispNumber -> LispNumber -> m LispNumber) -> LispNumber -> [LispVal] -> m LispVal
foldNumOp op identity ls =
  let mNumbers = Op.preview (GSum._As @"Number") <$> ls
   in case List.elemIndices Nothing mNumbers of
        [] -> Number <$> foldlM op identity (Maybe.catMaybes mNumbers)
        indices -> Exception.throw $ TypeMismatch "Number" . (ls !!) <$> indices

binNumOp :: Exception.MonadThrow m => (LispNumber -> LispNumber -> m LispNumber) -> [LispVal] -> m LispVal
binNumOp op [Number n, Number m] = Number <$> op n m
binNumOp op args@[_, _] =
  Exception.throw $
    TypeMismatch "Number"
      <$> List.filter (\case Number _ -> True; _ -> False) args
binNumOp op args = Exception.throw [NumArgs 2 args]

numOpToLispNumOp ::
  Exception.MonadThrow m =>
  (forall a. Num a => a -> a -> a) ->
  (LispNumber -> LispNumber -> m LispNumber)
numOpToLispNumOp op (Integer x) (Integer y) = pure $ Integer $ op x y
numOpToLispNumOp op (Integer x) (Rational y) = pure $ Rational $ op (fromIntegral x) y
numOpToLispNumOp op (Integer x) (Real y) = pure $ Real $ op (fromIntegral x) y
numOpToLispNumOp op (Rational x) (Integer y) = pure $ Rational $ op x $ fromIntegral y
numOpToLispNumOp op (Rational x) (Rational y) = pure $ Rational $ op x y
numOpToLispNumOp op (Rational x) (Real y) = pure $ Real $ op (fromRational x) y
numOpToLispNumOp op (Real x) (Integer y) = pure $ Real $ op x $ fromIntegral y
numOpToLispNumOp op (Real x) (Rational y) = pure $ Real $ op x $ fromRational y
numOpToLispNumOp op (Real x) (Real y) = pure $ Real $ op x y
numOpToLispNumOp op x y = pure $ Complex $ op (lispNumToComplex x) (lispNumToComplex y)

fractionalOpToLispNumOp ::
  Exception.MonadThrow m =>
  (forall a. Fractional a => a -> a -> a) ->
  (LispNumber -> LispNumber -> m LispNumber)
fractionalOpToLispNumOp op (Integer x) (Integer y) = pure $ Rational $ op (fromIntegral x) (fromIntegral y)
fractionalOpToLispNumOp op (Integer x) (Rational y) = pure $ Rational $ op (fromIntegral x) y
fractionalOpToLispNumOp op (Integer x) (Real y) = pure $ Real $ op (fromIntegral x) y
fractionalOpToLispNumOp op (Rational x) (Integer y) = pure $ Rational $ op x $ fromIntegral y
fractionalOpToLispNumOp op (Rational x) (Rational y) = pure $ Rational $ op x y
fractionalOpToLispNumOp op (Rational x) (Real y) = pure $ Real $ op (fromRational x) y
fractionalOpToLispNumOp op (Real x) (Integer y) = pure $ Real $ op x $ fromIntegral y
fractionalOpToLispNumOp op (Real x) (Rational y) = pure $ Real $ op x $ fromRational y
fractionalOpToLispNumOp op (Real x) (Real y) = pure $ Real $ op x y
fractionalOpToLispNumOp op x y = pure $ Complex $ op (lispNumToComplex x) (lispNumToComplex y)

integerOpToLispNumOp ::
  Exception.MonadThrow m =>
  (forall a. Integral a => a -> a -> a) ->
  (LispNumber -> LispNumber -> m LispNumber)
integerOpToLispNumOp op (Integer x) (Integer y) = pure $ Integer $ op x y
integerOpToLispNumOp _ x y = Exception.throw $ TypeMismatch "Integer" . Number <$> [x, y]

primitivesEnv :: Exception.MonadThrow m => Map.Map Text ([LispVal] -> m LispVal)
primitivesEnv =
  Map.fromList
    [ ("+", foldNumOp (numOpToLispNumOp (+)) $ Integer 0),
      ("-", foldNumOp (numOpToLispNumOp (-)) $ Integer 0),
      ("*", foldNumOp (numOpToLispNumOp (*)) $ Integer 1),
      ("/", foldNumOp (fractionalOpToLispNumOp (/)) $ Integer 1),
      ("quotient", binNumOp $ integerOpToLispNumOp quot),
      ("remainder", binNumOp $ integerOpToLispNumOp rem),
      ("modulo", binNumOp $ integerOpToLispNumOp mod)
    ]

eval :: Exception.MonadThrow m => LispVal -> m LispVal
eval (List [Atom "quote", ls]) = pure ls
eval (List (Atom func : args)) = do
  args <- traverse eval args
  apply func args
eval val = pure val

apply :: forall m. Exception.MonadThrow m => Text -> [LispVal] -> m LispVal
apply funcName args =
  case primitivesEnv !? funcName :: Maybe ([LispVal] -> m LispVal) of
    Nothing -> Exception.throw [UnboundVar funcName]
    Just f -> f args

main :: IO ()
main = do
  print $
    (eval :: LispVal -> Either Exception.SomeException LispVal)
      <$> [ List [Atom "+", Number $ Integer 1, Number $ Integer 3],
            List [Atom "*", Number $ Integer 1, Number $ Real 3, Number $ Rational 1, Number $ Complex 3],
            List [Atom "/", Number $ Real 1, Number $ Real 3, Number $ Rational 2, Number $ Integer 3],
            List [Atom "remainder", Number $ Integer 12, Number $ Integer 5]
          ]


text:
{-# LANGUAGE RecordWildCards #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE OverloadedStrings #-}

module StatReport (statReport, showStatEntryValue) where

import Data.Fixed (showFixed)
import Data.Text (Text)
import Fmt

import QuoteData
import Statistics

instance Buildable Statistic where
  build Mean = "Mean"
  build Min = "Minimum"
  build Max = "Maximum"
  build Days = "Days between Min/Max"

showStatEntryValue :: StatEntry -> String
showStatEntryValue StatEntry {..} = showFixed (removeTrailing stat qfield) value
  where
    removeTrailing Days _ = True
    removeTrailing Min Volume = True
    removeTrailing Max Volume = True
    removeTrailing _ _ = False

instance Buildable StatEntry where
  build se@StatEntry {..} = ""+|stat|+": "+|showStatEntryValue se|+""

instance Buildable StatQFieldData where
  build (qf, stats) = nameF ("Statistics for " +||qf||+"") $ unlinesF stats

statReport :: StatInfo -> Text
statReport = fmt . unlinesF

max_stars_repo_name: Francososa/stockell

max_issues_repo_name: Francososa/stockell

max_forks_repo_name: Francososa/stockell


text:
{-# LANGUAGE DeriveDataTypeable, DeriveGeneric #-}
-- |
-- Module    : Statistics.Distribution.Exponential
-- Copyright : (c) 2009 Bryan O'Sullivan
-- License   : BSD3
--
-- Maintainer  : bos@serpentine.com
-- Stability   : experimental
-- Portability : portable
--
-- The exponential distribution.  This is the continunous probability
-- distribution of the times between events in a poisson process, in
-- which events occur continuously and independently at a constant
-- average rate.

module Statistics.Distribution.Exponential
    (
      ExponentialDistribution
    -- * Constructors
    , exponential
    , exponentialFromSample
    -- * Accessors
    , edLambda
    ) where

import Data.Aeson (FromJSON, ToJSON)
import Data.Binary (Binary)
import Data.Data (Data, Typeable)
import GHC.Generics (Generic)
import Numeric.MathFunctions.Constants (m_neg_inf)
import qualified Statistics.Distribution         as D
import qualified Statistics.Sample               as S
import qualified System.Random.MWC.Distributions as MWC
import Statistics.Types (Sample)
import Data.Binary (put, get)

newtype ExponentialDistribution = ED {
      edLambda :: Double
    } deriving (Eq, Read, Show, Typeable, Data, Generic)

instance FromJSON ExponentialDistribution
instance ToJSON ExponentialDistribution

instance Binary ExponentialDistribution where
    put = put . edLambda
    get = fmap ED get

instance D.Distribution ExponentialDistribution where
    cumulative      = cumulative
    complCumulative = complCumulative

instance D.ContDistr ExponentialDistribution where
    density (ED l) x
      | x < 0     = 0
      | otherwise = l * exp (-l * x)
    logDensity (ED l) x
      | x < 0     = m_neg_inf
      | otherwise = log l + (-l * x)
    quantile = quantile

instance D.Mean ExponentialDistribution where
    mean (ED l) = 1 / l

instance D.Variance ExponentialDistribution where
    variance (ED l) = 1 / (l * l)

instance D.MaybeMean ExponentialDistribution where
    maybeMean = Just . D.mean

instance D.MaybeVariance ExponentialDistribution where
    maybeStdDev   = Just . D.stdDev
    maybeVariance = Just . D.variance

instance D.Entropy ExponentialDistribution where
  entropy (ED l) = 1 - log l

instance D.MaybeEntropy ExponentialDistribution where
  maybeEntropy = Just . D.entropy

instance D.ContGen ExponentialDistribution where
  genContVar = MWC.exponential . edLambda

cumulative :: ExponentialDistribution -> Double -> Double
cumulative (ED l) x | x <= 0    = 0
                    | otherwise = 1 - exp (-l * x)

complCumulative :: ExponentialDistribution -> Double -> Double
complCumulative (ED l) x | x <= 0    = 1
                         | otherwise = exp (-l * x)

quantile :: ExponentialDistribution -> Double -> Double
quantile (ED l) p
  | p == 1          = 1 / 0
  | p >= 0 && p < 1 = -log (1 - p) / l
  | otherwise       =
    error $ "Statistics.Distribution.Exponential.quantile: p must be in [0,1] range. Got: "++show p

-- | Create an exponential distribution.
exponential :: Double            -- ^ λ (scale) parameter.
            -> ExponentialDistribution
exponential l
  | l <= 0 =
    error $ "Statistics.Distribution.Exponential.exponential: scale parameter must be positive. Got " ++ show l
  | otherwise = ED l

-- | Create exponential distribution from sample. No tests are made to
-- check whether it truly is exponential.
exponentialFromSample :: Sample -> ExponentialDistribution
exponentialFromSample = ED . S.mean

max_stars_repo_name: StefanHubner/statistics

max_issues_repo_name: StefanHubner/statistics

max_forks_repo_name: StefanHubner/statistics


text:
module RadialPhase where

import           Control.Monad           as M
import           Control.Monad.Parallel  as MP
import           Data.Array.Repa         as R
import           Data.Binary
import           Data.Complex
import           Data.List               as L
import           Data.Vector.Unboxed     as VU
import           DFT.Plan
import           Filter.Utils
import           FokkerPlanck.MonteCarlo
import           FokkerPlanck.Pinwheel
import           Graphics.Gnuplot.Simple
import           Image.IO
import           System.Directory
import           System.Environment
import           System.FilePath
import           Text.Printf
import           Types
import           Utils.Array

main = do
  args@(numPointStr:thetaSigmaStr:scaleSigmaStr:maxScaleStr:maxScaleStr1:taoStr:numTrailStr:maxTrailStr:theta0FreqsStr:thetaFreqsStr:scale0FreqsStr:scaleFreqsStr:histFilePath:alphaStr:deltaThetaStr:numThreadStr:_) <-
    getArgs
  print args
  let numPoint = read numPointStr :: Int
      thetaSigma = read thetaSigmaStr :: Double
      scaleSigma = read scaleSigmaStr :: Double
      maxScale = read maxScaleStr :: Double
      maxScale1 = read maxScaleStr1 :: Double
      tao = read taoStr :: Double
      numTrail = read numTrailStr :: Int
      maxTrail = read maxTrailStr :: Int
      theta0Freq = read theta0FreqsStr :: Double
      theta0Freqs = [-theta0Freq .. theta0Freq]
      thetaFreq = read thetaFreqsStr :: Double
      thetaFreqs = [-thetaFreq .. thetaFreq]
      scale0Freq = read scale0FreqsStr :: Double
      scaleFreq = read scaleFreqsStr :: Double
      scale0Freqs = [-scale0Freq .. scale0Freq]
      scaleFreqs = [-scaleFreq .. scaleFreq]
      alpha = read alphaStr :: Double
      deltaTheta = read deltaThetaStr :: Double
      numThread = read numThreadStr :: Int
      folderPath = "output/test/RadialPhase"
  createDirectoryIfMissing True folderPath
  removePathForcibly (folderPath </> "Magnitude")
  createDirectoryIfMissing True (folderPath </> "Magnitude")
  flag <- doesFileExist histFilePath
  radialArr <-
    if flag
      then R.map magnitude . getNormalizedHistogramArr <$>
           decodeFile histFilePath
      else solveMonteCarloR2Z2T0S0Radial
             numThread
             numTrail
             maxTrail
             numPoint
             numPoint
             thetaSigma
             scaleSigma
             maxScale
             tao
             theta0Freqs
             thetaFreqs
             scale0Freqs
             scaleFreqs
             histFilePath
             (emptyHistogram
                [ (round . sqrt . fromIntegral $ 2 * (div numPoint 2) ^ 2)
                , L.length scale0Freqs
                , L.length theta0Freqs
                , L.length scaleFreqs
                , L.length thetaFreqs
                ]
                0)
  let arr3d =
        R.slice
          radialArr
          (Z :. (L.length thetaFreqs - 1) :. All :. (L.length theta0Freqs - 1) :.
           All :.
           All)
      (Z :. _ :. _ :. _ :. _ :. r) = extent radialArr
      idx =
        [ (tf, sf, t0f, s0f)
        | tf <- L.zip [0 .. (L.length thetaFreqs) - 1] [-3,-2,-1] --[(-thetaFreq) .. 0]
        , sf <- L.zip [0 .. (L.length scaleFreqs) - 1] [(-scaleFreq) .. 0]
        , t0f <- L.zip [0 .. (L.length theta0Freqs) - 1] theta0Freqs
        , s0f <- L.zip [0 .. (L.length scale0Freqs) - 1] scale0Freqs
        ]
      idxs =
        L.groupBy
          (\((_, tf), (_, sf), (_, t0f), (_, s0f)) ((_, tf'), (_, sf'), (_, t0f'), (_, s0f')) ->
             (tf + t0f, sf + s0f) == (tf' + t0f', sf' + s0f')) .
        L.sortBy
          (\((_, tf), (_, sf), (_, t0f), (_, s0f)) ((_, tf'), (_, sf'), (_, t0f'), (_, s0f')) ->
             compare (tf + t0f, sf + s0f) (tf' + t0f', sf' + s0f')) $
        idx
      -- idxs =
      --   L.groupBy
      --     (\((_, tf), (_, sf), (_, t0f), (_, s0f)) ((_, tf'), (_, sf'), (_, t0f'), (_, s0f')) ->
      --        (abs $ tf + t0f, abs $ tf - t0f) == (abs $ tf' + t0f', abs $ tf' - t0f')) .
      --   L.sortBy
      --     (\((_, tf), (_, sf), (_, t0f), (_, s0f)) ((_, tf'), (_, sf'), (_, t0f'), (_, s0f')) ->
      --        compare (abs $ tf + t0f, abs $ tf - t0f) (abs $ tf' + t0f', abs $ tf' - t0f')) $
      --   idx
  -- MP.mapM_
  --   (\((i, sf), (j, s0f)) ->
  --      plotPathStyle
  --        [ PNG (folderPath </> printf "%d_%d.png" i j)
  --        , Title (printf "sf: %.0f, s0f: %.0f" sf s0f)
  --        -- , XRange (0, log (fromIntegral len))
  --        , YRange (0, 2 * pi)
  --        ]
  --        (defaultStyle
  --           {plotType = LinesPoints, lineSpec = CustomStyle [LineTitle ""]}) .
  --      L.zip (L.map log [1 .. fromIntegral len - 1]) .
  --      L.tail . R.toList . R.map (\x -> (phase x) + pi) . R.slice arr3d $
  --      (Z :. i :. j :. All))
  --   [ (a, b)
  --   | a <- L.zip [0 .. (L.length scaleFreqs) - 1] scaleFreqs
  --   , b <- L.zip [0 .. (L.length scale0Freqs) - 1] scale0Freqs
  --   ]
  -- let maxMag = VU.maximum . toUnboxed . computeS . R.map magnitude $ arr3d
  plotPathsStyle
    [ PNG (folderPath </> "Magnitude" </> "all.png")
    , Title ""
    , XRange (60, 120)
    , YRange (0.01, 0.015)
    ] (L.map
         (\((i, tf), (j, sf), (k, t0f), (l, s0f)) ->
            ( defaultStyle
                { plotType = LinesPoints
                , lineSpec =
                    CustomStyle
                      [ LineTitle
                          (printf "%.0f %.0f %.0f %.0f" tf t0f sf s0f)
                      ]
                }
            , L.drop 5 .
              L.zip ([0,5 .. 5*(fromIntegral r - 1)]) .
              R.toList . R.slice radialArr $
              (Z :. i :. j :. k :. l :. All)))
         idx)
  MP.mapM_
    (\xs@(((i, tf), (j, sf), (k, t0f), (l, s0f)):_) ->
       let t = tf + t0f
           s = sf + s0f
        in plotPathsStyle
             [ PNG (folderPath </> "Magnitude" </> printf "%.0f_%.0f.png" t s)
             , Title (printf "tSum: %.0f sSum: %.0f" t s)
             , YRange (0, 0.045)
             ]
             (L.map
                (\((i, tf), (j, sf), (k, t0f), (l, s0f)) ->
                   ( defaultStyle
                       { plotType = LinesPoints
                       , lineSpec =
                           CustomStyle
                             [ LineTitle
                                 (printf "%.0f %.0f %.0f %.0f" tf t0f sf s0f)
                             ]
                       }
                   , L.drop 5 .
                     L.zip ([0 .. fromIntegral r - 1]) .
                     R.toList . R.slice radialArr $
                     (Z :. i :. j :. k :. l :. All)))
                xs))
    idxs
  -- MP.mapM_
  --   (M.mapM_
  --      (\((i, tf), (j, sf), (k, t0f), (l, s0f)) -> do
  --         createDirectoryIfMissing
  --           True
  --           (folderPath </> "Magnitude" </>
  --            (printf "%.0f_%.0f" (tf + t0f) (sf + s0f)))
  --         plotPathStyle
  --           [ PNG
  --               (folderPath </> "Magnitude" </>
  --                (printf "%.0f_%.0f" (tf + t0f) (sf + s0f)) </>
  --                printf "%.0f_%.0f_%.0f_%.0f.png" tf sf t0f s0f)
  --           , Title ""
  --               -- (printf "tf: %.0f sf: %.0f t0f: %.0f s0f: %.0f" tf sf t0f s0f)
  --           ]
  --           (defaultStyle
  --              {plotType = LinesPoints, lineSpec = CustomStyle [LineTitle ""]}) .
  --           L.zip ([0 .. fromIntegral r - 1]) . R.toList . R.slice radialArr $
  --           (Z :. i :. j :. k :. l :. All)))
  --   idxs
  -- MP.mapM_
  --   (\((i, tf), (j, sf), (k, t0f), (l, s0f)) ->
  --      plotPathStyle
  --        [ PNG
  --            (folderPath </> "Magnitude" </>
  --             printf "%.0f_%.0f_%.0f_%.0f.png" tf sf t0f s0f)
  --        , Title (printf "tf: %.0f sf: %.0f t0f: %.0f s0f: %.0f" tf sf t0f s0f)
  --        -- , XRange (0, log (fromIntegral len))
  --        -- , YRange (0, maxMag)
  --        ]
  --        (defaultStyle
  --           {plotType = LinesPoints, lineSpec = CustomStyle [LineTitle ""]}) .
  --      L.zip ([0 .. fromIntegral r - 1]) .
  --      -- L.tail .
  --      R.toList . R.slice radialArr $
  --      (Z :. i :. j :. k :. l :. All))
  --   [ (tf, sf, t0f, s0f)
  --   | tf <- L.zip [0 .. (L.length thetaFreqs) - 1] thetaFreqs
  --   , sf <- L.zip [0 .. (L.length scaleFreqs) - 1] scaleFreqs
  --   , t0f <- L.zip [0 .. (L.length theta0Freqs) - 1] theta0Freqs
  --   , s0f <- L.zip [0 .. (L.length scale0Freqs) - 1] scale0Freqs
  --   ]
  -- MP.mapM_
  --   (\((i, sf), (j, s0f)) ->
  --      plotPathStyle
  --        [ PNG (folderPath </> "Pinwheel" </> printf "%d_%d.png" i j)
  --        , Title (printf "sf: %.0f, s0f: %.0f" sf s0f)
  --        -- , XRange (0, log (fromIntegral len))
  --        , YRange (0, 2 * pi)
  --        ]
  --        (defaultStyle
  --           {plotType = LinesPoints, lineSpec = CustomStyle [LineTitle ""]}) .
  --      L.zip (L.map log [1 .. fromIntegral len - 1]) .
  --      L.map
  --        (\r ->
  --           pi +
  --           (phase $
  --            exp
  --              (0 :+
  --               (s0f *
  --                (log r * 2 * pi / log maxScale1 + (deltaTheta / 180 * pi)) +
  --                sf * (log r * 2 * pi / log maxScale1 + (deltaTheta / 180 * pi)))))) $
  --      [1 .. fromIntegral len - 1])
  --   [ (a, b)
  --   | a <- L.zip [0 .. (L.length scaleFreqs) - 1] scaleFreqs
  --   , b <- L.zip [0 .. (L.length scale0Freqs) - 1] scale0Freqs
  --   ]

max_stars_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_issues_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_forks_repo_name: XinhuaZhang/Stochastic-Completion-Field


text:
module Main where

import Correlation (cacf)
import Data.Complex (Complex (..))
import Graphics.Rendering.Chart.Backend.Diagrams
import Graphics.Rendering.Chart.Easy
import Linear.Metric
import Linear.V3
import Linear.Vector
import MLS (mls)

j :: RealFloat a => Complex a
j = 0 :+ 1

main :: IO ()
main = do
  print "IDDQD"
  print code
  toFile def "example1_big.png" $ do
    layout_title .= "Amplitude Modulation"
    setColors [opaque blue, opaque red]
    plot (line "am" [signal' [0, (0.5) .. 400]])
    plot (points "am points" (signal' [0, 7 .. 400]))

signal' :: [Double] -> [(Double, Double)]
signal' xs = [(x, (sin (x * 3.14159 / 45) + 1) / 2 * (sin (x * 3.14159 / 5))) | x <- xs]

code = mls [True, False, True, False, False, True, False, False, True, False, False, True] [True, False, False, False, False, False, False, False, False, False, False]

base = length code

carrierFeq :: Double
carrierFeq = 1e10

codeInterval :: Double
codeInterval = 1e-6

-- Теорема Котельникова
df :: Double
df = (carrierFeq + fromIntegral base / codeInterval) * 3

signal :: Double -> Complex Double
signal t
  | t >= 0 && t < signalLength = a t * exp (- j * (2 * pi * carrierFeq * t :+ 0))
  | otherwise = 0
  where
    signalLength = codeInterval * fromIntegral base
    i = floor (t / codeInterval) `mod` base
    a t =
      if code !! i
        then 1
        else -1

echo :: (Double -> V3 Double) -> Double -> Complex Double
echo obj t = signal (t - d / c)
  where
    d = distance (carrier t) (obj t)
    c = 3e8

scene :: [V3 Double]
scene = [V3 0 10000 0]

zeroPoint :: V3 Double
zeroPoint = V3 0 0 0

vel :: V3 Double
vel = V3 1 0 0

carrier :: Double -> V3 Double
carrier t = zeroPoint + vel ^* t

max_stars_repo_name: nihlete/sar-signal

max_issues_repo_name: nihlete/sar-signal

max_forks_repo_name: nihlete/sar-signal


text:
{-# LANGUAGE GADTs #-}
{-# LANGUAGE KindSignatures #-}
{-# LANGUAGE GADTs, DataKinds, PolyKinds, ScopedTypeVariables #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}
{-|
Module      : Grenade.Core.WeightInitialization
Description : Defines the Weight Initialization methods of Grenade.
Copyright   : (c) Manuel Schneckenreither, 2018
License     : BSD2
Stability   : experimental

This module defines the weight initialization methods.

-}

module Grenade.Core.WeightInitialization
    ( getRandomVector
    , getRandomMatrix
    , WeightInitMethod (..)
    ) where

import           Data.Proxy
import Control.Monad

import           Data.Singletons.TypeLits

import           Control.Monad.Primitive         (PrimBase, PrimState)
import           System.Random.MWC
import           System.Random.MWC.Distributions

import           GHC.TypeNats (type (*))

import           Numeric.LinearAlgebra.Static

data WeightInitMethod = UniformInit       -- ^ W_l,i ~ U(-1/sqrt(n_l),1/sqrt(n_l))                   where n_l is the number of nodes in layer l
                      | Xavier            -- ^ W_l,i ~ U(-sqrt (6/n_l+n_{l+1}),sqrt (6/n_l+n_{l+1})) where n_l is the number of nodes in layer l
                      | HeEtAl            -- ^ W_l,i ~ N(0,sqrt(2/n_l))                              where n_l is the number of nodes in layer l

-- | Get a random vector initialized according to the specified method.
getRandomVector :: forall m n . (PrimBase m, KnownNat n) => Integer -> Integer -> WeightInitMethod -> Gen (PrimState m) -> m (R n)
getRandomVector i o method gen = do
  unifRands <- vector <$> replicateM n (uniformR (-1,1) gen)
  gaussRands <- vector <$> replicateM n (standard gen)

  return $ case method of
             UniformInit -> (1/sqrt (fromIntegral i)) * unifRands
             Xavier      -> (sqrt 6/sqrt (fromIntegral i + fromIntegral o)) * unifRands
             HeEtAl      -> sqrt (2/fromIntegral i) * gaussRands
  where n = fromIntegral $ natVal (Proxy :: Proxy n)

-- | Get a matrix with weights initialized according to the specified method.
getRandomMatrix :: forall m r n nr . (PrimBase m, KnownNat r, KnownNat n, KnownNat nr, nr ~ (n*r))
                => Integer -> Integer -> WeightInitMethod -> Gen (PrimState m) -> m (L r n)
getRandomMatrix i o method gen = do
  unifRands <- matrix <$> replicateM nr (uniformR (-1,1) gen)
  gaussRands <- matrix <$> replicateM nr (standard gen)

  return $ case method of
             UniformInit -> (1/sqrt (fromIntegral i)) * unifRands
             Xavier -> (sqrt 6/sqrt (fromIntegral i + fromIntegral o)) * unifRands
             HeEtAl -> (sqrt (2/fromIntegral i)) * gaussRands

  where nr = fromIntegral $ natVal (Proxy :: Proxy nr)

max_stars_repo_name: koenigmaximilian/grenade

max_issues_repo_name: koenigmaximilian/grenade

max_forks_repo_name: koenigmaximilian/grenade


text:
-----------------------------------------------------------------------------
-- |
-- Module      :  Numeric.IO
-- Copyright   :  (c) Alberto Ruiz 2010
-- License     :  GPL
--
-- Maintainer  :  Alberto Ruiz <aruiz@um.es>
-- Stability   :  provisional
-- Portability :  portable
--
-- Display, formatting and IO functions for numeric 'Vector' and 'Matrix'
--
-----------------------------------------------------------------------------

module Numeric.IO (
    dispf, disps, dispcf, vecdisp, latexFormat, format,
    loadMatrix, saveMatrix, fromFile, fileDimensions,
    readMatrix, fromArray2D,
    fscanfVector, fprintfVector, freadVector, fwriteVector
) where

import Data.Packed
import Data.Packed.Internal
import System.Process(readProcess)
import Text.Printf(printf)
import Data.List(intersperse)
import Data.Complex

{- | Creates a string from a matrix given a separator and a function to show each entry. Using
this function the user can easily define any desired display function:

@import Text.Printf(printf)@

@disp = putStr . format \"  \" (printf \"%.2f\")@

-}
format :: (Element t) => String -> (t -> String) -> Matrix t -> String
format sep f m = table sep . map (map f) . toLists $ m

{- | Show a matrix with \"autoscaling\" and a given number of decimal places.

@disp = putStr . disps 2

\> disp $ 120 * (3><4) [1..]
3x4  E3
 0.12  0.24  0.36  0.48
 0.60  0.72  0.84  0.96
 1.08  1.20  1.32  1.44
@
-}
disps :: Int -> Matrix Double -> String
disps d x = sdims x ++ "  " ++ formatScaled d x

{- | Show a matrix with a given number of decimal places.

@disp = putStr . dispf 3

\> disp (1/3 + ident 4)
4x4
1.333  0.333  0.333  0.333
0.333  1.333  0.333  0.333
0.333  0.333  1.333  0.333
0.333  0.333  0.333  1.333
@
-}
dispf :: Int -> Matrix Double -> String
dispf d x = sdims x ++ "\n" ++ formatFixed (if isInt x then 0 else d) x

sdims x = show (rows x) ++ "x" ++ show (cols x)

formatFixed d x = format "  " (printf ("%."++show d++"f")) $ x

isInt = all lookslikeInt . toList . flatten

formatScaled dec t = "E"++show o++"\n" ++ ss
    where ss = format " " (printf fmt. g) t
          g x | o >= 0    = x/10^(o::Int)
              | otherwise = x*10^(-o)
          o = floor $ maximum $ map (logBase 10 . abs) $ toList $ flatten t
          fmt = '%':show (dec+3) ++ '.':show dec ++"f"

{- | Show a vector using a function for showing matrices.

@disp = putStr . vecdisp ('dispf' 2)

\> disp ('linspace' 10 (0,1))
10 |> 0.00  0.11  0.22  0.33  0.44  0.56  0.67  0.78  0.89  1.00
@
-}
vecdisp :: (Element t) => (Matrix t -> String) -> Vector t -> String
vecdisp f v
    = ((show (dim v) ++ " |> ") ++) . (++"\n")
    . unwords . lines .  tail . dropWhile (not . (`elem` " \n"))
    . f . trans . reshape 1
    $ v

-- | Tool to display matrices with latex syntax.
latexFormat :: String -- ^ type of braces: \"matrix\", \"bmatrix\", \"pmatrix\", etc.
            -> String -- ^ Formatted matrix, with elements separated by spaces and newlines
            -> String
latexFormat del tab = "\\begin{"++del++"}\n" ++ f tab ++ "\\end{"++del++"}"
    where f = unlines . intersperse "\\\\" . map unwords . map (intersperse " & " . words) . tail . lines

-- | Pretty print a complex number with at most n decimal digits.
showComplex :: Int -> Complex Double -> String
showComplex d (a:+b)
    | isZero a && isZero b = "0"
    | isZero b = sa
    | isZero a && isOne b = s2++"i"
    | isZero a = sb++"i"
    | isOne b = sa++s3++"i"
    | otherwise = sa++s1++sb++"i"
  where
    sa = shcr d a
    sb = shcr d b
    s1 = if b<0 then "" else "+"
    s2 = if b<0 then "-" else ""
    s3 = if b<0 then "-" else "+"

shcr d a | lookslikeInt a = printf "%.0f" a
         | otherwise      = printf ("%."++show d++"f") a

lookslikeInt x = show (round x :: Int) ++".0" == shx || "-0.0" == shx
   where shx = show x

isZero x = show x `elem` ["0.0","-0.0"]
isOne  x = show x `elem` ["1.0","-1.0"]

-- | Pretty print a complex matrix with at most n decimal digits.
dispcf :: Int -> Matrix (Complex Double) -> String
dispcf d m = sdims m ++ "\n" ++ format "  " (showComplex d) m

--------------------------------------------------------------------

-- | reads a matrix from a string containing a table of numbers.
readMatrix :: String -> Matrix Double
readMatrix = fromLists . map (map read). map words . filter (not.null) . lines

{- |  obtains the number of rows and columns in an ASCII data file
      (provisionally using unix's wc).
-}
fileDimensions :: FilePath -> IO (Int,Int)
fileDimensions fname = do
    wcres <- readProcess "wc" ["-w",fname] ""
    contents <- readFile fname
    let tot = read . head . words $ wcres
        c   = length . head . dropWhile null . map words . lines $ contents
    if tot > 0
        then return (tot `div` c, c)
        else return (0,0)

-- | Loads a matrix from an ASCII file formatted as a 2D table.
loadMatrix :: FilePath -> IO (Matrix Double)
loadMatrix file = fromFile file =<< fileDimensions file

-- | Loads a matrix from an ASCII file (the number of rows and columns must be known in advance).
fromFile :: FilePath -> (Int,Int) -> IO (Matrix Double)
fromFile filename (r,c) = reshape c `fmap` fscanfVector filename (r*c)

max_stars_repo_name: curiousleo/liquidhaskell

max_issues_repo_name: curiousleo/liquidhaskell

max_forks_repo_name: curiousleo/liquidhaskell


text:
module Lib.Practice.P6
    (
        f6
    ) where

import Lib.Utility.Mnist
import Control.Monad.Trans.Maybe
import Control.Monad.IO.Class
import Numeric.LinearAlgebra.Static
import qualified Numeric.LinearAlgebra as LA

f6' :: MaybeT IO ()
f6' = do
    ((trD, trL), _) <- getMnistData
    MaybeT $ fmap Just $ print $ size trD
    let (l, d) = (extract trL LA.! 0, extract trD LA.! 0)
    liftIO $ print d
    liftIO $ print l
    -- TODO 画像の表示

f6 :: IO ()
f6 = do
    let (MaybeT result) = f6'
    result' <- result
    case result' of
        Nothing -> print "失敗"
        Just x -> return x


text:
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE RecordWildCards #-}
{-# LANGUAGE RankNTypes #-}
-----------------------------------------------------------------------------
--
-- Module      :  AI.Network.RMM.LSTM
-- Copyright   :  (c) JP Moresmau
-- License     :  BSD3
--
-- Maintainer  :  JP Moresmau <jp@moresmau.fr>
-- Stability   :  experimental
-- Portability :
--
-- | Long Short Term Memory networks
--
-----------------------------------------------------------------------------

module AI.Network.RNN.LSTM where

import Control.DeepSeq

import Data.List

import Numeric.LinearAlgebra.HMatrix

import AI.Network.RNN.Types
import AI.Network.RNN.Util
import AI.Network.RNN.Expr
import Debug.Trace
import qualified Data.IntMap as I

import Numeric.AD

-- | The LSTM data type, with the different weights and states
data LSTMNetwork = LSTMNetwork
    { lstmSize :: !Int
    , lstmWeightsW :: !(Matrix Double)
    , lstmWeightsU :: !(Matrix Double)
    , lstmBias :: !(Vector Double)
    , lstmState :: !(Vector Double)
    , lstmOutput :: !(Vector Double)
    } deriving (Show,Read,Eq)

--instance Show LSTMNetwork where
--     show s = let
--        a= show $ (rnnsize s,toVector s)
--        in trace (show $ length a) a
--
--instance Read LSTMNetwork where
--    readsPrec d s =
--        let [(ok,r)] = readsPrec d s
--            dat = uncurry fromVector ok
--        in [(dat,r)]

-- | Force evaluation instance
instance NFData LSTMNetwork where
    rnf LSTMNetwork{..} = rnf (lstmSize,lstmWeightsW,lstmWeightsU,lstmBias,lstmState,lstmOutput)

-- | Implement the network evaluation and conversions functions
instance RNNEval LSTMNetwork where
    type Size LSTMNetwork = Int
    evalStep n@LSTMNetwork{..} is =
        let z = (lstmWeightsW #> is) + (lstmWeightsU #> lstmOutput) + lstmBias
            [i,f,c1,o] = takesV (replicate 4 lstmSize) z
            c2 = cmap tanh c1
            ns = (c2 * cmap sigmoid i) + (cmap sigmoid f * lstmState)
            no = cmap sigmoid o * cmap tanh ns
        in force (n{lstmState=ns,lstmOutput=no},no)
    fromVector sz vs =
        let
            msize = sz * sz
            [v1,v2,v3,v4,v5] = takesV [msize * 4,msize * 4, sz * 4,sz,sz] vs
            m1 = reshape sz v1
            m2 = reshape sz v2
        in LSTMNetwork sz m1 m2 v3 v4 v5
    toVector LSTMNetwork{..} = force $ vjoin
        [ flatten lstmWeightsW, flatten lstmWeightsU, lstmBias, lstmState, lstmOutput]
    rnnsize = lstmSize
    fullSize = lstmFullSize .lstmSize

-- | Full size of a network
lstmFullSize :: FullSize LSTMNetwork
lstmFullSize sz = (sz * sz) * 8 + sz* 4 + sz + sz

data LSTMIO = LSTMIO
    { lioInput :: Matrix Double
    , lioLstms :: [LSTMNetwork]
    , lioOutput :: Matrix Double
    , lioSer   :: Maybe (Vector Double)
    } deriving (Eq)

instance NFData LSTMIO where
  rnf LSTMIO{..} = rnf (lioInput,lioLstms,lioOutput)

instance RNNEval LSTMIO  where
    type Size LSTMIO = (Int,Int,Int,Int)
    evalStep io@LSTMIO{..} is =
        let inps = cmap sigmoid $ lioInput #> (vjoin [is,scalar 1])
            (lstm2,outs) = evalStep lioLstms inps
        in (io{lioLstms=lstm2,lioSer=Nothing},lioOutput #> outs)
    fromVector (is,lnnSize,lnnNumber,os) vs =
        let [v1,v2,v3] = takesV [(is+1)* lnnSize ,lnnNumber * (lstmFullSize lnnSize),os * lnnSize ] vs
            m1 = reshape (is+1) v1
            m2 = reshape lnnSize v3
        in LSTMIO m1 (fromVector (replicate lnnNumber lnnSize) v2) m2 $ Just vs
    toVector LSTMIO{..} = case lioSer of
        Just vs -> vs
        Nothing -> vjoin [flatten lioInput,toVector lioLstms,flatten lioOutput]
    rnnsize LSTMIO{..} = (cols lioInput - 1, rnnsize $ head lioLstms, length lioLstms,rows lioOutput)
    fullSize =lstmioFullSize . rnnsize
--        let (r1,c1)=size lioInput
--            (r2,c2)=size lioOutput
--        in (r1*c1) + (fullSize lioLstms) + (r2*c2)

instance Show LSTMIO where
    show l = show (rnnsize l,toVector l)

instance Read LSTMIO where
    readsPrec d s =
        let [(ok,r)] = readsPrec d s
            dat = uncurry fromVector ok
        in [(dat,r)]

lstmioFullSize :: FullSize LSTMIO
lstmioFullSize (is,lnnSize,lnnNumber,os) = (is+1) * lnnSize + lnnNumber * (lstmFullSize lnnSize) + os * lnnSize

data LSTMList = LSTMList
    { llSize :: (Int,Int,Int,Int)
    , llData :: [Double]
    } deriving (Show, Read, Eq)

instance NFData LSTMList where
  rnf LSTMList{..} = rnf (llSize,llData)

instance RNNEval LSTMList where
    type Size LSTMList = (Int,Int,Int,Int)
    rnnsize = llSize
    fullSize =lstmioFullSize . rnnsize
    toDList = llData
    fromDList sz = LSTMList sz
    toVector = fromList . toDList
    fromVector sz = fromDList sz . toList
    evalStep l@LSTMList{..} is =
        let lstmio :: LSTMIO = fromDList llSize llData
            (io2,out) = evalStep lstmio is
        in (l{llData=toDList io2},out)

-- | Implementation of the LSTM evaluation step without explicit matrices and vectors
-- just using lists, so we can use AD on it
lstmList :: (Num b,Floating b) => Int -> [b] -> [b] -> ([b],[b])
lstmList sz lstm is = let
    msize = sz * sz
    [mW,mU,vB,vS,vO] = takes [msize * 4,msize * 4, sz * 4,sz,sz] lstm
    z = zipWith3 (\a b c->a+b+c) (listMProd mW is) (listMProd mU vO) vB
    [i,f,c1,o] = takes (replicate 4 sz) z
    c2 = map tanh c1
    ns = zipWith (+) (zipWith (*) c2 (map sigmoid i)) (zipWith (*) (map sigmoid f) vS)
    no = zipWith (*) (map sigmoid o) (map tanh ns)
    in (mW++mU++vB++ns++no,no)

-- | Cost calculation using list representation for AD
cost' :: (Num b,Floating b,Fractional b) => Int -> [[b]] -> [[b]] -> [b] -> b
cost' sz is os lstm = let
    (_,res) = mapAccumL (lstmList sz) lstm is
--    in - (calcMeanList $ three (last res) (last os))
--    in - (calcMeanList $ concat $ zipWith (three) res os)
    in sum $ zipWith err os res
    where
      err :: (Num b,Floating b) => [b] -> [b] -> b
      err a b  = sum (zipWith (\c d -> (c- d)**2 ) a b)
--      one i o = zipWith (*) o (map log i)
--      two i o = zipWith (*) (map (\x->1 - x) o) (map (\x->1 - log x) i)
--      three i o= zipWith (+) (one i o) (two i o)

-- | Gradient descent learning
-- The third parameter is a call back function to monitor progress and stop the learning process if needed
learnGradientDescent :: (Monad m) => LSTMNetwork -> TrainData a -> (LSTMNetwork -> TrainData a -> Int -> m Bool) -> m LSTMNetwork
learnGradientDescent lstm td progressF =  go (toList $ toVector lstm) 0
    -- go2 gds 0
    where
      go2 (ls:lss) gen =  do
        let rnn::LSTMNetwork = fromVector (tdRecSize td) (fromList ls)
        cont <- progressF rnn td gen
        if cont
            then go2 lss (gen+1)
            else return rnn
      go ls gen = do
        let rnn::LSTMNetwork = fromVector (tdRecSize td) (fromList ls)
        cont <- progressF rnn td gen
        if cont
            then do
                let
                    gs= gf ls -- gradients using AD
                    ls2 = zipWith (\o g->o-g*0.1) ls gs
                go ls2 (gen+1)
            else return rnn
      lis = map toList (tdInputs td)
      los = map toList (tdOutputs td)
      gf = grad (cost' (tdRecSize td) (map (map auto) lis) (map (map auto) los))
      gds = gradientDescent (cost' (tdRecSize td) (map (map auto) lis) (map (map auto) los)) (toList $ toVector lstm)

-- | Gradient descent learning using symbolic differentiation
--   The goal was to calculate the derivative once and then just close and eval the result using the current data
--   However it is much slower than normal automatic differentiation
learnGradientDescentSym :: (Monad m) => LSTMNetwork -> TrainData a -> (LSTMNetwork -> TrainData a -> Int -> m Bool) -> m LSTMNetwork
learnGradientDescentSym lstm td progressF =  go (toList $ toVector lstm) 0
    where
      go ls gen = do
        let rnn::LSTMNetwork = fromVector (tdRecSize td) (fromList ls)
        cont <- progressF rnn td gen
        if cont
            then do
                let
                    i = I.fromList $ zip [0..] ls
                    cexpr = map (\g-> close g i) gf
                    gs = map eval cexpr
                    ls2 = zipWith (\o g->o-g*0.1) ls gs
                go ls2 (gen+1)
            else return rnn
      lis = map toList (tdInputs td)
      los = map toList (tdOutputs td)
      ls0 = toList $ toVector lstm
      gf = map fullSimplify $ grad (cost'
            (tdRecSize td)
            (map (map (\x -> autoEval (Lit $ Lit x) I.empty)) lis)
            (map (map (\x -> autoEval (Lit $ Lit x) I.empty)) los))
            (zipWith (\_ i->Var i) ls0 [0..])

-- | RMSProp learning, as far as I can make out
-- The third parameter is a call back function to monitor progress and stop the learning process if needed
learnRMSProp :: (Monad m) => LSTMNetwork -> TrainData a -> (LSTMNetwork -> TrainData a -> Int -> m Bool) -> m LSTMNetwork
learnRMSProp lstm td progressF = go ls0 (replicate myl 0) (replicate myl 0) (replicate myl 0) 0
    where
      go ls rgs rgs2 ugs gen = do
        let rnn::LSTMNetwork = fromVector (tdRecSize td) (fromList ls)
        cont <- progressF rnn td gen
        if cont
            then do
                let
                    gs= gf ls -- gradients using AD
                    rgup = zipWith (\rg g-> 0.95 * rg + 0.05 * g) rgs ls
                    rg2up = zipWith (\rg2 g-> 0.95 * rg2 + 0.05 * (g ** 2)) rgs2 ls
                    ugup = zipWith4 (\ud zg rg rg2 -> 0.9 * ud - 1e-4 * zg / sqrt(rg2 - rg ** 2 + 1e-4)) ugs gs rgup rg2up
                    ls2 = zipWith (+) ls ugup
                go (force ls2) (force rgup) (force rg2up) (force ugup) (gen+1)
            else return rnn
      lis = map toList (tdInputs td)
      los = map toList (tdOutputs td)
      gf = grad (cost' (tdRecSize td) (map (map auto) lis) (map (map auto) los))
      ls0 = toList $ toVector lstm
      myl= length ls0

max_stars_repo_name: JPMoresmau/rnn

max_issues_repo_name: JPMoresmau/rnn

max_forks_repo_name: JPMoresmau/rnn


text:
module ComplexUtils (toComplex, toComplexM) where

import Data.Complex (Complex((:+)))

toComplex :: Double -> (Complex Double)
toComplex a = (a :+ 0)

toComplexM :: [Double] -> [(Complex Double)]
toComplexM = map toComplex

max_stars_repo_name: vined/fourier

max_issues_repo_name: vined/fourier

max_forks_repo_name: vined/fourier


text:
-----------------------------------------------------------------------------
-- |
-- Module      :  Numeric.Statistics.Histogram
-- Copyright   :  (c) A. V. H. McPhail 2010, 2014
-- License     :  BSD3
--
-- Maintainer  :  haskell.vivian.mcphail <at> gmail <dot> com
-- Stability   :  provisional
-- Portability :  portable
--
-- create histograms from density functions
--
-----------------------------------------------------------------------------

module Numeric.Statistics.Histogram (
                                     cumulativeToHistogram
                                    , gaussianHistogram
                                  ) where

-----------------------------------------------------------------------------

--import qualified Data.Packed.Vector as V
import qualified Data.Vector.Storable as V

import qualified Numeric.GSL.Histogram as H
--import qualified Numeric.GSL.Histogram2D as H2

import qualified Numeric.GSL.Distribution.Continuous as C

--import Numeric.LinearAlgebra.Algorithms
--import Numeric.LinearAlgebra.Interface()

-----------------------------------------------------------------------------

vectorToTuples = toTuples . V.toList
    where toTuples []         = error "need a minimum of two elements"
          toTuples [_]        = error "need a minimum of two elements"
          toTuples [x1,x2]    = [(x1,x2)]
          toTuples (x1:x2:xs) = (x1,x2) : (toTuples (x2:xs))

-----------------------------------------------------------------------------

cumulativeToHistogram :: (Double -> Double)   -- ^ the cumulative distribution function D(x <= X)
                   -> V.Vector Double        -- ^ the bins
                   -> H.Histogram            -- ^ the resulting histogram
cumulativeToHistogram f v = H.addListWeighted (H.emptyRanges v) $ map (\(x1,x2) -> ((x1 + x2) / 2.0,f x2 - f x1)) (vectorToTuples v)

gaussianHistogram :: Double               -- ^ mean
                  -> Double              -- ^ standard deviation
                  -> V.Vector Double     -- ^ the bins
                  -> H.Histogram         -- ^ the resulting histogram
gaussianHistogram u s = cumulativeToHistogram (\x -> C.density_1p C.Gaussian C.Lower s (x-u))

-----------------------------------------------------------------------------

max_stars_repo_name: amcphail/hstatistics

max_issues_repo_name: amcphail/hstatistics

max_forks_repo_name: amcphail/hstatistics


text:
{-# LANGUAGE AllowAmbiguousTypes       #-}
{-# LANGUAGE DataKinds                 #-}
{-# LANGUAGE FlexibleContexts          #-}
{-# LANGUAGE GADTs                     #-}
{-# LANGUAGE NoMonomorphismRestriction #-}
{-# LANGUAGE OverloadedStrings         #-}
{-# LANGUAGE ScopedTypeVariables       #-}
{-# LANGUAGE TypeApplications          #-}
{-# LANGUAGE TypeOperators             #-}
{-# LANGUAGE QuasiQuotes               #-}
{-# LANGUAGE PartialTypeSignatures #-}
module Main where

import qualified Control.Foldl                as FL
import           Control.Monad.IO.Class (MonadIO (..))
import qualified Data.List                    as List
import           Data.Maybe                   (fromMaybe)
import qualified Data.Map                     as M
import qualified Data.Text                    as T
import qualified Data.Text.IO                 as T
import qualified Data.Text.Lazy               as TL
import qualified Data.Vector.Storable         as V
import qualified Data.Vinyl                   as V
import qualified Frames                       as F

import qualified Numeric.LinearAlgebra        as LA
import           Numeric.LinearAlgebra        (R, Matrix)

import qualified Frames.Regression            as FR
import qualified Frames.VegaLite               as FV

import qualified Text.Blaze.Html5              as H
import qualified Text.Blaze.Html5.Attributes   as HA
import           Text.Blaze.Html              ((!))
import qualified Text.Blaze.Html.Renderer.Text as BH

import qualified Statistics.Types             as S

import qualified Knit.Report                  as K
import qualified Knit.Report.Other.Blaze      as KB

{-
import qualified Polysemy                      as P
import qualified Knit.Effects.Logger           as Log
import qualified Knit.Effects.Pandoc           as PE
import qualified Knit.Effects.PandocMonad           as PM
import qualified Knit.Report.Blaze             as RB
import qualified Knit.Report.Pandoc           as PR
-}

import Data.String.Here

templateVars = M.fromList
  [
    ("lang", "English")
  , ("author", "Adam Conner-Sax")
  , ("pagetitle", "Frame Regression Examples")
--  , ("tufte","True")
  ]

regressionNotesMD
  = [here|
## Regression Algorithm Comparison
* To get data for testing we choose xs evenly spaced between -0.5 and 0.5.  Then we compute ys via $y=1 + 2.2x$.  Then we add Gaussian noise to the ys.  We can optionally add Gaussian noise to the observed xs as well (that is, the ys are still calculated from the exact xs but we add noise to the xs before we regress) and add weights to the (x,y) pairs.
* For unweighted data we use ordinary least squares (OLS) and total least square (TLS).
* For weighted data we use weighted least squares (WOLS) and weighted total least squares (WTLS).
* Both TLS and WTLS are computed via the Singular Value Decomposition.
* For each comparison we vary the level of noise on the ys and/or xs, regress, then show the results of each regression in tabular form, a plot of the fits and prediction intervals together with a scatter of the noisy data, and a plot of the coefficients and confidence intervals. The shaded regions in each plot are the "prediction intervals" which take into account the uncertainty of the coefficients as well as the remainining noise in the data.
* The F-stat is computed for each fit via comparision to the model of intercept-only.  So the p-value of the overall fit is the probability that the data are better explained as noise around their (weighted) mean value than by the fit.
|]

main :: IO ()
main = asPandoc

asPandoc :: IO ()
asPandoc = do
  let pandocWriterConfig = K.PandocWriterConfig (Just "pandoc-templates/minWithVega-pandoc.html")  templateVars K.mindocOptionsF
  htmlAsTextE <- K.knitHtml (Just "FrameRegressions.Main") K.logAll pandocWriterConfig $ do
    K.addMarkDown regressionNotesMD
    testMany
  case htmlAsTextE of
    Right htmlAsText -> T.writeFile "examples/html/FrameRegressions.html" $ TL.toStrict  $ htmlAsText
    Left err -> putStrLn $ "pandoc error: " ++ show err

{-
asBlaze :: IO ()
asBlaze = do
  let runAllP = FR.runPandocAndLoggingToIO Log.logAll . Log.wrapPrefix "Main" . blazeToText
  htmlAsTextE <- runAllP $ do
    regressionNotesBlaze <- P.markDownTextToBlazeFragment regressionNotesMD
    blaze $ H.makeReportHtml "Frame Regression Examples" $ do
      H.placeTextSection $ regressionNotesBlaze
    testMany
  case htmlAsTextE of
    Right htmlAsText -> T.writeFile "examples/html/FrameRegressions.html" $ TL.toStrict  $ htmlAsText
    Left err -> putStrLn $ "pandoc error: " ++ show err
-}

-- regression tests

-- build some data for testing
-- uniformly distributed measurements, normally distributed noise.  Separate noise amplitudes for xs and ys
-- also allow building heteroscedastic data

type Y = "y" F.:-> Double
type X = "x" F.:-> Double
--type X2 = "x2" F.:-> Double
type Weight = "weight" F.:-> Double

buildRegressable :: [Double] -> Maybe (LA.Vector R) -> Double -> LA.Vector R -> Double -> IO (LA.Vector R, LA.Matrix R)
buildRegressable variances offsetsM noiseObs coeffs noiseMeas = do
  -- generate random measurements
  let d = LA.size coeffs
      nObs = List.length variances
      xsO = LA.asColumn (LA.fromList (List.replicate nObs 1)) LA.<> LA.asRow (fromMaybe (LA.fromList $ List.replicate (d-1) 0) offsetsM)
--  xs0 <- LA.cmap (\x-> x - 0.5) <$> LA.rand nObs (d-1) -- 0 centered uniformly distributed random numbers
  let xs0 :: Matrix R = LA.asColumn $ LA.fromList $ [-0.5 + (realToFrac i/realToFrac nObs) | i <- [0..(nObs-1)]]
  xNoise <- LA.randn nObs (d-1) -- 0 centered normally (sigma=1) distributed random numbers
  let xsC = 1 LA.||| (xs0 + xsO)
      xsN = 1 LA.||| (xs0 + xsO + LA.scale noiseMeas xNoise)
  let ys0 = xsC LA.<> LA.asColumn coeffs
  yNoise <- fmap (List.head . LA.toColumns) (LA.randn nObs 1) -- 0 centered normally (sigma=1) distributed random numbers
  let ys = ys0 + LA.asColumn (LA.scale noiseObs (V.zipWith (*) yNoise (LA.cmap sqrt (LA.fromList variances))))
  return (List.head (LA.toColumns ys), xsN)

type AllCols = [Y,X,Weight]

makeFrame :: [Double] -> (LA.Vector R, LA.Matrix R) -> IO (F.FrameRec AllCols)
makeFrame vars (ys, xs) = do
  if snd (LA.size xs) /= 2 then error ("Matrix of xs wrong size (" ++ show (LA.size xs) ++ ") for makeFrame") else return ()
  let rows = LA.size ys
      rowIndex = [0..(rows - 1)]
      makeRecord :: Int -> F.Record AllCols
      makeRecord n = ys `LA.atIndex` n F.&: xs `LA.atIndex` (n,1) F.&: vars List.!! n F.&: V.RNil -- skip bias column
  return $ F.toFrame $ makeRecord <$> rowIndex

unweighted :: Int -> [Double]
unweighted n = List.replicate n (1.0)

coneWeighted :: Int -> Double -> Double -> [Double]
coneWeighted n base increment =
  let w0 = [base + (realToFrac i) * increment | i <- [1..n]]
      s = realToFrac n/FL.fold FL.sum w0
  in fmap (*s) w0

showText :: Show a => a -> T.Text
showText = T.pack . show

coeffs :: LA.Vector R = LA.fromList [1.0, 2.2]
offsets :: LA.Vector R = LA.fromList [1]
vars = coneWeighted 100 1 0.1
varListToWeights = LA.cmap (\x -> 1/sqrt x) . LA.fromList
wgts = varListToWeights vars

testRegressions :: ( F.ColumnHeaders '[w]
                   , FR.BoolVal (FR.NonVoidField w)
                   , V.KnownField w
                   , Traversable f
                   , Foldable f
                   , K.Member K.ToPandoc effs
                   , K.PandocEffects effs
                   , MonadIO (K.Sem effs))
                => Double
                -> Double
                -> Bool
                -> Bool
                -> T.Text
                -> f (T.Text, (F.FrameRec AllCols -> K.Sem effs (FR.FrameRegressionResult Y True '[X] w AllCols)))
                -> K.Sem effs ()
testRegressions  yNoise xNoise weighted offset vizId keyedFs = do
  let title = "Sy=" <> showText yNoise <> " gaussian noise added to ys & "
               <> "Sx=" <> showText xNoise <> " gaussian noise added to xs"
               <> " (" <> (if weighted then "cone weights" else "unweighted") <> (if offset then ", w/offsets" else "") <> ")"
      vars = if weighted then coneWeighted 100 1 0.1 else unweighted 100
      wgts = varListToWeights vars
      offsetM = if offset then (Just offsets) else Nothing
      doOne dat (key, f) = do
        result <- f dat
        return (key, result)
  frame <- liftIO (buildRegressable vars offsetM yNoise coeffs xNoise >>= makeFrame vars)
  results <- traverse (doOne frame) keyedFs
  let header _ _ = title
  K.addMarkDown $ "\n## " <> title
  K.addBlaze $ do
    H.div ! HA.style "display: block-inline" $ do
      FR.prettyPrintRegressionResults id results S.cl95 FR.prettyPrintRegressionResultBlaze mempty
      KB.placeVisualization (vizId <> "_fits") $ FV.keyedLayeredFrameScatterWithFit title id results S.cl95 frame
      KB.placeVisualization (vizId <> "_regresssionCoeffs") $ FV.regressionCoefficientPlotMany id "Parameters" ["intercept","x"] (fmap (\(k,frr) -> (k, FR.regressionResult frr)) results) S.cl95

-- I can't test the weighted and unweighted on the same things because those algos return different types in their results.  Which, maybe, is a good point.
testMany :: ( K.Member K.ToPandoc effs
            , K.PandocEffects effs
            , MonadIO (K.Sem effs)) =>  K.Sem effs ()
testMany = K.wrapPrefix "Many" $ do
  let toTestUW :: _ -- that this is required is suspicious
      toTestUW =
        [
          ("OLS", FR.ordinaryLeastSquares)
        , ("TLS", FR.totalLeastSquares)
        ]
      toTestW :: _ -- that this is required is suspicious
      toTestW =
        [
          ("WOLS", FR.weightedLeastSquares)
        , ("WTLS", FR.weightedTLS)
        ]
  testRegressions 0.0 0.0 False False "many1" toTestUW
  testRegressions 0.3 0.0 False False "many2" toTestUW
  testRegressions 0.5 0.0 False False "many3" toTestUW
  testRegressions 0.3 0.1 False False "many4" toTestUW
  testRegressions 3 0.0 False False "many9" toTestUW
  testRegressions @Weight 0.3 0.0 True False "many6" toTestW
  testRegressions @Weight 0.5 0.0 True False "many7" toTestW
  testRegressions @Weight 0.3 0.1 True False "many8" toTestW

max_stars_repo_name: teto/Frames-utils

max_issues_repo_name: teto/Frames-utils

max_forks_repo_name: teto/Frames-utils


text:
module BondCalculus.StochPyExtraction
  (stochPyModel, generateStochPy, reaction, ReactionVect, reactions, ReactionSystem(..), extractReactionSystem, fromIVPToStochPyModel, simulateStochPy) where

import BondCalculus.Symbolic
import BondCalculus.Vector
import BondCalculus.AST
import BondCalculus.Base
import qualified BondCalculus.AST as AST
import BondCalculus.Transitions (MTS)
import Data.Bifunctor
import BondCalculus.Simulation (Trace)
-- import qualified Data.HashMap.Strict as H
-- import BondCalculus.Simulation (Trace)
-- import Data.String.Utils
import BondCalculus.ODEExtraction (sympyExpr, IVP(..), sympySimplify, runPython)
-- , runPython)

-- import qualified Control.Exception as X
import qualified Data.Map as M
import qualified Data.List as L
import Data.Maybe
import BondCalculus.Processes

--import qualified Numeric.LinearAlgebra as LA

-- import System.IO.Unsafe
-- import qualified System.Process as OS
-- import System.Process (readProcess)
-- import Data.Maybe

type ReactionVect k = Vect (Tensor (ProcessVect Conc) (ProcessVect Conc)) k

reaction :: (Vector k (DirectionVect k), Vector k (ReactionVect k), Nullable (DirectionVect k), Expression k) =>
            [DirectionVect k] -> ReactionVect k
reaction = multilinear react' <$> filter (/=vectZero)
  where react' xs = vect( sourceVect xs :* targetVect xs )
        source (spec :* _) = spec
        targetVect xs = embed $ concretify $ foldl (<|>) (mkAbsBase Nil) (map target xs)
        target (_ :* spec') = spec'
        sourceVect xs = foldl (+>) vectZero (map (embed.source) xs)

reactions :: (Vector k (ProcessVect k), Show k, DoubleExpression k, ExpressionOver Double k) => (Species -> MTS) -> ConcreteAffinityNetwork Double -> ProcessVect k -> ReactionVect k
reactions tr network p = simplify $ dPdt'' reaction tr network p

type ConcVect = Vect String Conc
newtype ReactionSystem = ReactionSystem ([String], [(SymbolicExpr Double, ConcVect, ConcVect)], [Double])
  deriving (Show, Eq, Ord)

extractReactionSystem :: AST.Env -> ConcreteAffinityNetwork Double -> P' Double -> [Double] -> ReactionSystem
extractReactionSystem env network p inits = ReactionSystem (vars, transitions, inits)
  where tr = tracesGivenNetwork network env
        vars = map (pretty.simplify.snd) (toList p)
        reacts = reactions tr network p
        transitions = [ (simplify coeff, toConcVect u, toConcVect v)
                      | (coeff, u :* v) <- toList reacts]
        toConcVect  = fromList . map (second (pretty.simplify)) . toList

fromIVPToStochPyModel :: IVP Double -> Either String String
fromIVPToStochPyModel ivp =
  if any isNothing exprs
  then Left "An ODE equation has an unbound variable"
  else Right $
  "Output_In_Conc: True\n" ++
  "Species_In_Conc: True\n\n" ++
  L.intercalate "\n" rules ++ "\n" ++
  L.intercalate "\n" initlines
    where
      species = ["S" ++ show i | i <- [(0::Integer)..]]
      varmap = M.fromList $ zip vars species
      IVP (vars,rhss,inits) = ivp
      mkgrowrule spec y =
        "grow" ++ spec ++ ":\n" ++
        "\t$pool > " ++ spec ++ "\n" ++
        "\t" ++ mkgrowrate y ++ "\n"
      mkgrowrate y = "abs(" ++ y ++ " + abs(" ++ y ++ "))/2"
      mkdecayrule spec y =
        "decay" ++ spec ++ ":\n" ++
        "\t" ++ spec ++ " > $pool\n" ++
        "\t" ++ mkdecayrate y ++ "\n"
      mkdecayrate y = "abs(-" ++ y ++ " + abs(" ++ y ++ "))/2"
      exprs = map (sympyExpr varmap) rhss
      growrules = zipWith mkgrowrule species (catMaybes exprs)
      decayrules = zipWith mkdecayrule species (catMaybes exprs)
      rules = zipWith (\x y -> x ++ "\n" ++ y) growrules decayrules
      initlines = [y ++ " = " ++ show y0 | (y,y0) <- zip species inits]

stochPyModel :: ReactionSystem -> Double -> Either String String
stochPyModel (ReactionSystem (vars, transitions, inits)) h =
  if any isNothing transitions'
  then Left "An ODE equation has an unbound variable"
  else Right $
  "Output_In_Conc: True\n" ++
  "Species_In_Conc: True\n\n" ++
  L.intercalate "\n" rules ++ "\n" ++
  "h = " ++ show h ++ "\n" ++
  L.intercalate "\n" initlines
    where
      species = ["S" ++ show i | i <- [(0::Integer)..]]
      reacts = ["R" ++ show i | i <- [(0::Integer)..]]
      varmap = M.fromList $ zip vars species
      varmaph = fmap (\s -> "(" ++ s ++ "*h)") varmap
      rules = [
            i ++ ":\n" ++
            "\t" ++ formatVector u ++ " > " ++ formatVector v ++ "\n" ++
            "\t" ++ rate ++ "\n"
          | (i, (rate, u, v)) <- reacts `zip` catMaybes transitions'
        ]
      transitions' = fmap (\(rate,u,v) -> do
                      rate' <- sympyExpr varmaph rate
                      u' <- vectorVarLookup u
                      v' <- vectorVarLookup v
                      return (sympySimplify ("(" ++ rate' ++ ")/h"),u',v')) transitions
      vectorVarLookup u =
        let us = map (second (`M.lookup` varmap)) $ toList u
        in if any (isNothing.snd) us
          then Nothing
          else Just (fromList [(k,v) | (k, Just v) <- us])
      formatVector u = if null us
                       then "$pool"
                       else L.intercalate " + "
                              ["{" ++ show k ++ "} " ++ i | (k,i) <- us]
        where us = toList u
      initlines = [y ++ " = " ++ show (y0/h) | (y,y0) <- zip species inits]

generateStochPy :: String -> AST.Env -> ConcreteAffinityNetwork Double -> P' Double -> Double -> [Double] -> IO String
generateStochPy filename env network p h inits = case stochPyModel (extractReactionSystem env network p inits) h of
  Right script -> do
    -- putStrLn $ "Python script:\n\n" ++ script
    writeFile filename script
    return script
  Left _ -> undefined

simulateStochPy :: String -> Int -> AST.Env -> ConcreteAffinityNetwork Double -> P' Double -> Double -> [Double] -> IO Trace
simulateStochPy method steps env network p h inits = case stochPyModel (extractReactionSystem env network p inits) h of
  Right stochpy -> do
    -- putStrLn $ "Python script:\n\n" ++ script
    writeFile "/tmp/stochpy.psc" stochpy
    let script = "import os, sys\n" ++
                 "f = open(os.devnull, 'w')\n" ++
                 "stdout = sys.stdout\n" ++
                 "sys.stdout = f\n" ++
                 "import stochpy\n" ++
                 "smod = stochpy.SSA()\n" ++
                 "smod.Model('/tmp/stochpy.psc')\n" ++
                 "smod.DoStochSim(method='" ++ method ++ "', end=" ++
                   show steps ++ ", mode='time')\n" ++
                 "sys.stdout=stdout\n" ++
                 "xss = smod.data_stochsim.getSimData(" ++ L.intercalate "," ["'S" ++ show i ++ "'" | (i,_) <- [0..] `zip` kis] ++ ")\n" ++
                 "xss[:,1:] *= " ++ show h ++ "\n" ++
                 "for xs in xss:\n" ++
                 "    print(' '.join(map('{:.18f}'.format, list(xs))))\n"
    writeFile "stochpy_script.py" script
    res <- runPython script
    let yss = (map (map read.words) $ lines res) :: [[Double]]
    return [ (head ys, toP (tail ys)) | ys <- yss ]
  Left _ -> undefined
  where toP ks = fromList [(k,i) | (k, (_,i)) <- ks `zip` kis]
        kis = toList p

max_stars_repo_name: twright/bondwb

max_issues_repo_name: twright/bondwb

max_forks_repo_name: twright/bondwb


text:
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE DeriveGeneric       #-}
{-# LANGUAGE GADTs               #-}
{-# LANGUAGE KindSignatures      #-}
{-# LANGUAGE RankNTypes          #-}
{-# LANGUAGE RecordWildCards     #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE StandaloneDeriving  #-}
{-# LANGUAGE TypeApplications    #-}
{-# LANGUAGE TypeInType          #-}
{-# LANGUAGE TypeOperators       #-}

-- |
-- Module      : Numeric.Hamilton
-- Description : Hamiltonian dynamics for physical systems on generalized
--               coordinates using automatic differentiation
-- Copyright   : (c) Justin Le 2016
-- License     : BSD-3
-- Maintainer  : justin@jle.im
-- Stability   : unstable
-- Portability : portable
--
-- Simulate physical systems on generalized/arbitrary coordinates using
-- Hamiltonian mechanics and automatic differentiation!
--
-- See the <https://github.com/mstksg/hamilton#readme> for more
-- information on usage!
--

module Numeric.Hamilton
  ( -- * Systems and states
    -- ** Systems
    System
  , mkSystem
  , mkSystem'
  , underlyingPos
    -- ** States
  , Config(..)
  , Phase(..)
  , toPhase
  , fromPhase
    -- * State functions
  , momenta
  , velocities
  , keC
  , keP
  , pe
  , lagrangian
  , hamiltonian
  , hamEqs
    -- * Simulating hamiltonian dynamics
    -- ** Over phase space
  , stepHam
  , evolveHam
  , evolveHam'
    -- ** Over configuration space
    -- | Convenience wrappers over the normal phase-space
    -- steppers/simulators that allow you to provide input and expect
    -- output in configuration space instead of in phase space.  Note that
    -- the simulation itself still runs in phase space, so these all
    -- require conversions to and from phase space under the hood.
  , stepHamC
  , evolveHamC
  , evolveHamC'
  ) where

import           Control.Monad
import           Data.Bifunctor
import           Data.Foldable
import           Data.Kind
import           Data.Maybe
import           Data.Proxy
import           Data.Type.Equality
import           GHC.Generics                        (Generic)
import           GHC.TypeLits
import           GHC.TypeLits.Compare
import           Numeric.AD
import           Numeric.GSL.ODE
import           Numeric.LinearAlgebra.Static        as H
import           Numeric.LinearAlgebra.Static.Vector
import qualified Data.Vector.Generic.Sized           as VG
import qualified Data.Vector.Sized                   as V
import qualified Numeric.LinearAlgebra               as LA

-- | Represents the full state of a system of @n@ generalized coordinates
-- in configuration space (informally, "positions and velocities")
--
-- A configuration space representaiton is more directly "physically
-- meaningful" and intuitive/understandable to humans than a phase space
-- representation.  However, it's much less mathematically ideal to work
-- with because of the lack of some neat underlying symmetries.
--
-- You can convert a @'Config' n@ into a @'Phase' n@ (convert from
-- configuration space to phase space) for a given system with 'toPhase'.
-- This allows you to state your system in configuration space and then
-- convert it to phase space before handing it off to the hamiltonian
-- machinery.
data Config :: Nat -> Type where
    Cfg :: { -- | The current values ("positions") of each of the @n@
             -- generalized coordinates
             cfgPositions :: !(R n)
             -- | The current rate of changes ("velocities") of each of the
             -- @n@ generalized coordinates
           , cfgVelocities :: !(R n)
           }
        -> Config n
  deriving (Generic)

deriving instance KnownNat n => Show (Config n)

-- | Represents the full state of a system of @n@ generalized coordinates
-- in phase space (informally, "positions and momentums").
--
-- Phase space representations are much nicer to work with mathematically
-- because of some neat underlying symmetries.  For one, positions and
-- momentums are "interchangeable" in a system; if you swap every
-- coordinate's positions with their momentums, and also swap them in the
-- equations of motions, you get the same system back.  This isn't the case
-- with configuration space representations.
--
-- A hamiltonian simulation basically describes the trajectory of each
-- coordinate through phase space, so this is the /state/ of the
-- simulation.  However, configuration space representations are much more
-- understandable to humans, so it might be useful to give an initial state
-- in configuration space using 'Config', and then convert it to a 'Phase'
-- with 'toPhase'.
data Phase :: Nat -> Type where
    Phs :: { -- | The current values ("positions") of each of the @n@
             -- generalized coordinates.
             phsPositions :: !(R n)
             -- | The current conjugate momenta ("momentums") to each of
             -- the @n@ generalized coordinates
           , phsMomenta :: !(R n)
           }
        -> Phase n
  deriving (Generic)

deriving instance KnownNat n => Show (Phase n)

-- | Represents a physical system in which physics happens.  A @'System'
-- m n@ is a system whose state described using @n@ generalized coordinates
-- (an "@n@-dimensional" system), where the underlying cartesian coordinate
-- space is @m@-dimensional.
--
-- For the most part, you are supposed to be able to ignore @m@.  @m@ is
-- only provided because it's useful when plotting/drawing the system with
-- a given state back in rectangular coordinates. (The only function that
-- use the @m@ at the moment is 'underlyingPos')
--
-- A @'System' m n@'s state is described using a @'Config' n@ (which
-- describes the system in configuration space) or a @'Phase' n@ (which
-- describes the system in phase space).
data System :: Nat -> Nat -> Type where
    Sys :: { _sysInertia       :: R m
           , _sysCoords        :: R n -> R m
           , _sysJacobian      :: R n -> L m n
           , _sysHessian       :: R n -> V.Vector n (L m n)
           , _sysPotential     :: R n -> Double
           , _sysPotentialGrad :: R n -> R n
           }
        -> System m n

-- | Converts the position of generalized coordinates of a system to the
-- coordinates of the system's underlying cartesian coordinate system.
-- Useful for plotting/drawing the system in cartesian space.
underlyingPos
    :: System m n
    -> R n
    -> R m
underlyingPos = _sysCoords

-- | The potential energy of a system, given the position in the
-- generalized coordinates of the system.
pe  :: System m n
    -> R n
    -> Double
pe = _sysPotential

vec2l
    :: (KnownNat m, KnownNat n)
    => V.Vector m (V.Vector n Double)
    -> L m n
vec2l = rowsL . fmap gvecR

-- | Create a system with @n@ generalized coordinates by describing its
-- coordinate space (by a function from the generalized coordinates to the
-- underlying cartesian coordinates), the inertia of each of those
-- underlying coordinates, and the pontential energy function.
--
-- The potential energy function is expressed in terms of the genearlized
-- coordinate space's positions.
mkSystem
    :: forall m n. (KnownNat m, KnownNat n)
    => R m      -- ^ The "inertia" of each of the @m@ coordinates
                -- in the underlying cartesian space of the system.  This
                -- should be mass for linear coordinates and rotational
                -- inertia for angular coordinates.
    -> (forall a. RealFloat a => V.Vector n a -> V.Vector m a)
                -- ^ Conversion function to convert points in the
                -- generalized coordinate space to the underlying cartesian
                -- space of the system.
    -> (forall a. RealFloat a => V.Vector n a -> a)
                -- ^ The potential energy of the system as a function of
                -- the generalized coordinate space's positions.
    -> System m n
mkSystem m f u = Sys
    { _sysInertia       =                     m
    , _sysCoords        = gvecR             . f           . grVec
    , _sysJacobian      = tr   . vec2l      . jacobianT f . grVec
    , _sysHessian       = tr2  . fmap vec2l . hessianF f  . grVec
    , _sysPotential     =                     u           . grVec
    , _sysPotentialGrad = gvecR             . grad u      . grVec
    }
  where
    tr2 :: forall o. (KnownNat n, KnownNat o)
        => V.Vector m (L n o)
        -> V.Vector n (L m o)
    tr2 = fmap rowsL . traverse lRows
    {-# INLINE tr2 #-}

-- | Convenience wrapper over 'mkSystem' that allows you to specify the
-- potential energy function in terms of the underlying cartesian
-- coordinate space.
mkSystem'
    :: forall m n. (KnownNat m, KnownNat n)
    => R m      -- ^ The "inertia" of each of the @m@ coordinates
                -- in the underlying cartesian space of the system.  This
                -- should be mass for linear coordinates and rotational
                -- inertia for angular coordinates.
    -> (forall a. RealFloat a => V.Vector n a -> V.Vector m a)
                -- ^ Conversion function to convert points in the
                -- generalized coordinate space to the underlying cartesian
                -- space of the system.
    -> (forall a. RealFloat a => V.Vector m a -> a)
                -- ^ The potential energy of the system as a function of
                -- the underlying cartesian coordinate space's positions.
    -> System m n
mkSystem' m f u = mkSystem m f (u . f)

-- | Compute the generalized momenta conjugate to each generalized
-- coordinate of a system by giving the configuration-space state of the
-- system.
--
-- Note that getting the momenta from a @'Phase' n@ involves just using
-- 'phsMomenta'.
momenta
    :: (KnownNat m, KnownNat n)
    => System m n
    -> Config n
    -> R n
momenta Sys{..} Cfg{..} = tr j #> diag _sysInertia #> j #> cfgVelocities
  where
    j = _sysJacobian cfgPositions

-- | Convert a configuration-space representaiton of the state of the
-- system to a phase-space representation.
--
-- Useful because the hamiltonian simulations use 'Phase' as its working
-- state, but 'Config' is a much more human-understandable and intuitive
-- representation.  This allows you to state your starting state in
-- configuration space and convert to phase space for your simulation to
-- use.
toPhase
    :: (KnownNat m, KnownNat n)
    => System m n
    -> Config n
    -> Phase n
toPhase s = Phs <$> cfgPositions <*> momenta s

-- | The kinetic energy of a system, given the system's state in
-- configuration space.
keC :: (KnownNat m, KnownNat n)
    => System m n
    -> Config n
    -> Double
keC s = do
    vs <- cfgVelocities
    ps <- momenta s
    return $ (vs <.> ps) / 2

-- | The Lagrangian of a system (the difference between the kinetic energy
-- and the potential energy), given the system's state in configuration
-- space.
lagrangian
    :: (KnownNat m, KnownNat n)
    => System m n
    -> Config n
    -> Double
lagrangian s = do
    t <- keC s
    u <- pe s . cfgPositions
    return (t - u)

-- | Compute the rate of change of each generalized coordinate by giving
-- the state of the system in phase space.
--
-- Note that getting the velocities from a @'Config' n@ involves just using
-- 'cfgVelocities'.
velocities
    :: (KnownNat m, KnownNat n)
    => System m n
    -> Phase n
    -> R n
velocities Sys{..} Phs{..} = inv jmj #> phsMomenta
  where
    j   = _sysJacobian phsPositions
    jmj = tr j H.<> diag _sysInertia H.<> j

-- | Invert 'toPhase' and convert a description of a system's state in
-- phase space to a description of the system's state in configuration
-- space.
--
-- Possibly useful for showing the phase space representation of a system's
-- state in a more human-readable/human-understandable way.
fromPhase
    :: (KnownNat m, KnownNat n)
    => System m n
    -> Phase n
    -> Config n
fromPhase s = Cfg <$> phsPositions <*> velocities s

-- | The kinetic energy of a system, given the system's state in
-- phase space.
keP :: (KnownNat m, KnownNat n)
    => System m n
    -> Phase n
    -> Double
keP s = do
    ps <- phsMomenta
    vs <- velocities s
    return $ (vs <.> ps) / 2

-- | The Hamiltonian of a system (the sum of kinetic energy and the
-- potential energy), given the system's state in phase space.
hamiltonian
    :: (KnownNat m, KnownNat n)
    => System m n
    -> Phase n
    -> Double
hamiltonian s = do
    t <- keP s
    u <- pe s . phsPositions
    return (t + u)

-- | The "hamiltonian equations" for a given system at a given state in
-- phase space.  Returns the rate of change of the positions and
-- conjugate momenta, which can be used to progress the simulation through
-- time.
--
-- Computed using the maths derived in
-- <https://blog.jle.im/entry/hamiltonian-dynamics-in-haskell.html>.
hamEqs
    :: (KnownNat m, KnownNat n)
    => System m n
    -> Phase n
    -> (R n, R n)
hamEqs Sys{..} Phs{..} = (dHdp, -dHdq)
  where
    mm   = diag _sysInertia
    j    = _sysJacobian phsPositions
    trj  = tr j
    jmj  = trj H.<> mm H.<> j
    ijmj = inv jmj
    dTdq = gvecR
         . flip fmap (_sysHessian phsPositions) $ \djdq ->
             -phsMomenta <.> ijmj #> trj #> mm #> djdq #> ijmj #> phsMomenta
    dHdp = ijmj #> phsMomenta
    dHdq = dTdq + _sysPotentialGrad phsPositions

-- | Step a system through phase space over over a single timestep.
stepHam
    :: forall m n. (KnownNat m, KnownNat n)
    => Double           -- ^ timestep to step through
    -> System m n       -- ^ system to simulate
    -> Phase n          -- ^ initial state, in phase space
    -> Phase n
stepHam r s p = evolveHam @m @n @2 s p (V.fromTuple (0, r))
                  `V.unsafeIndex` 1

-- | Evolve a system using a hamiltonian stepper, with the given initial
-- phase space state.
--
-- Desired solution times provided as a list instead of a sized 'V.Vector'.
-- The output list should be the same length as the input list.
evolveHam'
    :: forall m n. (KnownNat m, KnownNat n)
    => System m n  -- ^ system to simulate
    -> Phase n     -- ^ initial state, in phase space
    -> [Double]    -- ^ desired solution times
    -> [Phase n]
evolveHam' _ _ [] = []
evolveHam' s p0 ts = V.withSizedList (toList ts') $ \(v :: V.Vector s Double) ->
                       case Proxy @2 %<=? Proxy @s of
                         LE Refl -> (if l1 then tail else id)
                                  . toList
                                  $ evolveHam s p0 v
                         NLE{}   -> error "evolveHam': Internal error"
  where
    (l1, ts') = case ts of
      [x] -> (True , [0,x])
      _   -> (False, ts   )

-- | Evolve a system using a hamiltonian stepper, with the given initial
-- phase space state.
evolveHam
    :: forall m n s. (KnownNat m, KnownNat n, KnownNat s, 2 <= s)
    => System m n           -- ^ system to simulate
    -> Phase n              -- ^ initial state, in phase space
    -> V.Vector s Double    -- ^ desired solution times
    -> V.Vector s (Phase n)
evolveHam s p0 ts = fmap toPs . fromJust . V.fromList . LA.toRows
                  $ odeSolveV RKf45 hi eps eps (const f) (fromPs p0) ts'
  where
    hi  = (V.unsafeIndex ts 1 - V.unsafeIndex ts 0) / 100
    eps = 1.49012e-08
    f :: LA.Vector Double -> LA.Vector Double
    f   = uncurry (\p m -> LA.vjoin [p,m])
        . join bimap extract . hamEqs s . toPs
    ts' = VG.fromSized . VG.convert $ ts
    n = fromInteger $ natVal (Proxy @n)
    fromPs :: Phase n -> LA.Vector Double
    fromPs p = LA.vjoin . map extract $ [phsPositions p, phsMomenta p]
    toPs :: LA.Vector Double -> Phase n
    toPs v = Phs pP pM
      where
        Just [pP, pM] = traverse create . LA.takesV [n, n] $ v

-- | A convenience wrapper for 'evolveHam'' that works on configuration
-- space states instead of phase space states.
--
-- Note that the simulation itself still runs in phase space; this function
-- just abstracts over converting to and from phase space for the inputs
-- and outputs.
evolveHamC'
    :: forall m n. (KnownNat m, KnownNat n)
    => System m n       -- ^ system to simulate
    -> Config n         -- ^ initial state, in configuration space
    -> [Double]         -- ^ desired solution times
    -> [Config n]
evolveHamC' s c0 = fmap (fromPhase s) . evolveHam' s (toPhase s c0)

-- | A convenience wrapper for 'evolveHam' that works on configuration
-- space states instead of phase space states.
--
-- Note that the simulation itself still runs in phase space; this function
-- just abstracts over converting to and from phase space for the inputs
-- and outputs.
evolveHamC
    :: forall m n s. (KnownNat m, KnownNat n, KnownNat s, 2 <= s)
    => System m n           -- ^ system to simulate
    -> Config n             -- ^ initial state, in configuration space
    -> V.Vector s Double    -- ^ desired solution times
    -> V.Vector s (Config n)
evolveHamC s c0 = fmap (fromPhase s) . evolveHam s (toPhase s c0)

-- | Step a system through configuration space over over a single timestep.
--
-- Note that the simulation itself still runs in phase space; this function
-- just abstracts over converting to and from phase space for the input
-- and output.
stepHamC
    :: forall m n. (KnownNat m, KnownNat n)
    => Double           -- ^ timestep to step through
    -> System m n       -- ^ system to simulate
    -> Config n         -- ^ initial state, in phase space
    -> Config n
stepHamC r s = fromPhase s . stepHam r s . toPhase s

max_stars_repo_name: mstksg/hamilton

max_issues_repo_name: mstksg/hamilton

max_forks_repo_name: mstksg/hamilton