text:
"All converter variables"
function variable_dc_converter(pm::_PM.AbstractPowerModel; kwargs...)
    variable_conv_tranformer_flow(pm; kwargs...)
    variable_conv_reactor_flow(pm; kwargs...)

    variable_converter_active_power(pm; kwargs...)
    variable_converter_reactive_power(pm; kwargs...)
    variable_acside_current(pm; kwargs...)
    variable_dcside_power(pm; kwargs...)
    variable_converter_firing_angle(pm; kwargs...)

    variable_converter_filter_voltage(pm; kwargs...)
    variable_converter_internal_voltage(pm; kwargs...)

    variable_converter_to_grid_active_power(pm; kwargs...)
    variable_converter_to_grid_reactive_power(pm; kwargs...)
end

function variable_conv_tranformer_flow(pm::_PM.AbstractPowerModel; kwargs...)
    variable_conv_transformer_active_power_to(pm; kwargs...)
    variable_conv_transformer_reactive_power_to(pm; kwargs...)
end

function variable_conv_reactor_flow(pm::_PM.AbstractPowerModel; kwargs...)
    variable_conv_reactor_active_power_from(pm; kwargs...)
    variable_conv_reactor_reactive_power_from(pm; kwargs...)
end

"variable: `pconv_ac[j]` for `j` in `convdc`"
function variable_converter_active_power(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    pc = _PM.var(pm, nw)[:pconv_ac] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_pconv_ac",
    start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc, i), "P_g", 1.0)
    )

    if bounded
        for (c, convdc) in _PM.ref(pm, nw, :convdc)
            JuMP.set_lower_bound(pc[c],  convdc["Pacmin"])
            JuMP.set_upper_bound(pc[c],  convdc["Pacmax"])
        end
    end

    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc, :pconv, _PM.ids(pm, nw, :convdc), pc)
end

"variable: `qconv_ac[j]` for `j` in `convdc`"
function variable_converter_reactive_power(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    qc = _PM.var(pm, nw)[:qconv_ac] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_qconv_ac",
    start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc, i), "Q_g", 1.0)
    )

    if bounded
        for (c, convdc) in _PM.ref(pm, nw, :convdc)
            JuMP.set_lower_bound(qc[c],  convdc["Qacmin"])
            JuMP.set_upper_bound(qc[c],  convdc["Qacmax"])
        end
    end

    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc, :qconv, _PM.ids(pm, nw, :convdc), qc)
end

"variable: `pconv_grid_ac_to[j]` for `j` in `convdc`"
function variable_conv_transformer_active_power_to(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    bigM = 2;
    ptfto = _PM.var(pm, nw)[:pconv_tf_to] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_pconv_tf_to",
    start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc, i), "P_g", 1.0)
    )
    if bounded
        for (c, convdc) in _PM.ref(pm, nw, :convdc)
            JuMP.set_lower_bound(ptfto[c],  -convdc["Pacrated"] * bigM)
            JuMP.set_upper_bound(ptfto[c],   convdc["Pacrated"] * bigM)
        end
    end

    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc, :ptf_to, _PM.ids(pm, nw, :convdc), ptfto)
end

"variable: `qconv_grid_ac_to[j]` for `j` in `convdc`"
function variable_conv_transformer_reactive_power_to(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    bigM = 2;
    qtfto = _PM.var(pm, nw)[:qconv_tf_to] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_qconv_tf_to",
    start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc, i), "Q_g", 1.0)
    )
    if bounded
        for (c, convdc) in _PM.ref(pm, nw, :convdc)
            JuMP.set_lower_bound(qtfto[c],  -convdc["Qacrated"] * bigM)
            JuMP.set_upper_bound(qtfto[c],   convdc["Qacrated"] * bigM)
        end
    end

    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc, :qtf_to, _PM.ids(pm, nw, :convdc), qtfto)
end

"variable: `pconv_pr_from[j]` for `j` in `convdc`"
function variable_conv_reactor_active_power_from(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    bigM = 2;
    pprfr = _PM.var(pm, nw)[:pconv_pr_fr] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_pconv_pr_fr",
    start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc, i), "P_g", 1.0)
    )
    if bounded
        for (c, convdc) in _PM.ref(pm, nw, :convdc)
            JuMP.set_lower_bound(pprfr[c],  -convdc["Pacrated"] * bigM)
            JuMP.set_upper_bound(pprfr[c],   convdc["Pacrated"] * bigM)
        end
    end

    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc, :ppr_fr, _PM.ids(pm, nw, :convdc), pprfr)
end

"variable: `qconv_pr_from[j]` for `j` in `convdc`"
function variable_conv_reactor_reactive_power_from(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    bigM = 2;
    qprfr = _PM.var(pm, nw)[:qconv_pr_fr] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_qconv_pr_fr",
    start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc, i), "Q_g", 1.0)
    )
    if bounded
        for (c, convdc) in _PM.ref(pm, nw, :convdc)
            JuMP.set_lower_bound(qprfr[c],  -convdc["Qacrated"] * bigM)
            JuMP.set_upper_bound(qprfr[c],   convdc["Qacrated"] * bigM)
        end
    end

    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc, :qpr_fr, _PM.ids(pm, nw, :convdc), qprfr)
end

"variable: `pconv_grid_ac[j]` for `j` in `convdc`"
function variable_converter_to_grid_active_power(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    bigM = 2;
    ptffr = _PM.var(pm, nw)[:pconv_tf_fr] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_pconv_tf_fr",
    start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc, i), "P_g", 1.0)
    )
    if bounded
        for (c, convdc) in _PM.ref(pm, nw, :convdc)
            JuMP.set_lower_bound(ptffr[c],  -convdc["Pacrated"] * bigM)
            JuMP.set_upper_bound(ptffr[c],   convdc["Pacrated"] * bigM)
        end
    end

    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc, :pgrid, _PM.ids(pm, nw, :convdc), ptffr)
end

"variable: `qconv_grid_ac[j]` for `j` in `convdc`"
function variable_converter_to_grid_reactive_power(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    bigM = 2;
    qtffr = _PM.var(pm, nw)[:qconv_tf_fr] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_qconv_tf_fr",
    start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc, i), "Q_g", 1.0)
    )
    if bounded
        for (c, convdc) in _PM.ref(pm, nw, :convdc)
            JuMP.set_lower_bound(qtffr[c],  -convdc["Qacrated"] * bigM)
            JuMP.set_upper_bound(qtffr[c],   convdc["Qacrated"] * bigM)
        end
    end

    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc, :qgrid, _PM.ids(pm, nw, :convdc), qtffr)
end

"variable: `pconv_dc[j]` for `j` in `convdc`"
function variable_dcside_power(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    bigM = 1.2; # to account for losses, maximum losses to be derived
    pcdc = _PM.var(pm, nw)[:pconv_dc] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_pconv_dc",
    start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc, i), "Pdcset", 1.0)
    )
    if bounded
        for (c, convdc) in _PM.ref(pm, nw, :convdc)
            JuMP.set_lower_bound(pcdc[c],  -convdc["Pacrated"] * bigM)
            JuMP.set_upper_bound(pcdc[c],   convdc["Pacrated"] * bigM)
        end
    end

    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc, :pdc, _PM.ids(pm, nw, :convdc), pcdc)
end

"variable: `pconv_dc[j]` for `j` in `convdc`"
function variable_converter_firing_angle(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    phic = _PM.var(pm, nw)[:phiconv] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_phiconv",
    start = acos(_PM.comp_start_value(_PM.ref(pm, nw, :convdc, i), "Pdcset", 1.0) / sqrt((_PM.comp_start_value(_PM.ref(pm, nw, :convdc, i), "Pacrated", 1.0))^2 + (_PM.comp_start_value(_PM.ref(pm, nw, :convdc, i), "Qacrated", 1.0))^2))
    )
    if bounded
        for (c, convdc) in _PM.ref(pm, nw, :convdc)
            JuMP.set_lower_bound(phic[c],  0)
            JuMP.set_upper_bound(phic[c],  pi)
        end
    end

    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc, :phi, _PM.ids(pm, nw, :convdc), phic)
end

"variable: `iconv_ac[j]` for `j` in `convdc`"
function variable_acside_current(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    ic = _PM.var(pm, nw)[:iconv_ac] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_iconv_ac",
    start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc, i), "P_g", 1.0)
    )
    if bounded
        for (c, convdc) in _PM.ref(pm, nw, :convdc)
            JuMP.set_lower_bound(ic[c],  0)
            JuMP.set_upper_bound(ic[c],  convdc["Imax"])
        end
    end

    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc, :iconv, _PM.ids(pm, nw, :convdc), ic)
end

"variable: `iconv_ac[j]` and `iconv_ac_sq[j]` for `j` in `convdc`"
function variable_acside_current(pm::_PM.AbstractWModels; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    ic = _PM.var(pm, nw)[:iconv_ac] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_iconv_ac",
    start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc, i), "P_g", 1.0)
    )
    icsq = _PM.var(pm, nw)[:iconv_ac_sq] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_iconv_ac_sq",
    start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc, i), "P_g", 1.0)
    )
    if bounded
        for (c, convdc) in _PM.ref(pm, nw, :convdc)
            JuMP.set_lower_bound(ic[c],  0)
            JuMP.set_upper_bound(ic[c],  convdc["Imax"])
            JuMP.set_lower_bound(icsq[c],  0)
            JuMP.set_upper_bound(icsq[c],  convdc["Imax"]^2)
        end
    end

    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc, :iconv_ac, _PM.ids(pm, nw, :convdc), ic)
    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc, :iconv_ac_sq, _PM.ids(pm, nw, :convdc), icsq)
end

"variable: `itf_sq[j]` for `j` in `convdc`"
function variable_conv_transformer_current_sqr(pm::_PM.AbstractWModels; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    bigM = 2; #TODO derive exact bound
    itfsq = _PM.var(pm, nw)[:itf_sq] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_itf_sq",
    start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc, i), "P_g", 1.0)^2
    )
    if bounded
        for (c, convdc) in _PM.ref(pm, nw, :convdc)
            JuMP.set_lower_bound(itfsq[c],  0)
            JuMP.set_upper_bound(itfsq[c], (bigM * convdc["Imax"])^2)
        end
    end
    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc, :itf_sq, _PM.ids(pm, nw, :convdc), itfsq)
end

"variable: `irc_sq[j]` for `j` in `convdc`"
function variable_conv_reactor_current_sqr(pm::_PM.AbstractWModels; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    bigM = 2; #TODO derive exact bound
    iprsq = _PM.var(pm, nw)[:irc_sq] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_irc_sq",
    start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc, i), "P_g", 1.0)^2
    )
    if bounded
        for (c, convdc) in _PM.ref(pm, nw, :convdc)
            JuMP.set_lower_bound(iprsq[c],  0)
            JuMP.set_upper_bound(iprsq[c], (bigM * convdc["Imax"])^2)
        end
    end
    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc, :ipr_sq, _PM.ids(pm, nw, :convdc), iprsq)
end

function variable_converter_filter_voltage(pm::_PM.AbstractPowerModel; kwargs...)
    variable_converter_filter_voltage_magnitude(pm; kwargs...)
    variable_converter_filter_voltage_angle(pm; kwargs...)
end

"variable: `vmf[j]` for `j` in `convdc`"
function variable_converter_filter_voltage_magnitude(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    bigM = 1.2; # only internal converter voltage is strictly regulated
    vmf = _PM.var(pm, nw)[:vmf] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_vmf",
    start = _PM.ref(pm, nw, :convdc, i, "Vtar")
    )
    if bounded
        for (c, convdc) in _PM.ref(pm, nw, :convdc)
            JuMP.set_lower_bound(vmf[c], convdc["Vmmin"] / bigM)
            JuMP.set_upper_bound(vmf[c], convdc["Vmmax"] * bigM)
        end
    end
    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc, :vmfilt, _PM.ids(pm, nw, :convdc), vmf)
end

"variable: `vaf[j]` for `j` in `convdc`"
function variable_converter_filter_voltage_angle(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    bigM = 2*pi; #
    vaf = _PM.var(pm, nw)[:vaf] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_vaf",
    start = 0
    )
    if bounded
        for (c, convdc) in _PM.ref(pm, nw, :convdc)
            JuMP.set_lower_bound(vaf[c], -bigM)
            JuMP.set_upper_bound(vaf[c],  bigM)
        end
    end
    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc, :vafilt, _PM.ids(pm, nw, :convdc), vaf)
end

function variable_converter_internal_voltage(pm::_PM.AbstractPowerModel; kwargs...)
    variable_converter_internal_voltage_magnitude(pm; kwargs...)
    variable_converter_internal_voltage_angle(pm; kwargs...)
end

"variable: `vmc[j]` for `j` in `convdc`"
function variable_converter_internal_voltage_magnitude(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    vmc = _PM.var(pm, nw)[:vmc] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_vmc",
    start = _PM.ref(pm, nw, :convdc, i, "Vtar")
    )
    if bounded
        for (c, convdc) in _PM.ref(pm, nw, :convdc)
            JuMP.set_lower_bound(vmc[c], convdc["Vmmin"])
            JuMP.set_upper_bound(vmc[c], convdc["Vmmax"])
        end
    end
    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc, :vmconv, _PM.ids(pm, nw, :convdc), vmc)
end

"variable: `vac[j]` for `j` in `convdc`"
function variable_converter_internal_voltage_angle(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    bigM = 2*pi; #
    vac = _PM.var(pm, nw)[:vac] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_vac",
    start = 0
    )
    if bounded
        for (c, convdc) in _PM.ref(pm, nw, :convdc)
            JuMP.set_lower_bound(vac[c], -bigM)
            JuMP.set_upper_bound(vac[c],  bigM)
        end
    end
    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc, :vaconv, _PM.ids(pm, nw, :convdc), vac)
end

"variable: `wrf_ac[j]` and `wif_ac`  for `j` in `convdc`"
function variable_converter_filter_voltage_cross_products(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    bigM = 1.2; # only internal converter voltage is strictly regulated
    wrfac = _PM.var(pm, nw)[:wrf_ac] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_wrf_ac",
    start = _PM.ref(pm, nw, :convdc, i, "Vtar")^2
    )
    wifac = _PM.var(pm, nw)[:wif_ac] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_wif_ac",
    start = _PM.ref(pm, nw, :convdc, i, "Vtar")^2
    )
    if bounded
        for (c, convdc) in _PM.ref(pm, nw, :convdc)
            JuMP.set_lower_bound(wrfac[c],  0)
            JuMP.set_upper_bound(wrfac[c],  (convdc["Vmmax"] * bigM)^2)
            JuMP.set_lower_bound(wifac[c], -(convdc["Vmmax"] * bigM)^2)
            JuMP.set_upper_bound(wifac[c],  (convdc["Vmmax"] * bigM)^2)
        end
    end

    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc, :wrfilt, _PM.ids(pm, nw, :convdc), wrfac)
    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc, :wifilt, _PM.ids(pm, nw, :convdc), wifac)
end

"variable: `wf_ac` for `j` in `convdc`"
function variable_converter_filter_voltage_magnitude_sqr(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    bigM = 1.2; # only internal converter voltage is strictly regulated
    wfac = _PM.var(pm, nw)[:wf_ac] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_wf_ac",
    start = _PM.ref(pm, nw, :convdc, i, "Vtar")^2
    )
    if bounded
        for (c, convdc) in _PM.ref(pm, nw, :convdc)
            JuMP.set_lower_bound(wfac[c], (convdc["Vmmin"] / bigM)^2)
            JuMP.set_upper_bound(wfac[c], (convdc["Vmmax"] * bigM)^2)
        end
    end
    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc, :wfilt, _PM.ids(pm, nw, :convdc), wfac)
end

"variable: `wrc_ac[j]` and `wic_ac[j]`  for `j` in `convdc`"
function variable_converter_internal_voltage_cross_products(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    bigM = 1.2; # only internal converter voltage is strictly regulated
    wrcac = _PM.var(pm, nw)[:wrc_ac] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_wrc_ac",
    start = _PM.ref(pm, nw, :convdc, i, "Vtar")^2
    )
    wicac = _PM.var(pm, nw)[:wic_ac] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_wic_ac",
    start = _PM.ref(pm, nw, :convdc, i, "Vtar")^2
    )
    if bounded
        for (c, convdc) in _PM.ref(pm, nw, :convdc)
            JuMP.set_lower_bound(wrcac[c],  0)
            JuMP.set_upper_bound(wrcac[c],  (convdc["Vmmax"] * bigM)^2)
            JuMP.set_lower_bound(wicac[c], -(convdc["Vmmax"] * bigM)^2)
            JuMP.set_upper_bound(wicac[c],  (convdc["Vmmax"] * bigM)^2)
        end
    end

    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc, :wrconv, _PM.ids(pm, nw, :convdc), wrcac)
    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc, :wiconv, _PM.ids(pm, nw, :convdc), wicac)
end

"variable: `wc_ac[j]` for `j` in `convdc`"
function variable_converter_internal_voltage_magnitude_sqr(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    wcac = _PM.var(pm, nw)[:wc_ac] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc)], base_name="$(nw)_wc_ac",
    start = _PM.ref(pm, nw, :convdc, i, "Vtar")^2
    )
    if bounded
        for (c, convdc) in _PM.ref(pm, nw, :convdc)
            JuMP.set_lower_bound(wcac[c], (convdc["Vmmin"])^2)
            JuMP.set_upper_bound(wcac[c], (convdc["Vmmax"])^2)
        end
    end
    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc, :wconv, _PM.ids(pm, nw, :convdc), wcac)
end

function variable_cos_voltage(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    #only for lpac
end

##################################### TNEP variables ##################################
"All converter variables"
function variable_dc_converter_ne(pm::_PM.AbstractPowerModel; kwargs...)
    variable_conv_tranformer_flow_ne(pm; kwargs...)
    variable_conv_reactor_flow_ne(pm; kwargs...)
    variable_converter_ne(pm; kwargs...)

    variable_converter_active_power_ne(pm; kwargs...)
    variable_converter_reactive_power_ne(pm; kwargs...)
    variable_acside_current_ne(pm; kwargs...)
    variable_dcside_power_ne(pm; kwargs...)
    # variable_converter_firing_angle_ne(pm; kwargs...)

    variable_converter_filter_voltage_ne(pm; kwargs...)
    variable_converter_internal_voltage_ne(pm; kwargs...)
    #
    variable_converter_to_grid_active_power_ne(pm; kwargs...)
    variable_converter_to_grid_reactive_power_ne(pm; kwargs...)
end

function variable_conv_tranformer_flow_ne(pm::_PM.AbstractPowerModel; kwargs...)
    variable_conv_transformer_active_power_to_ne(pm; kwargs...)
    variable_conv_transformer_reactive_power_to_ne(pm; kwargs...)
end

function variable_conv_reactor_flow_ne(pm::_PM.AbstractPowerModel; kwargs...)
    variable_conv_reactor_active_power_from_ne(pm; kwargs...)
    variable_conv_reactor_reactive_power_from_ne(pm; kwargs...)
end

"variable: `0 <= convdc_ne[c] <= 1` for `c` in `candidate converters"
function variable_converter_ne(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, relax::Bool=false, report::Bool=true)
    if !relax
        Z_dc_conv_ne = _PM.var(pm, nw)[:conv_ne] = JuMP.@variable(pm.model,
        [i in _PM.ids(pm, nw, :convdc_ne)], base_name="$(nw)_conv_ne",
        binary = true,
        start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc_ne, i), "branchdc_tnep_start", 1.0)
        )
    else
        Z_dc_conv_ne = _PM.var(pm, nw)[:conv_ne] = JuMP.@variable(pm.model,
        [i in _PM.ids(pm, nw, :convdc_ne)], base_name="$(nw)_conv_ne",
        lower_bound = 0,
        upper_bound = 1,
        start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc_ne, i), "branchdc_tnep_start", 1.0)
        )
    end
    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc_ne, :isbuilt, _PM.ids(pm, nw, :convdc_ne), Z_dc_conv_ne)
 end

"variable: `pconv_ac_ne[j]` for `j` in `candidate convdc`"
function variable_converter_active_power_ne(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    pc_ne = _PM.var(pm, nw)[:pconv_ac_ne] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc_ne)], base_name="$(nw)_pconv_ac_ne",
    start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc_ne, i), "P_g", 1.0)
    )
    if bounded
        for (c, convdc) in _PM.ref(pm, nw, :convdc_ne)
            JuMP.set_lower_bound(pc_ne[c],  convdc["Pacmin"])
            JuMP.set_upper_bound(pc_ne[c],  convdc["Pacmax"])
        end
    end
    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc_ne, :pconv, _PM.ids(pm, nw, :convdc_ne), pc_ne)
 end

"variable: `qconv_ac_ne[j]` for `j` in `candidate convdc`"
function variable_converter_reactive_power_ne(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    qc_ne =  _PM.var(pm, nw)[:qconv_ac_ne] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc_ne)], base_name="$(nw)_qconv_ac_ne",
    start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc_ne, i), "Q_g", 1.0)
    )
    if bounded
        for (c, convdc) in _PM.ref(pm, nw, :convdc_ne)
            JuMP.set_lower_bound(qc_ne[c],  convdc["Qacmin"])
            JuMP.set_upper_bound(qc_ne[c],  convdc["Qacmax"])
        end
    end
    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc_ne, :qconv, _PM.ids(pm, nw, :convdc_ne), qc_ne)
end

"variable: `pconv_grid_ac_to_ne[j]` for `j` in `candidate convdc`"
function variable_conv_transformer_active_power_to_ne(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    bigM = 2;
    ptfto_ne = _PM.var(pm, nw)[:pconv_tf_to_ne] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc_ne)], base_name="$(nw)_pconv_tf_to_ne",
    start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc_ne, i), "P_g", 1.0)
    )
    if bounded
        for (c, convdc) in _PM.ref(pm, nw, :convdc_ne)
            JuMP.set_lower_bound(ptfto_ne[c],  -convdc["Pacrated"] * bigM)
            JuMP.set_upper_bound(ptfto_ne[c],   convdc["Pacrated"] * bigM)
        end
    end
    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc_ne, :ptf_to, _PM.ids(pm, nw, :convdc_ne), ptfto_ne)
 end

"variable: `qconv_grid_ac_to_ne[j]` for `j` in `candidate convdc`"
function variable_conv_transformer_reactive_power_to_ne(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    bigM = 2;
    qtfto_ne = _PM.var(pm, nw)[:qconv_tf_to_ne] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc_ne)], base_name="$(nw)_qconv_tf_to_ne",
    start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc_ne, i), "Q_g", 1.0)
    )
    if bounded
        for (c, convdc) in _PM.ref(pm, nw, :convdc_ne)
            JuMP.set_lower_bound(qtfto_ne[c],  -convdc["Qacrated"] * bigM)
            JuMP.set_upper_bound(qtfto_ne[c],   convdc["Qacrated"] * bigM)
        end
    end
    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc_ne, :qtf_to, _PM.ids(pm, nw, :convdc_ne), qtfto_ne)
end

"variable: `pconv_pr_from_ne[j]` for `j` in `candidate convdc`"
function variable_conv_reactor_active_power_from_ne(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    bigM = 2
    pprfr_ne = _PM.var(pm, nw)[:pconv_pr_fr_ne] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc_ne)], base_name="$(nw)_pconv_pr_from_ne",
    start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc_ne, i), "P_g", 1.0)
    )
    if bounded
        for (c, convdc) in _PM.ref(pm, nw, :convdc_ne)
            JuMP.set_lower_bound(pprfr_ne[c],  -convdc["Pacrated"] * bigM)
            JuMP.set_upper_bound(pprfr_ne[c],   convdc["Pacrated"] * bigM)
        end
    end
    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc_ne, :ppr_fr, _PM.ids(pm, nw, :convdc_ne), pprfr_ne)
end

"variable: `qconv_pr_from_ne[j]` for `j` in `candidate convdc`"
function variable_conv_reactor_reactive_power_from_ne(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    bigM = 2
    qprfr_ne = _PM.var(pm, nw)[:qconv_pr_fr_ne] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc_ne)], base_name="$(nw)_qconv_pr_from_ne",
    start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc_ne, i), "Q_g", 1.0)
    )
    if bounded
       for (c, convdc) in _PM.ref(pm, nw, :convdc_ne)
           JuMP.set_lower_bound(qprfr_ne[c],  -convdc["Qacrated"] * bigM)
           JuMP.set_upper_bound(qprfr_ne[c],   convdc["Qacrated"] * bigM)
       end
   end
   report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc_ne, :qpr_fr, _PM.ids(pm, nw, :convdc_ne), qprfr_ne)
end

"variable: `pconv_grid_ac_ne[j]` for `j` in `candidate convdc`"
function variable_converter_to_grid_active_power_ne(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    bigM = 2
    ptffr_ne = _PM.var(pm, nw)[:pconv_tf_fr_ne] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc_ne)], base_name="$(nw)_pconv_tf_fr_ne",
    start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc_ne, i), "P_g", 1.0)
    )
    if bounded
       for (c, convdc) in _PM.ref(pm, nw, :convdc_ne)
           JuMP.set_lower_bound(ptffr_ne[c],  -convdc["Pacrated"] * bigM)
           JuMP.set_upper_bound(ptffr_ne[c],   convdc["Pacrated"] * bigM)
       end
   end
   report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc_ne, :pgrid, _PM.ids(pm, nw, :convdc_ne), ptffr_ne)
end

"variable: `qconv_grid_ac_ne[j]` for `j` in `candidate convdc`"
function variable_converter_to_grid_reactive_power_ne(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    bigM = 2
    qtffr_ne = _PM.var(pm, nw)[:qconv_tf_fr_ne] = JuMP.@variable(pm.model,
        [i in _PM.ids(pm, nw, :convdc_ne)], base_name="$(nw)_qconv_tf_fr_ne",
        start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc_ne, i), "Q_g", 1.0)
        )
        if bounded
           for (c, convdc) in _PM.ref(pm, nw, :convdc_ne)
               JuMP.set_lower_bound(qtffr_ne[c],  -convdc["Qacrated"] * bigM)
               JuMP.set_upper_bound(qtffr_ne[c],   convdc["Qacrated"] * bigM)
           end
       end
       report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc_ne, :qgrid, _PM.ids(pm, nw, :convdc_ne), qtffr_ne)
end

"variable: `pconv_dc_ne[j]` for `j` in `candidate convdc`"
function variable_dcside_power_ne(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    bigM = 1.2; # to account for losses, maximum losses to be derived
    pcdc_ne = _PM.var(pm, nw)[:pconv_dc_ne] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc_ne)], base_name="$(nw)_pconv_dc_ne",
    start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc_ne, i), "Pdcset", 1.0)
    )
    if bounded
        for (c, convdc) in _PM.ref(pm, nw, :convdc_ne)
            JuMP.set_lower_bound(pcdc_ne[c],  -convdc["Pacrated"] * bigM)
            JuMP.set_upper_bound(pcdc_ne[c],   convdc["Pacrated"] * bigM)
        end
    end
    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc_ne, :pdc, _PM.ids(pm, nw, :convdc_ne), pcdc_ne)
 end

"variable: `pconv_dc_ne[j]` for `j` in `candidate convdc`"
function variable_converter_firing_angle_ne(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    phic_ne = _PM.var(pm, nw)[:phiconv_ne] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc_ne)], base_name="$(nw)_phiconv",
    lower_bound = 0,
    upper_bound =  pi,
    start = acos(_PM.comp_start_value(_PM.ref(pm, nw, :convdc_ne, i), "Pdcset", 1.0) / sqrt((_PM.comp_start_value(_PM.ref(pm, nw, :convdc_ne, i), "Pacrated", 1.0))^2 + (_PM.comp_start_value(_PM.ref(pm, nw, :convdc_ne, i), "Qacrated", 1.0))^2))
    )
    if bounded
        for (c, convdc) in _PM.ref(pm, nw, :convdc_ne)
            JuMP.set_lower_bound(phic_ne[c],  0)
            JuMP.set_upper_bound(phic_ne[c],  pi)
        end
    end
    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc_ne, :phi, _PM.ids(pm, nw, :convdc_ne), phic_ne)
end

"variable: `iconv_ac_ne[j]` for `j` in `candidate convdc`"
function variable_acside_current_ne(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    ic_ne = _PM.var(pm, nw)[:iconv_ac_ne] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc_ne)], base_name="$(nw)_iconv_ac_ne",
    start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc_ne, i), "P_g", 1.0)
    )
    if bounded
        for (c, convdc) in _PM.ref(pm, nw, :convdc_ne)
            JuMP.set_lower_bound(ic_ne[c],  0)
            JuMP.set_upper_bound(ic_ne[c],  convdc["Imax"])
        end
    end

    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc_ne, :iconv, _PM.ids(pm, nw, :convdc_ne), ic_ne)
 end

"variable: `iconv_ac_ne[j]` and `iconv_ac_sq_ne[j]` for `j` in `candidate convdc`"
function variable_acside_current_ne(pm::_PM.AbstractWModels; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    ic_ne = _PM.var(pm, nw)[:iconv_ac_ne] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc_ne)], base_name="$(nw)_iconv_ac_ne",
    start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc_ne, i), "P_g", 1.0)
    )
    icsq_ne = _PM.var(pm, nw)[:iconv_ac_sq_ne] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc_ne)], base_name="$(nw)_iconv_ac_sq_ne",
    start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc_ne, i), "P_g", 1.0)
    )
    if bounded
    for (c, convdc) in _PM.ref(pm, nw, :convdc_ne)
        JuMP.set_lower_bound(ic_ne[c],  0)
        JuMP.set_upper_bound(ic_ne[c],  convdc["Imax"])
        JuMP.set_lower_bound(icsq_ne[c],  0)
        JuMP.set_upper_bound(icsq_ne[c],  convdc["Imax"]^2)
    end
end
report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc_ne, :iconv_ac, _PM.ids(pm, nw, :convdc_ne), ic_ne)
report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc_ne, :iconv_ac_sq, _PM.ids(pm, nw, :convdc_ne), icsq_ne)
end

"variable: `itf_sq_ne[j]` for `j` in `candidate convdc`"
function variable_conv_transformer_current_sqr_ne(pm::_PM.AbstractWModels;  nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    bigM = 2 #TODO derive exact bound
    itfsq_ne = _PM.var(pm, nw)[:itf_sq_ne] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc_ne)], base_name="$(nw)_itf_sq_ne",
    start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc_ne, i), "P_g", 1.0)^2
    )
    if bounded
        for (c, convdc) in _PM.ref(pm, nw, :convdc_ne)
            JuMP.set_lower_bound(itfsq_ne[c],  0)
            JuMP.set_upper_bound(itfsq_ne[c], (bigM * convdc["Imax"])^2)
        end
    end
report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc_ne, :itf_sq, _PM.ids(pm, nw, :convdc_ne), itfsq_ne)
end

"variable: `irc_sq_ne[j]` for `j` in `candidate convdc`"
function variable_conv_reactor_current_sqr_ne(pm::_PM.AbstractWModels;  nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    bigM = 2 #TODO derive exact bound
    iprsq_ne = _PM.var(pm, nw)[:irc_sq_ne] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc_ne)], base_name="$(nw)_irc_sq_ne",
    start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc_ne, i), "P_g", 1.0)^2
    )
    if bounded
        for (c, convdc) in _PM.ref(pm, nw, :convdc_ne)
            JuMP.set_lower_bound(iprsq_ne[c],  0)
            JuMP.set_upper_bound(iprsq_ne[c], (bigM * convdc["Imax"])^2)
        end
    end
report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc_ne, :ipr_sq, _PM.ids(pm, nw, :convdc_ne), iprsq_ne)
end

function variable_converter_filter_voltage_ne(pm::_PM.AbstractPowerModel; kwargs...)
    variable_converter_filter_voltage_magnitude_ne(pm; kwargs...)
    variable_converter_filter_voltage_angle_ne(pm; kwargs...)
end

"variable: `vmf_ne[j]` for `j` in `candidate convdc`"
function variable_converter_filter_voltage_magnitude_ne(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    bigM = 1.2; # only internal converter voltage is strictly regulated
    vmf_ne = _PM.var(pm, nw)[:vmf_ne] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc_ne)], base_name="$(nw)_vmf_ne",
    start = _PM.ref(pm, nw, :convdc_ne, i, "Vtar")
    )
    if bounded
        for (c, convdc) in _PM.ref(pm, nw, :convdc_ne)
            JuMP.set_lower_bound(vmf_ne[c], convdc["Vmmin"] / bigM)
            JuMP.set_upper_bound(vmf_ne[c], convdc["Vmmax"] * bigM)
        end
    end
    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc_ne, :vmfilt, _PM.ids(pm, nw, :convdc_ne), vmf_ne)
end

"variable: `vaf_ne[j]` for `j` in `candidate convdc`"
function variable_converter_filter_voltage_angle_ne(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    bigM = 2*pi; #
    vaf_ne = _PM.var(pm, nw)[:vaf_ne] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc_ne)], base_name="$(nw)_vaf_ne",
    start = 0
    )
    if bounded
        for (c, convdc) in _PM.ref(pm, nw, :convdc_ne)
            JuMP.set_lower_bound(vaf_ne[c], -bigM)
            JuMP.set_upper_bound(vaf_ne[c],  bigM)
        end
    end
    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc_ne, :vafilt, _PM.ids(pm, nw, :convdc_ne), vaf_ne)
end

function variable_converter_internal_voltage_ne(pm::_PM.AbstractPowerModel; kwargs...)
    variable_converter_internal_voltage_magnitude_ne(pm; kwargs...)
    variable_converter_internal_voltage_angle_ne(pm; kwargs...)
end

"variable: `vmc_ne[j]` for `j` in `candidate convdc`"
function variable_converter_internal_voltage_magnitude_ne(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    vmc_ne = _PM.var(pm, nw)[:vmc_ne] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc_ne)], base_name="$(nw)_vmc_ne",
    start = _PM.ref(pm, nw, :convdc_ne, i, "Vtar")
    )
    if bounded
        for (c, convdc) in _PM.ref(pm, nw, :convdc)
            JuMP.set_lower_bound(vmc_ne[c], convdc["Vmmin"])
            JuMP.set_upper_bound(vmc_ne[c], convdc["Vmmax"])
        end
    end
    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc_ne, :vmconv, _PM.ids(pm, nw, :convdc_ne), vmc_ne)
end

"variable: `vac_ne[j]` for `j` in `candidate convdc`"
function variable_converter_internal_voltage_angle_ne(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    bigM = 2*pi; #
    vac_ne = _PM.var(pm, nw)[:vac_ne] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc_ne)], base_name="$(nw)_vac_ne",
    start = 0
    )
    if bounded
        for (c, convdc) in _PM.ref(pm, nw, :convdc_ne)
            JuMP.set_lower_bound(vac_ne[c], -bigM)
            JuMP.set_upper_bound(vac_ne[c],  bigM)
        end
    end
    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc_ne, :vaconv, _PM.ids(pm, nw, :convdc_ne), vac_ne)
end

"variable: `wrf_ac_ne[j]` and `wif_ac`  for `j` in `candidate convdc`"
function variable_converter_filter_voltage_cross_products_ne(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    bigM = 1.2; # only internal converter voltage is strictly regulated
    wrfac_ne = _PM.var(pm, nw)[:wrf_ac_ne] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc_ne)], base_name="$(nw)_wrf_ac_ne",
    start = _PM.ref(pm, nw, :convdc_ne, i, "Vtar")^2
    )
    wifac_ne = _PM.var(pm, nw)[:wif_ac_ne] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc_ne)], base_name="$(nw)_wif_ac_ne",
    start = _PM.ref(pm, nw, :convdc_ne, i, "Vtar")^2
    )
    if bounded
        for (c, convdc) in _PM.ref(pm, nw, :convdc_ne)
            JuMP.set_lower_bound(wrfac_ne[c],  0)
            JuMP.set_upper_bound(wrfac_ne[c],  (convdc["Vmmax"] * bigM)^2)
            JuMP.set_lower_bound(wifac_ne[c], -(convdc["Vmmax"] * bigM)^2)
            JuMP.set_upper_bound(wifac_ne[c],  (convdc["Vmmax"] * bigM)^2)
        end
    end
    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc_ne, :wrfilt, _PM.ids(pm, nw, :convdc_ne), wrfac_ne)
    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc_ne, :wifilt, _PM.ids(pm, nw, :convdc_ne), wifac_ne)
end

"variable: `wf_ac_ne` for `j` in `candidate convdc`"
function variable_converter_filter_voltage_magnitude_sqr_ne(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    bigM = 1.2; # only internal converter voltage is strictly regulated
    wfac_ne = _PM.var(pm, nw)[:wf_ac_ne] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc_ne)], base_name="$(nw)_wf_ac_ne",
    start = _PM.ref(pm, nw, :convdc_ne, i, "Vtar")^2
    )
    if bounded
        for (c, convdc) in _PM.ref(pm, nw, :convdc_ne)
            JuMP.set_lower_bound(wfac_ne[c], 0) # check
            # JuMP.set_lower_bound(wfac_ne[c], (convdc["Vmmin"] / bigM)^2)
            JuMP.set_upper_bound(wfac_ne[c], (convdc["Vmmax"] * bigM)^2)
        end
    end
    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc_ne, :wfilt, _PM.ids(pm, nw, :convdc_ne), wfac_ne)
end

"variable: `wrc_ac_ne[j]` and `wic_ac_ne[j]`  for `j` in `candidate convdc`"
function variable_converter_internal_voltage_cross_products_ne(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    bigM = 1.2; # only internal converter voltage is strictly regulated
    wrcac_ne = _PM.var(pm, nw)[:wrc_ac_ne] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc_ne)], base_name="$(nw)_wrc_ac_ne",
    start = _PM.ref(pm, nw, :convdc_ne, i, "Vtar")^2
    )
    wicac_ne = _PM.var(pm, nw)[:wic_ac_ne] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc_ne)], base_name="$(nw)_wic_ac_ne",
    start = _PM.ref(pm, nw, :convdc_ne, i, "Vtar")^2
    )
    if bounded
        for (c, convdc) in _PM.ref(pm, nw, :convdc_ne)
            JuMP.set_lower_bound(wrcac_ne[c],  0)
            JuMP.set_upper_bound(wrcac_ne[c],  (convdc["Vmmax"] * bigM)^2)
            JuMP.set_lower_bound(wicac_ne[c], -(convdc["Vmmax"] * bigM)^2)
            JuMP.set_upper_bound(wicac_ne[c],  (convdc["Vmmax"] * bigM)^2)
        end
    end

    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc_ne, :wrconv, _PM.ids(pm, nw, :convdc_ne), wrcac_ne)
    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc_ne, :wiconv, _PM.ids(pm, nw, :convdc_ne), wicac_ne)

end

"variable: `wc_ac_ne[j]` for `j` in `candidate convdc`"
function variable_converter_internal_voltage_magnitude_sqr_ne(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    wcac_ne = _PM.var(pm, nw)[:wc_ac_ne] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc_ne)], base_name="$(nw)_wc_ac_ne",
    start = _PM.ref(pm, nw, :convdc_ne, i, "Vtar")^2
    )
    if bounded
        for (c, convdc) in _PM.ref(pm, nw, :convdc_ne)
            JuMP.set_lower_bound(wcac_ne[c], 0) # check
            # JuMP.set_lower_bound(wcac[c], (convdc["Vmmin"])^2)
            JuMP.set_upper_bound(wcac_ne[c], (convdc["Vmmax"])^2)
        end
    end
    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc_ne, :wconv, _PM.ids(pm, nw, :convdc_ne), wcac_ne)
end

function variable_voltage_slack(pm::_PM.AbstractWModels; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=false)
    w_du = _PM.var(pm, nw)[:w_du] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :convdc_ne)], base_name="$(nw)_w_du",
    start = 0,  # check
    )
    if bounded
        for (c, convdc) in _PM.ref(pm, nw, :convdc_ne)
            JuMP.set_lower_bound(w_du[c], 0) # check
            JuMP.set_upper_bound(w_du[c], 2) # check
        end
    end
    report && _IM.sol_component_value(pm, _PM.pm_it_sym, nw, :convdc_ne, :wdu, _PM.ids(pm, nw, :convdc_ne), w_du)
end

function variable_cos_voltage_ne(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, bounded::Bool = true, report::Bool=true)
    #only for lpac
end


text:
import accrue.cryptoerase.runtime.Condition;

class C {
    static final Condition c = new Condition();
    int{L /c H} x = 0;
    int{L /c H} y = 0;

    public static void main(String[] args) {
	C o1 = new C();
	C o2 = new C();
	o1.x = 42;
	o2.y = o1.x;
    }
}

ext: jl

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase


text:
using Jecco, Jecco.AdS5_3_1
using Plots
pyplot()

dirname = "/home/mzilhao/scratch/Jecco/spinodal1D_x_model1_en0.9_run02/"

en_ts    = VEVTimeSeries(dirname, :energy)
xi_ts    = XiTimeSeries(dirname)
B1_c2_ts = BulkTimeSeries(dirname, :B1, 2)

xi = xi_ts[1500, :,3]
t, x, y = get_coords(xi_ts, 1500,:,3)

xi_plot = plot(
    x, xi,
    xlabel = "x",
    ylabel = "ξ",
    title  = "t = $t",
    marker = :hex,
    reuse=false,
)

en = en_ts[1:10:1500,:,3]
t, x, y = get_coords(en_ts, 1:10:1500,:,3)

en_plot = wireframe(
    t, x, en, transpose=true,
    xlabel = "t",
    ylabel = "x",
    title  = "energy",
    reuse=false,
)

# now using the infrastructure from AdS5_3_1/recipes.jl
B1_plot = wireframe(
    B1_c2_ts, 300, :, :, 3,
    xlabel = "u",
    ylabel = "x",
    title  = "B1 c=2",
    reuse=false,
)

display(xi_plot)
display(en_plot)
display(B1_plot)


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

baremodule Lib

using Base
import Base: call, (==), copy

using ..cartesianmapreduce, ..cartesianarray, ..sum, ..(.*), ..map, ..map!

type MT{T}
  val :: T
  idx :: Int
end

function ==(x::MT, y::MT)
  x.val == y.val && x.idx == y.idx
end

function copy(x::MT)
  MT(x.val, x.idx)
end

@inline function indmin{T<:Number}(A::DenseArray{T})
  m = MT(A[1], 1)
  cartesianmapreduce((length(A),),
    (x -> begin
            if m.val > x.val
              m.val = x.val
              m.idx = x.idx
            end
            return m
          end,
     m)) do i
    if A[i] < m.val
      m.val = A[i]
      m.idx = i
    end
    0
  end
  return m.idx
end

@inline function indmax{T<:Number}(A::DenseArray{T})
  m = MT(A[1], 1)
  cartesianmapreduce((length(A),),
    (x -> begin
            if m.val < x.val
              m.val = x.val
              m.idx = x.idx
            end
            return m
          end,
     m)) do i
    if A[i] > m.val
      m.val = A[i]
      m.idx = i
    end
    0
  end
  return m.idx
end

@inline function sumabs2(A::DenseArray)
  sum(A .* A)
end

@inline function diag(A::DenseMatrix)
  d::Int = min(size(A, 1), size(A, 2))
  cartesianarray(Tuple{eltype(A)}, (d,)) do i
     A[i, i]
  end
end

@inline function diagm(A::DenseVector)
  d::Int = size(A, 1)
  cartesianarray(Tuple{eltype(A)}, (d, d)) do i, j
    # the assignment below is a hack to avoid mutiple return in body
    v = i == j ? A[i] : zero(eltype(A))
  end
end

@inline function trace(A::DenseMatrix)
  sum(diag(A))
end

@inline function scale(A::DenseMatrix, b::DenseVector)
  m::Int, n::Int = size(A)
  cartesianarray(Tuple{eltype(A)}, (m, n)) do i, j
    A[i,j] * b[j]
  end
end

@inline function scale(b::DenseVector, A::DenseMatrix)
  m::Int, n::Int = size(A)
  cartesianarray(Tuple{eltype(A)}, (m, n)) do i, j
    b[i] * A[i,j]
  end
end

@inline function eye(m::Int, n::Int)
  cartesianarray(Tuple{Float64}, (m, n)) do i, j
    # the assignment below is a hack to avoid mutiple return in body
    v = i == j ? 1.0 : 0.0
  end
end

@inline function eye(m::Int)
  eye(m, m)
end

@inline function eye(A::DenseMatrix)
  eye(size(A, 1), size(A, 2))
end

@inline function repmat{T<:Number}(A::Array{T, 3}, m::Int, n::Int, l::Int)
  s::Int = size(A, 1)
  t::Int = size(A, 2)
  u::Int = size(A, 3)
  cartesianarray(Tuple{eltype(A)}, (m * s, n * t, l * u)) do i, j, k
    A[1 + rem(i - 1, s), 1 + rem(j - 1, t), 1 + rem(k - 1, u)]
  end
end

@inline function repmat(A::DenseMatrix, m::Int, n::Int, l::Int)
  s::Int = size(A, 1)
  t::Int = size(A, 2)
  cartesianarray(Tuple{eltype(A)}, (m * s, n * t, l)) do i, j, k
    A[1 + rem(i - 1, s), 1 + rem(j -1, t)]
  end
end

@inline function repmat(A::DenseVector, m::Int, n::Int, l::Int)
  s::Int = size(A, 1)
  cartesianarray(Tuple{eltype(A)}, (m * s, n, l)) do i, j, k
    A[1 + rem(i - 1, s)]
  end
end

@inline function repmat(A::DenseVector, m::Int, n::Int)
  s::Int = size(A, 1)
  cartesianarray(Tuple{eltype(A)}, (m * s, n)) do i, j
    A[1 + rem(i - 1, s)]
  end
end

@inline function repmat(A::DenseMatrix, m::Int, n::Int)
  s::Int, t::Int = size(A)
  cartesianarray(Tuple{eltype(A)}, (m * s, n * t)) do i, j
    A[1 + rem(i - 1, s), 1 + rem(j - 1, t)]
  end
end

@inline function repmat(A::DenseVector, m::Int)
  repmat(A, m, 1)
end

@inline function repmat(A::DenseMatrix, m::Int)
  repmat(A, m, 1)
end

baremodule NoInline

using Base
import Base: call, (==)

@noinline function rand(args...)
  Base.Random.rand(args...)
end

@noinline function randn(args...)
  Base.Random.rand(args...)
end

@noinline function rand!(args...)
  Base.Random.rand(args...)
end

@noinline function randn!(args...)
  Base.Random.rand(args...)
end

end
import .NoInline

@inline function rand(dims::Int...)
  arr = Array(Float64, dims...)
  map!(x -> NoInline.rand(Float64)::Float64, arr)
  return arr
end

@inline function randn(dims::Int...)
  arr = Array(Float64, dims...)
  map!(x -> NoInline.randn()::Float64, arr)
  return arr
end

@inline function rand(t::Tuple)
  arr = Array(Float64, t)
  map!(x -> NoInline.rand(Float64)::Float64, arr)
  return arr
end

@inline function randn(t::Tuple)
  arr = Array(Float64, t)
  map!(x -> NoInline.randn()::Float64, arr)
  return arr
end

@inline function rand(T::Type{Float32}, d::Int, dims::Int...)
  arr = Array(Float32, d, dims...)
  map!(x -> NoInline.rand(Float32)::Float32, arr)
  return arr
end

@inline function rand(T::Type, d::Int, dims::Int...)
  arr = Array(T, d, dims...)
  map!(x -> NoInline.rand(T)::T, arr)
  return arr
end

#@inline function randn(T::Type, d::Int, dims::Int...)
#  arr = Array(T, d, dims...)
#  map!(x -> NoInline.randn(T)::T, arr)
#  return arr
#end

@inline function rand(r::AbstractRNG, T::Type, d::Int, dims::Int...)
  arr = Array(T, d, dims...)
  map!(x -> NoInline.rand(r, T)::T, arr)
  return arr
end

#@inline function randn(r::AbstractRNG, T::Type, d::Int, dims::Int...)
#  arr = Array(T, d, dims...)
#  map!(x -> NoInline.randn(r, T)::T, arr)
#  return arr
#end

@inline function rand!{T}(arr::DenseArray{T})
  map!(x -> NoInline.rand(T)::T, arr)
end

#@inline function randn!{T}(arr::DenseArray{T})
#  map!(x -> NoInline.randn(T)::T, arr)
#end

@inline function rand!{T}(r::AbstractRNG, arr::DenseArray{T})
  map!(x -> NoInline.rand(r, T)::T, arr)
end

#@inline function randn!{T}(r::AbstractRNG, arr::DenseArray{T})
#  map!(x -> NoInline.randn(r, T)::T, arr)
#end

@inline function rand(args...)
  Base.Random.rand(args...)
end

@inline function randn(args...)
  Base.Random.rand(args...)
end

@inline function rand!(args...)
  Base.Random.rand(args...)
end

@inline function randn!(args...)
  Base.Random.rand(args...)
end

export indmin, indmax, sumabs2
export diag, diagm, trace, scale, eye, repmat, rand, randn, rand!, randn!

end


text:
__precompile__()

module SGSDist

    using Logging
    import DelimitedFiles: readdlm
    import Random: rand, randn, AbstractRNG, GLOBAL_RNG

    import DifferentialEquations: SDEProblem, solve, EulerHeun, EM, RKMil
    import Distributions: @distr_support, @check_args, ContinuousUnivariateDistribution, Normal, params, pdf, cdf, ccdf, logpdf, logcdf, logccdf, fit, mean, var, skewness, kurtosis
    import Interpolations: LinearInterpolation, Line
    import QuadGK: quadgk
    import SpecialFunctions: loggamma
    import StatsBase: var, skewness, kurtosis

    export SGS, params, mean, var, std, skewness, kurtosis, pdf, pdf_xmax, logN, logpdf, cdf, logcdf, ccdf, logccdf, fit, CAM1D, Hasselmann1D, RK4, Euler, fdt, rand, air, epsilon_condition, epsilon_constraint, ks_constraint, small_E_constraint, fit_mm, param_E, param_b, param_g, invitp

    include("constraints.jl")
    include("continuous.jl")
    include("fit.jl")
    include("markov.jl")
    include("sampledata.jl")

    # precompile hints
    # d = SGS()
    # rand(d, 1)

end


text:
# Operators
# ------------------------------------------------------------------------------

#=
    Number
=#

# Number--ParameterRef
Base.:(+)(lhs::C, rhs::ParameterRef) where C<:Number = PAE{C}(GAEv{C}(convert(Float64, lhs)), GAEp{C}(zero(C), rhs => +one(C)))
Base.:(-)(lhs::C, rhs::ParameterRef) where C<:Number = PAE{C}(GAEv{C}(convert(Float64, lhs)), GAEp{C}(zero(C), rhs => -one(C)))
Base.:(*)(lhs::C, rhs::ParameterRef) where C<:Number = PAE{C}(GAEv{C}(zero(C)), GAEp{C}(zero(C), rhs => lhs))

# Number--PAE
Base.:(+)(lhs::Number, rhs::PAE{C}) where C<:Number = PAE{C}(lhs+rhs.v, copy(rhs.p))
Base.:(-)(lhs::Number, rhs::PAE{C}) where C<:Number = PAE{C}(lhs-rhs.v,-rhs.p)
Base.:(*)(lhs::Number, rhs::PAE{C}) where C<:Number = PAE{C}(lhs*rhs.v, lhs*rhs.p)

#=
    ParameterRef
=#

# AbstractJuMPScalar
Base.:(-)(lhs::ParameterRef) = PAE{Float64}(GAEv{Float64}(0.0), GAEp{Float64}(0.0, lhs => -1.0))

# ParameterRef--Number
Base.:(+)(lhs::ParameterRef, rhs::Number) = (+)(rhs, lhs)
Base.:(-)(lhs::ParameterRef, rhs::Number) = (+)(-rhs, lhs)
Base.:(*)(lhs::ParameterRef, rhs::Number) = (*)(rhs, lhs)
Base.:(/)(lhs::ParameterRef, rhs::Number) = (*)(1.0 / rhs, lhs)

# ParameterRef--VariableRef
Base.:(+)(lhs::ParameterRef, rhs::JuMP.VariableRef) = PAE{Float64}(GAEv{Float64}(0.0, rhs => +1.0), GAEp{Float64}(0.0, lhs => 1.0))
Base.:(-)(lhs::ParameterRef, rhs::JuMP.VariableRef) = PAE{Float64}(GAEv{Float64}(0.0, rhs => -1.0), GAEp{Float64}(0.0, lhs => 1.0))

# ParameterRef--ParameterRef
Base.:(+)(lhs::ParameterRef, rhs::ParameterRef) = PAE{Float64}(GAEv{Float64}(0.0), GAEp{Float64}(0.0, lhs => 1.0, rhs => +1.0))
Base.:(-)(lhs::ParameterRef, rhs::ParameterRef) = PAE{Float64}(GAEv{Float64}(0.0), GAEp{Float64}(0.0, lhs => 1.0, rhs => -1.0))

# ParameterRef--GAEp
# this one is used internally only, becaus no other gunction returns a GAEp
Base.:(+)(lhs::ParameterRef, rhs::GAEp{C}) where C = (+)(GAEp{C}(zero(C), lhs => 1.0),  rhs)
Base.:(-)(lhs::ParameterRef, rhs::GAEp{C}) where C = (+)(GAEp{C}(zero(C), lhs => 1.0), -rhs)

# ParameterRef--GAEv/GenericAffExpr{C,VariableRef}
Base.:(+)(lhs::ParameterRef, rhs::GAEv{C}) where {C} = PAE{C}(copy(rhs),GAEp{C}(zero(C), lhs => 1.))
Base.:(-)(lhs::ParameterRef, rhs::GAEv{C}) where {C} = PAE{C}(-rhs,GAEp{C}(zero(C), lhs => 1.))

# ParameterRef--PAE{C}
Base.:(+)(lhs::ParameterRef, rhs::PAE{C}) where {C} = PAE{C}(copy(rhs.v),lhs+rhs.p)
Base.:(-)(lhs::ParameterRef, rhs::PAE{C}) where {C} = PAE{C}(-rhs.v,lhs-rhs.p)

#=
    VariableRef
=#

# VariableRef--ParameterRef
Base.:(+)(lhs::JuMP.VariableRef, rhs::ParameterRef) = PAE{Float64}(GAEv{Float64}(zero(Float64), lhs => 1.0),GAEp{Float64}(zero(Float64), rhs =>  1.0))
Base.:(-)(lhs::JuMP.VariableRef, rhs::ParameterRef) = PAE{Float64}(GAEv{Float64}(zero(Float64), lhs => 1.0),GAEp{Float64}(zero(Float64), rhs => -1.0))

# VariableRef--GenericAffExpr{C,ParameterRef}
Base.:(+)(lhs::JuMP.VariableRef, rhs::GAEp{C}) where {C} = PAE{C}(GAEv{C}(zero(C), lhs => 1.),copy(rhs))
Base.:(-)(lhs::JuMP.VariableRef, rhs::GAEp{C}) where {C} = PAE{C}(GAEv{C}(zero(C), lhs => 1.),-rhs)

# VariableRef--PAE{C}
Base.:(+)(lhs::JuMP.VariableRef, rhs::PAE{C}) where {C} = PAE{C}(lhs + rhs.v, copy(rhs.p))
Base.:(-)(lhs::JuMP.VariableRef, rhs::PAE{C}) where {C} = PAE{C}(lhs - rhs.v, -rhs.p)

#=
    GenericAffExpr{C,VariableRef}
=#

# GenericAffExpr{C,VariableRef}--ParameterRef
Base.:(+)(lhs::GAEv{C}, rhs::ParameterRef) where {C} = (+)(rhs,lhs)
Base.:(-)(lhs::GAEv{C}, rhs::ParameterRef) where {C} = (+)(-rhs,lhs)

# GenericAffExpr{C,VariableRef}--GenericAffExpr{C,ParameterRef}
Base.:(+)(lhs::GAEv{C}, rhs::GAEp{C}) where {C} = PAE{C}(copy(lhs),copy(rhs))
Base.:(-)(lhs::GAEv{C}, rhs::GAEp{C}) where {C} = PAE{C}(copy(lhs),-rhs)

# GenericAffExpr{C,VariableRef}--PAE{C}
Base.:(+)(lhs::GAEv{C}, rhs::PAE{C}) where {C} = PAE{C}(lhs+rhs.v,copy(rhs.p))
Base.:(-)(lhs::GAEv{C}, rhs::PAE{C}) where {C} = PAE{C}(lhs-rhs.v,-rhs.p)

#=
    GenericAffExpr{C,ParameterRef}/GAEp
=#

# GenericAffExpr{C,ParameterRef}--ParameterRef
# DONE in JuMP

# GenericAffExpr{C,ParameterRef}--VariableRef
Base.:(+)(lhs::GAEp{C}, rhs::JuMP.VariableRef) where {C} = (+)(rhs,lhs)
Base.:(-)(lhs::GAEp{C}, rhs::JuMP.VariableRef) where {C} = (+)(-rhs,lhs)

# GenericAffExpr{C,ParameterRef}--GenericAffExpr{C,VariableRef}
Base.:(+)(lhs::GAEp{C}, rhs::GAEv{C}) where {C} = (+)(rhs,lhs)
Base.:(-)(lhs::GAEp{C}, rhs::GAEv{C}) where {C} = (+)(-rhs,lhs)

# GenericAffExpr{C,ParameterRef}--GenericAffExpr{C,ParameterRef}
# DONE in JuMP

# GenericAffExpr{C,ParameterRef}--PAE{C}

#=
    PAE{C}
=#

Base.:(-)(lhs::PAE{C}) where C = PAE{C}(-lhs.v, -lhs.p)

# Number--PAE
Base.:(+)(lhs::PAE, rhs::Number) = (+)(rhs,lhs)
Base.:(-)(lhs::PAE, rhs::Number) = (+)(-rhs,lhs)
Base.:(*)(lhs::PAE, rhs::Number) = (*)(rhs,lhs)

# PAE{C}--ParameterRef
Base.:(+)(lhs::PAE{C}, rhs::ParameterRef) where {C} = (+)(rhs,lhs)
Base.:(-)(lhs::PAE{C}, rhs::ParameterRef) where {C} = (+)(-rhs,lhs)

# VariableRef--PAE{C}
Base.:(+)(lhs::PAE{C}, rhs::JuMP.VariableRef) where {C} = (+)(rhs,lhs)
Base.:(-)(lhs::PAE{C}, rhs::JuMP.VariableRef) where {C} = (+)(-rhs,lhs)

# PAE{C}--GenericAffExpr{C,VariableRef}
# PAE{C}--GenericAffExpr{C,ParameterRef}
Base.:(+)(lhs::PAE{C}, rhs::GAE{C,V}) where {C,V} = (+)(rhs,lhs)
Base.:(-)(lhs::PAE{C}, rhs::GAE{C,V}) where {C,V} = (+)(-rhs,lhs)

# PAE{C}--PAE{C}
Base.:(+)(lhs::PAE{C}, rhs::PAE{C}) where {C} = PAE{C}(lhs.v+rhs.v,lhs.p+rhs.p)
Base.:(-)(lhs::PAE{C}, rhs::PAE{C}) where {C} = PAE{C}(lhs.v-rhs.v,lhs.p-rhs.p)


text:
function sumby_sort(id6,v1)
  a = DataFrames.DataFrame(id6=id6,v1=v1)
  sort!(a, cols=:id6; alg=RadixSort)
  a
end

@time sumby_sort(id6,v1) #

@time sumby(id6,v1)


text:
using SchroedingersSmoke, CLArrays
using Colors, GPUArrays, GeometryTypes, GLAbstraction

ArrayType = CLArray

vol_size = (4,2,2)# box size
dims = (64, 32, 32) .* 2 # volume resolution
hbar = 0.1f0      # Planck constant
dt = 1f0/48f0     # time step

jet_velocity = (1f0, 0f0, 0f0)
nozzle_cen = Float32.((2-1.7, 1-0.034, 1+0.066))
nozzle_len = 0.5f0
nozzle_rad = 0.5f0
n_particles = 500   # number of particles

isf = ISF{ArrayType, UInt32, Float32}(vol_size, dims, hbar, dt);

# initialize psi
psi = ArrayType([(one(Complex64), one(Complex64) * 0.01f0) for i=1:dims[1], j=1:dims[2], k=1:dims[3]]);

psi .= normalize_psi.(psi);

kvec = jet_velocity ./ hbar;
omega = sum(jet_velocity.^2f0) / (2f0*hbar);
isjetarr = isjet.(isf.positions, (nozzle_cen,), (nozzle_len,), (nozzle_rad,))
# constrain velocity
for iter = 1:10
    restrict_velocity!(isf, psi, kvec, isjetarr, 0f0)
    pressure_project!(isf, psi)
end

particles = ArrayType(map(x-> (0f0, 0f0, 0f0), 1:(10^6)))

add_particles!(particles, 1:n_particles, nozzle_cen, nozzle_rad)

using GLVisualize;
w = glscreen(color = RGBA(0f0, 0f0, 0f0, 0f0), resolution = (1920, 1080));
lincolor = RGBA{Float32}(0.0f0,0.74736935f0,1.0f0,0.1f0)
offset = translationmatrix(Vec3f0(0))
dircolor = map(x-> RGBA{Float32}(x[1], x[2], x[3], 0.1), isf.velocity);
dircolorcpu = Array{RGBA{Float32}}(size(isf.velocity))
copy!(dircolorcpu, dircolor)
dircolorviz = visualize(
    dircolorcpu, :absorption, model = offset,
    dimensions = Vec3f0(vol_size)
).children[]
# first view sets up camera
_view(dircolorviz, camera = :perspective, position = Vec3f0(4, -7, 2.1), lookat = Vec3f0(4, 0, 2))

_view(visualize(
    AABB(Vec3f0(0), Vec3f0(vol_size)), :lines,
    color =	lincolor,
    model = offset,
), camera = :perspective)
_view(visualize(
    "Velocity as Colors",
    start_position = Point3f0(0), billboard = false,
    model = translationmatrix(Vec3f0(0.1, -0.1, -0.2)) * rotationmatrix_x(Float32(0.5pi)),
    relative_scale = 0.2f0, color = RGBA(1f0,1f0,1f0,1f0)
), camera = :perspective)
offset = translationmatrix(Vec3f0(0, 0, 2))

particlescpu = Array(particles)
particle_vis = visualize(
    (GeometryTypes.Circle(GeometryTypes.Point2f0(0), 0.006f0), reinterpret(Point3f0, particlescpu)),
    boundingbox = nothing, # don't waste time on bb computation
    model = offset,
    color = fill(RGBA{Float32}(0, 0, 0, 0.09), length(particles)),
    billboard = true
).children[]
_view(particle_vis, camera = :perspective)
_view(visualize(
    AABB(Vec3f0(0), Vec3f0(vol_size)), :lines,
    color =	lincolor,
    model = offset,
), camera = :perspective)

particle_vis[:color][1:n_particles] = map(1:n_particles) do i
    xx = (i / n_particles) * 2pi
    RGBA{Float32}((sin(xx) + 1) / 2, (cos(xx) + 1.0) / 2.0, 0.0, 0.1)
end
_view(visualize(
    "Smoke Particle",
    start_position = Point3f0(0), billboard = false,
    model = translationmatrix(Vec3f0(0.1, -0.1, 4.1)) * rotationmatrix_x(Float32(0.5pi)),
    relative_scale = 0.2f0, color = RGBA(1f0,1f0,1f0,1f0)
), camera = :perspective)

offset = translationmatrix(Vec3f0(4, 0, 0))
magnitude = map(x->sqrt(sum(x .* x)), isf.velocity);
magnitudecpu = Array(magnitude)
magnitudeviz = visualize(
    magnitudecpu, :mip,
    model = offset,
    color_map = [RGBA(0f0, 0f0, 0f0, 0f0), RGBA(0.2f0, 0f0, 0.9f0, 0.9f0), RGBA(0.9f0, 0.2f0, 0f0, 0f0)],
    dimensions = Vec3f0(vol_size)
).children[]
_view(magnitudeviz, camera = :perspective)
_view(visualize(
    AABB(Vec3f0(0), Vec3f0(vol_size)), :lines,
    model = offset,
    color =	lincolor,
), camera = :perspective)
_view(visualize(
    "Velocity Magnitude",
    start_position = Point3f0(0), billboard = false,
    model = translationmatrix(Vec3f0(4.1, -0.1, -0.2)) * rotationmatrix_x(Float32(0.5pi)),
    relative_scale = 0.2f0, color = RGBA(1f0,1f0,1f0,1f0)
), camera = :perspective)

offset = translationmatrix(Vec3f0(4, 0, 2))
directions = copy(isf.velocity)
directionscpu = Array(directions)
directionsviz = visualize(
    reinterpret(Vec3f0, directionscpu),
    model = offset,
    color_map = [RGBA(0f0, 0f0, 0f0, 0f0), RGBA(0.2f0, 0f0, 0.9f0, 0.9f0), RGBA(0.9f0, 0.2f0, 0f0, 0f0)],
    ranges = map(x-> 0:x, vol_size)
).children[]
_view(directionsviz, camera = :perspective)
_view(visualize(
    "Velocity as Vectorfield",
    start_position = Point3f0(0), billboard = false,
    model = translationmatrix(Vec3f0(4.1, -0.1, 4.1)) * rotationmatrix_x(Float32(0.5pi)),
    relative_scale = 0.2f0, color = RGBA(1f0,1f0,1f0,1f0)
), camera = :perspective)

_view(visualize(
    AABB(Vec3f0(0), Vec3f0(vol_size)), :lines,
    model = offset,
    color =	lincolor,
), camera = :perspective)

dt = isf.dt; d = isf.d
iter = 1

io, buffer = GLVisualize.create_video_stream("test.mkv", w)

# main simulation loop
while isopen(w)
    t = iter * dt
    # incompressible Schroedinger flow
    schroedinger_flow!(isf, psi)
    psi .= normalize_psi.(psi)
    pressure_project!(isf, psi)

    start = mod((iter - 1) * n_particles + 1, length(particles))
    stop = start + n_particles - 1
    add_particles!(particles, start:stop, nozzle_cen, nozzle_rad)
    particle_vis[:color][start:stop] = map(1:n_particles) do i
        xx = (i / n_particles) * 2pi
        RGBA{Float32}((sin(xx) + 1) / 2, (cos(xx) + 1.0) / 2.0, mod(iter, 100) / 100, 0.1)
    end
    #constrain velocity
    restrict_velocity!(isf, psi, kvec, isjetarr, omega*t)
    pressure_project!(isf, psi)
    velocity_one_form!(isf, psi, isf.hbar)
    # inplace StaggeredSharp
    dinv = inv.(isf.d)
    broadcast!((x, y)-> x .* y, isf.velocity, isf.velocity, (dinv,))
    staggered_advect!(particles, isf)

    copy!(particlescpu, particles)
    GLAbstraction.set_arg!(particle_vis, :position, reinterpret(Point3f0, particlescpu))

    dircolor .= (x-> RGBA{Float32}(x[1] * 10f0, x[2] * 10f0, x[3] * 10f0, 1.0)).(isf.velocity)
    copy!(dircolorcpu, dircolor)
    GLAbstraction.set_arg!(dircolorviz, :volumedata, dircolorcpu)
    magnitude .= (x->sqrt(sum(x .* x))).(isf.velocity)
    copy!(magnitudecpu, magnitude)
    GLAbstraction.set_arg!(magnitudeviz, :volumedata, magnitudecpu)

    copy!(directionscpu, isf.velocity)
    GLAbstraction.set_arg!(directionsviz, :rotation, vec(reinterpret(Vec3f0, directionscpu)))
    GLWindow.poll_glfw()
    GLWindow.reactive_run_till_now()
    GLWindow.render_frame(w)
    GLWindow.swapbuffers(w)
    GLVisualize.add_frame!(io, w, buffer)
    iter += 1
end
close(io)
GLWindow.destroy!(w)
pwd()


text: # Proximal newton and quasi newton methods


text:
function GLONASSNumericalPropagator(arg0::GLONASSNumericalPropagator_Builder)
    return GLONASSNumericalPropagator((GLONASSNumericalPropagator_Builder,), arg0)
end

function add_additional_equations(obj::AbstractIntegratedPropagator, arg0::AdditionalEquations)
    return jcall(obj, "addAdditionalEquations", void, (AdditionalEquations,), arg0)
end

function add_additional_state_provider(obj::AbstractPropagator, arg0::AdditionalStateProvider)
    return jcall(obj, "addAdditionalStateProvider", void, (AdditionalStateProvider,), arg0)
end

function add_event_detector(obj::AbstractIntegratedPropagator, arg0::EventDetector)
    return jcall(obj, "addEventDetector", void, (EventDetector,), arg0)
end

function clear_events_detectors(obj::AbstractIntegratedPropagator)
    return jcall(obj, "clearEventsDetectors", void, ())
end

function get_additional_state_providers(obj::AbstractPropagator)
    return jcall(obj, "getAdditionalStateProviders", List, ())
end

function get_attitude_provider(obj::AbstractPropagator)
    return jcall(obj, "getAttitudeProvider", AttitudeProvider, ())
end

function get_basic_dimension(obj::AbstractIntegratedPropagator)
    return jcall(obj, "getBasicDimension", jint, ())
end

function get_calls(obj::AbstractIntegratedPropagator)
    return jcall(obj, "getCalls", jint, ())
end

function get_events_detectors(obj::AbstractIntegratedPropagator)
    return jcall(obj, "getEventsDetectors", Collection, ())
end

function get_frame(obj::AbstractPropagator)
    return jcall(obj, "getFrame", Frame, ())
end

function get_generated_ephemeris(obj::AbstractIntegratedPropagator)
    return jcall(obj, "getGeneratedEphemeris", BoundedPropagator, ())
end

function get_glonass_orbital_elements(obj::GLONASSNumericalPropagator)
    return jcall(obj, "getGLONASSOrbitalElements", GLONASSOrbitalElements, ())
end

function get_initial_state(obj::AbstractPropagator)
    return jcall(obj, "getInitialState", SpacecraftState, ())
end

function get_managed_additional_states(obj::AbstractIntegratedPropagator)
    return jcall(obj, "getManagedAdditionalStates", Vector{JString}, ())
end

function get_mode(obj::AbstractPropagator)
    return jcall(obj, "getMode", jint, ())
end

function get_mu(obj::AbstractIntegratedPropagator)
    return jcall(obj, "getMu", jdouble, ())
end

function get_pv_coordinates(obj::AbstractPropagator, arg0::AbsoluteDate, arg1::Frame)
    return jcall(obj, "getPVCoordinates", TimeStampedPVCoordinates, (AbsoluteDate, Frame), arg0, arg1)
end

function get_pv_in_pz90(obj::GLONASSNumericalPropagator, arg0::SpacecraftState)
    return jcall(obj, "getPVInPZ90", PVCoordinates, (SpacecraftState,), arg0)
end

function is_additional_state_managed(obj::AbstractIntegratedPropagator, arg0::JString)
    return jcall(obj, "isAdditionalStateManaged", jboolean, (JString,), arg0)
end

function propagate(obj::AbstractIntegratedPropagator, arg0::AbsoluteDate, arg1::AbsoluteDate)
    return jcall(obj, "propagate", SpacecraftState, (AbsoluteDate, AbsoluteDate), arg0, arg1)
end

function propagate(obj::GLONASSNumericalPropagator, arg0::AbsoluteDate)
    return jcall(obj, "propagate", SpacecraftState, (AbsoluteDate,), arg0)
end

function reset_initial_state(obj::AbstractPropagator, arg0::SpacecraftState)
    return jcall(obj, "resetInitialState", void, (SpacecraftState,), arg0)
end

function set_attitude_provider(obj::AbstractIntegratedPropagator, arg0::AttitudeProvider)
    return jcall(obj, "setAttitudeProvider", void, (AttitudeProvider,), arg0)
end

function set_ephemeris_mode(obj::AbstractIntegratedPropagator)
    return jcall(obj, "setEphemerisMode", void, ())
end

function set_ephemeris_mode(obj::AbstractIntegratedPropagator, arg0::OrekitStepHandler)
    return jcall(obj, "setEphemerisMode", void, (OrekitStepHandler,), arg0)
end

function set_master_mode(obj::AbstractIntegratedPropagator, arg0::OrekitStepHandler)
    return jcall(obj, "setMasterMode", void, (OrekitStepHandler,), arg0)
end

function set_master_mode(obj::AbstractPropagator, arg0::jdouble, arg1::OrekitFixedStepHandler)
    return jcall(obj, "setMasterMode", void, (jdouble, OrekitFixedStepHandler), arg0, arg1)
end

function set_mu(obj::AbstractIntegratedPropagator, arg0::jdouble)
    return jcall(obj, "setMu", void, (jdouble,), arg0)
end

function set_reset_at_end(obj::AbstractIntegratedPropagator, arg0::jboolean)
    return jcall(obj, "setResetAtEnd", void, (jboolean,), arg0)
end

function set_slave_mode(obj::AbstractIntegratedPropagator)
    return jcall(obj, "setSlaveMode", void, ())
end


text:
import ..Datasets: dataset, datafile, reader, Mode, Test

sts_init() = register(DataDep(
    "GLUE-STS",
    """
     The Semantic Textual Similarity Benchmark (STS) task (GLUE version)
    """,
    "https://firebasestorage.googleapis.com/v0/b/mtl-sentence-representations.appspot.com/o/data%2FSTS-B.zip?alt=media&token=bddb94a7-8706-4e0d-a694-1109e12273b5",
    "e60a6393de5a8b5b9bac5020a1554b54e3691f9d600b775bd131e613ac179c85";
    post_fetch_method = fn -> begin
      mv(fn, "STS-B.zip")
      DataDeps.unpack("STS-B.zip")
      innerdir = "STS-B"
      innerfiles = readdir(innerdir)
      mv.(joinpath.(innerdir, innerfiles), innerfiles)
      rm(innerdir)
    end
))

struct STS <: Dataset end

function dataset(::Type{M}, d::STS) where M <: Mode
    ds = reader(datafile(M, d))
    header = split(take!(ds), '\t')
    field_num = length(header)
    needed_field = (8,9,10)
    rds = get_channels(String, length(needed_field))

    task = @async begin
        for s ∈ ds
            s = split(s, '\t')
            length(s) != field_num && continue
            for (i, j) ∈ enumerate(needed_field)
                put!(rds[i], s[j])
            end
        end
    end
    for rd ∈ rds
        bind(rd, task)
    end

    rds
end

function trainfile(::STS)
    datadep"GLUE-STS/train.tsv"
end

function devfile(::STS)
    datadep"GLUE-STS/dev.tsv"
end

function testfile(::STS)
    datadep"GLUE-STS/test.tsv"
end


text:
###
### Contract: mutation without explicit copy is undefined behavior.
###

const PRETTY_PRINT = Ref(true)
const CoeffType = Union{Rational{<:Integer},Integer}

###
### AbstractMonomial: isless, degree
###
abstract type AbstractMonomial <: Number end

Base.isone(x::AbstractMonomial) = iszero(degree(x))
Base.one(::Type{<:T}) where {T<:AbstractMonomial} = T()

Base.:-(y::AbstractMonomial) = Term(-1, y)
Base.:+(y::T, x::T) where {T <: AbstractMonomial} = Term(y) + Term(x)
Base.:*(c::CoeffType, m::AbstractMonomial) = Term(c, m)
Base.:*(m::AbstractMonomial, c::CoeffType) = c * m
Base.:^(y::T, n::Integer) where {T <: AbstractMonomial} = iszero(n) ? T() : Base.power_by_squaring(y, n)
monomialtype(x::AbstractMonomial) = typeof(x)

###
### AbstractTerm: `coeff` and `monomial`
###
abstract type AbstractTerm <: Number end
Base.promote_rule(::Type{C}, ::Type{M}) where {C<:CoeffType,M<:AbstractMonomial} = Term{C,M}
Base.promote_rule(t::Type{<:AbstractTerm}, ::Type{<:AbstractMonomial}) = t
Base.promote_rule(t::Type{<:AbstractTerm}, ::Type{<:CoeffType}) = t
emptyterm(T::Type{<:AbstractTerm}) = T[]

degree(x::AbstractTerm) = degree(monomial(x))
ismatch(x::T, y::T) where {T<:AbstractTerm} = monomial(x) == monomial(y)
Base.:(==)(x::AbstractTerm, y::AbstractTerm) = x === y || (coeff(x) == coeff(y) && monomial(x) == monomial(y))
Base.copy(x::T) where {T<:AbstractTerm} = T(copy(coeff(x)), copy(monomial(x)))
Base.isless(x::T, y::T) where {T<:AbstractTerm} = isless(monomial(x), monomial(y))
Base.iszero(x::AbstractTerm) = iszero(coeff(x))
Base.isone(x::AbstractTerm) = isone(coeff(x)) && isone(monomial(x))
Base.zero(t::T) where {T<:AbstractTerm} = parameterless_type(T)(zero(coeff(t)), one(monomial(t)))
Base.one(t::T) where {T<:AbstractTerm} = parameterless_type(T)(one(coeff(t)), one(monomial(t)))
Base.isinteger(x::AbstractTerm) = isinteger(coeff(x))

monomialtype(x::AbstractTerm) = monomialtype(monomial(x))

Base.:*(x::T, y::T) where {T<:AbstractTerm} = T(coeff(x) * coeff(y), monomial(x) * monomial(y))
# TODO
function Base.:(/)(y::T, x::T) where {T<:AbstractTerm}
    m, fail = monomial(y) / monomial(x)
    if fail
        T(coeff(y), m), fail
    else
        T(coeff(y)/coeff(x), m), fail
    end
end

function Base.:+(x::T, y::T) where {T<:AbstractTerm}
    if ismatch(x, y)
        c = coeff(x) + coeff(y)
        return iszero(c) ? MPoly(emptyterm(T)) : MPoly(T(c, monomial(x)))
    else
        if iszero(x) && iszero(y)
            return MPoly(emptyterm(T))
        elseif iszero(x)
            return MPoly(T[y])
        elseif iszero(y)
            return MPoly(T[x])
        else
            if x < y
                x, y = y, x
            end
            return MPoly(T[x, y])
        end
    end
end

Base.:-(x::T) where {T<:AbstractTerm} = T(-coeff(x), monomial(x))
function Base.:-(x::T, y::T) where {T<:AbstractTerm}
    if ismatch(x, y)
        c = coeff(x) - coeff(y)
        return iszero(c) ? MPoly(T[]) : MPoly(T(c, monomial(x)))
    else
        y = -y
        if x < y
            x, y = y, x
        end
        return MPoly(T[x, y])
    end
end

function Base.gcd(x::T, y::T) where {T<:AbstractTerm}
    #g, a, b = gcd(monomial(x), monomial(y))
    g = gcd(monomial(x), monomial(y))
    gr = gcd(x.coeff, y.coeff)
    return T(gr, g)#, Term(x.coeff / gr, a), Term(y.coeff / gr, b)
end

print_coeff(io::IO, coeff) = isinteger(coeff) ? print(io, Integer(coeff)) : print(io, coeff)
function Base.show(io::IO, x::AbstractTerm)
    printed = false
    if coeff(x) != 1
        if coeff(x) == -1
            print(io, '-')
        else
            print_coeff(io, coeff(x))
            printed = true
            PRETTY_PRINT[] || print(io, '*')
        end
    end
    m = monomial(x)
    if !(printed && isone(m))
        show(io, m)
    end
end

# default term type
struct Term{C,M<:AbstractMonomial} <: AbstractTerm
    coeff::C
    monomial::M
end
coeff(x::Term) = x.coeff
monomial(x::Term) = x.monomial
Term{M}(x) where {M<:AbstractMonomial} = Term(x, M())
Term(x::M) where {M<:AbstractMonomial} = Term(1, x)
Term{A,B}(x::M) where {A,B,M<:AbstractMonomial} = Term(1, x)
Term{A,B}(x) where {A,B} = Term(x, B())
#const EMPTY_TERM = Term[]

###
### AbstractPolynomial: terms, copy
###
abstract type AbstractPolynomial <: Number end

Base.promote_rule(p::Type{<:AbstractPolynomial}, ::Type{<:AbstractMonomial}) = p
Base.promote_rule(p::Type{<:AbstractPolynomial}, ::Type{<:AbstractTerm}) = p
Base.promote_rule(p::Type{<:AbstractPolynomial}, ::Type{<:CoeffType}) = p

Base.iszero(x::AbstractPolynomial) = isempty(terms(x))
Base.isone(x::AbstractPolynomial) = (ts = terms(x); length(ts) == 1 && isone(only(ts)))
Base.:(==)(x::T, y::T) where {T<:AbstractPolynomial} = x === y || (terms(x) == terms(y))
monomialtype(p::AbstractPolynomial) = monomialtype(lt(p))

function Base.show(io::IO, p::AbstractPolynomial)
    ts = terms(p)
    if isempty(ts)
        print(io, 0)
        return
    end
    n = length(ts)
    t1, tr = Iterators.peel(ts)
    show(io, t1)
    for t in tr
        if t.coeff < 0
            print(io, " - ")
            t = -t
        else
            print(io, " + ")
        end
        show(io, t)
    end
end


text:
# FileManager
FileManager() = FileManager(create_file_manager())

"""
    create_file_manager() -> CXFileManager
Return a pointer to a `clang::FileManager` object.

For now, `FileSystemOptions` is set to nothing and `llvm::vfs::FileSystem` defaults to the
"real" file system, as seen by the operating system.

TODO: support custom `FileSystemOptions` and `llvm::vfs::FileSystem`
"""
function create_file_manager()
    status = Ref{CXInit_Error}(CXInit_NoError)
    mgr = clang_FileManager_create(status)
    @assert status[] == CXInit_NoError
    return mgr
end

dispose(x::FileManager) = clang_FileManager_dispose(x)

function PrintStats(mgr::FileManager)
    @check_ptrs mgr
    return clang_FileManager_PrintStats(mgr)
end

"""
    getFileEntry(filemgr::FileManager, filename::AbstractString; open_file::Bool=false, cache_failure::Bool=true) -> FileEntry
Get a file entry from the file manager.

If `open_file` is true, the file will be opened.
If `cache_failure` is true, the failure that this file does not exist will be cached.
"""
function getFileEntry(filemgr::FileManager, filename::AbstractString; open_file::Bool=false,
                      cache_failure::Bool=true)
    @check_ptrs filemgr
    GC.@preserve filename begin
        ref = clang_FileManager_getFileRef(filemgr, filename, open_file, cache_failure)
        @assert ref != C_NULL "failed to create a FileRef to $filename."
        entry = clang_FileEntryRef_getFileEntry(ref)
        clang_FileEntryRef_dispose(ref)
    end
    return FileEntry(entry)
end


text:
# This file was generated, do not modify it. # hide
Random.seed!(523)
perm = randperm(length(y))
X = X[perm,:]
y = y[perm];


text:
using ModelingToolkit, OrdinaryDiffEq
using ModelingToolkitStandardLibrary.Blocks
using ModelingToolkitStandardLibrary.Blocks: Dₜ, t

#=
Testing strategy:
The general strategy is to test systems using simple intputs where the solution is known on closed form. For algebraic systems (without differential variables), an integrator with a constant input is often used together with the system under test.
=#

@testset "Gain" begin
    @info "Testing Gain"
    @named c = Gain(42)
    @named int = Integrator(; k=1)
    @named iosys = ODESystem([int.u~c.y, c.u~1], t, systems=[c, int])
    sys = structural_simplify(iosys)
    prob = ODEProblem(sys, Pair[], (0.0, 1.0))
    sol = solve(prob, Rosenbrock23(), saveat=0:0.1:1)
    @test sol[int.y] ≈ 42 .* (0:0.1:1)

    # Matrix gain
    @named c = Gain([2 0; 0 3])
    ints = [Integrator(; k=1, name=Symbol("int$i")) for i in 1:2]
    @named iosys = ODESystem([
        ints[1].u~c.y[1],
        ints[2].u~c.y[2],
        c.u[1]~1,
        c.u[2]~2,
        ], t, systems=[c; ints])
    sys = structural_simplify(iosys)
    prob = ODEProblem(sys, Pair[], (0.0, 1.0))
    sol = solve(prob, Rosenbrock23(), saveat=0:0.1:1)
    @test sol[ints[1].y] ≈ 2 .* (0:0.1:1) # 2 * 1
    @test sol[ints[2].y] ≈ 6 .* (0:0.1:1) # 3 * 2
end

@testset "Sum" begin
    @info "Testing Sum"
    @named s = Sum(2)
    ints = [Integrator(; k=1, name=Symbol("int$i")) for i in 1:2]
    @named iosys = ODESystem([
        ints[1].u~1,
        ints[2].u~2,
        ints[1].y~s.u[1],
        ints[2].y~s.u[2],
        ], t, systems=[s; ints])
    sys = structural_simplify(iosys)
    prob = ODEProblem(sys, Pair[], (0.0, 1.0))
    sol = solve(prob, Rosenbrock23(), saveat=0:0.1:1)
    @test sol[s.y] ≈ 3 .* (0:0.1:1)

    @named s = Sum([1, -2])
    ints = [Integrator(; k=1, name=Symbol("int$i")) for i in 1:2]
    @named iosys = ODESystem([
        ints[1].u~1,
        ints[2].u~1,
        ints[1].y~s.u[1],
        ints[2].y~s.u[2],
        ], t, systems=[s; ints])
    sys = structural_simplify(iosys)
    prob = ODEProblem(sys, Pair[], (0.0, 1.0))
    sol = solve(prob, Rosenbrock23(), saveat=0:0.1:1)
    @test sol[s.y] ≈ (1 + (-2)) .* (0:0.1:1)
end


text:
using StaticArrays, BenchmarkTools, ForwardDiff
BenchmarkTools.DEFAULT_PARAMETERS.seconds = 10
show_ind_bench = true
show_jac = false
N = 1000000
println("Benchmarking evolution of Full Discrete systems (with Jac and eom and state)")
println("The eom will always return an SVector, since this is the fastest for")
println("the Jacobian. Evolution for total of $N steps")
u = rand(3); un = rand(3); J = rand(3,3)
su = SVector{3}(u); sun = SVector{3}(un);
mu = MVector{3}(u); mun = MVector{3}(un)
@inline function eom_towel(x)
  x1, x2, x3 = x[1], x[2], x[3]
  SVector(3.8*x1*(1-x1)-0.05*(x2+0.35)*(1-2*x3),
  0.1*((x2+0.35)*(1-2*x3)-1 )*(1-1.9*x1),
  3.78*x3*(1-x3)+0.2*x2)
end

println("V1: immutable struct with state SVector (new system at every evolve)")
struct dds1{S<:Real, F, J} #<: DynamicalSystem
  state::SVector{3,S}
  eom::F
  jacob!::J
end
function dds1(u0, eom)
  @inline jac!(J, x) = ForwardDiff.jacobian!(J, eom, x)
  return dds1(u0, eom, jac!)
end
@inline function evolve1(s::dds1, N::Int)
  for i in 1:N
    s = (dds1(s.eom(s.state), s.eom))
  end
  return s
end
s1 = dds1(su, eom_towel)
s1 = evolve1(s1, 1000)
b1 = @benchmark (s = evolve1($s1, $N))
j1 = @benchmark ()
display(b1)

println("V2: immutable struct with state Vector (changed in-place)")
struct dds2{S<:Real, F, J} #<: DynamicalSystem
  state::Vector{S}
  eom::F
  jacob!::J
end
function dds2(u0, eom)
  @inline jac!(J, x) = ForwardDiff.jacobian!(J, eom, x)
  return dds2(u0, eom, jac!)
end
@inline function evolve2!(s::dds2, N::Int)
  eom = s.eom
  for i in 1:N
    s.state .= eom(s.state)
  end
end
s2 = dds2(u, eom_towel)
evolve2!(s2, 1000)
b2 = @benchmark evolve2!($s2, $N)
display(b2)

println("V3: immutable struct with state MVector (changed in-place)")
struct dds3{S<:Real, F, J} #<: DynamicalSystem
  state::MVector{3, S}
  eom::F
  jacob!::J
end
function dds3(u0, eom)
  @inline jac!(J, x) = ForwardDiff.jacobian!(J, eom, x)
  return dds3(u0, eom, jac!)
end
@inline function evolve3!(s::dds3, N::Int)
  eom = s.eom
  for i in 1:N
    s.state .= eom(s.state)
  end
end
s3 = dds3(mu, eom_towel)
evolve3!(s3, 1000)
b3 = @benchmark evolve3!($s3, $N)
display(b3)

beches = [minimum(b) for b in [b1,b2,b3]]
i = indmin(beches)
println("Conclusions:")
println("Minimum time of evolve-call for N=$(N): v$i")
println("Judged vs v$i:")
for (j,b) in enumerate(beches)
  j == i && continue
  println("v$i vs v$j")
  sleep(0.01)
  display(judge(beches[i],beches[j]))
  sleep(0.01)
end

#### This is broken in Base (julia;s fault.) The immutable should be faster normally.
using BenchmarkTools
println("One dimensional benchmarks")
struct ds1d1{S<:Real, F}
  state::S
  eom::F
end
function evolve1d1(s::ds1d1, N::Int)
  x = deepcopy(s.state)
  f = s.eom
  for i in 1:N
    x = f(x)
  end
  return ds1d1(x, s.eom)
end

mutable struct ds2d2{S<:Real, F}
  state::S
  eom::F
end
function evolve2d2(s::ds2d2, N::Int)
  x = deepcopy(s.state)
  f = s.eom
  for i in 1:N
    x = f(x)
  end
  s.state = x
end

function bench(N)
  x0 = rand()

  @inline eom_logistic(x) = 4*x*(1-x)
  d1 = ds1d1(x0, eom_logistic)
  println("V1: Immutable struct")
  b1 = @benchmark evolve1d1($d1, $N)
  display(b1)

  println("V2: mutable struct")
  d2 = ds2d2(x0, eom_logistic)
  b2 = @benchmark evolve2d2($d2, $N)
  display(b2)
end

println("For N = 1000")
bench(1000)
println("For N = 1000000")
bench(1000000)

## Profiling:
Profile.clear()
N = 1000000000
x0 = rand()
@inline eom_logistic(x) = 4*x*(1-x)
d1 = ds1d1(x0, eom_logistic)
evolve1d1(d1, N)
@profile evolve1d1(d1, N)
Profile.print()

Profile.clear()
N = 1000000000
x0 = rand()
@inline eom_logistic(x) = 4*x*(1-x)
d2 = ds2d2(x0, eom_logistic)
evolve2d2(d2, N)
@profile evolve2d2(d2, N)
Profile.print()

## Profiling Juno
# Profile.clear()
# N = 10000000
# x0 = rand()
# eom_logistic(x) = 4*x*(1-x)
# d1 = ds1d1(x0, eom_logistic)
# @profile evolve1d1(d1, N)
# Atom.Profiler.tree()
#
# N = 10000000
# x0 = rand()
# eom_logistic(x) = 4*x*(1-x)
# d1 = ds1d1(x0, eom_logistic)
# @profile evolve1d1(d1, N)
# Juno.profiletree()


text:
include("docs/make.jl")
mv("docs\\build\\assets\\siteinfo.js","docs\\build\\siteinfo.js")
cp("docs\\build","..\\Repos-Pages\\v0.1.0",force=true)


text:
module Joice

using DataFrames, Dates

include("JoiceTypeDefinitions.jl")
include("JoiceTypeConstructorFunctions.jl")
include("JoiceHelpFile.jl")
include("JoiceEstimationFunctions.jl")

# function ModelInfo(Name::AbstractString="Joice-Model", )
# end
end

max_stars_repo_name: crforsythe/Joice-Releases

max_issues_repo_name: crforsythe/Joice-Releases

max_forks_repo_name: crforsythe/Joice-Releases


text:
export iar

using IterativeSolvers
using LinearAlgebra
using Random
using Statistics

"""
    iar(nep,[maxit=30,][σ=0,][γ=1,][linsolvecreator=DefaultLinSolverCreator(),][tol=eps()*10000,][neigs=6,][errmeasure,][v=rand(size(nep,1),1),][logger=0,][check_error_every=1,][orthmethod=DGKS,][proj_solve=false,][inner_solver_method=DefaultInnerSolver(),][inner_logger=0])

Run the infinite Arnoldi method on the nonlinear eigenvalue problem stored in `nep`.

The target `σ` is the center around which eiganvalues are computed. The value `γ` corresponds
to scaling and specifying a shift and scaling is effectively the same as the transformation `λ=γs+σ` where `s`
is now the eigenvalue parameter. If you want eigenvalues in a disk centered, select `σ` as the center
of the disk and `γ` as the radius.
The vector `v` is the starting vector for constructing the Krylov space. The orthogonalization
method, used in contructing the orthogonal basis of the Krylov space, is specified by
`orthmethod`, see the package `IterativeSolvers.jl`.
The iteration is continued until `neigs` Ritz pairs have converged.
This function throws a `NoConvergenceException` if the wanted eigenpairs are not computed after `maxit` iterations.
However, if `neigs` is set to `Inf` the iteration is continued until `maxit` iterations without an error being thrown.
The parameter `proj_solve` determines if the Ritz paris are extracted using the Hessenberg matrix (false),
or as the solution to a projected problem (true). If `true`, the method is descided by `inner_solver_method`, and the
logging of the inner solvers are descided by `inner_logger`, which works in the same way as `logger`.

See [`augnewton`](@ref) for other parameters.

# Example
```julia-repl
julia> using NonlinearEigenproblems, LinearAlgebra
julia> nep=nep_gallery("dep0",100);
julia> v0=ones(size(nep,1));
julia> λ,v=iar(nep;v=v0,tol=1e-5,neigs=3);
julia> norm(compute_Mlincomb!(nep,λ[1],v[:,1])) # Is it an eigenvalue?
julia> λ    # print the computed eigenvalues
3-element Array{Complex{Float64},1}:
 -0.15606211475666945 - 0.12273439802763578im
 -0.15606211475666862 + 0.12273439802763489im
  0.23169243065648365 - 9.464790582509696e-17im
```

# References
* Algorithm 2 in Jarlebring, Michiels Meerbergen, A linear eigenvalue algorithm for the nonlinear eigenvalue problem, Numer. Math, 2012
"""
iar(nep::NEP;params...)=iar(ComplexF64,nep;params...)
function iar(
    ::Type{T},
    nep::NEP;
    orthmethod::Type{T_orth}=DGKS,
    maxit=30,
    linsolvercreator=DefaultLinSolverCreator(),
    tol=eps(real(T))*10000,
    neigs=6,
    errmeasure::ErrmeasureType = DefaultErrmeasure(nep),
    σ=zero(T),
    γ=one(T),
    v=randn(real(T),size(nep,1)),
    logger=0,
    check_error_every=1,
    proj_solve=false,
    inner_solver_method=DefaultInnerSolver(),
    inner_logger=0)where{T<:Number,T_orth<:IterativeSolvers.OrthogonalizationMethod}

    @parse_logger_param!(logger)
    @parse_logger_param!(inner_logger)

    # Ensure types σ and v are of type T
    σ=T(σ)
    v=Array{T,1}(v)

    n = size(nep,1);
    m = maxit;

    # initialization
    V = zeros(T,n*(m+1),m+1);
    H = zeros(T,m+1,m);
    y = zeros(T,n,m+1);
    α=Vector{T}(γ.^(0:m)); α[1]=zero(T);
    local M0inv::LinSolver=create_linsolver(linsolvercreator,nep,σ)

    err = NaN*ones(m,m);
    λ=zeros(T,m+1); Q=zeros(T,n,m+1);

    vv=view(V,1:1:n,1); # next vector V[:,k+1]
    vv[:]=v; vv[:]=vv[:]/norm(vv);
    k=1; conv_eig=0;
    local pnep::NEP;
    if (proj_solve)
        pnep=create_proj_NEP(nep);
    end

    while (k <= m) && (conv_eig<neigs)

        VV=view(V,1:1:n*(k+1),1:k); # extact subarrays, memory-CPU efficient
        vv=view(V,1:1:n*(k+1),k+1); # next vector V[:,k+1]

        y[:,2:k+1] = reshape(VV[1:1:n*k,k],n,k);
        broadcast!(/,view(y,:,2:k+1),view(y,:,2:k+1),(1:k)')
        y[:,1] = compute_Mlincomb!(nep,σ,y[:,1:k+1],α[1:k+1]);
        y[:,1] = -lin_solve(M0inv,y[:,1]);

        vv[:]=reshape(y[:,1:k+1],(k+1)*n,1);
        # orthogonalization
        H[k+1,k] = orthogonalize_and_normalize!(VV, vv, view(H,1:k,k), orthmethod)

        # compute Ritz pairs (every check_error_every iterations)
        if ((rem(k,check_error_every)==0)||(k==m))
            # Extract eigenvalues from Hessenberg matrix
            D,Z = eigen(H[1:k,1:k])

            VV = view(V,1:1:n,1:k)
            Q = VV*Z
            λ = σ .+ γ ./ D

            if (proj_solve)  # Projected solve to extract eigenvalues (otw hessenberg matrix)
                QQ,RR=qr(VV); # Project on this space
                QQ = Matrix(QQ)
                set_projectmatrices!(pnep,QQ,QQ);
                # Make a call to the inner solve method
                λproj,Qproj=inner_solve(inner_solver_method,T,pnep,
                                        V=RR*Z,λv=copy(λ),
                                        neigs=k,#size(λ,1)+3,
                                        σ=mean(λ),
                                        inner_logger=inner_logger);

                Q=QQ*Qproj;
                λ=λproj;
             end
            conv_eig=0;
            # compute the errors
            err[k,1:size(λ,1)]=
              map(s-> estimate_error(errmeasure,λ[s],Q[:,s]), 1:size(λ,1))
            # Log them and compute the converged
            push_iteration_info!(logger,2, k,err=err[k,1:size(λ,1)], λ=λ,
                                 continues=true);
            for s=1:size(λ,1)
                if err[k,s]<tol
                    conv_eig=conv_eig+1;
                    push_info!(logger,"+", continues=true);
                elseif err[k,s]<tol*10
                    push_info!(logger,"=", continues=true);
                else
                    push_info!(logger,"-", continues=true);
                end
            end
            push_info!(logger,"");
            # Sort the errors
            idx=sortperm(err[k,1:k]); # sort the error
            err[k,1:k]=err[k,idx];

            # extract the converged Ritzpairs
            if (k==m)||(conv_eig>=neigs)
                nrof_eigs = Int(min(length(λ),neigs))
                λ=λ[idx[1:nrof_eigs]]
                Q=Q[:,idx[1:length(λ)]]
            end

        end

        k=k+1;
    end
    k=k-1

    # NoConvergenceException
    if conv_eig<neigs && neigs != Inf
        # Sort the errors
        idx=sortperm(err[end,1:k]); err=err[end,idx];
        kk=length(λ);
        λ=λ[idx[1:kk]]; Q=Q[:,idx[1:kk]]; err=err[1:kk];
        msg="Number of iterations exceeded. maxit=$(maxit)."
        if conv_eig<3
            msg=string(msg, "Try to change the inner_solver_method for better performance.")
        end
        throw(NoConvergenceException(λ,Q,err,msg))
    end

    # extract the converged Ritzpairs
    λ=λ[1:min(length(λ),conv_eig)];
    Q=Q[:,1:min(size(Q,2),conv_eig)];
    return λ,Q,V[:,1:k]
end


text:
# make this into test script, spline proxy vs lorentzian.

include("../src/NMRSpectraSimulator.jl")
import .NMRSpectraSimulator

using LinearAlgebra
using FFTW
import PyPlot
#import BSON
#import JLD

#import Clustering
import Statistics

import Random
Random.seed!(25)

PyPlot.close("all")
fig_num = 1

PyPlot.matplotlib["rcParams"][:update](["font.size" => 22, "font.family" => "serif"])

### user inputs.

tol_coherence = 1e-2
α_relative_threshold = 0.05
λ0 = 3.4
Δcs_max = 0.2
κ_λ_lb = 0.5
κ_λ_ub = 2.5

molecule_names = ["DSS";]
#molecule_names = ["L-Serine"; "L-Phenylalanine"; "DSS";]
#molecule_names = ["D-(+)-Glucose"; "DSS"]

# # # get the all GISSMO entries.
# import GISSMOReader
# tmp = GISSMOReader.getGISSMOentriesall()
# GISSMO_entries = GISSMOReader.extractfields(tmp, "entry")
# molecule_names = GISSMOReader.extractfields(tmp, "molecule_name")

# machine values taken from the BMRB 700 MHz 20 mM glucose experiment.
fs = 14005.602240896402
SW = 20.0041938620844
ν_0ppm = 10656.011933076665

# # machine values for the BMRB 500 MHz glucose experiment.
# ν_0ppm = 6752.490995937095
# SW = 16.0196917451925
# fs = 9615.38461538462

# path to the GISSMO Julia storage folder.
H_params_path = "/home/roy/Documents/repo/NMRData/input/coupling_info"

### end inputs.

hz2ppmfunc = uu->(uu - ν_0ppm)*SW/fs
ppm2hzfunc = pp->(ν_0ppm + pp*fs/SW)

Δcs_max_mixture = collect( Δcs_max for i = 1:length(molecule_names))

println("Timing: setupmixtureproxies()")
dummy_SSFID = NMRSpectraSimulator.SpinSysFIDType2(0.0)
@time mixture_params = NMRSpectraSimulator.setupmixtureproxies(molecule_names,
    H_params_path, ppm2hzfunc, fs, SW, λ0,
    ν_0ppm, dummy_SSFID;
    tol_coherence = tol_coherence,
    α_relative_threshold = α_relative_threshold)
As = mixture_params

#@assert 1==43

u_min = ppm2hzfunc(-0.5)
u_max = ppm2hzfunc(4.0)

NMRSpectraSimulator.fitproxies!(As;
    κ_λ_lb = κ_λ_lb,
    κ_λ_ub = κ_λ_ub,
    u_min = u_min,
    u_max = u_max,
    Δr = 1.0,
    Δκ_λ = 0.05)

### plot.

# purposely distort the spectra by setting non-autophased FID values.
compound_select = 1
spin_sys_select = 1
Ag = As[compound_select]
Ag.ss_params.d[spin_sys_select] = rand(length(Ag.part_inds_compound[spin_sys_select]))
Ag.ss_params.κs_λ[spin_sys_select] = rand(length(Ag.part_inds_compound[spin_sys_select])) .+ 1
Ag.ss_params.κs_β = collect( rand(length(Ag.ss_params.κs_β[i])) .* (2*π) for i = 1:length(Ag.ss_params.κs_β) )

f = uu->NMRSpectraSimulator.evalmixture(uu, mixture_params)

# test params.
ΩS_ppm = collect( hz2ppmfunc.( NMRSpectraSimulator.combinevectors(A.Ωs) ./ (2*π) ) for A in mixture_params )
#ΩS_ppm_flat = NMRSpectraSimulator.combinevectors(ΩS_ppm)

P = LinRange(hz2ppmfunc(u_min), hz2ppmfunc(u_max), 50000)
U = ppm2hzfunc.(P)
U_rad = U .* (2*π)

## parameters that affect qs.
# A.d, A.κs_λ, A.κs_β
# A.d_singlets, A.αs_singlets, A.Ωs_singlets, A.β_singlets, A.λ0, A.κs_λ_singlets
q = uu->NMRSpectraSimulator.evalitpproxymixture(uu, mixture_params)

f_U = f.(U_rad)
q_U = q.(U_rad)

discrepancy = abs.(f_U-q_U)
max_val, ind = findmax(discrepancy)
println("relative discrepancy = ", norm(discrepancy)/norm(f_U))
println("max discrepancy: ", max_val)
println()

## visualize.
PyPlot.plot(P, real.(f_U), label = "f")
PyPlot.plot(P, real.(q_U), label = "q")

PyPlot.legend()
PyPlot.xlabel("ppm")
PyPlot.ylabel("real")
PyPlot.title("f vs q")


text:
function TricubicInterpolator()
    return TricubicInterpolator(())
end

function interpolate(obj::TricubicInterpolator, arg0::Vector{jdouble}, arg1::Vector{jdouble}, arg2::Vector{jdouble}, arg3::Vector{Vector{jdouble}})
    return jcall(obj, "interpolate", TricubicInterpolatingFunction, (Vector{jdouble}, Vector{jdouble}, Vector{jdouble}, Vector{Vector{jdouble}}), arg0, arg1, arg2, arg3)
end


text:
function MultivariateNormalMixtureExpectationMaximization(arg0::Vector{Vector{jdouble}})
    return MultivariateNormalMixtureExpectationMaximization((Vector{Vector{jdouble}},), arg0)
end

function estimate(::Type{MultivariateNormalMixtureExpectationMaximization}, arg0::Vector{Vector{jdouble}}, arg1::jint)
    return jcall(MultivariateNormalMixtureExpectationMaximization, "estimate", MixtureMultivariateNormalDistribution, (Vector{Vector{jdouble}}, jint), arg0, arg1)
end

function fit(obj::MultivariateNormalMixtureExpectationMaximization, arg0::MixtureMultivariateNormalDistribution)
    return jcall(obj, "fit", void, (MixtureMultivariateNormalDistribution,), arg0)
end

function fit(obj::MultivariateNormalMixtureExpectationMaximization, arg0::MixtureMultivariateNormalDistribution, arg1::jint, arg2::jdouble)
    return jcall(obj, "fit", void, (MixtureMultivariateNormalDistribution, jint, jdouble), arg0, arg1, arg2)
end

function get_fitted_model(obj::MultivariateNormalMixtureExpectationMaximization)
    return jcall(obj, "getFittedModel", MixtureMultivariateNormalDistribution, ())
end

function get_log_likelihood(obj::MultivariateNormalMixtureExpectationMaximization)
    return jcall(obj, "getLogLikelihood", jdouble, ())
end