diff --git a/apps/egraphGP/Main.hs b/apps/egraphGP/Main.hs --- a/apps/egraphGP/Main.hs +++ b/apps/egraphGP/Main.hs @@ -140,8 +140,8 @@ evaluateUnevaluated runEqSat myCost rewriteBasic2 1 - while (numberOfEvalClasses nEvals) 1 $ - \radius -> + while (( do --nEvs <- gets (FingerTree.size . _fitRangeDB . _eDB) nCls <- gets (IM.size . _eClass) @@ -179,7 +179,9 @@ do runEqSat myCost rewriteBasic2 1 cleanDB pure () - if b then pure (min 20 $ radius+1) else pure (max 1 $ radius-1) + let radius' = if b then (min 20 $ radius+1) else (max 1 $ radius-1) + nEvs' = nEvs + if upd then 1 else 0 + pure (radius', nEvs') eclasses <- gets (IntMap.toList . _eClass) -- forM_ eclasses $ \(_, v) -> (io.print) (Set.size (_eNodes v), Set.size (_parents v)) paretoFront @@ -271,7 +273,8 @@ insertRndExpr :: Int -> RndEGraph EClassId insertRndExpr maxSize = do grow <- rnd toss - t <- rnd $ Random.randomTree 2 8 maxSize rndTerm rndNonTerm2 grow + n <- rnd (randomFrom [3 .. maxSize]) + t <- rnd $ Random.randomTree 2 8 n rndTerm rndNonTerm2 grow fromTree myCost t >>= canonical insertBestExpr :: RndEGraph EClassId @@ -353,9 +356,8 @@ insertFitness c f p pure c -while p arg prog = do b <- p - when b do arg' <- prog arg - while p arg' prog +while p arg prog = do when (p arg) do arg' <- prog arg + while p arg' prog {- egraphGP :: SRMatrix -> PVector -> [Fix SRTree] -> Int -> RndEGraph (Fix SRTree, Double) diff --git a/apps/srtools/IO.hs b/apps/srtools/IO.hs --- a/apps/srtools/IO.hs +++ b/apps/srtools/IO.hs @@ -51,7 +51,7 @@ (Nothing, _) -> _expr basic (_, Nothing) -> _expr basic (Just xV, Just yV) -> _expr $ getBasicStats args' seed dset{_xTr = xV, _yTr = yV} tree theta0 ix - sseOrig = getSSE dset tree + sseOrig = getSSE dset t sseOpt = getSSE dset (_expr basic) info = getInfo args' dset (_expr basic) treeVal cis = getCI args' dset basic (alpha args') diff --git a/src/Algorithm/SRTree/AD.hs b/src/Algorithm/SRTree/AD.hs --- a/src/Algorithm/SRTree/AD.hs +++ b/src/Algorithm/SRTree/AD.hs @@ -21,10 +21,11 @@ ( forwardMode , forwardModeUnique , reverseModeUnique + , reverseModeUniqueArr , forwardModeUniqueJac ) where -import Control.Monad (forM_) +import Control.Monad (forM_, foldM) import Control.Monad.ST ( runST ) import Data.Bifunctor (bimap, first, second) import qualified Data.DList as DL @@ -41,6 +42,8 @@ import qualified Data.Vector as V import Debug.Trace (trace, traceShow) import GHC.IO (unsafePerformIO) +import qualified Data.IntMap.Strict as IntMap +import Data.List ( foldl' ) applyUni :: (Index ix, Source r e, Floating e, Floating b) => Function -> Either (Array r ix e) b -> Either (Array D ix e) b applyUni f (Left t) = @@ -292,6 +295,179 @@ combine j (Uni f gs) s = gs combine j (Bin op l r) s = l+r +-- | Same as above, but using reverse mode with the tree encoded as an array, that is even faster. +--reverseModeUniqueArr :: SRMatrix +-- -> PVector +-- -> SRVector +-- -> (SRVector -> SRVector) +-- -> Array S Ix1 (Int, Int, Int, Double) -- arity, opcode, ix, const val +-- -> (Array D Ix1 Double, Array S Ix1 Double) +reverseModeUniqueArr xss theta ys f t j2ix = + {-let fwd = forward + v = fwd IntMap.! 0 + err = f v - delay ys + partial = reverseMode fwd + in -} + unsafePerformIO $ do + fwd <- M.newMArray (Sz2 m n) 0 + partial <- M.newMArray (Sz2 m n) 0 + jacob <- M.newMArray (Sz p) 0 + fwd' <- UMA.unsafeFreeze (getComp xss) fwd + let v = fwd' M. IO () + forward fwd = forM_ (Prelude.reverse t) makeFwd + where + makeFwd (j, (0, 0, ix, _)) = do let j' = j2ix IntMap.! j + forM_ [0..m-1] $ \i -> do + let val = xss M.! (i :. ix) + UMA.unsafeWrite fwd (i :. j') val + makeFwd (j, (0, 1, ix, _)) = do let j' = j2ix IntMap.! j + v = theta M.! ix + forM_ [0..m-1] $ \i -> do + UMA.unsafeWrite fwd (i :. j') v + makeFwd (j, (0, 2, _, x)) = do let j' = j2ix IntMap.! j + forM_ [0..m-1] $ \i -> do + UMA.unsafeWrite fwd (i :. j') x + makeFwd (j, (1, f, _, _)) = do let j' = j2ix IntMap.! j + j2 = j2ix IntMap.! (2*j + 1) + forM_ [0..m-1] $ \i -> do + v <- UMA.unsafeRead fwd (i :. j2) + let val = evalFun (toEnum f) v + UMA.unsafeWrite fwd (i :. j') val + makeFwd (j, (2, op, _, _)) = do let j' = j2ix IntMap.! j + j2 = j2ix IntMap.! (2*j + 1) + j3 = j2ix IntMap.! (2*j + 2) + forM_ [0..m-1] $ \i -> do + l <- UMA.unsafeRead fwd (i :. j2) + r <- UMA.unsafeRead fwd (i :. j3) + let val = evalOp (toEnum op) l r + UMA.unsafeWrite fwd (i :. j') val + {- + forward = foldr (makeFwd) IntMap.empty (IntMap.toAscList t) + where + makeFwd (j, (0, 0, ix, _)) fwd = IntMap.insert j (xss M. MArray (PrimState IO) S Ix2 Double -> IO () + reverseMode fwd partial = do forM_ [0..m-1] $ \i -> UMA.unsafeWrite partial (i :. 0) 1 + forM_ t makeRev + where + makeRev (j, (1, f, _, _)) = do forM_ [0..m-1] $ \i -> do + let dxj = j2ix IntMap.! j + vj = j2ix IntMap.! (2*j + 1) + v <- UMA.unsafeRead fwd (i :. vj) + dx <- UMA.unsafeRead partial (i :. dxj) + let val = dx * derivative (toEnum f) v + UMA.unsafeWrite partial (i :. vj) val + makeRev (j, (2, op, _, _)) = do forM_ [0..m-1] $ \i -> do + let dxj = j2ix IntMap.! j + lj = j2ix IntMap.! (2*j + 1) + rj = j2ix IntMap.! (2*j + 2) + l <- UMA.unsafeRead fwd (i :. lj) + r <- UMA.unsafeRead fwd (i :. rj) + dx <- UMA.unsafeRead partial (i :. dxj) + let (dxl, dxr) = diff (toEnum op) dx l r + UMA.unsafeWrite partial (i :. lj) dxl + UMA.unsafeWrite partial (i :. rj) dxr + makeRev _ = pure () + {- + reverseMode fwd = foldr (makeRev) rev0 (IntMap.toDescList t) + where + rev0 = IntMap.insert 0 (M.replicate (getComp xss) (M.Sz m) 1) IntMap.empty + + + makeRev (j, (1, f, _, _)) rev = let v = fwd IntMap.! (2*j + 1) + dx = rev IntMap.! j + val = dx !*! (M.map (derivative (toEnum f)) v) + in IntMap.insert (2*j + 1) val rev + makeRev (j, (2, op, _, _)) rev = let l = fwd IntMap.! (2*j + 1) + r = fwd IntMap.! (2*j + 2) + dx = rev IntMap.! j + (dxl, dxr) = diff (toEnum op) dx l r + in IntMap.insert (2*j + 2) dxr $ IntMap.insert (2*j + 1) dxl rev + makeRev (j, _) rev = rev + -} + + -- dx is the current derivative so far + -- fx is the evaluation of the left branch + -- gx is the evaluation of the right branch + -- + -- this should return a tuple, where the left element is + -- dx * d op(f(x), g(x)) / d f(x) and + -- the right branch dx * d op (f(x), g(x)) / d g(x) + arr1 !**! arr2 = M.zipWith (**) arr1 arr2 + + diff Add dx fx gy = (dx, dx) + diff Sub dx fx gy = (dx, negate dx) + diff Mul dx fx gy = (dx * gy, dx * fx) + diff Div dx fx gy = (dx / gy, dx * (negate fx / (gy * gy))) + diff Power dx fx gy = let dxl = dx * (fx ** (gy-1)) + dv2 = fx * log fx + in (dxl * gy, dxl * dv2) + diff PowerAbs dx fx gy = let dxl = (gy * fx) * (fx ** abs (gy - 2)) + dxr = (log (abs fx)) * (fx ** abs gy) + in (dxl * dx, dxr * dx) + diff AQ dx fx gy = let dxl = recip ((sqrt . (+1)) (gy * gy)) + dxy = fx * gy * (dxl^3) -- / (sqrt (gy*gy + 1)) + in (dxl * dx, dxy * dx) + {- + diff Mul dx fx gy = (dx !*! gy, dx !*! fx) + diff Div dx fx gy = (dx !/! gy, dx !*! (M.map negate fx !/! (gy !*! gy))) + diff Power dx fx gy = let dxl = dx !*! (fx !**! (M.map (subtract 1) gy)) + dv2 = fx !*! M.map log fx + in (dxl !*! gy, dxl !*! dv2) + diff PowerAbs dx fx gy = let dxl = (gy !*! fx) !*! (fx !**! M.map abs (M.map (subtract 2) gy)) + dxr = (M.map log (M.map abs fx)) !*! (fx !**! M.map abs gy) + in (dxl !*! dx, dxr !*! dx) + diff AQ dx fx gy = let dxl = M.map recip (M.map (sqrt . (+1)) (gy !*! gy)) + dxy = fx !*! gy !*! (M.map (^3) dxl) -- / (sqrt (gy*gy + 1)) + in (dxl !*! dx, dxy !*! dx) + -} + + -- once we reach a leaf with a parameter, we return a singleton + -- with that derivative upwards until the root + combine :: MArray (PrimState IO) S Ix2 Double -> MArray (PrimState IO) S Ix1 Double -> Array S Ix1 Double -> IO () + combine partial jacob err = forM_ t makeJacob + where + makeJacob (j, (0, 1, ix, _)) = do let j' = j2ix IntMap.! j + addI a b acc = do let v1 = err M.! a + v2 <- UMA.unsafeRead partial (a :. b) + pure (v1*v2 + acc) + acc <- foldM (\a i -> addI i j' a) 0 [0..m-1] + UMA.unsafeWrite jacob ix acc + makeJacob _ = pure () + {- + combine :: IntMap.IntMap (Array D Ix1 Double) -> MArray (PrimState IO) S Ix1 Double -> Array D Ix1 Double -> IO () + combine partial jacob err = forM_ (IntMap.toAscList t) makeJacob + where + makeJacob (j, (0, 1, ix, _)) = do v <- dotM (partial IntMap.! j) err + UMA.unsafeWrite jacob ix v + makeJacob _ = pure () + -} -- | The function `forwardModeUnique` calculates the numerical gradient of the tree and evaluates the tree at the same time. It assumes that each parameter has a unique occurrence in the expression. This should be significantly faster than `forwardMode`. forwardModeUniqueJac :: SRMatrix -> PVector -> Fix SRTree -> [PVector] diff --git a/src/Algorithm/SRTree/Likelihoods.hs b/src/Algorithm/SRTree/Likelihoods.hs --- a/src/Algorithm/SRTree/Likelihoods.hs +++ b/src/Algorithm/SRTree/Likelihoods.hs @@ -22,6 +22,7 @@ , nll , predict , gradNLL + , gradNLLArr , gradNLLNonUnique , fisherNLL , getSErr @@ -29,13 +30,14 @@ ) where -import Algorithm.SRTree.AD ( forwardMode, reverseModeUnique ) -- ( reverseModeUnique ) +import Algorithm.SRTree.AD ( forwardMode, reverseModeUnique, reverseModeUniqueArr ) -- ( reverseModeUnique ) import Data.Massiv.Array hiding (all, map, read, replicate, tail, take, zip) import qualified Data.Massiv.Array as M import Data.Maybe (fromMaybe) import Data.SRTree (Fix (..), SRTree (..), floatConstsToParam, relabelParams) import Data.SRTree.Derivative (deriveByParam) import Data.SRTree.Eval (PVector, SRMatrix, SRVector, compMode, evalTree) +import qualified Data.IntMap.Strict as IntMap -- | Supported distributions for negative log-likelihood data Distribution = Gaussian | Bernoulli | Poisson @@ -170,6 +172,38 @@ | otherwise = (nll' Poisson 1.0 yhat ys, delay grad) where (yhat, grad) = reverseModeUnique xss theta ys exp tree + --err = exp yhat - ys + +-- | Gradient of the negative log-likelihood +--Array B Ix1 (Int, Int, Int, Double) +gradNLLArr :: Distribution -> Maybe Double -> SRMatrix -> PVector -> [(Int,(Int, Int, Int, Double))] -> IntMap.IntMap Int -> PVector -> (Double, SRVector) +gradNLLArr Gaussian msErr xss ys tree j2ix theta = + (nll' Gaussian sErr yhat ys', delay grad ./ (sErr * sErr)) + where + (Sz m) = M.size ys + (Sz p) = M.size theta + ys' = delay ys + (yhat, grad) = reverseModeUniqueArr xss theta ys' id tree j2ix + -- err = yhat - delay ys + --ssr = sse xss ys tree theta + est = sqrt $ fromIntegral (m - p) -- $ ssr / fromIntegral (m - p) + sErr = getSErr Gaussian est msErr + +gradNLLArr Bernoulli _ xss (delay -> ys) tree j2ix theta + | M.any (\x -> x /= 0 && x /= 1) ys = error "For Bernoulli distribution the output must be either 0 or 1." + | otherwise = (nll' Bernoulli 1.0 yhat ys, delay grad) + where + (yhat, grad) = reverseModeUniqueArr xss theta ys logistic tree j2ix + grad' = M.map nanTo0 grad + --err = logistic yhat - ys + nanTo0 x = if isNaN x then 0 else x + +gradNLLArr Poisson _ xss (delay -> ys) tree j2ix theta + | M.any (<0) ys = error "For Poisson distribution the output must be non-negative." + -- | M.any isNaN grad = error $ "NaN gradient " <> show grad + | otherwise = (nll' Poisson 1.0 yhat ys, delay grad) + where + (yhat, grad) = reverseModeUniqueArr xss theta ys exp tree j2ix --err = exp yhat - ys -- | Gradient of the negative log-likelihood diff --git a/src/Algorithm/SRTree/Opt.hs b/src/Algorithm/SRTree/Opt.hs --- a/src/Algorithm/SRTree/Opt.hs +++ b/src/Algorithm/SRTree/Opt.hs @@ -1,3 +1,4 @@ +{-# LANGUAGE BangPatterns #-} ----------------------------------------------------------------------------- -- | -- Module : Algorithm.SRTree.Opt @@ -17,10 +18,38 @@ import Algorithm.SRTree.NonlinearOpt import Data.Bifunctor (bimap, second) import Data.Massiv.Array -import Data.SRTree (Fix (..), SRTree (..), floatConstsToParam, relabelParams) +import Data.SRTree (Fix (..), SRTree (..), floatConstsToParam, relabelParams, countNodes) import Data.SRTree.Eval (evalTree, compMode) import qualified Data.Vector.Storable as VS +import qualified Data.IntMap.Strict as IntMap +import Data.SRTree.Recursion +import Debug.Trace + +tree2arr :: Fix SRTree -> IntMap.IntMap (Int, Int, Int, Double) +tree2arr tree = IntMap.fromList listTree + where + height = cata alg + where + alg (Var ix) = 1 + alg (Const x) = 1 + alg (Param ix) = 1 + alg (Uni _ t) = 1 + t + alg (Bin _ l r) = 1 + max l r + listTree = accu indexer convert tree 0 + + indexer (Var ix) iy = Var ix + indexer (Const x) iy = Const x + indexer (Param ix) iy = Param ix + indexer (Bin op l r) iy = Bin op (l, 2*iy+1) (r, 2*iy+2) + indexer (Uni f t) iy = Uni f (t, 2*iy+1) + + convert (Var ix) iy = [(iy, (0, 0, ix, -1))] + convert (Const x) iy = [(iy, (0, 2, -1, x))] + convert (Param ix) iy = [(iy, (0, 1, ix, -1))] + convert (Uni f t) iy = (iy, (1, fromEnum f, -1, -1)) : t + convert (Bin op l r) iy = (iy, (2, fromEnum op, -1, -1)) : (l <> r) + -- | minimizes the negative log-likelihood of the expression minimizeNLL :: Distribution -> Maybe Double -> Int -> SRMatrix -> PVector -> Fix SRTree -> PVector -> (PVector, Double) minimizeNLL dist msErr niter xss ys tree t0 @@ -30,13 +59,15 @@ where tree' = relabelParams tree -- $ fst $ floatConstsToParam tree t0' = toStorableVector t0 + treeArr = IntMap.toAscList $ tree2arr tree' + j2ix = IntMap.fromList $ Prelude.zip (Prelude.map fst treeArr) [0..] (Sz n) = size t0 (Sz m) = size ys - funAndGrad = second (toStorableVector . computeAs S) . gradNLL dist msErr xss ys tree' . fromStorableVector compMode - (f, _) = gradNLL dist msErr xss ys tree t0 -- if there's no parameter or no iterations + funAndGrad = second (toStorableVector . computeAs S) . gradNLLArr dist msErr xss ys treeArr j2ix . fromStorableVector compMode + (f, _) = gradNLLArr dist msErr xss ys treeArr j2ix t0 -- if there's no parameter or no iterations algorithm = LBFGS funAndGrad Nothing - stop = ObjectiveRelativeTolerance 1e-10 :| [MaximumEvaluations (fromIntegral niter)] + stop = ObjectiveRelativeTolerance 1e-6 :| [MaximumEvaluations (fromIntegral niter)] problem = LocalProblem (fromIntegral n) stop algorithm t_opt = case minimizeLocal problem t0' of Right sol -> solutionParams sol diff --git a/src/Data/SRTree/Eval.hs b/src/Data/SRTree/Eval.hs --- a/src/Data/SRTree/Eval.hs +++ b/src/Data/SRTree/Eval.hs @@ -41,7 +41,7 @@ type SRMatrix = M.Array S Ix2 Double compMode :: M.Comp -compMode = M.Par' +compMode = M.Seq -- Improve quality of life with Num and Floating instances for our matrices instance Index ix => Num (M.Array D ix Double) where diff --git a/src/Numeric/Optimization/NLOPT/Bindings.hs b/src/Numeric/Optimization/NLOPT/Bindings.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Optimization/NLOPT/Bindings.hs @@ -0,0 +1,1065 @@ +{-# OPTIONS_GHC -Wall #-} +{-# LANGUAGE ForeignFunctionInterface #-} +{-# LANGUAGE NoMonomorphismRestriction #-} + +{- | +Module : Numeric.Optimization.NLOPT.Bindings +Copyright : (c) Matthew Peddie 2017 +License : BSD3 +Maintainer : Matthew Peddie +Stability : provisional +Portability : GHC + +Low-level interface to the NLOPT library. Please see + +for detailed information; the Haskell functions in this module closely +follow the interface to the C library in @nlopt.h@. + +Differences between this module and the C interface are documented +here; functions with identical interfaces are not. In general: + + ['Opt'] corresponds to an @nlopt_opt@ object + + ['Result'] corresponds to @nlopt_result@ + + ['V.Vector' 'Double'] corresponds to a @const double *@ input or a + @double *@ output + + ['ScalarFunction'] corresponds to @nlopt_func@ + + ['VectorFunction'] corresponds to @nlopt_mfunc@ + + ['PreconditionerFunction'] corresponds to @nlopt_precond@ + +User data that is handled by @void *@ in the C bindings can be any +Haskell value. + +-} + +module Numeric.Optimization.NLOPT.Bindings ( + -- * C enums + Algorithm(..) + , algorithm_name + , Result(..) + , isSuccess + -- * Optimizer object + , Opt + , create + , destroy + , copy + -- * Random number generator seeding + , srand + , srand_time + -- * Metadata + , Version(..) + , version + , get_algorithm + , get_dimension + -- * Callbacks + , ScalarFunction + , VectorFunction + , PreconditionerFunction + -- * Running the optimizer + , Output(..) + , optimize + -- * Objective function configuration + , set_min_objective + , set_max_objective + , set_precond_min_objective + , set_precond_max_objective + -- * Bound configuration + , set_lower_bounds + , set_lower_bounds1 + , get_lower_bounds + , set_upper_bounds + , set_upper_bounds1 + , get_upper_bounds + -- * Constraint configuration + , remove_inequality_constraints + , add_inequality_constraint + , add_precond_inequality_constraint + , add_inequality_mconstraint + , remove_equality_constraints + , add_equality_constraint + , add_precond_equality_constraint + , add_equality_mconstraint + -- * Stopping criterion configuration + , set_stopval + , get_stopval + , set_ftol_rel + , get_ftol_rel + , set_ftol_abs + , get_ftol_abs + , set_xtol_rel + , get_xtol_rel + , set_xtol_abs1 + , set_xtol_abs + , get_xtol_abs + , set_maxeval + , get_maxeval + , set_maxtime + , get_maxtime + , force_stop + , set_force_stop + , get_force_stop + -- * Algorithm-specific configuration + , set_local_optimizer + , set_population + , get_population + , set_vector_storage + , get_vector_storage + , set_default_initial_step + , set_initial_step + , set_initial_step1 + , get_initial_step + ) where + +import Foreign hiding (void) +import Foreign.C.String +import Foreign.C.Types +import qualified Foreign.Concurrent as CFP + +import qualified Data.Vector.Storable.Mutable as MV +import qualified Data.Vector.Storable as V + +{- C enums -} + +-- | The NLOPT algorithm names, apart from the names of the actual +-- optimization methods, follow this scheme: +-- +-- [@G@] means a global method +-- [@L@] means a local method +-- [@D@] means a method that requires the derivative +-- [@N@] means a method that does not require the derivative +-- [@*_RAND@] means the algorithm involves some randomization. +-- [@*_NOSCAL@] means the algorithm is *not* scaled to a unit +-- hypercube (i.e. it is sensitive to the units of x) +data Algorithm + = GN_DIRECT -- ^ DIviding RECTangles + | GN_DIRECT_L -- ^ DIviding RECTangles, + -- locally-biased variant + | GN_DIRECT_L_RAND -- ^ DIviding RECTangles, "slightly + -- randomized" + | GN_DIRECT_NOSCAL -- ^ DIviding RECTangles, unscaled version + | GN_DIRECT_L_NOSCAL -- ^ DIviding RECTangles, + -- locally-biased and unscaled + | GN_DIRECT_L_RAND_NOSCAL -- ^ DIviding RECTangles, locally-biased, + -- unscaled and "slightly randomized" + | GN_ORIG_DIRECT -- ^ DIviding RECTangles, original FORTRAN + -- implementation + | GN_ORIG_DIRECT_L -- ^ DIviding RECTangles, + -- locally-biased, original FORTRAN + -- implementation + | GD_STOGO -- ^ Stochastic Global Optimization + | GD_STOGO_RAND -- ^ Stochastic Global Optimization, + -- randomized variant + | LD_LBFGS_NOCEDAL -- ^ Limited-memory BFGS + | LD_LBFGS -- ^ Limited-memory BFGS + | LN_PRAXIS -- ^ PRincipal AXIS gradient-free local + -- optimization + | LD_VAR2 -- ^ Shifted limited-memory + -- variable-metric, rank-2 + | LD_VAR1 -- ^ Shifted limited-memory + -- variable-metric, rank-1 + | LD_TNEWTON -- ^ Truncated Newton's method + | LD_TNEWTON_RESTART -- ^ Truncated Newton's method with + -- automatic restarting + | LD_TNEWTON_PRECOND -- ^ Preconditioned truncated Newton's + -- method + | LD_TNEWTON_PRECOND_RESTART -- ^ Preconditioned truncated Newton's + -- method with automatic restarting + | GN_CRS2_LM -- ^ Controlled Random Search with + -- Local Mutation + | GN_MLSL -- ^ Original Multi-Level + -- Single-Linkage + | GD_MLSL -- ^ Original Multi-Level + -- Single-Linkage, user-provided + -- derivative + | GN_MLSL_LDS -- ^ Multi-Level Single-Linkage with + -- Sobol Low-Discrepancy Sequence for + -- starting points + | GD_MLSL_LDS -- ^ Multi-Level Single-Linkage with + -- Sobol Low-Discrepancy Sequence for + -- starting points, user-provided + -- derivative + | LD_MMA -- ^ Method of moving averages + | LN_COBYLA -- ^ Constrained Optimization BY Linear + -- Approximations + | LN_NEWUOA -- ^ Powell's NEWUOA algorithm + | LN_NEWUOA_BOUND -- ^ Powell's NEWUOA algorithm with + -- bounds by SGJ + | LN_NELDERMEAD -- ^ Nelder-Mead Simplex gradient-free + -- method + | LN_SBPLX -- ^ NLOPT implementation of Rowan's + -- Subplex algorithm + | LN_AUGLAG -- ^ AUGmented LAGrangian + | LD_AUGLAG -- ^ AUGmented LAGrangian, + -- user-provided derivative + | LN_AUGLAG_EQ -- ^ AUGmented LAGrangian with penalty + -- functions only for equality + -- constraints + | LD_AUGLAG_EQ -- ^ AUGmented LAGrangian with + -- penalty functions only for equality + -- constraints, user-provided + -- derivative + | LN_BOBYQA -- ^ Bounded Optimization BY Quadratic + -- Approximations + | GN_ISRES -- ^ Improved Stochastic Ranking + -- Evolution Strategy + + | AUGLAG -- ^ AUGmented LAGrangian, requires + -- local_optimizer to be set + | AUGLAG_EQ -- ^ AUGmented LAGrangian with penalty + -- functions only for equality + -- constraints, requires + -- local_optimizer to be set + | G_MLSL -- ^ Original Multi-Level + -- Single-Linkage, user-provided + -- derivative, requires local_optimizer + -- to be set + | G_MLSL_LDS -- ^ Multi-Level Single-Linkage with + -- Sobol Low-Discrepancy Sequence for + -- starting points, requires + -- local_optimizer to be set + | LD_SLSQP -- ^ Sequential Least-SQuares Programming + | LD_CCSAQ -- ^ Conservative Convex Separable + -- Approximation + | GN_ESCH -- ^ Evolutionary Algorithm + deriving (Eq, Show, Read, Bounded) + +instance Enum Algorithm where + fromEnum GN_DIRECT = 0 + fromEnum GN_DIRECT_L = 1 + fromEnum GN_DIRECT_L_RAND = 2 + fromEnum GN_DIRECT_NOSCAL = 3 + fromEnum GN_DIRECT_L_NOSCAL = 4 + fromEnum GN_DIRECT_L_RAND_NOSCAL = 5 + fromEnum GN_ORIG_DIRECT = 6 + fromEnum GN_ORIG_DIRECT_L = 7 + fromEnum GD_STOGO = 8 + fromEnum GD_STOGO_RAND = 9 + fromEnum LD_LBFGS_NOCEDAL = 10 + fromEnum LD_LBFGS = 11 + fromEnum LN_PRAXIS = 12 + fromEnum LD_VAR2 = 13 + fromEnum LD_VAR1 = 14 + fromEnum LD_TNEWTON = 15 + fromEnum LD_TNEWTON_RESTART = 16 + fromEnum LD_TNEWTON_PRECOND = 17 + fromEnum LD_TNEWTON_PRECOND_RESTART = 18 + fromEnum GN_CRS2_LM = 19 + fromEnum GN_MLSL = 20 + fromEnum GD_MLSL = 21 + fromEnum GN_MLSL_LDS = 22 + fromEnum GD_MLSL_LDS = 23 + fromEnum LD_MMA = 24 + fromEnum LN_COBYLA = 25 + fromEnum LN_NEWUOA = 26 + fromEnum LN_NEWUOA_BOUND = 27 + fromEnum LN_NELDERMEAD = 28 + fromEnum LN_SBPLX = 29 + fromEnum LN_AUGLAG = 30 + fromEnum LD_AUGLAG = 31 + fromEnum LN_AUGLAG_EQ = 32 + fromEnum LD_AUGLAG_EQ = 33 + fromEnum LN_BOBYQA = 34 + fromEnum GN_ISRES = 35 + fromEnum AUGLAG = 36 + fromEnum AUGLAG_EQ = 37 + fromEnum G_MLSL = 38 + fromEnum G_MLSL_LDS = 39 + fromEnum LD_SLSQP = 40 + fromEnum LD_CCSAQ = 41 + fromEnum GN_ESCH = 42 + toEnum 0 = GN_DIRECT + toEnum 1 = GN_DIRECT_L + toEnum 2 = GN_DIRECT_L_RAND + toEnum 3 = GN_DIRECT_NOSCAL + toEnum 4 = GN_DIRECT_L_NOSCAL + toEnum 5 = GN_DIRECT_L_RAND_NOSCAL + toEnum 6 = GN_ORIG_DIRECT + toEnum 7 = GN_ORIG_DIRECT_L + toEnum 8 = GD_STOGO + toEnum 9 = GD_STOGO_RAND + toEnum 10 = LD_LBFGS_NOCEDAL + toEnum 11 = LD_LBFGS + toEnum 12 = LN_PRAXIS + toEnum 13 = LD_VAR2 + toEnum 14 = LD_VAR1 + toEnum 15 = LD_TNEWTON + toEnum 16 = LD_TNEWTON_RESTART + toEnum 17 = LD_TNEWTON_PRECOND + toEnum 18 = LD_TNEWTON_PRECOND_RESTART + toEnum 19 = GN_CRS2_LM + toEnum 20 = GN_MLSL + toEnum 21 = GD_MLSL + toEnum 22 = GN_MLSL_LDS + toEnum 23 = GD_MLSL_LDS + toEnum 24 = LD_MMA + toEnum 25 = LN_COBYLA + toEnum 26 = LN_NEWUOA + toEnum 27 = LN_NEWUOA_BOUND + toEnum 28 = LN_NELDERMEAD + toEnum 29 = LN_SBPLX + toEnum 30 = LN_AUGLAG + toEnum 31 = LD_AUGLAG + toEnum 32 = LN_AUGLAG_EQ + toEnum 33 = LD_AUGLAG_EQ + toEnum 34 = LN_BOBYQA + toEnum 35 = GN_ISRES + toEnum 36 = AUGLAG + toEnum 37 = AUGLAG_EQ + toEnum 38 = G_MLSL + toEnum 39 = G_MLSL_LDS + toEnum 40 = LD_SLSQP + toEnum 41 = LD_CCSAQ + toEnum 42 = GN_ESCH + toEnum e = error $ + "Algorithm.toEnum: invalid C value '" ++ show e ++ "' received." + +foreign import ccall "nlopt.h nlopt_algorithm_name" + nlopt_algorithm_name :: CInt -> CString + +algorithm_name :: Algorithm -> IO String +algorithm_name = peekCString . nlopt_algorithm_name . fromIntegral . fromEnum + +-- | Mostly self-explanatory. +data Result + = FAILURE -- ^ Generic failure code + | INVALID_ARGS + | OUT_OF_MEMORY + | ROUNDOFF_LIMITED + | FORCED_STOP + | SUCCESS -- ^ Generic success code + | STOPVAL_REACHED + | FTOL_REACHED + | XTOL_REACHED + | MAXEVAL_REACHED + | MAXTIME_REACHED + deriving (Eq, Read, Show, Bounded) + +instance Enum Result where + fromEnum FAILURE = -1 + fromEnum INVALID_ARGS = -2 + fromEnum OUT_OF_MEMORY = -3 + fromEnum ROUNDOFF_LIMITED = -4 + fromEnum FORCED_STOP = -5 + fromEnum SUCCESS = 1 + fromEnum STOPVAL_REACHED = 2 + fromEnum FTOL_REACHED = 3 + fromEnum XTOL_REACHED = 4 + fromEnum MAXEVAL_REACHED = 5 + fromEnum MAXTIME_REACHED = 6 + toEnum (-1) = FAILURE + toEnum (-2) = INVALID_ARGS + toEnum (-3) = OUT_OF_MEMORY + toEnum (-4) = ROUNDOFF_LIMITED + toEnum (-5) = FORCED_STOP + toEnum 1 = SUCCESS + toEnum 2 = STOPVAL_REACHED + toEnum 3 = FTOL_REACHED + toEnum 4 = XTOL_REACHED + toEnum 5 = MAXEVAL_REACHED + toEnum 6 = MAXTIME_REACHED + toEnum e = error $ + "Result.toEnum: invalid C value '" ++ show e ++ "' received." + +isSuccess :: Result -> Bool +isSuccess SUCCESS = True +isSuccess STOPVAL_REACHED = True +isSuccess FTOL_REACHED = True +isSuccess XTOL_REACHED = True +isSuccess MAXEVAL_REACHED = True +isSuccess MAXTIME_REACHED = True +isSuccess _ = False + +parseEnum :: (Integral a, Enum b) => a -> b +parseEnum = toEnum . fromIntegral + +{- NLOPT optimizer object -} + +type NloptOpt = Ptr () + +-- | An optimizer object which must be created, configured and then +-- passed to 'optimize' to solve a problem +newtype Opt = Opt { pointerFromOpt :: ForeignPtr () } + +withOpt :: Opt -> (NloptOpt -> IO a) -> IO a +withOpt (Opt p) f = do + ret <- withForeignPtr p f + touchForeignPtr p -- This is critical! Otherwise the GC might + -- think it's done with everything in the middle + -- of the problem. + return ret + +useOpt :: (NloptOpt -> IO a) -> Opt -> IO a +useOpt = flip withOpt + +-- Every time we make a "wrapper" call, the runtime allocates a new +-- function pointer and won't release it until we explicitly tell it +-- to. This doesn't mesh well with NLOPT's "object-oriented" design, +-- wherein we have to allocate an object and make a bunch of setup +-- calls before we run the problem, so what we do is add a finalizer +-- to the 'Opt' object's 'ForeignPtr' every time we need to create a +-- function pointer for C to use. +addFunPtrFinalizer :: Opt -> FunPtr a -> IO () +addFunPtrFinalizer (Opt p) funptr = + CFP.addForeignPtrFinalizer p (freeHaskellFunPtr funptr) + +foreign import ccall "nlopt.h nlopt_create" + nlopt_create :: CInt -> CUInt -> IO (NloptOpt) + +-- | Create a new 'Opt' object +create :: Algorithm -- ^ Choice of algorithm + -> Word -- ^ Parameter vector dimension + -> IO (Maybe Opt) -- ^ Optimizer object +create alg dimension = do + outp <- nlopt_create (fromIntegral $ fromEnum alg) (fromIntegral dimension) + if (outp == nullPtr) + then return Nothing + else Just . Opt <$> CFP.newForeignPtr outp (nlopt_destroy outp) + +foreign import ccall "nlopt.h nlopt_destroy" + nlopt_destroy :: NloptOpt -> IO () + +-- It shouldn't be strictly necessary to call this by hand since we've +-- already put a call to 'nlopt_destroy' into the 'ForeignPtr', but +-- it's available in the C interface. +destroy :: Opt -> IO () +destroy = finalizeForeignPtr . pointerFromOpt + +foreign import ccall "nlopt.h nlopt_copy" + nlopt_copy :: NloptOpt -> IO (NloptOpt) + +copy :: Opt -> IO Opt +copy = useOpt $ \inp -> do + outp <- nlopt_copy inp + Opt <$> CFP.newForeignPtr outp (nlopt_destroy outp) + +{- Random seeding functions -} + +foreign import ccall "nlopt.h nlopt_srand" + nlopt_srand :: CUInt -> IO () + +srand :: Integral a => a -> IO () +srand = nlopt_srand . fromIntegral + +foreign import ccall "nlopt.h nlopt_srand_time" + nlopt_srand_time :: IO () + +srand_time :: IO () +srand_time = nlopt_srand_time + +{- Metadata -} + +foreign import ccall "nlopt.h nlopt_version" + nlopt_version :: Ptr CInt -> Ptr CInt -> Ptr CInt -> IO () + +-- | NLOPT library version, e.g. @2.4.2@ +data Version = Version + { major :: Int + , minor :: Int + , bugfix :: Int + } deriving (Eq, Ord, Read, Show) + +version :: IO Version +version = + alloca $ \majptr -> + alloca $ \minptr -> + alloca $ \bfptr -> do + nlopt_version majptr minptr bfptr + Version <$> pk majptr <*> pk minptr <*> pk bfptr + where + pk = fmap fromIntegral . peek + +foreign import ccall "nlopt.h nlopt_get_algorithm" + nlopt_get_algorithm :: NloptOpt -> IO CInt + +get_algorithm :: Opt -> IO Algorithm +get_algorithm = useOpt $ fmap parseEnum . nlopt_get_algorithm + +foreign import ccall "nlopt.h nlopt_get_dimension" + nlopt_get_dimension :: NloptOpt -> IO CUInt + +get_dimension :: Opt -> IO Word +get_dimension = useOpt $ fmap fromIntegral . nlopt_get_dimension + +{- Callback functions -} + +asMVector :: CUInt -> Ptr CDouble -> IO (MV.IOVector Double) +asMVector dim ptr = + MV.unsafeCast . flip MV.unsafeFromForeignPtr0 (fromIntegral dim) <$> + newForeignPtr_ ptr + +asVector :: CUInt -> Ptr CDouble -> IO (V.Vector Double) +asVector dim ptr = + V.unsafeCast . flip V.unsafeFromForeignPtr0 (fromIntegral dim) <$> + newForeignPtr_ ptr + +type CFunc a = CUInt -> Ptr CDouble -> Ptr CDouble -> StablePtr a -> IO CDouble + +-- | This function type corresponds to @nlopt_func@ in C and is used +-- for scalar functions of the parameter vector. You may pass data of +-- any type @a@ to the functions in this module that take a +-- 'ScalarFunction' as an argument; this data will be supplied to your +-- your function when it is called. +type ScalarFunction a + = V.Vector Double -- ^ Parameter vector + -> Maybe (MV.IOVector Double) -- ^ Gradient vector to be filled in + -> a -- ^ User data + -> IO Double -- ^ Scalar result + +-- | This function type corresponds to @nlopt_mfunc@ in C and is used +-- for vector functions of the parameter vector. You may pass data of +-- any type @a@ to the functions in this module that take a +-- 'VectorFunction' as an argument; this data will be supplied to your +-- function when it is called. +type VectorFunction a + = V.Vector Double -- ^ Parameter vector + -> MV.IOVector Double -- ^ Output vector to be filled in + -> Maybe (MV.IOVector Double) -- ^ Gradient vector to be filled in + -> a -- ^ User data + -> IO () + +-- | This function type corresponds to @nlopt_precond@ in C and is +-- used for functions that precondition a vector at a given point in +-- the parameter space. You may pass data of any type @a@ to the +-- functions in this module that take a 'PreconditionerFunction' as an +-- argument; this data will be supplied to your function when it is +-- called. +type PreconditionerFunction a + = V.Vector Double -- ^ Parameter vector + -> V.Vector Double -- ^ Vector @v@ to precondition + -> MV.IOVector Double -- ^ Output vector @vpre@ to be filled in + -> a -- ^ User data + -> IO () + +wrapCFunction :: ScalarFunction a -> CFunc a +wrapCFunction cfunc dim stateptr gradientptr userptr = do + nloptgradient <- asMVector dim gradientptr + statevec <- asVector dim stateptr + userdata <- deRefStablePtr userptr + let + gradptr = if gradientptr /= nullPtr + then Just nloptgradient + else Nothing + realToFrac <$> cfunc statevec gradptr userdata + +foreign import ccall safe "wrapper" + mkCFunction :: CFunc a -> IO (FunPtr (CFunc a)) + +type CMFunc a = CUInt -> Ptr CDouble -> CUInt -> Ptr CDouble + -> Ptr CDouble -> StablePtr a -> IO () + +wrapMFunction :: VectorFunction a -> CMFunc a +wrapMFunction mfunc constrdim constrptr dim stateptr gradientptr userptr + = do + nloptgradient <- asMVector (dim * constrdim) gradientptr + nloptconstraint <- asMVector constrdim constrptr + statevec <- asVector dim stateptr + userdata <- deRefStablePtr userptr + let + gradptr = if gradientptr /= nullPtr + then Just nloptgradient + else Nothing + mfunc statevec nloptconstraint gradptr userdata + +foreign import ccall safe "wrapper" + mkMFunction :: CMFunc a -> IO (FunPtr (CMFunc a)) + +type CPrecond a = CUInt -> Ptr CDouble -> Ptr CDouble + -> Ptr CDouble -> StablePtr a -> IO () + +wrapPreconditioner :: PreconditionerFunction a -> CPrecond a +wrapPreconditioner prec dim stateptr vptr preptr userptr = do + nloptpre <- asMVector dim preptr + statevec <- asVector dim stateptr + vvec <- asVector dim vptr + userdata <- deRefStablePtr userptr + prec statevec vvec nloptpre userdata + +foreign import ccall safe "wrapper" + mkPreconditionerFunction :: CPrecond a -> IO (FunPtr (CPrecond a)) + +-- We have to do the same silly dance with our user-data 'StablePtr's +-- as we do with function pointer wrappers: because NLOPT expects +-- these pointers before the actual optimization run, we have to +-- attach finalizers for them to the 'Opt' object so that they get +-- cleaned up properly. +addStablePtrFinalizer :: Opt -> StablePtr a -> IO () +addStablePtrFinalizer (Opt p) sp = + CFP.addForeignPtrFinalizer p (freeStablePtr sp) + +getStablePtr :: Opt -> a -> IO (StablePtr a) +getStablePtr opt a = do + aptr <- newStablePtr a + addStablePtrFinalizer opt aptr + return aptr + +exportFunPtr :: (t1 -> IO (FunPtr a)) -> (t -> t1) -> t -> Opt -> IO (FunPtr a) +exportFunPtr mk wrap fun opt = do + funptr <- mk $ wrap fun + addFunPtrFinalizer opt funptr + return funptr + +{- Invoking the optimizer -} + +-- | The output of an NLOPT optimizer run. +data Output = Output + { resultCode :: Result -- ^ Return code + , resultCost :: Double -- ^ Minimum of the objective + -- function if optimization + -- succeeded + , resultParameters :: V.Vector Double -- ^ Parameters corresponding + -- to the minimum if + -- optimization succeeded + } + +foreign import ccall "nlopt.h nlopt_optimize" + nlopt_optimize :: NloptOpt -> Ptr CDouble -> Ptr CDouble -> IO CInt + +-- | This function is very similar to the C function @nlopt_optimize@, +-- but it does not use mutable vectors and returns an 'Output' +-- structure. +optimize :: Opt -- ^ Optimizer object set up to solve the problem + -> V.Vector Double -- ^ Initial-guess parameter vector + -> IO Output -- ^ Results of the optimization run +optimize optimizer x0 = withOpt optimizer $ \opt -> do + vmut <- V.thaw $ V.unsafeCast x0 + alloca $ \costPtr -> do + result <- MV.unsafeWith vmut $ \xptr -> + parseEnum <$> nlopt_optimize opt xptr costPtr + outputCost <- peek . castPtr $ costPtr + iceout <- V.unsafeFreeze (MV.unsafeCast vmut) + return $ Output result outputCost iceout + +{- Objective function setup -} + +foreign import ccall "nlopt.h nlopt_set_min_objective" + nlopt_set_min_objective :: NloptOpt -> FunPtr (CFunc a) + -> StablePtr a -> IO CInt + +foreign import ccall "nlopt.h nlopt_set_max_objective" + nlopt_set_max_objective :: NloptOpt -> FunPtr (CFunc a) + -> StablePtr a -> IO CInt + +set_min_objective :: Opt -> ScalarFunction a -> a -> IO Result +set_min_objective opt objf userdata = do + objfunptr <- exportFunPtr mkCFunction wrapCFunction objf opt + userptr <- getStablePtr opt userdata + withOpt opt $ \o -> + parseEnum <$> + nlopt_set_min_objective o objfunptr userptr + +set_max_objective :: Opt -> ScalarFunction a -> a -> IO Result +set_max_objective opt objf userdata = do + objfunptr <- exportFunPtr mkCFunction wrapCFunction objf opt + userptr <- getStablePtr opt userdata + withOpt opt $ \o -> + parseEnum <$> nlopt_set_max_objective o objfunptr userptr + +foreign import ccall "nlopt.h nlopt_set_precond_min_objective" + nlopt_set_precond_min_objective :: NloptOpt + -> FunPtr (CFunc a) + -> FunPtr (CPrecond a) + -> StablePtr a + -> IO CInt + +foreign import ccall "nlopt.h nlopt_set_precond_max_objective" + nlopt_set_precond_max_objective :: NloptOpt + -> FunPtr (CFunc a) + -> FunPtr (CPrecond a) + -> StablePtr a + -> IO CInt + +set_precond_min_objective :: Opt + -> ScalarFunction a + -> PreconditionerFunction a + -> a + -> IO Result +set_precond_min_objective opt objf pref userdata = do + objfunptr <- exportFunPtr mkCFunction wrapCFunction objf opt + prefunptr <- exportFunPtr mkPreconditionerFunction wrapPreconditioner pref opt + userptr <- getStablePtr opt userdata + withOpt opt $ \o -> parseEnum <$> + nlopt_set_precond_min_objective o objfunptr prefunptr userptr + +set_precond_max_objective :: Opt + -> ScalarFunction a + -> PreconditionerFunction a + -> a + -> IO Result +set_precond_max_objective opt objf pref userdata = do + objfunptr <- exportFunPtr mkCFunction wrapCFunction objf opt + prefunptr <- exportFunPtr mkPreconditionerFunction wrapPreconditioner pref opt + userptr <- getStablePtr opt userdata + withOpt opt $ \o -> parseEnum <$> + nlopt_set_precond_max_objective o objfunptr prefunptr userptr + +{- Working with bounds -} + +foreign import ccall "nlopt.h nlopt_set_lower_bounds" + nlopt_set_lower_bounds :: NloptOpt -> Ptr CDouble -> IO CInt + +foreign import ccall "nlopt.h nlopt_set_lower_bounds1" + nlopt_set_lower_bounds1 :: NloptOpt -> CDouble -> IO CInt + +foreign import ccall "nlopt.h nlopt_get_lower_bounds" + nlopt_get_lower_bounds :: NloptOpt -> Ptr CDouble -> IO CInt + +foreign import ccall "nlopt.h nlopt_set_upper_bounds" + nlopt_set_upper_bounds :: NloptOpt -> Ptr CDouble -> IO CInt + +foreign import ccall "nlopt.h nlopt_set_upper_bounds1" + nlopt_set_upper_bounds1 :: NloptOpt -> CDouble -> IO CInt + +foreign import ccall "nlopt.h nlopt_get_upper_bounds" + nlopt_get_upper_bounds :: NloptOpt -> Ptr CDouble -> IO CInt + +set_lower_bounds :: Opt -> V.Vector Double -> IO Result +set_lower_bounds opt bounds = + withForeignPtr (fst . V.unsafeToForeignPtr0 . V.unsafeCast $ bounds) $ + \bptr -> withOpt opt $ \o -> + parseEnum <$> nlopt_set_lower_bounds o bptr + +set_lower_bounds1 :: Opt -> Double -> IO Result +set_lower_bounds1 opt bound = + withOpt opt $ \o -> + parseEnum <$> nlopt_set_lower_bounds1 o (realToFrac bound) + +get_lower_bounds :: Opt -> IO (V.Vector Double, Result) +get_lower_bounds opt = do + v <- get_dimension opt >>= MV.new . fromIntegral + MV.unsafeWith (MV.unsafeCast v) $ \vptr -> withOpt opt $ \o -> do + result <- parseEnum <$> nlopt_get_lower_bounds o vptr + retv <- V.unsafeFreeze v + return (retv, result) + +set_upper_bounds :: Opt -> V.Vector Double -> IO Result +set_upper_bounds opt bounds = + withForeignPtr (fst . V.unsafeToForeignPtr0 . V.unsafeCast $ bounds) $ + \bptr -> withOpt opt $ \o -> + parseEnum <$> nlopt_set_upper_bounds o bptr + +set_upper_bounds1 :: Opt -> Double -> IO Result +set_upper_bounds1 opt bound = + withOpt opt $ \o -> + parseEnum <$> nlopt_set_upper_bounds1 o (realToFrac bound) + +get_upper_bounds :: Opt -> IO (V.Vector Double, Result) +get_upper_bounds opt = do + v <- get_dimension opt >>= MV.new . fromIntegral + MV.unsafeWith (MV.unsafeCast v) $ \vptr -> withOpt opt $ \o -> do + result <- parseEnum <$> nlopt_get_upper_bounds o vptr + retv <- V.unsafeFreeze v + return (retv, result) + +{- Working with constraints -} + +foreign import ccall "nlopt.h nlopt_remove_inequality_constraints" + nlopt_remove_inequality_constraints :: NloptOpt -> IO CInt + +foreign import ccall "nlopt.h nlopt_add_inequality_constraint" + nlopt_add_inequality_constraint :: NloptOpt -> FunPtr (CFunc a) + -> StablePtr a -> CDouble -> IO CInt + +foreign import ccall "nlopt.h nlopt_add_precond_inequality_constraint" + nlopt_add_precond_inequality_constraint :: NloptOpt -> FunPtr (CFunc a) + -> FunPtr (CPrecond a) -> StablePtr a + -> CDouble -> IO CInt + +foreign import ccall "nlopt.h nlopt_add_inequality_mconstraint" + nlopt_add_inequality_mconstraint :: NloptOpt -> CUInt -> FunPtr (CMFunc a) + -> StablePtr a -> CDouble -> IO CInt + +remove_inequality_constraints :: Opt -> IO Result +remove_inequality_constraints = + useOpt $ fmap parseEnum . nlopt_remove_inequality_constraints + +add_inequality_constraint :: Opt -> ScalarFunction a + -> a -> Double -> IO Result +add_inequality_constraint opt objfun userdata tol = do + objfunptr <- exportFunPtr mkCFunction wrapCFunction objfun opt + userptr <- getStablePtr opt userdata + withOpt opt $ \o -> + parseEnum <$> + nlopt_add_inequality_constraint o objfunptr userptr (realToFrac tol) + +add_precond_inequality_constraint :: Opt -> ScalarFunction a + -> PreconditionerFunction a -> a -> Double + -> IO Result +add_precond_inequality_constraint opt objfun precfun userdata tol = do + objfunptr <- exportFunPtr mkCFunction wrapCFunction objfun opt + precfunptr <- + exportFunPtr mkPreconditionerFunction wrapPreconditioner precfun opt + userptr <- getStablePtr opt userdata + withOpt opt $ \o -> + parseEnum <$> + nlopt_add_precond_inequality_constraint o objfunptr + precfunptr userptr (realToFrac tol) + +add_inequality_mconstraint :: Opt -> Word -> VectorFunction a -> a + -> Double -> IO Result +add_inequality_mconstraint opt constraintsize constrfun userdata tol = do + constrfunptr <- exportFunPtr mkMFunction wrapMFunction constrfun opt + userptr <- getStablePtr opt userdata + withOpt opt $ \o -> + parseEnum <$> + nlopt_add_inequality_mconstraint o (fromIntegral constraintsize) + constrfunptr userptr (realToFrac tol) + +foreign import ccall "nlopt.h nlopt_remove_equality_constraints" + nlopt_remove_equality_constraints :: NloptOpt -> IO CInt + +foreign import ccall "nlopt.h nlopt_add_equality_constraint" + nlopt_add_equality_constraint :: NloptOpt -> FunPtr (CFunc a) + -> StablePtr a -> CDouble -> IO CInt + +foreign import ccall "nlopt.h nlopt_add_precond_equality_constraint" + nlopt_add_precond_equality_constraint :: NloptOpt -> FunPtr (CFunc a) + -> FunPtr (CPrecond a) -> StablePtr a + -> CDouble -> IO CInt + +foreign import ccall "nlopt.h nlopt_add_equality_mconstraint" + nlopt_add_equality_mconstraint :: NloptOpt -> CUInt -> FunPtr (CMFunc a) + -> StablePtr a -> CDouble -> IO CInt + +remove_equality_constraints :: Opt -> IO Result +remove_equality_constraints = + useOpt $ fmap parseEnum . nlopt_remove_equality_constraints + +add_equality_constraint :: Opt -> ScalarFunction a + -> a -> Double -> IO Result +add_equality_constraint opt objfun userdata tol = do + objfunptr <- exportFunPtr mkCFunction wrapCFunction objfun opt + userptr <- getStablePtr opt userdata + withOpt opt $ \o -> + parseEnum <$> + nlopt_add_equality_constraint o objfunptr userptr (realToFrac tol) + +add_precond_equality_constraint :: Opt -> ScalarFunction a + -> PreconditionerFunction a -> a -> Double + -> IO Result +add_precond_equality_constraint opt objfun precfun userdata tol = do + objfunptr <- exportFunPtr mkCFunction wrapCFunction objfun opt + precfunptr <- + exportFunPtr mkPreconditionerFunction wrapPreconditioner precfun opt + userptr <- getStablePtr opt userdata + withOpt opt $ \o -> + parseEnum <$> + nlopt_add_precond_equality_constraint o objfunptr + precfunptr userptr (realToFrac tol) + +add_equality_mconstraint :: Opt -> Word -> VectorFunction a -> a + -> Double -> IO Result +add_equality_mconstraint opt constraintsize constrfun userdata tol = do + constrfunptr <- exportFunPtr mkMFunction wrapMFunction constrfun opt + userptr <- getStablePtr opt userdata + withOpt opt $ \o -> + parseEnum <$> + nlopt_add_equality_mconstraint o (fromIntegral constraintsize) + constrfunptr userptr (realToFrac tol) + +{- Stopping criteria -} + +withInputVector :: (Storable c, Storable a) + => V.Vector c -> (Ptr a -> IO b) -> IO b +withInputVector = withForeignPtr . fst . V.unsafeToForeignPtr0 . V.unsafeCast +withOutputVector :: (Storable c, Storable a) + => V.MVector s c -> (Ptr a -> IO b) -> IO b +withOutputVector = withForeignPtr . fst . MV.unsafeToForeignPtr0 . MV.unsafeCast + +setScalar :: (Enum a, Integral b) => (NloptOpt -> t1 -> IO b) + -> (t -> t1) -> Opt -> t -> IO a +setScalar setter conv opt val = withOpt opt $ \o -> + parseEnum <$> setter o (conv val) + +getScalar :: (NloptOpt -> IO b) -> (b -> a) -> Opt -> IO a +getScalar getter conv = useOpt $ fmap conv . getter + +foreign import ccall "nlopt.h nlopt_set_stopval" + nlopt_set_stopval :: NloptOpt -> CDouble -> IO CInt + +foreign import ccall "nlopt.h nlopt_get_stopval" + nlopt_get_stopval :: NloptOpt -> IO CDouble + +set_stopval :: Opt -> Double -> IO Result +set_stopval = setScalar nlopt_set_stopval realToFrac + +get_stopval :: Opt -> IO Double +get_stopval = getScalar nlopt_get_stopval realToFrac + +foreign import ccall "nlopt.h nlopt_set_ftol_rel" + nlopt_set_ftol_rel :: NloptOpt -> CDouble -> IO CInt + +foreign import ccall "nlopt.h nlopt_get_ftol_rel" + nlopt_get_ftol_rel :: NloptOpt -> IO CDouble + +set_ftol_rel :: Opt -> Double -> IO Result +set_ftol_rel = setScalar nlopt_set_ftol_rel realToFrac + +get_ftol_rel :: Opt -> IO Double +get_ftol_rel = getScalar nlopt_get_ftol_rel realToFrac + +foreign import ccall "nlopt.h nlopt_set_ftol_abs" + nlopt_set_ftol_abs :: NloptOpt -> CDouble -> IO CInt + +foreign import ccall "nlopt.h nlopt_get_ftol_abs" + nlopt_get_ftol_abs :: NloptOpt -> IO CDouble + +set_ftol_abs :: Opt -> Double -> IO Result +set_ftol_abs = setScalar nlopt_set_ftol_abs realToFrac + +get_ftol_abs :: Opt -> IO Double +get_ftol_abs = getScalar nlopt_get_ftol_abs realToFrac + +foreign import ccall "nlopt.h nlopt_set_xtol_rel" + nlopt_set_xtol_rel :: NloptOpt -> CDouble -> IO CInt + +foreign import ccall "nlopt.h nlopt_get_xtol_rel" + nlopt_get_xtol_rel :: NloptOpt -> IO CDouble + +set_xtol_rel :: Opt -> Double -> IO Result +set_xtol_rel = setScalar nlopt_set_xtol_rel realToFrac + +get_xtol_rel :: Opt -> IO Double +get_xtol_rel = getScalar nlopt_get_xtol_rel realToFrac + +foreign import ccall "nlopt.h nlopt_set_xtol_abs1" + nlopt_set_xtol_abs1 :: NloptOpt -> CDouble -> IO CInt + +set_xtol_abs1 :: Opt -> Double -> IO Result +set_xtol_abs1 = setScalar nlopt_set_xtol_abs1 realToFrac + +foreign import ccall "nlopt.h nlopt_set_xtol_abs" + nlopt_set_xtol_abs :: NloptOpt -> Ptr CDouble -> IO CInt + +foreign import ccall "nlopt.h nlopt_get_xtol_abs" + nlopt_get_xtol_abs :: NloptOpt -> Ptr CDouble -> IO CInt + +set_xtol_abs :: Opt -> V.Vector Double -> IO Result +set_xtol_abs opt tolvec = + withInputVector tolvec $ \tolptr -> + withOpt opt $ \o -> parseEnum <$> nlopt_set_xtol_abs o tolptr + +get_xtol_abs :: Opt -> IO (Result, V.Vector Double) +get_xtol_abs opt = do + mutv <- get_dimension opt >>= MV.new . fromIntegral + withOutputVector mutv $ \vecptr -> + withOpt opt $ \o -> do + result <- parseEnum <$> nlopt_get_xtol_abs o vecptr + outvec <- V.unsafeFreeze mutv + return (result, outvec) + +foreign import ccall "nlopt.h nlopt_set_maxeval" + nlopt_set_maxeval :: NloptOpt -> CInt -> IO CInt + +foreign import ccall "nlopt.h nlopt_get_maxeval" + nlopt_get_maxeval :: NloptOpt -> IO CInt + +set_maxeval :: Opt -> Word -> IO Result +set_maxeval = setScalar nlopt_set_maxeval fromIntegral + +get_maxeval :: Opt -> IO Word +get_maxeval = getScalar nlopt_get_maxeval fromIntegral + +foreign import ccall "nlopt.h nlopt_set_maxtime" + nlopt_set_maxtime :: NloptOpt -> CDouble -> IO CInt + +foreign import ccall "nlopt.h nlopt_get_maxtime" + nlopt_get_maxtime :: NloptOpt -> IO CDouble + +set_maxtime :: Opt -> Double -> IO Result +set_maxtime = setScalar nlopt_set_maxtime realToFrac + +get_maxtime :: Opt -> IO Double +get_maxtime = getScalar nlopt_get_maxtime realToFrac + +foreign import ccall "nlopt.h nlopt_force_stop" + nlopt_force_stop :: NloptOpt -> IO CInt + +force_stop :: Opt -> IO Result +force_stop = useOpt $ fmap parseEnum . nlopt_force_stop + +foreign import ccall "nlopt.h nlopt_set_force_stop" + nlopt_set_force_stop :: NloptOpt -> CInt -> IO CInt + +foreign import ccall "nlopt.h nlopt_get_force_stop" + nlopt_get_force_stop :: NloptOpt -> IO CInt + +set_force_stop :: Opt -> Word -> IO Result +set_force_stop = setScalar nlopt_set_force_stop fromIntegral + +get_force_stop :: Opt -> IO Word +get_force_stop = getScalar nlopt_get_force_stop fromIntegral + +{- Algorithm-specific configuration -} + +foreign import ccall "nlopt.h nlopt_set_local_optimizer" + nlopt_set_local_optimizer :: NloptOpt -> NloptOpt -> IO CInt + +set_local_optimizer :: Opt -- ^ Primary optimizer + -> Opt -- ^ Subsidiary (local) optimizer + -> IO Result +set_local_optimizer p s = + withOpt p $ \primary -> withOpt s $ \secondary -> + parseEnum <$> nlopt_set_local_optimizer primary secondary + +foreign import ccall "nlopt.h nlopt_set_population" + nlopt_set_population :: NloptOpt -> Word -> IO CInt + +foreign import ccall "nlopt.h nlopt_get_population" + nlopt_get_population :: NloptOpt -> IO Word + +set_population :: Opt -> Word -> IO Result +set_population = setScalar nlopt_set_population fromIntegral + +get_population :: Opt -> IO Word +get_population = getScalar nlopt_get_population fromIntegral + +foreign import ccall "nlopt.h nlopt_set_vector_storage" + nlopt_set_vector_storage :: NloptOpt -> Word -> IO CInt + +foreign import ccall "nlopt.h nlopt_get_vector_storage" + nlopt_get_vector_storage :: NloptOpt -> IO Word + +set_vector_storage :: Opt -> Word -> IO Result +set_vector_storage = setScalar nlopt_set_vector_storage fromIntegral + +get_vector_storage :: Opt -> IO Word +get_vector_storage = getScalar nlopt_get_vector_storage fromIntegral + +foreign import ccall "nlopt.h nlopt_set_default_initial_step" + nlopt_set_default_initial_step :: NloptOpt -> Ptr CDouble -> IO CInt + +foreign import ccall "nlopt.h nlopt_set_initial_step" + nlopt_set_initial_step :: NloptOpt -> Ptr CDouble -> IO CInt + +foreign import ccall "nlopt.h nlopt_set_initial_step1" + nlopt_set_initial_step1 :: NloptOpt -> CDouble -> IO CInt + +foreign import ccall "nlopt.h nlopt_get_initial_step" + nlopt_get_initial_step :: NloptOpt -> Ptr CDouble -> Ptr CDouble -> IO CInt + +set_default_initial_step :: Opt -> V.Vector Double -> IO Result +set_default_initial_step opt stepvec = + withInputVector stepvec $ \stepptr -> + withOpt opt $ \o -> parseEnum <$> nlopt_set_default_initial_step o stepptr + +set_initial_step :: Opt -> V.Vector Double -> IO Result +set_initial_step opt stepvec = + withInputVector stepvec $ \stepptr -> + withOpt opt $ \o -> parseEnum <$> nlopt_set_initial_step o stepptr + +set_initial_step1 :: Opt -> Double -> IO Result +set_initial_step1 = setScalar nlopt_set_initial_step1 realToFrac + +get_initial_step :: Opt -> V.Vector Double -> IO (Result, V.Vector Double) +get_initial_step opt xvec = do + mutv <- get_dimension opt >>= MV.new . fromIntegral + withOutputVector mutv $ \outptr -> + withInputVector xvec $ \inptr -> + withOpt opt $ \o -> do + result <- parseEnum <$> nlopt_get_initial_step o inptr outptr + outvec <- V.unsafeFreeze mutv + return (result, outvec) diff --git a/srtree.cabal b/srtree.cabal --- a/srtree.cabal +++ b/srtree.cabal @@ -5,7 +5,7 @@ -- see: https://github.com/sol/hpack name: srtree -version: 2.0.0.0 +version: 2.0.0.1 synopsis: A general library to work with Symbolic Regression expression trees. description: A Symbolic Regression Tree data structure to work with mathematical expressions with support to first order derivative and simplification; category: Math, Data, Data Structures @@ -49,6 +49,7 @@ Data.SRTree.Print Data.SRTree.Random Data.SRTree.Recursion + Numeric.Optimization.NLOPT.Bindings Text.ParseSR Text.ParseSR.IO other-modules: @@ -56,11 +57,15 @@ hs-source-dirs: src ghc-options: -fwarn-incomplete-patterns + extra-lib-dirs: + /usr/local/lib + extra-libraries: + nlopt build-depends: attoparsec >=0.14.4 && <0.15 , attoparsec-expr >=0.1.1.2 && <0.2 , base >=4.16 && <5 - , bytestring ==0.11.* + , bytestring >=0.11 && <0.13 , containers >=0.6.7 && <0.8 , dlist ==1.0.* , exceptions >=0.10.7 && <0.11 @@ -71,7 +76,6 @@ , list-shuffle >=1.0.0.1 && <1.1 , massiv >=1.0.4.0 && <1.1 , mtl >=2.2 && <2.4 - , nlopt-haskell >=0.1.3.0 && <0.2 , random ==1.2.* , split >=0.2.5 && <0.3 , statistics >=0.16.2.1 && <0.17 @@ -93,7 +97,7 @@ attoparsec >=0.14.4 && <0.15 , attoparsec-expr >=0.1.1.2 && <0.2 , base >=4.16 && <5 - , bytestring ==0.11.* + , bytestring >=0.11 && <0.13 , containers >=0.6.7 && <0.8 , dlist ==1.0.* , exceptions >=0.10.7 && <0.11 @@ -104,7 +108,6 @@ , list-shuffle >=1.0.0.1 && <1.1 , massiv >=1.0.4.0 && <1.1 , mtl >=2.2 && <2.4 - , nlopt-haskell >=0.1.3.0 && <0.2 , optparse-applicative >=0.17 && <0.19 , random ==1.2.* , split >=0.2.5 && <0.3 @@ -127,7 +130,7 @@ attoparsec >=0.14.4 && <0.15 , attoparsec-expr >=0.1.1.2 && <0.2 , base >=4.16 && <5 - , bytestring ==0.11.* + , bytestring >=0.11 && <0.13 , containers >=0.6.7 && <0.8 , dlist ==1.0.* , exceptions >=0.10.7 && <0.11 @@ -138,7 +141,6 @@ , list-shuffle >=1.0.0.1 && <1.1 , massiv >=1.0.4.0 && <1.1 , mtl >=2.2 && <2.4 - , nlopt-haskell >=0.1.3.0 && <0.2 , random ==1.2.* , split >=0.2.5 && <0.3 , srtree @@ -160,7 +162,7 @@ attoparsec >=0.14.4 && <0.15 , attoparsec-expr >=0.1.1.2 && <0.2 , base >=4.16 && <5 - , bytestring ==0.11.* + , bytestring >=0.11 && <0.13 , containers >=0.6.7 && <0.8 , dlist ==1.0.* , exceptions >=0.10.7 && <0.11 @@ -171,7 +173,6 @@ , list-shuffle >=1.0.0.1 && <1.1 , massiv >=1.0.4.0 && <1.1 , mtl >=2.2 && <2.4 - , nlopt-haskell >=0.1.3.0 && <0.2 , optparse-applicative >=0.17 && <0.19 , random ==1.2.* , split >=0.2.5 && <0.3 @@ -194,7 +195,7 @@ attoparsec >=0.14.4 && <0.15 , attoparsec-expr >=0.1.1.2 && <0.2 , base >=4.16 && <5 - , bytestring ==0.11.* + , bytestring >=0.11 && <0.13 , containers >=0.6.7 && <0.8 , dlist ==1.0.* , exceptions >=0.10.7 && <0.11 @@ -205,7 +206,6 @@ , list-shuffle >=1.0.0.1 && <1.1 , massiv >=1.0.4.0 && <1.1 , mtl >=2.2 && <2.4 - , nlopt-haskell >=0.1.3.0 && <0.2 , optparse-applicative >=0.17 && <0.19 , random ==1.2.* , split >=0.2.5 && <0.3 @@ -231,7 +231,7 @@ attoparsec >=0.14.4 && <0.15 , attoparsec-expr >=0.1.1.2 && <0.2 , base >=4.16 && <5 - , bytestring ==0.11.* + , bytestring >=0.11 && <0.13 , containers >=0.6.7 && <0.8 , dlist ==1.0.* , exceptions >=0.10.7 && <0.11 @@ -242,7 +242,6 @@ , list-shuffle >=1.0.0.1 && <1.1 , massiv >=1.0.4.0 && <1.1 , mtl >=2.2 && <2.4 - , nlopt-haskell >=0.1.3.0 && <0.2 , normaldistribution >=1.1.0.3 && <1.2 , optparse-applicative >=0.17 && <0.19 , random ==1.2.* @@ -268,7 +267,7 @@ attoparsec >=0.14.4 && <0.15 , attoparsec-expr >=0.1.1.2 && <0.2 , base >=4.16 && <5 - , bytestring ==0.11.* + , bytestring >=0.11 && <0.13 , containers >=0.6.7 && <0.8 , dlist ==1.0.* , exceptions >=0.10.7 && <0.11 @@ -279,7 +278,6 @@ , list-shuffle >=1.0.0.1 && <1.1 , massiv >=1.0.4.0 && <1.1 , mtl >=2.2 && <2.4 - , nlopt-haskell >=0.1.3.0 && <0.2 , optparse-applicative >=0.17 && <0.19 , random ==1.2.* , split >=0.2.5 && <0.3 @@ -305,7 +303,7 @@ , attoparsec >=0.14.4 && <0.15 , attoparsec-expr >=0.1.1.2 && <0.2 , base >=4.16 && <5 - , bytestring ==0.11.* + , bytestring >=0.11 && <0.13 , containers >=0.6.7 && <0.8 , dlist ==1.0.* , exceptions >=0.10.7 && <0.11 @@ -316,7 +314,6 @@ , list-shuffle >=1.0.0.1 && <1.1 , massiv >=1.0.4.0 && <1.1 , mtl >=2.2 && <2.4 - , nlopt-haskell >=0.1.3.0 && <0.2 , random ==1.2.* , split >=0.2.5 && <0.3 , srtree