text:
using Surrogates
using Test
using SafeTestsets
using Pkg

function dev_subpkg(subpkg)
	subpkg_path = joinpath(dirname(@__DIR__), "lib", subpkg)
	Pkg.develop(PackageSpec(path=subpkg_path))
end
for pkg in ["SurrogatesAbstractGPs", "SurrogatesFlux", "SurrogatesPolyChaos",
				"SurrogatesRandomForest", "SurrogatesSVM"]
	@time begin
		dev_subpkg(pkg)
		Pkg.test(pkg)
	end
end

@time @safetestset  "Radials.jl" begin include("Radials.jl") end
@time @safetestset  "Kriging.jl" begin include("Kriging.jl") end
@time @safetestset  "Sampling" begin include("sampling.jl") end
@time @safetestset  "Optimization" begin include("optimization.jl") end
@time @safetestset  "LinearSurrogate" begin include("linearSurrogate.jl") end
@time @safetestset  "Lobachevsky" begin include("lobachevsky.jl") end
@time @safetestset  "InverseDistanceSurrogate" begin include("inverseDistanceSurrogate.jl") end
@time @safetestset  "SecondOrderPolynomialSurrogate" begin include("secondOrderPolynomialSurrogate.jl") end
# @time @safetestset  "AD_Compatibility" begin include("AD_compatibility.jl") end
@time @safetestset  "Wendland" begin include("Wendland.jl") end
@time @safetestset  "VariableFidelity" begin include("VariableFidelity.jl") end
@time @safetestset  "Earth" begin include("earth.jl") end
@time @safetestset  "Gradient Enhanced Kriging" begin include("GEK.jl") end
@time @safetestset  "Section Samplers" begin include("SectionSampleTests.jl") end


text:
function act_1_solution()
  cur_T = 20.6

  cur_reactor = act_1_deck(cur_T)

  cur_reactor.n_bar = 1.3
  cur_reactor.I_P = 10.95

  cur_reactor.R_0 = 6.25
  cur_reactor.B_0 = 6.0

  cur_reactor.P_F = 1813

  cur_reactor.f_BS = 0.91
  cur_reactor.f_CD = 1.00 - cur_reactor.f_BS
  cur_reactor.f_IN = 0.00

  cur_reactor.beta_N = cur_reactor.max_beta_N
  cur_reactor.q_95 = 4.5
  cur_reactor.P_W = 2.45

  cur_reactor.volume = 582

  cur_reactor
end

ext: jl


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed 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.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

module ConfigMCTS

export defineMCTSParams

using AdaptiveStressTesting

function defineMCTSParams(;
                          d::Int64 = 51,
                          ec::Float64 = 100.0,
                          n::Int64 = 1000,
                          k::Float64 = 0.5,
                          alpha::Float64 = 0.85,
                          kp::Float64 = 1.0,
                          alphap::Float64 = 0.0,
                          clear_nodes::Bool = true,
                          maxtime_s::Float64 = realmax(Float64),
                          rng_seed::UInt64 = UInt64(0)
                          )

  return p = DPWParams(d, ec, n, k, alpha, kp, alphap, clear_nodes, maxtime_s, rng_seed)
end

end #module

max_stars_repo_name: sisl/RLESCAS

max_issues_repo_name: sisl/RLESCAS

max_forks_repo_name: sisl/RLESCAS


text:
using Multisets
using Test

A = Multiset(1:5)
@test Set(A) == A
B = A + A
@test length(B) == 2 * length(A)
@test A ⊆ B
@test B ⊇ A
B[1] = 0
@test Set(B) == Set(2:5)
@test length(A ∪ B) == 9
@test length(A ∩ B) == 4
@test length(A * B) == 40
@test length(A + B) == 13
@test length(A - B) == 1
@test length(B - A) == 4
@test collect(A) == collect(1:5)
@test hash(A) != hash(B)

@test Multiset(1, 2, 3, 1) == Multiset([1, 1, 2, 3])
@test eltype(Multiset(1, 2, 3)) == Int
@test Multiset(1.0, 2.0) == Multiset(1 + 0im, big(2))

A = Multiset(1, 1, 1, 2, 1, 3, 3)
@test sum(A) == 12
@test sort(collect(keys(A))) == [1, 2, 3]
@test sort(collect(values(A))) == [1, 2, 4]
@test sort(collect(pairs(A))) == [1 => 4, 2 => 1, 3 => 2]

B = Set(3:5)
@test A - B == Multiset(1, 1, 1, 1, 2, 3)
@test B - A == BitSet(4:5)

@test length(A + B) == length(A) + length(B) == length(B + A)
@test A ∩ B == Set(3) == B ∩ A
@test A ∪ B == Multiset(1, 1, 1, 1, 2, 3, 3, 4, 5) == B ∪ A

X = Set(A)

@test X ⊆ A
@test A ⊇ X
@test 1 ∈ A
@test 9 ∉ A

A = Multiset(1, 2, 3)
B = Multiset(1, 2, 3, 4)
B[4] = 0
@test A == B
@test hash(A) == hash(B)

set_key_value_show()
A = Multiset(1, 1, 1, 2, 1, 3, 3)

@test string(A) == "Multiset{Int64}(2 => 1, 3 => 2, 1 => 4)"

B = Multiset()

@test string(B) == "Multiset{Any}()"

A = Multiset(1:5)
B = copy(A)
@test A == B
empty!(A)
@test length(A) == 0
@test A != B

nothing


text:
@static if VERSION < v"1.7"
    function Base.pkgdir(m::Module, path_::String, paths...)
        joinpath(pkgdir(m), path_, paths...)
    end
end


text:
using CSV
using DelimitedFiles
using DSP
using Random
using LinearAlgebra
using DataFrames

function loadtestdata(
    testCase::String;
    dataPath::String = (@__DIR__) * "/data_new_testcases",
)
    #println(pwd()) # to debug github action
    data = readdlm(joinpath(dataPath, "$(testCase)_data.csv"), ',', Float64, '\n')
    data = dropdims(data, dims = 1) # convert to vector
    evts = CSV.read(joinpath(dataPath, "$(testCase)_events.csv"), DataFrame)
    return data, evts
end

function gridexpand(conditionA = [0, 1.0], continuousA = [-1.0, 0, 1.0])

    tmp = reshape(
        [[x, y] for x in conditionA, y in continuousA],
        length(conditionA) * length(continuousA),
    )
    evts_grid = DataFrame(hcat(tmp...)')
    rename!(evts_grid, ["conditionA", "continuousA"])
    return evts_grid
end

simulate_lmm(args...; kwargs...) = simulate_lmm(Random.GLOBAL_RNG, args...; kwargs...)

function simulate_lmm(
    rng::AbstractRNG,
    τ = 1.5,
    fs = 12;
    β = [0.0, -1.0],
    σs = [1, 1, 1],
    σ = 1,
    n_sub = 20,
    n_item = 30,
    noise_type = "AR-exponential",
)
    rng_copy = deepcopy(rng)

    subj_btwn = item_btwn = both_win = nothing
    #subj_btwn = Dict("age" => ["O", "Y"])

    # there are no between-item factors in this design so you can omit it or set it to nothing
    item_btwn = Dict("stimType" => ["I", "II"])

    # put within-subject/item factors in a Dict
    #both_win = Dict("condition" => ["A", "B"])

    # simulate data
    evt = DataFrame(
        simdat_crossed(
            n_sub,
            n_item,
            subj_btwn = subj_btwn,
            item_btwn = item_btwn,
            both_win = both_win,
        ),
    )

    #    f1 = @formula dv ~ 1 + age * condition  + (1+condition|item) + (1+condition|subj);
    f1 = @formula dv ~ 1 + stimType + (1 + stimType | subj) + (1 | item)
    m = MixedModels.fit(MixedModel, f1, evt)

    # set the random effects

    #gen_han(τ,fs,1)
    basis = gen_han(τ, fs, 2)

    epoch_dat = zeros(Int(τ * fs), size(evt, 1))

    #MixedModels doesnt really support σ==0 because all Ranef are scaled residual variance
    σ_lmm = 0.0001
    σs = σs ./ σ_lmm

    for t = 1:size(epoch_dat, 1)
        b = basis[t]

        MixedModelsSim.update!(m, create_re(b .* σs[1], b .* σs[2]), create_re(b .* σs[3]))
        simulate!(deepcopy(rng_copy), m, β = [b .* β[1], b .* β[2]], σ = σ_lmm)

        epoch_dat[t, :] = m.y
    end

    epoch_dat = reshape(epoch_dat, (1, size(epoch_dat)...))

    # add some noise
    if noise_type == "normal"
        epoch_dat = epoch_dat .+ randn(rng, size(epoch_dat)) .* σ
    elseif noise_type == "AR-exponential"
        epoch_dat = epoch_dat .+ gen_noise_exp(rng, size(epoch_dat)) .* σ

    end

    return evt, epoch_dat
end

function gen_han(τ, fs, peak)
    hanLen = Int(τ * fs / 3)
    han = hanning(hanLen, zerophase = false)
    sig = zeros(Int(τ * fs))
    sig[1+hanLen*(peak-1):hanLen*peak] .= han
    return sig
end

function circulant(x)
    # Author: Jaromil Frossard
    # returns a symmetric matrix where X was circ-shifted.
    lx = length(x)
    ids = [1:1:(lx-1);]
    a = Array{Float64,2}(undef, lx, lx)
    for i = 1:length(x)
        if i == 1
            a[i, :] = x
        else
            a[i, :] = vcat(x[i], a[i-1, ids])
        end
    end
    return Symmetric(a)
end

function exponentialCorrelation(x; nu = 1, length_ratio = 1)
    # Author: Jaromil Frossard
    # generate exponential function
    R = length(x) * length_ratio
    return exp.(-3 * (x / R) .^ nu)
end

function noise_exp(rng, n_t)
    Σ = circulant(exponentialCorrelation([0:1:(n_t-1);], nu = 1.5))
    Ut = LinearAlgebra.cholesky(Σ).U'
    return (randn(rng, n_t)'*Ut')[1, :]
end

function gen_noise_exp(rng, si)
    n = hcat(map(x -> noise_exp(rng, si[2]), 1:si[3])...)
    return reshape(n, si)
end

max_stars_repo_name: unfoldtoolbox/jlunfold

max_issues_repo_name: unfoldtoolbox/jlunfold

max_forks_repo_name: unfoldtoolbox/jlunfold


text:
using QuantumPropagators: propstep!, write_to_storage!, get_from_storage!
using QuantumControlBase: evalcontrols!
using QuantumControlBase: resetgradvec!
using QuantumControlBase.ConditionalThreads: @threadsif
using LinearAlgebra
using Printf

"""Optimize a control problem using GRAPE.

```julia
result = optimize_grape(problem)
```

optimizes the given
control [`problem`](@ref QuantumControlBase.ControlProblem),
returning a [`GrapeResult`](@ref).

!!! note

    It is recommended to call [`optimize`](@ref QuantumControlBase.optimize)
    with `method=:GRAPE` instead of calling `optimize_grape` directly.

Keyword arguments that control the optimization are taken from the keyword
arguments used in the instantiation of `problem`.

# Required problem keyword arguments

* `J_T`: A function `J_T(ϕ, objectives, τ=τ)` that evaluates the final time
  functional from a list `ϕ` of forward-propagated states and
  `problem.objectives`.
* `gradient`:  A function `gradient!(G, τ, ∇τ)` that stores the gradient of
  `J_T` in `G`.

# Optional problem keyword arguments

* `update_hook`
* `info_hook`
* `check_convergence`
* `x_tol`
* `f_tol`
* `g_tol`
* `show_trace`
* `extended_trace`
* `show_every`
* `allow_f_increases`
* `optimizer`
* `prop_method`/`fw_prop_method`/`bw_prop_method`: The propagation method to
  use for each objective, see below.
* `prop_method`/`fw_prop_method`/`grad_prop_method`: The propagation method to
  use for the extended gradient vector for each objective, see below.

The propagation method for the forward propagation of each objective is
determined by the first available item of the following:

* a `fw_prop_method` keyword argument
* a `prop_method` keyword argument
* a property `fw_prop_method` of the objective
* a property `prop_method` of the objective
* the value `:auto`

The propagation method for the backword propagation is determined similarly,
but with `bw_prop_method` instead of `fw_prop_method`. The propagation method
for the forward propagation of the extended gradient vector for each objective
is determined from `grad_prop_method`, `fw_prop_method`, `prop_method` in order
of precedence.
"""
function optimize_grape(problem)
    update_hook! = get(problem.kwargs, :update_hook, (args...) -> nothing)
    info_hook = get(problem.kwargs, :info_hook, print_table)
    check_convergence! = get(problem.kwargs, :check_convergence, res -> res)

    wrk = GrapeWrk(problem)

    χ = wrk.bw_states
    Ψ = wrk.fw_states
    Ψ̃ = wrk.fw_grad_states
    τ = wrk.result.tau_vals
    ∇τ = wrk.tau_grads
    N_T = length(wrk.result.tlist) - 1
    N = length(wrk.objectives)
    L = length(wrk.controls)
    X = wrk.bw_storage

    gradfunc! = wrk.kwargs[:gradient]
    J_T_func = wrk.kwargs[:J_T]

    # calculate the functional only
    function f(F, G, pulsevals)
        @assert !isnothing(F)
        @assert isnothing(G)
        wrk.result.f_calls += 1
        @threadsif wrk.use_threads for k = 1:N
            copyto!(Ψ[k], wrk.objectives[k].initial_state)
            for n = 1:N_T  # `n` is the index for the time interval
                local (G, dt) = _fw_gen(pulsevals, k, n, wrk)
                propstep!(Ψ[k], G, dt, wrk.fw_prop_wrk[k])
            end
            τ[k] = dot(wrk.objectives[k].target_state, Ψ[k])
        end
        return J_T_func(Ψ, wrk.objectives; τ=τ)
    end

    # calculate the functional and the gradient
    function fg!(F, G, pulsevals)
        if isnothing(G)  # functional only
            return f(F, G, pulsevals)
        end
        wrk.result.fg_calls += 1
        # backward propagation of states
        @threadsif wrk.use_threads for k = 1:N
            copyto!(χ[k], wrk.objectives[k].target_state)
            write_to_storage!(X[k], N_T+1, χ[k])
            for n = N_T:-1:1
                local (G, dt) = _bw_gen(pulsevals, k, n, wrk)
                propstep!(χ[k], G, dt, wrk.bw_prop_wrk[k])
                write_to_storage!(X[k], n, χ[k])
            end
        end
        # forward propagation of gradients
        @threadsif wrk.use_threads for k = 1:N
            resetgradvec!(Ψ̃[k], wrk.objectives[k].initial_state)
            for n = 1:N_T  # `n` is the index for the time interval
                local (G̃, dt) = _fw_gradgen(pulsevals, k, n, wrk)
                propstep!(Ψ̃[k], G̃, dt, wrk.grad_prop_wrk[k])
                get_from_storage!(χ[k], X[k], n+1)
                for l = 1:L
                    ∇τ[k][l, n] = dot(χ[k], Ψ̃[k].grad_states[l])
                end
                resetgradvec!(Ψ̃[k])
            end
            τ[k] = dot(wrk.objectives[k].target_state, Ψ̃[k].state)
        end
        gradfunc!(G, τ, ∇τ)
        # TODO: set wrk.result.states
        return J_T_func([Ψ̃[k].state for k in 1:N], wrk.objectives; τ=τ)
    end

    optimizer = get_optimizer(wrk)
    res = run_optimizer(optimizer, wrk, fg!, info_hook, check_convergence!)
    finalize_result!(wrk, res)
    return wrk.result

end

function get_optimizer_optim_lbfgs(wrk) # TODO: get rid of this (not called)
    kwargs = wrk.kwargs
    lbfgs_kwargs = Dict{Symbol, Any}()
    lbfgs_keys = (:memory_length, :alphaguess, :linesearch, :P, :precond,
                    :manifold, :scaleinvH0)
    # TODO: get if of optim.jl-specific keywords: we'll default to LBFGS, and
    # if you want to use Optim.jl, you'll have to pass in a fully initialized
    # optimizer
    for key in lbfgs_keys
        if key in keys(kwargs)
            if :optimizer in keys(kwargs)
                @warn "keyword argument $(String(key)) will be ignored because due to custom optimizer"
            end
            val = kwargs[key]
            (key == :memory_length) && (key = :m)
            lbfgs_kwargs[key] = val
        end
    end
    optimizer = get(kwargs, :optimizer, Optim.LBFGS(;lbfgs_kwargs...))
    return optimizer
end

function get_optimizer(wrk)
    n = length(wrk.pulsevals)
    m = 10 # TODO: kwarg for number of limited memory corrections
    optimizer = get(wrk.kwargs, :optimizer, LBFGSB.L_BFGS_B(n, m))
    return optimizer
end

# The dynamical generator for the forward propagation (functional evaluation)
function _fw_gen(pulse_vals, k, n, wrk)
    vals_dict = wrk.vals_dict[k]
    t = wrk.result.tlist
    L = length(wrk.controls)
    ϵ = reshape(pulse_vals, L, :)
    for (l, control) in enumerate(wrk.controls)
        vals_dict[control] = ϵ[l, n]
    end
    dt = t[n+1] - t[n]
    evalcontrols!(wrk.G[k], wrk.objectives[k].generator, vals_dict)
    return wrk.G[k], dt
end

# The dynamical generator for the gradient propagation
function _fw_gradgen(pulse_vals, k, n, wrk)
    vals_dict = wrk.vals_dict[k]
    t = wrk.result.tlist
    L = length(wrk.controls)
    ϵ = reshape(pulse_vals, L, :)
    for (l, control) in enumerate(wrk.controls)
        vals_dict[control] = ϵ[l, n]
    end
    dt = t[n+1] - t[n]
    evalcontrols!(wrk.gradG[k], wrk.TDgradG[k], vals_dict)
    return wrk.gradG[k], dt
end

# The dynamical generator for the backward propagation
function _bw_gen(pulse_vals, k, n, wrk)
    vals_dict = wrk.vals_dict[k]
    L = length(wrk.controls)
    ϵ = reshape(pulse_vals, L, :)
    t = wrk.result.tlist
    for (l, control) in enumerate(wrk.controls)
        vals_dict[control] = ϵ[l, n]
    end
    dt = t[n+1] - t[n]
    evalcontrols!(wrk.G[k], wrk.adjoint_objectives[k].generator, vals_dict)
    return wrk.G[k], -dt
end


text:
export heterozygosity, samplehet

"""
    counthet(geno::T, allele::Int) where T<:GenoArray
    counthet(geno::T, allele::AbstractVector{U}) where T<:GenoArray where U<:Integer
Given a `GenoArray`, count the number of times `allele` appears in the
heterozygous state.
"""
function counthet(geno::T, allele::U) where T<:GenoArray where U<:Integer
    mapreduce(i -> _ishet(i, allele), +, skipmissing(geno))
end

function counthet(geno::T, allele::AbstractVector{U}) where T<:GenoArray where U<:Integer
    tmp = skipmissing(geno)
    isempty(tmp) && return fill(0, length(allele))
    map(_allele -> mapreduce(i -> _ishet(i, _allele), +, tmp), allele)
end

"""
    counthom(geno::T, allele::Int) where T<:GenoArray
Given a `GenoArray`, count the number of times `allele` appears in the
homozygous state.
"""
function counthom(geno::T, allele::U) where T<:GenoArray where U <: Signed
    mapreduce(i -> _ishom(i, allele), +, skipmissing(geno))
end

function counthom(geno::T, allele::AbstractVector{U}) where T<:GenoArray where U<:Integer
    tmp = skipmissing(geno)
    isempty(tmp) && return fill(0, length(allele))
    map(_allele -> mapreduce(i -> _ishom(i, _allele), +, tmp), allele)
end

"""
    _genediversitynei87(het_exp::Union{Missing,AbstractFloat}, het_obs::Union{Missing,AbstractFloat}, n::Union{Integer, Float64}, corr::Bool = true)
Calculate overall gene diversity with the adjustment/correction, use `corr = false` to ignore sample-size correction `* n/(n-1)`.

Hₜ = 1 −sum(pbar²ᵢ + Hₛ/(ñ * np) − Het_obs/(2ñ*np))
- _ñ_ is the number of genotypes for a locus for a population
- _np_ is the number of genotypes of a locus across all populations
    - i.e. sum(_ñ_)
- _pbar²_ is the observed homozygosity of a locus for that population
- _Hₛ_ is the within population gene diversity given by:
    - Hₛ =  ñ/(ñ-1) * (1 - sum(pbar²ᵢ - Het_observed / 2ñ))

(Nei M. (1987) Molecular Evolutionary Genetics. Columbia University Press).
"""
@inline function _genediversitynei87(het_exp::T, het_obs::T, n::Union{Integer,T}; corr::Bool = true) where T<: AbstractFloat
    corr_val = corr ? n/(n-1.0) : 1.0
    return @fastmath (het_exp - (het_obs/n/2.0)) * corr_val
end

@inline function _genediversitynei87(het_exp::AbstractFloat, het_obs::Missing, n::Union{Integer,AbstractFloat}; corr::Bool = true)
    return missing
end

@inline function _genediversitynei87(het_exp::Missing, het_obs::AbstractFloat, n::Union{Integer,AbstractFloat}; corr::Bool = true)
    return missing
end

@inline function _genediversitynei87(het_exp::Missing, het_obs::Missing, n::Union{Integer,AbstractFloat}; corr::Bool = true)
    return missing
end

"""
    _hetero_obs(data::T) where T <: GenoArray
Returns observed heterozygosity as a mean of the number of heterozygous genotypes, defined
as genotypes returning `true` for `_ishet()`. This is numerically feasible because
`true` values are mathematically represented as `1`, whereas `false` are represented
as `0`.
"""
@inline function _hetero_obs(data::T) where T <: GenoArray
    adjusted_vector = _ishet(data) |> skipmissing
    isempty(adjusted_vector) ? missing : mean(adjusted_vector)
end

"""
    _hetero_exp(allele_freqs::Vector{T}) where T <: GenoArray
Returns the expected heterozygosity of an array of genotypes,
calculated as 1 - sum of the squared allele frequencies.
"""
@inline function _hetero_exp(data::T) where T <: GenoArray
    all(ismissing.(data)) == true ? missing : 1.0 - mapreduce(i -> i^2, + , allelefreq_vec(data))
end

"""
    heterozygosity(data::PopData; by::String = "locus")
Calculate observed and expected heterozygosity in a `PopData` object. For loci,
heterozygosity is calculated in the Nei fashion, such that heterozygosity is
calculated as the average over heterozygosity per locus per population.
### Modes
- `"locus"` or `"loci"` : heterozygosity per locus (default)
- `"sample"` or `"ind"` or `"individual"` : heterozygosity per individual/sample
- `"population"` or `"pop"` : heterozygosity per population
## Example
heterozygosity(@nancycats, by = "population" )
"""
function heterozygosity(data::PopData; by::String = "locus")
    if by ∈ ["locus", "loci"]
        tmp = DataFrames.combine(
                groupby(data.genodata, [:locus, :population]),
                :genotype => nonmissing => :n_tmp,
                :genotype => _hetero_obs => :het_pop_obs,
                :genotype => _hetero_exp => :het_pop_exp
            )
        return DataFrames.combine(
                groupby(tmp, :locus),
                :n_tmp => sum => :n,
                :het_pop_obs => (h_o -> mean(skipmissing(h_o))) => :het_obs,
                :het_pop_exp => (h_e -> mean(skipmissing(h_e))) => :het_exp
            )

    elseif lowercase(by) ∈  ["sample", "ind", "individual"]
        return DataFrames.combine(
                groupby(data.genodata, :name),
                :genotype => nonmissing => :n,
                :genotype => _hetero_obs => :het_obs
            )

    elseif lowercase(by) ∈  ["pop", "population"]
        return DataFrames.combine(
                groupby(data.genodata, :population),
                :genotype => nonmissing => :n,
                :genotype => _hetero_obs => :het_obs,
                :genotype => _hetero_exp => :het_exp
            )
    else
        error("please specify by = \"locus\", \"sample\", or \"population\"")
    end
end

#NOTE this is not intended to be performant. It's a convenience function.
"""
    samplehet(data::PopData, individual::String)
Calculate the observed heterozygosity for an individual in a `PopData` object.
"""
@inline function samplehet(data::PopData, individual::String)
    individual ∉ samplenames(data) && throw(ArgumentError("$individual not found in data"))
    data.genodata[data.genodata.name .== individual, :genotype] |> _hetero_obs
end


text:
"""
    known_first(::Type{T})

If `first` of an instance of type `T` is known at compile time, return it.
Otherwise, return `nothing`.

@test isnothing(known_first(typeof(1:4)))
@test isone(known_first(typeof(Base.OneTo(4))))
"""
known_first(x) = known_first(typeof(x))
known_first(::Type{T}) where {T} = nothing
known_first(::Type{Base.OneTo{T}}) where {T} = one(T)
known_first(::Type{T}) where {T<:Base.Slice} = known_first(parent_type(T))

"""
    known_last(::Type{T})

If `last` of an instance of type `T` is known at compile time, return it.
Otherwise, return `nothing`.

@test isnothing(known_last(typeof(1:4)))
using StaticArrays
@test known_last(typeof(SOneTo(4))) == 4
"""
known_last(x) = known_last(typeof(x))
known_last(::Type{T}) where {T} = nothing
known_last(::Type{T}) where {T<:Base.Slice} = known_last(parent_type(T))

"""
    known_step(::Type{T})

If `step` of an instance of type `T` is known at compile time, return it.
Otherwise, return `nothing`.

@test isnothing(known_step(typeof(1:0.2:4)))
@test isone(known_step(typeof(1:4)))
"""
known_step(x) = known_step(typeof(x))
known_step(::Type{T}) where {T} = nothing
known_step(::Type{<:AbstractUnitRange{T}}) where {T} = one(T)

# add methods to support ArrayInterface

"""
    OptionallyStaticUnitRange(start, stop) <: AbstractUnitRange{Int}

This range permits diverse representations of arrays to communicate common information
about their indices. Each field may be an integer or `Val(<:Integer)` if it is known
at compile time. An `OptionallyStaticUnitRange` is intended to be constructed internally
from other valid indices. Therefore, users should not expect the same checks are used
to ensure construction of a valid `OptionallyStaticUnitRange` as a `UnitRange`.
"""
struct OptionallyStaticUnitRange{F <: Integer, L <: Integer} <: AbstractUnitRange{Int}
  start::F
  stop::L

  function OptionallyStaticUnitRange(start, stop)
    if eltype(start) <: Int
      if eltype(stop) <: Int
        return new{typeof(start),typeof(stop)}(start, stop)
      else
        return OptionallyStaticUnitRange(start, Int(stop))
      end
    else
      return OptionallyStaticUnitRange(Int(start), stop)
    end
  end

  function OptionallyStaticUnitRange{F,L}(x::AbstractRange) where {F,L}
    if step(x) == 1
      return OptionallyStaticUnitRange(static_first(x), static_last(x))
    else
        throw(ArgumentError("step must be 1, got $(step(r))"))
    end
  end

  function OptionallyStaticUnitRange(x::AbstractRange)
    if step(x) == 1
      return OptionallyStaticUnitRange(static_first(x), static_last(x))
    else
        throw(ArgumentError("step must be 1, got $(step(r))"))
    end
  end
end

Base.first(r::OptionallyStaticUnitRange) = r.start
Base.step(::OptionallyStaticUnitRange) = StaticInt(1)
Base.last(r::OptionallyStaticUnitRange) = r.stop

known_first(::Type{<:OptionallyStaticUnitRange{StaticInt{F}}}) where {F} = F
known_step(::Type{<:OptionallyStaticUnitRange}) = 1
known_last(::Type{<:OptionallyStaticUnitRange{<:Any,StaticInt{L}}}) where {L} = L

"""
    OptionallyStaticStepRange(start, step, stop) <: OrdinalRange{Int,Int}

Similarly to [`OptionallyStaticUnitRange`](@ref), `OptionallyStaticStepRange` permits
a combination of static and standard primitive `Int`s to construct a range. It
specifically enables the use of ranges without a step size of 1. It may be constructed
through the use of `OptionallyStaticStepRange` directly or using static integers with
the range operator (i.e., `:`).

```julia
julia> using ArrayInterface

julia> x = ArrayInterface.StaticInt(2);

julia> x:x:10
ArrayInterface.StaticInt{2}():ArrayInterface.StaticInt{2}():10

julia> ArrayInterface.OptionallyStaticStepRange(x, x, 10)
ArrayInterface.StaticInt{2}():ArrayInterface.StaticInt{2}():10
```
"""
struct OptionallyStaticStepRange{F <: Integer, S <: Integer, L <: Integer} <: OrdinalRange{Int,Int}
  start::F
  step::S
  stop::L

  function OptionallyStaticStepRange(start, step, stop)
    if eltype(start) <: Int
      if eltype(stop) <: Int
        lst = _steprange_last(start, step, stop)
        return new{typeof(start),typeof(step),typeof(lst)}(start, step, lst)
      else
        return OptionallyStaticStepRange(start, step, Int(stop))
      end
    else
      return OptionallyStaticStepRange(Int(start), step, stop)
    end
  end

  function OptionallyStaticStepRange(x::AbstractRange)
    return OptionallyStaticStepRange(static_first(x), static_step(x), static_last(x))
  end
end

# to make StepRange constructor inlineable, so optimizer can see `step` value
@inline function _steprange_last(start::StaticInt, step::StaticInt, stop::StaticInt)
    return StaticInt(_steprange_last(Int(start), Int(step), Int(stop)))
end
@inline function _steprange_last(start::Integer, step::StaticInt, stop::StaticInt)
    if step === one(step)
        # we don't need to check the `stop` if we know it acts like a unit range
        return stop
    else
        return _steprange_last(start, Int(step), Int(stop))
    end
end
@inline function _steprange_last(start::Integer, step::Integer, stop::Integer)
  z = zero(step)
  if step === z
    throw(ArgumentError("step cannot be zero"))
  else
    if stop == start
        return Int(stop)
    else
      if step > z
        if stop > start
          return stop - Int(unsigned(stop - start) % step)
        else
          return Int(start - one(start))
        end
      else
        if stop > start
          return Int(start + one(start))
        else
          return stop + Int(unsigned(start - stop) % -step)
        end
      end
    end
  end
end
Base.first(r::OptionallyStaticStepRange) = r.start
Base.step(r::OptionallyStaticStepRange) = r.step
Base.last(r::OptionallyStaticStepRange) = r.stop

known_first(::Type{<:OptionallyStaticStepRange{StaticInt{F}}}) where {F} = F
known_step(::Type{<:OptionallyStaticStepRange{<:Any,StaticInt{S}}}) where {S} = S
known_last(::Type{<:OptionallyStaticStepRange{<:Any,<:Any,StaticInt{L}}}) where {L} = L

Base.:(:)(L::Integer, ::StaticInt{U}) where {U} = OptionallyStaticUnitRange(L, StaticInt(U))
Base.:(:)(::StaticInt{L}, U::Integer) where {L} = OptionallyStaticUnitRange(StaticInt(L), U)
Base.:(:)(::StaticInt{L}, ::StaticInt{U}) where {L,U} = OptionallyStaticUnitRange(StaticInt(L), StaticInt(U))
Base.:(:)(::StaticInt{F}, ::StaticInt{S}, ::StaticInt{L}) where {F,S,L} = OptionallyStaticStepRange(StaticInt(F), StaticInt(S), StaticInt(L))
Base.:(:)(start::Integer, ::StaticInt{S}, ::StaticInt{L}) where {S,L} = OptionallyStaticStepRange(start, StaticInt(S), StaticInt(L))
Base.:(:)(::StaticInt{F}, ::StaticInt{S}, stop::Integer) where {F,S} = OptionallyStaticStepRange(StaticInt(F), StaticInt(S), stop)
Base.:(:)(::StaticInt{F}, step::Integer, ::StaticInt{L}) where {F,L} = OptionallyStaticStepRange(StaticInt(F), step, StaticInt(L))
Base.:(:)(start::Integer, step::Integer, ::StaticInt{L}) where {L} = OptionallyStaticStepRange(start, step, StaticInt(L))
Base.:(:)(start::Integer, ::StaticInt{S}, stop::Integer) where {S} = OptionallyStaticStepRange(start, StaticInt(S), stop)
Base.:(:)(::StaticInt{F}, step::Integer, stop::Integer) where {F} = OptionallyStaticStepRange(StaticInt(F), step, stop)
Base.:(:)(::StaticInt{F}, ::StaticInt{1}, ::StaticInt{L}) where {F,L} = OptionallyStaticUnitRange(StaticInt(F), StaticInt(L))
Base.:(:)(start::Integer, ::StaticInt{1}, ::StaticInt{L}) where {L} = OptionallyStaticUnitRange(start, StaticInt(L))
Base.:(:)(::StaticInt{F}, ::StaticInt{1}, stop::Integer) where {F} = OptionallyStaticUnitRange(StaticInt(F), stop)
Base.:(:)(start::Integer, ::StaticInt{1}, stop::Integer) = OptionallyStaticUnitRange(start, stop)

function Base.isempty(r::OptionallyStaticUnitRange)
  if known_first(r) === oneunit(eltype(r))
    return unsafe_isempty_one_to(last(r))
  else
    return unsafe_isempty_unit_range(first(r), last(r))
  end
end

function Base.isempty(r::OptionallyStaticStepRange)
  return (r.start != r.stop) & ((r.step > zero(r.step)) != (r.stop > r.start))
end

unsafe_isempty_one_to(lst) = lst <= zero(lst)
unsafe_isempty_unit_range(fst, lst) = fst > lst

unsafe_length_one_to(lst::Int) = lst
unsafe_length_one_to(::StaticInt{L}) where {L} = L

@inline function unsafe_length_step_range(start::Int, step::Int, stop::Int)
    if step > 1
        return Base.checked_add(Int(div(unsigned(stop - start), step)), 1)
    elseif step < -1
        return Base.checked_add(Int(div(unsigned(start - stop), -step)), 1)
    elseif step > 0
        return Base.checked_add(Int(div(Base.checked_sub(stop, start), step)), 1)
    else
        return Base.checked_add(Int(div(Base.checked_sub(start, stop), -step)), 1)
    end
end

@propagate_inbounds function Base.getindex(r::OptionallyStaticUnitRange, s::AbstractUnitRange{<:Integer})
    @boundscheck checkbounds(r, s)
    f = static_first(r)
    fnew = f - one(f)
    return (fnew + static_first(s)):(fnew + static_last(s))
end

@propagate_inbounds function Base.getindex(r::OptionallyStaticUnitRange, i::Integer)
  if known_first(r) === oneunit(eltype(r))
    return get_index_one_to(r, i)
  else
    return get_index_unit_range(r, i)
  end
end

@inline function get_index_one_to(r, i)
  @boundscheck if ((i < 1) || (i > last(r)))
    throw(BoundsError(r, i))
  end
  return convert(eltype(r), i)
end

@inline function get_index_unit_range(r, i)
  val = first(r) + (i - 1)
  @boundscheck if (i < 1) || val > last(r)
    throw(BoundsError(r, i))
  end
  return convert(eltype(r), val)
end

@inline _try_static(::StaticInt{N}, ::StaticInt{N}) where {N} = StaticInt{N}()
@inline _try_static(::StaticInt{M}, ::StaticInt{N}) where {M, N} = @assert false "Unequal Indices: StaticInt{$M}() != StaticInt{$N}()"
@propagate_inbounds function _try_static(::StaticInt{N}, x) where {N}
  @boundscheck begin
     @assert N == x "Unequal Indices: StaticInt{$N}() != x == $x"
  end
  return StaticInt{N}()
end
@propagate_inbounds function _try_static(x, ::StaticInt{N}) where {N}
  @boundscheck begin
     @assert N == x "Unequal Indices: x == $x != StaticInt{$N}()"
  end
  return StaticInt{N}()
end
@propagate_inbounds function _try_static(x, y)
  @boundscheck begin
    @assert x == y "Unequal Indices: x == $x != $y == y"
  end
  return x
end

###
### length
###
@inline function known_length(::Type{T}) where {T<:OptionallyStaticUnitRange}
  fst = known_first(T)
  lst = known_last(T)
  if fst === nothing || lst === nothing
    return nothing
  else
    if fst === oneunit(eltype(T))
      return unsafe_length_one_to(lst)
    else
      return unsafe_length_unit_range(fst, lst)
    end
  end
end

@inline function known_length(::Type{T}) where {T<:OptionallyStaticStepRange}
  fst = known_first(T)
  stp = known_step(T)
  lst = known_last(T)
  if fst === nothing || stp === nothing || lst === nothing
    return nothing
  else
    if stp === 1
      if fst === oneunit(eltype(T))
        return unsafe_length_one_to(lst)
      else
        return unsafe_length_unit_range(fst, lst)
      end
    else
      return unsafe_length_step_range(fst, stp, lst)
    end
  end
end

function Base.length(r::OptionallyStaticUnitRange)
  if isempty(r)
    return 0
  else
    if known_first(r) === 1
      return unsafe_length_one_to(last(r))
    else
      return unsafe_length_unit_range(first(r), last(r))
    end
  end
end

function Base.length(r::OptionallyStaticStepRange)
  if isempty(r)
    return 0
  else
    if known_step(r) === 1
      if known_first(r) === 1
        return unsafe_length_one_to(last(r))
      else
        return unsafe_length_unit_range(first(r), last(r))
      end
    else
        return unsafe_length_step_range(Int(first(r)), Int(step(r)), Int(last(r)))
    end
  end
end

unsafe_length_unit_range(start::Integer, stop::Integer) = Int((stop - start) + 1)

function Base.AbstractUnitRange{T}(r::OptionallyStaticUnitRange) where {T}
    if T <: Int
        return r
    else
        if known_first(r) === 1 && T <: Integer
            return OneTo{T}(last(r))
        else
            return UnitRange{T}(first(r), last(r))
        end
    end
end

const OptionallyStaticRange = Union{<:OptionallyStaticUnitRange,<:OptionallyStaticStepRange}

Base.to_shape(x::OptionallyStaticRange) = length(x)
Base.to_shape(x::Slice{T}) where {T<:OptionallyStaticRange} = length(x)

function Base.axes(S::Slice{T}) where {T<:OptionallyStaticRange}
    if known_first(T) === 1 && known_step(T) === 1
        return (S.indices,)
    else
        return (Base.IdentityUnitRange(S.indices),)
    end
end

function Base.axes1(S::Slice{T}) where {T<:OptionallyStaticRange}
    if known_first(T) === 1 && known_step(T) === 1
        return S.indices
    else
        return Base.IdentityUnitRange(S.indices)
    end
end

Base.:(-)(r::OptionallyStaticRange) = -static_first(r):-static_step(r):-static_last(r)

"""
    indices(x[, d])

Given an array `x`, this returns the indices along dimension `d`. If `x` is a tuple
of arrays, then the indices corresponding to dimension `d` of all arrays in `x` are
returned. If any indices are not equal along dimension `d`, an error is thrown. A
tuple may be used to specify a different dimension for each array. If `d` is not
specified, then the indices for visiting each index of `x` are returned.
"""
@inline function indices(x)
  inds = eachindex(x)
  if inds isa AbstractUnitRange && eltype(inds) <: Integer
    return Base.Slice(OptionallyStaticUnitRange(inds))
  else
    return inds
  end
end

function indices(x::Tuple)
  inds = map(eachindex, x)
  return reduce(_pick_range, inds)
end

@inline indices(x, d) = indices(axes(x, d))

@inline function indices(x::Tuple{Vararg{Any,N}}, dim) where {N}
  inds = map(x_i -> indices(x_i, dim), x)
  return reduce(_pick_range, inds)
end

@inline function indices(x::Tuple{Vararg{Any,N}}, dim::Tuple{Vararg{Any,N}}) where {N}
  inds = map(indices, x, dim)
  return reduce(_pick_range, inds)
end

@inline function _pick_range(x, y)
  fst = _try_static(static_first(x), static_first(y))
  lst = _try_static(static_last(x), static_last(y))
  return Base.Slice(OptionallyStaticUnitRange(fst, lst))
end


text:
# Find link to SQLite shared library
let
    global sqlite3_lib
    local lib
    succeeded = false
    @linux_only begin
        lib_choices = ["libsqlite3"]
    end
    @windows_only begin
        if WORD_SIZE == 64
            lib_choices = [Pkg.dir("SQLite", "lib", "sqlite3-64")]
        else
            lib_choices = [Pkg.dir("SQLite", "lib", "sqlite3")]
        end
    end
    @osx_only begin
        lib_choices = [
            "/usr/local/lib/libsqlite3.dylib",
            "/usr/lib/libsqlite3.dylib",
            Pkg.dir("SQLite", "lib", "libsqlite3")
        ]
    end
    for lib in lib_choices
        try
            dlopen(lib)
            succeeded = true
            break
        end
    end
    succeeded || error("SQLite library not found")
    @eval const sqlite3_lib = $lib
end

# Return codes
const SQLITE_OK = 0 # Successful result
const SQLITE_ROW = 100 # sqlite3_step() has another row ready
const SQLITE_DONE = 101  # sqlite3_step() has finished executing

# Types
const SQLITE_INTEGER = 1
const SQLITE_FLOAT = 2
const SQLITE3_TEXT = 3
const SQLITE_BLOB = 4
const SQLITE_NULL = 5

# I/O settings
const SQLITE_OPEN_READWRITE = 0x00000002


text:
using Pkg
Pkg.activate(joinpath(Pkg.devdir(), "MLCourse"))
using MLCourse, CSV, MLJ, DataFrames, MLJLinearModels, Random, Flux, MLJFlux, MLJXGBoostInterface, Plots
include("./Data_Processing.jl")

#Load the data
drop = generate(option = "drop", std = "fale", valid = "false", test = "true");
med = generate(option = "med", std = "false", valid = "false", test = "true");
drop_std = generate(option = "drop", std = "true", valid = "false", test = "true");
med_std = generate(option = "med", std = "true", valid = "false", test = "true");

#First scheme of possible neural network on med data set and drop data set, standardized or not.
Random.seed!(2711);
NN_1 = machine(NeuralNetworkClassifier(builder = MLJFlux.Short(n_hidden = 200, dropout = 0.1, σ = relu),
                                                batch_size = 1000, epochs = 100),
                                                select(drop_std.train[:,:], Not(:precipitation_nextday)),
                                                drop_std.train.precipitation_nextday) |> fit!;

MLJ.auc(predict(NN_1, select(drop_std.test[:,:], Not(:precipitation_nextday))), drop_std.test.precipitation_nextday)
# AUC_DROP_STD = 0.912
# Relu seems to perform better than sigmoid function.
# Drop seems to perform better than med.
# No particular differences between standardized and normal data.

#Tuned model
#Use this section to try multiple model with changes in some hyper-parametrs in order to see if there is some increase in the AUC
Random.seed!(2711)
model = NeuralNetworkClassifier( builder = MLJFlux.Short(n_hidden = 100, dropout = 0.1, σ = relu),
    epochs = 100,
    batch_size = 100);

result=[];
for _ in 1:10 #Change the second number if you want to try some changes
    NN_2 = machine(model, select(drop_std.train, Not(:precipitation_nextday)), drop_std.train.precipitation_nextday)
    fit!(NN_2)
    AUC_drop = MLJ.auc(predict(NN_2, select(drop_std.test[:,:], Not(:precipitation_nextday))), drop_std.test.precipitation_nextday)
    push!(result, AUC_drop)
    model.batch_size = model.batch_size + 100
end
scatter(1:10, result, ylims = [0.9, 0.95])

# Best model found in the previous part
Random.seed!(2711)
best_model = NeuralNetworkClassifier( builder = MLJFlux.Short(n_hidden = 200, dropout = 0.1, σ = relu),
    epochs = 300,
    batch_size = 1000)
best_NN = machine(best_model, select(drop_std.train, Not(:precipitation_nextday)), drop_std.train.precipitation_nextday)
fit!(best_NN, verbosity = 2)
MLJ.auc(predict(best_NN, select(drop_std.test[:,:], Not(:precipitation_nextday))), drop_std.test.precipitation_nextday)

#Write in the submission file with a machine trained on all the data
train, test = generate(option = "drop", std = "true", valid = "false", test = "false");
best_mach = machine(best_model, select(train, Not(:precipitation_nextday)), train.precipitation_nextday);
fit!(best_mach)
sample = CSV.read(joinpath(@__DIR__, "data", "sample_submission.csv"), DataFrame);
pred = pdf.(predict(best_mach, test), true);
sample.precipitation_nextday = pred;
CSV.write(joinpath(@__DIR__, "results", "MLP_submission.csv"), sample);

max_stars_repo_name: clara-rossignol/Will-it-rain-tomorrow

max_issues_repo_name: clara-rossignol/Will-it-rain-tomorrow

max_forks_repo_name: clara-rossignol/Will-it-rain-tomorrow


text:
using Documenter, MathOptInterface

"""
Pass `julia docs/make.jl --fix` to rebuild the doctests.
"""
const _FIX = findfirst(isequal("--fix"), ARGS) !== nothing

makedocs(
    sitename = "MathOptInterface",
    format = Documenter.HTML(
        # See https://github.com/JuliaDocs/Documenter.jl/issues/868
        prettyurls = get(ENV, "CI", nothing) == "true",
        mathengine = Documenter.MathJax2(),
        collapselevel = 1,
    ),
    strict = true,
    modules = [MathOptInterface],
    checkdocs = :exports,
    doctest = _FIX ? :fix : true,
    pages = [
        "Introduction" => "index.md",
        "Background" => [
            "background/motivation.md",
            "background/duality.md",
        ],
        "Manual" => [
            "manual/standard_form.md",
            "manual/models.md",
            "manual/variables.md",
            "manual/constraints.md",
            "manual/solutions.md",
            "manual/modification.md",
            "manual/example.md",
            "manual/implementing.md",
        ],
        "API Reference" => [
            "reference/standard_form.md",
            "reference/models.md",
            "reference/variables.md",
            "reference/constraints.md",
            "reference/modification.md",
            "reference/nonlinear.md",
            "reference/callbacks.md",
            "reference/errors.md",
        ],
        "Submodules" => [
            "Benchmarks" => [
                "Overview" => "submodules/Benchmarks/overview.md",
                "API Reference" => "submodules/Benchmarks/reference.md",
            ],
            "Bridges" => [
                "Overview" => "submodules/Bridges/overview.md",
                "Implementation" => "submodules/Bridges/implementation.md",
                "API Reference" => "submodules/Bridges/reference.md",
            ],
            "FileFormats" => [
                "Overview" => "submodules/FileFormats/overview.md",
                "API Reference" => "submodules/FileFormats/reference.md",
            ],
            "Utilities" => [
                "Overview" => "submodules/Utilities/overview.md",
                "API Reference" => "submodules/Utilities/reference.md",
            ],
            "Test" => [
                "Overview" => "submodules/Test/overview.md",
            ],
        ],
    ],
)

deploydocs(
    push_preview = true,
    repo = "github.com/jump-dev/MathOptInterface.jl.git",
)


text:
#
# Correctness Tests
#

using LsqFit

my_tests = [ "curve_fit.jl"]

println("Running tests:")

for my_test in my_tests
    println(" * $(my_test)")
    include(my_test)
end


text:
module MyResultsManager

import MyTools;

include("Core/Core.jl");
include("Interface/Interface.jl");

end # module

max_stars_repo_name: josePereiro/MyResultsManager

max_issues_repo_name: josePereiro/MyResultsManager

max_forks_repo_name: josePereiro/MyResultsManager


text:
module BitInformation

    export shave, set_one, groom, halfshave,
        shave!, set_one!, groom!, halfshave!, round!

    export bittranspose, bitbacktranspose,
        xor_delta, unxor_delta, xor_delta!, unxor_delta!,
        signed_exponent

    export bitinformation, mutual_information, redundancy, bitpattern_entropy,
        bitcount, bitcount_entropy, bitpaircount, bit_condprobability,
        bit_condentropy

    import StatsBase.entropy

    include("which_uint.jl")
    include("bitstring.jl")
    include("bittranspose.jl")
    include("round_nearest.jl")
    include("shave_set_groom.jl")
    include("xor_delta.jl")
    include("signed_exponent.jl")
    include("bit_information.jl")
    include("mutual_information.jl")
    include("bitpattern_entropy.jl")
    include("binomial.jl")

end


text:
    global graphvizsetup_is_working = false
    @info("Testing Graphviz Setup")
    @testset verbose = true "GraphViz Setup" begin
        if haskey(ENV,"graphvizdot")
            #if key exists this takes precedence
            global graphvizexe = ENV["graphvizdot"]
            #=
                ENV["graphvizdot"]="c:\\program files\\Graphviz\\bin\\dot.exe"
            =#
        else
            global graphvizexe = "dot" #linux & appleOS assume it is in in path
        end

        gcmd = `$(graphvizexe) "-V"`
        try
            gres = run(gcmd)
            @test gres.exitcode == 0
            @show gres.exitcode
            global graphvizsetup_is_working = true
            @test true #graphviz setup works
        catch e
            if Sys.iswindows()
                try
                    global graphvizexe = "c:\\program files\\Graphviz\\bin\\dot.exe"
                    gcmd = `$(graphvizexe) "-V"`
                    gres = run(gcmd)
                    @test gres.exitcode == 0
                    @show gres.exitcode
                    global graphvizsetup_is_working = true
                    @test true #graphviz setup works
                catch eiwn
                    @warn("DTM: unable to invoke gcmd")
                    @show gcmd
                    print(ewin)
                    @test false
                end
            else
                @warn("DTM: unable to invoke gcmd")
                @show gcmd
                print(e)
                @test false
            end
        end
    end


text:
module TestCollect

using StructuredQueries
using Base.Test

type MyData
    count::Int
end
StructuredQueries._collect(src::Tuple{MyData}, q::StructuredQueries.Node) =
    MyData(first(src).count+1)
StructuredQueries.prepare(src::MyData) = src
Base.collect(src::MyData) = src

src = MyData(0)

res = collect(@with src select(i))
@test res.count == 1

res = collect(@with src begin
    filter(i * j > k)
    select(l)
end)
@test res.count == 2

res = collect(@with src begin
    filter(sepal_length > 5.0)
    groupby(species)
    summarize(avg = mean(petal_width))
end)
@test res.count == 3

end


text:
### Julia OpenStreetMap Package ###
### MIT License                 ###
### Copyright 2014              ###

### Find the nearest node to a given location ###
function nearestNode{T<:@compat(Union{ENU,ECEF})}(nodes::Dict{Int,T}, loc::T)
    min_dist = Inf
    best_ind = 0

    for (key, node) in nodes
        dist = distance(node, loc)
        if dist < min_dist
            min_dist = dist
            best_ind = key
        end
    end

    return best_ind
end

### Find nearest node in a list of nodes ###
function nearestNode{T<:@compat(Union{ENU,ECEF})}(nodes::Dict{Int,T},
                                         loc::T,
                                         node_list::Vector{Int})
    min_dist = Inf
    best_ind = 0

    for ind in node_list
        dist = distance(nodes[ind], loc)
        if dist < min_dist
            min_dist = dist
            best_ind = ind
        end
    end

    return best_ind
end

### Find nearest node serving as a vertex in a routing network ###
function nearestNode{T<:@compat(Union{ENU,ECEF})}( nodes::Dict{Int,T},
                                          loc::T,
                                          network::Network )
    return nearestNode(nodes,loc,collect(keys(network.v)))
end

### Find all nodes within range of a location ###
function nodesWithinRange{T<:@compat(Union{ENU,ECEF})}(nodes::Dict{Int,T},
                                             loc::T,
                                             range::Float64=Inf)
    if range == Inf
        return keys(nodes)
    end
    indices = Int[]
    for (key, node) in nodes
        dist = distance(node, loc)
        if dist < range
            push!(indices, key)
        end
    end
    return indices
end

### Find nodes within range of a location using a subset of nodes ###
function nodesWithinRange{T<:@compat(Union{ENU,ECEF})}(nodes::Dict{Int,T},
                                              loc::T,
                                              node_list::Vector{Int},
                                              range::Float64=Inf)
    if range == Inf
        return node_list
    end
    indices = Int[]
    for ind in node_list
        dist = distance(nodes[ind], loc)
        if dist < range
            push!(ind, key)
        end
    end
    return indices
end

### Find vertices of a routing network within range of a location ###
function nodesWithinRange{T <: @compat(Union{ENU,ECEF}) }(nodes::Dict{Int,T},
                                              loc::T,
                                              network::Network,
                                              range::Float64=Inf)
    return nodesWithinRange(nodes,loc,collect(keys(network.v)),range)
end

### Add a new node ###
function addNewNode!{T <: @compat(Union{LLA,ENU}) }(nodes::Dict{Int,T},
                                        loc::T,
                                        start_id::Int = reinterpret(@compat(Int), hash(loc)) )
    id = start_id
    while id <= typemax(Int)
        if !haskey(nodes, id)
            nodes[id] = loc
            return id
        end
        id += 1
    end

    msg = "Unable to add new node to map, $(typemax(Int)) nodes is the current limit."
    throw(error(msg))
end

### Compute centroid of list of nodes ###
function centroid{T<:@compat(Union{LLA,ENU})}(nodes::Dict{Int,T}, node_list::Vector{Int})
    sum_1 = 0
    sum_2 = 0
    sum_3 = 0

    for k = 1:length(node_list)
        if typeof(nodes) == Dict{Int,LLA}
            sum_1 += nodes[node_list[k]].lat
            sum_2 += nodes[node_list[k]].lon
            sum_3 += nodes[node_list[k]].alt
        else
            sum_1 += nodes[node_list[k]].east
            sum_2 += nodes[node_list[k]].north
            sum_3 += nodes[node_list[k]].up
        end
    end

    if typeof(nodes) == Dict{Int,LLA}
        return LLA(sum_1/length(node_list),sum_2/length(node_list),sum_3/length(node_list))
    else
        return ENU(sum_1/length(node_list),sum_2/length(node_list),sum_3/length(node_list))
    end
end


text:
## Example script for building surrogate models using VLM analysis, needs Surrogates and DataFrames installed
using Surrogates
using DataFrames
using AeroMDAO

## Wing
wing_foils = Foil.(fill(naca4((0,0,1,2)), 2))
wing       = Wing(foils     = wing_foils,
                  chords    = [1.0, 0.6],
                  twists    = [0.0, 0.0],
                  spans     = [5.0],
                  dihedrals = [11.3],
                  sweep_LEs = [2.29]);
wing_mac   = mean_aerodynamic_center(wing)
wing_pos   = [0., 0., 0.]
print_info(wing, "Wing")

# Horizontal tail
htail = Wing(chords    = [0.7, 0.42],
             twists    = [0.0, 0.0],
             spans     = [1.25],
             dihedrals = [0.],
             sweep_LEs = [6.39])
htail_mac = mean_aerodynamic_center(htail)
htail_pos = [5., 0., 0.]
α_h_i     = 0.
print_info(htail, "Horizontal Tail")

# Vertical tail
vtail_foil = Foil(naca4((0,0,0,9)))
vtail = HalfWing(foils     = fill(vtail_foil, 2),
                 chords    = [0.7, 0.42],
                 twists    = [0.0, 0.0],
                 spans     = [1.0],
                 dihedrals = [0.],
                 sweep_LEs = [7.97])
vtail_mac = mean_aerodynamic_center(vtail) # NEEDS FIXING FOR ROTATION
vtail_pos = [5., 0., 0.]
print_info(vtail, "Vertical Tail")

## Panelling and assembly
wing_panels  = panel_wing(wing, [20], 10,
                          position = wing_pos
                         );
htail_panels = panel_wing(htail, [12], 12;
                          position	= htail_pos,
                          angle 	= deg2rad(α_h_i),
                          axis 	  	= [0., 1., 0.]
                         )
vtail_panels = panel_wing(vtail, [12], 10;
                          position 	= vtail_pos,
                          angle 	= π/2,
                          axis 	 	= [1., 0., 0.]
                         )

aircraft     = Dict(
                    "Wing"            => wing_panels,
                    "Horizontal Tail" => htail_panels,
                    "Vertical Tail"   => vtail_panels
                   )

## VLM setup
function vlm_analysis(aircraft, fs, ρ, ref, S, b, c, print = false)
    # Evaluate case
    data =  solve_case(aircraft, fs;
                       rho_ref   = ρ,
                       r_ref     = ref,
                       area_ref  = S,
                       span_ref  = b,
                       chord_ref = c,
                       print     = print
                      );

    # Get data
    nf_coeffs, ff_coeffs, CFs, CMs, horseshoe_panels, normals, horseshoes, Γs = data["Aircraft"]

    # Filter relevant data
    [ ff_coeffs[1:3]; nf_coeffs[4:6] ]
end

## Evaluate one case
# Case
x_w     = wing_pos + [ wing_mac[1], 0, 0 ]
ac_name = "My Aircraft"
ρ       = 1.225
ref     = x_w
V, α, β = 1.0, 0.0, 0.0
Ω       = [0.0, 0.0, 0.0]
fs      = Freestream(V, α, β, Ω)
S, b, c = projected_area(wing), span(wing), mean_aerodynamic_chord(wing)

@time coeffs = vlm_analysis(aircraft, fs, ρ, ref, S, b, c, true)

## Looping
αs, βs   = float.(-4:4), float.(-3:3)
fses     = [ Freestream(V, α, β, Ω) for α in αs, β in βs ]

## Evaluate cases
results = vlm_analysis.(Ref(aircraft), fses, ρ, Ref(ref), S, b, c, true);

## Generate DataFrame
αβs  = [ (α, β) for α in αs, β in βs ]
data = DataFrame([ (α, β, CD, CY, CL, Cl, Cm, Cn) for ((α, β), (CD, CY, CL, Cl, Cm, Cn)) in zip(αβs, results) ][:])
rename!(data, [:α, :β, :CD, :CY, :CL, :Cl, :Cm, :Cn])

##
using Plots
plotly(size = (1280, 720), markersize = 2, dpi = 300)

## The labels don't work for some reason, need to fix that later...
CD_plot = scatter(data[:,1], data[:,2], data[:,3], label = "CD")
CY_plot = scatter(data[:,1], data[:,2], data[:,4], label = "CY")
CL_plot = scatter(data[:,1], data[:,2], data[:,5], label = "CL")
Cl_plot = scatter(data[:,1], data[:,2], data[:,6], label = "Cl")
Cm_plot = scatter(data[:,1], data[:,2], data[:,7], label = "Cm")
Cn_plot = scatter(data[:,1], data[:,2], data[:,8], label = "Cn")

plots = [ CD_plot CY_plot CL_plot Cl_plot Cm_plot Cn_plot ]
plot(plots..., layout = (2, 3), xlabel = "α, ᵒ", ylabel = "β, ᵒ", zlabel = :nothing)
plot!()

## Surrogate training data
lower_bound, upper_bound = minimum.(eachcol(data[:,3])), maximum.(eachcol(data[:,3]))
zs = data[:,3]

## Generate surrogate
surrogate_model = LobachevskySurrogate(αβs, zs, lower_bound, upper_bound);

## Evaluate surrogate
num     = 100
xs, ys  = range(-4, 4, length = num), range(-3, 3, length = num)
xys     = [ (x, y) for x in xs, y in ys ]
surr_zs = surrogate_model.(xys)

## Create DataFrame
surr_data = DataFrame((x, y, z) for ((x, y), z) in zip(xys, surr_zs))
rename!(surr_data, [:α, :β, :CD])

## Plot surrogate
scatter(data[:,1], data[:,2], data[:,3], xlabel = "α, ᵒ", ylabel = "β, ᵒ", zlabel = "CD", label = "CD")
mesh3d!(surr_data[:,1], surr_data[:,2], surr_data[:,3], label = "Surrogate CD")

max_forks_repo_name: GodotMisogi/AeroMDAO


text:
__precompile__()

module StatGeochem

    using Reexport
    @reexport using NaNStatistics
    @reexport using StatGeochemBase

    # AVX vectorziation tools
    using LoopVectorization

    # General requirements
    using Statistics, DelimitedFiles, SpecialFunctions, Random
    using StatsBase: percentile, mean, std, ProbabilityWeights
    using ProgressMeter: @showprogress, Progress, update!
    using Interpolations
    include("utilities/System.jl")
    include("utilities/Import.jl")
    include("utilities/Resampling.jl")
    include("utilities/Changepoint.jl")

    include("resources/Chemistry.jl")
    include("utilities/Geochronology.jl")
    include("utilities/Geochemistry.jl")
    include("utilities/GIS.jl")
    include("utilities/Etc.jl")

    # Resources
    resourcepath = joinpath(homedir(),"resources")
    moduleresourcepath = joinpath(Base.source_dir(),"resources")
    export resourcepath, moduleresourcepath

    using FileIO, ImageIO, HDF5
    using Colors: Color, RGBX, RGB, N0f8
    include("resources/tc1/tc1.jl")
    include("resources/Crust1.jl")
    include("resources/Elevation.jl")
    include("resources/Seafloorage.jl")
    include("resources/PartitionCoefficients/PartitionCoefficients.jl")

    # Custom pretty printing for some types
    include("utilities/Display.jl")

end # module


text:
#
# need to set JULIA_NUM_THREADS
#
function digits_4(x::Int)::String
    if x < 10
        return string("000",x)
    elseif x<100
        return string("00",x)
    elseif x<1000
        return string("0",x)
    else
        string(x)
    end
end

function digits_3(x::Int)::String
    if x < 10
        return string("00",x)
    elseif x<100
        return string("0",x)
    else
        string(x)
    end
end

function digits_2(x::Int)::String
    if x < 10
        return string("0",x)
    else
        string(x)
    end
end

function main()
    filenames = readdir(pwd())

    Threads.@threads for i = 1:30
        p_file_name = string("ev",digits_2(i))
        mkdir(p_file_name)
        for j in filenames
            if isfile(j)
                if occursin(p_file_name,j)
                    mv(j,string(p_file_name,"/",j))
                end
            end
        end
    end
    auto_delete_empty_folder = true
    if auto_delete_empty_folder == true
        dir = readdir()
        for i in dir
            if isdir(i)
                if readdir(i) == []
                    rm(i)
                end
            end
        end
    end
end
main()

max_stars_repo_name: ditys/WebpMergeDecompose

max_issues_repo_name: ditys/WebpMergeDecompose

max_forks_repo_name: ditys/WebpMergeDecompose


text:
using ArgParse

include("BBTNet/models/retinaface.jl")
include("./BBTNet/utils/ImageReader.jl")
include("./BBTNet/utils/draw.jl")
include("configs.jl")

function parse_cmd()
    s = ArgParseSettings()
    @add_arg_table s begin
        "--image_path"
            help = "Path of the image to predict."
            arg_type = String
            required = true
        "--load_dir"
            help = "If there is a pretrained model, then the path of it."
            arg_type = String
            required = true
        "--save_dir", "-s"
            help = "Path for saving a the result after prediction."
            arg_type = String
            default = "./data/results/predicted.png"
            required = false
        "--mode", "-m"
            help = "Training mode: 0 for only baseline until FPN, 1 for full model, 2 for no cascaded structure."
            arg_type = Int
            default = 1
            required = false
        "--laterals", "-l"
            help = "How many lateral connections will be processed, either 3 or 5. 5 is needed for full model."
            arg_type = Int
            default = 5
            required = false
    end
    return parse_args(s)
end

function main()
    parsed_args = parse_cmd()
    img_path = parsed_args["image_path"]
    mode = parsed_args["mode"]
    scale_cnt = parsed_args["laterals"]
    num_anchors = scale_cnt == 3 ? 2 : 3
    anchor_info = scale_cnt == 3 ? lat_3_anchors : lat_5_anchors

    load_path = parsed_args["load_dir"]; save_path = parsed_args["save_dir"];
    model = RetinaFace(
        mode=mode, num_anchors=num_anchors, anchor_info=anchor_info, load_path=load_path, dtype=atype
    )
    print("[INFO] Model is loaded!\n")

    img, _ = read_img(img_path); reversed_img = reshape(reverse(img, dims=1), (size(img)..., 1));
    imgs_arr = convert(atype, permutedims(reversed_img, (3,2,1,4)))
    c, b, l = predict_image(model, imgs_arr, mode=mode); b[b .< 0] .= 0; l[l .< 0] .= 0;

    init_img = (img .+ avg_img) ./ 255
    # return draw_boxes_and_landmarks(init_img, b, l; len=640)

    if save_path !== nothing
        init_img = (img .+ avg_img) ./ 255
        plotted = draw_boxes_and_landmarks(init_img, b, l; conf=c, len=640)
        png(plotted, save_path)
    end

    println("# of persons found:", size(c, 1))

    return c, b, l
end

main()

max_stars_repo_name: barisbatuhan/retinaface

max_issues_repo_name: barisbatuhan/retinaface

max_forks_repo_name: barisbatuhan/retinaface


text:
# Explore a 1D Gaussian PDF using a basic MCMC chain
using ReversibleJumpMCMC
using Distributions
using Plots

# setup a MH structure that gets passed around inside RJMMCMC
struct MHStruct
    σ::Float32 # metropolis hasting jump size
end
mhs = MHStruct(1f0)

# setup a state type
mutable struct state1D
    x::Float32
end

# setup proposal and acceptance functions
function mypropose(mhs::MHStruct, s::state1D)
    teststate = state1D(s.x + mhs.σ * randn(Float32))
    return teststate, 0
end
function myaccept(mhs::MHStruct, s::state1D, teststate::state1D, vararg)
    σ = 5f0
    d = Normal(0f0, σ)
    α = pdf(d, teststate.x) / pdf(d, s.x)
    return α
end

# setup the RJMCMCStruct
burnin = Int32(100)
iterations = Int32(10000)
njumptypes = 1 # number of jump types for RJMCMC
jumpprobability = [1] # vector
proposalfuns = [mypropose] # vector
acceptfuns = [myaccept] # vector
rjs = ReversibleJumpMCMC.RJMCMCStruct(burnin, iterations, njumptypes, jumpprobability, proposalfuns, acceptfuns)

# initial state
state0 = state1D(0f0)

# run chain
rjc = ReversibleJumpMCMC.buildchain(rjs, mhs, state0)

# extract states
xchain = zeros(Float32, iterations)
for ii = 1:iterations
    xchain[ii] = rjc.states[ii].x
end

# plot found distributions
gr()
plt = histogram(xchain, normalize=true, xlabel="θ", ylabel="pdf(θ",label="mcmc dist.")

# true pdf
xvec = (-20:20)
σ = 5f0
d = Normal(0f0, σ)
plot!(plt,xvec,pdf(d, xvec),label="true pdf")


text:
## mustache template for ipynb
ipynb_tpl_v3 = mt"""
{"metadata": {
"language": "Julia",
 "name": "{{{TITLE}}}"
  },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
     {{{CELLS}}}
    ],
   "metadata": {}
  }
]
}
"""

ipynb_tpl_v4 = mt"""
{
  "cells": [
     {{{CELLS}}}
    ],
 "metadata": {
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.3"
  },
 "kernelspec": {
   "display_name": "Julia 1.3",
   "language": "julia",
   "name": "julia-1.3"
  }

 },
 "nbformat": 4,
 "nbformat_minor": 2

}
"""

not_needed = """

"""

const ipynb_tpl = ipynb_tpl_v4

## graphs. Don't want to dispatch, as packages are loaded in module, not global..
function render_winston(img)
    ## need a cell
    out = Dict()
    out["metadata"] = Dict()
    out["output_type"] = "execute_result"
    out["png"] = stringmime("image/png", img)
    out
end

function render_gadfly(img)
    ## need a cell
    info("render gadfly")
    imgfile = tempname()
    open(imgfile, "w") do io
        draw(PNG(io, 5inch, inch), img)
    end
    data = base64encode(read(imgfile, String))

    out = Dict()
#    out["metadata"] = Dict()
    out["data"] = Dict()
#    x = sprint(io -> tohtml(io, img))
#    x = split(x, "\n")
#    x = map(a -> a*"\n", x)
    #out["html"] = x
    out["data"]["image/png"] = data
    out["data"]["text/plain"] = ["Plot(...)"]
    out
end

function render_plotly(img)
    ## need a cell
    info("render plotly")
    imgfile = tempname()
    png(img, imgfile)
    data = base64encode(read(imgfile, String))

    out = Dict()
    out["data"] = Dict()
    out["data"]["image/png"] = data
    out["data"]["text/plain"] = ["Plot(...)"]
    out
end

function render_pyplot(img)
    info("render pyplot")
    out = Dict()
    out["metadata"] = Dict()
    out["output_type"] = "execute_result"
    out["png"] = stringmime("image/png", img)
    img[:clear]()
    out
end

function render_gr(img)
     info("render pyplot")
    out = Dict()
    out["metadata"] = Dict()
    out["output_type"] = "execute_result"
    out["png"] = stringmime("image/png", img)
    img[:clear]()
    out
end

## Main function to take a jmd file and turn into a ipynb file
function markdownToPynb(fname::AbstractString)

    dirnm, basenm = dirname(fname), basename(fname)
    newnm = replace(fname, r"[.].*" => ".ipynb")
    out = mdToPynb(fname)

    io = open(newnm, "w")
    write(io, out)
    close(io)
end

"""

Parses a markdown file into a ipynb file

Tries to handle graphics, but isn't perfect:

* Winston graphics work as expected

* PyPlot graphics have idiosyncracies:

- basic usage requires a call of `gcf()` as last entry of  a cell. This will *also* call clear on the figure, so that
any subsequent figures are added to a new canvas

- for 3d usage, this is not the case. The 3d graphics use a different backend and the display is different.

* `Gadfly` graphics are not (yet) supported, though this should be addressed .

"""
function mdToPynb(fname::AbstractString)

    m = make_module()

    newblocks = Any[]
    added_gadfly_preamble = false

    #    process_block("using WeavePynb, LaTeXStrings, Plots; pyplot()", m)
    process_block("using WeavePynb, LaTeXStrings, Plots; gr()", m)
    safeeval(m, Meta.parse("macro q_str(x)  \"`\$x`\" end"))

    out = Markdown.parse_file(fname,  flavor=Markdown.julia)
    for i in 1:length(out.content)
        cell = Dict()
        cell["metadata"] = Dict()
##        cell["prompt_number"] = i

        if isa(out.content[i], Markdown.Code)
            println("==== Block Code ====")
            println(out.content[i])
            ## Code Blocks are evaluated and their last value is added to the output
            ## this is different from IJulia, but similar.

            ## For graphics, we *now* use `GR`. This requires a flagging
            ## of `figure` when opening a code block. Otherwise, we get
            ## world age issues, as of v0.6.0

            txt = out.content[i].code
            lang = out.content[i].language

            ## we need to set
            ## nocode, noeval, noout
            langs = map(lstrip, split(lang, ","))

            docode, doeval, doout = true, true, true
            if "nocode" in langs
                docode = false
            end
            if "verbatim" in langs || "noeval" in langs
                doeval, doout = false, false
            end
            if "noout" in langs
                doout = false
            end

            ## language is used to pass in arguments
            result = nothing
            if doeval
                    result = process_block(txt, m)
            end

            !docode && (txt = "")

            cell["cell_type"] = "code"
            cell["execution_count"] = 1
            #cell["collapsed"] = false
            #cell["language"] = "python"
            #            cell["input"] = txt
            cell["source"] = [txt]

            println("Type of result: $(typeof(result))")

            if result == nothing
                cell["outputs"] = []
            elseif isa(result, Question)
                println("Process a question...")
                if isa(result, WeavePynb.Radioq)
                    # shove in an empty cell
                    #                cell["input"] = ""
#                    cell["outputs"] = []
                    tmp = Mustache.render(md_templates["Radioq"], items=["❆ $i" for i in result.labels])
                    cell["cell_type"] = "markdown"
                    cell["source"] = tmp
                    delete!(cell, "execution_count")
                else
                    continue
                end

            elseif isa(result, Verbatim)
                "Do not execute input, show as is"
                #                cell["input"] = result.x
                #                cell["outputs"] = []
                cell["cell_type"] = "markdown"
                cell["source"] = "<pre>$(result.x)</pre>"
                delete!(cell, "execution_count")
                ## We should be able to do this, but instead we now get World Age issues
            elseif isa(result, Plots.Plot)
#                if false  # XXX this has issues?
                    tmp = tempname()*".png"
                    Base.invokelatest(png, result, tmp)

                    dpi = 120
                    cell["outputs"] = [Dict(
                                            "output_type" => "execute_result",
                                            "execution_count" => 1,
                                            "data" => Dict("text/plain" => "Plot(...)",
                                                           "image/png" => base64encode(read(tmp, String))
                                                           #"image/png" => out
                                                           ),
                                       "metadata" => Dict("image/png" => Dict("width"=>5*dpi, "height"=>4*dpi))
                    )]
#                else
#                    cell["outputs"] = []
#                end

                # elseif string(typeof(result)) == "FramedPlot"
            #     ## Winston graphics
            #     cell["outputs"] = [render_winston(result)]
            # elseif  isa(result, Plots.Plot)
            #     # if isa(result, Plots.Plot{Plots.GadflyBackend})
            #     #     ## Gadfly graphics
            #     #     if !added_gadfly_preamble
            #     #         ## XXX this is *not* working, needed to figure out preamble... XXX
            #     #         ## Seems like injecting <script> failes.
            #     #         const gadfly_preamble = joinpath(dirname(@__FILE__), "..", "tpl", "gadfly-preamble.js")
                #     #         script = "<script>$(readstring(gadfly_preamble))</script>"
            #     #         added_gadfly_preamble = true
            #     #     end
            #     #     cell["outputs"] = [render_gadfly(result)]
            #     # else
            #         if  isa(result, Plots.Plot{Plots.PlotlyBackend})
            #         cell["outputs"] = [render_plotly(result)]
            #     end
            # elseif string(typeof(result)) == "Figure"
            #     ## *basic* PyPlot graphics.
            #     "Must do gcf() for last line"
            #     cell["outputs"] = [render_pyplot(result)]
            #     if ismatch(r"gcf\(\)$", txt)
            #         cell["input"] = join(split(txt, "\n")[1:(end-1)], "\n") ## trim last line which is gcf()
            #     else
            #         cell["input"] = txt ## trim last line which is gcf()
            #     end
            else
                ## Catch all
                tmp = Dict()
                tmp["metadata"] =Dict()
                mtype =  bestmime(result)
                tmp["output_type"] = "execute_result"
                tmp["execution_count"] = 1
                #                outtype = ifelse(ismatch(r"latex", string(mtype)), "latex", "text")
                outtype = ifelse(occursin(r"latex", string(mtype)), "text/latex", "text/plain")
                    output = ""
                try
                    output =  [sprint(io -> Base.invokelatest(show, io, mtype, result))]
                catch e
                end
                tmp["data"] = Dict()
                tmp["data"][outtype] = collect(output)

                cell["outputs"] = [tmp]
            end

        else
            cell["cell_type"] = "markdown"
            BigHeader = Union{Markdown.Header{1},Markdown.Header{2}}
            if isa(out.content[i], Markdown.Header)
                d = Dict()
                d["internals"] = Dict()
                if isa(out.content[i], BigHeader)
                    d["internals"]["slide_helper"] = "subslide_end"
                end
                d["internals"]["slide_type"] = "subslide"
                d["slide_helper"]="slide_end"
                d["slideshow"] = Dict()
                d["slideshow"]["slide_type"] = isa(out.content[i], BigHeader) ? "slide" : "subslide"
                cell["metadata"] = d
            end

            result = out.content[i]

            cell["source"] = sprint(io -> Markdown.html(io, out.content[i]))
        end

println("xxxxx")
cell["source"] == String[""] && println("XXXXXXX")
#        println("Source is ", cell["source"])

        push!(newblocks, JSON.json(cell))
    end
println("Here it is.....")
@show ipynb_tpl
    ## return string
    Mustache.render(ipynb_tpl, Dict("TITLE" => "TITLE", "CELLS" => join(newblocks, ",\n")))

end


text:
using Millboard
using Test

Millboard.set_table_mode(:grid)
Base.eval(:(have_color = false))

board = Dict()
@test table(board) |> string == ""

board = Dict(:a=>:b, "a"=>"b")
@test table(board) |> string == """
+-----+-------+
| KEY | VALUE |
+=====+=======+
|  :a |    :b |
+-----+-------+
|   a |     b |
+-----+-------+"""

board = Dict("1x3"=>[1 2 3], "2x3"=>[1 2 3; 4 5 6], "3x1"=> [1; 2; 3])
@test table(board) |> string == """
+-----+-------+
| KEY | VALUE |
+=====+=======+
| 1x3 | 1 2 3 |
+-----+-------+
|     | 1 2 3 |
| 2x3 | 4 5 6 |
+-----+-------+
|     |     1 |
|     |     2 |
| 3x1 |     3 |
+-----+-------+"""

Base.eval(:(have_color = true))


text:
module FundamentalStats

# Write your package code here.
    using Statistics, FFTW, Distributions
    # about poisson process
    export RCDF
    export ECDF
    export PowerSpectrum
    export AutoCorrelation
    export Variogram
    export Binning
    export Smooth
    export Q_Quantile
    export GaussFilter
    export StepFilter

    function StepFilter(t::Vector=[nothing],s = 1.;
                            data::Vector = t, stepsize = s)

        x2 = [i for i in 0:(1/stepsize):t[end]]
        y2 = [sum( (data.<i+stepsize/2) .& (data.>i-stepsize/2) )/stepsize for i in x2]
        return x2,y2
    end

    function Gaussian(t::Vector=[nothing],σ = 1.;
                            data::Vector = t, sigma = σ)::Vector{Float64}
        return (1/(sqrt(2*π)*sigma) ) .* exp.(-(data.^2)./(2 *sigma^2))
    end

    function GaussFilter(t::Vector = [nothing],  σ = 1.;
                                    data::Vector = t, sigma = σ)
        result = zeros(length(data))
        for i in 1:1:length(data)
            result[i] = sum(Gaussian(data = data.-data[i],σ))
        end
        return result
    end

    function Q_Quantile(x = nothing,y = Normal(0,1), num = 100; data_x = x, data_y = y)::Taple{Vector{Float64},Vector{Float64}}
        t = collect(0:(1/num):1)
        Q1 = quantile(data_x,t)
        Q2 = quantile(data_y,t)
        return Q1, Q2
    end

    function MakeWindowFilter(t::Vector,start = 0., finish = nothing,dt = 1., filter=Gaussian, params::Vector = [1.])
        if finish == nothing
            finish = t[end]
        end
        T::Vector{Float64} = [i for i in start:dt:finish]
        N::Vector{Float64} = [ sum(filter(t.-i, params...)) for i in T]
        return T, N
    end

    function Binning(X::Vector=[nothing], Y::Vector=[nothing], N = sqrt(length(X)) ;
         x::Vector=X, y::Vector=Y, num = N )::Tuple{Vector{Float64},Vector{Float64},Vector{Float64} }
        p::Vector{Int64} = sortperm(x)
        x_::Vector{Float64} = x[p]
        y_::Vector{Float64} = y[p]
        step::Float64 = (x_[end] - x_[1])/(num)
        start::Float64 = x_[1] + step
        x_save::Vector{Float64} = zeros(Float64,num+1)
        y_save::Vector{Float64} = zeros(Float64,num+1)
        y_err::Vector{Float64}  = zeros(Float64,num+1)
        P::Vector{Bool} = zeros(Bool,num+1)
        n::Int64 = 1
        Y::Float64 = 0
        Y2::Float64 = 0
        N::Int64 = 0
        for i in 1:length(x)
            while x_[i] > start
                x_save[n] = start - step/2
                y_save[n] = Y/N
                y_err[n] = sqrt( Y2/(N)-(Y/N)^2 )/sqrt(N)
                P[n] = (N!=0)
                Y = 0
                Y2 = 0
                start += step
                n += 1
                N=0
            end
            Y2 += y_[i]^2
            Y += y_[i]
            N += 1
        end
        x_save[n] = start - step/2
        y_save[n] = Y/N
        y_err[n] = sqrt( Y2/(N)-(Y/N)^2 )/sqrt(N)
        P[n] = (N!=0)
        return x_save[P], y_save[P],y_err[P]
    end

    function PowerSpectrum(x::Vector = [nothing], f::Bool = false;
                                data::Vector = x, full::Bool = f)::Vector{Float64}
        X = zeros(length(data)*2-1)
        X[1:length(data)] .= (data .- mean(data))
        X = abs2.(fft(X))
        if full
            return X
        else
            return X[1:length(data)]
        end
    end

    function AutoCorrelation(x::Vector = [nothing], l::Int64=10;
                                data::Vector = x,lag::Int64 = l)::Tuple{Vector{Float64}, Vector{Float64}}
        X = PowerSpectrum(data = data, full = true)
        X = real.(ifft(X))[1:lag]
        X = X./X[1]
        return X
    end

    function Variogram(x::Vector = [nothing], l::Int64=10;
                                data::Vector = x,lag::Int64 = l)::Tuple{Vector{Float64}, Vector{Float64}}
        n::Float64 = length(data)
        X::Vector{Float64} = collect(range(1,lag, length=lag))
        Y::Vector{Float64} = [sum( (data[1:end-i]-data[i+1:end]).^2 )/(n-i) for i in 1:l]
        return X, Y
    end

    # basic statistical tools
    ```
    For analyses on the stochastic process, explanatory data analysis do have crucial role. We thus provide some statistical tools, such as estimation of empirical cumulative distributions and ranked cumulative distributions.
    ```

    function RCDF(x::Vector = [nothing],b::Bool = true;
                            data::Vector = x, complementary::Bool = b)::Tuple{Vector, Vector{Float64}}
        if complementary
            return sort(data), [i for i in 1:1:length(data)]
        else
            return sort(data), [i for i in length(data):-1:1]
        end
    end

    function ECDF(x::Vector = [nothing], b::Bool = true;
                            data::Vector = x, complementary::Bool = b)::Tuple{Vector{Float64}, Vector{Float64}}
        n::Int64 = length(data)
        a::Vector = diff(sort(data)) .!= 0
        y::Vector = [i for i in 1:length(data)-1]
        num::Vector = y[a]
        if complementary
            pushfirst!(num,0.0)
            num = n.-num
            return num./n,sort(unique(data))
        else
            push!(num,length(data))
            return num./n,sort(unique(data))
        end
    end

end

max_stars_repo_name: yukki-jpn/JuliaTSA

max_issues_repo_name: yukki-jpn/JuliaTSA

max_forks_repo_name: yukki-jpn/JuliaTSA


text:
# For the function that should be executed on each single core, use @everywhere macro
@everywhere function myf(x)
    return 4.0 / (1.0 + x * x)
end

function mypi_par(n)

a = 0.0
b = 1.0
dx = (b - a) / (n - 1)

# cur_range = range(a, b, step = dx)
# tmpres = pmap(myf, cur_range)
# s = sum(tmpres)

s = 0.0
s = @distributed (+) for i = 1:n
    x = (i - 1) * dx + a
    myf(x)
end

s = (s - 0.5 * (myf(a) + myf(b))) * dx
return s

end

N = 33554432
for i = 1 : 50
    @time tmp = mypi_par(N)
end


text:
using GaussianProcesses
using Random

Random.seed!(20140430)
# Training data
n=183;                          #number of training points
x = 2π * rand(n);              #predictors
y = sin.(x) + 0.05*randn(n);   #regressors

#Select mean and covariance function
mZero = MeanZero()                   #Zero mean function
kern = SE(0.0,0.0)                   #Sqaured exponential kernel (note that hyperparameters are on the log scale)
                    # log standard deviation of observation noise (this is optional)
gp = GP(x,y,mZero,kern)       #Fit the GP

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
"""
convert to AutomotiveDrivingModels.Scene
"""
function state2scene(mdp::CarMDP, s::CarMDPState, car_type::VehicleDef = mdp.car_type)
    scene = Scene()
    push!(scene, Vehicle(s.ego, mdp.ego_type, EGO_ID))
    push!(scene, Vehicle(s.car, mdp.car_type, EGO_ID+1))
    return scene
end

function animate_states(mdp::CarMDP, states::Vector{CarMDPState}, actions::Vector{CarMDPAction};
                        overlays=SceneOverlay[IDOverlay()],
                        cam=StaticCamera(VecE2(0, -5.), 17.))
    duration = length(states)*mdp.ΔT
    fps = Int(1/mdp.ΔT)
    function render_states(t, dt)
        frame_index = Int(floor(t/dt)) + 1
        scene = state2scene(mdp, states[frame_index])
        return AutoViz.render(scene,
                mdp.env,
              cat(1, overlays,   TextOverlay(text = ["v: $(get_ego(scene).state.v)"],
                                            font_size=20,
                                            pos=VecE2(mdp.env.params.x_min + 3.,6.),
                                            incameraframe=true),
                                 TextOverlay(text = ["Acc: $(actions[frame_index].acc)"],
                                            font_size=20,
                                            pos=VecE2(mdp.env.params.x_min + 3.,8.),
                                            incameraframe=true),
                                TextOverlay(text = ["step: $frame_index"],
                                            font_size=20,
                                            pos=VecE2(mdp.env.params.x_min + 3.,4.),
                                            incameraframe=true)),
                cam=cam,
                car_colors=get_colors(scene))
    end
    return duration, fps, render_states
end

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
using MarkElide
using Base.Test

using MarkElide

f(x) = @mark :a println("hi $x")

println("non elided call:")
f(2)

println("elided call:")
@elide :a f(3)


text:
# Model configuration and hyperparameters
using Flux

const config_large = Dict(                      # Description                 | Release param
    "model"                 => "xlnet-large",   #
    "ff_activation"         => gelu,            # hidden feed-forward act.    | ff_activation
    "attn_type"             => "bi",            #
    "d_head"                => 64,              #
    "d_inner"               => 4096,            # intermediate size           | d_inner
    "d_model"               => 1024,            # hidden size                 | d_model
    "hidden_dropout_prob"   => 0.1,             #
    "n_head"                => 16,              # number of attention-heads   | d_head
    "n_layer"               => 24,              # number of hidden layers     | n_head
    "summary_activation"    => "tanh",          #
    "summary_attn_dropout"  => 0.1,             #
    "summary_type"          => "last",          #
    "vocab_size"            => 32000,           # vocab_size                  | n_token
    "untie_r"               => true             # untie_r
)

const config_base = Dict(                       # Description                 | Release param
    "model"                 => "xlnet-base",    #
    "ff_activation"         => gelu,            # hidden feed-forward act.    | ff_activation
    "attn_type"             => "bi",            #
    "d_head"                => 64,              #
    "d_inner"               => 3072,            # intermediate size           | d_inner
    "d_model"               => 768,             # hidden size                 | d_model
    "hidden_dropout_prob"   => 0.1,             #
    "n_head"                => 12,              # number of attention-heads   | d_head
    "n_layer"               => 12,              # number of hidden layers     | n_head
    "summary_activation"    => "tanh",          #
    "summary_attn_dropout"  => 0.1,             #
    "summary_type"          => "last",          #
    "vocab_size"            => 32000,           # vocabulary_size             | n_token
    "untie_r"               => true             # untie_r
)

# WIP


text:
"""
Basic stabilizer codes
"""
module BasicModels

# imports
using ..Model
using ..PauliTools:to_bsf, to_pauli

# exports
export BasicCode, FiveQubitCode, SteaneCode

struct BasicCode <: StabilizerCode
    stabilizers::BitMatrix
    logical_xs::BitMatrix
    logical_zs::BitMatrix
    nkd::Tuple{Int,Int,Union{Int,Missing}}
    label::String
end
"""
    BasicCode <: StabilizerCode

    BasicCode(
        pauli_stabilizers::AbstractVector{<:AbstractString},
        pauli_logical_xs::AbstractVector{<:AbstractString},
        pauli_logical_zs::AbstractVector{<:AbstractString},
        nkd::Tuple{Integer,Integer,Union{Integer,Missing}}=nothing,
        label::AbstractString=nothing
    )

Construct a basic code from string representations of stabilizers and logical operators.

Paulis are expressed as strings of capitalized I, X, Y, Z characters, with one character per
physical qubit. Logical X and Z operators are in matching order, with one of each for each
logical qubit. Optional `nkd` defaults to `n` and `k` evaluated and `d` missing. Optional
`label` defaults to "Basic [n,k,d]".

# Examples
```jldoctest
julia> using Qecsim.BasicModels

julia> code = BasicCode(["ZZI", "IZZ"], ["XXX"], ["IIZ"])  # 3-qubit repetition
BasicCode(["ZZI", "IZZ"], ["XXX"], ["IIZ"], (3, 1, missing), "Basic [3,1,missing]")

julia> validate(code)  # no error indicates operators satisfy commutation relations

julia> nkd(code)  # default nkd
(3, 1, missing)

julia> label(code)  # default label
"Basic [3,1,missing]"
```
"""
function BasicCode(
    pauli_stabilizers::AbstractVector{<:AbstractString},
    pauli_logical_xs::AbstractVector{<:AbstractString},
    pauli_logical_zs::AbstractVector{<:AbstractString},
    nkd=nothing,
    label=nothing
)
    # derived defaults (convert derived nkd type for derived label)
    dnkd::Tuple{Int,Int,Union{Int,Missing}} = !isnothing(nkd) ? nkd : (
        length(pauli_stabilizers) > 0 ? length(pauli_stabilizers[1]) : 0,
        length(pauli_logical_xs),
        missing
    )
    dlabel = !isnothing(label) ? label : "Basic [$(dnkd[1]),$(dnkd[2]),$(dnkd[3])]"
    return BasicCode(to_bsf(pauli_stabilizers), to_bsf(pauli_logical_xs),
                     to_bsf(pauli_logical_zs), dnkd, dlabel)
end
Model.stabilizers(code::BasicCode) = code.stabilizers
Model.logical_xs(code::BasicCode) = code.logical_xs
Model.logical_zs(code::BasicCode) = code.logical_zs
Model.nkd(code::BasicCode) = code.nkd
Model.label(code::BasicCode) = code.label
function Base.show(io::IO, x::BasicCode)
    show(io, typeof(x)); print(io, '(')
    for f in (x.stabilizers, x.logical_xs, x.logical_zs)
        show(io, to_pauli(f)); print(io, ", ")
    end
    show(io, x.nkd); print(io, ", ")
    show(io, x.label); print(io, ')')
end

"""
    FiveQubitCode() -> BasicCode

Construct 5-qubit [5,1,3] code as a [`BasicCode`](@ref).
"""
function FiveQubitCode()
    return BasicCode(["XZZXI", "IXZZX", "XIXZZ", "ZXIXZ"], ["XXXXX"], ["ZZZZZ"], (5, 1, 3),
                      "5-qubit")
end

"""
    SteaneCode() -> BasicCode

Construct Steane [7,1,3] code as a [`BasicCode`](@ref).
"""
function SteaneCode()
    return BasicCode(["IIIXXXX", "IXXIIXX", "XIXIXIX", "IIIZZZZ", "IZZIIZZ", "ZIZIZIZ"],
                     ["XXXXXXX"], ["ZZZZZZZ"], (7, 1, 3), "Steane")
end

end


text:
###################################################
## network/solvers/heuristic.jl
## Heuristic from MIT Senseable City Lab
###################################################

"""
    heuristicTimes :
    optimize travel times using Senseable City Heuristic
"""
function heuristicTimes(s::IterativeState)
	et = getPathTimes(s.timings)
	g = s.data.network.graph
    paths = s.paths
    tripData = s.trips
    roads = s.data.network.roads

    # compute relative error
    relErr = 0.
    for trip in tripData
    	relErr += abs(trip.time - et[trip.orig, trip.dest])/trip.time
    end
    relErr = relErr/length(tripData)
    # compute offset for each road
    offset = Dict{Tuple{Int, Int}, Float64}()
    for (i,trip) in enumerate(tripData)
    	for edge = keys(paths[i][1])
    		key = (src(edge), dst(edge))
    		if haskey(offset, key)
    			offset[key] += (et[trip.orig, trip.dest] - trip.time) * trip.weight
    		else
    			offset[key] = (et[trip.orig, trip.dest] - trip.time) * trip.weight
    		end
    	end
    end
    # set initial value for k
    k = 1.2
    times = spzeros(nNodes(s.data.network), nNodes(s.data.network))
    # use offsets go update road times
    while true
    	# update roadTimes
    	for roadKey in keys(offset)
    		if offset[roadKey] < 0
	    		times[roadKey[1], roadKey[2]] = k * s.timings.times[roadKey[1], roadKey[2]]
	    	else
	    		times[roadKey[1], roadKey[2]] = s.timings.times[roadKey[1], roadKey[2]]/k
	    	end
    	end
    	# compute new error and see if it is better than the old one
    	newErr = 0.
    	for (i, trip) in enumerate(tripData)
    		estTime = 0.
    		for edge = keys(paths[i][1])
    			estTime += times[src(edge), dst(edge)]
    		end
    		newErr += abs(trip.time - estTime)/trip.time
    	end
        newErr = newErr/length(tripData)
    	if newErr < relErr
    		break
    	else
    		k = 1 + (k - 1) * 0.75
    		k < 1.0001 && break
    	end
    end
    if k < 1.0001
        println("Convergence obtained")
    end
    # compute times for roads that did not appear in any path
    println("We have $(length(roads)-length(offset)) roads that did not appear in any path.")
    unknownRoads = Dict{Tuple{Int,Int}, Int}()
    for (key, road) in roads
    	if !(haskey(offset, key))
            unknownRoads[key] = 0
        end
    end
    for key in keys(unknownRoads)
        for neighbor in RoutingNetworks.in_neighbors(g, key[1])
            if !((neighbor, key[1]) in keys(unknownRoads))
                unknownRoads[key] += 1
            end
        end
        for neighbor in RoutingNetworks.outneighbors(g, key[2])
            if !((key[2], neighbor) in keys(unknownRoads))
                unknownRoads[key] += 1
            end
        end
    end

    # find most connected road, update its time and repeat
    while length(unknownRoads) > 0
        maxValue = 0
        maxKey = 0
        for (key, value) in unknownRoads
            if value > maxValue
                maxValue = value
                maxKey = key
            end
    	end
        newRoad = maxKey
        neighborV = 0.
        count = 0
        for neighbor in RoutingNetworks.in_neighbors(g, newRoad[1])
            if (neighbor, newRoad[1]) in keys(unknownRoads)
                unknownRoads[(neighbor, newRoad[1])] += 1
            else
                neighborV += times[neighbor, newRoad[1]]/roads[(neighbor, newRoad[1])].distance
                count += 1
            end
        end
        for neighbor in RoutingNetworks.outneighbors(g, newRoad[2])
            if (newRoad[2], neighbor) in keys(unknownRoads)
                unknownRoads[(newRoad[2], neighbor)] += 1
            else
                neighborV += times[newRoad[2], neighbor]/roads[(newRoad[2], neighbor)].distance
                count += 1
            end
        end
        neighborV = neighborV/count
        times[newRoad[1], newRoad[2]] = roads[newRoad].distance/neighborV
        delete!(unknownRoads, newRoad)
    end
    return times
end


text:
@inline function velocity(x, v, pbest, gbest, parameters)
    r1 = parameters.C1 * rand()
    r2 = parameters.C2 * rand()
    parameters.ω * v + r1 * (pbest - x) + r2 * (gbest - x)
end