text:
################
## Easter Dates
################

# easter dates from http://www.maa.mhn.de/StarDate/publ_holidays.html
@test_throws ErrorException BusinessDays.easter_date(Dates.Year(1100))
@test BusinessDays.easter_date(Dates.Year(1901)) == Dates.Date(1901, 04, 07)
@test BusinessDays.easter_date(Dates.Year(1902)) == Dates.Date(1902, 03, 30)
@test BusinessDays.easter_date(Dates.Year(1903)) == Dates.Date(1903, 04, 12)
@test BusinessDays.easter_date(Dates.Year(1904)) == Dates.Date(1904, 04, 03)
@test BusinessDays.easter_date(Dates.Year(1905)) == Dates.Date(1905, 04, 23)
@test BusinessDays.easter_date(Dates.Year(1906)) == Dates.Date(1906, 04, 15)
@test BusinessDays.easter_date(Dates.Year(1907)) == Dates.Date(1907, 03, 31)
@test BusinessDays.easter_date(Dates.Year(1908)) == Dates.Date(1908, 04, 19)
@test BusinessDays.easter_date(Dates.Year(1909)) == Dates.Date(1909, 04, 11)
@test BusinessDays.easter_date(Dates.Year(1910)) == Dates.Date(1910, 03, 27)
@test BusinessDays.easter_date(Dates.Year(1911)) == Dates.Date(1911, 04, 16)
@test BusinessDays.easter_date(Dates.Year(1912)) == Dates.Date(1912, 04, 07)
@test BusinessDays.easter_date(Dates.Year(1913)) == Dates.Date(1913, 03, 23)
@test BusinessDays.easter_date(Dates.Year(1914)) == Dates.Date(1914, 04, 12)
@test BusinessDays.easter_date(Dates.Year(1915)) == Dates.Date(1915, 04, 04)
@test BusinessDays.easter_date(Dates.Year(1916)) == Dates.Date(1916, 04, 23)
@test BusinessDays.easter_date(Dates.Year(1917)) == Dates.Date(1917, 04, 08)
@test BusinessDays.easter_date(Dates.Year(1918)) == Dates.Date(1918, 03, 31)
@test BusinessDays.easter_date(Dates.Year(1919)) == Dates.Date(1919, 04, 20)
@test BusinessDays.easter_date(Dates.Year(1920)) == Dates.Date(1920, 04, 04)
@test BusinessDays.easter_date(Dates.Year(1921)) == Dates.Date(1921, 03, 27)
@test BusinessDays.easter_date(Dates.Year(1922)) == Dates.Date(1922, 04, 16)
@test BusinessDays.easter_date(Dates.Year(1923)) == Dates.Date(1923, 04, 01)
@test BusinessDays.easter_date(Dates.Year(1924)) == Dates.Date(1924, 04, 20)
@test BusinessDays.easter_date(Dates.Year(1925)) == Dates.Date(1925, 04, 12)
@test BusinessDays.easter_date(Dates.Year(1926)) == Dates.Date(1926, 04, 04)
@test BusinessDays.easter_date(Dates.Year(1927)) == Dates.Date(1927, 04, 17)
@test BusinessDays.easter_date(Dates.Year(1928)) == Dates.Date(1928, 04, 08)
@test BusinessDays.easter_date(Dates.Year(1929)) == Dates.Date(1929, 03, 31)
@test BusinessDays.easter_date(Dates.Year(1930)) == Dates.Date(1930, 04, 20)
@test BusinessDays.easter_date(Dates.Year(1931)) == Dates.Date(1931, 04, 05)
@test BusinessDays.easter_date(Dates.Year(1932)) == Dates.Date(1932, 03, 27)
@test BusinessDays.easter_date(Dates.Year(1933)) == Dates.Date(1933, 04, 16)
@test BusinessDays.easter_date(Dates.Year(1934)) == Dates.Date(1934, 04, 01)
@test BusinessDays.easter_date(Dates.Year(1935)) == Dates.Date(1935, 04, 21)
@test BusinessDays.easter_date(Dates.Year(1936)) == Dates.Date(1936, 04, 12)
@test BusinessDays.easter_date(Dates.Year(1937)) == Dates.Date(1937, 03, 28)
@test BusinessDays.easter_date(Dates.Year(1938)) == Dates.Date(1938, 04, 17)
@test BusinessDays.easter_date(Dates.Year(1939)) == Dates.Date(1939, 04, 09)
@test BusinessDays.easter_date(Dates.Year(1940)) == Dates.Date(1940, 03, 24)
@test BusinessDays.easter_date(Dates.Year(1941)) == Dates.Date(1941, 04, 13)
@test BusinessDays.easter_date(Dates.Year(1942)) == Dates.Date(1942, 04, 05)
@test BusinessDays.easter_date(Dates.Year(1943)) == Dates.Date(1943, 04, 25)
@test BusinessDays.easter_date(Dates.Year(1944)) == Dates.Date(1944, 04, 09)
@test BusinessDays.easter_date(Dates.Year(1945)) == Dates.Date(1945, 04, 01)
@test BusinessDays.easter_date(Dates.Year(1946)) == Dates.Date(1946, 04, 21)
@test BusinessDays.easter_date(Dates.Year(1947)) == Dates.Date(1947, 04, 06)
@test BusinessDays.easter_date(Dates.Year(1948)) == Dates.Date(1948, 03, 28)
@test BusinessDays.easter_date(Dates.Year(1949)) == Dates.Date(1949, 04, 17)
@test BusinessDays.easter_date(Dates.Year(1950)) == Dates.Date(1950, 04, 09)
@test BusinessDays.easter_date(Dates.Year(1951)) == Dates.Date(1951, 03, 25)
@test BusinessDays.easter_date(Dates.Year(1952)) == Dates.Date(1952, 04, 13)
@test BusinessDays.easter_date(Dates.Year(1953)) == Dates.Date(1953, 04, 05)
@test BusinessDays.easter_date(Dates.Year(1954)) == Dates.Date(1954, 04, 18)
@test BusinessDays.easter_date(Dates.Year(1955)) == Dates.Date(1955, 04, 10)
@test BusinessDays.easter_date(Dates.Year(1956)) == Dates.Date(1956, 04, 01)
@test BusinessDays.easter_date(Dates.Year(1957)) == Dates.Date(1957, 04, 21)
@test BusinessDays.easter_date(Dates.Year(1958)) == Dates.Date(1958, 04, 06)
@test BusinessDays.easter_date(Dates.Year(1959)) == Dates.Date(1959, 03, 29)
@test BusinessDays.easter_date(Dates.Year(1960)) == Dates.Date(1960, 04, 17)
@test BusinessDays.easter_date(Dates.Year(1961)) == Dates.Date(1961, 04, 02)
@test BusinessDays.easter_date(Dates.Year(1962)) == Dates.Date(1962, 04, 22)
@test BusinessDays.easter_date(Dates.Year(1963)) == Dates.Date(1963, 04, 14)
@test BusinessDays.easter_date(Dates.Year(1964)) == Dates.Date(1964, 03, 29)
@test BusinessDays.easter_date(Dates.Year(1965)) == Dates.Date(1965, 04, 18)
@test BusinessDays.easter_date(Dates.Year(1966)) == Dates.Date(1966, 04, 10)
@test BusinessDays.easter_date(Dates.Year(1967)) == Dates.Date(1967, 03, 26)
@test BusinessDays.easter_date(Dates.Year(1968)) == Dates.Date(1968, 04, 14)
@test BusinessDays.easter_date(Dates.Year(1969)) == Dates.Date(1969, 04, 06)
@test BusinessDays.easter_date(Dates.Year(1970)) == Dates.Date(1970, 03, 29)
@test BusinessDays.easter_date(Dates.Year(1971)) == Dates.Date(1971, 04, 11)
@test BusinessDays.easter_date(Dates.Year(1972)) == Dates.Date(1972, 04, 02)
@test BusinessDays.easter_date(Dates.Year(1973)) == Dates.Date(1973, 04, 22)
@test BusinessDays.easter_date(Dates.Year(1974)) == Dates.Date(1974, 04, 14)
@test BusinessDays.easter_date(Dates.Year(1975)) == Dates.Date(1975, 03, 30)
@test BusinessDays.easter_date(Dates.Year(1976)) == Dates.Date(1976, 04, 18)
@test BusinessDays.easter_date(Dates.Year(1977)) == Dates.Date(1977, 04, 10)
@test BusinessDays.easter_date(Dates.Year(1978)) == Dates.Date(1978, 03, 26)
@test BusinessDays.easter_date(Dates.Year(1979)) == Dates.Date(1979, 04, 15)
@test BusinessDays.easter_date(Dates.Year(1980)) == Dates.Date(1980, 04, 06)
@test BusinessDays.easter_date(Dates.Year(1981)) == Dates.Date(1981, 04, 19)
@test BusinessDays.easter_date(Dates.Year(1982)) == Dates.Date(1982, 04, 11)
@test BusinessDays.easter_date(Dates.Year(1983)) == Dates.Date(1983, 04, 03)
@test BusinessDays.easter_date(Dates.Year(1984)) == Dates.Date(1984, 04, 22)
@test BusinessDays.easter_date(Dates.Year(1985)) == Dates.Date(1985, 04, 07)
@test BusinessDays.easter_date(Dates.Year(1986)) == Dates.Date(1986, 03, 30)
@test BusinessDays.easter_date(Dates.Year(1987)) == Dates.Date(1987, 04, 19)
@test BusinessDays.easter_date(Dates.Year(1988)) == Dates.Date(1988, 04, 03)
@test BusinessDays.easter_date(Dates.Year(1989)) == Dates.Date(1989, 03, 26)
@test BusinessDays.easter_date(Dates.Year(1990)) == Dates.Date(1990, 04, 15)
@test BusinessDays.easter_date(Dates.Year(1991)) == Dates.Date(1991, 03, 31)
@test BusinessDays.easter_date(Dates.Year(1992)) == Dates.Date(1992, 04, 19)
@test BusinessDays.easter_date(Dates.Year(1993)) == Dates.Date(1993, 04, 11)
@test BusinessDays.easter_date(Dates.Year(1994)) == Dates.Date(1994, 04, 03)
@test BusinessDays.easter_date(Dates.Year(1995)) == Dates.Date(1995, 04, 16)
@test BusinessDays.easter_date(Dates.Year(1996)) == Dates.Date(1996, 04, 07)
@test BusinessDays.easter_date(Dates.Year(1997)) == Dates.Date(1997, 03, 30)
@test BusinessDays.easter_date(Dates.Year(1998)) == Dates.Date(1998, 04, 12)
@test BusinessDays.easter_date(Dates.Year(1999)) == Dates.Date(1999, 04, 04)
@test BusinessDays.easter_date(Dates.Year(2000)) == Dates.Date(2000, 04, 23)
@test BusinessDays.easter_date(Dates.Year(2001)) == Dates.Date(2001, 04, 15)
@test BusinessDays.easter_date(Dates.Year(2002)) == Dates.Date(2002, 03, 31)
@test BusinessDays.easter_date(Dates.Year(2003)) == Dates.Date(2003, 04, 20)
@test BusinessDays.easter_date(Dates.Year(2004)) == Dates.Date(2004, 04, 11)
@test BusinessDays.easter_date(Dates.Year(2005)) == Dates.Date(2005, 03, 27)
@test BusinessDays.easter_date(Dates.Year(2006)) == Dates.Date(2006, 04, 16)
@test BusinessDays.easter_date(Dates.Year(2007)) == Dates.Date(2007, 04, 08)
@test BusinessDays.easter_date(Dates.Year(2008)) == Dates.Date(2008, 03, 23)
@test BusinessDays.easter_date(Dates.Year(2009)) == Dates.Date(2009, 04, 12)
@test BusinessDays.easter_date(Dates.Year(2010)) == Dates.Date(2010, 04, 04)
@test BusinessDays.easter_date(Dates.Year(2011)) == Dates.Date(2011, 04, 24)
@test BusinessDays.easter_date(Dates.Year(2012)) == Dates.Date(2012, 04, 08)
@test BusinessDays.easter_date(Dates.Year(2013)) == Dates.Date(2013, 03, 31)
@test BusinessDays.easter_date(Dates.Year(2014)) == Dates.Date(2014, 04, 20)
@test BusinessDays.easter_date(Dates.Year(2015)) == Dates.Date(2015, 04, 05)
@test BusinessDays.easter_date(Dates.Year(2016)) == Dates.Date(2016, 03, 27)
@test BusinessDays.easter_date(Dates.Year(2017)) == Dates.Date(2017, 04, 16)
@test BusinessDays.easter_date(Dates.Year(2018)) == Dates.Date(2018, 04, 01)
@test BusinessDays.easter_date(Dates.Year(2019)) == Dates.Date(2019, 04, 21)
@test BusinessDays.easter_date(Dates.Year(2020)) == Dates.Date(2020, 04, 12)
@test BusinessDays.easter_date(Dates.Year(2021)) == Dates.Date(2021, 04, 04)
@test BusinessDays.easter_date(Dates.Year(2022)) == Dates.Date(2022, 04, 17)
@test BusinessDays.easter_date(Dates.Year(2023)) == Dates.Date(2023, 04, 09)
@test BusinessDays.easter_date(Dates.Year(2024)) == Dates.Date(2024, 03, 31)
@test BusinessDays.easter_date(Dates.Year(2025)) == Dates.Date(2025, 04, 20)
@test BusinessDays.easter_date(Dates.Year(2026)) == Dates.Date(2026, 04, 05)
@test BusinessDays.easter_date(Dates.Year(2027)) == Dates.Date(2027, 03, 28)
@test BusinessDays.easter_date(Dates.Year(2028)) == Dates.Date(2028, 04, 16)
@test BusinessDays.easter_date(Dates.Year(2029)) == Dates.Date(2029, 04, 01)
@test BusinessDays.easter_date(Dates.Year(2030)) == Dates.Date(2030, 04, 21)
@test BusinessDays.easter_date(Dates.Year(2031)) == Dates.Date(2031, 04, 13)
@test BusinessDays.easter_date(Dates.Year(2032)) == Dates.Date(2032, 03, 28)
@test BusinessDays.easter_date(Dates.Year(2033)) == Dates.Date(2033, 04, 17)
@test BusinessDays.easter_date(Dates.Year(2034)) == Dates.Date(2034, 04, 09)
@test BusinessDays.easter_date(Dates.Year(2035)) == Dates.Date(2035, 03, 25)
@test BusinessDays.easter_date(Dates.Year(2036)) == Dates.Date(2036, 04, 13)
@test BusinessDays.easter_date(Dates.Year(2037)) == Dates.Date(2037, 04, 05)
@test BusinessDays.easter_date(Dates.Year(2038)) == Dates.Date(2038, 04, 25)
@test BusinessDays.easter_date(Dates.Year(2039)) == Dates.Date(2039, 04, 10)
@test BusinessDays.easter_date(Dates.Year(2040)) == Dates.Date(2040, 04, 01)
@test BusinessDays.easter_date(Dates.Year(2041)) == Dates.Date(2041, 04, 21)
@test BusinessDays.easter_date(Dates.Year(2042)) == Dates.Date(2042, 04, 06)
@test BusinessDays.easter_date(Dates.Year(2043)) == Dates.Date(2043, 03, 29)
@test BusinessDays.easter_date(Dates.Year(2044)) == Dates.Date(2044, 04, 17)
@test BusinessDays.easter_date(Dates.Year(2045)) == Dates.Date(2045, 04, 09)
@test BusinessDays.easter_date(Dates.Year(2046)) == Dates.Date(2046, 03, 25)
@test BusinessDays.easter_date(Dates.Year(2047)) == Dates.Date(2047, 04, 14)
@test BusinessDays.easter_date(Dates.Year(2048)) == Dates.Date(2048, 04, 05)
@test BusinessDays.easter_date(Dates.Year(2049)) == Dates.Date(2049, 04, 18)
@test BusinessDays.easter_date(Dates.Year(2050)) == Dates.Date(2050, 04, 10)
@test BusinessDays.easter_date(Dates.Year(2051)) == Dates.Date(2051, 04, 02)
@test BusinessDays.easter_date(Dates.Year(2052)) == Dates.Date(2052, 04, 21)
@test BusinessDays.easter_date(Dates.Year(2053)) == Dates.Date(2053, 04, 06)
@test BusinessDays.easter_date(Dates.Year(2054)) == Dates.Date(2054, 03, 29)
@test BusinessDays.easter_date(Dates.Year(2055)) == Dates.Date(2055, 04, 18)
@test BusinessDays.easter_date(Dates.Year(2056)) == Dates.Date(2056, 04, 02)
@test BusinessDays.easter_date(Dates.Year(2057)) == Dates.Date(2057, 04, 22)
@test BusinessDays.easter_date(Dates.Year(2058)) == Dates.Date(2058, 04, 14)
@test BusinessDays.easter_date(Dates.Year(2059)) == Dates.Date(2059, 03, 30)
@test BusinessDays.easter_date(Dates.Year(2060)) == Dates.Date(2060, 04, 18)
@test BusinessDays.easter_date(Dates.Year(2061)) == Dates.Date(2061, 04, 10)
@test BusinessDays.easter_date(Dates.Year(2062)) == Dates.Date(2062, 03, 26)
@test BusinessDays.easter_date(Dates.Year(2063)) == Dates.Date(2063, 04, 15)
@test BusinessDays.easter_date(Dates.Year(2064)) == Dates.Date(2064, 04, 06)
@test BusinessDays.easter_date(Dates.Year(2065)) == Dates.Date(2065, 03, 29)
@test BusinessDays.easter_date(Dates.Year(2066)) == Dates.Date(2066, 04, 11)
@test BusinessDays.easter_date(Dates.Year(2067)) == Dates.Date(2067, 04, 03)
@test BusinessDays.easter_date(Dates.Year(2068)) == Dates.Date(2068, 04, 22)
@test BusinessDays.easter_date(Dates.Year(2069)) == Dates.Date(2069, 04, 14)
@test BusinessDays.easter_date(Dates.Year(2070)) == Dates.Date(2070, 03, 30)
@test BusinessDays.easter_date(Dates.Year(2071)) == Dates.Date(2071, 04, 19)
@test BusinessDays.easter_date(Dates.Year(2072)) == Dates.Date(2072, 04, 10)
@test BusinessDays.easter_date(Dates.Year(2073)) == Dates.Date(2073, 03, 26)
@test BusinessDays.easter_date(Dates.Year(2074)) == Dates.Date(2074, 04, 15)
@test BusinessDays.easter_date(Dates.Year(2075)) == Dates.Date(2075, 04, 07)
@test BusinessDays.easter_date(Dates.Year(2076)) == Dates.Date(2076, 04, 19)
@test BusinessDays.easter_date(Dates.Year(2077)) == Dates.Date(2077, 04, 11)
@test BusinessDays.easter_date(Dates.Year(2078)) == Dates.Date(2078, 04, 03)


text:
module plot_utils
using Statistics
using Plots

function save_log(x::AbstractArray, logs::Matrix{Float64}, path::String)
    (_, imax) = size(logs)
    q25 = zeros(imax)
    q50 = zeros(imax)
    q75 = zeros(imax)
    for i in 1:imax
        q25[i] = quantile(logs[:, i], 0.25)
        q50[i] = quantile(logs[:, i], 0.50)
        q75[i] = quantile(logs[:, i], 0.75)
    end
    open(path, "w") do file
        for i in 1:imax
            write(file, string(Float64(x[i])," ", q25[i], " ", q50[i], " ", q75[i], "\n"))
        end
    end
end

function plot_log(paths::String...; xlims = (-Inf,Inf), ylims = (-Inf,Inf), symbols = 'A':'Z')
    p = plot(legend = :topright,
                 #xlabel = "fitness evals",
                 xlabel = "time (s)",
                 ylabel = "objective value",
                 yaxis =:log,
                 ylims = (1e-16, 1e+4))
    colors = Plots.palette(:tab10)
    for i in 1:length(paths)
        path = paths[i]
        denom = length(paths) > 1 ? string(symbols[i],": ") : ""
        open(path, "r") do file
            q25 = Vector{Float64}()
            q50 = Vector{Float64}()
            q75 = Vector{Float64}()
            t = Vector{Float64}()
            while !eof(file)
                (_t, _q25, _q50, _q75) = parse.(Float64, split(readline(file)))
                push!(t, _t)
                push!(q25, _q25)
                push!(q50, _q50)
                push!(q75, _q75)
            end
            yticks = [10.0^(-2*n) for n in -3:7]
            plot!(t, q50, label = denom*"median", lw = 2.0,
                xlims = xlims, ylims = ylims, yticks = yticks, linecolor = colors[i])
            plot!(t, q25, fillrange = q75, alpha = 0.3, lw = 2.0,
                label = denom*"first-third quantile",
                xlims = xlims, ylims = ylims, yticks = yticks,
                linecolor = colors[i], fillcolor = colors[i])
        end
    end
    return p
end

export save_log
export plot_log
end

ext: jl

max_stars_repo_name: OndrejKincl/evaTeaching-python

max_issues_repo_name: OndrejKincl/evaTeaching-python

max_forks_repo_name: OndrejKincl/evaTeaching-python


text:
mutable struct PBorn{T}
    attrib_mod::T
    medium::Medium
    ageom::AGeom
    srcwav::SrcWav
    fc::NamedStack{Float64}
    ic::NamedStack{Int64}
    cc::NamedStack{ComplexF64}
    fnpow2grid::FFTW.Frequencies{Float64}
    w::NamedStack{Vector{ComplexF64}}
    m::Vector{NamedStack{ComplexF64}}
    data::Records
end


text:
module Structures
# This module handles all functionality related to Structure and Atom
using ..DFControl
using ..Utils
using LinearAlgebra, StructTypes, Parameters, StaticArrays
using spglib_jll
const SPGLIB = spglib_jll.libsymspg

using UnitfulAtomic.Unitful: angstrom, Length, @unit, FreeUnits, unit, 𝐋, FreeUnits,
                             Quantity, ustrip, uconvert
const Ang = angstrom

const ReciprocalType{T,A} = Quantity{T,𝐋^-1,FreeUnits{A,𝐋^-1,nothing}}

StructTypes.StructType(::Type{<:Quantity}) = StructTypes.Struct()
@inline function StaticArrays._inv(::StaticArrays.Size{(3, 3)},
                                   A::SMatrix{3,3,LT}) where {LT<:Union{Length,
                                                                        ReciprocalType}}
    @inbounds x0 = SVector{3}(A[1], A[2], A[3])
    @inbounds x1 = SVector{3}(A[4], A[5], A[6])
    @inbounds x2 = SVector{3}(A[7], A[8], A[9])

    y0 = cross(x1, x2)
    d  = StaticArrays.bilinear_vecdot(x0, y0)
    x0 = x0 / d
    y0 = y0 / d
    y1 = cross(x2, x0)
    y2 = cross(x0, x1)

    @inbounds return SMatrix{3,3}((y0[1], y1[1], y2[1], y0[2], y1[2], y2[2], y0[3], y1[3],
                                   y2[3]))
end

include("projections.jl")
include("atom.jl")
include("structure.jl")

export Projection, Atom, Element, Structure, Pseudo, DFTU
export angstrom, Ang

end


text:
abstract type Integrator{DT} end

function integrate_step! end

struct ExplicitEuler{DT} <: Integrator{DT}
    Δt::DT
    ẋ::Vector{DT}

    function ExplicitEuler(equ::Equation{DT}, Δt::DT) where {DT}
        new{DT}(Δt, zeros(DT, ndims(equ)))
    end
end

function integrate_step!(int::ExplicitEuler, equ::Equation, x₀::AbstractVector, x₁::AbstractVector)
    equ.f(int.ẋ, x₀)
    x₁ .= x₀ .+ int.Δt .* int.ẋ
end

max_stars_repo_name: michakraus/JuliaAtHomeAssignments

max_issues_repo_name: michakraus/JuliaAtHomeAssignments

max_forks_repo_name: michakraus/JuliaAtHomeAssignments


text:
const VI = MOI.VariableIndex
const CI = MOI.ConstraintIndex
const SV = MOI.SingleVariable
const SAT = MOI.ScalarAffineTerm{Float64}
const SAF = MOI.ScalarAffineFunction{Float64}
const ET = MOI.EqualTo{Float64}
const LT = MOI.LessThan{Float64}
const GT = MOI.GreaterThan{Float64}
const IN = MOI.Interval{Float64}
const ZO = MOI.ZeroOne

function _build_base_model()
    # Model domain: x in [-1.1, 1.1]^2 satisfying
    #               x_1 + x_2 <= 0.5 AND
    #               x_1 - x_2 <= 0.6 AND
    #               x_1 is integer.
    model = MOIU.Model{Float64}()
    x_v, x_c =
        MOI.add_constrained_variables(model, [IN(-1.1, 1.1) for _ in 1:2])
    x = [SV(vi) for vi in x_v]
    MOI.add_constraint(model, 1.0 * x[1] + 1.0 * x[2], LT(0.5))
    MOI.add_constraint(model, 1.0 * x[1] - 1.0 * x[2], LT(0.6))
    MOI.add_constraint(model, x[1], MOI.Integer())
    return model, x
end

function _is_equal(u::SAF, v::SAF)
    u_t = sort(
        u.terms,
        lt = (x, y) -> x.variable_index.value < y.variable_index.value,
    )
    v_t = sort(
        v.terms,
        lt = (x, y) -> x.variable_index.value < y.variable_index.value,
    )
    if !(u.constant ≈ v.constant)
        return false
    end
    u_vars = [term.variable_index for term in u_t]
    v_vars = [term.variable_index for term in v_t]
    if u_vars != v_vars
        return false
    end
    u_coeffs = [term.coefficient for term in u_t]
    v_coeffs = [term.coefficient for term in v_t]
    @assert length(u_coeffs) == length(v_coeffs)
    for i in 1:length(u_coeffs)
        if !(u_coeffs[i] ≈ v_coeffs[i])
            return false
        end
    end
    return true
end

@testset "naive_big_m_formulation!" begin
    for (activity_method, expected_big_ms) in [
        DisjunctiveConstraints.IntervalArithmetic() => (2.2, 2.2, -2.15, -2.15),
        DisjunctiveConstraints.LinearProgrammingRelaxation(Cbc.Optimizer) =>
            (0.5, 0.6, -2.15, -2.15),
        DisjunctiveConstraints.FullFormulation(Cbc.Optimizer) =>
            (0.5, 0.6, -2.05, -2.05),
    ]
        method = DisjunctiveConstraints.NaiveBigM(activity_method)
        model, x = _build_base_model()

        # Disjunction: (x_1 + x_2 <= 0) OR
        #              (x_1 - x_2 <= 0) OR
        #              (x_1 + 0.5x_2 >= 0.5 AND x_1 - 0.5x_2 >= 0.5)
        f_1 = 1.0 * x[1] + 1.0 * x[2]
        f_2 = 1.0 * x[1] - 1.0 * x[2]
        f_3 = 1.0 * x[1] + 0.5 * x[2]
        f_4 = 1.0 * x[1] - 0.5 * x[2]
        f = MOIU.vectorize([f_1, f_2, f_3, f_4])

        lbs = [
            -Inf -Inf -Inf
            -Inf -Inf -Inf
            -Inf -Inf 0.5
            -Inf -Inf 0.5
        ]
        ubs = [
            0.0 Inf Inf
            Inf 0.0 Inf
            Inf Inf Inf
            Inf Inf Inf
        ]
        s = DisjunctiveConstraints.DisjunctiveSet(lbs, ubs)
        disjunction = DisjunctiveConstraints.Disjunction(f, s)

        state = @inferred DisjunctiveConstraints.formulate!(
            model,
            method,
            disjunction,
        )
        @test length(state.z_vis) == 3
        @test MOI.is_valid(model, state.sum_ci)
        @test size(state.lt_cis) == (4, 3)
        @test size(state.gt_cis) == (4, 3)

        @test MOI.get(model, MOI.NumberOfVariables()) == 5
        for i in 1:2
            @test MOIU.get_bounds(model, Float64, VI(i)) == (-1.1, 1.1)
        end
        for i in 3:5
            @test MOI.is_valid(model, MOI.ConstraintIndex{SV,ZO}(i))
        end
        z = [SV(VI(i)) for i in 3:5]
        moi_et_cis = MOI.get(model, MOI.ListOfConstraintIndices{SAF,ET}())
        @test length(moi_et_cis) == 1
        @test state.sum_ci == moi_et_cis[1]
        let f = MOI.get(model, MOI.ConstraintFunction(), state.sum_ci),
            s = MOI.get(model, MOI.ConstraintSet(), state.sum_ci)

            @test _is_equal(f, 1.0 * z[1] + 1.0 * z[2] + 1.0 * z[3])
            @test s == ET(1.0)
        end
        moi_lt_cis = MOI.get(model, MOI.ListOfConstraintIndices{SAF,LT}())
        # two from base model, two from big-M formulation
        @test length(moi_lt_cis) == 2 + 2
        @test count(!isnothing, state.lt_cis) == 2
        let lt_ci = state.lt_cis[1, 1]
            @test !isnothing(lt_ci)
            @test lt_ci in moi_lt_cis
            f = MOI.get(model, MOI.ConstraintFunction(), lt_ci)
            s = MOI.get(model, MOI.ConstraintSet(), lt_ci)
            @test _is_equal(
                f,
                1.0 * x[1] + 1.0 * x[2] + expected_big_ms[1] * z[1],
            )
            @test s isa MOI.LessThan
            @test s.upper ≈ 0.0 + expected_big_ms[1]
        end

        let lt_ci = state.lt_cis[2, 2]
            @test !isnothing(lt_ci)
            @test lt_ci in moi_lt_cis
            f = MOI.get(model, MOI.ConstraintFunction(), lt_ci)
            s = MOI.get(model, MOI.ConstraintSet(), lt_ci)
            @test _is_equal(
                f,
                1.0 * x[1] - 1.0 * x[2] + expected_big_ms[2] * z[2],
            )
            @test s isa MOI.LessThan
            @test s.upper ≈ 0.0 + expected_big_ms[2]
        end

        moi_lt_cis = MOI.get(model, MOI.ListOfConstraintIndices{SAF,GT}())
        # None from base model, two from big-M formulation
        @test length(moi_lt_cis) == 2
        @test count(!isnothing, state.gt_cis) == 2
        let gt_ci = state.gt_cis[3, 3]
            @test !isnothing(gt_ci)
            @test gt_ci in moi_lt_cis
            f = MOI.get(model, MOI.ConstraintFunction(), gt_ci)
            s = MOI.get(model, MOI.ConstraintSet(), gt_ci)
            @test _is_equal(
                f,
                1.0 * x[1] + 0.5 * x[2] + expected_big_ms[3] * z[3],
            )
            @test s isa MOI.GreaterThan
            @test s.lower ≈ 0.5 + expected_big_ms[3]
        end

        let gt_ci = state.gt_cis[4, 3]
            @test !isnothing(gt_ci)
            @test gt_ci in moi_lt_cis
            f = MOI.get(model, MOI.ConstraintFunction(), gt_ci)
            s = MOI.get(model, MOI.ConstraintSet(), gt_ci)
            @test _is_equal(
                f,
                1.0 * x[1] - 0.5 * x[2] + expected_big_ms[3] * z[3],
            )
            @test s isa MOI.GreaterThan
            @test s.lower ≈ 0.5 + expected_big_ms[3]
        end
    end
end


text:
# Autogenerated wrapper script for LibBirch_jll for x86_64-apple-darwin
export libbirch, libbirch_debug

using CompilerSupportLibraries_jll
using Eigen_jll
JLLWrappers.@generate_wrapper_header("LibBirch")
JLLWrappers.@declare_library_product(libbirch, "@rpath/libbirch.0.dylib")
JLLWrappers.@declare_library_product(libbirch_debug, "@rpath/libbirch-debug.0.dylib")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll, Eigen_jll)
    JLLWrappers.@init_library_product(
        libbirch,
        "lib/libbirch.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libbirch_debug,
        "lib/libbirch-debug.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# test partial pose3 constraints and evaluation

using RoME, Distributions
using IncrementalInference, TransformUtils
using Base: Test

N=50
fg = initfg()

v1 = addNode!(fg,:x1,0.001*randn(6,N),N=N)

mu1 = [0.0;0.0; -10.0]
prpz = PartialPriorRollPitchZ(
  MvNormal( mu1[1:2], [[1e-4;0.0]'; [0.0;1e-4]'] ),
  Normal( mu1[3], 0.0281 )
)

mu2 = [20.0,5.0,pi/2]
xyy = PartialPose3XYYaw(
  MvNormal(
    mu2,
    [1e-4 0.0 0.0; 0.0 1e-4 0.0; 0.0 0.0 4e-6]
  )
)

v2 = addNode!(fg,:x2,randn(6,N),N=N)

f1 = addFactor!(fg, [:x2], prpz)
f2 = addFactor!(fg, [:x1;:x2], xyy)

println("test PartialPriorRollPitchZ evaluations")

X2pts = getVal(v2)
pts = evalFactor2(fg, f1, v2.index, N=N)

newdims = collect(getData(f1).fnc.usrfnc!.partial)
olddims = setdiff(collect(1:6), newdims)

@test size(pts, 1) == 6
@test size(pts, 2) == N
@test norm(X2pts[olddims,:] - pts[olddims,:]) < 1e-10
# @show Base.mean(pts,2)[newdims]
@test norm(Base.mean(pts,2)[newdims]-mu1[[3;1;2]]) < 0.3
# ensure the correct response from
@test norm(X2pts[newdims,:] - pts[newdims,:]) > 2.0
# memcheck
@test norm(X2pts - getVal(v2)) < 1e-10

println("test residual function of PartialPose3XYYaw")

res = zeros(3)
idx = 1
meas = getSample(xyy)
xi = zeros(6,1)
xja = zeros(6,1)
xyy(res, idx, meas, xi, xja)
@test abs(res[1]-mu2[1]) < 0.2
@test abs(res[2]-mu2[2]) < 0.2
@test abs(res[3]-mu2[3]) < 0.2

xjb = zeros(6,1)
xjb[collect(xyy.partial),1] = mu2
res = zeros(3)
xyy(res, idx, meas, xi, xjb)
@test 0.0 < norm(res) < 0.2

meas = getSample(xyy,100)
@test norm(Base.std(meas[1],2)- [0.01;0.01;0.002]) < 0.005

println("test PartialPose3XYYaw evaluations")

X2pts = getVal(v2)
pts = evalFactor2(fg, f2, v2.index, N=N)

newdims = collect(getData(f2).fnc.usrfnc!.partial)
olddims = setdiff(collect(1:6), newdims)

@test size(pts, 1) == 6
@test size(pts, 2) == N
@test norm(X2pts[olddims,:] - pts[olddims,:]) < 1e-10

# @show Base.mean(getVal(v1),2)[newdims]
# @show mu2
# @show Base.mean(pts,2)[newdims]
# @show Base.std(pts,2)[newdims]
for i in 1:N
  pts[6,i] = wrapRad(pts[6,i])
end

@test norm(Base.mean(pts[newdims,:],2)-mu2) < 0.3
# ensure the correct response from
@test norm(X2pts[newdims,:] - pts[newdims,:]) > 2.0
memcheck = getVal(v2)
@test norm(X2pts - memcheck) < 1e-10

println("test predictbelief with two functions")
val = predictbelief(fg, :x2, ls(fg, :x2), N=N)

for i in 1:N
  val[6,i] = wrapRad(val[6,i])
end

@test size(val, 1) == 6
@test size(val, 2) == N

@test norm(Base.mean(val[collect(getData(f1).fnc.usrfnc!.partial),:],2)-mu1[[3;1;2]]) < 0.3
@test norm(Base.mean(val[collect(getData(f2).fnc.usrfnc!.partial),:],2)-mu2) < 0.3

memcheck = getVal(v2)
@test 1e-10 < norm(val - memcheck)

#


text:
function NOD(a,b)
    m = a
    n = b
    while n != 0
        m, n = n, m % n
    end
    return m
end

function zerodivisors(m)
    answer=Vector{Int}[]
    for i in 2:m-1
        if NOD(i,m)!=1
            push!(answer,i)
        end
    end
    return answer
end

zerodivisors(::ResidueRing{m}) where m = ResidueRing{m}.zerodivisors(m)


text:
using MathPlay
using Test

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


text:
using Test, Convex, SCS, SchattenNorms

import Base.kron
import LinearAlgebra

eye(n::Integer) = Matrix{Float64}(LinearAlgebra.I, (n, n))

x,y,z = [0 1; 1 0], [0 -1im; 1im 0], [1 0; 0 -1]

liou(x) = kron(conj(x),x)
liou(x,y) = kron(transpose(y),x)

function rand_unitary(d)
    rm = randn(d,d)+1im*randn(d,d)
    return LinearAlgebra.svd(rm).U
end

@testset "maximal dnorm examples" begin
    let # variable scope for cnot
        cnot = float([1 0 0 0; 0 1 0 0; 0 0 0 1; 0 0 1 0])
        cnotL = liou(cnot)
        @test isapprox(2.0, dnorm(cnotL-eye(16)), atol=1e-5)
        @test isapprox(2.0, ddist(cnotL,eye(16)), atol=1e-5)
        @test isapprox(2.0, ddistu(cnot,eye(4)), atol=1e-5)
    end

    @test isapprox(2.0, dnorm(eye(4)-liou(x)), atol=1e-5)
    @test isapprox(2.0, ddist(eye(4),liou(x)), atol=1e-5)
    @test isapprox(2.0, ddistu(eye(2),x), atol=1e-5)

    @test isapprox(2.0, dnorm(eye(4)-liou(z)), atol=1e-5)
    @test isapprox(2.0, ddist(eye(4),liou(z)), atol=1e-5)
    @test isapprox(2.0, ddistu(eye(2),z), atol=1e-5)

    @test isapprox(2.0, dnorm(eye(4)-liou(y)), atol=1e-5)
    @test isapprox(2.0, ddist(eye(4),liou(y)), atol=1e-5)
    @test isapprox(2.0, ddistu(eye(2),y), atol=1e-5)

    @test isapprox(2.0, dnorm(liou(x)-liou(y)), atol=1e-5)
    @test isapprox(2.0, ddist(liou(x),liou(y)), atol=1e-5)
    @test isapprox(2.0, ddistu(x,y), atol=1e-5)

    @test isapprox(2.0, dnorm(liou(x)-liou(z)), atol=1e-5)
    @test isapprox(2.0, ddist(liou(x),liou(z)), atol=1e-5)
    @test isapprox(2.0, ddistu(x,z), atol=1e-5)

    @test isapprox(2.0, dnorm(liou(z)-liou(y)), atol=1e-5)
    @test isapprox(2.0, ddist(liou(z),liou(y)), atol=1e-5)
    @test isapprox(2.0, ddistu(z,y), atol=1e-5)
end

@testset "maximal dnorm examples, random" begin
    global u = rand_unitary(2)
    global uu = liou(u)
    global v = rand_unitary(2)
    global vv = liou(v)

    @test isapprox(2.0, dnorm(uu*eye(4)-uu*liou(x)), atol=1e-5)
    @test isapprox(2.0, ddist(uu*eye(4),uu*liou(x)), atol=1e-5)
    @test isapprox(2.0, ddistu(u*eye(2),u*x), atol=1e-5)

    @test isapprox(2.0, dnorm(uu*eye(4)-uu*liou(z)), atol=1e-5)
    @test isapprox(2.0, ddist(uu*eye(4),uu*liou(z)), atol=1e-5)
    @test isapprox(2.0, ddistu(u*eye(2),u*z), atol=1e-5)

    @test isapprox(2.0, dnorm(uu*eye(4)-uu*liou(y)), atol=1e-5)
    @test isapprox(2.0, ddist(uu*eye(4),uu*liou(y)), atol=1e-5)
    @test isapprox(2.0, ddistu(u*eye(2),u*y), atol=1e-5)

    @test isapprox(2.0, dnorm(uu*liou(x)-uu*liou(y)), atol=1e-5)
    @test isapprox(2.0, ddist(uu*liou(x),uu*liou(y)), atol=1e-5)
    @test isapprox(2.0, ddistu(u*x,u*y), atol=1e-5)

    @test isapprox(2.0, dnorm(uu*liou(x)-uu*liou(z)), atol=1e-5)
    @test isapprox(2.0, ddist(uu*liou(x),uu*liou(z)), atol=1e-5)
    @test isapprox(2.0, ddistu(u*x,u*z), atol=1e-5)

    @test isapprox(2.0, dnorm(uu*liou(z)-uu*liou(y)), atol=1e-5)
    @test isapprox(2.0, ddist(uu*liou(z),uu*liou(y)), atol=1e-5)
    @test isapprox(2.0, ddistu(u*z,u*y), atol=1e-5)
end

@testset "dnorm for Pauli channels" begin
pauli = Matrix{ComplexF64}[eye(2), x, y, z]
    for i in 1:10
        p1 = rand(4); p1 = p1/sum(p1)
        p2 = rand(4); p2 = p2/sum(p2)
        pauli_dnorm = norm(p1-p2,1)

        pc1 = sum([ p1[ii+1]*liou(pauli[ii+1]) for ii in 0:3])
        pc2 = sum([ p2[ii+1]*liou(pauli[ii+1]) for ii in 0:3])

        calc_dnorm1 = ddist(pc1,pc2)

        @test isapprox(pauli_dnorm, calc_dnorm1, atol=1e-5)
    end
end

@testset "unitary invariance of dnorm" begin
    global u = rand_unitary(2)
    global uu = liou(u)
    global v = rand_unitary(2)
    global vv = liou(v)

    global duv  = dnorm(uu-vv)
    @time begin
        global duv2 = dnorm(uu*vv'-eye(4))
        global duv3 = dnorm(vv'*uu-eye(4))
        global duv4 = dnorm(eye(4)-uu'*vv)
        global duv5 = dnorm(eye(4)-vv*uu')
    end

    @test isapprox(duv, duv2, atol=1e-5)
    @test isapprox(duv, duv3, atol=1e-5)
    @test isapprox(duv, duv4, atol=1e-5)
    @test isapprox(duv, duv5, atol=1e-5)
end

@testset "ddist difference of unitary transformations" begin

    for i in 1:20
        for d in [2,3,4]
            global u = rand_unitary(d)
            global uu = liou(u)
            global v = rand_unitary(d)
            global vv = liou(v)
            @test isapprox(ddistu(u,v), ddist(uu,vv), atol=1e-4)
        end
    end

    for i in 1:20
        for d in [2,3,4]
            global u = rand_unitary(d)
            global uu = liou(u)
            global v = rand_unitary(d)
            global vv = liou(v)

            global duv  = ddistu(u,v)
            global duv2 = ddistu(u*v',eye(d))
            global duv3 = ddistu(v'u,eye(d))
            global duv4 = ddistu(eye(d),u'*v)
            global duv5 = ddistu(eye(d),v*u')

            @test isapprox(duv, duv2, atol=1e-10)
            @test isapprox(duv, duv3, atol=1e-10)
            @test isapprox(duv, duv4, atol=1e-10)
            @test isapprox(duv, duv5, atol=1e-10)
        end
    end
end

#def = dnorm(ee-ff)
#@time begin
#    def2 = dnorm(u*(ee-ff)*v)
#    def3 = dnorm(u*(ee-ff))
#    def4 = dnorm((ee-ff)*v)
#    def5 = dnorm(u'*(ee-ff))
#end

#println(def)
#println(def2)
#println(def3)
#println(def4)
#println(def5)


text:
"""
module that holds functions needed to  react to scrolling
Generally first we need to pass the GLFW callback to the Rocket obeservable
code adapted from https://discourse.julialang.org/t/custom-subject-in-rocket-jl-for-mouse-events-from-glfw/65133/3
"""
module ReactToScroll
using ModernGL, ..DisplayWords,Rocket, GLFW,  ..ForDisplayStructs,  ..TextureManag,Logging,  ..DataStructs,  ..StructsManag

export reactToScroll
export registerMouseScrollFunctions

"""
configuting Rocket on Subscribe so we get custom handler of input as we see we still need to define actor
"""
function Rocket.on_subscribe!(handler::ScrollCallbackSubscribable, actor::SyncActor{Any, ActorWithOpenGlObjects})
    return subscribe!(handler.subject, actor)
end

"""
we define how handler should act on the subject - observable so it will pass event onto subject
If we will scroll fast number will change much  and we will skip some slices
"""
function (handler::ScrollCallbackSubscribable)(_, xoff, yoff)

        if(!handler.isBusy[]) # if program is ready to repsond
           handler.numberToSend=0
           next!(handler.subject,Int64(handler.numberToSend+yoff))#true if we scroll up
        else
            handler.numberToSend+=yoff
        end
end

"""
uploading data to given texture; of given types associated
returns subscription in order to enable unsubscribing in the end
window - GLFW window
stopListening - atomic boolean able to stop the event listening cycle
return scrollback - that holds boolean subject (observable) to which we can react by subscribing appropriate actor
"""
function registerMouseScrollFunctions(window::GLFW.Window
                                    ,stopListening::Base.Threads.Atomic{Bool}
                                    ,isBusy::Base.Threads.Atomic{Bool}                   )

stopListening[]=true # stoping event listening loop to free the GLFW context

scrollback = ScrollCallbackSubscribable( isBusy,0 ,Subject(Int64, scheduler = AsyncScheduler()))
GLFW.SetScrollCallback(window, (a, xoff, yoff) -> scrollback(a, xoff, yoff))

stopListening[]=false # reactivate event listening loop

return scrollback

end #registerMouseScrollFunctions

"""
captures information send from handler that scroll was executed by the
"""

"""
in case of the scroll p true will be send in case of down - false
in response to it it sets new screen int variable and changes displayed screen
toBeSavedForBack - just marks weather we wat to save the info how to undo latest action
 - false if we invoke it from undoing
"""
function reactToScroll(scrollNumb::Int64
                    ,actor::SyncActor{Any, ActorWithOpenGlObjects}
                    ,toBeSavedForBack::Bool = true)
    actor.actor.mainForDisplayObjects.stopListening[]=true
    current = actor.actor.currentDisplayedSlice
    old = current
    #when shift is pressed scrolling is 10 times faster
    if(!actor.actor.mainForDisplayObjects.isFastScroll)
        current+=scrollNumb
    else
        current+=scrollNumb*10
    end

    #isScrollUp ? current+=1 : current-=1

   # we do not want to move outside of possible range of slices
   lastSlice = actor.actor.onScrollData.slicesNumber
   if(lastSlice>1)

    actor.actor.isSliceChanged = true
    actor.actor.isBusy[] = true
        if(current<1) current=1 end
        if(lastSlice<1) lastSlice=1 end
        if(current>=lastSlice) current=lastSlice end
        #logic to change displayed screen
        #we select slice that we are intrested in
        singleSlDat= actor.actor.onScrollData.dataToScroll|>
        (scrDat)-> map(threeDimDat->threeToTwoDimm(threeDimDat.type,Int64(current),actor.actor.onScrollData.dimensionToScroll,threeDimDat ),scrDat) |>
        (twoDimList)-> SingleSliceDat(listOfDataAndImageNames=twoDimList
                                    ,sliceNumber=current
                                    ,textToDisp = getTextForCurrentSlice(actor.actor.onScrollData, Int32(current))  )

        updateImagesDisplayed(singleSlDat
                            ,actor.actor.mainForDisplayObjects
                            ,actor.actor.textDispObj
                            ,actor.actor.calcDimsStruct
                            ,actor.actor.valueForMasToSet      )

        actor.actor.currentlyDispDat=singleSlDat
        # updating the last mouse position so when we will change plane it will better show actual position
        currentDim =Int64(actor.actor.onScrollData.dataToScrollDims.dimensionToScroll)
        lastMouse = actor.actor.lastRecordedMousePosition
        locArr = [lastMouse[1],lastMouse[2],lastMouse[3]]
        locArr[currentDim]= current
        actor.actor.lastRecordedMousePosition=CartesianIndex(locArr[1],locArr[2],locArr[3])
            #saving information about current slice for future reference
        actor.actor.currentDisplayedSlice = current
        #enable undoing the action
        if(toBeSavedForBack)
            func = ()->  reactToScroll(old-=scrollNumb, actor,false )
            addToforUndoVector(actor,func  )
        end

    end#if
   actor.actor.isBusy[] = false
   actor.actor.mainForDisplayObjects.stopListening[]=false

end#reactToScroll

end #ReactToScroll


text:
"""
    function fxmat(F)
---
If the length of categorical vector `F` is `n`, which has `v` leveles,
this function returns a `n × v` fixed effect matrix.
`μ` is thus not fitted.

In this simulation, the fixed effects are just the generation effects.
As selection on the binary trait goes on, the genetic values are elavated.
The challenge, however, is still killing half of the challenge population.

I may later replace this part using `MixedModels`.
"""
function fxmat(F)
    levels = unique(F)
    dic = begin
        tmp = Dict()
        col = 1
        for l in levels
            tmp[l] = col
            col += 1
        end
        tmp
    end
    nlvl = length(levels)
    nobs = length(F)
    X = zeros(Int16, nobs, nlvl) # I bet levels are no more than 32,767
    for i in 1:length(F)
        X[i, dic[F[i]]] = 1
    end
    X
end

"""
    function snp_blup(g, p, h²; Q = [], F = [], dd=0)
---
The plain SNP-BLUP procedure.
Returns E(population) and  a vector of SNP effects.

**Warning**: I use Float32 to save memory.
This will be changed if to run on a bigger machine.
"""
function snp_blup(g, p, h²; Q = [], F = [], dd=0)
    nlc, nid = size(g)
    λ = (1. - h²)/h² * nlc + dd
    X = begin           # incident matrix for fixed effects, include μ
        if length(F) > 0
            fxmat(F)
        else
            ones(Int16, nid, 1)
        end
    end
    nf = size(X)[2]
    sl = nlc + nf               # size of lhs

    lhs = zeros(Float32, sl, sl)
    # -- upper left block
    ul = view(lhs, 1:nf, 1:nf)
    matmul!(ul, X', X)
    # -- lower left block
    ll = view(lhs, nf+1:sl, 1:nf)
    matmul!(ll, g, X)
    # -- upper right block
    ur = view(lhs, 1:nf, nf+1:sl)
    copyto!(ur, ll')
    # -- lower right block
    lr = view(lhs, nf+1:sl, nf+1:sl)
    matmul!(lr, g, g')
    for i in nf+1:sl
        lhs[i, i] += λ
    end
    # -- if some QTL as fixed effects
    if length(Q) > 0
        l = Q .+ nf
        for x in l
            lhs[x, x] -= λ
        end
    end

    rhs = zeros(Float32, sl)
    y = reshape(p, length(p), :)
    rhs[1:nf]     = matmul(X', y)
    rhs[nf+1:end] = matmul(g, y)

    # Solving
    LAPACK.posv!('L', lhs, rhs)

    # return fixed effects and SNP effects separately
    rhs[1:nf], rhs[nf+1:end]
end


text:
using Dates
using HTTP
using HDF5
using LinearAlgebra
using MPIFiles
using Statistics
using Test
using UUIDs
using Unitful
using Scratch
using LazyArtifacts

const datadir = joinpath(artifact"data", "data")
@info "The test data is located at $datadir."

const tmpdir  = @get_scratch!("tmp")
@info "If you want to check the output of the tests, please head to $tmpdir."

mkpath(joinpath(tmpdir,"mdf"))
mkpath(joinpath(tmpdir,"positions"))
mkpath(joinpath(tmpdir,"transferFunction"))

include("DatasetStore.jl")
include("General.jl")
include("Cartesian.jl")
include("Positions.jl")
include("MDFv1.jl")
include("MultiMPIFile.jl")
include("Reco.jl")
include("IMT.jl")
include("TransferFunction.jl")
include("CustomSFMeas.jl")
include("MDFInMemory.jl")

@info "The unit tests are done!"


text:
module CUDAint

using ..CUDA

using CUDA: CUDNN
include("curnn.jl")
include("cudnn.jl")

end


text:
using LazyArrays, LinearAlgebra,  Test
import LazyArrays: ArrayLdivArrayStyle
import Base.Broadcast: materialize

@testset "Ldiv" begin
    @testset "Float64 \\ *" begin
        A = randn(5,5)
        b = randn(5)
        M = Mul(Inv(A),b)

        @test size(M) == (5,)
        @test similar(M) isa Vector{Float64}
        @test materialize(M) isa Vector{Float64}
        @test all(materialize(M) .=== (A\b))

        @test Base.BroadcastStyle(typeof(Ldiv(A,b))) isa ArrayLdivArrayStyle
        @test all(copyto!(similar(b), Ldiv(A,b)) .===
                    (similar(b) .= Ldiv(A,b)) .=== Inv(A) * b .===
                    materialize(Ldiv(A,b)) .===
                  (A\b) .=== (b̃ =  copy(b); LAPACK.gesv!(copy(A), b̃); b̃))

        @test copyto!(similar(b), Ldiv(UpperTriangular(A) , b)) ≈ UpperTriangular(A) \ b
        @test all(copyto!(similar(b), Ldiv(UpperTriangular(A) , b)) .===
                    (similar(b) .= Ldiv(UpperTriangular(A),b)) .=== Inv(UpperTriangular(A))*b .===
                    BLAS.trsv('U', 'N', 'N', A, b) )

        @test copyto!(similar(b), Ldiv(UpperTriangular(A)' , b)) ≈ UpperTriangular(A)' \ b
        @test all(copyto!(similar(b), Ldiv(UpperTriangular(A)' , b)) .===
                    (similar(b) .= Ldiv(UpperTriangular(A)',b)) .===
                    copyto!(similar(b), Ldiv(transpose(UpperTriangular(A)) , b)) .===
                            (similar(b) .= Ldiv(transpose(UpperTriangular(A)),b)) .===
                    BLAS.trsv('U', 'T', 'N', A, b))

        b = randn(5) + im*randn(5)
        @test Inv(A) * b ≈ Matrix(A) \ b
    end

    @testset "ComplexF64 \\ *" begin
        T = ComplexF64
        A = randn(T,5,5)
        b = randn(T,5)
        @test Base.BroadcastStyle(typeof(Ldiv(A,b))) isa ArrayLdivArrayStyle
        @test all(copyto!(similar(b), Ldiv(A,b)) .===
                    (similar(b) .= Ldiv(A,b)) .===
                  (A\b) .=== (b̃ =  copy(b); LAPACK.gesv!(copy(A), b̃); b̃))

        @test copyto!(similar(b), Ldiv(UpperTriangular(A) , b)) ≈ UpperTriangular(A) \ b
        @test all(copyto!(similar(b), Ldiv(UpperTriangular(A) , b)) .===
                    (similar(b) .= Ldiv(UpperTriangular(A),b)) .===
                    BLAS.trsv('U', 'N', 'N', A, b) )

        @test copyto!(similar(b), Ldiv(UpperTriangular(A)' , b)) ≈ UpperTriangular(A)' \ b
        @test all(copyto!(similar(b), Ldiv(UpperTriangular(A)' , b)) .===
                    (similar(b) .= Ldiv(UpperTriangular(A)',b)) .===
                    BLAS.trsv('U', 'C', 'N', A, b))

        @test copyto!(similar(b), Ldiv(transpose(UpperTriangular(A)) , b)) ≈ transpose(UpperTriangular(A)) \ b
        @test all(copyto!(similar(b), Ldiv(transpose(UpperTriangular(A)) , b)) .===
                            (similar(b) .= Ldiv(transpose(UpperTriangular(A)),b)) .===
                    BLAS.trsv('U', 'T', 'N', A, b))

        b = randn(5)
        @test Inv(A) * b ≈ A \ b
    end

    @testset "Triangular \\ matrix" begin
        A = randn(5,5)
        b = randn(5,5)
        M =  Ldiv(UpperTriangular(A), b)
        @test Base.Broadcast.broadcastable(M) === M
        @test UpperTriangular(A) \ b ≈ copyto!(similar(b) , Ldiv(UpperTriangular(A), b)) ≈
            (b .= Ldiv(UpperTriangular(A), b))
    end

    @testset "Rectangle PInv" begin
        A = randn(5,3)
        A_orig = copy(A)
        b = randn(5)
        b_orig = copy(b)
        @test_throws DimensionMismatch Inv(A)
        @test PInv(A) * b == (A\b)
        @test all(b .=== b_orig)
    end

    @testset "inv and pinv"  begin
        A = randn(5,5)
        @test inv(Inv(A)) === inv(PInv(A)) === pinv(Inv(A)) === pinv(PInv(A)) === A

        Ai = PInvMatrix(A)
        @test Matrix(Ai) ≈ inv(A)
        Ai = InvMatrix(A)
        @test Matrix(Ai) ≈ inv(A)

        A = randn(5,3)
        @test pinv(PInv(A)) === A
        Ai = PInvMatrix(A)
        @test Ai ≈ pinv(A)
        @test_throws DimensionMismatch inv(PInv(A))
        @test_throws DimensionMismatch InvMatrix(A)

        A = randn(3,5)
        @test pinv(PInv(A)) === A
        Ai = PInvMatrix(A)
        @test Ai ≈ pinv(A)
        @test_throws DimensionMismatch inv(PInv(A))
        @test_throws DimensionMismatch InvMatrix(A)
    end
end


text:
using Zygote: Grads

@testset "Params" begin
  @testset "delete!" begin
    w = rand(2,3)
    b = rand(2)
    ps = Params([w,b])
    delete!(ps, b)
    @test length(ps.order) == length(ps.params) == 1
    @test first(ps.order) == first(ps.params) == w
  end

  @testset "copy!" begin
    x = [0,0]
    ps = Params([x])
    copy!(ps, [1, 2])
    @test x == [1, 2]

    x = [0,0]
    y = [0]
    ps = Params([x, y])
    copy!(ps, [1, 2, 3])
    @test x == [1, 2]
    @test y == [3]

    ps = Params([[1,2]])
    x = [0, 0]
    copy!(x, ps)
    @test x == [1, 2]

    ps = Params([[1,2], [3]])
    x = [0, 0, 0]
    copy!(x, ps)
    @test x == [1, 2, 3]
  end

  @testset "broadcast" begin
    x, y = [1,2], [1]
    ps = Params([x, y])
    @test length.(ps) == length.([x, y]) # 617
    @test all(Params([[1,1]]) .== Params([[1,1]]))
  end

  @testset "indexing" begin
    x, y = [1,2], [1]
    ps = Params([x, y])
    @test ps[1] === x
    @test ps[2] === y
    @test ps[1:2] == [x, y]
  end

  @testset "comparison" begin
    x, y = [1,2], [1]
    ps1 = Params([x, y])
    ps2 = Params([x, y])
    ps3 = Params([y, x])
    @test ps1 == ps2
    @test ps1 != ps3  # comparison is order dependent
  end

  @testset "set interface" begin
    x, y, z = [1,2], [1], [3,4]
    ps1 = Params([x, y])
    ps2 = Params([z])
    ps3 = Params([y, z])

    ps = union(ps1, ps2)
    @test ps isa Params
    @test issetequal(ps, Set([x,y,z]))
    ps = union(ps1, ps3)
    @test ps isa Params
    @test issetequal(ps, Set([x,y,z]))

    ps = intersect(ps1, ps2)
    @test ps isa Params
    @test issetequal(ps, Set())

    ps = intersect(ps1, ps3)
    @test ps isa Params
    @test issetequal(ps, Set([y]))
  end
end

@testset "Grads" begin
  @testset "algebra" begin
    w, b = rand(2), rand(2)
    x1, x2 = rand(2), rand(2)

    gs1 = gradient(() -> sum(w .* x1), Params([w]))
    gs2 = gradient(() -> sum(w .* x2), Params([w]))

    @test .- gs1 isa Grads
    @test gs1 .- gs2 isa Grads
    @test .+ gs1 isa Grads
    @test gs1 .+ gs2 isa Grads
    @test 2 .* gs1 isa Grads
    @test (2 .* gs1)[w] ≈ 2 * gs1[w]
    @test gs1 .* 2 isa Grads
    @test gs1 ./ 2 isa Grads
    @test (gs1 .+ gs2)[w] ≈ gs1[w] .+ gs2[w]

    gs12 = gs1 .+ gs2
    gs1 .+= gs2
    @test gs12[w] ≈ gs1[w]

    gs3 = gradient(() -> sum(w .* x1), Params([w, b])) # grad nothing with respect to b
    gs4 = gradient(() -> sum(w .* x2 .+ b), Params([w, b]))

    @test .- gs3 isa Grads
    @test gs3 .- gs4 isa Grads
    @test .+ gs3 isa Grads
    @test gs3 .+ gs4 isa Grads
    @test 2 .* gs3 isa Grads
    @test gs3 .* 2 isa Grads
    @test gs3 ./ 2 isa Grads
    @test (gs3 .+ gs4)[w] ≈ gs3[w] .+ gs4[w]
    @test (gs3 .+ gs4)[b] ≈ gs4[b]

    @test gs3 .+ Dict(w => similar(w), b => similar(b)) isa Grads
    gs3 .+= Dict(p => randn(size(p)) for p in keys(gs3))
    @test gs3 isa Grads

  @test_throws ArgumentError gs1 .+ gs4
  end

  @testset "map and broadcast" begin
    w = rand(2)
    x1 = rand(2)
    x2 = rand(2)

    gs1 = gradient(() -> sum(w .* x1), Params([w]))
    gs2 = gradient(() -> sum(w .* x2), Params([w]))

    @test map(x -> zeros(2), gs1) isa Grads

    gs11 = map(x -> clamp.(x, -1e-5, 1e-5), gs1)
    @test gs11 isa Grads
    @test all(abs.(gs11[w]) .<= 1e-5)

    @test (x -> zeros(2)).(gs1) isa Grads
  end

  @testset "dictionary interface" begin
    w, b, x = rand(2), rand(2), rand(2)
    ps = Params([w, b])
    gs = gradient(() -> sum(tanh.(w .* x .+ b)), ps)

    @test issetequal(keys(gs), ps)
    @test length(values(gs)) == 2
    @test length(pairs(gs)) == 2
    k, v = first(pairs(gs))
    @test k === first(ps)
    @test v === gs[first(ps)]
  end

  @testset "iteration" begin
    w, b, x = rand(2), rand(2), rand(2)
    ps = Params([w, b])
    gs = gradient(() -> sum(tanh.(w .* x .+ b)), ps)

    # value-based iteration
    foreach(x -> clamp!(x, -1e-5, 1e-5), gs)
    @test all(abs.(gs[w]) .<= 1e-5)
    @test all(abs.(gs[b]) .<= 1e-5)
  end
end


text:
"""
Mutable type for storing medium parameters.
```julia
mod=Medium(mgrid, names=[:vp,:vs,:rho]) # initiate elastic medium
mod=Medium(mgrid, names=[:vp,:rho]) # initiate acoustic medium
```
This initializes a subsurface model with `:vp`, `:rho` and `:vs` parameters. Print the
names of the medium parameters stored in `mod`.
```julia
names(mod)
```

## Indexing
* `mod.mgrid` : returns the spatial-grid bundle
* `mod[:vp]` : P-wave velocity
* `mod[:vs]` : S-wave velocity
* `mod[:rho]` : mass density
* `mod[:Zp]` : P-wave impedance
* `mod[:K]` : bulk modulus (unrelaxed when considering attenuation)
* `mod[:M]` : P-wave modulus
* `mod[:mu]` : shear modulus
* `mod[:Q]` : quality factor (relaxation times are optimized to be constant over all frequencies; see Robertsson, et. al, 1994)
* `mod.ref` : reference medium parameters
* `mod.bounds` : bounds of medium parameters

"""
mutable struct Medium{N} # N is 2 for 2D, and 3 for 3D
    mgrid::Vector{T} where {T<:StepRangeLen{Float64}}
    m::NamedStack{Array{Float64,N}}
    bounds::NamedArrays.NamedArray{
        Array{Float64,1},
        1,
        Array{Array{Float64,1},1},
        Tuple{OrderedCollections.OrderedDict{Symbol,Int64}},
    }
    ref::NamedArray{
        Float64,
        1,
        Array{Float64,1},
        Tuple{OrderedCollections.OrderedDict{Symbol,Int64}},
    }
    # store some floats specific to the medium
    fc::NamedArrays.NamedArray{
        Float64,
        1,
        Array{Float64,1},
        Tuple{OrderedCollections.OrderedDict{Symbol,Int64}},
    }
    # store some integer constants
    ic::NamedArrays.NamedArray{
        Int64,
        1,
        Array{Int64,1},
        Tuple{OrderedCollections.OrderedDict{Symbol,Int64}},
    }
    # store 3D arrays, if any (currently used for attenuation, when N=2)
    m3::NamedStack{Array{Float64,3}}
end

# default construction without fcs and ics, and m3
function Medium(mgrid, m::NamedStack{Array{Float64,N}}, bounds, ref) where {N}
    return Medium(
        map(StepRangeLen{Float64},(mgrid)),
        Array{Float64}.(m),
        Array{Float64}.(bounds),
        Float64.(ref),
        NamedArray(Float64.([0.0]), ([:o],)),
        NamedArray([0], ([:o],)),
        NamedArray([zeros(Float64,1, 1, 1)], ([:o],)),
    )
end

function NamedArrays.names(mod::Medium)
    return names(mod.m)
end

include("attenuation.jl")
include("base.jl")
include("gallery.jl")

#=

* `vp0::Vector{Float64}` : [vpmin, vpmax]
* `vs0::Vector{Float64}` : [vsmin, vsmax]
* `rho0::Vector{Float64}` : [rhomin, rhomax]
* `χvp::Array{Float64,2}` : two-dimensional contrast model (χ) for vp, for e.g., zeros(length(mgrid[1]), length(mgrid[2]))
* `χvs::Array{Float64}` : two-dimensional contrast model (χ) for vs, for e.g., zeros(length(mgrid[1]), length(mgrid[2]))
* `χrho::Array{Float64}` : two-dimensional contrast model (χ) for density, for e.g., zeros(length(mgrid[1]), length(mgrid[2]))
* `mgrid` : array of ranges to determine the dimensions of models
"""
#
#	vp0::Vector{Float64}
#	vs0::Vector{Float64}
#	rho0::Vector{Float64}
#	χvp::Array{Float64,2}
#	χvs::Array{Float64,2}
#	χrho::Array{Float64,2}
#	# all derived fields from previous
#	K0::Vector{Float64}
#	mu0::Vector{Float64}
#	KI0::Vector{Float64}
#	muI0::Vector{Float64}
#	rhoI0::Vector{Float64}
#	"adding conditions that are to be false while construction"
#	Medium(vp0,vs0,rho0,χvp,χvs,χrho,mgrid,K0,mu0,KI0,muI0,rhoI0,ref) =
#		any([
#		     any(vp0.<0.0),
#		     any(vs0.<0.0),
#		     any(rho0.<0.0),
#		     #all([all(vp0 .≠ 0.0), any(χ(χvp,vp0,-1) .< vp0[1])]), # check vp bounds
#		     #all([all(vp0 .≠ 0.0), any(χ(χvp,vp0,-1) .> vp0[2])]), # check vp bounds
#		     #all([all(vs0 .≠ 0.0), any(χ(χvs,vs0,-1) .< vs0[1])]), # check vs bounds
#		     #all([all(vs0 .≠ 0.0), any(χ(χvs,vs0,-1) .> vs0[2])]), # check vs bounds
#		     #all([all(rho0 .≠ 0.0), any(χ(χrho,rho0,-1) .< rho0[1])]), # check rho bounds
#		     #all([all(rho0 .≠ 0.0), any(χ(χrho,rho0,-1) .> rho0[2])]), # check rho bounds
#		     size(χvp) != (length(mgrid[1]), length(mgrid[2])), # dimension check
#		     size(χvs) != (length(mgrid[1]), length(mgrid[2])), # dimension check
#		     size(χrho) != (length(mgrid[1]), length(mgrid[2])) # dimension check
#		    ]) ?
#		       error("error in Medium construction") : new(vp0,vs0,rho0,χvp,χvs,χrho,mgrid,K0,mu0,KI0,muI0,rhoI0,ref)
#end
#
# check whether a seismic model is bounded
function isbounded(mod::Medium)
	@warn "need to be implemented"
	result=true
end

# method to create Medium object without using derived fields
function Medium(vp0, vs0, rho0, χvp, χvs, χrho,
		 mgrid::Vector{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}})
	mod=Medium(vp0, vs0, rho0, χvp, χvs, χrho, mgrid,
	    zeros(2), zeros(2), zeros(2), zeros(2), zeros(2),
	    Medium_ref(zeros(8)...))
	update_derived!(mod)
	return mod
end

=#

# update the derived fields of the Seisimc model
function update!(mod::Medium)
    "note that: mean(K0) ≠ 1/mean(KI0)"
    mb = deepcopy(mod.bounds)

    rho = mb[:rho]
    rhoI = reverse(inv.(rho), dims = 1)

    if (:vs ∈ names(mb)[1])
        K = (mb[:vp] .* mb[:vp] .- 4.0 / 3.0 .* mb[:vs] .* mb[:vs]) .* mb[:rho]
        KI = reverse(inv.(K), dims = 1)
        newnames = [:K, :KI, :rhoI, :mu, :muI]
        mu = mb[:vs] .* mb[:vs] .* mb[:rho]
        muI = reverse(inv.(mu), dims = 1)
        mbnew = NamedArray([K, KI, rhoI, mu, muI], (newnames,))
    else
        K = (mb[:vp] .* mb[:vp]) .* mb[:rho]
        KI = reverse(inv.(K), dims = 1)
        newnames = [:K, :KI, :rhoI]
        mbnew = NamedArray([K, KI, rhoI], (newnames,))
    end

    for s in newnames
        # println(typeof(mod.bounds), '\'NamedArray([[0.0,0.0]],[[s],]))
        (s ∉ names(mod.bounds)[1]) &&
            (mod.bounds = Array{Float64}.(vcat(mod.bounds, NamedArray([[0.0, 0.0]], [s]))))
    end

    for s in newnames
        copyto!(mod.bounds[s], mbnew[s])
    end

    for s in names(mod.bounds)[1]
        meanvalue = Statistics.mean(mod.bounds[s])
        if (s ∈ names(mod.ref)[1])
            mod.ref[s] = meanvalue
        else
            mod.ref = vcat(mod.ref, NamedArray([meanvalue], [s]))
        end
    end

    return nothing
end

"""
Return medium property bounds based on maximum and minimum values of the array and frac.
The bounds cannot be less than zero
"""
function bounds(mod::AbstractArray, frac::Real = 0.1)
    any(mod .< 0.0) && error("model values less than zero")
    bounds = zeros(2)
    bound = frac * Statistics.mean(mod)
    bounds[1] = ((minimum(mod) - bound) < 0.0) ? 0.0 : (minimum(mod) - bound)
    bounds[2] = maximum(mod) + bound
    return bounds
end

function update!(mod::Medium, fields::AbstractVector{Symbol}, bounds)
    @assert !(bounds === nothing)
    for name in names(mod.bounds)[1]
        if (name ∈ fields)
            i = findall(x -> x == name, fields)[1]
            copyto!(mod.bounds[name], Float64.(bounds[i]))
        end
    end
    update!(mod)
end

"""
Adjust the bounds and hence the reference values.
Since the reference values are adjust the χ fields should also be changed
"""
function update!(mod::Medium, frac::Real)# 0.1
    for name in names(mod.m)[1]
        b = bounds(mod.m[name], frac)
        (name ∈ names(mod.bounds)[1]) && copyto!(mod.bounds[name], Float64.(b))
    end
    update!(mod)
    return mod
end

"""
Return dimensionless contrast model parameter using the
reference value.

# Arguments
* `mod::Array{Float64}` : subsurface parameter
* `mod0::Vector{Float64}` : reference value is mean of this vector
* `flag::Int64=1` :
"""
function χ(mod::AbstractArray{T}, mod0::T, flag::Int64 = 1) where {T<:Real}
    mod_out = copy(mod)
    χ!(mod_out, mod0, flag)
    return mod_out
end
function χ!(mod::AbstractArray{T}, mod0::T, flag::Int64 = 1) where {T<:Real}
    m0 = mod0
    if (flag == 1)
        if (iszero(m0))
            return nothing
        else
            @inbounds for i in eachindex(mod)
                mod[i] = ((mod[i] - m0) * inv(m0))
            end
            return mod
        end
    elseif (flag == -1)
        @inbounds for i in eachindex(mod)
            mod[i] = mod[i] * m0 + m0
        end
        return mod
    end
end

function χ(m::T, m0::T, flag::Int64 = 1) where {T<:Real}
    if (flag == 1)
        return (m - m0) * inv(m0)
    else
        (flag == -1)
        return (m * m0 + m0)
    end
end

"""
Gradients
Return contrast model parameter using the
reference value.
"""
function χg(mod::AbstractArray{T}, mod0::T, flag::Int64 = 1) where {T<:Real}
    mod_out = copy(mod)
    χg!(mod_out, mod0, flag)
    return mod_out
end
function χg!(mod::AbstractArray{T}, mod0::T, flag::Int64 = 1) where {T<:Real}
    if (flag == 1)
        rmul!(mod, mod0)
    elseif (flag == -1)
        m0 = inv(mod0)
        rmul!(mod, m0)
    end
    return mod
end

function χg(m::T, m0::T, flag::Int64 = 1) where {T<:Real}
    if (flag == 1)
        return m * m0
    else
        (flag == -1)
        return m * inv(m0)
    end
end

include("addons.jl")

include("chainrule.jl")

"""
Interpolate, extrapolate or truncate model to a new grid.
Check for memory allocations, specially when using Interpolations library.
Does this library contain inplace methods?
"""
function update!(modex::Medium{N}, mod::Medium{N}) where {N}
    for m in names(mod.m)[1]
        itp = extrapolate(
            interpolate(Tuple(mod.mgrid), mod[m], Gridded(Linear())),
            Interpolations.Flat(),
        )
        copyto!(modex[m], itp[modex.mgrid...])
    end
    update!(modex, 0.1)
    return modex
end

"""
Similar to previous method, but with initialization.
"""
function update(
    mod::Medium,
    mgrid::Vector{T} where {T<:StepRangeLen{Float64}}
)
    modex = Medium(mgrid, names(mod)[1]) # initialize a new medium
    update!(modex, mod)
    return modex
end

function padarray!(modex::Medium{N}, mod::Medium{N}, npml, edges) where {N}
    nex = length.(modex.mgrid)
    ist = [
        any(edges .== Symbol(string(dim), "min")) ? -npml + 1 : 1 for
        dim in dim_names(ndims(mod))
    ]
    vw = [ist[i]:ist[i]+nex[i]-1 for i = 1:ndims(mod)]
    pad = ImageFiltering.Pad(:replicate, fill(npml, ndims(mod))...)
    for m in names(mod.m)[1]
        m1 = ImageFiltering.BorderArray(mod[m], pad)
        m11 = view(m1, vw...)
        copyto!(modex[m], m11)
    end
    update!(modex, 0.1)
    return modex
end

# edges contains :xmin, :xmax, ... etc
function padarray(mod::Medium, npml, edges)
    minpads = [
        any(edges .== Symbol(string(dim), "min")) ? npml : 0 for
        dim in dim_names(ndims(mod))
    ]
    maxpads = [
        any(edges .== Symbol(string(dim), "max")) ? npml : 0 for
        dim in dim_names(ndims(mod))
    ]
    mgrid_new = [
        range(
            mg[1] - minpads[i] * step(mg),
            stop = mg[end] + maxpads[i] * step(mg),
            length = length(mg) + minpads[i] + maxpads[i],
        ) for (i, mg) in enumerate(mod.mgrid)
    ]
    modex = Medium(mgrid_new, names(mod)[1]) # initialize a new medium
    padarray!(modex, mod, npml, edges)
end

"""
Apply smoothing to `Medium` using a Gaussian filter of zwidth and xwidth

# Arguments

* `mod::Medium` : argument that is modified
* `zperc::Real` : smoothing percentage in z-direction
* `xperc::Real=zperc` : smoothing percentage in x-direction

# Keyword Arguments

* `zmin::Real=mod.mgrid[1][1]` :
* `zmax::Real=mod.mgrid[1][end]` :
* `xmin::Real=mod.mgrid[2][1]` :
* `xmax::Real=mod.mgrid[2][end]` :
* `fields` : fields of seismic model that are to be smooth
"""
function Medium_smooth(
    mod::Medium,
    zperc::Real,
    xperc::Real = zperc;
    zmin::Real = mod.mgrid[1][1],
    zmax::Real = mod.mgrid[1][end],
    xmin::Real = mod.mgrid[2][1],
    xmax::Real = mod.mgrid[2][end],
    fields = [:vp, :rho],
)
    xwidth = Float64(xperc) * 0.01 * abs(mod.mgrid[2][end] - mod.mgrid[2][1])
    zwidth = Float64(zperc) * 0.01 * abs(mod.mgrid[1][end] - mod.mgrid[1][1])
    xnwin = Int(div(xwidth, step(mod.mgrid[2]) * 2.0))
    znwin = Int(div(zwidth, step(mod.mgrid[1]) * 2.0))

    izmin = Interpolation.indminn(mod.mgrid[1], Float64(zmin), 1)[1]
    izmax = Interpolation.indminn(mod.mgrid[1], Float64(zmax), 1)[1]
    ixmin = Interpolation.indminn(mod.mgrid[2], Float64(xmin), 1)[1]
    ixmax = Interpolation.indminn(mod.mgrid[2], Float64(xmax), 1)[1]

    # @warn "check this routine, smooth contrast values instead?"

    modg = deepcopy(mod)
    for (i, iff) in enumerate(fields)
        m = view(mod[iff], izmin:izmax, ixmin:ixmax)
        mg = view(modg[iff], izmin:izmax, ixmin:ixmax)
        imfilter!(mg, m, Kernel.gaussian([znwin, xnwin]))
    end
    return modg
end

"""
function to resample in the model domain

# Arguments
* `mod::Medium` : model
* `modi::Medium` : model after interpolation
"""
function interp_spray!(
    mod::Medium,
    modi::Medium,
    attrib::Symbol,
    Battrib::Symbol = :B2,
    fields = [:vp, :rho];
    pa = nothing,
)
    if (pa === nothing)
        pa = Interpolation.Kernel(
            [mod.mgrid[2], mod.mgrid[1]],
            [modi.mgrid[2], modi.mgrid[1]],
            Battrib,
        )
    end

    # :vs is missing here....
    "loop over fields in `Medium`, add vs later"
    for field in fields
        Interpolation.interp_spray!(mod[field], modi[field], pa, attrib)
    end
end

#
#macro inter

function save(mod::Medium, folder; N = 100)
    !(isdir(folder)) && error("invalid directory")
    error("need to be updated")

    nx = length(mod.mgrid[2])
    nz = length(mod.mgrid[1])
    n = max(nx, nz)
    fact = (n > N) ? round(Int, n / N) : 1
    #mgrid=resamp(mod.mgrid, step(mod.mgrid[2])*fact, step(mod.mgrid[1])*fact)
    x = mgrid[2]
    z = mgrid[1]
    nx = length(x)
    nz = length(z)
    modo = Medium_zeros(mgrid)
    adjust_bounds!(modo, mod)
    interp_spray!(mod, modo, :interp)

    for m in [:vp, :rho, :Zp]
        # save original gf
        file = joinpath(folder, string("im", m, ".csv"))
        CSV.write(
            file,
            DataFrame(hcat(repeat(z, outer = nx), repeat(x, inner = nz), vec(modo[m]))),
        )
    end
end


text:
module GMT

using Printf

# Need to know what GMT version is available or if none at all to warn users on how to
# install GMT.
try
	#write("gmtversion__.txt", read(`gmt --version`))
	#global const GMTver = Meta.parse(read("gmtversion__.txt", String)[1:3])
	#rm("gmtversion__.txt")
	global const GMTver = Meta.parse(readlines(`gmt --version`)[1][1:3])
	global foundGMT = true
catch
	global foundGMT = false
	global const GMTver = 5.0
end

global legend_type = nothing
global img_mem_layout = ""			# "TCP"	 For Images.jl. The default is "TRBa"
global grd_mem_layout = ""			# "BRP" is the default for GMT PS images.
global current_cpt  = nothing		# To store the current palette
global current_view = nothing		# To store the current viewpoint (-p)
global multi_col    = false			# To allow plottig multiple columns at once
global IamModern    = false			# To know if we are in modern mode
global IamSubplot   = false			# To know if we are in subplot mode
global usedConfPar  = false			# Hacky solution for the session's memory trouble
const FMT = "ps"
const def_fig_size  = "12c/8c"              # Default fig size for plot like programs
const def_fig_axes  = " -Baf -BWSen"        # Default fig axes for plot like programs
const def_fig_axes3 = " -Baf -Bza -BWSenZ"  #		"" but for 3D views

export
	GMTver, FMT, gmt,
	arrows, arrows!, bar, bar!, bar3, bar3!, lines, lines!, legend, legend!,
	basemap, basemap!, blockmean, blockmedian, blockmode, clip, clip!, coast, coast!, colorbar, colorbar!,
	contour, contour!, filter1d, fitcircle, gmt2kml,  gmtconnect, gmtconvert, gmtinfo, gmtregress,
	gmtread, gmtselect, gmtset, gmtsimplify, gmtspatial, gmtvector, gmtwrite, gmtwhich,
	grd2cpt, grd2kml, grd2xyz, grdblend, grdclip, grdcontour, grdcontour!, grdcut, grdedit, grdfft,
	grdfilter, grdgradient, grdhisteq, grdimage, grdimage!, grdinfo, grdfill, grdlandmask, grdmask, grdpaste,
	grdproject, grdsample, grdtrack, grdtrend, grdvector, grdvector!, grdview, grdview!, grdvolume, greenspline,
	mat2ds, mat2grid, mat2img, histogram, histogram!, image, image!, imshow, kml2gmt, logo, logo!,
	makecpt, mask, mask!, mapproject, nearneighbor, plot, plot!, plot3d, plot3d!, project,
	pscontour, pscontour!, psconvert, psbasemap, psbasemap!, psclip, psclip!, pscoast, pscoast!,
	pshistogram, pshistogram!, psimage, psimage!, pslegend, pslegend!, psmask, psmask!, psrose, psrose!,
	psscale, psscale!, pssolar, pssolar!, psternary, psternary!, pstext, pstext!, pswiggle, pswiggle!,
	psxy, psxy!, psxyz, psxyz!, regress, rose, rose!, sample1d, scatter, scatter!, scatter3, scatter3!,
	solar, solar!, spectrum1d, sphdistance, sphinterpolate, sphtriangulate, surface, ternary, ternary!,
	text, text!, text_record, trend1d, trend2d, triangulate, splitxyz,
	decorated, vector_attrib, wiggle, wiggle!, xyz2grd,
	gmtbegin, gmtend, subplot, gmtfig, inset,
	linspace, logspace, contains, fields, tic, toc

include("common_docs.jl")
include("libgmt_h.jl")
include("libgmt.jl")
include("gmt_main.jl")
include("common_options.jl")
include("gmtbegin.jl")
include("blocks.jl")
include("filter1d.jl")
include("fitcircle.jl")
include("gmt2kml.jl")
include("gmtconnect.jl")
include("gmtconvert.jl")
include("gmtinfo.jl")
include("gmtlogo.jl")
include("gmtreadwrite.jl")
include("gmtset.jl")
include("gmtselect.jl")
include("gmtsimplify.jl")
include("gmtspatial.jl")
include("gmtregress.jl")
include("gmtvector.jl")
include("gmtwich.jl")
include("grd2cpt.jl")
include("grd2kml.jl")
include("grd2xyz.jl")
include("grdblend.jl")
include("grdclip.jl")
include("grdcut.jl")
include("grdedit.jl")
include("grdcontour.jl")
include("grdfft.jl")
include("grdfill.jl")
include("grdfilter.jl")
include("grdhisteq.jl")
include("grdinfo.jl")
include("grdimage.jl")
include("grdgradient.jl")
include("grdlandmask.jl")
include("grdmask.jl")
include("grdpaste.jl")
include("grdproject.jl")
include("grdsample.jl")
include("grdtrack.jl")
include("grdtrend.jl")
include("grdvector.jl")
include("grdview.jl")
include("grdvolume.jl")
include("greenspline.jl")
include("imshow.jl")
include("kml2gmt.jl")
include("makecpt.jl")
include("mapproject.jl")
include("nearneighbor.jl")
include("plot.jl")
include("project.jl")
include("psbasemap.jl")
include("psclip.jl")
include("pscoast.jl")
include("pscontour.jl")
include("psconvert.jl")
include("pshistogram.jl")
include("psimage.jl")
include("pslegend.jl")
include("psmask.jl")
include("psscale.jl")
include("psrose.jl")
include("pssolar.jl")
include("pstext.jl")
include("psxy.jl")
include("pswiggle.jl")
include("sample1d.jl")
include("spectrum1d.jl")
include("sphdistance.jl")
include("sphinterpolate.jl")
include("sphtriangulate.jl")
include("splitxyz.jl")
include("surface.jl")
include("subplot.jl")
include("triangulate.jl")
include("trend1d.jl")
include("trend2d.jl")
include("xyz2grd.jl")

if (!foundGMT)
	println("\n\nYou don't seem to have GMT installed and I don't currently install it automatically,")
	println("so you will have to do it yourself. Please follow instructions bellow but please note")
	println("that since GMT is migrating to Github (https://github.com/GenericMappingTools/gmt),")
	println("some of the links may change in a near future.\n\n")
	if (Sys.iswindows() && Sys.WORD_SIZE == 64)
		println("1) Download and install the official version at (the '..._win64.exe':")
		println("\t\t https://github.com/GenericMappingTools/gmt/releases")
	elseif (Sys.iswindows() && Sys.WORD_SIZE == 32)
		println("Download and install the official version at (the '..._win32.exe':")
		println("\t\t https://github.com/GenericMappingTools/gmt/releases")
	else
		println("https://github.com/GenericMappingTools/gmt/blob/master/INSTALL.md")
	end
end

function __init__()
	ver = Meta.parse(readlines(`gmt --version`)[1][1:3])
	if (ver != GMTver)
		println("\n\tYou seem to have changed your installed GMT version. Current version")
		println("\treports to be $ver but previously it was $GMTver. To fix this you need to")
		println("\tmake sure that GMT.jl recompiles again, otherwise fatal errors will occur.")
		t = joinpath(@__DIR__, "deps", "GMT.jl")
		println("\tA dirty way to force that is to make an irrelevant change in (and revert it after)\n\t$t")
	end
end

end # module


text:
using FunctionalCollections

immutable ArrayFunction
    expressions :: Dict{Unknown, ArrayExpr}
    parameters :: Array{Unknown, 1}
    returns :: Array{Unknown, 1}
    mutates :: Dict{Unknown, Unknown}
end

function all_nonleaf_nodes(aexpr)
    if is_leaf(aexpr)
        []
    else
        vcat([aexpr], map(all_nonleaf_nodes, children(aexpr))...)
    end
end

function build_parent_map(aexprs)
    # Only including non-leaves here!
    # We don't need to know leaves' parents, we won't replace them by aliases anyway.
    parents = Dict{ArrayExpr, Array{ArrayExpr, 1}}()
    for aexpr in aexprs
        for node in all_nonleaf_nodes(aexpr)
            for child in children(node)
                if is_node(child)
                    # There may be identical nodes that will be removed at merge stage
                    # For now we treat them separately
                    if !haskey(parents, child)
                        parents[child] = [aexpr]
                    else
                        push!(parents[child], aexpr)
                    end
                end
            end
        end

        # Add the top node to the mapping
        parents[aexpr] = []
    end
    parents
end

function merge(graph)
    parents = graph.parents # obj_id => [node]

    # Mark the nodes that can be merged
    original_nodes = Dict() # node => node
    to_merge = Dict() # obj_id => node
    for node in depth_first(graph.expressions)

        # Construct a replacement node, so that we could compare it with original ones
        # and decide right away whether it will require merging as well.
        # Since it's a depth-first search, we've already processed all children.
        node_children = children(node)
        child_ids = map(object_id, node_children)
        to_replace = Dict(
            child => original_nodes[to_merge[object_id(child)]]
            for child in node_children
            if object_id(child) in to_merge)
        new_node = replace_expr(node, to_replace)

        if new_node in original_nodes
            # `new_node` is equal to the original node, but it's a different object.
            # We need the original object.
            to_merge[object_id(node)] = original_nodes[new_node]
        else
            original_nodes[node] = node
        end
    end

    # Replace the marked nodes
    for (node, replacement) in reverse(to_merge)

    end
end

function array_function(returns, parameters, mutates=[])
    # Ideally, we can do the following backend-independent actions before generating code:
    # - merge expressions, so that we don't calculate the same thing twice
    # - for every node, calculate the number of nodes it is used in
    #   (those that are used in more than one node will have to be assigned to temporary vars)

    return_aliases = [Unknown(expr_type(expr)) for expr in returns]

    mutate_aliases = [Unknown(expr_type(expr)) for (expr, target) in mutates]
    mutate_targets = [target for (expr, target) in mutates]
    mutate_expressions = [expr for (expr, target) in mutates]
    mutate_map = Dict(
        alias => target
        for (alias, target) in zip(mutate_aliases, mutate_targets))

    aliases = vcat(return_aliases, mutate_aliases)
    expressions = vcat(returns, mutate_expressions)
    expression_map = Dict(
        alias => expr
        for (alias, expr) in zip(aliases, expressions))

    ArrayFunction(expression_map, parameters, return_aliases, mutate_map)
end

function (afunc::ArrayFunction)(args...)
    @assert length(afunc.mutates) == 0 # For now, only pure functions can be nested
    replacements = Dict(param => arg for (param, arg) in zip(afunc.parameters, args))
    returns = [replace_expr(afunc.expressions[ret], replacements)[2] for ret in afunc.returns]
    tuple(returns...)
end

function _replace_expr(
        used_aliases::PersistentSet{Unknown}, aexpr::ArrayExpr, replacements,
        # `true` on first call, so that we do not replace an expression that already has an alias,
        # but only replace its subexpressions
        is_root_expr)
    if !is_root_expr && haskey(replacements, aexpr)
        alias = replacements[aexpr]
        used_aliases = conj(used_aliases, alias)
        return used_aliases, convert(ArrayExpr, alias)
    elseif is_node(aexpr)
        map_accum(used_aliases, aexpr) do ua, arg
            _replace_expr(ua, arg, replacements, false)
        end
    else
        return used_aliases, aexpr
    end
end

function replace_expr(aexpr::ArrayExpr, replacements)
    used_aliases, new_expr = _replace_expr(
        PersistentSet{Unknown}(), aexpr, replacements, true)
    Set(used_aliases), new_expr
end

function separate(afunc)

    parents = build_parent_map([expr for (alias, expr) in afunc.expressions])

    mp_nodes_to_aliases = Dict(
        (node => Unknown(expr_type(node)))
        for (node, parent_nodes) in parents
        if length(parent_nodes) > 1)

    mp_aliases_to_nodes = map(reverse, mp_nodes_to_aliases)

    res = Base.merge(afunc.expressions, mp_aliases_to_nodes)

    aliases = Dict(
        (alias => replace_expr(expr, mp_nodes_to_aliases))
        for (alias, expr) in res)

    aliases
end

function serialize(serialized_aliases, alias, aliases_map)

    used_aliases, expr = aliases_map[alias]

    serialized = []
    for used_alias in used_aliases
        if !in(used_alias, serialized_aliases)
            serialized_aliases, l = serialize(serialized_aliases, used_alias, aliases_map)
            append!(serialized, l)
        end
    end

    serialized_aliases = conj(serialized_aliases, alias)
    push!(serialized, (alias => expr))

    serialized_aliases, serialized
end

function serialize(afunc)
# returns:
# - a list of pairs `(Unknown => ArrayExpr)` denoting "assignments" or "aliases".
#   The LHS arrays can appear in other expressions later in the list.
# These two are already contained in `afunc`:
# - a list of `Unknown` objects corresponding to return values
#   (each is either one of the parameters or one of the introduced aliases)
# - a dict `(Unknown => Unknown)`, where RHS is either a parameter or an alias,
#   and LHS is a parameter to which the LHS will be stored
# `Unknown` can be an array or a random state
    aliases_map = separate(afunc)

    serialized_aliases = PersistentSet{Unknown}()
    serialized = []
    for alias in keys(afunc.expressions)
        serialized_aliases, l = serialize(serialized_aliases, alias, aliases_map)
        append!(serialized, l)
    end

    return serialized
end

function mutated_parameters_indices(afunc::ArrayFunction)
    Set(find(
        # Cannot use `afunc.parameters .== target` instead, because the broadcasting
        # cannot be redirected into default functions in this case.
        # See the broadcasting methods in operators.jl.
        [param == target for param in afunc.parameters][1]
        for target in values(afunc.mutates)))
end


text:
module NumbersFromText

include("myparse.jl") # updates to parsing routines to work with Array{UInt8}

include("spacetokenizer.jl")
export SpaceTokenizer

include("readers.jl")
export readmatrix
export readarray, readarray!, readarrays, readarrays!

include("parallel.jl")
export partition_buffer, allocate_buffers

end


text:
# Autogenerated wrapper script for OpenLibm_jll for x86_64-apple-darwin14
export libopenlibm

JLLWrappers.@generate_wrapper_header("OpenLibm")
JLLWrappers.@declare_library_product(libopenlibm, "@rpath/libopenlibm.3.dylib")
JLLWrappers.@declare_library_product(libnonexisting, "libnonexisting.0.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libopenlibm,
        "lib/libopenlibm.3.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )
    JLLWrappers.@init_library_product(
        libnonexisting,
        "lib/libnonexisting.0.0.dylib",
        nothing,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()

max_stars_repo_name: giordano/JLLWrappers

max_issues_repo_name: giordano/JLLWrappers

max_forks_repo_name: giordano/JLLWrappers