{- | Magico -} {-# 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