text:
# The goal of this visualization is to see how `f(z) = eᶻ` can be
# approximated through `g(z) = (1 + z/n)ⁿ` as `n → ∞`. We will
# use a slider to change the value of `n` and also place a 'test set'
# in the domain to see how it gets transformed.
using Cairo
using ComplexVisual
@ComplexVisual.import_huge

const CVG = try using ComplexVisualGtk; ComplexVisualGtk; catch nothing; end

const fontface = cv_fontface("DejaVu Sans")

# We start by creating the function to draw the axes on the left and right
function setup_axis(setup::CV_SceneSetupChain)
    label_style = cv_black → fontface → cv_fontsize(20)
    # these are the rules used to label the x-axis ticks
                                # main ticks on left, get labeled
    ticks1, ticks1h, ticks2 = "%.0f" ⇒ -6:0,
                                # secondary ticks on left, don't get labeled
                                "" ⇒ -5.5:0.5,
                                # ticks on right, get labeled
                                "%.0f" ⇒ -3:3
    # rules used to label the y-axis ticks
    # uses a Tuple to label specific irrational points
    ticks3 = ("-π" ⇒ -pi, "-π/2" ⇒ -pi/2, "0" ⇒ 0.0, "π/2" ⇒ pi/2, "π" ⇒ pi,)

    # main ticks are drawn slightly longer than secondary ticks on left
    app_l = CV_TickLabelAppearance(; label_style, tick_length=10)
    app_s = CV_TickLabelAppearance(; label_style, tick_length=6)

    # use the standard axis builder from `lrdomains.jl`
    return cv_setup_lr_axis(setup,
        # x-axis labels                     # y-axis labels
        (app_l ↦ ticks1, app_s ↦ ticks1h), (app_l ↦ ticks3,),
        (app_l ↦ ticks2,),                  (app_l ↦ ticks3,)
    )
end

# Now we need to create labels + style for our custom slider. We wish
# to show more detail in the slider around the first few numbers,
# and have a slider position "at infinity"
function get_slider_rulers()  # {{{
    # points where the large ticks and medium ticks are drawn
    major_values, minor_values = Float64.(10:10:50), Float64.(15:5:45)
    # all leftover integers between 1 and 50, will be drawn with small ticks
    mini_values = setdiff(setdiff(Float64.(10:1:50), minor_values), major_values)
    # These values will all be labeled
    intro_values = Float64.(1:1:9)

    # label the major values between 10 and 50, and also ∞
    major_ticks = (("%.0f" ⇒  major_values)..., "∞" ⇒ 55.0)
    # minor values do not get labeled
    minor_ticks, mini_ticks = "" ⇒ minor_values, "" ⇒ mini_values
    intro_ticks = "%.0f" ⇒ intro_values

    # at this point the style is set, using `CV_TickLabelAppearance`
    major_ruler = CV_TickLabelAppearance(; tick_length=10,
            label_style=cv_black → fontface → cv_fontsize(20)) ↦ major_ticks
    minor_ruler = CV_TickLabelAppearance(; tick_length=7) ↦ minor_ticks
    mini_ruler = CV_TickLabelAppearance(; tick_length=4) ↦ mini_ticks
    intro_ruler = CV_TickLabelAppearance(; tick_length=4,
            label_style=cv_black → fontface → cv_fontsize(10)) ↦ intro_ticks
    return (major_ruler, minor_ruler, mini_ruler, intro_ruler)
end # }}}

# We have the labels, so we can build the interactive slider
function create_slider(setup, slider_pos)

    # setup will be a `CV_SceneSetupChain` which has a field `layout`.
    # For more info see `lrdomains.jl`
    layout = setup.layout

    # internal getter functions to access the layouts for the
    # left and right portraits. Referred to as the `seen rectangles`
    can_domain_l = cv_get_can_domain_l(layout)
    can_codomain_l = cv_get_can_codomain_l(layout)

    # will be used to position ('anchor') the slider,
    # based on the size of the portraits
    width = can_codomain_l.rectangle.right - can_domain_l.rectangle.left
    mid = cv_half(can_codomain_l.rectangle.right + can_domain_l.rectangle.left)
    bottom = min(can_domain_l.rectangle.bottom, can_codomain_l.rectangle.bottom)
    height = 15

    # now we create the basic slider object.
    # for more info, see `slider.jl`
    slider_data = cv_create_hslider(
        width, height, 0.5, 55.5, get_slider_rulers();
        decoraction_with_layout_and_position_callback=(inner_layout, pos) ->
            cv_border(inner_layout, pos, 1))

    # this is the bounding box ('container') that contains the drawing info
    cont_slider = slider_data.slider_container
    cc_cont_slider = slider_data.container_context

    # draw the slider in our setup.layout,
    # in the center (referring to x-position),
    # 20pts above the bottom of the left-right layout (y-position)
    cont_slider_l = cv_add_canvas!(layout,
        cont_slider,
        cv_anchor(cont_slider, :slider_south, cont_slider, :south),
        (mid, bottom - Int32(20)))

    # interactive part of the slider: if the user clicks
    # to the right of the max value, snap to the rightmost
    # slider position (labeled ∞)
    set_slider_value = z -> begin
        slider_pos.value = (z > 53) ? 55.0 : max(1.0, min(z, 50.0))
        return CV_Response(;redraw_flag=true)
    end

    # a few style decisions
    bg_painter = cv_color(.8,.8,.8) ↦ CV_FillPainter()
    mark_painter = (cv_op_source → cv_color(0,0,1) → cv_linewidth(2)) ↦
        CV_ValueMarkPainter(slider_pos,
            0.0, imag(cont_slider.can_slider.corner_ul), false)

    # finally the slider is added to the `CV_SceneSetupChain`
    setup = cv_setup_hslider(setup, slider_data, cont_slider_l,
        bg_painter → mark_painter, set_slider_value)
    return setup
end

# Putting everything together: almost identical to the
# 'standard' left-right setup process (see `cv_scene_lr_std` in 'lrdomains.jl')
function create_scene(trafo,
        domain, codomain, slider_pos; cut_test=nothing, gap=80,
        axis_label_style=cv_color(0,0,0) →
                  cv_fontface("sans-serif") → cv_fontsize(20), padding=30) # {{{
    # cyclic state counter counts right-clicks
    layout = CV_StateLayout(CV_2DLayout(), CV_CyclicValue(2))
    # creates the standard left-right layout
    layout = cv_do_lr_layout(cv_add(layout, trafo, domain, codomain), gap)

    # create a basic `CV_SceneSetupChain` with a default left-right setup
    setup = cv_setup_cycle_state(CV_LRSetupChain(layout))

    # place the slider in the setup chain
    setup = create_slider(setup, slider_pos)

    # standard left-right setup, using our custom axes
    # all the functions used here are in `lrdomains.jl`
    setup = cv_setup_lr_painters(setup; cut_test)
    setup = setup_axis(setup)
    setup = cv_setup_lr_border(setup)
    padding > 0 && cv_add_padding!(setup.layout, padding)
    setup = cv_setup_domain_codomain_scene(setup)
    cv_scene_lr_start(setup)
    return setup.layout
end # }}}

# `CV_TranslateByOffset(Float64)` lets the user interact with something that
# only needs to move left and right
slider_pos = CV_TranslateByOffset(Float64)
slider_pos.value = 1.0

# There is a discontinuity in our goal function due to our use of Log
# To get a better visualization, we add the cut test
cut_func = cv_create_angle_cross_test(π, 0.0, Inf; δ=1e-2)

cut_test = (z, w) -> (slider_pos.value == round(slider_pos.value)) ?
    false : cut_func(z + slider_pos.value, w + slider_pos.value)

trafo = z -> (slider_pos.value == 55.0) ? exp(z) :
    exp(slider_pos.value * log(1+z/slider_pos.value))
domain   = CV_Math2DCanvas(-6.5 + 3.5im, 0.5 - 3.5im, 100)
codomain = CV_Math2DCanvas(-3.5 + (pi+0.1)im, 3.5 - (pi+0.1)im, 100)

scene = create_scene(trafo, domain, codomain, slider_pos; cut_test)

# finally, run the scene!
if CVG !== nothing
    handler = CVG.cvg_visualize(scene); CVG.cvg_wait_for_destroy(handler.window)

    CVG.cvg_close(handler); CVG.cv_destroy(scene)
end

# vim:syn=julia:cc=79:fdm=marker:sw=4:ts=4:


text:
#TODO have to make manager into a specific thing with components to optimize...
@inline @inbounds function set_optim_values!(m::AbstractLedger, vals)
    # if length(@entities_in(m[Optimize])) <= 11
    #     for i = 1:8
    #         m[Spin].data[i] = Spin(ϕ=vals[i])
    #     end
    #     for i = 1:2
    #         m[L].data[i] = L(ϕ=vals[i+8])
    #     end
    #     if length(@entities_in(m[Optimize])) == 11
    #         m[Pb].data[1] = Pb(vals[end])
    #     end
    # else
        for (i, e) in enumerate(@entities_in(m[Optimize]))

            if e ∈ m[Spin]
                m[e] = Spin(ϕ=vals[i])
            elseif e ∈ m[L]
                m[e] = L(ϕ=vals[i])
            else
                m[e] = Pb(vals[i])
            end
        end
    # end
end
@export function optim_value(m, e)
    if e ∈ m[Spin]
        return m[Spin][e].ϕ
    elseif e ∈ m[L]
        return m[L][e].ϕ
    else
        return m[Pb][e].pb
    end
end
@export optim_values(m) = map(x->optim_value(m, x), @entities_in(m[Optimize]))

@export function obj(x, m)
    set_optim_values!(m, x)
    # set_angles!(m, x)
    m[Etot][1].e = 0.0
    update(stage(m, :E_all), m)
    return m[Etot][1].e
end

function BlackBoxOptim.bboptimize(m::AbstractLedger, cs...)
    empty!(m[Optimize])
    comps_to_optimize = map(x->m[x], cs)
    for c in comps_to_optimize
        for e in @entities_in(c)
            m[e] = Optimize()
        end
    end
    vals = optim_values(m)
    orig_angles = copy(vals)
    # res = optimize(x->obj(x, m), vals, ParticleSwarm(fill(-π, length(vals)), fill(1.0π, length(vals)), 140), Optim.Options(iterations=80))
    obj2 = x -> begin
        set_optim_values!(m, x)

        res = optimize(x->obj(x, m), optim_values(m, cs), ConjugateGradient())
        return res.minimum
    end
    res = bboptimize(x->obj2(x),
                     SearchRange = (0.0, 2.0π),
                     NumDimensions=length(vals),
                     MaxTime=0.1, Method=:generating_set_search,TraceMode=:silent)
    return res
end

function Optim.optimize(m::AbstractLedger, cs...)
    empty!(m[Optimize])
    comps_to_optimize = map(x->m[x], cs)
    for c in comps_to_optimize
        for e in @entities_in(c)
            m[e] = Optimize()
        end
    end
    vals = optim_values(m)
    orig_angles = copy(vals)
    res = optimize(x->obj(x, m), orig_angles, method=ConjugateGradient(),g_tol=1e-5)
    return res
end
Optim.optimize(m::AbstractLedger) = optimize(m, Spin, L, Pb)

@export function optimize!(m::AbstractLedger, cs...)
    res = optimize(m, cs...)
    set_optim_values!(m, res.minimizer)
    return m
end

"Takes the states and spaces the angles such that they are uniformly distributed from the starting point to the end point."
function space_out_states!(ledgers)
    n = length(ledgers[1][Optimize])
    total_rotations = zeros(n)
    nl = length(ledgers)
    rotation_directions = Vector{Int}[]
    for i=2:nl
        rot = all_angles(ledgers[i]) .- all_angles(ledgers[i-1])
        total_rotations .+= abs.(rot)
        push!(rotation_directions, sign.(rot))
    end
    rotation_per_ledger = total_rotations./ (nl - 1)
    for i=2:length(ledgers)-1
        rot_directions = rotation_directions[i-1]
        set_optim_values!(ledgers[i], all_angles(ledgers[i-1]) .+ rot_directions .* rotation_per_ledger)
    end
end

"Takes a set of states and performs nudged elastic band optimization on them, returning the optimized images."
function neb(ledgers, c=0.1)
    opt_ledgers = deepcopy.(ledgers)

    n = length(ledgers[1][Optimize])
    #10 is because 10 angles
    nebrange = i -> n*(i-1)+1:n*i
    start_ϕs = zeros(n*(length(ledgers)-2))
    for (i, l) in enumerate(ledgers[2:end-1])
        all_angles!(view(start_ϕs, nebrange(i)), l)
    end

    angle_vec = zeros(n)

    res = optimize(start_ϕs, ConjugateGradient()) do x
        tot = 0.0
        for il in 1:div(length(x), n)
            l = opt_ledgers[il+1]
            tot += obj(view(x, nebrange(il)), l)
        end
        #elastic part

        tot += c*sum((all_angles!(angle_vec, ledgers[1]) .- view(x, 1:n)) .^ 2)
        for i in 2:div(length(x),n)
            tot += c * sum((view(x, nebrange(i)) .- view(x, nebrange(i-1))) .^ 2)
        end
        tot += c*sum((all_angles!(angle_vec, ledgers[end]) .- view(x, length(x)-(n-1):length(x))).^2)
        return tot
    end
    for i = 2:length(opt_ledgers)-1
        ix = i-1
        set_optim_values!(opt_ledgers[i], view(res.minimizer, nebrange(ix)))
    end
    return opt_ledgers
end

ext: jl


text:
using DistributedArrays

module Star

import DistributedArrays:localpart
localpart(x) = x
#=
scan takes as input a function f that should be applied on the processors,
the data 'data' that should be distributed among the processors and passed as
an argument to f,
and the rest of the arguments 'args' for f.
=#
function scan(f, data, args...)
  assert(is_tuple(args))
  @sync begin
    for i = 1 : nworkers()
      in = length(data) >= nworkers() ? data[i] : data
      @spawnat workers()[i] f(in..., Base.map(localpart, args)...)
    end
  end
end

#=
reduce takes as input a function f that should be applied on the processors,
and the arguments 'args' for f.
It collects the result of applying the function f at each processor into
an array(or arrays), and returns a tuple of the array(or arrays).
=#
function reduce(f, args...)
  assert(is_tuple(args))
  rref = [@spawnat p f(Base.map(localpart, args)...) for p in workers()]
  root_in = [fetch(r) for r in rref]
  sz = size(root_in[1], 1)
  #convert array of tuples into a 2D array
  in = Array(Any, nworkers(), sz)
  for i = 1 : nworkers()
    for j = 1 : sz
      in[i,j] = root_in[i][j]
    end
  end
  #dynamically generate the function arguments from the 2D array
  ex = "$(in[:, 1])"
  for j = 2 : sz
    ex = ex * ", $(in[:, j])"
  end
  return eval(parse(ex))
end

function is_tuple(x::Tuple)
  true
end

function is_tuple(x)
  false
end

end


text:
function g(x,y,z)
  x(y,z)
end
println(g(+,2,3)) # prints 5

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
import Term.markup: has_markup, extract_markup

@testset "\e[34mMarkup - detection" begin
    # with no markup
    strings = [
        "asdasdasdsad",
        "asdasda[asdasdasd",
        "asdasdad[[asdasdasdas]]asdasd",
        "sadas[asdasdasd\nasdsadassadsa[[asdsadas]]asdsdadas",
    ]

    for str in strings
        @test has_markup(str) == false
        @test length(extract_markup(str)) == 0
    end

    # with markup
    function ntags(tag)
        return if length(tag.inner_tags) > 0
            length(tag.inner_tags) + sum(ntags.(tag.inner_tags))
        else
            0
        end
    end
    tottags(tags) = length(tags) + sum(ntags.(tags))

    strings = [
        ("[red]asdasdsa[/red]", 1),
        ("[red]dasdasdsa[green]asdasda[/green]asdasdasp[/red]", 2),
        ("[blue]onasdasdsa[green]asdasdasp[bold]asdasda", 3),
    ]

    for (str, n) in strings
        @test has_markup(str) == true
        @test tottags(extract_markup(str)) == n
    end
end


text:
module PowerSystem

using Printf
using SparseArrays
using ..ExaPF: ParsePSSE, ParseMAT, IndexSet, Spmat

import Base: show

const PQ_BUS_TYPE = 1
const PV_BUS_TYPE = 2
const REF_BUS_TYPE  = 3

"""
    AbstractPowerSystem

First layer of the package. Store the topology of a given
transmission network, including:

- the power injection at each bus
- the admittance matrix
- the default voltage at each bus

Data are imported either from a matpower file, or a PSSE file.

"""
abstract type AbstractPowerSystem end

"""
    AbstractNetworkElement

Abstraction for all physical elements being parts of a `AbstractPowerSystem`.
Elements are divided in

- transmission lines (`Lines`)
- buses (`Buses`)
- generators (`Generators`)

"""
abstract type AbstractNetworkElement end

"""
    Buses <: AbstractNetworkElement

Buses of a transmission network.
"""
struct Buses <: AbstractNetworkElement end

"""
    Lines <: AbstractNetworkElement

Lines of a transmission network.
"""
struct Lines <: AbstractNetworkElement end

"""
    Generator <: AbstractElement

Generator in a transmission network
"""
struct Generator <: AbstractNetworkElement end

"""
    AbstractNetworkAttribute

Attribute of a `AbstractPowerSystem`.

"""
abstract type AbstractNetworkAttribute end

"""
    NumberOfBuses <: AbstractNetworkAttribute

Number of buses in a `AbstractPowerSystem`.
"""
struct NumberOfBuses <: AbstractNetworkAttribute end

"""
    NumberOfLines <: AbstractNetworkAttribute

Number of lines in a `AbstractPowerSystem`.
"""
struct NumberOfLines <: AbstractNetworkAttribute end

"""
    NumberOfGenerators <: AbstractNetworkAttribute

Number of generators in a `AbstractPowerSystem`.
"""
struct NumberOfGenerators <: AbstractNetworkAttribute end

"""
    NumberOfPVBuses <: AbstractNetworkAttribute

Number of PV buses in a `AbstractPowerSystem`.
"""
struct NumberOfPVBuses <: AbstractNetworkAttribute end

"""
    NumberOfPQBuses <: AbstractNetworkAttribute

Number of PQ buses in a `AbstractPowerSystem`.
"""
struct NumberOfPQBuses <: AbstractNetworkAttribute end

"""
    NumberOfSlackBuses <: AbstractNetworkAttribute

Number of slack buses in a `AbstractPowerSystem`.
"""
struct NumberOfSlackBuses <: AbstractNetworkAttribute end

"""
    BaseMVA <: AbstractNetworkAttribute

Base MVA of the network.
"""
struct BaseMVA <: AbstractNetworkAttribute end

abstract type AbstractIndexing <: AbstractNetworkAttribute end

"""
    PVIndexes <: AbstractIndexing

Indexes of the PV buses in a `AbstractPowerSystem`.
"""
struct PVIndexes <: AbstractIndexing end

"""
    PQIndexes <: AbstractIndexing

Indexes of the PQ buses in a `AbstractPowerSystem`.
"""
struct PQIndexes <: AbstractIndexing end

"""
    SlackIndexes <: AbstractIndexing

Indexes of the slack buses in a `AbstractPowerSystem`.
"""
struct SlackIndexes <: AbstractIndexing end

"""
    GeneratorIndexes <: AbstractIndexing

Indexes of the generators in a `AbstractPowerSystem`.
"""
struct GeneratorIndexes <: AbstractIndexing end

# TODO: replace const *_BUS_TYPE with this enum
@enum BusType begin
    PQ=1
    PV=2
    Slack=3
end

"""
    AbstractNetworkValues

Numerical values attached to the different attributes of the network.
"""
abstract type AbstractNetworkValues end

"""
    VoltageMagnitude <: AbstractNetworkValues

Magnitude `|v|` of the voltage `v = |v| exp(i θ)`.
"""
struct VoltageMagnitude <: AbstractNetworkValues end

"""
    VoltageAngle <: AbstractNetworkValues

Angle `θ` of the voltage `v = |v| exp(i θ)`.
"""
struct VoltageAngle <: AbstractNetworkValues end

"""
    ActivePower <: AbstractNetworkValues

Active power `P` of the complex power `S = P + iQ`.
"""
struct ActivePower <: AbstractNetworkValues end

"""
    ReactivePower <: AbstractNetworkValues

Reactive power `Q` of the complex power `S = P + iQ`.
"""
struct ReactivePower <: AbstractNetworkValues end

# Templating
"""
    get(pf::AbstractPowerSystem, attr::AbstractNetworkAttribute)

Return value of attribute `attr` in the `AbstractPowerSystem` object `pf`.

    get(pf::AbstractPowerSystem, attr::AbstractIndexing)

Return indexing corresponding to a subset of the buses.

## Examples

```julia
npq = get(pf, NumberOfPQBuses())
npv = get(pf, NumberOfPVBuses())
index_pv = get(pf, PVIndexes()())
index_gen = get(pf, GeneratorIndexes()())

```
"""
function get end

"""
    bounds(pf::AbstractPowerSystem, attr::AbstractNetworkAttribute, val::AbstractNetworkValues)

Return lower and upper bounds corresponding to the admissible values
of the `AbstractNetworkAttribute` `attr`.

## Examples

```julia
p_min, p_max = bounds(pf, Generator(), ActivePower())
v_min, v_max = bounds(pf, Buses(), VoltageMagnitude())

```
"""
function bounds end

include("topology.jl")
include("power_network.jl")

end


text:
insts = OrderedDict()

insts["test"] = [
    ((4, MatNegLog()),),
    ((4, MatNegLogDirect()),),
    ((6, MatNegRtdet()),),
    ((6, MatNegRtdetEFExp()),),
]

inversecovariance_insts(f::MatSpecExt, max_d::Int) = [((d, f),) for d in vcat(5, 5:max_d)]

insts["nat"] = inversecovariance_insts(MatNegLog(), 16)
insts["ext"] = inversecovariance_insts(MatNegLogDirect(), 16)

return (InverseCovariance, insts)


text:
"""
    plot_code(code::TensorNetworkCode; use_coords=true)

Return a plot, as returned by `GraphPlot.gplot`, of the code. If used in a Jupyter
notebook the plot is rendered after the cell.

Physical qubit nodes are coloured red and virtual tensors are coloured green or blue:
green if they have no logical qubits and blue if they do.

# Examples
```julia
julia> code = TensorNetworkCode(five_qubit_code());

julia> plot_code(code)
```
![code plot](./plotting_functions/plot_code.svg)
"""
function plot_code(code::TensorNetworkCode; use_coords=true)
    Nodes = nodes(code.code_graph)
    graph = _light_graph(code)

    # nodelabels = [(node_types(code,node) == "physical") ? node : " "
    #     for node in Nodes]

    # nodecolours = [(node_types(code,node) == "physical") ?
    #     "tomato" : "lightgreen"
    #     for node in Nodes]

    nodelabels = fill("",length(Nodes))
    nodecolours = fill("",length(Nodes))

    for α in 1:length(Nodes)
        node = Nodes[α]
        if node_types(code,node) == "physical"
            nodelabels[α] = string(node)
            nodecolours[α] = "tomato"
        else
            code_name = node_types(code,node)
            seed_code = code.seed_codes[code_name]
            if length(seed_code.logicals) == 0
                nodecolours[α] = "lightgreen"
            else
                nodecolours[α] = "dodgerblue"
            end
        end
    end

    locs_x = Float64[]
    locs_y = Float64[]
    for node in Nodes
        locs = coords(code,node)
        push!(locs_x,locs[1])
        push!(locs_y,locs[2])
    end

    if use_coords
        return gplot(graph,locs_x,locs_y,nodefillc=nodecolours,nodelabel=nodelabels)
        # display(gplot(graph,locs_x,locs_y,nodefillc=nodecolours,nodelabel=nodelabels))
    else
        return gplot(graph,nodefillc=nodecolours,nodelabel=nodelabels)
        # display(gplot(graph,nodefillc=nodecolours,nodelabel=nodelabels))
    end
end

"""
    plot_operator(code::TensorNetworkCode, operator::AbstractVector{Int}; use_coords=true)

Return a plot, as returned by `GraphPlot.gplot`, of the operator. If used in a Jupyter
notebook the plot is rendered after the cell.

# Examples
```julia
julia> code = TensorNetworkCode(steane_code());

julia> plot_operator(code, code.stabilizers[1])
```
![operator plot](./plotting_functions/plot_operator.svg)
"""
function plot_operator(
    code::TensorNetworkCode, operator::AbstractVector{Int};
    use_coords = true
)
    Nodes = nodes(code.code_graph)
    graph = _light_graph(code)

    nodelabels = []
    nodecolours = String[]
    for node in Nodes
        node_type = node_types(code,node)

        if node_type == "physical"
            pauli = pauli_rep_change(operator[node])
            push!(nodelabels,pauli)
            if pauli == 'I'
                push!(nodecolours,"goldenrod1")
            else
                push!(nodecolours,"tomato")
            end
        else
            seed_code = code.seed_codes[node_type]
            if length(seed_code.logicals) == 0
                push!(nodecolours,"lightgreen")
            else
                push!(nodecolours,"dodgerblue")
            end
            push!(nodelabels," ")
        end
    end

    locs_x = Float64[]
    locs_y = Float64[]
    for node in Nodes
        locs = coords(code,node)
        push!(locs_x,locs[1])
        push!(locs_y,locs[2])
    end

    if use_coords
        return gplot(graph,locs_x,locs_y,nodefillc=nodecolours,nodelabel=nodelabels)
        # display(gplot(graph,locs_x,locs_y,nodefillc=nodecolours,nodelabel=nodelabels))
    else
        return gplot(graph,nodefillc=nodecolours,nodelabel=nodelabels)
        # display(gplot(graph,nodefillc=nodecolours,nodelabel=nodelabels))
    end
end

"""
    _light_graph(code_graph::CodeGraph)

Given a `CodeGraph`, returns a `LightGraph` for plotting.
"""
function _light_graph(code_graph::CodeGraph)

    Nodes = nodes(code_graph)
    virtual_nodes = filter(x->x<0,Nodes)
    n_virtual = length(virtual_nodes)

    Edges = edges(code_graph)

    output = Graph(length(Nodes))
    for edge in Edges
        graph_edge = _edge_to_graph_edge(edge,n_virtual)
        add_edge!(output,graph_edge...)
    end

    return output
end
_light_graph(code::TensorNetworkCode) = _light_graph(code.code_graph)

"""
    _edge_to_graph_edge(edge::Set{Int64},n_virtual::Int64)

Because `LightGraph` indices start at 1, but `CodeGraph` indices
include negative numbers, this converts `CodeGraph` edges to have
positive values.
"""
function _edge_to_graph_edge(edge::Set{Int64},n_virtual::Int64)

    output = [edge...] .+ n_virtual
    for α in 1:2
        if output[α] < n_virtual
            output[α] += 1
        end
    end

    return output
end


text:
# # Phillips model of Baroclinic Instability
#
#md # This example can be run online via [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/multilayerqg_2layer.ipynb).
#md # Also, it can be viewed as a Jupyter notebook via [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/multilayerqg_2layer.ipynb).
#
# A simulation of the growth of barolinic instability in the Phillips 2-layer model
# when we impose a vertical mean flow shear as a difference $\Delta U$ in the
# imposed, domain-averaged, zonal flow at each layer.

using FourierFlows, Plots, Printf

using FFTW: rfft, irfft
import GeophysicalFlows.MultilayerQG
import GeophysicalFlows.MultilayerQG: energies

# ## Choosing a device: CPU or GPU

dev = CPU()     # Device (CPU/GPU)
nothing # hide

# ## Numerical parameters and time-stepping parameters

nx = 128        # 2D resolution = nx^2
ny = nx

stepper = "FilteredRK4"  # timestepper
     dt = 6e-3           # timestep
 nsteps = 7000           # total number of time-steps
 nsubs  = 25             # number of time-steps for plotting (nsteps must be multiple of nsubs)
nothing # hide

# ## Physical parameters
Lx = 2π         # domain size
 μ = 5e-2       # bottom drag
 β = 5          # the y-gradient of planetary PV

nlayers = 2      # number of layers
f₀, g = 1, 1     # Coriolis parameter and gravitational constant
 H = [0.2, 0.8]  # the rest depths of each layer
 ρ = [4.0, 5.0]  # the density of each layer

 U = zeros(nlayers) # the imposed mean zonal flow in each layer
 U[1] = 1.0
 U[2] = 0.0
nothing # hide

# ## Problem setup
# We initialize a `Problem` by providing a set of keyword arguments,
prob = MultilayerQG.Problem(nlayers, dev; nx=nx, Lx=Lx, f₀=f₀, g=g, H=H, ρ=ρ, U=U, dt=dt, stepper=stepper, μ=μ, β=β)
nothing # hide

# and define some shortcuts.
sol, cl, pr, vs, gr = prob.sol, prob.clock, prob.params, prob.vars, prob.grid
x, y = gr.x, gr.y
nothing # hide

# ## Setting initial conditions

# Our initial condition is some small amplitude random noise. We smooth our initial
# condidtion using the `timestepper`'s high-wavenumber `filter`.

q_i  = 4e-3randn((nx, ny, nlayers))
qh_i = prob.timestepper.filter .* rfft(q_i, (1, 2)) # only apply rfft in dims=1, 2
q_i  = irfft(qh_i, gr.nx, (1, 2)) # only apply irfft in dims=1, 2

MultilayerQG.set_q!(prob, q_i)
nothing # hide

# ## Diagnostics

# Create Diagnostics -- `energies` function is imported at the top.
E = Diagnostic(energies, prob; nsteps=nsteps)
diags = [E] # A list of Diagnostics types passed to "stepforward!" will  be updated every timestep.
nothing # hide

# ## Output

# We choose folder for outputing `.jld2` files and snapshots (`.png` files).
filepath = "."
plotpath = "./plots_2layer"
plotname = "snapshots"
filename = joinpath(filepath, "2layer.jld2")
nothing # hide

# Do some basic file management
if isfile(filename); rm(filename); end
if !isdir(plotpath); mkdir(plotpath); end
nothing # hide

# And then create Output
get_sol(prob) = sol # extracts the Fourier-transformed solution
function get_u(prob)
  @. v.qh = sol
  streamfunctionfrompv!(v.ψh, v.qh, p, g)
  @. v.uh = -im * g.l * v.ψh
  invtransform!(v.u, v.uh, p)
  return v.u
end
out = Output(prob, filename, (:sol, get_sol), (:u, get_u))
nothing #hide

# ## Visualizing the simulation

# We define a function that plots the potential vorticity field and the evolution
# of energy and enstrophy.

symlims(data) = maximum(abs.(extrema(data))) |> q -> (-q, q)

function plot_output(prob)

  l = @layout grid(2, 3)
  p = plot(layout=l, size = (1000, 600), dpi=150)

  for m in 1:nlayers
    heatmap!(p[(m-1)*3+1], x, y, vs.q[:, :, m],
         aspectratio = 1,
              legend = false,
                   c = :balance,
               xlims = (-gr.Lx/2, gr.Lx/2),
               ylims = (-gr.Ly/2, gr.Ly/2),
               clims = symlims,
              xticks = -3:3,
              yticks = -3:3,
              xlabel = "x",
              ylabel = "y",
               title = "q_"*string(m),
          framestyle = :box)

    contourf!(p[(m-1)*3+2], x, y, vs.ψ[:, :, m],
              levels = 8,
         aspectratio = 1,
              legend = false,
                   c = :viridis,
               xlims = (-gr.Lx/2, gr.Lx/2),
               ylims = (-gr.Ly/2, gr.Ly/2),
               clims = symlims,
              xticks = -3:3,
              yticks = -3:3,
              xlabel = "x",
              ylabel = "y",
               title = "ψ_"*string(m),
          framestyle = :box)
  end

  plot!(p[3], 2,
             label = ["KE1" "KE2"],
            legend = :bottomright,
         linewidth = 2,
             alpha = 0.7,
             xlims = (-0.1, 2.35),
             ylims = (5e-10, 1e0),
            yscale = :log10,
            yticks = 10.0.^(-9:2:0),
            xlabel = "μt")

  plot!(p[6], 1,
             label = "PE",
            legend = :bottomright,
         linecolor = :red,
         linewidth = 2,
             alpha = 0.7,
             xlims = (-0.1, 2.35),
             ylims = (1e-10, 1e0),
            yscale = :log10,
            yticks = 10.0.^(-10:2:0),
            xlabel = "μt")

end
nothing # hide

# ## Time-stepping the `Problem` forward

# Finally, we time-step the `Problem` forward in time.

p = plot_output(prob)

startwalltime = time()

anim = @animate for j=0:Int(nsteps/nsubs)

  cfl = cl.dt*maximum([maximum(vs.u)/gr.dx, maximum(vs.v)/gr.dy])

  log = @sprintf("step: %04d, t: %.1f, cfl: %.2f, KE1: %.4f, KE2: %.4f, PE: %.4f, walltime: %.2f min", cl.step, cl.t, cfl, E.data[E.i][1][1], E.data[E.i][1][2], E.data[E.i][2][1], (time()-startwalltime)/60)

  if j%(1000/nsubs)==0; println(log) end

  for m in 1:nlayers
    p[(m-1)*3+1][1][:z] = @. vs.q[:, :, m]
    p[(m-1)*3+2][1][:z] = @. vs.ψ[:, :, m]
  end

  push!(p[3][1], μ*E.t[E.i], E.data[E.i][1][1])
  push!(p[3][2], μ*E.t[E.i], E.data[E.i][1][2])
  push!(p[6][1], μ*E.t[E.i], E.data[E.i][2][1])

  stepforward!(prob, diags, nsubs)
  MultilayerQG.updatevars!(prob)
end

mp4(anim, "multilayerqg_2layer.mp4", fps=18)

# ## Save

# Finally save the last snapshot.
savename = @sprintf("%s_%09d.png", joinpath(plotpath, plotname), cl.step)
savefig(savename)


text:
module Domains

import ..Geometry
using IntervalSets
export RectangleDomain

abstract type AbstractDomain end

abstract type HorizontalDomain <: AbstractDomain end
abstract type VerticalDomain <: AbstractDomain end

const BCTagType = Union{Nothing, Tuple{Symbol, Symbol}}

struct IntervalDomain{CT, B} <: VerticalDomain where {
    CT <: Geometry.Abstract1DPoint{FT},
    B <: BCTagType,
} where {FT}
    coord_min::CT
    coord_max::CT
    boundary_tags::B
end

"""
    IntervalDomain(coord⁻, coord⁺; boundary_tags::Tuple{Symbol,Symbol})

Construct a `IntervalDomain`, the closed interval is given by `coord⁻`, `coord⁺` coordinate arguments.
Because `IntervalDomain` does not support periodic boundary conditions, the `boundary_tags` keyword arugment must be supplied.
"""
function IntervalDomain(
    coord⁻::Geometry.Abstract1DPoint,
    coord⁺::Geometry.Abstract1DPoint;
    boundary_tags::Tuple{Symbol, Symbol},
)
    coords = promote(coord⁻, coord⁺)
    IntervalDomain(first(coords), last(coords), boundary_tags)
end

"""
    IntervalDomain(coords::ClosedInterval; boundary_tags::Tuple{Symbol,Symbol})

Construct a `IntervalDomain`, over the closed coordinate interval `coords`
Because `IntervalDomain` does not support periodic boundary conditions, the `boundary_tags` keyword arugment must be supplied.
"""
IntervalDomain(coords::ClosedInterval; boundary_tags::Tuple{Symbol, Symbol}) =
    IntervalDomain(coords.left, coords.right, boundary_tags)

coordinate_type(::IntervalDomain{CT}) where {CT} = CT
Base.eltype(domain::IntervalDomain) = coordinate_type(domain)

function Base.show(io::IO, domain::IntervalDomain)
    print(
        io,
        "IntervalDomain($(domain.coord_min) .. $(domain.coord_max), boundary_tags = $(domain.boundary_tags))",
    )
end
# coordinates (x1,x2)

struct RectangleDomain{CT, B1, B2} <: HorizontalDomain where {
    CT <: Geometry.Abstract2DPoint{FT},
    B1 <: BCTagType,
    B2 <: BCTagType,
} where {FT}
    x1x2min::CT
    x1x2max::CT
    x1boundary::B1
    x2boundary::B2
end

"""
    RectangleDomain(x1::ClosedInterval, x2::ClosedInterval;
        x1boundary::Tuple{Symbol,Symbol},
        x2boundary::Tuple{Symbol,Symbol},
        x1periodic = false,
        x2periodic = false,
    )

Construct a `RectangularDomain` in the horizontal.
If a given x1 or x2 boundary is not periodic, then `x1boundary` or `x2boundary` boundary tag keyword arguments must be supplied.
"""
function RectangleDomain(
    x1::ClosedInterval{X1CT},
    x2::ClosedInterval{X2CT};
    x1periodic = false,
    x2periodic = false,
    x1boundary::BCTagType = nothing,
    x2boundary::BCTagType = nothing,
) where {X1CT <: Geometry.Abstract1DPoint, X2CT <: Geometry.Abstract1DPoint}
    UX1CT = Geometry.unionalltype(X1CT)
    UX2CT = Geometry.unionalltype(X2CT)
    if UX1CT === UX2CT
        throw(
            ArgumentError(
                "x1 and x2 domain axis coordinates cannot be the same type: `$(UX1CT)`",
            ),
        )
    end
    if !x1periodic && !(x1boundary isa Tuple{Symbol, Symbol})
        throw(
            ArgumentError(
                "if `x1periodic=false` then an `x1boundary::Tuple{Symbol,Symbol}` boundary tag keyword argument must be supplied, got: $(x1boundary)",
            ),
        )
    end
    if !x2periodic && !(x2boundary isa Tuple{Symbol, Symbol})
        throw(
            ArgumentError(
                "if `x2periodic=false` then an `x2boundary::Tuple{Symbol,Symbol}` boundary tag keyword argument must be supplied, got: $(x2boundary)",
            ),
        )
    end
    x1x2min, x1x2max = promote(
        Geometry.product_coordinates(x1.left, x2.left),
        Geometry.product_coordinates(x1.right, x2.right),
    )
    RectangleDomain(
        x1x2min,
        x1x2max,
        x1periodic ? nothing : x1boundary,
        x2periodic ? nothing : x2boundary,
    )
end

function Base.show(io::IO, domain::RectangleDomain{CT}) where {CT}
    x1min = Geometry.coordinate(domain.x1x2min, 1)
    x2min = Geometry.coordinate(domain.x1x2min, 2)
    x1max = Geometry.coordinate(domain.x1x2max, 1)
    x2max = Geometry.coordinate(domain.x1x2max, 2)
    print(io, "RectangleDomain($(x1min)..$(x1max), $(x2min)..$(x2max)")
    if domain.x1boundary === nothing
        print(io, ", x1periodic=true")
    else
        print(io, ", x1boundary=$(domain.x1boundary)")
    end
    if domain.x2boundary === nothing
        print(io, ", x2periodic=true")
    else
        print(io, ", x2boundary=$(domain.x2boundary)")
    end
    print(io, ")")
end
coordinate_type(::RectangleDomain{CT}) where {CT} = CT

# coordinates (-pi/2 < lat < pi/2, -pi < lon < pi)
struct SphereDomain{FT} <: HorizontalDomain where {FT <: AbstractFloat}
    radius::FT
end

coordinate_type(::SphereDomain{FT}) where {FT} = Geometry.Cartesian123Point{FT}

struct CubePanelDomain{FT} <: HorizontalDomain where {FT <: AbstractFloat} end

coordinate_type(::CubePanelDomain{FT}) where {FT} =
    Geometry.Cartesian123Point{FT}

end # module


text:
# This file is part of Maxima.jl. It is licensed under the MIT license
#   Copyright (c) 2016 Nathan Smith

export sum,
       integrate,
       risch,
       diff,
       limit,
       taylor,
       product,
       changevar,
       lbfgs,
       minimum,
       laplace,
       ilt

import Base: sum,
             minimum

function integrate(m::MExpr, s)
    MExpr("integrate($m, $s)") |> mcall
end

"""
    integrate(f::T, x) where T

Evaluate the indefinite integral
```math
\\int f(x) dx
```
## Examples

```julia
julia> integrate(:(sin(x)), :x)
:(-cos(x))
```

## See also

* [`risch`](@refs Maxima.risch)
"""
function integrate(expr::T, s) where T
    m = MExpr(expr)
    out = integrate(m, s)
    convert(T, out)
end

function integrate(m::MExpr, s, lower, upper)
    MExpr("integrate($m, $s, $lower, $upper)") |> mcall
end

"""
    integrate(f::T, x, a, b) where T

Evaluate the definite integral of ``f`` with respect to ``x`` from ``a`` to ``b``.

```math
    \\int_a^b f(x) dx
```
"""
function integrate(expr::T, s, lower, upper) where T
    m = MExpr(expr)
    out = integrate(m, s, lower, upper)
    convert(T, out)
end

function risch(m::MExpr, var)
    MExpr("risch($m, $var)") |> mcall
end

"""
    risch(f::T, x) where T

Compute the indefinite integral of ``f`` with respect to ``x`` using the Risch algorithm

## See also
* [`integrate`](@ref Maxima.integrate)
"""
function risch(exp::T, var) where T
    m = MExpr(exp)
    convert(T, risch(m, var))
end

function diff(m::MExpr, s)
    MExpr("diff($m, $s)") |> mcall
end

"""
    diff(f, x) where T

Take the derivative of ``f`` with respect to ``x``

```math
\\frac{\\partial f(x, ...)}{\\partial x}
```
"""
function diff(exp::T, s) where T
    m = MExpr(exp)
    out = diff(m, s)
    convert(T, out)
end

function diff(m::MExpr, s, order::Integer)
    if order < 0
        error("Order of derivative must be positive integer")
    end
    MExpr("diff($m, $s, $order)") |> mcall
end

"""
    diff(f, x, n) where T

Take the nth order derivative of ``f`` with respect to ``x``.

```math
\\frac{\\partial^n f(x,...)}{\\partial x^n}
```
"""
function diff(exp::T, s, order::Integer) where T
    m = MExpr(exp)
    out = diff(m, s, order)
    convert(T, out)
end

function limit(m::MExpr, x, a)
    MExpr("limit($m, $x, $a)") |> mcall
end

"""
    limit(f, x, a) where T

Take the limit as ``x`` approaches ``a`` of ``f(x)``

```math
lim_{x \\rightarrow a} f(x)
```
"""
function limit(exp::T, x, a) where T
    m = MExpr(exp)
    out = limit(m, x, a)
    convert(T, out)
end

"""
    limit(f, x, a, side) where T

Take the left or right sided limit as ``x`` approaches ``a`` of ``f(x)``.

`side` may be either `\"plus\"` or `\"minus\"` denoting the right and left sided limit respectively
"""
function limit(m::MExpr, x, a, side)
    if "$side" != "plus" && "$side" != "minus"
        error("Side of limit must be \"plus\" or \"minus\"")
    end
    MExpr("limit($m, $x, $a, $side)") |> mcall
end

function limit(exp::T, x, a, side) where T
    m = MExpr(exp)
    out = limit(m, x, a, side)
    convert(T, out)
end

function sum(m::MExpr, k, start, finish)
    MExpr("sum($m, $k, $start, $finish), simpsum") |> mcall
end

"""
    sum(f::T, k, start, finish) where T

Compute the sum,

```math
    \\sum_{k=start}^{finish} f(k),
```
simplifying the sum if possible.
## Examples
```julia
julia> sum(m"1/n^2", :n, 1, "inf")

                                        2
                                     %pi
                                     ----
                                      6

```
"""
function sum(exp::T, k, start, finish) where T
    sumexp = Meta.parse("sum($exp, $k, $start, $finish)")
    m = MExpr(sumexp)
    out = mcall(MExpr("$m, simpsum"))
    convert(T, out)
end

function taylor(m::MExpr, x, x0, order::Integer)
    if order < 0
        error("Order of taylor expansion must be ≥ 0")
    end
    MExpr("taylor($m, $x, $x0, $order)") |> mcall
end

"""
    taylor(f::T, x, x0, n) where T

Taylor expand ``f(x)`` around ``x_0`` to nth order
"""
function taylor(exp::T, x, x0, order::Integer) where T
    m = MExpr(exp)
    out = taylor(m, x, x0, order)
    convert(T, out)
end

function product(m::MExpr, k, start, finish)
    MExpr("product($m, $k, $start, $finish), simpsum") |> mcall
end

"""
    product(f::T, k, start, finish) where T

Compute the product,

```math
    \\prod_{k=start}^{finish} f(k),
```
simplifying the product if possible.
"""
function product(exp::T, k, start, finish) where T
    sumexp = parse("product($exp, $k, $start, $finish)")
    convert(T, sumexp |> MExpr |> mcall)
end

function changevar(integral::MExpr, transform::MExpr, new, old)
    "changevar($integral, $transform, $new, $old)" |> MExpr |> mcall
end

function lbfgs(func::MExpr, var, guess; xtol = 1e-8)
    temp = "lbfgs($func, [$var], [$guess], $xtol)" |> MExpr |> mcall
    "rhs($temp[1])" |> MExpr |> mcall
end

function minimum(func::MExpr, var, guess; xtol = 1e-8)
    lbfgs(func, var, guess; xtol = xtol)
end

function laplace(func::MExpr, oldvar, newvar)
    "laplace($func, $oldvar, $newvar)" |> MExpr |> mcall
end

"""
    laplace(f::T, t, s) where T

Compute the Laplace transform of ``f(t)`` where ``s`` is the new variable.

```math
\\mathcal{L}\\lbrace f \\rbrace (s) = \\int_0^\\infty f(t) e^{-st} dt
```
"""
function laplace(func::T, oldvar, newvar) where T
    m = MExpr(func)
    convert(T, "laplace($m, $oldvar, $newvar)" |> MExpr |> mcall)
end

function ilt(func::MExpr, oldvar, newvar)
    "ilt($func, $oldvar, $newvar" |> MExpr |> mcall
end

"""
    ilt(f::T, s, t) where T

Compute the inverse Laplace transform of ``f(s)``.

``t`` is the new variable and ``s`` is the old variable.
"""
function ilt(func::T, oldvar, newvar) where T
    m = MExpr(func)
    convert(T, "ilt($m, $oldvar, $newvar)" |> MExpr |> mcall)
end


text:
module RandomWalk

export randomwalk1d, rwtrap1d_lifetime, rwtrap1d_lifetime_exp

"""
1D Random Walk

# Arguments

-`time::Integer`: the number of steps of the random walk
-`p::Float64`: the probability of going right
-`walks::Integer`: the number of walks, defaults to `100*time`

# Returns

a vector containing the coordinates of each point and a vector counting how many times
the walker landed on each point ath the end of the walk (only non-zero points are returned)
"""
function randomwalk1d(time::Integer, p::Float64; walks::Integer=100*time)
    path = zeros(2 * time + 1)
    coordinates = collect(Integer, -time:time)

    for _ in 1:walks
        pos = time + 1
        for _ in 1:time
            if rand() <= p
                pos += 1
            else
                pos -= 1
            end
        end
        path[pos] += 1
    end

    nonzero = findall(x -> x != 0, path)
    return coordinates[nonzero], path[nonzero]
end

"""
Calculate the average lifetime of a 1D random walker inside trapping boundaries

# Arguments

-`len::Integer`: The length of the path. The coordinates are in the range [-len÷2, len÷2]
-`initpos::Integer`: the initial position of the walker
-`walks::Integer`: the number of walks, defaults to the value of `time`
-`maxtime::Integer`: The maximum number of steps before the walker is abandoned;
    in case the number of steps exceeds `maxtime`, the lifetime is taken to be `maxtime`.
    defaults to `100*len`.

# Returns

The average lifetime, which is the average time passed before the walker gets trapped
in the boundaries of the path
"""
function rwtrap1d_lifetime(len::Integer, initpos::Integer, walks::Integer;
        maxtime::Integer=100*len)
    lifetime_sum = 0

    for _ in 1:walks
        pos = initpos
        time = 0

        while time < maxtime
            time += 1
            if rand() <= 0.5
                pos += 1
                if pos > len ÷ 2
                    break
                end
            else
                pos -= 1
                if pos < -len ÷ 2
                    break
                end
            end
        end

        lifetime_sum += time
    end

    lifetime_avg = lifetime_sum / walks
    return lifetime_avg
end

"""
Calculate the expected lifetime of a 1D random walker inside trapping boundaries

# Arguments

-`len::Integer`: The length of the path. The coordinates are in the range [-len÷2, len÷2]
-`initpos::Integer`: the position of the boundaries (traps)
-`maxtime::Integer`: The maximum number of steps before the walker is abandoned;
    in case the number of steps exceeds `maxtime`, the lifetime is taken to be `maxtime`.

# Returns

The average lifetime, which is the average time passed before the walker gets trapped
in the boundaries of the path
"""
function rwtrap1d_lifetime_exp(len::Integer, initpos::Integer; maxtime::Integer=100*len)
    lifetime_avg = 0

    probs = zeros(len + 3)
    probs[initpos + len ÷ 2 + 2] = 1

    for _ in 1:maxtime
        newprobs = zeros(len + 3)
        for i in 2:len+2
            newprobs[i - 1] += probs[i] * 0.5
            newprobs[i + 1] += probs[i] * 0.5
        end
        probs = newprobs

        lifetime_avg += sum(probs[2:end-1])
    end

    return lifetime_avg
end

end

max_stars_repo_name: slhshamloo/comp-phys

max_issues_repo_name: slhshamloo/comp-phys

max_forks_repo_name: slhshamloo/comp-phys


text:
mutable struct BSONWriteBuffer <: DenseVector{UInt8}
    data::Vector{UInt8}
    len::Int

    BSONWriteBuffer() = new(UInt8[], 0)
end

function Base.sizehint!(buffer::BSONWriteBuffer, capacity::Integer)
    resize!(buffer.data, max(capacity, buffer.len))
    buffer
end

@inline function Base.empty!(buffer::BSONWriteBuffer)
    buffer.len = 0
    buffer
end

@inline Base.length(buffer::BSONWriteBuffer) = buffer.len

@inline Base.size(buffer::BSONWriteBuffer) = (buffer.len,)

@inline Base.sizeof(buffer::BSONWriteBuffer) = buffer.len

@inline function Base.getindex(buffer::BSONWriteBuffer, i::Integer)
    @boundscheck (i > 0 && i <= buffer.len) || throw(BoundsError(buffer, i))
    @inbounds buffer.data[i]
end

@inline function Base.resize!(buffer::BSONWriteBuffer, len::Int)
    if length(buffer.data) < len
        resize!(buffer.data, len * 2)
    end
    buffer.len = len
    buffer
end

@inline Base.pointer(buffer::BSONWriteBuffer) = pointer(buffer.data)

@inline function Base.push!(buffer::BSONWriteBuffer, x::UInt8)
    if length(buffer.data) == buffer.len
        resize!(buffer.data, max(buffer.len * 2, 10))
    end
    buffer.len += 1
    @inbounds buffer.data[buffer.len] = x
    buffer
end


text:
abstract type AbstractPokemonScene <: AbstractVsScene end

const PokemonContext{T} = VsContext{T, PokemonBattle} where T <: AbstractVsStrategy


text:
function align_block(file; delim=r"=", block_begin=r"@c?enum", block_end=r"^end$")
    lines = readlines(file)
    final_lines = String[]

    is_target_block = false
    block_lines = nothing

    for (i, l) ∈ enumerate(lines)

        if is_target_block
            if !isnothing(match(block_end, l))
                lengths = map(block_lines) do decl
                    first(findfirst(delim, decl)) - 1
                end

                nmax = maximum(lengths)
                j = nmax == first(lengths) ? 2 : 1

                block_lines[j] = rpad(first(block_lines[j], lengths[j]), nmax) * block_lines[j][lengths[j]+1:end]

                append!(final_lines, block_lines)
                push!(final_lines, l)

                is_target_block = false
            else
                push!(block_lines, l)
            end
        else
            push!(final_lines, l)
        end

        if !isnothing(match(block_begin, l))
            block_lines = String[]
            is_target_block = true
        end
    end

    open(file, "w") do io
        print(io, (join(final_lines, '\n')))
    end

    format_file(file, margin=92, align_pair_arrow=true, align_assignment=true, align_struct_field=true, align_conditional=true)
end

custom_print(io::IO, expr::Expr) = println(io, prettify(expr))

function format()
    inst_file = src_dir("generated", "instructions.jl")
    enum_file = src_dir("generated", "enums.jl")

    @info "Formatting."

    # initial formatting
    format_file(inst_file, margin=92)
    format_file(enum_file, margin=92)

    # align blocks
    align_block(inst_file)
    align_block(inst_file, block_begin=r"const class_printing = Dict\(", block_end=r"^\)$", delim=r"=>")
    align_block(inst_file, block_begin=r"const kind_to_category = Dict\(", block_end=r"^\)$", delim=r"=>")

    align_block(enum_file)

    true
end


text:
using NumericalMethodsforEngineers

a = [16. 4. 8.; 4. 5. -4.; 8. -4. 22.]
b = [4., 2., 5.]

f = cholesky(a)
f |> display

y = f.L \ b
c = f.U \ y
println("Solution Vector: \n $c")

@test c == a \ b

a = [16. 4. 8.; 4. 5. -4.; 8. -4. 22.]
b = [4., 2., 5.]

a[2, 2] += 1.0e20
b[2] = a[2, 2] * 5.0

F = cholesky(a)
upper = Matrix(sparse(F.U))
upper |> display

lower = F.L
lower |> display

y = lower \ b
c = upper \ y
println("Solution Vector: \n $c")

@test ≈(c, a \ b, atol=1e-5)


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

"""
    ymd2dn(date) -> number_of_days

### Purpose ###

Convert from a date to day of the year.

### Explanation ###

Returns the day of the year for `date` with January 1st being day 1.

### Arguments ###

* `date`: the date with `Date` type.  Can be a single date or an array of dates.

### Output ###

The day of the year for the given `date`.  If `date` is an array, returns an
array of days.

### Example ###

Find the days of the year for March 5 in the years 2015 and 2016 (this is a leap
year).

``` julia
ymd2dn([Date(2015, 3, 5), Date(2016, 3, 5)])
# => 2-element Array{Int64,1}:
#     64
#     65
```

### Note ###

`ydn2md` converts from year and day number of year to a date.
"""
const ymd2dn = Dates.dayofyear


text:
 # This is test of the code for computing a transit model integrated over a time step2
# Computes derivatives over the timestep.
using PyPlot

#include("../src/integrate_transit_cubature.jl")
include("../src/integrate_transit_simpson.jl")

# Test it out:

include("../test/loglinspace.jl")

t1 = -1.5; t2 = 1.5; nt = 10000; dt = 0.3
t = zeros(nt)
t .= linearspace(t1,t2,nt)
#t = [t[1235]] ; nt =1
# The following compares different tolerances and maxdepths:
r = 0.1; b0 = 0.5; u_n = [0.3,0.3]; nu = length(u_n)
favg0 = zeros(nt,5+nu)
favg1 = zeros(nt,5+nu)
neval1= zeros(Int64,nt)
favg2 = zeros(nt,5+nu)
neval2= zeros(Int64,nt)
favg3 = zeros(nt,5+nu)
neval3= zeros(Int64,nt)
trans = transit_init(r,b0,u_n,true)
param = [0.0,1.0,b0]   # [t_0,v,b_0]

# Now, carry out speed test:

function compute_lightcurve!(trans::Transit_Struct{T},param::Array{T,1},t::Array{T,1},favg0::Array{T,2},nt::Int64) where {T <: Real}
@inbounds for i=1:nt
  trans.b = sqrt(param[3]^2+(param[2]*(t[i]-param[1]))^2)
  favg0[i,1] =  transit_poly!(trans)
  favg0[i,2] =  trans.dfdrb[1]
  favg0[i,3] =  trans.dfdrb[2]/trans.b*param[2]^2*(param[1]-t[i])
  favg0[i,4] =  trans.dfdrb[2]*param[2]/trans.b*(t[i]-param[1])^2
  favg0[i,5] =  trans.dfdrb[2]*param[3]/trans.b
  favg0[i,6:5+trans.n] = trans.dfdu
end
return
end

compute_lightcurve!(trans,param,t,favg0,nt)
@time compute_lightcurve!(trans,param,t,favg0,nt)

function integrate_lightcurve!(trans::Transit_Struct{T},param::Array{T,1},t::Array{T,1},dt::T,favg1::Array{T,2},nt::Int64,tol::T,maxdepth::Int64,neval_t::Array{Int64,1}) where {T <: Real}
# First, find the points of contact:
b0 = abs(param[3]); r = trans.r; t0 = param[1]; v = param[2]
if b0 > (1.0+r)
  println("No transit")
  favg1[:,1] = one(T)
  favg1[:,2:5+trans.n] = zero(T)
  return
elseif (b0+r) > 1.0
  # Two points of contact:
  tc = [t0 - sqrt((1.0+r)^2-b0^2)/v,t0 + sqrt((1.0+r)^2-b0^2)/v]
  nc = 2
else
  # Four points of contact:
  tc = [t0 - sqrt((1.0+r)^2-b0^2)/v,t0 - sqrt((1.0-r)^2-b0^2)/v,t0 + sqrt((1.0-r)^2-b0^2)/v,t0 + sqrt((1.0+r)^2-b0^2)/v]
  nc = 4
end
dtinv = inv(dt)
fint  = zeros(T,6+trans.n)
@inbounds for i=1:nt
  neval_t[i] = 0
  t1 = t[i]-0.5*dt ; t2 = t[i]+0.5*dt
  if t2 < tc[1] || t1 > tc[nc]
    # No points lie within the transit:
    ftmp=zeros(T,6+trans.n)
    ftmp[1]=one(T)
  else
    # Loop over points of conjunction:
    tlim = [t1]
    for j=1:nc
      if t1 < tc[j] && tc[j] < t2
        push!(tlim,tc[j])
      end
    end
    push!(tlim,t2)
    ftmp=zeros(T,6+trans.n)
    for j=1:length(tlim)-1
#      neval_t[i] += integrate_timestep_gradient!(param,trans,tlim[j],tlim[j+1],tol*trans.r^2,maxdepth,fint)
       nevali,depthmaxi =  integrate_timestep_gradient!(param,trans,tlim[j],tlim[j+1],tol*trans.r^2,maxdepth,fint)
       neval_t[i] += nevali
      ftmp += fint
    end
    ftmp *= dtinv
  end
  favg1[i,1:5]=ftmp[1:5]
  # Convert from d_n to u_n derivatives:
  favg1[i,6:5+trans.n]=BLAS.gemv('T',1.0,trans.dgdu,ftmp[6:6+trans.n])
#  println("i: ",i," t: ",t[i]," result: ",ftmp)
end
return
end

integrate_lightcurve!(trans,param,t,dt,favg1,nt,1e-4,3,neval1)
@time integrate_lightcurve!(trans,param,t,dt,favg1,nt,1e-4,8,neval1)
@time integrate_lightcurve!(trans,param,t,dt,favg2,nt,1e-6,128,neval2)
@time integrate_lightcurve!(trans,param,t,dt,favg3,nt,1e-8,128,neval3)

#@time for i=1:nt
#        ftmp = integrate_timestep_gradient(param,trans,t[i],dt,1e-6*r^2,32)*dtinv
#        favg1[i,1:5]=ftmp[1:5]
#        # Convert from d_n to u_n derivatives:
#        favg1[i,6:5+nu]=BLAS.gemv('T',1.0,trans.dgdu,ftmp[6:6+nu])
#      end

# Now plot the results:
fig,axes = subplots(2,3,figsize=(12,8))
ax = axes[1]
ax[:plot](t,favg0[:,1],linestyle="--",linewidth=2,label=L"$\mathcal{F}$")
ax[:plot](t,favg1[:,1],linewidth=2,label=L"$\overline{\mathcal{F}}$")
ax[:plot]([-dt/2,-dt/2,-dt/2,dt/2,dt/2,dt/2],[0.9925,0.9935,0.993,0.993,0.9925,0.9935],color="b")
ax[:annotate](L"$\Delta t$",xy=[-0.08,0.991])
ax[:legend](loc="upper center")
ax[:axis]([-1.5,1.5,0.988,1.002])
ax[:set_xlabel]("Time")
ax[:set_ylabel]("Flux")
ax = axes[2]
ax[:plot](t,favg0[:,2]*100.,linestyle="--",linewidth=2,label=L"$\frac{\partial\mathcal{F}}{\partial r}$")
ax[:plot](t,favg1[:,2]*100.,linewidth=2,label=L"$\frac{\partial\overline{\mathcal{F}}}{\partial r}$")
ax[:legend](loc="upper center")
ax[:axis]([-1.5,1.5,-30,2])
ax[:set_xlabel]("Time")
ax[:set_ylabel]("Flux derivative [pct]")
ax = axes[3]
ax[:plot](t,favg0[:,3]*100.,linestyle="--",linewidth=2,label=L"$\frac{\partial\mathcal{F}}{\partial t_0}$")
ax[:plot](t,favg1[:,3]*100.,linewidth=2,label=L"$\frac{\partial\overline{\mathcal{F}}}{\partial t_0}$")
ax[:set_xlabel]("Time")
ax[:set_ylabel]("Flux derivative [pct]")
#ax[:axis]([-1.5,1.5,-1,6])
ax[:legend](loc="lower left")
ax = axes[4]
ax[:plot](t,favg0[:,4]*100.,linestyle="--",linewidth=2,label=L"$\frac{\partial\mathcal{F}}{\partial v}$")
ax[:plot](t,favg1[:,4]*100.,linewidth=2,label=L"$\frac{\partial\overline{\mathcal{F}}}{\partial v}$")
ax[:set_xlabel]("Time")
ax[:set_ylabel]("Flux derivative [pct]")
ax[:axis]([-1.5,1.5,-1,6])
ax[:legend](loc="upper center")
ax = axes[5]
ax[:plot](t,favg0[:,5]*100.,linestyle="--",linewidth=2,label=L"$\frac{\partial\mathcal{F}}{\partial b_0}$")
ax[:plot](t,favg1[:,5]*100.,linewidth=2,label=L"$\frac{\partial\overline{\mathcal{F}}}{\partial b_0}$")
ax[:set_xlabel]("Time")
ax[:set_ylabel]("Flux derivative [pct]")
ax[:legend](loc="upper center")
ax = axes[6]
ax[:plot](t,favg0[:,6]*1000.,linestyle="--",linewidth=2,label=L"$\frac{\partial\mathcal{F}}{\partial u_1}$")
ax[:plot](t,favg1[:,6]*1000.,linewidth=2,label=L"$\frac{\partial\overline{\mathcal{F}}}{\partial u_1}$")
ax[:plot](t,favg0[:,7]*1000.,linestyle="--",linewidth=2,label=L"$\frac{\partial\mathcal{F}}{\partial u_2}$")
ax[:plot](t,favg1[:,7]*1000.,linewidth=2,label=L"$\frac{\partial\overline{\mathcal{F}}}{\partial u_2}$")
ax[:set_xlabel]("Time")
ax[:set_ylabel]("Flux derivative [ppt]")
ax[:legend](loc="upper center")
#b = sqrt.(param[3]^2+(param[2]*(t-param[1])).^2)
#plot(dt^2*(d2fdb2.*(param[2]*(t-param[1])./b).^2+dfdb.*param[2]^2./b)/24.0)
savefig("integrate_transit_gradient.pdf", bbox_inches="tight")

read(STDIN,Char)

for i=1:7
#  clf(); plot(t,favg1[:,i]); plot(t,favg2[:,i]); plot(t,favg3[:,i]); read(STDIN,Char); clf(); plot(t,abs.(favg1[:,i]-favg3[:,i]));semilogy(t,abs.(favg2[:,i]-favg3[:,i])); read(STDIN,Char)
  clf(); plot(t,favg1[:,i]); plot(t,favg2[:,i]); plot(t,favg3[:,i]); read(STDIN,Char); clf(); plot(t,favg1[:,i]-favg3[:,i]);plot(t,favg2[:,i]-favg3[:,i]); read(STDIN,Char)
end


text:
#######################################################################
#
# An example of creating an Excel charts with a date axis using
# XlsxWriter.
#
# Original Python Copyright 2013-2016, John McNamara, jmcnamara@cpan.org
# https://github.com/jmcnamara/XlsxWriter

using Dates
using XlsxWriter

function test()

	wb = Workbook("chart_date_axis.xlsx")
	ws = add_worksheet!(wb)

	# Add a format for the headings.

	chart = add_chart!(wb, Dict("type"=> "line"))
	date_format = add_format!(wb, Dict("num_format"=> "dd/mm/yyyy"))

	# Widen the first column to display the dates.
	set_column!(ws, "A:A", 12)

	# Some data to be plotted in the worksheet.
	dates = [Date(2013, 1, 1),
			 Date(2013, 1, 2),
			 Date(2013, 1, 3),
			 Date(2013, 1, 4),
			 Date(2013, 1, 5),
			 Date(2013, 1, 6),
			 Date(2013, 1, 7),
			 Date(2013, 1, 8),
			 Date(2013, 1, 9),
			 Date(2013, 1, 10)]

	values = [10, 30, 20, 40, 20, 60, 50, 40, 30, 30]

	# Write the date to the worksheet.
	write_column!(ws, "A1", dates, fmt=date_format)
	write_column!(ws, "B1", values)

	# Add a series to the chart.
	add_series!(chart, Dict(
		"categories"=> "=Sheet1!\$A\$1:\$A\$10",
		"values"=> "=Sheet1!\$B\$1:\$B\$10",
	))

	# Configure the X axis as a Date axis and set the max and min limits.
	set_x_axis!(chart, Dict(
		"date_axis"=> true,
		"min"=> Date(2013, 1, 2),
		"max"=> Date(2013, 1, 9),
	))

	# Turn off the legend.
	set_legend!(chart, Dict("none"=> true))

	# Insert the chart into the worksheet.
	insert_chart!(ws, "D2", chart)

	close(wb)
	isfile("chart_date_axis.xlsx")
end

test()


text:
"""
    AsTable(cols)

A type used for selection operations to signal that the columns selected by the
wrapped selector should be passed as a `NamedTuple` to the function.
"""
struct AsTable
    cols
end

Base.broadcastable(x::AsTable) = Ref(x)

function make_unique!(names::Vector{Symbol}, src::AbstractVector{Symbol};
                      makeunique::Bool=false)
    if length(names) != length(src)
        throw(DimensionMismatch("Length of src doesn't match length of names."))
    end
    seen = Set{Symbol}()
    dups = Int[]
    for i in 1:length(names)
        name = src[i]
        if in(name, seen)
            push!(dups, i)
        else
            names[i] = src[i]
            push!(seen, name)
        end
    end

    if length(dups) > 0
        if !makeunique
            dupstr = join(string.(':', unique(src[dups])), ", ", " and ")
            msg = "Duplicate variable names: $dupstr. Pass makeunique=true " *
                  "to make them unique using a suffix automatically."
            throw(ArgumentError(msg))
        end
    end

    for i in dups
        nm = src[i]
        k = 1
        while true
            newnm = Symbol("$(nm)_$k")
            if !in(newnm, seen)
                names[i] = newnm
                push!(seen, newnm)
                break
            end
            k += 1
        end
    end

    return names
end

function make_unique(names::AbstractVector{Symbol}; makeunique::Bool=false)
    make_unique!(similar(names), names, makeunique=makeunique)
end

"""
    gennames(n::Integer)

Generate standardized names for columns of a DataFrame.
The first name will be `:x1`, the second `:x2`, etc.
"""
function gennames(n::Integer)
    res = Vector{Symbol}(undef, n)
    for i in 1:n
        res[i] = Symbol(@sprintf "x%d" i)
    end
    return res
end

function funname(f)
    n = nameof(f)
    String(n)[1] == '#' ? :function : n
end

if isdefined(Base, :ComposedFunction) # Julia >= 1.6.0-DEV.85
    using Base: ComposedFunction
else
    using Compat: ComposedFunction
end

funname(c::ComposedFunction) = Symbol(funname(c.outer), :_, funname(c.inner))


text:
# example code for image manipulation using GtkUtilties

using  Gtk.ShortNames, Images, GtkUtilities, Graphics, TestImages

A = testimage("lighthouse")

win = Window("image")                                                 # Create Window
c = Canvas()                                                          # Create Canvas
push!(win, c)                                                         # Place Canvas in Window

yx_ranges = map(x->(first(x),last(x)), indices(A))                    # returns: yx_ranges = ((min-Y, max-Y), (min-X, max-X)). Note that y is first, since it refers to the rows of A
panzoom(c, yx_ranges[2], yx_ranges[1])                                # Initialize panzoom, with arguments for x-y yx_ranges.
panzoom_mouse(c)                                                      # Initialize mouse functions for panning/zooming

@guarded draw(c) do widget
    xview, yview = guidata[widget, :xview], guidata[widget, :yview]   #
    xv = round(Int, xview.min):round(Int, xview.max)                  #
    yv = round(Int, yview.min):round(Int, yview.max)                  #
    region_of_interest = view(A, yv, xv)                              # select subset of A according to the indices yv and xv
    copy!(widget, region_of_interest)                                 # copy the region of interest into c (widget)

    set_coordinates( getgc(widget),                                        # set_coordinates "associates" the two input intervals with the canvas "widget" (in this case, allows double zooming)
                BoundingBox(0,width(widget),                          # BoundingBox describing total viewable area in canvas (changes if you scale)
                            0,height(widget)),                        #
                BoundingBox(xview.min, xview.max,                     # BoundingBox describing current image size
                            yview.min, yview.max) )                   #

end

showall(win)


text:
function execute(config::Config, cmd::Cmd; dir::AbstractString="", name::AbstractString="")
    init(config)

    pg = PlaygroundConfig(config, dir, name)
    set_envs(pg)

    run(cmd)
end


text:
# Partial Least Squares (PLS1 and PLS2 NIPALS version)
module PLSRegressor

using JLD

include("utils.jl")
include("types.jl")
include("pls1.jl")
include("pls2.jl")
include("kpls.jl")
include("method.jl")

dir(path...) = joinpath(dirname(dirname(@__FILE__)),path...)

end


text:
#lowerBoundsUnimodal

#Based on Theorem 3
function _lower_bound_unimodal(D, M)
    if D/M <= 1/3
        return 1 / (1 - D/M)
    end

    return 8D / M / (1 + D/M )^2
end

####
## Separators for the Unimodal Case
####
#min_{t in [l, u]} a * H(t) + bt
#assumes h is convex, H(t) = 1/(mode-t) Int_t^mode h
function _sep_LB_Unimodal_1(h, H, a, b, mode, l, u)
	#check endpoints.
	val_star, t_star = a * H(l) + b * l, l

	if a * H(u) + b * u < val_star
		val_star = a * H(u) + b * u
		t_star = u
	end

	if a > 0 #convex case
		#dirty hack to get around t == mode
		function deriv(t)
			step = 1e-8
			if t == mode  #analytically = h_prime(mode)/ 2
				return a/4/step * (h(mode + step) - h(mode - step)) + b
			end
			a/(mode - t) * (H(t) - h(t)) + b
		end
		val_crit, t_crit = Inf, l
		if  deriv(l) * deriv(u) < 0 #critical point in interval
			t_crit = find_zero(deriv, (l,u), Bisection())
			val_crit = a * H(t_crit) + b * t_crit
		end
		if val_crit < val_star
			t_star = t_crit
			val_star = val_crit
		end
	end

	return val_star, t_star
end

#min_{t in [l, u]} a t^2 + b t + c int_[t, mode] h
#specifically for MAD:  h(t) = M/2 * abs(t-1) - D
#Notice this is specialized ot MAD bc otherwise optimization might be tricky
function _sep_LB_Unimodal_MAD_2(M, D, a, b, c, mode, l, u)
	h(t) = M/2 * abs(t-1) - D
	H_un(t) = quadgk(h, t, mode)[1]

	#check endpoints first
	t_star = l
	val_star = a * t_star^2 + b * t_star + c * H_un(t_star)

	val_u = a * u^2 + b * u + c * H_un(u)
	if val_u < val_star
		val_star = val_u
		t_star = u
	end

	#the <=1 critical point
	t_less = (-2b - 2c * D + c * M) / (4a + c * M)
	val_less = +Inf
	if t_less <=1 && l <= t_less <= u
		val_less = a * t_less^2 + b * t_less + c * H_un(t_less)
	end
	if val_less < val_star
		val_star = val_less
		t_star = t_less
	end

	#the >= 1 critcical point
	t_greater = (-2b  - 2c * D - c * M) / (4a - c * M)
	val_greater = Inf
	if t_greater >= 1 && l <= t_greater <= u
		val_greater = a * t_greater^2 + b * t_greater + c * H_un(t_greater)
	end
	if val_greater < val_star
		val_star = val_greater
		t_star = t_greater
	end

	return val_star, t_star
end

#min_{t in [l, u]} a t^2 + b t + c int_[t, mode] h
#specifically for CV:  h(t) = M^2 (t-1)^2 - C^2
#Notice this is specialized to CV bc otherwise optimization might be tricky
function _sep_LB_Unimodal_CV_2(M, C, a, b, c, mode, l, u)
	#define the functions
	h(t) = M^2 * (t-1)^2 - C^2
	H_unc(t) = quadgk(h, t, mode)[1]
	f(t) = a * t^2 + b*t + c * H_unc(t)

	#check end points
	t_star = l
	val_star = f(t_star)

	if f(u) < val_star
		t_star = u
		val_star = f(u)
	end

	#compute critical points and check feasibility
	disc = a^2 + c * (-2a + b - C^2) * M^2
	if disc >= 0
		t1 = a - c*M^2 - sqrt(disc)
		t1 /= c * M^2
		if l <= t1 <= u && f(t1) < val_star
			val_star = f(t1)
			t_star = t1
		end

		t2 = a - c*M^2 + sqrt(disc)
		t2 /= c*M^2
		if l <= t2 <= u && f(t2) < val_star
			val_star = f(t2)
			t_star = t2
		end
	end

	return val_star, t_star
end

#min_{t in [l, u]} a t^2 + b t + c int_[t, mode] h
#specifically for GM:  h(t) = -log(M ( t-1) + 1) + log(Bbar)
#Notice this is specialized to GM bc otherwise optimization might be tricky
#Bbar  = B/mu
function _sep_LB_Unimodal_GM_2(M, Bbar, a, b, c, mode, l, u)
	#define the functions
	h(t) = -log(M * (t - 1) + 1 ) + log(Bbar)
	H_unc(t) = quadgk(h, t, mode)[1]
	f(t) = a * t^2 + b * t + c * H_unc(t)

	#check end points
	t_star, val_star = l, f(l)

	if f(u) < val_star
		t_star, val_star = u, f(u)
	end

	#compute critical points and check feasibility
	#Degenerate case
	if c == 0.
		tp = -b /2/a
		if l <= tp <= u && f(tp) < val_star
			t_star, val_star = tp, f(tp)
		end
		return val_star, t_star
	end

	#else real work to do
	inner = 2a/c/M * Bbar * exp(2a/c * (1/M - 1) - b/c)

	#no interior solution
	if inner < -1/MathConstants.e
		return val_star, t_star
	end

	if inner < 0
		#check the -1 branch
		tp = 1 - 1/M + c/2/a * lambertw(inner)
		if l <= tp <= u && f(tp) < val_star
			val_star, t_star = f(tp), tp
		end
	end
	#always check the principal branch
	tp = 1 - 1/M + c/2/a * lambertw0(inner)
	if l <= tp <= u && f(tp) < val_star
		val_star, t_star = f(tp), tp
	end

	return val_star, t_star
end

#####################
#Convenience call that integrates h and calls workhorse
#INtegration routines can introduce numerical error. Better ot use other method when integrals known analytically
function _vopp_lb_unimodal_(S, mode, h, delta, sep_1, sep_2, numCuts, TOL, print_trace, ps)
	#Integrated functions
	function H_un(t)
		if min(mode, t) <= 1 <= max(mode, t)
			return quadgk(h, t, 1, mode)[1]
		else
			return quadgk(h, t, mode)[1]
		end
	end
	H(t) = t == mode ? h(mode) : H_un(t) / (mode - t)
	return _vopp_lb_unimodal_(S, mode, h, H, H_un, delta, sep_1, sep_2, numCuts, TOL, print_trace, ps)
end

#computes the maximum over the geometric grid
#sep_1(a, b, l, u) solves  min_{t in [l, u]}  a H(t) + b t
#sep_2(a, b, c, l, u) solves min_{t in [l, u]} a t^2 + b t + c int_{t,m} h
#user must ensure sep_1, sep_2 and h are consistent.
#both return val, tstar
#Note that H is normalized in sep 1 and the integral is UNnormalized in sep2
function _vopp_lb_unimodal_(S, mode, h, H, H_un, delta, sep_1, sep_2, numCuts, TOL, print_trace, ps)
	rstar = -Inf
	pstar = 0.
	rj = 0.

	for pj in ps
		#dispatch to appropriate case
		print_trace && println("\n Price Pj:\t", pj)
		if pj < mode
			rj = _rev_p_lb_unimodal_case1(S, mode, h, sep_1, sep_2, H, H_un, pj, numCuts, print_trace, TOL)
			print_trace && println("\t \t", "Rj:\t", rj, "\n")
		elseif pj == mode
			rj = _rev_p_lb_unimodal_case2(S, mode, h, sep_1, sep_2, H, H_un, pj, numCuts, print_trace, TOL)
			print_trace && println("\t \t", "Rj:\t", rj, "\n")
		else
			rj = _rev_p_lb_unimodal_case3(S, mode, h, sep_1, sep_2, H, H_un, pj, numCuts, print_trace, TOL)
			print_trace && println("\t \t", "Rj:\t", rj, "\n")
		end
		if rj > rstar
			rstar = rj
			pstar = pj
		end
	end
	return 1/rstar, pstar
end

#called by _vopp_lb_unimodal_
function _rev_p_lb_unimodal_case1(S, mode, h, sep_1, sep_2, H, H_un, pj, numCuts, print_trace, TOL)
	@assert pj < mode "Wrong Case called for LB Unimodal: price point: $pj mode: $mode"

	m = Model(solver=GurobiSolver(OutputFlag=false))

	@variable(m, theta)
	@variable(m, lam1)
	@variable(m, lam2)
	@objective(m, Min, theta + (2 - mode) * lam2)

	#add endpoints to get started
	#t = 0
	@constraint(m, theta * mode  + lam1 * H_un(0) >= pj * (mode - pj))
	#t = pj
	@constraint(m, theta + lam1 * H(pj) + lam2 * pj >= pj)
	#t = S
	@constraint(m, theta + lam1 * H(S) + lam2 * S >= pj)

	#Add a silly constraint to prevent unboundedness
	@constraint(m, theta + (2 - mode) * lam2 >= -100)

	for iter = 1:numCuts
		status = solve(m)
		feasible = true

		print_trace && println("\n   Iteration: \t", iter, "\t Status:\t", status)

		#1st constraint corresponds to second separator
		a = -getvalue(lam2)
		b = getvalue(lam2) * mode - getvalue(theta)
		c = getvalue(lam1)
		rhs = pj * (mode - pj) - getvalue(theta) * mode

		val, tstar = sep_2(a, b, c, 0, pj)
		print_trace && println("\t Iter: $iter  \t Sep2: LHS-RHS >=? 0:  $(val-rhs) : tstar $tstar")

		if val < rhs - TOL
			@constraint(m, theta * (mode - tstar) +
							lam1 * H_un(tstar) +
							lam2 * tstar * (mode - tstar) >= pj * (mode - pj)
							)
			feasible = false
		end

		#2nd constraint corresponds to first separator
		val, tstar = sep_1(getvalue(lam1), getvalue(lam2), pj, S)
		rhs = pj - getvalue(theta)
		print_trace && println("\t Iter: $iter  \t Sep1: LHS-RHS >=? 0:  $(val-rhs) : tstar $tstar")

		if val < rhs - TOL
			###DEBUG
			###VG Confirm that the constraint is actually violated right now?
			temp_lhs = getvalue(lam1) * H(tstar) + getvalue(lam2) * tstar
			temp_rhs = pj - getvalue(theta)
			print_trace && println("Violation check: tstar $tstar \n",
					"LHS >=? RHS: \t", temp_lhs, "\t", temp_rhs, "\t", temp_lhs - temp_rhs)

			print_trace && println("Violation Check 2\n",
				"val:\t $val", "temp_lhs:\t", temp_lhs)
			### END DEBUG

			cnst = @constraint(m, theta + lam1 * H(tstar) + lam2 * tstar >= pj)
			feasible = false

		end

		if feasible
			if status != :Optimal
				println("Problem likely infeasible.  Dual Status:\t", status)
				return -1.0
			else #correctly solved
				return getobjectivevalue(m)
			end
		end
	end
	##Exhausted iteration count
	print_trace && println("Max Iteration Count Reached. \n",
							"theta:\t", getvalue(theta),
							"\t lam1:\t", getvalue(lam1),
							"\t lam2:\t", getvalue(lam2))
	throw("Max Iteration Count Reached")
	return -1.
end

function _rev_p_lb_unimodal_case2(S, mode, h, sep_1, sep_2, H, H_un, pj, numCuts, print_trace, TOL)
	@assert abs(pj - mode)<= 1e-10 "Wrong Case called for LB Unimodal: price point: $pj mode: $mode"

	m = Model(solver=GurobiSolver(OutputFlag=false))
	@variable(m, theta)
	@variable(m, lam1)
	@variable(m, lam2)
	@objective(m, Min, theta + (2 - mode) * lam2)

	#add endpoints to get started
	#t = 0
	@constraint(m, theta + lam1 * H(0) >= 0)
	#t = mode
	@constraint(m, theta + lam1 * H(mode) + lam2 * mode >= mode)
	#t = S
	@constraint(m, theta + lam1 * H(S) + lam2 * S >= mode)

	#A trivial (bad) constraint for unboundedness?
	@constraint(m, theta + (2-mode) * lam2 >= -100)

	for iter = 1:numCuts
		status = solve(m)
		feasible = true

		print_trace && println("\n   Iteration: \t", iter, "Status:\t", status)
		#1st constraint corresponds to first opt type
		val, tstar = sep_1(getvalue(lam1), getvalue(lam2), 0, mode)
		rhs = -1 * getvalue(theta)

		print_trace && println("\t Iter: $iter  \t Sep1: LHS -RHS >=? 0:  $(val-rhs) :\t tstar $tstar")

		if val < rhs - TOL
			@constraint(m, theta + lam1 * H(tstar) + lam2 * tstar >= 0)
			feasible = false
		end

		rhs = mode - getvalue(theta)
		val, tstar = sep_1(getvalue(lam1), getvalue(lam2), mode, S)

		print_trace && println("\t Iter: $iter \t Sep1: LHS - RHS ?>= 0:  $(val-rhs) :\t tstar $tstar")

		if val < rhs - TOL
			@constraint(m, theta + lam1 * H(tstar) + lam2 * tstar >= mode)
			feasible = false
		end

		if feasible
			if status != :Optimal
				println("Problem likely infeasible.  Dual Status:\t", status)
				return -1.0
			else #correctly solved

				###DEBUG
				if getobjectivevalue(m) < 0
					print_trace && println("Weird Objective\n",
						"theta:\t", getvalue(theta),
						"\t lam1:\t", getvalue(lam1),
						"\t lam2:\t", getvalue(lam2))
				end

				return getobjectivevalue(m)
			end
		end
	end
	##Exhausted iteration count
	print_trace && println("Max Iteration Count Reached. \n",
							"theta:\t", getvalue(theta),
							"\t lam1:\t", getvalue(lam1),
							"\t lam2:\t", getvalue(lam2))
	throw("Max Iteration Count Reached")
	return -1.
end

function _rev_p_lb_unimodal_case3(S, mode, h, sep_1, sep_2, H, H_un, pj, numCuts, print_trace, TOL)
	@assert pj > mode "Wrong Case called for LB Unimodal: price point: $pj mode: $mode"

	m = Model(solver=GurobiSolver(OutputFlag=false))
	@variable(m, theta)
	@variable(m, lam1)
	@variable(m, lam2)

	@objective(m, Min, theta + (2 - mode) * lam2)

	#add endpoints to get started
	#t = 0
	@constraint(m, theta  + lam1 * H(0) >= 0)
	#t = pj
	@constraint(m, theta + lam1 * H(pj) + lam2 * pj >= 0)
	#t = S
	@constraint(m, theta * (S - mode) - lam1 * H_un(S) + lam2 * S * (S - mode) >= pj * (S - pj))

	#VG Try to prevent unbounded errors.
	@constraint(m, theta + (2 - mode) * lam2 >= -10.)

	for iter = 1:numCuts
		status = solve(m)
		feasible = true

		print_trace && println("\nIteration: \t", iter, "\t Status: ", status)
		#1st constraint corresponds to first separator
		val, tstar = sep_1(getvalue(lam1), getvalue(lam2), 0, pj)
		rhs = -getvalue(theta)
		print_trace && println("\t Iter: $iter Sep1 val -rhs ?>= 0: ", val-rhs, " : ", tstar)

		if val < rhs - TOL
			cnst = @constraint(m, theta + lam1 * H(tstar) + lam2 * tstar >= 0)
			feasible = false
		end

		#2nd constraint corresponds to second separator
		a = getvalue(lam2)
		b = getvalue(theta) - getvalue(lam2) * mode - pj
		c = -getvalue(lam1)
		rhs = getvalue(theta) * mode - pj^2

		val, tstar = sep_2(a, b, c, pj, S)
		print_trace && println("\t Iter: $iter Sep 2 val-rhs ?>= 0: ", val-rhs, " : ", tstar)

		if val < rhs - TOL
			# print_trace && println("\t VG Extra Detailed Debug")
			# print_trace && println("\t", getvalue(theta), "\t", getvalue(lam1), "\t", getvalue(lam2) )
			# print_trace && println("\t tstar:\t", tstar, "\t H_un", H_un(tstar), "\t mode\t", mode)

			@constraint(m, theta * (tstar - mode) -
							lam1 * H_un(tstar) +
							lam2 * tstar * (tstar - mode) >= pj * (tstar - pj)
							)
			feasible = false
		end

		if feasible
			if status != :Optimal
				println("Problem likely infeasible.  Dual Status:\t", status)
				return -1.0
			else #correctly solved
				return getobjectivevalue(m)
			end
		end
	end
	##Exhausted iteration count
	throw("Max Iteration Count Reached")
	return -1.
end

#safe_fail forces a return of -1 if D seems infeasible
function vopp_lb_unimodal_MAD(mu, S, M, D, mode;
			method=:Opt, delta = S/100, numCuts=100, print_trace=false, TOL=1e-6, safe_fail=false, pj = -1)
	#standardize
	c = mu * (1 - M)
	Sc = vopp.comp_Sc(S, M)
	mode_c = (mode - c) / (mu - c)

	if D/M > .25
		println("Deviation $(D/M) > .25, problem likely infeasible")
		if safe_fail
			println("Safe_fail set to true. Aborting Calc")
			return -1.0
		end
	end

	if method == :Formula
		return _lower_bound_unimodal(D, M)
	end
	if method != :Opt
		throw("Method must be one of :Opt or :Formula")
	end

	#create the separators
	h(t) = M/2 * abs(t-1) - D
	H(t) = mode_c == t ? h(mode_c) : quadgk(h, t, mode_c)[1] / (mode_c - t)
	sep_1(a, b, l, u) = _sep_LB_Unimodal_1(h, H, a, b, mode_c, l, u)
	sep_2(a, b, c, l, u) = _sep_LB_Unimodal_MAD_2(M, D, a, b, c, mode_c, l, u)

	if pj > 0
		ps = [ (pj - c)/(mu - c) ]
	else
		#search over the a geometric price ladder
		ps = geom_price_ladder(S, delta)
	end

	#kick it to the workhorse
	_vopp_lb_unimodal_(Sc, mode_c, h, delta, sep_1, sep_2, numCuts, TOL, print_trace, ps)
end

function vopp_lb_unimodal_CV(mu, S, M, C, mode; delta = S/100, numCuts=100, print_trace=false, TOL=1e-6, pj = -1)
	#standardize
	c = mu * (1 - M)
	Sc = vopp.comp_Sc(S, M)
	mode_c = (mode - c) / (mu - c)

	if C > max_cv(S, M)
		println("Warning:  Desired Coeff of Variation might exceed max-possible for uniform distributions.")
	end

	if pj > 0
		ps = [ (pj - c)/(mu - c) ]
	else
		#search over the a geometric price ladder
		ps = geom_price_ladder(S, delta)
	end

	#create the separators
	h(t) = M^2 * (t - 1)^2 - C^2
	H(t) = mode_c == t ? h(mode_c) : quadgk(h, t, mode_c)[1] / (mode_c - t)
	sep_1(a, b, l, u) = _sep_LB_Unimodal_1(h, H, a, b, mode_c, l, u)
	sep_2(a, b, c, l, u) = _sep_LB_Unimodal_CV_2(M, C, a, b, c, mode_c, l, u)

	#kick it to the workhorse
	_vopp_lb_unimodal_(Sc, mode_c, h, delta, sep_1, sep_2, numCuts, TOL, print_trace, ps)
end

function vopp_lb_unimodal_GM(mu, S, M, B, mode; delta = S/100, numCuts=100, print_trace=false, TOL=1e-6, pj = -1)
	#standardize
	c = mu * (1 - M)
	Sc = vopp.comp_Sc(S, M)
	mode_c = (mode - c) / (mu - c)
	Bbar = B/mu

	if Bbar > 1 || B < min_gm(S, M, mu)
		println("Warning:  Desired Geometric Mean might be infeasible.  $B not in [$(min_gm(S, M, mu)), 1]")
	end

	if pj > 0
		ps = [ (pj - c)/(mu - c) ]
	else
		#search over the a geometric price ladder
		ps = geom_price_ladder(S, delta)
	end

	#create the separators
	h(t) = -log(M * (t - 1) + 1) + log(Bbar)
	#notice scaled
	H(t) = mode_c == t ? h(mode_c) : quadgk(h, t, mode_c)[1] / (mode_c - t)

	#First separator works because of convexity
	sep_1(a, b, l, u) = _sep_LB_Unimodal_1(h, H, a, b, mode_c, l, u)
	sep_2(a, b, c, l, u) = _sep_LB_Unimodal_GM_2(M, Bbar, a, b, c, mode_c, l, u)

	#kick it to the workhorse
	_vopp_lb_unimodal_(Sc, mode_c, h, delta, sep_1, sep_2, numCuts, TOL, print_trace, ps)
end

function vopp_lb_unimodal_IC(mu, S, M, phat, q, mode; delta = S/100, numCuts=100, print_trace=false, TOL=1e-6, pj = -1)
	#standardize
	c = mu * (1 - M)
	Sc = vopp.comp_Sc(S, M)
	mode_c = (mode - c) / (mu - c)
	v0 = (phat - 1)/ M + 1

	@assert c <= phat <= S "Price outside valuation bounds:\t $phat not in [$c , $S ]"
	@assert 0 <= q <= 1 "Market fraction not $q not in [0, 1]"

	if pj > 0
		ps = [(pj - c)/(mu - c)]
	else
		#search over the a geometric price ladder
		ps = geom_price_ladder(S, delta)
		ps = [ps; range(max(0, phat/2), stop=min(2 * phat, 1), length=100)]
		push!(ps, v0)
		push!(ps, mode_c)
		sort!(ps)
	end

	#create the separators
	h(t) = t >= v0 ? 1 - q : -q
	#notice scaled
	function H(t)  #G(v_0, mode_c, t) - q
    	if v0 > max(mode_c, t)
        	return -q
    	elseif v0 < min(mode_c, t)
        	return 1 - q
    	elseif mode_c == t
        	return mode_c >= v0 ? 1 - q : -q
    	end
    	#non-trivial answer only if v0 in interval [mode_c, t]
    	return (max(mode_c, t) - v0) / abs(mode_c - t) - q
	end

	#sep_1(a, b, l, u) solves  min_{t in [l, u]}  a H(t) + b t
	function sep_1(a, b, l, u)
		f(t) = a * H(t) + b * t
		#check endpoints
		val, tstar = f(l), l
		if f(u) < val
			val, tstar = f(u), u
		end

		#check mode if its feasible
		if l <= mode_c <= u && f(mode_c) < val
			val, tstar = f(mode_c), mode_c
		end

		#check v0 if its feasible
		if l <= v0 <= u && f(v0) < val
			val, tstar = f(v0), v0
		end

		if a * b < 0
			if v0 > mode_c
				t1 = mode_c + sqrt( a/b * (mode_c - v0) )
				if l <= t1 <= u && f(t1) < val
					val, tstar = f(t1), t1
				end
			end
			if v0 < mode_c
				t4 = mode_c - sqrt( a/b * (v0 - mode_c) )
				if l <= t4 <= u && f(t4) < val
					val, tstar = f(t4), t4
				end
			end
		end  #end check for internal optimizers
		return val, tstar
	end  #ends def of separator

	#sep_2(a, b, c, l, u) solves min_{t in [l, u]} a t^2 + b t + c int_{t,m} h
	#Notice this reverses rolls of a, b from document

	H_un(t) = H(t) * (mode_c - t)
	function sep_2(a, b, c, l, u)
		obj(t) = a * t^2 + b * t + c * H_un(t)

		#check endpoints
		tstar, val = l, obj(l)
		if obj(u) < val
			tstar, val = u, obj(u)
		end

		#handle degenerate case where c == 0
		if c == 0.
			tp = -b/2/a
			if l <= tp <= u && obj(tp) < val
					val, tstar = obj(tp), tp
			end
			return val, tstar
		end

		#else real work to be done
		#check the two possible critical points
		if a != 0
			t1 = c/2/a * (1 - q - b/c)
			if t1 >= v0 && l <= t1 < u && obj(t1) < val
					val, tstar = obj(t1), t1
			end

			t0 = -c/2/a * (q + b/c)
			if t0 < v0 && l <= t0 < u && obj(t0) < val
					val, tstar = obj(t0), t0
			end
		end
		return val, tstar
	end  #end def of separator

	#kick it to the workhorse
	_vopp_lb_unimodal_(Sc, mode_c, h, H, H_un, delta, sep_1, sep_2, numCuts, TOL, print_trace, ps)
end