text:
"""
    gcve_pdf(epochlength::Real,fs::Int;filter=nothing,itr=1000000,pruned=0.5, σ = 1)

Calculate the distribution of the coefficient of variation of the envelope (CVE) of discrete-time zero-mean white Gaussian noise using random samples from N(0,1).

# Arguments
* `epochlength`: The length of the epoch in seconds
* `fs`: The sampling frequency
* `filter`: A filter function
* `itr`: The number of epochs used to compute the distribution
* `pruned`: The number of seconds cut from both tails of the envelope before calculating the CVE
* `σ`: Standard deviation of the Normal distribution

# Notes
* Notice that the pdf for the CVE of white Gaussian noise can also be bootstrapped using Fourier transform phase randomization.

# Examples
```julia
julia> gcve_pdf(10,1000,itr=10)
10-element Array{Float64,1}:
 0.5246024892582911
 0.52295322056966
 0.5256115270351591
 0.5263961505330913
 ⋮
 0.5237027693791765
 0.5181673343579238
 0.5260079659619668
```
"""
function gcve_pdf(epochlength::Real,fs::Int; filter = nothing,itr = 1000000,pruned = 0.5,σ = 1)
    epochlength>2*pruned ? nothing : error("Epoch length too short")
	isinteger(epochlength*fs) ? nothing : error("Number of points must be integer")
	isinteger(pruned*fs) ? nothing : error("Prunning points must be integer")
    println("The effective epoch length after prunning is $(epochlength-2*pruned) seconds")

	points = Int(epochlength*fs)
	d = Normal(0,σ)
    gcve_expdistro = zeros(Float64,itr)

    if filter != nothing
    	println("\n","Calculating CVE pdf for filtered white Gaussian noise...")
        @inbounds for i in 1:itr
            println("Iteration $(i)/$(itr)")
            gaussian_sample = rand(d, points)
            gaussian_sample = filter(gaussian_sample,fs)
            gaussian_sample = abs.(DSP.hilbert(gaussian_sample))
            gaussian_sample = gaussian_sample[Int(fs*pruned)+1:end-Int(fs*pruned)]
            gcve_expdistro[i] = std(gaussian_sample)/mean(gaussian_sample)
        end
    else
	println("\n","Calculating CVE pdf for non-filtered white Gaussian noise...")
        @inbounds for i in 1:itr
            println("Iteration $(i)/$(itr)")
            gaussian_sample = rand(d, epochlength*fs)
            gaussian_sample = abs.(DSP.hilbert(gaussian_sample))
            gaussian_sample = gaussian_sample[Int(fs*pruned)+1:end-Int(fs*pruned)]
            gcve_expdistro[i] = std(gaussian_sample)/mean(gaussian_sample)
        end
    end
    return gcve_expdistro
end

"""
    rcv_pdf(points::Int; itr=1000000, σ = 1)

Calculate the distribution of the coefficient of variation of Rayleigh noise using random samples from Rayleigh(σ).

# Arguments
* `epochlength`: The length of the epoch in seconds
* `fs`: The sampling frequency
* `itr`: The number of epochs used to compute the distribution
* `σ`: Scale parameter of the Rayleigh distribution

# Examples
```julia
julia> rcv_pdf(10000,itr=10,σ = 2)
10-element Array{Float64,1}:
 0.5208939217095829
 0.523914916275298
 0.5242368643462019
 0.519126955367746
 ⋮
 0.5229001421212949
 0.5239581722102657
 0.5275819378236367
```
"""
function rcv_pdf(points::Int;itr=1000000, σ = 1)

	d = Rayleigh(σ)
    rcv_expdistro = zeros(Float64,itr)

    println("\n","Calculating CV pdf for Rayleigh noise...")
    @inbounds for i in 1:itr
            println("Iteration $(i)/$(itr)")
            rayleigh_sample = rand(d, points)
            rcv_expdistro[i] = std(rayleigh_sample)/mean(rayleigh_sample)
    end
    return rcv_expdistro
end

"""
    ucve_pdf(epochlength::Int,fs::Int;filter=nothing,itr=1000000,prunning=1)

Calculate the distribution of the coefficient of variation of the envelope (CVE) of discrete-time zero-mean uniform noise.

# Arguments
* `epochlength`: The length of the epoch in seconds
* `fs`: The sampling frequency
* `filter`: A filter function
* `itr`: The number of epochs used to compute the distribution
* `prunning`: The number of seconds cut from both tails of the envelope before calculating the CVE

# Notes
* Notice that the pdf for the CVE of white Gaussian noise can also be bootstrapped using Fourier transform phase randomization.

# Examples
```julia
julia> ucve_pdf(10,1000,itr=10)
10-element Array{Float64,1}:
 0.42808995848724163
 0.4353783310060235
 0.4341374928673676
 0.43346989277220405
 ⋮
 0.4317732242362664
 0.42726595981167065
 0.4327565016764864
```
"""
function ucve_pdf(epochlength::Int,fs::Int;filter=nothing,itr=1000000,prunning=1)
    epochlength>2*prunning ? nothing : error(" Epoch length too short")
    println("The effective epoch length after prunning is $(epochlength-prunning) seconds")
    d = Uniform(-1,1)
    ucve_expdistro = zeros(Float64,itr)

    if filter != nothing
    	println("\n","Calculating CVE pdf for filtered uniform noise...")
        @inbounds for i in 1:itr
            println("Iteration $(i)/$(itr)")
            uniform_sample = rand(d, epochlength*fs)
            uniform_sample = filter(uniform_sample,fs)
            uniform_sample = abs.(DSP.hilbert(uniform_sample))
            uniform_sample = uniform_sample[(fs*prunning)+1:end-(fs*prunning)]
            ucve_expdistro[i] = std(uniform_sample)/mean(uniform_sample)
        end
    else
	println("\n","Calculating CVE pdf for non-filtered uniform noise...")
        @inbounds for i in 1:itr
            println("Iteration $(i)/$(itr)")
			uniform_sample = rand(d, epochlength*fs)
            uniform_sample = abs.(DSP.hilbert(uniform_sample))
            uniform_sample = uniform_sample[(fs*prunning)+1:end-(fs*prunning)]
            ucve_expdistro[i] = std(uniform_sample)/mean(uniform_sample)
        end
    end
    return ucve_expdistro
end


text:
adv_eq_dir  = joinpath(output_root,"AdvectionEquation")

print_norms = false

@testset "Linear advection, ∂_t ϕ + ∇•(cϕ) = 0 ∈ ∂Ω, ϕ(t=0) = Gaussian(σ, μ), ConservativeForm, explicit Euler" begin
  dd = DomainDecomp(gm=1)
  dss = DomainSubSet(gm=true)
  σ, μ = .1, 0.5
  δz = 0.2
  Triangle(z) = μ-δz < z < μ+δz ? ( μ > z ? (z-(μ-δz))/δz : ((μ+δz)-z)/δz) : 0.0
  Square(z) = μ-δz < z < μ+δz ? 1 : 0.0
  Gaussian(z) = exp(-1/2*((z-μ)/σ)^2)
  for n_elems_real in (64, 128)
    grid = Grid(0.0, 1.0, n_elems_real)
    domain_range = over_elems_real(grid)
    x = [grid.zc[k] for k in domain_range]
    unknowns = ( (:ϕ, dss), )
    vars = ( (:ϕ_initial, dss), (:ϕ_error, dss), (:ϕ_analytic, dss), )
    q = StateVec(unknowns, grid, dd)
    tmp = StateVec(vars, grid, dd)
    rhs = deepcopy(q)
    CFL = 0.1
    Δt = CFL*grid.Δz
    T = 0.25
    maxiter = Int(T/Δt)
    for scheme in (
                   UpwindAdvective(),
                   UpwindCollocated(),
                   CenteredUnstable(),
                   )
      for distribution in (Triangle,
                           Square,
                           Gaussian)
        for wave_speed in (-1, 1)

          scheme_name = replace(string(scheme), "()"=>"")
          distribution_name = joinpath("LinearAdvection",string(distribution))
          directory = joinpath(adv_eq_dir, distribution_name, scheme_name)
          export_plots && mkpath(directory)
          print_norms && print("\n", directory, ", ", n_elems_real, ", ", wave_speed, ", ")

          for k in over_elems_real(grid)
            tmp[:ϕ_initial, k] = distribution(grid.zc[k])
            q[:ϕ, k] = tmp[:ϕ_initial, k]
          end
          amax_w = max([max(tmp[:ϕ_initial, k]) for k in over_elems_real(grid)]...)
          for i in 1:maxiter
            for k in over_elems_real(grid)
              ϕ = q[:ϕ, Cut(k)]
              ϕ_dual = q[:ϕ, Dual(k)]
              w = [wave_speed, wave_speed, wave_speed]
              w_dual = [wave_speed, wave_speed]
              rhs[:ϕ, k] = - advect_old(ϕ, ϕ_dual, w, w_dual, grid, scheme, Δt)
            end
            for k in over_elems(grid)
              q[:ϕ, k] += Δt*rhs[:ϕ, k]
            end
            apply_Dirichlet!(q, :ϕ, grid, 0.0, Zmax())
            apply_Dirichlet!(q, :ϕ, grid, 0.0, Zmin())
          end
          sol_analtyic = [distribution(grid.zc[k] - wave_speed*T) for k in over_elems(grid)]
          sol_error = [sol_analtyic[k] - q[:ϕ, k] for k in over_elems(grid)]
          L2_norm = sum(sol_error.^2)/length(sol_error)

          if !(scheme==CenteredUnstable())
            @test all(abs.(sol_error) .< 100*grid.Δz)
            @test all(L2_norm < 100*grid.Δz)
          end
          for k in over_elems(grid)
            tmp[:ϕ_error, k] = sol_error[k]
            tmp[:ϕ_analytic, k] = sol_analtyic[k]
          end

          print_norms && print("L2_norm(err) = ", L2_norm)

          name = string(n_elems_real)
          markershape = wave_speed==-1 ? :dtriangle : :utriangle
          if export_plots
            if wave_speed==-1
              y = [tmp[:ϕ_initial , k] for k in domain_range]; plot(y, x, label="initial",
                markercolor="blue", linecolor="blue", markershapes=markershape, markersize=2, legend=:topleft)
              y = [  q[:ϕ         , k] for k in domain_range]; plot!(y, x, label="numerical",
                markercolor="black", linecolor="black", markershapes=markershape, markersize=2, legend=:topleft)
              y = [tmp[:ϕ_error   , k] for k in domain_range]; plot!(y, x, label="error",
                markercolor="red", linecolor="red", markershapes=markershape, markersize=2, legend=:topleft)
              y = [tmp[:ϕ_analytic, k] for k in domain_range]; plot!(y, x, label="analytic",
                markercolor="green", linecolor="green", markershapes=markershape, markersize=2, legend=:topleft)
            else
              y = [  q[:ϕ         , k] for k in domain_range]; plot!(y, x, label="",
                markercolor="black", linecolor="black", markershapes=markershape, markersize=2, legend=:topleft)
              y = [tmp[:ϕ_error   , k] for k in domain_range]; plot!(y, x, label="",
                markercolor="red", linecolor="red", markershapes=markershape, markersize=2, legend=:topleft)
              y = [tmp[:ϕ_analytic, k] for k in domain_range]; plot!(y, x, label="",
                markercolor="green", linecolor="green", markershapes=markershape, markersize=2, legend=:topleft)
            end
            wave_speed==1 && savefig(joinpath(directory, name))
          end
        end
      end
    end
  end
end

@testset "Non-linear Bergers, ∂_t w + ∇•(ww) = 0 ∈ ∂Ω, u(t=0) = Gaussian(σ, μ), ConservativeForm, explicit Euler" begin
  dd = DomainDecomp(gm=1)
  dss = DomainSubSet(gm=true)
  σ, μ = .1, 0.5
  δz = 0.2
  Triangle(z, velocity_sign)   = μ-δz < z < μ+δz ? ( μ > z ? velocity_sign*(z-(μ-δz))/δz : velocity_sign*((μ+δz)-z)/δz  ) : 0.0
  Square(z, velocity_sign) = μ-δz < z < μ+δz ? velocity_sign : 0.0
  Gaussian(z, velocity_sign) = velocity_sign*exp(-1/2*((z-μ)/σ)^2)
  for n_elems_real in (64, 128)
    grid = Grid(0.0, 1.0, n_elems_real)
    domain_range = over_elems_real(grid)
    x = [grid.zc[k] for k in domain_range]
    unknowns = ( (:w, dss), )
    vars = ( (:w_initial, dss), (:w_error, dss), (:w_analytic, dss), )
    q = StateVec(unknowns, grid, dd)
    tmp = StateVec(vars, grid, dd)
    rhs = deepcopy(q)
    CFL = 0.1
    Δt = CFL*grid.Δz
    T = 0.25
    maxiter = Int(T/Δt)
    for scheme in (
                   UpwindAdvective(),
                   UpwindCollocated(),
                   CenteredUnstable(),
                   )
      for distribution in (Triangle,
                           Square,
                           Gaussian)
        for velocity_sign in (-1, 1)
          for k in over_elems_real(grid)
            tmp[:w_initial, k] = distribution(grid.zc[k], velocity_sign)
            q[:w, k] = tmp[:w_initial, k]
          end
          amax_w = max([max(tmp[:w_initial, k]) for k in over_elems_real(grid)]...)
          for i in 1:maxiter
            for k in over_elems_real(grid)
              w = q[:w, Cut(k)]
              w_dual = (w[1:2]+w[2:3])/2
              rhs[:w, k] = - advect_old(w, w_dual, w, w_dual, grid, scheme, Δt)
            end
            for k in over_elems(grid)
              q[:w, k] += Δt*rhs[:w, k]
            end
            apply_Dirichlet!(q, :w, grid, 0.0, Zmax())
            apply_Dirichlet!(q, :w, grid, 0.0, Zmin())
          end

          scheme_name = replace(string(scheme), "()"=>"")
          distribution_name = joinpath("BurgersEquation",string(distribution))
          directory = joinpath(adv_eq_dir, distribution_name, scheme_name)
          export_plots && mkpath(directory)
          name = string(n_elems_real)
          markershape = velocity_sign==-1 ? :dtriangle : :utriangle
          if export_plots
            if velocity_sign==-1
              y = [tmp[:w_initial , k] for k in domain_range]; plot(y, x, label="initial",
                markercolor="blue", linecolor="blue", markershapes=markershape, markersize=2, legend=:topleft)
              y = [  q[:w         , k] for k in domain_range]; plot!(y, x, label="numerical",
                markercolor="black", linecolor="black", markershapes=markershape, markersize=2, legend=:topleft)
            else
              y = [tmp[:w_initial , k] for k in domain_range]; plot!(y, x, label="",
                markercolor="blue", linecolor="blue", markershapes=markershape, markersize=2, legend=:topleft)
              y = [  q[:w         , k] for k in domain_range]; plot!(y, x, label="",
                markercolor="black", linecolor="black", markershapes=markershape, markersize=2, legend=:topleft)
            end
            velocity_sign==1 && savefig(joinpath(directory, name))
          end
        end
      end
    end
  end

  end

ext: jl

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
using DataFrames, DataFramesMeta, Statistics
using PyPlot
using Colors

# Functions/Helpers to process dataframe files (outputs of simulations)

# Retain results for only one combination of
# stochasticy and changeprobability values
function df_trim(df, svalue, changeprobabilityvalue)
    dftrimmed = []
    if in(:sigma, names(df))
        dftrimmed = @where(df, :changeprobability .== changeprobabilityvalue,
                            :sigma .== svalue)
    elseif in(:stochasticity, names(df))
        dftrimmed = @where(df, :changeprobability .== changeprobabilityvalue,
                            :stochasticity .== svalue)
    end
    dftrimmed
end
# Merge seeds
function df_mergeseeds_errors(df)
    if in(:sigma, names(df))
        df = by(df, [:learner, :sigma, :changeprobability],
                allmerrors = :merrors => x -> [collect(skipmissing(reduce(vcat, x)))])
    elseif in(:stochasticity, names(df))
        df = by(df, [:learner, :stochasticity, :changeprobability],
                allmerrors = :terrors => x -> [collect(skipmissing(reduce(vcat, x)))])
    end
    df
end
# Get the mean error across all time steps and all seeds
function df_mergeseeds_meanerrors(df)
    if in(:sigma, names(df))
        df = by(df, [:learner, :sigma, :changeprobability],
                meanerrors = :merrors => x -> mean(collect(skipmissing(reduce(vcat, x)))))
    elseif in(:stochasticity, names(df))
        df = by(df, [:learner, :stochasticity, :changeprobability],
                meanerrors = :terrors => x -> mean(collect(skipmissing(reduce(vcat, x)))))
    end
    df
end
# Get the mean error across time steps for a single seed
function df_meanerrors_perseed(df)
    if in(:sigma, names(df))
        df = by(df, [:learner, :sigma, :changeprobability, :seed],
                meanerrors_seed = :merrors => x -> mean(collect(skipmissing(reduce(vcat, x)))))
    elseif in(:stochasticity, names(df))
        df = by(df, [:learner, :stochasticity, :changeprobability, :seed],
                meanerrors_seed = :terrors => x -> mean(collect(skipmissing(reduce(vcat, x)))))
    end
    df
end
# Get the mean error across time steps for each seed
function df_meanerrors_acrosseeds(df)
    df = df_meanerrors_perseed(df)
    if in(:sigma, names(df))
        df = by(df, [:learner, :sigma, :changeprobability],
                allmeanerrors_seed = :meanerrors_seed => x -> [collect(skipmissing(reduce(vcat, x)))])
    elseif in(:stochasticity, names(df))
        df = by(df, [:learner, :stochasticity, :changeprobability],
                allmeanerrors_seed = :meanerrors_seed => x -> [collect(skipmissing(reduce(vcat, x)))])
    end
    df
end
# Get the mean, std, se across seeds (Call df_meanerrors_perseed(df) first)
function df_stats_acrossseeds(df_meanerrors_seed) #
    if in(:sigma, names(df_meanerrors_seed))
        df_meanerrors_seed = by(df_meanerrors_seed, [:learner, :sigma, :changeprobability],
                meanerrors_seed = :meanerrors_seed => x -> mean(collect(skipmissing(reduce(vcat, x)))),
                stderrors_seed = :meanerrors_seed => x -> std(collect(skipmissing(reduce(vcat, x)))),
                nseed = :meanerrors_seed => x -> length(collect(skipmissing(reduce(vcat, x)))),
                seerrors_seed = :meanerrors_seed => x -> std(collect(skipmissing(reduce(vcat, x)))) / sqrt(length(collect(skipmissing(reduce(vcat, x)))))
                )
    elseif in(:stochasticity, names(df_meanerrors_seed))
        df_meanerrors_seed = by(df_meanerrors_seed, [:learner, :stochasticity, :changeprobability],
                meanerrors_seed = :meanerrors_seed => x -> mean(collect(skipmissing(reduce(vcat, x)))),
                stderrors_seed = :meanerrors_seed => x -> std(collect(skipmissing(reduce(vcat, x)))),
                nseed = :meanerrors_seed => x -> length(collect(skipmissing(reduce(vcat, x)))),
                seerrors_seed = :meanerrors_seed => x -> std(collect(skipmissing(reduce(vcat, x)))) / sqrt(length(collect(skipmissing(reduce(vcat, x)))))
                )
    end
    df_meanerrors_seed
end
function formatparamlist(paramlist)
    tuple([isa(x, AbstractFloat) ? round(x, digits=3) : x for x in paramlist]...)
end
# Keep K number of steps after a change point
function keep_Ksteps(dataafterswitches, Ksteps)
    Nsteps = [size(x, 1) for x in dataafterswitches]
    toosmall = findall(Nsteps .<= Ksteps)
    deleteat!(dataafterswitches, toosmall)
    dataafterswitches = [dataafterswitches[i][1:Ksteps]
                    for i in 1:length(dataafterswitches)]
end
function df_getlearnerinfo(df)
    unique(df[:learner])[1], unique(df[:param])[1]
end
function calc_stats(dataarray) # dataarray : timepoints x instances
    m = vec(mean(dataarray, dims=2))
    s = vec(std(dataarray, dims=2))
    serror = s ./ sqrt(size(dataarray, 2))
    m, s, serror
end
nanmean(x) = mean(filter(!isnan,x))
nanmean(x,dims) = mapslices(nanmean, x, dims=dims)
nanstd(x) = std(filter(!isnan,x))
nanstd(x,dims) = mapslices(nanstd, x, dims=dims)
nanlength(x) = length(filter(!isnan,x))
nanlength(x,dims) = mapslices(nanlength, x, dims=dims)
function nanserror(x, dims)
    std = nanstd(x, dims)
    n = nanlength(x, dims)
    std ./ sqrt.(n)
end
function calc_stats_nan(dataarray) # dataarray : timepoints x instances
    m = vec(nanmean(dataarray, 2))
    s = vec(nanstd(dataarray, 2))
    serror = vec(nanserror(dataarray, 2))
    m, s, serror
end
# Get the mean error after change points (in time)
# For Figure 2 and Figure 5
function df_meanerrorsafterswitches(df)
    df2 = df_geterrorsafterswitches(df)
    erallswitches = vcat(df2[:errorafterswitches]...)
    erallswitches = hcat(erallswitches...)
    merafterswitches, stderafterswitches, stderrorerafterswitches = calc_stats(erallswitches)
end
# Get the error after change points (in time) for each row (enviroment instance)
# in the dataframe
function df_geterrorsafterswitches(df; Ksteps = 50)
    if in(:merrors, names(df))
    df2 = @byrow! df begin
        @newcol errorafterswitches::Array{Array{Array{Float64,1} ,1}, 1}
        :errorafterswitches = geterrorsafterswitches(:merrors, :switches,
                                Int(1. / :changeprobability))
        end
    elseif in(:terrors, names(df))
    df2 = @byrow! df begin
        @newcol errorafterswitches::Array{Array{Array{Float64,1} ,1}, 1}
        :errorafterswitches = geterrorsafterswitches(:terrors, :switches,
                            Int(1. / :changeprobability))
        end
    end
end
# Get the errors in time (whole simulation) and align them with respect to
# the change points
function geterrorsafterswitches(errors, switches, Ksteps)
    switchindices = find_switchpoints(errors, switches)
    #[a[b[i]:b[i+1]] for i in 1:length(b)-1]
    erafterswitches = [errors[switchindices[i]:switchindices[i+1]]
                        for i in 1:length(switchindices)-1]
    erafterswitches = keep_Ksteps(erafterswitches, Ksteps)
end
# Find the time steps of change points
function find_switchpoints(errors, switches)
    switchindices = findall(switches)
    switchindices = vcat(switchindices, length(errors))
end
# ------------------------------------------------------------------------------
# Function to plot (or create csv data in order to plot them)
# the lower left inset plots of Figure 5 (task example)
function plottask_M(dflearner, stochasticity, changeprobability, seed;
                    savepath = homedir())
    states = @where(dflearner,
        :stochasticity .== stochasticity,
        :changeprobability .== changeprobability,
        :seed .== seed
        )[:states]

    switches = @where(dflearner,
            :stochasticity .== stochasticity,
            :changeprobability .== changeprobability,
            :seed .== seed
            )[:switches]

    fig = figure(); ax = gca()
    data = states[1][1:100]
    switchdata = switches[1][1:100]
    iswitch = findall(switchdata)
    ax.plot(1:length(data), data)

    if !isempty(switchdata)
        for i in 1:length(iswitch)
            ax.axvline(iswitch[i], color="red", linestyle="--", linewidth=0.5)
        end
    end
    ax.grid(true)
    xlabel("Time steps")
    ylabel("Observed state")
    PyPlot.savefig(joinpath(savepath, "states.png"))

    dfresults = DataFrame()
    dfresults[:x] = collect(1:length(states[1]))
    dfresults[:states] = states[1]
    dfresults[:switches] = Int.(switches[1])

    CSV.write(joinpath(savepath, "task_multi_s"* string(stochasticity) * "_Pc" *string(changeprobability)*".csv"), dfresults,
        delim = " ");
end
# Function to plot (or create csv data in order to plot them)
# the lower left inset plot of Figure 2 (task example)
function plottask_G(dflearner, sigma, changeprobability, seed;
                        savepath = homedir())

    learnertype = GLeakyIntegrator
    parameterset = (0.9,)

    callback = runG!(learnertype, parameterset,
            sigma = sigma,
            changeprobability = changeprobability,
            nsteps = 10^3,
            iseedenvpolicy = seed,
            iseedlearner = seed)

    true_mean = callback.mu_n_history[2:end]
    # true_mean = callback.mu_n_history[2:300]
    samples = callback.state_history[2:end]
    #samples = callback.state_history[2:300]
    switches = Int.(callback.switchflag)
    #switches = Int.(callback.switchflag[1:299])
    fig = figure(); ax = gca()
    ax.plot(true_mean)
    ax.scatter(1:length(samples), samples, color = "r")
    #iswitch = findall(callback.switchflag[1:299])
    iswitch = findall(callback.switchflag)
    if !isempty(switches)
        for i in 1:length(iswitch)
            ax.axvline(iswitch[i], color="red", linestyle="--", linewidth=0.5)
        end
    end
    PyPlot.savefig(joinpath(savepath, "gaussian_states.png"))

    dfresults = DataFrame()
    dfresults[:x] = collect(1:length(samples))
    dfresults[:truemean] = true_mean
    dfresults[:samples] = samples
    dfresults[:switches] = switches

    CSV.write(joinpath(savepath, "task_gaussian_s"* string(sigma) * "_Pc" *string(changeprobability)*".csv"), dfresults,
        delim = " ");
end
# ------------------------------------------------------------------------------
# Helper functions to explore learners
function explorelearner(df, learnername, stochasticity, stay, seed)
    dflearner = @where(df,
        :learner .== learnername,
        :stochasticity .== stochasticity,
        :changeprobability .== changeprobability,
        :seed .== seed
        )
end
function explorelearner_terrors(df, learnername, stochasticity, changeprobability)
    terrors = @where(df,
        :learner .== learnername,
        :stochasticity .== stochasticity,
        :changeprobability .== changeprobability,
        )[:terrors]
    terrors = reduce(vcat, terrors);
end
function explorelearner_merrors(df, learnername, sigma, changeprobability)
    merrors = @where(df,
        :learner .== learnername,
        :sigma .== sigma,
        :changeprobability .== changeprobability,
        )[:merrors]
    merrors = reduce(vcat, merrors);
end
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# --------- Plot within julia: not used eventually for the paper ---------------
function plot_shadedline(m, s, labelstring, color; titlestring = "")
    fig = figure(); ax = gca()
    plot_shadedline_core(m, s, labelstring, color, ax)
    title(titlestring)
end
function plot_shadedlines(means, stderrors, labels, ax;
                            xlabelstring = "Time steps after switch",
                            ylabelstring = "",
                            titlestring = "",
                            xdata = [],
                            ylims = [0., 1.])
    clist = [cmap(i) for i in range(0, 1, length = length(labels)+1)]
    for i in 1:length(labels)
        plot_shadedline_core(means[i], stderrors[i], labels[i], clist[i], ax,
                            xdata = xdata)
    end
    xlabel(xlabelstring)
    ylabel(ylabelstring)
    ax.set_ylim(ylims[1], ylims[2])
    title(titlestring)
    ax.legend()

end
function plot_shadedline_core(m, s, labelstring, color, ax; xdata = [])
    if isempty(xdata)
        ax.plot(1:length(m), m, label = labelstring, color = color)
        ax.fill_between(1:length(m), m + s, m - s, facecolor = color, alpha=0.5)
    else
        ax.plot(xdata, m, label = labelstring, color = color)
        ax.fill_between(xdata, m + s, m - s, facecolor = color, alpha=0.5)
    end
    ax.grid(true)
end
function plot_subplots_shadedlines(means, stderrors, labels, i; titlestring = "")
    subplot(1,5,i)
    ax = gca()
    plot_shadedlines(means, stderrors, labels, ax, titlestring = titlestring)
end


text:
# Based on https://github.com/cloud-oak/Tqdm.jl by @cloud-oak under Mozilla Public License 2.0
# Modified for Knet by Deniz Yuret

using Printf
import Base: length, size, iterate, eltype, IteratorSize, IteratorEltype, haslength, @propagate_inbounds

"""
    progress(f, itr; steps, seconds)
    progress(itr; o...) do x; f(x); end
    progress(itr; o...)
    progress!(...)

Return an iterator which acts exactly like `itr`, but prints a progressbar:

    ┣█████████████████▎  ┫ [86.83%, 903/1040, 01:36/01:50, 9.42i/s] 3.87835

Here `86.83%` is the percentage completed, `903` is the number of iterations completed,
`1040` is the total number of iterations. `01:36` is elapsed time, `01:50` is the estimated
total time, `9.42i/s` is the average number of iterations completed per second. If the speed
is less than 1, the average number of seconds per iteration (s/i) is reported instead.  The
bar, percent, total iterations, and estimated total time are omitted for iterators whose
size is unknown. `3.87835` is the output of `f` applied to the last value generated by itr.

The progress bar is updated and `f` is called with the most recent value of itr every
`steps` iterations or every `seconds` seconds in addition to the first and the last
iteration. If neither `steps` nor `seconds` is specified the default is to update every
second.

`f` can be specified by the first two forms above, if not specified (the third form) nothing
gets printed at the end of the line.  The last form, `progress!(...)`, is equivalent to
`(for x in progress(...) end)`, i.e.  iterates over the object created by `progress(...)`
and returns `nothing`.

"""
progress, progress!

mutable struct Progress{I}
    func
    iter::I
    steps
    seconds
    starttime::UInt
    lasttime::UInt
    lastiter::UInt
    curriter::UInt
    currval
end

progress(func::Base.Callable, iter::I; steps=0, seconds=0) where {I} =
    Progress{I}(func,iter,steps,seconds,0,0,0,0,nothing)

progress(iter; o...)=progress((x)->"",iter; o...)
progress!(x...; o...)=(for _ in progress(x...; o...) end)

IteratorSize(::Type{Progress{I}}) where {I} = IteratorSize(I)
IteratorEltype(::Type{Progress{I}}) where {I} = Base.EltypeUnknown()
length(p::Progress) = length(p.iter)

@propagate_inbounds function iterate(p::Progress, s...)
    if p.starttime == 0
        p.starttime = p.lasttime = time_ns()
        p.steps == p.seconds == 0 && (p.seconds = 1)
    end
    next = iterate(p.iter, s...)
    if next !== nothing
        p.curriter += 1
        (p.currval, s) = next
    end
    if ((next === nothing && p.curriter > 0) ||
        p.curriter == 1 || p.curriter == p.steps ||
        (p.steps != 0 && p.curriter == p.lastiter + p.steps) ||
        (p.seconds != 0 && time_ns() > p.lasttime + p.seconds*1e9))
        progressbar(p, next)
    end
    return next
end

function progressbar(p::Progress, next)
    fval_string = string(p.func(p.currval))
    currtime = time_ns()
    seconds = (currtime - p.starttime) * 1e-9
    speed = (next === nothing ? p.curriter / seconds : (p.curriter - p.lastiter) / ((currtime - p.lasttime) * 1e-9))
    p.lastiter, p.lasttime = p.curriter, currtime

    if haslength(p)
        ETA = (length(p) - p.curriter) / (p.curriter / seconds)
        percentage_string = string(@sprintf("%.2f%%",p.curriter/length(p)*100))
        status_string = string("[", percentage_string,
                               ", ", p.curriter, "/", length(p),
                               ", ", format_time(seconds), "/", format_time(seconds+ETA),
                               ", ", format_speed(speed),
                               "] ")
    else
        status_string = string("[", p.curriter,
                               ", ", format_time(seconds),
                               ", ", format_speed(speed),
                               "] ")
    end

    print("\r")

    if (haslength(p))
        width = 20
        print("┣")
        cellvalue = length(p) / width
        full_cells, remain = divrem(p.curriter, cellvalue)
        full_cells = round(Int, full_cells)
        print(repeat("█", full_cells))
        if (full_cells < width)
	    part = floor(Int, 8 * remain / cellvalue)
	    print(EIGHTS[part])
            print(repeat(" ", width - full_cells - 1))
        end
        print("┫ ")
    end

    print(status_string)
    print(fval_string)
    next === nothing && println()
    flush(stdout)
    return next
end

function format_time(seconds)
    if seconds != Inf
        mins,s  = divrem(round(Int,seconds), 60)
        h, m    = divrem(mins, 60)
    else
        h=0;m=Inf;s=Inf
    end
    if h!=0
         return @sprintf("%02d:%02d:%02d",h,m,s)
    else
         return @sprintf("%02d:%02d",m,s)
    end
end

format_speed(s)=(s >= 1 ? @sprintf("%.2fi/s",s) : @sprintf("%.2fs/i",1/s))

EIGHTS = Dict(0 => ' ',
	      1 => '▏',
	      2 => '▎',
	      3 => '▍',
	      4 => '▌',
	      5 => '▋',
	      6 => '▊',
	      7 => '▉',
	      8 => '█')


text:
@testset "ess.jl" begin
    @testset "copy and split" begin
        # check a matrix with even number of rows
        x = rand(50, 20)

        # check incompatible sizes
        @test_throws DimensionMismatch MCMCDiagnosticTools.copyto_split!(
            similar(x, 25, 20), x
        )
        @test_throws DimensionMismatch MCMCDiagnosticTools.copyto_split!(
            similar(x, 50, 40), x
        )

        y = similar(x, 25, 40)
        MCMCDiagnosticTools.copyto_split!(y, x)
        @test reshape(y, size(x)) == x

        # check a matrix with odd number of rows
        x = rand(51, 20)

        # check incompatible sizes
        @test_throws DimensionMismatch MCMCDiagnosticTools.copyto_split!(
            similar(x, 25, 20), x
        )
        @test_throws DimensionMismatch MCMCDiagnosticTools.copyto_split!(
            similar(x, 51, 40), x
        )

        MCMCDiagnosticTools.copyto_split!(y, x)
        @test reshape(y, 50, 20) == x[vcat(1:25, 27:51), :]
    end

    @testset "ESS and R̂ (IID samples)" begin
        rawx = randn(10_000, 40, 10)

        # Repeat tests with different scales
        for scale in (1, 50, 100)
            x = scale * rawx

            ess_standard, rhat_standard = ess_rhat(x)
            ess_standard2, rhat_standard2 = ess_rhat(x; method=ESSMethod())
            ess_fft, rhat_fft = ess_rhat(x; method=FFTESSMethod())
            ess_bda, rhat_bda = ess_rhat(x; method=BDAESSMethod())

            # check that we get (roughly) the same results
            @test ess_standard == ess_standard2
            @test ess_standard ≈ ess_fft
            @test rhat_standard == rhat_standard2 == rhat_fft == rhat_bda

            # check that the estimates are reasonable
            @test all(x -> isapprox(x, 100_000; rtol=0.1), ess_standard)
            @test all(x -> isapprox(x, 100_000; rtol=0.1), ess_bda)
            @test all(x -> isapprox(x, 1; rtol=0.1), rhat_standard)

            # BDA method fluctuates more
            @test var(ess_standard) < var(ess_bda)
        end
    end

    @testset "ESS and R̂ (identical samples)" begin
        x = ones(10_000, 40, 10)

        ess_standard, rhat_standard = ess_rhat(x)
        ess_standard2, rhat_standard2 = ess_rhat(x; method=ESSMethod())
        ess_fft, rhat_fft = ess_rhat(x; method=FFTESSMethod())
        ess_bda, rhat_bda = ess_rhat(x; method=BDAESSMethod())

        # check that the estimates are all NaN
        for ess in (ess_standard, ess_standard2, ess_fft, ess_bda)
            @test all(isnan, ess)
        end
        for rhat in (rhat_standard, rhat_standard2, rhat_fft, rhat_bda)
            @test all(isnan, rhat)
        end
    end

    @testset "ESS and R̂ (single sample)" begin # check that issue #137 is fixed
        x = rand(1, 5, 3)

        for method in (ESSMethod(), FFTESSMethod(), BDAESSMethod())
            # analyze array
            ess_array, rhat_array = ess_rhat(x; method=method)

            @test length(ess_array) == size(x, 2)
            @test all(ismissing, ess_array) # since min(maxlag, niter - 1) = 0
            @test length(rhat_array) == size(x, 2)
            @test all(ismissing, rhat_array)
        end
    end
end


text:
#=##############################################################################
# DESCRIPTION
    Tools for post-processing, plotting, and visualizing acoustic outputs.

# AUTHORSHIP
  * Author    : Eduardo J. Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Jan 2021
  * License   : MIT
=###############################################################################

"""
    `read_data(dataset_infos; datasets_psw=def_datasets_psw,
                datasets_bpm=def_datasets_bpm,
                psw_fieldnames=["pressure", "spl_spectrum", "OASPLdB", "OASPLdBA"],
                bpm_fieldnames=["spl_spectrum", "splA_spectrum", "OASPLdB", "OASPLdBA"],
                optargs...)`

Read a series of datasets comprised of PSU-WOPWOP and BPM outputs and stores
them under the dictionaries `datasets_psw` and `datasets_bpm`.

# ARGUMENTS
* `dataset_infos`    : Array of tuples (see example).

# OPTIONAL ARGUMENTS
* `psw_fieldnames`   : Fields to read from PSU-WOPWOP.
* `bpm_fieldnames`   : Fields to read from BPM.
* `optargs...`       : Optional arguments passed to `fetch_pswdataset(...)` and
                        `fetch_bpmdataset(...)`.

# EXAMPLE

```julia
# Dataset to read and associated information
dataset_infos = [   # (label, PSW solution, BPM solution, ...)

                    ("FLOWUnsteady w/BEMT",
                        "temps/example_dji9443_00_psw00/runcase/",
                        "temps/example_dji9443_00_bpm00",
                        "-", "r"),

                    ("FLOWUnsteady w/VPM",
                        "temps/example_dji9443_01_psw00/runcase/",
                        "temps/example_dji9443_01_bpm00",
                        "-", "b")
                ]

datasets_psw = Dict()     # Stores PSW data in this dictionary
datasets_bpm = Dict()     # Stores BPM data in this dictionary

# Read datasets and stores them in dictionaries
read_data(dataset_infos; datasets_psw=datasets_psw, datasets_bpm=datasets_bpm)
```
"""
function read_data(dataset_infos;
                    datasets_psw=def_datasets_psw, datasets_bpm=def_datasets_bpm,
                    psw_fieldnames=["pressure", "spl_spectrum", "spl_octFilt_spectrum", "OASPLdB", "OASPLdBA"],
                    bpm_fieldnames=["spl_spectrum", "splA_spectrum", "OASPLdB", "OASPLdBA", "frequencies"],
                    optargs...)

    # Read each dataset
    for (lbl, psw_read_path, bpm_read_path) in dataset_infos

        fetch_pswdataset(psw_read_path; datasets=datasets_psw, fieldnames=psw_fieldnames, optargs...)
        fetch_bpmdataset(bpm_read_path; datasets=datasets_bpm, fieldnames=bpm_fieldnames, optargs...)

    end
end

"""
`plot_pressure(dataset_infos, microphones, RPM, sph_ntht, pangle;
fieldname="pressure", datasets_psw=def_datasets_psw)`

Plots the pressure waveform of the tonal component of noise calculated by
PSU-WOPWOP at the requested microphone positions of a 2D grid.

# ARGUMENTS
* `dataset_infos`     : Array of tuples containing datasets to plot (see
                            example).
* `microphones`       : Array of microphone positions to plot.
* `RPM`               : RPM of the simulation (used to scale the time).
* `sph_ntht`          : Number of azimuthal nodes in the grid of observers (used
                            to convert `microphones` into node indices).
* `pangle::Function`  : Function to convert from indices of grid nodes to
                            degrees (or any other metric represented by the
                            azimuthal dimension of the grid).

# OPTIONAL ARGUMENTS
* `datasets_psw`      : Datasets where to read data from.

# EXAMPLE

```julia
# Dataset to read and associated information
dataset_infos = [   # (label, PSW solution, BPM solution, line style, color)

                    ("FLOWUnsteady w/BEMT",
                        "temps/example_dji9443_00_psw00/runcase/",
                        "temps/example_dji9443_00_bpm00",
                        "-", "r"),

                    ("FLOWUnsteady w/VPM",
                        "temps/example_dji9443_01_psw00/runcase/",
                        "temps/example_dji9443_01_bpm00",
                        "-", "b")
                ]

datasets_psw = Dict()     # Stores PSW data in this dictionary
datasets_bpm = Dict()     # Stores BPM data in this dictionary

# Read datasets and stores them in dictionaries
read_data(dataset_infos; datasets_psw=datasets_psw, datasets_bpm=datasets_bpm)

# Plot pressure waveform
microphones  = [-45, -90]            # (deg) microphones to plot
RPM          = 5400                  # RPM of solution

plot_pressure(dataset_infos, microphones, RPM, sph_ntht, pangle; datasets_psw=datasets_psw)
```
"""
function plot_pressure(dataset_infos, microphones, RPM, sph_ntht, pangle::Function;
                                fieldname="pressure", datasets_psw=def_datasets_psw,
                                xlims=[0, 3])
    # microphones  = [-45, -90]            # (deg) microphones to plot
    # RPM          = 5400                  # RPM of solution

    # fieldname = "pressure"               # Field to plot
    hash = Dict((round(pangle(mici), digits=1), mici) for mici in 1:sph_ntht)
    # mics         = Int.((-microphones .+ 180) * sph_ntht/360 .+ 1)   # Index of every microphone
    mics         = [hash[deg] for deg in microphones]

    for mici in mics

        plt.figure("$(fieldname)-$(mici)", figsize=[7*1*1.4, 3*3]*0.7)
        plt.suptitle("$(Int(round(pangle(mici))))"*L"$^\circ$ microphone")

        for (ploti, (ylbl, dlbl)) in enumerate([
                                                ("Total pressure (Pa)", "TotalAcousticPressure"),
                                                ("Loading pressure (Pa)", "LoadingAcousticPressure"),
                                                ("Thickness pressure (Pa)", "ThicknessAcousticPressure"),
                                              ])
            plt.subplot(310+ploti)
            if ploti==3; plt.xlabel("Rotor revolution"); end;
            plt.xlim(xlims)
            plt.ylabel(ylbl)
            plt.grid(true, which="major", color="0.8", linestyle="--")
            plt.grid(true, which="minor", color="0.8", linestyle="--")

            for (di, (lbl, read_path, _, stl, clr)) in enumerate(dataset_infos)

                data = datasets_psw[read_path][fieldname]
                xi = data["hs"]["ObserverTimes"]
                yi = data["hs"][dlbl]

                plt.plot(data["field"][mici, 1, 2:end, xi]*(RPM/60),
                         data["field"][mici, 1, 2:end, yi],  stl, alpha=0.8, label=lbl, color=clr)

            end

            if ploti==1; plt.legend(loc="center left", bbox_to_anchor=(1, 0.5)); end;
            # if ploti==1; plt.legend(loc="best"); end;
        end

        plt.tight_layout(rect=[0, 0.03, 1, 0.95])

    end
end

"""
    `plot_spectrum_spl(dataset_infos, microphones, BPF, sph_ntht, pangle::Function;
                                datasets_psw=def_datasets_psw,
                                datasets_bpm=def_datasets_bpm,
                                onethirdoctave=false,
                                Aweighted=false,
                                plot_csv=nothing, csv_lbl="Experimental", csv_stl="k", csv_weight=true,
                                xBPF=true, BPF_lines=0,
                                xlims=[6*10.0^(-1.0), 1.75e2], ylims=[0, 60],
                                Pref=L" re \$20\\mu\$Pa"
                                )`

Plot the SPL spectrum in `dataset_infos` at the given microphones.

# ARGUMENTS
* `dataset_infos`     : Array of tuples containing datasets to plot (see
                            example).
* `microphones`       : Array of microphone positions to plot.
* `BPF`               : Blade passing frequency (used to scale the frequency).
* `sph_ntht`          : Number of azimuthal nodes in the grid of observers (used
                            to convert `microphones` into node indices).
* `pangle::Function`  : Function to convert from indices of grid nodes to
                            degrees (or any other metric represented by the
                            azimuthal dimension of the grid).

# OPTIONAL ARGUMENTS
* `datasets_psw`      : PSU-WOPWOP datasets where to read data from.
* `datasets_bpm`      : BPM datasets where to read data from.
* `onethirdoctave`    : Converts the plot to 1/3 octave bands, otherwise it'll
                            convert BPM's 1/3 bands to narrow bands of PSW.
* `Aweighted`         : Plot A-weighted SPL.
* `plot_csv`          : Give it the path to a CSV file to include in the plot.
* `xBPF`              : Scale frequencies by BPF if true.
* `BPF_lines`         : Number of lines of BPF multiples to draw in the plot.

# EXAMPLE
```julia
# Dataset to read and associated information
dataset_infos = [   # (label, PSW solution, BPM solution, BPM freq bins, line style, color)
                    ("FLOWUnsteady w/BEMT",
                        "temps/example_dji9443_00_psw00/runcase/",
                        "temps/example_dji9443_00_bpm00",
                        "-", "r")
                ]

datasets_psw = Dict()     # Stores PSW data in this dictionary
datasets_bpm = Dict()     # Stores BPM data in this dictionary

# Read datasets and stores them in dictionaries
FLOWNoise.read_data(dataset_infos; datasets_psw=datasets_psw, datasets_bpm=datasets_bpm)

microphones  = [-45]                 # (deg) microphone to plot
Aweighted    = false                 # Plot A-weighted SPL

# Experimental data from Zawodny et al., Fig. 9
filename = joinpath("data", "zawodnydigitization2", "zawodny_dji9443_spl_5400_01.csv")
plot_experimental = [(filename, "Experimental", "k", Aweighted, [])]

# Plot spectrum in Hz with blade-passing-frequency lines
FLOWNoise.plot_spectrum_spl(dataset_infos, microphones, BPF, sph_ntht, pangle;
                              datasets_psw=datasets_psw, datasets_bpm=datasets_bpm,
                              Aweighted=Aweighted,
                              plot_csv=plot_experimental,
                              xBPF=false, xlims=[100, 3e4], BPF_lines=21)
```

"""
function plot_spectrum_spl(dataset_infos, microphones, BPF, sph_ntht, pangle::Function;
                            datasets_psw=def_datasets_psw,
                            datasets_bpm=def_datasets_bpm,
                            onethirdoctave=false,
                            psworg_onethirdoctave=false,
                            smooth_narrow=true,
                            Aweighted=false,
                            components=false,
                            add_broadband=true,
                            plot_csv=[],
                            xBPF=true, BPF_lines=0,
                            figname="automatic",
                            xlims=[6*10.0^(-1.0), 1.75e2], ylims=[0, 60],
                            Pref=L" re $20\mu$Pa",
                            verbose=true
                            )

    # microphone   = -45                   # (deg) microphone to plot
    # RPM          = 5400                  # RPM of solution
    # BPF          = 2*RPM/60              # Blade passing frequency

    fieldname_psw = onethirdoctave && psworg_onethirdoctave ? "spl_octFilt_spectrum" : "spl_spectrum"    # Field to plot
    fieldname_bpm = "spl"*"A"^Aweighted*"_spectrum"        # Field to plot

    hash = Dict((round(pangle(mici), digits=1), mici) for mici in 1:sph_ntht)
    # micis         = Int.((-microphones .+ 180) * sph_ntht/360 .+ 1)    # Index of each microphone
    micis         = [hash[deg] for deg in microphones]

    xscaling = (1/BPF)^xBPF

    aux = "SPL"*(onethirdoctave ? L"_{1/3}" : "")
    ylbl = Aweighted ? "A-weighted "*aux*"\n(dBA"*Pref*")" : aux*" (dB"*Pref*")"

    for mici in micis

        _figname = figname=="automatic" ? "$(fieldname_bpm)-$(mici)" : figname
        plt.figure(_figname, figsize=[7, 3.5])

        for (csv_filename, lbl, stl, weight, optargs) in plot_csv
            data = CSV.read(csv_filename, DataFrames.DataFrame, datarow=1)
            spl_exp = weight ? aWeight.(data[!, 1], data[!, 2]) : data[!, 2]
            plt.plot(data[!, 1]*xscaling, spl_exp, stl; label=lbl, optargs...)
        end

        # Plot datasets
        for (lbl, read_path_psw, read_path_bpm, stl, clr) in dataset_infos

            # Fetch tonal noise
            data_psw = datasets_psw[read_path_psw][fieldname_psw]
            xi = data_psw["hs"]["Frequency"]
            yi = data_psw["hs"]["Total_dB"*"A"^Aweighted]
            freqs_psw = data_psw["field"][mici, 1, 2:end, xi]
            org_spl_psw = data_psw["field"][mici, 1, 2:end, yi]

            # Fetch broadband noise
            data_bpm = datasets_bpm[read_path_bpm][fieldname_bpm]
            freqs_bpm = datasets_bpm[read_path_bpm]["frequencies"]["field"][:, 1]
            org_spl_bpm = data_bpm["field"][:, mici]

            # Convert narrow -> 1/3 or 1/3 -> narrow band to make tonal and
            # broadband comparable
            if onethirdoctave
                if !psworg_onethirdoctave
                    freqs_psw, org_spl_psw = narrow2onethird_spl(freqs_psw, org_spl_psw)
                end
            else
                df = (freqs_psw[2]-freqs_psw[1])
                println("Converting BPM to narrow band of bin width $(round(df, digits=1)) Hz")
                freqs_bpm, org_spl_bpm = onethird2narrow_spl(freqs_bpm, df, org_spl_bpm; smooth=smooth_narrow)
            end

            # Define the range of frequency as the union of both components
            freqs = freqs_psw                       # Grab PSW frequency range
            if freqs_psw[end] < freqs_bpm[end]  # Add BPM range
                bpm_i = findfirst( f -> f > freqs_psw[end], freqs_bpm)
                freqs = vcat(freqs, freqs_bpm[bpm_i:end])
                spl_psw = vcat(org_spl_psw, [-30.0 for i in bpm_i:length(freqs_bpm)])
            else
                spl_psw = org_spl_psw
            end

            # Interpolate broadband data into the same frequencies than tonal data
            spl_bpm = math.akima(freqs_bpm, org_spl_bpm, freqs)

            # Add tonal and broadband SPL together
            spl = addSPL(spl_psw, spl_bpm)

            if components
                plt.plot(freqs*xscaling, spl, "-", alpha=0.75, label=lbl*" --- Total", color=clr)
                plt.plot(freqs_bpm*xscaling, org_spl_bpm, ":", alpha=0.75, label=lbl*" --- BPM", color=clr)
                plt.plot(freqs_psw*xscaling, org_spl_psw, "--", alpha=0.75, label=lbl*" --- FW-H", color=clr)
            elseif !add_broadband
                plt.plot(freqs_psw*xscaling, org_spl_psw, stl, alpha=0.75, label=lbl, color=clr)
            else
                plt.plot(freqs*xscaling, spl, ("o"^onethirdoctave)*stl, alpha=0.75, label=lbl, color=clr)
            end

        end

        if BPF_lines != 0
            for BPFi in 1:BPF_lines
                this_x = BPFi*BPF*xscaling
                plt.plot(this_x*ones(100), range(ylims[1], ylims[2], length=100),
                                                               ",k", alpha=0.25)
                if BPFi<=8
                    plt.annotate("$BPFi"*(BPFi==1 ? L"$^\mathrm{st}$" :
                                          BPFi==2 ? L"$^\mathrm{nd}$" :
                                          BPFi==3 ? L"$^\mathrm{rd}$" : L"$^\mathrm{th}$"),
                                  (this_x, ylims[2]), (this_x, (1-0.1-0.001*BPFi^2.35)*(ylims[2]-ylims[1])), alpha=0.25)
                end
            end
        end

        plt.title("$(Int(round(pangle(mici))))"*L"$^\circ$ microphone")
        plt.xlim(xlims)
        plt.xscale("log")
        plt.xlabel("Frequency"*(xBPF ? " / BPF" : " (Hz)"))
        plt.ylim(ylims)
        plt.ylabel(ylbl)
        if BPF_lines==0 && BPF_lines!=-1
            plt.grid(true, which="major", color="0.8", linestyle=":")
            plt.grid(true, which="minor", color="0.8", linestyle=":")
        end
        plt.legend(loc="best", frameon=false)

        plt.tight_layout(rect=[0, 0.03, 1, 0.95])
    end
end

"""
    `plot_directivity_splbpf(dataset_infos, BPFi, BPF, pangle::Function;
                                        datasets_psw=def_datasets_psw,
                                        datasets_bpm=def_datasets_bpm,
                                        plot_csv=[],
                                        thetalims=180*[1, -1],
                                        thetagrids=collect(-180:45/2:180)[[3, 5, 7, 8, 9, 10, 11, 13, 15, 17]],
                                        rlbl_pos=90,
                                        rticks=40:4:52, rlims=[40, 54], rorigin=36,
                                        Pref=L" re \$20\\mu\$Pa")`

Plot the directivity of the requested multiple `BPFi` of the blade-passing
frequency `BPF`.

# ARGUMENTS
* `dataset_infos`     : Array of tuples containing datasets to plot (see
                            example).
* `pangle::Function`  : Function to convert from indices of grid nodes to
                            degrees (or any other metric represented by the
                            azimuthal dimension of the grid) as `pangle(i)`.

# OPTIONAL ARGUMENTS
* `datasets_psw`      : PSU-WOPWOP datasets where to read data from.
* `datasets_bpm`      : BPM datasets where to read data from.
* `plot_csv`          : Give it the path to a CSV file to include in the plot.

# EXAMPLE
```julia

# Dataset to read and associated information
dataset_infos = [   # (label, PSW solution, BPM solution, BPM freq bins, line style, color)
                    ("FLOWUnsteady w/BEMT",
                        "temps/example_dji9443_00_psw00/runcase/",
                        "temps/example_dji9443_00_bpm00",
                        "-", "r")
                ]

datasets_psw = Dict()     # Stores PSW data in this dictionary
datasets_bpm = Dict()     # Stores BPM data in this dictionary

# Read datasets and stores them in dictionaries
FLOWNoise.read_data(dataset_infos; datasets_psw=datasets_psw, datasets_bpm=datasets_bpm)

BPFi = 1                 # Multiple of blade-passing frequency to plot

# Experimental and computational data reported by Zawodny et al., Fig. 14
filename1 = joinpath("data", "zawodnydigitization2", "zawodny_fig14_topright_exp00.csv")
filename2 = joinpath("data", "zawodnydigitization2", "zawodny_fig14_topright_of00.csv")
filename3 = joinpath("data", "zawodnydigitization2", "zawodny_fig14_topright_pas00.csv")

plot_experimental = [(filename1, "Experimental", "ok", Aweighted, []),
                     (filename2, "OVERFLOW2", "-c", Aweighted, [(:alpha, 0.8)]),
                     (filename3, "BEM+ANOPP", ":b", Aweighted, [(:alpha, 0.8)])]

# Plot SPL directivity of first blade-passing frequency
FLOWNoise.plot_directivity_splbpf(dataset_infos, BPFi, BPF, pangle;
                                    datasets_psw=datasets_psw,
                                    datasets_bpm=datasets_bpm,
                                    plot_csv=plot_experimental,
                                    rticks=40:4:52, rlims=[40, 54], rorigin=36)
```
"""
function plot_directivity_splbpf(dataset_infos, BPFi, BPF, pangle::Function;
                                    datasets_psw=def_datasets_psw,
                                    datasets_bpm=def_datasets_bpm,
                                    plot_csv=[],
                                    thetalims=180*[1, -1],
                                    thetagrids=collect(-180:45/2:180)[[3, 5, 7, 8, 9, 10, 11, 13, 15, 17]],
                                    rlbl_pos=90,
                                    rticks=40:4:52, rlims=[40, 54], rorigin=36,
                                    Pref=L" re $20\mu$Pa")

    # BPFi         = 1                     # BPF multiple to plot
    fieldname_psw = "spl_spectrum"         # Field to plot
    fieldname_bpm = "spl_spectrum"         # Field to plot

    plt.figure("$(fieldname_bpm)-$(BPFi)", figsize=[3, 3]*2)

    for (csv_filename, lbl, stl, weight, optargs) in plot_csv
        data = CSV.read(csv_filename, DataFrames.DataFrame, datarow=1)
        plt.polar(pi/180*data[!, 2], data[!, 1], stl; label=lbl, optargs...)
    end

    # Plot datasets
    for (lbl, read_path_psw, read_path_bpm, stl, clr) in dataset_infos

        # Fetch tonal noise
        data_psw = datasets_psw[read_path_psw][fieldname_psw]
        yi = data_psw["hs"]["Total_dB"]
        fi = data_psw["hs"]["Frequency"]

        df = data_psw["field"][1, 1, 2, fi] - data_psw["field"][1, 1, 1, fi] # Frequency step
        freqi = ceil(Int, BPFi*BPF/df + 1)                           # Frequency index
        freq = data_psw["field"][1, 1, freqi, fi]                    # Frequency
        # _lbl = " @ $(Int(round(freq, digits=0))) Hz"                 # Frequency string
        _lbl = ""

        spl_psw = data_psw["field"][:, 1, freqi, yi]

        # Fetch broadband noise
        data_bpm = datasets_bpm[read_path_bpm][fieldname_bpm]
        freqs_bpm = datasets_bpm[read_path_bpm]["frequencies"]["field"][:, 1]
        spl_bpm = [math.akima(freqs_bpm, data_bpm["field"][:, mici], freq) for mici in 1:length(spl_psw)]

        # Add tonal and broadband SPL together
        # spl = addSPL(spl_psw, spl_bpm)
        # NOTE: Here we ignore the broadband component since it is in a 1/3 octave band
        spl = spl_psw

        angles = pangle.(1:length(spl))
        plt.polar(pi/180*angles, spl, stl, label=lbl*_lbl, alpha=0.8, color=clr)
    end

    ax = plt.gca()
    ax.set_theta_zero_location("W")
    ax.set_theta_direction(-1)
    ax.set_thetalim(pi/180*thetalims)
    ax.set_rlabel_position(rlbl_pos)
    ax.set_rticks(rticks)
    ax.set_rlim(rlims)
    ax.set_rorigin(rorigin)
    ax.set_thetagrids(thetagrids)
    ax.grid(true)
    plt.xlabel("SPL (dB"*Pref*")\n @ $(BPFi)"*(BPFi==1 ? L"$^\mathrm{st}$" :
                           BPFi==2 ? L"$^\mathrm{nd}$" :
                           BPFi==3 ? L"$^\mathrm{rd}$" : L"$^\mathrm{th}$")*" BPF ($(ceil(Int, BPFi*BPF)) hz)")
    plt.legend(loc="center left", bbox_to_anchor=(1.1, 0.5))
    plt.tight_layout();
end

"""
    `plot_directivity_oaspl(dataset_infos, pangle::Function;
                                        datasets_psw=def_datasets_psw,
                                        datasets_bpm=def_datasets_bpm,
                                        Aweighted=false,
                                        plot_csv=[],
                                        thetalims=180*[1, -1],
                                        thetagrids=collect(-180:45/2:180)[[3, 5, 7, 8, 9, 10, 11, 13, 15, 17]],
                                        rlbl_pos=90,
                                        rticks=40:10:70, rlims=[40, 72], rorigin=30,
                                        Pref=L" re \$20\\mu\$Pa")`

Plot overall SPL directivity.

# ARGUMENTS
* `dataset_infos`     : Array of tuples containing datasets to plot (see
                            example).
* `pangle::Function`  : Function to convert from indices of grid nodes to
                            degrees (or any other metric represented by the
                            azimuthal dimension of the grid) as `pangle(i)`.

# OPTIONAL ARGUMENTS
* `datasets_psw`      : PSU-WOPWOP datasets where to read data from.
* `datasets_bpm`      : BPM datasets where to read data from.
* `Aweighted`         : Plot A-weighted SPL.
* `plot_csv`          : Give it the path to a CSV file to include in the plot.

# EXAMPLE
```julia

# Dataset to read and associated information
dataset_infos = [   # (label, PSW solution, BPM solution, BPM freq bins, line style, color)
                    ("FLOWUnsteady w/BEMT",
                        "temps/example_dji9443_00_psw00/runcase/",
                        "temps/example_dji9443_00_bpm00",
                        "-", "r")
                ]

datasets_psw = Dict()     # Stores PSW data in this dictionary
datasets_bpm = Dict()     # Stores BPM data in this dictionary

# Read datasets and stores them in dictionaries
FLOWNoise.read_data(dataset_infos; datasets_psw=datasets_psw, datasets_bpm=datasets_bpm)

Aweighted = false

# Experimental and computational data reported by Zawodny et al., Fig. 12
filename1 = joinpath("data", "zawodnydigitization2", "zawodny_fig12_left_5400_00.csv")

plot_experimental = [(filename1, "Experimental", "ok", Aweighted, [])]

# Plot OASPL directivity
FLOWNoise.plot_directivity_oaspl(dataset_infos, pangle;
                                    datasets_psw=datasets_psw,
                                    datasets_bpm=datasets_bpm,
                                    Aweighted=Aweighted,
                                    plot_csv=plot_experimental,
                                    rticks=40:10:70, rlims=[40, 72], rorigin=30)
```
"""
function plot_directivity_oaspl(dataset_infos, pangle::Function;
                                    datasets_psw=def_datasets_psw,
                                    datasets_bpm=def_datasets_bpm,
                                    Aweighted=false,
                                    add_broadband=true,
                                    plot_csv=[],
                                    thetalims=180*[1, -1],
                                    thetagrids=collect(-180:45/2:180)[[3, 5, 7, 8, 9, 10, 11, 13, 15, 17]],
                                    rlbl_pos=90,
                                    rticks=40:10:70, rlims=[40, 72], rorigin=30,
                                    Pref=L" re $20\mu$Pa")

    fieldname_psw = "OASPLdB"*"A"^Aweighted # Field to plot
    fieldname_bpm = fieldname_psw           # Field to plot

    plt.figure("$(fieldname_bpm)", figsize=[3, 3]*2)

    for (csv_filename, lbl, stl, weight, optargs) in plot_csv
        data = CSV.read(csv_filename, DataFrames.DataFrame, datarow=1)
        plt.polar(pi/180*data[!, 2], data[!, 1], stl; label=lbl, optargs...)
    end

    # Plot datasets
    for (lbl, read_path_psw, read_path_bpm, stl, clr) in dataset_infos

        # Fetch tonal noise
        data_psw = datasets_psw[read_path_psw][fieldname_psw]
        yi = data_psw["hs"]["TotalOASPLdB"*"A"^Aweighted]
        oaspl_psw = data_psw["field"][:, 1, 1, yi]

        # Fetch broadband noise
        data_bpm = datasets_bpm[read_path_bpm][fieldname_bpm]
        oaspl_bpm = data_bpm["field"][:]

        # Add tonal and broadband OASPL together
        oaspl = addSPL(oaspl_psw, oaspl_bpm)

        angles = pangle.(1:length(oaspl))
        if add_broadband
            plt.polar(pi/180*angles, oaspl, stl, label=lbl, alpha=0.8, color=clr)
        else
            plt.polar(pi/180*angles, oaspl_psw, stl, label=lbl, alpha=0.8, color=clr)
        end
    end

    ax = plt.gca()
    ax.set_theta_zero_location("W")
    ax.set_theta_direction(-1)
    ax.set_thetalim(pi/180*thetalims)
    ax.set_rlabel_position(rlbl_pos)
    ax.set_rticks(rticks)
    ax.set_rlim(rlims)
    ax.set_rorigin(rorigin)
    ax.set_thetagrids(thetagrids)
    ax.grid(true)
    plt.xlabel("A-weighted "^Aweighted*"OASPL\n"*"(dB"*"A"^Aweighted*"$(Pref))")
    plt.legend(loc="center left", bbox_to_anchor=(1.1, 0.5))
    plt.tight_layout();
end

max_stars_repo_name: byuflowlab/FLOWNoise

max_issues_repo_name: byuflowlab/FLOWNoise

max_forks_repo_name: byuflowlab/FLOWNoise


text:
# Valkó, P.P. and Abate, J.
# Comparison of Sequence Accelerators for the Gaver Method of Numerical Laplace Transform Inversion
# Computers and Mathematics with Application, Vol. 48 (Iss.3-40) 2004 pp. 629-636
# Gaver Wynn rho method

# type GWR <: AbstractILt
#     func::Function
#     Nterms::Int
# end
# GWR(func::Function) = GWR(func, gwr_default_num_terms)

"""
    gwr(func::Function, t::AbstractFloat, M::Integer=gwr_default_num_terms)

Evaluate the inverse Laplace transform of `func` at the point `t`. Use `M` terms in the algorithm.
For `typeof(t)` is `Float64`, the default for `M` is `16`. For `BigFloat` the default is `64`.

If `BigFloat` precision is larger than default, try increasing `M`.

# Example

```jldoctest
julia> InverseLaplace.gwr( s -> 1/s^3, 3.0)
4.499985907607361
```

!!! note

    This function uses the Gaver-Wynn rho method.
    It evaluates `func` only for real arguments.

"""
function gwr(func, t, M)
    Dt = typeof(t)
    bM = convert(Dt, M)
    tau = log(convert(Dt, 2)) / t
    broken = false
    Fi = Array{Dt}(undef, 2 * M)
@inbounds  for i in 1: 2 * M
        Fi[i] = func(i * tau)
    end
    M1 = M
    G0 = zeros(Dt, M1 + 1)
    for n in 1:M
        sm = zero(Dt)
        bn = big(n)
        for i in 0:n
            bi = convert(Dt, i)
            sm += binomial(big(n), big(i)) * (-1)^i * Fi[n + i]
        end
        G0[n] = tau * SpecialFunctions.factorial(2 * bn) /   # FIXME: Do this more efficiently
            (SpecialFunctions.factorial(bn) * SpecialFunctions.factorial(bn - 1)) * sm
    end
    Gm = zeros(Dt, M1 + 1)
    Gp = zeros(Dt, M1 + 1)
    best = G0[M1]
    for k in 0:M1-2
        for n in (M1 - 2 - k):-1:0
            expr = G0[n + 2] - G0[n + 1]
            if expr == 0
                broken = true
                break
            end
            expr = Gm[n + 2] + (k + 1) / expr
            Gp[n + 1] = expr
            if isodd(k) && n == M1 - 2 - k
                best = expr
            end
        end
        if broken break end
        for n in 0:(M1-k)
            Gm[n + 1] = G0[n + 1]
            G0[n + 1] = Gp[n + 1]
        end
    end
    best
end

gwr(func, t::Float64) = gwr(func, t, gwr_default_num_terms)
const gwr_BigFloat_default_num_terms = 64
gwr(func, t::BigFloat) = gwr(func, t, gwr_BigFloat_default_num_terms)
gwr(func, t::Integer, args...) = gwr(func, BigFloat(t), args...)
gwr(func, t::Rational, args...) = gwr(func, float(t), args...)


text:
function json2df(s::String)
    s_json = parse(s)
    ks = Symbol.(collect(keys(s_json[1])))
    df = DataFrame(map(k -> k => [], ks)...)
    foreach(x -> push!(df, x), s_json)
    return df
end


text:
"""
    set_scale_sketch(s::T, pixels::Int) where {T <: Length}
    set_scale_sketch(s::T, pixels::Int) where {T <: Force}
    set_scale_sketch(s::T, pixels::Int) where {V <: Velocity}

To reset default 20m per screen height HE:
    set_scale_sketch(m), or any other length unit
    set_scale_sketch(:Length)

To reset to default 70 m/s per screen height in pixels HE:
    set_scale_sketch(m/s), or any other velocity unit
    set_scale_sketch(:Velocity)

To reset to default 20kN per screen height in pixels HE:
    set_scale_sketch(N), or any other force unit
    set_scale_sketch(:Force)

To reset all scale defaults:
    set_scale_sketch()
"""
function set_scale_sketch(x::FreeUnits)
    if 1x isa Length
        set_scale_sketch(:Length)
    elseif 1x isa Force
        set_scale_sketch(:Force)
    elseif 1x isa Velocity
        set_scale_sketch(:Velocity)
    else
        throw(DomainError(" set_scale_sketch with one argument only responds to units of length, force or velocity"))
    end
end
function set_scale_sketch(x::Symbol)
    if x == :Length
        set_scale_sketch(20m, HE)
    elseif x == :Force
        set_scale_sketch(20kN, HE)
    elseif x == :Velocity
        set_scale_sketch(70m/s, HE)
    else
        throw(DomainError(" set_scale_sketch with one argument only responds to (units of) Length, Force or Velocity"))
    end
end
set_scale_sketch(x::T, pixels::Int) where {T <: Length} =  push!(SCALE, :DISTANCE => x / pixels)
set_scale_sketch(x::T, pixels::Int) where {T <: Force} = push!(SCALE, :FORCE => x / pixels)
set_scale_sketch(x::T, pixels::Int) where {T <: Velocity} = push!(SCALE, :VELOCITY => x / pixels)
function set_scale_sketch()
    # Default scale
    set_scale_sketch(:Length)
    set_scale_sketch(:Force)
    set_scale_sketch(:Velocity)
end

"""
    scale_to_pt(q::Quantity)::Float64

Return the number of pixels (or pts, this may be ambiguous) corresponding to a given quantity.

    scale_to_pt(u::FreeUnits)::Float64

Return the number of pixels (or pts, this may be ambiguous) corresponding to one unit of u.

    To reset default 20m per screen height HE:
    set_scale_sketch()
    To reset default 20 kN per screen height HE:
    set_scale_sketch()
    To reset default 70 m/s per screen height HE:
    set_scale_sketch()
    70 m/s
"""
scale_to_pt(q::Length)::Float64 = upreferred(q / scale_pt_to_unit(oneunit(q)))
scale_to_pt(q::Velocity)::Float64 = upreferred(q / scale_pt_to_unit(oneunit(q)))
scale_to_pt(q::Force)::Float64 = upreferred(q / scale_pt_to_unit(oneunit(q)))
scale_to_pt(x)::Float64 = x

"""
    scale_pt_to_unit(u::FreeUnits)
"""
scale_pt_to_unit(u::FreeUnits) = scale_pt_to_unit(1∙u)
function scale_pt_to_unit(q::Length)
    @assert oneunit(q) == q "Ambiguous input. The scaling concerns one unit of the argument. Multiply by the number of units elsewhere."
    (q / unit(q)) ∙ unit(q)(SCALE[:DISTANCE])
end
function scale_pt_to_unit(q::Velocity)
    @assert oneunit(q) == q "Ambiguous input. The scaling concerns one unit of the argument. Multiply by the number of units elsewhere."
    (q / unit(q)) ∙ unit(q)(SCALE[:VELOCITY])
end
function scale_pt_to_unit(q::Force)
    @assert oneunit(q) == q "Ambiguous input. The scaling concerns one unit of the argument. Multiply by the number of units elsewhere."
    (q / unit(q)) ∙ unit(q)(SCALE[:FORCE])
end

"""
    set_figure_height(h::Int)

Updates global HE. Take care if functions are defined with closures (captured globals)
"""
function set_figure_height(h::Int)
    global HE = h
end
"""
    set_figure_width(w::Int)

Updates global WI. Take care if functions are defined with closures (captured globals)
"""
function set_figure_width(w::Int)
    global WI = w
end

function rotate(a::Angle)
    rotate(-ustrip( a |> rad))
end
"""
    x_y_iterators_at_pixels( physwidth = 10.0m, physheight = 4.0m, centered = true)
    ->(xs, ys)

Create position iterators with units corresponding to screen pixels.
    if centered = false:  (xs, ys)  >= 0m
    default: centered == true,  centered on (0.0m, 0.0m)

Use this to evaluate functions for visualization.
When looping over combinations, use ys as the inner loop for faster matrix lookup.
For image matrices, x vary over rows, y vary over columns.

NOTE! Replace this functionality with 'coords_spatial(img)' (google it)
"""
function x_y_iterators_at_pixels(;physwidth = 10.0m, physheight = 4.0m, centered = true)
    # Resolution for interpolation nodes
    nx = round(Int64, scale_to_pt(physwidth))
    ny = round(Int64, scale_to_pt(physheight))
    # Bounding box
    xmin = -centered * physwidth / 2
    xmax = (1 - centered / 2) * physwidth
    ymin = -centered * physheight / 2
    ymax = (1 - centered / 2) * physheight
    # Position iterators
    xs = range(xmin, xmax, length = nx)
    ys = range(ymax, ymin, length = ny) # Reversed order for faster matrix indexing.
    xs, ys
end


text:
"""
polyinterp(t,y)

Return a callable polynomial interpolant through the points in
vectors `t`,`y`. Uses the barycentric interpolation formula.
"""
function polyinterp(t,y)
    n = length(t)-1
    C = (t[end]-t[1]) / 4       # scaling factor to ensure stability
    tc = t/C

    # Adding one node at a time, compute inverses of the weights.
    omega = ones(n+1)
    for m = 1:n
        d = tc[1:m] .- tc[m+1]      # vector of node differences
        omega[1:m] = omega[1:m].*d  # update previous
        omega[m+1] = prod( -d )     # compute the new one
    end
    w = 1 ./ omega                  # go from inverses to weights

    p = function (x)
        # Compute interpolant.
        terms = @. w / (x - t)
        if any(isinf.(terms))       # divided by zero here
            # Apply L'Hopital's Rule exactly.
            idx = findfirst(x.==t)
            f = y[idx]
        else
            f = sum(y.*terms) / sum(terms)
        end
    end
    return p
end

"""
triginterp(t,y)

Return trigonometric interpolant for points defined by vectors `t`
and `y`.
"""
function triginterp(t,y)
    N = length(t)

    function trigcardinal(x)
        if isodd(N)      # odd
            tau = sin(N*pi*x/2) / (N*sin(pi*x/2))
        else             # even
            tau = sin(N*pi*x/2) / (N*tan(pi*x/2))
        end
        if isnan(tau)
            tau = 1
        end
        return tau
    end

    p = function (x)
        sum( y[k]*trigcardinal(x-t[k]) for k=1:N )
    end
    return p
end

"""
ccint(f,n)

Perform Clenshaw-Curtis integration for the function `f` on `n`+1
nodes in [-1,1]. Return integral and a vector of the nodes used.
Note: `n` must be even.
"""
function ccint(f,n)
    # Find Chebyshev extreme nodes.
    theta = [ i*pi/n for i=0:n ]
    x = -cos.(theta)

    # Compute the C-C weights.
    c = zeros(n+1)
    c[[1,n+1]] .= 1/(n^2-1)
    v = 1 .- 2*sum( cos.(2*k*theta[2:n])/(4*k^2-1) for k=1:n/2-1 )
    v -= cos.(n*theta[2:n])/(n^2-1)
    c[2:n] = 2*v/n

    # Evaluate integrand and integral.
    I = dot(c,f.(x))   # use vector inner product
    return I,x
end

"""
glint(f,n)

Perform Gauss-Legendre integration for the function `f` on `n` nodes
in (-1,1). Return integral and a vector of the nodes used.
"""
function glint(f,n)
    # Nodes and weights are found via a tridiagonal eigenvalue problem.
    beta = @. 0.5/sqrt(1-(2*(1:n-1))^(-2))
    T = diagm(-1=>beta,1=>beta)
    lambda,V = eigen(T)
    p = sortperm(lambda)
    x = lambda[p]           # nodes
    c = @. 2*V[1,p]^2       # weights

    # Evaluate the integrand and compute the integral.
    I = dot(c,f.(x))      # vector inner product
    return I,x
end

"""
intde(f,h,M)

Perform doubly-exponential integration of function `f` over
(-Inf,Inf), using discretization size `h` and truncation point `M`.
Return integral and a vector of the nodes used.
"""
function intde(f,h,M)
    # Find where to truncate the trapezoid sum.
    K = ceil( log(4/pi*log(2*M))/h )

    # Integrate by trapezoids in a transformed variable t.
    t = h*(-K:K)
    x = @. sinh(pi/2*sinh(t))
    dxdt = @. pi/2*cosh(t)*cosh(pi/2*sinh(t))

    I = h*dot(f.(x),dxdt)
    return I,x
end

"""
intsing(f,h,delta)

Integrate function `f` (possibly singular at 1 and -1) over
[-1+`delta`,1-`delta`] using discretization size `h`. Return
integral and a vector of the nodes used.
"""
function intsing(f,h,delta)
    # Find where to truncate the trapezoid sum.
    K = ceil(log(-2/pi*log(delta/2))/h)

    # Integrate over a transformed variable.
    t = h*(-K:K)
    x = @. tanh(pi/2*sinh(t))
    dxdt = @. pi/2*cosh(t) / (cosh(pi/2*sinh(t))^2)

    I = h*dot(f.(x),dxdt)
    return I,x
end

max_stars_repo_name: tobydriscoll/fnc

max_issues_repo_name: tobydriscoll/fnc

max_forks_repo_name: tobydriscoll/fnc


text:
@testset "165.compare-version-numbers.jl" begin
    @test compare_version("1.01", "1.001") == 0
    @test compare_version("1.0", "1.0.0") == 0
    @test compare_version("0.1", "1.1") == -1
    @test compare_version("1.0.1", "1") == 1
    @test compare_version("7.5.2.4", "7.5.3") == -1
end


text:
# ******************************************************************************************
# Notices:
#
# Copyright © 2022 United States Government as represented by the Administrator of the
# National Aeronautics and Space Administration.  All Rights Reserved.
#
# Disclaimers
#
# No Warranty: THE SUBJECT SOFTWARE IS PROVIDED "AS IS" WITHOUT ANY WARRANTY OF ANY KIND,
# EITHER EXPRESSED, IMPLIED, OR STATUTORY, INCLUDING, BUT NOT LIMITED TO, ANY WARRANTY
# THAT THE SUBJECT SOFTWARE WILL CONFORM TO SPECIFICATIONS, ANY IMPLIED WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, OR FREEDOM FROM INFRINGEMENT, ANY
# WARRANTY THAT THE SUBJECT SOFTWARE WILL BE ERROR FREE, OR ANY WARRANTY THAT DOCUMENTATION,
# IF PROVIDED, WILL CONFORM TO THE SUBJECT SOFTWARE. THIS AGREEMENT DOES NOT, IN ANY MANNER,
# CONSTITUTE AN ENDORSEMENT BY GOVERNMENT AGENCY OR ANY PRIOR RECIPIENT OF ANY RESULTS,
# RESULTING DESIGNS, HARDWARE, SOFTWARE PRODUCTS OR ANY OTHER APPLICATIONS RESULTING FROM
# USE OF THE SUBJECT SOFTWARE.  FURTHER, GOVERNMENT AGENCY DISCLAIMS ALL WARRANTIES AND
# LIABILITIES REGARDING THIRD-PARTY SOFTWARE, IF PRESENT IN THE ORIGINAL SOFTWARE, AND
# DISTRIBUTES IT "AS IS."
#
# Waiver and Indemnity:  RECIPIENT AGREES TO WAIVE ANY AND ALL CLAIMS AGAINST THE UNITED
# STATES GOVERNMENT, ITS CONTRACTORS AND SUBCONTRACTORS, AS WELL AS ANY PRIOR RECIPIENT.
# IF RECIPIENT'S USE OF THE SUBJECT SOFTWARE RESULTS IN ANY LIABILITIES, DEMANDS, DAMAGES,
# EXPENSES OR LOSSES ARISING FROM SUCH USE, INCLUDING ANY DAMAGES FROM PRODUCTS BASED ON,
# OR RESULTING FROM, RECIPIENT'S USE OF THE SUBJECT SOFTWARE, RECIPIENT SHALL INDEMNIFY AND
# HOLD HARMLESS THE UNITED STATES GOVERNMENT, ITS CONTRACTORS AND SUBCONTRACTORS, AS WELL
# AS ANY PRIOR RECIPIENT, TO THE EXTENT PERMITTED BY LAW.  RECIPIENT'S SOLE REMEDY FOR ANY
# SUCH MATTER SHALL BE THE IMMEDIATE, UNILATERAL TERMINATION OF THIS AGREEMENT.
# ******************************************************************************************

"""
AST core objective function abstract type. Defines fundamental reward function by combining
log probability, event bonus, and distance heuristic. Additional contributions to the
reward should be added by implementing `reward` function in GrayBox or BlackBox interface.
"""
abstract type AbstractCoreObjective end

"""
Default objective function. Sums components at each timestep.
"""
Base.@kwdef mutable struct WeightedObjective <: AbstractCoreObjective
    wl::Float64 = 1.0
    we::Float64 = 1.0
    wh::Float64 = 1.0
end

function (rf::WeightedObjective)(logprob::Float64, event::Float64, heuristic::Float64)
    return rf.wl * logprob + rf.we * event + rf.wh * heuristic
end

"""
Vector reward function. Maintains separate components to facilitate post-analysis and
enhanced learning methods.
"""
struct VectorObjective <: AbstractCoreObjective end

function (::VectorObjective)(logprob::Float64, event::Float64, heuristic::Float64)
    return (logprob, event, heuristic)
end

"""
Standard AST reward structure.
"""
Base.@kwdef mutable struct Reward <: AbstractReward
    marginalize::Bool                      = true
    heuristic::AbstractDistanceHeuristic   = GradientHeuristic()
    event_bonus::Float64                   = 0.0
    reward_function::AbstractCoreObjective = WeightedObjective()
end

"""
Internal handling of custom reward.
"""
custom_reward(sim::GrayBox, a::SampleAction) = Functoid(reward(sim, a.sample))
custom_reward(sim::BlackBox, ::SeedAction) = Functoid(reward(sim))


text:
struct PriorityFrame
    stream_identifier::UInt32
    exclusive::Bool
    dependent_stream_identifier::UInt32
    weight::UInt8
end

==(a::PriorityFrame, b::PriorityFrame) =
    a.stream_identifier == b.stream_identifier &&
    a.exclusive == b.exclusive &&
    a.dependent_stream_identifier == b.dependent_stream_identifier &&
    a.weight == b.weight

function decode_priority(header, payload)
    dependent_stream_identifier = (UInt32(payload[1]) << 24 + UInt32(payload[2]) << 16 +
                                   UInt32(payload[3]) << 8 + UInt32(payload[4])) & 0x7fffffff
    exclusive = payload[1] & 0x80 == 0x80
    weight = payload[5]

    return PriorityFrame(header.stream_identifier, exclusive, dependent_stream_identifier,
                         weight)
end

function encode_priority(frame)
    payload::Array{UInt8, 1} = [ UInt8(frame.dependent_stream_identifier >> 24);
                                 UInt8(frame.dependent_stream_identifier >> 16 & 0x000000ff);
                                 UInt8(frame.dependent_stream_identifier >> 8 & 0x000000ff);
                                 UInt8(frame.dependent_stream_identifier & 0x000000ff) ]
    payload[1] = (!isnothing(frame.exclusive) && frame.exclusive) ? (payload[1] | 0x80) : payload[1]
    push!(payload, frame.weight)

    return wrap_payload(payload, PRIORITY, 0x0, frame.stream_identifier)
end


text:
using MultiwayNumberPartitioning, HiGHS
using Test

@testset "MultiwayNumberPartitioning.jl" begin
    @testset "Example" begin
        include("../example/example.jl")
    end
end


text:
using LearnBase
using StatsBase
using SparseArrays
using Test

@testset "Observation" begin
    include("observation.jl")
end

@testset "Observation" begin
    include("iteration.jl")
end

include("aggmode.jl")
include("costs.jl")
include("other.jl")


text:
using FileIO
using Printf
include("julia2ims.jl")

"""
20201216: Package Spinning disk tif stack to imaris
"""

#src_dir = "/datahub/rawdata/tandeng/mRNA_imaging/mRNA_confocal_again/20201214/sCMOS/HE7-1NLS-uw-PWM_1/Default"
src_dir = "/home/hf/cf/20201216/sCMOS/HE7-1NLS-uw-PWM-3_1/Default"
out_dir = "/datahub/rawdata/tandeng/mRNA_imaging/mRNA_confocal_again/20201216"
#src_dir = ARGS[1]
t_len = 120
z_len = 20
pos_len = 40

for pos in 0:pos_len-1
#for pos in pos_len-1:pos_len-1
    println(pos)
    buffer = zeros(N0f16, 1900, 1300, z_len, t_len)
    Threads.@threads for t in 0:t_len-1
        id0 = t*z_len*pos_len + pos*z_len
        for z in 0:z_len-1
            id = @sprintf("img_channel000_position000_time%09d_z000.tif", id0+z)
            buffer[:, :, z+1, t+1] = load("$src_dir/$id");
        end
    end
    imsname = @sprintf("%s/pos%02d.ims", out_dir, pos)
    save2ims(reinterpret(UInt16, buffer), imsname, dx=0.108, compression=0)
end


text:
using JuMP, MathProgBase, MathOptInterface
using GLPK, GLPKMathProgInterface
using Cbc, Clp, Gurobi

struct KPDSupply{T<:Integer}
    item::String
    weight::T
    value::T
    quant::T
end
KPDSupply(itm::AbstractString, w::T, v::T, q::T=one(T)) where T<: Integer = KPDSupply(itm, w, v, q)
Base.show(io::IO, kdps::KPDSupply) = print(io, kdps.quant, " ", kdps.item, " ($(kdps.weight) kg, $(kdps.value) €)")

function solveIt(gear::Vector{KPDSupply{T}}, capacity::Integer, solver::DataType) where T<:Integer
    w = getfield.(gear, :weight)
    v = getfield.(gear, :value)
    q = getfield.(gear, :quant)
    sol = mixintprog(-v, w', '<', capacity, :Int, 0, q, solver())
    sol.status == :Optimal || error("this problem could not be solved")

    if all(q .== 1) # simpler case
        return gear[sol.sol == 1.0]
    else
        pack = similar(gear, 0)
        s = round.(Int, sol.sol)
        for (i, g) in enumerate(gear)
            iszero(s[i]) && continue
            push!(pack, KPDSupply(g.item, g.weight, g.value, s[i]))
        end
        return pack
    end
end

gear = [KPDSupply("map", 9, 150, 1),
        KPDSupply("compass", 13, 35, 1),
        KPDSupply("water", 153, 200, 2),
        KPDSupply("sandwich", 50, 60, 2),
        KPDSupply("glucose", 15, 60, 2),
        KPDSupply("tin", 68, 45, 3),
        KPDSupply("banana", 27, 60, 3),
        KPDSupply("apple", 39, 40, 3),
        KPDSupply("cheese", 23, 30, 1),
        KPDSupply("beer", 52, 10, 3),
        KPDSupply("suntan cream", 11, 70, 1),
        KPDSupply("camera", 32, 30, 1),
        KPDSupply("T-shirt", 24, 15, 2),
        KPDSupply("trousers", 48, 10, 2),
        KPDSupply("umbrella", 73, 40, 1),
        KPDSupply("waterproof trousers", 42, 70, 1),
        KPDSupply("waterproof overclothes", 43, 75, 1),
        KPDSupply("note-case", 22, 80, 1),
        KPDSupply("sunglasses", 7, 20, 1),
        KPDSupply("towel", 18, 12, 2),
        KPDSupply("socks", 4, 50, 1),
        KPDSupply("book", 30, 10, 2)]
mySolver = CbcSolver
pack = solveIt(gear, 400, mySolver)
println("The hiker should pack: \n - ", join(pack, "\n - "))
println("\nPacked weight: ", sum(getfield.(pack, :weight)), " kg")
println("Packed value: ", sum(getfield.(pack, :value)), " €")

max_stars_repo_name: rpb/Pivotal-Optimization

max_issues_repo_name: rpb/Pivotal-Optimization

max_forks_repo_name: rpb/Pivotal-Optimization


text:
using Photometry.Aperture: bounds

@testset "Apertures" begin
    ap_circ = CircularAperture(50, 40, 10)
    @test bounds(ap_circ) == (40, 60, 30, 50)
    @test size(ap_circ) == map(length, axes(ap_circ)) == (21, 21)
    @test size(ap_circ, 1) == 21
    @test all(axes(ap_circ) .== (30:50, 40:60))
    @test eachindex(ap_circ) == CartesianIndex(30, 40):CartesianIndex(50, 60)

    @test CircularAperture([50, 40], 10) == ap_circ

    ap_ann = CircularAnnulus(50, 40, 5, 10)
    @test bounds(ap_ann) == (40, 60, 30, 50)
    @test size(ap_ann) == (21, 21)
    @test all(axes(ap_ann) .== (30:50, 40:60))
    @test eachindex(ap_ann) == CartesianIndex(30, 40):CartesianIndex(50, 60)

    @test CircularAnnulus([50, 40], 5, 10) == ap_ann
end

@testset "broadcasting" begin
    ap = CircularAperture(3, 3, 2.5)
    data = ones(5, 7)
    weighted = ap .* data
    @test weighted == data .* ap # commutative
    @test sum(weighted) == sum(ap) == photometry(ap, data).aperture_sum
    @test all(iszero, weighted[:, 6:7])

end

@testset "Circle Aperture" begin
    c0 = CircularAperture(0, 0, 0)

    @test sprint(show, c0) == "CircularAperture(0, 0, r=0)"

    c1 = CircularAperture(0, 0, 1)

    @test sprint(show, c1) == "CircularAperture(0, 0, r=1)"

    sub_c1 = Subpixel(c1, 5)
    @test sprint(show, sub_c1) == "Subpixel(CircularAperture(0, 0, r=1), 5)"
end

@testset "Circle Annulus" begin
    c0 = CircularAnnulus(0, 0, 0, 0)

    @test sprint(show, c0) == "CircularAnnulus(0, 0, r_in=0, r_out=0)"

    c1 = CircularAnnulus(0, 0, 0, 1)

    @test sprint(show, c1) == "CircularAnnulus(0, 0, r_in=0, r_out=1)"
end


text:
"""
    read_hurdat(path,flag)
Get a hurdat url or filepath and downloads dataset from said resource. It returns the content of the downloaded csv as an array of arrays
"""

function read_hurdat(path,flag)
    @info "Started downloading and parsing HURDAT data"
    if flag
        res = get(path)
        res = decode(res.body,"UTF-8")
        res = split(res,"\n")
        res = replace.(res," " => "")
        res = split.(res,",")
        #@info res
    else
    end
    @info "Finished downloading and parsing HURDAT data"
    return res
end

function hurdat_to_dictionary!(content,data_dict,override_basin=false)
    @info "Started processing HURDAT2 data"
    is_header(line) = first(line[1]) ∈ ('A','C','E')
    current_id = "" #This will be re_defined in the header
    basin = ""
    for line in content
        if length(line) < 2
            continue
        end
        if is_header(line)
            #@info line
            basin = "north_atlantic"
            if first(line[1]) == "C"
                basin = "east_pacific"
            elseif first(line[1]) == "E"
                basin = "east_pacific"
            end
            if override_basin == true
                basin = "all"
            end
            data_dict[line[1]] = Dict(
                "id"=>line[1],
                "operational_id"=>"",
                "name"=>line[2],
                "year"=>parse(Int,line[1][nextind(line[1],0,5):end]),
                "season"=>parse(Int,line[1][nextind(line[1],0,5):end]),
                "basin"=>basin,
                "source_info"=>"NHC Hurricane Database",
                "source" => "hurdat"
            )
            current_id = line[1]
            data_dict[current_id]["date"]      = Union{Missing,DateTime}[]
            data_dict[current_id]["extra_obs"] = Union{Missing,Bool}[]
            data_dict[current_id]["special"]   = Union{Missing,String}[]
            data_dict[current_id]["type"]      = Union{Missing,String}[]
            data_dict[current_id]["lat"]       = Union{Missing,Float64}[]
            data_dict[current_id]["lon"]       = Union{Missing,Float64}[]
            data_dict[current_id]["vmax"]      = Union{Missing,Int}[]
            data_dict[current_id]["mslp"]      = Union{Missing,Int}[]
            data_dict[current_id]["wmo_basin"] = Union{Missing,String}[]
            data_dict[current_id]["ace"]       = 0.0
        else

            yyyymmdd,HHMM,special,storm_type,lat,lon,vmax,mslp = line[begin:9]

            date = DateTime(string(yyyymmdd,HHMM),"yyyymmddHHMM")

            if occursin('N',lat)
                lat = parse(Float64,split(lat,"N")[1])
            elseif occursin('S',lat)
                lat = -1.0*parse(Float64,split(lat,"S")[1])
            end
            if occursin('W',lon)
                lon = -1.0*parse(Float64,split(lon,"W")[1])
            elseif occursin('E',lon)
                lon = parse(Float64,split(lon,"E")[1])
            end
            # Handle missing data
            vmax = parse(Int,vmax) >= 0 ? parse(Int,vmax) : missing
            mslp = parse(Int,mslp) >= 800 ? parse(Int,mslp) : missing
            #Handel off-hour obs
            push!(data_dict[current_id]["extra_obs"], HHMM ∉ ["0000","0600","1200","1800"])
            #Fix storm type for cross-dateline storms
            storm_type = replace(storm_type,"ST" => "HU")
            storm_type = replace(storm_type,"TY" => "HU")

            #Add data to data_dict

            push!(data_dict[current_id]["date"],date)
            push!(data_dict[current_id]["special"],special)
            push!(data_dict[current_id]["type"],storm_type)
            push!(data_dict[current_id]["lat"],lat)
            push!(data_dict[current_id]["lon"],lon)
            push!(data_dict[current_id]["vmax"],vmax)
            push!(data_dict[current_id]["mslp"],mslp)

            if basin == "north_atlantic"
                push!(data_dict[current_id]["wmo_basin"],"north_atlantic")
            elseif basin == "east_pacific" & lon > 0.0
                push!(data_dict[current_id]["wmo_basin"],"west_pacific")
            elseif basin == "east_pacific" & lon <= 0.0
                push!(data_dict[current_id]["wmo_basin"],"east_pacific")
            else
                push!(data_dict[current_id]["wmo_basin"],"west_pacific")
            end

            #Calculate ACE & append to storm total
            if !ismissing(vmax)
                ace = (10 ^ (-4)) * (vmax ^ 2)
                if (HHMM ∈ ["0000","0600","1200","1800"]) & (storm_type ∈ ["SS","TS","HU"])
                    data_dict[current_id]["ace"] += round(ace; digits=4)
                end
            end

        end
    end

    #Account for operationally unnamed storms
    current_year = 0
    current_year_id = 1

    for storm in keys(data_dict)
        storm_data = data_dict[storm]
        storm_name = storm_data["name"]
        storm_year = storm_data["year"]
        storm_vmax = storm_data["vmax"]
        storm_id = storm_data["id"]

        max_wnd = all(ismissing.(storm_vmax)) ? missing : maximum(skipmissing(storm_vmax))
                #Fix current year
        if current_year == 0
            current_year = storm_year
        else
            if storm_year != current_year
                current_year = storm_year
                current_year_id = 1
                #special fix for 1992 in the Atlantic
                if (current_year == 1992) & (data_dict[current_id]["basin"] == "north_atlantic")
                    current_year_id = 2
                end

            end
        end

        #Estimate operational storm ID (which sometimes differs from HURDAT2 ID) ## TODO check if type can be constrained
        blocked_list = Union{Missing,String}[]
        potential_tcs = Union{Missing,String}["AL102017"]
        increment_but_pass = Union{Missing,String}[]

        if storm_name == "UNNAMED" &&  !ismissing(max_wnd) && (max_wnd >= 34) && (storm_id ∉ blocked_list)
            if storm_id ∈ increment_but_pass
                current_year_id += 1
            end
            continue
        elseif storm_id[begin:nextind(storm_id,1)] == "CP"
            continue
        else
            #Skip potential TCs
            if string(storm_id[begin:nextind(storm_id,1)],current_year_id,storm_year) ∈ potential_tcs
                current_year_id += 1
            end
            data_dict[storm]["operational_id"] = string(storm_id[begin:nextind(storm_id,1)],current_year_id,storm_year)
            current_year_id += 1
        end

        #Swap operational storm IDs, if necessary
        swap_list = ["EP101994","EP111994"]
        swap_pair = ["EP111994","EP101994"]

        if data_dict[storm]["operational_id"] ∈ swap_list
            @info swapping
            swap_idx = findfirst(isequal(data_dict[storm]["operational_id"]),swap_list)
            data_dict[key]["operational_id"] = swap_pair[swap_idx]
        end
    end

    @info "Completed processing HURDAT2 data"

    return data_dict
end

function read_and_process_hurdat(atlantic_url ="https://www.nhc.noaa.gov/data/hurdat/hurdat2-1851-2019-042820.txt",flag = true,override_basin=false)
    data_dict = Dict{String,Dict{T,N} where T where N}()
    res = Tropical.read_hurdat(atlantic_url,flag)
    hurdat_data = hurdat_to_dictionary!(res,data_dict,override_basin)
end


text:
# developed with Julia 1.1.1
#
# tests for EMSx package

using EMSx
using CSV
using Distributed
using Test

const test_directory = @__DIR__

const package_directory = joinpath(test_directory, "..")
const test_data_directory = joinpath(test_directory, "data")
const test_metadata_directory = joinpath(test_directory, "metadata")

isdir(test_data_directory) && rm(test_data_directory, recursive=true)
mkpath(test_data_directory)

const test_data_test_directory = joinpath(test_data_directory,"test")
const test_data_train_directory = joinpath(test_data_directory,"train")
const test_data_save_directory = joinpath(test_data_directory,"save")

mkpath(test_data_save_directory)

const test_prices_path = joinpath(test_metadata_directory, "edf_prices.csv")
const test_metadata_path = joinpath(test_metadata_directory, "metadata.csv")
const test_metadata2_path = joinpath(test_metadata_directory, "metadata2.csv")
const test_periods_path = joinpath(test_metadata_directory, "test_periods.csv")

@testset "EMSx.jl test set" begin

    @testset "Testing downloading data from Schneider" begin

        if haskey(ENV, "SCHNEIDER_API_KEY")

            @test isnothing(EMSx.download_sites_data(test_data_directory,
                                                      69:69;
                                                      periods = [1,2],
                                                      progress = false,
                                                      source = :schneider,
                                                      compressed = true))
            @test isfile(joinpath(test_data_directory, "69.csv.gz"))

            @test isnothing(EMSx.initialize_data(test_data_directory,
                                                  test_periods_path,
                                                  delete_files=false,
                                                  compressed=true))
            @test isfile(joinpath(test_data_directory, "69.csv.gz"))
            @test isfile(joinpath(test_data_test_directory, "69.csv.gz"))
            @test isfile(joinpath(test_data_train_directory, "69.csv.gz"))
            @test isnothing(EMSx.initialize_data(test_data_directory,
                                                  test_periods_path,
                                                  delete_files=true,
                                                  compressed=false))
            @test !isfile(joinpath(test_data_directory, "69.csv.gz"))
            @test isfile(joinpath(test_data_test_directory, "69.csv"))
            @test isfile(joinpath(test_data_train_directory, "69.csv"))

        else
            @warn "SCHNEIDER_API_KEY missing, Schneider download functions not tested"
        end

    end

    @testset "Testing downloading data from Zenodo" begin

        @test isnothing(EMSx.download_sites_data(test_data_directory,
                                                    8;
                                                    progress = false,
                                                    source = :zenodo,
                                                    compressed = true))
        @test isfile(joinpath(test_data_directory, "8.csv.gz"))

        # download pv

        @test isnothing(EMSx.download_pv(test_data_directory))
        @test isfile(joinpath(test_data_directory, "pv.csv.gz"))
        rm(joinpath(test_data_directory, "pv.csv.gz"))

        # non-compressed download from zenodo not available

        @test isnothing(EMSx.initialize_data(test_data_directory,
                                                test_periods_path,
                                                delete_files=true,
                                                compressed=true))
        @test !isfile(joinpath(test_data_directory, "8.csv.gz"))
        @test isfile(joinpath(test_data_test_directory, "8.csv.gz"))
        @test isfile(joinpath(test_data_train_directory, "8.csv.gz"))
        rm(joinpath(test_data_test_directory, "8.csv.gz"))
        rm(joinpath(test_data_train_directory, "8.csv.gz"))
    end
    @testset "EMS simulator's body" begin

        controller = EMSx.DummyController()
        prices = EMSx.load_prices(test_prices_path)

        site = EMSx.load_sites(test_metadata_path,
                                test_data_test_directory,
                                nothing,
                                test_data_save_directory)[1]

        period = EMSx.Period("1",
                                EMSx.read_site_file(site.path_to_test_data_csv),
                                site)

        @testset "Simulation" begin

            net_demand = period.data[98, :actual_consumption] - period.data[98, :actual_pv]
            @test EMSx.apply_control(1, 672, prices, period, 0., 0.) == (net_demand*prices.buy[1], 0.)
            simulation =  EMSx.simulate_period(controller, period, prices)
            @test simulation.result.soc == zeros(672)

            @test EMSx.simulate_site(controller, site, prices) === nothing
            @test EMSx.simulate_sites(controller,
                                        test_data_save_directory,
                                        test_prices_path,
                                        test_metadata_path,
                                        test_data_test_directory,
                                        nothing) === nothing
        end

        @testset "Evaluation" begin

            average_costs = EMSx.average_cost_per_site(joinpath(test_data_save_directory, "score.jld2"))
            performance_metrics = EMSx.evaluate_model(joinpath(test_data_save_directory, "score.jld2"))

            @test average_costs[!, :dummy][1] ≈ 7444.801871465684
            @test performance_metrics[!, :gain][1] ≈ -1940.8319533558906

        end

        @testset "Parallel computing" begin

            @test EMSx.init_parallel(2) == workers()
            @test nworkers() == 2
            @test isnothing(@everywhere EMSx.DIR)

            sites = 1:3
            @test isnothing(EMSx.download_sites_data_parallel(test_data_directory,
                                                                sites;
                                                                max_threads = 2,
                                                                progress = true,
                                                                compressed = true))
            @test all(isfile, joinpath.(test_data_directory, string.(sites) .* ".csv.gz"))

            @test isnothing(EMSx.initialize_data_parallel(test_data_directory,
                                                            test_periods_path;
                                                            progress = true))
            @test all(isfile, joinpath.([test_data_test_directory test_data_train_directory],
                                        string.(sites) .* ".csv.gz"))

            @test EMSx.simulate_sites_parallel(controller,
                                                test_data_save_directory,
                                                test_prices_path,
                                                test_metadata2_path,
                                                test_data_test_directory,
                                                nothing) === nothing

            average_costs = EMSx.average_cost_per_site(joinpath(test_data_save_directory, "score.jld2"))
            @test length(average_costs[!, :dummy]) == 3

            @test EMSx.init_parallel(1) == workers() == [1]

        end

    end

end

rm(test_data_directory, recursive=true)