text:
using BenchmarkTools, Test, Knet
# @show uses print, disp uses display which prints out details of benchmark
macro disp(x); :(println($(string(x)));display($(esc(x)));println()); end
macro date(_x) :(println("$(now()) "*$(string(_x)));flush(STDOUT);@time $(esc(_x))) end
macro benchmarx(x); :(@benchmark ($x;sync())); end
macro gpu(x); if gpu()>=0; esc(x); end; end

using Knet: poolx, im2col!, col2im!, im2col_dims, conv4x, conv4w
@gpu sync()=Knet.cudaDeviceSynchronize()
o = []

ax = rand(100,100,100,10)
aw = rand(3,3,100,100)
@gpu kx = KnetArray(ax)
@gpu kw = KnetArray(aw)

@gpu @show @test isapprox(pool(ax;o...), pool(kx;o...))
@gpu @show @benchmarx pool($kx;$o...) # 0.766ms
@show @benchmark pool($ax;$o...) # 3.407ms, 22.362ms (with and without omp on ai-test)

ap = pool(ax;o...)
@gpu kp = pool(kx;o...)
@gpu @show @test isapprox(poolx(ax,ap,ap;o...), poolx(kx,kp,kp;o...))
@gpu @show @benchmarx poolx($kx,$kp,$kp;$o...) # 2.638ms
@show @benchmark poolx($ax,$ap,$ap;$o...) # 30.780ms, 67.446ms

@gpu @show @test isapprox(conv4(aw,ax;o...), conv4(kw,kx;o...))
@gpu @show @benchmarx conv4($kw,$kx;$o...) # 24.614ms
@show @benchmark conv4($aw,$ax;$o...) # 269.747ms, 346.614ms

ay = conv4(aw,ax;o...)
@gpu ky = conv4(kw,kx;o...)
@gpu @show @test isapprox(conv4w(aw,ax,ay;o...), conv4w(kw,kx,ky;o...))
@gpu @show @benchmarx conv4w($kw,$kx,$ky;$o...) # 28.509ms
@show @benchmark conv4w($aw,$ax,$ay;$o...) # 190.270ms, 287.182ms

@gpu @show @test isapprox(conv4x(aw,ax,ay;o...), conv4x(kw,kx,ky;o...))
@gpu @show @benchmarx conv4x($kw,$kx,$ky;$o...) # 45.799ms
@show @benchmark conv4x($aw,$ax,$ay;$o...) # 221.181ms, 245.694ms

x2 = similar(ax, im2col_dims(aw,ax,ay))
@show @benchmark im2col!($aw, $ax, $x2, 1, 0, 0, 1, 1, 0) # 0.326ms, 10.509ms
@show @benchmark col2im!($aw, $ax, $x2, 1, 0, 0, 1, 1, 0) # 3.367ms, 10.144ms

nothing

#=
# commit 359d3646 2018-09-22, julia 1.0.0 vs commit 4aa5f92f 2018-08-14, julia 0.6.4
# GPU:V100, CPU:Intel(R) Xeon(R) CPU E5-2680 v4 @ 2.40GHz

@test(isapprox(pool(ax; o...), pool(kx; o...))) = Test Passed
@benchmarx(pool($(Expr(:$, :kx)); $(Expr(:$, :o))...)) = Trial(143.338 μs) = Trial(154.221 μs)
@benchmark(pool($(Expr(:$, :ax)); $(Expr(:$, :o))...)) = Trial(3.107 ms) = Trial(3.123 ms)
@test(isapprox(poolx(ax, ap, ap; o...), poolx(kx, kp, kp; o...))) = Test Passed
@benchmarx(poolx($(Expr(:$, :kx)), $(Expr(:$, :kp)), $(Expr(:$, :kp)); $(Expr(:$, :o))...)) = Trial(434.826 μs) = Trial(433.798 μs)
@benchmark(poolx($(Expr(:$, :ax)), $(Expr(:$, :ap)), $(Expr(:$, :ap)); $(Expr(:$, :o))...)) = Trial(60.695 ms) = Trial(48.381 ms) -----
@test(isapprox(conv4(aw, ax; o...), conv4(kw, kx; o...))) = Test Passed
@benchmarx(conv4($(Expr(:$, :kw)), $(Expr(:$, :kx)); $(Expr(:$, :o))...)) = Trial(4.055 ms) = Trial(4.054 ms)
@benchmark(conv4($(Expr(:$, :aw)), $(Expr(:$, :ax)); $(Expr(:$, :o))...)) = Trial(223.892 ms) = Trial(302.637 ms) +++++
@test(isapprox(conv4w(aw, ax, ay; o...), conv4w(kw, kx, ky; o...))) = Test Passed
@benchmarx(conv4w($(Expr(:$, :kw)), $(Expr(:$, :kx)), $(Expr(:$, :ky)); $(Expr(:$, :o))...)) = Trial(4.699 ms) = Trial(4.603 ms)
@benchmark(conv4w($(Expr(:$, :aw)), $(Expr(:$, :ax)), $(Expr(:$, :ay)); $(Expr(:$, :o))...)) = Trial(181.501 ms) = Trial(195.529 ms)
@test(isapprox(conv4x(aw, ax, ay; o...), conv4x(kw, kx, ky; o...))) = Test Passed
@benchmarx(conv4x($(Expr(:$, :kw)), $(Expr(:$, :kx)), $(Expr(:$, :ky)); $(Expr(:$, :o))...)) = Trial(4.450 ms) = Trial(4.430 ms)
@benchmark(conv4x($(Expr(:$, :aw)), $(Expr(:$, :ax)), $(Expr(:$, :ay)); $(Expr(:$, :o))...)) = Trial(311.572 ms) = Trial(276.654 ms) ----
@benchmark(im2col!($(Expr(:$, :aw)), $(Expr(:$, :ax)), $(Expr(:$, :x2)), 1, 0, 0, 1, 1, 0)) = Trial(799.827 μs) = Trial(917.461 μs) ++++
@benchmark(col2im!($(Expr(:$, :aw)), $(Expr(:$, :ax)), $(Expr(:$, :x2)), 1, 0, 0, 1, 1, 0)) = Trial(3.571 ms) = Trial(3.554 ms)
=#


text:
export MIDIFile, readMIDIFile, writeMIDIFile
export BPM, ms_per_tick

"""
    MIDIFile <: Any
Type representing a file of MIDI data.

## Fields
* `format::UInt16` : The format of the file. Can be 0, 1 or 2.
* `tpq::Int16` : The time division of the track, ticks-per-quarter-note.
* `tracks::Array{MIDITrack, 1}` : The array of contained tracks.
"""
mutable struct MIDIFile
    format::UInt16 # Can be 0, 1 or 2
    tpq::Int16 # The time division of the track. Ticks per quarter note
    tracks::Vector{MIDITrack}
end
# Pretty print
function Base.show(io::IO, midi::MIDIFile) where {N}
	tnames = tracknames(midi)
	s = "MIDIFile (format=$(Int(midi.format)), tpq=$(midi.tpq)) "
	if any(!isequal(NOTRACKNAME), tnames) # we have tracknames
		s *= "with tracks:\n"
		for t in tnames
			s *= " "*t*"\n"
		end
	else # some track doesn't have a name
		s *= "with $(length(midi.tracks)) tracks"
	end
	print(io, s)
end

MIDIFile() = MIDIFile(0,960,MIDITrack[])

function readMIDIFileastype0(filename::AbstractString)
	MIDIFile = readMIDIFile(filename)
	if MIDIFile.format == 1
		type1totype0!(MIDIFile)
	end
	MIDIFile
end

"""
    readMIDIFile(filename::AbstractString)
Read a file into a `MIDIFile` data type.
"""
function readMIDIFile(filename::AbstractString)
    if length(filename) < 4 || filename[end-3:end] != ".mid"
		filename *= ".mid"
    end
    f = open(filename)

    midifile = MIDIFile()
    # Check that it's a valid MIDI file - first four bytes should spell MThd
    mthd = join(map(Char, read!(f, Array{UInt8}(undef, 4))))
    if mthd != MTHD
        error("Not a valid MIDI file. Expected first 4 bytes to spell 'MThd', got $(mthd)")
    end

    # Skip the next four bytes - this is the header size, and it's always equal to 6.
    skip(f, 4)

    # Read the format code. 0 = single track, 1 = multiple tracks, 2 = multiple songs
    # Remember - MIDI files store data in big-endian format, which is why ntoh is used
    midifile.format = ntoh(read(f, UInt16))

    # Get the number of tracks and time division
    numberoftracks = ntoh(read(f, UInt16))
    midifile.tpq = ntoh(read(f, Int16))
    midifile.tracks = [readtrack(f) for x in 1:numberoftracks]
    close(f)

    midifile
end

readMIDIFile() = readMIDIFile(testmidi())

"""
    writeMIDIFile(filename::AbstractString, data::MIDIFile)
Write a `MIDIFile` as a ".mid" file to the given filename.

    writeMIDIFile(filename::AbstractString, notes::Notes)
Create a `MIDIFile` directly from `notes`, using format 1.
"""
function writeMIDIFile(filename::AbstractString, data::MIDIFile)
    if length(filename) < 4 || filename[end-3:end] != ".mid" && filename[end-3:end] != ".MID"
      filename *= ".mid"
    end

    f = open(filename, "w")

    write(f, MTHD) # File identifier
    write(f, hton(convert(UInt32, 6))) # Header length
    write(f, hton(data.format))
    write(f, hton(convert(UInt16, length(data.tracks))))
    write(f, hton(data.tpq))

    map(track->writetrack(f, track), data.tracks)

    close(f)
    return data
end

function writeMIDIFile(filename::AbstractString, notes::Notes)
    if length(filename) < 4 || lowercase(filename[end-3:end]) != ".mid"
      filename *= ".mid"
    end

    track = MIDITrack()
    addnotes!(track, notes)
    midi = MIDIFile(1, notes.tpq, [track])
    writeMIDIFile(filename, midi)
    return midi
end

"""
    BPM(midi)
Return the BPM where the given `MIDIFile` was exported at.
"""
function BPM(t::MIDI.MIDIFile)
    # META-event list:
    tttttt = Vector{UInt32}()
    # Find the one that corresponds to Set-Time:
    # The event tttttt corresponds to the command
    # FF 51 03 tttttt Set Tempo (in microseconds per MIDI quarter-note)
    # See here (page 8):
    # http://www.cs.cmu.edu/~music/cmsip/readings/Standard-MIDI-file-format-updated.pdf
    for event in t.tracks[1].events
        if typeof(event) == MetaEvent
            if event.metatype == 0x51
                tttttt = deepcopy(event.data)
                break
            end
        end
    end
    # Ensure that tttttt is with correct form (first entry should be 0x00)
    if tttttt[1] != 0x00
        pushfirst!(tttttt, 0x00)
    else
        # Handle correctly "incorrect" cases where 0x00 has entered more than once
        tttttt = tttttt[findin(tttttt, 0x00)[end]:end]
    end

    # Get the microsecond number from tttttt
    u = ntoh(reinterpret(UInt32, tttttt)[1])
    μs = Int64(u)
    # BPM:
    bpm = 60000000/μs
end

"""
    ms_per_tick(tpq, bpm)
    ms_per_tick(midi::MIDIFile)
Return how many miliseconds is one tick, based
on the beats per minute `bpm` and ticks per quarter note `tpq`.
"""
ms_per_tick(midi::MIDI.MIDIFile, bpm = BPM(midi)) = ms_per_tick(midi.tpq, bpm)
ms_per_tick(tpq, bpm) = (1000*60)/(bpm*tpq)

getnotes(midi::MIDIFile, trackno = midi.format == 0 ? 1 : 2) =
getnotes(midi.tracks[trackno], midi.tpq)

ext: jl


text:
using LinearAlgebra, BesselK, BenchmarkTools, Printf, StaticArrays
include("shared.jl")

const GRIDN  = 24
const GRID1D = range(0.0, 1.0, length=GRIDN)
const PTS    = map(x->SVector{2,Float64}(x[1], x[2]),
                   vec(collect.(Iterators.product(GRID1D, GRID1D))))

function checkcovmat(fn, p, v)
  _p     = pv(1.0, p, v)
  M      = assemble_matrix(fn, PTS, _p)
  em     = eigmin(M)
  Mf     = cholesky!(M, check=false)
  if issuccess(Mf)
    ld = logdet(Mf)
  else
    ld = NaN
  end
  (issuccess(Mf) ? "S" : "F", em, ld)
end

const VRANGE = (0.4, 1.25, 3.5)
const PRANGE = (0.01, 1.0, 100.0)

# For now, I might even keep the timing in seconds.
for (j, (v, p)) in enumerate(Iterators.product(VRANGE, PRANGE))
  _p  = pv(1.0, p, v)
  # test 1: assembly time.
  t_u = @belapsed assemble_matrix(matern_us,   $PTS, $_p) samples=8
  t_a = @belapsed assemble_matrix(matern_amos, $PTS, $_p) samples=8
  # test 2: Cholesky success or failure:
  (s_u, em_u, ld_u) = checkcovmat(matern_us,   p, v)
  (s_a, em_a, ld_a) = checkcovmat(matern_amos, p, v)
  # PRINTING:
  # (v,p) pair
  # times (u,p)
  # eigmins (u,p)
  # eigmin difference
  # logdets (u.p)
  # logdet difference
  if j != length(VRANGE)*length(PRANGE)
    @printf "(%1.3f, %1.2f) & %1.1e & %1.1e & %1.2e & %1.2e & %1.2e & %1.2e & %1.2e & %1.2e \\\\\n" p v t_u t_a em_u em_a abs(em_u-em_a) ld_u ld_a abs(ld_u-ld_a)
  else
    @printf "(%1.3f, %1.2f) & %1.1e & %1.1e & %1.2e & %1.2e & %1.2e & %1.2e & %1.2e & %1.2e\n" p v t_u t_a em_u em_a abs(em_u-em_a) ld_u ld_a abs(ld_u-ld_a)
  end
end


text:
export sample_vtree

using Random
using LogicCircuits: Vtree, PlainVtreeLeafNode

#####################
# Sampling functions
#####################

"""Split `X` into two partitions `A` and `B`, where `A` is a Bernoulli sample of each element in
`X` with probability `p` and `B=X∖A`. Guarantees at least one element in `A` and `B`."""
function bernoulli_partition(X::Vector{Int}, p::Float64)::Tuple{Vector{Int}, Vector{Int}}
    n = length(X)
    a = rand(1:n)
    b = ((rand(0:n-2)+a)%n)+1
    A, B = Int[X[a]], Int[X[b]]
    for (i, x) ∈ enumerate(X)
        if i == a continue end
        if i == b continue end
        if rand() > p push!(B, x)
        else push!(A, x) end
    end
    return A, B
end

"Samples a Vtree with a right bias of `p`. If `p<0`, then uniformly sample vtrees."
function sample_vtree(n::Int, p::Float64)::Vtree
    passdown(x::Int)::Vtree = PlainVtreeLeafNode(x)
    function passdown(X::Vector{Int})::Vtree
        R, L = bernoulli_partition(X, p)
        return Vtree(passdown(length(L) == 1 ? L[1] : L), passdown(length(R) == 1 ? R[1] : R))
    end
    return p < 0 ? Vtree(n, :random) : passdown(shuffle!(collect(1:n)))
end


text:
"""
    add_intercept(A)

Insert an intercept column (column of ones) at the beginning of a 2d array.

# Arguments

- A = 2d array of floats

# Value

Returns A with an intercept column

"""
function add_intercept(A::AbstractArray{Float64,2})

    # Check that the first column isn't a column of ones
    if all(A[:,1].==1)
        print("The first column of this array is already a column of ones. Are you sure you wanted to add an intercept?")
    end

    # Concatenate a column of ones to the beginning of A
    return hcat(ones(size(A,1)), A)
end

"""
    remove_intercept(A)

Remove the intercept column, assumed to be the first column of a 2d array.

# Arguments

- A = 2d array of floats

# Value

Returns A without the intercept column

"""
function remove_intercept(A::AbstractArray{Float64,2})

    # Check that the first column is a column of ones
    if any(A[:,1].!=1)
        print("The first column of this array does not look like an itercept
               since it is not a column of ones. Are you sure you wanted to
               remove it?")
    end

    # Remove the first column of A
    return A[:,2:end]
end

"""
    shuffle_rows(A)

Shuffles rows of a 2d array

# Arguments

- A = 2d array of floats

# Value

Returns A with rows shuffled

"""
function shuffle_rows(A::AbstractArray{Float64,2})

    return A[Random.shuffle(1:size(A,1)),:]
end

"""
    shuffle_cols(A)

Shuffles columns of a 2d array

# Arguments

- A = 2d array of floats

# Value

Returns A with columns shuffled

"""
function shuffle_cols(A::AbstractArray{Float64,2})

    return A[:,Random.shuffle(1:size(A,2))]
end


text:
# using EffectiveWaves, LinearAlgebra
# import StaticArrays: SVector

import MultipleScattering: outgoing_translation_matrix

function discrete_system(ω::T, source::AbstractSource{T,Acoustic{T,Dim}}, material::Material{Dim,Sphere{T,Dim}}, ::WithoutSymmetry{Dim}; kws...) where {T,Dim}

    return discrete_system(ω, source, material, AzimuthalSymmetry{Dim}(); kws...)
end

"""
    discrete_system(k,r1_vec::Vector,r2,θ2; basis_order::Int = 3)

documentation
"""
function discrete_system(ω::T, source::AbstractSource{T,Acoustic{T,Dim}}, material::Material{Dim,Sphere{T,Dim}}, ::AbstractAzimuthalSymmetry{Dim};
        basis_order::Int = 1,
        basis_field_order::Int = 2,
        legendre_order::Int = basis_field_order + 1,
        mesh_points::Int = Int(round(sqrt(1.1 * (basis_field_order) * legendre_order ))) + 1,
        rtol::T = 1e-2,
        maxevals::Int = Int(2e4),
        pair_corr = hole_correction_pair_correlation
    ) where {T,Dim}

    if length(material.species) > 1
        @warn "discrete_system has only been implemented for 1 species for now. Will use only first specie."
    end

    s1 = material.species[1]
    scale_number_density = one(T) - one(T) / material.numberofparticles
    bar_numdensity = scale_number_density * number_density(s1)

    R = outer_radius(material.shape)

    gs = regular_spherical_coefficients(source)(basis_field_order,origin(material.shape),ω);

    v = regular_basis_function(source.medium,  ω)

    Uout = outgoing_translation_matrix(ω, source.medium, material;
        basis_order = basis_order, tol = rtol
    )

    t_matrices = get_t_matrices(source.medium, material.species, ω, basis_order)
    t_diags = diag.(t_matrices)

    rθφ2xyz = radial_to_cartesian_coordinates

    r1s = LinRange(0,R-outer_radius(s1), mesh_points)
    θ1s = LinRange(0,π, mesh_points)

    ls, ms = spherical_harmonics_indices(basis_order)
    azi_inds(m::Int) = lm_to_spherical_harmonic_index.(abs(m):basis_field_order,-m)

    len = basisorder_to_basislength(Acoustic{T,Dim}, basis_order)
    len_p = sum(legendre_order for nd in 1:len for i in azi_inds(ms[nd]))

    function incident_coefficients(r1s::AbstractVector{T},θ1s::AbstractVector{T})
        lm2n = lm_to_spherical_harmonic_index

        coefs = [
            begin
                vs = v(basis_order + basis_field_order, rθφ2xyz(SVector(r1,θ1,zero(T))))
                data = [
                    t_diags[1][lm2n(l,m)] *
                    sum(
                        gaunt_coefficient(dl,0,l,m,l1,-m) * vs[lm2n(l1,-m)] * gs[lm2n(dl,0)]
                    for dl in 0:basis_field_order for l1 in max(abs(m),abs(dl-l)):(dl+l))
                for l = 0:basis_order for m = -l:l]

                data
            end
        for r1 in r1s, θ1 in θ1s][:];

        return vcat(coefs...)
    end

    function field_basis(rθφ::AbstractVector{T})
        Ys = spherical_harmonics(basis_field_order, rθφ[2], rθφ[3]);
        P = Legendre{T}()

        # [P_0(cos(θ)), …, P_(legendre_order-1)(cos(θ))]
        Prs = P[2 * rθφ[1] / (R - outer_radius(s1)) - one(T), 1:legendre_order]

        return [Pr * Yθ for Pr in Prs, Yθ in Ys]
    end
    # function field_basis(rθ::AbstractVector{T})
    #     P = Legendre{T}()
    #
    #     # [P_0(cos(θ)), …, P_(legendre_order-1)(cos(θ))]
    #     Pθs = P[cos(rθ[2]), 1:legendre_order]
    #     Prs = P[2 * rθ[1] / (R - outer_radius(s1)) - one(T), 1:legendre_order]
    #
    #     # [Pr * Pθ for Pr in Prs, Pθ in Pθs][:]
    #     return (Prs * transpose(Pθs))[:]
    # end

    function kernel_function(rθ1::SVector{2,T})
        x1 = rθφ2xyz(SVector(rθ1[1],rθ1[2],zero(T)))
        Kzero = zeros(Complex{T},len,len_p)

        fun = function (rθφ::SVector{3,T})
            x2 = rθφ2xyz(rθφ)
            if pair_corr(x1,s1,x2,s1) ≈ zero(T)
                return Kzero
            end
            basis2 = field_basis(rθφ)

            U = Uout(x1 - x2)
            U = U .* (bar_numdensity * pair_corr(x1,s1,x2,s1) * sin(rθφ[2]) * rθφ[1]^2)

            return reshape(
                [
                    t_diags[1][n] * U[nd,n] * b2
                for nd in 1:len for b2 in basis2[:,azi_inds(ms[nd])][:] for n in 1:len],
            (len, :))
        end

        return fun
    end

    function δφj(rθφ1::AbstractVector{T})
        basis1 = field_basis(rθφ1)
        return reshape(
            [
                (nd == n) ? b1 : zero(Complex{T})
            for nd in 1:len for b1 in basis1[:,azi_inds(ms[nd])][:] for n in 1:len],
        (len, len_p))
    end

    δφj(rθ1::SVector{2,T}) = δφj(SVector(rθ1[1],rθ1[2],zero(T)))

    test_ker = kernel_function(SVector(mean(r1s),θ1s[1]))
    (I,E) = hcubature(test_ker, SVector(0.0,0.0,-π), SVector(R-outer_radius(s1),π,π);
        rtol=rtol, maxevals=maxevals
    );

    println("The estimated max coefficient of the integrated kernel is:")
    println("I = ", maximum(abs.(I)) )
    println("with an estimated error of: ")
    println("E = ", E)

    Ks = [
        begin
            rθ1 = SVector(r1,θ1)
            ker = kernel_function(rθ1)
            ker_integrated = hcubature(ker, SVector(0.0,0.0,-π), SVector(R-outer_radius(s1),π,π);
                rtol=rtol, maxevals=maxevals
            )[1]

            δφj(rθ1) - ker_integrated
        end
    for r1 in r1s, θ1 in θ1s][:];

    bigK = vcat(Ks...);

    bs = incident_coefficients(r1s,θ1s);

    Fs = bigK \ bs;

    ## Alternative:
    # as = inv(transpose(conj.(bigK)) * bigK) * transpose(conj.(bigK)) * bs;

    if norm(bigK * Fs - bs) / norm(bs) > rtol
        @warn "Numerical solution has a relative residual error of $(norm(bigK * Fs - bs) / norm(bs)), where the requested relative tolernance was: $rtol. This residual error can be decreased by increasing the basis_field_order (current value: $basis_field_order) for the field."
    end

    function scattered_field(xs::Vector{T})
        rθφ = cartesian_to_radial_coordinates(xs)
        return δφj(rθφ) * Fs
        # The factor exp(-im * m * φ) is due to azimuthal symmetry
        # azi_factor = exp.((-im*rθφ[3]) .* ms)
        # return azi_factor .* (as * field_basis(rθφ[1:2]))
    end

    return scattered_field
end

function discrete_system(ω::T, source::AbstractSource{T,Acoustic{T,Dim}}, material::Material{Dim,Sphere{T,Dim}}, ::RadialSymmetry{Dim};
        basis_order::Int = 1,
        basis_field_order::Int = 2,
        legendre_order::Int = basis_field_order + 1,
        mesh_points::Int = Int(round(sqrt(1.1 * (basis_field_order) * legendre_order ))) + 1,
        rtol::T = 1e-2,
        maxevals::Int = Int(2e4),
        pair_corr = hole_correction_pair_correlation
    ) where {T,Dim}

    if length(material.species) > 1
        @warn "discrete_system has only been implemented for 1 species for now. Will use only first specie."
    end

    s1 = material.species[1]
    scale_number_density = one(T) - one(T) / material.numberofparticles
    bar_numdensity = scale_number_density * number_density(s1)

    R = outer_radius(material.shape)

    gs = regular_spherical_coefficients(source)(basis_field_order,origin(material.shape),ω);

    v = regular_basis_function(source.medium,  ω)

    Uout = outgoing_translation_matrix(ω, source.medium, material;
        basis_order = basis_order, tol = rtol
    )

    ns_to_l1s = lm_to_spherical_harmonic_index.(0:basis_field_order,0)

    t_matrices = get_t_matrices(source.medium, material.species, ω, basis_order)
    t_diags = diag.(t_matrices)

    rθφ2xyz = radial_to_cartesian_coordinates

    r1s = LinRange(0,R-outer_radius(s1), mesh_points)
    θ1s = LinRange(0,π, mesh_points)

    len = basisorder_to_basislength(Acoustic{T,Dim}, basis_order)
    len_p = legendre_order

    function incident_coefficients(r1s::AbstractVector{T})
        lm2n = lm_to_spherical_harmonic_index

        coefs = [
            begin
                vs = v(basis_order + basis_field_order, rθφ2xyz(SVector(r1,zero(T),zero(T))))
                data = [
                    t_diags[1][lm2n(l,0)] *
                        sqrt(4pi) * (-1)^l * vs[lm2n(l,0)] * gs[lm2n(0,0)]
                for l = 0:basis_order]
            end
        for r1 in r1s][:];

        return vcat(coefs...)
    end

    function field_basis(rθφ::AbstractVector{T})
        Ys = spherical_harmonics(basis_field_order, rθφ[2], rθφ[3]);
        P = Legendre{T}()

        Prs = P[2 * rθφ[1] / (R - outer_radius(s1)) - one(T), 1:legendre_order]

        return [Pr * Yθ for Pr in Prs, Yθ in Ys]
    end

    function kernel_function(rθ1::SVector{2,T})
        x1 = rθφ2xyz(SVector(rθ1[1],rθ1[2],zero(T)))
        Kzero = zeros(Complex{T},len,len_p)

        fun = function (rθφ::SVector{3,T})
            x2 = rθφ2xyz(rθφ)
            if pair_corr(x1,s1,x2,s1) ≈ zero(T)
                return Kzero
            end
            basis2 = field_basis(rθφ)

            U = Uout(x1 - x2)
            U = U .* (bar_numdensity * pair_corr(x1,s1,x2,s1) * sin(rθφ[2]) * rθφ[1]^2)

            return reshape(
                [
                    t_diags[1][n] * U[nd,n] * b2
                for nd in 1:len for b2 in basis2[:,azi_inds(ms[nd])][:] for n in 1:len],
            (len, :))
        end

        return fun
    end

    function δφj(rθφ1::AbstractVector{T})
        basis1 = field_basis(rθφ1)
        return reshape(
            [
                (nd == n) ? b1 : zero(Complex{T})
            for nd in 1:len for b1 in basis1[:,azi_inds(ms[nd])][:] for n in 1:len],
        (len, len_p))
    end

    δφj(rθ1::SVector{2,T}) = δφj(SVector(rθ1[1],rθ1[2],zero(T)))

    test_ker = kernel_function(SVector(mean(r1s),θ1s[1]))
    (I,E) = hcubature(test_ker, SVector(0.0,0.0,-π), SVector(R-outer_radius(s1),π,π);
        rtol=rtol, maxevals=maxevals
    );

    println("The estimated max coefficient of the integrated kernel is:")
    println("I = ", maximum(abs.(I)) )
    println("with an estimated error of: ")
    println("E = ", E)

    Ks = [
        begin
            rθ1 = SVector(r1,θ1)
            ker = kernel_function(rθ1)
            ker_integrated = hcubature(ker, SVector(0.0,0.0,-π), SVector(R-outer_radius(s1),π,π);
                rtol=rtol, maxevals=maxevals
            )[1]

            δφj(rθ1) - ker_integrated
        end
    for r1 in r1s, θ1 in θ1s][:];

    bigK = vcat(Ks...);

    bs = incident_coefficients(r1s,θ1s);

    Fs = bigK \ bs;

    ## Alternative:
    # as = inv(transpose(conj.(bigK)) * bigK) * transpose(conj.(bigK)) * bs;

    if norm(bigK * Fs - bs) / norm(bs) > rtol
        @warn "Numerical solution has a relative residual error of $(norm(bigK * Fs - bs) / norm(bs)), where the requested relative tolernance was: $rtol. This residual error can be decreased by increasing the basis_field_order (current value: $basis_field_order) for the field."
    end

    function scattered_field(xs::Vector{T})
        rθφ = cartesian_to_radial_coordinates(xs)
        return δφj(rθφ) * Fs
        # The factor exp(-im * m * φ) is due to azimuthal symmetry
        # azi_factor = exp.((-im*rθφ[3]) .* ms)
        # return azi_factor .* (as * field_basis(rθφ[1:2]))
    end

    return scattered_field
end

"""
    outgoing_translation_matrix(ω, ::Acoustic, material::Material{Dim,Sphere{T,Dim}};
        basis_order = 2, tol = 1e-3)

    return a function U where U(X) for X ∈ material.shape gives the the outgoing translation matrix
"""
function outgoing_translation_matrix(ω::T, medium::Acoustic{T,Dim}, material::Material{Dim,Sphere{T,Dim}};
        basis_order::Int = 2, tol::T = 1e-3) where {T,Dim}

    L = basisorder_to_basislength(typeof(medium), basis_order);
    R = outer_radius(material.shape) - minimum(outer_radius.(material.species))
    a12 = T(2) * minimum(outer_radius(s) * s.exclusion_distance for s in material.species)

    rθφ2xyz = radial_to_cartesian_coordinates

    # Just to get the type Tinter.
        xs = LinRange(a12,2R,4);
        data = [zeros(Complex{T},L,L) for x in xs];
        inter = interpolate(
            [data[j][1] for j in CartesianIndices(data)],
            BSpline(Cubic(Line(OnGrid())))
        );
        Tinterx = typeof(scale(inter, xs))

        data = [ zeros(Complex{T},L,L) for x in xs, y in xs, z in xs];
        inter = interpolate(
            [data[j][1] for j in CartesianIndices(data)],
            BSpline(Cubic(Line(OnGrid())))
        );
        Tinterxyz = typeof(scale(inter, xs, xs, xs))

    # Determine the number of mesh points
        Ns = 10:4:66;
        Nθs = 10:4:62;

        xs = LinRange(a12,2R,Ns[end]);
        θs = LinRange(zero(T),T(pi),Nθs[end]);

        datax = [
            outgoing_translation_matrix(medium, basis_order, ω, [x,zero(T),zero(T)])
        for x in xs];

        dataθ = [
            outgoing_translation_matrix(medium, basis_order, ω, rθφ2xyz(SVector(a12,θ,zero(T))))
        for θ in θs];

        CI1s = CartesianIndices(datax[1]);
        inter_xs = Array{Tinterx}(undef, size(CI1s));
        inter_θs = Array{Tinterx}(undef, size(CI1s));

        for i in CI1s
            inter = interpolate(
                [datax[j][i] for j in CartesianIndices(datax)],
                BSpline(Cubic(Line(OnGrid())))
            );
            inter_xs[i] = scale(inter, xs)

            inter = interpolate(
                [dataθ[j][i] for j in CartesianIndices(dataθ)],
                BSpline(Cubic(Line(OnGrid())))
            );
            inter_θs[i] = scale(inter, θs)
        end

        errorsx = Array{T}(undef, length(Ns))
        for n in eachindex(Ns)
            xNs = LinRange(a12,2R,Ns[n]);
            inter_xNs = Array{Tinterx}(undef, size(CI1s));

            for i in CI1s
                inter = interpolate(inter_xs[i].(xNs),BSpline(Cubic(Line(OnGrid()))));
                inter_xNs[i] = scale(inter, xNs)
            end

            errorsx[n] = norm(norm(inter_xNs[i].(xs) - inter_xs[i].(xs)) for i in CI1s) / norm(norm(inter_xs[i].(xs)) for i in CI1s)
        end
        n = findfirst(errorsx .< tol)
        Nr = Ns[n]

        errorsθ = Array{T}(undef, length(Nθs))
        for n in eachindex(Nθs)
            θNs = LinRange(zero(T),T(pi),Nθs[n]);
            inter_θNs = Array{Tinterx}(undef, size(CI1s));

            for i in CI1s
                inter = interpolate(inter_θs[i].(θNs),BSpline(Cubic(Line(OnGrid()))));
                inter_θNs[i] = scale(inter, θNs)
            end

            errorsθ[n] = norm(norm(inter_θNs[i].(θs) - inter_θs[i].(θs)) for i in CI1s) / norm(norm(inter_θs[i].(θs)) for i in CI1s)
        end
        n = findfirst(errorsθ .< tol)
        Nθ = Nθs[n]

    # interpolate the whole spherical region based on the xs mesh size
    rs = LinRange(a12,2R,Nr);
    φs = LinRange(-T(pi),T(pi),Nθ);
    θs = LinRange(zero(T),T(pi),Nθ);

    data = [
            outgoing_translation_matrix(medium, basis_order, ω, rθφ2xyz(SVector(r,θ,φ)))
    for r in rs, θ in θs, φ in φs];

    # reorganise the data to interpolate each element of outgoing_translation_matrix in terms of x, y, z.
    CI1s = CartesianIndices(data[1]);
    inters = Array{Tinterxyz}(undef, size(CI1s));

    for i in CI1s
        inter = interpolate(
            [data[j][i] for j in CartesianIndices(data)],
            BSpline(Cubic(Line(OnGrid())))
        );
        inters[i] = scale(inter, rs, θs, φs)
    end

    function Uout(X::AbstractVector{T})
        if norm(X) < a12
            zeros(Complex{T},L,L)
        else
            rθφ = cartesian_to_radial_coordinates(X)
            map(CI1s) do i
                inters[i](rθφ[1], rθφ[2], rθφ[3])
            end
        end
    end

    return Uout
end


text:
using DiffEqFlux, OrdinaryDiffEq, Plots

# model + problem parameters
(aₐ, aₙ, bₙ, cₙ, dₙ, aₘ, bₘ, cₘ, dₘ) = Float32.([-0.3, 19.373, -31.023, -9.717, -1.948, 40.44, -64.015, 2.922, -11.803])
(m, I_yy, S, d, ωₐ, ζₐ) = Float32.([204.02, 247.439, 0.0409, 0.2286, 150.0, 0.7])
(g, ρ₀, H, γₐ, Rₐ, T₀, λ) = Float32.([9.8, 1.225, 8435.0, 1.4, 286.0, 288.15, 0.0065])
(a_max, α_max, δ_max, q_max, M_max, h_max) = Float32.([100.0, deg2rad(30), deg2rad(25), deg2rad(5), 4, 11E3])

(k_a, k_δ, k_R) = Float32.([10000.0, 0.01, 0.001])
t_span = (0.0f0, 3.0f0)
x₀ = Float32.([5E3; 900; zeros(7)])
a_z_cmd = 50.0f0

# cost function
L(a_z, a_z_cmd, δ_c) = k_a * ((a_z - a_z_cmd) / a_max)^2 + k_δ * (δ_c / δ_max)^2
# + k_q * (q / q_max)^2
ϕ(a_z, a_z_cmd) = ((a_z - a_z_cmd) / a_max)^2
R(p_NN) = k_R * p_NN' * p_NN

# NN construction
dim_hidden = 64
(K_lb, K_ub) = Float32.([-0.2, 0])
Π = Chain(
    Dense(3, dim_hidden, tanh),
    Dense(dim_hidden, dim_hidden, tanh),
    Dense(dim_hidden, 3)
    # u -> (K_ub - K_lb) * σ.(u) .+ K_lb
)
(p_NN, restruct) = Flux.destructure(Π)

# plant + controller dynamics
function fwd_dynamics_model!(dx, x, p_NN, t)
    (h, V, α, q, θ, δ, δ̇, x_c, _) = x

    ρ = ρ₀ * exp(-h / H)
    Vₛ = sqrt(γₐ * Rₐ * (T₀ - λ * h))
    M = V / Vₛ
    γ = θ - α
    Q = 0.5f0 * ρ * V^2

    C_A = aₐ
    C_N = aₙ * α^3 + bₙ * α * abs(α) + cₙ * (2.0f0 - M / 3.0f0) * α + dₙ * δ
    C_M = aₘ * α^3 + bₘ * α * abs(α) + cₘ * (-7.0f0 + 8.0f0 * M / 3.0f0) * α + dₘ * δ

    K_A, K_I, K_R = min.( restruct(p_NN)([α / α_max; M / M_max; h / h_max]) , K_ub )
    δ_c = -K_I * x_c - K_R * q

    α̇ = Q * S / m / V * (C_N * cos(α) - C_A * sin(α)) + g / V * cos(γ) + q
    a_z = V * (q - α̇)

    dx[1] = V * sin(γ)
    dx[2] = Q * S / m * (C_N * sin(α) + C_A * cos(α)) - g * sin(γ)
    dx[3] = α̇
    dx[4] = Q * S * d / I_yy * C_M
    dx[5] = q
    dx[6] = δ̇
    dx[7] = -ωₐ^2 * (δ - δ_c) - 2.0f0 * ζₐ * ωₐ * δ̇
    dx[8] = -K_A * (a_z - a_z_cmd) + q - a_z_cmd / V
    dx[9] = L(a_z, a_z_cmd, δ_c)
end

function eval_output(x)
    (h, V, α, q, θ, δ, _, _, _) = x

    ρ = ρ₀ * exp(-h / H)
    Vₛ = sqrt(γₐ * Rₐ * (T₀ - λ * h))
    M = V / Vₛ
    γ = θ - α
    Q = 0.5f0 * ρ * V^2

    C_A = aₐ
    C_N = aₙ * α^3 + bₙ * α * abs(α) + cₙ * (2.0f0 - M / 3.0f0) * α + dₙ * δ

    α̇ = Q * S / m / V * (C_N * cos(α) - C_A * sin(α)) + g / V * cos(γ) + q
    return V * (q - α̇)
end

function eval_tune_params(p_NN)
    return function (x)
        (h, V, α) = x[1:3]

        Vₛ = sqrt(γₐ * Rₐ * (T₀ - λ * h))
        M = V / Vₛ

        return min.( restruct(p_NN)([α / α_max; M / M_max; h / h_max]) , K_ub )
    end
end

# ODE problem
prob = ODEProblem(fwd_dynamics_model!, x₀, t_span)

# loss construction
function loss()
    fwd_sol = Array(solve(prob, Tsit5(), p = p_NN, reltol = 1e-4, abstol = 1e-8, saveat = 0.1; sensealg = InterpolatingAdjoint(autojacvec = ZygoteVJP())))
    # fwd_sol = Array(solve(prob, Tsit5(), p = p_NN, saveat = 0.1; sensealg = InterpolatingAdjoint(autojacvec = ZygoteVJP())))

    return ϕ(eval_output(fwd_sol[:, end]), a_z_cmd) + R(p_NN) + fwd_sol[end, end] #, fwd_sol
end

# learning progress callback setup
cb_progress = function ()
# cb_progress = function (p, loss_val, pred)
    loss_val = loss()
    println("Loss = $loss_val\n")
    # display(scatter([vec(mapslices(eval_output, pred, dims = 1)) rad2deg.(pred[3:4, :]')], ylabel = ["\$a_{z} [m/s^{2}]\$" "\$\\alpha [deg]\$" "\$q [deg/s]\$"], label = :false, layout = (3, 1)))
    return false
end

function view_result(p_NN)
    fwd_sol = solve(prob, Tsit5(), p = p_NN, reltol = 1e-8, abstol = 1e-8, saveat = 0.01)
    a_z = eval_output.(fwd_sol.u)
    K = eval_tune_params(p_NN).(fwd_sol.u)

    f_aq = plot(fwd_sol.t, [a_z rad2deg.(fwd_sol[4,:])], label = :false, xlabel = "\$t [s]\$", ylabel = ["\$a_{z} [m/s^{2}]\$" "\$q [deg/s]\$"], layout = (2, 1))
    f_K = plot(fwd_sol.t, hcat(K...)', xlabel = "\$t [s]\$", ylabel = ["\$K_{A}\$" "\$K_{I}\$" "\$K_{R}\$"], label = :false, layout = (3, 1))
    display(plot(f_aq, f_K, layout = (1, 2)))

    f_hVα = plot(fwd_sol.t, [fwd_sol[1:2,:]' rad2deg.(fwd_sol[3,:])], label = :false, xlabel = "\$t [s]\$", ylabel = ["\$h [m]\$" "\$V [m/s]\$" "\$\\alpha [deg]\$"], layout = (3, 1))
    display(f_hVα)
    return fwd_sol
end

# view pre-training result
view_result(p_NN)

## NN training
println("1st Phase\n")
data = Iterators.repeated((), 2000)
Flux.train!(loss, Flux.params(p_NN), data, ADAM(0.5); cb = cb_progress)

##
println("2nd Phase\n")
Flux.train!(loss, Flux.params(p_NN), data, ADAM(); cb = cb_progress)

##
println("3rd Phase\n")
Flux.train!(loss, Flux.params(p_NN), data, Nesterov(); cb = cb_progress)

# view post-training result
fwd_sol = view_result(p_NN)


text:
#define	E_OK			(0)
#define	E_SYSTEM_ERROR	 	(-1)
#define	E_BAD_ARGUMENT	 	(-2)
#define	E_POSTED	 	(-3)
#define	E_CONNECTED	 	(-4)
#define	E_BAD_STATE	 	(-5)
#define	E_NO_ROOM	 	(-6)
#define	E_NOT_POSTED		(-7)
#define	E_UNKNOWN_COMMAND	(-8)
#define	E_NO_MATCH		(-9)
#define	E_NOT_SELECTABLE	(-10)
#define	E_NOT_CONNECTED	        (-11)
#define	E_REQUEST_DENIED	(-12)
#define	E_INVALID_FIELD	        (-13)
#define	E_CURRENT		(-14)
# this eti header is for ncurses->forms library
#ifndef FORM_PRIV_H
# typedef void *FIELD_CELL;
# #endif
#
# #ifndef NCURSES_FIELD_INTERNALS
# #define NCURSES_FIELD_INTERNALS /* nothing */
# #endif
#
# typedef int Form_Options;
# typedef int Field_Options;
#
#   /**********
#   *  _PAGE  *
#   **********/
#
# typedef struct {
#   short pmin;		/* index of first field on page			*/
#   short pmax;		/* index of last field on page			*/
#   short smin;		/* index of top leftmost field on page		*/
#   short smax;		/* index of bottom rightmost field on page	*/
# } _PAGE;
#
#   /**********
#   *  FIELD  *
#   **********/
#
# typedef struct fieldnode {
#   unsigned short	status;		/* flags			*/
#   short			rows;		/* size in rows			*/
#   short			cols;		/* size in cols			*/
#   short			frow;		/* first row			*/
#   short			fcol;		/* first col			*/
#   int			drows;		/* dynamic rows			*/
#   int			dcols;		/* dynamic cols			*/
#   int			maxgrow;	/* maximum field growth		*/
#   int			nrow;		/* off-screen rows		*/
#   short			nbuf;		/* additional buffers		*/
#   short			just;		/* justification		*/
#   short			page;		/* page on form			*/
#   short			index;		/* into form -> field		*/
#   int			pad;		/* pad character		*/
#   chtype		fore;		/* foreground attribute		*/
#   chtype		back;		/* background attribute		*/
#   Field_Options		opts;		/* options			*/
#   struct fieldnode *	snext;		/* sorted order pointer		*/
#   struct fieldnode *	sprev;		/* sorted order pointer		*/
#   struct fieldnode *	link;		/* linked field chain		*/
#   struct formnode *	form;		/* containing form		*/
#   struct typenode *	type;		/* field type			*/
#   void *		arg;		/* argument for type		*/
#   FIELD_CELL *		buf;		/* field buffers		*/
#   void *		usrptr;		/* user pointer			*/
#   /*
#    * The wide-character configuration requires extra information.  Because
#    * there are existing applications that manipulate the members of FIELD
#    * directly, we cannot make the struct opaque.  Offsets of members up to
#    * this point are the same in the narrow- and wide-character configuration.
#    * But note that the type of buf depends on the configuration, and is made
#    * opaque for that reason.
#    */
#   NCURSES_FIELD_INTERNALS
# } FIELD;
#
#
#   /*********
#   *  FORM  *
#   *********/
#
# typedef struct formnode {
#   unsigned short	status;	  	/* flags			*/
#   short			rows;		/* size in rows			*/
#   short			cols;		/* size in cols			*/
#   int			currow;		/* current row in field window	*/
#   int			curcol;		/* current col in field window	*/
#   int			toprow;		/* in scrollable field window	*/
#   int			begincol;	/* in horiz. scrollable field	*/
#   short			maxfield;	/* number of fields		*/
#   short			maxpage;	/* number of pages		*/
#   short			curpage;	/* index into page		*/
#   Form_Options		opts;		/* options			*/
#   WINDOW *		win;		/* window			*/
#   WINDOW *		sub;		/* subwindow			*/
#   WINDOW *		w;		/* window for current field	*/
#   FIELD **		field;		/* field [maxfield]		*/
#   FIELD *		current;	/* current field		*/
#   _PAGE *		page;		/* page [maxpage]		*/
#   void *		usrptr;		/* user pointer			*/
#
#   void			(*forminit)(struct formnode *);
#   void			(*formterm)(struct formnode *);
#   void			(*fieldinit)(struct formnode *);
#   void			(*fieldterm)(struct formnode *);
#
# } FORM;
#
#
#   /**************
#   *  FIELDTYPE  *
#   **************/
#
# typedef struct typenode {
#   unsigned short	status;			/* flags		    */
#   long			ref;			/* reference count	    */
#   struct typenode *	left;			/* ptr to operand for |     */
#   struct typenode *	right;			/* ptr to operand for |     */
#
#   void* (*makearg)(va_list *);			/* make fieldtype arg	    */
#   void* (*copyarg)(const void *);		/* copy fieldtype arg 	    */
#   void	(*freearg)(void *);			/* free fieldtype arg	    */
#
# #if NCURSES_INTEROP_FUNCS
#   union {
#     bool (*ofcheck)(FIELD *,const void *);	/* field validation	    */
#     bool (*gfcheck)(FORM*,FIELD *,const void*);	/* generic field validation */
#   } fieldcheck;
#   union {
#     bool (*occheck)(int,const void *);		/* character validation     */
#     bool (*gccheck)(int,FORM*,
#         FIELD*,const void*);        /* generic char validation  */
#   } charcheck;
#   union {
#     bool (*onext)(FIELD *,const void *);        /* enumerate next value     */
#     bool (*gnext)(FORM*,FIELD*,const void*);    /* generic enumerate next   */
#   } enum_next;
#   union {
#     bool (*oprev)(FIELD *,const void *);	/* enumerate prev value     */
#     bool (*gprev)(FORM*,FIELD*,const void*);    /* generic enumerate prev   */
#   } enum_prev;
#   void* (*genericarg)(void*);                   /* Alternate Arg method     */
# #else
#   bool	(*fcheck)(FIELD *,const void *);	/* field validation	*/
#   bool	(*ccheck)(int,const void *);		/* character validation */
#
#   bool	(*next)(FIELD *,const void *);		/* enumerate next value */
#   bool	(*prev)(FIELD *,const void *);		/* enumerate prev value */
# #endif
# } FIELDTYPE;
#
# typedef void (*Form_Hook)(FORM *);
#
#   /***************************
#   *  miscellaneous #defines  *
#   ***************************/
#
# /* field justification */
# #define NO_JUSTIFICATION	(0)
# #define JUSTIFY_LEFT		(1)
# #define JUSTIFY_CENTER		(2)
# #define JUSTIFY_RIGHT		(3)
#
# /* field options */
# #define O_VISIBLE		(0x0001U)
# #define O_ACTIVE		(0x0002U)
# #define O_PUBLIC		(0x0004U)
# #define O_EDIT			(0x0008U)
# #define O_WRAP			(0x0010U)
# #define O_BLANK			(0x0020U)
# #define O_AUTOSKIP		(0x0040U)
# #define O_NULLOK		(0x0080U)
# #define O_PASSOK		(0x0100U)
# #define O_STATIC		(0x0200U)
# #define O_DYNAMIC_JUSTIFY	(0x0400U)	/* ncurses extension	*/
#
# /* form options */
# #define O_NL_OVERLOAD		(0x0001U)
# #define O_BS_OVERLOAD		(0x0002U)
#
# /* form driver commands */
# #define REQ_NEXT_PAGE	 (KEY_MAX + 1)	/* move to next page		*/
# #define REQ_PREV_PAGE	 (KEY_MAX + 2)	/* move to previous page	*/
# #define REQ_FIRST_PAGE	 (KEY_MAX + 3)	/* move to first page		*/
# #define REQ_LAST_PAGE	 (KEY_MAX + 4)	/* move to last page		*/
#
# #define REQ_NEXT_FIELD	 (KEY_MAX + 5)	/* move to next field		*/
# #define REQ_PREV_FIELD	 (KEY_MAX + 6)	/* move to previous field	*/
# #define REQ_FIRST_FIELD	 (KEY_MAX + 7)	/* move to first field		*/
# #define REQ_LAST_FIELD	 (KEY_MAX + 8)	/* move to last field		*/
# #define REQ_SNEXT_FIELD	 (KEY_MAX + 9)	/* move to sorted next field	*/
# #define REQ_SPREV_FIELD	 (KEY_MAX + 10)	/* move to sorted prev field	*/
# #define REQ_SFIRST_FIELD (KEY_MAX + 11)	/* move to sorted first field	*/
# #define REQ_SLAST_FIELD	 (KEY_MAX + 12)	/* move to sorted last field	*/
# #define REQ_LEFT_FIELD	 (KEY_MAX + 13)	/* move to left to field	*/
# #define REQ_RIGHT_FIELD	 (KEY_MAX + 14)	/* move to right to field	*/
# #define REQ_UP_FIELD	 (KEY_MAX + 15)	/* move to up to field		*/
# #define REQ_DOWN_FIELD	 (KEY_MAX + 16)	/* move to down to field	*/
#
# #define REQ_NEXT_CHAR	 (KEY_MAX + 17)	/* move to next char in field	*/
# #define REQ_PREV_CHAR	 (KEY_MAX + 18)	/* move to prev char in field	*/
# #define REQ_NEXT_LINE	 (KEY_MAX + 19)	/* move to next line in field	*/
# #define REQ_PREV_LINE	 (KEY_MAX + 20)	/* move to prev line in field	*/
# #define REQ_NEXT_WORD	 (KEY_MAX + 21)	/* move to next word in field	*/
# #define REQ_PREV_WORD	 (KEY_MAX + 22)	/* move to prev word in field	*/
# #define REQ_BEG_FIELD	 (KEY_MAX + 23)	/* move to first char in field	*/
# #define REQ_END_FIELD	 (KEY_MAX + 24)	/* move after last char in fld	*/
# #define REQ_BEG_LINE	 (KEY_MAX + 25)	/* move to beginning of line	*/
# #define REQ_END_LINE	 (KEY_MAX + 26)	/* move after last char in line	*/
# #define REQ_LEFT_CHAR	 (KEY_MAX + 27)	/* move left in field		*/
# #define REQ_RIGHT_CHAR	 (KEY_MAX + 28)	/* move right in field		*/
# #define REQ_UP_CHAR	 (KEY_MAX + 29)	/* move up in field		*/
# #define REQ_DOWN_CHAR	 (KEY_MAX + 30)	/* move down in field		*/
#
# #define REQ_NEW_LINE	 (KEY_MAX + 31)	/* insert/overlay new line	*/
# #define REQ_INS_CHAR	 (KEY_MAX + 32)	/* insert blank char at cursor	*/
# #define REQ_INS_LINE	 (KEY_MAX + 33)	/* insert blank line at cursor	*/
# #define REQ_DEL_CHAR	 (KEY_MAX + 34)	/* delete char at cursor	*/
# #define REQ_DEL_PREV	 (KEY_MAX + 35)	/* delete char before cursor	*/
# #define REQ_DEL_LINE	 (KEY_MAX + 36)	/* delete line at cursor	*/
# #define REQ_DEL_WORD	 (KEY_MAX + 37)	/* delete word at cursor	*/
# #define REQ_CLR_EOL	 (KEY_MAX + 38)	/* clear to end of line		*/
# #define REQ_CLR_EOF	 (KEY_MAX + 39)	/* clear to end of field	*/
# #define REQ_CLR_FIELD	 (KEY_MAX + 40)	/* clear entire field		*/
# #define REQ_OVL_MODE	 (KEY_MAX + 41)	/* begin overlay mode		*/
# #define REQ_INS_MODE	 (KEY_MAX + 42)	/* begin insert mode		*/
# #define REQ_SCR_FLINE	 (KEY_MAX + 43)	/* scroll field forward a line	*/
# #define REQ_SCR_BLINE	 (KEY_MAX + 44)	/* scroll field backward a line	*/
# #define REQ_SCR_FPAGE	 (KEY_MAX + 45)	/* scroll field forward a page	*/
# #define REQ_SCR_BPAGE	 (KEY_MAX + 46)	/* scroll field backward a page	*/
# #define REQ_SCR_FHPAGE	 (KEY_MAX + 47) /* scroll field forward	 half page */
# #define REQ_SCR_BHPAGE	 (KEY_MAX + 48) /* scroll field backward half page */
# #define REQ_SCR_FCHAR	 (KEY_MAX + 49) /* horizontal scroll char	*/
# #define REQ_SCR_BCHAR	 (KEY_MAX + 50) /* horizontal scroll char	*/
# #define REQ_SCR_HFLINE	 (KEY_MAX + 51) /* horizontal scroll line	*/
# #define REQ_SCR_HBLINE	 (KEY_MAX + 52) /* horizontal scroll line	*/
# #define REQ_SCR_HFHALF	 (KEY_MAX + 53) /* horizontal scroll half line	*/
# #define REQ_SCR_HBHALF	 (KEY_MAX + 54) /* horizontal scroll half line	*/
#
# #define REQ_VALIDATION	 (KEY_MAX + 55)	/* validate field		*/
# #define REQ_NEXT_CHOICE	 (KEY_MAX + 56)	/* display next field choice	*/
# #define REQ_PREV_CHOICE	 (KEY_MAX + 57)	/* display prev field choice	*/
#
# #define MIN_FORM_COMMAND (KEY_MAX + 1)	/* used by form_driver		*/
# #define MAX_FORM_COMMAND (KEY_MAX + 57)	/* used by form_driver		*/
#
# #if defined(MAX_COMMAND)
# #  if (MAX_FORM_COMMAND > MAX_COMMAND)
# #    error Something is wrong -- MAX_FORM_COMMAND is greater than MAX_COMMAND
# #  elif (MAX_COMMAND != (KEY_MAX + 128))
# #    error Something is wrong -- MAX_COMMAND is already inconsistently defined.
# #  endif
# #else
# #  define MAX_COMMAND (KEY_MAX + 128)
# #endif
#
#   /*************************
#   *  standard field types  *
#   *************************/
# extern NCURSES_EXPORT_VAR(FIELDTYPE *) TYPE_ALPHA;
# extern NCURSES_EXPORT_VAR(FIELDTYPE *) TYPE_ALNUM;
# extern NCURSES_EXPORT_VAR(FIELDTYPE *) TYPE_ENUM;
# extern NCURSES_EXPORT_VAR(FIELDTYPE *) TYPE_INTEGER;
# extern NCURSES_EXPORT_VAR(FIELDTYPE *) TYPE_NUMERIC;
# extern NCURSES_EXPORT_VAR(FIELDTYPE *) TYPE_REGEXP;
#
#   /************************************
#   *  built-in additional field types  *
#   *  They are not defined in SVr4     *
#   ************************************/
# extern NCURSES_EXPORT_VAR(FIELDTYPE *) TYPE_IPV4;      /* Internet IP Version 4 address */
#
#   /***********************
#   *  FIELDTYPE routines  *
#   ***********************/
# extern NCURSES_EXPORT(FIELDTYPE *) new_fieldtype (
#         bool (* const field_check)(FIELD *,const void *),
#         bool (* const char_check)(int,const void *));
# extern NCURSES_EXPORT(FIELDTYPE *) link_fieldtype(
#         FIELDTYPE *, FIELDTYPE *);
#
# extern NCURSES_EXPORT(int)	free_fieldtype (FIELDTYPE *);
# extern NCURSES_EXPORT(int)	set_fieldtype_arg (FIELDTYPE *,
#         void * (* const make_arg)(va_list *),
#         void * (* const copy_arg)(const void *),
#         void (* const free_arg)(void *));
# extern NCURSES_EXPORT(int)	 set_fieldtype_choice (FIELDTYPE *,
#         bool (* const next_choice)(FIELD *,const void *),
#               bool (* const prev_choice)(FIELD *,const void *));
#
#   /*******************
#   *  FIELD routines  *
#   *******************/
# extern NCURSES_EXPORT(FIELD *)	new_field (int,int,int,int,int,int);
# extern NCURSES_EXPORT(FIELD *)	dup_field (FIELD *,int,int);
# extern NCURSES_EXPORT(FIELD *)	link_field (FIELD *,int,int);
#
# extern NCURSES_EXPORT(int)	free_field (FIELD *);
# extern NCURSES_EXPORT(int)	field_info (const FIELD *,int *,int *,int *,int *,int *,int *);
# extern NCURSES_EXPORT(int)	dynamic_field_info (const FIELD *,int *,int *,int *);
# extern NCURSES_EXPORT(int)	set_max_field ( FIELD *,int);
# extern NCURSES_EXPORT(int)	move_field (FIELD *,int,int);
# extern NCURSES_EXPORT(int)	set_field_type (FIELD *,FIELDTYPE *,...);
# extern NCURSES_EXPORT(int)	set_new_page (FIELD *,bool);
# extern NCURSES_EXPORT(int)	set_field_just (FIELD *,int);
# extern NCURSES_EXPORT(int)	field_just (const FIELD *);
# extern NCURSES_EXPORT(int)	set_field_fore (FIELD *,chtype);
# extern NCURSES_EXPORT(int)	set_field_back (FIELD *,chtype);
# extern NCURSES_EXPORT(int)	set_field_pad (FIELD *,int);
# extern NCURSES_EXPORT(int)	field_pad (const FIELD *);
# extern NCURSES_EXPORT(int)	set_field_buffer (FIELD *,int,const char *);
# extern NCURSES_EXPORT(int)	set_field_status (FIELD *,bool);
# extern NCURSES_EXPORT(int)	set_field_userptr (FIELD *, void *);
# extern NCURSES_EXPORT(int)	set_field_opts (FIELD *,Field_Options);
# extern NCURSES_EXPORT(int)	field_opts_on (FIELD *,Field_Options);
# extern NCURSES_EXPORT(int)	field_opts_off (FIELD *,Field_Options);
#
# extern NCURSES_EXPORT(chtype)	field_fore (const FIELD *);
# extern NCURSES_EXPORT(chtype)	field_back (const FIELD *);
#
# extern NCURSES_EXPORT(bool)	new_page (const FIELD *);
# extern NCURSES_EXPORT(bool)	field_status (const FIELD *);
#
# extern NCURSES_EXPORT(void *)	field_arg (const FIELD *);
#
# extern NCURSES_EXPORT(void *)	field_userptr (const FIELD *);
#
# extern NCURSES_EXPORT(FIELDTYPE *)	field_type (const FIELD *);
#
# extern NCURSES_EXPORT(char *)	field_buffer (const FIELD *,int);
#
# extern NCURSES_EXPORT(Field_Options)	field_opts (const FIELD *);
#
#   /******************
#   *  FORM routines  *
#   ******************/
#
# extern NCURSES_EXPORT(FORM *)	new_form (FIELD **);
#
# extern NCURSES_EXPORT(FIELD **)	form_fields (const FORM *);
# extern NCURSES_EXPORT(FIELD *)	current_field (const FORM *);
#
# extern NCURSES_EXPORT(WINDOW *)	form_win (const FORM *);
# extern NCURSES_EXPORT(WINDOW *)	form_sub (const FORM *);
#
# extern NCURSES_EXPORT(Form_Hook)	form_init (const FORM *);
# extern NCURSES_EXPORT(Form_Hook)	form_term (const FORM *);
# extern NCURSES_EXPORT(Form_Hook)	field_init (const FORM *);
# extern NCURSES_EXPORT(Form_Hook)	field_term (const FORM *);
#
# extern NCURSES_EXPORT(int)	free_form (FORM *);
# extern NCURSES_EXPORT(int)	set_form_fields (FORM *,FIELD **);
# extern NCURSES_EXPORT(int)	field_count (const FORM *);
# extern NCURSES_EXPORT(int)	set_form_win (FORM *,WINDOW *);
# extern NCURSES_EXPORT(int)	set_form_sub (FORM *,WINDOW *);
# extern NCURSES_EXPORT(int)	set_current_field (FORM *,FIELD *);
# extern NCURSES_EXPORT(int)	field_index (const FIELD *);
# extern NCURSES_EXPORT(int)	set_form_page (FORM *,int);
# extern NCURSES_EXPORT(int)	form_page (const FORM *);
# extern NCURSES_EXPORT(int)	scale_form (const FORM *,int *,int *);
# extern NCURSES_EXPORT(int)	set_form_init (FORM *,Form_Hook);
# extern NCURSES_EXPORT(int)	set_form_term (FORM *,Form_Hook);
# extern NCURSES_EXPORT(int)	set_field_init (FORM *,Form_Hook);
# extern NCURSES_EXPORT(int)	set_field_term (FORM *,Form_Hook);
# extern NCURSES_EXPORT(int)	post_form (FORM *);
# extern NCURSES_EXPORT(int)	unpost_form (FORM *);
# extern NCURSES_EXPORT(int)	pos_form_cursor (FORM *);
# extern NCURSES_EXPORT(int)	form_driver (FORM *,int);
# # if NCURSES_WIDECHAR
# extern NCURSES_EXPORT(int)	form_driver_w (FORM *,int,wchar_t);
# # endif
# extern NCURSES_EXPORT(int)	set_form_userptr (FORM *,void *);
# extern NCURSES_EXPORT(int)	set_form_opts (FORM *,Form_Options);
# extern NCURSES_EXPORT(int)	form_opts_on (FORM *,Form_Options);
# extern NCURSES_EXPORT(int)	form_opts_off (FORM *,Form_Options);
# extern NCURSES_EXPORT(int)	form_request_by_name (const char *);
#
# extern NCURSES_EXPORT(const char *)	form_request_name (int);
#
# extern NCURSES_EXPORT(void *)	form_userptr (const FORM *);
#
# extern NCURSES_EXPORT(Form_Options)	form_opts (const FORM *);
#
# extern NCURSES_EXPORT(bool)	data_ahead (const FORM *);
# extern NCURSES_EXPORT(bool)	data_behind (const FORM *);
#
# #if NCURSES_SP_FUNCS
# extern NCURSES_EXPORT(FORM *)	NCURSES_SP_NAME(new_form) (SCREEN*, FIELD **);
# #endif
#
# #ifdef __cplusplus
#   }
# #endif
# /* *INDENT-ON*/


text:
# function that fills the main queues of each department with all new incidents
# happening during the current minute "mnt"
function fill_queues!(sim_data::DataFrame,
                        incidents::DataFrame,
                        location_district::Dict,
                        incident_queue::Array{Union{Missing,Int},3},
                        queue_used::Array{Int64,2},
                        queue_change::Array{Bool,2},
                        cin::Int64,
                        mnt::Int64,
                        card_incidents::Int64,
                        card_priorities::Int64,
                        card_districts::Int64,
                        max_queue::Int64)
    ## queue_used: array that saves whether a queue is in use
    #    queue_used .= 1
    ## sort the incident queues after the incidents with the longest waiting time
    for p = 1:card_priorities
        for i = 1:card_districts
            if queue_change[i,p] == true
                incident_queue[:,i,p] = sort!(@view(incident_queue[:,i,p]))
                queue_change[i,p] = false
            end
        end
    end
    ##  fill the incident queue with all new cases arriving during the new minute "mnt"
    while cin <= card_incidents && sim_data[cin,:incident_minute] == mnt
        responsible_district = location_district[sim_data[cin,:location_responsible]]
        incident_priority = incidents[cin,:priority]
        if queue_used[responsible_district,incident_priority] > max_queue
            error("incident queue is to short in district",responsible_district,"!")
        else
            incident_queue[queue_used[responsible_district,incident_priority],responsible_district,incident_priority] = cin
            cin += 1
            queue_used[responsible_district,incident_priority] += 1
        end
    end
    return cin
end


text:
using Agents

mutable struct SchellingAgent <: AbstractAgent
  id::Int # The identifier number of the agent
  pos::Tuple{Int,Int} # The x, y location of the agent
  mood::Bool # whether the agent is happy in its node. (true = happy)
  group::Int # The group of the agent,
             # determines mood as it interacts with neighbors
end

space = Space((10,10), moore = true)

properties = Dict(:min_to_be_happy => 3)
schelling = ABM(SchellingAgent, space; properties = properties)

function initialize(;numagents=320, griddims=(20, 20), min_to_be_happy=3)
    space = Space(griddims, moore = true) # make a Moore grid
    properties = Dict(:min_to_be_happy => 3)
    model = ABM(SchellingAgent, space; properties=properties,
     scheduler = random_activation)
    # populate the model with agents, adding equal amount of the two types of agents
    # at random positions in the model
    for n in 1:numagents
        agent = SchellingAgent(n, (1,1), false, n < numagents/2 ? 1 : 2)
        add_agent_single!(agent, model)
    end
    return model
end

function agent_step!(agent, model)
    agent.mood == true && return # do nothing if already happy
    minhappy = model.properties[:min_to_be_happy]
    neighbor_cells = node_neighbors(agent, model)
    count_neighbors_same_group = 0
    # For each neighbor, get group and compare to current agent's group
    # and increment count_neighbors_same_group as appropriately.
    for neighbor_cell in neighbor_cells
        node_contents = get_node_contents(neighbor_cell, model)
        # Skip iteration if the node is empty.
        length(node_contents) == 0 && continue
        # Otherwise, get the first agent in the node...
        agent_id = node_contents[1]
        # ...and increment count_neighbors_same_group if the neighbor's group is
        # the same.
        neighbor_agent_group = model.agents[agent_id].group
        if neighbor_agent_group == agent.group
            count_neighbors_same_group += 1
        end
    end
    # After counting the neighbors, decide whether or not to move the agent.
    # If count_neighbors_same_group is at least the min_to_be_happy, set the
    # mood to true. Otherwise, move the agent to a random node.
    if count_neighbors_same_group ≥ minhappy
        agent.mood = true
    else
        move_agent_single!(agent, model)
    end
    return
end

# initialize the model with 370 agents on a 20 by 20 grid.
model = initialize()
step!(model, agent_step!)     # run the model one step
step!(model, agent_step!, 3)  # run the model 3 steps.

model = initialize()
# An array of Symbols for the agent fields that are to be collected.
properties = [:pos, :mood, :group]
# Specifies at which steps data should be collected.
n = 5  # number of time steps to run the simulation
when = 1:n  # At which steps to collect data
# Use the step function to run the model and collect data into a DataFrame.
data = step!(model, agent_step!, n, properties, when=when)
data[1:10, :] # print only a few rows

model = initialize(numagents=370, griddims=(20,20), min_to_be_happy=3);
properties = Dict(:mood => [sum])
n = 5; when = 1:n
data = step!(model, agent_step!, 5, properties, when=when)

# AGAIN
model = initialize()
# An array of Symbols for the agent fields that are to be collected.
properties = [:pos, :mood, :group]
# Specifies at which steps data should be collected.
n = 5  # number of time steps to run the simulation
when = 1:n  # At which steps to collect data
# Use the step function to run the model and collect data into a DataFrame.
data = step!(model, agent_step!, n, properties, when=when)
data[1:10, :] # print only a few rows

# Use the plot2D function from AgentsPlots.jl to plot distribution of agents at any step.
using AgentsPlots

#p = plot2D(data, :group, t=1, nodesize=10)
#p = plot2D(data, :group, t=2, nodesize=10)

for i in 1:2
  p = plot2D(data, :group, t=i, nodesize=10)
end

using DataVoyager
v = Voyager(data)

model = initialize(numagents=370, griddims=(20,20), min_to_be_happy=3);
data = step!(model, agent_step!, 5, properties, when=when, replicates=5)

using Distributed
addprocs(2)

data = step!(model, agent_step!, 2, properties,
             when=when, replicates=5, parallel=true)

paramscan

happyperc(moods) = count(x -> x == true, moods)/length(moods)

properties= Dict(:mood=>[happyperc])
parameters = Dict(:min_to_be_happy=>collect(2:5), :numagents=>[200,300], :griddims=>(20,20))

data = paramscan(parameters, initialize;
       properties=properties, n = 3, agent_step! = agent_step!)

data = paramscan(parameters, initialize;
       properties=properties, n = 3, agent_step! = agent_step!,
       replicates=3)

using DataFrames: Not, select!
using Statistics: mean
data_mean = Agents.aggregate(data, [:step, :min_to_be_happy, :numagents],  mean);
select!(data_mean, Not(:replicate_mean))


text:
"This is a demo grafter that has no effect"
function recurse_only_grafter(reflection)
    body = reflection.ast
    MacroTools.postwalk(body) do expr
        if isexpr(expr, :call)
            # TODO: make this ref `graft` right when in module
            # TODO: make this do kwargs right
            Expr(:block,
             #   :(println("Calling ", $(expr.args[1]))),
                recurse!(expr, reflection)
            )
        else
            expr
        end
    end
end

function translate_strings_grafter(reflection)
    ast = reflection.ast
    MacroTools.postwalk(ast) do expr
        if isexpr(expr, :call)
            return recurse!(expr, reflection)
        elseif isexpr(expr, String) ## String literal
            # Stub for actually translating the text
            return uppercase(expr)
        else
            return expr
        end
    end
end

@inline function fma_grafter(reflection)
    ast = reflection.ast
    MacroTools.prewalk(ast) do expr
        if (
            @capture(expr, x_ * y_ + z_) ||
            @capture(expr, z_ + x_ * y_)
            )

            return :(fma($x, $y, $z))
        elseif isexpr(expr, :call)
            return Expr(:call,
                :graft, fma_grafter, expr.args...
            )
        end
        return expr
    end
end


text:
using jInv.Mesh
using jInv.Utils
using jInv.ForwardShare
using ShapeReconstructionPaLS.Utils
using ParamLevelSet
using ShapeReconstructionPaLS.DipTransform;
using Statistics
using Distributed
using LinearAlgebra
using Test

println("Testing getData for all the methods:")
n = [32,32,32];
Mesh = getRegularMesh([0.0;3.0;0.0;3.0;0.0;3.0],n);
alpha = [1.5;2.5;2.0;-1.0];
beta = [1.0;2.0;-1.6;1.5];
Xs = [1.5 1.5 1.5 ; 2.0 2.0 2.0; 1.2 2.3 1.5; 2.2 1.5 2.0];
m = wrapTheta(alpha,beta,Xs);
sigmaH = getDefaultHeaviside();
u, = ParamLevelSetModelFunc(Mesh,m;computeJacobian=0,sigma = sigmaH);
samplingBinning = 2;

theta_phi = deg2rad.([37.0 24.0 ; 27.0 14.0;11.0 11.0]);
b = 4.0*[1.0 2.0 3.0 ; -2.3 1.2 2.5;-1.3 1.2 1.1]*mean(Mesh.h);
b[:] .= 0.0;

# RBFBased        = "RBFs";

methods = [MATBased,MATFree,RBFBasedSimple1,RBFBasedSimple2,RBFBased];
D = Array{Array{Float32,2}}(undef,length(methods));

for k=1:2
	pFor = getDipParam(Mesh,theta_phi,b,samplingBinning,nworkers(),methods[k]);
	Dremote, = getData(u[:],pFor);
	D[k] = arrangeRemoteCallDataIntoLocalData(Dremote);
end

for k=3:4
	pFor = getDipParam(Mesh,theta_phi,b,samplingBinning,nworkers(),methods[k]);
	Dremote, = getData(m[:],pFor);
	D[k] = arrangeRemoteCallDataIntoLocalData(Dremote);
end
k=5
m_wraped = wrapRBFparamAndRotationsTranslations(m,theta_phi,b)
pFor = getDipParam(Mesh,theta_phi,b,samplingBinning,nworkers(),methods[k]);
Dremote, = getData(m_wraped[:],pFor);
D[k] = arrangeRemoteCallDataIntoLocalData(Dremote);

println("All these should be approx 0")
for k=3:5
	r = norm(D[k] - D[2])./norm(D[2]);
	if r > 5e-1
		println("testGetData: Get data is not the same in all methods");
	end
	@test r < 5e-1;
end

println("Testing Sensitivities:");

du = 0.1*randn(size(u));
for k=1:2
	pFor = getDipParam(Mesh,theta_phi,b,samplingBinning,1,methods[k]);
	Dremote,pFor = getData(u[:],pFor);
	Dk = arrangeRemoteCallDataIntoLocalData(Dremote);

	ui = u; # not sure why this is needed. Otherwise julia do not recognise u in the loop
	for ii = 1:3
		dm = (0.5^ii)*du;
		println("Testing sensitivities for Method ",methods[k])
		pFor_t = getDipParam(Mesh,theta_phi,b,samplingBinning,1,methods[k]);

		Dremote, = getData(ui[:] .+ dm[:],pFor_t);
		Dk_t = arrangeRemoteCallDataIntoLocalData(Dremote);
		println("norm(dt-d0): ",norm(Dk_t[:].-Dk[:]),", norm(dt - d0 - J0*dm): ",norm(Dk_t[:] .- Dk[:] .- getSensMatVec(dm[:],u[:],fetch(pFor[1]))[:]));

		println("----------- Transpose test --------------");

		ui = reshape(ui,tuple(n...));
		u1 = copy(ui);
		d1 = sampleSlices(ui,0.01,1);
		t1 = randn(size(d1));
		sampleSlicesT(t1,0.01,u1);
		@test abs(dot(t1,d1) - dot(u1,ui)) < 1e-5

		Vd = randn(size(Dk_t));
		vm = randn(size(dm));
		val1 = dot(getSensMatVec(vm,ui[:],fetch(pFor[1]))[:],Vd);
		val2 = dot(getSensTMatVec(Vd[:],ui[:],fetch(pFor[1]))[:],vm);
		println(methods[k])
		@test abs(val1-val2) < 1e-5
	end
end
dmm = 0.05*randn(size(m));
for k=3:4
	println("Testing sensitivities for Method ",methods[k])
	pFor = getDipParam(Mesh,theta_phi,b,samplingBinning,1,methods[k]);
	Dremote, = getData(m[:],pFor);
	Dk = arrangeRemoteCallDataIntoLocalData(Dremote);
	println("Done!")

	for ii = 1:8
		dm = (0.5^ii)*dmm;
		pFor_t = getDipParam(Mesh,theta_phi,b,samplingBinning,1,methods[k]);
		Dremote, = getData(m[:] + dm[:],pFor_t);
		Dk_t = arrangeRemoteCallDataIntoLocalData(Dremote);
		println("norm(dt-d0): ",norm(Dk_t[:]-Dk[:]),", norm(dt - d0 - J0*dm): ",norm(Dk_t[:] - Dk[:] - getSensMatVec(dm[:],m[:],fetch(pFor[1]))[:]));
	end
end
k=5
println("Testing sensitivities for Method ",methods[k])
m_wraped = wrapRBFparamAndRotationsTranslations(m,theta_phi,b)
dm_wraped = wrapRBFparamAndRotationsTranslations(dmm,0.05*randn(size(theta_phi)),0.05*randn(size(b)));
pFor = getDipParam(Mesh,theta_phi,b,samplingBinning,nworkers(),methods[k]);
Dremote, = getData(m_wraped[:],pFor);
Dk = arrangeRemoteCallDataIntoLocalData(Dremote);
V = 0;
println("Done!")

for ii = 1:8
	dm = (0.5^ii)*dm_wraped;
	pFor_t = getDipParam(Mesh,theta_phi,b,samplingBinning,1,methods[k]);
	Dremote, = getData(m_wraped[:] + dm[:],pFor_t);
	Dk_t = arrangeRemoteCallDataIntoLocalData(Dremote);
	V = getSensMatVec(dm[:],m_wraped[:],fetch(pFor[1]))[:];
	println("norm(dt-d0): ",norm(Dk_t[:]-Dk[:]),", norm(dt - d0 - J0*dm): ",norm(Dk_t[:] - Dk[:] - V));
end

max_stars_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_issues_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_forks_repo_name: BGUCompSci/ShapeReconstructionPaLS


text:
# ALPX
#
# PART OF BAZINGA.jl

using Base.Iterators
using ProximalAlgorithms.IterationTools
using LinearAlgebra: norm
using OptiMo
using Printf

###########################################################################
# Options
###########################################################################
struct ALPX_options{R<:Real}
    tol_cviol::R
    tol_cviol_sqrt::R
    tol_optim_sqrt::R
    dual_reset::Bool
    dual_grow::R
    y_min::R
    y_max::R
    y_min_min::R
    y_max_max::R
    μ_min::R
    μ_max::R
    μ_est::R
    μ_down::R
    μ_grow::R
    μ_slow::R
    θ::R
    ϵ_min::R
    ϵ_init::R
    ϵ_down::R
    ϵ_rel::R
    max_sub_iter::Int
    subsolver::Symbol
    subsolver_verbose::Bool
    subsolver_freq::Int
    subsolver_gamma_min::R
    subsolver_tau_min::R
    bigeps::R
    warnings::Bool
end

###########################################################################
# Iterable
###########################################################################
struct ALPX_iterable{R<:Real,Tx<:AbstractArray{R},TP,TO<:ALPX_options{R}}
    prob::TP                        # problem
    x0::Tx                          # initial primal point
    y0::Tx                          # initial dual point
    μ::Maybe{Tx}                    # penalty parameters
    opts::TO                        # options
end

Base.IteratorSize(::Type{<:ALPX_iterable}) = Base.IsInfinite()

###########################################################################
# State
###########################################################################
mutable struct ALPX_state{R<:Real,Tx<:AbstractArray{R},TS,TA<:AugLagOptiModel{R}}
    x::Tx                   # primal iterate
    fx::R                   # value of smooth term
    gx::R                   # value of proximable term
    objec::R                # value of objective function
    cx::Tx                  # value of constraints
    y::Tx                   # dual iterate
    μ::Tx                   # penalty parameter vector
    optim::Maybe{R}         # optimality/criticality
    cviol::R                # constraint violation
    cslack::R               # compl. slackness
    # additional allocation
    px::Tx                  # constraint violation vector
    w::Tx                   # compl. slackness vector : cx - proj_S( cx + μ .* y )
    z::Tx                   # dual step vector : w ./ μ
    μ_zero::R               # slowly vanishing penalty parameter
    y_min::R                # lower dual safeguard
    y_max::R                # upper dual safeguard
    ϵsub::R                 # inner tolerance
    sx::Tx                  # infeasibility stationarity vector
    infsta::R               # infeasibility stationarity
    cslack_old::Maybe{R}    # previous compl. slackness
    sub_iter::Int           # subsolver iterations
    sub_flag::Bool          # subproblem successful
    sub_iter_tot::Int       # overall subsolver iterations
    subsolver::TS           # subsolver
    alprob::TA              # augmented Lagrangian subproblem
end

###########################################################################
# Initialize
###########################################################################
function Base.iterate(iter::ALPX_iterable{R}) where {R}
    # initial point
    x = copy(iter.x0)
    y = copy(iter.y0)
    gx, x = objprox(iter.prob, x, iter.opts.bigeps)
    fx = obj(iter.prob, x)
    objec = fx + gx
    y_min = iter.opts.y_min
    y_max = iter.opts.y_max
    dual_safeguard!(y, y_min, y_max, iter.opts)
    # constraint violation
    cx = cons(iter.prob, x)
    w = proj(iter.prob, cx) # tmp
    px = cx - w
    cviol = norm(px, Inf)
    # allocation
    z = similar(cx)
    sx = similar(x)
    # infeasibility stationarity
    infeas_station!(iter.prob, x, px, iter.opts.bigeps, sx)
    infsta = norm(sx, Inf)
    # penalty parameter vector [BiM12]
    μ_zero = copy(iter.opts.μ_max)
    μ = iter.μ
    if μ === nothing
        μ = similar(px)
        penalty_estimate!(objec, px, μ, iter.opts.μ_est)
    end
    μ .= max.(iter.opts.μ_min, min.(μ, iter.opts.μ_max))
    # complementarity slackness
    compl_slackness!(iter.prob, cx, μ, y, w, z)
    cslack = norm(w, Inf)
    # subproblem tolerance
    ϵsub = max(iter.opts.ϵ_min, iter.opts.ϵ_init)
    # subproblem solver
    if iter.opts.subsolver == :zerofpr
        subsolver = Bazinga.ZEROFPR(
            max_iter = iter.opts.max_sub_iter,
            verbose = iter.opts.subsolver_verbose,
            freq = iter.opts.subsolver_freq,
            gamma_min = iter.opts.subsolver_gamma_min,
            tau_min = iter.opts.subsolver_tau_min,
        )
    else
        subsolver = Bazinga.PANOC(
            max_iter = iter.opts.max_sub_iter,
            verbose = iter.opts.subsolver_verbose,
            freq = iter.opts.subsolver_freq,
            gamma_min = iter.opts.subsolver_gamma_min,
            tau_min = iter.opts.subsolver_tau_min,
        )
    end

    # augmented Lagrangian subproblem
    alprob = AugLagOptiModel(iter.prob, μ, y)

    state = ALPX_state(
        x,
        fx,
        gx,
        objec,
        cx,
        y,
        μ,
        nothing,
        cviol,
        cslack,
        px,
        w,
        z,
        μ_zero,
        y_min,
        y_max,
        ϵsub,
        sx,
        infsta,
        nothing,
        -1,
        false,
        0,
        subsolver,
        alprob,
    )

    return state, state
end

###########################################################################
# Step
###########################################################################
function Base.iterate(iter::ALPX_iterable{R}, state::ALPX_state{R,Tx}) where {R,Tx}
    # update the augmented Lagrangian subproblem
    AugLagUpdate!(state.alprob, state.μ, state.y)
    # update and call subsolver
    state.subsolver.tol_optim = state.ϵsub
    sub_out = state.subsolver(state.alprob, x0 = state.x)
    state.sub_iter = sub_out.iterations
    state.sub_iter_tot += state.sub_iter
    state.sub_flag = (sub_out.optimality <= state.ϵsub)
    # primal update
    copyto!(state.x, sub_out.x)
    # objective
    state.gx = sub_out.solver[:gx]
    state.fx = obj(iter.prob, state.x)
    state.objec = state.fx + state.gx
    # (outer) optimality coincides with subproblem optimality because of the
    # first-order dual update [LeV20]
    state.optim = sub_out.optimality
    # constraint violations
    cons!(iter.prob, state.x, state.cx)
    proj!(iter.prob, state.cx, state.w) # tmp
    state.px .= state.cx .- state.w
    state.cviol = norm(state.px, Inf)
    # infeasibility stationarity
    infeas_station!(iter.prob, state.x, state.px, iter.opts.bigeps, state.sx)
    state.infsta = norm(state.sx, Inf)
    # safeguarded dual update
    # complementarity slackness
    compl_slackness!(iter.prob, state.cx, state.μ, state.y, state.w, state.z)
    state.cslack = norm(state.w, Inf)
    state.y .+= state.z
    # penalty parameters
    state.μ_zero *= iter.opts.μ_slow
    if state.cslack_old === nothing
        if state.cviol > 0
            penalty_estimate!(state.objec, state.px, state.μ, iter.opts.μ_est)
        end
    else
        if (state.cviol <= iter.opts.tol_cviol)
            state.w .= state.cx .+ state.μ .* state.y
            proj!(iter.prob, state.w, state.z)
            state.w .= state.cx .- state.z
            if norm(state.w, Inf) <= iter.opts.tol_cviol
                state.μ *= iter.opts.μ_grow
                state.μ_zero = min(state.μ_zero, maximum(state.μ))
            elseif (state.cslack <= iter.opts.θ * state.cslack_old)
                # do nothing
            else
                state.μ .*= iter.opts.μ_down
            end
        elseif (state.cslack <= iter.opts.θ * state.cslack_old)
            # do nothing
        else
            state.μ .*= iter.opts.μ_down
        end
    end
    state.μ .= max.(iter.opts.μ_min, min.(state.μ, state.μ_zero))
    state.cslack_old = copy(state.cslack)
    dual_safeguard!(state.y, state.y_min, state.y_max, iter.opts)
    # subproblem tolerance
    if state.sub_flag &&
       (state.optim <= iter.opts.tol_optim_sqrt) &&
       (state.cslack <= iter.opts.tol_cviol_sqrt) &&
       (state.cviol <= iter.opts.tol_cviol_sqrt)
        state.ϵsub *= iter.opts.ϵ_down
        state.ϵsub = max(iter.opts.ϵ_min, min(state.ϵsub, iter.opts.ϵ_rel * state.optim))
    end

    return state, state
end

"""
    penalty_estimate!( objec, cviol, μ, scale )
Estimates suitable penalty parameters balancing objective `objec` and constraint
violations `cviol` [BiM12].
"""
function penalty_estimate!(
    objec::R,
    cviol::AbstractVector{R},
    μ::AbstractVector{R},
    scale::R,
) where {R<:Real}
    μ .= scale .* max.(1.0, 0.5 .* abs2.(cviol)) ./ max(1.0, abs(objec))
end

"""
    infeas_station!( prob, x, px, bigeps, sx )
Computes stationarity of the infeasibility measure.
"""
function infeas_station!(
    prob::TP,
    x::AbstractVector{R},
    px::AbstractVector{R},
    bigeps::R,
    sx::AbstractVector{R},
) where {R<:Real,TP}
    jtprod!(prob, x, px, sx) # tmp
    sx .= x .- bigeps .* sx # tmp
    prox!(prob, sx, bigeps, sx) # tmp
    sx .= (x .- sx) ./ bigeps
end

"""
    compl_slackness!( prob, cx, μ, y, w, z )
"""
function compl_slackness!(
    prob::TP,
    cx::AbstractVector{R},
    μ::AbstractVector{R},
    y::AbstractVector{R},
    w::AbstractVector{R},
    z::AbstractVector{R},
) where {R<:Real,TP}
    w .= cx .+ μ .* y
    proj!(prob, w, z)
    w .= cx .- z
    z .= w ./ μ
end

"""
    dual_safeguard!( y, ymin, ymax, opts )
If an entry of the dual estimate `y` is out of bounds `ymin` and `ymax`, projects
it, or reset it to zero, and expand bounds, within their limit [BiM14].
"""
function dual_safeguard!(
    y::AbstractVector{R},
    ymin::R,
    ymax::R,
    opts::TO,
) where {R<:Real,TO<:ALPX_options}
    idy = (y .<= ymin)
    if any(idy)
        y[idy] .= opts.dual_reset ? R(0) : ymin
        ymin *= opts.dual_grow
        ymin = max(opts.y_min_min, ymin)
    end
    idy = (y .>= ymax)
    if any(idy)
        y[idy] .= opts.dual_reset ? R(0) : ymax
        ymax *= opts.dual_grow
        ymax = min(opts.y_max_max, ymax)
    end
end

###########################################################################
# Solver
###########################################################################
struct ALPX{R<:Real}
    max_iter::Int
    max_sub_iter::Int
    max_sub_iter_tot::Int
    tol_optim::R
    tol_cviol::R
    tol_cviol_sqrt::R
    tol_inf_mu::R
    tol_diverging::R
    bigeps::R
    verbose::Bool
    freq::Int
    opts::ALPX_options{R}

    function ALPX{R}(;
        max_iter::Int = 100,
        max_sub_iter::Int = 10000,
        max_sub_iter_tot::Int = 100000,
        tol_optim::R = R(1e-8),
        tol_cviol::R = R(1e-8),
        theta::R = R(0.5),
        mu_down::R = R(0.25),
        eps_down::R = R(0.1),
        verbose::Bool = false,
        freq::Int = 1,
        dual_reset::Bool = true,
        dual_grow::R = R(1.1),
        mu_min::R = R(1e-12),
        mu_max::R = R(1e6),
        mu_est::R = R(0.1),
        mu_grow::R = R(10),
        mu_slow::R = R(0.99),
        y_min::R = R(-1e3),
        y_max::R = R(1e3),
        y_min_min::R = R(-1e6),
        y_max_max::R = R(1e6),
        eps_min::R = min(tol_optim, tol_cviol),
        eps_init::R = sqrt(tol_optim),
        eps_rel::R = R(0.5),
        tol_diverging::R = R(1e20),
        subsolver::Symbol = :zerofpr,
        subsolver_verbose::Bool = false,
        subsolver_freq::Int = 10,
        subsolver_gamma_min::R = 1e-16,
        subsolver_tau_min::R = 1e-8,
        bigeps::R = R(1e-12),
        warnings::Bool = true,
    ) where {R}
        @assert 0 < max_iter
        @assert 0 < max_sub_iter
        @assert 0 < max_sub_iter_tot
        @assert 0 < bigeps
        @assert 0 < tol_optim < 1 < tol_diverging
        @assert 0 < tol_cviol < 1
        @assert 0 < theta < 1
        @assert 0 < mu_down < 1
        @assert 0 < freq
        @assert 1 <= dual_grow
        @assert 0 < mu_min < 1
        @assert mu_min < mu_max
        @assert 0 < mu_est
        @assert 0 < mu_slow < 1 <= mu_grow
        @assert y_min_min <= y_min <= 0 <= y_max <= y_max_max
        @assert 0 < eps_min <= tol_optim
        @assert eps_min <= eps_init
        @assert 0 < eps_rel < 1
        @assert 0 < eps_down < 1
        @assert subsolver ∈ [:zerofpr, :panoc]
        @assert 0 < subsolver_freq
        @assert 0 < subsolver_gamma_min
        @assert 0 < subsolver_tau_min
        tol_optim_sqrt = sqrt(tol_optim)
        tol_cviol_sqrt = sqrt(tol_cviol)
        tol_inf_mu = cbrt(mu_min^2)
        opts = ALPX_options(
            tol_cviol,
            tol_cviol_sqrt,
            tol_optim_sqrt,
            dual_reset,
            dual_grow,
            y_min,
            y_max,
            y_min_min,
            y_max_max,
            mu_min,
            mu_max,
            mu_est,
            mu_down,
            mu_grow,
            mu_slow,
            theta,
            eps_min,
            eps_init,
            eps_down,
            eps_rel,
            max_sub_iter,
            subsolver,
            subsolver_verbose,
            subsolver_freq,
            subsolver_gamma_min,
            subsolver_tau_min,
            bigeps,
            warnings,
        )

        new(
            max_iter,
            max_sub_iter,
            max_sub_iter_tot,
            tol_optim,
            tol_cviol,
            tol_cviol_sqrt,
            tol_inf_mu,
            tol_diverging,
            bigeps,
            verbose,
            freq,
            opts,
        )
    end
end

###########################################################################
# Solver call
###########################################################################
function (solver::ALPX{R})(
    prob::TP;
    x0::AbstractArray{R} = prob.meta.x0,
    y0::AbstractArray{R} = prob.meta.y0,
    mu::Maybe{AbstractArray{R}} = nothing,
) where {R,TP<:AbstractOptiModel}

    tstart = time()

    stop_solved(state::ALPX_state) =
        (state.optim === nothing ? false : state.optim <= solver.tol_optim) &&
        (state.cslack === nothing ? false : state.cslack <= solver.tol_cviol) &&
        (state.cviol <= solver.tol_cviol)
    stop_unbound(state::ALPX_state) =
        (state.cviol <= solver.tol_cviol) &&
        (norm(state.x, Inf) > solver.tol_diverging || state.objec < -solver.tol_diverging)
    is_illegal(x) = (x === nothing ? false : (any(isnan.(x)) || any(isinf.(x))))
    stop_illegal(state::ALPX_state) =
        is_illegal(state.x) ||
        is_illegal(state.y) ||
        is_illegal(state.optim) ||
        is_illegal(state.cx)
    stop_infeas(state::ALPX_state) =
        state.infsta <= solver.tol_optim &&
        state.cviol > solver.tol_cviol_sqrt &&
        maximum(state.μ) <= solver.tol_inf_mu
    stop(state::ALPX_state) =
        stop_solved(state) ||
        stop_unbound(state) ||
        stop_illegal(state) ||
        stop_infeas(state) ||
        (state.sub_iter_tot >= solver.max_sub_iter_tot)

    disp((i, state)) = begin
        @printf("%5d | %+.3e | ", i-1, state.objec)
        state.optim === nothing ? @printf("%9s  ", "") : @printf("%.3e  ", state.optim)
        state.cslack === nothing ? @printf("%9s  ", "") :
            @printf("%.3e  ", state.cslack)
        @printf("%.3e | ", state.cviol)
        @printf("%.1e  %.1e  %.1e | ", norm(state.y, Inf), minimum(state.μ), state.ϵsub)
        if (state.sub_iter >= 0)
            state.sub_flag ? @printf("(+) ") : @printf("(-) ")
            @printf("%5d | ", state.sub_iter)
            #state.idf ? @printf("F | ") : (state.idp ? @printf("P | ") : @printf("U | "))
        end
        @printf("\n")
    end

    iter = ALPX_iterable(prob, x0, y0, mu, solver.opts)
    iter = halt(iter, stop)
    iter = take(iter, solver.max_iter + 1)
    iter = enumerate(iter)
    if solver.verbose
        @printf(
            "%5s | %10s | %9s  %9s  %9s | %7s  %7s  %7s | \n",
            "iter",
            "objec",
            "optim",
            "cslack",
            "cviol",
            "|y|",
            "|μ|",
            "sub-ϵ"
        )
        iter = sample(iter, solver.freq)
        iter = tee(iter, disp)
    end

    num_iters, state = loop(iter)

    status = if stop_solved(state)
        :first_order
    elseif stop_infeas(state)
        :infeasible
    elseif stop_unbound(state)
        :unbounded
    elseif stop_illegal(state)
        :ieee_nan_inf
    elseif (num_iters > solver.max_iter) || (state.sub_iter_tot >= solver.max_sub_iter_tot)
        :max_iter
    else
        :unknown
    end

    return OptiOutput(
        status,
        x = state.x,
        y = state.y,
        objective = state.objec,
        optimality = state.optim,
        cviolation = state.cviol,
        iterations = num_iters - 1,
        time = time() - tstart,
        solver_name = "ALPX",
        solver = Dict(:cslackness => state.cslack,
                      :infstation => state.infsta,
                      :sub_iterations => state.sub_iter_tot),
    )

end

###########################################################################
# Outer constructors
###########################################################################
"""
    ALPX([max_iter, tol_optim, tol_cviol, verbose, freq, ...])

Instantiate the ALPX algorithm [DeM21] for solving optimization problems
of the form

    minimize    f(x) + g(x)
    subject to  c(x) ∈ S

where `f` and `c` are smooth, `g` is proximable, and `S` is projectable. If
`solver = ALPX(args...)`, then the above problem is solved with

    solver( prob, [x0, y0] )

Optional keyword arguments:

* `max_iter::Integer` (default: `1000`), maximum number of iterations to perform.
* `tol_optim::Real` (default: `1e-8`), absolute tolerance on the dual residual.
* `tol_cviol::Real` (default: `1e-8`), absolute tolerance on the primal residual.
* `verbose::Bool` (default: `true`), whether or not to print information during the iterations.
* `freq::Integer` (default: `10`), frequency of verbosity.
* `grad_scaling::Bool`
* `scaling_grad_max::Real`
* `scaling_min_value::Real`

References:

[BiM12] Birgin, Martinez, ``Augmented Lagrangian method with nonmonotone penalty
        parameters for constrained optimization``, Computational Optimization and
        Applications, 51(3), (2012).
        DOI: 10.1007/s10589-011-9396-0

[BiM14] Birgin, Martinez, ``Practical Augmented Lagrangian Methods for Constrained
        Optimization``, SIAM, Philadelphia, (2014).
        ISBN: 978-1611973358

[LeV20] Leyffer, Vanaret, ``An Augmented Lagrangian Filter Method``
        Mathematical Methods of Operations Research, (2020).
        DOI: 10.1007/s00186-020-00713-x

[SFP20] Sopasakis, Fresk, Patrinos, ``OpEn: Code Generation for Embedded Nonconvex
        Optimization``, 21st IFAC World Congress Proceedings, (2020).
        arXiv: 2003.00292v1

[DeM21] De Marchi, ``Augmented Lagrangian Proximal solver for constrained structured
        optimization``, (2021).
"""
ALPX(::Type{R}; kwargs...) where {R} = ALPX{R}(; kwargs...)
ALPX(; kwargs...) = ALPX(Float64; kwargs...)


text:
workspace()
using CUDA
using CUDNN

function bench()
  x = CuArray(rand(Float32,30,30,3,10))
  #desc = CUDNN.ActivationDesc(CUDNN_ACTIVATION_SIGMOID)
  for i = 1:1000
    y = similar(x)
    CUDNN.activation!(CUDNN_ACTIVATION_SIGMOID, x, y)
  end
  sync()
end
@time bench()

CUDA.cugc()
x = rand(Float32,5,4,3,2) |> CuArray
Array(x)
w = CuArray(rand(Float32,2,2,3,4))
Array(w)
desc = ConvolutionDesc(Float32, (0,0), (1,1))
y = CUDNN.convolution(x, w, desc)
Array(y)

y = randn(10,5,1,1) |> CuArray
y = softmax(CUDNN_SOFTMAX_ACCURATE, CUDNN_SOFTMAX_MODE_INSTANCE, x)
yy = Array(y)

function bench()
  x = rand(Float32,100,100,10,10) |> CudaArray
  for i = 1:1000
    td = CUDNN.TensorDescriptor(x)
  end
end

@time bench()


text:
using LinearAlgebra;
using Printf;
using Random;
using BenchmarkTools;
using DelimitedFiles;

using BivMatFun;
import BivMatFun.mp;

include("experiments_common.jl")

f = (x,y,i,j) -> 1 ./ (x-y) ./ sqrt(x+y);

Random.seed!(123);

function run_test()

    if is_ci_test()
        # We skip this test in the the CI pipeline, as it does perform
        # any sanity check, it's only meant to test performances.
        return true;
    end

    success = true;
    ntests = 7;
    data = zeros(ntests, 6);

    # Set this to false to disable timings
    check_times = true;

    for i = 1 : ntests
        n = 2^(i-1) * 64;

        A = randn(n, n) + 1im * randn(n, n)
        B = randn(n, n) + 1im * randn(n, n)
        C = randn(n, n) + 1im * randn(n, n)

        # A, B, C, FF = construct_random_benchmark(f, n)
        # @time FF = evaluate_reference_solution(f, A, B, C)

        F, info = fun2m(f, A, B, C, method = BivMatFun.Diag);

        if check_times
            thp = @benchmark fun2m($f, $A, $B, $C, method = $BivMatFun.Diag);
            thp = median(thp.times) / 1e9;
        else
            thp = 0.0;
        end

        F2, info2 = fun2m(f, A, B, C, method = BivMatFun.DiagNoHp);

        if check_times
            tnohp = @benchmark fun2m($f, $A, $B, $C, method = $BivMatFun.DiagNoHp);
            tnohp = median(tnohp.times) / 1e9;
        else
            tnohp = 0.0;
        end

        XD = BivMatFun.diag_fun(f, A, B, C);
        if check_times
            tdiag = @benchmark BivMatFun.diag_fun($f, $A, $B, $C);
            tdiag = median(tdiag.times) / 1e9
        else
            tdiag = 0.0
        end

        # @btime BivMatFun.diag_fun($f, $A, $B, $C);
        reshp = 0.0 # norm(F - FF) / norm(FF)
        res = norm(F - F2) / norm(F) # norm(FF - F2) / norm(FF);
        res2 = norm(F - XD) / norm(F) # norm(XD - FF) / norm(FF)

        @printf("n = %d\n", n)
        @printf("DIAG_HP: time = %f, res = %e, nblocks = %d, %d\n", thp, reshp, info.nblocksA, info.nblocksB)
        @printf("DIAG_NOHP time = %f, Residual: %e, nblocks = %d, %d\n", tnohp, res, info2.nblocksA, info2.nblocksB);
        @printf("DIAG time = %f, Residual: %e\n", tdiag, res2);

        # data[i, 1, 1:4] = [ thp, info.nblocksA, info.nblocksB, reshp ]
        # data[i, 2, 1:4] = [ tnohp, info2.nblocksA, info2.nblocksB, res ]
        # data[i, 3, 1:4] = [ tdiag, 0, 0, res2 ]

        data[i, 1:6] = [ n, thp, tnohp, tdiag, info.nblocksA, info.nblocksB ]

        writedlm("exp3.dat", data, '\t')
    end

    return success;

end

run_test();

max_stars_repo_name: numpi/BivMatFun

max_issues_repo_name: numpi/BivMatFun

max_forks_repo_name: numpi/BivMatFun


text:
## location information

export DILocation

@checked struct DILocation <: MDNode
    ref::API.LLVMMetadataRef
end
register(DILocation, API.LLVMDILocationMetadataKind)

line(location::DILocation) = Int(API.LLVMDILocationGetLine(location))
column(location::DILocation) = Int(API.LLVMDILocationGetColumn(location))

function scope(location::DILocation)
    ref = API.LLVMDILocationGetScope(location)
    ref == C_NULL ? nothing : Metadata(ref)::DIScope
end

function inlined_at(location::DILocation)
    ref = API.LLVMDILocationGetInlinedAt(location)
    ref == C_NULL ? nothing : Metadata(ref)::DILocation
end

## nodes

export DINode

abstract type DINode <: MDNode end

## variables

export DIVariable

abstract type DIVariable <: DINode end

for var in (:Local, :Global)
    var_name = Symbol("DI$(var)Variable")
    var_kind = Symbol("LLVM$(var_name)MetadataKind")
    @eval begin
        @checked struct $var_name <: DIVariable
            ref::API.LLVMMetadataRef
        end
        register($var_name, API.$var_kind)
    end
end

function file(var::DIVariable)
    ref = API.LLVMDIVariableGetFile(var)
    ref == C_NULL ? nothing : Metadata(ref)::DIFile
end

line(var::DIVariable) = Int(API.LLVMDIVariableGetLine(var))

## scopes

export DIScope, DILocalScope

abstract type DIScope <: DINode end

file(scope::DIScope) = DIFile(API.LLVMDIScopeGetFile(scope))

function name(scope::DIScope)
    len = Ref{Cuint}()
    data = API.LLVMExtraDIScopeGetName(scope, len)
    unsafe_string(convert(Ptr{Int8}, data), len[])
end

abstract type DILocalScope <: DIScope end

## file

export DIFile

@checked struct DIFile <: DIScope
    ref::API.LLVMMetadataRef
end
register(DIFile, API.LLVMDIFileMetadataKind)

function directory(file::DIFile)
    len = Ref{Cuint}()
    data = API.LLVMDIFileGetDirectory(file, len)
    unsafe_string(convert(Ptr{Int8}, data), len[])
end

function filename(file::DIFile)
    len = Ref{Cuint}()
    data = API.LLVMDIFileGetFilename(file, len)
    unsafe_string(convert(Ptr{Int8}, data), len[])
end

function source(file::DIFile)
    len = Ref{Cuint}()
    data = API.LLVMDIFileGetSource(file, len)
    unsafe_string(convert(Ptr{Int8}, data), len[])
end

## type

export DIType

abstract type DIType <: DIScope end

for typ in (:Basic, :Derived, :Composite, :Subroutine)
    typ_name = Symbol("DI$(typ)Type")
    typ_kind = Symbol("LLVM$(typ_name)MetadataKind")
    @eval begin
        @checked struct $typ_name <: DIType
            ref::API.LLVMMetadataRef
        end
        register($typ_name, API.$typ_kind)
    end
end

function name(typ::DIType)
    len = Ref{Csize_t}()
    data = API.LLVMDITypeGetName(typ, len)
    unsafe_string(convert(Ptr{Int8}, data), len[])
end

Base.sizeof(typ::DIType) = 8*Int(API.LLVMDITypeGetSizeInBits(typ))
offset(typ::DIType) = Int(API.LLVMDITypeGetOffsetInBits(typ))
line(typ::DIType) = Int(API.LLVMDITypeGetLine(typ))
flags(typ::DIType) = API.LLVMDITypeGetFlags(typ)

## subprogram

export DISubProgram

@checked struct DISubProgram <: DIScope
    ref::API.LLVMMetadataRef
end
register(DISubProgram, API.LLVMDISubprogramMetadataKind)

line(subprogram::DISubProgram) = Int(API.LLVMDISubprogramGetLine(subprogram))

## compile unit

export DICompileUnit

@checked struct DICompileUnit <: DIScope
    ref::API.LLVMMetadataRef
end
register(DICompileUnit, API.LLVMDICompileUnitMetadataKind)

## other

export DEBUG_METADATA_VERSION, strip_debuginfo!

DEBUG_METADATA_VERSION() = API.LLVMDebugMetadataVersion()

strip_debuginfo!(mod::Module) = API.LLVMStripModuleDebugInfo(mod)

function get_subprogram(func::Function)
    ref = API.LLVMGetSubprogram(func)
    ref==C_NULL ? nothing : Metadata(ref)::DISubProgram
end
set_subprogram!(func::Function, sp::DISubProgram) = API.LLVMSetSubprogram(func, sp)


text:
#=
focus-chart.jl
generate "focus chart" phantom from triangles
2021-09-12, Jeff Fessler, University of Michigan
=#

export focus_chart

using ImagePhantoms: Triangle, Object2d, RealU

"""
    focus_chart( ; radius=1, nspoke=30, value=1)

Generate `nspoke` Triangle phantom parameters for a focus chart.
Returns `Vector{Object2d{Triangle}` for passing to `phantom`.

# Options
- `radius::RealU = 1` radius of phantom
- `nspoke::Int = 60` # of spokes
- `value::Number = 1` alternate between 0 and this value
"""
focus_chart( ;
    nspoke::Int = 60,
    radius::RealU = 1,
    value::Number = 1,
    values::AbstractVector = value * (iseven.(1:nspoke)*2 .- 1), # alternating
) = focus_chart.(1:nspoke, nspoke, radius, values)

# ith of nspoke triangles
function focus_chart(i::Int, nspoke::Int, radius::RealU, value::Number)
    ϕ = (i-1)/nspoke * 2π
    (s, c) = sincos(ϕ)
    θ = π/nspoke
    (w,h) = sincos(θ) .* radius .* (2, 2/sqrt(3))
    offset = cos(θ) * radius + eps() # trick to avoid center pile up
    return Triangle((s, -c) .* offset, (w,h), ϕ, value)
end