text:
#using CSV
using PyPlot

include("../../SysSimExClusters/src/clusters.jl")
include("../../SysSimExClusters/src/planetary_catalog.jl")

##### To get the SNR of each planet:

snr_singles = Float64[]
snr_multis = Float64[]
snr_multis_weakest = Float64[]

#=
# No cuts at all (except only keeping confirmed and candidate planets):
###f = open("snrs_kepler_all.out", "w")
planets_keep = planet_catalog[(planet_catalog[!,:koi_disposition] .== "CONFIRMED") .| (planet_catalog[!,:koi_disposition] .== "CANDIDATE"),:]
KOI_systems = [x[1:6] for x in planets_keep[!,:kepoi_name]]
checked_bools = zeros(size(planets_keep,1))
for i in 1:length(KOI_systems)
    if checked_bools[i] == 0
        system_i = (1:length(KOI_systems))[KOI_systems .== KOI_systems[i]]
        checked_bools[system_i] .= 1
        system_snr = planets_keep[!,:koi_model_snr][system_i]
        ###println(f, system_snr)
        if length(system_snr) > 1
            append!(snr_multis, system_snr[.~ismissing.(system_snr)])
            append!(snr_multis_weakest, minimum(system_snr[.~ismissing.(system_snr)]))
        elseif length(system_snr) == 1
            append!(snr_singles, system_snr[.~ismissing.(system_snr)])
        end
    end
end
###close(f)
=#

#=
# Only applying our period and radius cuts (3-300 d, 0.5-10 R_earth):
# WARNING: actually the cuts in radius here are based on the old radius values (i.e. without revised stellar properties from Gaia DR2)
###f = open("snrs_kepler_PR_cuts.out", "w")
planets_keep = planet_catalog[(planet_catalog[!,:koi_disposition] .== "CONFIRMED") .| (planet_catalog[!,:koi_disposition] .== "CANDIDATE"),:]
planets_keep = planets_keep[(planets_keep[!,:koi_period] .> 3.) .& (planets_keep[!,:koi_period] .< 300.) .& (planets_keep[!,:koi_prad] .> 0.5) .& (planets_keep[!,:koi_prad] .< 10.) .& (.~ismissing.(planets_keep[!,:koi_prad])), :]
KOI_systems = [x[1:6] for x in planets_keep[!,:kepoi_name]]
checked_bools = zeros(size(planets_keep,1))
for i in 1:length(KOI_systems)
    if checked_bools[i] == 0
        system_i = (1:length(KOI_systems))[KOI_systems .== KOI_systems[i]]
        checked_bools[system_i] .= 1
        system_snr = planets_keep[!,:koi_model_snr][system_i]
        ###println(f, system_snr)
        if length(system_snr) > 1
            append!(snr_multis, system_snr[.~ismissing.(system_snr)])
            append!(snr_multis_weakest, minimum(system_snr[.~ismissing.(system_snr)]))
        elseif length(system_snr) == 1
            append!(snr_singles, system_snr[.~ismissing.(system_snr)])
        end
    end
end
###close(f)
=#

#=
# Only applying our stellar cuts (stellar sample from HFR2020b):
###f = open("snrs_kepler_HFR2020b_stars.out", "w")
add_param_fixed(sim_param,"min_period", 0.)
add_param_fixed(sim_param,"max_period", Inf)
add_param_fixed(sim_param,"min_radius", 0.)
add_param_fixed(sim_param,"max_radius", Inf)
@time planets_keep = keep_planet_candidates_given_sim_param(planet_catalog; sim_param=sim_param, stellar_catalog=stellar_catalog, recompute_radii=true)
KOI_systems = [x[1:6] for x in planets_keep[!,:kepoi_name]]
checked_bools = zeros(size(planets_keep,1))
for i in 1:length(KOI_systems)
    if checked_bools[i] == 0
        system_i = (1:length(KOI_systems))[KOI_systems .== KOI_systems[i]]
        checked_bools[system_i] .= 1
        system_snr = planets_keep[!,:koi_model_snr][system_i]
        ###println(f, system_snr)
        if length(system_snr) > 1
            append!(snr_multis, system_snr[.~ismissing.(system_snr)])
            append!(snr_multis_weakest, minimum(system_snr[.~ismissing.(system_snr)]))
        elseif length(system_snr) == 1
            append!(snr_singles, system_snr[.~ismissing.(system_snr)])
        end
    end
end
###close(f)
=#

#
# Using all our cuts in HFR2020b:
###f = open("snrs_kepler_HFR2020b_all_cuts.out", "w")
KOI_systems = [x[1:6] for x in planets_cleaned[!,:kepoi_name]]
checked_bools = zeros(size(planets_cleaned,1))
for i in 1:length(KOI_systems)
    if checked_bools[i] == 0
        system_i = (1:length(KOI_systems))[KOI_systems .== KOI_systems[i]]
        checked_bools[system_i] .= 1
        system_snr = planets_cleaned[!,:koi_model_snr][system_i]
        ###println(f, system_snr)
        if length(system_snr) > 1
            append!(snr_multis, system_snr[.~ismissing.(system_snr)])
            append!(snr_multis_weakest, minimum(system_snr[.~ismissing.(system_snr)]))
        elseif length(system_snr) == 1
            append!(snr_singles, system_snr[.~ismissing.(system_snr)])
        end
    end
end
###close(f)
#

##### To make some plots:

logrange(x1, x2, n) = collect(10^logx for logx in range(log10(x1), log10(x2), length=n))

# Histogram of mass ratios:
hist([snr_singles, snr_multis, snr_multis_weakest], bins=logrange(5., 1e4, 101), histtype="step", label=["Singles","Multis","Multis weakest"])
ax = gca()
ax.set_xscale("log")
xlabel("SNR", fontsize=20)
legend(fontsize=20)


text:
import JuliennedArrays
using Documenter: doctest

doctest(JuliennedArrays)

ext: jl


text:
@enum Criteria CROWDING HV # different criterium for SELECTION_binary_tournament

struct _MOMKP
    name :: String # name of the instance
    n  :: Int # number of items
    m  :: Int # number of dimensions
    P  :: Matrix{Int} # profit of items for the objectives, k=1..p, j=1..n
    W  :: Matrix{Int} # weight of items for the constraints, i=1..m, j=1..n
    ω  :: Vector{Int} # capacity of knapsacks, i=1..m
end

mutable struct Sol
    prob::_MOMKP

    x::Vector{Int} # length(x) = prob.n
    z::Vector{Int} # length(z) = 2
    w::Vector{Int}

    rank::Int
    crowding::Int
    hv_contrib::Int

    Sol(prob, x, z, w) = new(prob, x, z, w, 0, 0, 0) # only aload constructor
end

import Base.:(==), Base.:(!=)

Base.:(==)(s1::Sol, s2::Sol) = (s1.x == s1.x) && (s1.z == s2.z)
Base.:(!=)(s1::Sol, s2::Sol) = (s1.x != s1.x) || (s1.z != s2.z)

#------------ Configurations -------------#

struct Config
    size_pop_init::Int
    crossover::Function
    mutation_rate::Float64
    reparation!::Function
end


text:
#=
    This is a script to test the whole simulation chain with all standard detector types.
    The script also produces some output plots but only very basic ones, since it
    is only to test the core functionality of the package. (So not detector specific plots)
=#

outputdir = joinpath(ENV["HOME"], "tmp/test_solidstatedetectors.jl/")
mkpath(outputdir)
@info "Test output dir: $outputdir"

@info "Loading packages"
using Plots
using SolidStateDetectors; SSD = SolidStateDetectors
using LinearAlgebra

T = Float32
@info "Testing now for Float32:"

plot() # creates a plot so that the plots during the following loop pop up.

for key in [:InvertedCoax, :BEGe, :Coax, :CGD, :Hexagon]

    @info "Now test detector type: $key"

    sim = Simulation(SSD_examples[key]);
    S = SSD.get_coordinate_system(sim)

    apply_initial_state!(sim, ElectricPotential)
    p = if S == SSD.Cartesian
        plot(
            plot(sim.electric_potential, y = 0.002),
            plot(sim.q_eff_imp, y = 0.002),
            plot(sim.point_types, y = 0.002),
            size = (1000, 600), layout= (1, 3)
        )
    else
        plot(
            plot(sim.electric_potential),
            plot(sim.q_eff_imp),
            plot(sim.point_types),
            size = (1000, 600), layout= (1, 3)
        )
    end
    savefig(joinpath(outputdir, "$(key)_0_init_setup"))

    nrefs = if key in [:InvertedCoax, :Heaxgon, :CGD]
        0:3
    else
        0:1
    end
    for nref in nrefs
        SSD.update_till_convergence!(sim, ElectricPotential)
        p = if S == SSD.Cartesian
            plot(
                plot(sim.electric_potential, y = 0.002),
                plot(sim.q_eff_imp, y = 0.002),
                plot(sim.point_types, y = 0.002),
                size = (1000, 600), layout= (1, 3)
            )
        else
            plot(
                plot(sim.electric_potential),
                plot(sim.q_eff_imp),
                plot(sim.point_types),
                size = (1000, 600), layout= (1, 3)
            )
        end
        savefig(joinpath(outputdir, "$(key)_1_Electric_Potential_$(nref)_refinements"))
        if nref != nrefs[end]
            SSD.refine!(sim, ElectricPotential, (100, 100, 100), (1e-4, 1e-4, 1e-4), update_other_fields = true)
        end
        @show size(sim.electric_potential.grid)
    end

    for contact in sim.detector.contacts
        calculate_weighting_potential!(sim, contact.id, refinement_limits = key == :Coax ? missing : [0.2], verbose = true)
    end
    wpot_plots = if S != SSD.Cartesian
        [ plot(sim.weighting_potentials[contact.id]) for contact in sim.detector.contacts ]
    else
        [ plot(sim.weighting_potentials[contact.id], y = 0.002) for contact in sim.detector.contacts ]
    end
    plot( wpot_plots..., size = (1000, 1000))
    savefig(joinpath(outputdir, "$(key)_2_Weighting_Potentials"))

    calculate_electric_field!(sim)

    if S == SSD.Cylindrical
        plot(sim.electric_field, φ = 30)
        savefig(joinpath(outputdir, "$(key)_3_Electric_Field_strength"))
        plot_electric_fieldlines!(sim, φ = 30)
        savefig(joinpath(outputdir, "$(key)_3_1_Electric_Field_Lines"))
    else
        plot(sim.electric_field, y = 0)
        savefig(joinpath(outputdir, "$(key)_3_Electric_Field_strength"))
        plot_electric_fieldlines!(sim, y = 0)
        savefig(joinpath(outputdir, "$(key)_3_1_Electric_Field_Lines"))
    end

    sim.detector = SolidStateDetector(sim.detector, ADLChargeDriftModel(T = T))

    pos = if key == :InvertedCoax
        CylindricalPoint{T}[ CylindricalPoint{T}( 0.02, deg2rad(10), 0.025 ) ]
    elseif key == :CGD
        CartesianPoint{T}[ CartesianPoint{T}( 0.003, 0.003, 0.003 ) ]
    elseif key == :BEGe
        CylindricalPoint{T}[ CylindricalPoint{T}( 0.016, deg2rad(10), 0.015  ) ]
    elseif key == :Coax
        CylindricalPoint{T}[ CylindricalPoint{T}( 0.016, deg2rad(10), 0.005  ) ]
    elseif key == :Spherical
        CylindricalPoint{T}[ CylindricalPoint{T}( 0.00, deg2rad(0), 0.0  ) ]
    elseif key == :Hexagon
        CartesianPoint{T}[ CartesianPoint{T}( 0.0004, 0.0, 0.0 )]
    end
    energy_depos = T[1460]
    @assert in(pos[1], sim.detector) "Test point $(pos[1]) not inside the detector $(key)."

    event = Event(pos, energy_depos)
    simulate!(event, sim)

    plot(sim.detector)
    plot!(event.drift_paths)
    savefig(joinpath(outputdir, "$(key)_4_charge_drift"))

    plot([event.waveforms...], size = (1200, 600), lw = 1.5)
    savefig(joinpath(outputdir, "$(key)_5_induced_signals"))

end

@info "Finished testing."
@info "Test output saved in: $outputdir"


text:
module SpherePacking

using AbstractAlgebra
using SpecialFunctions
using MPMP
using WriteFilesSDPB

export testbound_sphere_packing, Nsphere_packing_2point

function spherevolume(n, r)
    sqrt(BigFloat(pi))^n / gamma(BigFloat(n)/2 + 1) * r^n
end

laguerre(k, alpha, x) = MPMP.laguerrebasis(k, alpha, x)[end]

function standard_basis(N, i, j, element = 1, symmetric = true)
    # symmetric: place element on both i,j and j,i (i = j => element on i,i)
    # non-symmetric: place element on i,j
    E = zeros(parent(element), N, N)
    E[i,j] = element
    if symmetric
        E[j,i] = element
    end
    return E
end

function Nsphere_packing_2point(n, d, r, N=2, file_path="", write_only=false; omega = BigFloat(10)^2, normalize_prep = false,kwargs...)
    if precision(BigFloat) == 256 #default, usually need at least 512
        setprecision(512)
    end
    #NOTE: we use 2d as max degree, and sort of 4d because we use u = ||x||^2 as variable
    #variables M, a_{ij,k} : for k = 0:d for i = 1:N for j=1:i
    # max -M
    # M - f_ii(0)>=0 for i=1:N
    # F(f)(0) - (vol B(r_i)^1/2 vol B(r_j)^1/2 )_{ij=1}^N >=0
    # F(f)(t) >=0 for t>=0
    # -f(w)_ij >= 0 for w>=r_i+r_j for i=1:N j=1:i (symmetric)
    # With f(x) = sum_k a_k k!/pi^k L_k^{n/2-1}(pi ||x||^2)
    # and thus F(f)(x) = sum_k a_k x^{k}
    # max -M
    # s.t.  - (vol B(r_i)^1/2 vol B(r_j)^1/2 )_{ij=1}^N + \sum_i \sum_j a_{ij,k} E_ij >= 0 (G={1}) (NxN)
    #       0 + sum_k sum_i sum_j a_{ij,k} E_{ij} x^k >= 0        (G = {1,x}) (NxN)
    #       0 - sum_k a_ijk k! pi^-k L_k^{n/2-1}(pi x) >=0          (G = {1,x - (r_i + r_j)^2}) (1x1)
    #       M - sum_k a_iik k!/pi^k L_k^{n/2-1}(0) >= 0 for i=1:N   G = {1}, 1x1

    #Initialize polynomial rings for the polynomial constraints
    F = RealField
    R, (x,) = PolynomialRing(F,["x"])
    SN = MatrixSpace(R,N,N)
    S1 = MatrixSpace(R,1,1) # for the constraints on individual elements

    # The constraint matrices
    #order: constant matrix, matrix for M, matrices for a_ijk: for k=0:2d for i=1:N for j=1:i
    # - (vol B(r_i)^1/2 vol B(r_j)^1/2 )_{ij=1}^N + \sum_i \sum_j a_{ij,k} E_ij >= 0 (G={1}) (NxN)
    M0 = [vcat(SN([R(-sqrt(spherevolume(n,r[i])*spherevolume(n,r[j]))) for i=1:N, j=1:N]), SN(R(0)),
            [ k == 0 ? SN(standard_basis(N,i,j,R(1),true)) : SN(R(0)) for k=0:2d for i=1:N for j=1:i])]
    #    0 + sum_k sum_i sum_j a_{ij,k} E_{ij} x^k >= 0        (G = {1,x}) (NxN)
    M1 = [vcat(SN(R(0)), SN(R(0)),[SN(standard_basis(N,i,j,x^k,true)) for k=0:2d for i=1:N for j=1:i])]
    #    0 - sum_k a_ijk k! pi^-k L_k^{n/2-1}(pi x) >=0       i=1:N, j=1:i   (G = {1,x - (r_i + r_j)^2}) (1x1)
    M2 = [[vcat(S1(R(0)), S1(R(0)),
        [ (r == i && s == j) ? S1(-factorial(big(k)) / BigFloat(pi)^k * laguerre(k, BigFloat(n)/2-1, BigFloat(pi)*x)) : S1(R(0)) for k=0:2d for r=1:N for s=1:r])] for i=1:N for j=1:i]
    #    0 + M - sum_k a_iik k!/pi^k L_k^{n/2-1}(0) >= 0 for i=1:N   G = {1}, 1x1
    M3 = [[vcat(S1(R(0)), S1(R(1)),
        [ r == s == i ? S1(-factorial(big(k)) / BigFloat(pi)^k * laguerre(k, BigFloat(n)/2-1, BigFloat(0))) : S1(R(0)) for k=0:2d for r=1:N for s=1:r])] for i=1:N]
    M = vcat(M0,M1,M2...,M3...) #M2 and M3 contain multiple constraints

    #sample points
    sample_points = vcat([BigFloat(0)],
      [BigFloat.(MPMP.create_sample_points_1d(2d))], #laguerre sample points
      [BigFloat.(MPMP.create_sample_points_1d(2d)).+ (r[i]+r[j])^2 for i=1:N for j=1:i],
      [[BigFloat(0)] for i=1:N])

    #polynomial weights
    G = vcat([[R(1)]],
        [[R(1),x]],
        [[R(1),x-(r[i]+r[j])^2] for i=1:N for j=1:i],
        [[R(1)] for i=1:N])

    # bases
    q = MPMP.laguerrebasis(d, BigFloat(n)/2-1, BigFloat(2*pi)*x)
    max_coef = [max(coeffs(q[i])...) for i=1:length(q)]
    q = [max_coef[i]^(-1)*q[i] for i=1:length(q)]

    #maximum degrees
    delta = vcat(0,2d,[2d for i=1:N for j=1:i],[0 for i=1:N])

    #objective = -M
    b = vcat(BigFloat(-1),[BigFloat(0) for k=0:2d for i=1:N for j=1:i])

    #get the constraints and blockinfo
    constraints = [prepareabc(M[j],G[j],q,sample_points[j],delta[j];normalize=normalize_prep) for j=1:length(G)]
    blockinfo = MPMP.get_block_info(constraints)

    if length(file_path)>0
        #if needed, write the constraints to files in SDPB format
        write_files(file_path,constraints,blockinfo,b)
    end
    if length(M) == 7#only change ordering if N = 2, and hence J = 7.
        #ordering of sdpb: pairs (3,6),(5,7), (4,1) and (2) on cores (if 4 cores).
        # In any case, not the polynomial constraints 2,3,4,5 on the same core if multiple cores
        ordering = [3,6,5,7,4,1,2]
        constraints = constraints[ordering]
        blockinfo = MPMP.get_block_info(constraints)
    end

    if !write_only
        #solve the MPMP and return the values.
        # default start value is 100. For large n,d this may need to be increased.
        x,X,y,Y,P,p,d,dual_gap,primal_obj,dual_obj = solverank1sdp(constraints, b, blockinfo;omega_p = omega,omega_d = omega, kwargs...);
    else
        return true
    end
end

function test_bound_sphere_packing(n=3,d=8)
    if precision(BigFloat) == 256
        setprecision(512)
    end

    #best lower bound known with these ratio is obtained from the crystal structure of NaCl
    cur_bound = Nsphere_packing_2point(n,d,[BigFloat(1), sqrt(BigFloat(2))-1],2)
    println(-cur_bound[end])

    NaCl_density = 0.793
    NaCl_bound = 0.813 #this is the bound reported by de Laat et al. in Upper bounds for packings of spheres of several radii.
    #The same bound is programmed by Nsphere_packing_2point, but then with sampling instead of coefficient matching.
    println("Compare to the density of NaCL: $NaCl_density (Current bound: $NaCl_bound)")
end

end # end module

max_stars_repo_name: nanleij/Clustered-Low-Rank-SDP-solver

max_issues_repo_name: nanleij/Clustered-Low-Rank-SDP-solver

max_forks_repo_name: nanleij/Clustered-Low-Rank-SDP-solver


text:
# # Spin-weighted spherical harmonics
# This example plays with analysis of:
# ```math
# f(r) = e^{{\rm i} k\cdot r},
# ```
# for some $k\in\mathbb{R}^3$ and where $r\in\mathbb{S}^2$, using spin-$0$ spherical harmonics.
#
# It applies ð, the spin-raising operator,
# both on the spin-$0$ coefficients as well as the original function,
# followed by a spin-$1$ analysis to compare coefficients.
#
# For the storage pattern of the arrays, please consult the
# [documentation](https://MikaelSlevinsky.github.io/FastTransforms).

using FastTransforms, LinearAlgebra

# The colatitudinal grid (mod $\pi$):
N = 10
θ = (0.5:N-0.5)/N

# The longitudinal grid (mod $\pi$):
M = 2*N-1
φ = (0:M-1)*2/M

# Our choice of $k$ and angular parametrization of $r$:
k = [2/7, 3/7, 6/7]
r = (θ,φ) -> [sinpi(θ)*cospi(φ), sinpi(θ)*sinpi(φ), cospi(θ)]

# On the tensor product grid, our function samples are:
F = [exp(im*(k⋅r(θ,φ))) for θ in θ, φ in φ]

# We precompute a spin-$0$ spherical harmonic--Fourier plan:
P = plan_spinsph2fourier(F, 0)

# And an FFTW Fourier analysis plan on $\mathbb{S}^2$:
PA = plan_spinsph_analysis(F, 0)

# Its spin-$0$ spherical harmonic coefficients are:
U⁰ = P\(PA*F)

# We can check its $L^2(\mathbb{S}^2)$ norm against an exact result:
norm(U⁰) ≈ sqrt(4π)

# Spin can be incremented by applying ð, either on the spin-$0$ coefficients:
U¹c = zero(U⁰)
for n in 1:N-1
    U¹c[n, 1] = sqrt(n*(n+1))*U⁰[n+1, 1]
end
for m in 1:M÷2
    for n in 0:N-1
        U¹c[n+1, 2m] = -sqrt((n+m)*(n+m+1))*U⁰[n+1, 2m]
        U¹c[n+1, 2m+1] = sqrt((n+m)*(n+m+1))*U⁰[n+1, 2m+1]
    end
end

# or on the original function through analysis with spin-$1$ spherical harmonics:
F = [-(k[1]*(im*cospi(θ)*cospi(φ) + sinpi(φ)) + k[2]*(im*cospi(θ)*sinpi(φ)-cospi(φ)) - im*k[3]*sinpi(θ))*exp(im*(k⋅r(θ,φ))) for θ in θ, φ in φ]

# We change plans with spin-$1$ now and reanalyze:
P = plan_spinsph2fourier(F, 1)
PA = plan_spinsph_analysis(F, 1)
U¹s = P\(PA*F)

# Finally, we check $L^2(\mathbb{S}^2)$ norms against another exact result:
norm(U¹c) ≈ norm(U¹s) ≈ sqrt(8π/3*(k⋅k))


text:
include("../src/numder_module.jl")
include("../../../../std/util/src/util_module.jl")

# =========================================================================== #
# Test
# =========================================================================== #
module numdertest
using ..numder
using ..util
using LinearAlgebra

using PyPlot
PyPlot.pygui(true)

function test_step_size()
    x0 = [0.0, 1e-15, 1e-16, 1, 2, 10]
    h = map( numder.step_size, x0 )
end

function numdiff_fwbw()

    f = sin
    x0 = collect( -pi:0.1:pi )
    g_theo = cos.(x0)

    g_fw = map(x -> numder.numdiff_fw(f, x), x0)
    g_bw = map(x -> numder.numdiff_bw(f, x), x0)

    g_error_fw = g_fw .- g_theo
    g_error_bw = g_bw .- g_theo

    PyPlot.figure()
    PyPlot.grid()
    PyPlot.plot(x0, g_fw, marker=:o)
    PyPlot.plot(x0, g_theo, marker=:*)

    PyPlot.figure()
    PyPlot.grid()
    PyPlot.plot(x0, g_error_fw, marker=:o)
    PyPlot.plot(x0, g_error_bw, marker=:o)

end

function numdiff_central()

    f = sin
    x0 = collect( -pi:0.1:pi )
    g_theo = cos.(x0)

    g_c = map(x -> numder.numdiff_central(f, x), x0)

    g_error_c = g_c .- g_theo

    PyPlot.figure()
    PyPlot.grid()
    PyPlot.plot(x0, g_c, marker=:o)
    PyPlot.plot(x0, g_theo, marker=:*)

    PyPlot.figure()
    PyPlot.grid()
    PyPlot.plot(x0, g_error_c, marker=:o)

end

function gradient_fw()

    f(x) = x[1].^2 + x[2].^3

    x_vec = collect( -pi:0.1:pi )
    y_vec = collect( -pi:0.1:pi )

    (x_mat, y_mat) = util.meshgrid( x_vec, y_vec )
    x0 = [x_mat[:] y_mat[:]]
    g_theo_fcn(x) = [2 * x[1], 3 * x[2].^2]

    g_theo = mapslices( g_theo_fcn, x0, dims = 2)
    g_fw = mapslices( x -> numder.gradient_fw(f, x), x0, dims = 2)

    g_error = mapslices( norm, g_fw - g_theo, dims = 2 )
    g_error_mat = reshape(g_error, (length(x_vec), length(y_vec)) )

    PyPlot.figure()
    PyPlot.grid()
    PyPlot.plot_surface(x_vec, y_vec, g_error_mat)
    # PyPlot.plot(g_error)

end

function jacobian_fw()

    # f(x) = [x[1].^2  x[2].^3]
    f(x) = [cos(x[1])  x[1].^2 + x[2].^2   sin(x[2])] # R^2 -> R^3

    x_vec = collect( -pi:0.1:pi )
    y_vec = collect( -pi:0.1:pi )
    (x_mat, y_mat) = util.meshgrid( x_vec, y_vec )
    x0 = [x_mat[:] y_mat[:]]

    idx = 2
    #=
    g_theo_fcn(x) =
        [
        2 * x[1]    0.0;
        0.0         3 * x[2].^2
        ]
    =#
    g_theo_fcn(x) =
        [
            -sin(x[1])      0.0
            2 * x[1]        2 * x[2]
            0.0             cos(x[2])
        ]
    g_theo = g_theo_fcn( x0[idx, :] )
    g_fw = numder.jacobian_fw(f, x0[idx, :] )
    g_error = g_fw - g_theo

end

function test_directional_diff_fw()

    f(x) = x[1]^2 + x[2]^2

    x0 = [1.0, 2]
    d = [3.0, 6]
    d = d ./ norm(d)

    dfdd = numder.directional_diff_fw(f, x0, d)

    g = numder.gradient_fw(f, x0)
    df_theo = sum(g .* d)

    println(dfdd)

    println(df_theo)

end

end # numdertest

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
using Test, ParametricOptInterface, MathOptInterface, GLPK, Ipopt

const POI = ParametricOptInterface
const MOI = MathOptInterface
const MOIU = MathOptInterface.Utilities

const ATOL = 1e-4

include("production_problem_test.jl")
include("basic_tests.jl")
include("quad_tests.jl")


text:
using Requests
using Blosc
using HttpCommon
Blosc.set_num_threads(div(Sys.CPU_CORES,2))
# import Requests: get

const BOSS_TOKEN = ENV["INTERN_TOKEN"]

headers = Dict("Authorization"          => "Token $(BOSS_TOKEN)",
                        "content-type"  => "application/blosc",
                        "Accept"        => "application/blosc",
                        "Date"          => Dates.format(now(Dates.UTC), Dates.RFC1123Format))

a = rand(UInt8, 16, 512,512)
data = Blosc.compress( a )

resp = Requests.post(URI("https://api.theboss.io/v0.8/cutout/jingpengw_test"*
                "/test/em/0/0:512/0:512/32:48/");
                data = data,
                headers = headers)

@show statuscode(resp)
@show resp

resp = Requests.get(URI("https://api.theboss.io/v0.8/cutout/jingpengw_test"*
                "/test/em/0/0:512/0:512/32:48/");
                headers = headers)

data = resp.data

b = Blosc.decompress(UInt8, data)
b = reshape(b, (16, 512,512))

n = 0
for i in eachindex(a)
    if a[i] != b[i]
        n+=1
    end
end

(a.==b)
@assert all(a.==b)
a==b
# data = reshape(data,(5,6,4))

length(data)

@show resp


text:
using KernelAbstractions, Test, CUDAapi
if has_cuda_gpu()
    using CuArrays, CUDAdrv
    CuArrays.allowscalar(false)
end

@testset "Error propagation" begin
    let event = Event(()->error(""))
        @test_throws TaskFailedException wait(event)
    end

    let event = Event(error, "")
        @test_throws CompositeException wait(MultiEvent(event))
    end

    let event = Event(error, "")
        event = Event(wait, MultiEvent(event))
        @test_throws TaskFailedException wait(event)
    end

    let event = Event(error, "")
        event = Event(()->nothing, dependencies=event)
        @test_throws TaskFailedException wait(event)
    end
end

if has_cuda_gpu()
    barrier = Base.Threads.Event()
    cpu_event = Event(wait, barrier)

    KernelAbstractions.unsafe_wait(CUDA(), cpu_event) # Event edge on CuDefaultStream
    gpu_event = Event(CUDA()) # Event on CuDefaultStream

    notify(barrier)
    wait(gpu_event)
end


text:
using Boltzmann
using Base.Test

X = rand(1000, 2000)

function rbm_smoke_test()
    model = BaseRBM(Units(1000), Units(500))
    fit(model, X)
end

function gaussian_rbm_smoke_test()
    model = BaseRBM(Units(1000, :gaussian), Units(500))
    fit(model, X)
end

rbm_smoke_test()
gaussian_rbm_smoke_test()


text:
import EventHistory
##reload("EventHistory")

##################################################
### EventClass
##################################################

## Right-censored event-time
stop   = [2,3,3];
status = [false,true,true];
e1     = EventHistory.Event(stop,status)

EventHistory.Time(e1)
EventHistory.Status(e1)

## Right-censored+left truncation
start = [0,1,2];
e2    = EventHistory.Event(start,stop,status)

EventHistory.Entry(e2)

## Competing risks model
cause = [0,2,1];
e3    = EventHistory.Event(start,stop,cause)

EventHistory.Status(e3)
EventHistory.Cause(e3)

## Interval censoring
right = [2,3,Inf];
left  = [1,-Inf,1];
e4    = EventHistory.Event(left,right,"interval")

## Formula syntax (see also below)
d = DataFrames.DataFrame(start=start,stop=stop,status=status);
EventHistory.Event([:stop,:status],d)
EventHistory.Event([:start,:stop,:status],d)

##################################################
### Cox regression
##################################################
using RDatasets
ovarian = dataset("survival", "ovarian");
ovarian[:Group] = ovarian[:Rx].-1;
S = EventHistory.Event([:FUTime,:FUStat], ovarian);
ovarian[:S] = EventHistory.Event([:FUTime,:FUStat], ovarian);

mm = EventHistory.phreg(@formula(S~Age+Group), ovarian)

## Prediction
EventHistory.predict(mm, surv=false, X=[0 0]); ## Baseline
s56 = EventHistory.predict(mm, X=[56 1], order=true); ## Survival probabilities age 40, group 1
EventHistory.predict(mm, X=[56 0], time=[100,400,600]); ## ... at time 100,400,600
EventHistory.predict(mm, X=[56 1; 56 0], time=[600,100,400]) ## ... both groups

using Winston
plot(s56[:,1], s56[:,2])

s = EventHistory.predict(mm, X=[56 1; 56 0], order=true)
pr = DataFrame(Time=[s[:,1];s[:,1]], S=[s[:,2];s[:,3]], Group=repeat(["Group1","Group2"], inner=size(s,1)))

using Gadfly
p = Gadfly.plot(pr, x="Time", y="S", color="Group",
                Geom.step, Geom.point,
                Guide.ylabel("Survival probability"), Guide.title("Age 56"))

draw(PNG("surv.png",7inch,7inch),p)

##################################################
### Cox regression with left truncation
##################################################
d = DataFrame(start=[1,2,5,2,1,7,3,4,8,8],
              stop=[2,3,6,7,8,9,9.5,10,14,17],
              event=[1,1,1,1,1,1,1,0,0,0],
              x=[1,0,0,1,0,1,1,1,0,0],
              z=[1.0,0,2.0,0,3.0,0,4.0,0,5.0,0]);
d[:S] = EventHistory.Event([:start,:stop,:event], d);
mm = EventHistory.phreg(@formula(S~x), d)

mm = EventHistory.phreg(@formula(S~x+z), d)


text:
# ---
# title: 1106. Parsing A Boolean Expression
# id: problem1106
# author: Indigo
# date: 2022-03-08
# difficulty: Hard
# categories: String
# link: <https://leetcode.com/problems/parsing-a-boolean-expression/description/>
# hidden: true
# ---
#
# Return the result of evaluating a given boolean `expression`, represented as a
# string.
#
# An expression can either be:
#
#   * `"t"`, evaluating to `True`;
#   * `"f"`, evaluating to `False`;
#   * `"!(expr)"`, evaluating to the logical NOT of the inner expression `expr`;
#   * `"&(expr1,expr2,...)"`, evaluating to the logical AND of 2 or more inner expressions `expr1, expr2, ...`;
#   * `"|(expr1,expr2,...)"`, evaluating to the logical OR of 2 or more inner expressions `expr1, expr2, ...`
#
#
#
# **Example 1:**
#
#
#
#     Input: expression = "!(f)"
#     Output: true
#
#
# **Example 2:**
#
#
#
#     Input: expression = "|(f,t)"
#     Output: true
#
#
# **Example 3:**
#
#
#
#     Input: expression = "&(t,f)"
#     Output: false
#
#
# **Example 4:**
#
#
#
#     Input: expression = "|(&(t,f,t),!(t))"
#     Output: false
#
#
#
#
# **Constraints:**
#
#   * `1 <= expression.length <= 20000`
#   * `expression[i]` consists of characters in `{'(', ')', '&', '|', '!', 't', 'f', ','}`.
#   * `expression` is a valid expression representing a boolean, as given in the description.
#
#
## @lc code=start
using LeetCode

function parse_bool_expr(expr::String)
    operations = Char[]
    operands = Char[]
    for ch in expr
        if ch == '!' || ch == '&' || ch == '|'
            push!(operations, ch)
        elseif ch == ')'
            cur_oprs = Char[]
            while operands[end] != '('
                push!(cur_oprs, pop!(operands))
            end
            pop!(operands)
            op = pop!(operations)
            rs = true
            if op == '|'
                rs = any(==('t'), cur_oprs)
            elseif op == '&'
                rs = all(==('t'), cur_oprs)
            else
                rs = cur_oprs[1] != 't'
            end
            push!(operands, rs ? 't' : 'f')
        elseif ch != ','
            push!(operands, ch)
        end
    end
    return operands[1] == 't' ? true : false
end
## @lc code=end


text:
include("$(pwd())/startup.jl")

fp = "/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/H3ORF400Vpp_final/400Vpp_all/" # All files in this path will be processed
fp = "/media/wiebke/Extreme SSD/PSM_vs_PTR3/Data/PTR3data/FR_NH4_bCP_2/"

fp = "/media/wiebke/Extreme SSD/PSM_vs_PTR3/Data/apiTOFdata/CLOUD10/run1734_02/"

filefilterRegexp = r"\.h5$"
rf = "$(fp)2018-05-08-18h05m06.h5"  # The mass scale from this file defines the mass scale of all
rf = "$(fp)2019-08-01-07h32m25_FR.h5"

rf = "$(fp)APi4_Data_2015.10.20-13h34m44s.h5"

#masslist = MasslistFunctions.loadMasslist("/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/H3ORF400Vpp_final/masslists2combine/combined.csv")
#masslist = MasslistFunctions.loadMasslist("$(fp)results/masslist-NOV2020.csv")

masslist = MasslistFunctions.loadMasslist("$(fp)MassList_NO3-_AP_03-01-20.csv")

#cr = [29 371]
cr = [29 76]

# alternatively: use an auto generated masslist
# masslistMasses, masslistElements, masslistElementsMasses, masslistCompositions = createMassList(C=0:20, O=0:20, N=0:1, allowRadicals=false) #
(masslistMasses, masslistElements, masslistElementsMasses, masslistCompositions) = masslist
s = (masslistMasses.>0) .& ( masslistMasses.<680)
masslistMasses = masslistMasses[s]
masslistCompositions = masslistCompositions[s,:]

####################### END OF SETTINGS ###############################################################

####################### Processing sequence ###########################################################

correctMassScaleAndExtractSumSpec(
    fp,
    masslistMasses,
    masslistElements,
    masslistElementsMasses,
    masslistCompositions,
    rf,
    cr,
    filefilterRegexp=filefilterRegexp,
    onlyUseAverages = false,
    plotControlMass = true,
    firstNFiles=0,
    lastNFiles = 0,
    filePrecaching = false,
    openWholeFile = true
    )

baselineAndPeakshape(
    fp,
    peakshapeRegions=8,
    peakshapeQuantileValue = 0.1,
    peakfindingNoiseThresholdValue = 25
    )
mtrx = deconvolute(
    fp,
    calcTransposed = true
    )


text:
using Base: operator_associativity

# Order of tests is intended to reflect dependency order of functionality

include("recursive.jl")

include("DataLayouts/data1d.jl")
include("DataLayouts/data2d.jl")
include("DataLayouts/data1dx.jl")
include("DataLayouts/data2dx.jl")

include("Geometry/geometry.jl")
include("Geometry/axistensors.jl")

include("Meshes/rectangle.jl")
include("Meshes/cubedsphere.jl")

include("Topologies/rectangle.jl")
include("Topologies/cubedsphere.jl")

include("quadrature.jl")
include("spaces.jl")
include("field.jl")
include("spectraloperators.jl")
include("spectralspaces_opt.jl")
include("diffusion2d.jl")
include("fdspaces.jl")
include("fdspaces_opt.jl")
include("fdspaces_opt_examples.jl")
include("fielddiffeq.jl")
include("hybrid2d.jl")
include("hybrid3d.jl")
include("hybrid_opt.jl")
include("remapping.jl")
#spheres.jl")
include("sphere_geometry.jl")
include("sphere_metric_terms.jl")
include("sphere_gradient.jl")
include("sphere_divergence.jl")
include("sphere_curl.jl")
include("sphere_diffusion.jl")
include("sphere_diffusion_vec.jl")
include("sphere_hyperdiffusion.jl")
include("sphere_hyperdiffusion_vec.jl")

if "CUDA" in ARGS
    include("gpu/cuda.jl")
    include("gpu/data.jl")
end


text:
const c = Center()
const f = Face()

const XBoundedGrid = AbstractGrid{<:Any, <:Bounded}
const YBoundedGrid = AbstractGrid{<:Any, <:Any, <:Bounded}
const ZBoundedGrid = AbstractGrid{<:Any, <:Any, <:Any, <:Bounded}

const XYBoundedGrid = AbstractGrid{<:Any, <:Bounded, <:Bounded}
const XZBoundedGrid = AbstractGrid{<:Any, <:Bounded, <:Any, <:Bounded}
const YZBoundedGrid = AbstractGrid{<:Any, <:Any, <:Bounded, <:Bounded}

const XYZBoundedGrid =  AbstractGrid{<:Any, <:Bounded, <:Bounded, <:Bounded}

#####
##### Exterior node and peripheral node
#####

"""
    inactive_cell(i, j, k, grid)

Return `true` when the tracer cell at `i, j, k` is lies outside the "active domain" of
the grid in `Bounded` directions. Otherwise, return `false`.
"""
@inline inactive_cell(i, j, k, grid) = false
@inline inactive_cell(i, j, k, grid::XBoundedGrid) = (i < 1) | (i > grid.Nx)
@inline inactive_cell(i, j, k, grid::YBoundedGrid) = (j < 1) | (j > grid.Ny)
@inline inactive_cell(i, j, k, grid::ZBoundedGrid) = (k < 1) | (k > grid.Nz)

@inline inactive_cell(i, j, k, grid::XYBoundedGrid) = (i < 1) | (i > grid.Nx) | (j < 1) | (j > grid.Ny)
@inline inactive_cell(i, j, k, grid::XZBoundedGrid) = (i < 1) | (i > grid.Nx) | (k < 1) | (k > grid.Nz)
@inline inactive_cell(i, j, k, grid::YZBoundedGrid) = (j < 1) | (j > grid.Ny) | (k < 1) | (k > grid.Nz)

@inline inactive_cell(i, j, k, grid::XYZBoundedGrid) = (i < 1) | (i > grid.Nx) |
                                                       (j < 1) | (j > grid.Ny) |
                                                       (k < 1) | (k > grid.Nz)

"""
    inactive_node(LX, LY, LZ, i, j, k, grid)

Return `true` when the location `(LX, LY, LZ)` is "inactive" and thus not directly
associated with an "active" cell.

For `Face` locations, this means the node is surrounded by `inactive_cell`s:
the interfaces of "active" cells are _not_ `inactive_node`.

For `Center` locations, this means the direction is `Bounded` and that the
cell or interface centered on the location is completely outside the active
region of the grid.
"""
@inline inactive_node(LX, LY, LZ, i, j, k, grid) = inactive_cell(i, j, k, grid)

@inline inactive_node(::Face, LY, LZ, i, j, k, grid) = inactive_cell(i, j, k, grid) & inactive_cell(i-1, j, k, grid)
@inline inactive_node(LX, ::Face, LZ, i, j, k, grid) = inactive_cell(i, j, k, grid) & inactive_cell(i, j-1, k, grid)
@inline inactive_node(LX, LY, ::Face, i, j, k, grid) = inactive_cell(i, j, k, grid) & inactive_cell(i, j, k-1, grid)

@inline inactive_node(::Face, ::Face, LZ, i, j, k, grid) = inactive_node(c, f, c, i, j, k, grid) & inactive_node(c, f, c, i-1, j, k, grid)
@inline inactive_node(::Face, LY, ::Face, i, j, k, grid) = inactive_node(c, c, f, i, j, k, grid) & inactive_node(c, c, f, i-1, j, k, grid)
@inline inactive_node(LX, ::Face, ::Face, i, j, k, grid) = inactive_node(c, f, c, i, j, k, grid) & inactive_node(c, f, c, i, j, k-1, grid)

@inline inactive_node(::Face, ::Face, ::Face, i, j, k, grid) = inactive_node(c, f, f, i, j, k, grid) & inactive_node(c, f, f, i-1, j, k, grid)

"""
    peripheral_node(LX, LY, LZ, i, j, k, grid)

Return `true` when the location `(LX, LY, LZ)`, is _either_ inactive or
lies on the boundary between inactive and active cells in a `Bounded` direction.
"""
@inline peripheral_node(LX, LY, LZ, i, j, k, grid) = inactive_cell(i, j, k, grid)

@inline peripheral_node(::Face, LY, LZ, i, j, k, grid) = inactive_cell(i, j, k, grid) | inactive_cell(i-1, j, k, grid)
@inline peripheral_node(LX, ::Face, LZ, i, j, k, grid) = inactive_cell(i, j, k, grid) | inactive_cell(i, j-1, k, grid)
@inline peripheral_node(LX, LY, ::Face, i, j, k, grid) = inactive_cell(i, j, k, grid) | inactive_cell(i, j, k-1, grid)

@inline peripheral_node(::Face, ::Face, LZ, i, j, k, grid) = peripheral_node(c, f, c, i, j, k, grid) | peripheral_node(c, f, c, i-1, j, k, grid)
@inline peripheral_node(::Face, LY, ::Face, i, j, k, grid) = peripheral_node(c, c, f, i, j, k, grid) | peripheral_node(c, c, f, i-1, j, k, grid)
@inline peripheral_node(LX, ::Face, ::Face, i, j, k, grid) = peripheral_node(c, f, c, i, j, k, grid) | peripheral_node(c, f, c, i, j, k-1, grid)

@inline peripheral_node(::Face, ::Face, ::Face, i, j, k, grid) = peripheral_node(c, f, f, i, j, k, grid) | peripheral_node(c, f, f, i-1, j, k, grid)

"""
    boundary_node(LX, LY, LZ, i, j, k, grid)

Return `true` when the location `(LX, LY, LZ)` lies on a boundary.
"""
@inline boundary_node(LX, LY, LZ, i, j, k, grid) = peripheral_node(LX, LY, LZ, i, j, k, grid) & !inactive_node(LX, LY, LZ, i, j, k, grid)


text:
## info to read data
# datafile
fname_dat = "GAB_3_lines.txt"
# electronics file
fname_specs_halt = "electronics_halt.jl"
# column numbers from hdr file
X, Y, Z = 2, 70, 1
fid = 3
linenum = 43
frame_height = 4
frame_dx = 71
frame_dy = 72
frame_dz = 73
HM_Z = [20, 42]
LM_Z = [52+1, 69] # first LM gate dropped, also noise and times in electronics_halt
relerror = false
units = 1e-12
##
using transD_GP
soundings = transD_GP.SkyTEM1DInversion.read_survey_files(fname_dat = fname_dat,
						             fname_specs_halt = fname_specs_halt,
						             LM_Z             = LM_Z,
									 HM_Z             = HM_Z,
									 frame_height     = frame_height,
									 frame_dz         = frame_dz,
									 frame_dy         = frame_dy,
									 frame_dx         = frame_dx,
									 X                = X,
									 Y                = Y,
                                     Z                = Z,
									 fid              = fid,
									 units            = units,
									 relerror         = relerror,
									 linenum          = linenum,
                                     startfrom        = 1,
									 skipevery        = 10,
									 dotillsounding   = nothing,
									 makesounding     = true)


text:
struct SnpBitMatrix{T} <: AbstractMatrix{T}
    B1::BitMatrix
    B2::BitMatrix
    model::Union{Val{1}, Val{2}, Val{3}}
    center::Bool
    scale::Bool
    μ::Vector{T}
    σinv::Vector{T}
    storagev1::Vector{T}
    storagev2::Vector{T}
end

AbstractSnpBitMatrix = Union{SnpBitMatrix, SubArray{T, 1, SnpBitMatrix{T}},
    SubArray{T, 2, SnpBitMatrix{T}}} where T

"""
    SnpBitMatrix{T}(s; model=ADDITIVE_MODEL, center=false, scale=false)

Store an `AbstractSnpArray` `s` as two `BitMatrix`es to simplify linear algebra.

# Arguments
- s: an `AbstractSnpArray`.
- model: one of `ADDITIVE_MODEL`(default), `DOMINANT_MODEL`, `RECESSIVE_MODEL`.
- center: whether to center (default: false).
- scale: whether to scale to standard deviation 1 (default: false).
"""
function SnpBitMatrix{T}(
    s::AbstractSnpArray;
    model = ADDITIVE_MODEL,
    center::Bool = false,
    scale::Bool = false) where T <: AbstractFloat
    Base.depwarn("`SnpBitMatrix{T}` is deprecated and will be removed soon. " *
        "Please use `SnpLinAlg{T}` instead.", nothing)
    if model == ADDITIVE_MODEL
        B1 = s .≥ 0x02
        B2 = s .≥ 0x03
        if center || scale
            μ = Vector{T}(undef, size(s, 2))
            μ[:] = mean(s; dims=1, model=ADDITIVE_MODEL)
        else
            μ = T[]
        end
        if scale
            σinv = Vector{T}(undef, size(s, 2))
            @inbounds @simd for j in 1:size(s, 2)
                σinv[j] = sqrt(μ[j] * (1 - μ[j] / 2))
                σinv[j] = σinv[j] > 0 ? inv(σinv[j]) : one(T)
            end
        else
            σinv = T[]
        end
    elseif model == DOMINANT_MODEL
        B1 = s .≥ 0x02
        B2 = falses(0, 0)
        if center || scale
            μ = Vector{T}(undef, size(s, 2))
            μ[:] = mean(s; dims=1, model=DOMINANT_MODEL)
        else
            μ = T[]
        end
        if scale
            σinv = Vector{T}(undef, size(s, 2))
            @inbounds @simd for j in 1:size(s, 2)
                σinv[j] = sqrt(μ[j] * (1 - μ[j]))
                σinv[j] = σinv[j] > 0 ? inv(σinv[j]) : one(T)
            end
        else
            σinv = T[]
        end
    elseif model == RECESSIVE_MODEL
        B1 = s .== 0x03
        B2 = falses(0, 0)
        if center || scale
            μ = Vector{T}(undef, size(s, 2))
            μ[:] = mean(s; dims=1, model=RECESSIVE_MODEL)
        else
            μ = T[]
        end
        if scale
            σinv = Vector{T}(undef, size(s, 2))
            @inbounds @simd for j in 1:size(s, 2)
                σinv[j] = sqrt(μ[j] * (1 - μ[j]))
                σinv[j] = σinv[j] > 0 ? inv(σinv[j]) : one(T)
            end
        else
            σinv = T[]
        end
    else
        throw(ArgumentError("unrecognized model $model"))
    end
    storagev1 = Vector{T}(undef, size(s, 1))
    storagev2 = Vector{T}(undef, size(s, 2))
    SnpBitMatrix{T}(B1, B2, model, center, scale, μ, σinv, storagev1, storagev2)
end

Base.size(bm::SnpBitMatrix) = size(bm.B1)
Base.size(bm::SnpBitMatrix, k::Integer) = size(bm.B1, k)

function Base.getindex(s::SnpBitMatrix{T}, i::Int, j::Int) where T
    x = s.model == ADDITIVE_MODEL ?
        T(getindex(s.B1, i, j) + getindex(s.B2, i, j)) : T(getindex(s.B1, i, j))
    s.center && (x -= s.μ[j])
    s.scale && (x *= s.σinv[j])
    return x
end

eltype(bm::SnpBitMatrix) = eltype(bm.μ)
issymmetric(bm::SnpBitMatrix) = issymmetric(bm.B2) && issymmetric(bm.B1)

"""
    LinearAlgebra.mul!(out, s::SnpBitMatrix, v)

In-place matrix-vector multiplication.
"""
function mul!(
    out::AbstractVector{T},
    s::SnpBitMatrix{T},
    v::AbstractVector{T}) where T <: AbstractFloat
    @assert length(out) == size(s, 1) && length(v) == size(s, 2)
    if s.scale
        s.storagev2 .= s.σinv .* v
        w = s.storagev2
    else
        w = v
    end
    if s.model == ADDITIVE_MODEL
        mul!(out, s.B1, w)
        mul!(s.storagev1, s.B2, w)
        out .+= s.storagev1
    else
        mul!(out, s.B1, w)
    end
    if s.center
        return out .-= dot(s.μ, w)
    else
        return out
    end
end

"""
    LinearAlgebra.mul!(out, s::Union{Transpose{T, SnpBitMatrix{T}}, Adjoint{T, SnpBitMatrix{T}}}, v)

In-place matrix-vector multiplication, with transposed `SnpBitMatrix`.
"""
function mul!(
    out::AbstractVector{T},
    st::Union{Transpose{T, SnpBitMatrix{T}}, Adjoint{T, SnpBitMatrix{T}}},
    v::AbstractVector{T}) where T <: AbstractFloat
    s = st.parent
    @assert length(out) == size(st, 1) && length(v) == size(st, 2)
    if s.model == ADDITIVE_MODEL
        mul!(out, transpose(s.B1), v)
        mul!(s.storagev2, transpose(s.B2), v)
        out .+= s.storagev2
    else
        mul!(out, transpose(s.B1), v)
    end
    if s.center
        out .-= sum(v) .* s.μ
    end
    if s.scale
        return out .*= s.σinv
    else
        return out
    end
end

"""
    Base.copyto!(v, s)

Copy SnpBitMatrix `s` to numeric vector or matrix `v`. If `s` is centered/scaled,
`v` will be centered/scaled using precomputed column mean `s.μ` and inverse std
`s.σinv`.
"""
function Base.copyto!(
    v::AbstractVecOrMat{T},
    s::AbstractSnpBitMatrix
    ) where T <: AbstractFloat
    m, n = size(s, 1), size(s, 2)
    @inbounds for j in 1:n
        @simd for i in 1:m
            v[i, j] = s[i, j]
        end
    end
    return v
end

"""
    Base.convert(t, s)

Convert a SnpBitMatrix `s` to a numeric vector or matrix of same shape as `s`.
If `s` is centered/scaled, `v` will be centered/scaled using precomputed column
mean `s.μ` and inverse std `s.σinv`.

# Arguments
- `t::Type{AbstractVecOrMat{T}}`: Vector or matrix type.
"""
Base.convert(::Type{T}, s::AbstractSnpBitMatrix) where T <: Array = T(s)
Array{T,N}(s::AbstractSnpBitMatrix) where {T,N} =
    copyto!(Array{T,N}(undef, size(s)), s)


text:
using BinaryProvider # requires BinaryProvider 0.3.0 or later

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(prefix, String["libmp3lame"], :libmp3lame),
]

# Download binaries from hosted location
bin_prefix = "https://github.com/SimonDanisch/LAMEBuilder/releases/download/3.100.0"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:aarch64, :glibc) => ("$bin_prefix/liblame.v3.100.0.aarch64-linux-gnu.tar.gz", "1be606f6b7cb9a015cfcd22c3cf58ef8752b46da7edc25a5412b6b9fc34697f4"),
    Linux(:aarch64, :musl) => ("$bin_prefix/liblame.v3.100.0.aarch64-linux-musl.tar.gz", "651340c2e3e2b0231d04f023bea996285c9e396660d4f47a2f2c9a8268b8e437"),
    Linux(:armv7l, :glibc, :eabihf) => ("$bin_prefix/liblame.v3.100.0.arm-linux-gnueabihf.tar.gz", "1729adb511dad94c2dc2e358438fdc90c213ad97655752c6ec66251ac063fa95"),
    Linux(:armv7l, :musl, :eabihf) => ("$bin_prefix/liblame.v3.100.0.arm-linux-musleabihf.tar.gz", "0c119cdb397603005936095a51534b4c52a5a3dcc717d781a52510d17c735511"),
    Linux(:i686, :glibc) => ("$bin_prefix/liblame.v3.100.0.i686-linux-gnu.tar.gz", "97e12758e2dfb72cf788a2f968ef82e7e1ed772564b0fa097087180bfeb3a370"),
    Linux(:i686, :musl) => ("$bin_prefix/liblame.v3.100.0.i686-linux-musl.tar.gz", "7a254261af6e6d3fe9ff3dcbad91a5d2aa22ab9b943e61171cc2d2291c616422"),
    Linux(:powerpc64le, :glibc) => ("$bin_prefix/liblame.v3.100.0.powerpc64le-linux-gnu.tar.gz", "b39355b30705e7332e3d4097ac1cd7a6c4d7d0297c589cdc31a804321fb5bee9"),
    Linux(:x86_64, :glibc) => ("$bin_prefix/liblame.v3.100.0.x86_64-linux-gnu.tar.gz", "879f3e79384d0a24408d4a49c46b1934be9d379eb21967bb2fb2786531681909"),
    Linux(:x86_64, :musl) => ("$bin_prefix/liblame.v3.100.0.x86_64-linux-musl.tar.gz", "7ab49f13b4606371f045efcb104fc439b6eb8d606ed8b0162a398138a780cba5"),
    FreeBSD(:x86_64) => ("$bin_prefix/liblame.v3.100.0.x86_64-unknown-freebsd11.1.tar.gz", "03cb5033e9984d7cd51f6169ad69b494248bdd9d94e40877a78d0e623aaea4ce"),
)

# Install unsatisfied or updated dependencies:
unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
if haskey(download_info, platform_key())
    url, tarball_hash = download_info[platform_key()]
    if unsatisfied || !isinstalled(url, tarball_hash; prefix=prefix)
        # Download and install binaries
        install(url, tarball_hash; prefix=prefix, force=true, verbose=verbose)
    end
elseif unsatisfied
    # If we don't have a BinaryProvider-compatible .tar.gz to download, complain.
    # Alternatively, you could attempt to install from a separate provider,
    # build from source or something even more ambitious here.
    error("Your platform $(triplet(platform_key())) is not supported by this package!")
end

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps_codec.jl"), products)


text:
function unrepeat{Array}(input::Array)
    holder=Set()
    repeats=Set()

    for i in input
        if i in holder
            push!(repeats, i)
        else
            push!(holder, i)
        end
    end
    println(collect(setdiff(holder, repeats))[1])
end

i1 = [2,3,4,2,3,5,4,6,4,6,9,10,9,8,7,8,10,7]
i2 = [2,'a','l',3,'l',4,'k',2,3,4,'a',6,'c',4,'m',6,'m','k',9,10,9,8,7,8,10,7]
unrepeat(i1)
unrepeat(i2)


text:
function assemble(op::Identity,
        testnfs::AbstractTimeBasisFunction,
        trialfns::AbstractTimeBasisFunction)

    tbf = convolve(testnfs, trialfns)
    has_zero_tail = all(tbf.polys[end].data .== 0)
    # @show has_zero_tail

    T = scalartype(tbf)
    if has_zero_tail
        z = zeros(T, numintervals(tbf)-1)
    else
        z = zeros(T, numfunctions(tbf))
    end

    Δt = timestep(tbf)
    #for i in eachindex(z)
    for i in 1:numintervals(tbf)-1
        p = tbf.polys[i]
        t = (i-1)*Δt
        z[i] = evaluate(p,t)
    end

    for i in numintervals(tbf):length(z)
        p = tbf.polys[end]
        t = (i-1)*Δt
        z[i] = evaluate(p,t)
    end

    return z
end

mutable struct TensorOperator <: Operator
    spatial_factor
    temporal_factor
end

⊗(A::AbstractOperator, B::AbstractOperator) = TensorOperator(A, B)
function scalartype(A::TensorOperator)
    promote_type(
        scalartype(A.spatial_factor),
        scalartype(A.temporal_factor))
end

# function assemble(operator::TensorOperator, testfns, trialfns)
#
#     test_spatial_basis  = spatialbasis(testfns)
#     trial_spatial_basis = spatialbasis(trialfns)
#
#     M = numfunctions(test_spatial_basis)
#     N = numfunctions(trial_spatial_basis)
#
#     tbf = convolve(
#         temporalbasis(testfns),
#         temporalbasis(trialfns))
#
#     kmax = numintervals(tbf) - 1
#     k0 = fill(1,    (M,N))
#     k1 = fill(kmax, (M,N))
#
#     T = promote_type(scalartype(operator), scalartype(testfns), scalartype(trialfns))
#     data = zeros(T, M, N, kmax)
#     Z = SparseND.Banded3D(k0, k1, data)
#
#     store(v, m, n, k) = (Z[m,n,k] += v)
#     assemble!(operator, testfns, trialfns, store)
#     return Z
#
# end

function allocatestorage(op::TensorOperator, test_functions, trial_functions,
    ::Type{Val{:bandedstorage}}, ::Type{LongDelays{:ignore}},)

    M = numfunctions(spatialbasis(test_functions))
    N = numfunctions(spatialbasis(trial_functions))

    time_basis_function = BEAST.convolve(
        temporalbasis(test_functions),
        temporalbasis(trial_functions))

    # has_zero_tail = all(time_basis_function.polys[end].data .== 0)
    # @show has_zero_tail

    # if has_zero_tail
        # Numintervals includes the semi-infinite interval stretching to +Inf
        # K = numintervals(time_basis_function)-1
    # else
    #     speedoflight = 1.0
    #     @warn "Assuming speed of light to be equal to 1!"
    #     Δt = timestep(time_basis_function)
    #     ct, hs = boundingbox(geometry(spatialbasis(trial_functions)).vertices)
    #     diam = 2 * sqrt(3) * hs
    #     K = ceil(Int, (numintervals(time_basis_function)-1) + diam/speedoflight/Δt)+1
    # end
    # @assert K > 0

    space_operator = op.spatial_factor
    A = assemble(space_operator, spatialbasis(test_functions), spatialbasis(trial_functions))

    K0 = Int.(A .!= 0)
    # K0 = ones(Int,M,N)
    bandwidth = numintervals(time_basis_function) - 1
    data = zeros(scalartype(op), bandwidth, M, N)
    maxk1 = bandwidth
    Z = SparseND.Banded3D(K0, data, maxk1)
    return Z, (v,m,n,k)->(Z[m,n,k] += v)
end

function allocatestorage(op::TensorOperator, test_functions, trial_functions)

    M = numfunctions(spatialbasis(test_functions))
    N = numfunctions(spatialbasis(trial_functions))

    time_basis_function = BEAST.convolve(
        temporalbasis(test_functions),
        temporalbasis(trial_functions))

    tbf = time_basis_function
    has_zero_tail = all(tbf.polys[end].data .== 0)
    @show has_zero_tail

    if has_zero_tail
        K = numintervals(time_basis_function)-1
    else
        speedoflight = 1.0
        @warn "Assuming speed of light to be equal to 1!"
        Δt = timestep(tbf)
        ct, hs = boundingbox(geometry(spatialbasis(trial_functions)).vertices)
        diam = 2 * sqrt(3) * hs
        K = ceil(Int, (numintervals(tbf)-1) + diam/speedoflight/Δt)+1
    end
    @assert K > 0

    Z = zeros(M, N, K)
    return MatrixConvolution(Z), (v,m,n,k)->(Z[m,n,k] += v)
end

function assemble!(operator::TensorOperator, testfns, trialfns, store)

    space_operator = operator.spatial_factor
    time_operator  = operator.temporal_factor

    space_testfns = spatialbasis(testfns)
    space_trialfns = spatialbasis(trialfns)

    time_testfns = temporalbasis(testfns)
    time_trialfns = temporalbasis(trialfns)

    zt = assemble(time_operator, time_testfns, time_trialfns)

    # Truncate to a reasonable size in case the tbf has a semi-infinite support
    tbf = convolve(time_testfns, time_trialfns)
    has_zero_tail = all(tbf.polys[end].data .== 0)
    if !has_zero_tail
        speedoflight = 1.0
        @warn "Assuming speed of light to be equal to 1!"
        Δt = timestep(tbf)
        ct, hs = boundingbox(geometry(space_trialfns).vertices)
        diam = 2 * sqrt(3) * hs
        kmax = ceil(Int, (numintervals(tbf)-1) + diam/speedoflight/Δt)+1
        zt = zt[1:kmax]
    end

    function store1(v,m,n)
        for (k,w) in enumerate(zt)
            store(w*v,m,n,k)
        end
    end
    assemble!(space_operator, space_testfns, space_trialfns, store1)

end

mutable struct TemporalDifferentiation <: Operator
    operator
end

derive(op::Operator) = TemporalDifferentiation(op)
scalartype(op::TemporalDifferentiation) = scalartype(op.operator)
Base.:*(a::Number, op::TemporalDifferentiation) = TemporalDifferentiation(a * op.operator)

function assemble!(operator::TemporalDifferentiation, testfns, trialfns, store)

    trial_time_fns  = temporalbasis(trialfns)
    trial_space_fns = spatialbasis(trialfns)

    trialfns = SpaceTimeBasis(
        trial_space_fns,
        derive(trial_time_fns)
    )

    assemble!(operator.operator, testfns, trialfns, store)

end

struct TemporalIntegration <: Operator
    operator::Operator
end

integrate(op::Operator) = TemporalIntegration(op)
scalartype(op::TemporalIntegration) = scalartype(op.operator)
Base.:*(a::Number, op::TemporalIntegration) = TemporalIntegration(a * op.operator)

function allocatestorage(op::TemporalIntegration, testfns, trialfns,
	::Type{Val{:bandedstorage}},
	::Type{LongDelays{:ignore}})

	trial_time_fns  = temporalbasis(trialfns)
	trial_space_fns = spatialbasis(trialfns)

	trialfns = SpaceTimeBasis(
		trial_space_fns,
		integrate(trial_time_fns)
	)

	Z, store = allocatestorage(op.operator, testfns, trialfns, Val{:bandedstorage}, LongDelays{:ignore})
	@show size(Z)
	return Z, store
end

function assemble!(operator::TemporalIntegration, testfns, trialfns, store)

    trial_time_fns  = temporalbasis(trialfns)
    trial_space_fns = spatialbasis(trialfns)

    trialfns = SpaceTimeBasis(
        trial_space_fns,
        integrate(trial_time_fns)
    )

    assemble!(operator.operator, testfns, trialfns, store)

end


text:
using SIUnits
using SIUnits.ShortUnits
using Base.Test

if VERSION <= v"0.3-"
    macro test_throws_compat(a,b)
        esc(:(@test_throws $b))
    end
else
    macro test_throws_compat(a,b)
        esc(:(@test_throws $a $b))
    end
end

# Basic arithmetic things
@test 1V + 2V == 3V
@test (1//2)V - 1V == (-1//2)V
@test_throws_compat ErrorException 1//2V - 1V
@test_throws_compat ErrorException 1V + 2s + 2kg

OneNewton = 1*(kg*m/s^2)
@test OneNewton*(1s)^2 == 1kg*m
@test OneNewton*s^2 == 1*kg*m

@test OneNewton/(kg*m) == 1Hertz^2

@test OneNewton^2 == 1kg^2*m^2/s^4

@test s == sqrt(1s^2)
@test sqrt(2)*m == sqrt(2m^2)

@test_throws_compat InexactError sqrt(1s)

@test 1/s == 1Hz

@test 1Hz*1s == 1
@test 1s/(1s) == 1

# Issue #2

immutable note{T<:Real}
    pitch::quantity(T,Hz)     #has dimensions of inverse time
    duration::quantity(T,s) #has dimensions of time
    sustained::Bool
end

note{Float64}(1.0Hz,1.0s,true)

@test_throws_compat ErrorException immutable foo{T}
    bar::quantity(T,2s)
end

# Non-SI Units
@test_approx_eq_eps 3mtorr 3*SIUnits.Milli*as(1torr,Pascal) eps()*Pa

# Unitful Powers
@test_throws_compat ErrorException 5.0^(3Pa)
@test_throws_compat ErrorException 5.0^(3torr)
@test 1^((5torr)/3Pa) == 1.0

# Interaction of si and non-si units
@test (1mm)^2*1torr == as(1torr,Pascal)*(1mm)^2
@test_throws_compat ErrorException as(1torr,s)

# Ranges (#4)
r1 = 1Hz:5Hz
@test length(r1) == 5

@test collect(1Hz:5Hz) == collect(1:5)Hz # Tests the iteration protocol

# Others

@test mod(2µm,4µm) == 2µm

# Issue #9
a = [1m 2N]
b = [1m 2N 3V]
@test a[1:2] == b[1:2]

# Issue #10
a = 1m
b = 2m
@test sqrt(a*b) == (a*b)^(1/2) == (a*b)^(1//2) == sqrt(2)*m

# Issue #11
au=[1m 2m 3m]
bu=[2N 3N 4N]
@test au*bu' == bu*au'
@test (au*bu')[1] == dot(vec(au),vec(bu))

# MathConsts
@test (3pi)*m == 3*(pi*m)
a = 1m
b = 2mm
@test (pi*a^2)/b == (pi/2)*(m^2/mm)

# Comparisons with numbers and different units
a = SIUnits.UnitQuantity{Float64}(3.0)
@test a < 4
@test 2.8 < a
@test_throws_compat MethodError 1m < 2kg


text:
export SkipSequential

"""
$(TYPEDEF)

Represents a sequential (feed-forward) neural net, with `layers` ordered from input
to output.

## Fields:
$(FIELDS)
"""
struct SkipSequential <: NeuralNet
    layers::Array{Layer,1}
    UUID::String
end

function Base.show(io::IO, p::SkipSequential)
    println(io, "skip sequential net $(p.UUID)")
    for (index, value) in enumerate(p.layers)
        println(io, "  ($index) $value")
    end
end

# TODO (vtjeng): Think about the types carefully.
function apply(p::SkipSequential, x::Array{<:JuMPReal})
    xs::Array{Array} = [x]
    for layer in p.layers
        if typeof(layer) <: SkipBlock
            output = layer(xs)
        else
            output = layer(last(xs))
        end
        push!(xs, output)
    end
    return last(xs)
end

(p::SkipSequential)(x::Array{<:JuMPReal}) = apply(p, x)


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule LibMPDec_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("LibMPDec")
JLLWrappers.@generate_main_file("LibMPDec", UUID("7106de7a-f406-5ef1-84f7-3345f7341bd2"))
end  # module LibMPDec_jll


text:
function JacobianFunction(arg0::MultivariateDifferentiableVectorFunction)
    return JacobianFunction((MultivariateDifferentiableVectorFunction,), arg0)
end

function value(obj::JacobianFunction, arg0::Vector{jdouble})
    return jcall(obj, "value", Vector{Vector{jdouble}}, (Vector{jdouble},), arg0)
end


text:
#=
Very heavily inspired by Calculus.jl, but with an emphasis on performance and DiffEq API convenience.
=#

#=
Compute the finite difference interval epsilon.
Reference: Numerical Recipes, chapter 5.7.
=#
@inline function compute_epsilon(::Val{:forward}, x::T, relstep::Real, absstep::Real, dir::Real) where T<:Number
    return max(relstep*abs(x), absstep)*dir
end

@inline function compute_epsilon(::Val{:central}, x::T, relstep::Real, absstep::Real, dir=nothing) where T<:Number
    return max(relstep*abs(x), absstep)
end

@inline function compute_epsilon(::Val{:hcentral}, x::T, relstep::Real, absstep::Real, dir=nothing) where T<:Number
    return max(relstep*abs(x), absstep)
end

@inline function compute_epsilon(::Val{:complex}, x::T, ::Union{Nothing,T}=nothing, ::Union{Nothing,T}=nothing, dir=nothing) where T<:Real
    return eps(T)
end

default_relstep(v::Type, T) = default_relstep(v(), T)
@inline function default_relstep(::Val{fdtype}, ::Type{T}) where {fdtype,T<:Number}
    if fdtype==:forward
        return sqrt(eps(real(T)))
    elseif fdtype==:central
        return cbrt(eps(real(T)))
    elseif fdtype==:hcentral
        eps(T)^(1/4)
    else
        return one(real(T))
    end
end

function fdtype_error(::Type{T}=Float64) where T
    if T<:Real
        error("Unrecognized fdtype: valid values are Val{:forward}, Val{:central} and Val{:complex}.")
    elseif T<:Complex
        error("Unrecognized fdtype: valid values are Val{:forward} or Val{:central}.")
    else
        error("Unrecognized returntype: should be a subtype of Real or Complex.")
    end
end


text:
using LittleScienceTools
using LittleScienceTools.Measuring
# using LittleScienceTools.Random
# using LittleScienceTools.Roots
# using LittleScienceTools.Ising
# using Compat
# using Erdos
using Test
using Statistics
using LinearAlgebra
using Random
using KahanSummation

@testset "Measuring" begin
    include("Measuring/observable.jl")
    include("Measuring/obstable.jl")
end

# @testset "Ising" begin
#     include("Ising/runtests.jl")
# end

# @testset "Roots" begin
#     include("Roots/Roots.jl")
# end

# @testset "Random" begin
# #     include("Random/parisi_rapuano.jl")
#     include("Random/discrete_distribution.jl")
#     include("Random/Random.jl")
# end

# @testset "utils" begin
#     include("utils.jl")
# end

max_stars_repo_name: JuliaTagBot/LittleScienceTools

max_issues_repo_name: JuliaTagBot/LittleScienceTools

max_forks_repo_name: JuliaTagBot/LittleScienceTools


text:
export ConjugateGradientMethod

"""
    ConjugateGradientMethod(; M=30)

# CG
This object represents an iterative Krylov method for solving a linear system.
The constructor parameter `M` is the number maximum number of iterations. This
uses the preconditioned Conjugate Gradient method of Barrett et al. (1994).

## References
    @book{barrett1994cg,
    title = {Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods},
    author = {
            Barrett, Richard and Berry, Michael and Chan, Tony F. and Demmel, James
            and Donato, June and Dongarra, Jack and Eijkhout, Victor and Pozo, Roldan
            and Romine, Charles and van der Vorst, Henk
        },
    pages={12-15},
    year = {1994},
    publisher = {SIAM},
    doi = {10.1137/1.9781611971538}
    }
"""
struct ConjugateGradientMethod <: AbstractKrylovMethod
    "Maximum number of CG iterations"
    M::Int
    function ConjugateGradientMethod(; M=30)
        return new(M)
    end
end

mutable struct CGCache{AT} <: AbstractLinearSolverCache
    r::AT
    z::AT
    p::AT
    Ap::AT
end

function cache(
    krylov_alg::ConjugateGradientMethod,
    Q::AT,
) where {AT}
    r = similar(Q)
    z = similar(Q)
    p = similar(Q)
    Ap = similar(Q)

    return CGCache{AT}(
        r,
        z,
        p,
        Ap,
    )
end

function LSinitialize!(
    ::ConjugateGradientMethod,
    solver::LinearSolver,
    Q,
    Qrhs,
    args...,
)
    linearoperator! = solver.linop!
    pc = solver.pc

    cache = solver.cache
    r = cache.r
    z = cache.z
    p = cache.p

    linearoperator!(r, Q, args...)
    r .= Qrhs .- r

    residual_norm = norm(r)
    threshold = solver.rtol * residual_norm

    converged = false
    if threshold < solver.atol
        converged = true
        return converged, threshold
    end

    isa(pc.pc_side, PCleft) || error("Only supports left preconditioning")
    PCapply!(pc, z, r, args...)
    p .= z

    converged, max(threshold, solver.atol)
end

function LSsolve!(
    krylov_alg::ConjugateGradientMethod,
    solver::LinearSolver,
    threshold,
    Q,
    Qrhs,
    args...,
)
    converged = false
    iter = 0
    residual_norm = typemax(eltype(Q))

    linearoperator! = solver.linop!
    pc = solver.pc
    cache = solver.cache
    r = cache.r
    z = cache.z
    p = cache.p
    Ap = cache.Ap

    ω0 = dot(r, z)
    while !converged && iter < krylov_alg.M
        linearoperator!(Ap, p, args...)
        α = ω0 / dot(p, Ap)
        Q .+= α .* p
        r .= r .- α .* Ap
        residual_norm = norm(r)
        if residual_norm < threshold
            converged = true
            break
        end
        # apply preconditioning
        PCapply!(pc, z, r, args...)
        ω1 = dot(r, z)
        β =  ω1 / ω0
        p .= z .+ β .* p
        ω0 = ω1
        iter += 1
    end

    (converged, iter, residual_norm)
end

max_stars_repo_name: CliMA/SystemSolvers

max_issues_repo_name: CliMA/SystemSolvers


text:
# TODO atexit(f) to store the job history.

SCHEDULER_TASK = @task scheduler()

function new_scheduler_task()
    global SCHEDULER_TASK
    SCHEDULER_TASK = @task scheduler()
end

"""
    scheduler_start()

Start the job scheduler.
"""
function scheduler_start(; verbose=true)
    global SCHEDULER_TASK

    if istaskfailed(SCHEDULER_TASK) || istaskdone(SCHEDULER_TASK)
        verbose && @warn "Scheduler was interrupted or done. Restart."
        new_scheduler_task()
        schedule(SCHEDULER_TASK)
    elseif istaskstarted(SCHEDULER_TASK) # if done, started is also true
        verbose && @warn "Scheduler is running."
    else
        verbose && @info "Scheduler starts."
        schedule(SCHEDULER_TASK)
    end
end

"""
    scheduler_stop()

Stop the job scheduler.
"""
function scheduler_stop(; verbose=true)
    global SCHEDULER_TASK
    if istaskfailed(SCHEDULER_TASK) || istaskdone(SCHEDULER_TASK)
        verbose && @warn "Scheduler is not running."
    elseif istaskstarted(SCHEDULER_TASK) # if done, started is also true
        wait_for_lock()
            schedule(SCHEDULER_TASK, InterruptException; error=true)
        release_lock()
        verbose && @info "Scheduler stops."
    else
        verbose && @warn "Scheduler is not running."
    end
end

"""
    scheduler_status() :: Symbol

Print the settings and status of job scheduler. Return `:not_running` or `:running`.
"""
function scheduler_status(; verbose=true)
    global SCHEDULER_TASK
    global SCHEDULER_MAX_CPU
    global SCHEDULER_MAX_MEM
    global SCHEDULER_UPDATE_SECOND
    global JOB_QUEUE_MAX_LENGTH
    global SCHEDULER_TASK
    if istaskfailed(SCHEDULER_TASK) || istaskdone(SCHEDULER_TASK)
        verbose && @info "Scheduler is not running." SCHEDULER_MAX_CPU SCHEDULER_MAX_MEM SCHEDULER_UPDATE_SECOND JOB_QUEUE_MAX_LENGTH SCHEDULER_TASK
        :not_running
    elseif istaskstarted(SCHEDULER_TASK)
        verbose && @info "Scheduler is running." SCHEDULER_MAX_CPU SCHEDULER_MAX_MEM SCHEDULER_UPDATE_SECOND JOB_QUEUE_MAX_LENGTH SCHEDULER_TASK
        :running
    else
        verbose && @info "Scheduler is not running." SCHEDULER_MAX_CPU SCHEDULER_MAX_MEM SCHEDULER_UPDATE_SECOND JOB_QUEUE_MAX_LENGTH SCHEDULER_TASK
        :not_running
    end
end

"""
    set_scheduler_update_second(s::Float64 = 0.6)

Set the update interval of scheduler.
"""
function set_scheduler_update_second(s::Float64 = 0.6)
    s <= 0.001 && error("schedular update interval cannot be less than 0.001.")
    global SCHEDULER_UPDATE_SECOND = s
end
set_scheduler_update_second(s) = set_scheduler_update_second(convert(Float64, s))

"""
    set_scheduler_max_cpu(ncpu::Int = nthreads()>1 ? nthreads()-1 : Sys.CPU_THREADS)
    set_scheduler_max_cpu(percent::Float64)

Set the maximum CPU (thread) the scheduler can use. If starting Julia with multi-threads, the maximum CPU is `nthreads() - 1`.

# Example
    set_scheduler_max_cpu()     # use all available CPUs
    set_scheduler_max_cpu(4)    # use 4 CPUs
    set_scheduler_max_cpu(0.5)  # use 50% of CPUs
"""
function set_scheduler_max_cpu(ncpu::Int = nthreads() > 1 ? nthreads()-1 : Sys.CPU_THREADS)
    ncpu < 1 && error("number of CPU cannot be less than 1.")
    if ncpu > Sys.CPU_THREADS
        @warn "Assigning number of CPU > total CPU."
    end
    # nthreads == 1 will not triger schedule at different jobs.
    if nthreads() > 1
        if ncpu > nthreads()-1
            error("Assigning number of CPU > total threads available - 1. Thread one is reserved for schedulers. Try to start Julia with sufficient threads. Help: https://docs.julialang.org/en/v1/manual/multi-threading/#Starting-Julia-with-multiple-threads")
        elseif ncpu < 2
            error("Assigning number of CPU < 2 is not allowed in multi-threaded Julia. Thread one is reserved for schedulers.")
        end
    end
    global SCHEDULER_MAX_CPU = ncpu
end
function set_scheduler_max_cpu(percent::Float64)
    if 0.0 < percent <= 1.0
        ncpu = round(Int, Sys.CPU_THREADS * percent)
        if ncpu > nthreads() - 1
            @warn "Assigning number of CPU > total threads available - 1 is not allowed. Set to total threads available - 1. Thread one is reserved for schedulers. To use more threads, try to start Julia with sufficient threads. Help: https://docs.julialang.org/en/v1/manual/multi-threading/#Starting-Julia-with-multiple-threads"
            ncpu = nthreads() - 1
        end
        set_scheduler_max_cpu(ncpu)
    else
        @error "Percent::Float64 should be between 0 and 1. Are you looking for set_scheduler_max_cpu(ncpu::Int) ?"
    end
end

"""
    set_scheduler_max_mem(mem::Int = round(Int, Sys.total_memory() * 0.8))
    set_scheduler_max_mem(percent::Float64)

Set the maximum RAM the scheduler can use.

# Example
    set_scheduler_max_mem()             # use 80% of total memory

    set_scheduler_max_mem(4GB)          # use 4GB memory
    set_scheduler_max_mem(4096MB)
    set_scheduler_max_mem(4194304KB)
    set_scheduler_max_mem(4294967296B)

    set_scheduler_max_mem(0.5)          # use 50% of total memory

"""
function set_scheduler_max_mem(mem::Int = round(Int, Sys.total_memory() * 0.8))
    mem < 1 && error("number of memory cannot be less than 1.")
    if mem > Sys.total_memory() * 0.9
        @warn "Assigning memory > 90% of total memory."
    end
    global SCHEDULER_MAX_MEM = mem
end
function set_scheduler_max_mem(percent::Float64)
    if 0.0 < percent < 1.0
        set_scheduler_max_mem(round(Int, Sys.total_memory() * percent))
    else
        @error "Percent::Float64 should be between 0 and 1. Are you looking for set_scheduler_max_mem(mem::Int) ?"
    end
end

"""
    set_scheduler_max_job(n_finished_jobs::Int = 10000)

Set the number of finished jobs.
"""
function set_scheduler_max_job(n_finished_jobs::Int = 10000)
    if n_finished_jobs < 10
        @error "Cannot set number of finished jobs < 10"
    else
        global JOB_QUEUE_MAX_LENGTH = n_finished_jobs
    end
end


text:
using PiecewiseDeterministicMarkovProcesses, Random
const PDMP = PiecewiseDeterministicMarkovProcesses
jp = PDMP.Jump(PDMP.Delta_dummy)
jp = PDMP.Jump(rand(Int64,2,2))

PDMP.finalize_dummy(1, 1, 1, 1, 1)
PDMP.PDMPResult(rand(2), rand(2), rand(2))
PDMP.PDMPResult(rand(2), rand(2), rand(2), rand(2), (false, false))

xc0 = rand(2)
xdot0 = rand(2,2)
PDMP.Phi_dummy(xdot0, xc0, 1, 1, 1)


text:
function get_num_active_bots(u::String="*"; base_path::String=Pinguu.base_dir)
    config_paths = glob(joinpath(base_path, "users/" * u * "/*/config.json"));
    num = 0
    for path in config_paths
        tmp = JSON.parsefile(path)
        num = tmp["status"] != "dead" ? num + 1 : num
    end
    config_paths = glob(joinpath(base_path, "grid_bots/user/" * u * "/*/config.json"));
    for path in config_paths
        tmp = JSON.parsefile(path)
        if tmp["testing"] == "0"
            num = tmp["status"] != "finished" ? num + 1 : num
        end
    end
    return num
end
export get_num_active_bots


text:
@inline add_fpfp_dd(x::T, y::T) where {T<:AbstractFloat} = add_2(x, y)
@inline sub_fpfp_dd(x::T, y::T) where {T<:AbstractFloat} = sub_2(x, y)
@inline mul_fpfp_dd(x::T, y::T) where {T<:AbstractFloat} = mul_2(x, y)
@inline dvi_fpfp_dd(x::T, y::T) where {T<:AbstractFloat} = dvi_2(x, y)
# seems to help
@inline add_fpfp_dd(x::T, y::T) where {T<:IEEEFloat} = add_2(x, y)
@inline sub_fpfp_dd(x::T, y::T) where {T<:IEEEFloat} = sub_2(x, y)
@inline mul_fpfp_dd(x::T, y::T) where {T<:IEEEFloat} = mul_2(x, y)
@inline dvi_fpfp_dd(x::T, y::T) where {T<:IEEEFloat} = dvi_2(x, y)


text:
import IterativeSolvers
import IterativeSolvers: lsqr

function pinv_iter(A::AbstractMatrix{T}, method::Symbol=:lsqr) where {T<:Number}
    iA = zeros(Float64, size(A'))
    try
        iA = pinv(A)
    catch
        for c = 1:size(iA, 2)
            b = zeros(Float64, size(A,1))
            b[c] = 1.
            iA[:,c] = IterativeSolvers.lsqr(A, b)
        end
    end
    return iA
end

max_stars_repo_name: yanncalec/FracFin

max_issues_repo_name: yanncalec/FracFin

max_forks_repo_name: yanncalec/FracFin


text:
using LanguageServer, LanguageServer.SymbolServer

unquote(s) = length(s) >= 2 && (s[1] == '"' && s[end] == '"') ? s[2:end-1] : s
envpath = unquote(Base.ARGS[1])

runserver(stdin, stdout, envpath)

max_stars_repo_name: machakann/vim-lsp-julia

max_issues_repo_name: machakann/vim-lsp-julia

max_forks_repo_name: machakann/vim-lsp-julia


text:
function AdamsMoultonIntegratorBuilder(arg0::jint, arg1::jdouble, arg2::jdouble, arg3::jdouble)
    return AdamsMoultonIntegratorBuilder((jint, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3)
end

function build_integrator(obj::AdamsMoultonIntegratorBuilder, arg0::Orbit, arg1::OrbitType)
    return jcall(obj, "buildIntegrator", AbstractIntegrator, (Orbit, OrbitType), arg0, arg1)
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Test, Random
using CRC32c

function test_crc32c(crc32c)
    # CRC32c checksum (test data generated from @andrewcooke's CRC.jl package)
    for (n,crc) in [(0,0x00000000),(1,0xa016d052),(2,0x03f89f52),(3,0xf130f21e),(4,0x29308cf4),(5,0x53518fab),(6,0x4f4dfbab),(7,0xbd3a64dc),(8,0x46891f81),(9,0x5a14b9f9),(10,0xb219db69),(11,0xd232a91f),(12,0x51a15563),(13,0x9f92de41),(14,0x4d8ae017),(15,0xc8b74611),(16,0xa0de6714),(17,0x672c992a),(18,0xe8206eb6),(19,0xc52fd285),(20,0x327b0397),(21,0x318263dd),(22,0x08485ccd),(23,0xea44d29e),(24,0xf6c0cb13),(25,0x3969bba2),(26,0x6a8810ec),(27,0x75b3d0df),(28,0x82d535b1),(29,0xbdf7fc12),(30,0x1f836b7d),(31,0xd29f33af),(32,0x8e4acb3e),(33,0x1cbee2d1),(34,0xb25f7132),(35,0xb0fa484c),(36,0xb9d262b4),(37,0x3207fe27),(38,0xa024d7ac),(39,0x49a2e7c5),(40,0x0e2c157f),(41,0x25f7427f),(42,0x368c6adc),(43,0x75efd4a5),(44,0xa84c5c31),(45,0x0fc817b2),(46,0x8d99a881),(47,0x5cc3c078),(48,0x9983d5e2),(49,0x9267c2db),(50,0xc96d4745),(51,0x058d8df3),(52,0x453f9cf3),(53,0xb714ade1),(54,0x55d3c2bc),(55,0x495710d0),(56,0x3bddf494),(57,0x4f2577d0),(58,0xdae0f604),(59,0x3c57c632),(60,0xfe39bbb0),(61,0x6f5d1d41),(62,0x7d996665),(63,0x68c738dc),(64,0x8dfea7ae)]
        @test crc32c(UInt8[1:n;]) == crc == crc32c(String(UInt8[1:n;]))
    end

    # test that crc parameter is equivalent to checksum of concatenated data,
    # and test crc of subarrays:
    a = UInt8[1:255;]
    crc_256 = crc32c(a)
    @views for n = 1:255
        @test crc32c(a[n+1:end], crc32c(a[1:n])) == crc_256
    end
    @test crc32c(IOBuffer(a)) == crc_256
    let buf = IOBuffer()
        write(buf, a[1:3])
        @test crc32c(seekstart(buf)) == crc32c(a[1:3])
        @test crc32c(buf) == 0x00000000
        @test crc32c(seek(buf, 1)) == crc32c(a[2:3])
        @test crc32c(seek(buf, 0), 2) == crc32c(a[1:2])
        @test crc32c(buf) == crc32c(a[3:3])
    end

    let f = tempname()
        try
            write(f, a)
            @test open(crc32c, f) == crc_256
            open(f, "r") do io
                @test crc32c(io, 16) == crc32c(a[1:16])
                @test crc32c(io, 16) == crc32c(a[17:32])
                @test crc32c(io) == crc32c(a[33:end])
                @test crc32c(io, 1000) == 0x00000000
            end
            a = rand(UInt8, 30000)
            write(f, a)
            @test open(crc32c, f) == crc32c(a) == open(io -> crc32c(io, 10^6), f)
        finally
            rm(f, force=true)
        end
    end
end
unsafe_crc32c_sw(a, n, crc) =
    ccall(:jl_crc32c_sw, UInt32, (UInt32, Ptr{UInt8}, Csize_t), crc, a, n)
crc32c_sw(a::Union{Array{UInt8},Base.FastContiguousSubArray{UInt8,N,<:Array{UInt8}} where N},
          crc::UInt32=0x00000000) = unsafe_crc32c_sw(a, length(a), crc)
crc32c_sw(s::String, crc::UInt32=0x00000000) = unsafe_crc32c_sw(s, sizeof(s), crc)
function crc32c_sw(io::IO, nb::Integer, crc::UInt32=0x00000000)
    nb < 0 && throw(ArgumentError("number of bytes to checksum must be ≥ 0"))
    buf = Vector{UInt8}(uninitialized, min(nb, 24576))
    while !eof(io) && nb > 24576
        n = readbytes!(io, buf)
        crc = unsafe_crc32c_sw(buf, n, crc)
        nb -= n
    end
    return unsafe_crc32c_sw(buf, readbytes!(io, buf, min(nb, length(buf))), crc)
end
crc32c_sw(io::IO, crc::UInt32=0x00000000) = crc32c_sw(io, typemax(Int64), crc)
test_crc32c(crc32c)
test_crc32c(crc32c_sw)

max_stars_repo_name: blegat/julia

max_issues_repo_name: blegat/julia

max_forks_repo_name: blegat/julia


text:
# The heat capacity of water

using Dierckx  # for interpolation

export HeatCapacity, GridHeatCapacity

# Types for different heat capacity behaviours

"An isobaric heat capacity (cₚ)"
abstract type HeatCapacity end
"A heat capacity that's not constant"
abstract type VaryingHeatCapacity <: HeatCapacity end
"A heat capacity that depends varies according to some functional form"
abstract type FunctionalHeatCapacity <: VaryingHeatCapacity end

"A constant heat capacity (no T/P variation)"
struct ConstantHeatCapacity{T<:Real} <: HeatCapacity
    value::T
end

"A heat capacity which is a function of temperature"
struct TFuncHeatCapacity{F} <: FunctionalHeatCapacity
    func::F
end

"A heat capacity which is a function of pressure and temperature"
struct PTFuncHeatCapacity{F} <: FunctionalHeatCapacity
    func::F
end

"A heat capacity interpolated from a log-linear grid"
struct GridHeatCapacity <: HeatCapacity
    P::Vector{Float64}
    T::Vector{Float64}
    spline::Spline2D

    function GridHeatCapacity(P, T, c_p)
        new(P, T, Spline2D(P, T, c_p, kx=1, ky=1))
    end
end

# Making heat capacities from files

"Generate a heat capacity from a file"
function GridHeatCapacity(filename, scale=1000)
    data = readdlm(filename, Float64)
    T = vec(data[2:end, 1])
    P = vec(data[1, 2:end])
    c_p = Matrix{Float64}(data[2:end, 2:end])' * scale  # note the transpose
    # default scaling factor is 1000, corresponding to the values in the file
    # being provided in g/cm^3/K (cgs) and converted to kg/m^3/K (SI base)

    GridHeatCapacity(P, T, c_p)
end

# Evaluating heat capacities

(cp::GridHeatCapacity)(P, T) = evaluate(cp.spline, P, T)
(cp::ConstantHeatCapacity)(P) = cp.value
(cp::ConstantHeatCapacity)(P, T) = cp.value
(cp::TFuncHeatCapacity)(T) = cp.func(T)
(cp::TFuncHeatCapacity)(P, T) = cp.func(T)
(cp::PTFuncHeatCapacity)(P, T) = cp.func(P, T)

# Make and save heat capacity

"Save the water heat capacity to a JLD file"
function save_heat_capacity!()
    jldopen("$(config.datadir)/heatcapacity.jld", "w") do file
        # make the heat capacity
        c_p = GridHeatCapacity("$(config.rawdata)/heatcap-h2o.dat")
        write(file, "heatcap_h2o", c_p)
    end
end


text:
# gaussianprocess.jl

"""
    abstract type AbstractGaussianProcess end

Subtypes: GaussianProcess and FaberSchauderExpansionWithGaussianCoefficients.
Supertype: Any

Both subtypes represent a Gaussian process as a series expansion
in a basis of functions, where the vector of functions is multivariate normally
distributed.

See also: GaussianProcess and FaberSchauderExpansionWithGaussianCoefficients.
"""
abstract type AbstractGaussianProcess end

"""
  struct GaussianProcess{S, T} <: AbstractGaussianProcess where {
      S<:AbstractVector{U} where U<:Function, T<:Union{AbstractMvNormal,GaussianVector}}
    basis::S
    distribution::T
  end

Implements a Gaussian process defined as an inner product of a Gaussian vector
with distribution 'distribution' of type GaussianVector or a subtype of AbstractMvNormal
(from the Distributions.jl package) and a vector of functions,
the basis of the function space. For the use of a Gaussian process expanded in
Faber-Schauder functions up to level j, we recommend using
FaberSchauderExpansionWithGaussianCoefficients, as this makes efficiently use of
the sparsity structure of Faber-Schauder functions.

See also: FaberSchauderExpansionWithGaussianCoefficients, AbstractGaussianProcess

# Example

```julia
using Distributions
distribution = MvNormal([k^(-1.0) for k in 1:10])
Π = GaussianProcess([fourier(k) for k in 1:10], distribution)
f = rand(Π)
```
"""
struct GaussianProcess{S, T} <: AbstractGaussianProcess where {
    S<:AbstractVector{U} where U<:Function, T<:Union{AbstractMvNormal,GaussianVector}}
  basis::S
  distribution::T
  function GaussianProcess(basis::S, distribution::T) where {
        S<:AbstractArray{<:Function}, T<:Union{AbstractMvNormal,GaussianVector}}
      length(basis)==length(distribution) || throw(DimensionMismatch(
        "The basis and distribution are not of equal length."))
      new{S, T}(basis, distribution)
  end
end

# Tests need to be written.
"""
    calculateboundssupport(higestlevel::Int)

Internal function, not exported.

Returns a tuple of two Float64 vectors of length 2^(higestlevel+1)
where the i-th element of the first vector is the left bound of the
support of the i-th Faber-Schauder basis function and the i-th
element of the second vector is the right bound of the support of
i-th Faber-Schauder basis function.
"""
function calculateboundssupport(higestlevel::Int)
    higestlevel ≥ 0 || throw(AssertionError("higestlevel should be positive."))
    lengthvector = 2^(higestlevel+1)
    leftboundsupport = Vector{Float64}(undef, lengthvector)
    rightboundsupport = Vector{Float64}(undef, lengthvector)
    leftboundsupport[1] = leftboundsupport[2] = 0.0
    rightboundsupport[1] = rightboundsupport[2] = 1.0
    for j in 1:higestlevel
        for k in 1:2^j
            i = 2^j+k
            leftboundsupport[i] = ldexp(float(k-1), -j)
            rightboundsupport[i] = ldexp(float(k), -j)
        end
    end
    return (leftboundsupport, rightboundsupport)
end

"""
    createFaberSchauderBasisUpToLevelHigestLevel(higestlevel::Int)

Internal function, not exported.

Creates Faber-Schauder basis up to level higestlevel. Returns a
Vector{Function} type of length 2^(higestlevel+1). The first element
of the vector is faberschauderone, and the 2^j+k element is
faberschauder(j,k), with 0≤j≤higestlevel and 1≤k≤2^j.
"""
function createFaberSchauderBasisUpToLevelHigestLevel(higestlevel::Int)
    higestlevel≥0 || throw(AssertionError("The levels are non-negative."))
    # There are 2^(higestlevel+1) basis functions.
    basis = Vector{Function}(undef, 2^(higestlevel+1))
    basis[1] = faberschauderone
    for j in 0:higestlevel
        for k in 1:2^j
            basis[2^j+k] = faberschauder(j,k)
        end
    end
    return basis
end

"""
    struct FaberSchauderExpansionWithGaussianCoefficients{T} <:
            AbstractGaussianProcess where T<:Union{AbstractMvNormal,GaussianVector}
        higestlevel::Int64
        basis::Vector{Function}
        leftboundssupport::Vector{Float64}
        rightboundssupport::Vector{Float64}
        distribution::T
    end

Implements a Gaussian Process with Faber-Schauder functions and optimally exploits
the sparsity structure of Faber-Schauder functions.

Constructors:
    FaberSchauderExpansionWithGaussianCoefficients(higestlevel, distribution)
    FaberSchauderExpansionWithGaussianCoefficients(standarddeviationsperlevel::Vector{Float64})

The user is not allowed to set a basis, which the constructor calculates from the input and
is by its very definition a Faber-Schauder basis. The length of distribution
should be equal to 2^(higestlevel+1). The second constructor defines a Gaussian
process with Faber-Schauder basis with independent coefficients with
length(standarddeviationsperlevel)-1 levels (there is also a level zero), so
2^(length(standarddeviationsperlevel)) number of basis functions, where the standard deviation
of the coefficients belonging to level k is standarddeviationsperlevel[k+1].

# Example

```julia
α = 0.5
Π = FaberSchauderExpansionWithGaussianCoefficients([2^(-j*α) for j in 0:5])
f = rand(Π)
```
"""
struct FaberSchauderExpansionWithGaussianCoefficients{T} <:
         AbstractGaussianProcess where T<:Union{AbstractMvNormal, GaussianVector}
     higestlevel::Int64
     basis::Vector{Function}
     leftboundssupport::Vector{Float64}
     rightboundssupport::Vector{Float64}
     distribution::T
     function FaberSchauderExpansionWithGaussianCoefficients(higestlevel::Int64,
             distribution::T) where {T<:Union{AbstractMvNormal,GaussianVector}}
         higestlevel ≥ 0 || throw(AssertionError("higestlevel should be non-negative."))
         length(distribution) == 2^(higestlevel+1) || throw(AssertionError("The length of the
         distribution is not equal to 2^(higestlevel+1)."))
         basis = createFaberSchauderBasisUpToLevelHigestLevel(higestlevel)
         leftboundssupport, rightboundssupport = calculateboundssupport(
             higestlevel)
         new{T}(higestlevel, basis,
             leftboundssupport, rightboundssupport, distribution)
     end
end

"""
    createvectorofstandarddeviationsfromstandarddeviationsperlevel(standarddeviationsperlevel::AbstractArray{Float64})

Internal function, not exported!

Calculated a vector of standard deviations for all basis functions, where basis function i gets standard deviation
standarddeviationsperlevel[ℓ+1] if ℓ is the level of i (we start at level 0).
"""
function createvectorofstandarddeviationsfromstandarddeviationsperlevel(standarddeviationsperlevel::AbstractArray{Float64})
    lenghtstandarddeviationsperlevel = length(standarddeviationsperlevel)
    lenghtstandarddeviationsperlevel > 0 || throw(AssertionError("standarddeviationsperlevel is of zero length."))
    # The distribution has length 2^lenghtstandarddeviationsperlevel=2^(higestlevel+1).
    vectorofstandarddeviations = Vector{Float64}(undef, 2^lenghtstandarddeviationsperlevel)
    # There are two functions of level zero.
    vectorofstandarddeviations[1:2] .= standarddeviationsperlevel[1]
    # and 2^k of level k, k=1,2,...
    for k in 1:lenghtstandarddeviationsperlevel-1
        vectorofstandarddeviations[2^k+1:2^(k+1)] .= standarddeviationsperlevel[k+1]
    end
    return vectorofstandarddeviations
end

# Constructor
function FaberSchauderExpansionWithGaussianCoefficients(
        standarddeviationsperlevel::AbstractArray{Float64})
    lenghtstandarddeviationsperlevel = length(standarddeviationsperlevel)
    lenghtstandarddeviationsperlevel == 0 && throw(AssertionError("standarddeviationsperlevel is of
        zero length."))
    # We start with level zero.
    higestlevel = lenghtstandarddeviationsperlevel - 1
    vectorofstandarddeviations = createvectorofstandarddeviationsfromstandarddeviationsperlevel(standarddeviationsperlevel)
    distribution = MvNormal(vectorofstandarddeviations)
    return FaberSchauderExpansionWithGaussianCoefficients(higestlevel,
        distribution)
end

#Extends Base.length to objects of a subtypes of AbstractGaussianProcess.
"""
    length(Π::AbstractGaussianProcess)
    length(Π::FaberSchauderExpansionWithGaussianCoefficients)

length returns the number of basis functions, which is equal to the length of the distribution of the
coefficients, which is also equal to the number of coefficients.
For a FaberSchauderExpansionWithGaussianCoefficients object, the length is equal to 2^(higestlevel+1).

# Example
```julia
Π = GaussianProcess([sin, cos], MvNormal([1.,1.]))
length(Π)
#
Π = FaberSchauderExpansionWithGaussianCoefficients([2.0^(-0.5*j) for j in 0:3])
length(Π) # == 2^4 == 16.
```
"""
length(Π::AbstractGaussianProcess)=length(Π.basis)
length(Π::FaberSchauderExpansionWithGaussianCoefficients) = 2^(Π.higestlevel+1)

"""
    sumoffunctions(vectoroffunctions::AbstractArray{<:Function},
        vectorofscalars::AbstractArray{<:Number})

Internal function, not exported.

sumoffunctions calculates the 'inner product' of the function vector and the scalar vector. In
other words, the sum of the functions weighted by vectorofscalars. It returns a
function.
"""
function sumoffunctions(vectoroffunctions::AbstractArray{<:Function},
    vectorofscalars::AbstractArray{<:Number})
  n = length(vectoroffunctions)
  n == length(vectorofscalars) || throw(AssertionError("The vector of functions and the vector
    of scalars should be of equal length"))
  return x -> sum(vectorofscalars[k]*vectoroffunctions[k](x) for k in 1:n)
end

"""
    rand(Π::AbstractGaussianProcess)

rand returns a random function, where the coefficients have distribution Π.distribution
and the basis functions are defined in Π.basis.

# Example

```julia
using Distributions
distribution = MvNormal([k^(-1.0) for k in 1:100])
Π = GaussianProcess([fourier(k) for  k in 1:100], distribution)
f = rand(Π)
```
"""
function rand(Π::AbstractGaussianProcess)
  Z = rand(Π.distribution)
  return sumoffunctions(Π.basis, Z)
end

# Extends Base.mean to GaussianProcess objects.
"""
    mean(Π::AbstractGaussianProcess)

mean calculates the mean of the Gaussian process. It returns a function.

#Examples

```julia
using Distributions
distribution = MvNormal([k^(-1.0) for k in 1:100])
Π = GaussianProcess([fourier(k) for  k in 1:100], distribution)
mean(Π)
```
"""
mean(Π::AbstractGaussianProcess) = sumoffunctions(Π.basis, mean(Π.distribution))


text:
module OptiMimi

using NLopt
using ForwardDiff
using MathProgBase
using Compat

import Mimi: Model, CertainScalarParameter, CertainArrayParameter, addparameter

export problem, solution, unaryobjective, objevals, setparameters, nameindexes

include("registerdiff.jl")
include("matrixconstraints.jl")
include("linproghouse.jl")
include("makeroom.jl")

allverbose = false
objevals = 0

type OptimizationProblem
    model::Model
    components::Vector{Symbol}
    names::Vector{Symbol}
    opt::Opt
    constraints::Vector{Function}
end

type LinprogOptimizationProblem{T}
    model::Model
    components::Vector{Symbol}
    names::Vector{Symbol}
    objective::Function
    objectiveconstraints::Vector{Function}
    matrixconstraints::Vector{MatrixConstraintSet}
    exlowers::Vector{T}
    exuppers::Vector{T}
end

"""Returns (ii, len, isscalar) with the index of each symbol and its length."""
function nameindexes(model::Model, names::Vector{Symbol})
    ii = 1
    for name in names
        if isa(model.parameters[name], CertainScalarParameter)
            produce((ii, 1, true))
        elseif isa(model.parameters[name], CertainArrayParameter)
            produce((ii, length(model.parameters[name].values), false))
        else
            error("Unknown parameter type for " + string(name))
        end
        ii += 1
    end
end

"""Set parameters in a model."""
function setparameters(model::Model, components::Vector{Symbol}, names::Vector{Symbol}, xx::Vector)
    startindex = 1
    for (ii, len, isscalar) in @task nameindexes(model, names)
        if isscalar
            model.components[components[ii]].Parameters.(names[ii]) = xx[startindex]
        else
            shape = size(model.components[components[ii]].Parameters.(names[ii]))
            model.components[components[ii]].Parameters.(names[ii]) = reshape(collect(Number, xx[startindex:(startindex+len - 1)]), shape)
        end
        startindex += len
    end
end

"""Generate the form of objective function used by the optimization, taking parameters rather than a model."""
function unaryobjective(model::Model, components::Vector{Symbol}, names::Vector{Symbol}, objective::Function)
    function my_objective(xx::Vector)
        if allverbose
            println(xx)
        end

        global objevals
        objevals += 1

        setparameters(model, components, names, xx)
        run(model)
        objective(model)
    end

    my_objective
end

"""Create an NLopt-style objective function which does not use its grad argument."""
function gradfreeobjective(model::Model, components::Vector{Symbol}, names::Vector{Symbol}, objective::Function)
    myunaryobjective = unaryobjective(model, components, names, objective)
    function myobjective(xx::Vector, grad::Vector)
        myunaryobjective(xx)
    end

    myobjective
end

"""Create an NLopt-style objective function which computes an autodiff gradient."""
function autodiffobjective(model::Model, components::Vector{Symbol}, names::Vector{Symbol}, objective::Function)
    myunaryobjective = unaryobjective(model, components, names, objective)
    if VERSION < v"0.4.0-dev"
        # Slower: doesn't use cache
        function myobjective(xx::Vector, gradout::Vector)
            gradual = myunaryobjective(GraDual(xx))
            copy!(gradout, grad(gradual))
            value(gradual)
        end
    else
        function myobjective(xx::Vector, grad::Vector)
            out = GradientResult(xx)
            ForwardDiff.gradient!(out, myunaryobjective, xx)
            if any(isnan(ForwardDiff.gradient(out)))
                error("objective gradient is NaN")
            end
            copy!(grad, ForwardDiff.gradient(out))
            ForwardDiff.value(out)
        end
    end

    myobjective
end

"""Create a 0 point."""
function make0(model::Model, names::Vector{Symbol})
    initial = Float64[]
    for (ii, len, isscalar) in @task nameindexes(model, names)
        append!(initial, [0. for jj in 1:len])
    end

    initial
end

"""Setup an optimization problem."""
function problem{T<:Real}(model::Model, components::Vector{Symbol}, names::Vector{Symbol}, lowers::Vector{T}, uppers::Vector{T}, objective::Function; constraints::Vector{Function}=Function[], algorithm::Symbol=:LN_COBYLA_OR_LD_MMA)
    my_lowers = T[]
    my_uppers = T[]

    ## Replace with eachname
    totalvars = 0
    for (ii, len, isscalar) in @task nameindexes(model, names)
        append!(my_lowers, [lowers[ii] for jj in 1:len])
        append!(my_uppers, [uppers[ii] for jj in 1:len])
        totalvars += len
    end

    if algorithm == :GUROBI_LINPROG
        # Make no changes to objective!
    elseif model.numberType == Number
        if algorithm == :LN_COBYLA_OR_LD_MMA
            algorithm = :LD_MMA
        end
        if string(algorithm)[2] == 'N'
            warn("Model is autodifferentiable, but optimizing using a derivative-free algorithm.")
            myobjective = gradfreeobjective(model, components, names, objective)
        else
            println("Using AutoDiff objective.")
            myobjective = autodiffobjective(model, components, names, objective)
        end
    else
        if algorithm == :LN_COBYLA_OR_LD_MMA
            algorithm = :LN_COBYLA
        elseif string(algorithm)[2] == 'D'
            warn("Model is non-differentiable, but requested a gradient algorithm; instead using LN_COBYLA.")
            algorithm = :LN_COBYLA
        end

        myobjective = gradfreeobjective(model, components, names, objective)
    end

    if algorithm == :GUROBI_LINPROG
        LinprogOptimizationProblem(model, components, names, objective, constraints, MatrixConstraintSet[], my_lowers, my_uppers)
    else
        opt = Opt(algorithm, totalvars)
        lower_bounds!(opt, my_lowers)
        upper_bounds!(opt, my_uppers)
        xtol_rel!(opt, minimum(1e-6 * (uppers - lowers)))

        max_objective!(opt, myobjective)

        for constraint in constraints
            let this_constraint = constraint
                function my_constraint(xx::Vector, grad::Vector)
                    setparameters(model, components, names, xx)
                    this_constraint(model)
                end

                inequality_constraint!(opt, my_constraint)
            end
        end

        OptimizationProblem(model, components, names, opt, constraints)
    end
end

"""Solve an optimization problem."""
function solution(optprob::OptimizationProblem, generator::Function; maxiter=Inf, verbose=false)
    global allverbose
    allverbose = verbose

    if verbose
        println("Selecting an initial point.")
    end

    attempts = 0
    initial = []
    valid = false
    while attempts < maxiter
        initial = generator()

        setparameters(optprob.model, optprob.components, optprob.names, initial)

        valid = true
        for constraint in optprob.constraints
            if constraint(optprob.model) >= 0
                valid = false
                break
            end
        end

        if valid
            break
        end

        attempts += 1
        if attempts % 1000 == 0
            println("Could not find initial point after $attempts attempts.")
        end
    end

    if !valid
        throw(DomainError("Could not find a valid initial value."))
    end

    if verbose
        println("Optimizing...")
    end
    (minf,minx,ret) = optimize(optprob.opt, initial)

    (minf, minx)
end

"""Setup an optimization problem."""
function problem{T<:Real}(model::Model, components::Vector{Symbol}, names::Vector{Symbol}, lowers::Vector{T}, uppers::Vector{T}, objective::Function, objectiveconstraints::Vector{Function}, matrixconstraints::Vector{MatrixConstraintSet})
    my_lowers = T[]
    my_uppers = T[]

    ## Replace with eachname
    totalvars = 0
    for (ii, len, isscalar) in @task nameindexes(model, names)
        append!(my_lowers, [lowers[ii] for jj in 1:len])
        append!(my_uppers, [uppers[ii] for jj in 1:len])
        totalvars += len
    end

    LinprogOptimizationProblem(model, components, names, objective, objectiveconstraints, matrixconstraints, my_lowers, my_uppers)
end

"""Solve an optimization problem."""
function solution(optprob::LinprogOptimizationProblem, verbose=false)
    global allverbose
    allverbose = verbose

    initial = make0(optprob.model, optprob.names)

    if verbose
        println("Optimizing...")
    end

    if optprob.model.numberType == Number
        myobjective = unaryobjective(optprob.model, optprob.components, optprob.names, optprob.objective)
        f, b, A = lpconstraints(optprob.model, optprob.components, optprob.names, myobjective, objectiveconstraints)
    else
        f, b, A = lpconstraints(optprob.model, optprob.components, optprob.names, optprob.objective, optprob.objectiveconstraints)
    end

    f, b, A = combineconstraints(f, b, A, optprob.model, optprob.components, optprob.names, optprob.matrixconstraints)
    exlowers, exuppers = combinelimits(optprob.exlowers, optprob.exuppers, optprob.model, optprob.components, optprob.names, optprob.matrixconstraints)

    # Use -f, because linprog *minimizes* objective
    @time sol = linprog(-f, A, '<', b, optprob.exlowers, optprob.exuppers)

    sol.sol
end

end # module