1
\$\begingroup\$

I am using the least_squares() function from the scipy.optimize module to calibrate a Canopy structural dynamic model (CSDM). The calibrated model is then used to predict leaf area index (lai) based on thermal time (tt) data. I tried two variants, the first did not use the "loss" parameter of the least_squares() function, while the second set this parameter to produce a robust model. Both models give me runtime warnings even though the optimization complete successfully.

With the simple least squares model I get these warnings:

__main__:24: RuntimeWarning: overflow encountered in exp
__main__:24: RuntimeWarning: divide by zero encountered in true_divide

With the robust least squares model I get these warnings:

__main__:24: RuntimeWarning: overflow encountered in exp
__main__:24: RuntimeWarning: overflow encountered in power
# Canopy structural dynamic model (CSDM) implementation using scipy.optimize

import numpy as np
from scipy.optimize import least_squares
import matplotlib.pyplot as plt

#independent variable training data
tt_train = np.array([299.135, 408.143, 736.124, 1023.94, 1088.47, 1227.22,
1313.94, 1392.93, 1482.83, 1581.96, 2064.27, 2277.95, 2394.62, 2519.23])
#dependent variable training data
lai_train = np.array([0.0304313, 0.0833402, 0.682014, 0.973261, 2.54978,
4.93747, 5.31949, 6.25236, 6.64175, 7.3717, 3.61623, 2.96673, 1.72345, 0.803591])

# The CSDM formula (Duveiller et al. 2011. Retrieving wheat Green Area Index during the growing season...)
# LAI = k * (1 / ((1 + Exp(-a * (tt - T0 - Ta))) ^ c) - Exp(b * (tt - T0 - Tb)))

# initial estimates of parameters
To = 50      # plant emergence (x[0])
Ta = 1000    # midgrowth (x[1])
Tb = 2000    # end of cenescence (x[2])
k = 6        # scaling factor (arox. max LAI) (x[3])
a = 0.01     # rate of growth (x[4])
b = 0.01     # rate of senescence (x[5])
c = 1        # parameter allowing some plasticity to the shape of the curv (x[6])
x0 = np.array([To, Ta, Tb, k, a, b, c])

def model(x, tt):
    return x[3] * (1 / ((1 + np.exp(-x[4] * (tt - x[0] - x[1]))) ** x[6]) - np.exp(x[5] * (tt - x[0] - x[2])))

#Define the function computing residuals for least-squares minimization
def fun(x, tt, lai):
    return model(x, tt) - lai

#simple model
res_lsq = least_squares(fun, x0, args=(tt_train, lai_train))
#robust model
res_robust = least_squares(fun, x0, loss='soft_l1', f_scale=1, args=(tt_train, lai_train)) 

# termal time data for full season
tt_test = np.array([11.7584,22.1838,34.0008,47.7174,64.3092,81.1832,90.1728,101.494,116.125,127.732,140.229,
154.381,170.5,185.707,201.368,217.642,233.593,249.703,266.233,283.074,299.135,314.386,327.024,337.58,344.699,
354.328,367.247,379.627,391.51,400.93,408.143,414.941,423.678,433.2,442.072,448.923,454.699,462.479,471.187,
481.93,492.389,499.845,508.979,522.702,533.663,540.178,547.342,553.534,560.451,569.112,574.813,580.323,
589.95,597.542,601.937,606.161,609.48,613.321,615.876,619.44,623.754,630,636.784,640.978,643.625,646.384,
650.608,657.538,664.192,670.672,673.271,674.191,679.735,685.526,694.327,700.824,710.817,714.799,717.233,
718.539,718.669,718.669,718.669,719.985,726.038,736.124,740.441,745.865,751.463,757.85,761.474,763.216,
769.154,772.596,778.288,782.517,785.868,791.79,798.324,803.554,806.697,809.536,813.457,817.2,817.902,
817.902,817.902,817.902,817.902,820.271,824.126,826.609,826.668,827.619,827.619,827.629,827.629,827.629,
827.629,827.629,833.344,841.854,849.289,854.49,859.806,871.709,878.918,882.926,885.63,888.126,892.953,
898.661,899.547,900.031,903.327,906.253,909.183,912.358,917.222,921.757,925.36,927.341,927.819,929.745,
930.731,930.949,932.384,932.384,932.384,932.384,932.384,932.384,932.384,933.757,933.757,933.757,936.283,
940.396,945.01,952.758,961.418,973.865,986.804,999.508,1012.5,1023.94,1034.92,1048.68,1052.39,1052.39,
1052.8,1053.73,1053.73,1053.73,1054.09,1054.31,1056.48,1061.43,1068.88,1076.67,1088.47,1104.89,1119.38,
1130.99,1141.1,1155.06,1171.19,1185.48,1199.21,1213.17,1227.22,1242.87,1260.89,1277.97,1295.61,1313.94,
1331.04,1346.59,1359.13,1375.4,1392.93,1408.89,1424.56,1442.76,1461.92,1482.83,1502.78,1523.67,1544.39,
1563.29,1581.96,1599.23,1619.32,1637.81,1656.31,1678.33,1700.06,1721.59,1741.63,1761.09,1779.76,1799.04,
1818.54,1836.93,1855.25,1871.02,1890.62,1909.82,1928.01,1946.39,1966.27,1983.82,2003.26,2023.74,2043.92,
2064.27,2085.74,2107.14,2127.92,2148.44,2167.92,2188.01,2208.63,2231.33,2254.54,2277.95,2301.32,2323.56,
2347.52,2370.52,2394.62,2419.89,2442.6,2466.69,2492.51,2519.23,2540.78,2563.82,2585.14,2607.89,2628.95,
2652.57,2676.55,2700.73,2724,2742.09,2759.06,2778.77,2798.12,2815.01,2834.76,2855.37,2878.56])

# apply the two models to the full season data    
lai_lsq = model(res_lsq.x, tt_test)
lai_robust = model(res_robust.x, tt_test)

# plot the two model fits
plt.plot(tt_train, lai_train, 'o', markersize=4, label='training data')
plt.plot(tt_test, lai_lsq, label='fitted lsq model')
plt.plot(tt_test, lai_robust, label='fitted robust model')
plt.xlabel("tt")
plt.ylabel("LAI")
plt.legend(loc='upper left')
plt.show()

Here is an image showing the fitted lines for the two models. The simple lsq model seems OK but this overflow warnings may indicate serious problem that compromise the algorithm (especially divide by zero). The robust model at the other hand is totally wrong.

enter image description here

\$\endgroup\$
  • \$\begingroup\$ Welcome to CodeReview@SE. I see your question close to the border between what is and isn't on topic here: please explicitly express one or more concerns about the code presented. \$\endgroup\$ – greybeard Jan 24 at 8:22
  • \$\begingroup\$ An overflow in np.exp(-x[4] * (tt - x[0] - x[1]))) ** x[6]) does not sound so unexpected, x[6] does not even need to be very large, depending on the other parameters, e.g. np.exp(2**10) = inf. And then 1./np.exp(-2**10) is a divide by zero. \$\endgroup\$ – Graipher Jan 24 at 13:48
  • \$\begingroup\$ @Graipher thanks for the explanation. This equation is popular in the literature in my field. I have had some success in finding the parameters with Excel Solver but using another data set. But implementation in Python seems to be quite difficult; or maybe impossible? \$\endgroup\$ – ABC Jan 24 at 14:37
  • \$\begingroup\$ I found that there were no overflow and other warnings when constrain the model parameters to vary in certain bounds. I set the bounds parameter of least_squares like this: bounds=([1, 500, 1500, 3, 0.001, 0.001, 0.5], [50, 2000, 3000, 9, 0.1, 0.1, 1.5]). \$\endgroup\$ – ABC Jan 27 at 9:01
  • \$\begingroup\$ Also, the data sets should be cast from float32 to float64. \$\endgroup\$ – ABC Jan 27 at 9:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.