text:
{-|
Module      : ODMatrix.ElementaryDecomposition.ShortestPath
Description : This module is responsible for computing the shortest path between two ODMs.

Contains the necesary functions to measure and compare ODMs.

-}
module ODMatrix.ElementaryDecomposition.ShortestPath (
      Measure
    , sePath
    , ev
    , ev3
  ) where

  import Control.Monad.State.Lazy
  --import Control.Monad.Loops (iterateUntilM)

  import Numeric.LinearAlgebra

  import qualified Data.HashMap.Lazy as M
  import qualified Data.Hashable as H
  import Data.Maybe

  import ODMatrix.ElementaryDecomposition (ElementaryMatrix(..), ODM,
                                       getChildren, applyElementary,
                                       applyPath,
                                       applicables, elementaryValue,
                                       opposite, isNullElementary)

  -- import Debug.Trace (trace)

  -- | Measure or Metric prototype
  type Measure = Int -> Int -> Double

  -- | Elementary Value.
  -- Given a metric m, a in route metric mr and a scale factor alpha returns a elementary value metric.
  ev :: Double    -- ^ Scale factor alpha
     -> Measure   -- ^ Metric m
     -> Measure   -- ^ In route metric mr
     -> Measure   -- ^ Elementary value metric
  ev a m mr i j
    | a < 0 = error "The ajustment parameter must be greater than 0."
    | otherwise = m i j + a * mr i j

  -- | Shortest Path
  -- Computes the sequence of elementary matrices between ODMs
  sePath :: Measure             -- ^ Elementary value metric
         -> ODM                 -- ^ Source ODM
         -> ODM                 -- ^ Target ODM
         -> [ElementaryMatrix]  -- ^ Sequence of elementary matrices in the shortest path
  sePath m a b = evalState _process (_initialSEPState m a b)

  -- | Map a Elementary Matrix to it related cells indexes
  emc :: ElementaryMatrix -> (Int, Int)
  emc (Elementary _ (cell1, _) (cell2, _)) = (cell1,cell2)

  -- | Compares two ODMs
  odmc :: (Matrix Double -> Double) -- ^ Function g to compute the remainder
       -> Measure                   -- ^ Elementary value metric
       -> Double                    -- ^ Scale factor
       -> ODM                       -- ^ ODM A
       -> ODM                       -- ^ ODM B
       -> Double                    -- ^ Difference between A and B
  odmc g ev beta a b = beta * g s + sum [ev i j | (i,j) <- sp]
    where p = sePath ev a b
          sp = map emc . filter (not . isNullElementary) $ p
          s = b - applyPath a p -- Remainder

  newtype HashMatrix = HashMatrix { hm :: ODM }

  instance Show HashMatrix where
    show (HashMatrix m) = show m

  instance Eq HashMatrix where
    (==) (HashMatrix a) (HashMatrix b) = a == b

  instance H.Hashable HashMatrix where
    hashWithSalt salt (HashMatrix m) = floor $ m `atIndex` (salt `mod` rows m, salt `mod` cols m)

  data SEP = SEP {
      dist      :: M.HashMap HashMatrix Double
    , prev      :: M.HashMap HashMatrix ElementaryMatrix
    , visited   :: [ODM]
    , target    :: ODM
    , measure   :: Measure
    }

  instance Show SEP where
    show = show . visited

  -- | Elementary Value 3 (See the paper).
  -- This is a special case of elementary value computed with only one metric.
  -- Only works for regular grids.
  ev3 :: Double -> Measure -> Measure
  ev3 a m i j
    | a <= 0 || a >= 1 = error "The ajustment parameter a must be 0 < a < 1."
    | otherwise = a * m i j + (1-a) * sum [ m k (k+1) | k <- [min i j .. max i j] ]

  -- TODO: Implement metrics.

  _initialSEPState :: Measure -> ODM -> ODM -> SEP
  _initialSEPState m src trg = SEP (M.singleton (HashMatrix src) 0) M.empty [] trg m

  -- | TODO: take out the harcoded 60.
  -- | MAYBE IF IN THE FINAL STEP ASK FOR NEIGHBORS EMPTY.
  _process :: State SEP [ElementaryMatrix]
  _process = do
    s <- get
    let (curr,d) = _getNext (dist s)
        neighbors = applicables (0,60) curr . getChildren $ (target s) - curr
        -- Filter out the visited neighbors
        unvisited = filter (\e -> applyElementary curr e `notElem` visited s) neighbors
    -- Calculate min distance to all neighbors (unvisited)
    modify $ minAndPrev curr d unvisited
    -- \t -> t {dist = foldr (mindist (measure t) curr d) (dist t) unvisited}
    -- Add current to the previous of unvisited

    -- Mark the current as visited.
    modify $ \t -> t {dist = M.delete (HashMatrix curr) (dist t),
                      visited = curr:visited t}
    sn <- get
    if curr == (target sn) || null neighbors  then return . fst $ _wrapUp curr sn else _process
    -- if curr == (target sn) then return . fst $ _wrapUp curr sn else _process

  -- Calculate the min distance and the corresponding previous.
  minAndPrev :: ODM -> Double -> [ElementaryMatrix] -> SEP -> SEP
  minAndPrev curr d unvisited t = foldr f t nds
    where nds = map (neighborsDists t curr d) unvisited
          f (md,k,e) s = s { dist = M.insert k md (dist s),
                             prev = M.insert k e (prev s) }

  neighborsDists :: SEP                           -- ^ Final state
                 -> ODM                           -- ^ Current
                 -> Double                        -- ^ Dist from origin
                 -> ElementaryMatrix              -- ^ Neighbor
                 -> (Double,HashMatrix,ElementaryMatrix)     -- ^ Min distance and Neighbor
  neighborsDists s m d e = (min newd oldd, k, e)
    where newd = d + elementaryValue (measure s) e
          oldd = M.lookupDefault newd k (dist s)
          k = HashMatrix $ applyElementary m e

  -- | Get next matrix to process and the shortest distance from src
  _getNext :: M.HashMap HashMatrix Double -- ^ List of distances by ODM
           -> (ODM,Double)                -- ^ ODM with min distance and the distance
  _getNext ds
    | null ds = error "Not destination found"
    | otherwise = (\(m,d) -> (hm m, d)) $ M.foldrWithKey f initial ds
    where initial = head . M.toList $ ds            -- (k,v)
          f k v mn@(_, minv) | v < minv = (k,v)
                             | otherwise = mn

  _wrapUp :: ODM -> SEP -> ([ElementaryMatrix],SEP)
  _wrapUp m s = (es, s)
    where p = M.lookup (HashMatrix m) (prev s)
          --op e = e {sign = - sign e}
          odm = (\x -> applyElementary m (opposite x)) <$> p
          es = maybe [] (\x -> fromJust p : fst (_wrapUp x s)) odm

  -- _testMatrix :: ODM
  -- _testMatrix = matrix 10 [0,0,0,1,0,0,0,0,0,0,0,2,1,0,1,3,0,0,0,0,0,0,2,0,0,4,1,0,1,0,0,0,0,1,0,0,0,0,2,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]

  -- _testSimMatrix :: ODM
  -- _testSimMatrix = matrix 10 [0,0,1,0,0,0,0,0,0,0,0,2,0,1,1,3,0,0,0,0,0,0,2,0,0,4,1,0,1,0,0,0,0,1,0,0,0,0,2,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]

  -- _testDifMatrix :: ODM
  -- _testDifMatrix = matrix 10 [0,0,0,1,0,0,0,0,0,0,0,2,1,0,1,3,0,0,0,0,0,0,2,0,0,4,2,0,0,0,0,0,0,1,0,0,0,0,2,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]

max_stars_repo_name: renecura/odmatrix

max_issues_repo_name: renecura/odmatrix

max_forks_repo_name: renecura/odmatrix


text:
module Parser where

import Data.Ratio
import Data.Complex

import Number
import LispData
import Text.ParserCombinators.Parsec hiding (spaces)
import Text.Parsec.Number

import LispData
--TODO Learn monads again and check what mapM is doing

-- Recognizes if a character is a valid scheme symbol
symbol :: Parser Char
symbol = oneOf "-.!#$%&|*+/:<=>?@^_~"
-- Skips one or more spaces
spaces1 :: Parser ()
spaces1 = skipMany1 space

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

-- Parses a string which starts with a " and ends with a"
-- TODO \\t \\n \\r \\ \"
parseString :: Parser LispVal
parseString = do
  _ <- char '"' --Starts with a "
  -- Stops at "
  x <- many $ escapedChars <|> noneOf "\""
  _ <- char '"' --ends with a "
  return $ String x

-- Parses a symbol
parseAtom :: Parser LispVal
parseAtom = do
  first <- (letter <|> symbol) <?> "I HATE PARSEC"
  --first <- choice symlpars
  --the following chars must be one of letter, digit or symbol
  rest <-  many (letter <|> digit <|> symbol)
  let atom = first :rest
  --catch special atoms
  case atom of "#t" -> return $ Bool True
               "#f" -> return $ Bool False
               ('-':x:_)   -> case x of
                   ' ' -> return $ Atom atom
                   _   -> parseNumber
                             --return $ (LispNumber . Integer . read) atom --TODO THIS IS NOT GOOD, try to parse -3o
               _    -> {-trace ("attom"++ show atom)-} (return $ Atom atom)

parseComplex :: Parser LispVal
parseComplex = do
  f <- parseNumber
  s <- char '+'
  n <- parseNumber
  i <- char 'i'
  let x = case f of
        ((LispNumber (Integer n))) -> (fromIntegral n :: Double)
        ((LispNumber (Rational n))) -> (realToFrac n)
        ((LispNumber (Real n))) -> n
        _ -> 0

  let y = case n of
        ((LispNumber (Integer p))) -> (fromIntegral p :: Double)
        ((LispNumber (Rational p))) -> (realToFrac p)
        ((LispNumber (Real p))) -> p
        _ -> 0
  return $ (LispNumber . Complex) (x :+ y)

parseNegFloat :: Parser LispVal
parseNegFloat = do
  s <-  sign
  beforeDot <- int
  _ <- char '.' <?> "Floating Point Parse Error: expecting ."
  afterDot <- int
  let d = s (read (show beforeDot ++ "." ++ (show afterDot)))
  return $ (LispNumber . Real) d

parseNegRational :: Parser LispVal
parseNegRational = do
  top <- int
  _ <- char '/'
  bottom <- int
  return $ (LispNumber . Rational) (top % bottom)

parseInteger :: Parser LispVal
parseInteger =  do
  int' <- many1 digit
  return $ (LispNumber . Integer . read) int'

parseNegInteger :: Parser LispVal
parseNegInteger =  do
  _ <-  char '-'
  int' <- many1 digit
  return $ (LispNumber . Integer . negate . read) int'

parseVector :: Parser LispVal
parseVector = do
  _ <- char '('
  vec <- sepBy parseExpr spaces1
  _ <- char ')'
  return $ Vector vec

parseList :: Parser LispVal
-- Parse lispExpr which hare seperated by one or more whitespace
parseList = List <$> sepBy parseExpr spaces1

parseDottedList :: Parser LispVal
parseDottedList = do
    -- what?
    head' <- endBy parseExpr spaces1
    -- parses a dot then exactly one space and saves the expression after the space
    tail' <- char '.' >> spaces1 >> parseExpr
    return $ DottedList head' tail'

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

parseNumber = try parseNegFloat
              <|> try parseNegRational
              <|> try parseNegInteger --TODO WENN MIR WAS UM DIE OHREN FLIEGT LIEGTS HIER DRAN
              <|> parseInteger -- etc

parseExpr :: Parser LispVal
parseExpr = try parseAtom --first try to parse a atom
         <|> try parseComplex
         <|> parseNumber
         <|> parseString -- if this fails try to parse a string
         <|> parseQuoted
         <|> do _ <- char '('
                -- parses a normal list until it encounter a dot, at which point it will go back and sstart to parse
                -- a dotted list
                x <- try parseList <|> parseDottedList
                _ <- char ')'
                return x

ext: hs

max_stars_repo_name: Unaimend/schemehs

max_issues_repo_name: Unaimend/schemehs

max_forks_repo_name: Unaimend/schemehs


text:
module School.Unit.RecLin
( recLin ) where

import Numeric.LinearAlgebra (R, cmap, cond)
import School.Unit.Unit (Unit(..))
import School.Unit.UnitActivation (UnitActivation(..))
import School.Unit.UnitGradient (UnitGradient(..))
import School.Unit.UnitParams (UnitParams(..))

recLin :: Unit R
recLin = Unit { apply = recLinApply
              , deriv = recLinDeriv
              }

recLinApply :: UnitParams R
            -> UnitActivation R
            -> UnitActivation R
recLinApply EmptyParams (BatchActivation inMatrix) =
  BatchActivation $ cmap (max 0) inMatrix
recLinApply _ f@(ApplyFail _) = f
recLinApply _ _ = ApplyFail "Failure when applying RecLin unit"

recLinDeriv :: UnitParams R
            -> UnitGradient R
            -> UnitActivation R
            -> ( UnitGradient R
               , UnitParams R)
recLinDeriv EmptyParams
            (BatchGradient inGrad)
            (BatchActivation input)
  = (BatchGradient grad, EmptyParams)
    where grad = cond input 0 0 0 inGrad
recLinDeriv EmptyParams f@(GradientFail _) _ = (f, EmptyParams)
recLinDeriv EmptyParams _ (ApplyFail msg) = (GradientFail msg, EmptyParams)
recLinDeriv _ _ _ = (GradientFail "Failure when differentiating RecLin unit" , EmptyParams)

max_stars_repo_name: jfulseca/School

max_issues_repo_name: jfulseca/School

max_forks_repo_name: jfulseca/School


text:
--
-- Copyright (c) 2009-2011, ERICSSON AB
-- All rights reserved.
--
-- Redistribution and use in source and binary forms, with or without
-- modification, are permitted provided that the following conditions are met:
--
--     * Redistributions of source code must retain the above copyright notice,
--       this list of conditions and the following disclaimer.
--     * Redistributions in binary form must reproduce the above copyright
--       notice, this list of conditions and the following disclaimer in the
--       documentation and/or other materials provided with the distribution.
--     * Neither the name of the ERICSSON AB nor the names of its contributors
--       may be used to endorse or promote products derived from this software
--       without specific prior written permission.
--
-- THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
-- AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
-- IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
-- DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
-- FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
-- DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
-- SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
-- CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
-- OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
-- OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
--

module Feldspar.Core.Frontend.Eq
where

import qualified Prelude as P

import Data.Int
import Data.Word
import Data.Complex

import Feldspar.Prelude
import Feldspar.Core.Types
import Feldspar.Core.Constructs
import Feldspar.Core.Constructs.Eq

infix 4 ==
infix 4 /=

-- | Redefinition of the standard 'P.Eq' class for Feldspar
class (Type a) => Eq a
  where
    (==) :: Data a -> Data a -> Data Bool
    (==) = sugarSymF Equal
    (/=) :: Data a -> Data a -> Data Bool
    (/=) = sugarSymF NotEqual

instance Eq ()
instance Eq Bool
instance Eq Float
instance Eq Double
instance Eq Word8
instance Eq Word16
instance Eq Word32
instance Eq Word64
instance Eq WordN
instance Eq Int8
instance Eq Int16
instance Eq Int32
instance Eq Int64
instance Eq IntN

instance (Eq a, Eq b)                               => Eq (a,b)
instance (Eq a, Eq b, Eq c)                         => Eq (a,b,c)
instance (Eq a, Eq b, Eq c, Eq d)                   => Eq (a,b,c,d)
instance (Eq a, Eq b, Eq c, Eq d, Eq e)             => Eq (a,b,c,d,e)
instance (Eq a, Eq b, Eq c, Eq d, Eq e, Eq f)       => Eq (a,b,c,d,e,f)
instance (Eq a, Eq b, Eq c, Eq d, Eq e, Eq f, Eq g) => Eq (a,b,c,d,e,f,g)

instance (Eq a, RealFloat a) => Eq (Complex a)

max_stars_repo_name: kffaxen/feldspar-language

max_issues_repo_name: kffaxen/feldspar-language

max_forks_repo_name: kffaxen/feldspar-language


text:
module Strategy.RandomWalkMetropolis where

import Math.Probably.Sampler
import Math.Probably.MCMC
import Numeric.LinearAlgebra
import Control.Monad
import Data.Maybe

rwm = VStrategy rwmTrans rwmIni

rwmTrans posterior xi (sigma,(i, iaccept)) mpi = do
   let pi = fromMaybe (posterior xi) mpi
   xstar <- fmap fromList $ mapM (\x -> normal x sigma) $ toList xi
   let pstar =  posterior xstar
   let ratio = exp $ pstar - pi
   u <- unit
   let accept = u < ratio

   --diminishing adaptation:
   let sigmaNext = if accept then sigma*(min 1.4 $ 1+5.0/i)^3
                             else sigma*(max 0.7143 $ 1-5.0/i)

   return $ if accept
               then ((xstar, (sigmaNext, (i+1, iaccept+1))), Just pstar)
               else ((xi, (sigmaNext, (i+1, iaccept))), Just pi)

rwmIni _ = (0.1,(1.0, 0.0))

max_stars_repo_name: glutamate/probably-baysig

max_issues_repo_name: silky/probably-baysig

max_forks_repo_name: silky/probably-baysig


text:
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE BangPatterns #-}

{-# OPTIONS_GHC -fno-warn-orphans #-} -- do not warn on local instance of Random for Complex

module BBrotCompute(compute, selectCells, inSet) where

import Control.Monad.Par
import Data.Aeson(encode)
import Data.Array
import qualified Data.ByteString.Lazy as BS
import Data.Complex
import Data.List
import Data.Maybe
import Data.Time
import System.Environment
import System.Console.CmdArgs(whenNormal, whenLoud)
import System.Random
import Text.Printf

import BBrotConf
import BBrotSelection

instance Random (Complex Double) where
  randomR (!(x0 :+ y0), !(x1 :+ y1)) g = (r :+ i, g2)
    where (!r, !g1) = randomR (x0, x1) g
          (!i, !g2) = randomR (y0, y1) g1
  random = randomR (0.0 :+ 0.0, 1.0 :+ 1.0)

data GridPoint = GridPoint { xbase :: !Double
                           , ybase :: !Double
                           , inMandelbrotSet :: !Bool
                           }
                 deriving Show

makeGrid :: Double -> Int -> Array (Int, Int) GridPoint
makeGrid step bailout = array ((0, 0), (xsteps, ysteps)) points
  where (xmin :+ ymin) = loCorner
        (xmax :+ ymax) = hiCorner
        xsteps = floor $ (xmax - xmin) / step :: Int
        ysteps = floor $ (ymax - ymin) / step :: Int
        points = [ ((i,j), g i j) | i <- [0..xsteps], j <- [0..ysteps] ]
        g i j = GridPoint x y p
          where x = xmin + fromIntegral i * step
                y = ymin + fromIntegral j * step
                p = not $ inSet 0 bailout (x :+ y)

selectCells :: Double -> Int -> [(Double, Double)]
selectCells step bailout = [ (xbase g, ybase g) | g <- cells ]
  where !grid = makeGrid step bailout
        ((x0, y0), (xsteps, ysteps)) = bounds grid
        !cells = [ grid!(i,j) | i <- [x0 + 1 .. xsteps - 1],
                                j <- [y0 + 1 .. ysteps - 1],
                                isNearBorder i j ]
        isNearBorder i j = someIn && someOut
            where neighbours = [ grid!(i+a,j+b) | (a, b) <- range ((-1,-1), (1,1)) ]
                  inMandel = map inMandelbrotSet neighbours
                  someIn = or inMandel
                  someOut = not $ and inMandel

iterations :: Int -> Double -> Double -> Double -> Double -> Int -> Int
iterations !maxK !x !y !x0 !y0 !k =
  let x2 = x * x
      y2 = y * y
  in
   if k == maxK || x2 + y2 > 4
   then k
   else iterations maxK (x2 - y2 + x0) (2 * x * y + y0) x0 y0 (k + 1)

inCardioBulb :: Double -> Double -> Bool
inCardioBulb !x !y = inCardio || inBulb
  where sqr a = a * a
        q = sqr (x - 1/4) + sqr y
        inCardio = q * (q + (x - 1/4)) < sqr y / 4
        inBulb = sqr (x + 1) + sqr y < 1 / 16

toMaybeBBPoint :: Int -> Int -> Complex Double -> Maybe BBPoint
toMaybeBBPoint !minK !maxK !z =
  let x = realPart z
      y = imagPart z
      k = iterations maxK x y x y 0
  in if not (inCardioBulb x y) && k >= minK && k < maxK
     then Just (BBPoint x y k)
     else Nothing

inSet :: Int -> Int -> Complex Double -> Bool
inSet !minK !maxK !z = p $ toMaybeBBPoint minK maxK z
  where p (Just _) = False
        p Nothing  = True

toUnit :: Int -> String
toUnit n = show (n `div` d) ++ unit
  where pow :: Int -> Int -> Int
        pow x y = x ^ y
        l = [ (10 `pow` 9, "G"), (10 `pow` 6, "M"), (10 `pow` 3, "K") ]
        (d, unit) = fromMaybe (1, "") $ find ((<= n) . fst) l

compute :: BBrotConf -> IO ()
compute conf = do
  -- Pick random points in the complex plane, test their orbits,
  -- save the interesting ones.

  whenNormal $ putStrLn $ printf "Sampling %s points..." $ toUnit $ samples conf

  initRandGen <- case seed conf of
    Just s -> return $ mkStdGen s
    Nothing -> getStdGen

  let step = gridStep conf
      cells = selectCells step 255
      pointsPerCell = samples conf `div` length cells
      generators = iterate (fst . split) initRandGen
      genPoints (cell, gen) = take pointsPerCell $ randomRs (corners cell) gen
      corners (x, y) = ((x :+ y) - cellDiag / 2, (x :+ y) + cellDiag / 2)
      cellDiag = step :+ step

  let pointLists = map genPoints $ zip cells generators

  let selected = concat $ concat $ runPar $ parMap (map f) pointLists
        where f = maybeToList . toMaybeBBPoint (minIters conf) (maxIters conf)

  let cachefile = fromMaybe defPath (ocachepath conf)
        where defPath = printf "/tmp/buddhabrot-%s-%s_%s.bbc"
                        (toUnit $ samples conf)
                        (toUnit $ minIters conf)
                        (toUnit $ maxIters conf)

  whenNormal $ do
    putStrLn $ "Selected cells: " ++ show (length cells)
    putStrLn $ "Writing cache " ++ cachefile ++ " ..."
  args <- getArgs
  currentTime <- getCurrentTime
  BS.writeFile cachefile $ encode $
    PointSelection selected (Just args) (Just $ show initRandGen) (Just $ show currentTime)
  whenLoud $ putStrLn $ "selected points: " ++ show (length selected)

max_stars_repo_name: saffroy/buddhabrot

max_issues_repo_name: saffroy/buddhabrot

max_forks_repo_name: saffroy/buddhabrot


text:
{-|
Module : Network
|-}
module Network
  ( Net(..)
  , createNet
  , saveNetwork
  , loadNetwork
  ) where

import Numeric.LinearAlgebra
import Specs
import System.Random

-- | data representing Network
data Net = Net {
    numLayers :: Int -- ^ Number of layers in Network
  , sizes :: [Int] -- ^ Sizes of each layer
  , biases :: [Matrix Double] -- ^ Bias Matrix's of each layer
  , weights :: [Matrix Double] -- ^ Weight Matrix's of each layer
  } deriving (Eq, Show)

-- | Generate random small numbers
smallRandoms :: Int -> [Double]
smallRandoms seed = map (/100) (randoms (mkStdGen seed))

-- | Create Matrix with random weigths
randomWeightMatrix :: Int -> Int -> Int -> Matrix Double
randomWeightMatrix numInputs numOutputs seed = (numOutputs><numInputs) weights
    where weights = take (numOutputs*numInputs) (smallRandoms seed)

-- | Create new Network for given sizes of layers
createNet :: [Int] -- ^ Sizes of each layer
          -> IO Net
createNet sizes = do
   let b = (\ x -> randomWeightMatrix 1 x 7) <$> tail sizes
   let w = (\ (x, y) -> randomWeightMatrix x y 7) <$> zip (init sizes) (tail sizes)
   let net = Net{numLayers = length sizes, sizes = sizes, biases = b, weights = w}
   return net

-- | Change Network to String for saving
toString :: Net -> String
toString net = nL ++ s ++ b ++ w
  where
    nL = show (numLayers net) ++ " \n"
    s = foldr ((\ a b -> a (' ' : b)) . shows) "\n" (sizes net)
    b = foldr1 (++) (foldr ((\ a b -> a (' ' : b)) . shows) "\n" <$> (foldr1 (++) <$> (toLists <$> biases net)))
    w = foldr1 (++) (foldr ((\ a b -> a (' ' : b)) . shows) "\n" <$> (foldr1 (++) <$> (toLists <$> weights net)))

-- | Save Network to file
saveNetwork :: Net -- ^ Network to save
            -> FilePath -- ^ FilePath for saving
            -> IO ()
saveNetwork net file = writeFile file (toString net)

-- | Change line from file to list of Doubles
lineToListDouble :: String -> [Double]
lineToListDouble line = (\ x -> read x :: Double) <$> words line

-- | Load Network from file
loadNetwork :: FilePath -- ^ FilePath to saved network
            -> IO Net
loadNetwork file = do
  content <- readFile file
  let contentLines = lines content
  let numLayers = (read (contentLines !! 0) :: Int)
  let sizes = (\x-> read x :: Int) <$> (words (contentLines !! 1))
  let ba = (contentLines !!) <$> [2..numLayers]
  let b = (\ (x, y) -> (y><1) x :: Matrix Double) <$> zip (lineToListDouble <$> ba) (tail sizes)
  let we = (contentLines !!) <$> [(numLayers+1)..(numLayers+numLayers-1)]
  let w  = (\ (x, y, z) -> (z><y) x :: Matrix Double) <$> zip3 (lineToListDouble <$> we) (init sizes) (tail sizes)
  return Net{numLayers = length sizes, sizes = sizes, biases = b, weights = w}

max_stars_repo_name: Malenczuk/HuskNet

max_issues_repo_name: Malenczuk/HuskNet

max_forks_repo_name: Malenczuk/HuskNet


text:
module LogisticRegressionExercise where

import qualified Graphics.Gnuplot.Advanced as GP

import qualified Graphics.Gnuplot.Graph as Graph

import qualified Graphics.Gnuplot.Plot.ThreeDimensional as Plot3D
import qualified Graphics.Gnuplot.Graph.ThreeDimensional as Graph3D

import qualified Graphics.Gnuplot.Plot.TwoDimensional as Plot2D
import qualified Graphics.Gnuplot.Graph.TwoDimensional as Graph2D
import qualified Graphics.Gnuplot.Terminal.Default as DefaultTerm
import Graphics.Gnuplot.Plot.TwoDimensional (linearScale, )

import qualified Data.ByteString.Lazy as B
import Data.Csv.HMatrix
import Data.Csv.Incremental
import HVX
import LogisticRegression
import Numeric.LinearAlgebra
import Prelude hiding ((<>))
import System.Directory
import Util

-- replace this with the path to your data.
exerciseTwoDataOne :: IO String
exerciseTwoDataOne = do
  home <- getHomeDirectory
  return (home ++ "/code/haskell/CourseraMachineLearning/app/ex2data1.txt")

exerciseTwoDataTwo :: IO String
exerciseTwoDataTwo = do
  home <- getHomeDirectory
  return (home ++ "/code/haskell/CourseraMachineLearning/app/ex2data2.txt")

costFunction :: Matrix Double -> Matrix Double -> Matrix Double -> (Double, Matrix Double)
costFunction theta x y = ((sumElements cost) / m, grad)
  where sig = sigmoid $ x <> theta
        m = fromIntegral $ fst (size y)
        n = snd (size x)
	-- Cost function for logistic regression. Faithful formula translation.
        cost = (((-1) * (tr y)) <> (log sig)) - ((1 - (tr y)) <> (log (1 - sig)))
	-- Gradient. Faithful translation of the formula.
        grad = (columnWise sumElements (x * repmat (sig - y) 1 n)) / (scalar m)

plotData :: Matrix Double -> IO (Plot2D.T Double Double)
plotData completeMatrix = do
  let matrixRows = toRows completeMatrix
  -- select all rows whose third column is 1.
  let pos = fromRows $ filter (((==) 1) . (flip atIndex 2)) matrixRows
  -- do the same for those whose third column is 0.
  let neg = fromRows $ filter (((==) 0) . (flip atIndex 2)) matrixRows
  -- Helper function: Convert column vectors to lists for plotting.
  let getPointLists p = map toList $ toColumns (p ?? (All, Take 2))
  let [xPos, yPos] = getPointLists pos
  let [xNeg, yNeg] = getPointLists neg
  return $ plotPoints xPos yPos `mappend` plotPoints xNeg yNeg

plotDecisionBoundary :: Matrix Double -> Matrix Double -> IO (Plot2D.T Double Double)
plotDecisionBoundary theta x = do
  let theta' = flatten theta
  -- We only need two points to plot a line.
  -- So we pick the smallest and largest.
  let lowerPoint = minimum ((toList . flatten) $ x ¿ [1]) - 2
  let upperPoint = maximum ((toList . flatten) $ x ¿ [1]) + 2
  let plotX = [lowerPoint, upperPoint]
  -- A helper function to calculate the y values given by our theta.
  let f val = (-1) / (atIndex theta' 2) * ((atIndex theta' 1) * val + (atIndex theta' 0))
  let plotY = [f lowerPoint, f upperPoint]
  return $ plotLine plotX plotY

exerciseTwo = do
  dataOne <- exerciseTwoDataOne
  contents <- B.readFile dataOne
  let completeMatrix = decodeMatrix NoHeader contents
  let x'' = completeMatrix ?? (All {- Rows -}, Take 2 {- Columns -})
  let y = completeMatrix ¿ [2]
  let (m, n) = size x''
  let x = ones m 1 ||| x''
  let initialTheta = zeros (n + 1) 1
  let (cost, grad) = costFunction initialTheta x y
  putStrLn $ "Cost at initial theta (zeros): " ++ (show cost)
  putStrLn $ "Expected cost (approx): 0.693"
  putStrLn $ "Gradient at initial theta (zeros): " ++ (show grad)
  putStrLn $ "Expected gradients (approx):\n -0.1000\n -12.0092\n -11.2628\n"
  -- Compute and display cost and gradient with non-zero theta.
  -- the learning rate was achieved through empirical tests.
  -- Gives theta closest to fminunc.
  -- TODO(mchav): Find haskell optimisation library.
  let (_, theta) = sgd x y initialTheta costFunction 0.0041 1000000
  decisionBoundary <- plotDecisionBoundary theta x
  points <- plotData completeMatrix
  GP.plotSync DefaultTerm.cons $ decisionBoundary `mappend` points

max_stars_repo_name: mchav/CourseraMachineLearning

max_issues_repo_name: mchav/CourseraMachineLearning

max_forks_repo_name: mchav/CourseraMachineLearning


text:
#!/usr/bin/env stack
-- stack --resolver lts-11.1 --install-ghc runghc --package random
{-# LANGUAGE TypeApplications #-}
import Data.Complex (Complex(..), magnitude)
import Control.Monad (replicateM, forM_)
import System.Random (randomIO)
import Text.Printf (printf)

approximate :: Int -> IO Double
approximate count = do
    let randomPoint = (:+) <$> randomIO <*> randomIO @Double
        randomPoints = flip replicateM randomPoint
    inCircle <- (length . filter (\z -> magnitude z <= 1)) <$> randomPoints count
    pure $ 4 * (fromIntegral inCircle / fromIntegral count)

main = forM_ [1 .. 100000] $ \n -> printf "%6d %f\n" n =<< approximate n

max_stars_repo_name: kmein/a-pi-roximate

max_issues_repo_name: kmein/a-pi-roximate

max_forks_repo_name: kmein/a-pi-roximate


text:
import Data.Either
import DistanceGeometry
import Generators
import ReadWrite
import Numeric.LinearAlgebra
import Numeric.LinearAlgebra.Data
import System.Directory

idir = "test/molecules"
odir = "test/results"

main :: IO ()
main = do
  putStrLn ""
  putStrLn "generateOneMolecule: Begin"
  generateOneMolecule idir odir "example" 16
  putStrLn "generateOneMolecule: End"
  putStrLn ""
  putStrLn "generateManyMolecules: Begin"
  generateManyMolecules idir odir "example" 16 10
  putStrLn "generateManyMolecules: End"

max_stars_repo_name: wurthel/distance-geometry

max_issues_repo_name: wurthel/distance-geometry

max_forks_repo_name: wurthel/distance-geometry


text:
{-# LANGUAGE NamedFieldPuns #-}

module Word2Vec ( trainWith
                , mostSimilar
                , groupOf
                ) where

import Vocabulary
import Huffman
import Train
import Persistence

import Prelude hiding (lookup)
import Numeric.LinearAlgebra.Data (Matrix, Vector, loadMatrix, saveMatrix, cmap, size)
import qualified Numeric.LinearAlgebra.HMatrix as M
import qualified Numeric.LinearAlgebra.Devel as DV
import qualified Data.HashMap.Strict as HM
import qualified Data.IntMap.Strict as IM
import Data.List (sortBy)
import Data.Ord (comparing)
import Data.Foldable (fold)
import Control.Arrow (first)
import System.Random (getStdGen)
import System.IO

i2hFile = "i2h.txt"
i2wFile = "i2w.txt"
w2iFile = "w2i.txt"

-- | Train input-hidden weight matrix with sentences from a file.
trainWith :: FilePath -> Int -> IO ()
trainWith path num = do
  sentences <- take num <$> readGiganticLine path

  randGen <- getStdGen
  let (vocabulary, i2w, vocSize) = buildModel sentences
      huffman = encodeTree $ buildTree $ map snd $ HM.toList vocabulary
      vectorSize = 100
      matrices = buildMatrices randGen vocSize vectorSize
      config = Configuration {vocabulary, huffman, learningRate=0.025, vectorSize}
      (i2h, h2o) = train config matrices sentences

  saveMatrix i2hFile "%g" i2h
  saveI2W i2wFile i2w
  saveW2I w2iFile vocabulary

mostSimilar :: [String] -> IO ()
mostSimilar [] = putStrLn "Specify words."
mostSimilar words = do
  i2h <- normalizeMatrix <$> loadMatrix i2hFile
  i2w <- loadI2W i2wFile
  w2i <- loadW2I w2iFile
  mapM_ print $ take 5 $ findMostSimilar i2w w2i i2h words

findMostSimilar :: IM.IntMap String
                -> HM.HashMap String Int
                -> Matrix Double
                -> [String]
                -> [(String, Double)]
findMostSimilar i2w w2i i2h words = filter ((`notElem` nakedWords) . fst) similarWords
  where weightedWords = map weightedWord words
        nakedWords = map fst weightedWords
        similarities = M.toList $ M.app i2h (wordsMean w2i i2h weightedWords)
        indexSimilarities = sortBy (flip $ comparing snd) $ zip [0..] similarities
        similarWords = map (first (i2w IM.!)) indexSimilarities

-- | Normalized mean vector for multiple words.
wordsMean :: HM.HashMap String Int -> Matrix Double -> [(String, Double)] -> Vector Double
wordsMean w2i i2h words = meanVector vecs
  where vecs = map (\(s, w) -> M.scale w $ wordVector w2i i2h s) words

-- | Weigh a word according to its prefix.
weightedWord :: String -> (String, Double)
weightedWord word = case word of
  '+':rest -> (rest, 1)
  '-':rest -> (rest, -1)
  word     -> (word, 1)

-- | Return a vector for a word.
wordVector :: HM.HashMap String Int -> Matrix Double -> String -> Vector Double
wordVector w2i i2h = (i2h M.!) . (w2i HM.!)

-- | Normalized mean vector.
meanVector :: [Vector Double] -> Vector Double
meanVector vecs = M.scale (1 / magnitude' mean) mean
  where mean = sum vecs / fromIntegral (length vecs)

-- | Divide matrix values by each row vector's magnitude.
normalizeMatrix :: Matrix Double -> Matrix Double
normalizeMatrix mat = DV.runSTMatrix $ do
  m <- DV.thawMatrix mat
  mapM_ (uncurry $ modifyRow m) $ zip [0..] magnitudes
  return m
    where magnitudes = map magnitude' $ M.toRows mat
          colSize = M.cols mat
          modifyRow mm i mag = mapM_ (flip (DV.modifyMatrix mm i) (/ mag)) [0 .. colSize - 1]

-- | Calculate magnitude of a 'Vector'.
magnitude :: Vector Double -> Double
magnitude = sqrt . M.sumElements . (** 2)

-- | Calculate magnitude of a 'Vector'. Return 1 if magnitude is 0 not to
-- divide by 0.
magnitude' :: Vector Double -> Double
magnitude' v = if mag == 0 then 1 else mag
  where mag = magnitude v

-- | Read sentences from a normalized file.
readSentences :: FilePath -> IO [Sentence]
readSentences path = do
  contents <- readFile path
  return $ map words $ lines contents

-- | Read a gigantic line and split it into sentences by each 1000 words.
readGiganticLine :: FilePath -> IO [Sentence]
readGiganticLine path = do
  contents <- readFile path
  return $ groupOf 100 $ words contents

groupOf :: Int -> [a] -> [[a]]
groupOf _ [] = []
groupOf n xs = (\(f, s) -> f : groupOf n s) $ splitAt n xs

saveI2W :: FilePath -> Index2Word -> IO ()
saveI2W = saveIntMap (\(i, w) -> show i ++ " " ++ w)

saveW2I :: FilePath -> Vocabulary -> IO ()
saveW2I = saveHashMap (\(w, (i, _)) -> w ++ " " ++ show i)

loadI2W :: FilePath -> IO Index2Word
loadI2W = loadIntMap (\(i:w:_) -> (read i, w))

loadW2I :: FilePath -> IO (HM.HashMap String Int)
loadW2I = loadHashMap (\(w:i:_) -> (w, read i))

max_stars_repo_name: shuhei/wv

max_issues_repo_name: shuhei/wv

max_forks_repo_name: shuhei/wv


text:
{-# language TypeFamilies #-}
module FFT where

import Data.Complex
import Data.Vector.Generic
import Numeric.FFT.Vector.Invertible

fft :: (Vector v a, a ~ Complex Double) => v a -> v a
fft = run dft

-- vector-fftw hace la normalizacion por si mismo
ifft :: (Vector v a, a ~ Complex Double) => v a -> v a
ifft = run idft

max_stars_repo_name: RubenAstudillo/proyecto-lena

max_issues_repo_name: RubenAstudillo/proyecto-lena

max_forks_repo_name: RubenAstudillo/proyecto-lena


text:
{-# LANGUAGE ViewPatterns #-}
-- |
-- Tests for quantile
module Tests.Quantile (tests) where

import Control.Exception
import qualified Data.Vector.Unboxed as U
import Test.Framework
import Test.Framework.Providers.HUnit
import Test.Framework.Providers.QuickCheck2
import Test.HUnit (Assertion,assertEqual,assertFailure)
import Test.QuickCheck hiding (sample)
import Numeric.MathFunctions.Comparison (ulpDelta,ulpDistance)
import Statistics.Quantile

tests :: Test
tests = testGroup "Quantiles"
  [ testCase "R alg. 4" $ compareWithR cadpw (0.00, 0.50, 2.50, 8.25, 10.00)
  , testCase "R alg. 5" $ compareWithR hazen (0.00, 1.00, 5.00, 9.00, 10.00)
  , testCase "R alg. 6" $ compareWithR spss  (0.00, 0.75, 5.00, 9.25, 10.00)
  , testCase "R alg. 7" $ compareWithR s     (0.000, 1.375, 5.000, 8.625,10.00)
  , testCase "R alg. 8" $ compareWithR medianUnbiased
      (0.0, 0.9166666666666667, 5.000000000000003, 9.083333333333334, 10.0)
  , testCase "R alg. 9" $ compareWithR normalUnbiased
      (0.0000, 0.9375, 5.0000, 9.0625, 10.0000)
  , testProperty "alg 7." propWeigtedAverage
    -- Test failures
  , testCase "weightedAvg should throw errors" $ do
      let xs  = U.fromList [1,2,3]
          xs0 = U.fromList []
      shouldError "Empty sample" $ weightedAvg 1 4 xs0
      shouldError "N=0"  $ weightedAvg 1 0 xs
      shouldError "N=1"  $ weightedAvg 1 1 xs
      shouldError "k<0"  $ weightedAvg (-1) 4 xs
      shouldError "k>N"  $ weightedAvg 5    4 xs
  , testCase "quantile should throw errors" $ do
      let xs  = U.fromList [1,2,3]
          xs0 = U.fromList []
      shouldError "Empty xs" $ quantile s 1 4 xs0
      shouldError "N=0"  $ quantile s 1 0 xs
      shouldError "N=1"  $ quantile s 1 1 xs
      shouldError "k<0"  $ quantile s (-1) 4 xs
      shouldError "k>N"  $ quantile s 5    4 xs
    --
  , testProperty "quantiles    are OK" propQuantiles
  , testProperty "quantilesVec are OK" propQuantilesVec
  ]

sample :: U.Vector Double
sample = U.fromList [0, 1, 2.5, 7.5, 9, 10]

-- Compare quantiles implementation with reference R implementation
compareWithR :: ContParam -> (Double,Double,Double,Double,Double) -> Assertion
compareWithR p (q0,q1,q2,q3,q4) = do
  assertEqual "Q 0" q0 $ quantile p 0 4 sample
  assertEqual "Q 1" q1 $ quantile p 1 4 sample
  assertEqual "Q 2" q2 $ quantile p 2 4 sample
  assertEqual "Q 3" q3 $ quantile p 3 4 sample
  assertEqual "Q 4" q4 $ quantile p 4 4 sample

propWeigtedAverage :: Positive Int -> Positive Int -> Property
propWeigtedAverage (Positive k) (Positive q) =
  (q >= 2 && k <= q) ==> let q1 = weightedAvg k q sample
                             q2 = quantile s k q sample
                         in counterexample ("weightedAvg   = " ++ show q1)
                          $ counterexample ("quantile      = " ++ show q2)
                          $ counterexample ("delta in ulps = " ++ show (ulpDelta q1 q2))
                          $ ulpDistance q1 q2 <= 16

propQuantiles :: Positive Int -> Int -> Int -> NonEmptyList Double -> Property
propQuantiles (Positive n)
              ((`mod` n) -> k1)
              ((`mod` n) -> k2)
              (NonEmpty xs)
  =   n >= 2
  ==> [x1,x2] == quantiles s [k1,k2] n rndXs
  where
    rndXs = U.fromList xs
    x1 = quantile s k1 n rndXs
    x2 = quantile s k2 n rndXs

propQuantilesVec :: Positive Int -> Int -> Int -> NonEmptyList Double -> Property
propQuantilesVec (Positive n)
                 ((`mod` n) -> k1)
                 ((`mod` n) -> k2)
                 (NonEmpty xs)
  =   n >= 2
  ==> U.fromList [x1,x2] == quantilesVec s (U.fromList [k1,k2]) n rndXs
  where
    rndXs = U.fromList xs
    x1 = quantile s k1 n rndXs
    x2 = quantile s k2 n rndXs

shouldError :: String -> a -> Assertion
shouldError nm x = do
  r <- try (evaluate x)
  case r of
    Left  (ErrorCall{}) -> return ()
    Right _             -> assertFailure ("Should call error: " ++ nm)

max_stars_repo_name: intricate/statistics

max_issues_repo_name: intricate/statistics

max_forks_repo_name: intricate/statistics


text:
{-# LANGUAGE FlexibleContexts, NamedFieldPuns #-}

module School.Unit.WeightDecay
( weightDecay ) where

import Conduit (ConduitM, mapMC)
import Numeric.LinearAlgebra (Container, Matrix, Vector,
                              cmap, scale, sumElements)
import School.Types.Slinky (Slinky(..), slinkyAppend)
import School.Unit.CostFunction (CostFunction(..), defSetupCost)
import School.Unit.CostParams (CostParams(..))
import School.Unit.UnitActivation (UnitActivation(..))
import School.Unit.UnitForward (ForwardStack)
import School.Unit.UnitGradient (UnitGradient(..))

square :: (Num a) => a -> a
square a = a * a

weight :: (Container Vector a, Num a)
       => a -> Matrix a -> a
weight coeff input =
  (*coeff) . sumElements $ cmap square input

weightDeriv :: (Container Vector a, Num a)
            => a -> Matrix a -> UnitGradient a
weightDeriv coeff input =
  BatchGradient . scale coeff $
    cmap (*2) input

errorMsg :: String
errorMsg = "Weight decay expects batch activation and no cost params"

compute :: (Container Vector a, Num a)
        => a
        -> UnitActivation a
        -> Slinky CostParams
        -> Either String a
compute coeff
        (BatchActivation input)
        (SNode NoCostParams _) =
  Right $ weight coeff input
compute coeff
        (BatchActivation input)
        SNil =
  Right $ weight coeff input
compute _ _ _ = Left errorMsg

deriv :: (Container Vector a, Num a)
      => a
      -> UnitActivation a
      -> Slinky CostParams
      -> Either String (UnitGradient a)
deriv coeff
      (BatchActivation input)
      (SNode NoCostParams _) =
  Right $ weightDeriv coeff input
deriv coeff
      (BatchActivation input)
      SNil =
  Right $ weightDeriv coeff input
deriv _ _ _ = Left errorMsg

prepare :: (Monad m)
        => ConduitM (ForwardStack a)
                    (ForwardStack a)
                    m
                    ()
prepare = mapMC $ \(activations, cParams) ->
  return (activations, slinkyAppend NoCostParams cParams)

weightDecay :: (Container Vector a, Num a, Monad m)
            => a
            -> CostFunction a m
weightDecay coeff =
  CostFunction { computeCost = compute coeff
               , derivCost = deriv coeff
               , prepareCost = prepare
               , setupCost = defSetupCost
               }


text:
module Main where

import Graphics.Rendering.Chart.Backend.Diagrams (toFile)
import qualified Graphics.Rendering.Chart.Easy as Plt
import Numeric.GSL.ODE
import Numeric.LinearAlgebra

dynamics :: Double -> [Double] -> [Double]
dynamics _ [_, x1, w] =
  [ x1,
    - mass * g + k * w + cd * x1 ** 2,
    (torque - c * w) / j
  ]
  where
    -- timestep
    dt = 0.01
    -- Coefficient of drag
    cd = 0.1
    -- mass
    mass = 22
    -- gravity
    g = 9.8
    -- inertia
    j = 0.05
    -- Omega to thrust
    k = 5
    -- Coefficient of friction motor
    c = 0.1
    -- Calculate a physics based torque
    torque =
      clamp (-10, 10) $
        mass * g * j / (k * dt)
          - x1 / dt ** 2
          - w * j / dt
          + c * w
    -- Clamp between lower bound and upper bound
    clamp (lb, ub) x = max lb (min ub x)

dumbCopter :: IO ()
dumbCopter = toFile Plt.def "DumpCopter.svg" $ do
  let totalTime = 5
      dt = 0.01
      ts = linspace (round $ totalTime / dt) (0, totalTime)
      sol = odeSolve dynamics [0, -10, 0] ts
      [pos, vel, omega] = toColumns sol
      ts' = toList ts
  Plt.plot (Plt.line "pos" [zip ts' (toList pos)])
  Plt.plot (Plt.line "vel" [zip ts' (toList vel)])
  Plt.plot (Plt.line "omega" [zip ts' (toList omega)])

main :: IO ()
main = dumbCopter


text:
{-# LANGUAGE BangPatterns, CPP #-}
import System.Environment
import Control.Monad
import Control.Seq
import Data.Complex
import System.IO
import Debug.Trace
import Control.DeepSeq
import Control.Monad.Par
import Control.Exception

import PortablePixmap
import Control.Monad.Par.AList as A

mandel :: Int -> Complex Double -> Int
mandel max_depth c = loop 0 0
  where
   fn = magnitude
   loop i !z
    | i == max_depth = i
    | fn(z) >= 2.0   = i
    | otherwise      = loop (i+1) (z*z + c)

threshold = 1

runMandel :: Double -> Double -> Double -> Double -> Int -> Int -> Int -> Par (AList [Int])
runMandel minX minY maxX maxY winX winY max_depth = do

  parBuildThreshM threshold (InclusiveRange 0 (winY-1)) $ \y ->
       do
          let l = [ mandelStep y x | x <- [0.. winX-1] ]
          deepseq l (return l)

  where
    mandelStep i j = mandel max_depth z
        where z = ((fromIntegral j * r_scale) / fromIntegral winY + minY) :+
                  ((fromIntegral i * c_scale) / fromIntegral winX + minX)

    r_scale  =  maxY - minY  :: Double
    c_scale =   maxX - minX  :: Double

makeImage :: Integer -> Integer -> Int -> AList [Int] -> PixMap
makeImage x y depth ls =
  createPixmap x y depth
   (map prettyRGB (concat (toList ls)))
 where
   prettyRGB :: Int -> (Int,Int,Int)
   prettyRGB s = let t = (depth - s) in (s,t,t)

simple x y depth =
  runMandel 0 0 x' y' x y depth
 where
  x' = fromIntegral x
  y' = fromIntegral y

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

-- A meaningless checksum.
mandelCheck :: AList [Int] -> Int -> Int -> Int
mandelCheck als max_col max_depth = loop 0 als 0
 where
 loop i als !sum | A.null als = sum
 loop i als !sum = loop (i+1) (A.tail als)
		        (loop2 i 0 (A.head als) sum)
 loop2 i j []    !sum = sum
 loop2 i j (h:t) !sum | h == max_depth = loop2 i (j+1) t (sum + i*max_col + j)
		      | otherwise      = loop2 i (j+1) t  sum

main = do args <- getArgs

          let (x,y,depth) =
		case args of
		 [] ->
    --                 runPar $ simple 3 3 3
		       (3,3,3)

		 [x,y,depth] ->
    --		   simple (read x) (read y) (read depth)
		       (read x, read y, read depth)

		 -- [minX,minY,maxX,maxY,winX,winY,depth] ->
		 --    runPar $
		 -- 	    runMandel (read minX) (read minY)
		 -- 		      (read maxX) (read maxY)
		 -- 		      (read winX) (read winY) (read depth)

          let ls = runPar$ simple x y depth
          when (False) $ do
	     hnd <- openFile "mandel_image.ppm" WriteMode
	     hSetBinaryMode hnd True
	     hPrint hnd (makeImage (fromIntegral x) (fromIntegral y) depth ls)
	     hClose hnd

          putStrLn$ "Spot check: " ++ show (mandelCheck ls y depth)


text:
-----------------------------------------------------------------------------
-- |
-- Module      :  Numeric.GSL.Vector
-- Copyright   :  (c) Alberto Ruiz 2007
-- License     :  GPL-style
--
-- Maintainer  :  Alberto Ruiz <aruiz@um.es>
-- Stability   :  provisional
-- Portability :  portable (uses FFI)
--
-- Low level interface to vector operations.
--
-----------------------------------------------------------------------------

module Numeric.GSL.Vector (
    sumF, sumR, sumQ, sumC,
    prodF, prodR, prodQ, prodC,
    dotF, dotR, dotQ, dotC,
    FunCodeS(..), toScalarR, toScalarF, toScalarC, toScalarQ,
    FunCodeV(..), vectorMapR, vectorMapC, vectorMapF, vectorMapQ,
    FunCodeSV(..), vectorMapValR, vectorMapValC, vectorMapValF, vectorMapValQ,
    FunCodeVV(..), vectorZipR, vectorZipC, vectorZipF, vectorZipQ,
    RandDist(..), randomVector
) where

import Data.Packed.Internal.Common
import Data.Packed.Internal.Signatures
import Data.Packed.Internal.Vector

import Data.Complex
import Foreign.Marshal.Alloc(free)
import Foreign.Marshal.Array(newArray)
import Foreign.Ptr(Ptr)
import Foreign.C.Types
import System.IO.Unsafe(unsafePerformIO)

fromei x = fromIntegral (fromEnum x) :: CInt

data FunCodeV = Sin
              | Cos
              | Tan
              | Abs
              | ASin
              | ACos
              | ATan
              | Sinh
              | Cosh
              | Tanh
              | ASinh
              | ACosh
              | ATanh
              | Exp
              | Log
              | Sign
              | Sqrt
              deriving Enum

data FunCodeSV = Scale
               | Recip
               | AddConstant
               | Negate
               | PowSV
               | PowVS
               deriving Enum

data FunCodeVV = Add
               | Sub
               | Mul
               | Div
               | Pow
               | ATan2
               deriving Enum

data FunCodeS = Norm2
              | AbsSum
              | MaxIdx
              | Max
              | MinIdx
              | Min
              deriving Enum

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

-- | sum of elements
sumF :: Vector Float -> Float
sumF x = unsafePerformIO $ do
           r <- createVector 1
           app2 c_sumF vec x vec r "sumF"
           return $ r @> 0

-- | sum of elements
sumR :: Vector Double -> Double
sumR x = unsafePerformIO $ do
           r <- createVector 1
           app2 c_sumR vec x vec r "sumR"
           return $ r @> 0

-- | sum of elements
sumQ :: Vector (Complex Float) -> Complex Float
sumQ x = unsafePerformIO $ do
           r <- createVector 1
           app2 c_sumQ vec x vec r "sumQ"
           return $ r @> 0

-- | sum of elements
sumC :: Vector (Complex Double) -> Complex Double
sumC x = unsafePerformIO $ do
           r <- createVector 1
           app2 c_sumC vec x vec r "sumC"
           return $ r @> 0

foreign import ccall unsafe "gsl-aux.h sumF" c_sumF :: TFF
foreign import ccall unsafe "gsl-aux.h sumR" c_sumR :: TVV
foreign import ccall unsafe "gsl-aux.h sumQ" c_sumQ :: TQVQV
foreign import ccall unsafe "gsl-aux.h sumC" c_sumC :: TCVCV

-- | product of elements
prodF :: Vector Float -> Float
prodF x = unsafePerformIO $ do
           r <- createVector 1
           app2 c_prodF vec x vec r "prodF"
           return $ r @> 0

-- | product of elements
prodR :: Vector Double -> Double
prodR x = unsafePerformIO $ do
           r <- createVector 1
           app2 c_prodR vec x vec r "prodR"
           return $ r @> 0

-- | product of elements
prodQ :: Vector (Complex Float) -> Complex Float
prodQ x = unsafePerformIO $ do
           r <- createVector 1
           app2 c_prodQ vec x vec r "prodQ"
           return $ r @> 0

-- | product of elements
prodC :: Vector (Complex Double) -> Complex Double
prodC x = unsafePerformIO $ do
           r <- createVector 1
           app2 c_prodC vec x vec r "prodC"
           return $ r @> 0

foreign import ccall unsafe "gsl-aux.h prodF" c_prodF :: TFF
foreign import ccall unsafe "gsl-aux.h prodR" c_prodR :: TVV
foreign import ccall unsafe "gsl-aux.h prodQ" c_prodQ :: TQVQV
foreign import ccall unsafe "gsl-aux.h prodC" c_prodC :: TCVCV

-- | dot product
dotF :: Vector Float -> Vector Float -> Float
dotF x y = unsafePerformIO $ do
           r <- createVector 1
           app3 c_dotF vec x vec y vec r "dotF"
           return $ r @> 0

-- | dot product
dotR :: Vector Double -> Vector Double -> Double
dotR x y = unsafePerformIO $ do
           r <- createVector 1
           app3 c_dotR vec x vec y vec r "dotR"
           return $ r @> 0

-- | dot product
dotQ :: Vector (Complex Float) -> Vector (Complex Float) -> Complex Float
dotQ x y = unsafePerformIO $ do
           r <- createVector 1
           app3 c_dotQ vec x vec y vec r "dotQ"
           return $ r @> 0

-- | dot product
dotC :: Vector (Complex Double) -> Vector (Complex Double) -> Complex Double
dotC x y = unsafePerformIO $ do
           r <- createVector 1
           app3 c_dotC vec x vec y vec r "dotC"
           return $ r @> 0

foreign import ccall unsafe "gsl-aux.h dotF" c_dotF :: TFFF
foreign import ccall unsafe "gsl-aux.h dotR" c_dotR :: TVVV
foreign import ccall unsafe "gsl-aux.h dotQ" c_dotQ :: TQVQVQV
foreign import ccall unsafe "gsl-aux.h dotC" c_dotC :: TCVCVCV

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

toScalarAux fun code v = unsafePerformIO $ do
    r <- createVector 1
    app2 (fun (fromei code)) vec v vec r "toScalarAux"
    return (r `at` 0)

vectorMapAux fun code v = unsafePerformIO $ do
    r <- createVector (dim v)
    app2 (fun (fromei code)) vec v vec r "vectorMapAux"
    return r

vectorMapValAux fun code val v = unsafePerformIO $ do
    r <- createVector (dim v)
    pval <- newArray [val]
    app2 (fun (fromei code) pval) vec v vec r "vectorMapValAux"
    free pval
    return r

vectorZipAux fun code u v = unsafePerformIO $ do
    r <- createVector (dim u)
    app3 (fun (fromei code)) vec u vec v vec r "vectorZipAux"
    return r

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

-- | obtains different functions of a vector: norm1, norm2, max, min, posmax, posmin, etc.
toScalarR :: FunCodeS -> Vector Double -> Double
toScalarR oper =  toScalarAux c_toScalarR (fromei oper)

foreign import ccall unsafe "gsl-aux.h toScalarR" c_toScalarR :: CInt -> TVV

-- | obtains different functions of a vector: norm1, norm2, max, min, posmax, posmin, etc.
toScalarF :: FunCodeS -> Vector Float -> Float
toScalarF oper =  toScalarAux c_toScalarF (fromei oper)

foreign import ccall unsafe "gsl-aux.h toScalarF" c_toScalarF :: CInt -> TFF

-- | obtains different functions of a vector: only norm1, norm2
toScalarC :: FunCodeS -> Vector (Complex Double) -> Double
toScalarC oper =  toScalarAux c_toScalarC (fromei oper)

foreign import ccall unsafe "gsl-aux.h toScalarC" c_toScalarC :: CInt -> TCVV

-- | obtains different functions of a vector: only norm1, norm2
toScalarQ :: FunCodeS -> Vector (Complex Float) -> Float
toScalarQ oper =  toScalarAux c_toScalarQ (fromei oper)

foreign import ccall unsafe "gsl-aux.h toScalarQ" c_toScalarQ :: CInt -> TQVF

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

-- | map of real vectors with given function
vectorMapR :: FunCodeV -> Vector Double -> Vector Double
vectorMapR = vectorMapAux c_vectorMapR

foreign import ccall unsafe "gsl-aux.h mapR" c_vectorMapR :: CInt -> TVV

-- | map of complex vectors with given function
vectorMapC :: FunCodeV -> Vector (Complex Double) -> Vector (Complex Double)
vectorMapC oper = vectorMapAux c_vectorMapC (fromei oper)

foreign import ccall unsafe "gsl-aux.h mapC" c_vectorMapC :: CInt -> TCVCV

-- | map of real vectors with given function
vectorMapF :: FunCodeV -> Vector Float -> Vector Float
vectorMapF = vectorMapAux c_vectorMapF

foreign import ccall unsafe "gsl-aux.h mapF" c_vectorMapF :: CInt -> TFF

-- | map of real vectors with given function
vectorMapQ :: FunCodeV -> Vector (Complex Float) -> Vector (Complex Float)
vectorMapQ = vectorMapAux c_vectorMapQ

foreign import ccall unsafe "gsl-aux.h mapQ" c_vectorMapQ :: CInt -> TQVQV

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

-- | map of real vectors with given function
vectorMapValR :: FunCodeSV -> Double -> Vector Double -> Vector Double
vectorMapValR oper = vectorMapValAux c_vectorMapValR (fromei oper)

foreign import ccall unsafe "gsl-aux.h mapValR" c_vectorMapValR :: CInt -> Ptr Double -> TVV

-- | map of complex vectors with given function
vectorMapValC :: FunCodeSV -> Complex Double -> Vector (Complex Double) -> Vector (Complex Double)
vectorMapValC = vectorMapValAux c_vectorMapValC

foreign import ccall unsafe "gsl-aux.h mapValC" c_vectorMapValC :: CInt -> Ptr (Complex Double) -> TCVCV

-- | map of real vectors with given function
vectorMapValF :: FunCodeSV -> Float -> Vector Float -> Vector Float
vectorMapValF oper = vectorMapValAux c_vectorMapValF (fromei oper)

foreign import ccall unsafe "gsl-aux.h mapValF" c_vectorMapValF :: CInt -> Ptr Float -> TFF

-- | map of complex vectors with given function
vectorMapValQ :: FunCodeSV -> Complex Float -> Vector (Complex Float) -> Vector (Complex Float)
vectorMapValQ oper = vectorMapValAux c_vectorMapValQ (fromei oper)

foreign import ccall unsafe "gsl-aux.h mapValQ" c_vectorMapValQ :: CInt -> Ptr (Complex Float) -> TQVQV

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

-- | elementwise operation on real vectors
vectorZipR :: FunCodeVV -> Vector Double -> Vector Double -> Vector Double
vectorZipR = vectorZipAux c_vectorZipR

foreign import ccall unsafe "gsl-aux.h zipR" c_vectorZipR :: CInt -> TVVV

-- | elementwise operation on complex vectors
vectorZipC :: FunCodeVV -> Vector (Complex Double) -> Vector (Complex Double) -> Vector (Complex Double)
vectorZipC = vectorZipAux c_vectorZipC

foreign import ccall unsafe "gsl-aux.h zipC" c_vectorZipC :: CInt -> TCVCVCV

-- | elementwise operation on real vectors
vectorZipF :: FunCodeVV -> Vector Float -> Vector Float -> Vector Float
vectorZipF = vectorZipAux c_vectorZipF

foreign import ccall unsafe "gsl-aux.h zipF" c_vectorZipF :: CInt -> TFFF

-- | elementwise operation on complex vectors
vectorZipQ :: FunCodeVV -> Vector (Complex Float) -> Vector (Complex Float) -> Vector (Complex Float)
vectorZipQ = vectorZipAux c_vectorZipQ

foreign import ccall unsafe "gsl-aux.h zipQ" c_vectorZipQ :: CInt -> TQVQVQV

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

data RandDist = Uniform  -- ^ uniform distribution in [0,1)
              | Gaussian -- ^ normal distribution with mean zero and standard deviation one
              deriving Enum

-- | Obtains a vector of pseudorandom elements from the the mt19937 generator in GSL, with a given seed. Use randomIO to get a random seed.
randomVector :: Int      -- ^ seed
             -> RandDist -- ^ distribution
             -> Int      -- ^ vector size
             -> Vector Double
randomVector seed dist n = unsafePerformIO $ do
    r <- createVector n
    app1 (c_random_vector (fi seed) ((fi.fromEnum) dist)) vec r "randomVector"
    return r

foreign import ccall unsafe "random_vector" c_random_vector :: CInt -> CInt -> TV

max_stars_repo_name: curiousleo/liquidhaskell

max_issues_repo_name: curiousleo/liquidhaskell

max_forks_repo_name: curiousleo/liquidhaskell


text:
{-# LANGUAGE CPP #-}
module Driver (
    E,

    Nat(..),
    Nat2(..),
    Pos(..),
    Pos2(..),
    Index(..),
    Index2(..),
    Assocs(..),
    Assocs2(..),
    BandedAssocs(..),

    module Control.Arrow,
    module Control.Monad,
    module Control.Monad.ST,

    AEq,
    (===),
    (~==),
    module Data.Function,
    module Data.Ix,
    module Data.List,
    module Data.Ord,

    module Debug.Trace,

    module Test.QuickCheck,

    module Text.Printf,

    field,

    module Test.Framework,
    module Test.Framework.Providers.QuickCheck2,
    ) where

import Control.Arrow
import Control.Monad
import Control.Monad.ST

import Data.AEq( AEq )
import qualified Data.AEq as AEq
import Data.Complex
import Data.Ix
import Data.Function
import Data.List
import Data.Ord

import Debug.Trace

import System.IO
import System.Random

import Test.Vector.Dense()
import Test.Matrix.Dense()
import Test.Matrix.Banded()

import Test.QuickCheck hiding ( vector )
import Test.QuickCheck.BLAS( Nat(..), Nat2(..), Pos(..), Pos2(..),
    Index(..), Index2(..), Assocs(..), Assocs2(..), BandedAssocs(..) )
import qualified Test.QuickCheck.BLAS as Test

import Test.Framework
import Test.Framework.Providers.QuickCheck2

import Text.Printf
import Text.Show.Functions

#ifdef COMPLEX
field = "Complex Double"
type E = Complex Double
#else
field = "Double"
type E = Double
#endif

infix 4 ===, ~==

x === y | (AEq.===) x y = True
        | otherwise = trace (printf "expected `%s', but got `%s'" (show y) (show x)) False

x ~== y | (AEq.~==) x y = True
        | otherwise = trace (printf "expected `%s', but got `%s'" (show y) (show x)) False

max_stars_repo_name: patperry/hs-linear-algebra

max_issues_repo_name: cartazio/hs-cblas

max_forks_repo_name: cartazio/hs-cblas


text:
import Statistics.Distribution
import Statistics.Distribution.Normal

nd :: NormalDistribution
nd = normalDistr 0.0 1.0

main = do
  putStrLn $ "P( 2< Z) =      " ++ show (complCumulative nd 2)
  putStrLn $ "P(-2<=Z<=2)=    " ++ show (cumulative nd 2 - cumulative nd (-2))
  putStrLn $ "P( 0<=Z<=1.73)= " ++ show (cumulative nd 1.73 - cumulative nd 0)

max_stars_repo_name: raydsameshima/Stat

max_issues_repo_name: raydsameshima/Stat

max_forks_repo_name: raydsameshima/Stat


text:
{-# LANGUAGE GADTs #-}
{-# LANGUAGE FlexibleContexts #-}

module HVX.Primitives
-- TODO(mh): Currently primitives that should generate implicit constraints are
-- unsupported. The should be supported shortily. (2014-06-04)
  ( apply
  , hadd
  , (+~)
  , hmul
  , (*~)
  , habs
  , neg
  , hlog
  , hexp
  , logsumexp
  , hmax
  , hmin
  , norm
  , berhu
  , huber
  , quadform
  , powBaseP0
  , powBaseP01
  , powBaseP1
  , powBaseP1InfEven
  , powBaseP1InfNotInt
  ) where

import Numeric.LinearAlgebra hiding (i)
import Numeric.LinearAlgebra.LAPACK

import HVX.Internal.DCP
import HVX.Internal.Matrix
import HVX.Internal.Primitives
import HVX.Internal.Util

apply :: (Vex vf, Mon mf, Vex ve, Mon me, Vex (ApplyVex vf mf ve me), Mon (ApplyMon mf me))
  => Fun vf mf -> Expr ve me -> Expr (ApplyVex vf mf ve me) (ApplyMon mf me)
apply = EFun

hadd :: (Vex v1, Mon m1, Vex v2, Mon m2, Vex (AddVex v1 v2), Mon (AddMon m1 m2))
  => Expr v1 m1 -> Expr v2 m2 -> Expr (AddVex v1 v2) (AddMon m1 m2)
hadd = EAdd

infixl 6 +~
(+~) :: (Vex v1, Mon m1, Vex v2, Mon m2, Vex (AddVex v1 v2), Mon (AddMon m1 m2))
  => Expr v1 m1 -> Expr v2 m2 -> Expr (AddVex v1 v2) (AddMon m1 m2)
(+~) = hadd

-- Constructors that enforce DCP constraints.
hmul :: (Vex v1, Mon m1, Vex (ApplyVex Affine Nonmon v1 m1), Mon (ApplyMon Nonmon m1))
  => Expr Affine Const -> Expr v1 m1 -> Expr (ApplyVex Affine Nonmon v1 m1) (ApplyMon Nonmon m1)
hmul (EConst a) e = apply (Mul a) e
hmul _ _ = error "the left argument of a multiply must be a constant"

infixl 7 *~
(*~) :: (Vex v1, Mon m1, Vex (ApplyVex Affine Nonmon v1 m1), Mon (ApplyMon Nonmon m1))
  => Expr Affine Const -> Expr v1 m1 -> Expr (ApplyVex Affine Nonmon v1 m1) (ApplyMon Nonmon m1)
(*~) (EConst a) e = apply (Mul a) e
(*~) _ _ = error "the left argument of a multiply must be a constant"

habs :: (Vex v1, Mon m1, Vex (ApplyVex Convex Nonmon v1 m1), Mon (ApplyMon Nonmon m1))
  => Expr v1 m1 -> Expr (ApplyVex Convex Nonmon v1 m1) (ApplyMon Nonmon m1)
habs = apply Abs

neg :: (Vex v1, Mon m1, Vex (ApplyVex Affine Noninc v1 m1), Mon (ApplyMon Noninc m1))
  => Expr v1 m1 -> Expr (ApplyVex Affine Noninc v1 m1) (ApplyMon Noninc m1)
neg = apply Neg

hlog :: (Vex v1, Mon m1, Vex (ApplyVex Concave Nondec v1 m1), Mon (ApplyMon Nondec m1))
  => Expr v1 m1 -> Expr (ApplyVex Concave Nondec v1 m1) (ApplyMon Nondec m1)
hlog = apply Log

hexp :: (Vex v1, Mon m1, Vex (ApplyVex Convex Nondec v1 m1), Mon (ApplyMon Nondec m1))
  => Expr v1 m1 -> Expr (ApplyVex Convex Nondec v1 m1) (ApplyMon Nondec m1)
hexp = apply Exp

logsumexp :: (Vex v1, Mon m1, Vex (ApplyVex Convex Nondec v1 m1), Mon (ApplyMon Nondec m1))
  => Expr v1 m1 -> Expr (ApplyVex Convex Nondec v1 m1) (ApplyMon Nondec m1)
logsumexp = apply LogSumExp

hmax :: (Vex v1, Mon m1, Vex (ApplyVex Convex Nondec v1 m1), Mon (ApplyMon Nondec m1))
  => Expr v1 m1 -> Expr (ApplyVex Convex Nondec v1 m1) (ApplyMon Nondec m1)
hmax = apply Max

hmin :: (Vex v1, Mon m1, Vex (ApplyVex Concave Nondec v1 m1), Mon (ApplyMon Nondec m1))
  => Expr v1 m1 -> Expr (ApplyVex Concave Nondec v1 m1) (ApplyMon Nondec m1)
hmin = apply Min

norm :: (Vex v1, Mon m1, Vex (ApplyVex Convex Nonmon v1 m1), Mon (ApplyMon Nonmon m1))
  => Double -> Expr v1 m1 -> Expr (ApplyVex Convex Nonmon v1 m1) (ApplyMon Nonmon m1)
norm p
  | p == infinity = error "Internal: Infinity norm should become max . abs."
  | 1 <= p = apply (Norm p)
  | otherwise = error "Internal: Norm only supports p >= 1."

berhu :: (Vex v1, Mon m1, Vex (ApplyVex Convex Nonmon v1 m1), Mon (ApplyMon Nonmon m1))
  => Double -> Expr v1 m1 -> Expr (ApplyVex Convex Nonmon v1 m1) (ApplyMon Nonmon m1)
berhu m
  | 0 < m = apply (Berhu m)
  | otherwise = error "Internal: Berhu only supports m >= 0."

huber :: (Vex v1, Mon m1, Vex (ApplyVex Convex Nonmon v1 m1), Mon (ApplyMon Nonmon m1))
  => Double -> Expr v1 m1 -> Expr (ApplyVex Convex Nonmon v1 m1) (ApplyMon Nonmon m1)
huber m
  | 0 < m = apply (Huber m)
  | otherwise = error "Internal: Huber only supports m >= 0."

quadform :: (Mon m1, Vex (ApplyVex Convex Nonmon Affine m1), Mon (ApplyMon Nonmon m1))
  => Expr Affine Const -> Expr Affine m1 -> Expr (ApplyVex Convex Nonmon Affine m1) (ApplyMon Nonmon m1)
quadform (EConst a) e
  | rows a == cols a
    && fpequalsMat a (trans a)
    && 0 <= maxElement (eigOnlyS a) = apply (Quadform a) e
  | otherwise = error "Matrices in quadratic forms must be positive semidefinite."
quadform _ _ = error "The matrix sandwitched by the quadratic form must be constant."

powBaseP0 :: (Vex v1, Mon m1, Vex (ApplyVex Affine Const v1 m1), Mon (ApplyMon Const m1))
  => Double -> Expr v1 m1 -> Expr (ApplyVex Affine Const v1 m1) (ApplyMon Const m1)
powBaseP0 p
  | p `fpequals` 0 = apply PowBaseP0
  | otherwise = error "Internal: PowBaseP0 only supports p == 0."

powBaseP01 :: (Vex v1, Mon m1, Vex (ApplyVex Concave Nondec v1 m1), Mon (ApplyMon Nondec m1))
  => Double -> Expr v1 m1 -> Expr (ApplyVex Concave Nondec v1 m1) (ApplyMon Nondec m1)
powBaseP01 p
  | 0 < p && p < 1 = apply (PowBaseP01 p)
  | otherwise = error "Internal: PowBaseP01 only supports 0 < p < 1."

powBaseP1 :: (Vex v1, Mon m1, Vex (ApplyVex Affine Nondec v1 m1), Mon (ApplyMon Nondec m1))
  => Double -> Expr v1 m1 -> Expr (ApplyVex Affine Nondec v1 m1) (ApplyMon Nondec m1)
powBaseP1 p
  | p `fpequals` 1 = apply PowBaseP1
  | otherwise = error "Internal: PowBaseP1 only supports p == 1."

powBaseP1InfEven :: (Vex v1, Mon m1, Vex (ApplyVex Convex Nonmon v1 m1), Mon (ApplyMon Nonmon m1))
  => Double -> Expr v1 m1 -> Expr (ApplyVex Convex Nonmon v1 m1) (ApplyMon Nonmon m1)
powBaseP1InfEven p
  | 1 < p && isInteger p && even intP = apply (PowBaseP1InfEven intP)
  | otherwise = error "Internal: PowBaseP1InfEven only supports even p > 1."
  where intP = round p :: Integer

powBaseP1InfNotInt :: (Vex v1, Mon m1, Vex (ApplyVex Convex Nondec v1 m1), Mon (ApplyMon Nondec m1))
  => Double -> Expr v1 m1 -> Expr (ApplyVex Convex Nondec v1 m1) (ApplyMon Nondec m1)
powBaseP1InfNotInt p
  | 1 < p && not (isInteger p) = apply (PowBaseP1InfNotInt p)
  | otherwise = error "Internal: PowBaseP1InfNotInt only supports non integral p > 1."

max_stars_repo_name: wellposed/hvx

max_issues_repo_name: wellposed/hvx

max_forks_repo_name: wellposed/hvx


text:
{-# LANGUAGE NoMonomorphismRestriction,
             ScopedTypeVariables
             #-}

module PFP where

import Debug.Trace

import Data.Traversable
import Control.Applicative
import Control.Monad

import Control.Monad.Random
import Control.Monad.State
import Control.Monad.Fix

import Numeric.SpecFunctions -- From cabal install spec-functions

--import Rec
import SplitEval
import Mem

bernoulli :: (Functor m, MonadRandom m) => Double -> m Bool
bernoulli p = (< p) <$> getRandom

-- From Alexey Radul's HSVenture
beta :: (MonadRandom m) => Double -> Double -> m Double
beta alpha beta = do
  -- Adapted from Statistics.Distribution.Beta; not reused because of
  -- funny randomness management convention.
  x <- getRandomR (0.0,1.0)
  return $ quantile x
    where
      quantile x | x == 0 = 0
                 | x == 1 = 1
                 | 0 < x && x < 1 = invIncompleteBeta alpha beta x
                 | otherwise = error $ "x must be in the range [0,1], got: " ++ show x

-- Makes a CRP sampler given the concentration parameter alpha via stick breaking.
makeCRP ::
  (RandomGen g, Functor m, MonadSplit g m, Functor f, MonadRandom f) =>
  Double -> m (f Int)
makeCRP alpha = do
  -- This is probably our only choice for this type
  -- without adding Proxy arguments. Oh well.
  let breakStick :: (RandomGen g) => Rand g Double
      breakStick = beta 1 alpha

  sticks <- repeatS breakStick

  let choose (p:ps) = do
        b <- bernoulli p
        if b
          then return 0
          else (1 +) <$> choose ps

  return $ choose sticks

-- This type was inferred by ghc, so should be the most general.
dpmem ::
  (MonadEval g n, MonadRandom f, MonadSplit g m, RandomGen g,
      Applicative n, Functor f, Functor m) =>
  Double -> n a -> m (f a)
dpmem alpha sampler = do
  crp <- makeCRP alpha
  tables <- memIntS (const sampler)
  return $ tables <$> crp

dpMemInt :: forall a g n m f.
  (MonadEval g n, MonadRandom f, MonadSplit g m, RandomGen g,
      Applicative n, Functor f, Functor m) =>
  Double -> (Int -> n a) -> m (Int -> f a)
dpMemInt alpha f = memIntS f'
  where -- This type needs to have both (MonadSplit g m) for dpmem
        -- and (MonadEval g n) for memIntS
        f' :: Int -> Rand g (f a)
        f' = (dpmem alpha) . f

dpBernoulli = dpmem 1.0 $ (bernoulli 0.5 :: (RandomGen g) => Rand g Bool)

randomWalkRec ::
  (Functor m, MonadRandom m) => Double -> (Int -> m Int) -> Int -> m Int
randomWalkRec _ _ 0 = return 0
randomWalkRec p f n = do
  x <- f (n-1)
  b <- bernoulli p
  let dx = if b then 1 else -1
  return $ x + dx

randomWalk
  :: (MonadSplit g m, RandomGen g, MonadFix m, Functor m) =>
     Double -> m (Int -> Int)
randomWalk p = mfix $ memIntS . f
  where f :: RandomGen g => (Int -> Int) -> (Int -> Rand g Int)
        f = (randomWalkRec p) . (return .)

dpRandomWalk :: forall g n m f.
  (MonadEval g n, MonadRandom n, MonadSplit g m, RandomGen g,
    Applicative n, MonadFix m, Functor m) =>
   Double -> Double -> m (Int -> n Int)
dpRandomWalk alpha p = mfix $ (dpMemInt alpha) . (randomWalkRec p)

dpRandomWalkIO alpha p = (evalRandIO .) <$> dpRandomWalk alpha p


text:
import Numeric.LinearAlgebra
import Graphics.Plot
import Text.Printf(printf)

expand :: Int -> Vector Double -> Matrix Double
expand n x = fromColumns $ map (x^) [0 .. n]

polynomialModel :: Vector Double -> Vector Double -> Int
                -> (Vector Double -> Vector Double)
polynomialModel x y n = f where
    f z = expand n z <> ws
    ws  = expand n x <\> y

main = do
    [x,y] <- (toColumns . readMatrix) `fmap` readFile "data.txt"
    let pol = polynomialModel x y
    let view = [x, y, pol 1 x, pol 2 x, pol 3 x]
    putStrLn $ "  x      y      p 1    p 2    p 3"
    putStrLn $ format "  " (printf "%.2f") $ fromColumns view
    mplot view


text:
-----------------------------------------------------------------------------
-- |
-- Module     : Numeric.LinearAlgebra.Matrix.Cholesky
-- Copyright  : Copyright (c) 2010, Patrick Perry <patperry@gmail.com>
-- License    : BSD3
-- Maintainer : Patrick Perry <patperry@gmail.com>
-- Stability  : experimental
--
-- Cholesky factorizations of symmetric (Hermitian) positive-definite
-- matrices.
--

module Numeric.LinearAlgebra.Matrix.Cholesky (
    -- * Immutable interface
    cholFactor,
    cholSolveVector,
    cholSolveMatrix,

    -- * Mutable interface
    cholFactorM,
    cholSolveVectorM_,
    cholSolveMatrixM_,

    ) where

import Control.Monad( when )
import Control.Monad.ST( ST, runST, unsafeIOToST )
import Text.Printf( printf )

import qualified Foreign.LAPACK as LAPACK

import Numeric.LinearAlgebra.Types

import Numeric.LinearAlgebra.Matrix.Base( Matrix )
import Numeric.LinearAlgebra.Matrix.STBase( RMatrix, STMatrix )
import qualified Numeric.LinearAlgebra.Matrix.STBase as M

import Numeric.LinearAlgebra.Vector( Vector, STVector )
import qualified Numeric.LinearAlgebra.Vector as V

-- | @cholFactor a@ tries to compute the Cholesky
-- factorization of @a@.  If @a@ is positive-definite then the routine
-- returns @Right@ with the factorization.  If the leading minor of order @i@
-- is not positive-definite then the routine returns @Left i@.
cholFactor :: (LAPACK e)
                 => Herm Matrix e
                 -> Either Int (Chol Matrix e)
cholFactor (Herm uplo a) = runST $ do
    ma <- M.newCopy a
    cholFactorM (Herm uplo ma)
        >>= either (return . Left) (\(Chol uplo' ma') -> do
                a' <- M.unsafeFreeze ma'
                return $ Right (Chol uplo' a')
                )

-- | @cholSolveVector a x@ returns @a \\ x@.
cholSolveVector :: (LAPACK e)
                      => Chol Matrix e
                      -> Vector e
                      -> Vector e
cholSolveVector a x = V.create $ do
    x' <- V.newCopy x
    cholSolveVectorM_ a x'
    return x'

-- | @cholSolveMatrix a b@ returns @a \\ b@.
cholSolveMatrix :: (LAPACK e)
                => Chol Matrix e
                -> Matrix e
                -> Matrix e
cholSolveMatrix a c = M.create $ do
    c' <- M.newCopy c
    cholSolveMatrixM_ a c'
    return c'

-- | @cholFactorM a@ tries to compute the Cholesky
-- factorization of @a@ in place.  If @a@ is positive-definite then the
-- routine returns @Right@ with the factorization, stored in the same
-- memory as @a@.  If the leading minor of order @i@ is not
-- positive-definite then the routine returns @Left i@.
-- In either case, the original storage of @a@ is destroyed.
cholFactorM :: (LAPACK e)
            => Herm (STMatrix s) e
            -> ST s (Either Int (Chol (STMatrix s) e))
cholFactorM (Herm uplo a) = do
    (ma,na) <- M.getDim a
    let n = na

    when (not $ (ma,na) == (n,n)) $ error $
        printf ("cholFactorM"
                ++ " (Herm _ <matrix with dim (%d,%d)>): nonsquare matrix"
               ) ma na

    unsafeIOToST $
        M.unsafeWith a $ \pa lda -> do
            info <- LAPACK.potrf uplo n pa lda
            return $ if info > 0 then Left info
                                 else Right (Chol uplo a)

-- | @cholSolveVectorM_ a x@ sets @x := a \\ x@.
cholSolveVectorM_ :: (LAPACK e, RMatrix m)
                  => Chol m e
                  -> STVector s e
                  -> ST s ()
cholSolveVectorM_ a x =
    M.withFromColM x $ \x' ->
        cholSolveMatrixM_ a x'

-- | @cholSolveMatrixM_ a b@ sets @b := a \\ b@.
cholSolveMatrixM_ :: (LAPACK e, RMatrix m)
                  => Chol m e
                  -> STMatrix s e
                  -> ST s ()
cholSolveMatrixM_ (Chol uplo a) b = do
    (ma,na) <- M.getDim a
    (mb,nb) <- M.getDim b
    let (n,nrhs) = (mb,nb)

    when ((not . and) [ (ma,na) == (n,n)
                      , (mb,nb) == (n,nrhs)
                      ]) $ error $
        printf ("cholSolveMatrixM_"
                ++ " (Chol _ <matrix with dim (%d,%d)>)"
                ++ " <matrix with dim (%d,%d)>"
                ++ ": dimension mismatch")
               ma na mb nb

    unsafeIOToST $
        M.unsafeWith a $ \pa lda ->
        M.unsafeWith b $ \pb ldb ->
            LAPACK.potrs uplo n nrhs pa lda pb ldb


text:
module Math.Probably.Circular where

import Math.Probably.FoldingStats
import Numeric.LinearAlgebra
import Control.Applicative

--approximate method first mentioned http://en.wikipedia.org/wiki/Mean_of_circular_quantities
circularMean :: Fold Double Double
circularMean = before (after vmean getAngle) toCartesian where
  vmean = pure (scale) <*> after realLengthF recip <*> sumF
  getAngle v = atan2 (v@>1) (v@>0)
  toCartesian alpha = fromList [cos alpha, sin alpha]

--http://webspace.ship.edu/pgmarr/Geo441/Lectures/Lec%2016%20-%20Directional%20Statistics.pdf
circularDispersion :: Fold Double Double
circularDispersion = after (both xs ys) getR where
  xs = before sumF cos
  ys = before sumF sin
  getR (x,y) = sqrt (x*x + y*y)

--example for running statistic

runCircularMean :: [Double] -> Double
runCircularMean = runStat circularMean


text:
-- (C) Copyright Chris Banks 2011-2012

-- This file is part of The Continuous Pi-calculus Workbench (BondCalculusWB).

--     BondCalculusWB is free software: you can redistribute it and/or modify
--     it under the terms of the GNU General Public License as published by
--     the Free Software Foundation, either version 3 of the License, or
--     (at your option) any later version.

--     BondCalculusWB is distributed in the hope that it will be useful,
--     but WITHOUT ANY WARRANTY; without even the implied warranty of
--     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
--     GNU General Public License for more details.

--     You should have received a copy of the GNU General Public License
--     along with BondCalculusWB.  If not, see <http://www.gnu.org/licenses/>.

{-# OPTIONS_GHC -fno-warn-incomplete-patterns -w #-}

module CPi.Tests where

import CPi.Lib
import CPi.Parser
import CPi.Semantics
import CPi.ODE
import CPi.Logic
import CPi.Plot
import CPi.Matlab

import Text.ParserCombinators.Parsec
import System.IO
import Data.List as L
import qualified Data.Map as Map

import qualified Numeric.GSL as GSL
import qualified Numeric.LinearAlgebra as LA
import qualified Graphics.Plot as Plot

-----------------------------------
-- Parser tests:
-----------------------------------
{-
-- Parser Test harnesses:
tParse :: (Pretty a) => Parser a -> String -> IO ()
tParse p input = case (parse p "" input) of
                      Left err -> do putStr "parse error at";
                                     print err
                      Right x -> print (pretty x)

tParse' :: (Show a) => Parser a -> String -> IO ()
tParse' p input = case (parse p "" input) of
                      Left err -> do putStr "parse error at";
                                     print err
                      Right x -> print x

-- -- Test file input:

tFile x = do f <- readFile x
             tParse' pDefinitionLines f
-}
-- Load an Env for testing:
tEnv x = do f <- readFile x
            case parseFile f of
              Left err -> return []
              Right x -> return x

tProc env x = maybe tEmptyProc id (lookupProcName env x)

tEmptyProc = Process [] (AffNet [])

-- tSpec x = case (parse pSpecies "" x) of
--             Left err -> error $ show err
--             Right x -> return x

-------------------------------
-- Tests for transition system:
-------------------------------

-- Test transition system for enzyme example:
{-tTrans = do file <- readFile "testEnzyme.cpi"
            let defns = (\(Right x) -> x)(parse pDefinitionLines "" file)
            putStrLn $ concat $ map (\x->(pretty x)++"\n") defns
            -- transitions of the species individually:
            let mts = trans defns (MTS []) (Def "S" ["s"])
            let mts' = trans defns mts (Def "P" [])
            let mts'' = trans defns mts' (Def "E" ["e"])
            -- let mts''' = trans defns mts'' (Def "I" ["i"])
            putStrLn "Initial species transtions:\n"
            putStrLn $ pretty mts''
            -- find resultant complexes (pseudoapps):
            let compxs = appls mts''
            putStrLn "Complexes formed (pseudoapplications):\n"
            putStrLn $ concat $ map (\x->(pretty x)++"\n") compxs
            -- find transitions for the complexes and add to MTS
            let finalmts = transs defns mts'' compxs
            putStrLn "Multi-transition system:\n"
            putStrLn $ pretty finalmts
            -- calculate the closure of the MTS
            let fixedmts = fixMTS defns finalmts
            putStrLn "Closed MTS:\n"
            putStrLn $ pretty fixedmts
-}

-- Some constants for playing with the Enzyme example:
tEnzDefs = do file <- readFile "testEnzyme.cpi"
              return ((\(Right x) -> x) (parseFile file))

tEnzPi = Process [(Def "S" ["s"],1.0),(Def "E" ["e"],0.1),(Def "P" [],0.0)] (AffNet [Aff (("e","s"),1.0)])

tEnzPi' = Process [(Def "S" ["s"],1.0),(Def "E" ["e"],0.5),(Def "P" [],0.0),((New (AffNet [Aff (("a","t"),1.0),Aff (("a","u"),0.5)]) (Par [Sum [(Comm "a" [] [],Def "E" ["e"])],Sum [(Comm "t" [] [],Def "P" []),(Comm "u" [] [],Def "S" ["s"])]])),0.0)] (AffNet [Aff (("e","s"),1.0)])

{-
-- Test get full MTS of a process:
tPTrans = do file <- readFile "testEnzyme.cpi"
             let defns = (\(Right x) -> x)(parse pDefinitionLines "" file)
             let pi = Process [(Def "S" ["s"],"1.0"),(Def "E" ["e"],"0.1"),(Def "P" [],"0.0")] (AffNet [Aff (("e","s"),"1.0")])
             let mts = processMTS defns pi
             putStrLn $ (prettys defns)++"\nTransitions:\n"++(pretty mts)

-- test recursive species
tTransRec = do let defns = (\(Right x)->x)(parse pDefinitionLines "" "species P() = tau<1>.P();")
               let mts = trans defns (MTS []) (Def "P" [])
               putStrLn $ pretty mts

-- test infinite species
tTransInf = do let defns = (\(Right x)->x)(parse pDefinitionLines "" "species P() = tau<1>.(P()|P());")
               let mts = trans defns (MTS []) (Def "P" [])
               putStrLn $ pretty mts

-- Test trans of Def/Nil
tTrans' = trans [tcSpecP0] (MTS []) tcP

-- Test trans of singleton Sum of Tau.P()
tTrans'2 = trans [tcSpecP0] (MTS []) tcSum1TauP

-- Test trans of Sum of Tau.P() + Tau.Q()
tTrans'3 = trans [tcSpecP0,tcSpecQ0] (MTS []) tcSum2TauPQ

tLookupDef = lookupDef [tcSpecP0] (tcP)

-- Test tensor:
tTensor = do env <- tEnzDefs
             let net = AffNet [Aff (("e","s"),"1.0")]
             let e = Process [(Def "E" ["e"],"0.1")] net
             let s = Process [(Def "S" ["s"],"1.0")] net
             print $ tensor env net (partial env e) (partial env s)
-}
-- Test tensor':
{-tTensor' = do env <- tEnzDefs
              let net = AffNet [Aff (("e","s"),"1.0")]
              let e = Process [(Def "E" ["e"],"0.1")] net
              let s = Process [(Def "S" ["s"],"1.0")] net
              putStrLn $ prettyODE env $ tensor' env net (partial' env e) (partial' env s)
-}

-- Test symbolic dPdt:
{-tdPdt = do env <- tEnzDefs
           let pi = tEnzPi'
           putStrLn $ prettyODE env $ dPdt' env pi
-}

------------------
-- Test constants:
------------------

-- tau@<0.5>.P():
tcSum1TauP = Sum [((Tau 0.5),tcP)]
--  tau@<0.5>.P() + tau@<0.5>.Q()
tcSum2TauPQ = Sum [((Tau 0.5),tcP),((Tau 0.6),tcQ)]
-- P()
tcP = Def "P" []
-- Q()
tcQ = Def "Q" []
-- species P() = 0
tcSpecP0 = SpeciesDef "P" [] Nil
-- species Q() = 0
tcSpecQ0 = SpeciesDef "Q" [] Nil

-- tests for some struct.cong. rules
tcSs = Def "S" ["s"]
tcSsa = Def "S" ["s","a"]
tcXx = Def "X" ["x"]
tcNet0 = AffNet []
tcNet1 = AffNet [Aff (("s","s'"),1)]
tcNet2 = AffNet [Aff (("a","b"),1)]
tcNet3 = AffNet [Aff (("s","s'"),1),Aff (("x","y"),1)]
tcNNS = New tcNet1 (New tcNet2 tcSs)
tcNNSsa = New tcNet1 (New tcNet2 tcSsa)
tcNSX = New tcNet1 (Par [tcXx,tcSs])
tcN3SX = New tcNet3 (Par [tcXx,tcSs])
tcN0SX = New tcNet0 (Par [tcXx,tcSs])

tcC1 = ConcBase (Par [tcXx,tcSs]) ["a"] ["s"]
tcC2 = ConcPar tcC1 [tcQ,tcP]

tcNestPar = Par [tcQ, tcP,Par [tcXx,Par [tcSs,Nil]]]

-----------------------------------
-- Tests for new fixmts
-----------------------------------
{-
tNPs = do env <- tEnv "models/ddos.cpi"
          let pi = tProc env "Pi"
              net' (Process _ net) = net
              net = net' pi
              ss' (Process ss _) = ss
              ss = ss' pi
              initmts = transs env (MTS []) (map fst ss)
          return $ newPrimes env initmts

tFixmts = do env <- tEnv "models/ddos.cpi"
             let pi = tProc env "Pi"
             return $ processMTS env pi
-}
-----------------------------------
-- ODE solver tests
-----------------------------------

-- xdot t [x,v] = [v, -0.95*x - 0.1*v]
txdot t [e,s,p,c] = [(-1.0)*1.1*e*s + (-1.0)*0.9*c + 0.5*c,
                    (-1.0)*1.1*e*s + 0.9*c,
                    (-1.0)*0.5*p + 0.5*c,
                    1.1*e*s + (-1.0)*0.9*c + (-1.0)*0.5*c]
txdot _ _ = undefined
ts = LA.linspace 250 (0,25)
-- sol = GSL.odeSolve xdot [10,0] ts
sol = GSL.odeSolve txdot [0.4,2.3,0,0] ts
tODE = Plot.mplot (ts : LA.toColumns sol)

{-tXdot = do env <- tEnzDefs
           let pi = tEnzPi'
           let x = xdot env (dPdt' env pi)
           -- return $ x ts [1.0,0,0.5,0]
           let sol' = GSL.odeSolve x [1.0,0,0.5,0] ts
           Plot.mplot (ts : LA.toColumns sol')
-}
{-
tSeries = do env <- tEnv "testEnzyme.cpi"
             let pi = tProc env "Pi"
             let mts = processMTS env pi
             let pi' = wholeProc env pi mts
             let dpdt = dPdt' env mts pi'
             let odes = xdot env dpdt
             let inits = initials env pi' dpdt
             let ts = timePoints 250 (0,25)
             let soln = solveODE env pi' dpdt (250,(0,25))
             let ss = speciesIn env dpdt
             return $ timeSeries ts soln ss
-}

-------------------------
-- Model checker tests:
-------------------------
{-
tModelCheck src p f trc = do env <- tEnv src
                             let pi = tProc env p
                             return $ modelCheck env solveODE trc pi mcts f

tModelCheckDP src p f trc = do env <- tEnv src
                               let pi = tProc env p
                               return $ modelCheckDP env solveODE trc pi mcts f

tModelCheckHy src p f trc = do env <- tEnv src
                               let pi = tProc env p
                               return $ modelCheckHy env solveODE trc pi mcts f

tModelCheckHy2 src p f trc = do env <- tEnv src
                                let pi = tProc env p
                                return $ modelCheckHy2 env solveODE trc pi mcts f
-}
mcts = (1000,(0,100))
--
-- Some contrived formulae for benchmarking:
-- (Use on the testGT.cpi model)

-- F(S<0.1)
tF1 = Pos (0,10) (ValLT (Conc (Def "S" ["s"])) (R 0.1))

-- F{5}(S<0.4 AND (F{5} (S<0.1)))
tF1b = Pos (0,5)
       (Conj
        (ValLT (Conc (Def "S" ["s"])) (R 0.4))
        (Pos (0,5)
         (ValLT (Conc (Def "S" ["s"])) (R 0.1))))

tF1c = Pos (0,10) (Nec (0,5)
                  (ValLT (Conc (Def "S" ["s"])) (R 0.1)))

-- test gaurantee (introduce an inhibitor)
tF2 = Gtee "In" (Neg tF1)

-- G(E>0.001) enzyme never runs out
tF3 = Nec (0,25) (ValGT (Conc (Def "E" ["e"])) (R 0.001))
tFn3 = Neg tF3

-- still true with inhibitor:
tF4 = Gtee "In" tF3

-- test nested guarantee (re-solves for every time-point):
tF5 = Nec (0,25) (Gtee "In" tF3)

-- test nested TL
-- G(F(S<0.1))
tF6 = Nec (0,25) tF1
-- G(F(G(E>0.001)))
tF7 = Nec (0,25) $ Pos (0,25) tF3
-- G(F(G(Inhib|>(E>0.001))))
tF8 = Nec (0,25) $ Pos (0,25) tF4

-- should be faster in Hy than in DP:
tF9 = Pos (0,20) (ValGT (Conc (Def "P" [])) (R 0.5))

tF10 = Nec (0,25) tF9

-- for checking sim times:
phi = (ValGT (Conc (Def "A" [])) (R 0.5))
psi = (ValGT (Conc (Def "B" [])) (R 0.5))
tF11 = Pos (0,10) (Nec (0,20) phi)
tF12 = Nec (0,10) (Conj (Pos (0,10) phi) (Pos (0,20) psi))
tF13 = Gtee "I" (Pos (0,10) phi)

tFR1 = (Disj (Until (0.0,3.3935742971887564) T (ValLT (Conc (Def "S" ["s"])) (R 0.1))) (Disj (Until (0.0,3.4939759036144586) T (ValLT (Conc (Def "S" ["s"])) (R 0.1))) (Disj (Until (0.0,3.594377510040161) T (ValLT (Conc (Def "S" ["s"])) (R 0.1))) (Disj (Until (0.0,3.694779116465864) T (ValLT (Conc (Def "S" ["s"])) (R 0.1))) (Disj (Until (0.0,3.7951807228915664) T (ValLT (Conc (Def "S" ["s"])) (R 0.1))) (Disj (Until (0.0,3.8955823293172696) T (ValLT (Conc (Def "S" ["s"])) (R 0.1))) (Disj (Until (0.0,3.9959839357429727) T (ValLT (Conc (Def "S" ["s"])) (R 0.1))) (Disj (Until (0.0,4.096385542168675) T (ValLT (Conc (Def "S" ["s"])) (R 0.1))) (Disj (Until (0.0,4.196787148594378) T (ValLT (Conc (Def "S" ["s"])) (R 0.1))) (Disj (Until (0.0,4.2971887550200805) T (ValLT (Conc (Def "S" ["s"])) (R 0.1))) (Disj (Until (0.0,4.397590361445784) T (ValLT (Conc (Def "S" ["s"])) (R 0.1))) (Disj (Until (0.0,4.497991967871486) T (ValLT (Conc (Def "S" ["s"])) (R 0.1))) (Until (0.0,4.598393574297189) T (ValLT (Conc (Def "S" ["s"])) (R 0.1)))))))))))))))

tFR2 = (Disj (Rels (0.0,2.3895582329317264) F (ValLT (Conc (Def "S" ["s"])) (R 0.1))) (Disj (Rels (0.0,2.4899598393574305) F (ValLT (Conc (Def "S" ["s"])) (R 0.1))) (Disj (Rels (0.0,2.5903614457831328) F (ValLT (Conc (Def "S" ["s"])) (R 0.1))) (Disj (Rels (0.0,2.690763052208835) F (ValLT (Conc (Def "S" ["s"])) (R 0.1))) (Disj (Rels (0.0,2.791164658634539) F (ValLT (Conc (Def "S" ["s"])) (R 0.1))) (Disj (Rels (0.0,2.8915662650602414) F (ValLT (Conc (Def "S" ["s"])) (R 0.1))) (Disj (Rels (0.0,2.9919678714859437) F (ValLT (Conc (Def "S" ["s"])) (R 0.1))) (Disj (Rels (0.0,3.092369477911646) F (ValLT (Conc (Def "S" ["s"])) (R 0.1))) (Disj (Rels (0.0,3.19277108433735) F (ValLT (Conc (Def "S" ["s"])) (R 0.1))) (Disj (Rels (0.0,3.2931726907630523) F (ValLT (Conc (Def "S" ["s"])) (R 0.1))) (Rels (0.0,3.3935742971887546) F (ValLT (Conc (Def "S" ["s"])) (R 0.1)))))))))))))

tFR3 = (Conj (Disj (Rels (0.0,24.49799196787149) F (ValGT (Conc (Def "E" ["e"])) (R 1.0e-3))) (Disj (Rels (0.0,24.59839357429719) F (ValGT (Conc (Def "E" ["e"])) (R 1.0e-3))) (Disj (Rels (0.0,24.698795180722893) F (ValGT (Conc (Def "E" ["e"])) (R 1.0e-3))) (Disj (Rels (0.0,24.799196787148595) F (ValGT (Conc (Def "E" ["e"])) (R 1.0e-3))) (Disj (Rels (0.0,24.899598393574298) F (ValGT (Conc (Def "E" ["e"])) (R 1.0e-3))) (Until (0.0,24.49799196787149) T (Rels (0.0,25.0) F (ValGT (Conc (Def "E" ["e"])) (R 1.0e-3))))))))) (Conj (Disj (Rels (0.0,24.59839357429719) F (ValGT (Conc (Def "E" ["e"])) (R 1.0e-3))) (Disj (Rels (0.0,24.698795180722893) F (ValGT (Conc (Def "E" ["e"])) (R 1.0e-3))) (Disj (Rels (0.0,24.799196787148595) F (ValGT (Conc (Def "E" ["e"])) (R 1.0e-3))) (Disj (Rels (0.0,24.899598393574298) F (ValGT (Conc (Def "E" ["e"])) (R 1.0e-3))) (Until (0.0,24.59839357429719) T (Rels (0.0,25.0) F (ValGT (Conc (Def "E" ["e"])) (R 1.0e-3)))))))) (Conj (Disj (Rels (0.0,24.698795180722893) F (ValGT (Conc (Def "E" ["e"])) (R 1.0e-3))) (Disj (Rels (0.0,24.799196787148595) F (ValGT (Conc (Def "E" ["e"])) (R 1.0e-3))) (Disj (Rels (0.0,24.899598393574298) F (ValGT (Conc (Def "E" ["e"])) (R 1.0e-3))) (Until (0.0,24.698795180722893) T (Rels (0.0,25.0) F (ValGT (Conc (Def "E" ["e"])) (R 1.0e-3))))))) (Conj (Disj (Rels (0.0,24.799196787148595) F (ValGT (Conc (Def "E" ["e"])) (R 1.0e-3))) (Disj (Rels (0.0,24.899598393574298) F (ValGT (Conc (Def "E" ["e"])) (R 1.0e-3))) (Until (0.0,24.799196787148595) T (Rels (0.0,25.0) F (ValGT (Conc (Def "E" ["e"])) (R 1.0e-3)))))) (Conj (Disj (Rels (0.0,24.899598393574298) F (ValGT (Conc (Def "E" ["e"])) (R 1.0e-3))) (Until (0.0,24.899598393574298) T (Rels (0.0,25.0) F (ValGT (Conc (Def "E" ["e"])) (R 1.0e-3))))) (Rels (0.0,24.49799196787149) F (Until (0.0,25.0) T (Rels (0.0,25.0) F (ValGT (Conc (Def "E" ["e"])) (R 1.0e-3))))))))))

tFR4 = (Conj (Disj (Rels (0.0,24.899598393574298) F (ValGT (Conc (Def "E" ["e"])) (R 1.0e-3))) (Until (0.0,23.895582329317275) T (Rels (0.0,25.0) F (ValGT (Conc (Def "E" ["e"])) (R 1.0e-3))))) (Conj (Disj (Rels (0.0,24.899598393574298) F (ValGT (Conc (Def "E" ["e"])) (R 1.0e-3))) (Until (0.0,23.995983935742977) T (Rels (0.0,25.0) F (ValGT (Conc (Def "E" ["e"])) (R 1.0e-3))))) (Conj (Disj (Rels (0.0,24.899598393574298) F (ValGT (Conc (Def "E" ["e"])) (R 1.0e-3))) (Until (0.0,24.09638554216868) T (Rels (0.0,25.0) F (ValGT (Conc (Def "E" ["e"])) (R 1.0e-3))))) (Conj (Disj (Rels (0.0,24.899598393574298) F (ValGT (Conc (Def "E" ["e"])) (R 1.0e-3))) (Until (0.0,24.19678714859438) T (Rels (0.0,25.0) F (ValGT (Conc (Def "E" ["e"])) (R 1.0e-3))))) (Conj (Disj (Rels (0.0,24.899598393574298) F (ValGT (Conc (Def "E" ["e"])) (R 1.0e-3))) (Until (0.0,24.297188755020084) T (Rels (0.0,25.0) F (ValGT (Conc (Def "E" ["e"])) (R 1.0e-3))))) (Conj (Disj (Rels (0.0,24.899598393574298) F (ValGT (Conc (Def "E" ["e"])) (R 1.0e-3))) (Until (0.0,24.397590361445786) T (Rels (0.0,25.0) F (ValGT (Conc (Def "E" ["e"])) (R 1.0e-3))))) (Conj (Disj (Rels (0.0,24.899598393574298) F (ValGT (Conc (Def "E" ["e"])) (R 1.0e-3))) (Until (0.0,24.49799196787149) T (Rels (0.0,25.0) F (ValGT (Conc (Def "E" ["e"])) (R 1.0e-3))))) (Conj (Disj (Rels (0.0,24.899598393574298) F (ValGT (Conc (Def "E" ["e"])) (R 1.0e-3))) (Until (0.0,24.59839357429719) T (Rels (0.0,25.0) F (ValGT (Conc (Def "E" ["e"])) (R 1.0e-3))))) (Conj (Disj (Rels (0.0,24.899598393574298) F (ValGT (Conc (Def "E" ["e"])) (R 1.0e-3))) (Until (0.0,24.698795180722893) T (Rels (0.0,25.0) F (ValGT (Conc (Def "E" ["e"])) (R 1.0e-3))))) (Conj (Disj (Rels (0.0,24.899598393574298) F (ValGT (Conc (Def "E" ["e"])) (R 1.0e-3))) (Until (0.0,24.799196787148595) T (Rels (0.0,25.0) F (ValGT (Conc (Def "E" ["e"])) (R 1.0e-3))))) (Conj (Disj (Rels (0.0,24.899598393574298) F (ValGT (Conc (Def "E" ["e"])) (R 1.0e-3))) (Until (0.0,24.899598393574298) T (Rels (0.0,25.0) F (ValGT (Conc (Def "E" ["e"])) (R 1.0e-3))))) (Rels (0.0,23.895582329317275) F (Until (0.0,25.0) T (Rels (0.0,25.0) F (ValGT (Conc (Def "E" ["e"])) (R 1.0e-3))))))))))))))))

tf z@(Disj (Until (0.0,t1) a b) (Disj (Until (0.0,t2) c d) e))
    | a==c&&b==d
        = (Disj (Until (0,(max t1 t2)) a b) e)
    | otherwise
        = z

--------------------------
-- Graph plotting tests:
--------------------------

{-
tPlot = plot ts dims
    where
      ts = [0.0,0.05..60.0] -- ::[Double]
      {-dims = [("Times 2",map (*2) ts),
              ("Times 3",map (*3) ts)
             ]-}
      dims = [("Sin x + " ++ show x, map sin (map (+x) ts)) | x<-[0..5]]

tPlotTimeSeries = do env <- tEnv "models/testEnzyme"
                     let pi = tProc env "Pi"
                     let mts = processMTS env pi
                     let pi' = wholeProc env pi mts
                     let dpdt = dPdt' env mts pi'
                     let odes = xdot env dpdt
                     let inits = initials env pi' dpdt
                     let ts = timePoints 250 (0,25)
                     let soln = solveODE env pi' dpdt (250,(0,25))
                     let ss = speciesIn env dpdt
                     plotTimeSeries ts soln ss
-}

---------------------------------
-- Testing solver with Jacobain:
---------------------------------

{-tSolveMAPKwithJac = do env <- tEnv "models/mapk.cpi"
                       let pi = tProc env "MAPK"
                           mts = processMTS env pi
                           pi' = wholeProc env pi mts
                           dpdt = dPdt' env mts pi'
                           odes = xdot env dpdt
                           jacob = jac env dpdt
                           inits = initials env pi' dpdt
                           ts = timePoints 800 (0,80)
                           soln = solveODE' odes jacob inits ts
                       return soln

tSolveMAPKwithoutJac = do env <- tEnv "models/mapk.cpi"
                          let pi = tProc env "MAPK"
                              mts = processMTS env pi
                              pi' = wholeProc env pi mts
                              dpdt = dPdt' env mts pi'
                              odes = xdot env dpdt
                              inits = initials env pi' dpdt
                              ts = timePoints 800 (0,80)
                              soln = solveODE odes inits ts
                          return soln

tPlotEnzwithJac = do env <- tEnv "models/testEnzyme.cpi"
                     let pi = tProc env "Pi"
                         mts = processMTS env pi
                         pi' = wholeProc env pi mts
                         dpdt = dPdt' env mts pi'
                         odes = xdot env dpdt
                         jacob = jac env dpdt
                         inits = initials env pi' dpdt
                         ts = timePoints 250 (0,25)
                         soln = solveODE' odes jacob inits ts
                         ss = speciesIn env dpdt
                     plotTimeSeries ts soln ss

tPlotEnzwithoutJac = do env <- tEnv "models/testEnzyme.cpi"
                        let pi = tProc env "Pi"
                            mts = processMTS env pi
                            pi' = wholeProc env pi mts
                            dpdt = dPdt' env mts pi'
                            odes = xdot env dpdt
                            inits = initials env pi' dpdt
                            ts = timePoints 250 (0,25)
                            soln = solveODE env pi' dpdt (250,(0,25)) ts
                            ss = speciesIn env dpdt
                        plotTimeSeries ts soln ss
-}

max_stars_repo_name: twright/bondwb

max_issues_repo_name: twright/bondwb

max_forks_repo_name: twright/bondwb