Use explicit finite difference instead of scipy.misc.derivative in pvsyst_temperature_coeff

[kandersolar]

• Closes scipy.misc.derivative is deprecated #1644
• I am familiar with the contributing guidelines
• [ ] Tests added
• [ ] Updates entries in docs/sphinx/source/reference for API changes.
• Adds description and name entries in the appropriate "what's new" file in docs/sphinx/source/whatsnew for all changes. Includes link to the GitHub Issue with :issue:`num` or this Pull Request with :pull:`num`. Includes contributor name and/or GitHub username (link with :ghuser:`user`).
• [ ] New code is fully documented. Includes numpydoc compliant docstrings, examples, and comments where necessary.
• Pull request is nearly complete and ready for detailed review.
• Maintainer: Appropriate GitHub Labels (including remote-data) and Milestone are assigned to the Pull Request and linked Issue.

pvsyst_temperature_coeff uses the default step size of 1 in its call to scipy.misc.derivative. I haven't checked the reference to see if that value is intentional, but if not, it seems a little large to me. I've reduced the step size to the rather arbitrary value of 1e-3 here, which changes the answer very slightly:

In [49]: kwargs = dict(alpha_sc=0.0054, gamma_ref=1.058, mu_gamma=0.0053585, I_L_ref=7.6626, I_o_ref=2.1003e-09, R_sh_ref=236.57, R_sh_0=886.22, R_s=0.25483, cells_in_series=36)

In [50]: pvsyst_temperature_coeff(**kwargs)  # this PR
Out[50]: 0.0014326198691227289

In [51]: pvlib.ivtools.sdm.pvsyst_temperature_coeff(**kwargs)  # v0.9.4
Out[51]: 0.0014326221939508096

[mikofski]

Thanks for this! Some ideas, tell me what you think?

1. Would you be open to moving this to a new function in tools.py or utils.py where ever appropriate so it could be reused?
2. My preference would be to define a constant DX at module top instead of hardcoding here. If defining a new function, then ‘dx’ could be an optional argument with default DX
3. Traditionally dx is the cube root of epsilon, the machine precision of the system, which you can get from numpy

[kandersolar]

a new function in tools.py

Originally I did implement it as a standalone function in tools.py, then thought it was a lot of docstring for almost no code and decided to just inline it into the pvsyst function instead. Either way is fine with me.

a constant DX at module top

I generally have the view that module-level variables are best treated as read-only. If we want dx to be configurable, I'd favor it being an optional argument to pvsyst_temperature_coeff.

Traditionally dx is the cube root of epsilon

Nice! Some googling reports that this choice is not only traditional, but optimal (not sure in what sense...). Thanks!

[mikofski]

I generally have the view that module-level variables are best treated as read-only. If we want dx to be configurable, I'd favor it being an optional argument

I agree 💯 and that's sort of what I meant like this...

# in tools.py

import numpy as np

EPS = np.finfo('float64').eps  # machine precision NumPy-1.20
DX = EPS**(1/3)  # optimal differential element


def _my_private_explicit_derivative_calc(f, x0, dx=DX):
    """no docstring required for private functions"""
    df = f(x0+dx) - f(x0-dx)
    return df / 2/ dx

I believe there are advantages to having epsilon defined here rather than inside a function. Mainly, it can be reused without recalculating.

a lot of docstring for almost no code

I think it's okay to make it private and put minimal docstring, maybe just define args, don't even really need formatting, could doc in comments

googling reports that this choice is not only traditional, but optimal

That's a relief I've just seen this in other code, and sort of followed it without knowing exactly why, or maybe I did at one point but don't remember? Anyway, I think it has to do with floating point precision or something, this may be the smallest value for dx that can still calculate an equally precise derivative (or gradient).

Anyway, I won't be a blocker. Whatever you decide I approve.

[mikofski]

lgtm! Thanks!

[adriesse]

That's a relief I've just seen this in other code, and sort of followed it without knowing exactly why, or maybe I did at one point but don't remember? Anyway, I think it has to do with floating point precision or something, this may be the smallest value for dx that can still calculate an equally precise derivative (or gradient).

If the function is very flat it seems to me you can run into problems choosing a small dx. The derivative function had a default of 1.0, which could be for this reason. Also, I often use float32 in large data sets...

[cwhanse]

Since bishop88 is vectorized it shouldn't cost any more to compute 4 values than 2. Why not x = [x0-2dx, x0-dx, x0+dx, x0+2dx], pmp = maxp(x, ...) and use the fourth order formula?

gamma_pdc = np.array([1, -8, 8, -1]) * pmp / (12 * dx)

https://www.dam.brown.edu/people/alcyew/handouts/numdiff.pdf

[cwhanse]

I didn't hit the button to submit this comment. No action required here, archiving the comment in case the issue arises.