text:
"""
branch_mostinfeasible(m, node, disc2var_idx)

Get the index of an integer variable which is currently continuous which is most unintegral.
(nearest to *.5)
"""
function branch_mostinfeasible(m, node, disc2var_idx)
    x = node.solution
    idx = 0
    max_diff = 0
    for i in disc2var_idx
        diff = abs(x[i]-round(x[i]))
        if diff > max_diff
            idx = i
            max_diff = diff
        end
    end
    return idx
end

"""
init_strong_restart!(node, opts, var_idx, disc_var_idx, l_nd, r_nd, reasonable_disc_vars, infeasible_disc_vars,
 left_node, right_node, strong_restart, max_gain_var, atol, var_type, disc2var_idx)

Tighten the bounds for the node and check if there are variables that need to be checked for a restart.
"""
function init_strong_restart!(node, opts, var_idx, disc_var_idx, l_nd, r_nd,
                                reasonable_disc_vars, infeasible_disc_vars,
                                left_node, right_node, strong_restart, max_gain_var,
                                atol, var_type, disc2var_idx)
    restart = false
    set_to_last_var = false
    need_to_resolve = false

    # set the bounds directly for the node
    # also update the best bound and the solution
    if !state_is_optimal(l_nd.relaxation_state; allow_almost=opts.allow_almost_solved)
        # the solution shouldn't be updated when the current max_gain_var is then type correct
        if max_gain_var == 0 || !is_type_correct(r_nd.solution[max_gain_var],var_type[max_gain_var],atol)
            # if all variables reasonable_disc_vars are type correct => just continue with the tree
            # no restarts and no change in bounds
            if !all_reasonable_type_correct(r_nd.solution, disc2var_idx, reasonable_disc_vars, atol)
                node.l_var[var_idx] = ceil(node.solution[var_idx])
                node.best_bound = r_nd.best_bound
                node.solution = r_nd.solution
                need_to_resolve = true
                push!(infeasible_disc_vars, disc_var_idx)
                strong_restart && (restart = true)
            end
        end
    else
        # the solution shouldn't be updated when the current max_gain_var is then type correct
        if max_gain_var == 0 || !is_type_correct(l_nd.solution[max_gain_var],var_type[max_gain_var],atol)
            if !all_reasonable_type_correct(l_nd.solution, disc2var_idx, reasonable_disc_vars, atol)
                node.u_var[var_idx] = floor(node.solution[var_idx])
                node.best_bound = l_nd.best_bound
                node.solution = l_nd.solution
                need_to_resolve = true
                push!(infeasible_disc_vars, disc_var_idx)
                strong_restart && (restart = true)
            end
        end
    end

    if length(reasonable_disc_vars) == length(infeasible_disc_vars)
        # basically branching on the last infeasible variable
        set_to_last_var = true
        restart = false
    end
    return need_to_resolve, restart, infeasible_disc_vars, set_to_last_var
end

"""
    branch_strong_on!(m,opts,step_obj,reasonable_disc_vars, disc2var_idx, strong_restart, counter)

Try to branch on a few different variables and choose the one with highest obj_gain.
Update obj_gain for the variables tried and average the other ones.
"""
function branch_strong_on!(m,opts,step_obj,
    reasonable_disc_vars, disc2var_idx, strong_restart, counter)

    function set_temp_gains!(gains, gain_l, gain_r, disc_var_idx)
        if !isinf(gain_l)
            gains.minus[disc_var_idx] = gain_l
            gains.minus_counter[disc_var_idx] = 1
        else
            gains.inf_counter[disc_var_idx] = 1
        end
        if !isinf(gain_r)
            gains.plus[disc_var_idx] = gain_r
            gains.plus_counter[disc_var_idx] = 1
        else
            gains.inf_counter[disc_var_idx] = 1
        end

        if !isinf(gain_l) && !isinf(gain_r)
            gains.inf_counter[disc_var_idx] -= 1
        end
    end

    function get_current_gains(opts, node, l_nd, r_nd)
        gain_l = sigma_minus(opts, node, l_nd, node.solution[node.var_idx])
        gain_r = sigma_plus( opts, node,  r_nd, node.solution[node.var_idx])

        if isinf(gain_l)
            gain = gain_r
        elseif isinf(gain_r)
            gain = gain_l
        else
            gain = (gain_l+gain_r)/2
        end
        if isnan(gain)
            gain = Inf
        end
        return gain_l, gain_r, gain
    end

    function add_strong_step(strong_step,l_nd,r_nd; restart=false)
        strong_step.l_relaxation_state = l_nd.relaxation_state
        strong_step.r_relaxation_state = r_nd.relaxation_state
        strong_step.init_restart = restart
        push!(step_obj.strong_branching, strong_step)
    end

    strong_time = time()

    node = step_obj.node
    strong_restarts = -1
    restart = true
    status = :Normal
    node = step_obj.node
    infeasible_disc_vars = zeros(Int64,0)
    gains = init_gains(m.num_disc_var)

    left_node, right_node = nothing, nothing
    atol = opts.atol

    need_to_resolve = false
    opts.debug && (step_obj.strong_branching = [])

    max_gain = -Inf # then one is definitely better
    max_gain_var = 0
    max_gain_disc_var = 0

    allow_almost = opts.allow_almost_solved
    while restart
        strong_restarts += 1 # is init with -1
        restart = false
        for disc_var_idx in reasonable_disc_vars
            # don't rerun if the variable has already one infeasible node
            if disc_var_idx in infeasible_disc_vars
                continue
            end
            var_idx = disc2var_idx[disc_var_idx]
            step_obj.var_idx = var_idx
            u_b, l_b = node.u_var[var_idx], node.l_var[var_idx]
            # don't rerun if bounds are exact or is type correct
            if isapprox(u_b,l_b; atol=atol) || is_type_correct(node.solution[var_idx],m.var_type[var_idx], atol)
                continue
            end

            # branch on the current variable and get the corresponding children
            l_nd,r_nd = branch!(m, opts, step_obj, counter, disc2var_idx; temp=true)
            opts.debug && (strong_step = new_default_strong_branch_step_obj(var_idx))

            # no current restart => we can set max_gain and variable
            gain_l, gain_r, gain = get_current_gains(opts, node, l_nd, r_nd)

            if !state_is_optimal(l_nd.relaxation_state; allow_almost=allow_almost) &&
                !state_is_optimal(r_nd.relaxation_state; allow_almost=allow_almost) && counter == 1
                # TODO: Might be Error/UserLimit instead of infeasible
                status = :Infeasible
                left_node = l_nd
                right_node = r_nd
                opts.debug && add_strong_step(strong_step,l_nd,r_nd)
                break
            end

            # check if one part is infeasible => update bounds & restart if strong restart is true
            if !state_is_optimal(l_nd.relaxation_state; allow_almost=allow_almost) ||
                !state_is_optimal(r_nd.relaxation_state; allow_almost=allow_almost)
                if !state_is_optimal(l_nd.relaxation_state; allow_almost=allow_almost) && !state_is_optimal(r_nd.relaxation_state; allow_almost=allow_almost)
                    # TODO: Might be Error/UserLimit instead of infeasible
                    status = :PartlyInfeasible
                    left_node = l_nd
                    right_node = r_nd
                    opts.debug && add_strong_step(strong_step,l_nd,r_nd)
                    break
                end
                need_to_resolve, restart, new_infeasible_disc_vars, set_to_last_var = init_strong_restart!(node, opts, var_idx, disc_var_idx, l_nd, r_nd, reasonable_disc_vars, infeasible_disc_vars, left_node, right_node, strong_restart, max_gain_var, opts.atol, m.var_type, disc2var_idx)
                infeasible_disc_vars = new_infeasible_disc_vars

                # only variables where one branch is infeasible => no restart and break
                if set_to_last_var
                    max_gain_var = var_idx
                    max_gain_disc_var = disc_var_idx
                    left_node = l_nd
                    right_node = r_nd
                    restart = false
                    need_to_resolve = false
                    gain_l, gain_r, gain = get_current_gains(opts, node, l_nd, r_nd)
                    max_gain = gain
                    set_temp_gains!(gains, gain_l, gain_r, disc_var_idx)
                    opts.debug && add_strong_step(strong_step,l_nd,r_nd)
                    break
                end

                if restart && time()-strong_time > opts.strong_branching_time_limit
                    restart = false
                end
            end

            # don't update maximum if restart
            if gain > max_gain && !restart
                max_gain = gain
                max_gain_var = var_idx
                max_gain_disc_var = disc_var_idx
                left_node = l_nd
                right_node = r_nd
                # we are changing the bounds if one branch is infeasible then the current
                # selection might not work anymore => we have to resolve it later
                need_to_resolve = false
            end
            set_temp_gains!(gains, gain_l, gain_r, disc_var_idx)
            opts.debug && add_strong_step(strong_step,l_nd,r_nd;restart=restart)
            if restart
                break
            end

            if time()-m.start_time >= opts.time_limit
                break
            end

            if time()-strong_time >= opts.strong_branching_time_limit && !isinf(max_gain)
                @info "Breaking out of strong branching as the time limit of $(opts.strong_branching_time_limit) seconds got reached."
                break
            end
        end
    end

    # check if need to resolve (not if Local/Global Infeasible)
    if need_to_resolve && status == :Normal
        status = :Resolve
    end

    return status, max_gain_var, left_node, right_node, gains, strong_restarts
end

"""
branch_strong!(m,opts,disc2var_idx,step_obj,counter)

Try to branch on a few different variables and choose the one with highest obj_gain.
Update obj_gain for the variables tried and average the other ones.
"""
function branch_strong!(m,opts,disc2var_idx,step_obj,counter)
    node = step_obj.node

    # generate an of variables to branch on
    num_strong_var = Int(round((opts.strong_branching_perc/100)*m.num_disc_var))
    # if smaller than 2 it doesn't make sense
    num_strong_var = num_strong_var < 2 ? 2 : num_strong_var

    # get reasonable candidates (not type correct and not already perfectly bounded)
    disc_vars = m.num_disc_var
    atol = opts.atol
    reasonable_disc_vars = get_reasonable_disc_vars(node, m.var_type, disc_vars,  disc2var_idx, atol)
    if num_strong_var < length(reasonable_disc_vars)
        shuffle!(reasonable_disc_vars)
        reasonable_disc_vars = reasonable_disc_vars[1:num_strong_var]
    end

    # compute the gain for each reasonable candidate and choose the highest
    status, max_gain_var, left_node, right_node, gains, strong_restarts = branch_strong_on!(m,opts,step_obj,
        reasonable_disc_vars, disc2var_idx, opts.strong_restart, counter)

    step_obj.obj_gain += gains

    if status != :Resolve
        step_obj.l_nd = left_node
        step_obj.r_nd = right_node

        # set the variable to branch (best gain)
        node.state = :Done
        node.var_idx = max_gain_var
    end

    @assert max_gain_var != 0 || status == :PartlyInfeasible || status == :Infeasible || node.state == :Infeasible
    return status, max_gain_var, strong_restarts
end

function branch_reliable!(m,opts,step_obj,disc2var_idx,gains,counter)
    idx = 0
    node = step_obj.node
    mu = opts.gain_mu
    reliability_param = opts.reliability_branching_threshold
    reliability_perc = opts.reliability_branching_perc
    num_strong_var = Int(round((reliability_perc/100)*m.num_disc_var))
    # if smaller than 2 it doesn't make sense
    num_strong_var = num_strong_var < 2 ? 2 : num_strong_var

    gmc_r = gains.minus_counter .< reliability_param
    gpc_r = gains.plus_counter  .< reliability_param

    strong_restarts = 0
    reasonable_disc_vars = []
    atol = opts.atol
    for i=1:length(gmc_r)
        if gmc_r[i] || gpc_r[i]
            idx = disc2var_idx[i]
            u_b = node.u_var[idx]
            l_b = node.l_var[idx]
            if isapprox(u_b,l_b; atol=atol) || is_type_correct(node.solution[idx],m.var_type[idx],atol)
                continue
            end
            push!(reasonable_disc_vars,i)
        end
    end
    if length(reasonable_disc_vars) > 0
        unrealiable_idx = sortperm(gains.minus_counter[reasonable_disc_vars])
        reasonable_disc_vars = reasonable_disc_vars[unrealiable_idx]
        num_reasonable = num_strong_var < length(reasonable_disc_vars) ? num_strong_var : length(reasonable_disc_vars)
        reasonable_disc_vars = reasonable_disc_vars[1:num_reasonable]

        status, max_gain_var, left_node, right_node, gains, strong_restarts = branch_strong_on!(m,opts,step_obj,
            reasonable_disc_vars, disc2var_idx, opts.strong_restart, counter)

        step_obj.obj_gain += gains

        step_obj.upd_gains = :GainsToTree
        new_gains = copy(step_obj.obj_gain)

        if status == :Infeasible
            return :Infeasible, 0, strong_restarts
        end
    else
        step_obj.upd_gains = :GuessAndUpdate
        new_gains = copy(gains)
    end
    idx = branch_pseudo(m, node, disc2var_idx, new_gains, mu, atol)
    return :Normal, idx, strong_restarts
end

function branch_pseudo(m, node, disc2var_idx, obj_gain, mu, atol)
    # use the one with highest obj_gain which is currently continuous
    idx = 0
    scores, sort_idx = sorted_score_idx(node.solution, obj_gain, disc2var_idx, mu)
    for l_idx in sort_idx
        var_idx = disc2var_idx[l_idx]
        if !is_type_correct(node.solution[var_idx], m.var_type[var_idx], atol)
            u_b = node.u_var[var_idx]
            l_b = node.l_var[var_idx]
            # if the upper bound is the lower bound => should be type correct
            @assert !isapprox(u_b, l_b; atol=atol)
            idx = var_idx
            break
        end
    end
    return idx
end

function sorted_score_idx(x, gains, i2v, mu)
    g_minus, g_minus_c = gains.minus, gains.minus_counter
    g_plus, g_plus_c = gains.plus, gains.plus_counter
    g_minus_c += map(i -> (i == 0) && (i = 1), g_minus_c)
    g_plus_c += map(i -> (i == 0) && (i = 1), g_plus_c)
    scores = [score(f_minus(x[i2v[i]])*g_minus[i]/g_minus_c[i],f_plus(x[i2v[i]])*g_plus[i]/g_plus_c[i],mu) for i=1:length(g_minus)]
    sortedidx = sortperm(scores; rev=true)
    infsortedidx = sortperm(gains.inf_counter[sortedidx]; rev=true)
    return scores,sortedidx[infsortedidx]
end

"""
    Score function from
    http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.92.7117&rep=rep1&type=pdf
"""
function score(q_m, q_p, mu)
    minq = q_m < q_p ? q_m : q_p
    maxq = q_m < q_p ? q_p : q_m
    return (1-mu)*minq+mu*maxq
end

function diff_obj(opts, node, cnode)
    if state_is_optimal(cnode.relaxation_state; allow_almost=opts.allow_almost_solved)
        return abs(node.best_bound - cnode.best_bound)
    else
        return Inf
    end
end

f_plus(x) = ceil(x)-x
f_minus(x) = x-floor(x)
sigma_plus(opts, node,r_nd,x) = diff_obj(opts, node,r_nd)/f_plus(x)
sigma_minus(opts, node,l_nd,x) = diff_obj(opts, node,l_nd)/f_minus(x)


text:
using Colors
using LazySets
using Luxor

# Auxiliary drawing functions
function make_drawing(width, height, img_path, bkg_color, origin_p)
    d = Drawing(width, height, img_path)
    background(bkg_color)
    origin(origin_p)
    return d
end

## Canvas constants
width = 500
height = 500
path = "jarvis_march.png"
color = RGBA(0, 0, 0)
op = [width / 2, height / 2]

# Setting up drawing
my_draw = make_drawing(width, height, path, color, Point(op...))

sethue(RGBA(1.0, 0.9411764705882353, 0.0))
pts_list = [[rand(-250:250), rand(-250:250)] for _ in 1:10]
hull = convex_hull(pts_list)

map2luxor.(pts_list) .|> cp -> circle(cp, 15; action = :fill)
hull = map2luxor.(hull)

sethue(RGBA(1.0, 1.0, 1.0))
poly(hull, action = :stroke; close = true)

finish()

ext: jl

max_stars_repo_name: TheCedarPrince/DataExplorers

max_issues_repo_name: TheCedarPrince/DataExplorers

max_forks_repo_name: TheCedarPrince/DataExplorers


text:
"""
# DESCRIPTION
    Database for managing airfoil polar information.

# AUTHORSHIP
  * Author    : Eduardo J. Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Jun 2020
  * License   : MIT
"""
module AirfoilDatabase

# ------------ GENERIC MODULES -------------------------------------------------
import DataStructures: OrderedDict
import DataFrames
import CSV
import JuliaDB
import Plots
import Interact

const jdb = JuliaDB
const plt = Plots
const ntrct = Interact

# ------------ FLOW CODES ------------------------------------------------------
# https://github.com/byuflowlab/AirfoilPrep.jl
import AirfoilPrep
const ap = AirfoilPrep

# ------------ GLOBAL VARIABLES ------------------------------------------------
const module_path = splitdir(@__FILE__)[1]                   # Path to this module
const def_database = joinpath(module_path, "../database/")   # Path to default database

# ------------ DATA STRUCTURE --------------------------------------------------
const DIR_XY = "xy"         # xy directory
const DIR_CL = "Cl"         # Cl directory
const DIR_CD = "Cd"         # Cd directory
const DIR_CM = "Cm"         # Cm directory
const DIR_XUP = "xupsep"    # xupsep directory
const DIR_XLO = "xlosep"    # xlosep directory
const DEF_NDXFL = "index.csv" # Default indexing file name
const RQRD = "required"     # Default value for requiered fields
const PTNL_S = " "          # Default value for optional String fields
const PTNL_R = 0            # Default value for optional Real fields
const PTNL_I = 0            # Default value for optional Int fields
const PTNL_B = true         # Default value for optional Boolean fields
                            # Labels of each database entry
const LBLS = OrderedDict(#  hash_key        => (csv_header, default_value, type, description)
                            :airfoilname    => ("Airfoil", RQRD, String, "Airfoil name"),
                            :re             => ("Re", PTNL_R, Real, "Reynolds number"),
                            :ma             => ("Mach", PTNL_R, Real, "Mach number"),
                            :npanels        => ("Panels", PTNL_I, Int, "Number of panels (XFOIL)"),
                            :ncrit          => ("Ncrit", PTNL_R, Real, "XFOIL Ncrit parameter"),
                            :deg            => ("Degrees", PTNL_B, Bool, "XFOIL Ncrit parameter"),
                            :xyfile         => ("xy file", RQRD, String, "Airfoil contour file"),
                            :clfile         => ("Cl file", PTNL_S, String, "Cl file"),
                            :cdfile         => ("Cd file", PTNL_S, String, "Cd file"),
                            :cmfile         => ("Cm file", PTNL_S, String, "Cm file"),
                            :xupsepfile     => ("xupsep file", PTNL_S, String, "xupsep file"),
                            :xlosepfile     => ("xlosep file", PTNL_S, String, "xlosep file"),
                            :diff           => ("Differentiator", PTNL_I, Int, "Differentiator number (for accepting duplicates)"),
                        )

const RQRD_FIELDS = [key for (key, val) in LBLS  if val[2]==RQRD]   # Required fields
const ID_FIELDS = [:airfoilname, :re, :ma, :npanels, :ncrit, :deg, :diff] # Fields that make an entry unique

const HEADER2FIELD = Dict((val[1], key) for (key, val) in LBLS)     # Hash CSV header to field
const FIELD2HEADER = Dict((key, val[1]) for (key, val) in LBLS)     # Hash field to CSV header
const HEADER_2FIELD = Dict((replace(val[1], " "=>"_"), key) for (key, val) in LBLS)  # Hash CSV header (with spaces as _) to field
const FIELD2_HEADER = Dict((key, replace(val[1], " "=>"_")) for (key, val) in LBLS)  # Hash field to CSV header (with spaces as _)

const DIRS = Dict(  :xyfile=>DIR_XY,                                # Hash field to directory
                    :clfile=>DIR_CL, :cdfile=>DIR_CD, :cmfile=>DIR_CM,
                    :xupsepfile=>DIR_XUP, :xlosepfile=>DIR_XLO)

"Convert LBLS into a DataFrame"
function _lbls2dataframe()
    df = DataFrames.DataFrame(;(
                                (Symbol("Field"), collect(keys(LBLS))),
                                (Symbol("Header"), [val[1] for val in values(LBLS)]),
                                (Symbol("Default value"), [val[2] for val in values(LBLS)]),
                                (Symbol("Value type"), [val[3] for val in values(LBLS)]),
                                (Symbol("Description"), [val[4] for val in values(LBLS)]),
                             )...)
    return df
end

"Convert LBLS into a pretty String"
function _lbls2prettystring()
    b = IOBuffer();
    t = TextDisplay(b);
    display(t, LBLS);
    s = String(take!(b))
    return s
end

# ------------ HEADERS ---------------------------------------------------------
# Load headers
for header_name in ["manipulation", "visualization"]
    include("AirfoilDatabase_"*header_name*".jl")
end

end # END OF MODULE


text:
module TFRecords
using ProtoBuf
using CRC32c: crc32c

export TFRecord, readexample, writeexample, iterate, parse_example, build_example

# run(ProtoBuf.protoc(`-I=~/src/tensorflow/ --julia_out=src/generated ~/src/tensorflow/tensorflow/core/example/example.proto`))
include("generated/tensorflow.jl")

Example = tensorflow.Example
Feature = tensorflow.Feature
Features = tensorflow.Features
BytesList = tensorflow.BytesList
FloatList = tensorflow.FloatList
Int64List = tensorflow.Int64List

struct TFRecord
    filename::String
end

# https://www.tensorflow.org/tutorials/load_data/tfrecord#tfrecords_format_details
masked_crc(crc::UInt32) = ((crc >> 15) | (crc<<17)) + 0xa282ead8

lengthcrc(l::UInt64) = crc32c([x for x in reinterpret(UInt8, [l])])

"""
    readexample(f)

Read one Example proto from an IOStream, verifying checksum along the way.
"""
function readexample(s::IO)
    length = read(s, UInt64)
    length_crc = read(s, UInt32)
    computed_length_crc = masked_crc(lengthcrc(length))
    @assert(computed_length_crc == length_crc)

    data = read(s, length)
    computed_data_crc = masked_crc(crc32c(data))

    data_crc = read(s, UInt32)
    @assert(computed_data_crc == data_crc)

    iob = PipeBuffer(data)
    readproto(iob, Example())
end

"""
    writeexample(s, e)

Write one Example proto to an IOStream in the tfrecord format.
"""
function writeexample(s::IO, e::Example)
    iob = PipeBuffer()
    writeproto(iob, e)
    length = UInt64(iob.size)
    w = write(s, length)
    w += write(s, masked_crc(lengthcrc(length)))
    w += write(s, iob.data)
    w += write(s, masked_crc(crc32c(iob.data)))
end

"""
    iterate(r)

Iterate over a tfrecord file and produce one example at a time.

# Example
```julia-repl
julia> r = TFRecord("data.tfrecord")
TFRecord("data.tfrecord")
julia> records = Base.Iterators.Stateful(r)
Base.Iterators.Stateful{ ...
julia> record = popfirst!(records)
TFRecords.tensorflow.Example(TFRecords.tensorflow.Features(Dict{AbstractString,TFRecords.tensorflow.Feature} ...
```
"""
function Base.iterate(R::TFRecord, state::Union{Nothing,IO}=nothing)
    if state == nothing
        state = open(R.filename)
    end
    try
        (readexample(state), state)
    catch err
        if isa(err, EOFError)
            nothing
        else
            rethrow()
        end
    end
end

"""
    woo(f)

Extract the defined alternative value from the protobuf.
"""
function woo(f::Feature)
    woo = which_oneof(f, :kind)
    if woo == :bytes_list
        f.bytes_list.value
    elseif woo == :int64_list
        f.int64_list.value
    elseif woo == :float_list
        f.float_list.value
    end
end

"""
    parse_example(e)

Convert the ProtoBuf Example into a Dict{String,Array{T, 1}} with an entry for
each defined feature
"""
function parse_example(e::Example)
    feature_keys = e.features.feature.keys
    defined_keys = [feature_keys[i] for i in 1:length(feature_keys) if isassigned(feature_keys, i)]
    Dict(k=>woo(e.features.feature[k]) for k in defined_keys)
end

make_feature(f::Array{Array{UInt8, 1},1}) = Feature(bytes_list=BytesList(value=f))
make_feature(f::Array{Int64, 1}) = Feature(int64_list=Int64List(value=f))
make_feature(f::Array{Float32, 1}) = Feature(float_list=FloatList(value=f))
PossibleFeatureType = Union{Int64, Float32, Array{UInt8, 1}}

"""
    build_example(d)

Convert a dict like the one returned by parse_example into an Example for
serialization to a file
"""
build_example(d::Dict{T, Array{U, 1}} where {T<:AbstractString, U<:PossibleFeatureType}) = Example(features=Features(feature=Dict{AbstractString, Feature}(k => make_feature(v) for (k, v) in d)))

end


text:
using Test
using Quadmath

@testset "conversion $T" for T in (Float64, Int32, Int64, BigFloat, BigInt)
    @test Float128(T(1)) + Float128(T(2)) == Float128(T(3))
    @test Float128(T(1)) + Float128(T(2)) <= Float128(T(3))
    @test Float128(T(1)) + Float128(T(2)) != Float128(T(4))
    @test Float128(T(1)) + Float128(T(2)) < Float128(T(4))
    if isbitstype(T)
        @test T(Float128(T(1)) + Float128(T(2))) === T(3)
    else
        @test T(Float128(T(1)) + Float128(T(2))) == T(3)
    end
end

@test Base.exponent_one(Float128) == reinterpret(UInt128, Float128(1.0))

@testset "BigFloat" begin
    x = parse(Float128, "0.1")
    y = parse(Float128, "0.2")
    @test Float64(x+y) == Float64(BigFloat(x) + BigFloat(y))
    @test x+y == Float128(BigFloat(x) + BigFloat(y))
end

@testset "BigInt" begin
    x = parse(Float128, "100.0")
    y = parse(Float128, "25.0")
    @test Float64(x+y) == Float64(BigInt(x) + BigInt(y))
    @test x+y == Float128(BigInt(x) + BigInt(y))
end

@testset "flipsign" begin
    x = Float128( 2.0)
    y = Float128(-2.0)
    @test x == flipsign(y, -one(Float128))
    @test y == flipsign(y,  1)
end

@testset "modf" begin
    x = Float128(pi)
    fpart, ipart = modf(x)
    @test x == ipart + fpart
    @test signbit(fpart) == signbit(ipart) == false

    y = Float128(-pi)
    fpart, ipart = modf(y)
    @test y == ipart + fpart
    @test signbit(fpart) == signbit(ipart) == true

    z = x^3
    fpart, ipart = modf(x) .+ modf(y)
    @test x+y == ipart+fpart
end


text:
function find_roots(f, a::Real, b::Real, args...;
                    no_pts::Int=101,
                    abstol::Real=10*eps(), reltol::Real=10*eps(), ## should be abstol, reltol as used.
                    kwargs...)

    @assert a != b
    @assert no_pts > 1

    ( a < b ) || ( (a, b) = (b, a) )

    work_f(x) = real(float(f(complex(x))))

    root_list = _find_recursive_roots(
      work_f, a, b, args...;
      no_pts=no_pts, abstol=abstol, reltol=reltol, kwargs...
    )

    if length(root_list) > Int(ceil( (1/4) * no_pts ))
      custom_log("Redoing root hunt on larger grid")
      return find_roots(
        f, a, b, args...;
        no_pts = 10*no_pts, abstol=abstol, reltol=reltol, kwargs...
      )
    end

    sort!(root_list)

    root_list
end

function _find_recursive_roots(f, a::Real, b::Real, args...;
                      cur_depth::Int=1,
                      no_pts::Int=101,
                      abstol::Real=10*eps(), reltol::Real=10*eps(), ## should be abstol, reltol as used.
                      kwargs...)

    cur_range = collect(linspace(a, b, no_pts))

    root_list = find_root_list(f, cur_range, abstol, reltol; kwargs...)

    isempty(root_list) && return root_list

    sub_no_pts = no_pts
    sub_no_pts /= 1 + length(root_list) / 3
    sub_no_pts /= 2 ^ ( cur_depth - 1 )
    sub_no_pts = Int(round(sub_no_pts))

    ( sub_no_pts > 1 ) || return root_list

    cur_intervals = zip(
        vcat(first(cur_range),root_list),
        vcat(root_list,last(cur_range))
    )

    for (cur_index, (cur_a, cur_b)) in enumerate(cur_intervals)
        loose_isapprox(cur_a, cur_b, abstol, reltol) && continue

        tmp_a = find_next_non_root(cur_a, f, cur_b, reltol, abstol)
        tmp_b = find_next_non_root(cur_b, f, cur_a, reltol, abstol)

        isnan(tmp_a) || isnan(tmp_b) && continue
        loose_isapprox(tmp_a, tmp_b, abstol, reltol) && continue
        ( tmp_a < tmp_b ) || continue

        if cur_depth < 5 && 1 < cur_index < length(cur_intervals)
          tmp_no_pts = no_pts
        else
          tmp_no_pts = sub_no_pts
        end

        cur_roots = _find_recursive_roots(
            f, tmp_a, tmp_b, args...;
            no_pts=tmp_no_pts, abstol=abstol, reltol=reltol,
            cur_depth=(cur_depth+1), kwargs...
        )

        append!(root_list, cur_roots)
    end

    root_list
end

function find_next_non_root(cur_root, cur_func, barrier_point, reltol, abstol)
  cur_direction = sign(barrier_point - cur_root)

  cur_adder = abstol * cur_direction
  cur_value = cur_root + cur_adder

  while sign(barrier_point - cur_value) == cur_direction
    if !loose_isapprox(cur_value, cur_root, abstol, reltol)
      cur_f = cur_func(cur_value)
      isapprox(cur_f, 0.0, atol=abstol) || break
    end

    cur_adder *= 2
    cur_value += cur_adder
  end

  is_in_range = ( sign(barrier_point - cur_value) == cur_direction )
  is_in_range || return NaN

  for tmp_value in linspace(cur_value - cur_adder, cur_value)
    ( sign(tmp_value - cur_root) == cur_direction ) || continue

    loose_isapprox(tmp_value, cur_root, abstol, reltol) && continue

    cur_f = cur_func(tmp_value)

    isapprox(cur_f, 0.0, atol=abstol) && continue

    cur_value = tmp_value
    break
  end

  ( cur_root == cur_value ) && return NaN

  cur_value
end

function find_root_list(f, cur_range::AbstractVector{T}, abstol::Real, reltol::Real; kwargs...) where T <: Real
  isempty(cur_range) && return Real[]

  cur_roots = find_bisection_roots(f, cur_range, abstol, reltol; kwargs...)
  isempty(cur_roots) || return cur_roots

  cur_root = find_order_root(f, cur_range, abstol, reltol; kwargs...)
  isnan(cur_root) && return Real[]

  cur_roots = [cur_root]
  cur_roots
end

function find_bisection_roots(f, cur_range::AbstractVector{T}, abstol::Real, reltol::Real; kwargs...) where T <: Real
  no_pts = length(cur_range)

  root_list = Real[]

  cur_values = f.(cur_range)

  filtered_cur_values = filter(isfinite, cur_values)

  isempty(filtered_cur_values) && return []

  cur_signs = map(sign, cur_values[1:end-1] .* cur_values[2:end])

  for (cur_index, cur_sign) in enumerate(cur_signs)
    isinf(cur_sign) && continue
    isnan(cur_sign) && continue

    cur_a, cur_b = cur_range[cur_index:cur_index+1]

    if iszero(cur_sign)
      iszero(cur_values[cur_index]) && push!(root_list, cur_a)

      ( cur_index == no_pts - 1 ) &&
        iszero(cur_values[cur_index+1]) && push!(root_list, cur_b)

      continue
    end

    ( cur_sign > 0 ) && continue

    cur_root = NaN
    try
      cur_root = find_zero(f, [cur_a, cur_b], FalsePosition(), maxevals=16)
      @assert isapprox(f(cur_root), 0.0, atol=abstol)
    catch
      cur_root = custom_bisection(f, cur_a, cur_b, f(cur_a), f(cur_b); abstol=abstol, kwargs...)
    end
    isnan(cur_root) || push!(root_list, cur_root)
  end

  isempty(root_list) || return root_list

  super_diff = diff(filtered_cur_values)
  cur_attempt_index = []
  if length(super_diff) == 0 || ( any(is_positive, super_diff) && any(!is_positive, super_diff) )
    cur_attempt_index = [
      findfirst(isfinite, cur_values),
      findlast(isfinite, cur_values)
    ]
  else
    if all(is_positive, filtered_cur_values) || all(!is_positive, filtered_cur_values)
      cur_is_increasing = is_positive(first(super_diff))

      first_finite_index = findfirst(isfinite, cur_values)
      last_finite_index = findlast(isfinite, cur_values)

      if cur_is_increasing
        if cur_values[first_finite_index] < 0
          cur_attempt_index = [last_finite_index]
        else
          cur_attempt_index = [first_finite_index]
        end
      else
        if cur_values[first_finite_index] < 0
          cur_attempt_index = [first_finite_index]
        else
          cur_attempt_index = [last_finite_index]
        end
      end
    else
      tmp_signs = map(sign, filtered_cur_values[1:end-1] .* filtered_cur_values[2:end])
      cur_attempt_index = findin(cur_values, filtered_cur_values[find(!is_positive, tmp_signs)])
    end
  end

  cur_min_index = indmin(map(cur_value -> isfinite(cur_value) ? abs(cur_value) : Inf, cur_values))

  push!(cur_attempt_index, cur_min_index)

  cur_attempt_index = unique(cur_attempt_index)

  for cur_index in cur_attempt_index
    cur_point = cur_range[cur_index]

    cur_diff_list = diff(cur_range)
    filter_approx!(cur_diff_list; atol=10*eps())

    @assert length(cur_diff_list) == 1
    cur_diff = cur_diff_list[1]

    if 1 < cur_index < length(cur_values)
      cur_c = cur_point
      cur_f_c = f(cur_c)

      cur_a = cur_point - cur_diff
      cur_b = cur_point + cur_diff

      cur_root = custom_bisection(f, cur_a, cur_c, f(cur_a), cur_f_c; abstol=abstol, kwargs...)
      isnan(cur_root) || return [cur_root]

      cur_root = custom_bisection(f, cur_c, cur_b, cur_f_c, f(cur_b); abstol=abstol, kwargs...)
      isnan(cur_root) || return [cur_root]
    end
  end

  []
end

function find_order_root(f, cur_range::AbstractVector{T}, abstol::Real, reltol::Real; kwargs...) where T <: Real
  min_value = first(cur_range)
  max_value = last(cur_range)

  cur_guesses = view( cur_range , (3:4:length(cur_range)-2) )

  cur_root = NaN

  for cur_guess in cur_guesses
    tmp_root = NaN

    try
      tmp_root = find_zero(
        f, cur_guess, Order8();
        maxevals=10, abstol=abstol, reltol=reltol
      )
    catch cur_error
        continue
    end

    ( min_value < tmp_root < max_value ) || continue

    isapprox(f(tmp_root), 0.0, atol=abstol) || continue

    cur_root = tmp_root

    break
  end

  cur_root
end

function loose_isapprox(first_value, second_value, atol, rtol)
  isapprox(first_value, second_value, atol=atol) && return true
  isapprox(first_value, second_value, rtol=rtol) && return true

  false
end

function custom_bisection(f, cur_a::Number, cur_b::Number, cur_f_a::Number, cur_f_b::Number, cur_flag::Bool=true; abstol::Number=10*eps(), cur_c::Number=NaN, cur_f_c::Number=NaN, cur_d::Number=NaN, cur_f_d::Number=NaN, bad_streak::Number=0, kwargs...)
  ( bad_streak > 6 ) && return NaN

  init_a = cur_a
  init_b = cur_b

  init_f_a = cur_f_a
  init_f_b = cur_f_b

  isnan(cur_c) && ( cur_c = middle(cur_a, cur_b) )
  isnan(cur_f_c) && ( cur_f_c = f(cur_c) )

  if abs(cur_f_a) < abs(cur_f_b)
    cur_a, cur_b = cur_b, cur_a
    cur_f_a, cur_f_b = cur_f_b, cur_f_a
  end

  isapprox(cur_f_a, 0.0, atol=abstol) && return cur_a
  isapprox(cur_f_b, 0.0, atol=abstol) && return cur_b
  isapprox(cur_f_c, 0.0, atol=abstol) && return cur_c
  isapprox(cur_f_d, 0.0, atol=abstol) && return cur_d

  isapprox(cur_f_a, cur_f_b, atol=abstol) && return NaN
  isapprox(cur_a, cur_b, atol=2*eps()) && return NaN

  is_bad_a = isinf(cur_f_a) || isnan(cur_f_a)
  is_bad_b = isinf(cur_f_b) || isnan(cur_f_b)
  is_bad_c = isinf(cur_f_c) || isnan(cur_f_c)
  is_bad_d = isinf(cur_f_d) || isnan(cur_f_d)

  ( !is_bad_a && !is_bad_b && cur_f_a * cur_f_b > 0 ) && return NaN

  ( is_bad_a && is_bad_b ) && return NaN

  if is_bad_a || is_bad_b || is_bad_c
    if !is_bad_d && ( isapprox(cur_c, cur_a, atol=abstol) || isapprox(cur_c, cur_b, atol=abstol) )
      cur_g = cur_d
      cur_f_g = cur_f_d
    else
      cur_g = cur_c
      cur_f_g = cur_f_c
    end

    cur_f_list = [cur_f_a, cur_f_b, cur_f_g]
    filter!(isfinite, cur_f_list)

    is_on_good_streak = isfinite(cur_f_g)
    is_on_good_streak &= abs(cur_f_g) == minimum(map(abs, cur_f_list))

    is_on_good_streak &= !isapprox(cur_f_g, cur_f_a, atol=abstol)
    is_on_good_streak &= !isapprox(cur_f_g, cur_f_b, atol=abstol)

    is_on_good_streak &= !isapprox(cur_g, cur_a, atol=abstol)
    is_on_good_streak &= !isapprox(cur_g, cur_b, atol=abstol)

    if is_on_good_streak
      bad_streak -= 1

      filter!(work_f -> work_f == cur_f_g, cur_f_list)

      second_min_f = minimum(map(abs, cur_f_list))

      ( second_min_f / abs(cur_f_g) > 8 ) &&
        ( bad_streak -= 1 )

      bad_streak = max(0, bad_streak)
    else
      if length(cur_f_list) == 1
        bad_streak += 1
      else
        bad_streak += 2
      end
    end

    cur_root = NaN

    !isapprox(cur_f_a, cur_f_g, atol=abstol) &&
      ( cur_root = custom_bisection(f, cur_a, cur_g, cur_f_a, cur_f_g; abstol=abstol, bad_streak=bad_streak, kwargs...) )

    ( isnan(cur_root) && !isapprox(cur_f_g, cur_f_b, atol=abstol) ) &&
      ( cur_root = custom_bisection(f, cur_g, cur_b, cur_f_g, cur_f_b; abstol=abstol, bad_streak=bad_streak, kwargs...) )

    return cur_root
  end

  cur_d = NaN
  cur_s = 0.0

  if !isapprox(cur_f_a, cur_f_c, atol=abstol) && !isapprox(cur_f_b, cur_f_c, atol=abstol)
    cur_term_1 = cur_f_b / ( cur_f_a - cur_f_c )
    cur_term_2 = cur_f_c / ( cur_f_a - cur_f_b )
    cur_term_3 = cur_f_a / ( cur_f_b - cur_f_c )

    cur_s += cur_a * cur_term_1 * cur_term_2
    cur_s -= cur_b * cur_term_2 * cur_term_3
    cur_s += cur_c * cur_term_3 * cur_term_1
  else
    cur_term = ( cur_b - cur_a )
    cur_term /= ( cur_f_b - cur_f_a )

    cur_s = cur_b
    cur_s -= cur_f_b * cur_term
  end

  cur_bool = ( cur_s < cur_b && cur_s < ( ( 3 * cur_a + cur_b ) / 4 ) )
  cur_bool |= ( cur_s > cur_b && cur_s > ( ( 3 * cur_a + cur_b ) / 4 ) )

  if cur_flag
    cur_diff = abs( cur_b - cur_c )
  else
    cur_diff = abs( cur_d - cur_c )
  end

  cur_bool |= abs( cur_s - cur_b ) >= ( cur_diff / 2 )
  cur_bool |= isapprox(cur_diff, 0.0, atol=abstol)

  cur_flag = cur_bool
  cur_bool && ( cur_s = middle(cur_a, cur_b) )

  cur_d = cur_c
  cur_c = cur_b

  cur_f_d = cur_f_c
  cur_f_c = cur_f_b

  cur_f_s = f(cur_s)

  if cur_f_a * cur_f_s < 0
    cur_b = cur_s
    cur_f_b = cur_f_s
  else
    cur_a = cur_s
    cur_f_a = cur_f_s
  end

  is_bad_f_a = isapprox(cur_f_a, init_f_a, atol=abstol)
  is_bad_f_b = isapprox(cur_f_b, init_f_b, atol=abstol)

  is_bad_f_a |= abs(cur_f_a) > abs(init_f_a)
  is_bad_f_b |= abs(cur_f_b) > abs(init_f_b)

  if ( is_bad_f_a && is_bad_f_b )
    bad_streak += 1
  else
    bad_streak -= 2
    bad_streak = max(0, bad_streak)
  end

  is_unchanged = isapprox(cur_a, init_a, atol=2*eps())
  is_unchanged &= isapprox(cur_b, init_b, atol=2*eps())

  if is_unchanged
    cur_c = middle(cur_a, cur_b)
    cur_f_c = f(cur_c)

    cur_root = custom_bisection(f, cur_a, cur_c, cur_f_a, cur_f_c; abstol=abstol, bad_streak=bad_streak, kwargs...)
    isnan(cur_root) || return cur_root

    cur_root = custom_bisection(f, cur_c, cur_b, cur_f_c, cur_f_b; abstol=abstol, bad_streak=bad_streak, kwargs...)
    return cur_root
  end

  custom_bisection(
    f, cur_a, cur_b, cur_f_a, cur_f_b, cur_flag;
    abstol=abstol, bad_streak=bad_streak,
    cur_c=cur_c, cur_f_c=cur_f_c,
    cur_d=cur_d, cur_f_d=cur_f_d,
    kwargs...
  )
end


text:
using DataFrames, PyPlot, Statistics, FeatherFiles

AREA = 16.0 # area of the solar panels

df = DataFrame(load("./data/solaryear.feather"))

plot(df.Power_Generation*AREA)
title("Solar power, 16 m², Gouda, NL")

mean(df.Power_Generation*AREA)


text:
function tracem(x::Array{T,2}) where T <: AbstractFloat

  # We require the number of rows to be greater than the number of columns, so that m is greater than one.

  trans = false

  if size(x,1) < size(x,2)
    x = x'
    trans = true
  end

  n = size(x,2)
  m = int(size(x,1)/n)

  y = zeros(m,1)  # We want this to be a 2-d array for subsequent matrix multiplication

  for i = 1:m

    y[i,1] = trace(x[(n*(i-1)+1):i*n,1:n])

  end

  if trans == true
    y = y'
  end

  return y

end

function permutation(n1::S,n2::S) where S <: Int

  nn = n1*n2
  p = zeros(S,nn,nn)

  i = 1
  j = 1
  for r = 1:nn

    if j > n2
      j = 1
      i += 1
    end
    c = i+n1*(j-1)
    p[r,c] = 1

    j += 1

  end

  return p

end

function convert_second_order(model::M) where M <: Lombardo_Sutherland_Soln

  ssh       = copy(model.ssh)
  hx        = copy(model.hx)
  hxx       = copy(model.hxx)
  ssg       = copy(model.ssg)
  gx        = copy(model.gx)
  gxx       = copy(model.gxx)
  sigma     = copy(model.sigma)
  grc       = copy(model.grc)
  soln_type = copy(model.soln_type)

  (n1,n2) = size(hx)
  (n3,n4) = size(gx)

  new_hxx = reshape(vech_to_vec(n2)*hxx[1,:]',n2,n2)
  for i = 2:n1

    new_hxx = [new_hxx; reshape(vech_to_vec(n2)*hxx[i,:]',n2,n2)]

  end

  (n5,n6) = size(new_hxx)

  j = 0
  for i = 1:n5

    j += 1
    new_hxx[i,j] = 2*new_hxx[i,j]
    if j == n6
      j = 0
    end

  end

  new_hxx = new_hxx/2

  new_gxx = reshape(vech_to_vec(n4)*gxx[1,:]',n4,n4)
  for i = 2:n3

    new_gxx = [new_gxx; reshape(vech_to_vec(n4)*gxx[i,:]',n4,n4)]

  end

  (n7,n8) = size(new_gxx)

  j = 0
  for i = 1:n7

    j += 1
    new_gxx[i,j] = 2*new_gxx[i,j]
    if j == n8
      j = 0
    end

  end

  new_gxx = new_gxx/2

  soln = Gomme_Klein_Soln(ssh,hx,new_hxx,ssg,gx,new_gxx,sigma,grc,soln_type)

  return soln

end


text:
module ITensorsBenchmarks

ENV["GKSwstype"]= "100"

using Dates
using DelimitedFiles
using ITensors
using LinearAlgebra
using Plots
using Random
using Suppressor

export runbenchmarks,
       plotbenchmarks

#
# OPTIONS
#

const DEFAULT_OPTIONS = Dict{String, Any}()

# A range of maximum dimensions for each benchmark
DEFAULT_OPTIONS["maxdims"] = Dict{String, StepRange{Int64, Int64}}()
const maxdims = DEFAULT_OPTIONS["maxdims"]

# Descriptions of each benchmark
DEFAULT_OPTIONS["descriptions"] = Dict{String, String}()
const descriptions = DEFAULT_OPTIONS["descriptions"]

# Get a list of the built-in benchmarks
_benchmarks = readdir(joinpath(pkgdir(@__MODULE__), "src"))
_benchmarks = filter!(file -> startswith(file, "bench_"), _benchmarks)
_benchmarks .= replace.(_benchmarks, "bench_" => "")

for benchmark in _benchmarks
  include("bench_$benchmark/runbenchmark.jl")
end

DEFAULT_OPTIONS["write_results"] = false
DEFAULT_OPTIONS["cpp_or_julia"] = nothing
DEFAULT_OPTIONS["blocksparse_num_threads"] = 1
DEFAULT_OPTIONS["blas_num_threads"] = 1
DEFAULT_OPTIONS["test"] = false
DEFAULT_OPTIONS["benchmarks"] = _benchmarks
DEFAULT_OPTIONS["julia_itensor_version"] = v"0.1.37" #ITensors.version()
DEFAULT_OPTIONS["cpp_itensor_version"] = v"3.1.6"
DEFAULT_OPTIONS["cpp_itensor_dir"] = joinpath(pkgdir(@__MODULE__), "deps", "itensor_v$(DEFAULT_OPTIONS["cpp_itensor_version"])")
DEFAULT_OPTIONS["splitblocks"] = false

const OPTIONS = deepcopy(DEFAULT_OPTIONS)
function reset_options!()
  empty!(OPTIONS)
  OPTIONS_COPY = deepcopy(DEFAULT_OPTIONS)
  for k in keys(OPTIONS)
    OPTIONS[k] = OPTIONS_COPY[k]
  end
end

#
# Files needed from ITensors.jl examples
#

examples_dir = joinpath(pkgdir(ITensors), "examples", "src")
include(joinpath(examples_dir, "2d_classical_ising.jl"))
include(joinpath(examples_dir, "trg.jl"))
include(joinpath(examples_dir, "ctmrg_isotropic.jl"))
include(joinpath(examples_dir, "electronk.jl"))
include(joinpath(examples_dir, "hubbard.jl"))

include("runbenchmarks.jl")
include("plotbenchmarks.jl")

end # module


text:
# This file defines the PrerequisiteGroup type, used to define a group of
#   prerequisites that all need to be satisfied.

# The PrerequisiteGroup type requires the Prerequisite type.
requiredTypes = [ "prerequisite" ]

for reqType in requiredTypes
    if !isdefined( Symbol( uppercase( string( reqType[ 1 ] ) ) * reqType[ 2:end ] ) )
        include( joinpath( typePath, reqType * ".jl" ) )
    end  # if !isdefined( Symbol( ...
end  # for reqType in requiredTypes

export PrerequisiteGroup
type PrerequisiteGroup
    # The list of prerequisites.
    prereqs::Vector{Prerequisite}

    # Constructor.
    function PrerequisiteGroup()
        prereqGroup = new()
        prereqGroup.prereqs = Vector{Prerequisite}()
        return prereqGroup
    end  # PrerequisiteGroup()
end  # type PrerequisiteGroup


text:
function convert_type(obj::OrbitType, arg0::FieldOrbit)
    return jcall(obj, "convertType", FieldOrbit, (FieldOrbit,), arg0)
end

function convert_type(obj::OrbitType, arg0::Orbit)
    return jcall(obj, "convertType", Orbit, (Orbit,), arg0)
end

function get_drivers(obj::OrbitType, arg0::jdouble, arg1::Orbit, arg2::PositionAngle)
    return jcall(obj, "getDrivers", ParameterDriversList, (jdouble, Orbit, PositionAngle), arg0, arg1, arg2)
end

function map_array_to_orbit(obj::OrbitType, arg0::Vector{RealFieldElement}, arg1::Vector{RealFieldElement}, arg2::PositionAngle, arg3::FieldAbsoluteDate, arg4::RealFieldElement, arg5::Frame)
    return jcall(obj, "mapArrayToOrbit", FieldOrbit, (Vector{RealFieldElement}, Vector{RealFieldElement}, PositionAngle, FieldAbsoluteDate, RealFieldElement, Frame), arg0, arg1, arg2, arg3, arg4, arg5)
end

function map_array_to_orbit(obj::OrbitType, arg0::Vector{jdouble}, arg1::Vector{jdouble}, arg2::PositionAngle, arg3::AbsoluteDate, arg4::jdouble, arg5::Frame)
    return jcall(obj, "mapArrayToOrbit", Orbit, (Vector{jdouble}, Vector{jdouble}, PositionAngle, AbsoluteDate, jdouble, Frame), arg0, arg1, arg2, arg3, arg4, arg5)
end

function map_orbit_to_array(obj::OrbitType, arg0::FieldOrbit, arg1::PositionAngle, arg2::Vector{RealFieldElement}, arg3::Vector{RealFieldElement})
    return jcall(obj, "mapOrbitToArray", void, (FieldOrbit, PositionAngle, Vector{RealFieldElement}, Vector{RealFieldElement}), arg0, arg1, arg2, arg3)
end

function map_orbit_to_array(obj::OrbitType, arg0::Orbit, arg1::PositionAngle, arg2::Vector{jdouble}, arg3::Vector{jdouble})
    return jcall(obj, "mapOrbitToArray", void, (Orbit, PositionAngle, Vector{jdouble}, Vector{jdouble}), arg0, arg1, arg2, arg3)
end

function value_of(::Type{OrbitType}, arg0::JString)
    return jcall(OrbitType, "valueOf", OrbitType, (JString,), arg0)
end

function values(::Type{OrbitType})
    return jcall(OrbitType, "values", Vector{OrbitType}, ())
end


text:
module Traffic

include("Traffic/Simple.jl")

end


text:
# test the April grid camera calibration functions

using Test
using AprilTags
using FileIO
using Optim

# using ImageView

##

@testset "April grid camera calibration test" begin

##

# load the image used for calibration
imgfile = joinpath(dirname(dirname(pathof(AprilTags))), "data/CameraCalibration/taggridphoto.jpg")
img = load(imgfile);

##

# construct the cost function
IMGS = Vector{typeof(img)}()
push!(IMGS, img)

detector = AprilTagDetector()
tags = detector.(IMGS);

##

taglength = 0.0315

# assume obj([f_width, f_height, c_width, c_height])
obj = (fc_wh) -> calcCalibResidualAprilTags!(IMGS, tags,
                                            taglength=taglength,
                                            f_width =fc_wh[1],
                                            f_height=fc_wh[2],
                                            c_width =fc_wh[3],
                                            c_height=fc_wh[4]  )
#

## Run the optimization

# f_width, c_width, c_height = size(img, 1), size(img, 1)/2, size(img, 2)/2
# f_height = f_width

# nearby calibration
f_width  = 3370.4878918701756 + 5
f_height = 3352.8348099534364 + 5
c_width  = 2005.641610450976  + 5
c_height = 1494.8282013012076 + 5
# the physical size of the tag
# taglength = 0.0315

result = optimize(obj, [f_width, f_height, c_width, c_height], BFGS(), Optim.Options(iterations = 10, x_tol=1e-6, show_trace=true))

result.minimizer

## draw the result

SEL = 1

cimg_ = deepcopy(IMGS[SEL])
calcCornerProjectionsAprilTags!( cimg_, tags[SEL],
                                  taglength=0.0315,
                                  f_width=f_width,
                                  f_height=f_height,
                                  c_width=c_width,
                                  c_height=c_height,
                                  dodraw=true )
#

# imshow(cimg_)

##

freeDetector!(detector)

##

end

#


text:
function get_data_don(; n=2048, Δsamples=2^3, grid_size=div(2^13, Δsamples))
    file = matopen(joinpath(datadep"Burgers", "burgers_data_R10.mat"))
    x_data = collect(read(file, "a")[1:n, 1:Δsamples:end])
    y_data = collect(read(file, "u")[1:n, 1:Δsamples:end])
    close(file)

    return x_data, y_data
end

function train_don(; n=300, cuda=true, learning_rate=0.001, epochs=400)
    if cuda && has_cuda()
        @info "Training on GPU"
        device = gpu
    else
        @info "Training on CPU"
        device = cpu
    end

    x, y = get_data_don(n=n)

    xtrain = x[1:280, :]'
    ytrain = y[1:280, :]

    xval = x[end-19:end, :]' |> device
    yval = y[end-19:end, :] |> device

    grid = collect(range(0, 1, length=1024)') |> device

    opt = ADAM(learning_rate)

    m = DeepONet((1024,1024,1024), (1,1024,1024), gelu, gelu) |> device

    loss(X, y, sensor) = Flux.Losses.mse(m(X, sensor), y)
    evalcb() = @show(loss(xval, yval, grid))

    data = [(xtrain, ytrain, grid)] |> device
    Flux.@epochs epochs Flux.train!(loss, params(m), data, opt, cb=evalcb)
    ỹ = m(xval |> device, grid |> device)

    diffvec = vec(abs.(cpu(yval) .- cpu(ỹ)))
    mean_diff = sum(diffvec)/length(diffvec)
    return mean_diff
end


text:
function gate(context::NamedTuple,data::DataFrame,embedding::NamedTuple,selections::NamedTuple;channelRange=range(-3,7,length=50))
    request = context.request
	try
		response = HTTP.Response(200)
        HTTP.Messages.setheader(response,"Content-type" => "text/csv")
        feature = JSON.parse(String(request.body))

        body = IOBuffer(response.body, write=true)
        t = @elapsed save( Stream(format"CSV",body), gate(feature["coordinates"],data,embedding,selections;channelRange=channelRange) )

		@info """$(request.method) $(replace(request.target, r"\?.+" => "")) | $t seconds"""
		return response

	catch exception
		printstyled("ERROR $(replace(request.target, r"\?.+" => "")) | $exception\n",color=:red)
		return HTTP.Response(500)
	end
end

gate(polygon::Vector{<:Any},data::DataFrame,embedding::NamedTuple,selections::NamedTuple;kwargs...) = gate(convert(Vector{SVector{2,Float64}},polygon),data::DataFrame,embedding::NamedTuple,selections::NamedTuple;kwargs...)
function gate(polygon::Vector{SVector{2,Float64}},data::DataFrame,embedding::NamedTuple,selections::NamedTuple;channelRange=range(-3,7,length=50))

    inside(x::SVector{2,Float64}) = inpolygon(x,polygon;in=true,on=false,out=false)
    counts = combine(

        ############## for selected events in embedded polygon
        view( data, selections.rows .& map( inside, embedding.coordinates ), : ),

        ############## histogram per channel
        map( channel -> channel => ( x->fit(Histogram,filter(xi->~ismissing(xi),x),channelRange).weights ) => channel, Base.names(data) ))

    counts.channelValue = ( channelRange[1:end-1] + channelRange[2:end] ) / 2
    return counts
end

max_stars_repo_name: gszep/flow-explorer

max_issues_repo_name: gszep/flow-explorer

max_forks_repo_name: gszep/flow-explorer


text:
using DashBootstrapComponents, DashHtmlComponents

modal = html_div([
    dbc_formgroup([
        dbc_label("Backdrop:"),
        dbc_radioitems(
            id="backdrop-selector",
            options=[
                Dict("label" => "True (default)", "value" => true),
                Dict("label" => "False", "value" => false),
                Dict("label" => "'static'", "value" => "static"),
            ],
            inline=true,
            value=true,
        ),
    ]),
    dbc_button("Open modal", id="open-backdrop", n_clicks=0),
    dbc_modal(
        [
            dbc_modalheader("Header"),
            dbc_modalbody("Change the backdrop of this modal with the radio buttons"),
            dbc_modalfooter(
                dbc_button(
                    "Close",
                    id="close-backdrop",
                    className="ml-auto",
                    n_clicks=0,
                ),
            ),
        ],
        id="modal-backdrop",
        is_open=false,
    ),
]);

callback!(
    app,
    Output("modal-backdrop", "backdrop"),
    Input("backdrop-selector", "value"),
) do backdrop
    return backdrop
end;

callback!(
    app,
    Output("modal-backdrop", "is_open"),
    Input("open-backdrop", "n_clicks"),
    Input("close-backdrop", "n_clicks"),
    State("modal-backdrop", "is_open"),
) do n1, n2, is_open
    return n1 > 0 || n2 > 0 ? is_open == 0 : is_open
end;

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
using NoRainbows
using Test

@testset "NoRainbows.jl" begin
    # Write your tests here.
end


text:
"""Soil grid types."""
module SoilGrid

using DataFrames

export AbstractSoilGrid, FVGrid, getprofile, setprofile!, addprofile!,
    delprofile!, validategrid

"""USDA soil texture classification."""
const textures = ("sand", "loamy sand",
                  "sandy loam", "silt loam", "loam",
                  "sandy clay loam", "silty clay loam", "clay loam",
                  "sandy clay", "silty clay", "clay")

"""Abstract supertype for all soil grids."""
abstract type AbstractSoilGrid end

"""A 1-D finite-volume soil grid."""
struct FVGrid <: AbstractSoilGrid
    """Number of grid levels."""
    level::Integer
    """Top boundary depth."""
    top::Real
    """Bottom boundary depth."""
    bottom::Real
    """Soil texture."""
    texture::String
    """Names of soil grid properties."""
    gridfields::Array{String,1}
    """Names of soil physical properties."""
    propfields::Array{String,1}
    """Names of tracers."""
    tracerfields::Array{String,1}
    """A DataFrame to store soil variable profiles."""
    profiles::DataFrames.DataFrame

    # define default values in the constructor
    function FVGrid(level=25, top=0., bottom=1., texture="none",
                    gridfields=Array{String,1}[],
                    propfields=Array{String,1}[],
                    tracerfields=Array{String,1}[])
        gridfields = union(
            gridfields, ["grid_node", "grid_top", "grid_bottom", "grid_size"])

        if level <= 0
            throw(DomainError("Number of grid levels must be positive!"))
        end

        if top == bottom
            throw(ArgumentError("Top and bottom depths cannot be the same!"))
        end

        # ensure top boundary depth < bottom boundary depth
        if top > bottom
            tmp = top
            top = bottom
            bottom = tmp
        end

        # convert field names from String to Symbol
        allnamedfields = Symbol.(
            [sort(gridfields); sort(propfields); sort(tracerfields)])

        profiles = DataFrames.DataFrame(
            repeat([Float64], outer=length(allnamedfields)),
            allnamedfields, level)

        self = new(level, top, bottom, texture, gridfields, propfields,
                   tracerfields, profiles)
        initgrid!(self)
        return self
    end
end

"""
    getprofile(g::AbstractSoilGrid, name::Symbol)

Access the vertical profile of a variable in a soil grid.
"""
function getprofile(g::AbstractSoilGrid, name::Symbol)
    if name in g.profiles.colindex.names
        return g.profiles[name]
    else
        throw(ArgumentError("The soil grid has no attribute of \"$name\"\!"))
    end
end

"""
    setprofile!(g::AbstractSoilGrid, name::Symbol, val)

Set the vertical profile of a variable in a soil grid.
"""
function setprofile!(g::AbstractSoilGrid, name::Symbol, val)
    if name in g.profiles.colindex.names
        g.profiles[name] = val
    else
        throw(ArgumentError("The soil grid has no attribute of \"$name\"\!"))
    end
end

"""
    addprofile!(g::AbstractSoilGrid, name::Symbol[, val=NaN])

Add the vertical profile of a variable in a soil grid.
"""
function addprofile!(g::AbstractSoilGrid, name::Symbol, val=NaN)
    if name in g.profiles.colindex.names
        throw(ArgumentError(
            "The attribute \"$name\" already exists in the soil grid!"))
    else
        g.profiles[name] = val
    end
end

"""
    delprofile!(g::AbstractSoilGrid, name::Symbol)

Remove the vertical profile of a variable in a soil grid.
"""
function delprofile!(g::AbstractSoilGrid, name::Symbol)
    if name in g.profiles.colindex.names
        delete!(g.profiles, name)
    else
        throw(ArgumentError("The soil grid has no attribute of \"$name\"\!"))
    end
end

"""Initialize the vertical levels of FVGrid (face-centered)."""
function initgrid!(g::FVGrid)
    grid_node = exp.(collect(1:g.level) * 5. / g.level - 5.) .*
        (g.bottom - g.top) + g.top
    grid_size = Array{Float64,1}(g.level)
    grid_bottom = Array{Float64,1}(g.level)

    grid_size[1] = 0.5 * (grid_node[1] + grid_node[2])
    grid_size[2:end-1] = 0.5 * (grid_node[3:end] - grid_node[1:end-2])
    grid_size[end] = grid_node[end] - grid_node[end-1]

    grid_bottom[1:end-1] = 0.5 * (grid_node[1:end-1] + grid_node[2:end])
    grid_bottom[end] = grid_node[end] + 0.5 * grid_size[end]
    grid_top = grid_bottom - grid_size

    setprofile!(g, :grid_node, grid_node)
    setprofile!(g, :grid_size, grid_size)
    setprofile!(g, :grid_bottom, grid_bottom)
    setprofile!(g, :grid_top, grid_top)
end

"""
    validategrid(g::AbstractSoilGrid)

Check if the soil grid specifications meet the requirements for calculation.
"""
function validategrid(g::AbstractSoilGrid)
    if !(g.texture in textures)
        throw(ArgumentError("Soil texture \"", g.texture, "\" not supported!"))
    end

    required_fields = (:temp, :moisture, :porosity)
    for col in required_fields
        if !(col in g.profiles.colindex.names)
            throw(ArgumentError("Missing required field \"$col\"!"))
        end
    end
end

end


text:
const _LS_REMOTE_OK_REGEX = Regex("(?<hash>[0-9a-f]{40})\\s+HEAD")
function _remote_HEAD_hash(url::String)
    out = _run("git ls-remote $(url) HEAD 2>&1"; verbose = false, ignorestatus = true)
    m = match(_LS_REMOTE_OK_REGEX, out)
    isnothing(m) ? "" : m[:hash]
end
_remote_HEAD_hash(gl::GitLink) = _remote_HEAD_hash(remote_url(gl))


text:
export Space, domainspace, rangespace, maxspace,Space,conversion_type, transform,
            itransform, SequenceSpace, ConstantSpace

# Space maps the Domain to the type R
# For example, we have
#   Chebyshev{Interval{Float64}} <: Space{Interval{Float64},Float64}
#   Laurent{PeriodicSegment{Float64}} <: Space{PeriodicSegment{Float64},ComplexF64}
#   Fourier{Circle{ComplexF64}} <: Space{Circle{ComplexF64},Float64}
# Note for now Space doesn't contain any information about the coefficients

abstract type Space{D,R} end

const RealSpace = Space{D,R} where {D,R<:Real}
const ComplexSpace = Space{D,R} where {D,R<:Complex}
const UnivariateSpace = Space{D,R} where {D<:Domain{<:Number},R}
const BivariateSpace = Space{D,R}  where {D<:EuclideanDomain{2},R}
const RealUnivariateSpace = RealSpace{D,R} where {D<:Domain{<:Number},R<:Real}

eltype(S::Space{T}) where {T} = error("Eltype has been changed to domaintype, rangetype or prectype")
eltype(::Type{Space{D,R}}) where {D,R} = error("Eltype has been changed to domaintype, rangetype or prectype")

domaintype(::Space{D,R}) where {D,R} = D
domaintype(::Type{Space{D,R}}) where {D,R} = D
domaintype(::Type{FT}) where {FT<:Space} = domaintype(supertype(FT))
rangetype(::Space{D,R}) where {D,R} = R
rangetype(::Type{Space{D,R}}) where {D,R} = R
rangetype(::Type{FT}) where {FT<:Space} = rangetype(supertype(FT))

domaindimension(sp::Space) = dimension(domain(sp))
dimension(::Space) = ∞  # We assume infinite-dimensional spaces

# add indexing for all spaces, not just DirectSumSpace
# mimicking scalar vs vector

#supports broadcasting, overloaded for ArraySpace
size(::Space) = ()

transpose(sp::Space) = sp  # default no-op

# the default is all spaces have one-coefficient blocks
blocklengths(S::Space) = Ones{Int}(dimension(S))
nblocks(S::Space) = length(blocklengths(S))
block(S::Space,k) = Block(k)

Space(s::Space) = s

abstract type AmbiguousSpace <: Space{AnyDomain,UnsetNumber} end
domain(::AmbiguousSpace) = AnyDomain()

function setdomain(sp::Space{D},d::D) where D<:Domain
    S = typeof(sp)
    @assert length(fieldnames(S))==1
    S(d)
end

# function setdomain(sp::Space,d::Domain)
#     S=typeof(sp)
#     @assert length(fieldnames(S))==1
#     # the domain is not compatible, but maybe we c
#     # can drop the space depence.  For example,
#     # CosSpace{Circle{Float64}} -> CosSpace
#     eval(Meta.parse(string(S.name.module)*"."*string(S.name)))(d)
# end

setcanonicaldomain(s) = setdomain(s,canonicaldomain(s))
reverseorientation(S::Space) = setdomain(S,reverseorientation(domain(S)))

# UnsetSpace dictates that an operator is not defined until
#   its domainspace is promoted
# NoSpace is used to indicate no space exists for, e.g.,
# conversion_type

struct UnsetSpace <: AmbiguousSpace end
struct NoSpace <: AmbiguousSpace end

isambiguous(_) = false
isambiguous(::Type{UnsetNumber}) = true
isambiguous(::Type{Array{T}}) where {T} = isambiguous(T)
isambiguous(sp::Space) = isambiguous(rangetype(sp))

#TODO: should it default to canonicalspace?
points(d::Space,n) = points(domain(d),n)
points(d::Space) = points(d, dimension(d))

canonicalspace(T) = T
canonicaldomain(S::Space) = canonicaldomain(domain(S))

# Check whether spaces are the same, override when you need to check parameters
# This is used in place of == to support AnyDomain
spacescompatible(f::D,g::D) where D<:Space = error("Override spacescompatible for "*string(D))
spacescompatible(::UnsetSpace,::UnsetSpace) = true
spacescompatible(::NoSpace,::NoSpace) = true
spacescompatible(f,g) = false
==(A::Space,B::Space) = spacescompatible(A,B) && domain(A) == domain(B)
spacesequal(A::Space,B::Space) = A==B

pointscompatible(f,g) = spacescompatible(f,g)

# check a list of spaces for compatibility
for OP in (:spacescompatible,:domainscompatible,:spacesequal),TYP in (:AbstractArray,:Tuple)
    @eval function $OP(v::$TYP)
        for k=1:length(v)-1
            if !$OP(v[k],v[k+1])
                return false
            end
        end
        true
    end
end

domain(A::Space) = A.domain # assume it has a field domain

for op in (:tocanonical,:fromcanonical,:tocanonicalD,:fromcanonicalD,:invfromcanonicalD)
    @eval ($op)(sp::Space,x...)=$op(domain(sp),x...)
end

mappoint(a::Space,b::Space,x)=mappoint(domain(a),domain(b),x)
mappoint(a::Space,b::Domain,x)=mappoint(domain(a),b,x)
mappoint(a::Domain,b::Space,x)=mappoint(a,domain(b),x)

for FUNC in (:conversion_rule,:maxspace_rule,:union_rule)
    @eval begin
        function $FUNC(a,b)::Any
            if spacescompatible(a,b)
                a
            else
                NoSpace()
            end
        end
    end
end

for FUNC in (:conversion_type,:maxspace)
    @eval begin
        $FUNC(a::UnsetSpace,b::UnsetSpace) = a
        $FUNC(a::UnsetSpace,b::Space) = b
        $FUNC(a::Space,b::UnsetSpace) = a
    end
end

# gives a space c that has a banded conversion operator TO a and b
function conversion_type(a,b)::Any
    if spacescompatible(a,b)
        a
    elseif !domainscompatible(a,b)
        NoSpace()  # this avoids having to check eachtime
    else
        cr=conversion_rule(a,b)
        cr==NoSpace() ? conversion_rule(b,a) : cr
    end
end

# gives a space c that has a banded conversion operator FROM a and b
maxspace(a,b) = NoSpace()  # TODO: this fixes weird bug with Nothing
function maxspace(a::Space, b::Space)::Any
    if spacescompatible(a,b)
        return a
    elseif !domainscompatible(a,b)
        return NoSpace()  # this avoids having to check eachtime
    end

    cr=maxspace_rule(a,b)
    if !isa(cr,NoSpace)
        return cr
    end

    cr=maxspace_rule(b,a)
    if !isa(cr,NoSpace)
        return cr
    end

    cr=conversion_type(a,b)
    if cr==a
        return b
    elseif cr ==b
        return a
    end

    # check if its banded through canonicalspace
    cspa=canonicalspace(a)
    if spacescompatible(cspa,b)
        # we can't call maxspace(cspa,a)
        # maxspace/conversion_type should be implemented for canonicalspace
        error("Override conversion_type or maxspace for "*string(a)*" and "*string(b))
    end
    if cspa != a && maxspace(cspa,a)==cspa
        return maxspace(b,cspa)
    end

    cspb=canonicalspace(b)
    if spacescompatible(cspb,a)
        # we can't call maxspace(cspb,b)
        error("Override conversion_type or maxspace for "*string(a)*" and "*string(b))
    end
    if cspb !=b && maxspace(cspb,b)==cspb
        return maxspace(a,cspb)
    end

    NoSpace()
end

# union combines two spaces
# this is used primarily for addition of two funs
# that may be incompatible
union(a::AmbiguousSpace, b::AmbiguousSpace) = b
union(a::AmbiguousSpace, b::Space) = b
union(a::Space, b::AmbiguousSpace) = a

function union_by_union_rule(a::Space,b::Space)::Any
    if spacescompatible(a,b)
        if isambiguous(domain(a))
            return b
        else
            return a
        end
    end

    cr = union_rule(a,b)
    cr isa NoSpace || return cr

    union_rule(b,a)
end

function union(a::Space, b::Space)::Any
    cr = union_by_union_rule(a,b)
    cr isa NoSpace || return cr

    cspa=canonicalspace(a)
    cspb=canonicalspace(b)
    if cspa!=a || cspb!=b
        cr = union_by_union_rule(cspa,cspb)
    end
    # TODO: Uncomment when Julia bug is fixed
    cr=maxspace(a,b)  #Max space since we can convert both to it
    if !isa(cr,NoSpace)
        return cr
    end

    a ⊕ b
end

union(a::Space) = a

union(a::Space,b::Space,c::Space) = union(union(a,b),c)
union(a::Space,b::Space,c::Space,d::Space...) =
    union(union(a,b),c,d...)

# tests whether a Conversion operator exists
hasconversion(a,b) = maxspace(a,b) == b

# tests whether a coefficients can be converted to b
isconvertible(a,b) = a == b || hasconversion(a,b)

## Conversion routines
#       coefficients(v::AbstractVector,a,b)
# converts from space a to space b
#       coefficients(v::Fun,a)
# is equivalent to coefficients(v.coefficients,v.space,a)
#       coefficients(v::AbstractVector,a,b,c)
# uses an intermediate space b

coefficients(f,sp1,sp2,sp3) = coefficients(coefficients(f,sp1,sp2),sp2,sp3)

coefficients(f::AbstractVector,::Type{T1},::Type{T2}) where {T1<:Space,T2<:Space} =
    coefficients(f,T1(),T2())
coefficients(f::AbstractVector,::Type{T1},sp2::Space) where {T1<:Space} = coefficients(f,T1(),sp2)
coefficients(f::AbstractVector,sp1::Space,::Type{T2}) where {T2<:Space} = coefficients(f,sp1,T2())

## coefficients defaults to calling Conversion, otherwise it tries to pipe through Chebyshev

function defaultcoefficients(f,a,b)
    ct=conversion_type(a,b) # gives a space that has a banded conversion to both a and b

    if spacescompatible(a,b)
        f
    elseif hasconversion(a,b)
        mul_coefficients(Conversion(a,b),f)
    elseif hasconversion(b,a)
        ldiv_coefficients(Conversion(b,a),f)
    else
        csp=canonicalspace(a)

        if spacescompatible(a,csp)# a is csp, so try b
            csp=canonicalspace(b)
        end
        if spacescompatible(a,csp) || spacescompatible(b,csp)
            # b is csp too, so we are stuck, try Fun constructor
            coefficients(default_Fun(Fun(a,f),b))
        else
            coefficients(f,a,csp,b)
        end
    end
end

coefficients(f,a,b) = defaultcoefficients(f,a,b)

## rand
# checkpoints is used to give a list of points to double check
# the expansion
rand(d::Space,k...) = rand(domain(d),k...)
checkpoints(d::Space) = checkpoints(domain(d))

## default transforms
abstract type AbstractTransformPlan{T} <: Plan{T} end

space(P::AbstractTransformPlan) = P.space

# These plans are use to wrap another plan
for Typ in (:TransformPlan,:ITransformPlan)
    @eval begin
        struct $Typ{T,SP,inplace,PL} <: AbstractTransformPlan{T}
            space::SP
            plan::PL
        end
        $Typ(space,plan,::Type{Val{inplace}}) where {inplace} =
            $Typ{eltype(plan),typeof(space),inplace,typeof(plan)}(space,plan)
        # *(P::$Typ, x::AbstractArray) = P.plan*x
    end
end

for Typ in (:CanonicalTransformPlan,:ICanonicalTransformPlan)
    @eval begin
        struct $Typ{T,SP,PL,CSP} <: AbstractTransformPlan{T}
            space::SP
            plan::PL
            canonicalspace::CSP
        end
        $Typ(space,plan,csp) =
            $Typ{eltype(plan),typeof(space),typeof(plan),typeof(csp)}(space,plan,csp)
        $Typ(::Type{T},space,plan,csp) where {T} =
            $Typ{T,typeof(space),typeof(plan),typeof(csp)}(space,plan,csp)
    end
end

# Canonical plan uses coefficients
function CanonicalTransformPlan(space,v)
    csp = canonicalspace(space)
    CanonicalTransformPlan(eltype(v),space,plan_transform(csp,v),csp)
end
function plan_transform(sp::Space,vals)
    csp = canonicalspace(sp)
    if sp == csp
        error("Override for $sp")
    end
    CanonicalTransformPlan(sp,plan_transform(csp,vals),csp)
end

function ICanonicalTransformPlan(space,v)
    csp = canonicalspace(space)
    cfs = coefficients(v,space,csp)
    ICanonicalTransformPlan(eltype(v),space,plan_itransform(csp,cfs),csp)
end
function plan_itransform(sp::Space,v)
    csp = canonicalspace(sp)
    if sp == csp
        error("Override for $sp")
    end
    cfs = coefficients(v,sp,csp)
    ICanonicalTransformPlan(sp,plan_itransform(csp,cfs),csp)
end

plan_transform!(sp::Space,vals) = error("Override for $sp")
plan_itransform!(sp::Space,cfs) = error("Override for $sp")

# transform converts from values at points(S,n) to coefficients
# itransform converts from coefficients to values at points(S,n)

transform(S::Space,vals) = plan_transform(S,vals)*vals
itransform(S::Space,cfs) = plan_itransform(S,cfs)*cfs

itransform!(S::Space,cfs) = plan_itransform!(S,cfs)*cfs
transform!(S::Space,cfs) = plan_transform!(S,cfs)*cfs

*(P::CanonicalTransformPlan,vals::AbstractVector) = coefficients(P.plan*vals,P.canonicalspace,P.space)
*(P::ICanonicalTransformPlan,cfs::AbstractVector) = P.plan*coefficients(cfs,P.space,P.canonicalspace)

for OP in (:plan_transform,:plan_itransform,:plan_transform!,:plan_itransform!)
    # plan transform expects a vector
    # this passes an empty Float64 array
    @eval begin
        $OP(S::Space,::Type{T},n::Integer) where {T} = $OP(S,Vector{T}(undef, n))
        $OP(S::Space,n::Integer) = $OP(S, Float64, n)
    end
end

## sorting
# we sort spaces lexigraphically by default

for OP in (:<,:(<=),:>,:(>=),:(isless))
    @eval $OP(a::Space,b::Space)=$OP(string(a),string(b))
end

## Important special spaces

struct ZeroSpace{DD,R} <: Space{DD,R}
    domain::DD
    ZeroSpace{DD,R}(d::DD) where {DD,R} = new(d)
    ZeroSpace{DD,R}(d::AnyDomain) where {DD,R} = new(convert(DD,d))
end

ZeroSpace(S::Space) = ZeroSpace{domaintype(S),rangetype(S)}(domain(S))
ZeroSpace() = ZeroSpace{AnyDomain,UnsetNumber}(AnyDomain())
domain(S::ZeroSpace) = S.domain

dimension(::ZeroSpace) = 0

spacescompatible(::ZeroSpace,::ZeroSpace) = true
for FUNC in (:conversion_type,:maxspace)
    @eval begin
        $FUNC(::ZeroSpace,::UnsetSpace) = UnsetSpace()
        $FUNC(::UnsetSpace,::ZeroSpace) = UnsetSpace()
    end
end

"""
`ConstantSpace` is the 1-dimensional scalar space.
"""
struct ConstantSpace{DD,R} <: Space{DD,R}
    domain::DD
    ConstantSpace{DD,R}(d::DD) where {DD,R} = new(d)
    ConstantSpace{DD,R}(d::AnyDomain) where {DD,R} = new(convert(DD,d))
end

ConstantSpace(d::Domain) = ConstantSpace{typeof(d),real(prectype(d))}(d)

ConstantSpace(::Type{N},d::Domain) where {N<:Number} = ConstantSpace{typeof(d),real(N)}(d)
ConstantSpace(::Type{N}) where {N<:Number} = ConstantSpace(N,AnyDomain())
ConstantSpace() = ConstantSpace(Float64)

convert(::Type{Space}, z::Number) = ConstantSpace(convert(Domain, z))  # Spaces
convert(::Type{ConstantSpace}, d::Domain) = ConstantSpace(d)
Space(z::Number) = convert(Space, z)

isconstspace(::ConstantSpace) = true

for pl in (:plan_transform,:plan_transform!,:plan_itransform,:plan_itransform!)
    @eval $pl(sp::ConstantSpace,vals::AbstractVector) = I
end

# we override maxspace instead of maxspace_rule to avoid
# domainscompatible check.
for OP in (:maxspace,:(union))
    @eval begin
        $OP(A::ConstantSpace{AnyDomain},B::ConstantSpace{AnyDomain}) = A
        $OP(A::ConstantSpace{AnyDomain},B::ConstantSpace) = B
        $OP(A::ConstantSpace,B::ConstantSpace{AnyDomain}) = A
        $OP(A::ConstantSpace,B::ConstantSpace) = ConstantSpace(domain(A) ∪ domain(B))
    end
end

space(x::Number) = ConstantSpace(typeof(x))
space(f::AbstractArray{T}) where T<:Number = ArraySpace(ConstantSpace{T}(), size(f)...)

setdomain(A::ConstantSpace{DD,R}, d) where {DD,R} = ConstantSpace{typeof(d),R}(d)

blocklengths(::ConstantSpace) = Vec(1)

# Range type is Nothing since function evaluation is not defined
struct SequenceSpace <: Space{PositiveIntegers,Nothing} end

"""
`SequenceSpace` is the space of all sequences, i.e., infinite vectors.
Also denoted ℓ⁰.
"""
SequenceSpace()

const ℓ⁰ = SequenceSpace()
dimension(::SequenceSpace) = ∞
domain(::SequenceSpace) = ℕ
spacescompatible(::SequenceSpace,::SequenceSpace) = true

## Boundary

boundary(S::Space) = boundary(domain(S))


text:
export StagedRocket

"""
A staged rocket is made up of multiple stages which are rockets themselves
"""
mutable struct StagedRocket <: Rocket
    nextstage::Rocket
    tank::Tank
    engine::Engine
end

function Rocket(payload::Rocket, tank::Tank, engine::Engine)
    StagedRocket(payload, tank, engine)
end

"""
    Rocket(tank, engine)
Create a single stage rocket
"""
function Rocket(tank::Tank, engine::Engine)
    StagedRocket(emptypayload, tank, engine)
end

"""
    mass(r::StagedRocket)
Total mass of all stages in a staged rocket
"""
function mass(r::StagedRocket)
    mass(r.nextstage) + mass(r.tank) + mass(r.engine)
end

"""
    update!(r::StagedRocket)
Consume propellant for current active stage (bottom stage and any active boosters).
"""
function update!(r::StagedRocket, t::Number, Δt::Number)
    mflow = mass_flow(thrust(r), Isp(r.engine))
	consume!(r.tank, mflow * Δt)
end

"""
    thrust(r::StagedRocket)
Total thrust of active rocket stage
"""
thrust(r::StagedRocket) = thrust(r.engine)


text:
@testset "prefix" begin
    version = VERSION
    text = repeat("foobar\n", 100)
    html = string(PlutoStaticHTML.State(text)) * text
    state = PlutoStaticHTML.State(html)
    @test state.input_sha == PlutoStaticHTML.sha(html)
    @test state.julia_version == string(VERSION)
end

function try_read(file)::String
    for i in 1:100
        try
            return read(file, String)
        catch
            sleep(0.2)
        end
    end
    # Throw original error if still no success.
    read(file, String)
end

function try_rm(file)
    for i in 1:100
        try
            return rm(file)
        catch
            sleep(0.2)
        end
    end
    # Throw original error if still no success.
    rm(file)
end

@testset "caching" begin
    dir = mktempdir()
    use_distributed = false

    cd(dir) do
        path(name) = joinpath(dir, "$name.txt")
        code = pluto_notebook_content("""
            begin
                path = "$(path('a'))"
                println("Writing to \$path")
                write(path, "a")
            end
            """)
        write("a.jl", code)

        code = pluto_notebook_content("""write("$(path('b'))", "b")""")
        write("b.jl", code)

        bo = BuildOptions(dir; use_distributed)
        build_notebooks(bo)

        # Without try_read, Pluto in another process may still have a lock on a txt file.
        @test try_read("a.txt") == "a"
        try_read("a.html")
        @test try_read("b.txt") == "b"
        try_read("b.html")

        try_rm("a.html")
        try_rm("a.txt")
        try_rm("b.txt")

        previous_dir = dir
        dir = mktempdir()

        cd(dir) do
            cp(joinpath(previous_dir, "a.jl"), joinpath(dir, "a.jl"))
            cp(joinpath(previous_dir, "b.jl"), joinpath(dir, "b.jl"))

            bo = BuildOptions(dir; use_distributed, previous_dir)
            build_notebooks(bo)

            # a was evaluated because "a.html" was removed.
            # note that pluto always writes txt files to the first dir.
            @test try_read(joinpath(previous_dir, "a.txt")) == "a"
            try_read("a.html")

            # b was not evaluated because "b.html" was used from the cache.
            @test !isfile(joinpath(previous_dir, "b.txt"))
            try_read("b.html")
        end
    end
end

@testset "extract_previous" begin
    text = """
            prefix

            <!-- PlutoStaticHTML.Begin -->
            <!-- PlutoStaticHTML.End -->
            suffix
            """
    prev = PlutoStaticHTML.Previous(text::String)
    @test startswith(prev.text, PlutoStaticHTML.BEGIN_IDENTIFIER)
    @test endswith(prev.text, PlutoStaticHTML.END_IDENTIFIER)
end

@testset "html_cache_output" begin
    # Test whether the HTML copied from the cache doesn't contain anything outside PlutoStaticHTML.Begin and End.
    dir = mktempdir()

    cd(dir) do
        path = joinpath(dir, "notebook.jl")
        code = pluto_notebook_content("x = 1")
        write(path, code)

        output_format = html_output
        use_distributed = false
        bo = BuildOptions(dir; output_format, use_distributed)
        build_notebooks(bo)

        output_path = joinpath(dir, "notebook.html")
        output = read(output_path, String)

        # This is similar to how Franklin, for example, adds a prefix and suffix.
        write(output_path, "prefix\n" * output * "\nsuffix")

        previous_dir = dir
        bo = BuildOptions(dir; output_format, use_distributed, previous_dir)
        build_notebooks(bo)

        output2 = read(output_path, String)

        @test output == output2
    end
end

@testset "franklin_markdown_cache_output" begin
    # Test whether the Franklin Markdown copied from the cache is correct.
    dir = mktempdir()

    cd(dir) do
        path = joinpath(dir, "notebook.jl")
        code = pluto_notebook_content("""
            md\"\"\"
            +++
            title = \"foo\"
            +++
            \"\"\"
            """)
        write(path, code)

        use_distributed = false
        output_format = franklin_output
        bo = BuildOptions(dir; use_distributed, output_format)
        build_notebooks(bo)

        output_path = joinpath(dir, "notebook.md")
        output = read(output_path, String)

        previous_dir = dir
        bo = BuildOptions(dir; output_format, use_distributed, previous_dir)
        build_notebooks(bo)

        output2 = read(output_path, String)

        @test output == output2
    end
end

@testset "documenter_cache_output" begin
    # Test whether the Franklin Markdown copied from the cache is correct.
    dir = mktempdir()

    cd(dir) do
        path = joinpath(dir, "notebook.jl")
        code = pluto_notebook_content("""
            x = 1
            """)
        write(path, code)

        use_distributed = false
        output_format = documenter_output
        bo = BuildOptions(dir; use_distributed, output_format)
        build_notebooks(bo)

        output_path = joinpath(dir, "notebook.md")
        output = read(output_path, String)

        previous_dir = dir
        bo = BuildOptions(dir; output_format, use_distributed, previous_dir)
        build_notebooks(bo)

        output2 = read(output_path, String)

        @test output == output2
    end
end


text:
mutable struct Haploid <: AbstractAgent
    id::Int
    trait::Float64
end

"""
``` julia
wright_fisher(;
    numagents = 100,
    selection = true
)
```
Same as in [Wright-Fisher model of evolution](@ref).
"""
function wright_fisher(; numagents = 100, selection = true)
    model = ABM(Haploid)
    for i in 1:numagents
        add_agent!(model, rand(model.rng))
    end
    !selection && return model, dummystep, wright_fisher_model_step_neutral!
    return model, dummystep, wright_fisher_model_step_selection!
end

wright_fisher_model_step_neutral!(model::ABM) = sample!(model, nagents(model))

wright_fisher_model_step_selection!(model::ABM) = sample!(model, nagents(model), :trait)


text:
function MannWhitneyUTest()
    return MannWhitneyUTest(())
end

function MannWhitneyUTest(arg0::NaNStrategy, arg1::TiesStrategy)
    return MannWhitneyUTest((NaNStrategy, TiesStrategy), arg0, arg1)
end

function mann_whitney_u(obj::MannWhitneyUTest, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(obj, "mannWhitneyU", jdouble, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function mann_whitney_u_test(obj::MannWhitneyUTest, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(obj, "mannWhitneyUTest", jdouble, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function mann_whitney_u_test(obj::MannWhitneyUTest, arg0::Vector{jdouble}, arg1::Vector{jdouble}, arg2::jboolean)
    return jcall(obj, "mannWhitneyUTest", jdouble, (Vector{jdouble}, Vector{jdouble}, jboolean), arg0, arg1, arg2)
end


text:
"""
    count_nucleotides(strand)

The frequency of each nucleotide within `strand` as a dictionary.

Invalid strands raise a `DomainError`.

"""
function count_nucleotides(strand)
    d = Dict([(ch, 0) for ch in "ACGT"])
    for ch in strand
        if ch in keys(d)
            d[ch] += 1
        else
            throw(DomainError(ch))
        end
    end
    return d
end

max_stars_repo_name: twolodzko/Learning

max_issues_repo_name: twolodzko/Learning

max_forks_repo_name: twolodzko/Learning


text:
# test_to_s
sandboxed_test("diff", "test patch") do test_repo, test_repo_path
    t1 = GitTree(test_repo[Oid("d70d245ed97ed2aa596dd1af6536e4bfdb047b69")])
    t2 = GitTree(test_repo[Oid("7a9e0b02e63179929fed24f0a3e0f19168114d10")])
    d  = diff(test_repo, t1, t2, {:context_lines => 0})
    s1 = string(patches(d)[1])
    p1 = replace("diff --git a/another.txt b/another.txt
index 3e5bcba..546c735 100644
--- a/another.txt
+++ b/another.txt
@@ -2 +2 @@ Git is fast. With Git, nearly all operations are performed locally, giving
-it a huge speed advantage on centralized systems that constantly have to
+it an huge speed advantage on centralized systems that constantly have to
@@ -11,4 +10,0 @@ from the start.
-Let's see how common operations stack up against Subversion, a common
-centralized version control system that is similar to CVS or
-Perforce. Smaller is faster.
-
@@ -34,0 +31,4 @@ SVN.
+Let's see how common operations stack up against Subversion, a common
+centralized version control system that is similar to CVS or
+Perforce. Smaller is faster.
+
", "\r\n", "\n")
    @test s1 == p1

    s2 = string(patches(d)[2])
    p2 = replace("diff --git a/readme.txt b/readme.txt
index 7b808f7..29ab705 100644
--- a/readme.txt
+++ b/readme.txt
@@ -1 +1 @@
-The Git feature that really makes it stand apart from nearly every other SCM
+The Git feature that r3ally mak3s it stand apart from n3arly 3v3ry other SCM
@@ -10,4 +9,0 @@ This means that you can do things like:
-Frictionless Context Switching. Create a branch to try out an idea, commit a
-few times, switch back to where you branched from, apply a patch, switch
-back to where you are experimenting, and merge it in.
-
@@ -27,3 +22,0 @@ Notably, when you push to a remote repository, you do not have to push all
-of your branches. You can choose to share just one of your branches, a few
-of them, or all of them. This tends to free people to try new ideas without
-worrying about having to plan how and when they are going to merge it in or
@@ -35 +28 @@ incredibly easy and it changes the way most developers work when they learn
-it.
+it.!
\\ No newline at end of file
", "\r\n", "\n")
    @test s2 == p2
end


text:
"Returns the time-scale for advection on a regular grid across a single grid cell
 for ShallowWaterModel."

import Oceananigans.Utils: cell_advection_timescale

function shallow_water_cell_advection_timescale(uh, vh, h, grid)
    uhmax = maximum(abs, uh)
    vhmax = maximum(abs, vh)
    hmin  = minimum(abs,  h)

    return min(grid.Δx / uhmax, grid.Δy / vhmax) * hmin
end

cell_advection_timescale(model::ShallowWaterModel) =
    shallow_water_cell_advection_timescale(
        model.solution.uh.data.parent,
        model.solution.vh.data.parent,
        model.solution.h.data.parent,
        model.grid
        )


text:
export run_berl
using YAML

function run_berl()
    runs::Int64=0
    # Algorithms
    algs = []
    for (k, v) in pairs(YAML.load_file("./run_config/algorithms.yaml"))
        if k=="runs"
            runs = v
        elseif v
            push!(algs, k)
        end
    end

    # Environments
    envs=[]
    for (k, v) in pairs(YAML.load_file("./run_config/environments.yaml"))
        if v
            push!(envs, k)
        end
    end

    println("Algorithms: ", algs)
    println("Environments: ", envs)
    for e in envs
        for a in algs
            ids::Array{String}=[]
            for r in 1:runs
                if e in gym_envs
                    d = start_berl(a, "gym"; gym_env=e)
                elseif e in atari_envs
                    d = start_berl(a, "atari_ram"; atari_game=e)
                else
                    d = start_berl(a, e)
                end
                d["run"]=r
                print_metrics(d)
                push!(ids, d["id"])
            end
            try
                BERLplot(ids; name="benchmark_$(a)_$e.svg")
            catch

            end
        end
    end
end


text:
@testset "assemble" begin
    dofs = [1, 3, 5, 7]

    # residual
    ge = rand(4)
    g = zeros(8)
    assemble!(g, dofs, ge)
    @test g[1] == ge[1]
    @test g[3] == ge[2]
    @test g[5] == ge[3]
    @test g[7] == ge[4]

    # stiffness
    a = start_assemble()
    Ke = rand(4, 4)
    assemble!(a, dofs, Ke)
    K = end_assemble(a)
    @test K[1,1] == Ke[1,1]
    @test K[1,5] == Ke[1,3]
    @test K[5,1] == Ke[3,1]
end


text:
# ###################################################################################################
# Simulation studies for one-dimensional Ising model
#
# author: weiya
# date: Dec 25, 2018
#
# refer to Section 2.4.2 & Section 10.4 of
# Liu, J. S. (2008). Monte Carlo strategies in scientific computing (2. ed). New York, NY: Springer.
#
# ###################################################################################################
using Distributions
using StatsBase
using Plots
using LaTeXStrings

function exactMethod(M = 2000; beta = 10, d = 100)
    V = ones(d)
    V[1] = exp(beta) + exp(-beta)
    for t = 2:d
        V[t] = V[1] * V[t-1]
    end
    Z = 2*V[d]
    x = ones(Int, d)
    x[d] = 2*rand(Bernoulli(1/2))-1
    for t = (d-1):1
#        p1 = V[t] * exp(x[d] * 1)
#        p2 = V[t] * exp(x[d] * (-1))
        p1 = exp(x[d])
        p2 = exp(-x[d])
        x[t] = 2*rand(Bernoulli(p1/(p1+p2)))-1
    end
    return x
end

function pdfIsing(x, T)
    n = length(x)
    y = sum(x[1:(n-1)] .* x[2:n])
    return exp(y/T)
end

function logpdfIsing(x, T)
    n = length(x)
    y = sum(x[1:(n-1)] .* x[2:n])
    return y/T
end

function sampleX!(x::Array{Int}, T::Float64)
    d = length(x)
    for i = 1:d
        if i == 1
            p1 = exp( x[i+1] / T )
            p2 = exp( -x[i+1] / T )
        elseif i == d
            p1 = exp( x[i-1] / T )
            p2 = exp( -x[i-1] / T )
        else
            p1 = exp( (x[i-1] + x[i+1]) / T )
            p2 = exp( (-x[i-1] - x[i+1]) / T )
        end
        x[i] = 2*rand(Bernoulli( p1 / (p1+p2) )) - 1
    end
end

function sampleX!(x::Array{Int,2}, T::Array{Float64})
    I = length(T)
    d = size(x, 1)
    for j = 1:I
        for i = 1:d
            if i == 1
                p1 = exp( x[i+1,j] / T[j] )
                p2 = exp( -x[i+1,j] / T[j] )
            elseif i == d
                p1 = exp( x[i-1,j] / T[j] )
                p2 = exp( -x[i-1,j] / T[j] )
            else
                p1 = exp( (x[i-1,j] + x[i+1,j]) / T[j] )
                p2 = exp( (-x[i-1,j] - x[i+1,j]) / T[j] )
            end
            x[i,j] = 2*rand(Bernoulli( p1 / (p1+p2) )) - 1
        end
    end
end

function singleMCMC(M = 2000; T = 0.1, d = 100)
    # intial x
    X = ones(Int, M, d)
    x = 2 * rand(Bernoulli(1/2), d) .- 1
    for m = 1:M
        sampleX!(x, T)
        X[m, :] .= x
    end
    return X
end
# result: always trap in one of the modes

function parallelTemper(M = 2000, T = [0.1, 0.2, 0.3, 0.4]; d = 100, alpha0 = T[1])
    I = length(T)
    # initial
    x = 2*rand(Bernoulli(1/2), d, I) .- 1
    num1 = zeros(3)
    num2 = zeros(3)
    res = zeros(M, 4)
    for m = 1:M
        if rand() < alpha0 # parallel step
            #for i = 1:I
                # !!! Doesn't work !!!
                # y1 = x[:,i]
                # sampleX!(x[:,i], T[i])
                # y2 = x[:,i]
                # println(sum(abs.(y2.-y1)))
            #end
            sampleX!(x, T)
        else
#            idx = sample(1:I, 2, replace = false) # not neigbor
            idx1 = sample(1:(I-1))
            num1[idx1] += 1
            idx = [idx1, idx1+1]
#            rho = pdfIsing(x[:,idx[2]], T[idx[1]]) * pdfIsing(x[:,idx[1]], T[idx[2]]) / (pdfIsing(x[:,idx[1]], T[idx[1]]) * pdfIsing(x[:,idx[2]], T[idx[2]]))
            rho = logpdfIsing(x[:,idx[2]], T[idx[1]]) + logpdfIsing(x[:,idx[1]], T[idx[2]]) - (logpdfIsing(x[:,idx[1]], T[idx[1]]) + logpdfIsing(x[:,idx[2]], T[idx[2]]))
            if log(rand()) < rho
                # swappin step
                num2[idx1] += 1
                tmp = copy(x[:, idx[1]])
                x[:,idx[1]] .= x[:, idx[2]]
                x[:,idx[2]] .= tmp
            end
        end
        res[m, :] = sum(x, dims=1)
    end
    return res, num2 ./ num1
end

# run
for i in 1:10
    res, ratio = parallelTemper(200000)

    p1 = histogram(res[190000:end, 1], nbins = 20, title = L"T=0.1")
    p2 = histogram(res[190000:end, 2], nbins = 20, title = L"T=0.2")
    p3 = histogram(res[190000:end, 3], nbins = 20, title = L"T=0.3")
    p4 = histogram(res[190000:end, 4], nbins = 20, title = L"T=0.4")
    plot(p1, p2, p3, p4, layout = (2, 2),legend = false)
    savefig("IsingHist-rep$i.png")

    plot(autocor(res), label = [L"T=0.1", L"T=0.2", L"T=0.3", L"T=0.4"], title = "ACF")
    savefig("IsingACF-rep$i.png")
end

res, = parallel(d=3)
X = singleMCMC(d=3)
p1 = histogram(res[:,1])
p2 = histogram(X[:,1], nbins = 2)
plot(p1, p2, layout = (1, 2), legend = false)

max_stars_repo_name: szcf-weiya/MCMC

max_issues_repo_name: szcf-weiya/MCMC

max_forks_repo_name: szcf-weiya/MCMC


text:
# module MySQL_Octo.SQL

function prepare(adapter::QueryBuilders.AbstractAdapter, from::QueryBuilders.FromClause)
    tablename = from.tablename
    mysql_stmt_prepare(adapter.handle, "SELECT * FROM $tablename LIMIT 0")
    mysql_metadata(adapter.handle)
end


text:
# ElastoPlasticModels.jl

# The Elasto-Plastic alpha function.
function ep_alpha(
    mdl::ElastoPlasticModel,
    z::Number
)
    zba = mdl.breakaway_displacement
    zmax = mdl.max_bristle_displacement
    if abs(z) < zba
        alpha = 0.0
    elseif  abs(z) >= zba && abs(z) <= zmax
        half = 0.5
        arg = pi * (z - ((zmax + zba) / 2.0)) / (zmax - zba)
        alpha = half * (1.0 + sin(arg))
    else
        alpha = 1.0
    end
    return alpha
end

# The Elasto-Plastic model Stribeck equation.
function ep_stribeck(
    mdl::ElastoPlasticModel,
    v::T,
    N::T
) where T <: Number
    fc = mdl.coulomb_coefficient * N
    fs = mdl.static_coefficient * N
    return fc + (fs - fc) * exp(-(v / mdl.stribeck_velocity)^2)
end

# The Elasto-Plastic state equation (bristle velocity)
function heuristic_state_equation(
    mdl::ElastoPlasticModel,
    t::T,
    nrm::T,
    pos::T,
    vel::T,
    z::Array{T}
) where T <: Number
    alpha = ep_alpha(mdl, z[1])
    fss = ep_stribeck(mdl, vel, nrm)
    zdot = vel * (1.0 -
        alpha * mdl.bristle_stiffness * sign(vel) * z[1] / fss)
    return [zdot]
end

# The Elasto-Plastic friction force equation
function heuristic_force_equation(
    mdl::ElastoPlasticModel,
    t::T,
    nrm::T,
    pos::T,
    vel::T,
    z::Array{T},
    dzdt::Array{T}
) where T <: Number
    return mdl.bristle_stiffness * z[1] + mdl.bristle_damping * dzdt[1] +
        mdl.viscous_damping * vel
end

# ------------------------------------------------------------------------------
# Required routines to support friction model fitting.

function model_from_array(mdl::ElastoPlasticModel, x::Array{T}) where T <: Number
    ElastoPlasticModel(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8])
end

function model_to_array(x::ElastoPlasticModel)
    [
        x.static_coefficient,
        x.coulomb_coefficient,
        x.stribeck_velocity,
        x.bristle_stiffness,
        x.bristle_damping,
        x.breakaway_displacement,
        x.max_bristle_displacement,
        x.viscous_damping
    ]
end