magico 0.0.1.2 → 0.0.2
raw patch · 2 files changed
+129/−133 lines, 2 filesdep +comfort-arraydep +lapackdep −hmatrix
Dependencies added: comfort-array, lapack
Dependencies removed: hmatrix
Files
- magico.cabal +4/−3
- src/Main.hs +125/−130
magico.cabal view
@@ -1,5 +1,5 @@ Name: magico-Version: 0.0.1.2+Version: 0.0.2 Synopsis: Compute solutions for Magico puzzle Description: Compute solutions for Magico puzzle:@@ -14,7 +14,7 @@ Cabal-Version: >=1.10 Source-Repository this- Tag: 0.0.1.2+ Tag: 0.0.2 Type: darcs Location: http://hub.darcs.net/thielema/magico @@ -25,7 +25,8 @@ Executable magico Main-is: Main.hs Build-Depends:- hmatrix >=0.16 && <0.17,+ lapack >=0.2.2 && <0.3,+ comfort-array >=0.3.1 && <0.4, transformers >=0.3 && <0.6, utility-ht >=0.0.11 && <0.1, base >=4.5 && <5
src/Main.hs view
@@ -2,147 +2,138 @@ Magico <http://www.spectra-verlag.de/> -}+{-# LANGUAGE TypeOperators #-} module Main where -import qualified Numeric.Container as NC-import qualified Data.Packed.Matrix as Matrix-import qualified Data.Packed.Vector as Vector-import Numeric.LinearAlgebra.HMatrix (null1, nullspace, rank)-import Numeric.LinearAlgebra ((<\>))-import Data.Packed.Matrix (Matrix)-import Data.Packed.Vector (Vector)+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 Text.Printf (printf)+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 qualified Control.Monad.Trans.State as MS-import Control.Applicative (Applicative, liftA2, liftA3, pure, (<*>))+import Foreign.Storable (Storable)+ import Control.Functor.HT (outerProduct) import qualified Data.Foldable as Fold-import qualified Data.List as List-import Data.Traversable (Traversable, traverse, sequenceA)-import Data.Foldable (Foldable, foldMap)-import Data.Monoid ((<>))+import Data.Tuple.HT (thd3) -data Triple a = Triple a a a deriving (Show)-data Sums a = Sums a a (Triple a) (Triple a) deriving (Show)-type Solution a = Triple (Triple a)+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) -instance Functor Triple where- fmap f (Triple x y z) = Triple (f x) (f y) (f z) -instance Applicative Triple where- pure x = Triple x x x- Triple fx fy fz <*> Triple x y z = Triple (fx x) (fy y) (fz z)+data Index = I0 | I1 | I2 deriving (Eq, Ord, Enum, Bounded, Show)+type Index2 = (Index, Index) -instance Foldable Triple where- foldMap f (Triple x y z) = f x <> f y <> f z+squareShape :: SquareShape+squareShape = (Shape.Enumeration, Shape.Enumeration) -instance Traversable Triple where- traverse f (Triple x y z) = liftA3 Triple (f x) (f y) (f z) +type SumsShape =+ (():+:()) :+: (Shape.Enumeration Index :+: Shape.Enumeration Index) -instance Functor Sums where- fmap f (Sums diag0 diag1 horiz vert) =- Sums (f diag0) (f diag1) (fmap f horiz) (fmap f vert)+sumsShape :: SumsShape+sumsShape = (():+:()) :+: (Shape.Enumeration :+: Shape.Enumeration) -instance Foldable Sums where- foldMap f (Sums diag0 diag1 horiz vert) =- f diag0 <> f diag1 <> foldMap f horiz <> foldMap f vert+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]+ x<+>y = Vector.append x y+ in (Vector.singleton diag0 <+> Vector.singleton diag1)+ <+>+ (vec3 horiz <+> vec3 vert) -getFromList :: MS.State [a] a-getFromList =- MS.state $ \(x:xs) -> (x,xs)--solFromList :: [a] -> Solution a-solFromList =- MS.evalState (sequenceA $ pure $ sequenceA $ pure getFromList)--data Index = I0 | I1 | I2 deriving (Eq, Ord, Enum, Show)-type Index2 = (Index, Index)- indices :: Triple Index-indices = Triple I0 I1 I2--index :: Index -> Triple a -> a-index i (Triple x y z) =- case i of- I0 -> x- I1 -> y- I2 -> z--index2 :: Index2 -> Triple (Triple a) -> a-index2 (i,j) = index j . index i+indices = BoxedArray.indices Shape.Enumeration -flattenIndex2 :: Index2 -> Int-flattenIndex2 (i,j) = 3 * fromEnum i + fromEnum j+triple :: a -> a -> a -> Triple a+triple x y z = BoxedArray.fromList Shape.Enumeration [x,y,z] sumIndices :: Sums (Triple Index2) sumIndices =- Sums- (Triple (I0,I0) (I1,I1) (I2,I2))- (Triple (I0,I2) (I1,I1) (I2,I0))- (outerProduct (,) indices indices)- (outerProduct (flip (,)) indices indices)--sumOnBoard :: (Num a) => Solution a -> Triple Index2 -> a-sumOnBoard sol is = Fold.sum $ fmap (flip index2 sol) is+ 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 +++ [] -sums :: (Num a) => Solution a -> Sums a-sums sol = fmap (sumOnBoard sol) sumIndices+sumsFromSolution :: (Num a, Storable a) => Solution a -> Sums a+sumsFromSolution sol = fmap (Fold.sum . fmap (sol Array.!)) sumIndices -fullMatrix :: Matrix Double+fullMatrix :: Matrix SumsShape SquareShape Double fullMatrix =- Matrix.fromLists $- outerProduct- (\is j -> if Fold.elem j is then 1 else 0)- (Fold.toList sumIndices)- (liftA2 (,) (Fold.toList indices) (Fold.toList indices))--removeAt :: Int -> [a] -> (a, [a])-removeAt k xs =- case splitAt k xs of- (_, []) -> error "removeAt: index too large"- (leftXs, pivot:rightXs) -> (pivot, leftXs++rightXs)--insertAt :: Int -> a -> [a] -> [a]-insertAt k x xs =- case splitAt k xs of- (leftXs, rightXs) -> leftXs ++ x : rightXs-+ Matrix.fromRowArray squareShape $+ fmap+ (\is -> Array.sample squareShape (\j -> if Fold.elem j is then 1 else 0))+ sumIndices -isIntegerVector :: Vector Double -> Bool-isIntegerVector =- all (\x -> abs (x - fromInteger (round x)) < 1e-7) . Vector.toList+isIntegerVector :: Solution Double -> Bool+isIntegerVector x =+ (Vector.normInf $ Vector.sub x $ Array.map (fromInteger . round) x) < 1e-7 --- | it must be @NC.sumElements offset == 0@-integerSolutions :: Int -> Vector Double -> Vector Double -> [Vector Double]+-- | it must be @Vector.sum offset == 0@+integerSolutions ::+ Int -> Solution Double -> Solution Double -> [Solution Double] integerSolutions maxN offset start =- let maxI = NC.maxIndex offset- startI = NC.atIndex start maxI- offsetI = NC.atIndex offset maxI+ let (maxI,offsetI) = Vector.argAbs1Maximum offset+ startI = start!maxI in filter isIntegerVector $ map (\x -> let c = (fromIntegral x - startI) / offsetI- in NC.add start $ NC.scale c offset)+ in Vector.mac c offset start) [0 .. maxN] -splittedMatrix :: Index2 -> ((Int, Vector Double), Matrix Double)-splittedMatrix givenI2 =- let splitPos = flattenIndex2 givenI2- (givenCol, remCols) = removeAt splitPos $ Matrix.toColumns fullMatrix- in ((splitPos, givenCol), Matrix.fromColumns remCols)+extendedMatrix :: Index2 -> ExtMatrix Double+extendedMatrix givenI2 =+ fullMatrix+ ===+ Matrix.singleRow MatrixShape.RowMajor (Vector.unit squareShape givenI2) -ranks :: Triple (Triple Int)-ranks = outerProduct (curry $ rank . snd . splittedMatrix) indices indices -masks :: Triple (Triple (Matrix Double))-masks = outerProduct (curry $ nullspace . snd . splittedMatrix) indices indices+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.ZeroInt 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.ZeroInt 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 @@ -162,59 +153,63 @@ In this case, 'solve' misses some solutions, because it checks only one dimension, not two. -}-solve :: (Index2, Int) -> Sums Int -> [Solution Int]+solve :: (Index2, Int) -> SumsVector Int -> [Solution Int] solve (givenI2, given) ss =- let vec = Vector.fromList $ Fold.toList $ fmap fromIntegral ss- ((splitPos, givenCol), mat) = splittedMatrix givenI2- in map (solFromList . insertAt splitPos given) $- filter (all (>=0)) $- map (map round . Vector.toList) $- integerSolutions (Fold.maximum ss) (null1 mat) $- mat <\> NC.sub vec (NC.scale (fromIntegral given) givenCol)+ let vec =+ Vector.append+ (Array.map fromIntegral ss)+ (Vector.singleton $ fromIntegral given)+ 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) vec printSolution :: Solution Int -> IO () printSolution sol = do- Fold.forM_ sol $- putStrLn . List.intercalate " " . Fold.toList . fmap (printf "%2d")+ Matrix.fromRowMajor (Array.map (fromIntegral :: Int -> Float) sol) ## "%2.f" putStrLn "" exampleA1 :: [Solution Int]-exampleA1 =- solve ((I0,I2), 2) (Sums 5 8 (Triple 5 8 8) (Triple 9 6 6))+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 (Triple 8 9 7) (Triple 6 9 9))+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 (Triple 6 8 8) (Triple 9 7 6))+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 (Triple 8 9 7) (Triple 9 9 6))+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 (Triple 9 6 6) (Triple 9 9 3))+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 (Triple 9 6 8) (Triple 8 6 9))+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 (Triple 10 7 9) (Triple 9 9 8))+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 (Triple 43 40 36) (Triple 37 49 33))+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 (Triple 46 38 39) (Triple 49 35 39))+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