改进参数优化方案

src/Optim/DataType.jl → src/DataType.jl

@@ -75,6 +75,18 @@ function to_df(res::output_PML{T}) where {T<:Real}
   DataFrame(data, names)
 end
 
+function map_df_tuple(fun::Function, lst::GroupedDataFrame{DataFrame}, args...; kw...)
+  n = length(lst)
+  _keys = keys(lst)
+  map(i -> begin
+      d = lst[i]
+      key = NamedTuple(_keys[i])
+      r = fun(d, args...; kw...)
+      (; key..., r...)
+    end, 1:n)
+end
+
 
 export interm_PML, output_PML
 export to_mat;
+export map_df_tuple

src/ModelCalib.jl

@@ -0,0 +1,95 @@
+export ModelCalib, ModelCalib_IGBPs
+export model_goal
+export model_gof
+
+include("DataType.jl")
+# include("model_gof.jl")
+
+function ModelCalib_IGBPs(data::AbstractDataFrame;
+  parNames=ParNames,
+  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, _, _ = ModelCalib(df, par0, parNames; IGBPcode, of_gof, maxn, verbose=false)
+      theta
+    end, IGBPs; progress=false)
+
+  # printstyled("[i=$i, IGBP = $IGBP] \n", bold=true, color=:green, underline=false)
+  res = map(i -> begin
+      IGBP = IGBPs[i]
+      theta = params[i]
+      df = data[data.IGBP.==IGBP, vars]
+      model_goal(df, theta, parNames; verbose=true) # print GOF
+
+      par = theta2par(theta, parNames)
+      r = PMLV2_sites(df; par)
+      cbind(df[:, [:site, :date, :ET_obs, :GPP_obs]], r)
+    end, eachindex(params))
+  df_out = melt_list(res; IGBP=IGBPs)
+
+  gof = model_gof(df_out)
+  (; param=theta2param(params, IGBPs), gof, output=df_out)
+end
+
+
+## 一个站点的率定
+"""
+    ModelCalib(df::AbstractDataFrame, par0::AbstractETParam; 
+        IGBPcode=nothing, maxn=2500, of_gof=:NSE, kw...)
+"""
+function ModelCalib(df::AbstractDataFrame, par0::AbstractETParam, parNames;
+  IGBPcode=nothing, maxn=2500, of_gof=:NSE, kw...)
+  lower, upper = get_bounds(parNames)
+  x0 = select_param(par0, parNames)
+
+  sceua(theta -> -model_goal(df, theta, parNames; IGBPcode, of_gof, kw...),
+    x0, lower, upper; maxn, kw...) # theta, goal, flag
+end
+
+
+## gof of ET and GPP for one IGBP
+# cal df_out by Forcing and Param
+function model_goal(df, theta, parNames; IGBPcode=nothing, of_gof=:NSE, verbose=false)
+  par = theta2par(theta, parNames)
+  :hc ∉ parNames && !isnothing(IGBPcode) && (par.hc = hc_raw[IGBPcode])
+
+  dobs = df[!, [:GPP_obs, :ET_obs]]
+  dsim = PMLV2_sites(df; par)
+  # dsim = PMLV2(df; par) # the exact MATLAB version
+
+  ## 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]]
+  weighted_mean(gof, [1.0, 1.0]) # goal
+end
+
+
+## gof of ET and GPP for ALL IGBP
+function model_gof(df_out::DataFrame; all=true)
+  fun_et(d) = GOF(d.ET_obs, d.ET)
+  fun_gpp(d) = GOF(d.GPP_obs, d.GPP)
+
+  lst = groupby(df_out, :IGBP)
+
+  res = map_df_tuple(fun_et, lst)
+  all && push!(res, (; IGBP="ALL", fun_et(df_out)...))
+  gof_et = DataFrame(res)
+
+  res = map_df_tuple(fun_gpp, lst)
+  all && push!(res, (; IGBP="ALL", fun_gpp(df_out)...))
+  gof_gpp = DataFrame(res)
+  (; ET=gof_et, GPP=gof_gpp)
+end

src/PML.jl

@@ -3,16 +3,22 @@ module PML
 export PMLV2, PMLV2_sites, 
   photosynthesis, cal_Ei_Dijk2021, 
   T_adjust_Vm25, f_VPD_Zhang2019
-export DataFrame, GOF
+
 export file_FLUXNET_CRO, file_FLUXNET_CRO_USTwt
+export DataFrame, GOF
+export fread, fwrite, melt_list
+
 
 using Parameters
 using FieldMetadata
 using DataFrames
 using Statistics
-import HydroTools: cal_Uz, Cp, atm, GOF, sceua
+using RTableTools
 using DocStringExtensions
 
+import HydroTools: cal_Uz, Cp, atm, GOF, sceua
+import Ipaper: par_map
+
 ## global data
 dir_proj = "$(@__DIR__)/.."
 file_FLUXNET_CRO = "$dir_proj/data/CRO/FLUXNET_CRO" |> abspath
@@ -21,18 +27,16 @@ file_FLUXNET_CRO_USTwt = "$dir_proj/data/CRO/FLUXNET_CRO_US-Twt" |> abspath
 
 include("main_Ipaper.jl")
 include("Params.jl")
-include("Optim/ModelCalib.jl")
+include("ModelCalib.jl")
 include("ET_helper.jl")
+include("water_constrain.jl")
+include("photosynthesis.jl")
 # include("PET_equilibrium.jl")
 # include("Ei_EvapIntercepted.jl")
 # include("Ec_CanopyTrans.jl")
 # include("Es_EvapSoil.jl")
-include("water_constrain.jl")
-include("photosynthesis.jl")
 # include("model_PMLV2.jl")
 
-
-
 """
     PMLV2 (Penman–Monteith–Leuning Version 2) Evapotranspiration model
 

src/Params.jl

@@ -58,26 +58,30 @@ function Base.collect(par::AbstractETParam)
   [getfield(par, f) for f in fieldnames(typeof(par))]
 end
 
-function get_bounds(par::AbstractETParam, parNames::Vector{Symbol}=ParNames)
+function get_bounds(parNames::Vector{Symbol}=ParNames)
   inds = match2(parNames, ParNames).I_y # 选择的参数
-  bs = bounds(par)[inds]
-  x0 = vcat([getfield(par, f) for f in parNames]...)
+  bs = bounds(Param_PMLV2)[inds]
   lower = vcat(map(x -> x[1], bs)...)
   upper = vcat(map(x -> x[2], bs)...)
-  x0, lower, upper
+  lower, upper
 end
 
 
 # theta2par(theta) = Param_PMLV2(theta...)
-function theta2par(theta::Vector, par::AbstractETParam, parNames::Vector{Symbol})
+function theta2par!(theta::Vector, par::AbstractETParam, parNames::Vector{Symbol})
   k = 0
   for key in parNames
     x = getfield(par, key)
     n = length(x)
-    inds = k+1:k+n
+    inds = n == 1 ? k+1 : k+1:k+n
     k += n
     setfield!(par, key, theta[inds])
   end
+  return par
+end
+
+function theta2par(theta::Vector, parNames::Vector{Symbol})
+  theta2par!(theta, deepcopy(par0), parNames)
 end
 
 function select_param(par::AbstractETParam, parNames)
@@ -88,11 +92,8 @@ end
 hc_raw = [10, 10, 10, 10, 10, 1, 1, 5, 5, 0.2, 1, 0.5, 10, 1, 0.01, 0.05, 0.01]
 
 FT = Float64
-
 par0 = Param_PMLV2()
 theta0 = collect(par0)
-# parNames = fieldnames(Param_PMLV2)
-# parRanges = get_bounds(par0)
 
 function theta2param(params::Vector{Vector{T}}, IGBPs) where {T<:Real}
   parNames = fieldnames(Param_PMLV2) |> collect
@@ -103,7 +104,7 @@ function theta2param(params::Vector{Vector{T}}, IGBPs) where {T<:Real}
 end
 
 
-export bounds, Param_PMLV2
-export parRanges, parNames, theta0, par0, hc_raw
-export get_bounds, select_param
+export bounds, get_bounds, select_param
+export Param_PMLV2, ParNames
+export theta0, par0, hc_raw
 export param_PML, theta2par, theta2param

test/test-PMLV2.jl

@@ -1,8 +1,9 @@
 using PML, Test, Ipaper, DataFrames
 
 df_out, df, _par = deserialize(file_FLUXNET_CRO_USTwt)
-par = Param_PMLV2(;_par..., hc=0.5)
-
+df.GPP_obs = df.GPPobs
+df.ET_obs = df.ETobs
+par = Param_PMLV2(; _par..., hc=0.5)
 r = PMLV2_sites(df; par)
 
 @testset "PMLV2 scalar" begin
@@ -32,18 +33,21 @@ end
 end
 
 
-# @testset "ModelCalib" begin
-#   df.GPP_obs = df.GPPobs
-#   df.ET_obs = df.ETobs
-#   @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
-
-# @testset "model_gof" begin
-#   d = DataFrame(; IGBP=["CRO", "CRO", "CRO"],
-#     ET_obs=[1.0, 2, 3], ET=[0.5, 0.6, 0.7],
-#     GPP_obs=[1.0, 2, 3], GPP=[1, 1.5, 2.0])
-#   @test size(model_gof(d).ET, 1) == 2
-#   @test size(model_gof(d; all=false).ET, 1) == 1
-# end
+@testset "model_gof" begin
+  d = DataFrame(; IGBP=["CRO", "CRO", "CRO"],
+    ET_obs=[1.0, 2, 3], ET=[0.5, 0.6, 0.7],
+    GPP_obs=[1.0, 2, 3], GPP=[1, 1.5, 2.0])
+  @test size(model_gof(d).ET, 1) == 2
+  @test size(model_gof(d; all=false).ET, 1) == 1
+end
+
+
+@testset "ModelCalib" begin
+  parNames = [
+    :α, :η, :g1, :Am_25, :VPDmin, :VPDmax, :D0, :kQ, :kA, :S_sls, :fER0 # :hc
+  ]
+
+  @time _theta, goal, flag = ModelCalib(df, par0, parNames; IGBPcode=df.IGBPcode[1], maxn=2500)
+  goal = model_goal(df, _theta, parNames; verbose=true)
+  @test goal > 0.55 # mean(KGE_GPP, KGE_ET)
+end