wuerfelschlange (empty) → 0.1
raw patch · 4 files changed
+448/−0 lines, 4 filesdep +basedep +comfort-arraydep +containers
Dependencies added: base, comfort-array, containers, lapack, netlib-ffi, utility-ht
Files
- LICENSE +27/−0
- src/Graph.hs +93/−0
- src/Matrix.hs +268/−0
- wuerfelschlange.cabal +60/−0
+ LICENSE view
@@ -0,0 +1,27 @@+Copyright (c) Henning Thielemann 2018++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions+are met:+1. Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.+2. Redistributions in binary form must reproduce the above copyright+ notice, this list of conditions and the following disclaimer in the+ documentation and/or other materials provided with the distribution.+3. Neither the name of the author nor the names of his contributors+ may be used to endorse or promote products derived from this software+ without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE+ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS+OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)+HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT+LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY+OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF+SUCH DAMAGE.
+ src/Graph.hs view
@@ -0,0 +1,93 @@+module Main where++import Text.Printf (printf)++import qualified Data.Map as Map+import Data.Map (Map)+++type Node = (Int, Int)++pattern :: Node -> [Bool]+pattern (a,b) = map (flip elem [a,b]) [0..5]++ident, label :: Node -> String+ident (a,b) =+ map (\x -> if x then 'A' else '_') $ pattern (a,b)++label (a,b) =+ map (\x -> if x then ':' else '.') $ pattern (a,b)++move :: [((Node, Node), String)]+move =+ (((0,0), (0,0)), "loop above") :+ [(((a,b),(a-1,b-1)), "") | a<-[1..5], b<-[a..5]]++roll :: [((Node, Node), String)]+roll =+ [(((0,1),(0,0)), "")]+ +++ [(((0,1),(0,1)), "loop above")]+ +++ [(((0,1),(0,a)), "out=-80,in=170") | a<-[2..5]]+ +++ [((from,to), Map.findWithDefault "" (from,to) arcs) |+ b<-[2..5], let b' = pred b, a<-[0..5],+ let from = (0,b); to = (min a b', max a b')]++arcs :: Map (Node, Node) String+arcs =+ Map.fromList $+ (((0,2),(1,3)), "out=-10,in=-170") :+ (((0,3),(2,3)), "out=45,in=-80") :+ (((0,3),(2,5)), "out=-20,in=-160") :+ (((0,4),(3,4)), "out=45,in=-80") :+ (((0,5),(3,4)), "out=90,in=-80") :+ (((0,5),(4,5)), "out=45,in=-80") :+ []++edge1 :: ((Node, Node), String) -> String+edge1 (((a,b),(c,d)), option) =+ printf "(%s) edge [%s,dashed] (%s)" (ident (a,b)) option (ident (c,d))++edge2 :: ((Node, Node), String) -> String+edge2 (((a,b),(c,d)), option) =+ printf "(%s) edge [%s] (%s)" (ident (a,b)) option (ident (c,d))++pos :: Node -> String+pos (a,b) =+ printf "\\node[state%s] (%s) [at={(%d/2+%d*2,1.5*(%d-%d))}] {%s};"+ (if a==0 && b/=0 then ",rolling" else "")+ (ident (a,b)) a b a b (label (a,b))++posCircularSector :: Node -> String+posCircularSector p@(0,0) =+ printf "\\node[state] (%s) [at={(0,0)}] {%s};" (ident p) (label p)+posCircularSector (a,b) =+ let af = fromIntegral a+ bf = fromIntegral b+ cos90, sin90 :: Float -> Float+ cos90 x = cos (x/180*pi)+ sin90 x = sin (x/180*pi)+ in printf "\\node[state] (%s) [at={(%f,%f)}] {%s};"+ (ident (a,b))+ (2*bf*cos90(60/bf*(af-bf)))+ (2*bf*sin90(60/bf*(af-bf)))+ (label (a,b))++allNodes :: [Node]+allNodes = [(a,b) | a<-[0..5],b<-[a..5]]++graph :: String+graph =+ unlines $+ "\\begin{tikzpicture}[auto,rolling/.style=double]" :+ map pos allNodes +++ "\\path[->]" :+ map edge1 move ++ map edge2 roll +++ ";" :+ "\\end{tikzpicture}" :+ []++main :: IO ()+main = putStr graph
+ src/Matrix.hs view
@@ -0,0 +1,268 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}+module Main where++import qualified Numeric.LAPACK.Orthogonal as Ortho+import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape+import qualified Numeric.LAPACK.Matrix.Layout as Layout+import qualified Numeric.LAPACK.Matrix.Square as Square+import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix+import qualified Numeric.LAPACK.Matrix as Matrix+import qualified Numeric.LAPACK.Vector as Vector+import qualified Numeric.LAPACK.Shape as ExtShape+import qualified Numeric.Netlib.Class as Class+import Numeric.LAPACK.Matrix+ (ShapeInt,+ (#-#), (#*#), (#\##), (#*##), (##*#), (#*|), (#*\), (\*#), (|||))++import qualified Data.Array.Comfort.Storable as Array+import qualified Data.Array.Comfort.Boxed as BoxedArray+import qualified Data.Array.Comfort.Shape as Shape+import Data.Array.Comfort.Shape ((::+))+import Data.Array.Comfort.Storable ((!))++import qualified Data.Complex as Complex+import Data.Complex (Complex)++import Text.Printf (printf)++import qualified Data.List.Key as Key+import qualified Data.List as List+import Data.Semigroup ((<>))+import Data.Tuple.HT (swap)+++type ShapeInt2 = ShapeInt ::+ ShapeInt++data Score = P1 | P2 | P3 | P4 | P5 | P6 deriving (Eq, Ord, Enum, Bounded, Show)+type ScoreShape = Shape.Enumeration Score+type TriShape = Shape.UpperTriangular ScoreShape+type TriShapeInt = ExtShape.IntIndexed TriShape+type Transition = Matrix.Square TriShape+type Transformation = Matrix.Square TriShapeInt+type TallMatrix = Matrix.Tall TriShape ShapeInt+type Matrix = Matrix.General TriShape TriShape+type Vector = Vector.Vector TriShape+type Eigenvalues = Vector.Vector TriShapeInt+++triShape :: TriShape+triShape = Shape.upperTriangular Shape.Enumeration++transition :: Transition Double+transition =+ Square.fromFull $ Matrix.scale (1/6) $+ Matrix.fromRowArray triShape $ BoxedArray.fromList triShape $+ map (Vector.fromList triShape) $+ [6,1,0,0,0,0,6,0,0,0,0,0,0,0,0,0,0,0,0,0,0] :+ [0,1,1,0,0,0,0,6,0,0,0,0,0,0,0,0,0,0,0,0,0] :+ [0,1,0,1,0,0,0,0,6,0,0,0,0,0,0,0,0,0,0,0,0] :+ [0,1,0,0,1,0,0,0,0,6,0,0,0,0,0,0,0,0,0,0,0] :+ [0,1,0,0,0,1,0,0,0,0,6,0,0,0,0,0,0,0,0,0,0] :+ [0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0] :+ [0,0,1,0,0,0,0,0,0,0,0,6,0,0,0,0,0,0,0,0,0] :+ [0,0,1,1,0,0,0,0,0,0,0,0,6,0,0,0,0,0,0,0,0] :+ [0,0,1,0,1,0,0,0,0,0,0,0,0,6,0,0,0,0,0,0,0] :+ [0,0,1,0,0,1,0,0,0,0,0,0,0,0,6,0,0,0,0,0,0] :+ [0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0] :+ [0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,6,0,0,0,0,0] :+ [0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,6,0,0,0,0] :+ [0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,6,0,0,0] :+ [0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0] :+ [0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,6,0,0] :+ [0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,6,0] :+ [0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0] :+ [0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,6] :+ [0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0] :+ [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0] :+ []++-- cf. Graph.move, roll+transitionAlt :: Transition Double+transitionAlt =+ Matrix.scale (1/6) $ Matrix.transpose $ ArrMatrix.fromVector $+ Array.fromAssociations 0 (Layout.square MatrixShape.RowMajor triShape) $+ ((((P1,P1),(P1,P1)), 6) :+ [(((a,b),(pred a, pred b)), 6) | a<-[P2 .. P6], b<-[a .. P6]] +++ [((from,to), 1) |+ b<-[P2 .. P6], let b' = pred b, a<-[P1 .. P6],+ let from = (P1,b); to = (min a b', max a b')])++initial1 :: Vector Double+initial1 =+ Vector.scale (1/6) $+ Vector.fromList triShape [1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]++initial2 :: Vector Double+initial2 =+ Vector.scale (1/36) $+ Vector.fromList triShape [1,2,2,2,2,2,1,2,2,2,2,1,2,2,2,1,2,2,1,2,1]++eigenSystem ::+ (a ~ Complex Double) =>+ (Transformation a, Eigenvalues a, Transformation a)+eigenSystem =+ let (vr,d,vlAdj) =+ Square.eigensystem $ Square.mapSize ExtShape.IntIndexed transition+ scal = Square.takeDiagonal $ vlAdj <> vr+ in (vr, d, Vector.recip scal \*# vlAdj)++similarityTransformation ::+ (Class.Floating a, Vector.RealOf a ~ ar, Class.Real ar) =>+ (Transformation a, Eigenvalues a, Transformation a) -> Transition a+similarityTransformation (trafo, diag, trafoInv) =+ Square.mapSize ExtShape.deconsIntIndexed $ trafo #*\ diag ##*# trafoInv++reconstruct :: Transition (Complex Double)+reconstruct = similarityTransformation eigenSystem++reconstructionError :: (Complex Double, ((Score,Score), (Score,Score)))+reconstructionError =+ swap $ Vector.argAbsMaximum $ ArrMatrix.toVector $+ Matrix.sub (Matrix.fromReal transition) reconstruct+++unwrapEigenvector :: Vector.Vector TriShapeInt a -> Vector a+unwrapEigenvector = Array.mapShape ExtShape.deconsIntIndexed++eigenvectors :: [(Complex Double, Vector (Complex Double))]+eigenvectors =+ let (m,v,_minv) = eigenSystem+ in Key.sort (negate . Complex.magnitude . fst) $+ zip (Vector.toList v) (map unwrapEigenvector $ Matrix.toColumns m)++tallMatrixFromColumns :: [Vector Double] -> TallMatrix Double+tallMatrixFromColumns =+ either id (error "matrix not tall") . Matrix.caseTallWide .+ Matrix.fromColumns (Matrix.height transition)++dominantEigenmatrix :: TallMatrix Double+dominantEigenmatrix =+ tallMatrixFromColumns $+ map (Array.map Complex.realPart . snd) $ take 2 eigenvectors++{-+This is like the first step of the Schur decomposition.++https://cs.stackexchange.com/questions/33068/partially-diagonalizing-a-matrix+https://math.stackexchange.com/questions/3967036/partial-eigendecomposition-of-a-positive-semi-definite-matrix+-}+partialEigenSystem ::+ (Matrix.LiberalSquare TriShape ShapeInt2 Double,+ Matrix.LiberalSquare TriShape ShapeInt2 Double)+partialEigenSystem =+ let (ml,mr) = dominantEigenmatrixPairs+ scal = Square.takeDiagonal $ Square.fromFull $ Matrix.adjoint ml #*# mr+ mlScal = Matrix.generalizeTall $ ml #*\ Vector.recip scal+ ortho = Matrix.generalizeTall $ Ortho.complement ml+ in (Square.liberalFromFull $+ mlScal ||| ortho #-# mlScal #*# (Matrix.adjoint mr #*# ortho),+ Square.liberalFromFull $ Matrix.generalizeTall mr ||| ortho)++partiallyDiagonalized :: Matrix.LiberalSquare ShapeInt2 ShapeInt2 Double+partiallyDiagonalized =+ case fromInteger 1 :: Int of+ 0 -> let mr = snd partialEigenSystem in mr #\## transition #*## mr+ _ ->+ let (ml,mr) = partialEigenSystem+ in Matrix.adjoint ml #*## transition #*## mr++{-+U are the given eigenvectors.+Uorth is the basis for their orthogonal complement.++U _|_ Uorth+M*U _|_ M*Uorth - does not apply++Thus, (U | Uorth) cannot diagonalize the transition matrix.+-}+++filterEig :: (Double -> Bool) -> Transition (Complex Double)+filterEig p =+ let (trafo, diag, trafoInv) = eigenSystem+ in similarityTransformation+ (trafo,+ Array.map (\x -> if p (Complex.magnitude x) then 1 else 0) diag,+ trafoInv)++powerCoeff :: Complex Double+powerCoeff =+ (filterEig (\x->0.99<x && x<1.01) #*| Vector.fromReal initial2) ! (P1,P1)+++outerProduct ::+ (Eq height, Eq width, Shape.C height, Shape.C width, Class.Floating a) =>+ Vector.Vector height a ->+ Vector.Vector width a ->+ Matrix.General height width a+outerProduct = Matrix.outer MatrixShape.RowMajor++rank2matrix :: Matrix.General ShapeInt ShapeInt Double+rank2matrix =+ Matrix.add+ (Vector.autoFromList [38,37,20,80,38]+ `outerProduct`+ Vector.autoFromList [97,30,19,49,37])+ (Vector.autoFromList [7,43,89,64,20]+ `outerProduct`+ Vector.autoFromList [4,96,44,49,80])++eigenvectorPairs ::+ [(Complex Double, (Vector (Complex Double), Vector (Complex Double)))]+eigenvectorPairs =+ Key.sort (negate . Complex.magnitude . fst) $+ (\(vr,d,vlAdj) ->+ zip (Vector.toList d) $+ zip (map unwrapEigenvector $ Matrix.toRows vlAdj)+ (map unwrapEigenvector $ Matrix.toColumns vr)) $+ Square.eigensystem $ Square.mapSize ExtShape.IntIndexed transition++dominantEigenmatrixPairs :: (TallMatrix Double, TallMatrix Double)+dominantEigenmatrixPairs =+ (tallMatrixFromColumns $ map (Array.map Complex.realPart . fst . snd) $+ take 2 eigenvectorPairs,+ tallMatrixFromColumns $ map (Array.map Complex.realPart . snd . snd) $+ take 2 eigenvectorPairs)++rankReduced :: Matrix (Complex Double)+rankReduced =+ foldl1 Matrix.add $+ map+ (\(l,(vl,vr)) ->+ Matrix.scale (l / Vector.inner vl vr) $ outerProduct vr vl)+ (take 2 eigenvectorPairs)++powerParameters :: Vector Double -> [(Complex Double, Complex Double)]+powerParameters initial =+ let initialC = Vector.fromReal initial+ in map+ (\(l,(vl,vr)) ->+ (l, (vr!(P1,P1)) * Vector.inner vl initialC / Vector.inner vl vr))+ eigenvectorPairs+++probabilities :: Vector Double -> [Double]+probabilities = map (! (P1,P1)) . iterate (transition #*|)++probabilitiesApprox1 :: [Double]+probabilitiesApprox1 =+ map (1-) $ iterate (0.9083578427694454*) 0.9237265161044081++probabilitiesApprox2 :: [Double]+probabilitiesApprox2 =+ map (1-) $ iterate (0.9083578427694454*) 1.025929108036849++formatProbabilities :: [String]+formatProbabilities =+ List.zipWith5+ (printf "%3d & %1.5f & %1.5f & %1.5f & %1.5f \\\\")+ [0::Int ..]+ (probabilities initial1)+ probabilitiesApprox1+ (probabilities initial2)+ probabilitiesApprox2+++main :: IO ()+main = mapM_ putStrLn $ take 101 formatProbabilities
+ wuerfelschlange.cabal view
@@ -0,0 +1,60 @@+Name: wuerfelschlange+Version: 0.1+License: BSD3+License-File: LICENSE+Author: Henning Thielemann <haskell@henning-thielemann.de>+Maintainer: Henning Thielemann <haskell@henning-thielemann.de>+Homepage: http://code.henning-thielemann.de/wuerfelschlange/+Category: Math+Synopsis: Code for the dice chain problem+Description:+ Code for the dice chain problem+ <https://de.wikipedia.org/wiki/Würfelschlange>.+ It is similar to the Kruskal Count.+Tested-With: GHC==7.4.2, GHC==7.8.4, GHC==8.6.5+Cabal-Version: 1.14+Build-Type: Simple++Flag buildGraph+ description: Build graph generator+ default: True++Flag buildMatrix+ description: Build matrix computations+ default: True++Source-Repository this+ Tag: 0.1+ Type: darcs+ Location: http://code.henning-thielemann.de/wuerfelschlange/++Source-Repository head+ Type: darcs+ Location: http://code.henning-thielemann.de/wuerfelschlange/++Executable wuerfelschlange-graph+ If flag(buildGraph)+ Build-Depends:+ containers >=0.4 && <0.7,+ base >=4.5 && <5+ Else+ Buildable: False++ Main-Is: src/Graph.hs+ GHC-Options: -Wall+ Default-Language: Haskell98++Executable wuerfelschlange-matrix+ If flag(buildMatrix)+ Build-Depends:+ lapack >=0.4 && <0.6,+ netlib-ffi >=0.1.1 && <0.2,+ comfort-array >=0.5 && <0.6,+ utility-ht >=0.0.13 && <0.1,+ base >=4.5 && <5+ Else+ Buildable: False++ Main-Is: src/Matrix.hs+ GHC-Options: -Wall+ Default-Language: Haskell98