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 file imported")
Показать_таблицу = true # @param {type:"boolean"}
if to show a table
    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 ) # loading the model
else
    m = engee.load( model_name * ".engee" )
end
logout = engee.run(m, verbose=false)
input = plot(exp_time_data, exp_in_data, title = "Entrance", xlabel="t, c", label = "experiment", legend=true)
output = plot(exp_time_data, exp_out_data, title = "Exit", xlabel="t, c", label = "experiment", legend=true)
plot!(logout["w"].time, logout["w"].value, label = "simulation", 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 # Estimation of model parameters
# @markdown ## To start evaluating the model parameters, enter:
# @markdown ### Model name:
# Opening the model
model_name = "parameter_estimation" # @param {type:"string"}
if model_name in [m.name for m in engee.get_all_models()] # Checking the condition for loading a model into the kernel
    model = engee.open(model_name) # Open the model
    model = engee.gcm() 
else
    @error "The $model_name model is not open."
end
# @markdown ### Login (FromWorkspace type block):
input_name = "From the workspace" # @param {type:"string"}
input_path = engee.find_system(model; depth=0, blockparams=["BlockType" => "FromWorkspace","BlockName"=>input_name])
# @markdown ### Input value:
input_data = [exp_time_data exp_in_data] # @param {type:"raw"}
typeof(input_data) === Matrix{Float64} || throw(ArgumentError("The input signal must be a matrix containing Float64 numbers."))
size(input_data,2) == 2 || throw(ArgumentError("The input signal should have 2 columns - time and signal values."))

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 (logable signal):
output_name = "w" # @param {type:"string"}
output_name in keys(engee.run(model, verbose=false)) || throw(ArgumentError("Сигнал $output_name не найден среди логируемых сигналов."))

# @markdown ### The value of the output signal:
output_data = [exp_time_data exp_out_data] # @param {type:"raw"}
typeof(output_data) === Matrix{Float64} || throw(ArgumentError("The output signal must be a matrix containing Float64 numbers."))
size(output_data,2) == 2 || throw(ArgumentError("The output signal should have 2 columns - time and signal values."))

# @markdown ### Names of evaluated parameters:
# Creating a vector of variables, their initial values and constraints
p_names = ["b","c"] # @param {type:"raw"}
p_names isa Vector{String} || throw(ArgumentError("The names of the parameters must be entered as a vector of strings."))
parametrs_number = length(p_names)
parametrs_number>0 || throw(ArgumentError("The number of parameters must be greater than zero."))

# @markdown ### Initial parameter values:
p_0 = [0.01, 0.7] # @param {type:"raw"}
p_0 isa Vector{Float64} || throw(ArgumentError("The initial values of the parmeters must be a vector containing Float64 numbers."))
parametrs_number==length(p_0) || throw(ArgumentError("The number of initial values and parmeters must match."))

# @markdown ### Lower bound of parameter values:
p_lower_bounds = [0.0, 0.0] # @param {type:"raw"}
p_lower_bounds isa Vector{Float64} || throw(ArgumentError("The lower bound of the parameter values should be a vector containing Float64 numbers."))
parametrs_number==length(p_lower_bounds) || throw(ArgumentError("The number of values of the lower limit and the parmeters must match."))

# @markdown ### Upper bound of parameter values:
p_upper_bounds = [1.0, 10.0] # @param {type:"raw"}
p_upper_bounds isa Vector{Float64} || throw(ArgumentError("The upper bound of the parameter values should be a vector containing Float64 numbers."))
parametrs_number==length(p_0) || throw(ArgumentError("The number of initial values and parmeters must match."))

# @markdown ## Optimization Settings:
# @markdown ### Optimization algorithm (the "G" at the beginning of the algorithm name indicates that it is global, the "L" is local)
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 ### Auxiliary algorithm (only for :G_MLSL)
local_algoritm = :LN_COBYLA # @param [:LN_COBYLA,:LN_BOBYQA,:LN_PRAXIS,:LN_NELDERMEAD,:LN_SBPLX] {type:"raw"}
# Configuring optimization parameters
# @markdown ### Minimum relative step change
xtol_rel = 1.0e-1 # @param {type:"number"}
xtol_rel>0 || throw(ArgumentError("The minimum relative step change should be positive."))
# @markdown ### Minimum absolute step change
xtol_abs = 1.0e-1 # @param {type:"number"}
xtol_rel>0 || throw(ArgumentError("The minimum absolute step change should be more than positive."))
# @markdown ### Minimal relative change in objective function
ftol_rel = 1.0e-1 # @param {type:"number"}
xtol_rel>0 || throw(ArgumentError("The minimum relative change in the objective function should be positive."))
# @markdown ### Minimum absolute change in objective function
ftol_abs = 1.0e-1 # @param {type:"number"}
xtol_rel>0 || throw(ArgumentError("The minimum absolute change in the objective function should be positive."))
# @markdown ### Maximum number of estimates of the objective function
maxeval = 30 # @param {type:"integer"}
xtol_rel>0 || throw(ArgumentError("The maximum number of estimates of the objective function should be greater than zero."))

# @markdown ## Displaying parameter evaluation results
# @markdown ### Graph of measured data and simulation results:
compare = true # @param {type:"boolean"}
# @markdown ### Error graph between measured data and simulation results:
reside = true # @param {type:"boolean"}
# @markdown ### Graph of parameter estimation values:
variables = true # @param {type:"boolean"}
# @markdown ### Graph of objective function values:
objective = true # @param {type:"boolean"}

# Graph comparing the output of the model and the measured data
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="simulation", 
        title = "Exit", 
        xlabel="t, with", 
        legend=:bottomright,
        linewidth=3
        )
    plot!(out_plot, tdata, ydata, label="experiment")
    in_plot = plot(tdata, udata, title = "Entrance", xlabel="t, c", legend = false)
    p = plot(out_plot,in_plot,layout = (2, 1), plot_title="Experiment and simulation")
end

# schedule of discrepancies
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 = "Mistakes", xlabel="t, with", legend = false)
end

# Schedule of change of variables
function variable_plot(p_names::Vector{String}, parametrs_values)
    p_trace = plot(title = "Parameter values", xlabel="Number of estimates of the objective function",)
    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

# Vectors for logging the values of the objective function and variables
loss_function_vector = Vector{Float64}()
parametrs_values = Vector{Vector{Float64}}()

# loss function
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) # a gradient-free algorithm
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 = "\Evaluation results:\n"
for (name, value) in zip(p_names, minx) 
    str = str*"$name = $value\n"
end
str = str*"\n Target function value = $minf.\n\Optimization is completed after $numevals of estimates of the target function with the status: $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 = "Target function", xlabel="Number of estimates of the objective function", legend = false, markers = :circle)
    )
println("")   
variables && display(variable_plot(p_names, parametrs_values))
println("")

print("Evaluation results: b = ", b, ", c = ", c)