Source code for pyretis.core.tis

# -*- coding: utf-8 -*-
# Copyright (c) 2021, PyRETIS Development Team.
# Distributed under the LGPLv2.1+ License. See LICENSE for more info.
"""This file contains functions used in TIS.

This module defines the functions needed to perform TIS simulations.
The algorithms are implemented as described by van Erp et al. [TIS]_
and Riccardi et al. [SS+WT]_.

Methods defined here
~~~~~~~~~~~~~~~~~~~~

make_tis_step (:py:func:`.make_tis_step`)
    A method that will perform a single TIS step.

make_tis_step_ensemble (:py:func:`.make_tis_step_ensemble`)
    A method to perform a TIS step for a path ensemble. It will handle
    adding of the path to a path ensemble object.

shoot (:py:func:`.shoot`)
    A method that will perform a shooting move.

select_shoot (:py:func:`.select_shoot`)
    A method that randomly selects the shooting point.

make_path (:py:func:`.make_path`)
    A method that will generate a full path by combining
    back and forward integration.

stone_skipping (:py:func:`.stone_skipping`)
    A method that will do a stone skipping shooting move.

web_throwing (:py:func:`.web_throwing`)
    A method that will do a web throwing shooting move.

extender (:py:func:`.extender`)
    A method that will extend a path to target interfaces.

time_reversal (:py:func:`.time_reversal`)
    A method for performing the time reversal move.

compute_weight_ss (:py:func:`.compute_weight_ss`)
    A method to compute the statistical weight of a path generated by a
    Stone Skipping move.

ss_wt_acceptance (:py:func:`.ss_wt_acceptance`)
    A method to check the acceptance of a newly generated path
    according to super detailed balance rules for Stone Skipping (basic and
    High Acceptance version) and Web Throwing.


References
~~~~~~~~~~
.. [SS+WT] Enrico Riccardi, Oda Dahlen, Titus S. van Erp,
   J. Phys. Chem. letters, 8, 18, 4456, (2017),
   https://doi.org/10.1021/acs.jpclett.7b01617

.. [TIS] Titus S. van Erp, Daniele Moroni and Peter G. Bolhuis,
   J. Chem. Phys. 118, 7762 (2003),
   https://dx.doi.org/10.1063%2F1.1562614

"""
import logging
from pyretis.core.path import Path, paste_paths
from pyretis.core.common import (segments_counter,
                                 counter,
                                 select_and_trim_a_segment,
                                 crossing_counter,
                                 crossing_finder)
from pyretis.core.montecarlo import metropolis_accept_reject
logger = logging.getLogger(__name__)  # pylint: disable=invalid-name
logger.addHandler(logging.NullHandler())


__all__ = ['make_tis_step_ensemble',
           'make_tis_step',
           'select_shoot',
           'shoot',
           'make_path',
           'stone_skipping',
           'web_throwing',
           'compute_weight_ss',
           'ss_wt_acceptance',
           'extender',
           'time_reversal']


[docs]def make_tis_step_ensemble(path_ensemble, order_function, engine, rgen, tis_settings, cycle): """Perform a TIS step in a given path ensemble. This function will run :py:func:`.make_tis_step` for the given path ensemble and handling adding of the path to the path ensemble. Parameters ---------- path_ensemble : object like :py:class:`.PathEnsemble` This is the path ensemble to perform the TIS step for. order_function : object like :py:class:`.OrderParameter` The class used for obtaining the order parameter(s). engine : object like :py:class:`.EngineBase` The engine to use for propagating a path. rgen : object like :py:class:`.RandomGenerator` This is the random generator that will be used. tis_settings : dict This dictionary contains the TIS settings. cycle : int The current cycle number. Returns ------- out[0] : boolean True if the new path can be accepted, False otherwise. out[1] : object like :py:class:`.PathBase` The generated path. out[2] : string The status of the path. """ start_cond = path_ensemble.start_condition logger.info('TIS move in: %s', path_ensemble.ensemble_name) engine.exe_dir = path_ensemble.directory['generate'] if 'shooting_moves' in tis_settings and tis_settings['shooting_moves']: ens_num = path_ensemble.ensemble_number shooting_move = tis_settings['shooting_moves'][ens_num] else: shooting_move = 'sh' accept, trial, status = make_tis_step(path_ensemble.last_path, order_function, path_ensemble.interfaces, engine, rgen, tis_settings, start_cond, shooting_move) if accept: logger.info('The move was accepted!') else: logger.info('The move was rejected! (%s)', status) path_ensemble.add_path_data(trial, status, cycle=cycle) return accept, trial, status
[docs]def make_tis_step(path, order_function, interfaces, engine, rgen, tis_settings, start_cond, shooting_move=None): """Perform a TIS step and generate a new path. The new path will be generated from the input path, either by performing a time-reversal move or by a shooting move. This is determined pseudo-randomly by drawing a random number from a uniform distribution using the given random generator. Parameters ---------- path : object like :py:class:`.PathBase` This is the input path which will be used for generating a new path. order_function : object like :py:class:`.OrderParameter` The class used for obtaining the order parameter(s). interfaces : list of floats These are the interface positions on ordered form ``[left, middle, right]``. engine : object like :py:class:`.EngineBase` The engine to use for propagating a path when integrating the equations of motion. rgen : object like :py:class:`.RandomGenerator` Random number generator used to determine what TIS move to perform. tis_settings : dict This dictionary contains the settings for the TIS method. Here we explicitly use: * `freq`: float, the frequency of how often we should do time reversal moves. start_cond : string The starting condition for the path. This is determined by the ensemble we are generating for - it is 'R'ight or 'L'eft. shooting_move : string The flag of the shooting move to operate. Returns ------- out[0] : boolean True if the new path can be accepted. out[1] : object like :py:class:`.PathBase` The generated path. out[2] : string The status of the path. """ if rgen.rand() < tis_settings['freq']: logger.info('Performing a time reversal move') accept, new_path, status = time_reversal(path, order_function, interfaces, start_cond) else: logger.info('Selecting the shooting move.') accept, new_path, status = select_shoot(path, order_function, interfaces, engine, rgen, tis_settings, shooting_move, start_cond) return accept, new_path, status
[docs]def select_shoot(path, order_function, interfaces, engine, rgen, tis_settings, shooting_move, start_cond): """Select the shooting move to generate a new path. The new path will be generated from the input path, either by performing a normal shooting or web-throwing. This is determined pseudo-randomly by drawing a random number from a uniform distribution using the given random generator. Parameters ---------- path : object like :py:class:`.PathBase` This is the input path which will be used for generating a new path. order_function : object like :py:class:`.OrderParameter` The class used for obtaining the order parameter(s). interfaces : list of floats These are the interface positions on ordered form ``[left, middle, right]``. engine : object like :py:class:`.EngineBase` The engine to use for propagating a path when integrating the equations of motion. rgen : object like :py:class:`.RandomGenerator` Random number generator used to determine what TIS move to perform. tis_settings : dict This dictionary contains the settings for the TIS method. Here we explicitly use: * `freq`: float, the frequency of how often we should do time reversal moves. shooting_move : string The flag of the shooting move to operate. start_cond : string The starting condition for the path. This is determined by the ensemble we are generating for - it is 'R'ight or 'L'eft. Returns ------- out[0] : boolean True if the new path can be accepted. out[1] : object like :py:class:`.PathBase` The generated path. out[2] : string The status of the path. """ if shooting_move == 'wt': logger.info('Performing a web throwing move') accept, new_path, status = web_throwing(path, order_function, interfaces, engine, rgen, tis_settings) elif shooting_move == 'ss': logger.info('Performing a stone skipping move') accept, new_path, status = stone_skipping(path, order_function, interfaces, engine, rgen, tis_settings, start_cond) else: logger.info('Performing a shooting move') accept, new_path, status = shoot(path, order_function, interfaces, engine, rgen, tis_settings, start_cond) return accept, new_path, status
[docs]def time_reversal(path, order_function, interfaces, start_condition): """Perform a time-reversal move. Parameters ---------- path : object like :py:class:`.PathBase` This is the input path which will be used for generating a new path. order_function : object like :py:class:`.OrderParameter` The class used for obtaining the order parameter(s). interfaces : list/tuple of floats These are the interface positions of the form ``[left, middle, right]``. start_condition : string The starting condition, 'L'eft or 'R'ight. Returns ------- out[0] : boolean True if the path can be accepted. out[1] : object like :py:class:`.PathBase` The time reversed path. out[2] : string Status of the path, this is one of the strings defined in `.path._STATUS`. """ new_path = path.reverse(order_function) start, _, _, _ = new_path.check_interfaces(interfaces) # Explicitly set how this was generated. Keep the old if loaded. if path.get_move() == 'ld': new_path.set_move('ld') else: new_path.generated = ('tr', 0, 0, 0) if start == start_condition: accept = True status = 'ACC' else: accept = False status = 'BWI' # Backward trajectory end at wrong interface. new_path.status = status return accept, new_path, status
def prepare_shooting_point(path, order_function, engine, rgen, tis_settings): """Select and modify velocities for a shooting move. This method will randomly select a shooting point from a given path and modify its velocities. Parameters ---------- path : object like :py:class:`.PathBase` This is the input path which will be used for generating a new path. order_function : object like :py:class:`.OrderParameter` The class used to calculate the order parameter. engine : object like :py:class:`.EngineBase` The engine to use for propagating a path. rgen : object like :py:class:`.RandomGenerator` This is the random generator that will be used. tis_settings : dict This contains the settings for TIS. Here, we use the settings which dictates how we modify the velocities. Returns ------- out[0] : object like :py:class:`.System` The shooting point with modified velocities. out[1] : integer The index of the shooting point in the original path. out[2] : float The change in kinetic energy when modifying the velocities. """ shooting_point, idx = path.get_shooting_point() orderp = shooting_point.order logger.info('Shooting from order parameter/index: %f, %d', orderp[0], idx) # Copy the shooting point, so that we can modify velocities without # altering the original path: system = shooting_point.copy() # Modify the velocities: dek, _, = engine.modify_velocities( system, rgen, sigma_v=tis_settings['sigma_v'], aimless=tis_settings['aimless'], momentum=tis_settings['zero_momentum'], rescale=tis_settings['rescale_energy'], ) orderp = engine.calculate_order(order_function, system) system.order = orderp system.idx = idx return system, idx, dek def one_step_crossing(shooting_point, order_function, interface, engine, rgen): """Create a path of one step and check the crossing with the interface. This function will do a single step to try to cross the interface. Note that a step might involve several substeps, depending on the input file selection/sampling strategy. This task is the Achilles` Heel of Stone Skipping. If not wisely done, a large number of attempts to cross the interface in one step will be done, destroying the sampling efficiency. Parameters ---------- shooting_point : object like :py:class:`.System` The shooting point with modified velocities. order_function : object like :py:class:`.OrderParameter` The class used to calculate the order parameter. interface : float This is the interface position to be crossed. engine : object like :py:class:`.EngineBase` The engine to use for propagating a path. rgen : object like :py:class:`.RandomGenerator` This is the random generator that will be used. Returns ------- out[0] : boolean True if the path can be accepted. out[1] : object like :py:class:`.PathBase` Returns the generated path. """ # The trial path we need to generate. Note, 1 step = 2 points trial_path = Path(rgen=rgen, maxlen=2) engine.propagate(trial_path, shooting_point, order_function, interfaces=[-float("inf"), float("inf"), float("inf")]) if crossing_counter(trial_path, interface) == 0: return False, trial_path return True, trial_path def check_kick(shooting_point, interfaces, trial_path, rgen, dek, tis_settings): """Check the modification of the shooting point. After generating velocities for a shooting point, we do some additional checking to see if the shooting point is acceptable. Parameters ---------- shooting_point : object like :py:class:`.System` The shooting point with modified velocities. interfaces : list of floats The interfaces used for TIS, in the format ``[left, middle, right]``. trial_path : object like :py:class:`.PathBase` The path we are currently generating. rgen : object like :py:class:`.RandomGenerator` This is the random generator that will be used to check if we accept the shooting point based on the change in kinetic energy. dek : float The change in kinetic energy when modifying the velocities. tis_settings : dict This contains the settings for TIS. Returns ------- out : boolean True if the kick was OK, False otherwise. """ # 1) Check if the kick was too violent: left, _, right = interfaces if not left <= shooting_point.order[0] < right: # Shooting point was velocity dependent and was kicked outside # of boundaries when modifying velocities. trial_path.append(shooting_point) trial_path.status = 'KOB' return False # 2) If the kick is not aimless, we check if we reject it or not: if not tis_settings['aimless']: accept_kick = metropolis_accept_reject(rgen, shooting_point, dek) # If one wish to implement a bias call, this can be done here. if not accept_kick: trial_path.append(shooting_point) trial_path.status = 'MCR' # Momenta Change Rejection. return False return True def shoot_backwards(path_back, trial_path, shooting_point, engine, order_function, interfaces, tis_settings, start_cond): """Shoot in the backward time direction. Parameters ---------- path_back : object like :py:class:`.PathBase` The path we will fill with phase points from the propagation. trial_path : object like :py:class:`.PathBase` The current trial path generated by the shooting. shooting_point : object like :py:class:`.System` The shooting point used as the starting point for the propagation. engine : object like :py:class:`.EngineBase` The engine to use for propagating a path. order_function : object like :py:class:`.OrderParameter` The class used to calculate the order parameter. interfaces : list/tuple of floats These are the interface positions of the form ``[left, middle, right]``. tis_settings : dict This contains the settings for TIS. start_cond : string The starting condition for the current ensemble, 'L'eft or 'R'ight. Returns ------- out : boolean True if the backward path was generated successfully, False otherwise. """ logger.debug('Propagating backwards for the shooting move.') path_back.time_origin = trial_path.time_origin success_back, _ = engine.propagate(path_back, shooting_point, order_function, interfaces, reverse=True) if not success_back: # Something went wrong, most probably the path length was exceeded. trial_path.status = 'BTL' # BTL = backward trajectory too long. # Add the failed path to trial path for analysis: trial_path += path_back if path_back.length >= tis_settings['maxlength'] - 1: # BTX is backward trajectory longer than maximum memory. trial_path.status = 'BTX' return False # Backward seems OK so far, check if the ending point is correct: left, _, right = interfaces if path_back.get_end_point(left, right) != start_cond: # Nope, backward trajectory end at wrong interface. trial_path.status = 'BWI' trial_path += path_back # Store path for analysis. return False return True
[docs]def shoot(path, order_function, interfaces, engine, rgen, tis_settings, start_cond): """Perform a shooting-move. This function will perform the shooting move from a randomly selected time-slice. Parameters ---------- path : object like :py:class:`.PathBase` This is the input path which will be used for generating a new path. order_function : object like :py:class:`.OrderParameter` The class used to calculate the order parameter. interfaces : list/tuple of floats These are the interface positions of the form ``[left, middle, right]``. engine : object like :py:class:`.EngineBase` The engine to use for propagating a path. rgen : object like :py:class:`.RandomGenerator` This is the random generator that will be used. tis_settings : dict This contains the settings for TIS. Keys used here: * `aimless`: boolean, is the shooting aimless or not? * `allowmaxlength`: boolean, should paths be allowed to reach maximum length? * `maxlength`: integer, maximum allowed length of paths. start_cond : string The starting condition for the current ensemble, 'L'eft or 'R'ight. Returns ------- out[0] : boolean True if the path can be accepted. out[1] : object like :py:class:`.PathBase` Returns the generated path. out[2] : string Status of the path, this is one of the strings defined in :py:const:`.path._STATUS`. """ trial_path = path.empty_path() # The trial path we will generate. shooting_point, _, dek = prepare_shooting_point( path, order_function, engine, rgen, tis_settings ) # Store info about this point, just in case we have to return # before completing a full new path: trial_path.generated = ('sh', shooting_point.order[0], shooting_point.idx, 0) trial_path.time_origin = path.time_origin + shooting_point.idx # We now check if the kick was OK or not: if not check_kick(shooting_point, interfaces, trial_path, rgen, dek, tis_settings): return False, trial_path, trial_path.status # OK: kick was either aimless or it was accepted by Metropolis # we should now generate trajectories, but first check how long # it should be (if the path comes from a load, it is assumed to not # respect the detail balance anyway): if path.get_move() == 'ld' or tis_settings['allowmaxlength']: maxlen = tis_settings['maxlength'] else: maxlen = min(int((path.length - 2) / rgen.rand()) + 2, tis_settings['maxlength']) success, trial_path, _ = make_path(trial_path, order_function, interfaces, engine, tis_settings, maxlen, shooting_point, start_cond) # Deal with the rejections in the backward paths if not success: # Deal with the rejections in the forward paths if trial_path.status == 'FTL': # If we reached this point, the backward path was successful, # but the forward was not. For the case where the forward was # also successful, the length of the trial path cannot exceed # the maximum length given in the TIS settings. Thus we only # need to check this here, i.e. when given that the backward # was successful and the forward not: if trial_path.length >= tis_settings['maxlength']: trial_path.status = 'FTX' # exceeds "memory". return False, trial_path, trial_path.status # Deal with the rejections for path properties. # Make sure we did not hit the left interface on {0-} # Which is the only ensemble that allows paths starting in R if trial_path.check_interfaces(interfaces)[:2] == ('R', 'L'): trial_path.status = '0-L' return False, trial_path, trial_path.status # Last check - Did we cross the middle interface? if not trial_path.check_interfaces(interfaces)[-1][1]: # No, we did not cross the middle interface: trial_path.status = 'NCR' return False, trial_path, trial_path.status trial_path.status = 'ACC' return True, trial_path, trial_path.status
[docs]def make_path(path, order_function, interfaces, engine, tis_settings, maxlen, shooting_point, start_cond): """Generate a full path from a shooting point. Forth and Back in time. This function will perform the shooting move from a randomly selected time-slice. Parameters ---------- path : object like :py:class:`.PathBase` This is the input path which will be used for generating a new path. order_function : object like :py:class:`.OrderParameter` The class used to calculate the order parameter. interfaces : list/tuple of floats These are the interface positions of the form ``[left, middle, right]``. engine : object like :py:class:`.EngineBase` The engine to use for propagating a path. tis_settings : dict This contains the settings for TIS. Keys used here: * `aimless`: boolean, is the shooting aimless or not? * `allowmaxlength`: boolean, should paths be allowed to reach maximum length? * `maxlength`: integer, maximum allowed length of paths. maxlen : integer The maximum length of the path in respect to further constrains (e.g. detail balance). shooting_point: object like :py:class:`.phasepoint` The frame from which the new path is going to be generated. start_cond : string The starting condition for the current ensemble, 'L'eft or 'R'ight. Returns ------- out[0] : boolean True if the path can be accepted. out[1] : object like :py:class:`.PathBase` Returns the generated path. out[2] : string Status of the path, this is one of the strings defined in :py:const:`.path._STATUS`. """ path_back = path.empty_path(maxlen=maxlen-1) # Since the forward path must be at least one step, the maximum # length for the backward path is maxlen-1. # Generate the backward path: if not shoot_backwards(path_back, path, shooting_point, engine, order_function, interfaces, tis_settings, start_cond): return False, path, path.status # Everything seems fine, now propagate forward. # Note that the length of the forward path is adjusted to # account for the fact that it shares a point with the backward # path (i.e. the shooting point). The duplicate point is just # counted once when the paths are merged by the method # `paste_paths` by setting `overlap=True` (which indicates that # the forward and backward paths share a point). path_forw = path.empty_path(maxlen=(maxlen - path_back.length + 1)) logger.debug('Propagating forwards for shooting move...') success_forw, _ = engine.propagate(path_forw, shooting_point, order_function, interfaces, reverse=False) path_forw.time_origin = path.time_origin # Now, the forward propagation could have failed by exceeding the # maximum length for the forward path. However, it could also fail # when we paste together so that the length is larger than the # allowed maximum. We paste first and ask later: trial_path = paste_paths(path_back, path_forw, overlap=True, maxlen=tis_settings['maxlength']) # Also update information about the shooting: trial_path.generated = ('sh', shooting_point.order[0], shooting_point.idx, path_back.length - 1) if not success_forw: trial_path.status = 'FTL' return False, trial_path, trial_path.status trial_path.status = 'ACC' return True, trial_path, trial_path.status
[docs]def stone_skipping(path_old, order_function, interfaces, engine, rgen, tis_settings, start_cond): """Perform a stone_skipping move. This function will perform the famous stone skipping move from an initial path. Parameters ---------- path_old : object like :py:class:`.PathBase` This is the input path which will be used for generating a new path. order_function : object like :py:class:`.OrderParameter` The class used to calculate the order parameter. interfaces : list/tuple of floats These are the interface positions of the form ``[left, middle, right]``. engine : object like :py:class:`.EngineBase` The engine to use for propagating a path. rgen : object like :py:class:`.RandomGenerator` This is the random generator that will be used. tis_settings : dict This contains the settings for TIS. Keys used here: * `aimless`: boolean, is the shooting aimless or not? * `allowmaxlength`: boolean, should paths be allowed to reach maximum length? * `maxlength`: integer, maximum allowed length of paths. * `high_accept`: boolean, the option for High Acceptance SS. start_cond : string The starting condition for the current ensemble, 'L'eft or 'R'ight. Returns ------- out[0] : boolean True if the path can be accepted. out[1] : object like :py:class:`.PathBase` Returns the generated path. out[2] : string Status of the path, this is one of the strings defined in :py:const:`.path._STATUS`. """ maxlen = tis_settings['maxlength'] ph_pt1, ph_pt2 = crossing_finder(path_old, interfaces[1]) if not ph_pt1: return False, path_old, 'NCR' ss_interfaces = [interfaces[1], interfaces[1], interfaces[2]] osc_try = 0 # One step crossing attempt counter for i in range(tis_settings['n_jumps']): logger.debug(f'Trying a new stone skipping move, jump {i}') # Here we choose between the two # possible shooting points that describe a crossing. sh_pt = ph_pt1[-1] if rgen.rand() >= 0.5 else ph_pt2[-1] engine.dump_phasepoint(sh_pt, str(counter()) + '_ss_shoot') # To continue, we must be sure that the new path # CROSSES the interface in ONLY ONE step. # Generate paths until it succeed. That is # what makes this version of the SS move useless for large systems. for j in range(maxlen): # This function can become actually fun to work on. # e.g. have a 50% chance to give random v for each particle # Modify the velocities: logger.debug(f'jump{i}, try {j}, start: {sh_pt.order[0]}') engine.modify_velocities( sh_pt, rgen, sigma_v=tis_settings['sigma_v'], aimless=tis_settings['aimless'], momentum=tis_settings['zero_momentum'], rescale=tis_settings['rescale_energy'], ) # A path of two frames is going to be generated. success, path = one_step_crossing(sh_pt, order_function, interfaces[1], engine, rgen) osc_try += 1 if success: break else: # In case we reached maxlength in jumps attempts. success = False path.status = 'NSS' trial_path = path break # Depending on the shooting point (before or after the interface), # a backward path or a continuation has to be generated. new_segment = path.empty_path(maxlen=maxlen-1) if path.get_end_point(interfaces[1], interfaces[2]) == start_cond: path = path.reverse(order_function) new_segment.append(path.phasepoints[0].copy()) success, _ = engine.propagate(new_segment, path.phasepoints[1].copy(), order_function, ss_interfaces) if not success: new_segment.status = 'XSS' trial_path = new_segment break ph_pt1, ph_pt2 = crossing_finder(new_segment, interfaces[1]) if success: success, trial_path, _ = extender(new_segment, order_function, interfaces, engine, tis_settings) trial_path.generated = ('ss', sh_pt.order[0], osc_try, trial_path.length) if success: success = ss_wt_acceptance(path_old, trial_path, interfaces, rgen, order_function, tis_settings) reverse = False # A to A trajectories can change orientation at the end. 50% chances. # If a path goes from A to B, it we shall not reverse the path. # B to B paths are already rejected. if success and rgen.rand() < 0.5: if tis_settings.get('high_accept', False): if trial_path.get_start_point(interfaces[0], interfaces[2])\ == trial_path.get_end_point(interfaces[0], interfaces[2]): reverse = True else: reverse = True if reverse: trial_path = trial_path.reverse(order_function) trial_path.generated = ('ss', sh_pt.order[0], osc_try, trial_path.length) if success and trial_path.get_start_point(interfaces[0], interfaces[2]) != start_cond: trial_path.status = 'BWI' success = False logger.debug(f'SS move {trial_path.status}') if not success: return False, trial_path, trial_path.status # Compute, eventually, the weights if tis_settings.get('high_accept', False): trial_path.weight = compute_weight_ss(trial_path, interfaces) else: trial_path.weight = 1. trial_path.status = 'ACC' return True, trial_path, trial_path.status
[docs]def compute_weight_ss(path, interfaces): """Compute the path weight after a Stone Skipping move. This function computes the weights that will be used in the computation of the P cross. This trick allows the use of the High Acceptance version of Stone Skipping, allowing the acceptance of B to A paths. The drawback is that swapping moves needs to account also for this different weights. Parameters ---------- path : object like :py:class:`.PathBase` This is the input path which will be checked. interfaces : list/tuple of floats These are the interface positions of the form ``[left, middle, right]``. Returns ------- out[0] : float The weight of the path. """ weight = 1.*crossing_counter(path, interfaces[1]) if path.get_end_point(interfaces[0], interfaces[2]) == 'L': weight *= 0.5 return weight
[docs]def ss_wt_acceptance(path_old, path_new, interfaces, rgen, order_function, tis_settings, start_cond='L'): """Accept or reject the path_new. Super detailed balance rule is used in the original version and in the High Acceptance one. In the regular version, P acc = min (1, Cold/Cnew), where for Stone Skipping C is crossing, for Web Throwing C is segment. In the High Acceptance version, P acc = 1 for Stone Skipping and it also allows to accept paths that go from B to A, by reversing it. To respect super detailed balance, the weights have to be changed accordingly. This is done elsewhere. Parameters ---------- path_old : object like :py:class:`.PathBase` This is the previous path. path_new : object like :py:class:`.PathBase` This is the new path that might get accepted. interfaces : list/tuple of floats These are the interface positions of the form ``[left, middle, right]``. rgen : object like :py:class:`.RandomGenerator` This is the random generator that will be used. order_function : object like :py:class:`.OrderParameter` The class used to calculate the order parameter. tis_settings : dict This contains the settings for TIS. Keys used here: * `high_accept` : boolean, the option for High Acceptance SS. start_cond : string, optional The starting condition for the current ensemble, 'L'eft or 'R'ight. Returns ------- out[0] : boolean True if the path can be accepted. """ if start_cond != path_new.get_start_point(interfaces[0], interfaces[2]) \ == path_new.get_end_point(interfaces[0], interfaces[2]): path_new.status = 'BWI' return False if path_new.generated[0] == 'ss': if tis_settings.get('high_accept', False): if path_new.get_end_point(interfaces[0], interfaces[2]) != start_cond: path_new = path_new.reverse(order_function) else: cr_old = crossing_counter(path_old, interfaces[1]) cr_new = crossing_counter(path_new, interfaces[1]) if rgen.rand() >= min(1.0, cr_old/cr_new): path_new.status = 'SSA' return False elif path_new.generated[0] == 'wt': sour_int = tis_settings['interface_sour'] cr_old = segments_counter(path_old, sour_int, interfaces[1]) cr_new = segments_counter(path_new, sour_int, interfaces[1]) if rgen.rand() >= min(1.0, cr_old / cr_new): path_new.status = 'WTA' return False path_new.status = 'ACC' return True
[docs]def web_throwing(path_old, order_function, interfaces, engine, rgen, tis_settings): """Perform a web_throwing move. This function performs the great web throwing move from an initial path. Parameters ---------- path_old : object like :py:class:`.PathBase` This is the input path which will be used for generating a new path. order_function : object like :py:class:`.OrderParameter` The class used to calculate the order parameter. interfaces : list/tuple of floats These are the interface positions of the form ``[left, middle, right]``. engine : object like :py:class:`.EngineBase` The engine to use for propagating a path. rgen : object like :py:class:`.RandomGenerator` This is the random generator that will be used. tis_settings : dict This contains the settings for TIS. Keys used here: * `aimless`: boolean, is the shooting aimless or not? * `allowmaxlength`: boolean, should paths be allowed to reach maximum length? * `maxlength`: integer, maximum allowed length of paths. * `high_accept`: boolean, the option for High Acceptance SS. Returns ------- out[0] : boolean True if the path can be accepted. out[1] : object like :py:class:`.PathBase` Returns the generated path. out[2] : string Status of the path, this is one of the strings defined in :py:const:`.path._STATUS`. """ if not interfaces[0] < tis_settings['interface_sour'] <= interfaces[1]: raise ValueError('SOUR interface is not correctly positioned') maxlen = tis_settings['maxlength'] ccnt = segments_counter(path_old, tis_settings['interface_sour'], interfaces[1]) if ccnt == 0: return False, path_old, 'NSG' segment_to_pick = int(rgen.rand() * ccnt) wt_interfaces = [tis_settings['interface_sour'], tis_settings['interface_sour'], interfaces[1]] source_segment = select_and_trim_a_segment(path_old, tis_settings['interface_sour'], interfaces[1], segment_to_pick) shoots, save_acc = [0], 0 key = rgen.rand() >= 0.5 # Start from a random side for _ in range(tis_settings['n_jumps']): if rgen.rand() >= 0.5: shoots[-1] += 1 # One more on the Same side else: shoots.append(1) # A move in the other side for n_virtual in shoots: key = not key # Change side, key controls also path reverse for _ in range(n_virtual): if key: pre_shooting_point = source_segment.phasepoints[-1] shooting_point = source_segment.phasepoints[-2] else: pre_shooting_point = source_segment.phasepoints[0] shooting_point = source_segment.phasepoints[1] prefix = str(counter()) engine.dump_phasepoint(pre_shooting_point, prefix + '_wt_pre_shoot') engine.dump_phasepoint(shooting_point, prefix + '_wt_shoot') new_segment = path_old.empty_path(maxlen=maxlen) new_segment.append(pre_shooting_point) logger.debug('Trying a new web') engine.propagate(new_segment, shooting_point, order_function, wt_interfaces, reverse=key) if segments_counter(new_segment, tis_settings['interface_sour'], interfaces[1], reverse=key) == 1: logger.debug('Web successful') source_segment = new_segment.reverse( order_function) if key else new_segment source_segment.status = 'ACC' save_acc += 1 break accept, trial_path, _ = extender(source_segment, order_function, interfaces, engine, tis_settings) trial_path.generated = ('wt', source_segment.phasepoints[1].order[0], save_acc, trial_path.length) # Check that we did not get a B to A or a B to B path. if accept: accept = ss_wt_acceptance(path_old, trial_path, interfaces, rgen, order_function, tis_settings) if not accept: trial_path.generated = ('wt', source_segment.phasepoints[1].order[0], save_acc, trial_path.length) return False, trial_path, trial_path.status logger.debug(f'WT move {trial_path.status}') # Set the path flags if path_old.get_move() == 'ld' and save_acc == 0: trial_path.generated = ('ld', source_segment.phasepoints[1].order[0], save_acc, trial_path.length) trial_path.weight = 1. trial_path.status = 'ACC' return True, trial_path, trial_path.status
[docs]def extender(source_segment, order_function, interfaces, engine, tis_settings): """Extend a path to the given interfaces. This function will perform the web throwing move from an initial path. Parameters ---------- source_segment : object like :py:class:`.PathBase` This is the input path which will be prolonged. order_function : object like :py:class:`.OrderParameter` The class used to calculate the order parameter. interfaces : list/tuple of floats These are the interface positions of the form ``[left, middle, right]``. engine : object like :py:class:`.EngineBase` The engine to use for propagating a path. tis_settings : dict This contains the settings for TIS. Keys used here: * `aimless`: boolean, is the shooting aimless or not? * `allowmaxlength`: boolean, should paths be allowed to reach maximum length? * `maxlength`: integer, maximum allowed length of paths. Returns ------- out[0] : boolean True if the path can be accepted. out[1] : object like :py:class:`.PathBase` Returns the generated path. out[2] : string Status of the path, this is one of the strings defined in :py:const:`.path._STATUS`. """ maxlen = tis_settings['maxlength'] # Extender shooting_point = source_segment.phasepoints[0] if interfaces[0] <= shooting_point.order[0] < interfaces[-1]: back_segment = source_segment.empty_path(maxlen=maxlen) logger.debug('Trying to extend backwards') engine.propagate(back_segment, shooting_point, order_function, interfaces, reverse=True) if back_segment.length >= tis_settings['maxlength']: back_segment.status = 'BTX' # exceeds "memory". return False, back_segment, back_segment.status trial_path = paste_paths(back_segment, source_segment, overlap=True, maxlen=tis_settings['maxlength']) else: trial_path = source_segment.copy() shooting_point = trial_path.phasepoints[-1] if interfaces[0] <= shooting_point.order[0] < interfaces[-1]: trial_path.delete(-1) engine.propagate(trial_path, shooting_point, order_function, interfaces, reverse=False) if trial_path.length >= tis_settings['maxlength']: trial_path.status = 'FTX' # exceeds "memory". return False, trial_path, trial_path.status trial_path.status = 'ACC' return True, trial_path, trial_path.status