> {-# LANGUAGE ViewPatterns #-} > module Main > where Immutable unboxed arrays will be used for storing our 2D reference coloring, Vectors for shuffling our row and column permutations, and mersenne-random for faster randomness: > import Data.Array > import qualified Data.Vector.Unboxed as V > import System.Random.Mersenne > import System.Environment > import Text.CSV > import Control.Monad.State > import Control.Arrow(first,second) > import Data.List > import Graphics.HGL TODO: -check out statistics, adjust threshold schedule!!!! (like picking a lock) -check results of running with initial threshold of zero = ladder climbing -consider adding some kind of statistical element to the threshold. I think it is not fine-grained enough for us: i.e. there is most progress made within a very small span of a single threshold -try using a threshold based on making as few cells as possible have overlaps (i.e. so one cell with 4 overlaps is better than 2 cells with 2 overlaps) this should improve our step below) REALLY? WHY?: -try rendering the current 22-color grid, then remove some redundant cells and add cells that are opened later -after obtaining a population of 73-74-colorings, use this script as a fitness test for an evolutionary approach to combining various 74- colorings -continue with search script, to use a better than random cell picking/ adjusting scheme. -------------- | DATA TYPES | -------------- The reference grid is the original optimal single-coloring organized as an Array of (x,y) coordinates to Bool values indicating colored-ness. We represent permutations of this grid as simply permutations of row and column numbers: > type ReferenceGrid = Array (Int,Int) Bool The state of the grid is made up of a set of four cell-colorings, each a different permutation of our reference grid. These permutations are represented as pairs of (row numbers sequence, column number sequence) of the reference coloring: > data Grid = Grid { > redCells :: SubGrid > , grnCells :: SubGrid > , bluCells :: SubGrid > , ylwCells :: SubGrid > } deriving Show Automorphisms of our reference graph are created by shuffling rows and columns, so we will store these permutations as sequences of row and column numbers: > type SubGrid = (Order, Order) -- (X order, Y order) > type Order = V.Vector Int Our state monad encapsulation type: > type Annealer a = State Comp a > > data Comp = Comp { > grid :: Grid, > -- store rendering mechanism as a section with our > -- reference grid already consumed: > renderCol :: Int -> Grid -> AllSlices, > renderRow :: Int -> Grid -> AllSlices, > -- sources of randomness: > randDims :: [Int], > randCols :: [Int], > randRows :: [Int], > -- each round will increment from 0 until > head schedule > timeout :: !Int, > schedule :: ![Int], > -- bookkeeping: > globalScore :: !Int, > thresholds :: ![Int] > } > threshold, roundLimit :: Comp -> Int > threshold = head . thresholds > roundLimit = head . schedule We score current and potential solutions based on the number of UNCOLORED cells in a row column based on the union of all four colors. We compare these scores before and after a potential swap, and make the swap if the difference is not above a given threshold: TLDR: this is the number of uncolored cells in this row/column, low score is good like golf. Score of zero is optimal: > type Score = Int ----------------------- | META-HEURISTIC CODE | ----------------------- We use a technique similar to "simulated annealing" called "threshold accepting", in which swap two rows or two columns randomly, accepting the change if the swap would result in a better solution or a solution that is no worse than some threshold. This threshold changes as a function of time, eventually reaching zero, meaning that only beneficial mutations are accepted. This finds the local minimum of our neighborhood. Here is the schedule for our rounds. I am guessing that we should make later rounds longer, since we want to fully explore the local search space... ? > mkRoundSchedule :: Int -> [Int] > --mkRoundSchedule i = [ i, i*2.. ] > mkRoundSchedule i = repeat i The initial threshold is apparently important. Here the threshold means "ANY SWAP WHERE THE NEW SCORE MINUS THE OLD SCORE IS GREATER THAN $THRESHOLD WILL NOT BE ACCEPTED" For now we make a wild guess at a good threshold for 17x17: > --thresholdSchedule :: [Int] > --thresholdSchedule = [4,3.. -1] > mkThresholdSchedule :: Int -> [Int] > mkThresholdSchedule i = [i, i-1 .. -1] And we guess at some good settings to get a decent initial shuffle: > shuffleIterations, shuffleThresh :: Int > shuffleIterations = 9999 > shuffleThresh = 99 Here is the code we use to initialize and monitor the calculation, and adjust the threshold as the calculation proceeds. It starts by running for a long period with a very high threshold, simply to get an initial shuffle, then proceeds with the annealing procedure. The initial shuffle would be redundant given a proper annealing schedule, but I'm not confident I know what that would be: > annealer :: Annealer [Score] > annealer = do > s <- get > let (lim:lims) = schedule s > to = timeout s > (_:ts) = thresholds s > -- if this round is over, > if to > lim > -- and if the round just completed was round zero > then if null ts > then return [] > -- if not last: reset clock, move on to next round > else do put s{schedule = lims, thresholds = ts, timeout = 0} > annealer > -- otherwise just do an attempted swap in this round > else do gs <- swap > s' <- get > put s'{ timeout = to+1 } > -- return a running list of global scores: > gss <- annealer > -- if we found a permutation that is a total cover: SUCCESS! > if gs <= 0 > then return [gs] > else return (gs:gss) > -------------------------- | PROGRAM INITIALIZATION | -------------------------- > main = do > > -- read the CSV file of our initial grid > [f,lengthFstRound,initialThreshold] <- getArgs > (Right rs) <- parseCSVFromFile f > -- initial global score is same as number of free cells in our input > -- single coloring: > let (scoreI,refGrid) = buildRefGrid rs > (xN,yN) = snd (bounds refGrid) > -- initially all four colors overlap in all cells: > ordI = (V.fromList [0..xN], V.fromList [0..yN]) > gridI = Grid ordI ordI ordI ordI > roundSchedule = mkRoundSchedule (read lengthFstRound) > thresholdSchedule = mkThresholdSchedule $ read initialThreshold > > -- our mersenne random streams, for choosing row/col swaps: > rs1 <- getStdGen >>= randoms > rs2 <- getStdGen >>= randoms > rs3 <- getStdGen >>= randoms > let rrows = map (`mod` (yN+1)) rs1 ::[Int] > rcols = map (`mod` (xN+1)) rs2 ::[Int] > rdims = map (`mod` 8) rs3 ::[Int] > > > -- -- Build an initial state for the shuffle phase -- -- > > putStrLn "Running initial shuffle phase..." > let stI = Comp { grid = gridI, > renderCol = renderColWith refGrid , > renderRow = renderRowWith refGrid, > randDims = rdims, randCols = rcols, randRows = rrows, > -- for tracking & controlling computation: > globalScore = scoreI, timeout = 0, > -- annealing meta-parameters: > schedule = [shuffleIterations], > thresholds = [shuffleThresh] } > > let stShuffled = execState annealer stI > > > -- -- Run a set of random permutations to get of avg change -- -- > > putStrLn "Getting the average change for random free swaps..." > let stAvg = stShuffled { > timeout = 0, > schedule = [shuffleIterations], > thresholds = [shuffleThresh] } > > let ((scr0:scores),stAveraged) = runState annealer stAvg > delt prev cur = (cur, abs$ prev - cur) > delts = snd $ mapAccumL delt scr0 scores > deltGroups = map (\g-> (head g,length g))$ group $ sort delts > avgDelta = fromIntegral (sum delts) / fromIntegral shuffleIterations > > -- we should do something with this to get a good annealing schedule, but > -- not sure how at the moment: > print deltGroups > print avgDelta > --error "STOP" > > -- -- Run actual annealing procedure -- -- > > let stAnn = stAveraged { > timeout = 0, > schedule = roundSchedule, > thresholds = thresholdSchedule } > > let (scoresAnn,stFinal) = runState annealer stAnn > > -- print the running list of global scores, then show the final > -- state of the grid (i.e. our solution): > mapM_ putStrLn $ scoresToLines scoresAnn > putStrLn "------- OUR FINAL PERMUTATIONS: --------" > print $ grid stFinal > > putStrLn "Rendering four colors:" > runGraphics$ prettyPicturesFrom (grid stFinal) refGrid Takes a list of scores and produces a list of output line strings > scoresToLines :: [Int] -> [String] > scoresToLines = map (\(x,y)-> (show x)++" "++(show y)) . zip [0..] > --scoresToLines = map show . zip [0..] ----------------------- | ANNEALING PROCEDURE | ----------------------- HOW WE USE RANDOMNESS: ====================== We need 3 random Int streams to choose random row or column swaps: the first will be taken modulo 8 and will correspond to: [ col(Red), row(Red), col(Grn) ... row(Ylw) ] The second and third streams will correspond to column and row numberings respectively (we could use one stream here if we knew we only cared about square grids). Here is an example step We find the number 2 at the head of our first stream, so we will be swapping two columns in the Green subset. We see that the second stream looks like: [12,3,7...] so we will swap Green columns '12' and '3' if the resulting change is acceptable according to our current temperature threshold/cost function, otherwise the two columns will remain where they were. OUR ANNEALING/THRESHOLD-ACCEPTANCE CODE: ======================================== Maybe swap two rows or columns randomly, returning the new global score after the swap: > swap :: Annealer Score > swap = do > s <- get > let (d:ds) = randDims s > (c1:c2:cs) = dnub $ randCols s -- lazy bindings > (r1:r2:rs) = dnub $ randRows s > if even d then put s{ randCols = cs, randDims = ds } > else put s{ randRows = rs, randDims = ds } > case d of > 0 -> do swapRed True c1 c2 > 1 -> do swapRed False r1 r2 > 2 -> do swapGreen True c1 c2 > 3 -> do swapGreen False r1 r2 > 4 -> do swapBlue True c1 c2 > 5 -> do swapBlue False r1 r2 > 6 -> do swapYellow True c1 c2 > 7 -> do swapYellow False r1 r2 > _ -> error "sad panda..." > gets globalScore > UGH sorry about the boilerplate > swapRed :: Bool -> Int -> Int -> Annealer () > swapRed doCols c1 c2 = do > s <- get > let g = grid s > vecs = redCells g > (vec,onColOrRow,render) = if doCols > then (fst vecs,first,renderCol s) > else (snd vecs,second,renderRow s) > -- we render the two rows as arrays of bools of each color: > -- (use lazy state monad to ensure these bindings are lazy) > (c1sliceR,c1sliceG,c1sliceB,c1sliceY) = render c1 g > (c2sliceR,c2sliceG,c2sliceB,c2sliceY) = render c2 g > -- Here we score the current local solution against the > -- proposed solution after swapping: > let gby12 = ((c1sliceG,c1sliceB,c1sliceY), > (c2sliceG,c2sliceB,c2sliceY)) > scrs = score (c1sliceR,c2sliceR) gby12 > -- if the swap was accepted, put the new vector back in: > mvec' <- maybeSwap (c1,c2,vec) scrs > s' <- get > let putFunc v = put s'{grid = g{redCells = onColOrRow (const v) vecs}} > maybe (return ()) putFunc mvec' > > swapGreen doCols c1 c2 = do > s <- get > let g = grid s > vecs = grnCells g > (vec,onColOrRow,render) = if doCols > then (fst vecs,first,renderCol s) > else (snd vecs,second,renderRow s) > (c1sliceR,c1sliceG,c1sliceB,c1sliceY) = render c1 g > (c2sliceR,c2sliceG,c2sliceB,c2sliceY) = render c2 g > rby12 = ((c1sliceR,c1sliceB,c1sliceY), > (c2sliceR,c2sliceB,c2sliceY)) > scrs = score (c1sliceG,c2sliceG) rby12 > mvec' <- maybeSwap (c1,c2,vec) scrs > s' <- get > let putFunc v = put s'{grid = g{ grnCells = onColOrRow (const v) vecs}} > maybe (return ()) putFunc mvec' > > swapBlue doCols c1 c2 = do > s <- get > let g = grid s > vecs = bluCells g > (vec,onColOrRow,render) = if doCols > then (fst vecs,first,renderCol s) > else (snd vecs,second,renderRow s) > (c1sliceR,c1sliceG,c1sliceB,c1sliceY) = render c1 g > (c2sliceR,c2sliceG,c2sliceB,c2sliceY) = render c2 g > rgy12 = ((c1sliceR,c1sliceG,c1sliceY), > (c2sliceR,c2sliceG,c2sliceY)) > scrs = score (c1sliceB,c2sliceB) rgy12 > mvec' <- maybeSwap (c1,c2,vec) scrs > s' <- get > let putFunc v = put s'{grid = g{ bluCells =onColOrRow (const v) vecs }} > maybe (return ()) putFunc mvec' > > swapYellow doCols c1 c2 = do > s <- get > let g = grid s > vecs = ylwCells g > (vec,onColOrRow,render) = if doCols > then (fst vecs,first,renderCol s) > else (snd vecs,second,renderRow s) > (c1sliceR,c1sliceG,c1sliceB,c1sliceY) = render c1 g > (c2sliceR,c2sliceG,c2sliceB,c2sliceY) = render c2 g > rgb12 = ((c1sliceR,c1sliceG,c1sliceB), > (c2sliceR,c2sliceG,c2sliceB)) > scrs = score (c1sliceY,c2sliceY) rgb12 > > mvec' <- maybeSwap (c1,c2,vec) scrs > s' <- get > let putFunc v = put s'{grid = g{ ylwCells = onColOrRow (const v) vecs }} > maybe (return ()) putFunc mvec' > Here we swap the two vector elements (corresponding to row/col numbers) if swapping them leads to a better solution or a worse solution that is not worse than our allowed threshold (see below for more on that). > maybeSwap :: (Int,Int,Order) -> (Score,Score) -> Annealer (Maybe Order) > maybeSwap (s1,s2,vec) (curScore,newScore) = do > s <- get > let t = threshold s > sc = globalScore s --lower score is better > t' = newScore - curScore --negative val here means improvement > swappedVec = vec V.// [(s1,vec V.! s2 ),(s2,vec V.! s1 )] > if t' > t then return Nothing > -- or update the global score and return the new vector: > else do put s{ globalScore = sc + t' } > return$ Just swappedVec SCORING: ======== It would be really sweet to use DPH and use parallel arrays here (and in our reference array) to do lookups, unions in parallel. Too bad I have barely one core :-( > type Slice = V.Vector Bool > type AllSlices = (Slice,Slice,Slice,Slice) > type ThreeSlices = (Slice,Slice,Slice) Compare fst with union of fst, snd with union of snd, and sum to get the current score. Then compare if fst and snd swapped. Return: (current score, new potential score) > score :: (Slice,Slice) -> (ThreeSlices,ThreeSlices) -> (Score,Score) > score (v1,v2) (vUnion-> vu1, vUnion-> vu2) = > (sc v1 vu1 + sc v2 vu2, sc v2 vu1 + sc v1 vu2) > -- a value of 1 is given where the union of all colors is False: > where sc s = V.sum . V.zipWith (\v vu-> if v || vu then 0 else 1) s > vUnion :: ThreeSlices -> Slice > vUnion (v1,v2,v3) = V.zipWith3 (\x y z->x || y || z ) v1 v2 v3 ---------------- | SOME HELPERS | ---------------- GETTING ARRAY SLICES ===================== Return either row or column 'n' of all colors. Build a row/column renderer which we initialize at start and carry in our State. Hopefully this is an optimization: > renderColWith,renderRowWith :: ReferenceGrid -> > Int -> Grid -> AllSlices > renderColWith ref x (Grid rd gn bl yw) = > let render ((V.!x)->x', ys) = V.map (\y-> ref!(x',y)) ys > in (render rd, render gn, render bl, render yw) > > renderRowWith ref y (Grid rd gn bl yw) = > let render (xs, (V.!y)->y') = V.map (\x-> ref!(x,y')) xs > in (render rd, render gn, render bl, render yw) GRAPHICS: ========= > prettyPicturesFrom :: Grid -> ReferenceGrid -> IO () > prettyPicturesFrom (Grid r g b y) ref = do > let rendGraphic c = gridToGraphic c . renderAllColored ref > (wd,ht) = (V.length$ fst r, V.length$ snd r) > w <- openWindow "RGBY" (wd*20, ht*20) > directDraw w $ (rendGraphic Red r) `overGraphic` > (rendGraphic Green g)`overGraphic` > (rendGraphic Blue b)`overGraphic` > (rendGraphic Yellow y) > gridToGraphic c = overGraphics . map toGraphic > where toGraphic ((*20)->x, (*20)->y) = > withColor c $ polygon [(x,y),(x,y+18),(x+18,y+18),(x+18,y)] > renderAllColored :: ReferenceGrid -> SubGrid -> [(Int,Int)] > renderAllColored ref (xs,ys) = > let xCs = [0.. V.length xs - 1] > yCs = [0.. V.length ys - 1] > in [ (x,y) | x <- xCs, y <- yCs, ref!(xs V.! x, ys V.! y) ] OTHER: ====== Allows us to avoid self-swapping without using 'nub' which I believe would be bad statistical juju: > dnub :: [Int] -> [Int] > dnub xss@(x1:x2:xs) > | x1 == x2 = dnub (x2:xs) > | otherwise = xss Build the initial reference single coloring from our CSV file. > buildRefGrid :: CSV -> (Score,ReferenceGrid) > buildRefGrid cs = > let xN = subtract 1 $ maximum $ map length cs > yN = subtract 1 $ length cs > > cellOrdering = [(x,y) | y <-[0..yN], x <-[0..xN]] > cs' = concatMap (take (xN+1) . (++ repeat "" )) cs > cells = map (not . null) cs' > > refGrid = array ((0,0),(xN,yN)) $ zip cellOrdering cells > initialGlobalScore = length $ filter not cells > > in (initialGlobalScore, refGrid)