diff --git a/modelling/doc/monte_carlo.rst b/modelling/doc/monte_carlo.rst index 5e63a0a6674098e4244a5e20deae8cb70c23a4c3..9a0b6f5eef70a0bbe18b2ec804d6c78aed932b36 100644 --- a/modelling/doc/monte_carlo.rst +++ b/modelling/doc/monte_carlo.rst @@ -21,46 +21,7 @@ Carlo sampling to the N-terminal part of crambin: .. literalinclude:: ../../../tests/doc/scripts/modelling_monte_carlo.py -.. method:: SampleMonteCarlo(sampler, closer, scorer, cooler, steps,\ - bb_list, initialize=True, seed=0,\ - lowest_energy_conformation=True) - - A convenient function to perform Monte Carlo sampling using a simulated - annealing scheme. In every iteration, a new loop conformation gets proposed by - the provided *sampler* and closed by the *closer*. Upon scoring, this new - conformation gets accepted/rejected using a metropolis criterion based on the - temperature given by the *cooler* - => acceptance probability: exp(-delta_score/T). - The result is stored in *bb_list* and is either the lowest energy conformation - ever encountered or the last accepted proposal. - - :param sampler: Sampler object capable of initializing and altering - conformations. - :param closer: Closer object to adapt a new conformation to - the environment. - :param scorer: Scorer object to score new loop conformations. - :param cooler: Cooler object to control the temperature of the - Monte Carlo trajectory. - :param steps: Number of Monte Carlo iterations to be performed. - :param bb_list: The chosen conformation gets stored here. - :param initialize: Whether a new bb_list should be generated as starting - point, based on the samplers Initialize function. - The input *bb_list* gets used otherwise. - :param seed: Seed for internal random number generator. - :param lowest_energy_conformation: If True, we choose the lowest scoring - conformation of the trajectory. Otherwise, - the last accepted proposal. - - :type sampler: :ref:`mc-sampler-object` - :type closer: :ref:`mc-closer-object` - :type scorer: :ref:`mc-scorer-object` - :type cooler: :ref:`mc-cooler-object` - :type steps: :class:`int` - :type bb_list: :class:`~promod3.loop.BackboneList` - :type initialize: :class:`bool` - :type seed: :class:`int` - :type lowest_energy_conformation: :class:`bool` - +.. autofunction:: SampleMonteCarlo .. _mc-sampler-object: diff --git a/modelling/pymod/CMakeLists.txt b/modelling/pymod/CMakeLists.txt index 3f6e6f13be2079f47e8a11c1ffb979971cfb5534..5e44247d6bd03f7f80cbaee99fe973700ce4dcda 100644 --- a/modelling/pymod/CMakeLists.txt +++ b/modelling/pymod/CMakeLists.txt @@ -22,6 +22,7 @@ set(MODELLING_PYMOD _denovo.py _fragger_handle.py _reconstruct_sidechains.py + _monte_carlo.py ) pymod(NAME modelling diff --git a/modelling/pymod/__init__.py b/modelling/pymod/__init__.py index cdd41d485e453d2089b76a5a34c710559cd86a87..70ccfbb88eb891f54d21a19f01f99425b49d7f49 100644 --- a/modelling/pymod/__init__.py +++ b/modelling/pymod/__init__.py @@ -7,3 +7,4 @@ from _reconstruct_sidechains import * from _ring_punches import * from _denovo import * from _fragger_handle import * +from _monte_carlo import * diff --git a/modelling/pymod/_monte_carlo.py b/modelling/pymod/_monte_carlo.py new file mode 100644 index 0000000000000000000000000000000000000000..7c284aed00a979441c76a938f8282a0484cfb4aa --- /dev/null +++ b/modelling/pymod/_monte_carlo.py @@ -0,0 +1,90 @@ +import random +import math + + +def SampleMonteCarlo(sampler, closer, scorer, cooler, steps, bb_list, + initialize, seed = 0, lowest_energy_conformation = True): + ''' A convenient function to perform Monte Carlo sampling using a simulated + annealing scheme. In every iteration, a new loop conformation gets proposed + by the provided *sampler* and closed by the *closer*. Upon scoring, this new + conformation gets accepted/rejected using a metropolis criterion based on + the temperature given by the *cooler* + => acceptance probability: exp(-delta_score/T). + The result is stored in *bb_list* (passed by reference, so NO return value) + and is either the lowest energy conformation ever encountered or the last + accepted proposal. + + :param sampler: Sampler object capable of initializing and altering + conformations. + :param closer: Closer object to adapt a new conformation to + the environment. + :param scorer: Scorer object to score new loop conformations. + :param cooler: Cooler object to control the temperature of the + Monte Carlo trajectory. + :param steps: Number of Monte Carlo iterations to be performed. + :param bb_list: The chosen conformation gets stored here. + :param initialize: Whether a new bb_list should be generated as starting + point, based on the samplers Initialize function. + The input *bb_list* gets used otherwise. + :param seed: Seed for internal random number generator. + :param lowest_energy_conformation: If True, we choose the lowest scoring + conformation of the trajectory. + Otherwise, the last accepted proposal. + + :type sampler: :ref:`mc-sampler-object` + :type closer: :ref:`mc-closer-object` + :type scorer: :ref:`mc-scorer-object` + :type cooler: :ref:`mc-cooler-object` + :type steps: :class:`int` + :type bb_list: :class:`~promod3.loop.BackboneList` + :type initialize: :class:`bool` + :type seed: :class:`int` + :type lowest_energy_conformation: :class:`bool` + ''' + + # Initialize from scratch if necessary + if initialize: + initialized = False + for i in range(100): + sampler.Initialize(bb_list) + if closer.Close(bb_list, bb_list): + initialized = True + break + if not initialized: + raise RuntimeError("Failed to initialize monte carlo protocol!") + + # setup + random.seed(seed) + score = scorer.GetScore(bb_list) + closed = False + min_score = score + min_score_bb_list = bb_list.Copy() + proposed_bb_list = bb_list.Copy() + + for i in range(steps): + + # try several proposals (we might not be able to close at first attempt) + closed = False + for j in range(3): + sampler.ProposeStep(bb_list, proposed_bb_list) + if closer.Close(proposed_bb_list, proposed_bb_list): + closed = True + break + + # check for success, try again in next step in case of failure + if not closed: + continue + + # accept / reject based on Metropolis criterion + temperature = cooler.GetTemperature() + new_score = scorer.GetScore(proposed_bb_list) + + if math.exp(-(new_score-score)/temperature) > random.random(): + positions = proposed_bb_list.Copy() + score = new_score + if lowest_energy_conformation and score < min_score: + min_score = score + min_score_bb_list = positions + + if lowest_energy_conformation: + bb_list = min_score_bb_list diff --git a/modelling/pymod/export_monte_carlo.cc b/modelling/pymod/export_monte_carlo.cc index 5935c7534efd62b055bf57da7252e4f038685976..1cbf251a6f3114dc5f3c0cddd9b509b0dc1eb987 100644 --- a/modelling/pymod/export_monte_carlo.cc +++ b/modelling/pymod/export_monte_carlo.cc @@ -332,10 +332,4 @@ void export_monte_carlo() { register_ptr_to_python<NTerminalCloserPtr>(); register_ptr_to_python<CTerminalCloserPtr>(); register_ptr_to_python<MonteCarloCloserPtr>(); - - // export monte carlo sampler - def("SampleMonteCarlo", &SampleMonteCarlo, - (arg("sampler"), arg("closer"), arg("scorer"), arg("cooler"), - arg("steps"), arg("backbone"), arg("initialize")=true, arg("seed")=0, - arg("lowest_energy_conformation")=true)); }