Uncertainty Unwrapped

Please read my disclaimer.

I'm proud to announce UncertaintyWrapper at the Cheese Shop. This work was supported by my employer, SunPower Corp. and is currently offered with a standard 3-clause BSD license. The documentation, source code and releases are also available our SunPower org GitHub page.

So what does uncertainty_wrapper do? Let's say you have a Python function, to calculate solar position, and the function uses a C/C++ library via the Python ctypes library. Or maybe you just have a really complicated set of calculations that you repeat 8760 times, and you want it to run super fast, so you don't want it to calculate derivatives and uncertainty at every internal step, just the final output. Oh and by the way, you want all 8760 calculations vectorized, _ie_: done concurrently as much as possible.

Heres an example using PVLIB of just the first 24 hours.

#
import numpy as np  # v1.11.0
import pandas as pd  # v0.18.0
import pytz  # v2016.1
import pvlib  # v0.3.2
from uncertainty_wrapper import unc_wrapper_args  # v0.4.1

PST = pytz.timezone('US/Pacific')  # Pacific Standard Time
times = pd.DatetimeIndex(start='2015/1/1', end='2015/1/2', freq='1h', tz=PST)  # date range

# arguments for the number of observations
# new in UncertaintyWrapper==0.4.1, jagged arrays are okay
latitude, longitude, pressure, altitude, temperature = 37., -122., 101325., 0., 22.

# standard deviation of 1% assuming normal distribution
covariance = np.tile(np.diag([0.0001] * 5), (times.size, 1, 1))  # tile this for the number of observations


@unc_wrapper_args(1, 2, 3, 4, 5)
# indices specify positions of independent variables:
# 1: latitude, 2: longitude, 3: pressure, 4: altitude, 5: temperature
def spa(times, latitude, longitude, pressure, altitude, temperature):
    # location class only used prior to pvlib-0.3
    dataframe = pvlib.solarposition.spa_c(times, latitude, longitude, pressure, altitude, temperature)
    retvals = dataframe.to_records()
    zenith = retvals['apparent_zenith']
    zenith = np.where(zenith<90, zenith, np.nan)
    azimuth = retvals['azimuth']
    return zenith, azimuth


ze, az, cov, jac = spa(times, latitude, longitude, pressure, altitude, temperature, __covariance__=covariance)
df = pd.DataFrame({'zenith': ze, 'az': az}, index=times)  # easier to view as dataframe
print df
#                                    az     zenith
# 2015-01-01 00:00:00-08:00  349.297715        NaN
# 2015-01-01 01:00:00-08:00   40.210628        NaN
# 2015-01-01 02:00:00-08:00   66.719304        NaN
# 2015-01-01 03:00:00-08:00   80.930185        NaN
# 2015-01-01 04:00:00-08:00   90.852887        NaN
# 2015-01-01 05:00:00-08:00   99.212426        NaN
# 2015-01-01 06:00:00-08:00  107.181217        NaN
# 2015-01-01 07:00:00-08:00  115.450451        NaN
# 2015-01-01 08:00:00-08:00  124.564183  84.113440
# 2015-01-01 09:00:00-08:00  135.023137  74.984664
# 2015-01-01 10:00:00-08:00  147.247403  67.475783
# 2015-01-01 11:00:00-08:00  161.371578  62.273878
# 2015-01-01 12:00:00-08:00  176.922804  60.008978
# 2015-01-01 13:00:00-08:00  192.742327  61.017538
# 2015-01-01 14:00:00-08:00  207.519768  65.144340
# 2015-01-01 15:00:00-08:00  220.494108  71.839001
# 2015-01-01 16:00:00-08:00  231.600910  80.422988
# 2015-01-01 17:00:00-08:00  241.184075  89.948123
# 2015-01-01 18:00:00-08:00  249.726361        NaN
# 2015-01-01 19:00:00-08:00  257.751550        NaN
# 2015-01-01 20:00:00-08:00  265.873170        NaN
# 2015-01-01 21:00:00-08:00  275.014534        NaN
# 2015-01-01 22:00:00-08:00  287.078877        NaN
# 2015-01-01 23:00:00-08:00  307.283646        NaN
# 2015-01-02 00:00:00-08:00  348.921385        NaN

# covariance at 8AM
idx = 8
print times[idx]
# Timestamp('2015-01-01 08:00:00-0800', tz='US/Pacific', offset='H')
nf = 2  # number of dependent variables: [ze, az]
print cov[idx]
# [[ 0.6617299  -0.6152971 ]
#  [-0.6152971   0.62483904]]

# standard deviation
print np.sqrt(cov[8].diagonal())
# [ 0.81346782,  0.79046761]

# Jacobian at 8AM
nargs = 5  # number of independent args
print jac[idx]
# [[  5.56456716e-01  -6.45065654e-01  -1.37538277e-06   0.00000000e+00    4.72409055e-04]
#  [  8.29163154e-02   6.47436098e-01   0.00000000e+00   0.00000000e+00    0.00000000e+00]]
#

First this tells us that the standard deviation of the zenith is 1% if the input has a standard deviation of 1%. That's reasonable. This also tells that zenith is more sensitive to latitude and longitude than pressure or temperature and more sensitive to latitude than azimuth is.