text:
# Tutorial 2.4.1. Band structure calculations
# ===========================================
#
# Physics background
# ------------------
#  band structure of a simple quantum wire in tight-binding approximation
#
# Kwant features highlighted
# --------------------------
#  - Computing the band structure of a finalized lead.

function make_lead(a=1, t=1.0, W=10)
    # Start with an empty lead with a single square lattice
    lat = kwant.lattice.square(a)

    sym_lead = kwant.TranslationalSymmetry((-a, 0))
    lead = kwant.Builder(sym_lead)

    # build up one unit cell of the lead, and add the hoppings
    # to the next unit cell
    for j in range(0,length=W)
        lead[lat(0, j)] = 4 * t

        if j > 0
            lead[lat(0, j), lat(0, j - 1)] = -t
        end
        lead[lat(1, j), lat(0, j)] = -t
    end

    return lead
end


text:
#utility function for checking checksums of materials...
function check_jld2(material::MaterialModels.AbstractMaterial, loading::Vector, filename::String;
                        OVERWRITE_JLD2::Bool = false, debug_print::Bool = false)

    state = initial_material_state(material)
    stresses = []; tangents = [];
    for (i,load) in enumerate(loading)
        stress, tangent, state = material_response(material, load, state)

        #For user to check that the material acutally enters plastic state
        # and is not only testing the elastic part...
        if debug_print
            println("Iteration: $i")
            for symbol in fieldnames(typeof(state))
                println("$symbol: $(getproperty(state, symbol))")
            end
        end
        push!(stresses, stress)
        push!(tangents, tangent)
    end

    jld2_file = joinpath(@__DIR__, "jld2_files", string(filename, ".jld2"))
    if OVERWRITE_JLD2
        jldsave(jld2_file; stresses_jld2=stresses, tangents_jld2=tangents)
    else
        stresses_jld2, tangents_jld2 = load(jld2_file, "stresses_jld2", "tangents_jld2")
        @test stresses ≈ stresses_jld2
        @test tangents ≈ tangents_jld2
    end

end

function check_tangents_AD(material::MaterialModels.AbstractMaterial,loading::Vector)
    state = initial_material_state(material)
    for load in loading
        stress, tangent, state = material_response(material, load, state)
        ∂²Ψ∂C², ∂Ψ∂C, _ =  hessian(C -> MaterialModels.ψ(material, C), load, :all)
        stress_AD = 2*∂Ψ∂C
        tangent_AD = 4*∂²Ψ∂C²
        @test stress ≈ stress_AD
        @test tangent ≈ tangent_AD
    end
end

ext: jl


text:
em = ErgodicManagerR2("double gaussian", K=5, bins=100)

x0 = [0.5,0.01,pi/4]
N = 40
h = 0.1

tm = TrajectoryManager(x0, h, N, ConstantInitializer([0.0000]))

# things needed for dynamics
tm.dynamics = DubinsDynamics(0.3,0.1)
tm.Qn = eye(3)
tm.R = 0.01 * eye(1)
tm.Rn = 1 * eye(1)

xd, ud = pto_trajectory(em, tm)

tm.barrier_cost = 1
xd, ud = pto_trajectory(em, tm)


text:
using HDF5
using JSON
savepath = "SARA/NatCom2020/outer/data/"
xrd_file = "Bi2O3_19F44_01_outer_xrd_gradients_input_noise_iSARA.h5"

# savefile = "Bi2O3_19F44_01_outer_xrd_gradients_no_input_noise.h5"
xrd_f = h5open(savepath * xrd_file, "r")
xrd_temperatures = read(xrd_f, "temperatures")
xrd_dwelltimes = read(xrd_f, "dwelltimes")
xrd_coefficients = read(xrd_f, "coefficients")
close(xrd_f)

xrd_f = h5open(savepath * "Bi2O3_19F44_01_outer_xrd_data.h5", "r")
xrd_positions = read(xrd_f, "positions")
close(xrd_f)

include("../inner/inner_load_data.jl")
path = "/Users/sebastianament/Documents/SEA/XRD Analysis/SARA/Bi2O3_19F44_01/"
# path = ".."
file = "Bi2O3_19F44_01_inner_loop_data.json"
f = JSON.parsefile(path * file)
optical_positions, optical_coefficients, rescaling_parameters, optical_temperatures, optical_dwelltimes = load_data(f)

using Plots
plotly()
for i in 1:8
    p = 10
    match = (T, τ) -> xrd_temperatures[end, i] ≈ T && xrd_dwelltimes[end, i] ≈ τ
    j = 1
    while true
        if match(optical_temperatures[j]-p, optical_dwelltimes[j])
            break
        end
        j += 1
    end
    # j = findfirst(match, [optical_temperatures, optical_dwelltimes])
    plot(yscale = :log10)
    # xrd = xrd_coefficients[1, :, i]
    # opt = optical_coefficients[j][1]

    # xrd = abs.(diff(xrd))
    # opt = abs.(diff(opt))
    xrd = sum(abs2, diff(xrd_coefficients[:, :, i], dims = 2), dims = 1)
    opt = [optical_coefficients[j][k1][k2] for k1 in 1:16, k2 in 1:151]
    opt = sum(abs2, diff(opt, dims = 2), dims = 1)
    xrd = vec(xrd)
    opt = vec(opt)
    xrd ./= maximum(xrd)
    opt ./= maximum(opt)
    xrd_pos = xrd_positions[:, i]
    opt_pos = optical_positions[j]
    xrd_pos = xrd_pos[1:end-1]
    opt_pos = opt_pos[1:end-1]
    plot!(xrd_pos, xrd, label = "xrd")
    plot!(opt_pos, opt, label = "opt")
    gui()
end


text:
module Mumbling
    export accum

    function accum(s::AbstractString)::String
        accu = IOBuffer()
        for (index, chr) in enumerate(s)
          if index > 1
            write(accu,'-')
          end
          write(accu,titlecase(chr^index))
        end
        return String(take!(accu))
    end

end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
module VCFTool

include("filter_vcf_table.jl")

include("get_genotype.jl")

include("print_and_run_cmd.jl")

include("make_vcf_dataframe.jl")

include("make_vcf_indexedtable.jl")

include("make_vcf_ndsparse.jl")

include("search_region_with_julia.jl")

include("tabix_regions_from_file.jl")

include("query_regions.jl")

end


text:
@testset "loglikelihoods.jl" begin
    @model function demo(xs, y)
        s ~ InverseGamma(2, 3)
        m ~ Normal(0, √s)
        for i in eachindex(xs)
            xs[i] ~ Normal(m, √s)
        end

        y ~ Normal(m, √s)
    end

    xs = randn(3)
    y = randn()
    model = demo(xs, y)
    chain = sample(model, MH(), MCMCThreads(), 100, 2)
    var_to_likelihoods = pointwise_loglikelihoods(model, chain)
    @test haskey(var_to_likelihoods, "xs[1]")
    @test haskey(var_to_likelihoods, "xs[2]")
    @test haskey(var_to_likelihoods, "xs[3]")
    @test haskey(var_to_likelihoods, "y")

    for chain_idx in MCMCChains.chains(chain)
        for (i, (s, m)) in enumerate(zip(chain[:, :s, chain_idx], chain[:, :m, chain_idx]))
            @test logpdf(Normal(m, √s), xs[1]) == var_to_likelihoods["xs[1]"][i, chain_idx]
            @test logpdf(Normal(m, √s), xs[2]) == var_to_likelihoods["xs[2]"][i, chain_idx]
            @test logpdf(Normal(m, √s), xs[3]) == var_to_likelihoods["xs[3]"][i, chain_idx]
            @test logpdf(Normal(m, √s), y) == var_to_likelihoods["y"][i, chain_idx]
        end
    end

    var_info = VarInfo(model)
    results = pointwise_loglikelihoods(model, var_info)
    var_to_likelihoods = Dict(string(vn) => ℓ for (vn, ℓ) in results)
    s, m = var_info[SampleFromPrior()]
    @test logpdf(Normal(m, √s), xs[1]) == var_to_likelihoods["xs[1]"]
    @test logpdf(Normal(m, √s), xs[2]) == var_to_likelihoods["xs[2]"]
    @test logpdf(Normal(m, √s), xs[3]) == var_to_likelihoods["xs[3]"]
    @test logpdf(Normal(m, √s), y) == var_to_likelihoods["y"]
end


text:
using Distributed
addprocs(1)
@everywhere using Revise
@everywhere using VersatileHDPMixtureModels

function generate_data(dim, groups_count,sample_count,var,α = 10, γ = 1)
    crf_prior = hdp_prior_crf_draws(sample_count,groups_count,α,γ)
    pts,labels = generate_grouped_gaussian_from_hdp_group_counts(crf_prior[2],dim,var)
    return pts, labels
end

function results_stats(pred_dict, gt_dict)
    avg_nmi = 0
    for i=1:length(pred_dict)
        nmi = mutualinfo(pred_dict[i],gt_dict[i])
        avg_nmi += nmi
    end
    return avg_nmi / length(pred_dict)
end

function run_and_compare(pts,labels,gdim,iters = 100)
     gprior, lprior = create_default_priors(gdim,0,:niw)
     model = hdp_fit(pts,10,1,gprior,iters)
     model_results = get_model_global_pred(model[1])
     return results_stats(labels,model_results)
 end

gdim = 3
pts,labels = generate_data(gdim,4,100,100.0)
run_and_compare(pts,labels,gdim)


text:
###############################
####  OGGM configuration  #####
###############################

@everywhere using ParallelDataTransfer

export oggm_config, init_gdirs, PARAMS, PATHS

"""
    oggm_config()

Retrieves the initial glacier geometry (bedrock + ice thickness) for a glacier.
"""
function oggm_config()
    @eval begin
    @everywhere begin
    @eval ODINN begin
    cfg.initialize() # initialize OGGM configuration

    global PATHS = PyDict(cfg."PATHS")  # OGGM PATHS
    PATHS["working_dir"] = joinpath(homedir(), "Python/OGGM_data")  # Choose own custom path for the OGGM data
    global PARAMS = PyDict(cfg."PARAMS")
    PARAMS["hydro_month_nh"]=1

    # Multiprocessing
    # PARAMS["prcp_scaling_factor"], PARAMS["ice_density"], PARAMS["continue_on_error"]
    PARAMS["use_multiprocessing"] = true # Let's use multiprocessing for OGGM
    end # @eval ODINN
    end # @everywhere
    end # @eval
end

"""
    init_gdirs(rgi_ids; force=false)

Initializes Glacier Directories using OGGM. Wrapper function calling `init_gdirs_scratch(rgi_ids)`.
"""
function init_gdirs(rgi_ids; force=false)
    # Try to retrieve glacier gdirs if they are available
    try
        if force
            gdirs = init_gdirs_scratch(rgi_ids)
        else
            gdirs = workflow.init_glacier_directories(rgi_ids)
        end
        return gdirs
    catch
        @warn "Cannot retrieve gdirs from disk."
        println("Generating gdirs from scratch.")
        global create_ref_dataset = true # we force the creation of the reference dataset
        # Generate all gdirs if needed
        gdirs = init_gdirs_scratch(rgi_ids)
        return gdirs
    end
end

"""
    init_gdirs_scratch(rgi_ids)

Initializes Glacier Directories from scratch using OGGM.
"""
function init_gdirs_scratch(rgi_ids)
    # Check if some of the gdirs is missing files
    gdirs = workflow.init_glacier_directories(rgi_ids, from_prepro_level=3, prepro_border=40)
    list_talks = [
        # tasks.glacier_masks,
        # tasks.compute_centerlines,
        # tasks.initialize_flowlines,
        # tasks.compute_downstream_line,
        tasks.catchment_area,
        tasks.gridded_attributes,
        tasks.gridded_mb_attributes,
        # tasks.prepare_for_inversion,  # This is a preprocessing task
        # tasks.mass_conservation_inversion,  # This gdirsdoes the actual job
        # tasks.filter_inversion_output,  # This smoothes the thicknesses at the tongue a little
        # tasks.distribute_thickness_per_altitude,
        bedtopo.add_consensus_thickness   # Use consensus ice thicknesses from Farinotti et al. (2019)
    ]
    for task in list_talks
        # The order matters!
        workflow.execute_entity_task(task, gdirs)
    end
    return gdirs
end


text:
# You'll need to run julia from the directory that this
# file and MyModule.jl both reside in, for this to work:

push!(LOAD_PATH, ".")

using MyModule

MyModule.gleep()
println(code_lowered(MyModule.foop))

max_stars_repo_name: catseye/Dipple

max_issues_repo_name: catseye/Dipple

max_forks_repo_name: catseye/Dipple


text:
using Manifolds, Manopt, Plots, Test

@testset "Test Gradient checks" begin
    M = Sphere(10)
    q = zeros(11)
    q[1] = 1.0
    p = zeros(11)
    p[1:4] .= 1 / sqrt(4)
    r = log(M, p, q)

    F(M, p) = 1 / 2 * distance(M, p, q)^2
    gradF(M, p) = -log(M, p, q)

    @test check_gradient(M, F, gradF, p, r)

    gradF2(M, p) = -0.5 * log(M, p, q)
    @test_throws ErrorException check_gradient(M, F, gradF2, p, r; error=true)
    @test !check_gradient(M, F, gradF2, p, r)

    check_gradient(M, F, gradF, p, r; plot=true)

    #test windowsize error
    @test_throws ErrorException Manopt.find_best_slope_window(zeros(2), zeros(2), 20)
    @test_throws ErrorException Manopt.find_best_slope_window(zeros(2), zeros(2), [2, 20])
end


text:
##Saturation Tests
using Distributions, Random, DelimitedFiles, LinearAlgebra
include("KPNormal.jl")
using .KP

mutable struct DefaultExp
	cs
	thetas
	vs
end

function sim!(o, muhat)
	muhat[:] = randn!(muhat) ./ sqrt.(o.vs) .+ o.thetas
end

#simulates random normal stuff, but via a dataset
#dat is assumed n, S
function sim!(o, muhat, dat)
	S = size(dat, 2)
	n = size(dat, 1)
	for i = 1:S
		dat[:, i] .= randn(n) ./ sqrt.(o.vs/S)  .+ o.thetas #notice scaling by S
	end
	muhat[:] = mean(dat, dims=2)
end
######################
####
#For increasing n, many simulations compare methods
#on metrics
#		#Value wrt theta
#		#Time to compute
#		#optimal value of tau0
#Key is that the cs are scaled so that
# the constraint is effectively
#      c^T x \leq min(alpha*n, 1) --> 1/ n max(1/alpha ,n) *  c^T x \leq 1
#Since opt value is O(1), we do not scale opt value by n for stability.
function test_harness(f, numRuns, o, n_grid, Gamma_min, Gamma_max, Gamma_step, alpha)
	n_max = maximum(n_grid)
	S = 10
	muhat = zeros(Float64, n_max)
	xs  = zeros(Float64, n_max)
	dat = zeros(Float64, n_max, S)
	lam_t = 0.

	#write a header
	writedlm(f,  ["Run" "n" "Method" "thetaVal" "time" "tau0"], ',')

	for iRun = 1:numRuns
		#generate the entire path up to n_max
		sim!(o, muhat, dat)

		for n in n_grid
			#Take Views on Everything to speed up garbage collection
			x_t = view(xs, 1:n)
			vs  = view(o.vs, 1:n)
			thetas = view(o.thetas, 1:n)
			muhat_t = view(muhat, 1:n)
			dat_t = view(dat, 1:n, 1:S)

			#rescale the c's so that they work properly n what follows.
			cs  = o.cs[1:n] * max(1/alpha, n)

			#Compute performance of each method
			#SAA
			t = @elapsed x_t[:] = x(cs, muhat_t)
			thetaval = dot(thetas, x_t)
			writedlm(f,  [iRun n "SAA" thetaval t 0.], ',')

			#fullInfo val
			t = @elapsed x_t[:] = x(cs, thetas)
			thetaval = dot(thetas, x_t)
			writedlm(f,  [iRun n "FullInfo" thetaval t 0.], ',')

			#Tau MLE
			t =
			  @elapsed tauMLE, x_t[:] = x_MLE(cs, muhat_t, vs)

			thetaval = dot(thetas, x_t)
			writedlm(f,  [iRun n "EB_MLE" thetaval t tauMLE], ',')

			#Tau MM
			t =
			  @elapsed tauMM, x_t[:] = x_MM(cs, muhat_t, vs)
			thetaval = dot(thetas, x_t)
			writedlm(f,  [iRun n "EB_MM" thetaval t tauMM], ',')

			#Oracle MSE
			t =
			  @elapsed x_t[:], tau_CV = x_OR_MSE(cs, muhat_t, thetas, vs)

			thetaval = dot(thetas, x_t)
			writedlm(f,  [iRun n "OR_MSE" thetaval t tau_CV], ',')

			#Sure MSE
			t =
			  @elapsed x_t[:], tau_CV = x_sure_MSE(cs, muhat_t, vs)

			thetaval = dot(thetas, x_t)
			writedlm(f,  [iRun n "SURE_MSE" thetaval t tau_CV], ',')

			#Dirac Stein
			t =
			  @elapsed x_t[:], vals, objs = x_stein_exact(cs, muhat_t, vs, thetas)

			thetaval = dot(thetas, x_t)
			writedlm(f,  [iRun n "DiracStein" thetaval t vals[argmax(objs)]], ',')

			#Box with the optimized rate, i.e. h_n = n^-1/6 and scaling, altKernel
			h = n^-.16666
			t =
			  @elapsed x_t[:], vals, objs = x_stein_box(cs, muhat_t, vs, h, tau_step = .05)

			thetaval = dot(thetas, x_t)
			writedlm(f,  [iRun n "BoxStein" thetaval t vals[argmax(objs)]], ',')

			#Oracle Value
			t =
			  @elapsed x_t[:], vals, objs = best_x_tau(cs, muhat_t, vs, thetas)

			thetaval = dot(thetas, x_t)
			writedlm(f,  [iRun n "OR" thetaval t vals[argmax(objs)]], ',')

			#Hold-Out validation
			t =
			  @elapsed x_t[:], tau_grid, objs = KP.x_kFoldCV(cs, vs/S, dat_t, 2; tau_grid=collect(0.:.05:5))

			thetaval = dot(thetas, x_t)
			writedlm(f,  [iRun n "HO_EB" thetaval t tau_grid[argmax(objs)]], ',')

			#5-Fold validation
			t =
			  @elapsed x_t[:], tau_grid, objs = KP.x_kFoldCV(cs, vs/S, dat_t, 5; tau_grid=collect(0.:.05:5))

			thetaval = dot(thetas, x_t)
			writedlm(f,  [iRun n "K5_EB" thetaval t tau_grid[argmax(objs)]], ',')

			#LOO validation
			t =
			  @elapsed x_t[:], tau_grid, objs = KP.x_kFoldCV(cs, vs/S, dat_t, S; tau_grid=collect(0.:.05:5))

			thetaval = dot(thetas, x_t)
			writedlm(f,  [iRun n "LOO_EB" thetaval t tau_grid[argmax(objs)]], ',')

			###Regularization Methods
			#Oracle Regularization
			t =
			  @elapsed x_t[:], Gamma_grid, objs = KP.x_l2reg_CV(cs, muhat_t, vs, thetas,
														Gamma_min=Gamma_min, Gamma_max=Gamma_max, Gamma_step=Gamma_step)

			Gammahat = Gamma_grid[argmax(objs)]
			thetaval = dot(thetas, x_t)
			writedlm(f,  [iRun n "OracleReg" thetaval t Gammahat], ',')

			# #Our Stein Approach to Regularization
			t =
			  @elapsed x_t[:], Gamma_grid, objs = KP.x_stein_reg(cs, muhat_t, vs,
										Gamma_min=Gamma_min, Gamma_max=Gamma_max, Gamma_step=Gamma_step)

			Gammahat = Gamma_grid[argmax(objs)]
			thetaval = dot(thetas, x_t)
			writedlm(f,  [iRun n "SteinReg" thetaval t Gammahat], ',')

			# #NEW RO method with new threshold
			thresh = sqrt(2*log(1/.1))
			t =
			  @elapsed x_t[:] = KP.x_robFW(cs, muhat_t, vs, thresh, TOL=1e-4)

			thetaval = dot(thetas, x_t)
			writedlm(f,  [iRun n "FWRO_Eps_.1" thetaval t thresh], ',')

			thresh = sqrt(2*log(1/.05))
			t =
			  @elapsed x_t[:] = KP.x_robFW(cs, muhat_t, vs, thresh, TOL=1e-4)

			thetaval = dot(thetas, x_t)
			writedlm(f,  [iRun n "FWRO_Eps_.05" thetaval t thresh], ',')

			thresh = sqrt(2*log(1/.01))
			t =
			  @elapsed x_t[:] = KP.x_robFW(cs, muhat_t, vs, thresh, TOL=1e-4)

			thetaval = dot(thetas, x_t)
			writedlm(f,  [iRun n "FWRO_Eps_.01" thetaval t thresh], ',')

			#Hold Out Validation choice
			t =
			  @elapsed x_t[:], Gamma_grid, objs = KP.x_l2reg_kFoldCV(cs, vs/S, dat_t, 2;
                    Gamma_step=Gamma_step, Gamma_min=Gamma_min, Gamma_max=Gamma_max)

			Gammahat = Gamma_grid[argmax(objs)]
			thetaval = dot(thetas, x_t)
			writedlm(f,  [iRun n "HO_Reg" thetaval t Gammahat], ',')

			# #5-Fold Validation
			t =
			  @elapsed x_t[:], Gamma_grid, objs = KP.x_l2reg_kFoldCV(cs, vs/S, dat_t, 5;
                    Gamma_step=Gamma_step, Gamma_min=Gamma_min, Gamma_max=Gamma_max)

			Gammahat = Gamma_grid[argmax(objs)]
			thetaval = dot(thetas, x_t)
			writedlm(f,  [iRun n "K5_Reg" thetaval t Gammahat], ',')

			#LOO Validation
			t =
			  @elapsed x_t[:], Gamma_grid, objs = KP.x_l2reg_kFoldCV(cs, vs/S, dat_t, S;
	                Gamma_step=Gamma_step, Gamma_min=Gamma_min, Gamma_max=Gamma_max)
			Gammahat = Gamma_grid[argmax(objs)]
			thetaval = dot(thetas, x_t)
			writedlm(f,  [iRun n "LOO_Reg" thetaval t Gammahat], ',')
		end
		flush(f)
	end
end

#######################
###  Reads in a theta/cs/vs specification
function test_saturation(file_out, numRuns, n_grid, seed, param_path;
						Gamma_min = 1., Gamma_max = 100., Gamma_step=.5, alpha=1.)
	Random.seed!(seed)
	n_max = maximum(n_grid)
	dat, header = readdlm(param_path, ',', header=true)

	#confirm that n_max works
	@assert n_max <= size(dat, 1) "Param file too short for n_max"
	cs = dat[1:n_max, 3]
	thetas = dat[1:n_max, 1]
	vs = dat[1:n_max, 2]

	o = DefaultExp(cs, thetas, vs)
	file_name = "$(file_out)_$(seed).csv"
	f = open(file_name, "w")
	test_harness(f, numRuns, o, n_grid, Gamma_min, Gamma_max, Gamma_step, alpha)
	close(f)
	return file_name
end

max_stars_repo_path: Normal Model/testSaturation_Harness.jl

max_issues_repo_path: Normal Model/testSaturation_Harness.jl

max_forks_repo_path: Normal Model/testSaturation_Harness.jl


text:
import Common.test_clone

function run_test_clone()
    test_clone(Molecule([Atom(index=5)]), Dict(:atoms => [Atom(index=6)]), Molecule(atoms=[Atom(index=6)]), "Molecule test_clone 1")
end


text:
# GRADIENT OPERATOR WITH DIMENSION 1/Δ

function GTx_nonperiodic!(dTdx::Matrix{Numtype},T::Matrix{Numtype})
    #= Calculates the gradient in x-direction on the T-grid.
    The result dTdx sits on the u-grid. =#

    dTdx[:,:] = one_over_Δ*(T[2:end,:] - T[1:end-1,:])
end

function GTx_periodic!(dTdx::Matrix{Numtype},T::Matrix{Numtype})
    #= Calculates the gradient in x-direction on the T-grid.
    The result dTdx sits on the u-grid. =#

    dTdx[1,:] = one_over_Δ*(T[1,:] - T[end,:])
    dTdx[2:end,:] = one_over_Δ*(T[2:end,:] - T[1:end-1,:])
end

function GTy!(dTdy::Matrix{Numtype},T::Matrix{Numtype})
    #= Calculates the gradient in y-direction on the T-grid.
    The result dTdy sits on the v-grid. =#

    dTdy[:,:] = one_over_Δ*(T[:,2:end] - T[:,1:end-1])
end

function Gux_nonperiodic!(dudx::Matrix{Numtype},u::Matrix{Numtype})
    #= Calculates the gradient in x-direction on the u-grid.
    The result dudx sits on the T-grid. =#

    dudx[1,:] = one_over_Δ*u[1,:]              # -0, kinematic bc
    dudx[2:end-1,:] = one_over_Δ*(u[2:end,:] - u[1:end-1,:])
    dudx[end,:] = (-one_over_Δ)*u[end,:]       # +0, kinematic bc
end

function Gux_periodic!(dudx::Matrix{Numtype},u::Matrix{Numtype})
    #= Calculates the gradient in x-direction on the u-grid.
    The result dudx sits on the T-grid. =#

    dudx[1:end-1,:] = one_over_Δ*(u[2:end,:]-u[1:end-1,:])
    dudx[end,:] = one_over_Δ*(u[1,:]-u[end,:])
end

function Gvy!(dvdy::Matrix{Numtype},v::Matrix{Numtype})
    #= Calculates the gradient in y-direction on the v-grid.
    The result dvdy sits on the T-grid. =#

    dvdy[:,2:end-1] = one_over_Δ*(v[:,2:end] - v[:,1:end-1])
    dvdy[:,1] = one_over_Δ*v[:,1]              # -0, kinematic bc
    dvdy[:,end] = (-one_over_Δ)*v[:,end]       # +0, kinematic bc
end

function Guy_nonperiodic!(dudy::Matrix{Numtype},u::Matrix{Numtype})
    #= Calculates the gradient in y-direction on the u-grid.
    The result dudy sits on the q-grid. =#

    dudy[2:end-1,2:end-1] = one_over_Δ*(u[:,2:end] - u[:,1:end-1])
    dudy[2:end-1,1] = α_over_Δ*u[:,1] #  α is the lateral boundary condition parameter
    dudy[2:end-1,end] = -α_over_Δ*u[:,end]
    dudy[1,:] .= zeero       #TODO probably redundant if initialised with zeros
    dudy[end,:] .= zeero
end

function Guy_periodic!(dudy::Matrix{Numtype},u::Matrix{Numtype})
    #= Calculates the gradient in y-direction on the u-grid.
    The result dudy sits on the q-grid. =#

    dudy[:,2:end-1] = one_over_Δ*(u[:,2:end] - u[:,1:end-1])
    dudy[:,1] = α_over_Δ*u[:,1] #  α is the lateral boundary condition parameter
    dudy[:,end] = -α_over_Δ*u[:,end]
end

function Gvx_nonperiodic!(dvdx::Matrix{Numtype},v::Matrix{Numtype})
    #= Calculates the gradient in x-direction on the v-grid.
    The result dvdx sits on the q-grid. =#

    dvdx[2:end-1,2:end-1] = one_over_Δ*(v[2:end,:] - v[1:end-1,:])
    dvdx[1,2:end-1] = α_over_Δ*v[1,:] #  α is the lateral boundary condition parameter
    dvdx[end,2:end-1] = -α_over_Δ*v[end,:]
    dvdx[:,1] .= zeero      #TODO redundant if initialised with zeros
    dvdx[:,end] .= zeero
end

function Gvx_periodic!(dvdx::Matrix{Numtype},v::Matrix{Numtype})
    #= Calculates the gradient in x-direction on the v-grid.
    The result dvdx sits on the q-grid. =#

    dvdx[2:end,2:end-1] = one_over_Δ*(v[2:end,:] - v[1:end-1,:])
    dvdx[1,2:end-1] = one_over_Δ*(v[1,:]-v[end,:])
    dvdx[:,1] .= zeero      #TODO redundant if initialised with zeros
    dvdx[:,end] .= zeero
end

function Gqx_nonperiodic!(dqdx::Matrix{Numtype},q::Matrix{Numtype})
    #= Calculates the gradient in x-direction on the q-grid.
    The result dqdx sits on the v-grid. =#

    dqdx[:,:] = one_over_Δ*(q[2:end,2:end-1] - q[1:end-1,2:end-1])
end

function Gqx_periodic!(dqdx::Matrix{Numtype},q::Matrix{Numtype})
    #= Calculates the gradient in x-direction on the q-grid.
    The result dqdx sits on the v-grid. =#

    dqdx[1:end-1,:] = one_over_Δ*(q[2:end,2:end-1] - q[1:end-1,2:end-1])
    dqdx[end,:] = one_over_Δ*(q[1,2:end-1] - q[end,2:end-1])
end

function Gqy_nonperiodic!(dqdy::Matrix{Numtype},q::Matrix{Numtype})
    #= Calculates the gradient in y-direction on the q-grid.
    The result dqdy sits on the u-grid. =#

    dqdy[:,:] = one_over_Δ*(q[2:end-1,2:end] - q[2:end-1,1:end-1])
end

function Gqy_periodic!(dqdy::Matrix{Numtype},q::Matrix{Numtype})
    #= Calculates the gradient in y-direction on the q-grid.
    The result dqdy sits on the u-grid. =#

    dqdy[:,:] = one_over_Δ*(q[:,2:end] - q[:,1:end-1])
end

max_stars_repo_name: eschnett/juls

max_issues_repo_name: eschnett/juls

max_forks_repo_name: eschnett/juls


text:
export filter

import Base: filter
import Base: show

"""
    filter(filterFn::F) where { F <: Function }

Creates a filter operator, which filters items of the source Observable by emitting only
those that satisfy a specified `filterFn` predicate.

# Producing

Stream of type `<: Subscribable{L}` where `L` refers to type of source stream

# Arguments
- `filterFn::F`: predicate function with `(data::T) -> Bool` signature

# Examples
```jldoctest
using Rocket

source = from([ 1, 2, 3, 4, 5, 6 ])
subscribe!(source |> filter((d) -> d % 2 == 0), logger())
;

# output

[LogActor] Data: 2
[LogActor] Data: 4
[LogActor] Data: 6
[LogActor] Completed
```

See also: [`AbstractOperator`](@ref), [`InferableOperator`](@ref), [`ProxyObservable`](@ref), [`logger`](@ref)
"""
filter(filterFn::F) where { F <: Function }  = FilterOperator{F}(filterFn)

struct FilterOperator{F} <: InferableOperator
    filterFn::F
end

function on_call!(::Type{L}, ::Type{L}, operator::FilterOperator{F}, source) where { L, F }
    return proxy(L, source, FilterProxy{F}(operator.filterFn))
end

operator_right(operator::FilterOperator, ::Type{L}) where L = L

struct FilterProxy{F} <: ActorProxy
    filterFn::F
end

actor_proxy!(::Type{L}, proxy::FilterProxy{F}, actor::A) where { L, A, F } = FilterActor{L, A, F}(proxy.filterFn, actor)

struct FilterActor{L, A, F} <: Actor{L}
    filterFn :: F
    actor    :: A
end

function on_next!(actor::FilterActor{L}, data::L) where L
    if actor.filterFn(data)
        next!(actor.actor, data)
    end
end

on_error!(actor::FilterActor, err) = error!(actor.actor, err)
on_complete!(actor::FilterActor)   = complete!(actor.actor)

Base.show(io::IO, ::FilterOperator)         = print(io, "FilterOperator()")
Base.show(io::IO, ::FilterProxy)            = print(io, "FilterProxy()")
Base.show(io::IO, ::FilterActor{L}) where L = print(io, "FilterActor($L)")


text:
using ChainRulesTestUtils
using FiniteDifferences
using Random
using TensorRules
using Test
using Zygote

import TensorRules: @fn∇, rhs_to_args, make_only_product

Random.seed!(42)

Zygote.refresh()

@testset "TensorRules" begin
    include("rhsparse.jl")
    include("rules.jl")
    include("zygote.jl")
end


text:
__precompile__()

module ConcreteInterfaces

export @interface

include("ComputedFieldTypes/src/ComputedFieldTypes.jl")
using .ComputedFieldTypes
using FunctionWrappers: FunctionWrapper
using MacroTools

abstract type AbstractInterface end

function argtype(expr)
    if @capture(expr, x_::T_)
        return T
    else
        return :(Any)
    end
end

function parse_wrapper_field(expr)
    @assert @capture(expr, f_(y__)::X_ = ex_)
    :($(f)::FunctionWrapper{$(esc(X)), $(esc(:Tuple)){$(esc.(argtype.(y))...)}})
end

baretype(s::Symbol) = s
function baretype(ex::Expr)
    @assert @capture(ex, X_ <: T_)
    return X
end

function build_constructor(name, T, captures)
    :($name{$(baretype.(T)...)}($(captures...)) where {$(T...)})
end

function parse_closure(expr)
    @assert @capture(expr, f_(args__)::T_ = body_)
    if isempty(args)
        :(() -> $(body))
    else
        :($(args...) -> $(body))
    end
end

function outer_method(name, typeparams, methodexpr)
    @assert @capture(methodexpr, f_(args__)::T_ = body_)
    argtype = if isempty(typeparams)
        name
    else
        :($(name){$(baretype.(typeparams)...)})
    end
    :($f(self::$(argtype), $(args...)) where {$(typeparams...)} = self.$f($(args...)))
end

function parse_name(constructor)
    if @capture(constructor, S_{T__}(args__))
        return S, T, args
    elseif @capture(constructor, S_(args__))
        return S, [], args
    else
        error("couldn't parse $(constructor)")
    end
end

add_any(x::Symbol) = :($x::Any)
add_any(x::Expr) = x

macro interface(constructor, block)
    method_args = [arg for arg in block.args if arg.head == :(=)]
    wrapper_fields = parse_wrapper_field.(method_args)
    name, typeparams, captures = parse_name(constructor)
    quote
        @computed struct $(esc(name)){$(esc.(typeparams)...)} <: AbstractInterface
            $(esc.(add_any.(captures))...)
            $(wrapper_fields...)

            $(esc(build_constructor(name, typeparams, captures))) = new($(esc.(captures)...), $(esc.(parse_closure.(method_args))...))
        end

        $(esc.(outer_method.(name, [typeparams], method_args))...)
    end
end

end


text:
bxd_geno = "http://genenetwork.org/api/v_pre1/genotypes/BXD.geno"
spleen_pheno="http://datafiles.genenetwork.org/download/GN283/GN283_MeanDataAnnotated_rev081815.txt"

download(bxd_geno, joinpath(@__DIR__,"..", "data", "rawdata", "bxd.geno"))
download(spleen_pheno, joinpath(@__DIR__,"..", "data", "rawdata", "bxdspleen.txt"))


text:
using Gnuplot
Θ = LinRange(0, 2π, 100) # 50
Φ = LinRange(0, π, 20)
r = 0.8
x = [r * cos(θ) * sin(ϕ)      for θ in Θ, ϕ in Φ]
y = [r * sin(θ) * sin(ϕ)      for θ in Θ, ϕ in Φ]
z = [r * cos(ϕ) for θ in Θ, ϕ in Φ]
@gsp "set pm3d depthorder" :-
@gsp :- x y z "w l" "set view equal xyz"
save(term="pngcairo size 600,600", output="sphere1.png") # hide

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
# Arvind Prasadan
# 2019 January
# DMF Package
# Evaluates the Autocorrelation Function

"""
	autocorrelation(x; circular = false, bias_correct = false)

Evaluates the Autocorrelation Function of a signal

# Arguments
- `x`: signal (vector)
- `circular`: Boolean whether to use the circular inner product or not
- `bias_correct`: Boolean whether to use n - k instead of n in normalization for a lag of k (theoretically undesirable); ignored if circular is true

# Outputs
- `gamma`: Autocorrelation function (lags 0 to n - 1)
"""

function autocorrelation(x; circular = false, bias_correct = false)

    x = x[:]
    n = length(x)

    x = x .- mean(x)

    gamma = zeros(n)
    if false == circular
        for k = 0:1:(n - 1)
            gamma[k + 1] = x[1:1:(n - k)]' * x[(1 + k):1:n]

            if true == bias_correct
                gamma[k + 1] = gamma[k + 1] / (n - k)
            end
        end

        if false == bias_correct
            gamma = gamma / n
        end
    elseif true == circular
        for k = 0:1:(n - 1)
            gamma[k + 1] = circshift(x, k)' * x
        end

        gamma = gamma / n
    end

    return gamma
end


text:
using BinaryBuilder

# Collection of sources required to build ERFABuilder
sources = [
    "https://github.com/liberfa/erfa/releases/download/v1.4.0/erfa-1.4.0.tar.gz" =>
    "035b7f0ad05c1191b8588191ba4b19ba0f31afa57ad561d33bd5417d9f23e460",
]

script = raw"""
cd $WORKSPACE/srcdir
cd erfa-1.4.0/
if [[ ${target} == i686-w64* ]] || [[ ${target} == x86_64-w64* ]]; then
    sed -i 's/LT_INIT/LT_INIT([win32-dll])/' configure.ac
    sed -i 's/liberfa_la_LDFLAGS = -version-info \$(VI_ALL)/liberfa_la_LDFLAGS = -no-undefined -version-info $(VI_ALL)/' src/Makefile.am
fi
autoreconf -fi
./configure --prefix=$prefix --host=$target
make -j${nproc}
make install
"""

# Build on all default platforms
platforms = supported_platforms()

# No dependencies
dependencies = []

products = prefix -> [
    LibraryProduct(prefix, "liberfa", :liberfa)
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, "liberfa", sources, script, platforms, products, dependencies)


text:
@testset "Utils" begin
    @testset "testSparseUtils.jl" begin
        include("testSparseUtils.jl")
    end
    @testset "testSortpermFast.jl" begin
        include("testSortpermFast.jl")
    end
    @testset "testUniqueIdx.jl" begin
        include("testUniqueIdx.jl")
    end
    @testset "testVariousUtils.jl" begin
        include("testVariousUtils.jl")
    end
    @testset "testTests.jl" begin
        include("testTests.jl")
    end
end


text:
using OptiMo, Bazinga
using CUTEst, NLPModelsIpopt
using CSV, Printf
using DataFrames, Query
using BenchmarkProfiles

problems = CUTEst.select( min_var=1, max_var=100, min_con=1, max_con=100 )
nprob = length( problems )

data_alpx = DataFrame()
data_ipopt = DataFrame()

TOL_OPTIM = 1e-4
TOL_CVIOL = TOL_OPTIM

alpx = Bazinga.ALPX(  verbose=false,
                      tol_optim=TOL_OPTIM,
                      tol_cviol=TOL_CVIOL,
                      subsolver=:zerofpr,
                      max_iter=50,
                      max_sub_iter=1000 )

problem = problems[2]
nlp = CUTEstModel( problem )
out = ipopt( nlp, tol=1e-3 )
prob = NLPOptiModel( nlp )
out = alpx( prob )
finalize( nlp )

for id in 1:nprob
    local problem, nlp, out, prob

    problem = problems[id]
    @printf( "********************  %s  ********************\n", problem )
    nlp = CUTEstModel( problem )

    ##########
    # IPOPT
    ##########
    out = ipopt( nlp, tol=TOL_OPTIM,print_level=0 )
    etime = out.elapsed_time
    iters = out.iter
    optim = out.dual_feas
    cviol = out.primal_feas
    solved = out.solver_specific[:internal_msg] == :Solve_Succeeded ? 1 : 0
    push!( data_ipopt, (id=id,
                  time=etime,
                  iter=iters,
                  optim=optim,
                  cviol=cviol,
                  solved=solved) )

    ##########
    # ALPX
    ##########
    prob = NLPOptiModel( nlp )
    out = alpx( prob )
    etime = out.time
    iters = out.iterations
    optim = out.optimality
    cviol = out.cviolation
    solved = out.status == :first_order ? 1 : 0
    push!( data_alpx, (id=id,
                  time=etime,
                  iter=iters,
                  optim=optim,
                  cviol=cviol,
                  solved=solved) )

    finalize( nlp )
end

n_tot = size( data_alpx, 1 )
@printf(" %d problems \n", n_tot)
data_solved_alpx = data_alpx |> @filter(_.solved == 1) |> DataFrame
n_solved_alpx = size( data_solved_alpx, 1 )
@printf(" ALPX     %4.1f/100 first order\n", 100*n_solved_alpx/n_tot)
data_solved_ipopt = data_ipopt |> @filter(_.solved == 1) |> DataFrame
n_solved_ipopt = size( data_solved_ipopt, 1 )
@printf(" IPOPT    %4.1f/100 first order\n", 100*n_solved_ipopt/n_tot)

# store results
inttol = -Int(log10(TOL_OPTIM))
foldername = "/home/alberto/Documents/"
filename = "cutest_" * "alpx" * "_$(inttol).csv"
CSV.write( foldername * "Bazinga.jl/test/data/" * filename, data_alpx, header=false )
filename = "cutest_" * "ipopt" * "_$(inttol).csv"
CSV.write( foldername * "Bazinga.jl/test/data/" * filename, data_ipopt, header=false )

# plot profiles
cost(data) = data.time + Inf .* (data.solved .== 0)
T = [cost(data_ipopt) cost(data_alpx)]
performance_profile(T, ["IPOPT","ALPX"],title="CUTEst, ϵ=$TOL_OPTIM")


text:
"""
```julia
parsedlm([T::Type], filepath::String, delimiter::Char)
```

Parse a delimited text file, given a `filepath` and `delimiter`, and return the
parsed contents as a `MallocMatrix{T}`, that is a 2D `MallocArray` containing
numbers of type `T`.

If not specified, the parse type `T` will default to `Float64`.

## Examples
```julia
julia> using StaticTools

julia> m = (1:10) * (1:10)';

julia> fp = fopen(c"testfile.tsv", c"w"); printf(fp, m); fclose(fp);

julia> parsedlm(Int32, c"testfile.tsv", '\t')
10×10 MallocMatrix{Int32}:
  1   2   3   4   5   6   7   8   9   10
  2   4   6   8  10  12  14  16  18   20
  3   6   9  12  15  18  21  24  27   30
  4   8  12  16  20  24  28  32  36   40
  5  10  15  20  25  30  35  40  45   50
  6  12  18  24  30  36  42  48  54   60
  7  14  21  28  35  42  49  56  63   70
  8  16  24  32  40  48  56  64  72   80
  9  18  27  36  45  54  63  72  81   90
 10  20  30  40  50  60  70  80  90  100
```
"""
@inline parsedlm(filepath, delimiter::Char) = parsedlm(Float64, filepath, delimiter::Char)
@inline function parsedlm(::Type{T}, filepath, delimiter::Char) where {T}

	# File to open
	fp = fopen(filepath, c"r")

	if fp == C_NULL
		error(c"File does not exist!")
		return MallocMatrix{T}(undef, 0, 0)
	end

	delim = delimiter % Int32
	chars, rows, columns = 0, 0, 0
	maxcolumns, maxchars = 0, 0

	# Determine maximum number of characters per row, delimiters per row, and rows
	@inbounds while (c = getc(fp)) > 0
		chars += 1
		c == delim && (columns += 1)
		if c == 10 #Int32('\n')
			rows += 1
			# If there is a trailing delimiter, don't add an extra column for it
			fseek(fp, -2)
			getc(fp) != delim && (columns += 1)
			fseek(fp, +1)
			# See if we have a new maximum, and reset the counters
			chars > maxchars && (maxchars = chars)
			columns > maxcolumns && (maxcolumns = columns)
			chars, columns = 0, 0
		end
	end
	# If the last line isn't blank, add one more to the row counter
	fseek(fp, -1)
	getc(fp) != Int32('\n') && (rows += 1)
	frewind(fp)

	# If file starts with \ufeff, ignore it (with warning)
	c1, c2 = getc(fp), getc(fp)
	if c1 == 0xfe && c2 == 0xff
		warn(c"Skipping hidden U+FEFF character at start of input file.\n")
	else
		frewind(fp)
	end

	# # if debug
	# printf(c"Maximum number of characters: %d\n", maxchars)
	# printf(c"Maximum number of delimiters: %d\n", maxcolumns)
	# printf(c"Number of rows: %d\n", rows)

	# Allocate space for the imported array
	importedMatrix = MallocMatrix{T}(undef, rows, maxcolumns)
	field = MallocVector{Int}(undef, maxcolumns+2)
	str = MallocString(undef, maxchars+2)

	# For each line,
	i = 1
	@inbounds while gets!(str,fp) != C_NULL

		# identify the delimited fields,
		field[1] = 0
		k, columns = 1, 1
		while str[k] != 0x00 #UInt8('\0')
			if str[k] == UInt8(delim)
				str[k] = '\0'
				columns += 1
				field[columns] = k
			elseif str[k] == UInt8('\n')
				str[k] = '\0'
			end
			k += 1
		end

		# and perform operations on each field
		for j = 1:maxcolumns
			importedMatrix[i,j] = parse(T, pointer(str) + field[j])
		end
		i += 1
	end
	fclose(fp)
	free(str)
	free(field)

	return importedMatrix
end


text:
@shorthands marginalscatter

@recipe function f(::Type{Val{:marginalscatter}}, plt::AbstractPlot; density = false)
    x, y = plotattributes[:x], plotattributes[:y]
    i = isfinite.(x) .& isfinite.(y)
    x, y = x[i], y[i]
    scale = get(plotattributes, :scale, :identity)
    xlims, ylims = Plots.widen(Plots.ignorenan_extrema(x)..., scale), Plots.widen(Plots.ignorenan_extrema(y)..., scale)

    # set up the subplots
    legend --> false
    link := :both
    grid --> false
    layout --> @layout [
        topscatter           _
        scatter2d{0.9w,0.9h} rightscatter
    ]

    # main scatter2d
    @series begin
        seriestype := :scatter
        right_margin --> 0mm
        top_margin --> 0mm
        subplot := 2
        xlims --> xlims
        ylims --> ylims
    end

    # these are common to both marginal scatter
    ticks := nothing
    xguide := ""
    yguide := ""
    fillcolor --> Plots.fg_color(plotattributes)
    linecolor --> Plots.fg_color(plotattributes)

    if density
        trim := true
        seriestype := :density
    else
        seriestype := :scatter
    end

    # upper scatter
    @series begin
        subplot := 1
        bottom_margin --> 0mm
        showaxis := :x
        x := x
        y := ones(y |> size)
        xlims --> xlims
        ylims --> (0.95, 1.05)
    end

    # right scatter
    @series begin
        orientation := :h
        showaxis := :y
        subplot := 3
        left_margin --> 0mm
        # bins := edges2
        y := y
        x := ones(x |> size)
    end
end

# # now you can plot like:
# marginalscatter(rand(1000), rand(1000))


text:
function midpoints(s,t)
    maxs, mins, exts = find_extrema(s)
    mpx = Float64[]
    mpy = Float64[]
    for i=1:(length(exts)-1)
        push!(mpx, (t[exts[i]] + t[exts[i+1]]) /2.0)
        push!(mpy, (s[exts[i]] + s[exts[i+1]]) /2.0)
    end

    mpx,mpy
end

function get_spline(s,t,sparam)
    mpt, mps = midpoints(s,t)
    spl = Spline1D(mpt,mps,bc="extrapolate",s=sparam)
    return spl(t)
end

function bsift(s,t;sparam=0.0)
    count = 1
    sd = 1.0
    hlast = copy(s)
    hnew = copy(s)
    while sd > 0.25
        spls = get_spline(hlast,t,sparam)
        @. hnew = hlast - spls
        sd = sum(abs2.(hlast-hnew))/sum(hlast.^2)
        @. hlast = hnew
        count +=1
    end
    return hnew
end

function bEMD(s,t;maximfs=5,sparam=0.0)
    sig = copy(s)
    N=length(s)
    imfs = ElasticArray{Float64}(undef,N,0)
    for i in 1:maximfs
        h = bsift(sig,t;sparam=sparam)
        append!(imfs, h)
        sig .-= h
    end
    return convert(Array,imfs)
end


text:
using CSV

SpinLookUp = Dict{Tuple{Int64,Int64,Int64},String}()

for row in CSV.File(joinpath(@__DIR__, "SpinWeightedSphericalHarmonicsTable.csv"))
    s = row.s; l=row.l; m=row.m; func = row.Func
    SpinLookUp[(s,l,m)] = row.Func
end

function sterlings(n)
    return sqrt(2*pi*n)*(n/exp(1))^n
end

function MakeSwitchCase(param,iterable,args)
    newstr = """if """
    for x ∈ iterable
       newstr *= """$(param) == $(x)
    $(getcode(args...,x))
elseif """
    end
    newstr[1:(end-7)]*" end"
end

function getcode(s,l,m)
    return "return begin $(SpinLookUp[(s,l,m)]) end"
end

getcode(s,l) = MakeSwitchCase("m",-l:l,(s,l))
getcode(s) = MakeSwitchCase("l",abs(s):20,(s,))
getcode() = MakeSwitchCase("s",-2:2,())

#print("again")
"""
function SpinWeightedSphericalCalculation(z,s,l,m)
    zm = 1-z; zp = 1+z;
    if l < max(abs(s),abs(m))
        return 0.0
    end
    #println("BEFORE WE HAVE",s," ", l," ", m)
    $(getcode())
    if ((l+m > 20) || (l+s > 20) || (l-m > 20) || (l-s > 20))
        term1 = (2*l+1)*(sterlings(l+m)/sterlings(l+s))*(sterlings(l-m)/sterlings(l-s))
    else
        term1 = (2*l+1)*(factorial(l+m)/factorial(l+s))*(factorial(l-m)/factorial(l-s))
    end
    #println("HERE WE HAVE",s," ", l," ", m)
    A = sqrt(term1/(4*pi))*(-1)^m

    #The sin term right outside the summation
    sinterm = ((1-z)/2)^l
    #The summation terms
    sumterms = 0.0
    for r = 0:(l-s)
        consts = binomial(l-s,r)*binomial(l+s,r+s-m)*(-1)^(l-r-s)
        cotterm = (sqrt((1+z)/(1-z)))^((2*r+s-m))
        sumterms += consts*cotterm
    end
    return A*sinterm*sumterms
end
""" |> Meta.parse |> eval


text:
using LinearAlgebra

@inline function exp_gen!(cache,A,::Val{3})
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=6
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = getmem(cache,2)
    memslots3 = getmem(cache,3)
    memslots4 = getmem(cache,4)
    memslots5 = getmem(cache,5)
    memslots6 = getmem(cache,6)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    # Computation order: A2 A4 A6 Ua U V Z X P
    # Computing A2 with operation: mult
    mul!(memslots2,memslots1,memslots1)
    # Computing A4 with operation: mult
    mul!(memslots3,memslots2,memslots2)
    # Computing A6 with operation: mult
    mul!(memslots4,memslots2,memslots3)
    # Computing Ua = x*I+x*A2+x*A4+x*A6
    coeff1=8.64864e6
    coeff2=277200.0
    coeff3=1512.0
    coeff4=1.0
    memslots5 .= coeff2.*memslots2 .+ coeff3.*memslots3 .+ coeff4.*memslots4
    inplace_add!(memslots5,I*coeff1)
    # Computing U with operation: mult
    mul!(memslots6,memslots5,memslots1)
    # Deallocating Ua in slot 5
    # Deallocating A in slot 1
    # Computing V = x*I+x*A2+x*A4+x*A6
    coeff1=1.729728e7
    coeff2=1.99584e6
    coeff3=25200.0
    coeff4=56.0
    # Smart lincomb recycle A2
    memslots2 .= coeff2.*memslots2 .+ coeff3.*memslots3 .+ coeff4.*memslots4
    inplace_add!(memslots2,I*coeff1)
    # Deallocating A4 in slot 3
    # Deallocating A6 in slot 4
    # Computing Z = x*V+x*U
    coeff1=1.0
    coeff2=-1.0
    memslots1 .= coeff1.*memslots2 .+ coeff2.*memslots6
    # Computing X = x*V+x*U
    coeff1=1.0
    coeff2=1.0
    # Smart lincomb recycle V
    memslots2 .= coeff1.*memslots2 .+ coeff2.*memslots6
    # Deallocating U in slot 6
    # Computing P with operation: ldiv
    ldiv_for_generated!(memslots3, memslots1, memslots2)
    # Deallocating Z in slot 1
    # Deallocating X in slot 2
    copyto!(A,memslots3) # Returning P
end

max_stars_repo_name: JuliaDiffEq/ExponentialUtilities

max_issues_repo_name: JuliaDiffEq/ExponentialUtilities

max_forks_repo_name: JuliaDiffEq/ExponentialUtilities


text:
using LOWDER

fmin = Array{Function}(undef, 11)

for i = 1:11

    y = i + 14.0

    fmin[i] = x -> 2.0 * x[1] ^ 2.0 + 0.008 * y ^ 3.0 - 3.2 * y * x[1] - 2.0 * y

end

n_points = 2
l = [3.0]
u = [5.5]
xinit = [0.0]

sol = LOWDER.lowder( fmin, xinit, l, u; m = n_points )


text:
using Wordle
using Test
using TOML
using SHA
using Suppressor

data = TOML.parsefile(joinpath(@__DIR__(), "Data.toml"))

@testset "Wordle" begin
    include("fetch.jl")
    include("game.jl")
    include("guess.jl")
end


text:
"""
Contains the data for a BlockMap
"""
mutable struct BlockMapData{GID <: Integer, PID <:Integer, LID <: Integer}
    comm::Comm{GID, PID, LID}
    directory::Union{Directory, Nothing}
    lid::Vector{LID}
    myGlobalElements::Vector{GID}
#    firstPointInElementList::Array{Integer}
#    elementSizeList::Array{Integer}
#    pointToElementList::Array{Integer}

    numGlobalElements::GID
    numMyElements::LID
#    elementSize::Integer
#    minMyElementSize::Integer
#    maxMyElementSize::Integer
#    minElementSize::Integer
#    maxElementSize::Integer
    minAllGID::GID
    maxAllGID::GID
    minMyGID::GID
    maxMyGID::GID
    minLID::LID
    maxLID::LID
#    numGlobalPoints::Integer
#    numMyPoints::Integer

#    constantElementSize::Bool
    linearMap::Bool
    distributedGlobal::Bool
    oneToOneIsDetermined::Bool
    oneToOne::Bool
    lastContiguousGID::GID
    lastContiguousGIDLoc::GID
    lidHash::Dict{GID, LID}
end

function BlockMapData(numGlobalElements::GID, comm::Comm{GID, PID, LID}) where GID <: Integer where PID <: Integer where LID <: Integer
    BlockMapData(
        comm,
        nothing,
        LID[],
        GID[],

        numGlobalElements,
        LID(0),
        GID(0),
        GID(0),
        GID(0),
        GID(0),
        LID(0),
        LID(0),

        false,
        false,
        false,
        false,
        GID(0),
        GID(0),
        Dict{GID, LID}()
    )
end


text:
module HiveLoader

# https://github.com/JuliaDatabases/Hive.jl v0.3.0
using Hive # HiveSession HiveAuth

using Octo: Repo, AdapterBase, DBMS, SQLElement, Structured
using .Repo: SQLKeyword, ExecuteResult

# db_dbname
function db_dbname(nt::NamedTuple)::String
    ""
end

# db_connect
function db_connect(; host::String="localhost", port::Integer=10000, auth::HiveAuth=HiveAuth(), tprotocol::Symbol=:binary)
    HiveSession(host, port, auth; tprotocol=tprotocol)
end

# db_disconnect
function db_disconnect(sess)
    if sess isa HiveSession
        Hive.close(sess)
    end
end

# query
function query(sess, sql::String)
    pending = Hive.execute(sess, sql)
    rs = Hive.result(pending)
    sch = Hive.schema(rs)
    column_names = tuple(Symbol.(getproperty.(sch.columns, :columnName))...)
    df = reduce(vcat, Hive.records(rs))
    nts = NamedTuple{column_names}.(df)
    Hive.close(rs)
    nts
end

function query(sess, prepared::String, vals::Vector) # throw UnsupportedError
    throw(UnsupportedError("needs to be implemented"))
end

# execute
function execute(sess, sql::String)::ExecuteResult
    result = Hive.execute(sess, sql)
    nothing
end

function execute(sess, prepared::String, vals::Vector)::ExecuteResult # throw UnsupportedError
    throw(UnsupportedError("needs to be implemented"))
end

function execute(sess, prepared::String, nts::Vector{<:NamedTuple})::ExecuteResult # throw UnsupportedError
    throw(UnsupportedError("needs to be implemented"))
end

# execute_result
function execute_result(sess, command::SQLKeyword)::NamedTuple
    NamedTuple()
end

function Base.show(io::IO, mime::MIME"text/plain", element::Union{E,Structured} where E<:SQLElement)
    dbms = DBMS.Hive()
    AdapterBase._show(io, mime, dbms, element)
end

end # module Octo.Backends.HiveLoader


text:
"""
    AdaptiveExtrapolationD

This module is  implementing an
adaptive extrapolation of the explicit midpoint rule according to Deuflhard.
"""
module AdaptiveExtrapolationD

"""
    (Δ,Δx,statisitic) = (mySolver::solver)(f::Function, x₀::Vector{T}, t₀::S, tEnd::S; <options>)
        where {T<:Number,S<:AbstractFloat}

Computes the grid function `xΔ` on the grid `Δ` approximating the solution of the initial value problem

        x′ = f(t,x),  x(t₀) = x₀

on the interval `[t₀,tEnd]` (or `[tEnd,t₀]` if `tEnd < t₀`). The integrator is the extrapolated explicit midpoint rule.
The algorithm adaptively controls stepsize and extrapolation order in the line of thought of Deuflhard.
Note that `solver` is a callable instance of the `struct solver`.

***
## Initialize `solver`

    solver(N::Integer, sequence::String)

Initializes and returns the an instance of `solver`.
`N ≧ 1` is the maximal order of extrapolation. `sequence` specifies the subdividing sequence
used. The options are
* `Harmonic`, that is `1, 2, 3, 4, 5, 6,...`
* `Romberg`, that is `1, 2, 4, 8, 16, 32,...`
* `Bulirsch`, that is `1, 2, 3, 4, 6, 8,... `
The structure contains all quantities for the adaptive extrapolation of the the explicit
midpoint rule for all orders in `1:N` that can be tabulated before the actual computation .
These are the *subdividing sequence* and the *weights* for the extrapolation based on the first barycentric formula.

***
## Options

The following options can be passed as keyword arguments.
* `tol::S > 0` is tolerance used for accuracy check. Default is `1e-3`.
* `relativeScaling::Vector{S}` and `absoluteScaling::Vector{S}` are vectors of the same length as `x₀`
    containing no negative elements.
    Default (in every component) is `relativeScaling[i] = 1.0` and `absoluteScaling[i] = 1e-3`.
    The accuracy check is passed iff the estimated error satisfies in every component:

        |relativeError[i]| ≦ relativeScaling[i]⋅tol & |absoluteError[i]| ≦ absoluteScaling[i]⋅tol

* `firstGuessStepsize::S` is the length of the first step. Its modulus must be greater than
    `eps(S)` Default is `1e-3`.
* `maximalSteplength::S > 0` is the maximal absolute stepsize used for the integration. Default is `Inf`.
* `minimalStepsizeScaling::S` and `maximalStepsizeScaling::S` are safety factors for the stepsize
    selection. Default are `0.25` and `4.0`. The new stepsize satisfies:

         0 ≦ minimalStepsizeScaling ≦ newStepsize / oldStepsize ≦ maximalStepsizeScaling

* `safetyStepsizeSelection::S>0`. Use `safetyStepsizeSelection⋅tol` instead of `tol`. Default is `0.25`.
* `minimumOrder::Int64 ≧ 1` and `maximumOrder::Int64 ≦ N` are the minimal and maximal order of extrapolation for the integration (`N` is the number that has been passed to the constructor). Default are `1` and `N`.
* `firstGuessOrder::Int64` is the extrapolation order for the first step. It must satisfy
    `minimumOrder ≦ firstGuessOrder ≦ maximumOrder`. Default is `round(Int64,(maximumOrder-minimumOrder)/2)`.
* `maximumReduction::Int64` is the upper bound of the number of successive reductions allowed per step.
    Exceeding this boundary leads to ending the integration prematurely (`Δ[end] ≠ tEnd`). Default is `10`.
* `maximumSteps::Int64 ≧ 0` is the maximal length of the array `Δ`. Exceeding this boundary leads to
    ending the integration prematurely (`Δ[end] ≠ tEnd`). Default is `10 000`.
* `rescaleWeights::Bool` If `true` Lagrange interpolation is used otherwise the frist barycentric fromula. Default is `false`
* `detailedStatistics::Bool` flag for the content of `statisitic`. Default is `false`
    + `true`: `statisitic[:,i]` contains [extrapolation order; stepsize; number of `f`-evaluations;
        number of reductions] for the ith step.
    + `false`: `statisitic` contains [total number of evaluations of the right side `f`, total number of reductions].
"""
struct solver
    """
    The structure contains all quantities for the adaptive extrapolation of the the explicit
    midpoint rule that can be tabulated before the acutal computation.
    These are for all orders in `1:N` (`N` is provided the constructor):
    The subdividing sequence, the weights and scaling used for
    the extrapolation based on the first barycentric formula.
    """
   subdividingSequence::Array{BigInt,1}

   # weights and scaling factors for extrapolation operators
   ω::Array{Rational{BigInt},2}
   ρ::Array{Rational{BigInt},1}

   # weights and scaling factors for internal extrapolation operators (used for error estimate)
   ω2::Array{Rational{BigInt},2}
   ρ2::Array{Rational{BigInt},1}

   # constructor
   function solver(N::Integer, sequence::String)
       # check input and initialize subdividing sequence subdividingSequence for barycentric weights
       if N < 0
           error("Order of extrapolation must not be negative. But I got N = $N")
       end
       if sequence == "Harmonic"
           subdividingSequence = [BigInt(n+1) for n = 0:N]
       elseif sequence == "Romberg"
           subdividingSequence = [BigInt(2)^n for n = 0:N]
       elseif sequence == "Bulirsch"
           subdividingSequence = [n==0 ? BigInt(1) : (isodd(n) ? BigInt(2)^Int64(n/2+0.5) : 3*BigInt(2^Int64(n/2-1))) for n = 0:N]
       else
           error("Name of subdividing sequence must be ''Harmonic'', ''Romberg'' or ''Bulirsch''. But I got ''$sequence''")
       end

       # compute nodes corresponding to the subdividing sequence subdividingSequence
       nodes = BigInt(1).// subdividingSequence.^2

       # compute barycentric weights for internal extrapolation operators
       ω2 = zeros(Rational{BigInt},N,N)
       ω2[1,:] = ones(Rational{BigInt},1,N)
       for n = 2:N
           distance = nodes[2:n] .- nodes[n+1]
           ω2[1:(n-1),n] = ω2[1:n-1,n-1] .// distance
           ω2[n,n] = 1 // prod(-distance)
       end

       # compute barycentric weights for extrapolation operators
       ω = zeros(Rational{BigInt},N+1,N+1)
       for n = 1:N
           ω[n+1,(n+1):(N+1)] = ω2[n,n:N] // (nodes[n+1]-nodes[1])
           ω[1,n] = 1 // prod(nodes[1].-nodes[2:n])
       end
       ω[1,N+1] = 1 // prod(nodes[1].-nodes[2:N+1])

       #rescale barycentric weights to obtain weights of 1. barycentric formula
       for m = 1:(N+1)
           ω[1:m,m] = - ω[1:m,m] .// nodes[1:m]
           if 2 <= m
               ω2[1:m-1,m-1] = - ω2[1:m-1,m-1] .// nodes[2:m]
           end
       end

       # compute scaling factors for internal extrapolation operators
       ρ2 = ones(Rational{BigInt},N)
       ρ2[1] = -nodes[2]
       for n = 1:(N-1)
           ρ2[n+1] = -ρ2[n]*nodes[n+2]
       end

       # compute scaling factors for extrapolation operators
       ρ = -nodes[1]*[BigInt(1); ρ2]

       # initialize structure
       new(subdividingSequence,ω,ρ,ω2,ρ2)
   end
end

function (mySolver::solver)(
            f::Function, x₀::Vector{T}, t₀::S, tEnd::S;
            tol::S = parse(S,"1e-3"),
            relativeScaling::Vector{S} =ones(S,size(x₀)),
            absoluteScaling::Vector{S} = fill(parse(S,"1e-3"),size(x₀)),
            firstGuessStepsize::S = parse(S,"1e-3"),
            maximalSteplength::S = parse(S,"Inf"),
            minimalStepsizeScaling::S = parse(S,"0.02"),
            maximalStepsizeScaling::S = parse(S,"4"),
            safetyStepsizeSelection::S =  parse(S,"0.25"),
            minimumOrder::Int64 = 1,
            maximumOrder::Int64 = length(mySolver.subdividingSequence)-1,
            firstGuessOrder = round(Int64,(maximumOrder-minimumOrder)/2),
            maximumReduction::Int64 = 10,
            maximumSteps::Int64 = 10000,
            rescaleWeights::Bool = false,
            detailedStatistics::Bool = false,
            kwargs...
            )where {T<:Number,S<:AbstractFloat}
    # initialize
    # 1. constants
    d = length(x₀) # systems dimension
    if 1<= minimumOrder <= maximumOrder <= length(mySolver.ρ)-1
        n_ex, N_ex = minimumOrder, maximumOrder
    else
        error("Minimal and maximal order must satisfy: 1 ≦ minimumOrder ≦ maximumOrder ≦ $(length(mySolver.ρ)-1)!")
    end

    if length(relativeScaling) == length(absoluteScaling) == d
        if all(el->el>=zero(S),relativeScaling ) && all(el->el>=zero(S),absoluteScaling)
            σᵣ, σₐ = relativeScaling,  absoluteScaling # scaling factors for the  error estimate
        else
            error("All elements of relativeScaling and absoluteScaling must be non-negative!")
        end
    else
        error("relativeScaling and absoluteScaling must be vectors of length $(d)!")
    end

    if zero(S) <= minimalStepsizeScaling <= maximalStepsizeScaling
        p1, p2 = minimalStepsizeScaling, maximalStepsizeScaling
    else
        error("Minimal and maximal stepsize scaling must satisfy: 0.0 <= minimalStepsizeScaling <= maximalStepsizeScaling!")
    end

    # 2. arrays
    if maximumSteps >= 0
        Δ = zeros(S, maximumSteps + 1) # time grid
        xΔ = zeros(T, d, maximumSteps + 1) # solution
    else
        error("maximumSteps must not be negative!")
    end
    X = zeros(T, d, N_ex + 1) # storage for the internal discretisations obtained by the explicit midpoint rule
    λ = zeros(S, N_ex - n_ex +1) # storage for scaling factors of stepsize. λ[k] contains the scalar for extr. order (k + n_ex - 1)
    τ_opt =zeros(S, N_ex - n_ex +1) # storage for optimal stepsize. τ_opt[k] contains the stepsize for extr. order (k + n_ex - 1)
    s = [2*sum(Int64.(mySolver.subdividingSequence[1:n+1])) - n for n in n_ex:N_ex] # s[k] is the number of stages for  extrapolation order (k + n_ex - 1)
    statistics = detailedStatistics ? zeros(S,4,maximumSteps) : zeros(Int64,2)

    # 3. initialization for integration loop
    Δ[1], xΔ[:,1] = t₀, x₀ # store inital values
    tRest = tEnd - Δ[1] # remaining size of the interval
    if maximalSteplength > zero(S)
        maximalSteplength = min(abs(tRest),maximalSteplength)
    else
        error("maximalSteplength must be positive!")
    end

    if abs(firstGuessStepsize) > eps(S)
        τ = flipsign(min(abs(firstGuessStepsize),maximalSteplength),tRest) # stepsize for first step, accomodate backward integration
    else
        error("Absolute value of firstGuessStepsize must be greater than eps($S) = $(eps(S))!")
    end
    if n_ex <= firstGuessOrder <= N_ex
        m = firstGuessOrder # extrapolation order for first step
    else
        error("firstGuessOrder must be between $n_ex and $(N_ex)!")
    end

    counterStep = 1
    counterReduction = 0
    counterEvaluation = 0
    converged = false

    finished = maximumSteps == 0

    x_n = ones(T,d) # storage for the latest solution
    ϵ_n = zero(S) # storage for the latest error estimate
    if safetyStepsizeSelection > 0
        tol = safetyStepsizeSelection*tol # the tolerance for accuracy checks
    else
        error("safetyStepsizeSelection must be positive!")
    end

    # integration loop
    while !finished
        # the nth step: the approximation for the time t₀ + τ is saved in xΔ[n+1].

        # check if the computation must be aborted
        if counterStep > maximumSteps
            println("Not finished. Exceeded maximal number of steps ($maximumSteps)")
            return (Δ[1:counterStep], xΔ[:,1:counterStep], detailedStatistics ? statistics[:,1:counterStep-1] : statistics)
        end
        if counterReduction > maximumReduction
            println("Not finished. Too many reductions in step number $(counterStep)!")
            return (Δ[1:counterStep], xΔ[:,1:counterStep], detailedStatistics ? statistics[:,1:counterStep-1] : statistics)
        end
        if abs(τ) <= eps(S)
            println("Stepsize is becoming too small in step $(counterStep)!")
            return (Δ[1:counterStep], xΔ[:,1:counterStep], detailedStatistics ? statistics[:,1:counterStep-1] : statistics)
        end

        # preperations for current step:
        n_win, N_win = max(n_ex, m - 1), min(N_ex, m + 1) # order window
        n = n_win # start with smalles order in the order window

        # compute all necessary internal approximations:
        f₀ = f(t₀,x₀) # the information for the Euler step
        counterEvaluation = counterEvaluation + 1

        for i = 0:n
            X[:,i+1] = explicitMidpointRule(t₀,x₀,f₀,f,τ,2*Int64(mySolver.subdividingSequence[i+1]))
            counterEvaluation = counterEvaluation + 2*Int64(mySolver.subdividingSequence[i+1]) - 1
        end

        #  compute all information relating to an extrapolation order ≦ n_win:
        for i = n_ex : n
            if rescaleWeights
                # rescaling weights is equivalent to Lagrange interpolation
                x_i= X[:, 1:(i+1)]*T.(broadcast(*,mySolver.ω,mySolver.ρ')[1:(i+1),(i+1)]) # discretisation of order i
                xx_i = X[:, 2:(i+1)]*T.(broadcast(*,mySolver.ω2,mySolver.ρ2')[1:i, i]) # its internal counterpart
            else
                x_i = T.(mySolver.ρ[i+1])*(X[:, 1:(i+1)]*T.(mySolver.ω[1:(i+1), (i+1)]))
                xx_i = T.(mySolver.ρ2[i])*(X[:, 2:(i+1)]*T.(mySolver.ω2[1:i, i]))
            end
            σ = max.(σₐ, σᵣ.*abs.(x_i)) # scaling for error estimate
            ϵ_i = sqrt(sum((abs.(x_i-xx_i)./σ).^2)/d) # scaled mean square norm of the estimated error
            λ[i-n_ex+1] = min(p2, max(p1, (tol/ϵ_i)^(1/(2i+1)))) # factor for optimal stepsize, including safety scaling
            if i == n
                x_n, ϵ_n = x_i, ϵ_i # save solution and error for the current order n
            end
        end

        # check if soltution with for n in the order window can be accepted
        while n <= N_win
            if ϵ_n <= tol
                # accept order n
                converged = true
                break
            elseif ϵ_n <= tol^(s[n-n_ex+1]/s[N_win-n_ex+1])
                # reject order n but pass convergence monitor
                n = n + 1
                # compute x_n, λ_n and ϵ_n for new n:
                X[:,n+1] = explicitMidpointRule(t₀,x₀,f₀,f,τ,2*Int64.(mySolver.subdividingSequence[n+1]))
                counterEvaluation = counterEvaluation + 2Int64(mySolver.subdividingSequence[n+1]) - 1
                if rescaleWeights
                    # rescaling weights is equivalent to Lagrange interpolation
                    x_n= X[:, 1:(n+1)]*T.(broadcast(*,mySolver.ω,mySolver.ρ')[1:(n+1),(n+1)])
                    xx_n = X[:, 2:(n+1)]*T.(broadcast(*,mySolver.ω2,mySolver.ρ2')[1:n, n])
                else
                    x_n = T.(mySolver.ρ[n+1])*(X[:, 1:(n+1)]*T.(mySolver.ω[1:(n+1), (n+1)])) # approximation of extrapolation order n
                    xx_n = T.(mySolver.ρ2[n])*(X[:, 2:(n+1)]*T.(mySolver.ω2[1:n, n])) # and its internal counterpart
                end
                σ = max.(σₐ, σᵣ.*abs.(x_n)) # scaling for error estimate
                ϵ_n = sqrt(sum((abs.(x_n-xx_n)./σ).^2)/d) # scaled mean square norm of the estimated error
                λ[n-n_ex+1] = min(p2,max(p1,(tol/ϵ_n)^(1/(2n+1)))) # factor for optimal stepsize, including safety scaling

            else
                # reject order n and not pass convergence monitor
                break
            end
        end

        if converged
            #  step is accepted, save x_n for the time t₀ + τ and update statistics
            Δ[counterStep+1] = t₀ + τ
            xΔ[:, counterStep+1] = x_n
            tRest = tEnd - Δ[counterStep+1]
            maximalSteplength = min(abs(tRest),maximalSteplength)
            if detailedStatistics
                statistics[:,counterStep] = [n, τ, counterEvaluation, counterReduction]
            else
                statistics = statistics + [counterEvaluation,counterReduction]
            end

            #  check if tEnd is already reached
            if abs(tRest) < eps(S)
                finished = true
            else
                # preperations for the next step:
                t₀, x₀  = Δ[counterStep+1], xΔ[:, counterStep+1] # new inital data

                ## compute optimal order m and stepsize τ for the next step acc. Deuflhard
                # first compute the optimal extrapolation order
                temp = (n_ex:n) .- n_ex .+ 1 # index range of computed quantities
                τ_opt[temp] = min.(abs(τ)*λ[temp], maximalSteplength) # safety boundary for τ_opt
                ω = s[temp]./(τ_opt[temp]) # work per step

                m = argmin(ω) + n_ex - 1 # optimal order

                # check if we may increase m
                if m == n < N_win
                    # compute scaling for order m+1, including safety scaling
                    λ[n-n_ex+2] = min(p2, max(p1, (tol^(s[n-n_ex+1]/s[n-n_ex+2])/ϵ_n)^(1/(2n+1))))
                    # check if work dereases from order m to  m+1
                    if ω[end] > s[n-n_ex+2]/(τ_opt[n-n_ex+2] = min(abs(τ)*λ[n-n_ex+2], maximalSteplength))
                        m = m + 1
                    end
                end
                τ = flipsign(τ_opt[m-n_ex+1],τ) # accomodate backward integration

                counterStep = counterStep + 1
                converged = false
                counterEvaluation = 0
                counterReduction = 0
                fill!(λ,zero(S))
                fill!(τ_opt,zero(S))
            end
        else
            counterReduction = counterReduction + 1
            # compute reduced stepsize τ for order m. Use latest error estimate ϵ_n if estimate of order m is not available
            if n < m
                 λ[m-n_ex+1] = min(p2,max(p1,(tol^(s[n-n_ex+1]/s[m-n_ex+1])/ϵ_n)^(1/(2n+1))))
            end
            τ_opt[m-n_ex+1] = min(abs(τ)*λ[m-n_ex+1],maximalSteplength) # safety boundary for τ_opt
            τ = flipsign(τ_opt[m-n_ex+1],τ) # accomodate backward integration
        end
    end # integration loop
    return (Δ[1:counterStep+1], xΔ[:,1:counterStep+1], detailedStatistics ? statistics[:,1:counterStep] : statistics)
end

"""
    explicitMidpointRule(t₀,x₀,f₀,f,τ,N)

Computes the value ``x_\\Delta(t_0 + \\tau)`` recursively. Here ``x_\\Delta`` is the grid-function generated
by the explicit midpoint rule on the equidistant grid `[t₀ + n/N*τ for n = 0:N]`
with the initial data `(t₀,x₀)` and the right side `f`.
In the starting step `f₀` is used instead of `f(t₀,x₀)`.
"""
function explicitMidpointRule(t₀::S,x₀::Vector{T},f₀::Vector{T},f::Function,τ::S,N) where{S<:AbstractFloat, T<:Number}
    τ_n = τ/N
    x = x₀ + τ_n*f₀
    for i = 0:N-2
        temp = x
        x = x₀+2*τ_n*f(t₀+(i+1)*τ_n,x)
        x₀ = temp
    end
    return x
end

export solver
end # module AdaptiveExtrapolationD

max_stars_repo_name: AlthausKonstantin/Extrapolation

max_issues_repo_name: AlthausKonstantin/Extrapolation

max_forks_repo_name: AlthausKonstantin/Extrapolation


text:
using Rankers.Storages # update_with
using LOLTools.SpectatorV4: SpectatorV4, http_action

api_key = get(ENV, "RIOT_TOKEN", "")
region = "kr"
read_from_file = true # false
featured = SpectatorV4.featured_games(api_key, region; action=update_with(region, http_action, read_from_file=read_from_file, debug=true))

n = 0
for gameinfo in featured.gameList
    for participant in gameinfo.participants
        print(participant.summonerName)
        print(" ---- ")
        global n += 1
    end
end
println()
println("count: $n")


text:
using FSimBase

using LinearAlgebra  # for I, e.g., Matrix(I, n, n)
using ComponentArrays
using UnPack
using Transducers
using Plots
using DifferentialEquations

struct MyEnv <: AbstractEnv  # AbstractEnv exported from FSimBase
    a
    b
end

"""
FlightSims recommends you to use closures for State and Dynamics!. For more details, see https://docs.julialang.org/en/v1/devdocs/functions/.
"""
function State(env::MyEnv)
    return function (x1::Number, x2::Number)
        ComponentArray(x1=x1, x2=x2)
    end
end

function Dynamics!(env::MyEnv)
    @unpack a, b = env  # @unpack is very useful!
    @Loggable function dynamics!(dx, x, p, t; u)  # `Loggable` makes it loggable via SimulationLogger.jl (imported in FSimBase)
        @unpack x1, x2 = x
        @log x1  # to log x1
        @log x2  # to log x2
        dx.x1 = a*x2
        dx.x2 = b*u
    end
end

function my_controller(x, p, t)
    @unpack x1, x2 = x
    -(x1+x2)
end

function main()
    n = 2
    m = 1
    a, b = 1, 1
    A = -Matrix(I, n, n)
    B = Matrix(I, m, m)
    env = MyEnv(a, b)
    tf = 10.0
    Δt = 0.01
    x10, x20 = 10.0, 0.0
    x0 = State(env)(x10, x20)
    # prob: DE problem, df: DataFrame
    @time prob, df = sim(
                         x0,  # initial condition
                         apply_inputs(Dynamics!(env); u=my_controller);  # dynamics!; apply_inputs is exported from FSimBase and is so useful for systems with inputs
                         solver=Tsit5(),
                         tf=10.0,
                         savestep=Δt,  # savestep is NOT simulation step
                        )  # sim is exported from FSimBase
    ts = df.time
    x1s = df.sol |> Map(datum -> datum.x1) |> collect
    x2s = df.sol |> Map(datum -> datum.x2) |> collect
    # plot
    p_x1 = plot(ts, x1s;
                label="x1",
               )
    p_x2 = plot(ts, x2s;
                label="x2",
               )
    p_x = plot(p_x1, p_x2, layout=(2, 1))
    # save
    dir_log = "figures"
    mkpath(dir_log)
    savefig(p_x, joinpath(dir_log, "toy_example.png"))
    display(p_x)
end


text:
@assert 1 == 1

#=
CREATE TABLE test.types (
    bigints bigint PRIMARY KEY,
    blobs blob,
    booleans boolean,
    dates date,
    decimals decimal,
    doubles double,
    floats float,
    inets inet,
    ints int,
    smallints smallint,
    texts text,
    timestamps timestamp,
    timeuuids timeuuid,
    times time,
    tinyints tinyint,
    uuids uuid,
    varchars text,
    varints varint
) WITH bloom_filter_fp_chance = 0.01
    AND caching = {'keys': 'ALL', 'rows_per_partition': 'ALL'}
    AND comment = ''
    AND compaction = {'class': 'SizeTieredCompactionStrategy'}
    AND compression = {'sstable_compression': 'org.apache.cassandra.io.compress.LZ4Compressor'}
    AND crc_check_chance = 1.0
    AND dclocal_read_repair_chance = 0.1
    AND default_time_to_live = 0
    AND gc_grace_seconds = 864000
    AND max_index_interval = 2048
    AND memtable_flush_period_in_ms = 0
    AND min_index_interval = 128
    AND read_repair_chance = 0.0
    AND speculative_retry = '99.0PERCENTILE';
=#


text:
mutable struct Controller3 <: Controller
    Δt::Real
    g::Real
    M::Real

    c::Real
    k::Real
    Fs::Real

    K::AbstractMatrix
    F::Function
    control!::Function

    function Controller3(F,Δt,g,M,m,l,c,k,Fs,Q,R)
        I=1/3*m*l^2
        v1=(M+m)/(I*(M+m)+(l^2*m*M))
        v2=(I+l^2*m)/(I*(M+m)+(l^2*m*M))

        A = [0 1 0 0
            0 -k*v2 ((l*m)^2*g*v2)/(I+m*l^2) (-l*m*c*v2)/(I+m*l^2)
            0 0 0 1
            0 -(l*m*k*v1)/(M+m) l*m*g*v1 -1*c*v1
        ]

        B = [
            0
            v2
            0
            l*m*v1/(M+m)
        ]

        P = ConstrainedControl.care(A,B,Q,R)
        K = R\B'*P

        new(Δt,g,M,c,k,Fs,K,F,control!)
    end
end

function control!(mechanism,controller::Controller3,kstep)
    Δt = controller.Δt
    g = controller.g
    M = controller.M
    c = controller.c
    k = controller.k
    Fs = controller.Fs
    K = controller.K
    F = controller.F
    cart = mechanism.bodies[1]
    pend = mechanism.bodies[2]

    x = cart.state.xc[2]
    v = cart.state.vc[2]
    q = pend.state.qc
    θ = (rotation_angle(q) * rotation_axis(q))[1]
    ω = pend.state.ωc[1]

    z = [x;v;θ;ω]

    Fcart = F(z,K,kstep*Δt)
    Fcart, Fpend = friction(Fcart, 0, v, ω, k, c, Fs)

    setForce!(mechanism, geteqconstraint(mechanism, 3), [Fcart])
    setForce!(mechanism, geteqconstraint(mechanism, 4), [Fpend])

    return
end

function run3!(str; vis::Bool=true)
    open("Files/E3.jl",create=true,write=true) do file
        write(file,"
        Δt = 0.004

        # Parameters from Table 1.1
        g = 9.81
        L1 = 0.28125
        l1 = L1/2
        l = l1 # (l1-0.135)/2 + l1
        M = 4.5
        m = 0.2585
        # I = 1/3*m*l^2
        c = 0.002
        k = 0.3
        Fs = 8

        # Joint axes
        ex = [1.0;0.0;0.0]
        ey = [0.0;1.0;0.0]

        p = [0.0;0.0;l1] # joint connection point

        # Desired orientation
        θ1 = 0

        # Links
        origin = Origin{Float64}()
        cart = Box(0.1, 0.5, 0.1, M)
        pend = Box(0.055, 0.055, L1, m)

        # Constraints
        joint1 = EqualityConstraint(Prismatic(origin, cart, ey))
        joint2 = EqualityConstraint(Revolute(cart, pend, ex; p2 = -p))

        links = [cart;pend]
        constraints = [joint1;joint2]

        mech = Mechanism(origin, links, constraints, Δt = Δt, g = -g)

        xinit = 0.0
        theta1init = pi
        Q = [
            1 0 0 0
            0 1 0 0
            0 0 1 0
            0 0 0 1
        ]
        R = 1
        F(z,K) = 0

        ### Begin Student Input

        ")
        write(file, str)
        write(file, "

        ### End Student Input

        setPosition!(origin,cart,Δx = [0;xinit;0])
        setPosition!(cart,pend,p2 = -p, Δq = UnitQuaternion(RotX(theta1init)))

        controller = Controller3(F,Δt,g,M,m,l,c,k,Fs,Q,R)

        steps = Base.OneTo(Int(10/Δt))
        storage = Storage{Float64}(steps,2)

        try
            simulate!(mech,storage,controller,record = true)
        catch
            println(\"Unstable behavior\")
        end")

        if vis
            write(file, "
            visualize(mech,storage)")
        end
    end

    try
        include("Files/E3.jl")
    catch
        println("Error. Bad code.")
    end

    return
end


text:
function StatsBase.denserank(W::T; kw...) where {T <: Womble}
    return StatsBase.denserank(W.m; kw...)
end


text:
using GaussianDistributions
using Random
using Test
using Statistics
using LinearAlgebra

n = 500
Random.seed!(2)
L1 = rand(3,3)
Σ = L1*L1'
L2 = rand(2,2)
R = L2*L2'
H = rand(2,3)
μ = rand(3)
X = L1*randn(3, n)
Y = H*X + L2*randn(2, n)

@test norm(cov(vcat(X,Y)') - [L1*L1' L1*L1'*H' ; H*L1*L1' (H*L1*L1'*H' + L2*L2')]) < 20/sqrt(n)

function correct2(x, (y, R), H)
    GaussianDistributions.conditional(
    Gaussian([x.μ; H*x.μ], [x.Σ x.Σ*H'; H*x.Σ (H*x.Σ*H' + R)]),
    [1,2,3],
    [4,5],
    y)
end
x = Gaussian(μ, Σ)
y = H*rand(x) + rand(Gaussian(zeros(2), R))
z1 = GaussianDistributions.correct(x, (y,R), H)[1]
z2 = correct2(x, Gaussian(y,R), H)
@test z1 ≈ z2


text:
const _GW_STAGING_TOUT = 60.0

function _upload_jltask(gw::GitWorker,
        tid::String, ex::Expr;
        desc = "", src = _expr_src(ex),
        readme = false,
        vtime = _GW_DFLT_TASK_VTIME
    )

    tout = _GW_STAGING_TOUT
    gl = gitlink(gw)

    return GitLinks.writewdir(gl; tout) do _

    end # writewdir
end

max_stars_repo_name: josePereiro/GitWorkers

max_issues_repo_name: josePereiro/GitWorkers

max_forks_repo_name: josePereiro/GitWorkers


text:
##
## Particular divergent series with peculiar partial sums behaviour
##
## Tomáš Kalvoda, FIT CTU, 2015
##

using Images, Colors, FixedPointNumbers

# summands
u(n) = n*log(n)*sqrt(2)
z(n) = exp(2im*pi*u(n))

# first maxn partial sums
function partial_sums(maxn)
    out = Array(Complex{Float64}, maxn)
    out[1] = z(1)
    for j = 2:maxn
        out[j] = out[j-1] + z(j)
    end
    return out
end

# draw a point at coordinate p in image
function draw_point(p, image)
    x,y = p
    image[y, x, :] = [0.,0.,0.]
    image[y+1, x, :] = [0.,0.,0.]
    image[y, x+1, :] = [0.,0.,0.]
    image[y-1, x, :] = [0.,0.,0.]
    image[y, x-1, :] = [0.,0.,0.]
end

# draw a line connecting p1 and p2
function draw_line(p1, p2, image)
    x1,y1 = p1; x2,y2 = p2;
    d = convert(Int64, ceil(sqrt((x1-x2)^2 + (y1-y2)^2)))
    for j = 0:d
        x = convert(Int64, floor(x1 + j/d*(x2-x1)))
        y = convert(Int64, floor(y1 + j/d*(y2-y1)))
        draw_point((x,y), image)
    end
end

# create the image
function draw_image(sums, width, height)
    real_parts = real(sums)
    imag_parts = imag(sums)
    remin, remax = (1.05*minimum(real_parts), 1.05*maximum(real_parts))
    immin, immax = (1.05*minimum(imag_parts), 1.05*maximum(imag_parts))
    redelta, imdelta = ((remax - remin) / width, (immax - immin) / height)
    tox(re) = convert(Int64, ceil((re-remin)/redelta))
    toy(im) = convert(Int64, ceil((im-immin)/imdelta))

    img = Array(Float64, (height,width,3))
    fill!(img, 1.0)

    for j in 2:length(sums)
        p1 = (tox(real(sums[j-1])), toy(imag(sums[j-1])))
        p2 = (tox(real(sums[j])), toy(imag(sums[j])))
        draw_line(p1, p2, img)
    end

    return flipud(img)
end

# run the computation
function draw_sums(maxn; width=400, height=400)
    s = partial_sums(maxn)
    h = maximum(abs(s))

    img = draw_image(s, width, height)

    return img
end

max_stars_repo_name: kalvotom/julia-examples

max_issues_repo_name: kalvotom/julia-examples

max_forks_repo_name: kalvotom/julia-examples


text:
typealias ClientPtr Ptr{Void}
typealias PortPtr Ptr{Void}
typealias CFunPtr Ptr{Void}
typealias NFrames UInt32
typealias JACKSample Cfloat

const JACK_DEFAULT_AUDIO_TYPE = "32 bit float mono audio"

# this mirrors the struct defined in ringbuffer.h
type RingBuf
    buf::Ptr{Cchar}
    write_ptr::Csize_t
    read_ptr::Csize_t
    size::Csize_t
    size_mask::Csize_t
    mlocked::Cint
end
# add typealias for consistency with other *Ptr types
typealias RingBufPtr Ptr{RingBuf}

@enum(Option,
    # Null value to use when no option bits are needed.
    NullOption = 0x00,
    # Do not automatically start the JACK server when it is not
    # already running.  This option is always selected if
    # `JACK_NO_START_SERVER` is defined in the calling process
    # environment.
    NoStartServer = 0x01,
    # Use the exact client name requested. Otherwise, JACK
    # automatically generates a unique one, if needed.
    UseExactName = 0x02,
    # Open with optional server_name parameter.
    ServerName = 0x04,
    # Load internal client from optional
    # load_name.  Otherwise use the client_name.
    LoadName = 0x08,
    # Pass optional load_init string to the
    # jack_initialize() entry point of an internal client.
    LoadInit = 0x10,
    # pass a SessionID Token this allows the sessionmanager to identify the client again.
    SessionID = 0x20)

# useful for OR'ing options together
@compat Base.:|(l::Option, r::Option) = UInt(l) | UInt(r)

@enum(PortFlag,
  PortIsInput = 0x01,
  PortIsOutput = 0x02,
  PortIsPhysical = 0x04,
  PortCanMonitor = 0x08,
  PortIsTerminal = 0x10)

@compat Base.:|(l::PortFlag, r::PortFlag) = UInt(l) | UInt(r)

# some functions also return a -1 on failure
@enum(Status,
    Success = 0x00,
    # Overall operation failed.
    Failure = 0x01,
    # The operation contained an invalid or unsupported option.
    InvalidOption = 0x02,
    # The desired client name was not unique.  With the
    # JACKUseExactName option this situation is fatal.  Otherwise,
    # the name was modified by appending a dash and a two-digit
    # number in the range "-01" to "-99".  The
    # jack_get_client_name() function will return the exact string
    # that was used.  If the specified client_name plus these
    # extra characters would be too long, the open fails instead.
    NameNotUnique = 0x04,
    # The JACK server was started as a result of this operation.
    # Otherwise, it was running already.  In either case the caller
    # is now connected to jackd, so there is no race condition.
    # When the server shuts down, the client will find out.
    ServerStarted = 0x08,
    # Unable to connect to the JACK server.
    ServerFailed = 0x10,
    # Communication error with the JACK server.
    ServerError = 0x20,
    # Requested client does not exist.
    NoSuchClient = 0x40,
    # Unable to load internal client
    LoadFailure = 0x80,
    # Unable to initialize client
    InitFailure = 0x100,
    # Unable to access shared memory
    ShmFailure = 0x200,
    # Client's protocol version does not match
    VersionError = 0x400,
    # Backend error
    BackendError = 0x800,
    # Client zombified failure
    ClientZombie = 0x1000)

status_str(status::Status) = string(status)

# use & syntax for checking a flag, but return a boolean
@compat Base.:&{T <: Integer}(val::T, status::Status) = val & T(status) != 0

function status_str(status::Integer)
    if status == -1
        "GeneralError"
    elseif status == Int(Success)
        string(Status(status))
    else
        selected = filter(flag -> status & UInt(flag) != 0, instances(Status))
        join(selected, ", ")
    end
end

# low-level libjack wrapper functions

jack_client_open(name, options, statusref) =
    ccall((:jack_client_open, :libjack), ClientPtr,
        (Cstring, Cint, Ref{Cint}),
        name, 0, statusref)

jack_client_close(client) =
    ccall((:jack_client_close, :libjack), Cint, (ClientPtr, ), client)

jack_get_client_name(client) =
    ccall((:jack_get_client_name, :libjack), Cstring, (ClientPtr, ), client)

jack_get_sample_rate(client) =
    ccall((:jack_get_sample_rate, :libjack), NFrames, (ClientPtr, ), client)

jack_set_process_callback(client, callback, userdata) =
    ccall((:jack_set_process_callback, :libjack), Cint,
        (ClientPtr, CFunPtr, Ptr{Void}),
        client, callback, userdata)

jack_on_shutdown(client, callback, userdata) =
    ccall((:jack_on_shutdown, :libjack), Cint,
        (ClientPtr, CFunPtr, Ptr{Void}),
        client, callback, userdata)

jack_get_ports(client, portname, typename, flags) =
    ccall((:jack_get_ports, :libjack), Ptr{Cstring},
        (ClientPtr, Cstring, Cstring, Culong),
        client, portname, typename, flags)

jack_connect(client, src, dest) =
    ccall((:jack_connect, :libjack), Cint, (ClientPtr, Cstring, Cstring),
    client, src, dest)

jack_free(ptr) = ccall((:jack_free, :libjack), Void, (Ptr{Void}, ), ptr)

jack_port_register(client, portname, porttype, flags, bufsize) =
    ccall((:jack_port_register, :libjack), PortPtr,
        (ClientPtr, Cstring, Cstring, Culong, Culong),
        client, portname, porttype, flags, bufsize)

jack_port_unregister(client, port) =
    ccall((:jack_port_unregister, :libjack), Cint, (ClientPtr, PortPtr),
        client, port)

jack_port_get_buffer(port, nframes) =
    ccall((:jack_port_get_buffer, :libjack), Ptr{JACKSample},
        (PortPtr, NFrames),
        port, nframes)

jack_ringbuffer_create(bytes) =
    ccall((:jack_ringbuffer_create, :libjack), Ptr{RingBuf}, (Csize_t, ), bytes)

jack_ringbuffer_free(buf) =
    ccall((:jack_ringbuffer_free, :libjack), Void, (Ptr{RingBuf}, ), buf)

jack_ringbuffer_read(ringbuf, dest, bytes) =
    ccall((:jack_ringbuffer_read, :libjack), Csize_t,
        (Ptr{RingBuf}, Ptr{Void}, Csize_t), ringbuf, dest, bytes)

jack_ringbuffer_read_advance(ringbuf, bytes) =
    ccall((:jack_ringbuffer_read_advance, :libjack), Void,
        (Ptr{RingBuf}, Csize_t), ringbuf, bytes)

jack_ringbuffer_read_space(ringbuf) =
    ccall((:jack_ringbuffer_read_space, :libjack), Csize_t,
        (Ptr{RingBuf}, ), ringbuf)

jack_ringbuffer_write(ringbuf, src, bytes) =
    ccall((:jack_ringbuffer_write, :libjack), Csize_t,
        (Ptr{RingBuf}, Ptr{Void}, Csize_t), ringbuf, src, bytes)

jack_ringbuffer_write_space(ringbuf) =
    ccall((:jack_ringbuffer_write_space, :libjack), Csize_t,
        (Ptr{RingBuf}, ), ringbuf)


text:
"""
write_append(path::String, content::String)

Create the file if not exist. If exist write to file.

# Argument
- path: path to write file
- content: The content to write in the file
"""
function write_append(path::String, content::String)
    msg = "Project added to log successfully."
    if !isfile(path)
        open(path, "w") do io
            print(io, "$content\n")
        end
        return msg
    end
    open(path, "a") do io
        print(io, "$content\n")
    end
    return msg
end

"""
Log a project after consolidation to the given project.

# Argument
- path: path of the log file.
- prj: The consolidated project dictionary
"""
function log_project(path::String, pj::Dict)
    log_file = "projects.log"
    log_path = joinpath(path, log_file)
    metainfo = (
        "$(pj["uid"])\t$(pj["name"])\t$(pj["full_name"])\t"
        *"$(pj["created_date"])\t$(pj["profile"])"
    )
    try
        msg = write_append(log_path, metainfo)
        println(msg)
        return msg
    catch ex
        throw(ex)
    end
end


text:
using SchumakerSpline
using LinearAlgebra

x_row = [0.1,0.3]
x_col = [0.2,0.3,0.5]
ygrid = [1.0 2.0 3.0;
         18.0 -2.0 3.3]
sp = Schumaker2d(x_row,x_col,ygrid; extrapolation = (Linear, Linear))
sp(0.4, 0.3)


text:
export yin

"""
    yin(sig::Vector, sr::Int; w_len::Int = 512, f_step::Int = 256, f0_min = 100, f0_max = 500,
        harmonic_threshold = 0.1,
        difference_function::Function = difference_function_original, kwargs...) -> F0s, frame_times

Estimates the fundamental frequency (F0) of the signal `sig` using the YIN algorithm [^1].
The signal `sig` is a vector of points uniformly sampled at a rate `sr`.

## Keyword arguments
* `w_len`: size of the analysis window [samples == number of points]
* `f_step`: size of the lag between two consecutive frames [samples == number of points]
* `f0_min`: Minimum fundamental frequency that can be detected [linear frequency]
* `f0_max`: Maximum fundamental frequency that can be detected [linear frequency]
* `harmonic_threshold`: Threshold of detection. The algorithm returns the first minimum of
  the CMNDF function below this threshold.
* `diffference_function`: The difference function to be used (by default
  [`differenceFunction_original`](@ref)).

## Description
The YIN algorithm [^CheveigneYIN2002] estimates the signal's fundamental frequency `F0` by basically
looking for the period `τ0`  which minimizes the signal's autocorrelation. This
autocorrelation is calculated for signal segments (frames), composed of two windows of
length `w_len`. Each window is separated by a distance `τ`, and the idea is that the
distance which minimizes the pairwise difference between each window is considered to be the
fundamental period `τ0` of that frame.

More precisely, the algorithm first computes the cumulative mean normalized difference
function (MNDF) between two windows of a frame for several candidate periods `τ` ranging
from `τ_min=sr/f0_max` to `τ_max=sr/f0_min`. The MNDF is defined as
```math
d_t^\\prime(\\tau) = \\begin{cases}
        1 & \\text{if} ~ \\tau=0 \\\\
        d_t(\\tau)/\\left[{(\\frac 1 \\tau) \\sum_{j=1}^{\\tau} d_{t}(j)}\\right] & \\text{otherwise}
        \\end{cases}
```
where `d_t` is the difference function:
```math
d_t(\\tau) = \\sum_{j=1}^W (x_j - x_{j+\\tau})^2
```

It then refines the local minima of the MNDF using parabolic (quadratic) interpolation. This
is done by taking each minima, along with their first neighbor points, and finding the
minimum of the corresponding interpolated parabola. The MNDF minima are substituted by the
interpolation minima. Finally, the algorithm chooses the minimum with the smallest period
and with a corresponding MNDF below the `harmonic threshold`. If this doesn't exist, it
chooses the period corresponding to the global minimum. It repeats this for frames starting
at the first signal point, and separated by a distance `f_step` (frames can overlap), and
returns the vector of frequencies `F0=sr/τ0` for each frame, along with the start times of
each frame.

As a note, the physical unit of the frequency is 1/[time], where [time] is decided by the
sampling rate `sr`. If, for instance, the sampling rate is over seconds, then the frequency
is in Hertz.

[^CheveigneYIN2002]: De Cheveigné, A., & Kawahara, H. (2002). YIN, a fundamental frequency estimator for
speech and music. The Journal of the Acoustical Society of America, 111(4), 1917-1930.
"""
function yin(sig::Vector, sr::Int; w_len::Int = 512, f_step::Int = 256, f0_min = 100,
    f0_max = 500, harmonic_threshold = 0.1,
    difference_function::Function = difference_function_original, kwargs...)

    τ_min = floor(Int64, sr / f0_max)
    τ_max = floor(Int64, sr / f0_min)

    frame_times = range(1, length(sig) - w_len - τ_max, step=f_step)  # time values for start of  each analysis window
    # times = frame_times ./ eltype(sig)(sr)

    F0s = zeros(Float64, length(frame_times))

    for (i, t) in enumerate(frame_times)
        frame = sig[ (t) : (t + τ_max+ w_len-1) ]
        df = difference_function(frame, w_len, τ_max)
        cmdf = cumulative_mean_normalized_difference_function(df)
        y_refined, τ_refined, τ_indices = refine_local_minima(cmdf)
        idx_localminimum = absolute_threshold(y_refined, τ_min, τ_max, harmonic_threshold)
        τ0 = τ_refined[idx_localminimum]
        F0s[i] = sr / τ0

    end
    return  F0s, frame_times
end

"""
    difference_function_original(x, W, τmax) -> df
Computes the difference function of `x`. `W` is the window size, and `τmax` is the maximum
    period. This corresponds to equation (6) in [^1]:
```math
d_t(\\tau) = \\sum_{j=1}^W (x_j - x_{j+\\tau})^2
```
"""
function difference_function_original(x, W, τmax)
    df = zeros(eltype(x), τmax+1) #df corresponds to τ values from 0 to τ_max
    for τ in 1:τmax
        for j in 1:W
            df[τ+1] += (x[j] - x[j + τ]) ^ 2
        end
    end
    return df
end

"""
Compute cumulative mean normalized difference function (CMND), starting from the difference
    function `df`. This corresponds to equation (8) in [1]
"""
function cumulative_mean_normalized_difference_function(df)
    N = length(df)
    cmndf = df[2:end] .* range(1, N-1, step=1) ./ Float64.(cumsum(df[2:end]))
    return [1.0; cmndf]
end

"""
Returns the refined local minima of `y`, along with their refined `x` value and the
corresponding indices.  For each minimum, by the minimum of the parabola obtained by
interpolated the minimum with its first neighbors. Also returns the indices of `y`
(`x_nominal`) and the value `x` corresponding to each minima.
"""
function refine_local_minima(y)
    x_nominal = 0:length(y)-1 #nominal τ values: 0 to τ_max (cf. difference_function)
    idxs_local_minima = local_minima(y)
    if isempty(idxs_local_minima)
         @warn "No local minima found for the cumulative mean difference function. Adjusting
         the values of the minimum and maximum frequencies may fix this."
    end
    xv, yv = parabolic_interpolation(y, idxs_local_minima) #xv,yv are correction for the local minima
    y_locmin_refined = yv
    x_locmin_real= x_nominal[idxs_local_minima] .+ xv
    x_locmin_nominal = x_nominal[idxs_local_minima]
    return y_locmin_refined, x_locmin_real, x_locmin_nominal
end

"""
Applies the threshold step described in [^1]. It returns index (period) corresponding to the
first minimum below the `harmonic_threshod`. If that doesn't exist, it regurns the index of the
global minimum.
"""
function absolute_threshold(localminima, τ_min, τ_max, harmonic_threshold=0.1)
    for (idx, localminimum) in enumerate(localminima)
        if localminimum ≤ harmonic_threshold
            return idx
        end
    end
    return argmin(localminima)
end

"""
Calculates the parabolic (quadratic) interpolation for all the indixes of `y` given in
idxs_interpolate, and returns the coordinates of the respective parabola's minimum.
Assumes space all adjacent x values of `y` is 1.
"""
function parabolic_interpolation(y, idxs_interpolate)
    #separate x into triplets [x1,x2,x3].
    x1 = @view y[idxs_interpolate .- 1]
    x2 = @view y[idxs_interpolate]
    x3 = @view y[idxs_interpolate .+ 1]
    #calculate the vertix coordinates (xv, yv) for each triplet
    a = @. (x1 - 2x2 + x3)/2
    b = @. (x3 - x1) / 2
    xv = @. -b/2a
    yv = @. x2 - (b^2 /4a)
    return xv, yv
end

"""
Quickly written version for finding local minima by comparing first neighbors only.
#TODO: More detailed implementations available in Peaks.jl or Iamges.jl, maybe this should
be replaced by them.
"""
function local_minima(x)
    x1 = @view x[1:end - 2]
    x2 = @view x[2:end - 1]
    x3 = @view x[3:end]
    collect(1:length(x2))[x1 .> x2 .< x3] .+ 1
end


text:
using DemoPackagePM
using Test

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


text:
mutable struct Item{I<:Inventory, S<:Tuple, P} <: AbstractItem
    inventory::I
    sources::S
    policy::P
    position_log::Vector{Float64}
    name::String
end

function Item(inventory::Inventory, sources::Union{Assembly, Supplier}...; policy = sSPolicy(), name = "item")
    Item(inventory, sources, policy, zeros(0), name)
end

RLBase.state(e::Item) = [state(e.inventory); reduce(vcat, [state(source) for source in e.sources])]
state_size(e::Item) = state_size(e.inventory) + sum(state_size.(e.sources))
action_size(e::Item) = action_size(e.policy)*length(e.sources)
function print_state(e::Item)
    ps = [print_state(e.inventory); reduce(vcat, print_state.(e.sources))]
    return [e.name*" "*first(p) => last(p) for p in ps]
end
print_action(e::Item) = reduce(vcat, [e.name .* " " .* source.name .* " " .* print_action(e.policy) for source in e.sources])

function pull!(e::Item, quantity::Number, issuer)
    pull!(e.inventory, quantity, issuer)
end

function (e::Item)(Qs)
    for (source, Q) in zip(e.sources, Qs)
        pull!(source, Q, e.inventory)
    end
end

function dispatch!(e::Item)
    for source in e.sources
        dispatch!(source)
    end
    dispatch!(e.inventory)
end

function reward!(e::Item)
    r = 0.
    r += reward!(e.inventory)
    for source in e.sources
        r += reward!(source)
    end
    push!(e.position_log, inventory_position(e))
    return r
end

function RLBase.reset!(e::Item)
    empty!(e.position_log)
    reset!(e.inventory)
    for source in e.sources
        reset!(source)
    end
    return nothing
end

function inventory_position(e::Item)
    ip = inventory_position(e.inventory)
    for source in e.sources
        ip += inventory_position(source)
    end
    return ip
end

function children(e::Item)
    mapreduce(children, (x,y) -> tuple(x..., y...), e.sources)
end

function Base.show(io::IO, e::Item)
    print(io, "Item(", e.name,": ")
    show(io, e.inventory)
    print(io, ", ")
    for source in e.sources
        show(io, source)
        print(io, ", ")
    end
    show(io, e.policy)
    print(io, ")")
end


text:
using Documenter, ModelConstructors, SMC, Distributions

makedocs(modules = [SMC],
         clean = false,
         format = Documenter.HTML(),
         sitename = "SMC.jl",
         authors = "FRBNY-DSGE",
         linkcheck = false,
         strict = false,
         pages = Any[
                     "Home"                                   => "index.md",
                     "Using SMC"                              => "smc.md",
                     "Helper Functions"                       => "helpers.md",
                     "License"                                => "license.md"
         ],
         doctest = false # for now
)

deploydocs(
    repo = "github.com/FRBNY-DSGE/SMC.jl.git",
    target = "build",
    deps = nothing,
    devbranch = "master",
    branch = "gh-pages",
    # versions = "v#",
    # julia = "0.7",
    # osname = "osx",
    make = nothing
)


text:
using RepeatArrays
using Test

@testset "RepeatArrays.jl" begin
    a = [1. 2.
         3. 4.]

    @testset "Inner repetition" begin
        @test size(myrepeat(a, inner=(2,))) == (4, 2)
        @test size(myrepeat(a, inner=(1, 2))) == (2, 4)
        @test collect(myrepeat(a, inner=(2,))) == repeat(a, inner=(2, 1))
        @test collect(myrepeat(a, inner=(1, 2))) == repeat(a, inner=(1, 2))

        @testset "Inner rep of singleton treated as outer rep" begin
            r = myrepeat([1, 2], inner=(1, 3))
            @test r.repetitions[2] == (2, 1, 3)
        end
    end

    @testset "Outer repetition" begin
        @test size(myrepeat(a, outer=(2,))) == (4, 2)
        @test collect(myrepeat(a, outer=(2,))) == repeat(a, outer=(2, 1))
    end

    @testset "Mixed repetition" begin
        @test size(myrepeat(a, inner=(1, 3), outer=(2,))) == (4, 6)
    end

    @testset "Broadcasting" begin
        R = myrepeat([1, 2, 3], inner=(10,), outer=(2,))

        @testset "One dimension" begin
            a = myrepeat([1, 2], outer=(10,))
            b = myrepeat([.1, .2], inner=(10,))
            @test collect(a .+ b) == collect(a) .+ collect(b)

            S = myrepeat([.1, .2], inner=(5), outer=(6,))
            expected = repeat([1, 2, 3], inner=10, outer=2) .+ repeat([.1, .2], inner=5, outer=6)
            @test collect(R .+ S) == expected
        end

        @testset "Two dimensions" begin
            Q = myrepeat([1 2
                          3 4], inner=(5, 2), outer=(6, 3))
            expected = collect(R) .+ collect(Q)
            @test collect(R .+ Q) == expected
        end

        @testset "Promotion" begin
            a = myrepeat([1, 2], outer=(10,))
            b = myrepeat([1 2
                          3 4], inner=(10, 5))
            expected = collect(a) .* collect(b)
            check = myrepeat([1, 2], outer=(10, 10))
            @test collect(a .* b) == expected
            @test collect(a .* b) == collect(check) .* collect(b)

            c = myrepeat([4 5], outer=(1, 5))
            @test collect(b .* c) == collect(b) .* collect(c)
            @test collect(a .* c) == collect(a) .* collect(c)
        end

        @testset "With normal arrays" begin
            a = myrepeat([1, 2], outer=(10,))
            @test collect(a .* Array(1:20)) == collect(a) .* Array(1:20)
            @test collect(a .* Array(1:20)') == collect(a) .* Array(1:20)'
            @test collect(Array(1:20) .* a) == collect(a) .* Array(1:20)
            @test collect(Array(1:20)' .* a) == collect(a) .* Array(1:20)'
        end
    end

    @testset "Nested" begin
        a = myrepeat([1, 2], inner=(10,))
        b = myrepeat(a, outer=(2, 3))
        @test collect(b) == repeat(repeat([1, 2], inner=10), outer=(2, 3)) == repeat(collect(a), outer=(2, 3))
        c = myrepeat(b, inner=(4, 5))
        @test collect(c) == repeat(collect(b), inner=(4, 5))

        @testset "Nested broadcasting" begin
            d = myrepeat([4, 5], outer=(5,))
            e = myrepeat(d, inner=(4, 3))
            @test collect(e .+ b) == collect(e) .+ collect(b)
        end
    end
end

nothing


text: Please remove the contents of this file and replace it with your auth token provided by Quandl.com if you wish to have more than 10 downloads per day


text:
using GLMakie
xtmp = 10 .^(range(-1.5, stop=0.3,length=10))
x = vcat(-reverse(xtmp), xtmp)
y = vcat(-reverse(xtmp), xtmp)
f(z) = log(z)
fvals = [f(u+1im*v) for u in x, v in y]
fvalues = real.(fvals)
fargs = angle.(fvals)

set_theme!(theme_dark())
fig = Figure(resolution = (1200,600))
ax1 = Axis3(fig[1,1]; aspect = :data, perspectiveness = 0.5,
    elevation = 0.08, azimuth = -1.02)
ax2 = Axis3(fig[1,2]; aspect = :data, perspectiveness = 0.5,
    elevation = 0.11, azimuth = -0.45)
wireframe!(ax1, x, y, fvalues; color = :white, transparency = true)
surface!(ax2, x, y, fvalues; color = fargs,
    colormap = :diverging_rainbow_bgymr_45_85_c67_n256,
    colorrange = (-π,π), backlight = 1f0, highclip = :black)
save("wireframeSurface.png", fig)
display(fig)