text:
-- Copyright (c) 2015 James M. Lawrence
--
-- Permission is hereby granted, free of charge, to any person obtaining
-- a copy of this software and associated documentation files (the
-- "Software"), to deal in the Software without restriction, including
-- without limitation the rights to use, copy, modify, merge, publish,
-- distribute, sublicense, and/or sell copies of the Software, and to
-- permit persons to whom the Software is furnished to do so, subject to
-- the following conditions:
--
-- The above copyright notice and this permission notice shall be included
-- in all copies or substantial portions of the Software.
--
-- THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
-- EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
-- MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
-- IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
-- CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
-- TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
-- SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
--
-----------------------------------------------------------------------------
-- |
-- Module      :  Math.CayleyDickson
-- Copyright   :  (c) James M. Lawrence
-- License     :  MIT
--
-- Maintainer  :  James M. Lawrence <llmjjmll@gmail.com>
-- Stability   :  provisional
-- Portability :  portable
--
-- Cayley-Dickson constructions (complex numbers, quaternions,
-- octonions, sedenions, etc.) over general scalars without limit to
-- the number of dimensions.
--
-- An element of this structure is composed of an m-dimensional
-- /scalar part/ and an m*(2^n - 1)-dimensional /pure part/ (unrelated
-- to Haskell's uses of "pure"). An element whose scalar part is zero
-- is called a /pure/. Construction with real scalars yields the
-- Cayley-Dickson algebras, in which case the scalar part is also
-- called the /real part/. Other structures may be obtained by
-- considering general scalars, for instance the quaternions over
-- complex scalars.
-----------------------------------------------------------------------------

module Math.CayleyDickson (
    -- * Types
    Nion, Complex, Quaternion, Octonion, Sedenion,

    -- * Construction
    nion, fromScalar, complex, quaternion, octonion, sedenion,

    -- * Operations
    dot, cross, sqnorm, norm, polar,
    -- ** Operations with scalars
    --
    -- | The mnemonic is that the period (".") is on the side of the
    -- scalar.
    (**.),
    (.+), (+.), (.-), (-.), (.*), (*.), (/.),

    -- * Accessors
    coord, coords, setCoord, scalarPart, purePart,

    -- * Constants
    basisElement,

    -- * Classes
    Conjugable(..),

    -- ** Tags
    Tag(tagVal),
    Tag0,  Tag1,  Tag2,  Tag3,  Tag4,  Tag5,  Tag6,  Tag7,  Tag8,  Tag9,
    Tag10, Tag11, Tag12, Tag13, Tag14, Tag15, Tag16, Tag17, Tag18, Tag19,
    Tag20, Tag21, Tag22, Tag23, Tag24, Tag25, Tag26, Tag27, Tag28, Tag29,
    Tag30,

  ) where

----------------------------------------------------------
-- import

import Data.List (genericSplitAt, genericTake, genericReplicate, genericLength)
import Data.Bits (Bits, testBit)
import Data.Proxy (Proxy(Proxy))
import qualified Data.Int as Z
import qualified Data.Ratio as Q
import qualified Data.Complex as C
import qualified Data.Fixed as F
import qualified Data.Word as W

----------------------------------------------------------
-- infix

infix 6 :@

infix 6 .+
infix 6 +.
infix 6 .-
infix 6 -.

infix 7 .*
infix 7 *.
infix 7 /.

infixr 8 **.

----------------------------------------------------------
-- Nion

-- | General Cayley-Dickson construction producing \"N-ions\". The
-- first parameter is a 'Tag' instance that determines the dimension,
-- which is 2 raised to 'tagVal'. The second parameter is the scalar
-- type.
data Nion n a = Scalar a | Nion n a :@ Nion n a

----------------------------------------------------------
-- basic operations

-- | Equivalent to @'coord' x 0@.
scalarPart :: Nion n a -> a
scalarPart (Scalar x) = x
scalarPart (x :@ _) = scalarPart x

-- | Equivalent to @'setCoord' x 0 0@.
purePart :: Num a => Nion n a -> Nion n a
purePart (Scalar _) = Scalar 0
purePart (x :@ y) = purePart x :@ y

-- | Dot product. For real scalars, @(1 \`dot\`)@ is equivalent to
-- @'scalarPart'@.
dot :: Conjugable a => Nion n a -> Nion n a -> a
-- (y * conj x + x * conj y) / 2
Scalar x `dot` Scalar y = scalarPart' $ y * conj x
x@(Scalar _) `dot` (y1 :@ _) = x `dot` y1
(x1 :@ _) `dot` y@(Scalar _) = x1 `dot` y
(x1 :@ x2) `dot` (y1 :@ y2) = (x1 `dot` y1) + (x2 `dot` y2)

-- | Cross product. The cross product of two pures yields a pure that
-- is orthogonal to both operands. For real scalars, @(1 \`cross\`)@
-- is equivalent to @'purePart'@.
cross :: Conjugable a => Nion n a -> Nion n a -> Nion n a
-- (y * conj x - x * conj y) / 2
x `cross` y = y * conj x -. x `dot` y

-- | Squared norm: the dot product of an element with itself.
sqnorm :: Conjugable a => Nion n a -> a
sqnorm x = x `dot` x

-- | Square root of @sqnorm@.
norm :: (Conjugable a, Floating a) => Nion n a -> a
norm = sqrt . sqnorm

-- | Promote a scalar, returning an element whose scalar part is the
-- argument and whose pure part is zero. The element behaves as if it
-- were padded with zeros, but no actual padding is done.
fromScalar :: a -> Nion n a
fromScalar = Scalar

----------------------------------------------------------
-- operations with scalars

leftScalarOp :: (Nion n a -> Nion n a -> Nion n a) -> a -> Nion n a -> Nion n a
leftScalarOp f x y = f (Scalar x) y

rightScalarOp :: (Nion n a -> Nion n a -> Nion n a) -> Nion n a -> a -> Nion n a
rightScalarOp f x y = f x (Scalar y)

-- | Equivalent to @'fromScalar' x + y@.
(.+) :: Conjugable a => a -> Nion n a -> Nion n a
(.+) = leftScalarOp (+)

-- | Equivalent to @'fromScalar' x - y@.
(.-) :: Conjugable a => a -> Nion n a -> Nion n a
(.-) = leftScalarOp (-)

-- | Equivalent to @'fromScalar' x * y@.
(.*) :: Conjugable a => a -> Nion n a -> Nion n a
(.*) = leftScalarOp (*)

-- | Equivalent to @x + 'fromScalar' y@.
(+.) :: Conjugable a => Nion n a -> a -> Nion n a
(+.) = rightScalarOp (+)

-- | Equivalent to @x - 'fromScalar' y@.
(-.) :: Conjugable a => Nion n a -> a -> Nion n a
(-.) = rightScalarOp (-)

-- | Equivalent to @x * 'fromScalar' y@.
(*.) :: Conjugable a => Nion n a -> a -> Nion n a
(*.) = rightScalarOp (*)

-- | Equivalent to @x / 'fromScalar' y@.
(/.) :: (Conjugable a, Fractional a) => Nion n a -> a -> Nion n a
(/.) = rightScalarOp (/)

-- | Raise to a scalar power.
(**.) :: (Tag n, Conjugable a, RealFloat a) => Nion n a -> a -> Nion n a
Scalar x **. y = Scalar $ x ** y
x **. y = exp (y .* log x)

----------------------------------------------------------
-- polar form and complex function application

realPolar :: RealFloat a => Nion n a -> a -> (a, a, Nion n a)
realPolar sqrtMinus1 r
  | r >= 0 = (r, 0, sqrtMinus1)
  | otherwise = (-r, pi, sqrtMinus1)

polarUsing :: (Conjugable a, RealFloat a) =>
              Nion n a -> Nion n a -> (a, a, Nion n a)
polarUsing sqrtMinus1 (Scalar r) = realPolar sqrtMinus1 r
polarUsing sqrtMinus1 x
  | sqnormp == 0 = realPolar sqrtMinus1 r
  | otherwise = (normx, acos (r / normx), u)
  where
    p = purePart x
    sqnormp = sqnorm p
    u = p /. sqrt sqnormp
    r = scalarPart x
    normx = sqrt $ sqnormp + r * conj r

polar' :: (Tag n, Conjugable a, RealFloat a) =>
          Proxy n -> Nion n a -> (a, a, Nion n a)
polar' n x
  | tagVal n == 0 = error "polar: no polar form in 1 dimension"
  | otherwise = polarUsing basisElement1 x

-- | Return @(s, t, u)@ such that (approximately)
--
--     @x == s .* 'exp' (t .* u)@
--
-- where @s@ and @t@ are scalars, @s >= 0@, and @u@ is a unit pure.
--
-- If @x@ has no pure part then @u@ is arbitrarily chosen to be the
-- first pure basis element.
polar :: (Tag n, Conjugable a, RealFloat a) => Nion n a -> (a, a, Nion n a)
polar = polar' Proxy

applyUsing' :: (Tag n, Conjugable a, RealFloat a) =>
               Proxy n ->
               Nion n a -> (a -> a) -> (C.Complex a -> C.Complex a) ->
               Nion n a -> Nion n a
applyUsing' n sqrtMinus1 fr f z
  | tagVal n == 0 = Scalar $ fr $ scalarPart z
  | otherwise = x .+ u *. y
  where (s, t, u) = polarUsing sqrtMinus1 z
        -- handle special cases for a little more accuracy
        x C.:+ y | t == 0 = f $ c s 0
                 | t == pi = f $ c (-s) 0
                 | otherwise = f $ c s 0 * exp (c t 0 * c 0 1)
                 where c = (C.:+)

applyUsing :: (Tag n, Conjugable a, RealFloat a) =>
              Nion n a -> (a -> a) -> (C.Complex a -> C.Complex a) ->
              Nion n a -> Nion n a
applyUsing = applyUsing' Proxy

----------------------------------------------------------
-- constants

fill' :: Tag n => Proxy n -> a -> Nion n a
fill' n s = f $ tagVal n where
  f 0 = Scalar s
  f k = f k' :@ f k' where k' = k - 1

fill :: Tag n => a -> Nion n a
fill = fill' Proxy

paddedZero :: (Tag n, Num a) => Nion n a
paddedZero = fill 0

validIndex :: (Tag n, Num b, Ord b) => Proxy n -> b -> Bool
validIndex n i = i >= 0 && i < 2 ^ tagVal n

basisElement' :: (Tag n, Conjugable a, Bits i, Integral i) =>
                 Proxy n -> i -> Nion n a
basisElement' _ 0 = Scalar 1
basisElement' n index
  | validIndex n index = setCoord paddedZero index 1
  | otherwise = error "basisElement: out of range"

-- | The nth basis element.
basisElement :: (Tag n, Conjugable a, Bits i, Integral i) => i -> Nion n a
basisElement = basisElement' Proxy

basisElement1 :: (Tag n, Conjugable a) => Nion n a
basisElement1 = basisElement (1 :: Integer)

----------------------------------------------------------
-- util

-- Proper Functor and Foldable instances would have to translate
-- top-level scalars to their equivalent representation with padded
-- zeros. The machinations needed for this indirection are rather
-- cumbersome relative to the benefits of having the instances, whose
-- use would seem uncommon.

smap :: (a -> a) -> Nion n a -> Nion n a
smap f (Scalar s) = Scalar $ f s
smap f (x :@ y) = smap f x :@ smap f y

sfoldr :: (a -> b -> b) -> b -> Nion n a -> b
sfoldr f acc (Scalar x) = f x acc
sfoldr f acc (x :@ y) = sfoldr f (sfoldr f acc y) x

----------------------------------------------------------
-- accessors

coords' :: (Tag n, Num a) => Proxy n -> Nion n a -> [a]
coords' n (Scalar x) = x : genericReplicate k 0 where
                         k = 2 ^ tagVal n - 1 :: Integer
coords' _ x = sfoldr (:) [] x

-- | List of coordinates for this element.
coords :: (Tag n, Num a) => Nion n a -> [a]
coords = coords' Proxy

coord' :: (Tag n, Num a, Integral b, Bits b) => Proxy n -> Nion n a -> b -> a
coord' _ (Scalar x) 0 = x
coord' n elt index
  | validIndex n index = case elt of
                           Scalar _ -> 0
                           _ -> f elt $ fromInteger $ tagVal n - 1
  | otherwise = error "coord: out of range"
  where
    f (Scalar x) _ = x
    f (x :@ y) k = case testBit index k of
                     False -> f x k'
                     True  -> f y k'
                   where k' = k - 1

-- | Get the nth coordinate.
coord :: (Tag n, Num a, Integral b, Bits b) => Nion n a -> b -> a
coord = coord' Proxy

setCoord' :: (Tag n, Conjugable a, Integral b, Bits b) =>
             Proxy n -> Nion n a -> b -> a -> Nion n a
setCoord' _ (Scalar _) 0 value = Scalar value
setCoord' n elt index value
  | validIndex n index = case elt of
                           Scalar x -> setCoord (x .+ paddedZero) index value
                           _ -> f elt $ fromInteger $ tagVal n - 1
  | otherwise = error "setCoord: out of range"
  where
    f (Scalar _) _ = Scalar value
    f (x :@ y) k = case testBit index k of
                     False -> f x k' :@ y
                     True  -> x :@ f y k'
                   where k' = k - 1

-- | Set the nth coordinate, returning a new element.
setCoord :: (Tag n, Conjugable a, Integral b, Bits b) =>
            Nion n a -> b -> a -> Nion n a
setCoord = setCoord' Proxy

----------------------------------------------------------
-- construction

fromList :: Integer -> [a] -> Nion n a
fromList _ (x:[]) = Scalar x
fromList k xs = fromList k' l :@ fromList k' r where
                  k' = k `div` 2
                  (l, r) = genericSplitAt k' xs

nion' :: (Tag n, Num a) => Proxy n -> [a] -> Nion n a
nion' n elems = fromList d $ taken ++ padding where
                  d = 2 ^ tagVal n
                  taken = genericTake d elems
                  padding = genericReplicate (d - genericLength taken) 0

-- | Construct an element from a list of coordinates. If the list is
-- too small then the remaining coordinates are padded with zeros. If
-- the list is too large then the extra values are ignored.
nion :: (Tag n, Num a) => [a] -> Nion n a
nion = nion' Proxy

----------------------------------------------------------
-- instances

instance (Tag n, Show a, Num a) => Show (Nion n a) where
  show x = "nion " ++ show (coords x)

instance (Conjugable a, Eq a) => Eq (Nion n a) where
  Scalar x == Scalar y = x == y
  x@(Scalar _) == y1 :@ y2 = x == y1 && y2 == 0
  x1 :@ x2 == y@(Scalar _) = x1 == y && x2 == 0
  x1 :@ x2 == y1 :@ y2 = x1 == y1 && x2 == y2

instance Conjugable a => Num (Nion n a) where
  Scalar x + Scalar y = Scalar $ x + y
  x@(Scalar _) + (y1 :@ y2) = (x + y1) :@ y2
  (x1 :@ x2) + y@(Scalar _) = (x1 + y) :@ x2
  (x1 :@ y1) + (x2 :@ y2) = (x1 + x2) :@ (y1 + y2)

  Scalar x - Scalar y = Scalar $ x - y
  x@(Scalar _) - (y1 :@ y2) = (x - y1) :@ negate y2
  (x1 :@ x2) - y@(Scalar _) = (x1 - y) :@ x2
  (x1 :@ y1) - (x2 :@ y2) = (x1 - x2) :@ (y1 - y2)

  Scalar x * Scalar y = Scalar $ x * y
  x@(Scalar _) * (y1 :@ y2) = (x * y1) :@ (x * y2)
  (x1 :@ x2) * y@(Scalar _) = (x1 * y) :@ (x2 * y)
  (x1 :@ x2) * (y1 :@ y2) = (x1 * y1 - conj y2 * x2) :@ (y2 * x1 + x2 * conj y1)

  negate = smap negate
  fromInteger = fromScalar . fromInteger
  abs = doNotUse
  signum = doNotUse

instance (Conjugable a, Fractional a) => Fractional (Nion n a) where
  Scalar x / Scalar y = Scalar $ x / y
  x@(_ :@ _) / Scalar y = smap (/ y) x
  x / y = (x * conj y) /. sqnorm y

  recip x = conj x /. sqnorm x
  fromRational = fromScalar . fromRational

-- | The first pure basis element is arbitrarily chosen as sqrt (-1).
instance (Tag n, Conjugable a, RealFloat a) => Floating (Nion n a) where
  pi    = Scalar pi
  exp   = applyUsing basisElement1 exp exp
  log   = applyUsing basisElement1 log log
  sqrt  = applyUsing basisElement1 sqrt sqrt
  sin   = applyUsing basisElement1 sin sin
  cos   = applyUsing basisElement1 cos cos
  tan   = applyUsing basisElement1 tan tan
  asin  = applyUsing basisElement1 asin asin
  acos  = applyUsing basisElement1 acos acos
  atan  = applyUsing basisElement1 atan atan
  sinh  = applyUsing basisElement1 sinh sinh
  cosh  = applyUsing basisElement1 cosh cosh
  tanh  = applyUsing basisElement1 tanh tanh
  asinh = applyUsing basisElement1 asinh asinh
  acosh = applyUsing basisElement1 acosh acosh
  atanh = applyUsing basisElement1 atanh atanh

----------------------------------------------------------
-- convenience types

-- | Complex numbers, the 2^1-dimensional construction.
type Complex a = Nion Tag1 a

-- | Quaternions, the 2^2-dimensional construction.
type Quaternion a = Nion Tag2 a

-- | Octonions, the 2^3-dimensional construction.
type Octonion a = Nion Tag3 a

-- | Sedenions, the 2^4-dimensional construction.
type Sedenion a = Nion Tag4 a

-- | Construct a complex number.
complex :: a -> a -> Complex a
complex x y = (:@) (Scalar x) (Scalar y)

-- | Construct a quaternion.
quaternion :: a -> a -> a -> a -> Quaternion a
quaternion w x y z = (:@) ((:@) (Scalar w) (Scalar x))
                          ((:@) (Scalar y) (Scalar z))

-- | Construct an octonion.
octonion :: a -> a -> a -> a ->
            a -> a -> a -> a -> Octonion a
octonion s t u v
         w x y z = (:@) ((:@) ((:@) (Scalar s) (Scalar t))
                              ((:@) (Scalar u) (Scalar v)))
                        ((:@) ((:@) (Scalar w) (Scalar x))
                              ((:@) (Scalar y) (Scalar z)))

-- | Construct a sedenion.
sedenion :: a -> a -> a -> a ->
            a -> a -> a -> a ->
            a -> a -> a -> a ->
            a -> a -> a -> a -> Sedenion a
sedenion k l m n
         o p q r
         s t u v
         w x y z = (:@) ((:@) ((:@) ((:@) (Scalar k) (Scalar l))
                                    ((:@) (Scalar m) (Scalar n)))
                              ((:@) ((:@) (Scalar o) (Scalar p))
                                    ((:@) (Scalar q) (Scalar r))))
                        ((:@) ((:@) ((:@) (Scalar s) (Scalar t))
                                    ((:@) (Scalar u) (Scalar v)))
                              ((:@) ((:@) (Scalar w) (Scalar x))
                                    ((:@) (Scalar y) (Scalar z))))

----------------------------------------------------------
-- Conjugable

-- | The /conjugate/ of an element is obtained by negating the pure
-- part and conjugating the scalar part. The conjugate of a real
-- number is the identity ('id'), which is the default implementation.
class Num a => Conjugable a where
  -- | Conjugate.
  conj :: a -> a
  conj = id

  -- | Equivalent to @(x + conj x) / 2@.
  scalarPart' :: a -> a
  scalarPart' = id

instance Conjugable a => Conjugable (Nion n a) where
  conj (Scalar x) = Scalar $ conj x
  conj (x :@ y) = conj x :@ negate y

  scalarPart' (Scalar x) = Scalar $ scalarPart' x
  scalarPart' (x :@ _) = scalarPart' x

instance (Conjugable a, RealFloat a) => Conjugable (C.Complex a) where
  conj = C.conjugate
  scalarPart' (x C.:+ _) = scalarPart' x C.:+ 0

instance Conjugable Int
instance Conjugable Integer
instance Conjugable Float
instance Conjugable Double
instance Conjugable Z.Int8
instance Conjugable Z.Int16
instance Conjugable Z.Int32
instance Conjugable Z.Int64
instance Conjugable W.Word8
instance Conjugable W.Word16
instance Conjugable W.Word32
instance Conjugable W.Word64
instance Integral a => Conjugable (Q.Ratio a)
instance F.HasResolution a => Conjugable (F.Fixed a)

-----------------------------------------------------------------------------
-- doNotUse

rant :: String
rant = unlines $
  ["",
   "The Num class is a bit messed up, having tied (+), (-), and (*) to abs",
   "and signum. Number systems that have no appropriate definition for abs",
   "or signum must either invent their own operators for addition,",
   "subtraction, and multiplication, else break the contract with Num by",
   "raising an error such as this one when someone uses abs or signum.",
   "",
   "For some time I resisted hijacking Num, but eventually the replacement",
   "operators became too cumbersome and, coupled with the lack of numeric",
   "promotion, significantly detracted from the usability of the package.",
   "So here we are. Good luck, and stay away from abs and signum, which",
   "officially have cooties."]

doNotUse :: a -> a
doNotUse _ = error rant

----------------------------------------------------------
-- Tag

-- | Tags serve to determine a type's dimension, which is 2 raised to
-- `tagVal`. Tag instances are included for convenience only, as you
-- may create your own tag.
class Tag n where
  tagVal :: Proxy n -> Integer

data     Tag0
data     Tag1
data     Tag2
data     Tag3
data     Tag4
data     Tag5
data     Tag6
data     Tag7
data     Tag8
data     Tag9
data    Tag10
data    Tag11
data    Tag12
data    Tag13
data    Tag14
data    Tag15
data    Tag16
data    Tag17
data    Tag18
data    Tag19
data    Tag20
data    Tag21
data    Tag22
data    Tag23
data    Tag24
data    Tag25
data    Tag26
data    Tag27
data    Tag28
data    Tag29
data    Tag30

instance Tag    Tag0 where tagVal _ =    0
instance Tag    Tag1 where tagVal _ =    1
instance Tag    Tag2 where tagVal _ =    2
instance Tag    Tag3 where tagVal _ =    3
instance Tag    Tag4 where tagVal _ =    4
instance Tag    Tag5 where tagVal _ =    5
instance Tag    Tag6 where tagVal _ =    6
instance Tag    Tag7 where tagVal _ =    7
instance Tag    Tag8 where tagVal _ =    8
instance Tag    Tag9 where tagVal _ =    9
instance Tag   Tag10 where tagVal _ =   10
instance Tag   Tag11 where tagVal _ =   11
instance Tag   Tag12 where tagVal _ =   12
instance Tag   Tag13 where tagVal _ =   13
instance Tag   Tag14 where tagVal _ =   14
instance Tag   Tag15 where tagVal _ =   15
instance Tag   Tag16 where tagVal _ =   16
instance Tag   Tag17 where tagVal _ =   17
instance Tag   Tag18 where tagVal _ =   18
instance Tag   Tag19 where tagVal _ =   19
instance Tag   Tag20 where tagVal _ =   20
instance Tag   Tag21 where tagVal _ =   21
instance Tag   Tag22 where tagVal _ =   22
instance Tag   Tag23 where tagVal _ =   23
instance Tag   Tag24 where tagVal _ =   24
instance Tag   Tag25 where tagVal _ =   25
instance Tag   Tag26 where tagVal _ =   26
instance Tag   Tag27 where tagVal _ =   27
instance Tag   Tag28 where tagVal _ =   28
instance Tag   Tag29 where tagVal _ =   29
instance Tag   Tag30 where tagVal _ =   30

max_stars_repo_name: lmj/cayley-dickson

max_issues_repo_name: lmj/cayley-dickson

max_forks_repo_name: lmj/cayley-dickson


text:
module Fibon.Analyse.Statistics(
    computeSummary
  , Summary(..)
  , Estimate(..)
  , ConfidenceInterval(..)
)
where

import Statistics.Sample
import qualified Data.Vector.Unboxed as V

data Summary =
    Min
  | GeoMean
  | ArithMean
  | Max
  deriving(Read, Show)

data Estimate a = Estimate {
      ePoint  :: !a
    , eStddev :: !a
    , eSize   :: !Int
    , eCI     :: Maybe (ConfidenceInterval a)
  }
  deriving (Read, Show)

data ConfidenceInterval a = ConfidenceInterval {
      eLowerBound       :: !a
    , eUpperBound       :: !a
    , eConfidenceLevel  :: !Double
}
  deriving (Read, Show)

instance Functor Estimate where
  fmap f e = e {
      ePoint  = f (ePoint e)
    , eStddev = f (eStddev e)
    , eCI     = maybe Nothing (Just . fmap f) (eCI e)
  }

instance Functor ConfidenceInterval where
  fmap f c = c {
      eLowerBound = f (eLowerBound c)
    , eUpperBound = f (eUpperBound c)
  }

computeSummary :: Summary -> Sample -> Estimate Double
computeSummary summaryType vec =
  Estimate {
      ePoint  = sumF summaryType vec
    , eStddev = stdDev vec -- TODO: this is wrong stddev for geoMean
    , eSize   = V.length vec
    , eCI     = Nothing
  }
  where
  sumF ArithMean = mean
  sumF GeoMean   = geometricMean
  sumF Max       = V.maximum
  sumF Min       = V.minimum

ext: hs

max_stars_repo_name: dmpots/fibon

max_issues_repo_name: dmpots/fibon

max_forks_repo_name: dmpots/fibon


text:
module TestRotation where

import Cube (Cube)
import Numeric.LinearAlgebra ((<>))
import Rotation
import qualified Cube as Cube
import qualified Data.List as List

testRotations :: MatrixRotation -> Bool
testRotations m =
  let ms = take 4 $ iterate (<> m) m
      areIdentity = fmap (== identity) ms
      testAreIdentity = areIdentity == [False, False, False, True]
      testAllDifferent = (4 ==) . length $ List.nub ms
  in testAreIdentity && testAllDifferent

testTopToFront :: Bool
testTopToFront =
  let [top, frn, rgt, bck, lft, bot] = [Cube.Top ..]
      [w, g, r, b, o, y] = [Cube.White ..]
      expected = [ (top, 0, 0, b), (top, 0, 1, b), (top, 0, 2, b), (top, 1, 0, b), (top, 1, 1, b), (top, 1, 2, b), (top, 2, 0, b), (top, 2, 1, b), (top, 2, 2, b)
                 , (frn, 0, 0, w), (frn, 0, 1, w), (frn, 0, 2, w), (frn, 1, 0, w), (frn, 1, 1, w), (frn, 1, 2, w), (frn, 2, 0, w), (frn, 2, 1, w), (frn, 2, 2, w)
                 , (rgt, 0, 0, r), (rgt, 0, 1, r), (rgt, 0, 2, r), (rgt, 1, 0, r), (rgt, 1, 1, r), (rgt, 1, 2, r), (rgt, 2, 0, r), (rgt, 2, 1, r), (rgt, 2, 2, r)
                 , (bck, 0, 0, y), (bck, 0, 1, y), (bck, 0, 2, y), (bck, 1, 0, y), (bck, 1, 1, y), (bck, 1, 2, y), (bck, 2, 0, y), (bck, 2, 1, y), (bck, 2, 2, y)
                 , (lft, 0, 0, o), (lft, 0, 1, o), (lft, 0, 2, o), (lft, 1, 0, o), (lft, 1, 1, o), (lft, 1, 2, o), (lft, 2, 0, o), (lft, 2, 1, o), (lft, 2, 2, o)
                 , (bot, 0, 0, g), (bot, 0, 1, g), (bot, 0, 2, g), (bot, 1, 0, g), (bot, 1, 1, g), (bot, 1, 2, g), (bot, 2, 0, g), (bot, 2, 1, g), (bot, 2, 2, g) ]
      rotated = rotate topToFront Cube.solvedCube
  in expected == rotated
  && testRotations topToFront

getSideColors :: Cube -> [(Cube.Side, Cube.Color)]
getSideColors = fmap head . filter (not . null) . List.group . fmap (\(s,_,_,c) -> (s,c))

testTopToBack :: Bool
testTopToBack = testRotations topToBack

testTopToRight :: Bool
testTopToRight =
  let [top, frn, rgt, bck, lft, bot] = [Cube.Top ..]
      [w, g, r, b, o, y] = [Cube.White ..]
      expected = [ (top, 0, 0, o), (top, 0, 1, o), (top, 0, 2, o), (top, 1, 0, o), (top, 1, 1, o), (top, 1, 2, o), (top, 2, 0, o), (top, 2, 1, o), (top, 2, 2, o)
                 , (frn, 0, 0, g), (frn, 0, 1, g), (frn, 0, 2, g), (frn, 1, 0, g), (frn, 1, 1, g), (frn, 1, 2, g), (frn, 2, 0, g), (frn, 2, 1, g), (frn, 2, 2, g)
                 , (rgt, 0, 0, w), (rgt, 0, 1, w), (rgt, 0, 2, w), (rgt, 1, 0, w), (rgt, 1, 1, w), (rgt, 1, 2, w), (rgt, 2, 0, w), (rgt, 2, 1, w), (rgt, 2, 2, w)
                 , (bck, 0, 0, b), (bck, 0, 1, b), (bck, 0, 2, b), (bck, 1, 0, b), (bck, 1, 1, b), (bck, 1, 2, b), (bck, 2, 0, b), (bck, 2, 1, b), (bck, 2, 2, b)
                 , (lft, 0, 0, y), (lft, 0, 1, y), (lft, 0, 2, y), (lft, 1, 0, y), (lft, 1, 1, y), (lft, 1, 2, y), (lft, 2, 0, y), (lft, 2, 1, y), (lft, 2, 2, y)
                 , (bot, 0, 0, r), (bot, 0, 1, r), (bot, 0, 2, r), (bot, 1, 0, r), (bot, 1, 1, r), (bot, 1, 2, r), (bot, 2, 0, r), (bot, 2, 1, r), (bot, 2, 2, r) ]
      rotated = rotate topToRight Cube.solvedCube
  in expected == rotated
  && testRotations topToRight

testTopToLeft :: Bool
testTopToLeft = testRotations topToLeft

testTopTwistRight :: Bool
testTopTwistRight =
  let [top, frn, rgt, bck, lft, bot] = [Cube.Top ..]
      [w, g, r, b, o, y] = [Cube.White ..]
      expected = [ (top, 0, 0, w), (top, 0, 1, w), (top, 0, 2, w), (top, 1, 0, w), (top, 1, 1, w), (top, 1, 2, w), (top, 2, 0, w), (top, 2, 1, w), (top, 2, 2, w)
                 , (frn, 0, 0, r), (frn, 0, 1, r), (frn, 0, 2, r), (frn, 1, 0, r), (frn, 1, 1, r), (frn, 1, 2, r), (frn, 2, 0, r), (frn, 2, 1, r), (frn, 2, 2, r)
                 , (rgt, 0, 0, b), (rgt, 0, 1, b), (rgt, 0, 2, b), (rgt, 1, 0, b), (rgt, 1, 1, b), (rgt, 1, 2, b), (rgt, 2, 0, b), (rgt, 2, 1, b), (rgt, 2, 2, b)
                 , (bck, 0, 0, o), (bck, 0, 1, o), (bck, 0, 2, o), (bck, 1, 0, o), (bck, 1, 1, o), (bck, 1, 2, o), (bck, 2, 0, o), (bck, 2, 1, o), (bck, 2, 2, o)
                 , (lft, 0, 0, g), (lft, 0, 1, g), (lft, 0, 2, g), (lft, 1, 0, g), (lft, 1, 1, g), (lft, 1, 2, g), (lft, 2, 0, g), (lft, 2, 1, g), (lft, 2, 2, g)
                 , (bot, 0, 0, y), (bot, 0, 1, y), (bot, 0, 2, y), (bot, 1, 0, y), (bot, 1, 1, y), (bot, 1, 2, y), (bot, 2, 0, y), (bot, 2, 1, y), (bot, 2, 2, y) ]
      rotated = rotate topTwistRight Cube.solvedCube
  in expected == rotated
  && testRotations topTwistRight

testTopTwistLeft :: Bool
testTopTwistLeft = testRotations topTwistLeft

testRightL :: Bool
testRightL =
  let [top, frn, rgt, bck, lft, bot] = [Cube.Top ..]
      [w, g, r, b, o, y] = [Cube.White ..]
      expected = [ (top, 0, 0, w), (top, 0, 1, w), (top, 0, 2, b), (top, 1, 0, w), (top, 1, 1, w), (top, 1, 2, b), (top, 2, 0, w), (top, 2, 1, w), (top, 2, 2, b)
                 , (frn, 0, 0, g), (frn, 0, 1, g), (frn, 0, 2, w), (frn, 1, 0, g), (frn, 1, 1, g), (frn, 1, 2, w), (frn, 2, 0, g), (frn, 2, 1, g), (frn, 2, 2, w)
                 , (rgt, 0, 0, r), (rgt, 0, 1, r), (rgt, 0, 2, r), (rgt, 1, 0, r), (rgt, 1, 1, r), (rgt, 1, 2, r), (rgt, 2, 0, r), (rgt, 2, 1, r), (rgt, 2, 2, r)
                 , (bck, 0, 0, y), (bck, 0, 1, b), (bck, 0, 2, b), (bck, 1, 0, y), (bck, 1, 1, b), (bck, 1, 2, b), (bck, 2, 0, y), (bck, 2, 1, b), (bck, 2, 2, b)
                 , (lft, 0, 0, o), (lft, 0, 1, o), (lft, 0, 2, o), (lft, 1, 0, o), (lft, 1, 1, o), (lft, 1, 2, o), (lft, 2, 0, o), (lft, 2, 1, o), (lft, 2, 2, o)
                 , (bot, 0, 0, y), (bot, 0, 1, y), (bot, 0, 2, g), (bot, 1, 0, y), (bot, 1, 1, y), (bot, 1, 2, g), (bot, 2, 0, y), (bot, 2, 1, y), (bot, 2, 2, g) ]
      rotated = rotate rightL Cube.solvedCube
  in expected == rotated
  && testRotations rightL

testRightR :: Bool
testRightR = testRotations rightR

testLeftL :: Bool
testLeftL =
  let [top, frn, rgt, bck, lft, bot] = [Cube.Top ..]
      [w, g, r, b, o, y] = [Cube.White ..]
      expected = [ (top, 0, 0, g), (top, 0, 1, w), (top, 0, 2, w), (top, 1, 0, g), (top, 1, 1, w), (top, 1, 2, w), (top, 2, 0, g), (top, 2, 1, w), (top, 2, 2, w)
                 , (frn, 0, 0, y), (frn, 0, 1, g), (frn, 0, 2, g), (frn, 1, 0, y), (frn, 1, 1, g), (frn, 1, 2, g), (frn, 2, 0, y), (frn, 2, 1, g), (frn, 2, 2, g)
                 , (rgt, 0, 0, r), (rgt, 0, 1, r), (rgt, 0, 2, r), (rgt, 1, 0, r), (rgt, 1, 1, r), (rgt, 1, 2, r), (rgt, 2, 0, r), (rgt, 2, 1, r), (rgt, 2, 2, r)
                 , (bck, 0, 0, b), (bck, 0, 1, b), (bck, 0, 2, w), (bck, 1, 0, b), (bck, 1, 1, b), (bck, 1, 2, w), (bck, 2, 0, b), (bck, 2, 1, b), (bck, 2, 2, w)
                 , (lft, 0, 0, o), (lft, 0, 1, o), (lft, 0, 2, o), (lft, 1, 0, o), (lft, 1, 1, o), (lft, 1, 2, o), (lft, 2, 0, o), (lft, 2, 1, o), (lft, 2, 2, o)
                 , (bot, 0, 0, b), (bot, 0, 1, y), (bot, 0, 2, y), (bot, 1, 0, b), (bot, 1, 1, y), (bot, 1, 2, y), (bot, 2, 0, b), (bot, 2, 1, y), (bot, 2, 2, y) ]
      rotated = rotate leftL Cube.solvedCube
  in expected == rotated
  && testRotations leftL

testLeftR :: Bool
testLeftR = testRotations leftR

testTopL :: Bool
testTopL =
  let [top, frn, rgt, bck, lft, bot] = [Cube.Top ..]
      [w, g, r, b, o, y] = [Cube.White ..]
      expected = [ (top, 0, 0, w), (top, 0, 1, w), (top, 0, 2, w), (top, 1, 0, w), (top, 1, 1, w), (top, 1, 2, w), (top, 2, 0, w), (top, 2, 1, w), (top, 2, 2, w)
                 , (frn, 0, 0, o), (frn, 0, 1, o), (frn, 0, 2, o), (frn, 1, 0, g), (frn, 1, 1, g), (frn, 1, 2, g), (frn, 2, 0, g), (frn, 2, 1, g), (frn, 2, 2, g)
                 , (rgt, 0, 0, g), (rgt, 0, 1, g), (rgt, 0, 2, g), (rgt, 1, 0, r), (rgt, 1, 1, r), (rgt, 1, 2, r), (rgt, 2, 0, r), (rgt, 2, 1, r), (rgt, 2, 2, r)
                 , (bck, 0, 0, r), (bck, 0, 1, r), (bck, 0, 2, r), (bck, 1, 0, b), (bck, 1, 1, b), (bck, 1, 2, b), (bck, 2, 0, b), (bck, 2, 1, b), (bck, 2, 2, b)
                 , (lft, 0, 0, b), (lft, 0, 1, b), (lft, 0, 2, b), (lft, 1, 0, o), (lft, 1, 1, o), (lft, 1, 2, o), (lft, 2, 0, o), (lft, 2, 1, o), (lft, 2, 2, o)
                 , (bot, 0, 0, y), (bot, 0, 1, y), (bot, 0, 2, y), (bot, 1, 0, y), (bot, 1, 1, y), (bot, 1, 2, y), (bot, 2, 0, y), (bot, 2, 1, y), (bot, 2, 2, y) ]
      rotated = rotate topL Cube.solvedCube
  in expected == rotated
  && testRotations topL

testTopR :: Bool
testTopR = testRotations topR

testBottomL :: Bool
testBottomL =
  let [top, frn, rgt, bck, lft, bot] = [Cube.Top ..]
      [w, g, r, b, o, y] = [Cube.White ..]
      expected = [ (top, 0, 0, w), (top, 0, 1, w), (top, 0, 2, w), (top, 1, 0, w), (top, 1, 1, w), (top, 1, 2, w), (top, 2, 0, w), (top, 2, 1, w), (top, 2, 2, w)
                 , (frn, 0, 0, g), (frn, 0, 1, g), (frn, 0, 2, g), (frn, 1, 0, g), (frn, 1, 1, g), (frn, 1, 2, g), (frn, 2, 0, r), (frn, 2, 1, r), (frn, 2, 2, r)
                 , (rgt, 0, 0, r), (rgt, 0, 1, r), (rgt, 0, 2, r), (rgt, 1, 0, r), (rgt, 1, 1, r), (rgt, 1, 2, r), (rgt, 2, 0, b), (rgt, 2, 1, b), (rgt, 2, 2, b)
                 , (bck, 0, 0, b), (bck, 0, 1, b), (bck, 0, 2, b), (bck, 1, 0, b), (bck, 1, 1, b), (bck, 1, 2, b), (bck, 2, 0, o), (bck, 2, 1, o), (bck, 2, 2, o)
                 , (lft, 0, 0, o), (lft, 0, 1, o), (lft, 0, 2, o), (lft, 1, 0, o), (lft, 1, 1, o), (lft, 1, 2, o), (lft, 2, 0, g), (lft, 2, 1, g), (lft, 2, 2, g)
                 , (bot, 0, 0, y), (bot, 0, 1, y), (bot, 0, 2, y), (bot, 1, 0, y), (bot, 1, 1, y), (bot, 1, 2, y), (bot, 2, 0, y), (bot, 2, 1, y), (bot, 2, 2, y) ]
      rotated = rotate bottomL Cube.solvedCube
  in expected == rotated
  && testRotations bottomL

testBottomR :: Bool
testBottomR = testRotations bottomR

testFrontL :: Bool
testFrontL =
  let [top, frn, rgt, bck, lft, bot] = [Cube.Top ..]
      [w, g, r, b, o, y] = [Cube.White ..]
      expected = [ (top, 0, 0, w), (top, 0, 1, w), (top, 0, 2, w), (top, 1, 0, w), (top, 1, 1, w), (top, 1, 2, w), (top, 2, 0, r), (top, 2, 1, r), (top, 2, 2, r)
                 , (frn, 0, 0, g), (frn, 0, 1, g), (frn, 0, 2, g), (frn, 1, 0, g), (frn, 1, 1, g), (frn, 1, 2, g), (frn, 2, 0, g), (frn, 2, 1, g), (frn, 2, 2, g)
                 , (rgt, 0, 0, y), (rgt, 0, 1, r), (rgt, 0, 2, r), (rgt, 1, 0, y), (rgt, 1, 1, r), (rgt, 1, 2, r), (rgt, 2, 0, y), (rgt, 2, 1, r), (rgt, 2, 2, r)
                 , (bck, 0, 0, b), (bck, 0, 1, b), (bck, 0, 2, b), (bck, 1, 0, b), (bck, 1, 1, b), (bck, 1, 2, b), (bck, 2, 0, b), (bck, 2, 1, b), (bck, 2, 2, b)
                 , (lft, 0, 0, o), (lft, 0, 1, o), (lft, 0, 2, w), (lft, 1, 0, o), (lft, 1, 1, o), (lft, 1, 2, w), (lft, 2, 0, o), (lft, 2, 1, o), (lft, 2, 2, w)
                 , (bot, 0, 0, o), (bot, 0, 1, o), (bot, 0, 2, o), (bot, 1, 0, y), (bot, 1, 1, y), (bot, 1, 2, y), (bot, 2, 0, y), (bot, 2, 1, y), (bot, 2, 2, y) ]
      rotated = rotate frontL Cube.solvedCube
  in expected == rotated
  && testRotations frontL

testFrontR :: Bool
testFrontR = testRotations frontR

testBackL :: Bool
testBackL =
  let [top, frn, rgt, bck, lft, bot] = [Cube.Top ..]
      [w, g, r, b, o, y] = [Cube.White ..]
      expected = [ (top, 0, 0, o), (top, 0, 1, o), (top, 0, 2, o), (top, 1, 0, w), (top, 1, 1, w), (top, 1, 2, w), (top, 2, 0, w), (top, 2, 1, w), (top, 2, 2, w)
                 , (frn, 0, 0, g), (frn, 0, 1, g), (frn, 0, 2, g), (frn, 1, 0, g), (frn, 1, 1, g), (frn, 1, 2, g), (frn, 2, 0, g), (frn, 2, 1, g), (frn, 2, 2, g)
                 , (rgt, 0, 0, r), (rgt, 0, 1, r), (rgt, 0, 2, w), (rgt, 1, 0, r), (rgt, 1, 1, r), (rgt, 1, 2, w), (rgt, 2, 0, r), (rgt, 2, 1, r), (rgt, 2, 2, w)
                 , (bck, 0, 0, b), (bck, 0, 1, b), (bck, 0, 2, b), (bck, 1, 0, b), (bck, 1, 1, b), (bck, 1, 2, b), (bck, 2, 0, b), (bck, 2, 1, b), (bck, 2, 2, b)
                 , (lft, 0, 0, y), (lft, 0, 1, o), (lft, 0, 2, o), (lft, 1, 0, y), (lft, 1, 1, o), (lft, 1, 2, o), (lft, 2, 0, y), (lft, 2, 1, o), (lft, 2, 2, o)
                 , (bot, 0, 0, y), (bot, 0, 1, y), (bot, 0, 2, y), (bot, 1, 0, y), (bot, 1, 1, y), (bot, 1, 2, y), (bot, 2, 0, r), (bot, 2, 1, r), (bot, 2, 2, r) ]
      rotated = rotate backL Cube.solvedCube
  in expected == rotated
  && testRotations backL

testBackR :: Bool
testBackR = testRotations backR

testCombined :: Bool
testCombined =
  let [top, frn, rgt, bck, lft, bot] = [Cube.Top ..]
      [w, g, r, b, o, y] = [Cube.White ..]
      expected = [ (top, 0, 0, o), (top, 0, 1, o), (top, 0, 2, b), (top, 1, 0, o), (top, 1, 1, w), (top, 1, 2, b), (top, 2, 0, o), (top, 2, 1, r), (top, 2, 2, r)
                 , (frn, 0, 0, b), (frn, 0, 1, g), (frn, 0, 2, w), (frn, 1, 0, y), (frn, 1, 1, g), (frn, 1, 2, w), (frn, 2, 0, y), (frn, 2, 1, y), (frn, 2, 2, b)
                 , (rgt, 0, 0, g), (rgt, 0, 1, r), (rgt, 0, 2, w), (rgt, 1, 0, g), (rgt, 1, 1, r), (rgt, 1, 2, w), (rgt, 2, 0, y), (rgt, 2, 1, y), (rgt, 2, 2, o)
                 , (bck, 0, 0, r), (bck, 0, 1, b), (bck, 0, 2, g), (bck, 1, 0, r), (bck, 1, 1, b), (bck, 1, 2, w), (bck, 2, 0, y), (bck, 2, 1, y), (bck, 2, 2, g)
                 , (lft, 0, 0, w), (lft, 0, 1, w), (lft, 0, 2, w), (lft, 1, 0, b), (lft, 1, 1, o), (lft, 1, 2, o), (lft, 2, 0, y), (lft, 2, 1, g), (lft, 2, 2, r)
                 , (bot, 0, 0, g), (bot, 0, 1, g), (bot, 0, 2, r), (bot, 1, 0, o), (bot, 1, 1, y), (bot, 1, 2, r), (bot, 2, 0, o), (bot, 2, 1, b), (bot, 2, 2, b) ]
      rotation = bottomL <> leftL <> backL <> rightL <> frontL <> topL
      rotated = rotate rotation Cube.solvedCube
  in expected == rotated
  && testRotations backL

testNameMatrixRotation :: Bool
testNameMatrixRotation = and $ do
  (r, n) <- rotations
  let n' = nameMatrixRotation r
  return $ n == n'

testInverse :: Bool
testInverse = and $ do
  (r, _) <- rotations
  let (Just r') = inverse r
  let (Just r'') = inverse r'
  return $ r == r''

tests = [ ("Rotation.topToFront", testTopToFront)
        , ("Rotation.topToBack", testTopToBack)
        , ("Rotation.topToRight", testTopToRight)
        , ("Rotation.topToLeft", testTopToLeft)
        , ("Rotation.topTwistRight", testTopTwistRight)
        , ("Rotation.topTwistLeft", testTopTwistLeft)
        , ("Rotation.rightL", testRightL)
        , ("Rotation.rightR", testRightR)
        , ("Rotation.leftL", testLeftL)
        , ("Rotation.leftR", testLeftR)
        , ("Rotation.topL", testTopL)
        , ("Rotation.topR", testTopR)
        , ("Rotation.bottomL", testBottomL)
        , ("Rotation.bottomR", testBottomR)
        , ("Rotation.frontL", testFrontL)
        , ("Rotation.frontR", testFrontR)
        , ("Rotation.backL", testBackL)
        , ("Rotation.backR", testBackR)
        , ("combined", testCombined)
        , ("Rotation.nameMatrixRotation", testNameMatrixRotation)
        , ("Rotation.inverse", testInverse)
        ]

max_stars_repo_name: runjak/hRubiks

max_issues_repo_name: runjak/hRubiks

max_forks_repo_name: runjak/hRubiks


text:
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE Strict #-}
{-# LANGUAGE FlexibleContexts #-}

module CPVO.IO.DOS (
    showTotPDOS
  ) where

import CPVO.Numeric
import CPVO.IO
import CPVO.IO.Reader.Ecalj.Common
import CPVO.IO.Reader.Ecalj.DOS

import qualified Data.Text as T
import qualified Data.Text.IO as T
import Data.List
import Numeric.LinearAlgebra

showTotPDOS :: [String] -> IO ()
showTotPDOS allArgs = do
--  (invStat, ymax, xmin, xmax, ctrlAtoms, uniqAtoms, ctrlAtomicAOs,jdTable, cleanedJdHead, foldernya, tailer, colAlign, texFile) <- readHeaderData allArgs
  Right (invStat, _ , _ , _ , _ , _ , ctrlAtomicAOs,jdTable, cleanedJdHead, foldernya, tailer, _ , texFile) <- readHeaderData allArgs

  -------------------------------generating DOS data------------------------
  totalDOS <- readTotalDOSText tailer foldernya
  -------------------------------integrating DOS data------------------------
  let intgTot = map (\i -> integrateToZero $ totalDOS ¿ [0,i]) $ flipBy invStat [1,2]
  pdosAtomic <- readPDOS invStat tailer foldernya ctrlAtomicAOs
  let integratedAtPDOS = integrateAtomicPDOS pdosAtomic
  putStrLn $ show integratedAtPDOS
  let rIntgAll' =
        rendertable
        $ (:) cleanedJdHead
        $ nub
        $ (:) (concat [ ["Total"]
                      , map (showDouble (2::Integer)) $ (\[a,b] -> [a,b,a-b,a+b]) $ intgTot
                      ])
        $ map (\(iu,idn,(_,(j,_))) -> j:map (showDouble (2::Integer)) [iu,idn,iu-idn,iu+idn] ) integratedAtPDOS
  let rIntgAll = unlines  [
                          rIntgAll'
                          , jdTable
                          ]
  putStrLn rIntgAll
  {-
    -------------------------------generating LaTex Representation------------------------
  resIntAll' <- markdownToTex rIntgAll
  let resIntAll = T.replace "\\}" "}"
                $ T.replace "\\{" "{" $ T.pack
                $ unlines [
                          "\\begin{longtable}[]{" ++ colAlign ++ "}"
                          , unlines $ tail $ lines $ T.unpack resIntAll'
                          ]
  -}

  T.writeFile texFile $ T.pack rIntgAll
  putStrLn "====done===="

max_stars_repo_name: hasanalrasyid/cpvoh

max_issues_repo_name: hasanalrasyid/cpvoh

max_forks_repo_name: hasanalrasyid/cpvoh


text:
{-# LANGUAGE TemplateHaskell, GeneralizedNewtypeDeriving, ScopedTypeVariables #-}
module CNC.Figure(Pos, Figure, Path,
                  unPath, mapPath, mapPathM_, deltaPath,
                  solveFigure, declarePoints, figure,
                  xSize, ySize, len, xAngle, yAngle) where

import CNC.Geometry(RealT, eps)

import qualified Data.Map as M
import Data.Map((!))
import Data.Maybe(fromJust, isNothing)
import Data.Label
import Data.List(deleteBy,sort)
import qualified Data.Label.Monadic as LM
import qualified Control.Monad.State as S
import Data.Complex
import Text.Printf

newtype PointInd = PointInd Int deriving (Eq, Ord, Enum, Show)

type Pos = Complex RealT

type FigMap = M.Map PointInd (Maybe Pos)
data Figure = Figure { _fPoints ::  FigMap, _fNextPoint :: PointInd, _fConstraints :: [Constraint] } deriving Show
newtype Path = Path {unPath :: [Pos]} deriving Show

mkPos x y = x :+ y

-- applies a given transform to all points of path
mapPath f (Path p) = Path (map f p)
mapPathM_ f (Path p) = mapM_ f p
deltaPath (Path p) = last p - head p

zero_ind = PointInd 0
zero_point = mkPos 0 0

empty_figure = Figure M.empty zero_ind []
figNumPoints fig = case _fNextPoint fig of
  PointInd k -> k

instance (Ord a, RealFloat a) => Ord (Complex a) where
  compare x y = let c1 = compare (realPart x) (realPart y)
                in if c1 == EQ then compare (imagPart x) (imagPart y)
                   else c1

data UniConstraint pt = Disp {cPoint1 :: pt, cPoint2 :: pt, cPos :: Pos} -- A vector from first point to second
                      | Length {cPoint1 :: pt, cPoint2 :: pt, cLen :: RealT} -- distance between points
                      | Angle {cPoint1 :: pt, cPoint2 :: pt, cAngle :: RealT} deriving (Eq, Ord, Show)

instance Functor UniConstraint where
  fmap f c = case c of
    Disp p1 p2 d -> Disp (f p1) (f p2) d
    Length p1 p2 l -> Length (f p1) (f p2) l
    Angle p1 p2 a -> Angle (f p1) (f p2) a

type Constraint = UniConstraint PointInd

mkLabels [''Figure]

figure fig_m = let fig = S.execState fig_m empty_figure
               in solveFigure fig

-- Creates k new points
declarePoints k = do
  next <- LM.gets fNextPoint
  let pts = take k $ iterate succ next
      pairs = map (\p -> (p, Nothing)) pts
  LM.modify fPoints (\ps -> ps `M.union` M.fromList pairs)
  LM.puts fNextPoint (succ $ last pts)
  return pts

xSize p1 p2 d = LM.modify fConstraints (Disp p1 p2 (mkPos d 0) :) -- distance between horizontally aligned points
ySize p1 p2 d = LM.modify fConstraints (Disp p1 p2 (mkPos 0 d) :) -- distance between vertically aligned points
len p1 p2 d = LM.modify fConstraints (Length p1 p2 d :) -- distance between points
xAngle p1 p2 a = LM.modify fConstraints (Angle p1 p2 a :) -- angle between line and axe X
yAngle p1 p2 a = LM.modify fConstraints (Angle p1 p2 (pi/2 - a) :) -- angle between line and axe Y

-- solves system consisting of N declared points and user specified constraints
solveFigure :: Figure -> Path
solveFigure fig | figNumPoints fig == 0 = error "can't solve figure"
                | otherwise = let points = M.insert zero_ind (Just zero_point) $ get fPoints fig
                                  constraints = sort $ get fConstraints fig
                                  solution = solveConstraints points constraints [] False
                              in Path $ map fromJust $ M.elems $ checkSolution solution constraints

solveConstraints :: FigMap -> [Constraint] -> [Constraint] -> Bool -> FigMap
solveConstraints points [] [] _ = points -- no more constraints, finish here
solveConstraints points [] prev_cs True = solveConstraints points (reverse prev_cs) [] False -- iterating over constraints again
solveConstraints points [] prev_cs False = points -- no constraints was applied, finishing
solveConstraints points (c:cs) prev_cs shrinked = case c of -- process a single constraint
  Disp pi1 pi2 disp -> case (points ! pi1, points ! pi2) of
    (Nothing, Nothing) -> skipConstr
    (Just p1, Nothing) -> solveDisp p1 pi2 disp
    (Nothing, Just p2) -> solveDisp p2 pi1 (-disp)
    (Just p1, Just p2) -> checkConstr
  Length pi1 pi2 len -> case (points ! pi1, points ! pi2) of
    (Nothing, Nothing) -> skipConstr
    (Just p1, Nothing) -> solveLength p1 pi1 pi2 len
    (Nothing, Just p2) -> solveLength p2 pi2 pi1 len
    (Just p1, Just p2) -> checkConstr
  Angle pi1 pi2 a -> case (M.lookup pi1 points, M.lookup pi2 points) of -- only check these, they're used as part of Length constraints processing
    (Just p1, Just p2) -> checkConstr
    _ -> skipConstr
  where
    pointConstr = fmap (fromJust . (points !)) c
    checkConstr = checkEps (calcDiscrepancy pointConstr) (solveConstraints points cs prev_cs shrinked) ("disagreement while processing " ++ show c ++ " evaluated to " ++ show pointConstr)
    skipConstr = solveConstraints points cs (c:prev_cs) shrinked
    solveDisp p ind d = solveConstraints (M.insert ind (Just $ p + d) points) cs prev_cs True
    solveLength (p :: Pos) ind0 ind len = case findAngleConstr ind0 ind cs of
      Nothing -> skipConstr
      Just angle -> solveConstraints (M.insert ind (Just $ p + mkPolar len angle) points)
                                     (deleteAngleConstr ind0 ind cs) prev_cs True

checkEps :: RealT -> a -> String -> a
checkEps val cont msg = if abs val < eps then cont else error $ msg

calcDiscrepancy :: (UniConstraint Pos) -> RealT
calcDiscrepancy c =
  case c of
    Disp p1 p2 d -> magnitude $ p2 - p1 - d
    Length p1 p2 l -> magnitude (p2 - p1) - l
    Angle p1 p2 a -> phase (p2 - p1) - a

checkSolution points constrs | any isNothing $ M.elems points = error $ "couldn't find a complete solution, got " ++ show (M.elems points)
checkSolution points cs = checkConstrs points cs

checkConstrs :: FigMap -> [UniConstraint PointInd] -> FigMap
checkConstrs points [] = points
checkConstrs points (c:cs) = check $ calcDiscrepancy $ fmap (fromJust . (points !)) c
 where p1 = cPoint1 c
       p2 = cPoint2 c
       check x = checkEps x (checkConstrs points cs) $
                   printf "error checking constraint %s on points %s %s placed at %s %s"
                   (show c) (show p1) (show p2) (show $ points ! p1) (show $ points ! p2)

-- finds first angle constraint on points with specified indices
findAngleConstr :: PointInd -> PointInd -> [Constraint] -> Maybe RealT
findAngleConstr _ _ [] = Nothing
findAngleConstr ind1 ind2 (c:cs) = case c of
  Angle i1 i2 a | i1 == ind1 && i2 == ind2 -> Just a
                | i1 == ind2 && i2 == ind1 -> Just (-a)
  _ -> findAngleConstr ind1 ind2 cs

-- deletes first angle constraint
deleteAngleConstr ind1 ind2 = deleteBy findAngle
                              (Angle ind1 ind2 undefined)
  where findAngle (Angle d1 d2 _) c2 =
          case c2 of
            (Angle i1 i2 _) -> (i1, i2) == (d2,d1) || (i2, i1) == (d1,d2)
            _ -> False

max_stars_repo_name: akamaus/gcodec

max_issues_repo_name: akamaus/gcodec

max_forks_repo_name: akamaus/gcodec


text:
module CCEA.XOR
  (
    xorFitness
  , xorBreeding
  , xorEA
  , xorVars
  ) where

import NN.NeuralNetwork
import RandomUtil.Random
import AI.HNN.FF.Network

import qualified Data.Map.Strict as Map
import Numeric.LinearAlgebra.HMatrix
import Data.List (sortOn, partition, mapAccumL)
import System.Random
import System.Random.Shuffle
import Data.Tuple
import Data.Traversable

import CCEA.EA

xorVars :: NNVars
xorVars = nnVars

nnVars :: NNVars
nnVars = Map.fromList [("numberInputs", 2)
                      ,("numberHidden", 2)
                      ,("numberOutputs",1)
                      ,("timesToTrain",1)
                      ,("learningRate",0.25)
                      ,("sigmoidOrTanh",0)
                      ,("randomLowerBound", (-2))
                      ,("randomUpperBound", 2)
                      ,("gaussianMean", 0)
                      ,("gaussianStdDev", 1)
                      ,("mutationRate", 0.5)]

xorFitness :: NN n => FitnessFunction n
xorFitness = xorFitnessFunction

xorFitnessFunction :: NN n => n -> Double
xorFitnessFunction net = let x1  = (get net $ fromList [0.0,0.0]) ! 0
                             x2  = (get net $ fromList [0.0,1.0]) ! 0
                             x3  = (get net $ fromList [1.0,0.0]) ! 0
                             x4  = (get net $ fromList [1.0,1.0]) ! 0
                             x1e = 1 - abs (0 - x1)
                             x2e = 1 - abs (1 - x2)
                             x3e = 1 - abs (1 - x3)
                             x4e = 1 - abs (0 - x4)
                         in (x1e + x2e + x3e + x4e) / 4

xorBreeding :: (RandomGen g, NN n) => BreedingStrategy n g
xorBreeding g = notElitist xorBreedingStrategy xorBreedingStrategy g

xorBreedingStrategy ::  (RandomGen g, NN n) => BreedingStrategy n g
xorBreedingStrategy = mutateWeights

xorEA :: RandomGen g => IO (EA Net g)
xorEA = do
  pop <- sequence $ map (\_ -> create nnVars) [1..100]
  return $ EA pop (tournament 2) xorFitness xorBreeding

max_stars_repo_name: eklinkhammer/neural-algorithms

max_issues_repo_name: eklinkhammer/neural-algorithms

max_forks_repo_name: eklinkhammer/neural-algorithms


text:
{-# LANGUAGE DeriveFunctor, GeneralizedNewtypeDeriving #-}
{-# LANGUAGE DataKinds, TypeFamilies #-}
{-# LANGUAGE FlexibleContexts, FlexibleInstances #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE ConstraintKinds #-}

module DelayedSampling where

import Prelude hiding ((<>))

import Control.Exception.Base (assert)
import Control.Monad (forM_, when)
import Control.Monad.State (get, put)
import Control.Monad.Trans (MonadIO, liftIO)
import Control.Monad.Bayes.Class (MonadSample, MonadCond, MonadInfer)
import Control.Monad.Bayes.Weighted
import Control.Monad.Bayes.Sampler

import Data.Aeson (ToJSON (..), object, (.=))
import Data.List (delete)
import Data.Maybe (isJust)
import Data.Proxy
import Data.Random hiding (normal)

import Debug.Trace

import GHC.TypeLits

import Numeric.AD
import Numeric.AD.Rank1.Tower (Tower)

import Numeric.LinearAlgebra.Static hiding (M)

import Unsafe.Coerce (unsafeCoerce)

import Distributions
import MVDistributions
import Inference
import Util.Ref

data MarginalT = MGaussianT | MBetaT | MBernoulliT | forall (n :: Nat). MMVGaussianT (Proxy n)

type DelayedSampling = MonadState Heap
type DelayedSample m = (MonadState Heap m, MonadSample m)
type DelayedInfer m = (MonadState Heap m, MonadInfer m)

data SMarginalT (m :: MarginalT) where
  SMGaussianT :: SMarginalT MGaussianT
  SMBetaT :: SMarginalT MBetaT
  SMBernoulliT :: SMarginalT MBernoulliT
  SMMVGaussianT :: forall (n :: Nat) (p :: Proxy n). SMarginalT (MMVGaussianT p)

type family MType (m :: MarginalT) where
  MType MGaussianT = Double
  MType MBetaT = Double
  MType MBernoulliT = Bool
  MType (MMVGaussianT (p :: Proxy n)) = R n

data MDistr (m :: MarginalT) where
  MGaussian :: !Double -> !Double -> MDistr MGaussianT
  MBeta :: !Double -> !Double -> MDistr MBetaT
  MBernoulli :: !Double -> MDistr MBernoulliT
  MMVGaussian :: forall (n :: Nat) (p :: Proxy n). KnownNat n => !(R n) -> !(Sym n) -> MDistr (MMVGaussianT p)

deriving instance Show (MDistr m)

instance ToJSON (MDistr a) where
  toJSON (MMVGaussian mu cov) = object ["mean" .= mu, "cov" .= cov]
  toJSON (MGaussian mu cov) = object ["mean" .= mu, "cov" .= cov]
  toJSON (MBeta a b) = object ["alpha" .= a, "beta" .= b]
  toJSON (MBernoulli p) = object ["p" .= p]

data CDistr (m :: MarginalT) (m' :: MarginalT) where
  AffineMeanGaussian :: !Double -> !Double -> !Double -> CDistr MGaussianT MGaussianT
  CBernoulli :: CDistr MBetaT MBernoulliT
  MVAffineMeanGaussian :: forall (n :: Nat) (m :: Nat) (pn :: Proxy n) (pm :: Proxy m).
    (KnownNat n, KnownNat m) => (L m n) -> (R m) -> (Sym m) -> CDistr (MMVGaussianT pn) (MMVGaussianT pm)

data Result a where
  RConst :: a -> Result a
  RMarginal :: MDistr a -> Result (MType a)

deriving instance Show a => Show (Result a)

gaussianConditioning :: Double -> Double -> Double -> Double -> (Double, Double)
gaussianConditioning mu var obs obsvar = (mu', var')
  where
  ivar = recip var
  iobsvar = recip obsvar
  inf = ivar + iobsvar
  var' = recip inf
  mu' = (ivar * mu + iobsvar * obs) / inf

gaussianMeanGaussian :: Double -> CDistr MGaussianT MGaussianT
gaussianMeanGaussian = AffineMeanGaussian 1 0

ascend1 :: Floating d => d -> (d -> d) -> Int -> d -> d
ascend1 lr diff = go where
  go 0 x = x
  go n x = let d = diff x in go (n - 1) (x + (lr * d))

-- XXX: not sure that the factor amount is good
solveLaplaceGaussian :: MDistr MGaussianT
  -> (forall s. AD s (Tower Double) -> AD s (Tower Double))
  -> Double
  -> Int
  -> (MDistr MGaussianT, Double)
solveLaplaceGaussian (MGaussian mu var) likelihood learningRate numIters =
  (MGaussian mu' var', factorAmt - originalFactorAmt)
  where
  fmu' : _ : f''mu' : _ = diffs0 f mu'
  originalFactorAmt = log $ sqrt (pi / var) * exp (gaussian_ll' mu var mu)
  factorAmt = log $ sqrt (- 2 * pi / f''mu') * exp fmu'
  var' :: Double
  var' = -2 / f''mu'
  mu' ::  Double
  mu' = ascend1 learningRate ((!! 1) . diffs0 f) numIters mu
  f :: forall s. AD s (Tower Double) -> AD s (Tower Double)
  f x = gaussian_ll' (auto mu) (auto var) x + likelihood x

mdistrToDistr :: MDistr a -> Distr (MType a)
mdistrToDistr (MGaussian mu var) = normal mu var
mdistrToDistr (MMVGaussian mu var) = mvNormal mu var
mdistrToDistr (MBeta a b) = beta a b
mdistrToDistr (MBernoulli p) = bernoulli p

cdistrToDistr :: CDistr m m' -> MType m -> MDistr m'
cdistrToDistr (AffineMeanGaussian m b obsvar) mu = MGaussian (m * mu + b) obsvar
cdistrToDistr (MVAffineMeanGaussian f b obsvar) mu = MMVGaussian (f #> mu + b) obsvar
cdistrToDistr CBernoulli p = MBernoulli p

makeMarginal :: MDistr a -> CDistr a b -> MDistr b
makeMarginal (MGaussian mu var) (AffineMeanGaussian m b obsvar) =
  MGaussian (m * mu + b) (m^2 * var + obsvar)
makeMarginal (MMVGaussian mu var) (MVAffineMeanGaussian f b obsvar) =
  MMVGaussian (f #> mu + b) (conjugate f var + obsvar)
makeMarginal (MBeta a b) CBernoulli = MBernoulli (a / (a + b))
makeMarginal (MBernoulli _) _ = error "impossible"

makeConditional :: MDistr a -> CDistr a b -> MType b -> MDistr a
makeConditional (MGaussian mu var) (AffineMeanGaussian m b obsvar) obs =
  MGaussian mu' var'
  where (mu', var') = gaussianConditioning mu var ((obs - b) / m) (obsvar / m^2)
makeConditional (MMVGaussian mu var) (MVAffineMeanGaussian f b obsvar) obs =
  MMVGaussian mu' var'
  where (mu', var') = fst (kalmanUpdate obs f obsvar (mu, var))
makeConditional (MBeta a b) CBernoulli tf =
  if tf then MBeta (a + 1) b else MBeta a (b + 1)
makeConditional (MBernoulli _) _ _ = error "impossible"

data State b where
  Initialized :: State b
  Marginalized :: MDistr b -> State b

deriving instance Show (State MGaussianT)
deriving instance Show (State MBetaT)

isMarginalized :: State b -> Bool
isMarginalized (Marginalized _) = True
isMarginalized _ = False

data DSDistr a b where
   UDistr :: !(MDistr b) -> DSDistr a b
   CDistr :: !(Ref (Node z a)) -> !(CDistr a b) -> DSDistr a b

typeOfMDistr :: MDistr b -> SMarginalT b
typeOfMDistr (MGaussian _ _) = SMGaussianT
typeOfMDistr (MBeta _ _) = SMBetaT
typeOfMDistr (MBernoulli _) = SMBernoulliT
typeOfMDistr (MMVGaussian _ _) = SMMVGaussianT

typeOfCDistr :: CDistr m m' -> SMarginalT m'
typeOfCDistr (AffineMeanGaussian _ _ _) = SMGaussianT
typeOfCDistr CBernoulli = SMBernoulliT
typeOfCDistr (MVAffineMeanGaussian _ _ _) = SMMVGaussianT

typeOfDSDistr :: DSDistr a b -> SMarginalT b
typeOfDSDistr (UDistr d) = typeOfMDistr d
typeOfDSDistr (CDistr _ c) = typeOfCDistr c

data Node a b = Node
  { name :: !String
  , children :: [RefNodeFrom b]
  , state :: !(State b)
  , distr :: !(DSDistr a b)
  }
  | RealizedNode (MType b)

data RefNode a b where
  RefNode :: !(Ref (Node a b)) -> RefNode a b

data SomeRefNode where
  SomeRefNode :: !(Ref (Node a b)) -> SomeRefNode

instance Eq SomeRefNode where
  SomeRefNode x == SomeRefNode y = x == unsafeCoerce y

data Path e (a :: MarginalT) (b :: MarginalT) where
  PNil :: Path e a a
  PSnoc :: e b c -> Path e a b -> Path e a c

data PathTo e b where
  PathTo :: !(Path e a b) -> PathTo e b

data RefNodeTo b where
  RefNodeTo :: !(Ref (Node a b)) -> RefNodeTo b

data RefNodeFrom a where
  RefNodeFrom :: !(Ref (Node a b)) -> RefNodeFrom a

instance Eq (RefNodeFrom a) where
  RefNodeFrom x == RefNodeFrom y = x == unsafeCoerce y

instance Eq (RefNodeTo a) where
  RefNodeTo x == RefNodeTo y = x == unsafeCoerce y

ancestry :: DelayedSampling m => Node a b -> m (PathTo RefNode a)
ancestry n = case distr n of
  UDistr _ -> pure $ PathTo PNil
  CDistr par cdistr -> do
    PathTo ancestors <- ancestry =<< readRef par
    pure (PathTo (PSnoc (RefNode par) ancestors))

path_all :: Applicative f => (forall a b. e a b -> f Bool) -> Path e a b -> f Bool
path_all f PNil = pure True
path_all f (PSnoc e es) = (&&) <$> f e <*> path_all f es

isTerminal :: DelayedSampling m =>Node a b -> m Bool
isTerminal n = do
  PathTo ancestors <- ancestry n
  ancestors_all_marginalized <- path_all (\(RefNode n) -> isMarginalized . state <$> readRef n) ancestors
  pure $ isMarginalized (state n) && ancestors_all_marginalized

-- initialize without parent node
constant' :: String -> MDistr a -> Node z a
constant' n d = Node
  { name = n
  , children = []
  , state = Initialized
  , distr = UDistr d
  }

assumeConstant :: DelayedSampling m => String -> MDistr a -> m (Ref (Node z a))
assumeConstant n d = newRef (constant' n d)

-- initialize with parent node
newConditional' :: DelayedSampling m => String -> Ref (Node a b) -> CDistr b c -> m (Ref (Node b c))
newConditional' str par cdistr = do
  childRef <- newRef child
  modifyRef' par $ \n -> n { children = (RefNodeFrom childRef) : children n }
  --childRef' <- mkWeakRef childRef $ putStrLn ("deleted " ++ str)
  pure childRef
  where
  child = Node
    { name = str
    , children = []
    , state = Initialized
    , distr = CDistr par cdistr }

assumeConditional :: DelayedSampling m => String -> Ref (Node a b) -> CDistr b c -> m (Ref (Node b c))
assumeConditional str par cdistr = newConditional' str par cdistr

updater :: DelayedSampling m => (a -> m a) -> Ref a -> m ()
updater f ref = do
  x <- readRef ref
  y <- f x
  writeRef ref y

marginalize :: DelayedSampling m => MDistr a -> Ref (Node a b) -> m (MDistr b)
marginalize parMarginal nref = do
  n <- readRef nref
  writeLog ("marginalize " ++ name n)
  case (state n, distr n) of
    (Initialized, CDistr par cdistr) -> do
      let marginal' = makeMarginal parMarginal cdistr
      writeRef nref $ n { state = Marginalized marginal' }
      pure marginal'
    (state, _) -> error $  "marginalize': " ++ name n

marginalize2' :: MType a -> Node a b -> Node a b
marginalize2' value n =
  case (state n, distr n) of
    (Initialized, CDistr _ cdistr) ->
      let marg = cdistrToDistr cdistr value in
      n { state = Initialized, distr = UDistr marg }
    _ -> error "marginalize2'"

marginalize2 :: DelayedSampling m => MType a -> Ref (Node a b) -> m ()
marginalize2 a = updater (pure . marginalize2' a)

sample :: DelayedSampling m => MonadSample m => Ref (Node a b) -> m ()
sample nref = do
  n <- readRef nref
  ioAssert (isTerminal n)
  writeLog ("sample " ++ name n)
  case state n of
    Marginalized m -> do
      x <- Distributions.sample (mdistrToDistr m)
      realize x nref
    _ -> error "sample"

updateParent :: DelayedSampling m => Node b c -> (forall a. CDistr b c -> Node a b -> m (Node a b)) -> m ()
updateParent n f = case distr n of
  UDistr _ -> pure ()
  CDistr p cdistr -> updater (f cdistr) p

firstSatisfying :: Monad m => (a -> m Bool) -> [a] -> m (Maybe a)
firstSatisfying p [] = pure Nothing
firstSatisfying p (x : xs) = do
  px <- p x
  if px then pure (Just x) else firstSatisfying p xs

-- Invariant 2: A node always has at most one marginal Child
marginalChild :: DelayedSampling m => Node a b -> m (Maybe (RefNodeFrom b))
marginalChild n = firstSatisfying (\(RefNodeFrom x) -> isMarginalized . state <$> readRef x) (children n)

ioAssert :: Applicative m => m Bool -> m ()
ioAssert x = pure () -- do {b <- x; assert b (pure ())}

writeLog :: Applicative m => String -> m ()
writeLog = const $ pure () -- putStrLn

realize :: DelayedSampling m => MType b -> Ref (Node a b) -> m ()
realize val nref = do
  n <- readRef nref
  writeLog ("realize " ++ name n)
  ioAssert (isTerminal n)
  updateParent n $ \cdistr p -> do
    let Marginalized marg = state p
    let marg' = makeConditional marg cdistr val
    pure $ p { state = Marginalized marg', children = delete (RefNodeFrom nref) (children p) }
  forM_ (children n) $ \(RefNodeFrom c) -> do
    marginalize2 val c
  writeRef nref $ RealizedNode val

score' :: DelayedSampling m => MonadCond m => MType b -> Ref (Node a b) -> m Double
score' x nref = do
  n <- readRef nref
  ioAssert (isTerminal n)
  writeLog ("observe " ++ name n)
  pure $ case state n of
    Marginalized marg -> Distributions.score (mdistrToDistr marg) x
    _ -> error "observe'"

observe :: DelayedSampling m => MonadCond m => MType b -> Ref (Node a b) -> m ()
observe x nref = do
  ll <- score' x nref
  factor ll
  realize x nref

prune :: DelayedSampling m => MonadSample m => Ref (Node a b) -> m ()
prune nref = do
  n <- readRef nref
  writeLog ("prune " ++ name n)
  assert (isMarginalized (state n)) $ do
    maybeChild <- marginalChild n
    forM_ maybeChild $ \(RefNodeFrom c) ->
      prune c
    DelayedSampling.sample nref

-- turns `nref` into a terminal node
graft :: DelayedSampling m => MonadSample m => Ref (Node a b) -> m (MDistr b)
graft nref = do
  n <- readRef nref
  writeLog ("graft " ++ name n)
  case state n of
    Marginalized _ -> do
      maybeChild <- marginalChild n
      forM_ maybeChild $ \(RefNodeFrom c) -> prune c
      n'  <- readRef nref
      pure $ case state n' of
        Marginalized m -> m
        _ -> error "graft"
    Initialized -> case distr n of
      UDistr d -> do
        writeRef nref $ n { state = Marginalized d }
        pure d
      CDistr par cdistr -> do
        parMarg <- graft par
        marginalize parMarg nref

obs :: DelayedSampling m => MonadInfer m => MType b -> Ref (Node a b) -> m ()
obs x n = do
  graft n
  DelayedSampling.observe x n

score :: DelayedSampling m => MonadInfer m => MType b -> Ref (Node a b) -> m Double
score x n = withoutModifying $ do
  graft n
  score' x n

withoutModifying :: MonadState s m => m a -> m a
withoutModifying f = do
  s <- get
  x <- f
  put s
  return x

getValue :: DelayedSampling m => MonadSample m => Ref (Node a b) -> m (MType b)
getValue nref = do
  n <- readRef nref
  case n of
    RealizedNode x -> pure x
    _ -> do
      graft nref
      DelayedSampling.sample nref
      getValue nref

printValue :: Show (MType b) => DelayedSampling m => MonadSample m => MonadIO m => Ref (Node a b) -> m ()
printValue n = do
  x <- getValue n
  liftIO $ print x

printState :: Show (State b) => DelayedSampling m => MonadIO m => Ref (Node a b) -> m ()
printState nref = do
  n <- readRef nref
  liftIO $ print (state n)

forget :: DelayedSampling m => Ref (Node a b) -> m ()
forget nref = do
  n <- readRef nref
  case n of
    RealizedNode _ -> pure ()
    _ -> case state n of
      Initialized -> error "forget"
      Marginalized marg -> do
        forM_ (children n) $ \(RefNodeFrom cref) -> do
          modifyRef' cref $ \c ->
            c { distr = case distr c of
              UDistr d -> UDistr d
              CDistr _ _ -> case state c of
                Marginalized marg -> UDistr marg
                _ -> error "forget" }
        case distr n of
          UDistr d -> pure ()
          CDistr par cdistr -> do
            isStale <- stale par
            if isStale
              then error "can't forget: information to incorporate"
              else updateParent n $ \cdistr p -> pure $ p { children = delete (RefNodeFrom nref) (children p) }
        writeRef nref $ n { distr = UDistr marg }
  modifyRef' nref $ \n' -> n' { children = [] }
  freeRef nref

-- sweep :: [SomeRefNode] -> M ()
-- sweep liveVars = mapM_ (\(SomeRefNode n) -> forgetRelated n) liveVars where
--   forgetRelated :: Ref (Node a b) -> M ()
--   forgetRelated nref = undefined
--   forgettable :: Ref (Node a b) -> M Bool
--   forgettable nref =
--     if SomeRefNode nref `elem` liveVars
--       then pure False
--       else case state n of

stale :: DelayedSampling m => Ref (Node a b) -> m Bool
stale nref = do
  n <- readRef nref
  case n of
    RealizedNode _ -> pure False
    _ -> case state n of
      Initialized -> case distr n of
        UDistr _ -> pure False
        CDistr par _ -> stale par
      Marginalized d -> isJust <$> marginalChild n

-- Examples

delay_triplet :: DelayedInfer m => MonadIO m => Double -> m ()
delay_triplet zobs = do
  x <- assumeConstant "x" (MGaussian 0 1)
  writeLog "x created"
  y <- assumeConditional "y"  x (gaussianMeanGaussian 1)
  writeLog "y created"
  z <- assumeConditional "z" y (gaussianMeanGaussian 1)
  writeLog "z created"
  obs zobs z
  writeLog "z observed"
  printValue z
  printValue x
  printValue y

observeConditional :: DelayedInfer m => String -> Ref (Node a b) -> CDistr b c -> MType c -> m ()
observeConditional str nref cdistr observation = do
  y <- assumeConditional str nref cdistr
  obs observation y
  freeRef y

scoreConditional :: DelayedInfer m => String -> Ref (Node a b) -> CDistr b c -> MType c -> m Double
scoreConditional str nref cdistr observation = do
  y <- assumeConditional str nref cdistr
  ll <- DelayedSampling.score observation y
  modifyRef' nref $ \n -> n { children = tail (children n) }
  freeRef y
  pure ll

-- For some reason, likelihoods aren't exactly  the same
delay_iidADF :: DelayedInfer m => MonadIO m => Bool -> [Double] -> m ()
delay_iidADF useADF yobs = do
  x <- assumeConstant "x" (MGaussian 0 1)
  forM_ (zip [1..] yobs) $ \(t, obsyt) -> do
    if useADF
      then observeGeneric x $ \d -> solveLaplaceGaussian d (\mu -> gaussian_ll' mu 1 (auto obsyt)) 0.05 10000
      else observeConditional ("y" ++ show t) x (gaussianMeanGaussian 1) obsyt  -- normal(x, 1)
    printState x
  printValue x

observeGeneric :: DelayedInfer m => Ref (Node a b) -> (MDistr b -> (MDistr b, Double)) -> m ()
observeGeneric nref f = do
  graft nref
  n <- readRef nref
  case n of
    RealizedNode x -> error "observeGeneric: RealizedNode impossible"
    _ -> case state n of
      Marginalized d -> do
        let (d', factorAmt) = f d
        factor factorAmt
        writeRef nref $ n { state = Marginalized d' }
      _ -> error "observeGeneric"

delay_kalman :: DelayedInfer m => MonadIO m =>
  Bool -> Double -> Double -> [Double] -> m ()
delay_kalman shouldForget m b yobs = do
  x1 <- assumeConstant "x1" (MGaussian 0 1)
  observeConditional "y1" x1 (gaussianMeanGaussian 1) (head yobs)
  loop 2 x1 (tail yobs)
  where
  loop :: DelayedInfer m => MonadIO m =>
    Integer -> Ref (Node a MGaussianT) -> [Double] -> m ()
  loop t xpredt [] = pure ()
  loop t xpredt (y : ys) = do
    xt <- assumeConditional ("x" ++ show t) xpredt (AffineMeanGaussian m b 1)
    observeConditional ("y" ++ show t) xt (gaussianMeanGaussian 1) y
    when (t `mod` 100 == 0) $ do
      liftIO $ putStrLn ("t = " ++ show t)
      printState xt
    when shouldForget $ forget xpredt
    loop (t + 1) xt ys

max_stars_repo_name: psg-mit/probzelus-haskell

max_issues_repo_name: psg-mit/probzelus-haskell

max_forks_repo_name: psg-mit/probzelus-haskell


text:
module ReadWrite
  ( -- * Read
    readMolV2000
    -- * Write
  , writeXYZ
  ) where

import Prelude hiding (readFile)

import Control.Lens
import Control.Monad.State
import Numeric.LinearAlgebra.Data
import System.Directory
import System.IO
import System.IO.Unsafe
import Text.Printf (hPrintf)

import Types

-- * Read
-- | Read molecule at *.pdb format (V2000)
readMolV2000 :: FilePath -> (Molecule, [Bond])
readMolV2000 inf =
  let txt = (lines . unsafePerformIO . readFile) inf
      counts_line = words (txt !! 3)
      count_atoms = read (counts_line !! 0)
      count_bonds = read (counts_line !! 1)
      atoms_lines = (take count_atoms . drop 4) txt
      bonds_lines = (take count_bonds . drop (4 + count_atoms)) txt
   in (foldr addAtom molecule atoms_lines, foldr addBond [] bonds_lines)
  where
    addAtom a = over atoms ((:) (readline' a))
    addBond b = (:) (readline'' b)
    readline' l =
      execState
        (do let w = words l
            acoordin . x .= read (w !! 0)
            acoordin . y .= read (w !! 1)
            acoordin . z .= read (w !! 2)
            aelement .= w !! 3
            avdwrad .= vdwr (w !! 3))
        atom
    readline'' l =
      execState
        (do let w = words l
            bfid .= read (w !! 0) - 1
            bsid .= read (w !! 1) - 1
            btype .= read (w !! 2)
            bster .= read (w !! 3)
            btop .= read (w !! 4))
        bond

-- | Get VDW radius
vdwr :: Element -> Double
vdwr a =
  case a of
    "H" -> 0.500 -- 0.5 -- 1.000
    "O" -> 0.650 -- 0.5 -- 1.300
    "N" -> 0.700 -- 0.5 -- 1.400
    "C" -> 0.750 -- 0.5 -- 1.500
    "S" -> 0.950 -- 0.5 -- 1.900
    "Br"-> 0.950
    othrewise -> error $ "vdwr not found for: " ++ show a

-- * Write
-- | Write molecule in *.xyz format
writeXYZ :: FilePath -> String -> Molecule -> IO ()
writeXYZ ouf comment molecule = do
  (tmp_name, tmp_handle) <- openTempFile "." "temp"
  hPrint tmp_handle (views atoms length molecule)
  hPutStrLn tmp_handle comment
  mapM_ (writeData tmp_handle) (view atoms molecule)
  hClose tmp_handle
  renameFile tmp_name ouf
  where
    writeData hdl atom = do
      let e = view aelement atom
          (Point x y z) = view acoordin atom
      hPrintf hdl "%s\t%8.6f\t%8.6f\t%8.6f\n" e x y z

max_stars_repo_name: wurthel/distance-geometry

max_issues_repo_name: wurthel/distance-geometry

max_forks_repo_name: wurthel/distance-geometry


text:
{-# OPTIONS_GHC  -fno-warn-unused-binds -fno-warn-unused-matches -fno-warn-name-shadowing -fno-warn-missing-signatures #-}
{-# LANGUAGE FlexibleInstances, ConstraintKinds, ExistentialQuantification, GADTs, RankNTypes, MultiParamTypeClasses, RankNTypes, UndecidableInstances, FlexibleContexts, TypeSynonymInstances #-}

---------------------------------------------------------------------------------------------------
---------------------------------------------------------------------------------------------------
-- |
-- | Module : First attempt at approximate counting
-- | Creator: Xiao Ling
-- | Created: 12/08/2015
-- | TODO   : test standard deviation of alpha, beta, and final version
-- |
---------------------------------------------------------------------------------------------------
---------------------------------------------------------------------------------------------------

module Morris (

    morris
  , morris'

  , morris1
  , morris1'

  ) where

import Control.Monad.Random
import Control.Monad.Random.Class

import Data.Conduit
import Data.List.Split
import qualified Data.Conduit.List as Cl

import Core
import Statistics

{-----------------------------------------------------------------------------
    I. Morris Algorithm list of counter
------------------------------------------------------------------------------}

---- * Count the number of items in `as` to within `eps` of actual
---- * with confidence `delta`
morris' :: EpsDelta -> Batch a IO Counter
morris' ed = morris ed `using` evalRandIO

morris :: MonadRandom m => EpsDelta -> Streaming a m Counter
morris (ED e d) = medianOfMeans $ go cs
  where cs            = replicate (t*m) 0
        t             = round $ 1/(e^2*d)  :: Int
        m             = round . log $ 1/d
        go cs         = (\c -> 2^(round c) - 1) `ffmap` Cl.foldM (\cs _ -> traverse incr cs) cs
        ffmap         = fmap . fmap
        medianOfMeans = fmap median' . (fmap . fmap) mean' . fmap (chunksOf t)

{-----------------------------------------------------------------------------
    II. Morris Algorithm list of list of counter
------------------------------------------------------------------------------}

---- * Count the number of items in `as` to within `eps` of actual
---- * with confidence `delta`
morris1' :: EpsDelta -> Batch a IO Counter
morris1' ed = morris1 ed `using` evalRandIO

-- * Run on stream inputs `as` for t independent trials for `t = 1/eps^2 * d`,
-- * and `m` times in parallel, for `m = log(1/d)` and take the median
morris1 :: MonadRandom m => EpsDelta -> Streaming a m Counter
morris1 (ED e d) = medianOfMeans $ go ccs
      where
        medianOfMeans = fmap median' . (fmap . fmap) mean'
        ccs           = replicate m $ replicate t 0
        t             = round $ 1/(e^2*d)
        m             = round . log $ 1/d
        go ccs        = (\x -> 2^(round x) - 1) `fffmap` Cl.foldM (\xs _ -> incrs' xs) ccs
        fffmap        = fmap . fmap . fmap

-- * given a list of list of counters toss a coin for each counter and incr
-- * this can be flattened
incrs' :: MonadRandom m => [[Counter]] -> m [[Counter]]
incrs' = sequence . fmap (sequence . fmap incr)

{-----------------------------------------------------------------------------
    III. Utils
------------------------------------------------------------------------------}

-- * Increment a counter `x` with probability 1/2^x
incr :: MonadRandom m => Counter -> m Counter
incr x = do
  h <- toss . coin $ 0.5^(round x)
  return $ if isHead h then (seq () succ x) else seq () x

mean', median' :: (Floating a, Ord a, RealFrac a) => [a] -> Float
mean'   = fromIntegral . round . mean
median' = fromIntegral . round . median


text:
module NeuralNetwork.Types where

import Numeric.LinearAlgebra.Data (R, Vector)
import Numeric.Natural (Natural)

newtype Layer = Layer
  { numUnits :: Natural
  } deriving (Eq, Show)

data Network = Network
  { inputLayer :: Layer
  , hiddenLayer :: [Layer]
  , outputLayer :: Layer
  } deriving (Eq, Show)

max_stars_repo_name: sebashack/mlToolBox

max_issues_repo_name: sebashack/mlToolBox

max_forks_repo_name: sebashack/mlToolBox


text:
import Data.Map.Strict ((!))
import Numeric.LinearAlgebra hiding ((!))
import LeastSquares

leastSquaresEstimate :: Vector R
leastSquaresEstimate = result ! "x"
  where result = minimize (SumSquares (a*x-b)) [c*x:==:d]
        a = Mat $ (3><2) [1.0, 1.0, 1.0, 4.0, 2.0, -1.0]
        b = Vec $ vector [13.0, 27.0, 1.0]
        c = Mat $ (1><2) [3.0, 7.0]
        d = Vec $ vector [-1.0]
        x = Var "x" 2

max_stars_repo_name: tridao/leastsquares

max_issues_repo_name: tridao/leastsquares

max_forks_repo_name: tridao/leastsquares


text:
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE FunctionalDependencies #-}
{-# LANGUAGE TupleSections #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE ViewPatterns #-}
module Algorithm.Training where

import Prelude hiding (foldl, foldr)
-- import Numeric.AD
import Numeric.LinearAlgebra as LA hiding (fromList, (<>))
import Numeric.LinearAlgebra.Devel
import qualified Data.Vector.Storable as V hiding (foldr, modify)
import Data.Functor
import Data.Bifunctor
import Data.Monoid as M
import Data.Foldable as F
import Data.Ord
import Control.Monad
import Pipes
import Pipes.Prelude as P
import Control.DeepSeq
-- import Control.Monad.RWS
-- import Control.Monad.RWS as RWS
import Data.Network.Activation
import Data.Network
import Data.Network.Utils
import Algorithm.GradientDescent

-- type TrainingSet = [Sample]

-- type TestSet = [Sample]

-- type ResultSet = [(Output, Output)]

-- Loss function parameterized by sample type (a), output loss type(l) and parameter space (p)
data LossFunction a l = Evaluator { evaluate :: (a, a) -> l, evaluate' :: (a, a) -> l }

instance Differentiable (a,a) l (LossFunction a l) where
  derivate = evaluate'

-- The type of a loss/error evaluation
type Loss = Vector R
type LossBatch = Matrix R

-- y is target, y' is actual
mse :: (Floating a) => LossFunction a a
mse = let f (y,y') = let γ = y'-y in (/ 2) $ γ**2
          f' (y,y') = (y'-y)
      in  Evaluator f f'

mseBatch :: LossFunction Batch Batch
mseBatch = let fb (b, b')  = (/ 2) $ (**2) (b' - b)
               fb' (b, b') = (b' - b)
           in Evaluator fb fb'

type NetworkEvaluator a = LossFunction a a

type LearningRate = Double
-- Training parameters = (epochs, learning rate)
type TrainingParameters = (Int, Double)

class (Network i n, Differentiable (i, i) i lf, Monad m) => Backpropagation i lf n m | lf -> i, n -> i where
  -- given a neural network and an input, compute each layer's activation and derivative
  -- forward :: n -> i -> m [(i,i)]
  -- given a loss function, a network, an output and target pair and a list of activation derivatives, compute the deltas
  backward :: lf -> n -> (i, i) -> [i] -> m [i]

data TrainingState m  = TS {
  tsmodel :: m,      -- the model(network) state
  tseta   :: Double, -- the learning rate η
  tsiteration :: Int,-- the iteration number(one per gradient move)
  tsloss  :: Double,  -- the loss of the model over the last training batch/sample
  tssampleid :: Int
  } deriving Show

-- type DescentStopper i n = StopCondition (LossFunction i i) (n, (i,i))

class AsDouble a where
  toDouble :: a -> Double

instance AsDouble (Vector R) where
  toDouble = norm_2

instance AsDouble (Matrix R) where
  toDouble m = (F.sum . (fmap norm_2) . toRows) m / (fromIntegral $ rows m)

-- A pipe that takes a sample and makes one step in the parameter space
stepP :: (Monad m, AsDouble i, Network i n, GradientDescent (LossFunction i i) (n, (i,i)) m) =>
  n -> LossFunction i i -> Int -> (Int -> Double) -> (Int, (i,i)) -> Producer (TrainingState n) m r
stepP m lf n lrf (sid, (i, t)) = loop m n where
  loop m n = do
    let !η = lrf n -- compute learning rate for iteration
    (m', _) <- lift $ descend lf η (m, (i,t)) -- step
    let o = fst $ runNetwork m' i -- test
        loss = evaluate lf (t,o) -- evaluate
    m' `deepseq` yield $ TS m' η n (toDouble loss) sid -- yield new state
    loop m' (n+1)

-- A pipe that produce a reporting action on a training state
reportP :: (Monad m, AsDouble i, Network i n) => Int -> (TrainingState n -> m ()) -> Pipe (TrainingState n) (TrainingState n) m r
reportP interval action = P.mapM $ \ts -> do
  when (tsiteration ts `mod` interval == 0) (action ts)
  return ts

type Stopper n x m = (TrainingState n, TrainingState n) -> m (Maybe x)

-- Does gradient descent given a stop condition and a stream of "steps"(training states containing an updated model)
descendP :: (Monad m, AsDouble i, Network i n) =>
    Stopper n x m-> Consumer (TrainingState n) m x
descendP stopper = loop
  where loop = do
          ts <- await
          ts' <- await
          mx <- lift (stopper (ts, ts'))
          case mx of
            Just x -> return x
            Nothing -> loop

-- A stop condition: stop if the loss doesn't change significantly enough,
-- returning the state with the minimum loss
closeEnough :: (Network i n, Monad m) => Double -> Stopper n (TrainingState n) m
closeEnough tolerance (ts, ts') = do
  if abs (tsloss ts' - tsloss ts) < tolerance then
    return (Just $ minimumBy (comparing tsloss) [ts, ts']) else return Nothing

-- stops after a given number of iterations
longEnough :: (Network i n, Monad m) => Int -> Stopper n (TrainingState n) m
longEnough iteration (ts, ts') = if tsiteration ts' >= iteration then
  return (Just $ minimumBy (comparing tsloss) [ts, ts']) else return Nothing

-- stops when there's not enough improvement(i.e. loss reduction) between the two states
improvement :: (Network i n, Monad m) => Double -> Stopper n (TrainingState n) m
improvement tolerance (ts, ts') = if ((tsloss ts' - tsloss ts) > tolerance) then
  return (Just $ minimumBy (comparing tsloss) [ts, ts']) else return Nothing

max_stars_repo_name: DrPyser/NeuralNetwork

max_issues_repo_name: DrPyser/NeuralNetwork

max_forks_repo_name: DrPyser/NeuralNetwork


text:
module QIO.Circuit.Compile where

import qualified QIO.Circuit.Circuit as C

import Data.Complex
import QIO.QioSyn
import QIO.Qio

compileCircuit :: C.Circuit -> QIO [Bool]
compileCircuit circuit = do
  qs <- mkQbits (length $ C.qubits circuit) False
  applyU (gatesToU (C.gates circuit) qs)
  measQbits qs

gatesToU :: [C.Gate] -> [Qbit] -> U
gatesToU [] _      = mempty
gatesToU (g:gs) qs = (gateToU g qs) <> (gatesToU gs qs)

gateToU :: C.Gate -> [Qbit] -> U
gateToU (C.Had _ qis) qs         = condU uhad $ map (\i -> qs !! i) qis
gateToU (C.PX _ qis) qs          = condU unot $ map (\i -> qs !! i) qis
gateToU (C.PY _ qis) qs          = condU (\q -> rot q yRot) $ map (\i -> qs !! i) qis
                                   where
                                     yRot (False, True) = 0 :+ 1
                                     yRot (True, False) = 0 :+ (-1)
                                     yRot _             = 0
gateToU (C.PZ _ qis) qs          = condU (\q -> uphase q pi) $ map (\i -> qs !! i) qis
gateToU (C.Swap _ (q1:q2:[])) qs = swap (qs !! q1) (qs !! q2)
gateToU _ _                      = mempty

mkQbits :: Int -> Bool -> QIO [Qbit]
mkQbits n b = mkQbits' n b []
  where
    mkQbits' 0 _ qs = return qs
    mkQbits' n b qs = do
      q <- mkQbit b
      mkQbits' (n-1) b (qs ++ [q])

measQbits :: [Qbit] -> QIO [Bool]
measQbits qs = measQbits' qs []
  where
    measQbits' [] bs = return bs
    measQbits' (q:qs) bs = do
      b <- measQbit q
      measQbits' qs (bs ++ [b])

condU :: (Qbit -> U) -> [Qbit] -> U
condU u (ql:[]) = u ql
condU u (q:qs)  = cond q (\x -> if x then (condU u qs) else mempty)

max_stars_repo_name: psycjw/qio

max_issues_repo_name: psycjw/qio

max_forks_repo_name: psycjw/qio


text:
{-# LANGUAGE GADTs #-}
module Distrib (
    -- * Data Constructors
    Distrib(..)
    -- * Operations
  , compNormalDistr
  , pdf
  ) where

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

import           Classifier
import qualified Data.Vector.Unboxed            as VU
import           Statistics.Distribution
import           Statistics.Distribution.Normal
import qualified Statistics.Sample              as Statistics

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

data Distrib where
  Distrib :: ContDistr d => d -> Distrib

-- | We don't know much about our numeric variable
-- so we will assume normallity to simplify the implementation.
compNormalDistr :: [Double] -> Distrib
compNormalDistr xs =
  Distrib (normalDistr mean' stdDev')
    where
      v = VU.fromList xs
      mean' = Statistics.mean v
      stdDev'  = max 0.1 (Statistics.stdDev v) -- stdDev must be > 0

-- | Probability Density Function
pdf :: Distrib -> Double -> Probability
pdf (Distrib d) = Probability . density d

max_stars_repo_name: monadplus/naive-bayes-classifier

max_issues_repo_name: monadplus/naive-bayes-classifier

max_forks_repo_name: monadplus/naive-bayes-classifier


text:
{-# LANGUAGE RankNTypes, BangPatterns, ScopedTypeVariables, FlexibleInstances #-}

module Math.Probably.HamMC where

import Math.Probably.MCMC
import qualified Math.Probably.PDF as PDF
import Math.Probably.RandIO
import Math.Probably.FoldingStats
import Math.Probably.Sampler
import Control.Applicative

import Numeric.LinearAlgebra
import Numeric.AD
import Text.Printf
import System.IO
import Data.Maybe

import Statistics.Test.KolmogorovSmirnov
import Statistics.Test.MannWhitneyU
import qualified Data.Vector.Unboxed as U
import qualified Control.Monad.State.Strict as S

import Debug.Trace
import Data.IORef
import Control.Spoon

data HMCPar = HMCPar { hpXi :: !(Vector Double),
                       hpL :: !Int,
                       hpPi :: !Double,
                       hpEpsilon :: !Double,
                       hpCount :: !Int,
                       hpAccept :: !Int,
                       hpFreezeEps :: Bool }
    deriving Show

hmc1 :: (Vector Double -> (Double,Vector Double))
         -> HMCPar
         -> Sampler HMCPar
hmc1 postGrad (HMCPar current_q l _ epsMean count accept freeze) = do
   let dims = dim current_q
       grad_u = negate . snd . postGrad
       u = negate . fst . postGrad
   current_p <- fmap fromList $ gaussManyUnitD dims
   eps <- uniform (epsMean*0.8) (epsMean*1.2)
   --the initial half momentum step

   let p_half_step = current_p - scale (eps/2) (grad_u current_q)

   let step :: Int -> Vector Double -> Vector Double -> (Vector Double, Vector Double)
       step n p q  -- note that my l is Radford Neal's L+1
        | n == 0 = let qfinal = q + scale eps p -- only q not p update in last loop
                       pfinal = p - scale (eps/2) (grad_u qfinal) -- half momentum step at end
                   in (negate pfinal, -- negate momentum for symmetric proposal
                       qfinal)
        | otherwise =
            let q1 = q + scale eps p
                p1 = p - scale eps (grad_u q1)
            in step (n-1) p1 q1

   let (propose_p, propose_q) = step l p_half_step current_q
       current_U =  u current_q
       current_K =  (current_p `dot` current_p) / 2
       propose_U = u propose_q
       propose_K =  (propose_p `dot` propose_p) / 2
       ratio = exp $ current_U - propose_U + current_K - propose_K
       tr = max 1.0 $ realToFrac count

   u <- unitSample -- 0 to 1
   return $ if  u < ratio
               then HMCPar propose_q l propose_U (if freeze then epsMean else (min 2 $ 1+k_hmc/tr)*epsMean) (count+1) (accept+1) freeze
               else HMCPar current_q l current_U (if freeze then epsMean else (max 0.5 $ 1-k_hmc/tr)**1.8*epsMean) (count+1) (accept) freeze

traceAs s x = trace (s++" = "++show x) x

runHMC postgrad nsam hmc0 = go nsam hmc0 [] where
  go 0 hmcp xs = do io $ putStrLn $ "done with "++show hmcp
                    return (hmcp, xs)
  go n hmcp xs = do hmcnext <- sample $ hmc1 postgrad hmcp
                    if isNaN $ hpPi hmcnext
                       then return (hmcnext, [])
                       else do io $ do putStrLn $ show $ (hpCount hmcnext,
                                                          hpPi hmcnext,
                                                          hpEpsilon hmcnext)
                                       hFlush stdout
                               go (n-1) hmcnext $ hpXi hmcnext : xs

k_hmc = 2

max_stars_repo_name: glutamate/probably

max_issues_repo_name: glutamate/probably

max_forks_repo_name: glutamate/probably