Repository: StingraySoftware/stingray
Pull request #181 (2017-04-04 09:49:03)

Fix docs of fit_powerspectrum


<!-- Reviewable:start -->
This change is [<img src="https://reviewable.io/review_button.svg" height="34" align="absmiddle" alt="Reviewable"/>](https://reviewable.io/reviews/stingraysoftware/stingray/181)
<!-- Reviewable:end -->



Patch:
diff --git a/.gitignore b/.gitignore
index 20eb3d05..684ae3ba 100644
--- a/.gitignore
+++ b/.gitignore
@@ -50,6 +50,3 @@ distribute-*.tar.gz
 
 # Mac OSX
 .DS_Store
-
-# virtualenv
-.venv
diff --git a/stingray/modeling/scripts.py b/stingray/modeling/scripts.py
index df91c9be..1659dd39 100644
--- a/stingray/modeling/scripts.py
+++ b/stingray/modeling/scripts.py
@@ -9,7 +9,7 @@ __all__ = ["fit_powerspectrum", "fit_lorentzians"]
 
 
 def fit_powerspectrum(ps, model, starting_pars, max_post=False, priors=None,
-              fitmethod="L-BFGS-B"):
+                      fitmethod="L-BFGS-B"):
     """
     Fit a number of Lorentzians to a power spectrum, possibly including white
     noise. Each Lorentzian has three parameters (amplitude, centroid position,
@@ -29,8 +29,9 @@ def fit_powerspectrum(ps, model, starting_pars, max_post=False, priors=None,
     ps : Powerspectrum
         A Powerspectrum object with the data to be fit
 
-    nlor : int
-        The number of Lorentzians to fit
+    model: astropy.modeling.models class instance
+        The parametric model supposed to represent the data. For details
+        see the astropy.modeling documentation
 
     starting_pars : iterable
         The list of starting guesses for the optimizer. See explanation above
diff --git a/stingray/simulator/models.py b/stingray/simulator/models.py
index 23b876d4..875200da 100644
--- a/stingray/simulator/models.py
+++ b/stingray/simulator/models.py
@@ -1,55 +1,4 @@
 from __future__ import division, print_function, absolute_import
-from astropy.modeling.models import custom_model
-
-@custom_model
-def GeneralizedLorentz1D(x, x_0=1., fwhm=1., value=1., power_coeff=1.):
-    """
-    Generalized lorenzian function, 
-    implemented using astropy.modeling.models custom model
-
-    Parameters
-    ----------
-
-    x: numpy.ndarray
-        non-zero frequencies
-
-    x_0 : peak centeral frequency
-    fwhm : FWHM of the peak (gamma)
-    value : peak value at x=x0
-    power_coeff : power coefficient [n]
-
-    Returns
-    -------
-    model: astropy.modeling.Model
-        generalized Lorentzian psd model
-    """
-    assert power_coeff > 0., "The power coefficient should be greater than zero."
-    return value * (fwhm / 2)**power_coeff * 1./(abs(x - x_0)**power_coeff + (fwhm / 2)**power_coeff) 
-    
-@custom_model
-def SmoothBrokenPowerLaw(x, norm=1., gamma_low=1., gamma_high=1., break_freq=1.):
-    """
-    Smooth broken power law function,
-    implemented using astropy.modeling.models custom model
-
-    Parameters
-    ----------
-
-    x: numpy.ndarray
-        non-zero frequencies
-
-    norm: normalization frequency
-    gamma_low: power law index for f --> zero
-    gamma_high: power law index for f --> infinity
-    break_freq: break frequency
-
-    Returns
-    -------
-    model: astropy.modeling.Model
-        generalized smooth broken power law psd model
-    """
-    return norm * x**(-gamma_low) / (1. + (x / break_freq)**2)**(-(gamma_low - gamma_high) / 2)
-    
 
 def lorenzian(x, p):
     """
diff --git a/stingray/simulator/simulator.py b/stingray/simulator/simulator.py
index fbe3de82..3264180f 100644
--- a/stingray/simulator/simulator.py
+++ b/stingray/simulator/simulator.py
@@ -2,14 +2,13 @@ from __future__ import division, print_function, absolute_import
 import numpy as np
 import numbers
 from scipy import signal
-import astropy.modeling.models
 
 from stingray import Lightcurve, AveragedPowerspectrum, io, utils
 import stingray.simulator.models as models
 
 class Simulator(object):
 
-    def __init__(self, dt=1, N=1024, mean=0, rms=1, red_noise=1, random_state=None):
+    def __init__(self, dt=1, N=1024, mean=0, rms=1, red_noise=1, seed=None):
         """
         Methods to simulate and visualize light curves.
 
@@ -25,7 +24,7 @@ class Simulator(object):
             fractional rms of the simulated light curve,
             actual rms is calculated by mean*rms
         red_noise: int, default 1
-            multiple of real length of light curve, by
+            multiple of real length of light curve, by 
             which to simulate, to avoid red noise leakage
         seed: int, default None
             seed value for random processes
@@ -40,8 +39,9 @@ class Simulator(object):
 
         # Initialize a tuple of energy ranges with corresponding light curves
         self.channels = []
-
-        self.random_state = utils.get_random_state(random_state)
+        
+        if seed is not None:
+            np.random.seed(seed)
 
         assert rms<=1, 'Fractional rms must be less than 1.'
         assert dt>0, 'Time resolution must be greater than 0'
@@ -53,18 +53,20 @@ class Simulator(object):
 
         Examples
         --------
-        - x = simulate(beta)
+        - x = simulate(2)
             For generating a light curve using power law spectrum.
 
             Parameters
             ----------
-            beta: int
+            Beta: int
                 Defines the shape of spectrum
+            N: int
+                Number of samples
 
             Returns
             -------
             lightCurve: `LightCurve` object
-
+        
         - x = simulate(s)
             For generating a light curve from user-provided spectrum.
 
@@ -77,27 +79,17 @@ class Simulator(object):
             -------
             lightCurve: `LightCurve` object
 
-        - x = simulate(model)
+        - x = simulate('lorenzian', p)
             For generating a light curve from pre-defined model
 
             Parameters
             ----------
-            model: astropy.modeling.Model
-                the pre-defined model
-
-            Returns
-            -------
-            lightCurve: `LightCurve` object
-
-        - x = simulate('model', params)
-            For generating a light curve from pre-defined model
+            model: str
+                name of the pre-defined model
 
-            Parameters
-            ----------
-            model: string
-                the pre-defined model
-            params: list iterable or dict
-                the parameters for the pre-defined model
+            p: list iterable
+                model parameters. For details, see model definitions
+                in model.py
 
             Returns
             -------
@@ -141,18 +133,15 @@ class Simulator(object):
             lightCurve: `LightCurve` object
         """
 
-        if isinstance(args[0], numbers.Integral) and len(args) == 1:
+        if isinstance(args[0], numbers.Integral) and len(args) == 1: 
             return  self._simulate_power_law(args[0])
 
-        elif isinstance(args[0], astropy.modeling.Model) and len(args) == 1:
-            return self._simulate_model(args[0])
-
-        elif utils.is_string(args[0]) and len(args) == 2:
-            return self._simulate_model_string(args[0], args[1])
-
         elif len(args) == 1:
             return self._simulate_power_spectrum(args[0])
 
+        elif utils.is_string(args[0]) and len(args) == 2:
+            return self._simulate_model(args[0], args[1])
+
         elif len(args) == 2:
             return self._simulate_impulse_response(args[0], args[1])
 
@@ -162,10 +151,9 @@ class Simulator(object):
         else:
             raise ValueError("Length of arguments must be 1, 2 or 3.")
 
-
     def simulate_channel(self, channel, *args):
         """
-        Simulate a lightcurve and add it to corresponding energy
+        Simulate a lightcurve and add it to corresponding energy 
         channel.
 
         Parameters
@@ -173,7 +161,7 @@ class Simulator(object):
         channel: str
             range of energy channel (e.g., 3.5-4.5)
 
-        *args:
+        *args: 
             see description of simulate() for details
 
         Returns
@@ -199,7 +187,7 @@ class Simulator(object):
         """
         Get multiple light curves belonging to the energy channels.
         """
-
+        
         return [lc[1] for lc in self.channels if lc[0] in channels]
 
     def get_all_channels(self):
@@ -215,7 +203,7 @@ class Simulator(object):
         """
 
         channel = [lc for lc in self.channels if lc[0] == channel]
-
+        
         if len(channel) == 0:
             raise KeyError('This channel does not exist or has already been deleted.')
         else:
@@ -241,12 +229,12 @@ class Simulator(object):
         """
         Return total number of energy channels.
         """
-
+        
         return len(self.channels)
 
     def simple_ir(self, start=0, width=1000, intensity=1):
         """
-        Construct a simple impulse response using start time,
+        Construct a simple impulse response using start time, 
         width and scaling intensity.
         To create a delta impulse response, set width to 1.
 
@@ -314,21 +302,21 @@ class Simulator(object):
         # Create a rising exponential of user-provided slope
         x = np.linspace(t1/dt, t2/dt, (t2-t1)/dt)
         h_rise = np.exp(rise*x)
-
+        
         # Evaluate a factor for scaling exponential
         factor = np.max(h_rise)/(p2-p1)
         h_secondary = (h_rise/factor) + p1
 
         # Create a decaying exponential until the end time
         x = np.linspace(t2/dt, t3/dt, (t3-t2)/dt)
-        h_decay = (np.exp((-decay)*(x-4/dt)))
+        h_decay = (np.exp((-decay)*(x-4/dt))) 
 
         # Add the three responses
         h = np.append(h_primary, h_secondary)
         h = np.append(h, h_decay)
 
         return h
-
+        
     def _simulate_power_law(self, B):
         """
         Generate LightCurve from a power law spectrum.
@@ -345,10 +333,10 @@ class Simulator(object):
 
         # Define frequencies at which to compute PSD
         w = np.fft.rfftfreq(self.red_noise*self.N, d=self.dt)[1:]
-
+        
         # Draw two set of 'N' guassian distributed numbers
-        a1 = self.random_state.normal(size=len(w))
-        a2 = self.random_state.normal(size=len(w))
+        a1 = np.random.normal(size=len(w))
+        a2 = np.random.normal(size=len(w))
 
         # Multiply by (1/w)^B to get real and imaginary parts
         real = a1 * np.power((1/w),B/2)
@@ -380,108 +368,55 @@ class Simulator(object):
         self.red_noise = 1
 
         # Draw two set of 'N' guassian distributed numbers
-        a1 = self.random_state.normal(size=len(s))
-        a2 = self.random_state.normal(size=len(s))
+        a1 = np.random.normal(size=len(s))
+        a2 = np.random.normal(size=len(s))
 
         lc = self._find_inverse(a1*s, a2*s)
         lc = Lightcurve(self.time, self._extract_and_scale(lc))
 
         return lc
 
-
-    def _simulate_model(self, model):
+    def _simulate_model(self, mod, p):
         """
-        For generating a light curve from a pre-defined model
+        For generating a light curve from pre-defined model
 
         Parameters
         ----------
-        model: astropy.modeling.Model derived function
-            the pre-defined model
-            (library-based, available in astropy.modeling.models or custom-defined)
-
-        Returns
-        -------
-        lightCurve: `LightCurve` object
-        """
-
-        # Frequencies at which the PSD is to be computed
-        # (only positive frequencies, since the signal is real)
-        nbins = self.red_noise*self.N
-        simfreq = np.fft.rfftfreq(nbins, d=self.dt)[1:]
-
-        # Compute PSD from model
-        simpsd = model(simfreq)
-
-        fac = np.sqrt(simpsd/2.)
-        pos_real   = self.random_state.normal(size=nbins//2)*fac
-        pos_imag   = self.random_state.normal(size=nbins//2)*fac
-
-        pos_freq_transform = pos_real + 1j * pos_imag
-
-        # Simulate light curve from its Fourier transform
-        arg  = np.concatenate(([self.mean], pos_freq_transform))
-
-        # Inverse Fourier transform
-        long_lc = np.fft.irfft(arg)
-
-        lc = Lightcurve(self.time, self._extract_and_scale(long_lc))
-        return lc
-
-
-    def _simulate_model_string(self, model_str, params):
-        """
-        For generating a light curve from a pre-defined model
-
-        Parameters
-        ----------
-        model_str: string
+        mod: str
             name of the pre-defined model
-        params: list or dictionary
-            parameters of the pre-defined model
+
+        p: list iterable
+            model parameters. For details, see model definitions
+            in model.py
 
         Returns
         -------
         lightCurve: `LightCurve` object
         """
 
-        # Frequencies at which the PSD is to be computed
-        # (only positive frequencies, since the signal is real)
-        nbins = self.red_noise*self.N
-        simfreq = np.fft.rfftfreq(nbins, d=self.dt)[1:]
-
-        if model_str in dir(models):
-            if isinstance(params, dict):
-                model = eval('models.' + model_str + '(**params)')
-                # Compute PSD from model
-                simpsd = model(simfreq)
-            elif isinstance(params, list):
-                simpsd = eval('models.' + model_str + '(simfreq, params)')
-            else:
-                raise ValueError('Params should be list or dictionary!')
-
-            fac = np.sqrt(simpsd/2.)
-            pos_real   = self.random_state.normal(size=nbins//2)*fac
-            pos_imag   = self.random_state.normal(size=nbins//2)*fac
-
-            pos_freq_transform = pos_real + 1j * pos_imag
+        # Define frequencies at which to compute PSD
+        w = np.fft.rfftfreq(self.red_noise*self.N, d=self.dt)[1:]
 
-            # Simulate light curve from its Fourier transform
-            arg  = np.concatenate(([self.mean], pos_freq_transform))
+        if mod in dir(models):
+            mod = 'models.'+ mod
+            s = eval(mod +'(w, p)')
 
-            # Inverse Fourier transform
-            long_lc = np.fft.irfft(arg)
+            # Draw two set of 'N' guassian distributed numbers
+            a1 = np.random.normal(size=len(s))
+            a2 = np.random.normal(size=len(s))
 
+            long_lc = self._find_inverse(a1*s, a2*s)
             lc = Lightcurve(self.time, self._extract_and_scale(long_lc))
+
             return lc
+        
         else:
             raise ValueError('Model is not defined!')
 
-
-
     def _simulate_impulse_response(self, s, h, mode='same'):
         """
         Generate LightCurve from impulse response. To get
-        accurate results, binning intervals (dt) of variability
+        accurate results, binning intervals (dt) of variability 
         signal 's' and impulse response 'h' must be equal.
 
         Parameters
@@ -509,10 +444,10 @@ class Simulator(object):
             lc = lc[:-(len(h) - 1)]
         elif mode == 'filtered':
             lc = lc[(len(h) - 1):-(len(h) - 1)]
-
+        
         time = self.dt * np.arange(len(lc))
         return Lightcurve(time, lc)
-
+        
     def _find_inverse(self, real, imaginary):
         """
         Forms complex numbers corresponding to real and imaginary
@@ -536,7 +471,7 @@ class Simulator(object):
         f = [complex(r, i) for r,i in zip(real,imaginary)]
 
         f = np.hstack([self.mean, f])
-
+        
         # Obtain real valued time series
         f_conj = np.conjugate(np.array(f))
 
@@ -566,8 +501,8 @@ class Simulator(object):
             lc = long_lc
         else:
             # Make random cut and extract light curve of length 'N'
-            extract = self.random_state.randint(self.N-1, self.red_noise*self.N - self.N+1)
-            lc = np.take(long_lc, range(extract, extract + self.N))
+            extract = np.random.randint(self.N-1, self.red_noise*self.N - self.N+1)
+            lc = np.take(long_lc, range(extract, extract + self.N)) 
 
         avg = np.mean(lc)
         std = np.std(lc)
@@ -588,7 +523,7 @@ class Simulator(object):
         -------
         power: numpy.ndarray
             The array of normalized squared absolute values of Fourier
-            amplitudes
+            amplitudes    
 
         """
         if seg_size is None:
@@ -615,7 +550,7 @@ class Simulator(object):
         """
 
         if format_ == 'pickle':
-            data = io.read(filename, 'pickle')
+            data = io.read(filename, 'pickle')        
             return data
         else:
             raise KeyError("Format not supported.")
diff --git a/stingray/utils.py b/stingray/utils.py
index b28118ec..52a19833 100644
--- a/stingray/utils.py
+++ b/stingray/utils.py
@@ -2,11 +2,9 @@ from __future__ import (absolute_import, unicode_literals, division,
                         print_function)
 import sys
 import collections
-import numbers
 
 import warnings
 import numpy as np
-
 # If numba is installed, import jit. Otherwise, define an empty decorator with
 # the same name.
 
@@ -200,19 +198,3 @@ def contiguous_regions(condition):
     # Reshape the result into two columns
     idx.shape = (-1, 2)
     return idx
-
-def is_int(obj):
-    return isinstance(obj, (numbers.Integral, np.integer))
-
-def get_random_state(random_state = None):
-    if not random_state:
-        random_state = np.random.mtrand._rand
-    else:
-        if is_int(random_state):
-            random_state = np.random.RandomState(random_state)
-        elif not isinstance(random_state, np.random.RandomState):
-            raise ValueError("{value} can't be used to generate a numpy.random.RandomState".format(
-                value = random_state
-            ))
-
-    return random_state


Test patch:
diff --git a/stingray/simulator/tests/test_simulator.py b/stingray/simulator/tests/test_simulator.py
index 78d807f0..b1d09fd1 100644
--- a/stingray/simulator/tests/test_simulator.py
+++ b/stingray/simulator/tests/test_simulator.py
@@ -2,7 +2,6 @@ import numpy as np
 import os
 
 from astropy.tests.helper import pytest
-import astropy.modeling.models
 from stingray import Lightcurve, Crossspectrum, sampledata
 from stingray.simulator import simulator, models
 
@@ -40,12 +39,9 @@ class TestSimulator(object):
         """
         Simulate with a random seed value.
         """
-        self.simulator = simulator.Simulator(N=1024, random_state=12)
+        self.simulator = simulator.Simulator(N=1024, seed=12)
         assert len(self.simulator.simulate(2).counts), 1024
 
-    def test_simulate_with_random_state(self):
-        self.simulator = simulator.Simulator(N=1024, random_state=np.random.RandomState(12))
-
     def test_simulate_with_incorrect_arguments(self):
         with pytest.raises(ValueError):
             self.simulator.simulate(1, 2, 3, 4)
@@ -206,74 +202,6 @@ class TestSimulator(object):
         assert np.all(
             np.abs(actual_prob - simulated_prob) < 3*np.sqrt(actual_prob))
 
-    def test_simulate_GeneralizedLorentz1D_str(self):
-        """
-        Simulate a light curve using the GeneralizedLorentz1D model
-        called as a string
-        """
-        assert len(self.simulator.simulate('GeneralizedLorentz1D', {'x_0':10, 'fwhm':1., 'value':10., 'power_coeff':2})), 1024
-        
-    def test_simulate_GeneralizedLorentz1D(self):
-        """
-        Simulate a light curve using the GeneralizedLorentz1D model
-        called as a astropy.modeling.Model class
-        """
-        mod = models.GeneralizedLorentz1D(x_0=10, fwhm=1., value=10., power_coeff=2)
-        assert len(self.simulator.simulate(mod)), 1024
-        
-    def test_simulate_SmoothBrokenPowerLaw_str(self):
-        """
-        Simulate a light curve using SmoothBrokenPowerLaw model
-        called as a string
-        """
-        assert len(self.simulator.simulate('SmoothBrokenPowerLaw', {'norm':1., 'gamma_low':1., 'gamma_high':2., 'break_freq':1.})), 1024
-    
-    def test_simulate_SmoothBrokenPowerLaw(self):
-        """
-        Simulate a light curve using SmoothBrokenPowerLaw model
-        called as a astropy.modeling.Model class
-        """
-        mod = models.SmoothBrokenPowerLaw(norm=1., gamma_low=1., gamma_high=2., break_freq=1.)
-        assert len(self.simulator.simulate(mod)), 1024
-        
-        
-    def test_simulate_generic_model(self):
-        """
-        Simulate a light curve using a generic model
-        called as a astropy.modeling.Model class
-        """
-        mod = astropy.modeling.models.Gaussian1D(amplitude=10., mean=1., stddev=2.)
-        assert len(self.simulator.simulate(mod)), 1024
-        
-        
-    def test_compare_composite(self):
-        """
-        Compare the PSD of a light curve simulated using a composite model
-        (using SmoothBrokenPowerLaw plus GeneralizedLorentz1D)
-        with the actual model
-        """
-        N = 50000
-        dt = 0.01
-        m = 30000.
-        
-        self.simulator = simulator.Simulator(N=N, mean=m, dt=dt)
-        smoothbknpo = models.SmoothBrokenPowerLaw(norm=1., gamma_low=1., gamma_high=2., break_freq=1.)
-        lorentzian = models.GeneralizedLorentz1D(x_0=10, fwhm=1., value=10., power_coeff=2.)
-        myModel = smoothbknpo + lorentzian
-        
-        lc = [self.simulator.simulate(myModel) for i in range(1, 50)]
-
-        simulated = self.simulator.powerspectrum(lc, lc[0].tseg)
-
-        w = np.fft.rfftfreq(N, d=dt)[1:]
-        actual = myModel(w)[:-1]
-
-        actual_prob = actual/float(sum(actual))
-        simulated_prob = simulated/float(sum(simulated))
-
-        assert np.all(np.abs(actual_prob - simulated_prob) < 3*np.sqrt(actual_prob))
-        
-        
     def test_simulate_wrong_model(self):
         """
         Simulate with a model that does not exist.
diff --git a/stingray/tests/test_utils.py b/stingray/tests/test_utils.py
index 38e918bb..24b14e02 100644
--- a/stingray/tests/test_utils.py
+++ b/stingray/tests/test_utils.py
@@ -95,15 +95,3 @@ class TestUtils(object):
         cont = utils.contiguous_regions(array)
         assert np.all(cont == np.array([[1, 3], [4, 7]])), \
             'Contiguous region wrong'
-
-    def test_get_random_state(self):
-        # Life, Universe and Everything
-        lue = 42
-        random_state = np.random.RandomState(lue)
-
-        assert utils.get_random_state(None) is np.random.mtrand._rand
-        assert np.all(utils.get_random_state(lue).randn(lue) == np.random.RandomState(lue).randn(lue))
-        assert np.all(utils.get_random_state(np.random.RandomState(lue)).randn(lue) == np.random.RandomState(lue).randn(lue))
-
-        with pytest.raises(ValueError):
-            utils.get_random_state('foobar')