magico-0.0.2.1: src/Main.hs

{- |
Magico
<http://www.spectra-verlag.de/>
-}
{-# LANGUAGE TypeOperators #-}
module Main where

import qualified Numeric.LAPACK.Singular as Singular
import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape
import qualified Numeric.LAPACK.Matrix.Square as Square
import qualified Numeric.LAPACK.Matrix as Matrix
import qualified Numeric.LAPACK.Vector as Vector
import Numeric.LAPACK.Matrix ((===))
import Numeric.LAPACK.Vector (Vector, (+++))
import Numeric.LAPACK.Format ((##))

import qualified Data.Array.Comfort.Boxed as BoxedArray
import qualified Data.Array.Comfort.Storable as Array
import qualified Data.Array.Comfort.Shape as Shape
import Data.Array.Comfort.Storable ((!))
import Data.Array.Comfort.Shape ((:+:))

import Foreign.Storable (Storable)

import Control.Functor.HT (outerProduct)

import qualified Data.Foldable as Fold
import Data.Tuple.HT (thd3)


type Triple = BoxedArray.Array (Shape.Enumeration Index)
type Sums = BoxedArray.Array SumsShape
type Board = BoxedArray.Array SquareShape

type Matrix height width = Matrix.General height width
type ExtMatrix = Matrix (SumsShape:+:()) SquareShape
type SumsVector = Vector SumsShape
type Solution = Vector SquareShape
type SquareShape = (Shape.Enumeration Index, Shape.Enumeration Index)


data Index = I0 | I1 | I2 deriving (Eq, Ord, Enum, Bounded, Show)
type Index2 = (Index, Index)

squareShape :: SquareShape
squareShape = Shape.static


type SumsShape =
      (():+:()) :+: (Shape.Enumeration Index :+: Shape.Enumeration Index)

sumsShape :: SumsShape
sumsShape = Shape.static

sums :: (Storable a) => a -> a -> (a,a,a) -> (a,a,a) -> SumsVector a
sums diag0 diag1 horiz vert =
   let vec3 (x,y,z) = Vector.fromList Shape.Enumeration [x,y,z]
   in (Vector.singleton diag0 +++ Vector.singleton diag1)
      +++
      (vec3 horiz +++ vec3 vert)


indices :: Triple Index
indices = BoxedArray.indices Shape.Enumeration

triple :: a -> a -> a -> Triple a
triple x y z = BoxedArray.fromList Shape.Enumeration [x,y,z]

sumIndices :: Sums (Triple Index2)
sumIndices =
   BoxedArray.fromList sumsShape $
      triple (I0,I0) (I1,I1) (I2,I2) :
      triple (I0,I2) (I1,I1) (I2,I0) :
      outerProduct (,) (Shape.indices Shape.Enumeration) indices ++
      outerProduct (flip (,)) (Shape.indices Shape.Enumeration) indices ++
      []

sumsFromSolution :: (Num a, Storable a) => Solution a -> Sums a
sumsFromSolution sol = fmap (Fold.sum . fmap (sol Array.!)) sumIndices

fullMatrix :: Matrix SumsShape SquareShape Double
fullMatrix =
   Matrix.fromRowArray squareShape $
   fmap
      (Array.fromAssociations 0 squareShape . map (flip (,) 1) . Fold.toList)
      sumIndices

isIntegerVector :: Solution Double -> Bool
isIntegerVector x =
   (Vector.normInf $ Vector.sub x $ Array.map (fromInteger . round) x) < 1e-7

-- | it must be @Vector.sum offset == 0@
integerSolutions ::
   Int -> Solution Double -> Solution Double -> [Solution Double]
integerSolutions maxN offset start =
   let (maxI,offsetI) = Vector.argAbs1Maximum offset
       startI = start!maxI
   in  filter isIntegerVector $
       map
          (\x ->
             let c = (fromIntegral x - startI) / offsetI
             in  Vector.mac c offset start)
          [0 .. maxN]

extendedMatrix :: Index2 -> ExtMatrix Double
extendedMatrix givenI2 =
   fullMatrix
   ===
   Matrix.singleRow MatrixShape.RowMajor (Vector.unit squareShape givenI2)


rank :: ExtMatrix Double -> Int
rank = fst . Singular.pseudoInverseRCond 1e-7

ranks :: Board Int
ranks = BoxedArray.map (rank . extendedMatrix) $ BoxedArray.indices squareShape

nullspace ::
   ExtMatrix Double -> Matrix.General Matrix.ShapeInt SquareShape Double
nullspace mat =
   let (_vr,sigma,vl) = Singular.decompose mat
       rnk = length $ takeWhile (>1e-7) $ Vector.toList sigma
   in Matrix.dropRows rnk $
      Matrix.mapHeight (Shape.ZeroBased . Shape.size) $ Matrix.fromFull vl

masks :: Board (Matrix.General Matrix.ShapeInt SquareShape Double)
masks =
   BoxedArray.map (nullspace . extendedMatrix) $ BoxedArray.indices squareShape

lastRow :: Square.Square SquareShape Double -> Solution Double
lastRow mat = Matrix.takeRow mat (last $ Shape.indices $ Matrix.height mat)

null1 :: ExtMatrix Double -> Solution Double
null1 = lastRow . thd3 . Singular.decompose

{-
Vectors of the null space look like

 0  1 -1
-1  0  1
 1 -1  0

or

-1  2 -1
 0  0  0
 1 -2  1

The second one is the first one added to its horizontally flipped counterpart.

For givenI2 = (I1,I1) several solutions with a center zero can be combined.
In this case, 'solve' misses some solutions,
because it checks only one dimension, not two.
-}
solve :: (Index2, Int) -> SumsVector Int -> [Solution Int]
solve (givenI2, given) ss =
   let mat = extendedMatrix givenI2
   in  filter (all (>=0) . Vector.toList) $
       map (Array.map round) $
       integerSolutions (Fold.maximum $ Vector.toList ss) (null1 mat) $
       Matrix.unliftColumn MatrixShape.ColumnMajor
          (snd . Singular.leastSquaresMinimumNormRCond 1e-7 mat) $
       Array.map fromIntegral $ ss +++ Vector.singleton given

printSolution :: Solution Int -> IO ()
printSolution sol = do
   Matrix.fromRowMajor (Array.map (fromIntegral :: Int -> Float) sol) ## "%2.f"
   putStrLn ""


exampleA1 :: [Solution Int]
exampleA1 = solve ((I0,I2), 2) (sums 5 8 (5,8,8) (9,6,6))

exampleA2 :: [Solution Int]
exampleA2 = solve ((I2,I1), 3) (sums 9 6 (8,9,7) (6,9,9))

exampleA3 :: [Solution Int]
exampleA3 = solve ((I2,I1), 4) (sums 9 4 (6,8,8) (9,7,6))

exampleA4 :: [Solution Int]
exampleA4 = solve ((I1,I2), 4) (sums 9 6 (8,9,7) (9,9,6))

exampleA7 :: [Solution Int]
exampleA7 = solve ((I0,I0), 5) (sums 8 4 (9,6,6) (9,9,3))

exampleA20 :: [Solution Int]
exampleA20 = solve ((I1,I2), 4) (sums 7 7 (9,6,8) (8,6,9))

exampleB1 :: [Solution Int]
exampleB1 = solve ((I0,I1), 4) (sums 11 11 (10,7,9) (9,9,8))

exampleD7 :: [Solution Int]
exampleD7 = solve ((I0,I0), 20) (sums 44 34 (43,40,36) (37,49,33))

exampleD20 :: [Solution Int]
exampleD20 = solve ((I2,I2), 9) (sums 37 43 (46,38,39) (49,35,39))


analysis :: IO ()
analysis = do
   Fold.forM_ masks $ \maskBasis -> do
      Fold.forM_ (Matrix.toRows maskBasis) $ \mask -> do
         Matrix.fromRowMajor mask ## "%6.3f"
         putStrLn ""
      putStrLn ""
   Matrix.fromRowMajor
      (Array.map (fromIntegral :: Int -> Float) (Array.fromBoxed ranks))
         ## "%2.f"

main :: IO ()
main = mapM_ printSolution exampleD7