Source code for bowhead.velocity

from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import *
import numpy as np

[docs]class Model: """Wound velocity model. The model is based on Gaussian Process Regression (GPR) and numerical differentiation. The GPR uses scikit-learn GPR class ``sklearn.gaussian_process.GaussianProcessRegressor``. Parameters ---------- area_kwargs : dictionary The keyword arguments to the area GPR model. perimeter_kwargs : dictionary The keyword arguments to the perimeter GPR model. Attributes ----------- area : area model Instances of GaussianProcessRegressor from scikit-learn. perimeter : perimeter model Instances of GaussianProcessRegressor from scikit-learn. time : sequence Sequence of the input time points. Defined when fitted. Notes ----- The keyword argument ``kernel`` is specifying the covariance function of the GPR models handling the area and perimeter. The default kernel is a sum of a linear and squared dotproduct kernel. See `scikit-learn <>`_ for the different kernels that can be used. Note that the kernel hyperparameters are optimized during fitting. """ def __init__(self, area_kwargs = None, perimeter_kwargs = None): kernel = ConstantKernel() + 1*DotProduct()**2 kwargs = {'kernel': kernel, 'n_restarts_optimizer': 5, 'normalize_y': True} if area_kwargs is None: area_kwargs = kwargs if perimeter_kwargs is None: perimeter_kwargs = kwargs self.area = GaussianProcessRegressor(**area_kwargs) self.perimeter = GaussianProcessRegressor(**perimeter_kwargs) self.time = None
[docs] def fit(self, wounds): """Fitting the velocity model to wound data. Parameters ---------- wounds : sequence of dictionaries A sequence of wound dictionaries representing an experimental time series of the wound healing assay, as returned by :func:`detect`. The area, perimeter and time point of each wound is used to fit the overall velocity model. All these attributes should be positive scalars. Returns ------- self : Model The velocity model in a fitted state. """ wounds = sorted(wounds, key=lambda x: x['time']) self.time = np.array([x['time'] for x in wounds]) time = self.time.reshape(-1,1) self.area.alpha = np.array([x['area_variance'] for x in wounds]), np.array([x['area'] for x in wounds])) self.perimeter.alpha = np.array([x['perimeter_variance'] for x in wounds]), np.array([x['perimeter'] for x in wounds])) return self
[docs] def predict(self, time, dt=.1, return_std=False): """Predicting a velocity curve of a wound healing experiment. In addition to the mean of the predictive distribution, also its standard deviation (return_std=True) can be requested. Parameters ----------- time : sequence Desired time points of the velocity prediction. Should be positive. dt : scalar The time interval for calculating velocity. return_std : bool Whether to include the standard deviation of the prediction or not. Returns ------- velocity, [vel_std] : 1d arrays Velocity and (optionally) it's standard deviation. Returned either one array or (if return_std=True) as a tuple (velocity, vel_std). """ time = np.asarray(time, dtype=float).round(4) linear_time = np.arange(time.min(), time.max(), dt, dtype=float).round(4) if len(linear_time) < 10: raise ValueError('choose a smaller dt for numeric stability') ptime = np.union1d(time, linear_time) area, a_std = self.area.predict(ptime.reshape(-1,1), return_std=True) peri, p_std = self.perimeter.predict(ptime.reshape(-1,1), return_std=True) da = np.gradient(area) / np.gradient(ptime) resolution = np.mean(np.gradient(self.time)) da_std = np.sqrt(2)*a_std / resolution vel = -da / peri frac_q = (da_std/da)**2 + (p_std/peri)**2 vel_std = np.abs(vel) * np.sqrt(frac_q) mask = np.in1d(ptime, time) velocity = vel[mask][np.argsort(time)] v_std = vel_std[mask][np.argsort(time)] if return_std: return velocity, v_std return velocity