text:
{-# LANGUAGE BangPatterns #-}
module STCFourierSeries where

import           Control.DeepSeq
import           Control.Monad                      as M
import qualified Data.Array.Accelerate              as A
import qualified Data.Array.Accelerate.Data.Complex as A
import           Data.Array.Accelerate.LLVM.PTX
import           Data.Array.IArray
import           Data.Array.Repa                    as R
import           Data.Binary
import           Data.Complex
import           Data.List                          as L
import           Data.Maybe
import           Data.Vector.Generic                as VG
import           FokkerPlanck.BrownianMotion
import           FokkerPlanck.FourierSeries
import           FokkerPlanck.GreensFunction
import           FokkerPlanck.Histogram
import           FokkerPlanck.MonteCarlo
import           Foreign.CUDA.Driver                as CUDA
import           Image.IO
import           STC
import           System.Directory
import           System.Environment
import           System.FilePath
import           Text.Printf
import           Utils.Array
import           Utils.Time

main = do
  args@(gpuIDStr:numPointStr:deltaXStr:numOrientationStr:numScaleStr:thetaSigmaStr:scaleSigmaStr:maxScaleStr:deltaLogStr:taoStr:numTrailStr:maxTrailStr:phiFreqsStr:rhoFreqsStr:thetaFreqsStr:scaleFreqsStr:initDistStr:initScaleStr:histFilePath:stdStr:stdGStr:locationStr:numThreadStr:_) <-
    getArgs
  let gpuID = read gpuIDStr :: [Int]
      numPoint = read numPointStr :: Int
      deltaX = read deltaXStr :: Double
      numOrientation = read numOrientationStr :: Int
      numScale = read numScaleStr :: Int
      thetaSigma = read thetaSigmaStr :: Double
      scaleSigma = read scaleSigmaStr :: Double
      tao = read taoStr :: Double
      numTrail = read numTrailStr :: Int
      maxTrail = read maxTrailStr :: Int
      phiFreq = read phiFreqsStr :: Double
      phiFreqs = [-phiFreq .. phiFreq]
      rhoFreq = read rhoFreqsStr :: Double
      rhoFreqs = [-rhoFreq .. rhoFreq]
      thetaFreq = read thetaFreqsStr :: Double
      thetaFreqs = [-thetaFreq .. thetaFreq]
      scaleFreq = read scaleFreqsStr :: Double
      scaleFreqs = [-scaleFreq .. scaleFreq]
      initDist = read initDistStr :: [(Double, Double, Double, Double)]
      initScale = read initScaleStr :: Double
      initPoints = L.map (\(x, y, t, s) -> Point x y t s) initDist
      numThread = read numThreadStr :: Int
      folderPath = "output/test/STCFourierSeries"
      maxScale = read maxScaleStr :: Double
      halfLogPeriod = log maxScale
      deltaLog = read deltaLogStr :: Double
      std = read stdStr :: Double
      stdG = read stdGStr :: Double
      location = read locationStr :: (Int, Int)
  removePathForcibly folderPath
  createDirectoryIfMissing True folderPath
  flag <- doesFileExist histFilePath
  -- hist <-
  --   if flag
  --     then do
  --       printCurrentTime $ "read data from " L.++ histFilePath
  --       decodeFile histFilePath
  --     else runMonteCarloFourierCoefficientsGPU
  --            gpuID
  --            numThread
  --            numTrail
  --            maxTrail
  --            thetaSigma
  --            scaleSigma
  --            maxScale
  --            tao
  --            phiFreqs
  --            rhoFreqs
  --            thetaFreqs
  --            scaleFreqs
  --            deltaLog
  --            initScale
  --            histFilePath
  --            (emptyHistogram
  --               [ L.length phiFreqs
  --               , L.length rhoFreqs
  --               , L.length thetaFreqs
  --               , L.length scaleFreqs
  --               ]
  --               0)
  -- hist <-
  --   sampleLogpolar
  --     histFilePath
  --     (L.head gpuID)
  --     180
  --     180
  --     (fromIntegral $ div numPoint 2)
  --     initScale
  --     thetaSigma
  --     tao
  --     phiFreqs
  --     rhoFreqs
  --     thetaFreqs
  --     scaleFreqs
  hist <-
    sampleCartesian
      folderPath
      histFilePath
      gpuID
      (fromIntegral $ div numPoint 2)
      144
      deltaLog
      initScale
      thetaSigma
      tao
      phiFreqs
      rhoFreqs
      thetaFreqs
      scaleFreqs
  let !initSource =
        computeInitialDistribution'
          numPoint
          numPoint
          phiFreqs
          rhoFreqs
          -- thetaFreqs
          halfLogPeriod
          [L.head initPoints]
      !initSink =
        computeInitialDistribution'
          numPoint
          numPoint
          phiFreqs
          rhoFreqs
          -- thetaFreqs
          halfLogPeriod
          [L.last initPoints]
      !coefficients =
        normalizeFreqArr' std phiFreqs rhoFreqs . getNormalizedHistogramArr $
        hist
       -- = getNormalizedHistogramArr $ hist :: R.Array U DIM4 (Complex Double)
      !thetaRHarmonics =
        computeThetaRHarmonics
          numOrientation
          numScale
          thetaFreqs
          scaleFreqs
          halfLogPeriod
  plan <-
    makePlan
      folderPath
      emptyPlan
      numPoint
      numPoint
      (L.length thetaFreqs)
      (L.length scaleFreqs)
  printCurrentTime "harmonicsArray"
  -- gaussian <- gaussianFilter2D plan numPoint stdG
  harmonicsArray <-
    dftHarmonicsArray
      plan
      numPoint
      deltaX
      numPoint
      deltaX
      phiFreqs
      rhoFreqs
      thetaFreqs
      scaleFreqs
      halfLogPeriod
      (fromIntegral numPoint * sqrt 2)
      -- gaussian
  -- let harmonicsArray' =
  --       computeHarmonicsArray
  --         numPoint
  --         deltaX
  --         numPoint
  --         deltaX
  --         phiFreqs
  --         rhoFreqs
  --         thetaFreqs
  --         scaleFreqs
  --         halfLogPeriod
  --         maxScale
  --         -- (fromIntegral numPoint)
  -- harmonicsArrayGPU <-
  --   dftHarmonicsArrayGPU
  --     plan
  --     numPoint
  --     deltaX
  --     numPoint
  --     deltaX
  --     phiFreqs
  --     rhoFreqs
  --     thetaFreqs
  --     scaleFreqs
  --     halfLogPeriod
  --     maxScale
  -- initialise []
  -- dev <- device . L.head $ gpuID
  -- ctx <- CUDA.create dev []
  -- ptx <- createTargetFromContext ctx
  -- registerPinnedAllocatorWith ptx
  --Source
  printCurrentTime "Source"
  sourceArr <- convolve' Source plan coefficients harmonicsArray initSource
  -- sourceArr' <- convolve'' Source plan coefficients harmonicsArray' initSource
  -- printCurrentTime "CPU"
  -- sourceArr <-
  --   convolveGPU
  --     ptx
  --     Source
  --     plan
  --     (A.use .
  --      A.fromList
  --        (A.Z A.:. (L.length scaleFreqs) A.:. (L.length thetaFreqs) A.:.
  --         (L.length rhoFreqs) A.:.
  --         (L.length phiFreqs)) .
  --      R.toList $
  --      coefficients)
  --     (A.use .
  --      A.fromList
  --        (A.Z A.:. (L.length scaleFreqs) A.:. (L.length thetaFreqs) A.:.
  --         (L.length rhoFreqs) A.:.
  --         (L.length phiFreqs)) .
  --      R.toList $
  --      coefficients)
  --     harmonicsArrayGPU
  --     initSource
  -- printCurrentTime "GPU"
  plotDFTArrayPower
    (folderPath </> (printf "SourcePower.png"))
    numPoint
    numPoint $
    sourceArr
  plotDFTArrayThetaR
    (folderPath </> "Source.png")
    numPoint
    numPoint
    thetaRHarmonics
    sourceArr
  let sourceArr' =
        fromUnboxed (Z :. numOrientation :. numPoint :. numPoint) .
        VG.convert .
        VG.concat .
        computeFourierSeriesThetaR thetaRHarmonics . getDFTArrayVector $
        sourceArr
  plotThetaDimension folderPath "FourierSeries_" location . R.map magnitude $
    sourceArr'
  -- plotDFTArrayThetaR
  --   (folderPath </> "Source_test.png")
  --   numPoint
  --   numPoint
  --   thetaRHarmonics
  --   sourceArr'
  --Sink
  printCurrentTime "Sink"
  sinkArr <- convolve' Sink plan coefficients harmonicsArray initSink
  plotDFTArrayPower (folderPath </> (printf "SinkPower.png")) numPoint numPoint $
    sinkArr
  plotDFTArrayThetaR
    (folderPath </> "Sink.png")
    numPoint
    numPoint
    thetaRHarmonics
    sinkArr
  printCurrentTime "Completion"
  completionArr <- completionField' plan sourceArr sinkArr
  plotDFTArrayThetaR
    (folderPath </> "Completion.png")
    numPoint
    numPoint
    thetaRHarmonics
    completionArr
  plotDFTArrayPower
    (folderPath </> (printf "CompletionPower.png"))
    numPoint
    numPoint
    completionArr
  let a =
        computeFourierSeriesThetaR thetaRHarmonics . getDFTArrayVector $
        sourceArr
      b =
        computeFourierSeriesThetaR thetaRHarmonics . getDFTArrayVector $ sinkArr
  plotImageRepa (folderPath </> "CompletionR2S1.png") .
    ImageRepa 8 .
    fromUnboxed (Z :. (1 :: Int) :. numPoint :. numPoint) .
    VG.map (\x -> x ^ 2) . VG.convert . L.foldl1' (VG.zipWith (+)) $
    L.zipWith (VG.zipWith (\x y -> magnitude x * magnitude y)) a b
  printCurrentTime ""
  -- M.mapM_
  --   (\s -> do
  --      let !arr =
  --            getNormalizedHistogramArr hist :: R.Array U DIM4 (Complex Double)
  --          initArr =
  --            traverse4
  --              (fromListUnboxed (Z :. L.length phiFreqs) phiFreqs)
  --              (fromListUnboxed (Z :. L.length rhoFreqs) rhoFreqs)
  --              (fromListUnboxed (Z :. L.length thetaFreqs) thetaFreqs)
  --              (fromListUnboxed (Z :. L.length scaleFreqs) scaleFreqs)
  --              (\_ _ _ _ -> extent arr) $ \fPhi fRho fTheta fScale (Z :. scale :. theta :. rho :. phi) ->
  --              cis $
  --              (-pi * log s) * fRho (Z :. rho) / halfLogPeriod -
  --              (fPhi (Z :. phi) * t0 * pi / 180)
  --      printCurrentTime ""
  --      source <- computeUnboxedP $ arr *^ initArr
  --      plotImageRepaComplex (folderPath </> (printf "Source_%.2f.png" s)) .
  --        ImageRepa 8 .
  --        fromUnboxed (Z :. (1 :: Int) :. numPoint :. numPoint) .
  --        L.foldl1' (VU.zipWith (+)) .
  --        computeFourierSeriesThetaR
  --          numOrientation
  --          numScale
  --          thetaFreqs
  --          scaleFreqs
  --          halfLogPeriod .
  --        computeFourierSeriesR2
  --          numPoint
  --          deltaX
  --          numPoint
  --          deltaX
  --          phiFreqs
  --          rhoFreqs
  --          thetaFreqs
  --          scaleFreqs
  --          halfLogPeriod
  --          maxScale
  --          harmonicsArray $
  --        source)
  --   [1,1.25 .. s0]

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 NamedFieldPuns #-}

module ParticleProcessingTests
  ( testParticleHandling
  ) where

import Control.Applicative (liftA2)
import Data.Complex (Complex(..))
import Data.List (sort)
import qualified Data.Vector as V
import HABQTlib.Data
import HABQTlib.Data.Particle
import HABQTlib.RandomStates
import qualified Numeric.LinearAlgebra as LA
import StateGenTests (arbPropDensityMatrix)
import qualified Streaming.Prelude as S
import qualified System.Random.MWC as MWC
import qualified Test.QuickCheck as QC
import Test.QuickCheck
  ( Args(..)
  , Positive(..)
  , Property
  , (.&&.)
  , arbitrary
  , quickCheckWith
  , stdArgs
  , suchThat
  )
import TestHelpers

lik :: Double -> Double
lik p = p / (1 - p)

loglik :: Double -> Double
loglik = log . lik

likProp :: PureStateVector -> DensityMatrix -> Property
likProp sv dm =
  let l = pureStateLikelihood sv dm
   in QC.property (l >= 0) .&&. QC.property (l <= 1)

testLikelihood :: IO ()
testLikelihood = do
  putStrLn "Testing likelihood computation for density matrices:"
  quickCheckWith stdArgs {maxSuccess = 1000} likProp

expProp :: QC.Blind Particles -> Property
expProp (QC.Blind pts) =
  let dm = snd . getWDM . V.foldl1' (<+>) . ptsParticles $ pts
      dim = LA.rows . getDensityMatrix $ dm
      z = (dim LA.>< dim) $ repeat 0
      scaleDM :: Weight -> DensityMatrix -> LA.Matrix (Complex Double)
      scaleDM w (DensityMatrix dm'') = LA.scale (w :+ 0) dm''
      dm' = DensityMatrix $ foldOverPts id scaleDM (+) z pts
      matchProp = QC.property $ dm <==> dm'
      validProp = arbPropDensityMatrix dm'
   in matchProp .&&. validProp

testExpectation :: IO ()
testExpectation = do
  putStrLn "Testing particle expectation function:"
  quickCheckWith stdArgs {maxSuccess = 500} expProp

particleUpdateFidR :: PureStateVector -> Particles -> Double
particleUpdateFidR testState@(PureStateVector sv) particles =
  let testPure = DensityMatrix $ sv LA.<> LA.tr sv
      originalEstimate = snd . getWDM $ reduceParticlesToMean particles
      originalFid = fidelityDM testPure originalEstimate
      updatedParticles = updateParticles testState particles
      updatedEstimate = snd . getWDM $ reduceParticlesToMean updatedParticles
      updatedFid = fidelityDM testPure updatedEstimate
   in loglik updatedFid - loglik originalFid

updateProp :: Property
updateProp =
  let gen = do
        (dim, rank) <- arbDimRank
        Positive num <- arbitrary :: QC.Gen (Positive NumberOfParticles)
        stateVector <- QC.resize dim arbitrary
        particles <- arbParticles dim rank (fromIntegral num)
        return (stateVector, particles)
      prop (sv, pts) =
        let wtsf Particles {ptsParticles} =
              V.foldl1' (+) . V.map (fst . getWDM) $ ptsParticles
            ws = wtsf pts
            pts' = updateParticles sv pts
            ws' = wtsf pts'
            wsProp = abs (ws - 1) < 1e-10
            ws'Prop = abs (ws' - 1) < 1e-10
         in wsProp .&&. ws'Prop
      msg = "of updates move the estimate closer"
   in QC.forAll
        gen
        (\(sv, pts) ->
           QC.classify (particleUpdateFidR sv pts >= 0) msg (prop (sv, pts)))

ecdfProp :: Particles -> Property
ecdfProp Particles {ptsParticles} =
  let cdf = ecdf ptsParticles
      hp = V.head cdf > 0
      tp = abs (V.last cdf - 1) < 1e-10
      cdfSorted = V.fromList . sort . V.toList $ cdf
      sorted = cdf == cdfSorted
   in hp .&&. tp .&&. sorted

icdfProp :: Particles -> Property
icdfProp Particles {ptsParticles} =
  QC.ioProperty $ do
    gen <- MWC.createSystemRandom
    x <- MWC.uniform gen
    let cdf = ecdf ptsParticles
        idx = icdf cdf x
        Just idx' = V.findIndex (> x) cdf
    return (idx == idx')

hDistProp :: ParticleHierarchy -> PureStateVector -> Double
hDistProp ph sv =
  let dm = DensityMatrix $ getStateVector sv LA.<> LA.tr (getStateVector sv)
      originalEstimate = getMixedEstimate ph
      updatedEstimate = getMixedEstimate $ updateParticleHierarchy sv ph
      originalFid = fidelityDM dm originalEstimate
      updatedFid = fidelityDM dm updatedEstimate
   in loglik updatedFid - loglik originalFid

hierarchyBatchProp ::
     Int -> Positive Dim -> Positive NumberOfParticles -> IO Property
hierarchyBatchProp batchSize (Positive dim) (Positive num) = do
  ph <- initialiseParticleHierarchy dim num
  passed <-
    S.length_ . S.filter (>= 0) . S.replicateM batchSize $
    hDistProp ph <$> genPureSV dim
  return . QC.property $
    fromIntegral passed / (fromIntegral batchSize :: Double) >= 0.95

hierarchyUpdateProp :: Int -> Property
hierarchyUpdateProp batchSize =
  let gen = do
        dim <- suchThat arbitrary (\x -> x > 1 && x < 20)
        num <- suchThat arbitrary (> 100)
        return (Positive dim, Positive num)
   in QC.forAll
        gen
        (\(dim, num) -> QC.ioProperty $ hierarchyBatchProp batchSize dim num)

effectiveSampleSizeProp :: Property
effectiveSampleSizeProp =
  let genP = do
        (dim, rank) <- arbDimRank
        Positive num <- arbitrary :: QC.Gen (Positive NumberOfParticles)
        arbParticles dim rank (fromIntegral num)
      sizeP ps@(Particles _ _ n _) =
        QC.classify
          (es <= fromIntegral n / 2)
          "of effective sizes are less than half of total particle number"
          ((es <= fromIntegral n) && (es >= 0))
        where
          es = effectiveSize ps
   in QC.forAll genP sizeP

nudgeProp :: Property
nudgeProp =
  let gen = do
        (dim, rank) <- arbDimRank
        dm <- arbDM dim rank
        return (dim, rank, dm)
      nudged (dim, rank, dm) =
        truncateRank rank <$> nudgeParticle dim 1e-2 (mkWDM1 dm)
      prop (dim, rank, dm) =
        QC.ioProperty $ do
          WeighedDensityMatrix (w, ndm) <- nudged (dim, rank, dm)
          let validDM = arbPropDensityMatrix ndm
              preservesRank = getRank dm == getRank ndm
              preservesWeight = w == 1
          return $
            QC.classify
              (fidelityDM ndm dm > 0.99)
              "fidelities exceed 99%"
              (validDM .&&. preservesRank .&&. preservesWeight)
   in QC.forAll gen prop

type ResamplingFunIO = MWC.GenIO -> Particles -> IO Particles

resampleProp :: MWC.GenIO -> ResamplingFunIO -> Property
resampleProp gen resamplingFun =
  let genA = do
        dim <- suchThat (arbitrary :: QC.Gen Dim) (liftA2 (&&) (< 10) (> 1))
        Positive rank <-
          suchThat (arbitrary :: QC.Gen (Positive Rank)) (< Positive dim)
        num <- suchThat arbitrary (> 1000)
        QC.Blind <$> arbParticles dim rank num
      preservesMeanIO :: QC.Blind Particles -> IO Property
      preservesMeanIO (QC.Blind pts) = do
        rpts <- resamplingFun gen pts
        let WeighedDensityMatrix (w1, dm1) = reduceParticlesToMean pts
            WeighedDensityMatrix (w2, dm2) = reduceParticlesToMean rpts
            sumw wacc (WeighedDensityMatrix (w', _)) = wacc + w'
            getpw Particles {ptsParticles = ps} = V.foldl' sumw 0 ps
            w1' = getpw pts
            w2' = getpw rpts
            propWeight = QC.property $ 2 * abs (w1 - w2) / (w1 + w2) < 1e-12
            propWeight' =
              QC.property $ 2 * abs (w1' - w2') / (w1' + w2') < 1e-12
            Particles _ _ num _ = rpts
            propESabs =
              QC.property $ abs (effectiveSize rpts - fromIntegral num) < 1
            propESincrease =
              QC.property $ effectiveSize rpts > effectiveSize pts
            propAll =
              propWeight .&&. propWeight' .&&. propESincrease .&&. propESabs
        return $
          QC.classify
            (fidelityDM dm1 dm2 > 0.99)
            "fidelities between pre- and post-resampling means exceed 99%"
            propAll
   in QC.forAll genA (QC.ioProperty . preservesMeanIO)

testParticleResampling :: IO ()
testParticleResampling = do
  gen <- MWC.createSystemRandom
  putStrLn "Testing multinomial particle resampling without checks:"
  quickCheckWith stdArgs {maxSuccess = 100} $ resampleProp gen resampleMultinom

testParticleNudging :: IO ()
testParticleNudging = do
  putStrLn "Testing particle nudging:"
  quickCheckWith stdArgs {maxSuccess = 1000} nudgeProp

testEffectiveSampleSize :: IO ()
testEffectiveSampleSize = do
  putStrLn "Testing calculation of effective sample size:"
  quickCheckWith stdArgs {maxSuccess = 1000} effectiveSampleSizeProp

testParticleUpdate :: IO ()
testParticleUpdate = do
  putStrLn "Testing update of particles:"
  quickCheckWith stdArgs {maxSuccess = 1000} updateProp

testEcdf :: IO ()
testEcdf = do
  putStrLn "Testing ecdf calculation:"
  quickCheckWith stdArgs {maxSuccess = 10000} ecdfProp
  putStrLn "Testing iecdf calculation:"
  quickCheckWith stdArgs {maxSuccess = 10000} icdfProp

testHierarchyUpdate :: IO ()
testHierarchyUpdate = do
  putStrLn "Testing update of particle hierarchies:"
  quickCheckWith stdArgs {maxSuccess = 50} (hierarchyUpdateProp 100)

testParticleHandling :: IO ()
testParticleHandling = do
  testLikelihood
  testExpectation
  testParticleUpdate
  testEcdf
  testHierarchyUpdate
  testEffectiveSampleSize
  testParticleNudging
  testParticleResampling

ext: hs

max_stars_repo_name: Belinsky-L-V/HABQT

max_issues_repo_name: Belinsky-L-V/HABQT

max_forks_repo_name: Belinsky-L-V/HABQT


text:
{-# LANGUAGE ExistentialQuantification #-}
{-# LANGUAGE FlexibleInstances         #-}
{-# LANGUAGE MultiParamTypeClasses     #-}

module ComposableSDR.Types
  ( compose
  , addPipe
  , unPipe
  , toArray
  , elemSize
  , sampleFormat
  , Demodulator
  , Pipe(..)
  , ArrayPipe
  , AudioFormat(..)
  , SamplesIQCF32
  , SoapyException(..)
  , Array(..)
  ) where

import           Control.Category                           (Category (..))
import           Prelude                                    hiding ((.))

import           Data.Complex
import           Data.Typeable

import           Foreign.C.Types
import           Foreign.ForeignPtr                         (plusForeignPtr,
                                                             withForeignPtr)
import           Foreign.Ptr

import           Control.Exception                          (Exception)
import           Control.Monad
import qualified Control.Monad.Catch                        as MC
import           Control.Monad.State
import qualified Control.Monad.Trans.Control                as MTC

import qualified Sound.File.Sndfile                         as SF

import qualified Streamly.Internal.Data.Fold.Types          as FL
import           Streamly.Internal.Data.Stream.StreamK.Type (IsStream)
import qualified Streamly.Internal.Memory.Array.Types       as AT (Array (..))
import qualified Streamly.Memory.Array                      as A
import qualified Streamly.Prelude                           as S

data AudioFormat
  = AU
  | WAV
  deriving (Show, Read)

data Pipe m a b = forall r. Pipe
  { _start :: m r
  , _process :: r -> a -> m b
  , _done :: r -> m ()
  }

newtype Array a = Array
  { fromArray :: AT.Array a
  }

instance SF.Buffer Array Float where
  fromForeignPtr fp i n =
    withForeignPtr fp $ \p -> do
      let v =
            AT.Array
              { AT.aStart = fp `plusForeignPtr` i
              , AT.aEnd = p `plusPtr` (i + 4 * n)
              , AT.aBound = p `plusPtr` (i + 4 * n)
              }
      return $ toArray v
  toForeignPtr = return . (\x -> (AT.aStart x, 0, A.length x)) . fromArray

type ArrayPipe m a b = Pipe m (A.Array a) (A.Array b)
type Demodulator = ArrayPipe IO SamplesIQCF32 Float

data SoapyException =
  SoapyException
  deriving (Show, Typeable)

instance Exception SoapyException

type SamplesIQCF32 = Complex CFloat

elemSize :: Int
elemSize = 8

sampleFormat :: String
sampleFormat = "CF32"

toArray :: AT.Array a -> Array a
toArray = Array

compose :: Monad m => Pipe m b c -> Pipe m a b -> Pipe m a c
compose (Pipe start1 process1 done1) (Pipe start2 process2 done2) =
  Pipe start process done
  where
    start = (,) <$> start1 <*> start2
    done (r1, r2) = done2 r2 >> done1 r1
    process (r1, r2) = process2 r2 >=> process1 r1

instance Monad m => Category (Pipe m) where
  id = Pipe (return ()) (const return) (const $ return ())
  (.) = compose

instance Monad m => Functor (Pipe m a) where
  fmap f (Pipe start process done) =
    Pipe start (\r a -> fmap f (process r a)) done

unPipe :: (IsStream t, MonadIO m,
          MTC.MonadBaseControl IO m,
          MC.MonadThrow m) =>
          Pipe m a b -> m (t m a -> t m b, m ())
unPipe (Pipe creat process dest) = do
  r <- creat
  return (S.mapM (process r), dest r)

addPipe :: MonadIO m => Pipe m a b -> FL.Fold m b c -> FL.Fold m a c
addPipe (Pipe creat process dest) (FL.Fold step1 start1 done1) =
  FL.Fold step start done
  where
    start = do
      s <- start1
      r <- creat
      return (s, r, process r)
    done (s, r, _) = do
      dest r
      done1 s
    step (s, r, f) a = do
      b <- f a
      s' <- step1 s b
      return (s', r, f)

max_stars_repo_name: mryndzionek/composable-sdr

max_issues_repo_name: mryndzionek/composable-sdr

max_forks_repo_name: mryndzionek/composable-sdr


text:
{-# LANGUAGE BangPatterns #-}
{- |
   Description :  Functions to work with transition probability matrices on rooted trees
   Copyright   :  (c) Dominik Schrempf 2017
   License     :  GPLv3

   Maintainer  :  dominik.schrempf@gmail.com
   Stability   :  unstable
   Portability :  non-portable (not tested)

Calculate transition probability matrices, map rate matrices on trees, populate
a tree with states according to a stationary distribution, etc.

The implementation of the Markov process is more than basic and can be improved in a lot of ways.

TODO: Try https://hackage.haskell.org/package/Dist-0.4.1.0/docs/Numeric-Probability-Distribution.html

* Changelog

-}

module Transition where
import           Control.Monad.Random.Strict hiding (fromList)
import           Data.Distribution           hiding (toList)
import           Numeric.LinearAlgebra       hiding (fromList)
import           RateMatrix
import           RTree

-- | This may be moved to a different module ProbMatrix or alike.
type ProbMatrix   = Matrix R

-- | The important matrix that gives the probabilities to move from one state to
-- another in a specific time (branch length).
probMatrix :: RateMatrix -> BranchLn -> ProbMatrix
probMatrix m t = exponential
  where !exponential = expm $ scale t m

-- | Convert a tree with branch lengths into a tree that has transition
-- probability matrices assigned to each of its branches.
branchLengthsToTransitionProbs :: RateMatrix -> RTree a Double -> RTree a ProbMatrix
branchLengthsToTransitionProbs m = fmap (probMatrix m)

-- | Move from a given state to a new one according to a transition probability
-- matrix (for performance reasons this probability matrix needs to be given as
-- a list of generators, see
-- https://hackage.haskell.org/package/distribution-1.1.0.0/docs/Data-Distribution-Sample.html).
-- This function is the bottleneck of the simulator and takes up most of the
-- computation time. However, I was not able to find a faster implementation
-- than the one from Data.Distribution.
jump :: (MonadRandom m) => State -> [Generator State] -> m State
jump (State s) p = target
  where !target = getSample $ p !! s

-- | Perform N jumps from a given state and according to a transition
-- probability matrix transformed to a list of generators. This implementation
-- uses 'foldM' and I am not sure how to access or store the actual chain. This
-- could be done by an equivalent of 'scanl' for general monads, which I was
-- unable to find. This function is neat, but will most likely not be needed.
-- However, it is instructive and is left in place.
jumpN :: (MonadRandom m) => State -> [Generator State] -> Int -> m State
jumpN s p n = foldM jump s (replicate n p)

-- | This is the heart of the simulation. Take a tree (most probably with node
-- labels of type a, but this could be anything) and a root state. If there is a
-- split (i.e., if we have a node and not a leaf), jump down the left branch and
-- populate and flatten the tree and jump down the right branch and populate and
-- flatten the tree and append the two results. This is what 'liftM2' is doing,
-- it lifts the append function of lists (++) to the 'Rand g' monad. If we
-- encounter a leaf, just return the node label (or whatever the type a is) and
-- the state that we ended up at.
populateAndFlattenTree :: (MonadRandom m) => RTree a [Generator State] -> State -> m [(a, State)]
populateAndFlattenTree (Leaf a) s = return [(a, s)]
populateAndFlattenTree (Node _ lp lc rp rc) s = liftM2 (++) (jumpDownBranch lp lc) (jumpDownBranch rp rc)
  where jumpDownBranch p t = jump s p >>= populateAndFlattenTree t

-- | Convert a stationary distribution to a generator which speeds up random
-- picks. Cf. 'jump'.
stationaryDistToGenerator :: StationaryDist -> Generator State
stationaryDistToGenerator f = fG
  -- This is a little complicated. I need to convert the vector to a list to be
  -- able to create a distribution.
  where !fL = toList f
        !fD = fromList $ zip (map State [0..]) fL
        !fG = fromDistribution fD

-- | See 'jump'.
treeProbMatrixToTreeGenerator :: RTree a ProbMatrix -> RTree a [Generator State]
treeProbMatrixToTreeGenerator t = tG
  where
    -- Create a tree with the probability matrices as list of row vectors.
    !tL = fmap toLists t
    -- A complicated double map. We need to create generators for each branch on
    -- the tree (fmap) and for each target state on each branch (map).
    !tG = (fmap . map) (fromDistribution . fromList . zip (map State [0..])) tL

-- | Simulate data (states at the leaves) for a tree with transition
-- probabilities on its branches and with the stationary distribution of states
-- at the root. The state at the root is randomly chosen from the stationary
-- distribution and the states at the nodes and leaves are randomly chosen
-- according to the transition probabilities.
simulateSite :: (MonadRandom m) =>
                Generator State
             -> RTree a [Generator State]
             -> m [(a, State)]
simulateSite f t = do
  !rootState <- getSample f
  populateAndFlattenTree t rootState

-- | Randomly draw an index according to a given generator. Use the stationary
-- distribution and rooted tree at the drawn index to simulate a site. This is
-- useful for simulation, e.g., Gamma rate heterogeneity models.
simulateSiteGen :: (MonadRandom m) =>
                   Generator Int
                -> [Generator State]
                -> [RTree a [Generator State]]
                -> m [(a, State)]
simulateSiteGen gen fs trs = do
  !i <- getSample gen
  let !f = fs  !! i
      !t = trs !! i
  simulateSite f t

max_stars_repo_name: dschrempf/bmm-simulate

max_issues_repo_name: pomo-dev/bmm-simulate

max_forks_repo_name: pomo-dev/bmm-simulate


text:
{-# LANGUAGE OverloadedStrings #-}

import Control.Applicative ((<$>))
import Statistics.Sample.KernelDensity (kde)
import Text.Hastache (MuType(..), defaultConfig, hastacheFile)
import Text.Hastache.Context (mkStrContext)
import qualified Data.Attoparsec as B
import qualified Data.Attoparsec.Char8 as A
import qualified Data.ByteString as B
import qualified Data.ByteString.Lazy as L
import qualified Data.Vector.Unboxed as U

csv = do
  B.takeTill A.isEndOfLine
  (A.double `A.sepBy` A.char ',') `A.sepBy` A.endOfLine

main = do
  waits <- (either error (U.fromList . map last . filter (not.null)) .
            A.parseOnly csv) <$> B.readFile "data/faithful.csv"
  let xs = map (\(a,b) -> [a,b]) . U.toList . uncurry U.zip . kde 64 $ waits
      context "data" = MuVariable . show $ xs
  s <- hastacheFile defaultConfig "kde.tpl" (mkStrContext context)
  L.writeFile "kde.html" s

max_stars_repo_name: StefanHubner/statistics

max_issues_repo_name: StefanHubner/statistics

max_forks_repo_name: StefanHubner/statistics


text:
{- |
Description :  Default values and constants
Copyright   :  (c) Dominik Schrempf 2017
License     :  GPLv3

Maintainer  :  dominik.schrempf@gmail.com
Stability   :  unstable
Portability :  non-portable (not tested)

Default values and constants.

* Changelog

-}

module Defaults where

import           DNAModel              ( DNAModelSpec(..))
import           Numeric.LinearAlgebra
import qualified System.Random         as Rand

-- | Output file name.
outFileName :: String
outFileName = "Test"

-- | Virtual population size.
popSize :: Int
popSize = 9

-- | A kappa value of 6.25 corresponds to a transition to transversion ratio of
-- 3.0
dnaModelSpec :: DNAModelSpec
dnaModelSpec = HKY 6.25 (vector [0.3, 0.2, 0.2, 0.3])

-- | Heterozygosity value.
heterozygosity :: Double
heterozygosity = 0.0025

-- | Tree height.
treeHeight :: Double
treeHeight = 0.005

-- | Tree type. At the moment, ILS (incomplete lineage sorting) and Yule trees
-- are supported.
treeType :: String
treeType = "ILS"

-- | Number of sites to simulate.
nSites :: Int
nSites = 1000000

-- | By default, the seed is random.
seed :: IO Rand.StdGen
seed = Rand.getStdGen


text:
{-# OPTIONS_GHC -Wall #-}
{-# LANGUAGE CPP #-}

{-|

  Module      : AutoBench.Internal.Regression
  Description : Generating and fitting models for regression analysis.
  Copyright   : (c) 2018 Martin Handley
  License     : BSD-style
  Maintainer  : martin.handley@nottingham.ac.uk
  Stability   : Experimental
  Portability : GHC

  This module is responsible for generating and fitting regression models.
  See 'LinearType', 'LinearCandidate' and 'LinearFit' for more details.

  The system uses ridge regression because it has a closed-form solution so
  is fast. Previously used LASSO implemented using the HVX library but
  it was too computationally expensive, (see bottom of this file).

-}

{-
   ----------------------------------------------------------------------------
   <TO-DO>:
   ----------------------------------------------------------------------------
   - Too much white space;
   -
-}

module AutoBench.Internal.Regression
  (

  -- * Candidate generation
    generateLinearCandidate    -- Generate a 'LinearCandidate' from a 'LinearType' by filling
                               -- in the details of the model necessary to fit it to a given dataset.
  -- * Regression analysis
  , fitOLS                     -- Given a data set (a list of xy-'Coord's) and a 'LinearCandidate',
                               -- fit the model to the data using the OLS method.
  , fitRidgeRegress            -- Given a data set (a list of xy-'Coord's) and a 'LinearCandidate',
                               -- fit the model to the data using ridge regression.

  ) where

import           Data.List            (minimumBy, nub)
import           Data.Ord             (comparing)
import qualified Data.Vector.Storable as V

import Numeric.LinearAlgebra
  ( Element
  , Matrix
  , Vector
  , (<>)
  , (<\>)
  , (#>)
  , asColumn
  , cmap
  , diag
  , dot
  , fromColumns
  , fromRows
  , konst
  , norm_2
  , norm_Inf
  , qr
  , rows
  , scale
  , subVector
  , sumElements
  , takeRows
  , toColumns
  , tr'
  , vjoin
  )

import qualified AutoBench.Internal.Expr as E
import AutoBench.Internal.Types
  ( Coord
  , Exp
  , LinearCandidate(..)
  , LinearType(..)
  , numPredictors
  )

#ifdef MIN_VERSION_base
#if MIN_VERSION_base(4,11,0)
import Prelude hiding ((<>))
#endif
#endif

-- * 'LinearCandidate' generation

-- | Generate a 'LinearCandidate' from a 'LinearType' by filling in the
-- details of the model necessary to fit it to a given dataset.
-- Makes use of the generation helper functions 'generateLinearCandidateExpr'
-- and 'generateLinearCandidateYHat'.
generateLinearCandidate :: LinearType -> LinearCandidate
generateLinearCandidate t@(Poly _) =
  LinearCandidate
   {
     _lct   = t
   , _fxs   = id
   , _fex   = generateLinearCandidateExpr t
   , _fyhat = generateLinearCandidateYHat t
   }
generateLinearCandidate t@(Log b _) =
  LinearCandidate
   {
     _lct   = t
   , _fxs   = cmap $ logBase (fromIntegral b)
   , _fex   = generateLinearCandidateExpr t
   , _fyhat = generateLinearCandidateYHat t
   }
generateLinearCandidate t@(PolyLog b _) =
  LinearCandidate
   {
     _lct   = t
   , _fxs   = cmap (\x -> x * logBase (fromIntegral b) x)
   , _fex   = generateLinearCandidateExpr t
   , _fyhat = generateLinearCandidateYHat t
   }
generateLinearCandidate t@(Exp n) =
  LinearCandidate
   {
     _lct   = t
   , _fxs   = cmap (flip (**) $ fromIntegral n)
   , _fex   = generateLinearCandidateExpr t
   , _fyhat = generateLinearCandidateYHat t
   }

-- * Linear fitting

-- | Given a data set (a list of xy-'Coord's) and a 'LinearCandidate',
-- fit the model to the data using the ordinary least squares (OLS) method. The
-- result is the model's coefficients.
fitOLS :: [Coord] -> LinearCandidate -> Vector Double
fitOLS coords c = head $ toColumns $ (_preds <\> _resps)
  where (_preds, _resps) = genFitMatrices coords c

-- | Given a data set (a list of xy-'Coord's) and a 'LinearCandidate',
-- fit the model to the data using ridge regression. The result is the model's
-- coefficients.
fitRidgeRegress :: [Coord] -> LinearCandidate -> Vector Double
fitRidgeRegress coords c =
  -- Select the lambda that minimises mean squared error.
  fst $ minimumBy (comparing snd) $ zip fits errs
  where
    n        = length coords
    (xs, ys) = unzip  coords
    vXs      = V.fromList xs
    vYs      = V.fromList ys
    -- Fitting error is mean squared error.
    calcMse m vec = ((norm_2 $ V.zipWith (-) vYs vec) ** 2.0) / fromIntegral m
    -- Generate the matrices to solve.
    (_preds, _resps) = genFitMatrices coords c
    ---------------------------------------------------------------------------
    ---------------------------------------------------------------------------
    -- Using ridge regression with maximum tuning parameter:
    -- https://stats.stackexchange.com/questions/289075/what-is-the-
    -- smallest-lambda-that-gives-a-0-component-in-lasso?noredirect=1&lq=1
    ---------------------------------------------------------------------------
    maxLam   = (1.0 / fromIntegral n) * norm_Inf (tr' _preds <> _resps)
    -- Minimise 'mse' to select best tuning parameter.
    lamStep  = maxLam / 10
    lamRange = [lamStep, 2 * lamStep .. maxLam]
    fits     = nub $ fmap (head
                            . toColumns
                            . ridgeRegress _preds _resps
                          ) lamRange
    -- ** Do not use vFxs here, use vXs **
    errs     = fmap (calcMse n . flip (_fyhat c) vXs) fits
    ---------------------------------------------------------------------------
    ---------------------------------------------------------------------------
    -- Closed form solution of ridge regression.
    -- https://github.com/chris-taylor/aima-haskell/blob/master/src/AI/
    -- Learning/LinearRegression.hs
    ---------------------------------------------------------------------------
    ridgeRegress
      :: Matrix Double  -- X
      -> Matrix Double  -- Y
      -> Double         -- λ
      -> Matrix Double  -- (X'X + λI)^−1 X'Y.
    ridgeRegress x y lambda =
      let m      = rows x
          (_, r) = qr x
          rr     = takeRows m r
          ww     = diag $ vjoin [0, konst 1 (m - 1)]
      in (tr' rr <> rr + (lambda `scale` ww)) <\> (tr' x <> y)

-- * Helpers

-- | Generate an expression ('Exp') for a 'LinearCandidate' that will allow
-- its model's equation to be pretty printed at a later time. Important note:
-- the vectors are assumed to be of the correct form, /no/  checking is
-- performed currently.
generateLinearCandidateExpr
   :: LinearType      -- Model.
   -> Vector Double   -- Coefficients of the model.
   -> Exp             -- Model's equation as an expression for pretty printing.
--
-- Given coefficients:
--
--   a_0, a_1, ..., a_n
--
-- Generate:
--
--   a_0 + a_1 * x^1 + a_2 * x^2 + ... + a_n * x^n
--
generateLinearCandidateExpr (Poly n) coeffs =
  foldr1 E.Add . take (n + 1) $ zipWith ($) terms (V.toList coeffs)
  where terms = E.Num                                         -- a_0
                 : ( \a_k -> E.Mul (E.Num a_k) (E.Var "x")    -- a_1 * x
                   )
                 : fmap ( \p a_k -> E.Mul (E.Num a_k)         -- a_k * x^k, k > 1
                                          (E.Pow (E.Var "x")
                                                 (E.Num p)
                                          )
                        ) [2.0..]
--
-- Given coefficients:
--
-- a_0, a_1, ..., a_n
--
-- Generate:
--
-- a_0 + a_1 * log_b^1(x) + a_2 * log_b^2(x) + ... + a_n * log_b^n(x)
--
generateLinearCandidateExpr (Log b n) coeffs =
  foldr1 E.Add $ take (n + 1) $ zipWith ($) terms (V.toList coeffs)
  where terms = E.Num                                                -- a_0
                 : ( \a_k -> E.Mul (E.Num a_k)                       -- a_1 * log_b(x)
                                   (E.Log (E.Num $ fromIntegral b)
                                          (E.Var "x")
                                   )
                   )
                 : fmap ( \p a_k -> E.Mul (E.Num a_k)                -- a_k * log_b^k(x), k > 1
                                          (E.Pow (E.Log (E.Num $ fromIntegral b)
                                                        (E.Var "x")
                                                 )
                                                 (E.Num p)
                                          )
                        ) [2.0..]
--
-- Given coefficients:
--
-- a_0, a_1, ..., a_n
--
-- Generate:
--
-- a_0 + a_1 * x^1 * log_b^1(x) + a_2 * x^2 * log_b^2(x) + ... + a_n * x^n * log_b^n(x)
--
generateLinearCandidateExpr (PolyLog b n) coeffs =
  foldr1 E.Add $ take (n + 1) $ zipWith ($) terms (V.toList coeffs)
  where terms = E.Num                                                -- a_0
                 : ( \a_k -> E.Mul (E.Num a_k)                       -- a_1 * x * log_b(x)
                                   (E.Mul (E.Var "x")
                                          (E.Log (E.Num $ fromIntegral b)
                                                 (E.Var "x")
                                          )
                                   )
                   )
                 : fmap ( \p a_k -> E.Mul (E.Num a_k)                -- a_k * x^k * log_b^k(x), k > 1
                                          (E.Mul (E.Pow (E.Var "x")
                                                        (E.Num p)
                                                 )
                                                 (E.Pow (E.Log (E.Num $ fromIntegral b)
                                                               (E.Var "x")
                                                        )
                                                        (E.Num p)
                                                 )
                                          )
                        ) [2.0..]
--
-- Given coefficients:
--
-- a_0, a_1
--
-- Generate:
--
-- a_0 + a_1 * n^x
--
generateLinearCandidateExpr (Exp n) coeffs =
  foldr1 E.Add $ zipWith ($) terms (V.toList coeffs)
  where terms = [ E.Num  -- a_0
                , \a_1 -> E.Mul (E.Num a_1)   -- a_1 *
                                (E.Pow (E.Num $ fromIntegral n) -- n^x
                                        (E.Var "x")
                                )
                ]

-- | For a given set of x-coordinates, generate the y-coordinates that are
-- predicted by the model of a 'LinearCandidate'. Important note: the vectors
-- are assumed to be of the correct form, /no/ checking is performed currently.
generateLinearCandidateYHat
  :: LinearType
  -> Vector Double   -- Coefficients of the model.
  -> Vector Double   -- x-coordinates.
  -> Vector Double   -- Predicted y-coordinates.
--
-- Given coefficients and x-coords:
--
--   a_0, a_1, ..., a_n     x_0, x_1, ..., x_m
--
-- For each x_i, generate:
--
--   a_0 + a_1 * x_i^1 + a_2 * x_i^2 + ... + a_n * x_i^n
--
generateLinearCandidateYHat (Poly n) coeffs xs =
  (fromColumns $ V.fromList (replicate m 1.0) : xs') #> coeffs
  where
    m   = V.length xs
    xs' = zipWith vecPow [1.0..] (replicate n xs)
--
-- Given coefficients and x-coords:
--
-- a_0, a_1, ..., a_n     x_0, x_1, ..., x_m
--
-- For each x_i, generate:
--
-- a_0 + a_1 * log_b^1(x_i) + a_2 * log_b^2(x_i) + ... + a_n * log_b^n(x_i)
--
generateLinearCandidateYHat (Log b n) coeffs xs =
  (fromColumns $ V.fromList (replicate m 1.0) : xs') #> coeffs
  where
    m   = V.length xs
    xs' = zipWith (vecLogPow $ fromIntegral b) [1.0..] (replicate n xs)
--
-- Given coefficients and x-coords:
--
-- a_0, a_1, .., a_n     x_0, x_1, ..., x_m
--
-- For each x_i, generate:
--
-- a_0 + a_1 * x_i^1 * log_b^1(x_i) + a_2 * x_i^2 * log_b^2(x_i) + ... + a_n * x_i^n * log_b^n(x_i)
--
generateLinearCandidateYHat (PolyLog b n) coeffs xs =
  (fromColumns $ V.fromList (replicate m 1.0) : xs') #> coeffs
  where
    m   = V.length xs
    xs' = zipWith (vecPolyLogPow $ fromIntegral b) [1.0..] (replicate n xs)
--
-- Given coefficients and x-coords:
--
-- a_0, a_1      x_0, x_1, ..., x_m
--
-- For each x_i, generate:
--
-- a_0 + a_1 * n^x_i
--
generateLinearCandidateYHat (Exp n) coeffs xs =
  (fromColumns $ V.fromList (replicate m 1.0) : xs') #> coeffs
  where
    m   = V.length xs
    xs' = zipWith (vecBasePow $ fromIntegral n) [1.0..] (replicate 1 xs)

-- | Constructs predictor (@preds@) and responder (@resps@) matrices so that
-- regression equations can be solved in their matrix form.
genFitMatrices :: [Coord] -> LinearCandidate -> (Matrix Double, Matrix Double)
genFitMatrices coords c = gen (numPredictors (_lct c) - 1)
  where

    gen :: Int -> (Matrix Double, Matrix Double)
    gen k = let
               _xPows    = xPows k
               _xSumPows = xSumPows _xPows
               _ySums    = ySums _xPows k
               _preds    = preds _xSumPows k
               _resps    = resps _ySums
            in (_preds, asColumn _resps)

    ---------------------------------------------------------------------------
    --
    --
    (xs, ys) = unzip coords
    -- Transform the x-coords using the candidate's 'fxs' function. This allows
    -- us to fit different models to the data in a generic way. For example, a
    -- 'Log b n' candidate will first take 'LogBase b' of the x-coords. We can
    -- then use (a variant of) linear least squares to fit a straight line
    -- through the resulting data points. Although this is a straight line on
    -- the /resulting data/, it corresponds to a logarithmic function on the
    -- /initial/ data. This technique is standard.
    vFxs = (_fxs c) (V.fromList xs)
    vYs  = V.fromList ys
    -- Number of data points.
    n    = length coords
    -- As the x-'Coords' have been transformed by the candidate's 'fxs'
    -- function, this method works generically for any model and 'k'
    -- coefficients.
    ---------------------------------------------------------------------------
    --
    -- Helper functions to construct the matrices to solve using linear
    -- regression. Note: these functions are to be used in a particular way so
    -- their names reflect their use.
    --
    -- x_i^j for j = 1, 2, .., 2 * k
    xPows :: Int -> [Vector Double]
    xPows k = zipWith vecPow [1.0..] (replicate (2 * k) vFxs)
    -- SUM x_i^j for j = 1, 2, .., 2 * k
    xSumPows :: [Vector Double] -> [Double]
    xSumPows  = fmap sumElements
    -- SUM y_i, SUM x_i * y_i, SUM x_i^2 * y_i, .., SUM x_i^j * y_i
    ySums :: [Vector Double] -> Int -> [Double]
    ySums vs k = zipWith dot (V.fromList (replicate n 1.0) : vs) (replicate (k + 1) vYs)
    --
    -- Convert xSumPows into a matrix by 'chunking' off rows from the initial vector.
    --
    -- Rows are length k + 1, the initial rows are:
    --
    --   n      , SUM x_i  , SUM x_i^2, .., SUM x_i^k
    --
    --   SUM x_i, SUM x_i^2, SUM x_i^3, .., SUM x_i^(k + 1)
    --
    -- The last row is:
    --
    --   SUM x_i^j, SUM x_i^(j + 1), .., SUM x_i^(2 * k)
    --
    -- These values are the /predictors/.
    --
    preds :: [Double] -> Int -> Matrix Double
    preds ds k = vecToMatrix (V.fromList $ fromIntegral n : ds) (k + 1)
    -- 'ySums' are the /responders/.
    resps :: [Double] -> Vector Double
    resps  = V.fromList

-- | Raise each double in a vector of doubles to a given power.
vecPow :: Double -> Vector Double -> Vector Double
vecPow pow = cmap (** pow)

-- | For each double 'x' in a vector of doubles, apply the following
-- calculation:
--
-- > \x -> base ** (x * pow)
--
-- for a given 'base' and 'pow'.
vecBasePow :: Double -> Double -> Vector Double -> Vector Double
vecBasePow base pow = cmap (\x -> base ** (x * pow))

-- | For each double 'x' in a vector of doubles, apply the following
-- calculation:
--
-- > \x -> (x ** pow) * (logBase b x ** pow)
--
-- for a given 'b' and 'pow'.
vecPolyLogPow :: Double -> Double -> Vector Double -> Vector Double
vecPolyLogPow b pow = cmap (\x -> (x ** pow) * ((logBase b x) ** pow))

-- | For each double 'x' in a vector of doubles, apply the following
-- calculation:
--
-- > \x -> logBase b x ** pow
--
-- for a given 'b' and 'pow'.
vecLogPow :: Double -> Double -> Vector Double -> Vector Double
vecLogPow b pow = cmap (\x -> (logBase b x) ** pow)

-- | Create a matrix by 'chunking' off rows of a given length from a vector.
-- Works similarly to 'tails'.
vecToMatrix :: (V.Storable a, Element a) => Vector a -> Int -> Matrix a
vecToMatrix v j = fromRows [ subVector i j v | i <- [0..j - 1] ]

{-

LASSO has been retired for the time being due to computational expense.

{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE DataKinds        #-}

import           HVX.Internal.Constraints ((<=~))
import           HVX.Primitives           ((+~), (*~), neg, norm, quadform)
import           HVX.Internal.Solvers     (ellipsoidMinimize)
import qualified HVX.Internal.Primitives  as HVX

-- ** LASSO (Least Absolute Shrinkage and Selection Operator)

fitLASSO :: [Coord] -> LinearCandidate -> Vector Double
fitLASSO coords c =
  -- Select the lambda that minimises mean squared error.
  fst $ minimumBy (comparing snd) $ zip fits errs
  where
    n        = length coords
    (xs, ys) = unzip  coords
    vXs      = V.fromList xs
    vYs      = V.fromList ys
    -- Fitting error is mean squared error.
    calcMse m vec = ((norm_2 $ V.zipWith (-) vYs vec) ** 2.0) / fromIntegral m
    -- Generate the matrices to solve.
    (_preds, _resps) = genFitMatrices coords c
    ---------------------------------------------------------------------------
    ---------------------------------------------------------------------------
    initCoeffs = _preds <\> _resps
    maxT = V.sum $ cmap abs $ flatten initCoeffs
    -- Minimise 'mse' to select best tuning parameter.
    tStep  = maxT / 10
    tRange = [0 + tStep, 0 + (2 * tStep) .. maxT]
    fits   = filter (all (\m -> not (isNaN m || isInfinite m)) . V.toList)
              . nub
              $ fmap (head
                       . toColumns
                       . lassoRegress _preds _resps
                     ) tRange
    -- ** Do not use vFxs here, use vXs **
    errs   = fmap (calcMse n . flip (fyhat c) vXs) fits
    ---------------------------------------------------------------------------
    ---------------------------------------------------------------------------
    -- LASSO regression (convex equation using the subgradient method in HVX).
    lassoRegress
      :: Matrix Double      -- ^ n >< n matrix X.
      -> Matrix Double      -- ^ n >< 1 matrix Y.
      -> Double             -- ^ Constraint on sum of coefficients.
      -> Matrix Double      -- ^ Optimal variable values.
    lassoRegress x y t = snd . head $ res
      where
        (res, _, _) = ellipsoidMinimize expr [constr] [("a", m)] 1e-20 (2 * maxT)
        m      = rows x
        constr = norm 1 (HVX.EVar "a") <=~ (HVX.EConst $ konst t (1, 1))
        expr   = quadform (HVX.EConst $ ident m) $ ((HVX.EConst y) +~ neg (HVX.EConst x *~ HVX.EVar "a"))

-}

max_stars_repo_name: recursion-ninja/AutoBench

max_issues_repo_name: recursion-ninja/AutoBench

max_forks_repo_name: recursion-ninja/AutoBench


text:
{-# language GADTs                 #-}
{-# language StandaloneDeriving    #-}
{-# language TypeOperators         #-}
{-# language FlexibleInstances     #-}
{-# language FlexibleContexts      #-}
{-# language UndecidableInstances  #-}
{-# language MultiParamTypeClasses #-}
{-# language ConstraintKinds       #-}
{-# language TypeFamilies          #-}

{-# options_ghc -fwarn-incomplete-patterns #-}

module Feldspar.Software.Expression where

import Feldspar.Sugar
import Feldspar.Representation
import Feldspar.Software.Primitive
import Data.Struct

import Data.Complex (Complex)
import Data.Int
import Data.Word
import Data.List (genericTake)
import Data.Constraint
import Data.Typeable (Typeable)

-- syntactic.
import Language.Syntactic hiding (Signature, Args)
import Language.Syntactic.Functional hiding (Lam)
import Language.Syntactic.Functional.Tuple
import qualified Language.Syntactic as Syn

-- imperative-edsl.
import Language.Embedded.Expression

--------------------------------------------------------------------------------
-- * Software expressions.
--------------------------------------------------------------------------------

type instance Pred SoftwareDomain = SoftwarePrimType

--------------------------------------------------------------------------------
-- hmm...

type instance ExprOf   (SExp a) = SExp
type instance PredOf   SExp     = SoftwarePrimType
type instance DomainOf SExp     = SoftwareDomain
type instance RepresentationOf SoftwarePrimType = SoftwarePrimTypeRep

--------------------------------------------------------------------------------
-- ** Software types.

-- | Representation of supported software types.
type STypeRep  = TypeRep SoftwarePrimType SoftwarePrimTypeRep

-- | ...
type STypeRepF = TypeRepF SoftwarePrimType SoftwarePrimTypeRep

instance Type SoftwarePrimType Bool   where typeRep = Node BoolST
instance Type SoftwarePrimType Int8   where typeRep = Node Int8ST
instance Type SoftwarePrimType Int16  where typeRep = Node Int16ST
instance Type SoftwarePrimType Int32  where typeRep = Node Int32ST
instance Type SoftwarePrimType Int64  where typeRep = Node Int64ST
instance Type SoftwarePrimType Word8  where typeRep = Node Word8ST
instance Type SoftwarePrimType Word16 where typeRep = Node Word16ST
instance Type SoftwarePrimType Word32 where typeRep = Node Word32ST
instance Type SoftwarePrimType Word64 where typeRep = Node Word64ST
instance Type SoftwarePrimType Float  where typeRep = Node FloatST
instance Type SoftwarePrimType Double where typeRep = Node DoubleST
instance Type SoftwarePrimType (Complex Float)  where typeRep = Node ComplexFloatST
instance Type SoftwarePrimType (Complex Double) where typeRep = Node ComplexDoubleST

-- | Compare two software types for equality.
softwareTypeEq :: STypeRep a -> STypeRep b -> Maybe (Dict (a ~ b))
softwareTypeEq (Node t)       (Node u) = softwarePrimTypeEq t u
softwareTypeEq (Branch t1 u1) (Branch t2 u2) = do
  Dict <- softwareTypeEq t1 t2
  Dict <- softwareTypeEq u1 u2
  return Dict
softwareTypeEq _ _ = Nothing

-- | Construct the software type representation of 'a'.
softwareTypeRep :: Struct SoftwarePrimType c a -> STypeRep a
softwareTypeRep = mapStruct (const softwareRep)

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

-- | Short-hand for software types.
type SType    = Type SoftwarePrimType

-- | Short-hand for primitive software types.
type SType'   = PrimType SoftwarePrimType

--------------------------------------------------------------------------------
-- ** Software expression symbols.

data AssertionLabel =
    -- ^ Internal assertion to guarantee meaningful results.
    InternalAssertion
    -- ^ User made assertion.
  | UserAssertion String
  deriving (Eq, Show)

-- | Guard a value with an assertion.
data GuardVal sig
  where
    GuardVal :: AssertionLabel -> String -> GuardVal (Bool :-> a :-> Full a)

deriving instance Eq       (GuardVal a)
deriving instance Show     (GuardVal a)
deriving instance Typeable (GuardVal a)

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

-- | Hint that a value may appear in an invariant
data HintVal sig
  where
    HintVal :: SoftwarePrimType a => HintVal (a :-> b :-> Full b)

deriving instance Eq       (HintVal a)
deriving instance Show     (HintVal a)
deriving instance Typeable (HintVal a)

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

-- | For loop.
data ForLoop sig
  where
    ForLoop :: SType st =>
        ForLoop (Length :-> Length :-> st :-> (Index -> st -> st) :-> Full st)

deriving instance Eq       (ForLoop a)
deriving instance Show     (ForLoop a)
deriving instance Typeable (ForLoop a)

--------------------------------------------------------------------------------
{-
-- |
data Foreign sig
  where
    Foreign :: Signature sig => String -> Denotation sig -> Foreign sig
-}
--------------------------------------------------------------------------------

-- | Software symbols.
type SoftwareConstructs =
          BindingT
  Syn.:+: Let
  Syn.:+: Tuple
  Syn.:+: SoftwarePrimConstructs
  Syn.:+: Construct
  -- ^ Software specific symbol.
  Syn.:+: GuardVal
  Syn.:+: HintVal
  Syn.:+: ForLoop

-- | Software symbols tagged with their type representation.
type SoftwareDomain = SoftwareConstructs :&: TypeRepF SoftwarePrimType SoftwarePrimTypeRep

-- | Software expressions.
newtype SExp a = SExp { unSExp :: ASTF SoftwareDomain a }

-- | Evaluate a closed software expression.
eval :: (Syntactic a, Domain a ~ SoftwareDomain) => a -> Internal a
eval = evalClosed . desugar

-- | Sugar a software symbol as a smart constructor.
sugarSymSoftware
  :: ( Syn.Signature sig
       , fi  ~ SmartFun dom sig
       , sig ~ SmartSig fi
       , dom ~ SmartSym fi
       , dom ~ SoftwareDomain
       , SyntacticN f fi
       , sub :<: SoftwareConstructs
       , Type SoftwarePrimType (DenResult sig)
       )
    => sub sig -> f
sugarSymSoftware = sugarSymDecor $ ValT $ typeRep

-- | Sugar a software symbol as a primitive smart constructor.
sugarSymPrimSoftware
    :: ( Syn.Signature sig
       , fi  ~ SmartFun dom sig
       , sig ~ SmartSig fi
       , dom ~ SmartSym fi
       , dom ~ SoftwareDomain
       , SyntacticN f fi
       , sub :<: SoftwareConstructs
       , SoftwarePrimType (DenResult sig)
       )
    => sub sig -> f
sugarSymPrimSoftware = sugarSymDecor $ ValT $ Node softwareRep

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

instance Syntactic (SExp a)
  where
    type Domain   (SExp a) = SoftwareDomain
    type Internal (SExp a) = a

    desugar = unSExp
    sugar   = SExp

instance Syntactic (Struct SoftwarePrimType SExp a)
  where
    type Domain   (Struct SoftwarePrimType SExp a) = SoftwareDomain
    type Internal (Struct SoftwarePrimType SExp a) = a

    desugar (Node a)     = unSExp a
    desugar (Branch a b) = sugarSymDecor (ValT $ Branch ta tb) Pair a' b'
      where
        a' = desugar a
        b' = desugar b
        ValT ta = getDecor a'
        ValT tb = getDecor b'

    sugar a = case getDecor a of
      ValT (Node _)       -> Node $ SExp a
      ValT (Branch ta tb) -> Branch
        (sugarSymDecor (ValT ta) Fst a)
        (sugarSymDecor (ValT tb) Snd a)
      FunT _ _ -> error "Syntactic can't sugar a function."

instance Tuples SoftwareDomain
  where
    pair   = sugarSymSoftware Pair
    first  = sugarSymSoftware Fst
    second = sugarSymSoftware Snd

instance FreeExp SExp
  where
    type FreePred SExp = PrimType SoftwarePrimType
    constExp = sugarSymSoftware . Lit
    varExp   = sugarSymSoftware . FreeVar

--------------------------------------------------------------------------------
-- syntactic instances.

instance Eval GuardVal
  where
    evalSym (GuardVal InternalAssertion msg) = \cond a ->
      if cond then a else error $ "Internal assertion failure: " ++ show msg
    evalSym (GuardVal (UserAssertion ass) msg) = \cond a ->
      if cond then a else error $ "User assertion " ++ show ass ++ " failure: " ++ show msg

instance Symbol GuardVal
  where
    symSig (GuardVal _ _) = signature

instance Render GuardVal
  where
    renderSym  = show
    renderArgs = renderArgsSmart

instance EvalEnv GuardVal env

instance StringTree GuardVal

instance Equality GuardVal

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

instance Eval HintVal
  where
    evalSym (HintVal) = flip const

instance Symbol HintVal
  where
    symSig (HintVal) = signature

instance Render HintVal
  where
    renderSym  = show
    renderArgs = renderArgsSmart

instance EvalEnv    HintVal env
instance StringTree HintVal
instance Equality   HintVal

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

instance Eval ForLoop
  where
    evalSym ForLoop = \min max init body ->
      foldl (flip body) init [min..max]

instance Symbol ForLoop
  where
    symSig (ForLoop) = signature

instance Render ForLoop
  where
    renderSym  = show
    renderArgs = renderArgsSmart

instance EvalEnv    ForLoop env
instance StringTree ForLoop
instance Equality   ForLoop

--------------------------------------------------------------------------------
-- *** Temporary fix until GHC fixes their class resolution for DTC ***

instance {-# OVERLAPPING #-} Project sub SoftwareConstructs =>
    Project sub (AST SoftwareDomain)
  where
    prj (Sym s) = Syn.prj s
    prj _ = Nothing

instance {-# OVERLAPPING #-} Project sub SoftwareConstructs =>
    Project sub SoftwareDomain
  where
    prj (expr :&: info) = Syn.prj expr

instance {-# OVERLAPPING #-} Project BindingT SoftwareConstructs
  where
    prj (InjL a) = Just a
    prj _ = Nothing

instance {-# OVERLAPPING #-} Project Let SoftwareConstructs
  where
    prj (InjR (InjL a)) = Just a
    prj _ = Nothing

instance {-# OVERLAPPING #-} Project Tuple SoftwareConstructs
  where
    prj (InjR (InjR (InjL a))) = Just a
    prj _ = Nothing

instance {-# OVERLAPPING #-} Project SoftwarePrimConstructs SoftwareConstructs
  where
    prj (InjR (InjR (InjR (InjL a)))) = Just a
    prj _ = Nothing

instance {-# OVERLAPPING #-} Project Construct SoftwareConstructs
  where
    prj (InjR (InjR (InjR (InjR (InjL a))))) = Just a
    prj _ = Nothing

instance {-# OVERLAPPING #-} Project GuardVal SoftwareConstructs
  where
    prj (InjR (InjR (InjR (InjR (InjR (InjL a)))))) = Just a
    prj _ = Nothing

instance {-# OVERLAPPING #-} Project HintVal SoftwareConstructs
  where
    prj (InjR (InjR (InjR (InjR (InjR (InjR (InjL a))))))) = Just a
    prj _ = Nothing

instance {-# OVERLAPPING #-} Project ForLoop SoftwareConstructs
  where
    prj (InjR (InjR (InjR (InjR (InjR (InjR (InjR a))))))) = Just a
    prj _ = Nothing

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

max_stars_repo_name: markus-git/co-feldspar

max_issues_repo_name: markus-git/co-feldspar

max_forks_repo_name: markus-git/co-feldspar


text:
#!/usr/bin/env stack
-- stack runghc --package reanimate
{-# LANGUAGE OverloadedStrings #-}
module Main where

import           Codec.Picture.Types
import           Control.Exception
import           Control.Lens                    ()
import           Control.Monad
import           Data.Function
import           Data.List
import           Data.List.NonEmpty              (NonEmpty)
import qualified Data.List.NonEmpty              as NE
import           Data.Maybe
import           Data.Ratio
import qualified Data.Text                       as T
import           Data.Tuple
import qualified Data.Vector                     as V
import           Debug.Trace
import           Linear.Matrix                   hiding (trace)
import           Linear.Metric
import           Linear.V2
import           Linear.V3
import           Linear.Vector
import           Numeric.LinearAlgebra           hiding (polar, scale, (<>))
import qualified Numeric.LinearAlgebra           as Matrix
import           Numeric.LinearAlgebra.HMatrix   hiding (polar, scale, (<>))
import           Reanimate
import           Reanimate.Builtin.Documentation
import           Reanimate.Debug
import           Reanimate.Math.Balloon
import           Reanimate.Math.Common
import           Reanimate.Math.Compatible       (compatiblyTriangulateP)
import qualified Reanimate.Math.DCEL             as DCEL
import           Reanimate.Math.EarClip
import           Reanimate.Math.Polygon
import           Reanimate.Math.Render
import           Reanimate.Math.Smooth
import           Reanimate.Math.SSSP
import           Reanimate.Math.Visibility
import           Reanimate.Morph.Common
import           Reanimate.Morph.LeastDifference
import           Reanimate.Morph.LeastWork
import           Reanimate.Morph.Linear
import           Reanimate.Morph.Rigid
import           Reanimate.PolyShape             (svgToPolygons)
import           Text.Printf

-- p1 = centerPolygon $ shape2
--p1 = pScale 3.5 $ pAtCenter $ pAddPoints (0+2) (pSetOffset shape13 0)
-- p1 = setOffset (addPoints 2 shape13) 0
-- p2 = scalePolygon 0.5 $ centerPolygon shape20
-- p1 = centerPolygon shape2
--p2 = pScale 3.5 $ pAtCenter $ pAddPoints 0 (pSetOffset shape14 0)
-- p2 = setOffset shape14 0
p1 = pCopy $ pSetOffset (pCopy p1') 0
p2 = pCopy $ pSetOffset (pCopy p2') 0
(p1', p2') =
  -- normalizePolygons
  -- closestLinearCorrespondence
  (,)
  -- leastWork zeroStretchCosts defaultBendCosts
  -- leastWork defaultStretchCosts defaultBendCosts
  -- (pAtCenter $ unsafeSVGToPolygon 0.01 $ scale 6 $ latex "S")
  -- (pAtCenter $ unsafeSVGToPolygon 0.01 $ scale 6 $ latex "C")
  (pAtCenter $ unsafeSVGToPolygon 0.1 $ scale 6 $ latex "X")
  (pAtCenter $ unsafeSVGToPolygon 0.1 $ scale 6 $ latex "I")

p1_ = castPolygon p1
p2_ = castPolygon p2
polys = triangulate_ p1_ p2_
p1_circ = castPolygon (circumference (map fst polys') p1_)
p2_circ = castPolygon (circumference (map snd polys') p2_)

polys' = alignPolygons polys p1_ p2_
(p1s, p2s) = unzip $ compatTriagPairs [ (castPolygon a, castPolygon b) | (a,b) <- polys' ]

-- (p1s,p2s) = unzip (compatiblyTriangulateP p1 p2)

m2 = DCEL.buildMesh $ DCEL.polygonsMesh
        (map (fmap realToFrac) $ V.toList $ polygonPoints p2_circ)
        (map (map (fmap realToFrac) . V.toList . polygonPoints) p2s)

m1 = DCEL.buildMesh $ DCEL.polygonsMesh
        (map (fmap realToFrac) $ V.toList $ polygonPoints p1_circ)
        (map (map (fmap realToFrac) . V.toList . polygonPoints) p1s)

pipeline1 = last . take 20 . iterate
  ( uncurry DCEL.delaunayFlip .
    uncurry DCEL.splitInternalEdges .
   (\(a,b) -> (DCEL.meshSmoothPosition a, DCEL.meshSmoothPosition b)))
(m1final, m2final) = last $ take 30 $ iterate
  (pipeline1 .
    uncurry DCEL.splitLongestEdge .
    pipeline1
    ) (m1,m2)

-- (m1good, m2good) = last $ take 9 $ iterate
--   (pipeline1 .
--     -- uncurry DCEL.splitLongestEdge .
--     pipeline1
--     ) (m1,m2)

main :: IO ()
-- main = do
--   evaluate (last $ take 900 $ compatiblyTriangulateP p1 p2)
--   return ()
-- main = reanimate $ playTraces $ length $ take 29 p1s
-- main = reanimate $ playThenReverseA $ pauseAround 0.5 0.5 $ scene $ do

main = reanimate $ scene $ do
  newSpriteSVG_ $ mkBackgroundPixel rtfdBackgroundColor
  -- play
  --   $ playTraces
  --   $ traceA (mapA (withFillOpacity 0 . withStrokeColor "white") drawCircle)
  --   $ traceA (mapA (withFillOpacity 0 . withStrokeColor "white") drawBox)
  --   $ traceSVG ((withFillOpacity 0 . withStrokeColor "white") $ mkCircle 2)
  --   $ traceSVG ((withFillOpacity 0 . withStrokeColor "white") $ mkRect 2 2)
  --   $ 20
  -- play $ playTraces $ last $ take 35 $ compatiblyTriangulateP p1 p2
  -- fork $ newSpriteA $ drawCompatible p1 p2
  -- newSpriteSVG_ $ translate 5 0 $ lowerTransformations $ scale 2 $
  --   -- DCEL.renderMesh 0.02 m1final
  --   DCEL.renderMeshColored m1final
    -- renderMeshPair optMesh
  -- newSpriteSVG_ $ translate (-5) 0 $ lowerTransformations $ scale 2 $
  --   DCEL.renderMesh 0.02 m1good

  -- newSpriteSVG_ $
  --   translate (-3) 0 $ scale 1 $ mkGroup
  --   [ mkGroup []
  --   , withFillColor "grey" $ polygonShape p1
  --   , withFillColor "grey" $ polygonNumDots p1
  --   ]
  -- newSpriteSVG_ $
  --   translate 3 0 $ scale 1 $ mkGroup
  --   [ mkGroup []
  --   , withFillColor "grey" $ polygonShape p2
  --   , withFillColor "grey" $ polygonNumDots p2
  --   ]
  -- wait 1
  -- fork $ forM_ p1s $ \p1Piece -> do
  --   newSpriteSVG_ $
  --     translate (-3) 0 $ scale 1 $ mkGroup
  --     [ mkGroup []
  --     , withFillColor "white" $ polygonShape p1Piece
  --     -- , withFillColor "grey" $ polygonNumDots p1Piece
  --     ]
  --   wait (1/60)
  --   -- wait 1
  -- fork $ forM_ p2s $ \p1Piece -> do
  --   newSpriteSVG_ $
  --     translate 3 0 $ scale 1 $ mkGroup
  --     [ mkGroup []
  --     , withFillColor "white" $ polygonShape p1Piece
  --     -- , withFillColor "grey" $ polygonNumDots p1Piece
  --     ]
  --   wait (1/60)
  --   -- wait 1

  -- fork $ play $ mkAnimation 3 $ \t ->
  --   let points = interpolate prep t
  --   in  translate (-3) 0 $ scale 1 $ mkGroup
  --   [ mkGroup []
  --   , drawTrigs points (meshTriangles myMesh)
  --   -- , DCEL.renderMeshColored m1
  --   ]
  fork $ play $ mkAnimation 3 $ \t ->
    let points = interpolate optPrep t
    in  translate 3 0 $ lowerTransformations $ scale 1 $ mkGroup
      [ mkGroup []
      , drawTrigs points (meshTriangles optMesh)
      -- , DCEL.renderMeshColored m1final
      -- , withGroupOpacity 0.5 $ DCEL.renderMesh 0.02 m1final
      -- , withGroupOpacity 0.5 $ withFillOpacity 0 $ renderMeshPair optMesh
      -- , let V2 x y = V.head points in
      --   translate x y $ withFillColor "green" $
      --   mkCircle 0.03
      ]
        -- , scale 2 $ polygonNumDots p1
  -- play $ pauseAtEnd 1 $ mkAnimation 3 $ \t ->
  --   let points = interpolate bestPrep t
  --   in  translate 3 0 $ mkGroup [drawTrigs points (meshTriangles bestMesh)
  --       -- , scale 2 $ polygonNumDots p2
  --                                                                           ]
  -- play $ pauseAtEnd 1 $ mkAnimation 3 $ \t ->
  --   let points = interpolate prepRev t in
  --   mkGroup
  --   [ drawTrigs points (meshTriangles myMeshRev)
  --   -- , polygonNumDots shape2
  --   ]
  -- forM_ (myMesh : smoothMesh myMesh) $ \newMesh -> do
  --   let (minAngA, minAngB) = meshMinAngle newMesh
  --   txt <- newSpriteSVG $ withFillColor "white" $
  --     mkGroup
  --     [translate 5 2 $  center $
  --       latex $ T.pack $ printf "Min: %.1f" (minAngA/pi*180)
  --     ,translate 5 1 $  center $
  --       latex $ T.pack $ printf "Min: %.1f" (minAngB/pi*180) ]
  --   s <- newSpriteSVG $ renderAMesh newMesh
  --   wait (recip 60)
  --   destroySprite s
  --   destroySprite txt
  -- newSpriteSVG $ translate (-3) 0 $ scale 2 $ renderMesh p1 (map fst trigs)
  -- newSpriteSVG $ translate 3 0 $ scale 2 $ renderMesh p2 (map snd trigs)
  -- wait 1
 where
  prep      = prepare myMesh
  optPrep   = prepare optMesh
  -- bestPrep  = prepare bestMesh
  myMesh    = mkMesh p1 p2
  optMesh   = toRigidMesh  m1final m2final
  -- optMesh   = DCEL.toRigidMesh  m1good m2good
  -- bestMesh  = last $ smoothMesh myMesh
  prepRev   = prepare myMeshRev
  myMeshRev = reverseMesh myMesh

reverseMesh :: Mesh -> Mesh
reverseMesh mesh =
  mesh { meshPointsA = meshPointsB mesh, meshPointsB = meshPointsA mesh }

mkMesh :: Polygon -> Polygon -> Mesh
mkMesh a b = Mesh
  { meshPointsA   = V.map (fmap realToFrac) pointsA
  , meshPointsB   = V.map (fmap realToFrac) pointsB
  , meshOutline = V.map (fromJust . flip V.elemIndex pointsA) (polygonPoints a)
  , meshTriangles = relTrigs
  }
 where
  pointsA =
    V.fromList $ nub $ V.toList $ V.concat [ polygonPoints a | (a, b) <- trigs ]
  pointsB =
    V.fromList $ nub $ V.toList $ V.concat [ polygonPoints b | (a, b) <- trigs ]
  trigs = compatiblyTriangulateP a b
  mkRel arr p =
    ( fromJust $ V.elemIndex (pAccess p 0) arr
    , fromJust $ V.elemIndex (pAccess p 1) arr
    , fromJust $ V.elemIndex (pAccess p 2) arr
    )
  relTrigs = V.fromList [ mkRel pointsA a | (a, b) <- trigs ]
    -- data Mesh = Mesh (Vector P) (Vector (RelTrig, RelTrig))

testMesh2 :: Mesh
testMesh2 = mkMesh (pTranslate (V2 0 0) shape2) shape20

testPrep = prepare testMesh2

drawTrigs :: V.Vector (V2 Double) -> V.Vector RelTrig -> SVG
drawTrigs points trigs = mkGroup
  [ mkGroup
      [ withFillOpacity 1
        $ withStrokeWidth (defaultStrokeWidth * 0)
        $ withStrokeColor "grey"
        $ withFillColor "black"
        $ drawPolygon
        $ map (points V.!) [a, b, c]
      | (a, b, c) <- V.toList trigs
      ]
    -- , a == 0
  -- , mkGroup
  --   [ withFillColor "red" $ mkGroup
  --     [ drawPoint (points V.! a), drawPoint (points V.! b), drawPoint (points V.! c)]
  --   | (a,b,c) <- V.toList trigs
  --   -- , a == 0
  --   ]
  ]

drawTrigsLines :: V.Vector (V2 Double) -> V.Vector RelTrig -> SVG
drawTrigsLines points trigs = withFillOpacity 0 $ withStrokeColor "black" $ mkGroup
  [ mkLinePathClosed
    [ (aPx, aPy)
    , (bPx, bPy)
    , (cPx, cPy)]
  | (a,b,c) <- V.toList trigs
  -- , a==2 || b == 2 || c == 2
  , let V2 aPx aPy = points V.! a
        V2 bPx bPy = points V.! b
        V2 cPx cPy = points V.! c
  ]

drawPoint :: V2 Double -> SVG
drawPoint (V2 x y) = translate x y $ mkCircle 0.1

drawPolygon :: [V2 Double] -> SVG
drawPolygon lst = mkLinePathClosed [ (x, y) | V2 x y <- lst ]

drawCompatible :: Polygon -> Polygon -> Animation
drawCompatible a b = scene $ do
  let left  = -6
      right = 4
  newSpriteSVG $ translate left 0 $ mkGroup
    [ withFillColor "grey" $ polygonShape a
    , withFillColor "grey" $ polygonNumDots a
    ]
  newSpriteSVG $ translate right 0 $ mkGroup
    [ withFillColor "grey" $ polygonShape b
    , withFillColor "grey" $ polygonNumDots b
    ]
  return ()
  -- let compat = compatiblyTriangulateP a b
  -- forM_ compat $ \(l, r) -> do
  --   fork $ play $ staticFrame 1 $
  --     translate left 0 $ withStrokeColor "white" $ withStrokeWidth (defaultStrokeWidth*0.2) $
  --     withFillOpacity 0 $ polygonShape l
  --   fork $ play $ staticFrame 1 $
  --     translate right 0 $ withStrokeColor "white" $ withStrokeWidth (defaultStrokeWidth*0.2) $
  --     withFillOpacity 0 $ polygonShape r

max_stars_repo_name: cdodev/reanimate

max_issues_repo_name: cdodev/reanimate

max_forks_repo_name: cdodev/reanimate


text:
{-# LANGUAGE ScopedTypeVariables #-}

module NaiveBayes where

import Data.List (maximumBy)
import Data.Ord (comparing)
import Mnist (LabeledData, groupByLabel, dat)
import Data.Function (on)
import Numeric.LinearAlgebra (Vector, Matrix, R, size, fromList, toColumns, sumElements, fromRows, (#>), toList)

-- TODO: later expand this by receiving label type and data type
data NaiveBayesClassifier = NBC
    { p :: !(Vector R)
    , l :: !(Matrix R)
    }

-- represents a classifier that predicts the class from input sample
class Classifier a where
    -- given the prior and likelihood, with a sample data, predict the class
    predict :: a -> Vector R -> Int  -- TODO: parameterize label type

instance Classifier NaiveBayesClassifier where
    predict nbc sample = argmax $ (p nbc) * (l nbc #> sample)

-- divide two integers to make a fraction
fractionDiv :: Fractional b => Int -> Int -> b
fractionDiv = (/) `on` fromIntegral

-- calculates the prior probability of a dataset
-- the data are organized as rows having same label
-- P(y)
prior :: LabeledData -> Vector R
prior ld = fromList $ map (flip fractionDiv totalNum) numElemsPerLabel
  where
    grouped :: [Matrix R] = groupByLabel ld
    numElemsPerLabel :: [Int] = map (fst . size) grouped
    totalNum = sum numElemsPerLabel

-- used as prior distribution of likelihood
epsilon :: R = 0.001

-- conditional likelihood of each data attributes given the class
-- P(x | y)
likelihood :: LabeledData -> Matrix R  -- data_class * data_dimension shape
likelihood ld = fromRows $ map attrProb colSumByLabel
  where
    grouped :: [Matrix R] = groupByLabel ld
    colSum :: Matrix R -> [R] = map sumElements . toColumns
    colSumByLabel :: [[R]] = map colSum grouped
    totalSum :: [R] = colSum $ dat ld
    -- create a vector of size DATA_DIMENSION that contains probability of each attributes
    attrProb classSum = fromList $ map (\ a -> fst a / (epsilon + snd a)) $ zip classSum totalSum

-- find maximum index in the vector
argmax :: Vector R -> Int
argmax = fst . maximumBy (comparing snd) . zip [0..] . toList

max_stars_repo_name: deNsuh/hasml

max_issues_repo_name: deNsuh/hasml

max_forks_repo_name: deNsuh/hasml


text:
{-#Language TypeFamilies#-}
module Morse where

import Stratification as S hiding (shiftCeil)
import Data.Matrix (Matrix(..))
import Control.Monad
import Control.Arrow
import SymbolicImage (
    fromRad, toRad, fromCeil, toCeil, myNub, mySequence, (#=))
import Jac (jac, e)
import Data.Array as A
import Graph
import Data.Function
import Data.List as L
import Data.Map as M
import Data.IntMap as IM
import Data.Complex
import Data
import Point
import Data.Graph

morse :: F Point -> Imagination -> Int -> [(Double, Double, Int)]
morse fp im i = mrsElem <$> cyclics
    ((shiftCeil d1 d2 fp #=) <$> ls)
  where
    mrsElem :: [(Ceil3, [(Ceil3, Double)])] -> (Double, Double, Int)
    mrsElem a = (
                         minOptZ gr  (minBazeCircuit gr),
                negate $ minOptZ gr' (maxBazeCircuit gr'),
                length a
--                       maximum $ plat a,
--                       minimum $ plat a
            )
        where
            gr  = formGraph                    a
            gr' = formGraph $ modifyIns negate a

    (Stratification d1 d2 ls) = str fp im i

plat :: [(Ceil3, [(Ceil3, Double)])] -> [Double]
plat = concatMap (\x -> snd <$> snd x)

-- XXX Дзен
modifyIns :: (a -> b) -> [(c, [(d, a)])] -> [(c, [(d, b)])]
modifyIns f l = second (second f <$>) <$> l

-- выделяем области связности
--cyclics :: (Ord t0, Ord k0) => [(k0, t0, [(t0, b0)])] -> [[(k0, [(t0, b0)])]]
cyclics g = M.toList <$> do
    let g' = (\(a, b, l) -> (a, b, L.filter (\x -> fst x /= a) l)) <$> g
    CyclicSCC a <- stronglyConnComp
        [(a, b, fst <$> c) | (a, b, c) <- g']
    return $ elect a $ M.fromList [(a, b) | (a, _, b) <- g']

shiftCeil :: Diameter -> Diameter -> F Point -> Ceil3 -> [(Ceil3, Double)]
shiftCeil d1 d2 fp (Ceil3 c l) = myNub $ do
    let Point x y = fromCeil d1 c
    ret <$> mySequence x d1
        <*> mySequence y d1
        <*> mySequence (rad d2 l) (pi/2/d2)
  where
    ret x' y' t' = (
        func d1 d2 fp x' y' t',
        log $ x $ abs $ toPoint $ jacobian fp x' y' * e t')

elect :: Ord a => [a] -> F (Map a b)
elect a m = M.intersection m $ M.fromList $ (\a -> (a,a)) <$> a

max_stars_repo_name: vojiranto/dsa

max_issues_repo_name: vojiranto/dsa

max_forks_repo_name: vojiranto/dsa


text:
{-# OPTIONS_GHC -Wall #-}

import qualified Numeric.Sundials.ARKode.ODE as ARK
import qualified Numeric.Sundials.CVode.ODE  as CV
import           Numeric.LinearAlgebra

import           Plots as P
import qualified Diagrams.Prelude as D
import           Diagrams.Backend.Rasterific

import           Control.Lens

import           Test.Hspec

lorenz :: Double -> [Double] -> [Double]
lorenz _t u = [ sigma * (y - x)
              , x * (rho - z) - y
              , x * y - beta * z
              ]
  where
    rho = 28.0
    sigma = 10.0
    beta = 8.0 / 3.0
    x = u !! 0
    y = u !! 1
    z = u !! 2

_lorenzJac :: Double -> Vector Double -> Matrix Double
_lorenzJac _t u = (3><3) [ (-sigma), rho - z, y
                        , sigma   , -1.0   , x
                        , 0.0     , (-x)   , (-beta)
                        ]
  where
    rho = 28.0
    sigma = 10.0
    beta = 8.0 / 3.0
    x = u ! 0
    y = u ! 1
    z = u ! 2

brusselator :: Double -> [Double] -> [Double]
brusselator _t x = [ a - (w + 1) * u + v * u * u
                   , w * u - v * u * u
                   , (b - w) / eps - w * u
                   ]
  where
    a = 1.0
    b = 3.5
    eps = 5.0e-6
    u = x !! 0
    v = x !! 1
    w = x !! 2

_brussJac :: Double -> Vector Double -> Matrix Double
_brussJac _t x = (3><3) [ (-(w + 1.0)) + 2.0 * u * v, w - 2.0 * u * v, (-w)
                       , u * u                     , (-(u * u))     , 0.0
                       , (-u)                      , u              , (-1.0) / eps - u
                       ]
  where
    y = toList x
    u = y !! 0
    v = y !! 1
    w = y !! 2
    eps = 5.0e-6

stiffish :: Double -> [Double] -> [Double]
stiffish t v = [ lamda * u + 1.0 / (1.0 + t * t) - lamda * atan t ]
  where
    lamda = -100.0
    u = v !! 0

stiffishV :: Double -> Vector Double -> Vector Double
stiffishV t v = fromList [ lamda * u + 1.0 / (1.0 + t * t) - lamda * atan t ]
  where
    lamda = -100.0
    u = v ! 0

_stiffJac :: Double -> Vector Double -> Matrix Double
_stiffJac _t _v = (1><1) [ lamda ]
  where
    lamda = -100.0

predatorPrey :: Double -> [Double] -> [Double]
predatorPrey _t v = [ x * a - b * x * y
                    , d * x * y - c * y - e * y * z
                    , (-f) * z + g * y * z
                    ]
  where
    x = v!!0
    y = v!!1
    z = v!!2
    a = 1.0
    b = 1.0
    c = 1.0
    d = 1.0
    e = 1.0
    f = 1.0
    g = 1.0

lSaxis :: [[Double]] -> P.Axis B D.V2 Double
lSaxis xs = P.r2Axis &~ do
  let ts = xs!!0
      us = xs!!1
      vs = xs!!2
      ws = xs!!3
  P.linePlot' $ zip ts us
  P.linePlot' $ zip ts vs
  P.linePlot' $ zip ts ws

kSaxis :: [(Double, Double)] -> P.Axis B D.V2 Double
kSaxis xs = P.r2Axis &~ do
  P.linePlot' xs

main :: IO ()
main = do

  let res1 = ARK.odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0])
  renderRasterific "diagrams/brusselator.png"
                   (D.dims2D 500.0 500.0)
                   (renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res1))

  let res1a = ARK.odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0])
  renderRasterific "diagrams/brusselatorA.png"
                   (D.dims2D 500.0 500.0)
                   (renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res1a))

  let res2 = ARK.odeSolve stiffish [0.0] (fromList [0.0, 0.1 .. 10.0])
  renderRasterific "diagrams/stiffish.png"
                   (D.dims2D 500.0 500.0)
                   (renderAxis $ kSaxis $ zip [0.0, 0.1 .. 10.0] (concat $ toLists res2))

  let res2a = ARK.odeSolveV (ARK.SDIRK_5_3_4') Nothing 1e-6 1e-10 stiffishV (fromList [0.0]) (fromList [0.0, 0.1 .. 10.0])

  let res2b = ARK.odeSolveV (ARK.TRBDF2_3_3_2') Nothing 1e-6 1e-10 stiffishV (fromList [0.0]) (fromList [0.0, 0.1 .. 10.0])

  let maxDiffA = maximum $ map abs $
                 zipWith (-) ((toLists $ tr res2a)!!0) ((toLists $ tr res2b)!!0)

  let res2c = CV.odeSolveV (CV.BDF) Nothing 1e-6 1e-10 stiffishV (fromList [0.0]) (fromList [0.0, 0.1 .. 10.0])

  let maxDiffB = maximum $ map abs $
                 zipWith (-) ((toLists $ tr res2a)!!0) ((toLists $ tr res2c)!!0)

  let maxDiffC = maximum $ map abs $
                 zipWith (-) ((toLists $ tr res2b)!!0) ((toLists $ tr res2c)!!0)

  let res3 = ARK.odeSolve lorenz [-5.0, -5.0, 1.0] (fromList [0.0, 0.01 .. 10.0])

  renderRasterific "diagrams/lorenz.png"
                   (D.dims2D 500.0 500.0)
                   (renderAxis $ kSaxis $ zip ((toLists $ tr res3)!!0) ((toLists $ tr res3)!!1))

  renderRasterific "diagrams/lorenz1.png"
                   (D.dims2D 500.0 500.0)
                   (renderAxis $ kSaxis $ zip ((toLists $ tr res3)!!0) ((toLists $ tr res3)!!2))

  renderRasterific "diagrams/lorenz2.png"
                   (D.dims2D 500.0 500.0)
                   (renderAxis $ kSaxis $ zip ((toLists $ tr res3)!!1) ((toLists $ tr res3)!!2))

  let res4 = CV.odeSolve predatorPrey [0.5, 1.0, 2.0] (fromList [0.0, 0.01 .. 10.0])

  renderRasterific "diagrams/predatorPrey.png"
                   (D.dims2D 500.0 500.0)
                   (renderAxis $ kSaxis $ zip ((toLists $ tr res4)!!0) ((toLists $ tr res4)!!1))

  renderRasterific "diagrams/predatorPrey1.png"
                   (D.dims2D 500.0 500.0)
                   (renderAxis $ kSaxis $ zip ((toLists $ tr res4)!!0) ((toLists $ tr res4)!!2))

  renderRasterific "diagrams/predatorPrey2.png"
                   (D.dims2D 500.0 500.0)
                   (renderAxis $ kSaxis $ zip ((toLists $ tr res4)!!1) ((toLists $ tr res4)!!2))

  let res4a = ARK.odeSolve predatorPrey [0.5, 1.0, 2.0] (fromList [0.0, 0.01 .. 10.0])

  let maxDiffPpA = maximum $ map abs $
                   zipWith (-) ((toLists $ tr res4)!!0) ((toLists $ tr res4a)!!0)

  hspec $ describe "Compare results" $ do
    it "for SDIRK_5_3_4' and TRBDF2_3_3_2'" $ maxDiffA < 1.0e-6
    it "for SDIRK_5_3_4' and BDF" $ maxDiffB < 1.0e-6
    it "for TRBDF2_3_3_2' and BDF" $ maxDiffC < 1.0e-6
    it "for CV and ARK for the Predator Prey model" $ maxDiffPpA < 1.0e-3

max_stars_repo_name: steinitznavican/hmatrix-sundials

max_issues_repo_name: steinitznavican/hmatrix-sundials

max_forks_repo_name: steinitznavican/hmatrix-sundials


text:
module Main where
import Data.NeuralNetwork.Backend.BLASHS.Utils
import qualified Data.Vector.Storable as SV
import Numeric.LinearAlgebra
import Test.QuickCheck (generate)
import Test.Utils
import Test.Gen

main = t3

t1 = do
  ms <- test_corr2_arr 0 ks mt
  mapM_ (putStrLn . show) ms
  putStrLn "Good ones are:"
  let m4s = map (\k -> good_corr2 0 k mt) ks
  mapM_ (putStrLn . show) m4s

t2 = do
  let m1 = (3><3) [0, 1, 0,
                   1, 0, 1,
                   0 ,1, 0] :: Matrix Float
  let ma = (3><3) [1, 0, 1,
                   0, 1, 0,
                   1 ,0, 1] :: Matrix Float
  let mb = (3><3) [0, 0, 0,
                   0, 0, 0,
                   0 ,0, 1] :: Matrix Float
  let mc = (3><3) [0, 0, 0,
                   0, 1, 0,
                   0 ,0, 0] :: Matrix Float
  let md = (3><3) [1, 0, 0,
                   0, 0, 0,
                   0 ,0, 0] :: Matrix Float
  let m3 = (4><4) [4,3,2,1,
                   1,2,3,4,
                   2,4,6,8,
                   1,3,5,7] :: Matrix Float
  m3s <- test_corr2_arr 0 [m1, ma, mb, mc, md] m3
  mapM_ (putStrLn . show) m3s

t3 = do
   ks <- generate $ sequence $ replicate 4 $ squared_real_matrices 45
   mt <- generate $ squared_real_matrices 488
   test_corr2_arr 2 ks mt

ks = [(3><3)  [  -2.2182608,    -2.65496,   -4.452847
              , -0.50510013,   1.1137468,   0.5881784
              ,   -2.162154, -0.26882064, -0.21511263 ]
     ,(3><3)  [ -168.82103, -235.06078, 4.6464214
              ,  3.0136368,  0.9492494, 605.85254
              ,  -482.5887,  14.835919, -5.846858 ]
     ,(3><3)  [ -1.3539604,    8.014306, -4.851793
              , 0.84038895, -0.77802753, 1.4430395
              ,  0.5361183,  -1.2853559,   1.04196 ]
     ,(3><3)  [ -2.6642451, -2.149097, -4.7066493
              ,   1.040131,   1.52895,  -5.823082
              ,  -5.870729,  2.996159,  1.2088959 ]
     ,(3><3)  [  -2.5846646,   2.2002325, -2.7979004
              ,   0.5225487,   0.9918108, 0.85764366
              , -0.56938225, -0.14684846,  0.8047059 ]]
mt = (7><7)   [  3.4427254,    0.7416965,   1.2416595,  0.8233356,  -1.565333, 0.89987564,  0.8893491
              , -1.3306698,   -3.9727004,  -1.4404364, 8.25734e-2,  1.0835173, -6.5956297, 0.95600057
              ,  -5.312902, -9.969063e-2,  0.94423944, 0.25942376,  -1.177719, -219.11214,  0.7662355
              , -3.0743675,   0.16589901,   -1.338092,  1.2396933, -2.8903048, -0.2627482, 0.64401335
              ,  6.0397525,   -0.6807962, -0.36770278,   1.416498,  1.8280395, -0.5226003, 0.98356235
              , -1.1980977,   -0.5414966, -0.18184663, 0.20523632,  1.2170252,  1.2050351,  1.1030668
              ,   0.673801,   0.95085037,  -2.3031695,  0.4833055,  1.8391184,   1.117064, -1.7435362 ]

max_stars_repo_name: pierric/neural-network

max_issues_repo_name: pierric/neural-network

max_forks_repo_name: pierric/neural-network


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

{-# LANGUAGE BangPatterns        #-}
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE GADTs               #-}
{-# LANGUAGE RankNTypes          #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeOperators       #-}

module Grenade.Utils.LinearAlgebra
    ( bmean
    , bvar
    , bvar'
    , vsqrt
    , msqrt
    , sreshape
    , sflatten
    , vscale
    , vadd
    , vflatten
    , nsum
    , sumV
    , sumM
    , squareV
    , squareM
    , extractV
    , extractM2D
    , batchNormMean
    , batchNormVariance
    , vectorToList
    , listToVector
    ) where

import           Data.List                    (foldl')
import           Data.Maybe                   (fromJust)
import           Data.Proxy
import           Data.Singletons

import           GHC.TypeLits

import           Numeric.LinearAlgebra        hiding (R)
import           Numeric.LinearAlgebra.Static (L, R)
import qualified Numeric.LinearAlgebra.Static as H

import           Grenade.Core.Shape
import           Grenade.Types

sumV :: (KnownNat n) => R n -> RealNum
sumV v = v H.<.> 1

sumM :: (KnownNat m, KnownNat n) => L m n -> RealNum
sumM m = (m H.#> 1) H.<.> 1

squareV :: (KnownNat n) => R n -> R n
squareV = H.dvmap (^ (2 :: Int))

squareM :: (KnownNat m, KnownNat n) => L m n ->  L m n
squareM = H.dmmap (^ (2 :: Int))

-- | Helper function that sums the elements of a matrix
nsum :: S s -> RealNum
nsum (S1D x) = sumElements $ H.extract x
nsum (S2D x) = sumElements $ H.extract x
nsum (S3D x) = sumElements $ H.extract x
nsum (S4D x) = sumElements $ H.extract x

-- | Calculate elementwise mean of list of matrices
bmean :: forall s. SingI s => [S s] -> S s
bmean xs
  = let (!m, !l) = foldl' (\(m, l) x -> (m + x, l + 1)) (0 :: S s, 0) xs :: (S s, Int)
    in  case (m, l) of
          (S1D x, _) -> S1D $ H.dvmap (/ fromIntegral l) x
          (S2D x, _) -> S2D $ H.dmmap (/ fromIntegral l) x
          (S3D x, _) -> S3D $ H.dmmap (/ fromIntegral l) x
          (S4D x, _) -> S4D $ H.dmmap (/ fromIntegral l) x

-- | Calculate element wise variance
bvar :: forall s. SingI s => [S s] -> S s
bvar xs
  = let !m       = bmean xs
        (!v, !l) = foldl' (\(v, l) x -> (v + (x - m)**2, l + 1)) (0 :: S s, 0) xs :: (S s, Int)
    in  case (v, l) of
          (S1D x, _) -> S1D $ H.dvmap (/ fromIntegral l) x
          (S2D x, _) -> S2D $ H.dmmap (/ fromIntegral l) x
          (S3D x, _) -> S3D $ H.dmmap (/ fromIntegral l) x
          (S4D x, _) -> S4D $ H.dmmap (/ fromIntegral l) x

-- | Calculate element wise variance, with the mean precalculated
bvar' :: forall s. SingI s => S s -> [S s] -> S s
bvar' m xs
  = let (!v, !l) = foldl' (\(v, l) x -> (v + (x - m)**2, l + 1)) (0 :: S s, 0) xs :: (S s, Int)
    in  case (v, l) of
          (S1D x, _) -> S1D $ H.dvmap (/ fromIntegral l) x
          (S2D x, _) -> S2D $ H.dmmap (/ fromIntegral l) x
          (S3D x, _) -> S3D $ H.dmmap (/ fromIntegral l) x
          (S4D x, _) -> S4D $ H.dmmap (/ fromIntegral l) x

batchNormMean :: forall n. KnownNat n => [R n] -> RealNum
batchNormMean vs
  = let hw  = fromIntegral $ natVal (Proxy :: Proxy n)
        vs' = map (sumElements . H.extract) vs :: [RealNum]
        n   = fromIntegral $ length vs
    in  sum vs' / (hw * n)

batchNormVariance :: forall n. KnownNat n => [R n] -> RealNum
batchNormVariance vs
  = let mu  = batchNormMean vs
        hw  = fromIntegral $ natVal (Proxy :: Proxy n)
        vs' = map (sumElements . H.extract . H.dvmap (\a -> (a - mu) ^ (2 :: Int))) vs :: [RealNum]
        n   = fromIntegral $ length vs
    in  sum vs' / (hw * n)

extractV :: S ('D1 x) -> R x
extractV (S1D v) = v

extractM2D :: S ('D2 x y) -> L x y
extractM2D (S2D m) = m

vscale :: KnownNat n => RealNum -> R n -> R n
vscale = H.dvmap . (*)

vadd :: KnownNat n => RealNum -> R n -> R n
vadd = H.dvmap . (+)

vsqrt :: KnownNat n => R n -> R n
vsqrt = H.dvmap sqrt

msqrt :: (KnownNat n, KnownNat m) => L m n -> L m n
msqrt = H.dmmap sqrt

vflatten :: (KnownNat m, KnownNat n) => [R n] -> R m
vflatten xs = fromJust . H.create . flatten . fromRows $ map H.extract xs

sflatten :: (KnownNat m, KnownNat n) => L m n -> R (m * n)
sflatten = fromJust . H.create . flatten . H.extract

sreshape :: forall m n. (KnownNat m, KnownNat n) => R (m * n) -> L m n
sreshape v
  = let rows = fromIntegral $ natVal (Proxy :: Proxy m)
    in  fromJust . H.create . reshape rows . H.extract $ v

vectorToList :: KnownNat n => R n -> [RealNum]
vectorToList = toList . H.extract

listToVector :: KnownNat n => [RealNum] -> R n
listToVector = fromJust . H.create . fromList

-- bsqrt :: SingI s => [S s] -> [S s]
-- bsqrt xs = map msqrt xs

-- msqrt :: SingI s => S s -> S s
-- msqrt (S1D x) = S1D $ dvmap sqrt x
-- msqrt (S2D x) = S2D $ dmmap sqrt x
-- msqrt (S3D x) = S3D $ dmmap sqrt x

-- mscale :: SingI s => RealNum -> S s -> S s
-- mscale r (S1D x) = S1D $ dvmap (* r) x
-- mscale r (S2D x) = S2D $ dmmap (* r) x
-- mscale r (S3D x) = S3D $ dmmap (* r) x

-- madd :: SingI s => RealNum -> S s -> S s
-- madd r (S1D x) = S1D $ dvmap (+ r) x
-- madd r (S2D x) = S2D $ dmmap (+ r) x
-- madd r (S3D x) = S3D $ dmmap (+ r) x

max_stars_repo_name: th-char/grenade

max_issues_repo_name: th-char/grenade

max_forks_repo_name: th-char/grenade


text:
-- | see examples/
module DimMat (

   -- * Quasiquotes
   matD,
   blockD,

   -- * Data.Packed.Vector
   (@>),
   -- * Data.Packed.Matrix
   -- ** dimension
   cols, rows,
   colsNT, rowsNT,
   hasRows, hasCols,
   -- (><),
   trans,
   -- reshape, flatten, fromLists, toLists, buildMatrix,
   -- broken
   ToHLists(toHLists), toHList, FromHLists(fromHLists), fromHList,
   (@@>),

   -- asRow, asColumn, fromRows, toRows, fromColumns, toColumns
   -- fromBlocks
#if MIN_VERSION_hmatrix(0,15,0)
   diagBlock,
#endif
   -- toBlocks, toBlocksEvery, repmat, flipud, fliprl
   -- subMatrix, takeRows, dropRows, takeColumns, dropColumns,
   -- extractRows, diagRect, takeDiag, mapMatrix,
   -- mapMatrixWithIndexM, mapMatrixWithIndexM_, liftMatrix,
   -- liftMatrix2, liftMatrix2Auto, fromArray2D,

   ident, -- where to put this?
   -- * Numeric.Container
   -- constant, linspace,
   diag,
   ctrans,
   -- ** Container class
   scalar,
   conj,
   scale, scaleRecip,
   recipMat,
   addConstant,
   add,
   sub,
   mulMat, mulVec,
   divideMat, divideVec,
   equal,
   arctan2,
   hconcat,
   vconcat,
   cmap,
   konst,
   zeroes,
   -- build, atIndex, minIndex, maxIndex, minElement, maxElement,
   -- sumElements, prodElements, step, cond, find, assoc, accum,
   -- Convert
   -- ** Product class
   Dot(..),
   -- absSum, norm1, norm2, normInf,
   -- norm1, normInf,
   pnorm,
   -- optimiseMult, mXm, mXv, vXm, (<.>),
   -- (<>), (<\>), outer, kronecker,
   -- ** Random numbers
   -- ** Element conversion
   -- ** Input/Output
   -- ** Experimental

   -- * Numeric.LinearAlgebra.Algorithms
   -- | incomplete wrapper for "Numeric.LinearAlgebra.Algorithms"

   -- ** Linear Systems
   -- linearSolve, luSolve, cholSolve, linearSolveLS, linearSolveSVD,
   inv,
   PInv(pinv),
   pinvTol,
   det,
   -- invlndet,
   rank,
   -- rcond,
   -- ** Matrix factorizations

   -- *** Singular value decomposition
   -- *** Eigensystems
   -- eigs
   {-
   wrapEig, wrapEigOnly,
   EigV, EigE,
   -- **** eigenvalues and eigenvectors
   eig,
   eigC,
   eigH,
   eigH',
   eigR,
   eigS,
   eigS',
   eigSH,
   eigSH',

   -- **** eigenvalues
   eigOnlyC,
   eigOnlyH,
   eigOnlyR,
   eigOnlyS,
   eigenvalues,
   eigenvaluesSH,
   eigenvaluesSH',
   -}

   -- *** QR
   -- *** Cholesky
   -- *** Hessenberg
   -- *** Schur
   -- *** LU

   -- ** Matrix functions
   -- sqrtm, matFunc
   expm,

   -- ** Nullspace
   -- ** Norms
   -- ** Misc
   -- ** Util

   -- * Automatic Differentiation
   -- ad
#ifdef WITH_Ad
   diff,
#endif

    -- * todo arrange
    DotS,
    MultEq,

    -- * to keep types looking ok
    D,
    module Data.HList.CommonMain,
    Complex,

    -- ** "Numeric.NumType"
    Pos, Neg, Succ, Zero, Neg1,
    ) where

import Numeric.NumType
import DimMat.Internal
import DimMat.QQ
import Data.HList.CommonMain
import Data.Complex

max_stars_repo_name: aavogt/DimMat

max_issues_repo_name: aavogt/DimMat

max_forks_repo_name: aavogt/DimMat