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"}
if Показать_таблицу
    display(table)
end  

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

model_name = "parameter_estimation" # @param {type:"string"}
input_name = "Из рабочей области" # @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

df_in = DataFrame(time=exp_time_data, value=exp_in_data)
df_in = Float64.(df_in);
wa_in = WorkspaceArray(string(rand()), df_in)
engee.set_param!("$(model_name)/$input_name", "VariableName" =>"wa_in") 


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("Начальные значения парметров должны быть вектором, содержащим числа типа 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("Количество значений нижней границы и парметров должно совпадать."))

# @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Результаты оценки:\n"
for (name, value) in zip(p_names, minx) 
    str = str*"$name = $value\n"
end
str = str*"\nЗначение целевой функции = $minf.\n\nОптимизация завершена после $numevals оценок целевой функции cо статусом: $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)