# -*- coding: utf-8 -*-
# Copyright (c) 2019, 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,
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 shooting_move is None:
shooting_move = 'sh'
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_point1, ph_point2 = crossing_finder(path_old, interfaces[1])
if not ph_point1:
return False, path_old, 'NCR'
ss_interfaces = [interfaces[1], interfaces[1], interfaces[2]]
osc_try = 0 # One step crossing attempt counter
for _ in range(tis_settings['n_jumps']):
logger.debug('Trying a new stone skipping move')
# Here we choose between the two
# possible shooting points that describe a crossing.
if rgen.rand() >= 0.5:
shooting_point = ph_point1[-1]
else:
shooting_point = ph_point2[-1]
# 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 _ 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:
engine.modify_velocities(
shooting_point, 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(shooting_point,
order_function,
interfaces[1],
engine, rgen)
osc_try += 1
if success:
break
if not success: # In case we reached maxlen in jumps attempts.
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.
if path.get_end_point(interfaces[1]) == start_cond:
trial_path = path.empty_path(maxlen=maxlen-1)
success, _ = engine.propagate(trial_path,
shooting_point,
order_function,
ss_interfaces,
reverse=True) # Backward
new_segment = paste_paths(trial_path, path, overlap=True,
maxlen=maxlen)
else:
new_segment = path.empty_path(maxlen=maxlen)
new_segment.append(path.phasepoints[0])
success, _ = engine.propagate(new_segment,
path.phasepoints[1],
order_function,
ss_interfaces,
reverse=False) # Forward
if not success:
new_segment.status = 'XSS'
trial_path = new_segment
break
ph_point1, ph_point2 = crossing_finder(new_segment, interfaces[1])
if success:
success, trial_path, _ = extender(new_segment, order_function,
interfaces, engine, tis_settings)
if success:
success = ss_wt_acceptance(path_old, trial_path, interfaces, rgen,
order_function, tis_settings)
# 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 trial_path.get_start_point(interfaces[0], interfaces[2])\
== trial_path.get_end_point(interfaces[0], interfaces[2]):
if rgen.rand() < 0.5:
trial_path = trial_path.reverse(order_function)
trial_path.status = 'ACC'
trial_path.generated = ('ss', shooting_point.order[0],
osc_try, tis_settings['n_jumps'])
if not success:
return False, trial_path, trial_path.status
# Compute, eventually, the weights
if 'high_accept' in tis_settings and tis_settings['high_accept']:
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):
"""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.
Returns
-------
out[0] : boolean
True if the path can be accepted.
"""
if path_new.get_start_point(interfaces[0], interfaces[2]) == 'R' and \
path_new.get_end_point(interfaces[0], interfaces[2]) == 'R':
path_new.status = 'BWI'
return False
if path_old.generated[0] == 'ss':
if 'high_accept' in tis_settings and tis_settings['high_accept']:
if path_new.get_start_point(interfaces[0], interfaces[2]) == 'R' \
and \
path_new.get_end_point(interfaces[0], interfaces[2]) == 'L':
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_old.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.
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)
save_initial_op = source_segment.phasepoints[-2].order[0]
n_virtual, save_acc = 0, 0
side_picker0 = 'L' if rgen.rand() >= 0.5 else 'R'
side_picker = side_picker0
for i_j in range(tis_settings['n_jumps']):
# Choose a side
if side_picker != side_picker0 or i_j == n_virtual:
for _ in range(n_virtual):
if side_picker == 'L':
pre_shooting_point = source_segment.phasepoints[0]
shooting_point = source_segment.phasepoints[1]
reverse = False
else:
pre_shooting_point = source_segment.phasepoints[-1]
shooting_point = source_segment.phasepoints[-2]
reverse = True
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=reverse)
if segments_counter(new_segment,
tis_settings['interface_sour'],
interfaces[1]) == 1:
logger.debug('Web successful')
source_segment = new_segment.copy()
source_segment.status = 'ACC'
save_acc += 1
break
n_virtual = 0
side_picker0 = side_picker
else:
n_virtual += 1
side_picker = 'L' if rgen.rand() > 0.5 else 'R'
accept, trial_path, _ = extender(source_segment, order_function,
interfaces, engine, tis_settings)
# 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', save_initial_op,
segment_to_pick, save_acc)
return False, trial_path, trial_path.status
# Set the path flags
if path_old.get_move() == 'ld' and save_acc == 0:
trial_path.generated = ('ld', save_initial_op,
segment_to_pick, save_acc)
else:
trial_path.generated = ('wt', save_initial_op,
segment_to_pick, save_acc)
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