diff --git a/example/EpyMix with Septoria.ipynb b/example/EpyMix with Septoria.ipynb index 50d55d8..044330a 100644 --- a/example/EpyMix with Septoria.ipynb +++ b/example/EpyMix with Septoria.ipynb @@ -544,7 +544,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3 (ipykernel)", + "display_name": "Python 3", "language": "python", "name": "python3" }, @@ -558,7 +558,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.6" + "version": "3.9.2" } }, "nbformat": 4, diff --git a/example/app_epymix.ipynb b/example/app_epymix.ipynb new file mode 100644 index 0000000..c6d012a --- /dev/null +++ b/example/app_epymix.ipynb @@ -0,0 +1,306 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "b50cd883", + "metadata": {}, + "source": [ + "# Jupyter App" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "e05e4afa", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import pandas as pd\n", + "import matplotlib.pyplot as plt\n", + "import ipywidgets\n", + "\n", + "from epymix.rain import rain as _rain \n", + "from epymix.inoculum import inoculum \n", + "from epymix.configuration import configuration\n", + "from epymix.SEIR import SEIR \n", + "from epymix.dispersion_gradient import dispersion_kernel_rust, dispersion_kernel_septo\n", + "from epymix.growth_companion import growth_pois\n", + "from epymix.SEIR_practices_data import SEIR_practice\n", + "from epymix.app_parameter_generator import parameter_set" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "51d30dca", + "metadata": {}, + "outputs": [], + "source": [ + "#disease_par='septo_par'\n", + "delta_t0 = 10 # constant, the model has been parameterise such as t = 10 degree-day\n", + "delta_t = 10 # time step\n", + "n_season = 1 # number of season\n", + "season = 250*int(delta_t0/delta_t) # 2500 dd %% length of a cropping season\n", + "delta_companion = 0 # growth start lag of the companion crop (dd), negative: companion before, positive: companion after but for i" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "ipywidgets.interact(app, crop=['pure crop', 'intercropping'], experiment=[('fababean_N0',61),('fababean_N1',63),('pea_N0',13),('pea_N1',14)], disease_par=['septo_par', 'rust_par'], year=list(range(1993, 2016)), wheat_fraction=(0,1,0.05), \n", + "delta_companion=(-100,100,1))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0d0d90ec", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "epyland", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.2" + }, + "vscode": { + "interpreter": { + "hash": "cc495b77a574b98b57b7f5732d92c368da0bdafbfc1ca32d8caf81cd26f9b127" + } + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/example/app_epymix2.ipynb b/example/app_epymix2.ipynb new file mode 100644 index 0000000..7bc8db2 --- /dev/null +++ b/example/app_epymix2.ipynb @@ -0,0 +1,691 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "b50cd883", + "metadata": { + "extensions": { + "jupyter_dashboards": { + "version": 1, + "views": { + "default_view": { + "col": 0, + "height": 3, + "row": 0, + "width": 8 + } + } + } + } + }, + "source": [ + "# Epymix : mélange de cultures contre épidemies fongiques" + ] + }, + { + "cell_type": "code", + "execution_count": 481, + "id": "e05e4afa", + "metadata": { + "extensions": { + "jupyter_dashboards": { + "version": 1, + "views": { + "default_view": { + "hidden": true + } + } + } + } + }, + "outputs": [], + "source": [ + "#python jupyterflex.py app_epymix2.ipynb\n", + "import numpy as np\n", + "import pandas as pd\n", + "import matplotlib.pyplot as plt\n", + "import ipywidgets\n", + "from IPython.display import clear_output\n", + "\n", + "\n", + "from epymix.rain import rain as _rain \n", + "from epymix.inoculum import inoculum \n", + "from epymix.configuration import configuration\n", + "from epymix.SEIR import SEIR \n", + "from epymix.dispersion_gradient import dispersion_kernel_rust, dispersion_kernel_septo\n", + "from epymix.growth_companion import growth_pois\n", + "from epymix.SEIR_practices_data import SEIR_practice\n", + "from epymix.app_parameter_generator import parameter_set" + ] + }, + { + "cell_type": "code", + "execution_count": 482, + "id": "51d30dca", + "metadata": { + "extensions": { + "jupyter_dashboards": { + "version": 1, + "views": { + "default_view": { + "hidden": true + } + } + } + } + }, + "outputs": [], + "source": [ + "#disease_par='septo_par'\n", + "delta_t0 = 10 # constant, the model has been parameterise such as t = 10 degree-day\n", + "delta_t = 10 # time step\n", + "n_season = 1 # number of season\n", + "season = 250*int(delta_t0/delta_t) # 2500 dd %% length of a cropping season\n", + "delta_companion = 0 # growth start lag of the companion crop (dd), negative: companion before, positive: companion after but for i I + g_pre = np.concatenate([np.zeros((lambd-delta_ei-1)), np.arange( (1.0/(2*delta_ei+1)), 1+(1.0/(2*delta_ei+1)), (1.0/(2*delta_ei+1)) )]) + g = np.reshape(g_pre, (lambd+delta_ei, 1, 1)) * np.ones((lambd + delta_ei, 1, 1)) + + + ### MAIN FUNCTION + for i in range(0, t): + # Function initialisation + if np.mod(i, season) == 0: + ### ble = matrix (Lx, Ly), represents the wheat proportion in landscape fields + ### ex: [0 0 0 ; 1 1 1 ; 0.5 0.5 0.5] represents one landscape + ### 1st row of fields: cropped without wheat; 2nd: wheat only; 3rd: mix with 50% of wheat + #res_frac = arrangement[int((np.floor(i / season) % Lr)), :, :] # arrangement(:,:, year of arrangement) + #LAI_K_ble = res_frac * LAI_K + #LAI_K_pois = (1 - res_frac) * LAI_K # If we suppose that LAI of wheat and the companion species are equivalent + + ### Canopy initialisation + #sth0 = LAI_K_ble / 1000.0 + #pth0 = LAI_K_pois / 1000.0 + + ### INITIALISATION OF STATE VARIABLES + Eta = np.zeros((lambd + delta_ei, Lx, Ly)) + Nsp[i, :, :] = 0 * np.ones((Lx, Ly)) + Pth[i, :, :] = lai_comp_inc[i] * np.ones((Lx, Ly)) + Poi[i, :, :] = lai_comp_inc[i] * np.ones((Lx, Ly)) + Sth[i, :, :] = lai_wheat_inc[i] * np.ones((Lx, Ly)) + Sus[i, :, :] = lai_wheat_inc[i] * np.ones((Lx, Ly)) + Lat[i, :, :] = 0 * np.ones((Lx, Ly)) + Ifc[i, :, :] = 0 * np.ones((Lx, Ly)) + Ifv[i, :, :] = 0 * np.ones((Lx, Ly)) + Rem[i, :, :] = 0 * np.ones((Lx, Ly)) + LAI[i, :, :] = Sus[i, :, :] + Lat[i, :, :] + Ifc[i, :, :] + Rem[i, :, :] + Pth[i, :, :] + LAI_alive[i, :, :] = Sus[i, :, :] + Lat[i, :, :] + Ifc[i, :, :] + LAI_wheat[i, :, :] = Sus[i, :, :] + Lat[i, :, :] + Ifc[i, :, :] + Rem[i, :, :] + if i != 0 and theta != 0: # theta = 0.01 + if disease == 'rust': + Poo[i, :, :] = theta * (Poo[i - 1, :, :] + sigma * Ifc[i - 1, :, :]) + if disease == 'septo': + Poo[i, :, :] = theta * (Poo[i - 1, :, :] + (sigma + sigma_asco) * Ifc[i - 1, :, :]) + else: + Poo[i, :, :] = inoc_init * np.ones((Lx, Ly)) + if h_wheat == 0: + Eps[i, :, :] = 1 - np.exp(-ber_wheat*LAI_alive[i, :, :] - ber_companion*Poi[i, :, :]) + else: + DomFac_wheat = h_wheat / (h_wheat + h_companion) # wheat dominance factor + DomFac_comp = h_companion / (h_wheat + h_companion) # wheat dominance factor + F_wheat_d[i,:,:] = 1 - np.exp(-ber_wheat * LAI_alive[i,:,:]) # wheat dominant + F_comp_s[i,:,:] = np.exp(-ber_wheat * LAI_alive[i,:,:]) * (1 - np.exp(-ber_companion * Poi[i,:,:])) # companion under + F_comp_d[i,:,:] = 1 - np.exp(-ber_companion * Poi[i,:,:]) # companion dominant + F_wheat_s[i,:,:] = np.exp(-ber_companion * Poi[i,:,:]) * (1 - np.exp(-ber_wheat * LAI_alive[i,:,:])) # wheat under + F_wheat[i,:,:] = F_wheat_s[i,:,:] + DomFac_wheat * (F_wheat_d[i,:,:] - F_wheat_s[i,:,:]) + F_comp[i,:,:] = F_comp_s[i,:,:] + DomFac_comp * (F_comp_d[i,:,:] - F_comp_s[i,:,:]) + Eps[i,:,:] = F_wheat[i,:,:] + F_comp[i,:,:] + AUDPC[i, :, :] = 0 * np.ones((Lx, Ly)) + Scont[i, :, :] = 0 * np.ones((Lx, Ly)) + + # Main function after initialisation + else: + pluie = np.mod(i, season) in rain[:, np.mod(int(np.floor(i/season)), rain.shape[1])] #pour les évènement de dispersion de la septo + + pth = Pth[i-1, :, :] + poi = Poi[i-1, :, :] + sth = Sth[i-1, :, :] + sus = Sus[i-1, :, :] + lat = Lat[i-1, :, :] + ifc = Ifc[i-1, :, :] + ifv = Ifv[i-1, :, :] + rem = Rem[i-1, :, :] + poo = Poo[i-1, :, :] + lai = LAI[i-1, :, :] + lai_alive = LAI_alive[i-1, :, :] + lai_wheat = LAI_wheat[i-1, :, :] + eta = Eta + nsp = Nsp[i-1, :, :] + audpc = AUDPC[i-1, :, :] + + ### SPORE POOL: production and dispersion at the field-level + # (for septoriose, this takes also into account ascospores coming from the fields) + ng_field_asco = 0 + ng_field_pycnid = 0 + ng_field_ure = 0 + if disease == "septo": + #ascospores + if ((i % season) % (2*int(delta_t0/delta_t))) == 0: # ascospore dispersion every 20 dd (explaining the 2 modulus) + ng_field_asco = ndimage.convolve(ifc * sigma_asco, kernel_asco, mode='mirror') + #pycnidiospores + if pluie == 1: # pycnidiospore dispersion during rain events + ng_field_pycnid = ndimage.convolve(ifc * sigma, kernel_pycnid, mode='mirror') + elif disease == "rust": + #urediospores + if ((i % season) % (2*int(delta_t0/delta_t))) == 0: # urediospore dispersion every 20 dd (explaining the 2 modulus) + ng_field_ure = ndimage.convolve(ifc * sigma, kernel_ure, mode='mirror') + + ### SPORE POOL: Soil inoculum (Poo remobilisation) + if disease == "septo": + inoc = max(0,(1-np.mod(i,season)/(70*int(delta_t0/delta_t)))) * poo * (np.mod(i,season) >= inf_begin) * pluie + #spores_parcelle = sigma * ifc * pluie + elif disease == "rust": + inoc = poo * (np.mod(i, season) >= inf_begin) + #spores_parcelle = (1 - alpha) * sigma * ifc * (np.mod(np.mod(i, season), 2*int(delta_t0/delta_t)) == 0) + + ### SPORE POOL: External cloud infecting the field + if ((i % season) % (2 * int(delta_t0 / delta_t))) == 0: + #if isinstance(ng_ext0, int) == 1: + if disease == 'rust': + ng_ext = ng_ext0 * (np.mod(i, season) >= inf_begin and #i < season and + np.mod(i, season) < (20*int(delta_t0/delta_t) + inf_begin)) + if disease == 'septo': + ng_ext = ng_ext0 * (np.mod(i, season) >= inf_begin and + np.mod(i, season) < (20 * int(delta_t0 / delta_t) + inf_begin)) + # else: + # ng_ext = ng_ext0[np.mod(int(np.floor(i / season)), ng_ext0.shape[0]), :, :] * ( + # np.mod(i, season) >= inf_begin and + # np.mod(i, season) < (20*int(delta_t0/delta_t) + inf_begin)) + else: + ng_ext = 0 + + ### TOTAL SPORE POOL + if disease == 'septo': + nsp = inoc + ng_ext + ng_field_asco + ng_field_pycnid # + spores_parcelle + if disease == 'rust': + nsp = inoc + ng_ext + ng_field_ure # + spores_parcelle + + ### SPORE INTERCEPTION by healthy canopy (Poisson Law) + eps = Eps[i - 1, :, :] + if h_wheat == 0: #without ERIN + if ber_companion == 0: + fac_int_wheat = 1 + else: + fac_int_wheat = (ber_wheat * lai_alive) / (ber_wheat * lai_alive + ber_companion * poi) + fac_int_comp = 1 - fac_int_wheat + if 0 in lai_alive: + ex = np.zeros((Lx,Ly)) + for j in range(0,Lx): + for k in range (0, Ly): + if lai_alive[j,k]==0: + ex[j,k] = 1 + else: + ex[j,k] = np.exp(-(fac_int_wheat[j,k] * eps[j,k] * pi_inf0 * nsp[j,k] * s0/lai_alive[j,k])) + else: + ex = np.exp(-(fac_int_wheat * eps * pi_inf0 * nsp * s0 / lai_alive)) + else: # with ERIN + f_wheat = F_wheat[i - 1, :, :] + if 0 in lai_alive: + ex = np.zeros((Lx,Ly)) + for j in range(0,Lx): + for k in range (0, Ly): + if lai_alive[j,k]==0: + ex[j,k] = 1 + else: + ex[j,k] = np.exp(-(f_wheat[j,k] * pi_inf0 * nsp[j,k] * s0/lai_alive[j,k])) + else: + ex = np.exp(-(f_wheat * pi_inf0 * nsp * s0 / lai_alive)) + ### Healthy surface contaminated + Sc = (1 - ex) * sus + + + ### Crop growth + #crS = beta_wheat * (sth) * (1 - sth / (LAI_K_ble + (LAI_K_ble == 0))) * (np.mod(i, season) < end_wheat) + #crP = beta_companion * (pth) * (1 - pth / (LAI_K_pois + (LAI_K_pois == 0))) * (np.mod(i, season) < end_companion) + + ### COMPUTING STATE VARIABLES + # Pth[i,:,:] = pth + crP + # Poi[i,:,:] = poi + crP - mu * poi * (np.mod(i,season) >= end_companion) + # Sth[i,:,:] = sth + crS - mu_wheat * sth * (np.mod(i, season) >= end_wheat) + # Sus[i,:,:] = np.maximum(0, sus + crS * pow(((sus+lat)/(sth+(sth==0))),gamma) - Sc - mu_wheat*sus*(np.mod(i,season)>=end_wheat)) + + Pth[i,:,:] = pth + lai_comp_inc[i] + Poi[i,:,:] = poi + lai_comp_inc[i] + Sth[i,:,:] = sth + lai_wheat_inc[i] + Sus[i,:,:] = np.maximum(0, sus + lai_wheat_inc[i] * pow(((sus+lat)/(sth+(sth==0))),gamma) - Sc) + + Eta[1:lambd+delta_ei,:,:] = eta[0:-1,:,:] * (1 - g[0:-1,:,:]) + Eta[0,:,:] = Sc + Eta = (1 - lai_wheat_mu[i] * (lai_wheat_mu[i] > 0)) * Eta + Lat[i,:,:] = Eta.sum(axis=0) + + if disease == "septo": + Ifc[i, :, :] = ifc * (1 - psi * pluie) - nu * ifc + np.sum((eta * g), axis=0) + Rem[i,:,:] = rem + ifc * psi * pluie + np.sum((lai_wheat_mu[i] * (lai_wheat_mu[i] > 0) * Eta), axis=0) +\ + lai_wheat_mu[i] * (lai_wheat_mu[i] > 0) * sus + nu * ifc + elif disease == "rust": + Ifc[i,:,:] = ifc - nu * ifc + np.sum((eta * g), axis=0) + Rem[i,:,:] = rem + np.sum(lai_wheat_mu[i] * (lai_wheat_mu[i] > 0) * Eta, axis=0) +\ + lai_wheat_mu[i] * (lai_wheat_mu[i] > 0) * sus + nu * ifc + + Ifv[i,:,:] = ifv + np.sum((eta * g), axis=0) + LAI[i,:,:] = Sus[i,:,:] + Lat[i,:,:] + Ifc[i,:,:] + Rem[i,:,:] + Pth[i,:,:] + LAI_alive[i, :, :] = Sus[i, :, :] + Lat[i, :, :] + Ifc[i, :, :] + LAI_wheat[i, :, :] = Sus[i, :, :] + Lat[i, :, :] + Ifc[i, :, :] + Rem[i, :, :] + Poo[i, :, :] = poo - inoc + (1 - eps) * nsp + (1 - pi_inf0) * eps * nsp - rho * poo + AUDPC[i,:,:] = audpc + Ifv[i,:,:] + Nsp[i,:,:] = nsp + # if ber_companion == 0: + # Eps[i,:,:] = 1 - np.exp(-ber_wheat*LAI[i,:,:]) + # else: + if h_wheat == 0: + Eps[i,:,:] = 1 - np.exp((-ber_wheat*LAI_alive[i,:,:]) -(ber_companion*Poi[i,:,:])) + else: + DomFac_wheat = h_wheat / (h_wheat + h_companion) # wheat dominance factor + DomFac_comp = h_companion / (h_wheat + h_companion) # wheat dominance factor + F_wheat_d[i,:,:] = 1 - np.exp(-ber_wheat * LAI_alive[i,:,:]) # wheat dominant + F_comp_s[i,:,:] = np.exp(-ber_wheat * LAI_alive[i,:,:]) * (1 - np.exp(-ber_companion * Poi[i,:,:])) # companion under + F_comp_d[i,:,:] = 1 - np.exp(-ber_companion * Poi[i,:,:]) # companion dominant + F_wheat_s[i,:,:] = np.exp(-ber_companion * Poi[i,:,:]) * (1 - np.exp(-ber_wheat * LAI_alive[i,:,:])) # wheat under + F_wheat[i,:,:] = F_wheat_s[i,:,:] + DomFac_wheat * (F_wheat_d[i,:,:] - F_wheat_s[i,:,:]) + F_comp[i,:,:] = F_comp_s[i,:,:] + DomFac_comp * (F_comp_d[i,:,:] - F_comp_s[i,:,:]) + Eps[i,:,:] = F_wheat[i,:,:] + F_comp[i,:,:] + Scont[i,:,:] = Sc + + return Nsp, Pth, Poi, Sth, Sus, Lat, Ifc, Ifv, Rem, LAI, LAI_wheat, Poo, Eps, AUDPC, Scont diff --git a/src/epymix/app_parameter_generator.py b/src/epymix/app_parameter_generator.py new file mode 100644 index 0000000..a0dad86 --- /dev/null +++ b/src/epymix/app_parameter_generator.py @@ -0,0 +1,67 @@ +import numpy as np + +from epymix.rain import rain as rain ## f_rain +from epymix.inoculum import inoculum ## inoculum +from epymix.configuration import configuration +from epymix.SEIR import SEIR ## SEIR fonction principale +from epymix.dispersion_gradient import dispersion_kernel_rust, dispersion_kernel_septo +from epymix.growth_companion import growth_pois + +def parameter_set(disease_par, delta_t0, delta_t): + if disease_par == "septo_par": + ## MAIN INPUT PARAMETERS + disease = "septo" + ## E-I TRANSITION PARAMETERS + delta_ei = 5*int(delta_t0/delta_t) + lambd = 20*int(delta_t0/delta_t) # 200 dd %% latent duration + ## EPIDEMIC PARAMETERS + s0 = 0.0001 # 0.0001 %% lesion size (µm²) + pi_inf0 = 0.0002 # 0.0002 %% spore probability infection + rho = 0.002*delta_t/delta_t0 # 0.01 %% spore mortality in P reservoir + psi = 0.3 # 0.3 %% pycnide emptying rate by rain (only for septoriose) + gamma = 0 # 0 %% virulence parameter + theta = 0.15 # 0.015 %% spore intercropping survival rate + sigma = 45000000 # 50000000 %% spore production rate (septo: pycnidiospores by rain) + sigma_asco = 0.2 * sigma # 0.2*sigma %% (parameter for septoriose only, ascospores by wind) + inf_begin = 0*int(delta_t0/delta_t) # 1000 dd %% date of epidemic start (generally between 80 et 130 dd for rust) + ### INOCULUM PARAMETERS + inoc_init_abs = 20000000 # inoc_init_abs: initial inoculum intensity + ng_ext0_abs = int(20000*delta_t/delta_t0) # ng_ext0_abs: external cloud intensity + ### SPORE DISPERSION + day_length = (24/(2*delta_t))*60*60 + alpha_asco = 3 # coefficient of dispersal (n m-1;Fitt et al 1987) + radius_asco = 5 + alpha_pycnid = 0.2 * 0.0001 # coefficient of dispersal (m² s-1; Yang et al 1991) + radius_pycnid = 5 + alpha_ure = 3 + radius_ure = 5 + + elif disease_par == "rust_par": + ## MAIN INPUT PARAMETERS + disease = "rust" + ## E-I TRANSITION PARAMETERS + delta_ei = 5*int(delta_t0/delta_t) + lambd = 10*int(delta_t0/delta_t) # 100 dd %% latent duration + ## EPIDEMIC PARAMETERS + s0 = 0.0001 # 0.0001 %% lesion size (µm²) + pi_inf0 = 0.0002*delta_t/delta_t0 # 0.0002 %% spore probability infection + rho = 0.01*delta_t/delta_t0 # 0.01 %% spore mortality in P reservoir + psi = 0.3 # 0 %% pycnide emptying rate by rain (only for septoriose) + gamma = 0 # 0 %% virulence parameter + theta = 0.01 # 0.01 %% spore intercropping survival rate + sigma = 1500000 # 1500000 %% spore production rate (LAI-1; septo: pycnidiospores by rain)) + sigma_asco = 0.02 * sigma # 0.2*sigma %% (parameter for septoriose only, ascospores by wind) + inf_begin = 80*int(delta_t0/delta_t) # 1000 dd %% date of epidemic start (generally between 80 et 130 dd for rust) + ### INOCULUM PARAMETERS + inoc_init_abs = 500000 # inoc_init_abs: initial inoculum intensity + ng_ext0_abs = int(2000*delta_t/delta_t0) # ng_ext0_abs: external cloud intensity + ### SPORE DISPERSION + day_length = (24/(2*delta_t))*60*60 + alpha_asco = 3 # 3 coefficient of dispersal (n m-1;Fitt et al 1987) + radius_asco = 5 + alpha_pycnid = 0.2 * 0.0001 # coefficient of dispersal (m² s-1; Yang et al 1991) + radius_pycnid =5 + alpha_ure = 3 + radius_ure = 5 + + return disease, delta_ei, lambd, s0, pi_inf0, rho, psi, gamma, theta, sigma, sigma_asco, inf_begin, inoc_init_abs, ng_ext0_abs, day_length, alpha_asco, radius_asco, alpha_pycnid, radius_pycnid, alpha_ure, radius_ure diff --git a/src/epymix/app_practices_data_septo.py b/src/epymix/app_practices_data_septo.py new file mode 100644 index 0000000..2310b4f --- /dev/null +++ b/src/epymix/app_practices_data_septo.py @@ -0,0 +1,141 @@ +import os +import numpy as np +import math +import pandas as pd +import matplotlib.pyplot as plt + +from epymix.rain import rain as rain +from epymix.inoculum import inoculum +from epymix.configuration import configuration +from epymix.SEIR_practices_data import SEIR +from epymix.dispersion_gradient import dispersion_kernel_rust, dispersion_kernel_septo +from epymix.growth_companion import growth_pois + +### CHECK METEO TIME_STEP EQUALIZATION !!!! +### CHECK METEO TIME_STEP EQUALIZATION !!!! +### CHECK METEO TIME_STEP EQUALIZATION !!!! + +### OPEN DATA +data_lai_wheat = pd.read_csv('spline_coordinate_inc_wheat_10.csv', sep=';') +data_lai_comp = pd.read_csv('spline_coordinate_inc_comp_10.csv', sep=';') +data_lai_wheat_mu = pd.read_csv('spline_coordinate_coef_wheat_10.csv', sep=';') + +### BE CAREFUL +### CHECK Lx, Ly, RDA, DELTA, BER, SCENARIO + +## MAIN INPUT PARAMETERS +experiment = 5 +data_lai_wheat.columns.values[experiment] + +disease = "septo" +delta_t0 = 10 # constant, the model has been parameterise such as t = 10 degree-day +delta_t = 10 # time step +n_season = 1 # number of season +t = 136 +#t = len(data_lai_wheat) - data_lai_wheat.iloc[:, experiment].isna().sum() +season = t # 250*int(delta_t0/delta_t) # 2500 dd %% length of a cropping season + + +### ROTATION SCENARIO (Lx,Ly,Lr) +### f_rotation(Lr, Lx, Ly, scenario, frac_ble=1, melange=1) +### f_configuration(Lr, Lx, Ly, scenario, resistant_fraction) +Lr=1; Lx=1; Ly=1; scenario_rot='uniform'; wheat_fraction=0.5 +arrangement = configuration(Lr, Lx, Ly, scenario_rot, wheat_fraction) + +## GROWTH PARAMETERS +lai_wheat_inc = np.array(data_lai_wheat.iloc[:, experiment]) +lai_comp_inc = np.array(data_lai_comp.iloc[:, experiment]) +lai_wheat_mu = np.array(data_lai_wheat_mu.iloc[:, experiment]) +# mu_wheat = 0.003 * delta_t / delta_t0 # 0.03 %% mortality rate of S and E tissues (LAI/10dd) +nu = 0.01 * delta_t / delta_t0 # nu = mu %% mortality of I infectious tissues (LAI/10dd) +# mu_companion = 0.03 * delta_t / delta_t0 # 0.03 %% mortality rate of the companion species (LAI/10dd) +# beta_wheat = 0.09 * delta_t / delta_t0 # 0.09 %% wheat growth parameter (LAI/10dd) +# beta_companion = 0.09 * delta_t / delta_t0 # beta_companion = beta_wheat %% growth parameter of the companion crop (LAI/10dd) +# end_wheat = 140 * int(delta_t0 / delta_t) # 1400 dj %% date of wheat growth end +# end_companion = 140 * int(delta_t0 / delta_t) # end_companion = end_wheat %% date of growth end for the companion crop +# LAI_K = 6 # 6 %% carrying capacity (Maximum canopy LAI) +ber_wheat = 1 # 1 # wheat spore interception coefficient (the Beer-Lambert law) +ber_companion = 1 # 1 # spore interception coefficient of the companion species (the Beer-Lambert law) +h_wheat = 1 # wheat height +h_companion = 1 # companion height + +## E-I TRANSITION PARAMETERS +delta_ei = 5*int(delta_t0/delta_t) +lambd = 20*int(delta_t0/delta_t) # 200 dd %% latent duration + +## EPIDEMIC PARAMETERS +s0 = 0.0001 # 0.0001 %% lesion size (µm²) +pi_inf0 = 0.0002 # 0.0002 %% spore probability infection +rho = 0.002*delta_t/delta_t0 # 0.01 %% spore mortality in P reservoir +psi = 0.3 # 0.3 %% pycnide emptying rate by rain (only for septoriose) +gamma = 0 # 0 %% virulence parameter +theta = 0.15 # 0.015 %% spore intercropping survival rate +sigma = 45000000 # 50000000 %% spore production rate (septo: pycnidiospores by rain) +sigma_asco = 0.2 * sigma # 0.2*sigma %% (parameter for septoriose only, ascospores by wind) +inf_begin = 0*int(delta_t0/delta_t) # 1000 dd %% date of epidemic start (generally between 80 et 130 dd for rust) + +### RAIN PARAMETER +### f_rain, return rain +year=2000 #1995: défavorable; 1997: moyenne, 2000: très favorable +rain = rain(year, n_season, delta_t) + +### INOCULUM PARAMETERS +### inoculum(scenario, frac_inf, inoc_init_abs, ng_ext0_abs, rotation) +### return inoc_init, ng_ext0 +# inoc_init_abs = 20000000 # 20000000 %% spores initially present in reservoir P +# ng_ext0_abs = int(20000*delta_t/delta_t0) # 20000 %% spores coming from a cloud external to the landscape +scenario_ino = 'initial_inoculum'; frac_inf = 1; inoc_init_abs = 20000000; ng_ext0_abs = int(20000*delta_t/delta_t0) +inoc_init, ng_ext0 = inoculum(scenario_ino, Lx, Ly, frac_inf, inoc_init_abs, ng_ext0_abs) + +### SPORE DISPERSION +### dispersion_kernel(disease, rda, alpha) +### return Disp, C_Disp +day_length = (24/(2*delta_t))*60*60 +alpha_asco = 3 # coefficient of dispersal (n m-1;Fitt et al 1987) +radius_asco = 5 +alpha_pycnid = 0.2 * 0.0001 # coefficient of dispersal (m² s-1; Yang et al 1991) +radius_pycnid = 5 +alpha_ure = 3 +radius_ure = 5 +kernel_asco, C_Disp_asco, kernel_pycnid, C_Disp_pycnid = dispersion_kernel_septo(day_length, alpha_asco, radius_asco, alpha_pycnid, radius_pycnid) +kernel_ure, C_Disp_ure = dispersion_kernel_rust(day_length, alpha_ure, radius_ure) + + +### SEIR FUNCTION +Nsp, Pth, Poi, Sth, Sus, Lat, Ifc, Ifv, Rem, LAI, LAI_wheat, Poo, Eps, AUDPC, Scont = \ + SEIR(t=t , delta_t0=delta_t0, delta_t=delta_t, season=season, + disease=disease, rain=rain, arrangement=arrangement, inoc_init=inoc_init, ng_ext0=ng_ext0, + lai_wheat_inc=lai_wheat_inc, lai_comp_inc=lai_comp_inc, lai_wheat_mu=lai_wheat_mu, + nu=nu, ber_wheat=ber_wheat, ber_companion=ber_companion, h_wheat=h_wheat, h_companion=h_companion, + lambd=lambd, delta_ei=delta_ei, + s0=s0,pi_inf0=pi_inf0, rho=rho, psi=psi, gamma=gamma, theta=theta, sigma=sigma, sigma_asco=sigma_asco, inf_begin=inf_begin, + C_Disp_ure=C_Disp_ure, kernel_ure=kernel_ure, C_Disp_asco=C_Disp_asco, kernel_asco=kernel_asco, + C_Disp_pycnid=C_Disp_pycnid, kernel_pycnid=kernel_pycnid) + + +### GENERATING CSV +# ## CSV file +T = [*range(0,t)] +Sth = np.mean(Sth, axis=(1,2)) +Pth = np.mean(Pth, axis=(1,2)) +Poi = np.mean(Poi, axis=(1,2)) +Lat = np.mean(Lat, axis=(1,2)) +Ifc = np.mean(Ifc, axis=(1,2)) +Ifv = np.mean(Ifv, axis=(1,2)) +Sus = np.mean(Sus, axis=(1,2)) +Rem = np.mean(Rem, axis=(1,2)) +Nsp = np.mean(Nsp, axis=(1,2)) +LAI = np.mean(LAI, axis=(1,2)) +LAI_wheat = np.mean(LAI_wheat, axis=(1,2)) +Poo = np.mean(Poo, axis=(1,2)) +AUDPC = np.mean(AUDPC, axis=(1,2)) +Eps = np.mean(Eps, axis=(1,2)) +Scont = np.mean(Scont, axis=(1,2)) + +plt.plot(T, Sth, color='black') +plt.plot(T, Sus, color='green') +plt.plot(T, Pth, color='brown') +plt.plot(T, Ifv, color='red') +plt.plot(T, Ifc, color='red',linestyle="--") +plt.show() + diff --git a/src/epymix/rain.py b/src/epymix/rain.py index a88e6f1..0d494bc 100644 --- a/src/epymix/rain.py +++ b/src/epymix/rain.py @@ -1,7 +1,7 @@ import numpy as np from .data import meteo_path -def rain(years, delta_t): +def rain(year, n_season, delta_t): """ Define rainfall scenario. @@ -20,11 +20,13 @@ def rain(years, delta_t): numpy array rain """ + years = np.arange(year,year+n_season,1) #1995: bad; 1997: medium, 2000: good + years = years.tolist() if not isinstance(years, list): years = [years] maxLa = 0 - _rain=[] + rain=[] for year in years: path = meteo_path(year) with open(path) as f: @@ -32,8 +34,8 @@ def rain(years, delta_t): a = np.array([int(v) for v in a]) A = -100 * np.ones((1000)) A[0:len(a)] = a - _rain = np.append(_rain,A) + rain = np.append(rain,A) maxLa = max(maxLa, len(a)) - _rain = _rain.reshape([len(A), len(years)], order='F') - _rain = np.ceil(_rain[0:maxLa]/int(delta_t)) - return _rain + rain = rain.reshape([len(A), len(years)], order='F') + rain = np.ceil(rain[0:maxLa]/int(delta_t)) + return rain