{-# OPTIONS_GHC -Wall #-}
{-# LANGUAGE TypeFamilies #-}
module ToySolver.SDP
( dualize
, DualizeInfo (..)
) where
import qualified Data.Map.Strict as Map
import Data.Scientific (Scientific)
import ToySolver.Converter.Base
import qualified ToySolver.Text.SDPFile as SDPFile
dualize :: SDPFile.Problem -> (SDPFile.Problem, DualizeInfo)
dualize origProb =
( SDPFile.Problem
{ SDPFile.blockStruct = blockStruct
, SDPFile.costs = d
, SDPFile.matrices = a0 : as
}
, DualizeInfo m (SDPFile.blockStruct origProb)
)
where
m :: Int
m = SDPFile.mDim origProb
c :: [Scientific]
c = SDPFile.costs origProb
f0 :: SDPFile.Matrix
fs :: [SDPFile.Matrix]
f0:fs = SDPFile.matrices origProb
blockStruct :: [Int]
blockStruct = [-m, -m] ++ SDPFile.blockStruct origProb
a0 :: SDPFile.Matrix
a0 =
[ Map.fromList [((j,j), -cj) | (j,cj) <- zip [1..m] c, cj /= 0]
, Map.fromList [((j,j), cj) | (j,cj) <- zip [1..m] c, cj /= 0]
] ++
[ Map.empty | _ <- SDPFile.blockStruct origProb]
as :: [SDPFile.Matrix]
as =
[ [ Map.fromList [ ((k,k), - (if i == j then v else 2*v))
| (k,fk) <- zip [1..m] fs, let v = SDPFile.blockElem i j (fk!!b), v /= 0]
, Map.fromList [ ((k,k), (if i == j then v else 2*v))
| (k,fk) <- zip [1..m] fs, let v = SDPFile.blockElem i j (fk!!b), v /= 0]
] ++
[ if b /= b2 then
Map.empty
else if i == j then
Map.singleton (i,j) 1
else
Map.fromList [((i,j),1), ((j,i),1)]
| (b2, _) <- zip [0..] (SDPFile.blockStruct origProb)
]
| (b,block) <- zip [0..] (SDPFile.blockStruct origProb)
, (i,j) <- blockIndexes block
]
d =
[ - (if i == j then v else 2*v)
| (b,block) <- zip [0..] (SDPFile.blockStruct origProb)
, (i,j) <- blockIndexes block
, let v = SDPFile.blockElem i j (f0 !! b)
]
blockIndexes :: Int -> [(Int,Int)]
blockIndexes n = if n >= 0 then [(i,j) | i <- [1..n], j <- [1..i]] else [(i,i) | i <- [1..(-n)]]
blockIndexesLen :: Int -> Int
blockIndexesLen n = if n >= 0 then n*(n+1) `div` 2 else -n
data DualizeInfo = DualizeInfo Int [Int]
deriving (Eq, Show, Read)
instance Transformer DualizeInfo where
type Source DualizeInfo = SDPFile.Solution
type Target DualizeInfo = SDPFile.Solution
instance ForwardTransformer DualizeInfo where
transformForward (DualizeInfo _origM origBlockStruct)
SDPFile.Solution
{ SDPFile.primalVector = xV
, SDPFile.primalMatrix = xM
, SDPFile.dualMatrix = yM
} =
SDPFile.Solution
{ SDPFile.primalVector = zV
, SDPFile.primalMatrix = zM
, SDPFile.dualMatrix = wM
}
where
zV = concat [[SDPFile.blockElem i j block | (i,j) <- blockIndexes b] | (b,block) <- zip origBlockStruct yM]
zM = Map.empty : Map.empty : yM
wM =
[ Map.fromList $ zipWith (\i x -> ((i,i), if x >= 0 then x else 0)) [1..] xV
, Map.fromList $ zipWith (\i x -> ((i,i), if x <= 0 then -x else 0)) [1..] xV
] ++ xM
instance BackwardTransformer DualizeInfo where
transformBackward (DualizeInfo origM origBlockStruct)
SDPFile.Solution
{ SDPFile.primalVector = zV
, SDPFile.primalMatrix = _zM
, SDPFile.dualMatrix = wM
} =
case wM of
(xps:xns:xM) ->
SDPFile.Solution
{ SDPFile.primalVector = xV
, SDPFile.primalMatrix = xM
, SDPFile.dualMatrix = yM
}
where
xV = [SDPFile.blockElem i i xps - SDPFile.blockElem i i xns | i <- [1..origM]]
yM = f origBlockStruct zV
where
f [] _ = []
f (block : blocks) zV1 =
case splitAt (blockIndexesLen block) zV1 of
(vals, zV2) -> symblock (zip (blockIndexes block) vals) : f blocks zV2
_ -> error "ToySolver.SDP.transformSolutionBackward: invalid solution"
symblock :: [((Int,Int), Scientific)] -> SDPFile.Block
symblock es = Map.fromList $ do
e@((i,j),x) <- es
if x == 0 then
[]
else if i == j then
return e
else
[e, ((j,i),x)]