.. _python-extension-example-with-units:

Python Extension Example with Units
-----------------------------------

Often Python packages contain extensions in C/C++ which can't be tested using
automatic differentiation. The Numdidfftools is an alternative package that can
calculate derivatives more accurately than the central finite difference
approximation.

An example of using :func:`~uncertainty_wrapper.core.unc_wrapper_args` with a
Python extension written in C/C++ is in the tests called
:func:`~uncertainty_wrapper.tests.test_uncertainty_wrapper.test_solpos`. This
test using C/C++ code from NREL that is called using Python ``ctypes`` module.

This example also demonstrates using Pint's units wrapper. When using the units
wrapper, you must use :func:`~uncertainty_wrapper.core.unc_wrapper_args` and
specify the indices of the positional arguments which corresond to the covariance
matrix. Also, two additional ``None, None`` should be appended to the units
wrapper return values because otherwise Pint uses ``zip(out_units, retvals)``
and therefore the covariance and Jacobian matrices will get dropped. ::

    @UREG.wraps(('deg', 'deg', None, None),
                (None, 'deg', 'deg', 'Pa', 'm', 'degC'))
    @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):
        dataframe = pvlib.solarposition.spa_c(times, latitude, longitude, pressure,
                                              temperature)
        retvals = dataframe.to_records()
        zenith = retvals['apparent_zenith']
        zenith = np.where(zenith < 90, zenith, np.nan)
        azimuth = retvals['azimuth']
        return zenith, azimuth

Then test it out. ::

    times = pd.DatetimeIndex(start='2015/1/1', end='2015/1/2', freq='1h',
                             tz=PST).tz_convert(UTC)
    latitude, longitude = 37.0 * UREG.deg, -122.0 * UREG.deg
    pressure, temperature = 101325.0 * UREG.Pa, UREG.Quantity(22.0, UREG.degC)
    altitude = 0.0 * UREG.m
    # standard deviation of 1% assuming normal distribution
    covariance = np.diag([0.0001] * 5)
    ze, az, cov, jac = spa(times, latitude, longitude, pressure, altitude,
                           temperature, __covariance__=covariance)

The results are::

    >>> ze
    <Quantity([         nan          nan          nan          nan          nan
              nan          nan          nan  84.10855021  74.98258957
      67.47442104  62.27279883  60.00799371  61.01651321  65.14311785
      71.83729124  80.41979434  89.92923993          nan          nan
              nan          nan          nan          nan          nan], 'deg')>

    >>> az
    <Quantity([ 349.29771499   40.21062767   66.71930375   80.93018543   90.85288686
       99.21242575  107.18121735  115.45045069  124.56418347  135.02313717
      147.24740279  161.37157806  176.92280365  192.74232655  207.51976817
      220.49410796  231.60091006  241.18407504  249.7263611   257.75154961
      265.87317048  275.01453439  287.07887655  307.28364551  348.92138471], 'deg')>

Note that Pint corrects the ambient temperature from Kelvin to Celsius and also
converted Pascals to millibar. Finally Pint appends the specified units to the
return values.

    >>> idx = 8  # covariance at 8AM
    >>> 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[(nf * idx):(nf * (idx + 1)), (nf * idx):(nf * (idx + 1))]
    [[ 0.66082282, -0.61487518],
     [-0.61487518,  0.62483904]]

    >>> print np.sqrt(cov[(nf * idx), (nf * idx)]) / ze[idx]  # standard deviation
    0.0096710802029002577

This tells us that the standard deviation of the zenith is 1% if the input has a standard deviation
of 1%. That's reasonable.

    >>> nargs = 5  # number of independent args
    >>> print jac[nf*(idx-1):nf*idx, nargs*(idx-1):nargs*idx]  # Jacobian at 8AM
    [[  5.56075143e-01,  -6.44623321e-01,  -1.42364184e-06, 0.00000000e+00,   1.06672083e-10],
     [  8.29163154e-02,   6.47436098e-01,   0.00000000e+00, 0.00000000e+00,   0.00000000e+00]]

This also tells that zenith is more sensitive to latitude and longitude than pressure or temperature
and more sensitive to latitude than azimuth is.

Perhaps the most interesting outcome is the negative covariance between Zenith and Azimuth. From
`Wikipedia <https://en.wikipedia.org/wiki/Covariance>`_

    ... when the greater values of one variable mainly correspond to the lesser values of the other,
    the covariance is negative.

In other words when the error in Zenith increases, the error in Azimuth decreases. This is not
uncommon but it's not always intuitively obvious; we generally think that to get the largest error
we should choose the largest errors for all independent variables.