"""
Interface to the c++ core.
The data in c++ is dimensionless. The optimize function defined in this module
manages the conversion.
"""
from __future__ import print_function
import os
import pickle
import numpy as np
try:
import pydiffusioncma
except ImportError as e:
raise ImportError(
os.linesep*2 +
"Import error in dcma.optimize. "
"The C++ wrappers (and probably the binaries) "
"for the diffusioncma package (pydiffusioncma) "
"are not installed. Please check out "
"https://gitlab.com/Olllom/diffusioncma for instructions "
"on how to install diffusioncma properly."
)
import numpy as np
from dcma import tools, matrices, profiles, biasing
[docs]def optimize(transition_matrices, energy=None, diffusion=None,
dm="cos", em="cos", ds=6, es=9, cma_tolerance=None,
logfile="", bias=biasing.Bias()):
"""Optimize profiles using CMA-ES method.
Args:
transition_matrices (list of matrices.Transitions): The observed transitions
energy (np.array): Energy profile of shape (n_bins, ) if only D should be optimized (em="const")
diffusion (np.array): Diffusion profile of shape (nbins, ) if only F should be optimized (dm="const")
dm (str): Type of diffusion model (see dcma opt --help for available settings)
em (str): Type of energy model (see dcma opt --help for available settings)
ds (int): Number of coefficients for the diffusion profile
es (int): Number of coefficients for the energy profile
cma_tolerance (float): Tolerance for the optimization
logfile (str): Path to log file (default=no logging)
bias (biasing.Bias): A bias that was applied during simulation
Returns:
profiles.Profiles: the optimal profiles
"""
if logfile != "":
# write coefficient configuration for interpretation of coefficients
local_variables = locals()
with open(logfile + ".coco", 'wb') as fp:
pickle.dump(local_variables, fp)
if type(transition_matrices) is matrices.Transitions:
transition_matrices = [transition_matrices]
weights = np.array([t.weight for t in transition_matrices])
converter = tools.UnitConverter.from_transitions(transition_matrices)
if diffusion is not None:
diffusion = converter.normalize_diffusion(diffusion, apply_log=True)
if bias.get_biasing_type() is not "no":
print("Biasing with the following potential:")
print(bias.get_biasing_type(), bias.get_potential(converter))
kwargs = {
"tmats": [m.matrix for m in transition_matrices],
"lagtime": [m.lag_time for m in transition_matrices],
"ener": energy, "diff": diffusion,
"weight": weights,
"opt": "cma",
"dm": dm, "em": em, "ds": ds, "es": es, "cmat": cma_tolerance,
"profile_out": "",
"bias": bias.get_biasing_type(),
"biasing_potential": bias.get_potential(converter),
"logfile": logfile
}
# remove Nones
kwargs = {k: kwargs[k] for k in kwargs if kwargs[k] is not None}
solution = pydiffusioncma.optimize(**kwargs)
opt_free_energy = solution.getProfiles().getFreeEnergy()
opt_diffusion = solution.getProfiles().getDiffusion()
opt_diffusion = converter.denormalize_diffusion(opt_diffusion)
# this is a quick-fix: the optimization is only applied with interpolated diffusion profiles
# D_{i+0.5} = 0.5*(D_i + D_{i+1}), the optimal solution can thus have large oscillations.
# We apply a convolution to get rid of this effect.
opt_diffusion = 0.25*(
2*opt_diffusion
+ np.roll(opt_diffusion, 1)
+ np.roll(opt_diffusion, -1)
)
return profiles.Profiles(
converter.edges, opt_free_energy, opt_diffusion
)
[docs]class ProfileGenerator(object):
"""Generate Profiles from raw coefficients."""
def __init__(self, transition_matrices, energy=None, diffusion=None,
dm="cos", em="cos", ds=6, es=9):
if type(transition_matrices) is matrices.Transitions:
transition_matrices = [transition_matrices]
weights = np.array([t.weight for t in transition_matrices])
self.converter = tools.UnitConverter.from_transitions(transition_matrices)
if diffusion is not None:
diffusion = self.converter.normalize_diffusion(diffusion, apply_log=True)
kwargs = {
"tmats": [m.matrix for m in transition_matrices],
"lagtime": [m.lag_time for m in transition_matrices],
"ener": energy, "diff": diffusion,
"weight": weights,
"opt": "cma",
"dm": dm, "em": em, "ds": ds, "es": es, "cmat": None,
"profile_out": "",
"logfile": ""
}
# remove Nones
kwargs = {k: kwargs[k] for k in kwargs if kwargs[k] is not None}
self._generator = pydiffusioncma.make_profile_generator(**kwargs)
def __call__(self, coefficients):
raw_profile = self._generator.generate(coefficients)
opt_free_energy = raw_profile.getFreeEnergy()
opt_diffusion = raw_profile.getDiffusion()
opt_diffusion = self.converter.denormalize_diffusion(opt_diffusion)
return profiles.Profiles(
self.converter.edges, opt_free_energy, opt_diffusion
)
[docs] @staticmethod
def from_log(logfile):
"""
Create profile generator (and read data) from a log file.
Args:
logfile: logfile from a previous optimization run
Returns:
A tuple (generator, logs, best_coeff, best_f)
- generator: A ProfileGenerator instance
- logs (np.array): An array generated from the logs (time(s), objective, coeff1, coeff2, ...)
- best_coeff (np.array): The optimal coefficients
- best_f (float): The optimal objective function value
"""
coeff_config = os.path.expandvars(logfile+".coco")
assert os.path.isfile(logfile)
assert os.path.isfile(coeff_config)
with open(coeff_config, "rb") as fp:
coeff_config = pickle.load(fp)
coeff_config = {key: coeff_config[key]
for key in ["transition_matrices", "energy", "diffusion", "dm", "em", "ds", "es"]}
generator = ProfileGenerator(**coeff_config)
with open(logfile, "r") as fp:
lines = fp.readlines()
logs = np.genfromtxt(lines[1:-1])
best_coeff = np.genfromtxt(lines[-1].split()[2:])
best_f = float(lines[-1].split()[1])
return generator, logs, best_coeff, best_f