Source code for gdsctools.logistics

# -*- python -*-
# -*- coding utf-8 -*-
#  This file is part of GDSCTools software
#  Copyright (c) 2015 - Wellcome Trust Sanger Institute
#  All rights reserved
#  File author(s): Thomas Cokelaer <cokelaer@gmail.comWE HERE>
#  Distributed under the BSD 3-Clause License.
#  See accompanying file LICENSE.txt distributed with this software
#  website:
"""Sets of miscellaneous tools"""
import math
import warnings

import numpy as np
import pylab

__all__ = ['Logistic', 'LogisticMatchedFiltering']

def nextpower(x):
    """Returns next power of 2 times 2


        >>> nextpower(9)
        >>> nextpower(16)
        >>> nextpower(17)

    return int(2**(math.ceil(pylab.log2(x)) + 1 ))

[docs]class Logistic(object): """Simple logistic class to see the curve implied by xmid/scale parameters .. plot:: :include-source: from gdsctools.logistics import Logistic from pylab import legend tl = Logistic(2, 1) tl.plot() tl.scale = 4 tl.plot(hold=True) legend(['scale=1', 'scale=4']) The X values are set automatically given the number of data points :attr:`N` and the minimum and maximum X values. By default a sensible values for the minimun and maximum x values are guessed based on :attr:`scale` parameter but one can set the :attr:`xmin` and :attr:`xmax` .. plot:: :include-source: :width: 80% from gdsctools.logistics import Logistic from pylab import legend tl = Logistic(2,1) tl.plot() tl.xmin= -4 tl.xmax = 4 tl.N = 40 tl.plot(hold=True) legend(['default', 'user defined']) """ def __init__(self, xmid, scale, Asym=1, N=9, increase=False): r""".. rubric:: Constructor :param xmid: the first logistic function parameter :param scale: the second logistic function parameter :param Asym: Amplitude of the function :param N: number of data point :param increase: increasing or decreasing function. :param xmin: starting range of x-values :param xmax: ending range of x-values if :attr:`increase` is False: .. math:: L(x) = \frac{Asym}{1 + exp ((xmid-X)/scale)} if :attr:`increase` is True .. math:: L(x) = 1 - \frac{Asym}{1 + exp ((xmid-X)/scale)} Changing :attr:`xmin`, :attr:`xmax` or :attr:`N` does change the content of X. You can change :attr:`X` directly. """ self.xmid = xmid self._scale = None self.scale = scale self.Asym = Asym self._N = N self.increase = increase self._xmin = None self._xmax = None self._update_x() def _get_scale(self): return self._scale def _set_scale(self, scale): self._scale = scale scale = property(_get_scale, _set_scale) def _update_x(self): """Return a sensible linear space of X values """ dX = abs(self.scale) if self.xmin is None: self._xmin = self.xmid - dX * self.N if self.xmax is None: self._xmax = self.xmid + dX * self.N self._X = np.linspace(self.xmin, self.xmax, self.N)
[docs] def plot(self, hold=False): """Plot the logistic function""" if hold is False: pylab.clf() pylab.plot(self.X, self.Y, 'o-') pylab.grid(True) pylab.xlabel('X') pylab.ylabel('Y') pylab.axvline(self.xmid, color='r', linestyle='--')
def _get_xmin(self): return self._xmin def _set_xmin(self, xmin): self._xmin = xmin self._update_x() xmin = property(_get_xmin, _set_xmin, doc="set/get the minimum x range") def _get_xmax(self): return self._xmax def _set_xmax(self, xmax): self._xmax = xmax self._update_x() xmax = property(_get_xmax, _set_xmax, doc="set/get the maximum x range") def _get_N(self): return self._N def _set_N(self,N): self._N = N self._update_x() N = property(_get_N, _set_N, doc="set/get the number points") def _get_x(self): return self._X def _set_x(self, X): self._xmin = min(X) self._xmax = max(X) self._N = len(X) self._X = np.array(X) # no need to call _update X = property(_get_x, _set_x, doc="get/set of the x-values") def _get_y(self): with warnings.catch_warnings(): warnings.simplefilter("ignore") if self.increase is True: Y = self.Asym / (1. + np.exp((self.xmid - self.X)/self.scale)) else: Y = 1. - self.Asym / (1. + np.exp((self.xmid - self.X)/self.scale)) return Y Y = property(_get_y, doc='Getter for Y-values') def __str__(self): txt = "xmid: %s\n" % self.xmid txt += "scale: %s\n" % self.scale txt += "xmin: %(xmin)s\n" % {'xmin':self.xmin} txt += "xmax: %(xmax)s" % {'xmax':self.xmax} return txt
[docs]class LogisticMatchedFiltering(object): """Experimental class to identify parameters of a noisy Logistic function This class implements two methods to identify the parameters (xmid and scale) of a logistic function. Note One constraint is to define the x-values. .. plot:: :include-source: :width: 80% from gdsctools import logistics import pylab mf = logistics.LogisticMatchedFiltering(1,2) mf.scan(pylab.linspace(-3,3, 10), pylab.linspace(0,5,10)) """ def __init__(self, xmid, scale, N=9): """.. rubric:: constructor :param xmid: :param scale: :param N: """ # for book keeping and plotting self.signal = Logistic(xmid=xmid, scale=scale, N=N) # will be used = self.signal.Y # will be used for matchefd filtering self.template = Logistic(xmid=xmid, scale=scale, N=N) self.X = pylab.linspace(-8,8, N) def _get_N(self): return self.signal.N def _set_N(self, N): self.signal.N = N self.template.N = N = self.signal.Y N = property(_get_N, _set_N) def _get_xmid(self): return self.signal.xmid def _set_xmid(self, xmid): self.signal.xmid = xmid self.template.xmid = xmid = self.signal.Y xmid = property(_get_xmid, _set_xmid, doc="get/set the xmid value of the logistic signal") def _get_scale(self): return self.signal.scale def _set_scale(self, scale): self.signal.scale = scale self.template.scale = scale = xmid scale = property(_get_scale, _set_scale, doc="get/set the scale value of the logistic signal") def _get_X(self): return self.signal.X def _set_X(self, X): self.signal.X = X self.template.X = X = self.signal.Y X = property(_get_X, _set_X, doc="get/set the X values the logistic signal")
[docs] def set_noisy_data(self, sigma=0.1): = self.signal.Y + sigma * pylab.randn(self.N)[>1] = 1[<0] = 0
[docs] def get_snr(self, template, show=True): # first take fft of data and template # we double vector size to next power of 2 N = nextpower(self.N) data = np.concatenate([, [0] * (N-self.N+1)], axis=0) data_fft = np.fft.fft(data) #For this to work, we need the template and the data to be the # same length. So, we'll zero-pad the template before we take the FFT: template = np.concatenate([template, [0] * (N-self.N+1)], axis=0) template_fft = np.fft.fft(template) # we need an estimate of the noise power in each FFT bin. # we can assume white noise. # -- Calculate the PSD of the data #power_data, freq_psd = plt.psd(data[12*fs:], # Fs=fs, NFFT=fs, visible=False) # -- Interpolate to get the PSD values at the needed frequencies fs = 1 datafreq = np.fft.fftfreq(data.size)*fs #power_vec = np.interp(datafreq, freq_psd, power_data) power_vec = [1] * (N +1) ## 32 +f0 # we multiply the Fourier Space template and data, and divide # by the noise power in each frequency bin. Taking the Inverse # Fourier Transform (IFFT) of the filter output puts it back # in the time domain, so the result will be plotted as a function # of time off-set between the template and the data: # -- Calculate the matched filter output optimal = data_fft * template_fft.conjugate() / power_vec optimal_time = 2*np.fft.ifft(optimal) # Finally, we can normalize the matched filter output so that # we expect a value of 1 at times of just noise. Then, the peak # of the matched filter output will tell us the signal-to-noise # ratio (SNR) of the signal. # -- Normalize the matched filter output df = np.abs(datafreq[1] - datafreq[0]) sigmasq = 2*(template_fft * template_fft.conjugate() / power_vec).sum() * df sigma = np.sqrt(np.abs(sigmasq)) SNR = abs(optimal_time) / (sigma) # -- Plot the result if show:pylab.plot(SNR) return SNR
[docs] def scan(self, xmid_range, scale_range, method='mf', show=True): N = len(xmid_range) results = np.zeros((N, N)) coords = None best = 0 self.template.X = self.signal.X for i, xmid in enumerate(xmid_range): for j, scale in enumerate(scale_range): self.template.xmid = xmid self.template.scale = scale if method == 'corr': M = np.corrcoef(self.template.Y,[0,1] results[i, j] = M elif method == 'mf': M = self.get_snr(self.template.Y, show=False).max() results[i, j] = M if M > best: best = M coords = (xmid, scale) res = self.optimise(coords) Ex = abs(self.xmid-res.x[0])/self.xmid*100 Ey = abs(self.scale-res.x[1])/self.scale*100 if show is True: self.plot(results, coords, best, xmid_range, scale_range) pylab.figure(1) pylab.plot(self.xmid, self.scale, 'ok', markersize=10) pylab.plot(res.x[0], res.x[1], 'g', marker='D') pylab.title('Ex = {:.4} ; Ey={:.4}'.format(Ex, Ey)) results = {'results':results, 'coords':coords, 'best':best, 'Ex':Ex, 'Ey':Ey} return results
[docs] def plot(self, results, coords, best, xmid_range, scale_range): pylab.figure(1) pylab.clf() pylab.contourf(results.transpose(), 20, extent=[xmid_range[0], xmid_range[-1], scale_range[0], scale_range[-1]], aspect='auto') pylab.plot(coords[0], coords[1], 'ko') pylab.colorbar() pylab.figure(2) pylab.clf() pylab.plot(self.signal.X, self.signal.Y, '--k', label='input') pylab.plot(self.signal.X,, 'o', label='input') self.fitted = Logistic(coords[0], coords[1], N=self.N) self.fitted.X = self.X pylab.plot(self.fitted.X, self.fitted.Y, 'o-', label='input') pylab.legend(['signal','signal+noise','fitted']) pylab.grid()
[docs] def objective_function(self, xmid, scale): temp = Logistic(xmid, scale) temp.X = self.signal.X M = self.get_snr(temp.Y, show=False).max() return M
[docs] def optimise(self, guess=[1, 1]): from scipy.optimize import minimize def func(inputs): xmid = inputs[0] scale = inputs[1] M = self.objective_function(xmid, scale) return -M res = minimize(func, guess, method='nelder-mead', options={'xatol': 1e-8, 'disp': False}) return res