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

module MRP.WWC (post) where

import qualified Control.Foldl                 as FL
import qualified Data.List as L
import qualified Data.Set as S
import qualified Data.Map                      as M
import qualified Data.Text                     as T
import qualified Data.Serialize                as SE

import           Graphics.Vega.VegaLite.Configuration as FV
import qualified Graphics.Vega.VegaLite.Compat as FV
import qualified Frames as F
import qualified Data.Vinyl as V

import qualified Graphics.Vega.VegaLite        as GV
import qualified Knit.Report                   as K
import qualified Polysemy.Error                as P (mapError)
import           Text.Pandoc.Error             as PE
import qualified Text.Blaze.Colonnade          as BC

import           Data.String.Here               ( i )

import qualified Colonnade                     as C
import qualified Text.Blaze.Html5.Attributes   as BHA

import           BlueRipple.Configuration
import           BlueRipple.Utilities.KnitUtils
import           BlueRipple.Utilities.TableUtils

import qualified Numeric.GLM.ProblemTypes      as GLM
import qualified Numeric.GLM.ModelTypes      as GLM
import qualified Numeric.GLM.Bootstrap            as GLM
import qualified Numeric.GLM.Predict            as GLM
import qualified Numeric.GLM.Confidence            as GLM

import qualified Statistics.Types              as ST
import GHC.Generics (Generic)

import BlueRipple.Data.DataFrames
import qualified BlueRipple.Data.DemographicTypes as BR
import qualified BlueRipple.Data.ElectionTypes as ET
import qualified BlueRipple.Model.MRP as BR
import MRP.Common
import BlueRipple.Data.CCES as CCES

import qualified PreferenceModel.Common as PrefModel

brIntro :: T.Text
brIntro = [i|
In our research pieces so far, we've looked only at aggregate data, that
is data which comes from adding together a large number of people: census
counts or election results, for example.  In this piece we look at some
per-person data, namely the CCES survey, which surveys about 60,000 people
in every election cycle. For each survey response,
the CCES includes geographic and demographic information along with opinion about
various political questions, whether the person is registered to vote, and whether they
voted and who for in elections for Governor, House, Senate and President, whenever
each is applicable.

The geographic information makes this interesting.  This allows
us to start estimating a variety of things at the state level, something we couldn't do using
only aggregate data.  We do this using multi-level regression
(the "MR" of "MRP", which stands for Multi-level Regression with Post-stratification), a technique
we explain in more detail [here][MRP:Methods].

Following up on our earlier posts about
[changing voter preference][PrefModel:AcrossTime] and
the [white working class][PrefModel:WWCV] (WWC),
we look again at WWC voters in 2016 and
2018, this time estimating voter preference in each state and looking at
the changes, with a particular focus on the battleground states in the
mid-west: Indiana, Michigan, Ohio, Pennsylvania and Wisconsin.

[MRP:Methods]: <${brGithubUrl (postPath PostMethods)}>
[PrefModel:WWCV]: <${brGithubUrl (PrefModel.postPath PrefModel.PostWWCV)}>
[PrefModel:AcrossTime]: <${brGithubUrl (PrefModel.postPath PrefModel.PostAcrossTime)}>
|]

glmErrorToPandocError :: GLM.GLMError -> PE.PandocError
glmErrorToPandocError x = PE.PandocSomeError $ T.pack $ show x

--type BR.LocationCols = '[StateAbbreviation]
type CatCols = '[]
--type CCESGroup = Proxy (BR.GroupCols Location

post :: (K.KnitOne r, K.Member GLM.RandomFu r, K.Member GLM.Async r)
     => M.Map T.Text T.Text -- state names from state abbreviations
--     -> K.Cached es [F.Record CCES.CCES_MRP]
     -> K.Sem r ()
post stateNameByAbbreviation = P.mapError glmErrorToPandocError $ K.wrapPrefix "Intro" $ return () {- $ do

  K.logLE K.Info $ "Working on Intro post..."
  let isWWC r = (F.rgetField @BR.SimpleRaceC r == BR.White) && (F.rgetField @BR.CollegeGradC r == BR.NonGrad)
      countWWCDemPres2016VotesF = BR.weightedCountFold @ByCCESPredictors @CCES.CCES_MRP @'[CCES.Pres2016VoteParty,CCES.CCESWeightCumulative]
                                  (\r -> (F.rgetField @CCES.Turnout r == CCES.T_Voted)
                                         && (F.rgetField @Year r == 2016)
                                         && isWWC r
                                         && (F.rgetField @CCES.Pres2016VoteParty r `elem` [ET.Republican, ET.Democratic]))
                                  ((== ET.Democratic) . F.rgetField @CCES.Pres2016VoteParty)
                                  (F.rgetField @CCES.CCESWeightCumulative)
      countWWCDemHouseVotesF y = BR.weightedCountFold @ByCCESPredictors @CCES.CCES_MRP @'[CCES.HouseVoteParty,CCES.CCESWeightCumulative]
                                 (\r -> (F.rgetField @CCES.Turnout r == CCES.T_Voted)
                                        && (F.rgetField @Year r == y)
                                        && isWWC r
                                        && (F.rgetField @CCES.HouseVoteParty r `elem` [ET.Republican, ET.Democratic]))
                                 ((== ET.Democratic) . F.rgetField @CCES.HouseVoteParty)
                                 (F.rgetField @CCES.CCESWeightCumulative)
  let --makeTableRows :: K.Sem r [WWCTableRow]
      makeWWCTableRows = do
        ccesFrameAll <- ccesDataLoader -- F.toFrame <$> P.raise (K.useCached ccesRecordListAllCA)
        (mm2016p, rc2016p, ebg2016p, bu2016p, vb2016p, bs2016p) <- BR.inferMR @BR.LocationCols @CatCols @[BR.SimpleAgeC
                                                                                ,BR.SexC
                                                                                ,BR.CollegeGradC
                                                                                ,BR.SimpleRaceC
                                                                                ]
                                                                   countWWCDemPres2016VotesF
                                                                   [GLM.Intercept, GLM.Predictor CCES.P_WWC]
                                                                   ccesPredictor
                                                                   (fmap F.rcast ccesFrameAll)

        (mm2016, rc2016, ebg2016, bu2016, vb2016, bs2016) <- BR.inferMR @BR.LocationCols @CatCols @[BR.SimpleAgeC
                                                                                ,BR.SexC
                                                                                ,BR.CollegeGradC
                                                                                ,BR.SimpleRaceC
                                                                                ]
                                                             (countWWCDemHouseVotesF 2016)
                                                             [GLM.Intercept, GLM.Predictor P_WWC]
                                                             ccesPredictor
                                                             ccesFrameAll

        (mm2018, rc2018, ebg2018, bu2018, vb2018, bs2018) <- BR.inferMR @BR.LocationCols @CatCols @[BR.SimpleAgeC
                                                                                ,BR.SexC
                                                                                ,BR.CollegeGradC
                                                                                ,BR.SimpleRaceC
                                                                                ]
                                                             (countWWCDemHouseVotesF 2018)
                                                             [GLM.Intercept, GLM.Predictor P_WWC]
                                                             ccesPredictor
                                                             ccesFrameAll

        let states = FL.fold FL.set $ fmap (F.rgetField @StateAbbreviation) ccesFrameAll
            toPredict = [("National", M.empty)] <> (fmap (\s -> (s,M.singleton BR.RecordColsProxy (s F.&: V.RNil))) $ S.toList states)
            wwc = M.singleton P_WWC 1
            addPredictions (s, groupMap) = P.mapError glmErrorToPandocError $ do
              let cl = ST.mkCL 0.95
                  ct bs = GLM.BootstrapCI GLM.BCI_Accelerated bs
              p2016p <- GLM.predictWithCI mm2016p (flip M.lookup wwc) (flip M.lookup groupMap) rc2016p ebg2016p bu2016p vb2016p cl (ct bs2016p)
              p2016 <- GLM.predictWithCI mm2016 (flip M.lookup wwc) (flip M.lookup groupMap) rc2016 ebg2016 bu2016 vb2016 cl (ct bs2016)
              p2018 <- GLM.predictWithCI mm2018 (flip M.lookup wwc) (flip M.lookup groupMap) rc2018 ebg2018 bu2018 vb2018 cl (ct bs2018)
              return  $ WWCTableRow s p2016p p2016 p2018
        L.sortOn deltaHouse <$> traverse addPredictions toPredict
  forWWCTable <-  K.retrieveOrMake "mrp/intro/wwcPrefTable" makeWWCTableRows -- cache this so we don't need to redo
  forWWCChart <-  GLM.eitherToSem $ traverse (\tr -> do
                                            let sa = stateAbbr tr
                                            fullState <- maybe (Left $ "Couldn't find " <> sa) Right $ M.lookup sa stateNameByAbbreviation
                                            return (fullState, significantDeltaHouse tr)) $ filter (\tr -> (stateAbbr tr) /= "National") forWWCTable
  brAddMarkDown brIntro
  _ <- K.addHvega Nothing Nothing $ (vlPctStateChoropleth "Significant Change in WWC Dem Voter Preference 2016 to 2018" (FV.ViewConfig 800 400 10) forWWCChart)
  let mwBG = ["IA", "MI", "OH", "PA", "WI"]
  brAddRawHtmlTable
    "WWC Democratic Voter Preference"
    (BHA.class_ "brTable")
    (colWWC_Change $ emphasizeStates mwBG <> emphasizeNational) forWWCTable

  brAddMarkDown brReadMore
-}

data WWCTableRow = WWCTableRow { stateAbbr :: T.Text
                               , pref2016Pres :: (Double, (Double, Double))
                               , pref2016House :: (Double, (Double, Double))
                               , pref2018House :: (Double, (Double, Double))
                               } deriving (Generic)

instance SE.Serialize WWCTableRow

deltaHouse tr = fst (pref2018House tr) - fst (pref2016House tr)
inStates s tr = stateAbbr tr `elem` s

emphasizeStates s = CellStyle (\tr _ -> if inStates s tr then highlightCellBlue else "")
emphasizeNational = CellStyle (\tr _ -> if  stateAbbr tr == "National" then highlightCellPurple else "")

significantGivenCI :: (Double, Double) -> (Double, Double) -> Double
significantGivenCI (loA, hiA) (loB, hiB) =
  let maxLo = max loA loB
      minHi = min hiA hiB
      overlap = max 0 (minHi - maxLo)
      shortest = min (hiA - loA) (hiB - loB)
  in overlap / shortest

significantDeltaHouse tr = if signifWWCT tr < 0.25 then deltaHouse tr else 0

signifWWCT x = significantGivenCI (snd $ pref2016House x) (snd $ pref2018House x)

colWWC_Change :: CellStyle WWCTableRow T.Text -> C.Colonnade C.Headed WWCTableRow BC.Cell
colWWC_Change cas =
  C.headed "State" (toCell cas "State" "State" (textToStyledHtml . stateAbbr))
  <> C.headed "2016 Pres (%)" (toCell cas "2016 Pres" "2016 Pres" (numberToStyledHtml "%2.1f" . (*100) . fst . pref2016Pres))
  <> C.headed "2016 House (%)" (toCell cas "2016 House" "2016 House" (numberToStyledHtml "%2.1f" . (*100) . fst . pref2016House))
  <> C.headed "Low (95%)" (toCell cas "Low (95%)" "Low (95%)" (numberToStyledHtml "%2.1f" . (*100) . fst . snd . pref2016House))
  <> C.headed "High (95%)" (toCell cas "High (95%)" "High (95%)" (numberToStyledHtml "%2.1f" . (*100) . snd . snd . pref2016House))
  <> C.headed "2018 House (%)" (toCell cas "2018 D Pref" "2018 House" (numberToStyledHtml "%2.1f" . (*100) . fst . pref2018House))
  <> C.headed "2018-2016 House" (toCell cas "2018-2016 House" "2018-2016 House" (numberToStyledHtml "%2.1f" . (*100) . deltaHouse))
  <> C.headed "CI Overlap" (toCell cas "CI Overlap" "CI Overlap" (numberToStyledHtml "%2.1f" . (*100) . signifWWCT))

usStatesTopoJSONUrl = "https://cdn.jsdelivr.net/npm/us-atlas@3/states-10m.json"
usStatesAlbersTopoJSONUrl = "https://cdn.jsdelivr.net/npm/us-atlas@3/states-albers-10m.json"

-- Lessons learned
-- The data is the geography
-- transform with lookup (properties.name) to add data
-- to figure out the right lookup name you can just plot the geography (no transform, no encoding) and look at loaded data
-- in vega editor
vlPctStateChoropleth :: Foldable f => T.Text -> FV.ViewConfig -> f (T.Text, Double) -> GV.VegaLite
vlPctStateChoropleth title vc stateData =
  let datGeo = GV.dataFromUrl usStatesTopoJSONUrl [GV.TopojsonFeature "states"]
      datVal = GV.dataFromRows [] $ concat $ fmap (\(s,x) -> GV.dataRow [("state", GV.Str s),("value", GV.Number x)] []) $ FL.fold FL.list stateData
      dataSets = GV.datasets [("stateNums",datVal)]
      projection = GV.projection [GV.PrType GV.AlbersUsa]
      transform = GV.transform . GV.lookup "properties.name" datVal "state" (GV.LuFields ["state","value"])
      mark = GV.mark GV.Geoshape []
      colorEnc = GV.color [GV.MName "value", GV.MmType GV.Quantitative, GV.MScale [GV.SScheme "redyellowgreen" [], GV.SDomain (GV.DNumbers [-0.5,0.5])]]
      tooltip = GV.tooltips [[GV.TName "state", GV.TmType GV.Nominal],[GV.TName "value", GV.TmType GV.Quantitative, GV.TFormat ".0%"]]
      enc = GV.encoding .  colorEnc . tooltip
  in FV.configuredVegaLite vc [FV.title title, datGeo, mark, projection, transform [], enc []]

vlTest :: FV.ViewConfig -> GV.VegaLite
vlTest vc =
  let
    datGeo = GV.dataFromUrl usStatesTopoJSONUrl [GV.TopojsonFeature "states"]
    projection = GV.projection [GV.PrType GV.AlbersUsa]
    mark = GV.mark GV.Geoshape [GV.MFill "lightgrey" ]
  in FV.configuredVegaLite vc [datGeo, mark, projection]

max_stars_repo_name: blueripple/preference-model

max_issues_repo_name: blueripple/preference-model

max_forks_repo_name: blueripple/preference-model


text:
module FFTSerialSpec where

import qualified Clash.Prelude as Clash
import Clash.Prelude (Signal, Vec(..), BitVector, Index, Signed, Unsigned, SFixed, Bit, SNat(..),
                      simulate, simulate_lazy, listToVecTH, KnownNat, pack, unpack, (++#), mealy, mux, bundle, unbundle,
                      HiddenClockResetEnable, System)
import Test.Hspec
import Test.QuickCheck

import qualified Numeric.FFT as FFT
import qualified Data.Complex as C
import Clash.DSP.Complex (Complex, fromComplex, toComplex)
import Clash.DSP.FFT.Serial
import Clash.DSP.FFT.Example
import Clash.DSP.FFT.Twiddle (twiddleFactors)

--serial FFT
spec = describe "Serial FFTs" $ do
    specify "Decimation in time equals known good implementation"      $ property prop_fftSerialDIT
    specify "Decimation in frequency equals known good implementation" $ property prop_fftSerialDIF

approxEqual :: Double -> Double -> Bool
approxEqual x y = abs (x - y) < 0.0001

approxEqualComplex (a C.:+ b) (c C.:+ d) = approxEqual a c && approxEqual b d

twiddles4 :: Vec 4 (Complex Double)
twiddles4 = $(listToVecTH (twiddleFactors 4))

ditInputReorder :: Vec 8 a -> Vec 4 (a, a)
ditInputReorder (a :> b :> c :> d :> e :> f :> g :> h :> Nil) = (a, e) :> (c, g) :> (b, f) :> (d, h) :> Nil

ditOutputReorder :: [(a, a)] -> [a]
ditOutputReorder ((a, b) : (c, d) : (e, f) : (g, h) : _) = a : c : e : g : b : d : f : h : []

difInputReorder :: Vec 8 a -> Vec 4 (a, a)
difInputReorder (a :> b :> c :> d :> e :> f :> g :> h :> Nil) = (a, e) :> (b, f) :> (c, g) :> (d, h) :> Nil

difOutputReorder :: [(a, a)] -> [a]
difOutputReorder ((a, b) : (c, d) : (e, f) : (g, h) : _) = a : e : c : g : b : f : d : h : []

prop_fftSerialDIT :: Vec 8 (C.Complex Double) -> Bool
prop_fftSerialDIT vec = and $ zipWith approxEqualComplex (map toComplex result) (FFT.fft (Clash.toList vec))
    where
    result = ditOutputReorder $ drop 8 $ simulate_lazy @System (fftSerialDIT twiddles4 (pure True)) $ (Clash.toList (ditInputReorder (Clash.map fromComplex vec))) ++ repeat (0, 0)

prop_fftSerialDIF :: Vec 8 (C.Complex Double) -> Bool
prop_fftSerialDIF vec = and $ zipWith approxEqualComplex (map toComplex result) (FFT.fft (Clash.toList vec))
    where
    result = difOutputReorder $ drop 8 $ simulate_lazy @System (fftSerialDIF twiddles4 (pure True)) $ (Clash.toList (difInputReorder (Clash.map fromComplex vec))) ++ repeat (0, 0)

ext: hs

max_stars_repo_name: adamwalker/clash-utils

max_issues_repo_name: adamwalker/clash-utils

max_forks_repo_name: adamwalker/clash-utils


text:
module Universe
where

import Graphics.Rendering.OpenGL as OpenGL

import AObject
import Politics
import Statistics
import Tree
import OpenGLUtils

aobjs :: AObjTree
aobjs = Node (0.0, 0.0)
  [ Leaf $ AObject "Star"       0   (Color4 0.9 0.0 0.0 1.0) 6.0 1.0 0   glVector3Null Nothing
  , Leaf $ AObject "Murphy's"   10  (Color4 0.5 0.5 1.0 1.0) 2.0 1.0 28  glVector3Null (Just "Murphy")
  , Leaf $ AObject "Loki"       250 (Color4 0.0 0.4 0.5 1.0) 4.0 1.0 55  glVector3Null (Just "Harju")
  , Node (30, 115) $
       [Leaf $ AObject "Harju"         30  (Color4 0.6 0.6 0.6 1.0) 9.0 1.0 0  glVector3Null (Just "Harju")
      , Leaf $ AObject "Harju's Moon"  30  (Color4 0.2 0.9 0.6 1.0) 0.8 1.0 25 glVector3Null (Just "Harju")]
  , Leaf $ AObject "Riesenland" 80  (Color4 0.1 0.8 0.8 1.0) 2.0 1.0 230 glVector3Null (Just "Riesenland")
  , Leaf $ AObject "Riesenland" 80  (Color4 0.1 0.8 0.8 1.0) 2.0 1.0 230 glVector3Null (Just "Riesenland")
  , Node (180, 480) $
       [Leaf $ AObject "Natail"     180 (Color4 0.2 0.2 0.9 1.0) 1.0 1.5 60 glVector3Null (Just "Natail")
      , Leaf $ AObject "Mammoth"    0   (Color4 0.3 0.0 0.6 1.0) 1.5 1.0 40 glVector3Null (Just "Natail")]
  ]

relations = mkRelationshipMap relationsList

relationsList =
  [(("Murphy",     "Harju"),      (Peace, -5)),
   (("Murphy",     "Riesenland"), (Peace, 1)),
   (("Murphy",     "Natail"),     (Peace, -3)),
   (("Harju",      "Riesenland"), (Peace, -1)),
   (("Harju",      "Natail"),     (Peace, 2)),
   (("Riesenland", "Natail"),     (Peace, 0)),
   (("Murphy",     "Murphy"),     (Peace, 10)),
   (("Harju",      "Harju"),      (Peace, 10)),
   (("Riesenland", "Riesenland"), (Peace, 10)),
   (("Natail",     "Natail"),     (Peace, 10))]

randomAllegiance :: IO String
randomAllegiance = chooseIO allegiances

allegiances = ["Murphy", "Harju", "Riesenland", "Natail"]

plalleg :: String
plalleg = ""

initialAttitudes :: AttitudeMap
initialAttitudes = nullAttitudes allegiances


text:
module Data.Array.DynamicS.MatMul(matMul) where
import Prelude hiding ((<>))
import GHC.Stack
import Data.Array.DynamicS
import Numeric.LinearAlgebra as N

matMul :: (HasCallStack, N.Numeric a) => Array a -> Array a -> Array a
matMul x y =
  case (shapeL x, shapeL y) of
    ([m, n], [n', o]) | n == n' ->
      let x' = N.reshape n $ toVector x
          y' = N.reshape o $ toVector y
          xy' = x' N.<> y'
          xy = fromVector [m, o] $ N.flatten xy'
      in  xy
    sz -> error $ "matMul: expected two conforming two-dimensional arrays, got shape " ++ show sz

max_stars_repo_name: augustss/orthotope-hmatrix

max_issues_repo_name: augustss/orthotope-hmatrix

max_forks_repo_name: augustss/orthotope-hmatrix


text:
{-# LANGUAGE OverloadedStrings #-}

module Lib.Util where

import           Control.Arrow
import           Data.Aeson
import qualified Data.Complex           as C
import           Data.List              (foldl')
import qualified Data.Map.Strict        as Map
import qualified Data.Text              as T
import qualified Data.Vector            as V
import qualified Data.WAVE              as WAV
import           Database.Persist.Class
import           Database.Persist.Sql
import qualified Lib.STFT               as STFT
import qualified Lib.TWM                as TWM

-- Takes a vector of doubles, the sample frequency, and a name
-- and writes the audio file. Written in 32 bit.
writeWav :: STFT.Signal -> Int -> String -> IO ()
writeWav (STFT.MkSignal vec) sf name =
  let samples = fmap ((:[]) . WAV.doubleToSample . C.realPart) (V.toList vec)
      header = WAV.WAVEHeader 1 sf 32 Nothing
  in  WAV.putWAVEFile name (WAV.WAVE header samples)

-- Takes a Path and returns IO (sampling frequency, Vector signal).
readWav :: String -> IO (Int, V.Vector Double)
readWav path = do
  audio <- WAV.getWAVEFile path
  let header = WAV.waveHeader audio
      samples = WAV.waveSamples audio
      channels = WAV.waveNumChannels header
      sampRate = WAV.waveFrameRate header
  case channels of
    1 -> let sig = V.fromList $ fmap (WAV.sampleToDouble . Prelude.head) samples
          in return (sampRate, sig)
    _ -> error "Should be mono."

makeSignal :: (Int, V.Vector Double) -> (Int, STFT.Signal)
makeSignal (i, vec) = (i, makeSig vec)

makeSig :: V.Vector Double -> STFT.Signal
makeSig vec = STFT.MkSignal $ fmap (\x -> x C.:+ 0) vec

data Note = A | Bb | B | C | Db | D | Eb | E | F | Gb | G | Ab  deriving (Eq, Ord, Show, Enum)

fromChar :: Char -> Note
fromChar c = undefined

instance ToJSON Note where
    toJSON note = object ["note" .= (show note)]

instance PersistField Note where
  toPersistValue n = PersistText (T.pack $ show n)
  fromPersistValue (PersistText b) = undefined
  fromPersistValue x = Left $ "File.hs: When trying to deserialize a Note: expected PersistText, received: " <> (T.pack $ show x)

instance PersistFieldSql Note where
  sqlType _ = SqlString

data NoteCounter = MkCounter {
    counter :: Int,
    note    :: Note
}

twelthRoot :: Double
twelthRoot = 2 ** (1 / 12)

intToFreq :: Int -> TWM.Freq
intToFreq n = TWM.MkFreq $ 440 * twelthRoot ** fromIntegral n

intToNote :: Int -> Note
intToNote n = toEnum (abs $ n `mod` 12)

initCounter :: Note -> NoteCounter
initCounter n = MkCounter {
    counter = 0,
    note = n
}

closestFreq :: TWM.Freq -> [TWM.Freq] -> TWM.Freq
closestFreq f = snd . minimum . map (abs . subtract f &&& id)

type NoteMap = Map.Map TWM.Freq NoteCounter

-- |Look up table from frequency to note
-- Based on https://pages.mtu.edu/~suits/NoteFreqCalcs.html
freqToNoteTable :: NoteMap
freqToNoteTable = Map.fromList $ fmap (\i -> (intToFreq i, initCounter $ intToNote i)) [(-36)..24]

addToCounter :: NoteCounter -> NoteCounter
addToCounter noteCounter = MkCounter {
    counter = counter noteCounter + 1,
    note = note noteCounter
}

instance Eq NoteCounter where
    (MkCounter counter1 _) == (MkCounter counter2 _) = counter1 == counter2

instance Ord NoteCounter where
    (MkCounter counter1 _) `compare` (MkCounter counter2 _) = counter1 `compare` counter2

findf0 :: [TWM.Freq] -> Note
findf0 fs = note $ maximum $ Map.elems reduced
    where
        freqMap = flip closestFreq (Map.keys freqToNoteTable)
        hz = fmap freqMap (filter (/= 0) fs)
        adjust' = flip $ Map.adjust addToCounter
        reduced = foldl' adjust' freqToNoteTable hz

max_stars_repo_name: davlum/haskell-twm

max_issues_repo_name: davlum/haskell-twm

max_forks_repo_name: davlum/haskell-twm


text:
{-# LANGUAGE FlexibleContexts, ConstraintKinds, DataKinds #-}
import GHC.TypeLits
import Numeric.LinearAlgebra.Static
import qualified Numeric.LinearAlgebra as LA
import System.Environment
import Text.Printf (printf)

-- TODO: Learn about statically typed matrices.
-- TODO: Figure out how to avoid converting b/t lists and matrices.
--       (Will basically require slicing and mapping over rows/colmns).
--
--       Solution: extractRows and extractCols after the predictions line.
--       Could get rid of the toColumns entry as well and separately extract.
--       (Not really worth the time.)

-- Create design matrix for polynomial powers.
-- (or, just use Vector Double. But we want to keep generality.)
-- expand :: (Element t, Field t, Num (Vector t)) => Vector t -> Int -> Matrix t
expand :: (KnownNat n, KnownNat d) => R n -> L n d
expand x deg = fromColumns $ map (x^) [0..deg]

regress :: (Element t, Field t, Num (Vector t)) => Vector t -> Vector t -> Int -> Vector t
regress x y deg = (expand x deg) <\> y
--                                  ^ efficient pseudoinverse

evalPoly :: (Element t, Field t, Product t, Num (Vector t)) => Vector t -> Vector t -> Vector t
evalPoly coeff x = expand x (dim coeff - 1) <> coeff

-- sumVector :: (Element t, Field t) => Vector t -> t
squareDiff :: Vector Double -> Vector Double -> Vector Double
squareDiff x y = mapVector (flip (^) 2) (x - y)
sumVector :: Vector Double -> Double
sumVector x    = x <.> ones where ones = konst 1 (dim x)
sumSquareDiff x y = sumVector $ squareDiff x y

-- NB: This requires concrete types due to subroutines called:
-- readFile   :: FilePath -> IO String
-- readMatrix :: String -> Matrix Double
readData :: String -> IO [Vector Double]
readData filename = fmap (toColumns . readMatrix) $ readFile filename

maxDeg = 5
main = do
  args <- getArgs
  if null args
     then putStrLn "Usage: polyreg datafile"
     else do [x,y] <- readData "polyreg_data.txt"
             let models      = map (regress x y) [1..maxDeg]
                 predictions = map (flip evalPoly x) models
                 residuals   = map (squareDiff y) predictions
                 predsview   = fromColumns $ [x, y] ++ predictions
                 residsview  = fromColumns $ [x, y] ++ residuals
                 totalerrs   = map sumVector residuals
                 putFormat   = putStrLn . format "  " (printf "%.2f")

             putStrLn "Modeling predictions"
             putFormat predsview

             putStrLn "Modeling residuals"
             putFormat residsview

             putStrLn "Total errors"
             putFormat $ fromRows [fromList totalerrs]

max_stars_repo_name: kuitang/monad-learner

max_issues_repo_name: kuitang/monad-learner

max_forks_repo_name: kuitang/monad-learner


text:
module MnistLoader (loadData, loadDataWrapper, loadTestWrapper, investigate) where

import Data.IDX
import qualified Data.Vector.Unboxed as V
import Numeric.LinearAlgebra as L

vect n = (n><1)

loadData :: IO ([(Int, V.Vector Double)])
loadData = do
  mData <- decodeIDXFile "data/train-images-idx3-ubyte"
  let d = getSomething mData
  mLabels <- decodeIDXLabelsFile "data/train-labels-idx1-ubyte"
  let l = getSomething mLabels
  return $ getSomething $ labeledDoubleData l d

loadTest :: IO ([(Int, V.Vector Double)])
loadTest = do
  mData <- decodeIDXFile "data/t10k-images.idx3-ubyte"
  let d = getSomething mData
  mLabels <- decodeIDXLabelsFile "data/t10k-labels.idx1-ubyte"
  let l = getSomething mLabels
  return $ getSomething $ labeledDoubleData l d

-- Converts the label component to a vector representation where
-- n'th component is 1
loadDataWrapper :: IO ([(Matrix Float, Matrix Float)])
loadDataWrapper = do
  theData <- loadData
  return $ map (\(l, d) -> (mVect d, makeVector l)) theData
  where
    makeVector l = vect 10 $ (map (\i -> if i == l then 1 else 0) [0..10])
    mVect l = vect (V.length l) (map realToFrac $ V.toList l)

loadTestWrapper :: IO ([(Matrix Float, Matrix Float)])
loadTestWrapper = do
  theData <- loadTest
  return $ map (\(l, d) -> (mVect d, makeVector l)) theData
  where
    makeVector l = vect 10 $ (map (\i -> if i == l then 1 else 0) [0..10])
    mVect l = vect (V.length l) $ map realToFrac $ V.toList l

investigate :: IO ()
investigate = do
    putStrLn "Inspecting File"
    training_data_mb <- decodeIDXFile "data/train-images-idx3-ubyte"
    let td = getSomething training_data_mb

    putStrLn "Image Data"
    putStr "idxType: "
    putStrLn $ show $ idxType td

    putStr "idxDimensions: "
    putStrLn $ show $ idxDimensions td

    putStr "isIDXReal: "
    putStrLn $ show $ isIDXReal td

    putStr "isIDXIntegral: "
    putStrLn $ show $ isIDXIntegral td

    let training_vector = idxDoubleContent td
    putStrLn $ show $ V.length training_vector

    loaded_data <- loadData
    let (l, d) = loaded_data !! 1
    putStrLn $ show $ V.length $ d
    putStrLn $ show $ l

getSomething :: Maybe a -> a
getSomething (Just x) = x
getSomething Nothing = error "Bah"

max_stars_repo_name: madsbuch/haskell-nn

max_issues_repo_name: madsbuch/haskell-nn

max_forks_repo_name: madsbuch/haskell-nn


text:
{-# LANGUAGE BangPatterns, GADTs, StandaloneDeriving, OverloadedStrings #-}
-- |
-- Module      : Criterion
-- Copyright   : (c) 2009-2014 Bryan O'Sullivan
--
-- License     : BSD-style
-- Maintainer  : bos@serpentine.com
-- Stability   : experimental
-- Portability : GHC
--
-- 'Versus'-type bgroup support

module Criterion.Versus
       (
         vscsv
       , versusReports
       , VersusReport(..)
       ) where

import Criterion.Types
import Criterion.IO.Printf (writeCsv)
import Criterion.Monad (Criterion)
import Control.Monad
import Control.Monad.Trans (liftIO)
import Control.Monad.Reader (asks)
import Control.Arrow ((&&&))
import Data.Function (on)
import qualified Data.Csv as Csv
import Data.List (groupBy, sortBy)
import Statistics.Resampling.Bootstrap (Estimate(..))
import Data.Aeson (ToJSON(..), object, encode, (.=))

data VersusReport where
  VersusReport :: (Show l, Ord l) => {
    vsReportDescription :: String
  , vsReportDataPoints :: [l]
  , vsReportData :: [(String, [Estimate])]
  , vsReportIndices :: [((String, l), Int)]
  } -> VersusReport
deriving instance Show VersusReport

{- This instance is actually incomplete -}
instance ToJSON VersusReport where
  toJSON VersusReport{
      vsReportDescription = desc
    , vsReportDataPoints  = dp
    , vsReportData        = d
    } = object [ "name"       .= toJSON desc
               , "dataPoints" .= toJSON (map show dp)
               , "data"       .= toJSON (map mkArr d)]
        where
          mkArr (alg, p) = object [ "alg"  .= toJSON alg
                                  , "data" .= toJSON p]

vscsv :: [VersusReport] -> Criterion ()
vscsv = mapM_ f
  where f VersusReport{
            vsReportDescription = d
          , vsReportData        = r
          , vsReportDataPoints  = p} = do
          file <- asks vsCsvFile
          writeCsv file [d]
          writeCsv file $ "name":(map show p)
          forM_ r $ \(a, m) -> writeCsv file $ a:(map (show . estPoint) m)

versusReports :: [VersusReport] -> [Report] -> [VersusReport]
versusReports vrpts rpts = map (vsReport rpts) vrpts

vsReport :: [Report] -> VersusReport -> VersusReport
vsReport rpts vr@VersusReport{vsReportIndices = indices} =
  vr{vsReportData = map f l}
  where
    valg = fst . fst
    venv = snd . fst
    l    = groupBy' ((==) `on` valg) valg indices
    rpts' = sortBy (compare `on` reportNumber) rpts
    f (alg, idx) = (alg, [anMean . reportAnalysis $ rpts'!!i
                         | (_, i)<-sortBy (compare `on` venv) idx])

groupBy' :: (a -> a -> Bool) -> (a -> b) -> [a] -> [(b, [a])]
groupBy' f g = map (g . head &&& id) . groupBy f


text:
{-# LANGUAGE BangPatterns #-}
module MachineLearning.ClusteringTest
(
  tests
)

where

import Test.Framework (testGroup)
import Test.Framework.Providers.HUnit
import Test.HUnit
import Test.HUnit.Approx
import Test.HUnit.Plus

import MachineLearning.Types
import qualified Data.Vector as V
import qualified Numeric.LinearAlgebra as LA
import Numeric.LinearAlgebra ((><))
import qualified Control.Monad.Random as RndM
import MachineLearning.Clustering

x1 :: Matrix
x1 = (7><5) [ 1, 2, 3, 4, 5
            , 1, 2, 3, 4, 5
            , 1, 2, 3, 4, 5
            , 7, 4, 3, 2, 1
            , 7, 4, 3, 2, 1
            , 1, 2, 3, 4, 5
            , 1, 2, 3, 4, 5]

x2 :: Matrix
x2 = (7><5) [ 1.1, 2, 3, 4, 5
            , 1, 2, 4, 4, 5
            , 1, 2, 3, 4, 5
            , 5, 4, 3, 2, 1
            , 7, 4, 3, 2, 1
            , 1, 2, 3, 4, 5
            , 0.5, 2, 3, 4, 5]

testKmeans x k expectedK = do
  let gen = RndM.mkStdGen 10171
      clusters = RndM.evalRand (kmeans 10 x k) gen
  assertEqual "number of clusters" expectedK (V.length clusters)

isInDescendingOrder :: [Double] -> Bool
isInDescendingOrder lst = and . snd . unzip $ scanl (\(prev, _) current -> (current, prev-current > (-0.001))) (1/0, True) lst

testDescOrderOfCostValues = do
  let gen = RndM.mkStdGen 10171
      samples = V.fromList $ LA.toRows x1
      (clusters, js) = RndM.evalRand (kmeansIterM samples 3 1) gen
  assertBool "" (isInDescendingOrder js)

tests = [testGroup "kmeans" [
            testCase "perfect clusters, k = 2" $ testKmeans x1 2 2
            , testCase "perfect clusters, k = 3" $ testKmeans x1 3 2
            , testCase "good clusters, k = 2" $ testKmeans x2 2 2
            , testCase "good clusters, k = 3" $ testKmeans x2 3 3
            , testCase "good clusters, k = 4" $ testKmeans x2 4 4
            , testCase "descending order" testDescOrderOfCostValues
            ]
        ]

max_stars_repo_name: aligusnet/mltool

max_issues_repo_name: aligusnet/mltool

max_forks_repo_name: aligusnet/mltool


text:
module TwoMoon.Load where

  import Neural
  import Numeric.LinearAlgebra

  -- Parse two-moon dataset into a list of examples (pairs of a list of feature and a list of labels)
  loadExamples :: IO [Example]
  loadExamples = do
    twomoon <- loadMatrix "data/two_moon.txt"
    let x = vector `map` toLists (twomoon ¿ [0, 1])
    let gt = mkOneHot `map` toList (flatten (twomoon ¿ [2]))
    return $ zip x gt
      where
        mkOneHot :: Double -> Vector Double
        mkOneHot x = if x == 0.0 then vector [1.0, 0.0] else vector [0.0, 1.0]

max_stars_repo_name: ycjungSubhuman/hskmlp

max_issues_repo_name: ycjungSubhuman/hskmlp

max_forks_repo_name: ycjungSubhuman/hskmlp


text:
module Fractal.Utils.Field(
  generateField,
  generateRow,
  generateRowField,
  getFieldLength
)where

import Data.Complex

type Value = Complex Double
type Point = Value
type Row = [Point]
type Field = [Row]
type Step = Double

generateField :: Point -> Point -> Step -> Field
generateField start end step
  | (imagPart start) <= (imagPart end)     = (generateRow start end step):[]
  | otherwise                              = (generateRow start end step) : (generateField newStart end step)
      where newStart = (realPart start) :+ ((imagPart start) - step)

generateRow :: Point -> Point -> Step -> Row
generateRow start end step
  | (realPart start) >= (realPart end)  = [start]
  | otherwise                           = start : (generateRow newStart end step)
      where newStart = ((realPart start) + step) :+ (imagPart start)

generateRowField :: Point -> Point -> Step -> Row
generateRowField start end step = generateRowField' start start end step

generateRowField' :: Point -> Point -> Point -> Step -> Row
generateRowField' mark start end step
  | (imagPart start) <= (imagPart end) = endroutine
  | otherwise = fieldroutine
    where endroutine = if (realPart start) >= (realPart end) then [start] else start : (generateRowField' mark newStart end step)
          fieldroutine = if (realPart start) >= (realPart end) then start:(generateRowField' mark newRowStart end step) else start : (generateRowField' mark newStart end step)
          newStart = ((realPart start) + step) :+ (imagPart start)
          newRowStart = (realPart mark) :+ ((imagPart start) - step)

getFieldLength :: Point -> Point -> Step -> Int
getFieldLength a b step = count a b 0
  where count x y i = if (realPart x) >= (realPart y) then (i+1) else count (newX x) y (i+1)
        newX x = ((realPart x) + step) :+ (imagPart x)

max_stars_repo_name: sigmaticsMUC/FractalHaskell

max_issues_repo_name: sigmaticsMUC/FractalHaskell

max_forks_repo_name: sigmaticsMUC/FractalHaskell


text:
{-# LANGUAGE FlexibleContexts #-}

module Network.Visualizations
( networkHistogram
, weightList
, biasList
) where

import Numeric.LinearAlgebra
import Network.Neuron
import Network.Layer
import Network.Network
import Network.Network.FeedForwardNetwork

import Network.Trainer
import Data.Foldable (foldMap, )
import Graphics.Histogram
import GHC.Float

weightList :: FeedForwardNetwork -> [Double]
weightList n = concat $ map (toList . flatten . weightMatrix) (layers n)

biasList :: FeedForwardNetwork -> [Double]
biasList n = concat $ map (toList . biasVector) (layers n)

networkHistogram :: FilePath -> (FeedForwardNetwork -> [Double]) -> FeedForwardNetwork -> IO ()
networkHistogram filename listFunction n = do
  let hist = histogram binSturges (listFunction n)
  plot filename hist
  return ()

max_stars_repo_name: AkatsukiSirius/LambdaNet

max_issues_repo_name: world-admin/LambdaNet

max_forks_repo_name: world-admin/LambdaNet


text:
-- |
-- Module    : Statistics.Matrix.Mutable
-- Copyright : (c) 2014 Bryan O'Sullivan
-- License   : BSD3
--
-- Basic mutable matrix operations.

module Statistics.Matrix.Mutable
    (
      MMatrix(..)
    , MVector
    , replicate
    , thaw
    , bounds
    , unsafeFreeze
    , unsafeRead
    , unsafeWrite
    , unsafeModify
    , immutably
    , unsafeBounds
    ) where

import Control.Applicative ((<$>))
import Control.DeepSeq (NFData(..))
import Control.Monad.ST (ST)
import Statistics.Matrix.Types (Matrix(..), MMatrix(..), MVector)
import qualified Data.Vector.Unboxed as U
import qualified Data.Vector.Unboxed.Mutable as M
import Prelude hiding (replicate)

replicate :: Int -> Int -> Double -> ST s (MMatrix s)
replicate r c k = MMatrix r c 0 <$> M.replicate (r*c) k

thaw :: Matrix -> ST s (MMatrix s)
thaw (Matrix r c e v) = MMatrix r c e <$> U.thaw v

unsafeFreeze :: MMatrix s -> ST s Matrix
unsafeFreeze (MMatrix r c e mv) = Matrix r c e <$> U.unsafeFreeze mv

unsafeRead :: MMatrix s -> Int -> Int -> ST s Double
unsafeRead mat r c = unsafeBounds mat r c M.unsafeRead
{-# INLINE unsafeRead #-}

unsafeWrite :: MMatrix s -> Int -> Int -> Double -> ST s ()
unsafeWrite mat row col k = unsafeBounds mat row col $ \v i ->
  M.unsafeWrite v i k
{-# INLINE unsafeWrite #-}

unsafeModify :: MMatrix s -> Int -> Int -> (Double -> Double) -> ST s ()
unsafeModify mat row col f = unsafeBounds mat row col $ \v i -> do
  k <- M.unsafeRead v i
  M.unsafeWrite v i (f k)
{-# INLINE unsafeModify #-}

-- | Given row and column numbers, calculate the offset into the flat
-- row-major vector.
bounds :: MMatrix s -> Int -> Int -> (MVector s -> Int -> r) -> r
bounds (MMatrix rs cs _ mv) r c k
  | r < 0 || r >= rs = error "row out of bounds"
  | c < 0 || c >= cs = error "column out of bounds"
  | otherwise        = k mv $! r * cs + c
{-# INLINE bounds #-}

-- | Given row and column numbers, calculate the offset into the flat
-- row-major vector, without checking.
unsafeBounds :: MMatrix s -> Int -> Int -> (MVector s -> Int -> r) -> r
unsafeBounds (MMatrix _ cs _ mv) r c k = k mv $! r * cs + c
{-# INLINE unsafeBounds #-}

immutably :: NFData a => MMatrix s -> (Matrix -> a) -> ST s a
immutably mmat f = do
  k <- f <$> unsafeFreeze mmat
  rnf k `seq` return k
{-# INLINE immutably #-}

max_stars_repo_name: StefanHubner/statistics

max_issues_repo_name: StefanHubner/statistics

max_forks_repo_name: StefanHubner/statistics


text:
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE TemplateHaskell #-}
{-# LANGUAGE LambdaCase #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE TypeApplications #-}

module Sketch where

import Foreign.C.Types
import Codec.Picture.ColorQuant
import Codec.Picture.Gif
import Numeric.LinearAlgebra.HMatrix as H hiding (reshape)
import Codec.Picture.Types
import Data.Int
import Foreign
import System.Environment
import Data.Time
import qualified SDL
import Prelude hiding (init)
import Graphics.Rendering.OpenGL as GL
import SDL.Vect
import Control.Monad.State as S hiding (get)
import Control.Lens
import qualified Data.Vector.Storable as SV
import Data.Map.Strict as M

import GLCode
import SDLCode

data World = World {
        _windowSize :: (Int32, Int32),
        _shaderProgram :: M.Map String GL.Program,
        _mainWindow :: SDL.Window,
        _startTime :: UTCTime
    }

$(makeLenses ''World)

newtype Options = Options {
        _shaderDirectory :: FilePath
    } deriving Show

$(makeLenses ''Options)

type SketchMonad a = StateT World IO a

lookupProgram :: String -> SketchMonad Program
lookupProgram name = do
    programs <- use shaderProgram
    let Just program = M.lookup name programs
    return program

lookupEsp :: Either String Program -> SketchMonad Program
lookupEsp (Left name) = lookupProgram name
lookupEsp (Right program) = return program

-- Uniforms

class UniformType a where
    setUniform' :: Either String Program ->
                  (Either String Program -> SketchMonad ()) -> a

instance (a ~ ()) => UniformType (SketchMonad a) where
    setUniform' program cont = cont program

setUniform :: UniformType a => String -> a
setUniform name = setUniform' (Left name) $ \_ -> return ()

-- Uniform Float
instance (UniformType r) => UniformType (String -> Float -> r) where
    setUniform' esp0 cont attr value = setUniform' esp0 $ \esp1 -> do
        program <- lookupEsp esp1
        currentProgram $= Just program
        loc <- get $ uniformLocation program attr
        uniform loc GL.$= value
        cont esp1

-- Uniform vec2f
instance (UniformType r) => UniformType (String -> GL.Vertex2 Float -> r) where
    setUniform' esp0 cont attr value = setUniform' esp0 $ \esp1 -> do
        program <- lookupEsp esp1
        currentProgram $= Just program
        loc <- get $ uniformLocation program attr
        uniform loc GL.$= value
        cont esp1

-- Uniform mat4f
instance (UniformType r) => UniformType (String -> GL.GLmatrix Float -> r) where
    setUniform' esp0 cont attr value = setUniform' esp0 $ \esp1 -> do
        program <- lookupEsp esp1
        currentProgram $= Just program
        loc <- get $ uniformLocation program attr
        uniform loc GL.$= value
        cont esp1

-- Uniform mat4f
instance (UniformType r) => UniformType (String -> H.Matrix Float -> r) where
    setUniform' esp0 cont attr value = setUniform' esp0 $ \esp1 -> do
        program <- lookupEsp esp1
        currentProgram $= Just program
        loc <- get $ uniformLocation program attr
        matr <- io $ GL.newMatrix @GL.GLmatrix @Float GL.ColumnMajor $ concat $ H.toLists $ value
        uniform loc GL.$= matr
        cont esp1

-- Some duplication here XXX
init :: FilePath -> CInt -> CInt -> CInt -> IO World
init path width height samples = do
    window <- initWindow width height samples

    start <- getCurrentTime
    let world = World {
        _windowSize = (fromIntegral width, fromIntegral height),
        _shaderProgram = M.empty,
        _mainWindow = window,
        _startTime = start
    }
    programs <- installShaders path
    return $ world & shaderProgram .~ M.fromList programs

withProgram :: String -> (GL.Program -> StateT World IO ()) -> StateT World IO ()
withProgram name cmd = do
    programs <- use shaderProgram
    forM_ (M.lookup name programs) cmd

-- Going to be user responsibility to set window size
reshape :: (Int32, Int32) -> StateT World IO ()
reshape (w, h) = do
    -- withProgram "shader" $ \program -> io $ setShaderWindow program (w, h)
    GL.viewport GL.$= (GL.Position 0 0, GL.Size (i w) (i h))
    windowSize .= (w, h)

-- Going to be user resposibility to set mouse
mouse :: (Int32, Int32) -> StateT World IO ()
mouse _ = return ()

-- mouse (x, y) = withProgram $ \program -> do
--     (_, h) <- use windowSize
--     io $ setShaderMouse program (x, h-y)

handleKey :: SDL.Keysym -> StateT World IO Bool
handleKey SDL.Keysym {SDL.keysymScancode = SDL.ScancodeEscape} = return True
handleKey SDL.Keysym { } = return False

handlePayload :: SDL.EventPayload -> StateT World IO Bool
handlePayload (SDL.WindowResizedEvent
                            SDL.WindowResizedEventData { SDL.windowResizedEventSize = V2 w h }) =
                            reshape (w, h) >> return False
handlePayload (SDL.MouseMotionEvent
                            SDL.MouseMotionEventData { SDL.mouseMotionEventPos = P (V2 x y) }) =
        mouse (x, y) >> return False
handlePayload (SDL.KeyboardEvent
                            SDL.KeyboardEventData { SDL.keyboardEventKeyMotion = SDL.Pressed, SDL.keyboardEventKeysym = k }) =
        handleKey k
handlePayload SDL.QuitEvent = return True
handlePayload _ = return False

handleUIEvent :: SDL.Event -> StateT World IO Bool
handleUIEvent SDL.Event { SDL.eventPayload = payload} = handlePayload payload

parse :: Options -> [String] -> Options
parse options [] = options
parse options ("-d" : path : args) = parse (options { _shaderDirectory = path }) args
parse _ args = error ("Incomprehensible options " ++ unwords args)

{-
mainLoop :: (Float -> StateT World IO ()) -> IO ()
mainLoop render = do
    args <- getArgs
    let options = parse (Options { _shaderDirectory = "." }) args
    print options

    SDL.initialize [SDL.InitVideo]

    world <- init (_shaderDirectory options)

    evalStateT (loop options render) world

loop :: Options -> (Float -> StateT World IO ()) -> StateT World IO ()
loop options render = do
    window <- use mainWindow
    interval <- realToFrac <$> (diffUTCTime <$> io getCurrentTime <*> use startTime)
    -- user going to set time in shader
    -- withProgram $ \program -> io $ setShaderTime program interval
    render interval
    io $ SDL.glSwapWindow window
    events <- io SDL.pollEvents
    quit <- mapM handleUIEvent events
    unless (or quit) $ loop options render
-}

makeGif :: Int32 -> Int32 -> ForeignPtr (PixelBaseComponent PixelRGB8) -> Image PixelRGB8
makeGif width32 height32 pixelData =
    let width = fromIntegral width32
        height = fromIntegral height32
        array = SV.unsafeFromForeignPtr0 pixelData (width*height*3) :: SV.Vector (PixelBaseComponent PixelRGB8)
    in Image width height array

gifLoop :: String -> Float -> Int -> Int -> Options -> [(Palette, GifDelay, Image Pixel8)] -> (Float -> SketchMonad ()) -> SketchMonad ()
gifLoop filename _ j n _ frames _ | j >= n = do
    let Right zzz = writeGifImages filename LoopingForever frames
    io zzz
gifLoop filename fps j n options frames render = do
    window <- use mainWindow
    let interval = fromIntegral j
    render (interval/fps)
    (width, height) <- use windowSize
    pixelData <- io (mallocForeignPtrArray (fromIntegral $ width*height*3) :: IO (ForeignPtr (PixelBaseComponent PixelRGB8)))
    (im, pa) <- io $ withForeignPtr pixelData $ \ptr -> do
                GL.readPixels (GL.Position 0 0) (GL.Size width height) $ GL.PixelData GL.RGB GL.UnsignedByte ptr
                let gif = makeGif width height pixelData
                let (im, pa) = palettize defaultPaletteOptions gif
                return (im, pa)
    io $ SDL.glSwapWindow window
    events <- io SDL.pollEvents
    quit <- mapM handleUIEvent events
    unless (or quit) $ gifLoop filename fps (j+1) n options ((pa, 4, im) : frames) render

initSketch :: IO Options
initSketch = do
    SDL.initialize [SDL.InitVideo]
    parse Options { _shaderDirectory = "." } <$> getArgs

mainGifLoop :: String -> Float -> CInt -> CInt -> CInt -> Int -> Int -> (Float -> SketchMonad ()) -> IO ()
mainGifLoop filename fps width height samples start end render = do
    options <- initSketch
    world <- init (_shaderDirectory options) width height samples
    io $ evalStateT (gifLoop filename fps start end options [] render) world

mainGifLoopState :: String -> Float -> CInt -> CInt -> CInt -> Int -> Int -> a -> (Float -> StateT a (StateT World IO) ()) -> IO ()
mainGifLoopState filename fps width height samples start end initial render = do
    options <- initSketch
    world <- init (_shaderDirectory options) width height samples
    evalStateT (evalStateT (gifLoopState filename fps start end options [] render) initial) world

gifLoopState :: String -> Float -> Int -> Int -> Options -> [(Palette, GifDelay, Image Pixel8)] -> (Float -> StateT a (StateT World IO) ()) -> StateT a (StateT World IO) ()
gifLoopState filename _ j n _ frames _ | j >= n = do
    case writeGifImages filename LoopingForever frames of
        Right zzz -> io zzz
        Left zzz -> do
            io $ print zzz
            io $ forM_ frames $ \(p, d, f) -> do
                print (Codec.Picture.Types.imageWidth f, Codec.Picture.Types.imageHeight f)
gifLoopState filename fps j n options frames render = do
    window <- lift $ use mainWindow
    let interval = fromIntegral j
    render (interval/fps)
    (width, height) <- lift $ use windowSize
    pixelData <- io (mallocForeignPtrArray (fromIntegral $ width*height*3) :: IO (ForeignPtr (PixelBaseComponent PixelRGB8)))
    (im, pa) <- io $ withForeignPtr pixelData $ \ptr -> do
                GL.readPixels (GL.Position 0 0) (GL.Size width height) $ GL.PixelData GL.RGB GL.UnsignedByte ptr
                let gif = makeGif width height pixelData
                let (im, pa) = palettize defaultPaletteOptions gif
                return (im, pa)
    io $ SDL.glSwapWindow window
    events <- io SDL.pollEvents
    quit <- lift $ mapM handleUIEvent events
    unless (or quit) $ gifLoopState filename fps (j+1) n options ((pa, 4, im) : frames) render

loop :: Options -> (Float -> SketchMonad ()) -> SketchMonad ()
loop options render = do
    window <- use mainWindow
    interval <- realToFrac <$> (diffUTCTime <$> io getCurrentTime <*> use startTime)
    render interval
    io $ SDL.glSwapWindow window
    events <- io SDL.pollEvents
    quit <- mapM handleUIEvent events
    unless (or quit) $ loop options render

mainLoop :: CInt -> CInt -> CInt -> (Float -> SketchMonad ()) -> IO ()
mainLoop width height samples render = do
    options <- initSketch
    world <- init (_shaderDirectory options) width height samples
    io $ evalStateT (loop options render) world

loopState :: Options -> (Float -> StateT a (StateT World IO) ()) -> StateT a (StateT World IO) ()
loopState options render = do
    window <- lift $ use mainWindow
    now <- io getCurrentTime
    start <- lift (use startTime)
    let interval = realToFrac (diffUTCTime now start)
    render interval
    io $ SDL.glSwapWindow window
    events <- io SDL.pollEvents
    quit <- lift $ mapM handleUIEvent events
    unless (or quit) $ loopState options render

mainLoopState :: CInt -> CInt -> CInt -> a -> (Float -> StateT a (StateT World IO) ()) -> IO ()
mainLoopState width height samples initial render = do
    options <- initSketch
    world <- init (_shaderDirectory options) width height samples
    evalStateT (evalStateT (loopState options render) initial) world

max_stars_repo_name: dpiponi/SketchBox

max_issues_repo_name: dpiponi/SketchBox

max_forks_repo_name: dpiponi/SketchBox


text:
-----------------------------------------------------------------------------
-- |
-- Module      :  Numeric.Signal.EEG.BDF
-- Copyright   :  (c) Alexander Vivian Hugh McPhail 2010, 2014, 2015
-- License     :  BSD3
--
-- Maintainer  :  haskell.vivian.mcphail <at> gmail <dot> com
-- Stability   :  provisional
-- Portability :
--
-- EEG BDF data file functions
--
-----------------------------------------------------------------------------

module Numeric.Signal.EEG.BDF (
                               BDF(..)
                              , loadBDFRaw
                              , loadBDF
                          ) where

import qualified Data.ByteString as BS

--import Data.Char
--import Data.List
import Data.Word
import Data.Bits
--import Data.Array.Storable
import Numeric.LinearAlgebra.Data

import qualified Numeric.Signal.Multichannel as M

import Control.Monad hiding(join)
import Control.Monad.State hiding(join)

--import System.IO

--import Debug.Trace

{- BDF Header format: http://www.biosemi.com/faq/file_format.htm
Length in bytes  	BDF Header:  	EDF Header:  	Description
8 bytes 	Byte 1: "255" (non ascii) 	Byte 1: "0" (ASCII) 	Identification code
Bytes 2-8 : "BIOSEMI" (ASCII) 	Bytes 2-8 : " "(ASCII)
80 bytes
User text input (ASCII)
	Local subject identification
80 bytes
User text input (ASCII)
	Local recording identification
8 bytes
dd.mm.yy (ASCII)
	Startdate of recording
8 bytes
hh.mm.ss (ASCII)
	Starttime of recording
8 bytes
(ASCII)
	Number of bytes in header record
44 bytes
"24BIT" (ASCII)

"BIOSEMI" (ASCII)
	Version of data format.
8 bytes
(ASCII)
	Number of data records "-1" if unknown
8 bytes
e.g.: "1" (ASCII)
	Duration of a data record, in seconds
4 bytes
e.g.: "257" or "128" (ASCII)
	Number of channels (N) in data record
N x 16 bytes
e.g.: "Fp1", "Fpz", "Fp2", etc (ASCII)
	Labels of the channels
N x 80 bytes
e.g.: "active electrode", "respiration belt" (ASCII)
	Transducer type
N x 8 bytes
e.g.: "uV", "Ohm" (ASCII)
	Physical dimension of channels
N x 8 bytes 	e.g.: "-262144" (ASCII) 	e.g.: "-32768" (ASCII) 	Physical minimum in units of physical dimension
N x 8 bytes 	e.g.: "262143" (ASCII) 	e.g.: "32767" (ASCII) 	Physical maximum in units of physical dimension
N x 8 bytes 	e.g.: "-8388608" (ASCII) 	e.g.: "-32768" (ASCII) 	Digital minimum
N x 8 bytes 	e.g.: "8388607" (ASCII) 	e.g.: "32767" (ASCII) 	Digital maximum
N x 80 bytes 	e.g.: "HP:DC; LP:410" 	e.g.: "HP:0,16; LP:500" 	Prefiltering
N x 8 bytes
For example: "2048" (ASCII)

Number of samples in each data record
(Sample-rate if Duration of data record = "1")
N x 32 bytes
(ASCII)
	Reserved
-}

type BSM a = StateT BS.ByteString IO a

getN :: Int -> BSM BS.ByteString
getN n = do
         bs <- get
         let (v,bs') = BS.splitAt n bs
         put bs'
         return v

data Date = Date { day :: Int, month :: Int, year :: Int }
          deriving(Eq,Ord,Show)

data Time = Time { hour :: Int, minute :: Int, second :: Int }
          deriving(Eq,Ord,Show)

-- | BDF datatype.  NOTE: data has already been converted to Float
data BDF = BDF {
                id_ :: !Word8
                , type_ :: !String
                , subject :: !String
                , recording :: !String
                , date :: !Date
                , time :: !Time
                , head_bytes :: !Int
                , data_version :: !String
                , num_records :: !Int
                , duration :: !Int
                , channels :: !Int
                , chan_labels :: ![String]
                , tran_type :: ![String]
                , dimensions :: ![String]
                , phys_min :: ![Int]
                , phys_max :: ![Int]
                , dig_min :: ![Int]
                , dig_max :: ![Int]
                , prefilter :: ![String]
                , samples :: ![Int]
                , reserved :: ![String]
                , data_ :: ![Vector Float]
               } --deriving(Show)

getString :: Int -> BSM String
getString n = do
              bs <- getN n
              return $ (reverse . dropWhile (== ' ') . tail . reverse . tail . show) bs

getInt :: Int -> BSM Int
getInt n = do
           s <- getString n
           return $ read s

getDate = do
          s <- getString 8
          return $ strToDate' s

strToDate' :: String -> Date
strToDate' (d1:d2:'.':m1:m2:'.':y1:y2:[]) = Date (read [d1,d2]) (read [m1,m2]) (read [y1,y2])
strToDate' _                              = error "strToDate"

getTime = do
          s <- getString 8
          return $ strToTime' s

strToTime' :: String -> Time
strToTime' (h1:h2:'.':m1:m2:'.':s1:s2:[]) = Time (read [h1,h2]) (read [m1,m2]) (read [s1,s2])
strToTime' _                              = error "strToTime"

reverseBits :: Word8 -> Word8
--reverseBits w = foldl setBit 0 $ snd . unzip . filter ((== False) . fst) $ zip (map (testBit w) [1..8]) $ reverse [1..8]
--reverseBits = foldl setBit 0 . snd . unzip . filter ((False ==) . fst) . ($ reverse [0..7]) . zip . flip map [0..7] . testBit
--reverseBits w = foldl setBit 0 . foldl ((++) . \b -> if testBit w b then [9-b] else []) [] [1..8]
reverseBits w = foldr (\b r -> if not $ testBit w b then setBit r (7-b) else r) 0 [0..7]

get24Bit :: (Int -> Float) -> BSM Float
get24Bit f = do
             b1 <- getN 1
             b2 <- getN 1
             b3 <- getN 1
--           return $ (fromIntegral $ BS.head b1) `shiftL` 16 .|. (fromIntegral $ BS.head b2) `shiftL` 8 .|. (fromIntegral $ BS.head b3)
             return $ f $ (to32 b3) `shiftL` 16 .|. (to32 b2) `shiftL` 8 .|. (to32 b1)
    where to32 = fromIntegral {- . reverseBits -} . BS.head

readRecord :: (Int -> Float) -> Int -> BSM (Vector Float)
readRecord f s = do
                 m <- replicateM s $ get24Bit f
                 return $! fromList m

readRecordBlock :: (Int -> Int -> Float) -> [(Int,Int)] -> BSM [Vector Float]
readRecordBlock f ss = mapM (\(i,s) -> readRecord (f i) s) ss

convert :: [Int] -> [Int] -> [Int] -> [Int] -> Int -> Int -> Float
convert p_min p_max d_min d_max i = \x -> (fromIntegral x) - (fromIntegral (d_min !! i))*(fromIntegral ((p_max !! i) - (p_min !! i)))/(fromIntegral ((d_max !! i) - (d_min !! i)))

readData :: (Int -> Int -> Float) -> Int -> [(Int,Int)] -> BSM [Vector Float]
readData f rs ss = do
                   lift $ putStrLn $ "channels: " ++ (show $ length ss) ++ ", records: " ++ (show rs) ++ ", samples per record: " ++ (show $ snd $ head ss)
                   d <- mapM (\x -> do
                                    lift $ putStrLn $ "Record: " ++ show x
                                    readRecordBlock f ss) [1..rs]
                   -- let v = rotate d
                   -- lift $ putStrLn $ "vectors: " ++ (show $ length v)
                   -- lift $ putStrLn $ "slices: " ++ (show $ length $ head v)
                   return $! map vjoin $! rotate_ d
    where rotate_ []            = []
          rotate_ xs@((_:[]):_) = [concat xs]
          rotate_ ((x:xs):xss)  = (x : (map head xss)) : (rotate_ (xs : (map tail xss)))

{-
readData _ 0     _  = error "readData, zeroth record requested"
readData f 1     ss = do
                      lift $ putStrLn $ "data record: 1"
                      readRecordBlock f ss
readData f (n+1) ss = do
                      lift $ putStrLn $ "data record: " ++ show (n+1)
                      bs <- readRecordBlock f ss
                      ds <- readData f n ss
                      let result = zipWith (\x y -> join [x,y]) bs ds
                      return $ result
-}

readBDF :: BSM (Maybe BDF)
readBDF = do
          id_' <- getN 1
          type_' <- getString 7
          if (not $ BS.head id_' == 255 && type_' == "BIOSEMI")
               then do
                    lift $ putStrLn "Error: File not BDF Format"
                    return Nothing
               else do
                    subject' <- getString 80
                    recording' <- getString 80
                    date' <- getDate
                    time' <- getTime
                    head_bytes' <- getInt 8
                    data_version' <- getString 44
                    num_records' <- getInt 8
                    if (num_records' == -1)
                       then do
                            lift $ putStrLn "This file is probably a valid BDF file..."
                            lift $ putStrLn "   but this program cannot grok an unspecified"
                            lift $ putStrLn "   number of records"
                            lift $ putStrLn "So complain to the software author"
                            return Nothing
                       else do
                            duration' <- getInt 8
                            channels' <- getInt 4
                            chan_labels' <- replicateM channels' $ getString 16
                            tran_type' <- replicateM channels' $ getString 80
                            dimensions' <- replicateM channels' $ getString 8
                            phys_min' <- replicateM channels' $ getInt 8
                            phys_max' <- replicateM channels' $ getInt 8
                            dig_min' <- replicateM channels' $ getInt 8
                            dig_max' <- replicateM channels' $ getInt 8
                            prefilter' <- replicateM channels' $ getString 80
                            samples' <- replicateM channels' $ getInt 8
                            reserved' <- replicateM channels' $ getString 32
                            data_' <- readData (convert phys_min' phys_max' dig_min' dig_max') num_records' (zip [0..] samples')
                            return $ Just $ BDF {
                                                 id_ = BS.head id_'
                                                 , type_ = type_'
                                                 , subject = subject'
                                                 , recording = recording'
                                                 , date = date'
                                                 , time = time'
                                                 , head_bytes = head_bytes'
                                                 , data_version = data_version'
                                                 , num_records = num_records'
                                                 , duration = duration'
                                                 , channels = channels'
                                                 , chan_labels = chan_labels'
                                                 , tran_type = tran_type'
                                                 , dimensions = dimensions'
                                                 , phys_min = phys_min'
                                                 , phys_max = phys_max'
                                                 , dig_min = dig_min'
                                                 , dig_max = dig_max'
                                                 , prefilter = prefilter'
                                                 , samples = samples'
                                                 , reserved = reserved'
                                                 , data_ = data_'
                                                }

loadBDF :: FilePath -> IO (Maybe (M.Multichannel Float))
loadBDF fn = do
             bs <- BS.readFile fn
             (bdf,bs') <- runStateT readBDF bs
             m <- case bdf of
                           (Just b) -> do
                                       return $ Just $ M.createMultichannel (head $ samples b) 24 (data_ b)
                           _        -> do
                                       putStrLn "File not read"
                                       return Nothing
             when (not (BS.null bs')) $ do
                  putStrLn "data remaining..."
             return m

-- | Load a BDF data structure
--   NOTE: data has already been converted to Float
loadBDFRaw :: FilePath -> IO (Maybe BDF)
loadBDFRaw fn = do
  bs <- BS.readFile fn
  (bdf,bs') <- runStateT readBDF bs
  m <- case bdf of
        Just b -> return (Just b)
        _      -> do
                  putStrLn "File not read"
                  return Nothing
  when (not (BS.null bs')) $ do
    putStrLn "data remaining..."
  return m

{-
getDPConversion :: BDF -> [Word32 -> Float]
getDPConversion b = let chan = channels b
                        tup4 = zip4 (phys_min b) (phys_max b) (dig_min b) (dig_max b)
                    in map (\(a,b,c,d) -> \x -> ((fromIntegral x) - (fromIntegral c))*(fromIntegral (b-a))/(fromIntegral (d-c))) tup4

writeLine :: [Float] -> String
writeLine = ((++ "\n") . concat . intersperse " " . map show)

main = do
       bs <- BS.readFile "VivianMeditation.bdf"
       (bdf,bs') <- runStateT readBDF bs
       case bdf of
                (Just b) -> do
                            let conv = getDPConversion b
                            let b' = map (zipWith ($) conv) (data_ b)
                            putStrLn $ show $ length $ head b'
                            writeFile "VivianMeditation.txt" $ concat $ map writeLine b'
                _        -> putStrLn "File not read"
       when (not (BS.null bs')) $
            do
            putStrLn "data remaining..."
            return ()
       return ()
-}

max_stars_repo_name: amcphail/hsignal

max_issues_repo_name: amcphail/hsignal

max_forks_repo_name: amcphail/hsignal


text:
module Main where

import Statistics.Quantile.Bench
import Statistics.Quantile.Exact
import Statistics.Quantile.Util
import Statistics.Quantile.Types

import System.IO

main :: IO ()
main = do
  hSetBuffering stdin LineBuffering
  hSetBuffering stdout LineBuffering
  selectFromHandle median external stdin >>= print

max_stars_repo_name: fractalcat/slides

max_issues_repo_name: fractalcat/slides

max_forks_repo_name: fractalcat/slides


text:
{-# LANGUAGE DataKinds        #-}
{-# LANGUAGE DeriveGeneric    #-}
{-# LANGUAGE DeriveGeneric    #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE GADTs            #-}
{-# LANGUAGE TemplateHaskell  #-}
{-# LANGUAGE ViewPatterns     #-}

module Main where

-- Imports for network logic

import           Control.Lens hiding                   ((<.>))
import           Numeric.Backprop
import           Numeric.LinearAlgebra.Static.Backprop

-- Imports for running network on data

import           Control.DeepSeq
import           Control.Exception
import           Control.Monad.IO.Class
import           Control.Monad.Trans.Maybe
import           Control.Monad.Trans.State
import           Data.Bitraversable
import           Data.Foldable
import           Data.IDX
-- import           Data.List
import           Data.List.Split
import           Data.Traversable
import           Data.Tuple
import           GHC.Generics                    (Generic)
import           GHC.TypeLits
import           Numeric.OneLiner
import           System.Environment
import           System.FilePath hiding          ((<.>))
import           Text.Printf
import qualified Data.Vector                     as V
import qualified Data.Vector.Generic             as VG
import qualified Numeric.LinearAlgebra           as HM
import qualified Numeric.LinearAlgebra.Static    as H
import qualified System.Random.MWC               as MWC
import qualified System.Random.MWC.Distributions as MWC

data Net = N { _weights1 :: L 250 784
             , _bias1    :: R 250
             , _weights2 :: L 10 250
             , _bias2    :: R 10
             }
  deriving (Generic)
makeLenses ''Net

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

softMax
    :: Reifies s W
    => BVar s (R 10)
    -> BVar s (R 10)
softMax x = expx / konst (norm_1V expx)
  where
    expx = exp x

crossEntropy
    :: Reifies s W
    => BVar s (R 10)
    -> BVar s (R 10)
    -> BVar s (Double)
crossEntropy targ res = -(log res <.> targ)

runNet
    :: Reifies s W
    => BVar s Net
    -> BVar s (R 784)
    -> BVar s (R 10)
runNet n x = z
  where
    y = logistic $ (n ^^. weights1) #> x + (n ^^. bias1)
    z = softMax  $ (n ^^. weights2) #> y + (n ^^. bias2)

netErr
    :: Reifies s W
    => BVar s (R 784)
    -> BVar s (R 10)
    -> BVar s Net
    -> BVar s Double
netErr x targ n = crossEntropy targ (runNet n x)

stepNet :: R 784 -> R 10 -> Net -> Net
stepNet x targ net0 = net0 - 0.02 * gr
  where
    gr :: Net
    gr = gradBP (netErr (constVar x) (constVar targ)) net0

-- *********************************************
-- Plumbing for running the network on real data
-- *********************************************

main :: IO ()
main = MWC.withSystemRandom $ \g -> do
    datadir:_ <- getArgs
    Just train <- loadMNIST (datadir </> "train-images-idx3-ubyte")
                            (datadir </> "train-labels-idx1-ubyte")
    Just test  <- loadMNIST (datadir </> "t10k-images-idx3-ubyte")
                            (datadir </> "t10k-labels-idx1-ubyte")
    putStrLn "Loaded data."
    net0 <- MWC.uniformR (-0.5, 0.5) g
    flip evalStateT net0 . forM_ [1..] $ \e -> do
      train' <- liftIO . fmap V.toList $ MWC.uniformShuffle (V.fromList train) g
      liftIO $ printf "[Epoch %d]\n" (e :: Int)

      forM_ ([1..] `zip` chunksOf 5000 train') $ \(b, chnk) -> StateT $ \n0 -> do
        printf "(Batch %d)\n" (b :: Int)

        n' <- evaluate . force $ trainList chnk n0
        printf "Trained on %d points.\n" (length chnk)

        let trainScore = testNet chnk n'
            testScore  = testNet test n'
        printf "Training error:   %.2f%%\n" ((1 - trainScore) * 100)
        printf "Validation error: %.2f%%\n" ((1 - testScore ) * 100)

        return ((), n')

trainList :: [(R 784, R 10)] -> Net -> Net
trainList = flip $ foldl' (\n (x,y) -> stepNet x y n)

testNet :: [(R 784, R 10)] -> Net -> Double
testNet xs n = sum (map (uncurry test) xs) / fromIntegral (length xs)
  where
    test x (H.extract->t)
        | HM.maxIndex t == HM.maxIndex (H.extract r) = 1
        | otherwise                                  = 0
      where
        r = evalBP (`runNet` constVar x) n

loadMNIST
    :: FilePath
    -> FilePath
    -> IO (Maybe [(R 784, R 10)])
loadMNIST fpI fpL = runMaybeT $ do
    i <- MaybeT          $ decodeIDXFile       fpI
    l <- MaybeT          $ decodeIDXLabelsFile fpL
    d <- MaybeT . return $ labeledIntData l i
    MaybeT . return $ for d (bitraverse mkImage mkLabel . swap)
  where
    mkImage = H.create . VG.convert . VG.map (\i -> fromIntegral i / 255)
    mkLabel n = H.create $ HM.build 10 (\i -> if round i == n then 1 else 0)

instance Num Net where
    (+)         = gPlus
    (-)         = gMinus
    (*)         = gTimes
    negate      = gNegate
    abs         = gAbs
    signum      = gSignum
    fromInteger = gFromInteger

instance Fractional Net where
    (/)          = gDivide
    recip        = gRecip
    fromRational = gFromRational

instance KnownNat n => MWC.Variate (R n) where
    uniform g = H.randomVector <$> MWC.uniform g <*> pure H.Uniform
    uniformR (l, h) g = (\x -> x * (h - l) + l) <$> MWC.uniform g

instance (KnownNat m, KnownNat n) => MWC.Variate (L m n) where
    uniform g = H.uniformSample <$> MWC.uniform g <*> pure 0 <*> pure 1
    uniformR (l, h) g = (\x -> x * (h - l) + l) <$> MWC.uniform g

instance MWC.Variate Net where
    uniform g = N <$> MWC.uniform g
                  <*> MWC.uniform g
                  <*> MWC.uniform g
                  <*> MWC.uniform g
    uniformR (l, h) g = (\x -> x * (h - l) + l) <$> MWC.uniform g

instance NFData Net

max_stars_repo_name: joelburget/backprop-fun

max_issues_repo_name: joelburget/backprop-fun

max_forks_repo_name: joelburget/backprop-fun


text:
module Signal where

import Data.Complex as C
import FFT
import Utils

-- Hann window function.
-- scipy.signal.hann
-- https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.signal.hann.html
hann :: (Floating a) => Int -> [a]
hann n = [0.5 - 0.5 * cos ((2 * pi * i) / (n' - 1)) | i <- map fromIntegral [0..n-1]]
  where
    n' = fromIntegral n

-- Welch Power Spectral Density (PSD) Estimation.
-- scipy.signal.welch
-- https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.signal.welch.html
--
-- xs: Timeseries data
-- windowSize: Length of each window to compute FFT.
-- overlap: Number of points to overlap between windows.
-- Returns a tuple of sample frequencies and the power spectral density of xs
--
-- Unlike the scipy implementation, returns PSD corresponding to all frequencies
-- (including those of the negative frequencies, which are simply duplicates of
-- the positive frequencies for real data).
welch :: (RealFloat a) => [a] -> Int -> Int -> ([a], [a])
welch xs windowSize overlap = (fftfreq windowSize 1.0, welchPsd)
  where
    welchPsd = map (/ norm) $ sumFrequencies $ magnitudeSquared <$> rfft <$> applyHann <$> windows
    norm = hannWindowNorm * (fromIntegral . length $ windows)
    -- sum powers over all windows correponding to the frequency bins
    sumFrequencies = foldr (zipWith (+)) (repeat 0)
    -- compute square magnitude of frequencies computed by FFT for a given window
    magnitudeSquared fs = (^2) . C.magnitude <$> fs
    -- multiply by Hann windowing function for a given data window
    applyHann window = zipWith (*) hannWindow window
    -- create windows of size exactly `windowSize` with the specified overlap
    windows = filter (\w -> length w == windowSize) [slice start windowSize xs | start <- [0,stride..n-1]]
    slice start len = drop start . take (start + len)
    n = length xs
    stride = windowSize - overlap
    -- Hann window function of length `windowSize`
    hannWindow = hann windowSize
    hannWindowNorm = sum $ (^2) <$> hannWindow

max_stars_repo_name: viswanathgs/haskell-fft

max_issues_repo_name: viswanathgs/haskell-fft

max_forks_repo_name: viswanathgs/haskell-fft


text:
-- |
-- Module      :  Acme.Operators.Base
-- License     :  BSD-3-Clause (see the file LICENSE)
-- Maintainer  :  Oleg Grenrus <oleg.grenrus@iki.fi>
--
-- Have you ever been wondering about some magic lookin operator. Here are all operators in base package. Not so many.
--
-- The operators are in few bigger logical groups, subgrouped by package. Many of them are exported from 'Prelude'.
--
-- There is a Stack Overflow QA listing prononciations of some of the operators.
-- See: <http://stackoverflow.com/questions/7746894/are-there-pronounceable-names-for-common-haskell-operators>
{-# LANGUAGE CPP #-}
module Acme.Operators.Base (
  -- * Functions
  -- ** Data.Function
  ($),
  (&),
  -- * Functors
  -- | Infamous operators on 'Functor', 'Applicative' and 'Monad'. Also 'Alternative'.

  -- ** Data.Functor
  (<$>),
  (<$),
  ($>),
  -- ** Control.Applicative
  (<*>),
  (<*),
  (*>),
  (<**>),
  (<|>),
  -- ** Control.Monad
  (>>=),
  (>>),
  (=<<),
  (>=>),
  (<=<),
  (<$!>),
  -- * Category & Arrows
  -- ** Control.Category
  -- | Dot operator '.' clashes with one in 'Prelude': 'Prelude..'. You can use '<<<'.
  (.),
  (>>>),
  (<<<),
  -- ** Control.Arrow
  -- | Optionally one can use @Arrows@ language extensions. See: <https://downloads.haskell.org/~ghc/7.10.1/docs/html/users_guide/arrow-notation.html>.
  (***),
  (&&&),
  (^>>),
  (>>^),
  (<<^),
  (^<<),
  (+++),
  (|||),
  -- * Algebraic
  -- ** Data.Eq
  (==),
  (/=),
  -- ** Data.Ord
  (<=),
  (<),
  (>),
  (>=),
  -- * Data.Monoid
  (<>),
  -- * Data
  -- ** Data.Bits
  (.&.),
  (.|.),
  -- ** Data.Bool
  (&&),
  (||),
  -- ** Data.Complex
  Complex(..), -- :+
  -- ** Data.Ratio
  (%),
  -- ** Data.List
  (++),
  (!!),
  (\\),
  -- * Numeric
  (+),
  (*),
  (-),
  (/),
  -- ** Expontentials
  (^),
  (^^),
  (**),
  -- * Type level trickery
  -- ** Data.Type.Equality
  (:~:),
  ) where

import Control.Applicative
import Control.Arrow
import Control.Category
import Control.Monad
import Data.Bits
import Data.Complex (Complex(..))
import Data.Function hiding ((.))
import Data.Functor
import Data.List
import Data.Monoid
import Data.Ratio
import Data.Type.Equality
import Prelude hiding ((.))

#if !MIN_VERSION_base(4,8,0)

infixl 4 <$!>

-- | Strict version of 'Data.Functor.<$>'.
--
-- @since 4.8.0.0
(<$!>) :: Monad m => (a -> b) -> m a -> m b
{-# INLINE (<$!>) #-}
f <$!> m = do
  x <- m
  let z = f x
  z `seq` return z

infixl 1 &

-- | '&' is a reverse application operator.  This provides notational
-- convenience.  Its precedence is one higher than that of the forward
-- application operator '$', which allows '&' to be nested in '$'.
--
-- @since 4.8.0.0
(&) :: a -> (a -> b) -> b
x & f = f x

#endif

max_stars_repo_name: phadej/acme-operators

max_issues_repo_name: phadej/acme-operators

max_forks_repo_name: phadej/acme-operators


text:
import Data.List (intercalate)
import Numeric.LinearAlgebra.Data ((|>))
import Data.Array (elems, bounds)
import Raytracer.Geometry (cube)
import Raytracer.Camera (Camera(Camera), Point(Point), calculate_rays, fire_rays)

-- This makes a crude ASCII image
make_shitty_image array = intercalate "\n" $ fmap (\list -> "." ++ list ++ ".") $ takes (x+1) $ fmap pixel $ elems array
	where
	(x, y) = let ((Point x0 y0), (Point x1 y1)) = bounds array; in (x1 - x0, y1 - y0)
	pixel (Just _) = '#'
	pixel Nothing = ' '
	takes n list = map fst $ scanl (\acc v -> splitAt v $ snd acc) (splitAt n list) $ replicate y n

main = putStrLn $ make_shitty_image $ fire_rays test_rays test_cube
	where
	wres = 160
	hres = 30
	test_camera = Camera 4 4 wres hres Nothing (3 |> [-3, -6, -3]) (3 |> [1,1,1])
	test_rays = calculate_rays test_camera
	test_cube = cube (3 |> [2, 0, 0]) (3 |> [0, 2, 0]) (3 |> [0, 0, 2]) (3 |> [0, 0, 0])


text:
module Parser
    (
    readExpr
  , parseExpr
  , module X
    ) where
import Text.ParserCombinators.Parsec hiding (spaces)
import System.Environment
import Control.Monad
import Numeric
import Data.Ratio
import Data.Complex
import Data.Array

import Number as X
import LispType as X
import Eval as X

-- | @parse p filePath input@ runs a parser @p@ over Identity without user
-- state. The @filePath@ is only used in error messages and may be the
-- empty string. Returns either a 'ParseError' ('Left')
-- or a value of type @a@ ('Right').
--
-- >  main    = case (parse numbers "" "11, 2, 43") of
-- >             Left err  -> print err
-- >             Right xs  -> print (sum xs)
-- >
-- >  numbers = commaSep integer

-- parse :: (Stream s Identity t)
--       => Parsec s () a -> SourceName -> s -> Either ParseError a
-- parse p = runP p ()

-- |
-- >>> readExpr "hoge"
-- "Found value"
readExpr :: String -> Parser LispVal -> LispVal
readExpr input parser = case parse parser "lisp" input of
  Left  err -> String $ "No match: " ++ show err
  Right val -> val

parseExpr :: Parser LispVal
parseExpr = parseAtom
        <|> parseString
        <|> parseNumber
        <|> parseBool
        <|> parseCharacter
        <|> parseQuoted
        <|> parseQuasiQuated
        <|> parseUnQuote
        <|> parseVector
        <|> parseLists

-- 引数で与えられた文字列中のどれか一文字を 認識する
--
symbol :: Parser Char
symbol = oneOf "!$%&|*+-/:<=>?@^_~"

spaces :: Parser ()
spaces = skipMany1 space

-- lispのAtom(それ以上分割できないデータ型)をパースする
-- |
-- >>> readExpr "a!fdsafdsa" parseAtom
-- "Found value"
parseAtom :: Parser LispVal
parseAtom = do
  first <- letter <|> symbol
  rest  <- many (letter <|> digit <|> symbol)
  let atom = first : rest
  return $ Atom atom

escapeChars :: Parser Char
escapeChars = do
  _ <- char '\\'
  x <- oneOf "\\\"nrt"
  return $ case x of
    'n' -> '\n'
    'r' -> '\r'
    't' -> '\t'
    _   -> x

-- |
-- >>>readExpr "\"abc\"" parseString
-- "Found value"
parseString :: Parser LispVal
parseString = do
  _ <- char '"'
  x <- many $ escapeChars <|> noneOf("\"" ++ "\\")
  _ <- char '"'
  return $ String x

-- | TDOO : description
-- >>> readExpr "#\\j" parseCharacter
-- "Found value"
-- | TDOO : description
-- >>> readExpr "#\\newline" parseCharacter
-- "Found value"
-- | TDOO : description
-- >>> readExpr "#\\space" parseCharacter
-- "Found value"
parseCharacter :: Parser LispVal
parseCharacter = do
  _     <- try $ string "#\\"
  value <- try (string "newline" <|> string "space")
          <|> do {x <- anyChar; notFollowedBy alphaNum; return [x]}
  return $ Character $ case value of
    "space"   -> ' '
    "newline" -> '\n'
    _         -> head value

-- | TDOO : description
-- >>> readExpr "#t" parseBool
-- "Found value"
-- | TDOO : description
-- >>> readExpr "#f" parseBool
-- "Found value"
parseBool :: Parser LispVal
parseBool = do
  _ <- char '#'
  (char 't' >> return (Bool True)) <|> (char 'f' >> return (Bool False))

parseLists :: Parser LispVal
parseLists = do
  _ <- char '('
  x <- try parseList <|> parseDottedList
  _ <- char ')'
  return x

-- | TDOO : description
-- >>>  readExpr  "(a test)" parseExpr
-- "Found value"
-- | TDOO : description
-- >>>  readExpr  "(a (nested) test)" parseExpr
-- "Found value"
-- | TDOO : description
-- >>>  readExpr  "(a (dotteed . list) test)" parseExpr
-- "Found value"
parseList :: Parser LispVal
parseList = liftM List $ sepBy parseExpr spaces

parseDottedList :: Parser LispVal
parseDottedList = do
  head <- endBy parseExpr spaces
  tail <- char '.' >> spaces >> parseExpr
  return $ DottedList head tail

parseQuoted :: Parser LispVal
parseQuoted = do
  _ <- char '\''
  x <- parseExpr
  return $ List [Atom "quote", x]

-- | TDOO : description
-- >>> readExpr "`(list val val)" parseQuasiQuated
-- "Found value"
parseQuasiQuated :: Parser LispVal
parseQuasiQuated = do
  _ <- char '`'
  x <- parseExpr
  return $ List [Atom "quasiquote", x]

-- | TDOO : description
-- >>> readExpr ",(list val val)" parseUnQuote
-- "Found value"
parseUnQuote :: Parser LispVal
parseUnQuote = do
  _ <- char ','
  x <- parseExpr
  return $ List [Atom "unquote", x]

-- | TDOO : description
-- >>> readExpr "#(1 2 3)" parseVector
-- "Found value"
parseVector :: Parser LispVal
parseVector = do
  _ <- string "#("
  x <- try parseVector'
  _ <- char ')'
  return x

parseVector' :: Parser LispVal
parseVector' = do
  arrayValues <- sepBy parseExpr spaces
  return $ Vector $ listArray (0, length arrayValues - 1 ) arrayValues

max_stars_repo_name: bokuo-okubo/scheme-of-bko

max_issues_repo_name: bokuo-okubo/scheme-of-bko

max_forks_repo_name: bokuo-okubo/scheme-of-bko


text:
module Classifiers.QLearning.Prelude
  ( module Prelude
  , module Control.Monad.Primitive
  , module Statistics.Distribution
  , module Statistics.Distribution.Normal
  , module System.Random.MWC
  , module Data.Sequence
  , module GHC.Float
  , module Control.Exception
  , module Data.List
  , module Data.Vector
  ) where

import Control.Monad.Primitive
import Control.Exception (assert)
import Data.List (intercalate)
import Data.Vector (Vector, (!), (//))
import Data.Sequence (Seq, (<|), (|>))
import Statistics.Distribution
import Statistics.Distribution.Normal
import System.Random.MWC
import GHC.Float

import Prelude

max_stars_repo_name: stites/rl-sandbox

max_issues_repo_name: stites/rl-sandbox

max_forks_repo_name: stites/rl-sandbox


text:
{-# LANGUAGE RecordWildCards #-}
module Fish ( Fish (..)
            , Loc
            , FishKind
            , Tolerance
            , FishData
            , move
            ) where
import Data.Complex
import qualified Data.Map.Strict as M

import Consts

type Loc = Complex Double
type FishKind = Int
type Tolerance = Double
type FishData = M.Map (Int,Int) Tolerance

data Fish = Fish {_kind :: FishKind, _loc :: Loc} deriving (Eq, Show)

distSquare :: Fish -> Fish -> Double
distSquare (Fish {_loc=l0}) (Fish {_loc=l1}) =
        sq (realPart l0 - realPart l1) + sq (imagPart l0 - imagPart l1)

sq :: Num a => a -> a
sq = (^ (2::Int))

tolerance :: FishData -> Fish -> Fish -> Double
tolerance d f0 f1 = d M.! (_kind f0, _kind f1)

tooClose :: FishData -> Fish -> Fish -> Bool
tooClose d f0 f1 = distSquare f0 f1 < tolerance d f0 f1

move :: FishData -> [Fish] -> [(Fish, Bool)]
move dat fs = map moveOne fs where
    -- TODO : consider better implementation
    force :: Fish -> Fish -> Loc
    force f0 f1 = let
        v = _loc f0 - _loc f1 in
            (dt:+0) * ((tolerance dat f0 f1 :+ 0) - abs v) * signum v
    moveOne :: Fish -> (Fish, Bool)
    moveOne f@(Fish {..}) = let
        neighbours = filter (/= f) . filter (tooClose dat f) $ fs
        mean = (/ ((fromIntegral . length) neighbours :+ 0))
        totalForce = mean . sum . map (force f) $ neighbours
        in
            (f {_loc = _loc + totalForce}, toEnum . signum $ length neighbours)

max_stars_repo_name: lesguillemets/fishpackin

max_issues_repo_name: lesguillemets/fishpackin

max_forks_repo_name: lesguillemets/fishpackin


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

module Math.Probably.MALA where

import qualified Math.Probably.PDF as PDF
import Math.Probably.Sampler
import qualified Math.Probably.PDF as PDF
import Control.Applicative

import Numeric.LinearAlgebra
import Text.Printf
import System.IO
import Data.Maybe

import Statistics.Test.KolmogorovSmirnov
import Statistics.Test.MannWhitneyU
import qualified Data.Vector.Unboxed as U
import qualified Control.Monad.State.Strict as S
import qualified Data.Vector.Storable as VS

import Debug.Trace
import Data.IORef
import Control.Spoon

data MalaPar = MalaPar { mpXi :: !(Vector Double),
                         mpPi :: !Double,
                         mpGradLast :: !(Vector Double),
                         mpSigma :: !Double,
                         mpCount :: !Double,
                         mpAccept :: !Double,
                         mpFreezeSigma :: !Bool,
                         mpThin :: !Int,
                         mpMaxGradLen :: !Double } --,
--                         mpLastRatio :: !Double }
  deriving Show

trunc t v = scale (t/(max t (norm2 v))) v

mala1 :: Covariance a => a -> (Vector Double -> (Double,Vector Double))
         -> Bool
         -> MalaPar
         -> Sampler MalaPar
mala1 cov postgrad useCache
      (MalaPar xi piCached gradientiCached sigma tr tracc freeze thinn maxGrad) = do
  let (pi,gradienti) = if useCache
                          then (piCached, gradientiCached)
                          else let (p, gr_untr) = postgrad xi
                               in (p, trunc maxGrad gr_untr)
  let xstarMean = xi + scale (sigma/2) (cov `covMul` gradienti)
  xstar <- mvnSampler xstarMean sigma cov
  u <- unitSample
  let (!pstar, gradientStarUntrunc) = postgrad xstar
      gradientStar = trunc maxGrad gradientStarUntrunc
  let !revJumpMean = xstar + scale (sigma/2) (cov `covMul` gradientStar)
      ptop = mvnPDF revJumpMean sigma cov xi
      pbot = mvnPDF xstarMean sigma cov xstar
      ratio = exp $   pstar -pi + ptop - pbot
      tr' = max 1 tr
      sigmaNext = case () of
         _ | freeze -> sigma
      accept = tracc / tr
      freezeNext = freeze {-| freeze = True
                 | not freeze = tr > 100 && accept > 0.5 && accept < 0.6  -}
  if trace (show $ (tr, pstar ,pi , ratio, sigma)) $ u < ratio
     then return $ MalaPar xstar pstar (gradientStar)
                           (if freezeNext then sigma else (min 1.4 $ 1+kmala/tr')*sigma)
                           (tr+1) (tracc+1) freezeNext  thinn maxGrad
--                           sigma (tr+1) (tracc+1)
     else return $ MalaPar xi pi ( gradienti)
                           (if freezeNext then sigma else (max 0.7143 $ 1-kmala/tr')**1.3*sigma)
                           (tr+1) tracc freezeNext thinn maxGrad
--                           sigma (tr+1) tracc

runMalaMP :: Covariance a => a -> (Vector Double -> (Double,Vector Double))
         ->  Int -> MalaPar -> [(Double,Vector Double)] -> RIO (MalaPar, [(Double,Vector Double)])
runMalaMP cov  pdf nsam init xs0 = go nsam init xs0 where
  go 0 mpar xs = do io $ putStrLn $ "MALA accept = "++show (mpAccept mpar/mpCount mpar)
                    io $ putStrLn $ "MALA sigma = "++show (mpSigma mpar)
                    return (mpar, xs)
  go n y xs = do y1 <- sample $ mala1 cov pdf True y
--                 io $ do putStrLn $ show (mpCount y1, mpPi y1, mpSigma y1)
--                         hFlush stdout
                 let newChainRes = if mpThin y1 == 0 || round (mpCount y1) `mod` mpThin y1 ==0
                                      then let !xi = mpXi y1
                                               !pi = mpPi y1
                                               !more = (pi,xi)
                                           in more:xs
                                      else xs
                 go (n-1) y1 newChainRes

runMalaUntilBetter :: Covariance a => a ->(Vector Double -> (Double,Vector Double))
         ->  Int -> MalaPar -> RIO (MalaPar, [(Double,Vector Double)])
runMalaUntilBetter cov  pdf nsam init = go nsam init [] where
  pTarget = mpPi init
  go 0 mpar xs = do io $ putStrLn $ "MALA accept = "++show (mpAccept mpar/mpCount mpar)
                    io $ putStrLn $ "MALA sigma = "++show (mpSigma mpar)
                    return (mpar, xs)
  go n y xs = do y1 <- sample $ mala1 cov  pdf True y
                 io $ do putStrLn $ show (mpCount y1, mpPi y1, mpSigma y1)
                         hFlush stdout
                 let newChainRes = if mpThin y1 == 0 || round (mpCount y1) `mod` mpThin y1 ==0
                                      then (mpPi y1, mpXi y1):xs
                                      else xs
                 if mpPi y1 > pTarget
                    then return (y1, newChainRes)
                    else go (n-1) y1 $ newChainRes

{-runMalaMPaccept' :: Matrix Double -> (Vector Double -> (Double,Vector Double))
         ->  Int -> MalaPar -> RIO (MalaPar, [(Double,Vector Double)])
runMalaMPaccept'  cov pdf nsam init  = do
  stseed <- S.get
  res <- io $ do
   seedref <- newIORef stseed
   resRef <- newIORef Nil
   let go !mpar
        | (round $ mpAccept mpar) >= nsam
          = return mpar
        | otherwise
          = do seed <- readIORef seedref
               putStr ((show (round $ mpAccept mpar)) ++".") >> hFlush stdout
               let (!mpar1, !seed1) = unSam (mala1 cov pdf mpar) seed
               writeIORef seedref seed1
               modifyIORef resRef (Cons (Pair (mpPi mpar1) ( mpXi mpar1)))
               go mpar1
   mpar <- go init
   seed <- readIORef seedref
   reslist <- readIORef resRef
   return (seed, (mpar,reslist))
  S.put $ fst res
  let mp = fst $ snd res
  return $ (mp, map unPair $ unList $ snd $ snd res)

--  pi = pdf $ toList $ mpXi init -}

runMalaMPaccept :: Covariance a => a -> (Vector Double -> (Double,Vector Double))
         -> Int -> MalaPar ->  RIO (MalaPar, [(Double,Vector Double)])
runMalaMPaccept cov pdf nsam init = go init [] where
  go !mpar !xs
    | (round $ mpAccept mpar) >= nsam = do
         io $ putStrLn $ "MALA accept = "++show (mpAccept mpar/mpCount mpar)
         io $ putStrLn $ "MALA sigma = "++show (mpSigma mpar)
         return (mpar, xs)
    | mpAccept mpar < 0.5 && mpCount mpar > 20 =
         return (mpar, [])
    | otherwise = do
         !mpar1 <- sample $ mala1 cov  pdf True mpar
         io $ putStr ((show (round $ mpAccept mpar)) ++".") >> hFlush stdout
         go mpar1 $ (mpPi mpar1, mpXi mpar1):xs

kmala = 5

{-runMalaRioESS ::  Matrix Double -> (Vector Double -> (Double,Vector Double))
                  ->  Int -> Vector Double -> RIO [Vector Double]
runMalaRioESS cov pdf want_ess xi = do
    let (p0,grad0) = pdf xi
    let sigma0 = 1.5 / (realToFrac $ dim xi) -- determined empirically
    let mp0 = MalaPar xi p0 grad0 sigma0 0 0 False
        nsam0 = want_ess*20
    io $ putStrLn $ "initial sigma = "++show (mpSigma mp0)
    (mp1, xs1) <- runMalaMP cov pdf nsam0 mp0 []
    let have_ess = min (mpAccept mp1) $ calcESSprim $ map snd xs1
    if have_ess > realToFrac want_ess
       then return $ map snd xs1
       else do let need_ess =  max 1 $ realToFrac want_ess - have_ess
                   samples_per_es = realToFrac nsam0/have_ess
                   to_do = round $ samples_per_es * need_ess
               (mp2, xs2) <- runMalaMP cov pdf to_do mp1 xs1
               return  $ map snd  xs2 -}

runMalaRioCodaESS ::  Covariance a => a-> (Vector Double -> (Double,Vector Double))
                  ->  Int -> Double -> Vector Double -> RIO [Vector Double]
runMalaRioCodaESS cov pdf want_ess truncN xi = do
    let (p0,grad0) = pdf xi
    let sigma0 = 1.5 / (realToFrac $ dim xi) -- determined empirically
    let mp0 = MalaPar xi p0 (trunc truncN grad0) sigma0 0 0 False 0 truncN
        nsam0 = want_ess*1
    let converged mp  xs = do
         let have_ess = min (mpAccept mp) $ calcESSprim $ map snd xs
         io $ putStrLn $ "ESS=" ++show have_ess
         if have_ess > realToFrac want_ess
            then return $ map snd xs
            else do let need_ess =  max 1 $ realToFrac want_ess - have_ess
                        samples_per_es = realToFrac nsam0/have_ess
                        to_do = round $ samples_per_es * need_ess
                    io $ putStrLn $ "running converged for "++show to_do
                    (mp2, xs2) <- runMalaMP cov pdf to_do (mp {mpFreezeSigma = True}) xs
                    io $ putStrLn $ "All done"
                    return $ map snd xs2
    let go mp  n xs = do
            (mp2, xs2) <- runMalaMP cov  pdf n mp []
            let testres = mannWhitneyUtest TwoTailed 0.05 (U.fromList $ map fst xs2)
                                            (U.fromList $ map fst xs)

            if testres/= Just NotSignificant
               then do io$ putStrLn $ "not converged: "++show testres++" at "++show (mpPi mp2)
                       go mp2  (round $ realToFrac n*2) xs2
               else do io$ putStrLn "converged!"
                       converged mp2  (xs2++xs)

    io $ putStrLn $ "initial sigma = "++show (mpSigma mp0)
    (mp1, xs1) <- runMalaMPaccept cov  pdf nsam0 mp0
    case xs1 of
            [] -> return []
            _ -> go mp1  (round $ mpCount mp1) xs1

{-    case (spoon (invlndet cov), mbCholSH cov) of
       (Just (covInv, (lndt,_)), Just covChol) -> go_rest (covInv, covChol)
       _ -> let cov' = PDF.posdefify cov in
            case (spoon (invlndet cov'), mbCholSH cov') of
              (Just (covInv, (lndt,_)), Just covChol) -> go_rest (covInv, covChol)
              _ ->  do io $ putStrLn "non-invertible covariance matrix"
                       return [] -}

runMalaRioSimple ::  Covariance a => a -> (Vector Double -> (Double,Vector Double))
                  ->  Int -> Double -> Int -> Vector Double -> RIO [Vector Double]
runMalaRioSimple cov pdf samples maxGrad thinN xi = do
    let (p0,grad0) = pdf xi
    let sigma0 = 2.7e-4 --1.5 / (realToFrac $ dim xi) -- determined empirically
    let mp0 = MalaPar xi p0 (trunc maxGrad grad0) sigma0 0 0 False thinN maxGrad

    (mp2, xs2) <- runMalaMP cov  pdf samples mp0 []
    io $ putStrLn $ "All done"
    return $ map snd xs2

calcCovariance :: Vector Double ->
                  Vector Double ->
                  (Vector Double -> (Double,Vector Double)) ->
                  (Vector Double -> Double) ->
                  Either (Vector Double) (Matrix Double, Matrix Double, Matrix Double)
calcCovariance vinit vnear postgrad posterior = finalcov where
   ndim = dim vinit
   finalcov
     | ndim > 0 -- > 20000 --FIXME
        = Left $ calcFDindepVars vinit vnear posterior --Left $ iCov vinit --
     | otherwise
        = hessToCov (calcFDhess vinit vnear postgrad) Nothing

iCov v = VS.replicate (VS.length v) 1 -- in (m,m,m)

calcFDindepVars v v' post = trace ("FDVars = "++show (VS.take 10 vars))  $ vars where
   hv =  mapVector (*1e-4) v --mapVector (max 1e-9 .  abs) $ v - v'
   n = dim v
   postv = post v
   postPlus i = post $ v VS.// [(i,v @>i + hv @> i)]
   postMinus i = post $ v VS.// [(i,v @>i - hv@> i)]
   vars = buildVector n fvar
   fvar i = negate $ recip $ (postPlus i - 2*postv + postMinus i)/((hv @> i)*(hv @> i))

calcFDhess v v' postgrad = hess2 where
   grad =  snd . postgrad
   gradi i = (@>i) . grad
   hv = mapVector (max 1e-9 . abs) $ v - v'
   n = dim v
   gradv = grad v
   grads =  fromRows $ flip map [0..(n-1)] $ \i ->
               grad (v VS.// [(i,v @>i + hv @> i)])
   fhess (i,j) | i<j = 0
               | otherwise = (grads @@>(j,i) - (gradv @> i))
                                       /(2*(hv @> j)) +
                             (grads @@>(i,j)- (gradv @> j))
                                       /(2*(hv @> i))

--   hess3 = scale (recip $ realToFrac n) $ sum $ map outerSelf $ grads

   hess1 = buildMatrix n n fhess
   hess2 = buildMatrix n n $ \(i,j) -> if i>=j then hess1 @@> (i,j)
                                               else hess1 @@> (j,i)

outerSelf v= v `outer` v

hessToCov hess mOriginal =
   let mTryPosdefify = case mOriginal of
            Nothing -> hessToCov (PDF.posdefify hess) (Just hess)
            Just horiginal -> let ds =  mapVector (negate . recip)
                                         $ takeDiag horiginal
                              in trace ("USING DIAGS"++show (VS.take 10 ds)) $ Left $ ds
   in
   case spoon $ inv $ negate $ hess of
     Just cov -> case mbCholSH cov of
                   Just cholm ->  trace ("invert success:"++show (VS.take 10 $ takeDiag cov) ++ "det="++show (det cov)) $ Right (cov, negate hess, cholm)
                   Nothing -> trace ("chol fail") mTryPosdefify
     Nothing -> trace ("inv fail") mTryPosdefify
                  {- Just cov -> case mbCholSH cov of
                                Just cholm ->  Right (cov, negate hessToCov,cholm)
                                Nothing -> case mbCholSH $ PDF.posdefify cov of
                                            Just cholm -> Right (cov, negate hess,cholm)
                                            Nothing -> vars
                  Nothing -> vars -}

acceptSM ampar  | mpCount ampar == 0 = "0/0"
               | otherwise = printf "%.3g" (rate::Double) ++ " ("++show yes++"/"++show total++")" where
   rate = realToFrac (yes) / realToFrac (total)
   yes = mpAccept ampar
   total = mpCount ampar

max_stars_repo_name: glutamate/probably-baysig

max_issues_repo_name: silky/probably-baysig

max_forks_repo_name: silky/probably-baysig