major: change file dir

Author: kongdd

Diff for: docs/src/index.md

@@ -7,6 +7,6 @@ cal_Ei_Dijk2021
 photosynthesis
 aerodynamic_conductance
 
-model_calib
+ModelCalib
 Param_PMLV2
 ```

Diff for: examples/CalibOneSite.jl

@@ -11,7 +11,7 @@ df.ET_obs = df.ETobs
 # ## 模型参数率定
 
 #nb # %% A slide [markdown] {"slideshow": {"slide_type": "fragment"}}
-theta, goal, flag = model_calib(df, par0)
+theta, goal, flag = ModelCalib(df, par0)
 df_out = PMLV2_sites(df; par=theta2par(theta))
 df_out[1:10, :]
 

Diff for: examples/ex2_NorthChina_C3C4_calib.jl

@@ -11,7 +11,7 @@ df = data[data.site.=="固城", :]
 par = par0
 r = PMLV2(df; par)
 
-theta, goal, flag = model_calib(df, par0)
+theta, goal, flag = ModelCalib(df, par0)
 df_out = PMLV2_sites(df; par=theta2par(theta))
 # df_out[1:10, :]
 sites = unique(df.site)

Diff for: examples/s1_calibrate_model_params.qmd

@@ -7,8 +7,8 @@ include("data_PML_2019.jl")
 # data.LAI .= data.LAI_raw
 # data.LAI .= data.LAI_whit
 # data = @rename(data, LAI = LAI_whit) # LAI_raw, LAI_whit, LAI_sgfitw
-# @time r = calib_PMLV2(data; of_gof=:KGE, maxn=5000)
-@time r = calib_PMLV2(data; of_gof=:NSE, maxn=10000)
+# @time r = ModelCalib_IGBPs(data; of_gof=:KGE, maxn=5000)
+@time r = ModelCalib_IGBPs(data; of_gof=:NSE, maxn=10000)
 r.gof
 ```
 

Diff for: src/Utilize/DataType.jl renamed to src/Optim/DataType.jl

@@ -1,23 +1,21 @@
-using Parameters
-
 @with_kw mutable struct interm_PML{T}
-  ET::T = T(0)
-  GPP::T = T(0)
-  Ec::T = T(0)
-  Ecr::T = T(0)
-  Eca::T = T(0)
-
-  Ei::T = T(0)
-  Pi::T = T(0)
-  Es_eq::T = T(0)
-  Eeq::T = T(0)
-  ET_water::T = T(0)
-
-  Ga::T = T(0)
-  Gc_w::T = T(0)
-
-  fval_soil::T = T(0)
-  Es::T = T(0)
+  ET::T = 0.0
+  GPP::T = 0.0
+  Ec::T = 0.0
+  Ecr::T = 0.0
+  Eca::T = 0.0
+
+  Ei::T = 0.0
+  Pi::T = 0.0
+  Es_eq::T = 0.0
+  Eeq::T = 0.0
+  ET_water::T = 0.0
+
+  Ga::T = 0.0
+  Gc_w::T = 0.0
+
+  fval_soil::T = 0.0
+  Es::T = 0.0
 end
 
 @with_kw mutable struct output_PML{T}
@@ -79,4 +77,3 @@ end
 
 export interm_PML, output_PML
 export to_mat;
-

Diff for: src/Utilize/calib_PMLV2.jl renamed to src/Optim/ModelCalib.jl

@@ -1,15 +1,23 @@
+export ModelCalib, ModelCalib_IGBPs
+export fread, fwrite, melt_list
+
+using RTableTools
 import Ipaper: par_map
+include("DataType.jl")
+include("ModelCalib.jl")
+include("model_gof.jl")
+# include("PMLV2_sites.jl")
+
 
-# all sites
-function calib_PMLV2(data::AbstractDataFrame; of_gof=:KGE, maxn=2500)
+function ModelCalib_IGBPs(data::AbstractDataFrame; of_gof=:KGE, maxn=2500)
   vars = ["IGBPname", "IGBPcode", "site", "date", "GPP_obs", "ET_obs",
     "Prcp", "Tavg", "U2", "Rn", "Rs", "VPD", "LAI", "Pa", "Ca"]
   IGBPs = unique(data.IGBPname) |> sort
 
   @time params = par_map(IGBP -> begin
       df = data[data.IGBP.==IGBP, vars]
       IGBPcode = df.IGBPcode[1]
-      theta, _, _ = model_calib(df, par0; IGBPcode, of_gof, maxn, verbose=false)
+      theta, _, _ = ModelCalib(df, par0; IGBPcode, of_gof, maxn, verbose=false)
       theta
     end, IGBPs; progress=false)
 
@@ -30,4 +38,18 @@ function calib_PMLV2(data::AbstractDataFrame; of_gof=:KGE, maxn=2500)
   (; param=theta2param(params, IGBPs), gof, output=df_out)
 end
 
-export calib_PMLV2
+
+## 一个站点的率定
+"""
+    ModelCalib(df::AbstractDataFrame, par0::AbstractETParam; 
+        IGBPcode=nothing, maxn=2500, of_gof=:NSE, kw...)
+"""
+function ModelCalib(df::AbstractDataFrame, par0::AbstractETParam; IGBPcode=nothing, maxn=2500, of_gof=:NSE, kw...)
+  parRanges = get_bounds(par0)
+  lower = parRanges[:, 1]
+  upper = parRanges[:, 2]
+  theta0 = collect(par0) # must be vector
+
+  sceua(theta -> -model_goal(df, theta; IGBPcode, of_gof, kw...),
+    theta0, lower, upper; maxn, kw...) # theta, goal, flag
+end

Diff for: src/Utilize/model_gof.jl renamed to src/Optim/model_gof.jl

@@ -9,6 +9,7 @@ function map_df_tuple(fun::Function, lst::GroupedDataFrame{DataFrame}, args...;
     end, 1:n)
 end
 
+
 ## 统计每种植被类型的模拟效果
 function model_gof(df_out::DataFrame; all=true)
   fun_et(d) = GOF(d.ET_obs, d.ET)
@@ -26,4 +27,30 @@ function model_gof(df_out::DataFrame; all=true)
   (; ET=gof_et, GPP=gof_gpp)
 end
 
+
+function model_goal(df, theta; IGBPcode=nothing, of_gof=:NSE, verbose=false)
+  # IGBPcode !== nothing && (par.hc = hc_raw[IGBPcode])
+  IGBPcode !== nothing && (theta[end] = hc_raw[IGBPcode]) # the last one is hc
+  par = Param_PMLV2(theta...)
+
+  dobs = df[!, [:GPP_obs, :ET_obs]]
+  # dsim = PMLV2(df; par) # the exact MATLAB version
+  dsim = PMLV2_sites(df; par) 
+
+  ## 8-day, yearly, yearly_anom
+  info_GPP = GOF(dobs.GPP_obs, dsim.GPP)
+  info_ET = GOF(dobs.ET_obs, dsim.ET)
+
+  if verbose
+    indexes = [:KGE, :NSE, :R2, :RMSE, :bias, :bias_perc]
+    info_GPP = info_GPP[indexes] |> round2
+    info_ET = info_ET[indexes] |> round2
+  end
+  gof = [info_ET[of_gof], info_GPP[of_gof]]
+  goal = weighted_mean(gof, [1, 0.5])
+  goal
+end
+
+
+export model_goal
 export model_gof, map_df_tuple

Diff for: src/PML.jl

@@ -18,7 +18,7 @@ file_FLUXNET_CRO_USTwt = "$dir_proj/data/CRO/FLUXNET_CRO_US-Twt" |> abspath
 
 include("main_Ipaper.jl")
 include("Params.jl")
-include("Utilize/Utilize.jl")
+include("Optim/ModelCalib.jl")
 include("ET_helper.jl")
 # include("PET_equilibrium.jl")
 # include("Ei_EvapIntercepted.jl")

Diff for: src/Utilize/Utilize.jl

@@ -30,10 +30,10 @@ end
 end
 
 
-@testset "model_calib" begin
+@testset "ModelCalib" begin
   df.GPP_obs = df.GPPobs
   df.ET_obs = df.ETobs
-  @time _theta, goal, flag = model_calib(df, par0; IGBPcode=df.IGBPcode[1], maxn=2500)
+  @time _theta, goal, flag = ModelCalib(df, par0; IGBPcode=df.IGBPcode[1], maxn=2500)
   goal = model_goal(df, _theta; verbose=true)
   @test goal > 0.55 # mean(KGE_GPP, KGE_ET)
 end