Add C convolution routine for speedup

-Also changes to structure of ring and unit hydrographs.  Changes are internal and are backward compatable.
-Fixed bug in agg_tsteps logic

Author: Joe Hamman

rvic/c_convolve.c

@@ -0,0 +1,52 @@
+#include <stdio.h>
+#include <stdlib.h>
+
+void c_convolve(const int nsources,             /*scalar - number of sources*/
+                const int noutlets,             /*scalar - length of subset*/
+                const int subset_length,        /*scalar - length of subset*/
+                const int x_size,
+                const int* source2outlet_ind,   /*1d array - source to outlet mapping*/
+                const int* source_y_ind,        /*1d array - source y location*/
+                const int* source_x_ind,        /*1d array - source x location*/
+                const int* source_time_offset,  /*1d array - source time offset*/
+                const double* unit_hydrograph,  /*2d array[times][sources] - unit hydrographs*/
+                const double* aggrunin,         /*2d array[ysize][xsize] - vic runoff flux*/
+                double* ring)                   /*2d array[times][outlets] - convolution ring*/
+{
+    const double* runin;                      /*pointer to sources runoff flux*/
+    int s, i, j;                              /*counters*/
+    int y, x, offset, outlet;                 /*2d indicies*/
+    int xyind, rind, uhind;                   /*1d indicies*/
+
+    /*Loop through all sources*/
+    for (s = 0; s < nsources; s++) {
+
+        outlet = source2outlet_ind[s];
+        y = source_y_ind[s];
+        x = source_x_ind[s];
+        offset = source_time_offset[s];
+
+        //1d index location
+        //2d-->1d indexing goes like this:  ind = y*x_size + x
+        xyind = y*x_size + x;
+
+        runin = &aggrunin[xyind];
+
+        // if (*runin > 10000.0) {
+        //     printf("runin is big: %f\n", *runin);
+        //     printf("xyind: %i\n", xyind);
+        //     printf("y: %i\n", y);
+        //     printf("x: %i\n", x);
+
+        /* Do the convolution */
+        for (i = 0; i < subset_length; i++) {
+            j = i + offset;
+
+            //1d index locations
+            rind = j * noutlets + outlet;
+            uhind = i * nsources + s;
+
+            ring[rind] += unit_hydrograph[uhind] * *runin;
+        }
+    }
+}

rvic/convolution_wrapper.py

@@ -0,0 +1,101 @@
+"""
+convolution_wrapper.py
+
+ctypes wrapper for convolve.c
+
+gcc -shared -o c_convolve.so c_convolve.c
+"""
+
+import numpy as np
+import ctypes
+
+import os
+# os.system('gcc -shared -o c_convolve.so c_convolve.c')
+
+try:
+    path_to_file = os.path.split(__file__)[0]
+
+    _convolution = np.ctypeslib.load_library('c_convolve.so', path_to_file)
+    args = [ctypes.c_int,
+            ctypes.c_int,
+            ctypes.c_int,
+            ctypes.c_int,
+            np.ctypeslib.ndpointer(np.int32),
+            np.ctypeslib.ndpointer(np.int32),
+            np.ctypeslib.ndpointer(np.int32),
+            np.ctypeslib.ndpointer(np.int32),
+            np.ctypeslib.ndpointer(np.float64),
+            np.ctypeslib.ndpointer(np.float64),
+            np.ctypeslib.ndpointer(np.float64)]
+    _convolution.c_convolve.argtypes = args
+    _convolution.c_convolve.restype = None
+
+    def rvic_convolve(*args):
+        """args:
+
+           nsources,
+           noutlets,
+           subset_length,
+           xsize,
+           source2outlet_ind,
+           source_y_ind,
+           source_x_ind,
+           source_time_offset,
+           unit_hydrograph,
+           aggrunin,
+           ring
+        """
+        _convolution.c_convolve(*args)
+
+        return
+except:
+    print('Using brodcasting convolution method because there was a problem '
+          'loading c_convolve.c')
+
+    def rvic_convolve(nsources, noutlets, subset_length, xsize,
+                      source2outlet_ind, source_y_ind, source_x_ind,
+                      source_time_offset, unit_hydrograph, aggrunin, ring):
+        # numpy brodcasting method
+        for s, outlet in enumerate(source2outlet_ind):
+            y = source_y_ind[s]
+            x = source_x_ind[s]
+            i = source_time_offset[s]
+            j = i + subset_length
+            ring[i:j, outlet] += np.squeeze(unit_hydrograph[:, s] * aggrunin[y, x])
+
+        # # pure python convolution
+        # for s, outlet in enumerate(source2outlet_ind):
+        #     y = source_y_ind[s]
+        #     x = source_x_ind[s]
+        #     for i in xrange(subset_length):
+        #         j = i + source_time_offset[s]
+        #         ring[j, outlet] += (unit_hydrograph[i, s] * aggrunin[y, x])
+        return
+
+
+def test():
+    nsources = 20
+    subset_length = 10
+    full_time_length = 15
+    noutlets = 4
+    source2outlet_ind = np.linspace(0, noutlets, nsources,
+                                    endpoint=False).astype(np.int32)
+
+    ysize = 12
+    xsize = 15
+
+    source_y_ind = np.random.randint(0, ysize-1, nsources).astype(np.int32)
+    source_x_ind = np.random.randint(0, xsize-1, nsources).astype(np.int32)
+
+    source_time_offset = np.random.randint(0, 4, nsources).astype(np.int32)
+
+    aggrunin = np.random.uniform(size=(ysize, xsize))
+    unit_hydrograph = np.zeros((subset_length, nsources), dtype=np.float64)
+    unit_hydrograph[0:4, :] = 0.25
+    ring = np.zeros((full_time_length, noutlets), dtype=np.float64)
+
+    for i in xrange(10):
+        aggrunin = np.random.uniform(size=(ysize, xsize))
+        rvic_convolve(nsources, noutlets, subset_length, xsize,
+                      source2outlet_ind, source_y_ind, source_x_ind,
+                      source_time_offset, unit_hydrograph, aggrunin, ring)

rvic/history.py

@@ -3,10 +3,12 @@
 
 Summary:
     This is the core history file module for the rvic model.
-    The core of the module is the Tape class.  The basic procedure is as follows:
+    The core of the module is the Tape class.  The basic procedure is as
+    follows:
         - initialization:  sets tape options, determines filenames, etc.
         - update: method that incorporates new fluxes into the history tape.
-        - __next_update_out_data: method to determine when to update the outdata container
+        - __next_update_out_data: method to determine when to update the
+        outdata container
 
 """
 import os
@@ -17,9 +19,11 @@
 from logging import getLogger
 from log import LOG_NAME
 from share import SECSPERDAY, HOURSPERDAY, TIMEUNITS, NC_INT, NC_FLOAT
-from share import NC_DOUBLE, NC_CHAR, WATERDENSITY
+from share import NC_DOUBLE, WATERDENSITY
+# from share import NC_CHAR, RVIC_TRACERS
 import share
 
+
 # -------------------------------------------------------------------- #
 # create logger
 log = getLogger(LOG_NAME)
@@ -60,8 +64,8 @@ def __init__(self, time_ord, caseid, Rvar, tape_num=0,
         self._glob_ats = glob_ats
 
         self._out_data_i = 0            # position counter for out_data array
-        self._out_times = np.zeros(self._mfilt, type=np.float64)
-        self._out_time_bnds = np.zeros((self._mfilt, 2), type=np.float64)
+        self._out_times = np.zeros(self._mfilt, dtype=np.float64)
+        self._out_time_bnds = np.zeros((self._mfilt, 2), dtype=np.float64)
 
         self.__get_rvar(Rvar)           # Get the initial Rvar fields
         self._grid_shape = grid_area.shape
@@ -189,8 +193,9 @@ def update(self, data2tape, time_ord):
         # ------------------------------------------------------------ #
         # Update the fields
         for field in self._fincl:
+            tracer = 'LIQ'
             log.debug('updating {0}'.format(field))
-            fdata = data2tape[field]
+            fdata = data2tape[field][tracer]
             if self._avgflag == 'A':
                 self._temp_data[field] += fdata
             elif self._avgflag == 'I':
@@ -237,7 +242,9 @@ def __update_out_data(self):
             if self._outtype == 'grid':
                 # ---------------------------------------------------- #
                 # Grid the fields
-                self._out_data[field][self._out_data_i, self._outlet_y_ind, self._outlet_x_ind] = self._temp_data[field][:]
+                self._out_data[field][self._out_data_i,
+                                      self._outlet_y_ind,
+                                      self._outlet_x_ind] = self._temp_data[field][:]
                 # ---------------------------------------------------- #
             else:
                 self._out_data[field][self._out_data_i, :] = self._temp_data[field]
@@ -265,7 +272,7 @@ def __update_out_data(self):
     # ---------------------------------------------------------------- #
     # Get import rvar fields
     def __get_rvar(self, rvar):
-        """ Get the rvar Fields that are useful in writing output """
+        """ Get the rvar Fields that are useful for writing output """
         self._dt = rvar.unit_hydrograph_dt
         self._num_outlets = rvar.n_outlets
         self._outlet_decomp_ind = rvar.outlet_decomp_ind

rvic/read_forcing.py

@@ -249,18 +249,21 @@ def read(self, timestamp):
             if i == 0:
                 forcing['LIQ'] = self.current_fhdl.variables[fld][self.current_tind, :, :] * self.fld_mult[fld]
             else:
-                forcing['LIQ'] += self.current_fhdl.variables[fld][self.current_tind, :, :] * self.fld_mult[fld]
+                forcing['LIQ'] += self.current_fhdl.variables[fld][self.current_tind, :, :] * self.fld_mult[fld].astype(np.float64)
 
         for i, fld in enumerate(self.ice_flds):
             self.current_fhdl.variables[fld].set_auto_maskandscale(False)
             if i == 0:
-                forcing['ICE'] = self.current_fhdl.variables[fld][self.current_tind, :, :] * self.fld_mult[fld]
+                forcing['ICE'] = self.current_fhdl.variables[fld][self.current_tind, :, :] * self.fld_mult[fld].astype(np.float64)
             else:
-                forcing['ICE'] += self.current_fhdl.variables[fld][self.current_tind, :, :] * self.fld_mult[fld]
+                forcing['ICE'] += self.current_fhdl.variables[fld][self.current_tind, :, :] * self.fld_mult[fld].astype(np.float64)
 
         # move forward one step
         self.current_tind += 1
 
+        for field in forcing:
+            forcing[field] = forcing[field].astype(np.float64)
+
         return forcing
         # ------------------------------------------------------------ #
     # ---------------------------------------------------------------- #

rvic/share.py

@@ -53,7 +53,7 @@
 RPOINTER = 'rpointer'
 
 # tracers
-RVIC_TRACERS = ('LIQ',)
+RVIC_TRACERS = ('LIQ',)  # Before changing, update history module
 
 # Calendar key number for linking with CESM
 CALENDAR_KEYS = {0: ['None'],