text:
{-# LANGUAGE NamedFieldPuns, TemplateHaskell #-}

module School.Unit.Test.Affine
( affineTest ) where

import Control.Monad.IO.Class (liftIO)
import Numeric.LinearAlgebra ((><), (|>), IndexOf, Matrix, R, Vector, accum,
                              fromRows, ident, size)
import School.TestUtils (CostFunc, diffInput, jTest, matIndexes,
                         randomAffineParams, randomMatrix)
import School.Types.FloatEq ((~=))
import School.Unit.Affine
import School.Unit.Unit (Unit(..))
import School.Unit.UnitActivation (UnitActivation(..))
import School.Unit.UnitGradient (UnitGradient(..))
import School.Unit.UnitParams (UnitParams(..))
import School.Utils.LinearAlgebra (compareDoubleMatrix, compareDoubleVector,
                                   oneMatrix, zeroMatrix, zeroVector)
import Test.Tasty (TestTree)
import Test.Tasty.QuickCheck hiding ((><), scale)
import Test.Tasty.TH
import Test.QuickCheck.Monadic (assert, monadicIO, pre)

eps :: Double
eps = 1e-5

prec :: Double
prec = 5e-2

prop_trivial_apply :: (Positive Int) -> (Positive Int) -> Property
prop_trivial_apply (Positive bSize) (Positive fSize) = monadicIO $ do
  let affineBias = zeroVector fSize
  m <- liftIO $ randomMatrix bSize fSize
  let input = BatchActivation m
  let affineWeights = ident fSize
  let params = AffineParams { affineBias
                            , affineWeights }
  let result = apply affine params input
  assert $ result ~= input

randomSetup :: Int
            -> Int
            -> Int
            -> IO (UnitParams R, UnitActivation R)
randomSetup bSize fSize oSize = do
  input <- randomMatrix bSize fSize
  affineParams <- randomAffineParams fSize oSize
  return ( affineParams
         , BatchActivation input
         )

prop_overwrite_apply :: (Positive Int) -> (Positive Int) -> Positive Int -> Property
prop_overwrite_apply (Positive bSize) (Positive fSize) (Positive oSize) = monadicIO $ do
  let affineWeights = zeroMatrix oSize fSize
  (params, input) <- liftIO $ randomSetup bSize fSize oSize
  let params' = params { affineWeights }
  let result = apply affine params' input
  let check = fromRows . (replicate bSize) $ affineBias params
  assert $ result ~= BatchActivation check

prop_deriv_dims :: (Positive Int) -> (Positive Int) -> Positive Int -> Property
prop_deriv_dims (Positive bSize) (Positive fSize) (Positive oSize) = monadicIO $ do
  (params, input) <- liftIO $ randomSetup bSize fSize oSize
  gradIn <- liftIO $ fmap BatchGradient $ randomMatrix bSize oSize
  let (BatchGradient gradOut, AffineParams{ affineBias, affineWeights }) = deriv affine params gradIn input
  assert $ (size gradOut == (bSize, fSize))
        && (size affineBias == oSize)
        && (size affineWeights == (oSize, fSize))

type AlterParams a = Double
                  -> IndexOf a
                  -> UnitParams R
                  -> UnitParams R

diffBias :: UnitParams R
          -> UnitActivation R
          -> CostFunc
          -> IndexOf Vector
          -> Double
diffBias params input cost idx =
  (jAdd - jSub) / (2*eps) where
    outAdd = apply affine (alterBias eps idx params) input
    outSub = apply affine (alterBias (-eps) idx params) input
    jAdd = cost outAdd
    jSub = cost outSub

diffWeights :: UnitParams R
            -> UnitActivation R
            -> CostFunc
            -> IndexOf Matrix
            -> Double
diffWeights params input cost idx =
  (jAdd - jSub) / (2*eps) where
    outAdd = apply affine (alterWeights eps idx params) input
    outSub = apply affine (alterWeights (-eps) idx params) input
    jAdd = cost outAdd
    jSub = cost outSub

alterBias :: AlterParams Vector
alterBias change idx AffineParams{ affineBias, affineWeights } =
  let newBias = accum affineBias (+) [(idx, change)]
  in AffineParams { affineBias = newBias, affineWeights }
alterBias _ _ _ = EmptyParams

alterWeights :: AlterParams Matrix
alterWeights change idx AffineParams{ affineBias, affineWeights } =
  let newWeights = accum affineWeights (+) [(idx, change)]
  in AffineParams { affineBias, affineWeights = newWeights }
alterWeights _ _ _ = EmptyParams

vecIndexes :: Int -> [IndexOf Vector]
vecIndexes n = [0..n-1]

prop_numerical_gradient :: (Positive Int) -> (Positive Int) -> Positive Int -> Property
prop_numerical_gradient (Positive bSize) (Positive fSize) (Positive oSize) = monadicIO $ do
  pre $ bSize < 50 && fSize < 50 && oSize < 50
  (params, input) <- liftIO $ randomSetup bSize fSize oSize
  let idxs = matIndexes bSize fSize
  let num = map (\idx -> diffInput affine params input eps idx)
                idxs
  let numGrad = (bSize >< fSize) num
  let inGrad = BatchGradient $ oneMatrix bSize oSize
  let (BatchGradient result, _) = deriv affine params inGrad input
  assert $ compareDoubleMatrix prec result numGrad

prop_numerical_bias_deriv :: (Positive Int) -> (Positive Int) -> Positive Int -> Property
prop_numerical_bias_deriv (Positive bSize) (Positive fSize) (Positive oSize) = (mapSize (const 10)) . monadicIO $ do
  pre $ bSize < 50 && fSize < 50 && oSize < 50
  (params, input) <- liftIO $ randomSetup bSize fSize oSize
  let idxs = vecIndexes oSize
  let num = map (\idx -> diffBias params input jTest idx)
                idxs
  let numDeriv = oSize |> num
  let inGrad = BatchGradient $ oneMatrix bSize oSize
  let (_, AffineParams { affineBias }) = deriv affine params inGrad input
  assert $ compareDoubleVector prec affineBias numDeriv

prop_numerical_weights_deriv :: (Positive Int) -> (Positive Int) -> Positive Int -> Property
prop_numerical_weights_deriv (Positive bSize) (Positive fSize) (Positive oSize) = (mapSize (const 10)) . monadicIO $ do
  pre $ bSize < 50 && fSize < 50 && oSize < 50
  (params, input) <- liftIO $ randomSetup bSize fSize oSize
  let idxs = matIndexes oSize fSize
  let num = map (\idx -> diffWeights params input jTest idx)
                idxs
  let numDeriv = (oSize >< fSize) num
  let inGrad = BatchGradient $ oneMatrix bSize oSize
  let (_, AffineParams { affineWeights }) = deriv affine params inGrad input
  assert $ compareDoubleMatrix prec affineWeights numDeriv

affineTest :: TestTree
affineTest = $(testGroupGenerator)

max_stars_repo_name: jfulseca/School

max_issues_repo_name: jfulseca/School

max_forks_repo_name: jfulseca/School


text:
{-|
Module: MachineLearning.Optimization
Description: Optimization
Copyright: (c) Alexander Ignatyev, 2016
License: BSD-3
Stability: experimental
Portability: POSIX

Optimization module.
-}

module MachineLearning.Optimization
(
    MinimizeMethod(..)
  , minimize
  , checkGradient
)

where

import MachineLearning.Types (R, Vector, Matrix)
import MachineLearning.Model as Model
import MachineLearning.Regularization (Regularization)
import qualified MachineLearning.Optimization.GradientDescent as GD
import qualified MachineLearning.Optimization.MinibatchGradientDescent as MGD
import qualified Data.Vector.Storable as V
import qualified Numeric.LinearAlgebra as LA

import qualified Numeric.GSL.Minimization as Min

data MinimizeMethod = GradientDescent R       -- ^ Gradient descent, takes alpha. Requires feature normalization.
                    | MinibatchGradientDescent Int Int R  -- ^ Minibacth Gradietn Descent, takes seed, batch size and alpha
                    | ConjugateGradientFR R R -- ^ Fletcher-Reeves conjugate gradient algorithm,
                                              -- takes size of first trial step (0.1 is fine) and tol (0.1 is fine).
                    | ConjugateGradientPR R R -- ^ Polak-Ribiere conjugate gradient algorithm.
                                              -- takes size of first trial step (0.1 is fine) and tol (0.1 is fine).
                    | BFGS2 R R               -- ^ Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm,
                                              -- takes size of first trial step (0.1 is fine) and tol (0.1 is fine).

-- | Returns solution vector (theta) and optimization path.
-- Optimization path's row format:
-- [iter number, cost function value, theta values...]
minimize :: Model.Model a =>
            MinimizeMethod
         -> a       -- ^ model (Least Squares, Logistic Regression etc)
         -> R   -- ^ epsilon, desired precision of the solution
         -> Int     -- ^ maximum number of iterations allowed
         -> Regularization   -- ^ regularization parameter
         -> Matrix  -- ^ X
         -> Vector  -- ^ y
         -> Vector  -- ^ initial solution, theta
         -> (Vector, Matrix) -- ^ solution vector and optimization path

minimize (BFGS2 firstStepSize tol) = minimizeVD Min.VectorBFGS2 firstStepSize tol
minimize (ConjugateGradientFR firstStepSize tol) = minimizeVD Min.ConjugateFR firstStepSize tol
minimize (ConjugateGradientPR firstStepSize tol) = minimizeVD Min.ConjugatePR firstStepSize tol
minimize (GradientDescent alpha) = GD.gradientDescent alpha
minimize (MinibatchGradientDescent seed batchSize alpha) = MGD.minibatchGradientDescent seed batchSize alpha

minimizeVD method firstStepSize tol model epsilon niters reg x y initialTheta
  = Min.minimizeVD method epsilon niters firstStepSize tol costf gradientf initialTheta
  where costf = Model.cost model reg x y
        gradientf = Model.gradient model reg x y

-- | Gradient checking function.
-- Approximates the derivates of the Model's cost function
-- and calculates derivatives using the Model's gradient functions.
-- Returns norn_2 between 2 derivatives.
-- Takes model, regularization, X, y, theta and epsilon (used to approximate derivatives, 1e-4 is a good value).
checkGradient :: Model a => a -> Regularization -> Matrix -> Vector -> Vector -> R -> R
checkGradient model reg x y theta eps = LA.norm_2 $ grad1 - grad2
  where theta_m = LA.asColumn theta
        eps_v = V.replicate (V.length theta) eps
        eps_m = LA.diag eps_v
        -- +/- eps_m in case of zero theta
        theta_m1 = theta_m * (1 + eps_m) + eps_m
        theta_m2 = theta_m * (1 - eps_m) - eps_m
        cost1 = LA.vector $ map (cost model reg x y) $ LA.toColumns theta_m1
        cost2 = LA.vector $ map (cost model reg x y) $ LA.toColumns theta_m2
        eps_v1 = LA.takeDiag $ theta_m1 - theta_m
        eps_v2 = LA.takeDiag $ theta_m - theta_m2
        grad1 = (cost1 - cost2) / (eps_v1 + eps_v2)
        grad2 = gradient model reg x y theta

ext: hs

max_stars_repo_name: aligusnet/mltool

max_issues_repo_name: aligusnet/mltool

max_forks_repo_name: aligusnet/mltool


text:
module Instruction where

import           Protolude

import           Data.Complex (Complex (..))
import qualified Data.Matrix  as M

-- NOTE: could be newtyped for more safety
type Qubit = Int

type Qubits = [Qubit]

type Operation = M.Matrix (Complex Double)

data Instruction
  = Gate Operation Qubits
  | Measure
  | Halt
  deriving (Show)

max_stars_repo_name: colescott/quantum-interpreter

max_issues_repo_name: colescott/quantum-interpreter

max_forks_repo_name: colescott/quantum-interpreter


text:
{-
 - Name: Benjamin Kostiuk
 - Date: 9/25/2018
 -}
 -- Module for computing all roots of a polynomial of the third degree
module Assign_1_ExtraCredit where
import Data.Complex
-- see https://www.stackage.org/haddock/lts-8.24/base-4.9.1.0/Data-Complex.html

{- -----------------------------------------------------------------
 - cubicRoot
 - -----------------------------------------------------------------
 - Description: Computes the cubic root of a Double
 -}
cubicRoot :: Double -> Double
cubicRoot n
    | (n < 0) = - ((abs n) ** (1/3))
    | otherwise = n ** (1/3)

{- -----------------------------------------------------------------
 - cubicQ
 - -----------------------------------------------------------------
 - Description: Calculates the value of Q from coefficients a, b and c
 -}
cubicQ :: Double -> Double -> Double -> Double
cubicQ a b c = (3*a*c - b^2) / (9*a^2)

{- -----------------------------------------------------------------
 - cubicR
 - -----------------------------------------------------------------
 - Description: Calculates the value of R from coefficients a, b, c and d
 -}
cubicR :: Double -> Double -> Double -> Double -> Double
cubicR a b c d = (9*a*b*c - 27*a^2*d - 2*b^3) / (54*a^3)

{- -----------------------------------------------------------------
 - cubicDisc
 - -----------------------------------------------------------------
 - Description: Calculates the value of the discriminant from R and Q
 -}
cubicDisc :: Double -> Double -> Double
cubicDisc q r = q^3 + r^2

{- -----------------------------------------------------------------
 - cubicComplexS
 - -----------------------------------------------------------------
 - Description: Computes a complex double from the value of S from Q and R
 - The value of S is only calculated when the discriminant >= 0 because it is not needed for discriminant < 0
 -}
cubicComplexS :: Double -> Double -> Complex Double
cubicComplexS q r
    | (disc >= 0) = cubicRoot(r + sqrt(disc)) :+ 0
    | otherwise = 0 :+ 0 -- If disc < 0, not needed for calculations --
    where
        disc = cubicDisc q r

{- -----------------------------------------------------------------
 - cubicComplexT
 - -----------------------------------------------------------------
 - Description: Computes a complex double from the value of T from Q and R
 - The value of T is only calculated when the discriminant >= 0 because it is not needed for discriminant < 0
 -}
cubicComplexT :: Double -> Double -> Complex Double
cubicComplexT q r
    | (disc >= 0) = cubicRoot(r - sqrt(disc)) :+ 0
    | otherwise = 0 :+ 0 -- If disc < 0, not needed for calculations --
    where
        disc = cubicDisc q r

{- -----------------------------------------------------------------
 - cmplxRealSolutions
 - -----------------------------------------------------------------
 - Description: Computes a list of real roots from the angular components of a cubic equation
 - Recursively returns the list of real roots in complex double form
 -}
cmplxRealSolutions :: [Complex Double] -> Double -> Double -> Double -> Double -> Double -> [Complex Double]
cmplxRealSolutions xs n a b q theta
    | (length xs == 3) = xs
    | otherwise = cmplxRealSolutions (xs++[x]) (n+2) a b q theta    -- n is incremented by two to compute all angular solutions --
    where
        x = (2 * sqrt(-q) * cos((theta / 3) + (n*pi/3)) - (b / (3*a))) :+ 0

{- -----------------------------------------------------------------
 - cubicComplexSolutions
 - -----------------------------------------------------------------
 - Description: Computes a list of complex double roots of a cubic equation
 - If discriminant > 0, returns a list with one real and two imaginary roots
 - If discriminant equal to 0, returns a list with three roots, with x2 = x3
 - If discriminant < 0, returns a list with three distinct real roots
 -}
cubicComplexSolutions :: Double -> Double -> Double -> Double -> [Complex Double]
cubicComplexSolutions a b c d
    | ((abs disc) < 1e-9) = [x1 :+ 0, x2 :+ 0, x2 :+ 0]
    | (disc > 0)          = [x1 :+ 0, x2 :+ (((sqrt 3) / 2) * (s - t)), x2 :+ (-((sqrt 3) / 2) * (s - t))]
    | otherwise           = cmplxRealSolutions [] 0 a b q theta
    where
        q = cubicQ a b c
        r = cubicR a b c d
        s = realPart (cubicComplexS q r)
        t = realPart (cubicComplexT q r)
        disc = cubicDisc q r
        theta = acos(r / sqrt(-q^3)) -- Angle used in angular computation --
        x1 = ((s + t) - (b / (3*a)))
        x2 = (- ((s + t) / 2) - (b / (3*a)))

{- -----------------------------------------------------------------
 - Test Cases
 - -----------------------------------------------------------------
 -}
{-
-- See test cases from Assign_1.hs for any other function

--cubicComplexSolutions--
    Input: 3 4 2 1
    Expected: [-1.0 :+ 0.0, -0.1666667 :+ 0.55277, -0.1666667 :+ (-0.55277)]
    Actual: [(-1.0000000000000004) :+ 0.0,(-0.1666666666666664) :+ 0.5527707983925663,(-0.1666666666666664) :+ (-0.5527707983925663)]

    Input: 0.00033 0.00077 (-0.0002) 0.0001
    Expected: [(-2.6100218) :+ 0.0, 0.1383442795 :+ 0.31138954,  0.1383442795 :+ (-0.31138954)]
    Actual: [(-2.6100218924147884) :+ 0.0,0.1383442795407276 :+ 0.31138954838120664,0.1383442795407276 :+ (-0.31138954838120664)]

    Input: 1 (-3) 3 (-1)
    Expected: [1.0 :+ 0.0, 1.0 :+ 0.0, 1.0 :+ 0.0]
    Actual: [1.0 :+ 0.0,1.0 :+ 0.0,1.0 :+ 0.0]

    Input: 1 (-5) 8 (-4)
    Expected: [1.0 :+ 0.0, 2.0 :+ 0.0, 2.0 :+ 0.0]
    Actual: [1.0 :+ 0.0,2.0 :+ 0.0,2.0 :+ 0.0]

    Input: 1 0 (-3) 0
    Expected: [0.0 :+ 0.0, 1.7320508 :+ 0.0, (-1.7320508) :+ 0.0]
    Actual: [1.7320508075688774 :+ 0.0,(-1.732050807568877) :+ 0.0,(-3.673819061467132e-16) :+ 0.0]

    Input: 1 (-5) (-2) 24
    Expected: [4.0 :+ 0.0, (-2.0) :+ 0.0, 3.0 :+ 0.0]
    Actual: [4.000000000000001 :+ 0.0,(-1.9999999999999993) :+ 0.0,2.999999999999999 :+ 0.0]

    Input: 1 0 (-7) (-6)
    Expected: [3.0 :+ 0.0, (-2.0) :+ 0.0, (-1.0) :+ 0.0]
    Actual: [3.0000000000000004 :+ 0.0,(-1.9999999999999993) :+ 0.0,(-1.0000000000000002) :+ 0.0]

    Input: (-0.004) (-0.0122) 17 3
    Expected: [63.774967 :+ 0.0, (-66.6485177) :+ 0.0, (-0.17644953) :+ 0.0]
    Actual: [63.77496726734746 :+ 0.0,(-66.64851773002562) :+ 0.0,(-0.17644953732185886) :+ 0.0]
 -}

max_stars_repo_name: BenjaminKostiuk/haskell-slide-rule

max_issues_repo_name: BenjaminKostiuk/Code-Booklet

max_forks_repo_name: BenjaminKostiuk/Code-Booklet


text:
module LogisticRegression where

import Numeric.LinearAlgebra

sigmoid :: Matrix Double -> Matrix Double
sigmoid x = fromRows $ map (cmap sigmoid') (toRows x)

-- sigmoid applied to a single real number.
sigmoid' :: Double -> Double
sigmoid' x = 1 / (1 + (exp (-1 * x)))

max_stars_repo_name: mchav/CourseraMachineLearning

max_issues_repo_name: mchav/CourseraMachineLearning

max_forks_repo_name: mchav/CourseraMachineLearning


text:
{-# LANGUAGE DefaultSignatures   #-}
{-# LANGUAGE FlexibleContexts    #-}
{-# LANGUAGE InstanceSigs        #-}
{-# LANGUAGE KindSignatures      #-}
{-# LANGUAGE PolyKinds           #-}
{-# LANGUAGE RankNTypes          #-}
{-# LANGUAGE ScopedTypeVariables #-}
module Data.Roles
  ( Representational(rep)
  , Phantom(phantom)

  -- * Utility functions to define Representable easily
  , Proxy2(..)
  , wrapped
  , slash
  , coercionOf
  , eta
  ) where

import           Control.Applicative
import           Control.Category      ((>>>))
import           Control.Monad.ST      (ST)

import           Data.Complex
import           Data.Monoid
import           Data.Semigroup

import           Data.Proxy

import           Data.Functor.Identity
import           Data.Functor.Compose
import qualified Data.Functor.Product  as F
import qualified Data.Functor.Sum      as F

import           Data.IntMap
import           Data.Map

import Control.Monad.Reader      (ReaderT(..))
import Control.Monad.Except      (ExceptT(..), runExceptT)
import Control.Monad.Cont        (ContT(..))
import Control.Monad.Trans.Maybe (MaybeT(..))
import qualified Control.Monad.Writer.Lazy   as L (WriterT(..), runWriterT)
import qualified Control.Monad.Writer.Strict as S (WriterT(..), runWriterT)
import qualified Control.Monad.State.Lazy    as L (StateT(..), runStateT)
import qualified Control.Monad.State.Strict  as S (StateT(..), runStateT)
import qualified Control.Monad.RWS.Lazy      as L (RWST(..), runRWST)
import qualified Control.Monad.RWS.Strict    as S (RWST(..), runRWST)

import           Data.Coerce
import           Data.Type.Coercion
import           Unsafe.Coerce

class Representational (t :: k1 -> k2) where
  -- | An argument is representational if you can lift a coercion of the argument into one of the whole
  rep :: Coercion a b -> Coercion (t a) (t b)
  default rep :: Phantom t => Coercion a b -> Coercion (t a) (t b)
  rep _ = phantom

class Representational t => Phantom (t :: k1 -> k2) where
  -- | An argument is phantom if you can 'coerce' the whole ignoring the argument
  phantom :: Coercion (t a) (t b)
  default phantom :: Coercible (t a) (t b) => Coercion (t a) (t b)
  phantom = Coercion

wrapped :: forall x y proxy (f :: x -> y) (g :: x -> y) (a :: x) (b :: x).
           (Coercible f g) => proxy f g -> Coercion (f a) (f b) -> Coercion (g a) (g b)
wrapped _ c = (Coercion :: Coercion (g a) (f a)) >>> c >>> (Coercion :: Coercion (f b) (g b))

data Proxy2 a b = Proxy2

slash ::  proxy (f a) (g a) -> Proxy2 f g
slash _ = Proxy2

coercionOf :: Coercible a b => proxy a b -> Coercion a b
coercionOf _ = Coercion

eta :: forall f g a. Coercion f g -> Coercion (f a) (g a)
eta Coercion = Coercion

-- * Data.Proxy

instance Representational Proxy
instance Phantom Proxy

-- * Data.Functor.{Const, Identity, Sum, Product, Compose}

instance Representational Const where rep Coercion = Coercion
instance Representational (Const a)
instance Phantom (Const a)

instance Representational Identity where rep Coercion = Coercion

instance Representational F.Sum     where rep Coercion = Coercion
instance Representational (F.Sum f) where rep Coercion = Coercion

instance Representational F.Product where rep Coercion = Coercion
instance Representational (F.Product f) where rep Coercion = Coercion

instance Representational Compose where rep Coercion = Coercion
instance (Representational f,
          Representational g)
      => Representational (Compose f g) where
  rep c = coercionOf getCompose >>> rep (rep c) >>> coercionOf Compose

-- * Data.Type.Coercion

instance Representational Coercion     where rep = unsafeCoerce
instance Representational (Coercion a) where rep Coercion = Coercion

-- * Prelude

instance Representational (->)       where rep Coercion = Coercion
instance Representational ((->) a)   where rep Coercion = Coercion

instance Representational Either     where rep Coercion = Coercion
instance Representational (Either a) where rep Coercion = Coercion

instance Representational (,)     where rep Coercion = Coercion
instance Representational ((,) a) where rep Coercion = Coercion

instance Representational (,,)       where rep Coercion = Coercion
instance Representational ((,,) a)   where rep Coercion = Coercion
instance Representational ((,,) a b) where rep Coercion = Coercion

instance Representational (,,,)         where rep Coercion = Coercion
instance Representational ((,,,) a)     where rep Coercion = Coercion
instance Representational ((,,,) a b)   where rep Coercion = Coercion
instance Representational ((,,,) a b c) where rep Coercion = Coercion

instance Representational (,,,,)           where rep Coercion = Coercion
instance Representational ((,,,,) a)       where rep Coercion = Coercion
instance Representational ((,,,,) a b)     where rep Coercion = Coercion
instance Representational ((,,,,) a b c)   where rep Coercion = Coercion
instance Representational ((,,,,) a b c d) where rep Coercion = Coercion

instance Representational (,,,,,)             where rep Coercion = Coercion
instance Representational ((,,,,,) a)         where rep Coercion = Coercion
instance Representational ((,,,,,) a b)       where rep Coercion = Coercion
instance Representational ((,,,,,) a b c)     where rep Coercion = Coercion
instance Representational ((,,,,,) a b c d)   where rep Coercion = Coercion
instance Representational ((,,,,,) a b c d e) where rep Coercion = Coercion

instance Representational []      where rep Coercion = Coercion
instance Representational Maybe   where rep Coercion = Coercion
instance Representational IO      where rep Coercion = Coercion
instance Representational (ST s)  where rep Coercion = Coercion

-- * Control.Applicative

instance Representational ZipList where rep Coercion = Coercion

instance Representational WrappedMonad where rep Coercion = Coercion
instance Representational m => Representational (WrappedMonad m) where
  rep c = wrapped (slash WrapMonad) (rep c)

instance Representational WrappedArrow where rep Coercion = Coercion
instance (Representational p) => Representational (WrappedArrow p) where
  rep c = wrapped (slash (slash WrapArrow)) (rep c)
instance (Representational (p a)) => Representational (WrappedArrow p a) where
  rep c = wrapped (slash WrapArrow) (rep c)

-- * Data.Complex

instance Representational Complex where rep Coercion = Coercion

-- * Data.Semigroup

instance Representational Arg     where rep Coercion = Coercion
instance Representational (Arg a) where rep Coercion = Coercion
instance Representational Product where rep Coercion = Coercion
instance Representational Sum     where rep Coercion = Coercion
instance Representational Dual    where rep Coercion = Coercion
instance Representational Endo    where rep Coercion = Coercion
instance Representational Option  where rep Coercion = Coercion

instance Representational Data.Semigroup.First where rep Coercion = Coercion
instance Representational Data.Semigroup.Last  where rep Coercion = Coercion
instance Representational WrappedMonoid        where rep Coercion = Coercion

-- * Data.Monoid

instance Representational Ap                  where rep Coercion = Coercion
instance Representational f => Representational (Ap f) where
  rep c = wrapped (slash Ap) (rep c)

instance Representational Data.Monoid.First   where rep Coercion = Coercion
instance Representational Data.Monoid.Last    where rep Coercion = Coercion

-- * containers
instance Representational (Map k) where rep Coercion = Coercion
instance Representational IntMap  where rep Coercion = Coercion

-- * mtl, transformers
instance Representational ReaderT     where rep Coercion = Coercion
instance Representational (ReaderT r) where rep Coercion = Coercion
instance (Representational m)
      => Representational (ReaderT r m) where
  rep c = coercionOf runReaderT >>> rep (rep c) >>> coercionOf ReaderT

instance (Representational m)
      => Representational (ExceptT e m) where
  rep c = coercionOf runExceptT >>> rep (rep c) >>> coercionOf ExceptT

instance Representational (ContT r)   where rep Coercion = Coercion
instance Representational (ContT r m) where rep Coercion = Coercion

instance Representational MaybeT where rep Coercion = Coercion
instance (Representational m)
      => Representational (MaybeT m) where
  rep c = coercionOf runMaybeT >>> rep (rep c) >>> coercionOf MaybeT

firstR :: Coercion a b -> Coercion (a,c) (b,c)
firstR = eta . rep

instance (Representational m)
      => Representational (L.WriterT w m) where
  rep c = coercionOf L.runWriterT >>> rep (firstR c) >>> coercionOf L.WriterT

instance (Representational m)
      => Representational (S.WriterT w m) where
  rep c = coercionOf S.runWriterT >>> rep (firstR c) >>> coercionOf S.WriterT

instance (Representational m)
      => Representational (L.StateT s m) where
  rep c = coercionOf L.runStateT >>> rep (rep (firstR c)) >>> coercionOf L.StateT

instance (Representational m)
      => Representational (S.StateT s m) where
  rep c = coercionOf S.runStateT >>> rep (rep (firstR c)) >>> coercionOf S.StateT

instance Representational L.RWST where rep Coercion = Coercion
instance (Representational m)
      => Representational (L.RWST e w s m) where
  --                                     c :: a ~R b
  --                               rep $ c :: (,,) a ~R (,,) b
  --                   eta . eta . rep $ c :: (a,w,s) ~R (b,w,s)
  --             rep . eta . eta . rep $ c :: m (a,w,s) ~R m (b,w,s)
  -- rep . rep . rep . eta . eta . rep $ c :: e -> s -> m (a,w,s) ~R e -> s -> m (b,w,s)
  rep c = coercionOf L.runRWST >>> (rep . rep . rep . eta . eta . rep) c >>> coercionOf L.RWST

instance Representational S.RWST where rep Coercion = Coercion
instance (Representational m)
      => Representational (S.RWST e w s m) where
  rep c = coercionOf S.runRWST >>> (rep . rep . rep . eta . eta . rep) c >>> coercionOf S.RWST

max_stars_repo_name: viercc/roles

max_issues_repo_name: viercc/roles

max_forks_repo_name: viercc/roles


text:
{-# LANGUAGE DeriveGeneric, DeriveAnyClass, FlexibleContexts, FlexibleInstances, TypeSynonymInstances #-}

module Neural.FF where

import Control.Monad
import Control.Monad.State.Strict
import Control.Parallel.Strategies
import Data.List.Split
import qualified Data.Ix as I
import qualified Data.Vector as V
import GHC.Generics
import Numeric.LinearAlgebra
import System.Random

import Neural.Util
import Numeric.Extra

instance Show (a -> b) where show a = "function"

data FFLayerOptions = FFLayerOptions {
  -- # of slices
  layerDepth :: Int,

  -- z_j -> y_j
  layerActivationFn :: [Matrix Double] -> [Matrix Double],

  -- z_j -> yPrime_j
  layerActivationFnDeriv :: [Matrix Double] -> [Matrix Double],

  -- Pooling functions
  pooling :: (Int,Int) -> [Matrix Double] -> [Matrix Double],
  poolingBP :: (Int,Int) -> [Matrix Double] -> [Matrix Double] -> [Matrix Double],
  poolSize :: (Int,Int)
  }

data FFLayer = FFLayer {
  weights :: [Matrix Double],
  biases :: [Matrix Double],

  -- x_j -> (w_j, b_j) -> z_j
  zFn :: [Matrix Double] -> ([Matrix Double], [Matrix Double]) -> [Matrix Double],
  yFn :: [Matrix Double] -> [Matrix Double],

  -- nabla_b_j -> x_j -> nabla_w_j
  nablaW :: [Matrix Double] -> [Matrix Double] -> [Matrix Double],

  -- delta_j -> x_j -> nabla_b_j
  nablaB :: [Matrix Double] -> [Matrix Double] -> [Matrix Double],

  -- nabla_b_j -> z_i-> w_j -> delta_i
  backprop :: [Matrix Double] -> [Matrix Double] -> [Matrix Double] -> [Matrix Double]
    } deriving (Show, Generic, NFData)

data FFNetwork = FFNetwork {
  layers :: V.Vector(FFLayer)
  } deriving (Show, Generic, NFData)

-- descs is a list of [(spawner, dims, options)]
spawnFFNetwork :: [((Int, Int) -> (Int, Int) -> FFLayerOptions -> FFLayerOptions -> IO (FFLayer),
                   (Int, Int), (Int, Int), FFLayerOptions)] -> IO (FFNetwork)
spawnFFNetwork descs = liftM (FFNetwork . V.fromList)
                       $ zipWithM (\(ls_j, dims_i, dims_j, opts_j) opts_i ->
                                    ls_j dims_i dims_j opts_i opts_j)
                       (tail descs) (map (\(_,_,_,opts_) -> opts_) descs)

-- Parallel matrix multiplication; left operand is split into rows.
(<||#>) :: (Matrix Double) -> (Matrix Double) -> (Matrix Double)
x <||#> y = fromRows $ parMap rdeepseq (<# y) $ toRows x

-- Same as above, but right operand is split into columns
(<#||>) :: (Matrix Double) -> (Matrix Double) -> (Matrix Double)
x <#||> y = fromColumns $ parMap rdeepseq (x #>) $ toColumns y

-- Parallel dense matrix/vector product.
-- (m><n) #||> (n) -> (m)
(#||>) :: (Matrix Double) -> (Vector Double) -> (Vector Double)
x #||> y = fromList $ parMap rdeepseq (<.> y) $ toRows x

-- Idk what this product is called but is sure as hell is parallel
-- splits second operand into a list
-- I should really see if this is faster than <#||>
-- (m) <.||> (n) -> (m,n)
(<.||>) :: (Vector Double) -> (Vector Double) -> (Matrix Double)
x <.||> y = fromColumns $ parMap rdeepseq (\y_i -> scale y_i x) (toList y)

-- feeds forward, collecting layer activations.
-- note y, z are reversed.
-- net -> x_1 -> ([y_n, .. y_1, x_1], [z_n, .. z_1])
ff :: FFNetwork -> [Matrix Double] -> (V.Vector([Matrix Double]), V.Vector([Matrix Double]))
ff net x = V.foldl activator (V.singleton x, V.empty) $ layers net
  where activator (y, z) (FFLayer w b zf yf _ _ _) =
          (\(y_j, z_j) -> (y_j `V.cons` y, z_j `V.cons` z)) $ (\a -> (yf a, a)) $ zf (V.head y) (w, b)

-- ff that only returns final output of the net.
-- net -> x_1 -> y_n
ffOut :: FFNetwork -> [Matrix Double] -> [Matrix Double]
ffOut net x = V.head $ fst $ ff net x

-- feeds forward, then passes error backwards
bp :: FFNetwork -> [Matrix Double] -> [Matrix Double] -> [([Matrix Double], [Matrix Double])]
bp net expected x = V.foldl passBack [(nabla_w_n, nabla_b_n)] kit
  where passBack ((nabla_w_k, nabla_b_k) : nablas) (x_j, z_j, layer_k, layer_j) =
          (nabla_w_j, nabla_b_j) : (nabla_w_k, nabla_b_k) : nablas
          where nabla_w_j = (nablaW layer_j) nabla_b_j x_j
                nabla_b_j = (nablaB layer_j) delta_j z_j
                delta_j = (backprop layer_k) nabla_b_k z_j (weights layer_k)

        -- passBack sees these args:
        -- y[n-1 .. 1], z[n-1 .. 1], layers[n .. 2], layers[n-1 .. 1]
        kit = V.zip4 (V.tail $ V.tail y) (V.tail z) (V.reverse $ V.tail $ ls) (V.reverse $ V.init ls)
        nabla_w_n = (nablaW $ V.last $ ls) nabla_b_n (y V.! 1)
        nabla_b_n = zipWith (-) (y V.! 0) expected
        ls = layers net
        (y, z) = ff net x

-- Error/cost functions.

{-
crossEntropyDelta :: [Matrix Double] -> [Matrix Double] -> [Matrix Double]
crossEntropyDelta = zipWith (-)

crossEntropyError :: [Matrix Double] -> [Matrix Double] -> [Double]
crossEntropyError x y = map (\(x,y) -> sumElements $ - ((y * (cmap log x)) + (1.0 - y) * (cmap log $ 1.0 - x)))
                        $ zip x y
-}

-- Convert an output value to an actual output of the net
-- i.e. '5' -> [0,0,0,0,0,1,0,...]
buildOutput :: Int -> Int -> [Matrix Double]
buildOutput i len = [(len><1) (map (\x -> if x == i then 1.0 else 0.0) [0..len - 1])]

-- Convert a net's output to an integer output value
-- uses max value of the net.
-- i.e. [0,0.5,0.25,-0.1,0.999999,...] -> '4'
convertOutput :: [Matrix Double] -> Int
convertOutput output = (\(i,j) -> i) $ maxIndex $ head output

-- For a given test set, return the number the net gets correct
evaluate :: FFNetwork -> [((Matrix Double), Int)] -> Double -> Int
evaluate net testSet reg =
  sum $ parMap rdeepseq
  (\(input, expected) ->
    let output = ffOut net [input]
    in if convertOutput output == expected then 1 else 0) testSet

-- perform stochastic gradient descent on the net.
-- each sample data is a tuple of (input, expected output) where the output is an integer representing the desired result.
-- maxClassIndex is the max classification ID any output can produce
-- eta is the learning rate.
-- reg is the regularization parameter (for L2 reg)
sgd :: FFNetwork -> StdGen -> [((Matrix Double), Int)] -> Int -> Int -> Double -> Double -> ((FFNetwork, [((Matrix Double), Int)]), StdGen)
sgd net gen trainingSet maxClassIndex batchSize eta reg = ((newNet, newTrainingSet), newGen)
  where newNet = foldl applyBatch net $ chunksOf batchSize newTrainingSet
        applyBatch net batch =
          net { layers =
                  V.zipWith (\layer (nabla_w, nabla_b) ->
                             layer {
                                -- apply the changes to each layer of the net.
                                weights = zipWith adjustWeights (weights layer) nabla_w,
                                biases = zipWith adjustBiases (biases layer) nabla_b
                                }) (layers net)
                  -- backpropagate each member of the batch, summing the resulting change in the net.
                  -- this will yield [(nw, nb) ...] for each layer of the net.
                  $ V.fromList
                  $ (\(batch1 : batches) ->
                      foldl (zipWith (\(nw,nb) (dnw,dnb) -> (zipWith (+) nw dnw, zipWith (+) nb dnb)))
                      batch1 batches)
                  $ parMap rdeepseq
                  (\(input, expected) ->
                    bp net (buildOutput expected $ maxClassIndex + 1) [input]) batch
              }

        -- hmatrix gets a little confused when (Num a, Matrix b) => a * b happens
        -- so use `scale` for now.
        adjustWeights w nabla_w = (scale (1.0 - eta * reg / (fromIntegral $ length trainingSet)) w)
                                  - (scale (eta / (fromIntegral batchSize)) nabla_w)
        adjustBiases b nabla_b = b - (scale (eta / (fromIntegral batchSize)) nabla_b)

        -- cost = 0.5 * reg * l2_norm_sq * (fromIntegral batchSize) / (fromIntegral $ length trainingSet)
        -- l2_norm_sq = V.sum $ V.map (\l -> sum $ parMap rdeepseq (\w -> sumElements $ w * w) $ weights l) $ layers net

        (newTrainingSet, newGen) = shuffle' trainingSet gen


text:
{-# LANGUAGE TemplateHaskell, QuasiQuotes #-}

module Main where

import Language.Rust.Quote    ( ty )
import Language.Rust.Inline
import Foreign.Marshal.Unsafe ( unsafeLocalState )
import Foreign.Ptr            ( Ptr )
import Foreign.Storable       ( peek )
import Foreign.Marshal        ( with )

import Data.Complex           ( Complex(..) )

-- We add to our context the knowledge that a Rust pointer to a 'Complex' is
-- the same as a Haskell 'Ptr' to a 'Complex Float'.
extendContext basic
extendContext (singleton [ty| *mut Complex |] [t| Ptr (Complex Float) |])

setCrateRoot []

-- Rust representation of Haskell's 'Complex'
[rust|
#[repr(C)]
pub struct Complex {
  re: f32,
  im: f32,
}
|]

-- | Add two complex numbers in Rust
add :: Complex Float -> Complex Float -> Complex Float
add z1 z2 = unsafeLocalState $
  with z1 $ \ptr1 ->
    with z2 $ \ptr2 -> do
      ptr3 <- [rustIO| *mut Complex {
        let z1: *mut Complex = $(ptr1: *mut Complex);
        let z2: *mut Complex = $(ptr2: *mut Complex);

        let z3: Box<Complex> = Box::new(unsafe {
          Complex {
            re: (*z1).re + (*z2).re,
            im: (*z1).im + (*z2).im,
          }
        });
        Box::into_raw(z3)
      } |]
      z3 <- peek ptr3
      [rustIO| () {
        unsafe { Box::from_raw($(ptr3: *mut Complex)) };
      } |]
      pure z3

main = let z1 = 1.3 :+ 2.4
           z2 = 7.1 :+ 3.9
           z3 = add z1 z2
       in putStrLn $ unwords [ show z1, "+", show z2, "=", show z3 ]

max_stars_repo_name: hexresearch/inline-rust

max_issues_repo_name: hexresearch/inline-rust

max_forks_repo_name: hexresearch/inline-rust


text:
module Math.IRT.Model.ThreePLM
  ( ThreePLM (..)
  ) where

import Statistics.Distribution

import Math.IRT.Internal.Distribution
import Math.IRT.Internal.LogLikelihood
import Math.IRT.Model.FourPLM ( FourPLM(..) )
import Math.IRT.Model.Generic

data ThreePLM = ThreePLM { discrimination :: !Double
                         , difficulty     :: !Double
                         , pseudoGuessing :: !Double
                         } deriving (Show)

instance Distribution ThreePLM where
    cumulative = cumulative . toFourPLM

instance ContDistr ThreePLM where
    density    = density . toFourPLM
    quantile _ = error "This shouldn't be needed"

instance DensityDeriv ThreePLM where
    densityDeriv = densityDeriv . toFourPLM

instance GenericModel ThreePLM where
    fromRasch         b = ThreePLM 1.0 b 0.0
    fromOnePLM        b = ThreePLM 1.7 b 0.0
    fromTwoPLM      a b = ThreePLM   a b 0.0
    fromThreePLM        = ThreePLM
    fromFourPLM a b c _ = ThreePLM   a b   c

instance LogLikelihood ThreePLM where
    logLikelihood b = logLikelihood b . toFourPLM

toFourPLM :: ThreePLM -> FourPLM
toFourPLM (ThreePLM sa sb sc) = FourPLM sa sb sc 1.0

max_stars_repo_name: argiopetech/irt

max_issues_repo_name: argiopetech/irt

max_forks_repo_name: argiopetech/irt


text:
{-# LANGUAGE TypeSynonymInstances #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}

-- | This module defines the Syntax of the Quantum IO Monad, which is
-- an embedded language for writing quantum computations. It is an
-- alternative definition using the approach of defining F-Algebras.
module QIO.QioSynAlt where

import Data.Monoid as Monoid
import Data.Complex
import Control.Applicative (Applicative(..))
import Control.Monad       (liftM, ap)

-- | For Real numbers, we simply use the built in Double type
type RR = Double

-- | For Complex numbers, we use the built in Complex numbers, over our Real
-- number type (i.e. Double)
type CC = Complex RR

-- | The amplitude of a complex number is the magnitude squared.
amp :: CC -> RR
amp k = (magnitude k)*(magnitude k)

-- | The type of Qubits in QIO are simply integer references.
newtype Qbit = Qbit Int deriving (Num, Enum, Eq, Ord)

-- | We can display a qubit reference
instance Show Qbit where
    show (Qbit q) = "(Qbit:" ++ show q ++ ")"

-- | A rotation is in essence a two-by-two complex valued matrix
type Rotation = ((Bool,Bool) -> CC)

-- | We can display the matrix representation of a rotation
instance Show Rotation where
    show f = "(" ++ (show (f (False,False))) ++ "," ++ (show (f (False,True))) ++ "," ++ (show (f (True,False))) ++ "," ++ (show (f (True,True))) ++ ")"

-- | The identity rotation
rid :: Rotation
rid (x,y) = if x==y then 1 else 0

-- | The not rotation
rnot :: Rotation
rnot (x,y) = if x==y then 0 else 1

-- | The hadamard rotation
rhad :: Rotation
rhad (x,y) = if x && y then -h else h where h = (1/sqrt 2)

-- | The phase rotation
rphase :: RR -> Rotation
rphase _ (False,False)  = 1
rphase r (True,True)    = exp(0:+r)
rphase _ (_,_)          = 0

-- | Returns the inverse (or reverse) of the given rotation
rrev :: Rotation -> Rotation
rrev r (False,True)   = conjugate (r (True,False))
rrev r (True,False)   = conjugate (r (False,True))
rrev r xy             = conjugate (r xy)

-- | Rotations can be compared for equality.
-- They are equal if the define the same matrix.
instance Eq Rotation where
    f == g =    (f (False,False)  == g (False,False))
             && (f (False,True)   == g (False,True))
             && (f (True,False)   == g (True,False))
             && (f (True,True)    == g (True,True))
    f /= g =    (f (False,False)  /= g (False,False))
             || (f (False,True)   /= g (False,True))
             || (f (True,False)   /= g (True,False))
             || (f (True,True)    /= g (True,True))

-- | The non-recursive data type definition of a unitary operation
data UFunctor u = UReturn
                | Rot Qbit Rotation u
                | Swap Qbit Qbit u
                | Cond Qbit (Bool -> u) u
                | Ulet Bool (Qbit -> u) u

-- | In order to define an F-Algebra, 'UFunctor' must be a functor.
instance Functor UFunctor where
    fmap eval UReturn = UReturn
    fmap eval (Rot q r u) = Rot q r (eval u)
    fmap eval (Swap q1 q2 u) = Swap q1 q2 (eval u)
    fmap eval (Cond q f u) = Cond q (eval . f) (eval u)
    fmap eval (Ulet b f u) = Ulet b (eval . f) (eval u)

-- | The fix point type construtor.
newtype Fix f = Fx (f (Fix f))

-- | We can define the inverse of Fx
unFix :: Fix f -> f (Fix f)
unFix (Fx x) = x

-- | We fix the non-recursice data-type in order to get our type 'U'
-- of unitary operations.
type U = Fix UFunctor

-- | The type of an F-Algebra.
type Algebra f a = f a -> a

-- | The type of the initial algebra for UFunctor
type UInitialAlgebra = Algebra UFunctor U

-- | We can now define the initial algebra for U
uInitialAlgebra :: UInitialAlgebra
uInitialAlgebra = Fx

-- | We can use a catamorphism to abstract evaluation over a given
-- algebra
cata :: Functor f => Algebra f a -> Fix f -> a
cata algebra = algebra . fmap (cata algebra) . unFix

-- | The type "U" forms a Monoid.
instance Semigroup U where
    (Fx UReturn) <> u = u
    (Fx (Rot x a u)) <> u' = Fx $ Rot x a (u <> u')
    (Fx (Swap x y u)) <> u' = Fx $ Swap x y (u <> u')
    (Fx (Cond x br u')) <> u'' = Fx $ Cond x br (u' <> u'')
    (Fx (Ulet b f u)) <> u' = Fx $ Ulet b f (u <> u')

instance Monoid U where
    mempty = Fx UReturn

-- | Apply the given rotation to the given qubit
rot :: Qbit -> Rotation -> U
rot x r = Fx $ Rot x r mempty

-- | Swap the state of the two given qubits
swap :: Qbit -> Qbit -> U
swap x y = Fx $ Swap x y mempty

-- | Apply the conditional unitary, depending on the value of the given qubit
cond :: Qbit -> (Bool -> U) -> U
cond x br = Fx $ Cond x br mempty

-- | Introduce an Ancilla qubit in the given state, for use in the sub-unitary
ulet :: Bool -> (Qbit -> U) -> U
ulet b ux = Fx $ Ulet b ux mempty

-- | Apply a not rotation to the given qubit
unot :: Qbit -> U
unot x = rot x rnot

-- | Apply a hadamard rotation to the given qubit
uhad :: Qbit -> U
uhad x = rot x rhad

-- | Apply a phase rotation (of the given angle) to the given qubit
uphase :: Qbit -> RR -> U
uphase x r = rot x (rphase r)

-- | Returns the inverse (or reverse) of the given unitary operation,
-- using an F-Algebra
urev :: U -> U
urev = cata urev_algebra
 where
   urev_algebra :: UFunctor U -> U
   urev_algebra UReturn = Fx UReturn
   urev_algebra (Rot x r u) = u `mappend` rot x (rrev r)
   urev_algebra (Swap x y u) = u `mappend` swap x y
   urev_algebra (Cond x br u) = u `mappend` cond x br
   urev_algebra (Ulet b xu u) = u `mappend` ulet b xu

-- | We can display a representation of a unitary, using an F-Algebra
instance Show U where
  show = cata showU_algebra
    where
      showU_algebra :: UFunctor String -> String
      showU_algebra UReturn = ""
      showU_algebra (Rot q r u) =
       "Rotate " ++ show q ++ ":" ++ show r ++ "\n" ++ u
      showU_algebra (Swap q1 q2 u) =
       "Swap " ++ show q1 ++ " and " ++ show q2 ++ "\n" ++ u
      showU_algebra (Cond q br u) =
       "Cond " ++ show q ++ " \\b -> if b then (\n"
       ++ unlines (map (' ':) (lines $ br True))
       ++ ") else (\n"
       ++ unlines (map (' ':) (lines $ br False))
       ++ ")\n" ++ u
      showU_algebra (Ulet b xu u) =
       let fv = find_fv xu in
       "Ulet " ++ show fv ++ " = " ++ (if b then "1" else "0") ++ " in (\n"
       ++ unlines (map (' ':) (lines $ xu fv))
       ++ ")\n" ++ u
      -- this is currently a dummy function
      find_fv :: (Qbit -> String) -> Qbit
      find_fv _ = -1

-- | The non-recursive data type definition of a QIO computation
data QIOFunctor a q = QReturn a
                    | MkQbit Bool (Qbit -> q)
                    | ApplyU U q
                    | Meas Qbit (Bool -> q)

-- | In order to define an F-Algebra, 'UF' must be a functor.
instance Functor (QIOFunctor a) where
    fmap eval (QReturn a) = QReturn a
    fmap eval (MkQbit b f) = MkQbit b (eval . f)
    fmap eval (ApplyU u q) = ApplyU u (eval q)
    fmap eval (Meas q f) = Meas q (eval . f)

-- | We fix the non-recursice data-type in order to get our type 'U'
-- of unitary operations.
type QIOprim a = Fix (QIOFunctor a)

-- | The type of the initial algebra for UFunctor
type QIOInitialAlgebra a = Algebra (QIOFunctor a) (QIOprim a)

-- | We can now define the initial algebra for U
qioInitialAlgebra :: QIOInitialAlgebra a
qioInitialAlgebra = Fx

-- | The "QIO" type forms a Monad, by wrapping 'QIOprim'
data QIO a = Apply (Fix (QIOFunctor a))

-- | We can remove the wrapper.
primQIO :: QIO a -> QIOprim a
primQIO (Apply q) = q

instance Functor QIO where
    fmap = liftM

instance Applicative QIO  where
    pure  = Apply . Fx . QReturn
    (<*>) = ap

-- | The wrapper type 'ApplyFix' forms a Monad
instance Monad QIO where
    return = pure
    (Apply (Fx (QReturn a))) >>= f = f a
    (Apply (Fx (MkQbit b g))) >>= f = Apply . Fx $
      MkQbit b  (\q -> primQIO $ (Apply (g q)) >>= f)
    (Apply (Fx (ApplyU u q))) >>= f = Apply . Fx $
      ApplyU u $ primQIO (Apply q >>= f)
    (Apply (Fx (Meas x g))) >>= f = Apply . Fx $
      Meas x (\b -> primQIO $ (Apply (g b)) >>= f)

-- | Initialise a qubit in the given state (adding it to the overall quantum state)
mkQbit :: Bool -> QIO Qbit
mkQbit b = Apply . Fx $ MkQbit b (\q -> primQIO (return q))

-- | Apply the given unitary operation to the current quantum state
applyU :: U -> QIO ()
applyU u =  Apply . Fx $ ApplyU u $ primQIO (return ())

-- | Measure the given qubit, and return the measurement outcome (note that this
-- operation may affect the overall quantum state, as a measurement is destructive)
measQbit :: Qbit -> QIO Bool
measQbit x = Apply . Fx $ Meas x (\b -> primQIO (return b))

-- | We can show a QIO computation, using an F-Algebra
instance (Show a) => Show (QIO a) where
  show = (cata showQIO_algebra) . primQIO
   where
    showQIO_algebra :: (Show a) => Algebra (QIOFunctor a) String
    showQIO_algebra (QReturn a) =
       "Return: " ++ show a ++ "\n"
    showQIO_algebra (MkQbit b f) =
       "Init" ++ (if b then "1" else "0") ++ "\n"
       ++ f 0
    showQIO_algebra (ApplyU u qio) =
       "Apply Unitary: (\n"
       ++ unlines (map (' ':) (lines $ show u))
       ++ ")\n" ++ qio
    showQIO_algebra (Meas q f) =
       "Measure " ++ show q ++ " \\b -> if b then (\n"
       ++ unlines (map (' ':) (lines $ f True))
       ++ ") else (\n"
       ++ unlines (map (' ':) (lines $ f False))
       ++ ")\n"

-- | We can count the number of each primitive operation using an F-Algebra
count :: QIO a -> (Int,Int,Int)
count = (cata count_algebra) . primQIO
 where
  count_algebra :: Algebra (QIOFunctor a) (Int,Int,Int)
  count_algebra (QReturn _) = (0,0,0)
  count_algebra (MkQbit b f) = let (mk,ap,ms) = f 0 in
                               (mk+1,ap,ms)
  count_algebra (ApplyU _ (mk,ap,ms)) = (mk,ap+1,ms)
  count_algebra (Meas q f) = let (mk,ap,ms) = f False in
                             (mk,ap,ms+1)

toffoli :: Qbit -> Qbit -> Qbit -> U
toffoli q1 q2 q3 =
  cond q1 (\b1 -> if b1 then (
   cond q2 (\b2 -> if b2 then (unot q3)
    else mempty)) else mempty)

and :: Bool -> Bool -> QIO Bool
and a b = do
 q1 <- mkQbit a
 q2 <- mkQbit b
 q3 <- mkQbit False
 applyU (toffoli q1 q2 q3)
 measQbit q3

max_stars_repo_name: probprob/qio-haskell

max_issues_repo_name: probprob/qio-haskell

max_forks_repo_name: probprob/qio-haskell


text:
{-# LANGUAGE TypeApplications #-}
{-# LANGUAGE FlexibleContexts #-}

module Main where

-- import Prelude hidint (until)
import Control.Monad.State
-- import Control.Monad.Trans.Class
import Numeric.LinearAlgebra.HMatrix hiding (scale, (!))
import qualified Graphics.Rendering.OpenGL as GL
-- import Foreign.Storable
import GLCode
import Sketch
import Points
-- import Lines
-- import Triangles
import Data.Array
-- import Debug.Trace
import Turtle
import Triangles

ex5 :: IO ()
ex5 = do
    mainLoopState 512 512 16 (ident @Float 4) $ \time -> do

        put (ident @Float 4)

    --         let t = 5*realToFrac time :: Float

        GL.clearColor GL.$= GL.Color4 1.0 1.0 1.0 1
        GL.clearAccum GL.$= GL.Color4 1.0 1.0 1.0 1
        io $ GL.clear [GL.ColorBuffer, GL.AccumBuffer]

        GL.lineSmooth GL.$= GL.Enabled
        GL.lineWidth GL.$= 4
        GL.hint GL.LineSmooth GL.$= GL.Nicest
        GL.multisample GL.$= GL.Enabled
        GL.depthFunc GL.$= Just GL.Always

        let p = array ((0, 0), (19, 19)) [((ix, iy), 0.25*cos (-0.75*time+2.0*pi*fromIntegral (ix+iy)/10)+
                                                     0.25*sin (0.5*time+3.0*2*pi*(fromIntegral (ix-iy)/20))) |
                                          ix <- [0..19],
                                          iy <- [0..19]]

        translate (-1.0) (-1.0) 0.0
        drawDensityField (coolwarm (-0.45) 0.45) 20 20 0.1 0.1 p
        let (vx, vy) = grad p
--         let vel = fmap (\(x, y) -> (0.1*x, 0.1*y)) $ interpVelocity vx vy
--         drawVectorField 20 20 0.1 0.1 vx vy

--         lift $ drawPath [(x, y) | i <- [0..10], let x = 0.1*fromIntegral i, let y = x*x]

--         let f t = (t, 0.25*t*t)
--         let time' = 0.3*time
--         let d = 3*(time'-fromIntegral @Int (floor time'))
--         lift $ plotPath 100 f 0 d
--        let curve = integrate (\x y ->let (u,v)=interpVelocityField vx vy x y in (10*u,10*v)) 0 0 100 0.01
--         let curve = integrate (\x y ->let (u,v)=interpVelocityField vx vy x y in (u,v)) 1.0 1.0 300 0.01
        duringAndAfter time 5 10 $ \_ ->
            forM_ [-5.0, -3.0..5.0] $ \x -> do
                forM_ [-5.0, -3.0..5.0] $ \y -> do
                    let curve = integrate (interpVelocityField vx vy) (10+x) (10+y) 500 0.1
                    lift $ drawPath (0.0, 0.0, 0.0) [(0.1*u, 0.1*v) | (u, v) <- curve]

        return ()

image1 = do
    mainLoopState 768 768 16 (ident @Float 4) $ \time -> do
--     mainGifLoopState "euler_or_lagrange.gif" 12 768 768 16 1 100 (ident @Float 4) $ \time -> do

--         let time = fromIntegral 46/12

        put (ident @Float 4)

        GL.clearColor GL.$= GL.Color4 1.0 1.0 1.0 1
        io $ GL.clear [GL.ColorBuffer]

        GL.lineSmooth GL.$= GL.Enabled
        GL.lineWidth GL.$= 3
        GL.hint GL.LineSmooth GL.$= GL.Nicest
        GL.multisample GL.$= GL.Enabled

        translate (-0.95) (-0.90) 0.0
--         let vfield x y = (0.01*cos x+0.01*(y-cos 2*y+0.2*sin(2*x)), -0.02*(sin x-1+0.5*sin (5*y)))
        let vfield x y = (0.01*(y-1)+0.01*(sin y-1), -0.01*(x-1)-0.01*sin (3*x))
        beforeAndDuring time 14 16 $ \time ->
            drawVectorField' (arrow (1-time, 1-time, 1-time) 0.007 0.03 0.03) 13 13 0.15 0.15 $ \i j ->
                let x = 0.15*fromIntegral i
                    y = 0.15*fromIntegral j
                in vfield x y
--         let curve = integrate vfield 0.2 0.1 200 0.2
--         lift $ drawPath (0.0, 0.0, 0.0) (take 100 curve)
--         lift $ drawPath (0.0, 0.0, 0.0) (drop 99 curve)
        duringAndAfter time 5 12 $ \time' -> do
            forM_ [-0.5, -0.4..2.9] $ \starty -> do
                forM_ [-0.5, -0.4..2.9] $ \startx -> do
                    let curve = integrate vfield startx starty (floor (400*time')) 0.2
--                     beforeAndDuring time 14 16 $ \time ->
--                             lift $ drawPath (1-time, 1-time, 1-time) curve
                    when (length curve > 0) $ do
                        lift $ setUniform "turtle_point" "pointSize" (12.0 :: Float)
                        setTransformPoint
                        lift $ drawPoint "turtle_point"
                                   "vPosition" ((\(x, y) -> v2f x y) (last curve))
                                   "color" (v4f 0.0 0.2 0.9 1.0)

main :: IO ()
main = image1

max_stars_repo_name: dpiponi/SketchBox

max_issues_repo_name: dpiponi/SketchBox

max_forks_repo_name: dpiponi/SketchBox


text:
import System.Environment (getArgs)
import Control.Arrow (first)
import Graphics.Rendering.Chart.Easy
import Graphics.Rendering.Chart.Backend.Diagrams
import Data.Time
import qualified Statistics.Sample as S
import qualified Data.Vector as V
import System.Process

import Data.Text.Time (parseISODateTime)
import qualified Data.TimeSeries as TS

signal :: Num a => TS.Series a -> [(LocalTime, a)]
signal ts = map (first (utcToLocalTime utc)) (TS.toList ts)

main = do
    args <- getArgs
    let fn = head args
    ts <- TS.loadCSV TS.HasHeader parseISODateTime fn
    let y3 = TS.years 3
    let xs = TS.rolling y3 (S.mean . V.fromList) ts
    toFile def "dist/smooth.svg" $ do
        layout_title .= fn
        setColors [opaque blue, opaque red]
        plot (line "Orignal" [signal ts])
        plot (line "Smoothed" [signal xs])
    putStrLn "Plot saved to: dist/smooth.svg"
    createProcess (shell "firefox dist/smooth.svg")

max_stars_repo_name: klangner/timeseries

max_issues_repo_name: klangner/timeseries

max_forks_repo_name: klangner/time-series-lib


text:
{-# OPTIONS_HADDOCK hide #-}
-----------------------------------------------------------------------------
-- |
-- Module     : Numeric.LinearAlgebra.Types
-- Copyright  : Copyright (c) 2010, Patrick Perry <patperry@gmail.com>
-- License    : BSD3
-- Maintainer : Patrick Perry <patperry@gmail.com>
-- Stability  : experimental
--
-- Basic type classes and enums.
--

module Numeric.LinearAlgebra.Types (

    -- * Matrix views
    Herm(..),
    Tri(..),

    -- * Matrix factorization views
    Chol(..),

    -- * Vector math types
    VNum,
    VFractional,
    VFloating,

    -- * BLAS element types
    BLAS1,
    BLAS2,
    BLAS3,

    -- * LAPACK element types
    LAPACK,

    -- * Enums
    Trans(..),
    Uplo(..),
    Side(..),
    Diag(..),

    -- * Algorithm-specific parameters
    CovMethod(..),

    -- * Re-export of Complex from Data.Complex
    module Data.Complex,

    -- * Re-export of Storable from Foreign.Storable
    module Foreign.Storable,

    ) where

import Foreign.VMath( VNum, VFractional, VFloating )
import Foreign.BLAS( Trans(..), Uplo(..), Side(..), Diag(..), BLAS1, BLAS2, BLAS3 )
import Foreign.LAPACK( LAPACK )
import Data.Complex( Complex(..) )
import Foreign.Storable( Storable() )

-- | A hermitian view of an underlying matrix.  The view can either be
-- of the upper or lower triangular part of the matrix.  The type arguments
-- are as follows:
--
--     * @m@: the underlyting matrix type.
--
--     * @e@: the element type of the matrix.
--
data Herm m e = Herm Uplo (m e) deriving (Show)

-- | A triangular view of an underlying matrix.  The view can either be
-- upper or lower triangular, with a unit or non-unit diagonal.  The type
-- arguments are as follows:
--
--     * @m@: the underlyting matrix type.
--
--     * @e@: the element type of the matrix.
--
data Tri m e = Tri Uplo Diag (m e) deriving (Show)

-- | A Cholesky decomposition view of a matrix.
data Chol m e = Chol Uplo (m e) deriving (Show)

-- | The method of scaling the sample covariance matrix.
data CovMethod =
      UnbiasedCov -- ^ This is the default behavior. Corresponds to a
                  -- scaling of @n/(n-1)@ in the unweighed case, and
                  -- @1/(1 - \\sum w_i^2)@ in the weighted case, where @w_i@
                  -- is the normalized weight. Note the unweighted and
                  -- weighted cases agree when @w_i = 1/n@.

    | MLCov       -- ^ Returns the centered second moment matrix without
                  -- scaling the result.
    deriving (Eq, Show)

max_stars_repo_name: patperry/hs-linear-algebra

max_issues_repo_name: cartazio/hs-cblas

max_forks_repo_name: cartazio/hs-cblas


text:
module Lib.Practice.P3 where

import Control.Monad
import Numeric.LinearAlgebra

f3 = do
    let gates = [("AND", gAnd), ("NAND", gNand), ("OR", gOr)]
    mapM_ showOutput gates

gAnd = Neuron { weight = vector [0.5, 0.5], bias = -0.7 }
gNand = Neuron { weight = vector [-0.5, -0.5], bias = 0.7 }
gOr = Neuron { weight = vector [0.5, 0.5], bias = -0.2 }

inputs = replicateM 2 [NotFire, Fire]
showOutput (name, gate) = do
    putStrLn $ name ++ "ゲート"
    mapM_ (showOutput' gate) inputs
    putStrLn ""
    where
        showOutput' gate input = do
            let result = input --> gate
            putStrLn $ show input ++ " -> " ++ show result

data Neuron = Neuron {
    weight :: Vector R,
    bias :: R
} deriving Show

data Input = Fire | NotFire deriving Show
input2Int Fire = 1
input2Int NotFire = 0

(-->) = output
output :: [Input] -> Neuron -> Input
output inputs Neuron { weight = w, bias = b } = result
    where inputV    = fromList $ map input2Int inputs
          sum'      = sumElements $ inputV * w
          result'   = b + sum'
          result    = case () of
            _ | result' >= 0    -> Fire
            otherwise           -> NotFire


text:
module Section_2_3_3 where

import Numeric.LinearAlgebra

and' :: Double -> Double -> Double
and' x1 x2 =
  let x = vector [x1, x2]
      w = vector [0.5, 0.5]
      b = (-0.7)
      tmp = w <.> x + b
   in if tmp <= 0
         then 0
         else 1

nand :: Double -> Double -> Double
nand x1 x2 =
  let x = vector [x1, x2]
      w = vector [(-0.5), (-0.5)]
      b = 0.7
      tmp = w <.> x + b
   in if tmp <= 0
         then 0
         else 1

or' :: Double -> Double -> Double
or' x1 x2 =
  let x = vector [x1, x2]
      w = vector [0.5, 0.5]
      b = (-0.2)
      tmp = w <.> x + b
   in if tmp <= 0
         then 0
         else 1


text:
{-# LANGUAGE Arrows #-}

-- | Compute the value of pi by generating random number. Using o Storable vector thanks to Dunai MSF that can handle IO.
module Main where

import Data.Bool (bool)
import Data.Complex
import Data.Vector.Storable.Mutable
import Data.Word
import Dunai.Gloss
import Prelude hiding (read, replicate)

main :: IO ()
main = do
  bitmap <- replicate (s * s * 4) 255
  field <- generate (s * s) init
  field2 <- replicate (s * s) 0
  playDunaiIO (InWindow "Random Pi" (s, s) (100, 100)) white 30 (network field field2 bitmap)
  where
    init :: Int -> Complex Double
    init i =
      let x' = mod i s
          x = fromIntegral x' * dx
          y' = div i s
          y = fromIntegral y' * dx
          x0 = l / 5
          y0 = l / 2
          kx = l
          ky = 0
          r = l / 20
       in mkPolar (exp (- ((x - x0) ^ 2 + (y - y0) ^ 2) / r)) (kx * x + ky * y)

s :: Int
s = 200

l, dx, dt :: Double
l = 10
dx = l / fromIntegral s
dt = 3e-3

ii :: Complex Double
ii = 0 :+ 1

toRGBA ptr = bitmapOfForeignPtr s s (BitmapFormat TopToBottom PxRGBA) ptr False

network field field2 bitmap = proc (upd, _) -> do
  if upd > 0
    then do
      i <- count -< ()
      let (cur, nex) = if odd i then (field, field2) else (field2, field)
      arrM (uncurry update) -< (cur, nex)
      arrM (`render` bitmap) -< nex
    else returnA -< ()
  returnA -< img
  where
    img = toRGBA . fst $ unsafeToForeignPtr0 bitmap

update field nfield = go 0
  where
    dir x y
      | x < 0 = dir 0 y
      | x >= s = dir (s -1) y
      | y < 0 = dir x 0
      | y >= s = dir x (s -1)
      | otherwise = x + s * y
    go i
      | i >= s * s = return ()
      | otherwise = do
        let x = mod i s
            y = div i s
        c <- read field i
        cu <- read field (dir x (y -1))
        cd <- read field (dir x (y + 1))
        cl <- read field (dir (x -1) y)
        cr <- read field (dir (x + 1) y)
        let v = c + (dt :+ 0) * hv x y c cu cd cl cr
        write nfield i v
        go (i + 1)

hv :: Int -> Int -> Complex Double -> Complex Double -> Complex Double -> Complex Double -> Complex Double -> Complex Double
hv x y c cu cd cl cr = - ii * (- a * (cu + cd + cl + cr -4 * c))
  where
    a = 1 / dx * dx :+ 0

render field bitmap = go 0
  where
    go i
      | i >= s * s = return ()
      | otherwise = do
        c <- read field i
        let (re :+ im) = 255 * (c + (1.5 :+ 1.5)) / 3
            r = floor re
            g = 0
            b = floor im
        write bitmap (i * 4) r
        write bitmap (i * 4 + 1) g
        write bitmap (i * 4 + 2) b
        go (i + 1)


text:
{-# LANGUAGE DeriveFunctor             #-}
{-# LANGUAGE ExistentialQuantification #-}
{-# LANGUAGE OverloadedStrings         #-}
module ML.Gym.Value where

import           ML.Gym.DType
import           ML.Gym.Util

import           Control.Applicative
import qualified Control.Exception        as E
import           Control.Monad            (join, void, (>=>))
import qualified CPython                  as Py
import qualified CPython.Constants        as Py
import qualified CPython.Protocols.Object as Py
import qualified CPython.Types            as Py
import qualified CPython.Types.Exception  as Py
import qualified CPython.Types.Module     as Py
import qualified Data.ByteString          as BS
import           Data.Complex             (Complex)
import           Data.Foldable            (toList)
import           Data.Maybe               (fromMaybe)
import qualified Data.Text                as T
import           System.Environment       (getArgs, getProgName)
import           System.Exit
import           System.IO
import           System.Random            (randomRIO)

import           Debug.Trace

data GymValue a
  = GymScalar a
  | Gym1D [a]
  | Gym2D [[a]]
  | Gym3D [[[a]]]
  | GymTuple [GymValue a]
  deriving (Show, Eq, Functor)

gymValueTo1D :: GymValue a -> [a]
gymValueTo1D (GymScalar x) = [x]
gymValueTo1D (Gym1D xs)    = xs
gymValueTo1D (Gym2D xs)    = concat xs
gymValueTo1D (Gym3D xs)    = concat $ concat xs
gymValueTo1D (GymTuple xs) = concatMap gymValueTo1D xs

gymValueTo2D :: GymValue a -> [[a]]
gymValueTo2D (GymScalar x) = [[x]]
gymValueTo2D (Gym1D xs)    = [xs]
gymValueTo2D (Gym2D xs)    = xs
gymValueTo2D (Gym3D xs)    = concat xs
gymValueTo2D (GymTuple xs) = map gymValueTo1D xs

gymValueTo3D :: GymValue a -> [[[a]]]
gymValueTo3D (GymScalar x) = [[[x]]]
gymValueTo3D (Gym1D xs)    = [[xs]]
gymValueTo3D (Gym2D xs)    = [xs]
gymValueTo3D (Gym3D xs)    = xs
gymValueTo3D (GymTuple xs) = map gymValueTo2D xs

fromBool :: Num a => Bool -> a
fromBool True  = 1
fromBool False = 0

max_stars_repo_name: schnecki/gym-haskell

max_issues_repo_name: schnecki/gym-haskell

max_forks_repo_name: schnecki/gym-haskell


text:
{-# LANGUAGE DataKinds              #-}
{-# LANGUAGE FlexibleContexts       #-}
{-# LANGUAGE ScopedTypeVariables    #-}
{-# LANGUAGE TypeFamilies           #-}

module Numeric.Optimization.TwoPhase.Tableau.Pivot
    ( index
    , allCells
    , pivot
    ) where

import           Data.Finite                        (Finite)
import           GHC.TypeLits                       (KnownNat)
import           Numeric.LinearAlgebra              ((!))
import qualified Numeric.LinearAlgebra as LA        (atIndex)
import           Numeric.LinearAlgebra.Static       (L)
import qualified Numeric.LinearAlgebra.Static as LS

type Index rows cols = (Finite rows, Finite cols)

index :: (KnownNat rows, KnownNat cols)
    => Index rows cols
    -> L rows cols
    -> Double
index (i, j) x =
    LA.atIndex (LS.unwrap x) ix
  where
    ix = (fromIntegral i, fromIntegral j)

allCells :: (KnownNat rows, KnownNat cols)
    => (Double -> Bool)
    -> [Index rows cols]
    -> L rows cols
    -> Bool
allCells f is x =
    all f $ fmap (`index` x) is

pivot :: (KnownNat rows, KnownNat cols)
    => (Finite rows, Finite cols) -> L rows cols -> L rows cols
pivot (x', y') table =
    LS.build (newValue $ LS.unwrap table)
  where
    (x, y)          = (fromIntegral x', fromIntegral y')

    newValue mat i' j'
        | i == x    = mat!i!j / mat!x!y
        | j == y    = 0
        | otherwise = mat!x!y * mat!i!j - mat!x!j * mat!i!y
      where
        (i, j) = (round i', round j')

max_stars_repo_name: alex-mckenna/optimum

max_issues_repo_name: alex-mckenna/optimum

max_forks_repo_name: alex-mckenna/optimum


text:
-----------------------------------------------------------------------------
-- |
-- Module      :  Math.Tensor.Internal.LinearAlgebra
-- Copyright   :  (c) 2019 Tobias Reinhart and Nils Alex
-- License     :  MIT
-- Maintainer  :  tobi.reinhart@fau.de, nils.alex@fau.de
--
-- Gaussian elimination algorithm based on hmatrix.
-----------------------------------------------------------------------------
module Math.Tensor.Internal.LinearAlgebra (
-- * Gaussian Elimination
gaussianST,
gaussian,
-- * Linearly Independent Columns
independentColumns,
independentColumnsMat,
-- * Pivots
pivotsU,
findPivotMax)

where

import Numeric.LinearAlgebra
import Numeric.LinearAlgebra.Data
import Numeric.LinearAlgebra.Devel

import Data.List (maximumBy)

import Control.Monad
import Control.Monad.ST

-- | Returns the pivot columns of an upper triangular matrix.
--
-- @
-- λ let mat = (3 >< 4) [1, 0, 2, -3, 0, 0, 1, 0, 0, 0, 0, 0]
-- λ mat
-- (3><4)
--  [ 1.0, 0.0, 2.0, -3.0
--  , 0.0, 0.0, 1.0,  0.0
--  , 0.0, 0.0, 0.0,  0.0 ]
-- λ pivotsU mat
-- [0,2]
-- @
--

pivotsU :: Matrix Double -> [Int]
pivotsU mat = go (0,0)
  where
    go (i,j)
      = case findPivot mat e (i,j) of
          Nothing       -> []
          Just (i', j') -> j' : go (i'+1, j'+1)
    maxAbs = maximum $ map (maximum . map abs) $ toLists mat
    e = eps * maxAbs

eps :: Double
eps = 1e-12

-- find next pivot in upper triangular matrix

findPivot :: Matrix Double -> Double -> (Int, Int) -> Maybe (Int, Int)
findPivot mat e (i, j)
    | n == j = Nothing
    | m == i = Nothing
    | otherwise = case nonZeros of
                    []          -> if n == j+1
                                   then Nothing
                                   else findPivot mat e (i, j+1)
                    (pi, pj):_  -> Just (pi, pj+j)
    where
        m = rows mat
        n = cols mat
        col = mat ¿ [j]
        nonZeros = filter (\(i', _) -> i' >= i) $ find (not . (< e) . abs) col

-- | Find pivot element below position (i, j) with greatest absolute value in the ST monad.

findPivotMax :: Int -> Int -> Int -> Int -> STMatrix s Double -> ST s (Maybe (Int, Int))
findPivotMax m n i j mat
    | n == j = return Nothing
    | m == i = return Nothing
    | otherwise =
        do
          col      <- mapM (\i' -> do
                                    x <- readMatrix mat i' j
                                    return (i', abs x))
                      [i..m-1]
          let nonZeros = filter (not . (<eps) . abs . snd) col
          let (pi, _) = maximumBy (\(_, x) (_, y) -> x `compare` y) nonZeros
          case nonZeros of
            [] -> if n == j+1
                  then return Nothing
                  else findPivotMax m n i (j+1) mat
            _  -> return $ Just (pi, j)

-- gaussian elimination of sub matrix below position (i, j)

gaussian' :: Int -> Int -> Int -> Int -> STMatrix s Double -> ST s ()
gaussian' m n i j mat = do
    iPivot' <- findPivotMax m n i j mat
    case iPivot' of
        Nothing     -> return ()
        Just (r, p) -> do
                          rowOper (SWAP i r (FromCol j)) mat
                          pv <- readMatrix mat i p
                          mapM_ (reduce pv p) [i+1 .. m-1]
                          gaussian' m n (i+1) (p+1) mat
  where
    reduce pv p r = do
                      rv <- readMatrix mat r p
                      if abs rv < eps
                        then return ()
                        else
                         let frac = -rv / pv
                             op   = AXPY frac i r (FromCol p)
                         in do
                             rowOper op mat
                             mapM_ (\j' -> modifyMatrix mat r j' (\x -> if abs x < eps then 0 else x)) [p..n-1]

-- | Gaussian elimination perfomed in-place in the @'ST'@ monad.

gaussianST :: Int -> Int -> STMatrix s Double -> ST s ()
gaussianST m n = gaussian' m n 0 0

-- | Gaussian elimination as pure function. Involves a copy of the input matrix.
--
-- @
-- λ let mat = (3 >< 4) [1, 1, -2, 0, 0, 2, -6, -4, 3, 0, 3, 1]
-- λ mat
-- (3><4)
--  [ 1.0, 1.0, -2.0,  0.0
--  , 0.0, 2.0, -6.0, -4.0
--  , 3.0, 0.0,  3.0,  1.0 ]
-- λ gaussian mat
-- (3><4)
--  [ 3.0, 0.0,  3.0,                1.0
--  , 0.0, 2.0, -6.0,               -4.0
--  , 0.0, 0.0,  0.0, 1.6666666666666667 ]
-- @
--

gaussian :: Matrix Double -> Matrix Double
gaussian mat = runST $ do
    matST <- thawMatrix mat
    gaussianST m n matST
    freezeMatrix matST
  where
    m = rows mat
    n = cols mat

-- | Returns the indices of a maximal linearly independent subset of the columns
--   in the matrix.
--
-- @
-- λ let mat = (3 >< 4) [1, 1, -2, 0, 0, 2, -6, -4, 3, 0, 3, 1]
-- λ mat
-- (3><4)
--  [ 1.0, 1.0, -2.0,  0.0
--  , 0.0, 2.0, -6.0, -4.0
--  , 3.0, 0.0,  3.0,  1.0 ]
-- λ independentColumns mat
-- [0,1,3]
-- @
--

independentColumns :: Matrix Double -> [Int]
independentColumns mat = pivotsU mat'
    where
        mat' = gaussian mat

-- | Returns a sub matrix containing a maximal linearly independent subset of
--   the columns in the matrix.
--
-- @
-- λ let mat = (3 >< 4) [1, 1, -2, 0, 0, 2, -6, -4, 3, 0, 3, 1]
-- λ mat
-- (3><4)
--  [ 1.0, 1.0, -2.0,  0.0
--  , 0.0, 2.0, -6.0, -4.0
--  , 3.0, 0.0,  3.0,  1.0 ]
-- λ independentColumnsMat mat
-- (3><3)
--  [ 1.0, 1.0,  0.0
--  , 0.0, 2.0, -4.0
--  , 3.0, 0.0,  1.0 ]
-- @
--

independentColumnsMat :: Matrix Double -> Matrix Double
independentColumnsMat mat =
  case independentColumns mat of
    [] -> (rows mat >< 1) $ repeat 0
    cs -> mat ¿ cs