"""Parameter optimization."""
# Authors: Carolina Fernandez <cxf418@miami.edu>
# Nick Tolley <nicholas_tolley@brown.edu>
# Ryan Thorpe <ryan_thorpe@brown.edu>
# Mainak Jas <mjas@mgh.harvard.edu>
import numpy as np
from .objective_functions import _rmse_evoked, _maximize_psd
from ..externals.mne import _validate_type
[docs]
class Optimizer:
def __init__(
self,
initial_net,
tstop,
constraints,
set_params,
initial_params=None,
solver="bayesian",
obj_fun="dipole_rmse",
max_iter=200,
):
"""Parameter optimization.
Parameters
----------
initial_net : instance of Network
The network object.
tstop : float
The simulated dipole's duration.
constraints : dict
The user-defined constraints.
set_params : func
User-defined function that sets parameters in network drives.
``set_params(net, params) -> None``
where ``net`` is a Network object and ``params`` is a dictionary
of the parameters that will be set inside the function.
initial_params : dict, optional
Initial parameters for the objective function. Keys are parameter
names, values are initial parameters. The default is None.
If None, the parameters will be set to the midpoints of parameter ranges.
solver : str
The optimizer, 'bayesian' or 'cobyla'.
obj_fun : str | func
The objective function to be minimized. Can be 'dipole_rmse',
'maximize_psd', or a user-defined function. The default is
'dipole_rmse'.
max_iter : int, optional
The max number of calls to the objective function. The default is
200.
Attributes
----------
constraints : dict
The user-defined constraints.
initial_params : dict
Initial parameters for the objective function.
max_iter : int
The max number of calls to the objective function.
solver : func
The optimization function.
obj_fun : func
The objective function to be minimized.
obj_fun_name : str
The name of the template objective function.
tstop : float
The simulated dipole's duration.
net_ : instance of Network
The network object with optimized drives.
obj_ : list
The objective function values.
opt_params_ : list
The list of optimized parameter values.
"""
if initial_net.external_drives:
raise ValueError(
"The current Network instance has external "
+ "drives, provide a Network object with no "
+ "external drives."
)
self._initial_net = initial_net
self.constraints = constraints
self._set_params = set_params
self.max_iter = max_iter
# Optimizer method
if solver == "bayesian":
self.solver = "bayesian"
self._assemble_constraints = _assemble_constraints_bayesian
self._run_opt = _run_opt_bayesian
elif solver == "cobyla":
self.solver = "cobyla"
self._assemble_constraints = _assemble_constraints_cobyla
self._run_opt = _run_opt_cobyla
else:
raise ValueError("solver must be 'bayesian' or 'cobyla'")
# Response to be optimized
if obj_fun == "dipole_rmse":
self.obj_fun = _rmse_evoked
self.obj_fun_name = "dipole_rmse"
elif obj_fun == "maximize_psd":
self.obj_fun = _maximize_psd
self.obj_fun_name = "maximize_psd"
else:
self.obj_fun = obj_fun # user-defined function
self.obj_fun_name = None
# Initial weights for the objective function
if initial_params is not None:
# Check if initial_params is not a dict
_validate_type(initial_params, dict)
# Check if initial_params keys match constraints keys
if set(initial_params.keys()) != set(self.constraints.keys()):
raise ValueError(
"The keys of initial_params must match the keys of constraints."
)
# Check if initial_params values contain only int or float
for val in initial_params.values():
_validate_type(val, (int, float))
# Check if initial_params values are within constraints
for key in initial_params:
if not (
self.constraints[key][0]
<= initial_params[key]
<= self.constraints[key][1]
):
raise ValueError(
f"Initial parameter '{key}' with value {initial_params[key]} is out of bounds {self.constraints[key]}."
)
self.initial_params = initial_params
else:
self.initial_params = _get_initial_params(self.constraints)
self.tstop = tstop
self.net_ = None
self.obj_ = list()
self.opt_params_ = None
[docs]
def __repr__(self):
is_fit = False
if self.net_ is not None:
is_fit = True
name = self.__class__.__name__
return f"<{name}\nsolver={self.solver}\nfit={is_fit}>"
[docs]
def fit(self, **obj_fun_kwargs):
"""Runs optimization routine.
Parameters
----------
target : instance of Dipole (Required if obj_fun='dipole_rmse')
A dipole object with experimental data.
n_trials : int (Optional if obj_fun='dipole_rmse')
Number of trials to simulate and average.
f_bands : list of tuples (Required if obj_fun='maximize_psd')
Lower and higher limit for each frequency band.
relative_bandpower : list of float | float (Required if obj_fun='maximize_psd')
Weight for each frequency band in f_bands. If a single float is provided,
the same weight is applied to all frequency bands.
scale_factor : float, optional
The dipole scale factor.
smooth_window_len : float, optional
The smooth window length.
"""
if self.obj_fun_name == "dipole_rmse":
if "target" not in obj_fun_kwargs:
raise Exception("target must be specified")
elif "n_trials" not in obj_fun_kwargs:
obj_fun_kwargs["n_trials"] = 1
elif self.obj_fun_name == "maximize_psd" and (
"f_bands" not in obj_fun_kwargs
or "relative_bandpower" not in obj_fun_kwargs
):
raise Exception("f_bands and relative_bandpower must be specified")
constraints = self._assemble_constraints(self.constraints)
opt_params, obj, net_ = self._run_opt(
self._initial_net,
self.tstop,
constraints,
self._set_params,
self.initial_params,
self.obj_fun,
self.max_iter,
obj_fun_kwargs,
)
self.net_ = net_
self.obj_ = obj
self.opt_params_ = opt_params
[docs]
def plot_convergence(self, ax=None, show=True):
"""Convergence plot.
Parameters
----------
ax : instance of matplotlib figure, optional
The matplotlib axis. The default is None.
show : bool
If True, show the figure. The default is True.
Returns
-------
fig : instance of plt.fig
The matplotlib figure handle.
"""
import matplotlib as mpl
import matplotlib.pyplot as plt
if ax is None:
fig, ax = plt.subplots(constrained_layout=True)
axis = ax if isinstance(ax, mpl.axes._axes.Axes) else ax
x = list(range(1, self.max_iter + 1))
y_min = min(self.obj_) - 0.01
y_max = max(self.obj_) + 0.01
axis.plot(x, self.obj_, color="black")
axis.set_ylim([y_min, y_max])
axis.set_title("Convergence")
axis.set_xlabel("Number of calls")
axis.set_ylabel("Objective value")
axis.grid(visible=True)
if ax is None:
fig.show(show)
return axis.get_figure()
else:
return ax.get_figure()
def _get_initial_params(constraints):
"""Gets initial parameters as midpoints of parameter ranges.
Parameters
----------
constraints : dict
The user-defined constraints.
Returns
-------
initial_params : dict
Keys are parameter names, values are initial parameters.
"""
initial_params = dict()
for cons_key in constraints:
initial_params.update(
{cons_key: (constraints[cons_key][0] + constraints[cons_key][1]) / 2}
)
return initial_params
def _assemble_constraints_bayesian(constraints):
"""Assembles constraints in format required by gp_minimize.
Parameters
----------
constraints : dict
The user-defined constraints.
Returns
-------
cons_bayesian : list of tuples
Lower and higher limit for each parameter.
"""
# assemble constraints in solver-specific format
cons_bayesian = list(constraints.values())
return cons_bayesian
def _assemble_constraints_cobyla(constraints):
"""Assembles constraints in format required by fmin_cobyla.
Parameters
----------
constraints : dict
The user-defined constraints.
Returns
-------
cons_bayesian : dict
Set of functions.
"""
# assemble constraints in solver-specific format
cons_cobyla = list()
for idx, cons_key in enumerate(constraints):
cons_cobyla.append(lambda x, idx=idx: float(constraints[cons_key][1]) - x[idx])
cons_cobyla.append(lambda x, idx=idx: x[idx] - float(constraints[cons_key][0]))
return cons_cobyla
def _update_params(initial_params, predicted_params):
"""Update params with predicted parameters.
Parameters
----------
initial_params : dict
Keys are parameter names, values are initial parameters.
predicted_params : list
Parameters selected by the optimizer.
Returns
-------
params : dict
Keys are parameter names, values are parameters.
"""
params = dict()
for param_key, param_name in enumerate(initial_params):
params.update({param_name: predicted_params[param_key]})
return params
def _run_opt_bayesian(
initial_net,
tstop,
constraints,
set_params,
initial_params,
obj_fun,
max_iter,
obj_fun_kwargs,
):
"""Runs optimization routine with gp_minimize optimizer.
Parameters
----------
initial_net : instance of Network
The network object.
tstop : float
The simulated dipole's duration.
constraints : list of tuples
Parameter constraints in solver-specific format.
set_params : func
User-defined function that sets parameters in network drives.
initial_params : dict
Keys are parameter names, values are initial parameters.
obj_fun : func
The objective function.
max_iter : int
Number of calls the optimizer makes.
obj_fun_kwargs : dict
Additional arguments along with their respective values to be passed
to the objective function.
Returns
-------
opt_params : list
Optimized parameters.
obj : list
Objective values.
net_ : instance of Network
Optimized network object.
"""
from ..externals.bayesopt import bayes_opt, expected_improvement
obj_values = list()
def _obj_func(predicted_params):
return obj_fun(
initial_net=initial_net,
initial_params=initial_params,
set_params=set_params,
predicted_params=predicted_params,
update_params=_update_params,
obj_values=obj_values,
tstop=tstop,
obj_fun_kwargs=obj_fun_kwargs,
)
opt_results, _ = bayes_opt(
func=_obj_func,
x0=list(initial_params.values()),
cons=constraints,
acquisition=expected_improvement,
maxfun=max_iter,
)
# get optimized params
opt_params = opt_results
# get objective values
obj = [np.min(obj_values[:idx]) for idx in range(1, max_iter + 1)]
# get optimized net
params = _update_params(initial_params, opt_params)
net_ = initial_net.copy()
set_params(net_, params)
return opt_params, obj, net_
def _run_opt_cobyla(
initial_net,
tstop,
constraints,
set_params,
initial_params,
obj_fun,
max_iter,
obj_fun_kwargs,
):
"""Runs optimization routine with fmin_cobyla optimizer.
Parameters
----------
initial_net : instance of Network
The network object.
tstop : float
The simulated dipole's duration.
constraints : list of tuples
Parameter constraints in solver-specific format.
set_params : func
User-defined function that sets parameters in network drives.
initial_params : dict
Keys are parameter names, values are initial parameters.
obj_fun : func
The objective function.
max_iter : int
Number of calls the optimizer makes.
obj_fun_kwargs : dict
Additional arguments along with their respective values to be passed
to the objective function.
Returns
-------
opt_params : list
Optimized parameters.
obj : list
Objective values.
net_ : instance of Network
Optimized network object.
"""
from scipy.optimize import fmin_cobyla
obj_values = list()
def _obj_func(predicted_params):
return obj_fun(
initial_net=initial_net,
initial_params=initial_params,
set_params=set_params,
predicted_params=predicted_params,
update_params=_update_params,
obj_values=obj_values,
tstop=tstop,
obj_fun_kwargs=obj_fun_kwargs,
)
opt_results = fmin_cobyla(
_obj_func,
cons=constraints,
rhobeg=0.1,
rhoend=1e-4,
x0=list(initial_params.values()),
maxfun=max_iter,
catol=0.0,
)
# get optimized params
opt_params = opt_results
# get objective values
obj = [np.min(obj_values[:idx]) for idx in range(1, max_iter + 1)]
# get optimized net
params = _update_params(initial_params, opt_params)
net_ = initial_net.copy()
set_params(net_, params)
return opt_params, obj, net_
