text:
using Printf
using LinearAlgebra
using SparseArrays
using Random
using Flux

# Run FedDCD with logistic regression
function RunFedDCD(
    fileTrain::String,
    fileTest::String,
    lambda::Float64,
    participationRate::Float64,
    localLr::Float64,
    numRounds::Int64,
    writeFileName::String
    )
    numClients = 100
    numRounds = numRounds

    # read data
    Xtrain, Ytrain = read_libsvm(fileTrain)
    Xtest, Ytest = read_libsvm(fileTest)
    Ytrain, Ytest = labelTransform(Ytrain, Ytest)
    # Set Xtrain and Xtest same number of feature
    Itr, Jtr, Vtr = findnz(Xtrain)
    Ite, Jte, Vte = findnz(Xtest)
    d = max( size(Xtrain, 2), size(Xtest, 2) )
    Xtrain = sparse(Itr, Jtr, Vtr, size(Xtrain, 1), d)
    Xtest = sparse(Ite, Jte, Vte, size(Xtest, 1), d)

    numClasses = length( union( Set(Ytrain), Set(Ytest) ) )
    # Split data
    Xsplit, Ysplit = splitDataByRow(Xtrain, Ytrain, numClients)
    # Xsplit, Ysplit = splitDataByClass(Xtrain, Ytrain, numClients, numClasses)

    # Setup config, running FedAvg if mu=0.
    clientConfig = Dict(
        "num_classes" => numClasses,
        "lambda" => lambda,
        "learning_rate" => localLr,
        "participation_rate" => participationRate,
    )

    serverConfig = Dict(
        "num_classes" => numClasses,
        "num_clients" => numClients,
        "participation_rate" => participationRate,
        "learning_rate" => 0.99,
    )

    # Construct clients
    clients = Vector{FedDCDClient}(undef, numClients)
    for i = 1:numClients
        clients[i] = FedDCDClient(i, Xsplit[i], Ysplit[i], clientConfig, newton!)
    end
    # Construct server
    server = FedDCDServer(Xtest, Ytest, serverConfig)

    # Train
    W, objList, testAccList = fedDCD(server, clients, numRounds)

    writeToFile(
        fileTrain,
        "softmax classification",
        serverConfig,
        clientConfig,
        objList,
        testAccList,
        writeFileName
    )

    @printf("Finish training!\n")
end

# Run AccFedDCD with logistic regression
function RunAccFedDCD(
    fileTrain::String,
    fileTest::String,
    lambda::Float64,
    participationRate::Float64,
    localLr::Float64,
    numRounds::Int64,
    writeFileName::String
    )
    numClients = 100;
    numRounds = numRounds;
    # Read data
    Xtrain, Ytrain = read_libsvm(fileTrain);
    Xtest, Ytest = read_libsvm(fileTest);
    Ytrain, Ytest = labelTransform(Ytrain, Ytest);
    # Set Xtrain and Xtest same number of feature
    Itr, Jtr, Vtr = findnz(Xtrain);
    Ite, Jte, Vte = findnz(Xtest);
    d = max( size(Xtrain, 2), size(Xtest, 2) );
    Xtrain = sparse(Itr, Jtr, Vtr, size(Xtrain, 1), d);
    Xtest = sparse(Ite, Jte, Vte, size(Xtest, 1), d);

    numClasses = length( union( Set(Ytrain), Set(Ytest) ) );
    # Split data
    Xsplit, Ysplit = splitDataByRow(Xtrain, Ytrain, numClients)
    # Xsplit, Ysplit = splitDataByClass(Xtrain, Ytrain, numClients, numClasses)

    # Setup config, running FedAvg if mu=0.
    clientConfig = Dict(
        "num_classes" => numClasses,
        "lambda" => lambda,
        "learning_rate" => localLr,
        "participation_rate" => participationRate,
    )

    serverConfig = Dict(
        "num_classes" => numClasses,
        "num_clients" => numClients,
        "participation_rate" => participationRate,
        "learning_rate" => 0.99,
    )

    # Construct clients
    clients = Vector{AccFedDCDClient}(undef, numClients)
    for i = 1:numClients
        clients[i] = AccFedDCDClient(i, Xsplit[i], Ysplit[i], clientConfig, newton!)
    end
    # Construct server
    server = AccFedDCDServer(Xtrain, Ytrain, serverConfig)

    # Train
    W, objList, testAccList = accfedDCD(server, clients, numRounds, objMin = 0.0)

    writeToFile(
        fileTrain,
        "softmax classification",
        serverConfig,
        clientConfig,
        objList,
        testAccList,
        writeFileName
    )

    @printf("Finish training!\n")
end

# Run FedDCD with neural network
function TestFedDCDNN(
    fileTrain::String,
    fileTest::String,
    lambda::Float64,
    participationRate::Float64,
    localLr::Float64,
    numRounds::Int64,
    writeFileName::String
    )
    numClients = 100;
    numRounds = numRounds;
    # Read data
    Xtrain, Ytrain = read_libsvm(fileTrain);
    Xtest, Ytest = read_libsvm(fileTest);
    Ytrain, Ytest = labelTransform(Ytrain, Ytest);
    # Set Xtrain and Xtest same number of feature
    Itr, Jtr, Vtr = findnz(Xtrain);
    Ite, Jte, Vte = findnz(Xtest);
    d = max( size(Xtrain, 2), size(Xtest, 2) );
    Xtrain = sparse(Itr, Jtr, Vtr, size(Xtrain, 1), d);
    Xtest = sparse(Ite, Jte, Vte, size(Xtest, 1), d);

    numClasses = length( union( Set(Ytrain), Set(Ytest) ) );
    # Split data
    Xsplit, Ysplit = splitDataByRow(Xtrain, Ytrain, numClients)
    # Xsplit, Ysplit = splitDataByClass(Xtrain, Ytrain, numClients, numClasses)

    # Setup config, running FedAvg if mu=0.
    clientConfig = Dict(
        "num_classes" => numClasses,
        "lambda" => lambda,
    )

    serverConfig = Dict(
        "num_classes" => numClasses,
        "num_clients" => numClients,
        "participation_rate" => participationRate,
        "learning_rate" => localLr,
        "decay_rate" => 0.0,
    )

    # model structure
    dim = 32

    # Construct clients
    clients = Vector{FedDCDClientNN}(undef, numClients)
    for i = 1:numClients
        clients[i] = FedDCDClientNN(i, Xsplit[i], Ysplit[i], dim, clientConfig, adam!)
    end

    # Construct server
    server = FedDCDServerNN(Xtest, Ytest, dim, serverConfig)

    # # Train
    W, objList, testAccList = fedDCD(server, clients, numRounds)
    @printf("Finish training!\n")

    writeToFile(
        fileTrain,
        "softmax classification with MLP",
        serverConfig,
        clientConfig,
        objList,
        testAccList,
        writeFileName
    )
    return W
end

# test neural network
function TestNN(
    fileTrain::String,
    fileTest::String
    )
    # Read data
    Xtrain, Ytrain = read_libsvm(fileTrain)
    Xtest, Ytest = read_libsvm(fileTest)
    Ytrain, Ytest = labelTransform(Ytrain, Ytest)
    # Set Xtrain and Xtest same number of feature
    Itr, Jtr, Vtr = findnz(Xtrain)
    Ite, Jte, Vte = findnz(Xtest)
    d = max( size(Xtrain, 2), size(Xtest, 2) )
    Xtrain = sparse(Itr, Jtr, Vtr, size(Xtrain, 1), d)
    Xtest = sparse(Ite, Jte, Vte, size(Xtest, 1), d)
    XtrainT = copy(Xtrain')
    XtestT = copy(Xtest')
    Ytrain = Flux.onehotbatch(Ytrain, 1:10)
    Ytest =  Flux.onehotbatch(Ytest, 1:10)
    # build model
    model = Chain( Dense(780, 32, relu, bias=false), Dense(32, 10, bias=false), NNlib.softmax)
    # data loader
    data = Flux.Data.DataLoader((XtrainT, Ytrain), batchsize=128, shuffle=true)
    # optimizer
    opt = ADAM()
    # loss
    λ = 1e-2
    sqnorm(w) = sum(abs2, w)
    loss(x, l) = Flux.crossentropy(model(x), l) + (λ/2)*sum(sqnorm, params(model))
    #training
    num_epoches = 1000
    f = open("results/exp1/optval_MNIST_MLP_lambda1e-2.txt", "w")
    for t = 1:num_epoches
        # @printf "epoch: %d, training obj: %.4f, test accuracy: %.2f\n" t loss(XtrainT, Ytrain) accuracy(Xtest, Ytest, model)
        @printf(f, "epoch: %d, training obj: %.6f\n", t, loss(XtrainT, Ytrain))
        Flux.train!(loss, params(model), data, opt)
    end
    close(f)
    return model
end

# Run AccFedDCD with neural network
function RunAccFedDCDNN(
    fileTrain::String,
    fileTest::String,
    lambda::Float64,
    participationRate::Float64,
    localLr::Float64,
    numRounds::Int64,
    writeFileName::String
    )
    numClients = 100;
    numRounds = numRounds;
    # Read data
    Xtrain, Ytrain = read_libsvm(fileTrain);
    Xtest, Ytest = read_libsvm(fileTest);
    Ytrain, Ytest = labelTransform(Ytrain, Ytest);
    # Set Xtrain and Xtest same number of feature
    Itr, Jtr, Vtr = findnz(Xtrain);
    Ite, Jte, Vte = findnz(Xtest);
    d = max( size(Xtrain, 2), size(Xtest, 2) );
    Xtrain = sparse(Itr, Jtr, Vtr, size(Xtrain, 1), d);
    Xtest = sparse(Ite, Jte, Vte, size(Xtest, 1), d);

    numClasses = length( union( Set(Ytrain), Set(Ytest) ) );
    # Split data
    # Xsplit, Ysplit = splitDataByRow(Xtrain, Ytrain, numClients)
    Xsplit, Ysplit = splitDataByClass(Xtrain, Ytrain, numClients, numClasses)

    # Setup config
    clientConfig = Dict(
        "num_classes" => numClasses,
        "lambda" => lambda,
        "participation_rate" => participationRate,
    )

    serverConfig = Dict(
        "num_classes" => numClasses,
        "num_clients" => numClients,
        "participation_rate" => participationRate,
        "learning_rate" => localLr,
    )

    # model structure
    dim = 32

    # Construct clients
    clients = Vector{AccFedDCDClientNN}(undef, numClients)
    for i = 1:numClients
        clients[i] = AccFedDCDClientNN(i, Xsplit[i], Ysplit[i], dim, clientConfig, adam!)
    end

    # Construct server
    server = AccFedDCDServerNN(Xtest, Ytest, dim, serverConfig)

    # # Train
    W, objList, testAccList = accfedDCD(server, clients, numRounds)
    @printf("Test finished!\n")

    writeToFile(
        fileTrain,
        "softmax classification with MLP",
        serverConfig,
        clientConfig,
        objList,
        testAccList,
        writeFileName
    )
    return W
end


text:
using PyPlot
using PyCall
import Base64
@pyimport matplotlib.animation as anim

function plot_wigner(ψ)
    xlim = sqrt(size(ψ.data)[1])
    xvec = LinRange(-xlim, xlim, 100) |> collect
    wig = wigner(ψ, xvec, xvec)
    maxwig = maximum(wig)
    imshow(permutedims(wig), vmax=maxwig, vmin=-maxwig, cmap="RdBu",
           extent=[-xlim, xlim, -xlim, xlim], origin="lower")
    colorbar()
end

function plot_wigner(ψ,xlim::Number)
    xvec = LinRange(-xlim, xlim, 100) |> collect
    wig = wigner(ψ, xvec, xvec)
    maxwig = maximum(wig)
    imshow(permutedims(wig), vmax=maxwig, vmin=-maxwig, cmap="RdBu",
           extent=[-xlim, xlim, -xlim, xlim], origin="lower")
    colorbar()
end

function plot_wigner_without_cb(ψ,xlim::Number)
    xvec = LinRange(-xlim, xlim, 100) |> collect
    wig = wigner(ψ, xvec, xvec)
    maxwig = maximum(wig)
    imshow(permutedims(wig), vmax=maxwig, vmin=-maxwig, cmap="RdBu",
           extent=[-xlim, xlim, -xlim, xlim], origin="lower")
end

function showanim(filename)
    base64_video = Base64.base64encode(open(filename))
    display("text/html", """<video controls src="data:video/x-m4v;base64,$base64_video">""")
end

function draw_anim_two_wigners(filename,val_frames::Int,val_interval::Int,ψs,ϕs,xlim::Number)

    function make_frame(i)
        frame_number = i * (Int(floor(length(ψs)/val_frames))-1)+1
        clf()
        subplot(121)
        plot_wigner(ψs[frame_number],xlim)
        subplot(122)
        plot_wigner(ϕs[frame_number],xlim)
    end

    fig = figure(figsize=(12,6))
    withfig(fig) do
        myanim = anim.FuncAnimation(fig, make_frame, frames=val_frames, interval=val_interval)
        myanim[:save](filename, bitrate=-1, extra_args=["-vcodec", "libx264", "-pix_fmt", "yuv420p"])
    end
    showanim(filename)
end

function draw_anim_wigner(filename,val_frames::Int,val_interval::Int,ψs,xlim::Number)

    function make_frame(i)
        frame_number = i+1
        clf()
        plot_wigner(ψs[frame_number],xlim)
    end

    fig = figure(figsize=(6,6))
    withfig(fig) do
        myanim = anim.FuncAnimation(fig, make_frame, frames=val_frames, interval=val_interval)
        myanim[:save](filename, bitrate=-1, extra_args=["-vcodec", "libx264", "-pix_fmt", "yuv420p"])
    end
    showanim(filename)
end

function draw_anim_wigner(filename,val_frames::Int,val_interval::Int,ψs,xlim::Number,cb::Bool)

    function make_frame(i)
        frame_number = i+1
        clf()
        if cb == true
            plot_wigner(ψs[frame_number],xlim)
        else
            plot_wigner_without_cb(ψs[frame_number],xlim)
        end

    end

    fig = figure(figsize=(6,6))
    withfig(fig) do
        myanim = anim.FuncAnimation(fig, make_frame, frames=val_frames, interval=val_interval)
        myanim[:save](filename, bitrate=-1, extra_args=["-vcodec", "libx264", "-pix_fmt", "yuv420p"])
    end
    showanim(filename)
end

ext: jl

max_stars_repo_name: ash-y/QuantumCOMFrame

max_issues_repo_name: ash-y/QuantumCOMFrame

max_forks_repo_name: ash-y/QuantumCOMFrame


text:
using VortexDistributions, Plots

N=1000
θ = LinRange(0,4*pi,N)
psi = exp.(im*θ) + 0.1*(randn(N)+im*randn(N))
ϕ = angle.(psi)
ϕu = unwrap(ϕ)
ϕui = similar(ϕu)
unwrap!(ϕui,ϕ)

@test ϕu == ϕui

# plot(ϕ,label="phi",c=:red,alpha=0.4,lw=6)
# plot!(ϕu,label="phiu",lw=2,c=:black)
# plot!(ϕui,label="phiu",lw=2,c=:green)


text:
using WeightedNorms

@static if VERSION < v"0.7-"
    using Compat: range
else
    using Printf
end

## Norm benchmarking
const nsamples = 100
maxNormNameLength = maximum(length(string(s)) for s in WeightedNorms.exportedNormSymbols)
nExportedNorms = length(WeightedNorms.exportedNormSymbols)
for len in [2^n for n in 6:10]
  x = rand(len)
  grid = range(0, stop=1, length=len)
  tau = step(grid)
  tt1 = ones(nExportedNorms)
  tt2 = ones(nExportedNorms)
  for (k, func) in enumerate(WeightedNorms.exportedNorms)
    func(x, grid)
    func(x, tau)
    gc()
    tt1[k] = @elapsed for i in 1:nsamples; func(x, grid); end
    gc()
    tt2[k] = @elapsed for i in 1:nsamples; func(x, tau); end
  end
  for (i, (t1, t2)) in enumerate(zip(tt1, tt2))
      funci = rpad(string(WeightedNorms.exportedNormSymbols[i]), maxNormNameLength, " ")
      @printf "n: %d, func: %s, time 1: %e, time 2: %e, grid is: %f x slower\n" len funci t1 t2 t1/t2
  end
end


text:
function save_params!(g, params::ParamsB1)
	g["N"] = params.N
	g["m"] = params.m
	g["β"] = params.β
	g["u"] = params.u
	g["u₁"] = params.u₁
end
function save_params!(g, params::ParamsB1_pinned)
	g["N"] = params.N
	g["m"] = params.m
	g["β"] = params.β
	g["u"] = params.u
	g["u₁"] = params.u₁
	g["i_pinned"] = params.i_pinned
end
function save_params!(g, params::ParamsNoB1)
	g["N"] = params.N
	g["m"] = params.m
	g["β"] = params.β
	g["u"] = params.u
end
function save_params!(g, params::ParamsNoB1_pinned)
	g["N"] = params.N
	g["m"] = params.m
	g["β"] = params.β
	g["u"] = params.u
	g["i_pinned"] = params.i_pinned
end

function save_dispersion!(g, rdisp::ReducedDispersion)
	g["α"] = rdisp.α
	g["P"] = rdisp.P
	g["μ"] = rdisp.μ
	g["Λ"] = rdisp.Λ
end

isB1(params::AbstractParamsB1) = true
isB1(params::AbstractParamsNoB1) = false

struct Saver
	save_file::String
	log_file::String
	function Saver(filestring::AbstractString, params, rdisp::ReducedDispersion, Nₚ, max_iter)
		save_file = filestring * ".h5"
		init_file(save_file, params, rdisp, Nₚ, max_iter)
		log_file = filestring * ".log"
		new(save_file, log_file)
	end
end
@inline Saver(fname::Tuple{AbstractString, AbstractString}, args...) = (mkpath(fname[1]); Saver(fname[1]*fname[2], args...))

function init_file(save_file::AbstractString, params::AbstractParams, rdisp::ReducedDispersion, Nₚ::Int64, max_iter::Int64)
	max_iter = max_iter + 1
	h5open(save_file, "w") do file
        g = g_create(file, "parameters")
        g1 = g_create(g, "general")
        g2 = g_create(g, "dispersion")
        save_params!(g1, params)
        save_dispersion!(g2, rdisp)
        g["Nₚ"] = Nₚ
        g["isB1"] = isB1(params)
		#####
		N′ = get_length(params)
		file["bss"] = zeros(max_iter, N′)
		file["∂F"] = zeros(N′)
		file["∂²F"] = zeros(N′, N′)
		file["ΔFs"] = zeros(max_iter)
		file["ϵs"] = zeros(max_iter)
	end
end

function save_data(saver_o::Saver, bs, ΔF, ∂F, ∂²F, ϵ, params, rdisp, Nₚ, iter)
	iter = iter+1
	h5open(saver_o.save_file, "r+") do file
		file["bss"][iter, :] = bs
		file["∂F"][:] = ∂F
		file["∂²F"][:,:] = ∂²F
		file["ΔFs"][iter] = ΔF
		file["ϵs"][iter] = ϵ
	end
end
function log_data(saver_o::Saver, i, t, ϵ)
	io = open(saver_o.log_file, "a")
	print(io, @sprintf("%6d %12.3E %12.2E\n", i, t, ϵ))
	close(io)
end
function output(saver_o, bs, ΔF, ∂F, ∂²F, ϵ, t, i, params, rdisp, Nₚ, iter)
	save_data(saver_o, bs, ΔF, ∂F, ∂²F, ϵ, params, rdisp, Nₚ, iter)
	log_data(saver_o, i, t, ϵ)
end

function seek_minimum(params::AbstractParams, rdisp::ReducedDispersion, seed_func, Nₚ, fname, ϵ₀::Float64 = 1e-10, max_iter = 30)
	print(max_iter, "\n")
	saver_o = Saver(fname, params, rdisp, Nₚ, max_iter)
	bs = seed_func(params)
	∂Fs, ∂²Fs = generate_shared_cash(params)
	psamples = generate_psamples(rdisp, Nₚ)
	ΔF = 0.0
	ϵ = 0.0
	t = @elapsed data = precompute_newton_step!(∂Fs, ∂²Fs, bs, params, psamples)
	print(t, " ", data[1], " ", ϵ, " ",  0, "\n")
	output(saver_o, bs, data..., 0.0, t, 0, params, rdisp, Nₚ, 0)
	for i in Base.OneTo(max_iter)
		t1 = @elapsed data1 = newton_step!(bs, ∂Fs, ∂²Fs, params, psamples)
		print(t1, " ", data1[1], " ", data1[4], " ", i, "\n")
		output(saver_o, bs, data1..., t1, i, params, rdisp, Nₚ, i)
		ΔF, ϵ = data1[1], data1[4]
		if ϵ < ϵ₀
			return bs, data1...
		end
	end
	return bs, ΔF, ∂Fs[1], ∂²Fs[1], ϵ
end

function compute_ΔFs(bs::AbstractVector{Float64}, params::AbstractParams, rdisp::ReducedDispersion, Nₚ::Int64)
	psamples = generate_psamples(rdisp, Nₚ)
	ΔFs_chunked = [SharedVector{Float64}(length(psamples[wid])) for wid in eachindex(psamples)]
    @sync begin
        for wid in true_workers()
			@spawnat wid process_ΔF_chunk!(ΔFs_chunked[wid], bs, params, psamples[wid])
        end
		if nworkers()>1
			@spawnat 1 process_ΔF_chunk!(ΔFs_chunked[1], bs, params, psamples[1])
		end
    end
	return vcat(ΔFs_chunked...)
end

function compute_ΔFs_sequential(bs::AbstractVector{Float64}, params::AbstractParams, rdisp::ReducedDispersion, Nₚ::Int64)
	psamples = get_psamples(rdisp, Nₚ)
	ΔFs = Vector{Float64}(undef, length(psamples))
    process_ΔF_chunk!(ΔFs, bs, params, psamples)
	return ΔFs
end

sum_free_energy(bs::AbstractVector{Float64}, ΔFs::AbstractVector{Float64}, params::AbstractParams) = final_part(bs, params) - sum(ΔFs)

function final_part(bs, params::AbstractParamsB1)
	N, m, β, u, u₁ = params.N, params.m, params.β, params.u, params.u₁
	mΔτ = β/N
	return mΔτ*(u*sum(abs2, bs[1:N]) + u₁*sum(abs2, bs[(N+1):2N]))
end

function final_part(bs, params::AbstractParamsNoB1)
	N, m, β, u = params.N, params.m, params.β, params.u
	mΔτ = β/N
	return mΔτ*u*sum(abs2, bs)
end

compute_free_energy(bs::AbstractVector{Float64}, params::AbstractParams, rdisp::ReducedDispersion, Nₚ::Int64) = sum_free_energy(bs, compute_ΔFs_sequential(bs, params, rdisp, Nₚ), params)


text:
"""
The `extrap_prep` function is used by `getindex_impl` to generate the body of
the `getindex` function for extrapolation objects.

The methods of `extrap_prep` work in "layers", iteratively working out the exact
expression needed.

The first layer takes a specification of the extrapolation scheme(s) to be used
and a `Val` object that specifies the dimensionality of the extrapolation object:
`extrap_prep{T,N}(::Type{T}, Val{N})`. These methods only dispatch to the second
layer, and need not be extended for new schemes.

The second layer also takes a `Val` object that specifies a single dimension on
which to work: `extrap_prep{T,N,d}(::Type{T}, ::Val{N}, ::Val{d}). The methods
with this signature in src/extrapolation/extrap_prep.jl simply expand into a
block with sub-expressions for handling too-low and too-high values separately
(the third layer), but specific interpolation schemes can provide more specific
methods for this layer that handle both ends simultaneously. For example, the
`Flat` scheme has a layer-2 method that uses `clamp` to restrict the coordinate
when used in both directions, but uses `min` and `max` when handling each end
separately.

The third layer, to which the second dispatches if no scheme-specific method is
found, adds a final `Val` object with a symbol `:lo` or `:hi`:
`extrap_prep{T,N,d,l}(::Type{T}, ::Val{N}, ::Val{d}, ::Val{l})`. These methods
must be specified for each extrapolation scheme. However, the general framework
takes care of expanding all possible tuple combinations, so individual schemes
need only care about e.g. `T==Flat`.

In addition to these methods, there is a similar three-layer method hierarchy
for gradient evaluation, in which a `Val{:gradient}` is prepended to the other
arguments:
`extrap_prep{T,N,d,l}(::Val{:gradient}`, ::Type{T}, ::Val{N}, ::Val{d}, ::Val{l})`
If nothing else is specified for the individual schemes, these methods forward
to the same methods without the `:gradient` argument, i.e. the same behavior as
for value extrapolation. This works well with all schemes that are simple
coordinate transformations, but for anything else methods for the low- and high-
value cases need to be implemented for each scheme.
""" extrap_prep

extrap_prep{T}(::Type{T}, n::Val{1}) = extrap_prep(T, n, Val{1}())
extrap_prep{T}(::Type{Tuple{T}}, n::Val{1}) = extrap_prep(T, n)
extrap_prep{T}(::Type{Tuple{T,T}}, n::Val{1}) = extrap_prep(T, n)
extrap_prep{T}(::Type{Tuple{Tuple{T,T}}}, n::Val{1}) = extrap_prep(T, n)
function extrap_prep{S,T}(::Type{Tuple{S,T}}, n::Val{1})
    quote
        $(extrap_prep(S, n, Val{1}(), Val{:lo}()))
        $(extrap_prep(T, n, Val{1}(), Val{:hi}()))
    end
end
extrap_prep{S,T}(::Type{Tuple{Tuple{S,T}}}, n::Val{1}) = extrap_prep(Tuple{S,T}, n)

# needed for ambiguity resolution
extrap_prep{T<:Tuple}(::Type{T}, ::Val{1}) = :(throw(ArgumentError("The 1-dimensional extrap configuration $T is not supported")))

function extrap_prep{T,N}(::Type{T}, n::Val{N})
    exprs = Expr[]
    for d in 1:N
        push!(exprs, extrap_prep(T, n, Val{d}()))
    end
    return Expr(:block, exprs...)
end
function extrap_prep{N,T<:Tuple}(::Type{T}, n::Val{N})
    length(T.parameters) == N || return :(throw(ArgumentError("The $N-dimensional extrap configuration $T is not supported - must be a tuple of length $N (was length $(lenght(T.parameters)))")))
    exprs = Expr[]
    for d in 1:N
        Tdim = T.parameters[d]
        if Tdim <: Tuple
            length(Tdim.parameters) == 2 || return :(throw(ArgumentError("The extrap configuration $Tdim for dimension $d is not supported - must be a tuple of length 2 or a simple configuration type")))
            if Tdim.parameters[1] != Tdim.parameters[2]
                push!(exprs, extrap_prep(Tdim, n, Val{d}()))
            else
                push!(exprs, extrap_prep(Tdim.parameters[1], n, Val{d}()))
            end
        else
            push!(exprs, extrap_prep(Tdim, n, Val{d}()))
        end
    end
    return Expr(:block, exprs...)
end
extrap_prep{T,N,d}(::Type{T}, n::Val{N}, dim::Val{d}) = extrap_prep(Tuple{T,T}, n, dim)
function extrap_prep{S,T,N,d}(::Type{Tuple{S,T}}, n::Val{N}, dim::Val{d})
    quote
        $(extrap_prep(S, n, dim, Val{:lo}()))
        $(extrap_prep(T, n, dim, Val{:hi}()))
    end
end


text:
# CheckParameters.jl

function checitemparameters(pars::IRTmodel)
    checkitem(pars)
end

function checkitem(pars::graded)
    p = take(pars)
    for p in take(pars)
        if isnan.(p) || isnan.(p)
            @error "Any parameter is $(typeof(p)). Estimation seems to fail."
            return true
        end
    end
    return false
end

function checkitem(pars::guessing)
    p = take(pars)
    for p in take(pars)
        if isnan.(p) || isnan.(p)
            @error "Any parameter is $p. Estimation seems to fail."
            return true
        end
    end
    if !(0.0 ≤ p[end] ≤ 1.0)
        @warn "Guessing parameter is out of range 0 to 1. Estimation seems to fail"
        return true
    end
    return false
end

function checkitem(pars::LPE)
    p = take(pars)
    for p in take(pars)
        if isnan.(p)
            @error "Any parameter is $(typeof(p)). Estimation seems to fail."
            @show pars
            return true
        end
    end
    return false
end


text:
using DataStructures

PTDFValues = SortedDict{String,SortedDict{String, Float64}}
PTDFDict = SortedDict{String, PTDFValues}

BASECASE = "BASECASE"

struct Network
    name::String
    # Direct access containers
    buses::SortedDict{String, Bus}
    branches::SortedDict{String, Branch}

    generators::SortedDict{String, Generator}

    #FIXME : create ptdf structure ? to indicate reference bus, distributed or not, cut branches, original file...
    # usecase => PTDF( Branch_id => Bus_id => value)
    ptdf::PTDFDict

    # Constructor
    function Network(name::String)
        return new( name
                  , SortedDict{String, Bus}()
                  , SortedDict{String, Branch}()
                  , SortedDict{String, Generator}()
                  , SortedDict{String,SortedDict{String, Float64}}()
                  )
    end
end

function Network()
    return Network("empty_network")
end

#########
## Bus ##
#########

function add_new_bus!(network::Network, bus_id::String)
    if haskey(network.buses,bus_id)
        throw( error("Bus with id ", bus_id, " already exists in Network ", network.name) )
    end
    bus = Bus(bus_id)
    network.buses[bus_id] = bus
end

function add_new_buses!(network::Network, buses_id::Vector{String})
    for bus_id in buses_id
        add_new_bus!(network, bus_id)
    end
end

function add_bus!(network::Network, bus::Bus)
    bus_id = bus.id
    if haskey(network.buses,bus_id)
        throw( error("Bus with id ", bus_id, " already exists in Network ", network.name) )
    end
    network.buses[bus_id] = bus
end

function add_buses!(network::Network, buses::Vector{Bus})
    for bus in buses
        add_bus!(network, bus)
    end
end

function get_bus(network::Network, bus_id::String)::Union{Missing, Bus}
    if haskey(network.buses,bus_id)
        return network.buses[bus_id]
    else
        return missing
    end
end

function safeget_bus(network::Network, bus_id::String)::Bus
    bus::Union{Bus, Missing} = get_bus(network, bus_id)
    if !isequal(bus, missing)
        return bus
    else
        throw( error("Bus with id ", bus_id, " does not exist in Network ", network.name) )
    end
end

function get_buses(network::Network)
    return collect(values(network.buses))
end

############
## Branch ##
############

function add_new_branch!(network::Network, id_branch::String, limit::Float64)
    try
        branch = Branch(id_branch, limit)
        network.branches[id_branch] = branch
    catch
        rethrow()
    end
end

function add_branch!(network::Network, branch::Branch)
    try
        network.branches[branch.id] = branch
    catch
        rethrow()
    end
end

function add_branches!(network::Network, branches::Vector{Branch})
    for branch in branches
        add_branch!(network, branch)
    end
end

function get_branch(network::Network, branch_id::String)::Union{Missing, Branch}
    # Get branch
    if haskey(network.branches, branch_id)
        return network.branches[branch_id]
    else
        return missing
    end
end

function safeget_branch(network::Network, branch_id::String)::Branch
    branch::Union{Missing, Branch} = get_branch(network, branch_id)
    if !isequal(branch, missing)
        return branch
    else
        throw( error("Branch ", branch_id, " does not exist in Network ", network.name) )
    end
end

function get_branches(network::Network, bus::Union{Int, Bus})::Vector{Tuple{Branch, Bool}}
    bus_id::Int = typeof(bus) == Int ? bus : bus.id
    branches = Vector{Tuple{Branch, Bool}}()
    if haskey(network.branches, bus_id)
        for branch in values(network.branches[bus_id]) # il y a surement moyen de faire plus efficace
            push!(branches, branch)
        end
    end
    return branches
end

function get_branches(network::Network)
    return values(network.branches)
end

########################
##      Generator     ##
########################

function add_new_generator_to_bus!(network::Network, bus::Union{Bus, String}, generator_id::String,
                                type::GeneratorType, p_min::Float64, p_max::Float64,
                                start_cost::Float64, prop_cost::Float64, dmo::Dates.Second, dp::Dates.Second)
    bus::Bus = typeof(bus) == String ? safeget_bus(network, bus) : bus
    if generator_id in map(bus->get_id(bus), get_buses(network))
        throw( error("Generators and buses must have distinct names :
                    a bus having the name ", generator_id, " already exists.") )
    end
    generator::Generator = add_new_generator!(bus, generator_id,
                                            type, p_min, p_max, start_cost, prop_cost, dmo, dp)
    # Store it in direct access containers
    network.generators[generator_id] = generator
end

function get_generator(network::Network, generator_id::String)
    if haskey(network.generators,generator_id)
        return network.generators[generator_id]
    else
        return missing
    end
end

function safeget_generator(network::Network, generator_id::String)
    generator::Union{Generator, Missing} = get_generator(network, generator_id)
    if !isequal(generator, missing)
        return generator
    else
        throw( error("Generator with id ", generator_id, " does not exist in Network ", network.name) )
    end
end

function get_generators(network::Network)
    return values(network.generators)
end

function get_generators_of_type(generators, type::GeneratorType)
    selected_generators = Vector{Generator}()
    for generator in generators
        if get_type(generator) == type
            push!(selected_generators, generator)
        end
    end
    return selected_generators
end
function get_generators_of_type(bus::Bus, type::GeneratorType)
    return get_generators_of_type(get_generators(bus), type)
end
function get_generators_of_type(network::Network, type::GeneratorType)
    return get_generators_of_type(get_generators(network), type)
end

########################
##        PTDF        ##
########################

function get_cases(network::Network)
    return keys(network.ptdf)
end

"""
# Arguments
    - `ptdf::PTDFDict` : PTDF data for different cases (N and N-1)
    - `case::String` : BASECASE for the basecase PTDF, name of the branch we cut in N-1 cases
"""
function safeget_ptdf(ptdf::PTDFDict, case::String=BASECASE)::PTDFValues
    if haskey(ptdf,case)
        return ptdf[case]
    else
        msg_l = @sprintf("Usecase %s is missing in PTDF data.", case)
        error(msg_l)
    end
end

function get_ptdf!(ptdf::PTDFDict, case::String=BASECASE)::PTDFValues
    if !haskey(ptdf,case)
        ptdf[case] = PTDFValues()
    end
    return ptdf[case]
end

function get_ptdf_component(ptdf_values::PTDFValues, branch_id::String)::Union{Missing,SortedDict{String, Float64}}
    if haskey(ptdf_values,branch_id)
        return ptdf_values[branch_id]
    else
        return missing
    end
end
function get_ptdf_component(ptdf::PTDFDict, branch_id::String, case::String=BASECASE)::Union{Missing,SortedDict{String, Float64}}
    ptdf_values = safeget_ptdf(ptdf, case)
    return get_ptdf_component(ptdf_values, branch_id)
end
function get_ptdf_component(network::Network, branch_id::String, case::String=BASECASE)::Union{Missing,SortedDict{String, Float64}}
    ptdf_values = safeget_ptdf(network.ptdf, case)
    return get_ptdf_component(ptdf_values, branch_id, case)
end

function safeget_ptdf_component(ptdf_values::PTDFValues, branch_id::String)::Union{Missing,SortedDict{String, Float64}}
    if haskey(ptdf_values,branch_id)
        return ptdf_values[branch_id]
    else
        throw( error("PTDF value for branch ", branch_id, " does not exist in PTDF") )
    end
end
function safeget_ptdf_component(ptdf::PTDFDict, branch_id::String, case::String=BASECASE)::Union{Missing,SortedDict{String, Float64}}
    ptdf_values = safeget_ptdf(ptdf, case)
    return safeget_ptdf_component(ptdf_values, branch_id)
end
function safeget_ptdf_component(network::Network, branch_id::String, case::String=BASECASE)::Union{Missing,SortedDict{String, Float64}}
    ptdf_values = safeget_ptdf(network.ptdf, case)
    return safeget_ptdf_component(ptdf_values, branch_id)
end

function get_ptdf_elt(ptdf_values::PTDFValues, branch_id::String, bus_id::String)::Union{Missing,Float64}
    if haskey(ptdf_values,branch_id)
        if haskey(ptdf_values[branch_id], bus_id)
            return ptdf_values[branch_id][bus_id]
        else
            return missing
        end
    else
        return missing
    end
end
function get_ptdf_elt(ptdf::PTDFDict, branch_id::String, bus_id::String, case::String=BASECASE)::Union{Missing,Float64}
    ptdf_values = safeget_ptdf(ptdf, case)
    return get_ptdf_elt(ptdf_values, branch_id, bus_id)
end
function get_ptdf_elt(network::Network, branch_id::String, bus_id::String, case::String=BASECASE)::Union{Missing,Float64}
    ptdf_values = safeget_ptdf(network.ptdf, case)
    return get_ptdf_elt(ptdf_values, branch_id, bus_id)
end

function safeget_ptdf_elt(ptdf_values::PTDFValues, branch_id::String, bus_id::String)::Float64
    if haskey(ptdf_values,branch_id)
        if haskey(ptdf_values[branch_id], bus_id)
            return ptdf_values[branch_id][bus_id]
        else
            throw( error("PTDF value for branch ", branch_id, " and bus ", bus_id, " does not exist.") )
        end
    else
        throw( error("PTDF value for branch ", branch_id, " does not exist.") )
    end
end
function safeget_ptdf_elt(ptdf::PTDFDict, branch_id::String, bus_id::String, case::String=BASECASE)::Float64
    ptdf_values = safeget_ptdf(ptdf, case)
    return safeget_ptdf_elt(ptdf_values, branch_id, bus_id)
end
function safeget_ptdf_elt(network::Network, branch_id::String, bus_id::String, case::String=BASECASE)::Float64
    ptdf_values = safeget_ptdf(network.ptdf, case)
    return safeget_ptdf_elt(ptdf_values, branch_id, bus_id)
end

function add_ptdf_elt!(ptdf_values::PTDFValues, branch_id::String, bus_id::String, ptdf_value::Float64)
    ptdf_component = get!(ptdf_values, branch_id, SortedDict{String, Float64}())
    if haskey(ptdf_component, bus_id)
        msg_l = @sprintf("PTDF already contains value for branch %s and bus %s", branch_id, bus_id)
        error(msg_l)
    else
        ptdf_component[bus_id] = ptdf_value
    end
end
function add_ptdf_elt!(ptdf::PTDFDict, branch_id::String, bus_id::String, ptdf_value::Float64, case::String=BASECASE)::Float64
    ptdf_values = get_ptdf!(ptdf, case)
    return add_ptdf_elt!(ptdf_values, branch_id, bus_id, ptdf_value)
end
function add_ptdf_elt!(network::Network, branch_id::String, bus_id::String, ptdf_value::Float64, case::String=BASECASE)::Float64
    ptdf_values = get_ptdf!(network.ptdf, case)
    return add_ptdf_elt!(ptdf_values, branch_id, bus_id, ptdf_value)
end

########################
##        utils       ##
########################

function get_generator_or_bus(network::Network, id::String)
    generator::Union{Generator, Missing} = get_generator(network, id)
    if !isequal(generator, missing)
        return generator
    else
        return get_bus(network, id) #bus or missing
    end
end

function safeget_generator_or_bus(network::Network, id::String)
    generator_or_bus::Union{Generator, Bus, Missing} = get_generator_or_bus(network, id)
    if ismissing(generator_or_bus)
        throw( error("No Generator or Bus with id `", id, "` exists in Network ", network.name) )
    else
        return generator_or_bus
    end
end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
# Definitions for converting between types.

Base.convert(::Type{CartesianTuple{N,S,U}}, t::CartesianTuple{N,T,U}) where {N,S,T,U} = CartesianTuple{N,S,U}(t.coordinates...)
Base.convert(::Type{CartesianTuple{N,S,U}}, t::CartesianTuple{N,T,V}) where {N,S,T, U<:VecLike, V<:VecLike} = CartesianTuple{N,S,U}(t.coordinates...)
Base.convert(::Type{Vect{N,S}}, p::Point{N,T}) where {N,S,T} = Vect{N,S}(p.coordinates...)
Base.convert(::Type{CartesianPair{S,U}}, p::Point3{T}) where {S,T,U} = CartesianPair{S,U}(p[1:2]...)
Base.convert(::Type{Point2{S}}, v::Vector2{T}) where {S,T} = Point2{S}(v.coordinates...)

# Promotion rules for Cartesian Tuples.

Base.promote_rule(::Type{CartesianTuple{N,S,U}}, ::Type{CartesianTuple{N,T,U}}) where {N,S,T,U} = CartesianTuple{N, promote_type(S,T), U}
Base.promote_rule(::Type{CartesianTuple{N,S,U}}, ::Type{T}) where {N,S,T,U} = CartesianTuple{N, promote_type(S,T), U}

# Convenience constructors.

for type ∈ (:Vector2, :Point2, :Vector3, :Point3, :Normal3)
    @eval begin
        $(type)(args::Vararg{T, size($(type){T})}) where T = $(type){T}(args...)
        $(type)(args::Vararg{Any, size($(type){T})}) where T = $(type)(promote(args...)...)
        $(Symbol("$(type)f"))(args::Vararg{T, size($(type){T})}) where T<:AbstractFloat = $(type){T}(args...)
        $(Symbol("$(type)f"))(args::Vararg{AbstractFloat, size($(type){T})}) where T = $(type)(promote(args...)...)
        $(Symbol("$(type)f"))(args::Vararg{Any, size($(type){T})}) where T = $(type)(convert.(AbstractFloat, args)...)
    end
end


text:
# The square root of 2 can be written as an infinite continued fraction.
#
# √2 = 1 + 1 / (2 + 1 / (2 + 1 / (2 + ...)))
#
# The infinite continued fraction can be written, √2 = [1;(2)], (2) indicates
# that 2 repeats ad infinitum. In a similar way, √23 = [4;(1,3,1,8)].
#
# It turns out that the sequence of partial values of continued fractions for
# square roots provide the best rational approximations. Let us consider the
# convergents for 2–√.
#
# 1 + 1 / 2 = 3 / 2
# 1 + 1 / (2 + 1 / 2) = 7 / 5
# 1 + 1 / (2 + 1 / (2 + 1 / 2)) = 17/12
# 1 + 1 / (2 + 1 / (2 + 1 / (2 + 1 / 2))) = 41/29
#
# Hence the sequence of the first ten convergents for √2 are:
#
# 1,3/2,7/5,17/12,41/29,99/70,239/169,577/408,1393/985,3363/2378,...
#
# What is most surprising is that the important mathematical constant,
# e = [2;1,2,1,1,4,1,1,6,1,...,1,2k,1,...].
#
# The first ten terms in the sequence of convergents for e are:
#
# 2,3,8/3,11/4,19/7,87/32,106/39,193/71,1264/465,1457/536,...
#
# The sum of digits in the numerator of the 10th convergent is 1+4+5+7=17.
#
# Find the sum of digits in the numerator of the 100th convergent of the
# continued fraction for e.

using ProjectEulerSolutions

# Simple if we use Julia's BigInts and rational built-ins.
function p065solution_rationals(n::Integer=10)::Integer

    coefs = Array{Integer,1}([2, 1, 2])
    for k in 2:div(n+2,3)
        append!(coefs, [1, 1, k*2])
    end

    start = BigInt(coefs[n])
    for i in n:-1:2
       start = coefs[i-1] + 1 // start
    end

    return sum(digits(start.num))
end

# Explicitly calculates the numerator and denominator.
function p065solution_explicit(n::Integer=10)::Integer

    coefs = Array{Integer,1}([2, 1, 2])
    for k in 2:div(n+2,3)
        append!(coefs, [1, 1, k*2])
    end

    num = BigInt(1)
    denom = BigInt(coefs[n])
    count = 0
    for i in n:-1:2
        num, denom = denom, num + coefs[i-1] * denom
    end

    return sum(digits(denom))
end

p065 = Problems.Problem(Dict("Rationals" => p065solution_rationals,
                             "Explicit" => p065solution_explicit))

Problems.benchmark(p065, 100)

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
using Test
using Pkg.Artifacts
using NCDatasets

using RRTMGP
using RRTMGP.Device: array_type, array_device
using RRTMGP.Vmrs
using RRTMGP.LookUpTables
using RRTMGP.AtmosphericStates
using RRTMGP.Optics
using RRTMGP.Sources
using RRTMGP.BCs
using RRTMGP.Fluxes
using RRTMGP.AngularDiscretizations
using RRTMGP.RTE
using RRTMGP.RTESolver

using CLIMAParameters
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()
# overriding some parameters to match with RRTMGP FORTRAN code
CLIMAParameters.Planet.grav(::EarthParameterSet) = 9.80665
CLIMAParameters.Planet.molmass_dryair(::EarthParameterSet) = 0.028964
CLIMAParameters.Planet.molmass_water(::EarthParameterSet) = 0.018016

include("reference_files.jl")
include("read_all_sky.jl")

DA = array_type()

function all_sky(
    ::Type{OPC},
    ::Type{FT},
    ::Type{I},
    ::Type{DA};
    use_lut::Bool = true,
) where {FT<:AbstractFloat,I<:Int,DA,OPC}
    opc = Symbol(OPC)
    FTA1D = DA{FT,1}
    FTA2D = DA{FT,2}
    max_threads = Int(256)
    n_gauss_angles = 1
    ncol = 128 # repeats col#1 128 time per RRTMGP example

    lw_file = get_ref_filename(:lookup_tables, :clearsky, λ = :lw)             # lw lookup tables for gas optics
    lw_cld_file = get_ref_filename(:lookup_tables, :cloudysky, λ = :lw)        # lw cloud lookup tables
    sw_file = get_ref_filename(:lookup_tables, :clearsky, λ = :sw)             # sw lookup tables for gas optics
    sw_cld_file = get_ref_filename(:lookup_tables, :cloudysky, λ = :sw)        # lw cloud lookup tables

    input_file = get_ref_filename(:atmos_state, :cloudysky)                    # all-sky atmos state
    if use_lut
        flux_file = get_ref_filename(:comparison, :cloudysky, lut_pade = :lut) # flux files for comparison (LUT)
    else
        flux_file = get_ref_filename(:comparison, :cloudysky, lut_pade = :pade) # flux files for comparison (pade)
    end

    #reading longwave gas optics lookup data
    ds_lw = Dataset(lw_file, "r")
    lookup_lw, idx_gases = LookUpLW(ds_lw, I, FT, DA)
    close(ds_lw)
    # reading longwave cloud lookup data
    ds_lw_cld = Dataset(lw_cld_file, "r")
    lookup_lw_cld = LookUpCld(ds_lw_cld, I, FT, DA, use_lut)
    close(ds_lw_cld)
    #reading shortwave gas optics lookup data
    ds_sw = Dataset(sw_file, "r")
    lookup_sw, idx_gases = LookUpSW(ds_sw, I, FT, DA)
    close(ds_sw)
    # reading longwave cloud lookup data
    ds_sw_cld = Dataset(sw_cld_file, "r")
    lookup_sw_cld = LookUpCld(ds_sw_cld, I, FT, DA, use_lut)
    close(ds_sw_cld)
    # reading input file
    ds_in = Dataset(input_file, "r")

    as, sfc_emis, sfc_alb_direct, sfc_alb_diffuse, zenith, toa_flux, bot_at_1 =
        setup_allsky_as(
            ds_in,
            idx_gases,
            lookup_lw,
            lookup_sw,
            lookup_lw_cld,
            lookup_sw_cld,
            use_lut,
            ncol,
            FT,
            DA,
            max_threads,
        )
    close(ds_in)
    ncol, nlay = as.ncol, as.nlay
    nlev = nlay + 1
    #---reading comparison files -----------------------------------
    if bot_at_1
        flip_ind = 1:nlev
    else
        flip_ind = nlev:-1:1
    end

    # Setting up longwave problem---------------------------------------
    ngpt_lw = lookup_lw.n_gpt
    op = OPC(FT, ncol, nlay, DA) # allocating optical properties object
    src_lw = source_func_longwave(FT, ncol, nlay, opc, DA)   # allocating longwave source function object
    bcs_lw = LwBCs{FT,typeof(sfc_emis),Nothing}(sfc_emis, nothing)    # setting up boundary conditions
    fluxb_lw = FluxLW(ncol, nlay, FT, DA)                             # flux storage for bandwise calculations
    flux_lw = FluxLW(ncol, nlay, FT, DA)                              # longwave fluxes

    # Setting up shortwave problem---------------------------------------
    ngpt_sw = lookup_sw.n_gpt
    src_sw = source_func_shortwave(FT, ncol, nlay, opc, DA)        # allocating longwave source function object

    # setting up boundary conditions
    inc_flux_diffuse = nothing

    bcs_sw = SwBCs(
        zenith,
        toa_flux,
        sfc_alb_direct,
        inc_flux_diffuse,
        sfc_alb_diffuse,
    )

    fluxb_sw = FluxSW(ncol, nlay, FT, DA) # flux storage for bandwise calculations
    flux_sw = FluxSW(ncol, nlay, FT, DA)  # shortwave fluxes for band calculations
    #-------------------------------------------------------------------
    # initializing RTE solver
    slv = Solver(
        as,
        op,
        src_lw,
        src_sw,
        bcs_lw,
        bcs_sw,
        fluxb_lw,
        fluxb_sw,
        flux_lw,
        flux_sw,
    )
    #------solvers
    solve_lw!(slv, max_threads, lookup_lw, lookup_lw_cld)
    solve_sw!(slv, max_threads, lookup_sw, lookup_sw_cld)

    #for i = 1:10
    #    @time solve_lw!(slv, max_threads, lookup_lw, lookup_lw_cld)
    #    @time solve_sw!(slv, max_threads, lookup_sw, lookup_sw_cld)
    #    println("------------------------------------------------")
    #end
    #-------------
    # comparison
    method = use_lut ? "Lookup Table Interpolation method" : "PADE method"
    ds_comp = Dataset(flux_file, "r")
    comp_flux_up_lw =
        Array(transpose(ds_comp["lw_flux_up"][:][1:ncol, flip_ind]))
    comp_flux_dn_lw =
        Array(transpose(ds_comp["lw_flux_dn"][:][1:ncol, flip_ind]))
    comp_flux_up_sw =
        Array(transpose(ds_comp["sw_flux_up"][:][1:ncol, flip_ind]))
    comp_flux_dn_sw =
        Array(transpose(ds_comp["sw_flux_dn"][:][1:ncol, flip_ind]))
    close(ds_comp)

    flux_up_lw = Array(slv.flux_lw.flux_up)
    flux_dn_lw = Array(slv.flux_lw.flux_dn)

    max_err_flux_up_lw = maximum(abs.(flux_up_lw .- comp_flux_up_lw))
    max_err_flux_dn_lw = maximum(abs.(flux_dn_lw .- comp_flux_dn_lw))

    println("=======================================")
    println("Cloudy-sky longwave test - $opc")
    println(method)
    println("max_err_flux_up_lw = $max_err_flux_up_lw")
    println("max_err_flux_dn_lw = $max_err_flux_dn_lw")

    flux_up_sw = Array(slv.flux_sw.flux_up)
    flux_dn_sw = Array(slv.flux_sw.flux_dn)
    flux_dn_dir_sw = Array(slv.flux_sw.flux_dn_dir)

    max_err_flux_up_sw = maximum(abs.(flux_up_sw .- comp_flux_up_sw))
    max_err_flux_dn_sw = maximum(abs.(flux_dn_sw .- comp_flux_dn_sw))

    println("Cloudy-sky shortwave test - $opc")
    println(method)
    println("max_err_flux_up_sw = $max_err_flux_up_sw")
    println("max_err_flux_dn_sw = $max_err_flux_dn_sw")
    println("=======================================")

    toler = FT(1e-5)

    @test max_err_flux_up_lw < toler
    @test max_err_flux_dn_lw < toler

    @test max_err_flux_up_sw < toler
    @test max_err_flux_dn_sw < toler
end
println("running cloudy lookup table version")
all_sky(TwoStream, Float64, Int, DA, use_lut = true)
println("running cloudy pade version")
all_sky(TwoStream, Float64, Int, DA, use_lut = false)


text:
#Convenience for extracting positioner mod and mon with error checks
function pos_mod_mon(coms, recs)
    pos = getpositioners(coms)
    if length(pos) > 1 error("Provide only one positioner command.") end
    pos_mon = getpositionermonitors(recs)
    if length(pos_mon) > 1 error("Provide only one positioner monitor.") end
    return pos[1], pos_mon[1]
end

#take the mean MON cycle and calculate the number of sampls that MON is shifted from MOD
function mon_lag_nsamps(pos, pos_mon, ncycs_ignore)
    pos_cyc_name = first(sequence_names(pos))
    pos_cycle = get_samples(pos, pos_cyc_name)
    @assert all(sequence_names(pos).==pos_cyc_name)
    nsamps_cycle = Base.length(pos_cycle) #should always be even
    mon_samps = get_samples(pos_mon)
    mon_samps = view(mon_samps, (ncycs_ignore*nsamps_cycle+1):Base.length(mon_samps))
    cycs_mon = ImagineAnalyses.get_cycles(mon_samps, nsamps_cycle)
    mean_cyc = mean(ustrip.(cycs_mon); dims=2)[:] * unit(mon_samps[1])
    return ImagineAnalyses.mon_delay(pos_cycle, mean_cyc)
end

#function local_minima(mms)
#    min_is = Int[]
#    if mms[2] > mms[1]
#        push!(mms, 1)
#    end
#    for i = 2:(length(mms)-1)
#        if mms[i-1] > mms[i] < mms[i+1]
#            push!(min_is, i)
#        end
#    end
#    if mms[end] < mms[end-1]
#        push!(min_is,length(mms))
#        warn("Found local minimum at the last temporal offset.  You should use a wider range of offsets")
#    end
#    return min_is
#end
#
##Instead of just taking the lowest mismatch, find all local minima and take the highest-indexed minimum that satisfies thresh_fac
#function ind_best_offset(mms; thresh_fac = 0.1)
#    min_is = local_minima(mms)
#    min_mm = minimum(mms)
#    #thresh = min_mm + (maximum(mms) - min_mm) * thresh_fac
#    thresh = min_mm + min_mm * thresh_fac
#    for i = length(min_is):-1:1
#        if mms[min_is[i]] <= thresh
#            return min_is[i]
#        end
#    end
#end

#imageseq is a 3D image array where each 2D slice is an image acquired for calibration
#imageseq should be either forward or reverse slice trials, but not both (use exclude_fwd to separate)
#The first (slow) stack should not be included in imgseq (use exclude_template_stack to get that one)
get_slice_trials(imgseq, sliceidx, ntrials, ntoffsets) = view(imgseq, :, :, (sliceidx-1)*(ntrials*ntoffsets)+1:sliceidx*(ntrials*ntoffsets))
exclude_template_stack(imgseq, nslices) = view(imgseq, :, :, (nslices+1):size(imgseq,3))
get_template_stack(imgseq, nslices) = imgseq[:,:,1:nslices]
#assumes template stack is excluded
get_fwd_imgs(imgseq) = view(imgseq, :, :, 1:2:size(imgseq,3))
get_bck_imgs(imgseq) = view(imgseq, :, :, 2:2:size(imgseq,3))

#NOTE: the timings returned here are additional offsets relative to any lag already applied with a Calibration when
#generating the command signals.

function slicetimings(img, toffsets, nslices::Int; allow_shifts=true, allow_rotations=false, subpixel=true)
    timings, mms, tfms = timings_mms_tfms(img, toffsets, nslices; allow_shifts=allow_shifts, allow_rotations=allow_rotations, subpixel=subpixel)
    return timings
end

function timings_mms_tfms(img, toffsets, nslices; kwargs...)
    @assert ndims(img) == 3 # should be stored by Imagine as a 2D timeseries
    #first set of images is always the target / template image
    target = get_template_stack(img, nslices)
    testimgs = exclude_template_stack(img, nslices)
    fwdimgs = get_fwd_imgs(testimgs)
    backimgs = get_bck_imgs(testimgs)
    timings_mms_tfms(target, fwdimgs, backimgs, toffsets, nslices; kwargs...)
end

function safe_fetch(rr)
    rslt = fetch(rr)
    if isa(rslt, RemoteException)
        throw(rslt)
    else
        return rslt
    end
end

#Setting record_all will return arrays-of-arrays of mismatches and transforms (one for each condition) rather than just the best
function timings_mms_tfms(target, fwdimgs, backimgs, toffsets, nslices; allow_shifts=true, subpixel=false, allow_rotations=false, record_all=false)
    ntoffsets =length(toffsets)
    ntrials = convert(Int, size(fwdimgs,3) / (ntoffsets*nslices))

    fwd_timings = eltype(toffsets)[]
    back_timings = eltype(toffsets)[]
    fwd_tfms = []; back_tfms = []; fwd_mms = []; back_mms = []

    i_start = 1
    i_stop = 0
    nw = nworkers()
    ws = workers()
    if isodd(nw) && nw > 1
        nw-=1
        ws = ws[1:end-1]
        warn("Using one less worker (must be even number)")
    end
    print("Computing with $(nw) worker processes\n")
    while i_stop < nslices
        i_start = i_stop+1
        i_stop = min(nslices, i_stop+max(1, div(nw,2)))
        idxs = i_start:i_stop
        fwdrrs = []; bckrrs = [];
        for i = 1:length(idxs)
            widx_odd = 2*i-1
            fixed = view(target, :, :, idxs[i])
            movfwd = get_slice_trials(fwdimgs, idxs[i], ntrials, ntoffsets)
            #launch on odd-numbered worker
            push!(fwdrrs, remotecall(toffset_fits, ws[widx_odd], fixed, movfwd, toffsets; allow_shifts=allow_shifts, allow_rotations=allow_rotations, subpixel=subpixel))
            movback = get_slice_trials(backimgs, idxs[i], ntrials, ntoffsets)
            #launch on even-numbered worker if we have more than one worker
            if nw > 1
                push!(bckrrs, remotecall(toffset_fits, ws[widx_odd+1], fixed, movback, toffsets; allow_shifts=allow_shifts, allow_rotations=allow_rotations, subpixel=subpixel))
            else
                push!(bckrrs, remotecall(toffset_fits, ws[widx_odd], fixed, movback, toffsets; allow_shifts=allow_shifts, allow_rotations=allow_rotations, subpixel=subpixel))
            end
        end
        #fetch them all and store important results
        for i = 1:length(idxs)
            ftfms, fmms = safe_fetch(fwdrrs[i])
            #fwdi = ind_best_offset(fmms)
            fwdi = argmin(fmms)
            @show fwd_toffset = toffsets[fwdi]
            @show fwd_mm = fmms[fwdi]
            push!(fwd_timings, fwd_toffset)

            btfms, bmms = safe_fetch(bckrrs[i])
            #backi = ind_best_offset(bmms)
            backi = argmin(bmms)
            @show back_toffset = toffsets[backi]
            @show back_mm = bmms[backi]
            push!(back_timings, back_toffset)

            if record_all
                push!(fwd_tfms, ftfms)
                push!(back_tfms, btfms)
                push!(fwd_mms, fmms)
                push!(back_mms, bmms)
            else
                push!(fwd_tfms, ftfms[fwdi])
                push!(back_tfms, btfms[backi])
                push!(fwd_mms, fwd_mm)
                push!(back_mms, back_mm)
            end
            #print("Slice $(idxs[i]) complete\n")
        end
        print("------------------Finished $i_stop of $nslices slices-----------------------\n")
    end
    return (fwd_timings, back_timings), (fwd_mms, back_mms), (fwd_tfms, back_tfms)
end

function preprocess(img::AbstractArray{Float64,2}; sigmas=(1.0,1.0,), sqrt_tfm=true, bias=Float64(reinterpret(N0f16, UInt16(100))))
    img = imfilter(img, KernelFactors.IIRGaussian(Float64, sigmas))
    output = similar(img)
    for i in eachindex(img)
        temp = max(0.0, img[i] - bias)
        output[i] = ifelse(sqrt_tfm, sqrt(temp), temp)
    end
    return output
end

align2d(fixed, moving; kwargs...) = align2d(fixed, Float64.(moving); kwargs...)
function align2d(fixed, moving::AbstractMatrix{Float64}; thresh_fac=0.9, sigmas=(1.0,1.0,), allow_shifts = true, allow_rotations =false, subpixel=false)
    mxshift = (8,8,)
    moving = preprocess(moving; sigmas=sigmas, sqrt_tfm=false)
    thresh = (1-thresh_fac) * sum(abs2.(fixed[.!(isnan.(fixed))]))
    if allow_rotations
        maxradians = pi/180
        tfm, mm = RegisterQD.qd_rigid(fixed, moving, mxshift, [maxradians], [pi/10000]; thresh=thresh, initial_tfm=IdentityTransformation())
        return tfm, mm
    elseif allow_shifts
        if !subpixel
            shft, mm = RegisterQD.best_shift(fixed, moving, mxshift, thresh; normalization=:intensity, initial_tfm=IdentityTransformation())
            return Translation(shft...), mm
        else
            return RegisterQD.qd_translate(fixed, moving, mxshift; thresh=thresh, rtol=0.0, fvalue=0.0, maxevals=200)
        end
    else
        mm = RegisterMismatch.mismatch0(fixed, moving; normalization=:intensity)
        return Translation(0.0,0.0), ratio(mm, 0.0)
    end
end

#get_offset_trials returns the set of trials at a particular temporal offset
#Assumes that img_seq:
#   -does not include the slow "target" stack
#   -includes only forward or reverse images, not both
#   -includes data for only one slice
#(This means img_seq has size ntrials * ntoffsets in dimension 3)
get_toffset_trials(img_seq, itoffset, ntrials) = view(img_seq, :, :, (itoffset-1)*ntrials+1:itoffset*ntrials)

#consistent per-slice.  It depends on so many factors that we just do it emprically.
function toffset_fits(target, testimgs, toffsets::Vector; sigmas=(1.0,1.0), allow_shifts=true, subpixel=false, allow_rotations = false)
    ntoffsets = length(toffsets)::Int
    ntrials = convert(Int, size(testimgs,3) / ntoffsets)
    tfms = []
    mms = Float64[]
    target = preprocess(Float64.(target); sigmas=sigmas, sqrt_tfm=false)
    #Threads.@threads for i = 1:length(toffsets)  #Currently segfaults on julia 0.6.2
    for i = 1:length(toffsets)
        moving = dropdims(mean(Float64.(get_toffset_trials(testimgs, i, ntrials)); dims=3); dims=3)
        tfm, mm = align2d(target, moving; sigmas=sigmas, allow_shifts=allow_shifts, subpixel=subpixel, allow_rotations=allow_rotations)
        #@show mm
        push!(tfms, tfm)
        push!(mms, mm)
    end
    return tfms, mms
end

function find_bidi_centers(cyc::AbstractVector, sampr::HasInverseTimeUnits)
    @show pmin = minimum(cyc)
    @show pmax = maximum(cyc)
    pctr = pmin + (pmax-pmin)/2 #we will place the flash in the center
    find_bidi_locations(cyc, sampr, pctr)
end

function find_bidi_locations(cyc::AbstractVector, sampr::HasInverseTimeUnits, ploc::HasLengthUnits)
    stack_rate = 2 * (sampr / length(cyc))
    pad_secs = (1/stack_rate)/200 #wait 0.5% of a cycle before accepting threshold crossings (helps with noise)
    pad_nsamps = round(Int, pad_secs * sampr)
    mon_is = ImagineAnalyses.find_circular(cyc, [ploc;], pad_nsamps)[1]
    @assert Base.length(mon_is) == 2 #should have just two crossings, one fwd one back
    fwd_i = mon_is[1]
    back_i = mon_is[2]
    return fwd_i, back_i
end


text:
using FFTW
using LaTeXStrings
using Printf
using PyPlot

function dsdemo2()
    println("Discrete-Time Simulation")

    fig = figure(1)
    fig.clf()
    fig.canvas.set_window_title("Modulator Input & Output")

    OSR = 32
    H = synthesizeNTF(5, OSR, 1)
    N = 8192
    fB = ceil(Int, N/(2*OSR));
    ftest = floor(Int, 2/3*fB);
    u = 0.5 * sin.(2*pi*ftest/N * (0:N-1))
    v, = simulateDSM(u, H)
    t = 0:100

    step(t, u[t.+1], color="r", linewidth=1)
    step(t, v[t.+1], color="g", linewidth=1)
    axis([0, 100, -1.2, 1.2])
    xlabel("Sample Number")
    ylabel("u, v")

    print("paused")
    readline()

    fig = figure(2)
    fig.clf()
    fig.canvas.set_window_title("Output Spectrum")

    f = range(0, 0.5; length=N÷2+1)
    spec = fft(v .* ds_hann(N) / (N/4))

    plot(f, dbv(spec[1:N÷2+1]), color="b", linewidth=1)
    figureMagic([0, 0.5], 0.05, 2, [-120, 0], 20, 1)
    xlabel("Normalized Frequency")
    ylabel("dBFS")

    print("paused")
    readline()

    snr = calculateSNR(spec[3:fB+1], ftest-2)

    text(0.05, -10, @sprintf("SNR = %4.1fdB @ OSR = %d", snr, OSR))

    print("paused")
    readline()

    NBW = 1.5 / N
    Sqq = 4 * evalTF(H, exp.(2*im*pi*f)).^2 / 3

    plot(f, dbp(Sqq*NBW), color="m", linewidth=2)
    text(0.5, -90, string(@sprintf("NBW = %4.1E × ", NBW), L"f_s"), ha="right")
    legend(["Simulation", "Expected PSD"])

    print("paused")
    readline()

    fig = figure(3)
    fig.clf()
    fig.canvas.set_window_title("SQNR")

    snr_pred, amp_pred = predictSNR(H, OSR)
    snr, amp = simulateSNR(H, OSR)
    pk_snr, pk_amp = peakSNR(snr, amp)

    plot(amp_pred, snr_pred, linewidth=1)
    plot(amp, snr, color="g", linestyle="None", marker="o", markerfacecolor="none", markeredgewidth=0.5)
    text(-25, 85, @sprintf("peak SNR = %4.1fdB\n@ OSR = %d\n", pk_snr, OSR), ha="right")
    figureMagic([-100, 0], 10, 1, [0, 100], 10, 1)
    xlabel("Input Level (dBFS)")
    ylabel("SQNR (dB)")

    println("paused")
    readline()

    println("Bandpass Modulator")

    fig = figure(1)
    fig.clf()

    f0 = 1/8
    OSR = 64
    H = synthesizeNTF(8, OSR, 1, 1.5, f0)
    fB = ceil(Int, N/(2*OSR))
    ftest = round(Int, f0*N + 1/3*fB)
    u = 0.5 * sin.(2*pi*ftest/N * (0:N-1))   # half-scale sine-wave input
    v, = simulateDSM(u, H)
    t = 0:100

    step(t, u[t.+1], color="r", linewidth=1)
    step(t, v[t.+1], color="g", linewidth=1)
    axis([0, 100, -1.2, 1.2])
    xlabel("Sample Number")
    ylabel("u, v")

    fig.canvas.manager.window.raise_()

    print("paused")
    readline()

    fig = figure(2)
    fig.clf()

    f = range(0, 0.5; length=N÷2+1)
    spec = fft(v .* ds_hann(N) / (N/4))
    f1 = round(Int, (f0-0.25/OSR)*N)
    f2 = round(Int, (f0+0.25/OSR)*N)
    snr = calculateSNR(spec[f1:f2], ftest-f1+1)

    plot(f, dbv(spec[1:N÷2+1]), color="b", linewidth=1)
    text(0.15, -10, @sprintf("SNR = %4.1fdB @ OSR = %d", snr, OSR))
    figureMagic([0, 0.5], 0.05, 2, [-140, 0], 20, 1)
    xlabel("Normalized Frequency")
    ylabel("dBFS")

    fig.canvas.manager.window.raise_()

    print("paused")
    readline()

    NBW = 1.5/N
    Sqq = 4 * evalTF(H, exp.(2*im*pi*f)).^2 / 3

    plot(f, dbp(Sqq*NBW), color="m", linewidth=2)
    text(0.5, -90, string(@sprintf("NBW = %4.1E × ", NBW), L"f_s"), ha="right")
    legend(["Simulation", "Expected PSD"])

    print("paused")
    readline()

    fig = figure(3)
    fig.clf()

    snr_pred, amp_pred = predictSNR(H, OSR, [-120:10:-20; -15; -10:0], f0)
    snr, amp = simulateSNR(H, OSR, [-120:10:-20; -15; -10:0], f0)
    pk_snr, pk_amp = peakSNR(snr, amp)

    plot(amp_pred, snr_pred, linewidth=1)
    plot(amp, snr, color="g", linestyle="None", marker="o", markerfacecolor="none", markeredgewidth=0.5)
    text(-20, 95, @sprintf("peak SNR = %4.1fdB\n@ OSR = %d\n", pk_snr, OSR), ha="right")
    figureMagic([-110, 0], 10, 1, [0, 110], 10, 1)
    xlabel("Input Level (dBFS)")
    ylabel("SQNR (dB)")

    fig.canvas.manager.window.raise_()

    println("paused")
    readline()

    fig = figure(3)
    fig.clf()
    fig = figure(2)
    fig.clf()
    fig = figure(1)
    fig.clear()
    fig.canvas.set_window_title("Input & Output")
    fig.subplots_adjust(hspace=0.5)

    colors = ["m", "b"]
    Hinf_list = [2, 8]
    for i = 1:2
        Hinf = Hinf_list[i]
        col = colors[i]
        println("15-step 7th-order Lowpass")
        println("Hinf = $Hinf")

        fig = figure(1)

        OSR = 8
        M = 16
        H = synthesizeNTF(7, OSR, 1, Hinf)
        N = 8192
        fB = ceil(Int, N/(2*OSR))
        ftest = floor(Int, 2/7*fB)
        u = 0.5*M*sin.(2*pi*ftest/N*(0:N-1))
        v, = simulateDSM(u, H, M+1)
        t = 0:100

        subplot(2, 1, i)

        step(t, u[t.+1], color="g", linewidth=1)
        step(t, v[t.+1], color=col, linewidth=1)
        figureMagic([0, 100], 10, 2, [-M, M], 2, 4)
        xlabel("Sample Number")
        ylabel("u, v")

        fig.canvas.manager.window.raise_()

        print("paused")
        readline()

        fig = figure(2)

        f = range(0, 0.5; length=N÷2+1)
        spec = fft(v .* ds_hann(N) / (M*N/4))
        snr = calculateSNR(spec[3:fB+1], ftest-2)
        NBW = 1.5/N
        Sqq = 4 * evalTF(H, exp.(2*im*pi*f)).^2 / (3*M^2)

        plot(f, dbv(spec[1:N÷2+1]), color=col)
        plot(f, dbp(Sqq*NBW), color="c", linewidth=2)
        text(0.1, 10*(i-3), @sprintf("SNR = %4.1fdB @ OSR = %d", snr, OSR), color=col)
        if i == 1
            text(0.5, -110, string(@sprintf("NBW = %4.1E × ", NBW), L"f_s"), ha="right")
        end
        figureMagic([0, 0.5], 0.05, 2, [-160, 0], 20, 1)
        xlabel("Normalized Frequency")
        ylabel("dBFS")

        fig.canvas.manager.window.raise_()

        println("paused")
        readline()

        fig = figure(3)

        snr, amp = simulateSNR(H, OSR, vcat(-120:10:-20, -15, -10:0), 0, M+1)
        pk_snr, pk_amp = peakSNR(snr, amp)

        plot(amp, snr, color=col, linestyle="None", marker="o", markerfacecolor="none", markeredgewidth=0.5)
        plot(amp, snr, color=col, linestyle="--", linewidth=1)
        text(-13, pk_snr, @sprintf("peak SNR = %4.1fdB\n@ OSR = %d\n", pk_snr, OSR), fontsize=8, color=col, ha="right", va="top")
        figureMagic([-120, 0], 10, 2, [0, 120], 10, 2)
        xlabel("Input Level (dBFS)")
        ylabel("SQNR (dB)")

        fig.canvas.manager.window.raise_()
    end

    return nothing
end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
using ElasticStabilityConditions
using Test

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


text:
using BlackBoxOptim

# We also need the (not yet export) num_func_evals
using BlackBoxOptim: num_func_evals

function rosenbrock2d(x)
    return (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
end

fitness_progress_history = Array{Tuple{Int, Float64},1}()
callback = oc -> push!(fitness_progress_history, (num_func_evals(oc), best_fitness(oc)))

# A callback interval of 0.0 ensures we callback every step of the optimization.
res = bboptimize(rosenbrock2d;
        SearchRange=(-5.0,5.0), NumDimensions = 2,
        PopulationSize=10, MaxFuncEvals=200,
        CallbackFunction = callback, CallbackInterval = 0.0);

# Note that fitness_progress_history is NOT guaranteed to have all values for the
# number of function evaluations since it depends on the optimization algorithm how
# many function evaluations they do in each step. The callback happen only between
# steps. The number of func evals might even be stochastic or change per step.
println(fitness_progress_history)


text:
include("unrealizedNormals.jl")

meansteps = 10000
tol = .2 #kinda loose, but not too loose; ideally this would be calculated analytically but this'll do

function approxEq(a, b, tol = tol)
    a-tol < b < a+tol
end

function compareModels(integratedModel, explicitModel, args, constraints, locus, meansteps = meansteps)
    tr, w = generate(integratedModel, args, constraints)
    tr2s, logws, w2 = importance_sampling(explicitModel, args, constraints, meansteps)
    ws = AnalyticWeights(exp.(logws[1:meansteps]))
    locusvals = [tr2s[i][locus] for i in 1:meansteps]
    m = mean(locusvals, ws)
    vv = var(locusvals, ws, corrected=true)
    @info "Comparing models:" locus w w2 m tr[locus].currentMean vv tr[locus].currentCondSD^2
    @assert approxEq(w, w2) "Estimated weights unequal"
    @assert approxEq(m, tr[locus].currentMean) "Estimated means at locus unequal"
    @assert approxEq(vv, tr[locus].currentCondSD^2) "Estimated variances at locus unequal"
end

mul2 = 2.

@gen function testItExplicit()

    w ~ normal(0., 1.)
    x ~ normal(0., 1.)
    y ~ normal(x+w, 2.)
    z ~ normal(x-w, 2.)
    y2 ~ normal(x+w*mul2, 2.)
    z2 ~ normal(x-w*mul2, 2.)

end

@gen function testItIntegrated()

    w ~ unrealizedIndyNormal(0., 1.)
    x ~ unrealizedIndyNormal(0., 1.)
    y ~ observedNormal(x+w, 2.)
    z ~ observedNormal(x-w, 2.)
    y2 ~ observedNormal(x+w*mul2, 2.)
    z2 ~ observedNormal(x-w*mul2, 2.)
    #Note: must manually ensure that w and x are independent conditional on y, z, y2, z2.
end

if @isdefined doTestUnrealizedNormals
    for v in 0.:.1:1.
        c = Gen.choicemap()
        c[:y] = v
        c[:z] = v
        c[:y2] = v
        c[:z2] = v
        compareModels(testItIntegrated, testItExplicit, (), c, :w)
    end
end

max_stars_repo_name: jamesonquinn/Gen

max_issues_repo_name: jamesonquinn/Gen

max_forks_repo_name: jamesonquinn/Gen


text:
function apply_laser_over_atoms(laser::Laser{Gaussian3D}, atoms)
    @debug "start: apply_laser_over_atoms - Laser{Gaussian3D}"

    w₀, s, Δ = laser.pump.w₀, laser.s, laser.Δ

    E = zeros(Complex{eltype(atoms.r)}, atoms.N)
    E₀ = estimate_E₀(laser)
    v_atoms = view(Float64.(atoms.r), :, :)

    Threads.@threads for n in 1:(atoms.N)
        oneAtom = v_atoms[:,n]
        @inbounds E[n] = _core_Gaussian3D(oneAtom, E₀, w₀)
    end

    @debug "end  : apply_laser_over_atoms - Laser{Gaussian3D}"
    return E
end

function apply_laser_over_oneAtom(laser::Laser{Gaussian3D}, position::AbstractArray)
    w₀ = laser.pump.w₀

    E₀ = estimate_E₀(laser)
    E = _core_Gaussian3D(position, E₀, w₀)

    return E
end

function _core_Gaussian3D(oneAtom, E₀, w₀)
    x, y, z = oneAtom[1], oneAtom[2], oneAtom[3]

    ## This formula is stable for z==0
    # Ref: Eq 3.11 from "CHAPTER 3. PROPAGATION AND FOCUSING OF OPTICAL FIELDS"
    denominator_factor = 1 .+ 2im .* z / (k₀ * w₀^2)
    Eᵢ = E₀ * cis(+k₀ * z)
    Eᵢ = Eᵢ .* exp.(-(x .^ 2 + y .^ 2) ./ (denominator_factor .* w₀^2))
    Eᵢ = Eᵢ ./ denominator_factor

    return Eᵢ
end


text:
"""
    struct DuffyQuadrature{D,T} <: Quadrature{D,T}
      coordinates::Vector{Point{D,T}}
      weights::Vector{T}
    end

Duffy quadrature for simplices in [0,1]^D
"""
struct DuffyQuadrature{D,T} <: Quadrature{D,T}
  coordinates::Vector{Point{D,T}}
  weights::Vector{T}
end

get_coordinates(q::DuffyQuadrature) = q.coordinates

get_weights(q::DuffyQuadrature) = q.weights

"""
    DuffyQuadrature{D}(degree::Integer) where D
"""
function DuffyQuadrature{D}(order::Integer) where D
  x,w = _duffy_quad_data(order,D)
  DuffyQuadrature(x,w)
end

function _duffy_quad_data(order::Integer,D::Int)

  beta = 0
  dim_to_quad_1d = [
    _gauss_jacobi_in_0_to_1(order,(D-1)-(d-1),beta) for d in 1:(D-1) ]

  quad_1d = _gauss_legendre_in_0_to_1(order)
  push!(dim_to_quad_1d,quad_1d)

  x_pos = 1
  w_pos = 2
  dim_to_xs_1d = [quad_1d[x_pos] for quad_1d in dim_to_quad_1d]
  dim_to_ws_1d = [quad_1d[w_pos] for quad_1d in dim_to_quad_1d]

  a = 0.5
  for d in (D-1):-1:1
    ws_1d = dim_to_ws_1d[d]
    ws_1d[:] *= a
    a *= 0.5
  end

  x,w = _tensor_product_duffy(dim_to_xs_1d,dim_to_ws_1d)

  (_duffy_map.(x),w)

end

# Duffy map from the n-cube in [0,1]^d to the n-simplex in [0,1]^d
function _duffy_map(q::Point{D,T}) where {D,T}
  m = zero(mutable(Point{D,T}))
  m[1] = q[1]
  a = one(T)
  for i in 2:D
    a *= (1-q[i-1])
    m[i] = a*q[i]
  end
  Point{D,T}(m)
end

_duffy_map(q::Point{1,T}) where T = q

function _gauss_jacobi_in_0_to_1(order,alpha,beta)
  n = _npoints_from_order(order)
  x,w = gaussjacobi(n,alpha,beta)
  _map_to(0,1,x,w)
end

function _gauss_legendre_in_0_to_1(order)
  n = _npoints_from_order(order)
  x,w = gausslegendre(n)
  _map_to(0,1,x,w)
end

# Transforms a 1-D quadrature from `[-1,1]` to `[a,b]`, with `a<b`.
function _map_to(a,b,points,weights)
  points_ab = 0.5*(b-a)*points .+ 0.5*(a+b)
  weights_ab = 0.5*(b-a)*weights
  (points_ab, weights_ab)
end

function _npoints_from_order(order)
  ceil(Int, (order + 1.0) / 2.0 )
end

function _tensor_product_duffy(
  dim_to_xs_1d::Vector{Vector{T}},
  dim_to_ws_1d::Vector{Vector{W}}) where {T,W}

  D = length(dim_to_ws_1d)
  @assert D == length(dim_to_xs_1d)
  dim_to_n = [length(ws_1d) for ws_1d in dim_to_ws_1d]
  n = prod(dim_to_n)
  xs = zeros(Point{D,T},n)
  ws = zeros(W,n)
  cis = CartesianIndices(tuple(dim_to_n...))
  m = zero(mutable(Point{D,T}))
  _tensor_product_duffy!(xs,ws,dim_to_xs_1d,dim_to_ws_1d,cis,m)
  (xs,ws)
end

function _tensor_product_duffy!(
  xs,ws,dim_to_xs_1d,dim_to_ws_1d,cis::CartesianIndices{D},m) where D
  k = 1
  for ci in cis
    w = 1.0
    for d in 1:D
      xs_1d = dim_to_xs_1d[d]
      ws_1d = dim_to_ws_1d[d]
      i = ci[d]
      xi = xs_1d[i]
      wi = ws_1d[i]
      w *= wi
      m[d] = xi
    end
    xs[k] = m
    ws[k] = w
    k += 1
  end
end