{-# 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
  { forall a. Problem a -> Int
quboNumVars :: !Int
    -- ^ Number of variables N. Variables are numbered from 0 to N-1.
  , forall a. Problem a -> IntMap (IntMap a)
quboMatrix :: IntMap (IntMap a)
    -- ^ Upper triangular matrix Q
  }
  deriving (Problem a -> Problem a -> Bool
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
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 :: forall a b. (a -> b) -> Problem a -> Problem b
fmap a -> b
f Problem a
prob =
    Problem
    { quboNumVars :: Int
quboNumVars = forall a. Problem a -> Int
quboNumVars Problem a
prob
    , quboMatrix :: IntMap (IntMap b)
quboMatrix = forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap a -> b
f) (forall a. Problem a -> IntMap (IntMap a)
quboMatrix Problem a
prob)
    }

parseProblem :: (BS.ByteString -> a) -> BS.ByteString -> Either String (Problem a)
parseProblem :: forall a.
(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 forall a. Eq a => a -> a -> Bool
/= ByteString
"qubo" then
        forall a b. a -> Either a b
Left forall a b. (a -> b) -> a -> b
$ String
"unknown filetype: " forall a. [a] -> [a] -> [a]
++ ByteString -> String
BS.unpack ByteString
filetype
      else if ByteString
topology forall a. Eq a => a -> a -> Bool
/= ByteString
"0" Bool -> Bool -> Bool
&& ByteString
topology forall a. Eq a => a -> a -> Bool
/= ByteString
"unconstrained" then
        forall a b. a -> Either a b
Left forall a b. (a -> b) -> a -> b
$ String
"unknown topology: " forall a. [a] -> [a] -> [a]
++ ByteString -> String
BS.unpack ByteString
topology
      else
        forall a b. b -> Either a b
Right forall a b. (a -> b) -> a -> b
$ Problem
        { quboNumVars :: Int
quboNumVars = forall a. Read a => String -> a
read (ByteString -> String
BS.unpack ByteString
maxNodes)
        , quboMatrix :: IntMap (IntMap a)
quboMatrix =
            forall (f :: * -> *) a.
Foldable f =>
(a -> a -> a) -> f (IntMap a) -> IntMap a
IntMap.unionsWith forall a. IntMap a -> IntMap a -> IntMap a
IntMap.union 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] -> forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ forall a. Int -> a -> IntMap a
IntMap.singleton (forall a. Read a => String -> a
read (ByteString -> String
BS.unpack ByteString
i)) forall a b. (a -> b) -> a -> b
$ forall a. Int -> a -> IntMap a
IntMap.singleton (forall a. Read a => String -> a
read (ByteString -> String
BS.unpack ByteString
j)) forall a b. (a -> b) -> a -> b
$ ByteString -> a
f ByteString
v
        }
  where
    (ByteString
l:[ByteString]
ls) = forall a. (a -> Bool) -> [a] -> [a]
filter (Bool -> Bool
not 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 :: forall a. (a -> Builder) -> Problem a -> Builder
renderProblem a -> Builder
f Problem a
prob = Builder
header
    forall a. Semigroup a => a -> a -> a
<> forall a. Monoid a => [a] -> a
mconcat [ Int -> Builder
intDec Int
i forall a. Semigroup a => a -> a -> a
<> Char -> Builder
char7 Char
' ' forall a. Semigroup a => a -> a -> a
<> Int -> Builder
intDec Int
i forall a. Semigroup a => a -> a -> a
<> Char -> Builder
char7 Char
' ' forall a. Semigroup a => a -> a -> a
<> a -> Builder
f a
val forall a. Semigroup a => a -> a -> a
<> Char -> Builder
char7 Char
'\n'
               | (Int
i,a
val) <- forall a. IntMap a -> [(Int, a)]
IntMap.toList IntMap a
nodes
               ]
    forall a. Semigroup a => a -> a -> a
<> forall a. Monoid a => [a] -> a
mconcat [Int -> Builder
intDec Int
i forall a. Semigroup a => a -> a -> a
<> Char -> Builder
char7 Char
' ' forall a. Semigroup a => a -> a -> a
<> Int -> Builder
intDec Int
j forall a. Semigroup a => a -> a -> a
<> Char -> Builder
char7 Char
' ' forall a. Semigroup a => a -> a -> a
<> a -> Builder
f a
val forall a. Semigroup a => a -> a -> a
<> Char -> Builder
char7 Char
'\n'
               | (Int
i,IntMap a
row) <- forall a. IntMap a -> [(Int, a)]
IntMap.toList IntMap (IntMap a)
couplers, (Int
j,a
val) <- forall a. IntMap a -> [(Int, a)]
IntMap.toList IntMap a
row
               ]
  where
    nodes :: IntMap a
nodes = forall a b. (Int -> a -> Maybe b) -> IntMap a -> IntMap b
IntMap.mapMaybeWithKey forall a. Int -> IntMap a -> Maybe a
IntMap.lookup (forall a. Problem a -> IntMap (IntMap a)
quboMatrix Problem a
prob)
    nNodes :: Int
nNodes = forall a. IntMap a -> Int
IntMap.size IntMap a
nodes
    couplers :: IntMap (IntMap a)
couplers = forall a b. (Int -> a -> b) -> IntMap a -> IntMap b
IntMap.mapWithKey forall a. Int -> IntMap a -> IntMap a
IntMap.delete (forall a. Problem a -> IntMap (IntMap a)
quboMatrix Problem a
prob)
    nCouplers :: Int
nCouplers = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [forall a. IntMap a -> Int
IntMap.size IntMap a
row | IntMap a
row <- forall a. IntMap a -> [a]
IntMap.elems IntMap (IntMap a)
couplers]
    header :: Builder
header = forall a. Monoid a => [a] -> a
mconcat
      [Builder
"p qubo 0 "
      , Int -> Builder
intDec (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 = forall a.
(ByteString -> a) -> ByteString -> Either String (Problem a)
parseProblem (forall a. Read a => String -> a
read forall b c a. (b -> c) -> (a -> b) -> a -> c
. ByteString -> String
BS.unpack)
  render :: Problem Scientific -> Builder
render = forall a. (a -> Builder) -> Problem a -> Builder
renderProblem Scientific -> Builder
scientificBuilder


type Solution = UArray Int Bool

eval :: Num a => Solution -> Problem a -> a
eval :: forall a. Num a => Solution -> Problem a -> a
eval Solution
sol Problem a
prob = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum forall a b. (a -> b) -> a -> b
$ do
  (Int
x1, IntMap a
row) <- forall a. IntMap a -> [(Int, a)]
IntMap.toList forall a b. (a -> b) -> a -> b
$ forall a. Problem a -> IntMap (IntMap a)
quboMatrix Problem a
prob
  forall (f :: * -> *). Alternative f => Bool -> f ()
guard forall a b. (a -> b) -> a -> b
$ Solution
sol forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
a i e -> i -> e
! Int
x1
  (Int
x2, a
c) <- forall a. IntMap a -> [(Int, a)]
IntMap.toList IntMap a
row
  forall (f :: * -> *). Alternative f => Bool -> f ()
guard forall a b. (a -> b) -> a -> b
$ Solution
sol forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
a i e -> i -> e
! Int
x2
  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
  { forall a. IsingModel a -> Int
isingNumVars :: !Int
    -- ^ Number of variables N. Variables are numbered from 0 to N-1.
  , forall a. IsingModel a -> IntMap (IntMap a)
isingInteraction :: IntMap (IntMap a)
    -- ^ Interaction \(J_{i,j}\) with \(i < j\).
  , forall a. IsingModel a -> IntMap a
isingExternalMagneticField :: IntMap a
    -- ^ External magnetic field \(h_j\).
  }
  deriving (IsingModel a -> IsingModel a -> Bool
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
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 :: forall a. Num a => Solution -> IsingModel a -> a
evalIsingModel Solution
sol IsingModel a
m
  = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [ a
jj_ij forall a. Num a => a -> a -> a
* forall {a}. Num a => Int -> a
sigma Int
i forall a. Num a => a -> a -> a
*  forall {a}. Num a => Int -> a
sigma Int
j
        | (Int
i, IntMap a
row) <- forall a. IntMap a -> [(Int, a)]
IntMap.toList forall a b. (a -> b) -> a -> b
$ forall a. IsingModel a -> IntMap (IntMap a)
isingInteraction IsingModel a
m, (Int
j, a
jj_ij) <- forall a. IntMap a -> [(Int, a)]
IntMap.toList IntMap a
row
        ]
  forall a. Num a => a -> a -> a
+ forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [ a
h_i forall a. Num a => a -> a -> a
* forall {a}. Num a => Int -> a
sigma Int
i | (Int
i, a
h_i) <- forall a. IntMap a -> [(Int, a)]
IntMap.toList forall a b. (a -> b) -> a -> b
$ forall a. IsingModel a -> IntMap a
isingExternalMagneticField IsingModel a
m ]
  where
    sigma :: Int -> a
sigma Int
i = if Solution
sol forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
a i e -> i -> e
! Int
i then a
1 else -a
1