Commit 0d85a952 authored by eric pellegrini's avatar eric pellegrini
Browse files

Modified API for interpolation order configurator and related widget

docstringed InterpolationOrder configurator
parent 0d5930cd
......@@ -34,8 +34,6 @@ import wx
from MDANSE.Externals.pubsub import pub
from MDANSE.Mathematics.Signal import INTERPOLATION_ORDER
from MDANSE.App.GUI import DATA_CONTROLLER
from MDANSE.App.GUI.Framework.Widgets.IWidget import IWidget
......@@ -46,44 +44,24 @@ class InterpolationOrderWidget(IWidget):
def initialize(self):
pass
def add_widgets(self):
sizer = wx.BoxSizer(wx.HORIZONTAL)
label = wx.StaticText(self._widgetPanel, wx.ID_ANY, label="interpolation order")
self._interpolationOrder = wx.SpinCtrl(self._widgetPanel, min=-1, max=len(INTERPOLATION_ORDER)-1, initial=0, style=wx.SP_WRAP|wx.SP_ARROW_KEYS)
self._info = wx.StaticText(self._widgetPanel, wx.ID_ANY, label="")
self._interpolationOrder = wx.Choice(self._widgetPanel, wx.ID_ANY, choices=self.configurator.orders)
self._interpolationOrder.SetStringSelection(InterpolationOrderWidget.orders[self.configurator.default])
sizer.Add(label, 0, wx.ALL|wx.ALIGN_CENTER_VERTICAL, 5)
sizer.Add(self._interpolationOrder, 0, wx.ALL, 5)
sizer.Add(self._info, 0, wx.ALL|wx.ALIGN_CENTER_VERTICAL, 5)
self._interpolationOrder.Bind(wx.EVT_SPIN, self.on_select_interpolation_order)
pub.subscribe(self.set_wigdet_value, ("set_trajectory"))
self.select_interpolation_order(-1)
self._interpolationOrder.SetValue(self.configurator.orders[0])
return sizer
def select_interpolation_order(self, order):
self._interpolationOrder.SetValue(order)
if order == -1:
self._info.SetLabel("The velocities will be taken directly from the trajectory.")
else:
self._info.SetLabel("The velocities will be derived from coordinates using order-%d interpolation." % order)
def on_select_interpolation_order(self, event):
order = event.GetPosition()
self.select_interpolation_order(order)
def set_wigdet_value(self, message):
window, filename = message
......@@ -94,19 +72,12 @@ class InterpolationOrderWidget(IWidget):
trajectory = DATA_CONTROLLER[filename]
if "velocities" in trajectory.data.variables():
order = -1
self._interpolationOrder.SetStringSelection(self.configurator.orders)
else:
order = 0
self._interpolationOrder.SetRange(order,len(INTERPOLATION_ORDER)-1)
self.select_interpolation_order(order)
def get_widget_value(self):
value = self._interpolationOrder.GetValue()
self._interpolationOrder.SetStringSelection(self.configurator.orders[1:])
def get_widget_value(self):
value = self._interpolationOrder.GetStringSelection()
return value
def set_widget_value(self, value):
pass
......@@ -31,29 +31,46 @@ Created on Mar 30, 2015
'''
from MDANSE.Framework.Configurators.IConfigurator import ConfiguratorError
from MDANSE.Framework.Configurators.IntegerConfigurator import IntegerConfigurator
from MDANSE.Mathematics.Signal import INTERPOLATION_ORDER
from MDANSE.Framework.Configurators.SingleChoiceConfigurator import SingleChoiceConfigurator
class InterpolationOrderConfigurator(IntegerConfigurator):
class InterpolationOrderConfigurator(SingleChoiceConfigurator):
"""
This configurator allows to set as input the order of the interpolation apply when deriving velocities
from atomic coordinates to the atomic trajectories.
This configurator allows to input the interpolation order to be applied when deriving velocities
from atomic coordinates.
The value should be higher than 0 if velocities are not provided with the trajectory.
:note: this configurator depends on 'trajectory' configurator to be configured
"""
type = "interpolation_order"
_default = 0
_default = "no interpolation"
orders = ["no interpolation","1st order","2nd order","3rd order","4th order","5th order"]
def __init__(self, name, **kwargs):
'''
Initializes the configurator.
:param name: the name of the configurator as it will be appear in the configuration.
:type name: str.
'''
IntegerConfigurator.__init__(self, name, choices=range(-1,len(INTERPOLATION_ORDER)), **kwargs)
SingleChoiceConfigurator.__init__(self, name, choices=InterpolationOrderConfigurator.orders, **kwargs)
def configure(self, configuration, value):
'''
Configure the input interpolation order.
:param configuration: the current configuration.
:type configuration: a MDANSE.Framework.Configurable.Configurable object.
:param value: the interpolation order to be configured.
:type value: str one of 'no interpolation','1st order','2nd order','3rd order','4th order' or '5th order'.
'''
IntegerConfigurator.configure(self, configuration, value)
SingleChoiceConfigurator.configure(self, configuration, value)
if value == -1:
if value == "no interpolation":
trajConfig = configuration[self._dependencies['trajectory']]
......@@ -64,11 +81,4 @@ class InterpolationOrderConfigurator(IntegerConfigurator):
else:
self["variable"] = "configuration"
def get_information(self):
if self["value"] == -1:
return "No velocities interpolation"
else:
return "Velocities interpolated from atomic coordinates"
\ No newline at end of file
self["variable"] = "configuration"
\ No newline at end of file
......@@ -112,12 +112,12 @@ class DensityOfStates(IJob):
last=self.configuration['frames']['last']+1,
step=self.configuration['frames']['step'],
variable=self.configuration['interpolation_order']["variable"])
order = self.configuration["interpolation_order"]["value"]
if order != -1:
val = self.configuration["interpolation_order"]["value"]
if val != "no interpolation":
for axis in range(3):
series[:,axis] = differentiate(series[:,axis], order=order, dt=self.configuration['frames']['time_step'])
series[:,axis] = differentiate(series[:,axis], order=val, dt=self.configuration['frames']['time_step'])
if self.configuration["projection"]["projector"] is not None:
series = self.configuration['projection']["projector"](series)
......
......@@ -31,9 +31,7 @@ Created on Apr 10, 2015
'''
import collections
import os
from MDANSE import PREFERENCES, REGISTRY
from MDANSE.Framework.Jobs.IJob import IJob
from MDANSE.Mathematics.Signal import get_spectrum
......@@ -88,8 +86,8 @@ class StructureFactorFromScatteringFunction(IJob):
self._outputData["s(q,f)_%s" % suffix][:] = get_spectrum(v.getValue(),
self.configuration["instrument_resolution"]["time_window"],
self.configuration["instrument_resolution"]["time_step"],
axis=1)
axis=1)
def run_step(self, index):
"""
Runs a single step of the job.\n
......@@ -100,10 +98,8 @@ class StructureFactorFromScatteringFunction(IJob):
#. index (int): The index of the step.
"""
return index, None
def combine(self, index, x):
"""
Combines returned results of run_step.\n
......@@ -113,7 +109,6 @@ class StructureFactorFromScatteringFunction(IJob):
"""
pass
def finalize(self):
"""
Finalizes the calculations (e.g. averaging the total term, output files creations ...).
......@@ -122,5 +117,3 @@ class StructureFactorFromScatteringFunction(IJob):
self._outputData.write(self.configuration['output_files']['root'], self.configuration['output_files']['formats'], self.header)
self.configuration['netcdf_input_file']['instance'].close()
#MDANSE : Molecular Dynamics Analysis for Neutron Scattering Experiments
#------------------------------------------------------------------------------------------
#Copyright (C)
#2015- Eric C. Pellegrini Institut Laue-Langevin
#BP 156
#6, rue Jules Horowitz
#38042 Grenoble Cedex 9
#France
#pellegrini[at]ill.fr
#goret[at]ill.fr
#aoun[at]ill.fr
#
#This library is free software; you can redistribute it and/or
#modify it under the terms of the GNU Lesser General Public
#License as published by the Free Software Foundation; either
#version 2.1 of the License, or (at your option) any later version.
#
#This library is distributed in the hope that it will be useful,
#but WITHOUT ANY WARRANTY; without even the implied warranty of
#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
#Lesser General Public License for more details.
#
#You should have received a copy of the GNU Lesser General Public
#License along with this library; if not, write to the Free Software
#Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
'''
Created on May 21, 2015
@author: Eric C. Pellegrini
'''
import collections
from MDANSE import ELEMENTS
from MDANSE.Framework.Jobs.IJob import IJob
from MDANSE.Mathematics.Arithmetic import weight
from MDANSE.Mathematics.Signal import correlation, differentiate, normalize
from MDANSE.MolecularDynamics.Trajectory import read_atoms_trajectory
class VelocityAutoCorrelationFunction(IJob):
"""
The Velocity AutoCorrelation Function (VACF) is a property describing
the dynamics of a molecular system. It reveals the underlying nature of the forces acting
on the system by computing the cartesian density of states for a set of atoms.
In a molecular system that would be made of non interacting particles, the velocities would
be constant at any time triggering the VACF to be a constant value. Now, if we think about a
system with small interactions such as in a gas-phase, the magnitude and direction of the velocity
of a particle will change gradually over time due to its collision with the other particles of the
molecular system. In such a system, the VACF will be represented by a decaying exponential.
In the case of solid phase, the interaction are much stronger and, as a results, the atoms
are bound to a given position from which they will move backwards and forwards oscillating
between positive and negative values of their velocity. The oscillations will not be of equal
magnitude however, but will decay in time, because there are still perturbative forces acting on
the atoms to disrupt the perfection of their oscillatory motion. So, in that case the VACF will
look like a damped harmonic motion.
Finally, in the case of liquid phase, the atoms have more freedom than in solid phase and
because of the diffusion process, the oscillatory motion seen in solid phase will be cancelled quite
rapidly depending on the density of the system. So, the VACF will just have one very damped
oscillation before decaying to zero. This decaying time can be considered as the average time
for a collision between two atoms to occur before they diffuse away.
"""
type = 'vacf'
label = "Velocity AutoCorrelation Function"
category = ('Dynamics',)
ancestor = "mmtk_trajectory"
settings = collections.OrderedDict()
settings['trajectory'] = ('mmtk_trajectory',{})
settings['frames'] = ('frames', {'dependencies':{'trajectory':'trajectory'}})
settings['interpolation_order'] = ('interpolation_order', {'label':"velocities",
'dependencies':{'trajectory':'trajectory'}})
settings['projection'] = ('projection', {'label':"project coordinates"})
settings['normalize'] = ('boolean', {'default':False})
settings['atom_selection'] = ('atom_selection',{'dependencies':{'trajectory':'trajectory',
'grouping_level':'grouping_level'}})
settings['grouping_level'] = ('grouping_level',{})
settings['atom_transmutation'] = ('atom_transmutation',{'dependencies':{'trajectory':'trajectory',
'atom_selection':'atom_selection'}})
settings['weights'] = ('weights',{})
settings['output_files'] = ('output_files', {'formats':["netcdf","ascii"]})
settings['running_mode'] = ('running_mode',{})
def initialize(self):
"""
Initialize the input parameters and analysis self variables
"""
self.numberOfSteps = self.configuration['atom_selection']['n_groups']
# Will store the time.
self._outputData.add("time","line", self.configuration['frames']['time'], units='ps')
# Will store the mean square displacement evolution.
for element in self.configuration['atom_selection']['contents'].keys():
self._outputData.add("vacf_%s" % element,"line", (self.configuration['frames']['number'],), axis="time", units="nm2/ps2")
self._outputData.add("vacf_total","line", (self.configuration['frames']['number'],), axis="time", units="nm2/ps2")
def run_step(self, index):
"""
Runs a single step of the job.\n
:Parameters:
#. index (int): The index of the step.
:Returns:
#. index (int): The index of the step.
#. atomicDOS (numpy.array): The calculated density of state for atom of index=index
#. atomicVACF (numpy.array): The calculated velocity auto-correlation function for atom of index=index
"""
# get atom index
indexes = self.configuration['atom_selection']["groups"][index]
series = read_atoms_trajectory(self.configuration["trajectory"]["instance"],
indexes,
first=self.configuration['frames']['first'],
last=self.configuration['frames']['last']+1,
step=self.configuration['frames']['step'],
variable=self.configuration['interpolation_order']["variable"])
val = self.configuration["interpolation_order"]["value"]
if val != "no interpolation":
for axis in range(3):
series[:,axis] = differentiate(series[:,axis], order=val, dt=self.configuration['frames']['time_step'])
series = self.configuration['projection']["projector"](series)
atomicVACF = correlation(series,axis=0,reduce=1)
return index, atomicVACF
def combine(self, index, x):
"""
Combines returned results of run_step.\n
:Parameters:
#. index (int): The index of the step.\n
#. x (any): The returned result(s) of run_step
"""
# The symbol of the atom.
element = self.configuration['atom_selection']['elements'][index][0]
self._outputData["vacf_%s" % element] += x
def finalize(self):
"""
Finalizes the calculations (e.g. averaging the total term, output files creations ...).
"""
# The MSDs per element are averaged.
for element, number in self.configuration['atom_selection']['n_atoms_per_element'].items():
self._outputData["vacf_%s" % element] /= number
if self.configuration['normalize']["value"]:
for element in self.configuration['atom_selection']['n_atoms_per_element'].keys():
self._outputData["vacf_%s" % element] = normalize(self._outputData["vacf_%s" % element], axis=0)
props = dict([[k,ELEMENTS[k,self.configuration["weights"]["property"]]] for k in self.configuration['atom_selection']['n_atoms_per_element'].keys()])
vacfTotal = weight(props,
self._outputData,
self.configuration['atom_selection']['n_atoms_per_element'],
1,
"vacf_%s")
self._outputData["vacf_total"][:] = vacfTotal
self._outputData.write(self.configuration['output_files']['root'], self.configuration['output_files']['formats'], self.header)
self.configuration['trajectory']['instance'].close()
\ No newline at end of file
......@@ -7,34 +7,34 @@ from MDANSE.Core.Error import Error
class SignalError(Error):
pass
INTERPOLATION_ORDER = []
INTERPOLATION_ORDER = {}
INTERPOLATION_ORDER.append(numpy.array([[-3., 4.,-1.],
[-1., 0., 1.],
[ 1.,-4., 3.]], dtype=numpy.float64))
INTERPOLATION_ORDER["1st order"] = numpy.array([[-3., 4.,-1.],
[-1., 0., 1.],
[ 1.,-4., 3.]], dtype=numpy.float64)
INTERPOLATION_ORDER.append(numpy.array([[-3., 4.,-1.],
[-1., 0., 1.],
[ 1.,-4., 3.]], dtype=numpy.float64))
INTERPOLATION_ORDER["2nd order"] = numpy.array([[-3., 4.,-1.],
[-1., 0., 1.],
[ 1.,-4., 3.]], dtype=numpy.float64)
INTERPOLATION_ORDER.append(numpy.array([[-11., 18., -9., 2.],
[ -2., -3., 6.,-1.],
[ 1., -6., 3., 2.],
[ -2., 9.,-18.,11.]], dtype=numpy.float64))
INTERPOLATION_ORDER["3rd order"] = numpy.array([[-11., 18., -9., 2.],
[ -2., -3., 6.,-1.],
[ 1., -6., 3., 2.],
[ -2., 9.,-18.,11.]], dtype=numpy.float64)
INTERPOLATION_ORDER.append(numpy.array([[-50., 96.,-72., 32.,-6.],
[ -6.,-20., 36.,-12., 2.],
[ 2.,-16., 0., 16.,-2.],
[ -2., 12.,-36., 20., 6.],
[ 6.,-32., 72.,-96.,50.]], dtype=numpy.float64))
INTERPOLATION_ORDER["4th order"] = numpy.array([[-50., 96.,-72., 32.,-6.],
[ -6.,-20., 36.,-12., 2.],
[ 2.,-16., 0., 16.,-2.],
[ -2., 12.,-36., 20., 6.],
[ 6.,-32., 72.,-96.,50.]], dtype=numpy.float64)
INTERPOLATION_ORDER.append(numpy.array([[-274., 600.,-600., 400.,-150., 24.],
[ -24.,-130., 240.,-120., 40., -6.],
[ 6., -60., -40., 120., -30., 4.],
[ -4., 30.,-120., 40., 60., -6.],
[ 6., -40., 120.,-240., 130., 24.],
[ -24., 150.,-400., 600.,-600.,274.]], dtype=numpy.float64))
INTERPOLATION_ORDER["5th order"] = numpy.array([[-274., 600.,-600., 400.,-150., 24.],
[ -24.,-130., 240.,-120., 40., -6.],
[ 6., -60., -40., 120., -30., 4.],
[ -4., 30.,-120., 40., 60., -6.],
[ 6., -40., 120.,-240., 130., 24.],
[ -24., 150.,-400., 600.,-600.,274.]], dtype=numpy.float64)
def correlation(x, y=None, axis=0, sumOverAxis=None, average=None):
......@@ -96,19 +96,19 @@ def normalize(x, axis=0):
return nx
def differentiate(a, dt=1.0, order=0):
def differentiate(a, dt=1.0, order="1st order"):
try:
coefs = INTERPOLATION_ORDER[order]
except TypeError:
raise SignalError("Invalid type for differentiation order")
except IndexError:
if not INTERPOLATION_ORDER.has_key(order):
raise SignalError("Invalid differentiation order")
coefs = INTERPOLATION_ORDER[order]
# outputSeries is the output resulting from the differentiation
ts = numpy.zeros(a.shape, dtype=numpy.float)
fact = 1.0/dt
if order == 0:
if order == "1st order":
ts[0] = numpy.add.reduce(coefs[0,:]*a[:3])
ts[-1] = numpy.add.reduce(coefs[2,:]*a[-3:])
......@@ -116,8 +116,10 @@ def differentiate(a, dt=1.0, order=0):
gj = a[1:] - a[:-1]
ts[1:-1] = (gj[1:]+gj[:-1])
fact /= 2.0
# Case of the order 2
elif order == 1:
elif order == "2nd order":
ts[0] = numpy.add.reduce(coefs[0,:]*a[:3])
ts[-1] = numpy.add.reduce(coefs[2,:]*a[-3:])
......@@ -128,8 +130,10 @@ def differentiate(a, dt=1.0, order=0):
gj[:,2] = coefs[1,2]*a[2:]
ts[1:-1] = numpy.add.reduce(gj,-1)
fact /= 2.0
# Case of the order 3
elif order == 2:
elif order == "3rd order":
# Special case for the first and last elements
ts[0] = numpy.add.reduce(coefs[0,:]*a[:4])
......@@ -144,8 +148,10 @@ def differentiate(a, dt=1.0, order=0):
gj[:,3] = coefs[2,3]*a[3:]
ts[2:-1] = numpy.add.reduce(gj,-1)
fact /= 6.0
# Case of the order 4
elif order == 3:
elif order == "4th order":
# Special case for the first and last elements
ts[0] = numpy.add.reduce(coefs[0,:]*a[:5])
......@@ -162,8 +168,10 @@ def differentiate(a, dt=1.0, order=0):
gj[:,4] = coefs[2,4]*a[4:]
ts[2:-2] = numpy.add.reduce(gj,-1)
fact /= 24.0
# Case of the order 5
elif order == 4:
elif order == "5th order":
# Special case for the first and last elements
ts[0] = numpy.add.reduce(coefs[0,:]*a[:6])
......@@ -182,9 +190,8 @@ def differentiate(a, dt=1.0, order=0):
gj[:,5] = coefs[3,5]*a[5:]
ts[3:-2] = numpy.add.reduce(gj,-1)
fact /= 120.0
fact = 1.0/(dt*math.factorial(max(order+1,2)))
ts *= fact
return ts
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment