text:
{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}

{-# LANGUAGE CPP                  #-}
{-# LANGUAGE DataKinds            #-}
{-# LANGUAGE FlexibleContexts     #-}
{-# LANGUAGE GADTs                #-}
{-# LANGUAGE RankNTypes           #-}
{-# LANGUAGE ScopedTypeVariables  #-}
{-# LANGUAGE StandaloneDeriving   #-}
{-# LANGUAGE TypeFamilies         #-}
{-# LANGUAGE TypeOperators        #-}
{-# LANGUAGE UndecidableInstances #-}
{-|
Module      : Grenade.Core.Shape
Description : Dependently typed shapes of data which are passed between layers of a network
Copyright   : (c) Huw Campbell, 2016-2017
License     : BSD2
Stability   : experimental

-}
module Grenade.Core.Shape (
    S (..)
  , Shape (..)
#if MIN_VERSION_singletons(2,6,0)
  , SShape (..)
#else
  , Sing (..)
#endif

  , randomOfShape
  , fromStorable
  , fromStorableMatrix
  , nk
  , visualise2D
  , splitChannels
  , combineChannels
  ) where

#if MIN_VERSION_singletons(2,6,0)
import           Data.Kind                    (Type)
#endif

import           Control.DeepSeq              (NFData (..))
import           Data.List.Split              (chunksOf)
import           Data.Maybe                   (fromJust)
import           Data.Proxy
import           Data.Serialize
import           Data.Singletons
import           Data.Singletons.TypeLits
import           Data.Vector.Storable         (Vector)
import qualified Data.Vector.Storable         as V
import           GHC.TypeLits                 hiding (natVal)
import           Numeric.LinearAlgebra        (Matrix)
import qualified Numeric.LinearAlgebra        as NLA
import qualified Numeric.LinearAlgebra.Data   as NLAD
import           Numeric.LinearAlgebra.Static
import qualified Numeric.LinearAlgebra.Static as H

import           System.Random.MWC

import           Grenade.Types

-- | The current shapes we accept.
--   at the moment this is just one, two, and three dimensional
--   Vectors/Matricies.
--
--   These are only used with DataKinds, as Kind `Shape`, with Types 'D1, 'D2, 'D3.
data Shape
  = D1 Nat
  -- ^ One dimensional vector
  | D2 Nat Nat
  -- ^ Two dimensional matrix. Row, Column.
  | D3 Nat Nat Nat
  -- ^ Three dimensional matrix. Row, Column, Channels.
  | D4 Nat Nat Nat Nat
  -- ^ Four dimensional matrix. Depth, Channels, Row, Column

-- | Concrete data structures for a Shape.
--
--   All shapes are held in contiguous memory.
--   3D is held in a matrix (usually row oriented) which has height depth * rows.
data S (n :: Shape) where
  S1D :: ( KnownNat len )
      => R len
      -> S ('D1 len)

  S2D :: ( KnownNat rows, KnownNat columns )
      => L rows columns
      -> S ('D2 rows columns)

  S3D :: ( KnownNat rows
         , KnownNat columns
         , KnownNat channels
         , KnownNat (rows * channels) )
      => L (rows * channels) columns
      -> S ('D3 rows columns channels)

  S4D :: ( KnownNat rows
         , KnownNat columns
         , KnownNat channels
         , KnownNat depth)
      => L (depth * channels * rows) columns
      -> S ('D4 depth channels rows columns)

deriving instance Show (S n)

-- Singleton instances.
--
-- These could probably be derived with template haskell, but this seems
-- clear and makes adding the KnownNat constraints simple.
-- We can also keep our code TH free, which is great.
#if MIN_VERSION_singletons(2,6,0)
-- In singletons 2.6 Sing switched from a data family to a type family.
type instance Sing = SShape

-- | Datatype to allow "pattern matching" on shapes at the type level through
--   the use the sing function from Singletons.
data SShape :: Shape -> Type where
  D1Sing :: Sing a -> SShape ('D1 a)
  D2Sing :: Sing a -> Sing b -> SShape ('D2 a b)
  D3Sing :: KnownNat (a * c) => Sing a -> Sing b -> Sing c -> SShape ('D3 a b c)
  D4Sing :: KnownNat (a * c * d) => Sing a -> Sing b -> Sing c -> Sing d -> SShape ('D4 a b c d)
#else
data instance Sing (n :: Shape) where
  D1Sing :: Sing a -> Sing ('D1 a)
  D2Sing :: Sing a -> Sing b -> Sing ('D2 a b)
  D3Sing :: KnownNat (a * c) => Sing a -> Sing b -> Sing c -> Sing ('D3 a b c)
  D4Sing :: KnownNat (a * c * d) => Sing a -> Sing b -> Sing c -> Sing d -> Sing ('D4 a b c d)
#endif

instance KnownNat a => SingI ('D1 a) where
  sing = D1Sing sing
instance (KnownNat a, KnownNat b) => SingI ('D2 a b) where
  sing = D2Sing sing sing
instance (KnownNat a, KnownNat b, KnownNat c, KnownNat (a * c)) => SingI ('D3 a b c) where
  sing = D3Sing sing sing sing
instance (KnownNat a, KnownNat b, KnownNat c, KnownNat d) => SingI ('D4 a b c d) where
  sing = D4Sing sing sing sing sing

instance SingI x => Num (S x) where
  (+) = n2 (+)
  (-) = n2 (-)
  (*) = n2 (*)
  abs = n1 abs
  signum = n1 signum
  fromInteger x = nk (fromInteger x)

instance SingI x => Fractional (S x) where
  (/) = n2 (/)
  recip = n1 recip
  fromRational x = nk (fromRational x)

instance SingI x => Floating (S x) where
  pi = nk pi
  exp = n1 exp
  log = n1 log
  sqrt = n1 sqrt
  (**) = n2 (**)
  logBase = n2 logBase
  sin = n1 sin
  cos = n1 cos
  tan = n1 tan
  asin = n1 asin
  acos = n1 acos
  atan = n1 atan
  sinh = n1 sinh
  cosh = n1 cosh
  tanh = n1 tanh
  asinh = n1 asinh
  acosh = n1 acosh
  atanh = n1 atanh

--
-- I haven't made shapes strict, as sometimes they're not needed
-- (the last input gradient back for instance)
--
instance NFData (S x) where
  rnf (S1D x) = rnf x
  rnf (S2D x) = rnf x
  rnf (S3D x) = rnf x
  rnf (S4D x) = rnf x

-- | Generate random data of the desired shape
randomOfShape :: forall x . (SingI x) => IO (S x)
randomOfShape = do
  seed :: Int <- withSystemRandom . asGenST $ \gen -> uniform gen
  return $ case (sing :: Sing x) of
    D1Sing SNat ->
        S1D (H.randomVector seed H.Uniform * 2 - 1)

    D2Sing SNat SNat ->
        S2D (H.uniformSample seed (-1) 1)

    D3Sing SNat SNat SNat ->
        S3D (H.uniformSample seed (-1) 1)

    D4Sing SNat SNat SNat SNat ->
        S4D (H.uniformSample seed (-1) 1)

-- | Generate a shape from a Storable Vector.
--
--   Returns Nothing if the vector is of the wrong size.
fromStorable :: forall x. SingI x => Vector RealNum -> Maybe (S x)
fromStorable xs = case sing :: Sing x of
    D1Sing SNat ->
      S1D <$> H.create xs

    D2Sing SNat SNat ->
      S2D <$> mkL xs

    D3Sing SNat SNat SNat ->
      S3D <$> mkL xs

    D4Sing SNat SNat SNat SNat ->
      S4D <$> mkL xs
  where
    mkL :: forall rows columns. (KnownNat rows, KnownNat columns)
        => Vector RealNum -> Maybe (L rows columns)
    mkL v =
      let rows    = fromIntegral $ natVal (Proxy :: Proxy rows)
          columns = fromIntegral $ natVal (Proxy :: Proxy columns)
      in  if rows * columns == V.length v
             then H.create $ NLA.reshape columns v
             else Nothing

-- | Generate a shape from a Storable Matrix.
--
--   Returns Nothing if the matrix is of the wrong size (or one dimensional)
fromStorableMatrix :: forall x. SingI x => Matrix RealNum -> Maybe (S x)
fromStorableMatrix xs = case sing :: Sing x of
    D1Sing SNat ->
      Nothing

    D2Sing SNat SNat ->
      S2D <$> mkL xs

    D3Sing SNat SNat SNat ->
      S3D <$> mkL xs

    D4Sing SNat SNat SNat SNat ->
      S4D <$> mkL xs
  where
    mkL :: forall rows columns. (KnownNat rows, KnownNat columns)
        => Matrix RealNum -> Maybe (L rows columns)
    mkL m =
      let rows    = fromIntegral $ natVal (Proxy :: Proxy rows)
          columns = fromIntegral $ natVal (Proxy :: Proxy columns)
          (h, w)  = NLA.size m
      in  if rows == h && columns == w
             then H.create m
             else Nothing

instance SingI x => Serialize (S x) where
  put i = (case i of
            (S1D x) -> putListOf put . NLA.toList . H.extract $ x
            (S2D x) -> putListOf put . NLA.toList . NLA.flatten . H.extract $ x
            (S3D x) -> putListOf put . NLA.toList . NLA.flatten . H.extract $ x
            (S4D x) -> putListOf put . NLA.toList . NLA.flatten . H.extract $ x
          ) :: PutM ()

  get = do
    Just i <- fromStorable . V.fromList <$> getListOf get
    return i

-- Helper function for creating the number instances
n1 :: ( forall a. Floating a => a -> a ) -> S x -> S x
n1 f (S1D x) = S1D (f x)
n1 f (S2D x) = S2D (f x)
n1 f (S3D x) = S3D (f x)
n1 f (S4D x) = S4D (f x)

-- Helper function for creating the number instances
n2 :: ( forall a. Floating a => a -> a -> a ) -> S x -> S x -> S x
n2 f (S1D x) (S1D y) = S1D (f x y)
n2 f (S2D x) (S2D y) = S2D (f x y)
n2 f (S3D x) (S3D y) = S3D (f x y)
n2 f (S4D x) (S4D y) = S4D (f x y)

-- | Helper function for creating the number instances
nk :: forall x. SingI x => RealNum -> S x
nk x = case (sing :: Sing x) of
  D1Sing SNat ->
    S1D (H.konst x)

  D2Sing SNat SNat ->
    S2D (H.konst x)

  D3Sing SNat SNat SNat ->
    S3D (H.konst x)

  D4Sing SNat SNat SNat SNat ->
    S4D (H.konst x)

-- | Prints out the contents of a matrix with entries between zero and the max value
visualise2D :: S ('D2 a b)  -- ^ input matrix
            -> RealNum      -- ^ maximum element value
            -> String       -- ^ pretty printed matrix
visualise2D (S2D mm) max =
  let m  = H.extract mm
      ms = NLAD.toLists m
      render n' | n' <= 0.2 * max  = ' '
                | n' <= 0.4 * max  = '.'
                | n' <= 0.6 * max  = '-'
                | n' <= 0.8 * max  = '='
                | otherwise =  '#'
      px = (fmap . fmap) render ms
  in unlines px

-- | TODO Theo
splitChannels :: forall rows columns channels.
                 (KnownNat rows, KnownNat columns)
                 => S ('D3 rows columns channels) -> [S ('D2 rows columns)]
splitChannels (S3D x)
  = let r   = fromIntegral $ natVal (Proxy :: Proxy rows)
        rs  = NLA.toRows $ H.extract x
        rs' = chunksOf r rs
        ms  = map (S2D . fromJust . H.create . NLA.fromRows) rs' :: [S ('D2 rows columns)]
    in ms

-- | TODO Theo
combineChannels :: forall rows columns channels.
                 (KnownNat rows, KnownNat columns, KnownNat channels)
                 => [S ('D2 rows columns)] -> S ('D3 rows columns channels)
combineChannels xs
  = let xs' = map (\(S2D x) -> NLA.toRows $ H.extract x) xs :: [[Vector RealNum]]
        xs'' = concat xs'
    in  S3D $ fromJust . H.create . NLA.fromRows $ xs''

max_stars_repo_name: th-char/grenade

max_issues_repo_name: th-char/grenade

max_forks_repo_name: th-char/grenade


text:
module Math.IRT.MLE.Fenced
  ( DF (..)
  , MLEResult (..)
  , mleEst
  ) where

import Data.Default.Class

import Statistics.Distribution

import Math.IRT.Internal.Distribution
import Math.IRT.Internal.LogLikelihood
import Math.IRT.MLE.Internal.Generic
import Math.IRT.Model.Generic

data DF = DF { steps                :: !Int
             , thetaEstimate        :: !Double
             , lower_fence          :: !Double
             , upper_fence          :: !Double
             , fence_discrimination :: !Double }

instance Default DF where
  def = DF 10 0.0 (-3.5) 3.5 3.0

mleEst :: (ContDistr d, DensityDeriv d, LogLikelihood d, GenericModel d) => DF -> [Bool] -> [d] -> MLEResult
mleEst (DF n vTheta lf uf fd) rs params =
    let resp = True : False : rs
        pars = fromThreePLM fd lf 0 : fromThreePLM fd uf 0 : params
    in generic_mleEst resp pars n vTheta

ext: hs

max_stars_repo_name: argiopetech/irt

max_issues_repo_name: argiopetech/irt

max_forks_repo_name: argiopetech/irt


text:
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE ScopedTypeVariables #-}
module Main where

-- import Data.Array.Repa (Array, U, DIM2)
import qualified Data.Array.Repa as R
import qualified Data.IntMap as IM
import Control.Monad
import Data.Complex
import qualified Data.Vector as V
import qualified Linear.V as V
import LSI
import LSI.Examples
import LSI.FrequencyResponse
import LSI.Grid
import LSI.RationalFunction
import LSI.TransferFunction
import Unsafe.Coerce (unsafeCoerce)

graph :: SystemGraph 2 String Double
graph = toGraph $ deriche_ydiff 1.0

tfs :: IM.IntMap (RationalFunction 2 Double)
tfs = IM.map simplify $ transferFunctions graph

frequencyResponses :: IM.IntMap (V.V 2 Double -> Complex Double)
frequencyResponses = IM.map computeFR tfs

n :: Num a => a
n = 50

grid :: R.Array R.D R.DIM2 (Double, Double)
grid = frequencyGrid2D n

nPoints :: Num a => Double
nPoints = (2*n)^2

plotFR :: (V.V 2 Double -> Complex Double) -> R.Array R.D R.DIM2 (Complex Double)
plotFR fr = R.map (fr . V.V . V.fromList . \(w1,w2) -> [w1,w2]) grid

l2Norm :: R.Array R.D R.DIM2 (Complex Double) -> IO Double
l2Norm plot = fmap (/ nPoints) $ R.foldAllP (+) 0 squares
  where squares = R.map (\x -> realPart $ (abs x)^2) plot

plots :: IM.IntMap (R.Array R.D R.DIM2 (Complex Double))
plots = IM.map plotFR $ unsafeCoerce frequencyResponses -- TODO: Kind mismatch.

main :: IO ()
main = do
  -- print grid
  -- let plot = plots IM.! 13
  -- forM_ [0..2*n+1] $ \(i::Int) ->
  --   forM_ [0..2*n+1] $ \(j::Int) -> do
  --     let (x,y,z) = (i,j, (^ 2) . realPart . abs $ plot R.! (R.Z R.:. i R.:. j))
  --     putStrLn $ show i ++ ", " ++ show j ++ ", " ++ show z
  -- print $ toMatlabFunction $ tfs IM.! 13
  -- norm <- l2Norm plot
  -- print norm

  let ps = IM.elems plots
  norms <- forM ps $ \p -> l2Norm p
  print norms

max_stars_repo_name: gdeest/lsi-systems

max_issues_repo_name: gdeest/lsi-systems

max_forks_repo_name: gdeest/lsi-systems


text:
{-# LANGUAGE CPP                    #-}
{-# LANGUAGE FlexibleInstances      #-}
{-# LANGUAGE FunctionalDependencies #-}
{-# LANGUAGE GADTs                  #-}
{-# LANGUAGE RecordWildCards        #-}
-- |
--
-- @regression-simple@ provides (hopefully) simple regression functions.
--
-- The @'linear' :: Foldable f => (a -> (Double, Double)) -> f a -> 'V2'@
-- is the simplest one.
--
-- There are variants with weights, y-errors, and x and y-errors.
-- In addition, package includes Levenberg–Marquardt algorithm implementation
-- to fit arbitrary functions (with one, two or three parameters),
-- as long as you can give their partial derivatives as well (@ad@ package is handy for that).
--
-- For multiple independent variable ordinary least squares
-- or Levenberg-Marquard with functions with \> 3 parameter you should look elsewhere.
--
-- Package has been tested to return similar results as @fit@ functionality in @gnuplot@
-- (L-M doesn't always converge to exactly the same points in parameter space).
--
module Math.Regression.Simple (
    -- * Linear regression
    linear,
    linearFit,
    linearWithWeights,
    linearWithYerrors,
    linearWithXYerrors,
    -- ** Step-by-step interface
    linearFit',
    LinRegAcc (..),
    zeroLinRegAcc,
    addLinReg,
    addLinRegW,
    -- * Quadratic regression
    quadratic,
    quadraticFit,
    quadraticWithWeights,
    quadraticWithYerrors,
    quadraticWithXYerrors,
    -- ** Step-by-step interface
    quadraticFit',
    QuadRegAcc (..),
    zeroQuadRegAcc,
    addQuadReg,
    addQuadRegW,
    quadRegAccToLin,
    -- * Levenberg–Marquardt algorithm
    -- ** One parameter
    levenbergMarquardt1,
    levenbergMarquardt1WithWeights,
    levenbergMarquardt1WithYerrors,
    levenbergMarquardt1WithXYerrors,
    -- ** Two parameters
    levenbergMarquardt2,
    levenbergMarquardt2WithWeights,
    levenbergMarquardt2WithYerrors,
    levenbergMarquardt2WithXYerrors,
    -- ** Three parameters
    levenbergMarquardt3,
    levenbergMarquardt3WithWeights,
    levenbergMarquardt3WithYerrors,
    levenbergMarquardt3WithXYerrors,
    -- * Auxiliary types
    Fit (..),
    V2 (..),
    V3 (..),
) where

import Control.DeepSeq (NFData (..))

import qualified Data.Foldable      as F
import qualified Data.List.NonEmpty as NE

import Math.Regression.Simple.LinAlg
import Numeric.KBN

-- $setup
-- >>> :set -XDeriveFunctor -XDeriveFoldable -XDeriveTraversable
--
-- >>> import Numeric (showFFloat)
-- >>> import Data.List.NonEmpty (NonEmpty (..))
-- >>> import qualified Data.List.NonEmpty as NE
--
-- Don't show too much decimal digits
--
-- >>> let showDouble x = showFFloat (Just (min 5 (5 - ceiling (logBase 10 x)))) (x :: Double)
-- >>> showDouble 123.456 ""
-- "123.46"
--
-- >>> showDouble 1234567 ""
-- "1234567"
--
-- >>> showDouble 123.4567890123456789 ""
-- "123.46"
--
-- >>> showDouble 0.0000000000000012345 ""
-- "0.00000"
--
-- >>> newtype PP a = PP a
-- >>> class Show' a where showsPrec' :: Int -> a -> ShowS
-- >>> instance Show' a => Show (PP a) where showsPrec d (PP x) = showsPrec' d x
-- >>> instance Show' Double where showsPrec' d x = if x < 0 then showParen (d > 6) (showChar '-' . showDouble (negate x)) else showDouble x
-- >>> instance Show' Int where showsPrec' = showsPrec
-- >>> instance (Show' a, Show' b) => Show' (a, b) where showsPrec' d (x, y) = showParen True $ showsPrec' 0 x . showString ", " . showsPrec' 0 y
-- >>> instance Show' v => Show' (Fit v) where showsPrec' d (Fit p e ndf wssr) = showParen (d > 10) $ showString "Fit " . showsPrec' 11 p . showChar ' ' . showsPrec' 11 e . showChar ' ' . showsPrec' 11 ndf . showChar ' ' . showsPrec' 11 wssr
-- >>> instance Show' V2 where showsPrec' d (V2 x y)   = showParen (d > 10) $ showString "V2 " . showsPrec' 11 x . showChar ' ' . showsPrec' 11 y
-- >>> instance Show' V3 where showsPrec' d (V3 x y z) = showParen (d > 10) $ showString "V3 " . showsPrec' 11 x . showChar ' ' . showsPrec' 11 y . showChar ' ' . showsPrec' 11 z
-- >>> instance Show' a => Show' [a] where showsPrec' _ [] = id; showsPrec' _ [x] = showsPrec' 0 x; showsPrec' _ (x:xs) = showsPrec' 0 x . showChar '\n' . showsPrec' 0 xs
-- >>> instance Show' a => Show' (NonEmpty a) where showsPrec' _ (x :| []) = showsPrec' 0 x; showsPrec' _ (x:|xs) = showsPrec' 0 x . showChar '\n' . showsPrec' 0 xs
-- >>>
--
-- Inputs:
-- >>> let input1 = [(0, 1), (1, 3), (2, 5)]
-- >>> let input2 = [(0.1, 1.2), (1.3, 3.1), (1.9, 4.9), (3.0, 7.1), (4.1, 9.0)]
-- >>> let input3 = [(0, 2), (1, 3), (2, 6), (3, 11)]
--
-- >>> let sq z = z * z
--

-------------------------------------------------------------------------------
-- Linear
-------------------------------------------------------------------------------

-- | Linear regression.
--
-- >>> let input1 = [(0, 1), (1, 3), (2, 5)]
-- >>> PP $ linear id input1
-- V2 2.0000 1.00000
--
-- >>> let input2 = [(0.1, 1.2), (1.3, 3.1), (1.9, 4.9), (3.0, 7.1), (4.1, 9.0)]
-- >>> PP $ linear id input2
-- V2 2.0063 0.88685
--
linear :: F.Foldable f => (a -> (Double, Double)) -> f a -> V2
linear f = fitParams . linearFit f

-- | Like 'linear' but returns complete 'Fit'.
--
-- To get confidence intervals you should multiply the errors
-- by @quantile (studentT (n - 2)) ci'@ from @statistics@ package
-- or similar.
-- For big @n@ using value 1 gives 68% interval and using value 2 gives 95% confidence interval.
-- See https://en.wikipedia.org/wiki/Student%27s_t-distribution#Table_of_selected_values
-- (@quantile@ calculates one-sided values, you need two-sided, thus adjust @ci@ value).
--
-- The first input is perfect fit:
--
-- >>> let fit = linearFit id input1
-- >>> PP fit
-- Fit (V2 2.0000 1.00000) (V2 0.00000 0.00000) 1 0.00000
--
-- The second input is quite good:
--
-- >>> PP $ linearFit id input2
-- Fit (V2 2.0063 0.88685) (V2 0.09550 0.23826) 3 0.25962
--
-- But the third input isn't so much,
-- standard error of a slope parameter is 20%.
--
-- >>> let input3 = [(0, 2), (1, 3), (2, 6), (3, 11)]
-- >>> PP $ linearFit id input3
-- Fit (V2 3.0000 1.00000) (V2 0.63246 1.1832) 2 4.0000
--
linearFit :: F.Foldable f => (a -> (Double, Double)) -> f a -> Fit V2
linearFit f = linearFit' . linRegAcc f

-- | Weighted linear regression.
--
-- >>> let input2 = [(0.1, 1.2), (1.3, 3.1), (1.9, 4.9), (3.0, 7.1), (4.1, 9.0)]
-- >>> PP $ linearFit id input2
-- Fit (V2 2.0063 0.88685) (V2 0.09550 0.23826) 3 0.25962
--
-- >>> let input2w = [(0.1, 1.2, 1), (1.3, 3.1, 1), (1.9, 4.9, 1), (3.0, 7.1, 1/4), (4.1, 9.0, 1/4)]
-- >>> PP $ linearWithWeights id input2w
-- Fit (V2 2.0060 0.86993) (V2 0.12926 0.23696) 3 0.22074
--
linearWithWeights :: F.Foldable f => (a -> (Double, Double, Double)) -> f a -> Fit V2
linearWithWeights f = linearFit' . linRegAccW f

-- | Linear regression with y-errors.
--
-- >>> let input2y = [(0.1, 1.2, 0.12), (1.3, 3.1, 0.31), (1.9, 4.9, 0.49), (3.0, 7.1, 0.71), (4.1, 9.0, 1.9)]
-- >>> let fit = linearWithYerrors id input2y
-- >>> PP fit
-- Fit (V2 1.9104 0.98302) (V2 0.13006 0.10462) 3 2.0930
--
-- When we know actual y-errors, we can calculate the Q-value using @statistics@ package:
--
-- >>> import qualified Statistics.Distribution            as S
-- >>> import qualified Statistics.Distribution.ChiSquared as S
-- >>> S.cumulative (S.chiSquared (fitNDF fit)) (fitWSSR fit)
-- 0.446669639443138
--
-- or using @math-functions@
--
-- >>> import Numeric.SpecFunctions (incompleteGamma)
-- >>> incompleteGamma (fromIntegral (fitNDF fit) / 2) (fitWSSR fit / 2)
-- 0.446669639443138
--
-- It is not uncommon to deem acceptable on equal terms any models with, say, Q > 0.001.
-- If Q is too large, too near to 1 is most likely caused by overestimating
-- the y-errors.
--
--
linearWithYerrors :: F.Foldable f => (a -> (Double, Double, Double)) -> f a -> Fit V2
linearWithYerrors f = linearWithWeights f' where
    f' a = case f a of
        (x, y, dy) -> (x, y, recip (dy * dy))

-- | Iterative linear regression with x and y errors.
--
-- /Orear, J. (1982). Least squares when both variables have uncertainties. American Journal of Physics, 50(10), 912–916. doi:10.1119\/1.12972/
--
-- >>> let input2xy = [(0.1, 1.2, 0.01, 0.12), (1.3, 3.1, 0.13, 0.31), (1.9, 4.9, 0.19, 0.49), (3.0, 7.1, 0.3, 0.71), (4.1, 9.0, 0.41, 1.9)]
-- >>> let fit :| fits = linearWithXYerrors id input2xy
--
-- First fit is done using 'linearWithYerrors':
--
-- >>> PP fit
-- Fit (V2 1.9104 0.98302) (V2 0.13006 0.10462) 3 2.0930
--
-- After that the effective variance is used to refine the fit,
-- just a few iterations is often enough:
--
-- >>> PP $ take 3 fits
-- Fit (V2 1.9092 0.99251) (V2 0.12417 0.08412) 3 1.2992
-- Fit (V2 1.9092 0.99250) (V2 0.12418 0.08414) 3 1.2998
-- Fit (V2 1.9092 0.99250) (V2 0.12418 0.08414) 3 1.2998
--
linearWithXYerrors
    :: F.Foldable f
    => (a -> (Double, Double, Double, Double))  -- ^ \(x_i, y_i, \delta x_i, \delta y_i\)
    -> f a                                      -- ^ data
    -> NE.NonEmpty (Fit V2)
linearWithXYerrors f xs = iterate1 go fit0 where
    fit0   = linearWithYerrors (\a -> case f a of (x,y,_,dy) -> (x,y,dy)) xs
    go fit = linearWithWeights (\a -> case f a of (x,y,dx,dy) -> (x,y,recip $ sq (param1 * dx) + sq dy)) xs where
        V2 param1 _ = fitParams fit

-- >>> import qualified Numeric.AD.Mode.Reverse.Double as AD
-- >>> data H3 a = H3 a a a deriving (Functor, Foldable, Traversable)
-- >>> let linearF (H3 a b x) = a * x + b
-- >>> let lin' (V2 a b) (x, y, dx, dy) = case AD.grad' linearF (H3 a b x) of (f, H3 da db f') -> (y, f, V2 da db, recip $ sq (f' * dx) + sq dy)
--
-- >>> PP $ NE.last $ levenbergMarquardt2WithWeights lin' (V2 1 1) input2xy
-- Fit (V2 1.9092 0.99250) (V2 0.12418 0.08414) 3 1.2998

-- | Calculate linear fit from 'LinRegAcc'.
linearFit' :: LinRegAcc -> Fit V2
linearFit' LinRegAcc {..} = Fit params errors ndf wssr where
    matrix@(SM22 a11 _ a22) = inv (SM22 x2 x w)
    params@(V2 a b)         = mult matrix (V2 xy y)

    errors = V2 sa sb

    -- ensure that error is always non-negative.
    -- Due rounding errors, in perfect fit situations it can be slightly negative.
    wssr = max 0 (y2 - a * xy - b * y)
    ndf  = lra_n - 2
    ndf' = fromIntegral ndf :: Double

    sa   = sqrt (a11 * wssr / ndf')
    sb   = sqrt (a22 * wssr / ndf')

    w  = getKBN lra_w
    x  = getKBN lra_x
    x2 = getKBN lra_x2
    y  = getKBN lra_y
    xy = getKBN lra_xy
    y2 = getKBN lra_y2

    -- is it useful?
    -- r2 = (n * xy - x*y) / (sqrt (n * x2 - x*x) * sqrt (n * y2 - y*y))

linRegAcc :: F.Foldable f => (a -> (Double, Double)) -> f a -> LinRegAcc
linRegAcc f = F.foldl' (\acc a -> case f a of (x,y) -> addLinReg acc x y) zeroLinRegAcc

linRegAccW :: F.Foldable f => (a -> (Double, Double, Double)) -> f a -> LinRegAcc
linRegAccW f = F.foldl' (\acc a -> case f a of (x,y,w) -> addLinRegW acc x y w) zeroLinRegAcc

-------------------------------------------------------------------------------
-- Quadractic
-------------------------------------------------------------------------------

-- | Quadratic regression.
--
-- >>> let input1 = [(0, 1), (1, 3), (2, 5)]
-- >>> quadratic id input1
-- V3 0.0 2.0 1.0
--
-- >>> let input2 = [(0.1, 1.2), (1.3, 3.1), (1.9, 4.9), (3.0, 7.1), (4.1, 9.0)]
-- >>> PP $ quadratic id input2
-- V3 (-0.00589) 2.0313 0.87155
--
-- >>> let input3 = [(0, 2), (1, 3), (2, 6), (3, 11)]
-- >>> PP $ quadratic id input3
-- V3 1.00000 0.00000 2.0000
--
quadratic :: F.Foldable f => (a -> (Double, Double)) -> f a -> V3
quadratic f = fitParams . quadraticFit f

-- | Like 'quadratic' but returns complete 'Fit'.
--
-- >>> PP $ quadraticFit id input2
-- Fit (V3 (-0.00589) 2.0313 0.87155) (V3 0.09281 0.41070 0.37841) 2 0.25910
--
-- >>> PP $ quadraticFit id input3
-- Fit (V3 1.00000 0.00000 2.0000) (V3 0.00000 0.00000 0.00000) 1 0.00000
--
quadraticFit :: F.Foldable f => (a -> (Double, Double)) -> f a -> Fit V3
quadraticFit f = quadraticFit' . quadRegAcc f

-- | Weighted quadratic regression.
--
-- >>> let input2w = [(0.1, 1.2, 1), (1.3, 3.1, 1), (1.9, 4.9, 1), (3.0, 7.1, 1/4), (4.1, 9.0, 1/4)]
-- >>> PP $ quadraticWithWeights id input2w
-- Fit (V3 0.02524 1.9144 0.91792) (V3 0.10775 0.42106 0.35207) 2 0.21484
--
quadraticWithWeights :: F.Foldable f => (a -> (Double, Double, Double)) -> f a -> Fit V3
quadraticWithWeights f = quadraticFit' . quadRegAccW f

-- | Quadratic regression with y-errors.
--
-- >>> let input2y = [(0.1, 1.2, 0.12), (1.3, 3.1, 0.31), (1.9, 4.9, 0.49), (3.0, 7.1, 0.71), (4.1, 9.0, 0.9)]
-- >>> PP $ quadraticWithYerrors id input2y
-- Fit (V3 0.08776 1.6667 1.0228) (V3 0.10131 0.31829 0.11917) 2 1.5398
--
quadraticWithYerrors :: F.Foldable f => (a -> (Double, Double, Double)) -> f a -> Fit V3
quadraticWithYerrors f = quadraticWithWeights f' where
    f' a = case f a of
        (x, y, dy) -> (x, y, recip (dy * dy))

-- | Iterative quadratic regression with x and y errors.
--
-- /Orear, J. (1982). Least squares when both variables have uncertainties. American Journal of Physics, 50(10), 912–916. doi:10.1119\/1.12972/
--
quadraticWithXYerrors
    :: F.Foldable f
    => (a -> (Double, Double, Double, Double))  -- ^ \(x_i, y_i, \delta x_i, \delta y_i\)
    -> f a                                      -- ^ data
    -> NE.NonEmpty (Fit V3)
quadraticWithXYerrors f xs = iterate1 go fit0 where
    fit0   = quadraticWithYerrors (\a -> case f a of (x,y,_,dy) -> (x,y,dy)) xs
    go fit = quadraticWithWeights (\a -> case f a of (x,y,dx,dy) -> (x,y,recip $ sq ((2 * p1 * x + p2) * dx) + sq dy)) xs where
        V3 p1 p2 _ = fitParams fit

-- | Calculate quadratic fit from 'QuadRegAcc'.
quadraticFit' :: QuadRegAcc -> Fit V3
quadraticFit' QuadRegAcc {..} = Fit params errors ndf wssr where
    matrix@(SM33 a11
                 _   a22
                 _   _   a33) = inv (SM33 x4
                                          x3 x2
                                          x2  x w)

    params@(V3 a b c) = mult matrix (V3 x2y xy y)

    errors = V3 sa sb sc

    wssr = max 0 (y2 - a * x2y - b * xy - c * y)
    ndf  = qra_n - 3
    ndf' = fromIntegral ndf :: Double

    sa  = sqrt (a11 * wssr / ndf')
    sb  = sqrt (a22 * wssr / ndf')
    sc  = sqrt (a33 * wssr / ndf')

    w   = getKBN qra_w
    x   = getKBN qra_x
    x2  = getKBN qra_x2
    x3  = getKBN qra_x3
    x4  = getKBN qra_x4
    y   = getKBN qra_y
    xy  = getKBN qra_xy
    x2y = getKBN qra_x2y
    y2  = getKBN qra_y2

    -- is it useful?
    -- r2 = (n * xy - x*y) / (sqrt (n * x2 - x*x) * sqrt (n * y2 - y*y))

quadRegAcc :: F.Foldable f => (a -> (Double, Double)) -> f a -> QuadRegAcc
quadRegAcc f = F.foldl' (\acc a -> case f a of (x,y) -> addQuadReg acc x y) zeroQuadRegAcc

quadRegAccW :: F.Foldable f => (a -> (Double, Double, Double)) -> f a -> QuadRegAcc
quadRegAccW f = F.foldl' (\acc a -> case f a of (x,y,w) -> addQuadRegW acc x y w) zeroQuadRegAcc

-------------------------------------------------------------------------------
-- Levenberg–Marquardt 1
-------------------------------------------------------------------------------

-- | Levenberg–Marquardt for functions with one parameter.
--
-- See 'levenbergMarquardt2' for examples, this is very similar.
--
-- For example we can fit \(f = x \mapsto \beta x + 1\), its derivative is \(\partial_\beta f = x \mapsto x\).
--
-- >>> let scale a (x, y) = (y, a * x + 1, x)
-- >>> PP $ NE.last $ levenbergMarquardt1 scale 1 input2
-- Fit 1.9685 0.04735 4 0.27914
--
-- Not bad, but worse then linear fit which fits the intercept point too.
--
levenbergMarquardt1
    :: F.Foldable f
    => (Double -> a -> (Double, Double, Double))  -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \partial_\beta f(\beta, x)\)
    -> Double                                     -- ^ initial parameter, \(\beta_0\)
    -> f a                                        -- ^ data, \(d\)
    -> NE.NonEmpty (Fit Double)                   -- ^ non-empty list of iteration results
levenbergMarquardt1 f b0 xs = loop lambda0 b0 acc0 where
    acc0 = calcAcc b0

    lambda0 = c11
      where
        c11 = getKBN $ lm1_c11 acc0

    calcAcc beta = F.foldl' (\acc p -> case f beta p of (y, g, d) -> addLM1Acc acc y g d) zeroLM1Acc xs

    loop lambda beta acc
        | lmStop lambda wssr wssr'
        = Fit beta errors ndf wssr NE.:| []

        | wssr' >= wssr
        = loop (lambda * 10) beta acc

        | otherwise
        = Fit beta errors ndf wssr `NE.cons` loop (lambda / 10) beta' acc'

      where
        lambda1 z = (1 + lambda) * z
        matrix    = inv (lambda1 c11)
        delta     = mult matrix z1

        beta' = add beta delta
        acc'  = calcAcc beta'

        a11    = inv c11
        errors = sa

        wssr  = max 0 $ getKBN $ lm1_wssr acc
        wssr' =         getKBN $ lm1_wssr acc'

        ndf   = lm1_n acc - 1
        ndf'  = fromIntegral ndf :: Double

        sa   = sqrt (a11 * wssr / ndf')

        c11 = getKBN $ lm1_c11 acc
        z1  = getKBN $ lm1_z1 acc

-- | 'levenbergMarquardt1' with weights.
levenbergMarquardt1WithWeights
    :: F.Foldable f
    => (Double -> a -> (Double, Double, Double, Double))  -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \partial_\beta f(\beta, x), w_i\)
    -> Double                                             -- ^ initial parameter, \(\beta_0\)
    -> f a                                                -- ^ data, \(d\)
    -> NE.NonEmpty (Fit Double)                           -- ^ non-empty list of iteration results
levenbergMarquardt1WithWeights f b0 xs = loop lambda0 b0 acc0 where
    acc0 = calcAcc b0

    lambda0 = c11
      where
        c11 = getKBN $ lm1_c11 acc0

    calcAcc beta = F.foldl' (\acc p -> case f beta p of (y, g, d, w) -> addLM1AccW acc y g d w) zeroLM1Acc xs

    loop lambda beta acc
        | lmStop lambda wssr wssr'
        = Fit beta errors ndf wssr NE.:| []

        | wssr' >= wssr
        = loop (lambda * 10) beta acc

        | otherwise
        = Fit beta errors ndf wssr `NE.cons` loop (lambda / 10) beta' acc'

      where
        lambda1 z = (1 + lambda) * z
        matrix    = inv (lambda1 c11)
        delta     = mult matrix z1

        beta' = add beta delta
        acc'  = calcAcc beta'

        a11    = inv c11
        errors = sa

        wssr  = max 0 $ getKBN $ lm1_wssr acc
        wssr' =         getKBN $ lm1_wssr acc'

        ndf   = lm1_n acc - 1
        ndf'  = fromIntegral ndf :: Double

        sa   = sqrt (a11 * wssr / ndf')

        c11 = getKBN $ lm1_c11 acc
        z1  = getKBN $ lm1_z1 acc

-- | 'levenbergMarquardt1' with Y-errors.
levenbergMarquardt1WithYerrors
    :: F.Foldable f
    => (Double -> a -> (Double, Double, Double, Double))  -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \partial_\beta f(\beta, x), \delta y_i\)
    -> Double                                             -- ^ initial parameter, \(\beta_0\)
    -> f a                                                -- ^ data, \(d\)
    -> NE.NonEmpty (Fit Double)                           -- ^ non-empty list of iteration results
levenbergMarquardt1WithYerrors f = levenbergMarquardt1WithWeights f' where
    f' beta x = case f beta x of (y, fbetax, grad, dy) -> (y, fbetax, grad, recip $ sq dy)

-- | 'levenbergMarquardt1' with XY-errors.
levenbergMarquardt1WithXYerrors
    :: F.Foldable f
    => (Double -> a -> (Double, Double, Double, Double, Double, Double))  -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \partial_\beta f(\beta, x), \partial_x f(\beta, x), \delta x_i, \delta y_i\)
    -> Double                                                             -- ^ initial parameter, \(\beta_0\)
    -> f a                                                                -- ^ data, \(d\)
    -> NE.NonEmpty (Fit Double)                                           -- ^ non-empty list of iteration results
levenbergMarquardt1WithXYerrors g = levenbergMarquardt1WithWeights g' where
    g' beta x = case g beta x of (y, fbetax, grad, f', dx, dy) -> (y, fbetax, grad, recip $ sq (f' * dx) + sq dy)

data LM1Acc = LM1Acc
    { lm1_n    :: !Int
    , lm1_c11  :: !KBN
    , lm1_z1   :: !KBN
    , lm1_wssr :: !KBN
    }
  deriving Show

zeroLM1Acc :: LM1Acc
zeroLM1Acc = LM1Acc 0 zeroKBN zeroKBN zeroKBN

addLM1Acc :: LM1Acc -> Double -> Double -> Double -> LM1Acc
addLM1Acc LM1Acc {..} y f d1 = LM1Acc
    { lm1_n    = lm1_n + 1
    , lm1_c11  = addKBN lm1_c11  (d1 * d1)
    , lm1_z1   = addKBN lm1_z1   (d1 * res)
    , lm1_wssr = addKBN lm1_wssr (res * res)
    }
  where
    res = y - f

addLM1AccW :: LM1Acc -> Double -> Double -> Double -> Double -> LM1Acc
addLM1AccW LM1Acc {..} y f d1 w = LM1Acc
    { lm1_n    = lm1_n + 1
    , lm1_c11  = addKBN lm1_c11  (w * d1 * d1)
    , lm1_z1   = addKBN lm1_z1   (w * d1 * res)
    , lm1_wssr = addKBN lm1_wssr (w * res * res)
    }
  where
    res = y - f

-------------------------------------------------------------------------------
-- Levenberg–Marquardt 2
-------------------------------------------------------------------------------

-- | Levenberg–Marquardt for functions with two parameters.
--
-- You can use this sledgehammer to do a a linear fit:
--
-- >>> let lin (V2 a b) (x, y) = (y, a * x + b, V2 x 1)
--
-- We can then use 'levenbergMarquardt2' to find a fit:
--
-- >>> PP $ levenbergMarquardt2 lin (V2 1 1) input2
-- Fit (V2 1.00000 1.00000) (V2 1.0175 2.5385) 3 29.470
-- Fit (V2 1.0181 1.0368) (V2 0.98615 2.4602) 3 27.681
-- Fit (V2 1.1557 1.2988) (V2 0.75758 1.8900) 3 16.336
-- Fit (V2 1.5463 1.6577) (V2 0.29278 0.73043) 3 2.4400
-- Fit (V2 1.9129 1.1096) (V2 0.11033 0.27524) 3 0.34645
-- Fit (V2 2.0036 0.89372) (V2 0.09552 0.23830) 3 0.25970
-- Fit (V2 2.0063 0.88687) (V2 0.09550 0.23826) 3 0.25962
-- Fit (V2 2.0063 0.88685) (V2 0.09550 0.23826) 3 0.25962
--
-- This is the same result what 'linearFit' returns:
--
-- >>> PP $ linearFit id input2
-- Fit (V2 2.0063 0.88685) (V2 0.09550 0.23826) 3 0.25962
--
-- == Using AD
--
-- You can use @ad@ to calculate derivatives for you.
--
-- >>> import qualified Numeric.AD.Mode.Reverse.Double as AD
--
-- We need a ('Traversable') homogenic triple to represent the two parameters and @x@:
--
-- >>> data H3 a = H3 a a a deriving (Functor, Foldable, Traversable)
--
-- Then we define a function @ad@ can operate with:
--
-- >>> let linearF (H3 a b x) = a * x + b
--
-- which we can use to fit the curve in generic way:
--
-- >>> let lin' (V2 a b) (x, y) = case AD.grad' linearF (H3 a b x) of (f, H3 da db _f') -> (y, f, V2 da db)
-- >>> PP $ levenbergMarquardt2 lin' (V2 1 1) input2
-- Fit (V2 1.00000 1.00000) (V2 1.0175 2.5385) 3 29.470
-- Fit (V2 1.0181 1.0368) (V2 0.98615 2.4602) 3 27.681
-- Fit (V2 1.1557 1.2988) (V2 0.75758 1.8900) 3 16.336
-- Fit (V2 1.5463 1.6577) (V2 0.29278 0.73043) 3 2.4400
-- Fit (V2 1.9129 1.1096) (V2 0.11033 0.27524) 3 0.34645
-- Fit (V2 2.0036 0.89372) (V2 0.09552 0.23830) 3 0.25970
-- Fit (V2 2.0063 0.88687) (V2 0.09550 0.23826) 3 0.25962
-- Fit (V2 2.0063 0.88685) (V2 0.09550 0.23826) 3 0.25962
--
-- == Non-polynomial example
--
-- We can fit other curves too, for example an example from Wikipedia
-- https://en.wikipedia.org/wiki/Gauss%E2%80%93Newton_algorithm#Example
--
-- >>> let rateF (H3 vmax km s) = (vmax * s) / (km + s)
-- >>> let rateF' (V2 vmax km) (x, y) = case AD.grad' rateF (H3 vmax km x) of (f, H3 vmax' km' _) -> (y, f, V2 vmax' km')
-- >>> let input = zip [0.038,0.194,0.425,0.626,1.253,2.500,3.740] [0.050,0.127,0.094,0.2122,0.2729,0.2665,0.3317]
-- >>> PP $ levenbergMarquardt2 rateF' (V2 0.9 0.2) input
-- Fit (V2 0.90000 0.20000) (V2 0.43304 0.43936) 5 1.4455
-- Fit (V2 0.83306 0.25278) (V2 0.39164 0.49729) 5 1.0055
-- Fit (V2 0.59437 0.43508) (V2 0.21158 0.53403) 5 0.18832
-- Fit (V2 0.39687 0.56324) (V2 0.05723 0.25666) 5 0.01062
-- Fit (V2 0.36289 0.56104) (V2 0.04908 0.24007) 5 0.00784
-- Fit (V2 0.36190 0.55662) (V2 0.04887 0.23843) 5 0.00784
-- Fit (V2 0.36184 0.55629) (V2 0.04885 0.23830) 5 0.00784
-- Fit (V2 0.36184 0.55627) (V2 0.04885 0.23829) 5 0.00784
--
-- We get the same result as in the article: 0.362 and 0.556
--
-- The algorithm terminates when a scaling parameter \(\lambda\) becomes larger than 1e20 or smaller than 1e-20, or relative WSSR change is smaller than 1e-10, or sum-of-squared-residuals candidate becomes @NaN@ (i.e. when it would start to produce garbage).
-- You may want to terminate sooner, Numerical Recipes suggest to stop when WSSR decreases by a neglible amount absolutely or fractionally.
--
levenbergMarquardt2
    :: F.Foldable f
    => (V2 -> a -> (Double, Double, V2))  -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \nabla_\beta f(\beta, x)\)
    -> V2                                 -- ^ initial parameters, \(\beta_0\)
    -> f a                                -- ^ data, \(d\)
    -> NE.NonEmpty (Fit V2)               -- ^ non-empty list of iteration results
levenbergMarquardt2 f b0 xs = loop lambda0 b0 acc0 where
    acc0 = calcAcc b0

    calcAcc beta = F.foldl' (\acc p -> case f beta p of (y, g, d) -> addLM2Acc acc y g d) zeroLM2Acc xs

    lambda0 = max l1 l2
      where
        V2 l1 l2 = eigenSM22 c11 c12 c22
        c11 = getKBN $ lm2_c11 acc0
        c12 = getKBN $ lm2_c11 acc0
        c22 = getKBN $ lm2_c22 acc0

    loop lambda beta acc
        | lmStop lambda wssr wssr'
        = Fit beta errors ndf wssr NE.:| []

        | wssr' >= wssr
        = loop (lambda * 10) beta acc

        | otherwise
        = Fit beta errors ndf wssr `NE.cons` loop (lambda / 10) beta' acc'

      where
        lambda1 z = (1 + lambda) * z
        matrix    = inv (SM22 (lambda1 c11) c12 (lambda1 c22))
        delta     = mult matrix (V2 z1 z2)

        beta' = add beta delta
        acc'  = calcAcc beta'

        SM22 a11 _ a22 = inv (SM22 c11 c12 c22)
        errors = V2 sa sb

        wssr  = max 0 $ getKBN $ lm2_wssr acc
        wssr' =         getKBN $ lm2_wssr acc'

        ndf   = lm2_n acc - 2
        ndf'  = fromIntegral ndf :: Double

        sa   = sqrt (a11 * wssr / ndf')
        sb   = sqrt (a22 * wssr / ndf')

        c11 = getKBN $ lm2_c11 acc
        c12 = getKBN $ lm2_c12 acc
        c22 = getKBN $ lm2_c22 acc
        z1  = getKBN $ lm2_z1 acc
        z2  = getKBN $ lm2_z2 acc

-- | 'levenbergMarquardt2' with weights.
--
-- Because 'levenbergMarquardt2' is an iterative algorithm,
-- not only we can use it to fit curves with known y-errors ('levenbergMarquardt2WithYerrors'),
-- but also with both x and y-errors ('levenbergMarquardt2WithXYerrors').
--
levenbergMarquardt2WithWeights
    :: F.Foldable f
    => (V2 -> a -> (Double, Double, V2, Double))  -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \nabla_\beta f(\beta, x), w_i\)
    -> V2                                         -- ^ initial parameters, \(\beta_0\)
    -> f a                                        -- ^ data, \(d\)
    -> NE.NonEmpty (Fit V2)                       -- ^ non-empty list of iteration results
levenbergMarquardt2WithWeights f b0 xs = loop lambda0 b0 acc0 where
    acc0 = calcAcc b0

    lambda0 = max l1 l2
      where
        V2 l1 l2 = eigenSM22 c11 c12 c22
        c11 = getKBN $ lm2_c11 acc0
        c12 = getKBN $ lm2_c11 acc0
        c22 = getKBN $ lm2_c22 acc0

    calcAcc beta = F.foldl' (\acc p -> case f beta p of (y, g, d, w) -> addLM2AccW acc y g d w) zeroLM2Acc xs

    loop lambda beta acc
        | lmStop lambda wssr wssr'
        = Fit beta errors ndf wssr NE.:| []

        | wssr' >= wssr
        = loop (lambda * 10) beta acc

        | otherwise
        = Fit beta errors ndf wssr `NE.cons` loop (lambda / 10) beta' acc'

      where
        lambda1 z = (1 + lambda) * z
        matrix    = inv (SM22 (lambda1 c11) c12 (lambda1 c22))
        delta     = mult matrix (V2 z1 z2)

        beta' = add beta delta
        acc'  = calcAcc beta'

        SM22 a11 _ a22 = inv (SM22 c11 c12 c22)
        errors = V2 sa sb

        wssr  = max 0 $ getKBN $ lm2_wssr acc
        wssr' =         getKBN $ lm2_wssr acc'

        ndf   = lm2_n acc - 2
        ndf'  = fromIntegral ndf :: Double

        sa   = sqrt (a11 * wssr / ndf')
        sb   = sqrt (a22 * wssr / ndf')

        c11 = getKBN $ lm2_c11 acc
        c12 = getKBN $ lm2_c12 acc
        c22 = getKBN $ lm2_c22 acc
        z1  = getKBN $ lm2_z1 acc
        z2  = getKBN $ lm2_z2 acc

-- | 'levenbergMarquardt2' with Y-errors.
levenbergMarquardt2WithYerrors
    :: F.Foldable f
    => (V2 -> a -> (Double, Double, V2, Double))  -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \nabla_\beta f(\beta, x), \delta y_i\)
    -> V2                                         -- ^ initial parameters, \(\beta_0\)
    -> f a                                        -- ^ data, \(d\)
    -> NE.NonEmpty (Fit V2)                       -- ^ non-empty list of iteration results
levenbergMarquardt2WithYerrors f = levenbergMarquardt2WithWeights f' where
    f' beta x = case f beta x of (y, fbetax, grad, dy) -> (y, fbetax, grad, recip $ sq dy)

-- | 'levenbergMarquardt2' with XY-errors.
levenbergMarquardt2WithXYerrors
    :: F.Foldable f
    => (V2 -> a -> (Double, Double, V2, Double, Double, Double))  -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \nabla_\beta f(\beta, x), \partial_x f(\beta, x), \delta x_i, \delta y_i\)
    -> V2                                                         -- ^ initial parameters, \(\beta_0\)
    -> f a                                                        -- ^ data, \(d\)
    -> NE.NonEmpty (Fit V2)                                       -- ^ non-empty list of iteration results
levenbergMarquardt2WithXYerrors g = levenbergMarquardt2WithWeights g' where
    g' beta x = case g beta x of (y, fbetax, grad, f', dx, dy) -> (y, fbetax, grad, recip $ sq (f' * dx) + sq dy)

data LM2Acc = LM2Acc
    { lm2_n    :: !Int
    , lm2_c11  :: !KBN
    , lm2_c12  :: !KBN
    , lm2_c22  :: !KBN
    , lm2_z1   :: !KBN
    , lm2_z2   :: !KBN
    , lm2_wssr :: !KBN
    }
  deriving Show

zeroLM2Acc :: LM2Acc
zeroLM2Acc = LM2Acc 0 zeroKBN zeroKBN zeroKBN zeroKBN zeroKBN zeroKBN

addLM2Acc :: LM2Acc -> Double -> Double -> V2 -> LM2Acc
addLM2Acc LM2Acc {..} y f (V2 d1 d2) = LM2Acc
    { lm2_n    = lm2_n + 1
    , lm2_c11  = addKBN lm2_c11  (d1 * d1)
    , lm2_c12  = addKBN lm2_c12  (d1 * d2)
    , lm2_c22  = addKBN lm2_c22  (d2 * d2)
    , lm2_z1   = addKBN lm2_z1   (d1 * res)
    , lm2_z2   = addKBN lm2_z2   (d2 * res)
    , lm2_wssr = addKBN lm2_wssr (res * res)
    }
  where
    res = y - f

addLM2AccW :: LM2Acc -> Double -> Double -> V2 -> Double -> LM2Acc
addLM2AccW LM2Acc {..} y f (V2 d1 d2) w = LM2Acc
    { lm2_n    = lm2_n + 1
    , lm2_c11  = addKBN lm2_c11  (w * d1 * d1)
    , lm2_c12  = addKBN lm2_c12  (w * d1 * d2)
    , lm2_c22  = addKBN lm2_c22  (w * d2 * d2)
    , lm2_z1   = addKBN lm2_z1   (w * d1 * res)
    , lm2_z2   = addKBN lm2_z2   (w * d2 * res)
    , lm2_wssr = addKBN lm2_wssr (w * res * res)
    }
  where
    res = y - f

-------------------------------------------------------------------------------
-- Levenberg–Marquardt 3
-------------------------------------------------------------------------------

-- | Levenberg–Marquardt for functions with three parameters.
--
-- See 'levenbergMarquardt2' for examples, this is very similar.
--
-- >>> let quad (V3 a b c) (x, y) = (y, a * x * x + b * x + c, V3 (x * x) x 1)
-- >>> PP $ NE.last $ levenbergMarquardt3 quad (V3 2 2 2) input3
-- Fit (V3 1.00000 0.00000 2.0000) (V3 0.00000 0.00000 0.00000) 1 0.00000
--
-- Same as quadratic fit, just less direct:
--
-- >>> PP $ quadraticFit id input3
-- Fit (V3 1.00000 0.00000 2.0000) (V3 0.00000 0.00000 0.00000) 1 0.00000
--
levenbergMarquardt3
    :: F.Foldable f
    => (V3 -> a -> (Double, Double, V3))  -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \nabla_\beta f(\beta, x)\)
    -> V3                                 -- ^ initial parameters, \(\beta_0\)
    -> f a                                -- ^ data, \(d\)
    -> NE.NonEmpty (Fit V3)               -- ^ non-empty list of iteration results
levenbergMarquardt3 f b0 xs = loop lambda0 b0 acc0 where
    acc0 = calcAcc b0

    calcAcc beta = F.foldl' (\acc p -> case f beta p of (y, g, d) -> addLM3Acc acc y g d) zeroLM3Acc xs

    -- frobenius norm is larger than largest eigen value.
    -- calculating the eigen values for 3x3 (symmetric) matrix is becoming complicated.
    lambda0 = sqrt $ sumKBN [ sq c11
                            , 2 * sq c12, sq c22
                            , 2 * sq c13, 2 * sq c23, sq c33]
      where
        c11 = getKBN $ lm3_c11 acc0
        c12 = getKBN $ lm3_c12 acc0
        c13 = getKBN $ lm3_c13 acc0
        c22 = getKBN $ lm3_c22 acc0
        c23 = getKBN $ lm3_c23 acc0
        c33 = getKBN $ lm3_c33 acc0

    loop lambda beta acc
        | lmStop lambda wssr wssr'
        = Fit beta errors ndf wssr NE.:| []

        | wssr' >= wssr
        = loop (lambda * 10) beta acc

        | otherwise
        = Fit beta errors ndf wssr `NE.cons` loop (lambda / 10) beta' acc'

      where
        lambda1 z = (1 + lambda) * z
        matrix    = inv (SM33 (lambda1 c11)
                               c12          (lambda1 c22)
                               c13          c23           (lambda1 c33))
        delta     = mult matrix (V3 z1 z2 z3)

        beta' = add beta delta
        acc'  = calcAcc beta'

        SM33 a11
             _  a22
             _  _   a33 = inv (SM33 c11
                                    c12 c22
                                    c13 c23 c33)
        errors = V3 sa sb sc

        wssr  = max 0 $ getKBN $ lm3_wssr acc
        wssr' =         getKBN $ lm3_wssr acc'

        ndf   = lm3_n acc - 3
        ndf'  = fromIntegral ndf :: Double

        sa   = sqrt (a11 * wssr / ndf')
        sb   = sqrt (a22 * wssr / ndf')
        sc   = sqrt (a33 * wssr / ndf')

        c11 = getKBN $ lm3_c11 acc
        c12 = getKBN $ lm3_c12 acc
        c13 = getKBN $ lm3_c13 acc
        c22 = getKBN $ lm3_c22 acc
        c23 = getKBN $ lm3_c23 acc
        c33 = getKBN $ lm3_c33 acc
        z1  = getKBN $ lm3_z1 acc
        z2  = getKBN $ lm3_z2 acc
        z3  = getKBN $ lm3_z3 acc

-- | 'levenbergMarquardt3' with weights.
levenbergMarquardt3WithWeights
    :: F.Foldable f
    => (V3 -> a -> (Double, Double, V3, Double))  -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \nabla_\beta f(\beta, x), w_i\)
    -> V3                                         -- ^ initial parameters, \(\beta_0\)
    -> f a                                        -- ^ data, \(d\)
    -> NE.NonEmpty (Fit V3)                       -- ^ non-empty list of iteration results
levenbergMarquardt3WithWeights f b0 xs = loop lambda0 b0 acc0 where
    acc0 = calcAcc b0

    lambda0 = sqrt $ sumKBN [ sq c11
                            , 2 * sq c12, sq c22
                            , 2 * sq c13, 2 * sq c23, sq c33]
      where
        c11 = getKBN $ lm3_c11 acc0
        c12 = getKBN $ lm3_c12 acc0
        c13 = getKBN $ lm3_c13 acc0
        c22 = getKBN $ lm3_c22 acc0
        c23 = getKBN $ lm3_c23 acc0
        c33 = getKBN $ lm3_c33 acc0

    calcAcc beta = F.foldl' (\acc p -> case f beta p of (y, g, d, w) -> addLM3AccW acc y g d w) zeroLM3Acc xs

    loop lambda beta acc
        | lmStop lambda wssr wssr'
        = Fit beta errors ndf wssr NE.:| []

        | wssr' >= wssr
        = loop (lambda * 10) beta acc

        | otherwise
        = Fit beta errors ndf wssr `NE.cons` loop (lambda / 10) beta' acc'

      where
        lambda1 z = (1 + lambda) * z
        matrix    = inv (SM33 (lambda1 c11)
                               c12          (lambda1 c22)
                               c13          c23           (lambda1 c33))
        delta     = mult matrix (V3 z1 z2 z3)

        beta' = add beta delta
        acc'  = calcAcc beta'

        SM33 a11 _ _ a22 _ a33 = inv (SM33 c11 c12 c13 c22 c23 c33)
        errors = V3 sa sb sc

        wssr  = max 0 $ getKBN $ lm3_wssr acc
        wssr' =         getKBN $ lm3_wssr acc'

        ndf   = lm3_n acc - 3
        ndf'  = fromIntegral ndf :: Double

        sa   = sqrt (a11 * wssr / ndf')
        sb   = sqrt (a22 * wssr / ndf')
        sc   = sqrt (a33 * wssr / ndf')

        c11 = getKBN $ lm3_c11 acc
        c12 = getKBN $ lm3_c12 acc
        c13 = getKBN $ lm3_c13 acc
        c22 = getKBN $ lm3_c22 acc
        c23 = getKBN $ lm3_c23 acc
        c33 = getKBN $ lm3_c33 acc
        z1  = getKBN $ lm3_z1 acc
        z2  = getKBN $ lm3_z2 acc
        z3  = getKBN $ lm3_z3 acc

-- | 'levenbergMarquardt3' with Y-errors.
levenbergMarquardt3WithYerrors
    :: F.Foldable f
    => (V3 -> a -> (Double, Double, V3, Double))  -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \nabla_\beta f(\beta, x), \delta y_i\)
    -> V3                                         -- ^ initial parameters, \(\beta_0\)
    -> f a                                        -- ^ data, \(d\)
    -> NE.NonEmpty (Fit V3)                       -- ^ non-empty list of iteration results
levenbergMarquardt3WithYerrors f = levenbergMarquardt3WithWeights f' where
    f' beta x = case f beta x of (y, fbetax, grad, dy) -> (y, fbetax, grad, recip $ sq dy)

-- | 'levenbergMarquardt3' with XY-errors.
levenbergMarquardt3WithXYerrors
    :: F.Foldable f
    => (V3 -> a -> (Double, Double, V3, Double, Double, Double))  -- ^ \(\beta, d_i \mapsto y_i, f(\beta, x_i), \nabla_\beta f(\beta, x), \partial_x f(\beta, x), \delta x_i, \delta y_i\)
    -> V3                                                         -- ^ initial parameters, \(\beta_0\)
    -> f a                                                        -- ^ data, \(d\)
    -> NE.NonEmpty (Fit V3)                                       -- ^ non-empty list of iteration results
levenbergMarquardt3WithXYerrors g = levenbergMarquardt3WithWeights g' where
    g' beta x = case g beta x of (y, fbetax, grad, f', dx, dy) -> (y, fbetax, grad, recip $ sq (f' * dx) + sq dy)

data LM3Acc = LM3Acc
    { lm3_n    :: !Int
    , lm3_c11  :: !KBN
    , lm3_c12  :: !KBN
    , lm3_c13  :: !KBN
    , lm3_c22  :: !KBN
    , lm3_c23  :: !KBN
    , lm3_c33  :: !KBN
    , lm3_z1   :: !KBN
    , lm3_z2   :: !KBN
    , lm3_z3   :: !KBN
    , lm3_wssr :: !KBN
    }
  deriving Show

zeroLM3Acc :: LM3Acc
zeroLM3Acc = LM3Acc 0 zeroKBN zeroKBN zeroKBN zeroKBN zeroKBN zeroKBN zeroKBN zeroKBN zeroKBN zeroKBN

addLM3Acc :: LM3Acc -> Double -> Double -> V3 -> LM3Acc
addLM3Acc LM3Acc {..} y f (V3 d1 d2 d3) = LM3Acc
    { lm3_n    = lm3_n + 1
    , lm3_c11  = addKBN lm3_c11  (d1 * d1)
    , lm3_c12  = addKBN lm3_c12  (d1 * d2)
    , lm3_c13  = addKBN lm3_c12  (d1 * d3)
    , lm3_c22  = addKBN lm3_c22  (d2 * d2)
    , lm3_c23  = addKBN lm3_c22  (d2 * d3)
    , lm3_c33  = addKBN lm3_c22  (d3 * d3)
    , lm3_z1   = addKBN lm3_z1   (d1 * res)
    , lm3_z2   = addKBN lm3_z2   (d2 * res)
    , lm3_z3   = addKBN lm3_z3   (d3 * res)
    , lm3_wssr = addKBN lm3_wssr (res * res)
    }
  where
    res = y - f

addLM3AccW :: LM3Acc -> Double -> Double -> V3 -> Double -> LM3Acc
addLM3AccW LM3Acc {..} y f (V3 d1 d2 d3) w = LM3Acc
    { lm3_n    = lm3_n + 1
    , lm3_c11  = addKBN lm3_c11  (w * d1 * d1)
    , lm3_c12  = addKBN lm3_c12  (w * d1 * d2)
    , lm3_c13  = addKBN lm3_c12  (w * d1 * d3)
    , lm3_c22  = addKBN lm3_c22  (w * d2 * d2)
    , lm3_c23  = addKBN lm3_c22  (w * d2 * d3)
    , lm3_c33  = addKBN lm3_c22  (w * d3 * d3)
    , lm3_z1   = addKBN lm3_z1   (w * d1 * res)
    , lm3_z2   = addKBN lm3_z2   (w * d2 * res)
    , lm3_z3   = addKBN lm3_z3   (w * d3 * res)
    , lm3_wssr = addKBN lm3_wssr (w * res * res)
    }
  where
    res = y - f

-------------------------------------------------------------------------------
-- Output
-------------------------------------------------------------------------------

-- | Result of a curve fit.
data Fit v = Fit
    { fitParams :: !v       -- ^ fit parameters
    , fitErrors :: !v       -- ^ asympotic standard errors, /assuming a good fit/
    , fitNDF    :: !Int     -- ^ number of degrees of freedom
    , fitWSSR   :: !Double  -- ^ sum of squares of residuals
    }
  deriving Show

instance (NFData v) => NFData (Fit v) where
    rnf (Fit p e _ _) = rnf p `seq` rnf e

-------------------------------------------------------------------------------
-- LinRegAcc
-------------------------------------------------------------------------------

-- | Linear regression accumulator.
data LinRegAcc = LinRegAcc
    { lra_n  :: {-# UNPACK #-} !Int  -- ^ \(n\)
    , lra_w  :: {-# UNPACK #-} !KBN  -- ^ \(\sum w_i\)
    , lra_x  :: {-# UNPACK #-} !KBN  -- ^ \(\sum x_i \)
    , lra_x2 :: {-# UNPACK #-} !KBN  -- ^ \(\sum x_i^2 \)
    , lra_y  :: {-# UNPACK #-} !KBN  -- ^ \(\sum y_i \)
    , lra_xy :: {-# UNPACK #-} !KBN  -- ^ \(\sum x_i y_i \)
    , lra_y2 :: {-# UNPACK #-} !KBN  -- ^ \(\sum y_i^2 \)
    }
  deriving Show

instance NFData LinRegAcc where
    rnf LinRegAcc {} = ()

-- | All-zeroes 'LinRegAcc'.
zeroLinRegAcc :: LinRegAcc
zeroLinRegAcc = LinRegAcc 0 zeroKBN zeroKBN zeroKBN zeroKBN zeroKBN zeroKBN

-- | Add a point to linreg accumulator.
addLinReg
    :: LinRegAcc
    -> Double   -- ^ x
    -> Double   -- ^ y
    -> LinRegAcc
addLinReg LinRegAcc {..} x y = LinRegAcc
    { lra_n  = lra_n + 1
    , lra_w  = addKBN lra_w  1
    , lra_x  = addKBN lra_x  x
    , lra_x2 = addKBN lra_x2 (x * x)
    , lra_y  = addKBN lra_y  y
    , lra_xy = addKBN lra_xy (x * y)
    , lra_y2 = addKBN lra_y2 (y * y)
    }

-- | Add a weighted point to linreg accumulator.
addLinRegW
    :: LinRegAcc
    -> Double   -- ^ x
    -> Double   -- ^ y
    -> Double   -- ^ w
    -> LinRegAcc
addLinRegW LinRegAcc {..} x y w = LinRegAcc
    { lra_n  = lra_n + 1
    , lra_w  = addKBN lra_w  w
    , lra_x  = addKBN lra_x  (w * x)
    , lra_x2 = addKBN lra_x2 (w * x * x)
    , lra_y  = addKBN lra_y  (w * y)
    , lra_xy = addKBN lra_xy (w * x * y)
    , lra_y2 = addKBN lra_y2 (w * y * y)
    }

-------------------------------------------------------------------------------
-- QuadRegAcc
-------------------------------------------------------------------------------

-- | Quadratic regression accumulator.
data QuadRegAcc = QuadRegAcc
    { qra_n   :: {-# UNPACK #-} !Int  -- ^ \(n\)
    , qra_w   :: {-# UNPACK #-} !KBN  -- ^ \(\sum w_i\)
    , qra_x   :: {-# UNPACK #-} !KBN  -- ^ \(\sum x_i \)
    , qra_x2  :: {-# UNPACK #-} !KBN  -- ^ \(\sum x_i^2 \)
    , qra_x3  :: {-# UNPACK #-} !KBN  -- ^ \(\sum x_i^3 \)
    , qra_x4  :: {-# UNPACK #-} !KBN  -- ^ \(\sum x_i^4 \)
    , qra_y   :: {-# UNPACK #-} !KBN  -- ^ \(\sum y_i \)
    , qra_xy  :: {-# UNPACK #-} !KBN  -- ^ \(\sum x_i y_i \)
    , qra_x2y :: {-# UNPACK #-} !KBN  -- ^ \(\sum x_i^2 y_i \)
    , qra_y2  :: {-# UNPACK #-} !KBN  -- ^ \(\sum y_i^2 \)
    }
  deriving Show

instance NFData QuadRegAcc where
    rnf QuadRegAcc {} = ()

-- | All-zeroes 'QuadRegAcc'.
zeroQuadRegAcc :: QuadRegAcc
zeroQuadRegAcc = QuadRegAcc 0 zeroKBN zeroKBN zeroKBN zeroKBN zeroKBN zeroKBN zeroKBN zeroKBN zeroKBN

-- | Add a point to quadreg accumulator.
addQuadReg
    :: QuadRegAcc
    -> Double  -- ^ x
    -> Double  -- ^ y
    -> QuadRegAcc
addQuadReg QuadRegAcc {..} x y = QuadRegAcc
    { qra_n    = qra_n + 1
    , qra_w    = addKBN qra_w   1
    , qra_x    = addKBN qra_x   x
    , qra_x2   = addKBN qra_x2  x2
    , qra_x3   = addKBN qra_x3  (x * x2)
    , qra_x4   = addKBN qra_x4  (x2 * x2)
    , qra_y    = addKBN qra_y   y
    , qra_xy   = addKBN qra_xy  (x * y)
    , qra_x2y  = addKBN qra_x2y (x2 * y)
    , qra_y2   = addKBN qra_y2  (y * y)
    }
  where
    x2 = x * x

-- | Add a weighted point to quadreg accumulator.
addQuadRegW
    :: QuadRegAcc
    -> Double  -- ^ x
    -> Double  -- ^ y
    -> Double  -- ^ w
    -> QuadRegAcc
addQuadRegW QuadRegAcc {..} x y w = QuadRegAcc
    { qra_n    = qra_n + 1
    , qra_w    = addKBN qra_w   w
    , qra_x    = addKBN qra_x   (w * x)
    , qra_x2   = addKBN qra_x2  (w * x2)
    , qra_x3   = addKBN qra_x3  (w * x * x2)
    , qra_x4   = addKBN qra_x4  (w * x2 * x2)
    , qra_y    = addKBN qra_y   (w * y)
    , qra_xy   = addKBN qra_xy  (w * x * y)
    , qra_x2y  = addKBN qra_x2y (w * x2 * y)
    , qra_y2   = addKBN qra_y2  (w * y * y)
    }
  where
    x2 = x * x

-- | Convert 'QuadRegAcc' to 'LinRegAcc'.
--
-- Using this we can try quadratic and linear fits with a single data scan.
--
quadRegAccToLin :: QuadRegAcc -> LinRegAcc
quadRegAccToLin QuadRegAcc {..} = LinRegAcc
    { lra_n  = qra_n
    , lra_w  = qra_w
    , lra_x  = qra_x
    , lra_x2 = qra_x2
    , lra_y  = qra_y
    , lra_xy = qra_xy
    , lra_y2 = qra_y2
    }

-------------------------------------------------------------------------------
-- utils
-------------------------------------------------------------------------------

sq :: Num a => a -> a
sq x = x * x
{-# INLINE sq #-}

iterate1 :: (b -> b) -> b -> NE.NonEmpty b
iterate1 g x = NE.cons x (iterate1 g (g x))

eigenSM22 :: Double -> Double -> Double -> V2
eigenSM22 a b c = V2 ((ac + discr) / 2) ((ac - discr) / 2)
  where
    ac = a + c
    discr = sqrt (sq (a - c) + 4 * sq b)

-- | Levenberg-Marquard stop condition
lmStop :: Double -> Double -> Double -> Bool
lmStop lambda wssr wssr' =
    lambda < 1e-20 || lambda > 1e20 || isNaN wssr' || relDiff < 1e-10
  where
    relDiff = abs (wssr' - wssr) / wssr

max_stars_repo_name: phadej/regression-simple

max_issues_repo_name: phadej/regression-simple

max_forks_repo_name: phadej/regression-simple


text:
import Data.Complex
import Data.List (partition)
import Data.Map ((!))
import qualified Data.Map as M
import Data.Ratio

dft :: [Complex Double] -> [Complex Double]
dft x = matMult dftMat x
  where
    n = length x
    w = exp $ (-2) * pi * (0 :+ 1) / fromIntegral n
    dftMat = [[w ^ (j * k) | j <- [0 .. n - 1]] | k <- [0 .. n - 1]]
    matMult m x = map (sum . zipWith (*) x) m

fft :: [Complex Double] -> [Complex Double]
fft x = fft' x
  where
    n = length x
    w0 = exp ((-2) * pi * (0 :+ 1) / fromIntegral n)
    w = M.fromList [(k % n, w0 ^ k) | k <- [0 .. n - 1]]
    fft' [x] = [x]
    fft' x =
      let (evens, odds) = partition (even . fst) $ zip [0 ..] x
          e = fft' $ map snd evens
          o = fft' $ map snd odds
          x1 = zipWith3 (\e o k -> e + o * w ! (k %n)) e o [0 ..]
          x2 = zipWith3 (\e o k -> e - o * w ! (k %n)) e o [0 ..]
       in x1 ++ x2

main = do
  print $ dft [0, 1, 2, 3]
  print $ fft [0, 1, 2, 3]


text:
-----------------------------------------------------------------------------
-- |
-- Module     : Numeric.LinearAlgebra.Vector.Statistics
-- Copyright  : Copyright (c) 2010, Patrick Perry <patperry@gmail.com>
-- License    : BSD3
-- Maintainer : Patrick Perry <patperry@gmail.com>
-- Stability  : experimental
--
-- Basic multivariate statistics.
--

module Numeric.LinearAlgebra.Vector.Statistics (
    -- * Immutable interface
    sum,
    mean,
    weightedSum,
    weightedMean,

    -- * Mutable interface
    addSumTo,
    meanTo,
    addWeightedSumTo,
    weightedMeanTo,

    ) where

import Prelude hiding ( sum )

import Control.Monad( forM_ )
import Control.Monad.ST( ST )

import Numeric.LinearAlgebra.Types

import Numeric.LinearAlgebra.Vector.Base( Vector )
import Numeric.LinearAlgebra.Vector.STBase( RVector, STVector )
import qualified Numeric.LinearAlgebra.Vector.STBase as V

-- | Returns the sum of the vectors.  The first argument gives the dimension
-- of the vectors.
sum :: (BLAS1 e) => Int -> [Vector e] -> Vector e
sum p xs = V.create $ do
    s <- V.new_ p
    V.clear s
    addSumTo s xs
    return s

-- | Returns the mean of the vectors.  The first argument gives the dimension
-- of the vectors.
mean :: (BLAS1 e) => Int -> [Vector e] -> Vector e
mean p xs = V.create $ do
      m <- V.new_ p
      meanTo m xs
      return m

-- | Returns the weighted sum of the vectors.  The first argument gives the
-- dimension of the vectors.
weightedSum :: (BLAS1 e) => Int -> [(e, Vector e)] -> Vector e
weightedSum p wxs = V.create $ do
    s <- V.new_ p
    V.clear s
    addWeightedSumTo s wxs
    return s

-- | Returns the weighted mean of the vectors.  The first argument gives the
-- dimension of the vectors.
weightedMean :: (BLAS1 e)
             => Int -> [(Double, Vector e)] -> Vector e
weightedMean p wxs = V.create $ do
    m <- V.new_ p
    weightedMeanTo m wxs
    return m

-- | Adds the sum of the vectors to the target vector.
addSumTo :: (RVector v, BLAS1 e) => STVector s e -> [v e] -> ST s ()
addSumTo dst = addWeightedSumTo dst . zip (repeat 1)

-- | Sets the target vector to the mean of the vectors.
meanTo :: (RVector v, BLAS1 e)
       => STVector s e -> [v e] -> ST s()
meanTo dst = weightedMeanTo dst . zip (repeat 1)

-- | Adds the weigthed sum of the vectors to the target vector.
addWeightedSumTo :: (RVector v, BLAS1 e)
                 => STVector s e ->  [(e, v e)] -> ST s ()
addWeightedSumTo s wxs = do
    n <- V.getDim s
    err <- V.new n 0
    old_s <- V.new_ n
    diff <- V.new_ n
    val <- V.new_ n

    forM_ wxs $ \(w,x) -> do
        V.unsafeCopyTo old_s s -- old_s := s

        V.unsafeCopyTo val x   -- val := w * x
        V.scaleM_ w val

        V.addTo err err val    -- err := err + val
        V.addTo s s err        -- s := s + err

        V.subTo diff old_s s   -- diff := old_s - s
        V.addTo err diff val   -- err := diff + val

-- | Sets the target vector to the weighted mean of the vectors.
weightedMeanTo :: (RVector v, BLAS1 e)
               =>  STVector s e -> [(Double, v e)] -> ST s ()
weightedMeanTo m wxs = let
    go _ _ [] = return ()
    go diff w_sum ((w,x):wxs') | w == 0    = go diff w_sum wxs'
                               | otherwise = let w_sum' = w_sum + w
                                             in do
                                    V.subTo diff x m
                                    V.addWithScaleM_
                                        (realToFrac $ w/w_sum') diff m
                                    go diff w_sum' wxs'
    in do
        n <- V.getDim m
        diff <- V.new_ n
        V.clear m
        go diff 0 wxs

max_stars_repo_name: patperry/hs-linear-algebra

max_issues_repo_name: cartazio/hs-cblas

max_forks_repo_name: cartazio/hs-cblas


text:
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeApplications    #-}
module ArrayFire.DataSpec where

import Control.Exception
import Data.Complex
import Data.Word
import Foreign.C.Types
import GHC.Int
import Test.Hspec

import ArrayFire

spec :: Spec
spec =
  describe "Data tests" $ do
    it "Should create constant Array" $ do
      constant @Float [1] 1 `shouldBe` 1
      constant @Double [1] 1 `shouldBe` 1
      constant @Int16 [1] 1 `shouldBe` 1
      constant @Int32 [1] 1 `shouldBe` 1
      constant @Int64 [1] 1 `shouldBe` 1
      constant @Int [1] 1 `shouldBe` 1
      constant @Word16 [1] 1 `shouldBe` 1
      constant @Word32 [1] 1 `shouldBe` 1
      constant @Word64 [1] 1 `shouldBe` 1
      constant @Word [1] 1 `shouldBe` 1
      constant @CBool [1] 1 `shouldBe` 1
      constant @(Complex Double) [1] (1.0 :+ 1.0)
        `shouldBe`
          constant @(Complex Double) [1] (1.0 :+ 1.0)
      constant @(Complex Float) [1] (1.0 :+ 1.0)
        `shouldBe`
          constant @(Complex Float) [1] (1.0 :+ 1.0)

max_stars_repo_name: danielkroeni/arrayfire-haskell

max_issues_repo_name: danielkroeni/arrayfire-haskell

max_forks_repo_name: danielkroeni/arrayfire-haskell


text:
module Test.Grenade.Layers.Internal.Reference where

import           Numeric.LinearAlgebra

im2col :: Int -> Int -> Int -> Int -> Matrix Double -> Matrix Double
im2col nrows ncols srows scols m =
  let starts = fittingStarts (rows m) nrows srows (cols m) ncols scols
  in  im2colFit starts nrows ncols m

vid2col :: Int -> Int -> Int -> Int -> Int -> Int -> [Matrix Double] -> Matrix Double
vid2col nrows ncols srows scols inputrows inputcols ms =
  let starts = fittingStarts inputrows nrows srows inputcols ncols scols
      subs   = fmap (im2colFit starts nrows ncols) ms
  in  foldl1 (|||) subs

im2colFit :: [(Int,Int)] -> Int -> Int -> Matrix Double -> Matrix Double
im2colFit starts nrows ncols m =
  let imRows = fmap (\start -> flatten $ subMatrix start (nrows, ncols) m) starts
  in  fromRows imRows

col2vid :: Int -> Int -> Int -> Int -> Int -> Int -> Matrix Double -> [Matrix Double]
col2vid nrows ncols srows scols drows dcols m =
  let starts = fittingStart (cols m) (nrows * ncols) (nrows * ncols)
      r      = rows m
      mats   = fmap (\s -> subMatrix (0,s) (r, nrows * ncols) m) starts
      colSts = fittingStarts drows nrows srows dcols ncols scols
  in  fmap (col2imfit colSts nrows ncols drows dcols) mats

col2im :: Int -> Int -> Int -> Int -> Int -> Int -> Matrix Double -> Matrix Double
col2im krows kcols srows scols drows dcols m =
  let starts     = fittingStarts drows krows srows dcols kcols scols
  in  col2imfit starts krows kcols drows dcols m

col2imfit :: [(Int,Int)] -> Int -> Int -> Int -> Int -> Matrix Double -> Matrix Double
col2imfit starts krows kcols drows dcols m =
  let indicies   = (\[a,b] -> (a,b)) <$> sequence [[0..(krows-1)], [0..(kcols-1)]]
      convs      = fmap (zip indicies . toList) . toRows $ m
      pairs      = zip convs starts
      accums     = concatMap (\(conv',(stx',sty')) -> fmap (\((ix,iy), val) -> ((ix + stx', iy + sty'), val)) conv') pairs
  in  accum (konst 0 (drows, dcols)) (+) accums

poolForward :: Int -> Int -> Int -> Int -> Int -> Int -> Matrix Double -> Matrix Double
poolForward nrows ncols srows scols outputRows outputCols m =
  let starts = fittingStarts (rows m) nrows srows (cols m) ncols scols
  in  poolForwardFit starts nrows ncols outputRows outputCols m

poolForwardList :: Functor f => Int -> Int -> Int -> Int -> Int -> Int -> Int -> Int -> f (Matrix Double) -> f (Matrix Double)
poolForwardList nrows ncols srows scols inRows inCols outputRows outputCols ms =
  let starts = fittingStarts inRows nrows srows inCols ncols scols
  in  poolForwardFit starts nrows ncols outputRows outputCols <$> ms

poolForwardFit :: [(Int,Int)] -> Int -> Int -> Int -> Int -> Matrix Double -> Matrix Double
poolForwardFit starts nrows ncols _ outputCols m =
  let els    = fmap (\start -> maxElement $ subMatrix start (nrows, ncols) m) starts
  in  matrix outputCols els

poolBackward :: Int -> Int -> Int -> Int -> Matrix Double -> Matrix Double -> Matrix Double
poolBackward krows kcols srows scols inputMatrix gradientMatrix =
  let inRows     = rows inputMatrix
      inCols     = cols inputMatrix
      starts     = fittingStarts inRows krows srows inCols kcols scols
  in  poolBackwardFit starts krows kcols inputMatrix gradientMatrix

poolBackwardList :: Functor f => Int -> Int -> Int -> Int -> Int -> Int -> f (Matrix Double, Matrix Double) -> f (Matrix Double)
poolBackwardList krows kcols srows scols inRows inCols inputMatrices =
  let starts     = fittingStarts inRows krows srows inCols kcols scols
  in  uncurry (poolBackwardFit starts krows kcols) <$> inputMatrices

poolBackwardFit :: [(Int,Int)] -> Int -> Int -> Matrix Double -> Matrix Double -> Matrix Double
poolBackwardFit starts krows kcols inputMatrix gradientMatrix =
  let inRows     = rows inputMatrix
      inCols     = cols inputMatrix
      inds       = fmap (\start -> maxIndex $ subMatrix start (krows, kcols) inputMatrix) starts
      grads      = toList $ flatten gradientMatrix
      grads'     = zip3 starts grads inds
      accums     = fmap (\((stx',sty'),grad,(inx, iny)) -> ((stx' + inx, sty' + iny), grad)) grads'
  in  accum (konst 0 (inRows, inCols)) (+) accums

-- | These functions are not even remotely safe, but it's only called from the statically typed
--   commands, so we should be good ?!?!?
--   Returns the starting sub matrix locations which fit inside the larger matrix for the
--   convolution. Takes into account the stride and kernel size.
fittingStarts :: Int -> Int -> Int -> Int -> Int -> Int -> [(Int,Int)]
fittingStarts nrows kernelrows steprows ncols kernelcols stepcolsh =
  let rs = fittingStart nrows kernelrows steprows
      cs = fittingStart ncols kernelcols stepcolsh
      ls = sequence [rs, cs]
  in  fmap (\[a,b] -> (a,b)) ls

-- | Returns the starting sub vector which fit inside the larger vector for the
--   convolution. Takes into account the stride and kernel size.
fittingStart :: Int -> Int -> Int -> [Int]
fittingStart width kernel steps =
  let go left | left + kernel < width
              = left : go (left + steps)
              | left + kernel == width
              = [left]
              | otherwise
              = []
  in go 0

max_issues_repo_name: Alien-Inc/grenade

max_forks_repo_name: Alien-Inc/grenade


text:
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE Strict           #-}
{-# LANGUAGE StrictData       #-}
module Pinwheel.Transform where

import           Control.DeepSeq
import           Data.Array.IArray          as IA
import           Data.Complex
import           Data.List                  as L
import           Data.Vector.Generic        as VG
import           Data.Vector.Unboxed        as VU
import           FokkerPlanck.FourierSeries (getHarmonics)
import           FokkerPlanck.Histogram
import           Pinwheel.List
import           Sparse.Vector              as SV
import           Utils.List
import           Utils.Parallel             hiding (dot)

data PinwheelTransformData a b =
  PinwheelTransformData a
                        a
                        b

instance (NFData a, NFData b) => NFData (PinwheelTransformData a b) where
  rnf (PinwheelTransformData x y z) = x `seq` y `seq` z `seq` ()

{-# INLINE projectVec #-}
projectVec ::
     ( VG.Vector vector e
     , VG.Vector vector (Complex e)
     , VG.Vector vector Int
     , Num e
     , RealFloat e
     )
  => Int
  -> Int
  -> e
  -> vector (Complex e)
  -> PinwheelTransformData e (SparseVector vector (Complex e))
  -> Complex e
projectVec angularFreq radialFreq sigma pinwheel (PinwheelTransformData logR theta sparseVec) =
  ((exp ((sigma - 1) * logR)) :+ 0) *
  (cis $
   (-1) * (theta * fromIntegral angularFreq + logR * fromIntegral radialFreq)) *
  (sparseVec `dot` pinwheel)

{-# INLINE pinwheelTransform #-}
pinwheelTransform ::
     ( VG.Vector vector e
     , VG.Vector vector (Complex e)
     , VG.Vector vector Int
     , NFData e
     , RealFloat e
     , Unbox e
     )
  => Int
  -> Int
  -> Int
  -> Int
  -> Int
  -> e
  -> IA.Array (Int, Int) (vector (Complex e))
  -> [PinwheelTransformData e (SparseVector vector (Complex e))]
  -> Histogram (Complex e)
pinwheelTransform numThread maxPhiFreq maxRhoFreq maxThetaFreq maxRFreq sigma pinwheelArr xs =
  let coef =
        parMapChunk
          (ParallelParams numThread undefined)
          rdeepseq
          (\(rFreq, thetaFreq, rhoFreq, phiFreq) ->
             let pinwheel =
                   getHarmonics
                     pinwheelArr
                     (fromIntegral phiFreq)
                     (fromIntegral rhoFreq)
                     (fromIntegral thetaFreq)
                     (fromIntegral rFreq)
             in L.foldl'
                  (\s x -> s + projectVec thetaFreq rFreq sigma pinwheel x)
                  0
                  xs) $
        [ (rFreq, thetaFreq, rhoFreq, phiFreq)
        | rFreq <- [-maxRFreq .. maxRFreq]
        , thetaFreq <- [-maxThetaFreq .. maxThetaFreq]
        , rhoFreq <- [-maxRhoFreq .. maxRhoFreq]
        , phiFreq <- [-maxPhiFreq .. maxPhiFreq]
        ]
  in Histogram
       [ 2 * maxPhiFreq + 1
       , 2 * maxRhoFreq + 1
       , 2 * maxThetaFreq + 1
       , 2 * maxRFreq + 1
       ]
       1 .
     VU.fromList $
     coef

max_stars_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_issues_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_forks_repo_name: XinhuaZhang/Stochastic-Completion-Field


text:
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE ViewPatterns     #-}
-- Tests for Statistics.Test.NonParametric
module Tests.NonParametric (tests) where

import Statistics.Distribution.Normal (standard)
import Statistics.Test.KolmogorovSmirnov
import Statistics.Test.MannWhitneyU
import Statistics.Test.KruskalWallis
import Statistics.Test.WilcoxonT
import Statistics.Types (PValue,pValue,mkPValue)

import Test.Tasty                (testGroup)
import Test.Tasty.HUnit
import Tests.ApproxEq            (eq)
import Tests.Helpers             (testAssertion, testEquality)
import Tests.NonParametric.Table (tableKSD, tableKS2D)
import qualified Test.Tasty          as Tst
import qualified Data.Vector.Unboxed as U

tests :: Tst.TestTree
tests = testGroup "Nonparametric tests"
        $ concat [ mannWhitneyTests
                 , wilcoxonSumTests
                 , wilcoxonPairTests
                 , kruskalWallisRankTests
                 , kruskalWallisTests
                 , kolmogorovSmirnovDTest
                 ]

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

mannWhitneyTests :: [Tst.TestTree]
mannWhitneyTests = zipWith test [(0::Int)..] testData ++
  [ testEquality "Mann-Whitney U Critical Values, m=1"
      (replicate (20*3) Nothing)
      [mannWhitneyUCriticalValue (1,x) (mkPValue p) | x <- [1..20], p <- [0.005,0.01,0.025]]
  , testEquality "Mann-Whitney U Critical Values, m=2, p=0.025"
      (replicate 7 Nothing ++ map Just [0,0,0,0,1,1,1,1,1,2,2,2,2])
      [mannWhitneyUCriticalValue (2,x) (mkPValue 0.025) | x <- [1..20]]
  , testEquality "Mann-Whitney U Critical Values, m=6, p=0.05"
      (replicate 1 Nothing ++ map Just [0, 2,3,5,7,8,10,12,14,16,17,19,21,23,25,26,28,30,32])
      [mannWhitneyUCriticalValue (6,x) (mkPValue 0.05) | x <- [1..20]]
  , testEquality "Mann-Whitney U Critical Values, m=20, p=0.025"
      (replicate 1 Nothing ++ map Just [2,8,14,20,27,34,41,48,55,62,69,76,83,90,98,105,112,119,127])
      [mannWhitneyUCriticalValue (20,x) (mkPValue 0.025) | x <- [1..20]]
  ]
  where
    test n (a, b, c, d)
      = testCase "Mann-Whitney" $ do
          assertEqual ("Mann-Whitney U "     ++ show n) c us
          assertEqual ("Mann-Whitney U Sig " ++ show n) d ss
      where
        us = mannWhitneyU (U.fromList a) (U.fromList b)
        ss = mannWhitneyUSignificant SamplesDiffer (length a, length b) p005 us
    -- List of (Sample A, Sample B, (Positive Rank, Negative Rank))
    testData :: [([Double], [Double], (Double, Double), Maybe TestResult)]
    testData = [ ( [3,4,2,6,2,5]
                 , [9,7,5,10,6,8]
                 , (2, 34)
                 , Just Significant
                 )
               , ( [540,480,600,590,605]
                 , [760,890,1105,595,940]
                 , (2, 23)
                 , Just Significant
                 )
               , ( [19,22,16,29,24]
                 , [20,11,17,12]
                 , (17, 3)
                 , Just NotSignificant
                 )
               , ( [126,148,85,61, 179,93, 45,189,85,93]
                 , [194,128,69,135,171,149,89,248,79,137]
                 , (35,65)
                 , Just NotSignificant
                 )
               , ( [1..30]
                 , [1..30]
                 , (450,450)
                 , Just NotSignificant
                 )
               , ( [1 .. 30]
                 , [11.5 .. 40 ]
                 , (190.0,710.0)
                 , Just Significant
                 )
               ]

wilcoxonSumTests :: [Tst.TestTree]
wilcoxonSumTests = zipWith test [(0::Int)..] testData
  where
    test n (a, b, c) = testCase "Wilcoxon Sum"
                     $ assertEqual ("Wilcoxon Sum " ++ show n) c (wilcoxonRankSums (U.fromList a) (U.fromList b))
    -- List of (Sample A, Sample B, (Positive Rank, Negative Rank))
    testData :: [([Double], [Double], (Double, Double))]
    testData = [ ( [8.50,9.48,8.65,8.16,8.83,7.76,8.63]
                 , [8.27,8.20,8.25,8.14,9.00,8.10,7.20,8.32,7.70]
                 , (75, 61)
                 )
               , ( [0.45,0.50,0.61,0.63,0.75,0.85,0.93]
                 , [0.44,0.45,0.52,0.53,0.56,0.58,0.58,0.65,0.79]
                 , (71.5, 64.5)
                 )
               ]

wilcoxonPairTests :: [Tst.TestTree]
wilcoxonPairTests = zipWith test [(0::Int)..] testData ++
  -- Taken from the Mitic paper:
  [ testAssertion "Sig 16, 35" (to4dp 0.0467 $ wilcoxonMatchedPairSignificance 16 35)
  , testAssertion "Sig 16, 36" (to4dp 0.0523 $ wilcoxonMatchedPairSignificance 16 36)
  , testEquality   "Wilcoxon critical values, p=0.05"
      (replicate 4 Nothing ++ map Just [0,2,3,5,8,10,13,17,21,25,30,35,41,47,53,60,67,75,83,91,100,110,119])
      [wilcoxonMatchedPairCriticalValue x (mkPValue 0.05) | x <- [1..27]]
  , testEquality "Wilcoxon critical values, p=0.025"
      (replicate 5 Nothing ++ map Just [0,2,3,5,8,10,13,17,21,25,29,34,40,46,52,58,65,73,81,89,98,107])
      [wilcoxonMatchedPairCriticalValue x (mkPValue 0.025) | x <- [1..27]]
  , testEquality "Wilcoxon critical values, p=0.01"
      (replicate 6 Nothing ++ map Just [0,1,3,5,7,9,12,15,19,23,27,32,37,43,49,55,62,69,76,84,92])
      [wilcoxonMatchedPairCriticalValue x (mkPValue 0.01) | x <- [1..27]]
  , testEquality "Wilcoxon critical values, p=0.005"
      (replicate 7 Nothing ++ map Just [0,1,3,5,7,9,12,15,19,23,27,32,37,42,48,54,61,68,75,83])
      [wilcoxonMatchedPairCriticalValue x (mkPValue 0.005) | x <- [1..27]]
  ]
  where
    test n (a, b, c) = testEquality ("Wilcoxon Paired " ++ show n) c res
      where res = wilcoxonMatchedPairSignedRank (U.zip (U.fromList a) (U.fromList b))

    -- List of (Sample A, Sample B, (Positive Rank, Negative Rank))
    testData :: [([Double], [Double], (Int,Double, Double))]
    testData = [ ([1..10], [1..10], (0, 0, 0     ))
               , ([1..5],  [6..10], (5, 0, 5*(-3)))
               -- Worked example from the Internet:
               , ( [125,115,130,140,140,115,140,125,140,135]
                 , [110,122,125,120,140,124,123,137,135,145]
                 , ( 9
                   , sum $ filter (> 0) [7,-3,1.5,9,0,-4,8,-6,1.5,-5]
                   , sum $ filter (< 0) [7,-3,1.5,9,0,-4,8,-6,1.5,-5]
                   )
                 )
               -- Worked examples from books/papers:
               , ( [2.4,1.9,2.3,1.9,2.4,2.5]
                 , [2.0,2.1,2.0,2.0,1.8,2.0]
                 , (6, 18, -3)
                 )
               , ( [130,170,125,170,130,130,145,160]
                 , [120,163,120,135,143,136,144,120]
                 , (8, 27, -9)
                 )
               , ( [540,580,600,680,430,740,600,690,605,520]
                 , [760,710,1105,880,500,990,1050,640,595,520]
                 , (9, 3, -42)
                 )
               ]
    to4dp tgt (pValue -> x) = x >= tgt - 0.00005 && x < tgt + 0.00005

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

kruskalWallisRankTests :: [Tst.TestTree]
kruskalWallisRankTests = zipWith test [(0::Int)..] testData
  where
    test n (a, b) = testCase "Kruskal-Wallis Ranking"
                  $ assertEqual ("Kruskal-Wallis " ++ show n) (map U.fromList b) (kruskalWallisRank $ map U.fromList a)
    testData :: [([[Int]],[[Double]])]
    testData = [ ( [ [68,93,123,83,108,122]
                   , [119,116,101,103,113,84]
                   , [70,68,54,73,81,68]
                   , [61,54,59,67,59,70]
                   ]
                 , [ [8.0,14.0,16.0,19.0,23.0,24.0]
                   , [15.0,17.0,18.0,20.0,21.0,22.0]
                   , [1.5,8.0,8.0,10.5,12.0,13.0]
                   , [1.5,3.5,3.5,5.0,6.0,10.5]
                   ]
                 )
               ]

kruskalWallisTests :: [Tst.TestTree]
kruskalWallisTests = zipWith test [(0::Int)..] testData
  where
    test n (a, b, c) = testCase "Kruskal-Wallis" $ do
        assertEqual ("Kruskal-Wallis " ++ show n) (round100 b) (round100 kw)
        assertEqual ("Kruskal-Wallis Sig " ++ show n) c kwt
      where
        kw  = kruskalWallis $ map U.fromList a
        kwt = isSignificant p005 `fmap` kruskalWallisTest (map U.fromList a)
        round100 :: Double -> Integer
        round100 = round . (*100)

    testData :: [([[Double]], Double, Maybe TestResult)]
    testData = [ ( [ [68,93,123,83,108,122]
                   , [119,116,101,103,113,84]
                   , [70,68,54,73,81,68]
                   , [61,54,59,67,59,70]
                   ]
                 , 16.03
                 , Just Significant
                 )
               , ( [ [5,5,3,5,5,5,5]
                   , [5,5,5,5,7,5,5]
                   , [5,5,6,5,5,5,5]
                   , [4,5,5,5,6,5,5]
                   ]
               , 2.24
               , Just NotSignificant
               )
               , ( [ [36,48,5,67,53]
                   , [49,33,60,2,55]
                   , [71,31,140,59,42]
                   ]
                 , 1.22
                 , Just NotSignificant
                 )
               , ( [ [6,38,3,17,11,30,15,16,25,5]
                   , [34,28,42,13,40,31,9,32,39,27]
                   , [13,35,19,4,29,0,7,33,18,24]
                   ]
                 , 6.10
                 , Just Significant
                 )
               ]

----------------------------------------------------------------
-- K-S test
----------------------------------------------------------------

kolmogorovSmirnovDTest :: [Tst.TestTree]
kolmogorovSmirnovDTest =
  [ testAssertion "K-S D statistics" $
    and [ eq 1e-6 (kolmogorovSmirnovD standard (toU sample)) reference
        | (reference,sample) <- tableKSD
        ]
  , testAssertion "K-S 2-sample statistics" $
    and [ eq 1e-6 (kolmogorovSmirnov2D (toU xs) (toU ys)) reference
        | (reference,xs,ys) <- tableKS2D
        ]
  , testAssertion "K-S probability" $
    and [ eq 1e-14 (kolmogorovSmirnovProbability n d) p
        | (d,n,p) <- testData
        ]
  ]
  where
    toU = U.fromList
    -- Test data for the calculation of cumulative probability
    -- P(D[n] < d).
    --
    -- Test data is:
    --    (D[n], n, p)
    -- Table is generated using sample program from paper
    testData :: [(Double,Int,Double)]
    testData =
      [ (0.09           ,    3, 0                   )
      , (0.2            ,    3, 0.00177777777777778 )
      , (0.301          ,    3, 0.116357025777778   )
      , (0.392          ,    3, 0.383127210666667   )
      , (0.5003         ,    3, 0.667366306558667   )
      , (0.604          ,    3, 0.861569877333333   )
      , (0.699          ,    3, 0.945458198         )
      , (0.802          ,    3, 0.984475216         )
      , (0.9            ,    3, 0.998               )
      , (0.09           ,    5, 0                   )
      , (0.2            ,    5, 0.0384              )
      , (0.301          ,    5, 0.33993786080016    )
      , (0.392          ,    5, 0.66931908083712    )
      , (0.5003         ,    5, 0.888397260183794   )
      , (0.604          ,    5, 0.971609957879808   )
      , (0.699          ,    5, 0.994331075994008   )
      , (0.802          ,    5, 0.999391366368064   )
      , (0.9            ,    5, 0.99998             )
      , (0.09           ,    8, 3.37615237575e-06   )
      , (0.2            ,    8, 0.151622071801758   )
      , (0.301          ,    8, 0.613891042670582   )
      , (0.392          ,    8, 0.871491561427005   )
      , (0.5003         ,    8, 0.977534089199071   )
      , (0.604          ,    8, 0.997473116268255   )
      , (0.699          ,    8, 0.999806082005123   )
      , (0.802          ,    8, 0.999995133786947   )
      , (0.9            ,    8, 0.99999998          )
      , (0.09           ,   10, 3.89639433093119e-05)
      , (0.2            ,   10, 0.25128096          )
      , (0.301          ,   10, 0.732913126355935   )
      , (0.392          ,   10, 0.932185254518767   )
      , (0.5003         ,   10, 0.992276179340446   )
      , (0.604          ,   10, 0.999495533516769   )
      , (0.699          ,   10, 0.999979691783985   )
      , (0.802          ,   10, 0.999999801409237   )
      , (0.09           ,   20, 0.00794502217168886 )
      , (0.2            ,   20, 0.647279826376584   )
      , (0.301          ,   20, 0.958017466965765   )
      , (0.392          ,   20, 0.997206424843499   )
      , (0.5003         ,   20, 0.999962641414228   )
      , (0.09           ,   30, 0.0498147538075168  )
      , (0.2            ,   30, 0.842030838984526   )
      , (0.301          ,   30, 0.993403560017612   )
      , (0.392          ,   30, 0.99988478803318    )
      , (0.09           ,  100, 0.629367974413669   )
      ]

p005 :: PValue Double
p005 = mkPValue 0.05


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

module AI.DemoLayer
( LayerDefinition(..)
, Layer(..)

, Weights
, Biases
, LayerWidth
, Connectivity

--, createLayer
) where

import AI.DemoNeuron

import Numeric.LinearAlgebra
import Numeric.LinearAlgebra.Data
import Data.Binary
import Data.Binary.Put
import qualified Data.ByteString.Lazy as B
import System.Random

data LayerDefinition a = LayerDefinition { neuronType   :: a
                                         , width        :: LayerWidth
                                         , connectivity :: Connectivity
                                         }

data Layer a = Layer { neuron :: a
                     , weights :: Weights
                     , biases :: Biases
                     } deriving (Show)

type Weights = Matrix Double
type Biases = Vector Double
type LayerWidth = Int
type Connectivity = LayerWidth -> LayerWidth -> Weights

instance (Neuron a) => Binary (Layer a) where
  put Layer{..} = do put weights; put biases
  get = do weights <- get; biases <- get; return Layer{..}

-- createLayer :: (Neuron a, RandomGen g)
--               => RandomTransform -> g -> LayerDefinition a
--               -> LayerDefinition a -> Layer
-- createLayer t g li lj = Layer (neuron lj)

max_stars_repo_name: jbarrow/LambdaNet

max_issues_repo_name: jbarrow/LambdaNet

max_forks_repo_name: jbarrow/LambdaNet


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

module Mandelbrot
    ( mandelbrot
    ) where

import Data.Complex
import Debug.Trace

mag2 :: RealFloat a => Complex a -> a
mag2 z = (realPart z)^2 + (imagPart z)^2
{-# INLINE mag2 #-}

-- For given point and maxIterations returns iteration number when diverged
mandelbrot :: forall a. RealFloat a => Int -> Complex a -> Int
mandelbrot maxIter c = mandelbrot_ 0 0
    where
        mandelbrot_ :: Int -> Complex a -> Int
        mandelbrot_ !iter !z0 =
            if iter >= maxIter || (mag2 z0) > 4
            then iter
            else mandelbrot_ (iter+1) (z0*z0 + c)
        {-# INLINE mandelbrot_ #-}

{-# SPECIALIZE mandelbrot :: Int -> Complex Double -> Int #-}

max_stars_repo_name: hator/handelbrot

max_issues_repo_name: hator/handelbrot

max_forks_repo_name: hator/handelbrot


text:
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE TupleSections #-}
{-# LANGUAGE ScopedTypeVariables #-}
module Day12 where

import Control.Lens
import Control.Lens.TH

import qualified Data.ByteString as BS
import Data.Attoparsec.ByteString (Parser)
import qualified Data.Attoparsec.ByteString as DAB
import qualified Data.Attoparsec.ByteString.Char8 as DABC8
import Data.Attoparsec.ByteString.Char8 (decimal, inClass, digit, endOfLine, char, anyChar, letter_ascii, notChar)

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

import qualified Data.List as DL
import qualified Data.List.Split as DLS

import Debug.Trace

import Control.Arrow

import Data.Array.IArray (Array)
import qualified Data.Array.IArray as IA

import Data.Set (Set)
import qualified Data.Set as Set
import Data.Either (isRight)
import Data.Maybe (mapMaybe, catMaybes)
import Data.Bool (bool)

import Data.Vector (Vector, (!))
import qualified Data.Vector as Vec

import Data.Function (on)

import Control.Monad.State.Strict (State)
import qualified Control.Monad.State.Strict as State

import qualified Control.Monad as CM

-- import qualified Control.Foldl as L

import Data.Bifunctor (bimap)

import Data.Complex

--type Input = (IA.Array Int (IA.Array Int Char))
type Input = (Char, Int)

-- >>> solver [('F', 10),

solver input =
	uncurry (+) $
	(abs . realPart &&& abs . imagPart) $ snd $
	DL.foldl' (\x y -> traceShow (x, y) (navigate y x)) (10 :+ 1, 0 :+ 0) input
	where
	navigate :: Input -> (Complex Float, Complex Float) -> (Complex Float, Complex Float)
	navigate ('N', n) = first (+ imag n)
	navigate ('S', n) = first (dec (imag n))
	navigate ('E', n) = first (+ real n)
	navigate ('W', n) = first (dec (real n))

	navigate ('L', 90) = first (* (0 :+ 1))
	navigate ('L', 180) = first $ fmap round' . (* (0 :+ 1) ** 2)
	navigate ('L', 270) = first $ fmap round' . (* (0 :+ 1) ** 3)

	navigate ('R', 90) = first $ fmap round' . (* (0 :+ (-1)))
	navigate ('R', 180) = first $ fmap round' . (* (0 :+ (-1)) ** 2)
	navigate ('R', 270) = first $ fmap round' . (* (0 :+ (-1)) ** 3)
	navigate ('F', n) = uncurry (CM.liftM2 (.) (,) ((+) . (real n *)))

	navigate (c, n) = error $ show (c, n)

	norm :: Complex Float -> Complex Float
	norm x = x / abs x

	real n = fromIntegral n :+ 0
	imag n = 0 :+ fromIntegral n
	dec n x = x - n

	round' :: Float -> Float
	round' = fromIntegral . round

parseInput :: Parser [Input]
parseInput =
	((,) <$> anyChar <*> decimal) `DAB.sepBy1'` endOfLine

runSolution :: FilePath -> IO ()
runSolution filePath = do
	contents <- BS.readFile filePath
	let parseResult = DAB.parseOnly parseInput contents
	case parseResult of
		Left err -> putStrLn err
		Right input -> do
			print $ solver input

-- 1526 wrong
-- part ii
-- 1463 wrong


text:
module Network.Mux.DeltaQ.TraceStatsSupport where

-- This module is making use of hackage statisitical libraries. They
-- are not the most efficicient approches for this particlular use
-- case, and they may increase the package dependencies for the final
-- binaries (they have a lot of dependencies).
--
-- It may well be worthwhile constructing specialsed version for the
-- specific use case, but building those and creating the associated
-- test suite was not deemed a good use of time (at the time of
-- creation).
--
-- Definite space/time optimisation task here.

import           Network.Mux.DeltaQ.TraceTypes

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

estimateGS :: [(Int, SISec)] -> (Double, Double, Double)
estimateGS xys
  = let (xs', ys') = unzip xys
        xs = V.fromList $ map fromIntegral xs'
        ys = V.fromList $ map (\(S x) -> fromRational . toRational $ x) ys'
    in linearRegressionRSqr xs ys

max_stars_repo_name: RyanGlScott/ouroboros-network

max_issues_repo_name: RyanGlScott/ouroboros-network

max_forks_repo_name: RyanGlScott/ouroboros-network


text:
import Graphics.UI.GLUT
import Data.Complex

iterations = 400

chunk :: Int -> [a] -> [[a]]
chunk n [] = []
chunk n xs = let (xs', rest) = splitAt n xs in xs' : chunk n rest

vertices :: IO [Vertex2 GLfloat]
vertices = get windowSize >>= \(Size w h) ->
           return $ [pixel (Size w h) (x, y) | x <- [0..w-1], y <- [0..h-1]]
  where
    pixel (Size w h) (x, y) = Vertex2 ((fromIntegral 3 / fromIntegral w * fromIntegral x) - 2.0) ((fromIntegral 2 / fromIntegral h * fromIntegral y) - 1.0)

getcol :: Int -> Color3 Float
getcol iter | iter == iterations = Color3 0 0 0
              | otherwise          = Color3 (amt*0.5) amt (amt*0.5)
  where
    amt = fromIntegral iter / fromIntegral iterations

mandelbrot :: RealFloat a => Vertex2 a -> Int
mandelbrot (Vertex2 r i) = length . takeWhile (\z -> magnitude z <= 2) . take iterations $ iterate (\z -> z^2 + (r :+ i)) 0

drawVert :: (RealFloat a, VertexComponent a) => Vertex2 a -> IO ()
drawVert v = do color . getcol $ mandelbrot v
                vertex v

display' :: [[Vertex2 GLfloat]] -> IO ()
display' chunks = do mapM_ (\vs -> do renderPrimitive Points $ mapM_ drawVert vs
                                      flush) chunks
                     displayCallback $= display

display :: IO ()
display = do clear [ ColorBuffer ]
             displayCallback $= (vertices >>= display' . chunk 256)
             get currentWindow >>= postRedisplay

main :: IO ()
main = do
   getArgsAndInitialize
   initialDisplayMode $= [ SingleBuffered, RGBMode ]
   initialWindowSize $= Size 1200 1024
   initialWindowPosition $= Position 100 100
   createWindow "Mandelbrot"
   clearColor $= Color4 0 0 0 0
   matrixMode $= Projection
   loadIdentity
   ortho (-2) 1 (-1) 1 (-1) 1
   displayCallback $= display
   mainLoop

max_stars_repo_name: atlas-oss/mandelbrot

max_issues_repo_name: atlas-oss/mandelbrot

max_forks_repo_name: atlas-oss/mandelbrot


text:
module Main where
import System.Environment (getArgs)
import Data.Complex
import Data.List.Split (chunksOf)
import Text.Read (readMaybe)
import Data.Maybe (isJust)

type Point = Complex Double
data MandelbrotPoint = MandelbrotPoint Int {-n-} Point {-z_n-} Point {-c-}

instance Show MandelbrotPoint where
  show (MandelbrotPoint n _ _)
    | n == 100  = " "
    | n > 10    = "#"
    | otherwise = [".',;\"oO%8@" !! (n - 1)]

data Mandelbrot = Mandelbrot [MandelbrotPoint] {-state-} Int {-width-}

instance Show Mandelbrot where
  show (Mandelbrot state width) =
    unlines $ chunksOf width $ concatMap show state

data Arguments = Arguments Point {-topleft-} Point {-bottomright-} Int {-width-} Int {-height-}

main :: IO ()
main = do
  args <- getArgs
  putStr $ show $ mandelbrot 100 $ initialise $ parseArguments args

usage :: IO ()
usage = do
  putStrLn "Usage: mandelbrot [TOPLEFT BOTTOMRIGHT WIDTH HEIGHT]"
  putStrLn "  TOPLEFT      Top-left corner of complex subplane"
  putStrLn "  BOTTOMRIGHT  Bottom-right corner of complex subplane"
  putStrLn "  WIDTH        Output width (characters)"
  putStrLn "  HEIGHT       Output height (characters)"

parseArguments :: [String] -> Arguments
parseArguments [topLeft, bottomRight, width, height] =
  if (all isJust planeCoords) && (all isJust imageCoords)
    -- TODO Sanity check values
    then Arguments (read topLeft) (read bottomRight) (read width) (read height)
    else error "Bad arguments" -- TODO: Print usage and exit 1
  where planeCoords = map readMaybe [topLeft, bottomRight] :: [Maybe Point]
        imageCoords = map readMaybe [width, height]        :: [Maybe Int]

-- Default values
parseArguments _ = Arguments ((-2) :+ 1) (1 :+ (-1)) 99 33

initialise :: Arguments -> Mandelbrot
initialise (Arguments topLeft bottomRight width height) =
  Mandelbrot (map (MandelbrotPoint 0 0) field) width
  where [w, h]    = map fromIntegral [width, height]
        diagonal  = bottomRight - topLeft
        scaleReal = realPart diagonal / (w - 1)
        scaleImag = imagPart diagonal / (h - 1)
        field     = [ ((x * scaleReal) :+ (y * scaleImag)) + topLeft
                    | y <- [0..(h - 1)], x <- [0..(w - 1)]]

mandelbrot :: Int -> Mandelbrot -> Mandelbrot
mandelbrot 0 state = state
mandelbrot x (Mandelbrot state width) =
  mandelbrot (x - 1) (Mandelbrot (map step state) width)
    where step point@(MandelbrotPoint n z_n c)
            | magnitude z_n > 2 = point
            | otherwise         = MandelbrotPoint (n + 1) ((z_n * z_n) + c) c

max_stars_repo_name: Xophmeister/mandelbrot

max_issues_repo_name: Xophmeister/mandelbrot

max_forks_repo_name: Xophmeister/mandelbrot


text:
{-# LANGUAGE FlexibleInstances #-}

module TheField.Plot where

import Data.Complex
import Data.List as L
import Graphics.Rendering.Chart.Easy

class Plotable a where
  plotable :: RealFloat b => a -> (b, b)

instance RealFloat a => Plotable (Complex a) where
  plotable c = (realToFrac . realPart $ c, realToFrac . imagPart $ c)

instance RealFloat a => Plotable (a, a) where
  plotable (a, b) = (realToFrac a, realToFrac b)

plot' :: (Plotable a, RealFloat b, Show b) => [a] -> b -> b -> EC (Layout b b) ()
plot' s w h = do
  setColors [opaque red]
  plot (points "s" $ L.map plotable s)
  vectorPlot w h

plotComplexAddition :: Complex Double -> Complex Double -> EC (Layout Double Double) ()
plotComplexAddition a b = do
  vectorPlot (maxReal * 1.1) (maxImag * 1.1)
  plot . liftEC $ do
    plot_lines_style . line_color .= opaque red
    plot_lines_style . line_width .= 3
    plot_lines_values .= L.map (L.map plotable) [[0, a], [a, a + b]]
  plot . liftEC $ do
    plot_lines_style . line_color .= opaque green
    plot_lines_style . line_width .= 2
    plot_lines_style . line_dashes .= [5, 5]
    plot_lines_values .= [ L.map plotable [0, a + b] ]
  where
    maxReal = maximum [ realPart a, realPart b, realPart (a + b) ]
    maxImag = maximum [ imagPart a, imagPart b, imagPart (a + b) ]

vectorPlot :: (RealFloat a, Show a) => a -> a -> EC (Layout a a) ()
vectorPlot w h = do
  layout_x_axis . laxis_generate .= scaledAxis def (-w,w)
  layout_y_axis . laxis_generate .= scaledAxis def (-h,h)
  layout_bottom_axis_visibility . axis_show_line .= False
  layout_bottom_axis_visibility . axis_show_ticks .= False
  layout_left_axis_visibility . axis_show_line .= False
  layout_left_axis_visibility . axis_show_ticks .= False
  layout_legend .= Nothing
  setColors [opaque black, opaque black]
  plot (line "x" [x])
  plot (line "y" [y])
    where
      x = [(-w*2, 0), (w*2, 0)]
      y = [(0, -h*2), (0, h*2)]

max_stars_repo_name: danbroooks/coding-the-matrix

max_issues_repo_name: danbroooks/coding-the-matrix

max_forks_repo_name: danbroooks/coding-the-matrix