text:
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE TypeOperators    #-}
module FokkerPlanck.Pinwheel
  ( module Filter.Pinwheel
  , module FokkerPlanck.Pinwheel
  ) where

import           Data.Array.Repa               as R
import           Data.Complex
import           Data.List                     as L
import           Data.Vector.Unboxed           as VU
import           Filter.Pinwheel
import           FokkerPlanck.Interpolation
import           Numeric.LinearAlgebra.Data    as NL
import           Numeric.LinearAlgebra.HMatrix
import           Types
import           Utils.Array
import           Utils.Coordinates
import           Utils.Parallel

{-# INLINE computeR2Z1T0ArrayRadial #-}
computeR2Z1T0ArrayRadial ::
     (R.Source r Double)
  => R.Array r DIM3 Double
  -> Int
  -> Int
  -> Double
  -> [Double]
  -> [Double]
  -> R.Array D DIM4 (Complex Double)
computeR2Z1T0ArrayRadial radialArr xLen yLen scaleFactor thetaFreqs theta0Freqs =
  let pinwheelArr =
        traverse2
          (fromListUnboxed (Z :. L.length thetaFreqs) thetaFreqs)
          (fromListUnboxed (Z :. L.length theta0Freqs) theta0Freqs)
          (\(Z :. numThetaFreq) (Z :. numTheta0Freq) ->
             (Z :. numThetaFreq :. numTheta0Freq :. xLen :. yLen)) $ \f f0 (Z :. k :. l :. i :. j) ->
          pinwheel
            (f (Z :. k) - f0 (Z :. l))
            0
            (exp 1)
            0
            (i - center xLen)
            (j - center yLen)
   in radialCubicInterpolation radialArr scaleFactor pinwheelArr

{-# INLINE computeR2Z2T0S0ArrayRadial #-}
computeR2Z2T0S0ArrayRadial ::
     (R.Source r Double)
  => PinwheelType
  -> R.Array r DIM5 Double
  -> Int
  -> Int
  -> Double
  -> Double
  -> [Double]
  -> [Double]
  -> [Double]
  -> [Double]
  -> R.Array D DIM6 (Complex Double)
computeR2Z2T0S0ArrayRadial pinwheelType radialArr xLen yLen scaleFactor rMax thetaFreqs scaleFreqs theta0Freqs scale0Freqs =
  let pinwheelArr =
        traverse4
          (fromListUnboxed (Z :. L.length thetaFreqs) thetaFreqs)
          (fromListUnboxed (Z :. L.length scaleFreqs) scaleFreqs)
          (fromListUnboxed (Z :. L.length theta0Freqs) theta0Freqs)
          (fromListUnboxed (Z :. L.length scale0Freqs) scale0Freqs)
          (\(Z :. numThetaFreq) (Z :. numScaleFreq) (Z :. numTheta0Freq) (Z :. numScale0Freq) ->
             (Z :. numThetaFreq :. numScaleFreq :. numTheta0Freq :.
              numScale0Freq :.
              xLen :.
              yLen)) $ \ft fs ft0 fs0 (Z :. t :. s :. t0 :. s0 :. i :. j) ->
          pinwheelFunc
            pinwheelType
            (ft (Z :. t) - ft0 (Z :. t0))
            (fs (Z :. s) - fs0 (Z :. s0))
            rMax
            0
            (i - center xLen)
            (j - center yLen)
  in radialCubicInterpolation radialArr scaleFactor pinwheelArr

{-# INLINE computeR2Z2T0S0ArrayRadial' #-}
computeR2Z2T0S0ArrayRadial' ::
     Int
  -> Int
  -> Double
  -> Double
  -> [Double]
  -> [Double]
  -> [Double]
  -> [Double]
  -> R.Array D DIM6 (Complex Double)
computeR2Z2T0S0ArrayRadial' xLen yLen scaleFactor rMax thetaFreqs scaleFreqs theta0Freqs scale0Freqs =
  let pinwheelArr =
        traverse4
          (fromListUnboxed (Z :. L.length thetaFreqs) thetaFreqs)
          (fromListUnboxed (Z :. L.length scaleFreqs) scaleFreqs)
          (fromListUnboxed (Z :. L.length theta0Freqs) theta0Freqs)
          (fromListUnboxed (Z :. L.length scale0Freqs) scale0Freqs)
          (\(Z :. numThetaFreq) (Z :. numScaleFreq) (Z :. numTheta0Freq) (Z :. numScale0Freq) ->
             (Z :. numThetaFreq :. numScaleFreq :. numTheta0Freq :.
              numScale0Freq :.
              xLen :.
              yLen)) $ \ft fs ft0 fs0 (Z :. t :. s :. t0 :. s0 :. i :. j) ->
          pinwheel
            (ft (Z :. t) - ft0 (Z :. t0))
            (fs (Z :. s) - fs0 (Z :. s0))
            rMax
            0
            (i - center xLen)
            (j - center yLen)
   in pinwheelArr

{-# INLINE cutoff #-}
cutoff :: (R.Source r Double) => Int -> R.Array r DIM5 Double -> R.Array D DIM5 Double
cutoff r arr =
  R.traverse arr id $ \f idx@(Z :. _ :. _ :. _ :. _ :. e) ->
    if e > r
      then 0
      else f idx

{-# INLINE computeLocalEigenVector #-}
computeLocalEigenVector ::
     (R.Source r Double)
  => ParallelParams
  -> (Double -> Double -> Double -> Double -> Int -> Int -> Complex Double)
  -> R.Array r DIM5 Double
  -> Int
  -> Int
  -> Double
  -> [Double]
  -> [Double]
  -> R.Array U DIM4 (Complex Double)
computeLocalEigenVector parallelParams pinwheelFunc radialArr xLen yLen rMax thetaFreqs scaleFreqs =
  let (Z :. numThetaFreq :. numScaleFreq :. _ :. _ :. _) = extent radialArr
      len = numThetaFreq * numScaleFreq
   in computeS .
      rotate4D .
      rotate4D .
      fromUnboxed (Z :. xLen :. yLen :. numThetaFreq :. numScaleFreq) .
      VU.concat .
      parMapChunk
        parallelParams
        rdeepseq
        (\(x, y) ->
           let r =
                 sqrt $
                 (fromIntegral $ x - center xLen) ^ 2 +
                 (fromIntegral $ y - center yLen) ^ 2
               interpolatedArray = cubicInterpolation radialArr r
               interpolatedPinwheel =
                 R.traverse3
                   interpolatedArray
                   (fromListUnboxed (Z :. numThetaFreq) thetaFreqs)
                   (fromListUnboxed (Z :. numScaleFreq) scaleFreqs)
                   (\sh _ _ -> sh) $ \f ft fs idx@(Z :. t :. s :. t0 :. s0) ->
                   (f idx :+ 0) *
                   pinwheelFunc
                     (ft (Z :. t) - ft (Z :. t0))
                     (fs (Z :. s) + fs (Z :. s0))
                     rMax
                     0
                     (x - center xLen)
                     (y - center yLen)
               (eigVal, eigVec) =
                 eig . (len >< len) . R.toList $ interpolatedPinwheel
               (val, vec) =
                 L.head . L.reverse . L.sortOn (magnitude . fst) $
                 L.zip (NL.toList eigVal) (toColumns $ eigVec)
               dominantVec = VU.fromList . L.map (* val) . NL.toList $ vec
            in -- VU.zipWith
               --   (*)
               --   (VU.fromList
               --      [ exp
               --        (0 :+
               --         (tf) *
               --         (angleFunctionRad
               --            (fromIntegral $ x - center xLen)
               --            (fromIntegral $ y - center yLen)))
               --      | tf <- thetaFreqs
               --      , sf <- scaleFreqs
               --      ])
               --   dominantVec
               dominantVec
        ) $
      [(x, y) | x <- [0 .. xLen - 1], y <- [0 .. yLen - 1]]

{-# INLINE computeLocalEigenVectorSink #-}
computeLocalEigenVectorSink ::
     (R.Source r Double)
  => ParallelParams
  -> (Double -> Double -> Double -> Double -> Int -> Int -> Complex Double)
  -> R.Array r DIM5 Double
  -> Int
  -> Int
  -> Double
  -> [Double]
  -> [Double]
  -> R.Array U DIM4 (Complex Double)
computeLocalEigenVectorSink parallelParams pinwheelFunc radialArr xLen yLen rMax thetaFreqs scaleFreqs =
  let (Z :. numThetaFreq :. numScaleFreq :. _ :. _ :. _) = extent radialArr
      len = numThetaFreq * numScaleFreq
   in computeS .
      rotate4D .
      rotate4D .
      fromUnboxed (Z :. xLen :. yLen :. numThetaFreq :. numScaleFreq) .
      VU.concat .
      parMapChunk
        parallelParams
        rdeepseq
        (\(x, y) ->
           let r =
                 sqrt $
                 (fromIntegral $ x - center xLen) ^ 2 +
                 (fromIntegral $ y - center yLen) ^ 2
               interpolatedArray = cubicInterpolation radialArr r
               interpolatedPinwheel =
                 R.traverse3
                   interpolatedArray
                   (fromListUnboxed (Z :. numThetaFreq) thetaFreqs)
                   (fromListUnboxed (Z :. numScaleFreq) scaleFreqs)
                   (\sh _ _ -> sh) $ \f ft fs idx@(Z :. t :. s :. t0 :. s0) ->
                   (f idx :+ 0) * (exp (0 :+ (ft (Z :. t) - ft (Z :. t0)) * pi)) *
                   pinwheelFunc
                     (ft (Z :. t) - ft (Z :. t0))
                     (fs (Z :. s) + fs (Z :. s0))
                     rMax
                     0
                     (x - center xLen)
                     (y - center yLen)
               (eigVal, eigVec) =
                 eig . (len >< len) . R.toList $ interpolatedPinwheel
               (val, vec) =
                 L.head . L.reverse . L.sortOn (realPart . fst) $
                 L.zip (NL.toList eigVal) (toColumns $ eigVec)
            in VU.fromList . L.map (* val) . NL.toList $ vec) $
      [(x, y) | x <- [0 .. xLen - 1], y <- [0 .. yLen - 1]]

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:
module Util where

import Data.List
import Data.Ord

import Numeric.Container
import Numeric.LinearAlgebra ()

-- | Short function for logarithmus dualis
ld :: Floating a => a -> a
ld = logBase 2

-- | Returns a subsequnce of length k starting at index i.
cut :: Int -> Int -> [a] -> [a]
cut k i = take k . drop (i - 1)

-- | Returns all subsequences of length k.
cuts :: Int -> [a] -> [[a]]
cuts k s = map (take k) . zipWith const (iterate tail s) $ drop (k-1) s

-- | Checks two consecutive for convergence. If a convergence is found, the list
-- will be truncated.
convergenceOn :: (a -> a -> Double) -> Double -> [a] -> [a]
convergenceOn f c ls = map fst . takeWhile snd $ zip ls bs
  where
    bs = True : True : zipWith condition ls (tail ls)
    condition x y = f x y > c

-- | Returns the column indices for the biggest element in each row.
maxCols :: Matrix Double -> [Int]
maxCols = map row . toRows
  where
    row = fst . maximumBy (comparing snd) . zip [1 :: Int ..] . toList

-- | Normalizes the elements of a vector so that they sum up to 1.
normalizeVector :: Vector Double -> Vector Double
normalizeVector v = scale (recip $ sumElements v) v

-- | Normalizes the columns, so that each of column sums is 1.
normalizeColumns :: Matrix Double -> Matrix Double
normalizeColumns = fromColumns . map normalizeVector . toColumns

-- | Normalizes the row, so that each of row sums is 1.
normalizeRows :: Matrix Double -> Matrix Double
normalizeRows = fromRows . map normalizeVector . toRows

-- | Calculates the mean squared error for two matrices (elementwise + summing up).
meanSquaredError :: Matrix Double -> Matrix Double -> Double
meanSquaredError m m' = (/ n) . sumElements $ (m' - m)^^(2 :: Int)
  where
    n = fromIntegral $ rows m * cols m

-- vim: set ts=4 sts=4 sw=4 tw=80 et:

ext: hs

max_stars_repo_name: fgrsnau/emgene

max_issues_repo_name: fgrsnau/emgene

max_forks_repo_name: fgrsnau/emgene


text:
{-# LANGUAGE OverloadedStrings, LambdaCase, GeneralisedNewtypeDeriving #-}
-- | Module Simulation
--
-- Simulate races using the performance model.
module Simulation
    ( SimPip(..)
    , SimEntry(..)
    , toSimPips
    , simulateSingleRace
    , runExperimentFull
    , runExperimentPip
    , runExperimentProbe
    , simModelStrength
    , simAveragePositions
    , testSeed1
    ) where

import Types
import Orbital (orbitalDistr, OrbitalDistribution)
import Util.Lone (surroundL)

import System.Random.MWC
import Statistics.Distribution

import Data.Bifunctor
import qualified Data.Map.Strict as Map
import Data.Map.Strict (Map)
import Data.List
import Data.Ord
import Pipes
import qualified Pipes.Prelude as P
import Data.Default.Class
import Control.Monad.State.Strict
import qualified Data.Vector as V (fromList)

-- | Racer identification. The type is augmented by testing probes.
data SimPip = Probe !Int | SimPip !PipId
    deriving (Eq, Ord, Show)

-- | Racer result with rank generated by the simulation.
data SimEntry = SimEntry
    { simPip :: !SimPip
    , laptimeDelta :: !Double
    , pipRank :: !Int
    }
    deriving (Eq, Ord, Show)

-- | Prepare ratings from the engine to feed the simulation.
toSimPips
    :: Int                              -- ^ Gamma shape parameter.
    -> [(Int, Double)]                  -- ^ Identifications and ratings of the probes.
    -> Ratings                          -- ^ Ratings map from the engine.
    -> [(SimPip, OrbitalDistribution)]  -- ^ Identifications and performance models.
toSimPips gsh probes rtgs = map (second (orbitalDistr gsh))
    $ map (first Probe) probes
    ++ map (bimap SimPip rating) (Map.assocs rtgs)

-- | Infinite stream of simulated race results. If the number of runs that
-- will be used is known in advance, it is slightly less efficient than
-- 'simulatorN', as it has to save the generator state after each run.
simulator
    :: [(SimPip, OrbitalDistribution)]   -- ^ Identifications and performance models.
    -> Producer [SimEntry] SimM ()       -- ^ Identifications and laptime deltas.
simulator pips = do
    seed <- get
    g <- liftIO $ restore seed
    P.repeatM $ do
        results <- liftIO $ traverse (surroundL (flip genContVar g)) pips
        newSeed <- liftIO $ save g
        put newSeed
        return $ zipWith (uncurry SimEntry) (sortBy (comparing snd) results) [1..]

-- | Stream of a finite number of simulated results. The updated generator
-- seed is used as the return value to make it harder to discard it
-- accidentally by interrupting the stream. For interrupting the stream, use
-- 'void' to discard the return value, or switch to 'simulator' if it is
-- important to update the generator state.
simulatorN
    :: Int                               -- ^ Number of simulations
    -> [(SimPip, OrbitalDistribution)]   -- ^ Identifications and performance models.
    -> Producer [SimEntry] SimM Seed     -- ^ Identifications and laptime deltas.
simulatorN nRuns pips = do
    seed <- get
    g <- liftIO $ restore seed
    P.replicateM nRuns $ do
        results <- liftIO $ traverse (surroundL (flip genContVar g)) pips
        return $ zipWith (uncurry SimEntry) (sortBy (comparing snd) results) [1..]
    newSeed <- liftIO $ save g
    put newSeed
    return newSeed

-- | Generates results of a single race.
simulateSingleRace
    :: [(SimPip, OrbitalDistribution)]  -- ^ Identifications and performance models.
    -> SimM [SimEntry]                  -- ^ Identifications and laptime deltas.
simulateSingleRace pips = do
    P.head (void $ simulatorN 1 pips) >>= \case
        Just res -> return res
        Nothing -> error "Analysis.Simulation: simulator exhausted"

-- | Generates results of many races and tallies how often racers reached each
-- position.
runExperimentFull
    :: Int                              -- ^ Number of runs.
    -> [(SimPip, OrbitalDistribution)]  -- ^ Identifications and performance models.
    -> SimM (Map (SimPip, Int) Int)     -- ^ Count of racer-rank pairs.
runExperimentFull nRuns pips = fst <$> P.fold' updateCount counters id entries
    where
    counters = Map.empty
    entries = simulatorN nRuns pips >-> P.mapFoldable id
    updateCount ctrs se =
        Map.alter (maybe (Just 1) (Just . (+ 1))) (simPip se, pipRank se) ctrs

-- | Generates results of many races and tallies how often a specific racer
-- reached each position.
runExperimentPip
    :: Int                               -- ^ Number of runs.
    -> SimPip                            -- ^ Selected racer.
    -> [(SimPip, OrbitalDistribution)]   -- ^ Identifications and performance models.
    -> SimM (Map Int Int)                -- ^ Count of racer-rank pairs.
runExperimentPip nRuns selPip pips =
    fst <$> P.fold' updateCount counters id entries
    where
    counters = Map.empty
    entries = simulatorN nRuns pips >-> P.mapFoldable id
        >-> P.filter ((selPip ==) . simPip)
    updateCount ctrs se =
        Map.alter (maybe (Just 1) (Just . (+ 1))) (pipRank se) ctrs
    -- Note that this doesn't seem to be noticeably faster than
    -- runExperimentFull.

-- | Generates results of many races with an added probe and tallies how often
-- the probe reached each position.
runExperimentProbe
    :: EloOptions                        -- ^ Number of runs, seed, and probe.
    -> [(SimPip, OrbitalDistribution)]   -- ^ Identifications and performance models.
    -> SimM (Map Int Int)                -- ^ Count of racer-rank pairs.
runExperimentProbe eopts pips =
    runExperimentPip nRuns probeId ((probeId, probeModel) : pips)
    where
    nRuns = simRuns eopts
    probeModel = orbitalDistr (eloGammaShape eopts) (simProbeRating eopts)
    probeId = Probe 0

simModelStrength :: EloOptions -> Ratings -> SimM Double
simModelStrength eopts = fmap ((fromIntegral nRuns /) . fromIntegral
        . Map.foldl' (+) 0 . Map.filterWithKey (\key _ -> isTargetReached key))
    . runExperimentProbe eopts . toSimPips (eloGammaShape eopts) []
    where
    nRuns = simRuns eopts
    isTargetReached = (<= simTarget eopts)

-- | Average positions attained by racers in the simulations.
simAveragePositions :: EloOptions -> Ratings -> SimM (Map SimPip Double)
simAveragePositions eopts = fmap (fmap (/ fromIntegral nRuns)
        . Map.mapKeysWith (+) fst
        . Map.mapWithKey (\(p, i) m -> fromIntegral i * fromIntegral m))
    . runExperimentFull nRuns . toSimPips (eloGammaShape eopts) []
    where
    nRuns = simRuns eopts

example = bimap SimPip (orbitalDistr 3)
    <$> [("HAM", 500), ("BOT", 286), ("VER", 430), ("VET", 380), ("STR", 240)]

(testSeed1, testSeed2) =
    ( (toSeed . V.fromList)
    [2983729257,1613295996,1536938243,3382368331,999233052,3403717367,1505805409
    ,457896502,2924789643,126275855,1602675974,2774079579,2608979963
    ,1285505154,2190748158,2008670062,2323612213,1474206092,1096337431
    ,288815539,2167832047,89478764,3770564965,1831684145,4254467873
    ,4134720386,890644500,541497865,2034198580,1508666107,283455713
    ,3816390539,1145530897,2308506764,842610843,3446610465,2194641952
    ,1640187121,2991202145,2011144362,3678556400,2037948716,1876884772
    ,2860215653,1643680172,2482193874,3509845458,3561221336,4142605894
    ,2783016144,3695356134,2149339665,3027240920,3779754543,1091699906
    ,1296315652,4201840705,466463827,2771261985,3609434098,3174518634
    ,1883778791,1338358979,2103001057,1010315523,1850466367,1704534333
    ,1127116119,494022984,4131810412,1989435256,2189717985,1195790672
    ,546149896,3956437595,3430065006,996894469,652151680,2444801896
    ,3277167510,3898937104,1816281959,3887115177,756915671,1965221836
    ,3217297445,2524193384,1617990254,17611352,3883903320,654722342
    ,3035693226,2959956776,2480087566,1416574980,2332603833,1114892984
    ,3972168584,3190304332,2279966641,3984143279,506682112,708954088
    ,782831114,3580552292,1328043039,1445323362,545935803,3511911499
    ,3532637593,3092002176,1293437661,866180493,470170908,2213099343
    ,2610774100,770498881,1522960121,362223202,368376677,2696502499
    ,3735383147,959672244,3597809455,1480371635,1938962987,3574654602
    ,3123103398,3899686612,111986794,555794392,433322464,3512122633
    ,888439676,4013356910,1537488858,2998004672,1503517187,1676041891
    ,4284517682,4187913039,526411011,3844452821,2843102260,1863455660
    ,2745546381,3458778584,986126662,623225133,484127313,564056677
    ,3412008153,2637822297,2049738317,801004491,3845660219,598766095
    ,280702759,3220956212,2559119027,1387002130,1860669793,3878600592
    ,1646571175,2650709733,287304878,2708941441,3725941329,1923856846
    ,1348165256,360727946,734818281,3696343351,2757132999,3994319050
    ,821963309,4053106446,3956385121,3738530775,934656047,3471283180
    ,2977709497,1234997742,1733246858,1772131303,3541096545,1687511831
    ,628693315,3546632502,2171179846,2127187761,2008195960,537150990
    ,2347439733,368798333,1489117236,1030208402,1793664579,3357536366
    ,1379818798,1445263925,4136224347,2107945916,3629803435,2568860630
    ,4168840089,2151673581,2471786285,2480203426,3333035688,705116558
    ,3745612939,1227314979,1292798725,612448134,3396651774,975841030
    ,624512980,4190767123,1233459964,3008972087,601532817,3383745831
    ,773671035,1983882269,146329547,2740348494,1651745877,728960915
    ,3716272563,1763505775,3322684912,3890606420,213721570,223357770
    ,2416381159,1421939056,189957726,776510638,3327921427,3756638545
    ,3518106256,1713504596,1695176588,2195699141,1866690002,2307523501
    ,3633434759,832509879,2504451685,2104731942,2078007247,608050460
    ,1686619141,3763034973,1595031540,255,362436]
    , (toSeed . V.fromList)
    [3514361481,3842361471,795014533,1784838522,526088448,646126946,2812735612
    ,922018427,3378711992,782018649,1464785869,2983406718,2419633236
    ,1853276290,309755073,1617742846,4109644868,2445006221,1147110543
    ,2791681101,708108117,1877561500,86076843,3979878655,2149817768
    ,927290670,2660510741,976560117,3170810766,926610827,3509847976
    ,1108252766,601807489,4115751470,3120488071,3072205774,3290517234
    ,1820511216,3050864435,571828010,272780034,1808615920,2680969109
    ,196720231,2217492074,3498903999,978244655,2904543907,462407936
    ,2767797777,1805729825,471152103,2285843480,2591886531,1835500868
    ,2217326776,3533974798,1542713276,3796979153,924020557,1821142792
    ,2849580282,3686348589,1449495310,2253563057,2165272312,2064015893
    ,3905773535,2412868076,450944229,3801900089,2996821856,3902541163
    ,1622237188,1727880413,938683952,3672638018,1520966709,856986815
    ,1262553468,4209303792,1423083963,1481088361,4220947014,980218971
    ,2211614421,2023491475,873192432,258917836,2590338598,2264586590
    ,2935412575,3510677396,2083579827,1902702209,1248330234,763375492
    ,255297030,2164301963,2217557972,3368312704,1113206599,3368064395
    ,3721694346,4269798761,4101962784,1234564983,1936596601,1270230244
    ,1248522168,3808964907,3440012765,2846369544,1495023958,2839177494
    ,2358176524,1639191498,106020455,2934873951,3137110660,426604379
    ,4081972914,1956316804,427508886,710455103,3967527362,3510393578
    ,791432587,2989263762,4159249582,3959193725,82625198,1773816659
    ,1187237560,177169928,4088291379,1342962495,375846631,34078710
    ,4259092681,2238989091,2401696247,3212319020,2081634154,644007076
    ,1450607183,25503484,3904269573,839898787,1938851346,296165563
    ,472281688,3885688307,1744171857,107901069,3687678329,3225343202
    ,4274663955,2691973102,2944086281,560196343,417575873,746385354
    ,2817008926,2007973603,2800956071,4171484727,3441556823,788429435
    ,773478323,2946141746,4117575557,1835102299,1247443170,3285737716
    ,80965588,528551827,928476088,652810025,1728639523,3172597791
    ,3195498619,1668001640,1090898405,3635992671,724319477,3085777054
    ,3197077886,1892186618,3464705939,3263299682,2045281371,4056683611
    ,3422245534,281426931,2906982745,2362266837,2311470939,4024086179
    ,296904488,3549811521,3590869187,30877516,2190023718,1129803354
    ,2329929231,3208052139,3152957852,1071654726,2578550933,2962388832
    ,3992806217,3387397289,3784382569,3386603360,1782091423,133382899
    ,3644533203,3133883996,3443143560,2901080115,218653140,2133128666
    ,2660840034,3437312141,992717558,1663297758,2357872362,2910300196
    ,899348458,3304500258,612813003,3888003996,736806268,3876075907
    ,2968262906,2686602214,1787165222,1401421921,4129502730,464353724
    ,2707726113,4072910765,3141835839,2016856979,1313588609,3290206241
    ,3656476319,2831158121,4267011120,925306905,282140482,3869026817
    ,1821837849,538991518,915864998,221,539427143])

-- >$> :set -XOverloadedStrings
-- $>
-- >$> :set +s
-- $>
-- >$> evalSimM Nothing $ simulateSingleRace Nothing example
--
-- >$> evalSimM Nothing $ runExperimentFull 10000 example

max_stars_repo_name: duplode/elo-zs

max_issues_repo_name: duplode/elo-zs

max_forks_repo_name: duplode/elo-zs


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

import Graphics.UI.GLUT
import Control.Monad
import Data.Int
import Data.Complex
import System.Exit ( exitWith, ExitCode(ExitSuccess) )
import qualified Data.Map as Map

start = 5
iterations = 7
maxIterations = 30

x // y = fromIntegral x / fromIntegral y

-- Divides [a] into [[a], [a], ...] with each sublist of length n,
-- except the last sublist which has length <= n.
chunkify n [] = []
chunkify n xs = let (xs', rest) = splitAt n xs
                in xs' : chunkify n rest

-- Converts a coordinate in screen space to a vertex.
pix2vert (Size w h) (x, y) =
  Vertex2 ((3 // w * fromIntegral x) - 2.0) ((2 // h * fromIntegral y) - 1.0)

-- List of all of the vertices that represent screen pixels.
vertices :: IO [Vertex2 GLfloat]
vertices = get windowSize >>= (\x -> return $ vertices' x)

vertices' (Size w h) = [pix2vert (Size w h) (x, y) | x <- [0..w-1], y <- [0..h-1]]

-- Gets the color for a number of iterations.
color3 r g b = Color3 r g b

getcolor :: Int -> Color3 Float
getcolor numberOfIterations
  | inMandelbrotSet = color3 0 0 0
  | otherwise = color3 (escapeNumber*0.5) (escapeNumber * 0.5) escapeNumber
  where escapeNumber = numberOfIterations // iterations
        inMandelbrotSet = numberOfIterations == iterations

-- Returns the number of iterations <= the maximum iterations of the
-- Mandelbrot set at the given vertex.
mandel numberOfIterations (Vertex2 r i) =
  length . takeWhile (\z -> magnitude z <= 2) .
  take numberOfIterations $ iterate (\z -> z^2 + (r :+ i)) 0

mandel' numberOfIterations (Vertex2 r i) =
  length . takeWhile (\z -> magnitude z <= 2) .
  take numberOfIterations $ iterate (\z -> z^2 + (r :+ i)) 0

mandelMap :: (RealFloat a1, Fractional a2) => Vertex2 a1 -> Map.Map Int (Color3 a2)
mandelMap (Vertex2 r i) = Map.fromList iterList
  where mandelNumbers = iterate (\z -> z^2 + (r :+ i)) 0

        escapeVelocity numIter = length . takeWhile (\z -> magnitude z <= 2) .
                                 take numIter $ mandelNumbers

        --iterList = makeMap (\iter -> getColor (escapeVelocity iter) iter) [1..maxIterations]
        iterList = makeMap (\iter -> getColor (escapeVelocity iter) iter) [start..iterations]

        getColor escapedAt totalNumIters
          | inMandelbrotSet = color3 0 0 0
          | otherwise = color3 (escapeNumber*0.5) (escapeNumber * 0.5) escapeNumber
          where escapeNumber = escapedAt // totalNumIters
                inMandelbrotSet = escapedAt == totalNumIters

colorPoint :: RealFloat a => Int -> Vertex2 a -> Color3 Float
colorPoint numIterations = getcolor . (mandel numIterations)

-- plots one point.
-- drawVert v = do
--   color . getcolor $ mandel v
--   vertex v

drawVert' (point, color_) = do
  color color_
  vertex point

-- draws all the vertices in slices (to update the display while drawing).
-- display' chunks = do
  -- mapM_ render chunks
  -- displayCallback $= display
  -- where render points = do
          -- renderPrimitive Points (mapM_ drawVert points)
          -- flush

getColor :: (Map.Map (Vertex2 GLfloat) (Map.Map Int (Color3 Double))) -> Int -> Vertex2 GLfloat -> Color3 Double
getColor pointMap numIter vertex =
  case Map.lookup vertex pointMap of
    Just colorMap -> case Map.lookup numIter colorMap of
                       Just r -> r
                       Nothing -> Color3 0 0 0
    Nothing -> Color3 0 0 0

drawAll :: Int ->  Map.Map (Vertex2 GLfloat) (Map.Map Int (Color3 Double)) -> [Vertex2 GLfloat] -> IO()
drawAll count pointMap points = do
  render
  displayCallback $= (display (count + 1) pointMap)
  where
        render = do
          renderPrimitive Points (mapM_ drawVert' toDraw)
          flush

        toDraw :: [(Vertex2 GLfloat, Color3 Double)]
        toDraw = makeMap (getColor pointMap count) points

-- draws the whole fractal
display :: Int -> Map.Map (Vertex2 GLfloat) (Map.Map Int (Color3 Double)) -> IO ()
display 200 _ = do exitWith ExitSuccess
display num pointMap = do
  clear [ ColorBuffer ]
  -- displayCallback $= (vertices >>= display' . chunkify 256)
  displayCallback $= (vertices >>= (drawAll num pointMap))
  get currentWindow >>= postRedisplay

keyboard :: KeyboardMouseCallback
keyboard (Char keyPressed) Down _ _ = case keyPressed of
   ' '   -> postRedisplay Nothing
   '\27' -> exitWith ExitSuccess
   _     -> return ()
keyboard _ _ _ _ = return ()

mandelMain = do
  getArgsAndInitialize
  initialDisplayMode $= [ SingleBuffered, RGBMode]
  initialWindowSize $= Size 1200 1024
  initialWindowPosition $= Position 100 100
  createWindow "Mandelbrot"
  --fullScreen
  -- renderString Roman "Test"
  clearColor $= Color4 0 0 0 0
  matrixMode $= Projection
  loadIdentity
  ortho (-2) 1 (-1) 1 (-1) 1
  keyboardMouseCallback $= Just keyboard
  points <- vertices
  putStrLn "PreLoading.."
  let !preload = Map.fromList $ makeMap mandelMap points
  putStrLn "Loaded"
  displayCallback $= (display start preload)
  mainLoop

makeMap :: (a -> b) -> [a] -> [(a, b)]
makeMap fn lst = zip lst (map fn lst)

max_stars_repo_name: jordangedney/mandelbrot-explorer

max_issues_repo_name: jordangedney/mandelbrot-explorer

max_forks_repo_name: jordangedney/mandelbrot-explorer


text:
{-# LANGUAGE ViewPatterns        #-}
--------------------------------------------------------------------------------
-- |
-- Module      : ArrayFire.BLAS
-- Copyright   : David Johnson (c) 2019-2020
-- License     : BSD3
-- Maintainer  : David Johnson <djohnson.m@gmail.com>
-- Stability   : Experimental
-- Portability : GHC
--
-- Basic Linear Algebra Subprograms (BLAS) API
--
-- @
-- main :: IO ()
-- main = print (matmul x y xProp yProp)
--  where
--     x,y :: Array Double
--     x = matrix (2,3) [[1,2],[3,4],[5,6]]
--     y = matrix (3,2) [[1,2,3],[4,5,6]]
--
--     xProp, yProp :: MatProp
--     xProp = None
--     yProp = None
-- @
-- @
--  ArrayFire Array
--  [2 2 1 1]
--     22.0000    49.0000
--     28.0000    64.0000
-- @
--------------------------------------------------------------------------------
module ArrayFire.BLAS where

import Data.Complex

import ArrayFire.FFI
import ArrayFire.Internal.BLAS
import ArrayFire.Internal.Types

-- | The following applies for Sparse-Dense matrix multiplication.
--
-- This function can be used with one sparse input. The sparse input must always be the lhs and the dense matrix must be rhs.
--
-- The sparse array can only be of 'CSR' format.
--
-- The returned array is always dense.
--
-- optLhs an only be one of AF_MAT_NONE, AF_MAT_TRANS, AF_MAT_CTRANS.
--
-- optRhs can only be AF_MAT_NONE.
--
-- >>> matmul (matrix @Double (2,2) [[1,2],[3,4]]) (matrix @Double (2,2) [[1,2],[3,4]]) None None
-- ArrayFire Array
-- [2 2 1 1]
--     7.0000    15.0000
--    10.0000    22.0000
matmul
  :: AFType a
  => Array a
  -- ^ 2D matrix of Array a, left-hand side
  -> Array a
  -- ^ 2D matrix of Array a, right-hand side
  -> MatProp
  -- ^ Left hand side matrix options
  -> MatProp
  -- ^ Right hand side matrix options
  -> Array a
  -- ^ Output of 'matmul'
matmul arr1 arr2 prop1 prop2 = do
  op2 arr1 arr2 (\p a b -> af_matmul p a b (toMatProp prop1) (toMatProp prop2))

-- | Scalar dot product between two vectors. Also referred to as the inner product.
--
-- >>> dot (vector @Double 10 [1..]) (vector @Double 10 [1..]) None None
-- ArrayFire Array
-- [1 1 1 1]
--   385.0000
dot
  :: AFType a
  => Array a
  -- ^ Left-hand side input
  -> Array a
  -- ^ Right-hand side input
  -> MatProp
  -- ^ Options for left-hand side. Currently only AF_MAT_NONE and AF_MAT_CONJ are supported.
  -> MatProp
  -- ^ Options for right-hand side. Currently only AF_MAT_NONE and AF_MAT_CONJ are supported.
  -> Array a
  -- ^ Output of 'dot'
dot arr1 arr2 prop1 prop2 =
  op2 arr1 arr2 (\p a b -> af_dot p a b (toMatProp prop1) (toMatProp prop2))

-- | Scalar dot product between two vectors. Also referred to as the inner product. Returns the result as a host scalar.
--
-- >>> dotAll (vector @Double 10 [1..]) (vector @Double 10 [1..]) None None
-- 385.0 :+ 0.0
dotAll
  :: AFType a
  => Array a
  -- ^ Left-hand side array
  -> Array a
  -- ^ Right-hand side array
  -> MatProp
  -- ^ Options for left-hand side. Currently only AF_MAT_NONE and AF_MAT_CONJ are supported.
  -> MatProp
  -- ^ Options for right-hand side. Currently only AF_MAT_NONE and AF_MAT_CONJ are supported.
  -> Complex Double
  -- ^ Real and imaginary component result
dotAll arr1 arr2 prop1 prop2 = do
  let (real,imag) =
        infoFromArray22 arr1 arr2 $ \a b c d ->
          af_dot_all a b c d (toMatProp prop1) (toMatProp prop2)
  real :+ imag

-- | Transposes a matrix.
--
-- >>> array = matrix @Double (2,3) [[2,3],[3,4],[5,6]]
-- >>> array
-- ArrayFire Array
-- [2 3 1 1]
--     2.0000     3.0000     5.0000
--     3.0000     4.0000     6.0000
--
-- >>> transpose array True
-- ArrayFire Array
-- [3 2 1 1]
--     2.0000     3.0000
--     3.0000     4.0000
--     5.0000     6.0000
--
transpose
  :: AFType a
  => Array a
  -- ^ Input matrix to be transposed
  -> Bool
  -- ^ Should perform conjugate transposition
  -> Array a
  -- ^ The transposed matrix
transpose arr1 (fromIntegral . fromEnum -> b) =
  arr1 `op1` (\x y -> af_transpose x y b)

-- | Transposes a matrix.
--
-- * Warning: This function mutates an array in-place, all subsequent references will be changed. Use carefully.
--
-- >>> array = matrix @Double (2,2) [[1,2],[3,4]]
-- >>> array
-- ArrayFire Array
-- [3 2 1 1]
--    1.0000     4.0000
--    2.0000     5.0000
--    3.0000     6.0000
--
-- >>> transposeInPlace array False
-- >>> array
-- ArrayFire Array
-- [2 2 1 1]
--    1.0000     2.0000
--    3.0000     4.0000
--
transposeInPlace
  :: AFType a
  => Array a
  -- ^ Input matrix to be transposed
  -> Bool
  -- ^ Should perform conjugate transposition
  -> IO ()
transposeInPlace arr (fromIntegral . fromEnum -> b) =
  arr `inPlace` (`af_transpose_inplace` b)

max_stars_repo_name: arrayfire/arrayfire-haskell

max_issues_repo_name: arrayfire/arrayfire-haskell

max_forks_repo_name: arrayfire/arrayfire-haskell


text:
-- Haskell code for navigating the interface between lists of edges and
-- adjacency matrices.

module Adjacency where

import Numeric.LinearAlgebra hiding (fromList)
import qualified Numeric.LinearAlgebra as LA (fromList)
import Numeric.LinearAlgebra.Data
import qualified Data.Map.Strict as M
import qualified Data.List as L
import Data.Ord
import Data.KMeans

-- Data structure for graphs and networks. Just a list of edges.
data Graph a = Graph { nodes :: [a]
                     , edges :: [(a, a)]
                     } deriving (Show, Eq)

data Graph' a = Graph' { nodes' :: [a]
                       , edges' :: M.Map a [a]
                       } deriving (Show, Eq)

data Graph'' a = Graph'' { nodes'' :: [a]
                         , edges'' :: Matrix Double
                         } deriving (Show, Eq)

-- Converts the edge list to a map structure. Takes up more space, but it's
-- faster for looking up what else is connected to a particular node.
toMap :: (Ord a, Eq a) => Graph a -> Graph' a
toMap (Graph ns es) = Graph' ns es'
  where
    fInsert o o' m x = M.insertWith (++) (o x) [o' x] m
    iMap = M.fromList $ map (\x -> (x, [])) ns
    es' = L.foldl' (fInsert snd fst) (L.foldl' (fInsert fst snd) iMap es) es

-- Converts a graph (edge list) to an adjacency matrix structure.
toAdjacency :: (Ord a, Eq a) => Graph a -> Graph'' a
toAdjacency gr@(Graph ns es) = Graph'' ns m
  where
    n = length ns
    m = (n >< n) . concat . map (to01List ns) $ map ((M.!) . edges' $ toMap gr) ns

-- Converts a graph (edge list) to an adjacency matrix structure. About 1/3
-- slower than `toAdjacency`.
-- toAdjacency' :: (Ord a, Eq a) => Graph a -> Graph'' a
-- toAdjacency' gr@(Graph ns es) = Graph'' ns m
--   where
--     n = length ns
--     m = (n >< n) . M.foldl' (++) [] $ M.map (to01List ns) (edges' $ toMap es)

-- Given a list of edges which are all connected to a node, generates
-- a list of 0's and 1's denoting membership.
to01List ::  (Eq a) => [a] -> [a] -> [Double]
to01List ns es = map (boolTo01 . flip elem es) ns

-- Converts False/True to 0/1.
boolTo01 :: Bool -> Double
boolTo01 b
  | b    = 1
  | True = 0

-- SVD approximation to a graph (uses k-element SVD of the adjacency
-- matrix). In reality we need the SVD of XY^T, but that will come later.
svdGraph :: Int -> Graph'' a -> Graph'' a
svdGraph k (Graph'' ns es) = Graph'' ns $ compress k es

-- This function was written by Nicolas Favre-Felix. It returns an
-- approximation of the given matrix, by taking only the k largest entries
-- of Sigma.
compress :: Int -> Matrix Double -> Matrix Double
compress k m = u_k <> sigma_k <> v_k where
  (u,sigma,v) = svd m  -- get SVD
  sigma_k = takeBoth k $ diag sigma -- keep k values of Σ
  u_k = takeColumns k u    -- keep k columns of U
  v_k = takeRows k $ tr v   -- keep k rows of v

-- Constructs the graph Laplacian.
laplacian :: (Ord a, Eq a) => Graph a -> Matrix Double
laplacian gr = (diag $ LA.fromList counts) - (edges'' $ toAdjacency gr)
  where
    counts = M.foldl (\x y -> x ++ [fromIntegral $ length y]) [] . edges'
             $ toMap gr

-- A little wrapper to simplify taking the upper submatrix.
takeBoth :: (Element t) => Int -> Matrix t -> Matrix t
takeBoth d = takeColumns d . takeRows d

-- Spanning tree counter; uses Kirchhoff's Theorem.
spanTreeCount :: (Ord a, Eq a) => Graph a -> Integer
spanTreeCount gr@(Graph ns es) = round . det . takeBoth n' $ laplacian gr
  where
    n' = length ns - 1

-- spectralCluster performs spectral clustering to estimate group membership.
spectralCluster :: (Ord a, Eq a) => Int -> Graph a -> [[[Double]]]
spectralCluster k gr@(Graph ns es) = kmeans k u
  where
    n = length ns
    -- u = map toList . drop (n-k) . toColumns . snd . eigSH . sym $ laplacian gr
    u = map toList . drop (n-k) . toColumns . snd . eigSH $ laplacian gr

-- spectralCluster performs spectral clustering to estimate group
-- membership. Estimates the best choice of k using consecutive
-- differences of the eigenvalues.
spectralCluster' :: (Ord a, Eq a) => Graph a -> [[[Double]]]
spectralCluster' gr@(Graph ns es) = kmeans k u
  where
    n = length ns
    -- e = eigSH . sym $ laplacian gr
    e = eigSH $ laplacian gr
    k = (+1) . fst . L.minimumBy (comparing snd) . zip [0..] . reverse
        . diff . toList $ fst e
    u = map toList . drop (n-k) . toColumns $ snd e

lambdaSelect :: (Ord a, Eq a) => Graph a -> Int
lambdaSelect gr@(Graph ns es) = k
  where
    n = length ns
    -- e = eigSH . sym $ laplacian gr
    e = eigSH $ laplacian gr
    k = (+1) . fst . L.minimumBy (comparing snd) . zip [0..] . reverse
        . diff . toList $ fst e
    u = map toList . drop (n-k) . toColumns $ snd e

-- Consecutive differences.
diff :: (Num a) => [a] -> [a]
diff xs = zipWith (-) (drop 1 xs) xs

max_stars_repo_name: kimolas/causal-haskell

max_issues_repo_name: kimolas/causal-haskell

max_forks_repo_name: kimolas/causal-haskell


text:
{-|
Module      : Julia
Description : Computing Julia sets
Copyright   : (c) Jason Chen, 2021
License     : MIT
Maintainer  : chezbgone (at) gmail dot com
Stability   : experimental
-}
{-# LANGUAGE NamedFieldPuns #-}
module Julia where

import Data.Complex
import Data.List (iterate')
import Codec.Picture

import Interpolation

data Resolution = Resolution { width :: Int, height :: Int }

data Bounds = Bounds { top    :: Double
                     , bottom :: Double
                     , left   :: Double
                     , right  :: Double
                     }

data ComputationProps = ComputationProps
    { maxIterations :: Int
    , escapeRadius :: Double
    }

data RenderProps = RenderProps
    { resolution :: Resolution
    , bounds     :: Bounds
    }

data JuliaProps = JuliaProps
    { computationProps :: ComputationProps
    , renderProps      :: RenderProps
    }

escapeTime :: (Complex Double -> Complex Double) -- f
           -> ComputationProps
           -> Complex Double -- init
           -> Int
escapeTime f ComputationProps{ maxIterations, escapeRadius } init =
    length $ takeWhile (\z -> magnitude z < escapeRadius)
           $ take maxIterations
           $ iterate' f init

escapePixel :: Int -> PixelRGB8
escapePixel = undefined

gridToComplex :: RenderProps -> Int -> Int -> Complex Double
gridToComplex RenderProps{ resolution, bounds } x y
  = (left + xStep * fromIntegral x) :+ (top + yStep * fromIntegral y)
      where Resolution {width, height} = resolution
            Bounds {top, bottom, left, right} = bounds
            xStep = (right - left) / fromIntegral width
            yStep = (top - bottom) / fromIntegral height

generateJulia :: JuliaProps
              -> (Complex Double -> Complex Double) -- f
              -> FilePath
              -> IO ()
generateJulia JuliaProps{ computationProps, renderProps }
              f filepath = writePng filepath image
    where image = generateImage compute width height
          Resolution{ width, height } = resolution renderProps
          compute :: Int -> Int -> PixelRGB8
          compute x y = escapePixel $ escapeTime f computationProps
                                    $ gridToComplex renderProps x y

max_stars_repo_name: chezbgone/chaos-renderer

max_issues_repo_name: chezbgone/chaos-renderer

max_forks_repo_name: chezbgone/chaos-renderer


text:
{-# LANGUAGE NamedFieldPuns #-}

module School.Unit.Affine
( affine ) where

import Numeric.LinearAlgebra ((<>), R, add, tr')
import School.Unit.Unit (Unit(..))
import School.Unit.UnitActivation (UnitActivation(..))
import School.Unit.UnitParams (UnitParams(..))
import School.Unit.UnitGradient (UnitGradient(..))
import School.Utils.LinearAlgebra (mapRows, sumRows)

affine :: Unit Double
affine = Unit
  { apply = affineApply
  , deriv = affineDeriv
  }

affineApply :: UnitParams R
            -> UnitActivation R
            -> UnitActivation R
affineApply AffineParams { affineBias, affineWeights }
            (BatchActivation inMatrix) =
  BatchActivation $ mapRows (add affineBias) (inMatrix <> tr' affineWeights)
affineApply _ (ApplyFail msg) = ApplyFail msg
affineApply _ _ = ApplyFail "Failure when applying affine unit"

failure :: String
        -> ( UnitGradient R
           , UnitParams R
           )
failure msg = (GradientFail msg, EmptyParams)

affineDeriv :: UnitParams R
            -> UnitGradient R
            -> UnitActivation R
            -> ( UnitGradient R
               , UnitParams R
               )
affineDeriv AffineParams { affineWeights }
            (BatchGradient inGrad)
            (BatchActivation input) =
  (outGrad, paramDerivs) where
    outGrad = BatchGradient $ inGrad <> affineWeights
    bias = sumRows inGrad
    weights = tr' inGrad <> input
    paramDerivs = AffineParams { affineBias = bias
                               , affineWeights = weights
                               }
affineDeriv _ (GradientFail msg) _ = failure msg
affineDeriv _ _ (ApplyFail msg) = failure msg
affineDeriv _ _ _ = failure "Failure when differentiating affine unit"

max_stars_repo_name: jfulseca/School

max_issues_repo_name: jfulseca/School

max_forks_repo_name: jfulseca/School


text:
module Statistics.MLE (
	maxLikelihood,
	probit
	) where

import Numeric.GSL.Minimization
import Statistics.Distribution
import Statistics.Distribution.Normal
import Statistics.Matrix
import qualified Data.Vector.Generic as G
import qualified Data.Vector.Unboxed as U

data MaxControl = MaxControl { maxit :: Int, epsilon :: Double } deriving (Show)

type ObjectiveF = [Double] -> Double
type GradientF = [Double] -> [Double]
type LikelihoodModel = Matrix -> Vector -> (ObjectiveF, GradientF)

defaultMaxCtrl :: MaxControl
defaultMaxCtrl = MaxControl 100 10e-2

maxLikelihood :: LikelihoodModel -> MaxControl -> Matrix -> Vector -> Vector
maxLikelihood model ctrl x y = U.fromList . fst $ bhat
	where
		bhat = uncurry (minimizeD VectorBFGS2 (epsilon ctrl) (maxit ctrl) 1.0 0.1) (model x y) sval
		sval = replicate (snd (dimension x)) (epsilon ctrl)

llprobit :: LikelihoodModel
llprobit x y = (loglik, grad)
	where
		loglik b = negate $ G.sum $ G.map contr (G.zip y (x `multiplyV` U.fromList b))
		grad b = fmap negate $ U.toList $ transpose x `multiplyV` G.map gradcontr (G.zip y (x `multiplyV` U.fromList b))
		contr (y,xb) = y * log (cdf xb) + (1.0 - y) * log (1.0 - cdf xb)
		gradcontr (y,xb) = y * pdf xb / cdf xb - (1 - y) * pdf xb / (1 - cdf xb)
		pdf = density (normalDistr 0 1)
		cdf x = min (max (cumulative (normalDistr 0 1) x) (epsilon defaultMaxCtrl)) (1 - epsilon defaultMaxCtrl) -- BUG in cumulative: range not [0,1]

probit :: Matrix -> Vector -> Vector
probit = maxLikelihood llprobit defaultMaxCtrl

max_stars_repo_name: StefanHubner/statistics

max_issues_repo_name: StefanHubner/statistics

max_forks_repo_name: StefanHubner/statistics


text:
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE TypeSynonymInstances #-}
module Statistics.BBVI.Distribution.Normal
  ( NormalDist(..)
  , normalDistr
  , diffableNormalLogProb
  )
where

import qualified Data.Vector                   as V
import           Numeric.MathFunctions.Constants
                                                ( m_sqrt_2_pi )
import           Statistics.BBVI.Class
import           Statistics.BBVI.Propagator     ( SampleDouble )
import           System.Random.MWC.Distributions
                                                ( normal
                                                , standard
                                                )
-- | datatype for normal distribution
data NormalDist = ND {mean :: Double, stdDev :: Double} deriving (Show, Eq, Ord, Read)

-- | build a normal distribution
normalDistr
  :: Double -- ^ mean
  -> Double -- ^ standard deviation
  -> NormalDist
normalDistr mu std = ND mu (max 1e-10 std)

instance DistUtil NormalDist where
  fromParamVector xs = normalDistr (xs V.! 0) (xs V.! 1)
  toParamVector d = V.fromList [mean d, stdDev d]
  nParams _d = 2

instance Dist NormalDist SampleDouble where
  resample d gen = normal (mean d) (stdDev d) gen
  logProb d x = (-xm * xm / (2 * sd * sd)) - ndPdfDenom
   where
    xm         = x - mean d
    sd         = stdDev d
    ndPdfDenom = log $ m_sqrt_2_pi * sd
  -- |
  -- >>> paramGradOfLogQ (normalDistr 0.0 3.0) (2.0 :: Double)
  -- [0.2222222222222222,-0.18518518518518517]
  paramGradOfLogQ d x = V.fromList
    [ xm / std ^ (2 :: Int)
    , 1 / std ^ (3 :: Int) * (xm ^ (2 :: Int) - std ^ (2 :: Int))
    ]
   where
    std = stdDev d
    xm  = x - mean d

-- | log density of a normal distribution, with a type that allows us
-- to use automatic differentiation from the "ad" library.
diffableNormalLogProb
  :: Floating a
  => a -- ^ mean
  -> a -- ^ standard deviation
  -> a -- ^ observation
  -> a -- ^ log density
{-# INLINE diffableNormalLogProb #-}
diffableNormalLogProb mu sd x = -xm * xm / (2 * sd * sd) - ndPdfDenom
 where
  xm         = x - mu
  ndPdfDenom = log $ sqrt (2 * pi) * sd

instance Differentiable NormalDist SampleDouble where
  -- |
  -- >>> (transform (normalDistr 0.0 2.0) 1.0 :: Double)
  -- 2.0
  transform d samp = mean d + stdDev d * samp
  epsilon _d gen = standard gen
  sampleGradOfLogQ d z = -(z - mean d) / (stdDev d ** 2)
  -- |
  -- >>> gradTransform (normalDistr 0.0 1.0) (2.0 :: Double)
  -- [1.0,2.0]
  gradTransform _d samp = V.fromList [1.0, samp] -- ND $ V.fromList [1.0 , eps * stdDev d] -- --

max_stars_repo_name: massma/propagator-bbvi

max_issues_repo_name: massma/propagator-bbvi

max_forks_repo_name: massma/propagator-bbvi


text:
{-# LANGUAGE ForeignFunctionInterface #-}
module Grenade.Layers.Internal.Pad (
    pad
  , crop
  ) where

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

import           Foreign ( mallocForeignPtrArray, withForeignPtr )
import           Foreign.Ptr ( Ptr )

import           Numeric.LinearAlgebra ( flatten, Matrix )
import qualified Numeric.LinearAlgebra.Devel as U

import           System.IO.Unsafe ( unsafePerformIO )

pad :: Int -> Int -> Int -> Int -> Int -> Int -> Int -> Int -> Int -> Matrix Double -> Matrix Double
pad channels padLeft padTop padRight padBottom rows cols rows' cols' m
 = let outMatSize      = rows' * cols' * channels
       vec             = flatten m
   in unsafePerformIO $ do
     outPtr        <- mallocForeignPtrArray outMatSize
     let (inPtr, _) = U.unsafeToForeignPtr0 vec

     withForeignPtr inPtr $ \inPtr' ->
       withForeignPtr outPtr $ \outPtr' ->
         pad_cpu inPtr' channels rows cols padLeft padTop padRight padBottom outPtr'

     let matVec = U.unsafeFromForeignPtr0 outPtr outMatSize
     return (U.matrixFromVector U.RowMajor (rows' * channels) cols' matVec)
{-# INLINE pad #-}

foreign import ccall unsafe
    pad_cpu
      :: Ptr Double -> Int -> Int -> Int -> Int -> Int -> Int -> Int -> Ptr Double -> IO ()

crop :: Int -> Int -> Int -> Int -> Int -> Int -> Int -> Int -> Int -> Matrix Double -> Matrix Double
crop channels padLeft padTop padRight padBottom rows cols _ _ m
 = let outMatSize      = rows * cols * channels
       vec             = flatten m
   in unsafePerformIO $ do
     outPtr        <- mallocForeignPtrArray outMatSize
     let (inPtr, _) = U.unsafeToForeignPtr0 vec

     withForeignPtr inPtr $ \inPtr' ->
       withForeignPtr outPtr $ \outPtr' ->
         crop_cpu inPtr' channels rows cols padLeft padTop padRight padBottom outPtr'

     let matVec = U.unsafeFromForeignPtr0 outPtr outMatSize
     return (U.matrixFromVector U.RowMajor (rows * channels) cols matVec)

foreign import ccall unsafe
    crop_cpu
      :: Ptr Double -> Int -> Int -> Int -> Int -> Int -> Int -> Int -> Ptr Double -> IO ()

max_issues_repo_name: Alien-Inc/grenade

max_forks_repo_name: Alien-Inc/grenade


text:
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE TypeSynonymInstances #-}
{-# OPTIONS_GHC -fno-warn-orphans #-}

module IHaskell.Display.Spiral where

import Data.Complex.Cyclotomic (Cyclotomic)
import qualified Data.Text as T
import IHaskell.Display
import Text.LaTeX

import Spiral.Array
import Spiral.Exp
import Spiral.SPL
import Spiral.Util.Pretty.LaTeX

displayLaTeX :: Pretty a => a -> IO Display
displayLaTeX = display . IHaskell.Display.latex . T.unpack . render . math . ppr

instance IHaskellDisplay Cyclotomic where
  display = displayLaTeX

instance IHaskellDisplay (Const a) where
  display = displayLaTeX

instance IHaskellDisplay (Exp a) where
  display = displayLaTeX

instance (Pretty a, IArray r DIM2 a) => IHaskellDisplay (Matrix r a) where
  display = displayLaTeX

instance (Num a, Pretty a) => IHaskellDisplay (SPL a) where
  display = displayLaTeX

max_stars_repo_name: mainland/hspiral

max_issues_repo_name: mainland/hspiral

max_forks_repo_name: mainland/hspiral


text:
{-|
Module      : Qubism.QGate
Description : Types and functions for quantum gates
Copyright   : (c) Keith Pearson, 2018
License     : MIT
Maintainer  : keith@qubitrot.org
-}

{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE KindSignatures      #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeOperators       #-}

module Qubism.QGate
  ( QGate
  , (#>)
  , gate
  , ident
  , pauliX
  , pauliY
  , pauliZ
  , hadamard
  , unitary
  , cnot
  , controlled
  , ifBit
  , kronecker
  , onJust
  , onEvery
  , onRange
  ) where

-- For dependent typing
import GHC.TypeLits
import Data.Singletons
import Data.Singletons.TypeLits
import Data.Finite

import           Data.Complex
import           Data.Monoid
import           Control.Monad.Trans.State.Strict
import           Numeric.LinearAlgebra ((><))
import qualified Numeric.LinearAlgebra as LA

import Qubism.Algebra
import Qubism.StateVec
import Qubism.CReg

newtype QGate (n :: Nat) =
  UnsafeMkQGate (LA.Matrix C) -- must be a (2^n) x (2^n) unitary matrix
  deriving (Show)

-- | "Close enough" equality testing
instance Eq (QGate n) where
  (UnsafeMkQGate a) == (UnsafeMkQGate b) =
    LA.norm_2 (a - b) < 0.000001

instance KnownNat n => Semigroup (QGate n) where
  (UnsafeMkQGate a) <> (UnsafeMkQGate b) = UnsafeMkQGate $ a LA.<> b

instance KnownNat n => Monoid (QGate n) where
  mempty = ident

instance KnownNat n => VectorSpace (QGate n) where
  zero = UnsafeMkQGate . (l><l) $ repeat 0 where l = internalLen (sing :: Sing n)
  z                 .: (UnsafeMkQGate b) = UnsafeMkQGate $ LA.scalar z * b
  (UnsafeMkQGate a) +: (UnsafeMkQGate b) = UnsafeMkQGate $ a + b
  neg (UnsafeMkQGate a) = UnsafeMkQGate $ -a

instance KnownNat n => Algebra (QGate n) where
  (*:) = (<>)

internalLen :: (KnownNat n, Num a) => Sing n -> a
internalLen = (2 ^) . fromIntegral . fromSing

-- | Apply a quantum gate to a quantum register. Note that this
-- operator conflicts with the one from Numeric.LinearAlgebra.
infixr 5 #>
(#>) :: QGate n -> StateVec n -> StateVec n
(#>) (UnsafeMkQGate m) (UnsafeMkStateVec v) = UnsafeMkStateVec $ m LA.#> v

-- | Helper function to make state computations nicer
gate :: Monad m => QGate n -> StateT (StateVec n) m ()
gate g = state $ \qr -> ((), g #> qr)

ident :: forall n . KnownNat n => QGate n
ident = UnsafeMkQGate . LA.ident $ internalLen (sing :: Sing n)

-- | Also known as a NOT gate
pauliX :: QGate 1
pauliX = UnsafeMkQGate $
  (2><2) [0 :+ 0, 1 :+ 0,
          1 :+ 0, 0 :+ 0]

pauliY :: QGate 1
pauliY = UnsafeMkQGate $
  (2><2) [0 :+ 0, 0 :+ (-1),
          0 :+ 1, 0 :+ 0]

pauliZ :: QGate 1
pauliZ = UnsafeMkQGate $
  (2><2) [1 :+ 0, 0 :+ 0,
          0 :+ 0, (-1) :+ 0]

hadamard :: QGate 1
hadamard = UnsafeMkQGate $ 1 / sqrt 2 *
  (2><2) [1 :+ 0, 1 :+ 0,
          1 :+ 0, (-1) :+ 0]

-- | An arbitrary element of SU(2),
-- U(theta, phi, lambda)
unitary :: Double -> Double -> Double -> QGate 1
unitary theta phi lambda = UnsafeMkQGate $
  (2><2) [a, b, c, d]
  where a =  cis (phi+lambda/2) * ( cos (theta/2) :+ 0 )
        b = -cis (phi-lambda/2) * ( sin (theta/2) :+ 0 )
        c =  cis (phi-lambda/2) * ( sin (theta/2) :+ 0 )
        d =  cis (phi+lambda/2) * ( cos (theta/2) :+ 0 )

-- | Controlled not between control and target qubits
cnot :: KnownNat n => Finite n -> Finite n -> QGate n
cnot c t = controlled c . onJust t $ pauliX

-- | Transform an arbitrary QGate into a controled one depending on qubit i
controlled :: forall n . KnownNat n => Finite n -> QGate n -> QGate n
controlled finite (UnsafeMkQGate m) =
  UnsafeMkQGate $ (m <> projection) + LA.ident (2^n) - projection
  where
    projection = LA.diag $ (2^n) LA.|> fmap f [0..]
    f j = fromIntegral $ j `quot` 2^(n-i-1) `mod` 2
    n   = fromIntegral $ fromSing (sing :: Sing n)
    i   = fromIntegral $ getFinite finite

-- | If the given Bit is One then apply the gate, otherwise
--  do nothing
ifBit :: KnownNat n => Bit -> QGate n -> QGate n
ifBit b g = if (b == One) then g else ident

-- | The tensor product of a QGate's A and B is a QGate that acts as A on the
-- first n qubits and B on the rest. In the computational basis this is simply
-- the kronecker product of matricies.
kronecker :: QGate n -> QGate m -> QGate (m+n)
kronecker (UnsafeMkQGate a) (UnsafeMkQGate b) =
  UnsafeMkQGate $ LA.kronecker a b

-- | Promote a 1-qubit gate to an n-qubit gate with the original gate acting
-- on qubit i. Other qubits are unaffected.
onJust :: forall n . KnownNat n => Finite n -> QGate 1 -> QGate n
onJust i (UnsafeMkQGate m) = UnsafeMkQGate $  -- it should be possible to do
  pre `LA.kronecker` m `LA.kronecker` post    -- this with just QGate's
  where pre  = LA.ident $ 2^j                 -- kronecker, but I can't get
        post = LA.ident $ 2^(n-j-1)           -- the types to check.
        j    = fromIntegral $ getFinite i
        n    = fromIntegral $ fromSing (sing :: Sing n)

-- | Promote a 1-qubit gate to an n-qubit gate which acts on every qubit
-- identically.
onEvery :: forall n . KnownNat n => QGate 1 -> QGate n
onEvery (UnsafeMkQGate m) = UnsafeMkQGate $ iterate (LA.kronecker m) m !! (n-1)
  where n = fromIntegral $ fromSing (sing :: Sing n)

-- | Promote a 1-qubit gate to an n-qubit gate which acts on a range of
-- qubits, leaving others untouched
onRange :: forall n . KnownNat n => Finite n -> Finite n -> QGate 1 -> QGate n
onRange f l m = mconcat $ map (\i -> onJust i m) [f..l]

max_stars_repo_name: qubitrot/qubism

max_issues_repo_name: qubitrot/qubism

max_forks_repo_name: qubitrot/qubism


text:
{-# LANGUAGE DataKinds                 #-}
{-# LANGUAGE FlexibleContexts          #-}
{-# LANGUAGE GADTs                     #-}
{-# LANGUAGE OverloadedStrings         #-}
{-# LANGUAGE PolyKinds                 #-}
{-# LANGUAGE ScopedTypeVariables       #-}
{-# LANGUAGE TypeApplications          #-}
{-# LANGUAGE TypeOperators             #-}
{-# LANGUAGE QuasiQuotes               #-}
{-# OPTIONS_GHC  -fplugin=Polysemy.Plugin  #-}
--{-# LANGUAGE AllowAmbiguousTypes       #-}

import qualified Control.Foldl                 as FL

import qualified Data.List                     as L
import qualified Data.Map                      as M
import qualified Data.Vector                   as VB
import qualified Data.Text                     as T

import qualified Text.Printf                   as PF

import qualified Data.Vinyl                    as V
import qualified Frames                        as F
import qualified Frames.CSV                    as F

import qualified Pipes                         as P
import qualified Pipes.Prelude                 as P

import qualified Numeric.LinearAlgebra         as LA

import qualified Frames.Visualization.VegaLite.Data
                                               as FV
import qualified Frames.Visualization.VegaLite.StackedArea
                                               as FV
import qualified Frames.Visualization.VegaLite.LineVsTime
                                               as FV
import qualified Frames.Visualization.VegaLite.ParameterPlots
                                               as FV
import qualified Frames.Visualization.VegaLite.Correlation
                                               as FV

import qualified Frames.Transform              as FT
import qualified Frames.Folds                  as FF
import qualified Frames.MapReduce              as MR
import qualified Frames.Enumerations           as FE
import qualified Frames.Utils                  as FU

import qualified Knit.Report                   as K
import           Polysemy.Error                 ( Error )

import           Data.String.Here               ( here )

import           BlueRipple.Configuration
import           BlueRippleUtilitie.KnitUtils
import           BlueRipple.Data.DataFrames
import           BlueRipple.Data.MRP
import qualified BlueRipple.Model.TurnoutAdjustment
                                               as TA

templateVars = M.fromList
  [ ("lang"     , "English")
  , ("author"   , "Adam Conner-Sax & Frank David")
  , ("pagetitle", "Preference Model & Predictions")
--  , ("tufte","True")
  ]

{- TODO
1. Why are rows still being dropped?  Which col is missing?  In general, write a diagnostic for showing what is missing...
Some answers:
The boring:  CountyFIPS.  I don't want this anyway.  Dropped
The weird:  Missing voted_rep_party.  These are non-voters or voters who didn't vote in house race.  A bit more than 1/3 of
survey responses.  Maybe that's not weird??
-}

main :: IO ()
main = do
  let template = K.FromIncludedTemplateDir "mindoc-pandoc-KH.html"
--  let template = K.FullySpecifiedTemplatePath "pandoc-templates/minWithVega-pandoc.html"
  pandocWriterConfig <- K.mkPandocWriterConfig template
                                               templateVars
                                               K.mindocOptionsF
  eitherDocs <-
    K.knitHtmls (Just "MRP_Basics.Main") K.logAll pandocWriterConfig $ do
      K.logLE K.Info "Loading data..."
      let csvParserOptions =
            F.defaultParser { F.quotingMode = F.RFC4180Quoting ' ' }
          tsvParserOptions = csvParserOptions { F.columnSeparator = "\t" }
          preFilterYears   = FU.filterOnMaybeField @CCESYear (`L.elem` [2016])
      ccesMaybeRecs <- loadToMaybeRecs @CCES_MRP_Raw @(F.RecordColumns CCES)
        tsvParserOptions
        preFilterYears
        ccesTSV
      ccesFrame <-
        fmap transformCCESRow
          <$> maybeRecsToFrame fixCCESRow (const True) ccesMaybeRecs
      let firstFew = take 4 $ FL.fold FL.list ccesFrame
      K.logLE K.Diagnostic
        $  "ccesFrame (first 4 rows):\n"
        <> (T.pack $ show firstFew)
      K.logLE K.Info "Inferring..."
      K.logLE K.Info "Knitting docs..."
  case eitherDocs of
    Right namedDocs ->
      K.writeAllPandocResultsWithInfoAsHtml "reports/html/MRP_Basics" namedDocs
    Left err -> putStrLn $ "pandoc error: " ++ show err

max_stars_repo_name: blueripple/preference-model

max_issues_repo_name: blueripple/preference-model

max_forks_repo_name: blueripple/preference-model


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

module BayesStack.Models.Topic.CitationInfluence
  ( -- * Primitives
    NetData(..)
  , netData
  , HyperParams(..)
  , MState(..)
  , CitedUpdateUnit
  , CitingUpdateUnit
  , ItemSource(..)
  , CitedNode(..), CitedNodeItem(..)
  , CitingNode(..), CitingNodeItem(..)
  , Citing(..), Cited(..)
  , Item(..), Topic(..), NodeItem(..), Node(..), Arc(..)
  , setupNodeItems
    -- * Initialization
  , verifyNetData, cleanNetData
  , ModelInit
  , randomInitialize
  , model
  , updateUnits
    -- * Diagnostics
  , modelLikelihood
  , influence
  , arcTopicMixture
  ) where

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

import           Prelude hiding (mapM_, sum)

import           Data.Set (Set)
import qualified Data.Set as S

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

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

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

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

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

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

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

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

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

data HyperParams = HyperParams
                   { alphaPsi         :: Double
                   , alphaLambda      :: Double
                   , alphaPhi         :: Double
                   , alphaOmega       :: Double
                   , alphaGammaShared :: Double
                   , alphaGammaOwn    :: Double
                   }
                 deriving (Show, Eq, Generic)
instance Binary HyperParams

data NetData = NetData { dHypers             :: !(HyperParams)
                       , dArcs               :: !(Set Arc)
                       , dItems              :: !(Set Item)
                       , dTopics             :: !(Set Topic)
                       , dNodeItems          :: !(Map NodeItem (Node, Item))
                       , dCitingNodes        :: !(Map CitingNode (Set CitedNode))
                         -- ^ Maps each citing node to the set of nodes cited by it
                       , dCitedNodes         :: !(Map CitedNode (Set CitingNode))
                         -- ^ Maps each cited node to the set of nodes citing it
                       }
              deriving (Show, Eq, Generic)
instance Binary NetData

netData :: HyperParams -> Set Arc -> Map NodeItem (Node,Item) -> Set Topic -> NetData
netData hypers arcs nodeItems topics =
    NetData { dHypers       = hypers
            , dArcs         = arcs
            , dItems        = foldMap (S.singleton . snd) $ M.elems nodeItems
            , dTopics       = topics
            , dNodeItems    = nodeItems
            , dCitingNodes  = M.unionsWith S.union
                              $ map (\(Arc a b)->M.singleton a $ S.singleton b)
                              $ S.toList arcs
            , dCitedNodes   = M.unionsWith S.union
                              $ map (\(Arc a b)->M.singleton b $ S.singleton a)
                              $ S.toList arcs
            }

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

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

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

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

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

type CitedModelInit = Map CitedNodeItem (Setting CitedUpdateUnit)
type CitingModelInit = Map CitingNodeItem (Setting CitingUpdateUnit)
data ModelInit = ModelInit CitedModelInit CitingModelInit
               deriving (Show)

randomInitializeCited :: NetData -> CitedModelInit -> RVar CitedModelInit
randomInitializeCited d init = execStateT doInit init
    where doInit = let unset = citedNodeItems `S.difference` M.keysSet init
                       citedNodeItems = S.map Cited (M.keysSet (dNodeItems d))
                   in mapM_ (randomInitCitedUU d) (S.toList unset)

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

randomInitCitedUU :: NetData -> CitedNodeItem -> StateT CitedModelInit RVar ()
randomInitCitedUU d ni = do
    t' <- lift $ randomElement $ toList $ dTopics d
    modify' $ M.insert ni t'

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

randomInitCitingUU :: NetData -> CitingNodeItem -> StateT CitingModelInit RVar ()
randomInitCitingUU d cni@(Citing ni) =
    let (n,_) = dNodeItems d M.! ni
    in case dCitingNodes d M.! Citing n of
           a | S.null a -> do
               t <- lift $ randomElement $ toList $ dTopics d
               modify' $ M.insert cni $ OwnSetting t

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

randomInitialize :: NetData -> RVar ModelInit
randomInitialize d =
    ModelInit <$> randomInitializeCited d M.empty <*> randomInitializeCiting d M.empty

model :: NetData -> ModelInit -> MState
model d (ModelInit citedInit citingInit) =
    let citingNodes = M.keys $ dCitingNodes d
        HyperParams {..} = dHypers d
        s = MState { -- Citing model
                     stPsis = let dist n = case toList $ dCitingNodes d M.! n of
                                               []    -> M.empty
                                               nodes -> M.singleton n
                                                        $ symDirMulti alphaPsi nodes
                              in foldMap dist citingNodes
                   , stPhis = let dist = symDirMulti alphaPhi (toList $ dItems d)
                              in foldMap (\t->M.singleton t dist) $ dTopics d
                   , stGammas = let dist = multinom [ (Shared, alphaGammaShared)
                                                    , (Own, alphaGammaOwn) ]
                                in foldMap (\t->M.singleton t dist) citingNodes
                   , stOmegas = let dist = symDirMulti alphaOmega (toList $ dTopics d)
                                in foldMap (\t->M.singleton t dist) citingNodes
                   , stCiting = M.empty

                   -- Cited model
                   , stLambdas = let dist = symDirMulti alphaLambda (toList $ dTopics d)
                                 in foldMap (\t->M.singleton t dist) $ M.keys $ dCitedNodes d
                   , stT' = M.empty
                   }

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

        initCitedUU :: CitedUpdateUnit -> State MState ()
        initCitedUU uu = do
            let err = error $ "CitationInference: Initial value for "++show uu++" not given"
                s = maybe err id $ M.lookup (uuNI' uu) citedInit
            modify' $ setCitedUU uu (Just s)

    in execState (do mapM_ initCitingUU $ citingUpdateUnits d
                     mapM_ initCitedUU $ citedUpdateUnits d
                 ) s

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

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

data MState = MState { -- Citing model state
                       stGammas   :: !(Map CitingNode (Multinom Int ItemSource))
                     , stOmegas   :: !(Map CitingNode (Multinom Int Topic))
                     , stPsis     :: !(Map CitingNode (Multinom Int CitedNode))
                     , stPhis     :: !(Map Topic (Multinom Int Item))

                     , stCiting   :: !(Map CitingNodeItem CitingSetting)

                     -- Cited model state
                     , stLambdas  :: !(Map CitedNode (Multinom Int Topic))

                     , stT'       :: !(Map CitedNodeItem Topic)
                     }
            deriving (Show, Generic)
instance Binary MState

modelLikelihood :: MState -> Probability
modelLikelihood model =
    product $ map likelihood (M.elems $ stGammas model)
           ++ map likelihood (M.elems $ stPhis model)
           ++ map likelihood (M.elems $ stLambdas model)
           ++ map likelihood (M.elems $ stOmegas model)
           ++ map likelihood (M.elems $ stPsis model)

-- | Mixture of the topics of an edge
arcTopicMixture :: NetData -> MState -> Arc -> Topic -> Probability
arcTopicMixture nd m (Arc d c) t =
    let itemObs = zip (itemsOfCitingNode nd d) (repeat 1)
        phi = stPhis m M.! t
        lambda = stLambdas m M.! c
    in realToFrac (sampleProb lambda t) * obsProb phi itemObs

-- | Geometric mean
geomMean :: V.Vector (Log Double) -> Log Double
geomMean = Exp . mean . V.map ln

-- | The geometric mean of the probabilities of a collection of items under a
-- given topic mixture.
topicCompatibility :: MState -> [Item] -> Multinom Int Topic -> Probability
topicCompatibility m items lambda =
    geomMean $ V.fromList
            $ do t <- toList $ dmDomain lambda
                 let phi = stPhis m M.! t
                     itemObs = zip items (repeat 1)
                 return $ realToFrac (sampleProb lambda t) * obsProb phi itemObs

-- | Normalized compatibilities with a set of items
topicCompatibilities :: (Functor f, Foldable f)
                     => MState -> [Item] -> f (Multinom Int Topic) -> f Probability
topicCompatibilities m items topics =
    let scores = fmap (topicCompatibility m items) topics
    in fmap (/sum scores) scores

-- | The influence of cited nodes on a citing node.
influence :: NetData -> MState -> CitingNode -> Map CitedNode Probability
influence d m u =
    let lambdas = foldMap (\f->M.singleton f $ stLambdas m M.! f)
                 $ S.toList $ dCitingNodes d M.! u
    in topicCompatibilities m (itemsOfCitingNode d u) lambdas

-- Cited update unit (LDA-like)
data CitedUpdateUnit = CitedUpdateUnit { uuNI' :: !CitedNodeItem
                                       , uuN'  :: !CitedNode
                                       , uuX' :: !Item
                                       }
                     deriving (Show, Generic)
instance Binary CitedUpdateUnit

instance UpdateUnit CitedUpdateUnit where
    type ModelState CitedUpdateUnit = MState
    type Setting CitedUpdateUnit = Topic
    fetchSetting uu ms = maybe (error $ "Update unit "++show uu++" has no setting") id
                         $ M.lookup (uuNI' uu) (stT' ms)
    evolveSetting ms uu = categorical $ citedFullCond (setCitedUU uu Nothing ms) uu
    updateSetting uu _ s' = setCitedUU uu (Just s') . setCitedUU uu Nothing

citedProb :: MState -> CitedUpdateUnit -> Setting CitedUpdateUnit -> Double
citedProb st (CitedUpdateUnit {uuN'=n', uuX'=x'}) t =
    let lambda = stLambdas st M.! n'
        phi = stPhis st M.! t
    in realToFrac $ sampleProb lambda t * sampleProb phi x'

citedUpdateUnits :: NetData -> [CitedUpdateUnit]
citedUpdateUnits d =
    map (\(ni',(n',x'))->CitedUpdateUnit { uuNI'      = Cited ni'
                                         , uuN'       = Cited n'
                                         , uuX'       = x'
                                         }
        ) $ M.assocs
          $ M.filter (\(n,i)->Cited n `M.member` dCitedNodes d)
          $ dNodeItems d

setCitedUU :: CitedUpdateUnit -> Maybe Topic -> MState -> MState
setCitedUU uu@(CitedUpdateUnit {uuN'=n', uuNI'=ni', uuX'=x'}) setting ms =
    let t' = maybe (fetchSetting uu ms) id setting
        set = maybe Unset (const Set) setting
    in ms { stLambdas = M.adjust (setMultinom set t') n' (stLambdas ms)
          , stPhis = M.adjust (setMultinom set x') t' (stPhis ms)
          , stT' = case setting of Just _  -> M.insert ni' t' $ stT' ms
                                   Nothing -> stT' ms
          }

citedFullCond ::MState -> CitedUpdateUnit -> [(Double, Topic)]
citedFullCond ms uu = do
    t <- M.keys $ stPhis ms
    return (citedProb ms uu t, t)

-- Citing Update unit (Shared Taste-like)
data CitingUpdateUnit = CitingUpdateUnit { uuNI    :: !CitingNodeItem
                                         , uuN     :: !CitingNode
                                         , uuX     :: !Item
                                         , uuCites :: !(Set CitedNode)
                                         }
                      deriving (Show, Generic)
instance Binary CitingUpdateUnit

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

citingUpdateUnits :: NetData -> [CitingUpdateUnit]
citingUpdateUnits d =
    map (\(ni,(n,x))->CitingUpdateUnit { uuNI      = ni
                                       , uuN       = n
                                       , uuX       = x
                                       , uuCites   = dCitingNodes d M.! n
                                       }
        ) $ M.assocs $ dCitingNodeItems d

citingProb :: MState -> CitingUpdateUnit -> Setting CitingUpdateUnit -> Double
citingProb st (CitingUpdateUnit {uuN=n, uuX=x}) setting =
    let gamma = stGammas st M.! n
        omega = stOmegas st M.! n
        psi = stPsis st M.! n
    in case setting of
        SharedSetting t c -> let phi = stPhis st M.! t
                                 lambda = stLambdas st M.! c
                             in sampleProb gamma Shared
                              * sampleProb psi c
                              * sampleProb lambda t
                              * sampleProb phi x
        OwnSetting t      -> let phi = stPhis st M.! t
                             in sampleProb gamma Own
                              * sampleProb omega t
                              * sampleProb phi x

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

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

setCitingUU :: CitingUpdateUnit -> Maybe (Setting CitingUpdateUnit) -> MState -> MState
setCitingUU uu@(CitingUpdateUnit {uuNI=ni, uuN=n, uuX=x}) setting ms =
    let set = maybe Unset (const Set) setting
        ms' = case maybe (fetchSetting uu ms) id  setting of
            SharedSetting t c -> ms { stPsis = M.adjust (setMultinom set c) n $ stPsis ms
                                    , stLambdas = M.adjust (setMultinom set t) c $ stLambdas ms
                                    , stPhis = M.adjust (setMultinom set x) t $ stPhis ms
                                    , stGammas = M.adjust (setMultinom set Shared) n $ stGammas ms
                                    }
            OwnSetting t      -> ms { stOmegas = M.adjust (setMultinom set t) n $ stOmegas ms
                                    , stPhis = M.adjust (setMultinom set x) t $ stPhis ms
                                    , stGammas = M.adjust (setMultinom set Own) n $ stGammas ms
                                    }
    in ms' { stCiting = M.alter (const setting) ni $ stCiting ms' }

max_stars_repo_name: bgamari/bayes-stack

max_issues_repo_name: bgamari/bayes-stack

max_forks_repo_name: bgamari/bayes-stack


text:
module PeekyBlinders.Prelude
  ( module Exports,
    showBin,
  )
where

import Control.Applicative as Exports hiding (WrappedArrow (..))
import Control.Arrow as Exports hiding (first, second)
import Control.Category as Exports
import Control.Concurrent as Exports
import Control.Exception as Exports
import Control.Monad as Exports hiding (fail, forM, forM_, mapM, mapM_, msum, sequence, sequence_)
import Control.Monad.Fail as Exports
import Control.Monad.Fix as Exports hiding (fix)
import Control.Monad.IO.Class as Exports
import Control.Monad.ST as Exports
import Data.Bifunctor as Exports
import Data.Bits as Exports
import Data.Bool as Exports
import Data.ByteString as Exports (ByteString)
import Data.ByteString.Short as Exports (ShortByteString)
import Data.Char as Exports
import Data.Coerce as Exports
import Data.Complex as Exports
import Data.Data as Exports
import Data.Dynamic as Exports
import Data.Either as Exports
import Data.Fixed as Exports
import Data.Foldable as Exports hiding (toList)
import Data.Function as Exports hiding (id, (.))
import Data.Functor as Exports
import Data.Functor.Compose as Exports
import Data.Functor.Contravariant as Exports
import Data.IORef as Exports
import Data.Int as Exports
import Data.Ix as Exports
import Data.List as Exports hiding (all, and, any, concat, concatMap, elem, find, foldl, foldl', foldl1, foldr, foldr1, isSubsequenceOf, mapAccumL, mapAccumR, maximum, maximumBy, minimum, minimumBy, notElem, or, product, sortOn, sum, uncons)
import Data.List.NonEmpty as Exports (NonEmpty (..))
import Data.Maybe as Exports
import Data.Monoid as Exports hiding (Alt, (<>))
import Data.Ord as Exports
import Data.Proxy as Exports
import Data.Ratio as Exports
import Data.STRef as Exports
import Data.Semigroup as Exports hiding (First (..), Last (..))
import Data.String as Exports
import Data.Text as Exports (Text)
import Data.Traversable as Exports
import Data.Tuple as Exports
import Data.Unique as Exports
import Data.Version as Exports
import Data.Void as Exports
import Data.Word as Exports
import Debug.Trace as Exports
import Foreign as Exports hiding (void)
import GHC.Conc as Exports hiding (orElse, threadWaitRead, threadWaitReadSTM, threadWaitWrite, threadWaitWriteSTM, withMVar)
import GHC.Exts as Exports (IsList (..), groupWith, inline, lazy, sortWith)
import GHC.Generics as Exports (Generic)
import GHC.IO.Exception as Exports
import GHC.OverloadedLabels as Exports
import Numeric as Exports
import System.Environment as Exports
import System.Exit as Exports
import System.IO as Exports
import System.IO.Error as Exports
import System.IO.Unsafe as Exports
import System.Mem as Exports
import System.Mem.StableName as Exports
import System.Timeout as Exports
import Text.ParserCombinators.ReadP as Exports (ReadP, readP_to_S, readS_to_P)
import Text.ParserCombinators.ReadPrec as Exports (ReadPrec, readP_to_Prec, readPrec_to_P, readPrec_to_S, readS_to_Prec)
import Text.Printf as Exports (hPrintf, printf)
import Text.Read as Exports (Read (..), readEither, readMaybe)
import Unsafe.Coerce as Exports
import Prelude as Exports hiding (all, and, any, concat, concatMap, elem, fail, foldl, foldl1, foldr, foldr1, id, mapM, mapM_, maximum, minimum, notElem, or, product, sequence, sequence_, sum, (.))

showBin :: (Integral a, Show a) => a -> String
showBin a =
  showIntAtBase 2 intToDigit a ""

max_stars_repo_name: nikita-volkov/peeky-blinders

max_issues_repo_name: nikita-volkov/peeky-blinders

max_forks_repo_name: nikita-volkov/peeky-blinders


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

-- |
-- Module      :  Test.Opcount
-- Copyright   :  (c) 2016-2020 Drexel University
-- License     :  BSD-style
-- Maintainer  :  mainland@drexel.edu

module Test.Opcount (
    opCountTests,
    splitRadixOpcountTests,
    difOpcountTests,
    opcountRegressionTests
  ) where

import Control.Monad (mzero)
import Data.Complex
import Data.Typeable (Typeable)
import Test.HUnit (Assertion,
                   (@?=),
                   assertEqual)
import Test.Hspec

import qualified Data.FlagSet as FS

import Spiral (Config(..))
import Spiral.Config
import Spiral.Driver
import Spiral.Exp
import Spiral.FFT.CooleyTukey
import Spiral.Monad
import Spiral.OpCount
import Spiral.SPL
import Spiral.SPL.Run
import Spiral.Search
import Spiral.Search.FFTBreakdowns
import Spiral.Search.OpCount

import Test.Instances ()

opCountTests :: Spec
opCountTests =
    describe "Opcount" $ do
    splitRadixOpcountTests
    difOpcountTests

-- The number of multiplies and additions for spit radix decomposition of size n
-- when using three-multiply form of complex multiply. Taken from Table II of
-- Heideman and Burrus.
splitRadixOpcounts :: [(Int, Int, Int)]
splitRadixOpcounts = [ (4,    0,     16)
                     , (8,    4,     52)
                     , (16,   20,    148)
                     , (32,   68,    388)
                     , (64,   196,   964)
                     , (128,  516,   2308)
                     , (256,  1284,  5380)
                     --, (512,  3076,  12292)
                     --, (1024, 7172,  27652)
                     --, (2048, 16388, 61444)
                     --, (4096, 36868, 135172)
                     ]

splitRadixOpcountTests :: Spec
splitRadixOpcountTests =
    describe "Split radix opcounts" $
    sequence_ [mkTest n (muls + adds) | (n, muls, adds) <- splitRadixOpcounts]
  where
    mkTest :: Int -> Int -> Spec
    mkTest n nops =
      mkOpCountTest ("Split radix " ++ show n) fs nops (runSearch () f (DFT n))
      where
        fs :: [DynFlag]
        fs = [ StoreIntermediate
             , SplitComplex
             , CSE
             , Rewrite
             ]

    f :: (Typeable a, Typed a, MonadSpiral m)
      => SPL (Exp a)
      -> S s m (SPL (Exp a))
    f (F n w) = splitRadixBreakdown n w
    f _       = mzero

-- The DIF form should have the same operation count as split radix when the
-- DifRewrite flag is enabled.
difOpcountTests :: Spec
difOpcountTests =
    describe "DIF opcounts" $
    sequence_ [mkTest n (muls + adds) | (n, muls, adds) <- splitRadixOpcounts]
  where
    mkTest :: Int -> Int -> Spec
    mkTest n nops =
      mkOpCountTest ("DIF " ++ show n) fs nops (return $ dif n)
      where
        fs :: [DynFlag]
        fs = [ StoreIntermediate
             , SplitComplex
             , CSE
             , Rewrite
             , DifRewrite
             ]

mkOpCountTest :: String
              -> [DynFlag]
              -> Int
              -> Spiral (SPL (Exp (Complex Double)))
              -> Spec
mkOpCountTest desc fs nops gen =
    it desc $ do
    ops <- runSpiralWith mempty $ withOpcountFlags fs $ do
           f   <- gen
           toProgram "f" (Re f) >>= countProgramOps
    return $ mulOps ops + addOps ops @?= nops

withOpcountFlags :: MonadConfig m => [DynFlag] -> m a -> m a
withOpcountFlags fs =
    localConfig $ \env -> env { dynFlags  = FS.fromList fs
                              , maxUnroll = 256
                              }

opcountRegressionTests :: Int -> Spec
opcountRegressionTests max_size = do
    text <- runIO $ readFile rEGRESSION_FILE
    let opcounts = (map parseCSV . drop 1 . lines) text
    mapM_ test [(n,totalOps,mulOps,addOps) | [n,totalOps,mulOps,addOps] <- opcounts, n <= max_size]
  where
    rEGRESSION_FILE :: FilePath
    rEGRESSION_FILE = "benchmark/data/search-opcount.csv"

    parseCSV :: String -> [Int]
    parseCSV = map read . splitOn (== ',')

    test :: (Int,Int,Int,Int) -> Spec
    test (n, allOps0, mulOps0, addOps0) = it ("Regression(" ++ show n ++ ")") $ do
        ops <- runSpiralWith config $ do
               e    <- searchOpCount (Re (DFT n) :: SPL (Exp Double))
               prog <- toProgram "dft" e
               countProgramOps prog
        return $ checkOps (allOps ops, mulOps ops, addOps ops) (allOps0, mulOps0, addOps0)

    checkOps :: (Int, Int, Int) -> (Int, Int, Int) -> Assertion
    checkOps ops@(total, _, _) ops'@(total', _, _)
      | total < total' = assertEqual "IMPROVED!" ops' ops
      | otherwise      = ops @?= ops'

    config :: Config
    config = mempty { dynFlags  = FS.fromList fs
                    , maxUnroll = 256
                    }

    fs :: [DynFlag]
    fs = [ StoreIntermediate
         , SplitComplex
         , CSE
         , Rewrite
         ]

splitOn :: (Char -> Bool) -> String -> [String]
splitOn p s = case dropWhile p s of
                  "" -> []
                  s' -> w : splitOn p s''
                        where (w, s'') = break p s'


text:
module TB.Statistics.Examples (
) where

import Data.Vector
import Statistics.Sample

v1 = fromList [1,1,1,4,8,10,12,9,7,5,3,1,1,1,2,1]
weights = fromList [0.1,0.1,0.1,0.4,0.8,1,1.2,0.9,0.7,0.5,0.3,0.1,0.1,0.1,0.1,0.2,0.1]

stats =
  [
    ("stdDev", stdDev v1)
  , ("skewness", skewness v1)
  , ("range", range v1)
  , ("mean", mean v1)
  , ("welfordMean", welfordMean v1)
  , ("meanWeighted", meanWeighted $ Data.Vector.zip v1 weights)
  , ("harmonicMean", harmonicMean v1)
  , ("geometricMean", geometricMean v1)
  , ("kurtosis", kurtosis v1)
  , ("variance", variance v1)
  , ("varianceUnbiased", varianceUnbiased v1)
  , ("meanVariance fst", fst mv)
  , ("meanVariance snd", snd mv)
  ]
  where
    mv = meanVariance v1

max_stars_repo_name: adarqui/toybox

max_issues_repo_name: adarqui/toybox

max_forks_repo_name: adarqui/toybox


text:
module Numeric.DataFrame
  ( module Numeric.DataFrame.Type
    -- * Simplified type aliases
  , module Numeric.Scalar
  , module Numeric.Vector
  , module Numeric.Matrix
    -- * Functionality
  , module Numeric.DataFrame.SubSpace
  , module Numeric.DataFrame.Contraction
  , module Numeric.Basics
  , module Numeric.Subroutine.Sort
  ) where

import Numeric.DataFrame.Internal.Backend ()

import Numeric.DataFrame.Contraction
import Numeric.DataFrame.SubSpace
import Numeric.DataFrame.Type

import Numeric.Matrix
import Numeric.Scalar
import Numeric.Vector

import Numeric.Basics

import Numeric.Subroutine.Sort (sort, sortBy)

max_stars_repo_name: achirkin/fasttensor

max_issues_repo_name: achirkin/fasttensor

max_forks_repo_name: achirkin/fasttensor