text:
{-# LANGUAGE CPP                  #-}
{-# LANGUAGE DeriveDataTypeable   #-}
{-# LANGUAGE FlexibleContexts     #-}
{-# LANGUAGE NoImplicitPrelude    #-}
{-# LANGUAGE ScopedTypeVariables  #-}
{-# LANGUAGE StandaloneDeriving   #-}
{-# LANGUAGE UndecidableInstances #-}

--------------------------------------------------------------------------------
-- |
-- Module:     Numeric.LevMar
-- Copyright:  (c) 2009 - 2014 Roel van Dijk & Bas van Dijk
-- License:    BSD-style (see the file LICENSE)
-- Maintainer: Roel van Dijk <vandijk.roel@gmail.com>
--             Bas van Dijk <v.dijk.bas@gmail.com>
-- Stability:  Experimental
--
-- For additional documentation see the
-- <http://www.ics.forth.gr/~lourakis/levmar/ documentation of the levmar C>
-- library which this library is based on:
--
--------------------------------------------------------------------------------

module Numeric.LevMar
    ( -- * Model & Jacobian.
      Params, Samples
    , Model, Jacobian

      -- * Levenberg-Marquardt algorithm.
    , LevMarable(levmar)

      -- * Minimization options.
    , Options(..), defaultOpts

      -- * Constraints
    , Constraints(..), LinearConstraints

      -- * Output
    , Info(..), StopReason(..), LevMarError(..)
    ) where

-------------------------------------------------------------------------------
-- Imports
-------------------------------------------------------------------------------

-- from base:
import Control.Monad         ( return, mplus )
import Control.Exception     ( Exception )
import Data.Bool             ( (&&), (||), otherwise )
import Data.Data             ( Data )
import Data.Typeable         ( Typeable )
import Data.Either           ( Either(Left, Right) )
import Data.Eq               ( Eq, (==), (/=) )
import Data.Function         ( (.), ($) )
import Data.Functor          ( (<$>) )
import Data.Int              ( Int )
import Data.List             ( lookup, (++) )
import Data.Maybe            ( Maybe(Nothing, Just), isJust, fromJust, fromMaybe )
import Data.Monoid           ( Monoid, mempty, mappend )
import Data.Ord              ( Ord, (<) )
import Foreign.C.Types       ( CInt )
import Foreign.Marshal.Array ( allocaArray, withArray, peekArray, copyArray )
import Foreign.Ptr           ( Ptr, nullPtr )
import Foreign.ForeignPtr    ( ForeignPtr, newForeignPtr_, withForeignPtr )
import Foreign.Storable      ( Storable )
import Prelude               ( Num, Enum, Fractional, RealFrac, Float, Double
                             , fromIntegral, toEnum, (-), (*), error, floor
                             )
import System.IO             ( IO )
import System.IO.Unsafe      ( unsafePerformIO )
import Text.Read             ( Read )
import Text.Show             ( Show, show )

#if __GLASGOW_HASKELL__ >= 605
import GHC.ForeignPtr        ( mallocPlainForeignPtrBytes )
import Prelude               ( undefined )
import Foreign.Storable      ( sizeOf )
#else
import Foreign.ForeignPtr    ( mallocForeignPtrArray )
#endif

#if __GLASGOW_HASKELL__ < 700
import Prelude               ( fromInteger, (>>=), (>>), fail )
#endif

-- from hmatrix:
#if MIN_VERSION_hmatrix(0,17,0)
import Numeric.LinearAlgebra.Data ( Matrix, flatten, rows, reshape )
import Numeric.LinearAlgebra      ( Container, Element )
#else
import Data.Packed.Matrix         ( Matrix, Element, flatten, rows, reshape )
import Numeric.Container          ( Container )
import Numeric.LinearAlgebra      ( {- Instances for Matrix -} )
#endif

-- from vector:
import           Data.Vector.Storable       ( Vector )
import qualified Data.Vector.Storable as VS ( unsafeWith, length
                                            , unsafeFromForeignPtr
                                            , length
                                            )

-- from bindings-levmar:
import Bindings.LevMar ( c'LM_INFO_SZ

                       , withModel
                       , withJacobian

                       , c'LM_ERROR
                       , c'LM_ERROR_LAPACK_ERROR
                       , c'LM_ERROR_FAILED_BOX_CHECK
                       , c'LM_ERROR_MEMORY_ALLOCATION_FAILURE
                       , c'LM_ERROR_CONSTRAINT_MATRIX_ROWS_GT_COLS
                       , c'LM_ERROR_CONSTRAINT_MATRIX_NOT_FULL_ROW_RANK
                       , c'LM_ERROR_TOO_FEW_MEASUREMENTS
                       , c'LM_ERROR_SINGULAR_MATRIX
                       , c'LM_ERROR_SUM_OF_SQUARES_NOT_FINITE

                       , c'LM_INIT_MU
                       , c'LM_STOP_THRESH
                       , c'LM_DIFF_DELTA
                       )
import qualified Bindings.LevMar ( Model, Jacobian )

-- from levmar (this package):
import Bindings.LevMar.CurryFriendly ( LevMarDer,     LevMarDif
                                     , LevMarBCDer,   LevMarBCDif
                                     , LevMarLecDer,  LevMarLecDif
                                     , LevMarBLecDer, LevMarBLecDif

                                     , dlevmar_der,      slevmar_der
                                     , dlevmar_dif,      slevmar_dif
                                     , dlevmar_bc_der,   slevmar_bc_der
                                     , dlevmar_bc_dif,   slevmar_bc_dif
                                     , dlevmar_lec_der,  slevmar_lec_der
                                     , dlevmar_lec_dif,  slevmar_lec_dif
                                     , dlevmar_blec_der, slevmar_blec_der
                                     , dlevmar_blec_dif, slevmar_blec_dif
                                     )

--------------------------------------------------------------------------------
-- Model & Jacobian.
--------------------------------------------------------------------------------

-- | Parameter vector of length @m@.
--
-- Ensure that @m <= n@ where @n@ is the length of the 'Samples' vector.
type Params r = Vector r

-- | Sample vector of length @n@.
--
-- Ensure that @n >= m@ where @m@ is the length of the 'Params' vector.
type Samples r = Vector r

{-| A functional relation describing measurements represented as a function
from a vector of parameters to a vector of expected samples.

 * Ensure that the length @m@ of the parameter vector equals the length of the
   initial parameter vector in 'levmar'.

 * Ensure that the length @n@ of the output sample vector equals the length of
   the sample vector in 'levmar'.

 * Ensure that the length @n@ of the output sample vector vector is bigger than or
   equal to the length @m@ of the parameter vector.
-}
type Model r = Params r -> Samples r

{-| The <http://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant jacobian>
of the 'Model' function. Expressed as a function from a vector of
parameters to a matrix which for each expected sample describes the
partial derivatives of the parameters.

 * Ensure that the length @m@ of the parameter vector equals the length of the
   initial parameter vector in 'levmar'.

 * Ensure that the output matrix has the dimension @n><m@ where @n@ is the
   number of samples and @m@ is the number of parameters.
-}
type Jacobian r = Params r -> Matrix r

--------------------------------------------------------------------------------
-- Levenberg-Marquardt algorithm.
--------------------------------------------------------------------------------

-- | The Levenberg-Marquardt algorithm is overloaded to work on 'Double' and 'Float'.
class LevMarable r where

    -- | The Levenberg-Marquardt algorithm.
    --
    -- Returns a triple of the found parameters, a structure containing
    -- information about the minimization and the covariance matrix
    -- corresponding to LS solution.
    --
    -- Ensure that @n >= m@.
    levmar :: Model r            -- ^ Model
           -> Maybe (Jacobian r) -- ^ Optional jacobian
           -> Params r           -- ^ Initial parameters of length @m@
           -> Samples r          -- ^ Sample vector of length @n@
           -> Int                -- ^ Maximum iterations
           -> Options r          -- ^ Minimization options
           -> Constraints r      -- ^ Constraints
           -> Either LevMarError (Params r, Info r, Matrix r)

instance LevMarable Float where
    levmar = gen_levmar slevmar_der
                        slevmar_dif
                        slevmar_bc_der
                        slevmar_bc_dif
                        slevmar_lec_der
                        slevmar_lec_dif
                        slevmar_blec_der
                        slevmar_blec_dif

instance LevMarable Double where
    levmar = gen_levmar dlevmar_der
                        dlevmar_dif
                        dlevmar_bc_der
                        dlevmar_bc_dif
                        dlevmar_lec_der
                        dlevmar_lec_dif
                        dlevmar_blec_der
                        dlevmar_blec_dif

{-| @gen_levmar@ takes the low-level C functions as arguments and
executes one of them depending on the optional jacobian and constraints.

Preconditions:

@
  length ys >= length ps

     isJust mLowBs && length (fromJust mLowBs) == length ps
  && isJust mUpBs  && length (fromJust mUpBs)  == length ps

  boxConstrained && (all $ zipWith (<=) (fromJust mLowBs) (fromJust mUpBs))
@
-}
gen_levmar :: forall r. (RealFrac r, Element r)
           => LevMarDer r
           -> LevMarDif r
           -> LevMarBCDer r
           -> LevMarBCDif r
           -> LevMarLecDer r
           -> LevMarLecDif r
           -> LevMarBLecDer r
           -> LevMarBLecDif r

           -> Model r            -- ^ Model
           -> Maybe (Jacobian r) -- ^ Optional jacobian
           -> Params r           -- ^ Initial parameters
           -> Samples r          -- ^ Samples
           -> Int                -- ^ Maximum iterations
           -> Options r          -- ^ Options
           -> Constraints r      -- ^ Constraints
           -> Either LevMarError (Params r, Info r, Matrix r)
gen_levmar f_der
           f_dif
           f_bc_der
           f_bc_dif
           f_lec_der
           f_lec_dif
           f_blec_der
           f_blec_dif
           model mJac ps ys itMax opts (Constraints mLowBs mUpBs mWeights mLinC)
               | m == 0    = Left LevMarError -- LAPACK will crash otherwise!
               | otherwise =
  -- All effects are contained, so we can safely perform:
  unsafePerformIO $ do

    -- We need to pass the initial parameters 'ps' to the C function.
    -- However, we can't just pass a pointer to them because the C function
    -- will modify the parameters during execution which will violate
    -- referential transparanency. Instead we allocate new space
    -- and copy the parameters to it.
    --
    -- Note that, in the end, the array is returned from this function.
    -- This means that the only way to guarantee its finalisation
    -- is to allocate it using a ForeignPtr:
    psFP <- fastMallocForeignPtrArray m
    withForeignPtr psFP $ \psPtr -> do
      VS.unsafeWith ps $ \psPtrInp ->
        copyArray psPtr psPtrInp m

      -- Retrieve the (read-only) pointer 'ysPtr' to the samples vector 'ys'
      -- so we can pass it to the C function:
      VS.unsafeWith ys $ \ysPtr ->

        -- Convert the Options 'opts' to a list and then to an array
        -- so we can pass the (read-only) pointer 'optsPtr' to the C function:
        withArray (optsToList opts) $ \optsPtr ->

          -- Allocate space for the info array
          -- so we can pass it to the C function.
          -- Note that, in the end, this array is converted to an Info value
          -- and returned from this function.
          allocaArray c'LM_INFO_SZ $ \infoPtr -> do

            -- Allocate space for the covariance matrix
            -- so we can pass it to the C function.
            -- Like the parameters array the matrix
            -- needs to be returned from this function.
            -- So we also allocate it using a ForeignPtr:
            covarFP <- fastMallocForeignPtrArray mm
            withForeignPtr covarFP $ \covarPtr ->

              -- 'cmodel' is the low-level model function which is converted
              -- to the FunPtr 'modelFunPtr' and passed to the C function.
              -- 'cmodel' will first convert the parameters pointer 'parPtr'
              -- into a Vector after converting it into a ForeignPtr
              -- (without a finalizer).
              -- Then it will apply the high-level 'model' function
              -- to this parameter vector. The resulting vector is then copied
              -- to the output buffer 'hxPtr':
              let cmodel :: Bindings.LevMar.Model r
                  cmodel parPtr hxPtr _ _ _ = do
                    parFP <- newForeignPtr_ parPtr
                    let psV = VS.unsafeFromForeignPtr parFP 0 m
                        vector = model psV
                    VS.unsafeWith vector $ \p -> copyArray hxPtr p (VS.length vector)
              in withModel cmodel $ \modelFunPtr -> do

                 -- All the low-level C functions share a common set of arguments.
                 -- 'runDif' applies these arguments to the given C function 'f':
                 let runDif :: LevMarDif r -> IO CInt
                     runDif f = f modelFunPtr
                                  psPtr
                                  ysPtr
                                  (fromIntegral m)
                                  (fromIntegral n)
                                  (fromIntegral itMax)
                                  optsPtr
                                  infoPtr
                                  nullPtr
                                  covarPtr
                                  nullPtr

                 err <- case mJac of
                   Nothing -> if boxConstrained
                              then if linConstrained
                                   then withBoxConstraints
                                          (withLinConstraints $ withWeights runDif)
                                          f_blec_dif
                                   else withBoxConstraints runDif f_bc_dif
                              else if linConstrained
                                   then withLinConstraints runDif f_lec_dif
                                   else runDif f_dif

                   Just jac ->
                     let cjacobian :: Bindings.LevMar.Jacobian r
                         cjacobian parPtr jPtr _ _ _ = do
                           parFP <- newForeignPtr_ parPtr
                           let psV    = VS.unsafeFromForeignPtr parFP 0 m
                               matrix = jac psV
                               vector = flatten matrix
                           VS.unsafeWith vector $ \p ->
                             copyArray jPtr p (VS.length vector)
                     in withJacobian cjacobian $ \jacobPtr ->

                       let runDer :: LevMarDer r -> IO CInt
                           runDer f = runDif $ f jacobPtr
                       in if boxConstrained
                          then if linConstrained
                               then withBoxConstraints
                                      (withLinConstraints $ withWeights runDer)
                                      f_blec_der
                               else withBoxConstraints runDer f_bc_der
                          else if linConstrained
                               then withLinConstraints runDer f_lec_der
                               else runDer f_der

                 -- Handling errors:
                 if err < 0
                    -- we don't treat the following two as an error:
                    && err /= c'LM_ERROR_SINGULAR_MATRIX
                    && err /= c'LM_ERROR_SUM_OF_SQUARES_NOT_FINITE
                   then return $ Left $ convertLevMarError err

                   else do -- Converting results:
                     info <- listToInfo <$> peekArray c'LM_INFO_SZ infoPtr
                     let psV = VS.unsafeFromForeignPtr psFP 0 m
                     let covarM = reshape m $ VS.unsafeFromForeignPtr covarFP 0 mm

                     return $ Right (psV, info, covarM)
  where
    m  = VS.length ps
    n  = VS.length ys
    mm = m*m

    -- Whether the parameters are constrained by a linear equation.
    linConstrained = isJust   mLinC
    (cMat, rhcVec) = fromJust mLinC

    -- Whether the parameters are constrained by a bounding box.
    boxConstrained = isJust mLowBs || isJust mUpBs

    withBoxConstraints f g =
        maybeWithArray mLowBs $ \lBsPtr ->
          maybeWithArray mUpBs $ \uBsPtr ->
            f $ g lBsPtr uBsPtr

    withLinConstraints f g =
        VS.unsafeWith (flatten cMat) $ \cMatPtr ->
          VS.unsafeWith rhcVec $ \rhcVecPtr ->
            f . g cMatPtr rhcVecPtr $ fromIntegral $ rows cMat

    withWeights f g = maybeWithArray mWeights $ f . g

maybeWithArray :: (Storable a) => Maybe (Vector a) -> (Ptr a -> IO β) -> IO β
maybeWithArray Nothing  f = f nullPtr
maybeWithArray (Just v) f = VS.unsafeWith v f

#if __GLASGOW_HASKELL__ >= 605
{-# INLINE fastMallocForeignPtrArray #-}
fastMallocForeignPtrArray :: forall a. Storable a => Int -> IO (ForeignPtr a)
fastMallocForeignPtrArray n = mallocPlainForeignPtrBytes
                                (n * sizeOf (undefined :: a))
#else
fastMallocForeignPtrArray :: forall a. Storable a => Int -> IO (ForeignPtr a)
fastMallocForeignPtrArray = mallocForeignPtrArray
#endif

--------------------------------------------------------------------------------
-- Minimization options.
--------------------------------------------------------------------------------

-- | Minimization options
data Options r =
    Opts { optScaleInitMu      :: !r -- ^ Scale factor for initial /mu/.
         , optStopNormInfJacTe :: !r -- ^ Stopping thresholds for @||J^T e||_inf@.
         , optStopNorm2Dp      :: !r -- ^ Stopping thresholds for @||Dp||_2@.
         , optStopNorm2E       :: !r -- ^ Stopping thresholds for @||e||_2@.
         , optDelta            :: !r -- ^ Step used in the difference
                                     -- approximation to the Jacobian. If
                                     -- @optDelta<0@, the Jacobian is approximated
                                     -- with central differences which are more
                                     -- accurate (but slower!)  compared to the
                                     -- forward differences employed by default.
         } deriving (Eq, Ord, Read, Show, Data, Typeable)

-- | Default minimization options
defaultOpts :: Fractional r => Options r
defaultOpts = Opts { optScaleInitMu      = c'LM_INIT_MU
                   , optStopNormInfJacTe = c'LM_STOP_THRESH
                   , optStopNorm2Dp      = c'LM_STOP_THRESH
                   , optStopNorm2E       = c'LM_STOP_THRESH
                   , optDelta            = c'LM_DIFF_DELTA
                   }

optsToList :: Options r -> [r]
optsToList (Opts mu  eps1  eps2  eps3  delta) =
                [mu, eps1, eps2, eps3, delta]

--------------------------------------------------------------------------------
-- Constraints
--------------------------------------------------------------------------------

-- | Ensure that these vectors have the same length as the number of parameters.
data Constraints r = Constraints
    { lowerBounds       :: !(Maybe (Params r))            -- ^ Optional lower bounds
    , upperBounds       :: !(Maybe (Params r))            -- ^ Optional upper bounds
    , weights           :: !(Maybe (Params r))            -- ^ Optional weights
    , linearConstraints :: !(Maybe (LinearConstraints r)) -- ^ Optional linear constraints
    } deriving (Read, Show, Typeable)

deriving instance (Eq r, Container Vector r, Num r) => Eq (Constraints r)

-- | Linear constraints consisting of a constraints matrix, @k><m@ and
--   a right hand constraints vector, of length @k@ where @m@ is the number of
--   parameters and @k@ is the number of constraints.
type LinearConstraints r = (Matrix r, Vector r)

-- | * 'mempty' is defined as a 'Constraints' where all fields are 'Nothing'.
--
--   * 'mappend' merges two 'Constraints' by taking the first non-'Nothing' value
--     for each field.
instance Monoid (Constraints r) where
    mempty = Constraints Nothing Nothing Nothing Nothing
    mappend (Constraints lb1 ub1 w1 l1)
            (Constraints lb2 ub2 w2 l2) = Constraints (lb1 `mplus` lb2)
                                                      (ub1 `mplus` ub2)
                                                      (w1  `mplus` w2)
                                                      (l1  `mplus` l2)

--------------------------------------------------------------------------------
-- Output
--------------------------------------------------------------------------------

-- | Information regarding the minimization.
data Info r = Info
  { infNorm2initE      :: !r          -- ^ @||e||_2@             at initial parameters.
  , infNorm2E          :: !r          -- ^ @||e||_2@             at estimated parameters.
  , infNormInfJacTe    :: !r          -- ^ @||J^T e||_inf@       at estimated parameters.
  , infNorm2Dp         :: !r          -- ^ @||Dp||_2@            at estimated parameters.
  , infMuDivMax        :: !r          -- ^ @\mu/max[J^T J]_ii ]@ at estimated parameters.
  , infNumIter         :: !Int        -- ^ Number of iterations.
  , infStopReason      :: !StopReason -- ^ Reason for terminating.
  , infNumFuncEvals    :: !Int        -- ^ Number of function evaluations.
  , infNumJacobEvals   :: !Int        -- ^ Number of jacobian evaluations.
  , infNumLinSysSolved :: !Int        -- ^ Number of linear systems solved,
                                      --   i.e. attempts for reducing error.
  } deriving (Eq, Ord, Read, Show, Data, Typeable)

listToInfo :: (RealFrac r) => [r] -> Info r
listToInfo [a,b,c,d,e,f,g,h,i,j] =
    Info { infNorm2initE      = a
         , infNorm2E          = b
         , infNormInfJacTe    = c
         , infNorm2Dp         = d
         , infMuDivMax        = e
         , infNumIter         = floor f
         , infStopReason      = toEnum $ floor g - 1
         , infNumFuncEvals    = floor h
         , infNumJacobEvals   = floor i
         , infNumLinSysSolved = floor j
         }
listToInfo _ = error "liftToInfo: wrong list length"

-- | Reason for terminating.
data StopReason
  = SmallGradient  -- ^ Stopped because of small gradient @J^T e@.
  | SmallDp        -- ^ Stopped because of small Dp.
  | MaxIterations  -- ^ Stopped because maximum iterations was reached.
  | SingularMatrix -- ^ Stopped because of singular matrix. Restart from current
                   --   estimated parameters with increased 'optScaleInitMu'.
  | SmallestError  -- ^ Stopped because no further error reduction is
                   --   possible. Restart with increased 'optScaleInitMu'.
  | SmallNorm2E    -- ^ Stopped because of small @||e||_2@.
  | InvalidValues  -- ^ Stopped because model function returned invalid values
                   --   (i.e. NaN or Inf). This is a user error.
    deriving (Eq, Ord, Read, Show, Data, Typeable, Enum)

--------------------------------------------------------------------------------
-- Error
--------------------------------------------------------------------------------

data LevMarError
    = LevMarError                    -- ^ Generic error (not one of the others)
    | LapackError                    -- ^ A call to a lapack subroutine failed
                                     --   in the underlying C levmar library.
    | FailedBoxCheck                 -- ^ At least one lower bound exceeds the
                                     --   upper one.
    | MemoryAllocationFailure        -- ^ A call to @malloc@ failed in the
                                     --   underlying C levmar library.
    | ConstraintMatrixRowsGtCols     -- ^ The matrix of constraints cannot have
                                     --   more rows than columns.
    | ConstraintMatrixNotFullRowRank -- ^ Constraints matrix is not of full row
                                     --   rank.
    | TooFewMeasurements             -- ^ Cannot solve a problem with fewer
                                     --   measurements than unknowns.  In case
                                     --   linear constraints are provided, this
                                     --   error is also returned when the number
                                     --   of measurements is smaller than the
                                     --   number of unknowns minus the number of
                                     --   equality constraints.
      deriving (Eq, Ord, Read, Show, Data, Typeable)

-- Handy in case you want to thow a LevMarError as an exception:
instance Exception LevMarError

levmarCErrorToLevMarError :: [(CInt, LevMarError)]
levmarCErrorToLevMarError =
    [ (c'LM_ERROR,                                     LevMarError)
    , (c'LM_ERROR_LAPACK_ERROR,                        LapackError)
  --, (c'LM_ERROR_NO_JACOBIAN,                         can never happen)
  --, (c'LM_ERROR_NO_BOX_CONSTRAINTS,                  can never happen)
    , (c'LM_ERROR_FAILED_BOX_CHECK,                    FailedBoxCheck)
    , (c'LM_ERROR_MEMORY_ALLOCATION_FAILURE,           MemoryAllocationFailure)
    , (c'LM_ERROR_CONSTRAINT_MATRIX_ROWS_GT_COLS,      ConstraintMatrixRowsGtCols)
    , (c'LM_ERROR_CONSTRAINT_MATRIX_NOT_FULL_ROW_RANK, ConstraintMatrixNotFullRowRank)
    , (c'LM_ERROR_TOO_FEW_MEASUREMENTS,                TooFewMeasurements)
  --, (c'LM_ERROR_SINGULAR_MATRIX,                     we don't treat this as an error)
  --, (c'LM_ERROR_SUM_OF_SQUARES_NOT_FINITE,           we don't treat this as an error)
    ]

convertLevMarError :: CInt -> LevMarError
convertLevMarError err = fromMaybe (error $ "Unknown levmar error: " ++ show err)
                                   (lookup err levmarCErrorToLevMarError)

max_stars_repo_name: ann-yang-leapyear/levmar

max_issues_repo_name: ann-yang-leapyear/levmar

max_forks_repo_name: basvandijk/levmar


text:
module STCR2Z2T0S0EdgeBinary where

import           Control.Arrow
import           Control.Monad
import           Data.Array.Repa         as R
import           Data.Binary             (encodeFile, decodeFile)
import           Data.Complex
import           Data.List               as L
import           DFT.Plan
import           Filter.Utils
import           FokkerPlanck.MonteCarlo
import           FokkerPlanck.Pinwheel
import           Image.Edge
import           Image.IO
import           STC
import           System.Directory
import           System.Environment
import           System.FilePath
import           Text.Printf
import           Types
import           Utils.Array
import           Utils.Parallel

main = do
  args@(numPointStr:numOrientationStr:numScaleStr:thetaSigmaStr:scaleSigmaStr:maxScaleStr:taoStr:numTrailStr:maxTrailStr:thetaFreqsStr:scaleFreqsStr:hollowRadiusStr:cutoffRadiusStr:histFilePath:filterFileFolder:numIterationStr:writeSourceFlagStr:saveEdgeDataFlagStr:loadEdgeDataFlagStr:edgeFilePath:numNoisePointStr:scaleFactorStr:useFFTWWisdomFlagStr:fftwWisdomFileName:batchSizeStr:numThreadStr:_) <-
    getArgs
  print args
  let numPoint = read numPointStr :: Int
      numOrientation = read numOrientationStr :: Int
      numScale = read numScaleStr :: Int
      thetaSigma = read thetaSigmaStr :: Double
      scaleSigma = read scaleSigmaStr :: Double
      maxScale = read maxScaleStr :: Double
      tao = read taoStr :: Double
      numTrail = read numTrailStr :: Int
      maxTrail = read maxTrailStr :: Int
      thetaFreq = read thetaFreqsStr :: Double
      thetaFreqs = [-thetaFreq .. thetaFreq]
      scaleFreq = read scaleFreqsStr :: Double
      scaleFreqs = [-scaleFreq .. scaleFreq]
      hollowRadius = read hollowRadiusStr :: Double
      cutoffRadius = read cutoffRadiusStr :: Double
      numIteration = read numIterationStr :: Int
      writeSourceFlag = read writeSourceFlagStr :: Bool
      saveEdgeDataFlag = read saveEdgeDataFlagStr :: Bool
      loadEdgeDataFlag = read loadEdgeDataFlagStr :: Bool
      numNoisePoint = read numNoisePointStr :: Int
      scaleFactor = read scaleFactorStr :: Double
      useFFTWWisdomFlag = read useFFTWWisdomFlagStr :: Bool
      batchSize = read batchSizeStr :: Int
      numThread = read numThreadStr :: Int
      parallelParams = ParallelParams numThread batchSize
      folderPath = "output/test/STCR2Z2T0S0EdgeBinary"
      fftwWisdomFilePath = folderPath </> fftwWisdomFileName
      filterFileName =
        printf
          "Filter_%.0f_%.0f_%s"
          hollowRadius
          cutoffRadius
          (takeFileName histFilePath)
      filterFilePath = filterFileFolder </> filterFileName
      biasFilePath =
        printf
          "%s/Bias_%s_%d.dat"
          filterFileFolder
          (takeBaseName edgeFilePath)
          numNoisePoint
      eigenVecFilePath =
        printf
          "%s/EigenVec_%s_%d.dat"
          filterFileFolder
          (takeBaseName edgeFilePath)
          numNoisePoint
  createDirectoryIfMissing True folderPath
  createDirectoryIfMissing True filterFileFolder
  plan <-
    makePlanBinary
      emptyPlan
      useFFTWWisdomFlag
      fftwWisdomFilePath
      (L.length thetaFreqs)
      (L.length scaleFreqs)
      numPoint
      numPoint
  filterFlag <- doesFileExist filterFilePath
  flag <-
    if filterFlag
      then do
        size <- getFileSize filterFilePath
        return $
          if size == 0
            then False
            else True
      else return False
  unless
    flag
    (do histFlag <- doesFileExist histFilePath
        radialArr <-
          if histFlag
            then R.map magnitude . getNormalizedHistogramArr <$>
                 decodeFile histFilePath
            else do
              putStrLn
                "Couldn't find a Green's function data. Start simulation..."
              solveMonteCarloR2Z2T0S0Radial
                numThread
                numTrail
                maxTrail
                numPoint
                numPoint
                thetaSigma
                scaleSigma
                maxScale
                tao
                thetaFreqs
                thetaFreqs
                scaleFreqs
                scaleFreqs
                histFilePath
                (emptyHistogram
                   [ (round . sqrt . fromIntegral $ 2 * (div numPoint 2) ^ 2)
                   , L.length scaleFreqs
                   , L.length thetaFreqs
                   , L.length scaleFreqs
                   , L.length thetaFreqs
                   ]
                   0)
        putStrLn "Write filter to disk..."
        writeDFTPinwheel
          parallelParams
          plan
          radialArr
          hollowRadius
          cutoffRadius
          thetaFreqs
          scaleFreqs
          maxScale
          (numPoint, numPoint)
          filterFilePath)
  (bias, eigenVec) <-
    if loadEdgeDataFlag
      then do
        bias <- readRepaArray biasFilePath
        eigenVec <- readRepaArray eigenVecFilePath
        return (bias, eigenVec)
      else do
        xs <- parseEdgeFile edgeFilePath
        randomPonintSet <-
          generateRandomPointSet numNoisePoint numPoint numPoint
        let (centerX, centerY) =
              join (***) (\x -> div x . L.length $ xs) .
              L.foldl'
                (\(a, b) (R2S1RPPoint (c, d, _, _)) -> (a + c, b + d))
                (0, 0) $
              xs
            ys =
              L.map
                (\(R2S1RPPoint (a, b, c, d)) ->
                   (R2S1RPPoint
                      ( round $ (fromIntegral $ a - centerX) / scaleFactor
                      , round $ (fromIntegral $ b - centerX) / scaleFactor
                      , c
                      , d)))
                xs
            zs =
              L.map
                (\(R2S1RPPoint (a, b, c, d)) ->
                   (R2S1RPPoint (a - center numPoint, b - center numPoint, c, d)))
                randomPonintSet
            points = ys L.++ zs
        let bias =
              computeBiasR2T0S0 numPoint numPoint thetaFreqs scaleFreqs points
            eigenVec =
              computeInitialEigenVectorBinary
                numPoint
                numPoint
                thetaFreqs
                scaleFreqs
                points
        when saveEdgeDataFlag (writeRepaArray biasFilePath bias)
        return (bias, eigenVec)
  powerMethodBinary
    parallelParams
    plan
    folderPath
    numPoint
    numPoint
    numOrientation
    thetaFreqs
    numScale
    scaleFreqs
    maxScale
    filterFilePath
    numIteration
    writeSourceFlag
    (printf
       "_%.2f_%.2f_%d_%d_%d_%d"
       thetaSigma
       scaleSigma
       (round maxScale :: Int)
       (round tao :: Int)
       (round thetaFreq :: Int)
       (round scaleFreq :: Int))
    saveEdgeDataFlag
    eigenVecFilePath
    bias
    eigenVec

ext: hs

max_stars_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_issues_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_forks_repo_name: XinhuaZhang/Stochastic-Completion-Field


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

import MachineLearning.Utils.Config
import MachineLearning.Utils.Report
import MachineLearning.Utils.Data
import MachineLearning.TIR
import MachineLearning.TIR.Crossover
import MachineLearning.TIR.Mutation
import MachineLearning.Model.Fitness
import MachineLearning.Model.Measure
import Control.Evolution
import Data.SRTree
import Data.SRTree.Print

import System.Environment (getArgs)
import Control.Monad.State.Strict
import Data.Bifunctor             (second)
import Data.List.Split            (splitOn)
import Data.Maybe                 (fromJust, fromMaybe)
import Numeric.LinearAlgebra      (size)
import System.Random              (mkStdGen, getStdGen)
import System.IO                  (Handle)
import System.Clock               (Clock(..), getTime, sec)
import Numeric.ModalInterval
import qualified Data.Vector.Storable  as VS
import qualified Data.Vector  as V
import qualified Numeric.LinearAlgebra as LA

import Data.List (intercalate)
import Prelude hiding (null)
import Algorithm.ShapeConstraint (getViolationFun)

filterImage :: Kaucher Double -> [Function] -> [Function]
filterImage image = filter (isNotNullNaNInf . (`evalFun` image) . inverseFunc)
{-# INLINE filterImage #-}

isNotNullNaNInf :: Kaucher Double -> Bool
isNotNullNaNInf xs = not (isEmpty xs || isNaN xl || isNaN xu || isInfinite xl || isInfinite xu)
  where
    xl = fromMaybe (-1/0) $ inf xs
    xu = fromMaybe (1/0) $ sup xs
{-# INLINE isNotNullNaNInf #-}

instance EvoClass Individual where
  data Crossover Individual = OnePoint | UniformCX deriving (Show, Read)
  data Mutation  Individual = GroupMutation deriving (Show, Read)

parseCLI :: Task -> [String] -> IO ()
parseCLI Regression [expminP, expmaxP, tfuncsP, ytfuncsP, errorMetric, nGensP, nPopP, pcP, pmP, seedP, penalty, trainname] = do
  let mutCfg = dfltMutCfg { _kRange = (read expminP, read expmaxP)
                          , _yfuns  = map read $ splitOn "," ytfuncsP
                          , _funs   = map read $ splitOn "," tfuncsP
                          , _vars   = []
                          }
      algCfg = dfltAlgCfg { _gens     = read nGensP
                          , _nPop     = read nPopP
                          , _pm       = read pmP
                          , _pc       = read pcP
                          , _seed     = let s = read seedP in if s < 0 then Nothing else Just s
                          , _measures = [toMeasure errorMetric]
                          }
      ioCfg = IOCfg trainname trainname Screen
      pn    = read penalty
      cnst  = if pn == 0.0 then dfltCnstrCfg  else dfltCnstrCfg{ _penaltyType  = Len pn }
      cfg   = Conf mutCfg ioCfg algCfg cnst
  (champion, _, _) <- runGP cfg
  let bias = V.head $ VS.convert $ head $ _weights champion
  print $ (showPython . assembleTree bias . _chromo) champion <> ";" <> (show . _len) champion <> ";" <> (show . _fit) champion

parseCLI task [expminP, expmaxP, tfuncsP, ytfuncsP, errorMetric, nGensP, nPopP, pcP, pmP, seedP, penalty, niter, trainname] = do
  let mutCfg = dfltMutCfg { _kRange = (read expminP, read expmaxP)
                          , _yfuns  = map read $ splitOn "," ytfuncsP
                          , _funs   = map read $ splitOn "," tfuncsP
                          , _vars   = []
                          }
      algCfg = dfltAlgCfg { _gens     = read nGensP
                          , _nPop     = read nPopP
                          , _pm       = read pmP
                          , _pc       = read pcP
                          , _seed     = let s = read seedP in if s < 0 then Nothing else Just s
                          , _measures = [toMeasure errorMetric]
                          , _task     = case task of
                                          RegressionNL _ -> RegressionNL $ read niter
                                          Classification _ -> Classification $ read niter
                                          ClassMult _ -> ClassMult $ read niter
                          }
      ioCfg = IOCfg trainname trainname Screen
      pn    = read penalty
      cnst  = if pn == 0.0 then dfltCnstrCfg  else dfltCnstrCfg{ _penaltyType  = Len pn }
      cfg   = Conf mutCfg ioCfg algCfg cnst
  (champion, _, _) <- runGP cfg

  let
    tir    = _chromo champion
    trees  = intercalate "#" . map (showPython . getTree tir) $ _weights champion
  print $ trees <> ";" <> (show . _len) champion <> ";" <> (show . _fit) champion
  where
    getTree :: TIR -> LA.Vector Double -> SRTree Int Double
    getTree tir w = let bias   = V.head $ VS.convert w
                        consts = V.tail $ VS.convert w
                        sigm z = 1 / (1+exp(-z))

                in  case task of
                      Classification _ -> sigm $ assembleTree bias $ replaceConsts tir consts
                      ClassMult      _ -> sigm $ assembleTree bias $ replaceConsts tir consts
                      _                -> assembleTree bias $ replaceConsts tir consts

parseCLI _ _ = putStrLn "Usage: ./tir cli expmin expmax tfuncs ytfuncs errorMetric nGens nPop pc pm seed penalty trainname"

runWithCfg :: [FilePath] -> IO ()
runWithCfg [fname] = do
  cfg                         <- readConfig fname
  t0                          <- getTime Realtime
  (champion, mh, fitnessTest) <- runGP cfg
  putStrLn "Best expression (training set):\n"
  putStrLn $ prettyPrintsolution champion
  t1 <- getTime Realtime
  putStr "Fitness on the test set: "
  print $ (head . fromJust . fitnessTest) champion
  closeIfJust mh
  putStrLn $ "Total time: " ++ show (sec t1 - sec t0) ++ " secs."
  writeChampionStats cfg fitnessTest (sec t1 - sec t0) champion
runWithCfg _ = putStrLn "Usage: ./tir config filename"

runGP :: Config -> IO (Individual, Maybe Handle, Individual -> Maybe [Double])
runGP cfg@(Conf mutCfg _ algCfg cnstCfg) = do
  g <- case (_seed . _algorithmCfg) cfg of
         Nothing -> getStdGen
         Just x  -> return $ mkStdGen x
  (train, val, alldata, test, domains, image, nvars) <- processData cfg
  let
      budget       = max 5 $ min 15 $ size (snd train) `div` 10
      mutCfg'      = mutCfg{ _yfuns = filterImage image (_yfuns mutCfg), _vars = [0 .. nvars-1], _budget=budget }
      task         = _task algCfg
      measures     = _measures algCfg
      penalty      = _penaltyType cnstCfg
      cnstr        = case _evaluator cnstCfg of
                       Nothing -> const 0.0
                       Just e  -> getViolationFun e (_shapes cnstCfg) (_domains cnstCfg)
      fitnessTrain = evalTrain task False measures cnstr penalty domains (fst train) (snd train) (fst val) (snd val)
      fitnessAll   = evalTrain task True measures cnstr penalty domains (fst alldata) (snd alldata) (fst alldata) (snd alldata)
      fitnessTest  = evalTest task measures (fst test) (snd test)

      myCX OnePoint  = onepoint
      myCX UniformCX = uniformCx
      myMut GroupMutation = multiMut mutCfg'
      --myRep Generational = generational
      --mySelect (Tournament n) = tournament n
      --myFilter _ = id
      myCreate = createIndividual <$> randomTIR mutCfg'
      -- myFitness x = liftIO (print (_chromo x)) >> liftIO (print (tirToMatrix (fst train) $ selectValidTerms (_chromo x) domains)) >> pure (fitnessTrain x)
      myFitness = pure . fitnessTrain

      interpret :: Interpreter Individual
      interpret = Funs myCX myMut myCreate myFitness

      gp           = Reproduce Generational
                       [Cross OnePoint 2 (_pc algCfg) (Tournament 2) :> Mutate GroupMutation (_pm algCfg) :> Done]
      fi           = Reproduce Merge
                       [ With Feasible :> Cross OnePoint 2 (_pc algCfg) (Tournament 2) :> Mutate GroupMutation (_pm algCfg) :> Done
                       , With Infeasible :> Cross OnePoint 2 (_pc algCfg) (Tournament 2) :> Mutate GroupMutation (_pm algCfg) :> Done
                       ]

      alg          = case _algorithm algCfg of
                       GPTIR -> gp
                       SCTIR -> fi
  (logger, mh)  <- makeLogger cfg fitnessTest
  (_, champion) <- runEvolution (_gens algCfg) (_nPop algCfg) logger alg g interpret
  return (fitnessAll champion, mh, fitnessTest)

main :: IO ()
main = do
  args <- getArgs
  case args of
    ("regress":opts)    -> parseCLI Regression opts
    ("regressNL":opts)  -> parseCLI (RegressionNL 0) opts
    ("class":opts)      -> parseCLI (Classification 0) opts
    ("multiclass":opts) -> parseCLI (ClassMult 0) opts
    ("config":opts)     -> runWithCfg opts
    _                   -> putStrLn "Usage: ./tir {regress/regressNL/class/multiclass/config} args"

max_stars_repo_name: folivetti/tir

max_issues_repo_name: folivetti/tir

max_forks_repo_name: folivetti/tir


text:
module Math.AlgebraicNumbers.Overloads where
import Math.AlgebraicNumbers.Prelude
import Math.AlgebraicNumbers.Class
import Math.AlgebraicNumbers.AReal
import Math.AlgebraicNumbers.AComplex
import qualified Prelude
import qualified Data.Complex

-- Prelude functions generalized from Integral a to GCDDomain a / EuclideanDomain a

gcd :: (GCDDomain a) => a -> a -> a
gcd = gcdD

lcm :: (GCDDomain a) => a -> a -> a
lcm = lcmD

divMod :: (EuclideanDomain a) => a -> a -> (a, a)
divMod = divModD

infixl 7 `div`, `mod`

div, mod :: (EuclideanDomain a) => a -> a -> a
div = divD
mod = modD

class Sqrt a where
  sqrt :: a -> a
  maybeSqrt :: a -> Maybe a

instance Sqrt Prelude.Float where
  sqrt = Prelude.sqrt
  maybeSqrt x | x >= 0 = Just $ Prelude.sqrt x
              | otherwise = Nothing

instance Sqrt Prelude.Double where
  sqrt = Prelude.sqrt
  maybeSqrt x | x >= 0 = Just $ Prelude.sqrt x
              | otherwise = Nothing

instance (RealFloat a) => Sqrt (Data.Complex.Complex a) where
  sqrt = Prelude.sqrt
  maybeSqrt = Just . Prelude.sqrt

instance Sqrt AReal where
  sqrt = sqrtA
  maybeSqrt x | x >= 0 = Just $ sqrtA x
              | otherwise = Nothing

instance Sqrt AComplex where
  sqrt = sqrtAN
  maybeSqrt = Just . sqrtAN

-- toRational?

max_stars_repo_name: minoki/algebraic-numbers

max_issues_repo_name: minoki/algebraic-numbers

max_forks_repo_name: minoki/algebraic-numbers


text:
{-# LANGUAGE FlexibleContexts    #-}
{-# LANGUAGE OverloadedLabels    #-}
{-# LANGUAGE OverloadedStrings   #-}
{-# LANGUAGE PolyKinds           #-}
{-# LANGUAGE RankNTypes          #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TupleSections       #-}
{-# LANGUAGE CPP #-}

module Grenade.Onnx.Utils
  ( hasType
  , loadFailure
  , loadFailureReason
  , loadFailureAttr
  , loadFailureExpecting
  , guardOnnx
  , guardAttr
  , readInitializerMatrix
  , readInitializerVector
  , readInitializerTensorIntoMatrix
  , readInitializerTensorIntoVector
  , readIntAttribute
  , filterIntAttribute
  , readIntsAttribute
  , filterIntsAttribute
  , readFloatAttributeToRealNum
  , filterRealNumAttribute
  , doesNotHaveAttribute
  , hasSupportedGroup
  , hasSupportedDilations
  , hasMatchingShape
  , hasMatchingInt
  , hasMatchingRealNum
  , hasCorrectPadding
  , onnxIncorrectNumberOfInputs
  ) where

import           Control.Monad                (void)
import qualified Data.ByteString.Internal     as BS
import qualified Data.Vector.Storable         as V
import           Data.Either.Combinators
import           Data.List                    (find)
import qualified Data.Map.Strict              as Map
import           Data.Proxy

import           Foreign.ForeignPtr

import           Grenade.Onnx.OnnxLoadFailure
import           Grenade.Onnx.OnnxOperator

import           Numeric.LinearAlgebra        hiding (R, find, (<>))
import           Numeric.LinearAlgebra.Static (L, R)
import qualified Numeric.LinearAlgebra.Static as H

#if USE_FLOAT
#else
import           Data.Word                    (Word32)
import           GHC.Float                    (castWord32ToFloat, float2Double)
#endif
import           GHC.TypeLits

import           Data.ProtoLens.Labels        ()
import qualified Data.Text                    as T
import           Lens.Micro                   ((^.))
import qualified Proto.Onnx                   as P

import           Grenade.Types

hasType :: forall a. OnnxOperator a => P.NodeProto -> Proxy a -> Either OnnxLoadFailure ()
hasType node a
  | node ^. #opType `elem` expectedTypes = return ()
  | otherwise = void loadFail
  where
    nodeType = node ^. #opType
    expectedTypes = onnxOpTypeNames a
    loadFail :: Either OnnxLoadFailure (a, b, c)
    loadFail = loadFailureExpecting ("Unexpected node type: " ++ show nodeType) (Just node)

loadFailure :: String -> Maybe P.NodeProto -> Maybe P.NodeProto -> [String] -> Either OnnxLoadFailure a
loadFailure reason last curr types = Left (OnnxLoadFailure reason last curr types)

loadFailureReason :: String -> Either OnnxLoadFailure a
loadFailureReason reason = Left (OnnxLoadFailure reason Nothing Nothing [])

loadFailureAttr :: String -> P.NodeProto -> Either OnnxLoadFailure a
loadFailureAttr reason curr = loadFailure reason Nothing (Just curr) []

loadFailureExpecting :: forall a b c. OnnxOperator a => String -> Maybe P.NodeProto -> Either OnnxLoadFailure (a, b, c)
loadFailureExpecting reason curr = loadFailure reason Nothing curr (show <$> onnxOpTypeNames (Proxy :: Proxy a))

guardOnnx :: String -> Bool -> Either OnnxLoadFailure ()
guardOnnx reason cond
  | cond = return ()
  | otherwise = loadFailureReason reason

guardAttr :: T.Text -> Bool -> Either OnnxLoadFailure ()
guardAttr attribute = guardOnnx ("Attribute '" ++ show attribute ++ "' did not pass predicate")

maybeLoadFailureReason :: String -> Maybe a -> Either OnnxLoadFailure a
maybeLoadFailureReason reason = maybe (loadFailureReason reason) Right

readAttribute :: T.Text -> P.NodeProto -> Either OnnxLoadFailure P.AttributeProto
readAttribute attribute node = maybeLoadFailureReason ("Attribute '" ++ show attribute ++ "' was not found")
                                 (find ((== attribute) . (^. #name)) (node ^. #attribute))

readIntAttribute :: T.Text -> P.NodeProto -> Either OnnxLoadFailure Int
readIntAttribute attributeName node = readAttribute attributeName node >>= retrieve
  where
    retrieve attribute =
      case attribute ^. #type' of
        P.AttributeProto'INT -> Right $ fromIntegral $ attribute ^. #i
        _                    -> loadFailureReason $ "Type of attribute '" ++ show attributeName ++ "' is not int"

filterIntAttribute :: T.Text -> (Int -> Bool) -> P.NodeProto -> Either OnnxLoadFailure ()
filterIntAttribute attribute pred node = readIntAttribute attribute node >>= guardAttr attribute . pred

readFloatAttributeToRealNum :: T.Text -> P.NodeProto -> Either OnnxLoadFailure RealNum
readFloatAttributeToRealNum attributeName node = readAttribute attributeName node >>= retrieve
  where
    retrieve attribute = case (attribute ^. #type') of
#if USE_FLOAT
                           P.AttributeProto'FLOAT -> Right $ attribute ^. #f
#else
                           P.AttributeProto'FLOAT -> Right . float2Double $ attribute ^. #f
#endif
                           _                      -> loadFailureReason $
                                                       "Type of attribute '" ++ show attributeName ++ "' is not float"

filterRealNumAttribute :: T.Text -> (RealNum -> Bool) -> P.NodeProto -> Either OnnxLoadFailure ()
filterRealNumAttribute attribute pred node = readFloatAttributeToRealNum attribute node >>= guardAttr attribute . pred

readIntsAttribute :: T.Text -> P.NodeProto -> Either OnnxLoadFailure [Int]
readIntsAttribute attributeName node = readAttribute attributeName node >>= retrieve
  where
    retrieve attribute = case (attribute ^. #type') of
                           P.AttributeProto'INTS -> Right . map fromIntegral $ attribute ^. #ints
                           _                     -> loadFailureReason $
                                                      "Type of attribute '" ++ show attributeName ++ "' is not ints"

filterIntsAttribute :: T.Text -> ([Int] -> Bool) -> P.NodeProto -> Either OnnxLoadFailure ()
filterIntsAttribute attribute pred node = readIntsAttribute attribute node >>= guardAttr attribute . pred

doesNotHaveAttribute :: P.NodeProto -> T.Text -> Either OnnxLoadFailure ()
doesNotHaveAttribute node attribute = guardAttr attribute $ isLeft $ readAttribute attribute node

readInitializer :: Map.Map T.Text P.TensorProto -> T.Text -> Either OnnxLoadFailure ([Int], Vector RealNum)
readInitializer inits name = maybeLoadFailureReason ("Initializer '" ++ show name ++ "' does not exist")
                               (Map.lookup name inits) >>= retrieve
  where
    rawDataToVector :: P.TensorProto-> Maybe (Vector RealNum)
    rawDataToVector tensor
      = let (ptr, off, len)    = BS.toForeignPtr (tensor ^. #rawData)
            vec                = V.unsafeFromForeignPtr (castForeignPtr ptr) (div off 4) (div len 4)
#if USE_FLOAT
        in  if len == 0 then Nothing else Just vec
#else
        in if len == 0 then Nothing else Just $ V.map (float2Double . castWord32ToFloat) vec
#endif

    readFloatData :: P.TensorProto -> Either OnnxLoadFailure (Vector RealNum)
    readFloatData tensor
      = let floatData = tensor ^. #floatData
        in  if null floatData then loadFailureReason ("Failed reading raw_data field of initializer " ++ show name)
#if USE_FLOAT
                              else Right $ vector floatData
#else
                              else Right $ vector $ map float2Double floatData
#endif

    getData :: P.TensorProto -> Either OnnxLoadFailure (Vector RealNum)
    getData tensor = case rawDataToVector tensor of
      Just vec -> Right vec
      Nothing   -> readFloatData tensor

    retrieve tensor = let vec = getData tensor
                      in  case toEnum (fromIntegral (tensor ^. #dataType)) of
                            P.TensorProto'FLOAT -> (map fromIntegral $ tensor ^. #dims, ) <$> vec
                            _                   -> loadFailureReason $ "Type of attribute '" ++ show name ++ "' is not float"

readInitializerMatrix :: (KnownNat r, KnownNat c) => Map.Map T.Text P.TensorProto -> T.Text -> Either OnnxLoadFailure (L r c)
readInitializerMatrix inits name = do
  initializer <- readInitializer inits name
  maybeLoadFailureReason ("Failed to read initializer '" ++ show name ++ "' as matrix") (readMatrix initializer)

readInitializerTensorIntoMatrix :: (KnownNat r, KnownNat c) => Map.Map T.Text P.TensorProto -> T.Text -> Either OnnxLoadFailure (L r c)
readInitializerTensorIntoMatrix inits name = do
  initializer <- readInitializer inits name
  maybeLoadFailureReason ("Failed to read initializer tensor '" ++ show name ++ "' as matrix") (readTensorIntoMatrix initializer)

readInitializerTensorIntoVector :: (KnownNat c) => Map.Map T.Text P.TensorProto -> T.Text -> Either OnnxLoadFailure (R c)
readInitializerTensorIntoVector inits name = do
  initializer <- readInitializer inits name
  maybeLoadFailureReason ("Failed to read initializer tensor '" ++ show name ++ "' as matrix into vector") (readMatrixToVector initializer)

readInitializerVector :: KnownNat r => Map.Map T.Text P.TensorProto -> T.Text -> Either OnnxLoadFailure (R r)
readInitializerVector inits name = do
  initializer <- readInitializer inits name
  maybeLoadFailureReason ("Failed to read initializer '" ++ show name ++ "' as vector") (readVector initializer)

vector' :: forall r . (KnownNat r) => Vector RealNum -> Maybe (R r)
vector' = H.create
{-# INLINE vector' #-}

matrix' :: forall r c . (KnownNat r, KnownNat c) => Vector RealNum -> Maybe (L r c)
matrix' vec = let cols = fromIntegral $ natVal (Proxy :: Proxy c) in H.create $ reshape cols vec
{-# INLINE matrix' #-}

readMatrix :: forall r c . (KnownNat r, KnownNat c) => ([Int], Vector RealNum) -> Maybe (L r c)
readMatrix ([rows, cols], vals)
  | neededRows == rows && neededCols == cols = matrix' vals
  where
    neededRows = fromIntegral $ natVal (Proxy :: Proxy r)
    neededCols = fromIntegral $ natVal (Proxy :: Proxy c)
readMatrix _ = Nothing

readTensorIntoMatrix :: forall r c . (KnownNat r, KnownNat c) => ([Int], Vector RealNum) -> Maybe (L r c)
readTensorIntoMatrix (rows : cols, vals)
  | neededRows == rows && neededCols == product cols = matrix' vals
  where
    neededRows = fromIntegral $ natVal (Proxy :: Proxy r)
    neededCols = fromIntegral $ natVal (Proxy :: Proxy c)
readTensorIntoMatrix _ = Nothing

readVector :: forall r . KnownNat r => ([Int], Vector RealNum) -> Maybe (R r)
readVector ([rows], vals)
  | rows == neededRows = vector' vals
  where
    neededRows = fromIntegral $ natVal (Proxy :: Proxy r)
readVector _ = Nothing

readMatrixToVector :: forall r . KnownNat r => ([Int], Vector RealNum) -> Maybe (R r)
readMatrixToVector (dims, vals)
  | product dims == neededRows = vector' vals
  where
    neededRows = fromIntegral $ natVal (Proxy :: Proxy r)
readMatrixToVector _ = Nothing

hasSupportedGroup :: P.NodeProto -> Either OnnxLoadFailure ()
hasSupportedGroup node = (node `doesNotHaveAttribute` "group")
                      <> (filterIntAttribute "group" (== 1) node)

hasSupportedDilations :: P.NodeProto -> Either OnnxLoadFailure ()
hasSupportedDilations node = (node `doesNotHaveAttribute` "dilations")
                          <> (filterIntsAttribute "dilations" (all (==1)) node)

hasMatchingShape :: P.NodeProto -> T.Text -> [Integer] -> Either OnnxLoadFailure ()
hasMatchingShape node attribute dims = filterIntsAttribute attribute (== map fromIntegral dims) node

hasCorrectPadding :: (KnownNat padLeft, KnownNat padRight, KnownNat padTop, KnownNat padBottom)
                  => P.NodeProto -> Proxy padLeft -> Proxy padRight -> Proxy padTop -> Proxy padBottom -> Either OnnxLoadFailure ()
hasCorrectPadding node ppl ppr ppt ppb
  = let left   = fromIntegral $ natVal ppl
        right  = fromIntegral $ natVal ppr
        top    = fromIntegral $ natVal ppt
        bottom = fromIntegral $ natVal ppb
        paddings = case readIntsAttribute "pads" node of
                     Right ps -> ps
                     Left  _  -> [0, 0, 0, 0]
     in case paddings of
          [left', top', right', bottom'] -> guardOnnx "Padding dimensions mismatch" (left == left' && top == top' && right == right' && bottom == bottom')
          _                              -> loadFailureReason "Incorrect padding dimensions"

hasMatchingRealNum :: KnownSymbol a => P.NodeProto -> Proxy a -> T.Text -> Either OnnxLoadFailure ()
hasMatchingRealNum node a x = do
  let a' = (read $ symbolVal a) :: RealNum
  x' <- readFloatAttributeToRealNum x node
#if USE_FLOAT
  guardOnnx ("Value mismatch for attribute " ++ T.unpack x ++ " of type float") (x' == a')
#else
  guardOnnx ("Value mismatch for attribute " ++ T.unpack x ++ " of type double") (x' == a')
#endif

hasMatchingInt :: KnownNat a => P.NodeProto -> Proxy a -> T.Text -> Either OnnxLoadFailure ()
hasMatchingInt node a x = do
  let a' = fromIntegral $ natVal a
  x' <- readIntAttribute x node
  guardOnnx ("Value mismatch for attribute " ++ T.unpack x ++ " of type int") (x' == a')

onnxIncorrectNumberOfInputs :: Either OnnxLoadFailure a
onnxIncorrectNumberOfInputs = loadFailureReason "Incorrect number of inputs"

max_stars_repo_name: th-char/grenade

max_issues_repo_name: th-char/grenade

max_forks_repo_name: th-char/grenade


text:
{-# LANGUAGE TypeOperators #-}

--------------------------------------------------------------------------------
{- |
Module      :  Numeric.LinearAlgebra.Data
Copyright   :  (c) Alberto Ruiz 2015
License     :  BSD3
Maintainer  :  Alberto Ruiz
Stability   :  provisional

This module provides functions for creation and manipulation of vectors and matrices, IO, and other utilities.

-}
--------------------------------------------------------------------------------

module Numeric.LinearAlgebra.Data(

    -- * Elements
    R,C,I,Z,type(./.),

    -- * Vector
    {- | 1D arrays are storable vectors directly reexported from the vector package.
    -}

    fromList, toList, (|>), vector, range, idxs,

    -- * Matrix

    {- | The main data type of hmatrix is a 2D dense array defined on top of
         a storable vector. The internal representation is suitable for direct
         interface with standard numeric libraries.
    -}

    (><), matrix, tr, tr',

    -- * Dimensions

    size, rows, cols,

    -- * Conversion from\/to lists

    fromLists, toLists,
    row, col,

    -- * Conversions vector\/matrix

    flatten, reshape, asRow, asColumn,
    fromRows, toRows, fromColumns, toColumns,

    -- * Indexing

    atIndex,
    Indexable(..),

    -- * Construction
    scalar, Konst(..), Build(..), assoc, accum, linspace,  -- ones, zeros,

    -- * Diagonal
    ident, diag, diagl, diagRect, takeDiag,

    -- * Vector extraction
    subVector, takesV, vjoin,

    -- * Matrix extraction
    Extractor(..), (??),

    (?), (¿), fliprl, flipud,

    subMatrix, takeRows, dropRows, takeColumns, dropColumns,

    remap,

    -- * Block matrix
    fromBlocks, (|||), (===), diagBlock, repmat, toBlocks, toBlocksEvery,

    -- * Mapping functions
    conj, cmap, cmod,

    step, cond,

    -- * Find elements
    find, maxIndex, minIndex, maxElement, minElement,
    sortVector, sortIndex,

    -- * Sparse
    AssocMatrix, toDense,
    mkSparse, mkDiagR, mkDense,

    -- * IO
    disp,
    loadMatrix, loadMatrix', saveMatrix,
    latexFormat,
    dispf, disps, dispcf, format,
    dispDots, dispBlanks, dispShort,
-- * Element conversion
    Convert(..),
    roundVector,
    fromInt,toInt,fromZ,toZ,
    -- * Misc
    arctan2,
    separable,
    fromArray2D,
    module Data.Complex,
    Mod,
    Vector, Matrix, GMatrix, nRows, nCols

) where

import Internal.Vector
import Internal.Vectorized
import Internal.Matrix hiding (size)
import Internal.Element
import Internal.IO
import Internal.Numeric
import Internal.Container
import Internal.Util hiding ((&))
import Data.Complex
import Internal.Sparse
import Internal.Modular

max_stars_repo_name: schnecki/hmatrix-float

max_issues_repo_name: schnecki/hmatrix-float

max_forks_repo_name: schnecki/hmatrix-float


text:
module Y2016.M12.D14.Solution where

import Data.Array (elems, listArray)
import Data.Complex

-- import available via 1HaskellADay git repository

import Data.Matrix

{--
I'm thinking about quantum computation. IBM has released a 5-qubit computer
for public experimentation. So, let's experiment.

One way to go about that is to dive right in, so, yes, if you wish: dive
right in.

Another approach is to comprehend the maths behind quantum computation.

So, let's look at that.

I was going to bewail that Shor's prime factors algorithm needs 7 qubits to
work, but, NEWSFLASH! IBM has added Shor's algorithm to their API, so ...

CANCEL BEWAILMENT.

*ahem*

Moving on.

First, let's look at qubits. Qubits are 'bra-ket'ted numbers (ket numbers)
with the representation

|0> =  | 1 |   or    |1> = | 0 |
       | 0 |               | 1 |

OOH! MATRICES!

exercise 1. Represent ket0 and ket1 states as matrices in Haskell
--}

data Qubit = Qbit (Matrix (Complex Float))  --- where Float is 1 or 0 but okay

instance Show Qubit where
   show (Qbit mat) = '|':complexify mat ++ ">"

complexify :: Matrix (Complex Float) -> String
complexify = showC . determinant . cross (fromLists [[0, 1]])

showC :: (Eq a, Show a, Num a) => Complex a -> String
showC (x :+ 0) = show x
showC (0 :+ x) = (if x == 1 then "" else show x) ++ "i"
showC c = show c

ket0, ket1 :: Qubit
ket0 = Qbit (fromLists [[1],[0]])
ket1 = Qbit (fromLists [[0],[1]])

{--
It MAY be helpful to have a show-instance of a qubit that abbreviates the
complex number to something more presentable. Your choice.

A qubit state is most-times in a super-position of |0> or |1> and we represent
that as

|ψ> = α|0> + β|1>

And we KNOW that |α|² + |β|² = 1

YAY! Okay. Whatever.

So, we have a qubit at |0>-state and we want to flip it to |1>-state, or vice
versa. How do we do that?

We put it through a Pauli X gate

The Pauli X operator is =  | 0 1 |
                           | 1 0 |

That is to say, zero goes to 1 and 1 goes to zero.

excercise 2: represent the Pauli X, Y, and Z operators
--}

data PauliOperator = POp (Matrix (Complex Float))

data Cnum = C (Complex Float)

instance Show Cnum where show (C c) = showC c

printPauli :: PauliOperator -> IO ()
printPauli (POp mat) =
   pprint (M (listArray ((1,1),(2,2)) (map C (elems (matrix mat)))))

pauliX, pauliY, pauliZ :: PauliOperator
pauliX = POp (fromLists [[0,1],[1,0]])
pauliY = POp (fromLists [[0,0 :+ (-1)],[0 :+ 1,0]])
pauliZ = POp (fromLists [[1,0],[0,-1]])

{--
*Y2016.M12.D14.Solution> printPauli pauliX
Matrix 2x2
| 0.0 1.0 |
| 1.0 0.0 |
--}

-- exercise 3: rotate the qubits ket0 and ket1 through the pauliX operator
-- (figure out what that means). The intended result is:

-- X|0> = |1> and X|1> = |0>

-- what are your results?

rotate :: PauliOperator -> Qubit -> Qubit
rotate (POp p) (Qbit q) = Qbit (cross p q)

{--
*Y2016.M12.D14.Solution> rotate pauliX ket0 ~> |1.0>
*Y2016.M12.D14.Solution> rotate pauliX ket1 ~> |0.0>

BAM!
--}


text:
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE ScopedTypeVariables   #-}
module Bio.ChromVAR.Utils
    ( weight
    , Whitening(..)
    , covariance
    , whiten
    ) where

import qualified Data.Vector.Unboxed as U
import qualified Data.Vector.Storable as S
import Statistics.Sample (mean)
import Data.Matrix.Static.LinearAlgebra
import qualified Data.Matrix.Static.Dense as D
import qualified Data.Matrix.Static.Generic as D
import Data.Singletons (SingI)

weight :: (Double, Double) -> (Double, Double) -> Double
weight (x1, y1) (x2, y2) = gaussian $ sqrt $ (x1 - x2)**2 + (y1 - y2)**2
  where
    gaussian d = let sigma = 0.1 in exp $ negate $ d**2 / (2 * sigma**2)
    --gaussian d = d
{-# INLINE weight #-}

data Whitening = ZCA
               | Cholesky

-- | Rows are samples, columns are features.
whiten :: (SingI n, SingI m) => Whitening -> Matrix n m Double -> Matrix n m Double
whiten method mat = case method of
    Cholesky -> D.transpose $ inverse (cholesky cov) @@ D.transpose mat
  where
    cov = covariance mat

covariance :: forall n m. (SingI n, SingI m) => Matrix n m Double -> Matrix m m Double
covariance mat = D.map (/n) $ D.transpose cs @@ cs
  where
    n = fromIntegral $ D.rows mat - 1
    cs = D.fromColumns $ map f $ D.toColumns mat :: Matrix n m Double
    f x = let m = mean x in S.map (subtract m) x
{-# INLINE covariance #-}

-- | Mahalanobis or ZCA whitening
zca :: U.Vector (Double, Double) -> U.Vector (Double, Double)
zca points = U.map (`dot'` w) $ points -- U.zip xs' ys'
  where
    w = t $ evec `dot` d `dot` t evec
    d = (1 / sqrt (eval1 + eps), 0, 0, 1 / sqrt (eval2 + eps))
    (evec, (eval1, eval2)) = eigendecomposition (sx, sxy, sxy, sy)
    eps = 10e-6
    sx = U.sum (U.zipWith (*) xs' xs') / n
    sy = U.sum (U.zipWith (*) ys' ys') / n
    sxy = U.sum (U.zipWith (*) xs' ys') / n
    xs' = U.map (\x -> x - mx) xs
    ys' = U.map (\x -> x - my) ys
    mx = mean xs
    my = mean ys
    n = fromIntegral $ U.length points - 1
    (xs, ys) = U.unzip points
{-# INLINE zca #-}

type Mat2x2 = (Double, Double, Double, Double)

-- | Eigendecomposition of 2 x 2 symetric matrix
eigendecomposition :: Mat2x2
                   -> ( Mat2x2   -- ^ Eigenvectors in the columns
                      , (Double, Double) ) -- ^ Eigenvalues
eigendecomposition (a,b,_,d) =
    let (x1, y1) = evector eval1
        (x2, y2) = evector eval2
    in ((x1, x2, y1, y2), (eval1, eval2))
  where
    eval1 = (a + d + s) / 2
    eval2 = (a + d - s) / 2
    s = sqrt $ (a - d)**2 + 4 * b**2
    evector lambda =
        let denominator = sqrt $ b**2 + (lambda - d)**2
        in ((lambda - d) / denominator, b / denominator)
{-# INLINE eigendecomposition #-}

dot :: Mat2x2 -> Mat2x2 -> Mat2x2
dot (a,b,c,d) (a',b',c',d') =
    ( a * a' + b * c'
    , a * b' + b * d'
    , c * a' + d * c'
    , c * b' + d * d' )
{-# INLINE dot #-}

dot' :: (Double, Double) -> Mat2x2 -> (Double, Double)
dot' (x,y) (a,b,c,d) = (x * a + y * c, x * b + y * d)
{-# INLINE dot' #-}

t :: Mat2x2 -> Mat2x2
t (a,b,c,d) = (a,c,b,d)
{-# INLINE t #-}

max_stars_repo_name: kaizhang/ChromVAR

max_issues_repo_name: kaizhang/ChromVAR

max_forks_repo_name: kaizhang/ChromVAR


text:
{-# LANGUAGE RecordWildCards #-}

module Main where

import           Control.DeepSeq                  (force)
import           Control.Exception                (evaluate)
import           Control.Monad                    (replicateM)
import           Data.List                        (minimumBy)
import           Data.Ord                         (comparing)
import qualified Data.Text                        as T
import qualified Data.Vector.Unboxed              as V
import           GHC.Base                         (String)
import           Options.Applicative              (Parser (..), auto,
                                                   execParser, fullDesc, help,
                                                   helper, info, long, metavar,
                                                   option, progDesc,
                                                   showDefault, strOption,
                                                   value, (<**>))
import qualified Statistics.Distribution.StudentT as S
import qualified Statistics.Function              as S
import qualified Statistics.Sample                as S
import qualified Statistics.Test.StudentT         as S
import qualified Statistics.Test.Types            as ST
import           System.IO                        (hGetContents)
import           System.Process                   (CreateProcess (..),
                                                   StdStream (..), shell,
                                                   waitForProcess,
                                                   withCreateProcess)
import           Text.Printf                      (printf)

type Samples = V.Vector Double

data Args = Args
  { old   :: T.Text
  , new   :: T.Text
  , forks :: Int
  }

data Diff = Diff
  { meanPct :: Double
  , oldMean :: Double
  , oldMin  :: Double
  , oldMax  :: Double
  , newMin  :: Double
  , newMax  :: Double
  , newMean :: Double
  , test    :: Maybe (ST.Test S.StudentT) }
  deriving (Show)

args :: Parser Args
args = Args
    <$> strOption
        ( long "old"
        <> metavar "OLD"
        <> help "Command that returns samples for \"old\"" )
    <*> strOption
        ( long "new"
        <> metavar "NEW"
        <> help "Command that returns samples for \"new\"" )
    <*> option auto
        ( long "forks"
        <> help "How often the commands will be invoked to gather samples"
        <> showDefault
        <> value 1
        <> metavar "FORKS" )

execProc :: String -> IO Samples
execProc cmd =
  withCreateProcess createProc $ \_ (Just out) _ ph -> do
    samples <- textToSamples <$> hGetContents out
    evaluate $ force samples
    ex <- waitForProcess ph
    pure samples
  where
    textToSamples = V.fromList . map read . lines
    createProc = (shell cmd) { std_out = CreatePipe }

calcDiff :: Samples -> Samples -> Diff
calcDiff oldSamples newSamples = Diff{..}
  where
    (oldMin, oldMax) = S.minMax oldSamples
    (newMin, newMax) = S.minMax newSamples
    oldMean = S.mean oldSamples
    newMean = S.mean newSamples
    oldStdDev = S.stdDev oldSamples
    newStdDev = S.stdDev newSamples
    meanPct = percDiff oldMean newMean
    test = S.welchTTest S.BGreater oldSamples newSamples

percDiff :: Double -> Double -> Double
percDiff x y = (abs (x - y) / ((x + y) / 2)) * 100

cmpPerf :: Args -> IO ()
cmpPerf (Args old new forks) = do
  listOfOldSamples <- replicateM forks $ do
    putStrLn ("Gathering metrics for: " <> T.unpack old)
    execProc $ T.unpack old
  listOfNewSamples <- replicateM forks $ do
    putStrLn ("Gathering metrics for: " <> T.unpack new)
    execProc $ T.unpack new
  let
    bestOldSamples = minimumBy (comparing S.mean) listOfOldSamples
    bestNewSamples = minimumBy (comparing S.mean) listOfNewSamples
    diff = calcDiff bestOldSamples bestNewSamples
  printf   "%.2f%% mean difference\n" (meanPct diff)
  putStrLn " Cmd |       Mean |        Min |        Max |"
  printf   " Old | %10.3f | %10.3f | %10.3f |\n" (oldMean diff) (oldMin diff) (oldMax diff)
  printf   " New | %10.3f | %10.3f | %10.3f |\n" (newMean diff) (newMin diff) (newMax diff)

main :: IO ()
main = cmpPerf =<< execParser opts
  where
    opts = info (args <**> helper)
      ( fullDesc
      <> progDesc "Compare the (runtime) values returned by two different programs" )

max_stars_repo_name: mfussenegger/cmp-perf

max_issues_repo_name: mfussenegger/cmp-perf

max_forks_repo_name: mfussenegger/cmp-perf


text:
{-# LANGUAGE OverloadedStrings #-}

module ML.Gym.Space where

import           ML.Gym.DType
import           ML.Gym.Range
import           ML.Gym.Shape
import           ML.Gym.Util
import           ML.Gym.Value

import qualified Control.Exception        as E
import           Control.Monad            (void)
import qualified CPython                  as Py
import qualified CPython.Constants        as Py
import qualified CPython.Protocols.Object as Py
import qualified CPython.Types            as Py
import qualified CPython.Types.Exception  as Py
import qualified CPython.Types.Module     as Py
import qualified Data.ByteString          as BS
import           Data.Complex             (Complex)
import           Data.Maybe               (fromMaybe)
import qualified Data.Text                as T
import           System.Environment       (getArgs, getProgName)
import           System.Exit
import           System.IO
import           System.Random            (randomRIO)

import           Debug.Trace

class SpaceSize a where
  spaceSize :: a -> Integer     -- Rename to something better

data GymSpace = GymSpaceDiscrete Integer        -- ^ Range: 0...(x-1)
              | GymSpaceBox GymRange GymShape   -- ^ 1D, 2D, or 3D matrix. Range sizes and shape must coincide.
              | GymSpaceTuple [GymSpace]        -- ^ Tuple of Spaces.
              deriving (Show, Eq)

instance SpaceSize GymSpace where
  spaceSize (GymSpaceDiscrete nr)                 = nr
  spaceSize (GymSpaceBox _ (GymShape1D x1))       = x1
  spaceSize (GymSpaceBox _ (GymShape2D x1 x2))    = x1 * x2
  spaceSize (GymSpaceBox _ (GymShape3D x1 x2 x3)) = x1 * x2 * x3
  spaceSize (GymSpaceTuple spaces)                = product (map spaceSize spaces)

toAction :: GymSpace -> Integer -> IO Py.SomeObject
toAction (GymSpaceDiscrete nr) idx
  | idx < 0 || idx > (nr - 1) = error "Action index out of range"
  | otherwise = Py.toObject <$> Py.toInteger idx
toAction (GymSpaceBox _ (GymShape1D x1)) idx = error "Box (probably continuous) actions currently not supported"
toAction (GymSpaceBox _ (GymShape2D x1 x2)) idx = error "Box (probably continuous) actions currently not supported"
toAction (GymSpaceBox _ (GymShape3D x1 x2 x3)) idx = error "Box (probably continuous) actions currently not supported"
toAction (GymSpaceTuple spaces) idx
  | idx >= fromIntegral (length cats) = error "action index is out of range"
  | otherwise = Py.toObject <$> (mapM (fmap Py.toObject . Py.toInteger) (cats !! fromIntegral idx) >>= Py.toList)
  where
    sizes = map spaceSize spaces
    cats = fromCats sizes

fromCats :: [Integer] -> [[Integer]]
fromCats []     = []
fromCats [x]    = map return [0..x-1]
fromCats (x:xs) = concatMap (\cat -> map (cat ++) (fromCats xs)) (fromCats [x])

fromIdx :: [Integer] -> Integer -> [Integer]
fromIdx cats idx =
  -- trace ("rev: " ++ show (reverse $ scanl1 (*) cats))
  zipWith toCat cats (reverse $ scanl1 (*) cats)
  where toCat i c =
          -- trace ("nr: " ++ show nr)
          -- trace ("c: " ++ show c)
          -- trace ("i: " ++ show i)
          (idx `div` (max 1 (c `div` 2))) `mod` i
        nr = product cats

test :: [Integer] -> Bool
test cats = do

  all (\i ->(fromCats cats !! i) == fromIdx cats (fromIntegral i)) [0..product (map fromIntegral cats)-1]

getGymRangeFromSpace :: GymSpace -> GymRange
getGymRangeFromSpace (GymSpaceDiscrete nr)  = GymRange GymInteger (GymScalar 0) (GymScalar $ fromIntegral nr-1)
getGymRangeFromSpace (GymSpaceBox rng _)    = rng
getGymRangeFromSpace (GymSpaceTuple spaces) =
  trace ("from: " ++ show spaces)
  trace ("to: " ++ show (combineRanges $ map getGymRangeFromSpace spaces))
  combineRanges $ map getGymRangeFromSpace spaces

toGymSpace :: Py.Object self => self -> IO GymSpace
toGymSpace space = mkSpaceRepr >>= toGymSpace'
  where
    mkSpaceRepr = fromMaybe (error "could not get gym action space") <$> python (Py.callMethodArgs space "__repr__" [] >>= pyToText)
    toGymSpace' txt
      | T.isPrefixOf "Discrete" txt = do
        n <- fromMaybe (error "Gym discrete space 'n' not recognized!") <$> python (Py.toUnicode "n" >>= Py.getAttribute space >>= pyToInteger)
        return $ GymSpaceDiscrete n
      | T.isPrefixOf "Box" txt = do
        dType <- python $ Py.toUnicode "dtype" >>= Py.getAttribute space >>= getGymDType
        range <- fromMaybe (error $ "could not get gym range: " ++ show dType ++ ", " ++ T.unpack txt) <$> python (getGymRange dType space)
        shape <- fromMaybe (error "could not get gym shape") <$> python (getGymShape space)
        return $ GymSpaceBox range shape
      | T.isPrefixOf "Tuple" txt = do
          list <- fromMaybe (error "Gym discrete space 'n' not recognized!") <$> python (Py.toUnicode "spaces" >>= Py.getAttribute space >>= pyToList)
          spaces <- mapM toGymSpace list
          return $ GymSpaceTuple spaces
      | otherwise = do
        putStrLn "Sample space: "
        python (Py.callMethodArgs space "sample" [] >>= flip Py.print stdout)
        error $ "ERROR: Gym space '" ++ T.unpack txt ++ "' not supported yet!"

max_stars_repo_name: schnecki/gym-haskell

max_issues_repo_name: schnecki/gym-haskell

max_forks_repo_name: schnecki/gym-haskell


text:
import Prelude
import Control.Applicative
import Control.Arrow
import Control.Monad
import Data.Array
import Data.Bits
import Data.Bool
import Data.Char
import Data.Complex
import Data.Dynamic
import Data.Either
import Data.Eq
import Data.Fixed
import Data.Function
import Data.Graph
import Data.Int
import Data.Ix
import Data.List
import Data.Maybe
import Data.Monoid
import Data.Ord
import Data.Ratio
import Data.Tree
import Data.Tuple
import Data.Typeable
import Data.Word
import Text.Printf

' '?y=y
x?' '=x
_?_='~'
(x:c)??(y:d)=x?y:c??d
x??y=x++y
""#[]=[(0,"")]
x#y=[(c+1,a:d)|a:b<-[x],a/='~',(c,d)<-b#y]++[c|a:b<-[y],c<-x??a#b]
g a=snd$minimum$(#)""=<<permutations a

test input ans=
  let output = g input
  in do
    print ("output: " ++ output)
    print ("answer: " ++ ans)

main=do
  test ["hac m", "k e"] "hackme"
  test ["b", "a"] "ab"
  test [] ""
  test ["rst r", "rsy rsytsr", "a"] "rst rrsyarsytsr"

max_stars_repo_path: GoogleCTF/2019 Quals/code_golf/pos.hs

max_stars_repo_name: PurpEth/solved-hacking-problem

max_issues_repo_path: GoogleCTF/2019 Quals/code_golf/pos.hs

max_issues_repo_name: PurpEth/solved-hacking-problem

max_forks_repo_path: GoogleCTF/2019 Quals/code_golf/pos.hs

max_forks_repo_name: PurpEth/solved-hacking-problem


text:
module Forward (
    forwardClassification,
    forwardRegression,
    forward
  ) where

import Numeric.LinearAlgebra
import Common
import ActivFunc

forwardClassification :: (Matrix R -> Matrix R) -> [Matrix R] -> Matrix R -> Matrix R
forwardClassification f ws v = softmaxC $ forwardRegression f ws v

forwardRegression :: (Matrix R -> Matrix R) -> [Matrix R] -> Matrix R -> Matrix R
forwardRegression f ws v = forward id (last ws) nv
  where nv = forwards f (init ws) v

forwards :: (Matrix R -> Matrix R) -> [Matrix R] -> Matrix R -> Matrix R
forwards f ws v = foldl (flip $ forward f) v ws

forward :: (Matrix R -> Matrix R) -> Matrix R -> Matrix R -> Matrix R
forward f w v = f $ w <> inputWithBias v

max_stars_repo_name: pupuu/deep-neuralnet

max_issues_repo_name: pupuu/deep-neuralnet

max_forks_repo_name: pupuu/deep-neuralnet