module SelfAssembly where
-- Self-assembly of my initials, a midterm project for Teuscher's ECE 510.
-- Max OrHai, Portland State University, May 2010
-- Executive summary: to interact with this program, load it into GHCi
-- and type "demo" or "try n" where n is an integer.
import Data.Array (accumArray, elems)
import Random (getStdRandom, getStdGen, mkStdGen, randomR, randomRs)
import GHC.IO.Handle.FD (stdin)
import GHC.IO.Handle (hWaitForInput)
demo = try 5 -- Complete in 12817 steps, but only for worldSize = (12,9)
type BindSites = (Int,Int,Int,Int) -- a tuple
-- ( N E S W )
type Location = (Int,Int)
-- ( x y )
type Particle = (Location,BindSites)
type Part = [Particle] -- a list
type World = [Part]
worldSize = (12,9) :: (Int,Int) -- (horizontal, vertical)
-- This seems to be about the right size for two letters.
-- ... bigger Worlds (like 40 x 80) are fun too.
xMax = snd worldSize - 1 -- takes care of zero indexing
yMax = fst worldSize - 1 -- and transposition (see below)
rotateParticleL :: Particle -> Particle
-- Rotate a Particle counter-clockwise
rotateParticleL (loc,(n,e,s,w)) = (loc,(e,s,w,n))
rotateParticleR :: Particle -> Particle
-- Rotate a Particle clockwise
rotateParticleR (loc,(n,e,s,w)) = (loc,(w,n,e,s))
moveParticle :: Location -> Particle -> Particle
-- move a Particle by an increment
moveParticle (dx,dy) ((x,y),b) = ((x + dx, y + dy),b)
-- Move all constituent Particles in a Part together.
-- Due to implicit currying, "move (dx,dy)" has type (Part -> Part)
move :: Location -> Part -> Part
move incr parts = map (moveParticle incr) parts
-- All the things that may be done to a Part in one time step:
opList = [movN, movE, movS, movW, rotL, rotR]
-- The directions are mapped to make sense with printWorld.
convertToScreenCoords (x,y) = (-y,x)
movN = move (-1, 0)
movE = move ( 0, 1)
movS = move ( 1, 0)
movW = move ( 0,-1)
-- Rotating bound parts properly:
-- Find the center of the structure, offset each part relative to it (localize), rotate each, convert back to global coordinates. If the length or width of the structure is even, then the structure will shift toward (0,0) by 1/2 cell, to fit the grid.
rotL = rotatePart True
rotR = rotatePart False
-- l is a Bool: True means counterclockwise, False means clockwise
rotatePart l pt = map (if l then rotateParticleL else rotateParticleR)
(move back ((if l then rotL' else rotR')
(move toCenter pt)))
where back = xyCenter pt
toCenter = (\(x,y) -> (-x,-y)) back
rotL' = map (\((x,y),b) -> ((-y,x),b))
rotR' = map (\((x,y),b) -> ((y,-x),b))
center ns = (least ns + greatest ns) `div` 2 -- let's stick with integers
xyCenter bps = (center xs, center ys) where (xs,ys) = unzip (fst (unzip bps))
most rel ns = foldl (\a b -> if a `rel` b then a else b) (head ns) (tail ns)
least = most (<)
greatest = most (>)
-- Does this Part fit into this World?
fitsInto :: Part -> World -> Bool
fitsInto part world = clear theseXYs otherXYs
where theseXYs = fst $ unzip part
otherXYs = concat $ map (fst . unzip) world
outOfBounds (x,y) = x < 0 || x > xMax || y < 0 || y > yMax
clear [] _ = True
clear (p:ps) w | p `elem` w = False -- collision
| outOfBounds p = False -- past world bounds
| otherwise = clear ps w
-- If possible in this World (which must NOT contain the given Part), apply this (movement or rotation) operation to this Part. Return only the Part, not the whole World.
safe :: (Part -> Part) -> Part -> World -> Part
safe op part w = let movedPart = op part
in if movedPart `fitsInto` w then movedPart else part
-- Can these two Particles be bound together? Keep in mind the screen coordinate system.
bindableP :: Particle -> Particle -> Bool
bindableP ((x1,y1),(n1,e1,s1,w1)) ((x2,y2),(n2,e2,s2,w2))
| dist == (-1, 0) && s2 `matches` n1 = True -- particle 2 is N of p1
| dist == ( 0, 1) && w2 `matches` e1 = True -- p2 E of p1
| dist == ( 1, 0) && n2 `matches` s1 = True -- p2 S of p1
| dist == ( 0,-1) && e2 `matches` w1 = True -- p2 W of p1
| otherwise = False
where dist = (x2 - x1, y2 - y1) -- directed distance from p1 to p2
matches a b = a /= 0 && a + b == 0 -- 0's don't bind
-- How about these two Parts? Check all possible combinations of Particles.
bindable :: Part -> Part -> Bool
bindable pt1 pt2 =
any good [(thisParticle,thatParticle) | thisParticle <- pt1, thatParticle <- pt2]
where good (p1,p2) = bindableP p1 p2
-- If possible, bind these Part to any others available in this World, and return a new World.
safeBindParts p w = safeBindParts' p [] w
-- Needs an auxilary World to hold the checked, unbindable Parts
safeBindParts' :: Part -> World -> World -> World
safeBindParts' part checked [] = (part:checked)
safeBindParts' thisPart checked (thatPart:unchecked) =
if bindable thisPart thatPart
then safeBindParts' (thisPart ++ thatPart) [] (checked ++ unchecked) -- start over again
else safeBindParts' thisPart (thatPart:checked) unchecked
-- The time step function: given an index to some Part within a World, a movement or rotation operation, and the World itself, apply the operation safely to the Part and bind if possible, returning the new World
step :: Int -> (Part -> Part) -> World -> World
step index op w = let selectedPart = w !! index
wMinus = (take index w) ++ (drop (index+1) w)
in safeBindParts (safe op selectedPart wMinus) wMinus
randElem :: [a] -> IO a
randElem xs = do randIndex <- getStdRandom (randomR (0, length xs - 1))
return (xs !! randIndex)
randStep :: World -> IO World
randStep w = do gen <- getStdGen -- irreproducible!
return (fst (randStep' w gen))
-- randStep' :: World -> g -> (World,g) -- reproducible.
randStep' world gen = (newWorld, newGen) where
(randPartIndex, newGen) = randomR (0, length world - 1) gen
(randOpIndex, _) = randomR (0, length opList -1) gen
newWorld = step randPartIndex (opList !! randOpIndex) world
randSteps n w = if n < 0 then return w
else do nw <- randStep w
printWorld nw
randSteps (n - 1) nw
--------- Visualization and user interface -----------------------
run = randStepUntilReducedTo -- an alias to save typing during experiments
-- Use your favorite number as a seed to try and assemble one M and one O:
try n = run 2 (wMO 1 1 n) n
randStepUntilReducedTo :: Int -> World -> Int -> IO ()
-- Supplying the same generator seed and initial conditions should always produce
-- the same results. This is the only significant imperative part of this program.
randStepUntilReducedTo n w genSeed =
do putStrLn ("Initial conditions: (" ++ show (length w) ++ " Parts)"); printWorld w
putStrLn "Attempting self-assembly. \
\ \nPress to see a snapshot, or any other key to cancel."
rSUR' w 0 (mkStdGen genSeed)
where
rSUR' w numSteps gen =
let stats = show (length w) ++ " Parts at step " ++ show numSteps
in if length w > n
then do let (w', newGen) = randStep' w gen
putStr ('\r':stats)
keyDown <- hWaitForInput stdin 0
if keyDown then do key <- getChar
if key == ' '
then do putStrLn ""; printWorld w'
rSUR' w' (numSteps + 1) newGen
else putStrLn ("\rCancelled: " ++ stats)
else rSUR' w' (numSteps + 1) newGen
else do putStrLn ""; printWorld w
putStrLn ("Reduced to " ++ show (length w) ++
" parts in " ++ show numSteps ++ " steps.")
-- Single character tags for Particles. The 'b' stands for BindSites.
debugging = False -- set to True for more meaningful (and confusing) display output.
b2char b = if debugging then db2char b else b2char' b
where
-- ...these are for debugging the movement, rotation, and binding functions:
db2char (1,0,0,0) = 'N'
db2char (0,1,0,0) = 'E'
db2char (0,0,1,0) = 'S'
db2char (0,0,0,1) = 'W'
db2char (-1,0,0,0) = 'n'
db2char (0,-1,0,0) = 'e'
db2char (0,0,-1,0) = 's'
db2char (0,0,0,-1) = 'w'
-- ...these are for debugging the actual letter Particles: uppercase is for O, lower for M.
db2char b | is oA = 'A' | is oB = 'B' | is oC = 'C' | is oD = 'D'
| is mA = 'a' | is mB = 'b' | is mC = 'c' | is mD = 'd' | is mE = 'e'
| is mF = 'f' | is mG = 'g' | is mH = 'h' | is mI = 'i' | is mJ = 'j'
| is mK = 'k' | otherwise = '?'
is = isRotationOf b
-- ... and these are for the normal display mode.
b2char' b | any is [oA,oB,oC,oD] = 'O'
| any is [mA,mB,mC,mD,mE,mF,mG,mH,mI,mJ,mK] = 'M'
| otherwise = '?' -- the catch-all clause for any unrecognized Particles.
isRotationOf a (b1,b2,b3,b4)
| (b1,b2,b3,b4) == a = True
| (b2,b3,b4,b1) == a = True
| (b3,b4,b1,b2) == a = True
| (b4,b1,b2,b3) == a = True
| otherwise = False
part2chars pt = map (\(loc,b) -> (loc,b2char b)) pt
accf _ x = x -- Needed for accumArray. Replacement is a kind of accumulation... (-;
maybeWrap = map (\bps -> map (\(xy,b) -> (xy,Just b)) bps) -- takes a World implicitly
worldArr w = accumArray accf Nothing ((0,0),(xMax,yMax)) (concat $ maybeWrap w)
mkString es = hzFrame ++ ('\n':(lineBrk $ map unMaybe es)) ++ hzFrame
where n = yMax + 1
hzFrame = ('+':replicate n '-' ++ "+")
unMaybe Nothing = ' '
unMaybe (Just b) = b2char b
lineBrk [] = []
lineBrk chars = ('|':(take n chars)) ++ "|\n" ++ (lineBrk $ drop n chars)
printWorld w = putStrLn $ mkString $ elems $ worldArr w
{- Note transposition of output, a side effect of printWorld because of the ordering
produced by 'elems' and the fact that things print top to bottom. I figured it
doesn't really matter since Parts can rotate anyway, and sort of hacked around it.
(x,y) --> (-y,x)
(1,0) moves S on printout
(0,1) moves E
(-1,0) moves N
(0,-1) moves W -}
------- testing ---------------
aa, bb :: Particle
aa = ((0,0),(1,0,0,0))
bb = ((2,3),(0,0,-1,0))
w0 = [[aa],[bb]] -- a very simple World for debugging purposes.
{------------- Design of the letters: ---------------
# # ### My initials, in a 5 by 5 pixel font.
## ## ## ## Each letter requires 16 particles.
##### # # Designing an algorithm to decompose
# # # ## ## them automatically would be an
# # ### interesting task. I did these by hand.
DAB The O is the easier one, so I did it first.
BC CD As you can see, it can be made with only four
A A each of four distinct kinds of particles.
DC CB Symmetry is nice.
BAD -}
oA = (1,0,-4,0)
oB = (0,2,-1,0)
oC = (3,0,0,-2)
oD = (0,4,-3,0)
bindSitesForOneO = repList 4 [oA,oB,oC,oD]
repList n xs = take (n * length xs) (cycle xs)
{-
A A The M appears to require 11 types of particle.
CD EF Most of the symmetry in this letter is
GHIJK reflective rather than rotational; I could
B A B have used fewer particle types if the
A A particles could flip as well as rotate. -}
mA = (10,0,0,0)
mB = (-10,0,11,0)
mC = (-10,12,14,0)
mD = (15,-12,0,0)
mE = (13,-19,0,0)
mF = (20,-13,-10,0)
mG = (-14,16,-11,0)
mH = (-16,-15,17,0)
mI = (18,-10,-17,0)
mJ = (-18,19,18,0)
mK = (-11,-21,-20,0)
bindSitesForOneM = [mA,mA,mA,mA,mA,mB,mB,mC,mD,mE,mF,mG,mH,mI,mJ,mK]
singletons xs = [[x] | x <- xs]
uniqueRandLocs genSeed = uRLs' [] (randomRs (0,(max xMax yMax)) (mkStdGen genSeed)) where
uRLs' locs (r1:r2:rs) | length locs == (xMax + 1) * (yMax + 1) = locs -- overflow!
| r1 > xMax || r2 > yMax = uRLs' locs (r2:rs) -- out of bounds
| (r1,r2) `notElem` locs = uRLs' ((r1,r2):locs) rs -- got one
| otherwise = uRLs' locs (r2:rs) -- keep trying
randomWorldWith bs gS = whirl (singletons (zip (uniqueRandLocs gS) bs)) gS
whirl wrld gS = zipWith (\p r -> (iterate rotL p) !! r) wrld (randomRs (0,4) (mkStdGen gS))
-- Takes a number of M's and O's and a generator seed;
-- returns a fresh scrambled World with the unbound Particles to generate the M's and O's.
wMO nMs nOs genSeed = randomWorldWith ((repList nMs bindSitesForOneM) ++
(repList nOs bindSitesForOneO)) genSeed