Duplicate SampleMonteCarlo Function to also have it in Python

Advantage:
Makes it possible to provide sampler/closer/scorer/cooler objects implemented
in both, Python and C++
Disadvantage:
Code duplication as we also need that functionality directly in C++

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:
 

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

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 *

modelling/pymod/_monte_carlo.py

+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

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));
 }