text:
{-# LANGUAGE FlexibleInstances #-}

module Isumi.Math.AppxEq
  ( AppxEq(..)
  , (~/=)
  ) where

import           Data.Complex
import qualified Data.Vector.Unboxed as UV

class AppxEq a where
  (~=) :: a -> a -> Bool

infix 4 ~=

infix 4 ~/=
(~/=) :: AppxEq a => a -> a -> Bool
x ~/= y = not $ x ~= y

instance AppxEq Int where
  x ~= y = x == y

instance AppxEq Double where
  x ~= y = abs (x - y) < 0.000001

instance AppxEq a => AppxEq (Complex a) where
  (r :+ v) ~= (r2 :+ v2) = (r ~= r2) && (v ~= v2)

instance AppxEq a => AppxEq (Maybe a) where
  Just x ~= Just y = x ~= y
  Nothing ~= Nothing = True
  _ ~= _ = False

instance AppxEq a => AppxEq [a] where
  xs ~= ys = and $ zipWith (~=) xs ys

instance (UV.Unbox a, AppxEq a) => AppxEq (UV.Vector a) where
  xs ~= ys = UV.toList xs ~= UV.toList ys

instance (AppxEq a, AppxEq b) => AppxEq (a, b) where
  (x1, y1) ~= (x2, y2) = x1 ~= x2 && y1 ~= y2

ext: hs

max_stars_repo_name: IsumiF/fft-bruun

max_issues_repo_name: IsumiF/fft-bruun

max_forks_repo_name: IsumiF/fft-bruun


text:
{-# LANGUAGE TypeFamilies, GeneralizedNewtypeDeriving, DeriveGeneric, TupleSections, RecordWildCards, TemplateHaskell, RankNTypes, FlexibleContexts #-}

module BayesStack.Models.Topic.CitationInfluenceNoTopics
  ( -- * Primitives
    NetData
  , dHypers, dArcs, dItems, dNodeItems, dCitingNodes, dCitedNodes
  , netData
  , HyperParams(..)
  , MState
  , stGammas, stOmegas, stPsis, stCiting, stLambdas
  , CitingUpdateUnit
  , ItemSource(..)
  , CitedNode(..), CitedNodeItem(..)
  , CitingNode(..), CitingNodeItem(..)
  , Citing(..), Cited(..)
  , Item(..), Topic(..), NodeItem(..), Node(..), Arc(..)
  , setupNodeItems
    -- * Initialization
  , verifyNetData, cleanNetData
  , ModelInit
  , randomInitialize
  , model
  , updateUnits
    -- * Diagnostics
  , modelLikelihood
  ) where

import qualified Data.Vector as V
import Statistics.Sample (mean)

import           Prelude hiding (mapM_, sum)
import           Data.Maybe (fromMaybe)

import           Control.Lens hiding (Setting)
import           Data.Set (Set)
import qualified Data.Set as S

import           Data.Map.Strict (Map)
import qualified Data.Map.Strict as M

import           Data.Foldable hiding (product)
import           Control.Applicative ((<$>), (<*>))
import           Control.Monad (when)
import           Control.Monad.Trans.State.Strict
import           Control.Monad.Trans.Writer.Strict

import           Data.Random
import           Data.Random.Lift (lift)
import           Data.Random.Distribution.Categorical (categorical)
import           Numeric.Log hiding (sum)

import           BayesStack.Types
import           BayesStack.Gibbs
import           BayesStack.DirMulti
import           BayesStack.TupleEnum ()
import           BayesStack.Models.Topic.Types

import           GHC.Generics (Generic)
import           Data.Binary (Binary)
import           Control.DeepSeq

at' :: At m => Index m -> IndexedLens' (Index m) m (IxValue m)
at' i = at i . _fromMaybe
  where _fromMaybe = iso (fromMaybe $ error "at': Unexpected Nothing") Just

data ItemSource = Shared | Own deriving (Show, Eq, Enum, Ord, Generic)
instance Binary ItemSource
instance NFData ItemSource

newtype Citing a = Citing a deriving (Show, Eq, Enum, Ord, Generic, NFData)
newtype Cited a = Cited a deriving (Show, Eq, Enum, Ord, Generic, NFData)
instance Binary a => Binary (Citing a)
instance Binary a => Binary (Cited a)

type CitingNode = Citing Node
type CitedNode = Cited Node
type CitingNodeItem = Citing NodeItem
type CitedNodeItem = Cited NodeItem

-- ^ A directed edge
data Arc = Arc { citingNode :: !CitingNode, citedNode :: !CitedNode }
            deriving (Show, Eq, Ord, Generic)
instance Binary Arc

data HyperParams = HyperParams
                   { _alphaPsi         :: Double
                   , _alphaLambda      :: Double
                   , _alphaOmega       :: Double
                   , _alphaGammaShared :: Double
                   , _alphaGammaOwn    :: Double
                   , _alphaBetaFG      :: Double
                   , _alphaBetaBG      :: Double
                   }
                 deriving (Show, Eq, Generic)
instance Binary HyperParams
makeLenses ''HyperParams

data NetData = NetData { _dHypers             :: !(HyperParams)
                       , _dArcs               :: !(Set Arc)
                       , _dItems              :: !(Map Item Double)
                       , _dNodeItems          :: !(Map NodeItem (Node, Item))
                       , _dCitingNodes        :: !(Map CitingNode (Set CitedNode))
                         -- ^ Maps each citing node to the set of nodes cited by it
                       , _dCitedNodes         :: !(Map CitedNode (Set CitingNode))
                         -- ^ Maps each cited node to the set of nodes citing it
                       }
              deriving (Show, Eq, Generic)
instance Binary NetData
makeLenses ''NetData

netData :: HyperParams -> Set Arc -> Map Item Double -> Map NodeItem (Node,Item) -> NetData
netData hypers arcs items nodeItems =
    NetData { _dHypers       = hypers
            , _dArcs         = arcs
            , _dItems        = items
            , _dNodeItems    = nodeItems
            , _dCitingNodes  = M.unionsWith S.union
                               $ map (\(Arc a b)->M.singleton a $ S.singleton b)
                               $ S.toList arcs
            , _dCitedNodes   = M.unionsWith S.union
                               $ map (\(Arc a b)->M.singleton b $ S.singleton a)
                               $ S.toList arcs
            }

dCitingNodeItems :: NetData -> Map CitingNodeItem (CitingNode, Item)
dCitingNodeItems nd =
    M.mapKeys Citing
    $ M.map (\(n,i)->(Citing n, i))
    $ M.filter (\(n,i)->Citing n `M.member` (nd^.dCitingNodes))
    $ nd^.dNodeItems

itemsOfCitingNode :: NetData -> CitingNode -> [Item]
itemsOfCitingNode d (Citing u) =
    map snd $ M.elems $ M.filter (\(n,_)->n==u) $ d^.dNodeItems

connectedNodes :: Set Arc -> Set Node
connectedNodes arcs =
    S.map ((\(Cited n)->n) . citedNode) arcs `S.union` S.map ((\(Citing n)->n) . citingNode) arcs

cleanNetData :: NetData -> NetData
cleanNetData d =
    let nodesWithItems = S.fromList $ map fst $ M.elems $ d^.dNodeItems
        nodesWithArcs = connectedNodes $ d^.dArcs
        keptNodes = nodesWithItems `S.intersection` nodesWithArcs
        keepArc (Arc (Citing citing) (Cited cited)) =
            citing `S.member` keptNodes && cited `S.member` keptNodes
        go = do dArcs %= S.filter keepArc
                dNodeItems %= M.filter (\(n,i)->n `S.member` keptNodes)
    in execState go d

verifyNetData :: (Node -> String) -> NetData -> [String]
verifyNetData showNode d = execWriter $ do
    let nodesWithItems = S.fromList $ map fst $ M.elems $ d^.dNodeItems
    forM_ (d^.dArcs) $ \(Arc (Citing citing) (Cited cited))->do
        when (cited `S.notMember` nodesWithItems)
            $ tell [showNode cited++" has arc yet has no items"]
        when (citing `S.notMember` nodesWithItems)
            $ tell [showNode citing++" has arc yet has no items"]

-- Citing Update unit (Shared Taste-like)
data CitingUpdateUnit = CitingUpdateUnit { _uuNI    :: CitingNodeItem
                                         , _uuN     :: CitingNode
                                         , _uuX     :: Item
                                         , _uuCites :: Set CitedNode
                                         , _uuItemWeight :: Double
                                         }
                      deriving (Show, Generic)
instance Binary CitingUpdateUnit
makeLenses ''CitingUpdateUnit

citingUpdateUnits :: NetData -> [CitingUpdateUnit]
citingUpdateUnits d =
    map (\(ni,(n,x))->CitingUpdateUnit { _uuNI      = ni
                                       , _uuN       = n
                                       , _uuX       = x
                                       , _uuCites   = d^.dCitingNodes . at' n
                                       , _uuItemWeight = (d ^. dItems . at' x)
                                       }
        ) $ M.assocs $ dCitingNodeItems d

updateUnits :: NetData -> [WrappedUpdateUnit MState]
updateUnits d = map WrappedUU (citingUpdateUnits d)

-- | Model State
data CitingSetting = OwnSetting
                   | SharedSetting !CitedNode
                   deriving (Show, Eq, Generic)
instance Binary CitingSetting
instance NFData CitingSetting where
    rnf (OwnSetting)      = ()
    rnf (SharedSetting c) = rnf c `seq` ()

data MState = MState { -- Citing model state
                       _stGammas   :: !(Map CitingNode (Multinom Int ItemSource))
                     , _stOmegas   :: !(Map CitingNode (Multinom Int Item))
                     , _stPsis     :: !(Map CitingNode (Multinom Int CitedNode))

                     , _stCiting   :: !(Map CitingNodeItem CitingSetting)

                     -- Cited model state
                     , _stLambdas  :: !(Map CitedNode (Multinom Int Item))
                     }
            deriving (Show, Generic)
instance Binary MState
makeLenses ''MState

-- | Model initialization
type ModelInit = Map CitingNodeItem (Setting CitingUpdateUnit)

modify' :: Monad m => (a -> a) -> StateT a m ()
modify' f = do x <- get
               put $! f x

randomInitializeCiting :: NetData -> ModelInit -> RVar ModelInit
randomInitializeCiting d init = execStateT doInit init
    where doInit :: StateT ModelInit RVar ()
          doInit = let unset = M.keysSet (dCitingNodeItems d) `S.difference` M.keysSet init
                   in mapM_ (randomInitCitingUU d) (S.toList unset)

randomInitCitingUU :: NetData -> CitingNodeItem -> StateT ModelInit RVar ()
randomInitCitingUU d cni@(Citing ni) =
    let (n,_) = d ^. dNodeItems . at' ni
    in case d ^. dCitingNodes . at' (Citing n) of
           a | S.null a -> do
               modify' $ M.insert cni OwnSetting

           citedNodes -> do
               s <- lift $ randomElement [Shared, Own]
               c <- lift $ randomElement $ toList citedNodes
               modify' $ M.insert cni $
                   case s of Shared -> SharedSetting c
                             Own    -> OwnSetting

randomInitialize :: NetData -> RVar ModelInit
randomInitialize d = randomInitializeCiting d M.empty

model :: NetData -> ModelInit -> MState
model d citingInit =
    let citingNodes = M.keys $ d^.dCitingNodes
        hp = d^.dHypers
        s = MState { -- Citing model
                     _stPsis = let dist n = case d ^. dCitingNodes . at' n . to toList of
                                                []    -> M.empty
                                                nodes -> M.singleton n
                                                         $ symDirMulti (hp^.alphaPsi) nodes
                               in foldMap dist citingNodes
                   , _stGammas = let dist = multinom [ (Shared, hp^.alphaGammaShared)
                                                     , (Own, hp^.alphaGammaOwn) ]
                                 in foldMap (\t->M.singleton t dist) citingNodes
                   , _stOmegas = let dist = symDirMulti (hp^.alphaOmega) (M.keys $ d^.dItems)
                                 in foldMap (\t->M.singleton t dist) citingNodes
                   , _stCiting = M.empty

                   -- Cited model
                   , _stLambdas = let dist = symDirMulti (hp^.alphaLambda) (M.keys $ d^.dItems)
                                      lambdas0 = foldMap (\n->M.singleton n dist) $ M.keys $ d^.dCitedNodes
                                  in foldl' (\dms (n,x)->M.adjust (incMultinom x) (Cited n) dms) lambdas0 (M.elems $ d^.dNodeItems)
                   }

        initCitingUU :: CitingUpdateUnit -> State MState ()
        initCitingUU uu = do
            let err = error $ "CitationInference: Initial value for "++show uu++" not given\n"
                s = maybe err id $ M.lookup (uu^.uuNI) citingInit
            modify' $ setCitingUU uu (Just s)

    in execState (mapM_ initCitingUU $ citingUpdateUnits d) s

modelLikelihood :: MState -> Probability
modelLikelihood model =
    product $ (model ^.. stGammas  . folded . to likelihood)
           ++ (model ^.. stLambdas . folded . to likelihood)
           ++ (model ^.. stOmegas  . folded . to likelihood)
           ++ (model ^.. stPsis    . folded . to likelihood)

instance UpdateUnit CitingUpdateUnit where
    type ModelState CitingUpdateUnit = MState
    type Setting CitingUpdateUnit = CitingSetting
    fetchSetting uu ms = ms ^. stCiting . at' (uu^.uuNI)
    evolveSetting ms uu = categorical $ citingFullCond (setCitingUU uu Nothing ms) uu
    updateSetting uu _ s' = setCitingUU uu (Just s') . setCitingUU uu Nothing

citingProb :: MState -> CitingUpdateUnit -> Setting CitingUpdateUnit -> Double
citingProb st (CitingUpdateUnit {_uuN=n, _uuX=x}) setting =
    let gamma = st ^. stGammas . at' n
        omega = st ^. stOmegas . at' n
        psi = st ^. stPsis . at' n
    in case setting of
        SharedSetting c   -> let lambda = st ^. stLambdas . at' c
                             in sampleProb gamma Shared
                              * sampleProb psi c
                              * sampleProb lambda x
        OwnSetting        ->  sampleProb gamma Own
                            * sampleProb omega x

citingFullCond :: MState -> CitingUpdateUnit -> [(Double, Setting CitingUpdateUnit)]
citingFullCond ms uu = map (\s->(citingProb ms uu s, s)) $ citingDomain ms uu

citingDomain :: MState -> CitingUpdateUnit -> [Setting CitingUpdateUnit]
citingDomain ms uu = do
    s <- [Own, Shared]
    case s of
        Shared -> do c <- uu ^. uuCites . to S.toList
                     return $ SharedSetting c
        Own    -> do return $ OwnSetting

setCitingUU :: CitingUpdateUnit -> Maybe (Setting CitingUpdateUnit) -> MState -> MState
setCitingUU uu@(CitingUpdateUnit {_uuNI=ni, _uuN=n, _uuX=x}) setting ms = execState go ms
  where
    set = maybe Unset (const Set) setting
    go = case maybe (fetchSetting uu ms) id setting of
           SharedSetting c    -> do stPsis .    at' n %= setMultinom set c
                                    stLambdas . at' c %= setMultinom set x
                                    stGammas .  at' n %= setMultinom set Shared
                                    stCiting .  at ni .= setting

           OwnSetting         -> do stOmegas .  at' n %= setMultinom set x
                                    stGammas .  at' n %= setMultinom set Own
                                    stCiting .  at ni .= setting

max_stars_repo_name: bgamari/bayes-stack

max_issues_repo_name: bgamari/bayes-stack

max_forks_repo_name: bgamari/bayes-stack


text:
module Gosadenban
( main3

) where

-------------------------------------------------------------
--module
-------------------------------------------------------------

import Numeric.LinearAlgebra
import System.Random
import qualified Data.ByteString.Lazy as B
import Data.List
import Data.Time

import Tools
import Durstenfeld
import ForMNIST

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

data ActivationFunction = ActivationFunction
    { func  :: Double -> Double ,
      func' :: Double -> Double ,
      description :: String
    }

data Params = Params
    { weight :: Matrix Double ,
      bias :: Matrix Double
    } deriving (Show)

data DParams = DParams
    { dweight :: Matrix Double ,
      dbias :: Matrix Double
    } deriving (Show)

data LayerProp = LayerProp
    { params :: Params ,
      af :: ActivationFunction
    }

data HyperProp = HyperProp
    { vw0 :: Matrix Double ,
      vb0 :: Matrix Double
    }

data Hyperparameter = Hyperparameter
    { learningrate :: Double ,
      alpha0 :: Double ,
      beta0  :: Double ,
      beta1  :: Double
    }

data Network = Network
    { layers :: [LayerProp] ,
      hprop  :: [HyperProp] ,
      hpara  :: Hyperparameter
    }

data PropagateLayer
    = PropagateLayer
        { inData  :: Matrix Double ,
          outData :: Matrix Double ,
          pAf :: ActivationFunction ,
          pParams :: Params
        }
    | PropagateSensorLayer
        { outData :: Matrix Double
        }

data BackPropagateLayer
    = BackPropagateLayer
        { inGrad  :: Matrix Double ,
          outGrad :: Matrix Double ,
          bpAf :: ActivationFunction ,
          bpParams :: Params ,
          dParams :: DParams
        }
    | BackPropagateLastLayer
        { outGrad :: Matrix Double
        }

data ErrorLayer
    = ErrorLayer
        { lastData :: Matrix Double ,
          teachData :: Matrix Double ,
          errorValue :: Double ,
          grad :: Matrix Double
        } deriving (Show)

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

propagate :: PropagateLayer -> LayerProp -> PropagateLayer
propagate dataJ layerK = PropagateLayer
    { inData  = x ,
      outData = y ,
      pAf     = functions ,
      pParams = p
    }
    where x = outData dataJ
          p = params layerK
          w = weight p
          b = bias p
          a = add b $ x Numeric.LinearAlgebra.<> w
          functions =  af layerK
          h = func functions
          y = cmap h a

propagateNetwork :: Matrix Double -> Network -> [PropagateLayer]
propagateNetwork input net = calcs
    where calcs = scanl propagate layer0 (layers net)
          layer0 = PropagateSensorLayer { outData = input}

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

backpropagate :: PropagateLayer -> BackPropagateLayer -> BackPropagateLayer
backpropagate pLayerJ bpLayerK = BackPropagateLayer
    { inGrad = indx ,
      outGrad = dx ,
      bpAf = af ,
      bpParams  = p ,
      dParams = dP
    }
    where af = pAf pLayerJ
          p = pParams pLayerJ
          h' = func' af
          indx = outGrad bpLayerK
          y = outData pLayerJ
          dL' = cmap h' y
          dL = indx * dL'
          db = dL
          x = inData pLayerJ
          w = weight p
          dw = (tr x) Numeric.LinearAlgebra.<> dL
          dx = dL Numeric.LinearAlgebra.<> (tr w)
          dP = DParams { dweight = dw, dbias = db}

backpropagateNetwork :: ErrorLayer -> [PropagateLayer] ->[BackPropagateLayer]
backpropagateNetwork el prL = tail ans1
    where ans1 = scanr backpropagate layerFin prL
          layerFin = BackPropagateLastLayer {outGrad = grad0}
          grad0 = grad el

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

--Softmax & cross entropy error
errorNetwork :: Matrix Double -> Matrix Double -> ErrorLayer
errorNetwork ld td = ErrorLayer
    { lastData = dataFromSoftmaxFunction ,
      teachData = td ,
      errorValue = e / len ,
      grad = g
    }
    where dataFromSoftmaxFunction = softmaxBatch' ld
          g = add dataFromSoftmaxFunction $ scale (-1) td
          y = map asRow $ toRows dataFromSoftmaxFunction
          t = map asRow $ toRows td
          e' = zipWith cee t y
          e = sum e'
          len = length' t

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

inference :: Matrix Double -> Int -> Double
inference x t = ans
    where (r, c) = maxIndex x
          ans = if c == t then 1.0 else 0.0

accuracyNet :: Network -> Matrix Double -> Int -> Double
accuracyNet net x t = a
    where ans' = propagateNetwork x net
          ans  = last $ map outData ans'
          a    = inference ans t

accuracyNetwork :: Network -> [Matrix Double] -> [Int] -> Double
accuracyNetwork net xs ts = sumans / lenans
    where f x t = accuracyNet net x t
          ans = zipWith f xs ts
          sumans = sum ans
          lenans = length' xs

----------------------------------------------------------------------
sgd :: Matrix Double -> Matrix Double -> Hyperparameter -> Matrix Double
sgd x dx hp = x - scale eta dx
    where eta = learningrate hp

update :: Hyperparameter -> Params -> DParams -> Params
update hp p dp = newParam
    where b  = bias p
          w  = weight p
          db = dbias dp
          dw = dweight dp
          new_b = sgd b db hp
          new_w = sgd w dw hp
          newParam = Params {bias = new_b, weight = new_w}

updateNetowrk :: Network -> [DParams] -> Network
updateNetowrk net dp = new_Network
    where layerList = layers net
          hp = hpara net
          hpp = hprop net
          paramList = map params layerList
          afList = map af layerList
          f = update hp
          new_Para = zipWith f paramList dp
          new_Layer = zipWith makeLayerProp new_Para afList
          new_Network = Network{layers = new_Layer, hpara = hp, hprop = hpp  }

momentum :: Matrix Double -> Matrix Double -> Hyperparameter -> Matrix Double -> ( Matrix Double, Matrix Double)
momentum x dx hp v = (x + new_v, new_v)
    where eta = learningrate hp
          alpha = alpha0 hp
          new_v = (scale alpha v) - scale eta dx

adagrad :: Matrix Double -> Matrix Double -> Hyperparameter -> Matrix Double -> ( Matrix Double, Matrix Double)
adagrad x dx hp v = (new_x, new_v')
    where eta = learningrate hp
          alpha = alpha0 hp
          new_v' = v + dx*dx
          new_v = cmap sqrt new_v'
          new_v'' = cmap devdev new_v
          new_x = x - scale eta ( dx/new_v'' )

devdev :: Double -> Double
devdev x =
    if x == 0
        then 1.0
        else x

update' :: Hyperparameter -> (HyperProp , Params) -> DParams -> (Params, HyperProp)
update' hp (hpp, p) dp = (newParam, newHpp)
    where b  = bias p
          w  = weight p
          db = dbias dp
          dw = dweight dp
          vw = vw0 hpp
          vb = vb0 hpp
          (new_w, new_vw) = momentum w dw hp vw
          (new_b, new_vb) = momentum b db hp vb
          newParam = Params {bias = new_b, weight = new_w}
          newHpp = HyperProp {vw0 = new_vw, vb0 = new_vb}

updateNetowrk' :: Network -> [DParams] -> Network
updateNetowrk' net dp = new_Network
    where layerList = layers net
          hp = hpara net
          hpp = hprop net
          paramList = map params layerList
          afList = map af layerList
          f = update' hp
          hppp = zip hpp paramList
          new_Para_And_Hpp = zipWith f hppp dp
          new_Para = map fst new_Para_And_Hpp
          new_Hpp  = map snd new_Para_And_Hpp
          new_Layer = zipWith makeLayerProp new_Para afList
          new_Network = Network{layers = new_Layer, hpara = hp, hprop = new_Hpp  }

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

makeLayerProp :: Params -> ActivationFunction -> LayerProp
makeLayerProp p aF = LayerProp {params = p, af = aF }

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

loadL :: B.ByteString -> [Matrix Double]
loadL xs = tL4
    where tLabel  = loadLabel xs
          f x = pickUpNum x tLabel
          tL1 = map f batch_List
          tL2 = map ( map oneHotEncoder ) tL1
          tL3 = map ( map flatten ) tL2
          tL4 = map fromRows tL3

loadI :: B.ByteString -> [Matrix Double]
loadI xs = tI4
    where tImage  = loadImage xs
          g x = pickUpNum x tImage
          tI1 = map g batch_List
          tI2 = map ( map (scale (1/255.0) )) tI1
          tI3 = map ( map flatten ) tI2
          tI4 = map fromRows tI3

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

finPara :: Matrix Double -> Matrix Double -> Params
finPara b w = Params {bias = b, weight = w}

finLayer :: Params -> ActivationFunction -> LayerProp
finLayer p a = LayerProp {params = p, af = a }

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

sigmoidFunction :: Floating a => a -> a
sigmoidFunction x = 1 / (1 + exp (-x) )

sigmoidFunctionBack :: Floating a => a -> a
sigmoidFunctionBack  x = x*(1-x)

sigmoidFunc :: ActivationFunction
sigmoidFunc = ActivationFunction
    { func  = sigmoidFunction,
      func' = sigmoidFunctionBack,
      description = "Sigmoid Function"
    }

reluFunction :: Double -> Double
reluFunction x
    | x <= 0 = 0.0
    | x > 0  = x

reluFunctionBack :: Double -> Double
reluFunctionBack x
    | x <= 0 = 0.0
    | x > 0  = 1.0

reluFunc :: ActivationFunction
reluFunc = ActivationFunction
    { func  = reluFunction,
      func' = reluFunctionBack,
      description = "Relu Function"
    }

identity :: ActivationFunction
identity = ActivationFunction
    { func  = id,
      func' = backOne,
      description = "id"
    }

backOne :: Double -> Double
backOne x = 1.0

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

softmax :: (Linear t c, Floating t, Container c t) => c t -> c t
softmax m = scale  (1/bunbo) newM
    where maxValue = maxElement m
          newM     = cmap  (\x -> exp $ x - maxValue ) m
          bunbo    = sumElements newM

softmaxBatch' :: Matrix Double -> Matrix Double
softmaxBatch' m  =  m1
    where list  = toRows m
          list1 = fmap softmax list
          m1 = fromRows list1

cee :: Matrix Double -> Matrix Double -> Double
cee t y = (-1.0)*e
    where delta = 0.0000001
          y1 = cmap (\x -> log (x + delta) ) y
          y2 = flatten y1
          y3 = diag y2
          e = sumElements (t Numeric.LinearAlgebra.<> y3)

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

checkIt (e0, net) (x ,t) = (e, newNetwork)
    where ans' = propagateNetwork x net
          ans  = map outData ans'
          answer  = last ans
          answer2 = errorNetwork answer t
          e = errorValue answer2
          finalAnswer' = backpropagateNetwork answer2 ans'
          finalA = map dParams $ init finalAnswer'
          newNetwork = updateNetowrk' net finalA

checking net0 lt = ans
    where ans = scanl checkIt net0 lt

checking' net0 lt = ans
    where ans = foldl checkIt net0 lt

-------------------------------------------------------------
--Parameter
-------------------------------------------------------------

data_size = 10000
batch_size = 100
trials_num = 2000

input_size  = 784
hidden_size = 50
output_size = 10

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

-------------------------------------------------------------
--Initial Condition (?)
-------------------------------------------------------------

nw1 = randomMatrix input_size  hidden_size (mkStdGen 101) :: (Matrix Double)
nw2 = randomMatrix hidden_size output_size (mkStdGen 102) :: (Matrix Double)

nb1 = matrix hidden_size $ take (batch_size * hidden_size) $ repeat 0.0
nb2 = matrix output_size $ take (batch_size * output_size) $ repeat 0.0

nv_w1 =  matrix hidden_size $ take (input_size  * hidden_size) $ repeat 0.0
nv_w2 =  matrix output_size $ take (hidden_size * output_size) $ repeat 0.0

nv_b1 =  matrix hidden_size $ take (batch_size * hidden_size) $ repeat 0.0
nv_b2 =  matrix output_size $ take (batch_size * output_size) $ repeat 0.0

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

param1 = Params {bias = nb1, weight = nw1}
layer1 = LayerProp {params = param1, af = reluFunc}
hprop1 = HyperProp { vw0 = nv_w1, vb0 = nv_b1 }

param2 = Params {bias = nb2, weight = nw2}
layer2 = LayerProp {params = param2, af = identity }
hprop2  = HyperProp { vw0 = nv_w2, vb0 = nv_b2}

hyperParam = Hyperparameter {learningrate = 0.0001, alpha0 = 0.9, beta0 = 0.0, beta1 = 0.0}

networkA = Network {layers = [layer1,layer2], hpara = hyperParam, hprop = [hprop1,hprop2]}

--makeBatchList gen num_of_trials data_size batch_size
batch_List = makeBatchList 100000 trials_num data_size batch_size

main3 = do

    z <- getZonedTime
    print z

    x <- getCurrentTime

    contents1 <- B.readFile "MNIST/train-labels.idx1-ubyte"
    contents2 <- B.readFile "MNIST/train-images.idx3-ubyte"

    let tL4 = loadL contents1
        --tL4 = [[Batch X 10],[Batch X 10],...] :: [Matrix Double]

        tI4 = loadI contents2
        --tI4 = [[Batch X data_size],[Batch X data_size],...] :: [Matrix Double]

        test_Image_and_Label = zip tI4 tL4

        --e1 = checking' (0.0,networkA) test_Image_and_Label
        --e2 =  fst e1

        e1' =  checking (0.0,networkA) test_Image_and_Label
        e2' =  tail $ map fst e1'
        net1' = last $ map snd e1'

    --print e2
    writeFile "logs/e_mmt-0.0001.dat" $ encorderForDat2 e2'

    contents3 <- B.readFile "MNIST/t10k-labels.idx1-ubyte"
    contents4 <- B.readFile "MNIST/t10k-images.idx3-ubyte"

    let xs = loadImageTest contents4
        ts = loadLabelTest contents3
        --a = accuracyNet net1' xs ts
        --acc = accuracyNetwork net1' xs ts
        bs = map toRows $ map bias $ map params $ layers net1'
        bs' = map (foldl1 add) bs
        bsl = map recip $ map length' bs
        newbs = map asRow $ zipWith scale bsl bs'
        w = map weight $ map params $ layers net1'
        para = zipWith finPara newbs w
        acfunc = map af $ layers net1'
        newLayer = zipWith finLayer para acfunc
        finNetwork = Network {layers = newLayer, hpara = Hyperparameter {learningrate = 1.0, alpha0 = 0.0, beta0 = 0.0, beta1 = 0.0},hprop = [hprop1,hprop2]}

        acc = accuracyNetwork finNetwork xs ts

    print acc

    y <- getCurrentTime
    print$ diffUTCTime y x

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

loopCheck x net t = (e, finalA, newNetwork)
    where ans' = propagateNetwork x net
          ans  = map outData ans'
          answer  = last ans
          answer2 = errorNetwork answer t
          e = errorValue answer2
          finalAnswer' = backpropagateNetwork answer2 ans'
          finalA = map dParams $ init finalAnswer'
          newNetwork = updateNetowrk net finalA

loop _  _  _ 0 = []
loop x network t n = e : loop x newNetwork t (n-1)
    where  (e, finalA, newNetwork) = loopCheck x network t

randomList = randomChoice 10 10 (mkStdGen 10 )

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

max_stars_repo_name: llbxg/Fukami

max_issues_repo_name: llbxg/Fukami

max_forks_repo_name: llbxg/Fukami


text:
module OidTableSVG where

import Data.Word
import qualified Data.ByteString.Lazy as LB
import Control.Monad hiding (forM_)
import Data.Foldable (forM_)
import Text.Printf
import Control.Arrow ((***))
import Data.String
import Data.Complex

import qualified Text.Blaze.Svg11 as S
import qualified Text.Blaze.Svg11.Attributes as A
import Text.Blaze.Svg11 ((!))
import Text.Blaze.Svg.Renderer.Utf8

import OidCode
import KnownCodes
import Types
import Utils

type Point = Complex Double

oidTableSvg :: Conf -> Bool -> String -> [(String, Word16)] -> LB.ByteString
oidTableSvg conf usePNG title entries
    | entriesPerPage < 1 = error "OID codes too large to fit on a single page"
    | otherwise = renderSvg $
        S.docTypeSvg  ! A.version (S.toValue "1.1")
                      ! A.width  (S.toValue (printf "%fmm" (a4w/mm) :: String))
                      ! A.height (S.toValue (printf "%fmm" (a4h/mm) :: String))
                      ! A.viewbox (S.toValue (printf "0 0 %f %f" a4w a4h :: String))
                      ! A.fontFamily (S.toValue "sans-serif")
                      $ do
    let patid d c | d == show c = d
                  | otherwise   = printf "%s-%d" d c

    -- Create patterns for the codes
    S.defs $ forM_ entries $ \(d,c) ->
        case code2RawCode c of
            Nothing -> return ()
            Just rc -> oidSVGPattern conf usePNG (patid d c) rc

    -- For SVG, we put all on one page (and exceed the page if it is too big)
    let chunks = [entries]
    let totalPages = length chunks

    forM_ (zip [1::Int ..] chunks) $ \(pageNum, thisPage) -> do
            S.text_ ! A.x (S.toValue (a4w / 2))
                    ! A.y (S.toValue (padTop + titleHeight))
                    ! A.textAnchor (S.toValue "middle")
                    ! A.stroke (S.toValue "black")
                    ! A.fontSize (S.toValue (printf "%f" (12*pt) :: String))
                    $ fromString title

            S.text_ ! A.x (S.toValue padLeft)
                    ! A.y (S.toValue (a4h - padBot))
                    ! A.textAnchor (S.toValue "left")
                    ! A.stroke (S.toValue "black")
                    ! A.fontSize (S.toValue (printf "%f" (8*pt) :: String))
                    $ fromString $ "Created by tttool-" ++ tttoolVersion

            forM_ (zip thisPage positions) $ \((d,c), x :+ y) -> do
                S.rect ! A.width (S.toValue imageWidth)
                       ! A.height (S.toValue imageHeight)
                       ! A.x (S.toValue x)
                       ! A.y (S.toValue y)
                       ! A.fill (S.toValue $ "url(#"++patid d c++")")

                S.text_ ! A.x (S.toValue x)
                        ! A.y (S.toValue (y + imageHeight + subtitleSep + subtitleHeight))
                        ! A.textAnchor (S.toValue "left")
                        ! A.stroke (S.toValue "black")
                        ! A.fontSize (S.toValue (printf "%f" (8*pt) :: String))
                        $ fromString d
  where
    -- Configure-dependent dimensions (all in pt)
    (imageWidth,imageHeight) = (*mm) *** (*mm) $ fromIntegral *** fromIntegral $cCodeDim conf

    -- Static dimensions (all in pt)

    -- Page paddings
    padTop, padLeft, padBot, padRight :: Double
    padTop   = 1*cm
    padBot   = 1*cm
    padLeft  = 2*cm
    padRight = 2*cm

    titleHeight  = 1*cm
    titleSep     = 0.5*cm
    footerHeight = 0.5*cm
    footerSep    = 0.5*cm

    imageSepH = 0.4*cm
    imageSepV = 0.2*cm

    subtitleHeight = 0.4*cm
    subtitleSep    = 0.2*cm

    -- Derived dimensions (all in pt)
    {-
    titleRect = Rectangle
        (padLeft          :+ (a4h - padTop - titleHeight))
        ((a4w - padRight) :+ (a4h - padTop))
    titleFont = Font (PDFFont Helvetica 12) black black

    footerRect = Rectangle
        (padLeft          :+ padBot)
        ((a4w - padRight) :+ (padBot + footerHeight))
    footerFont = Font (PDFFont Helvetica 8) black black

    bodyFont = Font (PDFFont Helvetica 8) black black
    -}

    bodyWidth  = a4w - padLeft - padRight
    bodyHeight = a4h - padTop - titleHeight - titleSep - footerSep - footerHeight - padBot

    positions = map (+(padLeft :+ (padTop + titleHeight + titleSep))) $
        calcPositions bodyWidth  bodyHeight
                      imageWidth (imageHeight + subtitleSep + subtitleHeight)
                      imageSepH  imageSepV
    entriesPerPage = length positions

    -- Derived dimensions (all in pixels)
    imageWidthPx = floor (imageWidth * px)
    imageHeightPx = floor (imageHeight * px)

    -- config-dependent conversion factors
    px :: Double
    px = fromIntegral (cDPI conf) / 72

    {-
    -- Makes sure the given point is at a coordinate that is a multiple
    -- of an pixel
    align :: Point -> Point
    align pos = alignToPx (realPart pos) :+ (a4h - alignToPx (a4h - imagPart pos))

    -- Makes sure the given distance is an interal mulitple of a pixel
    alignToPx :: Double -> Double
    alignToPx x = fromIntegral (floor (x * px)) / px
    -}

calcPositions
    :: Double -- ^ total width
    -> Double -- ^ total height
    -> Double -- ^ entry width
    -> Double -- ^ entry height
    -> Double -- ^ pad width
    -> Double -- ^ pad height
    -> [Point]
calcPositions tw th ew eh pw ph = [ x :+ ({-(th - -} y) | y <- ys , x <- xs]
  where
    xs = [0,ew+pw..tw-ew]
    ys = [0,eh+ph..th-eh]

pt :: Double
pt = 48/2.83465

-- Conversation factor
cm :: Double
cm = 10 * mm

mm :: Double
mm = 2.83465 * pt

-- A4 dimensions
a4w, a4h :: Double
a4w = 595 * pt
a4h = 842 * pt

max_stars_repo_name: dalijolijo/tip-toi-reveng

max_issues_repo_name: dalijolijo/tip-toi-reveng

max_forks_repo_name: dalijolijo/tip-toi-reveng


text:
{-# LANGUAGE FlexibleContexts     #-}
{-# LANGUAGE InstanceSigs         #-}
{-# LANGUAGE TypeFamilies         #-}
{-# LANGUAGE UndecidableInstances #-}

module AI.Network.SOM
( SOM(..)
, MapDefinition(..)

, randomNeuron
, makeVectors
, reshapeList
, distance
) where

import           AI.Layer
import           AI.Network
import           AI.Neuron

import           Numeric.LinearAlgebra
import           System.Random

-- | The SOM definition is simple, it only contains a 2-dimensional list of weights
data SOM = SOM { neuronMap :: [[Vector Double]] }

-- | A definitution type for the SOM, it contains the dimensions of each layer (x, y)
--   and the dimension of the input vector (dim)
data MapDefinition = MapDefinition { x        :: Int
                                   , y        :: Int
                                   , inputDim :: Int
                                   }

instance Network SOM where
  type Parameters SOM g = MapDefinition

  predict :: Vector Double -> SOM -> Vector Double
  predict inputs network = inputs

  -- | Create a SOM and initialize it with given weights
  createNetwork :: (RandomGen g) => RandomTransform -> g -> Parameters SOM g -> SOM
  createNetwork transformation g def = SOM randomVectors
    where randomVectors = reshapeList (x def) $
                          makeVectors transformation g (inputDim def)
                          (x def * y def)

-- | A helper function to reshape a 1D list into a 2D list
reshapeList :: Int -> [a] -> [[a]]
reshapeList x [] = [[]]
reshapeList x lst = h : reshapeList x t
  where (h, t) = splitAt x lst

-- | Create a random set of weights for a given neuron
randomNeuron :: (RandomGen g) => RandomTransform -> g -> Int -> Vector Double
randomNeuron transform g inputDim = inputDim |> randomList transform g

-- | Make a 1D list of vectors to be used by the SOM in creating a map of weights
makeVectors :: (RandomGen g) => RandomTransform -> g -> Int -> Int -> [Vector Double]
makeVectors transform g inputDim 0 = []
makeVectors transform g inputDim num = randomNeuron transform g' inputDim :
                                       makeVectors transform g'' inputDim (num - 1)
  where (g', g'') = split g

-- | Calculate the distance between a SOM neuron and an input
distance :: Vector Double -> Vector Double -> Double
distance a b = sqrt $ sum $ map (^2) $ zipWith (-) (toList a) (toList b)

max_stars_repo_name: jbarrow/LambdaNet

max_issues_repo_name: jbarrow/LambdaNet

max_forks_repo_name: jbarrow/LambdaNet


text:
{-# LANGUAGE ConstraintKinds #-}

module Projection
( projectL2
)
where

import           Basis
import           Element
import           Node
import qualified Numeric.LinearAlgebra as LinAl
import           Quadrature
import           ShapeFcns

-- Synonym for numeric/fractional constraints
type FrElNuFi a = (Fractional a,LinAl.Element a,LinAl.Numeric a,LinAl.Field a)

-- Element mass integrand
elemMassIntegrand :: (Element e,ShapeFcn s,Basis b,FrElNuFi a) => e a -> s b -> [a] -> LinAl.Matrix a
elemMassIntegrand elem shpFcn xi = LinAl.fromLists [[
  ShapeFcns.n shpFcn xi i [0] * ShapeFcns.n shpFcn xi j [0] * Element.computeJacobianDet elem shpFcn xi
  | i <- idxRange] | j <- idxRange]
  where
    idxRange = [0..(ShapeFcns.getShapeFcnOrder shpFcn)]

-- This function computes the element load integrand when given a function
-- to integrate over an element.
elemLoadIntegrand :: (Element e,ShapeFcn s,Basis b,FrElNuFi a) => e a -> s b -> ([a] -> a) -> [a] -> LinAl.Matrix a
elemLoadIntegrand elem shpFcn f xi = LinAl.scale (computeJacobianDet elem shpFcn xi) (LinAl.fromLists $ map (:[]) [f physCoords * n shpFcn xi i [0] | i <- [0..(getShapeFcnOrder shpFcn)]])
  where
    physCoords = [sum $ zipWith (*) (map (head . nodeCoordinates) (getElementNodes elem)) (map (\j -> n shpFcn xi j [0]) [0..(getShapeFcnOrder shpFcn)])]

-- Function for projecting a function onto a basis over a single element.
projectL2 :: (Element e,ShapeFcn s,Basis b,FrElNuFi a) => e a -> s b -> Int -> ([a] -> a) -> Maybe (LinAl.Matrix a)
projectL2 elem shpFcn ngpts f = LinAl.linearSolve (integrate 1 ngpts (elemMassIntegrand elem shpFcn)) (integrate 1 ngpts (elemLoadIntegrand elem shpFcn f))


text:
{-# LANGUAGE DeriveDataTypeable #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE QuasiQuotes #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE UndecidableInstances #-}

module Main where

import Prelude hiding (map)
import qualified Prelude as P

import Control.Monad
import Criterion.Config
import Criterion.Main
import Criterion.Monad
import Criterion.Plot
import Data.Array.Vector
import Nikola
import Statistics.Function (minMax)
import Statistics.KernelDensity (epanechnikovPDF)
import Statistics.Types
import System.Environment (getArgs)

main :: IO ()
main = do
    (cfg, _) <- parseArgs defaultConfig defaultOptions =<< System.Environment.getArgs
    withNewContext $ \_ -> do
    reify f >>= compileTopFun "f" >>= print
    samples <- withCompiledFunction f $ \f -> do
               replicateM 10000 $ do
               timeKernel f $ \f -> do
               call f [1..1024]
    withConfig cfg $
        plotAll [("map", (toU . P.map (realToFrac . fst)) samples)]
  where
    f :: Exp [Int] -> Exp [Int]
    f = map (\x -> x + 1)

plotAll :: [(String, Sample)] -> Criterion ()
plotAll descTimes = forM_ descTimes $ \(desc,times) -> do
  plotWith Timing $ \o -> plotTiming o desc times
  plotWith KernelDensity $ \o -> uncurry (plotKDE o desc extremes)
                                     (epanechnikovPDF 100 times)
  where
    extremes = case descTimes of
                 (_:_:_) -> toJust . minMax . concatU . P.map snd $ descTimes
                 _       -> Nothing
    toJust r@(lo :*: hi)
        | lo == infinity || hi == -infinity = Nothing
        | otherwise                         = Just r
        where infinity                      = 1/0

max_stars_repo_name: mainland/nikola

max_issues_repo_name: mainland/nikola

max_forks_repo_name: mainland/nikola


text:
module SimilaritySpec where

import Test.Hspec
import Test.QuickCheck
import qualified Numeric.LinearAlgebra as NL

import Data.List

import Eris.Compute.Similarity
import ErisTestUtility

spec :: Spec
spec = do
    describe "L1-norm based. Property tests:" $ do
      it "SAD is L1 norm of the difference of two vectors:" $ property $
        \v1 v2 -> (length v1 /= length v2) || sumAbsoluteDifference v1 v2 == sum ( abs <$> zipWith (-) v1 v2)
      it "SAD is symmetric"$ property $ -- update this property
        \l1 l2 -> (length l1 /= length l2)  || (sumAbsoluteDifference l2 l1 == sum (fmap abs (zipWith (-) l1 l2)))
      it "manhattanDistance is SAD" $ property $
        \v1 v2 -> (length v1 /= length v2) || manhattanDistance v1 v2 == sumAbsoluteDifference v1 v2
      it "taxicab distance is SAD" $ property $
        \v1 v2 -> (length v1 /= length v2) || taxicab v1 v2 == sumAbsoluteDifference v1 v2
      it "mean absolute distance is SAD average on number of dimension" $ property $
        \v1 v2 -> (length v1 /= length v2) || meanAbsoluteDifference v1 v2 == average (abs <$> zipWith (-) v1 v2 :: [Double])
    describe "mean absolute Difference " $do
      it "naive test: based on Sum Absolute Difference" $
        meanAbsoluteDifference tl1 tl2 `shouldBe` sumAbsoluteDifference tl1 tl2 / (fromIntegral . length $ tl1)
      it "Pro test: symmetric" $ property $
        \l1 l2 -> (length l1 /= length l2)  || (meanAbsoluteDifference l2 l1 == sum (fmap abs (zipWith (-) l1 l2)) / fromIntegral (length l2))

    describe "based on l2-norm" $ do
      it "euclidean distance is l2-norm" $
        euclideanDistance tl1 tl2 `shouldBe` NL.norm_2 (NL.vector tl1 - NL.vector tl2)

    describe "mse distance: " $ do
      it "It is 0 with it self. " $
        (fromInteger. round) (meanSquaredDistance [1.0,1.0,1.0] [1.0,1.0,1.0]) `shouldBe` (0:: Double)
      it "It is squared L2-norm normal by number of dimension." $ do
        meanSquaredDistance [4.0, 0.0] [0.0, 5.0] `shouldBe` (4^2+5^2)/2
        meanSquaredDistance [1.0,1.0] [1.0,-1.0] `shouldBe` (0+4)/2

    describe "Cosine Similarity: " $ do
      it "The vector and itself is one." $
        (fromInteger. round) (cosineSimilarity [1.0,1.0,1.0] [1.0,1.0,1.0]) `shouldBe` (1.0:: Double)
      it "The similarity of two linear independent vetors is not 1 " $
        diff5Dec (cosineSimilarity [1,2] [1,3]) 1.0 > 0 `shouldBe` True
      it "The similarity of two linear dependent vetors is 1 " $
        diff5Dec (cosineSimilarity [1,2] [2,4]) 1.0 `shouldBe` 0
      it " It is 1 for linear dependencies: " $ do
        diff5Dec (cosineSimilarity [-1,2,-3,4,-5] [-4,5,-10,11,-16]) 1.0 `shouldBe` 0

      it "The vector and any perpendicular vector is 0." $ do
        cosineSimilarity [1.0, 0.0] [0.0, 1.0] `shouldBe` 0
        cosineSimilarity [1.0,1.0] [1.0,-1.0] `shouldBe` 0

    describe "pearson correlation coefficient: " $ do
      it " It is 1 for linear dependencies: " $ do
          (pearsonCC' [-1,2,-3,400,-30] [-4,5,-10,1199,-91]) `shouldBe` (1.0:: Double)

    describe "Pearson CC is a special Case of consine similarity: " $ do
        let l1 = [-3..1]
            l2 = fmap (\x -> x*3+120) l1
        it " The value of PCC and CosineSi: " $
          (pearsonCC l1 l2) `shouldBe` (cosineSimilarity l1 l2)
        it "pearson CC base on z score: " $
          (pearsonCC l1 l2) `shouldBe` (pearsonCC' l1 l2)
    where
        tl1 = [1,-1,2,-2]
        tl2 = [-1,1,-2,2]
        diff5Dec = relativeEq 5

max_stars_repo_name: emmettng/eris

max_issues_repo_name: emmettng/eris

max_forks_repo_name: emmettng/eris


text:
{-# LANGUAGE DeriveDataTypeable, DeriveGeneric, OverloadedStrings #-}
-- |
-- Module    : Statistics.Distribution.DiscreteUniform
-- Copyright : (c) 2016 André Szabolcs Szelp
-- License   : BSD3
--
-- Maintainer  : a.sz.szelp@gmail.com
-- Stability   : experimental
-- Portability : portable
--
-- The discrete uniform distribution. There are two parametrizations of
-- this distribution. First is the probability distribution on an
-- inclusive interval {1, ..., n}. This is parametrized with n only,
-- where p_1, ..., p_n = 1/n. ('discreteUniform').
--
-- The second parametrization is the uniform distribution on {a, ..., b} with
-- probabilities p_a, ..., p_b = 1/(a-b+1). This is parametrized with
-- /a/ and /b/. ('discreteUniformAB')

module Statistics.Distribution.DiscreteUniform
    (
      DiscreteUniform
    -- * Constructors
    , discreteUniform
    , discreteUniformAB
    -- * Accessors
    , rangeFrom
    , rangeTo
    ) where

import Control.Applicative (empty)
import Data.Aeson   (FromJSON(..), ToJSON, Value(..), (.:))
import Data.Binary  (Binary(..))
import Data.Data    (Data, Typeable)
import System.Random.Stateful (uniformRM)
import GHC.Generics (Generic)

import qualified Statistics.Distribution as D
import Statistics.Internal

-- | The discrete uniform distribution.
data DiscreteUniform = U {
      rangeFrom  :: {-# UNPACK #-} !Int
    -- ^ /a/, the lower bound of the support {a, ..., b}
    , rangeTo    :: {-# UNPACK #-} !Int
    -- ^ /b/, the upper bound of the support {a, ..., b}
    } deriving (Eq, Typeable, Data, Generic)

instance Show DiscreteUniform where
  showsPrec i (U a b) = defaultShow2 "discreteUniformAB" a b i
instance Read DiscreteUniform where
  readPrec = defaultReadPrecM2 "discreteUniformAB" (\a b -> Just (discreteUniformAB a b))

instance ToJSON   DiscreteUniform
instance FromJSON DiscreteUniform where
  parseJSON (Object v) = do
    a <- v .: "uniformA"
    b <- v .: "uniformB"
    return $ discreteUniformAB a b
  parseJSON _ = empty

instance Binary DiscreteUniform where
  put (U a b) = put a >> put b
  get         = discreteUniformAB <$> get <*> get

instance D.Distribution DiscreteUniform where
  cumulative (U a b) x
    | x < fromIntegral a = 0
    | x > fromIntegral b = 1
    | otherwise = fromIntegral (floor x - a + 1) / fromIntegral (b - a + 1)

instance D.DiscreteDistr DiscreteUniform where
  probability (U a b) k
    | k >= a && k <= b = 1 / fromIntegral (b - a + 1)
    | otherwise        = 0

instance D.Mean DiscreteUniform where
  mean (U a b) = fromIntegral (a+b)/2

instance D.Variance DiscreteUniform where
  variance (U a b) = (fromIntegral (b - a + 1)^(2::Int) - 1) / 12

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

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

instance D.Entropy DiscreteUniform where
  entropy (U a b) = log $ fromIntegral $ b - a + 1

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

instance D.ContGen DiscreteUniform where
  genContVar d = fmap fromIntegral . D.genDiscreteVar d

instance D.DiscreteGen DiscreteUniform where
  genDiscreteVar (U a b) = uniformRM (a,b)

-- | Construct discrete uniform distribution on support {1, ..., n}.
--   Range /n/ must be >0.
discreteUniform :: Int             -- ^ Range
                -> DiscreteUniform
discreteUniform n
  | n < 1     = error $ msg ++ "range must be > 0. Got " ++ show n
  | otherwise = U 1 n
  where msg = "Statistics.Distribution.DiscreteUniform.discreteUniform: "

-- | Construct discrete uniform distribution on support {a, ..., b}.
discreteUniformAB :: Int             -- ^ Lower boundary (inclusive)
                  -> Int             -- ^ Upper boundary (inclusive)
                  -> DiscreteUniform
discreteUniformAB a b
  | b < a     = U b a
  | otherwise = U a b

max_stars_repo_name: vaerksted/statistics

max_issues_repo_name: vaerksted/statistics

max_forks_repo_name: vaerksted/statistics


text:
module Ch03.Matrix where

import Prelude hiding (and, or)
import Numeric.LinearAlgebra

import Control.Exception

-- Matrix
as :: Matrix R
as = matrix 4 [1,2,3,4]
--  (1><4)
-- [ 1.0, 2.0, 3.0, 4.0 ]

as' :: Vector R
as' = vector [1,2,3,4]
-- [1.0,2.0,3.0,4.0]

bs :: Matrix R
bs = (3><2) [1..6]
-- (3><2)
--  [ 1.0, 2.0
--  , 3.0, 4.0
--  , 5.0, 6.0 ]

-- Dimensions
ndim  :: Matrix t -> Int
ndim     = rows

shape :: Matrix t -> (Int, Int)
shape xs = (rows xs, cols xs)

t1 = assert (shape as == (1,4)) "OK"
t2 = assert (shape bs == (3,2)) "OK"

-- Inner Product
as1 :: Matrix R
as1 = (2><2) [1..4]
bs1 :: Matrix R
bs1 = (2><2) [5..8]

as1bs1 :: Matrix R
as1bs1 = as1 <> bs1
-- (2><2)
--  [ 19.0, 22.0
--  , 43.0, 50.0 ]

as2 :: Matrix R
as2 = (2><3) [1..6]
bs2 :: Matrix R
bs2 = (3><2) [1..6]

as2bs2 :: Matrix R
as2bs2 = as2 <> bs2
-- (2><2)
--  [ 22.0, 28.0
--  , 49.0, 64.0 ]

c :: Matrix R
c = (2><2) [1..4]

-- Inconsistent Dimensions
as2c :: Matrix R
as2c = as2 <> c
-- *** Exception: inconsistent dimensions in matrix product (2,3) x (2,2)
-- CallStack (from HasCallStack):
--   error, called at src/Internal/LAPACK.hs:56:31 in hmatrix-0.17.0.2-HoBfYSJ1dMG9PoVkmbcnW2:Internal.LAPACK


text:
module Budget.Budgeter(
  budgeter
  )where

import Numeric.LinearAlgebra
import Budget.Types
import Budget.Parse
import Budget.Printer
import Text.Parsec
import Finance.FutureValue

getNeededExpense :: BudgetItem -> Expense
getNeededExpense bs = case needed bs of
  True -> Expense (name bs) (minPrice bs)
  False -> Expense (name bs) 0.0

reduceBs :: [BudgetItem] -> [Expense] -> [BudgetItem]
reduceBs bs xs = map removeMin $ zip bs xs

removeMin :: (BudgetItem, Expense) -> BudgetItem
removeMin ((BudgetItem s i l h n), Expense _ m) = BudgetItem s i (l-m) (h-m) n

aDiag :: BudgetItem -> Double
aDiag (BudgetItem _ i low high _)
  | low == high = 100000
  | otherwise = (-2.0) * i * k**2
  where k = 1.0 / (high-low)

aDiags :: [BudgetItem] -> [Double]
aDiags = map aDiag

type CurrentIndex = Integer
type Rank = Integer
aList :: CurrentIndex -> Rank -> [Double] -> [Double]
aList c r (x:xs)
  |c `mod` (r+1) == 0 = x : aList (c+1) r xs
aList c r xs
  | c == r^2-1 = [0]
  | ((c+1) `mod` r == 0) && (c /= r^2-1) = -1 : aList (c+1) r xs
  | (c `mod` (r+1) /= 0 && c >= r^2-r) = 1 : aList (c+1) r xs
  | otherwise = 0 : aList (c+1) r xs

availableCosts :: Double -> Double -> Double -> Double -> Double -> Double
availableCosts cpm y v ipm apr
  | cpm <= 0 = 0
  | retirement > sustainableRetirement = cpm
  | otherwise = availableCosts (cpm-10) y v ipm apr
  where retirement = futureValue numberMonths v (ipm-cpm) ratePerMonth
        sustainableRetirement = costsPerYear / apr
        numberMonths = floor (12 * y)
        ratePerMonth = apr/12
        costsPerYear = cpm * 12

aMatrix :: [BudgetItem] -> Matrix Double
aMatrix bs = let as = aDiags bs
                 r = (length as) + 1
             in (r><r) $ aList 0 (toInteger r) as

bElem :: BudgetItem -> Double
bElem (BudgetItem _ i low high _)
  | low == high = 100000
  | otherwise = (-2.0) * i * k**2 * high
  where k = 1.0 / (high-low)

bElems :: [BudgetItem] -> [Double]
bElems = map bElem

type TotalCost = Double
bVector :: [BudgetItem] -> TotalCost -> Vector Double
bVector bs tc = let vecElems = bElems bs ++ [tc]
                in vector vecElems

removeNeeded :: BudgetInfo -> (BudgetInfo, [Expense])
removeNeeded (BudgetInfo i s y a bs) = (BudgetInfo (i - totalNeeded) s y a newBs, needs)
  where needs = map getNeededExpense bs
        newBs = reduceBs bs needs
        totalNeeded = sum [x | (Expense _ x) <- needs]

type ItemNames = [String]
calculateBudget :: BudgetInfo -> [Expense]
calculateBudget (BudgetInfo i s y a bs)
  | i > totalMaxCosts = expenses
  | otherwise = calculateBudget' [] (BudgetInfo i s y a bs)
  where maxCosts = map maxPrice bs
        names = map name bs
        expenses = [Expense n m | (n, m) <- zip names maxCosts]
        totalMaxCosts = sum maxCosts

inefficientNames :: [BudgetItem] -> Costs -> ItemNames
inefficientNames bs costs  = [(name b)| (b,c) <- (zip bs costs), c < 0]

optimizedCosts :: Double -> [BudgetItem] -> [Double]
optimizedCosts i bs = init $ toList $ (aMatrix bs) <\> (bVector bs (realToFrac i))

calculateBudget' :: ItemNames -> BudgetInfo -> [Expense]
calculateBudget' unusedItems (BudgetInfo i s y a bs)
  | validCosts = expenses
  | otherwise = calculateBudget' newUnusedNames cleanedBI
  where costs = optimizedCosts i bs
        validCosts = length [c | c <- costs, c < 0] == 0
        allNames = (map name bs) ++ unusedItems
        newUnusedNames = unusedItems ++ inefficientNames bs costs
        cleanedBI = BudgetInfo i s y a [b | (b,c) <- (zip bs costs), c > 0]
        appendedCosts = costs ++ (replicate ((length allNames) - (length costs)) 0)
        expenses = [Expense n m | (n, m) <- zip allNames appendedCosts]

storeSavings (BudgetInfo ipm v ytr apr bs) = BudgetInfo cpm v ytr apr bs
  where cpm = availableCosts ipm ytr v ipm apr

unsafe :: Either ParseError BudgetInfo -> BudgetInfo
unsafe (Right x) = x

yearlySavings (BudgetInfo ipm _ _ _ _) cs = ipm - (sum cs)

combineExpenses :: [Expense] -> [Expense] -> [Expense]
combineExpenses e1s e2s = [Expense n1 (cost1 + cost2) | (Expense n1 cost1) <- e1s, (Expense n2 cost2) <- e2s,
                           n1 == n2]

getBudget :: BudgetInfo -> String
getBudget bi = let saved_bi = storeSavings bi
                   (available_bi, mandatory_costs) = removeNeeded saved_bi
                   distro = calculateBudget available_bi
                   finalExpenses = combineExpenses mandatory_costs distro
               in fullPrint bi finalExpenses

budgeter :: String -> String
budgeter s = let bi = parse getBudgetInfo "Fetching budget info" s
             in case fmap getBudget bi of
                Left parsecError -> "Encountered Input Parsing errror:\n" ++ show parsecError
                Right finalString -> finalString

max_stars_repo_name: PotentialParadox/budgeter

max_issues_repo_name: PotentialParadox/budgeter

max_forks_repo_name: PotentialParadox/budgeter


text:
{-# LANGUAGE FlexibleInstances #-}

import Debug.Trace
-- import Control.Monad.State.Lazy

import System.Environment (getArgs)
import System.IO (readFile)
import Data.Map.Strict (Map, (!), insert, elems, fromList, toList, findWithDefault, size, empty, member, findMin, findMax, singleton, filter)
import qualified Data.Map.Strict as M
import Data.Set (Set)
import qualified Data.Set as S
--import qualified Data.Array as A
import Data.List (find, intercalate, intersperse, permutations, inits, tails, isPrefixOf)
import Data.List.Split (splitOn)
-- import Data.Complex (Complex((:+)), realPart, imagPart) -- define my own complex

-- import UI.NCurses
-- import Data.Char (chr, ord)
-- -- import Data.Complex (Complex((:+)))
import Data.Maybe (fromJust, fromMaybe, isJust)
type Instructions = Map Integer Integer

data ComputerState = Ready | Running | Blocked | Done
  deriving (Show, Eq)

data Computer = Computer {
  state :: ComputerState,
  memory :: Map Integer Integer,
  iptr :: Integer, -- instruction pointer
  base :: Integer, -- base offset
  input :: [Integer],
  output :: [Integer],
  time :: Integer}

computer0 = Computer { state = Ready, memory = empty, iptr = 0, base = 0, input = [], output = [], time = 0 }

run1 :: Computer -> Computer -- (Integer, Integer) -> Instructions -> [Integer] -> [Integer]
run1 c
  | state c == Ready = run1 $ c {state = Running}
  | state c == Done = c
  | state c == Blocked = if null (input c) then c else run1 $ c {state = Running}
  | otherwise = -- running
    --traceShow (iptr c, base c, memory c) $
    case instr `mod` 100 of
      1 -> -- add
        addTime $ c {iptr = i+4, memory = insert (addr 3) (arg 1 + arg 2) instructions}
      2 -> -- multiply
        addTime $ c {iptr = i+4, memory = insert (addr 3) (arg 1 * arg 2) instructions}
      3 -> -- read input
        if (null $ input c)
        then c {state = Blocked}
        else addTime $ c {iptr = i+2, memory = insert (addr 1) (head $ input c) instructions, input = tail $ input c}
      4 -> -- output
        --traceShow (arg 1) $
        addTime $ c {iptr = i+2, output = output c ++ [arg 1]}
      5 -> -- jump-if-true
        addTime $ c {iptr = if arg 1 == 0 then i+3 else arg 2}
      6 -> -- jump-if-false
        addTime $ c {iptr = if arg 1 == 0 then arg 2 else i+3}
      7 -> -- less than
        addTime $ c {iptr = i+4, memory = insert (addr 3) (if arg 1 < arg 2 then 1 else 0) instructions}
      8 -> -- equals
        addTime $ c {iptr = i+4, memory = insert (addr 3) (if arg 1 == arg 2 then 1 else 0) instructions}
      9 -> -- set relative base
        addTime $ c {iptr = i+2, base = base c + arg 1}
      99 -> -- halt
        addTime $ c {state = Done}
      _ -> error "unknown opcode"
  where instructions = memory c
        i = iptr c
        instr = instructions!i
        ii x = findWithDefault 0 x instructions
        arg :: Integer -> Integer
        arg n = case (instr `mod` (100*10^n)) `div` (10*10^n) of
          0 -> ii $ ii $ i+n
          1 -> ii $ i+n
          2 -> ii ((ii $ i+n) + base c)
          _ -> error "bad argument mode"
        addr n = case (instr `mod` (100*10^n)) `div` (10*10^n) of
          0 -> ii $ i+n
          1 -> error "address in mode 1"
          2 -> --trace "address in mode 2" $
            ii (i+n) + base c
          _ -> error $ "address in unknown mode"
        addTime c = c{time = time c + 1}

run :: Computer -> Computer
run c | state c == Done = c
      | state c == Blocked = if null (input c) then c else run $ run1 c
      | length (output c) == 3 = c
      | otherwise = run $ run1 c

type NAT = (Maybe (Integer, Integer), Map Int Int, Maybe Integer)
  -- (packet, blocked computers, last y sent to computer0)

chainif :: a -> [(Bool, a->a)] -> a
-- chainif runs a value through a series of conditional transforms. If
-- a meets a test, it is transformed and continues in the testing
-- chain.  This makes the code look clean, by not having to explicitly
-- give names to intermediate results.
chainif x cfs = foldl (\x (cond, f)-> if cond then f x else x) x cfs

runNetwork :: Map Int Computer -> Set (Integer, Int) -> NAT-> Integer
runNetwork computers heap (xy, blocked, y0)
  | isDone = fromJust y0
  | otherwise = runNetwork computers' heap'' (xy', blocked', y0')
  where
    ((_, i), heap') = S.deleteFindMin heap -- find the computer with lowest time
    c = computers!i
    hasOutput = 3 <= length (output c)
    [addr_, x, y] = take 3 $ output c -- laziness means this won't be executed unless needed
    addr = fromIntegral addr_
    isBlocked = state c == Blocked && null (input c)
    isAllBlocked = isBlocked && isJust xy && all (> 0) (M.elems blocked)
    isDone = isAllBlocked && case (xy, y0) of
      (Nothing, _) -> False
      (_, Nothing) -> False
      (Just (x,y), Just y') -> y == y'
    c0' = (computers!0){input = (\(x,y)->[x,y]) (fromJust xy)}
    c' = chainif c
      [(hasOutput, \c -> c{output = drop 3 $ output c})
      ,(isBlocked, \c -> c{input = [-1]})
      ,(isAllBlocked && i == 0, \c -> c0')
      ,(True, run)
      ]
    c2 = computers!addr
    computers' = chainif computers
      [(hasOutput && addr /= 255, M.insert addr $ c2{input = input c2 ++ [x,y]})
      ,(True, M.insert i c')
      ,(isAllBlocked && i /= 0, M.insert 0 c0')
      ]
    heap'' = S.insert (time c', i) heap'
    xy' = if hasOutput && addr == 255 then Just (x, y) else xy
    blocked' = if not isBlocked then M.insert i 0 blocked
               else let m = M.adjust (+ 1) i blocked
                    in if not isAllBlocked then m
                       else M.insert 0 0 m
    y0' = if isAllBlocked then xy >>= return . snd else y0

main = do
  -- [instructionFile] <- getArgs
  instructionStrings <- readFile "23.input.txt" -- instructionFile
  let instructions = fromList . zip [0 ..] $ map read $ splitOn "," instructionStrings

  putStrLn "Part 2"
  let addrs = [0 .. 49]
      computers = M.fromList [(i, computer0{memory=instructions, input=[fromIntegral i]}) | i <- addrs]
      heap = S.fromList [(0,i) | i <- addrs]
      blocked0 = M.fromList [(i, 0) | i <- addrs]
  print $ runNetwork computers heap (Nothing, blocked0, Nothing)


text:
module NeuralNet.Activation (
  Activation (Sigmoid, Tanh, ReLU),
  forward,
  backward
) where

import Numeric.LinearAlgebra

data Activation = Sigmoid | Tanh | ReLU
  deriving (Eq, Show)

e :: Double
e = exp 1

forward :: Activation -> Matrix Double -> Matrix Double
forward Sigmoid = cmap (\x -> 1 / (1 + (e ** (-x))))
forward Tanh = cmap (\x -> ((e ** x) - (e ** (-x))) / (e ** x + e ** (-x)))
forward ReLU = cmap (max 0)

backward :: Activation -> Matrix Double -> Matrix Double -> Matrix Double
backward Sigmoid dA z = dA * (s * (1 - s))
  where s = cmap (\x -> 1 / (1 + exp x)) (-z)
backward Tanh _ _ = error "Unsure of impl"
backward ReLU dA z = dA * step z

max_stars_repo_name: danielholmes/neural-net

max_issues_repo_name: danielholmes/neural-net

max_forks_repo_name: danielholmes/neural-net


text:
{-# LANGUAGE FlexibleContexts #-}

module Lib.Utility.MBool where

import Numeric.LinearAlgebra

infix  4  .==., ./=., .<., .<=., .>=., .>.
infixr 3  .&&.
infixr 2  .||.

-- specialized for Int result
cond'
    :: (Element t, Ord t, Container c I, Container c t)
    => c t -> c t -> c I -> c I -> c I -> c I
cond' = cond

a .<.  b = cond' a b 1 0 0
a .<=. b = cond' a b 1 1 0
a .==. b = cond' a b 0 1 0
a ./=. b = cond' a b 1 0 1
a .>=. b = cond' a b 0 1 1
a .>.  b = cond' a b 0 0 1

a .&&. b  = step (a*b)
a .||. b  = step (a+b)
no a      = 1-a
xor a b   = a ./=. b
equiv a b = a .==. b
imp a b   = no a .||. b

taut x = minElement x == 1

minEvery a b = cond a b a a b
maxEvery a b = cond a b b b a


text:
{-# LANGUAGE ForeignFunctionInterface #-}
module Grenade.Layers.Internal.Convolution (
    im2col
  , col2im
  , col2vid
  , vid2col
  ) where

import qualified Data.Vector.Storable        as U (unsafeFromForeignPtr0,
                                                   unsafeToForeignPtr0)

import           Foreign                     (mallocForeignPtrArray, withForeignPtr)
import           Foreign.Ptr                 (Ptr)
import           Numeric.LinearAlgebra       (Matrix, cols, flatten, rows)
import qualified Numeric.LinearAlgebra.Devel as U
import           System.IO.Unsafe            (unsafePerformIO)

import           Grenade.Types

col2vid :: Int -> Int -> Int -> Int -> Int -> Int -> Matrix RealNum -> Matrix RealNum
col2vid kernelRows kernelColumns strideRows strideColumns height width dataCol =
  let channels = cols dataCol `div` (kernelRows * kernelColumns)
  in  col2im_c channels height width kernelRows kernelColumns strideRows strideColumns dataCol

col2im :: Int -> Int -> Int -> Int -> Int -> Int -> Matrix RealNum -> Matrix RealNum
col2im kernelRows kernelColumns strideRows strideColumns height width dataCol =
  let channels = 1
  in  col2im_c channels height width kernelRows kernelColumns strideRows strideColumns dataCol

col2im_c :: Int -> Int -> Int -> Int -> Int -> Int -> Int -> Matrix RealNum -> Matrix RealNum
col2im_c channels height width kernelRows kernelColumns strideRows strideColumns dataCol =
  let vec = flatten dataCol
  in unsafePerformIO $ do
    outPtr <- mallocForeignPtrArray (height * width * channels)
    let (inPtr, _) = U.unsafeToForeignPtr0 vec

    withForeignPtr inPtr $ \inPtr' ->
      withForeignPtr outPtr $ \outPtr' ->
        col2im_cpu inPtr' channels height width kernelRows kernelColumns strideRows strideColumns outPtr'

    let matVec = U.unsafeFromForeignPtr0 outPtr (height * width * channels)
    return $ U.matrixFromVector U.RowMajor (height * channels) width matVec

foreign import ccall unsafe
    col2im_cpu
      :: Ptr RealNum -> Int -> Int -> Int -> Int -> Int -> Int -> Int -> Ptr RealNum -> IO ()

vid2col :: Int -> Int -> Int -> Int -> Int -> Int -> Matrix RealNum -> Matrix RealNum
vid2col kernelRows kernelColumns strideRows strideColumns height width dataVid =
  let channels = rows dataVid `div` height
  in  im2col_c channels height width kernelRows kernelColumns strideRows strideColumns dataVid

im2col :: Int -> Int -> Int -> Int -> Matrix RealNum -> Matrix RealNum
im2col kernelRows kernelColumns strideRows strideColumns dataIm =
  let channels = 1
      height = rows dataIm
      width  = cols dataIm
  in  im2col_c channels height width kernelRows kernelColumns strideRows strideColumns dataIm

im2col_c :: Int -> Int -> Int -> Int -> Int -> Int -> Int -> Matrix RealNum -> Matrix RealNum
im2col_c channels height width kernelRows kernelColumns strideRows strideColumns dataIm =
  let vec    = flatten dataIm
      rowOut = (height - kernelRows) `div` strideRows + 1
      colOut = (width - kernelColumns) `div` strideColumns + 1
      kernelSize      = kernelRows * kernelColumns
      numberOfPatches = rowOut * colOut
  in unsafePerformIO $ do
    outPtr <- mallocForeignPtrArray (numberOfPatches * kernelSize * channels)
    let (inPtr, _) = U.unsafeToForeignPtr0 vec

    withForeignPtr inPtr $ \inPtr' ->
      withForeignPtr outPtr $ \outPtr' ->
        im2col_cpu inPtr' channels height width kernelRows kernelColumns strideRows strideColumns outPtr'

    let matVec = U.unsafeFromForeignPtr0 outPtr (numberOfPatches * kernelSize * channels)
    return $ U.matrixFromVector U.RowMajor numberOfPatches (kernelSize * channels) matVec

foreign import ccall unsafe
    im2col_cpu
      :: Ptr RealNum -> Int -> Int -> Int -> Int -> Int -> Int -> Int -> Ptr RealNum -> IO ()

max_stars_repo_name: schnecki/grenade

max_issues_repo_name: schnecki/grenade

max_forks_repo_name: schnecki/grenade


text:
{-# LANGUAGE TypeApplications #-}
{-# OPTIONS_GHC -fno-warn-type-defaults #-}

module ScalarWaveOld where

import Control.Exception (assert)
import Control.Monad.Loops
import qualified Data.Vector.Unboxed as U
import Data.VectorSpace
import Dual
import EquationsOld
import Foreign
import qualified Numeric.LinearAlgebra as L
import SSSFN
import System.CPUTime

default (Int)

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

bndw :: Floating a => Coord a -> a
bndw (Coord u v) = swamp * cos (pi / 4 * swk * t) * cos (pi / 4 * swk * r)
  where
    t = v + u
    r = v - u
    swk = 1
    swamp = 0.1

guess1 :: Floating a => Coord a -> State a
-- guess1 c = State 1 0 (bndw c)
guess1 c = State 1 0 0

guess :: (L.Container L.Vector a, Floating a) => Grid (State a)
guess = gsample (\(u, v) -> guess1 (Coord u v))

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

eqnInt1 :: Fractional a => Eqn -> Coord a -> State a -> State a -> State a -> State a -> State a -> State a -> a
eqnInt1 eqn c s su sv suv suu svv =
  case eqn of
    EqnG01 -> eqnG01 c s su sv suv
    EqnG22 -> eqnG22 c s su sv suv
    EqnSW2 -> eqnSW2 c s su sv suv
    EqnG00 -> eqnG00 c s su suu
    EqnG11 -> eqnG11 c s sv svv

eqnAxis1 :: Fractional a => Eqn -> Coord a -> State a -> State a -> State a -> State a -> State a -> State a -> a
eqnAxis1 eqn c@(Coord u v) s@(State g q w) su@(State gu qu wu) sv@(State gv qv wv) suv@(State guv quv wuv) suu@(State guu quu wuu) svv@(State gvv qvv wvv) =
  case eqn of
    EqnG01 -> eqnAxisG01 c s su sv suv
    EqnG22 -> eqnAxisG22 c s su sv suv
    EqnSW2 -> eqnAxisSW2 c s su sv suv
    EqnG00 -> eqnAxisG00 c s su sv suu
    EqnG11 -> eqnAxisG11 c s su sv svv

eqnOrig1 :: Floating a => Eqn -> Coord a -> State a -> State a -> State a -> State a -> State a -> State a -> a
eqnOrig1 eqn c s@(State g q w) su sv suv suu svv =
  case eqn of
    EqnG01 -> g - 1
    EqnG22 -> q - 0
    EqnSW2 -> w - bndw c
    EqnG00 -> 0
    EqnG11 -> 0

-- i == 0, u == -1
eqnBndU1 :: (Eq a, Floating a) => Eqn -> Coord a -> State a -> State a -> State a -> State a -> State a -> State a -> a
eqnBndU1 eqn c@(Coord u v) s@(State g q w) _ sv _ _ svv =
  assert (u == -1) $
    case eqn of
      EqnG01 -> eqnG11 c s sv svv
      EqnG22 -> q - 0
      EqnSW2 -> w - bndw c
      EqnG00 -> 0
      EqnG11 -> 0

-- j == 0, v == -1
eqnBndV1 :: (Eq a, Floating a) => Eqn -> Coord a -> State a -> State a -> State a -> State a -> State a -> State a -> a
eqnBndV1 eqn c@(Coord u v) s@(State g q w) su _ _ suu _ =
  assert (v == -1) $
    case eqn of
      EqnG01 -> eqnG00 c s su suu
      EqnG22 -> q - 0
      EqnSW2 -> w - bndw c
      EqnG00 -> 0
      EqnG11 -> 0

eqnEndU1 :: Floating a => Eqn -> Coord a -> State a -> State a -> State a -> State a -> State a -> State a -> a
eqnEndU1 eqn c@(Coord u v) s@(State g q w) _ _ _ _ _ =
  case eqn of
    EqnG01 -> g - 1
    EqnG22 -> q - 0
    EqnSW2 -> w - bndw c
    EqnG00 -> 0
    EqnG11 -> 0

eqnEndV1 :: Floating a => Eqn -> Coord a -> State a -> State a -> State a -> State a -> State a -> State a -> a
eqnEndV1 eqn c@(Coord u v) (State g q w) _ _ _ _ _ =
  case eqn of
    EqnG01 -> g - 1
    EqnG22 -> q - 0
    EqnSW2 -> w - bndw c
    EqnG00 -> 0
    EqnG11 -> 0

equation1 :: (Eq a, Floating a) => Coord Int -> Eqn -> Coord a -> State a -> State a -> State a -> State a -> State a -> State a -> a
equation1 (Coord i j) =
  if
    | i == 0 && j == 0 -> eqnOrig1
    | i == 0 && j == np - 1 -> eqnEndU1
    | i == np - 1 && j == 0 -> eqnEndV1
    | i == 0 -> eqnBndU1
    | j == 0 -> eqnBndV1
    | i == j -> eqnAxis1
    | True -> eqnInt1

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

g2s :: Storable a => Grid a -> State (Grid a)
g2s (Grid xs) =
  let [g, q, w] = L.takesV (replicate 3 (np * np)) xs
   in State (Grid g) (Grid q) (Grid w)

s2g :: Storable a => State (Grid a) -> Grid a
s2g (State (Grid g) (Grid q) (Grid w)) = Grid (L.vjoin [g, q, w])

eqn1 ::
  forall a.
  (Eq a, L.Field a, U.Unbox a) =>
  Eqn ->
  Coord (Grid a) ->
  State (Grid a) ->
  Grid a
eqn1 eqn c s =
  Grid $
    L.fromList
      [ eqnRow i j
        | i <- [0 .. np - 1],
          j <- [0 .. np - 1]
      ]
  where
    su = (derivU #>) <$> s
    sv = (derivV #>) <$> s
    suv = (derivUV #>) <$> s
    suu = (derivUU #>) <$> s
    svv = (derivVV #>) <$> s
    eqnRow :: Int -> Int -> a
    eqnRow i j =
      let at :: Grid a -> a
          at x = x ! (i, j)
       in equation1
            (Coord i j)
            eqn
            (at <$> c)
            (at <$> s)
            (at <$> su)
            (at <$> sv)
            (at <$> suv)
            (at <$> suu)
            (at <$> svv)

eqns ::
  forall a.
  (Eq a, L.Field a, U.Unbox a) =>
  Coord (Grid a) ->
  State (Grid a) ->
  Grid a
eqns c s =
  Grid $
    L.fromList
      [ eqnRow (i, j, eqn)
        | eqn <- [EqnG01, EqnG22, EqnSW2],
          i <- [0 .. np - 1],
          j <- [0 .. np - 1]
      ]
  where
    su = (derivU #>) <$> s
    sv = (derivV #>) <$> s
    suv = (derivUV #>) <$> s
    suu = (derivUU #>) <$> s
    svv = (derivVV #>) <$> s
    eqnRow :: (Int, Int, Eqn) -> a
    eqnRow (i, j, eqn) =
      let at :: Grid a -> a
          at x = x ! (i, j)
       in equation1
            (Coord i j)
            eqn
            (at <$> c)
            (at <$> s)
            (at <$> su)
            (at <$> sv)
            (at <$> suv)
            (at <$> suu)
            (at <$> svv)

eqns5 ::
  forall a.
  (Eq a, L.Field a, U.Unbox a) =>
  Coord (Grid a) ->
  State (Grid a) ->
  Grid a
eqns5 c s =
  Grid $
    L.fromList
      [ eqnRow (i, j, eqn)
        | eqn <- [EqnG01, EqnG22, EqnSW2, EqnG00, EqnG11],
          i <- [0 .. np - 1],
          j <- [0 .. np - 1]
      ]
  where
    su = (derivU #>) <$> s
    sv = (derivV #>) <$> s
    suv = (derivUV #>) <$> s
    suu = (derivUU #>) <$> s
    svv = (derivVV #>) <$> s
    eqnRow :: (Int, Int, Eqn) -> a
    eqnRow (i, j, eqn) =
      let at :: Grid a -> a
          at x = x ! (i, j)
       in equation1
            (Coord i j)
            eqn
            (at <$> c)
            (at <$> s)
            (at <$> su)
            (at <$> sv)
            (at <$> suv)
            (at <$> suu)
            (at <$> svv)

op ::
  forall a.
  (L.Container L.Vector a, Eq a, L.Field a, U.Unbox a) =>
  Coord (Grid a) ->
  State (Grid a) ->
  Op a
op (Coord u v) s =
  Op $
    L.fromRows
      [ jacRow (i, j, eqn)
        | eqn <- [EqnG01, EqnG22, EqnSW2],
          i <- [0 .. np - 1],
          j <- [0 .. np - 1]
      ]
  where
    su = (derivU #>) <$> s
    sv = (derivV #>) <$> s
    suv = (derivUV #>) <$> s
    suu = (derivUU #>) <$> s
    svv = (derivVV #>) <$> s
    jacRow :: (Int, Int, Eqn) -> L.Vector a
    jacRow (i, j, eqn) =
      L.vjoin $
        [ let con :: Grid a -> Dual a
              der :: Grid a -> Dual a
              con x = Dual (x ! (i, j)) 0
              der x = Dual (x ! (i, j)) 1
              scon :: State (Grid a) -> State (Dual a)
              sder :: State (Grid a) -> State (Dual a)
              scon (State g q w) = State (con g) (con q) (con w)
              sder (State g q w) = case var of
                VarG -> State (der g) (con q) (con w)
                VarQ -> State (con g) (der q) (con w)
                VarW -> State (con g) (con q) (der w)
           in ( ( dual $
                    equation1
                      (Coord i j)
                      eqn
                      (Coord (con u) (con v))
                      (sder s)
                      (scon su)
                      (scon sv)
                      (scon suv)
                      (scon suu)
                      (scon svv)
                )
                  `L.scale` (getOp delta L.! (i * np + j))
              )
                `L.add` ( ( dual $
                              equation1
                                (Coord i j)
                                eqn
                                (Coord (con u) (con v))
                                (scon s)
                                (sder su)
                                (scon sv)
                                (scon suv)
                                (scon suu)
                                (scon svv)
                          )
                            `L.scale` (getOp derivU L.! (i * np + j))
                        )
                `L.add` ( ( dual $
                              equation1
                                (Coord i j)
                                eqn
                                (Coord (con u) (con v))
                                (scon s)
                                (scon su)
                                (sder sv)
                                (scon suv)
                                (scon suu)
                                (scon svv)
                          )
                            `L.scale` (getOp derivV L.! (i * np + j))
                        )
                `L.add` ( ( dual $
                              equation1
                                (Coord i j)
                                eqn
                                (Coord (con u) (con v))
                                (scon s)
                                (scon su)
                                (scon sv)
                                (sder suv)
                                (scon suu)
                                (scon svv)
                          )
                            `L.scale` (getOp derivUV L.! (i * np + j))
                        )
                `L.add` ( ( dual $
                              equation1
                                (Coord i j)
                                eqn
                                (Coord (con u) (con v))
                                (scon s)
                                (scon su)
                                (scon sv)
                                (scon suv)
                                (sder suu)
                                (scon svv)
                          )
                            `L.scale` (getOp derivUU L.! (i * np + j))
                        )
                `L.add` ( ( dual $
                              equation1
                                (Coord i j)
                                eqn
                                (Coord (con u) (con v))
                                (scon s)
                                (scon su)
                                (scon sv)
                                (scon suv)
                                (scon suu)
                                (sder svv)
                          )
                            `L.scale` (getOp derivVV L.! (i * np + j))
                        )
          | var <- [VarG, VarQ, VarW]
        ]

op5 ::
  forall a.
  (L.Container L.Vector a, Eq a, L.Field a, U.Unbox a) =>
  Coord (Grid a) ->
  State (Grid a) ->
  Op a
op5 (Coord u v) s =
  Op $
    L.fromRows
      [ jacRow (i, j, eqn)
        | eqn <- [EqnG01, EqnG22, EqnSW2, EqnG00, EqnG11],
          i <- [0 .. np - 1],
          j <- [0 .. np - 1]
      ]
  where
    su = (derivU #>) <$> s
    sv = (derivV #>) <$> s
    suv = (derivUV #>) <$> s
    suu = (derivUU #>) <$> s
    svv = (derivVV #>) <$> s
    jacRow :: (Int, Int, Eqn) -> L.Vector a
    jacRow (i, j, eqn) =
      L.vjoin $
        [ let con :: Grid a -> Dual a
              der :: Grid a -> Dual a
              con x = Dual (x ! (i, j)) 0
              der x = Dual (x ! (i, j)) 1
              scon :: State (Grid a) -> State (Dual a)
              sder :: State (Grid a) -> State (Dual a)
              scon (State g q w) = State (con g) (con q) (con w)
              sder (State g q w) = case var of
                VarG -> State (der g) (con q) (con w)
                VarQ -> State (con g) (der q) (con w)
                VarW -> State (con g) (con q) (der w)
           in ( ( dual $
                    equation1
                      (Coord i j)
                      eqn
                      (Coord (con u) (con v))
                      (sder s)
                      (scon su)
                      (scon sv)
                      (scon suv)
                      (scon suu)
                      (scon svv)
                )
                  `L.scale` (getOp delta L.! (i * np + j))
              )
                `L.add` ( ( dual $
                              equation1
                                (Coord i j)
                                eqn
                                (Coord (con u) (con v))
                                (scon s)
                                (sder su)
                                (scon sv)
                                (scon suv)
                                (scon suu)
                                (scon svv)
                          )
                            `L.scale` (getOp derivU L.! (i * np + j))
                        )
                `L.add` ( ( dual $
                              equation1
                                (Coord i j)
                                eqn
                                (Coord (con u) (con v))
                                (scon s)
                                (scon su)
                                (sder sv)
                                (scon suv)
                                (scon suu)
                                (scon svv)
                          )
                            `L.scale` (getOp derivV L.! (i * np + j))
                        )
                `L.add` ( ( dual $
                              equation1
                                (Coord i j)
                                eqn
                                (Coord (con u) (con v))
                                (scon s)
                                (scon su)
                                (scon sv)
                                (sder suv)
                                (scon suu)
                                (scon svv)
                          )
                            `L.scale` (getOp derivUV L.! (i * np + j))
                        )
                `L.add` ( ( dual $
                              equation1
                                (Coord i j)
                                eqn
                                (Coord (con u) (con v))
                                (scon s)
                                (scon su)
                                (scon sv)
                                (scon suv)
                                (sder suu)
                                (scon svv)
                          )
                            `L.scale` (getOp derivUU L.! (i * np + j))
                        )
                `L.add` ( ( dual $
                              equation1
                                (Coord i j)
                                eqn
                                (Coord (con u) (con v))
                                (scon s)
                                (scon su)
                                (scon sv)
                                (scon suv)
                                (scon suu)
                                (sder svv)
                          )
                            `L.scale` (getOp derivVV L.! (i * np + j))
                        )
          | var <- [VarG, VarQ, VarW]
        ]

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

getTime :: IO Double
getTime = do
  t <- getCPUTime
  return $ fromInteger t / 1.0e+12

main :: IO ()
main =
  do
    putStrLn "Spherically Symmetric Scalar Field in double-Null coordinates"
    let coords = Coord coordU (coordV @Double)
    let ini = State (gmap g guess) (gmap q guess) (gmap w guess)
    vars0 <- return (s2g ini)
    let done (n, vars) =
          let res = eqns coords (g2s vars)
           in -- let res = eqns5 coords (g2s vars)
              gmaxabs res < 1.0e-12 || n >= 5
    let iter (n, vars) = do
          let res = eqns coords (g2s vars)
          -- let res = eqns5 coords (g2s vars)
          t <- getTime
          putStrLn $
            "iter " ++ show n ++ " time " ++ show (round t) ++ " sec |res| "
              ++ show (gmaxabs res)
          let jac = op coords (g2s vars)
          let dvars = solve jac res
          -- let jac = op5 coords (g2s vars)
          -- let dvars = solveLS jac res
          return (n + 1, vars ^-^ dvars)
    (n1, vars1) <- iterateUntilM done iter (0, vars0)
    putStrLn $ "var.g " ++ show (gmap approx $ g $ g2s vars1)
    putStrLn $ "var.q " ++ show (gmap approx $ q $ g2s vars1)
    putStrLn $ "var.w " ++ show (gmap approx $ w $ g2s vars1)
    putStrLn $ "res.G01 " ++ show (gmap approx $ eqn1 EqnG01 coords (g2s vars1))
    putStrLn $ "res.G22 " ++ show (gmap approx $ eqn1 EqnG22 coords (g2s vars1))
    putStrLn $ "res.SW2 " ++ show (gmap approx $ eqn1 EqnSW2 coords (g2s vars1))
    putStrLn $ "res.G00 " ++ show (gmap approx $ eqn1 EqnG00 coords (g2s vars1))
    putStrLn $ "res.G11 " ++ show (gmap approx $ eqn1 EqnG11 coords (g2s vars1))
    let res1 = eqns coords (g2s vars1)
    t <- getTime
    putStrLn $
      "iter " ++ show n1 ++ " time " ++ show (round t) ++ " sec |res| "
        ++ show (gmaxabs res1)
    let res5 = eqns5 coords (g2s vars1)
    putStrLn $ "|res5| " ++ show (gmaxabs res5)
    -- let h = 2 / (fromIntegral np - 1) :: Double
    -- let alpha = 0.01 * h ^ 3
    -- let smooth = delta - alpha *^ (derivUU * derivUU + derivVV * derivVV)
    -- let vars1s = s2g $ (expfilter #>) <$> g2s vars1
    -- let res5s = eqns5 coords (g2s vars1s)
    -- putStrLn $ "|res5s| " ++ show (gmaxabs res5s)
    putStrLn "Done."
  where
    approx :: RealFrac a => a -> a
    approx x = fromInteger (round (1.0e+10 * x)) / 1.0e+10

max_stars_repo_name: eschnett/sssfn

max_issues_repo_name: eschnett/sssfn

max_forks_repo_name: eschnett/sssfn


text:
module HLearn.Optimization.Amoeba
    where

import Prelude
import Control.Monad
import Control.Monad.ST
import Data.List
import Debug.Trace
import qualified Data.Vector as V
import qualified Data.Vector.Mutable as VM
import qualified Data.Vector.Storable as VS
import qualified Data.Vector.Storable.Mutable as VSM
import qualified Data.Vector.Generic as VG
import qualified Data.Vector.Generic.Mutable as VGM
import qualified Data.Vector.Algorithms.Intro as Intro
import Numeric.LinearAlgebra hiding ((<>))
import qualified Numeric.LinearAlgebra as LA

findMinAmoeba f x0 = runST $ do

    -- initialize simplex
    vec <- VM.new (VG.length x0+1)
    VGM.write vec 0 (f x0,x0)
    forM [1..VGM.length vec-1] $ \i -> do
        e_i <- VGM.replicate (VG.length x0) 0
        VGM.write e_i (i-1) 1
        e_i' <- VG.freeze e_i
        let x_i = x0 `LA.add` e_i'
        VGM.write vec i (f x_i,x_i)

    -- iterate
    vec' <- itrM 1000 (stepAmoeba f) vec

    -- return
    (_,ret) <- VGM.read vec 0
    return ret

stepAmoeba f vec = stepAmoebaRaw 1 2 (-1/2) (1/2) f vec

stepAmoebaRaw ::
    ( Fractional b
    , Ord b
    , VGM.MVector vec (b,a)
--     , a ~ LA.Matrix b
    , a ~ LA.Vector b
    , Field b
    , vec ~ VM.MVector
    , b ~ Double
    ) => b
      -> b
      -> b
      -> b
      -> (a -> b)
      -> vec s (b,a)
      -> ST s (vec s (b,a))
stepAmoebaRaw alpha gamma ro sigma f vec = do

    Intro.sortBy (\a b -> compare (fst a) (fst b)) vec

    (f_1,x_1) <- VGM.read vec 0
    (f_2,x_2) <- VGM.read vec 1
    (f_n1,x_n1) <- VGM.read vec $ VGM.length vec -1

    x_0 <- liftM ( scale (1/fromIntegral (VGM.length vec-1))
                 . foldl1' (LA.add)
                 . init
                 . map snd
                 . V.toList
                 ) $ VG.unsafeFreeze vec

    let x_r = x_0 `LA.add` (scale alpha $ x_0 `LA.sub` x_n1)
        f_r = f x_r

        x_e = x_0 `LA.add` (scale gamma $ x_0 `LA.sub` x_n1)
        f_e = f x_e

        x_c = x_0 `LA.add` (scale ro $ x_0 `LA.sub` x_n1)
        f_c = f x_c

    -- check reflection
    if f_1 <= f_r && f_r < f_1
        then VGM.write vec (VGM.length vec-1) (f_r,x_r)

        -- check expansion
        else if f_r < f_1
            then if f_e < f_r
                then VGM.write vec (VGM.length vec-1) (f_e,x_e)
                else VGM.write vec (VGM.length vec-1) (f_r,x_r)

            -- check contraction
            else if f_c < f_n1
                then VGM.write vec (VGM.length vec-1) (f_c,x_c)

                -- reduction
                else forM_ [1..VGM.length vec-1] $ \i -> do
                    (f_i,x_i) <- VGM.read vec i
                    let x_i' = x_1 `LA.add` (scale sigma $ x_i `LA.sub` x_1)
                        f_i' = f x_i'
                    VGM.write vec i (f_i',x_i')

    return vec
--     refMinVal <- newSTRef (-infinity)
--     refMinIndex <- newSTRef 0
--     forM [0..VGM.length vec-1] $ \i -> do
--         ival <- VGM.read vec i
--         minVal <- readSTRef refMinVal
--         if minVal < fst ival
--             then return ()
--             else do
--                 writeSTRef refMinVal $ fst ival
--                 writeSTRef refMinIndex i
--     undefined

itrM :: Monad m => Int -> (a -> m a) -> a -> m a
itrM 0 f a = return a
itrM i f a = do
    a' <- f a
--     if a' == a
--         then trace ("no movement\n   a="++show a++"\n   a'="++show a') $ return ()
--         else trace ("itrM i="++show i++"; f a ="++show a'++"; a="++show a) $ return ()
    itrM (i-1) f a'

max_stars_repo_name: Heather/HLearn

max_issues_repo_name: Heather/HLearn

max_forks_repo_name: Heather/HLearn