test/test_Signal.py

   1 #!/usr/bin/env python
   2 # Time-stamp: <2019-12-18 17:02:57 taoliu>
   3
   4 """Module Description: Test functions for Signal.pyx
   5
   6 This code is free software; you can redistribute it and/or modify it
   7 under the terms of the BSD License (see the file LICENSE included with
   8 the distribution).
   9 """
  10
  11 import unittest
  12 import pytest
  13
  14 from math import log10
  15 import numpy as np
  16 from MACS2.Signal import maxima, savitzky_golay, savitzky_golay_order2_deriv1
  17
  18 # ------------------------------------
  19 # Main function
  20 # ------------------------------------
  21
  22 class Test_maxima(unittest.TestCase):
  23
  24     def setUp(self):
  25         self.signal = np.array( [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 4, 4, 4,
  26                                   4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5,
  27                                   5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
  28                                   5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
  29                                   6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7,
  30                                   7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
  31                                   7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7,
  32                                   7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
  33                                   7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8,
  34                                   8, 8, 8, 8, 8, 8, 7, 7, 7, 7, 7, 7, 7, 7, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
  35                                   6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 4, 4, 4, 4, 4,
  36                                   4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
  37                                   4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
  38                                   4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
  39                                   3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
  40                                   3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
  41                                   3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
  42                                   3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
  43                                   3, 3, 3, 3, 3, 3 ], dtype = "float32" )
  44         self.windowsize = 253
  45         self.summit = 161       # this is based on 1-deriv smoothed data
  46         self.smoothed162 = -2.98155597e-18 # the value is based on python3.7+numpy1.17.4; from smoothing func with np.convolve
  47
  48     # this test uses only numpy functions. we aim to capture strange behavior under specific py+numpy
  49     @pytest.mark.skip( reason="it fails under some combinations of py+np with unknown reason." )
  50     def test_implement_smooth_here ( self ):
  51         signal = self.signal
  52         window_size = self.windowsize
  53         half_window = (window_size - 1) // 2
  54         # precompute coefficients
  55         b = np.array([[1, k, k**2] for k in range(-half_window, half_window+1)], dtype='int64')
  56         m = np.linalg.pinv(b)[1]
  57         # pad the signal at the extremes with
  58         # values taken from the signal itself
  59         firstvals = signal[0] - np.abs(signal[1:half_window+1][::-1] - signal[0])
  60         lastvals = signal[-1] + np.abs(signal[-half_window-1:-1][::-1] - signal[-1])
  61         signal = np.concatenate((firstvals, signal, lastvals))
  62         ret = np.convolve( m[::-1], signal.astype("float64"), mode='valid').astype("float32") # convolve function seems to have odd ret with spec py+numpy
  63         p = ret[162]
  64         print ("calculated step by step:\n", p)
  65         print ("expected:\n", self.smoothed162)
  66         self.assertAlmostEqual( p, self.smoothed162, places = 16 )
  67         self.assertEqual( np.sign(p), np.sign(self.smoothed162) )
  68
  69     def test_maxima(self):
  70         expect = self.summit
  71         result = maxima( self.signal, self.windowsize )[0]
  72         self.assertEqual( result, expect, msg=f"Not equal: result: {result}, expected: {expect}" )
  73
  74     def assertEqual_nparray1d ( self, a, b, places = 7 ):
  75         self.assertEqual( a.shape[0], b.shape[0] )
  76         l = a.shape[0]
  77         for i in range( l ):
  78             self.assertAlmostEqual( a[i], b[i], places = places, msg=f"Not equal at {i} {a[i]} {b[i]}" )
  79