Update APIs; add PEAQ Auto ITC test data

bayesitc/experiments.py

@@ -328,7 +328,7 @@ def _parse_heats(heats_file, unit):
 
         assert isinstance(heats_file, str)
         # Need python engine for skip_footer
-        dataframe = pd.read_table(heats_file, skip_footer=1, engine='python', sep='\s+', header=0)
+        dataframe = pd.read_table(heats_file, skipfooter=1, engine='python', sep='\s+', header=0)
         heats = numpy.array(dataframe['DH'])
         return Quantity(heats, unit)
 
@@ -575,10 +575,15 @@ def _plot_gaussian_baseline(self, full_x, full_y, sigma, x, y, y_pred):
 
         figure = Figure()
         canvas = FigureCanvas(figure)
-        axes = figure.add_subplot(1, 1, 1, axisbg='whitesmoke')
+        try:
+            # old API
+            axes = figure.add_subplot(1, 1, 1, axisbg='whitesmoke')
+        except AttributeError:
+            # new API
+            axes = figure.add_subplot(1, 1, 1, facecolor='whitesmoke')
 
         # Adds a 95% confidence interval to the plot
-        ExperimentMicroCal._plot_confidence_interval(axes, full_x, sigma, y_pred)
+        #ExperimentMicroCal._plot_confidence_interval(axes, full_x, sigma, y_pred) # DEBUG: Fix this
         # Entire set of data
         axes.plot(full_x, full_y, 'o', markersize=2, lw=1, color='deepskyblue', alpha=.5, label='Raw data')
         # Points for fit
@@ -615,7 +620,12 @@ def _plot_baseline_subtracted(self, x, y, raw=True, baseline=True):
 
         figure = Figure()
         canvas = FigureCanvas(figure)
-        axes1 = figure.add_subplot(1, 1, 1, axisbg='whitesmoke')
+        try:
+            # old API
+            axes1 = figure.add_subplot(1, 1, 1, axisbg='whitesmoke')
+        except AttributeError:
+            # new API
+            axes1 = figure.add_subplot(1, 1, 1, facecolor='whitesmoke')
 
         # Points for fit
         axes1.plot(x, y, 'o', color='deepskyblue', markersize=2, alpha=1, label='Baseline-subtracted data')
@@ -690,18 +700,32 @@ def fit_gaussian_process_baseline(self, fit_fraction=0.2, theta0=5.0, nugget=1.0
         full_y = numpy.array(self.differential_power).T
         y = numpy.array(y).T
 
-        # TODO look into GaussianProcessRegressor http://bit.ly/2kpUs0b
-        # current API will be deprecated as of scikit learn 0.8
-        gp = gaussian_process.GaussianProcess(regr='quadratic',
-                                              corr='squared_exponential',
-                                              theta0=theta0,
-                                              nugget=nugget,
-                                              random_start=100)
-
-        # Fit only based on the reduced set of the data
-        gp.fit(x, y)
-        y_pred, mean_squared_error = gp.predict(full_x, eval_MSE=True)
-        sigma = numpy.sqrt(mean_squared_error)
+        try:
+            # current API will be deprecated as of scikit learn 0.8
+            gp = gaussian_process.GaussianProcess(regr='quadratic',
+                                                  corr='squared_exponential',
+                                                  theta0=theta0,
+                                                  nugget=nugget,
+                                                  random_start=100)
+            # Fit only based on the reduced set of the data
+            gp.fit(x, y)
+            y_pred, mean_squared_error = gp.predict(full_x, eval_MSE=True)
+            sigma = numpy.sqrt(mean_squared_error)
+        except AttributeError:
+            # newer scikit-learn
+            kernel = gaussian_process.kernels.RBF(length_scale=theta0)
+            #from sklearn.gaussian_process.kernels import ConstantKernel, RBF
+            #kernel = y.std() * RBF(length_scale=theta0)
+            #gp = gaussian_process.GaussianProcessRegressor(kernel=kernel,
+            #                                               alpha=nugget*1e-6,
+            #                                               n_restarts_optimizer=100)
+            gp = gaussian_process.GaussianProcessRegressor(kernel=kernel,
+                                                           alpha=1e-3 * y.std(),
+                                                           n_restarts_optimizer=100)
+            # Fit only based on the reduced set of the data
+            gp.fit(x, y)
+            y_pred, mean_squared_error = gp.predict(full_x, return_std=True)
+            sigma = numpy.sqrt(mean_squared_error)
 
         self.baseline_power = Quantity(y_pred, 'microcalories per second')
         self.baseline_fit_data = {'x': full_x, 'y': y_pred, 'indices': fit_indices}

bayesitc/report.py

@@ -183,20 +183,29 @@ def plot_experiment(name, experiment):
     pylab.subplot(411)
 
     # plot baseline fit
-    pylab.hold(True)
+    try:
+        pylab.hold(True) # old API
+    except AttributeError:
+        pass # new API
     pylab.plot(experiment.filter_period_end_time / ureg.second, experiment.baseline_power / (ureg.microcalorie/ureg.second), 'g-')  # plot baseline fit
 
     # differential power
     pylab.plot(experiment.filter_period_end_time / ureg.second, experiment.differential_power / (ureg.microcalorie/ureg.second), 'k.', markersize=1)
     # plot red vertical line to mark injection times
-    pylab.hold(True)
+    try:
+        pylab.hold(True) # old API
+    except AttributeError:
+        pass # new API
     [xmin, xmax, ymin, ymax] = pylab.axis()
 
     for injection in experiment.injections:
         last_index = injection.first_index # timepoint at start of syringe injection
         t = experiment.filter_period_end_time[last_index] / ureg.second
         pylab.plot([t, t], [ymin, ymax], 'r-')
-    pylab.hold(False)
+    try:
+        pylab.hold(True) # old API
+    except AttributeError:
+        pass # new API
     xlabel = pylab.xlabel('time / s')
     xlabel.set_fontsize(fontsize)
     ylabel = pylab.ylabel('differential power / ucal/s')
@@ -214,7 +223,10 @@ def plot_experiment(name, experiment):
 
     # Plot enthalpogram.
     pylab.subplot(412)
-    pylab.hold(True)
+    try:
+        pylab.hold(True) # old API
+    except AttributeError:
+        pass # new API
     for injection in experiment.injections:
         # determine initial and final samples for injection i
         first_index = injection.first_index # index of timepoint for first filtered differential power measurement
@@ -229,14 +241,20 @@ def plot_experiment(name, experiment):
     # adjust axes to match first plot
     [xmin_new, xmax_new, ymin, ymax] = pylab.axis()
     pylab.axis([xmin, xmax, ymin, ymax])
-    pylab.hold(False)
+    try:
+        pylab.hold(True) # old API
+    except AttributeError:
+        pass # new API
     #pylab.title('evolved heat per injection')
     xlabel = pylab.xlabel('time / s')
     xlabel.set_fontsize(fontsize)
     ylabel = pylab.ylabel('evolved heat / ucal')
     ylabel.set_fontsize(fontsize)
     # plot zero
-    pylab.hold(True)
+    try:
+        pylab.hold(True) # old API
+    except AttributeError:
+        pass # new API
     pylab.plot(experiment.filter_period_end_time / ureg.second, 0.0*experiment.filter_period_end_time / (ureg.microcalorie/ureg.second), 'g-') # plot zero line
 
 
@@ -248,12 +266,18 @@ def plot_experiment(name, experiment):
 
     # Plot cell temperature.
     pylab.subplot(413)
-    pylab.hold(True)
+    try:
+        pylab.hold(True) # old API
+    except AttributeError:
+        pass # new API
     pylab.plot(experiment.filter_period_end_time / ureg.second, experiment.cell_temperature/ureg.kelvin - 273.15, 'r.', markersize=1)
     # adjust axes to match first plot
     [xmin_new, xmax_new, ymin, ymax] = pylab.axis()
     pylab.axis([xmin, xmax, ymin, ymax])
-    pylab.hold(False)
+    try:
+        pylab.hold(True) # old API
+    except AttributeError:
+        pass # new API
     #pylab.title('evolved heat per injection')
     xlabel = pylab.xlabel('time / s')
     xlabel.set_fontsize(fontsize)
@@ -268,12 +292,18 @@ def plot_experiment(name, experiment):
 
     # Plot adiabatic jacket temperature.
     pylab.subplot(414)
-    pylab.hold(True)
+    try:
+        pylab.hold(True) # old API
+    except AttributeError:
+        pass # new API
     pylab.plot(experiment.filter_period_end_time / ureg.second, experiment.jacket_temperature/ureg.kelvin - 273.15, 'b.', markersize=1)
     # adjust axes to match first plot
     [xmin_new, xmax_new, ymin, ymax] = pylab.axis()
     pylab.axis([xmin, xmax, ymin, ymax])
-    pylab.hold(False)
+    try:
+        pylab.hold(True) # old API
+    except AttributeError:
+        pass # new API
     #pylab.title('evolved heat per injection')
     xlabel = pylab.xlabel('time / s')
     xlabel.set_fontsize(fontsize)