-
Notifications
You must be signed in to change notification settings - Fork 1
/
Main.hs
170 lines (151 loc) · 5.3 KB
/
Main.hs
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
{-# LANGUAGE RecordWildCards #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE TemplateHaskell #-}
{-# LANGUAGE ViewPatterns #-}
{-# LANGUAGE LambdaCase #-}
import Bio.Motions.Prototype as Prototype
import Control.Monad.State.Strict
import Control.Monad.Except
import Options.Applicative
import System.IO
import System.Random
import Control.Lens
import Control.Monad.Random
import qualified Data.Serialize as Serialize
import qualified Data.ByteString as BS
data Settings
= Initialize InitializeSettings
| Simulate SimulateSettings
data InitializeSettings = InitializeSettings
{ settingsLaminsFile :: FilePath
, settingsBindersFile :: FilePath
, settingsRadius :: Int
, settingsNumBinders :: Int
, settingsChainLength :: Int
, settingsInitOutputFile :: FilePath
}
data SimulateSettings = SimulateSettings
{ settingsPDBFile :: FilePath
, settingsInputFile :: FilePath
, settingsOutputFile :: FilePath
, settingsNumSteps :: Int
, settingsWriteIntermediateStates :: Bool
}
data Counter = Counter
{ _counterNumSteps :: !Int
, _counterNumFrames :: !Int
}
makeLenses ''Counter
initializeParser :: Parser InitializeSettings
initializeParser = InitializeSettings
<$> strArgument
(metavar "LAMIN-BSITES"
<> help "File containing the lamin binding sites")
<*> strArgument
(metavar "BINDER-BSITES"
<> help "File containing the regular binding sites")
<*> option auto
(long "radius"
<> short 'r'
<> metavar "RADIUS"
<> help "Radius of the bounding sphere"
<> value 10)
<*> option auto
(long "num-binders"
<> short 'b'
<> metavar "BINDERS"
<> help "Number of binders"
<> value 256)
<*> option auto
(long "chain-length"
<> short 'l'
<> help "Chain length"
<> value 512)
<*> strOption
(long "output"
<> short 'o'
<> metavar "OUTPUT-FILE"
<> help "Output file")
simulateParser :: Parser SimulateSettings
simulateParser = SimulateSettings
<$> strOption
(long "pdb"
<> short 'p'
<> metavar "PDB-FILE"
<> help "PDB output file")
<*> strOption
(long "input"
<> short 'i'
<> metavar "INPUT-FILE"
<> help "Input file")
<*> strOption
(long "output"
<> short 'o'
<> metavar "OUTPUT-FILE"
<> help "Output file")
<*> option auto
(long "steps"
<> short 's'
<> metavar "STEPS"
<> help "Number of simulaton steps"
<> value 100000)
<*> switch
(long "intermediate-states"
<> short 'm'
<> help "Whether to write the intermediate states to the output file")
parser :: Parser Settings
parser = subparser
$ command "init" (info (helper <*> (Initialize <$> initializeParser))
(progDesc "Initialize the simulation"))
<> command "run" (info (helper <*> (Simulate <$> simulateParser))
(progDesc "Run the simulation"))
makeInput :: InitializeSettings -> IO Input
makeInput InitializeSettings{..} = do
let inputChainLength = settingsChainLength
inputRadius = settingsRadius
inputNumBinders = settingsNumBinders
inputLamins <- map read . lines <$> readFile settingsLaminsFile
inputBinders <- map read . lines <$> readFile settingsBindersFile
return Input{..}
runAndWrite :: (MonadSimulation m, MonadIO m) => Maybe Handle -> StateT Counter m ()
runAndWrite handle = do
oldEnergy <- lift $ gets energy
lift simulateStep
newEnergy <- lift $ gets energy
when (oldEnergy /= newEnergy) $ pushPDB handle
counterNumSteps += 1
pushPDB :: (MonadSimulation m, MonadIO m) => Maybe Handle -> StateT Counter m ()
pushPDB Nothing = return ()
pushPDB (Just handle) = do
st@SimulationState{..} <- lift get
headerSequenceNumber <- use counterNumFrames
headerStep <- use counterNumSteps
let headerTitle = "chromosome;bonds=" ++ show energy
liftIO $ do
putStrLn $ "gyration radius: " ++ show gyrationRadius
putStrLn $ "energy: " ++ show energy
writePDB handle Header{..} st
hPutStrLn handle "END"
counterNumFrames += 1
serialize :: FilePath -> SimulationState -> StdGen -> IO ()
serialize file st gen = BS.writeFile file $ Serialize.encode (st, show gen)
run :: Settings -> IO ()
run (Initialize settings) = do
input <- makeInput settings
(est, gen) <- newStdGen >>= runRandT (runExceptT $ initialize input)
case est of
Left err -> hPutStrLn stderr err
Right st -> serialize (settingsInitOutputFile settings) st gen
run (Simulate SimulateSettings{..}) = withFile settingsPDBFile WriteMode $ \pdbFile ->
Serialize.decode <$> BS.readFile settingsInputFile >>= \case
Left err -> hPutStrLn stderr err
Right (st, read -> gen) -> do
(st', gen') <- flip runRandT gen $ flip execStateT st $ flip execStateT (Counter 0 0) $ do
pushPDB $ Just pdbFile
replicateM_ settingsNumSteps $ runAndWrite $
guard settingsWriteIntermediateStates >> Just pdbFile
serialize settingsOutputFile st' gen'
main :: IO ()
main = run =<< execParser (info (helper <*> parser)
(fullDesc
<> progDesc "Perform a MCMC simulation of chromatine movements"))