text:
module LinearModels

using RobotDynamics
using LinearAlgebra
import RobotZoo: LinearModel, DiscreteLinearModel

"""
    DoubleIntegrator(D=1)

Create a LinearModel of a `D`-dimensional double integrator, such that the state dimension
is `2*D` (`D` positions and `D` velocities). The model is continuous and non-affine.

This is, in practice, identical to the `DoubleIntegrator` type in the base `RobotZoo`;
however the other is a child of `ContinuousDynamics` whereas this one is explicitly a
`LinearModel`. Both are included for reference and performance comparisons.
"""
function DoubleIntegrator(D=1)
    A = zeros(2D,2D)
    B = zeros(2D,D)
    for i = 1:D
        A[i, D+i] = 1
        B[D+i, i] = 1
    end
    return LinearModel(A, B)
end

function DiscreteDoubleIntegrator(h,D=1)
    b = h^2 / 2
    A = zeros(2D, 2D)
    B = zeros(2D, D)
    for i = 1:2D
        A[i,i] = 1.0
    end
    for i = 1:D
        A[i, i+D] = h
        B[i, i] = b
        B[i+D, i] = h
    end
    DiscreteLinearModel(h, A, B)
end

function FlexibleSatellite()
    eye(n) = Matrix(I,n,n)

    # inertia matrix
    J = diagm([1;2;3])

    # reaction wheel jacobian
    B_sc = diagm(ones(3))

    # linear momentum coupling matrix
    phi = [0 1 0;
           1 0 0;
           0 .2 -.8];

    # angular momentum coupling matrix
    delta = [0 0 1;
             0 1 0;
            -.7 .1 .1]

    # store this matrix for faster computations
    T = inv(J-delta'*delta)

    j = 3; # 3 modes

    # damping and stiffness
    zeta = [.001;.001;.001]
    Delta = [.05; .2; .125] * (2*pi)

    # damping and stiffness matrices
    C = zeros(j,j)
    K = zeros(j,j)
    for i =1:j
        C[i,i] = 2*zeta[i]*Delta[i];
        K[i,i] = Delta[i]^2;
    end

               #   mrp        w                  n                       ndot
    pdot_row = [zeros(3,3) .25*eye(3)       zeros(3,j)                 zeros(3,j)];
    wdot_row = [zeros(3,3) zeros(3,3)     T*delta'*K                  T*delta'*C];
    ndot_row = [zeros(j,3) zeros(j,3)     zeros(j,j)                  eye(j)];
    nddot_row = [zeros(j,3) zeros(j,3) (-K - delta*T*delta'*K)    (-C - delta*T*delta'*C)];

    # analytical A
    A_analytical = [pdot_row;wdot_row;ndot_row;nddot_row];

    # analytical B
    B_analytical = [zeros(3,3);
              -T*B_sc;
              zeros(j,3);
              delta*T*B_sc];

    LinearModel(A_analytical, B_analytical)
end

end

ext: jl


text:
export TrustRegion, convert_TR, convert_TR!, extract

"Exception type raised in case of error."
mutable struct TrustRegionException <: Exception
    msg::String
end

mutable struct TrustRegion{T}
    α₀::T
    α::T
    max_α::T
    acceptance_threshold::T
    increase_threshold::T
    reduce_threshold::T
    increase_factor::T
    decrease_factor::T
    max_unsuccinarow::Int

    function TrustRegion(
        α₀::T;
        max_α::T = T(1.0) / sqrt(eps(T)),
        acceptance_threshold::T = T(0.1),
        increase_threshold::T = T(0.75),
        reduce_threshold::T = T(0.1),
        increase_factor::T = T(5.0),
        decrease_factor::T = T(0.1),
        max_unsuccinarow::Int = 30,
    ) where {T}

        α₀ > T(0.0) || (α₀ = T(1.0))
        max_α > α₀ || throw(TrustRegionException("Invalid α₀"))
        (T(0.0) < acceptance_threshold < increase_threshold < T(1.0)) ||
            throw(TrustRegionException("Invalid thresholds"))
        (T(0.0) < decrease_factor < T(1.0) < increase_factor) ||
            throw(TrustRegionException("Invalid decrease/increase factors"))

        return new{T}(
            α₀,
            α₀,
            max_α,
            acceptance_threshold,
            increase_threshold,
            reduce_threshold,
            increase_factor,
            decrease_factor,
            max_unsuccinarow,
        )
    end
end

"""Compute the actual vs. predicted reduction radio ∆f/Δm, where
Δf = f - f_trial is the actual reduction is an objective/merit/penalty function,
Δq = q - q_trial is the reduction predicted by the model q of f.
We assume that q is being minimized, and therefore that Δq > 0.
"""
function compute_r(nlp, f, Δf, Δq, slope, d, xnext, gnext, robust)
    # If ever the next gradient is computed for round off errors reason, signal it.
    # # printstyled("dans compute r eltype(f) = $(eltype(f)) \n", color = :yellow)
    # # @show eltype(f)
    T = eltype(f)
    good_grad = false
    if robust & ((Δq < 10000 * eps(T)) | (abs(Δf) < 10000 * eps(T) * abs(f)))
        # trap potential cancellation errors
        grad!(nlp, xnext, gnext)
        good_grad = true
        slope_next = dot(gnext, d)

        Δf = -(slope_next + slope) / 2.0
    end
    r = Δf / Δq
    if isnan(r)
        r = 0.0
    end

    return r, good_grad, gnext
end

stop_norm(x) = norm(x, Inf)

function convert_TR(T, TR_init::TrustRegion)
    max_α_T = T(TR_init.max_α)
    acceptance_threshold_T = T(TR_init.acceptance_threshold)
    increase_threshold_T = T(TR_init.increase_threshold)
    reduce_threshold_T = T(TR_init.reduce_threshold)
    increase_threshold_T = T(TR_init.increase_threshold)
    decrease_factor_T = T(TR_init.decrease_factor)

    TRT = TrustRegion(
        T(TR_init.α),
        max_α = max_α_T,
        acceptance_threshold = acceptance_threshold_T,
        increase_threshold = increase_threshold_T,
        decrease_factor = decrease_factor_T,
    )

    return TRT
end


text:
using AbstractPlotting:
    NoConversion,
    convert_arguments,
    conversion_trait

@testset "Conversions" begin

    # NoConversion
    struct NoConversionTestType end
    conversion_trait(::NoConversionTestType) = NoConversion()

    let nctt = NoConversionTestType(),
        ncttt = conversion_trait(nctt)
        @test convert_arguments(ncttt, 1, 2, 3) == (1, 2, 3)
    end

end

@testset "functions" begin
    x = -pi..pi
    s = convert_arguments(Lines, x, sin)
    xy = s.args[1]
    @test xy[1][1] ≈ -pi
    @test xy[end][1] ≈ pi
    for (val, fval) in xy
        @test fval ≈ sin(val) atol=1f-6
    end

    x = range(-pi, stop=pi, length=100)
    s = convert_arguments(Lines, x, sin)
    xy = s.args[1]
    @test xy[1][1] ≈ -pi
    @test xy[end][1] ≈ pi
    for (val, fval) in xy
        @test fval ≈ sin(val) atol=1f-6
    end

    pts = [Point(1, 2), Point(4,5), Point(10, 8), Point(1, 2)]
    ls=LineString(pts)
    p = convert_arguments(AbstractPlotting.PointBased(), ls)
    @test p[1] == pts

    pts1 = [Point(5, 2), Point(4,8), Point(2, 8), Point(5, 2)]
    ls1 = LineString(pts1)
    lsa = [ls, ls1]
    p1 = convert_arguments(AbstractPlotting.PointBased(), lsa)
    @test p1[1][1:4] == pts
    @test p1[1][6:9] == pts1

    mls = MultiLineString(lsa)
    p2 = convert_arguments(AbstractPlotting.PointBased(), mls)
    @test p2[1][1:4] == pts
    @test p2[1][6:9] == pts1

    pol_e = Polygon(ls)
    p3_e = convert_arguments(AbstractPlotting.PointBased(), pol_e)
    @test p3_e[1] == pts

    pol = Polygon(ls, [ls1])
    p3 = convert_arguments(AbstractPlotting.PointBased(), pol)
    @test p3[1][1:4] == pts
    @test p3[1][6:9] == pts1

    pts2 = Point{2, Int}[(5, 1), (3, 3), (4, 8), (1, 2), (5, 1)]
    pts3 = Point{2, Int}[(2, 2), (2, 3),(3, 4), (2, 2)]
    pts4 = Point{2, Int}[(2, 2), (3, 8),(5, 6), (3, 4), (2, 2)]
    ls2 = LineString(pts2)
    ls3 = LineString(pts3)
    ls4 = LineString(pts4)
    pol1 = Polygon(ls2, [ls3, ls4])
    apol = [pol, pol1]
    p4 = convert_arguments(AbstractPlotting.PointBased(), apol)
    mpol = MultiPolygon([pol, pol1])
    @test p4[1][1:4] == pts
    @test p4[1][6:9] == pts1
    @test p4[1][11:15] == pts2
    @test p4[1][17:20] == pts3
    @test p4[1][22:26] == pts4
end

@testset "Categorical values" begin
    # AbstractPlotting.jl#345
    a = Any[Int64(1), Int32(1), Int128(2)] # vector of categorical values of different types
    ilabels = AbstractPlotting.categoric_labels(a)
    @test ilabels == [1, 2]
    @test AbstractPlotting.categoric_position.(a, Ref(ilabels)) == [1, 1, 2]
end

using AbstractPlotting: check_line_pattern, line_diff_pattern

@testset "Linetype" begin
    @test isnothing(check_line_pattern("-."))
    @test isnothing(check_line_pattern("--"))
    @test_throws ArgumentError check_line_pattern("-.*")

    # for readability, the length of dash and dot
    dash, dot = 3.0, 1.0

    @test line_diff_pattern(:dash)             ==
          line_diff_pattern("-",   :normal)    == [dash, 3.0]
    @test line_diff_pattern(:dot)              ==
          line_diff_pattern(".",   :normal)    == [dot, 2.0]
    @test line_diff_pattern(:dashdot)          ==
          line_diff_pattern("-.",  :normal)    == [dash, 3.0, dot, 3.0]
    @test line_diff_pattern(:dashdotdot)       ==
          line_diff_pattern("-..", :normal)    == [dash, 3.0, dot, 2.0, dot, 3.0]

    @test line_diff_pattern(:dash, :loose)     == [dash, 6.0]
    @test line_diff_pattern(:dot,  :loose)     == [dot, 4.0]
    @test line_diff_pattern("-",   :dense)     == [dash, 2.0]
    @test line_diff_pattern(".",   :dense)     == [dot, 1.0]
    @test line_diff_pattern(:dash, 0.5)        == [dash, 0.5]
    @test line_diff_pattern(:dot,  0.5)        == [dot, 0.5]
    @test line_diff_pattern("-",   (0.4, 0.6)) == [dash, 0.6]
    @test line_diff_pattern(:dot,  (0.4, 0.6)) == [dot, 0.4]
    @test line_diff_pattern("-..", (0.4, 0.6)) == [dash, 0.6, dot, 0.4, dot, 0.6]

    # gaps must be Symbol, a number, or two numbers
    @test_throws ArgumentError line_diff_pattern(:dash, :NORMAL)
    @test_throws ArgumentError line_diff_pattern(:dash, ())
    @test_throws ArgumentError line_diff_pattern(:dash, (1, 2, 3))
end


text:
abstract type PlaybackLayer{Next <: Layer} <: Layer{Next} end

const NoQuery = Dict{SubString{String}, SubString{String}}

function before(::Type{<:PlaybackLayer}, storage, path)
    data = load(storage, path)
    state = get_state()
    append!(state.responses, data[:responses])
end

function HTTP.request(::Type{<:PlaybackLayer}, m, u, h=Header[], b=nobody; kwargs...)
    method = string(m)
    uri = request_uri(u, get(kwargs, :query, nothing))
    headers = mkheaders(get(kwargs, :headers, h))
    body = get(kwargs, :body, b)
    state = get_state()
    isempty(state.responses) && error("No responses remaining in the data file")
    response = popfirst!(state.responses)
    request = response.request
    check_body(request, body)
    check_method(request, method)
    check_headers(request, headers; ignore=state.ignore_headers)
    check_uri(request, uri; ignore=state.ignore_query)
    return response
end

function after(::Type{<:PlaybackLayer}, storage, path)
    state = get_state()
    isempty(state.responses) || error("Found unused responses")
end

parse_query(uri) = isempty(uri.query) ? NoQuery() : Dict(split.(split(uri.query, '&'), '='))

function check_body(request, body)
    if request.body == body_was_streamed
        @warn "Can't verify streamed request body"
    else
        request.body == Vector{UInt8}(body) || error("Request body does not match")
    end
end

check_method(request, method) =
    request.method == method || error("Expected $(request.method) request, got $method")

function check_headers(request, headers; ignore)
    observed = filter(drop_keys(ignore), headers)
    # We only check for subset here because the `request` pipeline adds some headers
    # that the user didn't explicitly provide (Content-Length, Host, User-Agent).
    issubset(observed, request.headers) || error("Request headers do not match")
end

function check_uri(request, uri; ignore)
    # Check host.
    host = header(request, "Host")
    host == uri.host || error("Expected request to $host, got $(uri.host)")
    # Check path.
    expected = URI(request.target)
    path = isempty(uri.path) ? "/" : uri.path
    expected.path == path || error("Expected request to $(expected.path), got $(uri.path)")
    # Check query string parameters.
    expected_q = parse_query(expected)
    observed_q = filter(drop_keys(ignore), parse_query(uri))
    expected_q == observed_q || error("Query string parameters do not match")
end


text:
using Turing, ContinuousBenchmarks
using Turing: VarName, @VarName
using ForwardDiff: Dual

optRes = ">>> Runtime of ? for 1000 times <<<\n"
optRes *= "Pure function with symbol    : "

function f1()
    x = rand(Normal(0, 1))
end

t1 = @elapsed for i = 1:10000 f1() end

optRes *= "$t1\n"

optRes *= "Low-level API with symbol    : "

@model m2() = begin
    x = Turing.assume(nothing, Normal(0, 1), VarName(:x,:x,"",1), vi)[1]
end

f2 = m2()

t2 = @elapsed for i = 1:10000 f2() end

optRes *= "$t2\n"

optRes *= "~ notation with symbol       : "

@model m3() = begin
    x ~ Normal(0, 1)
end

f3 = m3()

t3 = @elapsed for i = 1:10000 f3() end

optRes *= "$t3\n"

optRes *= "~ notation with expression   : "

@model m4() = begin
    x = tzeros(1)
    x[1] ~ Normal(0, 1)
end

f4 = m4()

t4 = @elapsed for i = 1:10000 f4() end

optRes *= "$t4\n"

optRes *= "Low-level API with expression: "

@model m5() = begin
    x = tzeros(1)
    x[1] = Turing.assume(nothing, Normal(0, 1), VarName(:x,:x,"",1), vi)[1]
end

f5 = m5()

t5 = @elapsed for i = 1:10000 f5() end

optRes *= "$t5\n"

optRes *= "Pure function with expression: "

f6() = begin
    x = Vector{Real}(undef, 1)
    x[1] = rand(Normal(0, 1))
end

t6 = @elapsed for i = 1:10000 f6() end

optRes *= "$t6\n"

turnprogress(false)

optRes *= ">>> Pre-allocation? <<<\n"

run_total = 5

#= #broken

optRes *= "LDA (normal): "

include(joinpath(ContinuousBenchmarks.STAN_DATA_DIR, "lda-stan.data.jl"))
include(joinpath(ContinuousBenchmarks.STAN_MODELS_DIR, "lda.model.jl"))

data = map(ind -> ldastandata[1][ind], ["K", "V", "M", "N", "w", "doc", "beta", "alpha"])
ts = Float64[]
for run = 1:run_total
    ok = false
    chn = nothing
    counter = 0
    chn = sample(ldamodel_vec(data...), HMC(2000, 0.01, 5))
    t = sum(chn[:elapsed])
    push!(ts, t)
end
optRes *= "$ts, mean=$(mean(ts)), var=$(var(ts))\n"
=#

optRes *= "MoC (normal): "

include(joinpath(ContinuousBenchmarks.STAN_DATA_DIR, "MoC-stan.data.jl"))
include(joinpath(ContinuousBenchmarks.STAN_MODELS_DIR, "MoC.model.jl"))

ts2 = Float64[]

data = map(ind -> nbstandata[1][ind], ["K", "V", "M", "N", "z", "w", "doc", "alpha", "beta"])
for run = 1:run_total
    ok = false
    chn = nothing
    chn = sample(nbmodel(data...), HMC(100, 0.01, 5))
    t = sum(chn[:elapsed])
    push!(ts2, t)
end

optRes *= "$ts2, mean=$(mean(ts2)), var=$(var(ts2))\n"

#= #broken
optRes *= "LDA (pre-alloc): "

global theta = Matrix{Real}(undef, 2, 25)
global phi = Matrix{Real}(undef, 5, 2)

data = map(ind -> ldastandata[1][ind], ["K", "V", "M", "N", "w", "doc", "beta", "alpha"])
ts3 = Float64[]
for run = 1:run_total
    ok = false
    chn = nothing
    chn = sample(ldamodel_vec(data...), HMC(100, 0.01, 5))
    ok = true
    t = sum(chn[:elapsed])
    push!(ts3, t)
end

optRes *= "$ts3, mean=$(mean(ts3)), var=$(var(ts3))\n"
=#

optRes *= "MoC (pre-alloc): "

theta = Vector{Real}(undef, 4)
phi = Vector{Vector{Real}}(undef, 4)
for k = 1:4
    phi[k] = Vector{Real}(undef, 10)
end
log_theta = Vector{Real}(undef, 4)
log_phi = Vector{Real}(undef, 4)

@model nbmodel(K, V, M, N, z, w, doc, alpha, beta) = begin
    theta ~ Dirichlet(alpha)
    # phi = Array{Any}(K)
    for k = 1:K
        phi[k] ~ Dirichlet(beta)
    end

    log_theta = log.(theta)
    Turing.acclogp!(vi, sum(log_theta[z]))

    log_phi = map(x->log.(x), phi)
    lp = mapreduce(n->log_phi[z[doc[n]]][w[n]], +, 1:N)
    Turing.acclogp!(vi, lp)

    phi
end

ts4 = Float64[]
data = map(ind -> nbstandata[1][ind], ["K", "V", "M", "N", "z", "w", "doc", "alpha", "beta"])
for run = 1:run_total
    chn = sample(nbmodel(data...), HMC(100, 0.01, 5))
    t = sum(chn[:elapsed])
    push!(ts4, t)
end

optRes *= "$ts4, mean=$(mean(ts4)), var=$(var(ts4))\n"

optRes *= ">>> Gen VarName <<<\n"

t = @elapsed for i = 1:1000
    sym, idcs, csym = @VarName x[i]
    str = reduce(*, map(idx -> string(idx), idcs), init="")
    VarName(csym, sym, str, 1)
end

optRes *= "1st Run: \n"
optRes *= "Time : $t\n"

t = @elapsed for i = 1:1000
    sym, idcs, csym = @VarName x[i]
    str = reduce(*, map(idx -> string(idx), idcs), init="")
    VarName(csym, sym, str, 1)
end

optRes *= "2nd Run: \n"
optRes *= "Time : $t\n"

using ForwardDiff

optRes *= "realpart(): \n"

ds = [Dual{Nothing,Float64,10}(rand()) for i = 1:1000]

t_map = @elapsed for i = 1:1000 map(d -> d.value, ds) end
t_list = @elapsed for i = 1:1000 Float64[ds[i].value for i = 1:length(ds)] end

optRes *= "Map realpart: $t_map\n"
optRes *= "List realpart: $t_list\n"

optRes *= "Constructing Dual numbers: \n"

SEEDS = ForwardDiff.construct_seeds(ForwardDiff.Partials{44,Float64})
dps = SEEDS[11]

t_dualnumbers = @elapsed for _ = 1:(44*2000*5) ForwardDiff.Dual{Nothing, Float64, 44}(1.1, dps) end

optRes *= "44*2000*5 times: $t_dualnumbers\n"


text:
using OrdinaryDiffEq, DiffEqDevTools, DiffEqBase, Test, Random
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_2Dlinear, prob_ode_linear

Random.seed!(100)
## Convergence Testing
dts = (1/2) .^ (7:-1:4)
testTol = 0.2

prob = prob_ode_linear
sol = solve(prob,OwrenZen3())
@test length(sol) < 20
sol = solve(prob,OwrenZen4())
@test length(sol) < 20
sol = solve(prob,OwrenZen5())
@test length(sol) < 20

sim = test_convergence(dts,prob,OwrenZen3(),dense_errors=true)
@test sim.𝒪est[:final] ≈ 3 atol=testTol
@test sim.𝒪est[:L2] ≈ 3 atol=testTol
sim = test_convergence(dts,prob,OwrenZen4(),dense_errors=true)
@test sim.𝒪est[:final] ≈ 4 atol=testTol
@test sim.𝒪est[:L2] ≈ 4 atol=testTol
sim = test_convergence(dts,prob,OwrenZen5(),dense_errors=true)
@test sim.𝒪est[:final] ≈ 5 atol=testTol
@test sim.𝒪est[:L2] ≈ 5 atol=testTol

prob = prob_ode_2Dlinear
sol = solve(prob,OwrenZen3())
@test length(sol) < 20
sol = solve(prob,OwrenZen4())
@test length(sol) < 20
sol = solve(prob,OwrenZen5())
@test length(sol) < 20

sim = test_convergence(dts,prob,OwrenZen3(),dense_errors=true)
@test sim.𝒪est[:final] ≈ 3 atol=testTol
@test sim.𝒪est[:L2] ≈ 3 atol=testTol
sim = test_convergence(dts,prob,OwrenZen4(),dense_errors=true)
@test sim.𝒪est[:final] ≈ 4 atol=testTol
@test sim.𝒪est[:L2] ≈ 4 atol=testTol
sim = test_convergence(dts,prob,OwrenZen5(),dense_errors=true)
@test sim.𝒪est[:final] ≈ 5 atol=testTol
@test sim.𝒪est[:L2] ≈ 5 atol=testTol


text:
export AlleleMask, AlleleString, AlleleLinks, AlleleContribs

@doc """AlleleMask UInt128

A bit mask that identifies one or several alleles within a genotype.
"""
typealias AlleleMask UInt128

@doc """AlleleString UInt128

A bit string that represents the alleles of a particular genotype.
"""
typealias AlleleString UInt128

@doc """AlleleLinks Vector{AlleleMask}

A vector of masks specifying the epistatic links for each allele. The ith
element of the vector stores a mask of alleles that are epistatically linked
to the allele in the ith least significant position.

Each mask should have a 1 in the position that corresponds to the locus
discussed above, and in the K positions that correspond to the loci to which it
is epistatically linked.
"""
typealias AlleleLinks Vector{AlleleMask}

@doc """AlleleContribs Vector{Dict{AlleleString, Float64}}

A data structure in which fitness contributions may be stored. The ith element
of the vector stores contributions of the allele in the ith least significant
position. For example, the 4th element of the vector corresponds to the locus
which contains a 1 bit in the bit string 1000.

Each dictionary then maps a bit string into a fitness contribution. The bit string should contain all zeros except for the following:

  * The locus that corresponds to the dictionary should contain its actual
    value, either 1 or 0
  * The K loci to which the locus above is linked should contain their actual
    values, either 1 or 0
"""
typealias AlleleContribs Vector{Dict{AlleleString, Float64}}

newAlleleContribs(n::Int64) = [Dict{AlleleString, Float64}() for _ = 1:n]


text:
function triangularloss(d::Vector{T},y::Vector{I}) where {T<:Number,I<:Integer}
	l = length(y)
	e = zero(T)
	for i in 1:l
		for j in setdiff(findall( y.== y[i]),i)
			for k in findall( y.!= y[i])
				e += max(0, 1 + abs(d[i] - d[j]) - abs(d[i] - d[k]))
			end
		end
	end
	e
end

function triangularloss_back(d, y, Δ)
	l = length(y)
	e = 0.0
	g = zero(d)
	for i in 1:l
		for j in setdiff(findall( y.== y[i]),i)
			for k in findall( y.!= y[i])
				if 1 + abs(d[i] - d[j]) - abs(d[i] - d[k]) > 0
					g[i] += sign(d[i] - d[j]) - sign(d[i] - d[k])
					g[j] += -sign(d[i] - d[j])
					g[k] += sign(d[i] - d[k])
				end
			end
		end
	end
	g.*Δ
end

triangularloss(d::Flux.Tracker.TrackedArray, y) = Flux.Tracker.track(triangularloss, d, y)
Flux.Tracker.@grad function triangularloss(d, y)
  return(triangularloss(Flux.data(d), y), Δ -> (triangularloss_back(Flux.data(d), y, Δ),nothing))
end

max_stars_repo_name: pevnak/FluxExtensions

max_issues_repo_name: pevnak/FluxExtensions

max_forks_repo_name: pevnak/FluxExtensions


text:
#using Unitful, ImagineInterface, ImagineAnalyses, Images, ProgressMeter
#import Unitful:s, μm, V
#
##calculate lag mapping the piezo mon signal to the hardware response
##Procedure:
##   The idea is to use the camera as an independent measurement device instead of the piezo sensor signal
##   Delivere several paired-exposure pulse trials with bead sample:
##        first pulse always occurs in the middle of the forward stroke of the piezo (bidirectional motion)
##        second pulse always occurs in the middle of the reverse stroke of the piezo, offset by a set of lag times.
##        ssd on pairs of images to determine at what lag time ssd is minimized.  That is the monitor lag.
#
##lags is which lag values to check
#function mon_lag_sigs0(rig::AbstractString, pos_name::AbstractString, pos_start = 540μm, pos_stop = 630μm, stack_rate = 10*inv(Unitful.s))
#    sample_rate = 100000 * inv(Unitful.s)
#    rig_sigs = rigtemplate(rig; sample_rate = sample_rate)
#    posfwd, posback = gen_bidi_pos(pos_start, pos_stop, 1/stack_rate, sample_rate)
#    pos = getname(rig_sigs, pos_name)
#    add_sequence!(pos, "bidi_cycle", vcat(posfwd, posback))
#    cycle_nsamps = length(posfwd) + length(posback)
#    cam = first(getcameras(rig_sigs))
#    las = first(getlasers(rig_sigs))
#    add_sequence!(cam, "quiet", falses(cycle_nsamps)) #just to pass validation
#    ncycles = 300
#    for i = 1:ncycles
#        append!(pos, "bidi_cycle")
#        append!(cam, "quiet")
#        append!(las, "quiet")
#    end
#    return [pos;cam;las]
#end
#
##take the mean MON cycle and calculate the number of sampls that MON is shifted from MOD
#function mon_lag_nsamps(pos, pos_mon, ncycs_ignore)
#    pos_cyc_name = first(sequence_names(pos))
#    pos_cycle = get_samples(pos, pos_cyc_name)
#    @assert all(sequence_names(pos).==pos_cyc_name)
#    nsamps_cycle = Base.length(pos_cycle) #should always be even
#    mon_samps = get_samples(pos_mon)
#    mon_samps = view(mon_samps, (ncycs_ignore*nsamps_cycle+1):Base.length(mon_samps))
#    cycs_mon = ImagineAnalyses.get_cycles(mon_samps, nsamps_cycle)
#    mean_cyc = mean(ustrip.(cycs_mon), 2)[:] * unit(mon_samps[1])
#    return ImagineAnalyses.mon_delay(pos_cycle, mean_cyc)
#end
#
#function gen_mon_lag_procedure(coms, recs, las_name::AbstractString, cam_name::AbstractString)
#    pos = getpositioners(coms)
#    if length(pos) > 1
#        error("Provide only one positioner command.")
#    end
#    pos = pos[1]
#    pos_mon = getpositionermonitors(recs)
#    if length(pos_mon) > 1
#        error("Provide only one positioner monitor.")
#    end
#    pos_mon = pos_mon[1]
#    nsamps_cycle = full_length(first(sequences(pos)))
#    ncycs_ignore = 50
#    ntrials=10
#    @show nsamps_offset = mon_lag_nsamps(pos, pos_mon, ncycs_ignore)
#    mon_samps = get_samples(pos_mon)
#    mean_cyc = mean_cycle(mon_samps, nsamps_cycle; start_idx = ncycs_ignore*nsamps_cycle+1)
#    mean_cyc = circshift(mean_cyc, -nsamps_offset)
#    #@show lags = [0.0002s:0.00005s:0.0007s...]  #These worked to find the lag on the analog piezo of OCPI2
#    @show lags_pre = [0.0005s:0.00005s:0.002s...] #These worked to find the lag on the digital piezo of OCPI2
#    gen_desc = "Generate a set camera and laser commands to measure positioner monitor lag given a completed dynamic positioner recording"
#    sig_gen = SignalGenerator(gen_desc, mon_lag_sigs, (pos, mean_cyc, las_name, cam_name, lags_pre, ncycs_ignore, ntrials, nsamps_offset))
#    rig = rig_name(pos)
#    mon_name = name(pos_mon)
#    desc = "Calculates lag mapping the $rig $mon_name sensor signal to its true hardware response.  Returns a Calibration."
#    analyze_f = (coms, img) -> calc_lag(coms, img, lags_pre, ncycs_ignore, mean_cyc, nsamps_offset)
#    return ImagineProcedure(desc, sig_gen, analyze_f)
#end
#
#function calc_lag(coms, img, lags, ncycs_ignore::Int, mean_cyc, nsamps_offset::Int)
#    #mean_cyc has already been shifted by offset so that it aligns well with the command
#    sampr = samprate(first(coms))
#    fwd_ctr, back_ctr = find_bidi_centers(mean_cyc, sampr)
#    lag_pre = select_lag_pre(coms, img, lags, ncycs_ignore)
#    #NOTE: the lag chosen here isn't the signal lag, it needs more processing to find that.
#    #The lag chosen here is the temporal displacement from the MON measurement at the exposure moment to the moment when the
#    #forward and reverse sweeps cross the same location (that crossing takes place somewhere between the corrected z location
#    #and the original target z location
#    #so the true value is roughly half of this measurement
#    #a more precise way to find the true value is to advance two pointers until they reach the same z value:
#    #the first pointer begins at the fwd flash location and steps forward in time, the second
#    #pointer begins at the corrected second flash location and steps forward in time.
#    #The (equal) time interval traveled by the pointers when they meet is the lag
#    fwd_i = fwd_ctr
#    back_i = back_ctr - ImagineInterface.calc_num_samps(lag_pre, sampr)
#    lag_nsamps = 0
#    while back_i < length(mean_cyc)
#        if mean_cyc[fwd_i] < mean_cyc[back_i]
#            lag_nsamps+=1
#            fwd_i+=1
#            back_i+=1
#        else
#            return Calibration(upreferred(lag_nsamps / sampr), 0.0V, 1.0)
#        end
#    end
#    error("Couldn't find a lag value that makes sense.  Check inputs.")
#end
#
#
##Analyzes pairs of images to find the lag (in lags) at which image difference is minimized
##The lags vector is assumed to be ordered in the same way that trials were ordered when acquiring imgs
#function select_lag_pre(coms, img, lags, ncycs_ignore)
#    pos = first(getpositioners(coms))
#    ncycs = length(sequence_names(pos))
#    npairs = size(img,3) / 2
#    if !isinteger(npairs)
#        error("Expected an odd number of images (each trial requires an image pair).")
#    elseif npairs < 1
#        error("Not enough cycles.")
#    else
#        npairs = Int(npairs)
#    end
#    if abs(ncycs_ignore - abs(ncycs - npairs)) > 1 #allow a difference of one because sometimes an extra cycle gets added when the MON signal lags a lot
#        error("The number of cycles in the command signal disagrees with the image count and ncycs_ignore")
#    end
#    nlags = length(lags)
#    ntrials = npairs/nlags
#    if !isinteger(ntrials)
#        error("The number of image pairs should be evenly divisible by the number of lag conditions.")
#    else
#        ntrials = Int(ntrials)
#    end
#    sz_spatial = size(img)[1:2]
#    img = reshape(img, sz_spatial..., 2, ntrials, nlags)
#    ssds = zeros(ntrials, nlags)
#    #TODO: should also do a rigid 2D registration when comparing images because there is lateral wobble in fwd vs reverse images
#    @showprogress for l = 1:nlags
#        for t = 1:ntrials
#            ssds[t,l] = ssd(img[:,:,1,t,l], img[:,:,2,t,l])
#        end
#    end
#    @show mean_ssds = mean(ssds, 1)[:]
#    return lags[indmin(mean_ssds)]
#end
#
#function find_bidi_centers(cyc::AbstractVector, sampr::HasInverseTimeUnits)
#    @show pmin = minimum(cyc)
#    @show pmax = maximum(cyc)
#    prange = pmax - pmin
#    pctr = pmin + prange/2 #we will place the flash in the center
#    stack_rate = 2 * (sampr / length(cyc))
#    pos_avg_speed = stack_rate * prange #in microns per second
#    pad_secs = (1/stack_rate)/2 #wait half a cycle before accepting threshold crossings (helps with noise)
#    pad_nsamps = round(Int, pad_secs * sampr)
#    mon_ctr_is = ImagineAnalyses.find_circular(cyc, [pctr;], pad_nsamps)[1]
#    @assert Base.length(mon_ctr_is) == 2 #should have just two crossings, one fwd one back
#    fwd_ctr = mon_ctr_is[1]
#    back_ctr = mon_ctr_is[2]
#    return fwd_ctr, back_ctr
#end
#
##ncycs_ignore determines how many to skip before starting exposure pulses
##...and also how many cycles of pos_mon to skip when calculating the mean cycle waveform
#function mon_lag_sigs(pos::T1, mean_cyc, las_name::AbstractString, cam_name::AbstractString, lags, ncycs_ignore, ntrials, nsamps_offset::Int) where {T1<:ImagineSignal}
#    exp_time = 0.011s
#    flash_time = 0.001s
#    pos = deepcopy(pos)
#    rig = rig_name(pos)
#    #We want to deliver the laser pulse during the global shutter time (tglobal in PCO camera manual)
#    #If the user images with full chip then the global exposure segment doesn't begin until _10ms_ after the exposure starts.
#    #Therefore with a 1ms laser pulse we need an 11ms exposure time, and the pulse is placed in the last 1ms
#    rig_sigs = rigtemplate(rig; sample_rate = samprate(pos))
#    cam = getname(rig_sigs, cam_name)
#    las = -1
#    high_las = -1
#    all_las = -1
#    if rig == "ocpi-2"
#        if las_name != "all lasers"
#            las = getname(rig_sigs, "all lasers")
#            high_las = getname(rig_sigs, las_name)
#            all_las = [las;high_las]
#        else
#            error("Please specify a specific laser, not all lasers")
#        end
#    else
#        las = getname(rig_sigs, las_name)
#        all_las = [las;]
#    end
#
#    if !ImagineInterface.ispos(pos)
#        error("Only positioner signals are currently supported.")
#    end
#    if !ImagineInterface.iscam(cam)
#        error("$cam_name is not a camera trigger signal.")
#    end
#    if !ImagineInterface.islas(las)
#        error("$las_name is not a laser trigger signal.")
#    end
#
#    pos_cyc_name = first(sequence_names(pos))
#    pos_cycle = get_samples(pos, pos_cyc_name)
#    @assert all(sequence_names(pos).==pos_cyc_name)
#    empty!(pos)
#
#    nsamps_cycle = Base.length(pos_cycle) #should always be even
#    add_sequence!(cam, "wait_to_start", falses(nsamps_cycle))
#    if rig == "ocpi-2"
#        add_sequence!(high_las, "laser_high", trues(nsamps_cycle))
#    end
#    for i = 1:ncycs_ignore
#        append!(pos, pos_cyc_name)
#        append!(cam, "wait_to_start")
#        append!(las, "wait_to_start")
#        append!(high_las, "laser_high")
#    end
#    #After the ignored cycles each trial has this structure:
#    #pos: pos_cyc_name
#    #cam: "fwd_exp" -> "back_exp_$i" where i is an index into the lags vector
#    #las: "fwd_flash" -> "back_flash_$i"
#    #in order to know when/where to place the camera and laser pulses we must look at the pos_mon signal
#
#    fwd_ctr, back_ctr = find_bidi_centers(mean_cyc, samprate(pos))
#
#    nsamps_flash = ImagineInterface.calc_num_samps(flash_time, samprate(pos))
#    if iseven(nsamps_flash)
#	nsamps_flash -= 1 #force odd so that fwd and reverse flashes are easier to align
#    end
#    nsamps_from_ctr = div(nsamps_flash,2)
#    #forward flash and exposure are always timed the same, so we can calculate them just once
#    flash_itv_fwd = ClosedInterval(fwd_ctr-nsamps_from_ctr, fwd_ctr+nsamps_from_ctr)
#    nsamps_exp = ImagineInterface.calc_num_samps(exp_time, samprate(pos))
#    #This aligns the end of the exposure with the end of the flash (allows global exposure)
#    exp_itv_fwd = ClosedInterval(fwd_ctr+nsamps_from_ctr-nsamps_exp-1, fwd_ctr+nsamps_from_ctr)
#    nsamps_sweep = div(nsamps_cycle,2)
#    fwd_flash = ImagineInterface.gen_pulses(nsamps_sweep, [flash_itv_fwd;])
#    fwd_exp = ImagineInterface.gen_pulses(nsamps_sweep, [exp_itv_fwd;])
#    if nsamps_offset < 0
#        error("Cannot handle negative monitor lag")
#    elseif nsamps_offset > 0
#        append!(las, "delay", fill(false, nsamps_offset))
#        append!(cam, "delay")
#        append!(high_las, "delayhigh", fill(true, nsamps_offset))
#    end
#    add_sequence!(las, "fwd_flash", fwd_flash)
#    add_sequence!(cam, "fwd_exp", fwd_exp)
#    for l = 1:Base.length(lags)
#        #calculate flash and exposure samples for backward sweep
#        back_exp_name = "back_exp_$l"
#        back_flash_name = "back_flash_$l"
#        #A lag of 0.0s would mean that the laser pulse is centered on the measured stack center (back_ctr)
#        #Positive lags mean that the pulse centers precede back_ctr
#        nsamps_lag = ImagineInterface.calc_num_samps(lags[l], samprate(pos))
#        flash_itv_back = ClosedInterval(back_ctr-nsamps_lag-nsamps_from_ctr-nsamps_sweep, back_ctr-nsamps_lag+nsamps_from_ctr-nsamps_sweep)
#        exp_itv_back = ClosedInterval(back_ctr-nsamps_lag+nsamps_from_ctr-nsamps_exp-1-nsamps_sweep, back_ctr-nsamps_lag+nsamps_from_ctr-nsamps_sweep)
#        back_flash = ImagineInterface.gen_pulses(nsamps_sweep, [flash_itv_back;])
#        back_exp = ImagineInterface.gen_pulses(nsamps_sweep, [exp_itv_back;])
#        add_sequence!(cam, back_exp_name, back_exp)
#        add_sequence!(las, back_flash_name, back_flash)
#        for t = 1:ntrials
#            append!(pos, pos_cyc_name)
#            append!(cam, "fwd_exp")
#            append!(cam, back_exp_name)
#            append!(las, "fwd_flash")
#            append!(las, back_flash_name)
#            if rig == "ocpi-2"
#                append!(high_las, "laser_high")
#            end
#        end
#    end
#    if nsamps_offset > 0 #add another positioner cycle because we have too many cam and las samples
#        append!(pos, pos_cyc_name)
#        append!(cam, "extra_cyc", fill(false, nsamps_cycle-nsamps_offset))
#        append!(las, "extra_cyc")
#        append!(high_las, "extra_cychigh", fill(true, nsamps_cycle-nsamps_offset))
#    end
#    return [pos; cam; all_las...]
#end
#
##create procedures for each positioner on each rig
#rigs = ImagineInterface.RIGS
#rigs = filter(x->x!="dummy-6002", rigs)
#for rig in rigs
#    rig_sigs = rigtemplate(rig)
#    pos_mons = getpositionermonitors(rig_sigs)
#    for mon in pos_mons
#        mon_name = name(mon)
#        mod_name = ImagineInterface.actuator_name(mon)
#        desc = "Produces a procedure for measuring a lag mapping the $rig $mon_name sensor signal to its true hardware response.  Returns an ImagineProcedure."
#        sig_gen = SignalGenerator("Generate piezo command on which monitor lag calculation will be based.", mon_lag_sigs0, (rig, mod_name))
#        analyze_f = gen_mon_lag_procedure
#        p = ImagineProcedure(desc, sig_gen, analyze_f)
#        #now store it
#        store_calibration_procedure(rig, daq_channel(mon), p)
#    end
#end


text:
function createPoint2d(n::Int64)
    points,weights  = computeXYZandWeightsdim2(n)
	return points,weights
end

function computeXYZandWeightsdim2(_norder::Int64)
        n = _norder
        pointsxyz = zeros(n,3); # Even though we only need the x and y coordinate
        weights = zeros(n);

        # around z axis equidistant
		phi = [k/n*2*pi for k=0:n-1]

		for i=1:n
			pointsxyz[i,1] = cos(phi[i])
			pointsxyz[i,2] = sin(phi[i])
			pointsxyz[i,3] = 0.0

			weights[i] = 4*pi/n
		end

        return pointsxyz, weights
    end

max_stars_repo_name: vavrines/CSN

max_issues_repo_name: vavrines/CSN

max_forks_repo_name: vavrines/CSN


text:
export
    object_rts_low

"""
    object_rts_low(jd, ra, decl, h0)

Compute rise, set, and transit times with low accuracy.

Requires equatorial coordinates for the object in question for today at 0hr UT.

Return the Julian Days of the rise/transit/set times of an object.

Parameters:
    jd      : Julian Day number of the day in question, at 0 hr UT
    raList  : right ascension value, in radians, for jd
    decList : declination value, in radians, for jd
    h0      : the standard altitude in radians; the 'standard' altitude is
              deg2rad(-0.5667) for stars/planets, deg2rad(-0.8333) for the sun
    rise/set: boolean flag to calculate Rise (true) or Set (false)

Returns:
    tuple of Julian Day rise/transit/set times

This is the low accuracy part of Meeus' method from chapter 15.

The high accuracy interpolation routine is still to be done.
"""
function object_rts_low(jd, ra, decl, h0)
    global seconds_per_day, latitude, longitude
    THETA0 = mean_sidereal_time_greenwich(jd)
    deltaT_days = deltaT_seconds(jd) / seconds_per_day
    cosH0 = (sin(h0) - sin(latitude) * sin(decl)) / (cos(latitude) * cos(decl))

    # TODO object_rts_low: return some indicator when the
    # object is circumpolar or always below the horizon.

    if cosH0 < -1.0 # circumpolar
        println("circumpolar")
        return (false, false, false)
    end
    if cosH0 > 1.0  # never rises
        println("never rises")
        return (false, false, false)
    end

    H0 = acos(cosH0)
    # transit
    m0 = (ra + longitude - THETA0) / 2pi
    # must be between 0 and 1
    # rise
    m0 = mod(m0, 1)
    m1 = mod(m0 - (H0 / 2pi),1)
    # set
    m2 = mod(m0 + (H0 / 2pi),1)
    return(m1, m0, m2)
end


text:
using TensorNetworkCodes
using TensorNetworkCodes.TNDecode
using Test

@testset "tn_decode" begin
    code = rotated_surface_code(3)
    p = 0.2;  # error_probability
    error = pauli_random_operator(num_qubits(code), p)
    syndrome = find_syndrome(code, error);
    # default_contract
    recovery, success_prob = tn_decode(code, syndrome, p)
    @test find_syndrome(code, recovery) == syndrome
    @test 0 <= success_prob <= 1
    # basic_contract (same as default)
    recovery, success_prob = tn_decode(code, syndrome, p; contract_fn=basic_contract())
    @test find_syndrome(code, recovery) == syndrome
    @test 0 <= success_prob <= 1
    # mps_contract
    recovery, success_prob = tn_decode(code, syndrome, p; contract_fn=mps_contract(8))
    @test find_syndrome(code, recovery) == syndrome
    @test 0 <= success_prob <= 1
end


text:
using SunAsAStar
using Documenter

DocMeta.setdocmeta!(SunAsAStar, :DocTestSetup, :(using SunAsAStar); recursive=true)

makedocs(;
    modules=[SunAsAStar],
    authors="Eric Ford, Andrea Lin, Shubham Kanodia",
    repo="https://github.com/RvSpectML/SunAsAStar.jl/blob/{commit}{path}#{line}",
    sitename="SunAsAStar.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://RvSpectML.github.io/SunAsAStar.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/RvSpectML/SunAsAStar.jl",
    devbranch="main",
)


text:
using Lens

struct Loop end
struct LoopEnd end
function bubblesort(a::AbstractArray{T,1}) where T
  b = copy(a)
  isordered = false
  span = length(b)
  i = 0
  while !isordered && span > 1
    lens(Loop, (b = b, i = i)) # <--- lens here!!
    isordered = true
    for i in 2:span
      if b[i] < b[i-1]
        t = b[i]
        b[i] = b[i-1]
        b[i-1] = t
        isordered = false
      end
    end
    span -= 1
    i += 1
  end
  lens(LoopEnd, (sorteddata = b, niters = i)) # <--- and here!!
  return b
end

function bubblesortnolens(a::AbstractArray{T,1}) where T
  b = copy(a)
  isordered = false
  span = length(b)
  i = 0
  while !isordered && span > 1
    isordered = true
    for i in 2:span
      if b[i] < b[i-1]
        t = b[i]
        b[i] = b[i-1]
        b[i-1] = t
        isordered = false
      end
    end
    span -= 1
    i += 1
  end
  return b
end

function test()
  @leval bubblesort([1, 4, 2]) (after_loop = (println,),)
end

test()


text:
## readme script
using Pkg;Pkg.activate("./")
using VortexDistributions, Plots
gr(xlabel="x",ylabel="y",transpose=true,legend=false)

## make a simple 2D test field
Nx = 400; Ny = Nx
Lx = 200; Ly = Lx
x = LinRange(-Lx / 2, Ly / 2, Nx); y = x
psi0 = one.(x*y') |> complex

## doubly periodic boundary conditions
psi = Torus(psi0,x,y)

# make a point vortex
pv = PointVortex(30.0,70.3,-1)

# make a scalar GPE vortex with exact core
spv = ScalarVortex(pv)
vortex!(psi,spv)

## make some more random vortices
vort = rand_vortex(10,psi)
vortex!(psi,vort)

using BenchmarkTools
@btime findvortices(psi)


text:
module PyPSANetwork

using AxisArrays

using PyCall
import Pandas
const PD = Pandas
# const PD = Base.require(Base.PkgId(Base.UUID("eadc2687-ae89-51f9-a5d9-86b5a6373a9c"), "Pandas"))

using DataFrames
using Dates

using ..EnergyModels
const EM = EnergyModels

using ..EnergyModels:
    AbstractData, Data, pushaxes!, resolve, @consense, Component, Device,
    DeviceFormulation, attributes, components, naming, astype, typenames

import ..EnergyModels: store_results!

using JuMP: value, dual, ConstraintRef, VariableRef, GenericAffExpr, LowerBoundRef, UpperBoundRef

export PypsaNetworkData

function EM.load(network::PyObject)::AbstractData
    pypsa = pyimport("pypsa")
    if pyisinstance(network, pypsa.Network)
        Data(fallback=PypsaNetworkData(network))
    else
        error("Python object is no PyPSA Network")
    end
end

abstract type AbstractAttrDesc end

struct AttrDesc{T,S} <: AbstractAttrDesc
    data::S
    indices::Union{Colon,UnitRange{Int64},StepRange{Int64,Int64},Vector{Int64}}
    default::Union{T,Missing}

    AttrDesc{T}(data::S, indices, default) where {T,S} = new{T,S}(data, indices, default)
end

struct SingleAttrDesc{T} <: AbstractAttrDesc
    value::T
end

struct ComponentDesc
    name::Symbol
    axis::Axis
    componenttype
    attributes::Dict{Symbol,AbstractAttrDesc}
    typeparams::Union{Nothing,Dict{Symbol,Any}}
end

struct PypsaNetworkData <: AbstractData
    network::PyObject
    components::Dict{Symbol}{ComponentDesc}
    axes::Dict{Symbol}{Axis}

    function PypsaNetworkData(network, components)
        snapshots = PD.values(PD.Index(network.snapshots))
        if pyisinstance(network.snapshots, pyimport("pandas").DatetimeIndex)
            snapshots = DateTime.(snapshots)
        end
        axes = Dict(:snapshots => Axis{:snapshots}(snapshots),
                    (name => cd.axis for (name, cd) in components)...)
        new(network, components, axes)
    end
end

splitgroup(df, group) = [group]
function splitgroup(df, group, attr, attrs...)
    gdf = groupby(DataFrame(attr=>df[group.idx,attr], :idx=>group.idx, copycols=false), attr);

    if length(gdf) == 1
        ng = attr == :carrier ? (idx=group.idx, name=(group.name..., first(parent(gdf)[!,attr]))) : group
        groups = splitgroup(df, ng, attrs...)
    else
        groups = Any[]
        for sdf in gdf
            ng = (idx=collect(sdf.idx),
                  name=(group.name..., first(sdf[!,attr])))
            append!(groups, splitgroup(df, ng, attrs...))
        end
    end

    groups
end

function AttrDesc(::Type{Val{:single}}, attr, name, df, pnl, idx, axis)
    val = @consense(PD.values(df[attr[:PyPSA]])[idx],
                    "$(attr[:PyPSA]) for group $name must be a single value")
    SingleAttrDesc(astype(attr[:dtype], val))
end

AttrDesc(::Type{Val{:static}}, attr, name, df, pnl, idx, axis) =
    AttrDesc{typenames[attr[:dtype]]}(df[attr[:PyPSA]], idx, attr[:default])

function AttrDesc(::Type{Val{:series}}, attr, name, df, pnl, idx, axis)
    df_t = get(pnl, attr[:PyPSA], nothing)
    indices_t = !isnothing(df_t) ? findall(in(Symbol.(PD.values(PD.columns(df_t)))), axis) : []
    if length(indices_t) == 0
        AttrDesc(Val{:static}, attr, name, df, pnl, idx, axis)
    else
        @assert(length(indices_t) == length(axis),
                "$name has static and time-dependent $(attr[:PyPSA])")
        AttrDesc{typenames[attr[:dtype]]}(df_t, indices_t, attr[:default])
    end
end

attributedescriptions(name, df, pnl, idx, axis; attributes=nothing) =
    Dict(attr[:attribute] => AttrDesc(Val{attr[:dimensions]}, attr, name, df, pnl, idx, axis)
         for attr in eachrow(attributes)
         if attr[:PyPSA] in PD.values(PD.columns(df)))

function typeparamdescriptions(name, df, indices; types=nothing)
    isnothing(types) && return nothing
    typename = @consense(PD.values(df.type)[indices], "$name may not have more than one type")
    typename == "" && return nothing
    typ = types[types.name .== typename, :]
    Dict(c => typ[1, c] for c in names(typ) if c != :name)
end

"""
    PypsaNetworkData

"""
function PypsaNetworkData(net)
    # Determine defined components and split them into classes
    components = vcat(
        getcomponents(net, :buses),
        getcomponents(net, :generators, withname=false, exp=:p_nom_extendable, uc=:committable),
        getcomponents(net, :storageunits, withname=false, pypsaname=:storage_units, exp=:p_nom_extendable),
        getcomponents(net, :stores, carrier=:bus, exp=:e_nom_extendable),
        getcomponents(net, :loads, carrier=:bus),
        getcomponents(net, :links, carrier=(:bus0, :bus1), exp=:p_nom_extendable),
        getcomponents(net, :lines, carrier=:bus0, exp=:s_nom_extendable, types=:line_types),
        getcomponents(net, :transformers, carrier=:bus0, exp=:s_nom_extendable, types=:transformer_types)
    )

    compdict = Dict(cd.name => cd for cd in components)

    PypsaNetworkData(net, compdict)
end

function _maybe_as_range(x::AbstractArray{T,1}) where T<:Number
    if length(x) < 3 return Vector(x) end
    d = diff(x)
    s = d[1]
    if all(d .== s)
        s == 1 ? (x[1]:x[end]) : (x[1]:s:x[end])
    else
        Vector(x)
    end
end

_join(a...; delim="_") = join(a, delim)

function getcomponents(net, name; pypsaname=false, withname=true, carrier=true, exp=false, uc=false, types=nothing)
    if pypsaname === false
        pypsaname = name
    end

    df = PD.DataFrame(getproperty(net, pypsaname))
    pnl = Dict(k=>PD.DataFrame(v) for (k,v) in getproperty(net, Symbol(pypsaname, :_t)))

    components = ComponentDesc[]

    if isempty(df)
        return components
    end

    initialname = withname ? (name,) : ()

    attrs = DataFrame()
    if carrier === true
        attrs[!, :carrier] = PD.values(df.carrier)
    elseif carrier !== false
        buses_carrier = PD.Series(net.buses.carrier)
        attrs[!, :carrier] =
            if carrier isa Tuple
                _join.((map(buses_carrier, df[field]) for field in carrier)...)
            else
                PD.values(map(buses_carrier, df[carrier]))
            end
    end

    if exp !== false
        attrs[!, exp] = recode(Bool.(PD.values(df[exp])), false=>:fix, true=>:ext)
    end
    if uc !== false
        attrs[!, uc] = recode(Bool.(PD.values(df[uc])), false=>:lin, true=>:uc)
    end

    if !isnothing(types)
        types = DataFrame(PD.reset_index(PD.DataFrame(getproperty(net, types))))
    end

    componenttypeunion = resolve(Component, name)
    componentattributes = copy(attributes(componenttypeunion))
    componentattributes = componentattributes[componentattributes.quantitytype .!= "Variable", :]
    componentattributes[!, :PyPSA] = string.(componentattributes.attribute)
    componentattributes[!, :dimensions] = recode(length.(componentattributes.dimensions),
                                                 0=>:single, 1=>:static, 2=>:series,
                                                 3:20=>missing)

    names = PD.values(PD.index(df))
    for (timeattr, df) in pnl
        if !isempty(df)
            label = Symbol(:t_, timeattr)
            attrs[!, label] = recode(in.(names, Ref(PD.values(PD.columns(df)))),
                                     false=>missing, true=>label)
        end
    end

    groups =
        if isempty(attrs)
            [(idx=1:length(names),
              name=Symbol(join(skipmissing(initialname), "_")),
              attrs=Dict{Symbol}{Any}())]
        else
            [(idx=_maybe_as_range(group.idx),
              name=Symbol(join(skipmissing(group.name), "_")),
              attrs=Dict(pairs(attrs[first(group.idx),:])...))
             for group in splitgroup(attrs, (idx=1:size(attrs, 1), name=initialname), Base.names(attrs)...)]
        end

    for g in groups
        axis = Axis{g.name}(Symbol.(names[g.idx]))

        formulation = join(skipmissing([
            exp !== false ? ((coalesce(g.attrs[exp], :fix) == :ext) ? "linexp" : missing) : missing,
            uc !== false ?  ((coalesce(g.attrs[uc], :lin) == :uc)   ? "ucdisp" : "lindisp") : "lindisp"
        ]), "_")

        componenttype =
            if componenttypeunion <: Device && !isempty(formulation)
                componenttypeunion{resolve(DeviceFormulation, Symbol(formulation))}
            else
                componenttypeunion
            end

        attributes = attributedescriptions(g.name, df, pnl, g.idx, axis; attributes=componentattributes)
        typeparams = typeparamdescriptions(g.name, df, g.idx; types=types)

        push!(components, ComponentDesc(g.name, axis, componenttype, attributes, typeparams))
    end

    components
end

_value(a::AxisArray) = a
_value(a::AxisArray{<:ConstraintRef}) = AxisArray(dual.(a), AxisArrays.axes(a)...)
_value(a::AxisArray{<:VariableRef}) = AxisArray(value.(a), AxisArrays.axes(a)...)
_value(a::AxisArray{<:GenericAffExpr}) = AxisArray(value.(a), AxisArrays.axes(a)...)

function _set!(df, args...)
    @assert(length(args) % 2 == 0, "Needs an even number of arguments")
    arrays = map(_value, args[2:2:end])
    axes = AxisArrays.axes(first(arrays))
    factors = [EM.WrappedArray(f, axes...) for f in args[1:2:end]]
    agg = AxisArray(Array{Float64}(undef, size(first(arrays))), axes)
    for I in eachindex(IndexCartesian(), agg)
        agg[I] = sum(f[Tuple(I)...] * v[I] for (f,v) in zip(factors, arrays))
    end
    _set!(df, agg)
end

function _set!(ser::PD.Series, data)
    ax, = AxisArrays.axes(data)
    PD.loc(ser)[ax.val] = _value(data)
end

function _set!(df::PD.DataFrame, data)
    B, T = AxisArrays.axes(data)
    for (idx, b) in zip(PD.columns(df).pyo.get_indexer(B), B)
        if idx == -1
            df[b] = PD.Series(_value(data[b, :]), index=T.val)
        else
            PD.loc(df)[T.val, idx] = _value(data[b, :])
        end
    end
end

function _set!(df::PyObject, data)
    pd = pyimport("pandas")
    if pyisinstance(df, pd.DataFrame)
        _set!(PD.DataFrame(df), data)
    elseif pyisinstance(df, pd.Series)
        _set!(PD.Series(df), data)
    else
        error("Needs pandas series or dataframe")
    end
end

function store_results!(data::PypsaNetworkData, g::EM.Generator)
    _set!(data.network.generators_t["p"], g[:p])
    _set!(data.network.generators.p_nom_opt, g[:p_nom])
end

function store_results!(data::PypsaNetworkData, s::EM.StorageUnit)
    _set!(data.network.storage_units_t["p"], 1, s[:p_dispatch], -1, s[:p_store])
    _set!(data.network.storage_units_t["state_of_charge"], s[:state_of_charge])
    _set!(data.network.storage_units_t["spill"], s[:spill])
    _set!(data.network.storage_units.p_nom_opt, s[:p_nom])
end

function store_results!(data::PypsaNetworkData, s::EM.Store)
    _set!(data.network.stores_t["p"], s[:p])
    _set!(data.network.stores_t["e"], s[:e])
    _set!(data.network.stores.e_nom_opt, s[:e_nom])
end

function store_results!(data::PypsaNetworkData, l::EM.Load)
    m = EM.model(l)
    T = axis(m, :snapshots)
    L = axis(m, l)
    # TODO: FIX
    # _set!(data.network.loads_t["p"], EM.get(l, :p_set, L, T))
end

function _pypsaname(::Type{T}) where T<:EM.Component
    S = T
    while !(S isa Symbol)
        S = S.name
    end
    string(S)
end

_mu_lower(b::EM.Branch{DF}) where DF <: EM.ExpansionForm = b[:p_lower]
_mu_lower(b::EM.Branch{DF}) where DF <: EM.DispatchForm = AxisArray(LowerBoundRef.(b[:p]), AxisArrays.axes(b[:p])...)
_mu_upper(b::EM.Branch{DF}) where DF <: EM.ExpansionForm = l[:p_upper]
_mu_upper(b::EM.Branch{DF}) where DF <: EM.DispatchForm = AxisArray(UpperBoundRef.(b[:p]), AxisArrays.axes(b[:p])...)

function store_results!(data::PypsaNetworkData, b::T) where T <: EM.PassiveBranch
    pnl = data.network.pnl(_pypsaname(T))
    _set!(pnl["p0"], b[:p])
    _set!(pnl["p1"], -1, b[:p])
    _set!(pnl["mu_lower"], _mu_lower(b))
    _set!(pnl["mu_upper"], -1, _mu_upper(b))
    _set!(data.network.df(_pypsaname(T)).s_nom_opt, b[:s_nom])
end

function store_results!(data::PypsaNetworkData, l::EM.Link)
    pnl = data.network.links_t
    _set!(pnl["p0"], l[:p])
    _set!(pnl["p1"], AxisArray(- l[:efficiency], AxisArrays.axes(l[:efficiency])...), l[:p])
    _set!(pnl["mu_lower"], _mu_lower(l))
    _set!(pnl["mu_upper"], -1, _mu_upper(l))
end

function store_results!(data::PypsaNetworkData, b::EM.Bus)
    pnl = data.network.buses_t
    _set!(pnl["marginal_price"], b[:p_balance])
    _set!(pnl["p"], EM.expression(EM.NodalActivePower(), EM.model(b), b))
    # TODO: v_ang and v_mag_pu needs to be calculated from B
end

function store_results!(data::PypsaNetworkData, c)
    @warn "Storing nothing for $(naming(c))"
end

function store_results!(data::PypsaNetworkData, model::EM.EnergyModel)
    map(b->store_results!(data, b), values(model.buses))
    for d in devices(model)
        store_results!(data, d)
    end
end

Base.get(::PypsaNetworkData, attrdesc::SingleAttrDesc, _) = attrdesc.value

Base.get(data::PypsaNetworkData, attrdesc::AttrDesc{Symbol,S}, ax) where S<:PD.Series =
    AxisArray(Symbol.(PD.values(attrdesc.data)[attrdesc.indices]), ax)

Base.get(data::PypsaNetworkData, attrdesc::AttrDesc{T,S}, ax) where {T, S<:PD.Series} =
    AxisArray(PD.values(attrdesc.data)[attrdesc.indices], ax)

Base.get(data::PypsaNetworkData, attrdesc::AttrDesc{T,S}, ax) where {T, S<:PD.DataFrame} =
    AxisArray(attrdesc.data.pyo.values[:,attrdesc.indices], data.axes[:snapshots], ax)

function Base.get(data::PypsaNetworkData, component::Component, param::Symbol)
    componentdesc = data.components[naming(component)]
    attrdesc = get(componentdesc.attributes, param, nothing)
    !isnothing(attrdesc) ? get(data, componentdesc.attributes[param], componentdesc.axis) : nothing
end

EM.gettypeparams(data::PypsaNetworkData, device::Device) =
    data.components[naming(device)].typeparams

EM.axis(data::PypsaNetworkData, c::Component) = data.axes[naming(c)]
EM.axis(data::PypsaNetworkData, name::Symbol) = data.axes[name]

EM.components(data::PypsaNetworkData) = EM._components(data)

end

max_stars_repo_name: coroa/EnergyModels

max_issues_repo_name: coroa/EnergyModels

max_forks_repo_name: coroa/EnergyModels


text:
function mcmcmodel(panel::AbstractDataFrame,  zs::ZSpec,::Val{:randomeffects},
    ::Type{T}=PARAM[:mcmctype],
    ::Type{TV}=PARAM[:mcmcgpu] ? CUDA.CuVector{T} : Vector{T},
    ::Type{TM}=PARAM[:mcmcgpu] ? CUDA.CuMatrix{T} : Matrix{T},
    prior::Symbol=PARAM[:mcmcprior],
    Nburn::Int=PARAM[:mcmcnburn],
    Ndraws::Int=PARAM[:mcmcndraws]
    ) where {
      T<:Real, TV<:AbstractVector{T}, TM<:AbstractMatrix{T}}
  @assert PARAM[:predictiontype] ∈ [:leveldemean, :level]
  (T == Float64) || @warn("T=$T, which may lead to numerical instability.
    Is this really what you want?")
  CUDA.allowscalar(false)

  Zms::MeasureSpec = zs.Zms
  Fw::Vector{Symbol} = Fweights(Zms, :Fw)
  FRLw::Vector{Symbol} = Fweights(Zms, :FRLw)
  essentialfields = [Fw; FRLw; Zms.Fvol; ]

  #form the data into a more usable format and expand the arrays
  p = modelparts(panel, zs,T,TV,TM,
    weightdata = (Fw=Fw, FRLw=FRLw),
    Fvol=Zms.Fvol)
  dims = p.dims
  tidx::Dict = p.tidx
  cleanpanel::SubDataFrame = p.cleanpanel
  Ncleanpanel::Int = nrow(cleanpanel)
  ws::TM = p.dat[:Fw]
  RLws::TM = p.dat[:FRLw]
  ts::Vector{Int} = p.ts
  Ascale::Vector{T} = p.Ascale

  v::TV = p.dat[:Fvol]

  #println("ts: $ts")
  Xv = XYIter(ws, RLws, v, ts, Ascale)

  #construct the main MCMC object for book keeping and parameters
  hyper=PARAM[prior]
  initial=(
    A₀ = ones(T, dims.K),
    G = ones(T, dims.K, dims.T-1),
    σ² = ones(T, dims.T-1),
    ω = zeros(T, dims.T-1)
    )
  cc = CapacityChains(T, Nburn=Nburn, Ndraws=Ndraws, paramstemplate=initial, hyper=hyper)
  #construct the parameter
  Θ = VolumePartsMCMC(TM, TV, T; initial, cc, Xv, ts)
  #printmln(Θ.Ψ)
  createmcmc(Θ)

  #for debugging purposes- will trhow an error at some point if true
  PARAM[:runstacktraceonmcmc] #=&& solveA!(iterloss,Θ,Xv,_Xv, SolveType)

  #main optimization routine
  try
    λ⁺ = solveA!(iterloss,Θ,Xv,_Xv, SolveType)
  catch err
    errmessage = "$err"
    if length(errmessage) > 20_000 #sometimes we get really long error messages
      errmessage = errmessage[1:20_000]
    end
    @warn "Solver failed with error $errmessage"
  finally
    results::DataFrame = formresults(panel, Zms, Θ, Val{:iter}())

    #recompute the loss if needed
    if λ⁺<0
      updateA₁!(Θ.g|>deepcopy, iterloss, Θ, Xv, _Xv, RegressionType, Val{:identity}())
      updateAfromGAₜ!(Θ,Xv, _Xv,1) #update to the final values of A
      λ⁺ = currentloss(Θ, Xv, _Xv, projectvolumefromA(Θ,Xv,_Xv))
    end

    if Ncleanpanel < 500_000 && PARAM[:testiterresults]
      Aexpanded = view(Θ.A', Θ.ts, :)
      LAexpanded = view(Θ.A', Θ.tsM1, :)
      need2demean = genabsμₖ.(Aexpanded, LAexpanded, Xv.ws, Xv.RLws)
      factored = Xv.M * need2demean

      Fcoefs::Vector{Symbol} = Fcoefficients0(zs.Zms)
      Fwcoefs::Vector{Symbol} = (s->Symbol(:w,s)).(Fcoefs)
      FRLwcoefs::Vector{Symbol} = (s->Symbol(:RLw,s)).(Fcoefs)
      FAcoefs::Vector{Symbol} = (s->Symbol(:A,s)).(Fcoefs)
      FLAcoefs::Vector{Symbol} = (s->Symbol(:LA,s)).(Fcoefs)
      Ffcoefs::Vector{Symbol} = (s->Symbol(:f,s)).(Fcoefs)

      dfdebug = hcat(
        DataFrame(ts=Θ.ts, tsM1=Θ.tsM1, v=Xv.v, vtilde=Xv.ṽ),
        DataFrame(Xv.ws, Fwcoefs),
        DataFrame(Xv.RLws, FRLwcoefs),
        DataFrame(Aexpanded, FAcoefs),
        DataFrame(LAexpanded, FLAcoefs),
        DataFrame(factored, Ffcoefs))

      dfdebug |> CSV.write(
        "$(PARAM[:testpath])\\$(PARAM[:crspdatalabel])_opt$(iterresultsname(λ⁺|>round|>Int)).csv")
    end

    return Int(round(min(λ⁺, typemax(Int)÷10))), results
  end=#

  @assert false #should never reach here due to the finalizer
end

function updateΨ!(Θ::AbstractVolumePartsMCMC{TM, TV, T}, ::Val{:randomeffects}) where {TM,TV,T}
  #@info "TM:$TM, TV:$TV, T:$T"
  updateA₀!(Θ)
end

#this will create the mcmc chain
function createmcmc(Θ::AbstractVolumePartsMCMC; stateupdater::Val = Val{PARAM[:mcmcstateupdater]}())
  @unpack cc, Ψ = Θ
  @unpack Nburn = cc

  Niter = Nburn + size(cc.value,1)

  for i ∈ 1:Niter
    (i==Nburn+1) && @assert(cc.idx==1)
    updateΨ!(Θ, stateupdater)
    (i<Niter) && updateandincrement!(cc, Ψ, verify=true)
  end
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
module Ogg

## FIXME within the module, remove the Ogg prefix in the type names
immutable IOVec
    iov_base::Ptr{Void}
    iov_len::Csize_t
end

immutable PackBuffer
    endbyte::Clong
    endbit::Cint

    buffer::Ptr{Cuchar}
    ptr::Ptr{Cuchar}
    storage::Clong
end

@doc "Encapsulate the data in one Ogg bitstream page"->
immutable Page
    header::Ptr{Cuchar}
    header_len::Clong
    body::Ptr{Cuchar}
    body_len::Clong
end

@doc """
Contains the current encode/decode state of a logical Ogg bitstream.
"""->
immutable StreamState
    body_data::Ptr{Cuchar}         # bytes from packet bodies
    body_storage::Clong            # elements allocated
    body_fill::Clong               # elements stored; fill mark
    body_returned::Clong           # elements of fill returned

    lacing_vals::Ptr{Cint}         # values to go in segment table
    granule_vals::Ptr{Int64}       # granulepos values for headers
    lacing_storage::Clong
    lacing_fill::Clong
    lacing_packet::Clong
    lacing_returned::Clong

    header::NTuple{282,Cuchar}     # working space for header encoding
    e_o_s::Cint                    # set on end of stream
    b_o_s::Cint                    # set after writing initial page

    serialno::Clong
    pageno::Clong
    packetno::Int64
    granulepos::Int64
end

@doc """
Encapsulate the data and metadata belonging to a single
raw Ogg/Vorbis packet
"""->
immutable Packet
    packet::Ptr{Cuchar}
    bytes::Clong
    b_o_s::Clong
    e_o_s::Clong

    granulepos::Int64

    packetno::Int64
end

immutable SyncState
    data::Ptr{Cuchar}
    storage::Cint
    fill::Cint
    returned::Cint

    unsynced::Cint
    headerbytes::Cint
    bodybytes::Cint
end

## Ogg Bitstream Primitives: bitstream
## I skipped the C functions beginning with oggpackB_ b/c they
## seem to be the same as the corresponding oggpack_ C function.

for (rtyp,nm) in ((:Void,:writeinit),
                  (:Cint,:writecheck),
                  (:Void,:writealign),
                  (:Void,:reset),
                  (:Void,:writeclear),
                  (:Clong,:look1),
                  (:Void,:adv1),
                  (:Clong,:read1),
                  (:Clong,:bytes),
                  (:Clong,:bits),
                  ("Ptr{Cuchar}",:get_buffer))
    @eval begin
        $nm(b::PackBuffer) =
            ccall(($(string("oggpack_",nm)),ogg),$rtyp,
                  (Ref{PackBuffer},),b)
    end
end

writetrunc(b::PackBuffer,bits) =
    ccall((:oggpack_writetrunc,ogg),Void,
          (Ref{PackBuffer},Clong),b,bits)

writecopy(b::PackBuffer,source,bits) =
    ccall((:oggpack_writecopy),Void,
          (Ref{Packbuffer},Ptr{Void},Clong),b,source,bits)

readinit(b::PackBuffer,buf,bytes) =
    ccall((:oggpack_readinit,ogg),Void,
          (Ref{PackBuffer},Ptr{Cuchar},Cint),b,buf,bytes)

Base.write(b::PackBuffer,value,bits) =
    ccall((:oggpack_write,ogg),Void,
          (Ref{PackBuffer},Culong,Cint),b,value,bits)

Base.read(b::PackBuffer,bits) =
    ccall((:oggpack_read,ogg),Culong,
          (Ref{PackBuffer},Cint),b,bits)

for (nm,rtyp) in ((:look,:Clong),(:adv,:Void))
    @eval begin
        $nm(b::PackBuffer,bits) =
            ccall(($(string("oggpack_",nm)),ogg),Clong,
                  (Ref{PackBuffer},Cint),b,bits)
    end
end

## Ogg Bitstream primitives: encoding

packetin(os::StreamState,op::Packet) =
    ccall((:ogg_stream_packetin,ogg),Cint,
          (Ref{StreamState},Ref{Packet}),os,op)

iovecin(os::StreamState,iov::IOVec,count,e_o_s,granulepos) =
    ccall((:ogg_stream_iovecin,ogg),Cint,
          (Ref{StreamState},Ref{IOVec},Cint,Clong,Int64),
          os,iov,count,e_o_s,granulepos)

for nm in (:pageout,:flush)
    @eval begin
        $nm(os::StreamState,og::Page) =
            ccall(($(string("ogg_stream_",nm)),ogg),Cint,
                  (Ref{StreamState},Ref{Page}),os,og)
        $nm(os::StreamState,og::Page,nfill) =
            ccall(($(string("ogg_stream_",nm,"_fill")),ogg),Cint,
                  (Ref{StreamState},Ref{Page},Cint),os,og,nfill)
    end
end

## Ogg Bitstream primitives: decoding

for nm in (:init,:clear,:reset,:destroy,:check)
    @eval begin
        $nm(oy::SyncState) =
            ccall(($(string("ogg_sync_",nm)),ogg),Cint,
                  (Ref{SyncState},),oy)
    end
end

buffer(oy::SyncState,sz) =
    ccall((:ogg_sync_buffer,ogg),Ptr{Cuchar},
          (Ref{SyncState},Clong),oy,sz)

wrote(oy::SyncState,bytes) =
    ccall((:ogg_sync_wrote,ogg),Cint,
          (Ref{SyncState},Clong),oy,bytes)

for (nm,rtyp) in ((:pageseek,:Clong),(:pageout,:Cint))
    @eval begin
        $nm(oy::SyncState,op::Page) =
            ccall(($(string("ogg_sync_",nm)),ogg),$rtyp,
                  (Ref{SyncState},Ref{Page}),oy,op)
    end
end

pagein(os::StreamState,og::Page) =
    ccall((:ogg_stream_pagein,ogg),Cint,
          (Ref{StreamState},Ref{Page}),os,og)

packetout(os::StreamState,op::Packet) =
    ccall((:ogg_stream_packetout,ogg),Cint,
          (Ref{StreamState},Ref{Packet}),os,op)

packetpeek(os::StreamState,op::Packet) =
    ccall((:ogg_stream_packetpeek,ogg),Cint,
          (Ref{StreamState},Ref{Packet}),os,op)

## Ogg Bitstream primitives: general

for nm in (:clear,:reset,:destroy,:check,:eos)
    @eval begin
        $nm(os::StreamState) =
            ccall(($(string("ogg_stream_",nm)),ogg),Cint,
                  (Ref{StreamState},),os)
    end
end

for nm in (:init,:reset_serialno)
    @eval begin
        $nm(os::StreamState,sn) =
            ccall(($(string("ogg_stream_",nm)),ogg),Cint,
                  (Ref{StreamState},Cint),os,sn)
    end
end

for (nm,rtyp) in ((:checksum_set,:Void),
                  (:version,:Cint),
                  (:bos,:Cint),
                  (:eos,:Cint),
                  (:granulepos,:Int64),
                  (:serialno,:Cint),
                  (:pageno,:Clong),
                  (:packets,:Cint))
    @eval begin
        $nm(og::Page) =
            ccall(($(string("ogg_page_",nm)),ogg),$rtyp,
                  (Ref{Page},),og)
    end
end

clear(op::Packet) =
    ccall((:ogg_packet_clear,ogg),Void,(Ref{Packet},),op)

end


text:
#import Base: size, getindex, setindex!, IndexStyle, randn, similar

#using .Domain
#include("Domain.jl")

################################################################################
################################################################################
#Field type
struct Field{T, N, Dom} <: AbstractArray{T, N} where Dom <: AbstractDomain
    domain::Dom
    val::Array{T, N}
    Field(domain, val) = begin
        @assert size(domain) == size(val)
        new{eltype(val), dims(domain), typeof(domain)}(domain, val)
    end
end

struct MultiField
    domain::Union{Tuple, NamedTuple}
    val::Union{Tuple, NamedTuple}
end

function MultiField(fs::Union{Field, MultiField}...)
    domains = Tuple(f.domain for f in fs)
    vals = Tuple(f.val for f in fs)
    return MultiField(domains, vals)
end

#Implement correct broadcasting behaviour
#Pointwise operations on a field should always default to preserving the domain
#Binary broadcasting is only well defined for fields with matching domains
struct FieldStyle{N, Dom} <: Broadcast.AbstractArrayStyle{N} end
Base.BroadcastStyle(::Type{<:Field{T, N, Dom}}) where {T, N, Dom} = FieldStyle{N, Dom}()
Base.BroadcastStyle(::FieldStyle{N, Dom}, ::FieldStyle{N, Dom}) where {N, Dom} = FieldStyle{N, Dom}()
Base.BroadcastStyle(::FieldStyle, ::FieldStyle) = error("Domain mismatch")
#Base.BroadcastStyle(::Broadcast.Style{<:Number}, ::FieldStyle{N, Dom}) where {N, Dom} = FieldStyle{N, Dom}()
#Base.BroadcastStyle(::FieldStyle{N, Dom}, ::Broadcast.Style{<:Number}) where {N, Dom} = FieldStyle{N, Dom}()

similar(f::Field) = similar(f, eltype(f), size(f))
similar(f::Field, ::Type{T}, ::Tuple{Vararg{Int64, N}}) where {T, N} = Field(f.domain, similar(f.val, T))
function similar(bc::Broadcast.Broadcasted{FieldStyle{N, Dom}}, ::Type{ElType}) where {N, Dom, ElType}
    domain = getdomain(bc)
    return Field(domain, similar(Array{ElType}, axes(bc)))
end

getdomain(bc::Base.Broadcast.Broadcasted) = getdomain(bc.args...)
#TODO assure, that all args have the same domain?
getdomain(args::Vararg{Field}) = getdomain(args[1])
getdomain(f::Field) = f.domain

shape(f::Field) = shape(f.domain)
size(f::Field) = Int.(size(f.val))
IndexStyle(::Type{<:Field}) = IndexLinear()

getindex(f::Field, i::Integer) = getindex(f.val, i)
getindex(f::Field, I::Vararg{Integer, N}) where N = getindex(f.val, I...)
setindex!(f::Field, v, i::Integer) = setindex!(f.val, v, i)
setindex!(f::Field, v, I::Vararg{Integer, N}) where N = setindex!(f.val, v, I...)

randn(d::AbstractDomain) = Field(d, randn(size(d)...))
randn(::Type{T}, d::AbstractDomain) where T = Field(d, randn(T, size(d)...))
randn(::Type{T}, d::AbstractDomain) where T <: Integer = Field(d, round.(T, randn(size(d)...)))

pindices(ps::Field{T, N, PD}) where {T, N, PD <: PowerDomain} = ps.domain._pindices

#MultiField: named tuples of fields


text:
### A Pluto.jl notebook ###
# v0.16.3

using Markdown
using InteractiveUtils

# ╔═╡ e56ad662-0c01-48ca-addb-615e9350c6d9
1

# ╔═╡ 43e5dfe4-90b2-4f7b-b467-3107b78520f4
x = 1 + 1

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.0-rc1"
manifest_format = "2.0"

[deps]
"""

# ╔═╡ Cell order:
# ╠═e56ad662-0c01-48ca-addb-615e9350c6d9
# ╠═43e5dfe4-90b2-4f7b-b467-3107b78520f4
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
@testset "869.reordered-power-of-2.jl" begin
    @test reordered_power_of2(1)
    @test !reordered_power_of2(10)
    @test !reordered_power_of2(24)
    @test reordered_power_of2(46)
    @test reordered_power_of2(16)
end


text:
# Use typeof to find the types of types. Try Symbol or Int64. Can you find the types of operators?
typeof(Symbol)
typeof(Int64)
typeof(+)

# Create a typed dictionary with keys that are symbols and values that are floats. What happens when you add :thisis => :notanumber to the dictionary?
dict = Dict(:a => 1.1)
dict[:thisis] = :notanumber

# Create a 5×5×5 array where each 5×5 block in the first two dimensions is a single number but that number increases for each block. For example, magic[:,:,1] would have all elements equal to 1, and magic[:,:,2] would have all elements equal to 2.
magic = [z for x in 1:5, y in 1:5, z in 1:5]

# Run some arrays of various types through functions like sin and round. What happens?
sin(fill(1.5,(3,3)))
round(fill(1.5, (3,3)))
sin(fill(:a, (3,3))
round(fill(:a, (3,3)))

# Create a matrix and multiply it by its inverse. Hint: inv computes the inverse of a matrix, but not all matrices are invertable.
m = [1 2; 3 4]
m * inv(m)

# Create two dictionaries and merge them. Hint: Look up merge in the manual.
d1 = Dict(:a => 1)
d2 = Dict(:b => 2)
merge(d1, d2)

# sort and sort! both operate on arrays. What is the difference between them?
a = [2, 3, 4, 1, 5]
sort(a)
a
sort!(a)

max_stars_repo_name: hemslo/seven-in-seven

max_issues_repo_name: hemslo/seven-in-seven

max_forks_repo_name: hemslo/seven-in-seven


text:
using LogicCircuits
using LogicCircuits: JuiceTransformer, dimacs_comments, zoo_version

using LazyArtifacts
using Lerche: Lerche, Lark, Transformer, @rule, @inline_rule

include("circuit_io.jl")

# if no logic circuit file format is given on read, infer file format from extension
function file2pcformat(file)
    if endswith(file,".circuit")
        DiscriminativeFormat()
    else
        # try a logic circuit format
        LogicCircuits.file2logicformat(file)
    end
end

"""
    Base.read(file::AbstractString, ::Type{C}) where C <: LogisticCircuit

Reads circuit from file; uses extension to detect format type.
"""
Base.read(file::AbstractString, ::Type{C}) where C <: LogisticCircuit =
    read(file, C, file2pcformat(file))

Base.read(files::Tuple{AbstractString,AbstractString}, ::Type{C}) where C <: LogisticCircuit =
    read(files, C, (file2pcformat(files[1]), VtreeFormat()))

"""
    Base.write(file::AbstractString, circuit::LogisticCircuit)

Writes circuit to file; uses file name extention to detect file format.
"""
Base.write(file::AbstractString, circuit::LogisticCircuit) =
    write(file, circuit, file2pcformat(file))

"""
    Base.write(files::Tuple{AbstractString,AbstractString}, circuit::LogisticCircuit)

Saves circuit and vtree to file.
"""
Base.write(files::Tuple{AbstractString,AbstractString},
           circuit::LogisticCircuit) =
    write(files, circuit, (file2pcformat(files[1]), VtreeFormat()))

# # copy read/write API for tuples of files
function Base.read(files::Tuple{AbstractString, AbstractString}, ::Type{C}, args...) where C <: LogisticCircuit
    open(files[1]) do io1
        open(files[2]) do io2
            read((io1, io2), C, args...)
        end
    end
end

function Base.write(files::Tuple{AbstractString,AbstractString},
                    circuit::LogisticCircuit, args...)
    open(files[1], "w") do io1
        open(files[2], "w") do io2
            write((io1, io2), circuit, args...)
        end
    end
end

max_stars_repo_name: Juice-jl/DiscriminativeCircuits

max_issues_repo_name: Juice-jl/DiscriminativeCircuits

max_forks_repo_name: Juice-jl/DiscriminativeCircuits


text:
#=
Variations on Filon-Clenshaw-Curtis quadrature.
=#
module FCCQuad
export fccquad, fccquad_cc, fccquad_cs, fccquad_sc, fccquad_ss

include("Fct.jl")
include("Cheb.jl")
include("chirps.jl")
include("chebweights.jl")
include("Jets.jl")

using LinearAlgebra, Scratch

const AV = AbstractVector
const AA = AbstractArray

#=
Filon-Clenshaw-Curtis quadrature for f(x) on [-1,1].
Assuming chebfun is the coefficient vector of
a Chebyshev expansion of degree at least N of some function f(x),
returns f(x)*exp(im*w*x)*dx integrated over [-1,1]
using degree N, and discrepancy between that and using reduced_degree(N).
WARNING: N and baseN must be even. An @assert enforces.
=#
function fccquadUnit!(chebfun::AV,freq::Real,N::Integer,
                      weightmethod::Symbol,weights::AV,weights_workspace::AA,
                      baseN=reduced_degree(N))
    getweights!(weights,weights_workspace,length(weights)-1,freq,weightmethod)
    fcc_core(chebfun,weights,N,baseN)
end
reduced_degree(N)=div(3N,4)
function fcc_core(chebfun::AV,weights::Vector,N::Integer,baseN::Integer)
    @assert iseven(N) && iseven(baseN)
    a=chebfun
    i=weights
    base=(sum(a[n]*i[n] for n in baseN+1:-2:1) +
          sum(a[n]*i[n] for n in baseN:-2:2)*im)
    diff=(sum(a[n]*i[n] for n in N+1:-2:baseN+3) +
          sum(a[n]*i[n] for n in N:-2:baseN+2)*im)
    full = base + diff
    full,diff
end
#=
Alternative, more flexible error estimation method
that compares the results of using two user-provided Chebyshev expansions
instead of one Chebyshev expansion and its truncation.
WARNING: N and baseN must be even. An @assert enforces.
=#
function fcc_alt(chebfun::AV,basefun::AV,
                 weights::Vector,N::Integer,baseN::Integer)
    @assert iseven(N) && iseven(baseN)
    a,b,i=chebfun,basefun,weights
    full=(sum(a[n]*i[n] for n in N+1:-2:1) +
          sum(a[n]*i[n] for n in N:-2:2)*im)
    base=(sum(b[n]*i[n] for n in baseN+1:-2:1) +
          sum(b[n]*i[n] for n in baseN:-2:2)*im)
    diff = full - base
    full,diff
end

#=
Filon-Clenshaw-Curtis quadrature for f(x) with arbitrary finite interval domain.
Assuming chebfun is the coefficient vector of
a Chebyshev expansion of degree at least N of some function g(x)=f(c+rx),
returns the integral over [-1,1] of f(c+rx)*exp(im*(w*(c+rx)+v*rx))*d(c+rx)
using degree N, and discrepancy between that and using degree baseN.
Above, w=freq, c=center, r=radius, and v=doppler.
WARNING: N and baseN must be even. An @assert enforces.
=#
function fccquadSampled!(chebfun::AV,freq::Real,N::Integer,
                         center::Real,radius::Real,doppler::Real,
                         weightmethod::Symbol,weights::AV,weights_workspace::AA,
                         baseN::Integer=reduced_degree(N))
    w = freq*radius + doppler
    interval_transform(fccquadUnit!(chebfun,w,N,weightmethod,weights,weights_workspace,baseN)...,
                       freq,center,radius)
end
function interval_transform(full,diff,freq,center,radius)
    si,co=sincos(freq*center)
    mult = radius*complex(co,si)
    full*mult, diff*mult
end

#Integrates f1(x)*f2(x)*exp(im*w*x)*dx over [xmin,xmax] for each w in freqs
#using degree N and discrepancy between that and using reduced_degree(N).
#WARNING: log2N must be at least 3. An @assert enforces.
function fccquadBatch(f1::Function,f2::Function,freqs::AA,log2N::Integer;
                      weightmethod::Symbol=:thomas,T::Type=Complex{Float64},
                      xmin::Real=-one(real(T)),xmax::Real=one(real(T)))
    @assert 3 <= log2N <= 62
    xminT,xmaxT = real(T)(xmin),real(T)(xmax)
    center = real(T)(0.5)*(xmaxT + xminT)
    radius = real(T)(0.5)*(xmaxT - xminT)
    output,workspaces = fccquad_alloc(freqs,log2N,weightmethod,T)
    fccquadBatch!(output,workspaces,center,radius,f1,f2,freqs,log2N,weightmethod)
    output,1+(1<<log2N)
end
function fccquadBatch!(output::AA,workspaces,center::Real,radius::Real,
                       f1::Function,f2::Function,freqs::AA,log2N::Integer,
                       weightmethod::Symbol,baseN=reduced_degree(1<<log2N))
    N=1<<log2N
    samples,fct_workspaces = workspaces[1]
    a=Fct.chebcoeffs!(Fct.shiftedchebsample!(center,radius,f1,f2,samples,N),fct_workspaces...)
    M=length(freqs)
    weights,weights_workspace = workspaces[2]
    for m in 1:M
        output[:,m] = collect(fccquadSampled!(a,freqs[m],N,center,radius,zero(real(eltype(output))),
                                              weightmethod,weights,weights_workspace,baseN))
    end
    output
end
function fccquad_alloc(freqs::AA,log2N::Integer,weightmethod::Symbol,T::Type=Complex{Float64})
    output=Matrix{T}(undef,2,length(freqs))
    N=1<<log2N
    output,fccquad_workspaces(N,weightmethod,T)
end
function fccquad_workspaces(N::Integer,weightmethod::Symbol,T::Type=Complex{Float64})
    samples = Vector{T}(undef,1+N)
    (samples,Fct.fct_alloc(samples,T)),weights_alloc(N,weightmethod,real(T))[1:2]
end

#=
Degree-adaptive FCC quadrature implementation.
Integrates f1(c+rx)*f2(c+rx)*exp(im*(w*(c+rx)+v*rx))*d(c+rx) over [-1,1]
for each w in freqs, where c=center, r=radius, and v=doppler.

Expects output,workspaces = fccquad_alloc(freqs,maxlog2degree,weightmethod,T)
OR workspaces==nothing and size(output)==(2,length(freqs)).

WARNING: log2N is assumed to be at least 3.
=#
function adaptdegree!(f1::Function,f2::Function,freqs::AA,
                      center::Real,radius::Real,doppler::Real,
                      weightmethod::Symbol,output::AA,workspaces,
                      minlog2degree::Integer,maxlog2degree::Integer,
                      reltol::Real,abstol::Real,vectornorm::Function)
    N = 1 << minlog2degree
    T = eltype(output)
    if workspaces == nothing
        samples = Fct.shiftedchebsample(center,radius,f1,f2,N;T=T)
    else
        samples,fct_workspaces = workspaces[1]
        ifft_in,ifft_out,ifft_work = fct_workspaces
        weights_workspaces = workspaces[2]
        Fct.shiftedchebsample!(center,radius,f1,f2,samples,N)
    end

    R = real(T)
    base = zero(R)
    success = false
    maxN = 1 << maxlog2degree
    while true
        if workspaces == nothing
            ifft_out=Fct.chebcoeffs(samples)
        else
            Fct.chebcoeffs!(samples,ifft_in,ifft_out,ifft_work,N)
        end
        for m in 1:length(freqs)
            if workspaces == nothing
                weights_workspaces = weights_alloc(N,weightmethod,R)[1:2]
            end
            output[:,m] = collect(fccquadSampled!(ifft_out,freqs[m],N,center,radius,doppler,
                                                  weightmethod,weights_workspaces...))
        end
        base=vectornorm(view(output,1,:))
        delta=vectornorm(view(output,2,:))
        success = delta <= base * reltol || delta <= abstol
        if success || N >= maxN
            break
        end
        if workspaces == nothing
            samples = Fct.shifteddoublesample(center,radius,f1,f2,samples)
        else
            Fct.shifteddoublesample!(center,radius,f1,f2,samples,N)
        end
        N<<=1
    end
    success,N+1,base
end

#=
Degree-adaptive FCC quadrature with tone removal.
Expects output,workspaces = fccquad_alloc(freqs,maxlog2N,weightmethod,T).
Integrates prefactor(x)*oscillator(x)*exp(im*w*x)*dx over [c-r,c+r]
for each w in freqs, where c=center and r=radius.
WARNING: arg(oscillator(x)) is assumed to be smooth at x=c.
WARNING: log2N is assumed to be at least 3.
=#
function tonequad!(output::AA,workspaces,center::Real,radius::Real,
                   prefactor::Function,oscillator::Function,freqs::AA,
                   minlog2N::Integer,maxlog2N::Integer,weightmethod::Symbol,
                   T::Type,reltol::Real,abstol::Real,vectornorm::Function)
    cfreq = Jets.phase_velocity(oscillator(Jets.Jet(center,radius)))
    if !isfinite(cfreq)
        cfreq = zero(cfreq)
    end

    samples,fct_workspaces = workspaces[1]
    ifft_in,ifft_out,ifft_work = fct_workspaces
    weights,weights_work = workspaces[2]

    N = 1 << minlog2N
    maxN = 1 << maxlog2N
    Fct.tonechebsample!(cfreq,center,radius,prefactor,oscillator,samples,N)
    success = false
    base = 0
    while true
        Fct.chebcoeffs!(samples,ifft_in,ifft_out,ifft_work,N)
        for m in 1:length(freqs)
            output[:,m] = collect(fccquadSampled!(ifft_out,freqs[m],N,center,radius,cfreq,
                                                  weightmethod,weights,weights_work))
        end
        base=vectornorm(view(output,1,:))
        delta=vectornorm(view(output,2,:))
        success = delta <= base * reltol || delta <= abstol
        if success || N >= maxN
            break
        end
        Fct.tonedoublesample!(cfreq,center,radius,prefactor,oscillator,samples,N)
        N<<=1
    end
    success,N+1,base
end

#=
Degree-adaptive FCC quadrature with chirp removal.
Integrates prefactor(x)*oscillator(x)*exp(im*w*x)*dx over [c-r,c+r]
for each w in freqs, where c=center and r=radius.
WARNING: arg(oscillator(x)) is assumed to be smooth at x=c.
WARNING: log2N is assumed to be at least 3.
=#
function chirpquad!(failfast::Bool,output::AA,center::Real,radius::Real,
                    prefactor::Function,oscillator::Function,freqs::AA,
                    minlog2N::Integer,maxlog2N::Integer,weightmethod::Symbol,
                    T::Type,reltol::Real,abstol::Real,vectornorm::Function,
                    maxchirp::Integer)
    N = 1 << minlog2N
    maxN = 1 << maxlog2N
    R = real(T)

    jet = oscillator(Jets.Jet(center,radius))
    cfreq::Real = Jets.phase_velocity(jet)
    if !isfinite(cfreq)
        cfreq = zero(cfreq)
    end
    chirp::Real = Jets.phase_acceleration(jet)
    if !isfinite(chirp)
        chirp = zero(chirp)
    end
    if Chirps.rates[maxchirp] < abs(chirp)
        failfast && return false,1,zero(R)
        workspaces=fccquad_workspaces(maxN,weightmethod,T)
        return tonequad!(output,workspaces,center,radius,
                         prefactor,oscillator,freqs,
                         minlog2N,maxlog2N,weightmethod,
                         T,reltol,abstol,vectornorm)
    end

    samples = Fct.chirpchebsample(cfreq,chirp,center,radius,prefactor,oscillator,N;T=T)
    success = false
    base = zero(R)
    while true
        a=Cheb.ChebSeries(Fct.chebcoeffs(samples))

        index=findfirst(x->abs(chirp)<=x,view(Chirps.rates,1:maxchirp))::Integer
        faster=Chirps.rates[index]
        b0 = Cheb.EvenChebSeries(Chirps.getchirp(index))
        b=Cheb.dilate(b0,sqrt(abs(chirp)/faster))
        if chirp < zero(chirp)
            b = conj(b)
        end
        ab,ab2 = a*b, Cheb.truncate(a,reduced_degree(N))*b
        deg = Cheb.degree(Cheb.autotruncate(ab))
        deg2 = Cheb.degree(Cheb.autotruncate(ab2))
        @assert deg >=2 && deg2 >= 2
        deg1 = max(1, deg >> 2) << 2 #truncate to multiple of 4
        deg2 = max(1, deg2 >> 2) << 2 #truncate to multiple of 4
        deg3 = max(deg,deg2)

        weights,workspace,deg4=weights_alloc(deg3,weightmethod,R)
        for m in 1:length(freqs)
            w = freqs[m]*radius + cfreq
            getweights!(weights,workspace,deg4,w,weightmethod)
            step1 = fcc_alt(ab.coeff,ab2.coeff,weights,deg1,deg2)
            step2 = interval_transform(step1...,freqs[m],center,radius)
            output[:,m] = collect(step2)
        end

        base=vectornorm(view(output,1,:))
        delta=vectornorm(view(output,2,:))
        success = delta <= base * reltol || delta <= abstol
        if success || N >= maxN
            break
        end
        samples = Fct.chirpdoublesample(cfreq,chirp,center,radius,prefactor,oscillator,samples)
        N<<=1
    end
    success,N+1,base
end

#constants for fccquad()
const interval_methods = (:plain, :tone, :chirp)
const all_methods = (interval_methods..., :degree, :nonadaptive)
const supported_types = Dict(
    (m=>(Complex{T} for T in (Float32,Float64,BigFloat)) for m in (:plain, :tone))...,
    (m=>(Complex{T} for T in (Float64,BigFloat)) for m in (:degree, :nonadaptive))...,
    :chirp=>(Complex{Float64},))
const default_reltol = Dict(Float32=>1e-6,Float64=>1e-8,BigFloat=>1e-39)
const weightmethod=:thomas

"""
    ```
    fccquad(prefactor::Function,
            oscillator::Function,
            freqs::AbstractArray{<:Real};
            kwargs...)
    ```

Integrates `prefactor(x)*oscillator(x)*exp(im*w*x)*dx`
over a finite interval, for all `w` in `freqs`,
using a variant of Filon-Chenshaw-Curtis (FCC) quadrature.

Use the `xmin` and `xmax` keyword arguments to specify
a finite domain of integration [xmin,xmax]. Default: [-1,1].

Required: `oscillator` should be nonzero on [xmin,xmax].

Recommended: `prefactor` should not be highly oscillatory on [xmin,xmax].
The amplitude of `oscillator` and phase velocity of `oscillator`
should also not be highly oscillatory.

Outputs a tuple `(results, function_evaluation_count)` where `results`
is a 2xL complex matrix, where `L=length(freqs)`,
with row 1 containing the estimated integrals
and row 2 the discrepancies between row 1 and estimates made
using a reduced Chebyshev polynomial interpolation degree.

The `reltol` and `abstol` keyword arguments specify
the relative and absolute error goals used by the adaptive methods.
Success means satisfying both goals.
Defaults: `abstol=0` and `reltol=default_reltol[real(T)]`.

The `method` keyword argument specifies the FCC variant to be used.
`FCQuad.all_methods` lists the supported methods.
The default is `:tone`.

* Method `:nonadaptive` is FCC quadrature with the base-2 logarithm of
  the interpolation degree specified by `nonadaptivelog2degree` (default: 10).

* Method `:degree` is a degree-adaptive FCC quadrature with the base-2 logarithm of
  the interpolation degree starting at `minlog2degree` (default: 3)
  and adaptively increasing but not beyond `globalmaxlog2degree` (default: 20).

* Methods `:plain`, `:tone`, and `:chirp` are
  hybrid interval-adaptive degree-adaptive FCC quadrature variants
  with the base-2 logarithm of the per-subinterval interpolation degree
  starting at `minlog2degree` (default: 3) and adaptively increasing
  but not beyond `localmaxlog2degree` (default: 6).
  If a subinterval's accuracy goals are not met even with the maximum degree,
  then the subinterval is divided into `branching`-many equal subintervals.
  (default: 4).

* Methods `:tone` and `:chirp` use automatic differentiation
  to factor out a per-subinterval tone or chirp (respectively) from `oscillator`
  before the polynomial interpolation step. To use these methods,
  the function `oscillator` must be generic enough to take input
  of type `FCCQuad.Jets.Jet{real(T)}`, which is a subtype of `Number`.

The `T` keyword argument specifies the working complex scalar datatype,.
(Default: `Complex{Float64}`.)
`FCQuad.supported_types` lists which datatypes are supported
for which methods.
The functions `prefactor` and `oscillator` should output
something of type `T` or `real(T)` when the input has type `real(T)`.

The adaptive methods use stopping rules based on vectors of
discrepancies involving all frequencies in `freqs`.
The `vectornorm` keyword argument specifies the norm used
to compare vectors.
Default: `LinearAlgebra.norm` (Euclidean norm).

See also: `fccquad_cc`, `fccquad_cs`, `fccquad_sc`, `fccquad_ss`.
"""
function fccquad(prefactor::Function,oscillator::Function,freqs::AA{<:Real};
                 T::Type=Complex{Float64},xmin::Real=-one(real(T)),xmax::Real=one(real(T)),
                 reltol::Real=default_reltol[real(T)],abstol::Real=zero(real(T)),
                 method::Symbol=:tone,vectornorm=LinearAlgebra.norm,
                 branching::Integer=4,maxdepth::Integer=10,maxchirp::Integer=6,
                 minlog2degree::Integer=3,
                 localmaxlog2degree::Integer=6,
                 nonadaptivelog2degree::Integer=10,
                 globalmaxlog2degree::Integer=20)
    @assert xmin < xmax && isfinite(xmin) && isfinite(xmax)
    @assert 3 <= minlog2degree <= localmaxlog2degree <= 62
    @assert 3 <= minlog2degree <= globalmaxlog2degree <= 62
    @assert 3 <= nonadaptivelog2degree <= 62
    @assert method in all_methods
    @assert T in supported_types[method]
    @assert 2 <= branching
    @assert 1 <= maxchirp <= length(Chirps.rates)

    if method == :nonadaptive
        return fccquadBatch(prefactor,oscillator,freqs,nonadaptivelog2degree;
                            T=T,xmin=xmin,xmax=xmax,weightmethod=weightmethod)
    end

    xminT,xmaxT = real(T)(xmin),real(T)(xmax)
    center = real(T)(0.5)*(xmaxT + xminT)
    radius = real(T)(0.5)*(xmaxT - xminT)
    if method == :degree
        output=Matrix{T}(undef,2,length(freqs))
        success,evals,base=adaptdegree!(prefactor,oscillator,freqs,
                                        center,radius,zero(real(T)),
                                        weightmethod,output,nothing,
                                        minlog2degree,globalmaxlog2degree,
                                        reltol,abstol,vectornorm)
        return output,evals
    end

    output=zeros(T,2,length(freqs)) #initialize accumulator to zero
    if method == :chirp
        subintegrals,workspaces = Matrix{T}(undef,2,length(freqs)),nothing
    else # :plain, :tone
        subintegrals,workspaces = fccquad_alloc(freqs,localmaxlog2degree,weightmethod,T)
    end
    interval_adaptive!(output,subintegrals,workspaces,center,radius,
                       prefactor,oscillator,freqs,branching,maxdepth,0,
                       minlog2degree,localmaxlog2degree,reltol,abstol,
                       method,weightmethod,T,vectornorm,maxchirp)
end
#adds results in place to output
function interval_adaptive!(output::AA,subintegrals::AA,workspaces,center::Real,radius::Real,
                            prefactor::Function,oscillator::Function,freqs::AA,
                            branching::Integer,maxdepth::Integer,depth::Integer,
                            minlog2N::Integer,maxlog2N::Integer,reltol::Real,abstol::Real,
                            method::Symbol,weightmethod::Symbol,T::Type,vectornorm::Function,
                            maxchirp::Integer)
    @assert method in interval_methods
    if method == :chirp #degree-adaptive FCC quadrature with linear chirp removal
        success,evals,base=chirpquad!(depth<maxdepth,subintegrals,center,radius,
                                      prefactor,oscillator,freqs,minlog2N,maxlog2N,
                                      weightmethod,T,reltol,abstol,vectornorm,maxchirp)
    elseif method == :tone #degree-adaptive FCC quadrature with tone removal
        success,evals,base=tonequad!(subintegrals,workspaces,center,radius,
                                     prefactor,oscillator,freqs,
                                     minlog2N,maxlog2N,weightmethod,
                                     T,reltol,abstol,vectornorm)
    elseif method == :plain #degree-adaptive FCC quadrature
        success,evals,base=adaptdegree!(prefactor,oscillator,freqs,center,radius,zero(real(T)),
                                        weightmethod,subintegrals,workspaces,
                                        minlog2N,maxlog2N,reltol,abstol,vectornorm)
    end
    if success || depth >= maxdepth
        output .+= subintegrals
    else
        shrink = inv(real(T)(branching))
        r = shrink*radius
        abstol = shrink*max(abstol, base*reltol)
        for t in 1-branching:2:branching-1
            evals += interval_adaptive!(output,subintegrals,workspaces,r*t+center,r,
                                        prefactor,oscillator,freqs,branching,maxdepth,depth+1,
                                        minlog2N,maxlog2N,reltol,abstol,
                                        method,weightmethod,T,vectornorm,maxchirp)[2]
        end
    end
    output,evals
end

"""
    ```
    fccquad_cc(amplitude::Function,
               angle::Function,
               freqs::AbstractArray{<:Real};
               kwargs...)
    ```

Integrates `amplitude(x)*cos(angle(x))*cos(w*x)*dx`
over a finite interval, for all `w` in `freqs`,
using a variant of Filon-Chenshaw-Curtis quadrature.

Use the `xmin` and `xmax` keyword arguments to specify
a finite domain of integration [xmin,xmax]. Default: [-1,1].

Required: `amplitude` and `angle` should be real-valued on [xmin,xmax].

Recommended: `amplitude` should not be highly oscillatory on [xmin,xmax].
The first derivative of `angle` should also not be highly oscillatory.

Outputs a tuple (results, function_evaluation_count) where results
is a 2xL real matrix, where L=length(freqs),
with row 1 containing the estimated integrals
and row 2 the discrepancies between row 1 and estimates made
using a reduced Chebyshev interpolation degree.

The `T` keyword argument specifies the working real scalar datatype.
(Default: `Float64`.)
The functions `amplitude` and `angle` should output
something of type `T` when the input has type `T`.

See `fccquad` for documentation of all other keyword arguments.

See also: `fccquad_cs`, `fccquad_sc`, `fccquad_ss`, `fccquad`.
"""
function fccquad_cc(amplitude::Function,angle::Function,freqs::AA{<:Real};
                    kwargs...)
    fccquad_trig(coscos!,amplitude,angle,freqs;kwargs...)
end

"Like fccquad_cc(), but integrates `amplitude(x)*cos(angle(x))*sin(w*x)*dx`"
function fccquad_cs(amplitude::Function,angle::Function,freqs::AA{<:Real};
                    kwargs...)
    fccquad_trig(cossin!,amplitude,angle,freqs;kwargs...)
end

"Like fccquad_cc(), but integrates `amplitude(x)*sin(angle(x))*cos(w*x)*dx`"
function fccquad_sc(amplitude::Function,angle::Function,freqs::AA{<:Real};
                    kwargs...)
    fccquad_trig(sincos!,amplitude,angle,freqs;kwargs...)
end

"Like fccquad_cc(), but integrates `amplitude(x)*sin(angle(x))*sin(w*x)*dx`"
function fccquad_ss(amplitude::Function,angle::Function,freqs::AA{<:Real};
                    kwargs...)
    fccquad_trig(sinsin!,amplitude,angle,freqs;kwargs...)
end

function fccquad_trig(transform!,amplitude::Function,angle::Function,freqs::AA;
                      T::Type=Float64,vectornorm=LinearAlgebra.norm,kwargs...)
    symmetricfreqs = Vector{eltype(freqs)}(undef,2length(freqs))
    for i in 1:length(freqs)
        w = freqs[i]
        symmetricfreqs[2i-1],symmetricfreqs[2i] = w,-w
    end
    workspace = Vector{T}(undef,length(freqs))
    pseudonorm(v) = trignorm!(v, workspace, transform!, vectornorm)

    oscillator(z) = exp(im*angle(z))
    cis,evals = fccquad(amplitude,oscillator,symmetricfreqs;T=Complex{T},vectornorm=pseudonorm,kwargs...)
    output = Matrix{T}(undef,2,length(freqs))
    for i in 1:2
        transform!(view(cis,i,:),view(output,i,:))
    end
    output,evals
end

function trignorm!(v, workspace, transform!, vectornorm)
    transform!(v, workspace)
    vectornorm(workspace)
end

function coscos!(exps, trigs)
    for i in 1:length(trigs)
        #2cos(a)cos(b)==cos(a+b)+cos(a-b)
        trigs[i] = 0.5(real(exps[2i-1])+real(exps[2i]))
    end
end

function cossin!(exps, trigs)
    for i in 1:length(trigs)
        #2cos(a)sin(b)==sin(a+b)-sin(a-b)
        trigs[i] = 0.5(imag(exps[2i-1])-imag(exps[2i]))
    end
end

function sincos!(exps, trigs)
    for i in 1:length(trigs)
        #2sin(a)cos(b)==sin(a+b)+sin(a-b)
        trigs[i] = 0.5(imag(exps[2i-1])+imag(exps[2i]))
    end
end

function sinsin!(exps, trigs)
    for i in 1:length(trigs)
        #2sin(a)sin(b)==cos(a-b)-cos(a+b)
        trigs[i] = 0.5(real(exps[2i])-real(exps[2i-1]))
    end
end

function __init__()
    cache_dir = @get_scratch!("precomputed_chebyshev_weights")
    cache_file = joinpath(cache_dir,"chirps.h5")
    isfile(cache_file) || Chirps.store(cache_file)
    Chirps.load(cache_file)
end

end #module


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    Base.precompile(Tuple{typeof(hcat),Matrix{Int64},Matrix{Float64}})   # time: 0.0145908
end


text:
export General
export get_info, print_info, get_info_deep
# General Matching
# ----------------

# * All matching distances find optimal COMPLETE matching by default
# * To instead find general optimal matching one can construct the generalised distance

struct General{T<:CompleteMatchingDistance}
    d::T
end

function get_cost_matrix_fixed(
    d_gen::General{S},
    X::Vector{T}, Y::Vector{T}
    ) where {S<:MatchDist, T}

    d_g = d_gen.d.ground_dist
    N, M = (length(X), length(Y))
    K = N + M
    C = zeros(K, K)
    for j in 1:K
        for i in 1:K
            if (i ≤ N) & (j ≤ M)
                C[i,j] = d_g(X[i],Y[j])
            elseif (i ≤ N) & (j > M)
                C[i,j] = d_g(X[i],nothing) # Penalty for non-inclusion of X entry
            elseif (i > N) & (j ≤ M)
                C[i,j] = d_g(Y[j],nothing) # Penalty for non-inclusion of Y entry
            end
        end
    end
    return C
end

function get_cost_matrix_fixed(
    d_gen::General{S},
    X::Vector{T}, Y::Vector{T}
    ) where {S<:FixPenMatchDist, T}

    d_g = d_gen.d.ground_dist
    penalty = d_gen.d.penalty
    N, M = (length(X), length(Y))
    K = N + M
    C = fill(penalty, K, K)
    for j in 1:K
        for i in 1:K
            if (i ≤ N) & (j ≤ M)
                C[i,j] = d_g(X[i],Y[j])
            elseif (i > N) & (j > M)
                C[i,j] = 0.0
            end
        end
    end
    return C
end

function get_cost_matrix_fixed(
    d_gen::General{S},
    X::Vector{T}, Y::Vector{T}
    ) where {S<:AvgSizeMatchDist, T}

    d_g = d_gen.d.ground_dist
    penalty = d_gen.d.penalty
    N, M = (length(X), length(Y))
    K = N + M
    C = zeros(K, K)
    x_mean_size, y_mean_size = (
        mean(x->d_g(x,nothing), X),
        mean(y->d_g(y,nothing), Y)
    )
    for j in 1:K
        for i in 1:K
            if (i ≤ N) & (j ≤ M)
                C[i,j] = d_g(X[i],Y[j])
            elseif (i ≤ N) & (j > M)
                C[i,j] = abs(d_g(X[i],nothing)-y_mean_size) + penalty # Penalty for non-inclusion of X entry
            elseif (i > N) & (j ≤ M)
                C[i,j] = abs(d_g(Y[j],nothing)-x_mean_size) + penalty # Penalty for non-inclusion of Y entry
            end
        end
    end
    return C

end

function get_cost_matrix_fixed(
    d_gen::General{S},
    X::Vector{T}, Y::Vector{T}
    ) where {S<:MinDistMatchDist, T}

    d_g = d_gen.d.ground_dist
    penalty = d_gen.d.penalty
    N, M = (length(X), length(Y))
    # C_v = view(C, 1:N, 1:M)
    C = pairwise_inbounds(d_g, X, Y)
    x_min_vec, y_min_vec = (
        map(minimum, eachrow(C)),
        map(minimum, eachcol(C))
    )
    C = [
        C [x + penalty for x ∈ x_min_vec, i in 1:N] ;
        [y + penalty for i in 1:M, y ∈ y_min_vec] zeros(M, N)
    ]
    return C

end