{-# OPTIONS_GHC -Wall #-}
{-# OPTIONS_HADDOCK show-extensions #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE TypeFamilies #-}
-----------------------------------------------------------------------------
-- |
-- Module      :  ToySolver.SDP
-- Copyright   :  (c) Masahiro Sakai 2017
-- License     :  BSD-style
--
-- Maintainer  :  masahiro.sakai@gmail.com
-- Stability   :  exprimental
-- Portability :  non-portable
--
-- References:
--
-- * Convert Semidefinite program forms - Mathematics Stack Exchange
--   <https://math.stackexchange.com/questions/732658/convert-semidefinite-program-forms>
--
-----------------------------------------------------------------------------

module ToySolver.SDP
  ( dualize
  , DualizeInfo (..)
  ) where

import qualified Data.Aeson as J
import Data.Aeson ((.=), (.:))
import qualified Data.Map.Strict as Map
import Data.Scientific (Scientific)
import ToySolver.Converter.Base
import ToySolver.Internal.JSON (withTypedObject)
import qualified ToySolver.Text.SDPFile as SDPFile

-- | Given a primal-dual pair (P), (D), it returns another primal-dual pair (P'), (D')
-- such that (P) is equivalent to (D') and (D) is equivalent to (P').
dualize :: SDPFile.Problem -> (SDPFile.Problem, DualizeInfo)
dualize :: Problem -> (Problem, DualizeInfo)
dualize Problem
origProb =
  ( SDPFile.Problem
    { blockStruct :: [Int]
SDPFile.blockStruct = [Int]
blockStruct
    , costs :: [Scientific]
SDPFile.costs = [Scientific]
d
    , matrices :: [Matrix]
SDPFile.matrices = Matrix
a0 Matrix -> [Matrix] -> [Matrix]
forall a. a -> [a] -> [a]
: [Matrix]
as
    }
  , Int -> [Int] -> DualizeInfo
DualizeInfo Int
m (Problem -> [Int]
SDPFile.blockStruct Problem
origProb)
  )
  where
    {- original:
       (P)
           min Σ_i=1^m c_i x_i
           s.t.
             X = Σ_i=1^m F_i x_i - F_0
             X ⪰ 0
       (D)
           max F_0 • Y
           s.t.
             F_i • Y = c_i  for i ∈ {1,…,m}
             Y ⪰ 0
       where
         x : variable over R^m
         c ∈ R^m
         F_0, F_1, … , F_m ∈ R^(n × n)
    -}
    m :: Int
    m :: Int
m = Problem -> Int
SDPFile.mDim Problem
origProb
    c :: [Scientific]
    c :: [Scientific]
c = Problem -> [Scientific]
SDPFile.costs Problem
origProb
    f0 :: SDPFile.Matrix
    fs :: [SDPFile.Matrix]
    Matrix
f0:[Matrix]
fs = Problem -> [Matrix]
SDPFile.matrices Problem
origProb

    {- transformed
       (P')
           min d^T・z
           s.t.
             Z = Σ_i=1^n Σ_j=1^i A_ij z_ij - A_0
             Z ⪰ 0
       (D')
           max A_0 • W
           s.t.
             A_ij • W = d_ij  for i ∈ {1,…,n}, j ∈ {1,…,i}
             W ⪰ 0
       where
         z : variable over R^{n(n+1)/2}
         d_ij ∈ R  for i ∈ {1,…,n}, j ∈ {1,…,i}
         d_ij = - F0 [i,j]              if i=j
              = - (F0 [i,j] + F0 [j,i]) otherwise
         A_0  ∈ R^((2m+n)×(2m+n))
         A_0  = diag(-c, c, 0_{n×n})
         A_ij ∈ R^((2m+n)×(2m+n))  for i ∈ {1,…,n}, j ∈ {1,…,i}
         A_ij [   k,   k] = - (if i=j then F_k [i,j] else F_k [i,j] + F_k [j,i]) for k∈{1,…,m}
         A_ij [ m+k, m+k] =   (if i=j then F_k [i,j] else F_k [i,j] + F_k [j,i]) for k∈{1,…,m}
         A_ij [2m+i,2m+j] = 1
         A_ij [2m+j,2m+i] = 1
         A_ij [  _ ,  _ ] = 0

       correspondence:
         W = diag(x+, x-, X)
         Y [i,j] = z_ij if j≤i
                 = z_ji otherwise
         Z = diag(0, 0, Y)
    -}
    blockStruct :: [Int]
    blockStruct :: [Int]
blockStruct = [-Int
m, -Int
m] [Int] -> [Int] -> [Int]
forall a. [a] -> [a] -> [a]
++ Problem -> [Int]
SDPFile.blockStruct Problem
origProb
    a0 :: SDPFile.Matrix
    a0 :: Matrix
a0 =
      [ [((Int, Int), Scientific)] -> Map (Int, Int) Scientific
forall k a. Ord k => [(k, a)] -> Map k a
Map.fromList [((Int
j,Int
j), -Scientific
cj) | (Int
j,Scientific
cj) <- [Int] -> [Scientific] -> [(Int, Scientific)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
1..Int
m] [Scientific]
c, Scientific
cj Scientific -> Scientific -> Bool
forall a. Eq a => a -> a -> Bool
/= Scientific
0]
      , [((Int, Int), Scientific)] -> Map (Int, Int) Scientific
forall k a. Ord k => [(k, a)] -> Map k a
Map.fromList [((Int
j,Int
j),  Scientific
cj) | (Int
j,Scientific
cj) <- [Int] -> [Scientific] -> [(Int, Scientific)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
1..Int
m] [Scientific]
c, Scientific
cj Scientific -> Scientific -> Bool
forall a. Eq a => a -> a -> Bool
/= Scientific
0]
      ] Matrix -> Matrix -> Matrix
forall a. [a] -> [a] -> [a]
++
      [ Map (Int, Int) Scientific
forall k a. Map k a
Map.empty | Int
_ <- Problem -> [Int]
SDPFile.blockStruct Problem
origProb]
    as :: [SDPFile.Matrix]
    as :: [Matrix]
as =
      [ [ [((Int, Int), Scientific)] -> Map (Int, Int) Scientific
forall k a. Ord k => [(k, a)] -> Map k a
Map.fromList [ ((Int
k,Int
k), - (if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j then Scientific
v else Scientific
2Scientific -> Scientific -> Scientific
forall a. Num a => a -> a -> a
*Scientific
v))
                       | (Int
k,Matrix
fk) <- [Int] -> [Matrix] -> [(Int, Matrix)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
1..Int
m] [Matrix]
fs, let v :: Scientific
v = Int -> Int -> Map (Int, Int) Scientific -> Scientific
SDPFile.blockElem Int
i Int
j (Matrix
fkMatrix -> Int -> Map (Int, Int) Scientific
forall a. HasCallStack => [a] -> Int -> a
!!Int
b), Scientific
v Scientific -> Scientific -> Bool
forall a. Eq a => a -> a -> Bool
/= Scientific
0]
        , [((Int, Int), Scientific)] -> Map (Int, Int) Scientific
forall k a. Ord k => [(k, a)] -> Map k a
Map.fromList [ ((Int
k,Int
k),   (if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j then Scientific
v else Scientific
2Scientific -> Scientific -> Scientific
forall a. Num a => a -> a -> a
*Scientific
v))
                       | (Int
k,Matrix
fk) <- [Int] -> [Matrix] -> [(Int, Matrix)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
1..Int
m] [Matrix]
fs, let v :: Scientific
v = Int -> Int -> Map (Int, Int) Scientific -> Scientific
SDPFile.blockElem Int
i Int
j (Matrix
fkMatrix -> Int -> Map (Int, Int) Scientific
forall a. HasCallStack => [a] -> Int -> a
!!Int
b), Scientific
v Scientific -> Scientific -> Bool
forall a. Eq a => a -> a -> Bool
/= Scientific
0]
        ] Matrix -> Matrix -> Matrix
forall a. [a] -> [a] -> [a]
++
        [ if Int
b Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
b2 then
            Map (Int, Int) Scientific
forall k a. Map k a
Map.empty
          else if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j then
            (Int, Int) -> Scientific -> Map (Int, Int) Scientific
forall k a. k -> a -> Map k a
Map.singleton (Int
i,Int
j) Scientific
1
          else
            [((Int, Int), Scientific)] -> Map (Int, Int) Scientific
forall k a. Ord k => [(k, a)] -> Map k a
Map.fromList [((Int
i,Int
j),Scientific
1), ((Int
j,Int
i),Scientific
1)]
        | (Int
b2, Int
_) <- [Int] -> [Int] -> [(Int, Int)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0..] (Problem -> [Int]
SDPFile.blockStruct Problem
origProb)
        ]
      | (Int
b,Int
block) <- [Int] -> [Int] -> [(Int, Int)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0..] (Problem -> [Int]
SDPFile.blockStruct Problem
origProb)
      , (Int
i,Int
j) <- Int -> [(Int, Int)]
blockIndexes Int
block
      ]
    d :: [Scientific]
d =
      [ - (if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j then Scientific
v else Scientific
2Scientific -> Scientific -> Scientific
forall a. Num a => a -> a -> a
*Scientific
v)
      | (Int
b,Int
block) <- [Int] -> [Int] -> [(Int, Int)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0..] (Problem -> [Int]
SDPFile.blockStruct Problem
origProb)
      , (Int
i,Int
j) <- Int -> [(Int, Int)]
blockIndexes Int
block
      , let v :: Scientific
v = Int -> Int -> Map (Int, Int) Scientific -> Scientific
SDPFile.blockElem Int
i Int
j (Matrix
f0 Matrix -> Int -> Map (Int, Int) Scientific
forall a. HasCallStack => [a] -> Int -> a
!! Int
b)
      ]

blockIndexes :: Int -> [(Int,Int)]
blockIndexes :: Int -> [(Int, Int)]
blockIndexes Int
n = if Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 then [(Int
i,Int
j) | Int
i <- [Int
1..Int
n], Int
j <- [Int
1..Int
i]] else [(Int
i,Int
i) | Int
i <- [Int
1..(-Int
n)]]

blockIndexesLen :: Int -> Int
blockIndexesLen :: Int -> Int
blockIndexesLen Int
n = if Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 then Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
*(Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2 else -Int
n


data DualizeInfo = DualizeInfo Int [Int]
  deriving (DualizeInfo -> DualizeInfo -> Bool
(DualizeInfo -> DualizeInfo -> Bool)
-> (DualizeInfo -> DualizeInfo -> Bool) -> Eq DualizeInfo
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: DualizeInfo -> DualizeInfo -> Bool
== :: DualizeInfo -> DualizeInfo -> Bool
$c/= :: DualizeInfo -> DualizeInfo -> Bool
/= :: DualizeInfo -> DualizeInfo -> Bool
Eq, Int -> DualizeInfo -> ShowS
[DualizeInfo] -> ShowS
DualizeInfo -> String
(Int -> DualizeInfo -> ShowS)
-> (DualizeInfo -> String)
-> ([DualizeInfo] -> ShowS)
-> Show DualizeInfo
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> DualizeInfo -> ShowS
showsPrec :: Int -> DualizeInfo -> ShowS
$cshow :: DualizeInfo -> String
show :: DualizeInfo -> String
$cshowList :: [DualizeInfo] -> ShowS
showList :: [DualizeInfo] -> ShowS
Show, ReadPrec [DualizeInfo]
ReadPrec DualizeInfo
Int -> ReadS DualizeInfo
ReadS [DualizeInfo]
(Int -> ReadS DualizeInfo)
-> ReadS [DualizeInfo]
-> ReadPrec DualizeInfo
-> ReadPrec [DualizeInfo]
-> Read DualizeInfo
forall a.
(Int -> ReadS a)
-> ReadS [a] -> ReadPrec a -> ReadPrec [a] -> Read a
$creadsPrec :: Int -> ReadS DualizeInfo
readsPrec :: Int -> ReadS DualizeInfo
$creadList :: ReadS [DualizeInfo]
readList :: ReadS [DualizeInfo]
$creadPrec :: ReadPrec DualizeInfo
readPrec :: ReadPrec DualizeInfo
$creadListPrec :: ReadPrec [DualizeInfo]
readListPrec :: ReadPrec [DualizeInfo]
Read)

instance Transformer DualizeInfo where
  type Source DualizeInfo = SDPFile.Solution
  type Target DualizeInfo = SDPFile.Solution

instance ForwardTransformer DualizeInfo where
  transformForward :: DualizeInfo -> Source DualizeInfo -> Target DualizeInfo
transformForward (DualizeInfo Int
_origM [Int]
origBlockStruct)
    SDPFile.Solution
    { primalVector :: Solution -> [Scientific]
SDPFile.primalVector = [Scientific]
xV
    , primalMatrix :: Solution -> Matrix
SDPFile.primalMatrix = Matrix
xM
    , dualMatrix :: Solution -> Matrix
SDPFile.dualMatrix   = Matrix
yM
    } =
    SDPFile.Solution
    { primalVector :: [Scientific]
SDPFile.primalVector = [Scientific]
zV
    , primalMatrix :: Matrix
SDPFile.primalMatrix = Matrix
zM
    , dualMatrix :: Matrix
SDPFile.dualMatrix   = Matrix
wM
    }
    where
      zV :: [Scientific]
zV = [[Scientific]] -> [Scientific]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [[Int -> Int -> Map (Int, Int) Scientific -> Scientific
SDPFile.blockElem Int
i Int
j Map (Int, Int) Scientific
block | (Int
i,Int
j) <- Int -> [(Int, Int)]
blockIndexes Int
b] | (Int
b,Map (Int, Int) Scientific
block) <- [Int] -> Matrix -> [(Int, Map (Int, Int) Scientific)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int]
origBlockStruct Matrix
yM]
      zM :: Matrix
zM = Map (Int, Int) Scientific
forall k a. Map k a
Map.empty Map (Int, Int) Scientific -> Matrix -> Matrix
forall a. a -> [a] -> [a]
: Map (Int, Int) Scientific
forall k a. Map k a
Map.empty Map (Int, Int) Scientific -> Matrix -> Matrix
forall a. a -> [a] -> [a]
: Matrix
yM
      wM :: Matrix
wM =
        [ [((Int, Int), Scientific)] -> Map (Int, Int) Scientific
forall k a. Ord k => [(k, a)] -> Map k a
Map.fromList ([((Int, Int), Scientific)] -> Map (Int, Int) Scientific)
-> [((Int, Int), Scientific)] -> Map (Int, Int) Scientific
forall a b. (a -> b) -> a -> b
$ (Int -> Scientific -> ((Int, Int), Scientific))
-> [Int] -> [Scientific] -> [((Int, Int), Scientific)]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Int
i Scientific
x -> ((Int
i,Int
i), if Scientific
x Scientific -> Scientific -> Bool
forall a. Ord a => a -> a -> Bool
>= Scientific
0 then   Scientific
x else Scientific
0)) [Int
1..] [Scientific]
xV
        , [((Int, Int), Scientific)] -> Map (Int, Int) Scientific
forall k a. Ord k => [(k, a)] -> Map k a
Map.fromList ([((Int, Int), Scientific)] -> Map (Int, Int) Scientific)
-> [((Int, Int), Scientific)] -> Map (Int, Int) Scientific
forall a b. (a -> b) -> a -> b
$ (Int -> Scientific -> ((Int, Int), Scientific))
-> [Int] -> [Scientific] -> [((Int, Int), Scientific)]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Int
i Scientific
x -> ((Int
i,Int
i), if Scientific
x Scientific -> Scientific -> Bool
forall a. Ord a => a -> a -> Bool
<= Scientific
0 then  -Scientific
x else Scientific
0)) [Int
1..] [Scientific]
xV
        ] Matrix -> Matrix -> Matrix
forall a. [a] -> [a] -> [a]
++ Matrix
xM

instance BackwardTransformer DualizeInfo where
  transformBackward :: DualizeInfo -> Target DualizeInfo -> Source DualizeInfo
transformBackward (DualizeInfo Int
origM [Int]
origBlockStruct)
    SDPFile.Solution
    { primalVector :: Solution -> [Scientific]
SDPFile.primalVector = [Scientific]
zV
    , primalMatrix :: Solution -> Matrix
SDPFile.primalMatrix = Matrix
_zM
    , dualMatrix :: Solution -> Matrix
SDPFile.dualMatrix   = Matrix
wM
    } =
    case Matrix
wM of
      (Map (Int, Int) Scientific
xps:Map (Int, Int) Scientific
xns:Matrix
xM) ->
        SDPFile.Solution
        { primalVector :: [Scientific]
SDPFile.primalVector = [Scientific]
xV
        , primalMatrix :: Matrix
SDPFile.primalMatrix = Matrix
xM
        , dualMatrix :: Matrix
SDPFile.dualMatrix   = Matrix
yM
        }
        where
          xV :: [Scientific]
xV = [Int -> Int -> Map (Int, Int) Scientific -> Scientific
SDPFile.blockElem Int
i Int
i Map (Int, Int) Scientific
xps Scientific -> Scientific -> Scientific
forall a. Num a => a -> a -> a
- Int -> Int -> Map (Int, Int) Scientific -> Scientific
SDPFile.blockElem Int
i Int
i Map (Int, Int) Scientific
xns | Int
i <- [Int
1..Int
origM]]
          yM :: Matrix
yM = [Int] -> [Scientific] -> Matrix
f [Int]
origBlockStruct [Scientific]
zV
            where
              f :: [Int] -> [Scientific] -> Matrix
f [] [Scientific]
_ = []
              f (Int
block : [Int]
blocks) [Scientific]
zV1 =
                case Int -> [Scientific] -> ([Scientific], [Scientific])
forall a. Int -> [a] -> ([a], [a])
splitAt (Int -> Int
blockIndexesLen Int
block) [Scientific]
zV1 of
                  ([Scientific]
vals, [Scientific]
zV2) -> [((Int, Int), Scientific)] -> Map (Int, Int) Scientific
symblock ([(Int, Int)] -> [Scientific] -> [((Int, Int), Scientific)]
forall a b. [a] -> [b] -> [(a, b)]
zip (Int -> [(Int, Int)]
blockIndexes Int
block) [Scientific]
vals) Map (Int, Int) Scientific -> Matrix -> Matrix
forall a. a -> [a] -> [a]
: [Int] -> [Scientific] -> Matrix
f [Int]
blocks [Scientific]
zV2
      Matrix
_ -> String -> Solution
forall a. HasCallStack => String -> a
error String
"ToySolver.SDP.transformSolutionBackward: invalid solution"

instance J.ToJSON DualizeInfo where
  toJSON :: DualizeInfo -> Value
toJSON (DualizeInfo Int
origM [Int]
origBlockStruct) =
    [Pair] -> Value
J.object
    [ Key
"type" Key -> Value -> Pair
forall v. ToJSON v => Key -> v -> Pair
forall e kv v. (KeyValue e kv, ToJSON v) => Key -> v -> kv
.= (Value
"DualizeInfo" :: J.Value)
    , Key
"num_original_matrices" Key -> Int -> Pair
forall v. ToJSON v => Key -> v -> Pair
forall e kv v. (KeyValue e kv, ToJSON v) => Key -> v -> kv
.= Int
origM
    , Key
"original_block_structure" Key -> [Int] -> Pair
forall v. ToJSON v => Key -> v -> Pair
forall e kv v. (KeyValue e kv, ToJSON v) => Key -> v -> kv
.= [Int]
origBlockStruct
    ]

instance J.FromJSON DualizeInfo where
  parseJSON :: Value -> Parser DualizeInfo
parseJSON =
    String
-> (Object -> Parser DualizeInfo) -> Value -> Parser DualizeInfo
forall a. String -> (Object -> Parser a) -> Value -> Parser a
withTypedObject String
"DualizeInfo" ((Object -> Parser DualizeInfo) -> Value -> Parser DualizeInfo)
-> (Object -> Parser DualizeInfo) -> Value -> Parser DualizeInfo
forall a b. (a -> b) -> a -> b
$ \Object
obj ->
      Int -> [Int] -> DualizeInfo
DualizeInfo
        (Int -> [Int] -> DualizeInfo)
-> Parser Int -> Parser ([Int] -> DualizeInfo)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Object
obj Object -> Key -> Parser Int
forall a. FromJSON a => Object -> Key -> Parser a
.: Key
"num_original_matrices"
        Parser ([Int] -> DualizeInfo) -> Parser [Int] -> Parser DualizeInfo
forall a b. Parser (a -> b) -> Parser a -> Parser b
forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*> Object
obj Object -> Key -> Parser [Int]
forall a. FromJSON a => Object -> Key -> Parser a
.: Key
"original_block_structure"

symblock :: [((Int,Int), Scientific)] -> SDPFile.Block
symblock :: [((Int, Int), Scientific)] -> Map (Int, Int) Scientific
symblock [((Int, Int), Scientific)]
es = [((Int, Int), Scientific)] -> Map (Int, Int) Scientific
forall k a. Ord k => [(k, a)] -> Map k a
Map.fromList ([((Int, Int), Scientific)] -> Map (Int, Int) Scientific)
-> [((Int, Int), Scientific)] -> Map (Int, Int) Scientific
forall a b. (a -> b) -> a -> b
$ do
  e :: ((Int, Int), Scientific)
e@((Int
i,Int
j),Scientific
x) <- [((Int, Int), Scientific)]
es
  if Scientific
x Scientific -> Scientific -> Bool
forall a. Eq a => a -> a -> Bool
== Scientific
0 then
    []
  else if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j then
    ((Int, Int), Scientific) -> [((Int, Int), Scientific)]
forall a. a -> [a]
forall (m :: * -> *) a. Monad m => a -> m a
return ((Int, Int), Scientific)
e
  else
    [((Int, Int), Scientific)
e, ((Int
j,Int
i),Scientific
x)]