text:
#
"""
Solve boundary value problem
A * u = f

where A is a linear operator
we get

M * A * uh = M * ( B * f - A * ub)

where u = uh + ub (Dirichlet data)

learn to add neumann, robin data, and solve BVP

plug in to SciMLBase.BVProblem
"""
function makeRHS(prob::BoundaryValueProblem)

    rhs = b - applyBC()
end

function SciMLBase.solve(prob::BoundaryValueProblem)
    BoundaryValueCache()
end

"""
lhs(u) = b

move boundary data to right hand side, and solve BVP A \ b
using LinearSolve.jl
"""
struct BoundaryValuePDEProblem
    u0::Tu0
    bc::Tbc

    lhs::Tlhs # op or func

    rhs_func::Trhs
    b

    space
end

struct BoundaryValuePDECache
    u0::Tu
    bc::Tbc

    lhs::Tlhs

    rhs_func
    b::Tu

    space

    linprob
end

function SciMLBase.init(prob::BoundaryValuePDEProblem, alg::SciMLBase.AbstractBoundaryValueAlgorithm)
end

#


text:
# ---
# title: 477. Total Hamming Distance
# id: problem477
# author: Indigo
# date: 2021-06-22
# difficulty: Medium
# categories: Bit Manipulation
# link: <https://leetcode.com/problems/total-hamming-distance/description/>
# hidden: true
# ---
#
# The [Hamming distance](https://en.wikipedia.org/wiki/Hamming_distance) between
# two integers is the number of positions at which the corresponding bits are
# different.
#
# Now your job is to find the total Hamming distance between all pairs of the
# given numbers.
#
# **Example:**
#
#
#
#     Input: 4, 14, 2
#
#     Output: 6
#
#     Explanation: In binary representation, the 4 is 0100, 14 is 1110, and 2 is 0010 (just
#     showing the four bits relevant in this case). So the answer will be:
#     HammingDistance(4, 14) + HammingDistance(4, 2) + HammingDistance(14, 2) = 2 + 2 + 2 = 6.
#
#
# **Note:**
#
#   1. Elements of the given array are in the range of `0 ` to `10^9`
#   2. Length of the array will not exceed `10^4`.
#
#
## @lc code=start
using LeetCode

function total_hamming_distance(nums)
    res, n = 0, length(nums)
    for i in 0:30
        c = 0
        for num in nums
            c += (num >> i) & 1
        end
        res += c * (n - c)
    end
    return res
end
## @lc code=end

ext: jl


text:
using Test
import MathOptInterface
const MOI = MathOptInterface
const MOIT = MOI.Test
const MOIU = MOI.Utilities

const DoubleDicts = MathOptInterface.Utilities.DoubleDicts

include("../dummy.jl")

remove_moi(x::String) = replace(x, "MathOptInterface." => "")
function compare_without_moi(x::String, y::String)
    @test remove_moi(x) == remove_moi(y)
end

@testset "IndexMap" begin
    map = MOIU.IndexMap()
    @test length(map) == 0
    compare_without_moi(sprint(show, map), "Utilities.IndexMap()")
    x = MOI.VariableIndex(1)
    y = MOI.VariableIndex(2)
    cx = MOI.ConstraintIndex{MOI.SingleVariable, MOI.Integer}(1)
    cy = MOI.ConstraintIndex{MOI.SingleVariable, MOI.Integer}(2)
    map = MOIU.IndexMap(Dict(x => y), DoubleDicts.IndexDoubleDict())
    map.conmap[cx] = cy
    @test length(map) == 2
    # `x=>y` in Julia <= 1.1 and `x => y` in Julia >= 1.2
    if VERSION < v"1.2"
        x_y = string(x) *  "=>" * string(y)
        c_x_y = string(cx) *  "=>" * string(cy)
    else
        x_y = string(x => y)
        c_x_y = string(cx => cy)
    end
    compare_without_moi(sprint(show, map), "Utilities.IndexMap($x_y,$c_x_y)")
end

@testset "AUTOMATIC" begin
    src = DummyModel()
    dest = DummyModel()
    @test_throws ErrorException MOIU.automatic_copy_to(dest, src)
    try
        @test_throws ErrorException MOIU.automatic_copy_to(dest, src)
    catch err
        @test sprint(showerror, err) == "Model DummyModel does not" *
        " support copy with names."
    end
end
@testset "Default" begin
    @test !MOIU.supports_default_copy_to(DummyModel(), false)
    @test !MOIU.supports_default_copy_to(DummyModel(), true)
    model = MOIU.Model{Float64}()
    MOIT.failcopytestc(model)
    MOIT.failcopytestia(model)
    MOIT.failcopytestva(model)
    MOIT.failcopytestca(model)
    MOIT.copytest(model, MOIU.Model{Float64}())
end
@testset "Allocate-Load" begin
    @test !MOIU.supports_allocate_load(DummyModel(), false)
    @test !MOIU.supports_allocate_load(DummyModel(), true)
    mock = MOIU.MockOptimizer(MOIU.Model{Float64}(), needs_allocate_load=true)
    MOIT.failcopytestc(mock)
    MOIT.failcopytestia(mock)
    MOIT.failcopytestva(mock)
    MOIT.failcopytestca(mock)
    MOIT.copytest(mock, MOIU.Model{Float64}())
end

struct DummyEvaluator <: MOI.AbstractNLPEvaluator end

@testset "Create variables in same ordering when NLPBlock is used (#849)" begin
    model = MOIU.UniversalFallback(MOIU.Model{Float64}())
    a = MOI.add_variable(model)
    b, c = MOI.add_variables(model, 2)
    x, cx = MOI.add_constrained_variable(model, MOI.GreaterThan(0.0))
    y, cy = MOI.add_constrained_variables(model, MOI.Nonnegatives(1))
    nlp_data = MOI.NLPBlockData(
        [MOI.NLPBoundsPair(0.0, 1.0) for i in 1:5],
        DummyEvaluator(), false)
    MOI.set(model, MOI.NLPBlock(), nlp_data)
    copy = MOIU.UniversalFallback(MOIU.Model{Float64}())
    index_map = MOIU.default_copy_to(copy, model, true)
    for vi in [a, b, c, x, y[1]]
        @test index_map[vi] == vi
    end
    for ci in [cx, cy]
        @test index_map[ci] == ci
    end
end

# Model that distinguish variables created constrained
struct ConstrainedVariablesModel <: MOI.ModelLike
    allocate_load::Bool
    added_constrained::Vector{Bool}
end
function ConstrainedVariablesModel(allocate_load::Bool)
    return ConstrainedVariablesModel(allocate_load, Bool[])
end
function MOI.empty!(model::ConstrainedVariablesModel)
    empty!(model.added_constrained)
end

function MOIU.supports_default_copy_to(model::ConstrainedVariablesModel, ::Bool)
    return !model.allocate_load
end
function MOIU.supports_allocate_load(model::ConstrainedVariablesModel, ::Bool)
    return model.allocate_load
end
function MOI.copy_to(dest::ConstrainedVariablesModel, src::MOI.ModelLike; kws...)
    MOIU.automatic_copy_to(dest, src; kws...)
end

function MOI.add_variables(model::ConstrainedVariablesModel, n)
    m = length(model.added_constrained)
    for i in 1:n
        push!(model.added_constrained, false)
    end
    return MOI.VariableIndex.(m .+ (1:n))
end
function MOIU.allocate_variables(model::ConstrainedVariablesModel, n)
    return MOI.add_variables(model, n)
end
function MOIU.load_variables(model::ConstrainedVariablesModel, n)
end
function MOI.add_constrained_variables(model::ConstrainedVariablesModel,
                                       set::MOI.AbstractVectorSet)
    m = length(model.added_constrained)
    for i in 1:MOI.dimension(set)
        push!(model.added_constrained, true)
    end
    ci = MOI.ConstraintIndex{MOI.VectorOfVariables, typeof(set)}(m + 1)
    return MOI.VariableIndex.(m .+ (1:MOI.dimension(set))), ci
end
function MOIU.allocate_constrained_variables(model::ConstrainedVariablesModel, set::MOI.AbstractVectorSet)
    return MOI.add_constrained_variables(model, set)
end
function MOIU.load_constrained_variables(model::ConstrainedVariablesModel, vis::Vector{MOI.VariableIndex}, ci::MOI.ConstraintIndex{MOI.VectorOfVariables}, set::MOI.AbstractVectorSet)
end
function MOI.add_constraint(model::ConstrainedVariablesModel,
                            func::MOI.VectorOfVariables,
                            set::MOI.AbstractVectorSet)
    return MOI.ConstraintIndex{typeof(func), typeof(set)}(func.variables[1].value)
end
function MOIU.allocate_constraint(model::ConstrainedVariablesModel, func::MOI.VectorOfVariables, set::MOI.AbstractVectorSet)
    return MOI.add_constraint(model, func, set)
end
function MOIU.load_constraint(model::ConstrainedVariablesModel, ::MOI.ConstraintIndex{MOI.VectorOfVariables}, ::MOI.VectorOfVariables, ::MOI.AbstractVectorSet)
end

@testset "Duplicates in VectorOfVariables: $allocate_load" for allocate_load in [false, true]
    src = MOIU.Model{Int}()
    x = MOI.add_variables(src, 3)
    cx = MOI.add_constraint(src, [x[1], x[3], x[1], x[2]], MOI.Nonnegatives(4))
    y, cy = MOI.add_constrained_variables(src, MOI.Nonpositives(3))
    dest = ConstrainedVariablesModel(allocate_load)
    idxmap = MOI.copy_to(dest, src)
    for vi in x
        @test !dest.added_constrained[idxmap[vi].value]
    end
    for vi in y
        @test dest.added_constrained[idxmap[vi].value]
    end
end

abstract type AbstractConstrainedVariablesModel <: MOI.ModelLike end
mutable struct OrderConstrainedVariablesModel <: AbstractConstrainedVariablesModel
    constraintIndices       ::Array{MOI.ConstraintIndex}
    inner                   ::MOIU.Model{Float64}
    OrderConstrainedVariablesModel() = new(MOI.ConstraintIndex[], MOIU.Model{Float64}())
end
mutable struct ReverseOrderConstrainedVariablesModel <: AbstractConstrainedVariablesModel
    constraintIndices       ::Array{MOI.ConstraintIndex}
    inner                   ::MOIU.Model{Float64}
    ReverseOrderConstrainedVariablesModel() = new(MOI.ConstraintIndex[], MOIU.Model{Float64}())
end

MOI.add_variables(model::AbstractConstrainedVariablesModel, n) = MOI.add_variables(model.inner, n)
MOI.add_variable(model::AbstractConstrainedVariablesModel) = MOI.add_variable(model.inner)

function MOI.add_constraint(model::AbstractConstrainedVariablesModel, f::F, s::S) where {F<:MOI.AbstractFunction, S<:MOI.AbstractSet}
    ci = MOI.add_constraint(model.inner, f, s)
    push!(model.constraintIndices, ci)
    return ci
end

function MOI.copy_to(dest::AbstractConstrainedVariablesModel, src::MOI.ModelLike; kws...)
    MOIU.automatic_copy_to(dest, src; kws...)
end

MOIU.supports_default_copy_to(model::AbstractConstrainedVariablesModel, ::Bool) = true

function MOI.empty!(model::AbstractConstrainedVariablesModel)
    model.constraintIndices = MOI.ConstraintIndex[]
    MOI.empty!(model.inner)
end

MOI.supports_constraint(::OrderConstrainedVariablesModel, ::Type{MOI.VectorOfVariables}, ::Type{MOI.Nonnegatives}) = false
MOI.supports_add_constrained_variables(::OrderConstrainedVariablesModel, ::Type{MOI.Nonnegatives}) = true
MOI.supports_constraint(::OrderConstrainedVariablesModel, ::Type{MOI.VectorOfVariables}, ::Type{MOI.Nonnegatives}) = true
MOI.supports_add_constrained_variables(::OrderConstrainedVariablesModel, ::Type{MOI.Nonpositives}) = false

MOI.supports_constraint(::ReverseOrderConstrainedVariablesModel, ::Type{MOI.VectorOfVariables}, ::Type{MOI.Nonnegatives}) = true
MOI.supports_add_constrained_variables(::ReverseOrderConstrainedVariablesModel, ::Type{MOI.Nonnegatives}) = false
MOI.supports_constraint(::ReverseOrderConstrainedVariablesModel, ::Type{MOI.VectorOfVariables}, ::Type{MOI.Nonnegatives}) = false
MOI.supports_add_constrained_variables(::ReverseOrderConstrainedVariablesModel, ::Type{MOI.Nonpositives}) = true

MOI.supports_constraint(::OrderConstrainedVariablesModel, ::Type{MOI.SingleVariable}, ::Type{<:MOI.GreaterThan}) = true
MOI.supports_add_constrained_variable(::OrderConstrainedVariablesModel, ::Type{<:MOI.GreaterThan}) = false
MOI.supports_constraint(::OrderConstrainedVariablesModel, ::Type{MOI.SingleVariable}, ::Type{<:MOI.LessThan}) = false
MOI.supports_add_constrained_variable(::OrderConstrainedVariablesModel, ::Type{<:MOI.LessThan}) = true

MOI.supports_constraint(::ReverseOrderConstrainedVariablesModel, ::Type{MOI.SingleVariable}, ::Type{<:MOI.GreaterThan}) = false
MOI.supports_add_constrained_variable(::ReverseOrderConstrainedVariablesModel, ::Type{<:MOI.GreaterThan}) = true
MOI.supports_constraint(::ReverseOrderConstrainedVariablesModel, ::Type{MOI.SingleVariable}, ::Type{<:MOI.LessThan}) = true
MOI.supports_add_constrained_variable(::ReverseOrderConstrainedVariablesModel, ::Type{<:MOI.LessThan}) = false

@testset "Create variables using supports_add_constrained_variable(s) (#987)" begin
    # With vectors
    src = MOIU.Model{Float64}()
    a, c1 = MOI.add_constrained_variables(src, MOI.Nonpositives(3))
    c2 = MOI.add_constraint(src, a, MOI.Nonnegatives(3))

    dest = OrderConstrainedVariablesModel()
    index_map = MOI.copy_to(dest, src)
    @test typeof(c1) == typeof(dest.constraintIndices[2])
    @test typeof(c2) == typeof(dest.constraintIndices[1])

    dest = ReverseOrderConstrainedVariablesModel()
    index_map = MOI.copy_to(dest, src)
    @test typeof(c1) == typeof(dest.constraintIndices[1])
    @test typeof(c2) == typeof(dest.constraintIndices[2])

    # With single variables
    src = MOIU.Model{Float64}()
    a, c1 = MOI.add_constrained_variable(src, MOI.GreaterThan{Float64}(5.0))
    c2 = MOI.add_constraint(src, a, MOI.LessThan{Float64}(1.0))

    dest = OrderConstrainedVariablesModel()
    index_map = MOI.copy_to(dest, src)
    @test typeof(c1) == typeof(dest.constraintIndices[2])
    @test typeof(c2) == typeof(dest.constraintIndices[1])

    dest = ReverseOrderConstrainedVariablesModel()
    index_map = MOI.copy_to(dest, src)
    @test typeof(c1) == typeof(dest.constraintIndices[1])
    @test typeof(c2) == typeof(dest.constraintIndices[2])
end

@testset "Filtering copy: check based on index" begin
    # Create a basic model.
    src = MOIU.Model{Float64}()
    x = MOI.add_variable(src)
    c1 = MOI.add_constraint(src, x, MOI.LessThan{Float64}(1.0))
    c2 = MOI.add_constraint(src, x, MOI.GreaterThan{Float64}(0.0))

    # Filtering function: the default case where this function always returns
    # true is already well-tested by the above cases.
    # Only keep the constraint c1.
    f = (cidx) -> cidx == c1

    # Perform the copy.
    dst = OrderConstrainedVariablesModel()
    index_map = MOI.copy_to(dst, src, filter_constraints=f)

    @test typeof(c1) == typeof(dst.constraintIndices[1])
    @test length(dst.constraintIndices) == 1
end

mutable struct BoundModel <: MOI.ModelLike
    # Type of model that only supports ≤ bound constraints. In particular, it does not support integrality constraints.
    inner::MOIU.Model{Float64}
    BoundModel() = new(MOIU.Model{Float64}())
end

MOI.add_variable(model::BoundModel) = MOI.add_variable(model.inner)
MOI.add_constraint(model::BoundModel, f::MOI.AbstractFunction, s::MOI.LessThan{Float64}) = MOI.add_constraint(model.inner, f, s)
MOI.supports_constraint(::BoundModel, ::Type{MOI.SingleVariable}, ::MOI.LessThan{Float64}) = true

MOIU.supports_default_copy_to(::BoundModel, ::Bool) = true
MOI.copy_to(dest::BoundModel, src::MOI.ModelLike; kws...) = MOIU.automatic_copy_to(dest, src; kws...)
MOI.empty!(model::BoundModel) = MOI.empty!(model.inner)

MOI.supports(::BoundModel, ::Type{MOI.NumberOfConstraints}) = true
MOI.get(model::BoundModel, attr::MOI.NumberOfConstraints) = MOI.get(model.inner, attr)

@testset "Filtering copy: check based on constraint type" begin
    # Create a basic model.
    src = MOIU.Model{Float64}()
    x = MOI.add_variable(src)
    c1 = MOI.add_constraint(src, x, MOI.LessThan{Float64}(10.0))
    c2 = MOI.add_constraint(src, x, MOI.Integer())

    # Filtering function: get rid of integrality constraint.
    f = (cidx) -> MOI.get(src, MOI.ConstraintSet(), cidx) != MOI.Integer()

    # Perform the unfiltered copy. This should throw an error (i.e. the implementation of BoundModel
    # should be correct).
    dst = BoundModel()
    @test_throws MOI.UnsupportedConstraint{MOI.SingleVariable, MOI.Integer} MOI.copy_to(dst, src)

    # Perform the filtered copy. This should not throw an error.
    dst = BoundModel()
    MOI.copy_to(dst, src, filter_constraints=f)
    @test MOI.get(dst, MOI.NumberOfConstraints{MOI.SingleVariable, MOI.LessThan{Float64}}()) == 1
    @test MOI.get(dst, MOI.NumberOfConstraints{MOI.SingleVariable, MOI.Integer}()) == 0
end


text:
if occursin(r"RvSpectMLEcoSystem$", pwd())   cd("NeidSolarScripts")   end
 if occursin(r"NeidSolarScripts$", pwd())   cd("test")   end
 using Pkg
 Pkg.activate(".")

using NeidSolarScripts
using CSV, DataFrames

HARPSN_df = CSV.read("../data/Sun_harpsn_qualflag.csv",DataFrame)
jd = HARPSN_df.JD .+ 2400000
HARPSN_df.Delta = HARPSN_df.FWHMobs.^2 .- HARPSN_df.FWHMsid.^2

# Test scalar
Delta = SolarRotation.CalculateFWHMDifference_SolarRotation_from_obs(jd[1], obs=:HARPSN)

# Test on array of times
Delta = SolarRotation.CalculateFWHMDifference_SolarRotation_from_obs.(jd[1:50:end],obs=:HARPSN)

using Plots
plot(jd,HARPSN_df.Delta, label="ACC")
plot!(jd[1:50:end],Delta, label="Julia")
xlabel!("Time (d)")
ylabel!("FWHM_obs^2 - FWHM_sid^2")

plot(jd[1:50:end],Delta.-HARPSN_df.Delta[1:50:end], label="Julia-ACC")
xlabel!("Time (d)")
ylabel!("Residuals in Δ(FWHM^2)")

#=
plot(jd[1:50:end],Delta./HARPSN_df.Delta[1:50:end], label="Julia/ACC")
xlabel!("Time (d)")
ylabel!("Ratio of Δ(FWHM^2)'s")
=#


text:
# https://github.com/jheinen/GR.jl

# significant contributions by @jheinen

supported_args(::GRBackend) = merge_with_base_supported([
    :annotations,
    :background_color_legend, :background_color_inside, :background_color_outside,
    :foreground_color_legend, :foreground_color_grid, :foreground_color_axis,
    :foreground_color_text, :foreground_color_border,
    :label,
    :seriescolor, :seriesalpha,
    :linecolor, :linestyle, :linewidth, :linealpha,
    :markershape, :markercolor, :markersize, :markeralpha,
    :markerstrokewidth, :markerstrokecolor, :markerstrokealpha,
    :fillrange, :fillcolor, :fillalpha,
    :bins,
    :layout,
    :title, :window_title,
    :guide, :lims, :ticks, :scale, :flip,
    :tickfont, :guidefont, :legendfont,
    :grid, :legend, :colorbar,
    :marker_z, :levels,
    :ribbon, :quiver,
    :orientation,
    :overwrite_figure,
    :polar,
    :aspect_ratio,
    :normalize, :weights,
    :inset_subplots,
])
supported_types(::GRBackend) = [
    :path, :scatter,
    :heatmap, :pie, :image,
    :contour, :path3d, :scatter3d, :surface, :wireframe,
    :shape
]
supported_styles(::GRBackend) = [:auto, :solid, :dash, :dot, :dashdot, :dashdotdot]
supported_markers(::GRBackend) = vcat(_allMarkers, Shape)
supported_scales(::GRBackend) = [:identity, :log10]
is_subplot_supported(::GRBackend) = true

function _initialize_backend(::GRBackend; kw...)
    @eval begin
        import GR
        export GR
    end
end

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

const gr_linetype = KW(
    :auto => 1,
    :solid => 1,
    :dash => 2,
    :dot => 3,
    :dashdot => 4,
    :dashdotdot => -1
)

const gr_markertype = KW(
    :auto => 1,
    :none => -1,
    :circle => -1,
    :rect => -7,
    :diamond => -13,
    :utriangle => -3,
    :dtriangle => -5,
    :pentagon => -21,
    :hexagon => -22,
    :heptagon => -23,
    :octagon => -24,
    :cross => 2,
    :xcross => 5,
    :star4 => -25,
    :star5 => -26,
    :star6 => -27,
    :star7 => -28,
    :star8 => -29,
    :vline => -30,
    :hline => -31
)

const gr_halign = KW(
    :left => 1,
    :hcenter => 2,
    :right => 3
)

const gr_valign = KW(
    :top => 1,
    :vcenter => 3,
    :bottom => 5
)

const gr_font_family = Dict(
    "times" => 1,
    "helvetica" => 5,
    "courier" => 9,
    "bookman" => 14,
    "newcenturyschlbk" => 18,
    "avantgarde" => 22,
    "palatino" => 26
)

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

function gr_getcolorind(v, a =nothing)
    c = getColor(v)
    idx = convert(Int, GR.inqcolorfromrgb(c.r, c.g, c.b))
    GR.settransparency(float(a==nothing ? alpha(c) : a))
    idx
end

gr_set_linecolor(c, a=nothing) = GR.setlinecolorind(gr_getcolorind(c, a))
gr_set_fillcolor(c, a=nothing) = GR.setfillcolorind(gr_getcolorind(c, a))
gr_set_markercolor(c, a=nothing) = GR.setmarkercolorind(gr_getcolorind(c, a))
gr_set_textcolor(c, a=nothing) = GR.settextcolorind(gr_getcolorind(c, a))

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

# draw line segments, splitting x/y into contiguous/finite segments
# note: this can be used for shapes by passing func `GR.fillarea`
function gr_polyline(x, y, func = GR.polyline)
    iend = 0
    n = length(x)
    while iend < n-1
        # set istart to the first index that is finite
        istart = -1
        for j = iend+1:n
            if isfinite(x[j]) && isfinite(y[j])
                istart = j
                break
            end
        end

        if istart > 0
            # iend is the last finite index
            iend = -1
            for j = istart+1:n
                if isfinite(x[j]) && isfinite(y[j])
                    iend = j
                else
                    break
                end
            end
        end

        # if we found a start and end, draw the line segment, otherwise we're done
        if istart > 0 && iend > 0
            func(x[istart:iend], y[istart:iend])
        else
            break
        end
    end
end

function gr_polaraxes(rmin, rmax)
    GR.savestate()
    GR.setlinetype(GR.LINETYPE_SOLID)
    GR.setlinecolorind(88)
    tick = 0.5 * GR.tick(rmin, rmax)
    n = round(Int, (rmax - rmin) / tick + 0.5)
    for i in 0:n
        r = float(i) / n
        if i % 2 == 0
            GR.setlinecolorind(88)
            if i > 0
                GR.drawarc(-r, r, -r, r, 0, 359)
            end
            GR.settextalign(GR.TEXT_HALIGN_LEFT, GR.TEXT_VALIGN_HALF)
            x, y = GR.wctondc(0.05, r)
            GR.text(x, y, string(signif(rmin + i * tick, 12)))
        else
            GR.setlinecolorind(90)
            GR.drawarc(-r, r, -r, r, 0, 359)
        end
    end
    for alpha in 0:45:315
        a = alpha + 90
        sinf = sin(a * pi / 180)
        cosf = cos(a * pi / 180)
        GR.polyline([sinf, 0], [cosf, 0])
        GR.settextalign(GR.TEXT_HALIGN_CENTER, GR.TEXT_VALIGN_HALF)
        x, y = GR.wctondc(1.1 * sinf, 1.1 * cosf)
        GR.textext(x, y, string(alpha, "^o"))
    end
    GR.restorestate()
end

# using the axis extrema and limit overrides, return the min/max value for this axis
gr_x_axislims(sp::Subplot) = axis_limits(sp[:xaxis])
gr_y_axislims(sp::Subplot) = axis_limits(sp[:yaxis])
gr_z_axislims(sp::Subplot) = axis_limits(sp[:zaxis])
gr_xy_axislims(sp::Subplot) = gr_x_axislims(sp)..., gr_y_axislims(sp)...

function gr_lims(axis::Axis, adjust::Bool, expand = nothing)
    if expand != nothing
        expand_extrema!(axis, expand)
    end
    lims = axis_limits(axis)
    if adjust
        GR.adjustrange(lims...)
    else
        lims
    end
end

function gr_fill_viewport(vp::AVec{Float64}, c)
    GR.savestate()
    GR.selntran(0)
    GR.setfillintstyle(GR.INTSTYLE_SOLID)
    gr_set_fillcolor(c)
    GR.fillrect(vp...)
    GR.selntran(1)
    GR.restorestate()
end

normalize_zvals(zv::Void) = zv
function normalize_zvals(zv::AVec)
    vmin, vmax = extrema(zv)
    if vmin == vmax
        zeros(length(zv))
    else
        (zv - vmin) ./ (vmax - vmin)
    end
end

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

# draw ONE Shape
function gr_draw_marker(xi, yi, msize, shape::Shape)
    sx, sy = shape_coords(shape)
    GR.selntran(0)
    xi, yi = GR.wctondc(xi, yi)
    GR.fillarea(xi + sx * 0.0015msize,
                yi + sy * 0.0015msize)
    GR.selntran(1)
end

# draw ONE symbol marker
function gr_draw_marker(xi, yi, msize::Number, shape::Symbol)
    GR.setmarkertype(gr_markertype[shape])
    GR.setmarkersize(0.3msize)
    GR.polymarker([xi], [yi])
end

# draw the markers, one at a time
function gr_draw_markers(d::KW, x, y, msize, mz)
    shape = d[:markershape]
    if shape != :none
        for i=1:length(x)
            msize = cycle(msize, i)
            cfunc = isa(shape, Shape) ? gr_set_fillcolor : gr_set_markercolor
            cfuncind = isa(shape, Shape) ? GR.setfillcolorind : GR.setmarkercolorind

            # draw a filled in shape, slightly bigger, to estimate a stroke
            cfunc(d[:markerstrokecolor], d[:markerstrokealpha])
            gr_draw_marker(x[i], y[i], msize*1.2, shape, )

            # draw the shape
            if mz == nothing
                cfunc(d[:markercolor], d[:markeralpha])
            else
                # pick a color from the pre-loaded gradient
                ci = round(Int, 1000 + cycle(mz, i) * 255)
                cfuncind(ci)
            end
            gr_draw_marker(x[i], y[i], msize, shape)
        end
    end
end

function gr_draw_markers(series::Series, x, y)
    isempty(x) && return
    d = series.d
    mz = normalize_zvals(d[:marker_z])
    GR.setfillintstyle(GR.INTSTYLE_SOLID)
    gr_draw_markers(d, x, y, d[:markersize], mz)
    if mz != nothing
        gr_colorbar(d[:subplot])
    end
end

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

function gr_set_line(w, style, c, a)
    GR.setlinetype(gr_linetype[style])
    GR.setlinewidth(w)
    gr_set_linecolor(c, a)
end

function gr_set_fill(c, a)
    gr_set_fillcolor(c, a)
    GR.setfillintstyle(GR.INTSTYLE_SOLID)
end

# this stores the conversion from a font pointsize to "percentage of window height" (which is what GR uses)
const _gr_point_mult = zeros(1)

# set the font attributes... assumes _gr_point_mult has been populated already
function gr_set_font(f::Font)
    family = lowercase(f.family)
    GR.setcharheight(_gr_point_mult[1] * f.pointsize)
    GR.setcharup(sin(f.rotation), cos(f.rotation))
    if haskey(gr_font_family, family)
        GR.settextfontprec(100 + gr_font_family[family], GR.TEXT_PRECISION_STRING)
    end
    gr_set_textcolor(f.color)
    GR.settextalign(gr_halign[f.halign], gr_valign[f.valign])
end

# --------------------------------------------------------------------------------------
# viewport plot area

# this stays constant for a given subplot while displaying that subplot.
# values are [xmin, xmax, ymin, ymax].  they range [0,1].
const viewport_plotarea = zeros(4)

function gr_viewport_from_bbox(bb::BoundingBox, w, h, viewport_canvas)
    viewport = zeros(4)
    viewport[1] = viewport_canvas[2] * (left(bb) / w)
    viewport[2] = viewport_canvas[2] * (right(bb) / w)
    viewport[3] = viewport_canvas[4] * (1.0 - bottom(bb) / h)
    viewport[4] = viewport_canvas[4] * (1.0 - top(bb) / h)
    viewport
end

# change so we're focused on the viewport area
function gr_set_viewport_cmap(sp::Subplot)
    GR.setviewport(
        viewport_plotarea[2] + (is3d(sp) ? 0.04 : 0.02),
        viewport_plotarea[2] + (is3d(sp) ? 0.07 : 0.05),
        viewport_plotarea[3],
        viewport_plotarea[4]
    )
end

# reset the viewport to the plot area
function gr_set_viewport_plotarea()
    GR.setviewport(
        viewport_plotarea[1],
        viewport_plotarea[2],
        viewport_plotarea[3],
        viewport_plotarea[4]
    )
end

function gr_set_viewport_polar()
    xmin, xmax, ymin, ymax = viewport_plotarea
    ymax -= 0.05 * (xmax - xmin)
    xcenter = 0.5 * (xmin + xmax)
    ycenter = 0.5 * (ymin + ymax)
    r = 0.5 * min(xmax - xmin, ymax - ymin)
    GR.setviewport(xcenter -r, xcenter + r, ycenter - r, ycenter + r)
    GR.setwindow(-1, 1, -1, 1)
    r
end

# add the colorbar
function gr_colorbar(sp::Subplot)
    if sp[:colorbar] != :none
        gr_set_viewport_cmap(sp)
        GR.colormap()
        gr_set_viewport_plotarea()
    end
end

gr_view_xcenter() = 0.5 * (viewport_plotarea[1] + viewport_plotarea[2])
gr_view_ycenter() = 0.5 * (viewport_plotarea[3] + viewport_plotarea[4])
gr_view_xdiff() = viewport_plotarea[2] - viewport_plotarea[1]
gr_view_ydiff() = viewport_plotarea[4] - viewport_plotarea[3]

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

function gr_set_gradient(c, a)
    grad = isa(c, ColorGradient) ? c : default_gradient()
    grad = ColorGradient(grad, alpha=a)
    for (i,z) in enumerate(linspace(0, 1, 256))
        c = getColorZ(grad, z)
        GR.setcolorrep(999+i, red(c), green(c), blue(c))
    end
    grad
end

# this is our new display func... set up the viewport_canvas, compute bounding boxes, and display each subplot
function gr_display(plt::Plot)
    GR.clearws()

    # collect some monitor/display sizes in meters and pixels
    display_width_meters, display_height_meters, display_width_px, display_height_px = GR.inqdspsize()
    display_width_ratio = display_width_meters / display_width_px
    display_height_ratio = display_height_meters / display_height_px

    # compute the viewport_canvas, normalized to the larger dimension
    viewport_canvas = Float64[0,1,0,1]
    w, h = plt[:size]
    if w > h
        ratio = float(h) / w
        msize = display_width_ratio * w
        GR.setwsviewport(0, msize, 0, msize * ratio)
        GR.setwswindow(0, 1, 0, ratio)
        viewport_canvas[3] *= ratio
        viewport_canvas[4] *= ratio
    else
        ratio = float(w) / h
        msize = display_height_ratio * h
        GR.setwsviewport(0, msize * ratio, 0, msize)
        GR.setwswindow(0, ratio, 0, 1)
        viewport_canvas[1] *= ratio
        viewport_canvas[2] *= ratio
    end

    # fill in the viewport_canvas background
    gr_fill_viewport(viewport_canvas, plt[:background_color_outside])

    # update point mult
    px_per_pt = px / pt
    _gr_point_mult[1] = px_per_pt / h

    # subplots:
    for sp in plt.subplots
        gr_display(sp, w*px, h*px, viewport_canvas)
    end

    GR.updatews()
end

function gr_display(sp::Subplot{GRBackend}, w, h, viewport_canvas)
    # the viewports for this subplot
    viewport_subplot = gr_viewport_from_bbox(bbox(sp), w, h, viewport_canvas)
    viewport_plotarea[:] = gr_viewport_from_bbox(plotarea(sp), w, h, viewport_canvas)

    # fill in the plot area background
    bg = getColor(sp[:background_color_inside])
    gr_fill_viewport(viewport_plotarea, bg)

    # reduced from before... set some flags based on the series in this subplot
    # TODO: can these be generic flags?
    outside_ticks = false
    cmap = false
    draw_axes = true
    # axes_2d = true
    for series in series_list(sp)
        st = series.d[:seriestype]
        if st in (:contour, :surface, :heatmap) || series.d[:marker_z] != nothing
            cmap = true
        end
        if st == :pie
            draw_axes = false
        end
        if st == :heatmap
            outside_ticks = true
        end
    end

    if cmap && sp[:colorbar] != :none
        # note: add extra midpadding on the right for the colorbar
        viewport_plotarea[2] -= 0.1
    end

    # set our plot area view
    gr_set_viewport_plotarea()

    # these are the Axis objects, which hold scale, lims, etc
    xaxis = sp[:xaxis]
    yaxis = sp[:yaxis]
    zaxis = sp[:zaxis]

    # get data limits and set the scale flags and window
    data_lims = gr_xy_axislims(sp)
    xmin, xmax, ymin, ymax = data_lims
    scale = 0
    xtick, ytick = 1, 1
    if xmax > xmin && ymax > ymin
        # NOTE: for log axes, the major_x and major_y - if non-zero (omit labels) - control the minor grid lines (1 = draw 9 minor grid lines, 2 = no minor grid lines)
        # NOTE: for log axes, the x_tick and y_tick - if non-zero (omit axes) - only affect the output appearance (1 = nomal, 2 = scientiic notation)
        xaxis[:scale] == :log10 && (scale |= GR.OPTION_X_LOG)
        yaxis[:scale] == :log10 && (scale |= GR.OPTION_Y_LOG)
        xaxis[:flip]            && (scale |= GR.OPTION_FLIP_X)
        yaxis[:flip]            && (scale |= GR.OPTION_FLIP_Y)
        if scale & GR.OPTION_X_LOG == 0
            majorx = 1 #5
            xtick = GR.tick(xmin, xmax) / majorx
        else
            # log axis
            xtick = 2  # scientific notation
            majorx = 2 # no minor grid lines
        end
        if scale & GR.OPTION_Y_LOG == 0
            majory = 1 #5
            ytick = GR.tick(ymin, ymax) / majory
        else
            # log axis
            ytick = 2  # scientific notation
            majory = 2 # no minor grid lines
        end

        # NOTE: setwindow sets the "data coordinate" limits of the current "viewport"
        GR.setwindow(xmin, xmax, ymin, ymax)
        GR.setscale(scale)
    end

    # draw the axes
    gr_set_font(xaxis[:tickfont])
    gr_set_textcolor(xaxis[:foreground_color_text])
    GR.setlinewidth(1)

    if is3d(sp)
        zmin, zmax = gr_lims(zaxis, true)
        GR.setspace(zmin, zmax, 40, 70)
        xtick = GR.tick(xmin, xmax) / 2
        ytick = GR.tick(ymin, ymax) / 2
        ztick = GR.tick(zmin, zmax) / 2
        ticksize = 0.01 * (viewport_plotarea[2] - viewport_plotarea[1])

        # GR.setlinetype(GR.LINETYPE_DOTTED)
        if sp[:grid]
            GR.grid3d(xtick, 0, ztick, xmin, ymin, zmin, 2, 0, 2)
            GR.grid3d(0, ytick, 0, xmax, ymin, zmin, 0, 2, 0)
        end
        GR.axes3d(xtick, 0, ztick, xmin, ymin, zmin, 2, 0, 2, -ticksize)
        GR.axes3d(0, ytick, 0, xmax, ymin, zmin, 0, 2, 0, ticksize)

    elseif ispolar(sp)
        r = gr_set_viewport_polar()
        rmin, rmax = GR.adjustrange(minimum(r), maximum(r))
        gr_polaraxes(rmin, rmax)

    elseif draw_axes
        # draw the grid lines
        # TODO: control line style/width
        # GR.setlinetype(GR.LINETYPE_DOTTED)
        if sp[:grid]
            gr_set_linecolor(sp[:foreground_color_grid])
            GR.grid(xtick, ytick, 0, 0, majorx, majory)
        end

        window_diag = sqrt(gr_view_xdiff()^2 + gr_view_ydiff()^2)
        ticksize = 0.0075 * window_diag
        if outside_ticks
            ticksize = -ticksize
        end
        # TODO: this should be done for each axis separately
        gr_set_linecolor(xaxis[:foreground_color_axis])

        x1, x2 = xaxis[:flip] ? (xmax,xmin) : (xmin,xmax)
        y1, y2 = yaxis[:flip] ? (ymax,ymin) : (ymin,ymax)
        GR.axes(xtick, ytick, x1, y1, 1, 1, ticksize)
        GR.axes(xtick, ytick, x2, y2, -1, -1, -ticksize)
    end
    # end

    # add the guides
    GR.savestate()
    if sp[:title] != ""
        gr_set_font(sp[:titlefont])
        GR.settextalign(GR.TEXT_HALIGN_CENTER, GR.TEXT_VALIGN_TOP)
        gr_set_textcolor(sp[:foreground_color_title])
        GR.text(gr_view_xcenter(), viewport_subplot[4], sp[:title])
    end

    if xaxis[:guide] != ""
        gr_set_font(xaxis[:guidefont])
        GR.settextalign(GR.TEXT_HALIGN_CENTER, GR.TEXT_VALIGN_BOTTOM)
        gr_set_textcolor(xaxis[:foreground_color_guide])
        GR.text(gr_view_xcenter(), viewport_subplot[3], xaxis[:guide])
    end

    if yaxis[:guide] != ""
        gr_set_font(yaxis[:guidefont])
        GR.settextalign(GR.TEXT_HALIGN_CENTER, GR.TEXT_VALIGN_TOP)
        GR.setcharup(-1, 0)
        gr_set_textcolor(yaxis[:foreground_color_guide])
        GR.text(viewport_subplot[1], gr_view_ycenter(), yaxis[:guide])
    end
    GR.restorestate()

    # TODO: can we remove?
    gr_set_font(xaxis[:tickfont])
    GR.setcolormap(1000 + GR.COLORMAP_COOLWARM)

    for (idx, series) in enumerate(series_list(sp))
        d = series.d
        st = d[:seriestype]

        # update the current stored gradient
        if st in (:contour, :surface, :wireframe, :heatmap)
            gr_set_gradient(d[:fillcolor], d[:fillalpha])
        elseif d[:marker_z] != nothing
            d[:markercolor] = gr_set_gradient(d[:markercolor], d[:markeralpha])
        end

        GR.savestate()

        # update the bounding window
        if ispolar(sp)
            gr_set_viewport_polar()
        else
            xmin, xmax, ymin, ymax = data_lims
            if xmax > xmin && ymax > ymin
                GR.setwindow(xmin, xmax, ymin, ymax)
            end
        end

        x, y, z = d[:x], d[:y], d[:z]
        frng = d[:fillrange]

        # recompute data
        if st in (:contour, :surface, :wireframe)
            z = vec(transpose_z(d, z.surf, false))
        elseif ispolar(sp)
            if frng != nothing
                _, frng = convert_to_polar(x, frng, (rmin, rmax))
            end
            x, y = convert_to_polar(x, y, (rmin, rmax))
        end

        if st in (:path, :scatter)
            if length(x) > 1

                # do area fill
                if frng != nothing
                    gr_set_fillcolor(d[:fillcolor], d[:fillalpha])
                    GR.setfillintstyle(GR.INTSTYLE_SOLID)
                    frng = isa(frng, Number) ? Float64[frng] : frng
                    nx, ny, nf = length(x), length(y), length(frng)
                    n = max(nx, ny)
                    fx, fy = zeros(2n), zeros(2n)
                    for i=1:n
                        fx[i] = fx[end-i+1] = cycle(x,i)
                        fy[i] = cycle(y,i)
                        fy[end-i+1] = cycle(frng,i)
                    end
                    GR.fillarea(fx, fy)
                end

                # draw the line(s)
                if st == :path
                    gr_set_line(d[:linewidth], d[:linestyle], d[:linecolor], d[:linealpha])
                    gr_polyline(x, y)
                end
            end

            if d[:markershape] != :none
                gr_draw_markers(series, x, y)
            end

        elseif st == :contour
            zmin, zmax = gr_lims(zaxis, false)
            if typeof(d[:levels]) <: Array
                h = d[:levels]
            else
                h = linspace(zmin, zmax, d[:levels])
            end
            GR.setspace(zmin, zmax, 0, 90)
            if d[:fillrange] != nothing
                GR.surface(x, y, z, GR.OPTION_CELL_ARRAY)
            else
                GR.contour(x, y, h, z, 1000)
            end

            # create the colorbar of contour levels
            if sp[:colorbar] != :none
                gr_set_viewport_cmap(sp)
                l = round(Int32, 1000 + (h - minimum(h)) / (maximum(h) - minimum(h)) * 255)
                GR.setwindow(xmin, xmax, zmin, zmax)
                GR.cellarray(xmin, xmax, zmax, zmin, 1, length(l), l)
                ztick = 0.5 * GR.tick(zmin, zmax)
                GR.axes(0, ztick, xmax, zmin, 0, 1, 0.005)
                gr_set_viewport_plotarea()
            end

        elseif st in [:surface, :wireframe]
            if st == :surface
                GR.gr3.surface(x, y, z, GR.OPTION_COLORED_MESH)
            else
                GR.setfillcolorind(0)
                GR.surface(x, y, z, GR.OPTION_FILLED_MESH)
            end
            cmap && gr_colorbar(sp)

        elseif st == :heatmap
            z = vec(transpose_z(d, z.surf, false))
            zmin, zmax = gr_lims(zaxis, true)
            GR.setspace(zmin, zmax, 0, 90)
            GR.surface(x, y, z, GR.OPTION_COLORED_MESH)
            cmap && gr_colorbar(sp)

        elseif st in (:path3d, :scatter3d)
            # draw path
            if st == :path3d
                if length(x) > 1
                    gr_set_line(d[:linewidth], d[:linestyle], d[:linecolor], d[:linealpha])
                    GR.polyline3d(x, y, z)
                end
            end

            # draw markers
            if st == :scatter3d || d[:markershape] != :none
                x2, y2 = unzip(map(GR.wc3towc, x, y, z))
                gr_draw_markers(series, x2, y2)
            end

        # TODO: replace with pie recipe
        elseif st == :pie
            GR.selntran(0)
            GR.setfillintstyle(GR.INTSTYLE_SOLID)
            xmin, xmax, ymin, ymax = viewport_plotarea
            ymax -= 0.05 * (xmax - xmin)
            xcenter = 0.5 * (xmin + xmax)
            ycenter = 0.5 * (ymin + ymax)
            if xmax - xmin > ymax - ymin
                r = 0.5 * (ymax - ymin)
                xmin, xmax = xcenter - r, xcenter + r
            else
                r = 0.5 * (xmax - xmin)
                ymin, ymax = ycenter - r, ycenter + r
            end
            labels = pie_labels(sp, series)
            slices = d[:y]
            numslices = length(slices)
            total = sum(slices)
            a1 = 0
            x = zeros(3)
            y = zeros(3)
            for i in 1:numslices
                a2 = round(Int, a1 + (slices[i] / total) * 360.0)
                GR.setfillcolorind(980 + (i-1) % 20)
                GR.fillarc(xmin, xmax, ymin, ymax, a1, a2)
                alpha = 0.5 * (a1 + a2)
                cosf = r * cos(alpha * pi / 180)
                sinf = r * sin(alpha * pi / 180)
                x[1] = xcenter + cosf
                y[1] = ycenter + sinf
                x[2] = x[1] + 0.1 * cosf
                y[2] = y[1] + 0.1 * sinf
                y[3] = y[2]
                if 90 <= alpha < 270
                    x[3] = x[2] - 0.05
                    GR.settextalign(GR.TEXT_HALIGN_RIGHT, GR.TEXT_VALIGN_HALF)
                    GR.text(x[3] - 0.01, y[3], string(labels[i]))
                else
                    x[3] = x[2] + 0.05
                    GR.settextalign(GR.TEXT_HALIGN_LEFT, GR.TEXT_VALIGN_HALF)
                    GR.text(x[3] + 0.01, y[3], string(labels[i]))
                end
                gr_polyline(x, y)
                a1 = a2
            end
            GR.selntran(1)

        elseif st == :shape
            # draw the shapes
            gr_set_line(d[:markerstrokewidth], :solid, d[:markerstrokecolor], d[:markerstrokealpha])
            gr_polyline(d[:x], d[:y])

            # draw the interior
            gr_set_fill(d[:markercolor], d[:markeralpha])
            gr_polyline(d[:x], d[:y], GR.fillarea)

        elseif st == :image
            img = d[:z].surf
            w, h = size(img)
            if eltype(img) <: Colors.AbstractGray
                grey = round(UInt8, float(img) * 255)
                rgba = map(c -> UInt32( 0xff000000 + Int(c)<<16 + Int(c)<<8 + Int(c) ), grey)
            else
                rgba = map(c -> UInt32( round(Int, alpha(c) * 255) << 24 +
                                        round(Int,  blue(c) * 255) << 16 +
                                        round(Int, green(c) * 255) << 8  +
                                        round(Int,   red(c) * 255) ), img)
            end
            GR.drawimage(xmin, xmax, ymax, ymin, w, h, rgba)
        end

        GR.restorestate()
    end

    # add the legend
    if sp[:legend] != :none
        GR.savestate()
        GR.selntran(0)
        GR.setscale(0)
        gr_set_font(sp[:legendfont])
        w = 0
        i = 0
        n = 0
        for series in series_list(sp)
            should_add_to_legend(series) || continue
            n += 1
            if typeof(series.d[:label]) <: Array
                i += 1
                lab = series.d[:label][i]
            else
                lab = series.d[:label]
            end
            tbx, tby = GR.inqtext(0, 0, lab)
            w = max(w, tbx[3] - tbx[1])
        end
        if w > 0
            xpos = viewport_plotarea[2] - 0.05 - w
            ypos = viewport_plotarea[4] - 0.06
            dy = _gr_point_mult[1] * sp[:legendfont].pointsize * 1.75
            GR.setfillintstyle(GR.INTSTYLE_SOLID)
            gr_set_fillcolor(sp[:background_color_legend])
            GR.fillrect(xpos - 0.08, xpos + w + 0.02, ypos + dy, ypos - dy * n)
            GR.setlinetype(1)
            GR.setlinewidth(1)
            GR.drawrect(xpos - 0.08, xpos + w + 0.02, ypos + dy, ypos - dy * n)
            i = 0
            for series in series_list(sp)
                should_add_to_legend(series) || continue
                d = series.d
                st = d[:seriestype]
                GR.setlinewidth(d[:linewidth])
                if st == :path
                    gr_set_linecolor(d[:linecolor], d[:linealpha])
                    GR.setlinetype(gr_linetype[d[:linestyle]])
                    GR.polyline([xpos - 0.07, xpos - 0.01], [ypos, ypos])
                end

                gr_draw_markers(d, xpos-[0.06,0.02], [ypos,ypos], 10, nothing)
                # shape = d[:markershape]
                # if shape != :none #st == :scatter || d[:markershape] != :none
                #     msize = 10
                #     for xoff in [0.06,0.02]
                #         gr_set_markercolor(d[:markerstrokecolor], d[:markerstrokealpha])
                #         gr_draw_marker(xpos-xoff, ypos, msize*1.1, shape)
                #         gr_set_markercolor(d[:markercolor], d[:markeralpha])
                #         gr_draw_marker(xpos-xoff, ypos, msize, shape)
                #     end
                #     # gr_setmarkershape(d)
                #     # if st == :path
                #     #     gr_polymarker(d, [xpos - 0.06, xpos - 0.02], [ypos, ypos])
                #     # else
                #     #     gr_polymarker(d, [xpos - 0.06, xpos - 0.04, xpos - 0.02], [ypos, ypos, ypos])
                #     # end
                # end
                if typeof(d[:label]) <: Array
                    i += 1
                    lab = d[:label][i]
                else
                    lab = d[:label]
                end
                GR.settextalign(GR.TEXT_HALIGN_LEFT, GR.TEXT_VALIGN_HALF)
                gr_set_textcolor(sp[:foreground_color_legend])
                GR.text(xpos, ypos, lab)
                ypos -= dy
            end
        end
        GR.selntran(1)
        GR.restorestate()
    end

    # add annotations
    GR.savestate()
    for ann in sp[:annotations]
        x, y, val = ann
        x, y = GR.wctondc(x, y)
        gr_set_font(val.font)
        GR.text(x, y, val.str)
    end
    GR.restorestate()
end

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

const _gr_mimeformats = Dict(
    "application/pdf"         => "pdf",
    "image/png"               => "png",
    "application/postscript"  => "ps",
    "image/svg+xml"           => "svg",
)

for (mime, fmt) in _gr_mimeformats
    @eval function _writemime(io::IO, ::MIME{Symbol($mime)}, plt::Plot{GRBackend})
        GR.emergencyclosegks()
        wstype = haskey(ENV, "GKS_WSTYPE") ? ENV["GKS_WSTYPE"] : "0"
        filepath = tempname() * "." * $fmt
        ENV["GKS_WSTYPE"] = $fmt
        ENV["GKS_FILEPATH"] = filepath
        gr_display(plt)
        GR.emergencyclosegks()
        write(io, readall(filepath))
        ENV["GKS_WSTYPE"] = wstype
        rm(filepath)
    end
end

function _display(plt::Plot{GRBackend})
    gr_display(plt)
end


text:
fails = 0
@context "bar is bar"
f = jtest(1==1,
          2==2,
          3==3)
fails += f

@context "bar is not qux"
f = jtest( 1==1,
           2==2,
           3==3)
fails += f


text:
__precompile__()

#=
RequestsCache
A cache mechanism for [Request.jl](https://github.com/JuliaWeb/Requests.jl)

Based on JLD https://github.com/JuliaLang/JLD.jl

Inspired by [requests-cache](http://requests-cache.readthedocs.org/).

=#

module RequestsCache

    export get, CachedSession, Session

    #import Dates
    import Base: read
    import Requests: do_request, do_stream_request, ResponseStream
    import URIParser: URI
    import HttpCommon: Response
    import JLD: jldopen, write

    struct PreparedQuery
        verb::String
        uri::URI
        args::Array{Any,1}
    end
    Base.hash(q::PreparedQuery, h::UInt) = hash(string(q.verb), hash(string(q.uri), hash(q.args)))

    function create_query(verb::String, uri::URI; args...)
        return PreparedQuery(verb, uri, args)
    end
    create_query(verb::String, uri::String; args...) = create_query(verb, URI(uri); args...)

    struct CachedSessionType
        cache_name::String
        backend::String
        expire_after
    end

    function CachedSession(; cache_name="cache.jld", backend="jld", expire_after=Base.Dates.Day(1))
        if backend == "jld"
            CachedSessionType(cache_name, backend, expire_after)
        else
            error("'$(backend)' is not a supported backend")
        end
    end

    function Session()
        CachedSessionType("", "", Base.Dates.Day(0))
    end

    struct CachedResponse
        dt_stored::DateTime
        response::Response
    end

    function UTCnow()
        Dates.now(Dates.UTC)
    end

    #=
    struct CachedResponseStream{T}
        dt_stored::DateTime
        response::ResponseStream{T}
    end

    function write(session::CachedSessionType, prepared_query::PreparedQuery, response::ResponseStream{TCPSocket})
        backend = lowercase(session.backend)
        filename = session.cache_name
        key = string(hash(prepared_query))
        cached_response = CachedResponseStream(UTCnow(), response)
        if backend == "jld"
            jldopen(filename, "w") do file
                println("Write $cached_response with key='$key' to '$filename'")
                write(file, key, cached_response)
            end
        else
            error("'$(backend)' is not a supported backend for writing")
        end
    end
    =#

    function write(session::CachedSessionType, prepared_query::PreparedQuery, response::Response)
        backend = lowercase(session.backend)
        filename = session.cache_name
        key = string(hash(prepared_query))
        cached_response = CachedResponse(UTCnow(), response)
        if backend == "jld"
            jldopen(filename, "w") do file
                println("Write $cached_response with key='$key' to '$filename'")
                write(file, key, cached_response)
            end
        else
            error("'$(backend)' is not a supported backend for writing")
        end
    end

    function read(session::CachedSessionType, prepared_query::PreparedQuery)
        backend = lowercase(session.backend)
        filename = session.cache_name
        key = string(hash(prepared_query))
        if backend == "jld"
            retrieved_response = jldopen(filename, "r") do file
                println("Read key='$key' from '$filename'")
                read(file, key)
            end
            return retrieved_response
        else
            error("'$(backend)' is not a supported backend for reading")
        end
    end

    function execute_remote(prepared_query::PreparedQuery)
        println("execute_remote $(prepared_query.verb) $(prepared_query.uri) $(prepared_query.args)")
        #prepared_query.verb(string(prepared_query.uri); prepared_query.args...)
        verb = uppercase(string(prepared_query.verb))
        if !contains(verb, "_STREAMING")
            do_request(prepared_query.uri, verb; prepared_query.args...)
        else
            do_stream_request(prepared_query.uri, verb; prepared_query.args...)
        end
    end

    function execute_local(session::CachedSessionType, prepared_query::PreparedQuery)
        println("execute_local")
        try
            retrieved_response = read(session, prepared_query)
            dt_expiration = retrieved_response.dt_stored + session.expire_after
            if dt_expiration > UTCnow()
                println("Not expired")
                return retrieved_response.response
            else
                println("Cache expired - update is necessary")
                response = execute_remote(prepared_query)
                write(session, prepared_query, response)
            end
        catch LoadError
            println("LoadError $session")
            response = execute_remote(prepared_query)
            write(session, prepared_query, response)
            println("Write to $session")
            return response
        end
    end

    function execute(prepared_query::PreparedQuery; session=Session())
        println(session)
        if session.backend == ""
            execute_remote(prepared_query)
        else
            execute_local(session, prepared_query)
        end
    end

    for f in [:get, :post, :put, :delete, :head,
              :trace, :options, :patch, :connect]
        f_str = uppercase(string(f))
        #f_stream = symbol(string(f, "_streaming"))
        @eval begin
            function ($f)(session::CachedSessionType, uri::URI, data::String; headers::Dict=Dict())
                #do_request(uri, $f_str; data=data, headers=headers)
                prepared_query = create_query($f_str, uri; data=data, headers=headers)
                response = execute(prepared_query; session=session)
            end

            ($f)(session::CachedSessionType, uri::String; args...) = ($f)(session, URI(uri); args...)
            function ($f)(session::CachedSessionType, uri::URI; args...)
                #do_request(uri, $f_str; args...)
                prepared_query = create_query($f_str, uri; args...)
                response = execute(prepared_query; session=session)
            end

            #function ($f_stream)(uri::URI, data::String; headers::Dict=Dict())
            #    do_stream_request(uri, $f_str; data=data, headers=headers)
            #end
            #($f_stream)(uri::String; args...) = ($f_stream)(URI(uri); args...)
            #($f_stream)(uri::URI; args...) = do_stream_request(uri, $f_str; args...)
        end
    end

end


text:
## lsl benchmark

using MultipleTesting

p = rand(1000);
m = 1000;

@time for i in 1:m
    pi0 = MultipleTesting.lsl_pi0(p)
end

@time for i in 1:m
    pi0 = MultipleTesting.lsl_pi0_vec(p)
end


text: print("Hello world")

max_stars_repo_path: 1_introduction/Hello world.jl

max_stars_repo_name: F-Delpech/evolution

max_issues_repo_path: 1_introduction/Hello world.jl

max_issues_repo_name: F-Delpech/evolution

max_forks_repo_path: 1_introduction/Hello world.jl

max_forks_repo_name: F-Delpech/evolution


text:
# # Approximate Inference with NUTS-HMC

# Loading the necessary packages and setting seed.

using AbstractGPs, Plots, Random
Random.seed!(1234);

# Load toy regression
# [dataset](https://github.com/GPflow/docs/blob/master/doc/source/notebooks/basics/data/regression_1D.csv)
# taken from GPFlow examples.

x = [0.8658165855998895, 0.6661700880180962, 0.8049218148148531, 0.7714303440386239,
    0.14790478354654835, 0.8666105548197428, 0.007044577166530286, 0.026331737288148638,
    0.17188596617099916, 0.8897812990554013, 0.24323574561119998, 0.028590102134105955];
y = [1.5255314337144372, 3.6434202968230003, 3.010885733911661, 3.774442382979625,
    3.3687639483798324, 1.5506452040608503, 3.790447985799683, 3.8689707574953,
    3.4933565751758713, 1.4284538820635841, 3.8715350915692364, 3.7045949061144983];
scatter(x, y; xlabel="x", ylabel="y")

# Split the observations into train and test set.

(x_train, y_train) = (x[begin:8], y[begin:8]);
(x_test, y_test) = (x[9:end], y[9:end]);

# Instantiate the kernel.

k = Matern52Kernel()

# Instantiate a Gaussian Process with the given kernel `k`.

f = GP(k)

# Instantiate a `FiniteGP`, a finite dimentional projection at the inputs of the dataset
# observed under Gaussian Noise with $\sigma = 0.001$ .

fx = f(x_train, 0.001)

# Data's log-likelihood w.r.t prior `GP`.

logpdf(fx, y_train)

# Calculating the exact posterior over `f` given `y`. The GP's kernel currently has
# some arbitrary fixed parameters.

p_fx = posterior(fx, y_train)

# Data's log-likelihood under the posterior `GP`. We see that it drastically increases.

logpdf(p_fx(x_test), y_test)

# Plot the posterior `p_fx` along with the observations.

plt = scatter(x_train, y_train, label="Train Data")
scatter!(plt, x_test, y_test, label="Test Data")
plot!(plt, p_fx, 0:0.001:1; label="Posterior")

# # Hamiltonian Monte Carlo Sampler

# Previously we computed the log likelihood of the untuned kernel parameters of the GP,
# $-1.285$. We now also perform approximate inference over said kernel parameters using
# the No-U-Turn Sampler - Hamiltonian Monte Carlo(NUTS-HMC) Sampler
# ([paper](http://www.jmlr.org/papers/volume15/hoffman14a/hoffman14a.pdf))
# provided by [AdvancedHMC.jl](https://github.com/TuringLang/AdvancedHMC.jl/).
# We start of by loading necessary packages.

using AdvancedHMC, Distributions, ForwardDiff

# We define a function which returns log-probability of the data under the
# GP / log-likelihood of the parameters of the GP.

function logp(params; x=x_train, y=y_train)
    kernel = ScaledKernel(
        transform(
            Matern52Kernel(),
            ScaleTransform(exp(params[1]))
        ),
        exp(params[2])
    )
    f = GP(kernel)
    fx = f(x, 0.01)
    return logpdf(fx, y)
end

# Set the number of samples to draw and warmup iterations

n_samples, n_adapts = 2_000, 1_000
initial_params = rand(2)

# Define a Hamiltonian system

metric = DiagEuclideanMetric(2)
hamiltonian = Hamiltonian(metric, logp, ForwardDiff)

# Define a leapfrog solver, with initial step size chosen heuristically

initial_ϵ = find_good_stepsize(hamiltonian, initial_params)
integrator = Leapfrog(initial_ϵ)

# Define an HMC sampler, with the following components
#   - multinomial sampling scheme,
#   - generalised No-U-Turn criteria, and
#   - windowed adaption for step-size and diagonal mass matrix

proposal = NUTS{MultinomialTS, GeneralisedNoUTurn}(integrator)
adaptor = StanHMCAdaptor(MassMatrixAdaptor(metric), StepSizeAdaptor(0.8, integrator))

# Run the sampler to draw samples from the specified Gaussian, where
#   - `samples` will store the samples
#   - `stats` will store diagnostic statistics for each sample

samples, stats = sample(
    hamiltonian,
    proposal,
    initial_params,
    n_samples,
    adaptor,
    n_adapts;
    progress=false
)
samples_mat = reduce(hcat, samples)';

# Mean of samples of both the parameters.

mean_params = mean(samples_mat; dims=1)

# Plotting a histogram of the samples for the two parameters.
# The vertical line in each graph indicates the mean of the samples.

plt = histogram(samples_mat; layout=2, labels= "Param")
vline!(plt, mean_params; layout=2, label="Mean")

# Average log-marginal-probability of data with posterior kernel parameter samples
# sampled using ESS. We can observe that there is significant improvement over exact
# posterior with default kernel parameters.

mean(logp(param; x=x_test, y=y_test) for param in samples)

# Plotting sampled functions from posterior with tuned parameters

plt = scatter(x_train, y_train, label="Train Data")
scatter!(plt, x_test, y_test, label="Test Data")
for params in @view(samples[(end-100):end,:])
    opt_kernel = ScaledKernel(
        transform(
            Matern52Kernel(),
            ScaleTransform(exp(params[1]))
        ),
        exp(params[2])
    )
    f = GP(opt_kernel)
    p_fx = posterior(f(x_train, 0.1), y_train)
    sampleplot!(plt, p_fx(collect(0:0.02:1)), 1)
end
plt


text:
# Modified from Shreyas Kowshik's implementation.

# NOTE: TRPO will not work with categorical policies as nested AD is not currently defined for `softmax`
@with_kw mutable struct TRPOSolver
    episode_length::Int64 = 100 # Length of each episode
    terminate_horizon::Bool = false # Zero reward at end of horizon
    resume::Bool = false # Load saved policy and resume.

    # Optimization parameters
    η::Float64 = 1e-3 # learning rate η
    optimizer = ADAM(η) # Optimizer for the value function neural network
    δ = 0.01 # KL-divergence constraint
    value_iterations = 5 # Number of iterations to train the value function network

    # Generalized advantage estimate (GAE) parameters
    gae_γ::Float64 = 0.99 # GAE-Gamma
    gae_λ::Float64 = 0.95 # GAE-Lambda

    # Training parameters
    num_episodes::Int64 = 10 # Number of epochs of interaction (equal to number of policy updates)
    batch_size::Int64 = episode_length

    # Policy parameters
    hidden_layer_size::Int64 = 30 # Policy and value hidden layer size
    policy_type::Symbol = :continuous # :discrete or :continuous (TODO. infer)
    log_std::Float64 = log(1) # Initial log standard deviation used in the log-STD network

    # Environment parameters
    state_size::Int64 = 1 # Size of the state/observation space.
    action_size::Int64 = 1 # Size of the action space.

    # Frequencies and verbosity
    save::Bool = false
    save_frequency::Int64 = 10
    verbose_frequency::Int64 = 5
    verbose::Bool = false
    show_progress::Bool = true
end

mutable struct TRPOPlanner{P<:Union{MDP,POMDP}}
    solver::TRPOSolver
    mdp::P
    env::Union{MDPEnvironment,POMDPEnvironment}
    policy::Union{CategoricalPolicy,DiagonalGaussianPolicy,Nothing}
end

# Uses RLInterface environment wrapper
POMDPs.solve(solver::TRPOSolver, mdp::MDP) = TRPOPlanner(set_action_size!(solver, mdp), mdp, MDPEnvironment(mdp), nothing)
POMDPs.solve(solver::TRPOSolver, mdp::POMDP) = TRPOPlanner(set_action_size!(solver, mdp), mdp, POMDPEnvironment(mdp), nothing)

# NOTE: action and playback shared between TRPO and PPO (see PolicyOptimization.jl)

function policy_loss(policy, states::Array, actions::Array, advantages::Array, old_log_probs::Array)
    # Surrogate loss computation
    new_log_probs = log_prob(policy, states, actions)
    ratio = exp.(new_log_probs .- old_log_probs)
    π_loss = mean(ratio .* advantages)
    return π_loss
end

kl_loss(policy, states::Array, kl_vars) = mean(kl_divergence(policy, kl_vars, states))

value_loss(policy, states::Array, returns::Array) = mean((policy.value_net(states) .- returns).^2)

function linesearch!(solver, policy, step_dir, states, actions, advantages, old_log_probs, kl_vars, old_params, num_steps=10; α=0.5)
    old_loss = policy_loss(policy, states, actions, advantages, old_log_probs).data

    for i in 1:num_steps
        # Obtain new parameters
        new_params = old_params .+ (α^i .* step_dir)

        # Set the new parameters to the policy
        set_flat_params!(new_params, get_policy_net(policy))

        # Compute surrogate loss
        new_loss = policy_loss(policy, states, actions, advantages, old_log_probs).data

        # Compute kl divergence
        kl_div = kl_loss(policy, states, kl_vars).data

        # Output Statistics
        solver.verbose ? println("Old Loss : $old_loss") : nothing
        solver.verbose ? println("New Loss : $new_loss") : nothing
        solver.verbose ? println("KL Div : $kl_div") : nothing

        if new_loss >= old_loss && (kl_div <= solver.δ)
            solver.verbose ? println("Success.") : nothing
            set_flat_params!(new_params, get_policy_net(policy))
            return nothing
        end
    end

    set_flat_params!(old_params, get_policy_net(policy))
end

function trpo_update!(solver, policy, states, actions, advantages, returns, log_probs, kl_vars, old_params)
    model_params = get_policy_params(policy)
    policy_grads = Tracker.gradient(() -> policy_loss(policy, states, actions, advantages, log_probs), model_params)
    flat_policy_grads = get_flat_grads(policy_grads, get_policy_net(policy)).data

    x = conjugate_gradients(policy, states, kl_vars, flat_policy_grads, 10)
    solver.verbose ? println(minimum(x' * Hvp(policy, states, kl_vars, x))) : nothing

    step_dir = nothing
    try
        step_dir = sqrt.((2solver.δ) ./ (x' * Hvp(policy, states, kl_vars, x))) .* x
    catch err
        if err isa DomainError
            @info "Square root of a negative number received...Skipping update"
            return
        else
            throw(err)
        end
    end

    # Do a line search and update the parameters
    linesearch!(solver, policy, step_dir, states, actions, advantages, log_probs, kl_vars, old_params)

    # Update value function
    for _ in 1:solver.value_iterations
        value_params = get_value_params(policy)
        gs = Tracker.gradient(() -> value_loss(policy, states, returns), value_params)
        update!(solver.optimizer, value_params, gs)
    end
end


text:
using DiffEqBase, DiffEqOperators, ForwardDiff, LinearAlgebra, Test
const A = rand(300,300)
f(du,u) = mul!(du,A,u)
f(u) = A*u
x = rand(300)
v = rand(300)
du = similar(x)

cache1 = ForwardDiff.Dual{DiffEqOperators.JacVecTag}.(x, v)
cache2 = ForwardDiff.Dual{DiffEqOperators.JacVecTag}.(x, v)
@test DiffEqOperators.num_jacvec!(du, f, x, v) ≈ ForwardDiff.jacobian(f,similar(x),x)*v rtol=1e-6
@test DiffEqOperators.num_jacvec!(du, f, x, v, similar(v), similar(v)) ≈ ForwardDiff.jacobian(f,similar(x),x)*v rtol=1e-6
@test DiffEqOperators.num_jacvec(f, x, v) ≈ ForwardDiff.jacobian(f,similar(x),x)*v rtol=1e-6

@test DiffEqOperators.auto_jacvec!(du, f, x, v) ≈ ForwardDiff.jacobian(f,similar(x),x)*v
@test DiffEqOperators.auto_jacvec!(du, f, x, v, cache1, cache2) ≈ ForwardDiff.jacobian(f,similar(x),x)*v
@test DiffEqOperators.auto_jacvec(f, x, v) ≈ ForwardDiff.jacobian(f,similar(x),x)*v

f(du,u,p,t) = mul!(du,A,u)
f(u,p,t) = A*u
L = JacVecOperator(f,x)
@test L*x ≈ DiffEqOperators.auto_jacvec(f, x, x)
@test L*v ≈ DiffEqOperators.auto_jacvec(f, x, v)
@test mul!(du,L,v) ≈ DiffEqOperators.auto_jacvec(f, x, v)
DiffEqBase.update_coefficients!(L,v,nothing,nothing)
@test mul!(du,L,v) ≈ DiffEqOperators.auto_jacvec(f, v, v)

L = JacVecOperator(f,x,autodiff=false)
DiffEqBase.update_coefficients!(L,x,nothing,nothing)
@test L*x ≈ DiffEqOperators.num_jacvec(f, x, x)
@test L*v ≈ DiffEqOperators.num_jacvec(f, x, v)
@test mul!(du,L,v) ≈ DiffEqOperators.num_jacvec(f, x, v) rtol=1e-6
DiffEqBase.update_coefficients!(L,v,nothing,nothing)
@test mul!(du,L,v) ≈ DiffEqOperators.num_jacvec(f, v, v) rtol=1e-6

L2 = JacVecOperator{Float64}(f)
DiffEqBase.update_coefficients!(L2,x,nothing,nothing)
@test L2*x ≈ DiffEqOperators.auto_jacvec(f, x, x)
@test L2*v ≈ DiffEqOperators.auto_jacvec(f, v, v)
@test mul!(du,L2,x) ≈ DiffEqOperators.auto_jacvec(f, x, x)

L2 = JacVecOperator{Float64}(f,autodiff=false)
DiffEqBase.update_coefficients!(L2,x,nothing,nothing)
@test L2*x ≈ DiffEqOperators.num_jacvec(f, x, x)
@test L2*v ≈ DiffEqOperators.num_jacvec(f, v, v) rtol=1e-6

using OrdinaryDiffEq
function lorenz(du,u,p,t)
 du[1] = 10.0(u[2]-u[1])
 du[2] = u[1]*(28.0-u[3]) - u[2]
 du[3] = u[1]*u[2] - (8/3)*u[3]
end
u0 = [1.0;0.0;0.0]
tspan = (0.0,100.0)
ff1 = ODEFunction(lorenz,jac_prototype=JacVecOperator{Float64}(lorenz,u0))
ff2 = ODEFunction(lorenz,jac_prototype=JacVecOperator{Float64}(lorenz,u0,autodiff=false))

for ff in [ff1, ff2]
  prob = ODEProblem(ff,u0,tspan)
  @test solve(prob,TRBDF2()).retcode == :Success
  @test solve(prob,TRBDF2(linsolve=LinSolveGMRES())).retcode == :Success
  @test solve(prob,Exprb32()).retcode == :Success
  @test_broken sol = solve(prob,Rosenbrock23())
  @test_broken sol = solve(prob,Rosenbrock23(linsolve=LinSolveGMRES(tol=1e-10)))
end


text:
export up, down, Bit, Spin, Half, Clock, Potts, Continuous

abstract type AbstractSite{T} end

Base.show(io::IO, s::AbstractSite) = show(io, value(s))

Base.:(==)(lhs::AbstractSite, rhs::Number) = value(lhs) == rhs
Base.:(==)(lhs::Number, rhs::AbstractSite) = lhs == value(rhs)

abstract type IntegerSite{T} <: AbstractSite{T} end
abstract type BinarySite{T} <: IntegerSite{T} end
abstract type ContinuousSite{T} <: AbstractSite{T} end

"""
    Bit{T} <: BinarySite{T}
"""
struct Bit{T} <: BinarySite{T}
    value::T
end

struct Spin{T} <: BinarySite{T}
    value::T
end

struct Half{T} <: BinarySite{T}
    value::T
end

struct Clock{T, q} <: IntegerSite{T}
    value::T
end

struct Potts{T, q} <: IntegerSite{T}
    value::T
end

struct Continuous{T, d} <: ContinuousSite{T}
    value::T
end

"""
    up(site) -> site
    up(site_type) -> site

Up (highest value) tag for this label. e.g. `1` for `Bit`, `0.5` for `Half`.
"""
up(::ST) where {ST<: IntegerSite} = up(ST)
up(::Type{ST}) where {ST <: IntegerSite} = values(ST)[end]

"""
    down(site) -> site
    down(site_type) -> site

Down (lowest value) tag for this label. e.g. `0` for `Bit`, `-0.5` for `Half`.
"""
down(::ST) where {ST<: IntegerSite} = down(ST)
down(::Type{ST}) where {ST <: IntegerSite} = values(ST)[1]

"""
    values(site) -> site
    values(site_type) -> site

Returns a tuple of all possible values of the site type
"""
Base.values(::ST) where {ST <: AbstractSite} = values(ST)
Base.values(::Type{ST}) where {ST <: AbstractSite} = map(ST, _values(ST))
Base.values(::ST, i::Int) where {ST <: AbstractSite} = _values(ST)[i]

_values(::Type{Bit{T}}) where T = (zero(T), one(T))
_values(::Type{Spin{T}}) where T = (-one(T), one(T))
_values(::Type{Half{T}}) where T = (-0.5, 0.5)
_values(::Type{Clock{T, q}}) where {T, q} = Base.OneTo(q)
_values(::Type{Potts{T, q}}) where {T, q} = Tuple(-q:q)
# _values(::Type{Continuous{T}}) where {T} = # TODO

value(S::AbstractSite) = S.value

Base.length(::AbstractSite) = 1
Base.iterate(x::AbstractSite) = (x, nothing)
Base.iterate(x::AbstractSite, state) = nothing

Random.rand(rng::Random.AbstractRNG, sp::Random.SamplerType{Bit{T}}) where T = Bit{T}(rand(rng, Bool))
Random.rand(rng::Random.AbstractRNG, sp::Random.SamplerType{Spin{T}}) where T = Spin{T}(2 * rand(rng, Bool) - 1)
Random.rand(rng::Random.AbstractRNG, sp::Random.SamplerType{Half{T}}) where T = Half{T}(rand(rng, Bool) - 0.5)
Random.rand(rng::Random.AbstractRNG, sp::Random.SamplerType{Clock{T, q}}) where {T, q} = Clock{T, q}(rand(rng, 1:q))
Random.rand(rng::Random.AbstractRNG, sp::Random.SamplerType{Potts{T, q}}) where {T, q} = Potts{T, q}(rand(rng, -q:q))

Base.to_index(A::AbstractArray, i::Bit) = Int(value(i) + 1)
Base.to_index(A::AbstractArray, i::Spin) = Int(0.5 * (value(i) + 1) + 1)
Base.to_index(A::AbstractArray, i::Half) = Int(value(i) + 1.5)
Base.to_index(A::AbstractArray, i::Clock) = Int(value(i))
Base.to_index(A::AbstractArray, i::Potts) = Int(value(i) + q + 1)

Base.to_index(A::AbstractArray, I::AbstractArray{Bit{T}, N}) where {T, N} = Int(I) + 1

Base.getindex(t::Tuple, i::Bit) = getindex(t, value(i) + 1)
Base.getindex(t::Tuple, i::Spin) = getindex(t, Int(div(value(i) + 1, 2) + 1))

include("roundings.jl")
include("conversions.jl")
include("arraymath.jl")


text:
module PhysProps

using EzXML
using DataFrames

export
    # read from database, could decide to not export these at later point
      getDBRoot,
      name,
      value,
      units,
      parse_compound,
      parse_property

include("purecompbase.jl")
include("propsbase.jl")
include("properties.jl")
include("readDB.jl")

end # module


text:
function get_friendships_no_retweets_ids(auth::TwAuth, options::Dict)
  Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/friendships/no_retweets/ids.json", options))
end

function get_friendships_incoming(auth::TwAuth, options::Dict)
  Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/friendships/incoming.json", options))
end

function get_friendships_outgoing(auth::TwAuth, options::Dict)
  Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/friendships/outgoing.json", options))
end

function get_friendships_show(auth::TwAuth, options::Dict)
  Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/friendships/show.json", options))
end

function get_friendships_lookup(auth::TwAuth, options::Dict)
  Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/friendships/lookup.json", options))
end

max_stars_repo_name: nick-paul/JuliaTwitter

max_issues_repo_name: nick-paul/JuliaTwitter

max_forks_repo_name: nick-paul/JuliaTwitter


text:
# This code presents an example to illustrate how NLboxsolve.jl can be used

using NLboxsolve

function test_nlboxsolve_outplace()

    function test6(x)

        f = similar(x)

        for i in eachindex(x)
            f[i] = cos(x[i]) - 1.0
        end

        return f

    end

    # Solve test6 model

    n = 20
    x0 = ones(n)*0.5
    l = ones(n)*-0.5
    u = ones(n)*0.7

    soln_newton     = constrained_newton(test6,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_lm         = constrained_levenberg_marquardt(test6,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_lm_kyf     = constrained_levenberg_marquardt_kyf(test6,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_lm_ar      = constrained_levenberg_marquardt_ar(test6,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_dogleg     = constrained_dogleg_solver(test6,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_dogleg_bmp = constrained_dogleg_bmp_solver(test6,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_nk         = constrained_newton_krylov(test6,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_jfnk       = constrained_jacobian_free_newton_krylov(test6,x0,l,u,ftol=1e-15,xtol=1e-15)

    soln_newton_sparse     = constrained_newton_sparse(test6,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_lm_sparse         = constrained_levenberg_marquardt_sparse(test6,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_lm_kyf_sparse     = constrained_levenberg_marquardt_kyf_sparse(test6,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_lm_ar_sparse      = constrained_levenberg_marquardt_ar_sparse(test6,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_dogleg_sparse     = constrained_dogleg_solver_sparse(test6,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_dogleg_bmp_sparse = constrained_dogleg_bmp_solver_sparse(test6,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_nk_sparse         = constrained_newton_krylov_sparse(test6,x0,l,u,ftol=1e-15,xtol=1e-15)

    soln2 = nlboxsolve(test6,x0,l,u,xtol=1e-15,ftol=1e-15)

    test_one   = maximum(abs,soln_newton.zero       - soln2.zero) < 1e-8
    test_two   = maximum(abs,soln_lm.zero           - soln2.zero) < 1e-8
    test_three = maximum(abs,soln_lm_kyf.zero       - soln2.zero) < 1e-8
    test_four  = maximum(abs,soln_lm_ar.zero        - soln2.zero) < 1e-8
    test_five  = maximum(abs,soln_dogleg.zero       - soln2.zero) < 1e-8
    test_six   = maximum(abs,soln_dogleg_bmp.zero   - soln2.zero) < 1e-8
    test_seven = maximum(abs,soln_nk.zero           - soln2.zero) < 1e-8
    test_eight = maximum(abs,soln_jfnk.zero         - soln2.zero) < 1e-8

    test_nine     = maximum(abs,soln_newton_sparse.zero       - soln2.zero) < 1e-8
    test_ten      = maximum(abs,soln_lm_sparse.zero           - soln2.zero) < 1e-8
    test_eleven   = maximum(abs,soln_lm_kyf_sparse.zero       - soln2.zero) < 1e-8
    test_twelve   = maximum(abs,soln_lm_ar_sparse.zero        - soln2.zero) < 1e-8
    test_thirteen = maximum(abs,soln_dogleg_sparse.zero       - soln2.zero) < 1e-8
    test_fourteen = maximum(abs,soln_dogleg_bmp_sparse.zero   - soln2.zero) < 1e-8
    test_fifteen  = maximum(abs,soln_nk_sparse.zero           - soln2.zero) < 1e-8

    return test_one, test_two, test_three, test_four, test_five, test_six, test_seven, test_eight, test_nine, test_ten, test_eleven,
           test_twelve, test_thirteen, test_fourteen, test_fifteen

end

function test_nlboxsolve_inplace()

    function test6(f,x)

        for i in eachindex(x)
            f[i] = cos(x[i]) - 1.0
        end

    end

    # Solve test6 model

    n = 20
    x0 = ones(n)*0.5
    l = ones(n)*-0.5
    u = ones(n)*0.7

    soln_newton     = constrained_newton(test6,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_lm         = constrained_levenberg_marquardt(test6,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_lm_kyf     = constrained_levenberg_marquardt_kyf(test6,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_lm_ar      = constrained_levenberg_marquardt_ar(test6,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_dogleg     = constrained_dogleg_solver(test6,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_dogleg_bmp = constrained_dogleg_bmp_solver(test6,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_nk         = constrained_newton_krylov(test6,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_jfnk       = constrained_jacobian_free_newton_krylov(test6,x0,l,u,ftol=1e-15,xtol=1e-15)

    soln_newton_sparse     = constrained_newton_sparse(test6,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_lm_sparse         = constrained_levenberg_marquardt_sparse(test6,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_lm_kyf_sparse     = constrained_levenberg_marquardt_kyf_sparse(test6,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_lm_ar_sparse      = constrained_levenberg_marquardt_ar_sparse(test6,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_dogleg_sparse     = constrained_dogleg_solver_sparse(test6,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_dogleg_bmp_sparse = constrained_dogleg_bmp_solver_sparse(test6,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_nk_sparse         = constrained_newton_krylov_sparse(test6,x0,l,u,ftol=1e-15,xtol=1e-15)

    soln2 = nlboxsolve(test6,x0,l,u,xtol=1e-15,ftol=1e-15)

    test_one   = maximum(abs,soln_newton.zero       - soln2.zero) < 1e-8
    test_two   = maximum(abs,soln_lm.zero           - soln2.zero) < 1e-8
    test_three = maximum(abs,soln_lm_kyf.zero       - soln2.zero) < 1e-8
    test_four  = maximum(abs,soln_lm_ar.zero        - soln2.zero) < 1e-8
    test_five  = maximum(abs,soln_dogleg.zero       - soln2.zero) < 1e-8
    test_six   = maximum(abs,soln_dogleg_bmp.zero   - soln2.zero) < 1e-8
    test_seven = maximum(abs,soln_nk.zero           - soln2.zero) < 1e-8
    test_eight = maximum(abs,soln_jfnk.zero         - soln2.zero) < 1e-8

    test_nine     = maximum(abs,soln_newton_sparse.zero       - soln2.zero) < 1e-8
    test_ten      = maximum(abs,soln_lm_sparse.zero           - soln2.zero) < 1e-8
    test_eleven   = maximum(abs,soln_lm_kyf_sparse.zero       - soln2.zero) < 1e-8
    test_twelve   = maximum(abs,soln_lm_ar_sparse.zero        - soln2.zero) < 1e-8
    test_thirteen = maximum(abs,soln_dogleg_sparse.zero       - soln2.zero) < 1e-8
    test_fourteen = maximum(abs,soln_dogleg_bmp_sparse.zero   - soln2.zero) < 1e-8
    test_fifteen  = maximum(abs,soln_nk_sparse.zero           - soln2.zero) < 1e-8

    return test_one, test_two, test_three, test_four, test_five, test_six, test_seven, test_eight, test_nine, test_ten, test_eleven,
           test_twelve, test_thirteen, test_fourteen, test_fifteen

end

test_nlboxsolve_outplace()
test_nlboxsolve_inplace()


text:
## Algebraic identities

function mult_zeros(a, b::AbstractMatrix)
    size(a, 2) ≠ size(b, 1) &&
        throw(DimensionMismatch("Incompatible matrix multiplication dimensions"))
    return Zeros{promote_type(eltype(a), eltype(b))}(size(a, 1), size(b, 2))
end
function mult_zeros(a, b::AbstractVector)
    size(a, 2) ≠ size(b, 1) &&
        throw(DimensionMismatch("Incompatible matrix multiplication dimensions"))
    return Zeros{promote_type(eltype(a), eltype(b))}(size(a, 1))
end

const ZerosVecOrMat{T} = Union{Zeros{T,1}, Zeros{T,2}}
*(a::ZerosVecOrMat, b::AbstractMatrix) = mult_zeros(a, b)
*(a::AbstractMatrix, b::ZerosVecOrMat) = mult_zeros(a, b)
*(a::ZerosVecOrMat, b::AbstractVector) = mult_zeros(a, b)
*(a::AbstractVector, b::ZerosVecOrMat) = mult_zeros(a, b)
*(a::ZerosVecOrMat, b::ZerosVecOrMat) = mult_zeros(a, b)

function *(a::Adjoint{T, <:AbstractVector{T}}, b::Zeros{S, 1}) where {T, S}
    la, lb = length(a), length(b)
    if la ≠ lb
        throw(DimensionMismatch("dot product arguments have lengths $la and $lb"))
    end
    return zero(promote_type(T, S))
end
*(a::Adjoint{T, <:AbstractMatrix{T}} where T, b::Zeros{<:Any, 1}) = mult_zeros(a, b)

function *(a::Transpose{T, <:AbstractVector{T}}, b::Zeros{T, 1}) where T<:Real
    la, lb = length(a), length(b)
    if la ≠ lb
        throw(DimensionMismatch("dot product arguments have lengths $la and $lb"))
    end
    return zero(T)
end
*(a::Transpose{T, <:AbstractMatrix{T}}, b::Zeros{T, 1}) where T<:Real = mult_zeros(a, b)

+(a::Zeros) = a
-(a::Zeros) = a

# Zeros +/- Zeros
function +(a::Zeros{T}, b::Zeros{V}) where {T, V}
    size(a) ≠ size(b) && throw(DimensionMismatch("dimensions must match."))
    return Zeros{promote_type(T,V)}(size(a)...)
end
-(a::Zeros, b::Zeros) = -(a + b)

# Zeros +/- Fill and Fill +/- Zeros
function +(a::AbstractFill{T}, b::Zeros{V}) where {T, V}
    size(a) ≠ size(b) && throw(DimensionMismatch("dimensions must match."))
    return convert(AbstractFill{promote_type(T, V)}, a)
end
+(a::Zeros, b::AbstractFill) = b + a
-(a::AbstractFill, b::Zeros) = a + b
-(a::Zeros, b::AbstractFill) = a + (-b)

# Zeros +/- Array and Array +/- Zeros
function +(a::Zeros{T, N}, b::AbstractArray{V, N}) where {T, V, N}
    size(a) ≠ size(b) && throw(DimensionMismatch("dimensions must match."))
    return AbstractArray{promote_type(T,V),N}(b)
end
function +(a::Array{T, N}, b::Zeros{V, N}) where {T, V, N}
    size(a) ≠ size(b) && throw(DimensionMismatch("dimensions must match."))
    return AbstractArray{promote_type(T,V),N}(a)
end

function -(a::Zeros{T, N}, b::AbstractArray{V, N}) where {T, V, N}
    size(a) ≠ size(b) && throw(DimensionMismatch("dimensions must match."))
    return -b + a
end
-(a::Array{T, N}, b::Zeros{V, N}) where {T, V, N} = a + b

+(a::AbstractRange, b::Zeros) = b + a

function +(a::Zeros{T, 1}, b::AbstractRange) where {T}
    size(a) ≠ size(b) && throw(DimensionMismatch("dimensions must match."))
    Tout = promote_type(T, eltype(b))
    return convert(Tout, first(b)):convert(Tout, step(b)):convert(Tout, last(b))
end
function +(a::Zeros{T, 1}, b::UnitRange) where {T}
    size(a) ≠ size(b) && throw(DimensionMismatch("dimensions must match."))
    Tout = promote_type(T, eltype(b))
    return convert(Tout, first(b)):convert(Tout, last(b))
end

function -(a::Zeros{T, 1}, b::AbstractRange{V}) where {T, V}
    size(a) ≠ size(b) && throw(DimensionMismatch("dimensions must match."))
    return -b + a
end
-(a::AbstractRange{T}, b::Zeros{V, 1}) where {T, V} = a + b


text:
module LP

import ..MathOptFormat

import MathOptInterface
const MOI = MathOptInterface

MOI.Utilities.@model(InnerModel,
    (MOI.ZeroOne, MOI.Integer),
    (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan, MOI.Interval),
    (),
    (),
    (),
    (MOI.ScalarAffineFunction,),
    (),
    ()
)

struct Options
    maximum_length::Int
    warn::Bool
    warn_once::Bool
    warned_start::Set{Char}
    warned_illegal::Set{Char}
end

function get_options(m::InnerModel)
    default_options = Options(255, false, false, Set{Char}(), Set{Char}())
    return get(m.ext, :LP_OPTIONS, default_options)
end

"""
    Model(; kwargs...)

Create an empty instance of MathOptFormat.LP.Model.

Keyword arguments are:

 - `maximum_length::Int=255`: the maximum length for the name of a variable.
   lp_solve 5.0 allows only 16 characters, while CPLEX 12.5+ allow 255.

 - `warn::Bool=false`: print a warning when variables or constraints are renamed.

 - `warn_once::Bool=false`: print a warning when variables or constraints are
   renamed, but only once per kind of replacement (e.g., once per illegal
   character).
"""
function Model(;
    maximum_length::Int = 255, warn::Bool = false, warn_once::Bool = false
)
    model = InnerModel{Float64}()
    options = Options(maximum_length, warn, warn_once, Set{Char}(), Set{Char}())
    model.ext[:LP_OPTIONS] = options
    return model
end

function Base.show(io::IO, ::InnerModel)
    print(io, "A .LP-file model")
    return
end

# ==============================================================================
#
#   MOI.write_to_file
#
# ==============================================================================

const START_REG = r"^([\.0-9eE])"
const NAME_REG = r"([^a-zA-Z0-9\!\"\#\$\%\&\(\)\/\,\.\;\?\@\_\`\'\{\}\|\~])"

function sanitized_name(name::String, options::Options)
    m = match(START_REG, name)
    if m !== nothing
        plural = length(m.match) > 1

        if options.warn || (options.warn_once && !(m.match[1] in options.warned_start))
            @warn("Name $(name) cannot start with a period, a number, e, or E. " *
                  "Prepending an underscore to name.")
            push!(options.warned_start, m.match[1])
        end

        return sanitized_name("_" * name, options)
    end

    m = match(NAME_REG, name)
    if m !== nothing
        plural = length(m.match) > 1

        if options.warn || (options.warn_once && !(m.match[1] in options.warned_illegal))
            @warn("Name $(name) contains $(ifelse(plural, "", "an "))" *
                  "illegal character$(ifelse(plural, "s", "")): " *
                  "\"$(m.match)\". Removing the offending " *
                  "character$(ifelse(plural, "s", "")) from name.")
            push!(options.warned_illegal, m.match[1])
        end

        return sanitized_name(replace(name, NAME_REG => s"_"), options)
    end

    # Truncate at the end to fit as many characters as possible.
    if length(name) > options.maximum_length
        @warn("Name $(name) too long (length: $(length(name)); " *
              "maximum: $(options.maximum_length)). Truncating.")
        return sanitized_name(String(name[1:options.maximum_length]), options)
    end

    return name
end

function write_function(io::IO, model::InnerModel, func::MOI.SingleVariable, sanitized_names::Dict{MOI.VariableIndex, String})
    print(io, sanitized_names[func.variable])
    return
end

function write_function(io::IO, model::InnerModel, func::MOI.ScalarAffineFunction{Float64}, sanitized_names::Dict{MOI.VariableIndex, String})
    is_first_item = true
    if !(func.constant ≈ 0.0)
        Base.Grisu.print_shortest(io, func.constant)
        is_first_item = false
    end
    for term in func.terms
        if !(term.coefficient ≈ 0.0)
            if is_first_item
                Base.Grisu.print_shortest(io, term.coefficient)
                is_first_item = false
            else
                print(io, term.coefficient < 0 ? " - " : " + ")
                Base.Grisu.print_shortest(io, abs(term.coefficient))
            end

            print(io, " ", sanitized_names[term.variable_index])
        end
    end
    return
end

function write_constraint_suffix(io::IO, set::MOI.LessThan)
    print(io, " <= ", )
    Base.Grisu.print_shortest(io, set.upper)
    println(io)
    return
end

function write_constraint_suffix(io::IO, set::MOI.GreaterThan)
    print(io, " >= ", )
    Base.Grisu.print_shortest(io, set.lower)
    println(io)
    return
end

function write_constraint_suffix(io::IO, set::MOI.EqualTo)
    print(io, " = ", )
    Base.Grisu.print_shortest(io, set.value)
    println(io)
    return
end

function write_constraint_suffix(io::IO, set::MOI.Interval)
    print(io, " <= ", )
    Base.Grisu.print_shortest(io, set.upper)
    println(io)
    return
end

function write_constraint_prefix(io::IO, set::MOI.Interval)
    Base.Grisu.print_shortest(io, set.lower)
    print(io, " <= ")
    return
end

write_constraint_prefix(io::IO, set) = nothing

function write_constraint(io::IO, model::InnerModel, index, sanitized_names::Dict{MOI.VariableIndex, String}; write_name::Bool = true)
    func = MOI.get(model, MOI.ConstraintFunction(), index)
    set = MOI.get(model, MOI.ConstraintSet(), index)
    if write_name
        print(io, MOI.get(model, MOI.ConstraintName(), index), ": ")
    end
    write_constraint_prefix(io, set)
    write_function(io, model, func, sanitized_names)
    write_constraint_suffix(io, set)
end

const SCALAR_SETS = (
    MOI.LessThan{Float64}, MOI.GreaterThan{Float64}, MOI.EqualTo{Float64},
    MOI.Interval{Float64}
)

function write_sense(io::IO, model::InnerModel)
    if MOI.get(model, MOI.ObjectiveSense()) == MOI.MAX_SENSE
        println(io, "maximize")
    else
        println(io, "minimize")
    end
    return
end

function write_objective(io::IO, model::InnerModel, sanitized_names::Dict{MOI.VariableIndex, String})
    print(io, "obj: ")
    obj_func_type = MOI.get(model, MOI.ObjectiveFunctionType())
    obj_func = MOI.get(model, MOI.ObjectiveFunction{obj_func_type}())
    write_function(io, model, obj_func, sanitized_names)
    println(io)
    return
end

function MOI.write_to_file(model::InnerModel, io::IO)
    options = get_options(model)
    max_length = options.maximum_length
    # Ensure each variable has a unique name that does not infringe LP constraints.
    MathOptFormat.create_unique_names(model, warn = options.warn)
    sanitized_names = Dict{MOI.VariableIndex, String}()
    sanitized_names_set = Set{String}()
    for v in MOI.get(model, MOI.ListOfVariableIndices())
        proposed_sanitized_name = sanitized_name(MOI.get(model, MOI.VariableName(), v), options)
        # In case of duplicate names after sanitization, add a number at the end.
        if proposed_sanitized_name in sanitized_names_set
            # If the name is already too long, make some space for the suffix.
            if length(proposed_sanitized_name) >= max_length
                proposed_sanitized_name = String(proposed_sanitized_name[1:max_length - 2])
            end
            i = 1
            while proposed_sanitized_name * '_' * string(i) in sanitized_names_set
                i += 1

                # If the maximum length constraint would be broken with the *next* i,
                # truncate a bit more the proposed_sanitized_name.
                if length(proposed_sanitized_name * '_' * string(i)) > max_length
                    proposed_sanitized_name = String(proposed_sanitized_name[1:length(proposed_sanitized_name) - 1])
                end
            end
            proposed_sanitized_name *= '_' * string(i)
        end
        push!(sanitized_names_set, proposed_sanitized_name)
        sanitized_names[v] = proposed_sanitized_name
    end

    write_sense(io, model)
    write_objective(io, model, sanitized_names)
    println(io, "subject to")
    for S in SCALAR_SETS
        for index in MOI.get(model, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{Float64}, S}())
            write_constraint(io, model, index, sanitized_names; write_name = true)
        end
    end

    println(io, "Bounds")
    for S in SCALAR_SETS
        for index in MOI.get(model, MOI.ListOfConstraintIndices{MOI.SingleVariable, S}())
            write_constraint(io, model, index, sanitized_names; write_name = false)
        end
    end

    for (S, str_S) in [(MOI.Integer, "General"), (MOI.ZeroOne, "Binary")]
        indices = MOI.get(model, MOI.ListOfConstraintIndices{MOI.SingleVariable, S}())
        if length(indices) > 0
            println(io, str_S)
            for index in indices
                write_function(io, model, MOI.get(model, MOI.ConstraintFunction(), index), sanitized_names)
                println(io)
            end
        end
    end

    println(io, "End")

    return
end

# ==============================================================================
#
#   MOI.read_from_to_file
#
# ==============================================================================

function MOI.read_from_file(model::InnerModel, io::IO)
    error("Read from file is not implemented for LP files.")
end

end


text:
module TestDifferentiableGivens
# test
using Test
using LinearAlgebra
using LinearAlgebra: givensAlgorithm
using ForwardDiff: gradient, derivative

i = 3
g = x->givensAlgorithm(x, 2x)[i]
nexp = 16
@testset "differentiating givensAlgorithm" begin
    for _ in 1:nexp
        x = randn()
        gx = g(x)
        h = 1e-4
        y = x + h
        gy = g(y)
        gd_fd = (gy - gx) / h

        gd = derivative(g, x)

        tol = 1e-4
        @test abs(gd_fd - gd) < tol
    end
end

end


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate((c+d*x)^3*sech(a+b*x), x) == :(2*b^-1*(c+d*x)^3*arctan(exp(a+b*x))+-6*im*b^-4*d^3*Polylogarithms.polylog(4,-1*im*exp(a+b*x))+6*im*b^-4*d^3*Polylogarithms.polylog(4,im*exp(a+b*x))+-6*im*b^-3*d^2*(c+d*x)*Polylogarithms.polylog(3,im*exp(a+b*x))+-3*im*d*b^-2*(c+d*x)^2*Polylogarithms.polylog(2,-1*im*exp(a+b*x))+3*im*d*b^-2*(c+d*x)^2*Polylogarithms.polylog(2,im*exp(a+b*x))+6*im*b^-3*d^2*(c+d*x)*Polylogarithms.polylog(3,-1*im*exp(a+b*x)))
@test integrate((c+d*x)^2*sech(a+b*x), x) == :(2*b^-1*(c+d*x)^2*arctan(exp(a+b*x))+-2*im*b^-3*d^2*Polylogarithms.polylog(3,im*exp(a+b*x))+2*im*b^-3*d^2*Polylogarithms.polylog(3,-1*im*exp(a+b*x))+-2*im*d*b^-2*(c+d*x)*Polylogarithms.polylog(2,-1*im*exp(a+b*x))+2*im*d*b^-2*(c+d*x)*Polylogarithms.polylog(2,im*exp(a+b*x)))
@test integrate((c+d*x)*sech(a+b*x), x) == :(b^-1*(2c+2*d*x)*arctan(exp(a+b*x))+im*d*b^-2*Polylogarithms.polylog(2,im*exp(a+b*x))+-1*im*d*b^-2*Polylogarithms.polylog(2,-1*im*exp(a+b*x)))
@test integrate((c+d*x)^3*sech(a+b*x)^2, x) == :(b^-1*(c+d*x)^3+b^-1*(c+d*x)^3*tanh(a+b*x)+3//2*b^-4*d^3*Polylogarithms.polylog(3,-1*exp(2a+2*b*x))+-3*d*b^-2*(c+d*x)^2*log(1+exp(2a+2*b*x))+-3*b^-3*d^2*(c+d*x)*Polylogarithms.polylog(2,-1*exp(2a+2*b*x)))
@test integrate((c+d*x)^2*sech(a+b*x)^2, x) == :(b^-1*(c+d*x)^2+b^-1*(c+d*x)^2*tanh(a+b*x)+-1*b^-3*d^2*Polylogarithms.polylog(2,-1*exp(2a+2*b*x))+-2*d*b^-2*(c+d*x)*log(1+exp(2a+2*b*x)))
@test integrate(sech(a+b*x)^2*(c+d*x), x) == :(b^-1*(c+d*x)*tanh(a+b*x)+-1*d*b^-2*log(cosh(a+b*x)))
@test integrate((c+d*x)^3*sech(a+b*x)^3, x) == :(b^-1*(c+d*x)^3*arctan(exp(a+b*x))+(1/2)*b^-1*(c+d*x)^3*sech(a+b*x)*tanh(a+b*x)+-6*b^-3*d^2*(c+d*x)*arctan(exp(a+b*x))+-3*im*b^-4*d^3*Polylogarithms.polylog(2,im*exp(a+b*x))+-3*im*b^-4*d^3*Polylogarithms.polylog(4,-1*im*exp(a+b*x))+3*im*b^-4*d^3*Polylogarithms.polylog(2,-1*im*exp(a+b*x))+3*im*b^-4*d^3*Polylogarithms.polylog(4,im*exp(a+b*x))+3//2*d*b^-2*(c+d*x)^2*sech(a+b*x)+-3*im*b^-3*d^2*(c+d*x)*Polylogarithms.polylog(3,im*exp(a+b*x))+3*im*b^-3*d^2*(c+d*x)*Polylogarithms.polylog(3,-1*im*exp(a+b*x))+-3//2*im*d*b^-2*(c+d*x)^2*Polylogarithms.polylog(2,-1*im*exp(a+b*x))+3//2*im*d*b^-2*(c+d*x)^2*Polylogarithms.polylog(2,im*exp(a+b*x)))
@test integrate((c+d*x)^2*sech(a+b*x)^3, x) == :(b^-1*(c+d*x)^2*arctan(exp(a+b*x))+-1*b^-3*d^2*arctan(sinh(a+b*x))+im*b^-3*d^2*Polylogarithms.polylog(3,-1*im*exp(a+b*x))+d*b^-2*(c+d*x)*sech(a+b*x)+(1/2)*b^-1*(c+d*x)^2*sech(a+b*x)*tanh(a+b*x)+-1*im*b^-3*d^2*Polylogarithms.polylog(3,im*exp(a+b*x))+im*d*b^-2*(c+d*x)*Polylogarithms.polylog(2,im*exp(a+b*x))+-1*im*d*b^-2*(c+d*x)*Polylogarithms.polylog(2,-1*im*exp(a+b*x)))
@test integrate(sech(a+b*x)^3*(c+d*x), x) == :(b^-1*(c+d*x)*arctan(exp(a+b*x))+(1/2)*d*b^-2*sech(a+b*x)+(1/2)*im*d*b^-2*Polylogarithms.polylog(2,im*exp(a+b*x))+(1/2)*b^-1*(c+d*x)*sech(a+b*x)*tanh(a+b*x)+-1//2*im*d*b^-2*Polylogarithms.polylog(2,-1*im*exp(a+b*x)))
@test integrate(x*sech(x)^-3//2+-1//3*x*sech(x)^(1/2), x) == :(-4//9*sech(x)^-3//2+2//3*x*sech(x)^-1//2*sinh(x))
@test integrate(x*sech(x)^-5//2+-3//5*x*sech(x)^-1//2, x) == :(-4//25*sech(x)^-5//2+2//5*x*sech(x)^-3//2*sinh(x))
@test integrate(x*sech(x)^-7//2+-5//21*x*sech(x)^(1/2), x) == :(-20//63*sech(x)^-3//2+-4//49*sech(x)^-7//2+2//7*x*sech(x)^-5//2*sinh(x)+10//21*x*sech(x)^-1//2*sinh(x))
@test integrate(x^2*sech(x)^-3//2+-1//3*x^2*sech(x)^(1/2), x) == :(-8//9*x*sech(x)^-3//2+16//27*sech(x)^-1//2*sinh(x)+2//3*x^2*sech(x)^-1//2*sinh(x)+-16//27*im*cosh(x)^(1/2)*sech(x)^(1/2)*Elliptic.F((1/2)*im*x,2))

max_stars_repo_path: test/6 Hyperbolic functions/6.5 Hyperbolic secant/6.5.1 (c+d x)^m (a+b sech)^n.jl


text:
function hi!(x::Int)
    #this doesn't change input because Ints are immutable
    x = x + 2
    return x
end

function hi!(x::Array{T,1}) where T <: Real
    x .= x .+ 2
    return x
end


text:
using .ARPACK

## eigs

function eigs(A;nev::Integer=6, ncv::Integer=20, which::ASCIIString="LM",
          tol=0.0, maxiter::Integer=1000, sigma=nothing, v0::Vector=zeros((0,)),
          ritzvec::Bool=true, complexOP::Bool=false)
    n = chksquare(A)
    (n <= 6) && (nev = min(n-1, nev))
    ncv = blas_int(min(max(2*nev+2, ncv), n))

    sym   = issym(A)
    T     = eltype(A)
    cmplx = T<:Complex
    bmat  = "I"

    if !isempty(v0)
        length(v0)==n || throw(DimensionMismatch(""))
        eltype(v0)==T || error("Starting vector must have eltype $T")
    else
        v0=zeros(T,(0,))
    end

    if sigma===nothing # normal iteration
        mode = 1
        sigma = 0
        linop(x) = A * x
    else # inverse iteration with level shift
        C = lufact(sigma==0 ? A : A - sigma*eye(A))
        if cmplx
            mode = 3
            linop(x) = C\x
        else
            if !complexOP
                mode = 3
                linop(x) = real(C\x)
            else
                mode = 4
                linop(x) = imag(C\x)
            end
        end
    end

    # Compute the Ritz values and Ritz vectors
    (resid, v, ldv, iparam, ipntr, workd, workl, lworkl, rwork, TOL) =
       ARPACK.aupd_wrapper(T, linop, n, sym, cmplx, bmat, nev, ncv, which, tol, maxiter, mode, v0)

    # Postprocessing to get eigenvalues and eigenvectors
    ARPACK.eupd_wrapper(T, n, sym, cmplx, bmat, nev, which, ritzvec, TOL,
        resid, ncv, v, ldv, sigma, iparam, ipntr, workd, workl, lworkl, rwork)

end

max_stars_repo_name: tobig/julia

max_issues_repo_name: tobig/julia

max_forks_repo_name: tobig/julia