text:
export Particle, Disc, collision_time,
        Plane, FinitePlane, BilliardTable, AbstractPlane, AbstractParticle,
        Sinai_billiard,
        calculate_next_collision, billiard_dynamics, initial_condition

# exports of lattice functionality:
export ParticleOnLattice, CellBoundary, billiard_dynamics_on_lattice,
        calculate_next_collision_on_lattice

export continuous_time

# particle types:

abstract AbstractParticle{T}

type Particle{T} <: AbstractParticle{T}
    x::Vector2D{T}
    v::Vector2D{T}
end

Particle(x::Vector, v::Vector) = Particle(Vector2D(x), Vector2D(v))

type ParticleOnLattice{T} <: AbstractParticle{T}
    x::Vector2D{T}
    v::Vector2D{T}
    lattice_vector::Vector2D{Int}
end

# obstacle types:

abstract Obstacle{T}

immutable Disc{T} <: Obstacle{T}
    centre::Vector2D{T}
    radius::T
end

Disc(v::Vector, r) = Disc(Vector2D(v), r)

abstract AbstractPlane{T} <: Obstacle{T}

immutable Plane{T} <: AbstractPlane{T}
    point_on_plane::Vector2D{T}  # an arbitrary point on the plane
    normal::Vector2D{T}
end

Plane(v1::Vector, v2::Vector) = Plane(Vector2D(v1), Vector2D(v2))

make_perpendicular_vector{T}(v::Vector2D{T}) = Vector2D(v[2], -v[1])

immutable FinitePlane{T} <: AbstractPlane{T}
    point_on_plane::Vector2D{T}  # an arbitrary point on the plane
    normal::Vector2D{T}
    start_pos::Vector2D{T}  # where the finite piece starts and ends
    end_pos::Vector2D{T}
    director::Vector2D{T}  # vector pointing along
end

function FinitePlane(start_pos::Vector, end_pos::Vector)

    end_pos = Vector2D(end_pos)
    start_pos = Vector2D(start_pos)

    disp = end_pos - start_pos

    normal_vec = normalize(make_perpendicular_vector(disp))

    director = disp / normsq(disp)

    FinitePlane(start_pos, normal_vec, start_pos, end_pos, director)
end

@doc """A `CellBoundary` is a boundary between neighbouring cells on a lattice.
        It is assumed to be planar.  The constructor creates an incomplete
        object, so the `other_side` field must be added later
        (once the corresponding `CellBoundary` on the other side of the next cell has been constructed).

        `lattice_increment` is a vector which says which new lattice cell a particle enters when it crosses this
        `CellBoundary`.
        """ ->
type CellBoundary{T} <: AbstractPlane{T}  # immutable not (currently?) allowed due to the partial construction
    point_on_plane::Vector2D{T}
    normal::Vector2D{T}
    lattice_increment::Vector2D{Int}
    other_side::CellBoundary{T}

    CellBoundary(point, normal, lattice_increment) = new(point, normal, lattice_increment)  # partial constructor leaving other_side undefined
end

CellBoundary{T}(point::Vector{T}, normal::Vector{T}, lattice_increment::Vector) =
    CellBoundary{T}(Vector2D(point), Vector2D(normal), Vector2D(lattice_increment))

# BilliardTable type.

@doc """A `BilliardTable` is just a list of obstacles.""" ->
immutable BilliardTable{T}
  obstacles::Vector{Obstacle{T}}
end

@doc """`get_lattice_increment` returns the lattice_increment of an obstacle. This is a zero vector
except for `CellBoundary`s.""" ->

get_lattice_increment(o::Obstacle) = Vector2D(0, 0)  # non-boundary obstacles do not change the lattice cell
get_lattice_increment(o::CellBoundary) = o.lattice_increment

@doc """Compute *time of collision* of a particle and a disc,
        assuming the particle starts outside the disc and the
        particle speed is one (using the quadratic formula).

        Returns -Inf if no collision.
        """ ->
function collision_time{T}(p::AbstractParticle{T}, disc::Disc{T})

    disp = p.x - disc.centre

    B = dot(p.v, disp)
    C = dot(disp, disp) - disc.radius^2

    discriminant = B^2 - C

    if discriminant < zero(T)
        return -convert(T, Inf)
    end

    return convert(T, -B - √discriminant)  # NB: this supposes that the particle is *outside* the disc

end

@doc """Compute normal vector to disc boundary, at a point x that must lie *on the boundary*""" ->
function normal{T}(disc::Disc{T}, x::Vector2D{T})
    return (x - disc.centre) / disc.radius
end

@doc """Check whether a given position is inside a disc.""" ->
isoutside{T}(x::Vector2D{T}, disc::Disc{T}) = norm(x - disc.centre) > disc.radius

@doc """Calculate collison time of particle with plane.""" ->
collision_time{T}(p::AbstractParticle{T}, plane::AbstractPlane{T}) =
    dot(plane.point_on_plane - p.x, plane.normal) / dot(p.v, plane.normal)

normal(plane::AbstractPlane, x) = plane.normal  # same normal vector for all positions x on plane

function collision_time{T}(p::AbstractParticle{T}, plane::FinitePlane{T})
    t_collision = dot(plane.point_on_plane - p.x, plane.normal) / dot(p.v, plane.normal)
    new_pos = p.x + t_collision*p.v

    distance = (new_pos - plane.start_pos) ⋅ plane.director

    if zero(T) ≤ distance ≤ one(T)   # if distance between 0 and 1 then current point lies on plane
        t_collision
    else
        -convert(T, Inf)
    end

end

@doc """Convention: the normal points towards the allowed part of the billiard table""" ->
isoutside(x, plane::AbstractPlane) = dot(x - plane.point_on_plane, plane.normal) > 0.0

@doc """`isoutside` finds if the particle is in the allowed part of the billiard table,
i.e. is "outside" the billiard obstacles.""" ->
function isoutside(x, table::BilliardTable)
    for obstacle in table.obstacles
        if !isoutside(x, obstacle)
            return false
        end
    end

    true
end

@doc """The two versions of `calculate_next_collision` return the information of what
the result is of doing the next collision, *without* updating the particle.""" ->
function calculate_next_collision_on_lattice{T}(p::ParticleOnLattice{T},
                                             billiard_table::BilliardTable{T},
                                             previous_obstacle_hit::Obstacle{T})

    obstacles = billiard_table.obstacles

    first_collision_time = convert(T, Inf)
    artificial_obstacle = Plane{T}([zero(T), zero(T)], [zero(T), zero(T)])
    which_obstacle_hit = artificial_obstacle
    # artificial non-existent obstacle to avoid type instability

    for obstacle in obstacles
        obstacle === previous_obstacle_hit && continue

        t = collision_time(p, obstacle)

        if 0.0 < t < first_collision_time
            which_obstacle_hit, first_collision_time = obstacle, t
        end
    end

    #@assert which_obstacle_hit != artificial_obstacle  # hit a real obstacle

    x_collision = p.x + p.v*first_collision_time
    lattice_increment = get_lattice_increment(which_obstacle_hit)

    x_new, v_new, which_obstacle_hit = collide(x_collision, p.v, which_obstacle_hit)

    return x_new, v_new, first_collision_time, which_obstacle_hit, lattice_increment
end

function calculate_next_collision(p::Particle, billiard_table, previous_obstacle_hit)

    obstacles = billiard_table.obstacles

    first_collision_time = Inf
    which_obstacle_hit = nothing

    for obstacle in obstacles
        obstacle === previous_obstacle_hit && continue

        t = collision_time(p, obstacle)

        if t > 0.0 && t < first_collision_time
            which_obstacle_hit, first_collision_time = obstacle, t
        end
    end

    x_collision = p.x + p.v*first_collision_time

    x_collision, v_new, which_obstacle_hit = collide(x_collision, p.v, which_obstacle_hit)
    # Actually only need to return new velocity

    return x_collision, v_new, first_collision_time, which_obstacle_hit
end

@doc """`collide` *implements* an elastic collision""" ->
function collide{T}(x_collision::Vector2D{T}, v::Vector2D{T}, obstacle::Obstacle{T})
    n = normal(obstacle, x_collision)

    v_new = v - 2.0*dot(n,v)*n  # reflejar solo si es disco; si es plano, pasar a traves
    v_new /= norm(v_new)
    return x_collision, v_new, obstacle
end

@doc """This version of `collide` implements a "collision" with a `CellBoundary`.
        This just moves the particle to the opposite boundary (updating its position).
        This assumes that the distance between opposite
        CellBoundary objects is 1 and the normal is a unit normal
        pointing inwards (i.e. towards the available space on the billiard table).""" ->
function collide{T}(x_collision::Vector2D{T}, v::Vector2D{T}, boundary::CellBoundary{T})

    x_new = x_collision + boundary.normal

    return x_new, v, boundary.other_side  # velocity stays the same
end

@doc """Generate a random initial condition, uniformly in the allowed region of the billiard table
("outside") the billiard obstacles. Here, "outside" for a plane is taken to mean that
the allowed region lies *in the direction of the plane's normal vector*.

Velocities are generated uniformly on the unit circle S^1, so that in total the initial condition
is generated uniformly with respect to Liouville measure.
""" ->
function initial_condition(table, xmin, xmax, ymin, ymax)
    x, y = xmin, ymin

    while true
        x = xmin + rand()*(xmax-xmin)
        y = ymin + rand()*(ymax-ymin)

        if isoutside(Vector2D(x,y), table)  # valid if outside all discs
            break
        end
    end

    # velocity in 2D:
    theta = rand() * 2*pi

    xx = Vector2D(x, y)
    vv = Vector2D(cos(theta), sin(theta))

    return xx, vv
end

@doc """Create a Sinai billiard table (a square with a disc in the centre).
It may be periodic in the x and/or y directions.""" ->
function Sinai_billiard{T}(radius::T, periodic_x=false, periodic_y=false)

    obstacles = Obstacle{T}[]

    push!(obstacles, Disc{T}([0., 0.], radius) )

    if periodic_x
        # create a pair of opposite `CellBoundary`s:
        right = CellBoundary([0.5, -0.5], [-1., 0.], [1, 0])
        left  = CellBoundary([-0.5, 0.5], [1., 0.], [-1, 0])
        right.other_side = left
        left.other_side = right
        push!(obstacles, right, left)
    else
        push!(obstacles, Plane([0.5, -0.5], [-1., 0.]) )
        push!(obstacles, Plane([-0.5, 0.5], [1., 0.]) )
    end

    if periodic_y
        top = CellBoundary([0.5, 0.5], [0., -1.], [0, 1])
        bottom  = CellBoundary([-0.5, -0.5], [0., 1.], [0, -1])
        top.other_side = bottom
        bottom.other_side = top
        push!(obstacles, top, bottom)
    else
        push!(obstacles, Plane([-0.5, -0.5], [0., 1.]) )
        push!(obstacles, Plane([0.5, 0.5], [0., -1.]) )
    end

    return BilliardTable(obstacles)
end

@doc """Simulate a single particle p on a billiard table for a given number of collisions."""->
function billiard_dynamics{T}(p::Particle{T}, table::BilliardTable{T}, num_collisions::Integer)

    xs = [p.x]

    which_obstacle_hit = nothing

    for t in 1:num_collisions
        x_collision, v_new, first_collision_time, which_obstacle_hit =
            calculate_next_collision(p, table, which_obstacle_hit)

        p.x, p.v = x_collision, v_new

        push!(xs, p.x)

    end

    xs

end

@doc """Simulate discrete-time (collision) dynamics for a single particle
        p on a billiard table for a given number of obstacle collisions
        (excluding "collisions" on boundaries) on a *periodic* billiard table.\

        p contains the initial conditions.

        Returns the positions in the unit cells (*xs*) and the corresponding lattice cells (*ls*),
        as well as the free paths."""->
function billiard_dynamics_on_lattice{T}(p::ParticleOnLattice{T}, table::BilliardTable{T}, num_collisions::Integer)

    xs = [p.x]
    lattice_vectors = [p.lattice_vector]

    which_obstacle_hit = Plane([-Inf, -Inf], [0., 0.])
    # which_obstacle_hit = nothing

    free_path_length = 0.0
    free_paths = Float64[]
    num_obstacle_collisions = 0

    first_collision_time = Inf

    #for t in 1:num_collisions
    while num_obstacle_collisions < num_collisions
        x_new, v_new, first_collision_time, which_obstacle_hit, lattice_increment =
            calculate_next_collision_on_lattice(p, table, which_obstacle_hit)

        # version to draw periodised orbits:
        # outputs the position twice, once before and once after the collision

        ## Position where it hit the obstacle at the collision: (This is a bit of a hack)
        # x_collision = p.x + p.v*first_collision_time
        # push!(xs,x_collision)
        # push!(lattice_vectors, p.lattice_vector)

        # Periodise:
        p.x, p.v = x_new, v_new
        p.lattice_vector += lattice_increment

        # push!(xs, p.x)
        # push!(lattice_vectors, p.lattice_vector)

        ## version which takes account only of collisions with objects that are not CellBoundary

        free_path_length += first_collision_time

        #@show typeof(which_obstacle_hit)

        if !isa(which_obstacle_hit, CellBoundary)

            push!(xs, p.x)
            push!(lattice_vectors, p.lattice_vector)
            push!(free_paths, free_path_length)
            free_path_length = 0.0

            num_obstacle_collisions += 1

        end
    end

    return xs, lattice_vectors, free_paths

end

@doc """Extract the particle position for given continuous times that are multiples of delta_t
        from a discrete-time trajectory.

        Use delta_t that are representable floating-point numbers
        e.g. 0.125 instead of 0.1""" ->
function continuous_time(xs, lattice_vectors, free_paths, delta_t)
    t = 0.0
    collision_positions = xs + lattice_vectors
    collision_times = cumsum([0; free_paths])

    t_final = collision_times[end]
    push!(collision_times, t_final)  # duplicate -- necessary?

    collision = 1
    time_counter = 1

    positions = [xs[1];]  # initial position at time 0
    times = [0.0;]

    next_output_time = delta_t

    time_left_for_step = delta_t

    # velocity:
    direction = collision_positions[collision+1] - collision_positions[collision]
    direction /= norm(direction)

    while next_output_time < t_final

        if next_output_time < collision_times[collision+1]

            flight_time = next_output_time - collision_times[collision]

            push!(positions, collision_positions[collision] + flight_time * direction)
            push!(times, next_output_time)
            time_counter += 1
            next_output_time = time_counter * delta_t
            # more accurate than adding delta_t if delta_t is not representable

        else

            collision += 1

            direction = collision_positions[collision+1] - collision_positions[collision]
            direction /= norm(direction)

        end
    end

    positions, times

end

max_stars_repo_name: dpsanders/BilliardModel

max_issues_repo_name: dpsanders/BilliardModel

max_forks_repo_name: dpsanders/BilliardModel


text:
using Printf
using LinearAlgebra
using Random

include("../common/PWGrid_v03.jl")
include("../common/ortho_gram_schmidt.jl")
include("../common/wrappers_fft.jl")
include("../common/gen_lattice_pwscf.jl")

include("op_K.jl")
include("op_Vpot.jl")
include("op_H.jl")
include("calc_rho.jl")
include("calc_Etot.jl")
include("Kprec.jl")
include("diag_lobpcg.jl")
include("diag_davidson.jl")
include("calc_grad.jl")
include("schsolve_Emin_sd.jl")
include("schsolve_Emin_cg.jl")
include("structure_factor.jl")

include("gen_rho.jl")
include("gen_dr.jl")
include("calc_ewald.jl")

function test_main( ecutwfc_Ry::Float64 )

    LatVecs = gen_lattice_sc(16.0)

    pw = PWGrid( 0.5*ecutwfc_Ry, LatVecs )

    Ω  = pw.Ω
    r  = pw.r
    G  = pw.gvec.G
    G2 = pw.gvec.G2
    Ns = pw.Ns
    Npoints = prod(Ns)
    Ngwx = pw.gvecw.Ngwx

    @printf("ecutwfc (Ha) = %f\n", ecutwfc_Ry*0.5)
    @printf("Ns   = (%d,%d,%d)\n", Ns[1], Ns[2], Ns[3])
    @printf("Ngwx = %d\n", Ngwx)

    actual = Npoints/Ngwx
    theor = 1/(4*pi*0.25^3/3)
    @printf("Compression: actual, theor: %f , %f\n", actual, theor)

    Xpos = reshape( [0.0, 0.0, 0.0], (3,1) )

    Sf = structure_factor( Xpos, G )

    E_nn = calc_ewald( pw, Xpos, Sf )

    Vg = zeros(ComplexF64,Npoints)
    prefactor = -4*pi/Ω
    for ig=2:Npoints
        Vg[ig] = prefactor/G2[ig]
    end
    Vpot = real( G_to_R(Ns, Vg .* Sf) ) * Npoints

    Nstates = 1
    srand(2222)
    psi = randn(Ngwx,Nstates) + im*randn(Ngwx,Nstates)
    psi = ortho_gram_schmidt(psi)
    #
    #evals, evecs = diag_lobpcg( pw, Vpot, psi, verbose=true, tol_avg=1e-10 )
    #Etot = sum(evals)

    evals, evecs = diag_davidson( pw, Vpot, psi, verbose=true, tol_avg=1e-10 )
    Etot = sum(evals)

    #psi, Etot = schsolve_Emin_sd( pw, Vpot, psi, NiterMax=10 )
    #psi, Etot = schsolve_Emin_cg( pw, Vpot, psi, NiterMax=1000 )
    #
    #Y = ortho_gram_schmidt(psi)
    #mu = Y' * op_H( pw, Vpot, Y )
    #evals, evecs = eigen(mu)
    #psi = Y*evecs

    for ist = 1:Nstates
        @printf("State %d, Energy = %18.10f\n", ist, real(evals[ist]))
    end

    @printf("E_nn    = %18.10f\n", E_nn)
    @printf("E total = %18.10f\n", E_nn + Etot)

end

@time test_main( 60.0 )

ext: jl


text:
"""
    read_metadata()

Read metadata stored in excel sheet (Serier til Avrenningskart 2016.xlsx)
"""
function read_metadata()

    file_name = joinpath(dirname(@__FILE__), "../raw/metadata.csv")

    df_meta = CSV.read(file_name; delim = ';', missingstring = "NULL")

    filter!(row -> !ismissing(row[:drainage_basin_key]), df_meta)

    regine_area = df_meta[:regine_area]
    main_no = df_meta[:main_no]

    df_meta[:regine_main] = ["$(regine_area[i]).$(main_no[i])" for i in eachindex(regine_area)]

    return df_meta

end

"""
    subset_metadata(df_meta, stat_list)

Select metadata using station selection
"""
function subset_metadata(df_meta, stat_list)

    df_meta = filter(row -> row[:regine_main] in stat_list, df_meta)

    return df_meta

end

"""
    read_dbk_ind()

Read file with information about drainage basin key and indices indicating
the location of the watersheds in the senorge grid (SeNorge.txt)
"""
function read_dbk_ind()

    file_name = joinpath(dirname(@__FILE__), "../raw/SeNorge.txt")

    lines = readlines(file_name)

    dbk_ind = OrderedDict{Float64, Array{Int64, 1}}()

    for line in lines
        str_splitted = split(line)
        key = parse(Float64, str_splitted[2])
        values = map(x -> parse(Int64,x) + 1, str_splitted[3:end])
        dbk_ind[key] = values
    end

    return dbk_ind

end

"""
    subset_dbk_ind(dict_in, df_meta)

Select dictionary linking senorge indices with drainage basin key using metadata table
"""
function subset_dbk_ind(dict_in, df_meta)

    dbk_sel = df_meta[:drainage_basin_key]

    dict_out = OrderedDict{Float64, Array{Int, 1}}()

    for dbk in dbk_sel

        dict_out[dbk] = dict_in[dbk]

    end

    return dict_out

end

"""
    watershed_info(dbk_ind)

Collect information (elevation, landuse, ...) about gridcells in a watershed.
"""
function watershed_info(dbk_ind, elev_breaks = 0:200:4000)

    # Initilize dictonary

    df_all = OrderedDict{Float64, DataFrame}()

    for (dbk, ind) in dbk_ind
        df_all[dbk] = DataFrame(senorge_ind = ind, area = 1.0)
    end

    # Add landuse and elevation information to dataframes

    raster_names = ["lus_unclassified", "lus_glacier", "lus_agriculture", "lus_bog",
    "lus_lake", "lus_forest", "lus_bare_mountain", "lus_urban", "elevation"]

    for raster in raster_names

        file_name = joinpath(dirname(@__FILE__), "../raw/$(raster).asc")

        data_raster = read_esri_raster(file_name)

        data_vec = raster2vec(data_raster)

        for (dbk, senorge_ind) in dbk_ind

            df_all[dbk][Meta.parse(raster)] = data_vec[senorge_ind]

        end

    end

    # Remove watershed if elevation data is missing

    for (dbk, df_wsh) in df_all

        if any(isnan.(df_wsh[:elevation]))
            delete!(df_all, dbk)
            info("Removed watershed with drainage basin key $dbk because elevation data is missing!")
        end

    end

    # Add indicies for elevation bands to dataframes

    for (dbk, df_wsh) in df_all

        elevation = df_wsh[:elevation]

        elev_ind = similar(elevation)

        elev_band = 1

        for ielev in 1:length(elev_breaks)-1

            ind = findall(elev_breaks[ielev] .<= elevation .< elev_breaks[ielev+1])

            if !isempty(ind)
                elev_ind[ind] .= elev_band
                elev_band += 1
            end

        end

        df_all[dbk][:elev_ind] = elev_ind

    end

    return df_all

end

"""
    aggregate_elevations!(metdata_elev, itime, df_geo, senorge_raw, scale_func; agg_func = mean)

Average data over elevation bands.
"""
function aggregate_elevations!(metdata_elev, itime, df_geo, senorge_raw, scale_func; agg_func = mean)

    for dbk in keys(df_geo)

        senorge_ind = convert(Array{Int64, 1}, df_geo[dbk][:senorge_ind])

        elev_ind = convert(Array{Int64, 1}, df_geo[dbk][:elev_ind])

        elev_iter = sort(unique(elev_ind))

        for ielev in elev_iter

            iagg = senorge_ind[elev_ind .== ielev]

            metdata_elev[dbk][itime, ielev] = agg_func(scale_func(senorge_raw[iagg]))

        end

    end

    return nothing

end

"""
    read_senorge_data(time_vec, met_var, df_geo; senorge = "V2.0", scale_func = x -> x)

Read senorge meteorological data.
"""
function read_senorge_data(time_vec, met_var, df_geo; senorge = "V2.0", scale_func = x -> x)

    # Allocate memory for outputs

    metdata_elev = OrderedDict{Float64, Array{Float64, 2}}()

    for (dbk, df_tmp) in df_geo

        nrows = length(time_vec)
        ncols = length(unique(df_tmp[:elev_ind]))

        metdata_elev[dbk] = zeros(nrows, ncols)

    end

    senorge_raw = zeros(Float64, 1550*1195)

    # Loop over time

    for itime = 1:length(time_vec)

        date_str1 = Dates.format(time_vec[itime], "yyyy")
        date_str2 = Dates.format(time_vec[itime], "yyyy_mm_dd")

        if senorge == "V2.0"

            path_met = "//hdata/grid/metdata/met_obs_v2.0"

            file_name = joinpath(path_met, "$met_var", "$date_str1", "$(met_var)_$date_str2.bil")

        end

        print("Processing data for $date_str2 \n")

        try

            read_bil!(senorge_raw, file_name)

            aggregate_elevations!(metdata_elev, itime, df_geo, senorge_raw, scale_func)

        catch

            error("No $met_var data available from $(time_vec[i]) and onwards")

        end

    end

    return metdata_elev, time_vec

end

"""
    write_met_data(save_folder, file_desc, time_vec, met_data, df_meta)

Write meteorological data to files
"""
function write_met_data(save_folder, file_desc, time_vec, met_data, df_meta)

    for row in eachrow(df_meta)

        try

            dbk = row[:drainage_basin_key]

            regine_main = row[:regine_main]

            file_path = joinpath(save_folder, "$(replace(regine_main, "." => "_"))_data")

            mkpath(file_path)

            file_name = joinpath(file_path, "$file_desc.txt")
            data = hcat(Dates.format.(time_vec, "yyyy-mm-dd HH:MM"), round.(met_data[dbk], digits = 2))

            f = open(file_name, "w")
            writedlm(f, data, ";", quotes=false)
            close(f)

        catch

            warn("Unable to write $file_desc data for station $stat_name")

        end

    end

end

"""
    write_update(save_folder, file_desc, time_old, met_old, time_new, met_new, df_meta)

Combine old and new senorge data and write to file
"""
function write_update(save_folder, file_desc, time_old, met_old, time_new, met_new, df_meta)

    for row in eachrow(df_meta)

        try

            @info "Write station $(row[:station_name])"

            regine_main = row[:regine_main]

            dbk = row[:drainage_basin_key]

            file_path = joinpath(save_folder, "$(replace(regine_main, "." => "_"))_data")

            mkpath(file_path)

            met_data = vcat(met_old[regine_main], met_new[dbk])
            time_vec = vcat(time_old, time_new)

            file_name = joinpath(file_path, "$file_desc.txt")
            data = hcat(Dates.format.(time_vec, "yyyy-mm-dd HH:MM"), round.(met_data, digits = 2))

            f = open(file_name, "w")
            writedlm(f, data, ";", quotes=false)
            close(f)

        catch

            @warn "Unable to write $file_desc data for station $regine_main"

        end

    end

end

"""
    metadata_elevbands(save_folder, df_geo, df_meta)

Compute average coverage of the different landuse classes, average elevation
and total area for the individual elevation zones
"""
function metadata_elevbands(save_folder, df_geo, df_meta)

    regine_main = convert(Array, df_meta[:regine_main])

    for (dbk, df_tmp) in df_geo

        # Aggregate data for elevation bands

        df_agg = aggregate(df_tmp, :elev_ind, [mean, sum])

        colnames = [Symbol(replace(string(name), "Statistics." => "")) for name in names(df_agg)]

        names!(df_agg, colnames)

        keep = [:elev_ind, :lus_unclassified_mean, :lus_glacier_mean, :lus_agriculture_mean,
                :lus_bog_mean, :lus_lake_mean, :lus_forest_mean, :lus_bare_mountain_mean,
                :lus_urban_mean, :elevation_mean, :area_sum]

        df_save = df_agg[keep]

        # Save data to files

        irow = findall(df_meta[:drainage_basin_key] .== dbk)

        regine_main = df_meta[:regine_main][irow][1]

        file_path = joinpath(save_folder, "$(replace(regine_main, "." => "_"))_data")

        mkpath(file_path)

        file_name = joinpath(file_path, "metadata.txt")

        CSV.write(file_name, df_save; delim = ';')

    end

end

"""
    process_runoff_data(stat_list, time_vec, df_meta, path_runoff, save_folder)

Process runoff data obtained using shell scripts.
"""
function process_runoff_data(stat_list, time_vec, df_meta, path_runoff, save_folder)

    # File names for all available stations

    regine_area = convert(Array{Int64}, df_meta[:regine_area])
    main_no = convert(Array{Int64}, df_meta[:main_no])
    point_no = convert(Array{Int64}, df_meta[:point_no])
    param_key = convert(Array{Int64}, df_meta[:param_key])
    version_no_end = convert(Array{Int64}, df_meta[:version_no_end])

    stat_all = ["$(regine_area[i]).$(main_no[i])" for i in eachindex(regine_area)]

    # Catchment area

    area_total = df_meta[:area_total]

    # Loop over stations in input list

    for stat_name in stat_list

        @info "Save runoff data for station: $(stat_name)"

        # Load data for one station

        istat = findall(stat_all .== stat_name)
        istat = istat[1]

        file_name = "$(regine_area[istat]).$(main_no[istat]).$(point_no[istat]).$(param_key[istat]).$(version_no_end[istat])"

        data = readdlm(joinpath(path_runoff, file_name), skipstart = 2)

        # Match desired time period using data frames and joins

        df_runoff = DataFrame(Time = map(DateTime, data[:,1]), Runoff = data[:,2])

        df_ref = DataFrame(Time = time_vec)

        df_final = join(df_ref, df_runoff, on = :Time, kind = :left)

	    sort!(df_final, [:Time])

        # Save data to file

        try

            time_save = df_final[:Time]

            runoff_save = df_final[:Runoff]

            area = tryparse(Float64, replace(area_total[istat], "," => "."))

            runoff_save = (runoff_save * 86400 * 1000) / (area * 1e6)  # Convert from m3/s to mm/day

            runoff_save = convert(Array{Union{Missing, Float64}, 1}, runoff_save)

            runoff_save[ismissing.(runoff_save)] .= -999.0

            data = hcat(Dates.format.(time_save, "yyyy-mm-dd HH:MM"), round.(runoff_save; digits = 2))

            file_path = joinpath(save_folder, "$(replace(stat_name, "." => "_"))_data")
            file_name = joinpath(file_path, "runoff.txt")

            mkpath(file_path)

            f = open(file_name, "w")
            writedlm(f, data, ";", quotes=false)
            close(f)

        catch

            @warn "Unable to write runoff data for station $stat_name"

        end

    end

end

"""
     read_old_data(save_folder, met_var; crop = 10)

Read existing data.
"""
function read_old_data(save_folder, met_var; crop = 10)

    stat_list = readdir(save_folder)
    stat_list = [replace(x, "_data" => "") for x in stat_list]
    stat_list = [replace(x, "_" => ".") for x in stat_list]

    data = OrderedDict()

    time = []

    for stat_name in stat_list

        file_path = joinpath(save_folder, "$(replace(stat_name, "." => "_"))_data", "$met_var.txt")

        str  = readline(file_path)
        nsep = length(collect((m.match for m = eachmatch(r";", str))))
        tmp  = CSV.read(file_path, delim = ";", header = false,
                        dateformat="yyyy-mm-dd HH:MM", allowmissing=:none,
                        types = vcat(DateTime, repeat([Float64], nsep)))

        time = tmp[1:end-crop, 1]
        tmp  = convert(Array{Float64}, tmp[1:end-crop, 2:end])

        data[stat_name] = tmp

    end

    return data, time

end

"""
    find_last_senorge(senorge_version)

Find last available senorge data.
"""
function find_last_senorge(senorge_version)

    time_start = floor(now(), Dates.Day)
    time_stop = time_start + Dates.Day(10)

    time_final = []

    for time_check = time_start:Day(1):time_stop

        date_str1 = Dates.format(time_check, "yyyy")
        date_str2 = Dates.format(time_check, "yyyy_mm_dd")

        if senorge_version == "V2.0"
            path_met = "//hdata/grid/metdata/met_obs_v2.0"
            file_name = joinpath(path_met, "rr", "$date_str1", "rr_$date_str2.bil")
        end

        if isfile(file_name)
            time_final = time_check
        else
            break
        end

    end

    return time_final

end

max_forks_repo_name: jmgnve/NveData


text:
"""
    AdiabaticModels

All models defined within this module have only a single electronic state
and return potentials as scalars and derivatives as simple arrays.

The central abstract type is the [`AdiabaticModel`](@ref), which all models should subtype.
"""
module AdiabaticModels

using ..NQCModels: NQCModels
using Requires: Requires
using Parameters: Parameters

using Unitful: @u_str, ustrip
using UnitfulAtomic: austrip, auconvert

"""
    AdiabaticModel <: Model

`AdiabaticModel`s represent the potentials from classical molecular dynamics
where the potential is a function of the position.

# Implementation

`AdiabaticModel`s should implement:
* `potential(model, R)`
* `derivative!(model, D, R)` (this is the derivative of the potential energy with respect to the positions)
* `ndofs(model)` (these are the degrees of freedom)

# Example

This example creates a 2 dimensional adiabatic model `MyModel`.
We implement the 3 compulsory functions then evaluate the potential.
Here, the argument `R` is an `AbstractMatrix` since this is a 2D model
that can accept multiple atoms.

```jldoctest
struct MyModel{P} <: NQCModels.AdiabaticModels.AdiabaticModel
    param::P
end

NQCModels.ndofs(::MyModel) = 2

NQCModels.potential(model::MyModel, R::AbstractMatrix) = model.param*sum(R.^2)
NQCModels.derivative!(model::MyModel, D, R::AbstractMatrix) = D .= model.param*2R

model = MyModel(10)

NQCModels.potential(model, [1 2; 3 4])

# output

300

```
"""
abstract type AdiabaticModel <: NQCModels.Model end

NQCModels.nstates(::AdiabaticModel) = 1

NQCModels.zero_derivative(::AdiabaticModel, R) = zero(R)

include("free.jl")
export Free
include("morse.jl")
export Morse
include("harmonic.jl")
export Harmonic
include("diatomic_harmonic.jl")
export DiatomicHarmonic
include("darling_holloway_elbow.jl")
export DarlingHollowayElbow
include("ase_interface.jl")
export AdiabaticASEModel
include("averaged_potential.jl")
export AveragedPotential

function __init__()
    Requires.@require JuLIP="945c410c-986d-556a-acb1-167a618e0462" @eval include("julip.jl")
end

end # module


text:
@testset "pooling" begin

    # Test for agreement between CPU NNlib and CuDNN versions, across a variety of kwargs
    for num_spatial_dims in (1, 2, 3)
        # Initialize data we'll run our tests over
        C_in = 3
        batch_size = 1
        x = rand(Float64, fill(8, num_spatial_dims)..., C_in, batch_size)

        # Test that pooling is equivalent across GPU/CPU
        pdims = PoolDims(x, 2)
        y = maxpool(x, pdims)
        dy = ones(size(y))
        gputest(x -> maxpool(x, pdims), x)
        gputest((dy, y, x) -> ∇maxpool(dy, y, x, pdims), dy, y, x, checkgrad=false)
        gputest(x -> maxpool(x, pdims), x)
        gputest((dy, y, x) -> ∇maxpool(dy, y, x, pdims), dy, y, x, checkgrad=false)

        # Test the compatibility shims for pooling
        cx,cy,cdy = CuArray{Float32}.((x,y,dy))
        win,pad=2,1
        maxpool!(similar(cy), cx, win; pad=pad, stride=win) ≈ maxpool!(similar(cy), cx, PoolDims(cx, win; padding=pad, stride=win))
        meanpool!(similar(cy), cx, win; pad=pad, stride=win) ≈ meanpool!(similar(cy), cx, PoolDims(cx, win; padding=pad, stride=win))
    end
end


text:
"""
    Metaheuristic

Abstract type for metaheuristic algorithms, this makes a clear
distinction between different classifications of metaheuristic algorithms.
"""
abstract type Metaheuristic end

"""
    PopulationBase <: Metaheuristic

Type for population-based algorithms that employ [`Population`](@ref), i.e.
subroutines that _mutate_ an array of possible candidates in-place.
An example of this type is [`PSO`](@ref).
"""
abstract type PopulationBase <: Metaheuristic end

"""
    TrajectoryBase <: Metaheuristic

Type for trajectory-based algorithms that only uses specific heuristics to guide the
solution to a global minimum.
An example of this type is [`SimulatedAnnealing`](@ref).
"""
abstract type TrajectoryBase <: Metaheuristic end

"""
    OptimizationResults{T, U}

Type that formats the output of [`Metaheuristic`](@ref) to get better information
from it.

# Fields
- `x::T`: Stores the _solution_ array from the solver, i.e. the solution that minimizes
    the cost function.
- `min::U`: Stores the value obtained from evaluating the cost function with
    `x`, i.e. the minima found.
- `impl::AbstractString`: Stores the name of the `Metaheuristic` used, i.e. the name or identifier of the
    optimization algorithm.
- `iterations::Integer`: Stores the number of maximum iterations that the solver was run.
"""
mutable struct OptimizationResults{T,U}
    x::T
    min::U
    impl::AbstractString
    iterations::Integer
end

function Base.show(io::IO, r::OptimizationResults)
    println("Results from Optimization")
    @printf(io, "\tAlgorithm: %s\n", r.impl)
    @printf(io, "\tSolution: [%s]\n", join(r.x, ", "))
    @printf(io, "\tMinimum: %.4f\n", r.min)
    @printf(io, "\tMaximum iterations: %d\n", r.iterations)
end


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

struct BatchIterator
    bytes::Vector{UInt8}
    startpos::Int
end

"""
    Arrow.Stream(io::IO; convert::Bool=true)
    Arrow.Stream(file::String; convert::Bool=true)
    Arrow.Stream(bytes::Vector{UInt8}, pos=1, len=nothing; convert::Bool=true)

Start reading an arrow formatted table, from:
 * `io`, bytes will be read all at once via `read(io)`
 * `file`, bytes will be read via `Mmap.mmap(file)`
 * `bytes`, a byte vector directly, optionally allowing specifying the starting byte position `pos` and `len`

Reads the initial schema message from the arrow stream/file, then returns an `Arrow.Stream` object
which will iterate over record batch messages, producing an [`Arrow.Table`](@ref) on each iteration.

By iterating [`Arrow.Table`](@ref), `Arrow.Stream` satisfies the `Tables.partitions` interface, and as such can
be passed to Tables.jl-compatible sink functions.

This allows iterating over extremely large "arrow tables" in chunks represented as record batches.

Supports the `convert` keyword argument which controls whether certain arrow primitive types will be
lazily converted to more friendly Julia defaults; by default, `convert=true`.
"""
struct Stream
    batchiterator::BatchIterator
    pos::Int
    names::Vector{Symbol}
    schema::Meta.Schema
    dictencodings::Dict{Int64, DictEncoding} # dictionary id => DictEncoding
    dictencoded::Dict{Int64, Meta.Field} # dictionary id => field
    convert::Bool
end

Tables.partitions(x::Stream) = x

Stream(io::IO, pos::Integer=1, len=nothing; convert::Bool=true) = Stream(Base.read(io), pos, len; convert=convert)
Stream(str::String, pos::Integer=1, len=nothing; convert::Bool=true) = isfile(str) ? Stream(Mmap.mmap(str), pos, len; convert=convert) :
    throw(ArgumentError("$str is not a valid arrow file"))

# will detect whether we're reading a Stream from a file or stream
function Stream(bytes::Vector{UInt8}, off::Integer=1, tlen::Union{Integer, Nothing}=nothing; convert::Bool=true)
    len = something(tlen, length(bytes))
    if len > 24 &&
        _startswith(bytes, off, FILE_FORMAT_MAGIC_BYTES) &&
        _endswith(bytes, off + len - 1, FILE_FORMAT_MAGIC_BYTES)
        off += 8 # skip past magic bytes + padding
    end
    dictencodings = Dict{Int64, DictEncoding}() # dictionary id => DictEncoding
    dictencoded = Dict{Int64, Meta.Field}() # dictionary id => field
    batchiterator = BatchIterator(bytes, off)
    state = iterate(batchiterator)
    state === nothing && throw(ArgumentError("no arrow ipc messages found in provided input"))
    batch, (pos, id) = state
    schema = batch.msg.header
    schema isa Meta.Schema || throw(ArgumentError("first arrow ipc message MUST be a schema message"))
    # assert endianness?
    # store custom_metadata?
    names = Symbol[]
    for (i, field) in enumerate(schema.fields)
        push!(names, Symbol(field.name))
        # recursively find any dictionaries for any fields
        getdictionaries!(dictencoded, field)
        @debug 1 "parsed column from schema: field = $field"
    end
    return Stream(batchiterator, pos, names, schema, dictencodings, dictencoded, convert)
end

function Base.iterate(x::Stream, (pos, id)=(x.pos, 1))
    columns = AbstractVector[]
    while true
        state = iterate(x.batchiterator, (pos, id))
        state === nothing && return nothing
        batch, (pos, id) = state
        header = batch.msg.header
        if header isa Meta.DictionaryBatch
            id = header.id
            recordbatch = header.data
            @debug 1 "parsing dictionary batch message: id = $id, compression = $(recordbatch.compression)"
            if haskey(x.dictencodings, id) && header.isDelta
                # delta
                field = x.dictencoded[id]
                values, _, _ = build(field, field.type, batch, recordbatch, x.dictencodings, Int64(1), Int64(1), x.convert)
                dictencoding = x.dictencodings[id]
                append!(dictencoding.data, values)
                continue
            end
            # new dictencoding or replace
            field = x.dictencoded[id]
            values, _, _ = build(field, field.type, batch, recordbatch, x.dictencodings, Int64(1), Int64(1), x.convert)
            A = ChainedVector([values])
            x.dictencodings[id] = DictEncoding{eltype(A), typeof(A)}(id, A, field.dictionary.isOrdered)
            @debug 1 "parsed dictionary batch message: id=$id, data=$values\n"
        elseif header isa Meta.RecordBatch
            @debug 1 "parsing record batch message: compression = $(header.compression)"
            for vec in VectorIterator(x.schema, batch, x.dictencodings, x.convert)
                push!(columns, vec)
            end
            break
        else
            throw(ArgumentError("unsupported arrow message type: $(typeof(header))"))
        end
    end
    lookup = Dict{Symbol, AbstractVector}()
    types = Type[]
    for (nm, col) in zip(x.names, columns)
        lookup[nm] = col
        push!(types, eltype(col))
    end
    return Table(x.names, types, columns, lookup, Ref(x.schema)), (pos, id)
end

"""
    Arrow.Table(io::IO; convert::Bool=true)
    Arrow.Table(file::String; convert::Bool=true)
    Arrow.Table(bytes::Vector{UInt8}, pos=1, len=nothing; convert::Bool=true)

Read an arrow formatted table, from:
 * `io`, bytes will be read all at once via `read(io)`
 * `file`, bytes will be read via `Mmap.mmap(file)`
 * `bytes`, a byte vector directly, optionally allowing specifying the starting byte position `pos` and `len`

Returns a `Arrow.Table` object that allows column access via `table.col1`, `table[:col1]`, or `table[1]`.

NOTE: the columns in an `Arrow.Table` are views into the original arrow memory, and hence are not easily
modifiable (with e.g. `push!`, `append!`, etc.). To mutate arrow columns, call `copy(x)` to materialize
the arrow data as a normal Julia array.

`Arrow.Table` also satisfies the [Tables.jl](https://github.com/JuliaData/Tables.jl) interface, and so can easily be materialied via any supporting
sink function: e.g. `DataFrame(Arrow.Table(file))`, `SQLite.load!(db, "table", Arrow.Table(file))`, etc.

Supports the `convert` keyword argument which controls whether certain arrow primitive types will be
lazily converted to more friendly Julia defaults; by default, `convert=true`.
"""
struct Table <: Tables.AbstractColumns
    names::Vector{Symbol}
    types::Vector{Type}
    columns::Vector{AbstractVector}
    lookup::Dict{Symbol, AbstractVector}
    schema::Ref{Meta.Schema}
end

Table() = Table(Symbol[], Type[], AbstractVector[], Dict{Symbol, AbstractVector}(), Ref{Meta.Schema}())

names(t::Table) = getfield(t, :names)
types(t::Table) = getfield(t, :types)
columns(t::Table) = getfield(t, :columns)
lookup(t::Table) = getfield(t, :lookup)
schema(t::Table) = getfield(t, :schema)

Tables.istable(::Table) = true
Tables.columnaccess(::Table) = true
Tables.columns(t::Table) = Tables.CopiedColumns(t)
Tables.schema(t::Table) = Tables.Schema(names(t), types(t))
Tables.columnnames(t::Table) = names(t)
Tables.getcolumn(t::Table, i::Int) = columns(t)[i]
Tables.getcolumn(t::Table, nm::Symbol) = lookup(t)[nm]

# high-level user API functions
Table(io::IO, pos::Integer=1, len=nothing; convert::Bool=true) = Table(Base.read(io), pos, len; convert=convert)
Table(str::String, pos::Integer=1, len=nothing; convert::Bool=true) = isfile(str) ? Table(Mmap.mmap(str), pos, len; convert=convert) :
    throw(ArgumentError("$str is not a valid arrow file"))

# will detect whether we're reading a Table from a file or stream
function Table(bytes::Vector{UInt8}, off::Integer=1, tlen::Union{Integer, Nothing}=nothing; convert::Bool=true)
    len = something(tlen, length(bytes))
    if len > 24 &&
        _startswith(bytes, off, FILE_FORMAT_MAGIC_BYTES) &&
        _endswith(bytes, off + len - 1, FILE_FORMAT_MAGIC_BYTES)
        off += 8 # skip past magic bytes + padding
    end
    t = Table()
    sch = nothing
    dictencodings = Dict{Int64, DictEncoding}() # dictionary id => DictEncoding
    dictencoded = Dict{Int64, Meta.Field}() # dictionary id => field
    tsks = Channel{Task}(Inf)
    tsk = Threads.@spawn begin
        i = 1
        for tsk in tsks
            cols = fetch(tsk)
            if i == 1
                foreach(x -> push!(columns(t), x), cols)
            elseif i == 2
                foreach(1:length(cols)) do i
                    columns(t)[i] = ChainedVector([columns(t)[i], cols[i]])
                end
            else
                foreach(1:length(cols)) do i
                    append!(columns(t)[i], cols[i])
                end
            end
            i += 1
        end
    end
    for batch in BatchIterator(bytes, off)
        # store custom_metadata of batch.msg?
        header = batch.msg.header
        if header isa Meta.Schema
            @debug 1 "parsing schema message"
            # assert endianness?
            # store custom_metadata?
            for (i, field) in enumerate(header.fields)
                push!(names(t), Symbol(field.name))
                # recursively find any dictionaries for any fields
                getdictionaries!(dictencoded, field)
                @debug 1 "parsed column from schema: field = $field"
            end
            sch = header
            schema(t)[] = sch
        elseif header isa Meta.DictionaryBatch
            id = header.id
            recordbatch = header.data
            @debug 1 "parsing dictionary batch message: id = $id, compression = $(recordbatch.compression)"
            if haskey(dictencodings, id) && header.isDelta
                # delta
                field = dictencoded[id]
                values, _, _ = build(field, field.type, batch, recordbatch, dictencodings, Int64(1), Int64(1), convert)
                dictencoding = dictencodings[id]
                append!(dictencoding.data, values)
                continue
            end
            # new dictencoding or replace
            field = dictencoded[id]
            values, _, _ = build(field, field.type, batch, recordbatch, dictencodings, Int64(1), Int64(1), convert)
            A = ChainedVector([values])
            dictencodings[id] = DictEncoding{eltype(A), typeof(A)}(id, A, field.dictionary.isOrdered)
            @debug 1 "parsed dictionary batch message: id=$id, data=$values\n"
        elseif header isa Meta.RecordBatch
            @debug 1 "parsing record batch message: compression = $(header.compression)"
            put!(tsks, Threads.@spawn begin
                collect(VectorIterator(sch, batch, dictencodings, convert))
            end)
        else
            throw(ArgumentError("unsupported arrow message type: $(typeof(header))"))
        end
    end
    close(tsks)
    wait(tsk)
    lu = lookup(t)
    ty = types(t)
    for (nm, col) in zip(names(t), columns(t))
        lu[nm] = col
        push!(ty, eltype(col))
    end
    meta = sch !== nothing ? sch.custom_metadata : nothing
    if meta !== nothing
        setmetadata!(t, Dict(String(kv.key) => String(kv.value) for kv in meta))
    end
    return t
end

function getdictionaries!(dictencoded, field)
    d = field.dictionary
    if d !== nothing
        dictencoded[d.id] = field
    end
    for child in field.children
        getdictionaries!(dictencoded, child)
    end
    return
end

struct Batch
    msg::Meta.Message
    bytes::Vector{UInt8}
    pos::Int
    id::Int
end

function Base.iterate(x::BatchIterator, (pos, id)=(x.startpos, 0))
    @debug 1 "checking for next arrow message: pos = $pos"
    if pos + 3 > length(x.bytes)
        @debug 1 "not enough bytes left for another batch message"
        return nothing
    end
    if readbuffer(x.bytes, pos, UInt32) != CONTINUATION_INDICATOR_BYTES
        @debug 1 "didn't find continuation byte to keep parsing messages: $(readbuffer(x.bytes, pos, UInt32))"
        return nothing
    end
    pos += 4
    if pos + 3 > length(x.bytes)
        @debug 1 "not enough bytes left to read length of another batch message"
        return nothing
    end
    msglen = readbuffer(x.bytes, pos, Int32)
    if msglen == 0
        @debug 1 "message has 0 length; terminating message parsing"
        return nothing
    end
    pos += 4
    msg = FlatBuffers.getrootas(Meta.Message, x.bytes, pos-1)
    pos += msglen
    # pos now points to message body
    @debug 1 "parsing message: pos = $pos, msglen = $msglen, bodyLength = $(msg.bodyLength)"
    return Batch(msg, x.bytes, pos, id), (pos + msg.bodyLength, id + 1)
end

struct VectorIterator
    schema::Meta.Schema
    batch::Batch # batch.msg.header MUST BE RecordBatch
    dictencodings::Dict{Int64, DictEncoding}
    convert::Bool
end

buildmetadata(f::Meta.Field) = buildmetadata(f.custom_metadata)
buildmetadata(meta) = Dict(String(kv.key) => String(kv.value) for kv in meta)
buildmetadata(::Nothing) = nothing

function Base.iterate(x::VectorIterator, (columnidx, nodeidx, bufferidx)=(Int64(1), Int64(1), Int64(1)))
    columnidx > length(x.schema.fields) && return nothing
    field = x.schema.fields[columnidx]
    @debug 2 "building top-level column: field = $(field), columnidx = $columnidx, nodeidx = $nodeidx, bufferidx = $bufferidx"
    A, nodeidx, bufferidx = build(field, x.batch, x.batch.msg.header, x.dictencodings, nodeidx, bufferidx, x.convert)
    @debug 2 "built top-level column: A = $(typeof(A)), columnidx = $columnidx, nodeidx = $nodeidx, bufferidx = $bufferidx"
    @debug 3 A
    return A, (columnidx + 1, nodeidx, bufferidx)
end

Base.length(x::VectorIterator) = length(x.schema.fields)

const ListTypes = Union{Meta.Utf8, Meta.LargeUtf8, Meta.Binary, Meta.LargeBinary, Meta.List, Meta.LargeList}
const LargeLists = Union{Meta.LargeUtf8, Meta.LargeBinary, Meta.LargeList}

function build(field::Meta.Field, batch, rb, de, nodeidx, bufferidx, convert)
    d = field.dictionary
    if d !== nothing
        validity = buildbitmap(batch, rb, nodeidx, bufferidx)
        bufferidx += 1
        buffer = rb.buffers[bufferidx]
        S = d.indexType === nothing ? Int32 : juliaeltype(field, d.indexType, false)
        bytes, indices = reinterp(S, batch, buffer, rb.compression)
        encoding = de[d.id]
        A = DictEncoded(bytes, validity, indices, encoding, buildmetadata(field.custom_metadata))
        nodeidx += 1
        bufferidx += 1
    else
        A, nodeidx, bufferidx = build(field, field.type, batch, rb, de, nodeidx, bufferidx, convert)
    end
    return A, nodeidx, bufferidx
end

function buildbitmap(batch, rb, nodeidx, bufferidx)
    buffer = rb.buffers[bufferidx]
    voff = batch.pos + buffer.offset
    node = rb.nodes[nodeidx]
    if rb.compression === nothing
        return ValidityBitmap(batch.bytes, voff, node.length, node.null_count)
    else
        # compressed
        ptr = pointer(batch.bytes, voff)
        _, decodedbytes = uncompress(ptr, buffer, rb.compression)
        return ValidityBitmap(decodedbytes, 1, node.length, node.null_count)
    end
end

function uncompress(ptr::Ptr{UInt8}, buffer, compression)
    if buffer.length == 0
        return 0, UInt8[]
    end
    len = unsafe_load(convert(Ptr{Int64}, ptr))
    ptr += 8 # skip past uncompressed length as Int64
    encodedbytes = unsafe_wrap(Array, ptr, buffer.length - 8)
    if compression.codec === Meta.CompressionType.LZ4_FRAME
        decodedbytes = transcode(LZ4FrameDecompressor, encodedbytes)
    elseif compression.codec === Meta.CompressionType.ZSTD
        decodedbytes = transcode(ZstdDecompressor, encodedbytes)
    else
        error("unsupported compression type when reading arrow buffers: $(typeof(compression.codec))")
    end
    return len, decodedbytes
end

function reinterp(::Type{T}, batch, buf, compression) where {T}
    ptr = pointer(batch.bytes, batch.pos + buf.offset)
    if compression === nothing
        return batch.bytes, unsafe_wrap(Array, convert(Ptr{T}, ptr), div(buf.length, sizeof(T)))
    else
        # compressed
        len, decodedbytes = uncompress(ptr, buf, compression)
        return decodedbytes, unsafe_wrap(Array, convert(Ptr{T}, pointer(decodedbytes)), div(len, sizeof(T)))
    end
end

function build(f::Meta.Field, L::ListTypes, batch, rb, de, nodeidx, bufferidx, convert)
    @debug 2 "building array: L = $L"
    validity = buildbitmap(batch, rb, nodeidx, bufferidx)
    bufferidx += 1
    buffer = rb.buffers[bufferidx]
    ooff = batch.pos + buffer.offset
    OT = L isa LargeLists ? Int64 : Int32
    bytes, offs = reinterp(OT, batch, buffer, rb.compression)
    offsets = Offsets(bytes, offs)
    bufferidx += 1
    len = rb.nodes[nodeidx].length
    nodeidx += 1
    if L isa Meta.Utf8 || L isa Meta.LargeUtf8 || L isa Meta.Binary || L isa Meta.LargeBinary
        buffer = rb.buffers[bufferidx]
        bytes, A = reinterp(UInt8, batch, buffer, rb.compression)
        bufferidx += 1
    else
        bytes = UInt8[]
        A, nodeidx, bufferidx = build(f.children[1], batch, rb, de, nodeidx, bufferidx, convert)
    end
    meta = buildmetadata(f.custom_metadata)
    T = juliaeltype(f, meta, convert)
    return List{T, OT, typeof(A)}(bytes, validity, offsets, A, len, meta), nodeidx, bufferidx
end

function build(f::Meta.Field, L::Union{Meta.FixedSizeBinary, Meta.FixedSizeList}, batch, rb, de, nodeidx, bufferidx, convert)
    @debug 2 "building array: L = $L"
    validity = buildbitmap(batch, rb, nodeidx, bufferidx)
    bufferidx += 1
    len = rb.nodes[nodeidx].length
    nodeidx += 1
    if L isa Meta.FixedSizeBinary
        buffer = rb.buffers[bufferidx]
        bytes, A = reinterp(UInt8, batch, buffer, rb.compression)
        bufferidx += 1
    else
        bytes = UInt8[]
        A, nodeidx, bufferidx = build(f.children[1], batch, rb, de, nodeidx, bufferidx, convert)
    end
    meta = buildmetadata(f.custom_metadata)
    T = juliaeltype(f, meta, convert)
    return FixedSizeList{T, typeof(A)}(bytes, validity, A, len, meta), nodeidx, bufferidx
end

function build(f::Meta.Field, L::Meta.Map, batch, rb, de, nodeidx, bufferidx, convert)
    @debug 2 "building array: L = $L"
    validity = buildbitmap(batch, rb, nodeidx, bufferidx)
    bufferidx += 1
    buffer = rb.buffers[bufferidx]
    ooff = batch.pos + buffer.offset
    OT = Int32
    bytes, offs = reinterp(OT, batch, buffer, rb.compression)
    offsets = Offsets(bytes, offs)
    bufferidx += 1
    len = rb.nodes[nodeidx].length
    nodeidx += 1
    A, nodeidx, bufferidx = build(f.children[1], batch, rb, de, nodeidx, bufferidx, convert)
    meta = buildmetadata(f.custom_metadata)
    T = juliaeltype(f, meta, convert)
    return Map{T, OT, typeof(A)}(validity, offsets, A, len, meta), nodeidx, bufferidx
end

function build(f::Meta.Field, L::Meta.Struct, batch, rb, de, nodeidx, bufferidx, convert)
    @debug 2 "building array: L = $L"
    validity = buildbitmap(batch, rb, nodeidx, bufferidx)
    bufferidx += 1
    len = rb.nodes[nodeidx].length
    vecs = []
    nodeidx += 1
    for child in f.children
        A, nodeidx, bufferidx = build(child, batch, rb, de, nodeidx, bufferidx, convert)
        push!(vecs, A)
    end
    data = Tuple(vecs)
    meta = buildmetadata(f.custom_metadata)
    T = juliaeltype(f, meta, convert)
    return Struct{T, typeof(data)}(validity, data, len, meta), nodeidx, bufferidx
end

function build(f::Meta.Field, L::Meta.Union, batch, rb, de, nodeidx, bufferidx, convert)
    @debug 2 "building array: L = $L"
    buffer = rb.buffers[bufferidx]
    bytes, typeIds = reinterp(UInt8, batch, buffer, rb.compression)
    bufferidx += 1
    if L.mode == Meta.UnionMode.Dense
        buffer = rb.buffers[bufferidx]
        bytes2, offsets = reinterp(Int32, batch, buffer, rb.compression)
        bufferidx += 1
    end
    vecs = []
    nodeidx += 1
    for child in f.children
        A, nodeidx, bufferidx = build(child, batch, rb, de, nodeidx, bufferidx, convert)
        push!(vecs, A)
    end
    data = Tuple(vecs)
    meta = buildmetadata(f.custom_metadata)
    T = juliaeltype(f, meta, convert)
    if L.mode == Meta.UnionMode.Dense
        B = DenseUnion{T, typeof(data)}(bytes, bytes2, typeIds, offsets, data, meta)
    else
        B = SparseUnion{T, typeof(data)}(bytes, typeIds, data, meta)
    end
    return B, nodeidx, bufferidx
end

function build(f::Meta.Field, L::Meta.Null, batch, rb, de, nodeidx, bufferidx, convert)
    @debug 2 "building array: L = $L"
    return MissingVector(rb.nodes[nodeidx].length), nodeidx + 1, bufferidx
end

# primitives
function build(f::Meta.Field, ::L, batch, rb, de, nodeidx, bufferidx, convert) where {L}
    @debug 2 "building array: L = $L"
    validity = buildbitmap(batch, rb, nodeidx, bufferidx)
    bufferidx += 1
    buffer = rb.buffers[bufferidx]
    meta = buildmetadata(f.custom_metadata)
    # get storage type (non-converted)
    T = juliaeltype(f, nothing, false)
    @debug 2 "storage type for primitive: T = $T"
    bytes, A = reinterp(Base.nonmissingtype(T), batch, buffer, rb.compression)
    len = rb.nodes[nodeidx].length
    T = juliaeltype(f, meta, convert)
    @debug 2 "final julia type for primitive: T = $T"
    return Primitive(T, bytes, validity, A, len, meta), nodeidx + 1, bufferidx + 1
end

function build(f::Meta.Field, L::Meta.Bool, batch, rb, de, nodeidx, bufferidx, convert)
    @debug 2 "building array: L = $L"
    validity = buildbitmap(batch, rb, nodeidx, bufferidx)
    bufferidx += 1
    buffer = rb.buffers[bufferidx]
    meta = buildmetadata(f.custom_metadata)
    # get storage type (non-converted)
    T = juliaeltype(f, nothing, false)
    @debug 2 "storage type for primitive: T = $T"
    buffer = rb.buffers[bufferidx]
    voff = batch.pos + buffer.offset
    node = rb.nodes[nodeidx]
    if rb.compression === nothing
        decodedbytes = batch.bytes
        pos = voff
        # return ValidityBitmap(batch.bytes, voff, node.length, node.null_count)
    else
        # compressed
        ptr = pointer(batch.bytes, voff)
        _, decodedbytes = uncompress(ptr, buffer, rb.compression)
        pos = 1
        # return ValidityBitmap(decodedbytes, 1, node.length, node.null_count)
    end
    len = rb.nodes[nodeidx].length
    T = juliaeltype(f, meta, convert)
    return BoolVector{T}(decodedbytes, pos, validity, len, meta), nodeidx + 1, bufferidx + 1
end

max_stars_repo_name: sweb/arrow

max_issues_repo_name: sweb/arrow

max_forks_repo_name: sweb/arrow


text:
@testset "BicycleGame" begin
    # Test BicycleGame
    p = 3
    model = BicycleGame(p=p,lr=1.0,lf=2.0)
    @test model.lr == 1.0
    @test model.lf == 2.0
    @test typeof(model) <: BicycleGame{model.n, model.m, model.p, model.lf, model.lr}
    @test model.n == 4*p
    @test model.m == 2p
    @test model.p == p
    @test model.ni == 4*ones(Int,p)
    @test model.mi == 2*ones(Int,p)
    @test model.pu[1] == SVector{2,Int}((1,4))
    @test model.pu[2] == SVector{2,Int}((2,5))
    @test model.pu[3] == SVector{2,Int}((3,6))

    @test model.px[1] == SVector{2,Int}((1,4))
    @test model.px[2] == SVector{2,Int}((2,5))
    @test model.px[3] == SVector{2,Int}((3,6))

    @test model.pz[1] == SVector{4,Int}((1,4,7,10))
    @test model.pz[2] == SVector{4,Int}((2,5,8,11))
    @test model.pz[3] == SVector{4,Int}((3,6,9,12))

    @test size(model) == (4*p, 2p, model.pu, p)

    T = Float64
    x = ones(SVector{model.n,T})
    u = 100*ones(SVector{model.m,T})
    @test (@ballocated dynamics($model, $x, $u)) == 0
    # @btime dynamics($model, $x, $u)

end


text:
"""
A script attempting to make oscilatory assemblies in an AssemblyCalculus.Brain

Notes:

First attempt

1. Create a0 in area 3
2. Use four random inputs and a0 to create assemblies a1, a3 in area 1 and a2, a4
in area 2
3. Reciprocal project w/o a0
      a1 <-> a2
      a3 <-> a2
      a3 <-> a4
      a1 <-> a4

4. Run with zero input to area1 and area2 and a0 for 100 timesteps

In most of the simulations I tried, the final result was not oscilation, but
rather a fixed point that encompassed a number of neurons from both possible
area assemblies, but not all.

When I made β very small, and the edge probability very high, I was able to
create some noisy form of oscialting assemblies, but it was rough and at most
only 25% of the assemblies I created were present at a time.

This behavior occurs if the iterations are not allowed to fully converge or if
they do fully converge. It also occurs without step 3!

Updated notes:

Dropped step 3. The simulation appears to create non-noise oscillation
if β is small enough so that the assemblies don't converge too quickly
and they are stopped early. Hypothesis is that β is a matter of resolution
and what really matters is stopping early in the convergence process.

This could be a chaotic attractor!
"""

using AssemblyCalculus
using AssemblyCalculus: overlap, freeze_synapses!
using Plots
using StatsBase

K = 100
N = K^2
MEAN_DEGREE = K
P = MEAN_DEGREE / N

MAX_ITERS = 1000
TOL = floor(Int, 0.15 * K)

brain = Brain(num_areas=3, n=N, k=K, β=0.001, p=P)
assems, sp, ds = simulate!([rand_stim(brain[3])], max_iters=MAX_ITERS)
a0 = assems[1]

# Project alternating pattern across areas 1 and 2
assems, sp, ds = simulate!([rand_stim(brain[1])], [a0], tol = TOL, max_iters=MAX_ITERS)
a1 = assems[1]
assems, sp, ds = simulate!([rand_stim(brain[2])], [a0], tol = TOL, max_iters=MAX_ITERS)
a2 = assems[1]
assems, sp, ds = simulate!([rand_stim(brain[1])], [a0], tol = TOL, max_iters=MAX_ITERS)
a3 = assems[1]
assems, sp, ds = simulate!([rand_stim(brain[2])], [a0], tol = TOL, max_iters=MAX_ITERS)
a4 = assems[1]

println("overlap a1, a3: $(overlap(a1, a3))")
println("overlap a2, a4: $(overlap(a2, a4))")

freeze_synapses!(brain)
T = 200
assems, sp, ds = simulate!(
   [zero_stim(brain[1]), zero_stim(brain[2])],
   [a0],
   max_iters=T,
   tol=0,
)

assems = [a1, a2, a3, a4]
# Repeat spike to match assemby areas
sp = hcat(sp, sp)
assem_presence = [mapslices(x -> overlap(assems[i].neurons, x), sp[i], dims=1) for i in 1:length(assems)]
assem_presence = map(x -> reshape(x, :), assem_presence)

colors = [:teal, :pink, :purple, :orange]
names = ["a1", "a2", "a3", "a4"]
p1 = plot()
for (i, c, name) in zip(1:4, colors, names)
           plot!(p1, assem_presence[i], color=c, label=name, linewidth=3)
end
ylabel!("Area Overlap With Assembly")
xlabel!("Timestep")

p2 = plot()
plot!(p2, ds[1, :], color=colors[1], linewidth=3, label="Area 1")
plot!(p2, ds[2, :], color=colors[2], linewidth=3, label="Area 2")
xlabel!("Timestep")
ylabel!("Change in Area Neurons")

ap1, ap2, ap3, ap4 = assem_presence

p3 = plot()
plot!(crosscor(ap1, ap2), label="a2 Cross Correlation", linewidth=3)
plot!(crosscor(ap1, ap3), label="a3 Cross Correlation", linewidth=3)
plot!(crosscor(ap1, ap4), label="a4 Cross Correlation", linewidth=3)
ylabel!("Cross correlation with a1")
xlabel!("Lags")

p4 = plot()
plot!(crosscor(ap2, ap1), label="a1 Cross Correlation", linewidth=3)
plot!(crosscor(ap2, ap3), label="a3 Cross Correlation", linewidth=3)
plot!(crosscor(ap2, ap4), label="a4 Cross Correlation", linewidth=3)
ylabel!("Cross correlation with a2")
xlabel!("Lags")

plot(p1, p2, p3, p4, size=(1200, 1000), layout=@layout [a b; c d])

# assems_r = [a1r, a2r, a3r, a4r]
# # Repeat spike to match assemby areas
# sp = hcat(sp, sp)
# assem_presence = [mapslices(x -> overlap(assems_r[i].neurons, x), sp[i], dims=1) for i in 1:length(assems)]
#
# colors = [:teal, :pink, :purple, :orange]
# names = ["a1r", "a2r", "a3r", "a4r"]
# p2 = plot()
# for (i, c, name) in zip(1:4, colors, names)
#            plot!(p2, reshape(assem_presence[i], :), color=c, label=name, linewidth=3)
# end
# ylabel!("Area Overlap With Assembly")
# xlabel!("Timestep")
#
# plot(p1, p2, size=(1200, 400), layout=@layout [a b])


text:
"""
bound(estimator, δ, n::Int)
Give the theoritical bound B such that for the estimator
    #TODO! correction math formula
```math
\\mathbb{P}( |X- \\mathbb{E}(X)|/\\sigma \\leq B) \\leq \\delta
```

"""
bound(δ::Real, n, M::EmpiricalMean) = sqrt(1 / (n * δ))
bound(δ::Real, n, M::MedianOfMean) = sqrt(32log(1 / δ) / n)
bound(δ::Real, n, M::TrimmedMean) = 19sqrt(2log(8 / δ) / n) # Typo in the review by Lugosie. It seems that the ouragous 19 does th job. It is probably improvable.
bound(δ::Real, n, M::Catoni) = sqrt(2log(2 / δ) / (n - 2log(2 / δ)))
bound(δ::Real, n, M::Huber) = 8sqrt(2log(4 / δ) / n)
""""
sqrt(2log(1 / δ) / n) # Multiply by a factor (1 + o(1)) where
o(1) = (1+O(sqrt(log(1/δ)/n)))*(1+log(log(1/δ))/log(1/δ)) goes to zeros with (log(1/δ)/n, δ)→ (0, 0)
see https://www.youtube.com/watch?v=Kr0Kl_sXsJM Q&A
"""
bound(δ::Real, n, M::LeeValiant) = 8sqrt(2log(4 / δ) / n)


text:
"""
    INumericalValued(nexid)

An interface for numerical valued widgets.
"""
struct INumericalValued <: AbstractINextion
    _nid::NexID

    T
    rn  # RangeNumber type

    function INumericalValued(nid::NexID, T, rn)
        new(nid, T, rn)
    end
end

"""
    obj.value = new_val

Set value from `new_val` to Nextion object `obj`.
"""
function setproperty!(obj::INumericalValued, property::Symbol, new_val)

    if property == :value
        setnexproperty!(NexID(obj), :val, obj.T(new_val))
    else
        setfield!(obj, property, new_val)
    end

end

"""
    obj.value

Get `value` property from Nextion object `obj`.
"""
function getproperty(obj::INumericalValued, property::Symbol)

    if property == :value
        getnexproperty(NexID(obj), :val, obj.T)
    else
        getfield(obj, property)
    end

end


text:
using FFD
using Base.Test

function not_isapprox(args...; kwargs...)
  return !isapprox(args...; kwargs...)
end

# write your own tests here
@test 1 == 1

include("test_new.jl")


text:
export ModelDesigner

using ...Layer

struct LayerDetails
    design::LayerDesign
    input_ptr::PrealloctionPointer
    output_ptr::PrealloctionPointer
    const_idx

    function LayerDetails(layer::LayerDesign, pb::ProgramBuilder, initial_constants)
        S = pb.device.storagetype
        initializer = layer.const_initializer !== nothing ? layer.const_initializer(S) : nothing
        const_idx = length(initial_constants) + 1

        for c in 1:prod(l.const_shape)
            push!(initial_constants, initializer(c))
        end

        new(layer, prealloc_stack!(pb, S, layer.input_shape...), prealloc_stack!(pb, S, layer.output_shape...), const_idx)
    end
end

mutable struct ModelDesigner{StorageType} <: AbstractObject
    layers::Vector{LayerDetails}
    initial_constants::Vector{StorageType}
    error_function::Error.Func
    constant_range
    postprocessor::Processor.Func
    last_output_shape
    memsett::MemorySettings
    pb::ProgramBuilder
    was_built::Bool

   function ModelDesigner(; input_shape, settings=MemorySettings(), error::Error.Func = Error.None(),
                        postprocessor::Processor.Wrapper=Processor.None(),
                        constant_range = -1.0:.001:1.0)

        return new{settings.storagetype}(
            Vector{Layers.LayerDetails}(undef, 0)
            error,
            constant_range,
            postprocessor(settings),
            input_shape,
            settings,
            ProgramBuilder(d.settings.device),
            false)
   end

    function (d::ModelDesigner{S})() where {S}
        (length(d.layers) == 0) && return nothing
        (d.was_built) && error("Cannot Rebuild from same designer. Please make another")
        d.was_built = true
        constants = prealloc_global!(d.pb, :Constants, S, length(d.initial_constants))
        program = pb()
        constants = constants()
        copy!(constants, d.initial_constants)
        fun = (inputs) -> copy!(layer_inputs, inputs)
        last_output_array = nothing

        for layer in d.layers
            layer_inputs = layer.input_ptr()
            layer_outputs = layer.input_ptr()
            layer_const_shape = layer.design.const_shape
            layer_constants = reshape(view(constants, layer.const_idx:(layer.const_idx + prod(layer_const_shape))), layer_const_shape...)
            layer_fun = layer(layer_inputs, layer_outputs, layer_constants)
            last_output_array = layer_outputs

            fun = function (inputs)
                    fun(inputs)
                    layer_fun()
                  end
        end

        proc = d.postprocessor
        fun = function (inputs)
             fun(inputs)
             return proc(layer_outputs)
        end

        return Model(d, fun, program, constants, d.constant_range, d.settings.arrayType{d.settings.inputType}, d.settings.outputType)
    end

    "Do dimension checking and constant initializing"
    function prepare_layer(designer::ModelDesigner{S}, layer::Layers.LayerDesign) where S
        is_reshape_layer = layer.sig.name == :ReshapeLayer
        ((designer.last_output_shape !=  layer.input_shape) || is_reshape_layer) || error("Incorrect Dimensions at Layer:$(length(designer.layers)).
                 Last Output Shape:$(designer.last_output_shape). Recieved:$(layer.input_shape)")
        is_reshape_layer || push!(designer.layers, LayerDetails(layer, designer.pb, designer.initial_constants))
        designer.last_output_shape = layer.output_shape
    end

    function Base.push!(designer::ModelDesigner, layer::LayerGenerator)
        prepare_layer(designer, layer(designer.pb, designer.memsett, designer.last_output_shape))
    end

    function Base.append!(designer::ModelDesigner, layers::LayerGenerator...)
        for layer in layers
            push!(designer, layer)
        end
    end

    function Base.string(x::ModelDesigner)
        layerz = join([string(layer) for layer in x.layers], ", ")
        const_count = length(x.init_constants)
        err_func = x.error_function
        return "ModelDesigner:\n\tLayers: $layerz\n\tConstant Count: $const_count\n\tError Function: $err_func\n"
    end
end


text:
#  github: part of RxNav.jl

"""
    getAllConcepts

/allconcepts	All RxTerms concepts
"""
function getAllConcepts()
    clist = Vector{String}[]
    try
        doc = getdoc("baseurl", "RxTerms/allconcepts")
        for x in findall("//minConceptGroup/minConcept", doc)
            push!(clist, [nodecontent(findfirst(s, x)) for s in ["fullName", "termType", "rxcui"]])
        end
    catch y
        @warn y
    end
    return ndclist
end

"""
    getAllRxTermInfo

/rxcui/rxcui/allinfo	RxTerms information for a specified RxNorm concept
Returns the (quite variable) properties data in XML form.
"""
function getAllRxTermInfo(rxcui::String)
    argstring = "RxTerms/rxcui/" * rxcui * "/allinfo"
    try
        doc = getdoc("baseurl", argstring)
        rxn = findfirst("//rxtermsProperties", doc)
        return string(rxn)
    catch y
        @warn y
        return ""
    end
end

"""
    getRxTermDisplayName

/rxcui/rxcui/name	RxTerms display name for a specified RxNorm concept
"""
function getRxTermDisplayName(rxcui::String)
    argstring = "RxTerms/rxcui/" * rxcui * "/name"
    try
        doc = getdoc("baseurl", argstring)
        x = findfirst("//displayGroup", doc)
        return nodecontent(findfirst("displayName", x))
    catch y
        @warn y
        return ""
    end
end

"""
getRxTermsVersion

/version	RxTerms version
"""
function getRxTermsVersion()
    try
        doc = getdoc("baseurl", "rxTerms/version")
        return nodecontent(findfirst("//rxTermsVersion", doc))
    catch y
        @warn y
        return ""
    end
end


text:
## Implements bitonic merge networks that can merge two sorted SIMD vectors.
## The generated function supports any possible type and size.

using SIMD

"""
    bitonic_merge(input_a::Vec{N,T}, input_b::Vec{N,T}) where {N,T}

Merges two `SIMD.Vec` objects of the same type and size using a bitonic sort network. The inputs are assumed to be sorted. Returns a pair of vectors with the first and second halves of the merged sequence.
"""
@generated function bitonic_merge(input_a::Vec{N,T}, input_b::Vec{N,T}) where {N,T}

    pat = Val{ntuple(x->N-x, N)}

    ex = [
        Expr(:meta, :inline),
        :(la_0 = input_a),
        :(lb_0 = shufflevector(input_b, $pat)),
        :(L_0 = min(la_0, lb_0)),
        :(H_0 = max(la_0, lb_0))
    ]

    p = mylog(N)
    for n in 1:p
        la = Symbol("la_", n)
        lb = Symbol("lb_", n)
        Lp = Symbol("L_", n-1)
        Hp = Symbol("H_", n-1)
        L = Symbol("L_", n)
        H = Symbol("H_", n)

        ih = inverse_shuffle(bitonic_step(2^(n), 2^(p-n+1)))
        sh = bitonic_step(2^(n+1), 2^(p-n))
        pat_a = Val{tuple((ih[sh[1:N]].-1)...)}
        pat_b = Val{tuple((ih[sh[(N+1):end]].-1)...)}

        append!(ex, [
            :($la = shufflevector($Lp, $Hp, $pat_a)),
            :($lb = shufflevector($Lp, $Hp, $pat_b)),
            :($L = min($la, $lb)),
            :($H = max($la, $lb))
        ])
    end

    la = Symbol("la_", p+1)
    lb = Symbol("lb_", p+1)
    Lp = Symbol("L_", p)
    Hp = Symbol("H_", p)

    ih = inverse_shuffle(bitonic_step(2N, 1))
    pat_a = Val{tuple((ih[1:N].-1)...)}
    pat_b = Val{tuple((ih[(N+1):end].-1)...)}

    append!(ex, [
        :($la = shufflevector($Lp, $Hp, $pat_a)),
        :($lb = shufflevector($Lp, $Hp, $pat_b)),
        :(($la, $lb))
    ])

    quote $(ex...) end
end

function inverse_shuffle(tt)
    out = copy(tt)
    for n in 1:length(tt)
        out[tt[n]] = n
    end
    out
end

function bitonic_step(G, Glen)
    N = G * Glen
    gg = reshape(1:N, Glen, G)
    t_lh = [gg[:,1:2:end][:]; gg[:,2:2:end][:]]
end

"""
    merge_vecs(input...)

This brave function merges 4, 8 vectors or even more. Relies on the great ability from SIMD.jl and LLVM to handle
large vectors. Performance not yet tested."""
@generated function merge_vecs(input::Vararg{Vec{N,T}, L})::Vec{N*L,T} where {L,N,T}

    ex = [Expr(:meta, :inline)]

    for l in 1:L
        m = Symbol("m_0_", l)
        push!(ex, :($m = input[$l]))
    end

    p = mylog(L)
    for l in 1:p
        for b in 1:2^(p-l)
            m = Symbol("m_", l, "_", b)
            a1 = Symbol("m_", l-1, "_", b*2-1)
            a2 = Symbol("m_", l-1, "_", b*2)
            push!(ex, :($m = bitonic_merge_concat($a1, $a2)))
        end
    end

    m = Symbol("m_", p, "_", 1)
    push!(ex, :(return $m))

    quote $(ex...) end
end

@inline function bitonic_merge_concat(input_a::Vec{N,T}, input_b::Vec{N,T})::Vec{2*N,T} where {N,T}
    m_a, m_b = bitonic_merge(input_a, input_b)
    concat(m_a, m_b)
end

"""
    bitonic_merge_interleaved(output_a, output_b, input_a, input_b, Val(N))

Merges K pairs of vectors using a bitonic sort network, with interleaved execution. The inputs are assumed to be sorted. Inputs and outputs are acessed directly in memory.
"""
@generated function bitonic_merge_interleaved(
    output::AbstractArray{T,3},
    input::Array{Ptr{T},1},
    ::Val{V},
    J,
    ::Val{K}
) where {T,V,K}

    allex = []
    for indy in 1:K
        iSymbol(x...) = Symbol("indy_", indy, "_", x...)
        pat = Val{ntuple(x->V-x, V)}

        n = 0
        la = iSymbol("la_", n)
        lb = iSymbol("lb_", n)
        L = iSymbol("L_", n)
        H = iSymbol("H_", n)
        ex = [
            :($la = vloada(Vec{V,T}, pointer(@view output[:,J,$indy]))),
            :($lb = shufflevector(vloada(Vec{V,T}, input[$indy]), $pat)),
            :($L = min($la, $lb)),
            :($H = max($la, $lb))
        ]

        p = mylog(V)
        for n in 1:p
            la = iSymbol("la_", n)
            lb = iSymbol("lb_", n)
            Lp = iSymbol("L_", n-1)
            Hp = iSymbol("H_", n-1)
            L = iSymbol("L_", n)
            H = iSymbol("H_", n)

            ih = inverse_shuffle(bitonic_step(2^(n), 2^(p-n+1)))
            sh = bitonic_step(2^(n+1), 2^(p-n))
            pat_a = Val{tuple((ih[sh[1:V]].-1)...)}
            pat_b = Val{tuple((ih[sh[(V+1):end]].-1)...)}

            append!(ex, [
                :($la = shufflevector($Lp, $Hp, $pat_a)),
                :($lb = shufflevector($Lp, $Hp, $pat_b)),
                :($L = min($la, $lb)),
                :($H = max($la, $lb))
            ])
        end

        la = iSymbol("la_", p+1)
        lb = iSymbol("lb_", p+1)
        Lp = iSymbol("L_", p)
        Hp = iSymbol("H_", p)

        ih = inverse_shuffle(bitonic_step(2V, 1))
        pat_a = Val{tuple((ih[1:V].-1)...)}
        pat_b = Val{tuple((ih[(V+1):end].-1)...)}

        append!(ex, [
            :(vstorea(shufflevector($Lp, $Hp, $pat_a), pointer(@view output[:,J,$indy]))),
            :(vstorea(shufflevector($Lp, $Hp, $pat_b), pointer(@view output[:,J+1,$indy])))
        ])
        push!(allex, ex)
    end

    inteleaved = [x for xx in zip(allex...) for x in xx]

    final_ex = [
        Expr(:meta, :inline),
        inteleaved...
    ]

    quote $(final_ex...) end
end


text:
using Test
using StaticArrays, IntervalSets
import ClimaCore.DataLayouts: IJFH
import ClimaCore:
    Fields, Domains, Meshes, Topologies, Spaces, Operators, Geometry
using StaticArrays, IntervalSets, LinearAlgebra

include("sphere_sphericalharmonics.jl")

# diffusion on vector: input covariant vecotor
function ∇²(u)
    scurl = Operators.Curl()
    sdiv = Operators.Divergence()
    wcurl = Operators.WeakCurl()
    wgrad = Operators.WeakGradient()

    χ = Spaces.weighted_dss!(
        @. wgrad(sdiv(u)) - Geometry.Covariant12Vector(
            wcurl(Geometry.Covariant3Vector(scurl(u))),
        )
    )
    return χ
end

@testset "diffusion on vectors" begin
    FT = Float64

    l = Int(7)
    m = Int(4)

    Ne = 7
    Nq = 6

    radius = FT(1)
    domain = Domains.SphereDomain(radius)
    mesh = Meshes.Mesh2D(domain, Meshes.EquiangularSphereWarp(), Ne)
    grid_topology = Topologies.Topology2D(mesh)

    quad = Spaces.Quadratures.GLL{Nq}()
    space = Spaces.SpectralElementSpace2D(grid_topology, quad)
    coords = Fields.coordinate_field(space)

    # compute vector spherical harmonics (VSH) as the gradient of scalar spherical harmonics
    # eigenfunction properties of VSH: ∇² VSH(l,m) = -l(l+1) VSH(l,m)
    VSH_local = map(coords) do coord
        Geometry.UVVector(VSH(l, m, coord.lat, coord.long)...)
    end
    VSH_cov = Geometry.transform.(Ref(Geometry.Covariant12Axis()), VSH_local)

    # diffusion
    ∇²VSH_cov = ∇²(VSH_cov)
    ∇²VSH_local = Geometry.transform.(Ref(Geometry.UVAxis()), ∇²VSH_cov)

    # exact solution
    ∇²VSH_exact = @. -l * (l + 1) * VSH_local

    @test ∇²VSH_local ≈ ∇²VSH_exact rtol = 2e-2
end

convergence_rate(err, Δh) =
    [log(err[i] / err[i - 1]) / log(Δh[i] / Δh[i - 1]) for i in 2:length(Δh)]

@testset "convergence tests for vector diffusions on the sphere" begin
    FT = Float64

    Nes = [4, 8, 16]
    Nqs = [3, 4, 5, 6]

    l = Int(7)
    m = Int(4)

    radius = FT(1)
    domain = Domains.SphereDomain(radius)
    for (Iq, Nq) in enumerate(Nqs)
        err_diffusion = Array{FT}(undef, length(Nes))
        Δh = Array{FT}(undef, length(Nes))

        for (Ie, Ne) in enumerate(Nes)
            mesh = Meshes.Mesh2D(domain, Meshes.EquiangularSphereWarp(), Ne)
            grid_topology = Topologies.Topology2D(mesh)

            quad = Spaces.Quadratures.GLL{Nq}()
            space = Spaces.SpectralElementSpace2D(grid_topology, quad)
            coords = Fields.coordinate_field(space)

            VSH_local = map(coords) do coord
                Geometry.UVVector(VSH(l, m, coord.lat, coord.long)...)
            end
            VSH_cov =
                Geometry.transform.(Ref(Geometry.Covariant12Axis()), VSH_local)
            ∇²VSH_cov = ∇²(VSH_cov)
            ∇²VSH_local = Geometry.transform.(Ref(Geometry.UVAxis()), ∇²VSH_cov)
            ∇²VSH_exact = @. -l * (l + 1) * VSH_local

            err_diffusion[Ie] = norm(∇²VSH_local .- ∇²VSH_exact)
            Δh[Ie] = 1 / Ne
        end

        convergence_rate_∇² = convergence_rate(err_diffusion, Δh)
        for Ie in range(1, length = length(Nes) - 1)
            @test convergence_rate_∇²[Ie] > (Nq - 2)
        end
    end
end


text:
# Include packages
using Interpolations
using Printf: @printf, @sprintf
using Polynomials: polyfit, polyval
using DSP: welch_pgram
using DelimitedFiles: readdlm
using Dates

####################
# function
####################
function loadcsvsample()
    # define the filepath & filename
    fdir = "./data"
    fname = "dat_climate.csv"
    dataorg = readdlm(joinpath(fdir,fname), ',', skipstart=1)
    return dataorg
end
####################

####################
## main
####################

# read csv data
dataorg = loadcsvsample()
torg = DateTime.(dataorg[:,1], "yyyy/mm/dd HH:MM")
tsecorg = convert.(Float64, Dates.value.((torg-torg[1])/1000)) # ms to s
V = convert.(Float64, dataorg[:,2])
dataorg = nothing

# Interpolate
t = collect(torg[1]:Day(1):torg[end])
tsec = convert.(Float64, Dates.value.((t-t[1])/1000)) # ms to s
nt = size(t,1)
itp = interpolate((tsecorg, ), V, Gridded(Linear()))
#itp = interpolate((tsecorg, ), V, Gridded(Constant()))
Vint = itp(tsec)

# regression
lin_p = polyfit(tsec, Vint, 1)

# power spectral
# (注)360*2は結果を合わせにいった値のため根拠なし
pdg = welch_pgram(Vint, 360*2, onesided=true; fs=1.0)
days = 1.0./convert.(Float64, pdg.freq)
PSD = pdg.power
maxval, Tc = findmax(PSD[2:end])

####################

####################
## plot
####################

# directory where figures are printed
figdir="./fig"
if !isdir(figdir); mkdir(figdir); end

using Plots
pyplot()
#gr()
#plotlyjs()

# ticklabel format for Figure 1
xt = collect(DateTime(2012):Year(1):DateTime(2017))
xtl = Dates.format.(xt, DateFormat("Y"))
xt = Dates.value.(xt)

########
## Figure 1
#  trend
# constant interval
plot(t, Vint, line=(:solid, 1), lab="Uniformed", size=(800,600),
     tickfont=12, legend=:topleft, legendfont=14,
     xticks=(xt, xtl), xlims=(xt[1]-(0.5*365*86400*1e3), xt[end]),
     xlabel="Year", guidefont=12,
     ylims=(62,72),
     )
# plot!(torg, V, line=(:dash, 1), lab="Raw")

# regression
plot!(t, polyval(lin_p, tsec), lab="Regression 1")

# save
savefig(joinpath(figdir,"uniformed_data.png"))
########

########
## Figure 2
#  periods
plot(days[2:end], PSD[2:end], xscale=:log10, tickfont=12, size=(800, 600),
     lab="Power", legend=:topleft, legendfont=14,
     xlabel="Days", guidefont=14,
     xlims=(1e+00, 1e3),
     )
# natral frequency
scatter!([days[Tc+1]], [PSD[Tc+1]], ms=10.,
         lab="Period: "*@sprintf("%.0f",days[Tc+1])*"days",
         )
         #ylims=(0., 550.),
# save
savefig(joinpath(figdir,"trend.png"))
########

####################


text:
module TestSolver

using AMGX, Defer, Test, JSON, CUDA, SparseArrays
using AMGX: Config, Resources, AMGXVector, AMGXMatrix, Solver, dDDI, dFFI

@scope @testset "Solver" begin
    mode = dDDI

    @scope @testset "solver defaults" begin
        c = @! Config("")
        r = @! Resources(c)
        M = @! AMGXMatrix(r, mode)
        x = @! AMGXVector(r, mode)
        b = @! AMGXVector(r, mode)
        s = @! Solver(r, mode, c)

        AMGX.upload!(M, Cint[0, 1, 2, 3], Cint[0, 1, 2], [1.0, 1.0, 1.0])
        AMGX.upload!(x, [0.0, 0.0, 0.0])
        AMGX.upload!(b, [1.0, 2.0, 4.0])

        AMGX.setup!(s, M)
        AMGX.solve!(x, s, b)
        x_h = AMGX.download(b)
        @test x_h ≈ [1.0, 2.0, 4.0]
        status = AMGX.get_status(s)
        @test status == AMGX.SUCCESS

        AMGX.upload!(x, [1.0, 2.0, 3.0])
        AMGX.solve!(x, s, b; zero_inital_guess=true)
        x_h = AMGX.download(x)
        @test x_h ≈ [1.0, 2.0, 4.0]
        status = AMGX.get_status(s)
        @test status == AMGX.SUCCESS
    end

    @scope @testset "get_status / get_iterations_number" begin
        c = @! Config(Dict("monitor_residual" => 1, "max_iters" => 0, "store_res_history" => 1))
        r = @! Resources(c)
        M = @! AMGXMatrix(r, mode)
        x = @! AMGXVector(r, mode)
        b = @! AMGXVector(r, mode)
        s = @! Solver(r, mode, c)

        AMGX.upload!(M, Cint[0, 1, 2, 3], Cint[0, 1, 2], [1.0, 1.0, 1.0])
        AMGX.upload!(x, [0.0, 0.0, 0.0])
        AMGX.upload!(b, [1.0, 2.0, 4.0])
        AMGX.setup!(s, M)
        AMGX.solve!(x, s, b; zero_inital_guess=true)
        @test AMGX.get_status(s) == AMGX.DIVERGED
        @test AMGX.get_iterations_number(s) == 0
        @test AMGX.get_iteration_residual(s) > 1.0

        close(M)
        @test_throws ErrorException AMGX.solve!(x, s, b; zero_inital_guess=true)
    end

    @scope @testset "get_iterations_number" begin
        c = @! Config(Dict("monitor_residual" => 1, "tolerance" => 1e-14, "store_res_history" => 1))
        r = @! Resources(c)
        M = @! AMGXMatrix(r, mode)
        x = @! AMGXVector(r, mode)
        b = @! AMGXVector(r, mode)
        s = @! Solver(r, mode, c)

        AMGX.upload!(M, Cint[0, 1, 2, 3], Cint[0, 1, 2], [1.0, 1.0, 1.0])
        AMGX.upload!(x, [0.0, 0.0, 0.0])
        AMGX.upload!(b, [1.0, 2.0, 4.0])
        AMGX.setup!(s, M)
        AMGX.solve!(x, s, b)
        @test AMGX.get_iteration_residual(s) < 1e-14
        niter = AMGX.get_iterations_number(s)
        @test AMGX.get_iteration_residual(s, niter) == AMGX.get_iteration_residual(s)
    end
end

end # module


text:
using Cameras
using Test

@testset "Iteration" begin
    period = 0.04 # seconds
    image_size = (2, 2)

    function produce_images!(image_source::Channel)
        while true
            put!(image_source, zeros(image_size))
        end
    end
    image_source = Channel(produce_images!)

    camera = SimulatedCamera(period, image_source)

    start!(camera)

    iterations = 5
    let i = 0
        t_0 = time()
        for img in camera
            i += 1
            if i >= iterations
                break
            end
        end
        t_1 = time()
        @test i == iterations
        @test t_1 - t_0 >= iterations * period
    end

    t_0 = time()
    for img in Iterators.take(camera, iterations)
    end
    t_1 = time()
    @test t_1 - t_0 >= iterations * period

    t_0 = time()
    images = [img for img in Iterators.take(camera, iterations)]
    t_1 = time()
    @test t_1 - t_0 >= iterations * period
end


text:
module PostEstim

export MonteCarloIncome
using TypeDef, Helper, Estim, Gmm, GMMInterface, Optim, UserFileIO
using Printf, DelimitedFiles, Random, Distributions, LinearAlgebra

function exportδ(
  θ₂,
  dat  ::datType,
  dim  ::dimType,
  fxp  ::fxpType,
  mat  ::matType,
  indx ::indxType,
  pred ::predType,
  diags::diagsType,
  fid  ::IOStream)

  setupMu!(fxp.eμ,pred.eμ₀,pred.adj, θ₂,mat.Xᵢ,mat.Dᵢ,mat.νᵢ,dim,indx)
  fixPoint!(fxp.eδ,fxp.∂δ∂θ, fxp.eμ,dat,dim,mat,indx,pred,diags,θ₂)
  unpackStackedArray!(fxp.eδ1d, fxp.eδ, dim.market, indx.mkt)
  for i=1:dim.obs
    @printf(fid,", %28.15f",fxp.eδ1d[i])
  end
  @printf(fid,"\n")

end

function MonteCarloIncome(
  model  ::Int64,
  mcid   ::Int64,
  nRep   ::Int64,
  iv     ::String,
  ver    ::String,
  impF   ::String,
  outp   ::Bool;
  lowB   ::Bool=true,
  simEmp ::Bool=false,
  nPerson::Int64=200,
  stage  ::Int64=0,
  mcid1st::Int64=0,
  update ::Bool=true,
  fileName="paramDraws.csv")

  # Import
  dat,dim,fxp,lab,mat,set,indx,pred,diags =
    initialize(model,iv,nPerson,ver,impF,simEmp)

  pathName = string("out/",impF,"/m",model,"s",stage,iv,nPerson)
  imp = readdlm(string(pathName,"/gmmStage2.csv"),',')
  θ̂   = convert(Vector{Float64},imp[mcid,2:(end-2)])
  θ̂₁  = θ̂[dim.θ₂+1:end]
  θ̂₂  = θ̂[1:dim.θ₂]

  seθ̂,θ̂₁,covθ̂ = stdError(θ̂₂,model,stage,iv,nPerson,impF,update,
                         dat,dim,fxp,mat,set,indx,pred,diags,mcid1st)

  cholθ̂ = cholesky(Hermitian(covθ̂))

  # Import bounded income draws (for MC sim. Turn off for income prediction)
  dat,dim,fxp,lab,mat,set,indx,pred,diags =
    initialize(model,iv,nPerson,ver,impF,simEmp;lowBound=lowB)

  # Monte Carlo
  rng = MersenneTwister(1234)
  inc = zeros(dim.brand,dim.prepost,nRep)
  fid = open(string(pathName,"/",fileName),"w")
  for r=1:nRep
    tic = time()
    θ₁,θ₂ = sampleθ(θ̂,cholθ̂,dat,dim,fxp,mat,indx,pred,diags,rng,fid,outp)
    y = simulateIncome(θ₁,θ₂,dat,dim,fxp,mat,indx,pred,rng)
    @. inc[:,:,r] = y
    @printf("Bootstrap %d complete in %.1f sec.\n",r,time()-tic)
  end
  close(fid)
  yAvg = mean(inc,dims=3)
  y025 = zeros(dim.brand,dim.prepost)
  y975 = zeros(dim.brand,dim.prepost)

  for k=1:dim.prepost
    for j=1:dim.brand
      y025[j,k] = quantile(inc[j,k,:],0.025)
      y975[j,k] = quantile(inc[j,k,:],0.975)
    end
  end
  if lowB==false
    @printf("Average Income by Subscribers of Product Lines\n")
    display("text/plain",yAvg)
  end
  return yAvg, y025, y975
end

function sampleθ(
  θ̂::Vector{Float64},
  cholθ̂,
  dat  ::datType,
  dim  ::dimType,
  fxp  ::fxpType,
  mat  ::matType,
  indx ::indxType,
  pred ::predType,
  diags::diagsType,
  rng  ::MersenneTwister,
  fid  ::IOStream,
  outp ::Bool=false)

  θ  = θ̂ + cholθ̂.L * randn(rng,dim.θ)
  θ₁ = θ[dim.θ₂+1:end]
  θ₂ = θ[1:dim.θ₂]
  while (minimum(θ₂[1:2])<=0)
    θ  .= θ̂ + cholθ̂.L * randn(rng,dim.θ)
    θ₁ .= θ[dim.θ₂+1:end]
    θ₂ .= θ[1:dim.θ₂]
    σₐ  = θ₂[1]
    if (minimum(θ₂[1:2])<=0)
      @printf("\nMonte Carlo draw rejected: %14.4f %14.4f \n",θ₂[1],θ₂[2])
    end
  end
  # Export parameters to csv file
  if outp==true # export θ & δ to csv file
    @printf(fid,"%22.15f",θ₂[1])
    for k=2:dim.θ₂
      @printf(fid,", %22.15f",θ₂[k])
    end
    for k=1:dim.θ₁
      @printf(fid,", %22.15f",θ₁[k])
    end
    exportδ(θ₂,dat,dim,fxp,mat,indx,pred,diags,fid)
  end

  return θ₁,θ₂
end

function simulateIncome(
  θ₁::Vector{Float64},
  θ₂::Vector{Float64},
  dat::datType,
  dim::dimType,
  fxp::fxpType,
  mat::matType,
  indx::indxType,
  pred::predType,
  rng::MersenneTwister)

  yMean  = zeros(dim.brand,dim.prepost)
  yCount = zeros(dim.brand,dim.prepost)

  μ = 0; σ = 1; ξ = 0;
  randϵ = GeneralizedExtremeValue(μ, σ, ξ)

  @inbounds for t=1:dim.market
    eδ = fxp.eδ[t]
    eμ = fxp.eμ[t]
    Dᵢ = view(mat.Dᵢ[t],1,1,:)
    num     = zeros(dim.prodt[t])
    mask    = view(indx.mkt,:,t)
    brand   = dat.brand[mask]
    prepost = dat.prepost[mask]
    util = zeros(dim.prodt[t],dim.person)
    for i=1:dim.person
      for j=1:dim.prodt[t]
        util[j,i] = log(eδ[j]*eμ[j,i])
      end
    end

    for i=1:dim.person
      μ₀ = log(pred.eμ₀[t][i])
      for r=1:100
        for j=1:dim.prodt[t]
          num[j] = util[j,i] + rand(randϵ)
        end
        jmax = argmax(num)
        if num[jmax] > μ₀
          yMean[brand[jmax],prepost[jmax]] += 100/Dᵢ[i]
          yCount[brand[jmax],prepost[jmax]] += 1
        end
      end
    end # end of i loop
  end # end of t loop
  yCount[yCount.==0] .= 1
  @. yMean = yMean / yCount
  return yMean
end

end

max_stars_repo_name: yutec/FightBrand

max_issues_repo_name: yutec/FightBrand

max_forks_repo_name: yutec/FightBrand