{-# OPTIONS_GHC -Wall #-}
{-# OPTIONS_HADDOCK show-extensions #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE OverloadedStrings #-}
-----------------------------------------------------------------------------
-- |
-- Module      :  ToySolver.QUBO
-- Copyright   :  (c) Masahiro Sakai 2018
-- License     :  BSD-style
--
-- Maintainer  :  masahiro.sakai@gmail.com
-- Stability   :  provisional
-- Portability :  non-portable
--
-----------------------------------------------------------------------------
module ToySolver.QUBO
  ( -- * QUBO (quadratic unconstrained boolean optimization)
    Problem (..)
  , Solution
  , eval

    -- * Ising Model
  , IsingModel (..)
  , evalIsingModel
  ) where

import Control.Monad
import Data.Array.Unboxed
import Data.ByteString.Builder
import Data.ByteString.Builder.Scientific
import qualified Data.ByteString.Lazy.Char8 as BS
import Data.IntMap.Strict (IntMap)
import qualified Data.IntMap.Strict as IntMap
import Data.Scientific
import ToySolver.FileFormat.Base

-- | QUBO (quadratic unconstrained boolean optimization) problem.
--
-- Minimize \(\sum_{i\le j} Q_{i,j} x_i x_j\) where \(x_i \in \{0,1\}\) for \(i \in \{0 \ldots N-1\}\).
--
-- In the `Solution` type. 0 and 1 are represented as @False@ and @True@ respectively.
data Problem a
  = Problem
  { Problem a -> Int
quboNumVars :: !Int
    -- ^ Number of variables N. Variables are numbered from 0 to N-1.
  , Problem a -> IntMap (IntMap a)
quboMatrix :: IntMap (IntMap a)
    -- ^ Upper triangular matrix Q
  }
  deriving (Problem a -> Problem a -> Bool
(Problem a -> Problem a -> Bool)
-> (Problem a -> Problem a -> Bool) -> Eq (Problem a)
forall a. Eq a => Problem a -> Problem a -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: Problem a -> Problem a -> Bool
$c/= :: forall a. Eq a => Problem a -> Problem a -> Bool
== :: Problem a -> Problem a -> Bool
$c== :: forall a. Eq a => Problem a -> Problem a -> Bool
Eq, Int -> Problem a -> ShowS
[Problem a] -> ShowS
Problem a -> String
(Int -> Problem a -> ShowS)
-> (Problem a -> String)
-> ([Problem a] -> ShowS)
-> Show (Problem a)
forall a. Show a => Int -> Problem a -> ShowS
forall a. Show a => [Problem a] -> ShowS
forall a. Show a => Problem a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [Problem a] -> ShowS
$cshowList :: forall a. Show a => [Problem a] -> ShowS
show :: Problem a -> String
$cshow :: forall a. Show a => Problem a -> String
showsPrec :: Int -> Problem a -> ShowS
$cshowsPrec :: forall a. Show a => Int -> Problem a -> ShowS
Show)

instance Functor Problem where
  fmap :: (a -> b) -> Problem a -> Problem b
fmap a -> b
f Problem a
prob =
    Problem :: forall a. Int -> IntMap (IntMap a) -> Problem a
Problem
    { quboNumVars :: Int
quboNumVars = Problem a -> Int
forall a. Problem a -> Int
quboNumVars Problem a
prob
    , quboMatrix :: IntMap (IntMap b)
quboMatrix = (IntMap a -> IntMap b) -> IntMap (IntMap a) -> IntMap (IntMap b)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap ((a -> b) -> IntMap a -> IntMap b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap a -> b
f) (Problem a -> IntMap (IntMap a)
forall a. Problem a -> IntMap (IntMap a)
quboMatrix Problem a
prob)
    }

parseProblem :: (BS.ByteString -> a) -> BS.ByteString -> Either String (Problem a)
parseProblem :: (ByteString -> a) -> ByteString -> Either String (Problem a)
parseProblem ByteString -> a
f ByteString
s =
  case ByteString -> [ByteString]
BS.words ByteString
l of
    [ByteString
"p", ByteString
filetype, ByteString
topology, ByteString
maxNodes, ByteString
_nNodes, ByteString
_nCouplers] ->
      if ByteString
filetype ByteString -> ByteString -> Bool
forall a. Eq a => a -> a -> Bool
/= ByteString
"qubo" then
        String -> Either String (Problem a)
forall a b. a -> Either a b
Left (String -> Either String (Problem a))
-> String -> Either String (Problem a)
forall a b. (a -> b) -> a -> b
$ String
"unknown filetype: " String -> ShowS
forall a. [a] -> [a] -> [a]
++ ByteString -> String
BS.unpack ByteString
filetype
      else if ByteString
topology ByteString -> ByteString -> Bool
forall a. Eq a => a -> a -> Bool
/= ByteString
"0" Bool -> Bool -> Bool
&& ByteString
topology ByteString -> ByteString -> Bool
forall a. Eq a => a -> a -> Bool
/= ByteString
"unconstrained" then
        String -> Either String (Problem a)
forall a b. a -> Either a b
Left (String -> Either String (Problem a))
-> String -> Either String (Problem a)
forall a b. (a -> b) -> a -> b
$ String
"unknown topology: " String -> ShowS
forall a. [a] -> [a] -> [a]
++ ByteString -> String
BS.unpack ByteString
topology
      else
        Problem a -> Either String (Problem a)
forall a b. b -> Either a b
Right (Problem a -> Either String (Problem a))
-> Problem a -> Either String (Problem a)
forall a b. (a -> b) -> a -> b
$ Problem :: forall a. Int -> IntMap (IntMap a) -> Problem a
Problem
        { quboNumVars :: Int
quboNumVars = String -> Int
forall a. Read a => String -> a
read (ByteString -> String
BS.unpack ByteString
maxNodes)
        , quboMatrix :: IntMap (IntMap a)
quboMatrix =
            (IntMap a -> IntMap a -> IntMap a)
-> [IntMap (IntMap a)] -> IntMap (IntMap a)
forall (f :: * -> *) a.
Foldable f =>
(a -> a -> a) -> f (IntMap a) -> IntMap a
IntMap.unionsWith IntMap a -> IntMap a -> IntMap a
forall a. IntMap a -> IntMap a -> IntMap a
IntMap.union ([IntMap (IntMap a)] -> IntMap (IntMap a))
-> [IntMap (IntMap a)] -> IntMap (IntMap a)
forall a b. (a -> b) -> a -> b
$ do
              ByteString
l <- [ByteString]
ls
              case ByteString -> [ByteString]
BS.words ByteString
l of
                [ByteString
i, ByteString
j, ByteString
v] -> IntMap (IntMap a) -> [IntMap (IntMap a)]
forall (m :: * -> *) a. Monad m => a -> m a
return (IntMap (IntMap a) -> [IntMap (IntMap a)])
-> IntMap (IntMap a) -> [IntMap (IntMap a)]
forall a b. (a -> b) -> a -> b
$ Int -> IntMap a -> IntMap (IntMap a)
forall a. Int -> a -> IntMap a
IntMap.singleton (String -> Int
forall a. Read a => String -> a
read (ByteString -> String
BS.unpack ByteString
i)) (IntMap a -> IntMap (IntMap a)) -> IntMap a -> IntMap (IntMap a)
forall a b. (a -> b) -> a -> b
$ Int -> a -> IntMap a
forall a. Int -> a -> IntMap a
IntMap.singleton (String -> Int
forall a. Read a => String -> a
read (ByteString -> String
BS.unpack ByteString
j)) (a -> IntMap a) -> a -> IntMap a
forall a b. (a -> b) -> a -> b
$ ByteString -> a
f ByteString
v
        }
  where
    (ByteString
l:[ByteString]
ls) = (ByteString -> Bool) -> [ByteString] -> [ByteString]
forall a. (a -> Bool) -> [a] -> [a]
filter (Bool -> Bool
not (Bool -> Bool) -> (ByteString -> Bool) -> ByteString -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ByteString -> Bool
isCommentBS) (ByteString -> [ByteString]
BS.lines ByteString
s)

    isCommentBS :: BS.ByteString -> Bool
    isCommentBS :: ByteString -> Bool
isCommentBS ByteString
s =
      case ByteString -> Maybe (Char, ByteString)
BS.uncons ByteString
s of
        Just (Char
'c',ByteString
_) ->  Bool
True
        Maybe (Char, ByteString)
_ -> Bool
False

renderProblem :: (a -> Builder) -> Problem a -> Builder
renderProblem :: (a -> Builder) -> Problem a -> Builder
renderProblem a -> Builder
f Problem a
prob = Builder
header
    Builder -> Builder -> Builder
forall a. Semigroup a => a -> a -> a
<> [Builder] -> Builder
forall a. Monoid a => [a] -> a
mconcat [ Int -> Builder
intDec Int
i Builder -> Builder -> Builder
forall a. Semigroup a => a -> a -> a
<> Char -> Builder
char7 Char
' ' Builder -> Builder -> Builder
forall a. Semigroup a => a -> a -> a
<> Int -> Builder
intDec Int
i Builder -> Builder -> Builder
forall a. Semigroup a => a -> a -> a
<> Char -> Builder
char7 Char
' ' Builder -> Builder -> Builder
forall a. Semigroup a => a -> a -> a
<> a -> Builder
f a
val Builder -> Builder -> Builder
forall a. Semigroup a => a -> a -> a
<> Char -> Builder
char7 Char
'\n'
               | (Int
i,a
val) <- IntMap a -> [(Int, a)]
forall a. IntMap a -> [(Int, a)]
IntMap.toList IntMap a
nodes
               ]
    Builder -> Builder -> Builder
forall a. Semigroup a => a -> a -> a
<> [Builder] -> Builder
forall a. Monoid a => [a] -> a
mconcat [Int -> Builder
intDec Int
i Builder -> Builder -> Builder
forall a. Semigroup a => a -> a -> a
<> Char -> Builder
char7 Char
' ' Builder -> Builder -> Builder
forall a. Semigroup a => a -> a -> a
<> Int -> Builder
intDec Int
j Builder -> Builder -> Builder
forall a. Semigroup a => a -> a -> a
<> Char -> Builder
char7 Char
' ' Builder -> Builder -> Builder
forall a. Semigroup a => a -> a -> a
<> a -> Builder
f a
val Builder -> Builder -> Builder
forall a. Semigroup a => a -> a -> a
<> Char -> Builder
char7 Char
'\n'
               | (Int
i,IntMap a
row) <- IntMap (IntMap a) -> [(Int, IntMap a)]
forall a. IntMap a -> [(Int, a)]
IntMap.toList IntMap (IntMap a)
couplers, (Int
j,a
val) <- IntMap a -> [(Int, a)]
forall a. IntMap a -> [(Int, a)]
IntMap.toList IntMap a
row
               ]
  where
    nodes :: IntMap a
nodes = (Int -> IntMap a -> Maybe a) -> IntMap (IntMap a) -> IntMap a
forall a b. (Int -> a -> Maybe b) -> IntMap a -> IntMap b
IntMap.mapMaybeWithKey Int -> IntMap a -> Maybe a
forall a. Int -> IntMap a -> Maybe a
IntMap.lookup (Problem a -> IntMap (IntMap a)
forall a. Problem a -> IntMap (IntMap a)
quboMatrix Problem a
prob)
    nNodes :: Int
nNodes = IntMap a -> Int
forall a. IntMap a -> Int
IntMap.size IntMap a
nodes
    couplers :: IntMap (IntMap a)
couplers = (Int -> IntMap a -> IntMap a)
-> IntMap (IntMap a) -> IntMap (IntMap a)
forall a b. (Int -> a -> b) -> IntMap a -> IntMap b
IntMap.mapWithKey Int -> IntMap a -> IntMap a
forall a. Int -> IntMap a -> IntMap a
IntMap.delete (Problem a -> IntMap (IntMap a)
forall a. Problem a -> IntMap (IntMap a)
quboMatrix Problem a
prob)
    nCouplers :: Int
nCouplers = [Int] -> Int
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [IntMap a -> Int
forall a. IntMap a -> Int
IntMap.size IntMap a
row | IntMap a
row <- IntMap (IntMap a) -> [IntMap a]
forall a. IntMap a -> [a]
IntMap.elems IntMap (IntMap a)
couplers]
    header :: Builder
header = [Builder] -> Builder
forall a. Monoid a => [a] -> a
mconcat
      [Builder
"p qubo 0 "
      , Int -> Builder
intDec (Problem a -> Int
forall a. Problem a -> Int
quboNumVars Problem a
prob), Char -> Builder
char7 Char
' '
      , Int -> Builder
intDec Int
nNodes, Char -> Builder
char7 Char
' '
      , Int -> Builder
intDec Int
nCouplers, Char -> Builder
char7 Char
'\n'
      ]

instance FileFormat (Problem Scientific) where
  parse :: ByteString -> Either String (Problem Scientific)
parse = (ByteString -> Scientific)
-> ByteString -> Either String (Problem Scientific)
forall a.
(ByteString -> a) -> ByteString -> Either String (Problem a)
parseProblem (String -> Scientific
forall a. Read a => String -> a
read (String -> Scientific)
-> (ByteString -> String) -> ByteString -> Scientific
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ByteString -> String
BS.unpack)
  render :: Problem Scientific -> Builder
render = (Scientific -> Builder) -> Problem Scientific -> Builder
forall a. (a -> Builder) -> Problem a -> Builder
renderProblem Scientific -> Builder
scientificBuilder


type Solution = UArray Int Bool

eval :: Num a => Solution -> Problem a -> a
eval :: Solution -> Problem a -> a
eval Solution
sol Problem a
prob = [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ([a] -> a) -> [a] -> a
forall a b. (a -> b) -> a -> b
$ do
  (Int
x1, IntMap a
row) <- IntMap (IntMap a) -> [(Int, IntMap a)]
forall a. IntMap a -> [(Int, a)]
IntMap.toList (IntMap (IntMap a) -> [(Int, IntMap a)])
-> IntMap (IntMap a) -> [(Int, IntMap a)]
forall a b. (a -> b) -> a -> b
$ Problem a -> IntMap (IntMap a)
forall a. Problem a -> IntMap (IntMap a)
quboMatrix Problem a
prob
  Bool -> [()]
forall (f :: * -> *). Alternative f => Bool -> f ()
guard (Bool -> [()]) -> Bool -> [()]
forall a b. (a -> b) -> a -> b
$ Solution
sol Solution -> Int -> Bool
forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
a i e -> i -> e
! Int
x1
  (Int
x2, a
c) <- IntMap a -> [(Int, a)]
forall a. IntMap a -> [(Int, a)]
IntMap.toList IntMap a
row
  Bool -> [()]
forall (f :: * -> *). Alternative f => Bool -> f ()
guard (Bool -> [()]) -> Bool -> [()]
forall a b. (a -> b) -> a -> b
$ Solution
sol Solution -> Int -> Bool
forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
a i e -> i -> e
! Int
x2
  a -> [a]
forall (m :: * -> *) a. Monad m => a -> m a
return a
c


-- | Ising model.
--
-- Minimize \(\sum_{i<j} J_{i,j} \sigma_i \sigma_j + \sum_i h_i \sigma_i\) where \(\sigma_i \in \{-1,+1\}\) for \(i \in \{0 \ldots N-1\}\).
--
-- In the `Solution` type. -1 and +1 are represented as @False@ and @True@ respectively.
data IsingModel a
  = IsingModel
  { IsingModel a -> Int
isingNumVars :: !Int
    -- ^ Number of variables N. Variables are numbered from 0 to N-1.
  , IsingModel a -> IntMap (IntMap a)
isingInteraction :: IntMap (IntMap a)
    -- ^ Interaction \(J_{i,j}\) with \(i < j\).
  , IsingModel a -> IntMap a
isingExternalMagneticField :: IntMap a
    -- ^ External magnetic field \(h_j\).
  }
  deriving (IsingModel a -> IsingModel a -> Bool
(IsingModel a -> IsingModel a -> Bool)
-> (IsingModel a -> IsingModel a -> Bool) -> Eq (IsingModel a)
forall a. Eq a => IsingModel a -> IsingModel a -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: IsingModel a -> IsingModel a -> Bool
$c/= :: forall a. Eq a => IsingModel a -> IsingModel a -> Bool
== :: IsingModel a -> IsingModel a -> Bool
$c== :: forall a. Eq a => IsingModel a -> IsingModel a -> Bool
Eq, Int -> IsingModel a -> ShowS
[IsingModel a] -> ShowS
IsingModel a -> String
(Int -> IsingModel a -> ShowS)
-> (IsingModel a -> String)
-> ([IsingModel a] -> ShowS)
-> Show (IsingModel a)
forall a. Show a => Int -> IsingModel a -> ShowS
forall a. Show a => [IsingModel a] -> ShowS
forall a. Show a => IsingModel a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [IsingModel a] -> ShowS
$cshowList :: forall a. Show a => [IsingModel a] -> ShowS
show :: IsingModel a -> String
$cshow :: forall a. Show a => IsingModel a -> String
showsPrec :: Int -> IsingModel a -> ShowS
$cshowsPrec :: forall a. Show a => Int -> IsingModel a -> ShowS
Show)

evalIsingModel :: Num a => Solution -> IsingModel a -> a
evalIsingModel :: Solution -> IsingModel a -> a
evalIsingModel Solution
sol IsingModel a
m
  = [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [ a
jj_ij a -> a -> a
forall a. Num a => a -> a -> a
* Int -> a
forall p. Num p => Int -> p
sigma Int
i a -> a -> a
forall a. Num a => a -> a -> a
*  Int -> a
forall p. Num p => Int -> p
sigma Int
j
        | (Int
i, IntMap a
row) <- IntMap (IntMap a) -> [(Int, IntMap a)]
forall a. IntMap a -> [(Int, a)]
IntMap.toList (IntMap (IntMap a) -> [(Int, IntMap a)])
-> IntMap (IntMap a) -> [(Int, IntMap a)]
forall a b. (a -> b) -> a -> b
$ IsingModel a -> IntMap (IntMap a)
forall a. IsingModel a -> IntMap (IntMap a)
isingInteraction IsingModel a
m, (Int
j, a
jj_ij) <- IntMap a -> [(Int, a)]
forall a. IntMap a -> [(Int, a)]
IntMap.toList IntMap a
row
        ]
  a -> a -> a
forall a. Num a => a -> a -> a
+ [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [ a
h_i a -> a -> a
forall a. Num a => a -> a -> a
* Int -> a
forall p. Num p => Int -> p
sigma Int
i | (Int
i, a
h_i) <- IntMap a -> [(Int, a)]
forall a. IntMap a -> [(Int, a)]
IntMap.toList (IntMap a -> [(Int, a)]) -> IntMap a -> [(Int, a)]
forall a b. (a -> b) -> a -> b
$ IsingModel a -> IntMap a
forall a. IsingModel a -> IntMap a
isingExternalMagneticField IsingModel a
m ]
  where
    sigma :: Int -> p
sigma Int
i = if Solution
sol Solution -> Int -> Bool
forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
a i e -> i -> e
! Int
i then p
1 else -p
1