module SelfRepair where
import Random (StdGen, random, mkStdGen, randomR, randomRs, next)
import GHC.IO.Handle.FD (stdin)
import GHC.IO.Handle (hWaitForInput)
-- This module requires version 4.2 or better of the GHC 'base' library.
-- (i.e. use GHC >= 6.12)
{- A DEFECT AND FAULT TOLERANT MOLECULAR SYSTEM SIMULATION
by Max OrHai, 7 June 2010
for Christof Teuscher's ECE 510 at Portland State University
Problem Statement:
A 12 by 12 grid is occupied by 50 particles, each configurable as a horizontal or vertical wire, or a logic gate. These particles are affected by random motion, but not rotation. The grid has a set number of randomly distributed 'defects', which do not move or have any other function. At each discrete motion step, and with a given probability, a randomly selected particle may experience a 'fault', resulting in its removal from the grid. The particles are aware of their four immediate N, E, S, and W neighbors, and may bind together. There may be other forces or gradients introduced into the simulation, but the system must use only local information and interactions to maintain connections between inputs on the N and S walls of the grid, a logic gate particle, and outputs along the E grid wall.
Analysis:
It is clearly desirable to have the shortest circuit path possible, as this minimizes the probability of the circuit being affected by faults. Because any circuit path must bridge the N and S walls, the shortest possible circuit path is directly along the E wall, unless blocked by a defect. If a redundant path is desired, the second shortest possible is directly to the W of the shortest, the third directly W of that, etc. Therefore, the particles should accumulate along the E wall.
My Solution:
First, I define symbolic data to describe the particles and the states they may be in. My Particles don't ever bind to each other, so they have no binding sites. For convenience, defects are implemented as inert, immobile Particles. -}
data Flow x = Still | NFlowing x | SFlowing x | EFlowing x deriving Eq
data Particle = Defect
| HWire (Flow Bool)
| VWire (Flow Bool)
| Gate (Flow Bool) deriving Eq
instance Show Particle where
show Defect = "#"
show (HWire Still) = "-"
show (HWire (EFlowing _)) = ">"
show (VWire Still) = "|"
show (VWire (NFlowing _)) = "^"
show (VWire (SFlowing _)) = "v"
show (Gate Still) = "G"
show (Gate (EFlowing _)) = "@"
show _ = "?" -- to catch errors
type Location = (Int,Int)
type Circuit = (Bool,Bool,[(Location, Particle)])
-- inputs N , S , a list of pairs
-- It's important to note that the grid itself, and the blank spaces it contains, are not explicitly represented in the Circuit. The list of location-particle pairs may be in any order; in fact, randomly reordering it is part of the process. Only when the Circuit is printed out by the visualization functions does it assume an explicit grid form.
-- Global parameters:
(xMax,yMax) = (12,12) :: (Int,Int)
initialParticleCount = 50
defectRate = 0.01 :: Float -- once upon initialization
faultRate = 0.001 :: Float -- once per step
-- The number of Defects in a grid is proportional to its area:
defectCount = round (defectRate * (fromIntegral xMax) * (fromIntegral yMax))
-- Shuffle a list by selecting and popping a random element from the list to the head, 3 times for each item in the list. I'm not convinced that this is the most effective way to shuffle, but it seems to work.
pop n xs = ((xs!!n):(take n xs)) ++ (drop (n + 1) xs)
randPop xs g = pop (fst $ randomR (0, length xs -1) g) xs
shuffle :: [a] -> Int -> [a]
shuffle xs seed = shuffle' xs (3*(length xs)) (mkStdGen seed)
shuffle' xs n gen | n <= 0 = xs
| otherwise = let (rand, newgen) = randomR (1, length xs - 1) gen
in shuffle' (pop rand xs) (n-1) newgen
randLocs :: Int -> [Location]
randLocs seed = shuffle [(x,y) | x <- [1..xMax], y <- [1..yMax]] seed
-- Make a new Circuit with randomly scattered particles:
randCircuit seed = (True,False,
(zip (randLocs seed)
((replicate defectCount Defect)
++ (replicate initialParticleCount (HWire Still)))))
-- Make a new Circuit with all particles already settled along the E wall:
-- ( saves about 1200 steps at startup )
demoCircuit seed = (True, False, defects ++ xyps)
where defects = zip (randLocs seed) (replicate defectCount Defect)
xyps = take initialParticleCount
[((x,y),VWire Still)
| x <- [xMax,xMax-1..1], y <- [1..yMax],
(x,y) `notElem` (fst $ unzip defects)]
------------------- Visualization and User Interface: --------------------
prc = putStrLn . cString -- Print an ASCII art diagram of a given Circuit.
-- Read an output signal from this circuit, if one is present:
output :: Circuit -> Char
output (_,_,xyps) = ch [p | p <- xyps, outputting p]
where outputting ((x,_),(HWire (EFlowing _))) | x == xMax = True
outputting ((x,_),(Gate (EFlowing _))) | x == xMax = True
outputting _ = False
ch [] = ':' -- no inputs
ch [((xMax,_),HWire (EFlowing b))] = boolToChar b
ch [((xMax,_),Gate (EFlowing b))] = boolToChar b
ch _ = '?' -- more than one input; a very transient state.
boolToChar b = if b then '1' else '0' -- Numeric Boolean values,
b0 = False -- for all you EE folks
b1 = True -- and C programmers (-:
-- Construct a grid representation as a list of characters:
cString :: Circuit -> String
cString c@(inN, inS, xyPs) =
(replicate xMax (boolToChar inN))
++ ('\n':(lineBrk $ map (cell) [(i,j)|j<-[yMax,yMax-1..1],i<-[1..xMax]]))
++ (replicate xMax (boolToChar inS))
where cell (i,j) = unMaybe $ lookup (i,j) xyPs
unMaybe Nothing = ' '
unMaybe (Just p) = head $ show p
lineBrk [] = []
lineBrk chars = (take xMax chars)
++ (' ':(output c):'\n':(lineBrk $ drop xMax chars))
-- Cheap animation by repeated printing and stepping, with a brief delay.
-- Press to pause, to exit, 'n' and 's' to toggle inputs.
run :: Circuit -> Int -> IO (Circuit)
run c seed = run' (step c seed) 0
where
run' (circ,gen) numSteps =
let (circ'@(inN,inS,xyps),gen') = step' (circ,gen)
continue c = run' (c,gen') (numSteps+1)
keyAction k | k == ' ' = waitForSpace
| k == '\n' = return circ'
| k == 'n' = continue (not inN,inS,xyps)
| k == 's' = continue (inN,not inS,xyps)
| otherwise = continue circ' -- ignore other keys
waitForSpace = do key <- getChar
if key == ' ' then continue circ'
else waitForSpace
in do putStrLn $ show (length xyps) ++ " Particles at step "
++ show numSteps ++ ":"
prc circ'; putStrLn "" ;
keyDown <- hWaitForInput stdin 50 -- pause (in ms) for animation timing
if keyDown then do { key <- getChar ; keyAction key }
else continue circ'
rrun seed = run (randCircuit seed) seed
drun seed = run (demoCircuit seed) seed
-- How long can this Circuit run before encountering an output gap of this length?
stress :: Circuit -> Int -> Int -> IO ()
stress c maxGap seed = stress' (step c seed) 0 0
where stress' (c@(_,_,xyps),g) numSteps thisGap
| thisGap > maxGap = putStr "\n"
| otherwise = do putStr ("Step " ++ (show numSteps) ++ ": "
++ (show $ length xyps)
++ " Particles, Gap length "
++ (show thisGap) ++ ".\r")
stress' (step' (c,g)) (numSteps + 1)
(if (output c) == ':' then (thisGap + 1)
else 0)
{- ------------- Movement and Interactions -----------------------
For simplicity, and to remain clearly within the assignment guidlines, I've kept (stochastic) movement and faulting just as specified in the assignment description, functionally separate from other (deterministic) interaction rules in the system. So, this 'step' is actually two steps, the second of which is made of lots of little sub-steps.
1. 'shake' picks a random particle and either deletes it with probability given by 'faultRate', or moves it in a biased random direction obtained by sampling 'dxdyDist', if allowed by the wall and neighborhood occupancy constraints. If the Particle is not free to move, nothing happens.
2. 'update' applies the first matching transformation rule to each Particle in a Circuit. It is essentially an update function for a 2D von Neumann neighborhood cellular automaton on a randomly mobile partial-grid substrate. In practice, on a sequential computer, this takes as many sub-steps as there are Particles in the Circuit.
As a consequence, the random motion induced by 'shake' has a rate considerably slower than the "parallel" state changes propagated by 'update'. So, it sometimes takes a while for new particles to find the holes left by faults, even when a particle is directly W of a hole. This algorithm would repair faults a little faster (and perhaps be more physically realistic as a simulation) if the stochastic motion step were also applied "in parallel". One way to do this would be to 'shake' every individual particle in the list once per step. In practice, this would be done in sequence, probably by passing an intermediate list; the implicit sequencing would take care of collisions. -}
between lo hi n = n >= lo && n <= hi -- Interval-checking, doubly inclusive.
-- The movement distribution is biased toward the E to simulate a physical force of some kind (e.g. magnetic or electrostatic) which tends to pull the particles toward the E wall of the circuit.
dxdyDist r | between 0.0 0.02 r = (-1, 0 ) -- move W : 2% chance
| between 0.02 0.5 r = ( 1, 0 ) -- E : 48%
| between 0.05 0.75 r = ( 0, 1 ) -- N : 25%
| between 0.75 1.0 r = ( 0,-1 ) -- S : 25%
step c seed = step' (c, mkStdGen seed)
steps n c seed = (iterate step' (step c seed)) !! n
psteps n c = prc $ fst $ steps n c 777
step' :: (Circuit, StdGen) -> (Circuit, StdGen)
step' (c@(iN,iS,xyps),g) = (update (iN,iS, shake rxyps), g')
where (randFrac, g') = random g :: (Float, StdGen)
-- randFrac is a Float in the interval [0,1)
rxyps = randPop xyps g -- pop a random element of the list
-- shake x = x -- uncomment this line to disable 'shake'
shake ((xy,Defect):rest) = xyps -- Defects don't move; try again
shake ps = if randFrac < faultRate then tail ps
else safeMove ps (dxdyDist randFrac)
-- The head of the list is the particle to move, if it...
safeMove xyps@(((x,y),p):rest) (dx,dy) =
let xy'@(x',y') = (x+dx,y+dy)
in if xy' `notElem` (fst $ unzip rest) && -- ... doesn't collide and
between 1 xMax x' && between 1 yMax y' -- stays within the bounds
then ((xy',p):rest) else xyps
-- Apply the first matching 't' rule to every Particle in the Circuit
update :: Circuit -> Circuit
-- update x = x -- uncomment this line to disable 'update'
update (iN,iS,ps) = (iN,iS, map t ps) where
-- Neighbors:
en (x,y) = if x == xMax then Just (VWire Still) else lookup (x+1,y) ps
-- the E wall is made of these ^^^^^^^^^^^
nn (x,y) = if y == yMax then Just (VWire (SFlowing iN)) else lookup (x,y+1) ps
-- the N wall is made of these ^^^^^^^^^^^^^^^^^^
wn (x,y) = lookup (x-1,y) ps -- there is no W wall needed
sn (x,y) = if y == 1 then Just (VWire (NFlowing iS)) else lookup (x,y-1) ps
-- the S wall is made of ^^^^^^^^^^^^^^^^^^
t xyp@(_,Defect) = xyp -- Defects don't transform
t (xy,p) = xyp' -- ... meaning the step's particle state at this (x,y)
(xy,p) (nn xy) (en xy) (sn xy) (wn xy)
-- self ^ ^ ^ ^ ... four neighbors: five args total
-- When N and S flowing signals collide, they form two (NAND) Gates:
xyp' (xy, _) (Just (VWire (SFlowing flowN))) _
(Just (VWire (NFlowing flowS))) _
= (xy,(Gate (EFlowing (not (flowN && flowS)))))
-- ... and on the next step the S Gate of the two becomes a VWire:
xyp' (xy,Gate _) (Just (Gate _)) _ (Just (VWire (NFlowing f))) _
= (xy,(VWire (NFlowing f)))
-- Incoming flow induces realignment, and E-ward flow has priority:
xyp' (xy,p) _ _ _ (Just (Gate (EFlowing f))) = (xy,HWire (EFlowing f))
xyp' (xy,p) _ _ _ (Just (HWire (EFlowing f))) = (xy,HWire (EFlowing f))
xyp' (xy,VWire _) (Just (VWire (SFlowing f))) _ _ _ = (xy,VWire (SFlowing f))
xyp' (xy,VWire _) _ _ (Just (VWire (NFlowing f))) _ = (xy,VWire (NFlowing f))
-- Lay flat along any obstacles to the E:
xyp' (xy,HWire Still) _ e _ _ | e /= Nothing = (xy,VWire Still)
-- If the incoming flow ceases, become still
xyp' (xy,Gate (EFlowing _)) _ _ _ _ = (xy,VWire Still)
xyp' (xy,HWire (EFlowing _)) _ _ _ _ = (xy,HWire Still)
xyp' (xy,VWire (NFlowing _)) _ _ _ _ = (xy,VWire Still)
xyp' (xy,VWire (SFlowing _)) _ _ _ _ = (xy,VWire Still)
xyp' xyp _ _ _ _ = xyp -- Default catch-all rule: no change.
{- -------------- Conclusions and potential improvements: -----------------
In summary, I believe I have found a fairly simple algorithm which satisfies the assignment criteria and, further, is able to assemble a circuit from random initial conditions. My model uses no particle binding or gradients, only a universally applied biased stochastic motion and a deterministic cellular automaton update rule. The model is able to make a circuit in about 250 random single-particle motion steps, and attains a theoretically optimum configuration in about 1200 steps, starting from uniformly distributed random initial conditions. It will tend to repair itself while running, as long as there are enough remaining particles to complete a circuit.
The 'stress' function allows for easy evaluation of circuit robustness under various parameters. With a fault rate of 0.01 (ten times the assignment criteria's assumed rate), and beginning from pre-settled initial conditions (which don't have a defect in the farthest-E columns, as in the "lucky" random seed '7' in the example below), the model is able to maintain connectivity for 1990 steps with no output gaps over 50 time steps. If we allow gaps of up to 100 steps, the circuit will run for 4733 steps with the random seed '7'.
*SelfRepair> stress (demoCircuit 7) 50 7
Step 1990: 34 Particles, Gap length 50.
*SelfRepair> stress (demoCircuit 7) 100 7
Step 4733: 18 Particles, Gap length 100.
There remains some room for improvement of this approach. The first change I would make would be to synchronize the rates of particle movement and CA updates, as detailed above in the "Movement" section. Moreover, the simple 'update' rules don't handle defects in the E-most two or three columns very effectively, as they are unable to rout around defects in any way besides filling in another column to the W of them, which is somewhat wasteful. I have some preliminary ideas for an extension of this model, using another fundamental Particle type which would rout signal Flows around right angle bends, but I have not yet implemented them.
My primary motivation toward imagining such an extension is my skepticism toward one of the assumptions of the assignment: that faulty particles simply disappear from the grid, rather than clog it with their useless "corpses". As it turns out, the built-in capacity to die neatly, shrivel up, and disappear is actually an important and near-universal "feature" of both uni- and multi-cellular life, called Programmed Cell Death (PCD). Identified forms of PCD include apoptosis, autophagy, and anoikis; however, the mechanisms employed in these processes are far from simple or well-understood, and are unlikely to be part of any (non-biological) physical implementation of this circuit that I can readily imagine.
A robust self-repairing circuit similar to this one but able to accomodate a more realistic particle failure mode would probably require substantially more sophisticated signal routing rules, and possibly three spatial dimensions rather than just two. However, I think that the present model could accomodate these extensions gracefully, although a more sophisticated graphical display function would be helpful for such a project.
Finally, I believe that the the use of a language like Haskell, which has built-in pattern matching and higher-order functions, has been instrumental to the success of this approach. Haskell's pattern-matching function definition semantics made it easy to write what amounts to a simple domain-specific language for cellular automaton rules, without sacrificing performance or increasing code complexity. Thanks, Haskell! -}