text:
{-# LANGUAGE FlexibleInstances #-}
{-# OPTIONS_GHC -Wall #-}

-- |
-- Module      :  Numeric.NLOPT
-- Copyright   :  (c) Matthew Peddie 2017
-- License     :  BSD3
-- Maintainer  :  Matthew Peddie <mpeddie@gmail.com>
-- Stability   :  provisional
-- Portability :  GHC
--
-- This module provides a high-level, @hmatrix@-compatible interface to
-- the <http://ab-initio.mit.edu/wiki/index.php/NLopt NLOPT> library by
-- Steven G. Johnson.
--
-- = Documentation
--
-- Most non-numerical details are documented, but for specific
-- information on what the optimization methods do, how constraints are
-- handled, etc., you should consult:
--
--   * The <http://ab-initio.mit.edu/wiki/index.php/NLopt_Introduction NLOPT introduction>
--
--   * The <http://ab-initio.mit.edu/wiki/index.php/NLopt_Reference NLOPT reference manual>
--
--   * The <http://ab-initio.mit.edu/wiki/index.php/NLopt_Algorithms NLOPT algorithm manual>
--
-- = Example program
--
-- The following interactive session example uses the Nelder-Mead simplex
-- algorithm, a derivative-free local optimizer, to minimize a trivial
-- function with a minimum of 22.0 at @(0, 0)@.
--
-- >>> import Numeric.LinearAlgebra ( dot, fromList )
-- >>> let objf x = x `dot` x + 22                         -- define objective
-- >>> let stop = ObjectiveRelativeTolerance 1e-6 :| []    -- define stopping criterion
-- >>> let algorithm = NELDERMEAD objf [] Nothing          -- specify algorithm
-- >>> let problem = LocalProblem 2 stop algorithm         -- specify problem
-- >>> let x0 = fromList [5, 10]                           -- specify initial guess
-- >>> minimizeLocal problem x0
-- Right (Solution {solutionCost = 22.0, solutionParams = [0.0,0.0], solutionResult = FTOL_REACHED})
module Numeric.NLOPT
  ( -- * Specifying the objective function
    Objective,
    ObjectiveD,
    Preconditioner,

    -- * Specifying the constraints

    -- ** Bound constraints
    Bounds (..),

    -- ** Nonlinear constraints

    --
    -- $nonlinearconstraints

    -- *** Constraint functions
    ScalarConstraint,
    ScalarConstraintD,
    VectorConstraint,
    VectorConstraintD,

    -- *** Constraint types
    Constraint (..),
    EqualityConstraint (..),
    InequalityConstraint (..),

    -- *** Collections of constraints
    EqualityConstraints,
    EqualityConstraintsD,
    InequalityConstraints,
    InequalityConstraintsD,

    -- * Stopping conditions

    --
    -- $nonempty
    StoppingCondition (..),
    NonEmpty (..),

    -- * Additional configuration
    RandomSeed (..),
    Population (..),
    VectorStorage (..),
    InitialStep (..),

    -- * Minimization problems

    -- ** Local minimization
    LocalAlgorithm (..),
    LocalProblem (..),
    minimizeLocal,

    -- ** Global minimization
    GlobalAlgorithm (..),
    GlobalProblem (..),
    minimizeGlobal,

    -- ** Minimization by augmented Lagrangian
    AugLagAlgorithm (..),
    AugLagProblem (..),
    minimizeAugLag,

    -- ** Results
    Solution (..),
    N.Result (..),
  )
where

import Control.Exception (Exception)
import qualified Control.Exception as Ex
import Data.Foldable (traverse_)
import Data.List.NonEmpty (NonEmpty (..))
import Data.Typeable (Typeable)
import qualified Data.Vector.Storable as V
import Numeric.LinearAlgebra as HM
import qualified Numeric.Optimization.NLOPT.Bindings as N
import System.IO.Unsafe (unsafePerformIO)

{- Function wrapping for the immutable HMatrix interface -}
wrapScalarFunction :: (Vector Double -> Double) -> N.ScalarFunction ()
wrapScalarFunction f params _ _ = return $ f params

wrapScalarFunctionD ::
  (Vector Double -> (Double, Vector Double)) ->
  N.ScalarFunction ()
wrapScalarFunctionD f params grad _ = do
  case grad of
    Nothing -> return ()
    Just g -> V.copy g usergrad
  return result
  where
    (result, usergrad) = f params

wrapVectorFunction ::
  (Vector Double -> Word -> Vector Double) ->
  Word ->
  N.VectorFunction ()
wrapVectorFunction f n params vout _ _ = V.copy vout $ f params n

wrapVectorFunctionD ::
  (Vector Double -> Word -> (Vector Double, Matrix Double)) ->
  Word ->
  N.VectorFunction ()
wrapVectorFunctionD f n params vout jac _ = do
  V.copy vout result
  case jac of
    Nothing -> return ()
    Just j -> V.copy j (HM.flatten userjac)
  where
    (result, userjac) = f params n

wrapPreconditionerFunction ::
  (Vector Double -> Vector Double -> Vector Double) ->
  N.PreconditionerFunction ()
wrapPreconditionerFunction f params v vpre _ = V.copy vpre (f params v)

{- Objective functions -}

-- | An objective function that calculates the objective value at the
-- given parameter vector.
type Objective =
  -- | Parameter vector
  Vector Double ->
  -- | Objective function value
  Double

-- | An objective function that calculates both the objective value
-- and the gradient of the objective with respect to the input
-- parameter vector, at the given parameter vector.
type ObjectiveD =
  -- | Parameter vector
  Vector Double ->
  -- | (Objective function value, gradient)
  (Double, Vector Double)

-- | A preconditioner function, which computes @vpre = H(x) v@, where
-- @H@ is the Hessian matrix: the positive semi-definite second
-- derivative at the given parameter vector @x@, or an approximation
-- thereof.
type Preconditioner =
  -- | Parameter vector @x@
  Vector Double ->
  -- | Vector @v@ to precondition at @x@
  Vector Double ->
  -- | Preconditioned vector @vpre@
  Vector Double

data ObjectiveFunction f
  = MinimumObjective f
  | PreconditionedMinimumObjective Preconditioner f

applyObjective :: N.Opt -> ObjectiveFunction Objective -> IO N.Result
applyObjective opt (MinimumObjective f) =
  N.set_min_objective opt (wrapScalarFunction f) ()
applyObjective opt (PreconditionedMinimumObjective p f) =
  N.set_precond_min_objective
    opt
    (wrapScalarFunction f)
    (wrapPreconditionerFunction p)
    ()

applyObjectiveD :: N.Opt -> ObjectiveFunction ObjectiveD -> IO N.Result
applyObjectiveD opt (MinimumObjective f) =
  N.set_min_objective opt (wrapScalarFunctionD f) ()
applyObjectiveD opt (PreconditionedMinimumObjective p f) =
  N.set_precond_min_objective
    opt
    (wrapScalarFunctionD f)
    (wrapPreconditionerFunction p)
    ()

{- Constraint functions -}

-- | A constraint function which returns @c(x)@ given the parameter
-- vector @x@.  The constraint will enforce that @c(x) == 0@ (equality
-- constraint) or @c(x) <= 0@ (inequality constraint).
type ScalarConstraint =
  -- | Parameter vector @x@
  Vector Double ->
  -- | Constraint violation (deviation from 0)
  Double

-- | A constraint function which returns @c(x)@ given the parameter
-- vector @x@ along with the gradient of @c(x)@ with respect to @x@ at
-- that point.  The constraint will enforce that @c(x) == 0@ (equality
-- constraint) or @c(x) <= 0@ (inequality constraint).
type ScalarConstraintD =
  -- | Parameter vector
  Vector Double ->
  -- | (Constraint violation, constraint gradient)
  (Double, Vector Double)

-- | A constraint function which returns a vector @c(x)@ given the
-- parameter vector @x@.  The constraint will enforce that @c(x) == 0@
-- (equality constraint) or @c(x) <= 0@ (inequality constraint).
type VectorConstraint =
  -- | Parameter vector
  Vector Double ->
  -- | Constraint vector size
  Word ->
  -- | Constraint violation vector
  Vector Double

-- | A constraint function which returns @c(x)@ given the parameter
-- vector @x@ along with the Jacobian (first derivative) matrix of
-- @c(x)@ with respect to @x@ at that point.  The constraint will
-- enforce that @c(x) == 0@ (equality constraint) or @c(x) <= 0@
-- (inequality constraint).
type VectorConstraintD =
  -- | Parameter vector
  Vector Double ->
  -- | Constraint vector size
  Word ->
  -- | (Constraint violation vector,
  -- constraint Jacobian)
  (Vector Double, Matrix Double)

-- $nonlinearconstraints
--
-- Note that most NLOPT algorithms do not support nonlinear
-- constraints natively; if you need to enforce nonlinear constraints,
-- you may want to use the 'AugLagAlgorithm' family of solvers, which
-- can add nonlinear constraints to some algorithm that does not
-- support them by a principled modification of the objective
-- function.
--
-- == Example program
--
-- The following interactive session example enforces a scalar
-- constraint on the problem given in the beginning of the module: the
-- parameters must always sum to 1.  The minimizer finds a constrained
-- minimum of 22.5 at @(0.5, 0.5)@.
--
-- >>> import Numeric.LinearAlgebra ( dot, fromList, toList )
-- >>> let objf x = x `dot` x + 22
-- >>> let stop = ObjectiveRelativeTolerance 1e-9 :| []
-- >>>          -- define constraint function:
-- >>> let constraintf x = sum (toList x) - 1.0
-- >>>          -- define constraint object to pass to the algorithm:
-- >>> let constraint = EqualityConstraint (Scalar constraintf) 1e-6
-- >>> let algorithm = COBYLA objf [] [] [constraint] Nothing
-- >>> let problem = LocalProblem 2 stop algorithm
-- >>> let x0 = fromList [5, 10]
-- >>> minimizeLocal problem x0
-- Right (Solution {solutionCost = 22.500000000013028, solutionParams = [0.5000025521533521,0.49999744784664796], solutionResult = FTOL_REACHED})

data Constraint s v
  = -- | A scalar constraint.
    Scalar s
  | -- | A vector constraint.
    Vector Word v
  | -- | A scalar constraint with an attached preconditioning function.
    Preconditioned Preconditioner s

-- | An equality constraint, comprised of both the constraint function
-- (or functions, if a preconditioner is used) along with the desired
-- tolerance.
data EqualityConstraint s v = EqualityConstraint
  { eqConstraintFunctions :: Constraint s v,
    eqConstraintTolerance :: Double
  }

-- | An inequality constraint, comprised of both the constraint
-- function (or functions, if a preconditioner is used) along with the
-- desired tolerance.
data InequalityConstraint s v = InequalityConstraint
  { ineqConstraintFunctions :: Constraint s v,
    ineqConstraintTolerance :: Double
  }

-- | A collection of equality constraints that do not supply
-- constraint derivatives.
type EqualityConstraints =
  [EqualityConstraint ScalarConstraint VectorConstraint]

-- | A collection of inequality constraints that do not supply
-- constraint derivatives.
type InequalityConstraints =
  [InequalityConstraint ScalarConstraint VectorConstraint]

-- | A collection of equality constraints that supply constraint
-- derivatives.
type EqualityConstraintsD = [EqualityConstraint ScalarConstraintD VectorConstraintD]

-- | A collection of inequality constraints that supply constraint
-- derivatives.
type InequalityConstraintsD = [InequalityConstraint ScalarConstraintD VectorConstraintD]

class ApplyConstraint constraint where
  applyConstraint :: N.Opt -> constraint -> IO N.Result

instance ApplyConstraint (EqualityConstraint ScalarConstraint VectorConstraint) where
  applyConstraint opt (EqualityConstraint ty tol) = case ty of
    Scalar s ->
      N.add_equality_constraint opt (wrapScalarFunction s) () tol
    Vector n v ->
      N.add_equality_mconstraint opt n (wrapVectorFunction v n) () tol
    Preconditioned p s ->
      N.add_precond_equality_constraint
        opt
        (wrapScalarFunction s)
        (wrapPreconditionerFunction p)
        ()
        tol

instance ApplyConstraint (InequalityConstraint ScalarConstraint VectorConstraint) where
  applyConstraint opt (InequalityConstraint ty tol) = case ty of
    Scalar s ->
      N.add_inequality_constraint opt (wrapScalarFunction s) () tol
    Vector n v ->
      N.add_inequality_mconstraint opt n (wrapVectorFunction v n) () tol
    Preconditioned p s ->
      N.add_precond_inequality_constraint
        opt
        (wrapScalarFunction s)
        (wrapPreconditionerFunction p)
        ()
        tol

instance ApplyConstraint (EqualityConstraint ScalarConstraintD VectorConstraintD) where
  applyConstraint opt (EqualityConstraint ty tol) = case ty of
    Scalar s ->
      N.add_equality_constraint opt (wrapScalarFunctionD s) () tol
    Vector n v ->
      N.add_equality_mconstraint opt n (wrapVectorFunctionD v n) () tol
    Preconditioned p s ->
      N.add_precond_equality_constraint
        opt
        (wrapScalarFunctionD s)
        (wrapPreconditionerFunction p)
        ()
        tol

instance ApplyConstraint (InequalityConstraint ScalarConstraintD VectorConstraintD) where
  applyConstraint opt (InequalityConstraint ty tol) = case ty of
    Scalar s ->
      N.add_inequality_constraint opt (wrapScalarFunctionD s) () tol
    Vector n v ->
      N.add_inequality_mconstraint opt n (wrapVectorFunctionD v n) () tol
    Preconditioned p s ->
      N.add_precond_inequality_constraint
        opt
        (wrapScalarFunctionD s)
        (wrapPreconditionerFunction p)
        ()
        tol

{- Bounds -}

-- | Bound constraints are specified by vectors of the same dimension
-- as the parameter space.
--
-- == Example program
--
-- The following interactive session example enforces lower bounds on
-- the example from the beginning of the module.  This prevents the
-- optimizer from locating the true minimum at @(0, 0)@; a slightly
-- higher constrained minimum at @(1, 1)@ is found.  Note that the
-- optimizer returns 'N.XTOL_REACHED' rather than 'N.FTOL_REACHED',
-- because the bound constraint is active at the final minimum.
--
-- >>> import Numeric.LinearAlgebra ( dot, fromList )
-- >>> let objf x = x `dot` x + 22                           -- define objective
-- >>> let stop = ObjectiveRelativeTolerance 1e-6 :| []      -- define stopping criterion
-- >>> let lowerbound = LowerBounds $ fromList [1, 1]        -- specify bounds
-- >>> let algorithm = NELDERMEAD objf [lowerbound] Nothing  -- specify algorithm
-- >>> let problem = LocalProblem 2 stop algorithm           -- specify problem
-- >>> let x0 = fromList [5, 10]                             -- specify initial guess
-- >>> minimizeLocal problem x0
-- Right (Solution {solutionCost = 24.0, solutionParams = [1.0,1.0], solutionResult = XTOL_REACHED})
data Bounds
  = -- | Lower bound vector @v@ means we want @x >= v@.
    LowerBounds (Vector Double)
  | -- | Upper bound vector @u@ means we want @x <= u@.
    UpperBounds (Vector Double)
  deriving (Eq, Show, Read)

applyBounds :: N.Opt -> Bounds -> IO N.Result
applyBounds opt (LowerBounds lbvec) = N.set_lower_bounds opt lbvec
applyBounds opt (UpperBounds ubvec) = N.set_upper_bounds opt ubvec

{- Stopping conditions -}

-- | A 'StoppingCondition' tells NLOPT when to stop working on a
-- minimization problem.  When multiple 'StoppingCondition's are
-- provided, the problem will stop when any one condition is met.
data StoppingCondition
  = -- | Stop minimizing when an objective value @J@ less than or equal
    -- to the provided value is found.
    MinimumValue Double
  | -- | Stop minimizing when an optimization step changes the objective
    -- value @J@ by less than the provided tolerance multiplied by @|J|@.
    ObjectiveRelativeTolerance Double
  | -- | Stop minimizing when an optimization step changes the objective
    -- value by less than the provided tolerance.
    ObjectiveAbsoluteTolerance Double
  | -- | Stop when an optimization step changes /every element/ of the
    -- parameter vector @x@ by less than @x@ scaled by the provided
    -- tolerance.
    ParameterRelativeTolerance Double
  | -- | Stop when an optimization step changes /every element/ of the
    -- parameter vector @x@ by less than the corresponding element in
    -- the provided vector of tolerances values.
    ParameterAbsoluteTolerance (Vector Double)
  | -- | Stop when the number of evaluations of the objective function
    -- exceeds the provided count.
    MaximumEvaluations Word
  | -- | Stop when the optimization time exceeds the provided time (in
    -- seconds).  This is not a precise limit.
    MaximumTime Double
  deriving (Eq, Show, Read)

-- $nonempty
--
-- The 'NonEmpty' data type from 'Data.List.NonEmpty' is re-exported
-- here, because it is used to ensure that you always specify at least
-- one stopping condition.

applyStoppingCondition :: N.Opt -> StoppingCondition -> IO N.Result
applyStoppingCondition opt (MinimumValue x) = N.set_stopval opt x
applyStoppingCondition opt (ObjectiveRelativeTolerance x) = N.set_ftol_rel opt x
applyStoppingCondition opt (ObjectiveAbsoluteTolerance x) = N.set_ftol_abs opt x
applyStoppingCondition opt (ParameterRelativeTolerance x) = N.set_xtol_rel opt x
applyStoppingCondition opt (ParameterAbsoluteTolerance v) = N.set_xtol_abs opt v
applyStoppingCondition opt (MaximumEvaluations n) = N.set_maxeval opt n
applyStoppingCondition opt (MaximumTime deltat) = N.set_maxtime opt deltat

{- Random seed control -}

-- | This specifies how to initialize the random number generator for
-- stochastic algorithms.
data RandomSeed
  = -- | Seed the RNG with the provided value.
    SeedValue Word
  | -- | Seed the RNG using the system clock.
    SeedFromTime
  | -- | Don't perform any explicit initialization of the RNG.
    Don'tSeed
  deriving (Eq, Show, Read)

applyRandomSeed :: RandomSeed -> IO ()
applyRandomSeed Don'tSeed = return ()
applyRandomSeed (SeedValue n) = N.srand n
applyRandomSeed SeedFromTime = N.srand_time

{- Random stuff -}

-- | This specifies the population size for algorithms that use a pool
-- of solutions.
newtype Population = Population Word deriving (Eq, Show, Read)

applyPopulation :: N.Opt -> Population -> IO N.Result
applyPopulation opt (Population n) = N.set_population opt n

-- | This specifies the memory size to be used by algorithms like
-- 'LBFGS' which store approximate Hessian or Jacobian matrices.
newtype VectorStorage = VectorStorage Word deriving (Eq, Show, Read)

applyVectorStorage :: N.Opt -> VectorStorage -> IO N.Result
applyVectorStorage opt (VectorStorage n) = N.set_vector_storage opt n

-- | This vector with the same dimension as the parameter vector @x@
-- specifies the initial step for the optimizer to take.  (This
-- applies to local gradient-free algorithms, which cannot use
-- gradients to estimate how big a step to take.)
newtype InitialStep = InitialStep (Vector Double) deriving (Eq, Show, Read)

applyInitialStep :: N.Opt -> InitialStep -> IO N.Result
applyInitialStep opt (InitialStep v) = N.set_initial_step opt v

{- Algorithms -}

data GlobalProblem = GlobalProblem
  { -- | Lower bounds for @x@
    lowerBounds :: Vector Double,
    -- | Upper bounds for @x@
    upperBounds :: Vector Double,
    -- | At least one stopping
    -- condition
    gstop :: NonEmpty StoppingCondition,
    -- | Algorithm specification
    galgorithm :: GlobalAlgorithm
  }

-- | These are the global minimization algorithms provided by NLOPT.  Please see
-- <http://ab-initio.mit.edu/wiki/index.php/NLopt_Algorithms the NLOPT algorithm manual>
-- for more details on how the methods work and how they relate to one another.
--
-- Optional parameters are wrapped in a 'Maybe'; for example, if you
-- see 'Maybe' 'Population', you can simply specify 'Nothing' to use
-- the default behavior.
--
-- N.B.: There may be at most n equality constraints for problems with n dimensions due to restrictions
-- in NLopt. This may be circumvented by combination of equality constraints according to:
-- h1(x)==0 && h2(x)==0 <=> h1(x)^2+h2(x)^2==0.
data GlobalAlgorithm
  = -- | DIviding RECTangles
    DIRECT Objective
  | -- | DIviding RECTangles, locally-biased variant
    DIRECT_L Objective
  | -- | DIviding RECTangles, "slightly randomized"
    DIRECT_L_RAND Objective RandomSeed
  | -- | DIviding RECTangles, unscaled version
    DIRECT_NOSCAL Objective
  | -- | DIviding RECTangles, locally-biased and unscaled
    DIRECT_L_NOSCAL Objective
  | -- | DIviding RECTangles, locally-biased, unscaled and "slightly
    -- randomized"
    DIRECT_L_RAND_NOSCAL Objective RandomSeed
  | -- | DIviding RECTangles, original FORTRAN implementation
    ORIG_DIRECT Objective InequalityConstraints
  | -- | DIviding RECTangles, locally-biased, original FORTRAN
    -- implementation
    ORIG_DIRECT_L Objective InequalityConstraints
  | -- | Stochastic Global Optimization.
    -- __This algorithm is only available if you have linked with @libnlopt_cxx@.__
    STOGO ObjectiveD
  | -- | Stochastic Global Optimization, randomized variant.
    -- __This algorithm is only available if you have linked with @libnlopt_cxx@.__
    STOGO_RAND ObjectiveD RandomSeed
  | -- | Controlled Random Search with Local Mutation
    CRS2_LM Objective RandomSeed (Maybe Population)
  | -- | Improved Stochastic Ranking Evolution Strategy, the population defaults to 20×(n+1) in n dimensions
    ISRES Objective InequalityConstraints EqualityConstraints RandomSeed (Maybe Population)
  | -- | Evolutionary Algorithm
    ESCH Objective
  | -- | Original Multi-Level Single-Linkage
    MLSL Objective LocalProblem (Maybe Population)
  | -- | Multi-Level Single-Linkage with Sobol Low-Discrepancy
    -- Sequence for starting points
    MLSL_LDS Objective LocalProblem (Maybe Population)

algorithmEnumOfGlobal :: GlobalAlgorithm -> N.Algorithm
algorithmEnumOfGlobal (DIRECT _) = N.GN_DIRECT
algorithmEnumOfGlobal (DIRECT_L _) = N.GN_DIRECT_L
algorithmEnumOfGlobal (DIRECT_L_RAND _ _) = N.GN_DIRECT_L_RAND
algorithmEnumOfGlobal (DIRECT_NOSCAL _) = N.GN_DIRECT_NOSCAL
algorithmEnumOfGlobal (DIRECT_L_NOSCAL _) = N.GN_DIRECT_L_NOSCAL
algorithmEnumOfGlobal (DIRECT_L_RAND_NOSCAL _ _) = N.GN_DIRECT_L_RAND_NOSCAL
algorithmEnumOfGlobal (ORIG_DIRECT _ _) = N.GN_ORIG_DIRECT
algorithmEnumOfGlobal (ORIG_DIRECT_L _ _) = N.GN_ORIG_DIRECT_L
algorithmEnumOfGlobal (STOGO _) = N.GD_STOGO
algorithmEnumOfGlobal (STOGO_RAND _ _) = N.GD_STOGO_RAND
algorithmEnumOfGlobal (CRS2_LM _ _ _) = N.GN_CRS2_LM
algorithmEnumOfGlobal (ISRES _ _ _ _ _) = N.GN_ISRES
algorithmEnumOfGlobal (ESCH _) = N.GN_ESCH
algorithmEnumOfGlobal (MLSL _ _ _) = N.G_MLSL
algorithmEnumOfGlobal (MLSL_LDS _ _ _) = N.G_MLSL_LDS

applyGlobalObjective :: N.Opt -> GlobalAlgorithm -> IO ()
applyGlobalObjective opt alg = go alg
  where
    obj = tryTo . applyObjective opt . MinimumObjective
    objD = tryTo . applyObjectiveD opt . MinimumObjective

    go (DIRECT o) = obj o
    go (DIRECT_L o) = obj o
    go (DIRECT_NOSCAL o) = obj o
    go (DIRECT_L_NOSCAL o) = obj o
    go (ESCH o) = obj o
    go (STOGO o) = objD o
    go (DIRECT_L_RAND o _) = obj o
    go (DIRECT_L_RAND_NOSCAL o _) = obj o
    go (ORIG_DIRECT o _) = obj o
    go (ORIG_DIRECT_L o _) = obj o
    go (STOGO_RAND o _) = objD o
    go (CRS2_LM o _ _) = obj o
    go (ISRES o _ _ _ _) = obj o
    go (MLSL o _ _) = obj o
    go (MLSL_LDS o _ _) = obj o

applyGlobalAlgorithm :: N.Opt -> GlobalAlgorithm -> IO ()
applyGlobalAlgorithm opt alg = do
  applyGlobalObjective opt alg
  go alg
  where
    seed = applyRandomSeed
    pop = maybe (return ()) (tryTo . applyPopulation opt)
    ic = traverse_ (tryTo . applyConstraint opt)
    ec = traverse_ (tryTo . applyConstraint opt)

    local lp = setupLocalProblem lp >>= N.set_local_optimizer opt

    go (DIRECT_L_RAND _ s) = seed s
    go (DIRECT_L_RAND_NOSCAL _ s) = seed s
    go (ORIG_DIRECT _ ineq) = ic ineq
    go (ORIG_DIRECT_L _ ineq) = ic ineq
    go (STOGO_RAND _ s) = seed s
    go (CRS2_LM _ s p) = seed s *> pop p
    go (ISRES _ ineq eq s p) = ic ineq *> ec eq *> seed s *> pop p
    go (MLSL _ lp p) = local lp *> pop p
    go (MLSL_LDS _ lp p) = local lp *> pop p
    go _ = return ()

tryTo :: IO N.Result -> IO ()
tryTo act = do
  result <- act
  if (N.isSuccess result)
    then return ()
    else Ex.throw $ NloptException result

data NloptException = NloptException N.Result deriving (Show, Typeable)

instance Exception NloptException

-- | Solve the specified global optimization problem.
--
-- = Example program
--
-- The following interactive session example uses the 'ISRES'
-- algorithm, a stochastic, derivative-free global optimizer, to
-- minimize a trivial function with a minimum of 22.0 at @(0, 0)@.
-- The search is conducted within a box from -10 to 10 in each
-- dimension.
--
-- >>> import Numeric.LinearAlgebra ( dot, fromList )
-- >>> let objf x = x `dot` x + 22                         -- define objective
-- >>> let stop = ObjectiveRelativeTolerance 1e-12 :| []   -- define stopping criterion
-- >>> let algorithm = ISRES objf [] [] (SeedValue 22)     -- specify algorithm
-- >>> let lowerbounds = fromList [-10, -10]               -- specify bounds
-- >>> let upperbounds = fromList [10, 10]                 -- specify bounds
-- >>> let problem = GlobalProblem lowerbounds upperbounds stop algorithm
-- >>> let x0 = fromList [5, 8]                            -- specify initial guess
-- >>> minimizeGlobal problem x0
-- Right (Solution {solutionCost = 22.000000000002807, solutionParams = [-1.660591102367038e-6,2.2407062393213684e-7], solutionResult = FTOL_REACHED})
minimizeGlobal ::
  -- | Problem specification
  GlobalProblem ->
  -- | Initial parameter guess
  Vector Double ->
  -- | Optimization results
  Either N.Result Solution
minimizeGlobal prob x0 =
  unsafePerformIO $ (Right <$> minimizeGlobal' prob x0) `Ex.catch` handler
  where
    handler :: NloptException -> IO (Either N.Result a)
    handler (NloptException retcode) = return $ Left retcode

applyGlobalProblem :: N.Opt -> GlobalProblem -> IO ()
applyGlobalProblem opt (GlobalProblem lb ub stop alg) = do
  tryTo $ applyBounds opt (LowerBounds lb)
  tryTo $ applyBounds opt (UpperBounds ub)
  traverse_ (tryTo . applyStoppingCondition opt) stop
  applyGlobalAlgorithm opt alg

newOpt :: N.Algorithm -> Word -> IO N.Opt
newOpt alg sz = do
  opt' <- N.create alg sz
  case opt' of
    Nothing -> Ex.throw $ NloptException N.FAILURE
    Just opt -> return opt

setupGlobalProblem :: GlobalProblem -> IO N.Opt
setupGlobalProblem gp@(GlobalProblem _ _ _ alg) = do
  opt <- newOpt (algorithmEnumOfGlobal alg) (problemSize gp)
  applyGlobalProblem opt gp
  return opt

solveProblem :: N.Opt -> Vector Double -> IO Solution
solveProblem opt x0 = do
  (N.Output outret outcost outx) <- N.optimize opt x0
  if (N.isSuccess outret)
    then return $ Solution outcost outx outret
    else Ex.throw $ NloptException outret

minimizeGlobal' :: GlobalProblem -> Vector Double -> IO Solution
minimizeGlobal' gp x0 = do
  opt <- setupGlobalProblem gp
  solveProblem opt x0

data LocalProblem = LocalProblem
  { -- | The dimension of the
    -- parameter vector.
    lsize :: Word,
    -- | At least one stopping
    -- condition
    lstop :: NonEmpty StoppingCondition,
    -- | Algorithm specification
    lalgorithm :: LocalAlgorithm
  }

-- | These are the local minimization algorithms provided by NLOPT.  Please see
-- <http://ab-initio.mit.edu/wiki/index.php/NLopt_Algorithms the NLOPT algorithm manual>
-- for more details on how the methods work and how they relate to one
-- another.  Note that some local methods require you provide
-- derivatives (gradients or Jacobians) for your objective function
-- and constraint functions.
--
-- Optional parameters are wrapped in a 'Maybe'; for example, if you
-- see 'Maybe' 'VectorStorage', you can simply specify 'Nothing' to
-- use the default behavior.
data LocalAlgorithm
  = -- | Limited-memory BFGS
    LBFGS_NOCEDAL ObjectiveD (Maybe VectorStorage)
  | -- | Limited-memory BFGS
    LBFGS ObjectiveD (Maybe VectorStorage)
  | -- | Shifted limited-memory variable-metric, rank-2
    VAR2 ObjectiveD (Maybe VectorStorage)
  | -- | Shifted limited-memory variable-metric, rank-1
    VAR1 ObjectiveD (Maybe VectorStorage)
  | -- | Truncated Newton's method
    TNEWTON ObjectiveD (Maybe VectorStorage)
  | -- | Truncated Newton's method with automatic restarting
    TNEWTON_RESTART ObjectiveD (Maybe VectorStorage)
  | -- | Preconditioned truncated Newton's method
    TNEWTON_PRECOND ObjectiveD (Maybe VectorStorage)
  | -- | Preconditioned truncated Newton's method with automatic
    -- restarting
    TNEWTON_PRECOND_RESTART ObjectiveD (Maybe VectorStorage)
  | -- | Method of moving averages
    MMA ObjectiveD InequalityConstraintsD
  | -- | Sequential Least-Squares Quadratic Programming
    SLSQP ObjectiveD [Bounds] InequalityConstraintsD EqualityConstraintsD
  | -- | Conservative Convex Separable Approximation
    CCSAQ ObjectiveD Preconditioner
  | -- | PRincipal AXIS gradient-free local optimization
    PRAXIS Objective [Bounds] (Maybe InitialStep)
  | -- | Constrained Optimization BY Linear Approximations
    COBYLA
      Objective
      [Bounds]
      InequalityConstraints
      EqualityConstraints
      (Maybe InitialStep)
  | -- | Powell's NEWUOA algorithm
    NEWUOA Objective (Maybe InitialStep)
  | -- | Powell's NEWUOA algorithm with bounds by SGJ
    NEWUOA_BOUND Objective [Bounds] (Maybe InitialStep)
  | -- | Nelder-Mead Simplex gradient-free method
    NELDERMEAD Objective [Bounds] (Maybe InitialStep)
  | -- | NLOPT implementation of Rowan's Subplex algorithm
    SBPLX Objective [Bounds] (Maybe InitialStep)
  | -- | Bounded Optimization BY Quadratic Approximations
    BOBYQA Objective [Bounds] (Maybe InitialStep)

algorithmEnumOfLocal :: LocalAlgorithm -> N.Algorithm
algorithmEnumOfLocal (LBFGS_NOCEDAL _ _) = N.LD_LBFGS_NOCEDAL
algorithmEnumOfLocal (LBFGS _ _) = N.LD_LBFGS
algorithmEnumOfLocal (VAR2 _ _) = N.LD_VAR2
algorithmEnumOfLocal (VAR1 _ _) = N.LD_VAR1
algorithmEnumOfLocal (TNEWTON _ _) = N.LD_TNEWTON
algorithmEnumOfLocal (TNEWTON_RESTART _ _) = N.LD_TNEWTON_RESTART
algorithmEnumOfLocal (TNEWTON_PRECOND _ _) = N.LD_TNEWTON_PRECOND
algorithmEnumOfLocal (TNEWTON_PRECOND_RESTART _ _) = N.LD_TNEWTON_PRECOND_RESTART
algorithmEnumOfLocal (MMA _ _) = N.LD_MMA
algorithmEnumOfLocal (SLSQP _ _ _ _) = N.LD_SLSQP
algorithmEnumOfLocal (CCSAQ _ _) = N.LD_CCSAQ
algorithmEnumOfLocal (PRAXIS _ _ _) = N.LN_PRAXIS
algorithmEnumOfLocal (COBYLA _ _ _ _ _) = N.LN_COBYLA
algorithmEnumOfLocal (NEWUOA _ _) = N.LN_NEWUOA
algorithmEnumOfLocal (NEWUOA_BOUND _ _ _) = N.LN_NEWUOA
algorithmEnumOfLocal (NELDERMEAD _ _ _) = N.LN_NELDERMEAD
algorithmEnumOfLocal (SBPLX _ _ _) = N.LN_SBPLX
algorithmEnumOfLocal (BOBYQA _ _ _) = N.LN_BOBYQA

applyLocalObjective :: N.Opt -> LocalAlgorithm -> IO ()
applyLocalObjective opt alg = go alg
  where
    obj = tryTo . applyObjective opt . MinimumObjective
    objD = tryTo . applyObjectiveD opt . MinimumObjective
    precond p = tryTo . applyObjectiveD opt . PreconditionedMinimumObjective p

    go (LBFGS_NOCEDAL o _) = objD o
    go (LBFGS o _) = objD o
    go (VAR2 o _) = objD o
    go (VAR1 o _) = objD o
    go (TNEWTON o _) = objD o
    go (TNEWTON_RESTART o _) = objD o
    go (TNEWTON_PRECOND o _) = objD o
    go (TNEWTON_PRECOND_RESTART o _) = objD o
    go (MMA o _) = objD o
    go (SLSQP o _ _ _) = objD o
    go (CCSAQ o prec) = precond prec o
    go (PRAXIS o _ _) = obj o
    go (COBYLA o _ _ _ _) = obj o
    go (NEWUOA o _) = obj o
    go (NEWUOA_BOUND o _ _) = obj o
    go (NELDERMEAD o _ _) = obj o
    go (SBPLX o _ _) = obj o
    go (BOBYQA o _ _) = obj o

applyLocalAlgorithm :: N.Opt -> LocalAlgorithm -> IO ()
applyLocalAlgorithm opt alg = do
  applyLocalObjective opt alg
  go alg
  where
    ic = traverse_ (tryTo . applyConstraint opt)
    icd = traverse_ (tryTo . applyConstraint opt)
    ec = traverse_ (tryTo . applyConstraint opt)
    ecd = traverse_ (tryTo . applyConstraint opt)
    store = maybe (return ()) (tryTo . applyVectorStorage opt)
    bound = traverse_ (tryTo . applyBounds opt)
    step0 = maybe (return ()) (tryTo . applyInitialStep opt)

    go (LBFGS_NOCEDAL _ vs) = store vs
    go (LBFGS _ vs) = store vs
    go (VAR2 _ vs) = store vs
    go (VAR1 _ vs) = store vs
    go (TNEWTON _ vs) = store vs
    go (TNEWTON_RESTART _ vs) = store vs
    go (TNEWTON_PRECOND _ vs) = store vs
    go (TNEWTON_PRECOND_RESTART _ vs) = store vs
    go (MMA _ ineqd) = icd ineqd
    go (SLSQP _ b ineqd eqd) =
      bound b *> icd ineqd *> ecd eqd
    go (CCSAQ _ _) = return ()
    go (PRAXIS _ b s) = bound b *> step0 s
    go (COBYLA _ b ineq eq s) =
      bound b *> ic ineq *> ec eq *> step0 s
    go (NEWUOA _ s) = step0 s
    go (NEWUOA_BOUND _ b s) = bound b *> step0 s
    go (NELDERMEAD _ b s) = bound b *> step0 s
    go (SBPLX _ b s) = bound b *> step0 s
    go (BOBYQA _ b s) = bound b *> step0 s

applyLocalProblem :: N.Opt -> LocalProblem -> IO ()
applyLocalProblem opt (LocalProblem _ stop alg) = do
  traverse_ (tryTo . applyStoppingCondition opt) stop
  applyLocalAlgorithm opt alg

setupLocalProblem :: LocalProblem -> IO N.Opt
setupLocalProblem lp@(LocalProblem sz _ alg) = do
  opt <- newOpt (algorithmEnumOfLocal alg) sz
  applyLocalProblem opt lp
  return opt

minimizeLocal' :: LocalProblem -> Vector Double -> IO Solution
minimizeLocal' lp x0 = do
  opt <- setupLocalProblem lp
  solveProblem opt x0

-- |
-- == Example program
--
-- The following interactive session example enforces the same scalar
-- constraint as the nonlinear constraint example, but this time it
-- uses the SLSQP solver to find the minimum.
--
-- >>> import Numeric.LinearAlgebra ( dot, fromList, toList, scale )
-- >>> let objf x = (x `dot` x + 22, 2 `scale` x)
-- >>> let stop = ObjectiveRelativeTolerance 1e-9 :| []
-- >>> let constraintf x = (sum (toList x) - 1.0, fromList [1, 1])
-- >>> let constraint = EqualityConstraint (Scalar constraintf) 1e-6
-- >>> let algorithm = SLSQP objf [] [] [constraint]
-- >>> let problem = LocalProblem 2 stop algorithm
-- >>> let x0 = fromList [5, 10]
-- >>> minimizeLocal problem x0
-- Right (Solution {solutionCost = 22.5, solutionParams = [0.4999999999999998,0.5000000000000002], solutionResult = FTOL_REACHED})
minimizeLocal :: LocalProblem -> Vector Double -> Either N.Result Solution
minimizeLocal prob x0 =
  unsafePerformIO $ (Right <$> minimizeLocal' prob x0) `Ex.catch` handler
  where
    handler :: NloptException -> IO (Either N.Result a)
    handler (NloptException retcode) = return $ Left retcode

class ProblemSize c where
  problemSize :: c -> Word

instance ProblemSize LocalProblem where
  problemSize = lsize

instance ProblemSize GlobalProblem where
  problemSize = fromIntegral . HM.size . lowerBounds

instance ProblemSize AugLagProblem where
  problemSize (AugLagProblem _ _ alg) = case alg of
    AUGLAG_LOCAL lp _ _ -> problemSize lp
    AUGLAG_EQ_LOCAL lp -> problemSize lp
    AUGLAG_GLOBAL gp _ _ -> problemSize gp
    AUGLAG_EQ_GLOBAL gp -> problemSize gp

-- | __IMPORTANT NOTE__
--
-- For augmented lagrangian problems, you, the user, are responsible
-- for providing the appropriate type of constraint.  If the
-- subsidiary problem requires an `ObjectiveD`, then you should
-- provide constraint functions with derivatives.  If the subsidiary
-- problem requires an `Objective`, you should provide constraint
-- functions without derivatives.  If you don't do this, you may get a
-- runtime error.
data AugLagProblem = AugLagProblem
  { -- | Possibly empty set of
    -- equality constraints
    alEquality :: EqualityConstraints,
    -- | Possibly empty set of
    -- equality constraints with
    -- derivatives
    alEqualityD :: EqualityConstraintsD,
    -- | Algorithm specification.
    alalgorithm :: AugLagAlgorithm
  }

-- | The Augmented Lagrangian solvers allow you to enforce nonlinear
-- constraints while using local or global algorithms that don't
-- natively support them.  The subsidiary problem is used to do the
-- minimization, but the @AUGLAG@ methods modify the objective to
-- enforce the constraints.  Please see
-- <http://ab-initio.mit.edu/wiki/index.php/NLopt_Algorithms the NLOPT algorithm manual>
-- for more details on how the methods work and how they relate to one another.
--
-- See the documentation for 'AugLagProblem' for an important note
-- about the constraint functions.
data AugLagAlgorithm
  = -- | AUGmented LAGrangian with a local subsidiary method
    AUGLAG_LOCAL LocalProblem InequalityConstraints InequalityConstraintsD
  | -- | AUGmented LAGrangian with a local subsidiary method and with
    -- penalty functions only for equality constraints
    AUGLAG_EQ_LOCAL LocalProblem
  | -- | AUGmented LAGrangian with a global subsidiary method
    AUGLAG_GLOBAL GlobalProblem InequalityConstraints InequalityConstraintsD
  | -- | AUGmented LAGrangian with a global subsidiary method and with
    -- penalty functions only for equality constraints.
    AUGLAG_EQ_GLOBAL GlobalProblem

algorithmEnumOfAugLag :: AugLagAlgorithm -> N.Algorithm
algorithmEnumOfAugLag (AUGLAG_LOCAL _ _ _) = N.AUGLAG
algorithmEnumOfAugLag (AUGLAG_EQ_LOCAL _) = N.AUGLAG_EQ
algorithmEnumOfAugLag (AUGLAG_GLOBAL _ _ _) = N.AUGLAG
algorithmEnumOfAugLag (AUGLAG_EQ_GLOBAL _) = N.AUGLAG_EQ

-- | This structure is returned in the event of a successful
-- optimization.
data Solution = Solution
  { -- | The objective function value
    -- at the minimum
    solutionCost :: Double,
    -- | The parameter vector which
    -- minimizes the objective
    solutionParams :: Vector Double,
    -- | Why the optimizer stopped
    solutionResult :: N.Result
  }
  deriving (Eq, Show, Read)

applyAugLagAlgorithm :: N.Opt -> AugLagAlgorithm -> IO ()
applyAugLagAlgorithm opt alg = go alg
  where
    ic = traverse_ (tryTo . applyConstraint opt)
    icd = traverse_ (tryTo . applyConstraint opt)
    -- AUGLAG won't work at all if you don't pass it the same
    -- objective as the subproblem -- here we pull out the subproblem
    -- objectives from the algorithm spec and set the same objective
    -- function so the user can't mess it up.
    local lp = tryTo $ do
      localopt <- setupLocalProblem lp
      applyLocalObjective opt (lalgorithm lp)
      N.set_local_optimizer opt localopt
    global gp = do
      tryTo $ setupGlobalProblem gp >>= N.set_local_optimizer opt
      applyGlobalObjective opt (galgorithm gp)

    go (AUGLAG_LOCAL lp ineq ineqd) = local lp *> ic ineq *> icd ineqd
    go (AUGLAG_EQ_LOCAL lp) = local lp
    go (AUGLAG_GLOBAL gp ineq ineqd) = global gp *> ic ineq *> icd ineqd
    go (AUGLAG_EQ_GLOBAL gp) = global gp

applyAugLagProblem :: N.Opt -> AugLagProblem -> IO ()
applyAugLagProblem opt (AugLagProblem eq eqd alg) = do
  traverse_ (tryTo . applyConstraint opt) eq
  traverse_ (tryTo . applyConstraint opt) eqd
  applyAugLagAlgorithm opt alg

minimizeAugLag' :: AugLagProblem -> Vector Double -> IO Solution
minimizeAugLag' ap@(AugLagProblem _ _ alg) x0 = do
  opt <- newOpt (algorithmEnumOfAugLag alg) (problemSize ap)
  applyAugLagProblem opt ap
  solveProblem opt x0

-- |
-- == Example program
--
-- The following interactive session example enforces the same scalar
-- constraint as the nonlinear constraint example, but this time it
-- uses the augmented Lagrangian method to enforce the constraint and
-- the 'SBPLX' algorithm, which does not support nonlinear constraints
-- itself, to perform the minimization.  As before, the parameters
-- must always sum to 1, and the minimizer finds the same constrained
-- minimum of 22.5 at @(0.5, 0.5)@.
-- >>> import Numeric.LinearAlgebra ( dot, fromList, toList )
-- >>> let objf x = x `dot` x + 22
-- >>> let stop = ObjectiveRelativeTolerance 1e-9 :| []
-- >>> let algorithm = SBPLX objf [] Nothing
-- >>> let subproblem = LocalProblem 2 stop algorithm
-- >>> let x0 = fromList [5, 10]
-- >>> minimizeLocal subproblem x0
-- Right (Solution {solutionCost = 22.0, solutionParams = [0.0,0.0], solutionResult = FTOL_REACHED})
-- >>>          -- define constraint function:
-- >>> let constraintf x = sum (toList x) - 1.0
-- >>>          -- define constraint object to pass to the algorithm:
-- >>> let constraint = EqualityConstraint (Scalar constraintf) 1e-6
-- >>> let problem = AugLagProblem [constraint] [] (AUGLAG_EQ_LOCAL subproblem)
-- >>> minimizeAugLag problem x0
-- Right (Solution {solutionCost = 22.500000015505844, solutionParams = [0.5000880506776678,0.4999119493223323], solutionResult = FTOL_REACHED})
minimizeAugLag :: AugLagProblem -> Vector Double -> Either N.Result Solution
minimizeAugLag prob x0 =
  unsafePerformIO $ (Right <$> minimizeAugLag' prob x0) `Ex.catch` handler
  where
    handler :: NloptException -> IO (Either N.Result a)
    handler (NloptException retcode) = return $ Left retcode


text:
{-# language OverloadedStrings #-}
module Report (generateReports) where

import Graphics.Vega.VegaLite hiding (Sum)
import Prelude hiding (filter, lookup, repeat)
import Data.Monoid (Sum(..))
import Numeric.LinearAlgebra ((#>), Vector, toList)
import Data.Maybe

import GA
import Fitness

-- | Plot the evolution of average fitness along the generations
plotEvo :: [Double] -> IO ()
plotEvo avgs =
  let gens = [1.0 .. fromIntegral (length avgs)]
      plotData = dataFromColumns []
               . dataColumn "generation" (Numbers gens)
               . dataColumn "avg" (Numbers avgs)
               $ []
      enc      = encoding
               . position X [ PName "generation", PmType Quantitative]
               . position Y [ PName "avg", PmType Quantitative]
               $ []
  in  toHtmlFile "evolution.html" $ toVegaLite [ plotData, mark Line [], enc, height 800, width 600 ]

plotPoly :: [[Double]] -> Vector Double -> Solution (Poly n) -> IO ()
plotPoly xss ys sol =
  let x0 = map head xss
      zss = decode xss (map _getPoly $ _chromo sol)
      ysHat = toList $ zss #> (fromJust . _coeffs) sol
      ys'   = toList ys
      clusters = replicate (length x0) "poly" ++ replicate (length x0) "real"
      plotData = dataFromColumns []
               . dataColumn "x0" (Numbers $ x0 ++ x0)
               . dataColumn "y" (Numbers $ ysHat ++ ys')
               . dataColumn "Cluster" (Strings clusters)
               $ []
      enc      = encoding
               . position X [ PName "x0", PmType Quantitative]
               . position Y [ PName "y", PmType Quantitative]
               . color [ MName "Cluster", MmType Nominal ]
               $ []
  in  toHtmlFile "poly.html" $ toVegaLite [ plotData, mark Line [], enc, height 800, width 600 ]

generateReports :: [[Double]] -> Vector Double -> [Double] -> Solution (Poly n) -> IO ()
generateReports xss ys avgs best = do plotEvo avgs
                                      plotPoly xss ys best

ext: hs

max_stars_repo_name: folivetti/gapoly

max_issues_repo_name: folivetti/gapoly

max_forks_repo_name: folivetti/gapoly


text:
{-|
Module      : SmithDecomposition
Description : Descomposición de matrices para la resoluciones de sistemas de ecuaciones diofánticas. Específico a la forma del sistema de ecuaciones asociado a la matriz de viajes.

Este módulo no provee funciones generales para encontrar la Forma Normal de Smith, sino que genera los resultados basado en el patron regular del sistema.
-}
module ODMatrix.SmithDecomposition (
    computeODM
  --
  , systemMatrix
  , snfR
  --
  , vectorToODM
  , odmToVector
  --
  , tVectorSize
  -- Publish from Dimensions
  , srows
  , scols
  , nFromN
  , nFromM
  --
  , generalSolution
  ) where

  import Prelude hiding ((<>))

  import Numeric.LinearAlgebra as L

  import ODMatrix (ODM, FeatVec, toFlatList, fromFlatList)
  import Control.Monad (join)

  import ODMatrix.SmithDecomposition.SparseMat (SparseMat(SparseMat), denseSparseMat, makeAssoc)
  import ODMatrix.SmithDecomposition.Dimensions

  -- * Solution computation

  -- | Given a input vector and a t vector, computes a unique solution of the system.
  computeODM :: FeatVec      -- ^ Input Vector
             -> FeatVec      -- ^ t Vector
             -> ODM          -- ^ Solution Matrix
  computeODM b t = vectorToODM $ r1b + r2 <> asColumn t
    where (r1b, r2) = generalSolution b

  -- | Given and input vector, computes the general solution of the system.
  -- This is a partial solution. The aplicaction of a t vector is needed to find a particular one.
  generalSolution :: FeatVec      -- ^ Input Vector
                  -> (ODM, ODM)   -- ^ Solution Matrix
  generalSolution b = (r1 <> asColumn b, r2)
    where (r1, r2) = splitR $ snfR n
          n = nFromN $ size b

  -- * Smith's Decomposition

  -- | Given the size of the route, computes the system matrix.
  systemMatrix :: Int -> Matrix Double
  systemMatrix n = denseSparseMat $ SparseMat r c ab
    where r = srows n
          c = scols n
          a = zip (join $ map (\x -> [x..n]) [1..n]) [1..c]
          b = zip (join $ map (\x -> replicate (n-x) (n+x)) [1..n-1]) [n+1..c]
          ab = zip (a++b) $ repeat 1.0

  -- | Given the size of the route, computes the R matrix.
  snfR :: Int -> Matrix Double
  snfR n = stepA n <> stepB n <> stepC n

  stepA :: Int -> Matrix Double
  stepA n = accum (ident m) (+) as
    where as = makeAssoc (-1.0) $ zip sr [n+1..]
          m = scols n
          sr = join [[k..n] | k <- [2..n]]

  stepB :: Int -> Matrix Double
  stepB n = accum (ident m) (+) $ as
    where m = scols n
          as = makeAssoc (-1.0) $ sbs
          sbs = join . map (\(p,k) -> zip [k,k..] [k+1..k+p]) . zip [n-2,n-3..1] $ ser n

  stepC :: Int -> Matrix Double
  stepC n = accum (ident m) (+) $ as ++ bs
    where m = scols n
          as = makeAssoc (-1.0) . join . map (\(a,b) -> [(a,a),(b,b)]) $ swp
          bs = makeAssoc ( 1.0) . join . map (\(a,b) -> [(a,b),(b,a)]) $ swp
          swp = tail . zip [n+1..] $ ser n

  ser :: Int -> [Int]
  ser n = foldl (\(a:ac) x -> (a-x):a:ac) [scols n] [2..n-1]

  -- * Computation of the t-vector

  -- | The size of the t vector given the size of the route.
  tVectorSize :: Int -> Int
  tVectorSize n = scols n - srows n

  -- * Support Functions

  -- | Break the right matrix R, for the computation of the solution
  splitR :: Matrix Double -> (Matrix Double, Matrix Double)
  splitR r = (takeColumns n r, dropColumns n r)
    where n = srows . nFromM $ cols r

  -- | Transform a vector in column matrix form to a ODM
  vectorToODM :: Matrix Double -> ODM
  vectorToODM = fromFlatList . toList . flatten

  -- | Transform a ODM to vector in column matrix form
  odmToVector :: ODM -> Matrix Double
  odmToVector = matrix 1 . toFlatList

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

  -- | Indices fila y columna de la matriz cuadrada triangular superior con la diagonal en cero, a partir del índice aplanado.
  --
  -- fromFlatIndex :: Int       -- ^ Tamaño de la matriz
  --               -> Int       -- ^ Índice plano
  --               -> (Int,Int) -- ^ (fila, columna)
  -- fromFlatIndex n k = as !! (k-1)
  --   where as = [1..n] >>= (\x -> [(x,j) | j <- [x+1..n]])

  -- | Indice aplanado de la matriz cuadrada triangular superior con la diagonal en cero.
  --
  -- toFlatIndex :: Int        -- ^ Tamaño de la matriz
  --             -> (Int, Int) -- ^ (fila, columna)
  --             -> Int        -- ^ Índice plano
  -- toFlatIndex' n (r,c)
  --   | r >= c = 0
  --   | otherwise = fr r - n + c
  --   where fr 0 = 0
  --         fr r = fr (r-1) + n - r

  -- toFlatIndex n (r,c)
  --   | r >= c = 0
  --   | otherwise = fr r - n + c
  --   where fr r = (r * n) - (r * (r + 1) `div` 2)

  -- * Construcción de matriz del sistema

  -- upperCond :: Int -> Int -> Int -> R
  -- upperCond n k c
  --   | (i <= k) && (k < j) = 1
  --   | otherwise           = 0
  --   where (i,j) = fromFlatIndex n c

  -- lowerCond :: Int -> Int -> Int -> R
  -- lowerCond n k c
  --   | (i == k+1) && (k+1 < j) = 1
  --   | otherwise               = 0
  --   where (i,j) = fromFlatIndex n c

  -- partialSystemMatrix :: (Int -> Int -> Int -> R) -> Int -> Matrix R
  -- partialSystemMatrix c n = (rs><cs) h
  --   where rs = n - 1
  --         cs = scols n
  --         l = [1..(scols n)]
  --         h = [ c n i j | i <- [1..rs], j <- l ]

  -- systemMatrix :: Int -> Matrix Double
  -- systemMatrix n = um === lm
  --   where um = partialSystemMatrix upperCond n
  --         lm = takeRows (n-2) $ partialSystemMatrix lowerCond n
  --
  --

  -- * System calculation

  -- | Fuck yeah!
  -- calculatorMachine :: Vector R             -- ^ Input vector
  --                   -> (Matrix R, Matrix R) -- ^ Two part calculator
  -- calculatorMachine b = (sub1, sub2)
  --   where n = nFromInputVector b
  --         (d,l,r) = snf n
  --         b' = asColumn b
  --         (r1,r2) = splitr n r
  --         sub1 = r1 <> l <> b'
  --         sub2 = r2

  -- calcODMatrix :: InputVector  -- ^ Input Vector
  --              -> TVector      -- ^ t Vector
  --              -> ODM          -- ^ Solution Vector
  -- calcODMatrix b t = s1 + s2 <> asColumn t
  --   where (s1, s2) = calculatorMachine b

  -- | Validate the solution of the system Ax = b and that all values of the solution are positive integers.
  -- solutionValidation :: Vector R  -- ^ Input Vector
  --                    -> Vector R  -- ^ Solution Vector
  --                    -> Bool      -- ^ True if validate.
  -- solutionValidation b x = p1 && p2
  --   where p1 = all (>= 0) $ toList x
  --         n = nFromInputVector b
  --         a = systemMatrix n
  --         p2 = a <> asColumn x == asColumn b

  -- tVectors :: Int         -- ^ N
  --          -> Int         -- ^ K
  --          -> [Vector R]  -- ^ All possible t vectors. Not necessarilly for a valid solution.
  -- tVectors n k = tVectors' n k v0
  --   where ss = scols n - srows n
  --         v0 = vector $ replicate ss 0

  -- tVectors' :: Int         -- ^ N
  --           -> Int         -- ^ K
  --           -> Vector R    -- ^ Initial vector for search.
  --           -> [Vector R]  -- ^ All possible t vectors. Not necessarilly for a valid solution.
  -- tVectors' n k v0 = vector <$> list
  --   where list = iterate (nextc k) $ toList v0
  --         nextc k ls =
  --           let (t:ts) = dropWhile (== fromIntegral k) ls
  --               s = length ls - length (t:ts)
  --               h = replicate s 0
  --           in h ++ [t+1] ++ ts

  -- tVector :: Vector R  -- ^ Input Vector
  --         -> Int       -- ^ K
  --         -> Vector R  -- ^ t Vector
  -- tVector b k = tVector' b k v0
  --   where n = srows' b
  --         v0 = vector $ replicate (tSize n) 0

  -- tVector' :: Vector R  -- ^ Input Vector
  --          -> Int       -- ^ K
  --          -> Vector R  -- ^ Initial vector for search
  --          -> Vector R  -- ^ t Vector
  -- tVector' b k v0 = head $ dropWhile cond $ tVectors' n k v0
  --   where n = nFromInputVector b
  --         cond t = validate k . solutionVectorToODMatrix . flatten $ calcODMatrix b t

  -- solutionVectorToODMatrix :: Vector R -> ODM
  -- solutionVectorToODMatrix = denseSparseMat . solutionVectorToSparseMat

  -- * System decomposition

  -- | Smith Normal Form
  -- Dado el tamaño N de una matriz de viajes A, retorna una tripla (D,L,R) donde D está en la Forma Normal de Smith y LAR = D.
  -- snf :: Int
  --     -> (Matrix R, Matrix R, Matrix R)
  -- snf n = (snfD n, snfL n, snfR n)
  -- --
  -- --
  -- snfD :: Int
  --      -> Matrix R
  -- snfD n = diagRect 0 (fromList $ replicate (srows n) 1) (srows n) (scols n)
  -- --
  -- --
  -- snfL :: Int
  --      -> Matrix R
  -- snfL = toDense . cells . sparseSnfL

  -- -- | Componente derecho de la descomposición de Smith.
  -- snfR :: Int
  --      -> Matrix R
  -- snfR = toDense . cells . sparseSnfR

max_stars_repo_name: renecura/odmatrix

max_issues_repo_name: renecura/odmatrix

max_forks_repo_name: renecura/odmatrix


text:
module Main where

import Numeric.LinearAlgebra.Data (loadMatrix)

import qualified Euler.Problem001
import qualified Euler.Problem002
import qualified Euler.Problem003
import qualified Euler.Problem004
import qualified Euler.Problem005
import qualified Euler.Problem006
import qualified Euler.Problem007
import qualified Euler.Problem008
import qualified Euler.Problem009
import qualified Euler.Problem010
import qualified Euler.Problem011
import qualified Euler.Problem012
import qualified Euler.Problem013
import qualified Euler.Problem014
import qualified Euler.Problem015
import qualified Euler.Problem016
import qualified Euler.Problem017
import qualified Euler.Problem018
import qualified Euler.Problem019
import qualified Euler.Problem020
import qualified Euler.Problem021
import qualified Euler.Problem022
import qualified Euler.Problem023
import qualified Euler.Problem024
import qualified Euler.Problem025
import qualified Euler.Problem026
import qualified Euler.Problem027
import qualified Euler.Problem028
import qualified Euler.Problem029
import qualified Euler.Problem030
import qualified Euler.Problem031

main :: IO ()
main = do
  label "001" $ Euler.Problem001.solution 1000
  label "002" $ Euler.Problem002.solution 4000000
  label "003" $ Euler.Problem003.solution 600851475143
  label "004" $ Euler.Problem004.solution 1000
  label "005" $ Euler.Problem005.solution 20
  label "006" $ Euler.Problem006.solution 100
  label "007" $ Euler.Problem007.solution 10001
  label "008" $ Euler.Problem008.solution 13
  label "009" $ Euler.Problem009.solution 1000
  label "010" $ Euler.Problem010.solution 2000000
  matrix011 <- loadMatrix "input/011.txt"
  label "011" $ Euler.Problem011.solution 4 matrix011
  label "012" $ Euler.Problem012.solution 500
  label "013" $ Euler.Problem013.solution 10
  label "014" $ Euler.Problem014.solution 1000000
  label "015" $ Euler.Problem015.solution 20
  label "016" $ Euler.Problem016.solution 1000
  label "017" $ Euler.Problem017.solution 1000
  file018 <- readFile "input/018.txt"
  label "018" $ Euler.Problem018.solution file018
  label "019" $ Euler.Problem019.solution (1901, 1, 1) (2000, 12, 31)
  label "020" $ Euler.Problem020.solution 100
  label "021" $ Euler.Problem021.solution 10000
  file022 <- readFile "input/022.txt"
  label "022" $ Euler.Problem022.solution file022
  label "023" $ Euler.Problem023.solution 28123
  label "024" $ Euler.Problem024.solution 1000000
  label "025" $ Euler.Problem025.solution 1000
  label "026" $ Euler.Problem026.solution (1000 :: Int)
  label "027" $ Euler.Problem027.solution (1000 :: Int)
  label "028" $ Euler.Problem028.solution (1001 :: Int)
  label "029" $ Euler.Problem029.solution (100 :: Int)
  label "030" $ Euler.Problem030.solution (5 :: Int)
  label "031" $ Euler.Problem031.solution (200 :: Int)

label :: Integral n => String -> n -> IO ()
label s i = do
  putStr $ "Problem " ++ s ++ " solution: "
  print $ toInteger i

max_stars_repo_name: whittle/euler

max_issues_repo_name: whittle/euler

max_forks_repo_name: whittle/euler


text:
{-# LANGUAGE DataKinds                                #-}
{-# LANGUAGE DeriveDataTypeable                       #-}
{-# LANGUAGE DeriveGeneric                            #-}
{-# LANGUAGE FlexibleInstances                        #-}
{-# LANGUAGE GADTs                                    #-}
{-# LANGUAGE KindSignatures                           #-}
{-# LANGUAGE MultiParamTypeClasses                    #-}
{-# LANGUAGE PatternSynonyms                          #-}
{-# LANGUAGE RankNTypes                               #-}
{-# LANGUAGE RecordWildCards                          #-}
{-# LANGUAGE ScopedTypeVariables                      #-}
{-# LANGUAGE TypeApplications                         #-}
{-# LANGUAGE TypeFamilies                             #-}
{-# LANGUAGE TypeInType                               #-}
{-# LANGUAGE TypeOperators                            #-}
{-# LANGUAGE UndecidableInstances                     #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.Normalise       #-}

module Backprop.Learn.Model.Neural.LSTM (
    -- * Basic LSTM
    LSTM, pattern LSTM
  , LSTMp(..), lstmForget, lstmInput, lstmUpdate, lstmOutput
  , LSTM'(..)
    -- * GRU
  , GRU, pattern GRU
  , GRUp(..), gruMemory, gruUpdate, gruOutput
  , GRU'(..)
  ) where

import           Backprop.Learn.Initialize
import           Backprop.Learn.Model.Class
import           Backprop.Learn.Model.Function
import           Backprop.Learn.Model.Neural
import           Backprop.Learn.Model.Regression
import           Backprop.Learn.Model.State
import           Control.DeepSeq
import           Data.Typeable
import           GHC.Generics                          (Generic)
import           GHC.TypeNats
import           Lens.Micro
import           Numeric.Backprop
import           Numeric.LinearAlgebra.Static.Backprop
import           Numeric.OneLiner
import           Numeric.Opto.Ref
import           Numeric.Opto.Update
import qualified Data.Binary                           as Bi
import qualified Numeric.LinearAlgebra.Static          as H

-- | "Base" for 'LSTM'.  An 'LSTM' is just an 'LSTM'' where the second half
-- of the input vector is the previous output.
--
-- This is mostly an implementation detail.  It is recommended that you use
-- 'LSTM' and its constructor, instead of directly using this.
data LSTM' (i :: Nat) (o :: Nat) = LSTM'
  deriving (Typeable)

-- TODO: allow parameterize internal activation function?
-- TODO: Peepholes

-- | 'LSTM' layer parmateters
data LSTMp (i :: Nat) (o :: Nat) =
    LSTMp { _lstmForget :: !(FCp (i + o) o)
          , _lstmInput  :: !(FCp (i + o) o)
          , _lstmUpdate :: !(FCp (i + o) o)
          , _lstmOutput :: !(FCp (i + o) o)
          }
  deriving (Generic, Typeable, Show)

instance NFData (LSTMp i o)
instance (KnownNat i, KnownNat o) => Additive (LSTMp i o)
instance (KnownNat i, KnownNat o) => Scaling Double (LSTMp i o)
instance (KnownNat i, KnownNat o) => Metric Double (LSTMp i o)
instance (KnownNat i, KnownNat o, Ref m (LSTMp i o) v) => AdditiveInPlace m v (LSTMp i o)
instance (KnownNat i, KnownNat o, Ref m (LSTMp i o) v) => ScalingInPlace m v Double (LSTMp i o)
instance (KnownNat i, KnownNat o) => Bi.Binary (LSTMp i o)
instance (KnownNat i, KnownNat o) => Backprop (LSTMp i o)

-- | Forget biases initialized to 1
instance (KnownNat i, KnownNat o) => Initialize (LSTMp i o) where
    initialize d g = LSTMp <$> set (mapped . fcBias) 1 (initialize d g)
                           <*> initialize d g
                           <*> initialize d g
                           <*> initialize d g

instance (KnownNat i, KnownNat o) => Num (LSTMp i o) where
    (+)         = gPlus
    (-)         = gMinus
    (*)         = gTimes
    negate      = gNegate
    abs         = gAbs
    signum      = gSignum
    fromInteger = gFromInteger

instance (KnownNat i, KnownNat o) => Fractional (LSTMp i o) where
    (/)          = gDivide
    recip        = gRecip
    fromRational = gFromRational

instance (KnownNat i, KnownNat o) => Floating (LSTMp i o) where
    pi    = gPi
    sqrt  = gSqrt
    exp   = gExp
    log   = gLog
    sin   = gSin
    cos   = gCos
    asin  = gAsin
    acos  = gAcos
    atan  = gAtan
    sinh  = gSinh
    cosh  = gCosh
    asinh = gAsinh
    acosh = gAcosh
    atanh = gAtanh

lstmForget :: Lens' (LSTMp i o) (FCp (i + o) o)
lstmForget f x = (\y -> x { _lstmForget = y }) <$> f (_lstmForget x)

lstmInput :: Lens' (LSTMp i o) (FCp (i + o) o)
lstmInput f x = (\y -> x { _lstmInput = y }) <$> f (_lstmInput x)

lstmUpdate :: Lens' (LSTMp i o) (FCp (i + o) o)
lstmUpdate f x = (\y -> x { _lstmUpdate = y }) <$> f (_lstmUpdate x)

lstmOutput :: Lens' (LSTMp i o) (FCp (i + o) o)
lstmOutput f x = (\y -> x { _lstmOutput = y }) <$> f (_lstmOutput x)

instance (KnownNat i, KnownNat o, KnownNat io, io ~ (i + o)) => Learn (R io) (R o) (LSTM' i o) where
    type LParamMaybe (LSTM' i o) = 'Just (LSTMp i o)
    type LStateMaybe (LSTM' i o) = 'Just (R o)

    runLearn LSTM' (J_ p) x (J_ s) = (h, J_ s')
      where
        forget = logistic $ runLRp (p ^^. lstmForget) x
        input  = logistic $ runLRp (p ^^. lstmInput ) x
        update = tanh     $ runLRp (p ^^. lstmUpdate) x
        s'     = forget * s + input * update
        o      = logistic $ runLRp (p ^^. lstmOutput) x
        h      = o * tanh s'

-- | Long-term short-term memory layer
--
-- <http://colah.github.io/posts/2015-08-Understanding-LSTMs/>
--
-- @
-- instance 'Learn' ('R' i) (R o) ('LSTM' i o) where
--     type 'LParamMaybe' (LSTM i o) = ''Just' ('LSTMp' i o)
--     type 'LStateMaybe' (LSTM i o) = 'Just ('T2' (R o) (R o))     -- cell state, hidden state
-- @
type LSTM i o = Recurrent (R (i + o)) (R i) (R o) (R o) (LSTM' i o)

-- | Construct an 'LSTM'
pattern LSTM :: (KnownNat i, KnownNat o) => LSTM i o
pattern LSTM <- Rec { _recLearn = LSTM' }
  where
    LSTM = Rec
      { _recSplit = H.split
      , _recJoin  = (H.#)
      , _recLoop  = id
      , _recLearn = LSTM'
      }

-- | "Base" for 'GRU'.  An 'GRU' is just an 'GRU'' where the second half
-- of the input vector is the previous output.
--
-- This is mostly an implementation detail.  It is recommended that you use
-- 'GRU' and its constructor, instead of directly using this.
data GRU' (i :: Nat) (o :: Nat) = GRU'
  deriving (Typeable)

-- | 'GRU' layer parmateters
data GRUp (i :: Nat) (o :: Nat) =
    GRUp { _gruMemory :: !(FCp (i + o) o)
         , _gruUpdate :: !(FCp (i + o) o)
         , _gruOutput :: !(FCp (i + o) o)
         }
  deriving (Generic, Typeable, Show)

instance NFData (GRUp i o)
instance (KnownNat i, KnownNat o) => Additive (GRUp i o)
instance (KnownNat i, KnownNat o) => Scaling Double (GRUp i o)
instance (KnownNat i, KnownNat o) => Metric Double (GRUp i o)
instance (KnownNat i, KnownNat o, Ref m (GRUp i o) v) => AdditiveInPlace m v (GRUp i o)
instance (KnownNat i, KnownNat o, Ref m (GRUp i o) v) => ScalingInPlace m v Double (GRUp i o)
instance (KnownNat i, KnownNat o) => Bi.Binary (GRUp i o)
instance (KnownNat i, KnownNat o) => Backprop (GRUp i o)

instance (KnownNat i, KnownNat o) => Initialize (GRUp i o) where
    initialize d g = GRUp <$> initialize d g
                          <*> initialize d g
                          <*> initialize d g

instance (KnownNat i, KnownNat o) => Num (GRUp i o) where
    (+)         = gPlus
    (-)         = gMinus
    (*)         = gTimes
    negate      = gNegate
    abs         = gAbs
    signum      = gSignum
    fromInteger = gFromInteger

instance (KnownNat i, KnownNat o) => Fractional (GRUp i o) where
    (/)          = gDivide
    recip        = gRecip
    fromRational = gFromRational

instance (KnownNat i, KnownNat o) => Floating (GRUp i o) where
    pi    = gPi
    sqrt  = gSqrt
    exp   = gExp
    log   = gLog
    sin   = gSin
    cos   = gCos
    asin  = gAsin
    acos  = gAcos
    atan  = gAtan
    sinh  = gSinh
    cosh  = gCosh
    asinh = gAsinh
    acosh = gAcosh
    atanh = gAtanh

gruMemory :: Lens' (GRUp i o) (FCp (i + o) o)
gruMemory f x = (\y -> x { _gruMemory = y }) <$> f (_gruMemory x)

gruUpdate :: Lens' (GRUp i o) (FCp (i + o) o)
gruUpdate f x = (\y -> x { _gruUpdate = y }) <$> f (_gruUpdate x)

gruOutput :: Lens' (GRUp i o) (FCp (i + o) o)
gruOutput f x = (\y -> x { _gruOutput = y }) <$> f (_gruOutput x)

instance (KnownNat i, KnownNat o, KnownNat io, io ~ (i + o)) => Learn (R io) (R o) (GRU' i o) where
    type LParamMaybe (GRU' i o) = 'Just (GRUp i o)
    type LStateMaybe (GRU' i o) = 'Nothing

    runLearn GRU' (J_ p) = stateless $ \x ->
        let z      = logistic $ runLRp (p ^^. gruMemory) x
            r      = logistic $ runLRp (p ^^. gruUpdate) x
            r'     = 1 # r
            h'     = tanh     $ runLRp (p ^^. gruOutput) (r' * x)
        in  (1 - z) * snd (split @i x) + z * h'

-- | Gated Recurrent Unit
--
-- <http://colah.github.io/posts/2015-08-Understanding-LSTMs/>
--
-- @
-- instance 'Learn' ('R' i) (R o) ('GRU' i o) where
--     type 'LParamMaybe' (GRU i o) = ''Just' ('GRUp' i o)
--     type 'LStateMaybe' (GRU i o) = 'Just (R o)
-- @
type GRU i o = Recurrent (R (i + o)) (R i) (R o) (R o) (GRU' i o)

-- | Construct an 'GRU'
pattern GRU :: (KnownNat i, KnownNat o) => GRU i o
pattern GRU <- Rec { _recLearn = GRU' }
  where
    GRU = Rec
      { _recSplit = H.split
      , _recJoin  = (H.#)
      , _recLoop  = id
      , _recLearn = GRU'
      }

max_stars_repo_name: mstksg/backprop-learn

max_issues_repo_name: mstksg/backprop-learn

max_forks_repo_name: mstksg/backprop-learn


text:
{-# LANGUAGE FunctionalDependencies, FlexibleInstances, FlexibleContexts, BangPatterns #-}

module AI.Learning.LinearRegression where

import Data.List (foldl')
import Numeric.LinearAlgebra
import Numeric.LinearAlgebra.Util (ones)

import AI.Util.Matrix

{-

A linear model should contain the following information:

* Regression coefficients
* Covariance matrix of regression coefficients
* Residuals
* Confidence intervals for regression coefficients
* Confidence intervals for residuals
* T-Statistic for regression coefficients
* P-Values for regression coefficients
* F-Statistic for regression
* P-Value for regression
* MSE
* SST
* SSE
* PRESS Statistic
* R squared
* adjusted r squared
*

-}

-------------
-- Options --
-------------

-- |Type recording what kind of linear model we want to build. Choose from
--  ordinary least squares (OLS), ridge regression with parameter lambda
--  and LASSO with parameter lambda.
data LMType = OLS
            | Ridge Double
            | LASSO Double
            deriving (Show)

---- |Data type for linear regression options.
data LMOpts = LMOpts { fitIntercept :: Bool
                     , standardizeRegressors :: Bool }
                     deriving (Show)

---- |Standard options for linear regression (used by 'regress').
stdLMOpts = LMOpts { fitIntercept = True
                   , standardizeRegressors = False }

---------------------
---- Linear Models --
---------------------

-- |Data type for a linear model. It consists of the coefficient vector,
--  together with options that specify how to transform any data that is used
--  to make new predictions.
data LinearModel = LM { coefs   :: Vector Double
                      , lmIntercept :: Bool
                      , lmStandardize :: Bool
                      , lmMean :: Vector Double
                      , lmStd  :: Vector Double }
                      deriving (Show)

-- |Statistics structure for a linear regression.
data LMStats = LMStats { covBeta :: Maybe (Matrix Double)
                       , betaCI :: Maybe (Matrix Double)
                       , sst :: Double
                       , sse :: Double
                       , mse :: Double
                       , rSquare :: Double
                       , tBeta :: Maybe (Vector Double)
                       , pBeta :: Maybe (Vector Double)
                       , fRegression :: Maybe Double
                       , pRegression :: Maybe Double }
                       deriving (Show)

type Predictor = Matrix Double
type Response = Vector Double

-- |Fit an ordinary least squares linear model to data.
lm :: Matrix Double -> Vector Double -> LinearModel
lm = lmWith OLS stdLMOpts

-- |Fit a ridge regression (least squares with quadratic penalty) to data.
lmRidge :: Matrix Double -> Vector Double -> Double -> LinearModel
lmRidge x y lambda = lmWith (Ridge lambda) stdLMOpts x y

-- |Fit a linear model to data. The exact model fitted is specified by the
--  'LMType' argument, and the options are specified in 'LMOpts'.
lmWith :: LMType -> LMOpts -> Matrix Double -> Vector Double -> LinearModel
lmWith kind opts x y = LM { coefs = beta
                          , lmIntercept = fitIntercept opts
                          , lmStandardize = standardizeRegressors opts
                          , lmMean = mu
                          , lmStd = sigma }
    where
        (xx,mu,sigma) = lmPrepare opts x
        beta = case kind of
            OLS     -> regress xx y
            Ridge a -> ridgeRegress xx y (fitIntercept opts) a
            LASSO a -> error "LASSO not implemented."

-- |Prepare data according to the specified options structure. This may involve
--  centering and standardizing the data, or adding a column of constants.
lmPrepare :: LMOpts
          -> Matrix Double
          -> (Matrix Double, Vector Double, Vector Double)
lmPrepare opts x = (x3,mu,sigma)
    where
        (x1,mu,sigma) = standardize x
        x2            = if standardizeRegressors opts then x1 else x
        x3            = if fitIntercept opts then addOnes x2 else x2

---------------------------------
---- Predict from linear model --
---------------------------------

-- |Make predictions from a linear model.
lmPredict :: LinearModel -> Matrix Double -> Vector Double
lmPredict model x = x2 <> beta
    where
        beta    = coefs model
        xbar    = lmMean model
        sigma   = lmStd model
        x1      = if lmStandardize model
                    then eachRow (\x -> (x - xbar) / sigma) x
                    else x
        x2      = if lmIntercept model
                    then addOnes x1
                    else x1

-- |Calculate statistics for a linear regression.
lmStats :: LinearModel -> Matrix Double -> Vector Double -> LMStats
lmStats model x y =
    let ybar      = constant (mean y) (dim y)
        yhat      = lmPredict model x
        residuals = y - yhat
        covBeta   = Nothing
        betaCI    = Nothing
        sst       = sumVector $ (y - ybar) ^ 2
        sse       = sumVector $ residuals ^ 2
        mse       = sse / fromIntegral (rows x)
        rSq       = 1 - sse / sst
        tBeta     = Nothing
        pBeta     = Nothing
        fReg      = Nothing
        pReg      = Nothing
    in LMStats covBeta betaCI sst sse mse rSq tBeta pBeta fReg pReg

mseEvalFun :: Response -> Response -> Double
mseEvalFun actual predicted = mean $ (actual - predicted) ^ 2

lmPredFun :: Predictor -> Response -> Predictor -> Response -> Double
lmPredFun xtrain ytrain xtest ytest = mseEvalFun ytest ypred
    where
        ypred = lmPredict model xtest
        model = lm xtrain ytrain

--------------------------
---- Perform Regression --
--------------------------

---- |Regress a vector y against a matrix of predictors x, with the specified
----  options.
regress :: Matrix Double -- X
        -> Vector Double -- y
        -> Vector Double -- beta
regress x y
    | rows x /= dim y = error "Inconsistent dimensions -- REGRESS"
    | otherwise = let (_,n) = size x
                      (_,r) = qr x
                      rr = takeRows n r
                   in (trans rr <> rr) <\> trans x <> y

-- |Ridge regression. This adds a penalty term to OLS regression, which
--  discourages large coefficient values to prevent overfitting.
ridgeRegress :: Matrix Double  -- X
             -> Vector Double  -- y
             -> Bool           -- useConst?
             -> Double         -- lambda
             -> Vector Double  -- beta
ridgeRegress x y useConst lambda
    | rows x /= dim y = error "Inconsistent dimensions -- RIDGEREGRESS"
    | otherwise = let (_,n) = size x
                      (_,r) = qr x
                      rr = takeRows n r
                      ww = if useConst
                            then diag $ join [0, constant 1 (n-1)]
                            else ident n
                   in (trans rr <> rr + lambda `scale` ww) <\> trans x <> y

-----------------
---- Utilities --
-----------------

-- |De-mean a sample.
demean :: Matrix Double -> Matrix Double
demean x = eachRow (subtract $ mean x) x

-- |Standardize a sample to have zero mean and unit variance. Returns a triple
--  consisting of the standardized sample, the sample mean and the sample
--  standard deviation.
standardize :: Matrix Double -> (Matrix Double, Vector Double, Vector Double)
standardize m = (eachRow (\x -> (x - mu) / sigma) m, mu, sigma)
    where mu    = mean m
          sigma = std m

-- |Standardize a sample to have zero mean and unit variance. This variant of
--  the function discards the mean and standard deviation vectors, only
--  returning the standardized sample.
standardize_ :: Matrix Double -> Matrix Double
standardize_ x = a where (a,_,_) = standardize x

-------------------------
-- Mean, Variance etc. --
-------------------------

class Mean a b | a -> b where
    mean :: a -> b

class Floating b => Variance a b | a -> b where
    var :: a -> b
    var x = y * y where y = std x

    std :: a -> b
    std x = sqrt $ var x

instance Mean (Vector Double) Double where
    mean v = sumVector v / fromIntegral (dim v)

instance Variance (Vector Double) Double where
    var v  = mean $ (v - constant vbar (dim v)) ^ 2
        where vbar = mean v

instance Mean (Matrix Double) (Vector Double) where
    mean m = fromList $ mapCols mean m

instance Variance (Matrix Double) (Vector Double) where
    var m  = fromList $ mapCols var m

max_stars_repo_name: cagix/aima-haskell


text:
{-# LANGUAGE ScopedTypeVariables, TypeFamilies, TupleSections, FlexibleContexts #-}

module Fuml.Unsupervised where

import Fuml.Core
import qualified Fuml.Base.PCA as PCA
import Fuml.Base.KNN (euclideanDistance)
import Numeric.LinearAlgebra
import qualified Data.Vector.Storable as VS
import qualified Data.Vector.Unboxed as VU
import qualified Data.Vector as V
import Data.List (foldl1', sortBy, nub)
import Data.Ord (comparing)
import Control.Monad.Identity

type Unsupervisor m p a = Supervisor m () p a

with :: b -> [a] -> [(a,b)]
with y = map (,y)

data PCA = PCA { pcaComponents :: Int, pcaStat :: PCA.Stat }

pca :: Int -> Unsupervisor Identity PCA.Stat (Vector Double)
pca ncomp  = Supervisor $ \_ theData ->
  let stat@(m,_,v) = PCA.statVs $ map fst theData
  in return $ Predict stat $ \x -> takeRows ncomp v #> (x - m)

findNearestCentroidIx :: [Vector Double] -> Vector Double -> Int
findNearestCentroidIx ctrs x =
  let vdixs = zip [0..] $ map (\c -> euclideanDistance c x) ctrs
  in fst $ head $ sortBy (comparing snd) vdixs

centroid :: [Vector Double] -> Vector Double
centroid vs =
  let vadd = VS.zipWith (+)
      n = realToFrac $ length vs
  in VS.map (/n) $ foldl1' vadd vs

{-kmeans :: Int -> Unsupervisor Identity [Vector Double] Int
kmeans nclus
  =  Supervisor $ \_ theData ->
        let clus = KM.kmeans (VU.convert) KM.euclidSq nclus (map fst theData)
            ctrs = map (centroid . KM.elements) $ V.toList clus
            pr v = findNearestCentroidIx ctrs v
        in return $ Predict ctrs pr

cluster :: Eq b => [a] -> (a-> Vector Double) -> Unsupervisor Identity p b -> [[a]]
cluster xs f unsup =
  let p = runIdentity $ runSupervisor unsup Nothing $ with () $ map (f) xs
      withClus x = (x,predict p $ f x)
      withCluss = map withClus xs
      cluss = nub $ map snd withCluss
      getElems clus = map fst $ filter ((==clus) . snd) withCluss
  in map getElems cluss -}

max_stars_repo_name: ekalosak/open

max_issues_repo_name: ekalosak/open

max_forks_repo_name: ekalosak/open


text:
{-# OPTIONS_GHC -Wall #-}

-- | A magma heirarchy for multiplication. The basic magma structure is repeated and prefixed with 'Multiplicative-'.
module Plankton.Multiplicative
  ( MultiplicativeMagma(..)
  , MultiplicativeUnital(..)
  , MultiplicativeAssociative
  , MultiplicativeCommutative
  , MultiplicativeInvertible(..)
  , product
  , Multiplicative(..)
  , MultiplicativeRightCancellative(..)
  , MultiplicativeLeftCancellative(..)
  , MultiplicativeGroup(..)
  ) where

import Data.Complex (Complex(..))
import Plankton.Additive
import qualified Protolude as P
import Protolude (Bool(..), Double, Float, Int, Integer)

-- | 'times' is used as the operator for the multiplicative magam to distinguish from '*' which, by convention, implies commutativity
--
-- > ∀ a,b ∈ A: a `times` b ∈ A
--
-- law is true by construction in Haskell
class MultiplicativeMagma a where
  times :: a -> a -> a

instance MultiplicativeMagma Double where
  times = (P.*)

instance MultiplicativeMagma Float where
  times = (P.*)

instance MultiplicativeMagma Int where
  times = (P.*)

instance MultiplicativeMagma Integer where
  times = (P.*)

instance MultiplicativeMagma Bool where
  times = (P.&&)

instance (MultiplicativeMagma a, AdditiveGroup a) =>
         MultiplicativeMagma (Complex a) where
  (rx :+ ix) `times` (ry :+ iy) =
    (rx `times` ry - ix `times` iy) :+ (ix `times` ry + iy `times` rx)

-- | Unital magma for multiplication.
--
-- > one `times` a == a
-- > a `times` one == a
class MultiplicativeMagma a =>
      MultiplicativeUnital a where
  one :: a

instance MultiplicativeUnital Double where
  one = 1

instance MultiplicativeUnital Float where
  one = 1

instance MultiplicativeUnital Int where
  one = 1

instance MultiplicativeUnital Integer where
  one = 1

instance MultiplicativeUnital Bool where
  one = True

instance (AdditiveUnital a, AdditiveGroup a, MultiplicativeUnital a) =>
         MultiplicativeUnital (Complex a) where
  one = one :+ zero

-- | Associative magma for multiplication.
--
-- > (a `times` b) `times` c == a `times` (b `times` c)
class MultiplicativeMagma a =>
      MultiplicativeAssociative a

instance MultiplicativeAssociative Double

instance MultiplicativeAssociative Float

instance MultiplicativeAssociative Int

instance MultiplicativeAssociative Integer

instance MultiplicativeAssociative Bool

instance (AdditiveGroup a, MultiplicativeAssociative a) =>
         MultiplicativeAssociative (Complex a)

-- | Commutative magma for multiplication.
--
-- > a `times` b == b `times` a
class MultiplicativeMagma a =>
      MultiplicativeCommutative a

instance MultiplicativeCommutative Double

instance MultiplicativeCommutative Float

instance MultiplicativeCommutative Int

instance MultiplicativeCommutative Integer

instance MultiplicativeCommutative Bool

instance (AdditiveGroup a, MultiplicativeCommutative a) =>
         MultiplicativeCommutative (Complex a)

-- | Invertible magma for multiplication.
--
-- > ∀ a ∈ A: recip a ∈ A
--
-- law is true by construction in Haskell
class MultiplicativeMagma a =>
      MultiplicativeInvertible a where
  recip :: a -> a

instance MultiplicativeInvertible Double where
  recip = P.recip

instance MultiplicativeInvertible Float where
  recip = P.recip

instance (AdditiveGroup a, MultiplicativeInvertible a) =>
         MultiplicativeInvertible (Complex a) where
  recip (rx :+ ix) = (rx `times` d) :+ (negate ix `times` d)
    where
      d = recip ((rx `times` rx) `plus` (ix `times` ix))

-- | Idempotent magma for multiplication.
--
-- > a `times` a == a
class MultiplicativeMagma a =>
      MultiplicativeIdempotent a

instance MultiplicativeIdempotent Bool

-- | product definition avoiding a clash with the Product monoid in base
--
product :: (Multiplicative a, P.Foldable f) => f a -> a
product = P.foldr (*) one

-- | Multiplicative is commutative, associative and unital under multiplication
--
-- > one * a == a
-- > a * one == a
-- > (a * b) * c == a * (b * c)
-- > a * b == b * a
class ( MultiplicativeCommutative a
      , MultiplicativeUnital a
      , MultiplicativeAssociative a
      ) =>
      Multiplicative a where
  infixl 7 *
  (*) :: a -> a -> a
  a * b = times a b

instance Multiplicative Double

instance Multiplicative Float

instance Multiplicative Int

instance Multiplicative Integer

instance Multiplicative Bool

instance (AdditiveGroup a, Multiplicative a) => Multiplicative (Complex a)

-- | Non-commutative left divide
--
-- > recip a `times` a = one
class ( MultiplicativeUnital a
      , MultiplicativeAssociative a
      , MultiplicativeInvertible a
      ) =>
      MultiplicativeLeftCancellative a where
  infixl 7 ~/
  (~/) :: a -> a -> a
  a ~/ b = recip b `times` a

-- | Non-commutative right divide
--
-- > a `times` recip a = one
class ( MultiplicativeUnital a
      , MultiplicativeAssociative a
      , MultiplicativeInvertible a
      ) =>
      MultiplicativeRightCancellative a where
  infixl 7 /~
  (/~) :: a -> a -> a
  a /~ b = a `times` recip b

-- | Divide ('/') is reserved for where both the left and right cancellative laws hold.  This then implies that the MultiplicativeGroup is also Abelian.
--
-- > a / a = one
-- > recip a = one / a
-- > recip a * a = one
-- > a * recip a = one
class (Multiplicative a, MultiplicativeInvertible a) =>
      MultiplicativeGroup a where
  infixl 7 /
  (/) :: a -> a -> a
  (/) a b = a `times` recip b

instance MultiplicativeGroup Double

instance MultiplicativeGroup Float

instance (AdditiveGroup a, MultiplicativeGroup a) =>
         MultiplicativeGroup (Complex a)

max_stars_repo_name: chessai/plankton

max_issues_repo_name: chessai/plankton

max_forks_repo_name: chessai/plankton


text:
{-# LANGUAGE ViewPatterns ,DeriveGeneric,DeriveAnyClass ,BangPatterns #-}
module Report where
import           Criterion.IO
import           Criterion.Types
import           Data.Bifunctor
import           Data.List
import qualified Data.Map                        as M
import           Data.Ord
import           Statistics.Resampling.Bootstrap
import           System.Directory
import           Text.Printf
import System.FilePath
import Control.DeepSeq
import GHC.Generics
import Control.Monad

type Measures = M.Map String Measure

data Measure = Measure {mTest::String,mValue::Double} deriving (Show,Read,Eq,Generic,NFData)

mKindPkg :: Measure -> (String, String)
mKindPkg m = let (k,_:p) = break (== '-') $ mTest m in (k,p)

mKind :: Measure -> String
mKind = fst . mKindPkg

mKindPkgVal :: Measure -> (String, (String, Double))
mKindPkgVal m = let (k,p) = mKindPkg m in (k,(p,mValue m))

byTestKind :: M.Map k Measure -> [(String, (String, Double))]
byTestKind = map mKindPkgVal . M.elems

toMeasures :: [Report] -> Measures
toMeasures = M.fromList . map (\r-> let m = toMeasure r in (mTest m,m))

toMeasure :: Report -> Measure
toMeasure r = Measure (reportName r) ((1000 *) . estPoint . anMean . reportAnalysis $ r)

deleteMeasures :: FilePath -> IO ()
deleteMeasures = removeFile . measuresFile

updateMeasures :: FilePath -> IO ()
updateMeasures = void . updateMeasures_

updateMeasures_ :: FilePath -> IO (Measures, Measures, Measures)
updateMeasures_ dir = do
  m' <- readCriterionMeasures dir
  addMeasures_ dir m'

addMeasures_
  :: FilePath
     -> Measures
     -> IO (Measures, Measures, Measures)
addMeasures_ dir m' = do
  !m <- readMeasures dir
  let m'' = M.union m' m
  writeMeasures dir m''
  return (m,m',m'')

addMeasures
  :: FilePath
     -> String
     -> [(String, Double)]
     -> IO () -- Measures, Measures, Measures)
addMeasures dir name ms = void $ addMeasures_ dir (M.fromList $ map (\(pkg,val) -> let n = concat [name,"-",pkg] in (n,Measure n val)) ms)

readMeasures :: FilePath -> IO Measures
readMeasures dir =  do
  let f = measuresFile dir
  fe <- doesFileExist f
  if not fe
    then return M.empty
    else (force . read) <$> readFile f

writeMeasures :: FilePath -> Measures -> IO ()
writeMeasures dir = writeFile (measuresFile dir) . show

readCriterionMeasures :: FilePath -> IO Measures
readCriterionMeasures dir = toMeasures <$> readReports (reportsFile dir)

measuresFile :: FilePath -> FilePath
measuresFile dir = dir </> "measures"

reportsFile :: FilePath -> FilePath
reportsFile dir = dir </> "report.json"

printMeasures :: FilePath -> IO ()
printMeasures dir = readMeasures dir >>= printMeasuresAll_ (const True)

printMeasuresCurrent :: (Measures,Measures,Measures) -> IO ()
printMeasuresCurrent (_,m',m'') =
  let currentKinds = allKinds m'
  in printMeasuresAll_ (\ts -> fst ts `elem` currentKinds) m''

printMeasuresAll_
  :: ((String, [(String, Double)]) -> Bool)
     -> Measures -> IO ()
printMeasuresAll_ f m'' = mapM_ (uncurry report) . filter f . allTests . byTestKind $ m''

allKinds :: Measures -> [String]
allKinds = nub . map mKind . M.elems

allTests
  :: [(String, (String, Double))] -> [(String, [(String, Double)])]
allTests = sort . map (\g -> (fst . head $ g,map snd g)) . groupBy (\a b -> fst a == fst b)

printMeasuresDiff :: (Measures,Measures,Measures) -> IO ()
printMeasuresDiff (m,m',_) =
--  mapM_ (\(n,d) -> putStrLn (unwords [n,show d,"%"])) . M.toList $ M.intersectionWithKey (\k a b ->(k,round ((mValue b/mValue a-1)*100))) m m'
  mapM_ (\(n,d) -> putStrLn (unwords [n,show d,"%"])) . M.toList $ M.intersectionWith (\a b -> round ((mValue b/mValue a-1)*100)::Int) m m'

readReports :: FilePath -> IO [Report]
readReports jsonReportFile = do
  fe <- doesFileExist jsonReportFile
  if not fe
    then return [] -- M.empty
    else do
      er <- readJSONReports jsonReportFile
      return $
        case er of
          Right (_, _, reports) -> reports
          _ -> [] -- M.empty

report :: String -> [(String,Double)] -> IO ()
report _ [] = return ()
report name rs = do
  -- print rs
  let (_,rss) = report_ rs
  let width = maximum . map (length . fst) $ rs
  -- putStrLn $ unwords [name,"ordered by",prop,"("++fst best++":",printInt (snd best),unit++")"]
  putStrLn $ unwords [name,"(best first)"] -- package: "++fst best++" with ",printInt (snd best)++")"]
  -- putStrLn name -- ,"(Best package: "++fst best++" with ",printInt (snd best)++")"]
  mapM_ (\(n,v) -> putStrLn $ unwords [printString width n,printDouble v]) rss
  putStrLn ""

report_ :: (Fractional c, Ord c) => [(a, c)] -> ((a, c), [(a, c)])
report_ rs =
  let
    rss = sortBy (comparing snd) rs
    best = snd . head $ rss
  in (head rss, map (second (/ best)) rss)

printDouble :: Double -> String
printDouble = printf "%5.1f"

printInt :: Double -> String
printInt = printf "%.0f"

printString :: Int -> String -> String
printString width = printf ("%-"++show width++"s")

max_stars_repo_name: nickkuk/flat

max_issues_repo_name: nickkuk/flat

max_forks_repo_name: nickkuk/flat


text:
{-# LANGUAGE TypeOperators #-}
module Vect
    ( module Vect
    ) where

import Control.Arrow (first, second)
import Control.Monad (ap)
import Data.Complex
import Data.List (intercalate, sortOn)
import GHC.TypeLits

import QBasis

type C = Complex Double

-- |Free vector space over field
newtype Vect k b = V [(b, k)] deriving (Eq, Ord)

instance (Show k, Show b) => Show (Vect k b) where
    show (V xs) = intercalate " + " $ map (\(e, c) -> show c ++ show e) xs

instance Functor (Vect k) where
    fmap f (V xs) = V $ map (first f) xs

instance Num k => Applicative (Vect k) where
    pure = return
    (<*>) = ap

instance Num k => Monad (Vect k) where
    return a = V [(a, 1)]
    V ts >>= f = V $ concat [ [(b, y*x) | let V us = f a, (b, y) <- us] | (a, x) <- ts]

-- |Convert an element of Vect k b into normal form. Normal form consists in having the basis elements in ascending order,
-- with no duplicates, and all coefficients non-zero
nf :: (Ord b, Eq k, Num k) => Vect k b -> Vect k b
nf (V ts) = V $ nf' $ sortOn fst ts where
    nf' ((b1,x1):(b2,x2):ts) =
        case compare b1 b2 of
        LT -> if x1 == 0 then nf' ((b2,x2):ts) else (b1,x1) : nf' ((b2,x2):ts)
        EQ -> if x1+x2 == 0 then nf' ts else nf' ((b1,x1+x2):ts)
        GT -> error "nf': not pre-sorted"
    nf' [(b,x)] = if x == 0 then [] else [(b,x)]
    nf' [] = []

-- |The zero vector
zero :: Vect k b
zero = V []

-- |Vector addition
add :: (Ord b, Num k, Eq k) => Vect k b -> Vect k b -> Vect k b
add (V ts) (V us) = V $ addmerge ts us
    where addmerge axs@((a,x):xs) bxs@((b,y):ys) =
            case compare a b of
                LT -> (a,x) : addmerge xs bxs
                EQ -> if x+y == 0
                        then addmerge xs ys
                        else (a, x+y) : addmerge xs ys
                GT -> (b,y) : addmerge axs ys
          addmerge xs [] = xs
          addmerge [] ys = ys

-- |(Left) Scalar multiplication
smult :: (Eq k, Num k) => k -> Vect k b -> Vect k b
smult 0 _ = zero
smult k (V xs) = V $ map (second (k*)) xs

-- |Linear application to basis elements
linear :: (Ord b, Eq k, Num k) => (a -> Vect k b) -> Vect k a -> Vect k b
linear f v = nf $ v >>= f

-- |Tensor product of two free vector spaces
te :: (Num k, Eq k) => Vect k (QBasis m) -> Vect k (QBasis n) -> Vect k (QBasis (m+n))
te (V u) (V v) = V [(a `tb` b, x*y) | (a,x) <- u, (b,y) <- v]

-- |Tensor product of two functions between free vector spaces
tf :: (Num k, Eq k, KnownNat m) =>
            (Vect k (QBasis m) -> Vect k (QBasis m)) ->
            (Vect k (QBasis n) -> Vect k (QBasis n)) ->
            Vect k (QBasis (m+n))-> Vect k (QBasis (m+n))
tf f g (V ts) = sum [let (a,b) = bt basis in smult x (f (return a) `te` g (return b)) | (basis,x) <- ts]
        where sum = foldl add zero

max_stars_repo_name: caldwellwg/Aqua

max_issues_repo_name: caldwellwg/Aqua

max_forks_repo_name: caldwellwg/Aqua


text:
module ML.NN.Tests (tests) where

import ML.NN (ActivationFunction, Layer(..), Network(..), computeActivations,
              feedForward, runNetwork)

import Numeric.LinearAlgebra ((><), konst, size, sumElements, vector)

import Test.HUnit (Assertion, (@?=))
import Test.Tasty (TestTree, testGroup)
import Test.Tasty.HUnit (testCase)

identityActivation :: ActivationFunction
identityActivation x = x

emptyNetwork :: Network
emptyNetwork = Network []

tests :: TestTree
tests = testGroup "ML.NN"
    [
        testCase
        "FeedForward_SumInputWeights_NoBias" testFeedForward_SumInputWeights_NoBias
    ,   testCase
        "RunNetwork_Empty" testRunNetwork_Empty
    ,   testCase
        "ComputeActivations_EmptyNetwork" testComputeActivations_EmptyNetwork
    ,   testCase
        "ComputeActivations_SingleLayerSingleNeuron" testComputeActivations_SingleLayerSingleNeuron
    ]

-- | Test a simple layer where each neuron sums the input.
testFeedForward_SumInputWeights_NoBias :: Assertion
testFeedForward_SumInputWeights_NoBias =
    feedForward identityActivation x sumLayer @?= expected
    where x          = vector [1,2,3]
          expected   = vector $ replicate numNeurons (sumElements x)
          numNeurons = 5
          sumLayer   = Layer (konst 0 numNeurons) (konst 1 (numNeurons, size x))

-- | Test that running an empty network on an input 'x' yields 'x'.
testRunNetwork_Empty :: Assertion
testRunNetwork_Empty =
    runNetwork emptyNetwork identityActivation x @?= x
    where x = vector [1,2,3]

-- | When there are no layers, there should be no z's and the only
--   activation should be the input.
testComputeZsAndAs_EmptyNetwork :: Assertion
testComputeZsAndAs_EmptyNetwork =
    computeZsAndAs identityActivation emptyNetwork x @?= ([], [x])
    where x = vector [1,2,3]

-- | When there is a single identity layer, the zs = [z] where z is the sum
--   of the input values, and the activations are z and the input x.
--
--   The order of the activations is reversed from the order they appear in the
--   network so they can easily be used by the backpropogation algorithm.
testComputeZsAndAs_SingleLayerSingleNeuron :: Assertion
testComputeZsAndAs_SingleLayerSingleNeuron =
    computeZsAndAs identityActivation net x @?= ([z], [z,x])
    where x = vector [1,2,3]
          z = vector [6]
          net = Network [Layer (vector [0]) ((1><3) [1,1,1])]

max_stars_repo_name: m-renaud/ML

max_issues_repo_name: m-renaud/ML

max_forks_repo_name: m-renaud/ML


text:
module Lib.Practice.P3_1 (
    f3_1
) where

import Lib.Practice.P3 (Input( Fire, NotFire ), input2Int, inputs)
import Numeric.LinearAlgebra

f3_1 = showOutput ("XOR", gXor)

gAnd = Single (SNeuron (vector [0.5, 0.5]) (-0.7))
gNand = Single (SNeuron (vector [-0.5, -0.5]) 0.7)
gOr = Single (SNeuron (vector [0.5, 0.5]) (-0.2))
gXor = Multi MNeuron { inputLayer = [gNand, gOr], intermediateLayer = gAnd }

showOutput (name, gate) = do
    putStrLn $ name ++ "ゲート"
    mapM_ (showOutput' gate) inputs
    putStrLn ""
    where
        showOutput' gate input = do
            let result = input --> gate
            putStrLn $ show input ++ " -> " ++ show result

data SNeuron = SNeuron {
    weight :: Vector R,
    bias :: R
} deriving Show

data MNeuron = MNeuron {
    inputLayer :: [Neuron],
    intermediateLayer :: Neuron
} deriving Show

data Neuron = Single SNeuron | Multi MNeuron deriving Show
(-->) :: [Input] -> Neuron -> Input
inputs --> Single SNeuron { weight = w, bias = b } = result
    where inputV    = fromList $ map input2Int inputs
          sum'      = sumElements $ inputV * w
          result'   = b + sum'
          result    = case () of
            _ | result' >= 0 -> Fire
            _                -> NotFire
inputs --> Multi MNeuron { inputLayer = inpLs, intermediateLayer = intL } = result
    where inputs'   = map (inputs -->) inpLs
          result    = inputs' --> intL


text:
module Main where

import Model
import InteractionNetwork
import Config

import Numeric.LinearAlgebra

import System.Environment

getInteractionNetwork :: [Model] -> InteractionNetwork
getInteractionNetwork ts@(m:_) = interactionNetwork (numAgents m) $ map getInteractions ts

getInteractionNetworkAfter :: Int -> [Model] -> InteractionNetwork
getInteractionNetworkAfter k ts = getInteractionNetwork $ drop k ts

main :: IO ()
main = do
  [configFile,dataDir] <- getArgs
  config <- loadConfig configFile

  let model = newModel config
      inet = getInteractionNetworkAfter 1000 (runModel 1.0 1500 model)
      temporalBST = allPairsBFS inet
  putStrLn $ "TCC:  " ++ (show $ tcc inet)
  putStrLn $ "CTPL: " ++ (show $ ctpl inet temporalBST)
  putStrLn $ "TGE:  " ++ (show $ tge inet temporalBST)
  putStrLn $ "1/max ρ(·): " ++ (show $ 1 / (spectralRadius inet))
  putStrLn $ "1/min ρ(·): " ++ (show $ 1 / (spectralRadius' inet))

  -- let cm  = fromColumns . (\(a,b) -> a:b:[]) $ centralities inet
  --     ocm = fromColumns . (\(a,b) -> a:b:[]) $ oneHopCentralities inet

  -- saveMatrix (dataDir ++ "centralities_"++motionType++"-"++n++"-"++arena++fileNameExtra) "%f" cm
  -- saveMatrix (dataDir ++ "oneHopCentralities_"++motionType++"-"++n++"-"++arena++fileNameExtra) "%f" ocm

max_stars_repo_name: wfvining/tvg-smallworld

max_issues_repo_name: wfvining/tvg-smallworld

max_forks_repo_name: wfvining/tvg-smallworld


text:
module Square where

import Genetics
import qualified Neural
import Utilities
import qualified Numeric.LinearAlgebra
    ( cmap )
import Numeric.LinearAlgebra
    ( (<>)
    , (><)
    )
import System.Random

newtype Phenotype = Network Double

trainingSet :: TrainingSet
trainingSet = [ ((1><3)[1, 1, 1], (6><1)[1, 1, 1, 0, 1, 1])
              , ((1><3)[0, 0, 1], (6><1)[0, 0, 0, 0, 0, 1])
              , ((1><3)[0, 1, 0], (6><1)[0, 0, 0, 1, 0, 0])
              , ((1><3)[0, 1, 1], (6><1)[0, 0, 1, 0, 0, 1])
              , ((1><3)[1, 0, 0], (6><1)[0, 1, 0, 0, 0, 0])
              ]

randompopulation :: IO [Genotype]
randompopulation = do
    g <- getStdGen
    return $ map Neural.flatten $ take 10 $ randomNetworks g [3, 6]

morphogenesis :: Genotype -> Phenotype
morphogenesis g = Neural.fromList [3, 6]

fit :: Genotype -> Double
fit = (fitness trainingSet) . morphogenesis

-- Iterates fitness and selection on a population, where the
-- selection function uses randomness.
iterations :: RandomGen g => g -> Population -> [Population]
iterations g population =
    map fst $ iterate (mySelection fit) (population, g)

maximums :: [Population] -> [Double]
maximums populations = map maximum $ (map.map) fit $ populations

-- example to see how well we are doing
example :: IO [Double]
example = do
    g <- newStdGen
    population <- randompopulation
    return $ maximums (iterations g population)

fitness :: Neural.TrainingSet -> Neural.Network -> Double
fitness training network =
    let xs = map fst training
        ys = concat (map (Neural.toList . Neural.flatten . snd) training)
        actual = concat (map (Neural.toList . Neural.flatten . (Neural.forward network)) xs)
        differences = zipWith (-) ys actual
    in (-1) * sqrt (sum (map (^2) differences))

max_stars_repo_name: dbeecham/NeuralFourInARow

max_issues_repo_name: dbeecham/NeuralFourInARow

max_forks_repo_name: dbeecham/NeuralFourInARow