Improve sinefit resiliancy

And adjust the period locator accordingly.
Fitting \sin is the same mathematically, it's just conceptually more
straightforward since we're locating zero crossings anyway.

README.txt

@@ -5,10 +5,10 @@ by Jim Paris <jim@jtan.com>
 Prerequisites:
 
   # Runtime and build environments
-  sudo apt-get install python2.7 python2.7-dev python-setuptools
+  sudo apt-get install python2.7 python2.7-dev python-setuptools python-pip
-  sudo apt-get install python-numpy python-scipy python-matplotlib
+  sudo apt-get install python-numpy python-scipy
 
-  nilmdb (1.5.0+)
+  nilmdb (1.6.3+)
 
 Install:
 

setup.py

@@ -61,10 +61,10 @@ setup(name='nilmtools',
       long_description = "NILM Database Tools",
       license = "Proprietary",
       author_email = 'jim@jtan.com',
-      install_requires = [ 'nilmdb >= 1.5.0',
+      install_requires = [ 'nilmdb >= 1.6.3',
                            'numpy',
                            'scipy',
-                           'matplotlib',
+                           #'matplotlib',
                            ],
       packages = [ 'nilmtools',
                    ],
@@ -79,6 +79,7 @@ setup(name='nilmtools',
               'nilm-copy-wildcard = nilmtools.copy_wildcard:main',
               'nilm-sinefit = nilmtools.sinefit:main',
               'nilm-cleanup = nilmtools.cleanup:main',
+              'nilm-median = nilmtools.median:main',
               ],
           },
       zip_safe = False,

src/cleanup.py

@@ -238,12 +238,15 @@ def main(argv = None):
                        timestamp_to_seconds(total)))
             continue
         printf("  removing data before %s\n", timestamp_to_human(remove_before))
-        if args.yes:
+        # Clean in reverse order.  Since we only use the primary stream and not
-            client.stream_remove(path, None, remove_before)
+        # the decimated streams to figure out which data to remove, removing
-        for ap in streams[path].also_clean_paths:
+        # the primary stream last means that we might recover more nicely if
-            printf("  also removing from %s\n", ap)
+        # we are interrupted and restarted.
+        clean_paths = list(reversed(streams[path].also_clean_paths)) + [ path ]
+        for p in clean_paths:
+            printf("  removing from %s\n", p)
             if args.yes:
-                client.stream_remove(ap, None, remove_before)
+                client.stream_remove(p, None, remove_before)
 
     # All done
     if not args.yes:

src/decimate_auto.py

@@ -4,15 +4,19 @@ import nilmtools.filter
 import nilmtools.decimate
 import nilmdb.client
 import argparse
+import fnmatch
 
 def main(argv = None):
     parser = argparse.ArgumentParser(
         formatter_class = argparse.RawDescriptionHelpFormatter,
-        version = "1.0",
+        version = nilmtools.__version__,
         description = """\
     Automatically create multiple decimations from a single source
     stream, continuing until the last decimated level contains fewer
     than 500 points total.
+
+    Wildcards and multiple paths are accepted.  Decimated paths are
+    ignored when matching wildcards.
     """)
     parser.add_argument("-u", "--url", action="store",
                         default="http://localhost/nilmdb/",
@@ -23,20 +27,36 @@ def main(argv = None):
                         default = False,
                         help="Force metadata changes if the dest "
                         "doesn't match")
-    parser.add_argument("path", action="store",
+    parser.add_argument("path", action="store", nargs='+',
                         help='Path of base stream')
     args = parser.parse_args(argv)
 
     # Pull out info about the base stream
     client = nilmdb.client.Client(args.url)
 
-    info = nilmtools.filter.get_stream_info(client, args.path)
+    # Find list of paths to process
-    if not info:
+    streams = [ unicode(s[0]) for s in client.stream_list() ]
-        raise Exception("path " + args.path + " not found")
+    streams = [ s for s in streams if "~decim-" not in s ]
+    paths = []
+    for path in args.path:
+        new = fnmatch.filter(streams, unicode(path))
+        if not new:
+            print "error: no stream matched path:", path
+            raise SystemExit(1)
+        paths.extend(new)
 
-    meta = client.stream_get_metadata(args.path)
+    for path in paths:
+        do_decimation(client, args, path)
+
+def do_decimation(client, args, path):
+    print "Decimating", path
+    info = nilmtools.filter.get_stream_info(client, path)
+    if not info:
+        raise Exception("path " + path + " not found")
+
+    meta = client.stream_get_metadata(path)
     if "decimate_source" in meta:
-        print "Stream", args.path, "was decimated from", meta["decimate_source"]
+        print "Stream", path, "was decimated from", meta["decimate_source"]
         print "You need to pass the base stream instead"
         raise SystemExit(1)
 
@@ -53,7 +73,7 @@ def main(argv = None):
         if info.rows <= 500:
             break
         factor *= args.factor
-        new_path = "%s~decim-%d" % (args.path, factor)
+        new_path = "%s~decim-%d" % (path, factor)
 
         # Create the stream if needed
         new_info = nilmtools.filter.get_stream_info(client, new_path)
@@ -72,5 +92,7 @@ def main(argv = None):
         # Update info using the newly decimated stream
         info = nilmtools.filter.get_stream_info(client, new_path)
 
+    return
+
 if __name__ == "__main__":
     main()

src/filter.py

@@ -67,7 +67,7 @@ def get_stream_info(client, path):
 
 class Filter(object):
 
-    def __init__(self):
+    def __init__(self, parser_description = None):
         self._parser = None
         self._client_src = None
         self._client_dest = None
@@ -78,6 +78,9 @@ class Filter(object):
         self.end = None
         self.interhost = False
         self.force_metadata = False
+        if parser_description is not None:
+            self.setup_parser(parser_description)
+            self.parse_args()
 
     @property
     def client_src(self):
@@ -233,8 +236,14 @@ class Filter(object):
         metadata = self._client_dest.stream_get_metadata(self.dest.path)
         if not self.force_metadata:
             for key in data:
-                wanted = str(data[key])
+                wanted = data[key]
+                if not isinstance(wanted, basestring):
+                    wanted = str(wanted)
                 val = metadata.get(key, wanted)
+                # Force UTF-8 encoding for comparison and display
+                wanted = wanted.encode('utf-8')
+                val = val.encode('utf-8')
+                key = key.encode('utf-8')
                 if val != wanted and self.dest.rows > 0:
                     m =  "Metadata in destination stream:\n"
                     m += "  %s = %s\n" % (key, val)
@@ -275,6 +284,10 @@ class Filter(object):
         Return value of 'function' is the number of data rows processed.
         Unprocessed data will be provided again in a subsequent call
         (unless 'final' is True).
+
+        If unprocessed data remains after 'final' is True, the interval
+        being inserted will be ended at the timestamp of the first
+        unprocessed data point.
         """
         if args is None:
             args = []
@@ -319,7 +332,13 @@ class Filter(object):
 
                 # Last call for this contiguous interval
                 if old_array.shape[0] != 0:
-                    function(old_array, interval, args, insert_function, True)
+                    processed = function(old_array, interval, args,
+                                         insert_function, True)
+                    if processed != old_array.shape[0]:
+                        # Truncate the interval we're inserting at the first
+                        # unprocessed data point.  This ensures that
+                        # we'll not miss any data when we run again later.
+                        insert_ctx.update_end(old_array[processed][0])
 
 def main(argv = None):
     # This is just a dummy function; actual filters can use the other

src/median.py

@@ -0,0 +1,43 @@
+#!/usr/bin/python
+import nilmtools.filter, scipy.signal
+
+def main(argv = None):
+    f = nilmtools.filter.Filter()
+    parser = f.setup_parser("Median Filter")
+    group = parser.add_argument_group("Median filter options")
+    group.add_argument("-z", "--size", action="store", type=int, default=25,
+                       help = "median filter size (default %(default)s)")
+    group.add_argument("-d", "--difference", action="store_true",
+                       help = "store difference rather than filtered values")
+
+    try:
+        args = f.parse_args(argv)
+    except nilmtools.filter.MissingDestination as e:
+        print "Source is %s (%s)" % (e.src.path, e.src.layout)
+        print "Destination %s doesn't exist" % (e.dest.path)
+        print "You could make it with a command like:"
+        print "  nilmtool -u %s create %s %s" % (e.dest.url,
+                                                 e.dest.path, e.src.layout)
+        raise SystemExit(1)
+
+    meta = f.client_src.stream_get_metadata(f.src.path)
+    f.check_dest_metadata({ "median_filter_source": f.src.path,
+                            "median_filter_size": args.size,
+                            "median_filter_difference": repr(args.difference) })
+
+    f.process_numpy(median_filter, args = (args.size, args.difference))
+
+def median_filter(data, interval, args, insert, final):
+    (size, diff) = args
+    (rows, cols) = data.shape
+    for i in range(cols - 1):
+        filtered = scipy.signal.medfilt(data[:, i+1], size)
+        if diff:
+            data[:, i+1] -= filtered
+        else:
+            data[:, i+1] = filtered
+    insert(data)
+    return rows
+
+if __name__ == "__main__":
+    main()

src/prep.py

@@ -3,6 +3,8 @@
 # Spectral envelope preprocessor.
 # Requires two streams as input: the original raw data, and sinefit data.
 
+from nilmdb.utils.printf import *
+from nilmdb.utils.time import timestamp_to_human
 import nilmtools.filter
 import nilmdb.client
 from numpy import *
@@ -77,7 +79,8 @@ def main(argv = None):
     # Check and set metadata in prep stream
     f.check_dest_metadata({ "prep_raw_source": f.src.path,
                             "prep_sinefit_source": sinefit.path,
-                            "prep_column": args.column })
+                            "prep_column": args.column,
+                            "prep_rotation": repr(rotation) })
 
     # Run the processing function on all data
     f.process_numpy(process, args = (client_sinefit, sinefit.path, args.column,
@@ -105,7 +108,6 @@ def process(data, interval, args, insert_function, final):
     # Pull out sinefit data for the entire time range of this block
     for sinefit_line in client.stream_extract(sinefit_path,
                                               data[0, 0], data[rows-1, 0]):
-
         def prep_period(t_min, t_max, rot):
             """
             Compute prep coefficients from time t_min to t_max, which
@@ -162,7 +164,15 @@ def process(data, interval, args, insert_function, final):
                 break
             processed = idx_max
 
-    print "Processed", processed, "of", rows, "rows"
+    # If we processed no data but there's lots in here, pretend we
+    # processed half of it.
+    if processed == 0 and rows > 10000:
+        processed = rows / 2
+        printf("%s: warning: no periods found; skipping %d rows\n",
+               timestamp_to_human(data[0][0]), processed)
+    else:
+        printf("%s: processed %d of %d rows\n",
+               timestamp_to_human(data[0][0]), processed, rows)
     return processed
 
 if __name__ == "__main__":

src/sinefit.py

@@ -2,12 +2,18 @@
 
 # Sine wave fitting.  This runs about 5x faster than realtime on raw data.
 
+from nilmdb.utils.printf import *
 import nilmtools.filter
 import nilmdb.client
+from nilmdb.utils.time import (timestamp_to_human,
+                               timestamp_to_seconds,
+                               seconds_to_timestamp)
+
 from numpy import *
 from scipy import *
 #import pylab as p
 import operator
+import sys
 
 def main(argv = None):
     f = nilmtools.filter.Filter()
@@ -18,6 +24,15 @@ def main(argv = None):
     group.add_argument('-f', '--frequency', action='store', type=float,
                        default=60.0,
                        help='Approximate frequency (default: %(default)s)')
+    group.add_argument('-m', '--min-freq', action='store', type=float,
+                       help='Minimum valid frequency '
+                       '(default: approximate frequency / 2))')
+    group.add_argument('-M', '--max-freq', action='store', type=float,
+                       help='Maximum valid frequency '
+                       '(default: approximate frequency * 2))')
+    group.add_argument('-a', '--min-amp', action='store', type=float,
+                       default=20.0,
+                       help='Minimum signal amplitude (default: %(default)s)')
 
     # Parse arguments
     try:
@@ -34,17 +49,56 @@ def main(argv = None):
         parser.error("need a column number >= 1")
     if args.frequency < 0.1:
         parser.error("frequency must be >= 0.1")
+    if args.min_freq is None:
+        args.min_freq = args.frequency / 2
+    if args.max_freq is None:
+        args.max_freq = args.frequency * 2
+    if (args.min_freq > args.max_freq or
+        args.min_freq > args.frequency or
+        args.max_freq < args.frequency):
+        parser.error("invalid min or max frequency")
+    if args.min_amp < 0:
+        parser.error("min amplitude must be >= 0")
 
     f.check_dest_metadata({ "sinefit_source": f.src.path,
                             "sinefit_column": args.column })
-    f.process_numpy(process, args = (args.column, args.frequency))
+    f.process_numpy(process, args = (args.column, args.frequency, args.min_amp,
+                                     args.min_freq, args.max_freq))
+
+class SuppressibleWarning(object):
+    def __init__(self, maxcount = 10, maxsuppress = 100):
+        self.maxcount = maxcount
+        self.maxsuppress = maxsuppress
+        self.count = 0
+        self.last_msg = ""
+
+    def _write(self, sec, msg):
+        if sec:
+            now = "[" + timestamp_to_human(seconds_to_timestamp(sec)) + "] "
+        else:
+            now = ""
+        sys.stderr.write(now + msg)
+
+    def warn(self, msg, seconds = None):
+        self.count += 1
+        if self.count <= self.maxcount:
+            self._write(seconds, msg)
+        if (self.count - self.maxcount) >= self.maxsuppress:
+            self.reset(seconds)
+
+    def reset(self, seconds = None):
+        if self.count > self.maxcount:
+            self._write(seconds, sprintf("(%d warnings suppressed)\n",
+                                         self.count - self.maxcount))
+        self.count = 0
 
 def process(data, interval, args, insert_function, final):
-    (column, f_expected) = args
+    (column, f_expected, a_min, f_min, f_max) = args
     rows = data.shape[0]
 
     # Estimate sampling frequency from timestamps
-    fs = 1e6 * (rows-1) / (data[-1][0] - data[0][0])
+    fs = (rows-1) / (timestamp_to_seconds(data[-1][0]) -
+                     timestamp_to_seconds(data[0][0]))
 
     # Pull out about 3.5 periods of data at once;
     # we'll expect to match 3 zero crossings in each window
@@ -54,30 +108,41 @@ def process(data, interval, args, insert_function, final):
     if rows < N:
         return 0
 
+    warn = SuppressibleWarning(3, 1000)
+
     # Process overlapping windows
     start = 0
     num_zc = 0
+    last_inserted_timestamp = None
     while start < (rows - N):
         this = data[start:start+N, column]
-        t_min = data[start, 0]/1e6
+        t_min = timestamp_to_seconds(data[start, 0])
-        t_max = data[start+N-1, 0]/1e6
+        t_max = timestamp_to_seconds(data[start+N-1, 0])
 
         # Do 4-parameter sine wave fit
         (A, f0, phi, C) = sfit4(this, fs)
 
         # Check bounds.  If frequency is too crazy, ignore this window
-        if f0 < (f_expected/2) or f0 > (f_expected*2):
+        if f0 < f_min or f0 > f_max:
-            print "frequency", f0, "too far from expected value", f_expected
+            warn.warn(sprintf("frequency %s outside valid range %s - %s\n",
+                              str(f0), str(f_min), str(f_max)), t_min)
+            start += N
+            continue
+
+        # If amplitude is too low, results are probably just noise
+        if A < a_min:
+            warn.warn(sprintf("amplitude %s below minimum threshold %s\n",
+                              str(A), str(a_min)), t_min)
             start += N
             continue
 
         #p.plot(arange(N), this)
-        #p.plot(arange(N), A * cos(f0/fs * 2 * pi * arange(N) + phi) + C, 'g')
+        #p.plot(arange(N), A * sin(f0/fs * 2 * pi * arange(N) + phi) + C, 'g')
 
-        # Period starts when the argument of cosine is 3*pi/2 degrees,
+        # Period starts when the argument of sine is 0 degrees,
         # so we're looking for sample number:
-        #     n = (3 * pi / 2 - phi) / (f0/fs * 2 * pi)
+        #     n = (0 - phi) / (f0/fs * 2 * pi)
-        zc_n = (3 * pi / 2 - phi) / (f0 / fs * 2 * pi)
+        zc_n = (0 - phi) / (f0 / fs * 2 * pi)
         period_n = fs/f0
 
         # Add periods to make N positive
@@ -90,7 +155,13 @@ def process(data, interval, args, insert_function, final):
         while zc_n < (N - period_n/2):
             #p.plot(zc_n, C, 'ro')
             t = t_min + zc_n / fs
-            insert_function([[t * 1e6, f0, A, C]])
+            if (last_inserted_timestamp is None or
+                t > last_inserted_timestamp):
+                insert_function([[seconds_to_timestamp(t), f0, A, C]])
+                last_inserted_timestamp = t
+                warn.reset(t)
+            else:
+                warn.warn("timestamp overlap\n", t)
             num_zc += 1
             last_zc = zc_n
             zc_n += period_n
@@ -108,6 +179,7 @@ def process(data, interval, args, insert_function, final):
         start = int(round(start + advance))
 
     # Return the number of rows we've processed
+    warn.reset(last_inserted_timestamp)
     print "Marked", num_zc, "zero-crossings in", start, "rows"
     return start
 
@@ -123,15 +195,15 @@ def sfit4(data, fs):
 
     Output:
       Parameters [A, f0,  phi, C] to fit the equation
-        x[n] = A * cos(f0/fs * 2 * pi * n + phi) + C
+        x[n] = A * sin(f0/fs * 2 * pi * n + phi) + C
       where n is sample number.  Or, as a function of time:
-        x(t) = A * cos(f0 * 2 * pi * t + phi) + C
+        x(t) = A * sin(f0 * 2 * pi * t + phi) + C
 
     by Jim Paris
     (Verified to match sfit4.m)
     """
     N = len(data)
-    t = linspace(0, (N-1) / fs, N)
+    t = linspace(0, (N-1) / float(fs), N)
 
     ## Estimate frequency using FFT (step b)
     Fc = fft(data)
@@ -156,32 +228,31 @@ def sfit4(data, fs):
     i = arccos((Z2*cos(ni2) - Z1*cos(ni1)) / (Z2-Z1)) / n
 
     # Convert to Hz
-    f0 = i * fs / N
+    f0 = i * float(fs) / N
 
-    ## Fit it
+    # Fit it.  We'll catch exceptions here and just returns zeros
-    # first guess for A0, B0 using 3-parameter fit (step c)
+    # if something fails with the least squares fit, etc.
-    w = 2*pi*f0
-    D = c_[cos(w*t), sin(w*t), ones(N)]
-    s = linalg.lstsq(D, data)[0]
-
-    # Now iterate 6 times (step i)
-    for idx in range(6):
-        D = c_[cos(w*t), sin(w*t), ones(N),
-              -s[0] * t * sin(w*t) + s[1] * t * cos(w*t) ] # eqn B.16
-        s = linalg.lstsq(D, data)[0] # eqn B.18
-        w = w + s[3]	# update frequency estimate
-
-    ## Extract results
-    A = sqrt(s[0]*s[0] + s[1]*s[1]) # eqn B.21
-    f0 = w / (2*pi)
     try:
-        phi = -arctan2(s[1], s[0]) # eqn B.22
+        # first guess for A0, B0 using 3-parameter fit (step c)
-    except TypeError:
+        s = zeros(3)
+        w = 2*pi*f0
+
+        # Now iterate 7 times (step b, plus 6 iterations of step i)
+        for idx in range(7):
+            D = c_[cos(w*t), sin(w*t), ones(N),
+                  -s[0] * t * sin(w*t) + s[1] * t * cos(w*t) ] # eqn B.16
+            s = linalg.lstsq(D, data)[0] # eqn B.18
+            w = w + s[3]	# update frequency estimate
+
+        ## Extract results
+        A = sqrt(s[0]*s[0] + s[1]*s[1]) # eqn B.21
+        f0 = w / (2*pi)
+        phi = arctan2(s[0], s[1]) # eqn B.22 (flipped for sin instead of cos)
+        C = s[2]
+        return (A, f0, phi, C)
+    except Exception as e:
         # something broke down, just return zeros
         return (0, 0, 0, 0)
-    C = s[2]
-
-    return (A, f0, phi, C)
 
 if __name__ == "__main__":
     main()