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) Ta = 1000 # midgrowth (x) Tb = 2000 # end of cenescence (x) k = 6 # scaling factor (arox. max LAI) (x) a = 0.01 # rate of growth (x) b = 0.01 # rate of senescence (x) c = 1 # parameter allowing some plasticity to the shape of the curv (x) x0 = np.array([To, Ta, Tb, k, a, b, c]) def model(x, tt): return x * (1 / ((1 + np.exp(-x * (tt - x - x))) ** x) - np.exp(x * (tt - x - x))) #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.