Repository: StingraySoftware/stingray
Pull request #217 (2017-07-23 08:33:36)

Fix the simulator for an odd number of bins
Fixes issue reported by @ricardovb in DAVE


Patch:
diff --git a/stingray/simulator/simulator.py b/stingray/simulator/simulator.py
index fb67e086..72b5788c 100644
--- a/stingray/simulator/simulator.py
+++ b/stingray/simulator/simulator.py
@@ -419,13 +419,7 @@ class Simulator(object):
         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)
+        long_lc = self._find_inverse(pos_real, pos_imag)
 
         lc = Lightcurve(self.time, self._extract_and_scale(long_lc),
                         err_dist='gauss', dt=self.dt)
@@ -467,13 +461,7 @@ class Simulator(object):
             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)
+            long_lc = self._find_inverse(pos_real, pos_imag)
 
             lc = Lightcurve(self.time, self._extract_and_scale(long_lc),
                             err_dist='gauss', dt=self.dt)
@@ -481,8 +469,6 @@ class Simulator(object):
         else:
             raise ValueError('Model is not defined!')
 
-
-
     def _simulate_impulse_response(self, s, h, mode='same'):
         """
         Generate LightCurve from impulse response. To get


Test patch:
diff --git a/stingray/simulator/tests/test_simulator.py b/stingray/simulator/tests/test_simulator.py
index 78d807f0..5130f3e1 100644
--- a/stingray/simulator/tests/test_simulator.py
+++ b/stingray/simulator/tests/test_simulator.py
@@ -19,6 +19,7 @@ class TestSimulator(object):
     @classmethod
     def setup_class(self):
         self.simulator = simulator.Simulator(N=1024, mean=0.5, dt=0.125)
+        self.simulator_odd = simulator.Simulator(N=2039, mean=0.5, dt=0.125)
 
     def calculate_lag(self, lc, h, delay):
         """
@@ -57,6 +58,13 @@ class TestSimulator(object):
         self.simulator.simulate_channel('3.5-4.5', 'lorenzian', [1, 2, 3, 4])
         self.simulator.delete_channel('3.5-4.5')
 
+    def test_simulate_channel_odd(self):
+        """
+        Simulate an energy channel.
+        """
+        self.simulator_odd.simulate_channel('3.5-4.5', 'lorenzian', [1, 2, 3, 4])
+        self.simulator_odd.delete_channel('3.5-4.5')
+
     def test_incorrect_simulate_channel(self):
         """Test simulating a channel that already exists."""
         self.simulator.simulate_channel('3.5-4.5', 2)
@@ -122,6 +130,12 @@ class TestSimulator(object):
         """
         assert len(self.simulator.simulate(2).counts), 1024
 
+    def test_simulate_powerlaw_odd(self):
+        """
+        Simulate light curve from power law spectrum.
+        """
+        assert len(self.simulator_odd.simulate(2).counts), 2039
+
     def test_compare_powerlaw(self):
         """
         Compare simulated power spectrum with actual one.
@@ -156,6 +170,13 @@ class TestSimulator(object):
         """
         assert len(self.simulator.simulate('lorenzian', [1, 2, 3, 4])), 1024
 
+    def test_simulate_lorenzian_odd(self):
+        """
+        Simulate light curve using lorenzian model.
+        """
+        assert len(self.simulator_odd.simulate('lorenzian',
+                                               [1, 2, 3, 4])), 1024
+
     def test_compare_lorenzian(self):
         """
         Compare simulated lorenzian spectrum with original spectrum.
@@ -211,8 +232,20 @@ class TestSimulator(object):
         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
+        assert len(self.simulator.simulate('GeneralizedLorentz1D',
+                                           {'x_0':10, 'fwhm':1., 'value':10.,
+                                            'power_coeff':2})), 1024
         
+    def test_simulate_GeneralizedLorentz1D_odd_str(self):
+        """
+        Simulate a light curve using the GeneralizedLorentz1D model
+        called as a string
+        """
+        assert len(self.simulator_odd.simulate('GeneralizedLorentz1D',
+                                               {'x_0':10, 'fwhm':1.,
+                                                'value':10., 'power_coeff':2}
+                                               )), 2039
+
     def test_simulate_GeneralizedLorentz1D(self):
         """
         Simulate a light curve using the GeneralizedLorentz1D model
@@ -244,8 +277,16 @@ class TestSimulator(object):
         """
         mod = astropy.modeling.models.Gaussian1D(amplitude=10., mean=1., stddev=2.)
         assert len(self.simulator.simulate(mod)), 1024
-        
-        
+
+    def test_simulate_generic_model_odd(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_odd.simulate(mod)), 2039
+
     def test_compare_composite(self):
         """
         Compare the PSD of a light curve simulated using a composite model
@@ -286,15 +327,34 @@ class TestSimulator(object):
         Construct simple impulse response.
         """
         t0, w = 100, 500
-        assert len(self.simulator.simple_ir(t0, w)), (t0+w)/self.simulator.dt
+        assert len(self.simulator.simple_ir(t0, w)) == \
+               (t0+w)/self.simulator.dt
+
+    def test_construct_simple_ir_odd(self):
+        """
+        Construct simple impulse response.
+        """
+        t0, w = 100, 500
+        assert len(self.simulator_odd.simple_ir(t0, w)) == \
+               (t0+w)/self.simulator.dt
 
     def test_construct_relativistic_ir(self):
         """
         Construct relativistic impulse response.
         """
         t1, t3 = 3, 10
-        assert len(self.simulator.relativistic_ir(t1=t1, t3=t3)),\
-                  (t1+t3)/self.simulator.dt
+        ir = self.simulator.relativistic_ir(t1=t1, t3=t3)
+        assert np.all(ir[:int(t1 / self.simulator.dt)] == 0)
+        assert ir[int(t1 / self.simulator.dt)] == 1
+
+    def test_construct_relativistic_ir_odd(self):
+        """
+        Construct relativistic impulse response.
+        """
+        t1, t3 = 3, 10
+        ir = self.simulator_odd.relativistic_ir(t1=t1, t3=t3)
+        assert np.all(ir[:int(t1 / self.simulator_odd.dt)] == 0)
+        assert ir[int(t1 / self.simulator_odd.dt)] == 1
 
     def test_simulate_simple_impulse(self):
         """
@@ -305,12 +365,28 @@ class TestSimulator(object):
         h = self.simulator.simple_ir(10, 1, 1)
         output = self.simulator.simulate(s, h)
 
+    def test_simulate_simple_impulse_odd(self):
+        """
+        Simulate light curve from simple impulse response.
+        """
+        lc = sampledata.sample_data()
+        s = lc.counts
+        h = self.simulator_odd.simple_ir(10, 1, 1)
+        output = self.simulator_odd.simulate(s, h)
+
     def test_powerspectrum(self):
         """
         Create a power spectrum from light curve.
         """
-        self.simulator.simulate(2)
-        self.simulator.powerspectrum(self.simulator.lc)
+        lc = self.simulator.simulate(2)
+        self.simulator.powerspectrum(lc)
+
+    def test_powerspectrum_odd(self):
+        """
+        Create a power spectrum from light curve.
+        """
+        lc = self.simulator_odd.simulate(2)
+        self.simulator_odd.powerspectrum(lc)
 
     def test_simulate_relativistic_impulse(self):
         """
@@ -332,6 +408,16 @@ class TestSimulator(object):
         h = self.simulator.simple_ir()
         output = self.simulator.simulate(s, h, 'filtered')
 
+    def test_filtered_simulate_odd(self):
+        """
+        Simulate light curve using 'filtered' mode.
+        """
+        lc = sampledata.sample_data()
+        s = lc.counts
+
+        h = self.simulator_odd.simple_ir()
+        output = self.simulator_odd.simulate(s, h, 'filtered')
+
     def test_simple_lag_spectrum(self):
         """
         Simulate light curve from simple impulse response and
@@ -419,13 +505,6 @@ class TestSimulator(object):
         assert np.abs(5-h_cutoffs[0]) < np.sqrt(5)
         assert np.abs(5-h_cutoffs[1]) < np.sqrt(5)
 
-    def test_powerspectrum(self):
-        """
-        Create a powerspectrum from light curve.
-        """
-        lc = self.simulator.simulate(2)
-        self.simulator.powerspectrum(lc)
-
     def test_io(self):
         sim = simulator.Simulator(N=1024)
         sim.write('sim.pickle')