using DataFrames, CSV;

Import_data = "/user/start/examples/controls/parameter_estimation_webinar/parameter_estimation.csv" # @param {type:"filepicker"}
table = CSV.read(Import_data, DataFrame);
print(".导入csv文件")
Показать_таблицу = true # @param {type:"boolean"}
如果显示一个表
    display(table)
end  

Столбец_времени = 1 # @param {type:"integer"}
Столбец_входа = 2 # @param {type:"integer"}
Столбец_выхода = 3 # @param {type:"integer"}
Начинать_со_строки = 1 # @param {type:"integer"}
exp_time_data=table[Start with_strings:end,Time Column_]; 
exp_in_data=table[Start with_strings:end,Entry Column_];
exp_out_data=table[Start with_strings:end,Exit Column_];

df_in = DataFrame(time=exp_time_data, value=exp_in_data)
df_in = Float64.(df_in);
wa_in = WorkspaceArray(string(rand()), df_in)

model_name = "parameter_estimation" # @param {type:"string"}
if model_name in [m.name for m in engee.get_all_models()]
    m = engee.open( model_name ) # 加载模型
else
    m = engee.load( model_name * ".engee" )
end
logout = engee.run(m, verbose=false)
input = plot(exp_time_data, exp_in_data, title = "入口处", xlabel="t, c", label = "实验", legend=true)
output = plot(exp_time_data, exp_out_data, title = "出口；出口", xlabel="t, c", label = "实验", legend=true)
plot!(logout["w"].time, logout["w"].value, label = "模拟仿真", linewidth=2)
plot(input, output, layout = (2, 1))

using DataFrames
try
    using NLopt
catch
    Pkg.add(["NLopt"])
    using NLopt
end
try
    using Interpolations
catch
    Pkg.add(["Interpolations"])
    using Interpolations
end

# @markdown#模型参数的估计
# @markdown##要开始评估模型参数，请输入:
# @markdown###模型名称:
# 打开模型
model_name = "parameter_estimation" # @param {type:"string"}
if model_name in [m.name for m in engee.get_all_models()] # 检查将模型加载到内核的条件
    model = engee.open(model_name) # 打开模型
    model = engee.gcm() 
else
    @error "$Model_name模型未打开。"
end
# @markdown###登录（FromWorkspace类型块）:
input_name = "从工作区" # @param {type:"string"}
input_path = engee.find_system(model; depth=0, blockparams=["BlockType" => "FromWorkspace","BlockName"=>input_name])
# @markdown###输入值:
input_data = [exp_time_data exp_in_data] # @param {type:"raw"}
typeof(input_data) === Matrix{Float64} || throw(ArgumentError("输入信号必须是包含Float64数字的矩阵。"))
size(input_data,2) == 2 || throw(ArgumentError("输入信号应该有2列-时间和信号值。"))

isempty(engee.find_system(model; depth=0, blockparams=["BlockType" => "FromWorkspace","BlockName"=>input_name])) && throw(ArgumentError("Блока $input_name не существует в корневой системе модели $(model.name)."))
pe_udata_workspace = WorkspaceArray("pe_udata_workspace", DataFrame(time = input_data[:,1], value = input_data[:,2]));
engee.set_param!("$(model.name)/$input_name", "VariableName" =>"pe_udata_workspace") 

# @markdown###输出（可记录信号）:
output_name = "w" # @param {type:"string"}
output_name in keys(engee.run(model, verbose=false)) || throw(ArgumentError("Сигнал $output_name не найден среди логируемых сигналов."))

# @markdown###输出信号的值:
output_data = [exp_time_data exp_out_data] # @param {type:"raw"}
typeof(output_data) === Matrix{Float64} || throw(ArgumentError("输出信号必须是包含Float64数字的矩阵。"))
size(output_data,2) == 2 || throw(ArgumentError("输出信号应该有2列-时间和信号值。"))

# @markdown###被评估参数的名称:
# 创建变量、其初始值和约束的向量
p_names = ["b","c"] # @param {type:"raw"}
p_names isa Vector{String} || throw(ArgumentError("参数的名称必须作为字符串的向量输入。"))
parametrs_number = length(p_names)
parametrs_number>0 || throw(ArgumentError("参数的数量必须大于零。"))

# @markdown###初始参数值:
p_0 = [0.01, 0.7] # @param {type:"raw"}
p_0 isa Vector{Float64} || throw(ArgumentError("Parmeters的初始值必须是包含Float64数字的向量。"))
parametrs_number==length(p_0) || throw(ArgumentError("初始值和参数的数量必须匹配。"))

# @markdown###参数值下界:
p_lower_bounds = [0.0, 0.0] # @param {type:"raw"}
p_lower_bounds isa Vector{Float64} || throw(ArgumentError("参数值的下界应该是包含Float64数字的向量。"))
parametrs_number==length(p_lower_bounds) || throw(ArgumentError("下限值和parmeters的值的数量必须匹配。"))

# @markdown###参数值的上限:
p_upper_bounds = [1.0, 10.0] # @param {type:"raw"}
p_upper_bounds isa Vector{Float64} || throw(ArgumentError("参数值的上限应该是一个包含Float64数字的向量。"))
parametrs_number==length(p_0) || throw(ArgumentError("初始值和参数的数量必须匹配。"))

# @markdown##优化设置:
# @markdown###优化算法（算法名称开头的"G"表示是全局的，"L"是局部的）
algoritm = :G_MLSL # @param [:GN_DIRECT,:GN_DIRECT_L,:GN_DIRECT_L_RAND,:GN_CRS2_LM,:G_MLSL,:GN_AGS,:GN_ISRES,:GN_ESCH,:LN_COBYLA,:LN_BOBYQA,:LN_PRAXIS,:LN_NELDERMEAD,:LN_SBPLX] {type:"raw"}
# @markdown###辅助算法（仅适用于：G_MLSL）
local_algoritm = :LN_COBYLA # @param [:LN_COBYLA,:LN_BOBYQA,:LN_PRAXIS,:LN_NELDERMEAD,:LN_SBPLX] {type:"raw"}
# 配置优化参数
# @markdown###最小相对步长变化
xtol_rel = 1.0e-1 # @param {type:"number"}
xtol_rel>0 || throw(ArgumentError("最小相对阶跃变化应为正。"))
# @markdown###最小绝对步长变化
xtol_abs = 1.0e-1 # @param {type:"number"}
xtol_rel>0 || throw(ArgumentError("最小绝对阶跃变化应大于正值。"))
# @markdown###目标函数的最小相对变化
ftol_rel = 1.0e-1 # @param {type:"number"}
xtol_rel>0 || throw(ArgumentError("目标函数的最小相对变化应该是正的。"))
# @markdown###目标函数的最小绝对变化
ftol_abs = 1.0e-1 # @param {type:"number"}
xtol_rel>0 || throw(ArgumentError("目标函数的最小绝对变化应该是正的。"))
# @markdown###目标函数的最大估计数
maxeval = 30 # @param {type:"integer"}
xtol_rel>0 || throw(ArgumentError("目标函数的最大估计数应大于零。"))

# @markdown##显示参数评估结果
# @markdown###实测数据和仿真结果图:
compare = true # @param {type:"boolean"}
# @markdown###实测数据与仿真结果之间的误差图:
reside = true # @param {type:"boolean"}
# @markdown###参数估计值图:
variables = true # @param {type:"boolean"}
# @markdown###目标函数值图:
objective = true # @param {type:"boolean"}

# 比较模型输出和测量数据的图形
function compare_plot(model, output_name::String, tdata::Vector{Float64}, udata::Vector{Float64}, ydata::Vector{Float64})
    logout = engee.run(model, verbose=false)
    out_plot = plot(
        logout[output_name].time, 
        logout[output_name].value, 
        label="模拟仿真", 
        title = "出口；出口", 
        xlabel="t，与", 
        legend=:bottomright,
        linewidth=3
        )
    plot!(out_plot, tdata, ydata, label="实验")
    in_plot = plot(tdata, udata, title = "入口处", xlabel="t, c", legend = false)
    p = plot(out_plot,in_plot,layout = (2, 1), plot_title="实验及模拟")
end

# 差异附表
function residual_plot(model, output_name::String, tdata::Vector{Float64}, ydata::Vector{Float64})
    logout = engee.run(model, verbose=false)
    nearest = interpolate((logout[output_name].time,), logout[output_name].value, Gridded(Constant()))
    residals = ydata .- nearest.(tdata)
    reside = plot(tdata, residals, title = "错误", xlabel="t，与", legend = false)
end

# 变量更改时间表
function variable_plot(p_names::Vector{String}, parametrs_values)
    p_trace = plot(title = "参数值", xlabel="目标函数的估计数",)
    parametrs_values_m = mapreduce(permutedims, vcat, parametrs_values)
    for i in 1:n
        plot!(p_trace, 1:length(parametrs_values_m[:,i]), parametrs_values_m[:,i], markers=:auto, label = p_names[i])
    end
    return p_trace
end

n = parametrs_number

# 用于记录目标函数和变量值的向量
loss_function_vector = Vector{Float64}()
parametrs_values = Vector{Vector{Float64}}()

# 损失函数
function myfunc(x::Vector, grad::Vector, output_name::String, measure_time, measure_value, parametrs_name...)
    str = ""
    iterations_value = Vector{Float64}()
    for (value,p) in zip(x, parametrs_name)
        global parametr = Symbol(p)
        @eval(($parametr)=$(value))
        str = str*"$p = $value, "
        push!(iterations_value, value)
    end
    push!(parametrs_values, iterations_value)
    simulation_out = engee.run(model, verbose=false) |> x->collect(x[output_name])
    if measure_time[end]>simulation_out.time[end]  
        push!(simulation_out,[measure_time[end], simulation_out.time[end]])
    end
    nearest = interpolate(
        (simulation_out.time,), 
        simulation_out.value, 
        Gridded(Constant())
    )
    l = sum(abs2, measure_value .- nearest.(measure_time))
    str = "$(length(parametrs_values)). "*str*"loss = $l"
    println(str)
    push!(loss_function_vector, l)
    return l
end


opt = Opt(algoritm, n) # 无梯度算法
algoritm === :G_MLSL && NLopt.local_optimizer!(opt, Opt(local_algoritm, n))
opt.lower_bounds = p_lower_bounds
opt.upper_bounds = p_upper_bounds
opt.xtol_rel = xtol_rel
opt.xtol_abs = xtol_abs
opt.ftol_rel = ftol_rel
opt.ftol_abs = ftol_rel
opt.maxeval =  maxeval

opt.min_objective = (x,g) -> myfunc(x,g,output_name,output_data[:,1],output_data[:,2],p_names...)

(minf,minx,ret) = optimize(opt,p_0)
numevals = opt.numevals 

str = "\评价结果：\n"
for (name, value) in zip(p_names, minx) 
    str = str*"$name = $value\n"
end
str = str*"\N目标函数值=$minf。\n\优化在目标函数的估计值的num numevals之后完成，状态为：$ret。\n"
print(str)

compare && display(compare_plot(model, output_name, input_data[:,1], input_data[:,2], output_data[:,2]))
println("")
reside && display(residual_plot(model, output_name, output_data[:,1], output_data[:,2]))
println("")
objective &&
    display(
        plot(1:numevals, loss_function_vector, title = "目标功能", xlabel="目标函数的估计数", legend = false, markers = :circle)
    )
println("")   
variables && display(variable_plot(p_names, parametrs_values))
println("")

print("评价结果：b= ", b, ", c = ", c)