text:
const TYPES_INT = Set([Int8, Int16, Int32, Int64, Int128])
const TYPES_UINT = Set([UInt8, UInt16, UInt32, UInt64, UInt128])
const ℤ = TYPES_INT ∪ TYPES_UINT
const ℝ = Set([Float16, Float32, Float64])

function assert_square_matrix_indices(i::Int, j::Int, N::Int)
  @assert i > 0 && j > 0 prerequisite("i > 0 & j > 0", ("i", i), ("j", j))
  @assert N > 0 prerequisite("N > 0", ("N", N))
  @assert i <= N && j <= N prerequisite("i ≤ N & j ≤ N", ("i", i), ("j", j), ("N", N))
end

function assert_elements_in_set(elements::Set, set::Set)
  if isempty(elements) return end
  @assert length(elements) <= length(set) prerequisite("|{elements}| ≤ |{set}|", ("elements", elements), ("set", set))
  for element in elements
    @assert element ∈ set prerequisite("element ∈ {set}", ("element", element), ("set", set))
  end
end

function prerequisite(text::AbstractString, pars::Tuple{Any, Any}...)
  details =  if length(pars) > 0
      "found: {$(join(map(p->"$(p[1])=$(p[2])", pars), ", "))}"
    else
      "please fix it."
    end
  "Prerequisite '$text', $details"
end

function exists(f::Function, collection::Set)
  return findfirst(f, collection) > 0
end

function exists(f::Function, collection::Vector)
  return findfirst(f, collection) > 0
end

function exists(f::Function, collection::Array)
  return findfirst(f, collection) > 0
end

function frexp10(value::Int)
  return frexp10(Float64(value))
end

function frexp10(value::Float64)
  v = value
  p = 0.0
  while v < 1.0
    v *= 10
    p -= 1.0
  end
  while v > 10.0
    v /= 10
    p += 1.0
  end
  return (v, p)
end

function weighted_mean_max(values::Vector{Float64})
  return weighted_mean(values; weightf=it->it[1])
end

function weighted_mean(values::Vector{Float64}; weightf::Function=it::Tuple{Int, Float64}->1.0)
  e_values = enumerate(values)
  return mapreduce(it->weightf(it)*it[2], +, e_values)/mapreduce(it->weightf(it), +, e_values)
end

function mean_base10(values::Float64...; meanf::Function = mean)
  vals = collect(values)
  return mean_base10(vals; meanf = meanf)
end

function mean_base10(values::Vector{Float64}; meanf::Function = mean)
  values_parts = frexp10.(values)
  s_values = map(v->v[1], values_parts); e_values = map(v->v[2], values_parts)
  sv = meanf(s_values); ev = meanf(e_values)
  return sv * 10^ev
end

function mean_base2(values::Float64...; meanf::Function = mean)
  vals = collect(values)
  return average_base2(vals; meanf = meanf)
end

function mean_base2(values::Vector{Float64}; meanf::Function = mean)
  s_values = map(significand, values); e_values = collect(Float64, map(exponent, values))
  sv = meanf(s_values); ev = meanf(e_values)
  return sv * 2^ev
end

function data_column_upper_matrix(i::Int, j::Int, N::Int)
  assert_square_matrix_indices(i, j, N)
  @assert i < j prerequisite("i < j", ("i", i), ("j", j))
  return convert(Int, (2*(N - 1) - i) * (i - 1) / 2 + (j - 1))
end

function data_column_lower_matrix(i::Int, j::Int)
  assert_square_matrix_indices(i, j, typemax(Int))
  @assert i > j prerequisite("i > j", ("i", i), ("j", j))
  return convert(Int, convert(Int, (1 + (i - 2)) * (i - 2) / 2 + j))
end

function data_size_symmetric_matrix(N::Int)
  @assert N > 0 prerequisite("N > 0", ("N", N))
  convert(Int, N*(N - 1) / 2)
end

"""
Convert vector element number to a pair of N×N matrix indices.
"""
function mpos(l::Int, N::Int)
  n, r = divrem(l, N)
  return (r == 0) ? n : n + 1, (r == 0) ? N : r
end

"""
Convert a pair of N×N matrix indices to a one-dimentional array index.
"""
function mvec(i::Int, j::Int, N::Int)
  return N*(i - 1) + j
end


text:
@testset "map_pair" begin

# test of map_pair
sc = SparkContext(master="local")

rdd = parallelize(sc, 1:3)
pairs = map_pair(rdd, x -> (x,1))
values = collect(pairs)

@test values == [(1,1), (2,1), (3,1)]

close(sc)

end

ext: jl


text:
"""
    FIRITools

`FIRITools` is a collection of convenience functions for working with
Faraday-International Reference Ionosphere (FIRI) model profiles.

The underlying FIRI-2018 model `DATA`, `HEADER`, and `ALTITUDE`'s can be accessed
as e.g. `FIRITools.HEADER`. The unit of `ALTITUDE` is meters and `DATA` is electrons per
cubic meter.

Only the function `firi` is exported, but also see `FIRITools.quantile` and
`FIRITools.extrapolate`.

# References

Friedrich, M., Pock, C., & Torkar, K. (2018).
FIRI-2018, an updated empirical model of the lower ionosphere.
Journal of Geophysical Research: Space Physics, 123, 6737-6751.
[doi: 10.1029/2018JA025437](https://doi.org/10.1029/2018JA025437)
"""
module FIRITools

using Artifacts, Statistics, DelimitedFiles
using TypedTables
using Interpolations

export firi

const MIN_DENSITY = 1e-4

const DF = readdlm(joinpath(artifact"firi", "firi2018.csv"), ',')

function parseheader()
    df = Table(
        code=convert(Vector{String}, DF[2,2:end]),
        month=convert(Vector{Int}, DF[3,2:end]),
        chi=convert(Vector{Int}, DF[5,2:end]),  # deg
        lat=convert(Vector{Int}, DF[6,2:end]),  # deg
        f10_7=convert(Vector{Int}, DF[7,2:end])
    )

    return df
end

const HEADER = parseheader()

# Start altitude at 60 km even though table has down to 55 km (60 km is stated lower limit)
const DATA = convert(Matrix{Float64}, DF[14:end-2, 2:end])  # last 2 rows have Missing
const ALTITUDE = convert(Vector{Int}, DF[14:end-2, 1])*1000  # convert to m
@assert length(ALTITUDE) == size(DATA, 1) "ALTITUDE does not match DATA"

"""
    values(s)

Return the FIRI model values of parameter `s`.
"""
function values(s)
    if s in ("f10_7", :f10_7)
        return [75, 130, 200]  # == unique(HEADER.f10_7)
    elseif s in ("chi", :chi)
        return [0, 30, 45, 60, 75, 80, 85, 90, 95, 100, 130]
    elseif s in ("lat", :lat)
        return [0, 15, 30, 45, 60]
    elseif s in ("month", :month)
        return [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]
    else
        throw(ArgumentError("$s is not an FIRI model parameter"))
    end
end

"""
    firi(;chi=(0, 130), lat=(0, 60), f10_7=(75, 200), month=(1, 12))

Return the average FIRI profile across model parameters selected by the keyword arguments.
"""
function firi(;chi=(0, 130), lat=(0, 60), f10_7=(75, 200), month=(1, 12))
    profiles = selectprofiles(chi=chi, lat=lat, f10_7=f10_7, month=month)
    profile = mean(profiles, dims=2)
    replace!(x -> x<MIN_DENSITY ? MIN_DENSITY : x, profile)
    return dropdims(profile, dims=2)
end

"""
    firi(chi, lat; f10_7=(75, 200), month=(1, 12))

Return the average FIRI profile across the model parameters selected by the keyword
arguments and interpolated at solar zenith angle `chi` and latitude `lat`.
"""
function firi(chi, lat; f10_7=(75, 200), month=(1, 12))
    profiles = selectprofiles(chi, lat, f10_7=f10_7, month=month)
    profile = mean(profiles, dims=2)
    replace!(x -> x<MIN_DENSITY ? MIN_DENSITY : x, profile)
    return dropdims(profile, dims=2)
end

"""
    selectprofiles(;chi=(0, 130), lat=(0, 60), f10_7=(75, 200), month=(1, 12))

Return matrix of selected model profiles.
"""
function selectprofiles(;chi=(0, 130), lat=(0, 60), f10_7=(75, 200), month=(1, 12))
    mask = trues(length(HEADER))

    mask .&= select(HEADER.chi, chi) .&
             select(HEADER.lat, lat) .&
             select(HEADER.f10_7, f10_7) .&
             select(HEADER.month, month)

    return DATA[:,mask]
end

"""
    selectprofiles(chi, lat; f10_7=(75, 200), month=(1, 12))

Return matrix of selected model profiles interpolated at solar zenith angle `chi` and
latitude `lat`.
"""
function selectprofiles(chi, lat; f10_7=(75, 200), month=(1, 12))
    issorted(unique(HEADER.chi)) || throw(DomainError("chi column of HEADER is not sorted"))
    issorted(unique(HEADER.lat)) || throw(DomainError("lat column of HEADER is not sorted"))

    if chi > 130
        @warn "`chi` greater than 130° uses `chi = 130°`" maxlog=3
        chi = oftype(chi, 130)
    end

    if chi < 0 || lat < 0
        throw(ArgumentError("`chi` or `lat` below 0° is not supported. See Friedrich et al., “FIRI-2018”."))
    end

    mask = trues(length(HEADER))

    mask .&= select(HEADER.f10_7, f10_7) .&
             select(HEADER.month, month)

    chibnds = twoclosest(HEADER.chi, chi)
    latbnds = twoclosest(HEADER.lat, lat)

    # Mask all but the 2 closest values of chi and lat to the target chi and lat
    mask .&= .!isnothing.(indexin(HEADER.chi, chibnds)) .&
             .!isnothing.(indexin(HEADER.lat, latbnds))

    N = count(mask)
    maskidxs = findall(mask)

    # We need to explicitly handle 1D vs 2D interpolation (if lat, or chi don't need interp)
    if isequal(chibnds...) && isequal(latbnds...)
        profiles = DATA[:, maskidxs]
    elseif isequal(chibnds...)
        profiles = Matrix{Float64}(undef, length(ALTITUDE), N÷2)

        j = 1
        for i = 1:2:N
            for a in axes(DATA, 1)  # altitude
                d = DATA[a,maskidxs[[i, i+1]]]
                itp = LinearInterpolation(latbnds, d, extrapolation_bc=Line())
                profiles[a,j] = itp(lat)
            end
            j += 1
        end
    elseif isequal(latbnds...)
        profiles = Matrix{Float64}(undef, length(ALTITUDE), N÷2)

        j = 1
        for i = 1:2:N
            for a in axes(DATA, 1)  # altitude
                d = DATA[a,maskidxs[[i, i+1]]]
                itp = LinearInterpolation(chibnds, d, extrapolation_bc=Line())
                profiles[a,j] = itp(chi)
            end
            j += 1
        end
    else
        # Although not a general solution, we know the format of Chi guarantees that occurences
        # of the two closest values will be in pairs and that both Chi and Lat are sorted
        # so every 4 indices in the mask form a group that covers
        # [(chi_low, lat_low) (chi_high, lat_low);
        #  (chi_low, lat_high) (chi_high, lat_high)]
        profiles = Matrix{Float64}(undef, length(ALTITUDE), N÷4)

        j = 1
        for i = 1:4:N
            for a in axes(DATA, 1)  # altitude
                d = [DATA[a,maskidxs[i]] DATA[a,maskidxs[i+2]];
                     DATA[a,maskidxs[i+1]] DATA[a,maskidxs[i+3]]]
                itp = LinearInterpolation((chibnds, latbnds), d, extrapolation_bc=Line())
                profiles[a,j] = itp(chi, lat)
            end
            j += 1
        end
    end

    return profiles
end

"""
    twoclosest(s, v)

Return:
    - the two closest values of vector `s` to `v` if `v` is outside of the range of
        values in `s` (extrapolating)
    - or the closest value on each side of `v` in column `s` if interpolating
    - or `[v, v]` if `v` is exactly in `s`.

The two values are sorted.
"""
function twoclosest(s, v)
    uvals = unique(s)

    if v in uvals
        return [v, v]
    elseif v >= maximum(uvals) || v <= minimum(uvals)
        # We're extrapolating or at a min/max
        I = sortperm(abs.(v .- uvals))
        return sort([uvals[I[1]], uvals[I[2]]])  # 2 closest matches
    else
        # We're interpolating, make sure we return the closest on each side of `v`
        udiffs = v .- uvals
        v1 = minimum(x->x > 0 ? x : Inf, udiffs)  # closest val with positive diff
        v2 = -minimum(x->x < 0 ? abs(x) : Inf, udiffs)  # closest val with negative diff

        I1 = findfirst(isequal(v1), udiffs)
        I2 = findfirst(isequal(v2), udiffs)
        return sort([uvals[I1], uvals[I2]])
    end
end

"""
    select(h, x::Tuple)

Return `x[1] .<= h .<= x[2]`.
"""
function select(h, x::T) where T<:Tuple
    first(x) .<= h .<= last(x)
end

"""
    select(h, x)

Return `x .== h` if `h` is in `x`, otherwise `nothing`.
"""
function select(h, x::T) where T<:Number
    x in h || @warn "$x is not in model parameters. Looking for interpolating form of `selectprofiles`?"

    return x .== h
end

"""
    quantile(p; chi=(0, 130), lat=(0, 60), f10_7=(75, 200), month=(1, 12))

Compute the quantile(s) `p` at each `ALTITUDE` with data columns selected by the keyword
arguments.
"""
function quantile(p; chi=(0, 130), lat=(0, 60), f10_7=(75, 200), month=(1, 12))
    profiles = selectprofiles(chi=chi, lat=lat, f10_7=f10_7, month=month)
    profile = _quantile(profiles, p)

    replace!(x -> x<MIN_DENSITY ? MIN_DENSITY : x, profile)

    if p isa Number
        # make `profile` a Vector for consistency with `firi`
        profile = dropdims(profile, dims=2)
    end

    return profile
end

"""
    quantile(chi, lat, p; f10_7=(75, 200), month=(1, 12))

Compute the quantile(s) `p` at each `ALTITUDE` with data columns selected by keyword
arguments and interpolated to solar zenith angle `chi` and latitude `lat`.
"""
function quantile(chi, lat, p; f10_7=(75, 200), month=(1, 12))
    profiles = selectprofiles(chi, lat, f10_7=f10_7, month=month)
    profile = _quantile(profiles, p)

    replace!(x -> x<MIN_DENSITY ? MIN_DENSITY : x, profile)

    if p isa Number
        # make `profile` a Vector for consistency with `firi`
        profile = dropdims(profile, dims=2)
    end

    return profile
end

function _quantile(profiles, p)
    profile = Array{Float64, 2}(undef, length(ALTITUDE), length(p))
    for j in eachindex(p)
        i = 1
        for r in eachrow(profiles)
            profile[i,j] = Statistics.quantile(r, p[j])
            i += 1
        end
    end
    return profile
end

"""
    extrapolate([z,] profile, newz)

Return `profile`, originally defined at sorted altitudes `z`, exponentially interpolated
and extrapolated to altitudes of `newz`.

If `z` is not specified, it is assumed that `profile` is defined at `FIRITools.ALTITUDE`.
"""
function extrapolate(z::AbstractRange, profile::AbstractVector, newz::AbstractRange)
    length(z) == length(profile) ||
        throw(ArgumentError("`z` and `profile` must be equal lengths."))

    itp = LinearInterpolation(z, log.(profile), extrapolation_bc=Line())
    itpprofile = itp(newz)

    return exp.(itpprofile)
end
extrapolate(profile, newz) = extrapolate(first(ALTITUDE):1000:last(ALTITUDE), profile, newz)

function extrapolate(z::AbstractRange, profile::AbstractMatrix, newz::AbstractRange)
    expprofile = Vector{eltype(profile)}()
    nrows = 0
    for i in axes(profile, 2)
        p = extrapolate(z, view(profile,:,i), newz)
        append!(expprofile, p)
        nrows = length(p)
    end
    return reshape(expprofile, nrows, size(profile, 2))
end

# NOTE: at this time, models() isn't very important because only 1 model is supported
"""
    models()

Return the available models.
"""
models() = first.(splitext.(filter(endswith(".csv"), readdir(artifact"firi"))))

end  # module


text:
##Functions to test create new initial, design collision and change_gauge. Run the corresponding test files.
Base.@kwdef struct analytic_parameters
    energy   = 1.0
    phi0     = 1.0
    u_AH     = 1.0
    a4_amp   = 0.01
    phi2_amp = 0.01
    b14_amp  = 0.01
    b24_amp  = 0.01
end

a4(x,y,a40,a4_amp,xmid,ymid,Lx,Ly) = a40*(1 + a4_amp*cos(2*pi*(x-xmid)/Lx)*cos(2*pi*(y-ymid)/Ly))
phi2(x,y,phi2_amp,xmid,ymid,Lx,Ly) = phi2_amp*cos(2*pi*(x-xmid)/Lx)*cos(2*pi*(y-ymid)/Ly)
b14(x,y,b14_amp,xmid,ymid,Lx,Ly)   = b14_amp*cos(2*pi*(x-xmid)/Lx)*cos(2*pi*(y-ymid)/Ly)
b24(x,y,b24_amp,xmid,ymid,Lx,Ly)   = b24_amp*cos(2*pi*(x-xmid)/Lx)*cos(2*pi*(y-ymid)/Ly)
B1(u,x,y,b14_amp,xmid,ymid,Lx,Ly)  = u^4*b14(x,y,b14_amp,xmid,ymid,Lx,Ly)
B2(u,x,y,b24_amp,xmid,ymid,Lx,Ly)  = u^4*b24(x,y,b24_amp,xmid,ymid,Lx,Ly)
phi(u,x,y,phi0,xi0,phi2_amp,xmid,ymid,Lx,Ly) = phi0*u-xi0*phi0*u^2+u^3*(phi2(x,y,phi2_amp,xmid,ymid,Lx,Ly)+xi0^2*phi0)

function test_output_analytic(grid::SpecCartGrid3D,prm::analytic_parameters, io::InOut, potential::Potential)
    dirname   = io.out_dir
    energy    = prm.energy
    phi0      = prm.phi0
    u_AH      = prm.u_AH
    a4_amp    = prm.a4_amp
    phi2_amp  = prm.phi2_amp
    b14_amp   = prm.b14_amp
    b24_amp   = prm.b24_amp

    a40 = -4/3*energy
    xi0 = (-a40)^0.25-1/u_AH

    Nx   = grid.x_nodes
    Ny   = grid.y_nodes
    dx   = (grid.x_max-grid.x_min)/Nx
    dy   = (grid.y_max-grid.y_min)/Ny
    Lx   = grid.x_max-grid.x_min
    Ly   = grid.y_max-grid.y_min
    xmid = (grid.x_max+grid.x_min)/2
    ymid = (grid.y_max+grid.y_min)/2

    aa4(x,y)    = a4(x,y,a40,a4_amp,xmid,ymid,Lx,Ly)
    bb14(x,y)   = b14(x,y,b14_amp,xmid,ymid,Lx,Ly)
    bb24(x,y)   = b24(x,y,b24_amp,xmid,ymid,Lx,Ly)
    pphi2(x,y)  = phi2(x,y,phi2_amp,xmid,ymid,Lx,Ly)
    BB1(u,x,y)  = B1(u,x,y,b14_amp,xmid,ymid,Lx,Ly)
    BB2(u,x,y)  = B2(u,x,y,b24_amp,xmid,ymid,Lx,Ly)
    pphi(u,x,y) = phi(u,x,y,phi0,xi0,phi2_amp,xmid,ymid,Lx,Ly)

    atlas     = Atlas(grid)
    systems   = SystemPartition(grid)
    Nsys      = length(systems)
    boundary  = Boundary(grid)
    gauge     = Gauge(grid)
    bulkevols = BulkEvolvedPartition(grid)
    tinfo     = Jecco.TimeInfo(0, 0.0, 0.0, 0.0)
    empty     = Cartesian{1}("u", 0.0, 0.0, 1)
    x         = systems[1].xcoord[:]
    y         = systems[1].ycoord[:]
    chart2D   = Chart(empty, systems[1].xcoord, systems[1].ycoord)

    fill!(boundary.fx2,0.0)
    fill!(boundary.fy2,0.0)
    fill!(gauge.xi,xi0)

    for k in 1:Ny
        for j in 1:Nx
            boundary.a4[1,j,k] = aa4(x[j],y[k])
        end
    end

    for l in 1:Nsys
        u  = systems[l].ucoord[:]
        Nu = length(u)
        for k in 1:Ny
            for j in 1:Nx
                if l == 1
                    bulkevols[1].B1[:,j,k] .= bb14(x[j],y[k])
                    bulkevols[1].B2[:,j,k] .= bb24(x[j],y[k])
                    bulkevols[1].G[:,j,k]  .= 0.0
                    bulkevols[1].phi[:,j,k].= pphi2(x[j],y[k])/phi0^3+xi0^2/phi0^2
                else
                    for i in 1:Nu
                        bulkevols[l].B1[i,j,k]  = BB1(u[i],x[j],y[k])
                        bulkevols[l].B2[i,j,k]  = BB2(u[i],x[j],y[k])
                        bulkevols[l].G[i,j,k]   = 0.0
                        bulkevols[l].phi[i,j,k] = pphi(u[i],x[j],y[k])
                    end
                end
            end
        end
    end

    evolvars = AdS5_3_1.EvolVars(boundary, gauge, bulkevols)
    AdS5_3_1.create_outputs(dirname, evolvars, chart2D, atlas.charts, io, potential, phi0)
end

function test_create_new_data(grid::SpecCartGrid3D,new_grid::SpecCartGrid3D,prm::analytic_parameters,new_parameters::AdS5_3_1.NewParameters,
                                 io::InOut, potential::Potential)

    test_output_analytic(grid,prm,io,potential)
    #We modify the configuration now

    energy    = prm.energy
    phi0      = prm.phi0
    u_AH      = prm.u_AH
    a4_amp    = prm.a4_amp
    phi2_amp  = prm.phi2_amp
    b14_amp   = prm.b14_amp
    b24_amp   = prm.b24_amp

    a40 = -4/3*energy
    xi0 = (-a40)^0.25-1/u_AH

    Nx   = grid.x_nodes
    Ny   = grid.y_nodes
    dx   = (grid.x_max-grid.x_min)/Nx
    dy   = (grid.y_max-grid.y_min)/Ny
    Lx   = grid.x_max-grid.x_min
    Ly   = grid.y_max-grid.y_min
    xmid = (grid.x_max+grid.x_min)/2
    ymid = (grid.y_max+grid.y_min)/2

    aa4(x,y)    = a4(x,y,a40,a4_amp,xmid,ymid,Lx,Ly)
    bb14(x,y)   = b14(x,y,b14_amp,xmid,ymid,Lx,Ly)
    bb24(x,y)   = b24(x,y,b24_amp,xmid,ymid,Lx,Ly)
    pphi2(x,y)  = phi2(x,y,phi2_amp,xmid,ymid,Lx,Ly)
    BB1(u,x,y)  = B1(u,x,y,b14_amp,xmid,ymid,Lx,Ly)
    BB2(u,x,y)  = B2(u,x,y,b24_amp,xmid,ymid,Lx,Ly)
    pphi(u,x,y) = phi(u,x,y,phi0,xi0,phi2_amp,xmid,ymid,Lx,Ly)

    atlas         = Atlas(new_grid)
    systems       = SystemPartition(new_grid)
    Nsys          = length(systems)
    boundary_new  = Boundary(new_grid)
    gauge_new     = Gauge(new_grid)
    bulkevols_new = BulkEvolvedPartition(new_grid)
    empty         = Cartesian{1}("u", 0.0, 0.0, 1)
    x             = systems[1].xcoord[:]
    y             = systems[1].ycoord[:]
    chart2D       = Chart(empty, systems[1].xcoord, systems[1].ycoord)

    Nx = length(x)
    Ny = length(y)

    fill!(boundary_new.fx2,0.0)
    fill!(boundary_new.fy2,0.0)
    fill!(gauge_new.xi,xi0)

    phi2_new = zeros(Nx,Ny)

    for i in 1:Nsys
        u  = systems[i].ucoord[:]
        Nu = length(u)
        for l in 1:Ny
            for k in 1:Nx
                if grid.x_min<=x[k]<=grid.x_max && grid.y_min<=y[l]<=grid.y_max
                    if i == 1
                        boundary_new.a4[1,k,l]      = aa4(x[k],y[l])
                        phi2_new[k,l]               = pphi2(x[k],y[l])
                        bulkevols_new[1].B1[:,k,l] .= bb14(x[k],y[l])
                        bulkevols_new[1].B2[:,k,l] .= bb24(x[k],y[l])
                        bulkevols_new[1].G[:,k,l]  .= 0.0
                        bulkevols_new[1].phi[:,k,l].= pphi2(x[k],y[l])/phi0^3+xi0^2/phi0^2
                    else
                        for j in 1:Nu
                            bulkevols_new[i].B1[j,k,l]  = BB1(u[j],x[k],y[l])
                            bulkevols_new[i].B2[j,k,l]  = BB2(u[j],x[k],y[l])
                            bulkevols_new[i].G[j,k,l]   = 0.0
                            bulkevols_new[i].phi[j,k,l] = pphi(u[j],x[k],y[l])
                        end
                    end
                elseif grid.x_min<=x[k]<=grid.x_max
                    if i == 1
                        boundary_new.a4[1,k,l]      = aa4(x[k],grid.y_min)
                        phi2_new[k,l]               = pphi2(x[k],grid.y_min)
                        bulkevols_new[1].B1[:,k,l] .= bb14(x[k],grid.y_min)
                        bulkevols_new[1].B2[:,k,l] .= bb24(x[k],grid.y_min)
                        bulkevols_new[1].G[:,k,l]  .= 0.0
                        bulkevols_new[1].phi[:,k,l].= pphi2(x[k],grid.y_min)/phi0^3+xi0^2/phi0^2
                    else
                        for j in 1:Nu
                            bulkevols_new[i].B1[j,k,l]  = BB1(u[j],x[k],grid.y_min)
                            bulkevols_new[i].B2[j,k,l]  = BB2(u[j],x[k],grid.y_min)
                            bulkevols_new[i].G[j,k,l]   = 0.0
                            bulkevols_new[i].phi[j,k,l] = pphi(u[j],x[k],grid.y_min)
                        end
                    end
                elseif grid.y_min<=y[l]<=grid.y_max
                    if i == 1
                        boundary_new.a4[1,k,l]      = aa4(grid.x_min,y[l])
                        phi2_new[k,l]               = pphi2(grid.x_min,y[l])
                        bulkevols_new[1].B1[:,k,l] .= bb14(grid.x_min,y[l])
                        bulkevols_new[1].B2[:,k,l] .= bb24(grid.x_min,y[l])
                        bulkevols_new[1].G[:,k,l]  .= 0.0
                        bulkevols_new[1].phi[:,k,l].= pphi2(grid.x_min,y[l])/phi0^3+xi0^2/phi0^2
                    else
                        for j in 1:Nu
                            bulkevols_new[i].B1[j,k,l]  = BB1(u[j],grid.x_min,y[l])
                            bulkevols_new[i].B2[j,k,l]  = BB2(u[j],grid.x_min,y[l])
                            bulkevols_new[i].G[j,k,l]   = 0.0
                            bulkevols_new[i].phi[j,k,l] = pphi(u[j],grid.x_min,y[l])
                        end
                    end
                else
                    if i == 1
                        boundary_new.a4[1,k,l]      = aa4(grid.x_min,grid.y_min)
                        phi2_new[k,l]               = pphi2(grid.x_min,grid.y_min)
                        bulkevols_new[1].B1[:,k,l] .= bb14(grid.x_min,grid.y_min)
                        bulkevols_new[1].B2[:,k,l] .= bb24(grid.x_min,grid.y_min)
                        bulkevols_new[1].G[:,k,l]  .= 0.0
                        bulkevols_new[1].phi[:,k,l].= pphi2(grid.x_min,grid.y_min)/phi0^3+xi0^2/phi0^2
                    else
                        for j in 1:Nu
                            bulkevols_new[i].B1[j,k,l]  = BB1(u[j],grid.x_min,grid.y_min)
                            bulkevols_new[i].B2[j,k,l]  = BB2(u[j],grid.x_min,grid.y_min)
                            bulkevols_new[i].G[j,k,l]   = 0.0
                            bulkevols_new[i].phi[j,k,l] = pphi(u[j],grid.x_min,grid.y_min)
                        end
                    end
                end
            end
        end
    end

    e       = AdS5_3_1.compute_energy.(boundary_new.a4[1,:,:], phi2_new, phi0, potential.oophiM2)
    e_new   = new_parameters.e_new
    fx20    = new_parameters.fx20
    fy20    = new_parameters.fy20
    #change of energy
    if e_new != -1.0
        plan    = plan_rfft(e)
        e0      = real(1/(Nx*Ny)*(plan*e)[1])
        a40_new = real(1/(Nx*Ny)*(plan*boundary_new.a4[1,:,:])[1])
        k       = (-4/3*(e_new-e0)*phi0^4-maximum(boundary_new.a4[1,:,:])+a40_new)/(a40_new-maximum(boundary_new.a4[1,:,:]))
        boundary_new.a4[1,:,:] = k*(boundary_new.a4[1,:,:].-maximum(boundary_new.a4[1,:,:])).+maximum(boundary_new.a4[1,:,:])
        #a4_new  = a4_new.+4/3*(e0-e_new)
    end
    #boosting, with static cold phase (edge).
    if new_parameters.boostx
        boundary_new.fx2[1,:,:] .= fx20.*(e.-minimum(e))
    end
    if new_parameters.boosty
        boundary_new.fy2[1,:,:] .= fy20.*(e.-minimum(e))
    end
    #We add the specified a4 cos perturbation
    ampx    = new_parameters.a4_ampx
    ampy    = new_parameters.a4_ampy
    kx      = new_parameters.a4_kx
    ky      = new_parameters.a4_ky
    xmin    = new_grid.x_min
    xmax    = new_grid.x_max
    ymin    = new_grid.y_min
    ymax    = new_grid.y_max
    xmid    = (xmax+xmin)/2
    ymid    = (ymax+ymin)/2

    plan    = plan_rfft(boundary_new.a4[1,:,:])
    a40_new = real(1/(Nx*Ny)*(plan*boundary_new.a4[1,:,:])[1])

    for j in 1:Ny
        for i in 1:Nx
            boundary_new.a4[1,i,j] -= a40_new*(ampx*cos(2*π*kx*(x[i]-xmid)/(xmax-xmin))
                                                  +ampy*cos(2*π*ky*(y[j]-ymid)/(ymax-ymin)))
        end
    end

    return boundary_new, gauge_new, bulkevols_new
end

function test_design_collision(grid::SpecCartGrid3D, new_grid::SpecCartGrid3D, prm::analytic_parameters,
                                           new_parameters_coll::AdS5_3_1.NewParameters, io::InOut, potential::Potential)

    test_output_analytic(grid,prm,io,potential)

    energy    = prm.energy
    phi0      = prm.phi0
    u_AH      = prm.u_AH
    a4_amp    = prm.a4_amp
    phi2_amp  = prm.phi2_amp
    b14_amp   = prm.b14_amp
    b24_amp   = prm.b24_amp

    a40 = -4/3*energy
    xi0 = (-a40)^0.25-1/u_AH

    Nx   = grid.x_nodes
    Ny   = grid.y_nodes
    dx1  = (grid.x_max-grid.x_min)/Nx
    dy1  = (grid.y_max-grid.y_min)/Ny
    dx2  = dx1
    dy2  = dy1
    Lx   = grid.x_max-grid.x_min
    Ly   = grid.y_max-grid.y_min
    xmid = (grid.x_max+grid.x_min)/2
    ymid = (grid.y_max+grid.y_min)/2
    x1mid = xmid
    y1mid = ymid
    x2mid = xmid
    y2mid = ymid

    aa4(x,y)    = a4(x,y,a40,a4_amp,xmid,ymid,Lx,Ly)
    bb14(x,y)   = b14(x,y,b14_amp,xmid,ymid,Lx,Ly)
    bb24(x,y)   = b24(x,y,b24_amp,xmid,ymid,Lx,Ly)
    pphi2(x,y)  = phi2(x,y,phi2_amp,xmid,ymid,Lx,Ly)
    BB1(u,x,y)  = B1(u,x,y,b14_amp,xmid,ymid,Lx,Ly)
    BB2(u,x,y)  = B2(u,x,y,b24_amp,xmid,ymid,Lx,Ly)
    pphi(u,x,y) = phi(u,x,y,phi0,xi0,phi2_amp,xmid,ymid,Lx,Ly)

    systems = SystemPartition(grid)
    x1      = systems[1].xcoord[:]
    y1      = systems[1].ycoord[:]
    x2      = x1
    y2      = y1

    atlas         = Atlas(new_grid)
    systems       = SystemPartition(new_grid)
    Nsys          = length(systems)
    boundary_new  = Boundary(new_grid)
    gauge_new     = Gauge(new_grid)
    bulkevols_new = BulkEvolvedPartition(new_grid)
    empty         = Cartesian{1}("u", 0.0, 0.0, 1)
    x_new         = systems[1].xcoord[:]
    y_new         = systems[1].ycoord[:]
    chart2D       = Chart(empty, systems[1].xcoord, systems[1].ycoord)

    Nx = length(x_new)
    Ny = length(y_new)

#Position of the center of the boxes inside the new box
    x1mid_new = new_parameters_coll.x1_center
    y1mid_new = new_parameters_coll.y1_center
    x2mid_new = new_parameters_coll.x2_center
    y2mid_new = new_parameters_coll.y2_center
#New bounds of the old boxwes inside the new box
    x1max_new = x1[end]-x1mid+x1mid_new
    y1max_new = y1[end]-y1mid+y1mid_new
    x2max_new = x2[end]-x2mid+x2mid_new
    y2max_new = y2[end]-y2mid+y2mid_new
    x1min_new = x1[1]-x1mid+x1mid_new
    y1min_new = y1[1]-y1mid+y1mid_new
    x2min_new = x2[1]-x2mid+x2mid_new
    y2min_new = y2[1]-y2mid+y2mid_new
    #Check that we are within bounds.
    if x1max_new+dx1>new_grid.x_max || x1min_new<x_new[1] || y1max_new+dy1>new_grid.y_max || y1min_new<y_new[1]
        println("Box 1 out of bounds.")
        return nothing
    elseif x2max_new+dx2>new_grid.x_max || x2min_new<x_new[1] || y2max_new+dy2>new_grid.y_max || y2min_new<y_new[1]
        println("Box 2 out of bounds.")
        return nothing
    end

    e = zeros(length(x1),length(y1))
    for j in 1:length(y1)
        for i in 1:length(x1)
            e[i,j] = AdS5_3_1.compute_energy(aa4(x1[i],y1[j]),pphi2(x1[i],y1[j]),phi0,potential.oophiM2)
        end
    end
    e_min = minimum(e)

    #Read the desired new momenta
    fx21 = new_parameters_coll.fx21
    fy21 = new_parameters_coll.fy21
    fx22 = new_parameters_coll.fx22
    fy22 = new_parameters_coll.fy22
    #Embedding in the new box
    phi2_new = zeros(Nx,Ny)
    for i in 1:Nsys
        u  = systems[i].ucoord[:]
        Nu = length(u)
        for l in 1:Ny
            for k in 1:Nx
                if x1min_new<=x_new[k]<=x1max_new && y1min_new<=y_new[l]<=y1max_new #we displace box 1
                    if i == 1
                        boundary_new.a4[1,k,l]      = aa4(x_new[k]-x1mid_new+x1mid,y_new[l]-y1mid_new+y1mid)
                        phi2_new[k,l]               = pphi2(x_new[k]-x1mid_new+x1mid,y_new[l]-y1mid_new+y1mid)
                        e_aux                       = AdS5_3_1.compute_energy(boundary_new.a4[1,k,l],phi2_new[k,l],phi0,potential.oophiM2)-e_min
                        boundary_new.fx2[1,k,l]     = fx21*e_aux
                        boundary_new.fy2[1,k,l]     = fy21*e_aux
                        bulkevols_new[1].B1[:,k,l] .= bb14(x_new[k]-x1mid_new+x1mid,y_new[l]-y1mid_new+y1mid)
                        bulkevols_new[1].B2[:,k,l] .= bb24(x_new[k]-x1mid_new+x1mid,y_new[l]-y1mid_new+y1mid)
                        bulkevols_new[1].G[:,k,l]  .= 0.0
                        bulkevols_new[1].phi[:,k,l].= pphi2(x_new[k]-x1mid_new+x1mid,y_new[l]-y1mid_new+y1mid)/phi0^3+xi0^2/phi0^2
                    else
                        for j in 1:Nu
                            bulkevols_new[i].B1[j,k,l]  = BB1(u[j],x_new[k]-x1mid_new+x1mid,y_new[l]-y1mid_new+y1mid)
                            bulkevols_new[i].B2[j,k,l]  = BB2(u[j],x_new[k]-x1mid_new+x1mid,y_new[l]-y1mid_new+y1mid)
                            bulkevols_new[i].G[j,k,l]   = 0.0
                            bulkevols_new[i].phi[j,k,l] = pphi(u[j],x_new[k]-x1mid_new+x1mid,y_new[l]-y1mid_new+y1mid)
                        end
                    end
                elseif x2min_new<=x_new[k]<=x2max_new && y2min_new<=y_new[l]<=y2max_new #we displace box 2
                    if i == 1
                        boundary_new.a4[1,k,l]      = aa4(x_new[k]-x2mid_new+x2mid,y_new[l]-y2mid_new+y2mid)
                        phi2_new[k,l]               = pphi2(x_new[k]-x2mid_new+x2mid,y_new[l]-y2mid_new+y2mid)
                        e_aux                       = AdS5_3_1.compute_energy(boundary_new.a4[1,k,l],phi2_new[k,l],phi0,potential.oophiM2)-e_min
                        boundary_new.fx2[1,k,l]     = fx22*e_aux
                        boundary_new.fy2[1,k,l]     = fy22*e_aux
                        bulkevols_new[1].B1[:,k,l] .= bb14(x_new[k]-x2mid_new+x2mid,y_new[l]-y2mid_new+y2mid)
                        bulkevols_new[1].B2[:,k,l] .= bb24(x_new[k]-x2mid_new+x2mid,y_new[l]-y2mid_new+y2mid)
                        bulkevols_new[1].G[:,k,l]  .= 0.0
                        bulkevols_new[1].phi[:,k,l].= pphi2(x_new[k]-x2mid_new+x2mid,y_new[l]-y2mid_new+y2mid)/phi0^3+xi0^2/phi0^2
                    else
                        for j in 1:Nu
                            bulkevols_new[i].B1[j,k,l]  = BB1(u[j],x_new[k]-x2mid_new+x2mid,y_new[l]-y2mid_new+y2mid)
                            bulkevols_new[i].B2[j,k,l]  = BB2(u[j],x_new[k]-x2mid_new+x2mid,y_new[l]-y2mid_new+y2mid)
                            bulkevols_new[i].G[j,k,l]   = 0.0
                            bulkevols_new[i].phi[j,k,l] = pphi(u[j],x_new[k]-x2mid_new+x2mid,y_new[l]-y2mid_new+y2mid)
                        end
                    end
                else
                    if i == 1
                        boundary_new.a4[1,k,l]      = aa4(grid.x_min,grid.y_min)
                        phi2_new[k,l]               = pphi2(grid.x_min,grid.y_min)
                        boundary_new.fx2[1,k,l]     = 0.0
                        boundary_new.fy2[1,k,l]     = 0.0
                        bulkevols_new[1].B1[:,k,l] .= bb14(grid.x_min,grid.y_min)
                        bulkevols_new[1].B2[:,k,l] .= bb24(grid.x_min,grid.y_min)
                        bulkevols_new[1].G[:,k,l]  .= 0.0
                        bulkevols_new[1].phi[:,k,l].= pphi2(grid.x_min,grid.y_min)/phi0^3+xi0^2/phi0^2
                    else
                        for j in 1:Nu
                            bulkevols_new[i].B1[j,k,l]  = BB1(u[j],grid.x_min,grid.y_min)
                            bulkevols_new[i].B2[j,k,l]  = BB2(u[j],grid.x_min,grid.y_min)
                            bulkevols_new[i].G[j,k,l]   = 0.0
                            bulkevols_new[i].phi[j,k,l] = pphi(u[j],grid.x_min,grid.y_min)
                        end
                    end
                end
            end
        end
    end

    fill!(gauge_new.xi,xi0)

    return boundary_new, gauge_new, bulkevols_new
end


text:
using LightXML

export parseGPX, readGPX

# Converts a file to XML Document.
function readGPX(filename::String)
    return parse_file(filename)
end

#= Converts a XMLDocument to array of longitude and latitude.

Args:
    gpxDoc(XMLDocument): XML document in GPX format.
Returns
    Array{Float64, 2}: GPS coordinates with longitude and latitude.
=#
function parseGPX(gpxDoc::XMLDocument)

    gpxPath = Array{Float64, 2}
    gpxLng = zeros(0)
    gpxLon = zeros(0)

    xroot = root(gpxDoc)
    trk = get_elements_by_tagname(xroot, "trk")[1]
    trkseg = get_elements_by_tagname(trk, "trkseg")[1]
    trkpt = get_elements_by_tagname(trkseg, "trkpt")

    for c in child_nodes(trkseg)
        if is_elementnode(c)
            ad = attributes_dict(XMLElement(c))

            append!(gpxLng, parse(Float64, (ad["lat"])))
            append!(gpxLon, parse(Float64, (ad["lon"])))
        end
    end

    return cat(dims=2, gpxLng, gpxLon)
end


text:
module JuliaPackageExample

greet() = print("Hello World!")

include("extra_file.jl")
export my_f, my_g

end # module


text:
include("src/SSPS.jl")

using .SSPS

max_stars_repo_name: gitter-lab/ssps

max_issues_repo_name: gitter-lab/ssps

max_forks_repo_name: gitter-lab/ssps


text:
# ---
# title: 1400. Construct K Palindrome Strings
# id: problem1400
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Greedy
# link: <https://leetcode.com/problems/construct-k-palindrome-strings/description/>
# hidden: true
# ---
#
# Given a string `s` and an integer `k`. You should construct `k` non-empty
# **palindrome** strings using **all the characters** in `s`.
#
# Return _**True**_ if you can use all the characters in `s` to construct `k`
# palindrome strings or _**False**_ otherwise.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "annabelle", k = 2
#     Output: true
#     Explanation: You can construct two palindromes using all characters in s.
#     Some possible constructions "anna" + "elble", "anbna" + "elle", "anellena" + "b"
#
#
# **Example 2:**
#
#
#
#     Input: s = "leetcode", k = 3
#     Output: false
#     Explanation: It is impossible to construct 3 palindromes using all the characters of s.
#
#
# **Example 3:**
#
#
#
#     Input: s = "true", k = 4
#     Output: true
#     Explanation: The only possible solution is to put each character in a separate string.
#
#
# **Example 4:**
#
#
#
#     Input: s = "yzyzyzyzyzyzyzy", k = 2
#     Output: true
#     Explanation: Simply you can put all z's in one string and all y's in the other string. Both strings will be palindrome.
#
#
# **Example 5:**
#
#
#
#     Input: s = "cr", k = 7
#     Output: false
#     Explanation: We don't have enough characters in s to construct 7 palindromes.
#
#
#
#
# **Constraints:**
#
#   * `1 <= s.length <= 10^5`
#   * All characters in `s` are lower-case English letters.
#   * `1 <= k <= 10^5`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
################################################################################
#
# RGradient comparison for the BG Compass model
#
################################################################################
### Module loading
using LinearAlgebra, JLD2, PyPlot, Images
using JUDI.TimeModeling
push!(LOAD_PATH, string(pwd(), "/src/"));
using TWRIdual

### Load synthetic data
@load "./data/GaussLens_data_tti.jld"

# Filter data and source
freq = 0.006f0
df = .001f0
dat2 = low_filter(dat, dat.geometry.dt[1]; fmin=1f3*(freq-df), fmax=1f3*(freq+df))
fsrc2 = low_filter(fsrc, fsrc.geometry.dt[1]; fmin=1f3*(freq-df), fmax=1f3*(freq+df))

### Some model parameters
idx_w = 0
var = 20
n = model_true.n
d = model_true.d
o = model_true.o

function compute_gradients(model0, fsrc, dat)
    m0 = model0.m
    # Time sampling
    dt_comp = get_dt(model0)
    nt_comp = get_computational_nt(fsrc.geometry, dat.geometry, model0)

    # Pre- and post-conditioning
    mask = BitArray(undef, model0.n)
    mask .= false
    mask[2:end-1, idx_w+1:end-1] .= true
    preproc(x) = contr2abs(x, mask, m0)
    postproc(g) = gradprec_contr2abs(g, mask, m0)

    # [FWI]
    inv_name = "FWI"
    fun!(F, G, x) = objFWI!(F, G, preproc(x), model0, fsrc, dat; gradmprec_fun=postproc)

    # [TWRIdual]
    inv_name = "TWRIdual"

    ε0 = 0.01f0
    ε = Array{Float32, 1}(undef, fsrc.nsrc)
    ε .= ε0
    grad_corr = false
    objfact = 1f0
    v_bg = sqrt(1/m0[1])
    freq_peak = 0.003f0
    δ = 1f0*R(sqrt(2)/2)*v_bg/freq_peak
    weight_fun_pars = ("srcfocus", δ)
    fun!(F, G, x, weight_fun_pars, objfact) = objTWRIdual!(F, G, preproc(x), model0, fsrc, dat, ε;
                                                        objfact = objfact, comp_alpha = true, gradmprec_fun=postproc,
                                                        grad_corr = grad_corr, weight_fun_pars = weight_fun_pars)

    fun_wri!(F, G, x) = fun!(F, G, x, weight_fun_pars, objfact)
    ### Computing gradients
    x0 = zeros(R, length(findall(mask .== true)))

    G_FWI = zeros(R, n)
    G = zeros(R, size(x0))
    fun!(true, G, x0)
    G_FWI[mask] .= G
    G_FWI = G_FWI/norm(G_FWI, Inf)

    ### Computing gradients
    x02 = zeros(R, length(findall(mask .== true)))

    G_WRI = zeros(R, n)
    G2 = zeros(R, size(x02))
    fun_wri!(true, G2, x02)
    G_WRI[mask] .= G2
    G_WRI = G_WRI/norm(G_WRI, Inf)

    return G_FWI, G_WRI
end

# True Thomsen parameters model
model0_tti = deepcopy(model_true)
model0_tti.m .= .25

g_fwi_tti, g_wri_tti = compute_gradients(model0_tti, fsrc2, dat2)
@save "./data/gl_tti_g.jld" g_fwi_tti g_wri_tti model0_tti model_true

# Erronous Thomsen parameters model
model0_tti_err = deepcopy(model_true)

model0_tti_err.m .= .25
model0_tti_err.epsilon[:, idx_w+1:end] = R.(imfilter(model0_tti_err.epsilon[:, idx_w+1:end], Kernel.gaussian(var)))
model0_tti_err.delta[:, idx_w+1:end] = R.(imfilter(model0_tti_err.delta[:, idx_w+1:end], Kernel.gaussian(var)))
model0_tti_err.theta[:, idx_w+1:end] = R.(imfilter(model0_tti_err.theta[:, idx_w+1:end], Kernel.gaussian(var)))

g_fwi_tti_w, g_wri_tti_w = compute_gradients(model0_tti_err, fsrc2, dat2)
@save "./data/gl_tti_g_w.jld" g_fwi_tti_w g_wri_tti_w model0_tti_err model_true

# Acoustic model
model0_acou = Model(n, d, o, model0_tti.m)

g_fwi_a, g_wri_a = compute_gradients(model0_acou, fsrc, dat)
@save "./data/gl_acou_g.jld" g_fwi_a g_wri_a model0_acou model_true


text:
#!/usr/local/bin/julia

# NOTE: This doesn't work. Window GUI components don't get shown
using RocketBuilder
using Gtk
using SQLite

db = opendb()

tanklist = GtkListStore(String)

for name in tanknames(db)
   push!(tanklist, (name,))
end

win = maketankeditor(tanklist, db)

if !isinteractive()
    cond = Condition()
    signal_connect(win, :destroy) do widget
        notify(cond)
        DBInterface.close!(db)
    end
    @async Gtk.gtk_main()
    showall(win)
    wait(cond)
end


text:
# Plot deviation by (gpa, yp, quality).
# One plot per quality. One subplot per parental.
# function plot_gpa_yp_qual(ds :: DataSettings,
# 	dev :: ModelParams.Deviation,  plotModel :: Bool,  filePath :: String)

# 	xStrV = gpa_labels(ds);
# 	yStrV = parental_labels(ds);
#     dataM = get_data_values(dev);
#     modelM = get_model_values(dev; matchData = true);
#     nc = size(dataM, 3)

#     for ic = 1 : nc
#         qualPath = make_fpath_by_quality(filePath, ic);
#         plot_model_data_xy(modelM[:,:,ic], dataM[:,:,ic],
#             xStrV,  yStrV,  plotModel,  qualPath);
#     end
# 	return nothing
# end

# function make_fpath_by_quality(fPathIn, ic)
#     fPath, fExt = splitext(fPathIn);
#     return fPath * "_qual$ic" * fExt
# end

## ---------------  Individual moments

# Fraction in each quality (not conditional on entry), by GPA, parental
# Matches `StatsByXYQ.fracEnter_xyqM`.
# function dm_qual_entry_gpa_yp()
#     target = :fracEnter_gpqM;
#     return DataMoment(target, model_stat(target),
#         file_name("fracEnterQual", [:gpa, :parental], ".dat"), nothing,
#         qual_entry_gpa_yp, plot_gpa_yp_qual)
# end

## Fraction by quality, given [gpa, yp]
# Not conditional on entry. As a matrix by [gpa, yp, quality] that matches
# `StatsByXYQ.fracEnter_xyqM`
# function qual_entry_gpa_yp(ds :: DataSettings)
# 	target = :fracEnter_gpqM;
# 	m = load_qual_entry_gpa_yp_all(ds; conditionalOnEntry = false);
# 	n = length(m);

# 	return Deviation{Double}(name = target, dataV = m, modelV = m,
# 		scalarWt = 8.0 ./ n,
# 		shortStr = string(target),
# 		longStr = "Frac quality by [gpa, parental], not conditional",
# 		showPath = "fracEnterQualByGpaYp.txt")
# end

# Load fraction by quality | [gpa, yp]
# Not conditional on entry. As a matrix by [gpa, yp, quality] that matches
# `StatsByXYQ.fracEnter_xyqM`
# probably best to call those directly from CollegeStrat?
# check std errors ++++++
function load_qual_entry_gpa_yp_all(ds :: DataSettings;
	conditionalOnEntry :: Bool = false)

	# Counts should depend on conditionalOnEntry +++
	_, _, cnt_xyM = load_entry_gpa_yp(ds);

	nc = n_colleges(ds);
	fracEnter_xyqM = zeros(Double, n_gpa(ds), n_parental(ds), nc);
	cnt_xyqM = zeros(Int, n_gpa(ds), n_parental(ds), nc);

	for ic = 1 : nc
		fracEnter_xyqM[:,:,ic] = load_qual_entry_gpa_yp(ds, ic;
			conditionalOnEntry = conditionalOnEntry, momentType = :mean);
		cnt_xyqM[:,:,ic] = cnt_xyM;
	end
	if conditionalOnEntry
		# Should sum to 1 across qualities
		qSumM = sum(fracEnter_xyqM, dims = 3);
		@assert check_float_array(qSumM, 0.998, 1.002)
	end

	ses_xyqM = ses_from_choice_probs(fracEnter_xyqM, cnt_xyqM);

	return fracEnter_xyqM, ses_xyqM, cnt_xyqM
end

# Load fraction by quality | [gpa, yp]. Conditional on entry or not.
function load_qual_entry_gpa_yp(ds :: DataSettings, iCollege;
	conditionalOnEntry :: Bool = false, momentType :: Symbol = :mean)
	fPath = data_file(
		raw_qual_entry_gpa_parental(ds, iCollege; momentType = momentType));
	m = read_by_gpa_yp(ds, fPath);
	if momentType == :mean
		@assert all(m .< 1.0)  &&  all(m .>= 0.0)
	end
	@assert size(m) == (n_gpa(ds), n_parental(ds))

	if !conditionalOnEntry
		entryM, _ = load_moment(ds, :fracEnter_gpM);
		m .*= entryM;
	end
	return m
end

# ----------------


text:
####
# Moste complex code , can probebly be improved
####

abstract type AbstractDesign{T} end

""""
    Single_Design{T} <: AbstractDesign{T}

A Single design is a structure that contains one special define type of construct.
They can be combined to form more complex design space.
"""
abstract type Single_Design{T} <: AbstractDesign{T} end

"""
    Ordered_Design{Tm <: AbstractMod,Tcon <: AbstractConstrains,Tspace <: AbstractSpace{T} where T} <: Single_Design{Tspace}

Ordered space is where the possition in the construct has a functional role, the construct [a,b,c] is different compeared to [c,b,a].
In the General case every module can be used infinted times.

# Arguments
- `mod::AbstractMod`: The used moduels in this design
- `len::Int` : The allowed length of the construct
- `con::AbstractConstrains`:: The constrains of the design space
- `space::AbstractSpace`:: Design space containing all constructs, if possible, these constructs are not explicitly generated.
"""
struct Ordered_Design{Tm <: AbstractMod,Tcon <: AbstractConstrains,Tspace <: AbstractSpace{T} where T} <: Single_Design{Tspace}
    mod::Tm
    len::Int
    con::Tcon
    space:: Tspace
end

"""
    make_ordered_space(mod::Group_Mod,len::Int)

The function generates all ordered constructs with given length and modules. The modules can be infinity redrawn.
The constructs aren't explicitly generated.
With the use of the Kronecker product, every construct can be made on the fly only on the moment it is required.
"""
make_ordered_space(mod::Group_Mod,len::Int) =  reshape(mod.m,length(mod.m),1) |> (y -> len > 1 ? ⊗(y,len) : mod )

"""
    Unordered_Design{Tm <: AbstractMod,Tcon <: AbstractConstrains,Tspace <: AbstractSpace{T} where T} <: Single_Design{Tspace}

Unordered space is where the position in the construct has a non-functional role, the construct [a,b,c] is seen equal to [c,b,a], so only one of both will be generated.
Al modules can only be used ones.
# Arguments
- `mod::AbstractMod`: The used moduels in this design
- `len::Int` : The allowed length of the construct
- `con::AbstractConstrains`:: The constrains of the design space
- `space::AbstractSpace`:: Design space containing all constructs, if possible, these constructs are not explicitly generated.
"""
struct Unordered_Design{Tm <: AbstractMod,Tcon <: AbstractConstrains,Tspace <: AbstractSpace{T} where T} <: Single_Design{Tspace}
    mod::Tm
    len::Int
    con::Tcon
    space:: Tspace
end

"""
    make_unordered_space(mod::Group_Mod,len::Int)

The function generates all unordered constructs with given length and modules. Every module can only be used one single time in a construct.
The constructs aren't explicitly generated.
With the use of the ``Combination`` structure, every construct can be made on the fly only on the moment it is required.
"""

make_unordered_space(mod::Group_Mod,len::Int) =  len > 1 ? Combination(mod.m,len) : mod

"""
    construct_design(mod::Group_Mod, len::Int ; con::No_Constrain = No_Constrain(nothing) , order::Bool = false)

Returns a correct filled `Unordered_Design` with given input modules and desired length.
If `order = true`, then an `Ordered_Design` is returned using the same input settings
See [`Unordered_Design`](@ref) or [`Ordered_Design`](@ref)
"""
function construct_design(mod::Group_Mod, len::Int ; order::Bool = false)
    if order == true
        make_ordered_space(mod,len) |> Full_Ordered_space |> (y -> Ordered_Design(mod,len,No_Constrain(nothing),y))
    else
        make_unordered_space(mod,len) |> Full_Unordered_space |> (y -> Unordered_Design(mod,len,No_Constrain(nothing),y))
    end
end

# precomuted full space all other field are ignored, currently not used , manby usefull in the future.
construct_design(Tspace::Computed_Space ; mod = Group_Mod(nothing), len = Len_Constrain(nothing), con = No_Constrain(nothing), order= false) =  order == true ?  Ordered_Design(mod,len,con,Tspace) : Unordered_Design(mod,len,con,Tspace)

"""
    construct_design(mod::Group_Mod, len::Int, con::Construct_Constrains ; order = false)

Returns a correct filled `Unordered_Design` with given input modules and desired length.
The constraints are attributed and evaluated if they are useful for the given constructed type.
Ordered constraints are only helpful if  `order = true`, then an `Ordered_Design` is returned
See [`Unordered_Design`](@ref) or [`Ordered_Design`](@ref)
"""
function construct_design(mod::Group_Mod, len::Int, con::Construct_Constrains{T} where T; order = false)
    if order == true
        make_ordered_space(mod,len) |> x -> Frame_Space(Full_Ordered_space(x),con) |> (y -> Ordered_Design(mod,len,con,y))
    else
        if isa(con,UnOrdered_Constrain) || promote_type(eltype(con),UnOrdered_Constrain) == UnOrdered_Constrain
            return make_unordered_space(mod,len) |> x -> Frame_Space(Full_Unordered_space(x),con) |> (y -> Unordered_Design(mod,len,con,y))
        else
            error("Ordered contrains can't be used in unorderd design set order to true or remove constrains")
        end
    end
end

"""
    getspace(Design::Single_Design; full = false)

If the input design has no constrains, then an efficient space is generated This space type allow most `Base` function that you would have on the explicit generated design space.
If constraints are added than no efficient space can be generated. For most functionalities explicitly, generation is needed. In some case one can continue using the whole design structure.
If `full = true` the whole space is generated explicitly, which isn't recommended for large designs. Still for constrained problems is sometimes the best or only option.
"""

function getspace(Design::Single_Design; full = false)
          if full== false
              if  isa(Design.space, Frame_Space)
                  @warn "given space is the closed effienct space that can be made, not all constrains are used"
              end
              return Design.space
          else
            return  created_space(Design.space,Design.con)
        end
end

### Create space can be removed in the future because Base.iterate{Frame_space} has same functionalities.
#was implented before and it works so kept for now.

"""
    created_space(space::Eff_Space)

Generated explicated all constructed for a space that can be handled efficiently
"""
function created_space(space::Eff_Space,::AbstractConstrains)
    new_space = [construct for construct in space]
    return Computed_Space(new_space)
end

"""
    created_space(space::Frame_space,con::Construct_costrains)

It generated the space that has one or multiple constraints.
The constructed that aren't allowed base on given constraints are filtered out.
"""
function created_space(space::Frame_Space,con::Construct_Constrains)
    new_space = eltype(space)[]
    for temp_construct in space.space
        #print(temp_construct)
            if filter_constrain(temp_construct,con) == false
                 push!(new_space,temp_construct)
            end
   end
    return Computed_Space(new_space)
end

"""
    created_space(space::Comuputed_space)

fall back if space was computed before.
"""
created_space(space::Computed_Space,::AbstractConstrains) = space.space

"""
    constrain_check(len_construct::Int,constrain)

quick check if the constrain make sense given the length of the constructs.
"""
constrain_check(len_construct::Int,constrain) = (len_construct >= length(constrain.combination)) && (len_construct >= maximum(constrain.pos))

""""
    filter_constrain(construct::Ordered_Construct, con::Ordered_Constrain)

Evaluated if the `construct` is allowed given  `con`.  Returns `true` when the construct isn't allowed.
Because the constrain is ordered, the given modules are only tested on the given position. Only with an exact match, it returns true,if no match is found it returns `false`.
"""
function filter_constrain(construct::Ordered_Construct, con::Ordered_Constrain)
    # check if constrains makes sence to given construct
    if constrain_check(length(construct),con) == false
        @warn "skiped constrain that didn't make sence"
        return false
    else
         for (pos,element) in  zip(con.pos,con.combination)
             if construct[pos] != element
                 return false
            end
        end
    end
    return true
end

""""
    filter_constrain(construct::AbstractConstruct, con::UnOrdered_Constrain)

Evaluated if the `construct` is allowed given  `con`.  Returns `true` when the construct isn't allowed.  If no match is found it returns `false`
Because the constrain is unordered only all modules in the constrain need to be present in construct.
The position is not evaluated.
"""

filter_constrain(construct::AbstractConstruct, con::UnOrdered_Constrain) = (Set(construct.c) ∩ Set(con.combination)) |> length |> (y -> y == length(con.combination))

""""
    filter_constrain(construct::AbstractConstruct, con::Compose_Construct_Constrains)

A wrapper to evaluated multiple constraints regarding the same constructs.
The function terminates if one of the constrain matches to the given constructs and returns `true`.
"""

function filter_constrain(construct::AbstractConstruct, con::Compose_Construct_Constrains)
    for temp_con in con.construct_con
        if filter_constrain(construct,temp_con)
            return true
        end
    end
    return false
end

"""
    Multi_Design{Tm <: Single_Design{T} where T} <: AbstractDesign{T}

Designs that are constructed out of different `Single_Design`.
The underlining Designs are stored  in a Vector.
Currently only used to allow multiple lengths of constructs for the same design rules ( modules and constraints ).
"""

struct Multi_Design{T} <: AbstractDesign{T}
    d::Vector{T}
end

"""
addition of differnt designs
"""
Base.eltype(::Multi_Design{T}) where T = T
Base.length(design::Multi_Design) = length(design.d)
Base.iterate(design::Multi_Design,state = 1) = length(design.d) >= state ? (design[state],state+1) : nothing
Base.getindex(design::Multi_Design,i) = design.d[i]
Base.size(design::Multi_Design) = (length(design), map((i -> length(i.space,nothing)),design) |> sum )

"""
    construct_design(mod::Group_Mod, min::Int, max::Int; order = false)

Returns Mulit_Design contianing the design space for ever length between `min` and `max`.
If `order = true`, then an `Ordered_Design` is constructed.
"""

construct_design(mod::Group_Mod, min::Int, max::Int; order = false) = Multi_Design([construct_design(mod,len,order = order) for len in min:max])

"""
    construct_design(mod::Group_Mod, min::Int, max::Int; order = false)

Returns Mulit_Design contianing the design space for ever length in the `lens` array.
If `order = true`, then an `Ordered_Design` is constructed.
"""

construct_design(mod::Group_Mod, lens::Array{Int}; order = false) = Multi_Design([construct_design(mod,len,order = order) for len in lens])

"""
    construct_design(mod::Group_Mod, min::Int, max::Int; order = false)

Returns Mulit_Design contianing the design space for ever length in the `lens` array.
If `order = true`, then an `Ordered_Design` is constructed.
"""

construct_design(mod::Group_Mod, lens::Array{Int}, con::Construct_Constrains{T} where T; order = false) = Multi_Design([construct_design(mod,len,con,order = order) for len in lens])

"""
    getspace(multiple_design::Multi_Design; full = false )

Simple wrapper to obtain all different space in the Multi_Design object using `getspace(Design::Single_Design; full = false)` iteratively.
type determined on the space of the first design
"""

function getspace(design::Multi_Design; full = false )

    if full
        all_spaces = collect(design[1].space)
        for i in 2:length(design)
            append!(all_spaces,collect(design[i].space))
        end

        return Computed_Space(all_spaces)

    else
        all_spaces = AbstractSpace[]
        for d in design
             push!(all_spaces,getspace(d))
        end
        return Multi_Space(all_spaces)
    end

end

# Combining of combinations if not all intermediate-length are desired.

Base.:+(des1::T,des2::T) where T = Multi_Design([des1,des2])
Base.:+(des1::Multi_Design{T} where T, des2::Single_Design) = Multi_Design([des1.design...,des2])
Base.:+(des1::Single_Design,des2::Multi_Design) = +(des2,des1)
Base.:+(des1::Multi_Design{T} where T, des2::Multi_Design{N} where N) = Multi_Design([des1.design...,des2.design...])


text:
version = v"13.0.0"

include("../common.jl")

build_tarballs(ARGS, configure_build(ARGS, version; experimental_platforms=true)...;
               preferred_gcc_version=v"7", preferred_llvm_version=v"9", julia_compat="1.8")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# Tests for SPMF-code

using NonlinearEigenproblemsTest
using NonlinearEigenproblems
using Test
using LinearAlgebra
using Random
using SparseArrays

@testset "SPMF" begin
    # Create an SPMF
    n=5;
    Random.seed!(1)
    A0=sparse(randn(n,n));
    A1=sparse(randn(n,n));
    t::Float64=3.0

    minusop = S -> -S
    oneop = S -> one(S)
    expmop = S -> exp(-t*S)
    fi=[minusop, oneop, expmop];

    J = SparseMatrixCSC(1.0I, n, n)
    nep1 = SPMF_NEP([J, A0, A1], fi)
    nep2 = SPMF_NEP([J, A0, A1], fi; Schur_fact=true)

    @bench @testset "compute_Mlincomb" begin
        # test that the function compute_Mlincomb does not overwrite the input
        λ =randn(); V=randn(n,3); W=copy(V); a=[1; 0; 2];
        z=compute_Mlincomb(nep1,λ,V,a);
        @test norm(V-W)<sqrt(eps())*100
        # test that the function compute_Mlincomb! overwrites the input
        z=compute_Mlincomb!(nep1,λ,V,a);
        @test norm(V-W)>sqrt(eps())*100
    end

    @bench @testset "compute_MM" begin

        S=randn(3,3);
        V=randn(n,3);

        # Check if prefactorize with Schur gives the same result
        @test opnorm(compute_MM(nep1,S,V)-compute_MM(nep2,S,V))<sqrt(eps())
        # Check compute_MM
        @test opnorm(compute_MM(nep1,S,V)-(-V*S+A0*V+A1*V*exp(-t*S)))<sqrt(eps())

        # Check if compute_MM is correct (by comparing against diagonalization of S).
        # This is done with the identity MM(V,S)=MM(V,W*D*inv(W))=MM(V*W,D)*inv(W)
        # and MM(V1,D) for diagonal D can be computed with compute_Mlincomb

        N1=compute_MM(nep1,S,V);
        # Diagonalize S
        d,W = eigen(S)
        D = diagm(0 => d)
        V1 = V*W
        #
        N2=hcat(compute_Mlincomb(nep1,d[1],V1[:,1]),
                compute_Mlincomb(nep1,d[2],V1[:,2]),
                compute_Mlincomb(nep1,d[3],V1[:,3]))*inv(W)
        @test opnorm(N1-N2)<sqrt(eps())

    end
#    @bench @testset "SPMF sparse" begin
#
#
#        Random.seed!(0)
#        B0=sprandn(5,5,0.2)
#        B1=sprandn(5,5,0.2)
#
#        fv=[S->one(S),S->sin(S)];
#        spmf1=SPMF_NEP([B0,B1],fv,align_sparsity_patterns=false)
#        spmf2=SPMF_NEP([copy(B0),copy(B1)],fv,align_sparsity_patterns=true)
#
#        spmf3=SPMF_NEP([Matrix(B0),Matrix(B1)],fv)
#
#        λ=1.0;
#        Ma1=compute_Mder(spmf1,λ);
#        Ma2=compute_Mder(spmf2,λ);
#        Ma3=compute_Mder(spmf3,λ);
#        @test Ma1 ≈ Ma2
#        @test Matrix(Ma1) ≈ Matrix(Ma3)
#
#
#        λ=1.0+3im;
#        Mb1=compute_Mder(spmf1,λ);
#        Mb2=compute_Mder(spmf2,λ);
#        Mb3=compute_Mder(spmf3,λ);
#        @test Matrix(Mb1) ≈ Matrix(Mb2)
#        @test Matrix(Mb1) ≈ Matrix(Mb3)
#
#
#    end
    @bench @testset "compute_Mder_from_MM" begin

        S=randn(3,3);
        V=randn(n,3);

        # Check compute_Mder_from_MM()
        λ=2
        T1=compute_Mder_from_MM(nep1,λ,1)
        T2 = -I - t*A1*exp(-t*λ)
        @test opnorm(T1-T2)<sqrt(eps())

        λ=2
        T3=compute_Mder_from_MM(nep1,λ,2)
        T4=t^2*A1*exp(-t*λ)
        @test opnorm(T3-T4)<sqrt(eps())

        # Check consistency of MM and Mder

        # Exact Mlincomb:
        Zexact = (-λ*I + A0 + A1*exp(-t*λ)) * V[:,1] +
            (-I - t*A1*exp(-t*λ)) * V[:,2] +
            (t^2*A1*exp(-t*λ)) * V[:,3]

        Z1=compute_Mlincomb_from_MM(nep1,λ,V,[1.0,1.0,1.0])
        @test norm(Z1-Zexact)<sqrt(eps())

        Z2=compute_Mlincomb_from_Mder(nep1,λ,V,[1.0,1.0,1.0])
        @test norm(Z2-Zexact)<sqrt(eps())

    end

    @bench @testset "Compute_Mlincomb" begin

        Random.seed!(10)
        S=randn(3,3)+120^2*I;
        V=randn(n,3);

        # Same nonlinearities as the GUN NLEVP problem
        minusop= S-> -S
        oneop = S -> one(S)
        sqrt1op = S -> 1im * sqrt(S)
        sqrt2op = S -> 1im * sqrt(S - 108.8774^2*one(S))

        A0=sparse(randn(n,n))+1im*sparse(randn(n,n));
        A1=sparse(randn(n,n))+1im*sparse(randn(n,n));
        A2=sparse(randn(n,n))+1im*sparse(randn(n,n));
        A3=sparse(randn(n,n))+1im*sparse(randn(n,n));
        nep2=SPMF_NEP([A0,A1,A2,A3],[oneop,minusop,sqrt1op,sqrt2op]);

        N1=compute_MM(nep2,S,V);
        # Diagonalize S
        d,W = eigen(S)
        D = diagm(0 => d)
        V1 = V*W
        #
        N2=hcat(compute_Mlincomb(nep2,d[1],V1[:,1]),
                compute_Mlincomb(nep2,d[2],V1[:,2]),
                compute_Mlincomb(nep2,d[3],V1[:,3]))*inv(W)
        @test norm(N1-N2)<sqrt(eps())*100

    end

    @bench @testset "PEP" begin
        Random.seed!(99)
        A0=randn(5,5)
        A1=randn(5,5)
        A2=randn(5,5)
        A3=randn(5,5)
        # Check that PEP generates the
        # same values as an SPMF generated from the PEP
        nep0=PEP([A0,A1,A2,A3])
        nep1=SPMF_NEP(get_Av(nep0),get_fv(nep0))
        for  λ in [3,10,-3,-100]
            M0=compute_Mder(nep0,λ)
            M1=compute_Mder(nep1,λ)
            @test M0 ≈ M1

            M0=compute_Mder(nep0,λ,2)
            M1 = 2 .* A2 + (λ*6) .* A3
            @test M0 ≈ M1
        end

    end

    @onlybench @testset "SPMF benchmark" begin
        # To check performance of SPMF-compute_MM function
        large_benchmark = false
        if large_benchmark
            n_values = (5, 10, 50)
            p_mm_values = (5, 50)
            p_mder_values = (1, 3, 10, 20)
            m_terms = 500
        else
            n_values = (5, 10)
            p_mm_values = (5, 10, 15, 100)
            p_mder_values = (1, 3, 10)
            m_terms = 20
        end

        for Schur_fact in (true, false)
            for n in n_values
                # small number of terms.
                fv=[S->S, S->cos(S)]
                A0=randn(n,n);
                A1=randn(n,n);
                spmf1=SPMF_NEP([A0,A1],fv,Schur_fact = Schur_fact,Ftype=Float64);
                @testset "two-term SPMF (n=$n,schur=$Schur_fact): MM S-matrix p x p: p=$p" for p in p_mm_values
                    V=randn(size(spmf1,1),p);
                    S=randn(p,p);
                    Z=compute_MM(spmf1,S,V);
                    @test eltype(Z) == promote_type(eltype(S),eltype(V),eltype(A0))
                end
                @testset "two-term SPMF (n=$n,schur=$Schur_fact): Mder($p)" for p in p_mder_values
                    λ=3.0+1im;
                    MM=compute_Mder(spmf1,λ,p);
                    @test eltype(MM) == promote_type(typeof(λ),eltype(A0))
                end

                # large number of terms
                m=m_terms;
                fv=Vector{Function}(undef,m);
                Av=Vector{Matrix{Float64}}(undef,m);
                for k=1:m
                    fv[k] = S-> inv(sqrt(k)*I-S); # 1/(I√k-S)
                    Av[k] = randn(n,n)
                end

                spmf2=SPMF_NEP(Av,fv,Schur_fact = Schur_fact, Ftype=Float64);
                @testset "$m-term SPMF (n=$n,schur=$Schur_fact): MM S-matrix p x p: p=$p" for p in p_mm_values
                    p=5;
                    V=randn(size(spmf2,1),p);
                    S=randn(p,p);
                    Z=compute_MM(spmf2,S,V);
                    @test eltype(Z) == promote_type(eltype(S),eltype(V),eltype(A0))
                end
                # Disabled until it is optimized.
                #@testset "$m-term SPMF (n=$n,schur=$Schur_fact): Mder($p)" for p in p_mder_values
                #    λ=3.0+1im;
                #    MM=compute_Mder(spmf2,λ,p);
                #    @test eltype(MM) == promote_type(typeof(λ),eltype(A0))
                #end
            end
        end
    end

    @bench @testset "REP" begin
        Random.seed!(10)
        A0=randn(5,5)
        A1=randn(5,5)
        A2=randn(5,5)
        # Check that REP generates the
        # same values as an SPMF generated from the REP
        d1=1;
        d2=2;
        d3=3.3;
        c1=1.5;
        c2=2.4;
        c3=9.0;
        c4=-9.9
        nep0=REP([A0,A1],[c1; c2; c3; c4], [d1; d2; d3])
        nep1=SPMF_NEP(get_Av(nep0),get_fv(nep0))
        for  λ in [3,10,-3,-100]
            M0=compute_Mder(nep0,λ)
            M1=compute_Mder(nep1,λ)
            @test M0 ≈ M1
        end

    end

    @bench @testset "DEP" begin
        Random.seed!(88)
        A1=randn(5,5)
        A2=randn(5,5)
        τ1::Float64 = 1.5
        τ2::Float64 = 3.75
        # Check that DEP generates the
        # same values as an SPMF generated from the DEP
        nep0=DEP([A1,A2],[τ1, τ2])
        nep1=SPMF_NEP(get_Av(nep0),get_fv(nep0))
        for  λ in [3,10,-3,-100]
            M0=compute_Mder(nep0,λ)
            M1=compute_Mder(nep1,λ)
            @test M0 ≈ M1
            M2 = -λ*I + A1 .* exp(-τ1*λ) + A2 .* exp(-τ2*λ)
            @test M0 ≈ M2
        end
    end
end


text:
include("../../wannierize_utils.jl")

#=
We look at the free Hamiltonian H = -Delta in the unit cell 2pi x [-.5,.5]^d, with BZ [0,1]^d
For every k, we therefore solve the Hamiltonian H_k (-i nabla + k)^2 on 2pi [-.5,.5]^d with periodic BC
This is discretized in the basis of K/G vectors e_K(r) = exp(iK.r), with K integers
In this basis, H_KK' = delta_KK' |K|^2
Therefore unk(r) = e^(i K_n r), where K_n corresponds to the n-th element of the K ordered by increasing value
if perm = sortperm(eig_k), then perm gives the mapping from n to K
Ki is limited to -Li:Li. By Fourier duality this gives a real-space grid of the unit cell

i,j,k are indices in the BZ (from 1 to N1,N2,N3), kind is the flat version of (i,j,k) (from 1 to N1*N2*N3).
The BZ loops are always on [0,1]^d, which is shifted by "corner" when computing λ, so that the actual BZ is corner + [0,1]^d
=#

using SpecialFunctions
nwannier = 2
# k-grid
N1 = 1
N2 = 1
N3 = 40

# G-grid
L1 = 0
L2 = 0
L3 = 2

dim = 3 - count(x->x==1, [N1,N2,N3])

#SCDM parameters
mu = 0.0;
sigma = 2.0;

do_shift = false
corner = [0.,0.,0.]
# corner = [0.,0.,0.] # do not use directly in the code, only the functions below

# all loops in the code refer to a k point grid like (0:N-1)/N. This transfers to the actual kpt grid
function k_red_to_real(k)::Vector{Float64} #uses non-constant globals, so use explicit type
    if do_shift
        shift = [1/N1/2, 1/N2/2, 1/N3/2]
        shift[1:3-dim] = 0 #don't shift singleton dimensions
        k+shift+corner
    else
        k+corner
    end
end
function k_real_to_red(k)::Vector{Float64}
    if do_shift
        shift = [1/N1/2, 1/N2/2, 1/N3/2]
        shift[1:3-dim] = 0 #don't shift singleton dimensions
        k-shift-corner
    else
        k-corner
    end
end

k = randn(3); @assert k ≈ k_red_to_real(k_real_to_red(k)) #sanity check

λ(k,Ks) = squeeze(sum(abs2, (Ks .+ k_red_to_real(k)), 1),1) #eigenvalues at k in [0,1]^d

nband = (2L1+1)*(2L2+1)*(2L3+1) #size of the Hamiltonian
Ks = zeros(Int,3,nband)
ind = 1
for l1=-L1:L1
    for l2=-L2:L2
        for l3=-L3:L3
            Ks[:,ind] = [l1, l2, l3]
            ind += 1
        end
    end
end

function write_files(N1, N2, N3, L1, L2, L3, nwannier)
    nband = (2L1+1)*(2L2+1)*(2L3+1) #size of the Hamiltonian
    Ks = zeros(Int,3,nband) #list of K vectors
    Ntot = N1*N2*N3
    Ns = [N1,N2,N3]
    R = zeros(Complex128,nband,3)
    Rfrac = zeros(Complex128,nband,3)
    ind = 1
    for l1=-L1:L1
        for l2=-L2:L2
            for l3=-L3:L3
                Ks[:,ind] = [l1, l2, l3]
                R[ind,:]  = 2*pi*[l1/(2*L1+1), l2/(2*L2+1), l3/(2*L3+1)]
                Rfrac[ind,:] = [(l1+L1)/(2*L1+1), (l2+L2)/(2*L2+1), (l3+L3)/(2*L3+1)]
                ind += 1
            end
        end
    end

    win = open("free.win","w")
    write(win,"""
num_bands = $nband
num_wann = $nwannier

begin unit_cell_cart
bohr
$(2pi) 0 0
0 $(2pi) 0
0 0 $(2pi)
end unit_cell_cart

mp_grid = $N1 $N2 $N3

begin kpoints
""")
    kind = 1
    kind_to_ijk = zeros(Int64,Ntot,3) #kind runs from 1:Ntot, i,j,k from 1:N1,N2,N3
    ijk_to_kind = zeros(Int64,N1,N2,N3)
    kpt_order = [1,2,3] #potentially allow different ordering
    for i=0:N1-1
        for j=0:N2-1
            for k=0:N3-1
                write(win, "$(i/N1) $(j/N2) $(k/N3)\n")

                kind += 1
                ijk = [i,j,k]
                big_offset = Ns[kpt_order[1]]*Ns[kpt_order[2]]
                med_offset = Ns[kpt_order[2]]
                kind = ijk[kpt_order[1]]*big_offset + ijk[kpt_order[2]]*med_offset + ijk[kpt_order[3]]+1
                ijk_to_kind[i+1,j+1,k+1] = kind
                kind_to_ijk[kind,:] = [i+1 j+1 k+1]

            end
        end
    end

    write(win,"end kpoints")
    close(win)

    neighbors =
        [1 0 0;
         0 1 0;
         0 0 1;
         -1 0 0;
         0 -1 0;
         0 0 -1]

    # overlap matrix Mmnkb between k and k+b, with k+b = displ_vec + kpb
    # <un(k)|um(kpb + displ_vec)> = <un(k)|e^-(displ_vec)x um(kpb)>
    # = delta_KK', but mn are sorted versions of K K'
    function get_ovl(k, kpb, displ_vec,Ks)
        eig_k = λ(k,Ks)
        eig_kpb = λ(kpb,Ks)

        sort_k = sortperm(eig_k)
        sort_kpb = sortperm(eig_kpb)
        A = zeros(Int64, size(Ks,2),size(Ks,2))
        for m=1:nband
            for n=1:nband
                Km = sort_k[m]
                @assert eig_k[Km] == sort(eig_k)[m]
                Kpbn = sort_kpb[n]
                A[m,n] = Ks[:,Km] == (Ks[:,Kpbn] .- displ_vec)
            end
        end
        return A
    end

    eig = open("free.eig","w")
    eig_allk = zeros(Complex128,nband,N1*N2*N3)
    ind = 1
    for i=0:N1-1
        for j=0:N2-1
            for k=0:N3-1
                ijk = [i,j,k]
                kind = ijk_to_kind[i+1,j+1,k+1]
                eig_k = sort(λ(ijk./Ns,Ks))
                eig_allk[:,ind] = eig_k
                ind += 1
                for iband=1:nband
                    write(eig,"$iband $kind $(eig_k[iband])\n")
                end
            end
        end
    end
    close(eig)

    mmn = open("free.mmn","w")
    write(mmn, "Created by free.jl $L1 $L2 $L3\n")
    nneighbors = count(x -> x != 1, Ns)*2
    write(mmn, "$nband $Ntot $nneighbors\n")
    for i=0:N1-1
        for j=0:N2-1
            for k=0:N3-1
                ijk = [i,j,k]
                for ib = 1:6
                    b = neighbors[ib,:]
                    ijkpb_recen = mod.(ijk+b,[N1; N2; N3]) #bring back to BZ
                    kind = ijk_to_kind[(ijk.+1)...]
                    kindb = ijk_to_kind[(ijkpb_recen.+1)...]
                    if kind == kindb
                        continue #ignore singleton dimensions
                    end
                    displ_vec = div.((ijk+b) - ijkpb_recen,Ns)
                    write(mmn,"$kind $kindb $(displ_vec[1]) $(displ_vec[2]) $(displ_vec[3])\n")
                    ovl = get_ovl(ijk./Ns,(ijkpb_recen)./Ns,displ_vec,Ks)
                    for n=1:nband
                        for m=1:nband
                            write(mmn,string(ovl[m,n]), " 0\n")
                        end
                    end
                end
            end
        end
    end
    close(mmn)

    # Unk and SCDM stuff
    # same for all kpoints
    Amn = zeros(Complex128,N1,N2,N3,nband,nwannier)

    gamma = [0.0,0.0,0.0]
    eig_k = λ(k_real_to_red(gamma),Ks) # λ is computed in reduced coordinates
    Unk = exp.(1im*R*Ks)
    occ = (1/2)*erfc.((real(eig_k)-mu)/sigma)

    Qtemp, Rtemp, piv = qr(diagm(occ)*Unk',Val{true})

    cols = piv[1:nwannier]
    min_S = Inf
    for i=0:N1-1
        for j=0:N2-1
            for k=0:N3-1
                kpt = [i,j,k]./Ns

                phase = exp.(1im*R[cols,:]*k_red_to_real(kpt))

                eig_k = λ(kpt,Ks)
                n_to_K = sortperm(eig_k)
                K_to_n = invperm(n_to_K)
                eig_k = eig_k[n_to_K]
                # Note that the order of Ks is kpt dependent
                Unk = zeros(Complex128,nband,nband)
                for m = 1 : nband
                    Km = n_to_K[m]
                    Unk[:,m] = exp.(1im*R*(Ks[:,Km]))
                end
                occ = (1/2)*erfc.((real(eig_k)-mu)/sigma)

                Y = diagm(phase)*Unk[cols,:]*diagm(occ)
                Y = Y'

                U,S,V = svd(Y)
                Amn[i+1,j+1,k+1,:,:] = U*V'
                min_S = min(min_S, minimum(S))
            end
        end
    end
    println("Min singular value: $min_S")

    for file in ("free.amn", "free.SCDM.amn")
        out = open(file,"w")
        write(out, "Created by free.jl ", string(now()), "\n")
        write(out, "$(nband) $(N1*N2*N3) $(nwannier)\n")
        for i=1:N1
            for j=1:N2
                for k=1:N3
                    kind = ijk_to_kind[i,j,k]
                    for n=1:nwannier
                        for m = 1:nband
                            coeff = Amn[i,j,k,m,n]
                            write(out, "$m $n $kind $(real(coeff)) $(imag(coeff))\n")
                        end
                    end
                end
            end
        end
        close(out)
    end
end

write_files(N1,N2,N3,L1,L2,L3,nwannier)

include("../../run_optim.jl")

loc = omega_loc(p,A)
# writedlm("data",loc)
println(maximum(loc))
println("Max loc = ", maximum(loc))
include("../../plot_free_wannier.jl")

max_stars_repo_name: VVitale/wannier

max_issues_repo_name: VVitale/wannier

max_forks_repo_name: VVitale/wannier


text:
export Reconstruction, addReco, getReco, getRecons, findReco

# TODO: make store a member and remove path
mutable struct Reconstruction
  path::String
  num::Int64
  params::Dict
end

function getReco(d::MDFDatasetStore, study::Study, exp::Experiment, recoNum::Int64)
  path = joinpath(d.path, "reconstructions", getMDFStudyFolderName(study), string(exp.num), string(recoNum))
  filename = path*".mdf"
  if !isfile(filename)
    filename = path*".hdf"
  end
  r = Reconstruction(filename, recoNum, Dict())
  loadParams(r)
  return r
end

function extendPath(d::MDFDatasetStore, p::AbstractString)
  if isabspath(p)
    return normpath(p)
  else
    A = joinpath(d.path, p)
    if isfile(A)
      return normpath(A)
    end
    B = joinpath(d.path, "..", p)
    if isfile(B)
      return normpath(B)
    end
    D = splitpath(p)
    if length(D) > 1
      C = joinpath(d.path, D[2:end]...)
      @info D
      @info d.path
      @info C
      if isfile(C)
        return normpath(C)
      end
    end
    error("extendPath could not generate valid path. Something is wrong with
           the dataset store.")
  end
end

function normalizePathsRecoDict!(d::MDFDatasetStore, recoParams::Dict)
  for key in [:measPath, :SFPath, :emptyMeasPath]
    if haskey(recoParams, key)
      if isa(recoParams[key], AbstractString)
        recoParams[key] = normpath(extendPath(d, recoParams[key]))
      elseif isa(recoParams[key], AbstractVector)
        for l=1:length(recoParams[key])
          recoParams[key][l] = normpath(extendPath(d, recoParams[key][l]))
        end
      end
    end
  end
end

function splitOffStorePath(path, foldername)
  s = splitpath(path)
  i = findfirst(s .== foldername)
  if i != nothing
    return joinpath(s[i:end]...)
  else
    return path
  end
end

function makeRelativePathsRecoDict!(recoParams::Dict)
  coll = Dict((:measPath => "measurements"), (:SFPath => "calibrations"),
              (:emptyMeasPath => "measurements"))

  for (key,foldername) in coll
    if haskey(recoParams, key)
      if isa(recoParams[key], AbstractString)
        recoParams[key] = splitOffStorePath(recoParams[key],foldername)
      elseif isa(recoParams[key], AbstractVector)
        for l=1:length(recoParams[key])
          recoParams[key][l] = splitOffStorePath(recoParams[key][l],foldername)
        end
      end
    end
  end
end

# This functions searches for recoparams and returns the corresponding recoNumber
# 0 indicates that the set of parameters was not found
function findReco(d::MDFDatasetStore, study::Study, exp::Experiment, recoParams::Dict)
  recoNum = 0

  recoParams_ = deepcopy(recoParams)
  normalizePathsRecoDict!(d, recoParams_)

  # We do not care if the reconstruction has been done by a different
  # user. Therefore, we remove the :reconstructor field
  if haskey(recoParams_, :reconstructor)
    delete!(recoParams_, :reconstructor)
  end
  recons = getRecons(d, study, exp)
  for reco in recons
    normalizePathsRecoDict!(d, reco.params)

    if haskey(reco.params, :reconstructor)
      delete!(reco.params, :reconstructor)
    end

    if recoParams_ == reco.params
      recoNum = reco.num
    end
  end

  return recoNum
end

function getReco(d::MDFDatasetStore, study::Study, exp::Experiment, recoParams::Dict)
  getReco(d, study, exp, findReco(d,study,exp,recoParams) )
end

# The following function is certainly not ideal when considering a "getReco" scenario
function addReco(d::MDFDatasetStore, study::Study, exp::Experiment, image)

  outputpath = joinpath(d.path, "reconstructions", getMDFStudyFolderName(study), string(exp.num))
  # create data directory
  mkpath(outputpath)
  try_chmod(outputpath, 0o777, recursive=true)

  recoNum = getNewNumInFolder(outputpath)

  filepath = joinpath(outputpath, string(recoNum))

  makeRelativePathsRecoDict!(image.recoParams)

  saveRecoData(filepath*".mdf", image)
  #save(filepath*".jld","recoParams",recoParams)
end

function remove(reco::Reconstruction)
  if isfile(reco.path)
    rm(reco.path)
  end
  visufile = getVisuPath(reco)
  if isfile(visufile)
    rm(visufile)
  end
end

function save(reco::Reconstruction)
  h5open(reco.path, "r+") do file
    if haskey(file, "/reconstruction/_parameters")
      delete_object(file, "/reconstruction/_parameters")
    end
    saveParams(file, "/reconstruction/_parameters", reco.params)
  end
end

function loadParams(reco::Reconstruction)

  if isfile(reco.path)
   h5open(reco.path, "r") do file
    g = file["/reconstruction"]
    if haskey(g, "_parameters") #new world order
      reco.params = loadParams(reco.path, "/reconstruction/_parameters")
    else #this needs to go
      @debug "opening legacy file"
      prefix, ext = splitext(reco.path)
        if isfile(prefix*".jld")
      	  reco.params = load(prefix*".jld","recoParams")
	else
	  error("unable to load reco params")
	end
    end
   end
  end
  nothing
end

function getVisuPath(reco::Reconstruction)
  prefix, ext = splitext(reco.path)
  return prefix*".visu"
end

function getRecons(d::MDFDatasetStore, study::Study, exp::Experiment)

  recons = Reconstruction[]

  datadir = joinpath(d.path, "reconstructions", getMDFStudyFolderName(study), string(exp.num))

  if isdir(datadir)
    files = readdir(datadir)
    for file in files
        prefix, ext = splitext(file)
        fullfile = joinpath(datadir,file)

        if ext == ".hdf"  || ext == ".mdf"

        filename = joinpath(datadir,prefix*".mdf")
        if !isfile(filename)
            filename = joinpath(datadir,prefix*".hdf")
        end

        num = parse(Int64,prefix)
        r = Reconstruction(filename, num,  Dict())
        loadParams(r)

        push!(recons, r )
        end
    end
  end

  return recons
end


text:
@testset "symbols" begin
  r = symbols(false)
  println(r)
  r.response::Vector{Any}
  r = symbol_details(false, r.response[begin])
  println(r)
  r.response::Dict{String, Any}
end

max_stars_repo_name: rory-linehan/Gemini

max_issues_repo_name: rory-linehan/Gemini

max_forks_repo_name: rory-linehan/Gemini


text:
function relax!(f, xf, s::Solver, d::Domain; verbose = false)
    ω = s.ω
    c = CartesianIndices(d)
    l = LinearIndices(d)
    Δ = 0.0
    for I = c
        δ = xf[l[I]] - f[I]
        Δ = max(abs(δ), abs(Δ))
        f[I] += ω*δ
    end
    set_Δ!(s, Δ)
    increment!(s)
    verbose && print_convergence_info(s)
    return f
end


text:
struct Spec
    id
    reward_threshold
    nondeterministic
    tags
    max_episode_steps
end


text:
mutable struct DPout
    S
    V
    π
    α
end

mutable struct TargetData

    y::Array{Int64, 1}
    N::Array{Int64, 1}
    VAF::Array{Float64, 1}
    DF::DataFrame

    TargetData(y, N, VAF) =
    new(y, N, VAF, DataFrame(y = y, N = N, VAF = VAF))
end

mutable struct DPresults
    DF::DataFrame
    weights
    nclones
    cloneweights
    clonefrequencies
    allcloneweights
    allclonefrequencies
    posterior::DPout
    data::TargetData
end


text:
include("imports.jl")

@testset ExtendedTestSet "`AdjustContrast`" begin
    item = Image(rand(RGB{N0f8}, 64, 64))
    tfm = AdjustContrast(0.2)
    @test_nowarn apply(tfm, item)
    titem = apply(tfm, item)
    @test eltype(itemdata(titem)) == RGB{N0f8}
    testapply(tfm, Image)
    testapply!(tfm, Image)
end

@testset ExtendedTestSet "`AdjustBrightness`" begin
    item = Image(rand(RGB{N0f8}, 64, 64))
    tfm = AdjustBrightness(0.2)
    @test_nowarn apply(tfm, item)
    titem = apply(tfm, item)
    @test eltype(itemdata(titem)) == RGB{N0f8}

    testapply(tfm, Image)
    testapply!(tfm, Image)
end


text:
# Test shapes

using NumericExtensions
using Base.Test

@test mapshape(()) == ()
@test mapshape((3,)) == (3,)
@test mapshape((3,4)) == (3,4)

@test mapshape((3,), ()) == (3,)
@test mapshape((), (3,)) == (3,)
@test mapshape((3,), (3,)) == (3,)
@test_throws mapshape((3,), (4,))

@test mapshape((3,1), (3,)) == (3,1)
@test mapshape((3,), (3,1)) == (3,1)
@test_throws mapshape((3,1), (2,))
@test_throws mapshape((3,1), (3,2))

@test mapshape((3,4,5), (3,4,5)) == (3,4,5)
@test mapshape((3,4,1), (3,4)) == (3,4,1)
@test mapshape((3,4), (3,4,1,1)) == (3,4,1,1)
@test_throws mapshape((3,4), (3,4,2))
@test_throws mapshape((3,5), (3,4,1))

@test mapshape((3,4), (3,4), (3,4)) == (3,4)
@test mapshape((3,4,1), (3,4), (3,4)) == (3,4,1)
@test mapshape((3,4), (3,4,1), (3,4)) == (3,4,1)
@test mapshape((3,4), (3,4), (3,4,1)) == (3,4,1)
@test mapshape((3,4,1), (3,4), (3,4,1)) == (3,4,1)
@test_throws mapshape((3,4), (3,4), (3,5))

@test mapshape((3,4), (3,4), (3,4), (3,4)) == (3,4)
@test mapshape((3,), (3,1), (3,1,1), (3,1,1,1)) == (3,1,1,1)


text:
@testset "Univariate functions" begin
    # bspline0 from Daisy benchmarks
    f = x -> (1 - x)^3 / 6.0
    dom = Interval(-4.5, -0.3)
    xref = Interval(0.36616666666666675, 27.729166666666668)

    for solver in available_solvers
        x = enclose(f, dom, solver)
        r = relative_precision(x, xref)
        if solver == :TaylorModels
            # for this example, TaylorModels cannot tightly approximate the lhs
            @test inf(r) ≤ 30 && sup(r) ≤ 1e-10
        elseif solver == :SumOfSquares
            # solver returns negative -1.06259e-06
            @test abs(inf(r)) ≤ 1e-5 && sup(r) ≤ 1e-5
        else
            @test inf(r) ≤ 1e-5 && sup(r) ≤ 1e-5
        end
    end
end

@testset "Univariate example the Quickstart Guide" begin
    f(x) = -x^3/6 + 5x
    dom = Interval(1, 4)

    x = enclose(f, dom, :IntervalArithmetic)
    xref = Interval(-5.66667, 19.8334)
    r = relative_precision(x, xref)
    @test inf(r) ≤ 1e-5 && sup(r) ≤ 1e-5

    x = enclose(f, dom, :TaylorModels, order=4)
    xref = Interval(-4.2783, 12.7084)
    r = relative_precision(x, xref)
    @test inf(r) ≤ 1e-5 && sup(r) ≤ 1e-5

    x = enclose(f, dom, :IntervalOptimisation)
    xref = Interval(4.83299, 10.5448)
    r = relative_precision(x, xref)
    @test inf(r) ≤ 1e-5 && sup(r) ≤ 1e-5

    x = enclose(f, dom, :SumOfSquares)
    xref = Interval(4.8333, 10.541)
    r = relative_precision(x, xref)
    @test inf(r) ≤ 1e-5 && sup(r) ≤ 1e-5
end

@testset "Test univariate polynomial input" begin
    @polyvar x
    p = -x^3/6 + 5x
    dom = Interval(1, 4)

    x = enclose(p, dom)
    xref = Interval(-5.66667, 19.8334)
    r = relative_precision(x, xref)
    @test inf(r) ≤ 1e-5 && sup(r) ≤ 1e-5
end


text:
# code from https://github.com/MikeInnes/BrainForth/blob/master/src/brainfuck.jl
using Sockets

mutable struct Tape
    count::Int
    pos::Int
    tape::Vector{Int}
end

Tape() = Tape(0, 1, [0])

mutable struct Printer
    sum1::Int
    sum2::Int
    quiet::Bool
end

Printer(quiet) = Printer(0, 0, quiet)

getChecksum(p::Printer) = (p.sum2 << 8) | p.sum1

Base.getindex(t::Tape) = t.tape[t.pos]
Base.setindex!(t::Tape, v) = t.tape[t.pos] = v

Base.:(==)(a::Tape, b::Tape) = a.tape == b.tape

function left!(t::Tape)
    t.pos == length(t.tape) && t.tape[end] == 0 && pop!(t.tape)
    t.pos == 1 ? pushfirst!(t.tape, 0) : (t.pos -= 1)
    return
end

function right!(t::Tape)
    t.pos == 1 && t.tape[1] == 0 && (popfirst!(t.tape); t.pos -= 1)
    t.pos == length(t.tape) && push!(t.tape, 0)
    t.pos += 1
    return
end

inc!(t::Tape) = t[] += 0x01
dec!(t::Tape) = t[] -= 0x01

function write!(n::Int, p::Printer)
    if p.quiet
        p.sum1 = (p.sum1 + n) % 255
        p.sum2 = (p.sum2 + p.sum1) % 255
    else
        write(stdout, n)
    end
end

function interpret(t::Tape, bf, p::Printer)
    loops = Int[]
    scan = 0
    ip = 1
    @inbounds while ip <= length(bf)
        t.count += 1
        op = bf[ip]
        if op == '['
            scan > 0 || t[] == 0 ? (scan += 1) : push!(loops, ip)
        elseif op == ']'
            scan > 0 ? (scan -= 1) : t[] == 0 ? pop!(loops) : (ip = loops[end])
        elseif scan == 0
            if op == '+'
                inc!(t)
            elseif op == '-'
                dec!(t)
            elseif op == '<'
                left!(t)
            elseif op == '>'
                right!(t)
            elseif op == '.'
                write!(t[], p)
            end
        end
        ip += 1
    end
    return t
end

function notify(msg)
    try
        socket = connect("localhost", 9001)
        write(socket, msg)
        close(socket)
    catch
        # standalone usage
    end
end

function verify()
    text = """++++++++[>++++[>++>+++>+++>+<<<<-]>+>+>->>+[<]<-]>>.>
        ---.+++++++..+++.>>.<-.<.+++.------.--------.>>+.>++."""
    pLeft = Printer(true)
    interpret(Tape(), collect(text), pLeft)
    left = getChecksum(pLeft)

    pRight = Printer(true)
    for c in "Hello World!\n"
        write!(convert(Int, c), pRight)
    end
    right = getChecksum(pRight)

    if left != right
        println(stderr, "$(left) != $(right)")
        exit(1)
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    verify()
    text = open(ARGS[1]) do file
        read(file, String)
    end
    p = Printer(haskey(ENV, "QUIET"))

    notify("Julia\t$(getpid())")
    interpret(Tape(), collect(text), p)
    notify("stop")

    if p.quiet
        println("Output checksum: $(getChecksum(p))")
    end
end

max_stars_repo_name: clemenswasser/benchmarks

max_issues_repo_name: clemenswasser/benchmarks

max_forks_repo_name: clemenswasser/benchmarks


text:
module TestUtilFuncs

using Base.Test
using TimeData
using Dates

println("Running utility function tests")

#############
## narowrm ##
#############

## get input
tn = TimeData.testcase(TimeData.Timenum, 5)

## get expected outcome
dats = [Date(2010, 1, 1),
        Date(2010, 1, 2),
        Date(2010, 1, 4),
        Date(2010, 1, 5)]
df = DataFrame()
df[:prices1] = @data([100, 120, 170, 200])
df[:prices2] = @data([110, 120, 130, 150])
expOut = TimeData.Timenum(df, dats)

## test actual outcome
actOut = TimeData.narowrm(tn)
@test isequal(actOut, expOut)

## get input
dats = [Date(2010, 1, 1),
        Date(2010, 1, 2)]
df = DataFrame()
df[:prices1] = @data([100, NA])
df[:prices2] = @data([110, NA])
inputTn = TimeData.Timenum(df, dats)

## get expected outcome
dats = [Date(2010, 1, 1)]
df = DataFrame()
df[:prices1] = @data([100])
df[:prices2] = @data([110])
expOut = TimeData.Timenum(df, dats)

## test actual outcome
actOut = TimeData.narowrm(inputTn)
@test isequal(actOut, expOut)

## get input
tn = TimeData.testcase(TimeData.Timenum, 1)

## get expected outcome
expOut = TimeData.testcase(TimeData.Timenum, 1)

## test actual outcome
actOut = TimeData.narowrm(tn)
@test isequal(actOut, expOut)

#############
## anyToDa ##
#############

## findstub_vector
##----------------

@test TimeData.findstub_vector(Any[NA, NA, 3.4, 2]) == 3.4
@test TimeData.findstub_vector(Any[3.4, NA, 2]) == 3.4

## anyToDa
##--------

xx = Any[3, 4.2, NA]
@test isequal(@data([3, 4.2, NA]), TimeData.anyToDa(xx))

xx = Any[3, "hello", NA]
@test isequal(@data([3, "hello", NA]), TimeData.anyToDa(xx))

## anyToDf
##--------

xx = Any[3 4; NA 2.4]
expOut = DataFrame(a = @data([3, NA]), b = [4, 2.4])
@test isequal(expOut, TimeData.anyToDf(xx, [:a, :b]))

## wrong number of column names
@test_throws Exception TimeData.anyToDf(xx, [:a, :b, :c])

end


text:
function read_atoms(fn::IO)
	sane_readline() = strip(strip_split(readline(fn), ":")[1], '\'')
	nspecies = parse(Int, sane_readline())
	species_counter = 0

	reading_atoms = false
	atname = :nothing
	natoms = 0
	atcounter = 0
	positions = Point3{Float64}[]
	bfcmts    = Vec3{Float64}[]
	atoms     = Tuple{Symbol, Point3{Float64}, Vec3{Float64}}[]
	while species_counter < nspecies
		l = sane_readline()
		if !reading_atoms
			atname = Symbol(split(l, ".")[1])
			natoms = parse(Int, sane_readline())
			reading_atoms = true
		elseif atcounter < natoms
			sline = split(l)
			push!(positions, Point3{Float64}(parse.(Float64, sline[1:3])))
			push!(bfcmts, Vec3{Float64}(parse.(Float64, sline[4:end])))
			atcounter += 1
		elseif isempty(l)
			for (p, mag) in zip(positions, bfcmts)
				push!(atoms, (atname, p, mag))
			end
			atname = :nothing
			reading_atoms   = false
			natoms          = 0
			atcounter       = 0
			positions       = Point3{Float64}[]
			bfcmts          = Vec3{Float64}[]
			species_counter += 1
		end
	end
	return atoms
end

function elk_parse_DFTU(atoms::Vector{<:Atom}, blocknames_flaglines::Dict{Symbol, Any})
	lines = pop!(blocknames_flaglines, Symbol("dft+u"))
	t_l       = split(lines[1])
	species_start = 3
	if length(t_l) == 1
		dftu_type = parse(Int, t_l[1])
		inpdftu   = parse(Int, lines[2])
	else
		dftu_type, inpdftu = parse.(Int,  t_l)
		species_start = 2
	end
	species_names = name.(unique(atoms))
	for line in lines[species_start:end]
		sline = split(line)
		species_id, l = parse.(Int, sline[1:2])
		U, J0         = parse.(Float64, sline[3:4])
		for at in filter(x -> name(x) == species_names[species_id], atoms)
			at.dftu = DFTU{Float64}(U=U, J0=J0, l=l)
		end
	end
	return dftu_type, inpdftu
end

"""
    elk_read_input(filename; execs=[Exec("elk")], run=true, structure_name="noname")

Reads an Elk input file. The `ELK_EXEC` inside execs gets used to find which flags are allowed in this input file, and convert the read values to the correct Types.
Returns a `DFInput{Elk}` and the `Structure` that is found in the input.
"""
function elk_read_input(fn::String; execs=[Exec("elk")], run=true, structure_name="noname")
	blocknames_flaglines = Dict{Symbol, Any}()
	atoms = Atom{Float64}[]
	inputs = DFInput[]
	dir, file = splitdir(fn)

	#TODO it's a bit messy that the atoms block is not handled like this, not sure if we want to clean this up
	#blocks that require specific parsing
	special_blocks = (:avec, :tasks, Symbol("dft+u"))

	open(fn, "r") do f
		sane_readline() = strip_split(readline(f), ':')[1]
		while !eof(f)
			line = sane_readline()
			isempty(line) || line[1] == '!' && continue
			if !isempty(line)
				blockname = Symbol(line)
				info = elk_block_info(blockname)
				if blockname == :atoms
					atoms = read_atoms(f)
				elseif blockname == :wannierExtra
					@warn "Please supply the .win file instead of wannierExtra block for extracting the Wannier DFInput."

					wflags, wdata, ab, cb, proj_block = wan_read_input(f)
					continue
					# push!(inputs, DFInput{Wannier90}("wannier", dir, wflags, wdata, execs, run))
				elseif info != nothing
					blocklines = String[]
					line = sane_readline()
					while !isempty(line)
						push!(blocklines, strip(replace(strip(split(line, ':')[1]), "'" => ""), '.'))
						line = sane_readline()
					end
					if !in(blockname, special_blocks)
						blocknames_flaglines[blockname] =
						    parse_block_flags(blockname, collect(Iterators.flatten(split.(filter(x->!isempty(x), blocklines)))))
					else
						blocknames_flaglines[blockname] = filter(x->!isempty(x), blocklines)
					end
				end
			end
		end
	end
	# for i in
	#cell
	scale = haskey(blocknames_flaglines, :scale) ? last(pop!(blocknames_flaglines[:scale])) : 1.0
	cell  = uconvert.(Ang, Mat3(reshape(scale .* parse.(Float64, collect(Iterators.flatten(split.(pop!(blocknames_flaglines, :avec))))), 3, 3).*1a₀))'

    newats = [Atom{Float64, eltype(cell)}(name=x[1], element=element(x[1]), position_cryst=x[2], position_cart=cell*x[2], magnetization=x[3]) for x in atoms]
	#structure #TODO make positions be in lattices coordinates
	structure = Structure(structure_name, cell, newats)

	#different tasks
	tasks_     = pop!(blocknames_flaglines, :tasks)
	should_run = map(x -> !occursin("!", x), tasks_)
	tasks      = map(x -> strip(strip(x, '!')), tasks_)

	#The wan tasks can prbably be skipped completely and added depending on the wannier input
	wan_tasks = String[]
	for (x, run) in zip(tasks, should_run)
		if x ∈ ["601", "602", "603", "604", "605"]
			push!(wan_tasks, x)
		else
			push!(inputs, input_from_task(x, blocknames_flaglines, dir, execs, run))
		end
	end
	if haskey(blocknames_flaglines, :wannier) && !isempty(wan_tasks)
		flags = pop!(blocknames_flaglines, :wannier)
		# flags[:elk2wan_tasks] = wan_tasks
		push!(inputs, DFInput{Elk}(name="elk2wannier", dir=dir, flags=flags, execs=execs, run=true))
	end
	for f in (:ngrid, :vkloff, :plot1d, :plot2d, :plot3d)
		haskey(blocknames_flaglines, f) && pop!(blocknames_flaglines, f)
	end

	#elk2wannier

	flags = SymAnyDict()

	if haskey(blocknames_flaglines, Symbol("dft+u"))
		flags[:dftu], flags[:inpdftu] = elk_parse_DFTU(structure.atoms, blocknames_flaglines)
	end

	for (k, v) in blocknames_flaglines
		merge!(flags, v)
	end

	for i in inputs
		if isempty(i.flags)
			i.flags = deepcopy(flags)
		end
	end
	return inputs, structure
end

function find_data(flags, blocknames_flaglines)
	data    = InputData[]
	for f in flags
		if haskey(blocknames_flaglines, f)
			push!(data, InputData(f, :nothing, blocknames_flaglines[f]))
		end
	end
	return data
end

function input_from_task(task, blocknames_flaglines, dir, execs, run)
	if task ∈ ["0", "1", "2", "3", "5"]
		data = find_data((:ngridk, :vkloff), blocknames_flaglines)
	elseif task ∈ ["20", "21"]
		data = find_data((:plot1d, :plot2d, :plot3d), blocknames_flaglines)
	end
	return DFInput{Elk}(task, dir, SymAnyDict(), data, execs, run)
end

function parse(::Type{UnitRange{Int}}, l::AbstractString)
	for delim in [":", ",", "-"]
		if occursin(delim, l)
			return UnitRange(parse.(Int, split(l, delim))...)
		end
	end
end

function parse_block_flags(blockname::Symbol, lines::Vector{<:AbstractString})
	flags = SymAnyDict()
	i = 1
	totlen = length(lines)
	for flaginfo in elk_block_info(blockname).flags
		flagname = flaginfo.name
		i > length(lines) && break
		typ = eltype(flaginfo)
		if typ == String || typ == Vector{String}
			flags[flagname] = lines[i]
			i += 1
		elseif typ == Symbol
			flags[flagname] = Symbol(lines[i])
			i += 1
		else
			eltyp = eltype(typ)
			# The assumption is that the infinite repeating values are at the end of a block
			if flaginfo == last(elk_block_info(blockname).flags)
				typelen = length(typ)
				parsed = parse.(eltyp, lines[i:end])
				flags[flagname] = typ[]
				@assert length(parsed) % typelen == 0
				for j = 1:typelen:length(parsed)
					push!(flags[flagname], typ(parsed[j:j+typelen-1]...))
				end
			elseif typ <: UnitRange
				flags[flagname] = parse(typ, lines[i])
				i += 1
			else
				typelen = length(typ)
				parsed = parse.(eltyp, lines[i:i + typelen - 1])
				flags[flagname] = length(parsed) == 1 ? parsed[1] : typ(parsed)
				i += typelen
			end
		end
	end
	return flags
end

function elk_write_structure(f, structure)
	write(f, "atoms\n\t")
	species = unique(map(x->x.name, atoms(structure)))
	nspecies = length(species)
	write(f, "$nspecies\n")
	for s in species
		write(f, "\t'$s.in'\n")
		ats = filter(x->x.name == s, atoms(structure))
		write(f, "\t$(length(ats))\n")
		for a in ats
			write(f, "\t")
			for i in position_cryst(a)
				write(f, "$(round(i, digits=8)) ")
			end
			for i in magnetization(a)
				write(f, "$i ")
			end
			write(f, "\n")
		end
		write(f, "\n")
	end
	write(f, "\n")
	write(f, "avec\n\t")

	for i in eachindex(structure.cell')
		c = ustrip(uconvert(a₀, structure.cell[i]))
		if i%3 == 0
			write(f, "$c\n\t")
		else
			write(f, "$c ")
		end
	end
	write(f, "\n")
end

function construct_flag_blocks(inputs::Vector{DFInput{Elk}})
	all_flags = SymAnyDict()
	for i in inputs
		merge!(all_flags, i.flags)
		for d in data(i)
			merge!(all_flags, d.data)
		end
	end
	block_flags = SymAnyDict()

	for b in ELK_CONTROLBLOCKS
		bvals = Any[]
		for f in b.flags
			if haskey(all_flags, f.name)
				push!(bvals, pop!(all_flags, f.name))
			end
		end
		if !isempty(bvals)
			block_flags[b.name] = bvals
		end
	end
	return block_flags
end

function elk_write_flag_val(f, val)
	write(f, "\t")
	if isa(val, UnitRange)
		write(f, "$(first(val))-$(last(val))")
	elseif isa(val, AbstractVector)
		for i in val
			write(f, "$i ")
		end
	else
		if isa(val, AbstractString)
			write(f, "'$val'")
		else
			write(f, "$val")
		end
	end
	write(f, "\n")
end

function elk_write_DFTU(f::IO, structure::AbstractStructure, dftu_vals::Vector)
	write(f, "dft+u\n")
	for v in dftu_vals
		elk_write_flag_val(f, v)
	end
	species = unique(atoms(structure))
	for (i, s) in enumerate(species)
		write(f, "\t$i $(string(Elk, dftu(s)))\n")
	end
	write(f, "\n")
end

function save(inputs::Vector{DFInput{Elk}}, structure::Structure)
	tasks = map(x -> x.run ? "$(x.name)" : "!$(x.name)", filter(x->x.name != "elk2wannier", inputs))
	elk2wan_input = getfirst(x->x.name == "elk2wannier", inputs)
	if elk2wan_input != nothing && haskey(elk2wan_input.flags, :elk2wan_tasks)
		append!(tasks, pop!(elk2wan_input.flags, :elk2wan_tasks))
	end
	flags = construct_flag_blocks(inputs)
	open(joinpath(inputs[1].dir, "elk.in"), "w") do f

		write(f, "tasks\n")
		for t in tasks
			write(f, "\t$t\n")
		end
		write(f, "\n")

		elk_write_structure(f, structure)

		haskey(flags, Symbol("dft+u")) && elk_write_DFTU(f, structure, pop!(flags, Symbol("dft+u")))

		#remaining generic blocks
		for (k, v) in flags
			write(f, "$k\n")
			for v_ in values(v)
				if isa(v_, Vector{<:AbstractVector})
					for v__ in v_
						elk_write_flag_val(f, v__)
					end
				else
					elk_write_flag_val(f, v_)
				end
			end
			write(f, "\n")
		end
	end
end