-- SPDX-License-Identifier: MPL-2.0

{- | Dense bounded reference constrained optimization. Inequalities use g(x)<=0.
LP uses nonnegative variables and Ax<=b, supports negative right-hand sides
through phase I, and reports infeasibility/unboundedness. Revised simplex has
an explicit feasible-basis contract. Convex QP enumerates active sets up to a
caller budget; it is intended for small problems, not a sparse QP backend.
-}
module HMath.Optimization.Constrained (
    Constraint (..),
    feasibility,
    PenaltyMethod (..),
    constrainedMinimize,
    LPResult (..),
    linearProgramming,
    revisedSimplex,
    quadraticProgramming,
    slsqp,
    cobyla,
    weightedObjective,
    paretoFront,
) where

import qualified Data.Foldable as Foldable (foldl')
import Data.List (nub, sort)
import HMath.Error
import HMath.Optimization.Common
import HMath.Optimization.Local (QuasiNewton (BFGS), minimizeQuasiNewton)

{- | Scalar Equality or LessEqual constraint over the parameter vector; inequality functions
use the solver's \<=0 convention.
-}
data Constraint = Equality Objective Gradient | LessEqual Objective Gradient

{- | Constrained-optimization strategy selecting quadratic penalty,log barrier,augmented
Lagrangian,or interior-point treatment.
-}
data PenaltyMethod = QuadraticPenalty | LogBarrier | InteriorPoint | AugmentedLagrangian
    deriving (PenaltyMethod -> PenaltyMethod -> Bool
(PenaltyMethod -> PenaltyMethod -> Bool)
-> (PenaltyMethod -> PenaltyMethod -> Bool) -> Eq PenaltyMethod
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: PenaltyMethod -> PenaltyMethod -> Bool
== :: PenaltyMethod -> PenaltyMethod -> Bool
$c/= :: PenaltyMethod -> PenaltyMethod -> Bool
/= :: PenaltyMethod -> PenaltyMethod -> Bool
Eq, Int -> PenaltyMethod -> ShowS
[PenaltyMethod] -> ShowS
PenaltyMethod -> String
(Int -> PenaltyMethod -> ShowS)
-> (PenaltyMethod -> String)
-> ([PenaltyMethod] -> ShowS)
-> Show PenaltyMethod
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> PenaltyMethod -> ShowS
showsPrec :: Int -> PenaltyMethod -> ShowS
$cshow :: PenaltyMethod -> String
show :: PenaltyMethod -> String
$cshowList :: [PenaltyMethod] -> ShowS
showList :: [PenaltyMethod] -> ShowS
Show)

constraintValues :: [Constraint] -> [Double] -> Either MathError [(Bool, Double, [Double])]
constraintValues :: [Constraint]
-> [Double] -> Either MathError [(Bool, Double, [Double])]
constraintValues [Constraint]
cs [Double]
x = (Constraint -> Either MathError (Bool, Double, [Double]))
-> [Constraint] -> Either MathError [(Bool, Double, [Double])]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
forall (m :: * -> *) a b. Monad m => (a -> m b) -> [a] -> m [b]
mapM Constraint -> Either MathError (Bool, Double, [Double])
value [Constraint]
cs
  where
    value :: Constraint -> Either MathError (Bool, Double, [Double])
value Constraint
c = case Constraint
c of
        Equality Objective
f Gradient
g -> do y <- Objective -> Objective
checkedValue Objective
f [Double]
x; gs <- checkedGradient g x; Right (True, y, gs)
        LessEqual Objective
f Gradient
g -> do y <- Objective -> Objective
checkedValue Objective
f [Double]
x; gs <- checkedGradient g x; Right (False, y, gs)

-- | Return maximum absolute equality or positive inequality violation, with checked values and gradients.
feasibility :: [Constraint] -> [Double] -> Either MathError Double
feasibility :: [Constraint] -> Objective
feasibility [Constraint]
cs [Double]
x = do
    values <- [Constraint]
-> [Double] -> Either MathError [(Bool, Double, [Double])]
constraintValues [Constraint]
cs [Double]
x
    Right (maximum (0 : [if eq then abs y else max 0 y | (eq, y, _) <- values]))

-- | Outer continuation or multiplier updates around the named constrained objective. Every inner callback charges its nested objective/constraint/gradient calls; barriers require strict feasibility.
constrainedMinimize ::
    PenaltyMethod ->
    OptimizeOptions ->
    Objective ->
    Gradient ->
    [Constraint] ->
    [Double] ->
    Either MathError OptimizeResult
constrainedMinimize :: PenaltyMethod
-> OptimizeOptions
-> Objective
-> Gradient
-> [Constraint]
-> [Double]
-> Either MathError OptimizeResult
constrainedMinimize PenaltyMethod
method OptimizeOptions
o Objective
f Gradient
g [Constraint]
constraints [Double]
x0 = do
    OptimizeOptions -> Either MathError ()
checkedOptions OptimizeOptions
o
    if PenaltyMethod
method PenaltyMethod -> PenaltyMethod -> Bool
forall a. Eq a => a -> a -> Bool
== PenaltyMethod
LogBarrier Bool -> Bool -> Bool
|| PenaltyMethod
method PenaltyMethod -> PenaltyMethod -> Bool
forall a. Eq a => a -> a -> Bool
== PenaltyMethod
InteriorPoint
        then do
            vals <- [Constraint]
-> [Double] -> Either MathError [(Bool, Double, [Double])]
constraintValues [Constraint]
constraints [Double]
x0
            if any (\(Bool
eq, Double
v, [Double]
_) -> Bool
eq Bool -> Bool -> Bool
|| Double
v Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0) vals
                then Left (DomainError "log barrier requires strictly feasible inequalities only")
                else Right ()
        else () -> Either MathError ()
forall a b. b -> Either a b
Right ()
    Int
-> Int
-> Double
-> [Double]
-> [Double]
-> Either MathError OptimizeResult
loop Int
0 Int
0 Double
1 (Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate ([Constraint] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Constraint]
constraints) Double
0) [Double]
x0
  where
    loop :: Int
-> Int
-> Double
-> [Double]
-> [Double]
-> Either MathError OptimizeResult
loop Int
outerIter Int
cost Double
strength [Double]
multipliers [Double]
x
        | Int
outerIter Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
40 (OptimizeOptions -> Int
optMaxIterations OptimizeOptions
o) = MathError -> Either MathError OptimizeResult
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence (PenaltyMethod -> String
forall a. Show a => a -> String
show PenaltyMethod
method) Int
outerIter)
        | Bool
otherwise = do
            let evaluate :: Objective
evaluate [Double]
z = do
                    base <- Objective -> Objective
checkedValue Objective
f [Double]
z
                    vals <- constraintValues constraints z
                    terms <-
                        mapM
                            ( \((Bool
eq, Double
v, [Double]
_), Double
lambda) -> case PenaltyMethod
method of
                                PenaltyMethod
LogBarrier -> if Double
v Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0 then MathError -> Either MathError Double
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"outside barrier") else Double -> Either MathError Double
forall a b. b -> Either a b
Right (-Double -> Double
forall a. Floating a => a -> a
log (-Double
v) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
strength)
                                PenaltyMethod
InteriorPoint -> if Double
v Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0 then MathError -> Either MathError Double
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"outside interior") else Double -> Either MathError Double
forall a b. b -> Either a b
Right (-Double -> Double
forall a. Floating a => a -> a
log (-Double
v) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
strength)
                                PenaltyMethod
QuadraticPenalty -> Double -> Either MathError Double
forall a b. b -> Either a b
Right (Double
strength Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (if Bool
eq then Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v else Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
0 Double
v Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int)))
                                PenaltyMethod
AugmentedLagrangian ->
                                    Double -> Either MathError Double
forall a b. b -> Either a b
Right
                                        ( if Bool
eq
                                            then Double
lambda Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
strength Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2
                                            else ((Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
0 (Double
lambda Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
strength Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v)) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
lambda Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
lambda) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
strength)
                                        )
                            )
                            (zip vals multipliers)
                    Right (base + sum terms)
                derivative :: Gradient
derivative [Double]
z = do
                    base <- Gradient -> Gradient
checkedGradient Gradient
g [Double]
z
                    vals <- constraintValues constraints z
                    terms <-
                        mapM
                            ( \((Bool
eq, Double
v, [Double]
gs), Double
lambda) -> do
                                coeff <- case PenaltyMethod
method of
                                    PenaltyMethod
LogBarrier -> if Double
v Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0 then MathError -> Either MathError Double
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"outside barrier") else Double -> Either MathError Double
forall a b. b -> Either a b
Right (-Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
strength Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v))
                                    PenaltyMethod
InteriorPoint -> if Double
v Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0 then MathError -> Either MathError Double
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"outside interior") else Double -> Either MathError Double
forall a b. b -> Either a b
Right (-Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
strength Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v))
                                    PenaltyMethod
QuadraticPenalty -> Double -> Either MathError Double
forall a b. b -> Either a b
Right (Double
strength Double -> Double -> Double
forall a. Num a => a -> a -> a
* (if Bool
eq then Double
v else Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
0 Double
v))
                                    PenaltyMethod
AugmentedLagrangian -> Double -> Either MathError Double
forall a b. b -> Either a b
Right (if Bool
eq then Double
lambda Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
strength Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v else Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
0 (Double
lambda Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
strength Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v))
                                Right (scale coeff gs)
                            )
                            (zip vals multipliers)
                    Right (Foldable.foldl' add base terms)
                inner :: OptimizeOptions
inner =
                    OptimizeOptions
o
                        { optMaxIterations = max 10 (optMaxIterations o `div` 2)
                        , optMaxEvaluations = (optMaxEvaluations o - cost) `div` (1 + 2 * length constraints)
                        , optTolerance = max (optTolerance o / 10) (1 / (strength * strength) * 1e-4)
                        }
            -- A feasible backtracking gradient subsolver avoids evaluating an
            -- unconstrained line search outside the open barrier domain.
            r <-
                if PenaltyMethod
method PenaltyMethod -> PenaltyMethod -> Bool
forall a. Eq a => a -> a -> Bool
== PenaltyMethod
LogBarrier Bool -> Bool -> Bool
|| PenaltyMethod
method PenaltyMethod -> PenaltyMethod -> Bool
forall a. Eq a => a -> a -> Bool
== PenaltyMethod
InteriorPoint
                    then OptimizeOptions
-> Objective
-> Gradient
-> [Double]
-> Either MathError OptimizeResult
barrierInner OptimizeOptions
inner Objective
evaluate Gradient
derivative [Double]
x
                    else QuasiNewton
-> OptimizeOptions
-> Objective
-> Gradient
-> [Double]
-> Either MathError OptimizeResult
minimizeQuasiNewton QuasiNewton
BFGS OptimizeOptions
inner Objective
evaluate Gradient
derivative [Double]
x
            let xn = OptimizeResult -> [Double]
optPoint OptimizeResult
r; cost' = Int
cost Int -> Int -> Int
forall a. Num a => a -> a -> a
+ OptimizeResult -> Int
optEvaluations OptimizeResult
r Int -> Int -> Int
forall a. Num a => a -> a -> a
* (Int
1 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* [Constraint] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Constraint]
constraints)
            checkBudget o (cost' + 2 + 2 * length constraints)
            vals <- constraintValues constraints xn
            gs <- checkedGradient g xn
            let nextLambda =
                    [ let base :: Double
base = if PenaltyMethod
method PenaltyMethod -> PenaltyMethod -> Bool
forall a. Eq a => a -> a -> Bool
== PenaltyMethod
QuadraticPenalty then Double
0 else Double
l
                       in if Bool
eq then Double
base Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
strength Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v else Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
0 (Double
base Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
strength Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v)
                    | ((Bool
eq, Double
v, [Double]
_), Double
l) <- [(Bool, Double, [Double])]
-> [Double] -> [((Bool, Double, [Double]), Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [(Bool, Double, [Double])]
vals [Double]
multipliers
                    ]
                lagGradient = ([Double] -> [Double] -> [Double])
-> [Double] -> [[Double]] -> [Double]
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
Foldable.foldl' [Double] -> [Double] -> [Double]
add [Double]
gs [Double -> [Double] -> [Double]
scale Double
l [Double]
d | ((Bool
_, Double
_, [Double]
d), Double
l) <- [(Bool, Double, [Double])]
-> [Double] -> [((Bool, Double, [Double]), Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [(Bool, Double, [Double])]
vals [Double]
nextLambda]
                violation = [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum (Double
0 Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: [if Bool
eq then Double -> Double
forall a. Num a => a -> a
abs Double
v else Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
0 Double
v | (Bool
eq, Double
v, [Double]
_) <- [(Bool, Double, [Double])]
vals])
                dual =
                    if PenaltyMethod
method PenaltyMethod -> PenaltyMethod -> Bool
forall a. Eq a => a -> a -> Bool
== PenaltyMethod
LogBarrier Bool -> Bool -> Bool
|| PenaltyMethod
method PenaltyMethod -> PenaltyMethod -> Bool
forall a. Eq a => a -> a -> Bool
== PenaltyMethod
InteriorPoint
                        then Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Constraint] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Constraint]
constraints) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
strength
                        else [Double] -> Double
norm [Double]
lagGradient
            value <- checkedValue f xn
            checkBudget o (cost' + 2 + 2 * length constraints)
            if violation <= optTolerance o && dual <= sqrt (optTolerance o)
                then
                    Right
                        ( OptimizeResult
                            xn
                            value
                            (max violation dual)
                            (outerIter + 1)
                            (cost' + 2 + 2 * length constraints)
                            True
                        )
                else loop (outerIter + 1) (cost' + 2 + 2 * length constraints) (strength * 10) nextLambda xn
    barrierInner :: OptimizeOptions
-> Objective
-> Gradient
-> [Double]
-> Either MathError OptimizeResult
barrierInner OptimizeOptions
opts Objective
obj Gradient
grad [Double]
start = Int -> Int -> [Double] -> Either MathError OptimizeResult
go Int
0 Int
0 [Double]
start
      where
        go :: Int -> Int -> [Double] -> Either MathError OptimizeResult
go Int
k Int
count [Double]
x = do
            OptimizeOptions -> Int -> Either MathError ()
checkBudget OptimizeOptions
opts (Int
count Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2)
            y <- Objective -> Objective
checkedValue Objective
obj [Double]
x
            gs <- checkedGradient grad x
            if norm gs <= optTolerance opts
                then Right (OptimizeResult x y (norm gs) k (count + 2) True)
                else
                    if k >= optMaxIterations opts
                        then Right (OptimizeResult x y (norm gs) k (count + 2) False)
                        else do
                            (next, used) <- backtrack (optMaxEvaluations opts - count - 2) 0 1 y gs x
                            go (k + 1) (count + 2 + used) next
        backtrack :: Int
-> Int
-> Double
-> Double
-> [Double]
-> [Double]
-> Either MathError ([Double], Int)
backtrack Int
remaining Int
k Double
rate Double
y [Double]
gs [Double]
x
            | Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
remaining = MathError -> Either MathError ([Double], Int)
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"barrier evaluation budget" Int
k)
            | Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
50 = MathError -> Either MathError ([Double], Int)
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"barrier feasible line search" Int
k)
            | Bool
otherwise =
                let next :: [Double]
next = [Double] -> [Double] -> [Double]
sub [Double]
x (Double -> [Double] -> [Double]
scale Double
rate [Double]
gs)
                 in case Objective -> Objective
checkedValue Objective
obj [Double]
next of
                        Right Double
yn | Double
yn Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1e-4 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
rate Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double] -> [Double] -> Double
dot [Double]
gs [Double]
gs -> ([Double], Int) -> Either MathError ([Double], Int)
forall a b. b -> Either a b
Right ([Double]
next, Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)
                        Left (DomainError String
_) -> Int
-> Int
-> Double
-> Double
-> [Double]
-> [Double]
-> Either MathError ([Double], Int)
backtrack Int
remaining (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) (Double
rate Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2) Double
y [Double]
gs [Double]
x
                        Left MathError
e -> MathError -> Either MathError ([Double], Int)
forall a b. a -> Either a b
Left MathError
e
                        Either MathError Double
_ -> Int
-> Int
-> Double
-> Double
-> [Double]
-> [Double]
-> Either MathError ([Double], Int)
backtrack Int
remaining (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) (Double
rate Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2) Double
y [Double]
gs [Double]
x

{- | Linear-program primal solution,objective value,retained basis indices,and iteration
count.
-}
data LPResult = LPResult
    { LPResult -> [Double]
lpPoint :: [Double]
    -- ^ Returned primal variable vector of the linear program.
    , LPResult -> Double
lpValue :: Double
    -- ^ Objective value at the returned LP point.
    , LPResult -> [Int]
lpBasis :: [Int]
    -- ^ Zero-based basis indices retained by the simplex solution.
    , LPResult -> Int
lpIterations :: Int
    -- ^ Number of simplex iterations performed.
    }
    deriving (LPResult -> LPResult -> Bool
(LPResult -> LPResult -> Bool)
-> (LPResult -> LPResult -> Bool) -> Eq LPResult
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: LPResult -> LPResult -> Bool
== :: LPResult -> LPResult -> Bool
$c/= :: LPResult -> LPResult -> Bool
/= :: LPResult -> LPResult -> Bool
Eq, Int -> LPResult -> ShowS
[LPResult] -> ShowS
LPResult -> String
(Int -> LPResult -> ShowS)
-> (LPResult -> String) -> ([LPResult] -> ShowS) -> Show LPResult
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> LPResult -> ShowS
showsPrec :: Int -> LPResult -> ShowS
$cshow :: LPResult -> String
show :: LPResult -> String
$cshowList :: [LPResult] -> ShowS
showList :: [LPResult] -> ShowS
Show)

{- | Primal two-phase simplex with Bland entering/leaving rules. Tolerance is
used for feasibility and reduced costs, not to round returned variable data.
-}
linearProgramming ::
    Int -> Double -> [Double] -> [[Double]] -> [Double] -> Either MathError LPResult
linearProgramming :: Int
-> Double
-> [Double]
-> [[Double]]
-> [Double]
-> Either MathError LPResult
linearProgramming Int
budget Double
tolerance [Double]
c [[Double]]
a [Double]
b = do
    if Int
budget Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
1 Bool -> Bool -> Bool
|| Double
tolerance Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
0
        then MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"LP tolerance/budget")
        else [Double] -> Either MathError ()
checkedPoint (Double
tolerance Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: [Double]
c [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [[Double]] -> [Double]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [[Double]]
a [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double]
b)
    let n :: Int
n = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
c
        m :: Int
m = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
b
        negative :: [Int]
negative = [Int
i | (Int
i, Double
bi) <- [Int] -> [Double] -> [(Int, Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 :: Int ..] [Double]
b, Double
bi Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0]
        q :: Int
q = [Int] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Int]
negative
        total :: Int
total = Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
q
    if [[Double]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Double]]
a Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
m Bool -> Bool -> Bool
|| ([Double] -> Bool) -> [[Double]] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any ((Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
n) (Int -> Bool) -> ([Double] -> Int) -> [Double] -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length) [[Double]]
a
        then MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"LP constraint shape")
        else () -> Either MathError ()
forall a b. b -> Either a b
Right ()
    let rows :: [[Double]]
rows =
            [ let sign :: Double
sign = if Double
bi Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0 then -Double
1 else Double
1
                  slack :: [Double]
slack = [if Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
i then Double
sign else Double
0 | Int
j <- [Int
0 .. Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
                  artificial :: [Double]
artificial = [if Int -> Maybe Int
forall a. a -> Maybe a
Just Int
k Maybe Int -> Maybe Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int -> [(Int, Int)] -> Maybe Int
forall a b. Eq a => a -> [(a, b)] -> Maybe b
lookup Int
i ([Int] -> [Int] -> [(Int, Int)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int]
negative [Int
0 ..]) then Double
1 else Double
0 | Int
k <- [Int
0 .. Int
q Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
               in Double -> [Double] -> [Double]
scale Double
sign [Double]
row [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double]
slack [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double]
artificial [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double -> Double
forall a. Num a => a -> a
abs Double
bi]
            | (Int
i, ([Double]
row, Double
bi)) <- [Int] -> [([Double], Double)] -> [(Int, ([Double], Double))]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 ..] ([[Double]] -> [Double] -> [([Double], Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [[Double]]
a [Double]
b)
            ]
        basis :: [Int]
basis = [if Double
bi Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0 then Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
+ [Int] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ((Int -> Bool) -> [Int] -> [Int]
forall a. (a -> Bool) -> [a] -> [a]
filter (Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
i) [Int]
negative) else Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
i | (Int
i, Double
bi) <- [Int] -> [Double] -> [(Int, Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 ..] [Double]
b]
        phaseCost :: [Double]
phaseCost = Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
m) Double
0 [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate Int
q Double
1
    (r1, b1, k1) <- Int
-> Double
-> [Double]
-> [[Double]]
-> [Int]
-> Int
-> Either MathError ([[Double]], [Int], Int)
simplexTable Int
budget Double
tolerance [Double]
phaseCost [[Double]]
rows [Int]
basis Int
0
    let phaseValue = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [[Double]
phaseCost [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j Double -> Double -> Double
forall a. Num a => a -> a -> a
* ([Double]
row [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
total) | ([Double]
row, Int
j) <- [[Double]] -> [Int] -> [([Double], Int)]
forall a b. [a] -> [b] -> [(a, b)]
zip [[Double]]
r1 [Int]
b1]
    if phaseValue > tolerance then Left (DomainError "infeasible linear program") else Right ()
    (r2, b2, cleanupIterations) <- removeArtificial total (n + m) r1 b1 k1
    let retained = ([Double] -> [Double]) -> [[Double]] -> [[Double]]
forall a b. (a -> b) -> [a] -> [b]
map (\[Double]
row -> Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
take (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
m) [Double]
row [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [[Double]
row [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
total]) [[Double]]
r2
    (r3, b3, k3) <- simplexTable budget tolerance (c ++ replicate m 0) retained b2 cleanupIterations
    let point = [[Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [[Double]
row [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
m) | ([Double]
row, Int
j) <- [[Double]] -> [Int] -> [([Double], Int)]
forall a b. [a] -> [b] -> [(a, b)]
zip [[Double]]
r3 [Int]
b3, Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
i] | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
    checkedPoint point
    Right (LPResult point (dot c point) b3 k3)
  where
    removeArtificial :: Int
-> Int
-> [[Double]]
-> [Int]
-> Int
-> Either MathError ([[Double]], [Int], Int)
removeArtificial Int
total Int
cutoff [[Double]]
rows [Int]
basics Int
pivots = case [(Int
i, Int
basic) | (Int
i, Int
basic) <- [Int] -> [Int] -> [(Int, Int)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 :: Int ..] [Int]
basics, Int
basic Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
cutoff] of
        [] -> ([[Double]], [Int], Int)
-> Either MathError ([[Double]], [Int], Int)
forall a b. b -> Either a b
Right ([[Double]]
rows, [Int]
basics, Int
pivots)
        (Int
index, Int
_) : [(Int, Int)]
_ ->
            let row :: [Double]
row = [[Double]]
rows [[Double]] -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
index
                available :: [Int]
available = [Int
j | Int
j <- [Int
0 .. Int
cutoff Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1], Double -> Double
forall a. Num a => a -> a
abs ([Double]
row [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
tolerance, Int
j Int -> [Int] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`notElem` [Int]
basics]
             in case [Int]
available of
                    [] ->
                        if Double -> Double
forall a. Num a => a -> a
abs ([Double]
row [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
total) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
tolerance
                            then
                                Int
-> Int
-> [[Double]]
-> [Int]
-> Int
-> Either MathError ([[Double]], [Int], Int)
removeArtificial
                                    Int
total
                                    Int
cutoff
                                    [[Double]
r | (Int
i, [Double]
r) <- [Int] -> [[Double]] -> [(Int, [Double])]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 ..] [[Double]]
rows, Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
index]
                                    [Int
b' | (Int
i, Int
b') <- [Int] -> [Int] -> [(Int, Int)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 ..] [Int]
basics, Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
index]
                                    Int
pivots
                            else MathError -> Either MathError ([[Double]], [Int], Int)
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"infeasible artificial row")
                    Int
col : [Int]
_
                        | Int
pivots Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
budget -> MathError -> Either MathError ([[Double]], [Int], Int)
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"simplex artificial cleanup pivot budget" Int
pivots)
                        | Bool
otherwise ->
                            let normalized :: [Double]
normalized = Double -> [Double] -> [Double]
scale (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ ([Double]
row [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
col)) [Double]
row
                                next :: [[Double]]
next =
                                    [if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
index then [Double]
normalized else [Double] -> [Double] -> [Double]
sub [Double]
r (Double -> [Double] -> [Double]
scale ([Double]
r [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
col) [Double]
normalized) | (Int
i, [Double]
r) <- [Int] -> [[Double]] -> [(Int, [Double])]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 ..] [[Double]]
rows]
                                basis' :: [Int]
basis' = [if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
index then Int
col else Int
b' | (Int
i, Int
b') <- [Int] -> [Int] -> [(Int, Int)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 ..] [Int]
basics]
                             in Int
-> Int
-> [[Double]]
-> [Int]
-> Int
-> Either MathError ([[Double]], [Int], Int)
removeArtificial Int
total Int
cutoff [[Double]]
next [Int]
basis' (Int
pivots Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)

simplexTable ::
    Int -> Double -> [Double] -> [[Double]] -> [Int] -> Int -> Either MathError ([[Double]], [Int], Int)
simplexTable :: Int
-> Double
-> [Double]
-> [[Double]]
-> [Int]
-> Int
-> Either MathError ([[Double]], [Int], Int)
simplexTable Int
budget Double
tol [Double]
costs = [[Double]]
-> [Int] -> Int -> Either MathError ([[Double]], [Int], Int)
go
  where
    n :: Int
n = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
costs
    go :: [[Double]]
-> [Int] -> Int -> Either MathError ([[Double]], [Int], Int)
go [[Double]]
rows [Int]
basis Int
k
        | Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
budget = MathError -> Either MathError ([[Double]], [Int], Int)
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"simplex pivots" Int
k)
        | Bool
otherwise =
            let reduced :: [Double]
reduced = [[Double]
costs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [[Double]
costs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
bj Double -> Double -> Double
forall a. Num a => a -> a -> a
* ([Double]
row [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j) | ([Double]
row, Int
bj) <- [[Double]] -> [Int] -> [([Double], Int)]
forall a b. [a] -> [b] -> [(a, b)]
zip [[Double]]
rows [Int]
basis] | Int
j <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
                entering :: [Int]
entering = [Int
j | (Int
j, Double
v) <- [Int] -> [Double] -> [(Int, Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 ..] [Double]
reduced, Double
v Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< -Double
tol, Int
j Int -> [Int] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`notElem` [Int]
basis]
             in case [Int]
entering of
                    [] -> ([[Double]], [Int], Int)
-> Either MathError ([[Double]], [Int], Int)
forall a b. b -> Either a b
Right ([[Double]]
rows, [Int]
basis, Int
k)
                    Int
col : [Int]
_
                        | Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
budget -> MathError -> Either MathError ([[Double]], [Int], Int)
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"simplex pivots" Int
k)
                        | Bool
otherwise -> case [(Double, Int, Int)] -> [(Double, Int, Int)]
forall a. Ord a => [a] -> [a]
sort
                            [ (Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
0 ([Double]
row [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
n) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ ([Double]
row [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
col), Int
bj, Int
i)
                            | (Int
i, ([Double]
row, Int
bj)) <- [Int] -> [([Double], Int)] -> [(Int, ([Double], Int))]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 :: Int ..] ([[Double]] -> [Int] -> [([Double], Int)]
forall a b. [a] -> [b] -> [(a, b)]
zip [[Double]]
rows [Int]
basis)
                            , [Double]
row [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
col Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
tol
                            ] of
                            [] -> MathError -> Either MathError ([[Double]], [Int], Int)
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"unbounded linear program")
                            (Double
_, Int
_, Int
leave) : [(Double, Int, Int)]
_ ->
                                let normalized :: [Double]
normalized = Double -> [Double] -> [Double]
scale (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (([[Double]]
rows [[Double]] -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
leave) [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
col)) ([[Double]]
rows [[Double]] -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
leave)
                                    next :: [[Double]]
next =
                                        [ if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
leave then [Double]
normalized else [Double] -> [Double] -> [Double]
sub [Double]
row (Double -> [Double] -> [Double]
scale ([Double]
row [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
col) [Double]
normalized)
                                        | (Int
i, [Double]
row) <- [Int] -> [[Double]] -> [(Int, [Double])]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 ..] [[Double]]
rows
                                        ]
                                    bs :: [Int]
bs = [if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
leave then Int
col else Int
bj | (Int
i, Int
bj) <- [Int] -> [Int] -> [(Int, Int)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 ..] [Int]
basis]
                                 in [[Double]]
-> [Int] -> Int -> Either MathError ([[Double]], [Int], Int)
go [[Double]]
next [Int]
bs (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)

{- | Revised primal simplex, A*x=b,x>=0. A basis must have m distinct columns
and B^-1*b >=0. B is solved afresh each pivot; no tableau is materialized.
-}
revisedSimplex ::
    Int -> Double -> [Double] -> [[Double]] -> [Double] -> [Int] -> Either MathError LPResult
revisedSimplex :: Int
-> Double
-> [Double]
-> [[Double]]
-> [Double]
-> [Int]
-> Either MathError LPResult
revisedSimplex Int
budget Double
tol [Double]
c [[Double]]
a [Double]
b [Int]
initial = do
    [Double] -> Either MathError ()
checkedPoint (Double
tol Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: [Double]
c [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [[Double]] -> [Double]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [[Double]]
a [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double]
b)
    let n :: Int
n = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
c; m :: Int
m = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
b
    if Int
budget Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
1
        Bool -> Bool -> Bool
|| Double
tol Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
0
        Bool -> Bool -> Bool
|| [[Double]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Double]]
a Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
m
        Bool -> Bool -> Bool
|| ([Double] -> Bool) -> [[Double]] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any ((Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
n) (Int -> Bool) -> ([Double] -> Int) -> [Double] -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length) [[Double]]
a
        Bool -> Bool -> Bool
|| [Int] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Int]
initial Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
m
        Bool -> Bool -> Bool
|| [Int] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ([Int] -> [Int]
forall a. Eq a => [a] -> [a]
nub [Int]
initial) Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
m
        Bool -> Bool -> Bool
|| (Int -> Bool) -> [Int] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (\Int
j -> Int
j Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 Bool -> Bool -> Bool
|| Int
j Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
n) [Int]
initial
        then MathError -> Either MathError LPResult
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"revised simplex shape/basis/budget")
        else Int -> [Int] -> Either MathError LPResult
go Int
0 [Int]
initial
  where
    go :: Int -> [Int] -> Either MathError LPResult
go Int
k [Int]
basis
        | Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
budget = MathError -> Either MathError LPResult
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"revised simplex" Int
k)
        | Bool
otherwise = do
            let bm :: [[Double]]
bm = [[[Double]
row [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j | Int
j <- [Int]
basis] | [Double]
row <- [[Double]]
a]
            xb <- [[Double]] -> Gradient
solveDense [[Double]]
bm [Double]
b
            if any (< -tol) xb then Left (DomainError "revised simplex basis is infeasible") else Right ()
            prices <- solveDense (transpose bm) [c !! j | j <- basis]
            let reduced = [(Int
j, [Double]
c [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double] -> [Double] -> Double
dot [Double]
prices [[Double]
row [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j | [Double]
row <- [[Double]]
a]) | Int
j <- [Int
0 .. [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
c Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1], Int
j Int -> [Int] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`notElem` [Int]
basis]
            case [j | (j, r) <- reduced, r < -tol] of
                [] ->
                    let x :: [Double]
x = [[Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
0 Double
v | (Int
j, Double
v) <- [Int] -> [Double] -> [(Int, Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int]
basis [Double]
xb, Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
i] | Int
i <- [Int
0 .. [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
c Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
                     in LPResult -> Either MathError LPResult
forall a b. b -> Either a b
Right ([Double] -> Double -> [Int] -> Int -> LPResult
LPResult [Double]
x ([Double] -> [Double] -> Double
dot [Double]
c [Double]
x) [Int]
basis Int
k)
                Int
entering : [Int]
_
                    | Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
budget -> MathError -> Either MathError LPResult
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"revised simplex" Int
k)
                    | Bool
otherwise -> do
                        direction <- [[Double]] -> Gradient
solveDense [[Double]]
bm [[Double]
row [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
entering | [Double]
row <- [[Double]]
a]
                        case sort
                            [(max 0 x / d, bj, i) | (i, (x, d, bj)) <- zip [0 :: Int ..] (zip3 xb direction basis), d > tol] of
                            [] -> MathError -> Either MathError LPResult
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"unbounded revised linear program")
                            (Double
_, Int
_, Int
leaving) : [(Double, Int, Int)]
_ -> Int -> [Int] -> Either MathError LPResult
go (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
leaving then Int
entering else Int
j | (Int
i, Int
j) <- [Int] -> [Int] -> [(Int, Int)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 ..] [Int]
basis]

{- | Strictly convex QP: min .5*x'Hx+c'x, A*x<=b. Active-set enumeration
provides a correctness reference and a hard subset budget. H must be SPD.
-}
quadraticProgramming ::
    Int -> Double -> [[Double]] -> [Double] -> [[Double]] -> [Double] -> Either MathError OptimizeResult
quadraticProgramming :: Int
-> Double
-> [[Double]]
-> [Double]
-> [[Double]]
-> [Double]
-> Either MathError OptimizeResult
quadraticProgramming Int
budget Double
tol [[Double]]
h [Double]
c [[Double]]
a [Double]
b = do
    let n :: Int
n = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
c; m :: Int
m = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
b
    [Double] -> Either MathError ()
checkedPoint (Double
tol Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: [Double]
c [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [[Double]] -> [Double]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [[Double]]
h [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [[Double]] -> [Double]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [[Double]]
a [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double]
b)
    if Int
budget Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
1
        Bool -> Bool -> Bool
|| Double
tol Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
0
        Bool -> Bool -> Bool
|| [[Double]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Double]]
h Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
n
        Bool -> Bool -> Bool
|| ([Double] -> Bool) -> [[Double]] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any ((Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
n) (Int -> Bool) -> ([Double] -> Int) -> [Double] -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length) [[Double]]
h
        Bool -> Bool -> Bool
|| [[Double]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Double]]
a Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
m
        Bool -> Bool -> Bool
|| ([Double] -> Bool) -> [[Double]] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any ((Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
n) (Int -> Bool) -> ([Double] -> Int) -> [Double] -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length) [[Double]]
a
        then MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"QP shapes/budget")
        else () -> Either MathError ()
forall a b. b -> Either a b
Right ()
    if ((Int, Int) -> Bool) -> [(Int, Int)] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any
        (\(Int
i, Int
j) -> Double -> Double
forall a. Num a => a -> a
abs (([[Double]]
h [[Double]] -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j Double -> Double -> Double
forall a. Num a => a -> a -> a
- ([[Double]]
h [[Double]] -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
j) [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
tol)
        [(Int
i, Int
j) | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1], Int
j <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
        then MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"QP Hessian must be symmetric")
        else () -> Either MathError ()
forall a b. b -> Either a b
Right ()
    _ <- [[Double]] -> Either MathError [ZonkAny 0]
forall {a} {a}.
(Ord a, Fractional a) =>
[[a]] -> Either MathError [a]
spd [[Double]]
h
    let sets = Int -> [[Int]] -> [[Int]]
forall a. Int -> [a] -> [a]
take (Int
budget Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) (([Int] -> Bool) -> [[Int]] -> [[Int]]
forall a. (a -> Bool) -> [a] -> [a]
filter ((Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
n) (Int -> Bool) -> ([Int] -> Int) -> [Int] -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Int] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length) ([Int] -> [[Int]]
forall {a}. [a] -> [[a]]
subsets [Int
0 .. Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]))
    if length sets > budget then Left (NonConvergence "QP active-set budget" budget) else Right ()
    candidates <- mapM solveActive sets
    case sort [(y, x, r) | (Just (y, x, r)) <- candidates] of
        [] -> MathError -> Either MathError OptimizeResult
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"infeasible convex quadratic program")
        (Double
y, [Double]
x, Double
r) : [(Double, [Double], Double)]
_ -> OptimizeResult -> Either MathError OptimizeResult
forall a b. b -> Either a b
Right ([Double]
-> Double -> Double -> Int -> Int -> Bool -> OptimizeResult
OptimizeResult [Double]
x Double
y Double
r ([[Int]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Int]]
sets) ([[Int]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Int]]
sets) Bool
True)
  where
    subsets :: [a] -> [[a]]
subsets [] = [[]]
    subsets (a
x : [a]
xs) = let rest :: [[a]]
rest = [a] -> [[a]]
subsets [a]
xs in [[a]]
rest [[a]] -> [[a]] -> [[a]]
forall a. [a] -> [a] -> [a]
++ ([a] -> [a]) -> [[a]] -> [[a]]
forall a b. (a -> b) -> [a] -> [b]
map (a
x a -> [a] -> [a]
forall a. a -> [a] -> [a]
:) [[a]]
rest
    spd :: [[a]] -> Either MathError [a]
spd [] = [a] -> Either MathError [a]
forall a b. b -> Either a b
Right []
    spd [[a]]
rows = case [[a]]
rows of
        (a
pivot : [a]
rs) : [[a]]
others | a
pivot a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> a
0 -> [[a]] -> Either MathError [a]
spd [[a
v a -> a -> a
forall a. Num a => a -> a -> a
- a
x a -> a -> a
forall a. Num a => a -> a -> a
* a
y a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
pivot | (a
v, a
y) <- [a] -> [a] -> [(a, a)]
forall a b. [a] -> [b] -> [(a, b)]
zip [a]
row [a]
rs] | (a
x : [a]
row) <- [[a]]
others]
        [[a]]
_ -> MathError -> Either MathError [a]
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"QP requires positive definite Hessian")
    solveActive :: [Int] -> Either MathError (Maybe (Double, [Double], Double))
solveActive [Int]
active = do
        let n :: Int
n = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
c
            ac :: [[Double]]
ac = [[[Double]]
a [[Double]] -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i | Int
i <- [Int]
active]
            bc :: [Double]
bc = [[Double]
b [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i | Int
i <- [Int]
active]
            q :: Int
q = [Int] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Int]
active
            kkt :: [[Double]]
kkt = [[Double]
row [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [[Double]
ar [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i | [Double]
ar <- [[Double]]
ac] | (Int
i, [Double]
row) <- [Int] -> [[Double]] -> [(Int, [Double])]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 ..] [[Double]]
h] [[Double]] -> [[Double]] -> [[Double]]
forall a. [a] -> [a] -> [a]
++ [[Double]
row [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate Int
q Double
0 | [Double]
row <- [[Double]]
ac]
            rhs :: [Double]
rhs = (Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Double -> Double
forall a. Num a => a -> a
negate [Double]
c [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double]
bc
        case [[Double]] -> Gradient
solveDense [[Double]]
kkt [Double]
rhs of
            Left (SingularMatrix Int
_) -> Maybe (Double, [Double], Double)
-> Either MathError (Maybe (Double, [Double], Double))
forall a b. b -> Either a b
Right Maybe (Double, [Double], Double)
forall a. Maybe a
Nothing
            Left (NearSingularMatrix Int
_) -> Maybe (Double, [Double], Double)
-> Either MathError (Maybe (Double, [Double], Double))
forall a b. b -> Either a b
Right Maybe (Double, [Double], Double)
forall a. Maybe a
Nothing
            Left MathError
err -> MathError -> Either MathError (Maybe (Double, [Double], Double))
forall a b. a -> Either a b
Left MathError
err
            Right [Double]
solution -> do
                let x :: [Double]
x = Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
take Int
n [Double]
solution
                    multipliers :: [Double]
multipliers = Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
n [Double]
solution
                    violation :: Double
violation = [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum (Double
0 Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: ([Double] -> Double -> Double)
-> [[Double]] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\[Double]
row Double
bound -> [Double] -> [Double] -> Double
dot [Double]
row [Double]
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
bound) [[Double]]
a [Double]
b)
                    dual :: [Double]
dual =
                        [Double] -> [Double] -> [Double]
sub
                            ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Double -> Double
forall a. Num a => a -> a
negate [Double]
c)
                            ([Double] -> [Double] -> [Double]
add ([[Double]] -> [Double] -> [Double]
matVec [[Double]]
h [Double]
x) ([[Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double
l Double -> Double -> Double
forall a. Num a => a -> a -> a
* ([Double]
row [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) | (Double
l, [Double]
row) <- [Double] -> [[Double]] -> [(Double, [Double])]
forall a b. [a] -> [b] -> [(a, b)]
zip [Double]
multipliers [[Double]]
ac] | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]))
                if Double
violation Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
tol Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= -Double
tol) [Double]
multipliers
                    then Maybe (Double, [Double], Double)
-> Either MathError (Maybe (Double, [Double], Double))
forall a b. b -> Either a b
Right ((Double, [Double], Double) -> Maybe (Double, [Double], Double)
forall a. a -> Maybe a
Just (Double
0.5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double] -> [Double] -> Double
dot [Double]
x ([[Double]] -> [Double] -> [Double]
matVec [[Double]]
h [Double]
x) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ [Double] -> [Double] -> Double
dot [Double]
c [Double]
x, [Double]
x, Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
violation ([Double] -> Double
norm [Double]
dual)))
                    else Maybe (Double, [Double], Double)
-> Either MathError (Maybe (Double, [Double], Double))
forall a b. b -> Either a b
Right Maybe (Double, [Double], Double)
forall a. Maybe a
Nothing

{- | Sequential least-squares QP with an analytic Hessian and exact linearized
equality/active inequality constraints. Merit backtracking globalizes steps.
-}
slsqp ::
    OptimizeOptions ->
    Objective ->
    Gradient ->
    Hessian ->
    [Constraint] ->
    [Double] ->
    Either MathError OptimizeResult
slsqp :: OptimizeOptions
-> Objective
-> Gradient
-> Hessian
-> [Constraint]
-> [Double]
-> Either MathError OptimizeResult
slsqp OptimizeOptions
o Objective
f Gradient
g Hessian
hess [Constraint]
constraints [Double]
x0 = do OptimizeOptions -> Either MathError ()
checkedOptions OptimizeOptions
o; Int -> Int -> [Double] -> Either MathError OptimizeResult
loop Int
0 Int
0 [Double]
x0
  where
    baseCost :: Int
baseCost = Int
3 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* [Constraint] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Constraint]
constraints
    meritCost :: Int
meritCost = Int
1 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* [Constraint] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Constraint]
constraints
    loop :: Int -> Int -> [Double] -> Either MathError OptimizeResult
loop Int
k Int
cost [Double]
x = do
        OptimizeOptions -> Int -> Either MathError ()
checkBudget OptimizeOptions
o (Int
cost Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
baseCost)
        y <- Objective -> Objective
checkedValue Objective
f [Double]
x
        gs <- checkedGradient g x
        h <- hess x
        vals <- constraintValues constraints x
        let initial = [(Bool
eq, Double
v, [Double]
d) | (Bool
eq, Double
v, [Double]
d) <- [(Bool, Double, [Double])]
vals, Bool
eq Bool -> Bool -> Bool
|| Double
v Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= -Double -> Double
forall a. Floating a => a -> a
sqrt (OptimizeOptions -> Double
optTolerance OptimizeOptions
o)]
            n = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
x
        if length h /= n || any ((/= n) . length) h
            then Left (InvalidArgument "SQP Hessian shape")
            else checkedPoint (concat h)
        (active, solution) <- activeSolve h gs initial
        let a = [[Double]
d | (Bool
_, Double
_, [Double]
d) <- [(Bool, Double, [Double])]
active]
            p = Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
take Int
n [Double]
solution
            lambda = Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
n [Double]
solution
            dual = [Double] -> Double
norm ([Double] -> [Double] -> [Double]
add [Double]
gs [[Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double
l Double -> Double -> Double
forall a. Num a => a -> a -> a
* ([Double]
row [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) | (Double
l, [Double]
row) <- [Double] -> [[Double]] -> [(Double, [Double])]
forall a b. [a] -> [b] -> [(a, b)]
zip [Double]
lambda [[Double]]
a] | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]])
            violation = [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum (Double
0 Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: [if Bool
eq then Double -> Double
forall a. Num a => a -> a
abs Double
v else Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
0 Double
v | (Bool
eq, Double
v, [Double]
_) <- [(Bool, Double, [Double])]
vals])
            used = Int
cost Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
baseCost
        if dual <= optTolerance o && violation <= optTolerance o
            then Right (OptimizeResult x y (max dual violation) k used True)
            else
                if k >= optMaxIterations o
                    then Left (NonConvergence "SLSQP" k)
                    else do
                        let rho = Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum (Double
0 Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: (Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Double -> Double
forall a. Num a => a -> a
abs [Double]
lambda)
                            merit [Double]
z = do
                                yz <- Objective -> Objective
checkedValue Objective
f [Double]
z
                                vs <- constraintValues constraints z
                                Right (yz + rho * sum [if eq then abs v else max 0 v | (eq, v, _) <- vs])
                        checkBudget o (used + meritCost)
                        current <- merit x
                        (next, count) <- backtrack (optMaxEvaluations o - used - meritCost) merit current p x 0 1
                        loop (k + 1) (used + meritCost + count) next
    activeSolve :: [[Double]]
-> [Double]
-> [(Bool, Double, [Double])]
-> Either MathError ([(Bool, Double, [Double])], [Double])
activeSolve [[Double]]
h [Double]
gs [(Bool, Double, [Double])]
active = do
        let n :: Int
n = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
gs
            q :: Int
q = [(Bool, Double, [Double])] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [(Bool, Double, [Double])]
active
            a :: [[Double]]
a = [[Double]
d | (Bool
_, Double
_, [Double]
d) <- [(Bool, Double, [Double])]
active]
            kkt :: [[Double]]
kkt = [[Double]
row [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [[Double]
ar [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i | [Double]
ar <- [[Double]]
a] | (Int
i, [Double]
row) <- [Int] -> [[Double]] -> [(Int, [Double])]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 ..] [[Double]]
h] [[Double]] -> [[Double]] -> [[Double]]
forall a. [a] -> [a] -> [a]
++ [[Double]
row [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate Int
q Double
0 | [Double]
row <- [[Double]]
a]
        solution <- [[Double]] -> Gradient
solveDense [[Double]]
kkt ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Double -> Double
forall a. Num a => a -> a
negate [Double]
gs [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [-Double
v | (Bool
_, Double
v, [Double]
_) <- [(Bool, Double, [Double])]
active])
        case [ i
             | (i, ((eq, _, _), lambda)) <- zip [0 :: Int ..] (zip active (drop n solution))
             , not eq
             , lambda < -optTolerance o
             ] of
            [] -> ([(Bool, Double, [Double])], [Double])
-> Either MathError ([(Bool, Double, [Double])], [Double])
forall a b. b -> Either a b
Right ([(Bool, Double, [Double])]
active, [Double]
solution)
            Int
index : [Int]
_ -> [[Double]]
-> [Double]
-> [(Bool, Double, [Double])]
-> Either MathError ([(Bool, Double, [Double])], [Double])
activeSolve [[Double]]
h [Double]
gs [(Bool, Double, [Double])
a' | (Int
i, (Bool, Double, [Double])
a') <- [Int]
-> [(Bool, Double, [Double])] -> [(Int, (Bool, Double, [Double]))]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 ..] [(Bool, Double, [Double])]
active, Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
index]
    backtrack :: Int
-> ([Double] -> Either MathError t)
-> t
-> [Double]
-> [Double]
-> Int
-> Double
-> Either MathError ([Double], Int)
backtrack Int
remaining [Double] -> Either MathError t
merit t
current [Double]
p [Double]
x Int
count Double
rate
        | Int
count Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
meritCost Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
remaining = MathError -> Either MathError ([Double], Int)
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"SQP merit evaluation budget" Int
count)
        | Int
count Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
50 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
meritCost = MathError -> Either MathError ([Double], Int)
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"SQP merit line search" Int
count)
        | Bool
otherwise = do
            let next :: [Double]
next = [Double] -> [Double] -> [Double]
add [Double]
x (Double -> [Double] -> [Double]
scale Double
rate [Double]
p)
            value <- [Double] -> Either MathError t
merit [Double]
next
            if value < current
                then Right (next, count + meritCost)
                else backtrack remaining merit current p x (count + meritCost) (rate / 2)

{- | Derivative-free constrained linear-approximation trust region. This bounded
CPU variant uses an L-infinity trust box and axis interpolation simplex.
It solves each linear model as an LP, rather than relabeling a penalty solver.
-}
cobyla ::
    Double -> OptimizeOptions -> Objective -> [Objective] -> [Double] -> Either MathError OptimizeResult
cobyla :: Double
-> OptimizeOptions
-> Objective
-> [Objective]
-> [Double]
-> Either MathError OptimizeResult
cobyla Double
initialRadius OptimizeOptions
o Objective
f [Objective]
constraints [Double]
x0 = do
    OptimizeOptions -> Either MathError ()
checkedOptions OptimizeOptions
o
    [Double] -> Either MathError ()
checkedPoint [Double
initialRadius]
    if Double
initialRadius Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
0 Bool -> Bool -> Bool
|| [Double] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Double]
x0
        then MathError -> Either MathError OptimizeResult
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"COBYLA trust radius/dimension")
        else Int -> Int -> Double -> [Double] -> Either MathError OptimizeResult
loop Int
0 Int
0 Double
initialRadius [Double]
x0
  where
    loop :: Int -> Int -> Double -> [Double] -> Either MathError OptimizeResult
loop Int
k Int
cost Double
radius [Double]
x = do
        let n :: Int
n = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
x; funcs :: [Objective]
funcs = Objective
f Objective -> [Objective] -> [Objective]
forall a. a -> [a] -> [a]
: [Objective]
constraints; calls :: Int
calls = (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2) Int -> Int -> Int
forall a. Num a => a -> a -> a
* [Objective] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Objective]
funcs
        OptimizeOptions -> Int -> Either MathError ()
checkBudget OptimizeOptions
o (Int
cost Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
calls Int -> Int -> Int
forall a. Num a => a -> a -> a
+ [Objective] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Objective]
funcs)
        values <- (Objective -> Either MathError Double)
-> [Objective] -> Either MathError [Double]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
forall (m :: * -> *) a b. Monad m => (a -> m b) -> [a] -> m [b]
mapM (Objective -> Objective
`checkedValue` [Double]
x) [Objective]
funcs
        model <-
            mapM
                ( \Objective
fn -> do
                    base <- Objective -> Objective
checkedValue Objective
fn [Double]
x
                    mapM
                        ( \Int
i -> do
                            y <- Objective -> Objective
checkedValue Objective
fn ([Double] -> [Double] -> [Double]
add [Double]
x [if Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
i then Double
radius else Double
0 | Int
j <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]])
                            Right ((y - base) / radius)
                        )
                        [0 .. n - 1]
                )
                funcs
        let c = [[Double]]
model [[Double]] -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
0
            rows = Int -> [[Double]] -> [[Double]]
forall a. Int -> [a] -> [a]
drop Int
1 [[Double]]
model
            rhs = (Double -> [Double] -> Double)
-> [Double] -> [[Double]] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Double
v [Double]
row -> -Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
radius Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double]
row) (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
1 [Double]
values) [[Double]]
rows
            bounds = Int -> [[Double]]
identity Int
n
            aa = [[Double]]
rows [[Double]] -> [[Double]] -> [[Double]]
forall a. [a] -> [a] -> [a]
++ [[Double]]
bounds
            bb = [Double]
rhs [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate Int
n (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
radius)
        step <- case linearProgramming (max 100 (optMaxIterations o)) (optTolerance o / 100) c aa bb of
            Right LPResult
lp -> Gradient
forall a b. b -> Either a b
Right ([Double] -> [Double] -> [Double]
sub (LPResult -> [Double]
lpPoint LPResult
lp) (Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate Int
n Double
radius))
            Left (DomainError String
_) -> Gradient
forall a b. b -> Either a b
Right (Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate Int
n Double
0)
            Left MathError
e -> MathError -> Either MathError [Double]
forall a b. a -> Either a b
Left MathError
e
        nextVals <- mapM (\Objective
fn -> Objective -> Objective
checkedValue Objective
fn ([Double] -> [Double] -> [Double]
add [Double]
x [Double]
step)) funcs
        let merit [a]
vs = [a]
vs [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
0 a -> a -> a
forall a. Num a => a -> a -> a
+ a
1000 a -> a -> a
forall a. Num a => a -> a -> a
* [a] -> a
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a -> a -> a
forall a. Ord a => a -> a -> a
max a
0) (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
drop Int
1 [a]
vs))
            improve = [Double] -> Double
forall {a}. (Num a, Ord a) => [a] -> a
merit [Double]
nextVals Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< [Double] -> Double
forall {a}. (Num a, Ord a) => [a] -> a
merit [Double]
values
            next = if Bool
improve then [Double] -> [Double] -> [Double]
add [Double]
x [Double]
step else [Double]
x
            violation = [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum (Double
0 Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: (Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
0) (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
1 (if Bool
improve then [Double]
nextVals else [Double]
values)))
        if radius <= optTolerance o && violation <= optTolerance o
            then
                Right
                    ( OptimizeResult
                        next
                        (if improve then nextVals !! 0 else values !! 0)
                        (max radius violation)
                        k
                        (cost + calls + length funcs)
                        True
                    )
            else
                if k >= optMaxIterations o
                    then Left (NonConvergence "COBYLA" k)
                    else
                        loop
                            (k + 1)
                            (cost + calls + length funcs)
                            (if improve && norm step > radius / 2 then radius else radius / 2)
                            next

-- | Compose nonnegative scalarization weights and checked objective values; at least one weight must be positive.
weightedObjective :: [Double] -> [Objective] -> Objective
weightedObjective :: [Double] -> [Objective] -> Objective
weightedObjective [Double]
weights [Objective]
objectives [Double]
x = do
    [Double] -> Either MathError ()
checkedPoint [Double]
weights
    if [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
weights Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= [Objective] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Objective]
objectives Bool -> Bool -> Bool
|| [Double] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Double]
weights Bool -> Bool -> Bool
|| (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0) [Double]
weights Bool -> Bool -> Bool
|| [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double]
weights Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
0
        then MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"multiobjective weights")
        else () -> Either MathError ()
forall a b. b -> Either a b
Right ()
    vals <- (Objective -> Either MathError Double)
-> [Objective] -> Either MathError [Double]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
forall (m :: * -> *) a b. Monad m => (a -> m b) -> [a] -> m [b]
mapM (Objective -> Objective
`checkedValue` [Double]
x) [Objective]
objectives
    Right (dot weights vals)

-- | Return original indices of nondominated objective rows under componentwise minimization; duplicates remain nondominated.
paretoFront :: [[Double]] -> Either MathError [Int]
paretoFront :: [[Double]] -> Either MathError [Int]
paretoFront [[Double]]
points = do
    [Double] -> Either MathError ()
checkedPoint ([[Double]] -> [Double]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [[Double]]
points)
    case [[Double]]
points of
        [] -> [Int] -> Either MathError [Int]
forall a b. b -> Either a b
Right []
        [Double]
first : [[Double]]
rest ->
            if ([Double] -> Bool) -> [[Double]] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any ((Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
first) (Int -> Bool) -> ([Double] -> Int) -> [Double] -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length) [[Double]]
rest
                then MathError -> Either MathError [Int]
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"Pareto objective shape")
                else [Int] -> Either MathError [Int]
forall a b. b -> Either a b
Right [Int
i | (Int
i, [Double]
x) <- [Int] -> [[Double]] -> [(Int, [Double])]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 ..] [[Double]]
points, Bool -> Bool
not (([Double] -> Bool) -> [[Double]] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any ([Double] -> [Double] -> Bool
forall {b}. Ord b => [b] -> [b] -> Bool
dominates [Double]
x) [[Double]]
points)]
  where
    dominates :: [b] -> [b] -> Bool
dominates [b]
x [b]
y = [Bool] -> Bool
forall (t :: * -> *). Foldable t => t Bool -> Bool
and ((b -> b -> Bool) -> [b] -> [b] -> [Bool]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith b -> b -> Bool
forall a. Ord a => a -> a -> Bool
(<=) [b]
y [b]
x) Bool -> Bool -> Bool
&& [Bool] -> Bool
forall (t :: * -> *). Foldable t => t Bool -> Bool
or ((b -> b -> Bool) -> [b] -> [b] -> [Bool]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith b -> b -> Bool
forall a. Ord a => a -> a -> Bool
(<) [b]
y [b]
x)