text:
module GML_Glauber_Dynamics

export learn, inverse_ising
export learn_glauber_dynamics
export learn_glauber_dynamics_sho, learn_glauber_dynamics_sho2

export GMLFormulation, RISE, logRISE, RPLE, RISEA, multiRISE
export GMLMethod, NLP

using JuMP
using Ipopt
using Printf

import LinearAlgebra
import LinearAlgebra: diag
import LinearAlgebra: diagm
import Statistics: mean

include("models.jl")

include("sampling.jl")

abstract type GMLFormulation end

mutable struct multiRISE <: GMLFormulation
    regularizer::Real
    symmetrization::Bool
    interaction_order::Integer
end
# default values
multiRISE() = multiRISE(0.4, true, 2)

mutable struct RISE <: GMLFormulation
    regularizer::Real
    symmetrization::Bool
end
# default values
RISE() = RISE(0.4, true)

mutable struct RISEA <: GMLFormulation
    regularizer::Real
    symmetrization::Bool
end
# default values
RISEA() = RISEA(0.4, true)

mutable struct logRISE <: GMLFormulation
    regularizer::Real
    symmetrization::Bool
end
# default values
logRISE() = logRISE(0.8, true)

mutable struct RPLE <: GMLFormulation
    regularizer::Real
    symmetrization::Bool
end
# default values
RPLE() = RPLE(0.2, true)

abstract type GMLMethod end

mutable struct NLP <: GMLMethod
    #solver::JuMP.OptimizerFactory
    # Hard coded the update from the package online
    solver::Any
end
# default values
NLP() = NLP(with_optimizer(Ipopt.Optimizer, print_level=0, tol=1e-12))

# Other Solvers
# Coordinate descent
mutable struct CD <: GMLMethod
    solver_tolerance::Real
end
# default value
CD() = CD(1e-8)

# default settings
learn(samples::Array{T,2}) where T <: Real = learn(samples, RISE(), NLP())
learn(samples::Array{T,2}, formulation::S) where {T <: Real, S <: GMLFormulation} = learn(samples, formulation, NLP())

#TODO add better support for Adjoints
learn(samples::LinearAlgebra.Adjoint, args...) = learn(copy(samples), args...)

function data_info(samples::Array{T,2}) where T <: Real
    (num_conf, num_row) = size(samples)
    num_spins = num_row - 1
    num_samples = sum(samples[1:num_conf,1])
    return num_conf, num_spins, num_samples
end

function data_info_glauber_dynamics(samples::Array{T,2}) where T <: Real
    (num_conf, num_row) = size(samples)
    num_spins = num_row - 2
    num_spins = Int(num_spins/2)
    num_samples = sum(samples[1:num_conf,1])
    return num_conf, num_spins, num_samples
end

function learn(samples::Array{T,2}, formulation::multiRISE, method::NLP) where T <: Real
    num_conf, num_spins, num_samples = data_info(samples)

    lambda = formulation.regularizer*sqrt(log((num_spins^2)/0.05)/num_samples)
    inter_order = formulation.interaction_order

    reconstruction = Dict{Tuple,Real}()

    for current_spin = 1:num_spins
        nodal_stat = Dict{Tuple,Array{Real,1}}()

        for p = 1:inter_order
                nodal_keys = Array{Tuple{},1}()
                neighbours = [i for i=1:num_spins if i!=current_spin]
                if p == 1
                    nodal_keys = [(current_spin,)]
                else
                    perm = permutations(neighbours, p - 1)
                    if length(perm) > 0
                        nodal_keys = [(current_spin, perm[i]...) for i=1:length(perm)]
                    end
                end

                for index = 1:length(nodal_keys)
                    nodal_stat[nodal_keys[index]] =  [ prod(samples[k, 1 + i] for i=nodal_keys[index]) for k=1:num_conf]
                end
        end

        model = Model(method.solver)

        @variable(model, x[keys(nodal_stat)])
        @variable(model, z[keys(nodal_stat)])

        @NLobjective(model, Min,
            sum((samples[k,1]/num_samples)*exp(-sum(x[inter]*stat[k] for (inter,stat) = nodal_stat)) for k=1:num_conf) +
            lambda*sum(z[inter] for inter = keys(nodal_stat) if length(inter)>1)
        )

        for inter in keys(nodal_stat)
            @constraint(model, z[inter] >=  x[inter]) #z_plus
            @constraint(model, z[inter] >= -x[inter]) #z_minus
        end

        JuMP.optimize!(model)
        @assert JuMP.termination_status(model) == JuMP.MOI.LOCALLY_SOLVED

        nodal_reconstruction = JuMP.value.(x)
        for inter = keys(nodal_stat)
            reconstruction[inter] = deepcopy(nodal_reconstruction[inter])
        end
    end

    if formulation.symmetrization
        reconstruction_list = Dict{Tuple,Vector{Real}}()
        for (k,v) in reconstruction
            key = tuple(sort([i for i in k])...)
            if !haskey(reconstruction_list, key)
                reconstruction_list[key] = Vector{Real}()
            end
            push!(reconstruction_list[key], v)
        end

        reconstruction = Dict{Tuple,Real}()
        for (k,v) in reconstruction_list
            reconstruction[k] = mean(v)
        end
    end

    return FactorGraph(inter_order, num_spins, :spin, reconstruction)
end

function learn(samples::Array{T,2}, formulation::RISE, method::NLP) where T <: Real
    num_conf, num_spins, num_samples = data_info(samples)

    lambda = formulation.regularizer*sqrt(log((num_spins^2)/0.05)/num_samples)

    reconstruction = Array{Float64}(undef, num_spins, num_spins)

    for current_spin = 1:num_spins
        nodal_stat  = [ samples[k, 1 + current_spin] * (i == current_spin ? 1 : samples[k, 1 + i]) for k=1:num_conf , i=1:num_spins]

        model = Model(method.solver)

        @variable(model, x[1:num_spins])
        @variable(model, z[1:num_spins])

        @NLobjective(model, Min,
            sum((samples[k,1]/num_samples)*exp(-sum(x[i]*nodal_stat[k,i] for i=1:num_spins)) for k=1:num_conf) +
            lambda*sum(z[j] for j=1:num_spins if current_spin!=j)
        )

        for j in 1:num_spins
            @constraint(model, z[j] >=  x[j]) #z_plus
            @constraint(model, z[j] >= -x[j]) #z_minus
        end

        JuMP.optimize!(model)
        @assert JuMP.termination_status(model) == JuMP.MOI.LOCALLY_SOLVED
        reconstruction[current_spin,1:num_spins] = deepcopy(JuMP.value.(x))
    end

    if formulation.symmetrization
        reconstruction = 0.5*(reconstruction + transpose(reconstruction))
    end

    return reconstruction
end

function risea_obj(var, stat, weight)
    (num_conf, num_spins) = size(stat)
    #chvar = cosh.(var)
    #shvar = sinh.(var)
    #return sum(weight[k]*prod(chvar[i] - shvar[i]*stat[k,i] for i=1:num_spins) for k=1:num_conf)
    return sum(weight[k]*exp(-sum(var[i]*stat[k,i] for i=1:num_spins)) for k=1:num_conf)
end

function grad_risea_obj(g, var, stat, weight)
    (num_conf, num_spins) = size(stat)
    #chvar = cosh.(var)
    #shvar = sinh.(var)
    #partial_obj = [- weight[k] * prod(chvar[i] - shvar[i]*stat[k,i] for i=1:num_spins) for k=1:num_conf]
    partial_obj = [- weight[k]*exp(-sum(var[i]*stat[k,i] for i=1:num_spins)) for k=1:num_conf]
    for i=1:num_spins
        g[i] = sum(stat[k,i]*partial_obj[k] for k=1:num_conf)
    end
end

function learn(samples::Array{T,2}, formulation::RISEA, method::NLP) where T <: Real
    num_conf, num_spins, num_samples = data_info(samples)

    lambda = formulation.regularizer*sqrt(log((num_spins^2)/0.05)/num_samples)

    reconstruction = Array{Float64}(undef, num_spins, num_spins)

    for current_spin = 1:num_spins
        nodal_stat  = [ samples[k, 1 + current_spin] * (i == current_spin ? 1 : samples[k, 1 + i]) for k=1:num_conf , i=1:num_spins]
        weight = samples[1:num_conf,1] / num_samples

        obj(x...) = risea_obj(x, nodal_stat, weight)
        function grad(g, x...)
            grad_risea_obj(g, x, nodal_stat, weight)
        end

        function l1norm(z...)
            lambda*sum(z[j] for j=1:num_spins if current_spin!=j)
        end

        model = Model(method.solver)

        JuMP.register(model, :obj, num_spins, obj, grad)
        JuMP.register(model, :l1norm, num_spins, l1norm, autodiff=true)

        @variable(model, x[1:num_spins])
        @variable(model, z[1:num_spins])

        JuMP.setNLobjective(model, :Min, Expr(:call, :+,
                                            Expr(:call, :obj, x...),
                                            Expr(:call, :l1norm, z...)
                                        )
                            )

        for j in 1:num_spins
            @constraint(model, z[j] >=  x[j]) #z_plus
            @constraint(model, z[j] >= -x[j]) #z_minus
        end

        JuMP.optimize!(model)
        @assert JuMP.termination_status(model) == JuMP.MOI.LOCALLY_SOLVED
        reconstruction[current_spin,1:num_spins] = deepcopy(JuMP.value.(x))
    end

    if formulation.symmetrization
        reconstruction = 0.5*(reconstruction + transpose(reconstruction))
    end

    return reconstruction
end

function learn(samples::Array{T,2}, formulation::logRISE, method::NLP) where T <: Real
    num_conf, num_spins, num_samples = data_info(samples)

    lambda = formulation.regularizer*sqrt(log((num_spins^2)/0.05)/num_samples)

    reconstruction = Array{Float64}(undef, num_spins, num_spins)

    for current_spin = 1:num_spins
        nodal_stat  = [ samples[k, 1 + current_spin] * (i == current_spin ? 1 : samples[k, 1 + i]) for k=1:num_conf , i=1:num_spins]

        model = Model(method.solver)

        @variable(model, x[1:num_spins])
        @variable(model, z[1:num_spins])

        @NLobjective(model, Min,
            log(sum((samples[k,1]/num_samples)*exp(-sum(x[i]*nodal_stat[k,i] for i=1:num_spins)) for k=1:num_conf)) +
            lambda*sum(z[j] for j=1:num_spins if current_spin!=j)
        )

        for j in 1:num_spins
            @constraint(model, z[j] >=  x[j]) #z_plus
            @constraint(model, z[j] >= -x[j]) #z_minus
        end

        JuMP.optimize!(model)
        @assert JuMP.termination_status(model) == JuMP.MOI.LOCALLY_SOLVED
        reconstruction[current_spin,1:num_spins] = deepcopy(JuMP.value.(x))
    end

    if formulation.symmetrization
        reconstruction = 0.5*(reconstruction + transpose(reconstruction))
    end

    return reconstruction
end

function learn(samples::Array{T,2}, formulation::RPLE, method::NLP) where T <: Real
    num_conf, num_spins, num_samples = data_info(samples)

    lambda = formulation.regularizer*sqrt(log((num_spins^2)/0.05)/num_samples)

    reconstruction = Array{Float64}(undef, num_spins, num_spins)

    for current_spin = 1:num_spins
        nodal_stat  = [ samples[k, 1 + current_spin] * (i == current_spin ? 1 : samples[k, 1 + i]) for k=1:num_conf , i=1:num_spins]

        model = Model(method.solver)

        @variable(model, x[1:num_spins])
        @variable(model, z[1:num_spins])

        @NLobjective(model, Min,
            sum((samples[k,1]/num_samples)*log(1 + exp(-2*sum(x[i]*nodal_stat[k,i] for i=1:num_spins))) for k=1:num_conf) +
            lambda*sum(z[j] for j=1:num_spins if current_spin!=j)
        )

        for j in 1:num_spins
            @constraint(model, z[j] >=  x[j]) #z_plus
            @constraint(model, z[j] >= -x[j]) #z_minus
        end

        JuMP.optimize!(model)
        @assert JuMP.termination_status(model) == JuMP.MOI.LOCALLY_SOLVED
        reconstruction[current_spin,1:num_spins] = deepcopy(JuMP.value.(x))
    end

    if formulation.symmetrization
        reconstruction = 0.5*(reconstruction + transpose(reconstruction))
    end

    return reconstruction
end

## Solvers for learning Glauber Dynamics

# default settings in general
learn_glauber_dynamics(samples::Array{T,2}) where T <: Real = learn_glauber_dynamics(samples, RISE(), NLP())
learn_glauber_dynamics(samples::Array{T,2}, formulation::S) where {T <: Real, S <: GMLFormulation} = learn_glauber_dynamics(samples, formulation, NLP())

# default setting for coordinate descent
learn_glauber_dynamics(samples::Array{T,2}, formulation::S, method::M) where {T <: Real, S <: GMLFormulation, M <: GMLMethod} = learn_glauber_dynamics(samples, RISE(), CD())

## Using the JuMP Solvers
function learn_glauber_dynamics(samples::Array{T,2}, formulation::RISE, method::NLP) where T <: Real
    @info("using JuMP for RISE to learn Glauber dynamics")

    # samples are in the form of [num_samples matching config, config = node selected for update at t, \sigma^{t}, \sigma^{t+1}]
    num_conf, num_spins, num_samples = data_info_glauber_dynamics(samples)

    lambda = formulation.regularizer*sqrt(log((num_spins^2)/0.05)/num_samples)

    reconstruction = Array{Float64}(undef, num_spins, num_spins)

    for current_spin = 1:num_spins
        samples_current_spin = samples[findall(isequal(current_spin),samples[:,2]),:]
        num_conf_current_spin = size(samples_current_spin,1)
        num_samples_current_spin = sum(samples_current_spin[:,1])

        #display(samples_current_spin)

        nodal_stat  = [ samples_current_spin[k, 2 + num_spins + current_spin] * (i == current_spin ? 1 : samples_current_spin[k, 2 + i]) for k=1:num_conf_current_spin , i=1:num_spins]

        m = Model(method.solver)

        @variable(m, x[1:num_spins])
        @variable(m, z[1:num_spins])

        @NLobjective(m, Min,
            sum((samples_current_spin[k,1]/num_samples_current_spin)*exp(-sum(x[i]*nodal_stat[k,i] for i=1:num_spins)) for k=1:num_conf_current_spin) +
            lambda*sum(z[j] for j=1:num_spins if current_spin!=j)
        )

        for j in 1:num_spins
            @constraint(m, z[j] >=  x[j]) #z_plus
            @constraint(m, z[j] >= -x[j]) #z_minus
        end

        JuMP.optimize!(m)
        @assert JuMP.termination_status(m) == JuMP.MOI.LOCALLY_SOLVED
        reconstruction[current_spin,1:num_spins] = deepcopy(JuMP.value.(x))
    end

    if formulation.symmetrization
        reconstruction = 0.5*(reconstruction + transpose(reconstruction))
    end

    return reconstruction
end

function learn_glauber_dynamics(samples::Array{T,2}, formulation::logRISE, method::NLP) where T <: Real
    @info("using JuMP for logRISE")

    # samples are in the form of [num_samples matching config, config = node selected for update at t, \sigma^{t}, \sigma^{t+1}]
    num_conf, num_spins, num_samples = data_info_glauber_dynamics(samples)

    lambda = formulation.regularizer*sqrt(log((num_spins^2)/0.05)/num_samples)

    reconstruction = Array{Float64}(undef, num_spins, num_spins)

    for current_spin = 1:num_spins
        samples_current_spin = samples[find(isequal(current_spin),samples[:,2]),:]
        num_conf_current_spin = size(samples_current_spin,1)
        num_samples_current_spin = sum(samples_current_spin[:,1])

        #display(samples_current_spin)

        nodal_stat  = [ samples_current_spin[k, 2 + num_spins + current_spin] * (i == current_spin ? 1 : samples_current_spin[k, 2 + i]) for k=1:num_conf_current_spin , i=1:num_spins]

        m = Model(solver = method.solver)

        @variable(m, x[1:num_spins])
        @variable(m, z[1:num_spins])

        @NLobjective(m, Min,
            log(sum((samples_current_spin[k,1]/num_samples_current_spin)*exp(-sum(x[i]*nodal_stat[k,i] for i=1:num_spins)) for k=1:num_conf_current_spin)) +
            lambda*sum(z[j] for j=1:num_spins if current_spin!=j)
        )

        for j in 1:num_spins
            @constraint(m, z[j] >=  x[j]) #z_plus
            @constraint(m, z[j] >= -x[j]) #z_minus
        end

        status = solve(m)
        @assert status == :Optimal
        reconstruction[current_spin,1:num_spins] = deepcopy(getvalue(x))
    end

    if formulation.symmetrization
        reconstruction = 0.5*(reconstruction + transpose(reconstruction))
    end

    return reconstruction
end

function learn_glauber_dynamics(samples::Array{T,2}, formulation::RPLE, method::NLP) where T <: Real
    @info("using JuMP for RPLE to learn Glauber dynamics")

    # samples are in the form of [num_samples matching config, config = node selected for update at t, \sigma^{t}, \sigma^{t+1}]
    num_conf, num_spins, num_samples = data_info_glauber_dynamics(samples)

    lambda = formulation.regularizer*sqrt(log((num_spins^2)/0.05)/num_samples)

    reconstruction = Array{Float64}(undef, num_spins, num_spins)

    for current_spin = 1:num_spins
        samples_current_spin = samples[findall(isequal(current_spin),samples[:,2]),:]
        num_conf_current_spin = size(samples_current_spin,1)
        num_samples_current_spin = sum(samples_current_spin[:,1])

        #display(samples_current_spin)

        nodal_stat  = [ samples_current_spin[k, 2 + num_spins + current_spin] * (i == current_spin ? 1 : samples_current_spin[k, 2 + i]) for k=1:num_conf_current_spin , i=1:num_spins]

        m = Model(method.solver)

        @variable(m, x[1:num_spins])
        @variable(m, z[1:num_spins])

        @NLobjective(m, Min,
            sum((samples_current_spin[k,1]/num_samples_current_spin)*log(1 + exp(-2*sum(x[i]*nodal_stat[k,i] for i=1:num_spins))) for k=1:num_conf_current_spin) +
            lambda*sum(z[j] for j=1:num_spins if current_spin!=j)
        )

        for j in 1:num_spins
            @constraint(m, z[j] >=  x[j]) #z_plus
            @constraint(m, z[j] >= -x[j]) #z_minus
        end

        JuMP.optimize!(m)
        @assert JuMP.termination_status(m) == JuMP.MOI.LOCALLY_SOLVED
        reconstruction[current_spin,1:num_spins] = deepcopy(JuMP.value.(x))
    end

    if formulation.symmetrization
        reconstruction = 0.5*(reconstruction + transpose(reconstruction))
    end

    return reconstruction
end

function learn_glauber_dynamics_regularization(samples::Array{T,2}, formulation::RPLE, method::NLP) where T <: Real
    @info("using JuMP for RPLE with right regularization to learn Glauber dynamics")

    # samples are in the form of [num_samples matching config, config = node selected for update at t, \sigma^{t}, \sigma^{t+1}]
    num_conf, num_spins, num_samples = data_info_glauber_dynamics(samples)

    reconstruction = Array{Float64}(undef, num_spins, num_spins)

    for current_spin = 1:num_spins
        samples_current_spin = samples[findall(isequal(current_spin),samples[:,2]),:]
        num_conf_current_spin = size(samples_current_spin,1)
        num_samples_current_spin = sum(samples_current_spin[:,1])

        lambda = formulation.regularizer*sqrt(log((num_spins^2)/0.05)/num_samples_current_spin)

        nodal_stat  = [ samples_current_spin[k, 2 + num_spins + current_spin] * (i == current_spin ? 1 : samples_current_spin[k, 2 + i]) for k=1:num_conf_current_spin , i=1:num_spins]

        m = Model(method.solver)

        @variable(m, x[1:num_spins])
        @variable(m, z[1:num_spins])

        @NLobjective(m, Min,
            sum((samples_current_spin[k,1]/num_samples_current_spin)*log(1 + exp(-2*sum(x[i]*nodal_stat[k,i] for i=1:num_spins))) for k=1:num_conf_current_spin) +
            lambda*sum(z[j] for j=1:num_spins if current_spin!=j)
        )

        for j in 1:num_spins
            @constraint(m, z[j] >=  x[j]) #z_plus
            @constraint(m, z[j] >= -x[j]) #z_minus
        end

        JuMP.optimize!(m)
        @assert JuMP.termination_status(m) == JuMP.MOI.LOCALLY_SOLVED
        reconstruction[current_spin,1:num_spins] = deepcopy(JuMP.value.(x))
    end

    if formulation.symmetrization
        reconstruction = 0.5*(reconstruction + transpose(reconstruction))
    end

    return reconstruction
end

function learn_glauber_dynamics_regularization(samples::Array{T,2}, formulation::RISE, method::NLP) where T <: Real
    @info("using JuMP for RISE with right regularization to learn Glauber dynamics")
    @printf("c=%f\n", formulation.regularizer)

    # samples are in the form of [num_samples matching config, config = node selected for update at t, \sigma^{t}, \sigma^{t+1}]
    num_conf, num_spins, num_samples = data_info_glauber_dynamics(samples)

    reconstruction = Array{Float64}(undef, num_spins, num_spins)

    for current_spin = 1:num_spins
        samples_current_spin = samples[findall(isequal(current_spin),samples[:,2]),:]
        num_conf_current_spin = size(samples_current_spin,1)
        num_samples_current_spin = sum(samples_current_spin[:,1])

        lambda = formulation.regularizer*sqrt(log((num_spins^2)/0.05)/num_samples_current_spin)

        nodal_stat  = [ samples_current_spin[k, 2 + num_spins + current_spin] * (i == current_spin ? 1 : samples_current_spin[k, 2 + i]) for k=1:num_conf_current_spin , i=1:num_spins]

        m = Model(method.solver)

        @variable(m, x[1:num_spins])
        @variable(m, z[1:num_spins])

        @NLobjective(m, Min,
            sum((samples_current_spin[k,1]/num_samples_current_spin)*exp(-sum(x[i]*nodal_stat[k,i] for i=1:num_spins)) for k=1:num_conf_current_spin) +
            lambda*sum(z[j] for j=1:num_spins if current_spin!=j)
        )

        for j in 1:num_spins
            @constraint(m, z[j] >=  x[j]) #z_plus
            @constraint(m, z[j] >= -x[j]) #z_minus
        end

        JuMP.optimize!(m)
        @assert JuMP.termination_status(m) == JuMP.MOI.LOCALLY_SOLVED
        reconstruction[current_spin,1:num_spins] = deepcopy(JuMP.value.(x))
    end

    if formulation.symmetrization
        reconstruction = 0.5*(reconstruction + transpose(reconstruction))
    end

    return reconstruction
end

## Coordinate Descent Solvers for learning Glauber Dynamics (using analytical solution)
function learn_glauber_dynamics(samples::Array{T,2}, formulation::RISE, method::CD) where T <: Real
    @info("using coordinate descent solver for RISE to learn Glauber dyamics")

    # samples are in the form of [num_samples matching config, config = node selected for update at t, \sigma^{t}, \sigma^{t+1}]
    num_conf, num_spins, num_samples = data_info_glauber_dynamics(samples)

    lambda = formulation.regularizer*sqrt(log((num_spins^2)/0.05)/num_samples)

    reconstruction = Array{Float64}(num_spins, num_spins)

    # ISO Objective function
    RISEobjective(x) = exp(-x)

    for current_spin = 1:num_spins
        samples_current_spin = samples[find(isequal(current_spin),samples[:,2]),:]
        num_conf_current_spin = size(samples_current_spin,1)
        num_samples_current_spin = sum(samples_current_spin[:,1])

        conf_weight     = [samples_current_spin[k,1]/num_samples_current_spin for k=1:num_conf_current_spin]
        nodal_stat  = [ samples_current_spin[k, 2 + num_spins + current_spin] * (i == current_spin ? 1 : samples_current_spin[k, 2 + i]) for k=1:num_conf_current_spin , i=1:num_spins]

        # Initialize
        x = zeros(num_spins)    # Solution vector
        z = zeros(num_conf)     # Current sum for each sample in histogram

        func_new = sum(conf_weight[k]*(RISEobjective(z[k]))  for k=1:num_conf_current_spin) + lambda*sum(x[j] for j=1:num_spins if current_spin!=j)
        func_old = copy(func_new)

        # Counter for number of iterations
        n_iter = 0

        # Coordinate Descent Method
        diff_func_val = 1
        while diff_func_val > method.solver_tolerance
            # Choosing a cyclic variant
            cycle_updates = randperm(num_spins)

            # Update function value on should
            func_old = copy(func_new)

            for ind_spin = 1:num_spins
                x_old = copy(x)
                func_old = copy(func_new)
                coord = cycle_updates[ind_spin]

                a = sum(conf_weight[k]*(RISEobjective(z[k])/RISEobjective(x[coord]*nodal_stat[k,coord]))  for k=1:num_conf_current_spin)
                b = sum(conf_weight[k]*(RISEobjective(z[k])/RISEobjective(x[coord]*nodal_stat[k,coord]))*nodal_stat[k,coord]  for k=1:num_conf_current_spin)

                epsilon = b/a
                mu = (coord==current_spin ? 0.0 : lambda/a)

                x_update = RISE_soft_thresolding(epsilon, mu)
                z = [ z[k] + (x_update - x[coord])*nodal_stat[k,coord] for k=1:num_conf_current_spin ]
                x[coord] = x_update
            end
            # Update counter and function value
            n_iter = n_iter + 1
            func_new = sum(conf_weight[k]*(RISEobjective(z[k]))  for k=1:num_conf_current_spin) + lambda*sum(x[j] for j=1:num_spins if current_spin!=j)

            diff_func_val = abs(func_new - func_old)
            #diff_func_val = norm(grad_k,1)
        end
        # Asserting solution is optimal?

        reconstruction[current_spin,1:num_spins] = deepcopy(x)
    end

    if formulation.symmetrization
        reconstruction = 0.5*(reconstruction + transpose(reconstruction))
    end

    return reconstruction
end

## Solvers for learning Glauber Dynamics from Spin-History-Only samples
# default settings
learn_glauber_dynamics_sho(samples::Dict{Int,Array{T,2}}, num_samples::Int) where T <: Real = learn_glauber_dynamics_sho(samples, num_samples, RISE(), NLP())
learn_glauber_dynamics_sho(samples::Dict{Int,Array{T,2}}, num_samples::Int, formulation::S) where {T <: Real, S <: GMLFormulation} = learn_glauber_dynamics_sho(samples, num_samples, formulation, NLP())

function learn_glauber_dynamics_sho(samples::Dict{Int,Array{T,2}}, num_samples::Int, formulation::RISE, method::NLP) where T <: Real
    @info("using JuMP for RISE to learn Glauber dynamics from SHO samples")

    # Dictonary over samples
    # samples[1] are in the form of [num_samples matching config, config = node selected for update at t, \sigma^{t}, \sigma^{t+1}]
    num_spins = length(keys(samples))

    lambda = formulation.regularizer*sqrt(log((num_spins^2)/0.05)/num_samples)

    reconstruction = Array{Float64}(undef, num_spins, num_spins)

    for current_spin = 1:num_spins
        samples_current_spin = samples[current_spin]
        samples_current_spin_flip = samples_current_spin[findall(isequal(1),samples_current_spin[:,2]),:]

        num_conf_current_spin_flip = size(samples_current_spin_flip,1)
        num_samples_current_spin_flip = sum(samples_current_spin_flip[:,1])

        #display(samples_current_spin)

        nodal_stat  = [ samples_current_spin_flip[k, 2 + num_spins + current_spin] * (i == current_spin ? 1 : samples_current_spin_flip[k, 2 + i]) for k=1:num_conf_current_spin_flip , i=1:num_spins]

        m = Model(method.solver)

        @variable(m, x[1:num_spins])
        @variable(m, z[1:num_spins])

        @NLobjective(m, Min,
            sum((samples_current_spin_flip[k,1]/num_samples_current_spin_flip)*exp(-sum(x[i]*nodal_stat[k,i] for i=1:num_spins)) for k=1:num_conf_current_spin_flip) +
            lambda*sum(z[j] for j=1:num_spins if current_spin!=j)
        )

        for j in 1:num_spins
            @constraint(m, z[j] >=  x[j]) #z_plus
            @constraint(m, z[j] >= -x[j]) #z_minus
        end

        JuMP.optimize!(m)
        @assert JuMP.termination_status(m) == JuMP.MOI.LOCALLY_SOLVED
        reconstruction[current_spin,1:num_spins] = deepcopy(JuMP.value.(x))
    end

    if formulation.symmetrization
        reconstruction = 0.5*(reconstruction + transpose(reconstruction))
    end

    return reconstruction
end

function learn_glauber_dynamics_sho(samples::Dict{Int,Array{T,2}}, num_samples::Int, formulation::RPLE, method::NLP) where T <: Real
    @info("using JuMP for RPLE to learn Glauber dynamics from SHO samples")

    # Dictonary over samples
    # samples[1] are in the form of [num_samples matching config, config = node selected for update at t, \sigma^{t}, \sigma^{t+1}]
    num_spins = length(keys(samples))

    lambda = formulation.regularizer*sqrt(log((num_spins^2)/0.05)/num_samples)

    reconstruction = Array{Float64}(undef, num_spins, num_spins)

    for current_spin = 1:num_spins
        # Using i in variable names to denote current spin
        samples_i = samples[current_spin]
        samples_i_flip = samples_i[findall(isequal(1),samples_i[:,2]),:]
        samples_i_no_flip = samples_i[findall(isequal(0),samples_i[:,2]),:]

        num_conf_i_flip = size(samples_i_flip,1)
        num_samples_i_flip = sum(samples_i_flip[:,1])

        num_conf_i_no_flip = size(samples_i_no_flip,1)
        num_samples_i_no_flip = sum(samples_i_no_flip[:,1])

        nodal_stat_flip  = [ samples_i_flip[k, 2 + num_spins + current_spin] * (i == current_spin ? 1 : samples_i_flip[k, 2 + i]) for k=1:num_conf_i_flip , i=1:num_spins]
        nodal_stat_no_flip  = [ samples_i_no_flip[k, 2 + num_spins + current_spin] * (i == current_spin ? 1 : samples_i_no_flip[k, 2 + i]) for k=1:num_conf_i_no_flip , i=1:num_spins]

        # Probability that i was chosen for updating
        gamma_i = 1.0/num_spins

        m = Model(method.solver)

        @variable(m, x[1:num_spins])
        @variable(m, z[1:num_spins])

        @NLobjective(m, Min,
            sum((samples_i_flip[k,1]/num_samples_i_flip)*log(1 + exp(-2*sum(x[i]*nodal_stat_flip[k,i] for i=1:num_spins))) for k=1:num_conf_i_flip) +
            lambda*sum(z[j] for j=1:num_spins if current_spin!=j)
        )

        for j in 1:num_spins
            @constraint(m, z[j] >=  x[j]) #z_plus
            @constraint(m, z[j] >= -x[j]) #z_minus
        end

        JuMP.optimize!(m)
        @assert JuMP.termination_status(m) == JuMP.MOI.LOCALLY_SOLVED
        reconstruction[current_spin,1:num_spins] = deepcopy(JuMP.value.(x))
    end

    if formulation.symmetrization
        reconstruction = 0.5*(reconstruction + transpose(reconstruction))
    end

    return reconstruction
end

learn_glauber_dynamics_sho2(samples::Dict{Int,Array{T,2}}, num_samples::Int) where T <: Real = learn_glauber_dynamics_sho2(samples, num_samples, RISE(), NLP())
learn_glauber_dynamics_sho2(samples::Dict{Int,Array{T,2}}, num_samples::Int, formulation::S) where {T <: Real, S <: GMLFormulation} = learn_glauber_dynamics_sho2(samples, num_samples, formulation, NLP())

function learn_glauber_dynamics_sho2(samples::Dict{Int,Array{T,2}}, num_samples::Int, formulation::RISE, method::NLP) where T <: Real
    @info("using JuMP for RISE to learn Glauber dynamics from SHO samples v2")

    # Dictonary over samples
    # samples[1] are in the form of [num_samples matching config, config = node selected for update at t, \sigma^{t}, \sigma^{t+1}]
    num_spins = length(keys(samples))

    lambda = formulation.regularizer*sqrt(log((num_spins^2)/0.05)/num_samples)

    reconstruction = Array{Float64}(undef, num_spins, num_spins)

    for current_spin = 1:num_spins
        # Using i in variable names to denote current spin
        samples_i = samples[current_spin]
        samples_i_flip = samples_i[findall(isequal(1),samples_i[:,2]),:]
        samples_i_no_flip = samples_i[findall(isequal(0),samples_i[:,2]),:]

        num_conf_i_flip = size(samples_i_flip,1)
        num_samples_i_flip = sum(samples_i_flip[:,1])

        num_conf_i_no_flip = size(samples_i_no_flip,1)
        num_samples_i_no_flip = sum(samples_i_no_flip[:,1])

        nodal_stat_flip  = [ samples_i_flip[k, 2 + num_spins + current_spin] * (i == current_spin ? 1 : samples_i_flip[k, 2 + i]) for k=1:num_conf_i_flip , i=1:num_spins]
        nodal_stat_no_flip  = [ samples_i_no_flip[k, 2 + num_spins + current_spin] * (i == current_spin ? 1 : samples_i_no_flip[k, 2 + i]) for k=1:num_conf_i_no_flip , i=1:num_spins]

        # Probability that i was chosen for updating
        gamma_i = 1.0/num_spins

        m = Model(method.solver)

        @variable(m, x[1:num_spins])
        @variable(m, z[1:num_spins])

        @NLobjective(m, Min,
            sum((samples_i_flip[k,1]/num_samples_i_flip)*exp(-sum(x[i]*nodal_stat_flip[k,i] for i=1:num_spins)) for k=1:num_conf_i_flip) +
            sum((samples_i_no_flip[k,1]/num_samples_i_no_flip)*(gamma_i*exp(-sum(x[i]*nodal_stat_no_flip[k,i] for i=1:num_spins)) + (1.0-gamma_i)) for k=1:num_conf_i_no_flip) +
            lambda*sum(z[j] for j=1:num_spins if current_spin!=j)
        )

        for j in 1:num_spins
            @constraint(m, z[j] >=  x[j]) #z_plus
            @constraint(m, z[j] >= -x[j]) #z_minus
        end

        JuMP.optimize!(m)
        @assert JuMP.termination_status(m) == JuMP.MOI.LOCALLY_SOLVED
        reconstruction[current_spin,1:num_spins] = deepcopy(JuMP.value.(x))
    end

    if formulation.symmetrization
        reconstruction = 0.5*(reconstruction + transpose(reconstruction))
    end

    return reconstruction
end

function learn_glauber_dynamics_sho2(samples::Dict{Int,Array{T,2}}, num_samples::Int, formulation::RPLE, method::NLP) where T <: Real
    @info("using JuMP for RPLE to learn Glauber dynamics from SHO samples v2")

    # Dictonary over samples
    # samples[1] are in the form of [num_samples matching config, config = node selected for update at t, \sigma^{t}, \sigma^{t+1}]
    num_spins = length(keys(samples))

    lambda = formulation.regularizer*sqrt(log((num_spins^2)/0.05)/num_samples)

    reconstruction = Array{Float64}(undef, num_spins, num_spins)

    for current_spin = 1:num_spins
        # Using i in variable names to denote current spin
        samples_i = samples[current_spin]
        #samples_i_flip = samples_i[find(isequal(1),samples_i[:,2]),:]
        samples_i_no_flip = samples_i[findall(isequal(0),samples_i[:,2]),:]

        #num_conf_i_flip = size(samples_i_flip,1)
        #num_samples_i_flip = sum(samples_i_flip[:,1])

        num_conf_i_no_flip = size(samples_i_no_flip,1)
        num_samples_i_no_flip = sum(samples_i_no_flip[:,1])

        num_conf_i = size(samples_i,1)
        num_samples_i = sum(samples_i[:,1])

        nodal_stat  = [ samples_i[k, 2 + num_spins + current_spin] * (i == current_spin ? 1 : samples_i[k, 2 + i]) for k=1:num_conf_i, i=1:num_spins]
        #nodal_stat_flip  = [ samples_i_flip[k, 2 + num_spins + current_spin] * (i == current_spin ? 1 : samples_i_flip[k, 2 + i]) for k=1:num_conf_i_flip , i=1:num_spins]
        nodal_stat_no_flip  = [ samples_i_no_flip[k, 2 + num_spins + current_spin] * (i == current_spin ? 1 : samples_i_no_flip[k, 2 + i]) for k=1:num_conf_i_no_flip , i=1:num_spins]

        # Probability that i was chosen for updating
        gamma_i = 1.0/num_spins

        m = Model(method.solver)

        @variable(m, x[1:num_spins])
        @variable(m, z[1:num_spins])

        @NLobjective(m, Min,
            sum((samples_i[k,1]/num_samples_i)*log(1 + exp(-2*sum(x[i]*nodal_stat[k,i] for i=1:num_spins))) for k=1:num_conf_i) -
            #sum((samples_i_no_flip[k,1]/num_samples_i_no_flip)*log(1 + (1.0-gamma_i)*exp(-2*sum(x[i]*nodal_stat_no_flip[k,i] for i=1:num_spins))) for k=1:num_conf_i_no_flip) +
            sum((samples_i_no_flip[k,1]/num_samples_i_no_flip)*(1.0-gamma_i)*log(1 + exp(-2*sum(x[i]*nodal_stat_no_flip[k,i] for i=1:num_spins))) for k=1:num_conf_i_no_flip) +
            lambda*sum(z[j] for j=1:num_spins if current_spin!=j)
        )

        for j in 1:num_spins
            @constraint(m, z[j] >=  x[j]) #z_plus
            @constraint(m, z[j] >= -x[j]) #z_minus
        end

        JuMP.optimize!(m)
        @assert JuMP.termination_status(m) == JuMP.MOI.LOCALLY_SOLVED
        reconstruction[current_spin,1:num_spins] = deepcopy(JuMP.value.(x))
    end

    if formulation.symmetrization
        reconstruction = 0.5*(reconstruction + transpose(reconstruction))
    end

    return reconstruction
end

end


text:
using Logging

# try
#   using PyPlot
# catch e
#   println("caught: ", e)
# end

function checkDerivative(f,df,x;doPlot::Bool=false)
    N = 17
    h = exp10.(2:-1:-N)
    v = randn(size(x)).-0.5
    fx = f(x)
    errquad = zeros(length(h))
    errlin  = zeros(length(h))
    for i=1:length(h)
        errlin[i]  = norm(fx             - f(x+h[i]*v))
        errquad[i] = norm(fx + h[i]*df*v - f(x+h[i]*v))
        s = @sprintf(" h: %5e  ||f(x+h*v)||: %5e   elin: %5e   equad: %5e \n",h[i], norm(f(x+h[i]*v)), errlin[i],errquad[i])
        Logging.info(s)
    end
    #
    # if(doPlot)
    #   PyPlot.rc("legend",fontsize=10)
    #   PyPlot.rc("xtick", labelsize=10, color="black", direction="in")
    #   PyPlot.rc("ytick", labelsize=10, color="black", direction="in")
    #   fig = figure(figsize=(6,5),facecolor="white")
    #   plot(h,errlin,label=L"$||f(x) - f(x+hv)||$")
    #   plot(h,errquad,".-",label=L"$||f(x) + h \nabla f v - f(x+hv)||$")
    #   yscale("log")
    #   xscale("log")
    #   xlabel(L"$h$")
    #   ylabel(L"$error$")
    #   legend(loc="upper left")
    # end
    return errlin, errquad
end

function checkErrorDecay(errquad::Vector)
  errorReduction = log10.(errquad[1:end-1]./errquad[2:end])
  errorReduction[isinf.(errorReduction)]=1
  return sum(errorReduction[errorReduction.>1.7])>6
end

ext: jl


text:
module GFlops
import Statistics
using BenchmarkTools:   @benchmark
using InteractiveUtils: peakflops
using Printf:           @printf

export @gflops, @count_ops

include("overdub.jl")
include("count_ops.jl")

end # module


text:
abstract type AbstractTransmissionNetwork end

const AbstractTN = AbstractTransmissionNetwork

function Base.sum(x::TN) where {TN <: AbstractTN}
  return sum(x.external) + sum(x.internal)
end

function Base.show(io::IO, x::TN) where {TN <: AbstractTN}
  return print(io, "$TN Σexternal = $(sum(x.external)) Σinternal = $(sum(x.internal))" )
end

function Base.copy(x::TN) where {TN <: AbstractTN}
  return TN(copy(x.external),
            copy(x.internal))
end

function individuals(x::TN) where {TN <: AbstractTN}
  return length(x.external)
end

function Base.sum(x::TN, i::Integer) where {TN <: AbstractTN}
  if i < 1 || i > individuals(x)
    @error "Invalid individual identifier" i
  end
  return x.external[i] + sum(x.internal[:, i])
end


text:
## broker/exchange communication

## IMPORTANT:
## careful error handling needed here not to crash live trading session!
## try-catch nearly everything!

### TODO how to plug-in function sets for different brokers/exchanges
### for now: need to overwrite ALL these functions
### to start a trading instance

"Backtesting helper function: check if limit-price is reached"
function plimitcheck(orde::Order, pnow::Float64)
  if (orde.side == :buy) && (pnow <= orde.price)
    return true
  end
  if (orde.side == :sell) && (pnow >= orde.price)
    return true
  end
  return false
end

"""
Get order status by order ID string.
Returns `Symbol` in line with `Order`-type options for status-slot.
"""
:query_orderstatus

"Order status: backtesting version based on current price `pnow`"
function query_orderstatus(orde::Order, pnow::Float64)
  if orde.ordertype == :market
    return :complete
  end
  if plimitcheck(orde, pnow)
    # limit-price reached
    return :complete
  end
  return :pending
end

"Order status: live version"
function query_orderstatus(ordid::String)
  ### TODO
  return resp_orderstatus
end

"Submit new order. Returns order ID string or `FAIL`-string"
:submit_ordernew

"New order submission: backtesting version."
function submit_ordernew(orde::Order, backtest::Bool)
  if !backtest
    return "FAIL"
  end
  #println("Submitting order $orde")
  return orde.id
end

"New order submission: live version"
function submit_ordernew(orde::Order, position_actual::Int64)
  ### TODO
  ### add short-sell logic here by checking position_actual
  return resp_orderid # if error return "FAIL"
end

"Cancel order request. Returns `Bool` request result."
:submit_ordercancel

"Cancel pending order backtest version"
function submit_ordercancel(orde::Order)
  if orde.status != :pending
    warn("Attempting to cancel a non-pending order")
    return false
  end
  #println("Cancelling pending order $orde")
  return true
end

"Cancel order live version: provide order ID string `ordid`"
function submit_ordercancel(ordid::String)
  ### TODO
  return resp_status
end


text:
const DEFAULT_SEED = 0

"""
    AST.ASTParams

Adaptive Stress Testing specific simulation parameters.
"""
@with_kw mutable struct ASTParams
    max_steps::Int64 = 0 # Maximum simulation time step (for runaways in simulation)
    seed::UInt32 = DEFAULT_SEED # Initial seed value
    reset_seed::Union{Nothing, Int64} = nothing # Reset to this seed value on initialize()
    top_k::Int64 = 0 # Number of top performing paths to save (defaults to 0, i.e. do not record)
    debug::Bool = false # Flag to indicate debugging mode (i.e. metrics collection, etc)
    episodic_rewards::Bool = false # decision making process with epsidic rewards
    give_intermediate_reward::Bool = false # give log-probability as reward during intermediate gen calls (used only if `episodic_rewards`)
    reward_bonus::Float64 = episodic_rewards ? 100 : 0 # reward received when event is found, multiplicative when using `episodic_rewards`
    use_potential_based_shaping::Bool = true # apply potential-based reward shaping to speed up learning
    pass_seed_action::Bool = false # pass the selected RNG seed to the GrayBox.transition! and BlackBox.evaluate! functions
    collect_data::Bool = false # flag to indicate supervised dataset collection of (𝐱=disturbances, y=isevent)
    discount::Float64 = 1.0 # discount factor (generally 1.0 to not discount later samples in the trajectory)
end
ASTParams(max_steps::Int64, seed::Int64) = ASTParams(max_steps=max_steps, seed=seed)
ASTParams(max_steps::Int64, seed::Int64, top_k::Int64, debug::Bool) = ASTParams(max_steps=max_steps, seed=seed, top_k=top_k, debug=debug)

"""
Abstract type for the AST action variants.
"""
abstract type ASTAction end

"""
    AST.ASTSeedAction

Random seed AST action.
"""
@with_kw mutable struct ASTSeedAction <: ASTAction
    seed::UInt32 = DEFAULT_SEED
end

import Base.string
string(a::ASTSeedAction) = "0x" * string(a.seed, base=16)

"""
    AST.ASTSampleAction

Random environment sample as the AST action.
"""
@with_kw mutable struct ASTSampleAction <: ASTAction
    sample::Union{GrayBox.EnvironmentSample, Nothing} = nothing
end

function string(a::ASTSampleAction)
    strings = ["$k => $(a.sample[k].value) ($(a.sample[k].logprob))" for k in keys(a.sample)]
    return join(strings, ',')
end

"""
    AST.ASTState

State of the AST MDP.
"""
@with_kw mutable struct ASTState
    t_index::Int64 = 0 # Confidence check that time corresponds
    parent::Union{Nothing,ASTState} = nothing # Parent state, `nothing` if root
    action::Union{Nothing,ASTAction} = nothing # Action taken from parent
    state::GrayBox.State = nothing # State of the gray-box simulation (not required/available for all problems)
    hash::UInt64 = hash(t_index, parent, action) # Hash simulation state to match with ASTState
    q_value::Float64 = 0.0 # Saved Q-value
    terminal::Bool = false # Indication of termination state
end

"""
    AST.ASTMetrics

Debugging metrics.
"""
@with_kw mutable struct ASTMetrics
    miss_distance = Real[]
    rate = Real[]
    logprob = Real[]
    prob = Real[]
    reward = Real[]
    intermediate_reward = Real[] # for computing returns at episode termination
    returns = Vector{Real}[]
    event = Bool[]
    terminal = Bool[]
end

"""
    AST.ASTMDP

Adaptive Stress Testing MDP problem formulation object.
"""
@with_kw mutable struct ASTMDP{Action<:ASTAction} <: MDP{ASTState, Action}
    params::ASTParams = ASTParams() # AST simulation parameters
    sim::GrayBox.Simulation # Black-box simulation struct
    sim_hash::UInt64 = hash(0) # Hash to keep simulations in sync
    t_index::Int64 = 1 # Simulation time
    current_seed::Union{Nothing, UInt32} = nothing # Rolling seed (`nothing` indicates initial seed comes from params.seed)

    top_paths::PriorityQueue{Any, Float64} = PriorityQueue{Any, Float64}(Base.Order.Forward) # Collection of best paths in the tree
    metrics::ASTMetrics = ASTMetrics() # Debugging metrics

    dataset::Vector = [] # (𝐱=disturbances, y=isevent) supervised dataset

    rate::Real = NaN # Stored current rate value

    predict::Union{Function, Nothing} = nothing # failure prediction function (negative = non-failures, positive = failures)
end

ASTMDP{Action}(sim::GrayBox.Simulation) where {Action<:ASTAction} = ASTMDP{Action}(sim=sim)
ASTMDP{Action}(params::ASTParams, sim::GrayBox.Simulation) where {Action<:ASTAction} = ASTMDP{Action}(params=params, sim=sim)

# Default to ASTSeedAction
ASTMDP(sim::GrayBox.Simulation) = ASTMDP{ASTSeedAction}(sim=sim)
ASTMDP(params::ASTParams, sim::GrayBox.Simulation) = ASTMDP{ASTSeedAction}(params=params, sim=sim)


text:
using PartialRejectionSampling
const PRS = PartialRejectionSampling

c, r = rand(2), 10

using Random
rng = Random.MersenneTwister(123)
@time pts = rand(rng, PRS.BallWindow(c, r), 10)
@time pts = rand(PRS.BallWindow(c, r), 10)

using Plots

Plots.scatter(pts[1, :], pts[2, :], aspect_ratio=:equal)


text:
module DXFutils

using DataStructures
using Printf
import Base.parse

# This package:
include("all_subtypes.jl")
include("dxf_groups.jl")
include("reader.jl")
include("parser.jl")

include("explore.jl")

end


text:
function _epidemic_curve_distribution(
  x::Vector{Events{S}},
  state::DiseaseState,
  times;
  credibleinterval::Float64=0.95) where {
    S <: DiseaseStateSequence}
  !(0.0 <= credibleinterval <= 1.0) && error("`credibleinterval` must be between 0 and 1")
  ci = (1-credibleinterval)/2

  lwr  = Array{Tuple{Float64,Float64}, 1}(undef, length(times))
  upr  = Array{Tuple{Float64,Float64}, 1}(undef, length(times))
  for t in eachindex(times)
    y = Array{Int64, 1}(undef, length(x))
    @simd for i in eachindex(x)
      y[i] = _count_by_state(x[i], state, times[t])
    end
    q = quantile(y, [ci, 1-ci])
    lwr[t] = (times[t], q[1])
    upr[t] = (times[t], q[2])
  end
  return [lwr; reverse(upr)]
end

@recipe function f(
  x::Vector{Events{S}},
  state::DiseaseState,
  times;
  credibleinterval=0.95,
  thin=1,
  burnin=0)  where {
    S <: DiseaseStateSequence}
  xguide --> "Time"
  yguide --> "N"
  label --> convert(Char, state)
  legend --> :topright
  fillcolor --> _state_color(state)
  linecolor --> _state_color(state)
  fillalpha --> 0.5
  linealpha --> 0.7
  seriestype := :shape
  _epidemic_curve_distribution(x[1+burnin:thin:end], state, times, credibleinterval=credibleinterval)
end

@recipe function f(
  x::Vector{Events{S}},
  times;
  credibleinterval=0.95,
  thin=1,
  burnin=0) where {
    S <: DiseaseStateSequence}
  for s in convert(Tuple, S)
    @series begin
      credibleinterval := credibleinterval
      burnin := burnin
      thin := thin
      x, s, times
    end
  end
end

@recipe function f(
  x::MCMC,
  state::DiseaseState,
  times;
  credibleinterval=0.95,
  thin=1,
  burnin=0,
  bychain=false)
  if !bychain
    y = x.markov_chains[1].events[1+burnin:thin:end]
    for i in 2:length(x.markov_chains)
      append!(y, x.markov_chains[i].events[1+burnin:thin:end])
    end
    @series begin
      xguide --> "Time"
      yguide --> "N"
      legend --> :topright
      label --> convert(Char, state)
      fillcolor --> _state_color(state)
      linecolor --> _state_color(state)
      fillalpha --> 0.5
      linealpha --> 0.7
      seriestype := :shape
      _epidemic_curve_distribution(y, state, times, credibleinterval=credibleinterval)
    end
  else
    label --> convert(Char, state)
    for i in eachindex(x.markov_chains)
      @series begin
        burnin := burnin
        thin := thin
        credibleinterval := credibleinterval
        label := :none
        x.markov_chains[i].events, state, times
      end
    end
  end
end

@recipe function f(
  x::MCMC{S, M},
  times;
  credibleinterval=0.95,
  thin=1,
  burnin=0,
  bychain=false) where {
    S <: DiseaseStateSequence,
    M <: ILM}
  for s in convert(Tuple, S)
    @series begin
      credibleinterval := credibleinterval
      burnin := burnin
      thin := thin
      bychain := bychain
      x, s, times
    end
  end
end


text:
# julia> include("/Users/ian/jd/10-19 Code/12 Julia/12.03 PE/pe3.jl")

function sieve_e(n)
    ints = 2:n; p = 2; i = 1
    while ints[end] != p
        ints = setdiff(ints,[j*p for j=2:Int(ceil(n/p))])
        p = ints[i + 1]; i += 1
    end
    return ints
end

function gpd(n)
    test = []
    p = sieve_e(Int(ceil(sqrt(n))))
    for i in p
        if n%i==0
            push!(test,i)
        end
    end
    return test[end]
end

println("ans = $(gpd(600851475143))")

max_stars_repo_path: learning exercises/pe/pe3.jl

max_stars_repo_name: ian-double-u/NNfromScratch

max_issues_repo_path: learning exercises/pe/pe3.jl

max_issues_repo_name: ian-double-u/NNfromScratch

max_forks_repo_path: learning exercises/pe/pe3.jl

max_forks_repo_name: ian-double-u/NNfromScratch


text:
using Images, ImageMagick, ImageFiltering
using Plots

file_name_in = joinpath("data", "turtle.jpg")
file_name_out = joinpath("data", "turtles.png")

img = load(file_name_in)

K1 = Float64.([1 1 1; 1 1 1; 1 1 1] ./ 9)
img1 = imfilter(img, K1)

K2 = Float64.([-1 -1 -1; -1 8 -1; -1 -1 -1])
img2 = imfilter(img, K2)

function fix_bounds!(img)
    for i in 1:size(img,1)
        for j in 1:size(img,2)
            a1 = min(max(img[i,j].r, 0), 1)
            a2 = min(max(img[i,j].g, 0), 1)
            a3 = min(max(img[i,j].b, 0), 1)
            img[i,j] = RGB(a1,a2,a3)
        end
    end
end

fix_bounds!(img1)
save("turtle1.jpg", img1)

fix_bounds!(img2)
save("turtle2.jpg", img2)

p = plot(
    axis = nothing,
    layout = @layout([a b c]),
    size=1.5.*(700,300)
)

plot!(p[1], img, ratio=1)
plot!(p[2], img1, ratio=1)
plot!(p[3], img2, ratio=1)

savefig(p, file_name_out)

max_stars_repo_name: rzehumat/JuliaCourse

max_issues_repo_name: rzehumat/JuliaCourse

max_forks_repo_name: rzehumat/JuliaCourse


text:
"""
    function fox()

# Example

```jldoctest
julia> cowsay("What does the fox say?", cow=Cowsay.fox)
 ________________________
< What does the fox say? >
 ------------------------
         \\     ,-.      .-,
          \\    |-.\\ __ /.-|
           \\   \\  `    `  /
                /_     _ \\
              <  _`q  p _  >
              <.._=/  \\=_. >
                 {`\\()/`}`\\
                 {      }  \\
                 |{    }    \\
                 \\ '--'   .- \\
                 |-      /    \\
                 | | | | |     ;
                 | | |.;.,..__ |
               .-"";`         `|
              /    |           /
              `-../____,..---'`
```
"""
function fox(;eyes="oo", tongue="  ", thoughts="\\")
the_cow =
"""
         $thoughts     ,-.      .-,
          $thoughts    |-.\\ __ /.-|
           $thoughts   \\  `    `  /
                /_     _ \\
              <  _`q  p _  >
              <.._=/  \\=_. >
                 {`\\()/`}`\\
                 {      }  \\
                 |{    }    \\
                 \\ '--'   .- \\
                 |-      /    \\
                 | | | | |     ;
                 | | |.;.,..__ |
               .-"";`         `|
              /    |           /
              `-../____,..---'`
"""
return the_cow
end


text:
abstract type AbstractDepolymerization end

@columns struct Depolymerization{X0,VE,AKIN,KEP} <: AbstractDepolymerization
    Χ_0::X0     | _         | "Number of C atoms per product molecule"
    V_E::VE     | 1/h       | "Specific hydrolysis rate"
    α_kin::AKIN | mol/mol   | "Enzyme binding sites per polymer molecule"
    K_EP::KEP   | mol/m^3   | "Depolymerization half-saturation constant"
end

for fn in fieldnames(Depolymerization)
    @eval $fn(p::Depolymerization) = p.$fn
end

depolymerization!(J_P::Vector{Float64}, p::Depolymerization, P::Vector{Float64}, E::Vector{Float64}) = begin
    n_substrates = size(P,1)
    n_consumers  = size(E,1)
    ECA  = ECA_kinetics!(zeros(n_substrates,n_consumers), P, E, K_EP(p), V_E(p), α_kin(p), "depoly")
    J_P = vcat(sum(ECA[1:n_substrates, :], dims=2)...)
    return J_P
end

max_stars_repo_name: giannamars/DEBMicroTrait

max_issues_repo_name: giannamars/DEBMicroTrait

max_forks_repo_name: giannamars/DEBMicroTrait


text:
function optimalAssignments(estimator::SBMEstimator, dataset::Dataset, sbm::SBM; time_limit::Float64=400.0)::Tuple{OptResults, Matrix{Int}}
    return optimalAssignments(dataset, sbm, time_limit=time_limit, verbose=estimator.verbose)
end

function optimalAssignments(dataset::Dataset, sbm::SBM; time_limit::Float64=400.0, verbose::Bool=false)::Tuple{OptResults, Matrix{Int}}
    """
    Solves the SBM optimization problem with its descriptive formulation,
    given the matrix of probabilities w.

    Parameters
    ----------
    dataset : Dataset
        Dataset representing an observed graph.
    sbm : SBM
        Stochastic Block Model

    Returns
    -------
    (obj_value, elapsed, x_opt) : type
        Value of objective, elapsed time, and optimal values for x
    """
    A = dataset.A
    q = dataset.n_communities
    n = dataset.n
    m = dataset.m
    k = dataset.k
    w = sbm.w

    if verbose
        println("q = $q groups\nn = $n nodes\nm = $m edges");
        println("w=")
        println(w)
        println("Time limit: $time_limit seconds")
    end

    # convert from seconds to milli-seconds
    time_limit = Int(round(time_limit*1000))

    w[w .== 0] .= 1e-30
    W = zeros(n,n,q,q);
    for i=1:n, j=1:n, g=1:q, h=1:q
        if A[i,j] != 0
            W[i,j,g,h] = 0.5*( -A[i,j]*log(w[g,h]) + ((k[i]*k[j])/(2*m))*w[g,h] )
        else
            W[i,j,g,h] = 0.5*((k[i]*k[j])/(2*m))*w[g,h]
        end
    end

    ############################
    # TODO: set time limit and verbose options
    model = Model(solver=GLPKSolverMIP(tm_lim=time_limit)) # version 0.18 of JuMP

    # Variables
    @variable(model, 0 <= y[1:n,1:n,1:q,1:q] <= 1)
    @variable(model, x[1:n,1:q], Bin);

    # Objective
    @objective(model, Min, sum(
            (W[i,j,g,h] * y[i,j,g,h])
            for g=1:q, h=1:q, i=1:n, j=1:n));

    # Constraints
    @constraint(model, con1[i=1:n, j=1:n, g=1:q, h=1:q], x[i,g] - y[i,j,g,h] >= 0)
    @constraint(model, con2[i=1:n, j=1:n, g=1:q, h=1:q], x[j,h] - y[i,j,g,h] >= 0)
    @constraint(model, con3[i=1:n, j=1:n, g=1:q, h=1:q], 1 - x[i,g] - x[j,h] + y[i,j,g,h] >= 0)

    # Assignment constraints
    @constraint(model, assign[i = 1:n], sum(x[i,g] for g=1:q) == 1); # assignment constraint

    # Symmetry breaking constraints
    @constraint(model, conObj1, x[1,1] == 1); # object 1 must be in cluster 1
    @constraint(model, conObj[j=2:n], (1 - x[j,2]) <= sum((1 - x[i,1]) for i=2:(j-1)) + x[j,1] ) # if objects 2,...,j-1 are in cluster 1 and object j is not, then object j must be in cluster 2

    # Solve
    start = time()
    status = solve(model) # version 0.18 of JuMP
    solvetime = time() - start
    obj_lb = getobjectivebound(model);
    if isnan(obj_lb)
        obj_lb = -Inf
    end
    obj_ub = getobjectivevalue(model);
    if isnan(obj_ub)
        obj_ub = Inf
    end
    iterations = nothing
    nodecount = nothing
    lazycount = nothing
    ############################
    # Recover variable values
    x = getvalue(x)
    x = round.(x[:,:])
    x[map(v -> isnan(v) , x)] .= 0 # replace NaN values with 0
    x = Int.(x)

    opt_results = OptResults(obj_lb, obj_ub, status, solvetime, iterations, nodecount, lazycount)
    if verbose
        displayResults(opt_results)
    end
    return opt_results, x
end

function MINLP(estimator::SBMEstimator, dataset::Dataset)
    throw("not implemented")
end


text:
using CUSPARSE
using CUDArt
using Base.Test

m = 25
n = 35
k = 10
blockdim = 5

@testset "cssv" begin
    @testset for elty in [Float32,Float64,Complex64,Complex128]
        @testset "csrsv" begin
            A = rand(elty,m,m)
            A = triu(A)
            x = rand(elty,m)
            alpha = rand(elty)
            d_x = CudaArray(x)
            d_A = CudaSparseMatrixCSR(sparse(A))
            d_y = CUSPARSE.sv('N','T','U',alpha,d_A,d_x,'O')
            h_y = to_host(d_y)
            y = A\(alpha * x)
            @test y ≈ h_y
            x = rand(elty,n)
            d_x = CudaArray(x)
            info = CUSPARSE.sv_analysis('N','T','U',d_A,'O')
            @test_throws DimensionMismatch CUSPARSE.sv_solve('N','U',alpha,d_A,d_x,info,'O')
            A = sparse(rand(elty,m,n))
            d_A = CudaSparseMatrixCSR(A)
            @test_throws DimensionMismatch CUSPARSE.sv_analysis('T','T','U',d_A,'O')
            CUSPARSE.cusparseDestroySolveAnalysisInfo(info)
        end

        @testset "cscsv" begin
            A = rand(elty,m,m)
            A = triu(A)
            x = rand(elty,m)
            alpha = rand(elty)
            d_x = CudaArray(x)
            d_A = CudaSparseMatrixCSC(sparse(A))
            d_y = CUSPARSE.sv('N','T','U',alpha,d_A,d_x,'O')
            h_y = to_host(d_y)
            y = A\(alpha * x)
            @test y ≈ h_y
            x = rand(elty,n)
            d_x = CudaArray(x)
            info = CUSPARSE.sv_analysis('N','T','U',d_A,'O')
            @test_throws DimensionMismatch CUSPARSE.sv_solve('N','U',alpha,d_A,d_x,info,'O')
            A = sparse(rand(elty,m,n))
            d_A = CudaSparseMatrixCSC(A)
            @test_throws DimensionMismatch CUSPARSE.sv_analysis('T','T','U',d_A,'O')
            CUSPARSE.cusparseDestroySolveAnalysisInfo(info)
        end

        @testset "csrsv2" begin
            A = rand(elty,m,m)
            A = triu(A)
            X = rand(elty,m)
            alpha = rand(elty)
            d_X = CudaArray(X)
            d_A = CudaSparseMatrixCSR(sparse(A))
            d_Y = CUSPARSE.sv2('N','U',alpha,d_A,d_X,'O')
            h_Y = to_host(d_Y)
            Y = A\(alpha * X)
            @test Y ≈ h_Y
            d_y = UpperTriangular(d_A)\d_X
            h_y = to_host(d_y)
            y = A\X
            @test y ≈ h_y
            A = sparse(rand(elty,m,n))
            d_A = CudaSparseMatrixCSR(A)
            @test_throws DimensionMismatch CUSPARSE.sv2('N','U',alpha,d_A,d_X,'O')
        end

        @testset "cscsv2" begin
            A = rand(elty,m,m)
            A = triu(A)
            X = rand(elty,m)
            alpha = rand(elty)
            d_X = CudaArray(X)
            d_A = CudaSparseMatrixCSC(sparse(A))
            d_Y = CUSPARSE.sv2('N','U',alpha,d_A,d_X,'O')
            h_Y = to_host(d_Y)
            Y = A\(alpha * X)
            @test Y ≈ h_Y
            d_y = UpperTriangular(d_A)\d_X
            h_y = to_host(d_y)
            y = A\X
            @test y ≈ h_y
            A = sparse(rand(elty,m,n))
            d_A = CudaSparseMatrixCSC(A)
            @test_throws DimensionMismatch CUSPARSE.sv2('N','U',alpha,d_A,d_X,'O')
        end
    end
end


text:
function iup_app(f::Function)
    ccall((:IupOpen, "libiup"), Int32, (Int32, Array{AbstractString}), 0, [])
    f()
    ccall((:IupClose, "libiup"), Void, ())
end

function say_hello()
    iup_app(
        function()
            ccall((:IupMessage, "libiup"), Void, (AbstractString, AbstractString), "Hello World 1", "Hello, World, from IUP/Julia!")
        end
    )
end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: rds_data
using AWS.Compat
using AWS.UUIDs

"""
    batch_execute_statement(resource_arn, secret_arn, sql)
    batch_execute_statement(resource_arn, secret_arn, sql, params::Dict{String,<:Any})

Runs a batch SQL statement over an array of data. You can run bulk update and insert
operations for multiple records using a DML statement with different parameter sets. Bulk
operations can provide a significant performance improvement over individual insert and
update operations.  If a call isn't part of a transaction because it doesn't include the
transactionID parameter, changes that result from the call are committed automatically.

# Arguments
- `resource_arn`: The Amazon Resource Name (ARN) of the Aurora Serverless DB cluster.
- `secret_arn`: The name or ARN of the secret that enables access to the DB cluster.
- `sql`: The SQL statement to run.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"database"`: The name of the database.
- `"parameterSets"`: The parameter set for the batch operation. The SQL statement is
  executed as many times as the number of parameter sets provided. To execute a SQL statement
  with no parameters, use one of the following options:   Specify one or more empty parameter
  sets.   Use the ExecuteStatement operation instead of the BatchExecuteStatement operation.
    Array parameters are not supported.
- `"schema"`: The name of the database schema.
- `"transactionId"`: The identifier of a transaction that was started by using the
  BeginTransaction operation. Specify the transaction ID of the transaction that you want to
  include the SQL statement in. If the SQL statement is not part of a transaction, don't set
  this parameter.
"""
batch_execute_statement(resourceArn, secretArn, sql; aws_config::AbstractAWSConfig=global_aws_config()) = rds_data("POST", "/BatchExecute", Dict{String, Any}("resourceArn"=>resourceArn, "secretArn"=>secretArn, "sql"=>sql); aws_config=aws_config)
batch_execute_statement(resourceArn, secretArn, sql, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = rds_data("POST", "/BatchExecute", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("resourceArn"=>resourceArn, "secretArn"=>secretArn, "sql"=>sql), params)); aws_config=aws_config)

"""
    begin_transaction(resource_arn, secret_arn)
    begin_transaction(resource_arn, secret_arn, params::Dict{String,<:Any})

Starts a SQL transaction.  <important>
A transaction can run for a maximum
of 24 hours. A transaction is terminated and rolled back automatically after 24
hours.

A transaction times out if no calls use its transaction ID in
three minutes. If a transaction times out before it's committed, it's rolled back
automatically.

DDL statements inside a transaction cause an implicit
commit. We recommend that you run each DDL statement in a separate
ExecuteStatement call with
continueAfterTimeout enabled.
 </important>

# Arguments
- `resource_arn`: The Amazon Resource Name (ARN) of the Aurora Serverless DB cluster.
- `secret_arn`: The name or ARN of the secret that enables access to the DB cluster.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"database"`: The name of the database.
- `"schema"`: The name of the database schema.
"""
begin_transaction(resourceArn, secretArn; aws_config::AbstractAWSConfig=global_aws_config()) = rds_data("POST", "/BeginTransaction", Dict{String, Any}("resourceArn"=>resourceArn, "secretArn"=>secretArn); aws_config=aws_config)
begin_transaction(resourceArn, secretArn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = rds_data("POST", "/BeginTransaction", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("resourceArn"=>resourceArn, "secretArn"=>secretArn), params)); aws_config=aws_config)

"""
    commit_transaction(resource_arn, secret_arn, transaction_id)
    commit_transaction(resource_arn, secret_arn, transaction_id, params::Dict{String,<:Any})

Ends a SQL transaction started with the BeginTransaction operation and commits the changes.

# Arguments
- `resource_arn`: The Amazon Resource Name (ARN) of the Aurora Serverless DB cluster.
- `secret_arn`: The name or ARN of the secret that enables access to the DB cluster.
- `transaction_id`: The identifier of the transaction to end and commit.

"""
commit_transaction(resourceArn, secretArn, transactionId; aws_config::AbstractAWSConfig=global_aws_config()) = rds_data("POST", "/CommitTransaction", Dict{String, Any}("resourceArn"=>resourceArn, "secretArn"=>secretArn, "transactionId"=>transactionId); aws_config=aws_config)
commit_transaction(resourceArn, secretArn, transactionId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = rds_data("POST", "/CommitTransaction", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("resourceArn"=>resourceArn, "secretArn"=>secretArn, "transactionId"=>transactionId), params)); aws_config=aws_config)

"""
    execute_sql(aws_secret_store_arn, db_cluster_or_instance_arn, sql_statements)
    execute_sql(aws_secret_store_arn, db_cluster_or_instance_arn, sql_statements, params::Dict{String,<:Any})

Runs one or more SQL statements.  This operation is deprecated. Use the
BatchExecuteStatement or ExecuteStatement operation.

# Arguments
- `aws_secret_store_arn`: The Amazon Resource Name (ARN) of the secret that enables access
  to the DB cluster.
- `db_cluster_or_instance_arn`: The ARN of the Aurora Serverless DB cluster.
- `sql_statements`: One or more SQL statements to run on the DB cluster. You can separate
  SQL statements from each other with a semicolon (;). Any valid SQL statement is permitted,
  including data definition, data manipulation, and commit statements.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"database"`: The name of the database.
- `"schema"`: The name of the database schema.
"""
execute_sql(awsSecretStoreArn, dbClusterOrInstanceArn, sqlStatements; aws_config::AbstractAWSConfig=global_aws_config()) = rds_data("POST", "/ExecuteSql", Dict{String, Any}("awsSecretStoreArn"=>awsSecretStoreArn, "dbClusterOrInstanceArn"=>dbClusterOrInstanceArn, "sqlStatements"=>sqlStatements); aws_config=aws_config)
execute_sql(awsSecretStoreArn, dbClusterOrInstanceArn, sqlStatements, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = rds_data("POST", "/ExecuteSql", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("awsSecretStoreArn"=>awsSecretStoreArn, "dbClusterOrInstanceArn"=>dbClusterOrInstanceArn, "sqlStatements"=>sqlStatements), params)); aws_config=aws_config)

"""
    execute_statement(resource_arn, secret_arn, sql)
    execute_statement(resource_arn, secret_arn, sql, params::Dict{String,<:Any})

Runs a SQL statement against a database.  If a call isn't part of a transaction because it
doesn't include the transactionID parameter, changes that result from the call are
committed automatically.  The response size limit is 1 MB. If the call returns more than 1
MB of response data, the call is terminated.

# Arguments
- `resource_arn`: The Amazon Resource Name (ARN) of the Aurora Serverless DB cluster.
- `secret_arn`: The name or ARN of the secret that enables access to the DB cluster.
- `sql`: The SQL statement to run.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"continueAfterTimeout"`: A value that indicates whether to continue running the
  statement after the call times out. By default, the statement stops running when the call
  times out.  For DDL statements, we recommend continuing to run the statement after the call
  times out. When a DDL statement terminates before it is finished running, it can result in
  errors and possibly corrupted data structures.
- `"database"`: The name of the database.
- `"includeResultMetadata"`: A value that indicates whether to include metadata in the
  results.
- `"parameters"`: The parameters for the SQL statement.  Array parameters are not
  supported.
- `"resultSetOptions"`: Options that control how the result set is returned.
- `"schema"`: The name of the database schema.  Currently, the schema parameter isn't
  supported.
- `"transactionId"`: The identifier of a transaction that was started by using the
  BeginTransaction operation. Specify the transaction ID of the transaction that you want to
  include the SQL statement in. If the SQL statement is not part of a transaction, don't set
  this parameter.
"""
execute_statement(resourceArn, secretArn, sql; aws_config::AbstractAWSConfig=global_aws_config()) = rds_data("POST", "/Execute", Dict{String, Any}("resourceArn"=>resourceArn, "secretArn"=>secretArn, "sql"=>sql); aws_config=aws_config)
execute_statement(resourceArn, secretArn, sql, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = rds_data("POST", "/Execute", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("resourceArn"=>resourceArn, "secretArn"=>secretArn, "sql"=>sql), params)); aws_config=aws_config)

"""
    rollback_transaction(resource_arn, secret_arn, transaction_id)
    rollback_transaction(resource_arn, secret_arn, transaction_id, params::Dict{String,<:Any})

Performs a rollback of a transaction. Rolling back a transaction cancels its changes.

# Arguments
- `resource_arn`: The Amazon Resource Name (ARN) of the Aurora Serverless DB cluster.
- `secret_arn`: The name or ARN of the secret that enables access to the DB cluster.
- `transaction_id`: The identifier of the transaction to roll back.

"""
rollback_transaction(resourceArn, secretArn, transactionId; aws_config::AbstractAWSConfig=global_aws_config()) = rds_data("POST", "/RollbackTransaction", Dict{String, Any}("resourceArn"=>resourceArn, "secretArn"=>secretArn, "transactionId"=>transactionId); aws_config=aws_config)
rollback_transaction(resourceArn, secretArn, transactionId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = rds_data("POST", "/RollbackTransaction", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("resourceArn"=>resourceArn, "secretArn"=>secretArn, "transactionId"=>transactionId), params)); aws_config=aws_config)


text:
# AUTO GENERATED FILE - DO NOT EDIT

export vtk_celldata

"""
    vtk_celldata(;kwargs...)
    vtk_celldata(children::Any;kwargs...)
    vtk_celldata(children_maker::Function;kwargs...)

A CellData component.
CellData is exposing a vtkCellData to a downstream element
Keyword arguments:
- `children` (Array of a list of or a singular dash component, string or numbers | a list of or a singular dash component, string or number; optional)
"""
function vtk_celldata(; kwargs...)
        available_props = Symbol[:children]
        wild_props = Symbol[]
        return Component("vtk_celldata", "CellData", "dash_vtk", available_props, wild_props; kwargs...)
end

vtk_celldata(children::Any; kwargs...) = vtk_celldata(;kwargs..., children = children)
vtk_celldata(children_maker::Function; kwargs...) = vtk_celldata(children_maker(); kwargs...)

max_stars_repo_name: rodrigomologni/dash-vtk

max_issues_repo_name: rodrigomologni/dash-vtk

max_forks_repo_name: rodrigomologni/dash-vtk


text:
@testset "RecursiveVec - singular values full" begin
    @testset for T in (Float32, Float64, ComplexF32, ComplexF64)
        @testset for orth in (cgs2, mgs2, cgsr, mgsr)
            A = rand(T, (n,n))
            v = rand(T, (n,))
            v2 = RecursiveVec((v, zero(v)))
            alg = Lanczos(orth = orth, krylovdim = 2*n, maxiter = 1, tol = precision(T))
            D, V, info = eigsolve(v2, n, :LR, alg) do x
                        x1, x2 = x
                        y1 = A*x2
                        y2 = A'*x1
                        return RecursiveVec(y1, y2)
            end
            @test info.converged >= n
            S = D[1:n]
            @test S ≈ svdvals(A)
            UU = hcat((sqrt(2*one(T))*v[1] for v in V[1:n])...)
            VV = hcat((sqrt(2*one(T))*v[2] for v in V[1:n])...)
            @test UU*Diagonal(S)*VV' ≈ A
        end
    end
end

@testset "RecursiveVec - singular values iteratively" begin
    @testset for T in (Float32, Float64, ComplexF32, ComplexF64)
        @testset for orth in (cgs2, mgs2, cgsr, mgsr)
            A = rand(T, (N,2*N))
            v = rand(T, (N,))
            w = rand(T, (2*N,))
            v2 = RecursiveVec((v, w))
            alg = Lanczos(orth = orth, krylovdim = n, maxiter = 300, tol = precision(T))
            n1 = div(n,2)
            D, V, info = eigsolve(v2, n1, :LR, alg) do x
                        x1, x2 = x
                        y1 = A*x2
                        y2 = A'*x1
                        return RecursiveVec(y1, y2)
            end
            @test info.converged >= n1
            S = D[1:n1]
            @test S ≈ svdvals(A)[1:n1]
            UU = hcat((sqrt(2*one(T))*v[1] for v in V[1:n1])...)
            VV = hcat((sqrt(2*one(T))*v[2] for v in V[1:n1])...)
            @test Diagonal(S) ≈ UU'*A*VV
        end
    end
end


text:
#=
NodeWiseTriangleCount:
- Julia version: 1.6.2
- Author: samuel
- Date: 2021-09-10
=#

using SuiteSparseGraphBLAS
using SparseArrays
using LinearAlgebra

#A is the input matrix. Output: triangles originated by each node
function nodewise_tc(A)

    #extracting lower matrix from input matrix
    partial = mul(A, A, Semirings.PLUS_TIMES, mask=A)
    t = reduce(Monoids.PLUS_MONOID[Float64], partial, dims=2)
    print(t/2)
    return t/2

end


text:
using Literate

# * If true, the code will be executed by Literate and the output will be captured.
# * If false, the code will be executed by Documenter as an example block.
files = Dict(
    "example1.jl" => true,
    "example2.jl" => false,
    "example3.jl" => true,
    "example4.jl" => true
)

function lit_to_md(file; execute=false)
    examples_path = joinpath("src", "examples")
    out_md_path = "src/"
    Literate.markdown(
        joinpath(examples_path, file),
        out_md_path;
        documenter=true,
        execute=execute
    )

    return nothing
end

function parse_files(some_files)
    for (k, v) in some_files
        lit_to_md(k; execute=v)
    end

    return nothing
end

parse_files(files)

max_stars_repo_name: gmagannaDevelop/Elysivm

max_issues_repo_name: gmagannaDevelop/Elysivm

max_forks_repo_name: gmagannaDevelop/Elysivm


text:
mutable struct LogEnvFieldsEveryNStep <: AbstractHook
    n::Int
    t::Int
    fields::Vector{Symbol}
    lg::TBLogger
end
function LogEnvFieldsEveryNStep(save_dir, fields; n=1, t=0)
    lg = TBLogger(joinpath(save_dir, "tb_log"), min_level = Logging.Info)
    LogEnvFieldsEveryNStep(n, t, fields, lg)
end
function (hook::LogEnvFieldsEveryNStep)(::PostActStage, agent, env)
    hook.t += 1
    hook.t % hook.n != 0 && return

    while typeof(env) <: AbstractEnvWrapper
        env = env.env
    end

    with_logger(hook.lg) do
        for field in hook.fields
            @info "env_fields" field = getfield(env, field) log_step_increment = 0
        end
        @info "increment" log_step_increment=hook.n
    end
end

# Helper methods to get to base env for logging
Base.getindex(env::MultiThreadEnv, ::typeof(!)) = env[1][!]

mutable struct LogEveryNStep <: AbstractHook
    n::Int
    t::Int
    lg::TBLogger{String, IOStream}
end
function LogEveryNStep(save_dir::String; n=1)
    lg = TBLogger(joinpath(save_dir, "tb_log"), min_level = Logging.Info)
    LogEveryNStep(n, 0, lg)
end
function (hook::LogEveryNStep)(::PostActStage, agent, env)
    hook.t += 1
    hook.t % hook.n != 0 && return

    with_logger(hook.lg) do
        log_agent(agent)
        log_env(env[!])
        @info "increment" log_step_increment=hook.n
    end
end

log_agent(::AbstractPolicy) = nothing
log_agent(a::Agent{<:AbstractPolicy}) = log_agent(a.policy)
function log_agent(policy::PPOPolicy)
    @info(
        "training",
        actor_loss = policy.actor_loss[end, end],
        critic_loss = policy.critic_loss[end, end],
        entropy_loss = policy.entropy_loss[end, end],
        total_loss = policy.loss[end, end],
        log_step_increment = 0
    )
end
function log_agent(policy::SACPolicy)
    @info(
        "training",
        reward_term = policy.reward_term,
        entropy_term = policy.entropy_term,
        alpha = policy.α,
        log_step_increment = 0
    )
end
function log_agent(agent::GuidedPolicy)
    @info(
        "training",
        learner_action = agent.actor_history[agent.history_ptr],
        log_step_increment = 0
    )
    log_agent(agent.learner)
    log_agent(agent.teacher)
end
function log_env(env::ServiceMeshEnv)
    for (i, ms) in enumerate(env.microservices)
        @info(
            "MS$(i)",
            scale = ms.running_nodes,
            queue = max(0, length(ms.queue) - ms.running_nodes),
            booting = env.booting,
            closing = env.closing,
            processed_jobs = env.processed_jobs[i],
            utilization = ms.running_nodes == 0 ? isempty(ms.queue) ? 1 : ms.max_queue : length(ms.queue) / ms.running_nodes,
            log_step_increment = 0
        )
    end
    for i in eachindex(env.jobtypes)
        @info(
            "Job$(i)",
            arrivals = env.arrivals[i],
            dropped = env.dropped_jobs[i],
            finished = env.finished_jobs[i],
            missed = env.missed_deadlines[i],
            log_step_increment = 0
        )
    end
    @info(
        "other",
        time = env.time,
        log_step_increment = 0
    )
    @info(
        "training",
        profit = reward(env[!]), # Reward of base env is pure profit, other stuff is added in wrappers
        reward = reward(env),
        log_step_increment = 0
    )
end

"""
    Keep track of a running average of the reward over a number of steps `n`.
"""
mutable struct RunningReward{F} <: AbstractHook
    rewards::Vector{Float64}
    pointer::Int
    envmap::F
end
function RunningReward(n; envmap=identity)
    RunningReward(zeros(n), 1, envmap)
end
function (hook::RunningReward)(::PostActStage, agent, env)
    hook.rewards[hook.pointer] = mean(reward(hook.envmap(env))) # mean in case of MultiThreadEnv
    hook.pointer %= length(hook.rewards)
    hook.pointer += 1
end

max_stars_repo_name: albheim/ServiceMeshControl

max_issues_repo_name: albheim/ServiceMeshControl

max_forks_repo_name: albheim/ServiceMeshControl


text:
module ObjectDetectionStats
    using DataFrames
    include("LocalizationMeasures.jl")
    export Box, ==, barycenterbox, LLURtoULLR, translate, translate!,
        clamp, clamp!, area, intersection_area, intersection_over_union

    include("LocalizationStatistics.jl")
    export HotClassLocalization, ColdClassLocalization, ObjectDetectionScore,
        add!, micro_precision, micro_recall, micro_F1, macro_precision,
        macro_recall, macro_F1, classwise_PR_estimate
end # module


text:
"""
    Creates `RDAIO` wrapper for `io` stream depending on its format
    specified by `formatcode`.
"""
function rdaio(io::IO, formatcode::AbstractString)
    if formatcode == "X" XDRIO(io)
    elseif formatcode == "A" ASCIIIO(io)
    elseif formatcode == "B" NativeIO(io)
    else error("Unrecognized RDA format \"$formatcode\"")
    end
end

if LONG_VECTOR_SUPPORT
    # reads the length of any data vector from a stream
    # from R's serialize.c
    function readlength(io::RDAIO)
        len = convert(RVecLength, readint32(io))
        if (len < -1) error("negative serialized length for vector")
        elseif (len >= 0)
            return len
        else # big vectors, the next 2 ints encode the length
            len1, len2 = convert(RVecLength, readint32(io)), convert(RVecLength, readint32(io))
            # sanity check for now
            if (len1 > 65536) error("invalid upper part of serialized vector length") end
            return (len1 << 32) + len2
        end
    end
else
    # reads the length of any data vector from a stream
    # fails when long (> 2^31-1) vector encountered
    # from R's serialize.c
    function readlength(io::RDAIO)
        len = convert(RVecLength, readint32(io))
        if (len >= 0)
            return len
        elseif (len < -1)
            error("negative serialized length for vector")
        else
            error("negative serialized vector length:\nperhaps long vector from 64-bit version of R?")
        end
    end
end

immutable CHARSXProps # RDA CHARSXP properties
  levs::UInt32       # level flags (encoding etc) TODO process
  nchar::Int32       # string length, -1 for NA strings
end

function readcharsxprops(io::RDAIO) # read character string encoding and length
    fl = readuint32(io)
    @assert sxtype(fl) == CHARSXP
    @assert !hasattr(fl)
### watch out for levs in here.  Generally it has the value 0x40 so that fl = 0x00040009 (262153)
### if levs == 0x00 then the next line should be -1 to indicate the NA_STRING
    CHARSXProps(fl >> 12, readint32(io))
end

function readcharacter(io::RDAIO)  # a single character string
    props = readcharsxprops(io)
    props.nchar==-1 ? "" : readnchars(io, props.nchar)
end

function readcharacter(io::RDAIO, n::RVecLength)  # a single character string
    res = fill("", n)
    na = falses(n)
    for i in 1:n
        props = readcharsxprops(io)
        if (props.nchar==-1) na[i] = true
        else res[i] = readnchars(io, props.nchar)
        end
    end
    return res, na
end


text:
using Test
@testset "Multigrid" begin
    include("Multigrid/runtests.jl")
	include("ParallelJuliaSolver/runtests.jl")
	include("DomainDecomposition/runtests.jl")
end


text:
#=
Abstraction over the GPU thread indexing functions.
Uses CUDA like names
=#
for sym in (:x, :y, :z)
    for f in (:blockidx, :blockdim, :threadidx, :griddim)
        fname = Symbol(string(f, '_', sym))
        @eval $fname(state)::UInt32 = error("Not implemented")
        @eval export $fname
    end
end

"""
     synchronize_threads(state)

in CUDA terms `__synchronize`
in OpenCL terms: `barrier(CLK_LOCAL_MEM_FENCE)`
"""
function synchronize_threads(state)
    error("Not implemented")
end

"""
    linear_index(state)

linear index corresponding to each kernel launch (in OpenCL equal to get_global_id).

"""
@inline function linear_index(state)
    UInt32((blockidx_x(state) - UInt32(1)) * blockdim_x(state) + threadidx_x(state))
end

"""
    linearidx(A, statesym = :state)

Macro form of `linear_index`, which calls return when out of bounds.
So it can be used like this:

    ```julia
    function kernel(state, A)
        idx = @linear_index A state
        # from here on it's save to index into A with idx
        @inbounds begin
            A[idx] = ...
        end
    end
    ```
"""
macro linearidx(A, statesym = :state)
    quote
        x1 = $(esc(A))
        i1 = linear_index($(esc(statesym)))
        i1 > length(x1) && return
        i1
    end
end

"""
    cartesianidx(A, statesym = :state)

Like [`@linearidx(A, statesym = :state)`](@ref), but returns an N-dimensional `NTuple{ndim(A), Cuint}` as index
"""
macro cartesianidx(A, statesym = :state)
    quote
        x = $(esc(A))
        i2 = @linearidx(x, $(esc(statesym)))
        gpu_ind2sub(x, i2)
    end
end

"""
    global_size(state)

Global size == blockdim * griddim == total number of kernel execution
"""
@inline function global_size(state)
    # TODO nd version
    griddim_x(state) * blockdim_x(state)
end

"""
    device(A::AbstractArray)

Gets the device associated to the Array `A`
"""
function device(A::AbstractArray)
    # fallback is a noop, for backends not needing synchronization. This
    # makes it easier to write generic code that also works for AbstractArrays
end

"""
    synchronize(A::AbstractArray)

Blocks until all operations are finished on `A`
"""
function synchronize(A::AbstractArray)
    # fallback is a noop, for backends not needing synchronization. This
    # makes it easier to write generic code that also works for AbstractArrays
end
#
# @inline function synchronize_threads(state)
#     CUDAnative.__syncthreads()
# end

"""
    gpu_call(kernel::Function, A::GPUArray, args::Tuple, configuration = length(A))

Calls function `kernel` on the GPU.
`A` must be an GPUArray and will help to dispatch to the correct GPU backend
and supplies queues and contexts.
Calls the kernel function with `kernel(state, args...)`, where state is dependant on the backend
and can be used for getting an index into `A` with `linear_index(state)`.
Optionally, a launch configuration can be supplied in the following way:

    1) A single integer, indicating how many work items (total number of threads) you want to launch.
        in this case `linear_index(state)` will be a number in the range `1:configuration`
    2) Pass a tuple of integer tuples to define blocks and threads per blocks!

"""
function gpu_call(kernel, A::GPUArray, args::Tuple, configuration = length(A))
    ITuple = NTuple{N, Integer} where N
    # If is a single integer, we assume it to be the global size / total number of threads one wants to launch
    thread_blocks = if isa(configuration, Integer)
        thread_blocks_heuristic(configuration)
    elseif isa(configuration, ITuple)
        # if a single integer ntuple, we assume it to configure the blocks
        configuration,  ntuple(x-> x == 1 ? 256 : 1, length(configuration))
    elseif isa(configuration, Tuple{ITuple, ITuple})
        # 2 dim tuple of ints == blocks + threads per block
        if any(x-> length(x) > 3 || length(x) < 1, configuration)
            error("blocks & threads must be 1-3 dimensional. Found: $configuration")
        end
        map(x-> Int.(x), configuration) # make sure it all has the same int type
    else
        error("""Please launch a gpu kernel with a valid configuration.
            Found: $configurations
            Configuration needs to be:
            1) A single integer, indicating how many work items (total number of threads) you want to launch.
                in this case `linear_index(state)` will be a number in the range 1:configuration
            2) Pass a tuple of integer tuples to define blocks and threads per blocks!
                `linear_index` will be inbetween 1:prod((blocks..., threads...))
        """)
    end
    _gpu_call(kernel, A, args, thread_blocks)
end

# Internal GPU call function, that needs to be overloaded by the backends.
_gpu_call(f, A, args, thread_blocks) = error("Not implemented")


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

mutable struct IoK8sApiCoreV1ConfigMapVolumeSource <: SwaggerModel
    defaultMode::Any # spec type: Union{ Nothing, Int32 } # spec name: defaultMode
    items::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1KeyToPath} } # spec name: items
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    optional::Any # spec type: Union{ Nothing, Bool } # spec name: optional

    function IoK8sApiCoreV1ConfigMapVolumeSource(;defaultMode=nothing, items=nothing, name=nothing, optional=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1ConfigMapVolumeSource, Symbol("defaultMode"), defaultMode)
        setfield!(o, Symbol("defaultMode"), defaultMode)
        validate_property(IoK8sApiCoreV1ConfigMapVolumeSource, Symbol("items"), items)
        setfield!(o, Symbol("items"), items)
        validate_property(IoK8sApiCoreV1ConfigMapVolumeSource, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApiCoreV1ConfigMapVolumeSource, Symbol("optional"), optional)
        setfield!(o, Symbol("optional"), optional)
        o
    end
end # type IoK8sApiCoreV1ConfigMapVolumeSource

const _property_map_IoK8sApiCoreV1ConfigMapVolumeSource = Dict{Symbol,Symbol}(Symbol("defaultMode")=>Symbol("defaultMode"), Symbol("items")=>Symbol("items"), Symbol("name")=>Symbol("name"), Symbol("optional")=>Symbol("optional"))
const _property_types_IoK8sApiCoreV1ConfigMapVolumeSource = Dict{Symbol,String}(Symbol("defaultMode")=>"Int32", Symbol("items")=>"Vector{IoK8sApiCoreV1KeyToPath}", Symbol("name")=>"String", Symbol("optional")=>"Bool")
Base.propertynames(::Type{ IoK8sApiCoreV1ConfigMapVolumeSource }) = collect(keys(_property_map_IoK8sApiCoreV1ConfigMapVolumeSource))
Swagger.property_type(::Type{ IoK8sApiCoreV1ConfigMapVolumeSource }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiCoreV1ConfigMapVolumeSource[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1ConfigMapVolumeSource }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1ConfigMapVolumeSource[property_name]

function check_required(o::IoK8sApiCoreV1ConfigMapVolumeSource)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1ConfigMapVolumeSource }, name::Symbol, val)
end


text:
using DrWatson, DataFrames, Plots, Printf, Parameters

include(srcdir("load.jl"))
include(srcdir("binning.jl"))

function binnedgif(df)
    @unpack drug, phase, trial, gw, gh, videopath = df

    basepath = projectdir("plots", drug, string(trial), @sprintf "%dx%d" gh gw)
    mkpath(basepath)

    greenchannel = datadir(videopath) |> frames |> fs -> crop(fs, 400, 400, :) |> green
    greenbinned = coarse(greenchannel, gh, gw) |> rescale |> bin

    anim = @animate for t in 1:size(greenbinned, 3)
        title = @sprintf "Binned %s %s, Trial %d; Frame: %d" titlecase(phase) drug trial t
        heatmap(greenbinned[:,:,t] .- 1,
                title = title,
                titleloc = :left,
                clim = (0,1))
    end

    giffile = joinpath(basepath, phase * "-binned.gif")
    gif(anim, giffile, fps=10)
end

function coarsegif(df)
    @unpack drug, phase, trial, gw, gh, videopath = df

    basepath = projectdir("plots", drug, string(trial), @sprintf "%dx%d" gh gw)
    mkpath(basepath)

    greenchannel = datadir(videopath) |> frames |> fs -> crop(fs, 400, 400, :) |> green
    greengrid = coarse(greenchannel, gh, gw) |> rescale

    anim = @animate for t in 1:size(greengrid, 3)
        title = @sprintf "%s %s, Trial %d; Frame: %d" titlecase(phase) drug trial t
        heatmap(greengrid[:,:,t],
                title = title,
                titleloc = :left,
                clim = (0,1))
    end

    giffile = joinpath(basepath, phase * ".gif")
    gif(anim, giffile, fps=10)
end

function migif(df)
    simdata = first(unique(select(df, :drug, :phase, :trial, :gw, :gh, :nperms)))
    @unpack drug, phase, trial, gw, gh, nperms = simdata

    grid = @sprintf "%dx%d" gh gw
    basepath = projectdir("plots", drug, string(trial), grid, string(nperms))
    mkpath(basepath)
    mkpath(joinpath(basepath, "frames", "mi"))

    mis = sort(select(df, :lag, :mi), :lag)

    miframes = []
    anim = @animate for row in eachrow(mis)
        @unpack lag, mi = row

        title = if lag == 0
            @sprintf "Mutual Info %s %s" titlecase(phase) drug
        elseif lag == 1
            @sprintf "Mutual Info %s %s; Lagged %d Frame" titlecase(phase) drug lag
        else
            @sprintf "Mutual Info %s %s; Lagged %d Frames" titlecase(phase) drug lag
        end

        p = heatmap(linearize(mi),
                    title=title,
                    titleloc=:left,
                    xlabel="target cell",
                    ylabel="source cell",
                    clim=(0,1))
        push!(miframes, p)
        p
    end

    for (i, p) in enumerate(miframes)
        framename = @sprintf "%d-%s.png" (i-1) phase
        framefile = joinpath(basepath, "frames", "mi", framename)
        savefig(p, framefile)
    end

    giffile = joinpath(basepath, phase * "-mi.gif")
    gif(anim, giffile, fps=0.5)
end

function tegif(df)
    simdata = first(unique(select(df, :drug, :phase, :trial, :gw, :gh, :nperms)))
    @unpack drug, phase, trial, gw, gh, nperms = simdata

    grid = @sprintf "%dx%d" gh gw
    basepath = projectdir("plots", drug, string(trial), grid, string(nperms))
    mkpath(basepath)
    mkpath(joinpath(basepath, "frames", "te"))

    tes = sort(select(df, :k, :te), :k)

    teframes = []
    anim = @animate for row in eachrow(tes)
        @unpack k, te = row

        title = @sprintf "Transfer Entropy %s %s; k=%d" titlecase(phase) drug k

        p = heatmap(linearize(te),
                    title=title,
                    titleloc=:left,
                    xlabel="target cell",
                    ylabel="source cell",
                    clim=(0,1))
        push!(teframes, p)
        p
    end

    for (i, p) in enumerate(teframes)
        framename = @sprintf "%d-%s.png" i phase
        framefile = joinpath(basepath, "frames", "te", framename)
        savefig(p, framefile)
    end

    giffile = joinpath(basepath, phase * "-te.gif")
    gif(anim, giffile, fps=0.5)
end

function main()
    df = collect_results(datadir("info", "mi"))
    for group in groupby(df, [:drug, :trial, :phase, :gh, :gw])
        migif(group)
    end
    df = collect_results(datadir("info", "te"))
    for group in groupby(df, [:drug, :trial, :phase, :gh, :gw])
        tegif(group)
    end
    for group in eachrow(unique(select(df, :drug, :trial, :phase, :gh, :gw, :videopath)))
        coarsegif(group)
        binnedgif(group)
    end
end

@time main()

max_stars_repo_name: ELIFE-ASU/dglmoore

max_issues_repo_name: ELIFE-ASU/dglmoore

max_forks_repo_name: ELIFE-ASU/dglmoore