text:
# NOTE: Heatmapping assumes Flux's WHCN convention (width, height, color channels, batch size).

const HEATMAPPING_PRESETS = Dict{Symbol,Tuple{ColorScheme,Symbol,Symbol}}(
    # Analyzer => (colorscheme, reduce, rangescale)
    :LRP => (ColorSchemes.bwr, :sum, :centered), # attribution
    :InputTimesGradient => (ColorSchemes.bwr, :sum, :centered), # attribution
    :Gradient => (ColorSchemes.grays, :norm, :extrema), # gradient
)

"""
    heatmap(expl::Explanation; kwargs...)
    heatmap(attr::AbstractArray; kwargs...)

    heatmap(input, analyzer::AbstractXAIMethod)
    heatmap(input, analyzer::AbstractXAIMethod, neuron_selection::Int)

Visualize explanation.
Assumes Flux's WHCN convention (width, height, color channels, batch size).

## Keyword arguments
- `cs::ColorScheme`: ColorScheme that is applied.
    When calling `heatmap` with an `Explanation` or analyzer, the method default is selected.
    When calling `heatmap` with an array, the default is `ColorSchemes.bwr`.
- `reduce::Symbol`: How the color channels are reduced to a single number to apply a colorscheme.
    The following methods can be selected, which are then applied over the color channels
    for each "pixel" in the attribution:
    - `:sum`: sum up color channels
    - `:norm`: compute 2-norm over the color channels
    - `:maxabs`: compute `maximum(abs, x)` over the color channels in
    When calling `heatmap` with an `Explanation` or analyzer, the method default is selected.
    When calling `heatmap` with an array, the default is `:sum`.
- `rangescale::Symbol`: How the color channel reduced heatmap is normalized before the colorscheme is applied.
    Can be either `:extrema` or `:centered`.
    When calling `heatmap` with an `Explanation` or analyzer, the method default is selected.
    When calling `heatmap` with an array, the default for use with the `bwr` colorscheme is `:centered`.
- `permute::Bool`: Whether to flip W&H input channels. Default is `true`.
- `unpack_singleton::Bool`: When heatmapping a batch with a single sample, setting `unpack_singleton=true`
    will return an image instead of an Vector containing a single image.

**Note:** these keyword arguments can't be used when calling `heatmap` with an analyzer.
"""
function heatmap(
    attr::AbstractArray{T,N};
    cs::ColorScheme=ColorSchemes.bwr,
    reduce::Symbol=:sum,
    rangescale::Symbol=:centered,
    permute::Bool=true,
    unpack_singleton::Bool=true,
) where {T,N}
    N != 4 && throw(
        DomainError(
            N,
            """heatmap assumes Flux's WHCN convention (width, height, color channels, batch size) for the input.
            Please reshape your attribution to match this format if your model doesn't adhere to this convention.""",
        ),
    )
    if unpack_singleton && size(attr, 4) == 1
        return _heatmap(attr[:, :, :, 1], cs, reduce, rangescale, permute)
    end
    return map(a -> _heatmap(a, cs, reduce, rangescale, permute), eachslice(attr; dims=4))
end

# Use HEATMAPPING_PRESETS for default kwargs when dispatching on Explanation
function heatmap(expl::Explanation; permute::Bool=true, kwargs...)
    _cs, _reduce, _rangescale = HEATMAPPING_PRESETS[expl.analyzer]
    return heatmap(
        expl.attribution;
        reduce=get(kwargs, :reduce, _reduce),
        rangescale=get(kwargs, :rangescale, _rangescale),
        cs=get(kwargs, :cs, _cs),
        permute=permute,
    )
end
# Analyze & heatmap in one go
function heatmap(input, analyzer::AbstractXAIMethod, args...; kwargs...)
    return heatmap(analyze(input, analyzer, args...; kwargs...))
end

# Lower level function that is mapped along batch dimension
function _heatmap(
    attr::AbstractArray{T,3},
    cs::ColorScheme,
    reduce::Symbol,
    rangescale::Symbol,
    permute::Bool,
) where {T<:Real}
    img = dropdims(_reduce(attr, reduce); dims=3)
    permute && (img = permutedims(img))
    return ColorSchemes.get(cs, img, rangescale)
end

# Reduce attributions across color channels into a single scalar – assumes WHCN convention
function _reduce(attr::AbstractArray{T,3}, method::Symbol) where {T}
    if size(attr, 3) == 1 # nothing to reduce
        return attr
    elseif method == :sum
        return reduce(+, attr; dims=3)
    elseif method == :maxabs
        return reduce((c...) -> maximum(abs.(c)), attr; dims=3, init=zero(T))
    elseif method == :norm
        return reduce((c...) -> sqrt(sum(c .^ 2)), attr; dims=3, init=zero(T))
    end
    throw(
        ArgumentError(
            "Color channel reducer :$method not supported, `reduce` should be :maxabs, :sum or :norm",
        ),
    )
end


text:
path = "/home/user/GraphEvolve.jl/src"
# path = "/home/perot/julia/GraphEvolve.jl/src";
push!(LOAD_PATH, path);
using GraphEvolve;
using Plots; pyplot(fmt="png");
PyPlot.matplotlib.rc("mathtext", fontset="cm");
PyPlot.matplotlib.rc("text", usetex=true);
PyPlot.matplotlib.rc("font", family="serif", size=12)
using LaTeXStrings
using Printf
using Dates

# %%

function plot_Network(Ns, evolution_method, savepath)
	colors = RGB[cgrad(:plasma)[z] for z=range(0, stop=1, length=length(Ns))]
	reverse!(colors)
	markers = filter((m->begin m in Plots.supported_markers() end), Plots._shape_keys)

	plot_C = plot(dpi=300,
		legend=:topleft,
		xaxis=(L"t/N", (0, 1.5), 0:0.5:1.5),
		yaxis=(latexstring("|C|/N"), (0, 1), 0:0.2:1))
	plot_S = plot(dpi=300,
		legend=false,
		xaxis=(L"t/N", (0, 1.5), 0:0.5:1.5),
		yaxis=(L"\overline{s}"), (0, 10))
	plot_H = plot(dpi=300,
		legend=false,
		xaxis=(L"t/N", (0, 1.5), 0:0.5:1.5),
		yaxis=(L"H/N"), (0, 0.004))

	for (i, N) in enumerate(Ns)
		t₀ = now()
		g = Network(N)
		evolution_method(g, Int(floor(1.5N)))

		n_string = @sprintf("%.1e", N)
		k = Int(log2(N))
		x = collect(0:g.t) ./ N

		scatter!(plot_C, x, g.observables.largest_cluster_size ./ N,
			label=latexstring("N = 2^{$(k)}"),
			marker=(2, colors[i], markers[i], 0.8, stroke(colors[i])))
		scatter!(plot_S, x, g.observables.avg_cluster_size,
			label=latexstring("N = 2^{$(k)}"),
			marker=(2, colors[i], markers[i], 0.8, stroke(colors[i])))
		scatter!(plot_H, x, g.observables.heterogeneity ./ N,
			label=latexstring("N = 2^{$(k)}"),
			marker=(2, colors[i], markers[i], 0.8, stroke(colors[i])))

		t₁ = now()
		runtime = Dates.value(t₁ - t₀) / 1000
		println("k = $(k), time to run: $(runtime)s")
	end

	savefig(plot_C, joinpath(savepath, "k_scaling_m.png"))

	l = @layout [a{0.6h}; [b{0.5w} c]]
	plot_ = plot(plot_C, plot_H, plot_S, layout=l, size=(800, 600))
	savefig(plot_, joinpath(savepath, "k_scaling_triple.png"))
end

function plot_Lattice2D(Ls, evolution_method, savepath)
	colors = RGB[cgrad(:plasma)[z] for z=range(0, stop=1, length=length(Ls))]
	reverse!(colors)
	markers = filter((m->begin m in Plots.supported_markers() end), Plots._shape_keys)

	plot_C = plot(dpi=300,
		legend=:topleft,
		xaxis=(L"t/N", (0, 1.5), 0:0.5:1.5),
		yaxis=(latexstring("|C|/N"), (0, 1), 0:0.2:1))
	plot_S = plot(dpi=300,
		legend=false,
		xaxis=(L"t/N", (0, 1.5), 0:0.5:1.5),
		yaxis=(L"\overline{s}"), (0, 10))
	plot_H = plot(dpi=300,
		legend=false,
		xaxis=(L"t/N", (0, 1.5), 0:0.5:1.5),
		yaxis=(L"H/N"), (0, 0.004))

	for (i, L) in enumerate(Ls)
		N = L^2
		g = Lattice2D(L)
		evolution_method(g, Int(floor(1.5N)))

		x = collect(0:g.t) ./ N
		scatter!(plot_C, x, g.observables.largest_cluster_size ./ N,
			label=latexstring("\$N = $(L)\$"),
			marker=(2, colors[i], markers[i], stroke(0)))
		scatter!(plot_S, x, g.observables.avg_cluster_size,
			label=latexstring("\$N = $(L)\$"),
			marker=(2, colors[i], markers[i], stroke(0)))
		scatter!(plot_H, x, g.observables.heterogeneity ./ N,
			label=latexstring("\$N = $(L)\$"),
			marker=(2, colors[i], markers[i], stroke(0)))
	end

	# savefig(plot_C, joinpath(savepath, "Lattice2D_C.png"))
	# savefig(plot_S, joinpath(savepath, "Lattice2D_S.png"))
	# savefig(plot_H, joinpath(savepath, "Lattice2D_H.png"))

	l = @layout [a{0.6h}; [b{0.5w} c]]
	plot_ = plot(plot_C, plot_H, plot_S, layout=l)
	savefig(plot_, joinpath(savepath, "Lattice2D_observables.png"))
end

function plot_Lattice3D(Ls, evolution_method, savepath)
	colors = RGB[cgrad(:plasma)[z] for z=range(0, stop=1, length=length(Ls))]
	reverse!(colors)
	markers = filter((m->begin m in Plots.supported_markers() end), Plots._shape_keys)

	plot_C = plot(dpi=300,
		legend=:topleft,
		xaxis=(L"t/N", (0, 1.5), 0:0.5:1.5),
		yaxis=(latexstring("|C|/N"), (0, 1), 0:0.2:1))
	plot_S = plot(dpi=300,
		legend=false,
		xaxis=(L"t/N", (0, 1.5), 0:0.5:1.5),
		yaxis=(L"\overline{s}"), (0, 10))
	plot_H = plot(dpi=300,
		legend=false,
		xaxis=(L"t/N", (0, 1.5), 0:0.5:1.5),
		yaxis=(L"H/N"), (0, 0.004))

	for (i, L) in enumerate(Ls)
		N = L^3
		g = Lattice3D(L)
		evolution_method(g, Int(floor(1.5N)))

		x = collect(0:g.t) ./ N
		scatter!(plot_C, x, g.observables.largest_cluster_size ./ N,
			label=latexstring("\$L = $(L)\$"),
			marker=(2, colors[i], markers[i], stroke(0)))
		scatter!(plot_S, x, g.observables.avg_cluster_size,
			label=latexstring("\$L = $(L)\$"),
			marker=(2, colors[i], markers[i], stroke(0)))
		scatter!(plot_H, x, g.observables.heterogeneity ./ N,
			label=latexstring("\$L = $(L)\$"),
			marker=(2, colors[i], markers[i], stroke(0)))
	end

	# savefig(plot_C, joinpath(savepath, "Lattice3D_C.png"))
	# savefig(plot_S, joinpath(savepath, "Lattice3D_S.png"))
	# savefig(plot_H, joinpath(savepath, "Lattice3D_H.png"))

	l = @layout [a{0.6h}; [b{0.5w} c]]
	plot_ = plot(plot_C, plot_H, plot_S, layout=l)
	savefig(plot_, joinpath(savepath, "Lattice3D_observables.png"))
end

# %%

evolution_method = stochastic_edge_acceptance!
Ns = 2 .^ collect(15:23);
Ls_2D = Int.(floor.(collect(range(200, stop=8000, length=11))));
Ls_3D = Int.(floor.(collect(range(40, stop=400, length=11))));
savepath = "/home/user/thesis/latex/images"
# savepath = "/home/perot/julia/scaling"

plot_Network(Ns, evolution_method, savepath)
# plot_Lattice2D(Ls_2D, evolution_method, savepath)
# plot_Lattice3D(Ls_3D, evolution_method, savepath)

ext: jl

max_stars_repo_name: cameronperot/explosive-percolation

max_issues_repo_name: cameronperot/explosive-percolation

max_forks_repo_name: cameronperot/explosive-percolation


text:
module AbstractSDRsFMReceiver

using AbstractSDRs;
using DSP;
using FFTW;

include("Audio.jl");
using .Audio

include("GetSpectrum.jl")
using .GetSpectrum

using Blink;        # For GUI in electron
using Interact;      # For Widget and interactions
using PlotlyJS;

FLAG = false;

"""
Apply FM demodulation to the input signal
It applies Arg(sig[n+1] sig*[n]) where * stands for complex conjugate
# --- Syntax
out = fmDemod(sig)
# --- Input parameters
- sig : Input signal (expected to be complex Vector)
# --- Output parameters
- out : Demodulated signal (Float64 Vector)
"""
function fmDemod(sig)
    out = zeros(Float64,length(sig));
    @inbounds @simd for n ∈ (1:length(sig)-1)
        out[n+1] = angle(sig[n+1]*conj(sig[n]));
    end
    return out;
end

"""
Init a Low pass filter (based on FIR filter synthesis) that cut @bandPass
# --- Syntax
h = genFilter(sizeL,bandPass,samplingRate)
# --- Input parameters
- sizeL : Desired filter size
- bandPass : Max frequency in bandpass
- samplingRate : sampling frequency
# --- Output parameters
- h : filter (Vector{Float64}(sizeL,1))
"""
function genFilter(sizeL,bandPass,samplingRate)
    # --- Init frequency response
    cuttOff = bandPass / samplingRate;
    fRep  = zeros(Complex{Float64},sizeL);
    # --- Set passband value
    fRep[1:Int(floor(sizeL*cuttOff))] .= 1;
    # --- Adding linear phase term
    pulsation   = 2*pi*(0:sizeL-1) ./sizeL;# --- Pulsation \Omega
    groupDelay  = - (sizeL-1)/2;# --- Group delay
    fRep        .= fRep .* exp.(1im * groupDelay .* pulsation);
    # --- Get to time domain
    # With a force to cast in real domain
    hRep  = real(ifft(fRep));
    # --- Generates window
    wind  = hamming(sizeL);
    # --- Apply windowing (apodisation)
    h  = hRep .* wind;
end

"""
Apply the filter by direct convoltuion and downsample
# --- Syntax
 out = applyFilter(sig,h,decim)
# --- Input parameters
- sig : Signal to be filtered and downsampled (of size N)
- h : Filter to apply
- decim : Desired decimation factor (should be Int)
# --- Output parameters
- out : filtered and downsampled data (size N÷decim)
"""
function applyFilter(x,h,decim)
    l = length(h)÷2;
    return conv(x,h)[1+l:decim:end-l];
end

"""
Compute Power Spectral Density (PSD) of signal and returns a tuple (frequency,power) ready to be computed
# --- Syntax
(fs,pow) = getSpectrum(f,sig,N)
# --- Input parameters
- f : sampling frequency (in Hz)
- sig : Signal to be analyzed
# --- Output parameters
- fs : Vector of frequencies
- pow: Vector of power
"""
function getSpectrum(fs,sig;N=nothing)
    if isnothing(N)
        N = length(sig);
    end
    freqAx = collect(((0:N-1)./N .- 0.5)*fs);
    ss	   = @view sig[1:N];
    y	   = 10*log10.(abs2.(fftshift(fft(ss))));
    return (freqAx,y);
end
getSpectrum(sig) = getSpectrum(1,sig);

function main(sdr,carrierFreq,gain=20,p=nothing;kwargs...)
    # --- Simulation parameters
    global FLAG = true;
    global newCarrierFreq	  = carrierFreq;		# --- Starting frequency
    samplingRate  			  = 48e3*4;
    # --- Duration and buffer size
    duration                  = 0.5;
    nbSamples                 = Int( duration * samplingRate);
    # --- Update radio configuration
    global radio			= openSDR(sdr,newCarrierFreq,samplingRate,gain;kwargs...);
    # --- Create FIR filter
    audioRendering = 48e3;
    decim          = Int(samplingRate ÷ audioRendering);
    h              = genFilter(128,48e3,samplingRate);
    # --- Audio rendering
    stream = Audio.createAudioBuffer();
    try
        while(true)
            # --- We get buffer
            sig = recv(radio,nbSamples);
            # --- FM demodulator
            out = fmDemod(sig);
            # --- Audio decimation
            audio = applyFilter(out,h,decim);
            # --- Audio rendering
            Audio.play(stream,audio);
            # --- Plot
            deletetraces!(p,0);
            (x,y) = getSpectrum(samplingRate,sig[1:128]);
            pl1	  = scatter(; x,y, name=" ");
            addtraces!(p,pl1);
            # Exit
            if FLAG == false;
                print("BREAK");
                break;
            end
        end
    catch exception;
        close(stream);
        close(radio);
        rethrow(exception);
    end
    close(stream);
    close(radio);
end

"""
Main call to the Graphical User Interace (GUI)
# --- Syntax
gui()
# --- Input parameters
-
# --- Output parameters
-
"""
function gui()
    global FLAG = false;
    # ----------------------------------------------------
    # --- Define widgets
    # ----------------------------------------------------
    # --- Creating widget for FM frequencies
    frequencies = widget(88:0.1:107, label="FM frequencies (MHz)");
    widgetGain        = widget(0:50, label="Radio Gain");
    # --- Creating some stuff related to args
    widgetArgs = textbox(hint=""; value="addr=192.168.10.13");
    # --- Widget for radio configuration
    sdrList = AbstractSDRs.getSupportedSDRs();
    options = Observable(sdrList);
    wdg = dropdown(options);
    # --- Start button
    startButton = button("Start !"; value=0)
    # ----------------------------------------------------
    # ---
    w = Window(Blink.@d(:width=>600, :height=>600));
    title(w,"FM receiver in Julia");
    # --- Horizontal axis
    layout = hbox(pad(1em,wdg),pad(1em,widgetArgs));
    layout = vbox(layout,frequencies);
    layout = vbox(layout,widgetGain);
    layout = hbox(layout,pad(3em,startButton));
    #
    plotLayout = Layout(;
        width = 600,
        height = 400,
        title=" Spectrum ",
        xaxis_title=" Frequency  ",
        yaxis_title=" Power ",
        xaxis_showgrid=true, yaxis_showgrid=true,
        yaxis_range =[-80,0]
        )
    N = 100;
    x = 0:N-1;
    y = zeros(N);
    pl1	  = scatter(; x,y, name=" ");
    p = plot(pl1,plotLayout)
    layout = vbox(layout,p);

    body!(w,layout);
    # ----------------------------------------------------
    # ---
    # ----------------------------------------------------
    # --- Update Frequency
    on(updateC,frequencies);
    # --- Update Gain
    on(updateG,widgetGain);
    # --- Starting routine
    h = on(startButton) do click
        startButton["is-loading"][]=true
        #@show click
        #if click == 1
        #    map!(fromStartToStop,startButton,12,16)
        #end
        sdr =  wdg[];
        carrierFreq = frequencies[]*1e6
        gain        = widgetGain[];
        kwargs = (;args=widgetArgs[]);
        @async main(sdr,carrierFreq,gain,p;kwargs...);
        # @async main(:e310,carrierFreq;args="addr=192.168.10.13");
    end
    success(w.shell.proc)
    global FLAG=false;
    close(radio);
    # We get the hand back and can close the radio
    return nothing
end

function updateC(val)
    if FLAG == true
        carrierFreq = val*1e6;
        updateCarrierFreq!(radio,carrierFreq);
    end
end
function updateG(val)
    if FLAG == true
        gain = val;
        updateGain!(radio,gain);
    end
end

    end # module


text:
module Extras

include("loopinfo.jl")
using .LoopInfo
export @unroll

end # module


text:
##Simple Plots of MAD UB bound
using DelimitedFiles
include("src.jl")  #load up vopp

function simplePlot(file_out_name;
					S=4, M=1, N=100)
	#Create file_out and header
	outPath = "$(file_out_name)__$(S)_$(M)_$(N)"
	f = open("$(outPath).tab", "w")
	writedlm(f, ["D" "Bound"])

	for  D in range(1e-10, stop=vopp.delta_h(S) - 1e-10, length=N)
		#add a dumb bound for ease in plotting
		bound = vopp.vopp_ub_MAD(S, 1., D)
		writedlm(f, [D bound])
	end
	close(f)
end

function contour_plot(file_out_name;
				Smin = 1.1, Smax=4, N=100, TOL = 1e-10)
	#Create file_out and header
	outPath = "$(file_out_name)__$(Smin)_$(Smax)_$(N)"
	f = open("$(outPath).tab", "w")
	writedlm(f, ["S" "D" "Bound"])

	S_grid = range(Smin, stop=Smax, length=N)
	D_grid = range(1e-10, stop=1, length=N)

	for S in S_grid
		for D in D_grid
			if D > vopp.delta_h(S) - TOL
				continue
			end
			bound = vopp.vopp_ub_MAD(S, 1., D)
			writedlm(f, [S D bound])
		end
		flush(f)
	end
	close(f)
end

function scale_bound_dist(file_out_name; S=4, N=100)
	#Create file_out and header
	outPath = "$(file_out_name)__$(S)_$(N)"
	f = open("$(outPath).tab", "w")
	writedlm(f, ["x" "Fbar"])

	x_grid = collect(range(0., stop=S, length=N))
	alpha = 1/vopp.vopp_ub_scale(S, 1)
	push!(x_grid, alpha)
	sort!(x_grid)

	for x in x_grid
		writedlm(f, [x vopp.tight_dist_ub_scale(x, S)])
	end
	close(f)
end

function scale_bound_vary_s(file_out_name; Smin=1, Smax=10, N=100)
#Create file_out and header
	outPath = "$(file_out_name)__$(Smin)_$(Smax)"
	f = open("$(outPath).tab", "w")
	writedlm(f, ["x" "Bound"])

	S_grid = range(Smin, stop=Smax, length=N)
	for S in S_grid
		writedlm(f, [S vopp.vopp_ub_scale(S, 1.)])
	end
	close(f)
end

function scale_bound_vary_M(file_out_name; Mmin=.1, Mmax=1, S=5, N=100)
#Create file_out and header
	outPath = "$(file_out_name)__$(Mmin)_$(Mmax)"
	f = open("$(outPath).tab", "w")
	writedlm(f, ["M" "Bound"])

	M_grid = range(Mmin, stop=Mmax, length=N)
	for M in M_grid
		Sc = vopp.comp_Sc(S, M)
		writedlm(f, [M vopp.vopp_ub_scale(Sc, 1.)])
	end
	close(f)
end

function scale_bound_vary_ratio(file_out_name; xmin=.1, xmax=1, N=100)
	#Create file_out and header
	outPath = "$(file_out_name)__$(xmin)_$(xmax)"
	f = open("$(outPath).tab", "w")
	writedlm(f, ["x" "Bound"])

	x_grid = range(xmin, stop=xmax, length=N)
	for x in x_grid
		writedlm(f, [x vopp.vopp_ub_scale(1/x, 1.)])
	end
	close(f)
end

function plot_tight_distns(file_out_name; S=4, N=100)
	#Create file_out and header
	outPath = "$(file_out_name)__$(S)_$(N)"
	f = open("$(outPath).tab", "w")
	writedlm(f, ["x" "Low" "Med" "High"])

	dl = vopp.delta_l(S)/2
	dm = (vopp.delta_m(S) + vopp.delta_l(S))/2
	dh = (vopp.delta_m(S) + vopp.delta_h(S))/2

	println("Devs $dl \t $dm \t $dh")
	#add the weird poitns to make plots better
	x_grid = collect(range(0, stop=S, length=N))
	push!(x_grid, 1.)
	push!(x_grid, 1e-10)
	push!(x_grid, 1/vopp.vopp_ub_MAD(S, 1., dl))
	push!(x_grid, exp(1) * S^(1- 1/dm))
	push!(x_grid, 1/vopp.vopp_ub_MAD(S, 1., dh) / (1 - dh))
	sort!(x_grid)

	for x in x_grid
		Fl = vopp.tight_dist_ub_MAD(x, S, dl)
		Fm = vopp.tight_dist_ub_MAD(x, S, dm)
		Fh = vopp.tight_dist_ub_MAD(x, S, dh)
		writedlm(f, [x Fl Fm Fh])
	end
	close(f)
end

if ARGS[1] == "Simple"
	simplePlot(ARGS[2])
elseif ARGS[1] == "Contour"
	contour_plot(ARGS[2])
elseif ARGS[1] == "scale_tight_dist"
	scale_bound_dist(ARGS[2])
elseif ARGS[1] == "scale_vary_s"
	scale_bound_vary_s(ARGS[2])
elseif ARGS[1] == "scale_vary_M"
	scale_bound_vary_M(ARGS[2])
elseif ARGS[1] == "scale_vary_ratio"
	scale_bound_vary_ratio(ARGS[2])
elseif ARGS[1] == "MAD_tight_dist"
	plot_tight_distns(ARGS[2])
end


text:
module Command

using JLD2, FileIO
using LinearAlgebra
using Statistics, StatsBase
using Plots

import BSON

include("../src/scrna.jl")
include("../src/pointcloud.jl")
include("../src/distance.jl")
include("../src/SeqSpace.jl")

using .SeqSpace, .PointCloud

# ------------------------------------------------------------------------
# variable inputs

mutable struct Parameters
    hyper  :: HyperParams
    dim    :: Int
    plot   :: Bool
    metric :: Function
end

const default = Parameters(
    HyperParams(),
    50,
    true,
    Distances.euclidean
)

Parameters(;
    hyper=default.hyper,
    dim=default.dim,
    plot=default.plot,
    metric=default.metric,
) = Parameters(hyper,dim,plot,metric)

include("utils.jl")
using .CommandUtility

# ------------------------------------------------------------------------
# plotting code

module Plot

using Plots

import ..PointCloud, ..Distances
using Statistics, StatsBase

function scaling(D, N; δ=2)
    ϕ, Rs = PointCloud.scaling(D, N)
    p = plot(Rs, ϕ[1:δ:end,:]',
           alpha=0.05,
           color=:red,
           xscale=:log10,
           yscale=:log10,
           label="",
    )

    plot!(Rs, size(ϕ,1)*(Rs ./ maximum(Rs)).^3,
           color=:black,
           linestyle=:dashdot,
           linewidth=2,
           label="3d scaling",
    )
    xaxis!("ball size")
    yaxis!("number of points enclosed", (20,size(D,1)))

    return p
end

function isomapfit(ξ,D)
    c = [ cor(PointCloud.upper_tri(Distances.euclidean(ξ[:,1:i]')),
              PointCloud.upper_tri(D)) for i ∈ 1:10 ]
    p = plot(c, label="", linewidth=2, addmarker=true)

    xaxis!("number of isomap dimensions")
    yaxis!("pairwise distance correlation")

    return p
end

function fitequilibrate(loss)
    p = plot(loss.train,label="training",linewidth=3,yscale=:log10)
    plot!(loss.valid,label="validatation",linewidth=3,yscale=:log10)
    xlabel!("epoch")
    ylabel!("loss")

    return p
end

function fitgoodness(x, y)
    c = cor(x[:], y[:])
    p = histogram2d(x,y,title="corr=$(c)")
    xlabel!("data")
    ylabel!("predicted")

    return p
end

end

# ------------------------------------------------------------------------
# bulk functionality

function mainconncomp(distances)
    adj = distances .!= Inf
    ncs = sum(adj,dims=1)
    ccs = Set(ncs)

    return vec(ncs .== maximum(ccs))
end

function autoencode(data, param::Parameters, figs::AbstractString)
    alert("---->projecting data sized $(size(data)) to $(param.dim) dimensional linear subspace")
    input = linearprojection(data, param.dim)

    alert("---->computing geodesics")
    D₀ = param.metric(input.projection)
    χ  = percentile(upper_tri(D₀), 5)
    D  = geodesics(input.projection, param.hyper.k; D=D₀, accept=(d)->d≤χ)

    ι = mainconncomp(D)
    D = D[ι,ι]
    input = (projection=input.projection[:,ι], embed=input.embed)
    alert("---->main connected component sized $(sum(ι))")

    param.plot && let
        savefig(Plot.scaling(D₀, 100), "$figs/pointcloud_scaling_base.png")
        savefig(Plot.scaling(D,  100), "$figs/pointcloud_scaling_geodesic.png")
    end

    alert("---->computing isomap coordinates")
    ξ = real.(mds(D.^2, 10))
    param.plot && savefig(Plot.isomapfit(ξ,D), "$figs/isomap_corr.png")

    alert("---->fitting autoencoder replicates")
    # TODO: add niter parameter
    result, metadata = fitmodel(input.projection, param.hyper; D²=D.^2, chatty=false)
    param.plot && let
        savefig(Plot.fitequilibrate(result.loss), "$figs/loss_dynamics.png")
        savefig(Plot.fitgoodness(input.projection, result.model.identity(input.projection)), "$figs/fit_residual.png")
    end

    return (
        input    = input,
        distance = D,
        isomap   = ξ,
        fit      = marshal(result),
        latent   = result.model.pullback(input.projection),
    )
end

# ------------------------------------------------------------------------
# main point of entry

if abspath(PROGRAM_FILE) == @__FILE__
    params = (
        p = Arg("inpath", loadparams, default),
        o = Arg("outpath", (path)->path, nothing)
    )

    args = argparse(ARGS,params)
    length(args) > 1 && error("too many input paths")

    input  = args[1]

    outdir = dirname(params.o.value)
    !isdir(outdir) && mkdir(outdir)

    figdir = "$(outdir)/figs/model"
    !isdir(figdir) && mkpath(figdir)

    ENV["GKSwstype"] = "nul"
    output = params.o.value
    jldopen(output, "w") do model; jldopen(input, "r") do norm
        for name in keys(norm)
            alert("-->processing group $(name)")
            result = autoencode(norm[name]["data"], params.p.value, figdir)

            model["$name/fitparams"] = params.p.value
            for (key,val) in pairs(result)
                model["$name/$key"] = val
            end
        end
    end; end
end

end

max_stars_repo_name: nnoll/seqspace

max_issues_repo_name: nnoll/seqspace

max_forks_repo_name: nnoll/seqspace


text:
@everywhere const roots3 = Complex64[1.0, exp(2im*pi/3), exp(4im*pi/3)]
@everywhere const colors3 = ([1.0,0.0,0.0],[0.0,1.0,0.0],[0.0,0.0,1.0])

@everywhere function cubic(z0, roots, colors, maxiter, eps)
	n = 1
	z = z0 - (z0*z0*z0 - 1)/(3.0*z0*z0)
	while n <= maxiter
		for j in 1:length(roots)
			if abs2(z - roots[j]) <= eps*eps
				return (1.0 - (n-1)/(maxiter-1))*colors[j]
			end
		end
		z -= (z*z*z - 1)/(3.0*z*z)
		n += 1
	end
	return [0.0,0.0,0.0]
end

max_stars_repo_name: kalvotom/julia-examples

max_issues_repo_name: kalvotom/julia-examples

max_forks_repo_name: kalvotom/julia-examples


text:
using SegyIO
using DelimitedFiles

function line2segy(infile)
    print(infile);
    inversion_output = readdlm(infile);
    line_number = unique(round.(Int, inversion_output[:,line_col]))[1];

    # our first step is to define our datum as 20 m above the highest point in survey

    elevations = inversion_output[:,elevation_col]
    ymax = ceil(maximum(elevations)/10.) *10.

    # create an elevation grid

    ymin = floor(minimum(elevations)/10.)*10. - max_depth

    grid_elevations = collect(range(ymax, ymin, length = Int(ceil((ymax - ymin)/yres))))

    # get conductivity
    sigma = transpose(inversion_output[:,conductivity_cols[1]:conductivity_cols[2]])
    thickness = transpose(inversion_output[:,thickness_cols[1]:thickness_cols[2]])
    #layer_top_elevation = inversion_output[:,layer_top_elevation_cols[1]:layer_top_elevation_cols[2]]

    ## get coordinates

    easting = inversion_output[:,easting_col]
    northing = inversion_output[:,northing_col]
    ## create our array of interpolated AEM conductivity, negative indicates free air
    interpolated = -1*ones(Float32, length(grid_elevations), length(easting))

    # iterate
    for i in 1:length(easting)

        layer_top_elev = elevations[i] * ones(size(thickness)[1])
        layer_top_elev[2:end] = layer_top_elev[2:end] - cumsum(thickness[1:end-1,i])
        #layer_top_elev = layer_top_elevation[i,:]
        sigma_trace = sigma[:,i]
        # iterate through our layers
        for j in 1:(length(layer_top_elev) - 1)
            interpolated[(grid_elevations .< layer_top_elev[j]) .& (grid_elevations .>= layer_top_elev[j + 1]),i] .= sigma_trace[j]
        end
    end

    ## try writing these out in depth space

    block = SeisBlock(interpolated)

    set_header!(block, "SourceX", round.(Int, easting))
    set_header!(block, "SourceY", round.(Int, northing))
    set_header!(block, "GroupX", round.(Int, easting))
    set_header!(block, "GroupY", round.(Int, northing))
    set_header!(block, "ElevationScalar", round.(Int, elevations))
    set_header!(block, "TraceNumWithinLine", Array(1:length(easting)))
    set_header!(block, "TraceNumWithinFile", Array(1:length(easting)))
    #set_header!(block, "TraceIDCode", round.(Int, fiducial * 10))
    set_header!(block, "RecGroupElevation", round.(Int, elevations))
    set_header!(block, "DelayRecordingTime", Int(ymax) * -1)
    set_header!(block, "dt", Int(yres * 1e3))

    ##
    segy_write(joinpath(outdir,string(line_number) * ".segy"), block)
end

##
# choose the vertical resampling interval
yres = 4.0;
max_depth = 500.;

##
indir = raw"C:\Users\u77932\Documents\LEB\data\AEM\Frome\galei\lines2";
outdir = raw"C:\Users\u77932\Documents\LEB\data\AEM\Frome\galei\segy";

cd(indir);
files = readdir();

# define column numbers
easting_col = 7;
northing_col = 8;
elevation_col = 9;
conductivity_cols = [24,53] #ranger for 30 layer model;
thickness_cols = [54,83];
#layer_top_elevation_cols = [14 43]
line_col = 5;

##
for fname in files

    line2segy(fname)
end


text:
#=
multiply:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-08-06
=#

module Multiplication
    function multiply(a, b)
        return a * b
    end

    function main()
        a = 15
        b = 12
        println("a = $a, b = $b")
        println()
        println("use *")
        println("a * b = $(multiply(a, b))")
        println()
        println("omit *")
        println("2a + b = $(2a + b)")
        println("(a + 2)b = ", (a + 2)b)
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    using .Multiplication
    Multiplication.main()
end

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
#
# export JULIA_DEPOT_PATH=`pwd`/.julia
#
using Pkg
Pkg.instantiate()

using ClimateMachine

max_stars_repo_name: christophernhill/clima-impldiff-testing

max_issues_repo_name: christophernhill/clima-impldiff-testing

max_forks_repo_name: christophernhill/clima-impldiff-testing


text:
#= ForwardDiff is removed from dependencies in this version.
import MechanicalSketch
import MechanicalSketch: color_with_lumin, O, WI, HE, EM, FS, finish
import MechanicalSketch: PALETTE, setfont, settext, background, empty_figure, sethue, move, do_action
import MechanicalSketch: dimension_aligned, noise, line, lenient_max
import MechanicalSketch: ComplexQuantity, generate_complex_potential_source, generate_complex_potential_vortex
import MechanicalSketch: Quantity, x_y_iterators_at_pixels, @import_expand
import MechanicalSketch: place_image, set_scale_sketch, draw_legend
import MechanicalSketch: ∙, ∇_rectangle, SVector, normalize_datarange, BinLegend
import MechanicalSketch: Greys_9
import Interpolations: interpolate, Gridded, Linear, Flat, extrapolate

#using BenchmarkTools

let
BACKCOLOR = color_with_lumin(PALETTE[7], 10)

include("test_functions_25.jl")
restart_25(BACKCOLOR)
include("test_functions_24.jl")

# Reuse the flow field from test_23.jl, a matrix of complex velocities with an element per pixel.
A = flowfield_23();
set_scale_sketch(PHYSWIDTH_23, round(Int, SCREEN_WIDTH_FRAC_23 * WI))

#setfont("DejaVu Sans", FS)
str = "This noise image has three different spectra.\nA noise frequency along a line is 1 / (wavelength [m])."
settext(str, O + (-WI/2 + 1.25EM, -0.5HE + 2EM), markup = true)

# We are going to use the velocity vector field in a lot of calculations,
# and interpolate between the calculated pixel values
xs, ys = x_y_iterators_at_pixels(;physwidth = PHYSWIDTH_23, physheight = PHYSHEIGHT_23)
fxy_inter = interpolate((xs, reverse(ys)), map( cmplx -> (real(cmplx), imag(cmplx)), transpose(A)[ : , end:-1:1]), Gridded(Linear()));
fxy = extrapolate(fxy_inter, Flat());

# Generate a noise image for experimentation. Display it at the top of the figure.
NO_1 = noisepic_1(A);
toplegend = BinLegend(;maxlegend = lenient_max(NO_1), noofbins = 16,
                       colorscheme = reverse(Greys_9), name = :Noise)
colmat = toplegend.(NO_1)
upleftpoint, lowrightpoint = place_image(O + (0.0, -0.25HE + EM), colmat)
legendpos = lowrightpoint + (EM, 0) + (0.0m, PHYSHEIGHT_23)
draw_legend(legendpos, toplegend)

# And make a function that interpolates between pixels:
nxy_inter = interpolate((xs, reverse(ys)), transpose(NO_1)[ : , end:-1:1], Gridded(Linear()));
nxy = extrapolate(nxy_inter, Flat());
nxy(0.1m, 0.02m)

# Centre of the bottom figure
OB_25 = O + (0.0, + 0.25HE + 0.75EM )

# Duration to trace the streamlines forward and back
DU_25 = 2.0s

# Number of sampling points along our streamline in each direction (the number of steps is one less )
global const NS_1 = 10

# Step length (duration)
h = DU_25 / (NS_1 -1 )

@time M = convolute_image_1(xs, reverse(ys), fxy, nxy, h, CUTOFF_23) # 29.268005 seconds (1.28 M allocations: 71.975 MiB, 1.78% compilation time)

botlegend = BinLegend(;maxlegend = lenient_max(M), noofbins = 256,
                       colorscheme = reverse(Greys_9), name = :LIC)
colmat = botlegend.(M)
upleftpoint, lowrightpoint = place_image(OB_25, colmat)
legendpos = lowrightpoint + (EM, 0) + (0.0m, PHYSHEIGHT_23)
draw_legend(legendpos, botlegend)

draw_streamlines_1(OB_25, xs, ys, fxy, h)

dimension_aligned(OB_25 + (-PHYSWIDTH_23  / 2, PHYSHEIGHT_23  / 2), OB_25 + (PHYSWIDTH_23  / 2, PHYSHEIGHT_23  / 2))
dimension_aligned(OB_25, OB_25 + complex(0.0, 1.0)m)
dimension_aligned(OB_25 + (-PHYSWIDTH_23  / 2, - PHYSHEIGHT_23  / 2 ),  OB_25 +  (-PHYSWIDTH_23  / 2, PHYSHEIGHT_23  / 2 ))

str = "The grey image below is a first attempt at LIC  flow visualization.\nStreamlines cover $DU_25 forward (blue) and back."
settext(str, O + (-WI / 2 + EM, + 1.8EM ), markup = true)

finish()
end # begin
=#


text:
# A quadratic triangle, defined in 2D.
struct Triangle6_2D <: Face_2D
    # The points are assumed to be ordered as follows
    # p₁ = vertex A
    # p₂ = vertex B
    # p₃ = vertex C
    # p₄ = point on the quadratic segment from A to B
    # p₅ = point on the quadratic segment from B to C
    # p₆ = point on the quadratic segment from C to A
    points::SVector{6, Point_2D}
end

# Constructors
# -------------------------------------------------------------------------------------------------
Triangle6_2D(p₁::Point_2D, p₂::Point_2D, p₃::Point_2D,
             p₄::Point_2D, p₅::Point_2D, p₆::Point_2D
            ) = Triangle6_2D(SVector(p₁, p₂, p₃, p₄, p₅, p₆))

# Methods
# -------------------------------------------------------------------------------------------------
# Interpolation
function (tri6::Triangle6_2D)(r::Real, s::Real)
    # See Fhe Visualization Toolkit: An Object-Oriented Approach to 3D Graphics, 4th Edition
    # Chapter 8, Advanced Data Representation, in the interpolation functions section
    rₜ = Float64(r); sₜ = Float64(s)
    return (1 - rₜ - sₜ)*(2(1 - rₜ - sₜ) - 1)*tri6[1] +
                                   rₜ*(2rₜ-1)*tri6[2] +
                                   sₜ*(2sₜ-1)*tri6[3] +
                            4rₜ*(1 - rₜ - sₜ)*tri6[4] +
                                     (4rₜ*sₜ)*tri6[5] +
                            4sₜ*(1 - rₜ - sₜ)*tri6[6]
end

# Interpolation using a point instead of (r,s)
function (tri6::Triangle6_2D)(p::Point_2D)
    r,s = p
    return (1 - r - s)*(2(1 - r - s) - 1)*tri6[1] +
                                 r*(2r-1)*tri6[2] +
                                 s*(2s-1)*tri6[3] +
                           4r*(1 - r - s)*tri6[4] +
                                   (4r*s)*tri6[5] +
                           4s*(1 - r - s)*tri6[6]
end

area(tri6::Triangle6_2D) = area(tri6, Val(3))
function area(tri6::Triangle6_2D, ::Val{N}) where {N}
    # Numerical integration required. Gauss-Legendre quadrature over a triangle is used.
    # Let F(r,s) be the interpolation function for tri6,
    #                             1 1-r                          N
    # A = ∬ ||∂F/∂r × ∂F/∂s||dA = ∫  ∫ ||∂F/∂r × ∂F/∂s|| ds dr = ∑ wᵢ||∂F/∂r(rᵢ,sᵢ) × ∂F/∂s(rᵢ,sᵢ)||
    #      D                      0  0                          i=1
    #
    # N is the number of points used in the quadrature.
    w, r, s = gauss_legendre_quadrature(tri6, Val(N))
    return sum(@. w * abs( ∇(tri6, r, s) |> x->x[1] × x[2] ))
end

centroid(tri6::Triangle6_2D) = centroid(tri6, Val(6))
function centroid(tri6::Triangle6_2D, ::Val{N}) where {N}
    # Numerical integration required. Gauss-Legendre quadrature over a triangle is used.
    # Let F(r,s) be the interpolation function for tri6,
    #                             1 1-r                          N
    # A = ∬ ||∂F/∂r × ∂F/∂s||dA = ∫  ∫ ||∂F/∂r × ∂F/∂s|| ds dr = ∑ wᵢ||∂F/∂r(rᵢ,sᵢ) × ∂F/∂s(rᵢ,sᵢ)||
    #      D                      0  0                          i=1
    #
    # C_x = (∫∫ x dA)/A, C_y = (∫∫ y dA)/A
    #         D                  D
    w, r, s = gauss_legendre_quadrature(tri6, Val(N))
    # We can reuse our computed weighted derivative cross products, since we need these
    # in the C_y, C_y, and A.
    weighted_vals = @. w * abs( ∇(tri6, r, s) |> x->x[1] × x[2] )
    points = tri6.(r, s)
    A = sum(weighted_vals)
    C_x = sum(getindex.(points, 1) .* weighted_vals)
    C_y = sum(getindex.(points, 2) .* weighted_vals)
    return Point_2D(C_x/A, C_y/A)
end

function gradient(tri6::Triangle6_2D, r::Real, s::Real)
    # Let F(r,s) be the interpolation function for tri6
    # Returns ∂F/∂r, ∂F/∂s
    rₜ = Float64(r); sₜ = Float64(s)
    ∂F_∂r = (4rₜ + 4sₜ - 3)*tri6[1] +
                  (4rₜ - 1)*tri6[2] +
            4(1 - 2rₜ - sₜ)*tri6[4] +
                      (4sₜ)*tri6[5] +
                     (-4sₜ)*tri6[6]

    ∂F_∂s = (4rₜ + 4sₜ - 3)*tri6[1] +
                  (4sₜ - 1)*tri6[3] +
                     (-4rₜ)*tri6[4] +
                      (4rₜ)*tri6[5] +
            4(1 - rₜ - 2sₜ)*tri6[6]
    return ∂F_∂r, ∂F_∂s
end

function in(p::Point_2D, tri6::Triangle6_2D)
    # If the point is to the left of every edge
    #  3<-----2
    #  |     ^
    #  | p  /
    #  |   /
    #  |  /
    #  v /
    #  1
    return isleft(p, QuadraticSegment_2D(tri6[1], tri6[2], tri6[4])) &&
           isleft(p, QuadraticSegment_2D(tri6[2], tri6[3], tri6[5])) &&
           isleft(p, QuadraticSegment_2D(tri6[3], tri6[1], tri6[6]))
end

# # Slower method than above.
# function in(p::Point_2D, tri6::Triangle6_2D)
#     return in(p, tri6, 30)
# end
#
# function in(p::Point_2D, tri6::Triangle6_2D, N::Int64)
#     # Determine if the point is in the triangle using the Newton-Raphson method
#     # N is the max number of iterations of the method.
#     p_rs = real_to_parametric(p, tri6, N)
#     ϵ = parametric_coordinate_ϵ
#     # Check that the r coordinate and s coordinate are in [-ϵ,  1 + ϵ] and
#     # r + s ≤ 1 + ϵ
#     # These are the conditions for a valid point in the triangle ± some ϵ
#     if (-ϵ ≤ p_rs[1] ≤ 1 + ϵ) &&
#        (-ϵ ≤ p_rs[2] ≤ 1 + ϵ) &&
#        (p_rs[1] + p_rs[2] ≤ 1 + ϵ)
#         return true
#     else
#         return false
#     end
# end

function intersect(l::LineSegment_2D, tri6::Triangle6_2D)
    # Create the 3 quadratic segments that make up the triangle and intersect each one
    edges = SVector(QuadraticSegment_2D(tri6[1], tri6[2], tri6[4]),
                    QuadraticSegment_2D(tri6[2], tri6[3], tri6[5]),
                    QuadraticSegment_2D(tri6[3], tri6[1], tri6[6]))
    ipoints = MVector(Point_2D(), Point_2D(), Point_2D(),
                      Point_2D(), Point_2D(), Point_2D())
    n_ipoints = 0x00000000
    # We need to account for 6 points returned
    for k ∈ 1:3
        npoints, points = l ∩ edges[k]
        for i ∈ 1:npoints
            n_ipoints += 0x00000001
            ipoints[n_ipoints] = points[i]
        end
    end
    return n_ipoints, SVector(ipoints)
end

function jacobian(tri6::Triangle6_2D, r::Real, s::Real)
    # Return the 2 x 2 Jacobian matrix
    ∂F_∂r, ∂F_∂s = ∇(tri6, r, s)
    return SMatrix{2, 2}(∂F_∂r.x, ∂F_∂r.y,
                         ∂F_∂s.x, ∂F_∂s.y)
end

function real_to_parametric(p::Point_2D, tri6::Triangle6_2D)
    return real_to_parametric(p, tri6, 30)
end

function real_to_parametric(p::Point_2D, tri6::Triangle6_2D, N::Int64)
    # Convert from real coordinates to the triangle's local parametric coordinates using the
    # the Newton-Raphson method. N is the max number of iterations
    # If a conversion doesn't exist, the minimizer is returned.
    r = 0.3333333333333333 # Initial guess at triangle centroid
    s = 0.3333333333333333
    for i ∈ 1:N
        err = p - tri6(r, s)
        # Inversion is faster for 2 by 2 than \
        Δr, Δs = inv(jacobian(tri6, r, s)) * err
        r += Δr
        s += Δs
        if abs(Δr) + abs(Δs) < 1e-6
            break
        end
    end
    return Point_2D(r, s)
end

function triangulate(tri6::Triangle6_2D, N::Int64)
    # N is the number of divisions of each edge
    triangles = Vector{Triangle_2D}(undef, (N+1)*(N+1))
    if N === 0
        triangles[1] = Triangle_2D(tri6[1], tri6[2], tri6[3])
    else
        i = 1
        for S = 1:N, R = 0:N-S
            triangles[i]   = Triangle_2D(tri6(    R/(N+1),     S/(N+1)),
                                         tri6((R+1)/(N+1),     S/(N+1)),
                                         tri6(    R/(N+1), (S+1)/(N+1)))
            triangles[i+1] = Triangle_2D(tri6(    R/(N+1),     S/(N+1)),
                                         tri6((R+1)/(N+1), (S-1)/(N+1)),
                                         tri6((R+1)/(N+1),     S/(N+1)))
            i += 2
        end
        j = (N+1)*N + 1
        for S = 0:0, R = 0:N-S
            triangles[j] = Triangle_2D(tri6(    R/(N+1),     S/(N+1)),
                                       tri6((R+1)/(N+1),     S/(N+1)),
                                       tri6(    R/(N+1), (S+1)/(N+1)))
            j += 1
        end
    end
    return triangles
end

# Plot
# -------------------------------------------------------------------------------------------------
if enable_visualization
    function convert_arguments(LS::Type{<:LineSegments}, tri6::Triangle6_2D)
        q₁ = QuadraticSegment_2D(tri6[1], tri6[2], tri6[4])
        q₂ = QuadraticSegment_2D(tri6[2], tri6[3], tri6[5])
        q₃ = QuadraticSegment_2D(tri6[3], tri6[1], tri6[6])
        qsegs = [q₁, q₂, q₃]
        return convert_arguments(LS, qsegs)
    end

    function convert_arguments(LS::Type{<:LineSegments}, T::Vector{Triangle6_2D})
        point_sets = [convert_arguments(LS, tri6) for tri6 ∈ T]
        return convert_arguments(LS, reduce(vcat, [pset[1] for pset ∈ point_sets]))
    end

    function convert_arguments(P::Type{<:Mesh}, tri6::Triangle6_2D)
        triangles = triangulate(tri6, 13)
        return convert_arguments(P, triangles)
    end

    function convert_arguments(M::Type{<:Mesh}, T::Vector{Triangle6_2D})
        triangles = reduce(vcat, triangulate.(T, 13))
        return convert_arguments(M, triangles)
    end
end

max_stars_repo_name: KyleVaughn/MOCNeutronTransport

max_issues_repo_name: KyleVaughn/MOCNeutronTransport

max_forks_repo_name: KyleVaughn/MOCNeutronTransport


text:
module RandSeis
using Random, SeisIO, SeisIO.Quake

include("RandSeis/imports.jl")
if VERSION <= v"1.1.0"
  include("RandSeis/constants_1.jl")
else
  include("RandSeis/constants.jl")
end
include("RandSeis/iccodes_and_units.jl")
include("RandSeis/utils.jl")
include("RandSeis/randSeisChannel.jl")
include("RandSeis/randSeisData.jl")
include("RandSeis/randSeisEvent.jl")
include("RandSeis/exports.jl")

end


text:
function _function(forward::Function, input::Variable)::Variable
    x = input.data
    y::Array = forward(x)
    output = Variable(y)
    return output
end

abstract type _Function end
struct Square <: _Function end
(::Square)(input::Variable)::Variable = _function(x -> x .^ 2, input)


text:
using Test

include("data.jl")

@testset "randmgn tests" begin

    for i=1:100
        input_size = 15
        perm = randperm(input_size)
        a_inds = perm[1:div(input_size, 2)]
        b_inds = perm[div(input_size, 2)+1:end]

        min = 0
        max = rand(1:100)*rand()
        ops = [+, -, *, \]
        op = rand(ops)
        val = randmgn(op, input_size, min, max, a_inds, b_inds)
        @test abs(sum(val[a_inds])) > min
        @test abs(sum(val[a_inds])) < max
        @test abs(sum(val[b_inds])) > min
        @test abs(sum(val[b_inds])) < max
        yabs = abs(op(sum(val[a_inds]), sum(val[b_inds])))
        @test yabs < max
        @test yabs > min
    end

end

@testset "fixed binary_data tests" begin
    input_size = 10
    perm = collect(1:10)

    for op in [+, -, *, /]
        (xs, ys), (_, a_inds, b_inds, abs) = binary_data(op, 100, input_size, 5, perm=perm)
        for i in 1:length(xs)
            @test sum(xs[i][1:5]) == abs[i][1]
            @test sum(xs[i][6:10]) == abs[i][2]
            @test op(sum(xs[i][1:5]), sum(xs[i][6:10])) == ys[i]
        end
    end
end

@testset "random binary_data tests" begin
    input_size = 20

    for op in [+, -, *, /]
        (xs, ys), (perm, a_inds, b_inds, abs) = binary_data(op, 100, input_size, 10)
        (xs2, ys2), (perm2, a_inds2, b_inds2, abs2) = binary_data(op, 100, input_size, 10, perm=perm)

        @test all(perm .== perm2)
        @test all(a_inds .== a_inds2)
        @test all(b_inds .== b_inds2)

        @testset "random perm check" begin
            for i in 1:length(xs)
                @test sum(xs[i][a_inds]) == abs[i][1]
                @test sum(xs[i][b_inds]) == abs[i][2]
                @test op(sum(xs[i][a_inds]), sum(xs[i][b_inds])) == ys[i]
            end
        end

        @testset "re-used perm check" begin
            for i in 1:length(xs2)
                @test sum(xs2[i][a_inds]) == abs2[i][1]
                @test sum(xs2[i][b_inds]) == abs2[i][2]
                @test op(sum(xs2[i][a_inds]), sum(xs2[i][b_inds])) == ys2[i]
            end

        end
    end

end

@testset "unary_data tests" begin
    input_size = 10
    perm = collect(1:10)

    sqr(x) = ^(x, 2)

    for op in [sqr, sqrt]
        (xs, ys), (_, a_inds, b_inds, abs) = unary_data(op, 100, input_size, 5, perm)
        for i in 1:length(xs)
            @test sum(xs[i][1:5]) == abs[i][1]
            @test sum(xs[i][6:10]) == abs[i][2]
            @test op(sum(xs[i][1:5])) == ys[i]
        end
    end
end

max_stars_repo_name: dillondaudert/neural-alu

max_issues_repo_name: dillondaudert/neural-alu

max_forks_repo_name: dillondaudert/neural-alu


text:
"""
    stalta(s::SACtr, t, ws, wl) -> ::SACFloat

Return the ratio of the short-term average (window length `ws` s) to long-term
average (length `wl` s) at time `t` s.

    stalta(s::SACtr, t1, t2, ws, wl) -> ::Vector{SACFloat}

Return STA/LTA for each sample between times `t1` and `t2` s.

    stalta(s::SACtr, ws, wl) -> ::SACtr

Return STA/LTA for whole trace as a `SACtr` object.

#### Note on window positions

In this implemetation, the ratio refers to that between a long time window *before*
time `t` to a short window *after* `t`, and here the short window includes the sample
*at* `t` (or the nearest thereto).
"""
stalta(s::SACtr, t, ws, wl) = window_average(s, t, ws)/window_average(s, t-wl-s.delta, wl)

function stalta(s::SACtr, t1, t2, ws, wl)
   n = round(Int, (t2 - t1)/s.delta) + 1
   ratio = Array{eltype(s.t)}(undef, n)
   for i in 1:n
       t = t1 + (i-1)*s.delta
       ratio[i] = stalta(s, t, ws, wl)
   end
   ratio
end

function stalta(s::SACtr, ws, wl)
    t1 = s.b + wl + 2s.delta
    t2 = s.e - ws - s.delta
    t = stalta(s, t1, t2, ws, wl)
    s2 = deepcopy(s)
    s2.t = t
    s2.npts = length(t)
    s2.b = t1
    update_headers!(s2)
end

"""
    window_average(s::SACtr, ref_time, window_length, centred=false)

Return the mean value of the absolute amplitude within a window
`window_length` s long.  `ref_time` in s may be the start of the window (default)
or the centre if `centred` is `true`.

Note that window times down to 0 are allowed, in which case the amplitude of
the nearest point is returned.
"""
function window_average(s::SACtr, time, window_length, centred=false)
    window_length >= 0 || error("`window_length` ($window_length) must be positive")
    n = round(Int, window_length/s.delta) + 1
    it1 = round(Int, (time - s.b - (centred ? window_length/2 : 0.0))/s.delta) + 1
    it2 = it1 + n - 1
    it1 > 0 && it2 <= s.npts ||
        error("Window of times $time - $(time+window_length) s lies outside " *
              "trace (trace indices $it1 - $it2; start time $(s.b) s; end time $(s.e)) s")
    ave = zero(SACFloat)
    @inbounds for i in it1:it2
        ave += abs(s.t[i])
    end
    ave/n
end


text:
## Simplex Benchmark
## =================
# This benchmark samples from a small region around corners of a n-dimensional
# Simplex. If the regions are of equal size then we should expect an equal
# number of samples in each region.

"""The vertex coordinates for an ndim simplex
From http://people.sc.fsu.edu/~jburkardt/m_src/simplex_coordinates/"""
function simplex_coordinates(n::Int)
  x = zeros(Float64, n, n+1)
  for i = 1 : n
    x[i,i] = sqrt(1.0 - sum(x[1:i-1,i].^2))
    for j = i+1: n+1
      b = (-1.0/n - (x[1:i-1,i]' * x[1:i-1,j])) / x[i,i]
      x[i,j] = b[1]
    end
  end
  return x
end

"Is point p in box with center region and halfwidth halfwidth"
function point_near_region(p, region, halfwidth)
  conds = Array(Any,length(p))
  for i = 1:length(p)
    conds[i] = (p[i] >= region[i] - halfwidth) & (p[i] <= region[i] + halfwidth)
  end
  (&)(conds...)
end

"Constrain a point to be within a region around the vertices of a simplex"
function simplex(n::Int, box::RandArray, holesize)
  vertexcoords = simplex_coordinates(n)
  conds = RandVar{Bool}[]
  for i = 1:n+1
    push!(conds,point_near_region(box,vertexcoords[:,i],holesize))
  end
  box,(|)(conds...)
end

"Ground truth histrogram"
groundtruth(ndims) = Dict(i::Int => (1/(ndims+1))::Float64 for i = 1:(ndims+1))

"Constrain a point to be within a region around the vertices of a simplex"
simplex(n::Int) = simplex(n,mvuniform(-2,2,n))

"From a set of samples compute discrete distribution over the vertices"
function vertex_distribution(samples,n,holesize)
  vertexcoords = simplex_coordinates(n)
  counts =  DefaultDict(Int, Float64, 0.0)
  for sample in samples
    for i = 1:n+1
      if point_near_region(sample,vertexcoords[:,i], holesize)
        counts[i] += 1
        break
      end
    end
  end
  Dict(counts)
end

## Benchmark
## =========
"""Sample from a prior over R^ndim constrained to cubes of size holesize around
vertices of ndim Simplex"""
immutable Simplex <: Problem
  ndims::Int
  capture::Vector{Symbol}
  nsamples::Int
  holesize::Float64
  prior::RandArray
end

function benchmark(a::Algorithm, b::Simplex)
  print("Starting Sampler")
  Sigma.restart_counter!()
  captures::Vector{Symbol} = vcat(a.capture,b.capture)
  groundtruth = Dict(i => 1/(b.ndims+1) for i = 1:(b.ndims+1))
  model, condition = simplex(b.ndims, b.prior, b.holesize)
  @show b.nsamples

  value, results = capture(()->sample(a,model,condition,b.nsamples), captures)
  @show length(value)
  results
end

## Algorithm-dependent implementations of sample
## ============================================
"Sample with SMT based solver"
function sample(a::SigmaSMT, model, condition, nsamples)
  @show model
  @show condition
  @show nsamples
  rand(model, condition, nsamples; preimage_sampler = a.preimage_sampler,
                                   RandVarType = a.randvartype,
                                   cores = a.ncores,
                                   split = a.split)
end

"Sample with rejection sampling"
function sample(a::SigmaRej, X::AllRandVars, Y::RandVar{Bool}, nsamples::Integer)
  box = Sigma.LazyOmega(Float64)
  for dim in dims(Y)
    box[dim]
  end

  Y_func = Sigma.lambda(Y)
  X_func = Sigma.lambda(X)

  samples = []
  ntries = 10000000
  for i = 1:nsamples
    for i = 1:ntries
      ω = rand(box)
      if Y_func(ω)
        println("Got sample with rejection")
        push!(samples,X_func(ω))
        break
      end
    end
  end

  @show length(samples), nsamples
  length(samples) != nsamples && error("Couldn't draw sample")
  return samples
end

"Sample using breadth first search"
function sample(a::SigmaBFS, X::AllRandVars, Y::RandVar{Bool}, nsamples::Integer)
  Sigma.rand_bfs(X,Y,nsamples)
end

## Runs
## ====
include("runs/kl.jl")
# include("runs/dimensions.jl")


text:
## For now, just make sure the code runs ##

@testset "Lightcurve" begin
    n = 8
    BJD = 7250.0; t0_ic = 7.0
    tmax = 2.0
    ic = setup_ICs(n, BJD, t0_ic);
    intr = setup_integrator(ic, tmax);
    tt = compute_transit_times(ic, intr);
    pd = compute_pd(ic, tt, intr);

    # Transit parameters
    k = get_radius_ratios_trappist(n);
    u_n = get_limdark_coeffs_trappist();
    rstar = 0.00465047 * 0.1192 # Trappist-1 (Rstar/AU)

    # Simulate lightcurve
    dt = 2 / 60 / 24; # 2 minutes in days
    obs_duration = maximum(tt.tt) + 5 * intr.h
    nobs = round(Int64, obs_duration * 30 * 24)
    tobs = collect(pd.times[1] - 10 * intr.h:obs_duration / nobs:obs_duration);
    lc = Lightcurve(dt, copy(tobs), copy(tobs), zeros(length(tobs)), u_n, k, rstar);
    lc_noint = Lightcurve(0.0, copy(tobs), copy(tobs), zeros(length(tobs)), u_n, k, rstar);
    dlc = Lightcurve(dt, copy(tobs), copy(tobs), zeros(length(tobs)), u_n, k, rstar, n*7);
    compute_lightcurve!(lc, pd, tol=1e-11, maxdepth=40)
    compute_lightcurve!(lc_noint, pd)

    @time compute_lightcurve!(lc, pd, tol=1e-11, maxdepth=40)
    @time compute_lightcurve!(lc_noint, pd)
    @time compute_lightcurve!(dlc, pd, tol=1e-11, maxdepth=40)
    @test 1 == 1
end


text:
"""
    xT,RT,ll = smooth(kf::KalmanFilter, u::Vector, y::Vector)
Returns smoothed estimates of state `x` and covariance `R` given all input output data `u,y`
"""
function smooth(kf::KalmanFilter, u::AbstractVector, y::AbstractVector)
    reset!(kf)
    T            = length(y)
    x,xt,R,Rt,ll = forward_trajectory(kf, u, y)
    xT           = similar(xt)
    RT           = similar(Rt)
    xT[end]      = xt[end]      |> copy
    RT[end]      = Rt[end]      |> copy
    for t = T-1:-1:1
        C     = Rt[t]*kf.A/R[t+1]
        xT[t] = xt[t] .+ C*(xT[t+1] .- x[t+1])
        RT[t] = Rt[t] .+ symmetrize(C*(RT[t+1] .- R[t+1])*C')
    end
    xT,RT,ll
end

function smooth(pf::AbstractParticleFilter, M, u, y)
    xf,wf,wef,ll = forward_trajectory(pf, u, y)
    smooth(pf::AbstractParticleFilter, xf, wf, wef, ll, M, u, y)
end

"""
    xb,ll = smooth(pf, M, u, y)

Perform particle smoothing using forward-filtering, backward simulation. Return smoothed particles and loglikelihood.
See also `smoothed_trajs`, `smoothed_mean`, `smoothed_cov`
"""
function smooth(pf::AbstractParticleFilter, xf, wf, wef, ll, M, u, y)
    T = length(y)
    N = num_particles(pf)
    f = dynamics(pf)
    df = dynamics_density(pf)
    @assert M <= N "Must extend cache size of bins and j to allow this"
    xb = Array{particletype(pf)}(undef,M,T)
    j = resample(ResampleSystematic, wef[:,T], M)
    # @show Set(j)
    for i = 1:M
        xb[i,T] = xf[j[i], T]
    end
    wb = Vector{Float64}(undef,N)
    @inbounds for t = T-1:-1:1
        # tset = Set{Int}()
        for m = 1:M
            for n = 1:N
                wb[n] = wf[n,t] + logpdf(df, xb[m,t+1], f(xf[n,t],u[t],t), t)
            end
            i = draw_one_categorical(pf,wb)
            # push!(tset, i)
            xb[m,t] = xf[i, t]
        end
        # @show tset
    end
    return xb,ll
end

function sse(f::AbstractFilter, u, y, λ=1)
    reset!(f)
    ll = sum(zip(u, y)) do (u,y)
        ll, e = f(u,y)
        dot(e, λ, e)
    end
end

"""
    ll = loglik(filter,u,y)
Calculate loglikelihood for entire sequences `u,y`
"""
function loglik(f::AbstractFilter,u,y)
    reset!(f)
    ll = sum(x->f(x[1],x[2])[1], zip(u, y))
end

function loglik(pf::AuxiliaryParticleFilter,u,y)
    reset!(pf)
    ll = sum(t->pf(u[t],y[t],y[t+1],t)[1], 1:length(u)-1)
    ll + pf.pf(u[end],y[end], length(u))[1]
end

"""
    ll(θ) = log_likelihood_fun(filter_from_parameters(θ::Vector)::Function, priors::Vector{Distribution}, u, y)

returns function θ -> p(y|θ)p(θ)
"""
function log_likelihood_fun(filter_from_parameters,priors::Vector{<:Distribution},u,y)
    n = numargs(filter_from_parameters)
    pf = nothing
    function (θ)
        pf === nothing && (pf = filter_from_parameters(θ))
        length(θ) == length(priors) || throw(ArgumentError("Input must have same length as priors"))
        ll = sum(i->logpdf(priors[i], θ[i]), eachindex(priors))
        isfinite(ll) || return -Inf
        pf = filter_from_parameters(θ,pf)
        ll += loglik(pf,u,y)
    end
end

function naive_sampler(θ₀)
    !any(iszero.(θ₀)) || throw(ArgumentError("Naive sampler does not work if initial parameter vector contains zeros (it was going to return θ -> θ .+ rand(MvNormal(0.1abs.(θ₀))), but that is not a good ideas if θ₀ is zero."))
    θ -> θ .+ rand(MvNormal(0.1abs.(θ₀)))
end

"""
    metropolis(ll::Function(θ), R::Int, θ₀::Vector, draw::Function(θ) = naive_sampler(θ₀))

Performs MCMC sampling using the marginal Metropolis (-Hastings) algorithm
`draw = θ -> θ'` samples a new parameter vector given an old parameter vector. The distribution must be symmetric, e.g., a Gaussian. `R` is the number of iterations.
See `log_likelihood_fun`

# Example:
```julia
filter_from_parameters(θ) = ParticleFilter(N, dynamics, measurement, MvNormal(n,exp(θ[1])), MvNormal(p,exp(θ[2])), d0)
priors = [Normal(0,0.1),Normal(0,0.1)]
ll     = log_likelihood_fun(filter_from_parameters,priors,u,y,1)
θ₀ = log.([1.,1.]) # Initial point
draw = θ -> θ .+ rand(MvNormal(0.1ones(2))) # Function that proposes new parameters (has to be symmetric)
burnin = 200 # If using threaded call, provide number of burnin iterations
# @time theta, lls = metropolis(ll, 2000, θ₀, draw) # Run single threaded
# thetam = reduce(hcat, theta)'
@time thetalls = LowLevelParticleFilters.metropolis_threaded(burnin, ll, 5000, θ₀, draw) # run on all threads, will provide (2000-burnin)*nthreads() samples
histogram(exp.(thetalls[:,1:2]), layout=3)
plot!(thetalls[:,3], subplot=3) # if threaded call, log likelihoods are in the last column
```
"""
function metropolis(ll, R, θ₀, draw = naive_sampler(θ₀))
    params    = Vector{typeof(θ₀)}(undef,R)
    lls       = Vector{Float64}(undef,R)
    params[1] = θ₀
    lls[1]    = ll(θ₀)
    for i = 2:R
        θ = draw(params[i-1])
        lli = ll(θ)
        if rand() < exp(lli-lls[i-1])
            params[i] = θ
            lls[i] = lli
        else
            params[i] = params[i-1]
            lls[i] = lls[i-1]
        end
    end
    params, lls
end

function metropolis_threaded(burnin, args...)
    res = []
    mtx = ReentrantLock()
    Threads.@threads for i = 1:Threads.nthreads()
        p,l = metropolis(args...)
        resi = [reduce(hcat,p)' l]
        resi = resi[burnin+1:end,:]
        lock(mtx)
        push!(res, resi)
        unlock(mtx)
    end
    reduce(vcat,res)
end

"""
    smoothed_mean(xb)
Helper function to calculate the mean of smoothed particle trajectories
"""
function smoothed_mean(xb)
    M,T = size(xb)
    n = length(xb[1])
    xbm = vec(mean(xb,dims=1))
    reduce(hcat, xbm)
end

"""
    smoothed_cov(xb)
Helper function to calculate the covariance of smoothed particle trajectories
"""
function smoothed_cov(xb)
    M,T = size(xb)
    n = length(xb[1])
    xbc = [cov(copy(reshape(reinterpret(Float64, xb[:,t]),n,M)),dims=2) for t = 1:T]
end

"""
    smoothed_trajs(xb)
Helper function to get particle trajectories as a 3-dimensions array (N,M,T) instead of matrix of vectors.
"""
function smoothed_trajs(xb)
    M,T = size(xb)
    n = length(xb[1])
    copy(reshape(reinterpret(Float64, xb), n,M,T))
end


text:
"""
# module Default

# Examples

```jldoctest
julia>
```
"""
module Default

export QE_PROGRAMS,
    ALLOWED_CARDS,
    ALLOWED_NAMELISTS

const QE_PROGRAMS = (
    "bands.x",
    "bgw2pw.x",
    "cp.x",
    "cppp.x",
    "dos.x",
    "turbo_eels.x",
    "turbo_lanczos.x",
    "ld1.x",
    "molecularpdos.x",
    "neb.x",
    "ph.x",
    "pp.x",
    "projwfc.x",
    "pw.x",
    "pw2bgw.x",
    "pw2gw.x",
    "pw_export.x",
    "pwcond.x",
    "turbo_spectrum.x"
)

const ALLOWED_NAMELISTS = Dict(
    "pw.x" => ("CONTROL", "SYSTEM", "ELECTRONS", "IONS", "CELL")
)

const ALLOWED_CARDS = Dict(
    "pw.x" => ("ATOMIC_SPECIES", "ATOMIC_POSITIONS", "K_POINTS", "CELL_PARAMETERS", "CONSTRAINTS", "OCCUPATIONS", "ATOMIC_FORCES")
)

end


text:
"""
**FUNÇÃO ÁREA TOTAL**
Sintaxe:
areatotal(raio, altura_cilindro)
"""

function areatotal(raio, altura)
    println("Área total:")
    return arealateral(raio, altura) + 2*areaBC(raio)
end


text:
import Polymer

@testset "utils.jl" begin
    @test unicodesymbol2string('ϕ') == "phi"
    @test unicodesymbol2string('β') == "beta"
    @test unicodesymbol2string('b') == "b"

    d = Dict("a"=>1, "b"=>2, "c"=>3)
    @test reverse_dict(d) == Dict(1=>"a", 2=>"b", 3=>"c")
    d = Dict("c"=>1, "a"=>1, "b"=>2)
    @test reverse_dict(d) == Dict(1=>"a", 2=>"b")
    d = Dict("a"=>1, "b"=>2, "c"=>1)
    @test reverse_dict(d) == Dict(1=>"a", 2=>"b")

    @test Polymer._sort_tuple2((2,1)) == (1,2)
    @test Polymer._sort_tuple2((1,2)) == (1,2)
end


text:
using Fmt
using Printf
using Formatting

fmt_print(out, x, y)        = print(out, f"({$x}, {$y})\n")
sprintf_print(out, x, y)    = print(out, @sprintf("(%d, %d)\n", x, y))
naive_print(out, x, y)      = print(out, '(', x, ", ", y, ")\n")
string_print(out, x, y)     = print(out, "($x, $y)\n")
const expr = FormatExpr("({1}, {2})\n")
formatting_print(out, x, y) = print(out, format(expr, x, y))

function benchmark(printer, out, x, y)
    @assert length(x) == length(y)
    for i in 1:length(x)
        printer(out, x[i], y[i])
    end
end

using Random
Random.seed!(1234)
x = rand(-999:999, 1_000_000)
y = rand(-999:999, 1_000_000)

using BenchmarkTools
for printer in [fmt_print, sprintf_print, naive_print,
                string_print, formatting_print]
    print(f"{$printer:>20}:")
    @btime benchmark($printer, $devnull, $x, $y)
end


text:
"""
    plot_near_field(k0, kin, P, sp::ScatteringProblem, ui::Einc;
                        opt::FMMoptions = FMMoptions(), method = "multipole",
                        x_points = 201, y_points = 201, border = find_border(sp),
                        normalize = 1.0)

Plots the total electric field as a result of a plane wave with incident
TM field `ui` scattering from the ScatteringProblem `sp`, using matplotlib's
`pcolormesh`. Can accept number of sampling points in each direction plus
bounding box or calculate automatically.

Uses the FMM options given by `opt` (FMM is disabled by default);
`method = "multipole"` dictates whether electric field is calculated using the
multipole/cylindrical harmonics, uses a faster but less accurate Hankel recurrence
formula (`"recurrence"`), or falls back on potential densities (`"density"`).
Either way, the multiple-scattering system is solved in the cylindrical
harmonics space. Normalizes all distances and sizes in plot (but not output) by
`normalize`.

Returns the calculated field in two formats:
1. `(points, Ez)` where `Ez[i]` is the total electric field at `points[i,:]`, and
2. `(xgrid,ygrid,zgrid)`, the format suitable for `pcolormesh`, where `zgrid[i,j]`
contains the field at `(mean(xgrid[i, j:j+1]), mean(ygrid[i:i+1, j]))`.
"""
function plot_near_field(k0, kin, P, sp::ScatteringProblem, ui::Einc;
                        opt::FMMoptions = FMMoptions(), method = "multipole",
                        x_points = 201, y_points = 201, border = find_border(sp),
                        normalize = 1.0)

    x_min, x_max, y_min, y_max = border

    # with pcolormesh, the field is calculated at the center of each rectangle.
    # thus we need two grids - the rectangle grid and the sampling grid.
    x = range(x_min, stop=x_max, length=x_points + 1)
    y = range(y_min, stop=y_max, length=y_points + 1)
    xgrid = repeat(transpose(x), y_points + 1)
    ygrid = repeat(y, 1, x_points + 1)
    dx = (x_max - x_min)/2/x_points
    dy = (y_max - y_min)/2/y_points
    points = cat(vec(xgrid[1:y_points, 1:x_points]) .+ dx,
                vec(ygrid[1:y_points, 1:x_points]) .+ dy, dims=2)

    Ez = calc_near_field(k0, kin, P, sp, points, ui,
            method = method, opt = opt)
    zgrid = reshape(Ez, y_points, x_points)
    figure()
    pcolormesh(xgrid/normalize, ygrid/normalize, abs.(zgrid))

    ax = gca()
    draw_shapes(sp.shapes, sp.ids, sp.centers, sp.φs, ax = ax, normalize = normalize)
    xlim([x_min/normalize;x_max/normalize])
    ylim([y_min/normalize;y_max/normalize])
    tight_layout()
    ax.set_aspect("equal", adjustable = "box")
    return (points,Ez),(xgrid,ygrid,zgrid)
end

"""
    plot_far_field(k0, kin, P, sp::ScatteringProblem, pw::PlaneWave;
                        opt::FMMoptions = FMMoptions(), method = "multipole",
                        plot_points = 200)

Plots and returns the echo width (radar cross section in two dimensions) for a
given scattering problem. `opt`, `method` are as in `plot_near_field`.
"""
function plot_far_field(k0, kin, P, sp::ScatteringProblem, pw::PlaneWave;
                    opt::FMMoptions = FMMoptions(), method = "multipole",
                    plot_points = 200)

    Rmax = maximum(s.R for s in sp.shapes)

    x_max,y_max = maximum(sp.centers, dims=1) + 2*Rmax
    x_min,y_min = minimum(sp.centers, dims=1) - 2*Rmax
    Raggregate = 0.5*max(x_max - x_min, y_max - y_min) #radius of bounding circle
    x_center = 0.5*(x_max + x_min)
    y_center = 0.5*(y_max + y_min)
    Rfar = Raggregate*1e6
    theta_far = range(0, stop=2π, length=plot_points)
    x_far = x_center + Rfar*cos.(theta_far)
    y_far = y_center + Rfar*sin.(theta_far)
    points = [x_far y_far]

    Ez = calc_far_field(k0, kin, P, points, sp, pw,
            method = method, opt = opt)
    Ez[:] = (k0*Rfar)*abs2.(Ez)
    #plot echo width
    figure()
    plot(theta_far/π, Ez)
    xlabel("\$ \\theta/\\pi \$")
    ylabel("\$ \\sigma/\\lambda_0 \$")
    title("Echo Width")
    tight_layout()
    xlim([0;2])
    return Ez
end

"""
    draw_shapes(shapes, ids, centers, φs; ax = gca(), normalize = 1.0)
    draw_shapes(sp; ax = gca(), normalize = 1.0)

Draws all of the shapes in a given scattering problem in the PyPlot axis 'ax'.
Parametrized shapes are drawn as polygons while circles are drawn using
matplotlib's `patch.Circle`. Divides all lengths by 'normalize'.
"""
function draw_shapes(shapes, ids, centers, φs; ax = gca(), normalize = 1.0)
    #draw shapes
    for ic = 1:size(centers,1)
        if typeof(shapes[ids[ic]]) == ShapeParams
            if φs[ic] == 0.0
                ft_rot = shapes[ids[ic]].ft .+ centers[ic,:]'
            else
                Rot = cartesianrotation(φs[ic])'
                ft_rot = shapes[ids[ic]].ft*Rot .+ centers[ic,:]'
            end
            ax.plot([ft_rot[:,1];ft_rot[1,1]]/normalize,
                    [ft_rot[:,2];ft_rot[1,2]]/normalize, "k", linewidth = 2)
        else
            ax.add_patch(patch.Circle((centers[ic,1]/normalize,
                            centers[ic,2]/normalize),
                            radius = shapes[ids[ic]].R/normalize,
                            edgecolor="k", facecolor="none", linewidth = 2))
        end
    end
end

draw_shapes(sp; ax = gca(), normalize = 1.0) = draw_shapes(sp.shapes,
                    sp.ids, sp.centers, sp.φs; ax = ax, normalize = normalize)

"""
    calc_near_field(k0, kin, P, sp::ScatteringProblem, points, ui::Einc;
                            opt::FMMoptions = FMMoptions(), method = "multipole",
                            verbose = true)

Calculates the total electric field as a result of a plane wave with incident
field `ui` scattering from the ScatteringProblem `sp`, at `points`.
Uses the FMM options given by `opt` (default behaviour is disabled FMM);
`method = "multipole"` dictates whether electric field is calculated using the
multipole/cylindrical harmonics, uses a faster but less accurate Hankel recurrence
formula (`"recurrence"`), or falls back on potential densities (`"density"`).
Either way, the multiple-scattering system is solved in the cylindrical
harmonics space, and the field by a particular scatterer inside its own scattering
discs is calculated by potential densities, as the cylindrical harmonics
approximation is not valid there.
"""
function calc_near_field(k0, kin, P, sp::ScatteringProblem, points, ui::Einc;
                            opt::FMMoptions = FMMoptions(), method = "multipole",
                            verbose = true)

    shapes = sp.shapes;	ids = sp.ids; centers = sp.centers; φs = sp.φs
    u = zeros(Complex{Float64},size(points,1))
    if opt.FMM
        result,sigma_mu =  solve_particle_scattering_FMM(k0, kin, P, sp,
                            ui, opt, verbose = verbose)
        if result[2].isconverged == false
            @warn("FMM process did not converge")
            return
        end
        beta = result[1]
    else
        beta, sigma_mu = solve_particle_scattering(k0, kin, P, sp, ui,
                            verbose = verbose)
    end

    #first, let's mark which points are in which shapes in tags:
    #0 denotes outside everything, +-i means inside shape i or between it and its "multipole disk"
    dt_tag = @elapsed begin
        tags = tagpoints(sp, points)
    end

    dt_in = @elapsed begin
        rng_in = zeros(Bool,size(points,1))
        rng_out = zeros(Bool,size(points,1))
        Rot = Array{Float64}(undef, 2,2)
        for ic = 1:size(sp)
            rng_in[:] = (tags .== ic)
            rng_out[:] = (tags .== -ic)
            (any(rng_in) || any(rng_out)) || continue
            if typeof(shapes[ids[ic]]) == ShapeParams
                if φs[ic] == 0.0
                    ft_rot = shapes[ids[ic]].ft .+ centers[ic,:]'
                    dft_rot = shapes[ids[ic]].dft
                else
                    Rot[:] = cartesianrotation(φs[ic])'
                    ft_rot = shapes[ids[ic]].ft*Rot .+ centers[ic,:]'
                    dft_rot = shapes[ids[ic]].dft*Rot
                end
            end
            #field inside shape
            if any(rng_in)
                if typeof(shapes[ids[ic]]) == ShapeParams
                    u[rng_in] += scatteredfield(sigma_mu[ic], kin, shapes[ids[ic]].t,
                                    ft_rot, dft_rot, points[rng_in,:])
                else
                    u[rng_in] += innerFieldCircle(kin, sigma_mu[ic], centers[ic,:],
                                    points[rng_in,:])
                end
            end
            #field between shape and multipole disk (impossible for circle)
            if any(rng_out)
                u[rng_out] += scatteredfield(sigma_mu[ic], k0, shapes[ids[ic]].t,
                                ft_rot, dft_rot, points[rng_out,:])
                u[rng_out] += uinc(k0, points[rng_out,:], ui)
                for ic2 = 1:size(sp)
                    ic == ic2 && continue
                    if method == "multipole"
                        scattered_field_multipole!(u, k0, beta, P, centers, ic2, points,
                            findall(rng_out))
                    elseif method == "recurrence"
                        scattered_field_multipole_recurrence!(u, k0, beta, P, centers, ic2, points,
                            findall(rng_out))
                    else
                        if φs[ic2] == 0.0
                            ft_rot2 = shapes[ids[ic2]].ft .+ centers[ic2,:]'
                            dft_rot2 = shapes[ids[ic2]].dft
                        else
                            Rot[:] = cartesianrotation(φs[ic2])'
                            ft_rot2 = shapes[ids[ic2]].ft*Rot .+ centers[ic2,:]'
                            dft_rot2 = shapes[ids[ic2]].dft*Rot
                        end
                        u[rng_out] += scatteredfield(sigma_mu[ic2], k0, shapes[ids[ic2]].t,
                                        ft_rot2, dft_rot2, points[rng_out,:])
                    end
                end
            end
        end
    end

    #now compute field outside all shapes
    dt_out = @elapsed begin
        rng = (tags .== 0)
        #incident field
        u[rng] = uinc(k0, points[rng,:], ui)
        if method == "multipole"
            scattered_field_multipole!(u, k0, beta, P, centers, 1:size(sp), points, findall(rng))
        elseif method == "recurrence"
            scattered_field_multipole_recurrence!(u, k0, beta, P, centers, 1:size(sp), points, findall(rng))
        else
            for ic = 1:size(centers,1)
                if typeof(shapes[ids[ic]]) == ShapeParams
                    if φs[ic] == 0.0
                        ft_rot = shapes[ids[ic]].ft .+ centers[ic,:]'
                        u[rng] += scatteredfield(sigma_mu[ic], k0, shapes[ids[ic]].t,
                                    ft_rot, shapes[ids[ic]].dft, points[rng,:])
                    else
                        Rot[:] = copy(transpose(cartesianrotation(φs[ic])))
                        ft_rot = shapes[ids[ic]].ft*Rot .+ centers[ic,:]'
                        dft_rot = shapes[ids[ic]].dft*Rot
                        u[rng] += scatteredfield(sigma_mu[ic], k0, shapes[ids[ic]].t,
                                    ft_rot, dft_rot, points[rng,:])
                    end
                else
                    scattered_field_multipole!(u, k0, beta, P, centers, ic, points,
                        findall(rng))
                end
            end
        end
    end
    if verbose
        println("Time spent calculating field:")
        println("Location tagging: $dt_tag")
        println("In/around scatterers: $dt_in")
        println("Outside scatterers: $dt_out")
    end
    return u
end

function calc_far_field(k0, kin, P, points, sp::ScatteringProblem, pw::PlaneWave;
                        opt::FMMoptions = FMMoptions(), method = "multipole")
    #calc only scattered field + assumes all points are outside shapes
    shapes = sp.shapes; centers = sp.centers; ids = sp.ids; φs = sp.φs
    if opt.FMM
        result,sigma_mu =  solve_particle_scattering_FMM(k0, kin, P, sp, pw, opt)
        if result[2].isconverged == false
            @warn("FMM process did not converge")
            return
        end
        beta = result[1]
    else
        beta, sigma_mu = solve_particle_scattering(k0, kin, P, sp, pw)
    end
    Ez = zeros(Complex{Float64}, size(points,1))

    if method == "multipole"
        scattered_field_multipole!(Ez, k0, beta, P, centers, 1:size(sp),
            points, 1:size(points,1))
    elseif method == "recurrence"
        scattered_field_multipole_recurrence!(Ez, k0, beta, P, centers,
            1:size(sp), points, 1:size(points,1))
    else
        for ic = 1:size(sp)
            if typeof(shapes[ids[ic]]) == ShapeParams
                if φs[ic] == 0.0
                    ft_rot = shapes[ids[ic]].ft .+ centers[ic,:]'
                    Ez[:] += scatteredfield(sigma_mu[ic], k0, shapes[ids[ic]].t,
                                ft_rot, shapes[ids[ic]].dft, points)
                else
                    Rot = copy(transpose(cartesianrotation(φs[ic])))
                    ft_rot = shapes[ids[ic]].ft*Rot .+ centers[ic,:]'
                    dft_rot = shapes[ids[ic]].dft*Rot
                    Ez[:] += scatteredfield(sigma_mu[ic], k0, shapes[ids[ic]].t,
                                ft_rot, dft_rot, points)
                end
            else
                warning("should be only subset of beta!")
                scattered_field_multipole!(Ez, k0, beta, P, centers, ic, points,
                    1:size(points,1))
            end
        end
    end
    return Ez
end

function tagpoints_old(sp, points)
    shapes = sp.shapes;	ids = sp.ids; centers = sp.centers; φs = sp.φs

    tags = zeros(Integer, size(points,1))
    X = Array{Float64}(undef, 2)
    for ix = 1:size(points,1)
        for ic = 1:size(sp)
            X .= points[ix,:] - centers[ic,:]
            if sum(abs2,X) <= shapes[ids[ic]].R^2
                if typeof(shapes[ids[ic]]) == ShapeParams
                    if φs[ic] != 0.0 #rotate point backwards instead of shape forwards
                        Rot = [cos(-φs[ic]) -sin(-φs[ic]);sin(-φs[ic]) cos(-φs[ic])]
                        X = Rot*X
                    end
                    tags[ix] = pInPolygon(X, shapes[ids[ic]].ft) ? ic : -ic
        			break #can't be in two shapes
                else #CircleParams
                    tags[ix] = ic
                    break #can't be in two shapes
                end
        	end
        end
    end
    tags
end

function tagpoints(sp, points)
    shapes = sp.shapes;	ids = sp.ids; centers = sp.centers; φs = sp.φs

    tags = zeros(Integer, size(points,1))
    X = Array{Float64}(undef, 2) #tmp arrays
    for ix = 1:size(points,1) #need two loops due to "break"
        for ic = 1:size(sp)
            X[1] = points[ix,1] - centers[ic,1]
            X[2] = points[ix,2] - centers[ic,2]
            if hypot(X[1], X[2]) ≤ shapes[ids[ic]].R
                if isa(shapes[ids[ic]], ShapeParams)
                    if φs[ic] != 0.0 #rotate point backwards instead of shape forwards
                        Rot = cartesianrotation(-φs[ic])
                        X = Rot*X
                    end
                    tags[ix] = pInPolygon(X, shapes[ids[ic]].ft) ? ic : -ic
                else #CircleParams
                    tags[ix] = ic
                end
                break #can't be in two shapes
        	end
        end
    end
    tags
end


text:
const pyspatial = PyNULL()
scipy = false
function __init__()
    copy!(pyspatial,  pyimport("scipy.spatial"))
    # is scipy available?
    global scipy = pyspatial == PyNULL() ? false : true
    return
end
#-------------------------------------------------------------------------------
function edgestoloop(el::Array{Int64, 2})
    elcol = reshape(el, length(el), 1)
    listep, I, J = uniquerows(elcol)
    J = reshape(J, size(el))
    ptoe = zeros(Int64, length(I), 2)
    for k = 1:size(el, 1)
        for l = 1:2
            if ptoe[J[k, l], 1] == 0
                ptoe[J[k, l], 1] = k
            else
                ptoe[J[k, l], 2] = k
            end
        end
    end
    if minimum(ptoe) == 0
        error("Pb I in edgestoloop")
    end
    loop = zeros(Int64, length(I)); loop[1:2] = J[1, :]; ec = 1 # courant edge
    for k = 3:length(I)
        ec = ec == ptoe[loop[k - 1], 1] ? ptoe[loop[k - 1], 2] : ptoe[loop[k - 1], 1]
        #ec = setdiff(ptoe[loop[k - 1], :], ec) # actualize courant edge
        # setdiff(J[ec, :], loop[k - 1])[1]
        loop[k] = J[ec, 1] == loop[k - 1] ? J[ec, 2] : J[ec, 1]
    end
    if minimum(loop) == 0
        error("Pb II in edgestoloop")
    end
    return elcol[I[loop]]
end
#-------------------------------------------------------------------------------
"""
    Tri, pyhull = convexhull(p::Array{Float64, 2},
                             qhull_option::String = "QJ Pp")
    convexhull()
Call spatial scipy library to compute convehull with qhull.

WARNING : there is perhaps a problem of orientation of triangles!
"""
function convexhull(p::Array{Float64, 2}, qhull_option::String = "QJ Pp")
    pyhull = pyspatial.ConvexHull(p, qhull_options = qhull_option)
    Tri = convert(Array{Int64, 2}, pyhull.simplices) .+ 1
    return Tri, pyhull
end
#-------------------------------------------------------------------------------
"""
    Au, Iu, Ju = uniquerows(A::AbstractArray{T}, preci::Int64 = 6)
Matlab like unique_rows function

N.B.
* use preci = 0 for integers
* maximum(abs.((A[Iu, :])[Ju, :] .- A)) == 0
"""
function uniquerows(A::AbstractArray{T}, preci::Int64 = 0) where T
    Ar = if preci > 0 # use preci = 0 for integers
        round.(Int64, A * 10 ^ preci)
    else
        A
    end
    Ju = groupslices(Ar, dims = 1)
    p = sortperm(Ju)
    Iu, c = [Ju[p[1]],], 1
    Jun = similar(Ju)
    Jun[p[1]] = c
    for k = 2:length(Ju)
        if Ju[p[k]] > Ju[p[k - 1]]
            push!(Iu, Ju[p[k]])
            c += 1
        end
        Jun[p[k]] = c
    end
    #@test maximum(abs.((Ar[Iu, :])[Jun, :] .- Ar)) == 0
    return A[Iu, :], Iu, Jun
end


text:
@testset "Differentials" begin
    @testset "Wirtinger" begin
        w = Wirtinger(1+1im, 2+2im)
        @test wirtinger_primal(w) == 1+1im
        @test wirtinger_conjugate(w) == 2+2im
        @test w + w == Wirtinger(2+2im, 4+4im)

        @test w + One() == w + 1 == w + Thunk(()->1) == Wirtinger(2+1im, 2+2im)
        @test w * One() == One() * w == w
        @test w * 2 == 2 * w == Wirtinger(2 + 2im, 4 + 4im)

        # TODO: other + methods stack overflow
        @test_throws ErrorException w*w
        @test_throws ArgumentError extern(w)
        for x in w
            @test x === w
        end
        @test broadcastable(w) == w
        @test_throws MethodError conj(w)
    end
    @testset "Zero" begin
        z = Zero()
        @test extern(z) === false
        @test z + z == z
        @test z + 1 == 1
        @test 1 + z == 1
        @test z * z == z
        @test z * 1 == z
        @test 1 * z == z
        for x in z
            @test x === z
        end
        @test broadcastable(z) isa Ref{Zero}
        @test conj(z) == z
    end
    @testset "One" begin
        o = One()
        @test extern(o) === true
        @test o + o == 2
        @test o + 1 == 2
        @test 1 + o == 2
        @test o * o == o
        @test o * 1 == 1
        @test 1 * o == 1
        for x in o
            @test x === o
        end
        @test broadcastable(o) isa Ref{One}
        @test conj(o) == o
    end

    @testset "No ambiguities in $f" for f in (+, *)
        # We don't use `Test.detect_ambiguities` as we are only interested in
        # the +, and * operations. We also would catch any that are unrelated
        # to this package. but that is not a problem. Since no such failings
        # occur in our dependencies.

        ambig_methods = [
            (m1, m2) for m1 in methods(f), m2 in methods(f) if Base.isambiguous(m1, m2)
        ]
        @test isempty(ambig_methods)
    end
end


text:
function solve!(solver::AugmentedLagrangianSolver{T,S}) where {T,S}
	initialize!(solver)
    c_max::T = Inf

	conSet = get_constraints(solver)
	solver_uncon = solver.solver_uncon::S

	# Calculate cost
    J_ = TO.get_J(get_objective(solver))
    J = sum(J_)

    # Terminal tolerances
    cost_tol = solver.opts.cost_tolerance
    grad_tol = solver.opts.gradient_tolerance

    for i = 1:solver.opts.iterations_outer
		set_tolerances!(solver, solver_uncon, i, cost_tol, grad_tol)

        # Solve the unconstrained problem
        solve!(solver.solver_uncon)

        # Check solver status
        status(solver) > SOLVE_SUCCEEDED && break

        # Record the updated information
        J = sum(J_)
        TO.max_violation!(conSet)
        c_max = maximum(conSet.c_max)
        record_iteration!(solver, J, c_max)

        # Check for convergence before doing the outer loop udpate
        converged = evaluate_convergence(solver)
        if converged
            break
        end

        # Outer loop update
        dual_update!(solver)
        penalty_update!(solver)

        # Reset verbosity level after it's modified
        set_verbosity!(solver)

        reset!(solver_uncon)

        if i == solver.opts.iterations_outer
            solver.stats.status = MAX_ITERATIONS
        end
    end
    solver.opts.cost_tolerance = cost_tol
    solver.opts.gradient_tolerance = grad_tol
    terminate!(solver)
    return solver
end

function initialize!(solver::AugmentedLagrangianSolver)
	set_verbosity!(solver)
	clear_cache!(solver)

	# Reset solver
    reset!(solver)

	# Calculate cost
    TO.cost!(get_objective(solver), get_trajectory(solver))
end

function step!(solver::AugmentedLagrangianSolver)

    # Solve the unconstrained problem
    solve!(solver.solver_uncon)

    # Outer loop update
    dual_update!(solver)
    penalty_update!(solver)
    TO.max_violation!(get_constraints(solver))

	# Reset verbosity level after it's modified
	set_verbosity!(solver)
end

function record_iteration!(solver::AugmentedLagrangianSolver{T,S}, J::T, c_max::T) where {T,S}

	conSet = get_constraints(solver)
	max_penalty!(conSet)
    max_penalty = maximum(conSet.μ_max)
    J_prev = solver.stats.cost[solver.stats.iterations]
    dJ = J_prev - J

    # Just update constraint violation and max penalty
    record_iteration!(solver.stats, c_max=c_max, penalty_max=max_penalty, is_outer=true)
    j = solver.stats.iterations_outer::Int

	if is_verbose(solver)
        @logmsg OuterLoop :iter value=j
        @logmsg OuterLoop :total value=solver.stats.iterations
        @logmsg OuterLoop :cost value=J
        @logmsg OuterLoop :c_max value=c_max
        @logmsg OuterLoop :penalty value=max_penalty
		print_level(OuterLoop, global_logger())
	end
end

function set_tolerances!(solver::AugmentedLagrangianSolver{T},
        solver_uncon::AbstractSolver{T}, i::Int,
        cost_tol=solver.opts.cost_tolerance,
        grad_tol=solver.opts.gradient_tolerance
    ) where T
    if i != solver.opts.iterations_outer
        solver_uncon.opts.cost_tolerance = solver.opts.cost_tolerance_intermediate
        solver_uncon.opts.gradient_tolerance = solver.opts.gradient_tolerance_intermediate
    else
        solver_uncon.opts.cost_tolerance = cost_tol
        solver_uncon.opts.gradient_tolerance = grad_tol
    end

    return nothing
end

function evaluate_convergence(solver::AugmentedLagrangianSolver)
	i = solver.stats.iterations
    solver.stats.c_max[i] < solver.opts.constraint_tolerance ||
		solver.stats.penalty_max[i] >= solver.opts.penalty_max
end

"General Dual Update"
function dual_update!(solver::AugmentedLagrangianSolver) where {T,Q,N,M,NM}
    conSet = get_constraints(solver)
	dual_update!(conSet)
end

"General Penalty Update"
function penalty_update!(solver::AugmentedLagrangianSolver)
    conSet = get_constraints(solver)
	penalty_update!(conSet)
end


text:
include("tests.jl")
include("test_utils.jl")


text:
"""
    n"Fe"

Implements compile time parsing of strings to produce Element, SubShell, AtomicSubShell,
Transition or CharXRay objects.  The only oddity is that to get SubShell(\"K\") you
must enter n\"K1\" to differentiate the sub-shell from potassium.
Examples:

    n"Fe" == elements[26]
    n"K" == elements[19]
    n"K1" == subshell("K")
    n"L3" == subshell("L3")
    n"Fe L3" == AtomicSubShell(elements[26],subshell("L3"))
    n"L3-M5" == Transition(SubShell("L3"),"SubShell("M5"))
    n"Fe L3-M5" == CharXRay(elements[26],Transition(SubShell("L3"),"SubShell("M5"))
"""
macro n_str(str)
    parsex(str)
end

"""
    enx"Fe K"
    enx"Fe L3-M5"

Implements compile time parsing of strings to produce the energy associated with AtomicSubShell,
or CharXRay objects.  The advantage of enx"Fe K" over `energy(n"Fe K")` is that the evaluation occurs
entirely at compile time with the former while the later still needs to call `energy(...)` at runtime.

Examples:

    enx"Fe L3" == energy(AtomicSubShell(element(26),subshell("L3")))
    enx"Fe L3-M5" == energy(CharXRay(element(26),Transition(SubShell("L3"),"SubShell("M5")))
"""

macro enx_str(str)
    energy(parsex(str))
end

Base.parse(::Type{Element}, str::AbstractString) = element(str)

"""
    element(str::AbstractString)

Parses a string to determine if the string represents an Element by atomic number, symbol or name.
"""
function element(str::AbstractString)::Element
    elm = get(PeriodicTable.elements, str,
            get(PeriodicTable.elements, Symbol(str), missing))
    if ismissing(elm)
        try
            elm = PeriodicTable.elements[parse(Int, str)]
        catch
            error("Unable to parse ", str, " as an Element.")
        end
    end
    return elm
end

"""
    subshell(name::AbstractString)::SubShell

Returns a SubShell structure from a string of the form "K", "L1", ...., "O11"
"""
function subshell(name::AbstractString)
    i = findfirst(ss->repr(ss)==name, allsubshells)
    return !isnothing(i) ? allsubshells[i] : error("Unknown sub-shell $(name)")
end

Base.parse(::Type{SubShell}, name::AbstractString) = subshell(name)

"""
    Base.parse(::Type{SubShell}, name::AbstractString)
    subshell(name::AbstractString)

Parse an SubShell from a string of the form "K", "L1", "L2", ...

    Base.parse(::Type{AtomicSubShell}, str::AbstractString)::AtomicSubShell
    atomicsubshell(str::AbstractString)::AtomicSubShell

Parse an AtomicSubShell from a string of the form \"Fe K\" or \"U M5\".
"""
function Base.parse(::Type{AtomicSubShell}, str::AbstractString)::AtomicSubShell
    sp = split(str, " ")
    if length(sp) == 2
        sh = subshell(sp[2])
        elm = element(sp[1])
        return AtomicSubShell(z(elm), sh)
    end
    error("Cannot parse $str as an AtomicSubShell like \"Fe L3\"")
end

"""
    atomicsubshell(str::AbstractString)::AtomicSubShell

Parse an AtomicSubShell from a string of the form \"Fe K\" or \"U M5\".
"""
atomicsubshell(str::AbstractString) = Base.parse(AtomicSubShell, str)

"""
    characteristic(str::AbstractString)::CharXRay

Create a CharXRay structure from a string like \"Fe K-L3\" or \"U L3-M5\".
"""
function characteristic(str::AbstractString)::CharXRay
    sp1 = split(str, " ")
    if length(sp1) == 2
        elm = element(sp1[1])
        if !ismissing(elm)
            sp2 = split(sp1[2], "-")
            if length(sp2) == 2
                inner = subshell(sp2[1])
                outer = subshell(sp2[2])
                if (!ismissing(inner)) && (!ismissing(outer))
                    return CharXRay(z(elm), transition(inner, outer))
                end
            end
        end
    end
    error("Unable to parse $str as a characteristic X-ray.")
end

Base.parse(::Type{CharXRay}, str::AbstractString)::CharXRay = characteristic(str)

"""
    parsex(str::AbstractString)::Union{Element, SubShell, AtomicSubShell, Transition, CharXRay}

Implements compile time parsing of strings to produce Element, SubShell, AtomicSubShell,
Transition or CharXRay objects. The only oddity is that to get SubShell(\"K\") you
must enter n\"K1\" to differentiate the sub-shell from potassium. `parsex(...)` provides
the implementation for `n"..."`.

Examples:

  * n"Fe" => Element
  * n"L3" => SubShell
  * n"K" => Element
  * n"K1" => SubShell ("K" would be the element potassium)
  * n"Fe K" => AtomicSubShell ("K" not ambiguous here...)
  * n"Fe L3" => AtomicSubShell
  * n"L3-M5" => Transition
  * n"Fe L3-M5" => CharXRay.
"""
function parsex(
    str::AbstractString,
)::Union{Element,SubShell,AtomicSubShell,Transition,CharXRay}
    function _impl_(str::AbstractString)
        sp1 = split(str, " ")
        if length(sp1) == 1  # Could be an Element, an SubShell or a Transition
            if length(split(sp1[1], "-")) == 2 # A transition like "L3-M5"
                try
                    return parse(Transition, str)
                catch
                    return nothing
                end
            else # An Element or a SubShell like "Fe" or "L3"
                try # Try and Element first
                    return parse(Element, str)
                catch
                    # Ignore and try to parse as a sub-shell
                end
                # If it isn't an element, it must be a SubShell
                # Note that since "K" is an element we use "K1" to denote the SubShell
                try
                    return str == "K1" ? subshell("K") : parse(SubShell, str)
                catch
                    return nothing
                end
            end
        else # Could be a AtomicSubShell or a CharXRay
            sp2 = split(sp1[2], "-")
            try
                if length(sp2) == 1  # Like "Fe L3"
                    return parse(AtomicSubShell, str)
                else  # Like "Fe L3-M5"
                    return parse(CharXRay, str)
                end
            catch
                return nothing
            end
        end
    end
    res = _impl_(str)
    isnothing(res) && throw(ArgumentError("Unable to parse \"$str\" as an Element, a SubShell, a Transition, an AtomicSubShell or a CharXRay."))
    return res
end


text:
using SafeTestsets

files = [
    "MINERVA_tests.jl",
    "REM_tests.jl"
]

for file in files
    include(file)
end