text:
"Above the center/spawn sensor"
function inspawnsensorx(x::Float64)::Bool
    return SENSOR_SPAWN_LEFT_X < x < SENSOR_SPAWN_RIGHT_X
end

function insecondsensorx(x::Float64)::Bool
    return SENSOR_SECOND_LEFT_X < x < SENSOR_SECOND_RIGHT_X
end

function inthirdsensorx(x::Float64)::Bool
    return SENSOR_THIRD_LEFT_X < x < SENSOR_THIRD_RIGHT_X
end

function infourthsensorx(x::Float64)::Bool
    return SENSOR_FOURTH_LEFT_X < x < SENSOR_FOURTH_RIGHT_X
end

function insensor(xy::Array{Float64, 1})
    x, y = xy
    withinx = inspawnsensorx(x) || insecondsensorx(x) || inthirdsensorx(x) || infourthsensorx(x)
    withiny = y <= WALL_Y
    return withinx && withiny
end

"Parylene walls in between sensors"
function inwalls(x::Float64, y::Float64)::Bool
    withinx = !inspawnsensorx(x) && !insecondsensorx(x) && !inthirdsensorx(x) && !infourthsensorx(x)
    withiny = y <= WALL_Y
    return withinx && withiny
end

"Designate zones to simplify identical four sensors into easier subproblems"
function inzoneone(x::Float64)::Bool
    return -1 * ESCAPE_X <= x <= ZONE_ONE_RIGHT_X
end

function inzonetwo(x::Float64)::Bool
    return ZONE_ONE_RIGHT_X < x <= ZONE_TWO_RIGHT_X
end

function inzonethree(x::Float64)::Bool
    return ZONE_TWO_RIGHT_X < x <= ZONE_THREE_RIGHT_X
end

function inzonefour(x::Float64)::Bool
    return ZONE_THREE_RIGHT_X < x <= ESCAPE_X
end

function adjustxy(xy::Array{Float64,1})
    copyxy = copy(xy)
    x = xy[1]
    if inzoneone(x)
        return copyxy, "no adjust"
    elseif inzonetwo(x)
        copyxy[1] -= SENSOR_SECOND_LEFT_X
        return copyxy, "second"
    elseif inzonethree(x)
        copyxy[1] -= SENSOR_THIRD_LEFT_X
        return copyxy, "third"
    elseif inzonefour(x)
        copyxy[1] -= SENSOR_FOURTH_LEFT_X
        return copyxy, "fourth"
    else
        println("adjust error")
    end
end

function revertxy(xy::Array{Float64}, whichsensor::String)
    copyxy = copy(xy)
    if whichsensor == "no adjust"
        return xy
    elseif whichsensor == "second"
        copyxy[1] += SENSOR_SECOND_LEFT_X
        return copyxy
    elseif whichsensor == "third"
        copyxy[1] += SENSOR_THIRD_LEFT_X
        return copyxy
    elseif whichsensor == "fourth"
        copyxy[1] += SENSOR_FOURTH_LEFT_X
        return copyxy
    else
        println("revert error")
    end
end

"Assume inputs already adjusted"
function inwater(adjustedxy::Array{Float64,1})::Bool
    x, y = adjustedxy
    left_enz = x < ENZYME_LEFT_X
    right_enz = x > ENZYME_RIGHT_X
    above_enz = y > ENZYME_MAX_Y
    return left_enz || right_enz || above_enz
end

function inenz(adjustedxy::Array{Float64,1})::Bool
    x, y = adjustedxy
    left_ppd = ENZYME_LEFT_X <= x < PPD_LEFT_X && WALL_Y < y <= ENZYME_MAX_Y
    right_ppd = PPD_RIGHT_X < x <= ENZYME_RIGHT_X && WALL_Y < y <= ENZYME_MAX_Y
    above_ppd = ENZYME_LEFT_X <= x <= ENZYME_RIGHT_X && PPD_MAX_Y < y <= ENZYME_MAX_Y
    return left_ppd || right_ppd || above_ppd
end

function inppd(adjustedxy::Array{Float64, 1})::Bool
    x, y = adjustedxy
    withinx = PPD_LEFT_X <= x <= PPD_RIGHT_X
    withiny = PPD_MIN_Y < y <= PPD_MAX_Y
    return withinx && withiny
end

"Assumed no flow in Quinto"
function flow_arlett(xy::Array{Float64,1})
    return 0.0
end

"""
    spawnrandompoint(ranfloat, randfloat) -> Array{Float64,1}

Initialize an [x,y] coordinate, uniformly random at water-enz interface, just inside enz.
"""
function spawnrandompoint(randfloat::Float64=rand(), randfloatstub::Float64=rand())
    spawncoordinate = SPAWN_LENGTH * randfloat + ENZYME_LEFT_X - SPAWN_ENZYME_MAX_Y
    if SPAWN_LEFT_X <= spawncoordinate <= SPAWN_RIGHT_X
        return [spawncoordinate, SPAWN_ENZYME_MAX_Y]
    elseif randfloat <= 0.5
        return [SPAWN_LEFT_X, SPAWN_LENGTH * randfloat]
    else
        return [SPAWN_RIGHT_X, SPAWN_LENGTH * (1 - randfloat)]
    end
end

"""
    whereoutsideenz(xy)

Assign a cardinal direction (W, NW, N, NE, or E) when crossing the water-enzyme boundary
"""
function whereoutsideenz(xy::Array{Float64,1})
    x, y = xy
    if x < ENZYME_LEFT_X
        if y > ENZYME_MAX_Y
            return "NW"
        else
            return "W"
        end
    elseif x > ENZYME_RIGHT_X
        if y > ENZYME_MAX_Y
            return "NE"
        else
            return "E"
        end
    elseif y > ENZYME_MAX_Y
        return "N"
    else
        return "whereoutsideenz Error"
    end
end

"""
    whereoutsideppd(xy)

Assign a cardinal direction (W, NW, N, NE, or E) when crossing the enzyme-ppd boundary
"""
function whereoutsideppd(xy::Array{Float64,1})
    x, y = xy
    if x < PPD_LEFT_X
        if y > PPD_MAX_Y
            return "NW"
        else
            return "W"
        end
    elseif x > PPD_RIGHT_X
        if y > PPD_MAX_Y
            return "NE"
        else
            return "E"
        end
    elseif y > PPD_MAX_Y
        return "N"
    else
        return "whereoutsideppd Error"
    end
end

max_stars_repo_name: s-y-wu/HMCResearchRandomWalks

max_issues_repo_name: s-y-wu/HMCResearchRandomWalks

max_forks_repo_name: s-y-wu/HMCResearchRandomWalks


text:
const Accessor = Function;

abstract type Averager end;

# TODO: consider creating an averager which uses arbitrary precision arithmetic
mutable struct StandardAverager{T,N} <: Averager
  value::T;
  normalizer::N;
  accessor::Accessor;
end

function StandardAverager(accessor::Accessor, chain::EAPChain)
  return StandardAverager(
                          0.0 * accessor(chain),
                          0.0,
                          accessor
                         );
end

function StandardAverager(dt::DataType, accessor::Accessor, chain::EAPChain)
  return StandardAverager(
                          dt(0.0 * accessor(chain)),
                          dt(0.0),
                          accessor
                         );
end

function StandardAverager(dt1::DataType, dt2::DataType, accessor::Accessor, chain::EAPChain)
  return StandardAverager(
                          dt2(0.0 * accessor(chain)),
                          dt1(0.0),
                          accessor
                         );
end

@inline get_avg(stdavg::StandardAverager) = stdavg.value / stdavg.normalizer;

function record!(stdavg::StandardAverager, chain::EAPChain)
  stdavg.value += stdavg.accessor(chain);
  stdavg.normalizer += 1;
end

function record!(stdavg::StandardAverager{<:Vector}, chain::EAPChain)
  stdavg.value[:] += stdavg.accessor(chain);
  stdavg.normalizer += 1;
end

abstract type WeightFunction end;

struct UmbrellaAverager{T,N} <: Averager
  stdavg::StandardAverager{T,N};
  weightFunction::WeightFunction;
end

function UmbrellaAverager(ac::Accessor, wf::WeightFunction, chain::EAPChain)
  return UmbrellaAverager(StandardAverager(ac, chain), wf);
end

@inline get_avg(ua::UmbrellaAverager) = get_avg(ua.stdavg);

function record!(ua::UmbrellaAverager, chain::EAPChain)
  expw = exp(ua.weightFunction(chain));
  ua.stdavg.value += ua.stdavg.accessor(chain) / expw;
  ua.stdavg.normalizer += 1.0 / expw;
end

function record!(ua::UmbrellaAverager{Vector{<:Real}}, chain::EAPChain)
  expw = exp(ua.weightFunction(chain));
  ua.stdavg.value[:] += ua.stdavg.accessor(chain) / expw;
  ua.stdavg.normalizer += 1.0 / expw;
end

function record!(ua::UmbrellaAverager{Dec128}, chain::EAPChain)
  expw = exp(Dec128(ua.weightFunction(chain)));
  ua.stdavg.value += ua.stdavg.accessor(chain) / expw;
  ua.stdavg.normalizer += 1.0 / expw;
end

function record!(ua::UmbrellaAverager{Vector{Dec128}}, chain::EAPChain)
  expw = exp(Dec128(ua.weightFunction(chain)));
  ua.stdavg.value[:] += ua.stdavg.accessor(chain) / expw;
  ua.stdavg.normalizer += 1.0 / expw;
end

function record!(ua::UmbrellaAverager{BigFloat}, chain::EAPChain)
  expw = exp(BigFloat(ua.weightFunction(chain)));
  ua.stdavg.value += ua.stdavg.accessor(chain) / expw;
  ua.stdavg.normalizer += 1.0 / expw;
end

function record!(ua::UmbrellaAverager{Vector{BigFloat}}, chain::EAPChain)
  expw = exp(BigFloat(ua.weightFunction(chain)));
  ua.stdavg.value[:] += ua.stdavg.accessor(chain) / expw;
  ua.stdavg.normalizer += 1.0 / expw;
end

struct WeightlessFunction <: WeightFunction
end

(wf::WeightlessFunction)(::EAPChain) = 1.0;

struct AntiDipoleWeightFunction <: WeightFunction
  log_gauge::Float64;
end

# forward call
@inline AntiDipoleWeightFunction(chain::EAPChain) = AntiDipoleWeightFunction(chain, chain.μ);
function AntiDipoleWeightFunction(chain::EAPChain, dr::DielectricResponse)
  return AntiDipoleWeightFunction(-(dr.K1 + 2 * dr.K2) * chain.E0*chain.E0 *
                                   n(chain) / (3 * chain.kT) + chain.Ω);
end

function AntiDipoleWeightFunction(chain::EAPChain, dr::PolarResponse)
  return AntiDipoleWeightFunction(-eigvals(dr.M)[end] * chain.E0 *
                                   n(chain) / (3 * chain.kT) + chain.Ω);
end

function (wf::AntiDipoleWeightFunction)(chain::EAPChain)
  return (sum(chain.us) / chain.kT *
          (0.2 + 0.8*(exp(-(chain.Fx^2+chain.Fz^2)/chain.kT)))
          - wf.log_gauge);
end

ext: jl

max_stars_repo_name: grasingerm/polymer-stats

max_issues_repo_name: grasingerm/polymer-stats

max_forks_repo_name: grasingerm/polymer-stats


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_semiinfinite/sub_problems.jl
# Defines utilities for generic SIP subroutines.
#############################################################################

"""
    build_model

Create the model and variables used with extension `t::EAGO.ExtensionType` in algorithm
`a::AbstractSIPAlgo` in subproblem `s::AbstractSubproblemType` via the command
`build_model(t::ExtensionType, a::AbstractSIPAlgo, s::AbstractSubproblemType, p::SIPProblem)`.
"""
function build_model(t::DefaultExt, a::A, s::S, p::SIPProblem) where {A <: AbstractSIPAlgo, S <: AbstractSubproblemType}
    model = Model(get_sip_optimizer(t,a,s))
    for (k,v) in get_sip_kwargs(s,p)
        MOI.set(model, MOI.RawParameter(String(k)), v)
    end
    vL, vU, nv = get_bnds(s,p)
    @variable(model, vL[i] <= v[i=1:nv] <= vU[i])
    return model, v
end
function build_model(t::ExtensionType, a::A, s::S, p::SIPProblem) where {A <: AbstractSIPAlgo, S <: AbstractSubproblemType}
    build_model(DefaultExt(), a, s, p)
end

###
### General set tolerances
###
function set_tolerance_inner!(t::DefaultExt, alg, s, m::JuMP.Model, abs_tol::Float64)
    optimizer_name = JuMP.solver_name(m)
    if optimizer_name === "EAGO: Easy Advanced Global Optimization"
        set_optimizer_attribute(m, "absolute_tolerance", abs_tol)
        #set_optimizer_attribute(m, "constr_viol_tol", c_tol)
        #set_optimizer_attribute(m, "relative_tolerance", rel_tol)
    elseif optimizer_name === "SCIP"
        set_optimizer_attribute(m, "limits/absgap", abs_tol)
    elseif optimizer_name === "Alpine"
        set_optimizer_attribute(m, "absgap", abs_tol)
    elseif optimizer_name === "BARON"
        set_optimizer_attribute(m, "EpsA", abs_tol)
    elseif optimizer_name === "GAMS"
        set_optimizer_attribute(m, "OptCA", abs_tol)
    else
        error("A custom set_tolerance! function for solver = $optimizer_name
               specified for use with subproblem = $s in algorithm = $alg with
               extension = $t has not been defined and the selected solver is
               not support by default. Please open an issue requesting this feature
               at the following link https://github.com/PSORLab/EAGO.jl/issues.
               Extending the EAGO.set_tolerance! with a custom extension type
               will resolve this issue.")
    end
    return nothing
end
function set_tolerance!(t::DefaultExt, alg::AbstractSIPAlgo, s::S, m::JuMP.Model,
                        sr::SIPSubResult, i::Int) where {S <: AbstractSubproblemType}
    return nothing
end
function set_tolerance!(t::ExtensionType, alg::A, s::S, m::JuMP.Model,
                        sr::SIPSubResult, i::Int) where {A <: AbstractSIPAlgo,
                                                         S <: AbstractSubproblemType}
    set_tolerance!(t, alg, s, m, sr, i)
end

###
### General get discretation set
###
function get_disc_set(t::DefaultExt, alg::AbstractSIPAlgo, s::S, sr::SIPProblem, i::Int) where {S <: AbstractSubproblemType}
    return Vector{Float64}[]
end
function get_disc_set(t::ExtensionType, alg::AbstractSIPAlgo, s::S, sr::SIPProblem, i::Int) where {S <: AbstractSubproblemType}
    get_disc_set(DefaultExt(), alg, s, p, i)
end

###
### Add constraints solely on uncertainty
###
function add_uncertainty_constraint!(m::JuMP.Model, prob::SIPProblem)
    #if !isnothing(model.polyhedral_uncertainty_set)
        #@constraint(model, )
    #end
    #if !isnothing(model.ellipsodial_uncertainty_set)
        #@constraint(model, )
    #end
    return nothing
end

get_xbar(t::DefaultExt, alg::AbstractSIPAlgo, s::LowerLevel1, sr::SIPSubResult) = sr.lbd.sol
get_xbar(t::DefaultExt, alg::AbstractSIPAlgo, s::LowerLevel2, sr::SIPSubResult) = sr.ubd.sol
get_xbar(t::DefaultExt, alg::AbstractSIPAlgo, s::LowerLevel3, sr::SIPSubResult) = sr.lbd.res

function llp_check(islocal::Bool, t::MOI.TerminationStatusCode, r::MOI.ResultStatusCode)
    flag = true
    if islocal && ((t != MOI.LOCALLY_SOLVED) && (t != MOI.ALMOST_LOCALLY_SOLVED))
        error("Lower problem did not solve to local optimality.")
    elseif t !== MOI.OPTIMAL
        error("Error in lower level problem. Termination status = $t, primal status = $r.")
    end
    return flag
end

"""
    sip_llp!

Solves the lower level problem for the ith-SIP used with extension `t::EAGO.ExtensionType`
in algorithm `a::AbstractSIPAlgo` in subproblem `s::AbstractSubproblemType` via
the command  `sip_llp!(t::ExtensionType, a::AbstractSIPAlgo, s::AbstractSubproblemType, ..., i, tol)`.
"""
function sip_llp!(t::DefaultExt, alg::A, s::S, result::SIPResult,
                  sr::SIPSubResult, prob::SIPProblem, cb::SIPCallback,
                  i::Int64, tol::Float64 = -Inf) where {A <: AbstractSIPAlgo, S <: AbstractSubproblemType}

    # build the model
    m, p = build_model(t, alg, s, prob)
    set_tolerance!(t, alg, s, m, sr, i)

    # define the objective
    xbar = get_xbar(t, alg, s, sr)
    g(p...) = cb.gSIP[i](xbar, p)
    register(m, :g, prob.np, g, autodiff=true)
    if isone(prob.np)
        nl_obj = :(g($(p[1])))
    else
        nl_obj = Expr(:call)
        push!(nl_obj.args, :g)
        for i in 1:prob.np
            push!(nl_obj.args, p[i])
        end
    end
    set_nonlinear_objective(m, MOI.MAX_SENSE, nl_obj)

    # add uncertainty constraints
    add_uncertainty_constraint!(m, prob)

    # optimize model and check status
    JuMP.optimize!(m)
    tstatus = JuMP.termination_status(m)
    rstatus = JuMP.primal_status(m)
    feas = llp_check(prob.local_solver, tstatus, rstatus)

    # fill buffer with subproblem result info
    psol = JuMP.value.(p)
    load!(s, sr, feas, JuMP.objective_value(m), JuMP.objective_bound(m), psol)
    result.solution_time += MOI.get(m, MOI.SolveTimeSec())

    return nothing
end
function sip_llp!(t::ExtensionType, alg::A, s::S, result::SIPResult,
                     sr::SIPSubResult, prob::SIPProblem, cb::SIPCallback,
                     i::Int64) where {A <: AbstractSIPAlgo, S <: AbstractSubproblemType}
    sip_llp!(DefaultSubproblem(), alg, s, result, sr, prob, cb, i)
end

"""
    sip_bnd!

Solves the bounding problem for the ith-SIP used with extension `t::EAGO.ExtensionType`
in algorithm `a::AbstractSIPAlgo` in subproblem `s::AbstractSubproblemType` via
the command `sip_bnd!(t::ExtensionType, a::AbstractSIPAlgo, s::AbstractSubproblemType, ..., i, tol)`.
"""
function sip_bnd!(t::ExtensionType, alg::A, s::S, sr::SIPSubResult, result::SIPResult,
                  prob::SIPProblem, cb::SIPCallback) where {A <: AbstractSIPAlgo,
                                                            S <: AbstractSubproblemType}
    sip_bnd!(DefaultExt(), alg, s, sr, result, prob, cb)
end
function sip_bnd!(t::DefaultExt, alg::A, s::S, sr::SIPSubResult, result::SIPResult,
                  prob::SIPProblem, cb::SIPCallback) where {A <: AbstractSIPAlgo,
                                                            S <: AbstractSubproblemType}

    # create JuMP model
    m, x = build_model(t, alg, s, prob)

    for i = 1:prob.nSIP
        ϵ_g = get_eps(s, sr, i)
        disc_set = get_disc_set(t, alg, s, sr, i)
        for j = 1:length(disc_set)
            gi = Symbol("g$i$j")
            g(x...) = cb.gSIP[i](x, disc_set[j])
            register(m, gi, prob.nx, g, autodiff=true)
            gic = Expr(:call)
            push!(gic.args, gi)
            for i in 1:prob.nx
                push!(gic.args, x[i])
            end
            JuMP.add_nonlinear_constraint(m, :($gic + $ϵ_g <= 0))
        end
    end

    # define the objective
    obj(x...) = cb.f(x)
    register(m, :obj, prob.nx, obj, autodiff=true)
    if isone(prob.nx)
        nl_obj = :(obj($(x[1])))
    else
        nl_obj = Expr(:call)
        push!(nl_obj.args, :obj)
        for i in 1:prob.nx
            push!(nl_obj.args, x[i])
        end
    end
    set_nonlinear_objective(m, MOI.MIN_SENSE, nl_obj)

    # optimize model and check status
    JuMP.optimize!(m)
    t_status = JuMP.termination_status(m)
    r_status = JuMP.primal_status(m)
    feas = bnd_check(prob.local_solver, t_status, r_status)

    # fill buffer with subproblem result info
    load!(s, sr, feas, JuMP.objective_value(m), JuMP.objective_bound(m), JuMP.value.(x))
    result.solution_time += MOI.get(m, MOI.SolveTimeSec())

    return nothing
end

"""
    sip_res!

Solves the restriction problem for extension `t::EAGO.ExtensionType`
in algorithm `a::AbstractSIPAlgo` in subproblem `s::AbstractSubproblemType` via
the command `sip_res!(t::ExtensionType, a::AbstractSIPAlgo, ...)`.
"""
function sip_res!(t::ExtensionType, alg::A, sr::SIPSubResult, result::SIPResult,
                  prob::SIPProblem, cb::SIPCallback) where {A <: AbstractSIPAlgo}
    sip_res!(DefaultExt(), alg, sr, result, prob, cb)
end
function sip_res!(t::DefaultExt, alg::A, sr::SIPSubResult, result::SIPResult,
                  prob::SIPProblem, cb::SIPCallback) where {A <: AbstractSIPAlgo}

    # create JuMP model & variables
    s = ResProblem()
    m, x = build_model(t, alg, s, prob)
    @variable(m, η)

    # add discretized semi-infinite constraint
    for i = 1:prob.nSIP
        disc_set = get_disc_set(t, alg, s, sr, i)
        for j = 1:length(disc_set)
            gi = Symbol("g$i$j")
            g(x...) = cb.gSIP[i](x, disc_set[j])
            register(m, gi, prob.nx, g, autodiff=true)
            gic = Expr(:call)
            push!(gic.args, gi)
            for i in 1:prob.nx
                push!(gic.args, x[i])
            end
            gip = Expr(:call)
            push!(gip.args, :+)
            push!(gip.args, gic)
            push!(gip.args, η)
            JuMP.add_nonlinear_constraint(m, :($gip <= 0))
        end
    end

    # add epigraph reformulated objective
    obj(x...) = cb.f(x)
    register(m, :f, prob.nx, obj, autodiff=true)
    if isfinite(sr.fRes)
        if isone(prob.nx)
            nl_obj = :(f($(x[1])))
        else
            nl_obj = Expr(:call)
            push!(nl_obj.args, :f)
            for i in 1:prob.nx
                push!(nl_obj.args, x[i])
            end
        end
        JuMP.add_nonlinear_constraint(m, :($nl_obj + $(sr.fRes) <= 0))
    end

    # define the objective
    @objective(m, Min, -η)

    # optimize model and check status
    JuMP.optimize!(m)
    t_status = JuMP.termination_status(m)
    r_status = JuMP.primal_status(m)
    feas = llp_check(prob.local_solver, t_status, r_status)

    # fill buffer with subproblem result info
    load!(s, sr, feas, JuMP.objective_value(m), JuMP.objective_bound(m), JuMP.value.(x))
    result.solution_time += MOI.get(m, MOI.SolveTimeSec())

    return nothing
end

# get optimizer for use in subproblems
function get_sip_optimizer(t::DefaultExt, alg::A, s::S) where {A<:AbstractSIPAlgo, S<:AbstractSubproblemType}
    return EAGO.Optimizer
end

"""
    get_sip_optimizer

Specifices the optimizer to be used in extension `t::EAGO.ExtensionType` with
algorithm `alg::AbstractSIPAlgo` in subproblem `s::AbstractSubproblemType` via the
command `get_sip_optimizer(t::ExtensionType, alg::AbstractSIPAlgo, s::AbstractSubproblemType)`.
"""
function get_sip_optimizer(t::ExtensionType, alg::A, s::AbstractSubproblemType) where {A <: AbstractSIPAlgo}
    return get_sip_optimizer(DefaultExt(), alg, s)
end

# Printing
function print_int!(verb::Int, prob::SIPProblem, result::SIPResult, r::Float64)

    k = result.iteration_number
    lbd = result.lower_bound
    ubd = result.upper_bound

    if (prob.verbosity == 1 || prob.verbosity == 2)

        # prints header line every hdr_intv times
        if (mod(k, prob.header_interval) == 0 || k == 1)
            println("| Iteration | Lower Bound | Upper Bound |   r   |  Gap  |  Ratio  |")
        end

        # prints iteration summary every prnt_intv times
        if mod(k, prob.print_interval) == 0

            print_str = "| "

            max_len = 15
            temp_str = string(k)
            len_str = length(temp_str)
            print_str *= (" "^(max_len - len_str))*temp_str*" | "

            max_len = 15
            lower_adj = lbd
            temp_str = string(round(lower_adj, digits = 6))
            len_str = length(temp_str)
            print_str *= (" "^(max_len - len_str))*temp_str*" | "

            max_len = 15
            upper_adj = ubd
            temp_str = string(round(upper_adj, digits = 6))
            len_str = length(temp_str)
            print_str *= (" "^(max_len - len_str))*temp_str*" | "

            max_len =15
            temp_str = string(round(r, digits = 5))
            len_str = length(temp_str)
            print_str *= (" "^(max_len - len_str))*temp_str*" | "

            max_len = 15
            temp_str = string(round((ubd-lbd), digits = 5))
            len_str = length(temp_str)
            print_str *= (" "^(max_len - len_str))*temp_str*" | "

            max_len = 18
            temp_str = string(round((ubd-lbd)/abs(ubd), digits = 6))
            len_str = length(temp_str)
            print_str *= (" "^(max_len - len_str))*temp_str*" |"

            println(print_str)
        end
    end
    return nothing
end

function summary_inner!(verb::Int64, val::Float64, x::Vector{Float64},
                        feas::Bool, desc::String)
    (verb == 2 || verb == 1) && println("solved $desc: ", val, " ", x, " ", feas)
    return nothing
end

for (typ, fd) in SUBPROB_SYM
    @eval function print_summary!(s::$typ, v::Int64, r::SIPSubResult, i::Int = 0)
        if i == 0
            summary_inner!(v, r.$fd.obj_val, r.$fd.sol, r.$fd.feas, "$s ")
        else
            summary_inner!(v, r.$fd.obj_val, r.$fd.sol, r.$fd.feas, "$s($i) ")
        end
    end
end

# Check convergence
function check_convergence(result::SIPResult, atol::Float64, verb::Int64)
    if abs(result.upper_bound -  result.lower_bound) < atol
        (verb == 2 || verb == 1) && println("Algorithm Converged")
        return true
    end
    return false
end

get_bnds(s::Union{LowerLevel1,LowerLevel2,LowerLevel3}, p::SIPProblem) = p.p_l, p.p_u, p.np
get_bnds(s::Union{LowerProblem,UpperProblem, ResProblem}, p::SIPProblem) = p.x_l, p.x_u, p.nx

function bnd_check(is_local::Bool, t::MOI.TerminationStatusCode, r::MOI.ResultStatusCode)
    if t === MOI.OPTIMAL
        return true
    elseif !is_local
        error("Lower problem did not solve to global optimality.
               Termination status = $t. Primal status = $r")
    elseif is_local && !((t == MOI.LOCALLY_SOLVED) || (t == MOI.ALMOST_LOCALLY_SOLVED))
        error("Lower problem did not solve to local optimality.")
    end
    return false
end

max_stars_repo_name: MatthewStuber/EAGO

max_issues_repo_name: MatthewStuber/EAGO

max_forks_repo_name: MatthewStuber/EAGO


text:
module Web
using HTTP
using JuliaDB
include("./JuliaDBQuery.jl")
using .JuliaDBQuery
using CSVFiles

function write_row(http, row)
    write(http,join(row,","))
    write(http,"\n")
end

function write_chunk(http, chunk)
    local rows = chunk |> collect
    foreach(row -> write_row(http, row), rows)
end

const BASE_PATH = r"/api/v1/.*"
const API_PATH = r"/api/v1/([a-zA-Z0-9_]+)"

"""
This handles the http get requests
The url is expected to be of the form /api/{table}
and it will return the contents of the table, after the roles have been applied
in the form the accept header has asked for
"""
function handle_get(http::HTTP.Stream)
    local message = http.message
    local path = match(API_PATH,http.message.target)
    if isnothing(path)
        local path = match(BASE_PATH,http.message.target)
        startwrite(http)
        if isnothing(path)
            write(http,"our api urls are of the form /api/v1/{table}\n")
        else
            write(http,"nice try, but sadly, we check our API paths for that kind of shenanigans\n")
        end
        closewrite(http)
        return
    end
    t = JuliaDBQuery.load_table(path[1])
    startwrite(http)
    write_row(http, colnames(t)) # write header
    foreach(chunk -> write_chunk(http, chunk),t.chunks)
    closewrite(http)
end

function handle_post(http::HTTP.Stream)
    local message = http.message
    local query = join(readlines(http))
    println("QUERY:" , query)
    local path = match(API_PATH,http.message.target)
    if isnothing(path)
        local path = match(BASE_PATH,http.message.target)
        startwrite(http)
        if isnothing(path)
            write(http,"our api urls are of the form /api/v1/{table}\n")
        else
            write(http,"nice try, but sadly, we check our API paths for that kind of shenanigans\n")
        end
        closewrite(http)
        return
    end
    local table = JuliaDBQuery.load_table(path[1])
    local result = JuliaDBQuery.run_basic_query(table,query)
    startwrite(http)
    write_row(http, colnames(result)) # write header
    foreach(chunk -> write_chunk(http, chunk),result.chunks)
    closewrite(http)
end

"""
since this is AFTER the proxy, we can just ask for roles in the header
the role list is space seperated
"user admin" means someone is both a user and an admin
"""
function auth_handler(http::HTTP.Stream)
    @show http.message
    @show http.message.target
    @show http.message.method
    local roles = split(HTTP.header(http, "roles", "anonymous"))
    local token = HTTP.header(http, "Authorization", "")
    if token != "qWohIIuX5Oc0DxPHUAyX"
        HTTP.setstatus(http, 404)
        startwrite(http)
        write(http, "Bad auth, contact data ventures for a token\n")
        closewrite(http)
        return
    end
    println("roles: $roles")
    task_local_storage("roles", roles) do
        if http.message.method=="POST"
            handle_post(http)
        else
            handle_get(http)
        end
    end
end

function start()
    println("starting")
    HTTP.serve(auth_handler,"0.0.0.0"; stream=true)
end

end


text:
#=
    --------------------------------------------------------------------------
    Implements functionality for learning a mean diffusion's trajectory during
    mcmc sampling. The main object is `Adaptation`.
    --------------------------------------------------------------------------
=#
import Base.resize!

"""
    Adaptation{TV,T}

Stores the history of imputed paths from which the mean trajectory can be
computed, as well as the scheme according to which the adaptation is supposed
to be performed. `TV` indicates whether any adaptation should be done at all.
`TV` set to `Val{False}` acts as an indicator that no adaptation is to be done.
"""
struct Adaptation{TV,T}
    X::Vector{Vector{Vector{T}}} # history of paths
    λs::Vector{Float64}          # ladder of weights that balance initial
                                 # auxiliary law and the adaptive law based on
                                 # the mean trajectory
    sizes::Vector{Int64}         # ladder indicating number of paths that are to
                                 # be used for computing mean trajectory
    skip::Int64                  # save 1 in every ... many sampled paths
    N::Vector{Int64}             # counter #1-current position on the ladder
                                 # #2-current index of the last saved path

    """
        Adaptation(::T, λs, sizes_of_path_coll, skip=1)

    Initialise adaptation. `λs` and `sizes_of_path_coll` are ladders that
    are traversed during sampling.
    ...
    # Arguments
    - `::T`: Data type of a diffusion
    - `λs`: ladder of weights that balance between initial choice of auxiliary
            law and the adaptive law based on the mean trajectory
    - `sizes_of_path_coll`: ladder giving the number of paths that are to be
                            used for computing the mean trajectory
    - `skip`: save 1 in every ... many sampled paths
    ...
    """
    function Adaptation(::T, λs, sizes_of_path_coll, skip=1) where T
        TV = Val{true}
        M = maximum(sizes_of_path_coll)
        X = [[zeros(T,0)] for i in 1:M]
        N = [1,1]
        new{TV,T}(X, λs, sizes_of_path_coll, skip, N)
    end

    """
        Adaptation{TV,T}()

    Empty constructor.
    """
    Adaptation{TV,T}() where {TV,T} = new{TV,T}()
end

"""
    NoAdaptation()

Helper function for constructing a flag saying that no adaptation is to be done
"""
NoAdaptation() = Adaptation{Val{false},Nothing}()

"""
    check_if_adapt(::Adaptation{Val{T}})

Check if any adaptation needs to be done
"""
check_if_adapt(::Adaptation{Val{T}}) where T = T

"""
    still_adapting(adpt::Adaptation{Val{true}})

If the adaptation has not been completed then do nothing, else return a flag
that no further adaptation is to be done
"""
function still_adapting(adpt::Adaptation{Val{true}})
    adpt.N[1] > length(adpt.sizes) ? NoAdaptation() : adpt
end

"""
    still_adapting(adpt::Adaptation{Val{false}})

There is nothing to be done for a flag indicating no adaptation
"""
still_adapting(adpt::Adaptation{Val{false}}) = adpt

"""
    resize!(adpt::Adaptation{TV,T}, m, ns::Vector{Int64})

Resize internal containers `X` with paths so that each storage unit consists of
`m` subunits and each of these subunits is a length `ns[i]` vector of type `T`
"""
function resize!(adpt::Adaptation{TV,T}, m, ns::Vector{Int64}) where {TV,T}
    K = length(adpt.X)
    for i in 1:K
        adpt.X[i] = [[zero(T) for _ in 1:ns[i]] for i in 1:m]
    end
end

"""
    add_path!(adpt::Adaptation{Val{true},T}, X::Vector{SamplePath{T}}, i)

Save path `X` into the history stored by `adpt` object. Do so only if the index
`i` of the current update step is not supposed to be skipped.
"""
function add_path!(adpt::Adaptation{Val{true},T}, X::Vector{SamplePath{T}}, i) where T
    if i % adpt.skip == 0
        m = length(X)
        for j in 1:m
            adpt.X[adpt.N[2]][j] .= X[j].yy
        end
    end
end

"""
    add_path!(adpt::Adaptation{Val{false}}, ::Any, ::Any)

Nothing to be done when no adaptation is to be performed
"""
add_path!(adpt::Adaptation{Val{false}}, ::Any, ::Any) = false

"""
    mean_trajectory(adpt::Adaptation{Val{true}})

Compute the mean trajectory from the history of accepted paths stored in
`adpt.X`
"""
function mean_trajectory(adpt::Adaptation{Val{true}})
    X = adpt.X
    num_paths = adpt.sizes[adpt.N[1]]
    num_segments = length(X[1])
    for i in 2:num_paths
        for j in 1:num_segments
            num_pts = length(X[i][j])
            for k in 1:num_pts
                X[1][j][k] += X[i][j][k]
            end
        end
    end
    for j in 1:num_segments
        num_pts = length(X[1][j])
        for k in 1:num_pts
            X[1][j][k] /= num_paths
        end
    end
    X[1]
end

"""
    print_adaptation_info(adpt::Adaptation{Val{false}}, ::Any, ::Any, ::Any)

Nothing to print for no adaptation
"""
print_adaptation_info(adpt::Adaptation{Val{false}}, ::Any, ::Any, ::Any) = nothing

"""
    print_adaptation_info(adpt::Adaptation{Val{true}}, acc_imp_counter,
                          acc_updt_counter, i)

Print some information regarding acceptance rate during adaptation to the
console
"""
function print_adaptation_info(adpt::Adaptation{Val{true}}, acc_imp_counter,
                               accUpdtCounter, i)
    if i % adpt.skip == 0 && adpt.N[2] == adpt.sizes[adpt.N[1]]
        print("--------------------------------------------------------\n")
        print(" Adapting...\n")
        print(" Using ", adpt.N[2], " many paths, thinned by ", adpt.skip, "\n")
        print(" Previous imputation acceptance rate: ", acc_imp_counter/i, "\n")
        print(" Previous param update acceptance rate: ", acc_updt_counter./i, "\n")
        print("--------------------------------------------------------\n")
    end
end


text:
#!/mnt/e/Julia-1.6.2
#########################################################################
# Disorder_Reader, Version 1.0.0
# 2021, ShanghaiTech University, Wang Lin
# Email: wanglin3@shanghaitech.edu.cn
#########################################################################

using Base: Float64, Float32
using BioStructures

function get_next_line_frag(chain, frag_num, frag_start::Int64, chain_res_end, frag_matrix, pdb_name, model_num, chain_id)
    truncate_domain = collectatoms(chain, res -> resnumber(res) >= frag_start, calphaselector)
    frag_checkpoints::Int64 = 0
    update_matrix = false
    first_res = true
    residue_num_last::Int64 = 0
    frag_end::Int64 = 0
    for residue in truncate_domain
        b_factor_p::Float32 = tempfactor(residue)
        residue_num_p = resnumber(residue)
        chain_res_end == residue_num_p && begin
            frag_end = residue_num_p
            frag_checkpoints > 0 && ( update_matrix = true )
            break
        end
        first_res == true || begin
            if residue_num_last != residue_num_p - 1
                frag_end = residue_num_last
                (frag_checkpoints > 0) && ( update_matrix = true )
                println("Note: Gap $residue_num_last - $residue_num_p in $pdb_name-$model_num-$chain_id")
                break
            end
        end
        residue_num_last = residue_num_p
        if b_factor_p >= b_factor_min && b_factor_p <= b_factor_max
            ( frag_checkpoints == 0 ) && ( frag_start = residue_num_p )
            ( frag_checkpoints <= order_len_max ) && ( frag_checkpoints = frag_checkpoints + 1 )
        else
            if frag_checkpoints > 0
                frag_checkpoints = frag_checkpoints - 1
                ( frag_checkpoints == 0 ) && begin
                    frag_end = residue_num_p
                    update_matrix = true
                    break
                end
            end
        end
        first_res = false
    end
    ( update_matrix == true ) && begin
        frag_len = frag_end - frag_start + 1
        frag_len >= domain_min_len && ( push!(frag_matrix,"$frag_num,$frag_start,$frag_end") )
    end
    return frag_matrix, frag_end
end

function get_terminal_resnum(chain::Chain)
    C_end::Int64 = 1
    N_start::Int64 = 9999999999
    for residue in chain
        residue_num_p::Int64 = resnumber(residue)
        ( C_end <= residue_num_p ) && ( C_end = residue_num_p )
        ( N_start >= residue_num_p ) && ( N_start = residue_num_p )
    end
    return N_start::Int64,C_end::Int64
end

function get_fasta(chain, frag_matrix, pdb_name, model_num, chain_id,index_dir,pdb)
    fasta = open("$index_dir/$pdb_name-$model_num-$chain_id.fasta","w")
    for frag in frag_matrix
        frag_start = parse(Int64, split(frag, ",")[2])
        frag_end = parse(Int64, split(frag, ",")[3])
        frag_str = collectresidues(chain, res -> frag_start <= resnumber(res) <= frag_end )
        fasta_seq = LongAminoAcidSeq(frag_str)
        write(fasta, ">$pdb_name.$model_num.$chain_id:$frag_start-$frag_end\n")
        write(fasta, "$fasta_seq\n")
    end
end

function main(pdb::String, index_dir::String)
    Structure = read(pdb, PDB)
    pdb_name = chop(basename(pdb), tail = 4)
    for model in Structure
        model_num = modelnumber(model)
        for chain in model
            chain_id = chainid(chain)
            chain_res_start::Int64, chain_res_end::Int64= get_terminal_resnum(chain)
            frag_start::Int64 = chain_res_start
            frag_end::Int64 = chain_res_start
            frag_num = 0
            frag_matrix = []
            while frag_end < chain_res_end
                frag_num+=1
                frag_start = frag_end
                #将 frag_num, frag_start, frag_end 存入 frag_matrix
                frag_matrix, frag_end= get_next_line_frag(chain, frag_num, frag_start, chain_res_end, frag_matrix, pdb_name, model_num, chain_id)
            end
            get_fasta(chain, frag_matrix, pdb_name, model_num, chain_id,index_dir,pdb)
        end
    end
end

println("Input Dir for protein pdb files:")
const input_dir=readline(stdin)
println("output Dir for fasta files:")
const index_dir=readline(stdin)
global b_factor_max = 50
global b_factor_min = 0
global order_len_max = 2
global domain_min_len = 5

Threads.@threads for pdb_file in readdir(input_dir)
    pdb = "$(input_dir)/$(pdb_file)"
    endswith(pdb, ".pdb") && begin
        try
            main(pdb, index_dir)
        catch e1
            println("$e1 Error: $pdb")
        end
    end
end


text:
#!/usr/bin/env julia

using Luxor, Test, Random

function test_polytriangulate(fname)

    Random.seed!(42)

    width, height = 800, 800
    Drawing(width, height, fname)
    origin()
    background("white")

    rawpts = ngon(O, 300, 8, vertices = true)

    for i in 1:12
        push!(rawpts, rand(BoundingBox(rawpts) * 0.6))
    end

    @test length(rawpts) == 20

    ptriangles = polytriangulate(rawpts)

    @test length(ptriangles) == 30

    fontsize(16)
    # visual inspection...
    setline(3)
    for (n, p) in enumerate(ptriangles)
        sethue([Luxor.julia_purple, Luxor.julia_blue,  Luxor.julia_red,  Luxor.julia_green][mod1(n, end)])
        pgon = Point[p[1], p[2], p[3]]
        poly(pgon,  :fillpreserve, close = true)
        sethue("white")
        strokepath()
        text(string(n), polycentroid(pgon), halign=:middle)
    end

    sethue("white")
    circle.(rawpts, 5, :fill)

    # can't think of how to test the triangulation, ideas welcome
    @test isempty(polyintersect(ptriangles[1], ptriangles[2])) == false
    @test isempty(polyintersect(ptriangles[1], ptriangles[6])) == true
    @test isempty(polyintersect(ptriangles[1], ptriangles[10])) == true
    @test isempty(polyintersect(ptriangles[1], ptriangles[20])) == true

    @test finish() == true
end

fname = "polytriangulate.png"
test_polytriangulate(fname)
println("...finished test: output in $(fname)")


text:
using Plots
using LaTeXStrings
using Base.Test

# provides c0, G, m_H, k_B, sigma_sb, a_rad, Msun, Rsun, Lsun
include("constants.jl")
include("physics.jl")
include("opacity.jl")
include("derivs.jl")
include("nr.jl")
include("shootf.jl")

logkappa_spl = opal_test_spl()

# solar values
mc = 0.001 * Msun
Pc = 2.4e17 # dyne / cm2
Tc = 1.6e7 # K
X = 0.73
Y = 0.24
Z = 1 - X - Y
mu = mu_from_composition(X, Y)

function test_nuclear(savefig=nothing)
    temps = logspace(6.5, 8)
    X = 0.7
    Y = 0.28
    Z = 0.02
    mu = mu_from_composition(X, Y)
    rho = rho_ideal(Pc, temps, mu)
    pp = epsilon_pp(rho, temps, X)
    cno = epsilon_cno(rho, temps, X, Y)

    print("pp/CNO at T6=18, rho=80: ", epsilon_pp(80, 18e6, X) / epsilon_cno(80, 18e6, X, Y))

    plot(temps, pp, label="pp",
         xlabel="Temperature [K]", xscale=:log,
         ylabel=L"Luminosity density [g cm$^{-3}$]", yscale=:log)
    plot!(temps, cno, label="CNO")
    plot!(temps, pp + cno, label="pp + CNO")
    if ~is(savefig, nothing)
        png(savefig)
    end # if
end

function test_dTdm()
    mu = 0.6
    rhoc = rho_ideal(Pc, Tc, mu)
    kappac = 10 ^ logkappa_spl(log10(rhoc), log10(Tc))
    rc, lc, _, __ = load1(Pc, Tc, X, Y)
    mc = 0.001 * Msun
    dTdm(mc, rc, lc, Pc, Tc, kappac)
end

function test_profile(profile; mass=5)
    X = 0.7
    Y = 0.28
    m = 0.1 * Msun
    M = mass * Msun
    mf = 0.8 * M
    Rs, Ls, Pc, Tc = init_guess(M, X, Y)
    Rs, Ls, Ps, Ts = load2(M, Rs, Ls, X, Y)
    println("surface: ", [Rs, Ls, Ps, Ts])
    Rc, Lc, Pc, Tc = load1(m, Pc, Tc, X, Y)
    println("center: ", [Rc, Lc, Pc, Tc])

    rhos = rho_ideal(Ps, Ts, mu)
    rhoc = rho_ideal(Pc, Tc, mu)
    kappac = opacity(logkappa_spl, rhoc, Tc)
    kappas = opacity(logkappa_spl, rhos, Ts)
    println("kappa at center: ", kappac)
    println("log R at center: ", log10(rhoc / (Tc * 1e-6) ^ 3))
    println("kappa at surface: ", kappas)
    println("log R at surface: ", log10(rhos / (Ts * 1e-6) ^ 3))

    x1, y1, x2, y2 = profiles(m, M, Rs, Ls, Pc, Tc, X, Y, mf)

    idx_dict = Dict("r" => 1, "l" => 2, "P" => 3, "T" => 4)
    idx = idx_dict[profile]
    if idx == 1
        norm = Rsun
    elseif idx == 2
        norm = Lsun
    else
        norm = 1
    end
    if all(y1[1:end, idx] .== y1[1, idx])
        println("outward flatlined")
    end
    if all(y2[1:end, idx] .== y2[1, idx])
        println("inward flatlined")
    end

    println(profile, "_out(mf) = ", y1[end, idx])
    println(profile, "_in(mf) = ", y2[end, idx])
    plot(x1 / Msun, y1[1:end, idx] / norm, label=profile * "out", yscale=:log10)
    plot!(x2 / Msun, y2[1:end, idx] / norm, label=profile * "in", yscale=:log10)
    # return x1, x2, y1[1:end, idx], y2[1:end, idx]
end

function test_load()
    # load1
    center = load1(mc, Pc, Tc, X, Y)
    println("center ", center)
    # load2
    surface = load2(Msun, Rsun, Lsun, X, Y)
    println("surface ", surface)
    Rs, Ls, Ps, Ts = surface
    rhos = rho_ideal(Ps, Ts, mu)
    println(rhos)
    kappa = 10 ^ logkappa_spl(log10(rhos), log10(Ts))
    # test that the opacity gives a pressure consistant to within 0.001 percent
    @test_approx_eq_eps Ps P_surface(Msun, Rsun, kappa) (Ps * 1e-5)
end

function test_deriv()
    rc, lc, _, _ = load1(mc, Pc, Tc, X, Y)
    println("deriv at center: ", deriv(mc, rc, lc, Pc, Tc, X, Y))
    rs, ls, Ps, Ts = load2(Msun, Rsun, Lsun, X, Y)
    println("deriv at surface: ", deriv(Msun, rs, ls, Ps, Ts, X, Y))
end

function test_opacity(M=5 * Msun)
    X = 0.7
    Y = 0.28
    mu = mu_from_composition(X, Y)

    m = 0.01 * Msun
    # M = 5 * Msun
    mf = 0.8 * M
    Rs, Ls, Pc, Tc = init_guess(M, X, Y)
    Rs, Ls, Ps, Ts = load2(M, Rs, Ls, X, Y)
    println("surface: ", [Rs, Ls, Ps, Ts])
    Rc, Lc, Pc, Tc = load1(m, Pc, Tc, X, Y)
    println("center: ", [Rc, Lc, Pc, Tc])
    rhos = rho_ideal(Ps, Ts, mu)
    rhoc = rho_ideal(Pc, Tc, mu)
    temps = logspace(Ts, Tc)
    dens = logspace(rhos, rhoc)
    kappac = opacity(logkappa_spl, rhoc, Tc)
    kappas = opacity(logkappa_spl, rhos, Ts)
    println("kappa at center: ", kappac)
    println("log R at center: ", log10(rhoc / (Tc * 1e-6) ^ 3))

    println("kappa at surface: ", kappas)
    println("log R at surface: ", log10(rhos / (Ts * 1e-6) ^ 3))
end


text:
function estimate_return(::IS, π, θ, τ::Trajectory)
    G = sum(τ.rewards)
    blogp = sum(τ.blogps)
    logp = sum(@. logpdf((π,), τ.states, τ.actions))
    ρ = exp(logp - blogp)
    return exp(logp - blogp) * G
end

function estimate_return(::WIS, (π,), τ::Trajectory)
    G = sum(τ.rewards)
    blogp = sum(τ.blogps)
    logp = sum(@. logpdf((π,), τ.states, τ.actions))
    ρ = exp(logp - blogp)
    return ρ * G, ρ
end

function iterativePDIS(π, τ::Trajectory{T}, t, lnρ) where {T}
    lnρ += logpdf(π, τ.states[t], τ.actions[t]) - τ.blogps[t]
    ρ = exp(lnρ)
    return ρ * τ.rewards[t], ρ, lnρ
end

function estimate_return(::PDIS, π, τ::Trajectory{T}) where {T}
    N = length(τ)
    lnρ = T(0.0)
    G = T(0.0)
    for t in 1:N
        ρR, _, lnρ = iterativePDIS(π, τ, t, lnρ)
        G += ρR
    end
    return G
end

function estimate_returns!(G, ism::CWPDIS, π, H)
    N = length(H)
    lnρs = zeros(N)
    ρtot = 0.0
    L = maximum(length.(H))
    for i in 1:N
        for t in 1:L
            ptot = 0.0
            if length(H[i]) ≤ t
                ρR, ρ, lnρ = iterativePDIS(π, H[i], t, lnρs[i])
                G[i] += ρR
                lnρs[i] += lnρ
                ptot += ρ
            else
                ptot += exp(lnρs[i])
            end
            @. G /= ptot
        end
    end
end

function estimate_returns(ism::CWPDIS, π, H)
    N = length(H)
    G = zeros(N)
    estimate_returns!(G, ism, π, H)
    return G
end


text:
module LibTreeSitterHighlight

using tree_sitter_highlight_jll: libtree_sitter_highlight
using CEnum: @cenum

@cenum TSHighlightError begin
    TSHighlightOk
    TSHighlightUnknownScope
    TSHighlightTimeout
    TSHighlightInvalidLanguage
    TSHighlightInvalidUtf8
    TSHighlightInvalidRegex
    TSHighlightInvalidQuery
end

struct TSHighlighter end
struct TSHighlightBuffer end
struct TSLanguage end

# TSHighlighter *ts_highlighter_new(
#   const char **highlight_names,
#   const char **attribute_strings,
#   uint32_t highlight_count
# );
function ts_highlighter_new(highlight_names, attribute_strings, highlight_count)
    @ccall libtree_sitter_highlight.ts_highlighter_new(
        highlight_names::Ptr{Ptr{Cchar}},
        attribute_strings::Ptr{Ptr{Cchar}},
        highlight_count::Cuint,
    )::Ptr{TSHighlighter}
end

# // Delete a syntax highlighter.
# void ts_highlighter_delete(TSHighlighter *);
function ts_highlighter_delete(highlighter)
    @ccall libtree_sitter_highlight.ts_highlighter_delete(
        highlighter::Ptr{TSHighlighter},
    )::Cvoid
end

# // Add a `TSLanguage` to a highlighter. The language is associated with a
# // scope name, which can be used later to select a language for syntax
# // highlighting. Along with the language, you must provide a JSON string
# // containing the compiled PropertySheet to use for syntax highlighting
# // with that language. You can also optionally provide an 'injection regex',
# // which is used to detect when this language has been embedded in a document
# // written in a different language.
# TSHighlightError ts_highlighter_add_language(
#   TSHighlighter *self,
#   const char *scope_name,
#   const char *injection_regex,
#   const TSLanguage *language,
#   const char *highlight_query,
#   const char *injection_query,
#   const char *locals_query,
#   uint32_t highlight_query_len,
#   uint32_t injection_query_len,
#   uint32_t locals_query_len
# );
function ts_highlighter_add_language(
    self,
    scope_name,
    injection_regex,
    language,
    highlight_query,
    injection_query,
    locals_query,
    highlight_query_len,
    injection_query_len,
    locals_query_len,
)
    @ccall libtree_sitter_highlight.ts_highlighter_add_language(
        self::Ptr{TSHighlighter},
        scope_name::Cstring,
        injection_regex::Cstring,
        language::Ptr{TSLanguage},
        highlight_query::Cstring,
        injection_query::Cstring,
        locals_query::Cstring,
        highlight_query_len::Cuint,
        injection_query_len::Cuint,
        locals_query_len::Cuint,
    )::TSHighlightError
end

# // Compute syntax highlighting for a given document. You must first
# // create a `TSHighlightBuffer` to hold the output.
# TSHighlightError ts_highlighter_highlight(
#   const TSHighlighter *self,
#   const char *scope_name,
#   const char *source_code,
#   uint32_t source_code_len,
#   TSHighlightBuffer *output,
#   const size_t *cancellation_flag
# );
function ts_highlighter_highlight(
    self,
    scope_name,
    source_code,
    source_code_len,
    output,
    cancellation_flag,
)
    @ccall libtree_sitter_highlight.ts_highlighter_highlight(
        self::Ptr{TSHighlighter},
        scope_name::Cstring,
        source_code::Cstring,
        source_code_len::Cuint,
        output::Ptr{TSHighlightBuffer},
        cancellation_flag::Ptr{Csize_t},
    )::TSHighlightError
end

# // TSHighlightBuffer: This struct stores the HTML output of syntax
# // highlighting. It can be reused for multiple highlighting calls.
# TSHighlightBuffer *ts_highlight_buffer_new();
function ts_highlight_buffer_new()
    @ccall libtree_sitter_highlight.ts_highlight_buffer_new()::Ptr{TSHighlightBuffer}
end

# // Delete a highlight buffer.
# void ts_highlight_buffer_delete(TSHighlightBuffer *);
function ts_highlight_buffer_delete(highlighter)
    @ccall libtree_sitter_highlight.ts_highlight_buffer_delete(
        highlighter::Ptr{TSHighlightBuffer},
    )::Cvoid
end

# // Access the HTML content of a highlight buffer.
# const uint8_t *ts_highlight_buffer_content(const TSHighlightBuffer *);
function ts_highlight_buffer_content(highlighter)
    @ccall libtree_sitter_highlight.ts_highlight_buffer_content(
        highlighter::Ptr{TSHighlightBuffer},
    )::Ptr{Cuint}
end

# const uint32_t *ts_highlight_buffer_line_offsets(const TSHighlightBuffer *);
function ts_highlight_buffer_line_offsets(highlight_buffer)
    @ccall libtree_sitter_highlight.ts_highlight_buffer_line_offsets(
        highlight_buffer::Ptr{TSHighlightBuffer},
    )::Ptr{Cuint}
end

# uint32_t ts_highlight_buffer_len(const TSHighlightBuffer *);
function ts_highlight_buffer_len(highlight_buffer)
    @ccall libtree_sitter_highlight.ts_highlight_buffer_len(
        highlight_buffer::Ptr{TSHighlightBuffer},
    )::Cuint
end

# uint32_t ts_highlight_buffer_line_count(const TSHighlightBuffer *);
function ts_highlight_buffer_line_count(highlight_buffer)
    @ccall libtree_sitter_highlight.ts_highlight_buffer_line_count(
        highlight_buffer::Ptr{TSHighlightBuffer},
    )::Cuint
end

end # module LibTreeSitter


text:
function vmul!(C, A, B)
    @debug "vmul!($(typeof(C)), $(typeof(A)), $(typeof(B)))"
    mul!(C, A, B)
end
function vmul!(C, A::Diagonal, B, range)
    @debug "vmul!($(typeof(C)), $(typeof(A)), $(typeof(B)), $(typeof(range)))"
    @views mul!(C, Diagonal(A.diag[range]), B)
end
function vmul!(C, A, B::Diagonal, range)
    @debug "vmul!($(typeof(C)), $(typeof(A)), $(typeof(B)), $(typeof(range)))"
    @views mul!(C, A, Diagonal(B.diag[range]))
end
function rvmul!(A, B)
    @debug "rvmul!($(typeof(A)), $(typeof(B)))"
    rmul!(A, B)
end
function lvmul!(A, B)
    @debug "lvmul!($(typeof(A)), $(typeof(B)))"
    lmul!(A, B)
end
function rvadd!(A, B)
    @debug "rvadd!($(typeof(A)), $(typeof(B)))"
    A .+= B
end
function vsub!(A, B, C)
    @debug "vsub!($(typeof(A)), $(typeof(B)), $(typeof(C)))"
    A .= B .- C
end
function vsub!(A, B, ::UniformScaling)
    @debug "vsub!($(typeof(A)), $(typeof(B)), ::UniformScaling)"
    T1 = one(eltype(A))
    T0 = zero(eltype(A))
    @inbounds for i in axes(A, 1), j in axes(A, 2)
        A[i, j] = B[i, j] - ifelse(i==j, T1, T0)
    end
    A
end

function rdivp!(A::Matrix{<: Complex}, T::Matrix{<: Complex}, O::Matrix{<: Complex}, pivot)
    @debug "rdivp!($(typeof(A)), $(typeof(T)), $(typeof(O)), $(typeof(pivot)))"
    # assume Diagonal is ±1!
    @inbounds begin
        N = size(A, 1)

        # Apply pivot
        for (j, p) in enumerate(pivot)
            for i in 1:N
                O[i, j] = A[i, p]
            end
        end

        # do the rdiv
        for j in 1:N
            for i in 1:N
                x = O[i, j]
                @simd for k in 1:j-1
                    x -= A[i, k] * T[k, j]
                end
                A[i, j] = x / T[j, j]
            end
        end
    end
    A
end


text:
using RenyiEntropyOfMagic
using OMEinsum
using SparseArrays
using LinearAlgebra
using Test

function naive_two_site(st)
    psi = reshape(kron(st, st, st, st), ntuple(_->2, 8))
    O = kron(Consts.O, Consts.O)
    P = ein"abcdefgh,acegbdfhijklmnop->ijklmnop"(psi,reshape(O, ntuple(_->2, 16)))
    return ein"abcdefgh,acegbdfh->"(psi, P)[]
end

function naive_three_site(st)
    psi = reshape(kron(st, st, st, st), ntuple(_->2, 12))
    O = kron(Consts.O, Consts.O, Consts.O)
    P = ein"abcdefghijkl,adgjbehkcfilmnopqrstuvwx->mnopqrstuvwx"(psi,reshape(O, ntuple(_->2, 24)))
    return ein"abcdefghijkl,adgjbehkcfil->"(psi, P)[]
end

@testset "random periodic boundary MPS" begin
    # random MPS
    A = rand(3, 3, 2)
    T = transfer_matrix(A)

    # tests
    # periodic boundary
    # 1 site
    st = ein"iij->j"(A)
    psi = kron(st, st, st, st)
    psi' * Consts.O * psi
    tr(T)
    @test psi' * Consts.O * psi ≈ tr(T)

    # 2 site
    st = vec(ein"ija,jib->ab"(A, A))
    @test naive_two_site(st) ≈ tr(T^2)

    A1 = rand(2, 3, 2)
    A2 = rand(3, 2, 2)
    T1 = transfer_matrix(A1)
    T2 = transfer_matrix(A2)
    st = vec(ein"ija,jib->ab"(A1, A2))
    @test naive_two_site(st) ≈ tr(T1 * T2)

    A1,A2,A3 = rand(2, 3, 2), rand(3, 3, 2), rand(3, 2, 2)
    st = vec(ein"ija,jkb,kic->abc"(A1, A2, A3))
    T1, T2, T3 = transfer_matrix(A1), transfer_matrix(A2), transfer_matrix(A3)
    @test naive_three_site(st) ≈ tr(T1 * T2 * T3)
end

@testset "random open boundary MPS" begin
    # 2 site
    lA, rA = rand(3, 2), rand(3, 2)
    lT, rT = boundary_vector(lA), boundary_vector(rA)

    st = vec(transpose(lA) * rA)
    expect = naive_two_site(st)
    @test expect ≈ dot(lT, rT)
    @test expect ≈ replicated_expectation([lA, rA])

    # 3 site
    lA, A, rA = rand(3, 2), rand(3, 3, 2), rand(3, 2)
    lT, T, rT = boundary_vector(lA), transfer_matrix(A), boundary_vector(rA)

    st = vec(ein"ia,ijb,jc->abc"(lA, A, rA))
    expect = naive_three_site(st)
    @test expect ≈ dot(lT, T * rT)
    @test expect ≈ replicated_expectation([lA, A, rA])
end

@testset "tucker decomposition" begin
    include("tucker.jl")
end


text:
mutable struct Params
    a; b
end

@testset "obstable" begin
import DataStructures

t = ObsTable(Params)
 # or as an equivalent alternative
t2 = ObsTable()
set_params_names!(t2, ["a", "b"])
@test params_names(t2) == params_names(t)
@test t2 == ObsTable(["a","b"]) == ObsTable((:a,"b"))

t2 = ObsTable()
set_params_names!(t2, :a, :b)
@test t2.par_names == DataStructures.OrderedSet(["a","b"])

t2[1,2] #create entry
@test t2[1,2] === t2[(1,2)] === t2[[1,2]]

range = zip(1.:10., 1.:10.)
for (x,y) in range
    par = Params(x,y)
    for i=1:1e3
        r1, r2 = [x,y] + randn(2)
        t[par]["sum"] &= r1 + r2
        t[par]["sum2"] &= r1^2 + r2^2
        t[par][:sum3] &= r1 + r2
    end
    @test t[par][:sum] == t[par][:sum3]
end

@test_throws AssertionError t2[:a]
@test_throws AssertionError t2[:b][:E] &= 1

for (x,y) in range
    par = Params(x,y)
    @test t[par] === t[(x,y)]
    ob = t[par]["sum"]
    @test isapprox(mean(ob), x+y, atol=5err(ob))
end
println(t)
open("test.dat","w") do f
    print(f, t)
end
t2 = ObsTable("test.dat")
println(t2)

t1 = ObsTable(Params)
t2 = ObsTable(Params)
p = Params(1,2)
t1[p]["mag"] &= 5
t2[p]["mag"] &= 6
t = merge(t1,t2)
@test mean(t[p]["mag"]) == 5.5

t2[p]["ene"] &= 2
t = merge(t1,t2)
@test length(t.par_names) == 2
@test length(t.data) == 1
@test length(t[p]) == 2
@test mean(t[p]["ene"]) == 2

p2 = Params(2,3)
t2[p2]["ene"] &= 3
t = merge(t1,t2)
@test length(t.par_names) == 2
@test length(t.data) == 2
@test length(t[p]) == 2
@test length(t[p2]) == 1
@test mean(t[p]["ene"]) == 2
@test mean(t[p2]["ene"]) == 3

@test mean(t[(p.a, p.b)][:ene]) == 2
@test mean(t[(p2.a, p2.b)][:ene]) == 3

println(t)

end #testset

max_stars_repo_name: JuliaTagBot/LittleScienceTools

max_issues_repo_name: JuliaTagBot/LittleScienceTools

max_forks_repo_name: JuliaTagBot/LittleScienceTools


text:
# # [Time-Domain Simulation of a Joined-Wing](@id dynamic-joined-wing)
#
# In this example we use the same joined-wing model as used in the [previous example](@ref
# static-joined-wing), but with the following time varying loads applied at the wingtip:
#  - A piecewise-linear load ``F_L`` in the x and y-directions defined as follows:
# ```math
# F_L(t) = \begin{cases}
#     t10^6 \text{ N} & 0 \leq t \leq 0.01 \\
#     (0.02-t)10^6 & 0.01 \leq t \leq 0.02 \\
#     0 & \text{otherwise}
# \end{cases}
# ```
# - A sinusoidal load ``F_S`` applied in the z-direction defined as follows:
# ```math
# F_S(t) = \begin{cases}
#     0 & t \lt 0 \\
#     5 \times 10^3 (1-\cos(\pi t /0.02)) \text{ N} & 0 \leq t \lt 0.02 \\
#     10^4 \text{ N} & 0.02 \leq t
# \end{cases}
# ```
#
# We will also use the same compliance and mass matrix for all beams, in order to simplify
# the problem definition.
#
# ![](../assets/static-joined-wing-drawing.png)
#
#-
#md # !!! tip
#md #     This example is also available as a Jupyter notebook:
#md #     [`dynamic-joined-wing.ipynb`](@__NBVIEWER_ROOT_URL__/examples/dynamic-joined-wing.ipynb).
#-

using GXBeam, LinearAlgebra

## Set endpoints of each beam
p1 = [0, 0, 0]
p2 = [-7.1726, -12, -3.21539]
p3 = [7.1726, -12,  3.21539]

Cab_1 = [
0.5         0.866025  0.0
0.836516    -0.482963  0.258819
0.224144     -0.12941   -0.965926
]

Cab_2 = [
0.5         0.866025  0.0
-0.836516    0.482963 0.258819
0.224144    -0.12941   0.965926
]

## beam 1
L_b1 = norm(p1-p2)
r_b1 = p2
nelem_b1 = 8
lengths_b1, xp_b1, xm_b1, Cab_b1 = discretize_beam(L_b1, r_b1, nelem_b1, frame=Cab_1)

## beam 2
L_b2 = norm(p3-p1)
r_b2 = p1
nelem_b2 = 8
lengths_b2, xp_b2, xm_b2, Cab_b2 = discretize_beam(L_b2, r_b2, nelem_b2, frame=Cab_2)

## combine elements and points into one array
nelem = nelem_b1 + nelem_b2
points = vcat(xp_b1, xp_b2[2:end])
start = 1:nelem
stop = 2:nelem + 1
lengths = vcat(lengths_b1, lengths_b2)
midpoints = vcat(xm_b1, xm_b2)
Cab = vcat(Cab_b1, Cab_b2)

## assign all beams the same compliance and mass matrix
compliance = fill(Diagonal([2.93944738387698e-10, 8.42991725049126e-10,
    3.38313996669689e-08, 4.69246721094557e-08, 6.79584100559513e-08,
    1.37068861370898e-09]), nelem)
mass = fill(Diagonal([4.86e-2, 4.86e-2, 4.86e-2, 1.0632465e-2, 2.10195e-4,
    1.042227e-2]), nelem)

## create assembly
assembly = Assembly(points, start, stop;
    compliance = compliance,
    mass = mass,
    frames = Cab,
    lengths = lengths,
    midpoints = midpoints)

F_L = (t) -> begin
    if 0.0 <= t < 0.01
        1e6*t
    elseif 0.01 <= t < 0.02
        -1e6*(t-0.02)
    else
        zero(t)
    end
end

F_S = (t) -> begin
    if t < 0.0
        zero(t)
    elseif 0.0 <= t < 0.02
        5e3*(1-cos(pi*t/0.02))
    else
        1e4
    end
end

## assign boundary conditions and point load
prescribed_conditions = (t) -> begin
    Dict(
        ## fixed endpoint on beam 1
        1 => PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_y=0, theta_z=0),
        ## force applied on point 4
        nelem_b1 + 1 => PrescribedConditions(Fx=F_L(t), Fy=F_L(t), Fz=F_S(t)),
        ## fixed endpoint on last beam
        nelem+1 => PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_y=0, theta_z=0),
    )
end

## time
t = range(0, 0.04, length=1001)

system, history, converged = time_domain_analysis(assembly, t;
    prescribed_conditions=prescribed_conditions)

#!jl nothing #hide

# We can visualize tip displacements and the resultant forces accessing the post-processed
# results for each time step contained in the variable `history`.  Note that the fore-root
# and rear-root resultant forces for this case are equal to the external forces/moments,
# but with opposite sign.

#md using Suppressor #hide

#md @suppress_err begin #hide

using Plots
pyplot()

point = vcat(fill(nelem_b1+1, 6), fill(1, 6))
field = [:u, :u, :u, :theta, :theta, :theta, :F, :F, :F, :M, :M, :M]
direction = [1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3]
ylabel = ["\$u_x\$ (\$m\$)", "\$u_y\$ (\$m\$)", "\$u_z\$ (\$m\$)",
    "Rodriguez Parameter \$\\theta_x\$", "Rodriguez Parameter \$\\theta_y\$",
    "Rodriguez Parameter \$\\theta_z\$", "\$F_x\$ at the forewing root (\$N\$)",
    "\$F_y\$ at the forewing root (\$N\$)", "\$F_z\$ at the forewing root (\$N\$)",
    "\$M_x\$ at the forewing root (\$Nm\$)", "\$M_y\$ at the forewing root (\$Nm\$)",
    "\$M_z\$ at the forewing root (\$N\$)"]
#nb ph = Vector{Any}(undef, 12)

for i = 1:12
#nb     ph[i] = plot(
#!nb     plot(
        xlim = (0, 0.04),
        xticks = 0:0.01:0.04,
        xlabel = "Time (s)",
        ylabel = ylabel[i],
        grid = false,
        overwrite_figure=false
        )
    y = [getproperty(state.points[point[i]], field[i])[direction[i]] for state in history]

    if field[i] == :theta
        ## convert to angle
        @. y = 4*atan(y/4)
    end

    if field[i] == :F || field[i] == :M
        y = -y
    end

    plot!(t, y, label="")
#!nb     plot!(show=true)
#md     savefig("../assets/dynamic-joined-wing-"*string(field[i])*string(direction[i])*".svg"); #hide
#md     closeall() #hide
end

#nb ph[1]
#nb #-
#nb ph[2]
#nb #-
#nb ph[3]
#nb #-
#nb ph[4]
#nb #-
#nb ph[5]
#nb #-
#nb ph[6]
#nb #-
#nb ph[7]
#nb #-
#nb ph[8]
#nb #-
#nb ph[9]
#nb #-
#nb ph[10]
#nb #-
#nb ph[11]
#nb #-
#nb ph[12]
#nb #-

#md end #hide
#md nothing #hide

#md # ![](../assets/dynamic-joined-wing-u1.svg)
#md # ![](../assets/dynamic-joined-wing-u2.svg)
#md # ![](../assets/dynamic-joined-wing-u3.svg)
#md # ![](../assets/dynamic-joined-wing-theta1.svg)
#md # ![](../assets/dynamic-joined-wing-theta2.svg)
#md # ![](../assets/dynamic-joined-wing-theta3.svg)
#md # ![](../assets/dynamic-joined-wing-F1.svg)
#md # ![](../assets/dynamic-joined-wing-F2.svg)
#md # ![](../assets/dynamic-joined-wing-F3.svg)
#md # ![](../assets/dynamic-joined-wing-M1.svg)
#md # ![](../assets/dynamic-joined-wing-M2.svg)
#md # ![](../assets/dynamic-joined-wing-M3.svg)

#md # ![](../assets/dynamic-joined-wing-u1.svg)
#md # ![](../assets/dynamic-joined-wing-u2.svg)
#md # ![](../assets/dynamic-joined-wing-u3.svg)
#md # ![](../assets/dynamic-joined-wing-theta1.svg)
#md # ![](../assets/dynamic-joined-wing-theta2.svg)
#md # ![](../assets/dynamic-joined-wing-theta3.svg)
#md # ![](../assets/dynamic-joined-wing-F1.svg)
#md # ![](../assets/dynamic-joined-wing-F2.svg)
#md # ![](../assets/dynamic-joined-wing-F3.svg)
#md # ![](../assets/dynamic-joined-wing-M1.svg)
#md # ![](../assets/dynamic-joined-wing-M2.svg)
#md # ![](../assets/dynamic-joined-wing-M3.svg)

#-

# These graphs are identical to those presented in "GEBT: A general-purpose nonlinear analysis tool for composite beams" by Wenbin Yu and Maxwell Blair.
#
# We can also visualize the time history of the system using ParaView.  In order to view the small deflections we'll scale all the deflections up by a couple orders of magnitude.  We'll also set the color gradient to match the magnitude of the deflections at each point.

airfoil  = [ #FX 60-100 airfoil
    0.0000000 0.0000000;
    0.0010700 0.0057400;
    0.0042800 0.0114400;
    0.0096100 0.0177500;
    0.0170400 0.0236800;
    0.0265300 0.0294800;
    0.0380600 0.0352300;
    0.0515600 0.0405600;
    0.0669900 0.0460900;
    0.0842700 0.0508600;
    0.1033200 0.0556900;
    0.1240800 0.0598900;
    0.1464500 0.0640400;
    0.1703300 0.0675400;
    0.1956200 0.0708100;
    0.2222100 0.0733900;
    0.2500000 0.0756500;
    0.2788600 0.0772000;
    0.3086600 0.0783800;
    0.3392800 0.0788800;
    0.3705900 0.0789800;
    0.4024500 0.0784500;
    0.4347400 0.0775000;
    0.4673000 0.0759600;
    0.5000000 0.0740900;
    0.5327000 0.0717400;
    0.5652600 0.0691100;
    0.5975500 0.0660800;
    0.6294100 0.0627500;
    0.6607200 0.0590500;
    0.6913400 0.0551100;
    0.7211400 0.0508900;
    0.7500000 0.0465200;
    0.7777900 0.0420000;
    0.8043801 0.0374700;
    0.8296700 0.0329800;
    0.8535500 0.0286400;
    0.8759201 0.0244700;
    0.8966800 0.0205300;
    0.9157300 0.0168100;
    0.9330100 0.0134200;
    0.9484400 0.0103500;
    0.9619400 0.0076600;
    0.9734700 0.0053400;
    0.9829600 0.0034100;
    0.9903900 0.0019300;
    0.9957200 0.0008600;
    0.9989300 0.0002300;
    1.0000000 0.0000000;
    0.9989300 0.0001500;
    0.9957200 0.0007000;
    0.9903900 0.0015100;
    0.9829600 0.00251;
    0.9734700 0.00377;
    0.9619400 0.00515;
    0.9484400 0.00659;
    0.9330100 0.00802;
    0.9157300 0.00941;
    0.8966800 0.01072;
    0.8759201 0.01186;
    0.8535500 0.0128;
    0.8296700 0.01347;
    0.8043801 0.01381;
    0.7777900 0.01373;
    0.7500000 0.01329;
    0.7211400 0.01241;
    0.6913400 0.01118;
    0.6607200 0.00951;
    0.6294100 0.00748;
    0.5975500 0.00496;
    0.5652600 0.00217;
    0.532700  -0.00092;
    0.500000  -0.00405;
    0.467300  -0.00731;
    0.434740  -0.01045;
    0.402450  -0.01357;
    0.370590  -0.01637;
    0.339280  -0.01895;
    0.308660  -0.021;
    0.278860  -0.02275;
    0.250000  -0.02389;
    0.222210  -0.02475;
    0.195620  -0.025;
    0.170330  -0.02503;
    0.146450  -0.02447;
    0.124080  -0.02377;
    0.103320  -0.02246;
    0.084270  -0.0211;
    0.066990  -0.01913;
    0.051560  -0.0173;
    0.038060  -0.01481;
    0.026530  -0.01247;
    0.017040  -0.0097;
    0.009610  -0.00691;
    0.004280  -0.00436;
    0.001070  -0.002;
    0.0        0.0;
]

section = zeros(3, size(airfoil, 1))
for ic = 1:size(airfoil, 1)
    section[1,ic] = airfoil[ic,1] - 0.5
    section[2,ic] = 0
    section[3,ic] = airfoil[ic,2]
end

write_vtk("dynamic-joined-wing-simulation", assembly, history, t, scaling=1e2;
    sections = section)

# ![](../assets/dynamic-joined-wing-simulation.gif)


text:
using JuMP, LinearAlgebra, SparseArrays
using SCS
using MosekTools
const MOI = JuMP.MathOptInterface

# formulating LMI constraint with sparse matrices
#   https://github.com/JuliaOpt/JuMP.jl/issues/1950
# We are trying to solve:
# min c' x
# s.t.  A0 + A1 * x[1] + A2 * x[2] ⪰ 0

m = 2
A0 = sparse([1. 0; 0 3])
A1 = sparse([0. 1; 1 -2])
A2 = sparse([5. -1; -1 0])
B = [A1, A2]
c = [1.; 1]

# A0 = [1. 0; 0 3]
# A1 = [0. 1; 1 -2]
# A2 = [5. -1; -1 0]
# B = [A1, A2]
# c = [1.; 1]

# model = JuMP.Model(with_optimizer(SCS.Optimizer, verbose = 0));
model = JuMP.Model(with_optimizer(Mosek.Optimizer, QUIET=false))

@variable(model, x[1:m]);
@objective(model, Min, c' * x);

# formulate LMI constraint A0 + A1 x[1] + A2 x[2] in PSDCone
# works
K = A0 + sum(B[k] .* x[k] for k in 1:m)
@SDconstraint(model, con1,  Symmetric(K) >= 0);
# @constraint(model, con1,  Symmetric(K) in JuMP.PSDCone());

# throws error
# @constraint(model, con2,  A0 + sum(B[k] .* x[k] for k in 1:m) in JuMP.PSDCone());
# @SDconstraint(model, con1,  K >= 0);

JuMP.optimize!(model)
@show JuMP.objective_value(model)

@show JuMP.termination_status(model) == MOI.OPTIMAL
@show JuMP.primal_status(model) == MOI.FEASIBLE_POINT

L = JuMP.value.(K)
# L = A0 + sum(B[k] .* obj_x[k] for k in 1:m)
@show isposdef(L)


text:
# ---
# title: 1672. Richest Customer Wealth
# id: problem1672
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/richest-customer-wealth/description/>
# hidden: true
# ---
#
# You are given an `m x n` integer grid `accounts` where `accounts[i][j]` is the
# amount of money the `ith` customer has in the
# `jth` bank. Return _the **wealth** that the richest customer
# has._
#
# A customer's **wealth** is the amount of money they have in all their bank
# accounts. The richest customer is the customer that has the maximum
# **wealth**.
#
#
#
# **Example 1:**
#
#
#
#     Input: accounts = [[1,2,3],[3,2,1]]
#     Output: 6
#     **Explanation** **:**
#     1st customer has wealth = 1 + 2 + 3 = 6
#     2nd customer has wealth = 3 + 2 + 1 = 6
#     Both customers are considered the richest with a wealth of 6 each, so return 6.
#
#
# **Example 2:**
#
#
#
#     Input: accounts = [[1,5],[7,3],[3,5]]
#     Output: 10
#     **Explanation** :
#     1st customer has wealth = 6
#     2nd customer has wealth = 10
#     3rd customer has wealth = 8
#     The 2nd customer is the richest with a wealth of 10.
#
# **Example 3:**
#
#
#
#     Input: accounts = [[2,8,7],[7,1,3],[1,9,5]]
#     Output: 17
#
#
#
#
# **Constraints:**
#
#   * `m == accounts.length`
#   * `n == accounts[i].length`
#   * `1 <= m, n <= 50`
#   * `1 <= accounts[i][j] <= 100`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
#!/usr/bin/julia

tic()

using ArgParse

using JLD
using Multilane
using POMDPs
using MCTS
using DataFrames
using DataFramesMeta

s = ArgParseSettings()

@add_arg_table s begin
    "objectfile"
        help = "jld file containing all of the objects"
        nargs = 1
    "listfile"
        help = "jld file containing the stats table of the sims to be run"
        nargs = 1
    "--progress"
        help = "print progress information"
        action = :store_true
end

args = parse_args(ARGS, s)

objects = load(first(args["objectfile"]))
listfile = first(args["listfile"])
list = load(listfile)
stats = list["stats"]
metrics = get(list, "metrics", [])

sims = run_simulations(stats, objects, progress=args["progress"])
if any(isa(r, RemoteException) || nothing !== (r.exception) for r in sims)
    print_with_color(:red, """
        Error in simulations. Run the following to debug:

        ARGS = $ARGS
        include("$(@__FILE__)")
        """)
        # julia -i $(@__FILE__) $(first(args["objectfile"])) $(first(args["listfile"]))
    println()
end
fill_stats!(stats, objects, sims)

objects["stats"] = stats

dir = dirname(listfile)
b = basename(listfile)
if contains(b, "list")
    b = replace(b, "list", "results")
else
    b = string("results_", b)
end
resultsfile = joinpath(dir, b)
objects["histories"]=nothing

save(resultsfile, objects)
println("Done running $(nrow(stats)) tests.")
toc()


text:
using LightXML

if length(ARGS) < 1 error("L argument missing") end

L=parse(Int, ARGS[1])
N=L^2

# create an empty XML document
xdoc = XMLDocument()

# create & attach a graph (root) node
graph = create_root(xdoc, "GRAPH")

function addvertex!(node, id, typ, coordinates)
    v = new_child(node, "VERTEX")
    set_attributes(v, Dict("id"=>id, "type"=>typ))

    vc = new_child(v, "COORDINATE")
    add_text(vc, "$(join(coordinates .- 1, " "))")
end

function addedge!(node, src, trg, id, typ, vector)
    edge = new_child(node, "EDGE")
    set_attributes(edge, Dict("id"=>id))
    set_attributes(edge, Dict("source"=>src, "target"=>trg))
    set_attributes(edge, Dict("type"=>typ, "vector"=>join(vector, " ")))
end

# nn
sql = reshape(1:N, L, L)
up = circshift(sql,(-1,0))
right = circshift(sql,(0,-1))
down = circshift(sql,(1,0))
left = circshift(sql,(0,1))
neighbors = vcat(up[:]',right[:]',down[:]',left[:]') # colidx = site, rowidx = up right down left

# next nn (xy, yx)
ur = circshift(sql,(-1,-1))
dr = circshift(sql,(1,-1))
dl = circshift(sql,(1,1))
ul = circshift(sql,(-1,1))
Nneighbors = vcat(ur[:]',dr[:]',dl[:]',ul[:]') # colidx = site, rowidx = ur dr dl ul

# next-next nn (xx, yy)
uu = circshift(sql,(-2,0))
rr = circshift(sql,(0,-2))
dd = circshift(sql,(2,0))
ll = circshift(sql,(0,2))
NNneighbors = vcat(uu[:]',rr[:]',dd[:]',ll[:]') # colidx = site, rowidx = uu rr dd ll

# vertices
vi = 1 # vertex counter
for x in 1:L
    for y in 1:L
        global vi
        addvertex!(graph, vi, 0, [x, y])
        vi += 1
    end
end

# edges
ei = 1 # edge counter
for site in 1:N
    global ei

    # nn
    up = neighbors[1,site]
    right = neighbors[2,site]
    addedge!(graph, site, up, ei, 0, [0, 1])
    ei += 1
    addedge!(graph, site, right, ei, 0, [1, 0])
    ei +=1

    # next nn
    ur = Nneighbors[1,site]
    dr = Nneighbors[2,site]
    addedge!(graph, site, ur, ei, 1, [1, 1])
    ei += 1
    addedge!(graph, site, dr, ei, 1, [1, -1])
    ei +=1

    # next-next nn
    uu = NNneighbors[1,site]
    rr = NNneighbors[2,site]
    addedge!(graph, site, uu, ei, 2, [0, 2])
    ei += 1
    addedge!(graph, site, rr, ei, 2, [2, 0])
    ei +=1

end

set_attributes(graph, Dict("dimension"=>2, "vertices"=>vi-1, "edges"=>ei-1))

# s = string(xdoc);
# print(xdoc)

save_file(xdoc, "NNsquare_L_$(L)_W_$(L).xml")

max_stars_repo_name: crstnbr/dqmc

max_issues_repo_name: crstnbr/dqmc

max_forks_repo_name: crstnbr/dqmc


text:
module Utils

export
    second, minute, hour, day, meter, kilometer,
    KiB, MiB, GiB, TiB,
    @loop_xyz, @loop_xy, @loop_xz, @loop_yz,
    launch_config,
    cell_advection_timescale,
    TimeStepWizard, update_Δt!,
    prettytime, pretty_filesize,
    tupleit, parenttuple, datatuple, datatuples,
    validate_interval, time_to_run,
    ordered_dict_show,
    with_tracers

include("adapt_structure.jl")
include("units.jl")
include("loop_macros.jl")
include("launch_config.jl")
include("cell_advection_timescale.jl")
include("pretty_time.jl")
include("pretty_filesize.jl")
include("tuple_utils.jl")
include("output_writer_diagnostic_utils.jl")
include("ordered_dict_show.jl")
include("with_tracers.jl")

end


text:
module PlotuDiff

abstract type StringDiffAlgorithm end

include("./diffs.jl")
include("./viewer.jl")
include("./myers.jl")

"""
Compute the difference between two input strings and return the result of diffs wrapped
"""
function stringdiff(string1::AbstractString, string2::AbstractString; alg=Myers(), viewer=SimpleDiffViewer)
	stringdiff(alg, string1, string2) |> viewer
end

stringdiff_without_viewer(args...; kwargs...) = stringdiff(args...; kwargs..., viewer=identity)

export stringdiff, merge_diffs, SimpleDiffViewer, LineViewer, SimpleLineDiffViewer

end


text:
function relaxed_initial_state_TRI(mdp::NoCrashProblem, LastState::MLState, MKZAction::MLAction, DSRC_interval::Float64
                             ,rng=MersenneTwister(rand(UInt32)))

    pp = mdp.dmodel.phys_param
    #sim = HistoryRecorder(max_steps=steps, rng=rng)
    #policy = solve(solver, mdp)
    #policy = RandomPolicy(mdp)
    CurrentState = simulate_TRI(mdp, MKZAction,LastState, DSRC_interval,rng)
    #println("simulate done")
    return CurrentState
end

function simulate_TRI{S,A}(mdp::MDP{S,A}, MKZAction::MLAction,LastState::MLState, DSRC_interval::Float64,rng::AbstractRNG)
    sp = generate_sTRI(mdp, LastState, MKZAction, DSRC_interval,rng)
    return sp
end

function generate_sTRI(mdp::NoCrashProblem, s::MLState, a::MLAction, DSRC_interval::Float64, rng::AbstractRNG)

    @if_debug dbg_rng = copy(rng)

    sp::MLState=create_state(mdp)
    pp = mdp.dmodel.phys_param
    dt = pp.dt
    nb_cars = length(s.cars)
    resize!(sp.cars, nb_cars)
    #sp.x = 0
    sp.terminal = s.terminal

    ## Calculate deltas ##
    #====================#

    dxs = Array{Float64}(nb_cars)
    dvs = Array{Float64}(nb_cars)
    dys = Array{Float64}(nb_cars)
    lcs = Array{Float64}(nb_cars)

    # agent
    dvs[1] = a.acc*dt
    dxs[1] = s.cars[1].vel*dt + a.acc*dt^2/2.
    lcs[1] = a.lane_change
    dys[1] = a.lane_change*dt

    for i in 2:nb_cars
        neighborhood = get_neighborhood(pp, s, i)

        behavior = s.cars[i].behavior

        acc = gen_accel(behavior, mdp.dmodel, s, neighborhood, i, rng)
        dvs[i] = dt*acc
        dxs[i] = (s.cars[i].vel + dvs[i]/2.)*dt

        lcs[i] = gen_lane_change(behavior, mdp.dmodel, s, neighborhood, i, rng)
        dys[i] = lcs[i] * dt
    end

    ## Consistency checking ##
    #========================#

    # first prevent lane changes into each other
    changers = IntSet()
    for i in 1:nb_cars
        if lcs[i] != 0
            push!(changers, i)
        end
    end
    sorted_changers = sort!(collect(changers), by=i->s.cars[i].x, rev=true) # this might be slow because anonymous functions are slow
    # from front to back

    if length(sorted_changers) >= 2 #something to compare
        # iterate through pairs
        iter_state = start(sorted_changers)
        j, iter_state = next(sorted_changers, iter_state)
        while !done(sorted_changers, iter_state)
            i = j
            j, iter_state = next(sorted_changers, iter_state)
            car_i = s.cars[i] # front
            car_j = s.cars[j] # back

            # check if they are both starting to change lanes on this step
            if isinteger(car_i.y) && isinteger(car_j.y)

                # check if they are near each other lanewise
                if abs(car_i.y - car_j.y) <= 2.0

                    # make sure there is a conflict longitudinally
                    # if car_i.x - car_j.x <= pp.l_car || car_i.x + dxs[i] - car_j.x + dxs[j] <= pp.l_car
                    # if car_i.x - car_j.x <= mdp.dmodel.appear_clearance # made more conservative on 8/19
                    # if car_i.x - car_j.x <= get_idm_s_star(car_j.behavior.p_idm, car_j.vel, car_j.vel-car_i.vel) # upgraded to sstar on 8/19
                    ivp = car_i.vel + dt*dvs[i]
                    jvp = car_j.vel + dt*dvs[j]
                    ixp = car_i.x + dt*(car_i.vel + ivp)/2.0
                    jxp = car_j.x + dt*(car_j.vel + jvp)/2.0
                    n_max_acc_p = nullable_max_safe_acc(ixp-jxp-pp.l_car, jvp, ivp, pp.brake_limit,dt)
                    if ixp - jxp <= pp.l_car || car_i.x - car_j.x <= pp.l_car || isnull(n_max_acc_p) || get(n_max_acc_p) < -pp.brake_limit

                        # check if they are moving towards each other
                        # if dys[i]*dys[j] < 0.0 && abs(car_i.y+dys[i] - car_j.y+dys[j]) < 2.0
                        if true # prevent lockstepping 8/19 (doesn't prevent it that well)

                            # make j stay in his lane
                            dys[j] = 0.0
                            lcs[j] = 0.0
                        end
                    end
                end
            end
        end
    end

    # second, prevent cars hitting each other due to noise
    sorted = sort!(collect(1:length(s.cars)), by=i->s.cars[i].x, rev=true)

    if length(sorted) >= 2 #something to compare
        # iterate through pairs
        iter_state = start(sorted)
        j, iter_state = next(sorted, iter_state)
        while !done(sorted, iter_state)
            i = j
            j, iter_state = next(sorted, iter_state)
            if j == 1
                continue # don't check for the ego since the ego does not have noise
            end
            car_i = s.cars[i]
            car_j = s.cars[j]

            # check if they overlap longitudinally
            if car_j.x + dxs[j] > car_i.x + dxs[i] - pp.l_car

                # check if they will be in the same lane
                if occupation_overlap(car_i.y + dys[i], car_j.y + dys[j])
                    # warn and nudge behind
                    @if_debug begin
                        println("Conflict because of noise: front:$i, back:$j")
                        Gallium.@enter generate_s(mdp, s, a, dbg_rng)
                    end
                    if i == 1
                         warn("Car nudged because noise would cause a crash (ego in front).")
                        #error("Car nudged because noise would cause a crash (ego in front).")
                    else
                         warn("Car nudged because noise would cause a crash.")
                        #error("Car nudged because noise would cause a crash.")
                    end
                    dxs[j] = car_i.x + dxs[i] - car_j.x - 1.01*pp.l_car
                    dvs[j] = 2.0*(dxs[j]/dt - car_j.vel)
                end
            end
        end
    end

    ## Dynamics and Exits ##
    #======================#

    exits = IntSet()
    for i in 1:nb_cars
        car = s.cars[i]
        xp = car.x + (dxs[i] - dxs[1])
        yp = car.y + dys[i]
        # velp = max(min(car.vel + dvs[i],pp.v_max), pp.v_min)
        velp = max(car.vel + dvs[i], 0.0) # removed speed limits on 8/13
        # note lane change is updated above

        if dvs[i]/dt < -mdp.dmodel.brake_terminate_thresh
            sp.terminal = Nullable{Symbol}(:brake)
        end

        # check if a lane was crossed and snap back to it
        if isinteger(car.y)
            # prevent a multi-lane change in a single timestep
            if abs(yp-car.y) > 1.
                yp = car.y + sign(dys[i])
            end
        else # car.y is not an integer
            if floor(yp) >= ceil(car.y)
                yp = ceil(car.y)
            end
            if ceil(yp) <= floor(car.y)
                yp = floor(car.y)
            end
        end

        # if yp < 1.0 || yp > pp.nb_lanes
        #     @show i
        #     @show yp
        #     println("mdp = $mdp")
        #     println("s = $s")
        #     println("a = $a")
        # end
        #print("yp")
        #println(yp)
        @assert yp >= 1.0 && yp <= pp.nb_lanes

        if xp < 0.0 || xp >= pp.lane_length
            push!(exits, i)
        else
            sp.cars[i] = CarState(xp, yp, velp, lcs[i], car.behavior, s.cars[i].id)
        end
    end

    next_id = maximum([c.id for c in s.cars]) + 1

    deleteat!(sp.cars, exits)
    nb_cars -= length(exits)

    ## Generate new cars ##
    #=====================#

    if nb_cars < mdp.dmodel.nb_cars && rand(rng) <= mdp.dmodel.p_appear

        behavior = rand(rng, mdp.dmodel.behaviors)
        vel = typical_velocity(behavior) + randn(rng)*mdp.dmodel.vel_sigma

        clearances = Array{Float64}(pp.nb_lanes)
        fill!(clearances, Inf)
        closest_cars = Array{Int}(pp.nb_lanes)
        fill!(closest_cars, 0)
        sstar_margins = Array{Float64}(pp.nb_lanes)
        if vel > s.cars[1].vel
            # put at back
            # sstar is the sstar of the new guy
            for i in 1:length(s.cars)
                lowlane, highlane = occupation_lanes(s.cars[i].y, lcs[i])
                back = s.cars[i].x - pp.l_car
                if back < clearances[lowlane]
                    clearances[lowlane] = back
                    closest_cars[lowlane] = i
                end
                if back < clearances[highlane]
                    clearances[highlane] = back
                    closest_cars[highlane] = i
                end
            end
            for j in 1:pp.nb_lanes
                other = closest_cars[j]
                if other == 0
                    sstar = 0
                else
                    sstar = get_idm_s_star(behavior.p_idm, vel, vel-s.cars[other].vel)
                end
                sstar_margins[j] = clearances[j] - sstar
            end
        else
            for i in 1:length(s.cars)
                lowlane, highlane = occupation_lanes(s.cars[i].y, lcs[i])
                front = pp.lane_length - (s.cars[i].x + pp.l_car) # l_car is half the length of the old car plus half the length of the new one
                if front < clearances[lowlane]
                    clearances[lowlane] = front
                    closest_cars[lowlane] = i
                end
                if front < clearances[highlane]
                    clearances[highlane] = front
                    closest_cars[highlane] = i
                end
            end
            for j in 1:pp.nb_lanes
                other = closest_cars[j]
                if other == 0
                    sstar = 0
                else
                    sstar = get_idm_s_star(s.cars[other].behavior.p_idm,
                                           s.cars[other].vel,
                                           s.cars[other].vel-vel)
                end
                sstar_margins[j] = clearances[j] - sstar
            end
        end

        margin, lane = findmax(sstar_margins)

        if margin > 0.0
            if vel > s.cars[1].vel
                # at back
                push!(sp.cars, CarState(0.0, lane, vel, 0.0, behavior, next_id))
            else
                push!(sp.cars, CarState(pp.lane_length, lane, vel, 0.0, behavior, next_id))
            end
        end
    end

    if mdp.dmodel.lane_terminate && sp.cars[1].y == mdp.rmodel.target_lane
        sp.terminal = Nullable{Any}(:lane)
    elseif sp.x > mdp.dmodel.max_dist
        sp.terminal = Nullable{Any}(:distance)
    end

    # sp.crashed = is_crash(mdp, s, sp, warning=false)
    # sp.crashed = false

    @assert sp.cars[1].x == s.cars[1].x # ego should not move
    sp.x = s.x + sp.cars[1].vel*0.75
    sp.t = s.t
    ca = control_action(mdp,s)
    if s.Control_Signal && (a == ca)
       sp.Control_Signal = true
    else
       sp.Control_Signal = false
    end

    if s.x > 2500 && isinteger(s.cars[1].y) && s.t == 0.0
        sp.Control_Signal = true
        sp.t = 1.0
    end

    #if s.cars[1].y == 4.0 && s.t == 0.0
    #    sp.Control_Signal = true
    #    sp.t = 1.0
    #end
    return sp
end

max_stars_repo_name: zhcao/HighwayExit

max_issues_repo_name: zhcao/HighwayExit

max_forks_repo_name: zhcao/HighwayExit


text:
# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Script with functions to create τ_hat_Z, τ_hat_F from tariff data
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

"""
    country_sector_tariffs(ctry_one::Vector, ctry_two::Vector, tariffs_one::DataFrame, tariffs_two::DataFrame, industries::Vector, N::Integer, S::Integer)

The function matches industry bilateral tariffs which can then be used to create a matrix with the function 'create_τ_hat'. Notice that there are
    two country and tariff vector, i.e. asymmetric bilateral trade costs are possible. If the model only allows symmetric bilateral trade costs
    simply use the same DataFrame for both the arguments 'tariffs_one' and 'tariffs_two'.

# Arguments
- `ctry_one::Vector`: a vector of countries (ISO3 codes) on which the tarrifs of vector 'tariffs_one' are levied on by countries in vector 'ctry_two' .
- `ctry_two::Vector`: a vector of countries (ISO3 codes) on which the tarrifs of vector 'tariffs_two' are levied on by countries in vector 'ctry_one' .
- `tariffs_one::DataFrame`: a DataFrame carrying information on industry tariffs levied on 'ctry_two' by 'ctry_one'.
- `tariffs_two::DataFrame`: a DataFrame carrying information on industry tariffs levied on 'ctry_one' by 'ctry_two'.
- `industries::Vector{String}`: a vector of industries to match tariff and country data.
- `N::Integer`: number of origin/destination countries.
- `S::Integer`: number of origin/destination industries.

# Output
- `τ_new::DataFrame`: a DataFrame containing all affected country-industry bilateral trade costs.

# Examples
```julia-repl
julia> τ_new = country_sector_tariffs(ctry_exiting_EU, ctry_remaining_EU, df_tariffs, df_tariffs, industries, N, S)
```
"""

function country_sector_tariffs(ctry_one::Vector, ctry_two::Vector, tariffs_one::DataFrame, tariffs_two::DataFrame, industries::Vector, N::Integer, S::Integer)

    # uses countries in vector "ctry_one" as exporters which pay tariffs "tariffs_one" to importer countries "ctry_two"
    # as the importing ctry levy the tariffs depending on the product (sector) of the exporters
    # we need to match with the exporter with the sector of the exporter!
    exporter_ctry = ctry_one
    importer_ctry = ctry_two

    exporter_ID = repeat(exporter_ctry, inner=S) .* "__" .* repeat(industries, outer=length(exporter_ctry))
    importer_ID = repeat(importer_ctry, inner=S) .* "__" .* repeat(industries, outer=length(importer_ctry))

    exporter = repeat(exporter_ID, inner=length(importer_ID))
    importer = repeat(importer_ID, outer=length(exporter_ID))
    exporter_sector = reduce(vcat, permutedims.(split.(exporter, "__")))[:, 2] # to match with tariff data (on sectoral level)
    importer_ctry_F = reduce(vcat, permutedims.(split.(importer, "__")))[:, 1] # to match with final demand country

    setup_one = DataFrame([exporter importer exporter_sector importer_ctry_F], [:exporter, :importer, :exporter_sector, :importer_ctry])
    df_one = rename(tariffs_one, :industry_WIOD_2016 => :exporter_sector)

    τ_one = leftjoin(setup_one, df_one, on=:exporter_sector)

    # -------

    # uses countries in vector "ctry_two" as exporters which pay tariffs "tariffs_two" to importer countries "ctry_one"
    exporter_ctry = ctry_two
    importer_ctry = ctry_one

    exporter_ID = repeat(exporter_ctry, inner=S) .* "__" .* repeat(industries, outer=length(exporter_ctry))
    importer_ID = repeat(importer_ctry, inner=S) .* "__" .* repeat(industries, outer=length(importer_ctry))

    exporter = repeat(exporter_ID, inner=length(importer_ID))
    importer = repeat(importer_ID, outer=length(exporter_ID))
    exporter_sector = reduce(vcat, permutedims.(split.(exporter, "__")))[:, 2] # to match with tariff data (on sectoral level)
    importer_ctry_F = reduce(vcat, permutedims.(split.(importer, "__")))[:, 1] # to match with final demand

    setup_two = DataFrame([exporter importer exporter_sector importer_ctry_F], [:exporter, :importer, :exporter_sector, :importer_ctry])
    df_two = rename(tariffs_two, :industry_WIOD_2016 => :exporter_sector)

    τ_two = leftjoin(setup_two, df_two, on=:exporter_sector)

    # -------

    τ_new = vcat(τ_one, τ_two)

    return τ_new

end

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

"""
    create_τ_hat(df::DataFrame, scenario::Symbol, countries::Vector, N::Integer, S::Integer)

The function maps the industry bilateral country-industry tariffs supplied by the function 'country_sector_tariffs' in the correct cells
    for to produce the NS×NS matrix τ_hat_Z, and NS×N matrix τ_hat_F used for computing a counterfactual. Notice bilateral trade costs can
    vary between intermediate inputs and final demand.

# Arguments
- `df::DataFrame`: a DataFrame providing the counterfactual in form of bilateral country-industry tariffs.
- `scenario_Z::Symbol`: specifies the tariff scenario for intermediate input trade costs.
- `scenario_F::Symbol`: specifies the tariff scenario for final demand trade costs.
- `countries::Vector`: a vector of all countries (ISO3 codes) used in the IO table.
- `N::Integer`: number of origin/destination countries.
- `S::Integer`: number of origin/destination industries.

# Output
- `τ_hat_Z::Matrix{Float64}`: NS×NS, origin country-industry destination country-industry intermediate input trade cost (changes) matrix τ_hat_Z.
- `τ_hat_F::Matrix{Float64}`: NS×N, origin country-industry destination country final demand trade cost (changes) matrix τ_hat_Z.

# Examples
```julia-repl
julia> τ_hat_Z, τ_hat_F = create_τ_hat(df_affected, :avg_tariff_Z_soft, :avg_tariff_F_soft, ctry_names_2013, N, S)
```
"""

function create_τ_hat(df::DataFrame, scenario_Z::Symbol, scenario_F::Symbol, countries::Vector, N::Integer, S::Integer)

    col_names = repeat(countries, inner=S) .* "__" .* repeat(industries, outer=N) # NS×1

    # τ_hat_Z
    df_Z = df[:, [:exporter, :importer, :exporter_sector, :importer_ctry, scenario_Z]]
    rename!(df_Z, scenario_Z => :tariffs)

    τ_hat_Z = DataFrame([col_names ones(N*S, N*S)], ["exporter"; col_names])
    τ_hat_Z = stack(τ_hat_Z, Not(:exporter))
    rename!(τ_hat_Z, :variable => :importer)

    τ_hat_Z = leftjoin(τ_hat_Z, df_Z, on=[:exporter, :importer])
    τ_hat_Z.tariffs = ifelse.(ismissing.(τ_hat_Z.tariffs), τ_hat_Z.value, τ_hat_Z.tariffs)
    sort!(τ_hat_Z, [:exporter, :importer]) # sort to have same sorting as IO Table (thats why we use "ZoR" instead of "RoW")
    τ_hat_Z = τ_hat_Z[:, [:exporter, :importer, :tariffs]]
    τ_hat_Z = unstack(τ_hat_Z, :importer, :tariffs, allowduplicates=true)

    τ_hat_Z = Matrix(convert.(Float64, τ_hat_Z[:, 2:end])) # NS×NS

    # -------

    # τ_hat_F
    df_F = df[:, [:exporter, :importer, :exporter_sector, :importer_ctry, scenario_F]]
    rename!(df_F, scenario_F => :tariffs)

    τ_hat_F = DataFrame([col_names ones(N*S, N)], ["exporter"; countries])
    τ_hat_F = stack(τ_hat_F, Not(:exporter))
    rename!(τ_hat_F, :variable => :importer_ctry)

    τ_hat_F = leftjoin(τ_hat_F, df_F, on=[:exporter, :importer_ctry])
    τ_hat_F.tariffs = ifelse.(ismissing.(τ_hat_F.tariffs), τ_hat_F.value, τ_hat_F.tariffs)
    sort!(τ_hat_F, [:exporter, :importer_ctry]) # sort to have same sorting as IO Table (thats why we use "ZoR" instead of "RoW")
    τ_hat_F = τ_hat_F[:, [:exporter, :importer_ctry, :tariffs]]
    τ_hat_F = unstack(τ_hat_F, :importer_ctry, :tariffs, allowduplicates=true) # for some reason it gives out duplicate error?

    τ_hat_F = Matrix(convert.(Float64, τ_hat_F[:, 2:end])) # NS×N

    return τ_hat_Z, τ_hat_F
end


text:
export RedisBitVector

struct RedisBitVector <: AbstractRedisCollection
    conn::AbstractRedisConnection
    key::String
end

function getindex(rbv::RedisBitVector, index::Int64)
    exec(rbv.conn, "getbit", rbv.key, index) |> Bool
end

function setindex!(rbv::RedisBitVector, value::Bool, key::Int64)
    exec(rbv.conn, "setbit", key, value ? "1" : "0")
end

function isnull(rbv::RedisBitVector)
    exec(rbv.conn, "exists", rbv.key) |> Bool
end

function length(rbv::RedisBitVector)
    8 * exec(rbv.conn, "strlen", rbv.key)
end

function sum(rbv::RedisBitVector)
    exec(rbv.conn, "bitcount", rbv.key)
    wait(rbv.conn)::Int64
end

"this methods should always be used in the form of `rbv1 |= rbv2`"
function (|)(rbv1::RedisBitVector, rbv2::RedisBitVector)
    rbv1.conn == rbv2.conn || throw(ArgumentError("tow RedisBitVectors must have the same RedisConnection"))
    exec(rbv1.conn, "bitop", "or", rbv1.key, rbv2.key)
    rbv1
end

"this methods should always be used in the form of `rbv1 &= rbv2`"
function (&)(rbv1::RedisBitVector, rbv2::RedisBitVector)
    rbv1.conn == rbv2.conn || throw(ArgumentError("tow RedisBitVectors must have the same RedisConnection"))
    exec(rbv1.conn, "bitop", "and", rbv1.key, rbv2.key)
    rbv1
end

"this methods should always be used in the form of `rbv1 \$= rbv2`"
function ($)(rbv1::RedisBitVector, rbv2::RedisBitVector)
    rbv1.conn == rbv2.conn || throw(ArgumentError("tow RedisBitVectors must have the same RedisConnection"))
    exec(rbv1.conn, "bitop", "xor", rbv1.key, rbv2.key)
    rbv1
end

"this methods should always be used in the form of `rbv = ~rbv`"
function (~)(rbv::RedisBitVector)
    exec(rbv.conn, "bitop", "not", rbv.key)
    rbv1
end

"Caution: Redis consider the right of the string as padded with zeros, so it may act unexpectly when finding `false`"
function findfirst(rbv::RedisBitVector, v::Bool=true)
    exec(rbv.conn, "bitpos", rbv.key, v ? "1" : "0")
    1 + wait(rbv.conn)::Int64
end


text:
# RESERVED_WORDS, identifier, makeidentifier were originally part of DataFrames v0.20
const RESERVED_WORDS = Set(["local", "global", "export", "let",
    "for", "struct", "while", "const", "continue", "import",
    "function", "if", "else", "try", "begin", "break", "catch",
    "return", "using", "baremodule", "macro", "finally",
    "module", "elseif", "end", "quote", "do"])

function identifier(s::AbstractString)
    s = Unicode.normalize(s)
    if !Base.isidentifier(s)
        s = makeidentifier(s)
    end
    Symbol(in(s, RESERVED_WORDS) ? "_"*s : s)
end

function makeidentifier(s::AbstractString)
    (iresult = iterate(s)) === nothing && return "x"

    res = IOBuffer(zeros(UInt8, sizeof(s)+1), write=true)

    (c, i) = iresult
    under = if Base.is_id_start_char(c)
        write(res, c)
        c == '_'
    elseif Base.is_id_char(c)
        write(res, 'x', c)
        false
    else
        write(res, '_')
        true
    end

    while (iresult = iterate(s, i)) !== nothing
        (c, i) = iresult
        if c != '_' && Base.is_id_char(c)
            write(res, c)
            under = false
        elseif !under
            write(res, '_')
            under = true
        end
    end

    return String(take!(res))
end


text:
@testset "005" begin
    @test tdpe.e005(10) == 2520
end

max_stars_repo_name: mithodin/ttpe

max_issues_repo_name: mithodin/ttpe

max_forks_repo_name: mithodin/ttpe


text:
using Turing
using Distributions
using Base.Test

# Test the @predict macro on a deterministic model.

@model test begin
  x = 0
  @predict x
  x = 1
  y = 2
  @predict y
end

s = SMC(10)
p = PG(2,5)

res = sample(test, s)

@test res[:x][1] == 0
@test res[:y][1] == 2

res = sample(test, p)

@test res[:x][1] == 0
@test res[:y][1] == 2


text:
function avg2(x)
if (typeof(x[1])==Float64) | (typeof(x[1])==Float32)
        return (x[1:2:end] + x[2:2:end])/2;
else
        return x[1:2:end]
end

end

function avg2_2d(x)

return (x[1:2:end,:] + x[2:2:end,:])/2;

end

function avg2_df(x)

ncol = Integer(size(x)[2]);
nrow = Integer(size(x)[1]/2);
ans = deepcopy(x[1:nrow,:]);

for j in 1:ncol
        ans[:,j] = avg2(x[:,j]);
end

return ans;

end

function convert_sim(x)
a = avg2_df(x[1]);
b = avg2_2d(x[2]);
c = avg2_2d(x[3]);
d = avg2_2d(x[4]);
d2 = avg2(x[5]);

return a,b,c,d,d2
end

function get_diurnal(x,navg)
	if (typeof(x[1])==Float64) | (typeof(x[1])==Float32)
			return mean(reshape(x,(navg,:)),dims=2)[:,1];
	else
			return x[1:navg]
	end
end

function get_diurnal_summer(x,navg)
	if (typeof(x[1])==Float64) | (typeof(x[1])==Float32)
		day_reshape = reshape(x,(navg,:));
		day_list = collect(1:(size(day_reshape)[2])) .% 365;
		summer_mask = (day_list .>= 150) .* (day_list .< 250);
		summer_only = day_reshape[:,summer_mask];
		return mean(summer_only,dims=2)[:,1];
	else
			return x[1:navg]
	end
end

function get_daily(x,navg)
	if (typeof(x[1])==Float64) | (typeof(x[1])==Float32)
			return mean(reshape(x,(navg,:)),dims=1)[1,:];
	else
			return x[1:navg:end]
	end
end

function get_daily2d(x,navg)
	return mean(reshape(x,(navg,:,size(x)[2])),dims=1)[1,:,:];
end