text:
include("KinematicModel.jl")

function vars_state(m::KinematicModel, FT)
    @vars begin
        ρ::FT
        ρu::SVector{3, FT}
        ρe::FT
        ρq_tot::FT
        ρq_liq::FT
        ρq_ice::FT
        ρq_rai::FT
    end
end

function vars_aux(m::KinematicModel, FT)
    @vars begin
        # defined in init_aux
        p::FT
        z::FT
        # defined in update_aux
        u::FT
        w::FT
        q_tot::FT
        q_vap::FT
        q_liq::FT
        q_ice::FT
        q_rai::FT
        e_tot::FT
        e_kin::FT
        e_pot::FT
        e_int::FT
        T::FT
        S::FT
        RH::FT
        rain_w::FT
        # uncomment below for more diagnostics
        #src_cloud_liq::FT
        #src_cloud_ice::FT
        #src_acnv::FT
        #src_accr::FT
        #src_rain_evap::FT
        #flag_rain::FT
        #flag_cloud_liq::FT
        #flag_cloud_ice::FT
    end
end

function init_kinematic_eddy!(eddy_model, state, aux, (x, y, z), t)
    FT = eltype(state)
    _grav::FT = grav(param_set)

    dc = eddy_model.data_config

    # density
    q_pt_0 = PhasePartition(dc.qt_0)
    R_m, cp_m, cv_m, γ = gas_constants(q_pt_0)
    T::FT = dc.θ_0 * (aux.p / dc.p_1000)^(R_m / cp_m)
    ρ::FT = aux.p / R_m / T
    state.ρ = ρ

    # moisture
    state.ρq_tot = ρ * dc.qt_0
    state.ρq_liq = ρ * q_pt_0.liq
    state.ρq_ice = ρ * q_pt_0.ice
    state.ρq_rai = ρ * FT(0)

    # velocity (derivative of streamfunction)
    ρu::FT =
        dc.wmax * dc.xmax / dc.zmax *
        cos(π * z / dc.zmax) *
        cos(2 * π * x / dc.xmax)
    ρw::FT = 2 * dc.wmax * sin(π * z / dc.zmax) * sin(2 * π * x / dc.xmax)
    state.ρu = SVector(ρu, FT(0), ρw)
    u::FT = ρu / ρ
    w::FT = ρw / ρ

    # energy
    e_kin::FT = 1 // 2 * (u^2 + w^2)
    e_pot::FT = _grav * z
    e_int::FT = internal_energy(T, q_pt_0)
    e_tot::FT = e_kin + e_pot + e_int
    state.ρe = ρ * e_tot

    return nothing
end

function kinematic_model_nodal_update_aux!(
    m::KinematicModel,
    state::Vars,
    aux::Vars,
    t::Real,
)
    FT = eltype(state)
    _grav::FT = grav(param_set)
    # velocity
    aux.u = state.ρu[1] / state.ρ
    aux.w = state.ρu[3] / state.ρ
    # water
    aux.q_tot = state.ρq_tot / state.ρ
    aux.q_liq = state.ρq_liq / state.ρ
    aux.q_ice = state.ρq_ice / state.ρ
    aux.q_rai = state.ρq_rai / state.ρ
    aux.q_vap = aux.q_tot - aux.q_liq - aux.q_ice
    # energy
    aux.e_tot = state.ρe / state.ρ
    aux.e_kin = 1 // 2 * (aux.u^2 + aux.w^2)
    aux.e_pot = _grav * aux.z
    aux.e_int = aux.e_tot - aux.e_kin - aux.e_pot
    # supersaturation
    q = PhasePartition(aux.q_tot, aux.q_liq, aux.q_ice)
    aux.T = air_temperature(aux.e_int, q)
    aux.S =
        max(0, aux.q_vap / q_vap_saturation(aux.T, state.ρ, q) - FT(1)) *
        FT(100)
    aux.RH = aux.q_vap / q_vap_saturation(aux.T, state.ρ, q) * FT(100)

    aux.rain_w = terminal_velocity(aux.q_rai, state.ρ)

    # uncomment below for more diagnostics
    #q_eq = PhasePartition_equil(aux.T, state.ρ, aux.q_tot)
    #aux.src_cloud_liq = conv_q_vap_to_q_liq(q_eq, q)
    #aux.src_cloud_ice = conv_q_vap_to_q_ice(q_eq, q)
    #aux.src_acnv = conv_q_liq_to_q_rai_acnv(aux.q_liq)
    #aux.src_accr = conv_q_liq_to_q_rai_accr(aux.q_liq, aux.q_rai, state.ρ)
    #aux.src_rain_evap = conv_q_rai_to_q_vap(aux.q_rai, q, aux.T, aux.p, state.ρ)
    #aux.flag_cloud_liq = FT(0)
    #aux.flag_cloud_ice = FT(0)
    #aux.flag_rain = FT(0)
    #if (aux.q_liq >= FT(0))
    #    aux.flag_cloud_liq = FT(1)
    #end
    #if (aux.q_ice >= FT(0))
    #    aux.flag_cloud_ice = FT(1)
    #end
    #if (aux.q_rai >= FT(0))
    #    aux.flag_rain = FT(1)
    #end
end

function boundary_state!(
    ::Rusanov,
    m::KinematicModel,
    state⁺,
    aux⁺,
    n,
    state⁻,
    aux⁻,
    bctype,
    t,
    args...,
)
    #state⁺.ρu -= 2 * dot(state⁻.ρu, n) .* SVector(n)
    state⁺.ρq_rai = -state⁻.ρq_rai
end

@inline function wavespeed(
    m::KinematicModel,
    nM,
    state::Vars,
    aux::Vars,
    t::Real,
)
    FT = eltype(state)
    u = state.ρu / state.ρ
    rain_w = terminal_velocity(state.ρq_rai / state.ρ, state.ρ)
    nu = nM[1] * u[1] + nM[3] * max(u[3], rain_w, u[3] - rain_w)

    return abs(nu)
end

@inline function flux_nondiffusive!(
    m::KinematicModel,
    flux::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    FT = eltype(state)
    rain_w = terminal_velocity(state.ρq_rai / state.ρ, state.ρ)

    # advect moisture ...
    flux.ρq_tot = SVector(
        state.ρu[1] * state.ρq_tot / state.ρ,
        FT(0),
        state.ρu[3] * state.ρq_tot / state.ρ,
    )
    flux.ρq_liq = SVector(
        state.ρu[1] * state.ρq_liq / state.ρ,
        FT(0),
        state.ρu[3] * state.ρq_liq / state.ρ,
    )
    flux.ρq_ice = SVector(
        state.ρu[1] * state.ρq_ice / state.ρ,
        FT(0),
        state.ρu[3] * state.ρq_ice / state.ρ,
    )
    flux.ρq_rai = SVector(
        state.ρu[1] * state.ρq_rai / state.ρ,
        FT(0),
        (state.ρu[3] / state.ρ - rain_w) * state.ρq_rai,
    )
    # ... energy ...
    flux.ρe = SVector(
        state.ρu[1] / state.ρ * (state.ρe + aux.p),
        FT(0),
        state.ρu[3] / state.ρ * (state.ρe + aux.p),
    )
    # ... and don't advect momentum (kinematic setup)
end

function source!(
    m::KinematicModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
)
    # TODO - ensure positive definite
    FT = eltype(state)
    _grav::FT = grav(param_set)
    _e_int_v0::FT = e_int_v0(param_set)
    _cv_v::FT = cv_v(param_set)
    _cv_d::FT = cv_d(param_set)
    _T_0::FT = T_0(param_set)

    e_tot = state.ρe / state.ρ
    q_tot = state.ρq_tot / state.ρ
    q_liq = state.ρq_liq / state.ρ
    q_ice = state.ρq_ice / state.ρ
    q_rai = state.ρq_rai / state.ρ
    u = state.ρu[1] / state.ρ
    w = state.ρu[3] / state.ρ
    e_int = e_tot - 1 // 2 * (u^2 + w^2) - _grav * aux.z

    q = PhasePartition(q_tot, q_liq, q_ice)
    T = air_temperature(e_int, q)
    # equilibrium state at current T
    q_eq = PhasePartition_equil(T, state.ρ, q_tot)

    # zero out the source terms
    source.ρq_tot = FT(0)
    source.ρq_liq = FT(0)
    source.ρq_ice = FT(0)
    source.ρq_rai = FT(0)
    source.ρe = FT(0)

    # cloud water and ice condensation/evaporation
    source.ρq_liq += state.ρ * conv_q_vap_to_q_liq(q_eq, q)
    source.ρq_ice += state.ρ * conv_q_vap_to_q_ice(q_eq, q)

    # tendencies from rain
    src_q_rai_acnv = conv_q_liq_to_q_rai_acnv(q_liq)
    src_q_rai_accr = conv_q_liq_to_q_rai_accr(q_liq, q_rai, state.ρ)
    src_q_rai_evap = conv_q_rai_to_q_vap(q_rai, q, T, aux.p, state.ρ)

    src_q_rai_tot = src_q_rai_acnv + src_q_rai_accr + src_q_rai_evap

    source.ρq_liq -= state.ρ * (src_q_rai_acnv + src_q_rai_accr)
    source.ρq_rai += state.ρ * src_q_rai_tot
    source.ρq_tot -= state.ρ * src_q_rai_tot
    source.ρe -=
        state.ρ * src_q_rai_tot * (_e_int_v0 - (_cv_v - _cv_d) * (T - _T_0))
end

function main()
    CLIMA.init()

    # Working precision
    FT = Float64
    # DG polynomial order
    N = 4
    # Domain resolution and size
    Δx = FT(20)
    Δy = FT(1)
    Δz = FT(20)
    resolution = (Δx, Δy, Δz)
    # Domain extents
    xmax = 1500
    ymax = 10
    zmax = 1500
    # initial configuration
    wmax = FT(0.6)  # max velocity of the eddy  [m/s]
    θ_0 = FT(289) # init. theta value (const) [K]
    p_0 = FT(101500) # surface pressure [Pa]
    p_1000 = FT(100000) # reference pressure in theta definition [Pa]
    qt_0 = FT(7.5 * 1e-3) # init. total water specific humidity (const) [kg/kg]
    z_0 = FT(0) # surface height

    # time stepping
    t_ini = FT(0)
    t_end = FT(30 * 60)
    dt = FT(5)
    #CFL = FT(1.75)
    filter_freq = 1
    output_freq = 72

    driver_config = config_kinematic_eddy(
        FT,
        N,
        resolution,
        xmax,
        ymax,
        zmax,
        wmax,
        θ_0,
        p_0,
        p_1000,
        qt_0,
        z_0,
    )
    solver_config = CLIMA.SolverConfiguration(
        t_ini,
        t_end,
        driver_config;
        ode_dt = dt,
        init_on_cpu = true,
        #Courant_number = CFL,
    )

    model = driver_config.bl

    mpicomm = MPI.COMM_WORLD

    # get state variables indices for filtering
    ρq_liq_ind = varsindex(vars_state(model, FT), :ρq_liq)
    ρq_ice_ind = varsindex(vars_state(model, FT), :ρq_ice)
    ρq_rai_ind = varsindex(vars_state(model, FT), :ρq_rai)
    # get aux variables indices for testing
    q_tot_ind = varsindex(vars_aux(model, FT), :q_tot)
    q_vap_ind = varsindex(vars_aux(model, FT), :q_vap)
    q_liq_ind = varsindex(vars_aux(model, FT), :q_liq)
    q_ice_ind = varsindex(vars_aux(model, FT), :q_ice)
    q_rai_ind = varsindex(vars_aux(model, FT), :q_rai)
    S_ind = varsindex(vars_aux(model, FT), :S)
    rain_w_ind = varsindex(vars_aux(model, FT), :rain_w)

    # filetr out negative values
    cb_tmar_filter =
        GenericCallbacks.EveryXSimulationSteps(filter_freq) do (init = false)
            Filters.apply!(
                solver_config.Q,
                (ρq_liq_ind[1], ρq_ice_ind[1], ρq_rai_ind[1]),
                solver_config.dg.grid,
                TMARFilter(),
            )
            nothing
        end

    # output for paraview
    step = [0]
    cb_vtk =
        GenericCallbacks.EveryXSimulationSteps(output_freq) do (init = false)
            mkpath("vtk/")
            outprefix = @sprintf(
                "vtk/new_ex_2_mpirank%04d_step%04d",
                MPI.Comm_rank(mpicomm),
                step[1]
            )
            @info "doing VTK output" outprefix
            writevtk(
                outprefix,
                solver_config.Q,
                solver_config.dg,
                flattenednames(vars_state(model, FT)),
                solver_config.dg.auxstate,
                flattenednames(vars_aux(model, FT)),
            )
            step[1] += 1
            nothing
        end

    # call solve! function for time-integrator
    result = CLIMA.invoke!(
        solver_config;
        user_callbacks = (cb_tmar_filter, cb_vtk),
        check_euclidean_distance = true,
    )

    # supersaturation in the model
    max_S = maximum(abs.(solver_config.dg.auxstate[:, S_ind, :]))
    @test max_S < FT(0.25)
    @test max_S > FT(0)

    # qt < reference number
    max_q_tot = maximum(abs.(solver_config.dg.auxstate[:, q_tot_ind, :]))
    @test max_q_tot < FT(0.0077)

    # no ice
    max_q_ice = maximum(abs.(solver_config.dg.auxstate[:, q_ice_ind, :]))
    @test isequal(max_q_ice, FT(0))

    # q_liq ∈ reference range
    max_q_liq = max(solver_config.dg.auxstate[:, q_liq_ind, :]...)
    min_q_liq = min(solver_config.dg.auxstate[:, q_liq_ind, :]...)
    @test max_q_liq < FT(1e-3)
    @test abs(min_q_liq) < FT(1e-5)

    # q_rai ∈ reference range
    max_q_rai = max(solver_config.dg.auxstate[:, q_rai_ind, :]...)
    min_q_rai = min(solver_config.dg.auxstate[:, q_rai_ind, :]...)
    @test max_q_rai < FT(3e-5)
    @test abs(min_q_rai) < FT(3e-8)

    # terminal velocity ∈ reference range
    max_rain_w = max(solver_config.dg.auxstate[:, rain_w_ind, :]...)
    min_rain_w = min(solver_config.dg.auxstate[:, rain_w_ind, :]...)
    @test max_rain_w < FT(4)
    @test isequal(min_rain_w, FT(0))
end

main()

max_stars_repo_name: trontrytel/CLIMA

max_issues_repo_name: trontrytel/CLIMA

max_forks_repo_name: trontrytel/CLIMA


text:
export
ruleVBGammaOut,
ruleSPGammaOutNPP,
ruleVBGammaA,
ruleVBGammaB

ruleSPGammaOutNPP(  msg_out::Nothing,
                    msg_a::Message{PointMass, Univariate},
                    msg_b::Message{PointMass, Univariate}) =
    Message(Univariate, Gamma, a=deepcopy(msg_a.dist.params[:m]), b=deepcopy(msg_b.dist.params[:m]))

ruleVBGammaOut( dist_out::Any,
                dist_a::ProbabilityDistribution{Univariate},
                dist_b::ProbabilityDistribution{Univariate}) =
    Message(Univariate, Gamma, a=unsafeMean(dist_a), b=unsafeMean(dist_b))

ruleVBGammaA( dist_out::ProbabilityDistribution{Univariate},
                dist_a::Any,
                dist_b::ProbabilityDistribution{Univariate}) =
    Message(Univariate, Function, log_pdf = (a)->a*unsafeLogMean(dist_b) + (a-1)*unsafeLogMean(dist_out) - logabsgamma(a)[1])

ruleVBGammaB( dist_out::ProbabilityDistribution{Univariate},
              dist_a::ProbabilityDistribution{Univariate},
              dist_b::Any) =
        Message(Univariate, Gamma, a=unsafeMean(dist_a)+1, b=unsafeMean(dist_out))

ext: jl


text:
print("Input an integer: ")
n=parse(Int,readline())
println(rand(1:n,n,n))
@show b=3
c=4

max_stars_repo_path: src/Module A - Short Julia Course/deploy2.jl

max_stars_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course

max_issues_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course

max_forks_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course


text:
struct Discretization{WT, GT <: AbstractRange}
    size::Tuple{Int,Int}
    gridx::GT
    gridy::GT
    window::WT
end
function Discretization(grid::AbstractRange, window)
    n = length(grid)
    Discretization((n, n), grid, grid, window)
end

function issquare(d::Discretization)
    step(d.gridx) == step(d.gridy)
end

function discretize(x, x0, sx, x1)
    x0 <= x <= x1 || error("Not all points are in the window.")
    xi = 1+floor(Int, (x-x0)/sx)
    if x==x1 xi -= 1 end
    xi
end
# function cellcounts(::Type{SparseArrays.SparseVector}, d::Discretization, xy::PointPattern)
#
# end
function cellcounts(::Type{Vector{Tuple{Int,Int}}}, d::Discretization, xy::PointPattern)
    counts = cellcounts(Matrix, d, xy)

    [(k,v) for (k,v) in enumerate(counts) if v != 0]
end

function cellcounts(::Type{Matrix}, d::Discretization, xy::PointPattern)
    x0, x1 = d.window.x
    y0, y1 = d.window.y
    sx = (x1-x0)/d.size[1]
    sy = (y1-y0)/d.size[2]
    counts = zeros(Int, d.size)
    for (x, y) in xy
        counts[discretize(x, x0, sx, x1), discretize(y, y0, sy, y1)] += 1
    end
    counts
end


text:
# # Operations problems with [PowerSimulations.jl](https://github.com/NREL-SIIP/PowerSimulations.jl)

# **Originally Contributed by**: Clayton Barrows

# ## Introduction

# PowerSimulations.jl supports the construction and solution of optimal power system
# scheduling problems (Operations Problems). Operations problems form the fundamental
# building blocks for [sequential simulations](../../notebook/3_PowerSimulations_examples/sequential_simulations.ipynb).
# This example shows how to specify and customize a the mathematics that will be applied to the data with
# an `OperationsProblemTemplate`, build and execute an `OperationsProblem`, and access the results.

# ## Dependencies
using SIIPExamples

# ### Modeling Packages
using PowerSystems
using PowerSimulations
const PSI = PowerSimulations
using D3TypeTrees

# ### Data management packages
using Dates
using DataFrames

# ### Optimization packages
using Cbc #solver

# ### Data
# This data depends upon the [RTS-GMLC](https://github.com/gridmod/rts-gmlc) dataset. Let's
# download and extract the data.

rts_dir = SIIPExamples.download("https://github.com/GridMod/RTS-GMLC")
rts_src_dir = joinpath(rts_dir, "RTS_Data", "SourceData")
rts_siip_dir = joinpath(rts_dir, "RTS_Data", "FormattedData", "SIIP");

# ### Create a `System` from RTS-GMLC data just like we did in the [parsing tabular data example.](../../notebook/2_PowerSystems_examples/parse_tabulardata.jl)
rawsys = PowerSystems.PowerSystemTableData(
    rts_src_dir,
    100.0,
    joinpath(rts_siip_dir, "user_descriptors.yaml"),
    timeseries_metadata_file = joinpath(rts_siip_dir, "timeseries_pointers.json"),
    generator_mapping_file = joinpath(rts_siip_dir, "generator_mapping.yaml"),
);
sys = System(rawsys; time_series_resolution = Dates.Hour(1));

# ## Define a problem specification with an `OpModelTemplate`
# The `DeviceModel` constructor is to create an assignment between PowerSystems device types
# and the subtypes of `AbstractDeviceFormulation`. PowerSimulations has a variety of different
# `AbstractDeviceFormulation` subtypes that can be applied to different PowerSystems device types,
# each dispatching to different methods for populating optimization problem objectives, variables,
# and constraints.

#nb TypeTree(PSI.AbstractDeviceFormulation, scopesep="\n")

# ### Branch Formulations
# Here is an example of relatively standard branch formulations. Other formulations allow
# for selective enforcement of transmission limits and greater control on transformer settings.

branches = Dict{Symbol, DeviceModel}(
    :L => DeviceModel(Line, StaticLine),
    :T => DeviceModel(Transformer2W, StaticTransformer),
    :TT => DeviceModel(TapTransformer, StaticTransformer),
)

# ### Injection Device Formulations
# Here we define dictionary entries for all devices that inject or withdraw power on the
# network. For each device type, we can define a distinct `AbstractDeviceFormulation`. In
# this case, we're defining a basic unit commitment model for thermal generators,
# curtailable renewable generators, and fixed dispatch (net-load reduction) formulations
# for `HydroFix` and `RenewableFix` devices. Additionally, we've enabled a simple load
# shedding demand response formulation for `InterruptableLoad` devices.

devices = Dict(
    :Generators => DeviceModel(ThermalStandard, ThermalStandardUnitCommitment),
    :Ren => DeviceModel(RenewableDispatch, RenewableFullDispatch),
    :Loads => DeviceModel(PowerLoad, StaticPowerLoad),
    :HydroROR => DeviceModel(HydroDispatch, FixedOutput),
    :Hydro => DeviceModel(HydroEnergyReservoir, HydroDispatchRunOfRiver),
    :RenFx => DeviceModel(RenewableFix, FixedOutput),
    :ILoads => DeviceModel(InterruptibleLoad, InterruptiblePowerLoad),
)

# ### Service Formulations
# We have two `VariableReserve` types, parameterized by their direction. So, similar to
# creating `DeviceModel`s, we can create `ServiceModel`s. The primary difference being
# that `DeviceModel` objects define how constraints get created, while `ServiceModel` objects
# define how constraints get modified.

services = Dict(
    :ReserveUp => ServiceModel(VariableReserve{ReserveUp}, RangeReserve),
    :ReserveDown => ServiceModel(VariableReserve{ReserveDown}, RangeReserve),
)

# ### Wrap it up into an `OperationsProblemTemplate`
template_uc = OperationsProblemTemplate(CopperPlatePowerModel, devices, branches, services);

# ## `OperationsProblem`
# Now that we have a `System` and an `OperationsProblemTemplate`, we can put the two together
# to create an `OperationsProblem` that we solve.

# ### Optimizer
# It's most convenient to define an optimizer instance upfront and pass it into the
# `OperationsProblem` constructor. For this example, we can use the free Cbc solver with a
# relatively relaxed MIP gap (`ratioGap`) setting to improve speed.

solver = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 1, "ratioGap" => 0.5)

# ### Build an `OperationsProblem`
# The construction of an `OperationsProblem` essentially applies an `OperationsProblemTemplate`
# to `System` data to create a JuMP model.

horizon = 24;
interval = Dates.Hour(24);
transform_single_time_series!(sys, horizon, interval)

op_problem = OperationsProblem(
    GenericOpProblem,
    template_uc,
    sys;
    optimizer = solver,
    horizon = horizon,
)

#nb # The principal component of the `OperationsProblem` is the JuMP model. For small problems,
#nb # you can inspect it by simply printing it to the screen:
#nb # ```julia
#nb # op_problem.psi_container.JuMPmodel
#nb # ```
#nb #
#nb # For anything of reasonable size, that will be unmanageable. But you can print to a file:
#nb # ```julia
#nb # f = open("testmodel.txt","w"); print(f,op_problem.psi_container.JuMPmodel); close(f)
#nb # ```
#nb #
#nb # In addition to the JuMP model, an `OperationsProblem` keeps track of a bunch of metadata
#nb # about the problem and some references to pretty names for constraints and variables.
#nb # All of these details are contained within the `psi_container` field.
#
#nb print_struct(typeof(op_problem.psi_container))
#
#nb # ### Solve an `OperationsProblem`
#
#nb res = solve!(op_problem);
#
#nb # ## Results Inspection
#nb # PowerSimulations collects the `OperationsProblem` results into a struct:
#
#nb print_struct(PSI.SimulationResults)
#
#nb # ### Optimizer Log
#nb # The optimizer summary is included
#
#nb get_optimizer_log(res)
#
#nb # ### Total Cost (objective function value)
#
#nb get_total_cost(res)
#
#nb # ### Variable Values
#nb # The solution value data frames for variable in the `op_problem.psi_container.variables`
#nb # dictionary is stored:
#
#nb variable_values = get_variables(res)
#
#nb # Note that the time stamps are missing from the dataframes in `variable_values`...
#nb #
#nb # The time stamps for each value in the time series used in the `OperationsProblem` is
#nb # included separately from the variable value results.
#
# ## Plotting
# Take a look at the examples in [the plotting folder.](../../notebook/3_PowerSimulations_examples/Plotting)


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

"""
    Crash

Collection of crash methods used to generate initial decisions in structured algorithms.

...
# Available crash methods
- [`None`](@ref)
- [`EVP`](@ref)
- [`Scenario`](@ref)
- [`Custom`](@ref)
...

## Examples

The following solves a stochastic program `sp` created in `StochasticPrograms.jl` using an L-shaped algorithm with trust-region and Clp as an `lpsolver` and by generating an initial decision with the `EVP` crash.

```jldoctest
julia> optimize!(sp, solver = LShapedSolver(GLPKSolverLP(), crash=Crash.EVP(), regularize = TrustRegion()))
L-Shaped Gap  Time: 0:00:00 (8 iterations)
  Objective:       -855.8333333333339
  Gap:             0.0
  Number of cuts:  4
  Iterations:      8
:Optimal
```
"""
module Crash

using StochasticPrograms
using StochasticPrograms: AbstractCrash, StochasticProgramOptimizer
using MathOptInterface

const MOI = MathOptInterface

"""
    None

Randomize the initial decision (default).

"""
struct None <: AbstractCrash end

function (::None)(stochasticprogram::StochasticProgram)
    return rand(num_decisions(stochasticprogram))
end

function (::None)(stochasticmodel::StochasticModel, sampler::AbstractSampler)
    # Get instance optimizer
    optimizer = MOI.get(stochasticmodel, InstanceOptimizer())
    sp = instantiate(stochasticmodel, sampler, 1; optimizer = optimizer)
    return rand(num_decisions(sp))
end

"""
    EVP

Solve the expected value problem corresponding to the stochastic program and use the expected value solution as initial decision.

"""
struct EVP <: AbstractCrash end

function (::EVP)(stochasticprogram::StochasticProgram)
    return expected_value_decision(stochasticprogram)
end

function (crash::EVP)(stochasticmodel::StochasticModel, sampler::AbstractSampler)
    # Get instance optimizer
    optimizer = MOI.get(stochasticmodel, InstanceOptimizer())
    sp = instantiate(stochasticmodel, sampler, 10; optimizer = optimizer)
    return crash(sp)
end

"""
    FeasiblePoint

Generate a feasible first-stage decision as initial decision.

"""
struct FeasiblePoint <: AbstractCrash end

function (::FeasiblePoint)(stochasticprogram::StochasticProgram)
    # Generate first-stage
    m = stage_one_model(stochasticprogram, optimizer = master_optimizer(stochasticprogram))
    # Solve feasibility problem
    @objective(m, MOI.FEASIBILITY_SENSE, 0)
    optimize!(m)
    # Return feasible point
    return JuMP.value.(all_decision_variables(m, 1))
end

function (crash::FeasiblePoint)(stochasticmodel::StochasticModel, sampler::AbstractSampler)
    # Get instance optimizer
    optimizer = MOI.get(stochasticmodel, InstanceOptimizer())
    sp = instantiate(stochasticmodel, sampler, 10; optimizer = optimizer)
    return crash(sp)
end

"""
    Scenario

Solve the wait-and-see problem corresponding a supplied scenario and use the optimal solution as initial decision.

"""
struct Scenario{S <: AbstractScenario} <: AbstractCrash
    scenario::S

    function Scenario(scenario::S) where S <: AbstractScenario
        return new{S}(scenario)
    end
end

function (crash::Scenario)(stochasticprogram::StochasticProgram)
    return wait_and_see_decision(stochasticprogram, crash.scenario)
end

function (crash::Scenario)(stochasticmodel::StochasticModel, sampler::AbstractSampler)
    # Get instance optimizer
    optimizer = MOI.get(stochasticmodel, InstanceOptimizer())
    sp = instantiate(stochasticmodel, sampler, 10; optimizer = optimizer)
    return crash(sp)
end

"""
    PreSolve

Run a provided (ideally suboptimal) optimization procedure and use the (sub)optimal solution as initial decision.

"""
struct PreSolve <: AbstractCrash
    optimizer::StochasticProgramOptimizer

    function PreSolve(optimizer_constructor)
        return new(StochasticProgramOptimizer(optimizer_constructor))
    end
end

function (crash::PreSolve)(stochasticprogram::StochasticProgram)
    optimizer = crash.optimizer.optimizer
    structure_ = structure(stochasticprogram)
    # Generate random starting point
    x₀ = None()(stochasticprogram)
    # Sanity checks
    if !supports_structure(optimizer, structure_)
        @warn "The provided presolve optimizer does not support the underlying structure of the stochastic program. Defaulting to random starting point."
        return x₀
    end
    check_loadable(optimizer, structure_)
    # Better printing
    MOI.set(optimizer, MOI.RawParameter("keep"), false)
    # Load structure
    load_structure!(optimizer, structure_, x₀)
    # Run presolve
    MOI.optimize!(optimizer)
    # Get solution
    x₀ = map(all_decision_variables(stochasticprogram, 1)) do dvar
        return MOI.get(optimizer, MOI.VariablePrimal(), index(dvar))
    end
    # Restore problem
    restore_structure!(optimizer)
    return x₀
end

function (crash::PreSolve)(stochasticmodel::StochasticModel, sampler::AbstractSampler)
    # Get instance optimizer
    sp = instantiate(stochasticmodel, sampler, 10; optimizer = crash.optimizer.optimizer_constructor)
    return crash(sp)
end

"""
    Custom(x₀)

Use the user-supplied `x₀` as initial decision.

"""
struct Custom{T <: AbstractFloat} <: AbstractCrash
    x₀::Vector{T}

    function Custom(x₀::Vector{T}) where T <: AbstractFloat
        return new{T}(x₀)
    end
end

function (crash::Custom)(stochasticprogram::StochasticProgram)
    return crash.x₀[1:num_decisions(stochasticprogram)]
end

function (crash::Custom)(stochasticmodel::StochasticModel, sampler::AbstractSampler)
    # Get instance optimizer
    optimizer = MOI.get(stochasticmodel, InstanceOptimizer())
    sp = instantiate(stochasticmodel, sampler, 1; optimizer = optimizer)
    return crash.x₀[1:num_decisions(sp)]
end

end


text:
# setup paths -
function _include_my_codes(path_to_src::String, fileNameArray::Array{String,1})

    try

        for file_name in fileNameArray

            path_to_load = joinpath(path_to_src, file_name)
            include(path_to_load)
        end

    catch error
        rethrow(error)
    end

end

# get files from dir -
searchdir(path,key) = filter(x -> contains(x, key), readdir(path))

# setup project paths -
const _PATH_TO_ROOT = dirname(pathof(@__MODULE__))
const _PATH_TO_BASE = joinpath(_PATH_TO_ROOT, "base")
const _PATH_TO_DATABASE = joinpath(_PATH_TO_ROOT, "database")

# load packages -
using DataFrames
using CSV
using SQLite
using UUIDs
using JSON

# load my codes -
my_code_array = searchdir(_PATH_TO_BASE, ".jl")
_include_my_codes(_PATH_TO_BASE, my_code_array)

# create a database connection -
const path_to_database_file = joinpath(_PATH_TO_DATABASE, "Parameters.db")


text:
module TrainView

export TracksAhead,TrackProperties, OpticalProperties, TrainCar
export VideoCamera, focal_length, focalθφ, track_image_u, camera_image, left_track_image_u, right_track_image_u
export camera_distoration, rail_uvs_to_distortion, Δu, left_track_Δu, right_track_Δu

# methods that project small changes in track curvature to image of track
export v_to_dudβ, v_to_dudα

import StaticArrays: SVector, SMatrix

using LinearAlgebra
using Statistics
using RecipesBase

include("camera_types.jl")
include("camera_distortion.jl")
include("../plot/tracksahead.jl")

end


text:
import JuMP: GenericAffExpr

init_neurons(model::Model, layers::Vector{Layer})     = init_variables(model, layers, include_input = true)
init_deltas(model::Model, layers::Vector{Layer})      = init_variables(model, layers, binary = true)
init_multipliers(model::Model, layers::Vector{Layer}) = init_variables(model, layers)
# Allow ::Network input also (NOTE for legacy purposes mostly...)
init_neurons(m,     network::Network) = init_neurons(m, network.layers)
init_deltas(m,      network::Network) = init_deltas(m,  network.layers)
init_multipliers(m, network::Network) = init_multipliers(m, network.layers)

function init_variables(model::Model, layers::Vector{Layer}; binary = false, include_input = false)
    # TODO: only neurons get offset array
    vars = Vector{Vector{VariableRef}}(undef, length(layers))
    all_layers_n = n_nodes.(layers)

    if include_input
        # input to the first layer also gets variables
        # essentially an input constraint
        input_layer_n = size(first(layers).weights, 2)
        prepend!(all_layers_n, input_layer_n)
        push!(vars, Vector{VariableRef}())        # expand vars by one to account
    end

    for (i, n) in enumerate(all_layers_n)
        vars[i] = @variable(model, [1:n], binary = binary, base_name = "z$i")
    end
    return vars
end

function symbolic_max(m::Model, a, b)
    aux = @variable(m)
    @constraint(m, aux >= a)
    @constraint(m, aux >= b)
    return aux
end
symbolic_max(a::VariableRef, b::VariableRef)                         = symbolic_max(a.model, a, b)
symbolic_max(a::GenericAffExpr, b::GenericAffExpr)                   = symbolic_max(first(first(a.terms)).model, a, b)
symbolic_max(a::Array{<:GenericAffExpr}, b::Array{<:GenericAffExpr}) = symbolic_max.(a, b)

function symbolic_abs(m::Model, v)
    aux = @variable(m) #get an anonymous variable
    @constraint(m, aux >= 0)
    @constraint(m, aux >= v)
    @constraint(m, aux >= -v)
    return aux
end
symbolic_abs(v::VariableRef)             = symbolic_abs(v.m, v)
symbolic_abs(v::GenericAffExpr)          = symbolic_abs(first(first(v.terms)).model, v)
symbolic_abs(v::Array{<:GenericAffExpr}) = symbolic_abs.(v)

function symbolic_infty_norm(m::Model, v::Array{<:GenericAffExpr})
    aux = @variable(m)
    @constraint(m, aux >= 0)
    @constraint(m, aux .>= v)
    @constraint(m, aux .>= -v)
    return aux
end
# # in general, default to symbolic_abs behavior:
# symbolic_infty_norm(v) = symbolic_abs(v)
# only Array{<:GenericAffExpr} is needed
symbolic_infty_norm(v::Array{<:GenericAffExpr}) = symbolic_infty_norm(first(first(first(v).terms)).model, v)


text:
function updateDatastream!(
    mb    :: ModelBlock,
    clock :: ModelClock,
)
    cdatam     = mb.co.cdatam
    datastream = mb.tmpfi.datastream
    cfg        = mb.ev.config

    if cdatam != nothing

        interpData!(cdatam, clock.time, datastream)

        if cfg["MLD_scheme"] == "datastream"
            mb.fi.HMXL .= datastream["HMXL"]
        end

        if cfg["Qflx"] == "on"
            mb.fi._QFLXX_[:, 1] .= reshape(datastream["QFLXT"], :)
            mb.fi._QFLXX_[:, 2] .= reshape(datastream["QFLXS"], :)
        end

    end

end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
using Test
using Solvent
using StaticArrays
using LinearAlgebra
using SparseArrays
using Random
using Printf

# this test setup is partly based on IterativeSolvers.jl, see e.g
# https://github.com/JuliaMath/IterativeSolvers.jl/blob/master/test/cg.jl
@testset "Solvent small full system: GMRES" begin
    n = 10

    for T in [Float32, Float64]
        Random.seed!(44)

        A = @MMatrix rand(T, n, n)
        A = A' * A + I
        b = @MVector rand(T, n)

        mulbyA!(y, x) = (y .= A * x)

        x = @MVector rand(T, n)

        tol = sqrt(eps(T))
        solver_type = GeneralizedMinimalResidualMethod(M = n, K = 1)
        linearsolver = LinearSolver(
            mulbyA!,
            solver_type,
            x;
            rtol = tol,
            atol = tol,
        )
        x0 = copy(x)
        linearsolve!(linearsolver, x, b)
        @test norm(A * x - b) / norm(A * x0 - b) <= tol

        # test for convergence in 0 iterations by
        # initializing with the exact solution
        x = A \ b
        iters = linearsolve!(linearsolver, x, b)
        @test iters == 0
        @test norm(A * x - b) <= 100 * eps(T)

        newtol = 1000tol
        settolerance!(linearsolver, newtol)
        settolerance!(linearsolver, newtol; relative = true)

        x = @MVector rand(T, n)
        x0 = copy(x)
        linearsolve!(linearsolver, x, b)

        @test norm(A * x - b) / norm(A * x0 - b) <= newtol
        @test norm(A * x - b) / norm(A * x0 - b) >= tol

    end
end

@testset "Solvent large sparse system: GMRES" begin
    n = 10000
    M = 20
    K = 10

    for T in [Float32, Float64]
        Random.seed!(44)

        α = 1f-2
        A = I + α * sprandn(T, n, n, 0.05)
        b = rand(T, n)

        mulbyA!(y, x) = (y .= A * x)

        x = rand(T, n)

        tol = sqrt(eps(T))
        solver_type = GeneralizedMinimalResidualMethod(M = M, K = K)
        preconditioner = Identity(pc_side=PCright())
        linearsolver = LinearSolver(
            mulbyA!,
            solver_type,
            x;
            pc_alg = preconditioner,
            rtol = tol,
            atol = tol,
        )

        x0 = copy(x)
        linearsolve!(linearsolver, x, b)
        @test norm(A * x - b) / norm(A * x0 - b) <= tol

        # test for convergence in 0 iterations by
        # initializing with the exact solution
        x = A \ b
        iters = linearsolve!(linearsolver, x, b)
        @test iters == 0
        @test norm(A * x - b) <= tol

        newtol = 1000tol
        settolerance!(linearsolver, newtol)
        settolerance!(linearsolver, newtol; relative=true)

        x = rand(T, n)
        x0 = copy(x)
        linearsolve!(linearsolver, x, b)

        @test norm(A * x - b) / norm(A * x0 - b) <= newtol
    end
end

max_stars_repo_name: CliMA/SystemSolvers

max_issues_repo_name: CliMA/SystemSolvers


text:
# simple example of loading Flux data and examining.
using Flux, BenchmarkTools, Statistics

labels      = Flux.Data.Iris.labels()
features    = Flux.Data.Iris.features()

display(features[1,1:2])

n_features = Flux.normalise(features,dims = 2)

display(n_features[1,1:2])

display(size(n_features))
klasses = sort(unique(labels))
display(klasses)
hlabels = Flux.onehotbatch(labels,klasses)

display(hlabels)
display(size(labels))
display(size(hlabels))
# now bring back the old labels from one-hot encoding
old_labels = Flux.onecold(hlabels,klasses)
@assert labels == old_labels


text:
using Pkg.Artifacts

using ArtifactWrappers

# Get dycoms dataset folder:
dycoms_dataset = ArtifactWrapper(
    @__DIR__,
    isempty(get(ENV, "CI", "")),
    "dycoms",
    ArtifactFile[ArtifactFile(
        url = "https://caltech.box.com/shared/static/bxau6i46y6ikxn2sy9krgz0sw5vuptfo.nc",
        filename = "test_data_PhaseEquil.nc",
    ),],
)
dycoms_dataset_path = get_data_folder(dycoms_dataset)

@testset "Data tests" begin
    FT = Float64

    data = joinpath(dycoms_dataset_path, "test_data_PhaseEquil.nc")
    ds_PhaseEquil = Dataset(data, "r")
    e_int = Array{FT}(ds_PhaseEquil["e_int"][:])
    ρ = Array{FT}(ds_PhaseEquil["ρ"][:])
    q_tot = Array{FT}(ds_PhaseEquil["q_tot"][:])

    ts = PhaseEquil_ρeq.(Ref(param_set), ρ, e_int, q_tot, 4)
    # ts = PhaseEquil_ρeq.(Ref(param_set), ρ, e_int, q_tot, 3) # Fails
end

@testset "pθq data-driven tests" begin
    #! format: off
    pθq_broken = [
        # (; p = , θ_liq_ice = , q_tot = ),
        # (; p = , θ_liq_ice = , q_tot = ),
        # (; p = , θ_liq_ice = , q_tot = ),
        # (; p = , θ_liq_ice = , q_tot = ),
        # (; p = , θ_liq_ice = , q_tot = ),
    ]
    #! format: on
    # config = (50, 1e-3, RootSolvers.RegulaFalsiMethod)
    # config = (50, 1e-3, RootSolvers.NewtonMethodAD)
    config = ()
    if !isempty(pθq_broken)
        p_arr = getproperty.(pθq_broken, :p)
        θ_liq_ice_arr = getproperty.(pθq_broken, :θ_liq_ice)
        q_tot_arr = getproperty.(pθq_broken, :q_tot)
        for (p, θ_liq_ice, q_tot) in zip(p_arr, θ_liq_ice_arr, q_tot_arr)
            @test_broken begin
                ts = PhaseEquil_pθq(param_set, p, θ_liq_ice, q_tot, config...)
                air_pressure(ts) == p
            end
        end
    end

    #! format: off
    pθq = [
        (; p = 82307.30319719888, θ_liq_ice = 296.7074342326652, q_tot = 0.01894019026929829,),
        (; p = 88357.42589002676, θ_liq_ice = 296.7074342326652, q_tot = 0.01894019026929829,),
        (; p = 81090.35731696963, θ_liq_ice = 296.7074342326652, q_tot = 0.01894019026929829,),
        (; p = 75461.41343839701, θ_liq_ice = 296.7074342326652, q_tot = 0.01894019026929829,),
        (; p = 71158.23329080557, θ_liq_ice = 296.7074342326652, q_tot = 0.01894019026929829,),
        (; p = 68032.73180723468, θ_liq_ice = 296.7074342326652, q_tot = 0.01894019026929829,),
        (; p = 59906.16044902247, θ_liq_ice = 296.7074342326652, q_tot = 0.01894019026929829,),
        (; p = 51740.77281055945, θ_liq_ice = 296.7074342326652, q_tot = 0.01894019026929829,),
        (; p = 50056.7894012541, θ_liq_ice = 296.7074342326652, q_tot = 0.01894019026929829,),
        (; p = 48806.70567178993, θ_liq_ice = 296.7074342326652, q_tot = 0.01894019026929829,),
        (; p = 27564.73889538213, θ_liq_ice = 296.7074342326652, q_tot = 0.01894019026929829,),
    ]
    #! format: on
    if !isempty(pθq)
        p_arr = getproperty.(pθq, :p)
        θ_liq_ice_arr = getproperty.(pθq, :θ_liq_ice)
        q_tot_arr = getproperty.(pθq, :q_tot)
        for (p, θ_liq_ice, q_tot) in zip(p_arr, θ_liq_ice_arr, q_tot_arr)
            @test begin
                ts = PhaseEquil_pθq(param_set, p, θ_liq_ice, q_tot, config...)
                air_pressure(ts) == p
            end
        end
    end
end


text:
module Clover
    using LinearAlgebra
    import ..Actions:FermiActionParam_WilsonClover
    import ..LieAlgebrafields:clear!
    import ..Gaugefields:GaugeFields_1d,GaugeFields,set_wing!,
                substitute!,gauge_shift!,lambdamul,
                SU3GaugeFields_1d,SU2GaugeFields_1d,
                SUNGaugeFields_1d,
                Loops,evaluate_loops!,evaluate_loops,SU
    import ..Wilsonloops:Wilson_loop,Wilson_loop_set

    function Make_CloverFμν(fparam::FermiActionParam_WilsonClover,U,temps::Array{T,1})  where T <: GaugeFields_1d
        NV = temps[1].NV
        NC = temps[1].NC
        CloverFμν = zeros(ComplexF64,NC,NC,NV,6)
        Make_CloverFμν!(CloverFμν,fparam,U,temps)
        return CloverFμν
    end

    function Make_CloverFμν!(CloverFμν,fparam::FermiActionParam_WilsonClover,U,temps::Array{T,1}) where T <: GaugeFields_1d
        work1 = temps[1]
        work2 = temps[2]
        work3 = temps[3]
        work4 = temps[4]

        #CloverFμν .= 0
        set_wing!(U)

        NV = work1.NV
        NC = work1.NC

        # ... Calculation of 4 leaves under the counter clock order.
        μν = 0
        for μ=1:3
            for ν=μ+1:4
                μν += 1
                if μν > 6
                    error("μν > 6 ?")
                end

                loopset = Loops(U,fparam._cloverloops[μ,ν],[work2,work3,work4])
                evaluate_loops!(work1,loopset,U)
                setFμν!(CloverFμν,μν,work1)

            end
        end

        for μν = 1:6
            cimaglink!(CloverFμν,μν,NV,NC)
        end

        coe  = im*0.125*fparam.hop*fparam.Clover_coefficient
        CloverFμν .*= coe

        #println("clover ",sum(abs.(fparam.CloverFμν)))

        return

    end

    function Make_CloverFμν!(fparam::FermiActionParam_WilsonClover,U,temps::Array{T,1}) where T <: GaugeFields_1d
        Make_CloverFμν!(fparam.CloverFμν,fparam,U,temps)
        return
    end

    """
     y = x - x_aj
    """
    function cimaglink!(Fμν,μν,NV,NC)
        if NC == 3
            for i=1:NV
                z11 = Fμν[1,1,i,μν] - conj(Fμν[1,1,i,μν] )
                z12 = Fμν[1,2,i,μν] - conj(Fμν[2,1,i,μν] )
                z13 = Fμν[1,3,i,μν] - conj(Fμν[3,1,i,μν] )

                z22 = Fμν[2,2,i,μν] - conj(Fμν[2,2,i,μν] )
                z23 = Fμν[2,3,i,μν] - conj(Fμν[3,2,i,μν] )

                z33 = Fμν[3,3,i,μν] - conj(Fμν[3,3,i,μν] )

                Fμν[1,1,i,μν] = z11
                Fμν[1,2,i,μν] = z12
                Fμν[1,3,i,μν] = z13

                Fμν[2,1,i,μν] = -conj(z12)
                Fμν[2,2,i,μν] = z22
                Fμν[2,3,i,μν] = z23

                Fμν[3,1,i,μν] = -conj(z13)
                Fμν[3,2,i,μν] = -conj(z23)
                Fμν[3,3,i,μν] = z33
            end
        elseif NC==2
            for i=1:NV
                z11 = Fμν[1,1,i,μν] - conj(Fμν[1,1,i,μν] )
                z12 = Fμν[1,2,i,μν] - conj(Fμν[2,1,i,μν] )

                z22 = Fμν[2,2,i,μν] - conj(Fμν[2,2,i,μν] )

                Fμν[1,1,i,μν] = z11
                Fμν[1,2,i,μν] = z12

                Fμν[2,1,i,μν] = -conj(z12)
                Fμν[2,2,i,μν] = z22

            end
        else
            for i=1:NV
                for k2=1:NC
                    for k1=k2:NC
                        z = Fμν[k1,k2,i,μν] - conj(Fμν[k2,k1,i,μν] )
                        Fμν[k1,k2,i,μν] = z
                        if k1 != k2
                            Fμν[k2,k1,i,μν] = -conj(z)
                        end
                    end
                end
            end
            #error("error!")
        end
        return
    end

    function setFμν!(Fμν,μν,v::GaugeFields_1d{SU{NC}}) where NC
        NV=v.NV
        for i=1:NV
            for k2=1:NC
                for k1=1:NC
                    Fμν[k1,k2,i,μν] = v[ k1,k2,i]
                end
            end
        end
        return
    end

    function addFμν!(Fμν,μν,v::SU3GaugeFields_1d)
        NV=v.NV

        for i=1:NV
            Fμν[1,1,i,μν] += v[ 1,1,i]
            Fμν[1,2,i,μν] += v[ 1,2,i]
            Fμν[1,3,i,μν] += v[ 1,3,i]
            Fμν[2,1,i,μν] += v[ 2,1,i]
            Fμν[2,2,i,μν] += v[ 2,2,i]
            Fμν[2,3,i,μν] += v[ 2,3,i]
            Fμν[3,1,i,μν] += v[ 3,1,i]
            Fμν[3,2,i,μν] += v[ 3,2,i]
            Fμν[3,3,i,μν] += v[ 3,3,i]
        end
        return
    end

    function addFμν!(Fμν,μν,v::SU2GaugeFields_1d)
        NV=v.NV

        for i=1:NV
            Fμν[1,1,i,μν] += v[ 1,1,i]
            Fμν[1,2,i,μν] += v[ 1,2,i]
            Fμν[2,1,i,μν] += v[ 2,1,i]
            Fμν[2,2,i,μν] += v[ 2,2,i]

        end
        return
    end

    function addFμν!(Fμν,μν,v::SUNGaugeFields_1d)
        NV=v.NV
        NC = v.NC

        for i=1:NV
            for k2=1:NC
                for k1=1:NC
                    Fμν[k1,k2,i,μν] += v[ k1,k2,i]
                end
            end
        end
        return
    end

    jf(I,L) = (I+1) % L
    jb(I,L) = (I-1+L) %L
    ic(ix,iy,iz,it,N1,N2,N3) = 1+ix+N1*(iy+N2*(iz+N3*it))

    function fprep3!(NX,NY,NZ,NT,inn)
        for it=0:NT-1
            for iz=0:NZ-1
                for iy=0:NY-1
                    for ix=0:NX-1
                        i = ic(ix,iy,iz,it,NX,NY,NZ)
                        inn[i,1,1] = ic(jf(ix,NX),iy       ,iz       ,it       ,NX,NY,NZ)
                        inn[i,2,1] = ic(ix       ,jf(iy,NY),iz       ,it       ,NX,NY,NZ)
                        inn[i,3,1] = ic(ix       ,iy       ,jf(iz,NZ),it       ,NX,NY,NZ)
                        inn[i,4,1] = ic(ix       ,iy       ,iz       ,jf(it,NT),NX,NY,NZ)
                        inn[i,1,2] = ic(jb(ix,NX),iy       ,iz       ,it       ,NX,NY,NZ)
                        inn[i,2,2] = ic(ix       ,jb(iy,NY),iz       ,it       ,NX,NY,NZ)
                        inn[i,3,2] = ic(ix       ,iy       ,jb(iz,NZ),it       ,NX,NY,NZ)
                        inn[i,4,2] = ic(ix       ,iy       ,iz       ,jb(it,NT),NX,NY,NZ)
                    end
                end
            end
        end
        return
    end

    """
    Calculate   dS_clover/dA_mu(x)
    """
    function dSclover!(z,μ,X,Y,U,fparam::FermiActionParam_WilsonClover,temps::Array{T,1}) where T <: GaugeFields_1d
        NX = X.NX
        NY = X.NY
        NZ = X.NZ
        NT = X.NT
        NC = X.NC
        NV = NX*NY*NZ*NT
        numbasis = NC^2-1
        #numbasis = ifelse(NC==3,8,3)

        work1 = temps[1]
        work2 = temps[2]
        work3 = temps[3]
        work4 = temps[4]

        gtmp1 = temps[5]
        gtmp2 = temps[6]
        gtmp3 = temps[7]
        gtmp4 = temps[8]

        dF1 = temps[8+1:8+numbasis]
        #dF1 = temps[9:16]
        dF2 = temps[8+numbasis+1:8+2numbasis]
        #dF2 = temps[17:24]

        if fparam.internal_flags[1] ==  false
            fprep3!(NX,NY,NZ,NT,fparam.inn_table)
            fparam.internal_flags[1] = true
        end

        clear!(z)
        veta = fparam._ftmp_vectors[1]
        vxi = fparam._ftmp_vectors[2]
        ftmp1 = fparam._ftmp_vectors[3]
        ftmp2 = fparam._ftmp_vectors[4]
        v1 = fparam._ftmp_vectors[5]
        v2 = fparam._ftmp_vectors[6]
        is1 = fparam._is1
        is2 = fparam._is2
        inn = fparam.inn_table

        for ialpha = 1:4
            for ic=1:NC
                is = 0
                for it=1:NT
                    for iz=1:NZ
                        for iy=1:NY
                            for ix=1:NX
                                is += 1
                                veta[ic,is,ialpha] = Y[ic,ix,iy,iz,it,ialpha]
                                vxi[ic,is,ialpha] = X[ic,ix,iy,iz,it,ialpha]
                            end
                        end
                    end
                end
            end
        end

        for ν=1:4
            if ν==μ
                continue
            end

            #=
             .....................
                Case 1 and 3
             .....................
            =#
            iflag = 1
            cal_dFμν!(dF1,dF2,
                        U,fparam,
                        work1,work2,work3,work4,
                        gtmp1,gtmp2,gtmp3,gtmp4,
                        μ,ν,iflag)

            # ... Case 1
            VxSigxV!(veta,vxi,dF1,z,ftmp1,ftmp2,μ,ν)
            #println("z1 = ",z*z)

            # ... Case 3
            for is=1:NV
                is1[is] = inn[is,μ,1]
            end

            for is=1:NV
                is2[is] = inn[is1[is],ν,1]
            end

            for ic=1:NC
                for ialpha=1:4
                    for is=1:NV
                        v1[ic,is,ialpha] = veta[ic,is2[is],ialpha]
                        v2[ic,is,ialpha] = vxi[ic,is2[is],ialpha]
                    end
                end
            end

            VxSigxV!(v1,v2,dF2,z,ftmp1,ftmp2,μ,ν)
            #println("z2 = ",z*z)

            #=
             .....................
                Case 2 and 4
             .....................
            =#

            iflag = 2

            cal_dFμν!(dF1,dF2,
                    U,fparam,
                    work1,work2,work3,work4,
                    gtmp1,gtmp2,gtmp3,gtmp4,
                    μ,ν,iflag)

            # ... Case 2
            for is=1:NV
                is1[is] = inn[is,μ,1]
            end

            for ic=1:NC
                for ialpha=1:4
                    for is=1:NV
                        v1[ic,is,ialpha] = veta[ic,is1[is],ialpha]
                        v2[ic,is,ialpha] = vxi[ic,is1[is],ialpha]
                    end
                end
            end

            VxSigxV!(v1,v2,dF1,z,ftmp1,ftmp2,μ,ν)
            #println("z3 = ",z*z)

            # ... Case 4
            for is=1:NV
                is1[is] = inn[is,ν,1]
            end

            for ic=1:NC
                for ialpha=1:4
                    for is=1:NV
                        v1[ic,is,ialpha] = veta[ic,is1[is],ialpha]
                        v2[ic,is,ialpha] = vxi[ic,is1[is],ialpha]
                    end
                end
            end

            VxSigxV!(v1,v2,dF2,z,ftmp1,ftmp2,μ,ν)
            #println("z4 = ",z*z)

            #=
             .....................
                Case 4' and 2'
             .....................
            =#

            iflag = 3

            cal_dFμν!(dF1,dF2,
                    U,fparam,
                    work1,work2,work3,work4,
                    gtmp1,gtmp2,gtmp3,gtmp4,
                    μ,ν,iflag)

            # ... Case 4'
            VxSigxV!(veta,vxi,dF1,z,ftmp1,ftmp2,μ,ν)
            #println("z5 = ",z*z)

            # ... Case 2'
            for is=1:NV
                is1[is] = inn[is,μ,1]
            end

            for is=1:NV
                is2[is] = inn[is1[is],ν,2]
            end

            for ic=1:NC
                for ialpha=1:4
                    for is=1:NV
                        v1[ic,is,ialpha] = veta[ic,is2[is],ialpha]
                        v2[ic,is,ialpha] = vxi[ic,is2[is],ialpha]
                    end
                end
            end

            VxSigxV!(v1,v2,dF2,z,ftmp1,ftmp2,μ,ν)
            #println("z6 = ",z*z)

            #=
             .....................
                Case 3' and 1'
             .....................
            =#
            iflag = 4

            cal_dFμν!(dF1,dF2,
                    U,fparam,
                    work1,work2,work3,work4,
                    gtmp1,gtmp2,gtmp3,gtmp4,
                    μ,ν,iflag)

            # ... Case 3'
            for is=1:NV
                is1[is] = inn[is,μ,1]
            end

            for ic=1:NC
                for ialpha=1:4
                    for is=1:NV
                        v1[ic,is,ialpha] = veta[ic,is1[is],ialpha]
                        v2[ic,is,ialpha] = vxi[ic,is1[is],ialpha]
                    end
                end
            end

            VxSigxV!(v1,v2,dF1,z,ftmp1,ftmp2,μ,ν)

            # ... Case 1'
            for is=1:NV
                is1[is] = inn[is,ν,2]
            end

            for ic=1:NC
                for ialpha=1:4
                    for is=1:NV
                        v1[ic,is,ialpha] = veta[ic,is1[is],ialpha]
                        v2[ic,is,ialpha] = vxi[ic,is1[is],ialpha]
                    end
                end
            end

            VxSigxV!(v1,v2,dF2,z,ftmp1,ftmp2,μ,ν)

        end

        for ia=1:numbasis
            is = 0
            for it=1:NT
                for iz=1:NZ
                    for iy=1:NY
                        for ix=1:NX
                            is += 1
                            z[ia,ix,iy,iz,it] = -2*z[ia,ix,iy,iz,it]
                        end
                    end
                end
            end

        end

    end

    function VxSigxV!(v1,v2,u,z,tmp1,tmp2,μ,ν)
        NC,NV,_ = size(v1)
        NX = z.NX
        NY = z.NY
        NZ = z.NZ
        NT = z.NT
        #numbasis = ifelse(NC==3,8,3)
        numbasis = NC^2-1

        if μ == ν
            error("""μ should not be equal to ν (VsSigV)
            μ,ν : $μ $ν
            """)
        end

        if μ < ν
            facμν = 1
            μ0 = μ
            ν0 = ν
        else
            facμν = -1
            μ0 = ν
            ν0 = μ
        end

        # ... σ_μν x v2 (in Dirac space)  ....
        if μ0 == 1
            if ν0 == 2
                for ic=1:NC
                    for is=1:NV
                        tmp1[ic,is,1] = -v2[ic,is,1]
                        tmp1[ic,is,2] = +v2[ic,is,2]
                        tmp1[ic,is,3] = -v2[ic,is,3]
                        tmp1[ic,is,4] = +v2[ic,is,4]
                    end
                end

            elseif ν0 == 3
                for ic=1:NC
                    for is=1:NV
                        tmp1[ic,is,1] = -im*v2[ic,is,2]
                        tmp1[ic,is,2] = +im*v2[ic,is,1]
                        tmp1[ic,is,3] = -im*v2[ic,is,4]
                        tmp1[ic,is,4] = +im*v2[ic,is,3]
                    end
                end
            elseif ν0 == 4
                for ic=1:NC
                    for is=1:NV
                        tmp1[ic,is,1] = -v2[ic,is,2]
                        tmp1[ic,is,2] = -v2[ic,is,1]
                        tmp1[ic,is,3] = +v2[ic,is,4]
                        tmp1[ic,is,4] = +v2[ic,is,3]
                    end
                end
            else
                error("""something is wrong in VsSigV
                μ,ν : $μ $ν
                """)
            end
        elseif μ0 == 2
            if ν0 == 3
                for ic=1:NC
                    for is=1:NV
                        tmp1[ic,is,1] = -v2[ic,is,2]
                        tmp1[ic,is,2] = -v2[ic,is,1]
                        tmp1[ic,is,3] = -v2[ic,is,4]
                        tmp1[ic,is,4] = -v2[ic,is,3]
                    end
                end
            elseif ν0 == 4
                for ic=1:NC
                    for is=1:NV
                        tmp1[ic,is,1] = +im*v2[ic,is,2]
                        tmp1[ic,is,2] = -im*v2[ic,is,1]
                        tmp1[ic,is,3] = -im*v2[ic,is,4]
                        tmp1[ic,is,4] = +im*v2[ic,is,3]
                    end
                end
            else
                error("""something is wrong in VsSigV
                μ,ν : $μ $ν
                """)
            end
        elseif μ0 == 3
            if ν0 == 4
                for ic=1:NC
                    for is=1:NV
                        tmp1[ic,is,1] = -v2[ic,is,1]
                        tmp1[ic,is,2] = +v2[ic,is,2]
                        tmp1[ic,is,3] = +v2[ic,is,3]
                        tmp1[ic,is,4] = -v2[ic,is,4]
                    end
                end

            else
                error("""something is wrong in VsSigV
                μ,ν : $μ $ν
                """)
            end
        else
            error("""something is wrong in VsSigV
                μ,ν : $μ $ν
                """)
        end

        for ic=1:NC
            for is=1:NV
                tmp2[ic,is,1] = facμν*tmp1[ic,is,1]
                tmp2[ic,is,2] = facμν*tmp1[ic,is,2]
                tmp2[ic,is,3] = facμν*tmp1[ic,is,3]
                tmp2[ic,is,4] = facμν*tmp1[ic,is,4]
            end
        end

        if NC == 3

            # ... tmp1 = u x tmp2 (in color space)
            for ia = 1:8
                for ialpha=1:4
                    for is=1:NV
                        s1 = + u[ia][1,1,is]*tmp2[1,is,ialpha] +
                                u[ia][1,2,is]*tmp2[2,is,ialpha] +
                            + u[ia][1,3,is]*tmp2[3,is,ialpha]
                        s2 = + u[ia][2,1,is]*tmp2[1,is,ialpha] +
                            u[ia][2,2,is]*tmp2[2,is,ialpha] +
                            u[ia][2,3,is]*tmp2[3,is,ialpha]
                        s3 = + u[ia][3,1,is]*tmp2[1,is,ialpha] +
                            + u[ia][3,2,is]*tmp2[2,is,ialpha] +
                            + u[ia][3,3,is]*tmp2[3,is,ialpha]

                        tmp1[1,is,ialpha] = s1
                        tmp1[2,is,ialpha] = s2
                        tmp1[3,is,ialpha] = s3
                    end
                end

                # ... v1^{\dagger} * tmp1
                is = 0
                for it=1:NT

                    for iz=1:NZ
                        for iy=1:NY
                            for ix=1:NX
                                is += 1
                                s1 = conj(v1[1,is,1]) * tmp1[1,is,1] +
                                        conj(v1[1,is,2]) * tmp1[1,is,2] +
                                        conj(v1[1,is,3]) * tmp1[1,is,3] +
                                        conj(v1[1,is,4]) * tmp1[1,is,4]

                                s2 = conj(v1[2,is,1]) * tmp1[2,is,1] +
                                        conj(v1[2,is,2]) * tmp1[2,is,2] +
                                        conj(v1[2,is,3]) * tmp1[2,is,3] +
                                        conj(v1[2,is,4]) * tmp1[2,is,4]

                                s3 = conj(v1[3,is,1]) * tmp1[3,is,1] +
                                        conj(v1[3,is,2]) * tmp1[3,is,2] +
                                        conj(v1[3,is,3]) * tmp1[3,is,3] +
                                        conj(v1[3,is,4]) * tmp1[3,is,4]

                                z[ia,ix,iy,iz,it] += real(s1+s2+s3)
                                #println(z[ia,ix,iy,iz,it])
                            end
                        end
                    end
                end

            end
        elseif NC == 2

            # ... tmp1 = u x tmp2 (in color space)
            for ia = 1:3
                for ialpha=1:4
                    for is=1:NV
                        s1 = + u[ia][1,1,is]*tmp2[1,is,ialpha] +
                                u[ia][1,2,is]*tmp2[2,is,ialpha]
                        s2 = + u[ia][2,1,is]*tmp2[1,is,ialpha] +
                            u[ia][2,2,is]*tmp2[2,is,ialpha]

                        tmp1[1,is,ialpha] = s1
                        tmp1[2,is,ialpha] = s2

                    end
                end

                # ... v1^{\dagger} * tmp1
                is = 0
                for it=1:NT
                    for iz=1:NZ
                        for iy=1:NY
                            for ix=1:NX
                                is += 1
                                s1 = conj(v1[1,is,1]) * tmp1[1,is,1] +
                                        conj(v1[1,is,2]) * tmp1[1,is,2] +
                                        conj(v1[1,is,3]) * tmp1[1,is,3] +
                                        conj(v1[1,is,4]) * tmp1[1,is,4]

                                s2 = conj(v1[2,is,1]) * tmp1[2,is,1] +
                                        conj(v1[2,is,2]) * tmp1[2,is,2] +
                                        conj(v1[2,is,3]) * tmp1[2,is,3] +
                                        conj(v1[2,is,4]) * tmp1[2,is,4]

                                z[ia,ix,iy,iz,it] += real(s1+s2)
                            end
                        end
                    end
                end

            end
        else
            # ... tmp1 = u x tmp2 (in color space)
            for ia = 1:NC^2-1
                for ialpha=1:4
                    for is=1:NV
                        for k1=1:NC
                            tmp1[k1,is,ialpha] = 0

                            for k2=1:NC
                                tmp1[k1,is,ialpha] += u[ia][k1,k2,is]*tmp2[k2,is,ialpha]
                            end

                        end
                    end
                end

                # ... v1^{\dagger} * tmp1
                is = 0
                for it=1:NT

                    for iz=1:NZ
                        for iy=1:NY
                            for ix=1:NX
                                is += 1

                                for k1=1:NC
                                    z[ia,ix,iy,iz,it] += real(conj(v1[k1,is,1]) * tmp1[k1,is,1] +
                                        conj(v1[k1,is,2]) * tmp1[k1,is,2] +
                                        conj(v1[k1,is,3]) * tmp1[k1,is,3] +
                                        conj(v1[k1,is,4]) * tmp1[k1,is,4])
                                end

                                #println(z[ia,ix,iy,iz,it])
                            end
                        end
                    end
                end

            end
        end

        #println("z = ",z*z)

    end

    function cal_dFμν!(dFμν1,dFμν2,U,fparam::FermiActionParam_WilsonClover,work1::T,work2,work3,work4,temp1,temp2,temp3,temp4,μ,ν,iflag) where T <: GaugeFields_1d
        coe = im*0.125*fparam.hop*fparam.Clover_coefficient
        NC = U[1].NC
        #numbasis = ifelse(NC == 3,8,3)
        numbasis = NC^2-1

        #=
            ...  A) the upper staple  ....................
                                       w3
               nu                 .-----------+
                                  |           |
               /|             w4  |           | w2
                |                 |           |
                |                 |           |
                +----> mu         .-----------.

        =#
        if iflag == 1 || iflag == 2
            substitute!(work1,U[μ])
            gauge_shift!(work2,μ,U[ν])
            gauge_shift!(work3,ν,U[μ])
            substitute!(work4,U[ν])
        end

        # ...  iflag = 1 (Case 1 and 3)
        if iflag == 1
            #=
                                           w3
                                      .-----------+
              Case 1                  |           |
                   nu                 |           |
                    |              w4 |           | w2
                    |                 |           |
                    +----> mu         o-----------.
                                      x    w1

                                           w3
                                      .<----------o
              Case 3                  |           |
                   nu                 |           |
                    |              w4 |           | w2
                    |                 |           |
                    +----> mu         .-----------.

            =#
            mul!(temp1,work1,work2)
            mul!(temp2,work4,work3)

            for ia=1:numbasis
                lambdamul(temp3,temp1,ia,fparam.SUNgenerator)
                mul!(temp4,temp3,temp2')
                mul!(temp4,im)
                cimaglnk!(temp4)
                substitute!(dFμν1[ia],coe,temp4)

                mul!(temp4,temp2',temp3)
                mul!(temp4,im)
                cimaglnk!(temp4)
                substitute!(dFμν2[ia],coe,temp4)

            end

        end
        # ...  iflag = 2 (Case 2 and 4)
        if iflag == 2
            #=
                                            w3
                                      .-----------+
              Case 2                  |           |
                   nu                 |           |
                    |             w4  |           | w2
                    |                 |           |
                    +----> mu         .---------->o
                                      x    w1

                                            w3
                                      o-----------+
              Case 4                  |           |
                   nu                 |           |
                    |              w4 |           | w2
                    |                 |           |
                    +----> mu         .-----------.
                                      x    w1
            =#
            mul!(temp1,work2,work3')

            for ia=1:numbasis
                lambdamul(temp2,work1,ia,fparam.SUNgenerator) # temp2=(lambda_a/2)*work1
                mul!(temp3,work4',temp2)
                mul!(temp4,temp1,temp3)
                mul!(temp4,im)
                cimaglnk!(temp4)
                substitute!(dFμν1[ia],coe,temp4)

                mul!(temp4,temp3,temp1)
                mul!(temp4,im)
                cimaglnk!(temp4)
                substitute!(dFμν2[ia],coe,temp4)
            end

        end

        #=
            ...  B) the lower staple  ....................

               nu
               /|
                |
                |
                |                 x    w1
                +----> mu         .-----------+
                                  |           |
                              w4  |           | w2
                                  |           |
                                  |           |
                                  .-----------.

        =#
        if iflag == 3 || iflag == 4
            substitute!(work1,U[μ])
            idir2= (μ,-ν)
            gauge_shift!(work2,idir2,U[ν])
            gauge_shift!(work3,-ν,U[μ])
            gauge_shift!(work4,-ν,U[ν])
            #println(work2[1,1,1])
            #exit()
        end
        #   ...  iflag = 3 (Case 4' and 2')
        if iflag == 3
            #=
              Case 4'
                   nu
                    |
                    |                 x    w1
                    +----> mu         o<----------.
                                      |           |
                                  w4  |           | w2
                                      |           |
                                      |           |
                                      .-----------.
                                           w3
              Case 2'
                   nu
                    |
                    |                 x    w1
                    +----> mu         .<----------.
                                      |           |
                                  w4  |           | w2
                                      |           |
                                      |           |
                                      .-----------o
                                           w3
            =#

            mul!(temp1,work4',work3)

            for ia=1:numbasis
                lambdamul(temp2,work1,ia,fparam.SUNgenerator) # temp2=(lambda_a/2)*work1
                mul!(temp3,work2,temp2')

                mul!(temp4,temp1,temp3)
                mul!(temp4,-im)
                cimaglnk!(temp4)
                substitute!(dFμν1[ia],coe,temp4)

                mul!(temp4,temp3,temp1)
                mul!(temp4,-im)
                cimaglnk!(temp4)
                substitute!(dFμν2[ia],coe,temp4)
            end
        end

        # ...  iflag = 4 (Case 3' and 1')
        if iflag == 4
            #=
              Case 3'
                   nu
                    |
                    |                 x    w1
                    +----> mu         .<----------o
                                      |           |
                                  w4  |           | w2
                                      |           |
                                      |           |
                                      .-----------.
                                           w3
              Case 1'
                   nu
                    |
                    |                 x    w1
                    +----> mu         .<----------.
                                      |           |
                                  w4  |           | w2
                                      |           |
                                      |           |
                                      o-----------.

            =#

            mul!(temp1,work3,work2)

            for ia=1:numbasis
                lambdamul(temp2,work1,ia,fparam.SUNgenerator) # temp2=(lambda_a/2)*work1
                mul!(temp3,work4,temp2)

                mul!(temp4,temp3',temp1)
                mul!(temp4,-im)
                cimaglnk!(temp4)
                #println(temp4.g)
                #exit()

                substitute!(dFμν1[ia],coe,temp4)

                mul!(temp4,temp1,temp3')
                mul!(temp4,-im)
                cimaglnk!(temp4)

                substitute!(dFμν2[ia],coe,temp4)
            end
        end
        return

    end

    """
        x <= x - x_aj
    """
    function cimaglnk!(x::SUNGaugeFields_1d)
        NV = x.NV
        NC = x.NC
        for i=1:NV
            for k2=1:NC
                for k1=k2:NC
                    xdiff = x[k1,k2,i] - conj(x[k2,k1,i])
                    x[k1,k2,i] = xdiff
                    if k1 != k2
                        x[k2,k1,i] = -conj(xdiff)
                    end
                end
            end
        end
    end

    function cimaglnk!(x::SU3GaugeFields_1d)
        NV = x.NV
        for i=1:NV
            z11 = x[1,1,i] - conj(x[1,1,i])
            z12 = x[1,2,i] - conj(x[2,1,i])
            z13 = x[1,3,i] - conj(x[3,1,i])

            z22 = x[2,2,i] - conj(x[2,2,i])
            z23 = x[2,3,i] - conj(x[3,2,i])

            z33 = x[3,3,i] - conj(x[3,3,i])

            x[1,1,i] = z11
            x[1,2,i] = z12
            x[1,3,i] = z13
            x[2,1,i] = -conj(z12)
            x[2,2,i] = z22
            x[2,3,i] = z23
            x[3,1,i] = -conj(z13)
            x[3,2,i] = -conj(z23)
            x[3,3,i] = z33
        end
    end

    function cimaglnk!(x::SU2GaugeFields_1d)
        NV = x.NV
        for i=1:NV
            z11 = x[1,1,i] - conj(x[1,1,i])
            z12 = x[1,2,i] - conj(x[2,1,i])

            z22 = x[2,2,i] - conj(x[2,2,i])

            x[1,1,i] = z11
            x[1,2,i] = z12
            x[2,1,i] = -conj(z12)
            x[2,2,i] = z22

        end
    end

end


text:
using sem, Arrow, ModelingToolkit,
    LinearAlgebra, SparseArrays, DataFrames,
    Optim, LineSearches, StatsBase, Distributions, Statistics, Random,
    Distances

Random.seed!(2848102)

################################################### read data ###################################################

cd("test")

## Observed Data
data = DataFrame(Arrow.Table("comparisons/data_dem.arrow"))
par_ml = DataFrame(Arrow.Table("comparisons/par_dem_ml.arrow"))
par_ls = DataFrame(Arrow.Table("comparisons/par_dem_ls.arrow"))

data =
    select(
        data,
        [:x1, :x2, :x3, :y1, :y2, :y3, :y4, :y5, :y6, :y7, :y8])

data = convert(Matrix{Float64}, data)

n = 100
n_obs = size(data, 1)
n_obs_gen = 50
max_iter = 200

################################################### draw bootstrap samples ###################################################

function bootstrap_samples(data, nobs, n_obs_gen, n)
    obs = 1:nobs
    rowind = [sample(obs, n_obs_gen) for i in 1:n]
    samples = [data[rowind[i], :] for i in 1:n]
    return samples
end

data_boot = bootstrap_samples(data, n_obs, n_obs_gen, n)

################################################### draw normal samples from observed distribution ###################################################

distribution_observed = MvNormal(vec(mean(data; dims = 1)), cov(data))
data_obs = [rand(distribution_observed, n_obs_gen) for i in 1:n]
data_obs = permutedims.(data_obs)

################################################### draw normal samples from implied distribution ###################################################
# observed
function get_imp_cov(x)
    S =[x[1]  0     0     0     0     0     0     0     0     0     0     0     0     0
    0     x[2]  0     0     0     0     0     0     0     0     0     0     0     0
    0     0     x[3]  0     0     0     0     0     0     0     0     0     0     0
    0     0     0     x[4]  0     0     0     x[15] 0     0     0     0     0     0
    0     0     0     0     x[5]  0     x[16] 0     x[17] 0     0     0     0     0
    0     0     0     0     0     x[6]  0     0     0     x[18] 0     0     0     0
    0     0     0     0     x[16] 0     x[7]  0     0     0     x[19] 0     0     0
    0     0     0     x[15] 0     0     0     x[8]  0     0     0     0     0     0
    0     0     0     0     x[17] 0     0     0     x[9]  0     x[20] 0     0     0
    0     0     0     0     0     x[18] 0     0     0     x[10] 0     0     0     0
    0     0     0     0     0     0     x[19] 0     x[20] 0     x[11] 0     0     0
    0     0     0     0     0     0     0     0     0     0     0     x[12] 0     0
    0     0     0     0     0     0     0     0     0     0     0     0     x[13] 0
    0     0     0     0     0     0     0     0     0     0     0     0     0     x[14]]

    F =[1.0 0 0 0 0 0 0 0 0 0 0 0 0 0
        0 1 0 0 0 0 0 0 0 0 0 0 0 0
        0 0 1 0 0 0 0 0 0 0 0 0 0 0
        0 0 0 1 0 0 0 0 0 0 0 0 0 0
        0 0 0 0 1 0 0 0 0 0 0 0 0 0
        0 0 0 0 0 1 0 0 0 0 0 0 0 0
        0 0 0 0 0 0 1 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 1 0 0 0 0 0 0
        0 0 0 0 0 0 0 0 1 0 0 0 0 0
        0 0 0 0 0 0 0 0 0 1 0 0 0 0
        0 0 0 0 0 0 0 0 0 0 1 0 0 0]

    A =[0  0  0  0  0  0  0  0  0  0  0     1     0     0
        0  0  0  0  0  0  0  0  0  0  0     x[21] 0     0
        0  0  0  0  0  0  0  0  0  0  0     x[22] 0     0
        0  0  0  0  0  0  0  0  0  0  0     0     1     0
        0  0  0  0  0  0  0  0  0  0  0     0     x[23] 0
        0  0  0  0  0  0  0  0  0  0  0     0     x[24] 0
        0  0  0  0  0  0  0  0  0  0  0     0     x[25] 0
        0  0  0  0  0  0  0  0  0  0  0     0     0     1
        0  0  0  0  0  0  0  0  0  0  0     0     0     x[26]
        0  0  0  0  0  0  0  0  0  0  0     0     0     x[27]
        0  0  0  0  0  0  0  0  0  0  0     0     0     x[28]
        0  0  0  0  0  0  0  0  0  0  0     0     0     0
        0  0  0  0  0  0  0  0  0  0  0     x[29] 0     0
        0  0  0  0  0  0  0  0  0  0  0     x[30] x[31] 0]
    invA = inv(I-A)
    imp_cov = F*invA*S*invA'*F'
    return imp_cov
end

par_order = [collect(21:34); collect(15:20); 2;3; 5;6;7; collect(9:14)]
imp_cov = get_imp_cov(convert(Vector{Float64}, par_ml.est[par_order]))

distribution_implied = MvNormal(vec(mean(data; dims = 1)), Symmetric(imp_cov))
data_imp = [rand(distribution_implied, n_obs_gen) for i in 1:n]
data_imp = permutedims.(data_imp)

############################## specify the models #########################
@ModelingToolkit.variables x[1:31]

S =[x[1]  0     0     0     0     0     0     0     0     0     0     0     0     0
    0     x[2]  0     0     0     0     0     0     0     0     0     0     0     0
    0     0     x[3]  0     0     0     0     0     0     0     0     0     0     0
    0     0     0     x[4]  0     0     0     x[15] 0     0     0     0     0     0
    0     0     0     0     x[5]  0     x[16] 0     x[17] 0     0     0     0     0
    0     0     0     0     0     x[6]  0     0     0     x[18] 0     0     0     0
    0     0     0     0     x[16] 0     x[7]  0     0     0     x[19] 0     0     0
    0     0     0     x[15] 0     0     0     x[8]  0     0     0     0     0     0
    0     0     0     0     x[17] 0     0     0     x[9]  0     x[20] 0     0     0
    0     0     0     0     0     x[18] 0     0     0     x[10] 0     0     0     0
    0     0     0     0     0     0     x[19] 0     x[20] 0     x[11] 0     0     0
    0     0     0     0     0     0     0     0     0     0     0     x[12] 0     0
    0     0     0     0     0     0     0     0     0     0     0     0     x[13] 0
    0     0     0     0     0     0     0     0     0     0     0     0     0     x[14]]

F =[1.0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 1 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 1 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 1 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 1 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 1 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 1 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 1 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 1 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 1 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 1 0 0 0]

A =[0  0  0  0  0  0  0  0  0  0  0     1     0     0
    0  0  0  0  0  0  0  0  0  0  0     x[21] 0     0
    0  0  0  0  0  0  0  0  0  0  0     x[22] 0     0
    0  0  0  0  0  0  0  0  0  0  0     0     1     0
    0  0  0  0  0  0  0  0  0  0  0     0     x[23] 0
    0  0  0  0  0  0  0  0  0  0  0     0     x[24] 0
    0  0  0  0  0  0  0  0  0  0  0     0     x[25] 0
    0  0  0  0  0  0  0  0  0  0  0     0     0     1
    0  0  0  0  0  0  0  0  0  0  0     0     0     x[26]
    0  0  0  0  0  0  0  0  0  0  0     0     0     x[27]
    0  0  0  0  0  0  0  0  0  0  0     0     0     x[28]
    0  0  0  0  0  0  0  0  0  0  0     0     0     0
    0  0  0  0  0  0  0  0  0  0  0     x[29] 0     0
    0  0  0  0  0  0  0  0  0  0  0     x[30] x[31] 0]

S = sparse(S)

#F
F = sparse(F)

#A
A = sparse(A)

par_order = [collect(21:34); collect(15:20); 2;3; 5;6;7; collect(9:14)]

start_val_ml = [fill(1.0, 14); fill(0, 6); fill(1, 8); fill(0, 3)]
start_val_ls = [fill(1.0, 14); fill(0, 6); fill(1, 8); fill(0, 3)]
start_val_snlls = [fill(1.0, 8); fill(0, 3)]
start_val_snlls_r = [fill(1.0, 8)-rand(8)/20; rand(3)/20]

function gen_model(data, opt_algo, start_val, A, S, F, lossfun, implyfun, parsubset; opt_options = nothing)
    # observed
    semobserved = SemObsCommon(data = data)
    # diff
    diff = SemFiniteDiff(opt_algo, opt_options)
    # loss
    loss = Loss([lossfun(semobserved, [0.0], similar(start_val))])
    # imply
    imply = implyfun(A, S, F, x[parsubset], start_val)
    # model
    model = Sem(semobserved, imply, loss, diff)
    return model
end

function check_approx(estimated, true_par; tol = 0.05)
    all(abs.(estimated .- true_par) .< tol*abs.(true))
end

function swap_observed(model, data, lossfun)
    observed = SemObsCommon(data = data)
    loss = Loss([lossfun(observed, [0.0], similar(model.imply.start_val))])
    new_model = Sem(deepcopy(observed), deepcopy(model.imply), deepcopy(loss), deepcopy(model.diff))
    return new_model
end

function com_dist_max(a, b)
    dist = maximum(abs.(a-b))
    return dist
end

function get_omega(mod, fit)
    omega = inv(mod.imply.G'*mod.loss.functions[1].V*mod.imply.G)*mod.imply.G'*mod.loss.functions[1].sᵀV'
    omega = omega[[1, 2, 3, 4, 6, 9, 11, 13, 14, 16, 17, 18, 19, 20, 5, 7, 8, 10, 12, 15]]
    return vcat(omega, fit.minimizer)
end

function analyze_fits(fits, models; true_pars = nothing, postcond = false)
    converged = Optim.converged.(fits)
    minimizer = Optim.minimizer.(fits)
    iterations = Optim.iterations.(fits)
    time = getfield.(fits, :time_run)
    if postcond
        minimizer = get_omega.(models, fits)
    end
    res = DataFrame(
        converged = converged,
        minimizer = minimizer,
        iterations = iterations,
        time = time)
    if !isnothing(true_pars)
        dist_euk = euclidean.(minimizer, [true_pars])
        res.dist_euk = dist_euk
        dist_max = com_dist_max.(minimizer, [true_pars])
        res.dist_max = dist_max
    end
    return res
end

function summarize_fits(fit_table)
    away = fit_table.dist_euk .> 20
    fit_table.converged[away] .= false
    n_converged = sum(fit_table.converged)
    mean_iter = mean(fit_table.iterations[fit_table.converged])
    mean_dist = mean(fit_table.dist_euk[fit_table.converged])
    mean_time = mean(fit_table.time[fit_table.converged])
    return DataFrame(n_converged = n_converged, mean_iter = mean_iter, mean_dist = mean_dist, mean_time = mean_time)
end

#=
mod = gen_model(data, BFGS(), start_val_snlls, A, S, F, sem.SemSWLS, sem.ImplySymbolicSWLS, 21:31; opt_options = Optim.Options())
fit = sem_fit(mod)
check_approx(get_omega(mod, fit), par_ls.est[par_order])
mod2 = swap_observed(mod, data, sem.SemSWLS)
fit2 = sem_fit(mod2)
check_approx(get_omega(mod2, fit), par_ls.est[par_order])

mod = gen_model(data, BFGS(), start_val_ml, A, S, F, sem.SemML, sem.ImplySymbolic, 1:31; opt_options = Optim.Options())
fit = sem_fit(mod)
check_approx(fit.minimizer, par_ml.est[par_order])
mod2 = swap_observed(mod, data, sem.SemML)
fit2 = sem_fit(mod2)
check_approx(fit.minimizer, par_ml.est[par_order])

mod = gen_model(data, BFGS(), start_val_ls, A, S, F, sem.SemWLS, sem.ImplySymbolicWLS, 1:31; opt_options = Optim.Options())
fit = sem_fit(mod)
check_approx(fit.minimizer, par_ls.est[par_order])
mod2 = swap_observed(mod, data, sem.SemWLS)
fit2 = sem_fit(mod2)
check_approx(fit.minimizer, par_ls.est[par_order])
 =#

function gen_model_vec(data_vec, opt_algo, start_val, A, S, F, lossfun, implyfun, parsubset; opt_options = nothing)
    model = gen_model(data_vec[1], opt_algo, start_val, A, S, F, lossfun, implyfun, parsubset; opt_options = opt_options)
    models = [swap_observed(model, data_vec[i], lossfun) for i in 1:size(data_vec, 1)]
    return models
end

########### get optimizers to work ###############
using BenchmarkTools

################### Optim.jl

mod = gen_model(data, BFGS(), start_val_ls, A, S, F, sem.SemML, sem.ImplySymbolic, 1:31; opt_options = Optim.Options())
fit = sem_fit(mod)
check_approx(fit.minimizer, par_ml.est[par_order])

mod = gen_model(data, BFGS(;linesearch = BackTracking(order=3), alphaguess = InitialHagerZhang()), start_val_ls, A, S, F, sem.SemML, sem.ImplySymbolic, 1:31; opt_options = Optim.Options())
fit = sem_fit(mod)
check_approx(fit.minimizer, par_ml.est[par_order])

mod = gen_model(data, LBFGS(;linesearch = BackTracking(order=3), alphaguess = InitialHagerZhang()), start_val_ls, A, S, F, sem.SemML, sem.ImplySymbolic, 1:31; opt_options = Optim.Options())
fit = sem_fit(mod)
check_approx(fit.minimizer, par_ml.est[par_order])

mod = gen_model(data, ParticleSwarm(;n_particles = 20), start_val_ls, A, S, F, sem.SemML, sem.ImplySymbolic, 1:31; opt_options = Optim.Options())
fit = sem_fit(mod)
check_approx(fit.minimizer, par_ml.est[par_order])

mod = gen_model(data, GradientDescent(), start_val_ls, A, S, F, sem.SemML, sem.ImplySymbolic, 1:31; opt_options = Optim.Options())
fit = sem_fit(mod)
check_approx(fit.minimizer, par_ml.est[par_order])
# does not work

mod = gen_model(data, ConjugateGradient(), start_val_ls, A, S, F, sem.SemML, sem.ImplySymbolic, 1:31; opt_options = Optim.Options())
fit = sem_fit(mod)
check_approx(fit.minimizer, par_ml.est[par_order])

mod = gen_model(data, Newton(), start_val_ls, A, S, F, sem.SemML, sem.ImplySymbolic, 1:31; opt_options = Optim.Options())
fit = sem_fit(mod)
check_approx(fit.minimizer, par_ml.est[par_order])

mod = gen_model(data, Newton(;linesearch = BackTracking(order=3), alphaguess = InitialHagerZhang()), start_val_ls, A, S, F, sem.SemML, sem.ImplySymbolic, 1:31; opt_options = Optim.Options())
fit = sem_fit(mod)
check_approx(fit.minimizer, par_ml.est[par_order])

mod = gen_model(data, NewtonTrustRegion(), start_val_ls, A, S, F, sem.SemML, sem.ImplySymbolic, 1:31; opt_options = Optim.Options())
fit = sem_fit(mod)
check_approx(fit.minimizer, par_ml.est[par_order])

####################### NLopt.jl
####################### Local, derivative-based

mod = gen_model(data, :LD_LBFGS, start_val_ls, A, S, F, sem.SemML, sem.ImplySymbolic, 1:31; opt_options = nothing)
fit = sem.sem_fit_nlopt(mod)
check_approx(fit[2], par_ml.est[par_order])

mod = gen_model(data, :LD_MMA, start_val_ls, A, S, F, sem.SemML, sem.ImplySymbolic, 1:31; opt_options = nothing)
fit = sem.sem_fit_nlopt(mod)
check_approx(fit[2], par_ml.est[par_order])

mod = gen_model(data, :LD_CCSAQ, start_val_ls, A, S, F, sem.SemML, sem.ImplySymbolic, 1:31; opt_options = nothing)
fit = sem.sem_fit_nlopt(mod)
check_approx(fit[2], par_ml.est[par_order])

mod = gen_model(data, :LD_SLSQP, start_val_ls, A, S, F, sem.SemML, sem.ImplySymbolic, 1:31; opt_options = nothing)
fit = sem.sem_fit_nlopt(mod)
check_approx(fit[2], par_ml.est[par_order])

mod = gen_model(data, :LD_TNEWTON_PRECOND_RESTART, start_val_ls, A, S, F, sem.SemML, sem.ImplySymbolic, 1:31; opt_options = nothing)
fit = sem.sem_fit_nlopt(mod)
check_approx(fit[2], par_ml.est[par_order])

mod = gen_model(data, :LD_TNEWTON_PRECOND, start_val_ls, A, S, F, sem.SemML, sem.ImplySymbolic, 1:31; opt_options = nothing)
fit = sem.sem_fit_nlopt(mod)
check_approx(fit[2], par_ml.est[par_order])

mod = gen_model(data, :LD_TNEWTON_RESTART, start_val_ls, A, S, F, sem.SemML, sem.ImplySymbolic, 1:31; opt_options = nothing)
fit = sem.sem_fit_nlopt(mod)
check_approx(fit[2], par_ml.est[par_order])

mod = gen_model(data, :LD_TNEWTON, start_val_ls, A, S, F, sem.SemML, sem.ImplySymbolic, 1:31; opt_options = nothing)
fit = sem.sem_fit_nlopt(mod)
check_approx(fit[2], par_ml.est[par_order])

mod = gen_model(data, :LD_VAR2, start_val_ls, A, S, F, sem.SemML, sem.ImplySymbolic, 1:31; opt_options = nothing)
fit = sem.sem_fit_nlopt(mod)
check_approx(fit[2], par_ml.est[par_order])

mod = gen_model(data, :LD_VAR1, start_val_ls, A, S, F, sem.SemML, sem.ImplySymbolic, 1:31; opt_options = nothing)
fit = sem.sem_fit_nlopt(mod)
check_approx(fit[2], par_ml.est[par_order])

####################### Global

mod = gen_model(data, :GN_CRS2_LM, start_val_ls, A, S, F, sem.SemML, sem.ImplySymbolic, 1:31; opt_options = nothing)
fit = sem.sem_fit_nlopt(mod; upper = fill(10.0, 31), lower = fill(-10.0, 31))#, local_algo = :LD_LBFGS)
check_approx(fit[2], par_ml.est[par_order])

####################### Local derivative free
mod = gen_model(data, :LN_BOBYQA, start_val_ls, A, S, F, sem.SemML, sem.ImplySymbolic, 1:31; opt_options = nothing)
fit = sem.sem_fit_nlopt(mod; xtol_rel = 1e-10, ftol_rel = 1e-12)
check_approx(fit[2], par_ml.est[par_order])
## --> does not converge

mod = gen_model(data, :LN_NEWUOA, start_val_ls, A, S, F, sem.SemML, sem.ImplySymbolic, 1:31; opt_options = nothing)
fit = sem.sem_fit_nlopt(mod; xtol_rel = 1e-10, maxeval = 10000)
check_approx(fit[2], par_ml.est[par_order])

mod = gen_model(data, :LN_PRAXIS, start_val_ls, A, S, F, sem.SemML, sem.ImplySymbolic, 1:31; opt_options = nothing)
fit = sem.sem_fit_nlopt(mod; maxeval = 10000)
check_approx(fit[2], par_ml.est[par_order])

mod = gen_model(data, :LN_NELDERMEAD, start_val_ls, A, S, F, sem.SemML, sem.ImplySymbolic, 1:31; opt_options = nothing)
fit = sem.sem_fit_nlopt(mod; maxeval = 200000, xtol_rel = 1e-40, ftol_rel = 1e-32)
check_approx(fit[2], par_ml.est[par_order])
## -> does not converge

mod = gen_model(data, :LN_SBPLX, start_val_ls, A, S, F, sem.SemML, sem.ImplySymbolic, 1:31; opt_options = nothing)
@benchmark fit = sem.sem_fit_nlopt(mod; maxeval = 100000)
check_approx(fit[2], par_ml.est[par_order])
########################################################## simulation implied ############################################################

models_imp_ls = gen_model_vec(data_imp, BFGS(), start_val_ls, A, S, F, sem.SemWLS, sem.ImplySymbolicWLS, 1:31; opt_options = Optim.Options(;iterations = max_iter))
fits_imp_ls = sem_fit.(models_imp_ls)
res_imp_ls = analyze_fits(fits_imp_ls, models_imp_ls; true_pars = par_ml.est[par_order])
summarize_fits(res_imp_ls)

models_imp_ml = gen_model_vec(data_imp, BFGS(), start_val_ml, A, S, F, SemML, ImplySymbolic, 1:31; opt_options = Optim.Options(;iterations = max_iter))
fits_imp_ml = sem_fit.(models_imp_ml)
res_imp_ml = analyze_fits(fits_imp_ml, models_imp_ml; true_pars = par_ml.est[par_order])
summarize_fits(res_imp_ml)

models_imp_snlls = gen_model_vec(data_imp, BFGS(), start_val_snlls, A, S, F, sem.SemSWLS, sem.ImplySymbolicSWLS, 21:31; opt_options = Optim.Options(;iterations = max_iter))
fits_imp_snlls = sem_fit.(models_imp_snlls)
res_imp_snlls = analyze_fits(fits_imp_snlls, models_imp_snlls; true_pars = par_ml.est[par_order], postcond = true)
summarize_fits(res_imp_snlls)

########################################################## simulation bootstrap ############################################################
models_imp_ls = gen_model_vec(data_boot,BFGS(;linesearch = BackTracking(order=3), alphaguess = InitialHagerZhang()), start_val_ls, A, S, F, sem.SemWLS, sem.ImplySymbolicWLS, 1:31; opt_options = Optim.Options(;iterations = max_iter))
fits_imp_ls = sem_fit.(models_imp_ls)
res_imp_ls = analyze_fits(fits_imp_ls, models_imp_ls; true_pars = par_ml.est[par_order])
summarize_fits(res_imp_ls)

models_imp_ml = gen_model_vec(data_boot, BFGS(;linesearch = BackTracking(order=3), alphaguess = InitialHagerZhang()), start_val_ml, A, S, F, SemML, ImplySymbolic, 1:31; opt_options = Optim.Options(;iterations = max_iter))
fits_imp_ml = sem_fit.(models_imp_ml)
res_imp_ml = analyze_fits(fits_imp_ml, models_imp_ml; true_pars = par_ml.est[par_order])
summarize_fits(res_imp_ml)

models_imp_snlls = gen_model_vec(data_boot, BFGS(;linesearch = BackTracking(order=3), alphaguess = InitialHagerZhang()), start_val_snlls, A, S, F, sem.SemSWLS, sem.ImplySymbolicSWLS, 21:31; opt_options = Optim.Options(;iterations = max_iter))
fits_imp_snlls = sem_fit.(models_imp_snlls)
res_imp_snlls = analyze_fits(fits_imp_snlls, models_imp_snlls; true_pars = par_ml.est[par_order], postcond = true)
summarize_fits(res_imp_snlls)

########################################################## simulation obs ############################################################
models_imp_ls = gen_model_vec(data_obs, BFGS(), start_val_ls, A, S, F, sem.SemWLS, sem.ImplySymbolicWLS, 1:31; opt_options = Optim.Options(;iterations = max_iter))
fits_imp_ls = sem_fit.(models_imp_ls)
res_imp_ls = analyze_fits(fits_imp_ls, models_imp_ls; true_pars = par_ml.est[par_order])
summarize_fits(res_imp_ls)

models_imp_ml = gen_model_vec(data_obs, BFGS(), start_val_ml, A, S, F, SemML, ImplySymbolic, 1:31; opt_options = Optim.Options(;iterations = max_iter))
fits_imp_ml = sem_fit.(models_imp_ml)
res_imp_ml = analyze_fits(fits_imp_ml, models_imp_ml; true_pars = par_ml.est[par_order])
summarize_fits(res_imp_ml)

models_imp_snlls = gen_model_vec(data_obs, BFGS(), start_val_snlls, A, S, F, sem.SemSWLS, sem.ImplySymbolicSWLS, 21:31; opt_options = Optim.Options(;iterations = max_iter))
fits_imp_snlls = sem_fit.(models_imp_snlls)
res_imp_snlls = analyze_fits(fits_imp_snlls, models_imp_snlls; true_pars = par_ml.est[par_order], postcond = true)
summarize_fits(res_imp_snlls)

########### speed up optimizer ###############
using BenchmarkTools

mod = gen_model(data, BFGS(), start_val_ls, A, S, F, sem.SemWLS, sem.ImplySymbolicWLS, 1:31; opt_options = Optim.Options())
fit = sem_fit(mod)
check_approx(fit.minimizer, par_ls.est[par_order])
### BFGS
@benchmark sem_fit(mod)

mod = gen_model(data, LBFGS(), start_val_ls, A, S, F, sem.SemWLS, sem.ImplySymbolicWLS, 1:31; opt_options = Optim.Options())
fit = sem_fit(mod)
check_approx(fit.minimizer, par_ls.est[par_order])
### LBFGS
@benchmark sem_fit(mod)

## --> BFGS wins (for this small model)

mod = gen_model(data, BFGS(;linesearch = BackTracking(order=3)), start_val_ls, A, S, F, sem.SemWLS, sem.ImplySymbolicWLS, 1:31; opt_options = Optim.Options())
fit = sem_fit(mod)
check_approx(fit.minimizer, par_ls.est[par_order])
### BFGS + Backtracking
@benchmark sem_fit(mod)

mod = gen_model(data, BFGS(;linesearch = HagerZhang(), alphaguess = InitialHagerZhang()), start_val_ls, A, S, F, sem.SemWLS, sem.ImplySymbolicWLS, 1:31; opt_options = Optim.Options())
fit = sem_fit(mod)
check_approx(fit.minimizer, par_ls.est[par_order])
### BFGS + HagerZang
@benchmark sem_fit(mod)

mod = gen_model(data, BFGS(;linesearch = BackTracking(order=3), alphaguess = InitialHagerZhang()), start_val_ls, A, S, F, sem.SemWLS, sem.ImplySymbolicWLS, 1:31; opt_options = Optim.Options())
fit = sem_fit(mod)
check_approx(fit.minimizer, par_ls.est[par_order])
### BFGS + Backtracking
@benchmark sem_fit(mod)

mod = gen_model(data, BFGS(;linesearch = MoreThuente()), start_val_ls, A, S, F, sem.SemWLS, sem.ImplySymbolicWLS, 1:31; opt_options = Optim.Options())
fit = sem_fit(mod)
check_approx(fit.minimizer, par_ls.est[par_order])
### BFGS + MoreThuente
@benchmark sem_fit(mod)

# --> similar performance

mod = gen_model(data, ParticleSwarm(;n_particles = 50), start_val_ls, A, S, F, sem.SemWLS, sem.ImplySymbolicWLS, 1:31; opt_options = Optim.Options())
fit = sem_fit(mod)
check_approx(fit.minimizer, par_ls.est[par_order])
### BFGS + Backtracking
@benchmark sem_fit(mod)

mod = gen_model(data, OACCEL(nlprecon=LBFGS(), wmax=5), start_val_ls, A, S, F, sem.SemWLS, sem.ImplySymbolicWLS, 1:31; opt_options = Optim.Options())
fit = sem_fit(mod)
check_approx(fit.minimizer, par_ls.est[par_order])
### BFGS + Backtracking
@benchmark sem_fit(mod)

################# NLopt
mod = gen_model(data, :LD_LBFGS, start_val_ls, A, S, F, sem.SemWLS, sem.ImplySymbolicWLS, 1:31; opt_options = nothing)
fit = sem.sem_fit_nlopt(mod)
check_approx(fit[2], par_ls.est[par_order])
### BFGS
@benchmark sem.sem_fit_nlopt(mod)

########################################################## varying sample sizes + bootstrap ############################################################
n_obs_gen = 40
n = 300
data_boot = bootstrap_samples(data, n_obs, n_obs_gen, n)

models_imp_ls = gen_model_vec(data_boot,BFGS(;linesearch = BackTracking(order=3), alphaguess = InitialHagerZhang()), start_val_ls, A, S, F, sem.SemWLS, sem.ImplySymbolicWLS, 1:31; opt_options = Optim.Options(;iterations = max_iter))
fits_imp_ls = sem_fit.(models_imp_ls)
res_imp_ls = analyze_fits(fits_imp_ls, models_imp_ls; true_pars = par_ml.est[par_order])
summarize_fits(res_imp_ls)

models_imp_ml = gen_model_vec(data_boot, BFGS(;linesearch = BackTracking(order=3), alphaguess = InitialHagerZhang()), start_val_ml, A, S, F, SemML, ImplySymbolic, 1:31; opt_options = Optim.Options(;iterations = max_iter))
fits_imp_ml = sem_fit.(models_imp_ml)
res_imp_ml = analyze_fits(fits_imp_ml, models_imp_ml; true_pars = par_ml.est[par_order])
summarize_fits(res_imp_ml)

models_imp_snlls = gen_model_vec(data_boot, BFGS(;linesearch = BackTracking(order=3), alphaguess = InitialHagerZhang()), start_val_snlls, A, S, F, sem.SemSWLS, sem.ImplySymbolicSWLS, 21:31; opt_options = Optim.Options(;iterations = max_iter))
fits_imp_snlls = sem_fit.(models_imp_snlls)
res_imp_snlls = analyze_fits(fits_imp_snlls, models_imp_snlls; true_pars = par_ml.est[par_order], postcond = true)
summarize_fits(res_imp_snlls)