text:
mutable struct VoidUpdater <: Updater end
POMDPs.update(::VoidUpdater, ::B, ::Any, ::Any, b=nothing) where B = nothing

mutable struct RandomPolicy{P <: Union{MDP, POMDP}} <: Policy
    rng::AbstractRNG
    problem::P
end
RandomPolicy(problem::Union{POMDP,MDP};
            rng=Random.GLOBAL_RNG) = RandomPolicy(rng, problem)

function POMDPs.action(policy::RandomPolicy, s)
    return rand(policy.rng, actions(policy.problem, s))
end

mutable struct RandomSolver <: Solver
    rng::AbstractRNG
end

RandomSolver(;rng=Base.GLOBAL_RNG) = RandomSolver(rng)
POMDPs.solve(solver::RandomSolver, problem::P) where {P<:Union{POMDP,MDP}} = RandomPolicy(solver.rng, problem)

let
    rng = MersenneTwister(7)

    mdp = LegacyGridWorld()
    POMDPs.DDNStructure(::Type{typeof(mdp)}) = DDNStructure(MDP) |> add_infonode
    @test :info in nodenames(DDNStructure(mdp))
    s = initialstate(mdp, rng)
    a = rand(rng, actions(mdp))
    sp, r, i = @inferred gen(DDNOut(:sp,:r,:info), mdp, s, a, rng)
    @test i === nothing

    pomdp = TigerPOMDP()
    POMDPs.DDNStructure(::Type{typeof(pomdp)}) = DDNStructure(POMDP) |> add_infonode
    @test :info in nodenames(DDNStructure(pomdp))
    s = initialstate(pomdp, rng)
    a = rand(rng, actions(pomdp))
    sp, o, r, i = @inferred gen(DDNOut(:sp,:o,:r,:info), pomdp, s, a, rng)
    @test i === nothing

    up = VoidUpdater()
    policy = RandomPolicy(rng, pomdp)
    @inferred action_info(policy, s)

    solver = RandomSolver(rng=rng)
    policy, sinfo = solve_info(solver, pomdp)
    @test isa(sinfo, Nothing)

    d = initialstate_distribution(pomdp)
    b = initialize_belief(up, d)
    a = action(policy, b)
    sp, o, r = gen(DDNOut(:sp,:o,:r), pomdp, rand(rng, d), a, rng)
    @inferred update_info(up, b, a, o)
end


text:
# InternalOptions
mutable struct InternalOptions
    wrap_structs::Bool
    ismutable::Bool
end
InternalOptions() = InternalOptions(true, false)

"""
    WrapContext
Store shared information about the wrapping session.
"""
mutable struct WrapContext
    index::Index
    headers::Vector{String}
    output_file::String
    common_file::String
    clang_includes::Vector{String}               # clang include paths
    clang_args::Vector{String}                   # additional {"-Arg", "value"} pairs for clang
    header_library::Function                     # called to determine shared library for given header
    common_buf::OrderedDict{Symbol, ExprUnit}    # output buffer for common items: typedefs, enums, etc.
    output_bufs::DefaultOrderedDict{String, Array{Any}}
    options::InternalOptions
    exclude_symbols::Set{String}
    only_select_symbols::Set{String}
    fields_align::Dict{Tuple{Symbol,Symbol},Int} # {struct_symbol,field_symbol} -> alignment
    rewriter::Function
end

### Convenience function to initialize wrapping context with defaults
function init(; headers::Vector{String}                     = String[],
                index::Union{Index,Nothing}                 = nothing,
                output_dir::String                          = "",
                output_file::String                         = "",
                common_file::String                         = "",
                clang_args::Vector{String}                  = String[],
                clang_includes::Vector{String}              = String[],
                clang_diagnostics::Bool                     = true,
                header_library                              = nothing,
                options                                     = InternalOptions(),
                exclude_symbols::Set{String}                = Set{String}(),
                only_select_symbols::Set{String}            = Set{String}(),
                fields_align::Dict{Tuple{Symbol,Symbol},Int} = Dict{Tuple{Symbol,Symbol},Int}(),
                rewriter                                    = x -> x)

    # Set up some optional args if they are not explicitly passed.
    index == nothing && (index = Index(clang_diagnostics);)

    if output_file == ""
        output_file = joinpath(output_dir, strip(splitext(basename(x))[1]) * ".jl")
    end

    if common_file == ""
        common_file = joinpath(output_dir, common_file)
    end

    if header_library == nothing
        header_library = strip(splitext(basename(x))[1])
    elseif isa(header_library, String)
        libname = copy(header_library)
        header_library = x->libname
    end

    # Instantiate and return the WrapContext
    global context = WrapContext(index,
                                 headers,
                                 output_file,
                                 common_file,
                                 clang_includes,
                                 clang_args,
                                 header_library,
                                 OrderedDict{Symbol,ExprUnit}(),
                                 DefaultOrderedDict{String, Array{Any}}(()->Any[]),
                                 options,
                                 exclude_symbols,
                                 only_select_symbols,
                                 fields_align,
                                 rewriter)
    return context
end

function Base.run(wc::WrapContext, generate_template=true)
    # parse headers
    ctx = DefaultContext(wc.index, wc.exclude_symbols, wc.only_select_symbols, wc.fields_align)
    parse_headers!(ctx, wc.headers, args=wc.clang_args, includes=wc.clang_includes, index=wc.index)
    ctx.options["is_function_strictly_typed"] = false
    # Helper to store file handles
    filehandles = Dict{String,IOStream}()
    getfile(f) = (f in keys(filehandles)) ? filehandles[f] : (filehandles[f] = open(f, "w"))

    # First pass on the AST, add declarations to be translated to the queue.
    # Declarations to be translated include:
    # - macro definitions
    # - enumerations
    # - structures
    # - unions
    # - functions
    # and consists only in declarations from selected headers.
    # Types on which these declarations depend (e.g. parameters and return types of functions,
    # types of fields in structures) will be added to the queue during the traversal.
    for trans_unit in ctx.trans_units
        root_cursor = getcursor(trans_unit)
        header = spelling(root_cursor)

        root_children = children(root_cursor)
        for (i, child) in enumerate(root_children)

            # Save next cursor.
            ctx.next_cursor[child] = length(root_children) < i+1 ? nothing : root_children[i+1]

            child_kind = kind(child)
            child_name = split(name(child), "(")[1]
            child_header = filename(child)

            # Skip irrelevant cursors.
            if child in ctx.visited || # already added to the queue
                (is_forward_declaration(child) && child_kind != CXCursor_FunctionDecl && child_kind != CXCursor_MacroDefinition) || # forward type declaration
                startswith(child_name, "__") || # compiler definitions
                child_header != header || # cursors from other headers
                child_name in ctx.exclude_symbols ||
                (length(ctx.only_select_symbols) > 0 && !(child_name in ctx.only_select_symbols)) ||
                (child_kind != CXCursor_EnumDecl &&
                 child_kind !=  CXCursor_StructDecl &&
                 child_kind != CXCursor_UnionDecl &&
                 child_kind != CXCursor_EnumConstantDecl &&
                 child_kind != CXCursor_FunctionDecl &&
                 child_kind != CXCursor_MacroDefinition)
                continue
            end

            enqueue!(ctx.queue, child)
            push!(ctx.visited, child)
        end
    end

    # Wrap all cursors until the queue is empty.
    # For each declaration wrapped, other declarations on which it depends will be added to the queue.
    while !isempty(ctx.queue)
        cursor = dequeue!(ctx.queue)
        ctx.libname = wc.header_library(filename(cursor))
        try
            wrap!(ctx, cursor)
        catch err
            @error "error thrown."
            rethrow(err)
        end
    end

    # apply user-supplied transformation
    ctx.api_buffer = wc.rewriter(ctx.api_buffer)

    # write output
    out_file = wc.output_file
    @info "writing $(out_file)"

    out_stream = getfile(out_file)
    println(out_stream, "# Julia wrapper automatically generated using Clang.jl\n")
    print_buffer(out_stream, ctx.api_buffer)
    ctx.api_buffer = Expr[]  # clean api_buffer for the next header

    common_buf = dump_to_buffer(ctx.common_buffer)

    # apply user-supplied transformation
    common_buf = wc.rewriter(common_buf)

    # write "common" definitions: types, typealiases, etc.
    open(wc.common_file, "w") do f
        println(f, "# Automatically generated using Clang.jl\n")
        print_buffer(f, common_buf)
    end

    map(close, values(filehandles))

    if generate_template
        copydeps(dirname(wc.common_file))
        print_template(joinpath(dirname(wc.common_file), "LibTemplate.jl"))
    end
end

ext: jl


text:
# # Internal wave example
#
# In this example, we initialize an internal wave packet in two-dimensions
# and watch it propagate. This example illustrates how to set up a two-dimensional
# model, set initial conditions, and how to use `BackgroundField`s.
#
# ## Install dependencies
#
# First let's make sure we have all required packages installed.

using Pkg
pkg"add Oceananigans, JLD2, Plots"

# ## The physical domain
#
# First, we pick a resolution and domain size. We use a two-dimensional domain
# that's periodic in ``(x, y, z)``:

using Oceananigans

grid = RegularCartesianGrid(size=(128, 1, 128), x=(-π, π), y=(-π, π), z=(-π, π),
                            topology=(Periodic, Periodic, Periodic))

# ## Internal wave parameters
#
# Inertia-gravity waves propagate in fluids that are both _(i)_ rotating, and
# _(ii)_ density-stratified. We use Oceananigans' Coriolis abstraction
# to implement a background rotation rate:

coriolis = FPlane(f=0.2)

# On an `FPlane`, the domain is idealized as rotating at a constant rate with
# rotation period `2π/f`. `coriolis` is passed to `IncompressibleModel` below.
# Our units are arbitrary.

# We use Oceananigans' `background_fields` abstraction to define a background
# buoyancy field `B(z) = N^2 * z`, where `z` is the vertical coordinate
# and `N` is the "buoyancy frequency". This means that the modeled buoyancy field
# perturbs the basic state `B(z)`.

using Oceananigans.Fields: BackgroundField

## Background fields are functions of `x, y, z, t`, and optional parameters.
## Here we have one parameter, the buoyancy frequency
B_func(x, y, z, t, N) = N^2 * z

N = 1 ## buoyancy frequency

B = BackgroundField(B_func, parameters=N)

# We are now ready to instantiate our model. We pass `grid`, `coriolis`,
# and `B` to the `IncompressibleModel` constructor.
# We add a small amount of `IsotropicDiffusivity` to keep the model stable
# during time-stepping, and specify that we're using a single tracer called
# `b` that we identify as buoyancy by setting `buoyancy=BuoyancyTracer()`.

using Oceananigans.Advection

model = IncompressibleModel(
                 grid = grid,
            advection = CenteredFourthOrder(),
          timestepper = :RungeKutta3,
              closure = IsotropicDiffusivity(ν=1e-6, κ=1e-6),
             coriolis = coriolis,
              tracers = :b,
    background_fields = (b=B,), # `background_fields` is a `NamedTuple`
             buoyancy = BuoyancyTracer()
)

# ## A Gaussian wavepacket
#
# Next, we set up an initial condition that excites an internal wave that propates
# through our rotating, stratified fluid. This internal wave has the pressure field
#
# ```math
# p(x, y, z, t) = a(x, z) \, \cos(kx + mz - ω t) \, .
# ```
#
# where ``m`` is the vertical wavenumber, ``k`` is the horizontal wavenumber,
# ``ω`` is the wave frequncy, and ``a(x, z)`` is a Gaussian envelope.
# The internal wave dispersion relation links the wave numbers ``k`` and ``m``,
# the Coriolis parameter ``f``, and the buoyancy frequency ``N``:

## Non-dimensional internal wave parameters
m = 16      # vertical wavenumber
k = 8       # horizontal wavenumber
f = coriolis.f

## Dispersion relation for inertia-gravity waves
ω² = (N^2 * k^2 + f^2 * m^2) / (k^2 + m^2)

ω = sqrt(ω²)
nothing # hide

# We define a Gaussian envelope for the wave packet so that we can
# observe wave propagation.

## Some Gaussian parameters
A = 1e-9
δ = grid.Lx / 15

## A Gaussian envelope centered at ``(x, z) = (0, 0)``.
a(x, z) = A * exp( -( x^2 + z^2 ) / 2δ^2 )
nothing # hide

# An inertia-gravity wave is a linear solution to the Boussinesq equations.
# In order that our initial condition excites an inertia-gravity wave, we
# initialize the velocity and buoyancy perturbation fields to be consistent
# with the pressure field ``p = a \, \cos(kx + mx - ωt)`` at ``t=0``.
# These relations are sometimes called the "polarization
# relations". At ``t=0``, the polarization relations yield

u₀(x, y, z) = a(x, z) * k * ω   / (ω^2 - f^2) * cos(k*x + m*z)
v₀(x, y, z) = a(x, z) * k * f   / (ω^2 - f^2) * sin(k*x + m*z)
w₀(x, y, z) = a(x, z) * m * ω   / (ω^2 - N^2) * cos(k*x + m*z)
b₀(x, y, z) = a(x, z) * m * N^2 / (ω^2 - N^2) * sin(k*x + m*z)

set!(model, u=u₀, v=v₀, w=w₀, b=b₀)

# Recall that the buoyancy `b` is a perturbation, so that the total buoyancy field
# is ``N^2 z + b``.

# ## A wave packet on the loose
#
# We're ready to release the packet. We build a simulation with a constant time-step,

simulation = Simulation(model, Δt = 0.1 * 2π/ω, stop_iteration = 15)

# and add an output writer that saves the vertical velocity field every two iterations:

using Oceananigans.OutputWriters: JLD2OutputWriter, IterationInterval

simulation.output_writers[:velocities] = JLD2OutputWriter(model, model.velocities,
                                                          schedule = IterationInterval(1),
                                                            prefix = "internal_wave",
                                                             force = true)

# With initial conditions set and an output writer at the ready, we run the simulation

run!(simulation)

# ## Animating a propagating packet
#
# To visualize the solution, we load snapshots of the data and use it to make contour
# plots of vertical velocity.

using JLD2, Printf, Plots, Oceananigans.Grids

# We use coordinate arrays appropriate for the vertical velocity field,

x, y, z = nodes(model.velocities.w)
nothing # hide

# open the jld2 file with the data,

file = jldopen(simulation.output_writers[:velocities].filepath)

## Extracts a vector of `iterations` at which data was saved.
iterations = parse.(Int, keys(file["timeseries/t"]))

# and makes an animation with Plots.jl:

anim = @animate for (i, iter) in enumerate(iterations)

    @info "Drawing frame $i from iteration $iter..."

    w = file["timeseries/w/$iter"][:, 1, :]
    t = file["timeseries/t/$iter"]

    contourf(x, z, w', title = @sprintf("ωt = %.2f", ω * t),
                      levels = range(-1e-8, stop=1e-8, length=10),
                       clims = (-1e-8, 1e-8),
                      xlabel = "x",
                      ylabel = "z",
                       xlims = (-π, π),
                       ylims = (-π, π),
                   linewidth = 0,
                       color = :balance,
                      legend = false,
                 aspectratio = :equal)
end

gif(anim, "internal_wave.gif", fps = 8) # hide


text:
using Statistics, StatsBase

export relayed, isi_efficacy, efficacy

const TS = Vector{Float64}

# ============================================================================ #
function relayed(pre::TS, post::TS, kbin::AbstractArray{Int,1},
    bin_size::Real, nbase::Integer, thr::Real)

    # NOTE: we are looking at the spiking of pre relative to post, so that
    # the indicies in <grp> refer to the the pre spike train
    xcm, grp = psth(pre, post, kbin, bin_size, Binwise)
    xc = vec(sum(xcm, dims=2))

    _, kmax = findmax(xc)
    baseline = xc[[1:nbase; end-(nbase-1):end]]
    threshold = mean(baseline) + (std(baseline) * thr)

    kstart = kmax - (findfirst(<(threshold), xc[kmax:-1:1]) - 1)
    kstart = 1 <= kstart < kmax ? kstart : kmax

    kend = kmax + (findfirst(<(threshold), xc[kmax:end]) - 1)
    kend = length(xc) >= kend > kmax ? kend : kmax

    return sort(unique(cat(grp[kstart:kend]..., dims=1)))
end
# ============================================================================ #
"""
`eff, lab = isi_efficacy(pre::TS, post::TS, bin_size::Float64,
    isi_bin_size::Float64, isimax::Float64=0.03, baseline::Float64=0.01,
    deadtime::Float64=0.02, max_delay::Float64=0.015)`
### Inputs:
* `pre` - timestamps of presynaptic spikes (in seconds)
* `post` - timestamps of postsynaptic spikes (in seconds)
* `bin_size` - psth/xcorr bin size (in seconds)
* `isi_bin_size` - bin size for isi discretization
* `isimax` - (0.03) max ISI length to include in the analysis (in seconds)
* `baseline` - (0.01) baseline window size (in seconds) used to construct threshold for identifying peak bins
* `deadtime` - (0.02) minimum deadtime preceeding first spike
* `max_delay` - (+/-0.015) the max synaptic delay to use to constrain search for xcorr peak (in seconds)

### Output:
* `eff` - a vector of % efficacy values
* `lab` - a vector of bin labels for <eff>

#### Notes:
* output `eff` will have NaN entries where the ISI distribution == 0
"""
function isi_efficacy(pre::TS, post::TS, bin_size::Float64,
    isi_bin_size::Float64, isimax::Float64=0.03, baseline::Float64=0.01,
    deadtime::Float64=0.02, max_delay::Float64=0.015)

    bin = round(Int64, (max_delay + baseline) / bin_size)
    kbin = -bin:+bin

    nbase = round(Int64, baseline / bin_size)

    # find all relayed spikes: by performing this operation on the full spike
    # train we minimize the liklihood of missing the true peak
    krel = relayed(pre, post, kbin, bin_size, nbase, 3.0)

    # preceeding isi (i.e. duration between each spike and the spike before it)
    # NOTE: the last spike cannot be the first spike of a pair (hence [1:end-1])
    pre_isi = [0.0; diff(pre[1:end-1])]

    # all pre spikes with a pre_isi >= deadtime (20ms in Usrey et al. 1998)
    first_spk = findall(x -> x >= deadtime, pre_isi)

    second_spk = first_spk .+ 1

    # relayed second spikes are simply the intersection of all second spikes and
    # all relayed spikes (indicies)
    krel_second = findall(in(krel), second_spk)

    # post isi (i.e. isi between first (<pre_isi> spikes) and second spikes)
    post_isi = pre[second_spk] - pre[first_spk]

    edges = 0.0:isi_bin_size:isimax+isi_bin_size

    # isi distribution of relayed second spikes divided by isi distribution of
    # all second spikes
    rel = fit(Histogram, post_isi[krel_second], edges)
    all = fit(Histogram, post_isi, edges)

    return rel.weights ./ all.weights, edges[1:end-1]
end
# ============================================================================ #
"""
`efficacy(pre, post, evt, dur; bin_size=0.001, baseline=0.01, max_delay=0.015)`

Trialwise efficacy calculation

### Inputs:
* pre - vector of timestamps from presynaptic cell
* post - vector of timestamps from postsynaptic cell
* evt - vector of event timestamps
* dur - trial duration in seconds
* bin_size - (0.001) bin size in seconds
* baseline - (0.01) baseline window size (in seconds)
* max_delay - (0.015) the max synaptic delay in seconds

### Outputs:
* ef - a vector of efficacy values, one for each trial
"""
function efficacy(pre::TS, post::TS, evt::TS, dur::Real,
    bin_size::Real=.001, baseline::Real=0.01, max_delay::Real=0.015)

    bin = round(Int64, (max_delay + baseline) / bin_size)
    kbin = -bin:+bin

    nbase = round(Int64, baseline / bin_size)

    # get indicies of all relayed spikes
    krel = relayed(pre, post, kbin, bin_size, nbase, 3.0)

    kbin = 0:round(Int64, dur / bin_size)

    pre_grp = psth(pre, evt, kbin, bin_size, Trialwise)[2]

    ef = zeros(length(pre_grp))
    @inbounds for k = 1:length(pre_grp)
        # ef[k] = length(findin(pre_grp[k], krel)) / length(pre_grp[k])
        ef[k] = length(findall(in(krel), pre_grp[k])) / length(pre_grp[k])
    end

    return ef
end
# ============================================================================ #


text:
#
# example of how you can get a const typed array into global scope,
# by wrapping it in a struct and relying on the fact that the array
# inself will always be mutable
#
module HHLDs

    using BenchmarkTools
    using InteractiveUtils

    struct HH{T<:Real}
      rent :: T
      rates :: T
    end

    function tenkhhlds(t::Type)::Vector
        a = Vector{HH{t}}(undef,10_000)
        for i in 1:10_000
            a[i] = HH{t}(i, i*2)
        end
        a
    end

    struct A
       m :: Vector{HH{Float64}}
     end

    const AA = A(tenkhhlds(Float64))
    BB = A(tenkhhlds(Float64))
    CC = tenkhhlds(Float64)

    function f(x::A)::Real
        s = 0.0
        for i in x.m
          s += i.rates
        end
        s
    end

    function f(x::Vector{HH{Float64}})::Real
        s = 0.0
        for i in x
          s += i.rates
        end
        s
    end

    function getAA(i::Int)::HH{Float64}
        AA.m[i]
    end

    function getBB(i::Int)::HH{Float64}
        BB.m[i]
    end

    function getCC(i::Int)::HH{Float64}
        CC[i]
    end

    function fgAA()::Real
        s = 0.0
        for i in 1:10_000
            hh = getAA(i)
            s += hh.rates
        end
        s
    end

    function fgBB()::Real
        s = 0.0
        for i in 1:10_000
            hh = getBB(i)
            s += hh.rates
        end
        s
    end

    function fgCC()::Real
        s = 0.0
        for i in 1:10_000
            hh = getCC(i)
            s += hh.rates
        end
        s
    end

    function do_all()
        hhbench = BenchmarkGroup()
        hhbench[:direct_constant] = @benchmark f(AA)
        hhbench[:direct_variable_struct] = @benchmark f(BB)
        hhbench[:direct_variable] = @benchmark f(CC)
        hhbench[:getter_struct_const] = @benchmark fgAA()
        hhbench[:getter_struct_var] = @benchmark fgBB()
        hhbench[:getter_array_variable] = @benchmark fgCC()

        hhbench
    end

end

# push!(HHLDs.AA.m,HHLDs.HH(10.0,10.0))
# or ..
# push!(HHLDs.AA.m,HHLDs.HH{Float64}(11,11))

HHLDs.do_all()

@code_warntype HHLDs.getCC(1)
# @code_warntype HHLDs.getAA(1)

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
# Speed up rotation search of TIV-GOGMA by using the ROCS alignment algorithm to produce a best rotation
function rocs_gogma_align(gmmx::AbstractGMM, gmmy::AbstractGMM; kwargs...)
    t = promote_type(numbertype(gmmx),numbertype(gmmy))
    p = t(π)

    # rotation alignment
    rocs_res = rocs_align(gmmx, gmmy)
    rotpos = rot_to_axis(rocs_res.tform.linear)

    # translation alignment
    trl_res = trl_gogma_align(gmmx, gmmy; rot=rotpos, kwargs...)

    # local alignment in the full transformation space
    pos = NTuple{6, t}([rotpos..., trl_res.tform_params...])
    trlim = translation_limit(gmmx, gmmy)
    localblock = Block(((-p,p), (-p,p), (-p,p), (-trlim,trlim), (-trlim,trlim), (-trlim,trlim)), pos, trl_res.upperbound, -Inf)
    min, bestpos = local_align(gmmx, gmmy, localblock)

    return GMAlignmentResult(gmmx, gmmy, min, trl_res.lowerbound, AffineMap(bestpos...), bestpos, trl_res.obj_calls, trl_res.num_splits, trl_res.num_blocks, trl_res.stagnant_splits)

end


text:
cd(@__DIR__); runapprox = true; allNGWPs = true; include("setups/sunflower.jl")
gr(dpi = 200)

## (a) barbara sunflower pants
heatmap(barbara, yflip=true, ratio=1, c=:greys)
scatter_gplot!(transform2D(X; s = 20, t = [280, 320]); ms = 2, c = :red)
sample_location_plt = plot!(cbar = false, frame = :none)
savefig(sample_location_plt, "../figs/barb_sunflower_pants.png")

## (b) barbara eye graph signal
f = matread("../datasets/sunflower_barbara_voronoi.mat")["f_trouser_voronoi"]
G_Sig.f = reshape(f, (N, 1))
scatter_gplot(X; marker = f, ms = LinRange(4.0, 14.0, N), c = :greys);
plt = plot!(xlim = [-1.2,1.2], ylim = [-1.2,1.2], frame = :none)
savefig(plt, "../figs/SunFlower_barbara_ftrouser.png")

## (c) barbara eye relative l2 approximation error by various methods
DVEC = getall_expansioncoeffs2(G_Sig, GP_dual, GP_dual_Lsym, VM_NGWP, PC_NGWP, LP_NGWP,
                               VM_NGWP_Lsym, PC_NGWP_Lsym, LP_NGWP_Lsym, 𝚽, 𝚽sym)
approx_error_plot2(DVEC);
plt = plot!(xguidefontsize = 14, yguidefontsize = 14, legendfontsize = 11, size = (600, 600))
savefig(plt, "../figs/SunFlower_barbara_ftrouser_DAG_approx.png")


text:
function SplitStr(s12, sep)
#Given an array of joint array with a separater for each string in the array,
#Split the array in two array with onside of the array on each array.
#Assuem that separator is always inside the elements of s12

    strsplit = [split(x,sep)[y] for x in s12, y in 1:2]
    return strsplit[:,1],strsplit[:,2]

end


text:
"""
	struct ExtractScalar{T}
		datatype::Type{T}
		c::T
		s::T
	end

	extract a scalar value and center. If `T` is of type number, it is centered by first subtracting c and then
	multiplying that with s.
"""
struct ExtractScalar{T<:Number,V} <: AbstractExtractor
	datatype::Type{T}
	c::V
	s::V
end

ExtractScalar(datatype) = ExtractScalar(datatype, zero(datatype), one(datatype))
extractsmatrix(s::ExtractScalar) = true

(s::ExtractScalar{T,V})(v) where {T<:Number,V}			 = ArrayNode(s.s .* (fill(s.datatype(v),1,1) .- s.c))
(s::ExtractScalar{T,V} where {V,T<:Number})(v::String)   = s((parse(s.datatype,v)))
(s::ExtractScalar{T,V})(v::S) where {T<:Number,V,S<:Nothing}= ArrayNode(fill(zero(T),(1,1)))
function Base.show(io::IO, m::ExtractScalar;pad = [], key::String="")
	c = COLORS[(length(pad)%length(COLORS))+1]
	key *= isempty(key) ? "" : ": ";
	paddedprint(io,"$(key)$(m.datatype)\n", color = c)
end

Base.hash(e::ExtractScalar, h::UInt) = hash((e.datatype, e.c, e.s), h)
Base.:(==)(e1::ExtractScalar, e2::ExtractScalar) = e1.datatype == e2.datatype && e1.c === e2.c && e1.s === e2.s


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Base: @propagate_inbounds
import Base: length, size, axes, IndexStyle, getindex, setindex!, parent, vec, convert, similar

### basic definitions (types, aliases, constructors, abstractarray interface, sundry similar)

# note that Adjoint and Transpose must be able to wrap not only vectors and matrices
# but also factorizations, rotations, and other linear algebra objects, including
# user-defined such objects. so do not restrict the wrapped type.
"""
    Adjoint

Lazy wrapper type for an adjoint view of the underlying linear algebra object,
usually an `AbstractVector`/`AbstractMatrix`, but also some `Factorization`, for instance.
Usually, the `Adjoint` constructor should not be called directly, use [`adjoint`](@ref)
instead. To materialize the view use [`copy`](@ref).

This type is intended for linear algebra usage - for general data manipulation see
[`permutedims`](@ref Base.permutedims).

# Examples
```jldoctest
julia> A = [3+2im 9+2im; 0 0]
2×2 Matrix{Complex{Int64}}:
 3+2im  9+2im
 0+0im  0+0im

julia> Adjoint(A)
2×2 adjoint(::Matrix{Complex{Int64}}) with eltype Complex{Int64}:
 3-2im  0+0im
 9-2im  0+0im
```
"""
struct Adjoint{T,S} <: AbstractMatrix{T}
    parent::S
end
"""
    Transpose

Lazy wrapper type for a transpose view of the underlying linear algebra object,
usually an `AbstractVector`/`AbstractMatrix`, but also some `Factorization`, for instance.
Usually, the `Transpose` constructor should not be called directly, use [`transpose`](@ref)
instead. To materialize the view use [`copy`](@ref).

This type is intended for linear algebra usage - for general data manipulation see
[`permutedims`](@ref Base.permutedims).

# Examples
```jldoctest
julia> A = [2 3; 0 0]
2×2 Matrix{Int64}:
 2  3
 0  0

julia> Transpose(A)
2×2 transpose(::Matrix{Int64}) with eltype Int64:
 2  0
 3  0
```
"""
struct Transpose{T,S} <: AbstractMatrix{T}
    parent::S
end

# basic outer constructors
Adjoint(A) = Adjoint{Base.promote_op(adjoint,eltype(A)),typeof(A)}(A)
Transpose(A) = Transpose{Base.promote_op(transpose,eltype(A)),typeof(A)}(A)

Base.dataids(A::Union{Adjoint, Transpose}) = Base.dataids(A.parent)
Base.unaliascopy(A::Union{Adjoint,Transpose}) = typeof(A)(Base.unaliascopy(A.parent))

# wrapping lowercase quasi-constructors
"""
    A'
    adjoint(A)

Lazy adjoint (conjugate transposition). Note that `adjoint` is applied recursively to
elements.

For number types, `adjoint` returns the complex conjugate, and therefore it is equivalent to
the identity function for real numbers.

This operation is intended for linear algebra usage - for general data manipulation see
[`permutedims`](@ref Base.permutedims).

# Examples
```jldoctest
julia> A = [3+2im 9+2im; 0  0]
2×2 Matrix{Complex{Int64}}:
 3+2im  9+2im
 0+0im  0+0im

julia> B = A' # equivalently adjoint(A)
2×2 adjoint(::Matrix{Complex{Int64}}) with eltype Complex{Int64}:
 3-2im  0+0im
 9-2im  0+0im

julia> B isa Adjoint
true

julia> adjoint(B) === A # the adjoint of an adjoint unwraps the parent
true

julia> Adjoint(B) # however, the constructor always wraps its argument
2×2 adjoint(adjoint(::Matrix{Complex{Int64}})) with eltype Complex{Int64}:
 3+2im  9+2im
 0+0im  0+0im

julia> B[1,2] = 4 + 5im; # modifying B will modify A automatically

julia> A
2×2 Matrix{Complex{Int64}}:
 3+2im  9+2im
 4-5im  0+0im
```

For real matrices, the `adjoint` operation is equivalent to a `transpose`.

```jldoctest
julia> A = reshape([x for x in 1:4], 2, 2)
2×2 Matrix{Int64}:
 1  3
 2  4

julia> A'
2×2 adjoint(::Matrix{Int64}) with eltype Int64:
 1  2
 3  4

julia> adjoint(A) == transpose(A)
true
```

The adjoint of an `AbstractVector` is a row-vector:
```jldoctest
julia> x = [3, 4im]
2-element Vector{Complex{Int64}}:
 3 + 0im
 0 + 4im

julia> x'
1×2 adjoint(::Vector{Complex{Int64}}) with eltype Complex{Int64}:
 3+0im  0-4im

julia> x'x # compute the dot product, equivalently x' * x
25 + 0im
```

For a matrix of matrices, the individual blocks are recursively operated on:
```jldoctest
julia> A = reshape([x + im*x for x in 1:4], 2, 2)
2×2 Matrix{Complex{Int64}}:
 1+1im  3+3im
 2+2im  4+4im

julia> C = reshape([A, 2A, 3A, 4A], 2, 2)
2×2 Matrix{Matrix{Complex{Int64}}}:
 [1+1im 3+3im; 2+2im 4+4im]  [3+3im 9+9im; 6+6im 12+12im]
 [2+2im 6+6im; 4+4im 8+8im]  [4+4im 12+12im; 8+8im 16+16im]

julia> C'
2×2 adjoint(::Matrix{Matrix{Complex{Int64}}}) with eltype Adjoint{Complex{Int64}, Matrix{Complex{Int64}}}:
 [1-1im 2-2im; 3-3im 4-4im]    [2-2im 4-4im; 6-6im 8-8im]
 [3-3im 6-6im; 9-9im 12-12im]  [4-4im 8-8im; 12-12im 16-16im]
```
"""
adjoint(A::AbstractVecOrMat) = Adjoint(A)

"""
    transpose(A)

Lazy transpose. Mutating the returned object should appropriately mutate `A`. Often,
but not always, yields `Transpose(A)`, where `Transpose` is a lazy transpose wrapper. Note
that this operation is recursive.

This operation is intended for linear algebra usage - for general data manipulation see
[`permutedims`](@ref Base.permutedims), which is non-recursive.

# Examples
```jldoctest
julia> A = [3 2; 0 0]
2×2 Matrix{Int64}:
 3  2
 0  0

julia> B = transpose(A)
2×2 transpose(::Matrix{Int64}) with eltype Int64:
 3  0
 2  0

julia> B isa Transpose
true

julia> transpose(B) === A # the transpose of a transpose unwraps the parent
true

julia> Transpose(B) # however, the constructor always wraps its argument
2×2 transpose(transpose(::Matrix{Int64})) with eltype Int64:
 3  2
 0  0

julia> B[1,2] = 4; # modifying B will modify A automatically

julia> A
2×2 Matrix{Int64}:
 3  2
 4  0
```

For complex matrices, the `adjoint` operation is equivalent to a conjugate-transpose.
```jldoctest
julia> A = reshape([Complex(x, x) for x in 1:4], 2, 2)
2×2 Matrix{Complex{Int64}}:
 1+1im  3+3im
 2+2im  4+4im

julia> adjoint(A) == conj(transpose(A))
true
```

The `transpose` of an `AbstractVector` is a row-vector:
```jldoctest
julia> v = [1,2,3]
3-element Vector{Int64}:
 1
 2
 3

julia> transpose(v) # returns a row-vector
1×3 transpose(::Vector{Int64}) with eltype Int64:
 1  2  3

julia> transpose(v) * v # compute the dot product
14
```

For a matrix of matrices, the individual blocks are recursively operated on:
```jldoctest
julia> C = reshape(1:4, 2, 2)
2×2 reshape(::UnitRange{Int64}, 2, 2) with eltype Int64:
 1  3
 2  4

julia> D = reshape([C, 2C, 3C, 4C], 2, 2) # construct a block matrix
2×2 Matrix{Matrix{Int64}}:
 [1 3; 2 4]  [3 9; 6 12]
 [2 6; 4 8]  [4 12; 8 16]

julia> transpose(D) # blocks are recursively transposed
2×2 transpose(::Matrix{Matrix{Int64}}) with eltype Transpose{Int64, Matrix{Int64}}:
 [1 2; 3 4]   [2 4; 6 8]
 [3 6; 9 12]  [4 8; 12 16]
```
"""
transpose(A::AbstractVecOrMat) = Transpose(A)

# unwrapping lowercase quasi-constructors
adjoint(A::Adjoint) = A.parent
transpose(A::Transpose) = A.parent
adjoint(A::Transpose{<:Real}) = A.parent
transpose(A::Adjoint{<:Real}) = A.parent

# printing
function Base.showarg(io::IO, v::Adjoint, toplevel)
    print(io, "adjoint(")
    Base.showarg(io, parent(v), false)
    print(io, ')')
    toplevel && print(io, " with eltype ", eltype(v))
    return nothing
end
function Base.showarg(io::IO, v::Transpose, toplevel)
    print(io, "transpose(")
    Base.showarg(io, parent(v), false)
    print(io, ')')
    toplevel && print(io, " with eltype ", eltype(v))
    return nothing
end

# some aliases for internal convenience use
const AdjOrTrans{T,S} = Union{Adjoint{T,S},Transpose{T,S}} where {T,S}
const AdjointAbsVec{T} = Adjoint{T,<:AbstractVector}
const AdjointAbsMat{T} = Adjoint{T,<:AbstractMatrix}
const TransposeAbsVec{T} = Transpose{T,<:AbstractVector}
const TransposeAbsMat{T} = Transpose{T,<:AbstractMatrix}
const AdjOrTransAbsVec{T} = AdjOrTrans{T,<:AbstractVector}
const AdjOrTransAbsMat{T} = AdjOrTrans{T,<:AbstractMatrix}

# for internal use below
wrapperop(_) = identity
wrapperop(::Adjoint) = adjoint
wrapperop(::Transpose) = transpose

# AbstractArray interface, basic definitions
length(A::AdjOrTrans) = length(A.parent)
size(v::AdjOrTransAbsVec) = (1, length(v.parent))
size(A::AdjOrTransAbsMat) = reverse(size(A.parent))
axes(v::AdjOrTransAbsVec) = (Base.OneTo(1), axes(v.parent)...)
axes(A::AdjOrTransAbsMat) = reverse(axes(A.parent))
IndexStyle(::Type{<:AdjOrTransAbsVec}) = IndexLinear()
IndexStyle(::Type{<:AdjOrTransAbsMat}) = IndexCartesian()
@propagate_inbounds getindex(v::AdjOrTransAbsVec{T}, i::Int) where {T} = wrapperop(v)(v.parent[i-1+first(axes(v.parent)[1])])::T
@propagate_inbounds getindex(A::AdjOrTransAbsMat{T}, i::Int, j::Int) where {T} = wrapperop(A)(A.parent[j, i])::T
@propagate_inbounds setindex!(v::AdjOrTransAbsVec, x, i::Int) = (setindex!(v.parent, wrapperop(v)(x), i-1+first(axes(v.parent)[1])); v)
@propagate_inbounds setindex!(A::AdjOrTransAbsMat, x, i::Int, j::Int) = (setindex!(A.parent, wrapperop(A)(x), j, i); A)
# AbstractArray interface, additional definitions to retain wrapper over vectors where appropriate
@propagate_inbounds getindex(v::AdjOrTransAbsVec, ::Colon, is::AbstractArray{Int}) = wrapperop(v)(v.parent[is])
@propagate_inbounds getindex(v::AdjOrTransAbsVec, ::Colon, ::Colon) = wrapperop(v)(v.parent[:])

# conversion of underlying storage
convert(::Type{Adjoint{T,S}}, A::Adjoint) where {T,S} = Adjoint{T,S}(convert(S, A.parent))
convert(::Type{Transpose{T,S}}, A::Transpose) where {T,S} = Transpose{T,S}(convert(S, A.parent))

# Strides and pointer for transposed strided arrays — but only if the elements are actually stored in memory
Base.strides(A::Adjoint{<:Real, <:AbstractVector}) = (stride(A.parent, 2), stride(A.parent, 1))
Base.strides(A::Transpose{<:Any, <:AbstractVector}) = (stride(A.parent, 2), stride(A.parent, 1))
# For matrices it's slightly faster to use reverse and avoid calling stride twice
Base.strides(A::Adjoint{<:Real, <:AbstractMatrix}) = reverse(strides(A.parent))
Base.strides(A::Transpose{<:Any, <:AbstractMatrix}) = reverse(strides(A.parent))

Base.unsafe_convert(::Type{Ptr{T}}, A::Adjoint{<:Real, <:AbstractVecOrMat}) where {T} = Base.unsafe_convert(Ptr{T}, A.parent)
Base.unsafe_convert(::Type{Ptr{T}}, A::Transpose{<:Any, <:AbstractVecOrMat}) where {T} = Base.unsafe_convert(Ptr{T}, A.parent)

Base.elsize(::Type{<:Adjoint{<:Real, P}}) where {P<:AbstractVecOrMat} = Base.elsize(P)
Base.elsize(::Type{<:Transpose{<:Any, P}}) where {P<:AbstractVecOrMat} = Base.elsize(P)

# for vectors, the semantics of the wrapped and unwrapped types differ
# so attempt to maintain both the parent and wrapper type insofar as possible
similar(A::AdjOrTransAbsVec) = wrapperop(A)(similar(A.parent))
similar(A::AdjOrTransAbsVec, ::Type{T}) where {T} = wrapperop(A)(similar(A.parent, Base.promote_op(wrapperop(A), T)))
# for matrices, the semantics of the wrapped and unwrapped types are generally the same
# and as you are allocating with similar anyway, you might as well get something unwrapped
similar(A::AdjOrTrans) = similar(A.parent, eltype(A), axes(A))
similar(A::AdjOrTrans, ::Type{T}) where {T} = similar(A.parent, T, axes(A))
similar(A::AdjOrTrans, ::Type{T}, dims::Dims{N}) where {T,N} = similar(A.parent, T, dims)

# AbstractMatrix{T} constructor for adjtrans vector: preserve wrapped type
AbstractMatrix{T}(A::AdjOrTransAbsVec) where {T} = wrapperop(A)(AbstractVector{T}(A.parent))

# sundry basic definitions
parent(A::AdjOrTrans) = A.parent
vec(v::TransposeAbsVec{<:Number}) = parent(v)
vec(v::AdjointAbsVec{<:Real}) = parent(v)

### concatenation
# preserve Adjoint/Transpose wrapper around vectors
# to retain the associated semantics post-concatenation
hcat(avs::Union{Number,AdjointAbsVec}...) = _adjoint_hcat(avs...)
hcat(tvs::Union{Number,TransposeAbsVec}...) = _transpose_hcat(tvs...)
_adjoint_hcat(avs::Union{Number,AdjointAbsVec}...) = adjoint(vcat(map(adjoint, avs)...))
_transpose_hcat(tvs::Union{Number,TransposeAbsVec}...) = transpose(vcat(map(transpose, tvs)...))
typed_hcat(::Type{T}, avs::Union{Number,AdjointAbsVec}...) where {T} = adjoint(typed_vcat(T, map(adjoint, avs)...))
typed_hcat(::Type{T}, tvs::Union{Number,TransposeAbsVec}...) where {T} = transpose(typed_vcat(T, map(transpose, tvs)...))
# otherwise-redundant definitions necessary to prevent hitting the concat methods in LinearAlgebra/special.jl
hcat(avs::Adjoint{<:Any,<:Vector}...) = _adjoint_hcat(avs...)
hcat(tvs::Transpose{<:Any,<:Vector}...) = _transpose_hcat(tvs...)
hcat(avs::Adjoint{T,Vector{T}}...) where {T} = _adjoint_hcat(avs...)
hcat(tvs::Transpose{T,Vector{T}}...) where {T} = _transpose_hcat(tvs...)
# TODO unify and allow mixed combinations

### higher order functions
# preserve Adjoint/Transpose wrapper around vectors
# to retain the associated semantics post-map/broadcast
#
# note that the caller's operation f operates in the domain of the wrapped vectors' entries.
# hence the adjoint->f->adjoint shenanigans applied to the parent vectors' entries.
map(f, avs::AdjointAbsVec...) = adjoint(map((xs...) -> adjoint(f(adjoint.(xs)...)), parent.(avs)...))
map(f, tvs::TransposeAbsVec...) = transpose(map((xs...) -> transpose(f(transpose.(xs)...)), parent.(tvs)...))
quasiparentt(x) = parent(x); quasiparentt(x::Number) = x # to handle numbers in the defs below
quasiparenta(x) = parent(x); quasiparenta(x::Number) = conj(x) # to handle numbers in the defs below
broadcast(f, avs::Union{Number,AdjointAbsVec}...) = adjoint(broadcast((xs...) -> adjoint(f(adjoint.(xs)...)), quasiparenta.(avs)...))
broadcast(f, tvs::Union{Number,TransposeAbsVec}...) = transpose(broadcast((xs...) -> transpose(f(transpose.(xs)...)), quasiparentt.(tvs)...))
# Hack to preserve behavior after #32122; this needs to be done with a broadcast style instead to support dotted fusion
Broadcast.broadcast_preserving_zero_d(f, avs::Union{Number,AdjointAbsVec}...) = adjoint(broadcast((xs...) -> adjoint(f(adjoint.(xs)...)), quasiparenta.(avs)...))
Broadcast.broadcast_preserving_zero_d(f, tvs::Union{Number,TransposeAbsVec}...) = transpose(broadcast((xs...) -> transpose(f(transpose.(xs)...)), quasiparentt.(tvs)...))
# TODO unify and allow mixed combinations with a broadcast style

### reductions
# faster to sum the Array than to work through the wrapper
Base._mapreduce_dim(f, op, init::Base._InitialValue, A::Transpose, dims::Colon) =
    transpose(Base._mapreduce_dim(_sandwich(transpose, f), _sandwich(transpose, op), init, parent(A), dims))
Base._mapreduce_dim(f, op, init::Base._InitialValue, A::Adjoint, dims::Colon) =
    adjoint(Base._mapreduce_dim(_sandwich(adjoint, f), _sandwich(adjoint, op), init, parent(A), dims))
# sum(A'; dims)
Base.mapreducedim!(f, op, B::AbstractArray, A::TransposeAbsMat) =
    transpose(Base.mapreducedim!(_sandwich(transpose, f), _sandwich(transpose, op), transpose(B), parent(A)))
Base.mapreducedim!(f, op, B::AbstractArray, A::AdjointAbsMat) =
    adjoint(Base.mapreducedim!(_sandwich(adjoint, f), _sandwich(adjoint, op), adjoint(B), parent(A)))

_sandwich(adj::Function, fun) = (xs...,) -> adj(fun(map(adj, xs)...))
for fun in [:identity, :add_sum, :mul_prod] #, :max, :min]
    @eval _sandwich(::Function, ::typeof(Base.$fun)) = Base.$fun
end

### linear algebra

(-)(A::Adjoint)   = Adjoint(  -A.parent)
(-)(A::Transpose) = Transpose(-A.parent)

## multiplication *

function _dot_nonrecursive(u, v)
    lu = length(u)
    if lu != length(v)
        throw(DimensionMismatch("first array has length $(lu) which does not match the length of the second, $(length(v))."))
    end
    if lu == 0
        zero(eltype(u)) * zero(eltype(v))
    else
        sum(uu*vv for (uu, vv) in zip(u, v))
    end
end

# Adjoint/Transpose-vector * vector
*(u::AdjointAbsVec{<:Number}, v::AbstractVector{<:Number}) = dot(u.parent, v)
*(u::TransposeAbsVec{T}, v::AbstractVector{T}) where {T<:Real} = dot(u.parent, v)
*(u::AdjOrTransAbsVec, v::AbstractVector) = _dot_nonrecursive(u, v)

# vector * Adjoint/Transpose-vector
*(u::AbstractVector, v::AdjOrTransAbsVec) = broadcast(*, u, v)
# Adjoint/Transpose-vector * Adjoint/Transpose-vector
# (necessary for disambiguation with fallback methods in linalg/matmul)
*(u::AdjointAbsVec, v::AdjointAbsVec) = throw(MethodError(*, (u, v)))
*(u::TransposeAbsVec, v::TransposeAbsVec) = throw(MethodError(*, (u, v)))

# AdjOrTransAbsVec{<:Any,<:AdjOrTransAbsVec} is a lazy conj vectors
# We need to expand the combinations to avoid ambiguities
(*)(u::TransposeAbsVec, v::AdjointAbsVec{<:Any,<:TransposeAbsVec}) = _dot_nonrecursive(u, v)
(*)(u::AdjointAbsVec,   v::AdjointAbsVec{<:Any,<:TransposeAbsVec}) = _dot_nonrecursive(u, v)
(*)(u::TransposeAbsVec, v::TransposeAbsVec{<:Any,<:AdjointAbsVec}) = _dot_nonrecursive(u, v)
(*)(u::AdjointAbsVec,   v::TransposeAbsVec{<:Any,<:AdjointAbsVec}) = _dot_nonrecursive(u, v)

## pseudoinversion
pinv(v::AdjointAbsVec, tol::Real = 0) = pinv(v.parent, tol).parent
pinv(v::TransposeAbsVec, tol::Real = 0) = pinv(conj(v.parent)).parent

## left-division \
\(u::AdjOrTransAbsVec, v::AdjOrTransAbsVec) = pinv(u) * v
\(u::AdjointAbsVec, y::Number) = adjoint(conj(y) / u.parent)
\(u::TransposeAbsVec, y::Number) = transpose(y / u.parent)

## right-division /
/(u::AdjointAbsVec, A::AbstractMatrix) = adjoint(adjoint(A) \ u.parent)
/(u::TransposeAbsVec, A::AbstractMatrix) = transpose(transpose(A) \ u.parent)
/(u::AdjointAbsVec, A::Transpose{<:Any,<:AbstractMatrix}) = adjoint(conj(A.parent) \ u.parent) # technically should be adjoint(copy(adjoint(copy(A))) \ u.parent)
/(u::TransposeAbsVec, A::Adjoint{<:Any,<:AbstractMatrix}) = transpose(conj(A.parent) \ u.parent) # technically should be transpose(copy(transpose(copy(A))) \ u.parent)

## complex conjugate
conj(A::Transpose) = adjoint(A.parent)
conj(A::Adjoint) = transpose(A.parent)

max_stars_repo_name: TechPenguineer/julia

max_issues_repo_name: TechPenguineer/julia

max_forks_repo_name: TechPenguineer/julia


text:
abstract type AbstractWorkCost end

@inline (mdl::AbstractWorkCost)(n_vertices, n_pins, k) = mdl(n_vertices, n_pins)

struct AffineWorkModel{Tv} <: AbstractWorkCost
    α::Tv
    β_vertex::Tv
    β_pin::Tv
end

function AffineWorkModel(; α = false, β_vertex = false, β_pin = false)
    AffineWorkModel(promote(α, β_vertex, β_pin)...)
end

@inline cost_type(::Type{AffineWorkModel{Tv}}) where {Tv} = Tv

(mdl::AffineWorkModel)(n_vertices, n_pins) = mdl.α + n_vertices * mdl.β_vertex + n_pins * mdl.β_pin

struct WorkOracle{Ti, Mdl <: AbstractWorkCost} <: AbstractOracleCost{Mdl}
    pos::Vector{Ti}
    mdl::Mdl
end

oracle_model(ocl::WorkOracle) = ocl.mdl

function oracle_stripe(hint::AbstractHint, mdl::AbstractWorkCost, A::SparseMatrixCSC; kwargs...)
    return WorkOracle(A.colptr, mdl)
end

@inline function (cst::WorkOracle{Ti, Mdl})(j::Ti, j′::Ti, k...) where {Ti, Mdl}
    @inbounds begin
        w = cst.pos[j′] - cst.pos[j]
        return cst.mdl(j′ - j, w, k...)
    end
end

bound_stripe(A::SparseMatrixCSC, K, ocl::WorkOracle{<:Any, <:AffineWorkModel}) =
    bound_stripe(A, K, oracle_model(ocl))
function bound_stripe(A::SparseMatrixCSC, K, mdl::AffineWorkModel)
    m, n = size(A)
    N = nnz(A)
    c_lo = mdl.α + fld(mdl.β_vertex * n + mdl.β_pin * N, K)
    if mdl.β_vertex ≥ 0 && mdl.β_pin ≥ 0
        c_hi = mdl.α + mdl.β_vertex * n + mdl.β_pin * N
    elseif mdl.β_vertex ≤ 0 && mdl.β_pin ≤ 0
        c_hi = mdl.α
    else
        @assert false
    end
    return (c_lo, c_hi)
end

function compute_objective(g::G, A::SparseMatrixCSC, Π::SplitPartition, mdl::AbstractWorkCost) where {G}
    cst = objective_identity(g, cost_type(mdl))
    for k = 1:Π.K
        j = Π.spl[k]
        j′ = Π.spl[k + 1]
        cst = g(cst, mdl(j′ - j, A.colptr[j′] - A.colptr[j], k))
    end
    return cst
end

function compute_objective(g::G, A::SparseMatrixCSC, Π::DomainPartition, mdl::AbstractWorkCost) where {G}
    cst = objective_identity(g, cost_type(mdl))
    for k = 1:Π.K
        s = Π.spl[k]
        s′ = Π.spl[k + 1]
        n_vertices = s′ - s
        n_pins = 0
        for _s = s : s′ - 1
            j = Π.prm[_s]
            n_pins += A.colptr[j + 1] - A.colptr[j]
        end
        cst = g(cst, mdl(n_vertices, n_pins, k))
    end
    return cst
end

function compute_objective(g, A::SparseMatrixCSC, Π::MapPartition, mdl::AbstractWorkCost)
    return compute_objective(g, A, convert(DomainPartition, Π), mdl)
end


text:
# source: https://github.com/dirk/HTTP.jl/blob/master/src/HTTP/Util.jl

escaped_regex = r"%([0-9a-fA-F]{2})"

function unescape(str)
    # def _unescape(str, regex) str.gsub(regex){ $1.hex.chr } end
    for m in eachmatch(escaped_regex, str)
        for capture in m.captures
            rep = string(char(parseint(capture, 16)))
            str = replace(str, "%"*capture, rep)
        end
    end
    return str
end

function unescape_form(str)
    str = replace(str, "+", " ")
    return unescape(str)
end

# Escaping
control_array = convert(Array{Uint8,1}, [vec(0:(parseint("1f", 16)))])
control = utf8(ascii(control_array)*"\x7f")
space = utf8(" ")
delims = utf8("%<>\"")
unwise   = utf8("{}|\\^`")
nonascii_array = convert(Array{Uint8,1}, [vec(parseint("80", 16):(parseint("ff", 16)))])

reserved = utf8(",;/?:@&=+\$![]'*#")
# Strings to be escaped
# (Delims goes first so '%' gets escaped first.)
unescaped = delims * reserved * control * space * unwise# * nonascii
unescaped_form = delims * reserved * control * unwise# * nonascii
# Escapes chars (in second string); also escapes all non-ASCII chars.
function escape_with(str, use)
    chars = split(use, "")

    for c in chars
        _char = c[1] # Character string as Char
        h = hex(int(_char))
        if length(h) < 2
            h = "0"*h
        end
        str = replace(str, c, "%" * h)
    end

    for i in nonascii_array
        str = replace(str, char(i), "%" * hex(i))
    end

    return str
end

function escape(str)
    return escape_with(str, unescaped)
end

function escape_form(str)
    str = escape_with(str, unescaped_form)
    return replace(str, " ", "+")
end


text:
import Base: iterate, elsize, size, getindex, setindex!, parent

abstract type AbstractDiagonalIterator{T, V}  <: AbstractVector{T} end

"""
    DiagonalIterator{T <: AbstractMatrix}

Iterator over the diagonal of an AbstractMatrix
"""
struct DiagonalIterator{T, V <:AbstractMatrix} <: AbstractDiagonalIterator{T, V}
    M::V
    length::Int
end

"""
    altdiagonal(M::AbstractMatrix)

An iterator over the diagonal of `M`. `altdiagonal` is an alternative
to the iterator `diagonal`.  An iterator is constructed much more quickly
with `altdiagonal` than with `diagonal`.
"""
altdiagonal(M::V) where V <: AbstractMatrix  = DiagonalIterator{eltype(M), V}(M, checksquare(M))

iterate(iter::DiagonalIterator, i = 1) = i > iter.length ? nothing : (iter.M[i,i], i + 1)
size(iter::DiagonalIterator) = (iter.length,)
elsize(::Type{DiagonalIterator{T, V}}) where {T, V} = sizeof(T)
getindex(iter::DiagonalIterator, i) = iter.M[i,i]
setindex!(iter::DiagonalIterator, v, i) = (iter.M[i,i] = v)
parent(iter::DiagonalIterator) = iter.M
Base.dataids(d::DiagonalIterator) = Base.dataids(parent(d))


text:
module ElementwiseTests

using ReverseDiff

using DiffRules
using DiffTests
using ForwardDiff
using LogExpFunctions
using NaNMath
using SpecialFunctions

using Test

include(joinpath(dirname(@__FILE__), "../utils.jl"))

x, y = rand(3, 3), rand(3, 3)
a, b = rand(3), rand(3)
n = rand()
tp = InstructionTape()
offset = 0.00001

function test_elementwise(f, fopt, x, tp)
    xt = track(copy(x), tp)

    ########################################

    y = map(f, x)

    # record
    yt = map(fopt, xt)
    @test yt == y
    @test length(tp) == 1

    # reverse
    out = similar(y, (length(x), length(x)))
    ReverseDiff.seeded_reverse_pass!(out, yt, xt, tp)
    test_approx(out, ForwardDiff.jacobian(z -> map(f, z), x))

    # forward
    x2 = x .- offset
    ReverseDiff.value!(xt, x2)
    ReverseDiff.forward_pass!(tp)
    @test yt == map(f, x2)
    ReverseDiff.value!(xt, x)

    empty!(tp)

    ########################################

    y = broadcast(f, x)

    # record
    yt = broadcast(fopt, xt)
    @test yt == y
    @test length(tp) == 1

    # reverse
    out = similar(y, (length(x), length(x)))
    ReverseDiff.seeded_reverse_pass!(out, yt, xt, tp)
    test_approx(out, ForwardDiff.jacobian(z -> broadcast(f, z), x))

    # forward
    x2 = x .- offset
    ReverseDiff.value!(xt, x2)
    ReverseDiff.forward_pass!(tp)
    @test yt == broadcast(f, x2)
    ReverseDiff.value!(xt, x)

    empty!(tp)
end

function test_map(f, fopt, a, b, tp)
    at, bt = track(copy(a), tp), track(copy(b), tp)
    c = map(f, a, b)

    ########################################

    # record
    ct = map(fopt, at, b)
    @test ct == c
    @test length(tp) == 1

    # reverse
    out = similar(c, (length(a), length(a)))
    ReverseDiff.seeded_reverse_pass!(out, ct, at, tp)
    test_approx(out, ForwardDiff.jacobian(x -> map(f, x, b), a))

    # forward
    a2 = a .- offset
    ReverseDiff.value!(at, a2)
    ReverseDiff.forward_pass!(tp)
    @test ct == map(f, a2, b)
    ReverseDiff.value!(at, a)

    empty!(tp)

    ########################################

    # record
    ct = map(fopt, a, bt)
    @test ct == c
    @test length(tp) == 1

    # reverse
    out = similar(c, (length(a), length(a)))
    ReverseDiff.seeded_reverse_pass!(out, ct, bt, tp)
    test_approx(out, ForwardDiff.jacobian(x -> map(f, a, x), b))

    # forward
    b2 = b .- offset
    ReverseDiff.value!(bt, b2)
    ReverseDiff.forward_pass!(tp)
    @test ct == map(f, a, b2)
    ReverseDiff.value!(bt, b)

    empty!(tp)

    ########################################

    # record
    ct = map(fopt, at, bt)
    @test ct == c
    @test length(tp) == 1

    # reverse
    out_a = similar(c, (length(a), length(a)))
    out_b = similar(c, (length(a), length(a)))
    ReverseDiff.seeded_reverse_pass!(out_a, ct, at, tp)
    ReverseDiff.seeded_reverse_pass!(out_b, ct, bt, tp)
    test_approx(out_a, ForwardDiff.jacobian(x -> map(f, x, b), a))
    test_approx(out_b, ForwardDiff.jacobian(x -> map(f, a, x), b))

    # forward
    a2, b2 = a .- offset, b .- offset
    ReverseDiff.value!(at, a2)
    ReverseDiff.value!(bt, b2)
    ReverseDiff.forward_pass!(tp)
    @test ct == map(f, a2, b2)
    ReverseDiff.value!(at, a)
    ReverseDiff.value!(bt, b)

    empty!(tp)
end

function test_broadcast(f, fopt, a::AbstractArray, b::AbstractArray, tp, builtin::Bool = false)
    at, bt = track(copy(a), tp), track(copy(b), tp)

    if builtin
        g = fopt
    else
        g = (x, y) -> broadcast(fopt, x, y)
    end

    c = g(a, b)

    ########################################

    # record
    ct = g(at, b)
    @test ct == c
    @test length(tp) == 1

    # reverse
    out = similar(c, (length(c), length(a)))
    ReverseDiff.seeded_reverse_pass!(out, ct, at, tp)
    test_approx(out, ForwardDiff.jacobian(x -> g(x, b), a))

    # forward
    a2 = a .- offset
    ReverseDiff.value!(at, a2)
    ReverseDiff.forward_pass!(tp)
    @test ct == g(a2, b)
    ReverseDiff.value!(at, a)

    empty!(tp)

    ########################################

    # record
    ct = g(a, bt)
    @test ct == c
    @test length(tp) == 1

    # reverse
    out = similar(c, (length(c), length(b)))
    ReverseDiff.seeded_reverse_pass!(out, ct, bt, tp)
    test_approx(out, ForwardDiff.jacobian(x -> g(a, x), b))

    # forward
    b2 = b .- offset
    ReverseDiff.value!(bt, b2)
    ReverseDiff.forward_pass!(tp)
    @test ct == g(a, b2)
    ReverseDiff.value!(bt, b)

    empty!(tp)

    ########################################

    # record
    ct = g(at, bt)
    @test ct == c
    @test length(tp) == 1

    # reverse
    out_a = similar(c, (length(c), length(a)))
    out_b = similar(c, (length(c), length(b)))
    ReverseDiff.seeded_reverse_pass!(out_a, ct, at, tp)
    ReverseDiff.seeded_reverse_pass!(out_b, ct, bt, tp)
    test_approx(out_a, ForwardDiff.jacobian(x -> g(x, b), a))
    test_approx(out_b, ForwardDiff.jacobian(x -> g(a, x), b))

    # forward
    a2, b2 = a .- offset, b .- offset
    ReverseDiff.value!(at, a2)
    ReverseDiff.value!(bt, b2)
    ReverseDiff.forward_pass!(tp)
    @test ct == g(a2, b2)
    ReverseDiff.value!(at, a)
    ReverseDiff.value!(bt, b)

    empty!(tp)
end

function test_broadcast(f, fopt, n::Number, x::AbstractArray, tp, builtin::Bool = false)
    nt, xt = track(copy(n), tp), track(copy(x), tp)

    if builtin
        g = fopt
    else
        g = (x, y) -> broadcast(fopt, x, y)
    end

    y = g(n, x)

    ########################################

    # record
    yt = g(nt, x)
    @test yt == y
    @test length(tp) == 1

    # reverse
    out = similar(y)
    ReverseDiff.seeded_reverse_pass!(out, yt, nt, tp)
    test_approx(out, ForwardDiff.derivative(z -> g(z, x), n))

    # forward
    n2 = n + offset
    ReverseDiff.value!(nt, n2)
    ReverseDiff.forward_pass!(tp)
    @test yt == g(n2, x)
    ReverseDiff.value!(nt, n)

    empty!(tp)

    ########################################

    # record
    yt = g(n, xt)
    @test yt == y
    @test length(tp) == 1

    # reverse
    out = similar(y, (length(y), length(x)))
    ReverseDiff.seeded_reverse_pass!(out, yt, xt, tp)
    test_approx(out, ForwardDiff.jacobian(z -> g(n, z), x))

    # forward
    x2 = x .- offset
    ReverseDiff.value!(xt, x2)
    ReverseDiff.forward_pass!(tp)
    @test yt == g(n, x2)
    ReverseDiff.value!(xt, x)

    empty!(tp)

    ########################################

    # record
    yt = g(nt, xt)
    @test yt == y
    @test length(tp) == 1

    # reverse
    out_n = similar(y)
    out_x = similar(y, (length(y), length(x)))
    ReverseDiff.seeded_reverse_pass!(out_n, yt, nt, tp)
    ReverseDiff.seeded_reverse_pass!(out_x, yt, xt, tp)
    test_approx(out_n, ForwardDiff.derivative(z -> g(z, x), n))
    test_approx(out_x, ForwardDiff.jacobian(z -> g(n, z), x))

    # forward
    n2, x2 = n + offset , x .- offset
    ReverseDiff.value!(nt, n2)
    ReverseDiff.value!(xt, x2)
    ReverseDiff.forward_pass!(tp)
    @test yt == g(n2, x2)
    ReverseDiff.value!(nt, n)
    ReverseDiff.value!(xt, x)

    empty!(tp)
end

function test_broadcast(f, fopt, x::AbstractArray, n::Number, tp, builtin::Bool = false)
    xt, nt = track(copy(x), tp), track(copy(n), tp)

    if builtin
        g = fopt
    else
        g = (x, y) -> broadcast(fopt, x, y)
    end

    y = g(x, n)

    ########################################

    # record
    yt = g(x, nt)
    @test yt == y
    @test length(tp) == 1

    # reverse
    out = similar(y)
    ReverseDiff.seeded_reverse_pass!(out, yt, nt, tp)
    test_approx(out, ForwardDiff.derivative(z -> g(x, z), n))

    # forward
    n2 = n + offset
    ReverseDiff.value!(nt, n2)
    ReverseDiff.forward_pass!(tp)
    @test yt == g(x, n2)
    ReverseDiff.value!(nt, n)

    empty!(tp)

    ########################################

    # record
    yt = g(xt, n)
    @test yt == y
    @test length(tp) == 1

    # reverse
    out = similar(y, (length(y), length(x)))
    ReverseDiff.seeded_reverse_pass!(out, yt, xt, tp)
    test_approx(out, ForwardDiff.jacobian(z -> g(z, n), x))

    # forward
    x2 = x .- offset
    ReverseDiff.value!(xt, x2)
    ReverseDiff.forward_pass!(tp)
    @test yt == g(x2, n)
    ReverseDiff.value!(xt, x)

    empty!(tp)

    ########################################

    # record
    yt = g(xt, nt)
    @test yt == y
    @test length(tp) == 1

    # reverse
    out_n = similar(y)
    out_x = similar(y, (length(y), length(x)))
    ReverseDiff.seeded_reverse_pass!(out_n, yt, nt, tp)
    ReverseDiff.seeded_reverse_pass!(out_x, yt, xt, tp)
    test_approx(out_n, ForwardDiff.derivative(z -> g(x, z), n))
    test_approx(out_x, ForwardDiff.jacobian(z -> g(z, n), x))

    # forward
    x2, n2 = x .- offset, n + offset
    ReverseDiff.value!(xt, x2)
    ReverseDiff.value!(nt, n2)
    ReverseDiff.forward_pass!(tp)
    @test yt == g(x2, n2)
    ReverseDiff.value!(xt, x)
    ReverseDiff.value!(nt, n)

    empty!(tp)
end

for f in DiffTests.NUMBER_TO_NUMBER_FUNCS
    test_println("DiffTests.NUMBER_TO_NUMBER_FUNCS", f)
    test_elementwise(f, ReverseDiff.@forward(f), x, tp)
    test_elementwise(f, ReverseDiff.@forward(f), a, tp)
end

for (M, fsym, arity) in DiffRules.diffrules(; filter_modules=nothing)
    # ensure that all rules can be tested
    if !(isdefined(@__MODULE__, M) && isdefined(getfield(@__MODULE__, M), fsym))
        error("$M.$fsym is not available")
    end
    fsym === :rem2pi && continue
    if arity == 1
        f = eval(:($M.$fsym))
        test_println("forward-mode unary scalar functions", f)
        test_elementwise(f, f, modify_input(fsym, x), tp)
        test_elementwise(f, f, modify_input(fsym, a), tp)
    elseif arity == 2
        in(fsym, SKIPPED_BINARY_SCALAR_TESTS) && continue
        f = eval(:($M.$fsym))
        test_println("forward-mode binary scalar functions", f)
        test_map(f, f, x, y, tp)
        test_map(f, f, a, b, tp)
        test_broadcast(f, f, x, y, tp)
        test_broadcast(f, f, a, b, tp)
        test_broadcast(f, f, x, a, tp)
        test_broadcast(f, f, a, x, tp)
        test_broadcast(f, f, n, x, tp)
        test_broadcast(f, f, x, n, tp)
        test_broadcast(f, f, n, a, tp)
        test_broadcast(f, f, a, n, tp)
    end
end

for f in DiffTests.BINARY_BROADCAST_OPS
    test_println("built-in broadcast operators", f)
    test_broadcast(f, f, x, y, tp, true)
    test_broadcast(f, f, a, b, tp, true)
    test_broadcast(f, f, x, a, tp, true)
    test_broadcast(f, f, a, x, tp, true)
    test_broadcast(f, f, n, x, tp, true)
    test_broadcast(f, f, x, n, tp, true)
    test_broadcast(f, f, n, a, tp, true)
    test_broadcast(f, f, a, n, tp, true)
end

end # module

max_stars_repo_name: jrevels/ReverseDiffPrototype

max_issues_repo_name: jrevels/ReverseDiffPrototype

max_forks_repo_name: jrevels/ReverseDiffPrototype


text:
mutable struct CompDotAcc{T}
    s :: T
    e :: T
end

compDotAcc(T) = CompDotAcc{T}(vzero(T), vzero(T))

@inline function add!(acc::CompDotAcc{T}, x::T, y::T) where {T}
    SIMDops.@explicit

    p, ep = two_prod(x, y)
    acc.s, es = two_sum(acc.s, p)
    acc.e += ep + es
end

@inline function add!(acc::A, x::A) where {A<:CompDotAcc}
    SIMDops.@explicit

    acc.s, e = two_sum(acc.s, x.s)
    acc.e += x.e + e
end

@inline function Base.sum(acc::CompDotAcc{T}) where {T<:Vec}
    acc_r = compDotAcc(fptype(T))
    acc_r.e = vsum(acc.e)
    for xi in acc.s
        acc_r.s, ei = two_sum(acc_r.s, xi.value)
        acc_r.e += ei
    end
    acc_r
end

@inline function Base.sum(acc::CompDotAcc{T}) where {T<:Real}
    acc.s + acc.e
end