text:
{- |
Module      :  Camfort.Specification.Units.Analysis
Description :  Helpers for units refactoring and analysis.
Copyright   :  (c) 2017, Dominic Orchard, Andrew Rice, Mistral Contrastin, Matthew Danish
License     :  Apache-2.0

Maintainer  :  dom.orchard@gmail.com
Stability   :  experimental
-}

{-# LANGUAGE OverloadedStrings #-}

module Camfort.Specification.Units.Analysis
  ( UnitAnalysis
  , compileUnits
  , initInference
  , runInference
  , runUnitAnalysis
    -- ** Helpers
  , puName
  , puSrcName
  ) where

import           Camfort.Analysis
import           Camfort.Analysis.Annotations (Annotation)
import           Camfort.Analysis.CommentAnnotator (annotateComments)
import           Camfort.Analysis.Logger (LogLevel(..))
import           Camfort.Analysis.ModFile (withCombinedEnvironment)
import qualified Camfort.Specification.Units.Annotation as UA
import           Camfort.Specification.Units.Environment
import           Camfort.Specification.Units.InferenceBackend
import qualified Camfort.Specification.Units.InferenceBackendFlint as Flint
import qualified Camfort.Specification.Units.InferenceBackendSBV as BackendSBV
import           Camfort.Specification.Units.ModFile
  (genUnitsModFile, initializeModFiles, runCompileUnits)
import           Camfort.Specification.Units.Monad
import           Camfort.Specification.Units.MonadTypes
import           Camfort.Specification.Units.Parser (unitParser)
import qualified Camfort.Specification.Units.Parser.Types as P
import           Control.Lens ((^?), _1)
import           Control.Monad
import           Control.Monad.Reader
import           Control.Monad.State.Strict
import           Control.Monad.Writer.Lazy
import qualified Data.Array as A
import           Data.Data (Data)
import           Data.Generics.Uniplate.Operations
import qualified Data.IntMap.Strict as IM
import           Data.List (nub, intercalate)
import qualified Data.Map.Strict as M
import           Data.Maybe (isJust, fromMaybe, mapMaybe)
import qualified Data.Set as S
import           Data.Text (Text)
import qualified Language.Fortran.AST as F
import           Language.Fortran.Analysis (constExp, varName, srcName)
import qualified Language.Fortran.Analysis as FA
import qualified Language.Fortran.Analysis.BBlocks as FAB
import qualified Language.Fortran.Analysis.DataFlow as FAD
import           Language.Fortran.Parser.Utils (readReal, readInteger)
import           Language.Fortran.Util.ModFile
import qualified Numeric.LinearAlgebra as H -- for debugging
import           Prelude hiding (mod)

-- | Prepare to run an inference function.
initInference :: UnitSolver ()
initInference = do
  pf <- getProgramFile

  -- Parse unit annotations found in comments and link to their
  -- corresponding statements in the AST.
  let (linkedPF, _) =
        runWriter $ annotateComments unitParser
        (\srcSpan err -> tell $ "Error " ++ show srcSpan ++ ": " ++ show err) pf
  modifyProgramFile $ const linkedPF

  -- The following insert* functions examine the AST and insert
  -- mappings into the tables stored in the UnitState.

  -- First, find all given unit annotations and insert them into our
  -- mappings.  Also obtain all unit alias definitions.
  insertGivenUnits

  -- For function or subroutine parameters (or return variables) that
  -- are not given explicit units, give them a parametric polymorphic
  -- unit.
  insertParametricUnits

  -- Any other variables get assigned a unique undetermined unit named
  -- after the variable. This assumes that all variables have unique
  -- names, which the renaming module already has assured.
  insertUndeterminedUnits

  -- Now take the information that we have gathered and annotate the
  -- variable expressions within the AST with it.
  annotateAllVariables

  -- Annotate the literals within the program based upon the
  -- Literals-mode option.
  annotateLiterals

  -- With the variable expressions annotated, we now propagate the
  -- information throughout the AST, giving units to as many
  -- expressions as possible, and also constraints wherever
  -- appropriate.
  propagateUnits

  -- Gather up all of the constraints that we identified in the AST.
  -- These constraints will include parametric polymorphic units that
  -- have not yet been instantiated into their particular uses.
  abstractCons <- extractConstraints
  dumpConsM "***abstractCons" abstractCons

  -- Eliminate all parametric polymorphic units by copying them for
  -- each specific use cases and substituting a unique call-site
  -- identifier that distinguishes each use-case from the others.
  cons <- applyTemplates abstractCons
  dumpConsM "***concreteCons" cons

  -- Remove any traces of CommentAnnotator, since the annotations can
  -- cause generic operations traversing the AST to get confused.
  modifyProgramFile UA.cleanLinks

  modifyConstraints (const cons)

  debugLogging

-- | Run a 'UnitSolver' analysis within a 'UnitsAnalysis'.
runInference :: UnitSolver a -> UnitAnalysis (a, UnitState)
runInference solver = do
  pf <- asks unitProgramFile
  mfs <- lift analysisModFiles

  let (pf', _, _) = withCombinedEnvironment mfs . fmap UA.mkUnitAnnotation $ pf
  let pvm = combinedParamVarMap mfs
  let pf'' = FAD.analyseConstExps . FAD.analyseParameterVars pvm . FAB.analyseBBlocks $ pf'
  runUnitSolver pf'' $ do
    initializeModFiles
    initInference
    solver

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

-- | Seek out any parameters to functions or subroutines that do not
-- already have units, and insert parametric units for them into the
-- map of variables to UnitInfo.
insertParametricUnits :: UnitSolver ()
insertParametricUnits = getProgramFile >>= (mapM_ paramPU . universeBi)
  where
    paramPU pu =
      forM_ (indexedParams pu) $ \ (i, param) ->
        -- Insert a parametric unit if the variable does not already have a unit.
        modifyVarUnitMap $ M.insertWith (curry snd) param (UnitParamPosAbs (fname, i))
      where
        fname = (puName pu, puSrcName pu)

-- | Return the list of parameters paired with its positional index.
indexedParams :: F.ProgramUnit UA -> [(Int, VV)]
indexedParams pu
  | F.PUFunction _ _ _ _ _ Nothing (Just r) _ _       <- pu = [(0, toVV r)]
  | F.PUFunction _ _ _ _ _ Nothing _ _ _              <- pu = [(0, (fname, sfname))]
  | F.PUFunction _ _ _ _ _ (Just paList) (Just r) _ _ <- pu = zip [0..] $ map toVV (r : F.aStrip paList)
  | F.PUFunction _ _ _ _ _ (Just paList) _ _ _        <- pu = zip [0..] $ (fname, sfname) : map toVV (F.aStrip paList)
  | F.PUSubroutine _ _ _ _ (Just paList) _ _          <- pu = zip [1..] $ map toVV (F.aStrip paList)
  | otherwise                                               = []
  where
    fname  = puName pu
    sfname = puSrcName pu
    toVV e = (varName e, srcName e)

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

-- | Any remaining variables with unknown units are given unit UnitVar
-- with a unique name (in this case, taken from the unique name of the
-- variable as provided by the Renamer), or UnitParamVarAbs if the
-- variables are inside of a function or subroutine.
insertUndeterminedUnits :: UnitSolver ()
insertUndeterminedUnits = do
  pf   <- getProgramFile
  dmap <- lift . lift $ M.union (extractDeclMap pf) . combinedDeclMap <$> analysisModFiles
  forM_ (universeBi pf :: [F.ProgramUnit UA]) $ \ pu ->
    modifyPUBlocksM (transformBiM (insertUndeterminedUnitVar dmap)) pu

-- Specifically handle variables
insertUndeterminedUnitVar :: DeclMap -> F.Expression UA -> UnitSolver (F.Expression UA)
insertUndeterminedUnitVar dmap v@(F.ExpValue _ _ (F.ValVariable _))
  | Just (FA.IDType { FA.idVType = Just bty }) <- FA.idType (F.getAnnotation v)
  , bty `elem` acceptableTypes = do
  let vname = varName v
  let sname = srcName v
  let unit  = toUnitVar dmap (vname, sname)
  modifyVarUnitMap $ M.insertWith (curry snd) (varName v, srcName v) unit
  pure v
insertUndeterminedUnitVar _ e = pure e

-- Choose UnitVar or UnitParamVarAbs depending upon how the variable was declared.
toUnitVar :: DeclMap -> VV -> UnitInfo
toUnitVar dmap (vname, sname) = unit
  where
    unit = case fst <$> M.lookup vname dmap of
      Just (DCFunction (F.Named fvname, F.Named fsname))   -> UnitParamVarAbs ((fvname, fsname), (vname, sname))
      Just (DCSubroutine (F.Named fvname, F.Named fsname)) -> UnitParamVarAbs ((fvname, fsname), (vname, sname))
      _                                                    -> UnitVar (vname, sname)

-- Insert undetermined units annotations on the following types of variables.
acceptableTypes :: [F.BaseType]
acceptableTypes = [F.TypeReal, F.TypeDoublePrecision, F.TypeComplex, F.TypeDoubleComplex, F.TypeInteger]

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

-- | Convert explicit polymorphic annotations such as (UnitName "'a")
-- into UnitParamEAPAbs with a 'context-unique-name' given by the
-- ProgramUnitName combined with the supplied unit name.
transformExplicitPolymorphism :: Maybe F.ProgramUnitName -> UnitInfo -> UnitInfo
transformExplicitPolymorphism (Just (F.Named f)) (UnitName a@('\'':_)) = UnitParamEAPAbs (a, f ++ "_" ++ a)
transformExplicitPolymorphism _ u                                      = u

-- | Any units provided by the programmer through comment annotations
-- will be incorporated into the VarUnitMap.
insertGivenUnits :: UnitSolver ()
insertGivenUnits = do
  pf <- getProgramFile
  mapM_ checkPU (universeBi pf)
  where
    -- Look through each Program Unit for the comments
    checkPU :: F.ProgramUnit UA -> UnitSolver ()
    checkPU (F.PUComment a _ _)
      -- Look at unit assignment between function return variable and spec.
      | Just (P.UnitAssignment (Just vars) unitsAST) <- mSpec
      , Just pu                                      <- mPU = insertPUUnitAssigns (toUnitInfo unitsAST) pu vars
      -- Add a new unit alias.
      | Just (P.UnitAlias name unitsAST) <- mSpec = modifyUnitAliasMap (M.insert name (toUnitInfo unitsAST))
      | otherwise                                 = pure ()
      where
        mSpec = UA.unitSpec (FA.prevAnnotation a)
        mPU   = UA.unitPU   (FA.prevAnnotation a)
    -- Other type of ProgramUnit (e.g. one with a body of blocks)
    checkPU pu = mapM_ (checkBlockComment getName) [ b | b@F.BlComment{} <- universeBi (F.programUnitBody pu) ]
      where
        getName = case pu of
          F.PUFunction {}   -> Just $ F.getName pu
          F.PUSubroutine {} -> Just $ F.getName pu
          _                 -> Nothing

    -- Look through each comment that has some kind of unit annotation within it.
    checkBlockComment :: Maybe F.ProgramUnitName -> F.Block UA -> UnitSolver ()
    checkBlockComment pname (F.BlComment a _ _)
      -- Look at unit assignment between variable and spec.
      | Just (P.UnitAssignment (Just vars) unitsAST) <- mSpec
      , Just b                                       <- mBlock = insertBlockUnitAssigns pname (toUnitInfo unitsAST) b vars
      -- Add a new unit alias.
      | Just (P.UnitAlias name unitsAST) <- mSpec  = modifyUnitAliasMap (M.insert name (toUnitInfo unitsAST))
      | otherwise                                  = pure ()
      where
        mSpec  = UA.unitSpec  (FA.prevAnnotation a)
        mBlock = UA.unitBlock (FA.prevAnnotation a)
    checkBlockComment _ _ = error "received non-comment in checkBlockComment"

    -- Figure out the unique names of the referenced variables and
    -- then insert unit info under each of those names.
    insertBlockUnitAssigns :: Maybe F.ProgramUnitName -> UnitInfo -> F.Block UA -> [String] -> UnitSolver ()
    insertBlockUnitAssigns pname info (F.BlStatement _ _ _ (F.StDeclaration _ _ _ _ decls)) varRealNames = do
      -- figure out the 'unique name' of the varRealName that was found in the comment
      -- FIXME: account for module renaming
      -- FIXME: might be more efficient to allow access to variable renaming environ at this program point
      let info' = transform (transformExplicitPolymorphism pname) info
      let m = M.fromList [ ((varName e, srcName e), info')
                         | e@(F.ExpValue _ _ (F.ValVariable _)) <- universeBi decls :: [F.Expression UA]
                         , varRealName <- varRealNames
                         , varRealName == srcName e ]
      modifyVarUnitMap $ M.unionWith const m
      modifyGivenVarSet . S.union . S.fromList . map fst . M.keys $ m
    insertBlockUnitAssigns _ _ _ _ = error "received non-statement/declaration in insertBlockUnitAssigns"

    -- Insert unit annotation for function return variable
    insertPUUnitAssigns :: UnitInfo -> F.ProgramUnit UA -> [String] -> UnitSolver ()
    insertPUUnitAssigns info pu@(F.PUFunction _ _ _ _ _ _ mret _ _) varRealNames
      | (retUniq, retSrc) <- case mret of Just ret -> (FA.varName ret, FA.srcName ret)
                                          Nothing  -> (puName pu, puSrcName pu)
      , retSrc `elem` varRealNames = do
          let pname = Just $ F.getName pu
          let info' = transform (transformExplicitPolymorphism pname) info
          let m = M.fromList [ ((retUniq, retSrc), info') ]
          modifyVarUnitMap $ M.unionWith const m
          modifyGivenVarSet . S.union . S.fromList . map fst . M.keys $ m
    insertPUUnitAssigns _ _ _ = error "received non-function in insertPUUnitAssigns"

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

-- | Take the unit information from the VarUnitMap and use it to
-- annotate every variable expression in the AST.
annotateAllVariables :: UnitSolver ()
annotateAllVariables = modifyProgramFileM $ \ pf -> do
  varUnitMap <- getVarUnitMap
  importedVariables <- getImportedVariables
  let varUnitMap' = M.unionWith (curry snd) varUnitMap importedVariables

  let annotateExp e@(F.ExpValue _ _ (F.ValVariable _))
        | Just info <- M.lookup (varName e, srcName e) varUnitMap' = UA.setUnitInfo info e
      -- may need to annotate intrinsics separately
      annotateExp e = e
  pure $ transformBi annotateExp pf

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

-- | Give units to literals based upon the rules of the Literals mode.
--
-- LitUnitless: All literals are unitless.
-- LitPoly:     All literals are polymorphic.
-- LitMixed:    The literal "0" or "0.0" is fully parametric polymorphic.
--              All other literals are monomorphic, possibly unitless.
annotateLiterals :: UnitSolver ()
annotateLiterals = modifyProgramFileM (transformBiM annotateLiteralsPU)

annotateLiteralsPU :: F.ProgramUnit UA -> UnitSolver (F.ProgramUnit UA)
annotateLiteralsPU pu = do
  mode <- asks (uoLiterals . unitOpts)
  case mode of
    LitUnitless -> modifyPUBlocksM (transformBiM expUnitless) pu
    LitPoly     -> modifyPUBlocksM (transformBiM (withLiterals genParamLit)) pu
    LitMixed    -> modifyPUBlocksM (transformBiM expMixed) pu
  where
    -- Follow the LitMixed rules.
    expMixed e = case e of
      F.ExpValue _ _ (F.ValInteger i) | readInteger i == Just 0 -> withLiterals genParamLit e
                                      | otherwise               -> withLiterals genUnitLiteral e
      F.ExpValue _ _ (F.ValReal i) | readReal i == Just 0       -> withLiterals genParamLit e
                                   | otherwise                  -> withLiterals genUnitLiteral e

      F.ExpBinary a s op e1 e2
        | op `elem` [F.Multiplication, F.Division] -> case () of
            -- leave it alone if they're both constants
            _ | Just _ <- constExp (F.getAnnotation e1)
              , Just _ <- constExp (F.getAnnotation e2)        -> pure e
            -- a constant multiplier is unitless
            _ | Just _ <- constExp (F.getAnnotation e1)
              , Just UnitLiteral{} <- UA.getUnitInfo e1         ->
                pure $ F.ExpBinary a s op (UA.setUnitInfo UnitlessLit e1) e2
            -- a constant multiplier is unitless
              | Just _ <- constExp (F.getAnnotation e2)
              , Just UnitLiteral{} <- UA.getUnitInfo e2         ->
                pure $ F.ExpBinary a s op e1 (UA.setUnitInfo UnitlessLit e2)
            _                                                   -> pure e

      _ | Just _ <- constExp (F.getAnnotation e)                -> case UA.getUnitInfo e of
            -- Treat constant expressions as if they were fresh
            -- literals, unless assigned units already.
            Just UnitLiteral{} -> genLit e
            Just UnitVar{}     -> genLit e
            _                  -> pure e
        | otherwise                                             -> pure e

    -- Set all literals to unitless.
    expUnitless e
      | isLiteral e = pure $ UA.setUnitInfo UnitlessLit e
      | otherwise   = pure e

    -- Set all literals to the result of given monadic computation.
    withLiterals m e
      | isLiteral e = flip UA.setUnitInfo e <$> m
      | otherwise   = pure e

    -- isPolyCtxt = case pu of F.PUFunction {} -> True; F.PUSubroutine {} -> True; _ -> False

    genLit e
      | isLiteralZero e = withLiterals genParamLit e
      | otherwise       = withLiterals genUnitLiteral e

-- | Is it a literal, literally?
isLiteral :: F.Expression UA -> Bool
isLiteral (F.ExpValue _ _ (F.ValReal _))    = True
isLiteral (F.ExpValue _ _ (F.ValInteger _)) = True
-- allow propagated constants to be interpreted as literals
isLiteral e                                 = isJust (constExp (F.getAnnotation e))

-- | Is expression a literal and is it non-zero?
isLiteralNonZero :: F.Expression UA -> Bool
isLiteralNonZero (F.ExpValue _ _ (F.ValInteger i)) = readInteger i /= Just 0
isLiteralNonZero (F.ExpValue _ _ (F.ValReal i))    = readReal i    /= Just 0
-- allow propagated constants to be interpreted as literals
isLiteralNonZero e = case constExp (F.getAnnotation e) of
  Just (FA.ConstInt i)          -> i /= 0
  Just (FA.ConstUninterpInt s)  -> readInteger s /= Just 0
  Just (FA.ConstUninterpReal s) -> readReal s /= Just 0
  _                             -> False

isLiteralZero :: F.Expression UA -> Bool
isLiteralZero x = isLiteral x && not (isLiteralNonZero x)

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

-- | Filter out redundant constraints.
cullRedundant :: Constraints -> Constraints
cullRedundant = nub . mapMaybe ( \ con -> case con of
  ConEq u1 u2 | u1 /= u2                              -> Just con
  ConConj cs | cs' <- cullRedundant cs, not (null cs) -> Just (ConConj cs')
  _                                                   -> Nothing
  )

-- | Convert all parametric templates into actual uses, via substitution.
applyTemplates :: Constraints -> UnitSolver Constraints
-- postcondition: returned constraints lack all Parametric constructors
applyTemplates cons = do
  dumpConsM "applyTemplates" cons
  -- Get a list of the instances of parametric polymorphism from the constraints.
  let instances = nub [ (name, i) | UnitParamPosUse ((name, _), _, i) <- universeBi cons ]

  -- Also generate a list of 'dummy' instances to ensure that every
  -- 'toplevel' function and subroutine is thoroughly expanded and
  -- analysed, even if it is not used in the current ProgramFile. (It
  -- might be part of a library module, for instance).
  pf <- getProgramFile
  dummies <- forM (topLevelFuncsAndSubs pf) $ \ pu -> do
    ident <- freshId
    pure (puName pu, ident)

  logDebug' pf $ ("instances: " <> describeShow instances)
  logDebug' pf $ ("dummies: " <> describeShow dummies)

  importedVariables <- getImportedVariables

  -- Prepare constraints for all variables imported via StUse.
  let importedCons = [ ConEq (UnitVar vv) units | (vv, units) <- M.toList importedVariables ]

  -- Work through the instances, expanding their templates, and
  -- substituting the callId into the abstract parameters.
  concreteCons <- cullRedundant <$>
                  liftM2 (++) (foldM (substInstance False []) importedCons instances)
                              (foldM (substInstance True []) [] dummies)
  dumpConsM "applyTemplates: concreteCons" concreteCons

  -- Also include aliases in the final set of constraints, where
  -- aliases are implemented by simply asserting that they are equal
  -- to their definition.
  aliasMap <- getUnitAliasMap
  let aliases = [ ConEq (UnitAlias name) def | (name, def) <- M.toList aliasMap ]
  let transAlias (UnitName a) | a `M.member` aliasMap = UnitAlias a
      transAlias u                                    = u

  dumpConsM "aliases" aliases
  pure . transformBi transAlias . cullRedundant $ cons ++ concreteCons ++ aliases

-- | Look up the Parametric templates for a given function or
-- subroutine, and do the substitutions. Process any additional
-- polymorphic calls that are uncovered, unless they are recursive
-- calls that have already been seen in the current call stack.
substInstance :: Bool -> [F.Name] -> Constraints -> (F.Name, Int) -> UnitSolver Constraints
substInstance isDummy callStack output (name, callId) = do
  tmap <- getTemplateMap

  -- Look up the templates associated with the given function or
  -- subroutine name. And then transform the templates by generating
  -- new callIds for any constraints created by function or subroutine
  -- calls contained within the templates.
  --
  -- The reason for this is because functions called by functions can
  -- be used in a parametric polymorphic way.

  -- npc <- nameParamConstraints name -- In case it is an imported function, use this.
  let npc = [] -- disabled for now
  template <- transformBiM callIdRemap $ npc `fromMaybe` M.lookup name tmap
  dumpConsM ("substInstance " ++ show isDummy ++ " " ++ show callStack ++ " " ++ show (name, callId) ++ " template lookup") template

  -- Reset the usCallIdRemap field so that it is ready for the next
  -- set of templates.
  modifyCallIdRemap (const IM.empty)

  -- If any new instances are discovered, also process them, unless recursive.
  let instances = nub [ (name', i) | UnitParamPosUse ((name', _), _, i) <- universeBi template ]
  template' <- if name `elem` callStack then
                 -- Detected recursion: we do not support polymorphic-unit recursion,
                 -- ergo all subsequent recursive calls are assumed to have the same
                 -- unit-assignments as the first call.
                 pure []
               else
                 foldM (substInstance False (name:callStack)) [] instances

  dumpConsM ("instantiating " ++ show (name, callId) ++ ": (output ++ template) is") (output ++ template)
  dumpConsM ("instantiating " ++ show (name, callId) ++ ": (template') is") template'

  -- Convert abstract parametric units into concrete ones.

  let output' = -- Do not instantiate explicitly annotated polymorphic
                -- variables from current context when looking at dummy (name, callId)
                (if isDummy then output ++ template
                            else instantiate callId (output ++ template)) ++

                -- Only instantiate explicitly annotated polymorphic
                -- variables from nested function/subroutine calls.
                instantiate callId template'

  dumpConsM ("final output for " ++ show (name, callId)) output'

  pure output'

-- -- | Generate constraints from a NameParamMap entry.
-- nameParamConstraints :: F.Name -> UnitSolver Constraints
-- nameParamConstraints fname = do
--   let filterForName (NPKParam (n, _) _) _ = n == fname
--       filterForName _ _                   = False
--   nlst <- (M.toList . M.filterWithKey filterForName) <$> getNameParamMap
--   pure [ ConEq (UnitParamPosAbs (n, pos)) (foldUnits units) | (NPKParam n pos, units) <- nlst ]

-- | If given a usage of a parametric unit, rewrite the callId field
-- to follow an existing mapping in the usCallIdRemap state field, or
-- generate a new callId and add it to the usCallIdRemap state field.
callIdRemap :: UnitInfo -> UnitSolver UnitInfo
callIdRemap info = modifyCallIdRemapM $ \ idMap -> case info of
    UnitParamPosUse (n, p, i)
      | Just i' <- IM.lookup i idMap -> pure (UnitParamPosUse (n, p, i'), idMap)
      | otherwise                    -> freshId >>= \ i' ->
                                          pure (UnitParamPosUse (n, p, i'), IM.insert i i' idMap)
    UnitParamVarUse (n, v, i)
      | Just i' <- IM.lookup i idMap -> pure (UnitParamVarUse (n, v, i'), idMap)
      | otherwise                    -> freshId >>= \ i' ->
                                          pure (UnitParamVarUse (n, v, i'), IM.insert i i' idMap)
    UnitParamLitUse (l, i)
      | Just i' <- IM.lookup i idMap -> pure (UnitParamLitUse (l, i'), idMap)
      | otherwise                    -> freshId >>= \ i' ->
                                          pure (UnitParamLitUse (l, i'), IM.insert i i' idMap)
    UnitParamEAPUse (v, i)
      | Just i' <- IM.lookup i idMap -> pure (UnitParamEAPUse (v, i'), idMap)
      | otherwise                    -> freshId >>= \ i' ->
                                          pure (UnitParamEAPUse (v, i'), IM.insert i i' idMap)

    _                                -> pure (info, idMap)

-- | Convert a parametric template into a particular use.
instantiate :: Data a => Int -> a -> a
instantiate callId = transformBi $ \ info -> case info of
  UnitParamPosAbs (name, position) -> UnitParamPosUse (name, position, callId)
  UnitParamLitAbs litId            -> UnitParamLitUse (litId, callId)
  UnitParamVarAbs (fname, vname)   -> UnitParamVarUse (fname, vname, callId)
  UnitParamEAPAbs vname            -> UnitParamEAPUse (vname, callId)
  _                                -> info

-- | Return a list of ProgramUnits that might be considered 'toplevel'
-- in the ProgramFile, e.g., possible exports. These must be analysed
-- independently of whether they are actually used in the same file,
-- because other files might use them.
topLevelFuncsAndSubs :: F.ProgramFile a -> [F.ProgramUnit a]
topLevelFuncsAndSubs (F.ProgramFile _ pus) = topLevel =<< pus
  where
    topLevel (F.PUModule _ _ _ _ (Just contains)) = topLevel =<< contains
    topLevel (F.PUMain _ _ _ _ (Just contains))   = topLevel =<< contains
    topLevel f@F.PUFunction{}                  = pure f
    topLevel s@F.PUSubroutine{}                = pure s
    topLevel _                                    = []

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

-- | Gather all constraints from the main blocks of the AST, as well as from the varUnitMap
extractConstraints :: UnitSolver Constraints
extractConstraints = do
  pf         <- getProgramFile
  dmap       <- lift . lift $ M.union (extractDeclMap pf) . combinedDeclMap <$> analysisModFiles
  varUnitMap <- getVarUnitMap
  pure $ [ con | b <- mainBlocks pf, con@ConEq{} <- universeBi b ] ++
           [ ConEq (toUnitVar dmap v) u | (v, u) <- M.toList varUnitMap ]

-- | A list of blocks considered to be part of the 'main' program.
mainBlocks :: F.ProgramFile UA -> [F.Block UA]
mainBlocks = concatMap getBlocks . universeBi
  where
    getBlocks (F.PUMain _ _ _ bs _)   = bs
    getBlocks (F.PUModule _ _ _ bs _) = bs
    getBlocks _                       = []

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

-- | Propagate* functions: decorate the AST with constraints, given
-- that variables have all been annotated.
propagateUnits :: UnitSolver ()
-- precondition: all variables have already been annotated
propagateUnits = modifyProgramFileM $ transformBiM propagateInterface <=<
                                      transformBiM propagatePU        <=<
                                      transformBiM propagateDoSpec    <=<
                                      transformBiM propagateStatement <=<
                                      transformBiM propagateExp

propagateExp :: F.Expression UA -> UnitSolver (F.Expression UA)
propagateExp e = case e of
  F.ExpValue{}                           -> pure e -- all values should already be annotated
  F.ExpBinary _ _ F.Multiplication e1 e2 -> setF2 UnitMul (UA.getUnitInfo e1) (UA.getUnitInfo e2)
  F.ExpBinary _ _ F.Division e1 e2       -> setF2 UnitMul (UA.getUnitInfo e1) (flip UnitPow (-1) <$> UA.getUnitInfo e2)
  F.ExpBinary _ _ F.Exponentiation e1 e2 -> setF2 UnitPow (UA.getUnitInfo e1) (constantExpression e2)
  F.ExpBinary _ _ o e1 e2 | isOp AddOp o -> setF2C ConEq  (UA.getUnitInfo e1) (UA.getUnitInfo e2)
                          | isOp RelOp o -> setF2C ConEq  (UA.getUnitInfo e1) (UA.getUnitInfo e2)
  F.ExpFunctionCall {}                   -> propagateFunctionCall e
  F.ExpSubscript _ _ e1 _                -> pure $ UA.maybeSetUnitInfo (UA.getUnitInfo e1) e
  F.ExpUnary _ _ _ e1                    -> pure $ UA.maybeSetUnitInfo (UA.getUnitInfo e1) e
  F.ExpInitialisation{}                  -> pure e
  _                                      -> do
    logDebug' e $ "progagateExp: unhandled " <> describeShow e
    pure e
  where
    -- Shorter names for convenience functions.
    setF2 f u1 u2  = pure $ UA.maybeSetUnitInfoF2 f u1 u2 e
    -- Remember, not only set a constraint, but also give a unit!
    setF2C f u1 u2 = pure . UA.maybeSetUnitInfo u1 $ UA.maybeSetUnitConstraintF2 f u1 u2 e

propagateFunctionCall :: F.Expression UA -> UnitSolver (F.Expression UA)
propagateFunctionCall (F.ExpFunctionCall a s f Nothing)                     = do
  (info, _) <- callHelper f []
  let cons = intrinsicHelper info f []
  pure . UA.setConstraint (ConConj cons) . UA.setUnitInfo info $ F.ExpFunctionCall a s f Nothing
propagateFunctionCall (F.ExpFunctionCall a s f (Just (F.AList a' s' args))) = do
  (info, args') <- callHelper f args
  let cons = intrinsicHelper info f args'
  pure . UA.setConstraint (ConConj cons) . UA.setUnitInfo info $ F.ExpFunctionCall a s f (Just (F.AList a' s' args'))
propagateFunctionCall _ = error "received non-function-call in propagateFunctionCall"

propagateDoSpec :: F.DoSpecification UA -> UnitSolver (F.DoSpecification UA)
propagateDoSpec ast@(F.DoSpecification _ _ (F.StExpressionAssign _ _ e1 _) e2 m_e3) = do
  -- express constraints between the iteration variable, the bounding
  -- expressions and the step expression, or treat the step expression
  -- as a literal 1 if not specified.
  pure . maybe ast (flip UA.setConstraint ast) $ ConConj <$> mconcat [
        -- units(e1) ~ units(e2)
        (:[]) <$> liftM2 ConEq (UA.getUnitInfo e1) (UA.getUnitInfo e2)

        -- units(e1) ~ units(e3) or if e3 not specified then units(e1) ~ 1 in a polymorphic context
        , do u1 <- UA.getUnitInfo e1
             u3 <- (UA.getUnitInfo =<< m_e3) `mplus` if isMonomorphic u1 then mzero else pure UnitlessVar
             pure [ConEq u1 u3]

        -- units(e2) ~ units(e3) or if e3 not specified then units(e2) ~ 1 in a polymorphic context
        , do u2 <- UA.getUnitInfo e1
             u3 <- (UA.getUnitInfo =<< m_e3) `mplus` if isMonomorphic u2 then mzero else pure UnitlessVar
             pure [ConEq u2 u3]
        ]
propagateDoSpec _ = error "propagateDoSpec: called on invalid DoSpec"

propagateStatement :: F.Statement UA -> UnitSolver (F.Statement UA)
propagateStatement stmt = case stmt of
  F.StExpressionAssign _ _ e1 e2               -> literalAssignmentSpecialCase e1 e2 stmt
  F.StCall a s sub (Just (F.AList a' s' args)) -> do
    (info, args') <- callHelper sub args
    let cons = intrinsicHelper info sub args'
    pure . UA.setConstraint (ConConj cons) $ F.StCall a s sub (Just (F.AList a' s' args'))
  F.StDeclaration {}                           -> transformBiM propagateDeclarator stmt
  _                                            -> pure stmt

propagateDeclarator :: F.Declarator UA -> UnitSolver (F.Declarator UA)
propagateDeclarator decl = case decl of
  F.DeclVariable _ _ e1 _ (Just e2) -> literalAssignmentSpecialCase e1 e2 decl
  F.DeclArray _ _ e1 _ _ (Just e2)  -> literalAssignmentSpecialCase e1 e2 decl
  _                                 -> pure decl

-- Allow literal assignment to overload the non-polymorphic
-- unit-assignment of the non-zero literal.
literalAssignmentSpecialCase :: (F.Annotated f)
                             => F.Expression UA -> F.Expression UA
                             -> f UA -> UnitSolver (f UA)
literalAssignmentSpecialCase e1 e2 ast
  | isLiteralZero e2 = pure ast
  | isLiteral e2
  , Just u1            <- UA.getUnitInfo e1
  , Just UnitLiteral{} <- UA.getUnitInfo e2
  , isMonomorphic u1 = pure ast
  -- otherwise express the constraint between LHS and RHS of assignment.
  | otherwise = pure $ UA.maybeSetUnitConstraintF2 ConEq (UA.getUnitInfo e1) (UA.getUnitInfo e2) ast

-- Generic Interface template mapping will be same as first module procedure.
propagateInterface :: F.Block UA -> UnitSolver (F.Block UA)
propagateInterface b@(F.BlInterface _ _ (Just e) _ _ bs) = do
  let iname = varName e
  case [ varName e1 | F.StModuleProcedure _ _ (F.AList _ _ (e1:_)) <- universeBi bs :: [F.Statement UA] ] of
    mpname:_ -> do
      -- translate any instance of mpname into iname within the template
      let trans = transformBi (\ x -> if x == mpname then iname else x)
      -- copy (translated) template from first module procedure to interface
      modifyTemplateMap $ \ m -> fromMaybe m ((\ t -> M.insert iname (trans t) m) <$> M.lookup mpname m)
    _        ->
      pure ()
  pure b

propagateInterface b = pure b

propagatePU :: F.ProgramUnit UA -> UnitSolver (F.ProgramUnit UA)
propagatePU pu = do
  let name     = puName pu
  let sname    = puSrcName pu
  let nn       = (name, sname)
  let bodyCons = [ con | con@ConEq{} <- universeBi pu ] -- Constraints within the PU.

  varMap <- getVarUnitMap

  -- If any of the function/subroutine parameters was given an
  -- explicit unit annotation, then create a constraint between that
  -- explicit unit and the UnitParamPosAbs corresponding to the
  -- parameter. This way all other uses of the parameter get linked to
  -- the explicit unit annotation as well.
  givenCons <- forM (indexedParams pu) $ \ (i, param) ->
    case M.lookup param varMap of
      Just UnitParamPosAbs{}    -> pure . ConEq (UnitParamVarAbs (nn, param)) $ UnitParamPosAbs (nn, i)
      Just u                    -> pure . ConEq u $ UnitParamPosAbs (nn, i)
      _                         -> pure . ConEq (UnitParamVarAbs (nn, param)) $ UnitParamPosAbs (nn, i)

  let cons = givenCons ++ bodyCons
  case pu of F.PUFunction {}   -> modifyTemplateMap (M.insert name cons)
             F.PUSubroutine {} -> modifyTemplateMap (M.insert name cons)
             _                 -> pure ()

  -- Set the unitInfo field of a function program unit to be the same
  -- as the unitInfo of its result.
  let pu' = case (pu, indexedParams pu) of
              (F.PUFunction {}, (0, res):_) -> UA.setUnitInfo (UnitParamPosAbs (nn, 0) `fromMaybe` M.lookup res varMap) pu
              _                             -> pu

  pure (UA.setConstraint (ConConj cons) pu')

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

-- | Coalesce various function and subroutine call common code.
callHelper :: F.Expression UA -> [F.Argument UA] -> UnitSolver (UnitInfo, [F.Argument UA])
callHelper nexp args = do
  let name = (varName nexp, srcName nexp)
  let ctyp = FA.idCType =<< FA.idType (F.getAnnotation nexp)
  callId <- case ctyp of
    Just FA.CTExternal -> pure 0  -- if external with no further info then no polymorphism
    _                  -> freshId -- every call-site gets its own unique identifier
  let eachArg i arg@(F.Argument _ _ _ e)
        -- add site-specific parametric constraints to each argument
        | Just u <- UA.getUnitInfo e = UA.setConstraint (ConEq u (UnitParamPosUse (name, i, callId))) arg
        | otherwise               = arg
  let args' = zipWith eachArg [1..] args
  -- build a site-specific parametric unit for use on a return variable, if any
  let info = UnitParamPosUse (name, 0, callId)
  pure (info, args')

-- FIXME: use this function to create a list of constraints on intrinsic call-sites...
intrinsicHelper :: Foldable t => UnitInfo -> F.Expression (FA.Analysis a) -> t b -> [Constraint]
intrinsicHelper (UnitParamPosUse (_, _, callId)) f@(F.ExpValue _ _ (F.ValIntrinsic _)) args
  | Just (retU, argUs) <- intrinsicLookup sname = zipWith eachArg [0..numArgs] (retU:argUs)
  where
    numArgs     = length args
    sname       = srcName f
    vname       = varName f
    eachArg i u = ConEq (UnitParamPosUse ((vname, sname), i, callId)) (instantiate callId u)
intrinsicHelper _ _ _ = []

-- | Get info about intrinsics by source name 'sname', taking into
-- account the special case of those with arbitrary number of
-- arguments.
intrinsicLookup :: F.Name -> Maybe (UnitInfo, [UnitInfo])
intrinsicLookup sname = do
  (retU, argUs) <- M.lookup sname intrinsicUnits
  return (retU, if sname `S.member` specialCaseArbitraryArgs then cycle argUs else argUs)

-- | Generate a unique identifier for a literal encountered in the code.
genUnitLiteral :: UnitSolver UnitInfo
genUnitLiteral = UnitLiteral <$> freshId

-- | Generate a unique identifier for a polymorphic literal encountered in the code.
genParamLit :: UnitSolver UnitInfo
genParamLit = UnitParamLitAbs <$> freshId

-- Operate only on the blocks of a program unit, not the contained sub-programunits.
modifyPUBlocksM :: Monad m => ([F.Block a] -> m [F.Block a]) -> F.ProgramUnit a -> m (F.ProgramUnit a)
modifyPUBlocksM f pu = case pu of
  F.PUMain a s n b pus                    -> flip fmap (f b) $ \ b' -> F.PUMain a s n b' pus
  F.PUModule a s n b pus                  -> flip fmap (f b) $ \ b' -> F.PUModule a s n b' pus
  F.PUSubroutine a s r n p b subs         -> flip fmap (f b) $ \ b' -> F.PUSubroutine a s r n p b' subs
  F.PUFunction   a s r rec n p res b subs -> flip fmap (f b) $ \ b' -> F.PUFunction a s r rec n p res b' subs
  F.PUBlockData  a s n b                  -> flip fmap (f b) $ \ b' -> F.PUBlockData  a s n b'
  F.PUComment {}                          -> pure pu -- no blocks

-- Fortran semantics for interpretation of constant expressions
-- involving numeric literals.
data FNum = FReal Double | FInt Integer

fnumToDouble :: FNum -> Double
fnumToDouble (FReal x) = x
fnumToDouble (FInt x)  = fromIntegral x

fAdd, fSub, fMul, fDiv, fPow :: FNum -> FNum -> FNum
fAdd (FReal x) fy      = FReal $ x + fnumToDouble fy
fAdd fx (FReal y)      = FReal $ fnumToDouble fx + y
fAdd (FInt x) (FInt y) = FInt  $ x + y
fSub (FReal x) fy      = FReal $ x - fnumToDouble fy
fSub fx (FReal y)      = FReal $ fnumToDouble fx - y
fSub (FInt x) (FInt y) = FInt  $ x - y
fMul (FReal x) fy      = FReal $ x * fnumToDouble fy
fMul fx (FReal y)      = FReal $ fnumToDouble fx * y
fMul (FInt x) (FInt y) = FInt  $ x * y
fDiv (FReal x) fy      = FReal $ x / fnumToDouble fy
fDiv fx (FReal y)      = FReal $ fnumToDouble fx / y
fDiv (FInt x) (FInt y) = FInt  $ x `quot` y  -- Haskell quot truncates towards zero, like Fortran
fPow (FReal x) fy      = FReal $ x ** fnumToDouble fy
fPow fx (FReal y)      = FReal $ fnumToDouble fx ** y
fPow (FInt x) (FInt y)
  | y >= 0             = FInt  $ x ^ y
  | otherwise          = FReal $ fromIntegral x ^^ y

fDivMaybe :: Maybe FNum -> Maybe FNum -> Maybe FNum
fDivMaybe mx my
  | Just y <- my,
    fnumToDouble y == 0.0 = Nothing
  | otherwise             = liftM2 fDiv mx my

-- | Statically computes if the expression is a constant value.
constantExpression :: F.Expression a -> Maybe Double
constantExpression expr = fnumToDouble <$> ce expr
  where
    ce e = case e of
      (F.ExpValue _ _ (F.ValInteger i))        -> FInt <$> readInteger i
      (F.ExpValue _ _ (F.ValReal r))           -> FReal <$> readReal r
      (F.ExpBinary _ _ F.Addition e1 e2)       -> liftM2 fAdd (ce e1) (ce e2)
      (F.ExpBinary _ _ F.Subtraction e1 e2)    -> liftM2 fSub (ce e1) (ce e2)
      (F.ExpBinary _ _ F.Multiplication e1 e2) -> liftM2 fMul (ce e1) (ce e2)
      (F.ExpBinary _ _ F.Division e1 e2)       -> fDivMaybe (ce e1) (ce e2)
      (F.ExpBinary _ _ F.Exponentiation e1 e2) -> liftM2 fPow (ce e1) (ce e2)
      -- FIXME: expand...
      _                                        -> Nothing

-- | Asks the question: is the operator within the given category?
isOp :: BinOpKind -> F.BinaryOp -> Bool
isOp cat = (== cat) . binOpKind

data BinOpKind = AddOp | MulOp | DivOp | PowerOp | LogicOp | RelOp deriving Eq
binOpKind :: F.BinaryOp -> BinOpKind
binOpKind F.Addition         = AddOp
binOpKind F.Subtraction      = AddOp
binOpKind F.Multiplication   = MulOp
binOpKind F.Division         = DivOp
binOpKind F.Exponentiation   = PowerOp
binOpKind F.Concatenation    = AddOp
binOpKind F.GT               = RelOp
binOpKind F.GTE              = RelOp
binOpKind F.LT               = RelOp
binOpKind F.LTE              = RelOp
binOpKind F.EQ               = RelOp
binOpKind F.NE               = RelOp
binOpKind F.Or               = LogicOp
binOpKind F.And              = LogicOp
binOpKind F.XOr              = LogicOp
binOpKind F.Equivalent       = RelOp
binOpKind F.NotEquivalent    = RelOp
binOpKind (F.BinCustom _)    = RelOp

-- | Get information about imported variables coming from mod files.
getImportedVariables :: UnitSolver (M.Map VV UnitInfo)
getImportedVariables = do
  pf <- getProgramFile
  nmap <- getNameParamMap
  -- Translate a Use AST node into a pair mapping unique name to 'local' source name in this program file.
  let useToPair (F.UseID _ _ e) = (varName e, srcName e)
      useToPair (F.UseRename _ _ e1 _) = (varName e1, srcName e1) -- (unique name, 'local' source name)
  -- A map of modules -> (maps of variables -> their unit info).
  let modnmaps = [ M.fromList (mapMaybe f (M.toList npkmap))
                 -- find all StUse statements and identify variables that need to be imported from nmap
                 | F.StUse _ _ e _ only alist <- universeBi pf :: [ F.Statement UA ]
                 , let mod = srcName e
                 , let uses = map useToPair (fromMaybe [] (F.aStrip <$> alist))
                 , Just npkmap <- [M.lookup (F.Named mod) nmap]
                 , let f (npk, ui) = case npk of
                         (NPKVariable (var, src))
                           -- import all variables from module -- apply any renames from uses
                           | only == F.Permissive -> Just (NPKVariable (var, src `fromMaybe` lookup var uses), ui)
                           -- only import variable mentioned in uses
                           | Just src' <- lookup var uses -> Just (NPKVariable (var, src'), ui)
                         _ -> Nothing
                 ]
  pure $ M.fromList [ (vv, foldUnits units) | (NPKVariable vv, units) <- M.toList (M.unions modnmaps) ]

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

logDebugNoOrigin :: Text -> UnitSolver ()
logDebugNoOrigin msg = do
  pf <- gets usProgramFile
  logDebug' pf msg

dumpConsM :: String -> Constraints -> UnitSolver ()
dumpConsM str = logDebugNoOrigin . describe . unlines . ([replicate 50 '-', str ++ ":"]++) . (++[replicate 50 '^']) . map f
  where
    f (ConEq u1 u2)  = show (flattenUnits u1) ++ " === " ++ show (flattenUnits u2)
    f (ConConj cons) = intercalate " && " (map f cons)

debugLogging :: UnitSolver ()
debugLogging = do
    (logDebugNoOrigin . describe . unlines . map (\ (ConEq u1 u2) -> "  ***AbsConstraint: " ++ show (flattenUnits u1) ++ " === " ++ show (flattenUnits u2))) =<< extractConstraints
    pf   <- getProgramFile
    cons <- getConstraints
    vum  <- getVarUnitMap
    logDebugNoOrigin . describe . unlines $ [ "  " ++ show info ++ " :: " ++ n | ((n, _), info) <- M.toList vum ]
    logDebugNoOrigin ""
    uam <- getUnitAliasMap
    logDebugNoOrigin . describe . unlines $ [ "  " ++ n ++ " = " ++ show info | (n, info) <- M.toList uam ]
    logDebugNoOrigin . describe . unlines $ map (\ (ConEq u1 u2) -> "  ***Constraint: " ++ show (flattenUnits u1) ++ " === " ++ show (flattenUnits u2)) cons
    logDebugNoOrigin $ describeShow cons <> "\n"
    forM_ (universeBi pf) $ \ pu -> case pu of
      F.PUFunction {}
        | Just (ConConj con) <- UA.getConstraint pu ->
          logDebugNoOrigin . describe . unlines $ (puName pu ++ ":"):map (\ (ConEq u1 u2) -> "    constraint: " ++ show (flattenUnits u1) ++ " === " ++ show (flattenUnits u2)) con
      F.PUSubroutine {}
        | Just (ConConj con) <- UA.getConstraint pu ->
          logDebugNoOrigin . describe . unlines $ (puName pu ++ ":"):map (\ (ConEq u1 u2) -> "    constraint: " ++ show (flattenUnits u1) ++ " === " ++ show (flattenUnits u2)) con
      _ -> pure ()
    let (lhsM, rhsM, _, lhsColA, rhsColA) = constraintsToMatrices cons
    logDebugNoOrigin "--------------------------------------------------\nLHS Cols:"
    logDebugNoOrigin $ describeShow lhsColA
    logDebugNoOrigin "--------------------------------------------------\nRHS Cols:"
    logDebugNoOrigin $ describeShow rhsColA
    logDebugNoOrigin "--------------------------------------------------\nLHS M:"
    logDebugNoOrigin $ describeShow lhsM
    logDebugNoOrigin "--------------------------------------------------\nRHS M:"
    logDebugNoOrigin $ describeShow rhsM
    logDebugNoOrigin "--------------------------------------------------\nAUG M:"
    let augM = if H.rows rhsM == 0 || H.cols rhsM == 0 then lhsM else H.fromBlocks [[lhsM, rhsM]]
    logDebugNoOrigin $ describeShow augM
    logDebugNoOrigin "--------------------------------------------------\nSolved (hnf) M:"
    let hnfM = Flint.hnf augM
    logDebugNoOrigin $ describeShow hnfM
    logDebugNoOrigin "--------------------------------------------------\nSolved (normHNF) M:"
    let (solvedM, newColIndices) = Flint.normHNF augM
    logDebugNoOrigin . describeShow $ solvedM
    logDebugNoOrigin $ "newColIndices = " <> describeShow newColIndices
    logDebugNoOrigin "--------------------------------------------------\nLHS Cols with newColIndices:"
    let lhsCols = A.elems lhsColA ++ map (lhsColA A.!) newColIndices
    logDebugNoOrigin $ describe . unlines . map show $ zip [(0::Int)..] lhsCols
    -- logDebugNoOrigin "--------------------------------------------------\nSolved (SVD) M:"
    -- logDebugNoOrigin $ show (H.linearSolveSVD lhsM rhsM)
    -- logDebugNoOrigin "--------------------------------------------------\nSingular Values:"
    -- logDebugNoOrigin $ show (H.singularValues lhsM)
    logDebugNoOrigin "--------------------------------------------------"
    logDebugNoOrigin $ "Rank LHS: " <> describeShow (H.rank lhsM)
    logDebugNoOrigin "--------------------------------------------------"
    let augA = if H.rows rhsM == 0 || H.cols rhsM == 0 then lhsM else H.fromBlocks [[lhsM, rhsM]]
    logDebugNoOrigin $ "Rank Augmented: " <> describeShow (H.rank augA)
    logDebugNoOrigin "--------------------------------------------------\nGenUnitAssignments:"
    let unitAssignments = genUnitAssignments cons
    logDebugNoOrigin . describe . unlines $ map (\ (u1s, u2) -> "  ***UnitAssignment: " ++ show u1s ++ " === " ++ show (flattenUnits u2) ++ "\n") unitAssignments
    logDebugNoOrigin "--------------------------------------------------"
    let unitAssignmentsSBV = BackendSBV.genUnitAssignments cons
    logDebugNoOrigin . describe . unlines $ map (\ (u1s, u2) -> "  ***UnitAssignmentSBV: " ++ show u1s ++ " === " ++ show (flattenUnits u2)) unitAssignmentsSBV
    logDebugNoOrigin "--------------------------------------------------\nProvenance:"
    let (augM', p) = provenance augM
    logDebugNoOrigin . describeShow $ augM'
    logDebugNoOrigin . describeShow $ p

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

-- convenience
puName :: F.ProgramUnit UA -> F.Name
puName pu
  | F.Named n <- FA.puName pu = n
  | otherwise                 = "_nameless"

puSrcName :: F.ProgramUnit UA -> F.Name
puSrcName pu
  | F.Named n <- FA.puSrcName pu = n
  | otherwise                    = "_nameless"

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

-- | Intrinics that take arbitrary number of arguments. Entry in table
-- 'intrinsicUnits' will contain a single item in the argument list,
-- corresponding to the template used for all arguments.
specialCaseArbitraryArgs :: S.Set F.Name
specialCaseArbitraryArgs = S.fromList [ "max", "max0", "amax1", "dmax1", "amax0", "max1"
                                      , "min", "min0", "amin1", "dmin1", "amin0", "min1" ]

-- | Intrinsics table: name => (return-unit, parameter-units). See also 'specialCaseArbitraryArgs'.
intrinsicUnits :: M.Map F.Name (UnitInfo, [UnitInfo])
intrinsicUnits =
  M.fromList
    [ ("transfer", (UnitParamEAPAbs ("'b", "'b"), [UnitParamEAPAbs ("'a", "'a"), UnitParamEAPAbs ("'b", "'b")]))
    , ("abs", (UnitParamEAPAbs ("'a", "'a"), [UnitParamEAPAbs ("'a", "'a")]))
    , ("iabs", (UnitParamEAPAbs ("'a", "'a"), [UnitParamEAPAbs ("'a", "'a")]))
    , ("dabs", (UnitParamEAPAbs ("'a", "'a"), [UnitParamEAPAbs ("'a", "'a")]))
    , ("cabs", (UnitParamEAPAbs ("'a", "'a"), [UnitParamEAPAbs ("'a", "'a")]))
    , ("aimag", (UnitParamEAPAbs ("'a", "'a"), [UnitParamEAPAbs ("'a", "'a")]))
    , ("aint", (UnitParamEAPAbs ("'a", "'a"), [UnitParamEAPAbs ("'a", "'a")]))
    , ("dint", (UnitParamEAPAbs ("'a", "'a"), [UnitParamEAPAbs ("'a", "'a")]))
    , ("anint", (UnitParamEAPAbs ("'a", "'a"), [UnitParamEAPAbs ("'a", "'a")]))
    , ("dnint", (UnitParamEAPAbs ("'a", "'a"), [UnitParamEAPAbs ("'a", "'a")]))
    , ("cmplx", (UnitParamEAPAbs ("'a", "'a"), [UnitParamEAPAbs ("'a", "'a")]))
    , ("conjg", (UnitParamEAPAbs ("'a", "'a"), [UnitParamEAPAbs ("'a", "'a")]))
    , ("dble", (UnitParamEAPAbs ("'a", "'a"), [UnitParamEAPAbs ("'a", "'a")]))
    , ("dim", (UnitParamEAPAbs ("'a", "'a"), [UnitParamEAPAbs ("'a", "'a"), UnitParamEAPAbs ("'a", "'a")]))
    , ("idim", (UnitParamEAPAbs ("'a", "'a"), [UnitParamEAPAbs ("'a", "'a"), UnitParamEAPAbs ("'a", "'a")]))
    , ("ddim", (UnitParamEAPAbs ("'a", "'a"), [UnitParamEAPAbs ("'a", "'a"), UnitParamEAPAbs ("'a", "'a")]))
    , ("dprod", (UnitParamEAPAbs ("'a", "'a"), [UnitParamEAPAbs ("'a", "'a")]))
    , ("ceiling", (UnitParamEAPAbs ("'a", "'a"), [UnitParamEAPAbs ("'a", "'a")]))
    , ("floor", (UnitParamEAPAbs ("'a", "'a"), [UnitParamEAPAbs ("'a", "'a")]))
    , ("int", (UnitParamEAPAbs ("'a", "'a"), [UnitParamEAPAbs ("'a", "'a")]))
    , ("ifix", (UnitParamEAPAbs ("'a", "'a"), [UnitParamEAPAbs ("'a", "'a")]))
    , ("idint", (UnitParamEAPAbs ("'a", "'a"), [UnitParamEAPAbs ("'a", "'a")]))
    , ("maxval", (UnitParamEAPAbs ("'a", "'a"), [UnitParamEAPAbs ("'a", "'a")]))
    , ("minval", (UnitParamEAPAbs ("'a", "'a"), [UnitParamEAPAbs ("'a", "'a")]))
    , ("max", (UnitParamEAPAbs ("'a", "'a"), [UnitParamEAPAbs ("'a", "'a")]))   -- special case: arbitrary # of parameters
    , ("min", (UnitParamEAPAbs ("'a", "'a"), [UnitParamEAPAbs ("'a", "'a")]))   -- special case: arbitrary # of parameters
    , ("min0", (UnitParamEAPAbs ("'a", "'a"), [UnitParamEAPAbs ("'a", "'a")]))  -- special case: arbitrary # of parameters
    , ("amin1", (UnitParamEAPAbs ("'a", "'a"), [UnitParamEAPAbs ("'a", "'a")])) -- special case: arbitrary # of parameters
    , ("dmin1", (UnitParamEAPAbs ("'a", "'a"), [UnitParamEAPAbs ("'a", "'a")])) -- special case: arbitrary # of parameters
    , ("amin0", (UnitParamEAPAbs ("'a", "'a"), [UnitParamEAPAbs ("'a", "'a")])) -- special case: arbitrary # of parameters
    , ("min1", (UnitParamEAPAbs ("'a", "'a"), [UnitParamEAPAbs ("'a", "'a")]))  -- special case: arbitrary # of parameters
    , ("mod", (UnitParamEAPAbs ("'a", "'a"), [UnitParamEAPAbs ("'a", "'a"), UnitParamEAPAbs ("'b", "'b")]))
    , ("modulo", (UnitParamEAPAbs ("'a", "'a"), [UnitParamEAPAbs ("'a", "'a"), UnitParamEAPAbs ("'b", "'b")]))
    , ("amod", (UnitParamEAPAbs ("'a", "'a"), [UnitParamEAPAbs ("'a", "'a"), UnitParamEAPAbs ("'b", "'b")]))
    , ("dmod", (UnitParamEAPAbs ("'a", "'a"), [UnitParamEAPAbs ("'a", "'a"), UnitParamEAPAbs ("'b", "'b")]))
    , ("nint", (UnitParamEAPAbs ("'a", "'a"), [UnitParamEAPAbs ("'a", "'a")]))
    , ("real", (UnitParamEAPAbs ("'a", "'a"), [UnitParamEAPAbs ("'a", "'a")]))
    , ("float", (UnitParamEAPAbs ("'a", "'a"), [UnitParamEAPAbs ("'a", "'a")]))
    , ("sngl", (UnitParamEAPAbs ("'a", "'a"), [UnitParamEAPAbs ("'a", "'a")]))
    , ("sign", (UnitParamEAPAbs ("'a", "'a"), [UnitParamEAPAbs ("'a", "'a"), UnitParamEAPAbs ("'b", "'b")]))
    , ("isign", (UnitParamEAPAbs ("'a", "'a"), [UnitParamEAPAbs ("'a", "'a"), UnitParamEAPAbs ("'b", "'b")]))
    , ("dsign", (UnitParamEAPAbs ("'a", "'a"), [UnitParamEAPAbs ("'a", "'a"), UnitParamEAPAbs ("'b", "'b")]))
    , ("present", (UnitParamEAPAbs ("'a", "'a"), [UnitlessVar]))
    , ("sqrt", (UnitParamEAPAbs ("'a", "'a"), [UnitPow (UnitParamEAPAbs ("'a", "'a")) 2]))
    , ("dsqrt", (UnitParamEAPAbs ("'a", "'a"), [UnitPow (UnitParamEAPAbs ("'a", "'a")) 2]))
    , ("csqrt", (UnitParamEAPAbs ("'a", "'a"), [UnitPow (UnitParamEAPAbs ("'a", "'a")) 2]))
    , ("exp", (UnitlessVar, [UnitlessVar]))
    , ("dexp", (UnitlessVar, [UnitlessVar]))
    , ("cexp", (UnitlessVar, [UnitlessVar]))
    , ("alog", (UnitlessVar, [UnitlessVar]))
    , ("dlog", (UnitlessVar, [UnitlessVar]))
    , ("clog", (UnitlessVar, [UnitlessVar]))
    , ("alog10", (UnitlessVar, [UnitlessVar]))
    , ("dlog10", (UnitlessVar, [UnitlessVar]))
    , ("sin", (UnitlessVar, [UnitlessVar]))
    , ("dsin", (UnitlessVar, [UnitlessVar]))
    , ("csin", (UnitlessVar, [UnitlessVar]))
    , ("cos", (UnitlessVar, [UnitlessVar]))
    , ("dcos", (UnitlessVar, [UnitlessVar]))
    , ("ccos", (UnitlessVar, [UnitlessVar]))
    , ("tan", (UnitlessVar, [UnitlessVar]))
    , ("dtan", (UnitlessVar, [UnitlessVar]))
    , ("asin", (UnitlessVar, [UnitlessVar]))
    , ("dasin", (UnitlessVar, [UnitlessVar]))
    , ("acos", (UnitlessVar, [UnitlessVar]))
    , ("dacos", (UnitlessVar, [UnitlessVar]))
    , ("atan", (UnitlessVar, [UnitlessVar]))
    , ("datan", (UnitlessVar, [UnitlessVar]))
    , ("atan2", (UnitlessVar, [UnitParamEAPAbs ("'a", "'a"), UnitParamEAPAbs ("'a", "'a")]))
    , ("datan2", (UnitlessVar, [UnitParamEAPAbs ("'a", "'a"), UnitParamEAPAbs ("'a", "'a")]))
    , ("sinh", (UnitlessVar, [UnitlessVar]))
    , ("dsinh", (UnitlessVar, [UnitlessVar]))
    , ("cosh", (UnitlessVar, [UnitlessVar]))
    , ("dcosh", (UnitlessVar, [UnitlessVar]))
    , ("tanh", (UnitlessVar, [UnitlessVar]))
    , ("dtanh", (UnitlessVar, [UnitlessVar]))
    , ("iand", (UnitParamEAPAbs ("'a", "'a"), [UnitParamEAPAbs ("'a", "'a"), UnitParamEAPAbs ("'a", "'a")]))
    ]

-- Others: reshape, merge need special handling

-- | Compile a program to a 'ModFile' containing units information.
compileUnits :: UnitOpts -> ModFiles -> F.ProgramFile Annotation -> IO ModFile
compileUnits uo mfs pf = do
  let (pf', _, _) = withCombinedEnvironment mfs . fmap UA.mkUnitAnnotation $ pf

  let analysis = runReaderT (runInference runCompileUnits) $
        UnitEnv
        { unitOpts = uo
        , unitProgramFile = pf
        }

  report <- runAnalysisT (F.pfGetFilename pf) (logOutputNone True) LogError mfs analysis

  case report ^? arResult . _ARSuccess . _1 of
    Just cu -> return (genUnitsModFile pf' cu)
    Nothing -> fail "compileUnits: units analysis failed"

max_stars_repo_name: apthorpe/camfort

max_issues_repo_name: apthorpe/camfort

max_forks_repo_name: apthorpe/camfort


text:
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE DeriveTraversable #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE PatternSynonyms #-}
{-# LANGUAGE StandaloneDeriving #-}
{-# LANGUAGE TypeFamilies #-}
-- |
-- Module : Graphics.Color.Model.LCH
module Graphics.Color.Model.LCH
  ( LCH
  -- * Constructors for an LCH color model.
  , pattern ColorLCH
  , pattern ColorLCHA
  , Color
  , ColorModel(..)
  , lch2lxy
  , lxy2lch
  ) where

import Data.Complex ( Complex(..), polar, mkPolar )
import Data.Fixed ( mod' )
import Foreign.Storable
import Graphics.Color.Model.Internal

-----------
--- LCH ---
-----------

-- | CIEL*C*H color model, representing a cylindrical reparameterization
--   of CIEL*a*b* or CIEL*u*v*.
data LCH

-- | `LCH` color model
newtype instance Color LCH e = LCH (V3 e)

-- | Constructor for @LCH@.
pattern ColorLCH :: e -> e -> e -> Color LCH e
pattern ColorLCH l c h = LCH (V3 l c h)
{-# COMPLETE ColorLCH #-}

-- | Constructor for @LCH@ with alpha channel.
pattern ColorLCHA :: e -> e -> e -> e -> Color (Alpha LCH) e
pattern ColorLCHA l c h a = Alpha (ColorLCH l c h) a
{-# COMPLETE ColorLCHA #-}

-- | `LCH` color model
deriving instance Eq e => Eq (Color LCH e)
-- | `LCH` color model
deriving instance Ord e => Ord (Color LCH e)
-- | `LCH` color model
deriving instance Functor (Color LCH)
-- | `LCH` color model
deriving instance Applicative (Color LCH)
-- | `LCH` color model
deriving instance Foldable (Color LCH)
-- | `LCH` color model
deriving instance Traversable (Color LCH)
-- | `LCH` color model
deriving instance Storable e => Storable (Color LCH e)

-- | `LCH` color model
instance Elevator e => Show (Color LCH e) where
  showsPrec _ = showsColorModel

-- | `LCH` color model
instance Elevator e => ColorModel LCH e where
  type Components LCH e = (e, e, e)
  toComponents (ColorLCH l c h) = (l, c, h)
  {-# INLINE toComponents #-}
  fromComponents (l, c, h) = ColorLCH l c h
  {-# INLINE fromComponents #-}

lch2lxy :: Color LCH Double -> Components LCH Double
lch2lxy (ColorLCH l c h) = (l, x, y)
 where
  !h' = h * pi / 180
  (x :+ y) = mkPolar c h'
{-# INLINE lch2lxy #-}

lxy2lch :: Components LCH Double -> Color LCH Double
lxy2lch (l, x, y) = ColorLCH l c h
 where
  (c,h') = polar (x :+ y)
  !h = (h' * 180 / pi) `mod'` 360
{-# INLINE lxy2lch #-}

ext: hs

max_stars_repo_name: TravisCardwell/Color

max_issues_repo_name: TravisCardwell/Color

max_forks_repo_name: TravisCardwell/Color


text:
{-# LANGUAGE CPP #-}
-- |
-- Module      : Test.LeanCheck.Basic
-- Copyright   : (c) 2015-2020 Rudy Matela
-- License     : 3-Clause BSD  (see the file LICENSE)
-- Maintainer  : Rudy Matela <rudy@matela.com.br>
--
-- This module is part of LeanCheck,
-- a simple enumerative property-based testing library.
--
-- This module exports "Test.LeanCheck.Core" along with:
--
--   * support for 'Listable' 6-tuples up to 12-tuples;
--   * 'tiers' constructors (@consN@) with arities from 6 up to 12;
--   * a 'Listable' 'Word' instance;
--   * a 'Listable' 'Ratio' instance (consequently 'Listable' 'Rational');
--   * a 'Listable' 'Complex' instance;
--   * 'Listable' 'Int8/16/32/64' instances;
--   * 'Listable' 'Word8/16/32/64' instances;
--   * 'Listable' instances for "Foreign.C" types;
--   * a 'Listable' 'ExitCode' instance;
--   * a 'Listable' 'GeneralCategory' instance;
--   * 'Listable' 'Buffer/IO/SeekMode' instances;
--   * the operators 'addWeight' and 'ofWeight'.
--
-- The above includes all types defined in the Haskell 2010 Report
-- with the exception of Array, IO, Handle, HandlePosn, IOErrorType.
--
-- "Test.LeanCheck" already exports everything from this module.
-- You are probably better off importing it.
--
-- You should /only/ import "Test.LeanCheck.Basic"
-- if you /only/ want the above basic functionality.
module Test.LeanCheck.Basic
  ( module Test.LeanCheck.Core

  , cons6
  , cons7
  , cons8
  , cons9
  , cons10
  , cons11
  , cons12

  , ofWeight
  , addWeight
  )
where

import Test.LeanCheck.Core
import Data.Ratio
import Data.Complex
import Data.Int
import Data.Word
import Data.Char (GeneralCategory)
import System.IO (IOMode (..), BufferMode (..), SeekMode (..))
import Foreign.C
import System.Exit

instance (Listable a, Listable b, Listable c,
          Listable d, Listable e, Listable f) =>
         Listable (a,b,c,d,e,f) where
  tiers  =  productWith (\x (y,z,w,v,u) -> (x,y,z,w,v,u)) tiers tiers

instance (Listable a, Listable b, Listable c, Listable d,
          Listable e, Listable f, Listable g) =>
         Listable (a,b,c,d,e,f,g) where
  tiers  =  productWith (\x (y,z,w,v,u,r) -> (x,y,z,w,v,u,r)) tiers tiers

instance (Listable a, Listable b, Listable c, Listable d,
          Listable e, Listable f, Listable g, Listable h) =>
         Listable (a,b,c,d,e,f,g,h) where
  tiers  =  productWith (\x (y,z,w,v,u,r,s) -> (x,y,z,w,v,u,r,s)) tiers tiers

instance (Listable a, Listable b, Listable c, Listable d, Listable e,
          Listable f, Listable g, Listable h, Listable i) =>
         Listable (a,b,c,d,e,f,g,h,i) where
  tiers  =  productWith (\x (y,z,w,v,u,r,s,t) -> (x,y,z,w,v,u,r,s,t))
                        tiers tiers

instance (Listable a, Listable b, Listable c, Listable d, Listable e,
          Listable f, Listable g, Listable h, Listable i, Listable j) =>
         Listable (a,b,c,d,e,f,g,h,i,j) where
  tiers  =  productWith (\x (y,z,w,v,u,r,s,t,o) -> (x,y,z,w,v,u,r,s,t,o))
                        tiers tiers

instance (Listable a, Listable b, Listable c, Listable d,
          Listable e, Listable f, Listable g, Listable h,
          Listable i, Listable j, Listable k) =>
         Listable (a,b,c,d,e,f,g,h,i,j,k) where
  tiers  =  productWith (\x (y,z,w,v,u,r,s,t,o,p) -> (x,y,z,w,v,u,r,s,t,o,p))
                        tiers tiers

instance (Listable a, Listable b, Listable c, Listable d,
          Listable e, Listable f, Listable g, Listable h,
          Listable i, Listable j, Listable k, Listable l) =>
         Listable (a,b,c,d,e,f,g,h,i,j,k,l) where
  tiers  =  productWith (\x (y,z,w,v,u,r,s,t,o,p,q) ->
                          (x,y,z,w,v,u,r,s,t,o,p,q))
                        tiers tiers

-- | Returns tiers of applications of a 6-argument constructor.
cons6 :: (Listable a, Listable b, Listable c, Listable d, Listable e, Listable f)
      => (a -> b -> c -> d -> e -> f -> g) -> [[g]]
cons6 f  =  delay $ mapT (uncurry6 f) tiers

-- | Returns tiers of applications of a 7-argument constructor.
cons7 :: (Listable a, Listable b, Listable c, Listable d,
          Listable e, Listable f, Listable g)
      => (a -> b -> c -> d -> e -> f -> g -> h) -> [[h]]
cons7 f  =  delay $ mapT (uncurry7 f) tiers

-- | Returns tiers of applications of a 8-argument constructor.
cons8 :: (Listable a, Listable b, Listable c, Listable d,
          Listable e, Listable f, Listable g, Listable h)
      => (a -> b -> c -> d -> e -> f -> g -> h -> i) -> [[i]]
cons8 f  =  delay $ mapT (uncurry8 f) tiers

-- | Returns tiers of applications of a 9-argument constructor.
cons9 :: (Listable a, Listable b, Listable c, Listable d, Listable e,
          Listable f, Listable g, Listable h, Listable i)
      => (a -> b -> c -> d -> e -> f -> g -> h -> i -> j) -> [[j]]
cons9 f  =  delay $ mapT (uncurry9 f) tiers

-- | Returns tiers of applications of a 10-argument constructor.
cons10 :: (Listable a, Listable b, Listable c, Listable d, Listable e,
           Listable f, Listable g, Listable h, Listable i, Listable j)
       => (a -> b -> c -> d -> e -> f -> g -> h -> i -> j -> k) -> [[k]]
cons10 f  =  delay $ mapT (uncurry10 f) tiers

-- | Returns tiers of applications of a 11-argument constructor.
cons11 :: (Listable a, Listable b, Listable c, Listable d,
           Listable e, Listable f, Listable g, Listable h,
           Listable i, Listable j, Listable k)
       => (a -> b -> c -> d -> e -> f -> g -> h -> i -> j -> k -> l) -> [[l]]
cons11 f  =  delay $ mapT (uncurry11 f) tiers

-- | Returns tiers of applications of a 12-argument constructor.
cons12 :: (Listable a, Listable b, Listable c, Listable d,
           Listable e, Listable f, Listable g, Listable h,
           Listable i, Listable j, Listable k, Listable l)
       => (a->b->c->d->e->f->g->h->i->j->k->l->m) -> [[m]]
cons12 f  =  delay $ mapT (uncurry12 f) tiers

uncurry6 :: (a->b->c->d->e->f->g) -> (a,b,c,d,e,f) -> g
uncurry6 f (x,y,z,w,v,u)  =  f x y z w v u

uncurry7 :: (a->b->c->d->e->f->g->h) -> (a,b,c,d,e,f,g) -> h
uncurry7 f (x,y,z,w,v,u,r)  =  f x y z w v u r

uncurry8 :: (a->b->c->d->e->f->g->h->i) -> (a,b,c,d,e,f,g,h) -> i
uncurry8 f (x,y,z,w,v,u,r,s)  =  f x y z w v u r s

uncurry9 :: (a->b->c->d->e->f->g->h->i->j) -> (a,b,c,d,e,f,g,h,i) -> j
uncurry9 f (x,y,z,w,v,u,r,s,t)  =  f x y z w v u r s t

uncurry10 :: (a->b->c->d->e->f->g->h->i->j->k) -> (a,b,c,d,e,f,g,h,i,j) -> k
uncurry10 f (x,y,z,w,v,u,r,s,t,o)  =  f x y z w v u r s t o

uncurry11 :: (a->b->c->d->e->f->g->h->i->j->k->l)
          -> (a,b,c,d,e,f,g,h,i,j,k) -> l
uncurry11 f (x,y,z,w,v,u,r,s,t,o,p)  =  f x y z w v u r s t o p

uncurry12 :: (a->b->c->d->e->f->g->h->i->j->k->l->m)
          -> (a,b,c,d,e,f,g,h,i,j,k,l) -> m
uncurry12 f (x,y,z,w,v,u,r,s,t,o,p,q)  =  f x y z w v u r s t o p q

-- | > list :: [Rational]  =
--   >   [   0  % 1
--   >   ,   1  % 1
--   >   , (-1) % 1
--   >   ,   1  % 2,   2  % 1
--   >   , (-1) % 2, (-2) % 1
--   >   ,   1  % 3,   3  % 1
--   >   , (-1) % 3, (-3) % 1
--   >   ,   1  % 4,   2  % 3,   3  % 2,   4  % 1
--   >   , (-1) % 4, (-2) % 3, (-3) % 2, (-4) % 1
--   >   ,   1  % 5,   5  % 1
--   >   , (-1) % 5, (-5) % 1
--   >   , ...
--   >   ]
instance (Integral a, Listable a) => Listable (Ratio a) where
  tiers  =  mapT (uncurry (%)) . reset
        $ tiers `suchThat` (\(n,d) -> d > 0 && n `gcd` d == 1)

instance (RealFloat a, Listable a) => Listable (Complex a) where
  tiers  =  cons2 (:+)

-- | > list :: [Word]  =  [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, ...]
instance Listable Word where
  list  =  listIntegral

-- | > list :: [Word8]  =  [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, ..., 255]
instance Listable Word8 where
  list  =  listIntegral

-- | > list :: [Word16]  =  [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, ..., 65535]
instance Listable Word16 where
  list  =  listIntegral

-- | > list :: [Word32]  =  [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, ...]
instance Listable Word32 where
  list  =  listIntegral

-- | > list :: [Word64]  =  [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, ...]
instance Listable Word64 where
  list  =  listIntegral

-- | > list :: [Int8]  =  [0, 1, -1, 2, -2, 3, -3, ..., 127, -127, -128]
instance Listable Int8 where
  list  =  listIntegral

-- | > list :: [Int16]  =  [0, 1, -1, 2, -2, ..., 32767, -32767, -32768]
instance Listable Int16 where
  list  =  listIntegral

-- | > list :: [Int32]  =  [0, 1, -1, 2, -2, 3, -3, 4, -4, 5, -5, 6, ...]
instance Listable Int32 where
  list  =  listIntegral

-- | > list :: [Int64]  =  [0, 1, -1, 2, -2, 3, -3, 4, -4, 5, -5, 6, ...]
instance Listable Int64 where
  list  =  listIntegral

instance Listable CChar      where  list  =  listIntegral
instance Listable CSChar     where  list  =  listIntegral
instance Listable CUChar     where  list  =  listIntegral
instance Listable CShort     where  list  =  listIntegral
instance Listable CUShort    where  list  =  listIntegral
instance Listable CInt       where  list  =  listIntegral
instance Listable CUInt      where  list  =  listIntegral
instance Listable CLong      where  list  =  listIntegral
instance Listable CULong     where  list  =  listIntegral
instance Listable CPtrdiff   where  list  =  listIntegral
instance Listable CSize      where  list  =  listIntegral
instance Listable CWchar     where  list  =  listIntegral
instance Listable CSigAtomic where  list  =  listIntegral
instance Listable CLLong     where  list  =  listIntegral
instance Listable CULLong    where  list  =  listIntegral
instance Listable CIntPtr    where  list  =  listIntegral
instance Listable CUIntPtr   where  list  =  listIntegral
instance Listable CIntMax    where  list  =  listIntegral
instance Listable CUIntMax   where  list  =  listIntegral
instance Listable CClock     where  list  =  listIntegral
instance Listable CTime      where  list  =  listIntegral
instance Listable CFloat     where  tiers  =  tiersFloating
instance Listable CDouble    where  tiers  =  tiersFloating
#if __GLASGOW_HASKELL__ >= 802
instance Listable CBool      where  list  =  listIntegral
#endif
#if __GLASGOW_HASKELL__
instance Listable CUSeconds  where  list  =  listIntegral
instance Listable CSUSeconds where  list  =  listIntegral
#endif

-- | Only includes valid POSIX exit codes
--
-- > > list :: [ExitCode]
-- > [ExitSuccess, ExitFailure 1, ExitFailure 2, ..., ExitFailure 255]
instance Listable ExitCode where  list  =  ExitSuccess : map ExitFailure [1..255]

instance Listable GeneralCategory where  list  =  [minBound..maxBound]

instance Listable IOMode where
  tiers  =  cons0 ReadMode
         \/ cons0 WriteMode
         \/ cons0 AppendMode
         \/ cons0 ReadWriteMode

instance Listable BufferMode where
  tiers  =  cons0 NoBuffering
         \/ cons0 LineBuffering
         \/ cons1 BlockBuffering

instance Listable SeekMode where
  tiers  =  cons0 AbsoluteSeek
         \/ cons0 RelativeSeek
         \/ cons0 SeekFromEnd

-- | Resets the weight of a constructor or tiers.
--
-- > > [ [], [], ..., xs, ys, zs, ... ] `ofWeight` 1
-- > [ [], xs, ys, zs, ... ]
--
-- > > [ xs, ys, zs, ... ] `ofWeight` 2
-- > [ [], [], xs, ys, zs, ... ]
--
-- > > [ [], xs, ys, zs, ... ] `ofWeight` 3
-- > [ [], [], [], xs, ys, zs, ... ]
--
-- Typically used as an infix operator when defining 'Listable' instances:
--
-- > instance Listable <Type> where
-- >   tiers  =  ...
-- >          \/ cons<N> <Cons>  `ofWeight`  <W>
-- >          \/ ...
--
-- /Warning:/ do not apply @ \`ofWeight\` 0 @ to recursive data structure
-- constructors.  In general this will make the list of size 0 infinite,
-- breaking the tier invariant (each tier must be finite).
--
-- @ \`ofWeight\` /n/ @ is equivalent to 'reset' followed
-- by @/n/@ applications of 'delay'.
ofWeight :: [[a]] -> Int -> [[a]]
ofWeight xss w  =  dropWhile null xss `addWeight` w

-- | Adds to the weight of a constructor or tiers.
--
-- > instance Listable <Type> where
-- >   tiers  =  ...
-- >          \/ cons<N> <Cons>  `addWeight`  <W>
-- >          \/ ...
--
-- Typically used as an infix operator when defining 'Listable' instances:
--
-- > > [ xs, ys, zs, ... ] `addWeight` 1
-- > [ [], xs, ys, zs, ... ]
--
-- > > [ xs, ys, zs, ... ] `addWeight` 2
-- > [ [], [], xs, ys, zs, ... ]
--
-- > > [ [], xs, ys, zs, ... ] `addWeight` 3
-- > [ [], [], [], [], xs, ys, zs, ... ]
--
-- @ \`addWeight\` /n/ @ is equivalent to @/n/@ applications of 'delay'.
addWeight :: [[a]] -> Int -> [[a]]
addWeight xss w  =  replicate w [] ++ xss

max_stars_repo_name: rudymatela/llcheck

max_issues_repo_name: rudymatela/llcheck

max_forks_repo_name: rudymatela/llcheck


text:
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE Strict, StrictData #-}
{-# LANGUAGE DeriveGeneric, DeriveAnyClass #-}

------------------------------------------------
-- |
-- Module    :  ArcGraph.EnhancedState
-- Copyright :  (c) Jun Yoshida 2019
-- License   :  BSD3
--
-- Implementing manipulations on enhanced states
--
------------------------------------------------

module ArcGraph.EnhancedState where

import GHC.Generics (Generic)

import Control.DeepSeq
import Control.Parallel.Strategies

import Control.Applicative
import Control.Monad
import Control.Monad.ST

import Data.Bifunctor

import Data.STRef
import Data.Maybe as M
import Data.Foldable
import qualified Data.List as L
import qualified Data.Vector as V
import qualified Data.Vector.Mutable as MV
import Data.Map.Strict (Map)
import qualified Data.Map.Strict as Map
import Data.IntMap.Strict (IntMap)
import qualified Data.IntMap.Strict as IMap
import qualified Data.Set as Set
import qualified Data.BitArray as BA

import qualified Numeric.LinearAlgebra as LA

import ArcGraph
import ArcGraph.Component
import ArcGraph.State

import Numeric.Matrix.Integral
import Numeric.Algebra.FreeModule as FM
import Numeric.Algebra.Frobenius as Frob
import Numeric.Algebra.Presentation
import Numeric.Algebra.Homology

{-- for debug
import Debug.Trace

traceCond :: Bool -> String -> a -> a
traceCond False _  = id
traceCond True msg = trace msg
--}

-----------------------
-- * Enhanced states --
-----------------------
class (DState ds, Ord e) => Enhancement ds e where
  listEnh :: (Alternative f) => Int -> ArcGraph -> ds -> f e
  diffEnh :: ArcGraph -> ds -> e -> FreeMod Int (ds,e)

listEStates :: (Enhancement ds e) => ArcGraph -> Int -> Int -> [(ds,e)]
listEStates ag i j = concatMap (\s -> (,) s <$> listEnh j ag s) $ listStates ag i

-----------------
-- * Instances --
-----------------
newtype MapEState pc = MEState (Map.Map pc SL2B)
  deriving (Eq, Ord, Show, Generic, NFData)

instance (DState ds, PComponent pc) => Enhancement ds (MapEState pc) where
  listEnh j ag st =
    let comps = getComponents (smoothing ag st)
        deg' = j - degree ag st + L.length comps
        subs = filter (\sub -> 2*L.length sub == deg') $ L.subsequences comps
        mapS sub = Map.fromList $! map (\c -> (c, if elem c sub then SLI else SLX)) comps
    in MEState <$> foldr' (\x xs -> pure (mapS x) <|> xs) empty subs

  diffEnh ag st (MEState mp) =
    FM.sumFM $ (diffState ag st :: V.Vector (Int,ds)) >>= \dstp -> do
      let (sign,dst) = dstp
      return $! ((,) dst . MEState . Map.fromList) FM.@$>% sign FM.@*% Frob.tqftZ (doIntersect ag) (Map.toList mp) (getComponents (smoothing ag dst))

----------------------
-- WORK IN PROGRESS --
----------------------{--
newtype BArray = BArray BA.BitArray
  deriving (Eq, Ord, Show, Generic)

instance NFData BArray where
  rnf (BArray x) = x `seq` ()

data BitEState pc = BEState pc BArray
  deriving (Eq, Show, Generic, NFData)
--}

data ArcGraphE ds e = AGraphE {
  arcGraph :: ArcGraph,
  state :: ds,
  enhancement :: e
  } deriving (Eq,Show,Generic,NFData)

enhancements :: (DState ds, PComponent pc, Alternative f) => Int -> ArcGraph -> ds -> f (Map.Map pc SL2B)
enhancements deg ag st =
  let comps = getComponents (smoothing ag st)
      deg' = deg + L.length comps
      subs = filter (\sub -> 2*L.length sub == deg') $ L.subsequences comps
      mapS sub = Map.fromList $! map (\c -> (c, if elem c sub then SLI else SLX)) comps
  in foldr' (\x xs -> pure (mapS x) <|> xs) empty subs

----------------------------------------------------------
-- The computation of (unnormalized) Khovanov homology
-----------------------------------------------------------
-- | The type to carry the data of Khovanov homologies
data KHData a ds e = KHData {
  subject :: ArcGraph,
  rank :: Int,
  tors :: [a],
  cycleV :: [FreeMod a (ds, e)],
  bndryV :: Maybe [FreeMod a (ds, e)] }
  deriving (Show,Eq,Generic, NFData)

vecToSum :: (Coefficient a, Ord b) => [b] -> Vec a -> FreeMod a b
vecToSum bs v = sumFM $! zipWith (@*@%) (vecToList v) bs

cohomologyToKH :: (Coefficient a, DState ds, Enhancement ds e) => ArcGraph -> [(ds, e)] -> Bool -> Homology a -> KHData a ds e
cohomologyToKH ag basis hasBndry hdata =
  let basisAGE = uncurry (AGraphE ag) <$> basis
  in KHData {
    subject = ag,
    rank = L.length (freeCycs hdata),
    tors = fmap fst (torsions hdata),
    cycleV = fmap (vecToSum basis) (freeCycs hdata ++ fmap snd (torsions hdata)),
    bndryV = if hasBndry
             then Just (fmap (vecToSum basis) (bndries hdata))
             else Nothing }

-- | Compute Khovanov homology for given range of cohomological degrees and a given quantum-degree
computeKhovanov :: (ChainEliminable a, NFData ds, NFData e, Show ds, Show e, Eq ds, Eq e, Enhancement ds e) => ArcGraph -> Int -> [ds] -> Bool -> IntMap (KHData a ds e)
computeKhovanov ag qdeg states hasBndry =
  let numCrs = countCross ag
      hdegs = [0..numCrs]
      slimAG = slimCross ag
      basis i = L.concatMap (\st -> (,) st <$> listEnh qdeg slimAG st) (filter ((==i). degree slimAG) states)
      !basisMap = force (IMap.fromAscList (map (\i -> (i,basis i)) hdegs))
      !diffs = flip (parMap rdeepseq) [0..numCrs-1] $ \i ->
        let sbasis = basisMap IMap.! i
            tbasis = basisMap IMap.! (i+1)
        in force $! present (uncurry (diffEnh slimAG)) sbasis tbasis
  in if numCrs == 0
     then -- The case where there is no crossing point;
       IMap.map (\v -> KHData slimAG (L.length v) [] (fmap (1@*@%) v) Nothing) basisMap
     else -- The case where there is at least one crossing point;
       let !hdata = force (homology diffs)
           !hdataMap = IMap.fromAscList $ filter (\hdti -> not (null (freeCycs (snd hdti)) && null (torsions (snd hdti))) ) $ zip hdegs hdata
       in flip IMap.mapWithKey hdataMap $ \i hdt -> cohomologyToKH slimAG (basisMap IMap.! i) hasBndry hdt

max_stars_repo_name: Junology/linkedgram

max_issues_repo_name: Junology/linkedgram

max_forks_repo_name: Junology/linkedgram


text:
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE DeriveDataTypeable, DeriveGeneric #-}
-----------------------------------------------------------------------------
-- |
-- Module      :  Statistics.Distribution.Beta
-- Copyright   :  (C) 2012 Edward Kmett,
-- License     :  BSD-style (see the file LICENSE)
--
-- Maintainer  :  Edward Kmett <ekmett@gmail.com>
-- Stability   :  provisional
-- Portability :  DeriveDataTypeable
--
----------------------------------------------------------------------------
module Statistics.Distribution.Beta
  ( BetaDistribution
    -- * Constructor
  , betaDistr
  , betaDistrE
  , improperBetaDistr
  , improperBetaDistrE
    -- * Accessors
  , bdAlpha
  , bdBeta
  ) where

import Control.Applicative
import Data.Aeson            (FromJSON(..), ToJSON, Value(..), (.:))
import Data.Binary           (Binary(..))
import Data.Data             (Data, Typeable)
import GHC.Generics          (Generic)
import Numeric.SpecFunctions (
  incompleteBeta, invIncompleteBeta, logBeta, digamma, log1p)
import Numeric.MathFunctions.Constants (m_NaN,m_neg_inf)
import qualified Statistics.Distribution as D
import Statistics.Internal

-- | The beta distribution
data BetaDistribution = BD
 { bdAlpha :: {-# UNPACK #-} !Double
   -- ^ Alpha shape parameter
 , bdBeta  :: {-# UNPACK #-} !Double
   -- ^ Beta shape parameter
 } deriving (Eq, Typeable, Data, Generic)

instance Show BetaDistribution where
  showsPrec n (BD a b) = defaultShow2 "improperBetaDistr" a b n
instance Read BetaDistribution where
  readPrec = defaultReadPrecM2 "improperBetaDistr" improperBetaDistrE

instance ToJSON BetaDistribution
instance FromJSON BetaDistribution where
  parseJSON (Object v) = do
    a <- v .: "bdAlpha"
    b <- v .: "bdBeta"
    maybe (fail $ errMsgI a b) return $ improperBetaDistrE a b
  parseJSON _ = empty

instance Binary BetaDistribution where
  put (BD a b) = put a >> put b
  get = do
    a <- get
    b <- get
    maybe (fail $ errMsgI a b) return $ improperBetaDistrE a b

-- | Create beta distribution. Both shape parameters must be positive.
betaDistr :: Double             -- ^ Shape parameter alpha
          -> Double             -- ^ Shape parameter beta
          -> BetaDistribution
betaDistr a b = maybe (error $ errMsg a b) id $ betaDistrE a b

-- | Create beta distribution. Both shape parameters must be positive.
betaDistrE :: Double             -- ^ Shape parameter alpha
          -> Double             -- ^ Shape parameter beta
          -> Maybe BetaDistribution
betaDistrE a b
  | a > 0 && b > 0 = Just (BD a b)
  | otherwise      = Nothing

errMsg :: Double -> Double -> String
errMsg a b = "Statistics.Distribution.Beta.betaDistr: "
          ++ "shape parameters must be positive. Got a = "
          ++ show a
          ++ " b = "
          ++ show b

-- | Create beta distribution. Both shape parameters must be
-- non-negative. So it allows to construct improper beta distribution
-- which could be used as improper prior.
improperBetaDistr :: Double             -- ^ Shape parameter alpha
                  -> Double             -- ^ Shape parameter beta
                  -> BetaDistribution
improperBetaDistr a b
  = maybe (error $ errMsgI a b) id $ improperBetaDistrE a b

-- | Create beta distribution. Both shape parameters must be
-- non-negative. So it allows to construct improper beta distribution
-- which could be used as improper prior.
improperBetaDistrE :: Double             -- ^ Shape parameter alpha
                   -> Double             -- ^ Shape parameter beta
                   -> Maybe BetaDistribution
improperBetaDistrE a b
  | a >= 0 && b >= 0 = Just (BD a b)
  | otherwise        = Nothing

errMsgI :: Double -> Double -> String
errMsgI a b
  =  "Statistics.Distribution.Beta.betaDistr: "
  ++ "shape parameters must be non-negative. Got a = " ++ show a
  ++ " b = " ++ show b

instance D.Distribution BetaDistribution where
  cumulative (BD a b) x
    | x <= 0    = 0
    | x >= 1    = 1
    | otherwise = incompleteBeta a b x
  complCumulative (BD a b) x
    | x <= 0    = 1
    | x >= 1    = 0
    -- For small x we use direct computation to avoid precision loss
    -- when computing (1-x)
    | x <  0.5  = 1 - incompleteBeta a b x
    -- Otherwise we use property of incomplete beta:
    --  > I(x,a,b) = 1 - I(1-x,b,a)
    | otherwise = incompleteBeta b a (1-x)

instance D.Mean BetaDistribution where
  mean (BD a b) = a / (a + b)

instance D.MaybeMean BetaDistribution where
  maybeMean = Just . D.mean

instance D.Variance BetaDistribution where
  variance (BD a b) = a*b / (apb*apb*(apb+1))
    where apb = a + b

instance D.MaybeVariance BetaDistribution where
  maybeVariance = Just . D.variance

instance D.Entropy BetaDistribution where
  entropy (BD a b) =
    logBeta a b
    - (a-1) * digamma a
    - (b-1) * digamma b
    + (a+b-2) * digamma (a+b)

instance D.MaybeEntropy BetaDistribution where
  maybeEntropy = Just . D.entropy

instance D.ContDistr BetaDistribution where
  density (BD a b) x
    | a <= 0 || b <= 0 = m_NaN
    | x <= 0 = 0
    | x >= 1 = 0
    | otherwise = exp $ (a-1)*log x + (b-1) * log1p (-x) - logBeta a b
  logDensity (BD a b) x
    | a <= 0 || b <= 0 = m_NaN
    | x <= 0 = m_neg_inf
    | x >= 1 = m_neg_inf
    | otherwise = (a-1)*log x + (b-1)*log1p (-x) - logBeta a b

  quantile (BD a b) p
    | p == 0         = 0
    | p == 1         = 1
    | p > 0 && p < 1 = invIncompleteBeta a b p
    | otherwise      =
        error $ "Statistics.Distribution.Gamma.quantile: p must be in [0,1] range. Got: "++show p

instance D.ContGen BetaDistribution where
  genContVar = D.genContinuous


text:
{-# LANGUAGE DeriveDataTypeable, DeriveGeneric #-}
-- |
-- Module    : Statistics.Distribution.Gamma
-- Copyright : (c) 2009, 2011 Bryan O'Sullivan
-- License   : BSD3
--
-- Maintainer  : bos@serpentine.com
-- Stability   : experimental
-- Portability : portable
--
-- The gamma distribution.  This is a continuous probability
-- distribution with two parameters, /k/ and ϑ. If /k/ is
-- integral, the distribution represents the sum of /k/ independent
-- exponentially distributed random variables, each of which has a
-- mean of ϑ.

module Statistics.Distribution.Gamma
    (
      GammaDistribution
    -- * Constructors
    , gammaDistr
    , improperGammaDistr
    -- * Accessors
    , gdShape
    , gdScale
    ) where

import Data.Aeson (FromJSON, ToJSON)
import Control.Applicative ((<$>), (<*>))
import Data.Binary (Binary)
import Data.Binary (put, get)
import Data.Data (Data, Typeable)
import GHC.Generics (Generic)
import Numeric.MathFunctions.Constants (m_pos_inf, m_NaN, m_neg_inf)
import Numeric.SpecFunctions (incompleteGamma, invIncompleteGamma, logGamma, digamma)
import Statistics.Distribution.Poisson.Internal as Poisson
import qualified Statistics.Distribution as D
import qualified System.Random.MWC.Distributions as MWC

-- | The gamma distribution.
data GammaDistribution = GD {
      gdShape :: {-# UNPACK #-} !Double -- ^ Shape parameter, /k/.
    , gdScale :: {-# UNPACK #-} !Double -- ^ Scale parameter, ϑ.
    } deriving (Eq, Read, Show, Typeable, Data, Generic)

instance FromJSON GammaDistribution
instance ToJSON GammaDistribution

instance Binary GammaDistribution where
    put (GD x y) = put x >> put y
    get = GD <$> get <*> get

-- | Create gamma distribution. Both shape and scale parameters must
-- be positive.
gammaDistr :: Double            -- ^ Shape parameter. /k/
           -> Double            -- ^ Scale parameter, ϑ.
           -> GammaDistribution
gammaDistr k theta
  | k     <= 0 = error $ msg ++ "shape must be positive. Got " ++ show k
  | theta <= 0 = error $ msg ++ "scale must be positive. Got " ++ show theta
  | otherwise  = improperGammaDistr k theta
    where msg = "Statistics.Distribution.Gamma.gammaDistr: "

-- | Create gamma distribution. This constructor do not check whether
--   parameters are valid
improperGammaDistr :: Double            -- ^ Shape parameter. /k/
                   -> Double            -- ^ Scale parameter, ϑ.
                   -> GammaDistribution
improperGammaDistr = GD

instance D.Distribution GammaDistribution where
    cumulative = cumulative

instance D.ContDistr GammaDistribution where
    density    = density
    logDensity (GD k theta) x
      | x <= 0    = m_neg_inf
      | otherwise = log x * (k - 1) - (x / theta) - logGamma k - log theta * k
    quantile   = quantile

instance D.Variance GammaDistribution where
    variance (GD a l) = a * l * l

instance D.Mean GammaDistribution where
    mean (GD a l) = a * l

instance D.MaybeMean GammaDistribution where
    maybeMean = Just . D.mean

instance D.MaybeVariance GammaDistribution where
    maybeStdDev   = Just . D.stdDev
    maybeVariance = Just . D.variance

instance D.MaybeEntropy GammaDistribution where
  maybeEntropy (GD a l)
    | a > 0 && l > 0 =
      Just $
      a
      + log l
      + logGamma a
      + (1-a) * digamma a
    | otherwise = Nothing

instance D.ContGen GammaDistribution where
    genContVar (GD a l) = MWC.gamma a l

density :: GammaDistribution -> Double -> Double
density (GD a l) x
  | a < 0 || l <= 0   = m_NaN
  | x <= 0            = 0
  | a == 0            = if x == 0 then m_pos_inf else 0
  | x == 0            = if a < 1 then m_pos_inf else if a > 1 then 0 else 1/l
  | a < 1             = Poisson.probability (x/l) a * a / x
  | otherwise         = Poisson.probability (x/l) (a-1) / l

cumulative :: GammaDistribution -> Double -> Double
cumulative (GD k l) x
  | x <= 0    = 0
  | otherwise = incompleteGamma k (x/l)

quantile :: GammaDistribution -> Double -> Double
quantile (GD k l) p
  | p == 0         = 0
  | p == 1         = 1/0
  | p > 0 && p < 1 = l * invIncompleteGamma k p
  | otherwise      =
    error $ "Statistics.Distribution.Gamma.quantile: p must be in [0,1] range. Got: "++show p

max_stars_repo_name: StefanHubner/statistics

max_issues_repo_name: StefanHubner/statistics

max_forks_repo_name: StefanHubner/statistics


text:
module Main where

import           Control.Concurrent
import           Control.Concurrent.Chan
import           Control.Monad.Writer
import           Data.Array
import           Data.Bits
import qualified Data.ByteString           as B
import           Data.Complex
import           Data.Function
import           Data.List
import           Data.Word
import           Network.Socket            hiding (sendTo)
import           Network.Socket.ByteString
import           System.IO

artWorkLeft :: Rythm
artWorkLeft = Rest 1 :+: Hit 1 1 :+: incrAmp  (incrDur artWorkLeft)

artWorkRight :: Rythm
artWorkRight = Hit 1 1 :+: Rest 1 :+: incrAmp (incrDur artWorkRight)

artWorkMiddle :: Rythm
artWorkMiddle = Hit 1 1 :+: Hit 1 1 :+: go 1 True
  where
    go :: Int -> Bool -> Rythm
    go dur up =
        Rest 1 :+: Hit 119 dur :+: go dur' up'
      where
        dur' =
            if up
                then dur + 1
                else dur - 1
        up' = up && dur < 119 || (not up && dur < 2)

decrAmpN 0 r = r
decrAmpN n r = r' :+: decrAmpN (n - 1) r'
  where
    r' = decrAmp r

incrAmpN 0 r = r
incrAmpN n r = r' :+: incrAmpN (n-1) r'
  where
    r' = incrAmp r

nTimes r 1 = r
nTimes r n = r :+: nTimes r (n-1)

ssht :: Rythm
ssht = Hit 12 1

incrDur (Hit a d) = Hit a (cap 1 120 (d + 1))
incrDur (r1 :+: r2) = incrDur r1 :+: incrDur r2
incrDur ( = r

decrAmp :: Rythm -> Rythm
decrAmp (Hit a d) | a > 1 = Hit (a - 1) d
                  | otherwise = Hit a d
decrAmp (r1 :+: r2) = decrAmp r1 :+: decrAmp r2
decrAmp r = r

incrAmp :: Rythm -> Rythm
incrAmp (Hit a d) | a < 120 = Hit (a + 1) d
                  | otherwise = Hit a d
incrAmp (r1 :+: r2) = incrAmp r1 :+: incrAmp r2
incrAmp r = r

data Action
    = Render Track
    | Stop
    deriving (Show)

ipLeft = "2001:67c:20a1:1095:ba27:ebff:feb9:db12"
ipMiddle = "2001:67c:20a1:1095:ba27:ebff:fe23:60d7"
ipRight = "2001:67c:20a1:1095:ba27:ebff:fe71:dd32"

-- small panel: 16x20
-- panel: 3 x 6 small: 48 x 120

data Panel = P Track (Chan Action)

main :: IO ()
main = do
    l <- startPanel 'L' artWorkLeft (mandelBrot 0 newImage) ipLeft
    m <- startPanel 'M' artWorkMiddle (mandelBrot 49 newImage) ipMiddle
    r <- startPanel 'R' artWorkRight (mandelBrot 97 newImage) ipRight
    hSetBuffering stdin NoBuffering
    hSetBuffering stdout NoBuffering
    xxxLoop l m r
    hSetBuffering stdin LineBuffering
    hSetBuffering stdout LineBuffering
    stopPanel l
    stopPanel m
    stopPanel r

xxxLoop :: Panel -> Panel -> Panel -> IO ()
xxxLoop l m r = do
    i <- readInput
    processInput i l m r

processInput :: Input -> Panel -> Panel -> Panel -> IO ()
processInput Quit _ _ _ = return ()
processInput Tick l m r = do
    let l' = updatePanel l
        m' = updatePanel m
        r' = updatePanel r
    renderToPanel l
    renderToPanel m
    renderToPanel r
    xxxLoop l' m' r'

updatePanel :: Panel -> Panel
updatePanel (P t c) = P (dropCurrent t) c
  where
    dropCurrent (_ :>>: t') = t'
    dropCurrent _ = Done

data Input
    = Quit
    | Tick

readInput :: IO Input
readInput = do
    c <- getChar
    return
        (case c of
             ' ' -> Tick
             _ -> Quit)

startPanel :: Char -> Rythm -> Image -> String -> IO Panel
startPanel name artWork initialImage host = do
    c <- newChan
    ai <- head <$> getAddrInfo Nothing (Just host) (Just "2323")
    s <- socket (addrFamily ai) Datagram defaultProtocol
    void $ forkIO (panelLoop name (addrAddress ai) c s)
    return (P (renderRythm artWork initialImage) c)

stopPanel :: Panel -> IO ()
stopPanel = sendAction Stop

renderToPanel :: Panel -> IO ()
renderToPanel p@(P t _) = sendAction (Render t) p

sendAction :: Action -> Panel -> IO ()
sendAction a (P _ c) = writeChan c a

panelLoop :: Char -> SockAddr -> Chan Action -> Socket -> IO ()
panelLoop name a c s = do
    putStr [name, '?']
    action <- readChan c
    case action of
        Stop -> do
            putStr [name, '$']
            return ()
        Render t -> do
            panelRenderNext t
            panelLoop name a c s
  where
    panelRenderNext t =
        case t of
            Done -> do
                putStr [name, 'X']
                return ()
            Sleep d -> do
                putStr [name, '<']
                threadDelay (durationToMicros d)
                putStr [name, '>']
            Draw i -> do
                putStr [name, '#']
                void $ sendTo s (convertToDisplay i) a
            t' :>>: _ -> panelRenderNext t'

-- * MUSIC API

type Amplitude = Int
type Duration = Int

durationToMicros :: Int -> Int
durationToMicros = (microsPerCol *)
  where
    microsPerCol = 1500000 `div` panelW

data Rythm = Hit Amplitude Duration | Rest Duration | Rythm :+: Rythm
    deriving (Show)

type Micros = Int
data Track = Done | Sleep Micros | Draw Image | Track :>>: Track
    deriving (Show)

instance Monoid Track where
    mempty = Done
    mappend Done t = t
    mappend t Done = t
    mappend (ll :>>: lr) r = ll :>>: mappend lr r
    mappend l r = l :>>: r

renderRythm :: Rythm -> Image -> Track
renderRythm r orig = execWriter (go r orig)
  where
    go :: Rythm -> Image -> Writer Track Image
    go (r1 :+: r2) i = do
        i' <- go r1 i
        go r2 i'
    go (Rest d) i = do
        tell (Sleep (cap 1 120 d))
        return i
    go (Hit a d) i = do
        let i' = flipBlock d a i
        tell (Draw i')
        tell (Sleep (cap 1 120 d))
        return i'

cap l u t | l <= t && t <= u = t
          | l <= t         = u
          | otherwise     = l

-- * Image API
type Image = Array (Int,Int) Bool

panelW, panelH :: Int
panelW = 48
panelH = 120

capCoords (x, y) = (cap 1 48 x, cap 1 120 y)

mandelBrot :: Int -> Image -> Image
mandelBrot offset img = foldr mandelPixel img (indices img)
  where
    mandelPixel (x,y) = setPixel x y (inMandelSet (0.0 :+ 0.0) 0)
      where
        maxIter = 100
        c :: Complex Double
        c =
            (fromIntegral (offset + x) * (4.0 / 144.0) - 2.0) :+
            (fromIntegral y * (4.0 / 120.0) - 2.0)
        inMandelSet z i
          | ((magnitude z) <= 2.0) && i < maxIter =
              inMandelSet (z * z + c) (i + 1)
          | otherwise = rem i 2 == 1

newImage :: Image
newImage =
    array
        ((1, 1), (panelW, panelH))
        [((x, y), False) | x <- [1 .. panelW]
                         , y <- [1 .. panelH]]

flipBlock width height img = foldr (flipCol height) img [1 .. width]

flipCol height col img = foldr flipPixel img [(col, y) | y <- [1..height]]

flipPixel (x,y) img = setPixel x y (not (getPixel x y img)) img

setPixel :: Int -> Int -> Bool -> Image -> Image
setPixel x y v img = img // [((x, y), v)]

getPixel :: Int -> Int -> Image -> Bool
getPixel x y img = img ! (x,y)

convertToDisplay :: Image -> B.ByteString
convertToDisplay img = pixels
  where
    pixels = B.pack $ toFlipDotFormat $ toBit <$> elems img
    toFlipDotFormat :: [Word8] -> [Word8]
    toFlipDotFormat x = go x []
      where
        go x acc =
            case take 8 x of
              [] -> reverse acc
              sx -> go (drop 8 x) (toByte (reverse sx) 1 0 : acc)
    toByte [] _ acc =  acc
    toByte (b:rest) p acc = toByte rest (p*2) (acc .|. b * p)
    toBit :: Bool -> Word8
    toBit = fromIntegral . fromEnum

    (_,(w,h)) = bounds img

toPixel :: Int -> Char
toPixel = toEnum

max_stars_repo_name: sheyll/rythm-flipdot

max_issues_repo_name: sheyll/rythm-flipdot

max_forks_repo_name: sheyll/rythm-flipdot


text:
-- ch1 ex 1.35

import qualified Data.Vector as V

import Statistics.Sample.Histogram (histogram)
import Statistics.Sample (mean, varianceUnbiased, stdDev, range)

samples :: V.Vector Double
samples = V.fromList
  [ 172, 140, 123, 130, 115
  , 148, 108, 129, 137, 161
  , 123, 152, 133, 128, 142
  ]

hist :: (V.Vector Double, V.Vector Int)
hist = histogram 8 samples

main = do
  putStrLn $ "the range/4:            " ++ show (range samples / 4)

  let av = mean samples
      sd = stdDev samples

  putStrLn $ "the standard deviation: " ++ show sd
  putStrLn $ "the mean value:         " ++ show av

  let (bins, hs) = hist
  print bins
  print hs

  let f k = V.length . V.filter (\x -> (av-k*sd)<x && x< (av+k*sd)) $ samples
  putStrLn $ "the sample size:    " ++ show (V.length samples)
  putStrLn $ "within one sigma:   " ++ show (f 1)
  putStrLn $ "within two sigma:   " ++ show (f 2)
  putStrLn $ "within three sigma: " ++ show (f 3)

max_stars_repo_name: raydsameshima/Stat

max_issues_repo_name: raydsameshima/Stat

max_forks_repo_name: raydsameshima/Stat


text:
module Statistics.Quantile.Bench where

import Control.Applicative

import Criterion

import Pipes
import qualified Pipes.Prelude as P

import Statistics.Quantile.Types
import Statistics.Quantile.Util

import System.IO

dummy :: Selector IO
dummy = Selector dummy'

dummy' :: Quantile -> Stream IO -> IO Double
dummy' _ (Stream src) = do
  runEffect $ src >-> P.map show >-> P.stdoutLn
  pure 0

benchSelector :: [Int]
              -> String
              -> Quantile
              -> Selector IO
              -> FilePath
              -> Benchmark
benchSelector ns name q (Selector select) fp =
  bgroup name $ benchSelector' <$> ns
  where benchSelector' n = bench' n fp (select q)

benchBufferingSelector :: (Int -> Int)
                       -> [Int]
                       -> String
                       -> Quantile
                       -> (Int -> Int -> Selector IO)
                       -> FilePath
                       -> Benchmark
benchBufferingSelector f ns name q select fp =
  bgroup name $ benchSelector' <$> ns
  where benchSelector' n = bench' n fp (select' n q)

        select' n = unSelector  (select (f n) n)

bench' :: Int -> FilePath -> (Stream IO -> IO Double) -> Benchmark
bench' n fp f = do
  bench (show n) $ nfIO $ do
    fh <- liftIO $ openFile fp ReadMode
    (f . takeStream n $ streamHandle fh)

bufSize :: (Int -> Int)
        -> Int
        -> Int
bufSize f n
  | n < 10000 = n
  | otherwise = f n

sqrtN :: Int -> Int
sqrtN = bufSize (floor . sqrt . fromIntegral)

logN :: Int -> Int
logN = bufSize (floor . logBase 2 . fromIntegral)

max_stars_repo_name: fractalcat/slides

max_issues_repo_name: fractalcat/slides

max_forks_repo_name: fractalcat/slides