text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct FirewallRulesApi <: SwaggerApi
    client::Swagger.Client
end

"""

Creates or updates the specified firewall rule. During update, the firewall rule with the specified name will be replaced with this new firewall rule.
Param: subscriptionId::String (required)
Param: resourceGroupName::String (required)
Param: accountName::String (required)
Param: firewallRuleName::String (required)
Param: parameters::CreateOrUpdateFirewallRuleParameters (required)
Param: api_version::String (required)
Return: FirewallRule
"""
function _swaggerinternal_firewallRulesCreateOrUpdate(_api::FirewallRulesApi, subscriptionId::String, resourceGroupName::String, accountName::String, firewallRuleName::String, parameters, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PUT", FirewallRule, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.DataLakeStore/accounts/{accountName}/firewallRules/{firewallRuleName}", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "accountName", accountName)  # type String
    Swagger.set_param(_ctx.path, "firewallRuleName", firewallRuleName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function firewallRulesCreateOrUpdate(_api::FirewallRulesApi, subscriptionId::String, resourceGroupName::String, accountName::String, firewallRuleName::String, parameters, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_firewallRulesCreateOrUpdate(_api, subscriptionId, resourceGroupName, accountName, firewallRuleName, parameters, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function firewallRulesCreateOrUpdate(_api::FirewallRulesApi, response_stream::Channel, subscriptionId::String, resourceGroupName::String, accountName::String, firewallRuleName::String, parameters, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_firewallRulesCreateOrUpdate(_api, subscriptionId, resourceGroupName, accountName, firewallRuleName, parameters, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Deletes the specified firewall rule from the specified Data Lake Store account.
Param: subscriptionId::String (required)
Param: resourceGroupName::String (required)
Param: accountName::String (required)
Param: firewallRuleName::String (required)
Param: api_version::String (required)
Return: Nothing
"""
function _swaggerinternal_firewallRulesDelete(_api::FirewallRulesApi, subscriptionId::String, resourceGroupName::String, accountName::String, firewallRuleName::String, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "DELETE", Nothing, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.DataLakeStore/accounts/{accountName}/firewallRules/{firewallRuleName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "accountName", accountName)  # type String
    Swagger.set_param(_ctx.path, "firewallRuleName", firewallRuleName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function firewallRulesDelete(_api::FirewallRulesApi, subscriptionId::String, resourceGroupName::String, accountName::String, firewallRuleName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_firewallRulesDelete(_api, subscriptionId, resourceGroupName, accountName, firewallRuleName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function firewallRulesDelete(_api::FirewallRulesApi, response_stream::Channel, subscriptionId::String, resourceGroupName::String, accountName::String, firewallRuleName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_firewallRulesDelete(_api, subscriptionId, resourceGroupName, accountName, firewallRuleName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Gets the specified Data Lake Store firewall rule.
Param: subscriptionId::String (required)
Param: resourceGroupName::String (required)
Param: accountName::String (required)
Param: firewallRuleName::String (required)
Param: api_version::String (required)
Return: FirewallRule
"""
function _swaggerinternal_firewallRulesGet(_api::FirewallRulesApi, subscriptionId::String, resourceGroupName::String, accountName::String, firewallRuleName::String, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", FirewallRule, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.DataLakeStore/accounts/{accountName}/firewallRules/{firewallRuleName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "accountName", accountName)  # type String
    Swagger.set_param(_ctx.path, "firewallRuleName", firewallRuleName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function firewallRulesGet(_api::FirewallRulesApi, subscriptionId::String, resourceGroupName::String, accountName::String, firewallRuleName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_firewallRulesGet(_api, subscriptionId, resourceGroupName, accountName, firewallRuleName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function firewallRulesGet(_api::FirewallRulesApi, response_stream::Channel, subscriptionId::String, resourceGroupName::String, accountName::String, firewallRuleName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_firewallRulesGet(_api, subscriptionId, resourceGroupName, accountName, firewallRuleName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Lists the Data Lake Store firewall rules within the specified Data Lake Store account.
Param: subscriptionId::String (required)
Param: resourceGroupName::String (required)
Param: accountName::String (required)
Param: api_version::String (required)
Return: FirewallRuleListResult
"""
function _swaggerinternal_firewallRulesListByAccount(_api::FirewallRulesApi, subscriptionId::String, resourceGroupName::String, accountName::String, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", FirewallRuleListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.DataLakeStore/accounts/{accountName}/firewallRules", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "accountName", accountName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function firewallRulesListByAccount(_api::FirewallRulesApi, subscriptionId::String, resourceGroupName::String, accountName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_firewallRulesListByAccount(_api, subscriptionId, resourceGroupName, accountName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function firewallRulesListByAccount(_api::FirewallRulesApi, response_stream::Channel, subscriptionId::String, resourceGroupName::String, accountName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_firewallRulesListByAccount(_api, subscriptionId, resourceGroupName, accountName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Updates the specified firewall rule.
Param: subscriptionId::String (required)
Param: resourceGroupName::String (required)
Param: accountName::String (required)
Param: firewallRuleName::String (required)
Param: api_version::String (required)
Param: parameters::UpdateFirewallRuleParameters
Return: FirewallRule
"""
function _swaggerinternal_firewallRulesUpdate(_api::FirewallRulesApi, subscriptionId::String, resourceGroupName::String, accountName::String, firewallRuleName::String, api_version::String; parameters=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PATCH", FirewallRule, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.DataLakeStore/accounts/{accountName}/firewallRules/{firewallRuleName}", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "accountName", accountName)  # type String
    Swagger.set_param(_ctx.path, "firewallRuleName", firewallRuleName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function firewallRulesUpdate(_api::FirewallRulesApi, subscriptionId::String, resourceGroupName::String, accountName::String, firewallRuleName::String, api_version::String; parameters=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_firewallRulesUpdate(_api, subscriptionId, resourceGroupName, accountName, firewallRuleName, api_version; parameters=parameters, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function firewallRulesUpdate(_api::FirewallRulesApi, response_stream::Channel, subscriptionId::String, resourceGroupName::String, accountName::String, firewallRuleName::String, api_version::String; parameters=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_firewallRulesUpdate(_api, subscriptionId, resourceGroupName, accountName, firewallRuleName, api_version; parameters=parameters, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export firewallRulesCreateOrUpdate, firewallRulesDelete, firewallRulesGet, firewallRulesListByAccount, firewallRulesUpdate


text:
function reshape_product(A::AbstractArray)
    k = 0
    v = Vector{Vector{Int64}}(undef, prod(size(A)))
    for i in eachindex(A)
        v[i] = [A[i]...]
    end
    return sort!(v)
end

function matchrow(a::AbstractVector, B::AbstractArray)
    findfirst(i -> all(j -> a[j] == B[i,j], 1:size(B,2)), 1:size(B,1))
end

function circlecoord(x::Float64, y::Float64, r::Float64)
    radians = -0.1:0.1:2pi
    cx = [cos(rad)r + x for rad in radians]
    cy = [sin(rad)r + y for rad in radians]
    cx, cy
end

function rerr_approx(u::Vector)
    lb = round(Int, 0.9 * length(u))
    u_approx = mean(u[lb:end])
    um = cumsum(u) ./ collect(1:T)
    u_rerr_approx = abs.(um .- u_approx) ./ abs.(u_approx)
    return u_rerr_approx
end

rerr(u::Number, u_true::Number) = abs(u - u_true) / abs(u_true)

tosecond(t::T) where {T} = t / convert(T, Dates.Second(1))
tominute(t::T) where {T} = t / convert(T, Dates.Minute(1))
tohour(t::T) where {T} = t / convert(T, Dates.Hour(1))

ext: jl

max_stars_repo_name: ajkeith/Cyber-Air-Defense

max_issues_repo_name: ajkeith/Cyber-Air-Defense


text:
## Projection methods for global solution for non-linear models
## using neural network to approximate the policy function
## Based on Computationa Macroeconomics for the open economy
## Lim and McNelis 2008
## Chapter 2: Flexible Price Model
using Parameters, Distributions, Optim, Plots
plotly()

type Parameter
    Θ_s::Array
    Θ_rf::Array
    Θ_ss::Array
end

type Variable
    endo::Array
    state_endo::Array
    euler_error::Array
    other::Array
    shock::Array
end

type Descriptive
    T::Int64
    K::Int64
    B::Int64
    State::Int64
    Euler::Int64
    Neuron::Int64
    Neuronx1::Int64
    Param::Int64
    Start::Int64
    Other::Int
end

type Model
    Θ::Parameter
    Var::Variable
    Des::Descriptive
    d::Distribution
end

function Model(Θ::Parameter,Des::Descriptive,d::Distribution)
    Var = initialize_system(Θ,Des,d)
    return Model(Θ::Parameter,Var::Variable, Des::Descriptive,d::Distribution)
end

function initialize_system(Θ::Parameter,Des::Descriptive,d::Distribution)
    @unpack T, K, B, State, Euler, Neuron, Neuronx1, Param, Start,Other = Des
    Θ_ss = Θ.Θ_ss
    (Blam_ss,C_ss,F_ss,G_ss,K_ss,L_ss,P_ss,Pk_ss,R_ss,S_ss,W_ss,X_ss,Y_ss,Z_ss) = Θ_ss
    # Endo
    endo = initialize_endo(T,K,Θ_ss)
    # State
    state_endo = initialize_state(T,K,State)
    # Errors
    euler_error = initialize_error(T,K,Euler)
    #Other
    other = initialize_other(T,K,Other)
    return Variable(endo,state_endo,euler_error,other,Θ,Des,d)
end

initialize_endo(T,K,Θ_ss) = [initialize_var(T,K,i)  for i in Θ_ss]

function initialize_var(T::Int64,K::Int64,x_ss::Float64)
    x = x_ss*ones(T,K)
end

initialize_state(T,K,State) = [zeros(T,K)  for i in 1:State]
initialize_error(T,K,Euler) = [zeros(T,K)  for i in 1:Euler+1]
initialize_other(T,K,Other) = [zeros(T,K)  for i in 1:Other]

function Variable(endo::Array,state_endo::Array,euler_error::Array,other::Array,Θ::Parameter,Des::Descriptive,d::Distribution)
    @unpack T, K, B, State, Euler, Neuron, Neuronx1, Param, Start,Other = Des
    @unpack Θ_s, Θ_rf, Θ_ss = Θ
    nshock = get_shock(Θ_rf)
    vStd = get_std(Θ_rf,nshock)
    shock = initialize_shock(T,K,vStd,d)
    return  Variable(endo::Array,state_endo::Array,euler_error::Array,other::Array,shock::Array)
end

get_shock(Θ_rf::Array) = Int(length(Θ_rf)/2)
get_std(Θ_rf::Array,nshock::Int64) = Θ_rf[nshock+1:end]

initialize_shock(T::Int64,K::Int64,vStd::Array,d::Distribution) = [rand(d,T,K).*σ for σ in vStd]

function demean(x::Float64,x_ss::Float64)
    x_dm = x - x_ss
end

function ar_process(x_l::Float64,x_ss::Float64,ρ::Float64,shock::Float64)
    x = ρ*log(x_l) + (1-ρ)*log(x_ss) + shock
    return exp(x)
end

function create_neuron(xstate,gamax,jj,kk,Euler,Neuron)
    neuron = zeros(1,Euler*Neuron)
    for nn = 1: Euler*Neuron
        neuron[1,nn] = 1./(1 + exp.(-xstate*gamax[jj[nn]:kk[nn]]))[1] - 0.5
    end
    return neuron
end

function create_endogenous!(m::Model,gamax::Array,jj,kk)
    @unpack T, K, B, State, Euler, Neuron, Neuronx1, Param, Start,Other = m.Des
    @unpack Θ_s, Θ_rf, Θ_ss = m.Θ
    @unpack endo, state_endo, euler_error, other, shock = m.Var
    (Blam_ss,C_ss,F_ss,G_ss,K_ss,L_ss,P_ss,Pk_ss,R_ss,S_ss,W_ss,X_ss,Y_ss,Z_ss) = Θ_ss
    (Rstar,PFstar,PXstar,eta,omega,beta,alpha1,kappa,chis,phi0,phi1)= Θ_s
    (ρ, σ) = Θ_rf
    (Blam,C,F,G,κ,L,P,Pk,R,S,W,X,Y,Z) = endo
    (ZZ, FF, RR) =  state_endo
    (error_C, error_S) = euler_error
    (Zrisk, Zinf, trade, Zder) = other
    zshock = shock[1]
    for j = 1:K
        for i = (Start+1):T
            # Create the log AR(1) shock process and exponenciate it
            Z[i,j] = ar_process(Z[i-1,j],Z_ss,ρ,zshock[i,j])
            # Demeaning the exogenous and endogenous state variables
            ZZ[i,j] = demean(Z[i,j],Z_ss)
            FF[i,j] = demean(F[i - 1,j],F_ss)
            RR[i,j] = demean(R[i - 1,j],R_ss)
            xstate = [ZZ[i,j] FF[i,j] RR[i,j]]
            # Setting the approximating functions
            neuron = create_neuron(xstate,gamax,jj,kk,Euler,Neuron)
            pea_C = neuron[1, 1: Neuron][1]
            pea_S = neuron[1, Neuron + 1:2*Neuron][1]
            C[i,j] = exp(pea_C)*C_ss # L: WHY EXPONENCIATING?
            S[i,j] = exp(pea_S)*S_ss
            # Generating the endogeous variables
            Y[i,j] = C[i,j] + G_ss + X_ss
            LL = 0.5*(1 - alpha1)*(Z[i,j]^kappa)*(Y[i,j]^(1 - kappa))*(C[i,j]^-eta)
            L[i,j] = LL^(1/(1 - kappa + omega))
            L[i,j] = real(L[i,j])
            KK = ((Y[i,j]/Z[i,j])^kappa) - (1 - alpha1)*L[i,j]^kappa
            κ[i,j] = (KK/alpha1)^(1/kappa)
            κ[i,j] = real(κ[i,j])
            mpl = (1 - alpha1)*(Z[i,j]^kappa)*(Y[i,j]/L[i,j])^(1 - kappa)
            mpk = (alpha1)*(Z[i,j]^kappa)*(Y[i,j]/κ[i,j])^(1 - kappa)
            mpl = real(mpl)
            mpk = real(mpk)
            Pk[i,j] = S[i,j]*PFstar
            P[i,j] = 2*Pk[i,j]/mpk
            W[i,j] = (L[i,j]^omega)*(C[i,j]^eta)*P[i,j]
            W[i,j] = real(W[i,j])
            Zinf[i,j] = 0.25*((P[i,j]/P[i - 4,j]) - 1)
            R[i,j] = phi0*R[i - 1,j] + (1 - phi0)*(Rstar + phi1*Zinf[i,j])
            trade[i,j] = P[i,j]*X_ss - S[i,j]*PFstar*κ[i,j]
            trade1 = trade[i,j]/S[i,j]
            F[i,j] = F[i - 1,j]*(1 + Rstar + Zrisk[i - 1,j]) - trade1
            Blam[i,j] = (C[i,j]^-eta)/P[i,j]
            Blam[i,j] = real(Blam[i,j])
            Zrisk[i,j] = sign(F[i - 1,j])*chis*(exp(abs(F[i - 1,j])) - 1)
            Zder[i,j] = chis*(exp(abs(F[i - 1,j])))
            # Obtaining the Euler errors
            MUC = Blam[i,j]*(beta*(1 + R[i - 1,j]))
            MUCLAG = Blam[i - 1,j]
            error_C[i,j] = (MUC/MUCLAG) - 1
            MUS = S[i,j]*(1 + Rstar + Zrisk[i - 1,j] + Zder[i - 1,j]*F[i - 1,j])
            MUSLAG = (1 + R[i - 1,j])*S[i - 1,j]
            error_S[i,j] = (MUS/MUSLAG) - 1
        end
    end
    endo = [Blam,C,F,G,κ,L,P,Pk,R,S,W,X,Y,Z]
    state_endo = [ZZ, FF, RR]
    euler_error = [error_C, error_S]
    other = [Zrisk, Zinf, trade, Zder]
    shock[1] = zshock
    Var = Variable(endo, state_endo, euler_error, other, shock)
    return m = Model(m.Θ,Var,m.Des,m.d)
end

minimize_euler_error!(gamax::Array) = minimize_euler_error!(gamax,m)

function minimize_euler_error!(gamax::Array,m::Model)
    @unpack T, K, B, State, Euler, Neuron, Neuronx1, Param, Start, Other = m.Des
    jk = Param
    jj = 1:State:Param
    kk = State:State:Param
    # Recursively create the endogenous variables
    m = create_endogenous!(m,gamax,jj,kk)
    # Defining the errors function to be minimized
    iError = create_euler_error(m)
end

function create_euler_error(m::Model)
    @unpack endo, state_endo, euler_error, other, shock = m.Var
    @unpack T, K, B, State, Euler, Neuron, Neuronx1, Param, Start,Other = m.Des
    mError = reshape_error(euler_error,T*K,1)
    iError = mean_squared_error(mError)
end

reshape_error(euler_error::Array,T::Int64,K::Int64) = hcat([reshape(error,T*K,1) for error in euler_error]...)

mean_squared_error(mError::Array) = mean([sum(mError[i,:])^2 for i in 1:size(mError,1)])

function evaluate_model!(gamax::Array,m::Model)
    @unpack T, K, B, State, Euler, Neuron, Neuronx1, Param, Start,Other = m.Des
    jk = Param
    jj = 1:State:Param
    kk = State:State:Param
    # Recursively create the endogenous variables
    m = create_endogenous!(m,gamax,jj,kk)
end

function burn_in!(endo::Array, state_endo::Array, euler_error::Array, other::Array, shock::Array, B::Int64)
    endo = [burn_in(i,B) for i in endo]
    state_endo = [burn_in(i,B) for i in state_endo]
    euler_error = [burn_in(i,B) for i in euler_error]
    other = [burn_in(i,B) for i in other]
    shock = [burn_in(i,B) for i in shock]
    return Var = Variable(endo, state_endo, euler_error, other, shock)
end

burn_in(x::Array,B::Int64) = x[B+1:end,:]

function plot_irf!(m::Model,magnitude::Float64,horizon::Int64,timing::Int64)
    compute_irf!(m,magnitude,horizon,timing)
    @unpack endo, state_endo, euler_error, other, shock = m.Var
    (Blam,C,F,G,κ,L,P,Pk,R,S,W,X,Y,Z) = endo
    plot_range = timing:timing+horizon-1
    # Plot
    irf = plot(layout = grid(2,5), legend = false)
    plot!(mean(Z[plot_range,:],2), linewidth = 2, subplot = 1, title = "Z")
    plot!(mean(C[plot_range,:],2), linewidth = 2, subplot = 2, title = "C")
    plot!(mean(S[plot_range,:],2), linewidth = 2, subplot = 3, title = "S")
    plot!(mean(Y[plot_range,:],2), linewidth = 2, subplot = 4, title = "Y")
    plot!(mean(κ[plot_range,:],2), linewidth = 2, subplot = 5, title = "K")
    plot!(mean(L[plot_range,:],2), linewidth = 2, subplot = 6, title = "L")
    plot!(mean(P[plot_range,:],2), linewidth = 2, subplot = 7, title = "P")
    plot!(mean(R[plot_range,:],2), linewidth = 2, subplot = 8, title = "R")
    plot!(mean(F[plot_range,:],2), linewidth = 2, subplot = 9, title = "F")
    plot!(mean(W[plot_range,:],2)./P[plot_range], linewidth = 2, subplot = 10, title = "W/P")
    gui()
end

function compute_irf!(m::Model,magnitude::Float64,horizon::Int64,timing::Int64)
    @unpack T, K, B, State, Euler, Neuron, Neuronx1, Param, Start,Other = m.Des
    vShock = one_time_shock(magnitude,timing,T,K)
    m.Var.shock = vShock
    m = evaluate_model!(vParam_star,m)
end

one_time_shock(magnitude::Float64,timing::Int64,T::Int64,K::Int64) = [[zeros(timing-1,K); magnitude*ones(1,K) ; zeros(T-timing,K)]]

function judd_gaspar_stat(m::Model)
    # I AM NOT SURE THIS IS CORRECT
    @unpack endo, state_endo, euler_error, other, shock = m.Var
    (Blam,C,F,G,κ,L,P,Pk,R,S,W,X,Y,Z) = endo
    (error_C, error_S) = euler_error
    jg_c = mean(abs.(error_C)./C,2)
    jg_s = mean(abs.(error_S)./S,2)
    jg = histogram(layout = grid(2,1), legend = false)
    histogram!(jg_c, subplot = 1, title = "Error C")
    histogram!(jg_s, subplot = 2, title = "Error S")
end

function haan_marcet_statistics(m::Model)
    # I AM NOT SURE THIS IS CORRECT
    @unpack T, K, B, State, Euler, Neuron, Neuronx1, Param, Start,Other = m.Des
    @unpack endo, state_endo, euler_error, other, shock = m.Var
    DM = zeros(K,Euler)
    for j in 1:Euler
        dm = zeros(K)
        for i in 1:K
            e = euler_error[j][:,i]
            x = hcat([state_endo[j][:,i] for j in 1:State]...)
            Q = (1/T)*e'*x
            A = (1/T)sum(x*x'*e.^2)
            dm[i] = T*Q*A^(-1)*Q'
        end
        DM[:,j] = dm
    end
    return sort(DM,1)
end

function plot_dm_stat(DM::Array)
    # I AM NOT SURE THIS IS CORRECT
    K = size(DM,2)
    dm = plot(layout = grid(K,1), legend = false, show = true)
    [plot!(DM[:,i], subplot = i) for i in 1:K]
end


text:
import AtmosphericProfilesLibrary
const APL = AtmosphericProfilesLibrary

function units(var::String)
    var == "θ_liq_ice" && return "[K]"
    var == "T" && return "[K]"
    var == "p" && return "[kPa]"
    var == "q_tot" && return "[kg/kg]"
    var == "u" && return "[m/s]"
    var == "v" && return "[m/s]"
    var == "ug" && return "[m/s]"
    var == "vg" && return "[m/s]"
    var == "u0" && return "[m/s]"
    var == "v0" && return "[m/s]"
    var == "shf" && return "[W/m^2]"
    var == "lhf" && return "[W/m^2]"
    var == "tke" && return "[m^2/s^2]"
    var == "dTdt" && return "[K/s]"
    var == "dqtdt" && return "[kg/(kg s)]"
    var == "area" && return "[1]"
    var == "w" && return "[m/s]"
    var == "RH" && return "[%?]"
    var == "subsidence" && return "[?]"
    error("No units found for variable $var")
end

include("define_save_plots.jl")

const zr_hi = range(0, stop = 2.5e4, length = 100);
const zr_lo = range(0, stop = 4000.0, length = 100);
const t_range = range(0, stop = 3600*15, length = 200);

#####
##### z-t profiles
#####
tz_profiles = [
    (; func = APL.TRMM_LBA_radiation, kwargs = (;z_range = zr_hi, t_range, label = "radiation")),
    (; func = APL.ARM_SGP_dTdt      , kwargs = (;z_range = zr_lo, t_range, label = "dTdt", units = "[K/s]")),
]

#####
##### z profiles
#####
z_profiles = [
    (; func = APL.Soares_q_tot         , kwargs = (;z_range = zr_lo, xlabel = "q_tot")),
    (; func = APL.Soares_θ_liq_ice     , kwargs = (;z_range = zr_lo, xlabel = "θ_liq_ice")),
    (; func = APL.Soares_u             , kwargs = (;z_range = zr_lo, xlabel = "u")),
    (; func = APL.Soares_tke           , kwargs = (;z_range = zr_lo, xlabel = "tke")),
    (; func = APL.Nieuwstadt_θ_liq_ice , kwargs = (;z_range = zr_lo, xlabel = "θ_liq_ice")),
    (; func = APL.Nieuwstadt_u         , kwargs = (;z_range = zr_lo, xlabel = "u")),
    (; func = APL.Nieuwstadt_tke       , kwargs = (;z_range = zr_lo, xlabel = "tke")),
    (; func = APL.Bomex_q_tot          , kwargs = (;z_range = zr_lo, xlabel = "q_tot")),
    (; func = APL.Bomex_θ_liq_ice      , kwargs = (;z_range = zr_lo, xlabel = "θ_liq_ice")),
    (; func = APL.Bomex_u              , kwargs = (;z_range = zr_lo, xlabel = "u")),
    (; func = APL.Bomex_tke            , kwargs = (;z_range = zr_lo, xlabel = "tke")),
    (; func = APL.Bomex_geostrophic_u  , kwargs = (;z_range = zr_lo, xlabel = "u")),
    (; func = APL.Bomex_geostrophic_v  , kwargs = (;z_range = zr_lo, xlabel = "v")),
    # (; func = APL.Bomex_dTdt           , kwargs = (;z_range = zr_lo, xlabel = "dTdt")), # Callable by (Π, z)
    (; func = APL.Bomex_dqtdt          , kwargs = (;z_range = zr_lo, xlabel = "dqtdt")),
    (; func = APL.Bomex_subsidence     , kwargs = (;z_range = zr_lo, xlabel = "subsidence")),

    (; func = APL.LifeCycleTan2018_θ_liq_ice  , kwargs = (;z_range = zr_lo, xlabel = "θ_liq_ice")),
    (; func = APL.LifeCycleTan2018_q_tot      , kwargs = (;z_range = zr_lo, xlabel = "q_tot")),
    (; func = APL.LifeCycleTan2018_u          , kwargs = (;z_range = zr_lo, xlabel = "u")),
    (; func = APL.LifeCycleTan2018_tke        , kwargs = (;z_range = zr_lo, xlabel = "tke")),
    (; func = APL.LifeCycleTan2018_geostrophic_u, kwargs = (;z_range = zr_lo, xlabel = "u")),
    (; func = APL.LifeCycleTan2018_geostrophic_v, kwargs = (;z_range = zr_lo, xlabel = "v")),
    # (; func = APL.LifeCycleTan2018_dTdt       , kwargs = (;z_range = zr_lo, xlabel = "dTdt")), # callable by (Π, z)
    (; func = APL.LifeCycleTan2018_dqtdt      , kwargs = (;z_range = zr_lo, xlabel = "dqtdt")),
    (; func = APL.LifeCycleTan2018_subsidence , kwargs = (;z_range = zr_lo, xlabel = "subsidence")),

    (; func = APL.Rico_u               , kwargs = (;z_range = zr_lo, xlabel = "u")),
    (; func = APL.Rico_v               , kwargs = (;z_range = zr_lo, xlabel = "v")),
    (; func = APL.Rico_θ_liq_ice       , kwargs = (;z_range = zr_lo, xlabel = "θ_liq_ice")),
    (; func = APL.Rico_q_tot           , kwargs = (;z_range = zr_lo, xlabel = "q_tot")),
    (; func = APL.Rico_geostrophic_ug  , kwargs = (;z_range = zr_lo, xlabel = "u")),
    (; func = APL.Rico_geostrophic_vg  , kwargs = (;z_range = zr_lo, xlabel = "v")),
    # (; func = APL.Rico_dTdt            , kwargs = (;z_range = zr_lo, xlabel = "dTdt")), # callable by (Π, z)
    (; func = APL.Rico_dqtdt           , kwargs = (;z_range = zr_lo, xlabel = "dqtdt")),
    (; func = APL.Rico_subsidence      , kwargs = (;z_range = zr_lo, xlabel = "v")),

    (; func = APL.TRMM_LBA_p           , kwargs = (;z_range = zr_hi, xlabel = "p")),
    (; func = APL.TRMM_LBA_T           , kwargs = (;z_range = zr_hi, xlabel = "T")),
    (; func = APL.TRMM_LBA_RH          , kwargs = (;z_range = zr_hi, xlabel = "RH")),
    (; func = APL.TRMM_LBA_u           , kwargs = (;z_range = zr_hi, xlabel = "u")),
    (; func = APL.TRMM_LBA_v           , kwargs = (;z_range = zr_hi, xlabel = "v")),
    (; func = APL.TRMM_LBA_tke         , kwargs = (;z_range = zr_hi, xlabel = "tke")),

    (; func = APL.ARM_SGP_u            , kwargs = (;z_range = zr_lo, xlabel = "u")),
    (; func = APL.ARM_SGP_θ_liq_ice    , kwargs = (;z_range = zr_lo, xlabel = "θ_liq_ice")),
    (; func = APL.ARM_SGP_q_tot        , kwargs = (;z_range = zr_lo, xlabel = "q_tot")),
    (; func = APL.ARM_SGP_tke          , kwargs = (;z_range = zr_lo, xlabel = "tke")),
    # (; func = APL.ARM_SGP_dqtdt        , kwargs = (;z_range = zr_lo, xlabel = "dqdt")), # callable by `(Π, t, z)`

    (; func = APL.GATE_III_q_tot       , kwargs = (;z_range = zr_lo, xlabel = "q_tot")),
    (; func = APL.GATE_III_u           , kwargs = (;z_range = zr_lo, xlabel = "u")),
    (; func = APL.GATE_III_T           , kwargs = (;z_range = zr_lo, xlabel = "T")),
    (; func = APL.GATE_III_tke         , kwargs = (;z_range = zr_lo, xlabel = "tke")),
    (; func = APL.GATE_III_dqtdt       , kwargs = (;z_range = zr_lo, xlabel = "dqtdt")),
    (; func = APL.GATE_III_dTdt        , kwargs = (;z_range = zr_lo, xlabel = "dTdt")),

    (; func = APL.Dycoms_RF01_θ_liq_ice  , kwargs = (;z_range = zr_lo, xlabel = "θ_liq_ice")),
    (; func = APL.Dycoms_RF01_q_tot      , kwargs = (;z_range = zr_lo, xlabel = "q_tot")),
    (; func = APL.Dycoms_RF01_u0         , kwargs = (;z_range = zr_lo, xlabel = "u0")),
    (; func = APL.Dycoms_RF01_v0         , kwargs = (;z_range = zr_lo, xlabel = "v0")),
    (; func = APL.Dycoms_RF01_tke        , kwargs = (;z_range = zr_lo, xlabel = "tke")),

    (; func = APL.Dycoms_RF02_θ_liq_ice  , kwargs = (;z_range = zr_lo, xlabel = "θ_liq_ice")),
    (; func = APL.Dycoms_RF02_q_tot      , kwargs = (;z_range = zr_lo, xlabel = "q_tot")),
    (; func = APL.Dycoms_RF02_u          , kwargs = (;z_range = zr_lo, xlabel = "u")),
    (; func = APL.Dycoms_RF02_v          , kwargs = (;z_range = zr_lo, xlabel = "v")),
    (; func = APL.Dycoms_RF02_tke        , kwargs = (;z_range = zr_lo, xlabel = "tke")),

    (; func = APL.GABLS_u              , kwargs = (;z_range = zr_lo, xlabel = "u")),
    (; func = APL.GABLS_v              , kwargs = (;z_range = zr_lo, xlabel = "v")),
    (; func = APL.GABLS_θ_liq_ice      , kwargs = (;z_range = zr_lo, xlabel = "θ_liq_ice")),
    (; func = APL.GABLS_q_tot          , kwargs = (;z_range = zr_lo, xlabel = "q_tot")),
    (; func = APL.GABLS_tke            , kwargs = (;z_range = zr_lo, xlabel = "tke")),
    (; func = APL.GABLS_geostrophic_ug , kwargs = (;z_range = zr_lo, xlabel = "u")),
    (; func = APL.GABLS_geostrophic_vg , kwargs = (;z_range = zr_lo, xlabel = "v")),

    (; func = APL.SP_u                 , kwargs = (;z_range = zr_lo, xlabel = "u")),
    (; func = APL.SP_v                 , kwargs = (;z_range = zr_lo, xlabel = "v")),
    (; func = APL.SP_θ_liq_ice         , kwargs = (;z_range = zr_lo, xlabel = "θ_liq_ice")),
    (; func = APL.SP_q_tot             , kwargs = (;z_range = zr_lo, xlabel = "q_tot")),
    (; func = APL.SP_geostrophic_u     , kwargs = (;z_range = zr_lo, xlabel = "ug")),
    (; func = APL.SP_geostrophic_v     , kwargs = (;z_range = zr_lo, xlabel = "vg")),

    (; func = APL.DryBubble_θ_liq_ice    , kwargs = (;z_range = zr_lo, xlabel = "θ_liq_ice")),
    (; func = APL.DryBubble_updrafts_θ_liq_ice , kwargs = (;z_range = zr_lo, xlabel = "θ_liq_ice")),
    (; func = APL.DryBubble_updrafts_area      , kwargs = (;z_range = zr_lo, xlabel = "area")),
    (; func = APL.DryBubble_updrafts_w         , kwargs = (;z_range = zr_lo, xlabel = "w")),
    (; func = APL.DryBubble_updrafts_T         , kwargs = (;z_range = zr_lo, xlabel = "T")),

]

#####
##### t profiles
#####
t_profiles = [
    (; func = APL.ARM_SGP_shf         , kwargs = (;t_range, label = "shf")),
    (; func = APL.ARM_SGP_lhf         , kwargs = (;t_range, label = "lhf")),
]

for profile in t_profiles
    # We don't have t-profiles yet.
    # @info "Saving t-profile $(nameof(profile.func))" # for debugging
    save_t_profile(profile.func; profile.kwargs...)
end
for profile in tz_profiles
    # @info "Saving tz-profile $(nameof(profile.func))" # for debugging
    save_tz_profile(profile.func; profile.kwargs...)
end
for profile in z_profiles
    # @info "Saving z-profile $(nameof(profile.func))" # for debugging
    save_z_profile(profile.func; profile.kwargs...)
end

nothing


text:
module SpineBasedRecordLinkage

export run_linkage, summarise_linkage_run, compare_linkage_runs

include("unexported/distances.jl")     # Independent
include("unexported/TableIndexes.jl")  # Independent
include("unexported/config.jl")        # Depends on distances
include("runlinkage.jl")               # Exports: run_linkage
include("reporting.jl")                # Exports: summarise_linkage_run, compare_linkage_runs

using .distances
using .TableIndexes
using .config
using .runlinkage
using .reporting

end


text:
module MyExample

using ForwardDiff

# Write your package code here.
include("extra_file.jl")

export my_f # from extra_file.jl
export derivative_of_my_f

end


text:
using SuperTB, PyPlot, LinearAlgebra
close("all")

# set up a model (amorphous system with AII-class hopping)
tb = TBModel(hopAgarwala)
tb.periodicity = trues(2) # periodic boundary conditions
tb.basisdim    = 2

# run over this parameter range
sites = 24^2
Nensembles = 20

# spans/preallocation
enespan = range(-5,5,length=150)
kx = range(-60*pi,60*pi,length=301); ky = kx
kspan = [(i,j) for i in kx, j in ky]
sdosvals = zeros(Float64, (size(kspan)...,size(enespan)...))

for i = 1:Nensembles
    tb.pos = rand(sites,2)
    # solve the system
    H = Hermitian(Matrix(setHopping!(tb)))
    F = eigen(H)

    # compute sdos
    sdosvals .+= sdos(F.values, F.vectors, tb, kspan, enespan)
end
sdosvals ./= Nensembles # normalize

figure()
pcolormesh(kx, enespan, sdosvals[:,151,:].',
       cmap=get_cmap("inferno"), vmin=0, vmax=maximum(sdosvals), linewidth=0)
colorbar()
#PyPlot.svg(true)
#savefig("sdos_ensembleaveraged")

max_issues_repo_name: thchr/super-tb

max_forks_repo_name: thchr/super-tb


text:
"""
    random_matrix(n, m)
Generates a correlation matrix from a `n` by `m` random matrix.
`n` is the number of culumns. `m` is the number of rows.

# Example
```julia
julia> random_matrix(10_000, 5)
5×5 Matrix{Float64}:
  1.0          0.00446344   0.00358672  -0.00662204   0.00685964
  0.00446344   1.0         -0.00526056   0.0165235    0.00962451
  0.00358672  -0.00526056   1.0          0.00227261  -0.0133972
 -0.00662204   0.0165235    0.00227261   1.0         -9.9025e-5
  0.00685964   0.00962451  -0.0133972   -9.9025e-5    1.0
```
"""
function random_matrix(n, m)
    M = randn(n, m)
    R = cor(M)
    return R
end

"""

# Example
```julia
julia> using StatsBase
julia> X = hcat(map(i -> sample([1,2,3,4,5], 100), 1:10)...)
100×10 Matrix{Int64}:
 3  1  5  1  2  5  3  3  1  3
 2  1  1  1  2  4  4  4  3  1
 5  4  3  5  2  5  3  3  3  1
 3  4  4  1  3  1  2  3  1  3
 2  3  1  4  1  2  2  1  5  3
 5  2  3  1  4  1  4  3  5  3
 2  2  4  3  4  2  4  3  1  3
 ⋮              ⋮
 2  5  3  1  1  4  1  4  5  5
 3  2  4  3  4  5  5  1  1  4
 5  1  5  5  2  1  2  1  5  1
 1  5  1  4  3  4  2  4  1  1
 1  2  1  1  2  1  1  1  2  5
 1  1  2  5  5  5  1  1  5  3

julia> random_sample(X)

```
"""
function random_sample(x)
	n = size(x, 1)
	hcat(map(i -> sample(i, n; replace = true), eachcol(x))...)
end


text:
"""
Функция сложения
"""
function summa(a::Float64,b::Float64,c::Float64)
    d = round(a+b+c)
end

function summa(a::Int64,b::Int64,c::Int64)
    d = a+b+c
end

function summa(a,b,c)
    return "Invalid input values"
end

"""
#Блабабла
Функция расчета квадратных корней
"""
function calcQuadRoots(a,b,c)
    D = b^2 - 4*a*c
    if D<0
        x1 = (-b+sqrt(Complex(D)))/2/a
        x2 = (-b-sqrt(Complex(D)))/2/a
    else
        x1 = (-b+sqrt(D))/2/a
        x2 = (-b-sqrt(D))/2/a
    end
    return x1,x2
end


text:
using Test

@testset "Test of reconstr function" begin

T    = 2 * pi

Ntaumm = 64

tau  = zeros(Float64, Ntaumm)
tau .= T * collect(0:Ntaumm-1) / Ntaumm

aa = exp.(cos.(tau)) * transpose(ones(4))

bb = reconstr(aa, 1, T, Ntaumm)

@test isapprox(vec(abs.(bb .- exp(cos(1)))), zeros(length(bb)), atol=1e-15 )

end


text:
module CommonTools

require("LoadEnvironment.jl")
require("DataTypes.jl")

using DataTypes

export tanh, tanh!, hard_tanh, ultra_fast_tanh
export square!, l2sq, l2sq_distance
export add!, subtract_vec_matrix!
export rmsPropAddGradScale!, pregrads_l2sq_distance, CVMIntermediates
export adaGradAddGradScale!
export cosine_distance_grad, pregrads_cosine_distance
export mult_by_0!, sigmoid, loopdot
export get_rand
export @nogc

using Corpora.Token

immutable CVMIntermediates
    ofull::Matrix{Float};
    result::Vector{Float};
    local1::Vector{Float};
end

CVMIntermediates(wvd::Int64, max_len::Int64) = CVMIntermediates(
                                 zeros(Float, wvd, max_len),
                                 zeros(Float, wvd),
                                 zeros(Float, wvd));

sigmoid(z::Float) = const_one/(const_one + exp(-z));

# Disable garbage collection in the block this macro is applied to
# and perform one large garbage collection at the end.
# Prevents the garbage collection from being triggered constantly
# and instead performs one big sweep.
macro nogc(ex)
    quote
        local val
        try
            gc_disable()
            val = $(esc(ex))
        finally
            gc_enable()
            gc()
        end
        val
    end
end

# Get random integer in range between 1 and l
function get_rand(l)
    return 1+int64(trunc(rand() * l))
end

# Dot product of two vectors
function loopdot(result::DenseArray{Float,1}, ip::DenseArray{Float,1})
    o = 0.0
    @simd for i=1:size(result,1)
        @inbounds o += result[i] * ip[i]
    end
    return o
end

# Multiply all entries of a matrix by zero
function mult_by_0!(x::Matrix{Float})
    @simd for i=1:length(x)
        @inbounds x[i] *= const_zero
    end
end

# Multiply all entries of a vector by zero
# TODO: write one generic function that works for all array-type objects
function mult_by_0!(x::Vector{Float})
    @simd for i=1:length(x)
        @inbounds x[i] *= const_zero
    end
end

# Apply tanh to a float
function tanh(x::Float)
    ex = exp(2*x)
    return (ex-1)/(ex+1)
end

# Apply tanh in place to all entries of a given matrix
function tanh!(x::DenseArray{Float,2})
    @simd for i=1:length(x)
        x[i] = tanh(x[i])
    end
end

# Crude approximation of tanh (inspired by Theano)
function hard_tanh(x::Float64)
    z = x
    if x < -1.0
        z = -0.995054
    elseif x > 1.0
        z = 0.995054
    end
    return z
end

# Apply approximation of tanh in place to all entries of given array
function hard_tanh!(x::DenseArray{Float})
    @simd for i=1:length(x)
        x[i] = hard_tanh(x[i])
    end
end

# Fast approximation of tanh, inspired by Theano
function ultra_fast_tanh(x::Float64)
    z = 0.0
    if x >= 0.0
        if x < 1.7
            z = (1.5 * x / (1.0 + x))
        elseif x < 3.0
            z = (0.935409070603099 + 0.0458812946797165 * (x - 1.7))
        else
            z = 0.99505475368673
        end
    else
        xx = -x
        if xx < 1.7
            z = (1.5 * xx / (1.0 + xx))
        elseif xx < 3
            z = (0.935409070603099 + 0.0458812946797165 * (xx - 1.7))
        else
            z = 0.99505475368673
        end
        z = -z
    end
    return z
end

# Apply approximation of tanh in place to all entries of given array
function ultra_fast_tanh!(x::DenseArray{Float})
    @simd for i=1:length(x)
        x[i] = ultra_fast_tanh(x[i])
    end
end

# Square all entries of given matrix in place
function square!(m::Matrix{Float})
    for y=1:size(m,2)
        @simd for x=1:size(m,1)
            m[x,y] *= m[x,y]
        end
    end
end

# Compute squared l2 norm of given vector
function l2sq(v::Vector{Float})
    result = 0.0
    @simd for i=1:size(v,1)
        result += v[i] * v[i]
    end
    return result
end

# Compute squared l2 norm of given matrix
# TODO: Write generic function that works for all array-type objects
function l2sq(m::Matrix{Float})
    result = 0.0
    @simd for i=1:length(m)
        @inbounds result += m[i] * m[i]
    end
    return result
end

# Computes cosine distance between two vectors
cosine_distance(w0,w1) = -(dot(w0,w1))/norm(w0)/norm(w1)
# Computes cosine distance gradient
function cosine_distance_grad(w0,w1)
    c = dot(w0,w1)
    s0 = norm(w0)
    s1 = norm(w1)
    sb = s0*s1
    sbsq = sb*sb

    grad_w0 = (w1 * sb - c * w0/s0 * s1) / sbsq * (-1)
    grad_w1 = (w0 * sb - c * w1/s1 * s0) / sbsq * (-1)

    return grad_w0, grad_w1
end

function pregrads_cosine_distance(v1::Vector{Float}, v2::Vector{Float},
                                vn::Vector{Float})
    # We are still below the margin
    ## precompute
    grad_v1_1, grad_v2 = cosine_distance_grad(v1, v2)
    grad_v1_2, grad_vn = cosine_distance_grad(v1, vn)

    # gradient of the contrastive criterion w.r.t.
    # the valid pair difference and the noise pair
    # difference
    return vn-v2, v2-v1, v1-vn
end

function l2sq_distance(v1::DenseArray{Float,1}, v2::DenseArray{Float,1})
    @assert length(v1) == length(v2)
    result = 0.0
    @simd for i=1:length(v1)
        @inbounds result += (v1[i]-v2[i])^2
    end
    return result
end

function pregrads_l2sq_distance(v1::DenseArray{Float,1}, v2::DenseArray{Float,1},
                                vn::DenseArray{Float,1})
    # gradient of the contrastive criterion w.r.t.
    # the valid pair difference and the noise pair
    # difference
    return vn-v2, v2-v1, v1-vn
end
# Adds two vectors v and o up and store the result in o
function add!(v::DenseArray{Float}, o::DenseArray{Float})
    """
        Inplace add function that adds up two vectors.
        :param v: vector 1
        :param o: vector 2, will be modified to contain the result
                  of the addition
    """
    @assert length(v) == length(o)
    @simd for i=1:length(v)
        @inbounds o[i] += v[i]
    end
end

# Subtract vector from given column in matrix
# (in place)
function subtract_vec_matrix!(v::DenseArray{Float,1}, idx::Int,
                                        m::DenseArray{Float,2})
    @simd for i=1:size(v,1)
        m[i,idx] -= v[i]
    end
end

end

max_stars_repo_name: ogh/binclusion

max_issues_repo_name: ogh/binclusion

max_forks_repo_name: ogh/binclusion


text:
import Base.show

# IPython message structure
type Msg
    idents::Vector{String}
    header::Dict
    content::Dict
    parent_header::Dict
    metadata::Dict
    function Msg(idents, header::Dict, content::Dict,
                 parent_header=Dict{String,Any}(), metadata=Dict{String,Any}())
        new(idents,header,content,parent_header,metadata)
    end
end

# PUB/broadcast messages use the msg_type as the ident, except for
# stream messages which use the stream name (e.g. "stdout").
# [According to minrk, "this isn't well defined, or even really part
# of the spec yet" and is in practice currently ignored since "all
# subscribers currently subscribe to all topics".]
msg_pub(m::Msg, msg_type, content, metadata=Dict{String,Any}()) =
  Msg([ msg_type == "stream" ? content["name"] : msg_type ],
      ["msg_id" => uuid4(),
       "username" => m.header["username"],
       "session" => m.header["session"],
       "msg_type" => msg_type],
      content, m.header, metadata)

msg_reply(m::Msg, msg_type, content, metadata=Dict{String,Any}()) =
  Msg(m.idents,
      ["msg_id" => uuid4(),
       "username" => m.header["username"],
       "session" => m.header["session"],
       "msg_type" => msg_type],
      content, m.header, metadata)

function show(io::IO, msg::Msg)
    print(io, "IPython Msg [ idents ")
    print_joined(io, msg.idents, ", ")
    print(io, " ] {\n  header = $(msg.header),\n  metadata = $(msg.metadata),\n  content = $(msg.content)\n}")
end

function send_ipython(socket, m::Msg)
    @vprintln("SENDING $m")
    for i in m.idents
        send(socket, i, SNDMORE)
    end
    send(socket, "<IDS|MSG>", SNDMORE)
    header = json(m.header)
    parent_header = json(m.parent_header)
    metadata = json(m.metadata)
    content = json(m.content)
    send(socket, hmac(header, parent_header, metadata, content), SNDMORE)
    send(socket, header, SNDMORE)
    send(socket, parent_header, SNDMORE)
    send(socket, metadata, SNDMORE)
    send(socket, content)
end

function recv_ipython(socket)
    msg = recv(socket)
    idents = String[]
    s = bytestring(msg)
    @vprintln("got msg part $s")
    while s != "<IDS|MSG>"
        push!(idents, s)
        msg = recv(socket)
        s = bytestring(msg)
        @vprintln("got msg part $s")
    end
    signature = bytestring(recv(socket))
    request = Dict{String,Any}()
    header = bytestring(recv(socket))
    parent_header = bytestring(recv(socket))
    metadata = bytestring(recv(socket))
    content = bytestring(recv(socket))
    if signature != hmac(header, parent_header, metadata, content)
        error("Invalid HMAC signature") # What should we do here?
    end
    m = Msg(idents, JSON.parse(header), JSON.parse(content), JSON.parse(parent_header), JSON.parse(metadata))
    @vprintln("RECEIVED $m")
    return m
end


text:
snoop_file = joinpath(@__DIR__, "check-snoop.txt")
ENV["PkgCompilerEnabled"] = "1"

using PkgCompiler
PkgCompiler.set_global_snoop_file(snoop_file)

include("code.jl")

PkgCompiler.finish_snooping()


text:
using StudentHousing, StatsBase, Distributions, JuMP, Gurobi
include("utils.jl")

srand(32)

# =============================================================================
# Toy problem data
# =============================================================================
nhouses = 50
nbedrooms_range = collect(1:3)
nbedrooms_frequency = Weights([0.5, 0.25, 0.25])
nbathrooms_range = collect(1:2)
nbathrooms_frequency = Weights([0.5, 0.5])
prices_range_pp = [800.0, 1000.0]
area_ranges = [0.0, 800.0]
market_data = StudentHousing.MarketData(nbedrooms_range, nbedrooms_frequency,
    nbathrooms_range, nbathrooms_frequency, prices_range_pp, area_ranges)
budget = 1e3
d = StudentHousingData(market_data, nhouses = nhouses, budget = budget,
                  demand_distribution = Uniform(0.0, 10.0), easy=false)
nhouses = length(d.houses)

# =============================================================================
# Solve with direct formulation
# =============================================================================
m1 = StudentHousing.ia_model(d)
@assert solve(m1, relaxation=true) == :Optimal
getobjectivevalue(m1)

# =============================================================================
# Solve problem with column generation algorithm
# =============================================================================
m, V_generated, λ_generated, house_choice = solve_ia_generation(d)

# =============================================================================
# Look at solution
# =============================================================================
println(getobjectivevalue(m))

# If not integer...
for i = 1:nhouses
    setcategory(m[:λ][i], :Bin)
end
for λnew in λ_generated
    setcategory(λnew, :Bin)
end
solve(m)

# Recover the solution

# In case one of our initial columns were in the optimal basis
results = zeros(nhouses, length(d.patterns_allow))
for k = 1:nhouses
    if getvalue(m[:λ][k]) ≈ 1.0
        results[k, StudentHousing.house_allowedby(k, d)[1]] = 1.0
    end
end
# Check all the other columns we added
for k = 1:length(V_generated)
    if getvalue(λ_generated[k]) ≈ 1.0
        # Find the assignment that we picked for our house
        patterns = find(V_generated[k] .≈ 1.0)
        # Check the house index
        h = house_choice[k]
        # Mark the assignment
        results[h, patterns] .= 1.0
    end
end

println("House 1 assigned to:", find(results[1,:] .≈ 1.0))
# House 1 assigned to:[6, 7]

# sanity check
StudentHousing.house_allowedby(1, d)
# 5-element Array{Int64,1}:
#   6
#   7
#  11
#  12
#  13

println("House 2 assigned to:", find(results[2,:] .≈ 1.0))
# House 2 assigned to:[3]
# sanity check
StudentHousing.house_allowedby(2, d)
# 11-element Array{Int64,1}:
#   3
#   4
#   5
#   6
#   7
#   8
#   9
#  10
#  11
#  12
#  13

for i = 1:nhouses
    println(println("House $i assigned to:", find(results[i,:] .≈ 1.0)))
end

# =============================================================================
# Compare with original single-stage MIP, 329 pattern example
# =============================================================================
nbedrooms_range = collect(1:2)
nbedrooms_frequency = Weights([0.75, 0.25])
nbathrooms_range = collect(1:2)
nbathrooms_frequency = Weights([0.75; 0.25])
prices_range_pp = [800.0, 1000.0]
area_ranges = [0.0, 700.0, 850.0]
budget = 1000.0
market_data = StudentHousing.MarketData(nbedrooms_range, nbedrooms_frequency,
    nbathrooms_range, nbathrooms_frequency, prices_range_pp, area_ranges)
nhouses = 50
srand(1)
d = StudentHousingData(market_data, nhouses = nhouses, budget = budget,
                  demand_distribution = Uniform(0.0, 100.0))

# =============================================================================
# Solve with direct formulation
# =============================================================================
m1 = StudentHousing.ia_model(d)
@assert solve(m1) == :Optimal
println("Shortage = ", sum(sum(d.demands)) - getobjectivevalue(m1))
# Shortage = 16616.0
@assert solve(m1, relaxation=true) == :Optimal
println("Shortage = ", sum(sum(d.demands)) - getobjectivevalue(m1))
# Shortage = 16615.18328390582
# These are the same objectives we had with the original single stage
# formulation

# =============================================================================
# Solve problem with column generation algorithm, house subproblems
# =============================================================================
tic()
m2, V_generated, λ_generated, house_choice = solve_ia_generation(d)
toc()
# elapsed time: 2.031915379 seconds

# =============================================================================
# Look at solution
# =============================================================================
println(getobjectivevalue(m2))
# 25.63353800440406
println("Shortage = ", sum(sum(d.demands)) - getobjectivevalue(m2))
# Shortage = 16615.366461995596
# a little stronger than the LP relaxation

# Not integer, so get integer solution with existing columns
r = recover_integer_soln(m2, d)
println("Integer shortage = ", sum(sum(d.demands)) - sum(sum(r)))
# Integer shortage = 16616.0
# ... which matches the MIP solution


text:
using Revise
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using Random
using Test
using SparseArrays
Random.seed!(233)

sparse_indexing = load_op_and_grad("build/libSparseIndexing", "sparse_indexing", multiple=true)
# function Base.:getindex(s::SparseTensor, i1::Union{PyObject,Array{S,1}},
#          i2::Union{PyObject,Array{T,1}}) where {S<:Real,T<:Real}
#     m_, n_ = length(i1), length(i2)
#     i1 = convert_to_tensor(i1, dtype=Int64)
#     i2 = convert_to_tensor(i2, dtype=Int64)
#     ii1, jj1, vv1 = find(s)
#     m = tf.convert_to_tensor(s.o.shape[1],dtype=tf.int64)
#     n = tf.convert_to_tensor(s.o.shape[2],dtype=tf.int64)
#     ii2, jj2, vv2 = sparse_indexing(ii1,jj1,vv1,m,n,i1,i2)
#     SparseTensor(ii2, jj2, vv2, m_, n_)
# end
################## End Load Operator ##################

i1 = unique(rand(1:20,3))
j1 = unique(rand(1:30,3))
A = sprand(20,30,0.3)
@show i1, j1
Ad = Array(A[i1, j1])
B = SparseTensor(A)
Bd = Array(B[i1, j1])
sess = Session()
init(sess)
Bd_ = run(sess, Bd)
@test Ad≈Bd_

# error()
# TODO: change your test parameter to `m`
# gradient check -- v
function scalar_function(m_)
    A = SparseTensor(ii1, jj1, m_, 20, 30)
    B = A[i1, j1]
    return sum(B)^2
end

# TODO: change `m_` and `v_` to appropriate values
m = 20
n = 30
ii1,jj1,vv1 = find(B)
m_ = constant(rand(length(B.o.values)))
v_ = rand(length(B.o.values))
y_ = scalar_function(m_)
dy_ = gradients(y_, m_)
ms_ = Array{Any}(undef, 5)
ys_ = Array{Any}(undef, 5)
s_ = Array{Any}(undef, 5)
w_ = Array{Any}(undef, 5)
gs_ =  @. 1 / 10^(1:5)

for i = 1:5
    g_ = gs_[i]
    ms_[i] = m_ + g_*v_
    ys_[i] = scalar_function(ms_[i])
    s_[i] = ys_[i] - y_
    w_[i] = s_[i] - g_*sum(v_.*dy_)
end

sess = Session()
init(sess)
sval_ = run(sess, s_)
wval_ = run(sess, w_)
close("all")
loglog(gs_, abs.(sval_), "*-", label="finite difference")
loglog(gs_, abs.(wval_), "+-", label="automatic differentiation")
loglog(gs_, gs_.^2 * 0.5*abs(wval_[1])/gs_[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
loglog(gs_, gs_ * 0.5*abs(sval_[1])/gs_[1], "--",label="\$\\mathcal{O}(\\gamma)\$")

plt.gca().invert_xaxis()
legend()
xlabel("\$\\gamma\$")
ylabel("Error")


text:
#=
 # In this section you can exclude folders and files you don't want to be index in the nav bar,
 # as well as words you want either to be casted to uppercase or not be formatted at all.
 # md files should be stated with their full name, while folders and other file types can also just
 # contain parts of names.
 # TODO: The file black does not filter very efficiently. Could be replaced by regex
=#
##########################################################################################################
file_black_list = ["404.md", "README.md", "config.md", "search.md", "impressum.md", "index.md", "internal.md", "index_content.md"]
folder_black_list = ["_", "node_modules", ".git", "LICENSE", ".idea"]

uppercase_list = ["eeg", "glmm"]
format_blacklist = ["porics"]
##########################################################################################################


text: Base.show(io::IO, self::Term) = print(io, "$(Meta.parse("($(args(self)[1])) $(data(self)) ($(args(self)[2]))"))")


text:
################################################################################
#
#   Smith normal form
#
################################################################################

@doc Markdown.doc"""
    diagonal_form(A::SMat{fmpz}) -> SMat{fmpz}

A matrix $D$ that is diagonal and obtained via unimodular row and column operations.
Like a snf without the divisibility condition.
"""
function diagonal_form(A::SMat{fmpz})
  s = 0
  while !isdiagonal(A)
    s += 1
    A = hnf(A')
  end
  if isodd(s)
    return A'
  else
    return A
  end
end

@doc Markdown.doc"""
    diagonal(A::SMat) -> fmpz[]

The diagonal elements of $A$ in an array.
"""
function diagonal(A::SMat)
  return [A[i,i] for i=1:min(nrows(A), ncols(A))]
end

@doc Markdown.doc"""
    isdiagonal(A::SMat) -> Bool

True iff only the i-th entry in the i-th row is non-zero.
"""
function isdiagonal(A::SMat)
  i = 1
  while i<= nrows(A) && length(A.rows[i].pos)>0
    if length(A.rows[i].pos) > 1
      return false
    end
    if A.rows[i].pos[1] != i
      return false
    end
    i += 1
  end
  if i > nrows(A)
    return true
  end
  return all(i -> length(A.rows[i].pos) == 0, i:nrows(A))
end

@doc Markdown.doc"""
    snf(A::SMat{fmpz})

The Smith normal form (snf) of $A$.
"""
function snf(A::SMat{fmpz})
  A = diagonal_form(A)
  e = elementary_divisors(A)
  for i=1:length(e)
    if iszero(e[i])
      return A
    end
    A.rows[i].values[1] = e[i]
  end
  return A
end

@doc Markdown.doc"""
    elementary_divisors(A::SMat{fmpz}) -> Vector{fmpz}

The elementary divisors of $A$, i.e. the diagonal elements of the Smith normal
form of $A$.
"""
function elementary_divisors(A::SMat{fmpz})
  A = diagonal_form(A)
  b = [x for x in diagonal(A) if !iszero(x)]
  c = coprime_base(b)
  e = [fmpz(1) for x = b]
  for p = c
    if isone(p)
      continue
    end
    v = [valuation(x, p) for x = b]
    sort!(v)
    for i=1:length(v)
      e[i] *= p^v[i]
    end
  end
  for i = length(e):min(nrows(A), ncols(A))
    push!(e, fmpz(0))
  end
  return e
end


text:
using SASLib

files = filter(x -> endswith(x, "sas7bdat"), Base.Filesystem.readdir())
for f in files
    println("=== $f ===")
    result = readsas(f)
end


text:
module ModelTools
using SemanticModels.Parsers

export model, callsites, structured, AbstractModel,
    pusharg!, setarg!, bodyblock, argslist, issome,
    head, isblock, isfunc, or, and, isexpr, iscall,
    isusing, isimport, funcarg, funcname,
    Edge, Edges, typegraph, @typegraph

# TODO Possible imports/exports: invoke

""" isexpr(x)

predicate for isa(x, Expr).
"""
isexpr(x) = isa(x, Expr)

"""    AbstractModel

a placeholder struct to dispatch on how to parse the expression tree into a model.
"""
abstract type AbstractModel end

""" model(::AbstractModel, expr::Expr)

dig into the expression that describes a model and break it down into components. This allows you to construct a structured representation of the modeling problem at the expression level. Just like how julia modeling frameworks build structured representations of the problems in data structures. This version builds them at the expression level.

The first argument is the type you want to construct, the second argument is the expression that you want to analyze. For example

```
model(ExpODEModel, expr)::ExpODEModel
```

"""
function model(::Type{T}, expr::Expr) where T<:AbstractModel
    error("NotImplemented: model(::$T,::Expr")
end

function invoke(m::AbstractModel, args...)
    Mod = eval(m.expr)
    Base.invokelatest(Mod.main, args...)
end

"""    callsites(expr::Expr, name::Symbol)

extract the location where the function `name` is called in `expr`.
"""
function callsites(expr::Expr, name::Symbol)
    matches = Expr[]
    f(x::Any) = x
    f(x::Expr) = begin
        if x.head == :call && x.args[1] == name
            push!(matches, x)
            return x
        else
            return walk(x, f, g)
        end
    end
    g(x) = x
    walk(expr, f, g)
    return matches
end

"""    structured(func, var::Symbol, assign=true)

extract the expressions that use structuring/destructuring assignment
to name the components of `var`
"""
function structured(func::Expr, var::Symbol, assign=true)
    body = func.args[2]
    paramvec = Expr[]
    for line in body.args
        if isa(line, LineNumberNode)
            continue
        end
        if  assign && line.head == :(=) && inexpr(line.args[2], var)
            push!(paramvec, line)
        end
        if !assign && line.head == :(=) && inexpr(line.args[1], var)
            push!(paramvec, line)
        end
    end
    return paramvec
end

lhs(x::Expr) = begin
    @show x
    x.head == :(=) || error("x is not an assignment")
    return x.args[1]
end

"""    bodyblock(expr::Expr)

get the array of args representing the body of a function definition.
"""
function bodyblock(expr::Expr)
    expr.head == :function || error("$expr is not a function definition")
    return expr.args[2].args
end

"""    funclines(expr::Expr, s::Symbol)

clean up the lines of a function definition for presentation
"""
function funclines(expr::Expr, s::Symbol)
    q = Expr(:block)
    q.args = (filter(isexpr, findfunc(expr, s))[end]
              |> bodyblock
              |> arr -> filter(x->!isa(x, LineNumberNode),arr))
    return q
end

"""    argslist(expr::Expr)

get the array of args representing the arguments of a defined function.
the first element of this list is the function name

See also [`bodyblock`](@ref), [`pusharg!`](@ref),
"""
function argslist(expr::Expr)
    expr.head == :function || error("$expr is not a function definition")
    return expr.args[1].args
end

"""    pusharg!(expr::Expr, s::Symbol)

push a new argument onto the definition of a function.

See also [`argslist`](@ref), [`setarg!`](@ref)
"""
function pusharg!(ex::Expr, s::Symbol)
    ex.head == :function || error("ex is not a function definition")
    push!(argslist(ex), s)
    return ex
end

"""    setarg!(expr::Expr, s::Symbol)

replace the argument in a function call.

See also [`argslist`](@ref), [`pusharg!`](@ref)
"""
function setarg!(ex::Expr, old, new)
    ex.head == :call || error("ex is not a function call")
    for (i, x) in enumerate(ex.args)
        if x == old
            ex.args[i] = new
        end
    end
    return ex
end

function funcarg(ex::Expr)
    return ex.args[1].args[2]
end

include("exprs.jl")
include("typegraphs.jl")
include("Transformations.jl")
include("SimpleModels.jl")
include("ExpODEModels.jl")
include("ExpStateModels.jl")
include("MonomialRegressionModels.jl")
include("CategoryTheory.jl")
include("PetriModels.jl")
include("OpenModels.jl")
include("OpenPetris.jl")
include("WiringDiagrams.jl")

end


text:
abstract type Cofils <: Persa.CFModel end

mutable struct GlobalCofils <: Cofils
    slmodel::Nullable{DecisionTree.Ensemble}
    lv::LatentVariable
    features::Int
    preferences::Persa.RatingPreferences
end

mutable struct UserCofils <: Cofils
    slmodel::Nullable{DecisionTree.Ensemble}
    lv::LatentVariable
    features::Int
    means::Array{Float64}
    preferences::Persa.RatingPreferences
end

mutable struct ItemCofils <: Cofils
    slmodel::Nullable{DecisionTree.Ensemble}
    lv::LatentVariable
    features::Int
    means::Array{Float64}
    preferences::Persa.RatingPreferences
end

userlv(model::Cofils, id::Int) = user(model.lv, id)
itemlv(model::Cofils, id::Int) = item(model.lv, id)

function Cofils(dataset::Persa.CFDatasetAbstract, features::Int; normalization::Symbol = :user, extractlv::Function = svd)
    @assert in(normalization, [:user, :item, :none]) "Incorrect value of normalization. Use :user, :item or :none."

    if normalization == :user
        means = Persa.means(dataset; mode = :user, α = 10)
        return UserCofils(Nullable{DecisionTree.Ensemble}(), extractlv(dataset, features), features, means, dataset.preferences)
    elseif normalization == :item
        means = Persa.means(dataset; mode = :item, α = 10)
        return ItemCofils(Nullable{DecisionTree.Ensemble}(), extractlv(dataset, features), features, means, dataset.preferences)
    end

    return GlobalCofils(Nullable{DecisionTree.Ensemble}(), extractlv(dataset, features), features, dataset.preferences)
end

function Persa.train!(model::Cofils, dataset::Persa.CFDatasetAbstract; nfeaturestrees::Int = 10, trees::Int = 20)
    if length(model.lv) < nfeaturestrees
        warn("RandomFlorest: # features trees is greeter than the features quantity.");
        nfeaturestrees = length(model.lv);
    end

    (attributes, labels) = convert2sl(model, dataset)

    model.slmodel = Nullable{DecisionTree.Ensemble}(build_forest(labels, attributes, nfeaturestrees, trees))

    return nothing
end

function convert2sl(model::GlobalCofils, dataset::Persa.CFDatasetAbstract)
    m, n = size(model.lv)

    attributes = Array{Float64,2}(length(dataset), m + n)
    labels = Array{Float64,1}(length(dataset))

    for i=1:length(dataset)
      (u, v, r) = dataset[i]

      attributes[i, :] = convert2sl(model, u, v)
      labels[i] = r
    end

    return attributes, labels
end

function convert2sl(model::UserCofils, dataset::Persa.CFDatasetAbstract)
    m, n = size(model.lv)

    attributes = Array{Float64,2}(length(dataset), m + n)
    labels = Array{Float64,1}(length(dataset))

    for i=1:length(dataset)
      (u, v, r) = dataset[i]

      attributes[i, :] = convert2sl(model, u, v)
      labels[i] = r - model.means[u]
    end

    return attributes, labels
end

function convert2sl(model::ItemCofils, dataset::Persa.CFDatasetAbstract)
    m, n = size(model.lv)

    attributes = Array{Float64,2}(length(dataset), m + n)
    labels = Array{Float64,1}(length(dataset))

    for i=1:length(dataset)
      (u, v, r) = dataset[i]

      attributes[i, :] = convert2sl(model, u, v)
      labels[i] = r - model.means[v]
    end

    return attributes, labels
end

function convert2sl(model::Cofils, user::Int, item::Int)
    m, n = size(model.lv)

    attributes = Array{Float64,2}(1, m + n)

    attributes[1,1:m] = userlv(model, user)
    attributes[1,(m+1):(m+n)] = itemlv(model, item)

    return attributes
end

function Persa.predict(model::GlobalCofils, user::Int64, item::Int64)
    attributes = convert2sl(model, user, item)

    p = apply_forest(get(model.slmodel), attributes)[1]

    return Persa.correct(p, model.preferences)
end

function Persa.predict(model::UserCofils, user::Int64, item::Int64)
    attributes = convert2sl(model, user, item)

    p = apply_forest(get(model.slmodel), attributes)[1] + model.means[user]

    return Persa.correct(p, model.preferences)
end

function Persa.predict(model::ItemCofils, user::Int64, item::Int64)
    attributes = convert2sl(model, user, item)

    p = apply_forest(get(model.slmodel), attributes)[1] + model.means[item]

    return Persa.correct(p, model.preferences)
end

Persa.canpredict(model::Cofils, user::Int, item::Int) = true


text:
using PCHIPInterpolation
using Test

function test_interpolation_is_piecewise_monotone(xs, ys, N=10000)
    itp = Interpolator(xs, ys)
    for (i,j) in monotone_intervals(ys)
        @test is_monotone([itp(x) for x ∈ range(xs[i], stop=xs[j], length=N)])
    end
end

function monotone_intervals(xs)
    N = length(xs)
    function _monotone_intervals(ch :: Channel)
        up_down = 0
        previous_i = 1
        for i = 1:N-1
            current_x, next_x = xs[i:i+1]
            if up_down != 0  &&  up_down*current_x > up_down*next_x
                put!(ch, (previous_i, i))
                previous_i = i
                up_down = 0
            end
            if up_down == 0  &&  current_x != next_x
                up_down = current_x < next_x ? 1 : -1
            end
        end
        put!(ch, (previous_i, N))
    end

    [interval for interval in Channel(_monotone_intervals)]
end

function is_monotone(zs)
    (all(is_approx_le(z, znext) for (z,znext) = zip(zs[1:end-1], zs[2:end]))
    || all(is_approx_le(znext, z) for (z,znext) = zip(zs[1:end-1], zs[2:end])))
end

function is_approx_le(a, b)
    a < b || a ≈ b
end

xs = [0.0, 1.0]
ys = [1.0, 1.0]
test_interpolation_is_piecewise_monotone(xs, ys)

xs = [0.0, 1.0]
ys = [1.0, 0.0]
test_interpolation_is_piecewise_monotone(xs, ys)

xs = [0.0, 1.0]
ys = [0.0, 1.0]
test_interpolation_is_piecewise_monotone(xs, ys)

xs = [0.0, 1.0, 2.0, 100.0]
ys = [0.0, 0.0, 1.0,   1.0]
test_interpolation_is_piecewise_monotone(xs, ys)

xs = [1.0, 2.0, 5.0, 6.0, 10.0, 22.0]
ys = [0.0, 1.0, 0.0, 1.0,  0.0,  1.0]
test_interpolation_is_piecewise_monotone(xs, ys)

xs = [   0.0, 10.0, 10000.0]
ys = [-100.0,  1.2,     1.1]
test_interpolation_is_piecewise_monotone(xs, ys)

xs = range(0, stop=2π, length=5000)
ys = sin.(xs)
test_interpolation_is_piecewise_monotone(collect(xs), ys)

# Example data from Fritsch and Carlson, SIAM J. NUMER. ANAL. 17 (1980) 238-246.
xs = [ 0.0,  2,  3,  5,  6,  8,  9,   11, 12, 14, 15]
ys = [10.0, 10, 10, 10, 10, 10, 10.5, 15, 50, 60, 85]
test_interpolation_is_piecewise_monotone(xs, ys)

# Test the test helper functions
@test is_monotone([0.0, 0.0, 0.0])
@test is_monotone([0.0, 1.0, 10.0])
@test is_monotone([10.0, 1.0, 0.0])
@test ! is_monotone([0.0, 1.0, 0.0])
@test ! is_monotone([1.0, 0.0, 1.0])

# Test internal functions

# Make sure the correct interval is identified
p = Interpolator([1.0, 2.0, 3.0], [0.0, 0.0, 0.0])
for search ∈ (PCHIPInterpolation._pchip_index_linear_search, PCHIPInterpolation._pchip_index_bisectional_search)
    @test search(p, 1.0) == 1
    @test search(p, 1.0 + eps(1.0)) == 1
    @test search(p, 2.0 - eps(2.0)) == 1
    @test search(p, 2.0) == 2
    @test search(p, 2.0 + eps(2.0)) == 2
    @test search(p, 3.0 - eps(3.0)) == 2
    @test search(p, 3.0) == 3
end

# Test integration
@test integrate(p, 1, 3) == 0
@test integrate(p, 3, 1) == 0
@test integrate(p, 1, 2) == 0
@test integrate(p, 2, 3) == 0
@test integrate(p, 1.25, 2) == 0
@test integrate(p, 1, 2.25) == 0
@test integrate(p, 1.25, 2.25) == 0

p = Interpolator([1.0, 2.0, 3.0, 4.0], [4.0, 3.0, 2.0, 1.0])
@test integrate(p, 1, 4) == 3*3/2 + 3
@test integrate(p, 1, 2) + integrate(p, 2, 4) == integrate(p, 1, 4)
@test integrate(p, 1, 2.75) + integrate(p, 2.75, 4) == integrate(p, 1, 4)
@test integrate(p, 3, 1) == -integrate(p, 1, 3)

# Test out of domain
p = Interpolator([1.0, 2.0, 3.0, 4.0], [4.0, 3.0, 2.0, 1.0])
@test p(1) == 4
@test p(4) == 1
@test_throws DomainError p(1 - 1e-6)
@test_throws DomainError p(4 + 1e-6)
@test_throws DomainError integrate(p, 1 - 1e-6, 4 + 1e-6)
@test_throws DomainError integrate(p, 1 - 1e-6, 3)
@test_throws DomainError integrate(p, 2, 4 + 1e-6)

# Test too short
@test_throws ArgumentError Interpolator([1.0], [4.0])
@test_throws ArgumentError Interpolator([], [])

# Test mismatched lengths
@test_throws DimensionMismatch Interpolator([1.0, 2.0, 3.0, 4.0], [4.0, 3.0, 2.0])
@test_throws DimensionMismatch Interpolator([1.0, 2.0, 3.0], [4.0, 3.0, 2.0, 1.0])

# Test non-increasing xs
@test_throws ArgumentError Interpolator([1.0, 2.0, 3.0, 3.0], [4.0, 3.0, 2.0, 1.0])


text:
export mkchans

function mkchans(chans::ChanSpec, S::GphysData; f::Symbol=:x, keepempty::Bool=false, keepirr::Bool=true)
  chan_list = (if chans == Int64[]
    Int64.(collect(1:S.n))
  elseif typeof(chans) <: UnitRange
    Int64.(collect(chans))
  elseif typeof(chans) <: Integer
    [Int64(chans)]
  else
    chans
  end)

  # added 2020-02-12; prevents processing an empty channel
  k = trues(length(chan_list))
  if keepempty == false
    F = getfield(S, f)
    for (j,i) in enumerate(chan_list)
      if isempty(F[i])
        k[j] = false
      end
    end
  end

  # added 2020-03-08: option to auto-delete irregular channels
  if keepirr == false
    for (j,i) in enumerate(chan_list)
      if S.fs[i] == 0.0
        k[j] = false
      end
    end
  end
  return chan_list[k]
end


text:
# signal.jl
##########################
fs=Int(1000);
Ts=1/fs;
N=Int(12000);
T=N/fs;
#t=(0:Int(N-1))*Ts;

x=zeros(Float64, N)
for m=1:N
    x[m]=sin(2*Ts*(m-1)*(Ts*(m-1)-3)*(Ts*(m-1)-6)*(Ts*(m-1)-9)*(Ts*(m-1)-12))
end

using Plots

plot(x,
    xaxis="n = 0, 1, 2, ..., N-1",
    yaxis="amplitude x[n]",
    title="Discrete time signal x(t)=sin(2t(t-3)(t-6)(t-9)(t-12))",
    grid=false,
    legend=false)
plot!(xticks = ([0 5000 10000;], [0, 5000, 10000]))
plot!(yticks = ([-1 1;], [-1, 1]))

plot(x,
    xaxis="time t [sec]",
    yaxis="amplitude",
    title="Audio Signal x(t)=sin(2t(t-3)(t-6)(t-9)(t-12))",
    grid=false,
    legend=false)
plot!(xticks = ([0 5000 10000;], [0, 5, 10]))
plot!(yticks = ([-1 1;], [-1, 1]))
savefig("signal.png")


text:
# The crabs don't seem interested in your proposed solution. Perhaps you
# misunderstand crab engineering?
#
# As it turns out, crab submarine engines don't burn fuel at a constant rate.
# Instead, each change of 1 step in horizontal position costs 1 more unit of
# fuel than the last: the first step costs 1, the second step costs 2, the
# third step costs 3, and so on.
#
# As each crab moves, moving further becomes more expensive. This changes the
# best horizontal position to align them all on; in the example above, this
# becomes 5:
#
# Move from 16 to 5: 66 fuel
# Move from 1 to 5: 10 fuel
# Move from 2 to 5: 6 fuel
# Move from 0 to 5: 15 fuel
# Move from 4 to 5: 1 fuel
# Move from 2 to 5: 6 fuel
# Move from 7 to 5: 3 fuel
# Move from 1 to 5: 10 fuel
# Move from 2 to 5: 6 fuel
# Move from 14 to 5: 45 fuel
#
# This costs a total of 168 fuel. This is the new cheapest possible outcome;
# the old alignment position (2) now costs 206 fuel instead.
#
# Determine the horizontal position that the crabs can align to using the least
# fuel possible so they can make you an escape route! How much fuel must they
# spend to align to that position?

# read input and sort it
positions = sort!(parse.(UInt16, split(readline("input.txt"), ',')))

# Similar to part 1 except that as the position moves right, the fuel cost to
# the left increases by the sum of the difference between the current position
# and the crab's starting position. The fuel to the right decreases by the sum
# of the difference between the crab's previous position and the previous
# "current" position.
fuel_left = 0
fuel_right = mapreduce(p -> (p * (p+1)) >> 1, +, positions)
right_idx = findfirst(>(0), positions)
cheapest = fuel_right
for position in 1:positions[end]
  # calculate new fuels
  left = @view positions[1:(right_idx - 1)]
  right = @view positions[right_idx:end]
  global (fuel_left, fuel_right) = (fuel_left + sum(position .- left), fuel_right - sum(right .- (position - 1)))
  if fuel_left + fuel_right < cheapest
    global cheapest = fuel_left + fuel_right
  end

  # potentially update right_idx
  if position == positions[right_idx]
    global right_idx = findnext(>(position), positions, right_idx)
    if isnothing(right_idx)
      # honestly, when this happens, we can break - continuing to move to the
      # right will only increase fuel cost because all crabs are to the left
      break
    end
  end
end

println("Cheapest fuel cost: ", cheapest)


text:
# DESCRIPTION
# Trim initialization.
#
# ------------------------------------------------------------------------------

# initial guess for state in trim solution
x0=zeros(12,1)
x0[1]=VXTRIM
x0[2]=VYTRIM
x0[3]=VZTRIM
if (VXTRIM>=60*1.688)
    if PSIDTRIM==0
        # Trim heading for zero bank angle above 60 kts
        x0[7]=0.0;
        TRIMVARS[7:8]=[8; 9]
    else
        x0[9]=atan2(VYTRIM,VXTRIM)
        TRIMVARS[7:8]=[7; 8]
        # initial guess for phi (using turn coord)
        x0[7]=atan(PSIDTRIM*sqrt(VXTRIM^2+VYTRIM^2+VZTRIM^2)/G)
        # initial guess for r
        x0[6]=G/sqrt(VXTRIM^2+VYTRIM^2+VZTRIM^2)*sin(x0[7])
        # initial guess for q
        x0[5]=x0[6]*tan(x0[7])
    end
else
    # otherwise set heading to flight path
    x0[9]=atan2(VYTRIM,VXTRIM)
    TRIMVARS[7:8]=[7; 8]
end
# initial guess for controls [deg]
u0=[50.0; 50.0; 50.0; 50.0]
# target values for state derivatives
xdot_targ=zeros(12)
xdot_targ[9]=PSIDTRIM
xdot_targ[10]=VXTRIM
xdot_targ[11]=VYTRIM
xdot_targ[12]=VZTRIM

max_stars_repo_name: umbertosaetti/J-SimpleHel

max_issues_repo_name: umbertosaetti/J-SimpleHel

max_forks_repo_name: umbertosaetti/J-SimpleHel


text:
module EAUtils

export code_escapes, @code_escapes, __clear_cache!

const CC = Core.Compiler
const EA = CC.EscapeAnalysis

# entries
# -------

import Base: unwrap_unionall, rewrap_unionall
import InteractiveUtils: gen_call_with_extracted_types_and_kwargs

"""
    @code_escapes [options...] f(args...)

Evaluates the arguments to the function call, determines its types, and then calls
[`code_escapes`](@ref) on the resulting expression.
As with `@code_typed` and its family, any of `code_escapes` keyword arguments can be given
as the optional arguments like `@code_escapes optimize=false myfunc(myargs...)`.
"""
macro code_escapes(ex0...)
    return gen_call_with_extracted_types_and_kwargs(__module__, :code_escapes, ex0)
end

"""
    code_escapes(f, argtypes=Tuple{}; [debuginfo::Symbol = :none], [optimize::Bool = true]) -> result::EscapeResult

Runs the escape analysis on optimized IR of a generic function call with the given type signature.

# Keyword Arguments

- `optimize::Bool = true`:
  if `true` returns escape information of post-inlining IR (used for local optimization),
  otherwise returns escape information of pre-inlining IR (used for interprocedural escape information generation)
- `debuginfo::Symbol = :none`:
  controls the amount of code metadata present in the output, possible options are `:none` or `:source`.
"""
function code_escapes(@nospecialize(f), @nospecialize(types=Base.default_tt(f));
                      world::UInt = get_world_counter(),
                      interp::Core.Compiler.AbstractInterpreter = Core.Compiler.NativeInterpreter(world),
                      debuginfo::Symbol = :none,
                      optimize::Bool = true)
    ft = Core.Typeof(f)
    if isa(types, Type)
        u = unwrap_unionall(types)
        tt = rewrap_unionall(Tuple{ft, u.parameters...}, types)
    else
        tt = Tuple{ft, types...}
    end
    interp = EscapeAnalyzer(interp, tt, optimize)
    results = Base.code_typed_by_type(tt; optimize=true, world, interp)
    isone(length(results)) || throw(ArgumentError("`code_escapes` only supports single analysis result"))
    return EscapeResult(interp.ir, interp.state, interp.linfo, debuginfo===:source)
end

# in order to run a whole analysis from ground zero (e.g. for benchmarking, etc.)
__clear_cache!() = empty!(GLOBAL_CODE_CACHE)

# AbstractInterpreter
# -------------------

# imports
import .CC:
    AbstractInterpreter, NativeInterpreter, WorldView, WorldRange,
    InferenceParams, OptimizationParams, get_world_counter, get_inference_cache, code_cache,
    lock_mi_inference, unlock_mi_inference, add_remark!,
    may_optimize, may_compress, may_discard_trees, verbose_stmt_info
# usings
import Core:
    CodeInstance, MethodInstance, CodeInfo
import .CC:
    InferenceResult, OptimizationState, IRCode, copy as cccopy,
    @timeit, convert_to_ircode, slot2reg, compact!, ssa_inlining_pass!, sroa_pass!,
    adce_pass!, type_lift_pass!, JLOptions, verify_ir, verify_linetable
import .EA: analyze_escapes, ArgEscapeCache, EscapeInfo, EscapeState, is_ipo_profitable

# when working outside of Core.Compiler,
# cache entire escape state for later inspection and debugging
struct EscapeCache
    cache::ArgEscapeCache
    state::EscapeState # preserved just for debugging purpose
    ir::IRCode         # preserved just for debugging purpose
end

mutable struct EscapeAnalyzer{State} <: AbstractInterpreter
    native::NativeInterpreter
    cache::IdDict{InferenceResult,EscapeCache}
    entry_tt
    optimize::Bool
    ir::IRCode
    state::State
    linfo::MethodInstance
    EscapeAnalyzer(native::NativeInterpreter, @nospecialize(tt), optimize::Bool) =
        new{EscapeState}(native, IdDict{InferenceResult,EscapeCache}(), tt, optimize)
end

CC.InferenceParams(interp::EscapeAnalyzer)    = InferenceParams(interp.native)
CC.OptimizationParams(interp::EscapeAnalyzer) = OptimizationParams(interp.native)
CC.get_world_counter(interp::EscapeAnalyzer)  = get_world_counter(interp.native)

CC.lock_mi_inference(::EscapeAnalyzer,   ::MethodInstance) = nothing
CC.unlock_mi_inference(::EscapeAnalyzer, ::MethodInstance) = nothing

CC.add_remark!(interp::EscapeAnalyzer, sv, s) = add_remark!(interp.native, sv, s)

CC.may_optimize(interp::EscapeAnalyzer)      = may_optimize(interp.native)
CC.may_compress(interp::EscapeAnalyzer)      = may_compress(interp.native)
CC.may_discard_trees(interp::EscapeAnalyzer) = may_discard_trees(interp.native)
CC.verbose_stmt_info(interp::EscapeAnalyzer) = verbose_stmt_info(interp.native)

CC.get_inference_cache(interp::EscapeAnalyzer) = get_inference_cache(interp.native)

const GLOBAL_CODE_CACHE = IdDict{MethodInstance,CodeInstance}()

function CC.code_cache(interp::EscapeAnalyzer)
    worlds = WorldRange(get_world_counter(interp))
    return WorldView(GlobalCache(), worlds)
end

struct GlobalCache end

CC.haskey(wvc::WorldView{GlobalCache}, mi::MethodInstance) = haskey(GLOBAL_CODE_CACHE, mi)

CC.get(wvc::WorldView{GlobalCache}, mi::MethodInstance, default) = get(GLOBAL_CODE_CACHE, mi, default)

CC.getindex(wvc::WorldView{GlobalCache}, mi::MethodInstance) = getindex(GLOBAL_CODE_CACHE, mi)

function CC.setindex!(wvc::WorldView{GlobalCache}, ci::CodeInstance, mi::MethodInstance)
    GLOBAL_CODE_CACHE[mi] = ci
    add_callback!(mi) # register the callback on invalidation
    return nothing
end

function add_callback!(linfo)
    if !isdefined(linfo, :callbacks)
        linfo.callbacks = Any[invalidate_cache!]
    else
        if !any(@nospecialize(cb)->cb===invalidate_cache!, linfo.callbacks)
            push!(linfo.callbacks, invalidate_cache!)
        end
    end
    return nothing
end

function invalidate_cache!(replaced, max_world, depth = 0)
    delete!(GLOBAL_CODE_CACHE, replaced)

    if isdefined(replaced, :backedges)
        for mi in replaced.backedges
            mi = mi::MethodInstance
            if !haskey(GLOBAL_CODE_CACHE, mi)
                continue # otherwise fall into infinite loop
            end
            invalidate_cache!(mi, max_world, depth+1)
        end
    end
    return nothing
end

function CC.optimize(interp::EscapeAnalyzer,
    opt::OptimizationState, params::OptimizationParams, caller::InferenceResult)
    ir = run_passes_with_ea(interp, opt.src, opt, caller)
    return CC.finish(interp, opt, params, ir, caller)
end

function CC.cache_result!(interp::EscapeAnalyzer, caller::InferenceResult)
    if haskey(interp.cache, caller)
        GLOBAL_ESCAPE_CACHE[caller.linfo] = interp.cache[caller]
    end
    return Base.@invoke CC.cache_result!(interp::AbstractInterpreter, caller::InferenceResult)
end

const GLOBAL_ESCAPE_CACHE = IdDict{MethodInstance,EscapeCache}()

"""
    cache_escapes!(caller::InferenceResult, estate::EscapeState, cacheir::IRCode)

Transforms escape information of call arguments of `caller`,
and then caches it into a global cache for later interprocedural propagation.
"""
function cache_escapes!(interp::EscapeAnalyzer,
    caller::InferenceResult, estate::EscapeState, cacheir::IRCode)
    cache = ArgEscapeCache(estate)
    ecache = EscapeCache(cache, estate, cacheir)
    interp.cache[caller] = ecache
    return cache
end

function get_escape_cache(interp::EscapeAnalyzer)
    return function (linfo::Union{InferenceResult,MethodInstance})
        if isa(linfo, InferenceResult)
            ecache = get(interp.cache, linfo, nothing)
        else
            ecache = get(GLOBAL_ESCAPE_CACHE, linfo, nothing)
        end
        return ecache !== nothing ? ecache.cache : nothing
    end
end

function run_passes_with_ea(interp::EscapeAnalyzer, ci::CodeInfo, sv::OptimizationState,
    caller::InferenceResult)
    @timeit "convert"   ir = convert_to_ircode(ci, sv)
    @timeit "slot2reg"  ir = slot2reg(ir, ci, sv)
    # TODO: Domsorting can produce an updated domtree - no need to recompute here
    @timeit "compact 1" ir = compact!(ir)
    nargs = let def = sv.linfo.def; isa(def, Method) ? Int(def.nargs) : 0; end
    local state
    if is_ipo_profitable(ir, nargs) || caller.linfo.specTypes === interp.entry_tt
        try
            @timeit "[IPO EA]" begin
                state = analyze_escapes(ir, nargs, false, get_escape_cache(interp))
                cache_escapes!(interp, caller, state, cccopy(ir))
            end
        catch err
            @error "error happened within [IPO EA], insepct `Main.ir` and `Main.nargs`"
            @eval Main (ir = $ir; nargs = $nargs)
            rethrow(err)
        end
    end
    if caller.linfo.specTypes === interp.entry_tt && !interp.optimize
        # return back the result
        interp.ir = cccopy(ir)
        interp.state = state
        interp.linfo = sv.linfo
    end
    @timeit "Inlining"  ir = ssa_inlining_pass!(ir, ir.linetable, sv.inlining, ci.propagate_inbounds)
    # @timeit "verify 2" verify_ir(ir)
    @timeit "compact 2" ir = compact!(ir)
    if caller.linfo.specTypes === interp.entry_tt && interp.optimize
        try
            @timeit "[Local EA]" state = analyze_escapes(ir, nargs, true, get_escape_cache(interp))
        catch err
            @error "error happened within [Local EA], insepct `Main.ir` and `Main.nargs`"
            @eval Main (ir = $ir; nargs = $nargs)
            rethrow(err)
        end
        # return back the result
        interp.ir = cccopy(ir)
        interp.state = state
        interp.linfo = sv.linfo
    end
    @timeit "SROA"      ir = sroa_pass!(ir)
    @timeit "ADCE"      ir = adce_pass!(ir)
    @timeit "type lift" ir = type_lift_pass!(ir)
    @timeit "compact 3" ir = compact!(ir)
    if JLOptions().debug_level == 2
        @timeit "verify 3" (verify_ir(ir); verify_linetable(ir.linetable))
    end
    return ir
end

# printing
# --------

import Core: Argument, SSAValue
import .CC: widenconst, singleton_type

Base.getindex(estate::EscapeState, @nospecialize(x)) = CC.getindex(estate, x)

function get_name_color(x::EscapeInfo, symbol::Bool = false)
    getname(x) = string(nameof(x))
    if x === EA.⊥
        name, color = (getname(EA.NotAnalyzed), "◌"), :plain
    elseif EA.has_no_escape(EA.ignore_argescape(x))
        if EA.has_arg_escape(x)
            name, color = (getname(EA.ArgEscape), "✓"), :cyan
        else
            name, color = (getname(EA.NoEscape), "✓"), :green
        end
    elseif EA.has_all_escape(x)
        name, color = (getname(EA.AllEscape), "X"), :red
    elseif EA.has_return_escape(x)
        name = (getname(EA.ReturnEscape), "↑")
        color = EA.has_thrown_escape(x) ? :yellow : :blue
    else
        name = (nothing, "*")
        color = EA.has_thrown_escape(x) ? :yellow : :bold
    end
    name = symbol ? last(name) : first(name)
    if name !== nothing && !isa(x.AliasInfo, Bool)
        name = string(name, "′")
    end
    return name, color
end

# pcs = sprint(show, collect(x.EscapeSites); context=:limit=>true)
function Base.show(io::IO, x::EscapeInfo)
    name, color = get_name_color(x)
    if isnothing(name)
        Base.@invoke show(io::IO, x::Any)
    else
        printstyled(io, name; color)
    end
end
function Base.show(io::IO, ::MIME"application/prs.juno.inline", x::EscapeInfo)
    name, color = get_name_color(x)
    if isnothing(name)
        return x # use fancy tree-view
    else
        printstyled(io, name; color)
    end
end

struct EscapeResult
    ir::IRCode
    state::EscapeState
    linfo::Union{Nothing,MethodInstance}
    source::Bool
    function EscapeResult(ir::IRCode, state::EscapeState,
        linfo::Union{Nothing,MethodInstance} = nothing,
        source::Bool=false)
        return new(ir, state, linfo, source)
    end
end
Base.show(io::IO, result::EscapeResult) = print_with_info(io, result)
@eval Base.iterate(res::EscapeResult, state=1) =
    return state > $(fieldcount(EscapeResult)) ? nothing : (getfield(res, state), state+1)

Base.show(io::IO, cached::EscapeCache) = show(io, EscapeResult(cached.ir, cached.state, nothing))

# adapted from https://github.com/JuliaDebug/LoweredCodeUtils.jl/blob/4612349432447e868cf9285f647108f43bd0a11c/src/codeedges.jl#L881-L897
function print_with_info(io::IO, (; ir, state, linfo, source)::EscapeResult)
    # print escape information on SSA values
    function preprint(io::IO)
        ft = ir.argtypes[1]
        f = singleton_type(ft)
        if f === nothing
            f = widenconst(ft)
        end
        print(io, f, '(')
        for i in 1:state.nargs
            arg = state[Argument(i)]
            i == 1 && continue
            c, color = get_name_color(arg, true)
            printstyled(io, c, ' ', '_', i, "::", ir.argtypes[i]; color)
            i ≠ state.nargs && print(io, ", ")
        end
        print(io, ')')
        if !isnothing(linfo)
            def = linfo.def
            printstyled(io, " in ", (isa(def, Module) ? (def,) : (def.module, " at ", def.file, ':', def.line))...; color=:bold)
        end
        println(io)
    end

    # print escape information on SSA values
    # nd = ndigits(length(ssavalues))
    function preprint(io::IO, idx::Int)
        c, color = get_name_color(state[SSAValue(idx)], true)
        # printstyled(io, lpad(idx, nd), ' ', c, ' '; color)
        printstyled(io, rpad(c, 2), ' '; color)
    end

    print_with_info(preprint, (args...)->nothing, io, ir, source)
end

function print_with_info(preprint, postprint, io::IO, ir::IRCode, source::Bool)
    io = IOContext(io, :displaysize=>displaysize(io))
    used = Base.IRShow.stmts_used(io, ir)
    if source
        line_info_preprinter = function (io::IO, indent::String, idx::Int)
            r = Base.IRShow.inline_linfo_printer(ir)(io, indent, idx)
            idx ≠ 0 && preprint(io, idx)
            return r
        end
    else
        line_info_preprinter = Base.IRShow.lineinfo_disabled
    end
    line_info_postprinter = Base.IRShow.default_expr_type_printer
    preprint(io)
    bb_idx_prev = bb_idx = 1
    for idx = 1:length(ir.stmts)
        preprint(io, idx)
        bb_idx = Base.IRShow.show_ir_stmt(io, ir, idx, line_info_preprinter, line_info_postprinter, used, ir.cfg, bb_idx)
        postprint(io, idx, bb_idx != bb_idx_prev)
        bb_idx_prev = bb_idx
    end
    max_bb_idx_size = ndigits(length(ir.cfg.blocks))
    line_info_preprinter(io, " "^(max_bb_idx_size + 2), 0)
    postprint(io)
    return nothing
end

end # module EAUtils

max_stars_repo_name: JuliaLabs/julia

max_issues_repo_name: JuliaLabs/julia

max_forks_repo_name: JuliaLabs/julia


text:
"""
    module Domains

for constructing domains to be fed to Simulation
"""
module Domains

export Domain

# using ...Defaults
using ..Boundaries
using ..DielectricFunctions
using ..Dispersions
using ..Lattices
using ..Points
using ..Shapes
# using LinearAlgebra
# using RecipesBase

"""
    struct Domain

    Domain(type,boundary,lattice,[dielectric,pump];align=false,[name]) -> domain
    Domain(type,lattice,boundary,[dielectric,pump];align=false,[name]) -> domain

combines `boundary` with `lattice` (in either order) to generate a list of sites.
`T` just labels the kinds of domain (e.g. Cavity, or Waveguide)
each site is labeled by being in the interior, or containing the interior (`domain.surface`),
or being a corner.

----------------------
    (::Domain)(args...;align=false) -> dom

construct a new domain with modified parameters in `args`
This is for "updating" non-geometric fields of the immutable `domain`.
For example, to change the polarity of the boundary layers, do
    `new_dom = old_dom(conj(old_dom.boundary))`

For geometric parameters it recomputes the whole domain, so this command is no
more efficient than explicitly constructing a new domain. It might save some typing, however.
"""
struct Domain{N,TBND,TLAT,TDF,TPF,TCHI}
    type::Symbol
    name::Symbol
    boundary::TBND
    lattice::TLAT
    dielectric::TDF
    pump::TPF

    x::Array{Point{N},1}
    indices::Array{CartesianIndex{N},1}
    ε::Array{ComplexF64,1}
    F::Array{Float64,1}
    χ::TCHI

    interior::BitArray{1}
    bulk::BitArray{1}
    surface::BitArray{1}
    corner::BitArray{1}

    nnm::NTuple{N,Array{Int,1}}
    nnp::NTuple{N,Array{Int,1}}

    function Domain(
        type::Symbol,
        name::Symbol,
        boundary::TBND,
        lattice::TLAT,
        dielectric::TDF,
        pump::TPF,
        x::Array{Point{N},1},
        indices::Array{CartesianIndex{N},1},
        ε::Array{ComplexF64,1},
        F::Array{Float64,1},
        χ::TC,
        interior::BitArray{1},
        bulk::BitArray{1},
        surface::BitArray{1},
        corner::BitArray{1},
        nnm::NTuple{N,Array{Int,1}},
        nnp::NTuple{N,Array{Int,1}}
        ) where {TBND<:Boundary,TLAT<:Lattice{N},TDF<:DielectricFunction,TPF<:PumpFunction,TC<:Tuple} where N

        new{N,TBND,TLAT,TDF,TPF,TC}(type,name,boundary,lattice,dielectric,pump,x,indices,ε,F,χ,interior,bulk,surface,corner,nnm,nnp)
    end

    # function (dom::Domain{N})(dielectric::DielectricFunction) where N
    #     return new{N,typeof(bnd.boundary),typeof(bnd.lattice),typeof(bnd.dielectric),typeof(bnd.pump)}(dom.type,
    #     dom.name,dom.boundary,dom.lattice,dielectric,dom.pump,
    #     dom.x,dom.indices,dom.ε,dom.F,dom.interior,dom.bulk,
    #     dom.surface,dom.corner,dom.nnm,dom.nnp)
    # end
    # function (dom::Domain{N})(pump::PumpFunction) where N
    #     return newnew{N,typeof(bnd.boundary),typeof(bnd.lattice),typeof(bnd.dielectric),typeof(bnd.pump)}(dom.type,
    #     dom.name,dom.boundary,dom.lattice,dom.dielectric,pump,
    #     dom.x,dom.indices,dom.ε,dom.F,dom.interior,dom.bulk,
    #     dom.surface,dom.corner,dom.nnm,dom.nnp)
    # end
    # function (dom::Domain{N})(boundary::Boundary) where N
    #     return newnew{N,typeof(bnd.boundary),typeof(bnd.lattice),typeof(bnd.dielectric),typeof(bnd.pump)}(dom.type,
    #     dom.name,boundary,dom.lattice,dom.dielectric,dom.pump,
    #     dom.x,dom.indices,dom.ε,dom.F,dom.interior,dom.bulk,
    #     dom.surface,dom.corner,dom.nnm,dom.nnp)
    # end
end

function Base.conj(d::Domain)
    Domain(
        d.type,
        d.name,
        conj(d.boundary),
        d.lattice,
        conj(d.dielectric),
        d.pump,
        d.x,
        d.indices,
        conj(d.ε),
        d.F,
        conj(d.χ),
        d.interior,
        d.bulk,
        d.surface,
        d.corner,
        d.nnm,
        d.nnp
        )
end

# # for modifying non-geometric parameters of an already existing domain
# function (dom::Domain)(args...)
#     for a ∈ args
#         dom = dom(a)
#     end
#     return dom
# end

include("1D/Domains.jl")
# include("2D/Domains.jl")
# include("3D/Domains.jl")

function Base.show(io::IO,dom::Domain)
    println(io,"Domain ",dom.type,": ",dom.name)
    println(io,"Number of sites: ",length(dom.x))
    println(io,"Number of bulk points: ",sum(dom.bulk))
    println(io,"Number of surface points: ",sum(dom.surface))
    println(io,"Number of corner points: ",sum(dom.corner))
    println(io)
    println(io,dom.dielectric)
    println(io)
    println(io,dom.pump)
    println(io)
    println(io,dom.boundary)
    println(io)
    print(io,dom.lattice)
end

# @recipe function f(d::Domain)
#     aspect_ratio --> 1
#     legend --> false
#     @series begin
#         seriestype --> :scatter
#         markersize --> MARKERSIZE_SCALE/sqrt(length(d.x))
#         markerstrokealpha --> 0
#         shape --> MARKERSHAPE
#         (d.x[d.bulk],d.y[d.bulk])
#     end
#     @series begin
#         seriestype --> :scatter
#         markersize --> MARKERSIZE_SCALE/sqrt(length(d.x))
#         markerstrokealpha --> 0
#         shape --> MARKERSHAPE
#         (d.x[d.surface],d.y[d.surface])
#     end
#     @series begin
#         seriestype --> :scatter
#         markersize --> MARKERSIZE_SCALE/sqrt(length(d.x))
#         markerstrokealpha --> 0
#         shape --> MARKERSHAPE
#         (d.x[d.corner],d.y[d.corner])
#     end
#     @series d.boundary
# end
#

end # module


text:
import Distributions

function logsumexp(arr)
    min_arr = maximum(arr)
    min_arr + log(sum(exp.(arr .- min_arr)))
end

"""
A scheme for sequential Monte Carlo inference in a state space model, without rejuvenation.
Concerete subtypes should implement the following methods:
    init(scheme::StateSpaceSMCScheme{H})
    init_score(scheme::StateSpaceSMCScheme{H}, state::H)
    forward(scheme::StateSpaceSMCScheme{H}, prev_state::H, t::Integer)
    forward_score(scheme::StateSpaceSMCScheme{H}, prev_state::H, state::H, t::Integer)
    get_num_steps(scheme::StateSpaceSMCScheme{H})
    get_num_particles(scheme::StateSpaceSMCScheme{H})
    get_ess_threshold(scheme::StateSpaceSMCScheme{H})
"""
abstract type StateSpaceSMCScheme{H} end

function init end
function init_score end
function forward end
function forward_score end
function get_num_steps end
function get_num_particles end
function get_ess_threshold end

struct StateSpaceSMCResult{H}
    states::Matrix{H}
    parents::Matrix{Int}
    log_weights::Vector{Float64}
    log_ml_estimate::Float64
    num_resamples::Int
end

function effective_sample_size(log_weights::Vector{Float64})
    # assumes weights are normalized
    log_ess = -logsumexp(2. * log_weights)
    exp(log_ess)
end

"""
Sequential Monte Carlo for state space models without rejuvenation and with multinomial resampling.
See algoritihm 3.1.1 of:
Del Moral, Pierre, Arnaud Doucet, and Ajay Jasra. "Sequential monte carlo samplers."
Journal of the Royal Statistical Society: Series B (Statistical Methodology) 68.3 (2006): 411-436.
Also see algorithm 2.2.1 of:
Andrieu, Christophe, Arnaud Doucet, and Roman Holenstein. "Particle markov chain monte carlo methods."
Journal of the Royal Statistical Society: Series B (Statistical Methodology) 72.3 (2010): 269-342.
"""
function smc(scheme::StateSpaceSMCScheme{H}) where {H}
    N = get_num_particles(scheme)
    T = get_num_steps(scheme)
    ess_threshold = get_ess_threshold(scheme)
    states = Matrix{H}(undef, N, T)
    parents = Matrix{Int}(undef, N, T-1)
    log_unnormalized_weights = Vector{Float64}(undef, N)
    log_ml_estimate = 0.
    for i=1:N
        (states[i, 1], log_unnormalized_weights[i]) = init(scheme)
    end

    num_resamples = 0
    for t=2:T
        log_total_weight = logsumexp(log_unnormalized_weights)
        log_normalized_weights = log_unnormalized_weights .- log_total_weight
        if effective_sample_size(log_normalized_weights) < ess_threshold
            weights = exp.(log_normalized_weights)
            parents[:, t-1] = rand(Distributions.Categorical(weights / sum(weights)), N)
            log_ml_estimate += log_total_weight - log(N)
            log_unnormalized_weights = zeros(N)
            num_resamples += 1
        else
            parents[:, t-1] = 1:N
        end
        for i=1:N
            parent = parents[i, t-1]
            (states[i, t], log_incremental_weight) = forward(scheme, states[parent, t-1], t)
            log_unnormalized_weights[i] += log_incremental_weight
        end
    end
    log_total_weight = logsumexp(log_unnormalized_weights)
    log_normalized_weights = log_unnormalized_weights .- log_total_weight
    log_ml_estimate += log_total_weight - log(N)
    StateSpaceSMCResult(states, parents, log_normalized_weights, log_ml_estimate, num_resamples)
end

# TODO modify conditional SMC below

const ONE = 1
"""
Conditional sequential Monte Carlo update for state space models without rejuvenation and with multinomial resampling.
See sections 2.4.3 and 4.3 of:
Andrieu, Christophe, Arnaud Doucet, and Roman Holenstein. "Particle markov chain monte carlo methods."
Journal of the Royal Statistical Society: Series B (Statistical Methodology) 72.3 (2010): 269-342.
"""

function conditional_smc(scheme::StateSpaceSMCScheme{H}, distinguished_particle::Vector{H}) where {H}
    N = get_num_particles(scheme)
    T = get_num_steps(scheme)
    ess_threshold = get_ess_threshold(scheme)
    if length(distinguished_particle) != T
        error("Expected particle length $T, actual length was $(length(distinguished_particle))")
    end
    states = Matrix{H}(undef, N, T)
    parents = Matrix{Int}(undef, N, T-1)
    log_unnormalized_weights = Vector{Float64}(undef, N)
    log_ml_estimate = 0.

    # Due to symmetries, the ancestral indices of the distinguished
    # particle do not matter, so we set them all to 1.
    # Note that this may not suffice for other resampling schemes.
    states[ONE, 1] = distinguished_particle[1]
    log_unnormalized_weights[ONE] = init_score(scheme, states[ONE, 1])

    for i=2:N
        (states[i, 1], log_unnormalized_weights[i]) = init(scheme)
    end

    num_resamples = 0
    for t=2:T
        log_total_weight = logsumexp(log_unnormalized_weights)
        log_normalized_weights = log_unnormalized_weights .- log_total_weight
        if effective_sample_size(log_normalized_weights) < ess_threshold
            weights = exp.(log_normalized_weights)
            parents[:, t-1] = rand(Distributions.Categorical(weights / sum(weights)), N)
            log_ml_estimate += log_total_weight - log(N)
            log_unnormalized_weights = zeros(N)
            num_resamples += 1
        else
            parents[:, t-1] = 1:N
        end

        # handle distinguished particle
        parents[ONE, t-1] = ONE
        states[ONE, t] = distinguished_particle[t]
        log_unnormalized_weights[ONE] += forward_score(scheme, states[ONE, t-1], states[ONE, t], t)

        for i=2:N
            parent = parents[i, t-1]
            (states[i, t], log_incremental_weight) = forward(scheme, states[parent, t-1], t)
            log_unnormalized_weights[i] += log_incremental_weight
        end
    end
    log_total_weight = logsumexp(log_unnormalized_weights)
    log_normalized_weights = log_unnormalized_weights .- log_total_weight
    log_ml_estimate += log_total_weight - log(N)
    StateSpaceSMCResult(states, parents, log_normalized_weights, log_ml_estimate, num_resamples)
end

function get_particle(result::StateSpaceSMCResult{H}, final_index::Integer) where {H}
    (N, T) = size(result.states)
    particle = Vector{H}(undef, T)
    particle[T] = result.states[final_index, T]
    current_index = final_index
    for t=T-1:-1:1
        current_index = result.parents[current_index, t]
        particle[t] = result.states[current_index, t]
    end
    return particle
end