text:
using MDPTools, NestedFixedPoint, Plots
srand(122)

n = [0, 2, 5]
next_i(i, k) = max(0, min(I, i+n[k]-1))

p = [0.6, 0.8, 1.] # FIX this
F_P = [0.15 0.0 0.85; 0.0 0.25 0.75; 0.10 0.10 0.80]
SP = CommonState(p, F_P)

I = 15
K = length(n)
F_i = [spzeros(I+1, I+1) for k = 1:K]
for i = 0:I, j = 1:K
        F_i[j][i+1, next_i(i, j)+1]=1
end

Si = State(0:I, F_i);

S = States(SP, Si)

delta, alpha = 3., 3.

Z1 = [[zeros(1); ones(I)] zeros(I+1) -next_i(0:I, 1);
      [zeros(1); ones(I)] zeros(I+1) -next_i(0:I, 1);
      [zeros(1); ones(I)] zeros(I+1) -next_i(0:I, 1)]
Z2 = [ones(3*(I+1)) -[p[1]*n[2]*ones(I+1);p[2]*n[2]*ones(I+1);p[3]*n[2]*ones(I+1)] -kron(ones(3), next_i(0:I, 2))]
Z3 = [ones(3*(I+1)) -[p[1]*n[3]*ones(I+1);p[2]*n[3]*ones(I+1);p[3]*n[3]*ones(I+1)] -kron(ones(3), next_i(0:I, 3))]
U = LinearUtility([Z1, Z2, Z3], 0.95, [delta; alpha; 0.05])

M, N, T = 1, 400, 3*12
D = simulate(U, S, M, N, T);

hitsch_nfxp = fit_nfxp(U, S, D)

theme(:dark)
plot(plot([U.P[i][1:16] for i = 1:3], labels = ["buy 0" "buy 2" "buy 5"], title = "Low price", xlabel = "inventory", ylabel = "probability", ylabel = "probability"),
            plot([U.P[i][17:32] for i = 1:3], labels = ["buy 0" "buy 2" "buy 5"], title = "Medium price", xlabel = "inventory", ylabel = "probability"),
            plot([U.P[i][33:end] for i = 1:3], labels = ["buy 0" "buy 2" "buy 5"], title = "High price", xlabel = "inventory", ylabel = "probability"), layout = (1,3), size = (900,300))


text:
export meas, viewmeas, state, viewstate

"""
    meas(X, Ny, Nx)

A function to extract the first `Ny` lines of the ensemble matrix `X`, typically storing the observations.
"""
meas(X, Ny, Nx) = X[1:Ny,:]

"""
    viewmeas(X, Ny, Nx)

Create a view of the first `Ny` lines of the ensemble matrix `X`, typically storing the observations.
"""
viewmeas(X, Ny, Nx) = view(X, 1:Ny,:)

"""
    state(X, Ny, Nx)

A function to extract the lines `Ny+1` to `Ny+Nx` of the ensemble matrix `X`, typically storing the state.
"""
state(X, Ny, Nx) = X[Ny+1:Nx+Ny,:]

"""
    viewstate(X, Ny, Nx)

Create a view of the lines `Ny+1` to `Ny+Nx` of the ensemble matrix `X`, typically storing the state.
"""
viewstate(X, Ny, Nx) = view(X, Ny+1:Nx+Ny,:)

ext: jl


text:
#
# This file is part of module
#   ModiaMedia (ModiaMedia/src/ModiaMedia.jl)
#

### Data structures that are additionally available for media <: PartialMixtureMedium -----------------------

"""
    state = ThermodynamicState_pTX(p,T,X)

Generate a `ThermodynamicState_pT <: MixtureThermodynamicState` object containing
pressure `p` [Pa], temperature `T` [K], and a vector of mass fractions `X` as states.
"""
mutable struct ThermodynamicState_pTX <: MixtureThermodynamicState
    p::Float64
    T::Float64
    X::Vector{Float64}
end

"R = gas_constant(medium,state) - return gas constant for `medium` from `state` in [J/mol.K]"
# should this be MixtureMedium?
# should this take MixtureThermodynamicState?
gas_constant(m::AbstractMedium, state::ThermodynamicState) = undefinedFunction("gas constant", m)

# moleToMassFractions: use inner constructor with fluid information.
# massToMoleFractions: use inner constructor with fluid information from abstract medium, since we have that type.


text:
function _prob_contact(contact_mean::Int, num_persons::Int, num_places::Int)
    return contact_mean / (num_persons / num_places) * 100
end

function _get_random_size(range::String, max::Int)
    intervals = split(range, '-')
    if contains(intervals[1], "+")
        intervals[1] = intervals[1][1:(length(intervals[1])-1)]
        intervals = parse.(Int, intervals)
        push!(intervals, max)
    else
        intervals = parse.(Int, intervals)
    end
    return rand(intervals[1]:intervals[2])
end

function _get_occupation_by_age(age::Int)
    p_quit_school = 0.17
    p_unoccupied = 0.17
    p_unoccupied_15_29 = 0.32

    if age <= 17
        return SCHOOL
    elseif 18 <= age <= 24
        if rand() > p_quit_school
            return SCHOOL
        else
            if rand() > p_unoccupied_15_29
                return WORKING
            else
                return UNEMPLOYED
            end
        end
    elseif 25 <= age <= 29
        if rand() > p_unoccupied_15_29
            return WORKING
        else
            return UNEMPLOYED
        end
    elseif age < 65
        if rand() > p_unoccupied
            return WORKING
        else
            return UNEMPLOYED
        end
    else
        return UNEMPLOYED
    end
end

function _get_school_level(age, levels)
    for level in levels
        intervals = split(level[1], '-')
        if contains(intervals[1], "+")
            intervals[1] = intervals[1][1:(length(intervals[1])-1)]
            intervals = parse.(Int, intervals)
            if age >= intervals[1]
                return level[2]
            end
        else
            intervals = parse.(Int, intervals)
            if age >= intervals[1] && age <= intervals[2]
                return level[2]
            end
        end
    end
end

function _get_school_work_place(
    age::Int,
    occupation::Int,
    companies::Set{Company},
    schools::Dict{String,Array{Int,1}},
    levels,
)
    if occupation == UNEMPLOYED
        return -1
    elseif occupation == SCHOOL
        school_level = _get_school_level(age, levels)
        return rand(schools[school_level])
    else
        company = rand(companies)
        company.current_size += 1
        return company.id
    end
end

function _add_person_to_right_set(person::Person, students, employeds, unemployeds)
    occupation = person.occupation
    if occupation == SCHOOL
        push!(students, person)
    elseif occupation == WORKING
        push!(employeds, person)
    elseif occupation == UNEMPLOYED
        push!(unemployeds, person)
    end
end

"""
Write an household vector as a CSV dataset
"""
function _write_households(households, output_path::String)
    output_path = output_path * "dataset.csv"
    table = DataFrame(
        household_id = Int[],
        person_id = Int[],
        age = Int[],
        occupation = Int[],
        work_school_place = Int[],
    )
    mutex = ReentrantLock()
    Threads.@threads for household in households
        people = household.workers ∪ household.students ∪ household.unemployeds
        for person in people
            lock(mutex) do
                push!(
                    table,
                    [
                        household.id,
                        person.id,
                        person.age,
                        person.occupation,
                        person.work_school_place,
                    ],
                )
            end
        end
    end
    table |> CSV.write(output_path)
end

function _write_dict(input,output_path::String,filename::String)
    output_path = output_path * filename *".json"
    open(output_path,"w+") do io
        JSON3.write(io,input)
    end
end

function _read_household_dataset(dataset_path::String)
    households = Array{Household,1}()
    df = CSV.File(dataset_path) |> DataFrame
    gdf = groupby(df, :household_id)
    mutex = ReentrantLock()
    Threads.@threads for group_idx = 1:length(gdf)
        group = gdf[group_idx]
        household = Household(group.household_id[1], 0)
        for record in eachrow(group)
            if record.occupation == WORKING
                push!(
                    household.workers,
                    Person(
                        record.person_id,
                        record.age,
                        record.occupation,
                        record.work_school_place,
                        -1,
                    ),
                )
            elseif record.occupation == SCHOOL
                push!(
                    household.students,
                    Person(
                        record.person_id,
                        record.age,
                        record.occupation,
                        record.work_school_place,
                        -1,
                    ),
                )
            else
                push!(
                    household.unemployeds,
                    Person(
                        record.person_id,
                        record.age,
                        record.occupation,
                        record.work_school_place,
                        -1,
                    ),
                )
            end
            household.num_components += 1
        end
        lock(mutex) do
            push!(households, household)
        end
    end
    return households
end

function generate_dataset(config_file::String, output_path::String)
    Random.seed!(0)
    config = JSON.parsefile(config_file)

    households = []
    companies = Set{Company}()
    companies_full = Set{Company}()
    schools = Dict{String,Array{Int,1}}()
    num_people = config["people"]
    num_companies = config["workplaces"]
    num_schools = sum([pair[2] for pair in config["schools"]])
	probs = Dict{Int,Float16}() # Probability dict
    info = Dict{String,Int}() # Useful information dict

    students = Set{Person}()
    employeds = Set{Person}()
    unemployeds = Set{Person}()

    min_workplace_size = 66.7 # Percentage
    id = 1
    push!(info,WORKPLACES_START_ID => id) # Save where the ids starts
    for company_class in config["companies_size"]
        num_company_to_generate = ceil(Int, num_companies * company_class[2] / 100)
        while num_company_to_generate > 0
            num_company_to_generate -= 1
            size = _get_random_size(company_class[1], config["companies_max_size"])
            company = Company(id, size)
            push!(companies, company)
			push!(probs,id => config["company_probs"][company_class[1]])
            id += 1
        end
    end

    # Get the ids for the schools
    push!(info,SCHOOLS_START_ID => id) # Save where the ids starts
    for school_class in config["schools"]
        cur_school_ids = Array{Int,1}()
        cur_school_num = school_class[2]
        for i = 1:cur_school_num
            push!(cur_school_ids, id)
			push!(probs,id => config["school_probs"][school_class[1]])
            id += 1
        end
        push!(schools, school_class[1] => cur_school_ids)
    end
    next_id = id

    id = 1
    for population_class in config["population_dist"]
        num_people_to_generate = floor(Int, num_people * population_class[2] / 100)
        println("$(population_class) -> $(num_people_to_generate)")
        while num_people_to_generate > 0
            num_people_to_generate -= 1
            age = _get_random_size(population_class[1], config["person_max_age"])
            occupation = _get_occupation_by_age(age)
            school_work_place = _get_school_work_place(
                age,
                occupation,
                companies,
                schools,
                config["students_school"],
            )

            if occupation == WORKING
                company = first(filter(c -> c.id == school_work_place, companies))
                if company.max_size == company.current_size
                    delete!(companies, company)
                    push!(companies_full, company)
                end
            end

            person = Person(id, age, occupation, school_work_place, -1)
            _add_person_to_right_set(person, students, employeds, unemployeds)

            id += 1
        end
    end

    id = next_id #Because household ids are also used in the simulation, they must be unique
    push!(info,HOUSEHOLDS_START_ID => id) # Save where the ids starts
    p_one_adult = 0.3
    p_more_adults = 0.2
    μ = 2.5 # Medium househol size
    σ = 0.5 # Minum size 1
    households_size_distribution = Normal(μ, σ)
    age_difference = 15 # max age difference between two person in a household
    remaining_people = num_people

    Adult = employeds ∪ unemployeds
    Children = students

    while remaining_people > 0
        person_added = 0
        household_size = max(1, round(Int, rand(households_size_distribution)))
        household = Household(id, household_size)

        adult_1 = rand(Adult)
        _add_person_to_right_set(
            adult_1,
            household.students,
            household.workers,
            household.unemployeds,
        )
        delete!(Adult, adult_1)
        person_added += 1

        if household_size > 1
            if !isempty(Adult) && rand() > p_one_adult
                Adult_close_age = filter(
                    person ->
                        max(CHILD_AGE, person.age - age_difference) <=
                        person.age <=
                        person.age + age_difference,
                    Adult,
                )
                adult_2 = rand(Adult_close_age)
                _add_person_to_right_set(
                    adult_2,
                    household.students,
                    household.workers,
                    household.unemployeds,
                )
                delete!(Adult, adult_2)
                person_added += 1
            end

            if !isempty(Children)
                while person_added < household_size && !isempty(Children)
                    child = rand(Children)
                    _add_person_to_right_set(
                        child,
                        household.students,
                        household.workers,
                        household.unemployeds,
                    )
                    delete!(Children, child)
                    person_added += 1
                end
            elseif !isempty(Adult)
                while person_added < household_size &&
                          rand() > p_more_adults &&
                          !isempty(Adult)
                    adult_2 = rand(Adult)
                    _add_person_to_right_set(
                        adult_2,
                        household.students,
                        household.workers,
                        household.unemployeds,
                    )
                    delete!(Adult, adult_2)
                    person_added += 1
                end
            end

        end
        remaining_people -= person_added
        household.num_components = person_added # When some sets are empty, these values can be different
        push!(probs,id => 1) # Contact prob inside household is always 1
        push!(households, household)
        id += 1
    end

    #Transports are placed at the end
    push!(info,TRANSPORTS_START_ID => id) # Save where the ids starts
    for i in id:id+config["transports"]
        push!(probs,i => config["transports_prob"])
    end

    _write_households(households, output_path)
    _write_dict(probs,output_path,"probs")
    _write_dict(info,output_path,"info")

    return households
end

max_stars_repo_name: Luigi-Crisci/Hypergraphs-Epidemics

max_issues_repo_name: Luigi-Crisci/Hypergraphs-Epidemics

max_forks_repo_name: Luigi-Crisci/Hypergraphs-Epidemics


text:
"""
A generated function, with the code that generated it
"""
struct GeneratedFunction{F}
    f::F
    code::Expr
end

#GeneratedFunction(code::Expr) = GeneratedFunction(eval(code), code)

(f::GeneratedFunction{F})(x...) where {F} = f.f(x...)

function make_tuple(args)

    if isa(args, Symbol)
        # args = [args]
        return args
    end

    length(args) == 1 && return args[1]

    return Expr(:tuple, args...)
end

function really_make_tuple(args)

    return Expr(:tuple, args...)
end

function emit_forward_code(a::Assignment)

    args = isa(a.args, Vector) ? a.args : [a.args]

    lhs = make_tuple(a.lhs)

    if a.op == :()  # empty
        args = make_tuple(args)
        :($lhs = $args)

    else
        :($lhs = $(a.op)($(args...) ) )
    end
end

function emit_backward_code(a::Assignment)

    args = isa(a.args, Vector) ? a.args : [a.args]

    return_args = [a.lhs, args...]
    rev_op = reverse_operations[a.op]  # find reverse operation

    if rev_op == :()   # empty
        args = make_tuple(args)
        lhs = make_tuple(a.lhs)
        return :($args = $lhs)
        #return :($(args...) = $(a.lhs))
    else
        rev_code = :($(rev_op)($(return_args...)))
    end

    # delete non-symbols in return args:
    for (i, arg) in enumerate(return_args)
        if !(isa(arg, Symbol))
            return_args[i] = :_
        end
    end

    return_tuple = Expr(:tuple, return_args...)  # make tuple out of array
# or: :($(return_args...),)

    return :($(return_tuple) = $(rev_code))

end

# TODO: Just pass intermediate as tuple between forward and backward for functions

function emit_forward_code(a::FunctionAssignment)
    f = a.f

    args = isa(a.args, Vector) ? a.args : [a.args]
    args_tuple = really_make_tuple(args)

    return_tuple = really_make_tuple(a.return_arguments)
    intermediate = really_make_tuple(a.intermediate)

    # Remove the following once https://github.com/JuliaLang/julia/issues/20524 fixed and replace with
    # :( ( $return_tuple, $intermediate ) = $(esc(f)).forward($args_tuple))

    temp1 = make_symbol(:z_tuple)
    temp2 = make_symbol(:z_tuple)

    return quote
        ($temp1, $temp2) = $(esc(f)).forward($args_tuple)
        $return_tuple = $temp1
        $intermediate = $temp2
    end
end

function emit_backward_code(a::FunctionAssignment)
    f = a.f

    args = isa(a.args, Vector) ? a.args : [a.args]
    args_tuple = really_make_tuple(args)

    intermediate = really_make_tuple(a.intermediate)

    return_tuple = really_make_tuple(a.return_arguments)

    :($args_tuple = $(esc(f)).backward($args_tuple, $return_tuple, $intermediate))
end

function emit_forward_code(code)  # code::Vector{Assignment})
    new_code = quote end
    new_code.args = vcat([emit_forward_code(line) for line in code])
    return new_code
end

function emit_backward_code(code) #::Vector{Assignment})
    new_code = quote end
    new_code.args = vcat([emit_backward_code(line) for line in reverse(code)])
    return new_code
end

function forward_backward(flatAST::FlatAST)

    # @show flatAST.input_variables
    # @show flatAST.intermediate
    input = collect(flatAST.input_variables)

    # @show flatAST.top

    if isa(flatAST.top, Symbol)
        output = [flatAST.top]

    elseif isa(flatAST.top, Expr) && flatAST.top.head == :tuple
        output = flatAST.top.args

    else
        output = flatAST.top
        # @show output
    end

    # @show input
    # @show flatAST.intermediate
    # @show output

    input = setdiff(input, flatAST.intermediate)  # remove local variables
    intermediate = setdiff(flatAST.intermediate, output)

    flatAST.variables = input

    forward_code = emit_forward_code(flatAST.code)
    forward = make_forward_function(input, output, intermediate, forward_code)

    # @show input
    # @show intermediate
    # @show output

    backward_code = emit_backward_code(flatAST.code)
    backward = make_backward_function(input, output, intermediate,
                                input, backward_code)

    # @show input
    # @show output
    # @show intermediate

    return (forward, backward)
end

"""
Generate code for an anonymous function with given
input arguments, output arguments, and code block.
"""
function make_forward_function(input_args, output_args, intermediate, code)

    input = really_make_tuple(input_args)  # make a tuple of the variables
    intermediate = really_make_tuple(intermediate)
    output = make_tuple(output_args)  # make a tuple of the variables

    quote
        t -> begin
                $input = t
                $code
                return ($output, $intermediate)
            end
    end
end

function make_backward_function(input1, input2, input3, output_args, code)

    input1 = really_make_tuple(input1)  # make a tuple of the variables
    input2 = really_make_tuple(input2)
    input3 = really_make_tuple(input3)
    output = really_make_tuple(output_args)  # make a tuple of the variables

    quote
        (t1, t2, t3) -> begin
            $input1 = t1
            $input2 = t2
            $input3 = t3
            $code
            return $output
        end
    end
end


text:
__precompile__()

module FastaIO

using Compat
using GZip

export
    FastaReader,
    readentry,
    rewind,
    readfasta,
    FastaWriter,
    writeentry,
    writefasta

import Base.start, Base.done, Base.next, Base.readall,
       Base.close, Base.show, Base.eof, Base.write

import Compat: String

const fasta_buffer_size = 4096

type FastaReader{T}
    # public but read-only
    num_parsed::Int        # number of parsed entries so far
    # private
    f::IO
    is_eof::Bool           # did we reach end of file?
    rbuffer::Vector{UInt8} # read buffer
    rbuf_sz::Int           # read buffer size
    rbuf_pos::Int          # read buffer cursor
    lbuffer::Vector{UInt8} # line buffer
    lbuf_sz::Int           # line buffer size
    mbuffer::Vector{UInt8} # multi-line buffer
    mbuf_sz::Int           # multi-line buffer size
    own_f::Bool
    function FastaReader(filename::AbstractString)
        fr = new(0, gzopen(filename), false,
            Array(UInt8, fasta_buffer_size), 0, 0,
            Array(UInt8, fasta_buffer_size), 0,
            Array(UInt8, fasta_buffer_size), 0,
            true)
        finalizer(fr, close)
        return fr
    end
    function FastaReader(io::IO)
        new(0, io, false,
            Array(UInt8, fasta_buffer_size), 0, 0,
            Array(UInt8, fasta_buffer_size), 0,
            Array(UInt8, fasta_buffer_size), 0,
            false)
    end
end

FastaReader(filename::AbstractString) = FastaReader{String}(filename)
FastaReader(io::IO) = FastaReader{String}(io)

function FastaReader(f::Function, filename::AbstractString, T::Type=String)
    fr = FastaReader{T}(filename)
    try
        f(fr)
    finally
        close(fr)
    end
end

close(fr::FastaReader) = fr.own_f && close(fr.f)

function rewind(fr::FastaReader)
    seek(fr.f, 0)
    fr.is_eof = false
    fr.num_parsed = 0
    fr.rbuf_sz = 0
    fr.rbuf_pos = 0
    fr.lbuf_sz = 0
    fr.mbuf_sz = 0
    return
end

read_chunk_ll(fr::FastaReader, s::GZipStream) = gzread(s, pointer(fr.rbuffer), fasta_buffer_size)
read_chunk_ll(fr::FastaReader, s::IOStream) =
    ccall(:ios_readall, UInt, (Ptr{Void}, Ptr{Void}, UInt), fr.f.ios, fr.rbuffer, fasta_buffer_size)
function read_chunk_ll(fr::FastaReader, s::IO)
    ret = 0
    while !eof(fr.f) && ret < fasta_buffer_size
        ret += 1
        fr.rbuffer[ret] = read(fr.f, UInt8)
    end
    return ret
end

function read_chunk(fr::FastaReader)
    if fr.is_eof
        return
    end
    ret = read_chunk_ll(fr, fr.f)
    ret == -1 && error("read failure")
    fr.rbuf_sz = ret
    fr.rbuf_pos = 1
    if ret == 0
        fr.is_eof = true
    end
    return
end

function readline(fr::FastaReader)
    fr.lbuf_sz = 0
    found = false
    while !fr.is_eof
        if fr.rbuf_pos == 0
            read_chunk(fr::FastaReader)
        end
        i = fr.rbuf_pos
        cr = false
        while i <= fr.rbuf_sz
            c = fr.rbuffer[i]
            if c == UInt8('\n')
                found = true
                break
            else
                cr = (c == UInt8('\r'))
            end
            i += 1
        end
        i -= 1 + cr
        chunk_len = i - fr.rbuf_pos + 1
        free_sbuf = length(fr.lbuffer) - fr.lbuf_sz
        gap = chunk_len - free_sbuf
        if gap > 0
            resize!(fr.lbuffer, length(fr.lbuffer) + gap)
        end

        #fr.lbuffer[fr.lbuf_sz + (1:chunk_len)] = fr.rbuffer[fr.rbuf_pos:i]
        copy!(fr.lbuffer, fr.lbuf_sz + 1, fr.rbuffer, fr.rbuf_pos, chunk_len)
        fr.lbuf_sz += chunk_len

        i += 2 + cr
        if i > fr.rbuf_sz
            i = 0
        end
        fr.rbuf_pos = i
        found && break
    end
    return
end

function start(fr::FastaReader)
    rewind(fr)
    readline(fr)
    if fr.lbuf_sz == 0
        error("empty FASTA file")
    end
    return
end
done(fr::FastaReader, x::Void) = fr.is_eof
function _next_step(fr::FastaReader)
    if fr.lbuffer[1] != UInt8('>')
        error("invalid FASTA file: description does not start with '>'")
    end
    if fr.lbuf_sz == 1
        error("invalid FASTA file: empty description")
    end
    name = String(fr.lbuffer[2:fr.lbuf_sz])
    isascii(name) || error("invalid non-ASCII description in FASTA file")
    fr.mbuf_sz = 0
    while true
        readline(fr)
        if fr.lbuf_sz == 0 || fr.lbuffer[1] == UInt8('>')
            break
        end
        gap = fr.lbuf_sz - (length(fr.mbuffer) - fr.mbuf_sz)
        if gap > 0
            resize!(fr.mbuffer, length(fr.mbuffer) + gap)
        end
        #fr.mbuffer[fr.mbuf_sz + (1:fr.lbuf_sz)] = fr.lbuffer[1:fr.lbuf_sz]
        copy!(fr.mbuffer, fr.mbuf_sz + 1, fr.lbuffer, 1, fr.lbuf_sz)
        fr.mbuf_sz += fr.lbuf_sz
    end
    return name
end
function _next(fr::FastaReader{Vector{UInt8}})
    name = _next_step(fr)
    fr.num_parsed += 1
    return (name, fr.mbuffer[1:fr.mbuf_sz])
end
function _next(fr::FastaReader{String})
    name = _next_step(fr)
    out_str = ccall(:jl_pchar_to_string, Ref{String}, (Ptr{UInt8},Int), fr.mbuffer, fr.mbuf_sz)
    fr.num_parsed += 1
    return (name, out_str)
end
function _next{T}(fr::FastaReader{T})
    name = _next_step(fr)
    fr.num_parsed += 1
    return (name, T(fr.mbuffer[1:fr.mbuf_sz]))
end

next(fr::FastaReader, x::Void) = (_next(fr), nothing)

function readall(fr::FastaReader)
    ret = Any[]
    for item in fr
        push!(ret, item)
    end
    return ret
end

function readentry(fr::FastaReader)
    fr.is_eof && throw(EOFError())
    if fr.num_parsed == 0
        readline(fr)
        if fr.lbuf_sz == 0
            error("empty FASTA file")
        end
    end
    item, _ = next(fr, nothing)
    return item
end

eof(fr::FastaReader) = fr.is_eof

function show{T}(io::IO, fr::FastaReader{T})
    print(io, "FastaReader(input=\"$(fr.f)\", out_type=$T, num_parsed=$(fr.num_parsed), eof=$(fr.is_eof))")
end

function readfasta(filename::AbstractString, T::Type=String)
    FastaReader(filename, T) do fr
        readall(fr)
    end
end
readfasta(io::IO, T::Type=String) = readall(FastaReader{T}(io))

type FastaWriter
    f::IO
    in_seq::Bool
    entry_chars::Int
    desc_chars::Int
    parsed_nl::Bool
    pos::Int
    entry::Int
    own_f::Bool
    at_start::Bool
    function FastaWriter(io::IO)
        fw = new(io, false, 0, 0, false, 0, 1, false, true)
        finalizer(fw, close)
        return fw
    end
    function FastaWriter(filename::AbstractString, mode::AbstractString = "w")
        if endswith(filename, ".gz")
            of = gzopen
        else
            of = open
        end
        fw = new(of(filename, mode), false, 0, 0, false, 0, 1, true, true)
        finalizer(fw, close)
        return fw
    end
end

FastaWriter() = FastaWriter(STDOUT)

function FastaWriter(f::Function, args...)
    fw = FastaWriter(args...)
    try
        f(fw)
    finally
        close(fw)
    end
end

function write(fw::FastaWriter, c)
    ch = convert(Char, c)
    isascii(ch) || error("invalid (non-ASCII) character: $c (entry $(fw.entry) of FASTA input)")
    if ch == '\n' && !fw.at_start
        fw.parsed_nl = true
        if !fw.in_seq
            fw.desc_chars == 1 && error("empty description (entry $(fw.entry) of FASTA input")
            write(fw.f, '\n')
            fw.pos = 0
            fw.in_seq = true
        end
    end
    isspace(ch) && (fw.at_start || fw.in_seq || fw.desc_chars <= 1) && return
    fw.at_start && ch != '>' && error("no desctiption given (entry $(fw.entry) of FASTA input")
    fw.at_start = false
    if fw.parsed_nl
        @assert fw.in_seq
        if ch == '>'
            fw.entry_chars > 0 || error("description must span a single line (entry $(fw.entry) of FASTA input)")
            write(fw.f, '\n')
            fw.in_seq = false
            fw.pos = 0
            fw.entry += 1
            fw.entry_chars = 0
            fw.desc_chars = 0
        end
    elseif fw.in_seq && ch == '>'
        error("character '>' not allowed in sequence data (entry $(fw.entry) of FASTA input)")
    end
    if fw.pos == 80
        if !fw.in_seq
            warn("description line longer than 80 characters (entry $(fw.entry) of FASTA input)")
        else
            write(fw.f, '\n')
            fw.pos = 0
        end
    end
    write(fw.f, ch)
    fw.pos += 1
    if fw.in_seq
        fw.entry_chars += 1
    else
        fw.desc_chars += 1
    end
    fw.parsed_nl = false
    return
end

function write(fw::FastaWriter, s::Vector)
    for c in s
        write(fw, c)
    end
end

function write(fw::FastaWriter, s::AbstractString)
    for c in s
        write(fw, c)
    end
    write(fw, '\n')
end

function writeentry(fw::FastaWriter, desc::AbstractString, seq)
    !fw.at_start && write(fw, '\n')
    desc = strip(String(desc))
    isascii(desc) || error("description must be ASCCII (entry $(fw.entry+1) of FASTA input)")
    if search(desc, '\n') != 0
        error("newlines are not allowed within description (entry $(fw.entry+1) of FASTA input)")
    end
    write(fw, '>')
    write(fw, desc)
    write(fw, '\n')
    #write(fw, seq)
    #write(fw, '\n')
    fw.entry_chars = writefastaseq(fw.f, seq, fw.entry, false)
    fw.in_seq = true
    fw.parsed_nl = false
    fw.pos = 0
    fw.entry_chars > 0 || error("empty sequence data (entry $(fw.entry) of FASTA input)")
    return
end

function close(fw::FastaWriter)
    try
        write(fw.f, '\n')
        flush(fw.f)
    catch err
        isa(err, EOFError) || rethrow(err)
    end
    fw.pos = 0
    fw.parsed_nl = true
    fw.own_f && close(fw.f)
    return
end

function show(io::IO, fw::FastaWriter)
    print(io, "FastaWriter(input=\"$(fw.f)\", entry=$(fw.entry)")
end

function writefastaseq(io::IO, seq, entry::Int, nl::Bool = true)
    i = 0
    entry_chars = 0
    for c in seq
        if i == 80
            write(io, '\n')
            i = 0
        end
        ch = convert(Char, c)
        isascii(ch) || error("invalid (non-ASCII) character: $c (entry $entry of FASTA input)")
        isspace(ch) && continue
        ch != '>' || error("character '>' not allowed in sequence data (entry $entry of FASTA input)")
        write(io, ch)
        i += 1
        entry_chars += 1
    end
    nl && write(io, '\n')
    return entry_chars
end

function writefasta(io::IO, data)
    entry = 0
    for (desc, seq) in data
        entry += 1
        desc = strip(String(desc))
        isascii(desc) || error("description must be ASCCII (entry $entry of FASTA input)")
        if isempty(desc)
            error("empty description (entry $entry of FASTA input")
        end
        if search(desc, '\n') != 0
            error("newlines are not allowed within description (entry $entry of FASTA input)")
        end
        if length(desc) > 79
            warn("description line longer than 80 characters (entry $entry of FASTA input)")
        end
        println(io, ">", desc)
        entry_chars = writefastaseq(io, seq, entry)
        entry_chars > 0 || error("empty sequence data (entry $entry of FASTA input)")
    end
end
writefasta(data) = writefasta(STDOUT, data)

function writefasta(filename::AbstractString, data, mode::AbstractString = "w")
    if endswith(filename, ".gz")
        of = gzopen
    else
        of = open
    end
    of(filename, mode) do f
        writefasta(f, data)
    end
end

end


text:
module DGMethods

using MPI
using StaticArrays
using DocStringExtensions
using KernelAbstractions
using KernelAbstractions.Extras: @unroll
using ..MPIStateArrays
using ..Mesh.Grids
using ..Mesh.Topologies
using ..VariableTemplates
using ..Courant
using ..BalanceLaws:
    BalanceLaw,
    AbstractStateType,
    Prognostic,
    Auxiliary,
    Gradient,
    GradientFlux,
    GradientLaplacian,
    Hyperdiffusive,
    UpwardIntegrals,
    DownwardIntegrals,
    vars_state,
    number_states

import ..BalanceLaws:
    BalanceLaw,
    init_state_prognostic!,
    init_state_auxiliary!,
    flux_first_order!,
    flux_second_order!,
    compute_gradient_flux!,
    compute_gradient_argument!,
    source!,
    transform_post_gradient_laplacian!,
    wavespeed,
    boundary_state!,
    update_auxiliary_state!,
    update_auxiliary_state_gradient!,
    nodal_update_auxiliary_state!,
    nodal_init_state_auxiliary!,
    integral_load_auxiliary_state!,
    integral_set_auxiliary_state!,
    indefinite_stack_integral!,
    reverse_indefinite_stack_integral!,
    reverse_integral_load_auxiliary_state!,
    reverse_integral_set_auxiliary_state!

export DGModel,
    init_ode_state,
    restart_ode_state,
    restart_auxiliary_state,
    basic_grid_info,
    init_state_auxiliary!,
    continuous_field_gradient!,
    courant

include("NumericalFluxes.jl")
include("DGModel.jl")
include("DGModel_kernels.jl")
include("create_states.jl")

"""
    calculate_dt(dg, model, Q, Courant_number, direction, t)

For a given model, compute a time step satisying the nondiffusive Courant number
`Courant_number`
"""
function calculate_dt(dg, model, Q, Courant_number, t, direction)
    Δt = one(eltype(Q))
    CFL = courant(nondiffusive_courant, dg, model, Q, Δt, t, direction)
    return Courant_number / CFL
end

end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
#miscellenious helper functions and constants

import Base.@math_const

Q95 = sqrt(2.)*erfinv(0.95)

range(x) = (minimum(x), maximum(x))

#95% normal quantile

eps2 = sqrt(eps())

intervalgaussian(U, a, b) = broadcast(intervalgaussian, U, a, b)

function intervalgaussian (U::Real,a::Real, b::Real)

    a, b = min(a,b), max(a,b)

    c =  erf(a/sqrt(2.))

    d =  erf(b/sqrt(2.))

    sqrt(2.)*erfinv(c + U.*(d-c))

end

function scalar(x)

assert(length(x) == 1)

x[1]

end

function norma!(x)

 minx = min(x)

 x = (x - minx) / (max(x)-minx)

 x

end

function cut(x, a, b)

    warn("cut(): better use clamp()")

    clamp(x, a,b)

end

#compute normal confidence interval for Monte Carlo estimate, precision aware rounding

function mc(M)

  m = mean(M)

  ste = std(M)/sqrt(length(M))

  round([m, Q95*ste], int(2-log(10,Q95*ste)))

end

#for sequential estimates, it suffices to provide the running sum Z, running sum of squares Z2 and the number of observations k

va(k, Z, Z2) =  Z2/k-(Z/k)^2
ste(k, Z, Z2) = sqrt(Z2/k^2-(Z/k)^2/k)

function mc2(k, Z::Real, Z2::Real, mode = :ci)
    m = Z/k
    if mode ==: ci
        F = Q95 else F = 1.
    end

    try
        ste = sqrt(Z2/k-m.^2)/sqrt(k)
        res = (m, F*ste)
        return res
    catch
        return (m, NaN.*m)
    end
end

function mc2(k, Z::Vector, Z2::Vector, mode = :ci)
   map((z,z2) -> mc2(k, z, z2, mode), Z, Z2)
end
function mc2str(k, Z, Z2)
    m, ste = mc2(k, Z, Z2, :se)
    m, ci = round([m, Q95*ste], int(2-log(10,Q95*ste)))
    "$m+-$(ci)ci "
end
function mc2str(k, Z::Vector, Z2::Vector)
   reduce(*,map((z,z2) -> mc2str(k, z, z2), Z, Z2))
end

function isignif_og(x, digits, base)

    if base == 10

        ifloor(log10(abs(x)) - digits + 1.)

    elseif base == 2

        ifloor(log2(abs(x)) - digits + 1.)

    else

        ifloor(log2(abs(x))/log2(base) - digits + 1.)

    end

end

function roundste(out::IOBuffer, m, ste::FloatingPoint)

    if

        !isfinite(m) print(out, m)

        return out

    end

    if isnan(ste)

        print(out, m, " ± NaN (se)")

        return out

    elseif isinf(ste)

        print(out, m, " ± Inf (se)")

        return out

    elseif ste < eps(m)

        print(out, m)

        return out

    end

    assert(ste >= 0.)

    og = max(isignif_og(m, 3, 10) - isignif_og(ste, 3, 10)+3,0)

    grisu(out, m, Base.Grisu.PRECISION, og)

    print(out, " ± ")

    grisu(out, ste, Base.Grisu.PRECISION, 3)

    print(out, " (se)")

out

end

function roundste(out::IOBuffer, r)

    m, ste = r

    roundste(out, m[1], ste[1])

    for i in 2:length(m)

        print(", ")

        roundste(out, m[i], ste[i])

    end

end

function strmc2(k, Z, Z2)

    s = IOBuffer(true, true)

    truncate(s,0)

    roundste(s, mc2(k, Z, Z2, false))

    takebuf_string(s)

end

#stdv(k, Z, Z2) = va(k, Z, Z2)/k

# X = [ZL, L]

# X2 = X * X'

function selfadjmc(k, X, X2)

    m = X[1]/X[2]

    v = (va(k, X[1], X2[1,1]) - 2m*(X2[1,2]/k - X[1]*X[2]/k^2) + m^2*va(k, X[2], X2[2,2]))/k

#   println([va(k, X[1], X2[1,1]), - 2m*(X2[1,2]/k - X[1]*X[2]/k^2),m^2*va(k, X[2], X2[2,2])])

    try

        ste = sqrt(v)

        res = [m, ste]

        return res

    catch

        return [m, NaN.*m]

    end

end

function grisu(io::IO, x::FloatingPoint, mode::Int32, n::Int)

    if isnan(x) return write(io, "NaN"); end

    if isinf(x)

        if x < 0 write(io,'-') end

        write(io, "Inf")

        return

    end

    Base.Grisu.@grisu_ccall(x, mode, n)

    pdigits = pointer(Base.Grisu.DIGITS)

    neg = Base.Grisu.NEG[1]

    len = int(Base.Grisu.LEN[1])

    pt  = int(Base.Grisu.POINT[1])

#    if mode == PRECISION

#        while len > 1 && DIGITS[len] == '0'

#            len -= 1

#        end

#    end

    if neg write(io,'-') end

    if pt <= -4 || pt > 6 # .00001 to 100000.

        # => #.#######e###

        write(io, pdigits, 1)

        write(io, '.')

        if len > 1

            write(io, pdigits+1, len-1)

        else

            write(io, '0')

        end

        write(io, 'e')

        write(io, dec(pt-1))

        return

    elseif pt <= 0

        # => 0.00########

        write(io, "0.")

        while pt < 0

            write(io, '0')

            pt += 1

        end

        write(io, pdigits, len)

    elseif pt >= len

        # => ########00.0

        write(io, pdigits, len)

        while pt > len

            write(io, '0')

            len += 1

        end

        write(io, ".0")

    else # => ####.####

        write(io, pdigits, pt)

        write(io, '.')

        write(io, pdigits+pt, len-pt)

    end

    nothing

end

#plot (2d) sample path

function pl(x::Array{ FloatingPoint,2})

    p = FramedPlot()

    #compute range

    R1 = range(x[1,:])

    R2 = range(x[2,:])

    #widen a bit

    R1 = (R1[1] -0.1*(R1[2]-R1[1]), R1[2] + 0.1*(R1[2]-R1[1]))

    R2 = (R2[1] -0.1*(R2[2]-R2[1]), R2[2] + 0.1*(R2[2]-R2[1]))

    setattr(p, "xrange", R1)

    setattr(p, "yrange", R2)

    add(p, Curve(x[1,:],x[2,:]))

    Winston.display(p)

    p

end


text:
module RFCVGT
__precompile__()
# package code goes here
using DataFrames,DecisionTree,MLBase

export build_forest,apply_forest,gridtune
export searchparm
export Unievalfun,Estfun,Evalfun,Gridsearch
export unievalfun,estfun,evalfun,gridsearch

include("estfun.jl")
include("unievalfun.jl")
include("evalfun.jl")
include("rfcvgtfun.jl")
end # module

max_stars_repo_name: wakakusa/RFCVGT

max_issues_repo_name: wakakusa/RFCVGT

max_forks_repo_name: wakakusa/RFCVGT


text:
# Reversed indices
struct ReverseIndices{I, Inds <: AbstractIndices{I}} <: AbstractIndices{I}
    inds::Inds
end

Base.parent(inds::ReverseIndices) = getfield(inds, :inds)
@propagate_inbounds Base.in(inds::ReverseIndices{I}, i::I) where {I} = in(parent(inds), i)
Base.length(inds::ReverseIndices) = length(parent(inds))
Base.IteratorSize(inds::ReverseIndices) = Base.IteratorSize(parent(inds))

istokenizable(inds::ReverseIndices) = istokenizable(parent(inds))
tokentype(inds::ReverseIndices) = tokentype(parent(inds))
@propagate_inbounds gettoken(inds::ReverseIndices{I}, i::I) where {I} = gettoken(parent(inds), i)
@propagate_inbounds gettokenvalue(inds::ReverseIndices, t) = gettokenvalue(parent(inds), t)

@propagate_inbounds iteratetoken(inds::ReverseIndices, s...) = iteratetoken_reverse(parent(inds), s...)
@propagate_inbounds iteratetoken_reverse(inds::ReverseIndices, s...) = iteratetoken(parent(inds), s...)

Iterators.reverse(inds::AbstractIndices) = ReverseIndices{eltype(inds), typeof(inds)}(inds)
Iterators.reverse(inds::ReverseIndices) = parent(inds)

Base.reverse(inds::AbstractIndices) = copy(Iterators.reverse(inds))

Base.empty(inds::ReverseIndices, ::Type{I}) where {I} = empty(parent(inds), I)
Base.empty(inds::ReverseIndices, ::Type{I}, ::Type{T}) where {I, T} = empty(parent(inds), I, T)

Base.similar(inds::ReverseIndices, ::Type{T}) where {T} = Iterators.reverse(similar(parent(inds), T))

# Reversed dictionary
struct ReverseDictionary{I, T, Dict <: AbstractDictionary{I, T}} <: AbstractDictionary{I, T}
    dict::Dict
end

Base.parent(dict::ReverseDictionary) = getfield(dict, :dict)

@inline function Base.keys(dict::ReverseDictionary)
    Iterators.reverse(keys(parent(dict)))
end

issetable(dict::ReverseDictionary) = issetable(parent(dict))

Base.isassigned(dict::ReverseDictionary{I}, i::I) where {I} = isassigned(parent(dict), i)
@propagate_inbounds Base.getindex(dict::ReverseDictionary{I}, i::I) where {I} = getindex(parent(dict), i)
@propagate_inbounds Base.setindex!(dict::ReverseDictionary{I, T}, i::I, v::T) where {I, T} = setindex!(parent(dict), v, i)

istokenassigned(dict::ReverseDictionary, t) = istokenassigned(parent(dict), t)
@propagate_inbounds gettokenvalue(dict::ReverseDictionary, t) = gettokenvalue(parent(dict), t)
@propagate_inbounds settokenvalue!(dict::ReverseDictionary{<:Any, T}, t, v::T) where {T} = settokenvalue!(parent(dict), t, v)

Iterators.reverse(dict::AbstractDictionary) = ReverseDictionary{keytype(dict), eltype(dict), typeof(dict)}(dict)
Iterators.reverse(dict::ReverseDictionary) = parent(dict)

function Base.reverse(dict::AbstractDictionary)
    out = similar(reverse(keys(dict)), eltype(dict))
    @inbounds copyto!(out, Iterators.reverse(dict))
    return out
end

Base.empty(inds::ReverseDictionary, ::Type{I}) where {I} = empty(parent(inds), I)
Base.empty(inds::ReverseDictionary, ::Type{I}, ::Type{T}) where {I, T} = empty(parent(inds), I, T)


text:
using ArgumentModes
using Test

@test Mode(:a, :b => 1, :c => "sdf") isa
  Mode[:a, :b => Int, :c => Union{Float64, AbstractString}]
@test !(Mode(:a, :b => 1, :c => 2) isa
  Mode[:a, :b => Int, :c => Union{Float64, AbstractString}])

@test repr(Mode(:a => 1., :b => "1", :c , :d => 2)) ==
  "Mode[==, :a => Float64, :b => String, :c, :d => Int64]"*
  "((a = 1.0, b = \"1\", c = nothing, d = 2))"
@test repr(Mode[:a => Int, :b => Float64, :c , :d => String]) ==
  "Mode[:a => Int64, :b => Float64, :c, :d => String]"
@test repr(Mode) == "Mode"

@test Mode[:a, :b => Int, :c => Float64] ==
  Mode[:c => Float64, :a, :b => Int] == Mode[:c => Float64, :b => Int, :a]

@test typeof(Mode(:a, :b => 1, :c => "sdf")) ==
  Mode[==, :b => Int, :c => String, :a] ==
  typeof(Mode(:b => 1, :a, :c => "sdf")) ==
  Mode[==, :a, :c => String,:b => Int]

f(x::Mode[:mode1]) = 1
f(x::Mode[:mode2]) = 2
f(x::Mode[:mode1, :mode2]) = 3
f(x::Mode[:mode1, :mode2, :mode3]) = 4

@test f(Mode(:mode1)) == 1
@test f(Mode(:mode2)) == 2
@test f(Mode(:mode1, :mode2)) == 3
@test f(Mode(:mode3)) == 4
@test f(Mode(:mode1, :mode3)) == 4

@test Mode(:a) == Mode(:a => nothing)
@test (Mode(:a) ~ Mode(:b)) == Mode(:a, :b)
@test (Mode(:a)=>3 ~ Mode(:b)) == Mode(:a=>3, :b)
@test_throws ArgumentError Mode(:a)=>3 ~ Mode(:a)

@test keys(Mode(:a, :b=> 3)) == (:a, :b)
@test values(Mode(:a, :b=> 3)) == (nothing, 3)
@test collect(pairs(Mode(:a, :b=> 3))) == [:a=> nothing, :b=> 3]

@test !(Mode(:a => 5) isa Mode[:a => Float64])
@test Mode(:a => 5) isa Mode[:a => Real]
@test Mode(:a => 5) isa Mode[:a => Real, :b => Int]

@test Mode(:a, :b=> 3)[:a, :b] == (nothing, 3)
@test Mode(:a, :b=> 3)[:b] == 3
@test_throws ArgumentError Mode(:a, :b=> 3)[]
@test Mode(:b=> 3)[] == 3

m = Mode(:a, :b=> 2, :c=> 3.)
x = 4
@test !checkmode(x, :a)
@test checkmode(m, :a, :b)
@test !checkmode(m, :a, :b, :d)
@test checkmode(m, &, :a, :b)
@test !checkmode(m, &, :a, :b, :d)
@test checkmode(m, |, :a, :b, :d)
@test checkmode(m, |, :a, :b)
@test !checkmode(m, |, :d, :e)
@test !checkmode(m, ==, :a, :b)
@test !checkmode(m, ==, :a, :b, :c, :d)
@test checkmode(m, ==, :a, :b, :c)
@test !checkmode(m, !, :a, :b, :c)
@test !checkmode(m, !, :a, :d, :e)
@test checkmode(m, !, :d, :e)

@test checkmode(m, :b, :c) do  x,y; x,y   end == (2,3.)
@test checkmode(m, :b, :d) do  x,y; x,y   end === nothing
@test checkmode(m, &, :b, :c) do  x,y; x,y   end == (2,3.)
@test checkmode(m, &, :b, :d) do  x,y; x,y   end === nothing
@test checkmode(m, :b) do  x;  x   end == 2
@test checkmode(m, !, :d, :e) do;  1   end == 1
@test checkmode(m, !, :c, :e) do;  1   end === nothing
@test checkmode(m, |, :c, :e) do;  1   end == 1
@test checkmode(m, |, :d, :e) do;  1   end === nothing
@test checkmode(m, ==, :a, :b, :c) do x,y,z;  x,y,z   end == (nothing, 2, 3.)
@test checkmode(m, ==, :a, :b, :c, :d) do x,y,z;  1   end === nothing
@test checkmode(m, ==, :a, :b, :d) do x,y,z;  1   end === nothing
@test checkmode(m, ==, :a, :b) do x,y,z;  1   end === nothing

f2(x::Mode[:m1, :m2 => Int]) = 1
f2(x::Mode[:m1, :m2 => String]) = 2

@test_throws MethodError f2(Mode(:m1) ~ Mode(:m2))
@test f2(Mode(:m1) ~ Mode(:m2)=> 1) == 1
@test f2(Mode(:m1) ~ Mode(:m2)=> "1") == 2
@test_throws MethodError f2(Mode(:m1)=> 1 ~ Mode(:m2)=> 1)
@test_throws MethodError f2(Mode(:m1) ~ Mode(:m3))


text:
# fitnessPath: Synapse id, file or folder (should point to MyProject/Metadata).
#              If it's a folder, all .csv files in the folder will be checked.
function getsamplefitness(syn, fitnessPath)
    fitness = getsamplemetadata(syn, fitnessPath) # hijack metadata loading which is very similar
    rename!(fitness, :MetadataID, :FitnessID)
end

function appendfitness!(metadata::DataFrame, fitnessTable::DataFrame; allowInexactMatches=true)
    ids = metadata[:SampleID]
    @assert length(unique(ids)) == length(ids)
    ids = dropna(fitnessTable[:SampleID])
    @assert length(unique(ids)) == length(ids)
    @assert eltype(metadata[:Passage])<:Integer
    @assert eltype(fitnessTable[:Passage])<:Integer

    fitnessTable = copy(fitnessTable) # create copy since we will modify the table
    fitness = hcat(fitnessTable[:A],fitnessTable[:B],fitnessTable[:C])  # create Nx3 DataArray
    fitness = map(i->mean(dropna(fitness[i,:][:])), 1:size(fitness,1)); # take mean of each row after dropping NAs
    fitnessTable[:Fitness] = fitness
    fitnessTable = fitnessTable[.~isnan.(fitnessTable[:Fitness]), :] # drop measurements where 3 fitness values for a sample are NAs

    # create new columns
    metadata[:Fitness] = 0.0
    metadata[:FitnessID] = ""

    # find exact SampleID matches between metadata and fitnessTable
    matches = indexin(metadata[:SampleID], map(string,fitnessTable[:SampleID])) # converting to string converts NA->"NA" which can be handled by indexin
    hasMatches = matches.!=0
    matches = matches[hasMatches]

    metadata[hasMatches, :Fitness]   = fitnessTable[matches, :Fitness]
    metadata[hasMatches, :FitnessID] = fitnessTable[matches, :FitnessID]

    if allowInexactMatches
        fitnessTable = fitnessTable[setdiff(1:end,matches),:] # remove those fitness entries that have been used
        used = falses(size(fitnessTable,1)) # but from now on, keep track of which ones we have used rather than removing rows in the loop

        for i in find(ismissing.(metadata[:Fitness])) # for those samples we don't have a fitness value for
            if ismissing(metadata[i,:Reference]) || ismissing(metadata[i,:Mutagen]) || ismissing(metadata[i,:Dose]) || ismissing(metadata[i,:Passage])
                continue # skip samples with incomplete metadata
            end

            mask = (metadata[i,:Reference] .== fitnessTable[:Virus]) .&
                   (metadata[i,:Mutagen]   .== fitnessTable[:Mutagen]) .&
                   (metadata[i,:Dose]      .== fitnessTable[:Dose]) .&
                   (metadata[i,:Passage]   .== fitnessTable[:Passage]) .&
                   .~used

            ind = findfirst(mask) # if there are multiple matches, use the first one and leave others to match other samples
            ind==0 && continue # no match

            metadata[i,:Fitness]   = fitnessTable[ind,:Fitness]
            metadata[i,:FitnessID] = fitnessTable[ind,:FitnessID]
            used[ind] = true
        end

        # fitnessTable = fitnessTable[.~used,:]
    end

    metadata
end


text:
using Test

@testset "Classes" begin
    include("test_classes.jl")
end


text:
function get_callback_trace(trace, trace_p, show_every = 1)
    function callback(k::Int, dx, dlogp, log_p, p)
        trace[:f][k] = log_p
        trace[:f_abs][k] = dlogp
        trace[:f_rel][k] = abs(dlogp)/abs(log_p)

        trace[:x_abs][k] = dx

        if mod(k,show_every) == 0
            time = Dates.format(Dates.now(), "yyyy-mm-dd HH:MM:SS")
            @printf("%20s %7d logp %.16e dlogp % .8e ||dx|| %.8e\n",
                    time, k, log_p, dlogp, trace[:x_abs][k])

            # C = 0.0
            # for k = 1:8
            #     C += tr(T[k+8])
            # end
            # C /= (3*8)
            # @printf(" C_norm %.8e\n", C)
            push!(trace_p[:r], p[:r])
            push!(trace_p[:T], p[:T])
            push!(trace_p[:b], p[:b])
            push!(trace_p[:k], k)
        end
    end
    return callback
end

r0_m = hcat(r0...) |> Array
T0_m = cat(T0...;dims = 3) |> Array
b0_m = hcat(b0...) |> Array
y_m = reshape(vcat((y .|> Array )...), :,size(y,2))
Q_inv_m = cat((Q_inv .|> Array)...;dims = (1,2))

show_every = 1
# N_bcd  = 10
tol_dlogp = 1.0e-9

trace = Dict(
    :f => zeros(N_bcd),
    :f_abs => zeros(N_bcd),
    :f_rel => zeros(N_bcd),
    :x_abs => zeros(N_bcd)
)
trace_p = Dict(:r => [], :T => [], :b => [], :k => [])
r_hat_m = copy(r0_m)
T_hat_m = copy(T0_m)
b_hat_m = copy(b0_m)

theta_m, eta_m, conv = bcd!(
    r_hat_m, T_hat_m, b_hat_m, y_m, Q_inv_m, Ns, g_mag,
    1.0e-12, 100, 1.0e-13, 200, tol_dlogp, 1.0-14, N_bcd;
    callback = get_callback_trace(
        trace, trace_p, show_every
    ),
    # verbose = true
)

# end


text:
using Pkg
Pkg.activate("../")

using StatsPlots, NearestNeighbors, SCARGC

function nearestNeighborStatic(dataset::Array{T, 2} where {T<:Number}, percentTraining::Float64)
    labeledData, labeledDataLabels, streamData, streamLabels, _ = SCARGC.fitData(dataset, percentTraining)

    tree = KDTree(permutedims(labeledData))
    accuracyVector = zeros(size(streamData, 1))

    for stream = 1:size(streamData, 1)
        index, _ = knn(tree, streamData[stream, :], 1, true)

        if (streamLabels[stream] == labeledDataLabels[index[1]])
            global accuracyVector[stream] = 1
        end
    end

    finalAccuracy = (sum(accuracyVector)/size(streamData, 1)) * 100

    println(finalAccuracy)

    return finalAccuracy
end

function nearestNeighborSliding(dataset::Array{T, 2} where {T<:Number}, percentTraining::Float64)
    labeledData, labeledDataLabels, streamData, streamLabels, _ = SCARGC.fitData(dataset, percentTraining)

    tree = KDTree(permutedims(labeledData))
    accuracyVector = zeros(size(streamData, 1))

    windowSize = 50
    data = zeros(50, size(streamData, 2) + 1)
    dataStreamIndex = 1

    for stream = 1:size(streamData, 1)
        index, _ = knn(tree, streamData[stream, :], 1, true)

        if (streamLabels[stream] == labeledDataLabels[index[1]])
            global accuracyVector[stream] = 1
        end

        data[dataStreamIndex, 1:end-1] = streamData[stream, :]
        data[dataStreamIndex, end] = streamLabels[stream]
        dataStreamIndex += 1

        if stream % windowSize == 0
            labeledData = data[:, 1:end-1]
            labeledDataLabels = data[:, end]
            data = zeros(50, size(streamData, 2) + 1)
            dataStreamIndex = 1
            tree = KDTree(permutedims(labeledData))
        end
    end

    finalAccuracy = (sum(accuracyVector)/size(streamData, 1)) * 100

    println(finalAccuracy)

    return finalAccuracy
end

function nearestNeighborLandmark(dataset::Array{T, 2} where {T<:Number}, percentTraining::Float64)
    labeledData, labeledDataLabels, streamData, streamLabels, _ = SCARGC.fitData(dataset, percentTraining)

    tree = KDTree(permutedims(labeledData))
    accuracyVector = zeros(size(streamData, 1))

    windowSize = 50
    maxWindowSize = 100

    data = zeros(100, size(streamData, 2) + 1)
    dataStreamIndex = 1

    for stream = 1:size(streamData, 1)
        index, _ = knn(tree, streamData[stream, :], 1, true)

        if (streamLabels[stream] == labeledDataLabels[index[1]])
            global accuracyVector[stream] = 1
        end

        data[dataStreamIndex, 1:end-1] = streamData[stream, :]
        data[dataStreamIndex, end] = streamLabels[stream]
        dataStreamIndex += 1

        if dataStreamIndex == maxWindowSize
            labeledData = data[1:windowSize, 1:end-1]
            labeledDataLabels = data[1:windowSize, end]
            tree = KDTree(permutedims(labeledData))

            data[1:windowSize, :] = data[windowSize+1:end, :]
            dataStreamIndex = 50
        end
    end

    finalAccuracy = (sum(accuracyVector)/size(streamData, 1)) * 100

    println(finalAccuracy)

    return finalAccuracy
end

path = "../src/datasets/synthetic/"

datasets = [
    "1CDT.txt", "2CDT.txt", "1CHT.txt", "2CHT.txt", "4CR.txt", "4CRE-V1.txt",
    "4CRE-V2.txt", "5CVT.txt", "1CSurr.txt", "4CE1CF.txt", "FG_2C_2D.txt", "UG_2C_2D.txt",
    "UG_2C_3D.txt", "UG_2C_5D.txt", "MG_2C_2D.txt", "GEARS_2C_2D.txt"
]

K_values = [2, 2, 2, 2, 4, 4, 4, 5, 4, 5, 4, 2, 2, 2, 4, 2]

rowCounts = [
    16000, 16000, 16000, 16000, 144400, 125000, 183000, 24000, 55283,
    173250, 200000, 100000, 200000, 200000, 200000, 200000
]

featureCount = [2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 5, 2, 2]

poolValues = [300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300]

trainingValues = [
    0.3125, 0.3125, 0.3125, 0.3125, 0.0346, 0.04, 0.0273, 0.125,
    0.0904, 0.0288, 0.05, 0.025, 0.025, 0.025, 0.025, 0.025
]

x = [
    "1CDT", "2CDT", "1CHT", "2CHT", "4CR", "4CRE-V1", "4CRE-V2", "5CVT",
    "1CSurr", "4CE1CF", "FG_2C_2D", "UG_2C_2D", "UG_2C_3D", "UG_2C_5D",
    "MG_2C_2D", "GEARS_2C_2D"
]

accuracies = zeros(16, 4)

for i = 1:4
    data = SCARGC.extractValuesFromFile(path * datasets[i], rowCounts[i], featureCount[i] + 1)
    _, accuracy = SCARGC.scargc_1NN(data, trainingValues[i], poolValues[i], K_values[i])
    ac = nearestNeighborStatic(data, trainingValues[i])
    ac2 = nearestNeighborSliding(data, trainingValues[i])
    ac3 = nearestNeighborLandmark(data, trainingValues[i])
    accuracies[i, 1] = accuracy
    accuracies[i, 2] = ac
    accuracies[i, 3] = ac2
    accuracies[i, 4] = ac3
end

for i = 5:7
    data = SCARGC.extractValuesFromFile(path * datasets[i], rowCounts[i], featureCount[i] + 1)
    _, accuracy = SCARGC.scargc_1NN(data, trainingValues[i], poolValues[i], K_values[i])
    ac = nearestNeighborStatic(data, trainingValues[i])
    ac2 = nearestNeighborSliding(data, trainingValues[i])
    ac3 = nearestNeighborLandmark(data, trainingValues[i])
    accuracies[i, 1] = accuracy
    accuracies[i, 2] = ac
    accuracies[i, 3] = ac2
    accuracies[i, 4] = ac3
end

for i = 8:10
    data = SCARGC.extractValuesFromFile(path * datasets[i], rowCounts[i], featureCount[i] + 1)
    _, accuracy = SCARGC.scargc_1NN(data, trainingValues[i], poolValues[i], K_values[i])
    ac = nearestNeighborStatic(data, trainingValues[i])
    ac2 = nearestNeighborSliding(data, trainingValues[i])
    ac3 = nearestNeighborLandmark(data, trainingValues[i])
    accuracies[i, 1] = accuracy
    accuracies[i, 2] = ac
    accuracies[i, 3] = ac2
    accuracies[i, 4] = ac3
end

for i = 11:13
    data = SCARGC.extractValuesFromFile(path * datasets[i], rowCounts[i], featureCount[i] + 1)
    _, accuracy = SCARGC.scargc_1NN(data, trainingValues[i], poolValues[i], K_values[i])
    ac = nearestNeighborStatic(data, trainingValues[i])
    ac2 = nearestNeighborSliding(data, trainingValues[i])
    ac3 = nearestNeighborLandmark(data, trainingValues[i])
    accuracies[i, 1] = accuracy
    accuracies[i, 2] = ac
    accuracies[i, 3] = ac2
    accuracies[i, 4] = ac3
end

for i = 14:16
    data = SCARGC.extractValuesFromFile(path * datasets[i], rowCounts[i], featureCount[i] + 1)
    _, accuracy = SCARGC.scargc_1NN(data, trainingValues[i], poolValues[i], K_values[i])
    ac = nearestNeighborStatic(data, trainingValues[i])
    ac2 = nearestNeighborSliding(data, trainingValues[i])
    ac3 = nearestNeighborLandmark(data, trainingValues[i])
    accuracies[i, 1] = accuracy
    accuracies[i, 2] = ac
    accuracies[i, 3] = ac2
    accuracies[i, 4] = ac3
end

gr(size=(900,400))
ctg = repeat(["SCARGC_1NN", "Static", "Sliding", "Landmark"], inner=16)
groupedbar(accuracies, group = ctg, xlabel = "Datasets",
    ylabel = "Accuracy", bar_position = :dodge, w = 0, framestyle = :box, bar_width = 0.67, xticks=(1:16, x), xrotation=30)


text:
## from module gaiaClustering
## Functions to deal with GAIA data and to normalize
##

############
## Query the GAIA data towards the coord with a radius conesearch
#####
function query_gaia(coord, radius, dump = false)
#####
    return(0)

end

function copy(s::Df)::Df
    c = Df(s.ndata, zeros(length(s.data[:,1]),s.ndata),zeros(length(s.raw[:,1]),s.ndata), zeros(length(s.err[:,1]),s.ndata),
    zeros(length(s.sourceid[:,1]),s.ndata))
    c.data[:,:] = s.data[:,:]
    c.raw[:,:]  = s.raw[:,:]
    c.err[:,:]  = s.err[:,:]
    c.sourceid[:,:]  = s.sourceid[:,:]

    return(c)
end

## dummy ...
function copy1(s::Df)::Df
    c = Df(s.ndata, zeros(length(s.data[:,1]),s.ndata),zeros(length(s.raw[:,1]),s.ndata), zeros(length(s.err[:,1]),s.ndata),
    zeros(length(s.sourceid[:,1]),s.ndata))
    c.data[:,:] = s.data[:,:]
    c.raw[:,:]  = s.raw[:,:]
    c.err[:,:]  = s.err[:,:]
    c.sourceid[:,:]  = s.sourceid[:,:]

    return(c)
end

#########
## function to create the df
######
function read_votable(voname::String)
######
    warnings = pyimport("warnings")
    warnings.filterwarnings("ignore")
    votable = pyimport("astropy.io.votable")
    vot = votable.parse(voname)
    data = vot.get_first_table()

    println("## Votable $voname read")

    ### return(data["array"]["data"])
    return(data.array.data)
end

#########
function filter_data(gaia, dist_range = [0., 2000], vra_range = [-250,250], vdec_range = [-250.,250], mag_range =[-1e9, 1e9])::Df
########
    ngaia = length(gaia)

    source_id = zeros(Int64,ngaia)
    lgal = zeros(ngaia)
    bgal = zeros(ngaia)
    ra = zeros(ngaia)
    dec= zeros(ngaia)
    distance = zeros(ngaia)
    pmra = zeros(ngaia)
    pmdec = zeros(ngaia)
    parallax = zeros(ngaia)
    vra = zeros(ngaia)
    vdec = zeros(ngaia)
    g = zeros(ngaia)
    rp = zeros(ngaia)
    bp = zeros(ngaia)
    parallax_error = zeros(ngaia)
    pmra_error     = zeros(ngaia)
    pmdec_error    = zeros(ngaia)
    radialvel      = zeros(ngaia)
    ## Galactic proper motion and velocities
    pml = zeros(ngaia)
    pmb = zeros(ngaia)
    vl = zeros(ngaia)
    vb = zeros(ngaia)

    ## Extinction A_G
    ag= zeros(ngaia)

    for i in 1:ngaia

        source_id[i]= convert(Int64, get(gaia,i-1).source_id)
        lgal[i]     = convert(Float64, get(gaia,i-1).l)
        bgal[i]     = convert(Float64, get(gaia,i-1).b)
        ra[i]       = convert(Float64, get(gaia,i-1).ra)
        dec[i]      = convert(Float64, get(gaia,i-1).dec)
        distance[i] = 1000. / convert(Float64, get(gaia,i-1).parallax)
        parallax[i] = convert(Float64, get(gaia,i-1).parallax)
        pmra[i]     = convert(Float64, get(gaia,i-1).pmra)
        pmdec[i]    = convert(Float64, get(gaia,i-1).pmdec)
        vra[i]      = 4.74e-3 * pmra[i]  * distance[i]
        vdec[i]     = 4.74e-3 * pmdec[i] * distance[i]

        ## Galactic Proper motions
        muG = PM_equatorial2galactic(pmra[i],pmdec[i]  , ra[i] , dec[i] , lgal[i])
        pml[i] = muG[1]
        pmb[i] = muG[2]
        vl[i]  = 4.74e-3 * pml[i]  * distance[i]
        vb[i]  = 4.74e-3 * pmb[i]  * distance[i]

        #fix for EDR3
        radialvel[i]    = convert(Float64, get(gaia,i-1).dr2_radial_velocity)

        ### errors.
        parallax_error[i]  = convert(Float64, get(gaia,i-1).parallax_error)
        pmra_error[i]  = convert(Float64, get(gaia,i-1).pmra_error)
        pmdec_error[i] = convert(Float64, get(gaia,i-1).pmdec_error)

        g[i]        = convert(Float64, get(gaia,i-1).phot_g_mean_mag)
        rp[i]       = convert(Float64, get(gaia,i-1).phot_rp_mean_mag)
        bp[i]       = convert(Float64, get(gaia,i-1).phot_bp_mean_mag)

        #extinction
        # fix fo EDR3, extinction not found
        ag[i]       = 99.  ##  convert(Float64, get(gaia,i-1).a_g_val)

    end

    ## Filtering ...
    i1 =  distance .> dist_range[1]
    i2 =  distance .< dist_range[2]
    i3 =  vra .> vra_range[1]
    i4 =  vra .< vra_range[2]
    i5 =  vdec .> vdec_range[1]
    i6 =  vdec .< vdec_range[2]
    i7 = g  .> mag_range[1]
    i8 = g  .< mag_range[2]
    i9  = rp  .> mag_range[1]
    i10 = rp  .< mag_range[2]
    i11  = bp  .> mag_range[1]
    i12  = bp  .< mag_range[2]

    ifinal = i1 .& i2 .& i3 .& i4 .& i5 .& i6 .& i7 .& i8 .& i9 .& i10 .& i11 .& i12

    ## G magnitude
    gbar =  g[ifinal] - 5 .* log10.(distance[ifinal]) .+ 17.

    ## Df of the filtered dat
    ndata = length(distance[ifinal])
    s = Df(ndata, zeros(8,ndata), zeros(15,ndata) , zeros(8,ndata) , zeros(1,ndata))

    s.data[1,:] = lgal[ifinal]
    s.data[2,:] = bgal[ifinal]
    s.data[3,:] = distance[ifinal]
    s.data[4,:] = vl[ifinal]
    s.data[5,:] = vb[ifinal]
    s.data[6,:] = gbar
    s.data[7,:] = g[ifinal] .- rp[ifinal]
    s.data[8,:] = bp[ifinal] .- g[ifinal]

    s.raw[1,:] = ra[ifinal]
    s.raw[2,:] = dec[ifinal]
    s.raw[3,:] = lgal[ifinal]
    s.raw[4,:] = bgal[ifinal]
    s.raw[5,:] = parallax[ifinal]
    s.raw[6,:] = pmra[ifinal]
    s.raw[7,:] = pmdec[ifinal]
    s.raw[8,:] = pml[ifinal]
    s.raw[9,:] = pmb[ifinal]
    s.raw[10,:] = g[ifinal]
    s.raw[11,:] = rp[ifinal]
    s.raw[12,:] = bp[ifinal]
    s.raw[13,:] = radialvel[ifinal]
    s.raw[14,:] = ag[ifinal]

    ## Errors ..
    s.err[1,:] = parallax_error[ifinal]
    s.err[4,:] = pmra_error[ifinal]
    s.err[5,:] = pmdec_error[ifinal]

    ## GAIA source id
    s.sourceid[1,:] = source_id[ifinal]

    println("## Filtering done ...")
    println("## Stars selected: $ndata")

    return(s)
end

function equatorial2galactic(α , δ)
    ## NGP coordinates
    αG = 192.85948
    δG = 27.12825
    lascend= 32.93

    b= asind(cosd(δ)*cosd(δG)*cosd(α-αG)+sind(δ)*sind(δG))

    x= cosd(δ)*cosd(δG)*sind(α-αG)
    y= sind(δ)-sind(b)*sind(δG)
    l= atand(y,x) + lascend
    println(x)
    println(y)

    if x>=0 && y<=0  l += 360.0 end
    if x<=0 && y<0   l += 360.0 end
    #if x<0  && y>=0  l -= 180.0 end

    return(l,b)
end

## angle between two points on a sphere
function angle4sphere(long1, lat1, long2, lat2)
    dLon= long2 - long1
    cosang= cosd(lat1)*cosd(lat2)*cosd(dLon)+sind(lat1)*sind(lat2)
    ang= acosd(cosang)
end

## Transform PM from equatorial to galactic system.
## See Poleski 1997 / arXiv
## PM,,corr is from Conrad (2015)
function PM_equatorial2galactic(μα , μδ , α , δ , l )
    ## NGP coordinates
    αG = 192.85948
    δG = 27.12825

    C1 = sind(δG)*cosd(δ) - cosd(δG)*sind(δ)*cosd(α - αG)
    C2 =  cosd(δG)*sind(α - αG)
    k = 1 / sqrt(C1^2 + C2^2)
    A = k * [C1 C2 ; -C2 C1]
    PMG = A * [μα ; μδ ]

    ## PM along gal. lat. corrected for differential velocity
    ## Oort constants. Not applied.
    #A , B = (14.5 , -13.) ./ 4.74
    # PMG[1] =  PMG[1] - (A*cosd(2l) + B)

    return(PMG)
end

## Compute the X,Y,Z galactic coordinates (centered on the Galactic Center)
## See Ellsworth-Bowers et al. (2013)
## Rgal was updated from Anderson et al. (2018)
## xg,yg,zg in pc
## Modified to be left-handed (yg is modified to -yg from original)
## see http://astro.utoronto.ca/~bovy/AST1420/notes/notebooks/A.-Coordinate-systems.html
##

function galXYZ(l,b,distance)
    Rgal= 8.34e3
    zsun= 25
    θ= asin(zsun/Rgal)

    xg= Rgal*cos(θ)-distance*(cosd(l)cosd(b)cos(θ)+sind(b)sin(θ))
    yg= distance*sind(l)cosd(b)
    zg= Rgal*sin(θ)-distance*(cosd(l)cosd(b)sin(θ)-sind(b)cos(θ))

    return(xg,yg,zg)
end

### Correction of the radial velocity
### Conrad (3025)

function RVEL_corr(rvel , distance , l)
  ## Oort's constant
  A = 14.5
  rv = rvel - A * 1e-3 * distance * sind(2l)
  return(rv)
end

### compute galactic U V W from idlastro gal_uvw.pro
##
## ra, dec: degrees
## distance; pc
## pmra,pmdec: milliarcsec/yr
## vrad: km/s

## UVW: km/s
#      U - Velocity (km/s) positive toward the Galactic *anti*center
#      V - Velocity (km/s) positive in the direction of Galactic rotation
#      W - Velocity (km/s) positive toward the North Galactic Pole
#
function galUVW(ra, dec, distance, pmra, pmdec, vrad ; LSR_vel=[-8.5 ; 13.38 ; 6.49])
    k = 4.74047     #Equivalent of 1 A.U/yr in km/s

    T=  [0.0548756   0.873437  0.483835 ;
          0.494109   -0.44483   0.746982 ;
         -0.867666   -0.198076  0.455984]

    A1= [ cosd(ra) sind(ra) 0 ; sind(ra) -cosd(ra) 0 ; 0 0 -1]
    A2= [ cosd(dec) 0 -sind(dec) ; 0 -1 0 ; -sind(dec) 0 -cosd(dec)]
    vec1 = vrad
    vec2 = k*pmra*1e-3*distance
    vec3 = k*pmdec*1e-3*distance
    v= [vec1 ; vec2 ; vec3]

    uvw= T*A1*A2*v + LSR_vel

    return(uvw)
end

######
function add_cartesian(s::Df, centering = true)::Df
######
    dfresult = copy(s)
    off = zeros(2)

    if centering
        off[1] = mean(s.data[1,:])
        off[2] = mean(s.data[2,:])
    end

    lgal = DEG2RAD .* (s.data[1,:] .- off[1])
    bgal = DEG2RAD .* (s.data[2,:] .- off[2])

    dfresult.data[1,:] = s.data[3,:] .* cos.(bgal) .* cos.(lgal)
    dfresult.data[2,:] = s.data[3,:] .* cos.(bgal) .* sin.(lgal)
    dfresult.data[3,:] = s.data[3,:] .* sin.(bgal)

    println("## Cartesian transformation done ...")

    return(dfresult)
end

######
function  normalization_PerBlock(s::Df, block , weightblock, norm , density = false , verbose = true)
######
    dfresult = copy(s)
    ndf = size(s.data)
    scale8d = zeros(ndf[1])
    vector8d = 0.

    ind = 1
    for aw in zip(block,weightblock)
        weight = aw[2]
        for ak in aw[1]
            normK = normalizationVector(norm, density, dfresult.data[ak,:])
            # normK[2] = normK[2] * totalWeight
            dfresult.data[ak,:]    =   weight .* (s.data[ak,:] .- normK[1] ) ./ normK[2]
            scale8d[ind] = weight / normK[2]
            vector8d += scale8d[ind] ^ 2
            ind += 1
        end
    end

    vector8d = sqrt(vector8d)
    scale8d[:] = scale8d[:] ./ vector8d
    dfresult.data[:,:] = dfresult.data[:,:] ./ vector8d

    if verbose
        println("## Normalization $norm done...")
        println("### [1pc,1pc,1pc,1km/s,1km/s,1mag,1mag,1mag] equivalent to $scale8d")
        println("##")
    end

    return(dfresult , scale8d)
end

######
function normalizationVector(norm, density, arr)
######
    vecNorm = [0.0, 1.0]

    if norm == "identity"
        vecNorm = [0.0, 1.0]

        elseif norm == "normal"
            stdArr = std(arr)
            meanArr = mean(arr)
            vecNorm  = [meanArr , stdArr]

        elseif norm == "minmax"
            minarr  = minimum(vcat(arr...))
            maxarr  = maximum(vcat(arr...))
            vecNorm = [minarr, maxarr-minarr]

    end

    if density
        vecNorm[2] = vecNorm[2] * length(arr)
    end

    return(vecNorm)
end

######
function subsetDf(df::Df, indx)::Df
######

  ndat = length(indx)
  subset = Df(ndat, df.data[:,indx], df.raw[:,indx], df.err[:,indx], df.sourceid[:,indx])

    return(subset)
end

####
# Create the DataFrame to save the cluster...
##
function export_df(votname, ocdir, df , dfcart, labels , labelmax)
    ra= df.raw[1, labels[labelmax]]
    dec= df.raw[2,labels[labelmax]]
    l= df.data[1, labels[labelmax]]
    b= df.data[2,labels[labelmax]]
    d= df.data[3,labels[labelmax]]
    pmra= df.raw[6, labels[labelmax]]
    pmdec= df.raw[7, labels[labelmax]]
    X= dfcart.data[1, labels[labelmax]]
    Y= dfcart.data[2, labels[labelmax]]
    Z= dfcart.data[3, labels[labelmax]]
    vl= df.data[4,labels[labelmax]]
    vb= df.data[5,labels[labelmax]]
    vrad= df.raw[13,labels[labelmax]]
    gbar= df.raw[10,labels[labelmax]]
    rp= df.raw[11,labels[labelmax]]
    bp= df.raw[12,labels[labelmax]]
    ag= df.raw[14,labels[labelmax]]

    source_id= df.sourceid[1,labels[labelmax]]

    # galactocentric coordinates...
    s=size(ra)
    xg= zeros(s[1])
    yg= zeros(s[1])
    zg= zeros(s[1])

    ind= 1
    for i in 1:s[1]
        xx, yy, zz= galXYZ.(l[i],b[i],d[i])
        xg[i]= xx
        yg[i]= yy
        zg[i]= zz
    end

    oc= DataFrame(sourceid=source_id,ra=ra,dec=dec,l=l,b=b, distance=d,pmra=pmra, pmdec=pmdec, X=X,Y=Y,Z=Z,vl=vl,
        vb=vb,vrad=vrad, Xg=xg,Yg=yg,Zg=zg,gbar=gbar,rp=rp,bp=bp, ag=ag)

    name= split(votname,".")
    infix= ""
    for iname in name
        if iname != "vot"
            infix *= iname*"."
        end
    end
    infix *= "oc.csv"
    filename= @sprintf("%s/%s",ocdir, infix)
    CSV.write(filename,oc,delim=';')
    @printf("### %s created \n",filename)
end
#######################################
## a built-in version of getdata
function get_data(m::GaiaClustering.meta)
    println("## Distance cut : $(m.mindist) $(m.maxdist) pc")

    data       = read_votable(m.votname)
    df         = filter_data(data, [m.mindist, m.maxdist])
    dfcart     = add_cartesian(df)
    blck       = [[1,2,3],[4,5], [6,7,8]]
    wghtblck   = [4.0,5.0,1.0]
    norm       = "identity"

    dfcartnorm , scale8 = normalization_PerBlock(dfcart, blck, wghtblck , norm, false)

    return(df, dfcart , dfcartnorm)
end

max_stars_repo_name: bosscha/gaia-shock

max_issues_repo_name: bosscha/gaia-shock

max_forks_repo_name: bosscha/gaia-shock


text:
"""
$(TYPEDEF)

A system of operators and caches for immersed layer problems. This is constructed
by [`__init`](@ref)
"""
struct ILMSystem{PT,BCT<:BasicILMCache,ECT<:Union{AbstractExtraILMCache,Nothing}}

  base_cache :: BCT
  extra_cache :: ECT

end

"""
    __init(prob::AbstractILMProblem)

Initialize `ILMSystem` with the given problem `prob` specification.
Depending on the type of problem, this sets up a base cache of scalar or
vector type, as well as an optional extra cache
"""
function __init(prob::AbstractILMProblem{DT,ST}) where {DT,ST}
    @unpack g, bodies, phys_params, bc, f_funcs = prob

    if typeof(prob) <: AbstractScalarILMProblem
        base_cache = SurfaceScalarCache(bodies,g,ddftype=DT,scaling=ST,phys_params=phys_params,bc=bc,f_funcs=f_funcs)
    elseif typeof(prob) <: AbstractVectorILMProblem
        base_cache = SurfaceVectorCache(bodies,g,ddftype=DT,scaling=ST,phys_params=phys_params,bc=bc,f_funcs=f_funcs)
    end

    extra_cache = prob_cache(prob,base_cache)

    return ILMSystem{typeof(prob),typeof(base_cache),typeof(extra_cache)}(base_cache,extra_cache)

end

# Create the basic solve function, to be extended
function solve(prob::AbstractILMProblem,sys::ILMSystem) end

## Extend functions on `BasicILMCache` type to `ILMSystem`
for f in [:zeros_surface,:zeros_grid,:zeros_gridcurl,:zeros_gridgrad,
          :similar_surface,:similar_grid,:similar_gridcurl,:similar_gridgrad,
          :ones_surface,:ones_grid,:ones_gridgrad,:ones_gridcurl,
          :x_grid,:y_grid,:x_gridcurl,:y_gridcurl,
          :normals,:areas,:points,
          :create_nRTRn,:create_GLinvD,:create_CLinvCT,:create_CL2invCT,
          :create_RTLinvR,
          :create_GLinvD_cross,:create_surface_filter]
   @eval $f(sys::ILMSystem) = $f(sys.base_cache)
end

for f in [:regularize!, :interpolate!, :regularize_normal!,
          :normal_interpolate!,:regularize_normal_cross!,
          :normal_cross_interpolate!,
          :surface_curl!,:surface_divergence!,:surface_grad!,
          :surface_curl_cross!,:surface_divergence_cross!,:surface_grad_cross!]
  @eval $f(a,b,sys::ILMSystem) = $f(a,b,sys.base_cache)
end

for f in [:norm,:integrate]
   @eval $f(a,sys::ILMSystem) = $f(a,sys.base_cache)
   @eval $f(a,sys::ILMSystem,i::Int) = $f(a,sys.base_cache,i)
end

for f in [:view]
  @eval $f(a::PointData,sys::ILMSystem,i::Int) = $f(a,sys.base_cache,i)
end

for f in [:dot,:copyto!]
   @eval $f(a,b,sys::ILMSystem) = $f(a,b,sys.base_cache)
   @eval $f(a,b,sys::ILMSystem,i::Int) = $f(a,b,sys.base_cache,i)
end

for f in [:RegularizationMatrix,:InterpolationMatrix]
  @eval CartesianGrids.$f(sys::ILMSystem,args...) = $f(sys.base_cache,args...)
end


text:
using CartesianCutCell
using Test

include("rectangular.jl")
#include("cylindrical.jl")


text:
@testset "Kronecker products" begin

    A = randn(4, 4)
    B = Array{Float64, 2}([1 2 3;
         4 5 6;
         7 -2 9])
    C = rand(5, 6)
    D = rand(4, 4)

    v = rand(12)

    K = A ⊗ B
    K3 = kronecker(A, B, C)

    X = kron(A, B)  # true result

    @testset "Types and basic properties" begin
        @test issquare(A)
        @test !issquare(C)

        @test getmatrices(A)[1] === A

        @test !issymmetric(A)

        @test issquare(K)
        @test !issymmetric(K)

        @test K ≈ X

        @test order(A) == 1
        @test order(K) == 2

        @test eltype(K) <: Float64

        @test K isa AbstractKroneckerProduct
        @test K isa AbstractKroneckerProduct{Float64}

        @test K isa AbstractMatrix{Float64}
    end

    @testset "Linear algebra" begin
        @test tr(K) ≈ tr(X)
        @test det(K) ≈ det(X)
        @test !isposdef(K)
        @test transpose(K) ≈ transpose(X)
        @test conj(K) ≈ conj(X)
        @test K' ≈ X'
        @test inv(K) ≈ inv(X)

        # test on pos def functions
        As = A' * A
        Bs = B * B'
        @test logdet(As ⊗ Bs) ≈ log(det(As ⊗ Bs)) ≈ log(det(kron(As, Bs)))
    end

    @testset "Mismatch errors" begin
        P, Q = rand(10, 4), rand(4, 5)
        Kns = P ⊗ Q
        @test_throws DimensionMismatch inv(Kns)
        @test_throws DimensionMismatch det(Kns)
        @test_throws DimensionMismatch Kns * [1, 2, 3]
    end

    @testset "Higher order" begin
        @test order(K ⊗ A) == 3
        @test order(K3) == 3
        @test collect(K3) ≈ kron(X, C)

    end

    @testset "Kronecker powers" begin
        Kpow = ⊗(A, 5)
        @test order(Kpow) == 5
        @test size(Kpow, 1) == 4^5
        @test Kpow[1,1] ≈ A[1,1]^5
    end

    @testset "kron" begin
        @test kron(A ⊗ B, C) ≈ kron(A, B, C)
        @test kron(A, B ⊗ C) ≈ kron(A, B, C)
        @test kron(A ⊗ B, C ⊗ D) ≈ kron(A, B, C, D)
        @test collect((A⊗B) ⊗ (C⊗D)) ≈ kron(A, B, C, D)
    end

    @testset "Mixed product" begin
        A = rand(5, 4)
        B = rand(2, 3)
        C = rand(4, 6)
        D = rand(3, 4)

        K1 = (A ⊗ B)
        K2 = (C ⊗ D)

        @test K1 * K2 ≈ collect(K1) * collect(K2)

        @test_throws DimensionMismatch (A ⊗ C) * (B ⊗ D)
        @test_throws DimensionMismatch (A ⊗ D) * (C ⊗ B)
    end

    @testset "Add to dense" begin
        @test K + X ≈ Matrix(K) + X
        @test X + K ≈ X + Matrix(K)
    end

    @testset "Scalar multiplication" begin
        @test 3.0K ≈ 3.0X
        @test K * 2 ≈ 2X
        @test π * K3 ≈ π * collect(K3)
        @test 3.0K isa AbstractKroneckerProduct
        @test K * 2 isa AbstractKroneckerProduct
        @test 2(K ⊗ K) isa AbstractKroneckerProduct
    end

    @testset "Inplace scalar multiplication" begin
        A = rand(2, 2)
        B = rand(3, 4)

        K = copy(A) ⊗ copy(B)
        lmul!(3, K)
        rmul!(K, 2)
        @test K.A ≈ 3A
        @test K.B ≈ 2B
        @test K ≈ 6kron(A, B)
    end

    @testset "Solving Linear Systems" begin
        n = 8
        m = 16
        # 1) testing square kronecker product with square A, B
        A = randn(n, n)
        B = randn(m, m)
        K = A ⊗ B
        x = randn(n * m)
        b = K*x
        @test K\b ≈ x
        @test (x'*K)/K ≈ x'

        # testing least squares solution
        function test_ls_solve(size_A, size_B)
            A = randn(size_A)
            B = randn(size_B)
            x = randn(size(A, 2) * size(B, 2))
            K = kronecker(A, B)
            b = K*x
            b .+= randn(size(b)) # this moves b out of range(K), necessitating least-squares
            xls = K\b
            return K'*(K*xls) ≈ K'b # test via normal equations
        end

        # 2) kronecker product is square, but A, B aren't
        size_A = (m, n)
        size_B = (n, m)
        @test test_ls_solve(size_A, size_B)

        # 3) kronecker product is not square, dimension of x is larger than b
        size_A = (n, n)
        size_B = (n, m)
        @test test_ls_solve(size_A, size_B)

        # 4) kronecker product is not square, dimension of x is smaller than b
        size_A = (m, n)
        size_B = (m, m)
        @test test_ls_solve(size_A, size_B)
    end
end


text:
# Include packages
using NetCDF
using Printf: @printf, @sprintf
using Dates

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

# Name of NetCDF file
ncfile="./data/wspd.mon.mean.nc"
if !isfile(ncfile)
    # download NetCDF file
    # These data were collected and distributed by
    # Research Institute for Sustainable Humanosphere,
    # Kyoto University (http://database.rish.kyoto-u.ac.jp/index-e.html).
    url="http://database.rish.kyoto-u.ac.jp/arch/ncep/data/ncep.reanalysis.derived/surface/wspd.mon.mean.nc"
    download(url, ncfile)
end

# Get information
# ncinfo(ncfile)

# Get variables
lon = convert.(Float64, ncread(ncfile,"lon"))
lat = convert.(Float64, ncread(ncfile,"lat"))
wspd = permutedims(ncread(ncfile,"wspd"), [2 1 3])
# convert the order of latitudes, 90:-90 to -90:90
lat = reverse(lat, dims=1)
wspd = reverse(wspd, dims=1)
torg = ncread(ncfile,"time")
nt = length(torg);
T = DateTime(1800,1,1)+Hour.(Int.(torg))

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

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

using Plots
pyplot()

####################
## functions
####################
function DrawSnapShot(k::Int, lon, lat, wspd, T::Array{DateTime,1})
    # setup for ticks
    xt = collect(Int64, 0:60:360)
    xtl = [collect(0:60:180);collect(-120:60:0)]
    xtl = [@sprintf("%d",xtl[i]) for i=1:length(xtl)]
    # clibrary
    clibrary(:misc)
    # filled contour
    ## (2018/03 メモ)pyplotバックエンドでは，clims=()の設定は反映されない。なぜ？
    contour(lon, lat, wspd[:,:,k],fill=true, clims=(0.,16.), color=(:rainbow),
             size=(800, 400), tickfont=12, axis_ratio=:equal,
             xlabel="Longitude", ylabel="Latitude", guidefont=12,
             colorbar_title="(m/s)",
             title=Dates.format(T[k], "yyyy/mm"), titlefont=14,
             xticks=(xt, xtl),
             )
end
####################

# DrawSnapShot(1, lon, lat, wspd, T) # 1st step
#if !isdir("./forgif"); mkdir("./forgif"); end

# animation
anim = @animate for k=1:24
    @printf("%d, ",k)
    DrawSnapShot(k, lon, lat, wspd, T)
end
gifname = joinpath(figdir,"monthlywspd.gif")
if isfile(gifname); rm(gifname); end
gif(anim, gifname, fps=4) #save the animation
####################


text:
module aoc_06
    using Test
    export read_input, part1, part2
    export test_part1, test_part2, test
    export input, demo

    const input_r = r""

    function read_input(filename = "input.txt")
        lines = readlines(filename)
        # for line in eachline(filename)
        groups = Vector{Vector{String}}()
        group = Vector{String}()
        for line in lines
            if isempty(line)
                push!(groups, copy(group))
                empty!(group)
            else
                push!(group, line)
            end
        end
        push!(groups, copy(group))
        # [parse_line(line) for line in eachline(filename)]
        groups
    end

    function parse_line(line)
    end

    const input = read_input("input.txt")
    const demo = read_input("demo.txt")

    function part1(input = input)
        sum(length.(unique.(Vector{Char}.(join.(input)))))
    end

    function test_part1()
        @test part1( demo ) == 11
        @test part1() == 6457
    end

    function part2(input = input)
        sum(length.([intersect.( group... ) for group in input]))
    end

    function test_part2()
        @test part2( demo ) == 6
        @test part2() == 3260
    end

    function test()
        test_part1()
        test_part2()
    end
end


text:
# module Bukdu.Plug

# follow the https://github.com/elixir-plug/plug/blob/master/lib/plug/conn.ex

"""
    Conn
"""
mutable struct Conn <: AbstractPlug
    request::Deps.Request

    # Fetchable fields
    # cookies           # the request cookies with the response cookies  FIXME
    body_params::Assoc  # Plug.Parsers
    query_params::Assoc # fetch_query_params
    path_params::Assoc
    params::Assoc       # merge(body_params, query_params, path_params)

    # Connection fields
    halted::Bool        # the boolean status on whether the pipeline was halted
end
Conn(request::Deps.Request) = Conn(request, Assoc(), Assoc(), Assoc(), Assoc(), false)

# module Bukdu.Plug

max_stars_repo_name: wookay/TestBukdu

max_issues_repo_name: wookay/TestBukdu

max_forks_repo_name: wookay/TestBukdu


text:
module IndexFunArrays

using LinearAlgebra
using Base.Iterators

export IndexFunArray
export selectsizes, get_scale

const DEFAULT_T = Float64
const IterType = Union{NTuple{N,Tuple} where N, Vector, Matrix, Base.Iterators.Repeated}

include("scaling_offset_types.jl")
include("utils.jl")
include("other_ifas.jl")
# include docstrings
include("docstrings.jl")
include("scalar_ifas.jl")
include("window_ifas.jl")
include("tuple_ifas.jl")
include("ifa_type.jl")
# include("derivative_support.jl")

# These are the type promotion rules, taken from float.jl but written in terms of types
# see also
promote_type()
default_type(::Type{Bool}, def_T)    = def_T
default_type(::Type{Int8}, def_T)    = def_T
default_type(::Type{Int16}, def_T)   = def_T
default_type(::Type{Int32}, def_T)   = def_T
default_type(::Type{Int64}, def_T)   = def_T # LOSSY
default_type(::Type{Int128}, def_T)  = def_T # LOSSY
default_type(::Type{UInt8}, def_T)   = def_T
default_type(::Type{UInt16}, def_T)  = def_T
default_type(::Type{UInt32}, def_T)  = def_T
default_type(::Type{UInt64}, def_T)  = def_T # LOSSY
default_type(::Type{UInt128}, def_T) = def_T # LOSSY
default_type(::Type{T}, def_T) where{T} = T # all other types remain to be the same

end # module


text:
push!(LOAD_PATH, "./")
using Rando
using Stiffness
function prSubDiag(m)

    for i = 0:5
        for j=0:i
            print(m[div(i*(i+1),2)+j+1]," ")
        end
        println()
    end
end
function RandomTriangle!(R::RandoData,x::Array{Float64,1},y::Array{Float64,1})
    for i in 1:3
        x[i]=Rando.fv!(R,10.)
    end
    for i in 1:3
        y[i]=Rando.fv!(R,10.)
    end
end
S=StiffnessData()
x=Float64[0.,1.,0.]
y=Float64[0.,0.,1.]
m=zeros(21)
print("\nVerify that, on the reference element, we are coherent with sage ")
println("(see ../sage/):\n")
Stiffness.op!(S,x,y,m)
prSubDiag(m)
print("\nWe must get the same result if we dilate the triangle:\n")
x*=2.0
y*=2.0
Stiffness.op!(S,x,y,m)
prSubDiag(m)

R=RandoData()

# be sure to run once before actually running the benchmark !
RandomTriangle!(R,x,y)
Stiffness.op!(S,x,y,m)
#

println("\nNow, start the benchmark:")
const ntri=1000000
println(ntri," triangles.")

t1 = time_ns()
for t=1:ntri
    RandomTriangle!(R,x,y)
    Stiffness.op!(S,x,y,m)
end
t1 = time_ns()-t1
println("first phase: ",float(t1)*10.0^(-9)," seconds.")

t11 = time_ns()
for t=1:ntri
    RandomTriangle!(R,x,y)
end
t11 = time_ns()-t11
println("second phase: ",float(t11)*10.0^(-9)," seconds.")
T=t1-t11
Tsecond=float(T)*10.0^(-9)
println("Total time: ",Tsecond)
byTr=Tsecond/ntri
println("Time by triangle: ",byTr," seconds.")
open("RunningOn"*gethostname(),"w") do f
    write(f,string(byTr),"\n")
end
println("end.")

# to profile the code,uncomment:
# ntri1=100000
# @profile for t=1:ntri1
#     RandomTriangle!(R,x,y)
#     Stiffness.op!(S,x,y,m)
#  end
#  Profile.print(format=:flat)

max_stars_repo_name: fluiddyn/BenchmarksPythonJuliaAndCo

max_issues_repo_name: fluiddyn/BenchmarksPythonJuliaAndCo

max_forks_repo_name: fluiddyn/BenchmarksPythonJuliaAndCo


text:
"""
    CNumber <: Number

Abstract type for all symbolic numbers, i.e. [`Parameter`](@ref), [`average`](@ref).
"""
abstract type CNumber <: Number end

"""
    Parameter <: CNumber

Type used as symbolic type in a `SymbolicUtils.Sym` variable to represent
a parameter.
"""
struct Parameter <: CNumber
    function Parameter(name; metadata=source_metadata(:Parameter, name))
        s = SymbolicUtils.Sym{Parameter, typeof(metadata)}(name, metadata)
        return SymbolicUtils.setmetadata(s, MTK.MTKParameterCtx, true)
    end
end

# Promoting to CNumber ensures we own the symtype; could be used to dispatch
# on Base methods (e.g. latex printing)
Base.promote_rule(::Type{<:CNumber},::Type{<:Number}) = CNumber

Base.one(::Type{Parameter}) = 1
Base.zero(::Type{Parameter}) = 0

"""
    @cnumbers(ps...)

Convenience macro to quickly define symbolic cnumbers.

Examples
========
```
julia> @cnumbers ω κ
(ω, κ)
```
"""
macro cnumbers(ps...)
    ex = Expr(:block)
    pnames = []
    for p in ps
        @assert p isa Symbol
        push!(pnames, p)
        d = source_metadata(:cnumbers, p)
        ex_ = Expr(:(=), esc(p), Expr(:call, :Parameter, Expr(:quote, p), Expr(:kw, :metadata, Expr(:quote, d))))
        push!(ex.args, ex_)
    end
    push!(ex.args, Expr(:tuple, map(esc, pnames)...))
    return ex
end

"""
    cnumbers(symbols::Symbol...)
    cnumbers(s::String)

Create symbolic cnumbers.

Expamples
=========
```
julia> ps = cnumbers(:a, :b)
(a, b)

julia> cnumbers("a b") == ps
true
```
"""
function cnumbers(syms::Symbol...)
    ps = Tuple(Parameter(s; metadata=source_metadata(:cnumbers, s)) for s in syms)
    return ps
end
function cnumbers(s::String)
    syms = [Symbol(p) for p in split(s, " ")]
    return cnumbers(syms...)
end


text:
# Bi-objective unidimensionnal 01 knapsack problem (biukp)
#
# Example 2 with 20 variables

# ---- Packages to use
using vOptGeneric, JuMP, GLPK, LinearAlgebra

# ---- Values of the instance to solve
p1 = [77,94,71,63,96,82,85,75,72,91,99,63,84,87,79,94,90,60,69,62]
p2 = [65,90,90,77,95,84,70,94,66,92,74,97,60,60,65,97,93,60,69,74]
w  = [80,87,68,72,66,77,99,85,70,93,98,72,100,89,67,86,91,79,71,99]
c  = 900
size = length(p1)

# ---- setting the model
m = vModel( GLPK.Optimizer ) ; JuMP.set_silent(m)
@variable( m, x[1:size], Bin )
@addobjective( m, Max, dot(x, p1) )
@addobjective( m, Max, dot(x, p2) )
@constraint( m, dot(x, w) <= c )

# ---- Invoking the solver (epsilon constraint method)
vSolve( m, method=:epsilon, step=0.5, verbose=true )

# ---- Querying the results
Y_N = getY_N(m)

# ---- Displaying the results (X_E and Y_N)
for n = 1:length(Y_N)
    X = value.(x, n)
    print(findall(elt -> elt ≈ 1, X))
    println("| z = ",Y_N[n])
end

# using PyPlot
# f1 = map(y -> y[1], Y_N)
# f2 = map(y -> y[2], Y_N)
# xlabel("z1")
# ylabel("z2")
# plot(f1,f2,"bx", markersize = "6")
# !isinteractive() && show()

max_stars_repo_name: vOptSolver/vOptGeneric

max_issues_repo_name: vOptSolver/vOptGeneric

max_forks_repo_name: vOptSolver/vOptGeneric