text:
{-|
Module      : Marvin.API.Algorithms.Internal.GradientDescent
Description : Optimization algorithm.
-}
module Marvin.API.Algorithms.Internal.GradientDescent where

import Marvin.API.Table.Internal as Table
import Marvin.API.Fallible

import Numeric.LinearAlgebra.HMatrix
import Numeric.LinearAlgebra.Data as LA

import qualified Data.Vector as Vec
import qualified Data.Vector.Unboxed as UVec

import Control.Arrow
import Control.Monad.Except
import Data.List (intercalate)

-- | Resulting model
data LinearModel = LM {
  featureColumnNames :: [Fallible ColumnName]
  , targetVariableName :: Fallible ColumnName
  , coefs :: LA.Vector Double
  , intercept_ :: Double
} deriving Eq

instance Show LinearModel where
  show model = intercalate " +\n" $ map ((++) "\t") $ show (intercept_ model) : columnStrings
    where
      columnStrings = map (\(idx, (coef, fallibleName)) ->
        let name = case fallibleName of
              Right n -> n
              Left _ -> "column #" ++ show idx
        in show coef ++ " * " ++ name) $
        zip [0..] $ zip (LA.toList (coefs model)) (featureColumnNames model)

fit' :: GradientDescent -> (NumericTable, NumericColumn) -> LinearModel
fit' gradDesc (nt, nc) = LM {
    featureColumnNames = columnNames nt
    , targetVariableName = columnName nc
    , intercept_ = intercept'
    , coefs = theta'
  }
  where
    featureMtx = tableToMatrix nt
    targetVar = columnToVector nc
    modelParams = gradientDescent gradDesc featureMtx targetVar
    theta = modelParams
    (intercept', theta') = if (addIntercept gradDesc)
      then (theta ! 0, LA.subVector 1 (LA.size theta - 1) theta)
      else (0, theta)

predict' :: LinearModel -> NumericTable -> Fallible NumericColumn
predict' model test = do
      dimCheck
      return $ unsafePredict model test
    where
      modelDim = LA.size $ coefs model
      testDim = numberOfColumns test
      dimCheck = if modelDim == testDim
          then return ()
          else throwError $
            RowLengthMismatch $
              "When attempting to make predictions based on a linear model. " ++
              "coefficients: " ++ (show modelDim) ++
              ", columns: " ++ (show testDim) ++ "."

-- | Predicts without checks.
unsafePredict :: LinearModel -> NumericTable -> NumericColumn
unsafePredict model nt = pred
  where
    pred = unsafeFromVector $ fromLinAlgVector $ batchPredict theta c x
    x = tableToMatrix nt
    theta = coefs model
    c = intercept_ model

-- | Predicts with hmatrix types.
batchPredict :: LA.Vector R -> R -> FeatureMatrix -> TargetVariable
batchPredict theta intercept x = cmap (+ intercept) (x #> theta)

-- * Type aliases
type TargetVariable   = Vector R
type FeatureVector    = Vector R

type FeatureMatrix    = Matrix R
type ModelParameters  = Vector R

type CostFunction = FeatureMatrix -> TargetVariable -> ModelParameters -> (R, Vector R)
type Regularization = R -> ModelParameters -> (R, Vector R)

-- * Parameters
data GradientDescent = GradientDescent {
  learningRate :: Double
  , addIntercept :: Bool
  , lambda :: Double
  , cost :: CostFunction
  , initModelParams :: FeatureMatrix -> TargetVariable -> ModelParameters
  , numIter :: Int
}

defaultGradientDescent = GradientDescent {
  learningRate = 0.1
  , initModelParams = \x y -> cols x |> [1,1..]
  , numIter = 10
  , lambda = 0.1
  , addIntercept = True
  , cost = linearRegressionCost
}
-- * Cost functions

logisticRegressionCost :: FeatureMatrix -> TargetVariable -> ModelParameters
                                                     -> (R, Vector R)
logisticRegressionCost x@features y@targetVariable theta@modelParameters =
  (cost, gradient)
  where
    cost            = (1 / m) * sumElements (-y * log h - (1-y) * log(1-h))
    gradient        = scale (1 / m) $ tr x #> error
    m               = fromIntegral (rows x) :: R
    error           = h - y
    h@hypothesis    = cmap sigmoid $ x #> theta

linearRegressionCost :: FeatureMatrix -> TargetVariable -> ModelParameters -> (R, Vector R)
linearRegressionCost x@features y@targetVariable theta@modelParameters =
  (cost, gradient)
  where
    cost            = (1 / (2 * m)) * sumElements (error ^ 2)
    gradient        = scale (1 / m) $ tr x #> error
    m               = fromIntegral (rows x) :: R
    error           = h - y
    h@hypothesis    = x #> theta

addRegularization :: R -> CostFunction -> Regularization -> CostFunction
addRegularization lambda costFunction regularization x y theta =
  (cost + regCost, gradient + regGradient)
  where
    (cost, gradient) = costFunction x y theta
    (regCost, regGradient) = regularization lambda theta

l2 :: R -> ModelParameters -> (R, Vector R)
l2 lambda@regularizationCoefficient theta@modelParameters =
  (cost, gradient)
  where
    cost = lambda * (norm_2 modelParameters)^2
    gradient = scale (2 * lambda) modelParameters

-- * Algorithm

gradientDescent :: GradientDescent -> FeatureMatrix -> TargetVariable ->
  ModelParameters
gradientDescent params featureMtx targetVariable =
  runIteration (numIter params) modelUpdater initModel
    where
      features = if addIntercept params
        then addInterceptFeature featureMtx
        else featureMtx
      lambda' = lambda params
      cost' = cost params
      costFunction = if lambda' /= 0
        then addRegularization lambda' cost' l2
        else cost'
      alpha                       = learningRate params
      initModel                   = initModelParams params features targetVariable
      gradient model              = snd $ costFunction features targetVariable model
      modelUpdater theta          = theta - scale alpha (gradient theta)

-- * Helpers

addInterceptFeature :: FeatureMatrix -> FeatureMatrix
addInterceptFeature features = LA.fromColumns $ constantFeature:toColumns x
  where
    x                   = features
    (m,_)               = size x
    constantFeature     = m |> repeat 1

sigmoid :: R -> R
sigmoid z = 1 / (1 + exp(-z))

runIteration :: Int -> (a -> a) -> a -> a
runIteration n f init = iterate f init !! n

-- * Conversions

columnToVector :: NumericColumn -> LA.Vector R
columnToVector = LA.fromList . UVec.toList . values

tableToMatrix :: NumericTable -> LA.Matrix R
tableToMatrix =
  columnsVec >>>
  Vec.toList >>>
  map columnToVector >>>
  LA.fromColumns

fromLinAlgVector :: LA.Vector R -> UVec.Vector R
fromLinAlgVector = LA.toList >>> UVec.fromList

toLinAlgVector :: Vec.Vector R -> LA.Vector R
toLinAlgVector = Vec.toList >>> LA.fromList

max_stars_repo_name: gaborhermann/marvin

max_issues_repo_name: gaborhermann/marvin

max_forks_repo_name: gaborhermann/marvin


text:
module Mandelbrot (
  vMandelbrot
  )where
import Data.Array
import Data.Complex

maxIterations = 100
minConvergeValue = 0.000001

vMandelbrotCurve :: Complex Double -> Complex Double
vMandelbrotCurve z = z * z - 1

vMandelbrot :: (Double,Double) -> Double -> Array (Double,Double) Integer
vMandelbrot center radius = undefined

converges :: (Complex Double -> Complex Double)
          -> Double -> Integer -> Complex Double -> Bool
converges f minVal iterations initVal =
  if magnitude (f initVal) < minVal
  then True
  else converges f minVal (iterations - 1) (f initVal)

ext: hs


text:
{-# LANGUAGE DataKinds                 #-}
{-# LANGUAGE FlexibleContexts          #-}
{-# LANGUAGE GADTs                     #-}
{-# LANGUAGE OverloadedStrings         #-}
{-# LANGUAGE ScopedTypeVariables       #-}
{-# LANGUAGE TypeApplications          #-}
{-# LANGUAGE TypeOperators             #-}
{-# LANGUAGE QuasiQuotes               #-}
{-# LANGUAGE AllowAmbiguousTypes       #-}
{-# LANGUAGE PolyKinds                 #-}
{-# LANGUAGE PartialTypeSignatures     #-}
module Main where

import           Control.Lens                    ((^.))
import qualified Control.Foldl                   as FL
import           Control.Monad.IO.Class          (MonadIO(liftIO))
import qualified Data.Map                        as M
import           Data.Maybe                      (catMaybes)
import qualified Data.Monoid                     as MO
import           Data.Proxy                      (Proxy(..))
import qualified Data.Text                       as T
import qualified Data.Text.IO                    as T
import qualified Data.Text.Lazy                  as TL
import qualified Data.Time.Calendar              as Time
import qualified Data.Vinyl                      as V
import qualified Data.Vinyl.TypeLevel            as V
import qualified Frames                          as F
import           Frames                          ((:->),(<+>),(&:))
import qualified Frames.CSV                      as F
import qualified Frames.InCore                   as F hiding (inCoreAoS)
import qualified Pipes                           as P
import qualified Pipes.Prelude                   as P
import qualified Statistics.Types               as S

import qualified Text.Blaze.Html.Renderer.Text  as BH

import qualified Frames.ParseableTypes as FP
import qualified Frames.VegaLite as FV
import qualified Frames.Transform as FT
import qualified Frames.Folds as FF
import qualified Frames.Regression as FR
import qualified Frames.MapReduce as MR
import qualified Frames.Table as Table

import qualified Frames.Visualization.VegaLite.Data        as FHV
import qualified Graphics.Vega.VegaLite.Configuration as VLC
import qualified Frames.Visualization.VegaLite.Histogram as FHV

import qualified Knit.Report                    as K
import qualified Knit.Report.Other.Blaze        as KB

import           Data.String.Here (here)

import           BlueRipple.Data.DataFrames

templateVars = M.fromList
  [
    ("lang", "English")
  , ("author", "Adam Conner-Sax")
  , ("pagetitle", "BlueRipple Adam Mission")
--  , ("tufte","True")
  ]

loadCSVToFrame :: forall rs effs. ( MonadIO (K.Sem effs)
                                  , K.LogWithPrefixesLE effs
                                  , F.ReadRec rs
                                  , F.RecVec rs
                                  , V.RMap rs)
               => F.ParserOptions -> FilePath -> (F.Record rs -> Bool) -> K.Sem effs (F.FrameRec rs)
loadCSVToFrame po fp filterF = do
  let --producer :: F.MonadSafe m => P.Producer rs m ()
      producer = F.readTableOpt po fp P.>-> P.filter filterF
  frame <- liftIO $ F.inCoreAoS producer
  let reportRows :: Foldable f => f x -> FilePath -> K.Sem effs ()
      reportRows f fn = K.logLE K.Diagnostic $ T.pack (show $ FL.fold FL.length f) <> " rows in " <> T.pack fn
  reportRows frame fp
  return frame

main :: IO ()
main = do
  let template = K.FromIncludedTemplateDir "mindoc-pandoc-KH.html"
  pandocWriterConfig <- K.mkPandocWriterConfig template
                                               templateVars
                                               K.mindocOptionsF
  let knitConfig = K.defaultKnitConfig
        { K.outerLogPrefix = Just "mission.Main"
        , K.logIf = K.logAll
        , K.pandocWriterConfig = pandocWriterConfig
        }
  eitherDocs <- K.knitHtmls knitConfig $ do
    -- load the data
    let parserOptions = F.defaultParser { F.quotingMode =  F.RFC4180Quoting ' ' }
    K.logLE K.Info "Loading data..."
    totalSpendingFrame :: F.Frame TotalSpending <- loadCSVToFrame parserOptions totalSpendingCSV (const True)
    totalSpendingBeforeFrame :: F.Frame TotalSpending <- loadCSVToFrame parserOptions totalSpendingBeforeCSV (const True)
    totalSpendingDuringFrame :: F.Frame TotalSpending <- loadCSVToFrame parserOptions totalSpendingDuringCSV (const True)
    forecastAndSpendingFrame :: F.Frame ForecastAndSpending <- loadCSVToFrame parserOptions forecastAndSpendingCSV (const True)
    electionResultsFrame :: F.Frame ElectionResults <- loadCSVToFrame parserOptions electionResultsCSV (const True)
    contextDemographicsFrame :: F.Frame ContextDemographics <- loadCSVToFrame parserOptions contextDemographicsCSV (const True)
    angryDemsFrame :: F.Frame AngryDems <- loadCSVToFrame parserOptions angryDemsCSV (const True)
--    reportRows angryDemsFrame "angryDems"
    K.logLE K.Info "Knitting..."
    K.newPandoc (K.PandocInfo "mission" M.empty) $ do
      totalSpendingHistograms totalSpendingFrame
      spendVsChangeInVoteShare totalSpendingDuringFrame totalSpendingFrame forecastAndSpendingFrame electionResultsFrame
      angryDemsAnalysis angryDemsFrame
  case eitherDocs of
    Right namedDocs ->   K.writeAllPandocResultsWithInfoAsHtml "mission/html" namedDocs --T.writeFile "mission/html/mission.html" $ TL.toStrict  $ htmlAsText
    Left err -> putStrLn $ "pandoc error: " ++ show err
--
{-
TODO
** Trends
** Registration (see Census)
* Repeat Runners
* local races?
* Ballot initiatives?
* GOTV
-}
--
angryDemsNotes :: T.Text
angryDemsNotes
  = [here|
## Angry Democrats Donations
* The [angry democrats post][AngryDemPost] generated a variety of donation size.
* As the histogram below shows, donations were mostly small with a few much larger ones.

[AngryDemPost]: <https://medium.com/@frank_s_david/angrydems-cc7e8caefe7b>
|]

angryDemsAnalysis :: K.KnitOne effs => F.Frame AngryDems -> K.Sem effs ()
angryDemsAnalysis angryDemsFrame = do
  -- aggregate by ReceiptID
  K.addMarkDown angryDemsNotes
  let byDonationFrame = FL.fold (MR.concatFold $ MR.mapReduceFold
                                 MR.noUnpack
                                 (MR.assignKeysAndData @'[ReceiptID] @'[Amount])
                                 (MR.foldAndAddKey (FF.foldAllMonoid @MO.Sum)))
                        angryDemsFrame
      vc = VLC.ViewConfig 800 400 50
  _ <- K.addHvega Nothing Nothing $
    FHV.singleHistogram @Amount "Angry Democrats Donations" (Just "# Donations") 10 VLC.Default False vc byDonationFrame
  _ <- K.addHvega Nothing Nothing $
    FHV.singleHistogram @Amount "Angry Democrats Donations (<$3000)" (Just "# Donations") 10 (VLC.GivenMinMax 0 3000) True vc byDonationFrame
  _ <- K.addHvega Nothing Nothing $
    FHV.singleHistogram @Amount "Angry Democrats Donations (<$200)" (Just "# Donations") 10 (VLC.GivenMinMax 0 200) True vc byDonationFrame
  return ()

--
spendFor r = (r ^. disbursement) + (r ^. indSupport) + realToFrac (r ^. partyExpenditures)
spendAgainst r = r ^. indOppose
proxyRace = Proxy :: Proxy '[StateAbbreviation, CongressionalDistrict]

type RaceTotalFor = "race_total_for" :-> Double
type RaceTotalAgainst = "race_total_against" :-> Double
type RaceTotalCands = "race_total_cands" :-> Int
type CandidateDiffSpend  = "candidate_diff_spend" :-> Double
type ResultVsForecast = "result_vs_forecast" :-> Double

differentialSpendNotes :: T.Text
differentialSpendNotes
  = [here|
## Differential Spending and changes in the race
* Below we look at the spending of each candidate from 8/1/2018 (the earliest date for which we have any polling/forecast data) through election day 2018.  Since we are interested in whether that spending affects the outcome of a race, we pair the spending data with an estimated change in vote share, calculated using the [538][538_house] forecast on 8/1 and the actual election result.
* Candidates, committees and parties spend money on behalf of candidates.  But Committees (Super-PACs) etc. can also spend money *against* candidates.  That makes it hard to determine how much was spent on behalf of each candidate.  We take the simplistic approach of assuming that money spent against a candidate is exactly like money spent for the opposing candidate(s).  If there is more than one opposing candidate, we divide it evenly among them.
* We divide each candidates total spending by the total spent on the entire race by all candidates, considering % of spending by each candidate rather than the actual dollars spent.  We do this so we can compare races with very different levels of spending.
* We visualize this data a couple of ways.  First, a histogram of the differential spending and then a scatter-plot of differential spending vs change in vote-share, along with a simple regression. If anything, this regression shows that extra spending has a slightly negative correlation with changes in vote-share; in line with [this][538_electionMoney] 538 story about how money affects elections.

[538_house]:<https://projects.fivethirtyeight.com/2018-midterm-election-forecast/house/>
[538_electionMoney]: <https://fivethirtyeight.com/features/money-and-elections-a-complicated-love-story/>
|]

setF :: V.KnownField t => V.Snd t -> V.ElField t
setF = V.Field

-- differential spend vs (result - 8/1 forecast)
spendVsChangeInVoteShare :: K.KnitOne effs
  => F.Frame TotalSpending -> F.Frame TotalSpending -> F.Frame ForecastAndSpending -> F.Frame ElectionResults -> K.Sem effs ()
spendVsChangeInVoteShare spendingDuringFrame totalSpendingFrame fcastAndSpendFrame eResultsFrame = K.wrapPrefix "spendVsChangeInVoteShare" $ do
  -- create a Frame with each candidate, candidate total, race total, candidate differential
  -- use aggregateFs by aggregating by [StateAbbreviation,CongressionalDistrict] and then use extract to turn those records into the ones we want
  K.logLE K.Info "Transforming spending before, forecasts and spending during, and election results..."
  let firstForecastFrame = fmap (F.rcast @[CandidateId, StateAbbreviation, CongressionalDistrict, CandidateParty, Voteshare]) $
        F.filterFrame (\r -> r ^. date == FP.FrameDay (Time.fromGregorian 2018 08 01)) fcastAndSpendFrame
  K.logLE K.Info "Filtered forecasts for first voteshare forecast on 8/1/2018 to get first forecast by candidateId"
  let proxyRace = Proxy :: Proxy '[StateAbbreviation, CongressionalDistrict]
      spendAgainst r = (r ^. indOppose)
      raceTotalsF = FF.sequenceRecFold (FF.recFieldF @RaceTotalCands FL.length id
                                        V.:& FF.recFieldF @RaceTotalFor FL.sum spendFor
                                        V.:& FF.recFieldF @RaceTotalAgainst FL.sum spendAgainst
                                        V.:& V.RNil)
      raceTotalFrameF = (MR.concatFold $ MR.mapReduceFold
                          MR.noUnpack
                          (MR.assignKeys @[StateAbbreviation, CongressionalDistrict])
                          (MR.ReduceFold $ \_ -> let f t c = fmap (V.rappend t) c in f <$> raceTotalsF <*> FL.list ))

      raceTotalFrame = F.toFrame $ fmap (F.rcast @[CandidateId,RaceTotalFor,RaceTotalAgainst,RaceTotalCands]) $ FL.fold raceTotalFrameF totalSpendingFrame

      retypeCols = FT.retypeColumn @RaceTotalFor @("race_during_for" :-> Double)
                   . FT.retypeColumn @RaceTotalAgainst @("race_during_against" :-> Double)
      raceDuringFrame = F.toFrame $ fmap (retypeCols . F.rcast @[CandidateId,RaceTotalFor,RaceTotalAgainst,Disbursement,IndSupport,IndOppose,PartyExpenditures])
                        $ FL.fold raceTotalFrameF spendingDuringFrame
  K.logLE K.Info "aggregated by race (state and district) to get total spending by race and then attached that to each candidateId"
  -- all ur joins r belong to us
  let selectResults = fmap (F.rcast @[CandidateId, FinalVoteshare])
      fRFrame = F.toFrame $ catMaybes $ fmap F.recMaybe $ F.leftJoin @'[CandidateId] firstForecastFrame (selectResults eResultsFrame)
      fRTFrame = F.toFrame $ catMaybes $ fmap F.recMaybe $ F.leftJoin @'[CandidateId] fRFrame raceTotalFrame
      fRTBFrame = F.toFrame $ catMaybes $ fmap F.recMaybe $ F.leftJoin @'[CandidateId] fRTFrame raceDuringFrame
  K.logLE K.Info $ "Did all the joins. Final frame has " <> T.pack (show $ FL.fold FL.length fRTBFrame) <> " rows."
  let addCandDiff r = FT.recordSingleton @CandidateDiffSpend $ 100*candSpend/totalSpend where
        allCandsFor = F.rgetField @RaceTotalFor r
        allCandsAgainst = F.rgetField @RaceTotalAgainst r
        allCandsAgainstDuring = F.rgetField @("race_during_against" :-> Double) r
        allCandsForDuring = F.rgetField @("race_during_for" :-> Double) r
        nCands = F.rgetField @RaceTotalCands r
        candForDuring = spendFor r
        candAgainstDuring = spendAgainst r
        candSpend = candForDuring + (allCandsAgainstDuring - realToFrac candAgainstDuring)/(realToFrac nCands - 1)
        totalSpend = allCandsForDuring + allCandsAgainstDuring
      addResultVsForecast r = FT.recordSingleton @ResultVsForecast $ (finalVS - forecastVS) where
        forecastVS = F.rgetField @Voteshare r
        finalVS = F.rgetField @FinalVoteshare r
      addAll r = addCandDiff r <+> addResultVsForecast r
      contestedRacesFrame = F.filterFrame ((==2) . F.rgetField @RaceTotalCands) $ fmap (FT.mutate addAll) fRTBFrame
  K.logLE K.Info $ "Added normalized differential spend and filtered to cheap and close races. " <> (T.pack $ show (FL.fold FL.length contestedRacesFrame)) <> " rows left."
  K.addMarkDown differentialSpendNotes
  K.addMarkDown "### All Races With 2 or more candidates"
  _ <- K.addHvega Nothing Nothing  $
    FHV.singleHistogram @CandidateDiffSpend "Distribution of Differential Spending (8/1/2018-11/6/2018)" (Just "# Candidates") 10 VLC.Default True (VLC.ViewConfig 800 400 50) contestedRacesFrame
  diffSpendVsdiffVs <- FR.ordinaryLeastSquares @_ @ResultVsForecast @True @'[CandidateDiffSpend] contestedRacesFrame
  _ <- K.addHvega Nothing Nothing $ FV.frameScatterWithFit "Differential Spending vs Change in Voteshare (8/1/208-11/6/2018)" (Just "regression") diffSpendVsdiffVs S.cl95 contestedRacesFrame
  K.addBlaze $ FR.prettyPrintRegressionResultBlaze (\y _ -> "Regression Details") diffSpendVsdiffVs S.cl95
  _ <- K.addHvega Nothing Nothing $ FV.regressionCoefficientPlot "Parameters" ["intercept","differential spend"] (FR.regressionResult diffSpendVsdiffVs) S.cl95
  let filterCloseAndCheap r =
        let vs = F.rgetField @Voteshare r
            rtf = F.rgetField @RaceTotalFor r
            rta = F.rgetField @RaceTotalAgainst r
            party = F.rgetField @CandidateParty r
        in (vs > 40) && (vs < 60) && ((rtf + rta) < 2000000) -- && (party == "Republican")
      cheapAndCloseFrame = F.filterFrame filterCloseAndCheap  contestedRacesFrame
  K.addMarkDown "### All Races With 2 or more candidates, total spending below $2,000,000, and first forecast closer than 60/40."
  _ <- K.addHvega Nothing Nothing  $
    FHV.singleHistogram @CandidateDiffSpend "Distribution of Differential Spending (8/1/2018-11/6/2018)" (Just "# Candidates") 10 VLC.Default True (VLC.ViewConfig 800 400 50) cheapAndCloseFrame
  diffSpendVsdiffVs <- FR.ordinaryLeastSquares @_ @ResultVsForecast @True @'[CandidateDiffSpend] cheapAndCloseFrame
  _ <- K.addHvega Nothing Nothing $ FV.frameScatterWithFit "Differential Spending vs Change in Voteshare (8/1/208-11/6/2018)" (Just "regression") diffSpendVsdiffVs S.cl95 cheapAndCloseFrame
  K.addBlaze $ FR.prettyPrintRegressionResultBlaze (\y _ -> "Regression Details") diffSpendVsdiffVs S.cl95
  _ <- K.addHvega Nothing Nothing $ FV.regressionCoefficientPlot "Parameters" ["intercept","differential spend"] (FR.regressionResult diffSpendVsdiffVs) S.cl95
  K.addMarkDown "### All Races With 2 or more candidates, total spending below $2,000,000, and first forecast closer than 60/40."

-- Spending histograms
spendingHistNotes :: T.Text
spendingHistNotes
  = [here|
## Spending By Party in 2018 House Races
* The bar charts below summarize the distribution of per-candidate spending in 2018 house races.  The top chart shows all the candidates, the one below zooms in on candidates who spent less than $1,000,000 and the one below that zooms in further to candidates who spent less than $100,000.
* These were produced with data from the [FEC][FECMain], using the open-data [API][FECAPI].
* This chart includes candidate expenditures, independent expenditures in support of the candidate and party expenditures on behalf of the candidate.  It ignores independent expenditures *against* the candidate, which can be considerable in some races. See the [FEC web-site for more information][FECDefs] about all of these categories.

[FECMain]: <https://www.fec.gov/>
[FECAPI]: <https://api.open.fec.gov/developers/>
[FECDefs]: <https://www.fec.gov/data/browse-data/?tab=spending>
|]

type AllSpending = "all_spending" F.:-> Double
sumSpending r =
  let db = realToFrac $ F.rgetField @Disbursement r
      is = realToFrac $ F.rgetField @IndSupport r
      pe = realToFrac $ F.rgetField @PartyExpenditures r
  in FT.recordSingleton @AllSpending (db + is + pe)

totalSpendingHistograms :: K.KnitOne effs
  => F.Frame TotalSpending -> K.Sem effs ()
totalSpendingHistograms tsFrame = do
  K.addMarkDown spendingHistNotes
  let frameWithSum = F.filterFrame ((>0). F.rgetField @AllSpending) $ fmap (FT.mutate sumSpending) tsFrame
      mergeOtherParties :: F.Record '[CandidateParty] -> F.Record '[CandidateParty]
      mergeOtherParties r =
        let p = F.rgetField @CandidateParty r
            np = case p of
              "Republican" -> "Republican"
              "Democrat" -> "Democrat"
              _ -> "_AllOthers"
        in FT.recordSingleton @CandidateParty np
      frameWithMergedOtherParties = fmap (FT.transform mergeOtherParties) frameWithSum
      vc = VLC.ViewConfig 800 400 50
--  K.log K.Diagnostic $ T.pack $ show $ fmap (show . F.rcast @[CandidateId, AllSpending]) $ FL.fold FL.list frameWithSum
  _ <- K.addHvega Nothing Nothing  $
    FHV.multiHistogram @AllSpending @CandidateParty "Distribution of Spending By Party" (Just "# Candidates") 10 (VLC.GivenMinMax 0 1e7) True FHV.AdjacentBar vc frameWithMergedOtherParties
  _ <- K.addHvega Nothing Nothing  $
    FHV.multiHistogram @AllSpending @CandidateParty "Distribution of Spending By Party (< $1,000,000)" (Just "# Candidates") 10 (VLC.GivenMinMax 0 1e6) False FHV.AdjacentBar vc frameWithMergedOtherParties
  _ <- K.addHvega Nothing Nothing  $
    FHV.multiHistogram @AllSpending @CandidateParty "Distribution of Spending By Party (< $100,000)" (Just "# Candidates") 10 (VLC.GivenMinMax 0 1e5) False FHV.AdjacentBar vc frameWithMergedOtherParties
  return ()

max_stars_repo_name: blueripple/preference-model

max_issues_repo_name: blueripple/preference-model

max_forks_repo_name: blueripple/preference-model


text:
{-# LANGUAGE TemplateHaskell #-}
{-# LANGUAGE RankNTypes #-}

module Main(main) where

import Graphics.Gloss (Display, Color, Picture(..),
    circleSolid, red, color, translate, dark,
    pictures, black, white, rotate)
import Graphics.Gloss.Data.ViewPort
import Graphics.Gloss.Interface.Pure.Game (
      Event(EventKey)
    , MouseButton(..)
    , Key(Char, MouseButton)
    , KeyState(Up, Down)
    , play
    , Display(InWindow)
    )

import qualified Data.Map.Strict as Map
import Data.List (elemIndex, mapAccumL, partition)
import Data.Tuple (swap)
import Data.Maybe (catMaybes)
import System.Random
import Control.Monad.State

import Control.Lens

{-} ###########
    # PHYSICS #
    ########### -}

import Data.Complex

type Vector = Complex Float
j :: Vector
j = 0:+1

normalize :: RealFloat a => Complex a -> Complex a
normalize (0:+0) = 0:+0
normalize a = a / (sqrt $ a * conjugate a)

toxy :: Complex a -> (a, a)
toxy (x :+ y) = (x, y)

toVector :: Float -> Vector
toVector x = x:+0

{-} ##############
    # STRUCTURES #
    ############## -}

data Ball    = Ball Vector Vector Float
data Timer   = Timer
    { _total  :: Float
    , _actual :: Float
    }
makeLenses ''Timer

data Body    = Body Ball Int

data Entity d = Entity {
      _pos    :: Vector
    , _vel    :: Vector
    , _rad    :: Float
    , _lifes  :: Int
    , _damage :: Int
    , _body :: d
    }
makeLenses ''Entity

type Bullet  = Entity ()
type Barrel  = Entity ()
type Player  = Entity ()
type Smoke   = Entity Timer

data Satelite = Satelite Float Float Float Float Float Float -- r w h rot v t

data MonsterData = MonsterData {
    _hull  :: [Satelite]
  , _shift :: Int
  }
makeLenses ''MonsterData

type Monster = Entity MonsterData

data Game = Game {
      _player   :: Player
    , _keyboard :: Map.Map Key KeyState
    , _bullets  :: [Bullet]
    , _barrels  :: [Barrel]
    , _clicks   :: [Vector]
    , _smoke    :: [Smoke]
    , _seed     :: StdGen
    , _ticks    :: Int
    , _monsters :: [Monster]
    }
makeLenses ''Game

{-} ##########
    # ENTITY #
    ########## -}

renderBody :: Entity e -> Picture
renderBody ent = translate x y $
                 color white $
                 circleSolid (ent^.rad)
                 where (x:+y) = ent^.pos

updateBody :: Float -> Entity e -> Entity e
updateBody time ent = ent & pos +~ (toVector time)*(ent^.vel)

collide :: Entity a -> Entity b -> Bool
collide a b = (realPart $ abs $ (a^.pos - b^.pos)) < (a^.rad + b^.rad)

isAlive :: Entity e -> Bool
isAlive ent = (ent^.lifes) > 0

{-} ##########
    #  GAME  #
    ########## -}

window :: Display
window = InWindow "JampaZen" (500, 500) (5, 5) -- w, h, offset

background :: Color
background = black

initialState :: Game
initialState = Game
    { _player = Entity
        { _pos = 0
        , _vel = 0
        , _rad = 10
        , _lifes = 10
        , _damage = 1
        , _body = ()
        }
    , _keyboard = Map.fromList $ zip (map Char "wasd") (repeat Up)
    , _bullets  = []
    , _barrels  = []
    , _monsters = []
    , _clicks   = []
    , _smoke    = []
    , _seed     = mkStdGen 1
    , _ticks    = 0
    }

fps :: Int
fps = 60

handleKeys :: Event -> Game -> Game
handleKeys (EventKey (Char key) keystate _ _) =
    keyboard %~ Map.adjust (const keystate) (Char key)
handleKeys (EventKey (MouseButton LeftButton) Down _ (mx, my)) =
    clicks %~ (:) (mx:+my)
handleKeys _ = id

renderGame :: Game -> Picture
renderGame game =
    let (mx:+my) = game^.player.pos
        canmap x = map renderBody (game^.x)
    in  translate (-mx) (-my) $
        pictures $ concat $ [
          [renderBody (game^.player)]
        , map renderBody (game^.bullets)
        , map renderSmoke (game^.smoke)
        , map renderMonster (game^.monsters)
        ]

updateGame :: Float -> Game -> Game
updateGame time game =
    game & ifPlayerAlive shot
         & removeSmoke
         & onceUpon 4 bulletSmoke
         & ifPlayerAlive (onceUpon 50 spawnMonster)
         & playerMove
         & collisions
         & player   %~ updateBody time
         & bullets  %~ map (updateBody time)
         & smoke    %~ map (updateSmoke time)
         & monsters %~ map (updateMonster time game)
         & ticks    +~ 1

main :: IO ()
main = play window background fps initialState renderGame handleKeys updateGame

{-} ########
    # SHOT #
    ######## -}

newBullet :: Vector -> Vector -> Bullet
newBullet p v = Entity
    { _pos=p, _vel=v, _rad=5, _lifes=1, _damage=1, _body=()}

shot :: Game -> Game
shot game = game & (clicks .~ [])
                 & (bullets %~ (++) newBullets)
      where newBullets = (newBullet (game^.player.pos) . (500*) . normalize) <$> game^.clicks

{-} ##############
    # EXPLOSIONS #
    ############## -}

randomSt :: (Float, Float) -> State StdGen Float
randomSt = state . randomR

generateExplosionSmoke :: Float -> Vector -> State StdGen Smoke
generateExplosionSmoke size pos = do
    pa <- randomSt (0, 2*pi)
    pr <- randomSt (0, 5*size)
    let p = mkPolar pr pa
    va <- randomSt (0, 2*pi)
    vr <- randomSt (10, 10*size)
    let v = mkPolar vr va
    r <- randomSt (10, 50)
    ttl <- randomSt (0.25, (sqrt size)/3)
    return Entity {_pos=p+pos, _vel=v, _rad=r,
                   _lifes=1, _damage=0, _body=newTimer ttl}

generateExplosion :: Float -> Vector -> Game -> Game
generateExplosion size pos game =
    game & seed  .~ seed'
         & smoke %~ (++) newsmoke
    where
        (newsmoke, seed') = flip runState (game^.seed) $
            sequence $ take 30 $ repeat $ generateExplosionSmoke size pos

{-} ##############
    # COLLISIONS #
    ############## -}

-- two entities collide
(<..>) :: Entity a -> Entity b -> (Entity a, Entity b)
a <..> b | collide a b && isAlive a && isAlive b = (,)
              (a & lifes %~ flip (-) (b^.damage))
              (b & lifes %~ flip (-) (a^.damage))
         | otherwise = (a, b)

(<.+>) :: Entity a -> [Entity b] -> (Entity a, [Entity b])
(<.+>) = mapAccumL (<..>)

(<+.>) :: [Entity a] -> Entity b -> ([Entity a], Entity b)
(<+.>) as b = swap $ b <.+> as

-- two list of entities collide
(<++>) :: [Entity a] -> [Entity b] -> ([Entity a], [Entity b])
(<++>) = mapAccumL (<+.>)

smash :: Lens' Game [Entity a]
      -> Lens' Game [Entity b]
      -> (Entity a -> Game -> Game)
      -> (Entity b -> Game -> Game)
      -> Game -> Game

--smash as bs destroyA destroyB b game =
smash as bs destroyA destroyB game =
    game & as .~ aliveAs
         & bs .~ aliveBs
         & foldr (.) id destructions
    where
        (as', bs') = (game^.as) <++> (game^.bs)
        (aliveAs, deadAs) = partition isAlive as'
        (aliveBs, deadBs) = partition isAlive bs'
        destructions = map destroyA deadAs ++ map destroyB deadBs

collisions :: Game -> Game
collisions =
    smash bullets monsters destroyBullet destroyMonster

{-} #########
    # TIMER #
    ######### -}

newTimer :: Float -> Timer
newTimer time = Timer {_total=time, _actual=time}

onceUpon :: Int -> (Game -> Game) -> Game -> Game
onceUpon n f game | game^.ticks `mod` n == 0 = f game
                  | otherwise                = game

{-} ##########
    # PLAYER #
    ########## -}

playerMove game = game & player.vel .~ 100*normalize (keyboardToVector $ game^.keyboard)

ifPlayerAlive :: (Game -> Game) -> Game -> Game
ifPlayerAlive f game = if isAlive (game^.player) then f game else game

stateToScalar :: KeyState -> Vector
stateToScalar Up = 1
stateToScalar Down = 0

keyToVector :: Key -> Vector
keyToVector (Char c) = maybe 0 (j^) (elemIndex c "asdw")

keyboardToVector :: Map.Map Key KeyState -> Vector
keyboardToVector = Map.foldrWithKey f (0:+0) where
    f k v a = a + (stateToScalar v)*(keyToVector k)

destroyPlayer :: Player -> Game -> Game
destroyPlayer player = generateExplosion 100 (player^.pos)

{-} ###########
    #  SMOKE  #
    ########### -}

bulletSmoke :: Game -> Game
bulletSmoke game = game & smoke %~ (++) (map bulletToSmoke $ game^.bullets)
      where bulletToSmoke bullet = Entity{
            _pos=bullet^.pos
          , _vel=0
          , _rad=bullet^.rad
          , _lifes=1
          , _damage= 0
          , _body=newTimer 0.25
          }

updateSmoke :: Float -> Smoke -> Smoke
updateSmoke time = (updateBody time) . (body.actual -~ time)

renderSmoke smoke = translate x y $
               color white $
               circleSolid r
    where (x:+y) = smoke^.pos
          r = smoke^.rad * smoke^.body.actual / smoke^.body.total

removeSmoke :: Game -> Game
removeSmoke = smoke %~ filter ((>0) . (^.body.actual))

{-} ##########
    # BULLET #
    ########## -}

destroyBullet :: Bullet -> Game -> Game
destroyBullet = const id

{-} ###########
    # MONSTER #
    ########### -}

generateSatelite :: State StdGen Satelite
generateSatelite = do
    r <- randomSt (10, 20)
    w <- randomSt (10, 50)
    h <- randomSt (10, 50)
    v <- randomSt (1, 5)
    t <- randomSt (0, 2*pi)
    rot <- randomSt (0, 2*pi)
    o <- signum <$> randomSt (-1, 1) -- orientation
    return $ Satelite r w h rot (o*v) t

generateMonster :: Vector -> Game -> Game
generateMonster p game =
    game & seed     .~ seed'
         & monsters %~ (:) Entity {
                _pos=p
              , _vel=0
              , _rad=50
              , _lifes=3
              , _damage=1
              , _body=MonsterData {
                    _hull=hull'
                  , _shift=game^.ticks
                  }
              }
    where (hull', seed') = flip runState (game^.seed) $
              sequence $ take 20 $ repeat $ generateSatelite

spawnMonster :: Game -> Game
spawnMonster game =
    game & generateMonster (mkPolar d a)
         & seed     .~ seed'
    where ((a, d), seed') = runState randomMonster (game^.seed)
          randomMonster = do
              a <- randomSt (0, 2*pi)
              d <- randomSt (1000, 2000)
              return (a, d)

renderSatelite :: Satelite -> Picture
renderSatelite (Satelite r w h rot v t) =
    rotate rot $
    translate (w*cos t) (h*sin t)$
    color white $
    circleSolid r

updateSatelite :: Float -> Satelite -> Satelite
updateSatelite time (Satelite r w h rot v t) = (Satelite r w h rot v (t + v*time))

updateMonster :: Float -> Game -> Monster -> Monster
updateMonster time game monster =
        monster & updateBody time
                & followPlayer
                & body.hull %~ map (updateSatelite time)
        where
            pp = game^.player.pos
            mp = monster^.pos
            followPlayer
              | ((game^.ticks) + (monster^.body.shift)) `mod` 120 /= 0 = id
              | otherwise = vel .~ 100*normalize (pp - mp)

renderMonster :: Monster -> Picture
renderMonster monster =
    translate x y $ pictures $ renderSatelite <$> (monster^.body.hull)
    where (x:+y) = monster^.pos

destroyMonster :: Monster -> Game -> Game
destroyMonster monster = generateExplosion 15 (monster^.pos)

max_stars_repo_name: jarys/jampazen

max_issues_repo_name: jarys/jampazen

max_forks_repo_name: jarys/jampazen


text:
module TextHistogram (
  printHistogram,
) where

import           Control.Monad (forM_, when)
import           Statistics.Test.ApproxRand
import           System.IO (hPutStrLn, stderr)
import           Text.Printf (printf)

import           Histogram

printHistogram :: Int -> TestResult -> IO ()
printHistogram bins result@(TestResult _ score _) =
  case histogram bins result of
    Left err   -> hPutStrLn stderr err
    Right hist ->
      let bucketSize = (fst $ hist !! 1) - (fst $ head hist)
          bucketHalf = bucketSize / 2
          charsPerDot = (maximum $ map snd hist) `div` 50 in
            forM_ hist $ \(label, freq) -> do
              let blocks = freq `div` charsPerDot
              when (blocks > 0) $ do
                let lower = label - bucketHalf
                let barChar = if score >= lower && score < lower + bucketSize then
                            '✣'
                          else
                            '█'
                putStr $ printf "%12.3e | " $ lower + bucketHalf
                putStrLn $ replicate blocks barChar

max_stars_repo_name: danieldk/approx-rand-test

max_issues_repo_name: danieldk/approx-rand-test

max_forks_repo_name: danieldk/approx-rand-test


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

module Data.SkyBox where

import qualified Graphics.Rendering.OpenGL.GL as GL
import           Graphics.Rendering.OpenGL.GL (($=))
import           Numeric.LinearAlgebra (flatten,ident)
import           Control.Monad
import           Control.Lens
import qualified Data.Vector.Storable as VS
import           Foreign (sizeOf)
import           Foreign.Ptr (intPtrToPtr)

import Data.Mesh
import Data.GameObject
import Data.Shaders

data SkyBox = SkyBox { _skyTex :: GL.TextureObject,
                       _skyVao :: GL.VertexArrayObject,
                       _skyVbo :: GL.BufferObject,
                       _shad   :: Shaders }

makeLenses ''SkyBox

drawSkyBox cam proj b = do
  GL.depthMask $= GL.Disabled
  GL.currentProgram $= Just (b^.shad.program)
  GL.bindVertexArrayObject $= Just (b^.skyVao)

  passMatrix "camera" $ flatten cam
  passMatrix "proj" proj

  passTexture "skytex" 3 $ Just $ b^.skyTex
  GL.drawArrays GL.Triangles 0 36
  GL.depthMask $= GL.Enabled
  where
    passMatrix = passMatrixToProg $ b^.shad
    passTexture = passTextureToProg GL.TextureCubeMap $ b^.shad

loadSkyBox = do
  let uni = ["camera", "proj", "skytex"] :: [String]
  shad <- loadShaders "shaders/cubevert.glsl" "shaders/cubefrag.glsl" uni
  tex <- loadTextures fnames
  createSkyBox shad tex

createSkyBox shad tex = do
  vao <- GL.genObjectName
  GL.bindVertexArrayObject $= Just vao
  vbo <- createAndBindArrayBuffer GL.ArrayBuffer skyboxVertices
  let glFloatSize = sizeOf (undefined :: GL.GLfloat)
  loc <- GL.get $ GL.attribLocation (shad^.program) "aPos"
  GL.vertexAttribPointer loc $=
    (GL.ToFloat,
     GL.VertexArrayDescriptor 3 GL.Float 0 $ intPtrToPtr $ fromIntegral 0)
  GL.vertexAttribArray loc $= GL.Enabled
  GL.bindVertexArrayObject $= Nothing
  return $ SkyBox tex vao vbo shad

texTargets =
  [ GL.TextureCubeMapPositiveX
  , GL.TextureCubeMapNegativeX
  , GL.TextureCubeMapPositiveY
  , GL.TextureCubeMapNegativeY
  , GL.TextureCubeMapPositiveZ
  , GL.TextureCubeMapNegativeZ ]

fnames :: [String]
fnames = -- map (("media/SkyBox/"++).(++".jpg")) . replicate 6 $ "Starscape"
  map ("media/sbox2/cgp_1/textures/"++) [ "xpos.png", "xneg.png", "ypos.png", "yneg.png", "zpos.png", "zneg.png"]

  -- [ "right"
  -- , "left"
  -- , "up"
  -- , "bottom"
  -- , "front"
  -- , "back" ]

loadTextures fns = do
  texName <- GL.genObjectName
  GL.textureBinding GL.TextureCubeMap $= Just texName
  zipWithM_ loadAndBindImage texTargets fns
  GL.textureFilter GL.TextureCubeMap $=
    ((GL.Linear', Nothing),GL.Linear')
  mapM_ setWrapping dirs
  return texName
  where setWrapping dir =
          GL.textureWrapMode GL.TextureCubeMap dir $=
            (GL.Repeated, GL.ClampToBorder)
        dirs = [GL.S, GL.T, GL.R] :: [GL.TextureCoordName]

skyboxVertices :: VS.Vector Float
skyboxVertices = [
    -1.0,  1.0, -1.0,
    -1.0, -1.0, -1.0,
     1.0, -1.0, -1.0,
     1.0, -1.0, -1.0,
     1.0,  1.0, -1.0,
    -1.0,  1.0, -1.0,

    -1.0, -1.0,  1.0,
    -1.0, -1.0, -1.0,
    -1.0,  1.0, -1.0,
    -1.0,  1.0, -1.0,
    -1.0,  1.0,  1.0,
    -1.0, -1.0,  1.0,

     1.0, -1.0, -1.0,
     1.0, -1.0,  1.0,
     1.0,  1.0,  1.0,
     1.0,  1.0,  1.0,
     1.0,  1.0, -1.0,
     1.0, -1.0, -1.0,

    -1.0, -1.0,  1.0,
    -1.0,  1.0,  1.0,
     1.0,  1.0,  1.0,
     1.0,  1.0,  1.0,
     1.0, -1.0,  1.0,
    -1.0, -1.0,  1.0,

    -1.0,  1.0, -1.0,
     1.0,  1.0, -1.0,
     1.0,  1.0,  1.0,
     1.0,  1.0,  1.0,
    -1.0,  1.0,  1.0,
    -1.0,  1.0, -1.0,

    -1.0, -1.0, -1.0,
    -1.0, -1.0,  1.0,
     1.0, -1.0, -1.0,
     1.0, -1.0, -1.0,
    -1.0, -1.0,  1.0,
     1.0, -1.0,  1.0 ]

max_stars_repo_name: Antystenes/CPG

max_issues_repo_name: Antystenes/CPG

max_forks_repo_name: Antystenes/CPG


text:
#!/usr/bin/env stack
-- stack --install-ghc runghc --resolver lts-10.4 --package hmatrix --package lens -- -Wall

{-# LANGUAGE DataKinds       #-}
{-# LANGUAGE DeriveGeneric   #-}
{-# LANGUAGE TemplateHaskell #-}

-- Imports for network logic

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

-- Imports for utility functions

import           GHC.Generics (Generic)
import           Numeric.OneLiner

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
    :: R 10
    -> R 10
softMax x = expx / konst (norm_1 expx)
  where
    expx = exp x

crossEntropy
    :: R 10
    -> R 10
    -> Double
crossEntropy targ res = -(log res <.> targ)

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

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

main :: IO ()
main = return ()

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

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

max_stars_repo_name: mstksg/inCode

max_issues_repo_name: mstksg/inCode

max_forks_repo_name: mstksg/inCode


text:
-- |
-- Module      : Foreign.CUDA.Solver.Internal.C2HS
-- Copyright   : [2017] Trevor L. McDonell
-- License     : BSD3
--
-- Maintainer  : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>
-- Stability   : experimental
-- Portability : non-portable (GHC extensions)
--

module Foreign.CUDA.Solver.Internal.C2HS (

  -- * Conversion between C and Haskell types
  cIntConv, cFloatConv, cToBool, cFromBool, cToEnum, cFromEnum,

  -- * Composite marshalling functions
  withComplex,
  peekIntConv, peekFloatConv,

) where

-- system
import Control.Monad
import Data.Complex
import Foreign
import Foreign.C
import Foreign.Storable.Complex ()

-- Conversions -----------------------------------------------------------------
--

-- | Integral conversion
--
{-# INLINE cIntConv #-}
cIntConv :: (Integral a, Integral b) => a -> b
cIntConv  = fromIntegral

-- | Floating conversion
--
{-# INLINE [1] cFloatConv #-}
cFloatConv :: (RealFloat a, RealFloat b) => a -> b
cFloatConv  = realToFrac
-- As this conversion by default goes via `Rational', it can be very slow...
{-# RULES
  "cFloatConv/Float->Float"    forall (x::Float).  cFloatConv x = x;
  "cFloatConv/Double->Double"  forall (x::Double). cFloatConv x = x;
  "cFloatConv/Float->CFloat"   forall (x::Float).  cFloatConv x = CFloat x;
  "cFloatConv/CFloat->Float"   forall (x::Float).  cFloatConv CFloat x = x;
  "cFloatConv/Double->CDouble" forall (x::Double). cFloatConv x = CDouble x;
  "cFloatConv/CDouble->Double" forall (x::Double). cFloatConv CDouble x = x
 #-}

-- | Obtain C value from Haskell 'Bool'.
--
{-# INLINE cFromBool #-}
cFromBool :: Num a => Bool -> a
cFromBool  = fromBool

-- | Obtain Haskell 'Bool' from C value.
--
{-# INLINE cToBool #-}
cToBool :: (Eq a, Num a) => a -> Bool
cToBool  = toBool

-- | Convert a C enumeration to Haskell.
--
{-# INLINE cToEnum #-}
cToEnum :: (Integral i, Enum e) => i -> e
cToEnum  = toEnum . cIntConv

-- | Convert a Haskell enumeration to C.
--
{-# INLINE cFromEnum #-}
cFromEnum :: (Enum e, Integral i) => e -> i
cFromEnum  = cIntConv . fromEnum

-- | Marshalling of complex numbers
--
{-# INLINE withComplex #-}
withComplex :: Storable a => Complex a -> (Ptr () -> IO b) -> IO b
withComplex c f = with c (f . castPtr)

-- | Marshalling of numerals
--
{-# INLINE peekIntConv #-}
peekIntConv :: (Storable a, Integral a, Integral b) => Ptr a -> IO b
peekIntConv = liftM cIntConv . peek

{-# INLINE peekFloatConv #-}
peekFloatConv :: (Storable a, RealFloat a, RealFloat b) => Ptr a -> IO b
peekFloatConv = liftM cFloatConv . peek

max_stars_repo_name: tmcdonell/cusolver

max_issues_repo_name: tmcdonell/cusolver

max_forks_repo_name: tmcdonell/cusolver


text:
module LSA
    ( compress
    ) where

import Data.List
import Data.Char
import Numeric.LinearAlgebra
--import Data.Packed.Matrix

title :: String
titles = [ "The Neatest Little Guide to Stock Market Investing"
         , "Investing For Dummies, 4th Edition"
         , "The Little Book of Common Sense Investing: The Only Way to Guarantee Your Fair Share of Stock Market Returns"
         , "The Little Book of Value Investing"
         , "Value Investing: From Graham to Buffett and Beyond"
         , "Rich Dad's Guide to Investing: What the Rich Invest in, That the Poor and the Middle Class Do Not!"
         , "Investing in Real Estate, 5th Edition"
         , "Stock Investing For Dummies"
         , "Rich Dad's Advisors: The ABC's of Real Estate Investing: The Secrets of Finding Hidden Profits Most Investors Miss"
         ]

stopwords :: String
stopwords = [ "and"
            , "edition"
            , "for"
            , "in"
            , "little"
            , "of"
            , "the"
            , "to"
            ]

-- Preps input
-- All letters to lowercase, remove all punctionation, remove stop words
docs :: [[String]]
docs = map (filter (not . (`elem` stopwords))) $
        map words $
          lines $
            filter (\x -> isAlpha x || isSpace x) $
              map toLower (unlines titles)

-- Finds word frequency over all documents
-- Removes words that appear only once, for some reason
tf :: [(String, Int)]
tf = filter (\(_,f) -> f>1) $ map (\l@(x:xs) -> (x,length l)) . group . sort $ concat docs

--Number of times term t is in document d
doc_freq :: Int -> String -> Int
doc_freq d t = length (filter (==t) (docs !! d))

--
mat :: Matrix Double
mat = buildMatrix (length tf) (length docs) ( \(term, doc) ->
          let occurances = fromIntegral $ doc_freq doc $ fst $ tf !! term -- occurance count
              docLength = length $ docs !! doc                            -- words per doc
              numDocs = length docs                                       -- number of docs
              commonness = fromIntegral $ snd $ tf !! term                -- number of docs this word occurs in
          in (occurances / docLength * log (numDocs / commonness))
      )

compress :: Int -> Matrix Double -> Matrix Double
compress k = u_k <> sigma_k <> v_k where
	(u,sigma,v) = full svd mat		                -- get SVD
	sigma_k = (takeColumns k . takeRows k) sigma	-- keep k values of Σ
	u_k = takeColumns k u				                  -- keep k columns of U
	v_k = takeRows k $ trans v			              -- keep k rows of v

max_stars_repo_name: Brightgalrs/latent-semantic-analysis

max_issues_repo_name: Brightgalrs/latent-semantic-analysis

max_forks_repo_name: Brightgalrs/latent-semantic-analysis


text:
{-# LANGUAGE FlexibleContexts #-}
-- | Pearson's chi squared test.
module Statistics.Test.ChiSquared (
    chi2test
  , chi2testCont
  , module Statistics.Test.Types
  ) where

import Prelude hiding (sum)

import Statistics.Distribution
import Statistics.Distribution.ChiSquared
import Statistics.Function        (square)
import Statistics.Sample.Internal (sum)
import Statistics.Test.Types
import Statistics.Types
import qualified Data.Vector as V
import qualified Data.Vector.Generic as G
import qualified Data.Vector.Unboxed as U

-- | Generic form of Pearson chi squared tests for binned data. Data
--   sample is supplied in form of tuples (observed quantity,
--   expected number of events). Both must be positive.
--
--   This test should be used only if all bins have expected values of
--   at least 5.
chi2test :: (G.Vector v (Int,Double), G.Vector v Double)
         => Int                 -- ^ Number of additional degrees of
                                --   freedom. One degree of freedom
                                --   is due to the fact that the are
                                --   N observation in total and
                                --   accounted for automatically.
         -> v (Int,Double)      -- ^ Observation and expectation.
         -> Maybe (Test ChiSquared)
chi2test ndf vec
  | ndf <  0  = error $ "Statistics.Test.ChiSquare.chi2test: negative NDF " ++ show ndf
  | n   > 0   = Just Test
              { testSignificance = mkPValue $ complCumulative d chi2
              , testStatistics   = chi2
              , testDistribution = chiSquared ndf
              }
  | otherwise = Nothing
  where
    n     = G.length vec - ndf - 1
    chi2  = sum $ G.map (\(o,e) -> square (fromIntegral o - e) / e) vec
    d     = chiSquared n
{-# INLINABLE  chi2test #-}
{-# SPECIALIZE
    chi2test :: Int -> U.Vector (Int,Double) -> Maybe (Test ChiSquared) #-}
{-# SPECIALIZE
    chi2test :: Int -> V.Vector (Int,Double) -> Maybe (Test ChiSquared) #-}

-- | Chi squared test for data with normal errors. Data is supplied in
--   form of pair (observation with error, and expectation).
chi2testCont
  :: (G.Vector v (Estimate NormalErr Double, Double), G.Vector v Double)
  => Int                                   -- ^ Number of additional
                                           --   degrees of freedom.
  -> v (Estimate NormalErr Double, Double) -- ^ Observation and expectation.
  -> Maybe (Test ChiSquared)
chi2testCont ndf vec
  | ndf < 0   = error $ "Statistics.Test.ChiSquare.chi2testCont: negative NDF " ++ show ndf
  | n   > 0   = Just Test
              { testSignificance = mkPValue $ complCumulative d chi2
              , testStatistics   = chi2
              , testDistribution = chiSquared ndf
              }
  | otherwise = Nothing
  where
    n     = G.length vec - ndf - 1
    chi2  = sum $ G.map (\(Estimate o (NormalErr s),e) -> square (o - e) / s) vec
    d     = chiSquared n


text:
{-# LANGUAGE FlexibleInstances #-}
module Marvin.Test.Metric where

import Marvin.API.Table.Internal
import Numeric.LinearAlgebra.HMatrix
import Numeric.LinearAlgebra.Data
import qualified Data.Vector as Vec
import qualified Data.Vector.Unboxed as UVec
import Data.Array

class Metric a where
    dist :: a -> a -> Double

instance Metric R where
    dist a b = abs $ a - b

instance Metric (Vector R) where
    dist a b = norm_Inf $ a - b

instance Metric (Matrix R) where
    dist a b = norm_Inf $ a - b

instance Metric (Vec.Vector R) where
    dist a b = Vec.foldl1 max $ Vec.zipWith (\x y -> abs (x - y)) a b

instance Metric (UVec.Vector R) where
    dist a b = UVec.foldl1 max $ UVec.zipWith (\x y -> abs (x - y)) a b

instance Metric NumericColumn where
    dist a b = dist (values a) (values b)

instance Metric NumericTable where
    dist a b = Vec.maximum $ Vec.zipWith dist (columnsVec a) (columnsVec b)

instance Metric [Double] where
    dist a b = maximum $ zipWith dist a b

instance Metric [[Double]] where
    dist a b = maximum $ zipWith dist a b

-- infinite norm
instance (Ix i, Metric a) => Metric (Array i a) where
    dist a b = maximum $ zipWith dist (elems a) (elems b)


text:
-- Copyright (c) Facebook, Inc. and its affiliates.

{-# LANGUAGE DefaultSignatures #-}
module Foreign.CPP.Addressable
  ( Addressable(..)
  , Addressable1(..)
  ) where

import Data.Complex (Complex)
import Data.Ratio (Ratio)
import Data.Functor.Const (Const)
import Data.Functor.Identity (Identity)
import Foreign hiding (sizeOf, alignment)
import Foreign.C.Types
import qualified Foreign.Storable as Storable
import GHC.Fingerprint (Fingerprint)
import System.Posix.Types

class Addressable a where
  -- | Computes the storage requirements (in bytes) of the argument.
  -- The value of the argument is not used.
  sizeOf :: a -> Int
  default sizeOf :: Storable a => a -> Int
  sizeOf = Storable.sizeOf

  -- | Computes the alignment constraint of the argument. An alignment
  -- constraint @x@ is fulfilled by any address divisible by @x@.
  -- The value of the argument is not used.
  alignment :: a -> Int
  default alignment :: Storable a => a -> Int
  alignment = Storable.alignment

class Addressable1 f where
  sizeOf1 :: f a -> Int
  alignment1 :: f a -> Int

instance Addressable Bool
instance Addressable Char
instance Addressable Double
instance Addressable Float
instance Addressable Int
instance Addressable Int8
instance Addressable Int16
instance Addressable Int32
instance Addressable Int64
instance Addressable Word
instance Addressable Word8
instance Addressable Word16
instance Addressable Word32
instance Addressable Word64
instance Addressable ()
instance Addressable Fingerprint
instance Addressable CUIntMax
instance Addressable CIntMax
instance Addressable CUIntPtr
instance Addressable CIntPtr
instance Addressable CSUSeconds
instance Addressable CUSeconds
instance Addressable CTime
instance Addressable CClock
instance Addressable CSigAtomic
instance Addressable CWchar
instance Addressable CSize
instance Addressable CPtrdiff
instance Addressable CDouble
instance Addressable CFloat
instance Addressable CULLong
instance Addressable CLLong
instance Addressable CULong
instance Addressable CLong
instance Addressable CUInt
instance Addressable CInt
instance Addressable CUShort
instance Addressable CShort
instance Addressable CUChar
instance Addressable CSChar
instance Addressable CChar
instance Addressable IntPtr
instance Addressable WordPtr
instance Addressable Fd
instance Addressable CRLim
instance Addressable CTcflag
instance Addressable CSpeed
instance Addressable CCc
instance Addressable CUid
instance Addressable CNlink
instance Addressable CGid
instance Addressable CSsize
instance Addressable CPid
instance Addressable COff
instance Addressable CMode
instance Addressable CIno
instance Addressable CDev
instance (Storable a, Integral a) => Addressable (Ratio a)
instance Addressable (StablePtr a)
instance Addressable (Ptr a)
instance Addressable (FunPtr a)
instance Storable a => Addressable (Complex a)
instance Storable a => Addressable (Identity a)
instance Storable a => Addressable (Const a b)

max_stars_repo_name: nhawkes/hsthrift

max_issues_repo_name: nhawkes/hsthrift

max_forks_repo_name: nhawkes/hsthrift


text:
module STCR2S1Edge where

import           Control.Monad           as M
import           Control.Monad.Parallel  as MP
import           Data.Array.Repa         as R
import           Data.Binary             (decodeFile)
import           Data.Complex
import           Data.List               as L
import           Data.Vector.Unboxed     as VU
import           DFT.Plan
import           Filter.Utils
import           FokkerPlanck.MonteCarlo
import           Image.Edge
import           Image.IO
import           STC.CompletionFieldR2S1
import           STC.PointSet
import           System.Directory
import           System.Environment
import           System.FilePath
import           Text.Printf
import           Types
import           Utils.Array

main = do
  args@(numPointStr:numOrientationStr:sigmaStr:taoStr:numTrailStr:maxTrailStr:initOriStr:initSpeedStr:writeFlagStr:numIterationStr:thresholdStr:histFilePath:edgeFilePath:numThreadStr:_) <-
    getArgs
  print args
  let numPoint = read numPointStr :: Int
      numOrientation = read numOrientationStr :: Int
      sigma = read sigmaStr :: Double
      tao = read taoStr :: Double
      numTrail = read numTrailStr :: Int
      maxTrail = read maxTrailStr :: Int
      initOri = read initOriStr :: Double
      initSpeed = read initSpeedStr :: Double
      writeFlag = read writeFlagStr :: Bool
      numIteration = read numIterationStr :: Int
      threshold = read thresholdStr :: Double
      numThread = read numThreadStr :: Int
      folderPath = "output/test/STCR2S1Edge"
  createDirectoryIfMissing True folderPath
  flag <- doesFileExist histFilePath
  arrG <-
    if flag
      then getNormalizedHistogramArr .
           mapHistogram
             (\x ->
                let y = x :: Int
                 in fromIntegral y) <$>
           decodeFile histFilePath
      else do
        putStrLn "Couldn't find a Green's function data. Start simulation..."
        computeS . R.map magnitude <$>
          solveMonteCarloR2S1
            numThread
            numTrail
            numPoint
            numPoint
            numOrientation
            sigma
            tao
            histFilePath
            (0, 0, 0, 0, initOri / 180 * pi, initSpeed)
  xs <- parseEdgeFile edgeFilePath
  randomPonintSet <- generateRandomPointSet 40 numPoint numPoint
  let ys =
        L.map
          (\(R2S1RPPoint (a, b, c, d)) -> (R2S1RPPoint (a - 150, b - 150, c, d)))
          xs
      zs =
        L.map
          (\(R2S1RPPoint (a, b, c, d)) ->
             (R2S1RPPoint (a - center numPoint, b - center numPoint, c, d)))
          randomPonintSet
      pointset = ys L.++ zs
      bias = computeBias numPoint numPoint numOrientation pointset
      initialEigenVec =
        computeInitialEigenVec numPoint numPoint numOrientation pointset
  powerMethod
    emptyPlan
    folderPath
    arrG
    numIteration
    writeFlag
    ""
    threshold
    bias
    initialEigenVec

max_stars_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_issues_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_forks_repo_name: XinhuaZhang/Stochastic-Completion-Field


text:
import Control.Monad.ST (runST)
import Criterion.Main
import Data.Complex
import Statistics.Sample
import Statistics.Transform
import Statistics.Correlation.Pearson
import System.Random.MWC
import qualified Data.Vector.Unboxed as U

-- Test sample
sample :: U.Vector Double
sample = runST $ flip uniformVector 10000 =<< create

-- Weighted test sample
sampleW :: U.Vector (Double,Double)
sampleW = U.zip sample (U.reverse sample)

-- Comlex vector for FFT tests
sampleC :: U.Vector (Complex Double)
sampleC = U.zipWith (:+) sample (U.reverse sample)

-- Simple benchmark for functions from Statistics.Sample
main :: IO ()
main =
  defaultMain
  [ bgroup "sample"
    [ bench "range"            $ nf (\x -> range x)            sample
      -- Mean
    , bench "mean"             $ nf (\x -> mean x)             sample
    , bench "meanWeighted"     $ nf (\x -> meanWeighted x)     sampleW
    , bench "harmonicMean"     $ nf (\x -> harmonicMean x)     sample
    , bench "geometricMean"    $ nf (\x -> geometricMean x)    sample
      -- Variance
    , bench "variance"         $ nf (\x -> variance x)         sample
    , bench "varianceUnbiased" $ nf (\x -> varianceUnbiased x) sample
    , bench "varianceWeighted" $ nf (\x -> varianceWeighted x) sampleW
      -- Correlation
    , bench "pearson"          $ nf (\x -> pearson (U.reverse sample) x) sample
    , bench "pearson'"          $ nf (\x -> pearson' (U.reverse sample) x) sample
    , bench "pearsonFast"      $ nf (\x -> pearsonFast (U.reverse sample) x) sample
      -- Other
    , bench "stdDev"           $ nf (\x -> stdDev x)           sample
    , bench "skewness"         $ nf (\x -> skewness x)         sample
    , bench "kurtosis"         $ nf (\x -> kurtosis x)         sample
      -- Central moments
    , bench "C.M. 2"           $ nf (\x -> centralMoment 2 x)  sample
    , bench "C.M. 3"           $ nf (\x -> centralMoment 3 x)  sample
    , bench "C.M. 4"           $ nf (\x -> centralMoment 4 x)  sample
    , bench "C.M. 5"           $ nf (\x -> centralMoment 5 x)  sample
    ]
  , bgroup "FFT"
    [ bgroup "fft"
      [ bench  (show n) $ whnf fft   (U.take n sampleC) | n <- fftSizes ]
    , bgroup "ifft"
      [ bench  (show n) $ whnf ifft  (U.take n sampleC) | n <- fftSizes ]
    , bgroup "dct"
      [ bench  (show n) $ whnf dct   (U.take n sample)  | n <- fftSizes ]
    , bgroup "dct_"
      [ bench  (show n) $ whnf dct_  (U.take n sampleC) | n <- fftSizes ]
    , bgroup "idct"
      [ bench  (show n) $ whnf idct  (U.take n sample)  | n <- fftSizes ]
    , bgroup "idct_"
      [ bench  (show n) $ whnf idct_ (U.take n sampleC) | n <- fftSizes ]
    ]
  ]

fftSizes :: [Int]
fftSizes = [32,128,512,2048]


text:
{-# LANGUAGE MagicHash, CPP, UnboxedTuples, BangPatterns, FlexibleContexts #-}
-----------------------------------------------------------------------------
-- |
-- Module      :  Data.Packed.Internal.Vector
-- Copyright   :  (c) Alberto Ruiz 2007
-- License     :  GPL-style
--
-- Maintainer  :  Alberto Ruiz <aruiz@um.es>
-- Stability   :  provisional
-- Portability :  portable (uses FFI)
--
-- Vector implementation
--
-----------------------------------------------------------------------------

module Data.Packed.Internal.Vector (
    Vector, dim,
    fromList, toList, (|>),
    join, (@>), safe, at, at', subVector, takesV,
    mapVector, mapVectorWithIndex, zipVectorWith, unzipVectorWith,
    mapVectorM, mapVectorM_, mapVectorWithIndexM, mapVectorWithIndexM_,
    foldVector, foldVectorG, foldLoop, foldVectorWithIndex,
    createVector, vec,
    asComplex, asReal, float2DoubleV, double2FloatV,
    stepF, stepD, condF, condD,
    conjugateQ, conjugateC,
    fwriteVector, freadVector, fprintfVector, fscanfVector,
    cloneVector,
    unsafeToForeignPtr,
    unsafeFromForeignPtr,
    unsafeWith
) where

import Data.Packed.Internal.Common
import Data.Packed.Internal.Signatures
import Foreign.Marshal.Alloc(free)
import Foreign.Marshal.Array(peekArray, pokeArray, copyArray, advancePtr)
import Foreign.ForeignPtr(ForeignPtr, castForeignPtr)
import Foreign.Ptr(Ptr)
import Foreign.Storable(Storable, peekElemOff, pokeElemOff, sizeOf)
import Foreign.C.String
import Foreign.C.Types
import Data.Complex
import Control.Monad(when)
import System.IO.Unsafe(unsafePerformIO)

#if __GLASGOW_HASKELL__ >= 605
import GHC.ForeignPtr           (mallocPlainForeignPtrBytes)
#else
import Foreign.ForeignPtr       (mallocForeignPtrBytes)
#endif

import GHC.Base
#if __GLASGOW_HASKELL__ < 612
import GHC.IOBase hiding (liftIO)
#endif

import qualified Data.Vector.Storable as Vector
import Data.Vector.Storable(Vector,
                            unsafeToForeignPtr,
                            unsafeFromForeignPtr,
                            unsafeWith)

-- | Number of elements
dim :: (Storable t) => Vector t -> Int
dim = Vector.length

-- C-Haskell vector adapter
-- vec :: Adapt (CInt -> Ptr t -> r) (Vector t) r
vec :: (Storable t) => Vector t -> (((CInt -> Ptr t -> t1) -> t1) -> IO b) -> IO b
vec x f = unsafeWith x $ \p -> do
    let v g = do
        g (fi $ dim x) p
    f v
{-# INLINE vec #-}

-- allocates memory for a new vector
createVector :: Storable a => Int -> IO (Vector a)
createVector n = do
    when (n <= 0) $ error ("trying to createVector of dim "++show n)
    fp <- doMalloc undefined
    return $ unsafeFromForeignPtr fp 0 n
  where
    --
    -- Use the much cheaper Haskell heap allocated storage
    -- for foreign pointer space we control
    --
    doMalloc :: Storable b => b -> IO (ForeignPtr b)
    doMalloc dummy = do
#if __GLASGOW_HASKELL__ >= 605
        mallocPlainForeignPtrBytes (n * sizeOf dummy)
#else
        mallocForeignPtrBytes      (n * sizeOf dummy)
#endif

{- | creates a Vector from a list:

@> fromList [2,3,5,7]
4 |> [2.0,3.0,5.0,7.0]@

-}
fromList :: Storable a => [a] -> Vector a
fromList l = unsafePerformIO $ do
    v <- createVector (length l)
    unsafeWith v $ \ p -> pokeArray p l
    return v

safeRead v = inlinePerformIO . unsafeWith v
{-# INLINE safeRead #-}

inlinePerformIO :: IO a -> a
inlinePerformIO (IO m) = case m realWorld# of (# _, r #) -> r
{-# INLINE inlinePerformIO #-}

{- | extracts the Vector elements to a list

@> toList (linspace 5 (1,10))
[1.0,3.25,5.5,7.75,10.0]@

-}
toList :: Storable a => Vector a -> [a]
toList v = safeRead v $ peekArray (dim v)

{- | An alternative to 'fromList' with explicit dimension. The input
     list is explicitly truncated if it is too long, so it may safely
     be used, for instance, with infinite lists.

     This is the format used in the instances for Show (Vector a).
-}
(|>) :: (Storable a) => Int -> [a] -> Vector a
infixl 9 |>
n |> l = if length l' == n
            then fromList l'
            else error "list too short for |>"
  where l' = take n l

-- | access to Vector elements without range checking
at' :: Storable a => Vector a -> Int -> a
at' v n = safeRead v $ flip peekElemOff n
{-# INLINE at' #-}

--
-- turn off bounds checking with -funsafe at configure time.
-- ghc will optimise away the salways true case at compile time.
--
#if defined(UNSAFE)
safe :: Bool
safe = False
#else
safe = True
#endif

-- | access to Vector elements with range checking.
at :: Storable a => Vector a -> Int -> a
at v n
    | safe      = if n >= 0 && n < dim v
                    then at' v n
                    else error "vector index out of range"
    | otherwise = at' v n
{-# INLINE at #-}

{- | takes a number of consecutive elements from a Vector

@> subVector 2 3 (fromList [1..10])
3 |> [3.0,4.0,5.0]@

-}
subVector :: Storable t => Int       -- ^ index of the starting element
                        -> Int       -- ^ number of elements to extract
                        -> Vector t  -- ^ source
                        -> Vector t  -- ^ result
subVector = Vector.slice

{- | Reads a vector position:

@> fromList [0..9] \@\> 7
7.0@

-}
(@>) :: Storable t => Vector t -> Int -> t
infixl 9 @>
(@>) = at

{- | creates a new Vector by joining a list of Vectors

@> join [fromList [1..5], constant 1 3]
8 |> [1.0,2.0,3.0,4.0,5.0,1.0,1.0,1.0]@

-}
join :: Storable t => [Vector t] -> Vector t
join [] = error "joining zero vectors"
join [v] = v
join as = unsafePerformIO $ do
    let tot = sum (map dim as)
    r <- createVector tot
    unsafeWith r $ \ptr ->
        joiner as tot ptr
    return r
  where joiner [] _ _ = return ()
        joiner (v:cs) _ p = do
            let n = dim v
            unsafeWith v $ \pb -> copyArray p pb n
            joiner cs 0 (advancePtr p n)

{- | Extract consecutive subvectors of the given sizes.

@> takesV [3,4] (linspace 10 (1,10))
[3 |> [1.0,2.0,3.0],4 |> [4.0,5.0,6.0,7.0]]@

-}
takesV :: Storable t => [Int] -> Vector t -> [Vector t]
takesV ms w | sum ms > dim w = error $ "takesV " ++ show ms ++ " on dim = " ++ (show $ dim w)
            | otherwise = go ms w
    where go [] _ = []
          go (n:ns) v = subVector 0 n v
                      : go ns (subVector n (dim v - n) v)

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

-- | transforms a complex vector into a real vector with alternating real and imaginary parts
asReal :: (RealFloat a, Storable a) => Vector (Complex a) -> Vector a
asReal v = unsafeFromForeignPtr (castForeignPtr fp) (2*i) (2*n)
    where (fp,i,n) = unsafeToForeignPtr v

-- | transforms a real vector into a complex vector with alternating real and imaginary parts
asComplex :: (RealFloat a, Storable a) => Vector a -> Vector (Complex a)
asComplex v = unsafeFromForeignPtr (castForeignPtr fp) (i `div` 2) (n `div` 2)
    where (fp,i,n) = unsafeToForeignPtr v

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

float2DoubleV :: Vector Float -> Vector Double
float2DoubleV v = unsafePerformIO $ do
    r <- createVector (dim v)
    app2 c_float2double vec v vec r "float2double"
    return r

double2FloatV :: Vector Double -> Vector Float
double2FloatV v = unsafePerformIO $ do
    r <- createVector (dim v)
    app2 c_double2float vec v vec r "double2float2"
    return r

foreign import ccall unsafe "float2double" c_float2double:: TFV
foreign import ccall unsafe "double2float" c_double2float:: TVF

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

stepF :: Vector Float -> Vector Float
stepF v = unsafePerformIO $ do
    r <- createVector (dim v)
    app2 c_stepF vec v vec r "stepF"
    return r

stepD :: Vector Double -> Vector Double
stepD v = unsafePerformIO $ do
    r <- createVector (dim v)
    app2 c_stepD vec v vec r "stepD"
    return r

foreign import ccall unsafe "stepF" c_stepF :: TFF
foreign import ccall unsafe "stepD" c_stepD :: TVV

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

condF :: Vector Float -> Vector Float -> Vector Float -> Vector Float -> Vector Float -> Vector Float
condF x y l e g = unsafePerformIO $ do
    r <- createVector (dim x)
    app6 c_condF vec x vec y vec l vec e vec g vec r "condF"
    return r

condD :: Vector Double -> Vector Double -> Vector Double -> Vector Double -> Vector Double -> Vector Double
condD x y l e g = unsafePerformIO $ do
    r <- createVector (dim x)
    app6 c_condD vec x vec y vec l vec e vec g vec r "condD"
    return r

foreign import ccall unsafe "condF" c_condF :: CInt -> PF -> CInt -> PF -> CInt -> PF -> TFFF
foreign import ccall unsafe "condD" c_condD :: CInt -> PD -> CInt -> PD -> CInt -> PD -> TVVV

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

conjugateAux fun x = unsafePerformIO $ do
    v <- createVector (dim x)
    app2 fun vec x vec v "conjugateAux"
    return v

conjugateQ :: Vector (Complex Float) -> Vector (Complex Float)
conjugateQ = conjugateAux c_conjugateQ
foreign import ccall unsafe "conjugateQ" c_conjugateQ :: TQVQV

conjugateC :: Vector (Complex Double) -> Vector (Complex Double)
conjugateC = conjugateAux c_conjugateC
foreign import ccall unsafe "conjugateC" c_conjugateC :: TCVCV

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

cloneVector :: Storable t => Vector t -> IO (Vector t)
cloneVector v = do
        let n = dim v
        r <- createVector n
        let f _ s _ d =  copyArray d s n >> return 0
        app2 f vec v vec r "cloneVector"
        return r

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

-- | map on Vectors
mapVector :: (Storable a, Storable b) => (a-> b) -> Vector a -> Vector b
mapVector f v = unsafePerformIO $ do
    w <- createVector (dim v)
    unsafeWith v $ \p ->
        unsafeWith w $ \q -> do
            let go (-1) = return ()
                go !k = do x <- peekElemOff p k
                           pokeElemOff      q k (f x)
                           go (k-1)
            go (dim v -1)
    return w
{-# INLINE mapVector #-}

-- | zipWith for Vectors
zipVectorWith :: (Storable a, Storable b, Storable c) => (a-> b -> c) -> Vector a -> Vector b -> Vector c
zipVectorWith f u v = unsafePerformIO $ do
    let n = min (dim u) (dim v)
    w <- createVector n
    unsafeWith u $ \pu ->
        unsafeWith v $ \pv ->
            unsafeWith w $ \pw -> do
                let go (-1) = return ()
                    go !k = do x <- peekElemOff pu k
                               y <- peekElemOff pv k
                               pokeElemOff      pw k (f x y)
                               go (k-1)
                go (n -1)
    return w
{-# INLINE zipVectorWith #-}

-- | unzipWith for Vectors
unzipVectorWith :: (Storable (a,b), Storable c, Storable d)
                   => ((a,b) -> (c,d)) -> Vector (a,b) -> (Vector c,Vector d)
unzipVectorWith f u = unsafePerformIO $ do
      let n = dim u
      v <- createVector n
      w <- createVector n
      unsafeWith u $ \pu ->
          unsafeWith v $ \pv ->
              unsafeWith w $ \pw -> do
                  let go (-1) = return ()
                      go !k   = do z <- peekElemOff pu k
                                   let (x,y) = f z
                                   pokeElemOff      pv k x
                                   pokeElemOff      pw k y
                                   go (k-1)
                  go (n-1)
      return (v,w)
{-# INLINE unzipVectorWith #-}

foldVector :: Storable a => (a -> b -> b) -> b -> Vector a -> b
foldVector f x v = unsafePerformIO $
    unsafeWith v $ \p -> do
        let go (-1) s = return s
            go !k !s = do y <- peekElemOff p k
                          go (k-1::Int) (f y s)
        go (dim v -1) x
{-# INLINE foldVector #-}

-- the zero-indexed index is passed to the folding function
foldVectorWithIndex :: Storable a => (Int -> a -> b -> b) -> b -> Vector a -> b
foldVectorWithIndex f x v = unsafePerformIO $
    unsafeWith v $ \p -> do
        let go (-1) s = return s
            go !k !s = do y <- peekElemOff p k
                          go (k-1::Int) (f k y s)
        go (dim v -1) x
{-# INLINE foldVectorWithIndex #-}

foldLoop f s0 d = go (d - 1) s0
     where
       go 0 s = f (0::Int) s
       go !j !s = go (j - 1) (f j s)

foldVectorG f s0 v = foldLoop g s0 (dim v)
    where g !k !s = f k (at' v) s
          {-# INLINE g #-} -- Thanks to Ryan Ingram (http://permalink.gmane.org/gmane.comp.lang.haskell.cafe/46479)
{-# INLINE foldVectorG #-}

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

-- | monadic map over Vectors
--    the monad @m@ must be strict
mapVectorM :: (Storable a, Storable b, Monad m) => (a -> m b) -> Vector a -> m (Vector b)
mapVectorM f v = do
    w <- return $! unsafePerformIO $! createVector (dim v)
    mapVectorM' w 0 (dim v -1)
    return w
    where mapVectorM' w' !k !t
              | k == t               = do
                                       x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
                                       y <- f x
                                       return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y
              | otherwise            = do
                                       x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
                                       y <- f x
                                       _ <- return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y
                                       mapVectorM' w' (k+1) t
{-# INLINE mapVectorM #-}

-- | monadic map over Vectors
mapVectorM_ :: (Storable a, Monad m) => (a -> m ()) -> Vector a -> m ()
mapVectorM_ f v = do
    mapVectorM' 0 (dim v -1)
    where mapVectorM' !k !t
              | k == t            = do
                                    x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
                                    f x
              | otherwise         = do
                                    x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
                                    _ <- f x
                                    mapVectorM' (k+1) t
{-# INLINE mapVectorM_ #-}

-- | monadic map over Vectors with the zero-indexed index passed to the mapping function
--    the monad @m@ must be strict
mapVectorWithIndexM :: (Storable a, Storable b, Monad m) => (Int -> a -> m b) -> Vector a -> m (Vector b)
mapVectorWithIndexM f v = do
    w <- return $! unsafePerformIO $! createVector (dim v)
    mapVectorM' w 0 (dim v -1)
    return w
    where mapVectorM' w' !k !t
              | k == t               = do
                                       x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
                                       y <- f k x
                                       return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y
              | otherwise            = do
                                       x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
                                       y <- f k x
                                       _ <- return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y
                                       mapVectorM' w' (k+1) t
{-# INLINE mapVectorWithIndexM #-}

-- | monadic map over Vectors with the zero-indexed index passed to the mapping function
mapVectorWithIndexM_ :: (Storable a, Monad m) => (Int -> a -> m ()) -> Vector a -> m ()
mapVectorWithIndexM_ f v = do
    mapVectorM' 0 (dim v -1)
    where mapVectorM' !k !t
              | k == t            = do
                                    x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
                                    f k x
              | otherwise         = do
                                    x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
                                    _ <- f k x
                                    mapVectorM' (k+1) t
{-# INLINE mapVectorWithIndexM_ #-}

mapVectorWithIndex :: (Storable a, Storable b) => (Int -> a -> b) -> Vector a -> Vector b
--mapVectorWithIndex g = head . mapVectorWithIndexM (\a b -> [g a b])
mapVectorWithIndex f v = unsafePerformIO $ do
    w <- createVector (dim v)
    unsafeWith v $ \p ->
        unsafeWith w $ \q -> do
            let go (-1) = return ()
                go !k = do x <- peekElemOff p k
                           pokeElemOff      q k (f k x)
                           go (k-1)
            go (dim v -1)
    return w
{-# INLINE mapVectorWithIndex #-}

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

-- | Loads a vector from an ASCII file (the number of elements must be known in advance).
fscanfVector :: FilePath -> Int -> IO (Vector Double)
fscanfVector filename n = do
    charname <- newCString filename
    res <- createVector n
    app1 (gsl_vector_fscanf charname) vec res "gsl_vector_fscanf"
    free charname
    return res

foreign import ccall unsafe "vector_fscanf" gsl_vector_fscanf:: Ptr CChar -> TV

-- | Saves the elements of a vector, with a given format (%f, %e, %g), to an ASCII file.
fprintfVector :: FilePath -> String -> Vector Double -> IO ()
fprintfVector filename fmt v = do
    charname <- newCString filename
    charfmt <- newCString fmt
    app1 (gsl_vector_fprintf charname charfmt) vec v "gsl_vector_fprintf"
    free charname
    free charfmt

foreign import ccall unsafe "vector_fprintf" gsl_vector_fprintf :: Ptr CChar -> Ptr CChar -> TV

-- | Loads a vector from a binary file (the number of elements must be known in advance).
freadVector :: FilePath -> Int -> IO (Vector Double)
freadVector filename n = do
    charname <- newCString filename
    res <- createVector n
    app1 (gsl_vector_fread charname) vec res "gsl_vector_fread"
    free charname
    return res

foreign import ccall unsafe "vector_fread" gsl_vector_fread:: Ptr CChar -> TV

-- | Saves the elements of a vector to a binary file.
fwriteVector :: FilePath -> Vector Double -> IO ()
fwriteVector filename v = do
    charname <- newCString filename
    app1 (gsl_vector_fwrite charname) vec v "gsl_vector_fwrite"
    free charname

foreign import ccall unsafe "vector_fwrite" gsl_vector_fwrite :: Ptr CChar -> TV

max_stars_repo_name: curiousleo/liquidhaskell

max_issues_repo_name: curiousleo/liquidhaskell

max_forks_repo_name: curiousleo/liquidhaskell


text:
{-# LANGUAGE OverloadedStrings #-}
import Test.HUnit
import Halite

import Data.Text (Text(..))
import Data.Attoparsec.Text

import Numeric.LinearAlgebra (vector)

main :: IO ()
main = runTestTT allTests >>= print

allTests = TestList [ parserTests, geometryTests ]

exampleInit = "\
   \1\n\
   \240 160\n\
   \ 2 0 3 0 60.0000 80.0000 255 0.0000 0.0000 0 0 0 0 1 60.0000 77.0000 255 0.0000 0.0000 0 0 0 0 2 60.0000 83.0000 255 0.0000 0.0000 0 0 0 0 1 3 3 180.0000 80.0000 255 0.0000 0.0000 0 0 0 0 4 180.0000 77.0000 255 0.0000 0.0000 0 0 0 0 5 180.0000 83.0000 255 0.0000 0.0000 0 0 0 0 12 0 130.7715 90.7715 1814 7.1166 3 0 1024 0 0 0 1 109.2285 90.7715 1814 7.1166 3 0 1024 0 0 0 2 109.2285 69.2285 1814 7.1166 3 0 1024 0 0 0 3 130.7715 69.2285 1814 7.1166 3 0 1024 0 0 0 4 85.6058 18.5976 3224 12.6469 5 0 1821 0 0 0 5 165.4537 33.5377 3224 12.6469 5 0 1821 0 0 0 6 154.3942 141.4024 3224 12.6469 5 0 1821 0 0 0 7 74.5463 126.4623 3224 12.6469 5 0 1821 0 0 0 8 134.7591 28.6182 1150 4.5129 2 0 649 0 0 0 9 214.8587 87.9945 1150 4.5129 2 0 649 0 0 0 10 105.2409 131.3818 1150 4.5129 2 0 649 0 0 0 11 25.1413 72.0055 1150 4.5129 2 0 649 0 0 0\n"

examplePlayer = "0 3 0 60.0000 80.0000 255 0.0000 0.0000 0 0 0 0 1 60.0000 77.0000 255 0.0000 0.0000 0 0 0 0 2 60.0000 83.0000 255 0.0000 0.0000 0 0 0 0\n"

assertDoneWith :: Show r => (r -> IO a) -> Parser r -> Text -> Test
assertDoneWith f p t = TestCase $ do
  case parse p t of
    Fail _ contexts err -> assertFailure err
    Partial f -> assertFailure $ "failed: partial parse: " ++ show (f "")
    Done _ r -> putStr " succeeded with result: " >> f r >> print r

assertDone :: Show r => Parser r -> Text -> Test
assertDone = assertDoneWith (const $ return ())

parserTests = TestLabel "Parser tests" $ TestList
  [ TestLabel "parse testUid" $ assertDone (lineOf parseUid) "0\n"

  , TestLabel "parse sepByN decimal" $
      assertDoneWith (assertEqual "result" [1,2,3,4,5])
        (lineOf $ sepByN 5 skipSpace decimal) "1 2 3 4 5\n"

  , TestLabel "parse sepByN empty" $
      assertDoneWith (assertEqual "result" [])
        (lineOf $ sepByN 0 skipSpace decimal) "\n"

  , TestLabel "parse list decimal" $
      assertDoneWith (assertEqual "result" (5,[1,2,3,4,5]))
        (lineOf $ list decimal) "5 1 2 3 4 5\n"

  , TestLabel "parse list empty" $
      assertDoneWith (assertEqual "result" (0, []))
        (lineOf $ list decimal) "0\n"

  , TestLabel "parse header" $ assertDone parseHeader "0\n240 160\n"

  , TestLabel "parse docking info Undocked" $
      assertDoneWith (assertEqual "pdi" Undocked) (lineOf parseDockingInfo) "0 0 0\n"

  , TestLabel "parse docking info Docking" $
      assertDoneWith (assertEqual "pdi" (Docking (Id 314) 0)) (lineOf parseDockingInfo) "1 314 0\n"

  , TestLabel "parse pos" $
      assertDoneWith (assertEqual "pos" $ vector [60.0, 80.0]) (lineOf parsePos) "60.0000 80.0000\n"

  , TestLabel "parse ship" $
      assertDone (lineOf ship) "0 60.0000 80.0000 255 0.0000 0.0000 0 0 0 0\n"

  , TestLabel "parse player" $
      assertDone (lineOf player) examplePlayer

  , TestLabel "parse planet" $
      assertDone (lineOf planet) "11 25.1413 72.0055 1150 4.5129 2 0 649 0 0 0\n"

  , TestLabel "parse init"    $ assertDone parseInit exampleInit
  ]

testLineIntersectsCircle expected x0 x1 c r =
  assertEqual "intersects" expected . maybe False (const True) $
    lineIntersectsCircle (Line x0 x1) (Circle c r)

geometryTests = TestLabel "Geometry tests" $ TestList
  [ TestLabel "lineIntersectsCircle 1" . TestCase $
      testLineIntersectsCircle True  (vector [0,0]) (vector [2,2]) (vector [1,1]) 0.5

  , TestLabel "lineIntersectsCircle 2" . TestCase $
      testLineIntersectsCircle False (vector [0,0]) (vector [1,0]) (vector [10,10]) 0.5

  -- Test unit circle intersects line y = 0 at (1,0) and (-1,0)
  , TestLabel "lineIntersectsCircle 3" . TestCase $
      assertEqual "unity" (Just (vector [1.0,0.0], vector[-1.0,0.0])) $ lineIntersectsCircle (Line (vector [0,0]) (vector [1,0]) ) (Circle (vector [0,0]) 1)
  ]


text:
module Main where

import           BuckinghamPi
import           FileOperations
import           Numeric.LinearAlgebra.HMatrix

rho = DimensionalVariable "rho" [(Mass   , 1.0)  , (Length , -3.0)]
vel = DimensionalVariable "v"   [(Length , 1.0)  , (Time   , -1.0)]
l   = DimensionalVariable "l"   [(Length , 1.0)]
mu  = DimensionalVariable "mu"  [(Mass   , 1.0)  , (Length , -1.0)  , (Time , -1.0)]
a   = DimensionalVariable "a"   [(Length , 1.0)  , (Time   , -1.0)]

dimVars = [rho, vel, l, mu, a]
funUnits = [Mass, Length, Time]

--writeDimensionlessGroups :: [DimensionlessVariable] -> FilePath -> IO ()
main :: IO ()
main = fmap (\d -> generatePiGroups d [0..(numFundamentalUnits dimVars - 1)]) (readVariablesFromFile "./app/dimensional_variables") >>= (\v -> writeDimensionlessGroups v "dimensionless_groups")


text:
{-# LANGUAGE KindSignatures, DataKinds #-}
{-# LANGUAGE FlexibleContexts, FlexibleInstances, MultiParamTypeClasses #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE TypeFamilies #-}

module DSProg (
  module DSProg,
  Result
) where

import Prelude hiding ((<>))

import Control.Monad (join, (<=<), forM_, void, replicateM)
import Control.Monad.State (StateT, MonadState (get, put), evalStateT)
import qualified Control.Monad.State as State
import Control.Monad.Bayes.Class (MonadSample, MonadInfer)
import Control.Monad.Bayes.Sampler
import Control.Monad.Bayes.Weighted
import Control.Monad.Trans (liftIO)
import Control.Monad.Fix

import Data.Aeson (ToJSON (toJSON))
import Data.Maybe (isNothing, fromJust)
import qualified Data.Map
import Data.Proxy (Proxy (Proxy))

import GHC.TypeLits

import Numeric.LinearAlgebra.Static hiding (Gaussian, M)

import Unsafe.Coerce (unsafeCoerce)

import qualified Util.AssocList as M
import Util.AssocList (Map)
import qualified Util.AssocListSet as S
import Util.AssocListSet (Set)
import qualified Util.MStream as MS
import Util.MStream (MStream)
import Distributions
import MVDistributions
import DelayedSampling
import Inference (particles, Weight, zunheap)
import Util.Numeric (average)
import Util.Ref

import qualified Util.ZStream as ZS

import Numeric.Log (Log (Exp))

instance ToJSON a => ToJSON (Result a) where
  toJSON (RConst x) = toJSON x
  toJSON (RMarginal d) = toJSON d

instance MonadState s m => MonadState s (Weighted m) where
  get = State.lift get
  put = State.lift . put
  state = State.lift . State.state

meanMDistr :: Double ~ MType a => MDistr a -> Double
meanMDistr (MGaussian mu var) = mu
meanMDistr (MBeta a b) = a / (a + b)

meanResult :: Result Double -> Double
meanResult (RConst x) = x
meanResult (RMarginal m) = meanMDistr m

addResult :: Num a => Result a -> Result a -> Maybe (Result a)
addResult (RConst x) (RConst y) = Just $ RConst (x + y)
addResult (RConst x) (RMarginal (MGaussian mu var)) = Just $ RMarginal (MGaussian (mu + x) var)
addResult (RMarginal (MGaussian mu var)) (RConst x) = addResult (RConst x) (RMarginal (MGaussian mu var))
addResult (RMarginal (MGaussian mu var)) (RMarginal (MGaussian mu' var')) = Nothing
  --  Just $ RMarginal (MGaussian (mu + mu') (var + var')) -- only valid when variables are independent
addResult _ _ = Nothing

negateResult :: Num a => Result a -> Maybe (Result a)
negateResult (RConst x) = Just $ RConst (- x)
negateResult (RMarginal (MGaussian mu var)) = Just $ RMarginal (MGaussian (- mu) var)
negateResult _ = Nothing

multResult :: Num a => Result a -> Result a -> Maybe (Result a)
multResult (RConst x) (RConst y) = Just $ RConst (x * y)
multResult (RConst x) (RMarginal (MGaussian mu var)) = Just $ RMarginal (MGaussian (x * mu) (x^2 * var))
multResult (RMarginal (MGaussian mu var)) (RConst x) = multResult (RConst x) (RMarginal (MGaussian mu var))
multResult _ _ = Nothing

data Expr' var a where
  Const :: a -> Expr' var a
  Var :: var a -> Expr' var a
  Plus :: Num a => Expr' var a -> Expr' var a -> Expr' var a
  Times :: Num a => Expr' var a -> Expr' var a -> Expr' var a
  MVMul :: (KnownNat m, KnownNat n) => Expr' var (L m n) -> Expr' var (R n) -> Expr' var (R m)

instance Show (RefNodeToT a) where
  show (RefNodeToT r) = show r

deriving instance (Show a) => Show (Expr' RefNodeToT a)

class DeepForce e where
  type Forced e :: *
  deepForce :: MonadState Heap m => MonadSample m => e -> m (Forced e)
  deepConst :: Forced e -> e

deepForce' :: MonadState Heap m => MonadSample m => DeepForce e => e -> m e
deepForce' = fmap deepConst . deepForce

instance DeepForce Double where
  type Forced Double = Double
  deepForce = pure
  deepConst = id

instance DeepForce Int where
  type Forced Int = Int
  deepForce = pure
  deepConst = id

instance DeepForce (R n) where
  type Forced (R n) = R n
  deepForce = pure
  deepConst = id

instance DeepForce (L m n) where
  type Forced (L m n) = L m n
  deepForce = pure
  deepConst = id

instance DeepForce (Expr a) where
  type Forced (Expr a) = a
  deepForce = force
  deepConst = Const

instance DeepForce e => DeepForce [e] where
  type Forced [e] = [Forced e]
  deepForce = mapM deepForce
  deepConst = map deepConst

instance DeepForce e => DeepForce (Data.Map.Map k e) where
  type Forced (Data.Map.Map k e) = Data.Map.Map k (Forced e)
  deepForce = mapM deepForce
  deepConst = fmap deepConst

instance (DeepForce a, DeepForce b) => DeepForce (a, b) where
  type Forced (a, b) = (Forced a, Forced b)
  deepForce (x, y) = (,) <$> deepForce x <*> deepForce y
  deepConst (x, y) = (deepConst x, deepConst y)

instance (DeepForce a, DeepForce b, DeepForce c) => DeepForce (a, b, c) where
  type Forced (a, b, c) = (Forced a, Forced b, Forced c)
  deepForce (x, y, z) = (,,) <$> deepForce x <*> deepForce y <*> deepForce z
  deepConst (x, y, z) = (deepConst x, deepConst y, deepConst z)

instance (DeepForce a, DeepForce b, DeepForce c, DeepForce d) => DeepForce (a, b, c, d) where
  type Forced (a, b, c, d) = (Forced a, Forced b, Forced c, Forced d)
  deepForce (x, y, z, w) = (,,,) <$> deepForce x <*> deepForce y <*> deepForce z <*> deepForce w
  deepConst (x, y, z, w) = (deepConst x, deepConst y, deepConst z, deepConst w)

type M = Weighted (StateT Heap SamplerIO)

evalM :: M a -> SamplerIO (Double, a)
evalM x = (\(x, Exp w) -> (w, x)) <$> evalStateT (runWeighted x) emptyHeap

forgettableVar :: MonadState Heap m => Ref (Node a b) -> m (Expr (MType b))
forgettableVar nref = do
  --  addFinalizer v (forget nref)
  pure v
  where
  v = Var (RefNodeToT nref)

type Expr = Expr' RefNodeToT

data RefNodeToT b where
  RefNodeToT :: Ref (Node a b) -> RefNodeToT (MType b)

instance Eq (RefNodeToT a) where
  RefNodeToT (Ref x) == RefNodeToT (Ref y) = x == unsafeCoerce y

instance Num a => Num (Expr' var a) where
  fromInteger = Const . fromInteger
  x + y = Plus x y
  x * y = Times x y

fmapExpr' :: forall f var var' a. Applicative f => (forall x. var x -> f (var' x)) -> Expr' var a -> f (Expr' var' a)
fmapExpr' varMap = f where
  f :: forall a. Expr' var a -> f (Expr' var' a)
  f (Const x) = pure $ Const x
  f (Var i) = Var <$> varMap i
  f (Plus x y) = Plus <$> f x <*> f y
  f (Times x y) = Times <$> f x <*> f y
  f (MVMul a x) = MVMul <$> f a <*> f x

read :: MonadState Heap m => Expr' RefNodeToT a -> m (Expr' NodeTo a)
read = fmapExpr' (\(RefNodeToT nref) -> NodeTo <$> readRef nref)

data AffineMVMult (m :: Nat) where
  AffineMVMult :: KnownNat n => RefNodeToT (R n) -> L m n -> AffineMVMult m

getAffineMV :: KnownNat m => Expr (R m) -> Maybe ([AffineMVMult m], R m)
getAffineMV (Const x) = Just ([], x)
getAffineMV (Var i) = Just ([AffineMVMult i eye], 0)
getAffineMV (Plus x y) = do
  (m1, b1) <- getAffineMV x
  (m2, b2) <- getAffineMV y
  Just (m1 ++ m2, b1 + b2)
getAffineMV (Times x y) = Nothing
getAffineMV (MVMul (Const f) x) = do
  (m, b) <- getAffineMV x
  Just ([ AffineMVMult z (f <> g) | AffineMVMult z g <- m ], f #> b)
getAffineMV (MVMul _ _) = Nothing

getAffineMV1 :: KnownNat m => Expr (R m) -> Maybe (R m, Maybe (AffineMVMult m))
getAffineMV1 e = getAffineMV e >>= toAffine
  where
  toAffine :: KnownNat m => ([AffineMVMult m], R m) -> Maybe (R m, Maybe (AffineMVMult m))
  toAffine (m, b) = do
    mx <- case m of
      [] -> Just Nothing
      [x] -> Just (Just x)
      _ -> Nothing
    pure (b, mx)

getAffine :: Num a => Expr a -> Maybe (Map (Maybe (RefNodeToT a)) a)
getAffine (Const x) = Just (M.singleton Nothing x)
getAffine (Var i) = Just (M.singleton (Just i) 1)
getAffine (Plus x y) = M.unionWith (+) <$> getAffine x <*> getAffine y
getAffine (Times x y) = do
  mx <- getAffine x
  my <- getAffine y
  terms <- sequence [ mult x1 x2 | x1 <- M.toList mx, x2 <- M.toList my ]
  Just (M.fromList terms)
  where
  mult (Nothing, x) (Nothing, y) = Just (Nothing, x * y)
  mult (Nothing, x) (Just i, y) = Just (Just i, x * y)
  mult (Just i, x) (Nothing, y) = Just (Just i, x * y)
  mult (Just i, x) (Just j, y) = Nothing
getAffine (MVMul f x) = Nothing

getAffine1 :: Num a => Expr a -> Maybe (a, Maybe (RefNodeToT a, a))
getAffine1 e = getAffine e >>= toAffine

toAffine :: Num v => Map (Maybe k) v -> Maybe (v, Maybe (k, v))
toAffine mp = let (mconst, mvar) = M.partition Data.Maybe.isNothing mp in do
  let b = case mconst of
            M.Map [] -> 0
            M.Map [(Nothing, v)] -> v
            _ -> error "impossible"
  mx <- case mvar of
    M.Map [] -> Just Nothing
    M.Map [(Just k, v)] -> Just (Just (k, v))
    _ -> Nothing
  pure (b, mx)

realizedToConst' :: forall var f a. Applicative f => (forall x. var x -> f (Maybe x)) -> Expr' var a -> f (Expr' var a)
realizedToConst' varMap = f where
  f :: forall a. Expr' var a -> f (Expr' var a)
  f (Const x) = pure $ Const x
  f (Var i) = flip fmap (varMap i) $ \mk -> case mk of
      Nothing -> Var i
      Just k -> Const k
  f (Plus x y) = Plus <$> f x <*> f y
  f (Times x y) = Times <$> f x <*> f y
  f (MVMul a x) = MVMul <$> f a <*> f x

realizedToConst :: MonadState Heap m => Expr a -> m (Expr a)
realizedToConst = realizedToConst' f where
  f :: MonadState Heap m => RefNodeToT b -> m (Maybe b)
  f (RefNodeToT nref) = do
    n <- readRef nref
    pure $ case n of
      RealizedNode x -> Just x
      _ -> Nothing

allVars :: Expr a -> Set (SomeRefNode)
allVars (Const x) = S.empty
allVars (Var (RefNodeToT i)) = S.singleton (SomeRefNode i)
allVars (Plus x y) = allVars x `S.union` allVars y
allVars (Times x y) = allVars x `S.union` allVars y
allVars (MVMul f x) = allVars f `S.union` allVars x

initializedMarginal :: MonadState Heap m => Ref (Node a b) -> m (MDistr b)
initializedMarginal nref = do
  n <- readRef nref
  case distr n of
    UDistr d -> pure d
    CDistr par cd -> do
      Just pard <- marginalRef par
      pure (makeMarginal pard cd)

marginalRef :: MonadState Heap m => Ref (Node a b) -> m (Maybe (MDistr b))
marginalRef nref = do
  n <- readRef nref
  case n of
    RealizedNode x -> error "shouldn't occur"
    _ -> case state n of
      Marginalized d -> do
        isStale <- stale nref
        pure $ if isStale then Nothing else Just d
      Initialized -> Just <$> initializedMarginal nref

-- This is currently unsound
-- Consider
-- x <- normal(0, 1)
-- pure (x + (- x))
marginal :: MonadState Heap m => Expr a -> m (Maybe (Result a))
marginal (Const x) = pure (Just (RConst x))
marginal (Var (RefNodeToT nref)) = do
  n <- readRef nref
  case n of
    RealizedNode x -> pure $ Just (RConst x)
    _ -> case state n of
      Marginalized d -> do
        isStale <- stale nref
        pure $ if isStale then Nothing else Just (RMarginal d)
      Initialized -> Just . RMarginal <$> initializedMarginal nref
marginal (Plus x y) = pure Nothing -- could do independence ananlysis, etc.
marginal (Times x y) = pure Nothing -- could do independence ananlysis, etc.
marginal (MVMul f x) = pure Nothing

marginal' :: MonadState Heap m => Expr a -> m (Result a)
marginal' = fmap fromJust . marginal

data NodeTo a where
  NodeTo :: Node z a -> NodeTo (MType a)

force :: MonadSample m => MonadState Heap m => Expr a -> m a
force (Const x) = pure x
force (Var (RefNodeToT nref)) = getValue nref
force (Plus x y) = (+) <$> force x <*> force y
force (Times x y) = (*) <$> force x <*> force y
force (MVMul a x) = (#>) <$> force a <*> force x

typeOfRefNodeToT :: MonadState Heap m => Ref (Node a b) -> m (SMarginalT b)
typeOfRefNodeToT nref = typeOfDSDistr . distr <$> readRef nref

forgetE :: MonadState Heap m => Expr a -> m ()
forgetE e = forM_ (allVars e) forgetSomeRefNode

forgetSomeRefNode :: MonadState Heap m => SomeRefNode -> m ()
forgetSomeRefNode (SomeRefNode nref) = forget nref

justRun :: Show p => MonadState Heap m => State.MonadIO m => MStream m (Expr p) (Expr a) -> m (Maybe (Result a))
justRun s = do
  res <- MS.runStream (\x -> marginal x >>= liftIO . print) s
  res' <- marginal res
  pure res'

runParticle :: Int -> MStream M (Maybe (Result Double)) (Maybe (Result Double)) -> StateT Heap SamplerIO Double
runParticle numParticles =
  fmap avg . MS.runStream (liftIO . print . avg)
  . Inference.particles numParticles
  where
  avg = average . map (meanResult . fromJust)
  fromJust (Just x) = x
  fromJust Nothing = RConst $ Prelude.read "NaN"

zdeepForce :: DeepForce b => MonadSample m => ZS.ZStream (StateT Heap m) a b -> ZS.ZStream m a (Forced b)
zdeepForce model = zunheap $ ZS.zconstM deepForce ZS.<<< model

zdeepForce' :: DeepForce b => MonadSample m => ZS.ZStream (StateT Heap m) a b -> ZS.ZStream m a b
zdeepForce' model = fmap deepConst (zdeepForce model)

max_stars_repo_name: psg-mit/probzelus-haskell

max_issues_repo_name: psg-mit/probzelus-haskell

max_forks_repo_name: psg-mit/probzelus-haskell