text:
"""
`optimalPLD_createRanking_rfn(n, revenues, numPermutations, orderings, lambda)`

Create a function that will compute the expected per-customer profit/revenue of a product line according to the first-choice model specified by the input arguments.

`n` = number of candidate products.

`revenues` = array of n numbers representing marginal revenue/profit of each of the n products. (These are the pi_i's in the Bertsimas and Mišić paper.)

`numPermutations` = number of rankings / customer types.

`orderings` = 2D array of integers encoding each ranking. orderings[k,j] contains the index of ranking k's jth top product; (n+1) encodes the no-purchase option. For example, if n = 5 and

	orderings[k,:] = [3,5,2,6,1,4]

then the customer prefers 3 < 5 < 2 < 0 (no-purchase option).

`lambda` = array specifying the probability of each ranking / customer type; should be nonnegative and sum to 1.

"""
function optimalPLD_createRanking_rfn(n, revenues, numPermutations, orderings, lambda)
	unpermed_orderings = zeros(Int64,size(orderings))
	for k=1:numPermutations
		unpermed_orderings[k,:] = sortperm(orderings[k,:][:])
	end
	revenue_fn = function ( s::Array{Int64} )
		s0 = [s;n+1];
		temp = 0;
		# @show revenues
		for k=1:numPermutations
			# @show s
			s_ranks = unpermed_orderings[k,s0[:]]
			# @show s_ranks
			ind = indmin(s_ranks)
			temp = temp + revenues[s0[ind]] * lambda[k]
			# @show k
			# @show
		end

		return temp
	end

	return revenue_fn

end

max_stars_repo_name: vvmisic/EFCPLO

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_name: vvmisic/EFCPLO

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_name: vvmisic/EFCPLO

max_forks_repo_licenses: MIT Unlicense


text:
using Compat
@compat readstring
# Markdown Report Generation #
#----------------------------#

function printreport(io::IO, results)
    # print report preface + job properties #
    #---------------------------------------#
    println(io, """
    # Benchmark Report
    The percentages accompanying time and memory values in the below table are noise tolerances. The "true"
    time/memory value for a given benchmark is expected to fall within this percentage of the reported value.
    """)

    # print result table #
    #--------------------#
    print(io, """
    | ID | time | GC time | memory | allocations |
    |----|------|---------|--------|-------------|
    """)

    entries = BenchmarkTools.leaves(median(results))

    try
        entries = entries[sortperm(map(x -> string(first(x)), entries))]
    end

    for (ids, t) in entries
        println(io, resultrow(ids, t))
    end

    println(io)

    # print list of executed benchmarks #
    #-----------------------------------#
    println(io, """
    ## Benchmark Group List
    Here's a list of all the benchmark groups executed by this job:
    """)

    for id in unique(map(pair -> pair[1][1:end-1], entries))
        println(io, "- `", idrepr(id), "`")
    end

    println(io)

    # print build version info #
    #--------------------------#
    print(io, """
    ## Version Info
    #### Module Commit ID
    ```
    $(chomp(readstring(`git rev-parse HEAD`)))
    ```
    #### Julia Version
    ```
    $(first(split(chomp(readstring(`julia -e 'versioninfo(true)'`)), "Environment")))
    ```
    """)
end

idrepr(id) = (str = repr(id); str[searchindex(str, '['):end])

function resultrow(ids, t::BenchmarkTools.TrialEstimate)
    t_tol = BenchmarkTools.prettypercent(BenchmarkTools.params(t).time_tolerance)
    m_tol = BenchmarkTools.prettypercent(BenchmarkTools.params(t).memory_tolerance)
    timestr = string(BenchmarkTools.prettytime(BenchmarkTools.time(t)), " (", t_tol, ")")
    memstr = string(BenchmarkTools.prettymemory(BenchmarkTools.memory(t)), " (", m_tol, ")")
    gcstr = BenchmarkTools.prettytime(BenchmarkTools.gctime(t))
    allocstr = string(BenchmarkTools.allocs(t))
    return "| `$(idrepr(ids))` | $(timestr) | $(gcstr) | $(memstr) | $(allocstr) |"
end

ext: jl


text:
# load pkgs
using Ipopt, HarmonicPowerModels, PowerModels

# pkg const
const _PMs = PowerModels
const _HPM = HarmonicPowerModels

# path to the data
path = joinpath(_HPM.BASE_DIR,"test/data/matpower/case_xfmr_YNyn0_simplified.m")

# transformer excitation data
xfmr = Dict("voltage_harmonics" => [1,3],
            "current_harmonics" => [1,3],
            "N" => 50,
            "current_type" => :rectangular,
            "excitation_type" => :sigmoid,
            "inom" => 0.4,
            "ψmax" => 0.5,
            "voltage_type" => :rectangular,
            "dv" => [0.1,0.1],
            "vmin" => [0.0,0.0],
            "vmax" => [1.1,1.1],
            "dθ" => [π/5,π/5],
            "θmin" => [0.0,0.0],
            "θmax" => [2π,2π])

# load data
data  = _PMs.parse_file(path)

vmmin = 0.8
for (b,branch) in data["branch"]
    branch["c_rating"] = branch["rate_a"]/vmmin
end
for (d,load) in data["load"]
    load["c_rating"] = abs(load["pd"] + im* load["qd"])/vmmin
end
for (g,gen) in data["gen"]
    gen["c_rating"] = abs(gen["pmax"] + im* gen["qmax"])/vmmin
end

hdata = _HPM.replicate(data)
# hdata = _HPM.replicate(data, xfmr_exc=xfmr)

# set the solver
solver = Ipopt.Optimizer

# solve the hopf
# run_hopf_iv(hdata, _PMs.IVRPowerModel, solver)
pm = _PMs.instantiate_model(hdata, _PMs.IVRPowerModel, _HPM.build_hopf_iv; ref_extensions=[_HPM.ref_add_xfmr!]);
result = optimize_model!(pm, optimizer=solver, solution_processors=[ _HPM.sol_data_model!])
_HPM.append_indicators!(result, hdata)

for (n,nw) in result["solution"]["nw"]
    for (i,bus) in nw["bus"]
            bus["vm"] =  abs(bus["vr"] +im*  bus["vi"])
            bus["va"] =  angle(bus["vr"] +im*  bus["vi"])*180/pi
    end
end
println("Harmonic 3")
_PMs.print_summary(result["solution"]["nw"]["2"])
println("Harmonic 1")
_PMs.print_summary(result["solution"]["nw"]["1"])
result["objective"]
result["termination_status"]

##
# Yy -> third harmonic of the excitation cannot flow -> current is inherently forced to 0
# this means the third harmonic voltage is 0
#


text:
# No aggregation
# ------------------------------------------------------------
"""
    NoAggregation

Empty functor object for running an L-shaped algorithm without aggregation (multi-cut L-shaped).

"""
struct NoAggregation <: AbstractAggregation end

function aggregate_cut!(lshaped::AbstractLShaped, ::NoAggregation, cut::HyperPlane)
    return add_cut!(lshaped, cut)
end

function aggregate_cut!(cutqueue::CutQueue, ::NoAggregation, ::MetaData, t::Integer, cut::HyperPlane, x::AbstractArray)
    put!(cutqueue, (t, cut))
    return nothing
end

function num_thetas(num_subproblems::Integer, ::NoAggregation)
    return num_subproblems
end

function num_thetas(num_subproblems::Integer, ::NoAggregation, ::AbstractScenarioProblems)
    return num_subproblems
end

function flush!(::AbstractLShaped, ::NoAggregation)
    return false
end

function flush!(::CutQueue, ::NoAggregation, ::MetaData, ::Integer, ::AbstractArray)
    return false
end

# API
# ------------------------------------------------------------
"""
    DontAggregate

Factory object for [`NoAggregation`](@ref). Passed by default to `aggregate` in the `LShapedSolver` factory function.

"""
struct DontAggregate <: AbstractAggregator end

function (::DontAggregate)(::Integer, ::Type{<:AbstractFloat})
    return NoAggregation()
end

function remote_aggregator(::NoAggregation, ::AbstractScenarioProblems, ::Integer)
    return DontAggregate()
end

function str(::DontAggregate)
    return "disaggregate cuts"
end


text:
using GraphBlasAlgorithmsFinal
using Test


text:
# include("../utils.jl")
include("RVRP_data.jl")
include("RVRP_solution.jl")
include("utils.jl")
include("parser_cvrplib.jl")
include("instance_check.jl")
include("parser_json.jl")
include("solution_check.jl")

include("mock_solver.jl")
include("vroom.jl")
include("jsprit.jl")

function unit_tests()
    @testset "RVRP_data.jl" begin
        rvrp_data_unit_tests()
    end
    @testset "RVRP_solution.jl" begin
        rvrp_solution_unit_tests()
    end
    @testset "utils.jl" begin
        utils_unit_tests()
    end
    @testset "parser_cvrplib.jl" begin
        parser_cvrplib_unit_tests()
    end
    @testset "instance_check.jl" begin
        instance_check_unit_tests()
    end
    @testset "parser_json.jl" begin
        parser_json_unit_tests()
    end
    @testset "solution_check.jl" begin
        solution_check_unit_tests()
    end

    @testset "vroom_cpp.jl" begin
        vroom_cpp_unit_tests()
    end
    @testset "jsprit_java.jl" begin
        jsprit_java_unit_tests()
    end
    @testset "mock_solver.jl" begin
        mock_solver_unit_tests()
    end
end


text:
cd(@__DIR__); include("setups/path512.jl");
gr(dpi = 200)

## Build VM-NGWP
@time VM_NGWP = vm_ngwp(𝚽, GP_dual)

#################### Fig.5
j = 5
for k in [1, 2, 5]
   WW = sort_wavelets(VM_NGWP[GP_dual.rs[k, j]:(GP_dual.rs[k + 1, j] - 1), j, :]')
   sc = 0.75
   if k == 1
      sc = 0.5
   end
   if k == 2
      WW[:, end] *= -1
   end
   plt = wiggle(WW; sc = sc)
   savefig(plt, "../figs/Path512_VM_NGWP_j$(j-1)k$(k-1).png")
end
current()


text:
module Polygons

using Gadfly
using Colors
using FixedSizeArrays
using CHull2d
import Base: *, .*, +, union

export
    # Types
    Polygon,

    # Methods
    acwOrder,
    ptsToDirs,
    dirsToPts,
    (*),
    (.*),
    (+),
    setSum,
    weightedSum,
    grahamScan,
    chull,
    crop,
    polyPlot,
    deeDoop,
    hausdorff

# TODO: Hausdorff distance

## Source files
include("common.jl")
include("simple.jl")
include("convex.jl")
include("utils.jl")
include("hausdorff.jl")
include("polyPlot.jl")

end # module


text:
using Random, LinearAlgebra, Krylov
using Plots
using RandomizedPreconditioners
const RP = RandomizedPreconditioners

n = 500
μ = 1e-4

# Build Spectrum (exp decay)
λs = max.(√n * exp.(-(0:n-1) ./ √n), μ)
plot(λs, lw=3, legend=false, title="Spectrum")

# Build Matrix
V = Array(qr(randn(n,n)).Q)
A = V*Diagonal(λs)*V'
A = 0.5(A + A')

# Build system
xtrue = randn(n)
b = A * xtrue

# Build preconditioner
r = 110
k = round(Int, r - 5)
Anys = NystromSketch(A, k, r)
P = RP.NystromPreconditionerInverse(Anys, μ)

# Solve system
_, stats_npc = cg(A, b; history=true)
_, stats_nys = cg(A, b; history=true, M = P)
_, stats_diag = cg(A, b; history=true, M=Diagonal(diag(A)))

nb = norm(b)
res_npc = prepend!(stats_npc.residuals, nb)
res_nys = prepend!(stats_nys.residuals, nb)
res_diag = prepend!(stats_diag.residuals, nb)

plt_cg = plot(
    res_npc,
    dpi=300,
    lw=2,
    label="No Preconditioner",
    ylabel="residual",
    xlabel="iteration",
    title="Convergence of CG",
    legend=:topright,
    yaxis=:log
)
plot!(plt_cg,
    res_nys,
    label="Nystrom Preconditioner",
    lw=2
)
plot!(plt_cg,
    res_diag,
    label="Diagonal Preconditioner",
    lw=2
)
display(plt_cg)


text:
using Test
using DynamicPipe

# basic tests
@testset "@>, @>> pipeing into one function" begin
    x = [1,2,3]

    y = x |>
        @>  sum
    @test y == sum(x)

    y = x |>
        @> sum()
    @test y == sum(x)

    y = x |>
        @> sum(_)
    @test y == sum(x)

   y = @>> x |> sum
   @test y == sum(x)

   y = @>> x |> sum()
   @test y == sum(x)

   y = @>> x |> sum(_)
   @test y == sum(x)

end

macro double(x)
    return :($(esc(:($x + $x))))
end

@testset "@>, @>> test pipeing into macros" begin
    x = [1,2,3]
    y = x + x

    t1 = x |>
        @>  @double()
    t2 = x |>
        @>  @double
    t3 = x |>
        @>  @double(_)

    @test y == t1
    @test y == t2
    @test y == t3

    t1 = @>> x |> @double
    t2 = @>> x |> @double()
    t3 = @>> x |> @double(_)

    @test y == t1
    @test y == t2
    @test y == t3
end

@testset "@>, @>> test first argument imputation" begin
    x = [1, 2, 3]

    y = copy(x) |>
        @>  append!(4)
    @test y == [1,2,3,4]

    y = @>> copy(x) |>
        append!(4)
    @test y == [1, 2, 3, 4]
end

@testset "@>, @>> test multiple pipes" begin
    x = [1, 2, 3]

    y = x |>
        @>  copy |>
            append!(4)
    @test y == [1,2,3,4]

    y = @>> x |>
        copy |>
        append!(4)

    @test y == [1, 2, 3, 4]
end

@testset "@>, @>> test vectorisation" begin
    x = [1,2,3]

    y = x |>
        @>  copy() |>
            .+(1)
    @test y == [2,3,4]

    y = @>> x |>
        copy() |>
        .+(1)

    @test y == [2,3,4]

    x = [1,2,3]

    y = x |>
        @>  copy .|>
            +(1)

    @test y == [2,3,4]

    y = @>> x |>
        copy .|>
        +(1)

    @test y == [2,3,4]

end

@testset "@>, @>> test pipes in pipes" begin
    x = [1,2,3]
    z = [1,1,1]
    y_ans = (
        (x + ((x .+ 1) .* 2)) + z .- 1
    )

    # TODO fix
    # make 2 versions of create create_anonymous_function
    # one that also acts on second part and pipe func
    # problem is this :(x |> @> _  .+ 1 |> _ .* 2).args
    # we are creating a function for the macro second part
    y = x |>
        @>  +(x |>
                @>  _ .+ 1 |>
                    _ .* 2) |>
            _ + z |>
            .-(1)

    @test y == y_ans

    y = @>> x |>
        +(x |>
            @>  _ .+ 1 |>
                _ .* 2) |>
        _ + z |>
        .-(1)

    @test y == y_ans

    y = @>> x |>
        +(@>> x |>
            _ .+ 1 |>
            _ .* 2) |>
        _ + z |>
        .-(1)

    @test y == y_ans

    y = @>> x |>
        +(@>> _ |>
            _ .+ 1 |>
            _ .* 2) |>
        _ + z |>
        .-(1)

    @test y == (
        (((x .+ 1) .* 2) + z) .- 1
    )

    y = @>> x |>
        _ + @>> _ |>
            _ .+ 1 |>
            _ .* 2 |>
        _ + z |>
        .-(1)

    @test y == y_ans

    y = @>> x |>
            +(_ |>
                @>  _ .+ 1 |>
                    _ .* 2) |>
            _ + z |>
            .-(1)

    @test y == (
        (((x .+ 1) .* 2) + z) .- 1
    )

    y = x |>
        @>  +(_ |>
                @>  _ .+ 1 |>
                    _ .* 2) |>
            _ + z |>
            .-(1)

    @test y == (
        (((x .+ 1) .* 2) + z) .- 1
    )

    y = x |>
        @>  _ + (_ |>
                    @>  _ .+ 1 |>
                        _ .* 2) |>
            _ + z |>
            .-(1)

    @test y == y_ans
end

@testset "@>, @>> test block input" begin
    x = [1,2,3]
    z = [1,1,1]
    y_ans = (
        (x + ((x .+ 1) .* 2)) + z .- 1
    )

    y = x |>
        @> begin
            +(x |>
                @>  _ .+ 1 |>
                    _ .* 2)
            _ + z
            .-(1)
        end

    @test y == y_ans

    y = x |>
        @> begin
            +(x |>
                @> begin
                    _ .+ 1
                    _ .* 2
                end)
            _ + z
            .-(1)
        end

    @test y == y_ans

    y = @>> begin
            x
            +(x |>
                @>  _ .+ 1 |>
                    _ .* 2)
            _ + z
            .-(1)
        end

    @test y == y_ans

    y = @>> begin
            x
            +(_, _ |>
                @>  _ .+ 1 |>
                    _ .* 2)
            _ + z
            .-(1)
        end

    @test y == y_ans

    y = @>> x begin
            +(x |>
                @>  _ .+ 1 |>
                    _ .* 2)
            _ + z
            .-(1)
        end

    @test y == y_ans

    y = @>> x begin
            +(_, _ |>
                @>  _ .+ 1 |>
                    _ .* 2)
            _ + z
            .-(1)
        end

    @test y == y_ans

    y = @>> x begin
        +(@>> begin x
            _ .+ 1
            _ .* 2
        end)
        _ + z
        .-(1)
    end

    @test y == y_ans

    y = @>> x begin
        +(@>> x begin
            _ .+ 1
            _ .* 2
        end)
        _ + z
        .-(1)
    end

    @test y == y_ans

    y = @>> x begin
        +(_, @>> begin
            _
            _ .+ 1
            _ .* 2
        end)
        _ + z
        .-(1)
    end

    @test y == y_ans

    y = @>> x begin
        +(_, @>> _ begin
            _ .+ 1
            _ .* 2
        end)
        _ + z
        .-(1)
    end

    @test y == y_ans
end

@testset "@>, @>> test pipes in pipes without macros calls" begin
    x = [1,2,3]
    z = [1,1,1]
    y_ans = (
        (x + ((x .+ 1) .* 2)) + z .- 1
    )

    y = x |>
        @> +(x |>
                _ .+ 1 |>
                _ .* 2) |>
            _ + z |>
            .-(1)

    @test y == y_ans

    y = @>> x |>
        +(x |>
            _ .+ 1 |>
            _ .* 2) |>
        _ + z |>
        .-(1)

    @test y == y_ans

    y = @>> x |>
        +(_ |>
            _ .+ 1 |>
            _ .* 2) |>
        _ + z |>
        .-(1)

    @test y == (
        (((x .+ 1) .* 2) + z) .- 1
    )

    y = @>> x |>
        _ + (_ |>
            _ .+ 1 |>
            _ .* 2 |>
        _ + z |>
        .-(1))

    @test y == y_ans

    y = @>> x |>
            +(_ |>
                _ .+ 1 |>
                _ .* 2) |>
            _ + z |>
            .-(1)

    @test y == (
        (((x .+ 1) .* 2) + z) .- 1
    )

    y = x |>
        @>  +(_ |>
                _ .+ 1 |>
                _ .* 2) |>
            _ + z |>
            .-(1)

    @test y == (
        (((x .+ 1) .* 2) + z) .- 1
    )

    x = [1,2,3]
    z = [1,1,1]
    y_ans = (
        (x + ((x .+ 1) .* 2)) + z .- 1
    )

    y = x |>
        @> begin
            +(x |>
                _ .+ 1 |>
                _ .* 2)
            _ + z
            .-(1)
        end

    @test y == y_ans

    y = @>> begin
            x
            +(x |>
                _ .+ 1 |>
                _ .* 2)
            _ + z
            .-(1)
        end

    @test y == y_ans

    y = @>> begin
            x
            +(_, _ |>
                _ .+ 1 |>
                _ .* 2)
            _ + z
            .-(1)
        end

    @test y == y_ans

    y = @>> x begin
            +(x |>
                _ .+ 1 |>
                _ .* 2)
            _ + z
            .-(1)
        end

    @test y == y_ans

    y = @>> x begin
            +(_, _ |>
                _ .+ 1 |>
                _ .* 2)
            _ + z
            .-(1)
        end

    @test y == y_ans

    y = x |>
        @>  _ + (_ .|> (_ + 1)) |>
            (_ .+ sum((_ |> .*(3)))) |>
            _ .+ (2 |> -(3))

    @test y == (x + (x .+ 1)) .+ sum((x + (x .+ 1)) .* 3) .+ (-1)

    y = x |>
        @>  _ .+ (2 |> -(3 |> +(6))) .|>
            _ + (7 |> -(1 |> +(2))) |>
            _ .- 2 .|>
            +(1)

    @test y == (
        x .+ (2 - 9) |>
        y -> y .+ (7 -3 ) |>
        y -> y .- 2 |>
        y -> y .+ 1
    )
end

max_stars_repo_name: ConallTimoney/DynamicPipe

max_issues_repo_name: ConallTimoney/DynamicPipes

max_forks_repo_name: ConallTimoney/DynamicPipes


text:
function first_term_eval(H::hypergraph,Z::Array{<:Integer,1}, Ω::IntensityFunction; α)
    """
    Not optimized, goal is to make this as quick and easy as
    possible using existing code.
    H: hypergraph
    Z: array storing cluster indices; Z[i] is the cluster node i is in
    kmax: maximum hyperedges size in H
    Ω: group interation function (e.g., planted partition)
    """

    kmin, kmax = minimum(keys(H.E)), maximum(keys(H.E))

    obj = 0
    for l = kmin:kmax
        El = H.E[l]
        for edge in keys(El)
            p = Ω.P(Z[edge])
#             a = Ω.aggregator(p)
            obj += El[edge]*log(Ω.ω(p, α))
        end
    end
    return obj
end

function modularity(H::hypergraph, Z::Array{<:Integer, 1}, Ω::IntensityFunction; α, bigInt::Bool=true)
    """
    Compute the modularity of a partition Z in a hypergraph H with interaction function Ω.
    H::hypergraph: the input hypergraph
    Z::array{Int64, 1}: array of group labels.
    Ω: group interaction function, as constructed by ΩFromDict(D)
    bigInt::Bool: whether to convert the degree sequence to an array of BigInt when evaluating the volume term in second_term_eval(). Recommended.
    return: Q::Float, the modularity term in the HSBM likelihood for H, Z, and Ω.
    """

    cut = first_term_eval(H, Z, Ω; α=α)
    vol = second_term_eval(H, Z, Ω; α=α, bigInt = bigInt)

    return cut - vol
end

function logLikelihood(H::hypergraph, Z::Array{<:Integer, 1}, Ω::IntensityFunction; α, bigInt::Bool=true)
    """
    Compute the HSBM log-likelihood of a partition Z in a hypergraph H with interaction function Ω.
    H::hypergraph: the input hypergraph
    Z::array{Int64, 1}: array of group labels.
    Ω: group interaction function, as constructed by ΩFromDict(D)
    bigInt::Bool: whether to convert the degree sequence to an array of BigInt when evaluating the volume term in second_term_eval(). Recommended.
    return: L::Float, the log-likelihood of H with parameters Z and Ω under the HSBM model.
    """
    Q = modularity(H, Z, Ω; α=α, bigInt=bigInt)

    D = H.D

    K, C = 0, 0

    logD = log.(D)

    kmax = maximum(keys(H.E))

    for ℓ = 1:kmax
        if haskey(H.E, ℓ)
            El = H.E[ℓ]
            for edge in keys(El)
                c = counting_coefficient(edge)
                weight = El[edge]
                K += weight*sum(logD[edge])
                C += weight*log(c) - log(factorial(weight)) # can maybe improve on performance here?
            end
        end
    end
    return Q, K, C
end

max_stars_repo_name: PhilChodrow/HypergraphModularity

max_issues_repo_name: PhilChodrow/HypergraphModularity

max_forks_repo_name: PhilChodrow/HypergraphModularity


text:
struct GeneralJONSWAP{K} <: WhittleLikelihoodInference.UnknownAcvTimeSeriesModel{1,Float64}
    α::Float64
    ωₚ::Float64
    γ::Float64
    r::Float64
    s::Float64

    r_over_s::Float64
    ωₚ²::Float64
    ωₚ³::Float64
    ωₚˢ::Float64
    ωₚˢ⁻¹::Float64
    logγ::Float64
    s⁻¹::Float64
    logωₚminuss⁻¹::Float64
    function GeneralJONSWAP{K}(α,ωₚ,γ,r,s) where {K}
        α > 0 || throw(ArgumentError("GeneralJONSWAP requires α > 0"))
        ωₚ > 0 || throw(ArgumentError("GeneralJONSWAP requires ωₚ > 0"))
        γ >= 1 || throw(ArgumentError("GeneralJONSWAP requires γ > 1"))
        r > 1 || throw(ArgumentError("GeneralJONSWAP requires r > 1"))
        s > 0 || throw(ArgumentError("GeneralJONSWAP requires s > 0"))
        new{K}(α,ωₚ,γ,r,s, r/s,ωₚ^2,ωₚ^3,ωₚ^s,ωₚ^(s-1),log(γ),inv(s),log(ωₚ)-inv(s))
    end
    function GeneralJONSWAP{K}(x::AbstractVector{Float64}) where {K}
        @boundscheck checkparameterlength(x,GeneralJONSWAP{K})
        @inbounds GeneralJONSWAP{K}(x[1], x[2], x[3], x[4], x[5])
    end
end

# functions to throw informative error if type parameter not provided
GeneralJONSWAP(x::AbstractVector{Float64}) = GeneralJONSWAP(1,1,1,1,1)
GeneralJONSWAP(α,ωₚ,γ,r,s) = error("GeneralJONSWAP process requires the ammount of aliasing specified as a type parameter. Use GeneralJONSWAP{K}() where K ∈ N.")

WhittleLikelihoodInference.npars(::Type{GeneralJONSWAP{K}}) where {K} = 5
WhittleLikelihoodInference.nalias(::GeneralJONSWAP{K}) where {K} = K

WhittleLikelihoodInference.lowerbounds(::Type{GeneralJONSWAP{K}}) where {K} = [0,0,1,1,0]
WhittleLikelihoodInference.upperbounds(::Type{GeneralJONSWAP{K}}) where {K} = [Inf,Inf,Inf,Inf,Inf]

@inline @fastmath function WhittleLikelihoodInference.sdf(model::GeneralJONSWAP{K}, ω::Real) where {K}
    α,ωₚ,γ,r,s = model.α,model.ωₚ,model.γ,model.r,model.s
    ω = abs(ω)
    if ω > 1e-10
        σ1² = 0.0049 + 0.0032 * (ω > ωₚ)
        δ = exp(-1 / (2σ1²) * (ω / ωₚ - 1)^2)
        ωₚˢ_over_ωˢ = model.ωₚˢ / (ω^s)
        return (α * ω^(-r) * exp(-(model.r_over_s) * ωₚˢ_over_ωˢ) * γ^δ) / 2
    else
        return 0.0
    end
end

@propagate_inbounds @fastmath function WhittleLikelihoodInference.grad_add_sdf!(out, model::GeneralJONSWAP{K}, ω::Real) where {K}
    @boundscheck checkbounds(out,1:5)
    @inbounds begin
        ω = abs(ω)
        if ω > 1e-10
            α,ωₚ,γ,r,s = model.α,model.ωₚ,model.γ,model.r,model.s
            σ1² = 0.0049 + 0.0032 * (ω > ωₚ)
            δ = exp(-1 / (2σ1²) * (ω / ωₚ - 1)^2)
            ω⁻ˢ = ω^(-s)
            ωₚˢ_over_ωˢ = model.ωₚˢ * ω⁻ˢ
            ∂S∂α = ω^(-r) * exp(-(model.r_over_s) * ωₚˢ_over_ωˢ) * γ^δ / 2
            sdf = α * ∂S∂α

            ∂S∂ωₚ = sdf * (δ*model.logγ * ω / σ1² *(ω-ωₚ) / model.ωₚ³ - r*model.ωₚˢ⁻¹*ω⁻ˢ)
            ∂S∂γ = sdf * δ / γ
            ∂S∂r = sdf * (-log(ω)-ωₚˢ_over_ωˢ*model.s⁻¹)
            ∂S∂s = sdf * model.r_over_s * ωₚˢ_over_ωˢ * (log(ω)-model.logωₚminuss⁻¹)

            out[1] += ∂S∂α
            out[2] += ∂S∂ωₚ
            out[3] += ∂S∂γ
            out[4] += ∂S∂r
            out[5] += ∂S∂s
        end # 0 otherwise
    end
    return nothing
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Center()

The transform that removes the mean of the variables.
"""
struct Center <: Colwise end

assertions(::Type{Center}) = [assert_continuous]

isrevertible(::Type{Center}) = true

colcache(::Center, x) = mean(x)

colapply(::Center, x, μ)  = @. x - μ

colrevert(::Center, y, μ) = @. y + μ


text:
using ExoplanetsSysSim

#=
    p_lim_arr_num = [0.5, 1.25, 2.5, 5., 10., 20., 40., 80., 160., 320.]
    r_lim_arr_num = [1., 1.25, 1.5, 1.75, 2.]
    rate_tab_init = [0.027  0.149 0.354 1.37  2.56  1.96 2.85  4.77 0.0001 ;
                     0.02   0.063 0.432 0.9   1.72  1.46 2.4   5.42 3.63 ;
                     0.0074 0.032 0.389 0.545 1.32  1.65 1.96  1.69 1.17 ;
                     0.0001 0.017 0.156 0.542 0.777 1.35 0.965 1.22 2.49]*0.01
=#

function set_test_param(sim_param_closure::SimParam)
    add_param_fixed(sim_param_closure,"num_targets_sim_pass_one",50000)  # For faster simulated catalogs

    p_lim_arr_num = [10., 20.]
    r_lim_arr_num = [1.0, 1.25]
    rate_tab_init = reshape([2.5]*0.01, (1,1))
    add_param_fixed(sim_param_closure, "p_lim_arr", p_lim_arr_num)
    add_param_fixed(sim_param_closure, "r_lim_arr", r_lim_arr_num*ExoplanetsSysSim.earth_radius)
    add_param_active(sim_param_closure, "obs_par", rate_tab_init)
    return sim_param_closure
end


text:
"""
	get_bin_content(inFile, h2dName, Emin, Emax)

Description of ```get_bin_content```
------------------------------
Provides bin content in the specified energy range.

Arguments:
	+ inFile 	- The ROOTFile of the 2d histogram
	+ h2dName 	- The name of the Th2F. Provided as string. (e.g. "h")
	+ Emin 		- Minimum energy of the deisred E range
	+ Emax      - Maximum energy of the deisred E range
-------------------------------
Example use:

We have TH2F file called "hist2D.root" opened in UnROOT as f = ROOTFile("hist2D.root"). The <NAME> of the histogram (generated in root as TH2F("<NAME>", "TITLE", ...)) is "h2".

To find bin content n at E ∈ (500, 600): n = ```get_bin_content(f, "h2", 500, 600)```.

FUNCTION ONLY WORKS IF THE HISTOGRAM IS SQUARE - meaning it has the same number of bins on x and y axes.
"""
function get_bin_content(inFile, h2dName, Emin, Emax)
	nbins = inFile[h2dName][:fXaxis_fNbins]
	binsize = inFile[h2dName][:fXaxis_fXmax] / nbins

	row = floor(Emax/binsize) - 1
	col = floor(Emin/binsize) + 1

	#=
	To get the index from the array :fN the algorithm is as follows.
	First 3 bins are 0 (for some unknown reason, that's how root is written)
	The next <nbins> number of the array are also zero - these correspond to all the underwflow xbins
	Afterwards <nbins> number of bins correspond to the first row in the histogram.
	2 bins are empty again - these are overflow bins for x and y (in one row).
	=#
	idx = convert(Int64, 3 + nbins + (nbins+2)*row + col)

	bincontent = inFile[h2dName][:fN][idx]
	return convert(Float64, bincontent)
end

"""
	`MP_heatmap( xTicks::Vector{<:Real}, yTicks::Vector{<:Real}, data::Vector{<:Real}, stepSize = 100 )``

Description of `MP_heatmap`
------------------------------
Returns a Plots.heatmap object when provided 3 arrays x,y,z. Where x,y are the axes and z is the color of the cell.

Arguments:

	xTicks 	- Vector x values associated with z values. Works with dataframe columns. Can be provided as eg. "df.EMins".
	yTicks 	- Vector y values associated with z values. Works with dataframe columns. Can be provided as eg. "df.EMaxs".
	data 		- Values for the z axis / color. Works with dataframe columns. Can be provided as eg. "df.NPassed".

-------------------------------
Example use:
hm = MP_heatmap( df.EMins, df.EMaxs, df.NPassed)

"""
function MP_heatmap( xTicks::Vector{<:Real}, yTicks::Vector{<:Real}, data::Vector{<:Real}, stepSize = 100; kwargs... )
    dataMatrix = zeros(
                        length(unique(yTicks)),
                        length(unique(xTicks))
                        )
    if ( unique(xTicks)[1] <= 0 )
        x  = unique(xTicks) .+ stepSize/2   # x specifies the midpoints of the bins
    else
        x  = unique(xTicks) .- stepSize/2   # x specifies the midpoints of the bins
    end

    if ( unique(yTicks)[1] <= 0 )
        y  = unique(yTicks) .+ stepSize/2   # y specifies the midpoints of the bins
    else
        y  = unique(yTicks) .- stepSize/2   # y specifies the midpoints of the bins
    end

    for d in 1:length(data)  # each data point is mapped onto the matrix, where row and col are specified
        if ( length(data) != length(xTicks) || length(data) != length(yTicks))
            error("Arrays must be the same size!")
        end

        if (yTicks[1] <= 0) # since Julia indexing starts from 1, we must add 1 if row was to be 0
            r = convert(Int, yTicks[d]/stepSize + 1) # gives row index
        else
            r = convert(Int, yTicks[d]/stepSize ) # gives row index
        end

        if (xTicks[1] <= 0) # since Julia indexing starts from 1, we must add 1 if col was to be 0
            c = convert(Int, xTicks[d]/stepSize + 1) # gives row index
        else
            c = convert(Int, xTicks[d]/stepSize ) # gives row index
        end

        dataMatrix[r,c] = data[d]
    end
    dataMatrix = replace!(dataMatrix, 0.0 => NaN)

    hm = Plots.heatmap(x,y,dataMatrix; kwargs...)

    return hm
end

"""
Function borrowed and slightly changed from UnROOT.jl

Returns 2d Heatmap when TH1D/TH2D object is passed.

Example use:
f =  ROOTFile("82Se_2ubb_SpectrumG0.root")
th = f["h2"] ; where h2 is the name of the TH2D object from root.
MP_heatmap(th).
"""
function MP_heatmap(th::Dict{Symbol, Any}; kwargs...)
    xmin = th[:fXaxis_fXmin]
    xmax = th[:fXaxis_fXmax]
    xnbins = th[:fXaxis_fNbins]
    xbins = isempty(th[:fXaxis_fXbins]) ? range(xmin, xmax, length=xnbins+1) : th[:fXaxis_fXbins];
    counts = th[:fN]

    if th[:fYaxis_fNbins] > 1
        ymin = th[:fYaxis_fXmin]
        ymax = th[:fYaxis_fXmax]
        ynbins = th[:fYaxis_fNbins]

        ybins = isempty(th[:fYaxis_fXbins]) ? range(ymin, ymax, length=ynbins+1) : th[:fYaxis_fXbins];

        counts = reshape(counts, (xnbins+2, ynbins+2))[2:end-1, 2:end-1]   # matrix with dimensions (xnbins, ynbins) and values of the individual bin counts
        edges = (xbins, ybins)
    else
        counts = counts[2:end-1]
        edges = (xbins,)
    end
    return Plots.heatmap(StatsBase.Histogram(edges, counts), nbins = (xbins, ybins), ; kwargs...)
end


text:
using MyPkgDemo
using Test

@testset "test for add_nums" begin
 # Write your tests here.
   @test add_nums(1, 3) == 4
   @test add_nums(4, 3) == 7
   @test add_nums(7, 3) == 10
   @test add_nums(8, 3) == 11
   @test add_nums(8, 3, 1) == 12
   @test add_nums(8, 3, 8) == 19
   @test add_nums(1, 4, 3) == 8
end

max_stars_repo_name: junedwx/MyPkgDemo

max_issues_repo_name: junedwx/MyPkgDemo

max_forks_repo_name: junedwx/MyPkgDemo


text:
# This file is a part of Julia, but is derived from
# https://github.com/floitsch/double-conversion which has the following license
#
# Copyright 2006-2014, the V8 project authors. All rights reserved.
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are
# met:
#
#     * Redistributions of source code must retain the above copyright
#       notice, this list of conditions and the following disclaimer.
#     * Redistributions in binary form must reproduce the above
#       copyright notice, this list of conditions and the following
#       disclaimer in the documentation and/or other materials provided
#       with the distribution.
#     * Neither the name of Google Inc. nor the names of its
#       contributors may be used to endorse or promote products derived
#       from this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
# A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
# OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
# LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
# THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

module Bignums

import Base: ==, <

export Bignum

const kMaxSignificantBits = 3584

const Chunk = UInt32
const DoubleChunk = UInt64

const kChunkSize = sizeof(Chunk) * 8
const kDoubleChunkSize = sizeof(DoubleChunk) * 8
# With bigit size of 28 we loose some bits, but a double still fits easily
# into two chunks, and more importantly we can use the Comba multiplication.
const kBigitSize = 28
const kBigitMask = Chunk((1 << kBigitSize) - 1)
# Every instance allocates kBigitLength chunks on the stack. Bignums cannot
# grow. There are no checks if the stack-allocated space is sufficient.
const kBigitCapacity = div(kMaxSignificantBits,kBigitSize)

mutable struct Bignum
    bigits::Array{UInt32,1}
    used_digits::Int32
    exponent::Int32
    function Bignum()
        bigits = Array{UInt32}(kBigitCapacity)
        @inbounds for i = 1:kBigitCapacity
            bigits[i] = 0
        end
        new(bigits,0,0)
    end
end

==(a::Bignum,b::Bignum) = compare(a,b) == 0
<(a::Bignum,b::Bignum) = compare(a,b) < 0

times10!(x::Bignum) = multiplybyuint32!(x,UInt32(10))

plusequal(a,b,c) = pluscompare(a,b,c) == 0
pluslessequal(a,b,c) = pluscompare(a,b,c) <= 0
plusless(a,b,c) = pluscompare(a,b,c) < 0
lessequal(a::Bignum,b::Bignum) = compare(a,b) <= 0
less(a::Bignum,b::Bignum) = compare(a,b) < 0

bigitlength(x::Bignum) = x.used_digits + x.exponent

bitsize(value) = 8 * sizeof(value)

function zero!(x::Bignum)
    for i = 1:x.used_digits
        @inbounds x.bigits[i] = 0
    end
    x.used_digits = 0
    x.exponent = 0
    return
end

function clamp!(x::Bignum)
    @inbounds while (x.used_digits > 0 && x.bigits[x.used_digits] == 0)
        x.used_digits -= 1
    end
    x.used_digits == 0 && (x.exponent = 0)
    return
end

isclamped(x::Bignum) = x.used_digits == 0 || x.bigits[x.used_digits] != 0

function align!(x::Bignum,other::Bignum)
    @inbounds if x.exponent > other.exponent
        zero_digits = x.exponent - other.exponent
        for i = x.used_digits:-1:1
            x.bigits[i + zero_digits] = x.bigits[i]
        end
        for i = 1:zero_digits
            x.bigits[i] = 0
        end
        x.used_digits += zero_digits
        x.exponent -= zero_digits
    end
    return
end

function bigitshiftleft!(x::Bignum,shift_amount)
    carry::UInt32 = 0
    @inbounds begin
    for i = 1:x.used_digits
        new_carry::Chunk = x.bigits[i] >> (kBigitSize - shift_amount)
        x.bigits[i] = ((x.bigits[i] << shift_amount) + carry) & kBigitMask
        carry = new_carry
    end
    if carry != 0
        x.bigits[x.used_digits+1] = carry
        x.used_digits += 1
    end
    end
    return
end

function subtracttimes!(x::Bignum,other::Bignum,factor)
    if factor < 3
        for i = 1:factor
            subtractbignum!(x,other)
        end
        return
    end
    borrow::Chunk = 0
    exponent_diff = other.exponent - x.exponent
    @inbounds begin
    for i = 1:other.used_digits
        product::DoubleChunk = DoubleChunk(factor) * other.bigits[i]
        remove::DoubleChunk = borrow + product
        difference::Chunk = (x.bigits[i+exponent_diff] - (remove & kBigitMask)) % Chunk
        x.bigits[i+exponent_diff] = difference & kBigitMask
        borrow = ((difference >> (kChunkSize - 1)) + (remove >> kBigitSize)) % Chunk
    end
    for i = (other.used_digits + exponent_diff + 1):x.used_digits
        borrow == 0 && return
        difference::Chunk = x.bigits[i] - borrow
        x.bigits[i] = difference & kBigitMask
        borrow = difference >> (kChunkSize - 1)
    end
    end
    clamp!(x)
end

function assignuint16!(x::Bignum,value::UInt16)
    zero!(x)
    value == 0 && return
    x.bigits[1] = value
    x.used_digits = 1
    return
end

const kUInt64Size = 64
function assignuint64!(x::Bignum,value::UInt64)
    zero!(x)
    value == 0 && return
    needed_bigits = div(kUInt64Size,kBigitSize) + 1
    @inbounds for i = 1:needed_bigits
        x.bigits[i] = value & kBigitMask
        value >>= kBigitSize
    end
    x.used_digits = needed_bigits
    clamp!(x)
end

function assignbignum!(x::Bignum,other::Bignum)
    x.exponent = other.exponent
    @inbounds begin
    for i = 1:other.used_digits
        x.bigits[i] = other.bigits[i]
    end
    for i = (other.used_digits+1):x.used_digits
        x.bigits[i] = 0
    end
    end
    x.used_digits = other.used_digits
    return
end

function adduint64!(x::Bignum,operand::UInt64)
    operand == 0 && return
    other = Bignum()
    assignuint64!(other,operand)
    addbignum!(x,other)
end

function addbignum!(x::Bignum,other::Bignum)
    align!(x,other)
    carry::Chunk = 0
    bigit_pos = other.exponent - x.exponent
    @inbounds for i = 1:other.used_digits
        sum::Chunk = x.bigits[bigit_pos+1] + other.bigits[i] + carry
        x.bigits[bigit_pos+1] = sum & kBigitMask
        carry = sum >> kBigitSize
        bigit_pos += 1
    end
    @inbounds while carry != 0
        sum = x.bigits[bigit_pos+1] + carry
        x.bigits[bigit_pos+1] = sum & kBigitMask
        carry = sum >> kBigitSize
        bigit_pos += 1
    end
    x.used_digits = max(bigit_pos,x.used_digits)
    return
end

function subtractbignum!(x::Bignum,other::Bignum)
    align!(x,other)
    offset = other.exponent - x.exponent
    borrow = Chunk(0)
    @inbounds begin
    for i = 1:other.used_digits
        difference = x.bigits[i+offset] - other.bigits[i] - borrow
        x.bigits[i+offset] = difference & kBigitMask
        borrow = difference >> (kChunkSize - 1)
    end
    i = other.used_digits+1
    while borrow != 0
        difference = x.bigits[i+offset] - borrow
        x.bigits[i+offset] = difference & kBigitMask
        borrow = difference >> (kChunkSize - 1)
        i += 1
    end
    end
    clamp!(x)
end

function shiftleft!(x::Bignum,shift_amount)
    x.used_digits == 0 && return
    x.exponent += div(shift_amount,kBigitSize)
    local_shift = shift_amount % kBigitSize
    bigitshiftleft!(x,local_shift)
end

function multiplybyuint32!(x::Bignum,factor::UInt32)
    factor == 1 && return
    if factor == 0
        zero!(x)
        return
    end
    x.used_digits == 0 && return
    carry::DoubleChunk = 0
    @inbounds begin
    for i = 1:x.used_digits
        product::DoubleChunk = (factor % DoubleChunk) * x.bigits[i] + carry
        x.bigits[i] = (product & kBigitMask) % Chunk
        carry = product >> kBigitSize
    end
    while carry != 0
        x.bigits[x.used_digits+1] = carry & kBigitMask
        x.used_digits += 1
        carry >>= kBigitSize
    end
    end
    return
end

function multiplybyuint64!(x::Bignum,factor::UInt64)
    factor == 1 && return
    if factor == 0
        zero!(x)
        return
    end
    carry::UInt64 = 0
    low::UInt64 = factor & 0xFFFFFFFF
    high::UInt64 = factor >> 32
    @inbounds begin
    for i = 1:x.used_digits
        product_low::UInt64 = low * x.bigits[i]
        product_high::UInt64 = high * x.bigits[i]
        tmp::UInt64 = (carry & kBigitMask) + product_low
        x.bigits[i] = tmp & kBigitMask
        carry = (carry >> kBigitSize) + (tmp >> kBigitSize) +
                (product_high << (32 - kBigitSize))
    end
    while carry != 0
        x.bigits[x.used_digits+1] = carry & kBigitMask
        x.used_digits += 1
        carry >>= kBigitSize
    end
    end
    return
end

const kFive27 = UInt64(0x6765c793fa10079d)
const kFive1 = UInt16(5)
const kFive2 = UInt16(kFive1 * 5)
const kFive3 = UInt16(kFive2 * 5)
const kFive4 = UInt16(kFive3 * 5)
const kFive5 = UInt16(kFive4 * 5)
const kFive6 = UInt16(kFive5 * 5)
const kFive7 = UInt32(kFive6 * 5)
const kFive8 = UInt32(kFive7 * 5)
const kFive9 = UInt32(kFive8 * 5)
const kFive10 = UInt32(kFive9 * 5)
const kFive11 = UInt32(kFive10 * 5)
const kFive12 = UInt32(kFive11 * 5)
const kFive13 = UInt32(kFive12 * 5)
const kFive1_to_12 = UInt32[kFive1, kFive2, kFive3, kFive4, kFive5, kFive6,
        kFive7, kFive8, kFive9, kFive10, kFive11, kFive12]
function multiplybypoweroften!(x::Bignum,exponent)
    exponent == 0 && return
    x.used_digits == 0 && return
    remaining_exponent = exponent
    while remaining_exponent >= 27
        multiplybyuint64!(x,kFive27)
        remaining_exponent -= 27
    end
    while remaining_exponent >= 13
        multiplybyuint32!(x,kFive13)
        remaining_exponent -= 13
    end
    remaining_exponent > 0 && multiplybyuint32!(x,
                            kFive1_to_12[remaining_exponent])
    shiftleft!(x,exponent)
end

function square!(x::Bignum)
    product_length = 2 * x.used_digits
    (1 << (2 * (kChunkSize - kBigitSize))) <= x.used_digits && error("unimplemented")
    accumulator::DoubleChunk = 0
    copy_offset = x.used_digits
    @inbounds begin
    for i = 1:x.used_digits
        x.bigits[copy_offset + i] = x.bigits[i]
    end
    for i = 1:x.used_digits
        bigit_index1 = i-1
        bigit_index2 = 0
        while bigit_index1 >= 0
            chunk1::Chunk = x.bigits[copy_offset + bigit_index1 + 1]
            chunk2::Chunk = x.bigits[copy_offset + bigit_index2 + 1]
            accumulator += (chunk1 % DoubleChunk) * chunk2
            bigit_index1 -= 1
            bigit_index2 += 1
        end
        x.bigits[i] = (accumulator % Chunk) & kBigitMask
        accumulator >>= kBigitSize
    end
    for i = x.used_digits+1:product_length
        bigit_index1 = x.used_digits - 1
        bigit_index2 = i - bigit_index1 - 1
        while bigit_index2 < x.used_digits
            chunk1::Chunk = x.bigits[copy_offset + bigit_index1 + 1]
            chunk2::Chunk = x.bigits[copy_offset + bigit_index2 + 1]
            accumulator += (chunk1 % DoubleChunk) * chunk2
            bigit_index1 -= 1
            bigit_index2 += 1
        end
        x.bigits[i] = (accumulator % Chunk) & kBigitMask
        accumulator >>= kBigitSize
    end
    end
    x.used_digits = product_length
    x.exponent *= 2
    clamp!(x)
end

function assignpoweruint16!(x::Bignum,base::UInt16,power_exponent::Int)
    if power_exponent == 0
        assignuint16!(x,UInt16(1))
        return
    end
    zero!(x)
    shifts::Int = 0
    while base & UInt16(1) == UInt16(0)
        base >>= UInt16(1)
        shifts += 1
    end
    bit_size::Int = 0
    tmp_base::Int= base
    while tmp_base != 0
        tmp_base >>= 1
        bit_size += 1
    end
    final_size = bit_size * power_exponent
    mask::Int = 1
    while power_exponent >= mask
        mask <<= 1
    end
    mask >>= 2
    this_value::UInt64 = base
    delayed_multiplication = false
    max_32bits::UInt64 = 0xFFFFFFFF
    while mask != 0 && this_value <= max_32bits
        this_value *= this_value
        if (power_exponent & mask) != 0
            base_bits_mask::UInt64 = ~(UInt64(1) << (64 - bit_size) - 1)
            high_bits_zero = (this_value & base_bits_mask) == 0
            if high_bits_zero
                this_value *= base
            else
                delayed_multiplication = true
            end
        end
        mask >>= 1
    end
    assignuint64!(x,this_value)
    delayed_multiplication && multiplybyuint32!(x,UInt32(base))
    while mask != 0
        square!(x)
        (power_exponent & mask) != 0 && multiplybyuint32!(x,UInt32(base))
        mask >>= 1
    end
    shiftleft!(x,shifts * power_exponent)
end

function dividemodulointbignum!(x::Bignum,other::Bignum)
    bigitlength(x) < bigitlength(other) && return UInt16(0)
    align!(x,other)
    result::UInt16 = 0
    @inbounds begin
    while bigitlength(x) > bigitlength(other)
        result += x.bigits[x.used_digits] % UInt16
        subtracttimes!(x,other,x.bigits[x.used_digits])
    end
    this_bigit::Chunk = x.bigits[x.used_digits]
    other_bigit::Chunk = other.bigits[other.used_digits]
    if other.used_digits == 1
        quotient = reinterpret(Int32,div(this_bigit,other_bigit))
        x.bigits[x.used_digits] = this_bigit - other_bigit * reinterpret(UInt32,quotient)
        result += quotient % UInt16
        clamp!(x)
        return result
    end
    end
    division_estimate = reinterpret(Int32,div(this_bigit,other_bigit+Chunk(1)))
    result += division_estimate % UInt16
    subtracttimes!(x,other,division_estimate)
    other_bigit * (division_estimate+1) > this_bigit && return result
    while lessequal(other, x)
        subtractbignum!(x,other)
        result += UInt16(1)
    end
    return result
end

function pluscompare(a::Bignum,b::Bignum,c::Bignum)
    bigitlength(a) < bigitlength(b) && return pluscompare(b,a,c)
    bigitlength(a) + 1 < bigitlength(c) && return -1
    bigitlength(a) > bigitlength(c) && return 1
    a.exponent >= bigitlength(b) && bigitlength(a) < bigitlength(c) && return -1
    borrow::Chunk = 0
    min_exponent = min(a.exponent,b.exponent,c.exponent)
    for i = (bigitlength(c)-1):-1:min_exponent
        chunk_a::Chunk = bigitat(a,i)
        chunk_b::Chunk = bigitat(b,i)
        chunk_c::Chunk = bigitat(c,i)
        sum::Chunk = chunk_a + chunk_b
        if sum > chunk_c + borrow
            return 1
        else
            borrow = chunk_c + borrow - sum
            borrow > 1 && return -1
            borrow <<= kBigitSize
        end
    end
    borrow == 0 && return 0
    return -1
end

function compare(a::Bignum,b::Bignum)
    bigit_length_a = bigitlength(a)
    bigit_length_b = bigitlength(b)
    bigit_length_a < bigit_length_b && return -1
    bigit_length_a > bigit_length_b && return 1
    for i = (bigit_length_a-1):-1:min(a.exponent,b.exponent)
        bigit_a::Chunk = bigitat(a,i)
        bigit_b::Chunk = bigitat(b,i)
        bigit_a < bigit_b && return -1
        bigit_a > bigit_b && return 1
    end
    return 0
end

function bigitat(x::Bignum,index)
    index >= bigitlength(x) && return Chunk(0)
    index < x.exponent && return Chunk(0)
    @inbounds ret = x.bigits[index - x.exponent+1]::Chunk
    return ret
end

end # module


text:
# This file was generated, do not modify it. # hide
df = DataFrame(a=1:3, b=4:6)
combine(df, :a => sum, nrow)


text:
#####
##### `AbstractTangent`
#####

"""
The subtypes of `AbstractTangent` define a custom \"algebra\" for chain
rule evaluation that attempts to factor various features like complex derivative
support, broadcast fusion, zero-elision, etc. into nicely separated parts.

In general a differential type is the type of a derivative of a value.
The type of the value is for contrast called the primal type.
Differential types correspond to primal types, although the relation is not one-to-one.
Subtypes of  `AbstractTangent` are not the only differential types.
In fact for the most common primal types, such as `Real` or `AbstractArray{Real}` the
the differential type is the same as the primal type.

In a circular definition: the most important property of a differential is that it should
be able to be added (by defining `+`) to another differential of the same primal type.
That allows for gradients to be accumulated.

It generally also should be able to be added to a primal to give back another primal, as
this facilitates gradient descent.

All subtypes of `AbstractTangent` implement the following operations:

 - `+(a, b)`: linearly combine differential `a` and differential `b`
 - `*(a, b)`: multiply the differential `b` by the scaling factor `a`
 -  `Base.zero(x) = ZeroTangent()`: a zero.

Further, they often implement other linear operators, such as `conj`, `adjoint`, `dot`.
Pullbacks/pushforwards are linear operators, and their inputs are often
`AbstractTangent` subtypes.
Pullbacks/pushforwards in-turn call other linear operators on those inputs.
Thus it is desirable to have all common linear operators work on `AbstractTangent`s.
"""
abstract type AbstractTangent end

Base.:+(x::AbstractTangent) = x

@inline Base.conj(x::AbstractTangent) = x


text:
####################################################
# Script to train a deep network with feedback alignment (SGD)
# patch connectivity is controled via "patch_size" variable

###################################################
# Include libraries and core

using Pkg; Pkg.activate("./../core/"); Pkg.instantiate()
include("./../core/src/ratenets/core.jl")

###################################################
# Define task and network

dataset = "CIFAR10" # dataset to be used: "MNIST" or "CIFAR10"
colorimage = dataset == "CIFAR10"
hidden_size = 1000 # number of neurons in hidden layer (0 for simple perceptron)
patch_size = 10 # (linear) patch size of receptive field in first layer p ∈ [1,28] or [1,32] for MNIST/CIFAR10 resp.
n_inits = 3 # number of initializations per net (for averaging)
iterations = 10^6 # iterations per initialization
learningrates = 5e-3*[1; 1] # learningrates per layer
nonlinearity = [relu!, relu!] # nonlinearities
nonlinearity_diff = [relu_diff!, relu_diff!]

###################################################
# Load and preprocess data

smallimgs, labels, smallimgstest, labelstest = import_data(dataset)
subtractmean!(smallimgs) # substract pixel-wise mean
subtractmean!(smallimgstest)

###################################################
# Train network

print(string("train network with 1 hidden layer with ",hidden_size," neurons...\n"))
for j in 1:n_inits
    net = Network_RF_sparseconnect([size(smallimgs)[1], hidden_size, 10])
    set_connectivity_matrix!(net,patch_size; mode = "randn", colorimage = colorimage) # set connectivity of patches
    @time loss = learn!(net, getsmallimg, getlabel, iterations, learningrates, 0.,
                          nonlinearity = nonlinearity, nonlinearity_diff = nonlinearity_diff)
    error_train = geterrors(net, smallimgs, labels, nonlinearity = nonlinearity)
    error_test = geterrors(net, smallimgstest, labelstest, nonlinearity = nonlinearity)
    print(string(100 * error_train, " % on training set\n"))
    print(string(100 * error_test, " % on test set\n"))
end


text:
using .Glimpse
const Gl = Glimpse

using ..Structures
using ..DFControl

function Gl.Diorama(str::Structure)
    dio = Diorama(; name = Symbol("DFControl Structure"))
    add_structure!(dio, str)
    dio[Gl.PointLight].data[1] = Gl.PointLight(; diffuse = 1.5f0, ambient = 1.5f0)
    Gl.center_cameras(dio)
    return dio
end
function add_structure!(dio::Diorama, str::Structure; polyhedra::Vector{Symbol} = Symbol[])
    Entity(dio,
           Gl.assemble_wire_axis_box(; position = zero(Point3f0), x = Structures.ustrip(Structures.a(str)),
                                     y = Structures.ustrip(Structures.b(str)), z = Structures.ustrip(Structures.c(str)),
                                     color = Gl.BLACK)...)
    for at in str.atoms
        origin = Point3f0(Structures.ustrip(at.position_cart)...)
        color = RGB{Float32}(at.element.color...)
        Entity(dio,
               Gl.assemble_sphere(origin; color = color,
                                  radius = Float32(0.2 + at.element.atomic_weight / 400.0f0))...,
               Gl.Text(; str = String(at.name), font_size = 2))
        if sum(at.magnetization) != 0
            Entity(dio,
                   Gl.assemble_arrow(origin, origin + Point3f0(at.magnetization...);
                                     color = color, radius_ratio = 1.5f0,
                                     thickness = 0.05f0)...)
        end
    end
    for el in map(x->x.element, polyhedra)
        for at in str[el]
            distances = sort(distance.((at,), str.atoms))
            count = length(filter(i -> distances[i+1] - distances[i] < 1e-1angstrom),
                           1:length(distances)-1)
            poly = polyhedron(str, count)
        end
    end
    return dio
end


text:
const FIG_DIR = joinpath(WATERFALL_DIR,"fig")

# Dimensions
WIDTH = 800
HEIGHT = 400
BORDER = 10
TOP_BORDER = 20
BOTTOM_BORDER = 30
LEFT_BORDER = 70
RIGHT_BORDER = 20
SEP = 10
FONTSIZE = 18

VMAX = 22.5
VMIN = 16

# Names
const SAMPLE_COL = :Sample
const VALUE_COL = :Value

# Colors
const HEX_LOSS = parse(Luxor.Colorant,"#ca0020")
const HEX_GAIN = parse(Luxor.Colorant,"#0571b0")

# Default input editing
const SEED = 1234
const DEFAULT_FUZZINESS = (0.01,0.3)
const DEFAULT_DISTRIBUTION = :normal      # linear, normal, triangular
const DEFAULT_NSAMPLE = 5
const DEFAULT_NCOR = true

const LEADING = 1.25

# COLORCYCLE = Dict(
#     # https://colorbrewer2.org/?type=qualitative&scheme=Paired&n=12
#     :bright => ["#a6cee3","#1f78b4","#b2df8a","#33a02c","#fb9a99","#e31a1c","#fdbf6f","#ff7f00","#cab2d6","#6a3d9a","#ffff99","#b15928"],
#     # https://colorbrewer2.org/?type=qualitative&scheme=Set3&n=12
#     :pastel => ["#8dd3c7","#ffffb3","#bebada","#fb8072","#80b1d3","#fdb462","#b3de69","#fccde5","#d9d9d9","#bc80bd","#ccebc5","#ffed6f"],
#     # https://arxiv.org/abs/2107.02270
#     # https://github.com/mpetroff/accessible-color-cycles
#     6 => ["#5790fc", "#f89c20", "#e42536", "#964a8b", "#9c9ca1", "#7a21dd"],
#     8 => ["#1845fb", "#ff5e02", "#c91f16", "#c849a9", "#adad7d", "#86c8dd", "#578dff", "#656364"],
#     10 => ["#3f90da", "#ffa90e", "#bd1f01", "#94a4a2", "#832db6", "#a96b59", "#e76300", "#b9ac70", "#717581", "#92dadd"],
# )
# const COLORCYCLE = Dict(k => parse.(Luxor.Colorant, v) for (k,v) in COLORCYCLE)
COLORCYCLE = ["#3f90da", "#ffa90e", "#bd1f01", "#94a4a2", "#832db6", "#a96b59", "#e76300", "#b9ac70", "#717581", "#92dadd","#92dadd","#92dadd"]
# COLORCYCLE = ["#8dd3c7","#ffffb3","#bebada","#fb8072","#80b1d3","#fdb462","#b3de69","#fccde5","#d9d9d9","#bc80bd","#ccebc5","#ffed6f"]
# COLORCYCLE = ["#a6cee3","#1f78b4","#b2df8a","#33a02c","#fb9a99","#e31a1c","#fdbf6f","#ff7f00","#cab2d6","#6a3d9a","#ffff99","#b15928"]


text:
#hoursperyear(yr::Integer) = isleap(yr) ? 8784 : 8760
isleap(yr::Integer) = yr % 4 == 0 && (yr % 100 != 0 || yr % 400 == 0)
flip(x) = permutedims(x, (2,1))

showall(x) = show(stdout, "text/plain", x)

@userplot StackedArea
# Bugfix of StatsPlots.areaplot(), i.e. reverse the legend order.
@recipe function f(a::StackedArea)
    data = cumsum(a.args[end], dims=2)
    n = size(data, 2)
    x = length(a.args) == 1 ? (1:size(data, 1)) : a.args[1]
    seriestype := :line
    labels = haskey(plotattributes, :label) ? plotattributes[:label] : ["y$i" for i = 1:n]

    for i in n:-1:1
        @series begin
            fillrange := i > 1 ? data[:,i-1] : 0
            fillcolor --> i
            label := labels[i]
            x, data[:,i]
        end
    end
end

@userplot StackedArea_Polygons
# A simple "recipe" for Plots.jl to get stacked area plots. This version draws a line around each color patch.
# (Calling the other function with linecolor=:black doesn't draw the color around the sides or bottom.)
# usage: stackedarea(xvector, datamatrix, plotsoptions)
@recipe function f(a::StackedArea_Polygons)
    data = cumsum(a.args[end], dims=2)
    nx, nseries = size(data)
    x = length(a.args) == 1 ? (1:nx) : a.args[1]
    labels = haskey(plotattributes, :label) ? plotattributes[:label] : ["y$i" for i = 1:n]
    seriestype := :shape

    # create a filled polygon for each item
    for i = nseries:-1:1
        sx = vcat(x, reverse(x))
        sy = vcat(data[:, i], i==1 ? zeros(nx) : reverse(data[:, i-1]))
        @series begin
            fillcolor --> i
            label := labels[i]
            (sx, sy)
        end
    end
end

@userplot StackedBar
# Bugfix of StatsPlots.groupedbar(), i.e. reverse the legend order. Also default to bar_position=:stack.
@recipe function f(g::StackedBar; spacing = 0)
    x, y = StatsPlots.grouped_xy(g.args...)

    nr, nc = size(y)

    isstack = pop!(plotattributes, :bar_position, :stack) == :stack

    # extract xnums and set default bar width.
    # might need to set xticks as well
    xnums = if eltype(x) <: Number
        xdiff = length(x) > 1 ? mean(diff(x)) : 1
        bar_width --> 0.8 * xdiff
        x
    else
        bar_width --> 0.8
        ux = unique(x)
        xnums = (1:length(ux)) .- 0.5
        xticks --> (xnums, ux)
        xnums
    end
    @assert length(xnums) == nr

    # compute the x centers.  for dodge, make a matrix for each column
    x = if isstack
        x
    else
        bws = plotattributes[:bar_width] / nc
        bar_width := bws * clamp(1 - spacing, 0, 1)
        xmat = zeros(nr,nc)
        for r=1:nr
            bw = StatsPlots._cycle(bws, r)
            farleft = xnums[r] - 0.5 * (bw * nc)
            for c=1:nc
                xmat[r,c] = farleft + 0.5bw + (c-1)*bw
            end
        end
        xmat
    end

    # compute fillrange
    fillrange := if isstack
        y, fr = StatsPlots.groupedbar_fillrange(reverse(y, dims=2))
        fr
    else
        get(plotattributes, :fillrange, nothing)
    end

    seriestype := :bar

    if isstack
        label := haskey(plotattributes, :label) ? reverse(plotattributes[:label],dims=2) : ["y$i" for j = 1:1, i = nc:-1:1]
        fillcolor --> permutedims(nc:-1:1)
    end

    x, y
end

#=
function printtable(title::String, setnames::Vector, datarows::Array, columnname::Symbol = :none)
    println(title)
    push!(setnames, :value)
    df = DataFrame([fill(Symbol, length(setnames)-1); Float64], setnames, 0)
    for row in datarows
        push!(df, row)
    end
    if columnname != :none
        df = unstack(df, columnname, :value)
    end
    show(IOContext(STDOUT, displaysize=(100,120)), "text/plain", df)
end

printtable(jarr::JuMP.JuMPArray{Float64,2}) = printtable("", jarr)
printtable(title::String, jarr::JuMP.JuMPArray{Float64,2}) =
    println("$title\n", DataFrame([jarr.indexsets[1] jarr.innerArray], [:_; jarr.indexsets[2]]))
=#

max_stars_repo_name: niclasmattsson/Supergrid

max_issues_repo_name: niclasmattsson/Supergrid

max_forks_repo_name: niclasmattsson/Supergrid


text:
using Test
using MotifSequenceGenerator

@testset "possiblesums" begin
    a = [1,1,3,4]
    mina = minimum(a); maxa = maximum(a)
    for n in 1:3
        allsum = all_possible_sums(a, n)
        m = length(a)
        for i in 1:binomial(m + n - 1, n)
            @test mina*n ≤ allsum[i][1] ≤ maxa*n
        end
    end
end

include("integer_length_tests.jl")
include("float_length_tests.jl")


text:
abstract type Perturbation end

"""
    mutable struct BranchTrip <: Perturbation
        time::Float64
        branch_type::Type{<:PowerSystems.ACBranch}
        branch_name::String
    end

A BranchTrip completely disconnects a branch from the system. Currently there is only support for static branches disconnection, `PowerSystems.Line` and `PowerSystems.Transformer2W`.
Future releases will provide support for a Dynamic Line disconnection.
**Note:** Islanding is currently not supported in `PowerSimulationsDynamics.jl`. If a `BranchTrip` isolates a generation unit, the system may diverge due to the isolated generator.

# Arguments:

- `time::Float64` : Defines when the Branch Trip will happen. This time should be inside the time span considered in the Simulation
- `branch_tipe::Type{<:PowerSystems.ACBranch}` : Type of branch disconnected
- `branch_name::String` : User defined name for identifying the branch
"""
mutable struct BranchTrip <: Perturbation
    time::Float64
    branch_type::Type{<:PSY.ACBranch}
    branch_name::String
end

"""
    mutable struct BranchImpedanceChange <: Perturbation
        time::Float64
        branch_type::Type{<:PSY.ACBranch}
        branch_name::String
        multiplier::Float64
    end

A BranchImpedanceChange change the impedance of a branch by a user defined multiplier. Currently there is only support for static branches disconnection, `PowerSystems.Line` and `PowerSystems.Transformer2W`.
Future releases will provide support for a Dynamic Line disconnection.

# Arguments:

- `time::Float64` : Defines when the Branch Impedance Change will happen. This time should be inside the time span considered in the Simulation
- `branch_tipe::Type{<:PowerSystems.ACBranch}` : Type of branch modified
- `branch_name::String` : User defined name for identifying the branch
- `multiplier::Float64` : User defined value for impedance multiplier.
"""
mutable struct BranchImpedanceChange <: Perturbation
    time::Float64
    branch_type::Type{<:PSY.ACBranch}
    branch_name::String
    multiplier::Float64
end

function _get_branch_for_perturbation(
    sys::PSY.System,
    perturbation::T,
) where {T <: Perturbation}
    if perturbation.branch_type == PSY.DynamicBranch
        error("DynamicBranch is not supported currently with perturbation $T")
    end
    branch = PSY.get_component(perturbation.branch_type, sys, perturbation.branch_name)
    if branch === nothing || !PSY.get_available(branch)
        throw(
            IS.ConflictingInputsError(
                "Branch $(perturbation.branch_type) $(perturbation.branch_name) not existent or unavailable in the system. Use PSY.show_components(system, $(perturbation.branch_type)) to check valid component names",
            ),
        )
    end
    return branch
end

function get_affect(::SimulationInputs, sys::PSY.System, pert::BranchImpedanceChange)
    branch = _get_branch_for_perturbation(sys, pert)
    mult = 0.0
    if pert.multiplier < 0.0
        throw(
            IS.ConflictingInputsError(
                "Negative multipliers are not allowed for BranchImpedanceChange perturbations",
            ),
        )
    elseif pert.multiplier < 1.0
        mult = 1.0 - pert.multiplier
    elseif pert.multiplier > 1.0
        mult = pert.multiplier - 1.0
    else
        @assert pert.multiplier == 1.0
    end

    return (integrator) -> begin
        @debug "Changing impedance line $(PSY.get_name(branch)) by a factor of $(pert.multiplier)"
        ybus_update!(integrator.p, branch, mult)
    end

    return
end

function get_affect(::SimulationInputs, sys::PSY.System, pert::BranchTrip)
    branch = _get_branch_for_perturbation(sys, pert)
    return (integrator) -> begin
        @debug "Tripping line $(PSY.get_name(branch))"
        ybus_update!(integrator.p, branch, -1.0)
    end
    return
end

function _record_change!(
    ybus::SparseArrays.SparseMatrixCSC{Float64, Int},
    bus_from_no::Int,
    bus_to_no::Int,
    n_buses::Int,
    Y11_real::Float64,
    Y11_imag::Float64,
    Y22_real::Float64,
    Y22_imag::Float64,
    Y12_real::Float64,
    Y21_real::Float64,
    Y12_imag::Float64,
    Y21_imag::Float64,
)

    # First Quadrant Real Part Changes
    ybus[bus_from_no, bus_from_no] += Y11_real
    ybus[bus_from_no, bus_to_no] += Y12_real
    ybus[bus_to_no, bus_from_no] += Y21_real
    ybus[bus_to_no, bus_to_no] += Y22_real

    # Second Quadrant Imag Part changes
    ybus[bus_from_no, bus_from_no + n_buses] += Y11_imag
    ybus[bus_from_no, bus_to_no + n_buses] += Y12_imag
    ybus[bus_to_no, bus_from_no + n_buses] += Y21_imag
    ybus[bus_to_no, bus_to_no + n_buses] += Y22_imag

    # Third Quadrant -1*Imag Part changes
    ybus[bus_from_no + n_buses, bus_from_no] -= Y11_imag
    ybus[bus_from_no + n_buses, bus_to_no] -= Y12_imag
    ybus[bus_to_no + n_buses, bus_from_no] -= Y21_imag
    ybus[bus_to_no + n_buses, bus_to_no] -= Y22_imag

    # Fourth Quadrant Real Part Changes
    ybus[bus_from_no + n_buses, bus_from_no + n_buses] += Y11_real
    ybus[bus_from_no + n_buses, bus_to_no + n_buses] += Y12_real
    ybus[bus_to_no + n_buses, bus_from_no + n_buses] += Y21_real
    ybus[bus_to_no + n_buses, bus_to_no + n_buses] += Y22_real
    return
end

function ybus_update!(
    ybus::SparseArrays.SparseMatrixCSC{Float64, Int},
    b::PSY.ACBranch,
    num_bus::Dict{Int, Int},
    mult::Float64,
)
    arc = PSY.get_arc(b)
    bus_from_no = num_bus[PSY.get_number(arc.from)]
    bus_to_no = num_bus[arc.to.number]
    n_buses = length(num_bus)

    Y_l = (1 / (PSY.get_r(b) + PSY.get_x(b) * 1im))

    Y11_real = mult * real(Y_l)
    Y11_imag = mult * (imag(Y_l) + PSY.get_b(b).from)
    Y22_real = mult * real(Y_l)
    Y22_imag = mult * (imag(Y_l) + PSY.get_b(b).to)
    Y12_real = Y21_real = -mult * real(Y_l)
    Y12_imag = Y21_imag = -mult * imag(Y_l)

    _record_change!(
        ybus,
        bus_from_no,
        bus_to_no,
        n_buses,
        Y11_real,
        Y11_imag,
        Y22_real,
        Y22_imag,
        Y12_real,
        Y21_real,
        Y12_imag,
        Y21_imag,
    )
    return
end

function ybus_update!(
    ybus::SparseArrays.SparseMatrixCSC{Float64, Int},
    b::PSY.TapTransformer,
    num_bus::Dict{Int, Int},
    mult::Float64,
)
    arc = PSY.get_arc(b)
    bus_from_no = num_bus[PSY.get_number(arc.from)]
    bus_to_no = num_bus[arc.to.number]

    Y_t = 1 / (PSY.get_r(b) + PSY.get_x(b) * 1im)
    c = 1 / PSY.get_tap(b)
    b = PSY.get_primary_shunt(b)

    Y11 = (Y_t * c^2)
    Y21 = Y12 = (-Y_t * c)
    Y22 = Y_t + (1im * b)

    Y11_real = mult * real(Y11)
    Y11_imag = mult * imag(Y12)
    Y22_real = mult * real(Y22)
    Y22_imag = mult * imag(Y22)
    Y12_real = -mult * real(Y12)
    Y21_real = -mult * real(Y21)
    Y12_imag = -mult * imag(Y12)
    Y21_imag = -mult * imag(Y21)

    _record_change!(
        ybus,
        bus_from_no,
        bus_to_no,
        n_buses,
        Y11_real,
        Y11_imag,
        Y22_real,
        Y22_imag,
        Y12_real,
        Y21_real,
        Y12_imag,
        Y21_imag,
    )
    return
end

function ybus_update!(
    ybus::SparseArrays.SparseMatrixCSC{Float64, Int},
    b::PSY.PhaseShiftingTransformer,
    num_bus::Dict{Int, Int},
    mult::Float64,
)
    arc = PSY.get_arc(b)
    bus_from_no = num_bus[PSY.get_number(arc.from)]
    bus_to_no = num_bus[arc.to.number]

    Y_t = 1 / (PSY.get_r(b) + PSY.get_x(b) * 1im)
    tap = (PSY.get_tap(b) * exp(PSY.get_α(b) * 1im))
    c_tap = (PSY.get_tap(b) * exp(-1 * PSY.get_α(b) * 1im))
    b = PSY.get_primary_shunt(b)

    Y11 = (Y_t / abs(tap)^2)
    Y12 = (-Y_t / c_tap)
    Y21 = (-Y_t / tap)
    Y22 = Y_t

    Y11_real = mult * real(Y11)
    Y11_imag = mult * imag(Y12)
    Y22_real = mult * real(Y22)
    Y22_imag = mult * imag(Y22)
    Y12_real = -mult * real(Y12)
    Y21_real = -mult * real(Y21)
    Y12_imag = -mult * imag(Y12)
    Y21_imag = -mult * imag(Y21)

    _record_change!(
        ybus,
        bus_from_no,
        bus_to_no,
        n_buses,
        Y11_real,
        Y11_imag,
        Y22_real,
        Y22_imag,
        Y12_real,
        Y21_real,
        Y12_imag,
        Y21_imag,
    )

    return
end

function ybus_update!(
    ybus::SparseArrays.SparseMatrixCSC{Float64, Int},
    fa::PSY.FixedAdmittance,
    num_bus::Dict{Int, Int},
    mult::Float64,
)
    bus = PSY.get_bus(fa)
    bus_ix = num_bus[PSY.get_number(bus)]
    ybus[bus_ix, bus_ix] += mult * fa.Y
    return
end

function ybus_update!(integrator_params, branch::PSY.ACBranch, mult::Float64)
    ybus_update!(integrator_params.ybus_rectangular, branch, integrator_params.lookup, mult)
    return
end

"""
    function NetworkSwitch(
        time::Float64,
        ybus::SparseArrays.SparseMatrixCSC{Complex{Float64}, Int},
    )

Allows to modify directly the admittance matrix, Ybus, used in the Simulation.
This allows the user to perform branch modifications, three phase faults (with impedance larger than zero) or branch trips, as long as the new Ybus provided captures that perturbation.

# Arguments:

- `time::Float64` : Defines when the Network Switch will happen. This time should be inside the time span considered in the Simulation
- `ybus::SparseArrays.SparseMatrixCSC{Complex{Float64}, Int}` : Complex admittance matrix
"""
mutable struct NetworkSwitch <: Perturbation
    time::Float64
    ybus_rectangular::SparseArrays.SparseMatrixCSC{Float64, Int}
    function NetworkSwitch(
        time::Float64,
        ybus::SparseArrays.SparseMatrixCSC{Complex{Float64}, Int},
    )
        n_bus = size(ybus)[1]
        if n_bus < 15_000
            I = PSY._goderya(ybus)
            if length(Set(I)) != n_bus
                error("The Ybus provided has islands and can't be used in the simulation")
            end
        else
            @warn(
                "The number of buses in the Ybus provided is $n_bus and is too large to be verified for connectivity"
            )
        end
        # TODO: Improve performance here
        ybus_rect = transform_ybus_to_rectangular(ybus)
        new(time, ybus_rect)
    end
end

function NetworkSwitch(time::Float64, ybus::PSY.Ybus)
    return NetworkSwitch(time, ybus.data)
end

function get_affect(::SimulationInputs, ::PSY.System, pert::NetworkSwitch)
    return (integrator) -> begin
        # TODO: This code can be more performant using SparseMatrix methods
        for (i, v) in enumerate(pert.ybus_rectangular)
            @debug "Changing Ybus network"
            integrator.p.ybus_rectangular[i] = v
        end
        return
    end
end

"""
    mutable struct ControlReferenceChange <: Perturbation
        time::Float64
        device::PowerSystems.DynamicInjection
        signal::Symbol
        ref_value::Float64
    end

A ControlReferenceChange allows to change the reference setpoint provided by a generator/inverter.

# Arguments:

- `time::Float64` : Defines when the Control Reference Change will happen. This time should be inside the time span considered in the Simulation
- `device::Type{<:PowerSystems.DynamicInjection}` : Dynamic device modified
- `signal::Symbol` : determines which reference setpoint will be modified. The accepted signals are:
    - `:P_ref`: Modifies the active power reference setpoint.
    - `:V_ref`: Modifies the voltage magnitude reference setpoint (if used).
    - `:Q_ref`: Modifies the reactive power reference setpoint (if used).
    - `:ω_ref`: Modifies the frequency setpoint.
- `ref_value::Float64` : User defined value for setpoint reference.
"""
mutable struct ControlReferenceChange <: Perturbation
    time::Float64
    device::PSY.DynamicInjection
    signal::Symbol
    ref_value::Float64

    function ControlReferenceChange(
        time::Float64,
        device::PSY.DynamicInjection,
        signal::Symbol,
        ref_value::Float64,
    )
        if signal ∉ ACCEPTED_CONTROL_REFS
            error("Signal $signal not accepted as a control reference change")
        end
        new(time, device, signal, ref_value)
    end
end

function _is_same_device(
    device1::DynamicWrapper{T},
    device2::T,
) where {T <: PSY.DynamicInjection}
    return PSY.get_name(device1) == PSY.get_name(device2)
end

function _is_same_device(
    ::DynamicWrapper{T},
    ::U,
) where {T <: PSY.DynamicInjection, U <: PSY.DynamicInjection}
    return false
end

function _find_device_index(inputs::SimulationInputs, device::PSY.DynamicInjection)
    wrapped_devices = get_dynamic_injectors(inputs)
    wrapped_device_ixs = findall(x -> _is_same_device(x, device), wrapped_devices)
    if isempty(wrapped_device_ixs)
        error(
            "Device $(typeof(device))-$(PSY.get_name(device)) not found in the simulation inputs",
        )
    end
    return wrapped_device_ixs[1]
end

function _is_same_device(
    device1::StaticWrapper{T},
    device2::U,
) where {T <: PSY.StaticInjection, U <: PSY.StaticInjection}
    if T != U
        return false
    end
    if PSY.get_name(device1) == PSY.get_name(device2)
        return true
    elseif PSY.get_name(device1) != PSY.get_name(device2)
        return false
    else
        error("comparison failed for $device1 and $device2")
    end
end

function _find_device_index(inputs::SimulationInputs, device::PSY.StaticInjection)
    wrapped_devices = get_static_injectors(inputs)
    wrapped_device_ixs = findall(x -> _is_same_device(x, device), wrapped_devices)
    if isempty(wrapped_device_ixs)
        error(
            "Device $(typeof(device))-$(PSY.get_name(device)) not found in the simulation inputs",
        )
    end
    return wrapped_device_ixs[1]
end

function get_affect(inputs::SimulationInputs, ::PSY.System, pert::ControlReferenceChange)
    wrapped_device_ix = _find_device_index(inputs, pert.device)
    return (integrator) -> begin
        wrapped_device = get_dynamic_injectors(integrator.p)[wrapped_device_ix]
        @debug "Changing $(PSY.get_name(wrapped_device)) $(pert.signal) to $(pert.ref_value)"
        getfield(wrapped_device, pert.signal)[] = pert.ref_value
        return
    end
end

"""
    mutable struct SourceBusVoltageChange <: Perturbation
        time::Float64
        device::PSY.Source
        signal_index::Int
        ref_value::Float64
    end

A `SourceBusVoltageChange` allows to change the reference setpoint provided by a voltage source.

# Arguments:

- `time::Float64` : Defines when the Control Reference Change will happen. This time should be inside the time span considered in the Simulation
- `device::Type{<:PowerSystems.Source}` : Device modified
- `signal::Int` : determines which reference setpoint will be modified. The accepted signals are:
    - `1`: Modifies the internal voltage magnitude reference setpoint.
    - `2`: Modifies the internal voltage angle reference setpoint.
- `ref_value::Float64` : User defined value for setpoint reference.
"""
mutable struct SourceBusVoltageChange <: Perturbation
    time::Float64
    device::PSY.Source
    signal_index::Int
    ref_value::Float64
end

function get_affect(inputs::SimulationInputs, ::PSY.System, pert::SourceBusVoltageChange)
    wrapped_device_ix = _find_device_index(inputs, pert.device)
    return (integrator) -> begin
        wrapped_device = get_static_injectors(integrator.p)[wrapped_device_ix]
        set_V_ref(wrapped_device, pert.ref_value)
        return
    end
end

"""
    mutable struct GeneratorTrip <: Perturbation
        time::Float64
        device::PowerSystems.DynamicInjection
    end

A `GeneratorTrip` allows to disconnect a Dynamic Generation unit from the system at a specified time.

# Arguments:

- `time::Float64` : Defines when the Generator Trip will happen. This time should be inside the time span considered in the Simulation
- `device::Type{<:PowerSystems.DynamicInjection}` : Device to be disconnected
"""
mutable struct GeneratorTrip <: Perturbation
    time::Float64
    device::PSY.DynamicInjection
end

function get_affect(inputs::SimulationInputs, ::PSY.System, pert::GeneratorTrip)
    wrapped_device_ix = _find_device_index(inputs, pert.device)
    return (integrator) -> begin
        wrapped_device = get_dynamic_injectors(integrator.p)[wrapped_device_ix]
        ix_range = get_ix_range(wrapped_device)
        @debug "Changing connection status $(PSY.get_name(wrapped_device)), setting states $ix_range to 0.0"
        if integrator.du !== nothing
            @debug "setting du $ix_range to 0.0"
            integrator.du[ix_range] .= 0.0
        end
        integrator.u[ix_range] .= 0.0
        set_connection_status(wrapped_device, 0)
    end
end

"""
    mutable struct LoadChange <: Perturbation
        time::Float64
        device::PowerSystems.ElectricLoad
        signal::Symbol
        ref_value::Float64
    end

A LoadChange allows to change the active or reactive power setpoint from a load.

# Arguments:

- `time::Float64` : Defines when the Load Change will happen. This time should be inside the time span considered in the Simulation
- `device::Type{<:PowerSystems.ElectricLoad}` : Dynamic device modified
- `signal::Symbol` : determines which reference setpoint will be modified. The accepted signals are:
    - `:P_ref`: Modifies the active power reference setpoint.
    - `:Q_ref`: Modifies the reactive power reference setpoint.
- `ref_value::Float64` : User defined value for setpoint reference.
"""
mutable struct LoadChange <: Perturbation
    time::Float64
    device::PSY.ElectricLoad
    signal::Symbol
    ref_value::Float64

    function LoadChange(
        time::Float64,
        device::PSY.ElectricLoad,
        signal::Symbol,
        ref_value::Float64,
    )
        if signal ∉ [:P_ref, :Q_ref]
            error("Signal $signal not accepted as a control reference change in Loads")
        end
        new(time, device, signal, ref_value)
    end
end

function _find_device_index(inputs::SimulationInputs, device::PSY.ElectricLoad)
    wrapped_devices = get_static_loads(inputs)
    wrapped_device_ixs = findall(x -> _is_same_device(x, device), wrapped_devices)
    if isempty(wrapped_device_ixs)
        error(
            "Device $(typeof(device))-$(PSY.get_name(device)) not found in the simulation inputs",
        )
    end
    return wrapped_device_ixs[1]
end

function get_affect(inputs::SimulationInputs, ::PSY.System, pert::LoadChange)
    wrapped_device_ix = _find_device_index(inputs, pert.device)
    return (integrator) -> begin
        wrapped_device = get_static_loads(integrator.p)[wrapped_device_ix]
        @debug "Changing $(PSY.get_name(wrapped_device)) $(pert.signal) to $(pert.ref_value)"
        getfield(wrapped_device, pert.signal)[] = pert.ref_value
        return
    end
end

"""
    mutable struct LoadTrip <: Perturbation
        time::Float64
        device::PowerSystems.ElectricLoad
    end

A `LoadTrip` allows the user to disconnect a load from the system.

# Arguments:

- `time::Float64` : Defines when the Generator Trip will happen. This time should be inside the time span considered in the Simulation
- `device::Type{<:PowerSystems.ElectricLoad}` : Device to be disconnected
"""
mutable struct LoadTrip <: Perturbation
    time::Float64
    device::PSY.ElectricLoad
end

function get_affect(inputs::SimulationInputs, ::PSY.System, pert::LoadTrip)
    wrapped_device_ix = _find_device_index(inputs, pert.device)
    return (integrator) -> begin
        wrapped_device = get_static_loads(integrator.p)[wrapped_device_ix]
        @info "Changing connection status $(PSY.get_name(wrapped_device))"
        set_connection_status(wrapped_device, 0)
        return
    end
end


text:
using ..Serde: VanillaPayloadSerde

"""
    VanillaConfig(base::BaseConfig) <: Configuration

Configuration for the baseline, where no compression is applied. The raw data
are sent over the network. This configuration instanciates a
[VanillaPayloadSerde](@ref).
"""
struct VanillaConfig <: Configuration
    base::BaseConfig
    payload_serde::VanillaPayloadSerde
    VanillaConfig(base::BaseConfig) = new(base, VanillaPayloadSerde())
end


text:
"""
    c = @sc_llvm f(args...)

LLVM IR explore.

```julia
c                  # view IR in the REPL
display(c)         # (ditto)
edit(c)            # open the IR in editor
print(c)           # print the IR
abspath(c)         # file path to the text containing the IR

c.native           # create native code explore
c.att              # (ditto)
c.intel            # create native code explore in intel syntax
eidt(c.native)
abspath(c.native)

c.cfg              # control-flow graph (CFG) visualizer
display(c.cfg)     # display CFG
edit(c.cfg.png)    # open PNG file in your editor
edit(c.cfg.svg)    # same for SVG
abspath(c.cfg.png) # file path to the PNG image
c.cfg_only
c.dom
# ... and so on; type `c.` + TAB to see all the list
```
"""
:(@sc_llvm)

"""
    ShowCode.from_llvm(ir::AbstractString)

Construct a LLVM IR explore from a snippet of LLVM IR.

For example,

```julia
a = @sc_llvm dump_module=true f(...)
b = ShowCode.from_llvm(string(a))
```

should be roughly equivalent.
"""
ShowCode.from_llvm(ir::AbstractString) =
    CodeLLVM(String(ir), true, ("f?", "Tuple{?}"), nothing)

struct CodeLLVM <: AbstractCode
    ir::String
    user_dump_module::Bool
    args::Any
    kwargs::Any
    cache::Dict{Symbol,Any}
    abspath::Base.RefValue{Union{Nothing,String}}
end

CodeLLVM(ir, user_dump_module, args, kwargs) = CodeLLVM(
    ir,
    user_dump_module,
    args,
    kwargs,
    Dict{Symbol,Any}(),
    Ref{Union{Nothing,String}}(nothing),
)

Base.string(llvm::CodeLLVM) = Fields(llvm).ir
Base.print(io::IO, llvm::CodeLLVM) = print(io, string(llvm))

function Base.abspath(llvm::CodeLLVM)
    p = Fields(llvm).abspath[]
    p === nothing || return p
    p = joinpath(mktempdir(prefix = "jl_codeviz_"), "code.ll")
    write(p, string(llvm))
    Fields(llvm).abspath[] = p
    return p
end

macro sc_llvm(args...)
    gen_call_with_extracted_types_and_kwargs(__module__, sc_llvm, args)
end

function sc_llvm(args...; dump_module = false, kwargs...)
    @nospecialize
    ir = sprint() do io
        @nospecialize
        code_llvm(io, args...; dump_module = true, kwargs...)
    end
    return CodeLLVM(ir, dump_module, args, kwargs)
end

function Base.summary(io::IO, llvm::CodeLLVM)
    f, t = Fields(llvm).args
    print(io, "CodeLLVM of ", f, " with ", t)
    return
end

function Base.show(io::IO, ::MIME"text/plain", llvm::CodeLLVM)
    @unpack ir, user_dump_module = Fields(llvm)
    summary(io, llvm)
    println(io)
    if !user_dump_module
        ir = sprint(print_main_llvm_ir, ir)
    end
    if get(io, :color, false)
        print_llvm(io, ir)
    else
        print(io, ir)
    end
    return
end

function print_main_llvm_ir(io, ir)
    input = IOBuffer(ir)
    for line in eachline(input)
        if startswith(line, "define")
            println(io, line)
            for line in eachline(input)
                println(io, line)
                if startswith(line, "}")
                    return
                end
            end
        end
    end
end

Base.propertynames(llvm::CodeLLVM) = (
    # conversions to `CodeNative`:
    :native,
    :intel,
    :att,
    # `-dot-*` passes https://llvm.org/docs/Passes.html
    :callgraph,
    :cfg,
    :cfg_only,
    :dom,
    :dom_only,
    :postdom,
    :postdom_only,
)

# TODO: make getproperty I/O-free; just return a lazy object
Base.getproperty(llvm::CodeLLVM, name::Symbol) =
    get!(Fields(llvm).cache, name) do
        if name in (:native, :intel, :att)
            syntax = name == :native ? :att : name
            CodeNative(llvm, syntax)
        else
            run_opt_dot(llvm, name)
        end
    end::Union{LLVMDot,CodeNative}

function run_opt_dot(llvm::CodeLLVM, name)
    ir = string(llvm)
    dotarg = "-dot-" * replace(string(name), "_" => "-")  # e.g., -dot-cfg
    dumpdir = mktempdir(prefix = "jl_codeviz_")
    cmd = getcmd(:opt)
    cmd = cmd0 = `$cmd $dotarg -`
    cmd = setenv(cmd0; dir = dumpdir)
    # TODO: don't hide handle error messages (if any)

    write_silently(cmd, ir)

    return LLVMDot(llvm, cmd, dumpdir)
end

abstract type AbstractLLVMDot <: AbstractCode end
AbstractLLVMDot(dot::AbstractLLVMDot) = dot
Base.getproperty(dot::AbstractLLVMDot, ext::Symbol) = LLVMDotImage(dot, ".$ext")
dotpath(dot) = abspath(AbstractLLVMDot(dot))
Base.basename(dot::AbstractLLVMDot) = basename(abspath(dot))

struct LLVMDot <: AbstractLLVMDot
    llvm::CodeLLVM
    cmd::Cmd
    dumpdir::String
end

Base.dirname(dot::LLVMDot) = Fields(dot).dumpdir

function Base.abspath(dot::LLVMDot)
    candidates = map(abspath, collect(dot))
    if length(candidates) > 1
        @debug "Ignoring some dot files." candidates[2:end]
    end
    return candidates[1]
end

struct SubLLVMDot <: AbstractLLVMDot
    dot::LLVMDot
    abspath::String
end

Base.abspath(dot::SubLLVMDot) = Fields(dot).abspath
dotpath(dot::SubLLVMDot) = abspath(dot)

function Base.collect(dot::LLVMDot)
    alldots = [
        SubLLVMDot(dot, p)
        for
        p in readdir(dirname(dot); join = true) if
        match(r".*\.jfptr_.*", basename(p)) === nothing && endswith(p, ".dot")
    ]
    sort!(alldots; by = length ∘ abspath)
    return alldots
end

Base.keys(dot::LLVMDot) = keys(pairs(dot))
Base.values(dot::LLVMDot) = values(pairs(dot))

function Base.pairs(dot::LLVMDot)
    defaultpath = abspath(dot)
    prefix, _ = splitext(basename(defaultpath))
    kvs = Dict{String,SubLLVMDot}()
    for dot in collect(dot)
        abspath(dot) == defaultpath && continue
        k, _ = splitext(basename(abspath(dot)))
        kvs[k] = dot
    end
    for k in collect(keys(kvs))
        if startswith(k, prefix)
            knew = k[length(prefix)+1:end]
            if !haskey(kvs, knew)
                kvs[knew] = pop!(kvs, k)
            end
        end
    end
    return kvs
end

stemname(x) = splitext(basename(x))[1]

function Base.getindex(dot::LLVMDot, needle::AbstractString)
    alldots = collect(dot)
    candidates = filter(d -> needle == stemname(d), alldots)
    length(candidates) == 1 && return candidates[1]

    defaultpath = abspath(dot)
    prefix, _ = splitext(basename(defaultpath))
    candidates = filter(d -> prefix * needle == stemname(d), alldots)
    length(candidates) == 1 && return candidates[1]

    candidates = filter(d ->　occursin(needle, stemname(d)), alldots)
    length(candidates) == 1 && return candidates[1]

    error("No unique match found:\n", map(stemname, candidates))
end

function Base.getindex(dot::LLVMDot, pattern::Regex)
    candidates = filter!(d ->　match(pattern, stemname(d)) !== nothing, collect(dot))
    length(candidates) == 1 && return candidates[1]
    error("No unique match found:\n", map(stemname, candidates))
end

struct LLVMDotImage <: AbstractCode
    dot::AbstractLLVMDot
    ext::String

    function LLVMDotImage(dot::AbstractLLVMDot, ext::String)
        @assert startswith(ext, ".")
        return new(dot, ext)
    end
end

AbstractLLVMDot(dotimg::LLVMDotImage) = Fields(dotimg).dot
Base.dirname(dotimg::LLVMDotImage) = dirname(AbstractLLVMDot(dotimg))

function _imagepath(dotimg::LLVMDotImage)
    @unpack dot, ext = Fields(dotimg)
    stem, _ = splitext(dotpath(dot))
    return stem * ext
end

function Base.abspath(dotimg::LLVMDotImage)
    imgpath = _imagepath(dotimg)
    ensure_dot_compile(dotpath(dotimg), imgpath)
    return imgpath
end

function Base.show(io::IO, ::MIME"text/plain", dotimg::LLVMDotImage)
    print(io, "LLVMDotImage at ")
    print(io, abspath(dotimg))
end

function dot_compile(dotpath, imgpath)
    # TODO: don't hide handle error messages (if any)
    _, ext = splitext(imgpath)
    fmt = lstrip(ext, '.')
    cmd = getcmd(:dot)
    cmd = `$cmd -o$imgpath -T$fmt $dotpath`
    @debug "Run: $cmd"
    run(pipeline(cmd, stdout = devnull, stderr = devnull, stdin = devnull))
end

function ensure_dot_compile(dotpath, imgpath)
    isfile(imgpath) || dot_compile(dotpath, imgpath)
    return
end

function showdot(io::IO, ext::AbstractString, dotpath::AbstractString)
    @assert startswith(ext, ".")
    stem, _ = splitext(dotpath)
    imgpath = stem * ext
    ensure_dot_compile(dotpath, imgpath)
    write(io, read(imgpath))
    return
end

Base.show(io::IO, ::MIME"image/png", dot::AbstractLLVMDot) =
    showdot(io, ".png", abspath(dot))
Base.show(io::IO, ::MIME"image/svg+xml", dot::AbstractLLVMDot) =
    showdot(io, ".svg", abspath(dot))


text:
f(x,y) = 2*x + 3y
f2(x,y) = 3*x + 2*y
dfdx(x,y) = ForwardDiff.derivative(x->f(x,y),x)


text:
# Michelle Medina
# June 20, 2020
# Tests Linear Interpolation Code

#TODO comment
x = [1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0]
u=sin.(x)

@test line(3.0,5.0,1.0,3.0,1.5) == 3.5

@test find_interval(x,3.5) == 3

@test interpolate(x,u,1.5) == 0.8753842058167891

# interpolateX(0.0, 1.0, 0.0, 1.0, 2.0, 3.0, 0.5)

interpolateY(0.0, 1.0, 1.0, 2.0, 0.5)

p1 = (0.0, 0.0, 0.0)
p2 = (0.0, 1.0, 1.0)
p3 = (1.0, 0.0, 2.0)
p4 = (1.0, 1.0, 3.0)
p5 = (0.5, 0.5)

p6 = (0.0, 0.0)
p7 = (0.0, 1.0)
p8 = (1.0, 0.0)
p9 = (1.0, 1.0)

@test interpolate(p1, p2, p3, p4, p5) == 1.5

# Test at corner points
@test interpolate(p1, p2, p3, p4, p6) == 0.0
@test interpolate(p1, p2, p3, p4, p7) == 1.0
@test interpolate(p1, p2, p3, p4, p8) == 2.0
@test interpolate(p1, p2, p3, p4, p9) == 3.0


text:
function FieldAttitude(arg0::Field, arg1::Attitude)
    return FieldAttitude((Field, Attitude), arg0, arg1)
end

function FieldAttitude(arg0::FieldAbsoluteDate, arg1::Frame, arg2::FieldAngularCoordinates)
    return FieldAttitude((FieldAbsoluteDate, Frame, FieldAngularCoordinates), arg0, arg1, arg2)
end

function FieldAttitude(arg0::FieldAbsoluteDate, arg1::Frame, arg2::FieldRotation, arg3::FieldVector3D, arg4::FieldVector3D)
    return FieldAttitude((FieldAbsoluteDate, Frame, FieldRotation, FieldVector3D, FieldVector3D), arg0, arg1, arg2, arg3, arg4)
end

function FieldAttitude(arg0::FieldAbsoluteDate, arg1::Frame, arg2::Rotation, arg3::Vector3D, arg4::Vector3D, arg5::Field)
    return FieldAttitude((FieldAbsoluteDate, Frame, Rotation, Vector3D, Vector3D, Field), arg0, arg1, arg2, arg3, arg4, arg5)
end

function FieldAttitude(arg0::Frame, arg1::TimeStampedFieldAngularCoordinates)
    return FieldAttitude((Frame, TimeStampedFieldAngularCoordinates), arg0, arg1)
end

function get_date(obj::FieldAttitude)
    return jcall(obj, "getDate", FieldAbsoluteDate, ())
end

function get_orientation(obj::FieldAttitude)
    return jcall(obj, "getOrientation", TimeStampedFieldAngularCoordinates, ())
end

function get_reference_frame(obj::FieldAttitude)
    return jcall(obj, "getReferenceFrame", Frame, ())
end

function get_rotation(obj::FieldAttitude)
    return jcall(obj, "getRotation", FieldRotation, ())
end

function get_rotation_acceleration(obj::FieldAttitude)
    return jcall(obj, "getRotationAcceleration", FieldVector3D, ())
end

function get_spin(obj::FieldAttitude)
    return jcall(obj, "getSpin", FieldVector3D, ())
end

function interpolate(obj::FieldAttitude, arg0::FieldAbsoluteDate, arg1::Stream)
    return jcall(obj, "interpolate", FieldAttitude, (FieldAbsoluteDate, Stream), arg0, arg1)
end

function interpolate(obj::FieldTimeInterpolable, arg0::FieldAbsoluteDate, arg1::Collection)
    return jcall(obj, "interpolate", FieldTimeInterpolable, (FieldAbsoluteDate, Collection), arg0, arg1)
end

function shifted_by(obj::FieldAttitude, arg0::RealFieldElement)
    return jcall(obj, "shiftedBy", FieldAttitude, (RealFieldElement,), arg0)
end

function shifted_by(obj::FieldAttitude, arg0::jdouble)
    return jcall(obj, "shiftedBy", FieldAttitude, (jdouble,), arg0)
end

function to_attitude(obj::FieldAttitude)
    return jcall(obj, "toAttitude", Attitude, ())
end

function with_reference_frame(obj::FieldAttitude, arg0::Frame)
    return jcall(obj, "withReferenceFrame", FieldAttitude, (Frame,), arg0)
end