text:
using JuMP, LinearAlgebra, Suppressor, StatsBase, CSV, DataFrames, Compat, Test, Random
include("exact_approaches.jl")
include("alt_min.jl")
include("convexpenalties.jl")

results_template = DataFrame(n=Int[], r=Int[], p=Real[], err_NN=Real[], rnk_NN=Int[],
solve_NN=Real[], err_GD=Real[], rnk_GD=Int[], solve_GD=Real[], err_AM=Real[], rnk_AM=Int[], solve_AM=Real[], err_OA=Real[], rnk_OA=Int[], solve_OA=Real[], nc_OA=Int[], gap_OA=Real[], cuts_OA=Int[])

results_exact=similar(results_template, 0)

for ARG in ARGS
    theArg=parse(Int, ARG)
    pvals=[0.0125, 0.025, 0.0375, 0.05, 0.0625, 0.075, 0.1, 0.125, 0.15, 0.175, 0.2, 0.225, 0.25, 0.275, 0.3, 0.325, 0.35, 0.375, 0.4]
    p=pvals[theArg]
    for iter in 1:1

        n=50
        r=1
        gamma=100/p
        M_L=randn(n, r)
        M_R=randn(n, r)
        M=M_L*M_R'

        @show maximum(M)
        numEntriesSampled=Int(round(p*n^2))
        coords=[repeat(1:n, inner=size(1:n,1)) repeat(1:n, outer=size(1:n,1))];
        entries=sample(1:n^2, Int(numEntriesSampled), replace=false)

        theIndices=zeros(n,n)
        for t=1:Int(numEntriesSampled)
            theIndices[coords[entries[t],1], coords[entries[t],2]]=1.0
        end

        k=r
        mu=1.0 #imposing constraints so mu doesn't matter
        t_NN=@elapsed X_NN=getMatrixNuclearNorm(coords, entries, M, n, numEntriesSampled, mu, true) #imposing constraints as noiseless so we want exact signal recovery asap
        u,sigma,v,=svd(X_NN)
        rnk_NN=sum(sigma.>1e-2);
        err_NN = norm(X_NN.-M)/norm(M);
        # Get SDP lower bound (for sake of comparison)
        t_GD=@elapsed X_GD,LB, UB, Y_relax=getMatrixFrobeniusNorm_greedy(coords, entries, M, n, numEntriesSampled, k, gamma, false)

        @show "SDP lower bound is: " LB
        u,sigma,v,=svd(X_GD)
        rnk_GD=sum(sigma.>1e-2);
        err_GD = norm(X_GD.-M)/norm(M);

        t_AM=@elapsed X_AM,=alternatingminimization_matrixcompletion(coords, entries, M, n, numEntriesSampled, k, gamma, 50, X_GD)#AM warm started with GD

        # Test solution on overall data.
        u,sigma,v,=svd(X_AM)
        @show rnk_AM=sum(sigma.>1e-2);
        @show err_AM = norm(X_AM.-M)/norm(M);

        u_0=u[:,1:k]
        v_0=(Diagonal(sigma[1:k])*v[:,1:k]')'

        # Apply the cutting-plane method (multi-tree version)
        t_OA=@elapsed X_OA, nc_OA, cuts_OA, gap_OA=getMatrixFrobeniusNorm_cuttingplanesmultitree_exact_inout(coords, entries, M, n, numEntriesSampled, k, gamma, Y_relax, X_AM, true)
        # Test solution on overall data.
        u,v,=svd(X_OA)
        @show rnk_OA=sum(v.>1e-2);
        @show err_OA = norm(X_OA.-M)/norm(M);

        push!(results_exact, [n, r, p, err_NN, rnk_NN, t_NN, err_GD, rnk_GD, t_GD, err_AM, rnk_AM, t_AM, err_OA, rnk_OA, t_OA, nc_OA, gap_OA, cuts_OA])

        end
    end

CSV.write("results_basispursuit_varyingp_n50_k1_gamma100onp.csv", results_exact, append=true)

max_stars_repo_path: Matrix Completion/createFigure4b.jl

max_stars_repo_name: ryancorywright/MixedProjectionSoftware

max_issues_repo_path: Matrix Completion/createFigure4b.jl

max_issues_repo_name: ryancorywright/MixedProjectionSoftware

max_forks_repo_path: Matrix Completion/createFigure4b.jl

max_forks_repo_name: ryancorywright/MixedProjectionSoftware


text:
function fieldtype(t, f) nothing end

"""
    Genser.gensertypefor(T::Type) :: Type{<: GenserDataType}

Get the Genser data model type for a Julia type.
"""
gensertypefor(T::Type{<:GenserDataType}; kwargs...) = T
gensertypefor(::Type{Any}; kwargs...) = GenserAny
gensertypefor(::Type{Nothing}; kwargs...) = GenserNull
gensertypefor(::Type{Int8}; kwargs...) = GenserInt8
gensertypefor(::Type{UInt8}; kwargs...) = GenserUInt8
gensertypefor(::Type{Int16}; kwargs...) = GenserInt16
gensertypefor(::Type{UInt16}; kwargs...) = GenserUInt16
gensertypefor(::Type{Int32}; kwargs...) = GenserInt32
gensertypefor(::Type{UInt32}; kwargs...) = GenserUInt32
gensertypefor(::Type{Int64}; kwargs...) = GenserInt64
gensertypefor(::Type{UInt64}; kwargs...) = GenserUInt64
gensertypefor(::Type{Int128}; kwargs...) = GenserInt128
gensertypefor(::Type{UInt128}; kwargs...) = GenserUInt128
gensertypefor(::Type{Char}; kwargs...) = GenserChar
gensertypefor(::Type{Bool}; kwargs...) = GenserBool
gensertypefor(::Type{BigInt}; kwargs...) = GenserBigInt
gensertypefor(::Type{Float16}; kwargs...) = GenserFloat16
gensertypefor(::Type{Float32}; kwargs...) = GenserFloat32
gensertypefor(::Type{Float64}; kwargs...) = GenserFloat64
gensertypefor(::Type{BigFloat}; kwargs...) = GenserBigFloat
gensertypefor(::Type{Rational}; kwargs...) where V = GenserRational
gensertypefor(::Type{Rational{V}}; kwargs...) where V = GenserRational

gensertypefor(::Type{String}; kwargs...) = GenserString
gensertypefor(T::Type{<: AbstractString}; kwargs...) = GenserString
gensertypefor(::Type{Vector{UInt8}}; kwargs...) = begin
    if haskey(kwargs, :containertype) && haskey(kwargs, :fieldkey)
        containertype = kwargs[:containertype]
        fieldname = kwargs[:fieldkey]
        encoding = fieldencoding(containertype, Val{fieldname})
        if encoding <: Encoding
            GenserBinaryValue{encoding}
        else
            GenserBinary
        end
    else
        GenserBinary
    end
end
gensertypefor(::Type{UUID}; kwargs...) = GenserUUID
gensertypefor(::Type{Symbol}; kwargs...) = GenserSymbol
gensertypefor(::Type{<: AbstractArray{V}}; kwargs...) where V = GenserSequence{gensertypefor(V)}
gensertypefor(::Type{<: AbstractSet{V}}; kwargs...) where V = GenserSet{gensertypefor(V)}
gensertypefor(T::Type{<: Tuple}; kwargs...) = begin
    types = Base.map(gensertypefor, [T.parameters...])
    GenserTuple{Tuple{types...}}
end
gensertypefor(::Type{<: AbstractDict{K,V}}; kwargs...) where K where V = GenserDict{gensertypefor(K),gensertypefor(V)}
gensertypefor(T::Type{<: NamedTuple}; kwargs...) = begin
    keys = fieldnames(T)
    vals = Base.map(gensertypefor, [T.types...])
    vals = Tuple{vals...}
    NewT = NamedTuple{keys, vals}
    GenserRecord{NewT}
end

function gensertypefor(T::Type; kwargs...)
    # Could be a struct or union
    ttype = typeof(T)

    if ttype == Union
        if Nothing <: T
            # optional type
            return GenserOptional{gensertypefor(unoptionalize(T); kwargs...)}
        end

        # variant type
        return GenserVariant{mapuniontype(T; kwargs...)}
    elseif ttype === DataType
        if T <: Enum
            return GenserString
        end
        # Assume a struct
        keys = fieldnames(T)
        vals = Base.map(zip(keys, T.types)) do (fieldkey, type)
            definedtype = fieldtype(T, Val{fieldkey})
            if definedtype isa DataType && definedtype <: GenserDataType && definedtype.isconcretetype
                definedtype
            else
                gensertypefor(type, fieldkey=fieldkey, containertype=T)
            end
        end
        vals = Tuple{vals...}
        NewT = NamedTuple{keys, vals}
        return GenserRecord{NewT}
    else
        throw(ArgumentError("unable to derive genser type for $T"))
    end
end

function unoptionalize(T::Type)
    @assert typeof(T) === Union
    @assert Nothing <: T
    subtypes = []
    while T.a !== Nothing
        push!(subtypes, T.a)
        T = T.b
    end
    push!(subtypes, T.b)

    return Union{subtypes...}
end

function mapuniontype(T::Type; kwargs...)
    subtypes = []
    while typeof(T) === Union
        push!(subtypes, gensertypefor(T.a; kwargs...))
        T = T.b
    end
    push!(subtypes, gensertypefor(T; kwargs...))

    Union{subtypes...}
end

ext: jl


text:
using Gadfly, RDatasets

set_default_plot_size(6inch, 2inch)

plot(dataset("Zelig", "macro"), x=:Year, y=:Unem, color=:Country, Geom.point)


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

@nospecialize

if Pair != Base.Pair
import Base: Base, IOContext, string, join, sprint
IOContext(io::IO, KV::Pair) = IOContext(io, Base.Pair(KV[1], KV[2]))
length(s::String) = Base.length(s)
^(s::String, i::Int) = Base.:^(s, i)
end

import Base: show_unquoted
using Base: printstyled, with_output_color, prec_decl

function Base.show(io::IO, cfg::CFG)
    for (idx, block) in enumerate(cfg.blocks)
        print(io, idx, "\t=>\t")
        join(io, block.succs, ", ")
        println(io)
    end
end

function print_stmt(io::IO, idx::Int, @nospecialize(stmt), used::BitSet, maxlength_idx::Int, color::Bool, show_type::Bool)
    indent = maxlength_idx + 4
    if idx in used
        pad = " "^(maxlength_idx - length(string(idx)) + 1)
        print(io, "%", idx, pad, "= ")
    else
        print(io, " "^indent)
    end
    if !color && stmt isa PiNode
        # when the outer context is already colored (yellow, for pending nodes), don't use the usual coloring printer
        print(io, "π (")
        show_unquoted(io, stmt.val, indent)
        print(io, ", ")
        print(io, stmt.typ)
        print(io, ")")
    elseif isexpr(stmt, :invoke)
        # TODO: why is this here, and not in Base.show_unquoted
        print(io, "invoke ")
        linfo = stmt.args[1]
        show_unquoted(io, stmt.args[2], indent)
        print(io, "(")
        # XXX: this is wrong if `sig` is not a concretetype method
        # more correct would be to use `fieldtype(sig, i)`, but that would obscure / discard Varargs information in show
        sig = linfo.specTypes == Tuple ? Core.svec() : Base.unwrap_unionall(linfo.specTypes).parameters::Core.SimpleVector
        print_arg(i) = sprint() do io
            show_unquoted(io, stmt.args[i], indent)
            if (i - 1) <= length(sig)
                print(io, "::", sig[i - 1])
            end
        end
        join(io, (print_arg(i) for i = 3:length(stmt.args)), ", ")
        print(io, ")")
    # given control flow information, we prefer to print these with the basic block #, instead of the ssa %
    elseif isexpr(stmt, :enter) && length(stmt.args) == 1 && stmt.args[1] isa Int
        print(io, "\$(Expr(:enter, #", stmt.args[1]::Int, "))")
    elseif stmt isa GotoNode
        print(io, "goto #", stmt.label)
    elseif stmt isa PhiNode
        show_unquoted_phinode(io, stmt, indent, "#")
    elseif stmt isa GotoIfNot
        show_unquoted_gotoifnot(io, stmt, indent, "#")
    # everything else in the IR, defer to the generic AST printer
    else
        show_unquoted(io, stmt, indent, show_type ? prec_decl : 0)
    end
    nothing
end

show_unquoted(io::IO, val::Argument, indent::Int, prec::Int) = show_unquoted(io, Core.SlotNumber(val.n), indent, prec)

show_unquoted(io::IO, stmt::PhiNode, indent::Int, ::Int) = show_unquoted_phinode(io, stmt, indent, "%")
function show_unquoted_phinode(io::IO, stmt::PhiNode, indent::Int, prefix::String)
    args = map(1:length(stmt.edges)) do i
        e = stmt.edges[i]
        v = !isassigned(stmt.values, i) ? "#undef" :
            sprint() do io′
                show_unquoted(io′, stmt.values[i], indent)
            end
        return "$prefix$e => $v"
    end
    print(io, "φ ", '(')
    join(io, args, ", ")
    print(io, ')')
end

function show_unquoted(io::IO, stmt::PhiCNode, indent::Int, ::Int)
    print(io, "φᶜ (")
    first = true
    for v in stmt.values
        first ? (first = false) : print(io, ", ")
        show_unquoted(io, v, indent)
    end
    print(io, ")")
end

function show_unquoted(io::IO, stmt::PiNode, indent::Int, ::Int)
    print(io, "π (")
    show_unquoted(io, stmt.val, indent)
    print(io, ", ")
    printstyled(io, stmt.typ, color=:cyan)
    print(io, ")")
end

function show_unquoted(io::IO, stmt::UpsilonNode, indent::Int, ::Int)
    print(io, "ϒ (")
    isdefined(stmt, :val) ?
        show_unquoted(io, stmt.val, indent) :
        print(io, "#undef")
    print(io, ")")
end

function show_unquoted(io::IO, stmt::ReturnNode, indent::Int, ::Int)
    if !isdefined(stmt, :val)
        print(io, "unreachable")
    else
        print(io, "return ")
        show_unquoted(io, stmt.val, indent)
    end
end

show_unquoted(io::IO, stmt::GotoIfNot, indent::Int, ::Int) = show_unquoted_gotoifnot(io, stmt, indent, "%")
function show_unquoted_gotoifnot(io::IO, stmt::GotoIfNot, indent::Int, prefix::String)
    print(io, "goto ", prefix, stmt.dest, " if not ")
    show_unquoted(io, stmt.cond, indent)
end

function compute_inlining_depth(linetable::Vector, iline::Int32)
    iline == 0 && return 1
    depth = -1
    while iline != 0
        depth += 1
        lineinfo = linetable[iline]::LineInfoNode
        iline = lineinfo.inlined_at
    end
    return depth
end

function should_print_ssa_type(@nospecialize node)
    if isa(node, Expr)
        return !(node.head in (:gc_preserve_begin, :gc_preserve_end, :meta, :return, :enter, :leave))
    end
    return !isa(node, PiNode)   && !isa(node, GotoIfNot) &&
           !isa(node, GotoNode) && !isa(node, ReturnNode) &&
           !isa(node, QuoteNode)
end

function default_expr_type_printer(io::IO, @nospecialize(typ), used::Bool)
    printstyled(io, "::", typ, color=(used ? :cyan : :light_black))
    nothing
end

# converts the linetable for line numbers
# into a list in the form:
#   1 outer-most-frame
#   2   inlined-frame
#   3     innermost-frame
function compute_loc_stack(linetable::Vector, line::Int32)
    stack = Int[]
    while line != 0
        entry = linetable[line]::LineInfoNode
        pushfirst!(stack, line)
        line = entry.inlined_at
    end
    return stack
end

"""
    Compute line number annotations for an IRCode

This functions compute three sets of annotations for each IR line. Take the following
example (taken from `@code_typed sin(1.0)`):

```
    **                                                    ***         **********
    35 6 ── %10  = :(Base.mul_float)(%%2, %%2)::Float64   │╻╷         sin_kernel
       │    %11  = :(Base.mul_float)(%10, %10)::Float64   ││╻          *
```

The three annotations are indicated with `*`. The first one is the line number of the
active function (printed once whenver the outer most line number changes). The second
is the inlining indicator. The number of lines indicate the level of nesting, with a
half-size line (╷) indicating the start of a scope and a full size line (│) indicating
a continuing scope. The last annotation is the most complicated one. It is a heuristic
way to print the name of the entered scope. What it attempts to do is print the outermost
scope that hasn't been printed before. Let's work a number of examples to see the impacts
and tradeoffs involved.

```
f() = leaf_function() # Delibarately not defined to end up in the IR verbatim
g() = f()
h() = g()
top_function() = h()
```

After inlining, we end up with:
```
1 1 ─ %1 = :(Main.leaf_function)()::Any   │╻╷╷ h
  └──      return %1                      │
```

We see that the only function printed is the outermost function. This certainly loses
some information, but the idea is that the outermost function would have the most
semantic meaning (in the context of the function we're looking at).

On the other hand, let's see what happens when we redefine f:
```
function f()
    leaf_function()
    leaf_function()
    leaf_function()
end
```

We get:
```
1 1 ─      :(Main.leaf_function)()::Any   │╻╷╷ h
  │        :(Main.leaf_function)()::Any   ││┃│  g
  │   %3 = :(Main.leaf_function)()::Any   │││┃   f
  └──      return %3                      │
```

Even though we were in the `f` scope since the first statement, it tooks us two statements
to catch up and print the intermediate scopes. Which scope is printed is indicated both
by the indentation of the method name and by an increased thickness of the appropriate
line for the scope.
"""
function compute_ir_line_annotations(code::Union{IRCode, CodeInfo})
    loc_annotations = String[]
    loc_methods = String[]
    loc_lineno = String[]
    cur_group = 1
    last_line = 0
    last_lineno = 0
    last_stack = []
    last_printed_depth = 0
    stmts = (code isa IRCode ? code.stmts : code.code)
    linetable = code.linetable
    lines = (code isa IRCode ? code.lines : code.codelocs)
    for idx in eachindex(stmts)
        buf = IOBuffer()
        line = lines[idx]
        depth = compute_inlining_depth(linetable, line)
        iline = line
        lineno = 0
        loc_method = ""
        print(buf, "│")
        if line != 0
            stack = compute_loc_stack(linetable, line)
            lineno = linetable[stack[1]].line
            x = min(length(last_stack), length(stack))
            if length(stack) != 0
                # Compute the last depth that was in common
                first_mismatch = findfirst(i->last_stack[i] != stack[i], 1:x)
                # If the first mismatch is the last stack frame, that might just
                # be a line number mismatch in inner most frame. Ignore those
                if length(last_stack) == length(stack) && first_mismatch == length(stack)
                    last_entry, entry = linetable[last_stack[end]], linetable[stack[end]]
                    if last_entry.method == entry.method && last_entry.file == entry.file
                        first_mismatch = nothing
                    end
                end
                last_depth = something(first_mismatch, x+1)-1
                if min(depth, last_depth) > last_printed_depth
                    printing_depth = min(depth, last_printed_depth + 1)
                    last_printed_depth = printing_depth
                elseif length(stack) > length(last_stack) || first_mismatch != nothing
                    printing_depth = min(depth, last_depth + 1)
                    last_printed_depth = printing_depth
                else
                    printing_depth = 0
                end
                stole_one = false
                if printing_depth != 0
                    for _ in 1:(printing_depth-1)
                        print(buf, "│")
                    end
                    if printing_depth <= last_depth-1 && first_mismatch === nothing
                        print(buf, "┃")
                        for _ in printing_depth+1:min(depth, last_depth)
                            print(buf, "│")
                        end
                    else
                        stole_one = true
                        print(buf, "╻")
                    end
                else
                    for _ in 1:min(depth, last_depth)
                        print(buf, "│")
                    end
                end
                print(buf, "╷"^max(0,depth-last_depth-stole_one))
                if printing_depth != 0
                    if length(stack) == printing_depth
                        loc_method = String(linetable[line].method)
                    else
                        loc_method = String(linetable[stack[printing_depth+1]].method)
                    end
                end
                loc_method = string(" "^printing_depth, loc_method)
            end
            last_stack = stack
            entry = linetable[line]
        end
        push!(loc_annotations, String(take!(buf)))
        push!(loc_lineno, (lineno != 0 && lineno != last_lineno) ? string(lineno) : "")
        push!(loc_methods, loc_method)
        last_line = line
        (lineno != 0) && (last_lineno = lineno)
    end
    return (loc_annotations, loc_methods, loc_lineno)
end

Base.show(io::IO, code::IRCode) = show_ir(io, code)

function show_ir(io::IO, code::IRCode, expr_type_printer=default_expr_type_printer; verbose_linetable=false)
    cols = displaysize(io)[2]
    used = BitSet()
    stmts = code.stmts
    types = code.types
    cfg = code.cfg
    max_bb_idx_size = length(string(length(cfg.blocks)))
    for stmt in stmts
        scan_ssa_use!(push!, used, stmt)
    end
    bb_idx = 1
    new_nodes = code.new_nodes
    if any(i -> !isassigned(code.new_nodes, i), 1:length(code.new_nodes))
        printstyled(io, :red, "ERROR: New node array has unset entry\n")
        new_nodes = new_nodes[filter(i -> isassigned(code.new_nodes, i), 1:length(code.new_nodes))]
    end
    for nn in new_nodes
        scan_ssa_use!(push!, used, nn.node)
    end
    perm = sortperm(new_nodes, by = x->x.pos)
    new_nodes_perm = Iterators.Stateful(perm)

    if isempty(used)
        maxlength_idx = 0
    else
        maxused = maximum(used)
        maxlength_idx = length(string(maxused))
    end
    if !verbose_linetable
        (loc_annotations, loc_methods, loc_lineno) = compute_ir_line_annotations(code)
        max_loc_width = maximum(length(str) for str in loc_annotations)
        max_lineno_width = maximum(length(str) for str in loc_lineno)
        max_method_width = maximum(length(str) for str in loc_methods)
    end
    max_depth = maximum(compute_inlining_depth(code.linetable, line) for line in code.lines)
    last_stack = []
    for idx in eachindex(stmts)
        if !isassigned(stmts, idx)
            # This is invalid, but do something useful rather
            # than erroring, to make debugging easier
            printstyled(io, :red, "#UNDEF\n")
            continue
        end
        stmt = stmts[idx]
        # Compute BB guard rail
        bbrange = cfg.blocks[bb_idx].stmts
        bbrange = bbrange.first:bbrange.last
        bb_idx_str = string(bb_idx)
        bb_pad = max_bb_idx_size - length(bb_idx_str)
        bb_type = length(cfg.blocks[bb_idx].preds) <= 1 ? "─" : "┄"
        bb_start_str = string(bb_idx_str, " ", bb_type, "─"^bb_pad, " ")
        bb_guard_rail_cont = string("│  ", " "^max_bb_idx_size)
        if idx == first(bbrange)
            bb_guard_rail = bb_start_str
        else
            bb_guard_rail = bb_guard_rail_cont
        end
        # Print linetable information
        if verbose_linetable
            stack = compute_loc_stack(code.linetable, code.lines[idx])
            # We need to print any stack frames that did not exist in the last stack
            ndepth = max(1, length(stack))
            rail = string(" "^(max_depth+1-ndepth), "│"^ndepth)
            start_column = cols - max_depth - 10
            for (i, x) in enumerate(stack)
                if i > length(last_stack) || last_stack[i] != x
                    entry = code.linetable[x]
                    printstyled(io, "\e[$(start_column)G$(rail)\e[1G", color = :light_black)
                    print(io, bb_guard_rail)
                    ssa_guard = " "^(maxlength_idx + 4 + (i - 1))
                    entry_label = "$(ssa_guard)$(entry.method) at $(entry.file):$(entry.line) "
                    hline = string("─"^(start_column-length(entry_label)-length(bb_guard_rail)+max_depth-i), "┐")
                    printstyled(io, string(entry_label, hline), "\n"; color=:light_black)
                    bb_guard_rail = bb_guard_rail_cont
                end
            end
            printstyled(io, "\e[$(start_column)G$(rail)\e[1G", color = :light_black)
            last_stack = stack
        else
            annotation = loc_annotations[idx]
            loc_method = loc_methods[idx]
            lineno = loc_lineno[idx]
            # Print location information right aligned. If the line below is too long, it'll overwrite this,
            # but that's what we want.
            if get(io, :color, false)
                method_start_column = cols - max_method_width - max_loc_width - 2
                filler = " "^(max_loc_width-length(annotation))
                printstyled(io, "\e[$(method_start_column)G$(annotation)$(filler)$(loc_method)\e[1G", color = :light_black)
            end
            printstyled(io, lineno, " "^(max_lineno_width - length(lineno) + 1); color = :light_black)
        end
        idx != last(bbrange) && print(io, bb_guard_rail)
        print_sep = false
        if idx == last(bbrange)
            print_sep = true
        end
        floop = true
        # print new nodes first in the right position
        while !isempty(new_nodes_perm) && new_nodes[Iterators.peek(new_nodes_perm)].pos == idx
            node_idx = popfirst!(new_nodes_perm)
            new_node = new_nodes[node_idx]
            node_idx += length(stmts)
            if !floop && !verbose_linetable
                print(io, " "^(max_lineno_width + 1))
            end
            if print_sep
                if idx == first(bbrange) && floop
                    print(io, bb_start_str)
                else
                    print(io, "│  ", " "^max_bb_idx_size)
                end
            end
            print_sep = true
            floop = false
            show_type = should_print_ssa_type(new_node.node)
            with_output_color(:yellow, io) do io′
                print_stmt(io′, node_idx, new_node.node, used, maxlength_idx, false, show_type)
            end
            if show_type
                expr_type_printer(io, new_node.typ, node_idx in used)
            end
            println(io)
        end
        if !floop && !verbose_linetable
            print(io, " "^(max_lineno_width + 1))
        end
        if print_sep
            if idx == first(bbrange) && floop
                print(io, bb_start_str)
            elseif idx == last(bbrange)
                print(io, "└", "─"^(1 + max_bb_idx_size), " ")
            else
                print(io, "│  ", " "^max_bb_idx_size)
            end
        end
        if idx == last(bbrange)
            bb_idx += 1
        end
        show_type = should_print_ssa_type(stmt)
        print_stmt(io, idx, stmt, used, maxlength_idx, true, show_type)
        if !isassigned(types, idx)
            # This is an error, but can happen if passes don't update their type information
            printstyled(io, "::#UNDEF", color=:red)
        elseif show_type
            typ = types[idx]
            expr_type_printer(io, typ, idx in used)
        end
        println(io)
    end
end

function show_ir(io::IO, code::CodeInfo, expr_type_printer=default_expr_type_printer; verbose_linetable=false)
    cols = displaysize(io)[2]
    used = BitSet()
    stmts = code.code
    types = code.ssavaluetypes
    cfg = compute_basic_blocks(stmts)
    max_bb_idx_size = length(string(length(cfg.blocks)))
    for stmt in stmts
        scan_ssa_use!(push!, used, stmt)
    end
    bb_idx = 1

    if isempty(used)
        maxlength_idx = 0
    else
        maxused = maximum(used)
        maxlength_idx = length(string(maxused))
    end
    if !verbose_linetable
        (loc_annotations, loc_methods, loc_lineno) = compute_ir_line_annotations(code)
        max_loc_width = maximum(length(str) for str in loc_annotations)
        max_lineno_width = maximum(length(str) for str in loc_lineno)
        max_method_width = maximum(length(str) for str in loc_methods)
    end
    max_depth = maximum(compute_inlining_depth(code.linetable, line) for line in code.codelocs)
    last_stack = []
    for idx in eachindex(stmts)
        if !isassigned(stmts, idx)
            # This is invalid, but do something useful rather
            # than erroring, to make debugging easier
            printstyled(io, :red, "#UNDEF\n")
            continue
        end
        stmt = stmts[idx]
        # Compute BB guard rail
        bbrange = cfg.blocks[bb_idx].stmts
        bbrange = bbrange.first:bbrange.last
        bb_idx_str = string(bb_idx)
        bb_pad = max_bb_idx_size - length(bb_idx_str)
        bb_type = length(cfg.blocks[bb_idx].preds) <= 1 ? "─" : "┄"
        bb_start_str = string(bb_idx_str, " ", bb_type, "─"^bb_pad, " ")
        bb_guard_rail_cont = string("│  ", " "^max_bb_idx_size)
        if idx == first(bbrange)
            bb_guard_rail = bb_start_str
        else
            bb_guard_rail = bb_guard_rail_cont
        end
        # Print linetable information
        if verbose_linetable
            stack = compute_loc_stack(code.linetable, code.codelocs[idx])
            # We need to print any stack frames that did not exist in the last stack
            ndepth = max(1, length(stack))
            rail = string(" "^(max_depth+1-ndepth), "│"^ndepth)
            start_column = cols - max_depth - 10
            for (i, x) in enumerate(stack)
                if i > length(last_stack) || last_stack[i] != x
                    entry = code.linetable[x]
                    printstyled(io, "\e[$(start_column)G$(rail)\e[1G", color = :light_black)
                    print(io, bb_guard_rail)
                    ssa_guard = " "^(maxlength_idx + 4 + (i - 1))
                    entry_label = "$(ssa_guard)$(entry.method) at $(entry.file):$(entry.line) "
                    hline = string("─"^(start_column-length(entry_label)-length(bb_guard_rail)+max_depth-i), "┐")
                    printstyled(io, string(entry_label, hline), "\n"; color=:light_black)
                    bb_guard_rail = bb_guard_rail_cont
                end
            end
            printstyled(io, "\e[$(start_column)G$(rail)\e[1G", color = :light_black)
            last_stack = stack
        else
            annotation = loc_annotations[idx]
            loc_method = loc_methods[idx]
            lineno = loc_lineno[idx]
            # Print location information right aligned. If the line below is too long, it'll overwrite this,
            # but that's what we want.
            if get(io, :color, false)
                method_start_column = cols - max_method_width - max_loc_width - 2
                filler = " "^(max_loc_width-length(annotation))
                printstyled(io, "\e[$(method_start_column)G$(annotation)$(filler)$(loc_method)\e[1G", color = :light_black)
            end
            printstyled(io, lineno, " "^(max_lineno_width-length(lineno)+1); color = :light_black)
        end
        idx != last(bbrange) && print(io, bb_guard_rail)
        if idx == last(bbrange) # print separator
            if idx == first(bbrange)
                print(io, bb_start_str)
            elseif idx == last(bbrange)
                print(io, "└", "─"^(1 + max_bb_idx_size), " ")
            else
                print(io, "│  ", " "^max_bb_idx_size)
            end
        end
        if idx == last(bbrange)
            bb_idx += 1
        end
        # convert statement index to labels, as expected by print_stmt
        if stmt isa Expr
            if stmt.head === :gotoifnot && length(stmt.args) == 2 && stmt.args[2] isa Int
                stmt = GotoIfNot(stmt.args[1], block_for_inst(cfg, stmt.args[2]::Int))
            elseif stmt.head === :enter && length(stmt.args) == 1 && stmt.args[1] isa Int
                stmt = Expr(:enter, block_for_inst(cfg, stmt.args[1]::Int))
            end
        elseif isa(stmt, GotoIfNot)
            stmt = GotoIfNot(stmt.cond, block_for_inst(cfg, stmt.dest))
        elseif stmt isa GotoNode
            stmt = GotoNode(block_for_inst(cfg, stmt.label))
        elseif stmt isa PhiNode
            e = stmt.edges
            stmt = PhiNode(Any[block_for_inst(cfg, e[i]) for i in 1:length(e)], stmt.values)
        end
        show_type = types isa Vector{Any} && should_print_ssa_type(stmt)
        print_stmt(io, idx, stmt, used, maxlength_idx, true, show_type)
        if types isa Vector{Any} # ignore types for pre-inference code
            if !isassigned(types, idx)
                # This is an error, but can happen if passes don't update their type information
                printstyled(io, "::#UNDEF", color=:red)
            elseif show_type
                typ = types[idx]
                expr_type_printer(io, typ, idx in used)
            end
        end
        println(io)
    end
end

@specialize

max_stars_repo_name: pochoi/julia

max_issues_repo_name: pochoi/julia

max_forks_repo_name: pochoi/julia


text:
"""
    serialize_model(model::MM, filename::String)::Serialized{MM} where {MM<:MetabolicModel}

Serialize the `model` to file `filename`, returning a [`Serialized`](@ref)
model that can be loaded back transparently by [`precache!`](@ref). The result
does _not_ contain the actual model data that are deferred to the disk; it may
thus be used to save memory, or send the model efficiently to remote workers
within distributed shared-storage environments.

The benefit of using this over "raw" `Serialization.serialize` is that the
resulting `Serialized` model will reload itself automatically with
[`precache!`](@ref) at first usage, which needs to be done manually when using
the `Serialization` package directly.
"""
function serialize_model(
    model::MM,
    filename::String,
)::Serialized{MM} where {MM<:MetabolicModel}
    open(f -> serialize(f, model), filename, "w")
    Serialized{MM}(filename)
end

"""
    serialize_model(model::Serialized, filename::String)::Serialized

Specialization of [`serialize_model`](@ref) that prevents nested serialization
of already-serialized models.
"""
function serialize_model(model::Serialized, filename::String)
    precache!(model)
    serialize_model(model.m, filename)
end


text:
using PolicyComplianceModel
using Test

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


text:
# hdf5utils.jl

using MIRTio: h5_get_keys, h5_get_attributes, h5_get_ismrmrd
using MIRTio: h5_get_ESC, h5_get_RSS, h5_get_kspace
using HDF5: h5open, h5read, attributes
using Test: @test

filename = tempname()

hdr = "ismrm header test"
T = Float32
esc = rand(T, 3,4,5)
rss = rand(T, 3,4,5)
ksp = complex.(rand(T, 3,4,5), rand(T, 3,4,5))
h5open(filename, "w") do file
    write(file, "ismrmrd_header", hdr)
    write(file, "reconstruction_esc", esc)
    write(file, "reconstruction_rss", rss)
    write(file, "kspace", ksp)
    attributes(file)["kspace"] = "spiral"
end

@test hdr == h5read(filename, "ismrmrd_header")
@test rss == h5read(filename, "reconstruction_rss")
@test ksp == h5read(filename, "kspace")

pdims = x -> permutedims(x, 3:-1:1)
@test h5_get_keys(filename) isa Array{String}
@test h5_get_keys(filename)[2] == "kspace"
@test h5_get_attributes(filename) isa Dict
@test h5_get_ismrmrd(filename) == hdr
@test h5_get_ESC(filename ; T=T) == pdims(esc)
@test h5_get_RSS(filename ; T=T) == pdims(rss)
@test h5_get_kspace(filename ; T=ComplexF32) == pdims(ksp)

#=
h5open(filename, "r") do fid
#   @test h5_getcomplextype(fid["kspace"]) == ComplexF32
    @test eltype(fid["kspace"]) == ComplexF32
end
=#


text:
function check_variables(z)
    q3, n, sϕ, b, sb = unpack(z)

    if any([ni <= 0.0 for ni in n])
        # println("n not in cone")
        return true
    end

    if any([si <= 0.0 for si in sϕ])
        # println("ϕ not in cone")
        return true
    end

    if any([!κ_so(bi)[2] for bi in b])
        # println("b not in cone")
        return true
    end

    if any([!κ_so(si)[2] for si in sb])
        # println("bs not in cone")
        return true
    end

    return false
end

function r(z, θ)
    q3, n, sϕ, b, sb = unpack(z)
    q2 = view(θ, nq .+ (1:nq))
    q1 = view(θ, 1:nq)
    u1 = view(θ, 2 * nq .+ (1:nu))
    h = θ[end-1]
    μ = θ[end]

    λ = vcat([[view(b[i], 2:3); n[i]] for i = 1:num_contacts]...) # contact forces
    ϕ = signed_distance(model, q3) # signed-distance function
    vT = P_func(model, q3) * (q3 - q2) ./ h

    e = [1.0; 0.0; 0.0]

    sb_stack = vcat([view(sb[i], 2:3) for i = 1:num_contacts]...)
    fc_stack = vcat([b[i][1] - model.μ * n[i] for i = 1:num_contacts]...)
    cp_stack = vcat([cone_product(sb[i], b[i]) - μ * e for i = 1:num_contacts]...)

    # action optimality conditions
    [dynamics(model, q1, q2, q3, u1, λ, h);
     sϕ - ϕ;
     n .* sϕ .- μ;

     # maximum dissipation optimality conditions
     vT - sb_stack;
     fc_stack;
     cp_stack]
end

# Jacobian
function Rz(z, θ)
    # differentiate r
    _r(w) = r(w, θ)
    _R = ForwardDiff.jacobian(_r, z)
    return _R
end

function Rθ(z, θ)
    # differentiate r
    _r(w) = r(z, w)
    _R = ForwardDiff.jacobian(_r, θ)

    return _R
end

"""
    step particle system
        solves 1-step feasibility problem
"""
function step(q1, q2, u1, h;
    r_tol = 1.0e-5,
    μ_tol = 1.0e-5,
    z_init = 1.0,
    μ_init = 1.0,
    max_iter = 100)
    # 1-step optimization problem:
    #     find z
    #     s.t. r(z) = 0
    #
    # z = (q, n, sϕ, b, sb)
    #     s are slack variables for convenience
    #     b[2:3] is the friction force

    z = initialize(q2, num_var, z_init = z_init)
    μ = μ_init # barrier parameter
    θ = [q1; q2; u1; h; μ] # problem data

    flag = false

    for k = 1:10
        for i = 1:max_iter
            # compute residual, residual Jacobian
            res = r(z, θ)
            if norm(res) < r_tol
                # println("   iter ($i) - norm: $(norm(res))")
                # println("     μ = $μ")
                # return z, true
                flag = true
                continue
            end

            jac = Rz(z, θ)

            # compute step
            Δ = jac \ res

            # line search the step direction
            α = 1.0

            iter = 0
            while check_variables(z - α * Δ) # backtrack inequalities
                α = 0.5 * α

                # println("   α = $α")
                iter += 1
                if iter > 50
                    @error "backtracking line search fail"
                    q3, n, sϕ, b, sb = unpack(z)
                    println("q3: $q3")
                    println("q2: $q2")
                    println("q1: $q1")
                    println("μ: $μ")
                    flag = false
                    return q2, zeros(num_contacts), zeros(3 * num_contacts), zeros(nq, nq), zeros(nq, nq), zeros(nq , nu), false
                end
            end

            while norm(r(z - α * Δ, θ), Inf)^2.0 >= (1.0 - 0.001 * α) * norm(res, Inf)^2.0
                α = 0.5 * α
                # println("   α = $α")

                iter += 1
                if iter > 50
                    @error "line search fail"
                    q3, n, sϕ, b, sb = unpack(z)
                    println("q3: $q3")
                    println("q2: $q2")
                    println("q1: $q1")
                    println("μ: $μ")

                    flag = false
                    return q2, zeros(num_contacts), zeros(3 * num_contacts), zeros(nq, nq), zeros(nq, nq), zeros(nq , nu), false
                end
            end

            # update
            z .-= α * Δ
        end

        if μ < μ_tol
            break
        else
            μ = 0.1 * μ
            θ[end] = μ
        end
    end

    Δz = -1.0 * Rz(z, θ) \ Rθ(z, θ)

    q3, n, sϕ, b, sb = unpack(z)
    Δq1 = view(Δz, 1:nq, 1:nq)
    Δq2 = view(Δz, 1:nq, nq .+ (1:nq))
    Δu1 = view(Δz, 1:nq, 2 * nq .+ (1:nu))

    return q3, n, b, Δq1, Δq2, Δu1, flag
end


text:
using Test, GpABC, DelimitedFiles
import GpABC.covariance, GpABC.get_hyperparameters_size

struct NoGradKernel <: AbstractGPKernel
    rbf_iso_kernel::GpABC.SquaredExponentialIsoKernel
end

function NoGradKernel()
    NoGradKernel(GpABC.SquaredExponentialIsoKernel())
end

function get_hyperparameters_size(ker::NoGradKernel, x::AbstractArray{Float64, 2})
    get_hyperparameters_size(ker.rbf_iso_kernel, x)
end

function covariance(ker::NoGradKernel, log_theta::AbstractArray{Float64, 1},
    x1::AbstractArray{Float64, 2}, x2::AbstractArray{Float64, 2})
    covariance(ker.rbf_iso_kernel, log_theta, x1, x2)
end

@testset "GP Reression Tests" begin
    training_x = readdlm("$(@__DIR__)/training_x1.csv")
    training_y = readdlm("$(@__DIR__)/training_y1.csv")
    test_x = readdlm("$(@__DIR__)/test_x1.csv")

    gpem = GPModel(training_x, training_y, test_x)

    theta_mle = gp_train(gpem)
    @test [58.7668, 2.94946, 5.84869] ≈ theta_mle rtol=1e-5
    @test -183.80873007668373 ≈ gp_loglikelihood(gpem)
    test_grad = gp_loglikelihood_grad(log.(gpem.gp_hyperparameters), gpem)
    @test all(abs.(test_grad - [46540815059614715e-22, 6451822871511581e-22, 20565826801544063e-22]) .< 1e-5)

    regression_mean = readdlm("$(@__DIR__)/test_y1_mean.csv")
    regression_var = readdlm("$(@__DIR__)/test_y1_var.csv")
    (r_mean, r_var) = gp_regression(gpem, batch_size=10)
    @test r_mean ≈ regression_mean
    @test r_var ≈ regression_var rtol=1e-3

    regression_cov = readdlm("$(@__DIR__)/test_y1_cov.csv", ',')
    (r_mean, r_cov) = gp_regression(gpem, batch_size=10, full_covariance_matrix=true)
    @test r_cov ≈ regression_cov rtol=1e-3

    gpem = GPModel(training_x, training_y, NoGradKernel())
    theta_mle = gp_train(gpem)
    @test [58.7668, 2.94946, 5.84869] ≈ theta_mle rtol=1e-3

    # Test that 1-d inputs work
    training_x = reshape(training_x, size(training_x, 1))
    training_y = reshape(training_y, size(training_y, 1))
    gpem = GPModel(training_x, training_y)
    theta_mle = gp_train(gpem)

end


text:
"""
    MessageOnlyLogger(stream=stderr, min_level=Info)

A simple logger for logging only the message field.
"""
struct MessageOnlyLogger <: AbstractLogger
    stream::IO
    min_level::LogLevel
end

MessageOnlyLogger(stream::IO=stderr, level=Logging.Info) = MessageOnlyLogger(stream, level)

Logging.shouldlog(logger::MessageOnlyLogger, level, _module, group, id) = true
Logging.min_enabled_level(logger::MessageOnlyLogger) = logger.min_level
Logging.catch_exceptions(logger::MessageOnlyLogger) = false

function Logging.handle_message(logger::MessageOnlyLogger,
                                level, message, _module, group, id, filepath, line;
                                kwargs...)
    println(logger.stream, message)
end


text:
include("LocalField/LocalField.jl")
include("LocalField/Elem.jl")
include("LocalField/pAdic.jl")
include("LocalField/qAdic.jl")
include("LocalField/Ring.jl")
include("LocalField/Conjugates.jl")
include("LocalField/Poly.jl")
include("LocalField/Completions.jl")


text:
module AffineM

export affine, coprime, affinebrute

using Main.UtilsM: removeExtras
using Main.ComparisonM: faComparison
using Main.FAM: fa

# Affine Substitute
function affine(text::String, alphabet::String, nMultiply::Int, nAdd::Int)

    return join([
        alphabet[(findfirst(isequal(letter), alphabet)*nMultiply+nAdd+1)%length(alphabet)+1]
        for letter in text
    ])

end # function

# Coprime Calculator
function coprime(n::Int)

    return [i for i in 2:n if gcd(i, n) == 1]

end # function

# Affine brute scoring mechanism
function affinebrute(text::String, alphabet::String)

    alphabet = uppercase(alphabet)
    text = removeExtras(uppercase(text), alphabet)

    bestScore = -10000
    bestDecrypt = ""
    for i in coprime(length(alphabet))
        for u = 1:length(alphabet)
            aff = affine(text, alphabet, i, u)
            if faComparison(fa(aff, alphabet)) > bestScore
                bestDecrypt = aff
                bestScore = faComparison(fa(aff, alphabet))
            end # if
        end # for
    end # for

    return bestDecrypt

end # function

end # module

max_stars_repo_name: Squalm/Cipher-Tools

max_issues_repo_name: Squalm/Cipher-Tools

max_forks_repo_name: Squalm/Cipher-Tools


text:
using SpinMC

a1 = (1.0, 0.0, 0.0)
a2 = (0.0, 1.0, 0.0)
a3 = (0.0, 0.0, 1.0)
uc = UnitCell(a1,a2,a3)

b = addBasisSite!(uc, (0.0, 0.0, 0.0))
M = [-1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 -1.0]
addInteraction!(uc, b, b, M, (1, 0, 0))
addInteraction!(uc, b, b, M, (0, 1, 0))
addInteraction!(uc, b, b, M, (0, 0, 1))

L = (8, 8, 8)
lattice = Lattice(uc, L)
thermalizationSweeps = 10000
measurementSweeps = 100000

tmin = 0.1
tmax = 10.0
N = 32
temperature = zeros(N)
heat = zeros(N)
dheat = zeros(N)
magnetization = zeros(N)
dmagnetization = zeros(N)
for i in 1:N
    println("Performing MC simulation ", i, "/", N)
    beta = 1.0 / (tmax * (tmin / tmax)^((i-1)/(N-1)))
    m = MonteCarlo(lattice, beta, thermalizationSweeps, measurementSweeps, reportInterval=20000)
    run!(m)
    temperature[i] = 1.0 / beta

    # Calculate specific heat
    # (Note that if file output is activated in the run!() call, the specific heat is automatically included.)
    c(e) = beta * beta * (e[2] - e[1] * e[1]) * length(m.lattice)
    ∇c(e) = [-2.0 * beta * beta * e[1] * length(m.lattice), beta * beta * length(m.lattice)]
    heat[i] = mean(m.observables.energy, c)
    dheat[i] = std_error(m.observables.energy, ∇c)
    # Calculate magnetization
    magnetization[i] = mean(m.observables.magnetization)
    dmagnetization[i] = std_error(m.observables.magnetization)
end

using Plots
display(plot(temperature, heat, yerror=dheat, xaxis=:log, yrange=(0.0,2.5), xrange=(0.1,10.0), xlabel="temperature", ylabel="specific heat", label="8*8*8 cubic lattice"))
display(plot(temperature, magnetization, yerror=dmagnetization, xaxis=:log, yrange=(0.0,1.0), xrange=(0.1,10.0), xlabel="temperature", ylabel="magnetization", label="8*8*8 cubic lattice"))


text:
##############################################################################
#
# Wishart Distribution
#
#  Parameters nu and S such that E(X) = nu * S
#  See the rwish and dwish implementation in R's MCMCPack
#  This parametrization differs from Bernardo & Smith p 435
#  in this way: (nu, S) = (2.0 * alpha, 0.5 * beta^-1)
#
##############################################################################

immutable Wishart <: ContinuousMatrixDistribution
    nu::Float64
    Schol::Cholesky{Float64}
    function Wishart(n::Real, Sc::Cholesky{Float64})
        if n > size(Sc, 1) - 1
            new(float64(n), Sc)
        else
            error("Wishart parameters must be df > p - 1")
        end
    end
end

Wishart(nu::Real, S::Matrix{Float64}) = Wishart(nu, cholfact(S))

function insupport(W::Wishart, X::Matrix{Float64})
    return size(X, 1) == size(X, 2) && isApproxSymmmetric(X) &&
           size(X, 1) == size(W.Schol, 1) && hasCholesky(X)
end
# This just checks if X could come from any Wishart
function insupport(::Type{Wishart}, X::Matrix{Float64})
    return size(X, 1) == size(X, 2) && isApproxSymmmetric(X) && hasCholesky(X)
end

mean(w::Wishart) = w.nu * w.Schol[:U]' * w.Schol[:U]

pdf(W::Wishart, X::Matrix{Float64}) = exp(logpdf(W, X))

size(W::Wishart) = size(W.Schol)
dim(W::Wishart) = size(W.Schol, 1)

function expected_logdet(W::Wishart)
    logd = 0.
    d = dim(W)

    for i=1:d
        logd += digamma(0.5 * (W.nu + 1 - i))
    end

    logd += d * log(2)
    logd += logdet(W.Schol)

    return logd
end

function lognorm(W::Wishart)
    d = dim(W)
    return (W.nu / 2) * logdet(W.Schol) + (d * W.nu / 2) * log(2) + lpgamma(d, W.nu / 2)
end

function logpdf(W::Wishart, X::Matrix{Float64})
    if !insupport(W, X)
        return -Inf
    else
        d = dim(W)
        logd = -lognorm(W)
        logd += 0.5 * (W.nu - d - 1.0) * logdet(X)
        logd -= 0.5 * trace(W.Schol \ X)
        return logd
    end
end

function rand(w::Wishart)
    p = size(w.Schol, 1)
    X = zeros(p, p)
    for ii in 1:p
        X[ii, ii] = sqrt(rand(Chisq(w.nu - ii + 1)))
    end
    if p > 1
        for col in 2:p
            for row in 1:(col - 1)
                X[row, col] = randn()
            end
        end
    end
    Z = X * w.Schol[:U]
    return Z' * Z
end

function entropy(W::Wishart)
    d = dim(W)
    return lognorm(W) - (W.nu - d - 1) / 2 * expected_logdet(W) + W.nu * d / 2
end

var(w::Wishart) = error("Not yet implemented")


text:
module DifferentialGeometry

export Δ, poisson

using SparseArrays, LinearAlgebra

"""the cotangent of the angle between two vectors"""
cotan( v1::T, v2::T ) where {T<:AbstractArray} = (v1⋅v2)/norm(v1×v2)

"""
    Δ(F, P)

Laplace-Beltrami operator built from list of triangle indices and points.

"""
function Δ( F::Vector{T}, P::Vector{PT} ) where {T<:Union{AbstractArray,NTuple{3}}, FT, PT<:AbstractVector{FT}}
  n = length(F)*9
  indexI = Vector{Int}(undef,n)
  indexJ = Vector{Int}(undef,n)
  element = Vector{FT}(undef,n)
  s = 1
  @inbounds for ijk in F
    i,j,k = ijk
    Pᵢ = P[i]
    Pⱼ = P[j]
    Pk = P[k]
    ij = Pⱼ-Pᵢ
    jk = Pk-Pⱼ
    ki = Pᵢ-Pk
    a = cotan(ij,-ki)*FT(0.5)
    b = cotan(jk,-ij)*FT(0.5)
    c = cotan(ki,-jk)*FT(0.5)
    #L[i,i] += b+c
    indexI[s] = i; indexJ[s] = i; element[s] = b+c; s+=1
    #L[i,j] -= c
    indexI[s] = i; indexJ[s] = j; element[s] = -c; s+=1
    #L[i,k] -= b
    indexI[s] = i; indexJ[s] = k; element[s] = -b; s+=1
    #L[j,i] -= c
    indexI[s] = j; indexJ[s] = i; element[s] = -c; s+=1
    #L[j,j] += c+a
    indexI[s] = j; indexJ[s] = j; element[s] = c+a; s+=1
    #L[j,k] -= a
    indexI[s] = j; indexJ[s] = k; element[s] = -a; s+=1
    #L[k,i] -= b
    indexI[s] = k; indexJ[s] = i; element[s] = -b; s+=1
    #L[k,j] -= a
    indexI[s] = k; indexJ[s] = j; element[s] = -a; s+=1
    #L[k,k] += a+b
    indexI[s] = k; indexJ[s] = k; element[s] = a+b; s+=1
  end
  sparse(indexI, indexJ, element)
end

Δ( F::Vector{IT}, P ) where {IT<:Integer} = Δ(Iterators.partition( F, 3 ) |> collect, P)

"""
    constrain( L, pin )

Add Dirichlet boundary conditions to Linear Operator L.
return tuple with reduced L and columns of L representing constrained dofs
"""
function constrain_system( L_in::S, pin ) where {S<:AbstractArray}
  L = copy(L_in)
  I_ = sparse(one(eltype(S))*I,size(L)...)
  L_constrained = L[:,pin]
  L[:,pin] = I_[:,pin]
  L[pin,:] = I_[pin,:]

  (L,L_constrained)
end

"""
  add boundary conditions to rhs of linear system
"""
function constrain_rhs(rhs, u, pin, L_constrained)
  b = rhs - L_constrained*u[pin,:]

  for ic in pin
    b[ic,:] = u[ic,:]
  end
  b
end

"""
    poisson(L, b, u_boundary, boundary_index)

Solve the Poisson problem Lu = b, where L is a Laplace-beltrami operator with Dirichlet boundary specified by u_boundary at locations boundary_index.
"""
function poisson(L, b, u_boundary, boundary_index)
  (Lₒ, Lpin) = constrain_system(L, boundary_index)

  rhs = constrain_rhs(b, u_boundary, boundary_index, Lpin)

  Lₒ\rhs
end

poisson(L, b::Vector{P}, u_boundary::Vector{P}, boundary_index) where P <: AbstractVector =
  map(P,
      eachrow(poisson(L,
                      reduce(vcat, adjoint.(b)),
                      reduce(vcat, adjoint.(u_boundary)),
                      boundary_index)))

end


text:
#### Path

# First just do it for the case where there we assume
# tainted gotoifnots do not go in a loop!
# TODO: write a thing to detect this! (overdub predicates only in tainted ifs)
# implement snapshotting function state as an optimization for branch exploration
mutable struct Path
    path::BitVector
    cursor::Int
end

Path() = Path([], 1)

function increment!(bitvec)
    for i=1:length(bitvec)
        if bitvec[i] === true
            bitvec[i] = false
        else
            bitvec[i] = true
            break
        end
    end
end

function reset!(p::Path)
    p.cursor=1
    increment!(p.path)
    nothing
end

function alldone(p::Path) # must be called at the end of the function!
    all(identity, p.path)
end

function current_predicate!(p::Path)
    #bt = backtrace()
    #Base.show_backtrace(stdout, bt)
    if p.cursor > length(p.path)
        push!(p.path, false)
    else
        p.path[p.cursor]
    end
    val = p.path[p.cursor]
    p.cursor+=1
    val
end

alldone(c) = alldone(c.metadata[2])
reset!(c) = reset!(c.metadata[2])
current_predicate!(c) = current_predicate!(c.metadata[2])

#=
julia> p=Path()
Path(Bool[], 1)

julia> alldone(p) # must be called at the end of a full run
true

julia> current_predicate!(p)
false

julia> alldone(p) # must be called at the end of a full run
false

julia> current_predicate!(p)
false

julia> p
Path(Bool[false, false], 3)

julia> alldone(p)
false

julia> reset!(p)

julia> p
Path(Bool[true, false], 1)

julia> current_predicate!(p)
true

julia> current_predicate!(p)
false

julia> alldone(p)
false

julia> reset!(p)

julia> p
Path(Bool[false, true], 1)

julia> current_predicate!(p)
false

julia> current_predicate!(p)
true

julia> reset!(p)

julia> current_predicate!(p)
true

julia> current_predicate!(p)
true

julia> alldone(p)
true
=#

"""
`abstract_run(g, ctx, overdubbed_fn, args...)`

First rewrites every if statement

```julia
if <expr>
  ...
end

as

```julia
cond = <expr>
if istainted(ctx, cond) ? @amb(true, false) : cond
  ...
end
```

Then runs `g(Cassette.overdub(ctx, overdubbed_fn, args...)`
as many times as there are available paths. i.e. `2^n` ways
where `n` is the number of tainted branch conditions.

# Examples:
```
meta = Any[]
abstract_run(ctx, f. args...) do result
    push!(meta, metadata(result, ctx))
end
# do something to merge metadata from all the runs
```
"""
function abstract_run(acc, ctx::Cassette.Context, overdub_fn, args...; verbose=true)
    path = Path()
    pass_ctx = Cassette.similarcontext(ctx, metadata=(ctx.metadata, path), pass=AbsintPass)

    while true
        acc(Cassette.recurse(pass_ctx, ()->overdub_fn(args...)))

        verbose && println("Explored path: ", path)
        alldone(path) && break
        reset!(path)
    end
end

"""
`istainted(ctx, cond)`

Does `cond` have any metadata?
"""
function istainted(ctx, cond)
    error("Method needed: istainted(::$(typeof(ctx)), ::Bool)." *
          " See docs for `istainted`.")
end

# Must return 7 exprs
function rewrite_branch(ctx, stmt, extraslot, i)
    # turn
    #   gotoifnot %p #g
    # into
    #   %t = istainted(%p)
    #   gotoifnot %t #orig
    #   %rec = @amb true false
    #   gotoifnot %rec #orig+1 (the next statement after gotoifnot)

    exprs = Any[]
    cond = stmt.args[1]        # already an SSAValue

    # insert a check to see if SSAValue(i) isa Tainted
    istainted_ssa = Core.SSAValue(i)
    push!(exprs, :($(Expr(:nooverdub, istainted))($(Expr(:contextslot)),
                              $cond)))

    # not tainted? jump to the penultimate statement
    push!(exprs, Expr(:gotoifnot, istainted_ssa, i+5))

    # tainted? then use current_predicate!(SSAValue(1))
    current_pred = i+2
    push!(exprs, :($(Expr(:nooverdub, current_predicate!))($(Expr(:contextslot)))))

    # Store the interpreter-provided predicate in the slot
    push!(exprs, Expr(:(=), extraslot, SSAValue(i+2)))

    push!(exprs, Core.GotoNode(i+6))

    push!(exprs, Expr(:(=), extraslot, cond))

    # here we put in the original code
    stmt1 = copy(stmt)
    stmt.args[1] = extraslot
    push!(exprs, stmt)

    exprs
end

function rewrite_ir(ctx, ref)
    # turn
    #   <val> ? t : f
    # into
    #   istainted(<val>) ? current_predicate!(p) : <val> ? t : f

    ir = ref.code_info
    ir = copy(ir)

    extraslot = gensym("tmp")
    push!(ir.slotnames, extraslot)
    push!(ir.slotflags, 0x00)
    extraslot = Core.SlotNumber(length(ir.slotnames))

    Cassette.insert_statements!(ir.code, ir.codelocs,
        (stmt, i) -> Base.Meta.isexpr(stmt, :gotoifnot) ? 7 : nothing,
        (stmt, i) -> rewrite_branch(ctx, stmt, extraslot, i))

    ir.ssavaluetypes = length(ir.code)
    # Core.Compiler.validate_code(ir)
    #@show ref.method
    #@show ir
    return ir
end

const AbsintPass = Cassette.@pass rewrite_ir


text:
function D3Trees.D3Tree(D::DESPOT; title="DESPOT Tree", kwargs...)
    lenb = length(D.children)
    lenba = length(D.ba_children)
    len = lenb + lenba
    K = length(D.scenarios[1])
    children = Vector{Vector{Int}}(undef, len)
    text = Vector{String}(undef, len)
    tt = fill("", len)
    link_style = fill("", len)
    L = calc_L(D)
    for b in 1:lenb
        children[b] = D.children[b] .+ lenb
        text[b] = @sprintf("""
                           o:%s (|Φ|:%3d)
                           L:%6.2f, U:%6.2f
                           l:%6.2f, μ:%6.2f, l₀:%6.2f""",
                           b==1 ? "<root>" : string(D.obs[b]),
                           length(D.scenarios[b]),
                           L[b],
                           D.U[b],
                           D.l[b],
                           D.mu[b],
                           D.l_0[b]
                          )
        tt[b] = """
                o: $(b==1 ? "<root>" : string(D.obs[b]))
                |Φ|: $(length(D.scenarios[b]))
                L: $(L[b])
                U: $(D.U[b])
                l: $(D.l[b])
                μ: $(D.mu[b])
                l₀: $(D.l_0[b])
                $(length(D.children[b])) children
                """
        link_width = 20.0*sqrt(length(D.scenarios[b])/K)
        link_style[b] = "stroke-width:$link_width"
        for ba in D.children[b]
            link_style[ba+lenb] = "stroke-width:$link_width"
        end

        for ba in D.children[b]
            weighted_sum_U = 0.0
            for bp in D.ba_children[ba]
                weighted_sum_U += length(D.scenarios[bp]) * D.U[bp]
            end
            U = (D.ba_Rsum[ba] + infer_discount(D) * weighted_sum_U) / length(D.scenarios[b])
            children[ba+lenb] = D.ba_children[ba]
            text[ba+lenb] = @sprintf("""
                                     a:%s (ρ:%6.2f)
                                     L:%6.2f, U:%6.2f,
                                     l:%6.2f, μ:%6.2f""",
                                     D.ba_action[ba], D.ba_rho[ba],
                                     L[ba+lenb], U,
                                     ba_l(D, ba), D.ba_mu[ba])
            tt[ba+lenb] = """
                          a: $(D.ba_action[ba])
                          ρ: $(D.ba_rho[ba])
                          L: $(L[ba+lenb])
                          U: $U
                          l: $(ba_l(D, ba))
                          μ: $(D.ba_mu[ba])
                          $(length(D.ba_children[ba])) children
                          """
        end
    end
    return D3Tree(children;
                  text=text,
                  tooltip=tt,
                  link_style=link_style,
                  title=title,
                  kwargs...
                 )
end

Base.show(io::IO, mime::MIME"text/html", D::DESPOT) = show(io, mime, D3Tree(D))
Base.show(io::IO, mime::MIME"text/plain", D::DESPOT) = show(io, mime, D3Tree(D))

"""
Return a vector of lower bounds L of length lenb+lenba, with b nodes first followed by ba nodes.
"""
function calc_L(D::DESPOT)
    lenb = length(D.children)
    lenba = length(D.ba_children)
    if lenb == 1
        @assert lenba == 0
        return [D.l_0[1]]
    end
    len = lenb + lenba
    cache = fill(NaN, len)
    disc = infer_discount(D)
    fill_L!(cache, D, 1, disc)
    return cache
end

function infer_discount(D::DESPOT)
    # @assert !isempty(D.children[1])
    # K = length(D.scenarios[0])
    # firstba = first(D.children[1])
    # lambda = D.ba_rsum[firstba]/K - D.ba_rho[firstba]
    disc = D._discount
    return disc
end

"""
Fill all the elements of the cache for b and children of b and return L[b]
"""
function fill_L!(cache::Vector{Float64}, D::DESPOT, b::Int, disc::Float64)
    K = length(D.scenarios[1])
    lenb = length(D.children)
    if isempty(D.children[b])
        L = D.l_0[b]*K/(length(D.scenarios[b])*disc^D.Delta[b])
        cache[b] = L
        return L
    else
        max_L = -Inf
        for ba in D.children[b]
            weighted_sum_L = 0.0
            for bp in D.ba_children[ba]
                weighted_sum_L += length(D.scenarios[bp]) * fill_L!(cache, D, bp, disc)
            end
            new_L = (D.ba_Rsum[ba] + disc * weighted_sum_L) / length(D.scenarios[b])
            cache[lenb+ba] = new_L
            max_L = max(max_L, new_L)
        end
        cache[b] = max_L
        return max_L
    end
end


text:
dir = joinpath(@__DIR__, "../config")
isdir(dir) || mkdir(dir)

if !isfile(joinpath(dir, "odbc.ini"))
    open(joinpath(dir, "odbc.ini"), "w") do io
        write(io, "[ODBC Data Sources]\n\n[ODBC]\nTrace=0\nTraceFile=stderr\n")
    end
end

if !isfile(joinpath(dir, "odbcinst.ini"))
    open(joinpath(dir, "odbcinst.ini"), "w") do io
        write(io, "[ODBC Drivers]\n")
    end
end


text:
macro import_lrdomains_huge()
    :(
        using ComplexVisual:
            cv_destroy,
            CV_DomainCodomainParts,
            cv_get_trafo, cv_get_can_domain, cv_get_can_codomain,
            cv_get_cc_can_domain, cv_get_cc_can_codomain,
            cv_add,
            CV_DomainPosLayout, CV_CodomainPosLayout,
            cv_get_can_domain_l, cv_get_can_codomain_l,
            cv_do_lr_layout,
            CV_DomainCodomainLayout,
            cv_get_can_layout, cv_get_cc_can_layout,
            cv_get_update_math_domains, cv_get_actionpixel_update,
            cv_get_statepixel_update,
            cv_setup_domain_codomain_scene,
            cv_setup_lr_axis, cv_setup_lr_painters, cv_setup_lr_border,
            cv_scene_lr_start, cv_scene_lr_std,
            CV_LRSetupChain, cv_combine
    )
end

import Base:show

"""
SetupChain with layout (type inferable) and a vector for every
callback-type where all the callback functions are stored (type inference
is lost). See `CV_SceneSetupChain` for the reasons of this tradeoff.
"""
struct CV_LRSetupChain{layoutT} <: CV_SceneSetupChain # {{{
    layout               :: layoutT
    draw_once_func       :: Vector{Any}
    actionpixel_update   :: Vector{Any}
    statepixel_update    :: Vector{Any}
    update_painter_func  :: Vector{Any}
    update_state_func    :: Vector{Any}
    redraw_func          :: Vector{Any}
end

function CV_LRSetupChain(layout)
    return CV_LRSetupChain(layout, Vector(), Vector(), Vector(),
        Vector(), Vector(), Vector())
end

function cv_combine(old::CV_LRSetupChain;
        layout=missing, draw_once_func=missing, update_painter_func=missing,
        update_state_func=missing, actionpixel_update=missing,
        statepixel_update=missing, redraw_func=missing)
    new = ismissing(layout) ? old : CV_LRSetupChain(layout,
        old.draw_once_func, old.actionpixel_update, old.statepixel_update,
        old.update_painter_func, old.update_state_func, old.redraw_func)
    !ismissing(draw_once_func) && push!(new.draw_once_func, draw_once_func)
    !ismissing(actionpixel_update) && push!(new.actionpixel_update, actionpixel_update)
    !ismissing(statepixel_update) && push!(new.statepixel_update, statepixel_update)
    !ismissing(update_painter_func) && push!(new.update_painter_func, update_painter_func)
    !ismissing(update_state_func) && push!(new.update_state_func, update_state_func)
    !ismissing(redraw_func) && push!(new.redraw_func, redraw_func)
    return new
end

function cv_setup_cycle_state(setup::CV_LRSetupChain)
    state_counter = cv_get_state_counter(setup.layout)

    update_state_func = z -> begin
        state_counter()
        return nothing
    end
    return cv_combine(setup; update_state_func)
end
# }}}

"""
Saves domain and codomain (and their contexts) together with a trafo in
(the) layout.
"""
struct CV_DomainCodomainParts{parentT<:CV_Abstract2DLayout,
                    trafoT, domainT<:CV_Canvas, codomainT<:CV_Canvas,
                    ccdT, cccT} <: CV_2DLayoutWrapper  # {{{
    parent_layout   :: parentT
    trafo           :: trafoT
    can_domain      :: domainT
    can_codomain    :: codomainT
    cc_can_domain   :: ccdT
    cc_can_codomain :: cccT
end

@layout_composition_getter(trafo,           CV_DomainCodomainParts)
@layout_composition_getter(can_domain,      CV_DomainCodomainParts)
@layout_composition_getter(can_codomain,    CV_DomainCodomainParts)
@layout_composition_getter(cc_can_domain,   CV_DomainCodomainParts)
@layout_composition_getter(cc_can_codomain, CV_DomainCodomainParts)

function cv_destroy(l::CV_DomainCodomainParts)
    cv_destroy(l.cc_can_domain)
    cv_destroy(l.cc_can_codomain)
    cv_destroy(l.can_domain)
    cv_destroy(l.can_codomain)
    cv_destroy(l.parent_layout)
    return nothing
end

"""
convenience function to combine domain and codomain (and their contexts)
together with a trafo to (the) layout.
"""
function cv_add(layout::CV_Abstract2DLayout, trafo,
        can_domain::domainT=CV_Math2DCanvas(-2.0+2.0im, 2.0-2.0im, 50),
        can_codomain::codomainT=CV_Math2DCanvas(-2.0+2.0im, 2.0-2.0im, 50),
        cc_can_domain::CV_CanvasContext=cv_create_context(can_domain),
        cc_can_codomain::CV_CanvasContext=cv_create_context(can_codomain)
        ) where {domainT<:CV_Canvas, codomainT<:CV_Canvas}

    return CV_DomainCodomainParts(
        layout, trafo, can_domain, can_codomain,
        cc_can_domain, cc_can_codomain)
end

# }}}

"""
A `CV_DomainPosLayout` with domain positioned.
"""
struct CV_DomainPosLayout{parentT<:CV_Abstract2DLayout, canT, dcbT, styleT
                          } <: CV_2DLayoutWrapper    # {{{
    parent_layout   :: parentT
    can_domain_l    :: CV_2DLayoutPosition{canT, dcbT, styleT}
end

@layout_composition_getter(can_domain_l, CV_DomainPosLayout)
# }}}

"""
A `CV_CodomainPosLayout` with codomain positioned.
"""
struct CV_CodomainPosLayout{parentT<:CV_Abstract2DLayout,canT, dcbT, styleT
                            } <: CV_2DLayoutWrapper   # {{{
    parent_layout   :: parentT
    can_codomain_l  :: CV_2DLayoutPosition{canT, dcbT, styleT}
end

@layout_composition_getter(can_codomain_l, CV_CodomainPosLayout)

# }}}

"""
Do a left-right layout of domain and codomain canvas with a gap between.
"""
function cv_do_lr_layout(layout::CV_Abstract2DLayout, gap::Integer=50) # {{{
    can_domain = cv_get_can_domain(layout)
    can_codomain = cv_get_can_codomain(layout)
    can_domain_l = cv_add_canvas!(layout, can_domain,
        cv_anchor(can_domain, :northwest), (0, 0))
    can_codomain_l = cv_add_canvas!(layout, can_codomain,
        cv_anchor(can_codomain, :west),
        cv_translate(cv_anchor(can_domain_l, :east), gap, 0))
    return CV_CodomainPosLayout(
        CV_DomainPosLayout(layout, can_domain_l), can_codomain_l)
end # }}}

"""
Scene with `update_math_domains` for domain and codomain.
"""
struct CV_DomainCodomainLayout{parentT, umdT} <: CV_2DLayoutWrapper # {{{
    parent_layout            :: parentT
    update_math_domains      :: umdT
end

@layout_composition_getter(update_math_domains,     CV_DomainCodomainLayout)

"""
create scene for given scene setup.
"""
function cv_setup_domain_codomain_scene(setup::CV_SceneSetupChain)
    layout = setup.layout

    can_domain_l = cv_get_can_domain_l(layout)
    can_codomain_l = cv_get_can_codomain_l(layout)

    redraw_func = (future_layout) -> begin
        cc_can_layout = cv_get_cc_can_layout(future_layout)
        can_domain_l(cc_can_layout)
        can_codomain_l(cc_can_layout)
        return nothing
    end

    update_math_domains = (z, future_layout) -> begin
        cc_can_layout = cv_get_cc_can_layout(future_layout)
        for func in setup.update_painter_func
            func(z)
        end
        can_domain_l(cc_can_layout)
        can_codomain_l(cc_can_layout)
        return nothing
    end

    new_layout = CV_DomainCodomainLayout(layout, update_math_domains)

    actionpixel_update = (px, py, future_layout) -> begin
        can_layout = cv_get_can_layout(future_layout)
        canvas = can_domain_l.canvas
        lx, ly = cv_global2local(can_layout, can_domain_l, px, py)
        if (lx < 0) || (ly < 0) ||
              (lx > canvas.pixel_width) || (ly > canvas.pixel_height)
            return nothing
        end
        x, y = cv_pixel2math(canvas, lx, ly)
        update_math_domains(x + y*1im, future_layout)
        return nothing
    end

    statepixel_update = (px, py, future_layout) -> begin
        can_layout = cv_get_can_layout(future_layout)
        canvas = can_domain_l.canvas
        lx, ly = cv_global2local(can_layout, can_domain_l, px, py)
        if (lx < 0) || (ly < 0) ||
              (lx > canvas.pixel_width) || (ly > canvas.pixel_height)
            return nothing
        end
        x, y = cv_pixel2math(canvas, lx, ly)
        for func in setup.update_state_func
            func(x + y*1im)
        end
        return nothing
    end

    setup = cv_combine(setup; layout=new_layout,
        actionpixel_update, statepixel_update, redraw_func)

    return cv_setup_2dminimal_scene(setup)
end
# }}}

function cv_setup_lr_axis(setup::CV_SceneSetupChain,
        domain_re_ruler::CV_Ruler, domain_im_ruler::CV_Ruler,
        codomain_re_ruler::CV_Ruler, codomain_im_ruler::CV_Ruler) # {{{
    return cv_setup_lr_axis(setup,
        (domain_re_ruler, ), (domain_im_ruler, ),
        (codomain_re_ruler, ), (codomain_im_ruler, ))
end # }}}

function cv_setup_lr_axis(setup::CV_SceneSetupChain,
        domain_re_rulers::NTuple{A, CV_Ruler},
        domain_im_rulers::NTuple{B, CV_Ruler},
        codomain_re_rulers::NTuple{C, CV_Ruler},
        codomain_im_rulers::NTuple{D, CV_Ruler}) where {A, B, C, D} # {{{

    layout = setup.layout
    can_domain_l = cv_get_can_domain_l(layout)
    can_codomain_l = cv_get_can_codomain_l(layout)

    domain_re_axis = cv_ticks_labels(layout, can_domain_l,
        cv_south, domain_re_rulers)
    domain_im_axis = cv_ticks_labels(layout, can_domain_l,
        cv_west, domain_im_rulers)
    codomain_re_axis = cv_ticks_labels(layout, can_codomain_l,
        cv_south, codomain_re_rulers)
    codomain_im_axis = cv_ticks_labels(layout, can_codomain_l,
        cv_west, codomain_im_rulers)

    draw_once_func = future_layout -> begin
        cc_can_layout = cv_get_cc_can_layout(future_layout)
        domain_re_axis(cc_can_layout)
        domain_im_axis(cc_can_layout)
        codomain_re_axis(cc_can_layout)
        codomain_im_axis(cc_can_layout)
        return nothing
    end
    return cv_combine(setup; draw_once_func)
end  # }}}

function cv_setup_lr_axis(setup::CV_SceneSetupChain;
        label_style::CV_ContextStyle = cv_black →
                cv_fontface("serif") → cv_fontsize(20)) where {A, B, C, D} # {{{

    layout = setup.layout
    get_range_mean = (a, b) -> (
        range(ceil(Int64, a); stop=floor(Int64, b)), a+0.5*(b-a))
    function format_ticks(a, b)
        ra = range(ceil(Int64, a); stop=floor(Int64, b))
        mean = a+0.5*(b-a)
        return isempty(ra) ? cv_format_ticks("%.2f", mean) :
                             cv_format_ticks("%.0f", ra...)
    end

    can_domain = cv_get_can_domain(layout)
    can_codomain = cv_get_can_codomain(layout)

    domain_re_ticks = format_ticks(
        real(can_domain.corner_ul), real(can_domain.corner_lr))
    domain_im_ticks = format_ticks(
        imag(can_domain.corner_lr), imag(can_domain.corner_ul))
    codomain_re_ticks = format_ticks(
        real(can_codomain.corner_ul), real(can_codomain.corner_lr))
    codomain_im_ticks = format_ticks(
        imag(can_codomain.corner_lr), imag(can_codomain.corner_ul))
    app = CV_TickLabelAppearance(; label_style)

    return cv_setup_lr_axis(setup,
        CV_Ruler(domain_re_ticks, app), CV_Ruler(domain_im_ticks, app),
        CV_Ruler(codomain_re_ticks, app), CV_Ruler(codomain_im_ticks, app))
end  # }}}

function cv_setup_lr_painters(setup::CV_SceneSetupChain,
        codomain_painter::Tuple{codos1P, codos2P},
        domain_painter::Tuple{dos1P, dos2P}, save_action_pos) where {
            codos1P<:Union{CV_Painter, Nothing},
            codos2P<:Union{CV_Painter, Nothing},
            dos1P<:Union{CV_Painter, Nothing},
            dos2P<:Union{CV_Painter, Nothing}} # {{{

    layout = setup.layout
    state_counter = cv_get_state_counter(layout)

    cc_can_domain = cv_get_cc_can_domain(layout)
    cc_can_codomain = cv_get_cc_can_codomain(layout)

    update_painter_func = z -> begin
        save_action_pos.value = z

        state = state_counter.value
        if state == 1
            if domain_painter[1] !== nothing
                cv_paint(cc_can_domain, domain_painter[1])
            end
            if codomain_painter[1] !== nothing
                cv_paint(cc_can_codomain, codomain_painter[1])
            end
        elseif state == 2
            if domain_painter[2] !== nothing
                cv_paint(cc_can_domain, domain_painter[2])
            end
            if codomain_painter[2] !== nothing
                cv_paint(cc_can_codomain, codomain_painter[2])
            end
        end

        return nothing
    end

    redraw_func = layout -> begin
        if domain_painter[1] !== nothing
            cv_clear_cache(domain_painter[1])
        end
        if domain_painter[2] !== nothing
            cv_clear_cache(domain_painter[2])
        end
        if codomain_painter[1] !== nothing
            cv_clear_cache(codomain_painter[1])
        end
        if codomain_painter[2] !== nothing
            cv_clear_cache(codomain_painter[2])
        end
        update_painter_func(save_action_pos.value)
        return nothing
    end

    return cv_combine(setup; update_painter_func, redraw_func)
end # }}}

const cv_setup_lr_painters_default_phs = cv_op_source →
                cv_antialias(Cairo.ANTIALIAS_BEST) →
                cv_linewidth(3) → cv_black
const cv_setup_lr_painters_default_pvs = cv_op_source →
                cv_antialias(Cairo.ANTIALIAS_BEST) →
                cv_linewidth(3) → cv_white
const cv_setup_lr_painters_default_imgps = cv_op_over →
                     cv_antialias(Cairo.ANTIALIAS_NONE)

"""
convenience function for `cv_setup_lr_painters` with all
arguments as keyword arguments.
"""
function cv_setup_lr_painters(setup::CV_SceneSetupChain;
        cut_test=nothing,
        canvas_test_img=cv_example_image_letter(),
        img_painter_style=cv_setup_lr_painters_default_imgps,
        img_painter_domain=missing,
        img_painter_codomain=missing,
        portrait_painter_domain=missing,
        portrait_painter_codomain=missing,
        parallel_hlines_style=cv_setup_lr_painters_default_phs,
        parallel_vlines_style=cv_setup_lr_painters_default_pvs,
        parallel_hlines=cv_parallel_lines(1.0+0.0im),
        parallel_vlines=cv_parallel_lines(0.0+1.0im),
        parallel_lines_painter_domain=missing,
        parallel_lines_painter_codomain=missing,
        action_pos::CV_TranslateByOffset{ComplexF64}=CV_TranslateByOffset(ComplexF64)
        ) # {{{

    trafo = cv_get_trafo(setup.layout)
    translate = action_pos
    trafo_translate = w -> trafo(translate(w))

    if ismissing(img_painter_domain)
        img_painter_domain =
            img_painter_style ↦ CV_Math2DCanvasPainter(translate, canvas_test_img)
    end
    if ismissing(img_painter_codomain)
        img_painter_codomain =
            img_painter_style ↦ CV_Math2DCanvasPainter(trafo_translate,
                canvas_test_img,
                cut_test === nothing ? nothing : translate, cut_test)
    end

    if ismissing(portrait_painter_domain)
        portrait_painter_domain = CV_PortraitPainter(trafo)
    end
    if ismissing(portrait_painter_codomain)
        portrait_painter_codomain = CV_PortraitPainter()
    end

    if ismissing(parallel_lines_painter_domain)
        parallel_lines_painter_domain = (
            (parallel_hlines_style ↦
                CV_LinePainter(translate, parallel_hlines)) →
            (parallel_vlines_style ↦
                CV_LinePainter(translate, parallel_vlines)))
    end
    if ismissing(parallel_lines_painter_codomain)
        parallel_lines_painter_codomain = (
            (parallel_hlines_style ↦
                CV_LinePainter(trafo_translate, parallel_hlines, false,
                    cut_test === nothing ? nothing : translate, cut_test)) →
            (parallel_vlines_style ↦
                CV_LinePainter(trafo_translate, parallel_vlines, false,
                    cut_test === nothing ? nothing : translate, cut_test)))
    end

    domain_state1_painter = portrait_painter_domain
    domain_state2_painter = portrait_painter_domain
    if img_painter_domain !== nothing
        domain_state1_painter = domain_state1_painter → img_painter_domain
    end
    if parallel_lines_painter_domain !== nothing
        domain_state2_painter = domain_state2_painter → parallel_lines_painter_domain
    end

    codomain_state1_painter = portrait_painter_codomain
    codomain_state2_painter = portrait_painter_codomain
    if img_painter_codomain !== nothing
        codomain_state1_painter = codomain_state1_painter → img_painter_codomain
    end
    if parallel_lines_painter_codomain !== nothing
        codomain_state2_painter = codomain_state2_painter → parallel_lines_painter_codomain
    end

    return cv_setup_lr_painters(setup,
        (codomain_state1_painter, codomain_state2_painter),
        (domain_state1_painter, domain_state2_painter),
        action_pos)

end # }}}

"""
creates borders for the domain and codomain canvas.
"""
function cv_setup_lr_border(setup::CV_SceneSetupChain; width::Integer=2,
        style=cv_black)
    layout = setup.layout
    domain_border = cv_border(layout, cv_get_can_domain_l(layout), width;
        style)
    codomain_border = cv_border(layout, cv_get_can_codomain_l(layout), width;
        style)

    draw_once_func = future_layout -> begin
        cc_can_layout = cv_get_cc_can_layout(future_layout)
        domain_border(cc_can_layout)
        codomain_border(cc_can_layout)
        return nothing
    end

    return cv_combine(setup; draw_once_func)
end

function cv_scene_lr_start(scene::CV_SceneSetupChain;
        z_start::Union{ComplexF64, Missing}=missing,
        state_start::Union{Int, Missing}=missing) # {{{

    layout = scene.layout
    z = 0.0+0.0im
    if ismissing(z_start)
        can_codomain = cv_get_can_codomain(layout)
        z = (can_codomain.corner_ul + can_codomain.corner_lr)/2
    else
        z = z_start :: ComplexF64
    end
    if !ismissing(state_start)
        cv_set_value!(cv_get_state_counter(layout), state_start)
    end
    cv_get_update_math_domains(layout)(z, layout)
    return nothing
end # }}}

"""
creates "standard" scene with domain and codomain in a left-right layout.
"""
function cv_scene_lr_std(trafo,
        domain, codomain; cut_test=nothing, gap=80,
        axis_label_style=cv_black →
                  cv_fontface("sans-serif", Cairo.FONT_WEIGHT_BOLD) →
                  cv_fontsize(20), padding=30,
        lr_painters_kwargs = Dict(:cut_test => cut_test),
        lr_start_kwargs = Dict()) # {{{
    layout = CV_StateLayout(CV_2DLayout(), CV_CyclicValue(2))
    layout = cv_do_lr_layout(cv_add(layout, trafo, domain, codomain), gap)

    setup = cv_setup_cycle_state(CV_LRSetupChain(layout))
    setup = cv_setup_lr_painters(setup; lr_painters_kwargs...)
    setup = cv_setup_lr_axis(setup; label_style=axis_label_style)
    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; lr_start_kwargs...)
    return setup.layout
end # }}}

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


text:
using Mimi, Distributions

include("MyModel.jl")
using .MyModel

m = construct_MyModel()

# Step 2: Define the Simulation

"""
@defsim macro is the first step in the process and returns a SimulationDef

allows users to define random variables as distributions, and associate model
parameters with the defined random variables
"""

sd = @defsim begin

    # Define random variables. The rv() is only required
    # when defining correlations or sharing an RV across parameters
    # otherwise, you can use the shortcut syntax
    # to assign a distribution to a parameter name
    rv(name1) = Normal(1, 0.2)
    rv(name2) = Uniform(0.75, 1.25)
    rv(name3) = LogNormal(20, 4)

    # Define the sampling strategy
    sampling(LHSData,corrlist = [(:name1, :name2, 0.7), (:name1, :name3, 0.5)])

    # assign RVs to model Parameters
    grosseconomy.share = Uniform(0.2, 0.8)

    # you can use the *= operator to replace the values in the parameter with the
    # product of the original value and the value of the RV for the current
    # trial (note that in both lines below, all indexed values will be mulitplied by the
    # same draw from the given random parameter (name2 or Uniform(0.8, 1.2))

    emissions.sigma[:,Region1] *= name2
    emissions.sigma[2020:5:2050, (Region2, Region3)] *= Uniform(0.8, 1.2)

    # For parameters that have a region dimension, you can assign an array of distributions,
    # keyed by region label, which must match the region labels in the model
    grosseconomy.depk = [Region1 => Uniform(0.7, .9),
                        Region2 => Uniform(0.8, 1.),
                        Region3 => Truncated(Normal(), 0, 1)]

    # Indicate which variables to save for each model run
    # The syntax is: component_name.variable_name
    save(grosseconomy.K, grosseconomy.YGROSS,
        emissions.E, emissions.E_Global)
end

# STEP 3: Run Simulation

# Run 100 trials and optionally save results to the indicated directories
si = run(sd, m, 100; trials_output_filename = "/tmp/trialdata.csv", results_output_dir = "/tmp/tutorial4")

# Explore the results saved in-memory by using getdataframe with the returned SimulationInstance.
# Values are saved from each trial for each variable or parameter specified by the call to "save()" at the end of the @defsim block.
K_results = getdataframe(si, :grosseconomy, :K)
E_results = getdataframe(si, :emissions, :E)

explore(si)

max_stars_repo_name: tlorans/SustainableFinance

max_issues_repo_name: tlorans/SustainableFinance

max_forks_repo_name: tlorans/SustainableFinance


text:
function mytest1()
x_data_a2 = [0.1, 0.102, 0.153, 0.2, 0.202,
    0.248, 0.297, 0.3, 0.348, 0.399, 0.4,
    0.5, 0.5, 0.6, 0.601, 0.699, 0.7, 0.794,
    0.8, 0.9]

a2_data_a2 = [0.9, 0.866488436, 0.76805011,
    0.63, 0.649506857, 0.524300393, 0.38397441,
    0.347, 0.256229362, 0.162509557, 0.149,
    0.073868064, 0.0665, 0.0281, 0.030561125,
    0.014647489, 0.0116, 0.006033884,
    0.0042, 0.0012]

T_data_a2 = [998., 1073., 1073., 998., 1073.,
    1073., 1073., 998., 1073., 1073., 998.,
    1073., 998., 998., 1073., 1073., 998.,
    1073., 998., 998.]

x_data_Tf = Float64[0, 0, 0, 0.049533151, 0.097097104, 0.1,
    0.149600304, 0.192848845, 0.2, 0.200826782, 0.25,
    0.250912166, 0.280680928, 0.295130904, 0.3, 0.304846515,
    0.306117172, 0.31, 0.312350413, 0.31546046, 0.32,
    0.324720543, 0.33, 0.331772771, 0.333, 0.34, 0.346155446,
    0.35, 0.350583994, 0.354739648, 0.369980613, 0.37,
    0.375348261, 0.391228458, 0.395210589, 0.4, 0.418663535,
    0.42, 0.432446439, 0.446446984, 0.45, 0.451241823,
    0.465687834, 0.47, 0.480656179, 0.485668569, 0.49637962,
    0.5, 0.5, 0.501330506, 0.516401957, 0.52, 0.526433927,
    0.530199064, 0.535751614, 0.546510668, 0.55, 0.55092091,
    0.555777151, 0.57, 0.577064846, 0.6, 0.611741422, 0.631046965,
    0.648471584, 0.65, 0.68, 0.7, 0.751358143, 0.8,
    0.848845459, 0.9, 0.947644142, 1, 1, 1]

Tf_data_Tf = Float64[1041, 1043.45, 1044.626691, 1023.27028,
    997.9162418, 983, 962.5594888, 913.7211707, 908,
    891.9362264, 833.5, 803.3174683, 755.19996, 738.1305908,
    716, 700.78878, 707.6685787, 700, 703.9170541,
    701.2203147, 699, 704.0910172, 701, 706.0509132,
    705.35, 702, 703.3818541, 713, 705.6608689,
    701.9512391, 711.9000928, 710, 713.7865827, 722.8688513,
    729.2058882, 740, 742.475943, 746, 748.6008856,
    750.7296288, 760, 750.7372196, 755.7544464, 761,
    759.1474261, 757.9121455, 758.2744076, 761.5,
    757.95, 761.5921327, 755.7726235, 761.5, 752.1919898,
    757.0571252, 748.6111, 743.5970352, 761.5, 751.9226969,
    745.7508941, 752, 761.5291335, 748, 751.051037,
    771.8169362, 792.59026, 747, 746, 868, 882.1841289,
    908.5, 932.9572053, 956.5, 974.5672105, 985, 987.05, 987.00111]

γ1_data_γ1 = [0.0048, 0.0042, 0.0046, 3.90E-03, 0.0056,
    0.009, 0.011, 0.029, 0.057, 0.073, 0.09, 1.13E-01,
    0.13, 3.00E-01, 0.29, 0.45, 6.14E-01, 0.66, 0.79,
    0.7875]

x_data_γ1 = [0.015, 0.02, 0.088, 0.1, 0.168, 0.2,
    0.252, 0.3, 0.325, 0.343, 0.398, 0.4,
    0.469, 0.5, 0.563, 0.6, 0.7, 0.71, 0.792,
    0.8]

T_data_γ1 = [1073., 1073., 1073., 1073., 1073.,
    1073., 1073., 1073., 1073., 1073., 1073.,
    1073., 1073., 1073., 1073., 1073., 1073.,
    1073., 1073., 1073.]

x_data_Hex = [0.0608, 0.1059, 0.2047,
    0.2602, 0.3524, 0.4495, 0.53,
    0.623, 0.6577, 0.8147, 0.8816,
    0.9415]

Hex_data_Hex = [-3313.728, -5614.928,
    -10447.448, -12991.32, -15627.24,
    -15627.24, -15008.008, -13167.048,
    -12811.408, -8393.104, -5640.032,
    -3016.664]

T_data_Hex = [1073., 1073., 1073., 1073.,
    1073., 1073., 1073., 1073., 1073.,
    1073., 1073., 1073.]

  # Pure Comp data from Table I of Pelton and Chartrand, Modified Quasichemical I
  KCls = MoltenSaltThermodynamics.PureComp([-436684.1, 82.55023, 40.01578, 25.46801, 3.64845, 0.,0.])
  KCll = MoltenSaltThermodynamics.PureComp([-421824.9, 86.5525, 73.59656, 0.,0.,0.,0.])
  MgCl2s = MoltenSaltThermodynamics.PureComp([-641616.0, 89.629, 54.58434, 21.42127, -11.12119, 399.1767, -2.356672])
                     MgCl2l = MoltenSaltThermodynamics.PureComp([-606887.4, 117.29708, 92.048, 0.,0.,0.,0.])
                     K2MgCl4s = MoltenSaltThermodynamics.PureComp([-1550013.0, 216.8, 263.05,0.,0.,0.,0.])
                     KMgCl3s = MoltenSaltThermodynamics.PureComp([-1100924.0, 162.3,157.65,0.,0.,0.,0.0])
  Temp = 1073.0
  mygab0 = -17947.0
  mygai = [-1026.0]
  mygbj = [-14801.0]
  mycoordparams = [6.0,6.0,3.0,6.0]
  KClMgCl2_model = BinaryMQCModel(mygab0, mygai, mygbj, mycoordparams)
  initguess = [0.990, 0.005, 0.005]
  my_mfracs = [[x, 1-x] for x in 0.999:-0.001:0.001]
  myfunction(mf, pfg=missing) = find_configuration2(KClMgCl2_model, mf,[KCll, MgCl2l], Temp, pfg)
  manyans = Vector{Any}(undef, length(my_mfracs))
  manyans[1] = myfunction(my_mfracs[1], initguess)
  max_nab = mycoordparams[1]/2
  for i in 2:length(manyans)
    manyans[i] = myfunction(my_mfracs[i], manyans[i-1].zero)
  end
  coordnums = [coordination(KClMgCl2_model, manyans[i].zero) for i in 1:length(my_mfracs)]
  totpairs = [sum(coordnums[i] .* my_mfracs[i])/2 for i in 1:length(my_mfracs)]
  Xminims = [x.zero for x in manyans]
  npairspace = Xminims .* totpairs
  xaa = [thing[1] for thing in Xminims]
  xbb = [thing[2] for thing in Xminims]
  xab = [thing[3] for thing in Xminims]
  naa = [thing[1] for thing in npairspace]
  nbb = [thing[2] for thing in npairspace]
  nab = [thing[3] for thing in npairspace]
  Xaspace = [my_mfracs[i][1] for i in 1:length(my_mfracs)]
  # Yes, I know that I don't have to use a different internal variable
  # for the following functions. I just want to.
  mygibbs(i) = gibbs_binary_solution(KClMgCl2_model, my_mfracs[i], Xminims[i], [KCll, MgCl2l], Temp)
  myent(m) = Sconfig(KClMgCl2_model, my_mfracs[m], Xminims[m])
  mygibbsex(j) = gibbs_excess_eq35(KClMgCl2_model, [naa[j], nbb[j], nab[j]])
  mygibbsex2(k) = gibbs_excess2(KClMgCl2_model, my_mfracs[k], Xminims[k], [KCll, MgCl2l], Temp)
  mygibbsexc1(k) = gibbs_excess_comp2(KClMgCl2_model,[KCll, MgCl2l], my_mfracs[k], npairspace[k],  Temp)
  mygibbsexc2(k) = gibbs_excess_comp2(KClMgCl2_model,[KCll, MgCl2l], my_mfracs[k], npairspace[k], Temp,2)
  myenthex2(l) = enthalpy_excess(KClMgCl2_model, my_mfracs[l], Xminims[l], [KCll, MgCl2l], Temp)
  myact(m) = activity(KClMgCl2_model, [KCll, MgCl2l], my_mfracs[m], npairspace[m], Temp)
  gibbsltot = [mygibbs(i) for i in 1:length(manyans)]
  entconfig = [myent(m) for m in 1:length(manyans)]
  gibbsl = [mygibbsex(j) for j in 1:length(manyans)]
  gibbsx2 = [mygibbsex2(k) for k in 1:length(manyans)]
  gibbsxc1 = [mygibbsexc1(k) for k in 1:length(manyans)]
  gibbsxc2 = [mygibbsexc2(k) for k in 1:length(manyans)]
  enthex = [myenthex2(l) for l in 1:length(manyans)]
  activkcl = [myact(m) for m in 1:length(manyans)]
  #println(activkcl)
  #println(length(gibbsl))
  sc1 = plot(xab, label="Xab")
  plot!(sc1, xaa, label="Xaa")
  plot!(sc1, xbb, label="Xbb")
  plot!(sc1, naa, label="naa")
  plot!(sc1, nab, label="nab")
  plot!(sc1, nbb, label="nbb")
  savefig(sc1, "/home/doomphoenix/ENTER/Research/MQCM/testfig1.png")
  display(sc1)
  #println(my_mfracs)
  #print(Xaspace)
  sc2 = plot(Xaspace ,gibbsl, label="GEx by Eqn 35")
  plot!(sc2, Xaspace, gibbsx2, label="GEx by definition")
  plot!(sc2, Xaspace, enthex, label="Hex")
  scatter!(sc2, x_data_Hex, Hex_data_Hex, label="Hex data")
  savefig(sc2, "/home/doomphoenix/ENTER/Research/MQCM/testfig2.png")
  ##display(sc2)
  #println("Gibbs of solid KCl at 1400: ", gibbsenergy(KCls, 1400.0))
  #println("Gibbs of solid MgCl2 at 1400: ", gibbsenergy(MgCl2s, 1400.0))
  sc3 = plot(Xaspace, gibbsltot, label="Total Gibbs liquid")
  savefig(sc3, "/home/doomphoenix/ENTER/Research/MQCM/testfig3.png")
  ##display(sc3)
  sc4 = plot(Xaspace, entconfig, label="Config. entropy")
  savefig(sc4, "/home/doomphoenix/ENTER/Research/MQCM/testfig4.png")
  ##display(sc4)
  sc5 = plot(Xaspace ,activkcl, label="Act. KCl calc")
  plot!(sc5, [0.0,1.0],[1.0,0.0])
  scatter!(sc5, x_data_a2, a2_data_a2, label="Activ. data KCl")
  savefig(sc5, "/home/doomphoenix/ENTER/Research/MQCM/testfig5.png")
  sc6 = plot(Xaspace ,gibbsxc1, label="GEx KCl")
  plot!(sc6, Xaspace, gibbsxc2, label="GEx MgCl2")
  savefig(sc6, "/home/doomphoenix/ENTER/Research/MQCM/testfig6.png")
  println("Gibbs of total liquid at mole frac ", Xaspace[500], " KCl: ", gibbsltot[500])
  println("Gibbs excess of total liquid at mole frac ", Xaspace[500], " KCl: ", gibbsx2[500])
  println("Entropy of total liquid at mole frac ", Xaspace[500], " KCl: ", gibbsltot[500])
  println("Config Entropy of total liquid at mole frac ", Xaspace[500], " KCl: ", entconfig[500])
end


text:
using JetPackWaveFD, Random, Test

function test_ginsu_sub(g::Ginsu{2}, x, xsub, interior, extend)
    N1 = size(x,1)
    N2 = size(x,2)

    I1 = lextents(g, interior=interior)[1]
    I2 = lextents(g, interior=interior)[2]

    for i2 = 1:N2, i1 = 1:N1
        if in(i1, I1) && in(i2, I2)
            @test x[i1,i2] ≈ xsub[i1-I1[1]+1,i2-I2[1]+1]
        end
    end
end

function test_ginsu_sub(g::Ginsu{3}, x, xsub, interior, extend)
    N1 = size(x,1)
    N2 = size(x,2)
    N3 = size(x,3)

    I1 = lextents(g, interior=interior)[1]
    I2 = lextents(g, interior=interior)[2]
    I3 = lextents(g, interior=interior)[3]

    for i3 = 1:N3, i2 = 1:N2, i1 = 1:N1
        if in(i1, I1) && in(i2, I2) && in(i3, I3)
            @test x[i1,i2,i3] ≈ xsub[i1-I1[1]+1,i2-I2[1]+1,i3-I3[1]+1]
        end
    end
end

function test_ginsu_super(g::Ginsu{2}, x, xsuper, interior)
    N1 = size(x,1)
    N2 = size(x,2)

    I1 = lextents(g, interior=interior)[1]
    I2 = lextents(g, interior=interior)[2]

    for i2 = 1:N2, i1 = 1:N1
        if in(i1,I1) && in(i2,I2)
            @test x[i1,i2] ≈ xsuper[i1,i2]
        else
            @test xsuper[i1,i2] ≈ 0
        end
    end
end

function test_ginsu_super(g::Ginsu{3}, x, xsuper, interior)
    N1 = size(x,1)
    N2 = size(x,2)
    N3 = size(x,3)

    I1 = lextents(g, interior=interior)[1]
    I2 = lextents(g, interior=interior)[2]
    I3 = lextents(g, interior=interior)[3]

    for i3 = 1:N3, i2 = 1:N2, i1 = 1:N1
        if in(i1,I1) && in(i2,I2) && in(i3,I3)
            @test x[i1,i2,i3] ≈ xsuper[i1,i2,i3]
        else
            @test xsuper[i1,i2,i3] ≈ 0
        end
    end
end

@testset "Ginsu 2D" begin
    nz, nx = 100, 200
    dz, dx = 10.0, 20.0
    z0, x0 = 0.0, 0.0

    sz, sx = [10.0], [2000.0]
    rz, rx = [0.0, 0.0], [0.0, 1990.0]
    padz, padx = 100.0, 200.0
    ndamp = 10

    ginsu = Ginsu((z0,x0), (dz,dx), (nz,nx), (sz,sx), (rz,rx), ((padz,padz),(padx,padx)), ((0,0),(ndamp,ndamp)), dims=(:z,:x), T=Float64)
    @test lextents(ginsu, interior=false) == (-9:110,-19:124)
    @test lextents(ginsu, interior=true) == (-9:110,-9:114)

    map(i->@test(pextents(ginsu, interior=false)[i] ≈ (-100.0:10.0:1090.0,-400.0:20.0:2460.0)[i]), 1:2)
    map(i->@test(pextents(ginsu, interior=true)[i] ≈ (-100.0:10.0:1090.0,-200.0:20.0:2260.0)[i]), 1:2)

    for interior in (false, true), extend in (false, true)
        x = rand(nz,nx)
        z = zeros(nz,nx)

        y = sub(ginsu, x, interior=interior, extend=extend)
        test_ginsu_sub(ginsu, x, y, interior, extend)

        super!(z, ginsu, y, interior=interior, accumulate=false)
        test_ginsu_super(ginsu, x, z, interior)

        super!(z, ginsu, y, interior=interior, accumulate=false)
        test_ginsu_super(ginsu, x, z, interior)

        super!(z, ginsu, y, interior=interior, accumulate=true)
        test_ginsu_super(ginsu, 2x, z, interior)
    end
end

@testset "Ginsu 2D Array Src&Rec" begin
    nz, nx = 100, 200
    dz, dx = 10.0, 20.0
    z0, x0 = 0.0, 0.0

    sz, sx = [[10.0]], [[2000.0]]
    rz, rx = [[0.0, 0.0]], [[0.0, 1990.0]]
    padz, padx = 100.0, 200.0
    ndamp = 10

    ginsu = Ginsu((z0,x0), (dz,dx), (nz,nx), (sz,sx), (rz,rx), ((padz,padz),(padx,padx)), ((0,0),(ndamp,ndamp)), dims=(:z,:x), T=Float64)
    @test lextents(ginsu, interior=false) == (-9:110,-19:124)
    @test lextents(ginsu, interior=true) == (-9:110,-9:114)

    map(i->@test(pextents(ginsu, interior=false)[i] ≈ (-100.0:10.0:1090.0,-400.0:20.0:2460.0)[i]), 1:2)
    map(i->@test(pextents(ginsu, interior=true)[i] ≈ (-100.0:10.0:1090.0,-200.0:20.0:2260.0)[i]), 1:2)

    for interior in (false, true), extend in (false, true)
        x = rand(nz,nx)
        z = zeros(nz,nx)

        y = sub(ginsu, x, interior=interior, extend=extend)
        test_ginsu_sub(ginsu, x, y, interior, extend)

        super!(z, ginsu, y, interior=interior, accumulate=false)
        test_ginsu_super(ginsu, x, z, interior)

        super!(z, ginsu, y, interior=interior, accumulate=false)
        test_ginsu_super(ginsu, x, z, interior)

        super!(z, ginsu, y, interior=interior, accumulate=true)
        test_ginsu_super(ginsu, 2x, z, interior)
    end
end

@testset "Ginsu 3D" begin
    nz, ny, nx = 100, 30, 200
    dz, dy, dx = 10.0, 10.0, 20.0
    z0, y0, x0 = 0.0, 0.0, 0.0

    sz, sy, sx = [10.0], [150.0], [2000.0]
    rz, ry, rx = [0.0, 0.0], [0.0,270.0], [0.0, 1990.0]
    padz, pady, padx = 100.0, 50.0, 200.0
    ndamp = 10

    ginsu = Ginsu((z0,y0,x0), (dz,dy,dx), (nz,ny,nx), (sz,sy,sx), (rz,ry,rx), ((padz,padz),(pady,pady),(padx,padx)), ((0,0),(ndamp,ndamp),(ndamp,ndamp)), dims=(:z,:y,:x), T=Float64)
    @test lextents(ginsu, interior=false) == (-9:110,-14:49,-19:124)
    @test lextents(ginsu, interior=true) == (-9:110,-4:39,-9:114)

    map(i->@test(pextents(ginsu, interior=false)[i] ≈ (-100.0:10.0:1090.0,-150.0:10.0:480.0,-400.0:20.0:2460.0)[i]), 1:2)
    map(i->@test(pextents(ginsu, interior=true)[i] ≈ (-100.0:10.0:1090.0,-50.0:10.0:380.0,-200.0:20.0:2260.0)[i]), 1:2)

    for interior in (false, true), extend in (false, true)
        x = rand(nz,ny,nx)
        z = zeros(nz,ny,nx)

        y = sub(ginsu, x, interior=interior, extend=extend)
        test_ginsu_sub(ginsu, x, y, interior, extend)

        super!(z, ginsu, y, interior=interior, accumulate=false)
        test_ginsu_super(ginsu, x, z, interior)

        super!(z, ginsu, y, interior=interior, accumulate=false)
        test_ginsu_super(ginsu, x, z, interior)

        super!(z, ginsu, y, interior=interior, accumulate=true)
        test_ginsu_super(ginsu, 2x, z, interior)
    end
end

@testset "Ginsu 3D Array Src&Rec" begin
    nz, ny, nx = 100, 30, 200
    dz, dy, dx = 10.0, 10.0, 20.0
    z0, y0, x0 = 0.0, 0.0, 0.0

    sz, sy, sx = [[10.0]], [[150.0]], [[2000.0]]
    rz, ry, rx = [[0.0, 0.0]], [[0.0,270.0]], [[0.0, 1990.0]]
    padz, pady, padx = 100.0, 50.0, 200.0
    ndamp = 10

    ginsu = Ginsu((z0,y0,x0), (dz,dy,dx), (nz,ny,nx), (sz,sy,sx), (rz,ry,rx), ((padz,padz),(pady,pady),(padx,padx)), ((0,0),(ndamp,ndamp),(ndamp,ndamp)), dims=(:z,:y,:x), T=Float64)
    @test lextents(ginsu, interior=false) == (-9:110,-14:49,-19:124)
    @test lextents(ginsu, interior=true) == (-9:110,-4:39,-9:114)

    map(i->@test(pextents(ginsu, interior=false)[i] ≈ (-100.0:10.0:1090.0,-150.0:10.0:480.0,-400.0:20.0:2460.0)[i]), 1:2)
    map(i->@test(pextents(ginsu, interior=true)[i] ≈ (-100.0:10.0:1090.0,-50.0:10.0:380.0,-200.0:20.0:2260.0)[i]), 1:2)

    for interior in (false, true), extend in (false, true)
        x = rand(nz,ny,nx)
        z = zeros(nz,ny,nx)

        y = sub(ginsu, x, interior=interior, extend=extend)
        test_ginsu_sub(ginsu, x, y, interior, extend)

        super!(z, ginsu, y, interior=interior, accumulate=false)
        test_ginsu_super(ginsu, x, z, interior)

        super!(z, ginsu, y, interior=interior, accumulate=false)
        test_ginsu_super(ginsu, x, z, interior)

        super!(z, ginsu, y, interior=interior, accumulate=true)
        test_ginsu_super(ginsu, 2x, z, interior)
    end
end

@testset "Ginsu Vector Width" begin
    nz, nx = 100, 200
    dz, dx = 10.0, 20.0
    z0, x0 = 0.0, 0.0

    sz, sx = [10.0], [2000.0]
    rz, rx = [0.0, 0.0], [0.0, 1990.0]
    padz, padx = 100.0, 200.0
    ndamp = 10
    for vector_width = (1,2,4,8,16,32)
        ginsu = Ginsu((z0,x0), (dz,dx), (nz,nx), (sz,sx), (rz,rx), ((padz,padz),(padx,padx)), ((0,0),(ndamp,ndamp)), dims=(:z,:x), T=Float64, vector_width=vector_width)
        @test size(ginsu) .% vector_width == (0,0)
    end
end

@testset "Ginsu dot product" begin
    rz = 750.0*ones(size(75:125))
    rx = 10.0*collect(75:125)
    ginsu = Ginsu((0.0,0.0), (10.0,10.0), (201,201), ([1000.0],[1000.0]), (rz,rx), ((500.0,500.0),(500.0,500.0)), ((50,50),(50,50)), T=Float64)
    m = rand(201,201)
    d = rand(size(ginsu)...)
    ds = sub(ginsu, m, extend=false)
    ms = zeros(201,201)
    super!(ms, ginsu, d)
    lhs = sum(d .* ds)
    rhs = sum(m .* ms)
    err = (lhs - rhs) / (lhs + rhs)
    @test isapprox(err, 0.0, atol=100*eps(Float32))

    rz = 10.0*ones(size(1:201))
    rx = 10.0*collect(1:201)
    ginsu = Ginsu((0.0,0.0), (10.0,10.0), (201,201), ([10.0],[1000.0]), (rz,rx), ((500.0,500.0),(500.0,500.0)), ((50,50),(50,50)), T=Float64)
    m = rand(201,201)
    d = rand(size(ginsu)...)
    ds = sub(ginsu, m, extend=false)
    ms = zeros(201,201)
    super!(ms, ginsu, d)
    lhs = sum(d .* ds)
    rhs = sum(m .* ms)
    @test lhs ≈ rhs rtol=1e-5
end

@testset "Ginsu, absolute constructor" begin
    ginsu = Ginsu((0.0,0.0), (10.0,10.0), (201,201), (-100.0:200.0,-50.0:150.0), ((50,50),(50,50)))
    @test lextents(ginsu, interior=false) == (-59:76,-54:73)
    @test lextents(ginsu, interior=true) == (-9:26,-4:23)

    map(i->@test(pextents(ginsu, interior=false)[i] ≈ (-600.0:10.0:750.0,-550.0:10.0:720.0)[i]), 1:2)
    map(i->@test(pextents(ginsu, interior=true)[i] ≈ (-100.0:10.0:250.0,-50.0:10.0:220.0)[i]), 1:2)
end

max_stars_repo_name: ChevronETC/JetPackWaveFD

max_issues_repo_name: ChevronETC/JetPackWaveFD

max_forks_repo_name: ChevronETC/JetPackWaveFD