text:
using ProtoBuf

using BondLengthDistributions
using SmuUtilities
using dataset_pb2

export MatchingParameters
export SmuMolecule
export set_initial_score_and_incrementer!
export generate_search_state
export place_bonds!

"""Parameters controlling settable parameters during matching.
"""
mutable struct MatchingParameters
  # If true, all atoms must have their full valence satisfied.
  must_match_all_bonds::Bool
end
MatchingParameters() = MatchingParameters(true)

"""Description of a molecule used for building different bonding patterns.
"""
mutable struct SmuMolecule
  # The starting bond_topology, optionally with Hydrogens, and other singly
  # bonded atoms attached.
  starting_bond_topology::BondTopology

  # For each atom, the maximum number of bonds that can be attached.
  max_bonds::Vector{Int32}

  # During enumeration of possible molecules, we start with the
  # Hydrogen atom attached form of the BondTopology. That causes
  # atoms to have an initial number of bonds attached.
  bonds_with_hydrogens_attached::Vector{Int32}

  # During building a molecule, how many bonds are currently attached
  # to each atom
  current_bonds_attached::Vector{Int32}

  # Scores can be implemented as either multiplicative or additive forms.
  # Each needs an initial value and a function to accumulate the score.
  initial_score::Float32
  accumulate_score::Function

  # The `bonds` and `scores` are copied from the keys and values of
  # the bonds_to_scores dictionary given to the constructor.
  # The important thing is that these arrays are in the same order.
  bonds::Vector{Tuple{Int32, Int32}}
  scores::Vector{Vector{Float32}}

  # If true, then a successful match is only reported if each atom has the
  # maximum number of bonds attached.
  must_match_all_bonds::Bool

  function SmuMolecule(hydrogens_attached::BondTopology,
              bonds_to_scores::Dict{Tuple{Int32,Int32}, Vector{Float32}},
              matching_parameters::MatchingParameters)
    mol = new();
    mol.starting_bond_topology = hydrogens_attached;
    # Delegate processing to a separate function, keep this small.
    smu_molecule_constructor!(mol);

    # separately extract the bonds and scores
    mol.bonds = collect(keys(bonds_to_scores));
    mol.scores = collect(values(bonds_to_scores));
    mol.must_match_all_bonds = matching_parameters.must_match_all_bonds
    mol
  end
end

"""Called from the SmuMolecule inner constructor as a regular function.
Initializes many of the attributes of `mol`.
"""
function smu_molecule_constructor!(mol::SmuMolecule)
  mol.accumulate_score = Base.:+
  mol.initial_score = 0.0
  natoms = length(mol.starting_bond_topology.atoms);
  mol.bonds_with_hydrogens_attached = zeros(Int32, natoms);
  mol.max_bonds = zeros(Int32, natoms);
  for i in 1:natoms
    mol.max_bonds[i] = SmuUtilities.smu_atom_type_to_max_con(mol.starting_bond_topology.atoms[i])
  end;
  # And during molecule building the current state.
  mol.current_bonds_attached = zeros(Int32, natoms);

  hasproperty(mol.starting_bond_topology, :bonds) || return
  for bond in mol.starting_bond_topology.bonds
    mol.bonds_with_hydrogens_attached[bond.atom_a + 1] += 1
    mol.bonds_with_hydrogens_attached[bond.atom_b + 1] += 1
  end
end

function set_initial_score_and_incrementer!(score::T, op, smu_molecule::SmuMolecule) where{T<:Number}
  smu_molecule.initial_score = score
  smu_molecule.accumulate_score = op
end

function accumulate_score(mol::SmuMolecule, existing_score::T, increment::T)::Float32 where {T<:Number}
  return mol.accumulate_score(existing_score, increment)
end

"""For each pair of atoms, return a list of plausible bond types.

The result will be used to enumerate all the possible bonding forms.
The resulting Vector has the same length as mol.bonds, and each item
in that vector is avector of the indices of the plausible bond types
for that connection.
Args:
  mol:
Returns:
  Vector of Vectors - one for each atom pair.
"""
function generate_search_state(mol::SmuMolecule)::Vector{Vector{Int32}}
  result = Vector{Vector{Int32}}()
  for ndx in 1:length(mol.bonds)
    # For each pair of atoms, the plausible bond types - non zero score.
    push!(result, findall(x->x > 0.0, mol.scores[ndx]))
  end

  return result
end

"""Possibly add a new bond to the current config.

If the bond can be placed, updates self._current_bonds_attached for
both `a`` and `a2`.
  Args:
    a1:
    a2:
    btype:
  Returns:
    true always - even if no bond is added. Perhaps change to return nothing.
print(f"Trying to place bond {btype} current {self._current_bonds_attached[a1]} and {self._current_bonds_attached[a2]}")
"""
function _place_bond!(a1::Int32, a2::Int32, btype, mol::SmuMolecule)::Bool
  btype > 0 || return true
  mol.current_bonds_attached[a1] + btype > mol.max_bonds[a1] && return false
  mol.current_bonds_attached[a2] + btype > mol.max_bonds[a2] && return false

  mol.current_bonds_attached[a1] += btype
  mol.current_bonds_attached[a2] += btype
  return true
end

"""Place bonds corresponding to `state`.

Args:
  state: for each pair of atoms, the kind of bond to be placed.
Returns:
  If successful, the score.
"""
function place_bonds!(state::Tuple,
                      mol::SmuMolecule)::Union{dataset_pb2.BondTopology, Nothing}
  result = dataset_pb2.BondTopology()  # To be returned.
  copy!(result, mol.starting_bond_topology)  # Only Hydrogens attached.
  result.score = mol.initial_score

  # Initialize state in mol to Hydrogens attached.
  mol.current_bonds_attached = copy(mol.bonds_with_hydrogens_attached)

  for i in 1:length(state)
    a1 = mol.bonds[i][1]
    a2 = mol.bonds[i][2]
    btype = state[i]
    _place_bond!(a1, a2, btype - 1, mol) || return nothing
    # If the bond is anything other than BOND_UNDEFINED, add it to result.
    (btype - 1) == BondTopology_BondType.BOND_UNDEFINED || add_bond!(a1 - 1, a2 - 1, btype - 1, result)
    result.score = accumulate_score(mol, result.score, mol.scores[i][btype])
  end

  # Optionally check whether all bonds have been matched
  mol.must_match_all_bonds || return result

  @debug("Cf bonds $(mol.current_bonds_attached) $(mol.max_bonds)")
  mol.current_bonds_attached == mol.max_bonds ? result : nothing
end

max_stars_repo_name: IanAWatson/smu

max_issues_repo_name: IanAWatson/smu

max_forks_repo_name: IanAWatson/smu


text:
module essentials
"""
Computes the Long run average cost and
long run marginal cost

Lrac = Lrtc/x

Lrcm = d* Lrtc/d*x
"""

end  # module essentials

ext: jl

max_stars_repo_name: adamwillisXanax/SolvencyPredictor

max_issues_repo_name: adamwillisXanax/SolvencyPredictor

max_forks_repo_name: adamwillisXanax/SolvencyPredictor


text:
function update!{T}(tr::OptimizationTrace{T},
                 iteration::Integer,
                 f_x::Real,
                 grnorm::Real,
                 dt::Dict,
                 store_trace::Bool,
                 show_trace::Bool,
                 show_every::Int = 1,
                 callback = nothing)
    os = OptimizationState{T}(iteration, f_x, grnorm, dt)
    if store_trace
        push!(tr, os)
    end
    if show_trace
        if iteration % show_every == 0
            show(os)
        end
    end
    if callback != nothing && (iteration % show_every == 0)
        if store_trace
            callback(tr)
        else
            callback(os)
        end
    end
    return
end


text:
const PKPLOTSTYLE = (
(:solid, :blue, :circle, :blue),
(:solid, :red, :utriangle, :red),
(:solid, :green, :diamond, :green),
(:solid, :magenta, :pentagon, :magenta),
(:solid, :purple, :heptagon, :purple),
(:solid, :indigo, :octagon, :indigo),
(:solid, :gold, :star, :gold),
(:solid, :yellow, :rect, :yellow),
(:solid, :gray, :xcross, :gray),
(:solid, :cyan, :cross, :cyan),
(:dot, :blue, :utriangle, :blue),
(:dot, :red, :circle, :red),
(:dot, :green, :rect, :green),
(:dot, :yellow, :diamond, :yellow),
(:dot, :gray, :cross, :gray),
(:dot, :cyan, :xcross, :cyan),
(:dot, :gold, :pentagon, :gold),
(:dot, :magenta, :star, :magenta),
(:dot, :purple, :octagon, :purple),
(:dot, :indigo, :heptagon, :indigo),
)

function plotstyle(n)
    if n >= 1 && n <= 20
        return PKPLOTSTYLE[n]
    elseif n > 1900
        return (:auto, :auto, :auto, :auto)
    end
    n -= 20
    linestyle   = (:dash, :dashdot, :dashdotdot)
    linecolor   = (:yellow, :gray, :cyan, :gold, :magenta, :purple, :indigo)
    markershape = (:star4, :star6, :star7, :star8, :rect, :star5, :diamond, :hexagon, :cross, :xcross, :utriangle, :dtriangle, :rtriangle, :ltriangle, :pentagon, :heptagon, :octagon, :vline, :hline)
    markercolor = (:gray, :gold, :magenta, :purple, :indigo)
    a = n % 3 + 1
    b = n % 7 + 1
    c = n % 19 + 1
    d = n % 5 + 1
    (linestyle[a], linecolor[b], markershape[c], markercolor[d])
end

@userplot PKPlot
@userplot PKElimpPlot

function luceil(x)
    fl = Int(floor(log10(x)))
    if fl < 0 fl = 0 end
    ceil(x/10^fl)*10^fl
end

@recipe function f(subj::PKPlot; lcd = :auto, tcd = :auto)
    x, y = subj.args

    if isa(lcd, Real)
        lc = luceil(maximum(x->isnan(x) ? -Inf : x, y)/lcd)
        yt = 0:lc:lc*lcd
    elseif isa(lcd, StepRange)
        yt = lcd
    elseif lcd == :all
        yt = y
    end

    if isa(tcd, Real)
        tc = luceil(maximum(x)/tcd)
        xt = 0:tc:tc*tcd
    elseif isa(tcd, StepRange)
        xt = tcd
    elseif tcd == :all
        xt = x
    end

    widen             --> true
    seriestype        --> :line
    xguide            --> "Time"
    link              --> :both
    legend            --> true
    grid              --> true
    gridstyle         --> :auto
    #ticks       := [nothing :auto nothing]
    #xlims             --> (minimum(x), maximum(x)*1.1)
    #ylims             --> (0, maximum(y)*1.1)

    if !isa(lcd, Symbol) || lcd != :auto
        yticks            --> yt
    end
    if !isa(tcd, Symbol) || tcd != :auto
        xticks            --> xt
    end

    seriescolor       --> :blue
    markershape       --> :circle
    markersize        --> 3
    markercolor       --> :match
    markerstrokealpha --> 0
    (x, y)
end

@recipe function f(subj::PKElimpPlot)
    x, y = subj.args
    seriestype        --> :line
    legend            --> false
    markersize        --> 0
    markerstrokealpha --> 0
    (x, y)
end

# Text label from ID
function plotlabel(d)
    title = ""
    if isnothing(d) return title end
    if length(d) > 0
        for (k, v) in d
            title *= "$(k) = $(v); "
        end
    end
    return title
end

"""
    pkplot(subj; ls = false, elim = false, xticksn = :auto, yticksn = :auto, kwargs...)

Plot for subject

* `ls` - concentration in log scale;
* `elim` - draw elimination curve;
* `xticksn` - number of ticks on x axis;
* `yticksn` - number of ticks on y axis/

"""
function pkplot(subj::AbstractSubject; ls = false, elim = false, xticksn = :auto, yticksn = :auto, kwargs...)
    time = subj.time
    obs  = subj.obs
    kwargs = Dict{Symbol, Any}(kwargs)
    k = keys(kwargs)

    if !(:plotstyle in k)
        kwargs[:linestyle], kwargs[:linecolor], kwargs[:markershape],  kwargs[:markercolor]  = PKPLOTSTYLE[1]
    else
        kwargs[:linestyle], kwargs[:linecolor], kwargs[:markershape],  kwargs[:markercolor]  = kwargs[:plotstyle]
    end
    if !(:title in k)
        kwargs[:title] = plotlabel(subj.id)
    end
    if !(:xlims in k)
        kwargs[:xlims] = (minimum(subj.time), maximum(subj.time)*1.1)
    end
    if !(:ylims in k)
        kwargs[:ylims] = (minimum(subj.obs), maximum(subj.obs)*1.15)
    end
    if !(:legend in k)
        kwargs[:legend] = true
    end
    if !(:ylabel in k)
        kwargs[:ylabel] = "Concentration"
    end
    if :yscale in k
        if kwargs[:yscale] in [:ln, :log, :log2, :log10] ls = true end
    end
    if ls == true
        inds = findall(x->x>0, subj.obs)
        time = subj.time[inds]
        obs = log.(subj.obs[inds])
        if (:ylims in k)
            kwargs[:ylims] = (0, log(kwargs[:ylims][2]))
        end
    end

    p = pkplot(time, obs;  lcd = yticksn, tcd = xticksn, kwargs...)
    if elim
        if length(subj.keldata) > 0
            rsq, rsqn = findmax(subj.keldata.ar)
            lz        = subj.keldata.a[rsqn]
            lzint     = subj.keldata.b[rsqn]
            ts        = subj.keldata.s[rsqn]
            te        = subj.keldata.e[rsqn]

            if ls true
                x = [ts, te]
                y = [lzint + lz * x[1], lzint + lz * x[2]]
            else
                x = collect(ts:(te-ts)/100:te)
                y = exp.(lzint .+ lz .* x)
            end
            pkelimpplot!(p, x, y)
        end
    end
    return p
end

function pkplot!(subj; ls = false, elim = false, xticksn = :auto, yticksn = :auto, kwargs...)
    time = subj.time
    obs = subj.obs
    kwargs = Dict{Symbol, Any}(kwargs)
    k = keys(kwargs)
    if !(:plotstyle in k)
        kwargs[:linestyle], kwargs[:linecolor], kwargs[:markershape],  kwargs[:markercolor]  = PKPLOTSTYLE[1]
    else
        kwargs[:linestyle], kwargs[:linecolor], kwargs[:markershape],  kwargs[:markercolor]  = kwargs[:plotstyle]
    end
    if !(:xlims in k)
        kwargs[:xlims] = (minimum(subj.time), maximum(subj.time)*1.1)
    end
    if !(:ylims in k)
        kwargs[:ylims] = (minimum(subj.obs), maximum(subj.obs)*1.15)
    end
    if !(:legend in k)
        kwargs[:legend] = true
    end
    if :yscale in k
        if kwargs[:yscale] in [:ln, :log, :log2, :log10] ls = true end
    end
    if ls == true
        inds = findall(x->x>0, subj.obs)
        time = subj.time[inds]
        obs = log.(subj.obs[inds])
        if (:ylims in k)
            kwargs[:ylims] = (0,log(kwargs[:ylims][2]))
        end
    end

    p = pkplot!(time, obs;  lcd = yticksn, tcd = xticksn, kwargs...)
    return p
end

function pageplot(data, id, ulist; kwargs...)
    kwargs = Dict{Symbol, Any}(kwargs)
    k = keys(kwargs)
    if !(:title in k)
        kwargs[:title] = plotlabel(id)
    end
    #=
    if !(:xticksn in k)
        kwargs[:xticksn] = 6
    end
    if !(:yticksn in k)
        kwargs[:yticksn] = 5
    end
    =#
    #utypes    = keys(styledict)
    fst       = true
    p         = nothing
    labvec    = Vector{Int}(undef, 0)
    # Make subdata by ID
    isnothing(id) ? subdata   = data : subdata   = subset(data, id)
    # Y lims
    if !(:ylims in k)
        kwargs[:ylims] = (findmin(x->minconc(x), getdata(subdata))[1], findmax(x->maxconc(x), getdata(subdata))[1]*1.15)
    end
    # Plotting subdata
    for subj in subdata
        #if isnothing(id) || id ⊆ subj.id
            if !isnothing(ulist)
                num = findfirst(x-> x ⊆ subj.id, ulist)
                style = plotstyle(num)
                if num ∈ labvec
                    kwargs[:label] = nothing
                else
                    kwargs[:label] = plotlabel(ulist[num])
                    push!(labvec, num)
                end
            else
                style = plotstyle(1)
            end
            if fst
                p = pkplot(subj; plotstyle = style, kwargs...)
                fst = false
            else
                pkplot!(subj; plotstyle = style, kwargs...)
            end

        #end
    end
    p
end

"""
    pkplot(data::DataSet{T};
    typesort::Union{Nothing, Symbol, AbstractVector{Symbol}} = nothing,
    pagesort::Union{Nothing, Symbol, AbstractVector{Symbol}} = nothing,
    sort::Union{Nothing, Dict{Symbol}} = nothing,
    kwargs...) where T <: AbstractSubject

PK plot for subject set.

* `typesort` - sort on page by this id key;
* `pagesort` - different pages by this id key;
* `sort` - use only subjects if sort ⊆ subject id.
"""
function pkplot(data::DataSet{T};
    typesort::Union{Nothing, Symbol, AbstractVector{Symbol}} = nothing,
    pagesort::Union{Nothing, Symbol, AbstractVector{Symbol}} = nothing,
    sort::Union{Nothing, Dict{Symbol}} = nothing,
    kwargs...) where T <: AbstractSubject

    kwargs = Dict{Symbol, Any}(kwargs)
    k = keys(kwargs)
    if !(:ls in k)
        kwargs[:ls] = false
    end
    if !(:elim in k)
        kwargs[:elim] = false
    end
    #=
    if !(:xticksn in k)
        kwargs[:xticksn] = :auto
    end
    if !(:yticksn in k)
        kwargs[:yticksn] = :auto
    end
    =#
    if !isnothing(sort) data = subset(data, sort) end

    if isnothing(typesort) && isnothing(pagesort)
        p = []
        printtitle = false
        if !(:title in k)
            printtitle = true
        end
        for subj in data
            if printtitle
                kwargs[:title] = plotlabel(subj.id)
            end
            if !(:legend in k)
                kwargs[:legend] = false
            end
            push!(p, pkplot(subj;  kwargs...))
        end
        return p
    end

    if !isnothing(typesort)
        if isa(typesort, Symbol) typesort = [typesort] end
        typelist = uniqueidlist(data, typesort)
    else
        typelist = nothing
        if !(:legend in k)
            kwargs[:legend] = false
        end
    end

    if !isnothing(pagesort)
        kwargs[:elim] = false
        if isa(pagesort, Symbol) pagesort = [pagesort] end
        p = []
        pagelist = uniqueidlist(data, pagesort)
        for id in pagelist
            push!(p, pageplot(data, id, typelist; kwargs...))
        end
        return p
    else
        if !(:title in k) && !isnothing(sort)
            kwargs[:title] = plotlabel(sort)
        end
        return pageplot(data, pagesort, typelist; kwargs...)
    end
end

function uniqueidlist(data::DataSet{T}, ::Nothing) where T <: AbstractIdData
    nothing
end


text:
# Generic concepts and functions

# a premetric is a function d that satisfies:
#
#   d(x, y) >= 0
#   d(x, x) = 0
#
abstract PreMetric

# a semimetric is a function d that satisfies:
#
#   d(x, y) >= 0
#   d(x, x) = 0
#   d(x, y) = d(y, x)
#
abstract SemiMetric <: PreMetric

# a metric is a semimetric that satisfies triangle inequality:
#
#   d(x, y) + d(y, z) >= d(x, z)
#
abstract Metric <: SemiMetric

# Generic functions

result_type(::PreMetric, ::AbstractArray, ::AbstractArray) = Float64

# Generic column-wise evaluation

function colwise!(r::AbstractArray, metric::PreMetric, a::AbstractVector, b::AbstractMatrix)
    n = size(b, 2)
    length(r) == n || throw(DimensionMismatch("Incorrect size of r."))
    for j = 1 : n
        @inbounds r[j] = evaluate(metric, a, view(b, :, j))
    end
    r
end

function colwise!(r::AbstractArray, metric::PreMetric, a::AbstractMatrix, b::AbstractVector)
    n = size(a, 2)
    length(r) == n || throw(DimensionMismatch("Incorrect size of r."))
    for j = 1 : n
        @inbounds r[j] = evaluate(metric, view(a, :, j), b)
    end
    r
end

function colwise!(r::AbstractArray, metric::PreMetric, a::AbstractMatrix, b::AbstractMatrix)
    n = get_common_ncols(a, b)
    length(r) == n || throw(DimensionMismatch("Incorrect size of r."))
    for j = 1 : n
        @inbounds r[j] = evaluate(metric, view(a, :, j), view(b, :, j))
    end
    r
end

function colwise!(r::AbstractArray, metric::SemiMetric, a::AbstractMatrix, b::AbstractVector)
    colwise!(r, metric, b, a)
end

function colwise(metric::PreMetric, a::AbstractMatrix, b::AbstractMatrix)
    n = get_common_ncols(a, b)
    r = Array(result_type(metric, a, b), n)
    colwise!(r, metric, a, b)
end

function colwise(metric::PreMetric, a::AbstractVector, b::AbstractMatrix)
    n = size(b, 2)
    r = Array(result_type(metric, a, b), n)
    colwise!(r, metric, a, b)
end

function colwise(metric::PreMetric, a::AbstractMatrix, b::AbstractVector)
    n = size(a, 2)
    r = Array(result_type(metric, a, b), n)
    colwise!(r, metric, a, b)
end

# Generic pairwise evaluation

function pairwise!(r::AbstractMatrix, metric::PreMetric, a::AbstractMatrix, b::AbstractMatrix)
    na = size(a, 2)
    nb = size(b, 2)
    size(r) == (na, nb) || throw(DimensionMismatch("Incorrect size of r."))
    for j = 1 : size(b, 2)
        bj = view(b,:,j)
        for i = 1 : size(a, 2)
            @inbounds r[i,j] = evaluate(metric, view(a,:,i), bj)
        end
    end
    r
end

function pairwise!(r::AbstractMatrix, metric::PreMetric, a::AbstractMatrix)
    pairwise!(r, metric, a, a)
end

function pairwise!(r::AbstractMatrix, metric::SemiMetric, a::AbstractMatrix)
    n = size(a, 2)
    size(r) == (n, n) || throw(DimensionMismatch("Incorrect size of r."))
    for j = 1 : n
        aj = view(a,:,j)
        for i = j+1 : n
            @inbounds r[i,j] = evaluate(metric, view(a,:,i), aj)
        end
        @inbounds r[j,j] = 0
        for i = 1 : j-1
            @inbounds r[i,j] = r[j,i]   # leveraging the symmetry of SemiMetric
        end
    end
    r
end

function pairwise(metric::PreMetric, a::AbstractMatrix, b::AbstractMatrix)
    m = size(a, 2)
    n = size(b, 2)
    r = Array(result_type(metric, a, b), (m, n))
    pairwise!(r, metric, a, b)
end

function pairwise(metric::PreMetric, a::AbstractMatrix)
    n = size(a, 2)
    r = Array(result_type(metric, a, a), (n, n))
    pairwise!(r, metric, a)
end

function pairwise(metric::SemiMetric, a::AbstractMatrix)
    n = size(a, 2)
    r = Array(result_type(metric, a, a), (n, n))
    pairwise!(r, metric, a)
end


text:
# ---
# jupyter:
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# +
# IJulia rocks! So does Plotly. Check it out

using Plotly
api_key = "" # visit https://plot.ly/api to generate an API username and password
username = ""

Plotly.signin(username, api_key)

# +
# Following data taken from http://julialang.org/ frontpage
benchmarks = ["fib", "parse_int", "quicksort3", "mandel", "pi_sum", "rand_mat_stat", "rand_mat_mul"]
platforms = ["Fortran", "Julia", "Python", "R", "Matlab", "Mathematica", "Javascript", "Go"]

data = {
        platforms[1] => [0.26, 5.03, 1.11, 0.86, 0.80, 0.64, 0.96],
        platforms[2] => [0.91, 1.60, 1.14, 0.85, 1.00, 1.66, 1.01],
        platforms[3] => [30.37, 13.95, 31.98, 14.19, 16.33, 13.52, 3.41 ],
        platforms[4] => [411.36, 59.40, 524.29, 106.97, 15.42, 10.84, 3.98  ],
        platforms[5] => [1992.00, 1463.16, 101.84, 64.58, 1.29, 6.61, 1.10   ],
        platforms[6] => [64.46, 29.54, 35.74, 6.07, 1.32, 4.52, 1.16 ],
        platforms[7] => [2.18, 2.43, 3.51, 3.49, 0.84, 3.28, 14.60],
        platforms[8] => [1.03, 4.79, 1.25, 2.36, 1.41, 8.12, 8.51]
        }

pdata = [ {"x"=>benchmarks,"y"=>data[k],"bardir"=>"h","type"=>"bar","name"=>k} for k = platforms ]

layout = {
          "title"=> "Julia benchmark comparison (smaller is better, C performance = 1.0)",
          "barmode"=> "group",
          "autosize"=> false,
          "width"=> 900,
          "height"=> 900,
          "titlefont"=>
          {
           "family"=> "Open Sans",
           "size"=> 18,
           "color"=> "rgb(84, 39, 143)"
           },
          "margin"=> {"l"=>160, "pad"=>0},
          "xaxis"=> {
                     "title"=> "Benchmark log-time",
                     "type"=> "log"
                     },
          "yaxis"=> {"title"=> "Benchmark Name"}
          }

response = Plotly.plot(pdata,["layout"=>layout])

# Embed in an iframe within IJulia
s = string("")
display("text/html", s)

# +
# checkout https://plot.ly/api/ for more Julia examples!
# But to show off some other Plotly features:
x = 1:1500
y1 = sin(2*pi*x/1500.) + rand(1500)-0.5
y2 = sin(2*pi*x/1500.)

fish = {"x"=>x,"y"=> y1,
	"type"=>"scatter","mode"=>"markers",
	"marker"=>{"color"=>"rgb(0, 0, 255)","opacity"=>0.5 } }

fit = {"x"=> x,"y"=> y2,
	"type"=>"scatter", "mode"=>"markers", "opacity"=>0.8,
	"marker"=>{"color"=>"rgb(255, 0, 0)"} }

layout = {"autosize"=> false,
    "width"=> 650, "height"=> 550,
    "title"=>"Fish School",
	"xaxis"=>{ "ticks"=> "",
        "gridcolor"=> "white",
        "zerolinecolor"=> "white",
        "linecolor"=> "white",
        "autorange"=> false,
        "range"=>[0,1500] },
	"yaxis"=>{ "ticks"=> "",
        "gridcolor"=> "white",
        "zerolinecolor"=> "white",
		"linecolor"=> "white",
        "autorange"=> false,
        "range"=>[-2.2,2.2] },
	"plot_bgcolor"=> "rgb(245,245,247)",
    "showlegend"=> false,
    "hovermode"=> "closest"}

response = Plotly.plot([fish, fit],["layout"=>layout])
s = string("")
display("text/html", s)

max_stars_repo_name: hwasiti/jupytext

max_issues_repo_name: royalosyin/jupytext

max_forks_repo_name: royalosyin/jupytext


text:
function hshoe_kernel(
	cx,
	cy,
	cz,
	x1,
	x2,
	x3,
	p1x,
	p1y,
	p1z,
	p2x,
	p2y,
	p2z;
	M∞::Real = 0.0
)

	β = sqrt(1.0 - M∞ ^ 2)

	ax = cx .- p1x'
	ay = (cy .- p1y') .* β
	az = (cz .- p1z') .* β

	bx = cx .- p2x'
	by = (cy .- p2y') .* β
	bz = (cz .- p2z') .* β

	na = @. sqrt(ax ^ 2 + ay ^ 2 + az ^ 2) + 1e-10
	nb = @. sqrt(bx ^ 2 + by ^ 2 + bz ^ 2) + 1e-10

	ab = @. ax * bx + ay * by + az * bz

	mult = @. (1.0 / na + 1.0 / nb) / (ab + na * nb)

	d2 = @. na - ax * x1 - ay * x2 - az * x3
	d3 = @. nb - bx * x1 - by * x2 - bz * x3

	mult2 = @. 1.0 / (na * d2)
	mult3 = @. 1.0 / (nb * d3)

	vx = @. (
		(ay * bz - az * by) * mult +
		(ay * x3 - az * x2) * mult2 -
		(by * x3 - bz * x2) * mult3
	) / (4 * π * β ^ 2)
	vy = @. (
		(az * bx - ax * bz) * mult +
		(az * x1 - ax * x3) * mult2 -
		(bz * x1 - bx * x3) * mult3
	) / (4 * π * β)
	vz = @. (
		(ax * by - ay * bx) * mult +
		(ax * x2 - ay * x1) * mult2 -
		(bx * x2 - by * x1) * mult3
	) / (4 * π * β)

	[vx, vy, vz]

end

function vlm_kernel(
	cx,
	cy,
	cz,
	p1x,
	p1y,
	p1z,
	p2x,
	p2y,
	p2z;
	M∞::Real = 0.0
)

	β = sqrt(1.0 - M∞ ^ 2)

	ax = cx .- p1x'
	ay = (cy .- p1y') .* β
	az = (cz .- p1z') .* β

	bx = cx .- p2x'
	by = (cy .- p2y') .* β
	bz = (cz .- p2z') .* β

	na = @. sqrt(ax ^ 2 + ay ^ 2 + az ^ 2) + 1e-10
	nb = @. sqrt(bx ^ 2 + by ^ 2 + bz ^ 2) + 1e-10

	ab = @. ax * bx + ay * by + az * bz

	mult = @. (1.0 / na + 1.0 / nb) / (ab + na * nb)

	vx = @. (
		ay * bz - az * by
	) * mult / (4 * π * β ^ 2)
	vy = @. (
		az * bx - ax * bz
	) * mult / (4 * π * β)
	vz = @. (
		ax * by - ay * bx
	) * mult / (4 * π * β)

	[vx, vy, vz]

end

max_stars_repo_name: Equipe-AeroDesign-ITA/PropBiology

max_issues_repo_name: Equipe-AeroDesign-ITA/PropBiology

max_forks_repo_name: Equipe-AeroDesign-ITA/PropBiology


text:
using Hecke

using Test

using LinearAlgebra

for scope in Hecke.ASSERT_SCOPE
  set_assert_level(scope, 3)
end

#for scope in Hecke.VERBOSE_SCOPE
#  set_verbose_level(scope, 3)
#end

@time include("AlgAss.jl")
@time include("AlgAssOrd.jl")
@time include("EllCrv.jl")
@time include("GrpAb.jl")
@time include("LinearAlgebra.jl")
@time include("Map.jl")
@time include("Misc.jl")
@time include("NfAbs.jl")
@time include("NfOrd.jl")
@time include("NfRel.jl")
@time include("RCF.jl")
@time include("Examples.jl")
@time include("Sparse.jl")

# x^5 + 514944*x^2 + 123904 test prime decomposition with this (2 is index divisor and only one prime ideal over 2)


text:
using Query
using DataFrames
using NamedTuples

# We want a DataFrame without Nullable's, so
# we need to construct a DataFrame with columns of Array
columns=[]
push!(columns, ["John", "Sally", "Kirk"])
push!(columns, [34., 23., 59.])
push!(columns, [3,5,2])
df = DataFrame(columns, [:name, :age, :children])

x = @from i in df begin
    @orderby i.age
    @select {Name=lowercase(i.name)}
    @collect DataFrame
end

println(x)

x = @from i in df begin
    @orderby descending(i.age)
    @select {Name=lowercase(i.name)}
    @collect DataFrame
end

println(x)

x = @from i in df begin
    @orderby ascending(i.age)
    @select {Name=lowercase(i.name)}
    @collect DataFrame
end

println(x)


text:
"""
    pole(sys)

Compute the poles of system `sys`."""
pole(sys::AbstractStateSpace) = eigvalsnosort(sys.A)
pole(sys::SisoTf) = error("pole is not implemented for type $(typeof(sys))")

# Seems to have a lot of rounding problems if we run the full thing with sisorational,
# converting to zpk before works better in the cases I have tested.
pole(sys::TransferFunction) = pole(zpk(sys))

function pole(sys::TransferFunction{<:TimeEvolution,SisoZpk{T,TR}}) where {T, TR}
    # With right TR, this code works for any SisoTf

    # Calculate least common denominator of the minors,
    # i.e. something like least common multiple of the pole-polynomials
    individualpoles = [map(pole, sys.matrix)...;]
    lcmpoles = TR[]
    for poles = minorpoles(sys.matrix)
        # Poles have to be equal to existing poles for the individual transfer functions and this
        # calculation probably is more precise than the full. Seems to work better at least.
        for i = 1:length(poles)
            idx = argmin(map(abs, individualpoles .- poles[i]))
            poles[i] = individualpoles[idx]
        end
        for pole = lcmpoles
            idx = findfirst(poles .≈ pole)
            if idx != nothing
                deleteat!(poles, idx)
            end
        end
        append!(lcmpoles, poles)
    end

    return lcmpoles
end

"""
    minorpoles(sys)

Compute the poles of all minors of the system."""
# TODO: Improve implementation, should be more efficient ways.
# Calculates the same minors several times in some cases.
function minorpoles(sys::Matrix{SisoZpk{T, TR}}) where {T, TR}
    minors = Array{TR,1}[]
    ny, nu = size(sys)
    if ny == nu == 1
        push!(minors, pole(sys[1, 1]))
    elseif ny == nu
        push!(minors, pole(det(sys)))
        for i = 1:ny
            for j = 1:nu
                newmat = sys[1:end .!=i, 1:end .!= j]
                append!(minors, minorpoles(newmat))
            end
        end
    elseif ny < nu
        for i = 1:nu
            newmat = sys[1:end, 1:end .!= i]
            append!(minors, minorpoles(newmat))
        end
    else
        for i = 1:ny
            newmat = sys[1:end .!= i, 1:end]
            append!(minors, minorpoles(newmat))
        end
    end
    return minors
end

"""
    det(sys)

Compute the determinant of the Matrix `sys` of SisoTf systems, returns a SisoTf system."""
# TODO: improve this implementation, should be more efficient ones
function det(sys::Matrix{S}) where {S<:SisoZpk}
    ny, nu = size(sys)
    @assert ny == nu "Matrix is not square"
    if ny == 1
        return sys[1, 1]
    end
    tot = zero(S)
    sign = -1
    for i = 1:ny
        sign = -sign
        tot += sign * sys[i, 1] * det(sys[1:end .!= i, 2:end])
    end
    return tot
end

"""
    dcgain(sys)

Compute the dcgain of system `sys`.

equal to G(0) for continuous-time systems and G(1) for discrete-time systems."""
function dcgain(sys::LTISystem)
    return iscontinuous(sys) ? evalfr(sys, 0) : evalfr(sys, 1)
end

"""
    markovparam(sys, n)

Compute the `n`th markov parameter of state-space system `sys`. This is defined
as the following:

`h(0) = D`

`h(n) = C*A^(n-1)*B`"""
function markovparam(sys::AbstractStateSpace, n::Integer)
    n < 0 && error("n must be >= 0")
    return n == 0 ? sys.D : sys.C * sys.A^(n-1) * sys.B
end

"""
    z, p, k = zpkdata(sys)

Compute the zeros, poles, and gains of system `sys`.

### Returns
`z` : Matrix{Vector{ComplexF64}}, (ny x nu)

`p` : Matrix{Vector{ComplexF64}}, (ny x nu)

`k` : Matrix{Float64}, (ny x nu)"""
function zpkdata(sys::LTISystem)
    G = convert(TransferFunction{typeof(timeevol(sys)),SisoZpk}, sys)

    zs = map(x -> x.z, G.matrix)
    ps = map(x -> x.p, G.matrix)
    ks = map(x -> x.k, G.matrix)

    return zs, ps, ks
end

"""
    Wn, zeta, ps = damp(sys)

Compute the natural frequencies, `Wn`, and damping ratios, `zeta`, of the
poles, `ps`, of `sys`"""
function damp(sys::LTISystem)
    ps = pole(sys)
    if isdiscrete(sys)
        ps = log.(ps)/sys.Ts
    end
    Wn = abs.(ps)
    order = sortperm(Wn; by=z->(abs(z), real(z), imag(z)))
    Wn = Wn[order]
    ps = ps[order]
    ζ = -cos.(angle.(ps))
    return Wn, ζ, ps
end

"""
    dampreport(sys)

Display a report of the poles, damping ratio, natural frequency, and time
constant of the system `sys`"""
function dampreport(io::IO, sys::LTISystem)
    Wn, zeta, ps = damp(sys)
    t_const = 1 ./ (Wn.*zeta)
    header =
    ("|     Pole      |   Damping     |   Frequency   | Time Constant |\n"*
    "|               |    Ratio      |   (rad/sec)   |     (sec)     |\n"*
    "+---------------+---------------+---------------+---------------+")
    println(io, header)
    if all(isreal, ps)
        for i=eachindex(ps)
            p, z, w, t = ps[i], zeta[i], Wn[i], t_const[i]
            Printf.@printf(io, "|  %-13.3e|  %-13.3e|  %-13.3e|  %-13.3e|\n", real(p), z, w, t)
        end
    else
        for i=eachindex(ps)
            p, z, w, t = ps[i], zeta[i], Wn[i], t_const[i]
            Printf.@printf(io, "|  %-13.3e|  %-13.3e|  %-13.3e|  %-13.3e|\n", real(p), z, w, t)
            Printf.@printf(io, "|  %-+11.3eim|               |               |               |\n", imag(p))
        end
    end
end
dampreport(sys::LTISystem) = dampreport(stdout, sys)

"""
    tzero(sys)

Compute the invariant zeros of the system `sys`. If `sys` is a minimal
realization, these are also the transmission zeros."""
function tzero(sys::TransferFunction)
    if issiso(sys)
        return tzero(sys.matrix[1,1])
    else
        return tzero(ss(sys))
    end
end

# Implements the algorithm described in:
# Emami-Naeini, A. and P. Van Dooren, "Computation of Zeros of Linear
# Multivariable Systems," Automatica, 18 (1982), pp. 415–430.
#
# Note that this returns either Vector{ComplexF32} or Vector{Float64}
tzero(sys::AbstractStateSpace) = tzero(sys.A, sys.B, sys.C, sys.D)
# Make sure everything is BlasFloat
function tzero(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, D::AbstractMatrix)
    T = promote_type(eltype(A), eltype(B), eltype(C), eltype(D))
    A2, B2, C2, D2, _ = promote(A,B,C,D, fill(zero(T)/one(T),0,0)) # If Int, we get Float64
    tzero(A2, B2, C2, D2)
end
function tzero(A::AbstractMatrix{T}, B::AbstractMatrix{T}, C::AbstractMatrix{T}, D::AbstractMatrix{T}) where {T <: Union{AbstractFloat,Complex{<:AbstractFloat}}#= For eps(T) =#}
    # Balance the system
    A, B, C = balance_statespace(A, B, C)

    # Compute a good tolerance
    meps = 10*eps(real(T))*norm([A B; C D])

    # Step 1:
    A_r, B_r, C_r, D_r = reduce_sys(A, B, C, D, meps)

    # Step 2: (conjugate transpose should be avoided since single complex zeros get conjugated)
    A_rc, B_rc, C_rc, D_rc = reduce_sys(transpose(A_r), transpose(C_r), transpose(B_r), transpose(D_r), meps)
    if isempty(A)   return complex(T)[]    end

    # Step 3:
    # Compress cols of [C D] to [0 Df]
    mat = [C_rc D_rc]
    # To ensure type-stability, we have to annote the type here, as qrfact
    # returns many different types.
    W = qr(mat').Q
    W = reverse(W, dims=2)
    mat = mat*W
    if fastrank(mat', meps) > 0
        nf = size(A_rc, 1)
        m = size(D_rc, 2)
        Af = ([A_rc B_rc] * W)[1:nf, 1:nf]
        Bf = ([Matrix{T}(I, nf, nf) zeros(nf, m)] * W)[1:nf, 1:nf]
        zs = eigvalsnosort(Af, Bf)
        _fix_conjugate_pairs!(zs) # Generalized eigvals does not return exact conj. pairs
    else
        zs = complex(T)[]
    end
    return zs
end

"""
Implements REDUCE in the Emami-Naeini & Van Dooren paper. Returns transformed
A, B, C, D matrices. These are empty if there are no zeros.
"""
function reduce_sys(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, D::AbstractMatrix, meps::AbstractFloat)
    T = promote_type(eltype(A), eltype(B), eltype(C), eltype(D))
    Cbar, Dbar = C, D
    if isempty(A)
        return A, B, C, D
    end
    while true
        # Compress rows of D
        U = qr(D).Q
        D = U'D
        C = U'C
        sigma = fastrank(D, meps)
        Cbar = C[1:sigma, :]
        Dbar = D[1:sigma, :]
        Ctilde = C[(1 + sigma):end, :]
        if sigma == size(D, 1)
            break
        end

        # Compress columns of Ctilde
        V = qr(Ctilde').Q
        V = reverse(V,dims=2)
        Sj = Ctilde*V
        rho = fastrank(Sj', meps)
        nu = size(Sj, 2) - rho

        if rho == 0
            break
        elseif nu == 0
            # System has no zeros, return empty matrices
            A = B = Cbar = Dbar = Matrix{T}(undef, 0,0)
            break
        end
        # Update System
        n, m = size(B)
        Vm = [V zeros(T, n, m); zeros(T, m, n) Matrix{T}(I, m, m)]
        if sigma > 0
            M = [A B; Cbar Dbar]
            Vs = [V' zeros(T, n, sigma) ; zeros(T, sigma, n) Matrix{T}(I, sigma, sigma)]
        else
            M = [A B]
            Vs = V'
        end
        sigma, rho, nu
        M = Vs * M * Vm
        A = M[1:nu, 1:nu]
        B = M[1:nu, (nu + rho + 1):end]
        C = M[(nu + 1):end, 1:nu]
        D = M[(nu + 1):end,  (nu + rho + 1):end]
    end
    return A, B, Cbar, Dbar
end

# Determine the number of non-zero rows, with meps as a tolerance. For an
# upper-triangular matrix, this is a good proxy for determining the row-rank.
function fastrank(A::AbstractMatrix, meps::Real)
    n, m = size(A)
    if n*m == 0     return 0    end
    norms = Vector{real(eltype(A))}(undef, n)
    for i = 1:n
        norms[i] = norm(A[i, :])
    end
    mrank = sum(norms .> meps)
    return mrank
end

"""
`ωgₘ, gₘ, ωϕₘ, ϕₘ = margin{S<:Real}(sys::LTISystem, w::AbstractVector{S}; full=false, allMargins=false)`

returns frequencies for gain margins, gain margins, frequencies for phase margins, phase margins

If `!allMargins`, return only the smallest margin

If `full` return also `fullPhase`

"""
function margin(sys::LTISystem, w::AbstractVector{<:Real}; full=false, allMargins=false)
    ny, nu = size(sys)

    if allMargins
        wgm         = Array{Array{numeric_type(sys),1}}(undef, ny,nu)
        gm          = Array{Array{numeric_type(sys),1}}(undef, ny,nu)
        wpm         = Array{Array{numeric_type(sys),1}}(undef, ny,nu)
        pm          = Array{Array{numeric_type(sys),1}}(undef, ny,nu)
        fullPhase   = Array{Array{numeric_type(sys),1}}(undef, ny,nu)
    else
        wgm         = Array{numeric_type(sys),2}(undef, ny, nu)
        gm          = Array{numeric_type(sys),2}(undef, ny, nu)
        wpm         = Array{numeric_type(sys),2}(undef, ny, nu)
        pm          = Array{numeric_type(sys),2}(undef, ny, nu)
        fullPhase   = Array{numeric_type(sys),2}(undef, ny, nu)
    end
    for j=1:nu
        for i=1:ny
            wgm[i,j], gm[i,j], wpm[i,j], pm[i,j], fullPhase[i,j] = sisomargin(sys[i,j], w, full=true, allMargins=allMargins)
        end
    end
    if full
        wgm, gm, wpm, pm, fullPhase
    else
        wgm, gm, wpm, pm
    end
end

"""
`ωgₘ, gₘ, ωϕₘ, ϕₘ = sisomargin{S<:Real}(sys::LTISystem, w::AbstractVector{S}; full=false, allMargins=false)`

returns frequencies for gain margins, gain margins, frequencies for phase margins, phase margins
"""
function sisomargin(sys::LTISystem, w::AbstractVector{<:Real}; full=false, allMargins=false)
    ny, nu = size(sys)
    if ny !=1 || nu != 1
        error("System must be SISO, use `margin` instead")
    end
    mag, phase, w = bode(sys, w)
    wgm, = _allPhaseCrossings(w, phase)
    gm = similar(wgm)
    for i = eachindex(wgm)
        gm[i] = 1 ./ abs(freqresp(sys,[wgm[i]])[1][1])
    end
    wpm, fi = _allGainCrossings(w, mag)
    pm = similar(wpm)
    for i = eachindex(wpm)
        pm[i] = mod(rad2deg(angle(freqresp(sys,[wpm[i]])[1][1])),360)-180
    end
    if !allMargins #Only output the smallest margins
        gm, idx = findmin([gm;Inf])
        wgm = [wgm;NaN][idx]
        fi = [fi;NaN][idx]
        pm, idx = findmin([abs.(pm);Inf])
        wpm = [wpm;NaN][idx]
        if full
            if !isnan(fi) #fi may be NaN, fullPhase is a scalar
                fullPhase = interpolate(fi, phase)
            else
                fullPhase = NaN
            end
        end
    else
        if full #We know that all values are defined and fullPhase is a vector
            fullPhase = interpolate(fi, phase)
        end
    end
    if full
        wgm, gm, wpm, pm, fullPhase
    else
        wgm, gm, wpm, pm
    end
end
margin(system::LTISystem; kwargs...) =
margin(system, _default_freq_vector(system, Val{:bode}()); kwargs...)
#margin(sys::LTISystem, args...) = margin(LTISystem[sys], args...)

# Interpolate the values in "list" given the floating point "index" fi
function interpolate(fi, list)
    fif = floor.(Integer, fi)
    fic = ceil.(Integer, fi)
    list[fif]+mod.(fi,1).*(list[fic]-list[fif])
end

function _allGainCrossings(w, mag)
    _findCrossings(w,mag.>1,mag.-1)
end

function _allPhaseCrossings(w, phase)
    #Calculate numer of times real axis is crossed on negative side
    n = fld.(phase.+180,360) #Nbr of crossed
    ph = mod.(phase,360) .- 180 #Residual
    _findCrossings(w, n, ph)
end

function _findCrossings(w, n, res)
    wcross = Vector{eltype(w)}()
    tcross = Vector{eltype(w)}()
    for i in 1:(length(w)-1)
        if res[i] == 0
            wcross = [wcross; w[i]]
            tcross = [tcross; i]
        elseif n[i] != n[i+1]
            #Interpolate to approximate crossing
            t = res[i]/(res[i]-res[i+1])
            tcross = [tcross; i+t]
            wt = w[i]+t*(w[i+1]-w[i])
            wcross = [wcross; wt]
        end
    end
    if res[end] == 0 #Special case if multiple points
        wcross = [wcross; w[end]]
        tcross = [tcross; length(w)]
    end
    wcross, tcross
end

"""
    dₘ = delaymargin(G::LTISystem)

Only supports SISO systems"""
function delaymargin(G::LTISystem)
    # Phase margin in radians divided by cross-over frequency in rad/s.
    if G.nu + G.ny > 2
        error("delaymargin only supports SISO systems")
    end
    m     = margin(G,allMargins=true)
    ϕₘ, i = findmin(m[4])
    ϕₘ   *= π/180
    ωϕₘ   = m[3][i]
    dₘ    = ϕₘ/ωϕₘ
    if isdiscrete(G)
        dₘ /= G.Ts # Give delay margin in number of sample times, as matlab does
    end
    dₘ
end

function robust_minreal(G, args...; kwargs...)
    try
        return minreal(G, args...; kwargs...)
    catch
        return G
    end
end

"""
    S,D,N,T = gangoffour(P,C; minimal=true)
    gangoffour(P::AbstractVector,C::AbstractVector; minimal=true)

Given a transfer function describing the Plant `P` and a transfer function describing the controller `C`, computes the four transfer functions in the Gang-of-Four.

- `minimal` determines whether or not to call `minreal` on the computed systems.
- `S = 1/(1+PC)` Sensitivity function
- `D = P/(1+PC)`
- `N = C/(1+PC)`
- `T = PC/(1+PC)` Complementary sensitivity function

Only supports SISO systems
"""
function gangoffour(P::LTISystem,C::LTISystem; minimal=true)
    if P.nu + P.ny + C.nu + C.ny > 4
        error("gangoffour only supports SISO systems")
    end
    minfun = minimal ? robust_minreal : identity
    S = (1/(1+P*C)) |> minfun
    D = (P*S)       |> minfun
    N = (C*S)       |> minfun
    T = (P*N)       |> minfun
    return S, D, N, T
end

"""
    S, D, N, T, RY, RU, RE = gangofseven(P,C,F)

Given transfer functions describing the Plant `P`, the controller `C` and a feed forward block `F`,
computes the four transfer functions in the Gang-of-Four and the transferfunctions corresponding to the feed forward.

`S = 1/(1+PC)` Sensitivity function

`D = P/(1+PC)`

`N = C/(1+PC)`

`T = PC/(1+PC)` Complementary sensitivity function

`RY = PCF/(1+PC)`

`RU = CF/(1+P*C)`

`RE = F/(1+P*C)`

Only supports SISO systems"""
function gangofseven(P::TransferFunction,C::TransferFunction,F::TransferFunction)
    if P.nu + P.ny + C.nu + C.ny + F.nu + F.ny > 6
        error("gof only supports SISO systems")
    end
    S, D, N, T = gangoffour(P,C)
    RY = T*F
    RU = N*F
    RE = S*F
    return S, D, N, T, RY, RU, RE
end


text:
# lhs = left hand side; rhs = right hand side
# src = source, dst = destination
# order is: lhs, rhs, src, dst
function Arrow(lhs, rhs, src, dst; edge_fun = DirectedEdge, cross = false)
    if cross
        edges = edge_fun(src, dst)
    else
        edges = edge_fun.(src, dst)
    end
    edges = isa(edges, VecOrMat) ? edges : [edges]
    edges = isa(edges, Matrix) ? vec(edges) : edges
    edges = isa(lhs, Arrow) ? vcat(lhs.edges..., edges...) : edges
    edges = isa(rhs, Arrow) ? vcat(edges..., rhs.edges...) : edges
    Arrow(edges, lhs, rhs)
end

"""
Single (broadcasting) arrow right (\\rightarrow, Alt Gr + i)

```jldoctest
julia> @StenoGraph a → b
a → b

julia> eltype(ans)
DirectedEdge{SimpleNode{Symbol}, SimpleNode{Symbol}}

julia> @StenoGraph [a b] → [c d]
a → c
b → d
```
"""
→(lhs, rhs) = Arrow(lhs, rhs, keep(lhs, Right, ModifiedNode), keep(rhs, Left, AbstractNode))

"""
Double (cross product) arrow right (\\Rightarrow)

```jldoctest
julia> @StenoGraph a ⇒ b
a → b

julia> eltype(ans)
DirectedEdge{SimpleNode{Symbol}, SimpleNode{Symbol}}

julia> @StenoGraph [a b] ⇒ [c d]
a → c
a → d
b → c
b → d
```
"""
⇒(lhs, rhs) = Arrow(lhs, rhs, keep(lhs, Right, ModifiedNode), keep(rhs, Left, AbstractNode); cross = true)

"""
Single (broadcasting) arrow left (\\leftarrow, Alt Gr + z)

```jldoctest
julia> @StenoGraph a ← b
b → a

julia> eltype(ans)
DirectedEdge{SimpleNode{Symbol}, SimpleNode{Symbol}}

julia> @StenoGraph [a b] ← [c d]
c → a
d → b
```
"""
←(lhs, rhs) = Arrow(lhs, rhs, keep(rhs, Left, ModifiedNode), keep(lhs, Right, AbstractNode))

"""
Double (cross product) arrow left (\\Leftarrow)

```jldoctest
julia> @StenoGraph a ⇐ b
b → a

julia> eltype(ans)
DirectedEdge{SimpleNode{Symbol}, SimpleNode{Symbol}}

julia> @StenoGraph [a b] ⇐ [c d]
c → a
c → b
d → a
d → b
```
"""
⇐(lhs, rhs) = Arrow(lhs, rhs, keep(rhs, Left, ModifiedNode), keep(lhs, Right, AbstractNode); cross = true)

"""
Single (broadcasting) arrow left-right (\\leftrightarrow)

```jldoctest
julia> @StenoGraph a ↔ b
a ↔ b

julia> eltype(ans)
UndirectedEdge{SimpleNode{Symbol}, SimpleNode{Symbol}}

julia> @StenoGraph [a b] ↔ [c d]
a ↔ c
b ↔ d
```
"""
↔(lhs, rhs) = Arrow(lhs, rhs, keep(lhs, Right, ModifiedNode), keep(rhs, Left, AbstractNode); edge_fun = UndirectedEdge)

"""
Double (cross product) arrow left-right (\\Leftrightarrow)

```jldoctest
julia> @StenoGraph a ⇔ b
a ↔ b

julia> eltype(ans)
UndirectedEdge{SimpleNode{Symbol}, SimpleNode{Symbol}}

julia> @StenoGraph [a b] ⇔ [c d]
a ↔ c
a ↔ d
b ↔ c
b ↔ d
```
"""
⇔(lhs, rhs) = Arrow(lhs, rhs, keep(lhs, Right, ModifiedNode), keep(rhs, Left, AbstractNode); edge_fun = UndirectedEdge, cross = true)

# if either lhs or rhs are a single node, broadcast arrow (barrow) behave like cross arrow (carrow)
barrows = [:→, :←, :↔]
carrows = [:⇒, :⇐, :⇔]
for (barrow, carrow) in zip(barrows,  carrows)
    @eval $barrow(lhs::Union{AbstractNode, Symbol}, rhs::Union{AbstractNode, Symbol}) = $carrow(lhs, rhs)
    @eval $barrow(lhs::Union{AbstractNode, Symbol}, rhs) = $carrow(lhs, rhs)
    @eval $barrow(lhs, rhs::Union{AbstractNode, Symbol}) = $carrow(lhs, rhs)
end


text:
using LinearAlgebra,
      SparseArrays

function normalize_matrix(A::AbstractMatrix, dim::Integer; f=LinearAlgebra.norm)
    if (size(A, dim) == 1)
        return A
    end
    return A./mapslices(f, A, dims=[dim])
end

normalize_rows(A::AbstractMatrix; f=LinearAlgebra.norm) = normalize_matrix(A, 2, f=f)
normalize_cols(A::AbstractMatrix; f=LinearAlgebra.norm) = normalize_matrix(A, 1, f=f)

function normalize_matrix!(A::AbstractMatrix, dim::Integer; f=LinearAlgebra.norm)
    if (size(A, dim) == 1)
        return A
    end
    A ./= mapslices(f, A, dims=[dim])
end
normalize_rows!(A::AbstractMatrix; f=LinearAlgebra.norm) = normalize_matrix!(A, 2, f=f)
normalize_cols!(A::AbstractMatrix; f=LinearAlgebra.norm) = normalize_matrix!(A, 1, f=f)
#=
# get the column sums of A
S = vec(sum(A,1))

# get the nonzero entries in A. ei is row index, ej is col index, ev is the value in A
ei,ej,ev = findnz(A)

# get the number or rows and columns in A
m,n = size(A)

# create a new normalized matrix. For each nonzero index (ei,ej), its new value will be
# the old value divided by the sum of that column, which can be obtained by S[ej]
A_normalized = sparse(ei,ej,ev./S[ej],m,n)
http://stackoverflow.com/questions/24296856/in-julia-how-can-i-column-normalize-a-sparse-matrix
=#
function normalize_cols(A::SparseMatrixCSC)
    sums = sum(A, 1) + eps()
    I, J, V = findnz(A)
    for idx in 1:length(V)
        V[idx] /= sums[J[idx]]
    end
    sparse(I, J, V)
end

function svd_whiten(X)
    U, s, Vt = svd(X)
    return U * Vt
end


text:
using Printf, Distributions

newgrid(p::Float64, r::Int, c::Int=r) = rand(Bernoulli(p), r, c)

function walkmaze!(grid::Matrix{Int}, r::Int, c::Int, indx::Int)
    NOT_CLUSTERED = 1 # const
    N, M = size(grid)
    dirs = [[1, 0], [-1, 0], [0, 1], [0, -1]]
    # fill cell
    grid[r, c] = indx
    # check for each direction
    for d in dirs
        rr, cc = (r, c) .+ d
        if checkbounds(Bool, grid, rr, cc) && grid[rr, cc] == NOT_CLUSTERED
            walkmaze!(grid, rr, cc, indx)
        end
    end
end

function clustercount!(grid::Matrix{Int})
    NOT_CLUSTERED = 1 # const
    walkind = 1
    for r in 1:size(grid, 1), c in 1:size(grid, 2)
        if grid[r, c] == NOT_CLUSTERED
            walkind += 1
            walkmaze!(grid, r, c, walkind)
        end
    end
    return walkind - 1
end
clusterdensity(p::Float64, n::Int) = clustercount!(newgrid(p, n)) / n ^ 2

function printgrid(G::Matrix{Int})
    LETTERS = vcat(' ', '#', 'A':'Z', 'a':'z')
    for r in 1:size(G, 1)
        println(r % 10, ") ", join(LETTERS[G[r, :] .+ 1], ' '))
    end
end

G = newgrid(0.5, 15)
@printf("Found %i clusters in this %i×%i grid\n\n", clustercount!(G), size(G, 1), size(G, 2))
printgrid(G)
println()

const nrange = 2 .^ (4:2:12)
const p = 0.5
const nrep = 5
for n in nrange
    sim = mean(clusterdensity(p, n) for _ in 1:nrep)
    @printf("nrep = %2i p = %.2f dim = %-13s sim = %.5f\n", nrep, p, "$n × $n", sim)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
export build_table

# bm_table = Matrix{Tuple{Int,Int}}(undef, 2, 2)
bm_table = Vector{Tuple{Int8,Int8}}(undef, 2)

bm_flag = false

@inline function bm_table_save(c::Int, d::Int, bm::Tuple{Int,Int})
    @inbounds bm_table[10_000*c+d+1] = bm
end

@inline function bm_table_fetch(c::Int, d::Int)::Tuple{Int,Int}
    @inbounds bm_table[10_000*c+d+1]
end

"""
`build_table()` creates a lookup table for `fast_bm_count()`
"""
function build_table(verbose::Bool = true)
    if !bm_flag
        global bm_flag = true
    else
        return
    end
    @info "Building lookup table"
    global bm_table = Vector{Tuple{Int8,Int8}}(undef, 10_000 * 10_000)

    if verbose
        PM = Progress(50005000)
    end
    for a = 0:9999
        for b = a:9999
            bm = slow_bm_count(a, b)
            bm_table_save(a, b, bm)
            bm_table_save(b, a, bm)
            if verbose
                next!(PM)
            end
        end
    end
end


text:
using UTFDetect
using Test

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

max_stars_repo_name: StefanKarpinski/UTFDetect

max_issues_repo_name: StefanKarpinski/UTFDetect

max_forks_repo_name: StefanKarpinski/UTFDetect


text:
# Generate Data
A_ssp = [-4.0, 3.0, -8.0]
A_ppp = [8.0, -6.0, 3]
Γ     = [9.0,  8.0, 7.0]
ω     = [2870, 2910, 2940]
a_ssp = [1.0, 1.0, 1.0,
        0.99, 0.98, 0.97,
        0.98, 0.96, 0.94]
a_pow = 0.9   # factor by which we power the a_ssp values to
               # to get the a_ppp values

a_ssp = reshape(a_ssp, (3,3))' # num time steps, num resonances
a_ppp = a_ssp .^ a_pow
wn    = range(2830, stop=2990, length=301) |> collect

sig_ssp = zeros(size(a_ssp,1), length(wn))
sig_ppp = similar(sig_ssp)

for i = 1:size(a_ssp,1), j = 1:length(wn)
    sig_ssp[i,j] = FiSUDiPPuS.sfspec(wn[j], a_ssp[i,:] .* A_ssp, ω, Γ)
    sig_ppp[i,j] = FiSUDiPPuS.sfspec(wn[j], a_ppp[i,:] .* A_ppp, ω, Γ)
end

# Save Data
savepath = joinpath(@__DIR__, "../data")
d_ssp = Dict(
    "wavenumber" => wn,
    "signal"     => sig_ssp,
    "dltime"     => range(-3, 5, length=length(a_ssp)÷3)
)
d_ppp = Dict(
    "wavenumber" => wn,
    "signal"     => sig_ppp,
    "dltime"     => range(-3, 5, length=length(a_ssp)÷3)
)

save(joinpath(savepath, "sample_ssp.jld2"), d_ssp)
save(joinpath(savepath, "sample_ppp.jld2"), d_ppp)


text:
module MNIST
export trainimages, testimages,
       trainlabels, testlabels,
       traindata,   testdata

include("getdata.jl")

end # module


text:
using KernelDensity, Distributions, Interpolations
import KernelDensity: kernel_dist
kernel_dist(::Type{Epanechnikov}, w::Float64) = Epanechnikov(0.0, w)
kernel(data) = kde(data; kernel=Epanechnikov)


text:
"""
What is the distribution of static parameters of the SD filter,  when used as DGP? is it well approximated by the MLE normal?
"""

#region import and models
begin

using ScoreDrivenERGM
import ScoreDrivenERGM:StaticNets, DynNets
import ScoreDrivenERGM.DynNets:SdErgm,SdErgmDirBin0Rec0, sample_ergm_sequence, stats_from_mat, array_2_vec_all_par, unrestrict_all_par, conf_bands_par_uncertainty, avg_grad_and_hess_obj_SD_filter_time_seq, conf_bands_par_uncertainty, number_ergm_par, estimate_filter_and_conf_bands, conf_bands_coverage, estimate, mle_distrib_filtered_par, plot_filtered_and_conf_bands
using ScoreDrivenERGM.Utilities

using PyPlot
pygui(true)

using ForwardDiff
using StatsBase
using LinearAlgebra
using Distributions
using Statistics
using JLD2

model_mle = DynNets.SdErgmDirBin0Rec0_mle()
model_pmle = DynNets.SdErgmDirBin0Rec0_pmle()
indTvPar = trues(2)

end
#endregion

#region
begin
T=300
N=300
nSample = 100
model = model_mle

quantilesVals = [[0.975, 0.025] ]

# sample dgp
dgpSetAR, ~, dgpSetSD = ScoreDrivenERGM.DynNets.list_example_dgp_settings(model_mle)

dgpSettings = dgpSetSD
parDgpT = DynNets.sample_time_var_par_from_dgp(model, dgpSettings.type, N, T;  dgpSettings.opt...)

A_T_dgp = DynNets.sample_ergm_sequence(model,N , parDgpT, 1)[:,:,:,1]

@show std(parDgpT, dims=2)

res_mle_pmle = DynNets.estimate_filter_and_conf_bands(model, A_T_dgp, quantilesVals; plotFlag =true, parDgpT = parDgpT, parUncMethod = "WHITE-MLE")

end

res_mle_pmle[1]

begin
T=300
N=300
nSample = 100
model = model_pmle
alpha = 0.2
beta = mean(DynNets.beta_min_max_from_alpha_min(0.1, N))
UM = [alpha, beta]
vUnPar, ~ = DynNets.starting_point_optim(model,  UM)
vResParDgp = DynNets.restrict_all_par(model, vUnPar)

BDgp = 0.96
ADgp = 0.2
vResParDgp[2:3:end] .= BDgp
vResParDgp[3:3:end] .= ADgp

#start optimization on the correct values
vParUnOptim_0, ARe_min = DynNets.starting_point_optim(model, UM; indTargPar = falses(number_ergm_par(model)))

vParOptim_0 = Real.(DynNets.restrict_all_par(model, vParUnOptim_0))

nStaticSdPar = length(vResParDgp)
vResParEstDistrib = zeros(length(vResParDgp), nSample)
vUnParEstDistrib = zeros(length(vResParDgp), nSample)
HessSum = zeros(nStaticSdPar, nStaticSdPar, nSample)
OpGradSum = zeros(nStaticSdPar, nStaticSdPar, nSample)

quantilesVals = [[0.975, 0.025] ]

for n=1:2
    # sample SD dgp
    global fVecTDgp, A_T_dgp, ~ = DynNets.score_driven_filter_or_dgp( model_mle,N,  vResParDgp, indTvPar; dgpNT = (N,T))

    #global res_mle_mle = DynNets.estimate_filter_and_conf_bands(model_mle, A_T_dgp, quantilesVals; plotFlag =true, parDgpT = fVecTDgp, parUncMethod = "WHITE-MLE")

    global res_mle_pmle = DynNets.estimate_filter_and_conf_bands(model_pmle, A_T_dgp, quantilesVals; plotFlag =true, parDgpT = fVecTDgp, parUncMethod = "WHITE-MLE")

   # global res_mle_boot = DynNets.estimate_filter_and_conf_bands(model, A_T_dgp, quantilesVals; plotFlag =true, parDgpT = fVecTDgp, parUncMethod = "PB-SAMPLE")

    # global res_mle_boot = DynNets.estimate_filter_and_conf_bands(model, A_T_dgp, quantilesVals; plotFlag =true, parDgpT = fVecTDgp, parUncMethod = "PB-COVMAT")

end
end

for n=1:nSample
    # sample SD dgp
    fVecTDgp, A_T_dgp, ~ = DynNets.score_driven_filter_or_dgp( model, N, vResParDgp, indTvPar; dgpNT = (N,T))

    obsT = [stats_from_mat(model, A_T_dgp[:,:,t]) for t in 1:T ]

    arrayAllParHat, conv_flag,UM, ftot_0 = estimate(model, N, obsT; indTvPar=indTvPar, indTargPar=falses(2), vParOptim_0=vParOptim_0)

    vResParEstDistrib[:,n] = DynNets.array_2_vec_all_par(model, arrayAllParHat, indTvPar)

    OpGradSum[:,:,n], HessSum[:,:,n] = DynNets.A0_B0_est_for_white_cov_mat_obj_SD_filter_time_seq(model, N, obsT, vResParEstDistrib[:,n], ftot_0)

    # estCovWhite[:,:,n], errorFlag, OpGradSum[:,:,n], HessSum[:,:,n] = DynNets.white_estimate_cov_mat_static_sd_par(model, obsT, indTvPar, ftot_0, vResParEstDistrib[:,n]; returnAllMats=true)

end

# end
#endregion
n=1
begin
figure()
fVecTDgp, A_T_dgp, ~ = DynNets.score_driven_filter_or_dgp( model, vResParDgp, indTvPar; dgpNT = (N,T))

obsT = [stats_from_mat(model, A_T_dgp[:,:,t]) for t in 1:T ]

arrayAllParHat, conv_flag,UM, ftot_0 = estimate(model, obsT; indTvPar=indTvPar, indTargPar=falses(2), vParOptim_0=vParOptim_0)

vResParEstDistrib[:,n] = DynNets.array_2_vec_all_par(model, arrayAllParHat, indTvPar)

fVecTFilt, A_T_dgp, ~ = DynNets.score_driven_filter_or_dgp( model, vResParEstDistrib[:,n], indTvPar, obsT=obsT, ftot_0=ftot_0)

plot(fVecTFilt[2,:])
plot(fVecTDgp[2,:])
end

begin
indPar = 3
fig, ax = subplots()
ax.hist(vResParEstDistrib[indPar,:], 30)
ylim = ax.get_ylim()
ax.vlines(vResParDgp[indPar], ylim[1], ylim[2], color = "r" )
ax.vlines(vResPar_0[indPar], ylim[1], ylim[2], color = "r" )
#region distribution of filtered parameters checks
end

#region distribution of white estimators
# 1. C'é un bias di campione finito. Lo stimatore, MLE o White, della covarianza degli stimatori dei parametri é ok?
begin
estCovWhite = zeros(nStaticSdPar, nStaticSdPar, nSample)
estCovWhitePos = zeros(nStaticSdPar, nStaticSdPar, nSample)
errFlag = falses(nSample)
for n=1:nSample
    opGrad = OpGradSum[:,:,n]
    hess = HessSum[:,:,n]

        parCovHat = pinv(hess) * opGrad * pinv(hess)

        parCovHatPosDef, minEigenVal = make_pos_def(parCovHat)

        estCovWhite[:,:,n] = parCovHat
        estCovWhitePos[:,:,n] = parCovHatPosDef
end

parBootCov = cov(vResParEstDistrib')
parNames = ["w_theta", "B_theta", "A_theta", "w_eta", "B_eta", "A_eta"]

fig, ax = subplots(3,2)
for i = 1:6
    j=i
    ax[i].hist(estCovWhite[i,j,:], range=quantile(estCovWhite[i,j,:], [0.01, 0.99]))
    ylim = ax[i].get_ylim()
    ax[i].vlines(parBootCov[i,j], ylim[1], ylim[2], color = "r" )
    ax[i].set_title(parNames[i])
end
suptitle("White Estimators of Diagonal Elements of Covariance Matrix \n Par Bootstrap Estimate in Red. nSample = $nSample ")

tight_layout()
end
#endregion

# Domande
# 1. C'é un bias di campione finito. Lo stimatore, MLE o White, della covarianza degli stimatori dei parametri é ok?
# 2. Visto che posso costruire la distribuzione dei parametri con parametric bootstrap, perché non usare quella?
# 3. Che relazione c'é tra la distribuzione da parametric bootstrap e la mia idea di usare non parametric bootstrap campionando le osservazioni (che portano in dote la storia fino a quel punto, ma considerando la likelihood della sola osservazione)?

#endregion

max_stars_repo_name: domenicodigangi/SDERGM

max_issues_repo_name: domenicodigangi/SDERGM

max_forks_repo_name: domenicodigangi/SDERGM


text:
# run from project directory

# use this line to test coverage
# include("test/rm_cov.jl");Pkg.test("PlanerTransformer",coverage=true)

cd("./src")
for (root, dirs, files) in walkdir(".")
  for file in files
    (a,b) = splitext(file)
    if b == ".cov"
      rm(file)
    end
  end
end
cd("..")

cd("./test")
for (root, dirs, files) in walkdir(".")
  for file in files
    (a,b) = splitext(file)
    if b == ".cov"
      rm(file)
    end
  end
end
cd("..")


text:
# imgpredict.jl
# This script takes a single image and runs the random forest predictor on it
# author: Timothy Evans
# date: 04/08/14

push!(LOAD_PATH, pwd())
using DecisionTree, JLD, Images, DataFrames, Winston.imagesc, Winston.figure
importall config, dataprocessing

#load paths from config.jl
specpath,labelledpath,trainpath,modelpath = loadpaths()

println("Loading Random Forest...")
model = jldopen(string(modelpath,"model.jld"),"r") do file
  read(file,"model")
end

spec = "BSR01_140107_0_123.jpg"#readdir(specpath)[1];
print("Loading image $spec")
img = Images.load(string(trainpath,spec));
imgmask = Images.load(string(labelledpath,spec[1:end-3],"png"));
recorderdata = readtable(string(modelpath,"recorderdata.csv"));
(pixels,freq,boxmu,boxvar,label) = imgfeature(img, imgmask,spec,recorderdata);
featuredata = hcat(pixels,freq,boxmu,boxvar);

#normalise features using saved mean and standard dev
normvalues = readcsv(string(modelpath,"musigma.csv"));
mu = normvalues[1,:];
sigma = normvalues[2,:];
(m,n) = size(featuredata);
X = (featuredata - ones(m)*mu) ./ (ones(m)*sigma);

y = apply_forest(model,X);
ytest = convertclasses(label); # make zeros and ones all zeros and twos into ones

#
gt = reshapepred(ytest,width(img),height(img));
predimg = convert(Array{Float64,2},reshapepred(y,width(img),height(img)));
predimg = convthresh(predimg)
figure()
imagesc(gt)
figure()
imagesc(predimg)

confm = confusmat(ytest,y);
acc = trace(confm)./length(ytest);
C = confm./sum(confm,2);
@printf "\nConfusion Matrix:\n\n\t[%.3f , %.3f\n\t%.3f , %.3f]\n\nConfusion Accuracies:\n\n\t[%d , %d\n\n\t %d , %d]\n\nTest Accuracy: %.3f" C[1,1] C[1,2] C[2,1] C[2,2] confm[1,1] confm[1,2] confm[2,1] confm[2,2] acc

max_stars_repo_name: kokakodetect/kokako-detection

max_issues_repo_name: kokakodetect/kokako-detection

max_forks_repo_name: kokakodetect/kokako-detection


text:
### CR3BP Functions ###
V_pseudo(r) = _V_pseudo(r,μ)
CR3BP_ode!(t, x, dx) = _CR3BP_ode!(t, x, dx, μ)
C_Jacobi(r,v) = _C_Jacobi(r,v,μ)
vy(r, vx; E=C_L1) = _vy(r, vx, μ, E=E)
### CR3BP params ###
L_1  =  1 - (μ/3)^(1/3)
L_2  =  1 + (μ/3)^(1/3)
L_3  =  -(1 + (5μ/12))
L_x4 = 1/2 * (1 - 2μ)
L_x5 = 1/2 * (1 - 2μ)

C_L1 = V_pseudo([L_1,0.0])
C_L2 = V_pseudo([L_2,0.0])
C_L3 = V_pseudo([L_3,0.0])
C_L4 = V_pseudo([L_x4,L_y4])
C_L5 = V_pseudo([L_x5,L_y5])

println("L_x4: $(round(L_x4,3))")
println("C_L1: $(round(C_L1,3))")

### Integration Parameters ###
abstol = 1e-7
ord_taylor = 20
ord_jet = 4

println("\nTaylorIntegration parameters:
abstol = $abstol
ord_taylor = $ord_taylor
ord_jet = $ord_jet")

### Independent configuration variables ###
δq1, δq2 = set_variables("x y",order=ord_jet)

println("\nIndependent Variables:
δq1: $δq1
δq2: $δq2")

println("\nμ: $μ")

max_stars_repo_name: blas-ko/tesis

max_issues_repo_name: blas-ko/tesis

max_forks_repo_name: blas-ko/tesis


text:
"""
**Represents an occurrence in the GBIF format**

This is currently a subset of all the fields. This `struct` is *not* mutable --
this ensures that the objects returned from the GBIF database are never modified
by the user.

The `taxon` field is a `GBIFTaxon` object, and can therefore be manipulated as
any other `GBIFTaxon`.
"""
struct GBIFRecord
    key::Integer
    datasetKey::String
    dataset::Union{Missing, String}
    publishingOrgKey::Union{Missing, AbstractString}
    publishingCountry::Union{Missing, AbstractString}
    institutionCode::Union{Missing, AbstractString}
    protocol::Union{Missing, AbstractString}
    crawled::Union{Missing, DateTime}
    parsed::Union{Missing, DateTime}
    modified::Union{Missing, DateTime}
    interpreted::Union{Missing, DateTime}
    countryCode::Union{Missing, AbstractString}
    country::Union{Missing, AbstractString}
    basisOfRecord::Symbol
    individualCount::Union{Missing, Integer}
    latitude::Union{Missing, AbstractFloat}
    longitude::Union{Missing, AbstractFloat}
    precision::Union{Missing, AbstractFloat}
    uncertainty::Union{Missing, AbstractFloat}
    geodetic::Union{Missing, AbstractString}
    date::Union{Missing, DateTime}
    issues::Array{Symbol,1}
    taxonKey::Union{Missing, Integer}
    rank::Union{Missing, String}
    taxon::GBIFTaxon
    generic::Union{Missing, AbstractString}
    epithet::Union{Missing, AbstractString}
    vernacular::Union{Missing, AbstractString}
    scientific::Union{Missing, AbstractString}
    observer::Union{Missing, AbstractString}
    license::Union{Missing, AbstractString}
end

"""
**Internal function to format dates in records**
"""
function format_date(o, d)
    t = get(o, d, missing)
    if t === nothing || ismissing(t)
        return missing
    else
        DateTime(t[1:19])
    end
end

"""
**Generates an occurrence from the JSON response of GBIF**
"""
function GBIFRecord(o::Dict{String, Any})
    # Prepare the taxon object
    levels = ["kingdom", "phylum", "class", "order", "family", "genus", "species"]
    r = Dict{Any,Any}()
    for l in levels
        if haskey(o, l)
            r[l] = Pair(o[l], o[l*"Key"])
        else
            r[l] = missing
        end
    end

    this_name = o["genericName"]
    if !ismissing(get(o, "specificEpithet", missing))
        this_name *= " "*get(o, "specificEpithet", missing)
    end

    this_record_taxon = GBIFTaxon(
        this_name,
        o["scientificName"],
        :ACCEPTED,
        :EXACT,
        r["kingdom"],
        r["phylum"],
        r["class"],
        r["order"],
        r["family"],
        r["genus"],
        r["species"],
        100,
        false
    )

    return GBIFRecord(
    o["key"],
    o["datasetKey"],
    get(o, "datasetName", missing),
    get(o, "publishingOrgKey", missing),
    get(o, "publishingCountry", missing),
    get(o, "institutionCode", missing),
    get(o, "protocol", missing),
    format_date(o, "lastCrawled"),
    format_date(o, "lastParsed"),
    format_date(o, "modified"),
    format_date(o, "lastInterpreted"),
    get(o, "countryCode", missing),
    get(o, "country", missing),
    Symbol(o["basisOfRecord"]),
    get(o, "individualCount", missing),
    get(o, "decimalLatitude", missing),
    get(o, "decimalLongitude", missing),
    get(o, "precision", missing),
    get(o, "coordinateUncertaintyInMeters", missing),
    get(o, "geodeticDatum", missing),
    format_date(o, "eventDate"),
    Symbol.(o["issues"]),
    get(o, "taxonKey", missing),
    get(o, "taxonRank", missing),
    this_record_taxon,
    get(o, "genericName", missing),
    get(o, "specificEpithet", missing),
    get(o, "vernacularName", missing),
    get(o, "scientificName", missing),
    get(o, "recordedBy", missing),
    get(o, "license", missing)
    )
end

"""
**List of occurrences and metadata**

This type has actually very few information, besides `offset` (the number of
records already retrieved) and `count` (the total number of records). The
`query` field stores the query parameters, and `show` is a vector of boolean
values to decide which of the `GBIFRecord` (stored in `occurrences`) will be
displayed.

This type is mutable and fully iterable.
"""
mutable struct GBIFRecords
    offset::Integer
    count::Integer
    query::Union{Dict{String,Any},Nothing}
    occurrences::Vector{GBIFRecord}
    show::Vector{Bool}
end


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.
# *****************************************************************************

using TensorFlow
using TensorFlow.API: gradients

type Circuit
    blks::Vector{AbstractMux}
    opnames::Vector{Vector{ASCIIString}}
end

function hardselect_tensor(ckt::Circuit)
    hardselects = [hardselect(blk) for blk in ckt.blks]
    hardtensor = transpose_(pack(Tensor(hardselects)))
    hardtensor
end

function softselect_tensor(ckt::Circuit)
    softsels = [softselect(blk) for blk in ckt.blks]
    softtensor = Tensor(softsels)
    softtensor
end

#"""
#Output a simple string concatenated by delimiter
#"""
using Debug
@debug function simplestring(sess::Session, ckt::Circuit, fd::FeedDict=FeedDict();
    order::Vector{Int64}=Int64[], delim::ASCIIString="_")

    hs = run(sess, hardselect_tensor(ckt), fd)
    n_examples, n_selectors = size(hs)

    #@bp

    if isempty(order)
        order = collect(1:n_selectors)
    end

    stringout = ASCIIString[]
    for i = 1:n_examples
        A = [ckt.opnames[j][hs[i,j]+1] for j in order] #x+1 for python encoding
        s = join(A, delim)
        push!(stringout, s)
    end
    stringout
end

function topstrings{S<:AbstractString}(stringout::Vector{S}, n_top::Int64)
    cmap = countmap(stringout)
    cmap1 = collect(cmap)
    tmp = sort(collect(cmap1), by=kv->kv[2], rev=true)
    n_top = min(n_top, length(tmp)) #don't exceed length
    topentries = tmp[1:n_top] #take top n
    topentries = map(x->(x...), topentries)
    topentries
end

"""
Executes each blk individually in a loop.  Good for pin pointing exceptions when debugging.
"""
function exec_blks(sess::Session, ckt::Circuit, fd::FeedDict=FeedDict();
    verbose::Bool=true)
    results = Any[]
    for i = 1:endof(ckt.blks)
        verbose && println("blk=$(ckt.blks[i])")
        result = run(sess, out(ckt.blks[i]), fd)
        push!(results, result)
    end
    results
end

function gradient_tensor(target::Tensor, vars::Variable...)
    Tensor(gradients(target, Tensor([vars...]))[1])
end

function softselect_by_example(sess::Session, ckt::Circuit, fd::FeedDict=FeedDict())
    softsel = run(sess, softselect_tensor(ckt), fd)
    byexample = Any[]
    n_examples = size(softsel[1], 1) #rows of the first blk
    for i = 1:n_examples
        A = [softsel[j][i,:] for j=1:length(ckt.blks)]
        push!(byexample, A)
    end
    byexample
end


text:
## --- Auto-formatted `printf` methods for Julia types

# Pick a printf format string depending on the type
@inline printfmt(::Type{<:AbstractFloat}) = c"%e"
@inline printfmt(::Type{<:Integer}) = c"%d"
@inline printfmt(::Type{<:Ptr})  = c"Ptr @0x%016x" # Assume 64-bit pointers
@inline printfmt(::Type{UInt64}) = c"0x%016x"
@inline printfmt(::Type{UInt32}) = c"0x%08x"
@inline printfmt(::Type{UInt16}) = c"0x%04x"
@inline printfmt(::Type{UInt8})  = c"0x%02x"
@inline printfmt(::Type{<:Union{MallocString, StaticString}}) = c"\"%s\"" # Can I offer you a string in this trying time?

# Top-level formats, single numbers
"""
```julia
printf([fp::Ptr{FILE}], [fmt], n::Number)
```
Libc `printf` function, accessed by direct `llvmcall`.

Prints a number `n` to a filestream specified by the file pointer `fp`,
defaulting to the current standard output `stdout` if not specified.

Optionally, a C-style format specifier string `fmt` may be provided as well.

Returns `0` on success.

## Examples
```julia
julia> printf(1)
1
0

julia> printf(1/3)
3.333333e-01
0

julia> printf(c"%f\n", 1/3)
0.333333
0
```
"""
@inline function printf(n::T) where T <: Union{Number, Ptr}
    printf(printfmt(T), n)
    newline()
end

# Print a vector
"""
```julia
printf([fp::Ptr{FILE}], a::AbstractArray{<:Number})
```
Print a matrix or vector of numbers `a` to a filestream specified by the file
pointer `fp`, defaulting to the current standard output `stdout` if not specified.

Returns `0` on success.

## Examples
```julia
julia> printf(rand(5,5))
5.500186e-02    8.425572e-01    3.871220e-01    5.442254e-01    5.990694e-02
5.848425e-01    6.714915e-01    5.616896e-01    6.668248e-01    2.643873e-01
9.156712e-01    1.276033e-01    3.350369e-01    6.513146e-01    9.999104e-01
3.301038e-01    6.027120e-01    5.139433e-01    2.219796e-01    4.057417e-01
2.821340e-01    9.258760e-01    7.950481e-01    1.152236e-01    7.949463e-01
0
```
"""
@inline function printf(v::AbstractVector{T}) where T <: Union{Number, Ptr, StaticString}
    fmt = printfmt(T)
    p = pointer(fmt)
    @inbounds GC.@preserve fmt for i ∈ eachindex(v)
        printf(p, v[i])
        newline()
    end
    return zero(Int32)
end

# Print a tuple
@inline function printf(v::NTuple{N, T} where N) where T <: Union{Number, Ptr, StaticString}
    fmt = printfmt(T)
    p = pointer(fmt)
    putchar(0x28) # open paren
    @inbounds GC.@preserve fmt for i ∈ eachindex(v)
        printf(p, v[i])
        putchar(0x2c) # comma
        putchar(0x20) # space
    end
    putchar(0x29) # close paren
    newline()
end

# Print a 2d matrix
@inline function printf(m::AbstractMatrix{T}) where T <: Union{Number, Ptr, StaticString}
    fmt = printfmt(T)
    p = pointer(fmt)
    @inbounds GC.@preserve fmt for i ∈ axes(m,1)
        for j ∈ axes(m,2)
            printf(p, m[i,j])
            putchar(0x09) # tab
        end
        newline()
    end
    return zero(Int32)
end

## --- Printing to file

# Top-level formats, single numbers
@inline function printf(fp::Ptr{FILE}, n::T) where T <: Union{Number, Ptr}
    printf(fp, printfmt(T), n)
    newline(fp)
end

# Print a vector
@inline function printf(fp::Ptr{FILE}, v::AbstractVector{T}) where T <: Union{Number, Ptr, StaticString}
    fmt = printfmt(T)
    p = pointer(fmt)
    @inbounds GC.@preserve fmt for i ∈ eachindex(v)
        printf(fp, p, v[i])
        newline(fp)
    end
    return zero(Int32)
end

# Print a tuple
@inline function printf(fp::Ptr{FILE}, v::NTuple{N, T} where N) where T <: Union{Number, Ptr, StaticString}
    fmt = printfmt(T)
    p = pointer(fmt)
    putchar(fp, 0x28) # open paren
    @inbounds GC.@preserve fmt for i ∈ eachindex(v)
        printf(fp, p, v[i])
        putchar(fp, 0x2c) # comma
        putchar(fp, 0x20) # space
    end
    putchar(fp, 0x29) # close paren
    newline(fp)
end

# Print a 2d matrix
@inline function printf(fp::Ptr{FILE}, m::AbstractMatrix{T}) where T <: Union{Number, Ptr, StaticString}
    fmt = printfmt(T)
    p = pointer(fmt)
    @inbounds GC.@preserve fmt for i ∈ axes(m,1)
        for j ∈ axes(m,2)
            printf(fp, p, m[i,j])
            putchar(fp, 0x09) # tab
        end
        newline(fp)
    end
    return zero(Int32)
end

## -- Print errors

@inline Base.error(s::Union{StaticString,MallocString}) = (perror(c"ERROR: "); perror(s))
@inline warn(s::Union{StaticString,MallocString}) = (perror(c"Warning: "); perror(s))

"""
```julia
perror(s)
```
Print the string `s` to the standard error filestream, `stderr`.

Returns `0` on success.

## Examples
```julia
julia> StaticTools.perror(c"ERROR: could not do thing\n")
ERROR: could not do thing
0
```
"""
@inline perror(s::MallocString) = perror(pointer(s))
@inline perror(s::StaticString) = GC.@preserve s perror(pointer(s))
@inline function perror(s::Ptr{UInt8})
    Base.llvmcall(("""
    ; External declaration of the fprintf function
    declare i32 @fprintf(i8*, i8*)

    define i32 @main(i64 %jlfp, i64 %jls) #0 {
    entry:
      %fp = inttoptr i64 %jlfp to i8*
      %str = inttoptr i64 %jls to i8*
      %status = call i32 (i8*, i8*) @fprintf(i8* %fp, i8* %str)
      ret i32 0
    }

    attributes #0 = { alwaysinline nounwind ssp uwtable }
    """, "main"), Int32, Tuple{Ptr{FILE}, Ptr{UInt8}}, stderrp(), s)
end


text:
"""
```
createInitialValueDistribution( attribute::Attribute )
```
This function initialises the `initValueDist` field of the personnel attribute `attribute`, using the information int he other fields of the attribute object.

This function returns `nothing`.
"""
function createInitialValueDistribution( attribute::Attribute )

    probs = attribute.initValueWeights / sum( attribute.initValueWeights )
    attribute.initValueDist = Categorical( probs )
    return

end  # createInitialValueDistribution( attribute )

"""
```
isAttributeValuePossible(
    attribute::Attribute,
    value::String )
```
This function checks if the attribute `attribute` van have the value `value`.

This function returns a `Bool`, the result of the check.
"""
isAttributeValuePossible( attribute::Attribute, value::String ) =
    value ∈ attribute.possibleValues


text:
# TODO: prevent motion outside of domain

abstract type MotionModel end

struct NoMotion <: MotionModel
end
move_target(sm::NoMotion, theta, L::Float64) = theta

struct RandomMotion <: MotionModel
    d::Float64
end
RandomMotion() = RandomMotion(5.0)

function move_target(rm::RandomMotion, theta, L::Float64)
    ang = rand() * 2.0 * pi
    new_x = theta[1] + cos(ang) * rm.d
    new_y = theta[2] + sin(ang) * rm.d
    return theta[1] + rm.d*rand(), theta[2] + rm.d*rand(), 0.0, 0.0
end

struct ConstantMotion <: MotionModel
    noise_vals::NTuple{4,Float64}   # std deviations for each dimension
end
ConstantMotion() = ConstantMotion((1.0,1.0,0.05,0.05))
#ConstantMotion(a::Float64, b::Float = ConstantMotion((a,b,c,d))

function move_target(cm::ConstantMotion, theta, L::Float64)
    #new_x = theta[1] + theta[3]# + cm.noise_vals[1] * randn()
    #new_y = theta[2] + theta[4]# + cm.noise_vals[2] * randn()
    #new_xdot = theta[3]# + cm.noise_vals[3] * randn()
    #new_ydot = theta[4]## + cm.noise_vals[4] * randn()

    new_x = theta[1] + theta[3] + cm.noise_vals[1] * randn()
    new_y = theta[2] + theta[4] + cm.noise_vals[2] * randn()
    new_xdot = theta[3] + cm.noise_vals[3] * randn()
    new_ydot = theta[4] + cm.noise_vals[4] * randn()

    return new_x, new_y, new_xdot, new_ydot
end


text:
T = Float64
fs = 100. # sampling rate
N = 90000 # number of points
expect_output_type = Vector{Complex{Float64}}

# Test compute_fft function
@testset "compute fft" begin
    data = rand(T, N)
    fft_result = compute_fft(data)
    @test size(fft_result, 1) == N
    @test typeof(fft_result) == expect_output_type
end


text:
abstract type CompositeMean <: Mean end

components(m::CompositeMean) = m.means

@deprecate submeans components

include("sum_mean.jl")       # Sum mean functions
include("prod_mean.jl")      # Product of mean functions


text:
using Printf
using DataFrames
include(joinpath("src", "data_loader.jl"))
include(joinpath("src", "model.jl"))

df = data_loader.read_csv("data_proj_414.csv")

train_df, test_df = data_loader.data_split(df)
#fruit_number, centers, esigma = model.fit(df)
fruit_number, centers, esigma = model.fit(train_df)
@printf("fruit number:\t%d\n", floor(fruit_number))
@printf("Sample tree number:\t%d\n", size(centers)[2])
@printf("Total tree number:\t%d\n", size(centers)[2]*3)
loc_arr = [train_df.X train_df.Y]'
train_pred_fruit = model.predict(loc_arr, fruit_number, centers, esigma)
train_mse = model.mse(train_pred_fruit, train_df.Close)
train_acc_soft = model.accuarcy_soft(train_pred_fruit, train_df.Close)*100
@printf("train mse:\t%2.4f\n", train_mse)
@printf("train soft accuarcy:\t%2.2f%%\n", train_acc_soft)
println()
loc_arr = [test_df.X test_df.Y]'
test_pred_fruit = model.predict(loc_arr, fruit_number, centers, esigma)
test_mse = model.mse(test_pred_fruit, test_df.Close)
test_acc_soft = model.accuarcy_soft(test_pred_fruit, test_df.Close)*100
@printf("test mse:\t%2.4f\n", test_mse)
@printf("test soft accuarcy:\t%2.2f%%\n", test_acc_soft)
train_df[:predict_close] =  train_pred_fruit
train_df[:diff_close] = map((x, y) -> abs((y-x)),
train_df[:, :predict_close], train_df[:, :Close])
data_loader.save_csv(train_df, "predict.csv")
center_df = DataFrame(X = centers[1, :], Y = centers[2, :])
data_loader.save_csv(center_df, "centers.csv")

# println(summary(df))
# println(summary(train_df))
# println(summary(test_df)

# main()


text:
function puzzle09(path::String=joinpath(@__DIR__, "..", "data/09.txt"))
    stream = readline(path)

    idx = 1

    total = 0
    noncanceled = 0

    arr = []
    groups = 0
    isgarbage = false

    while idx <= length(stream)
        if stream[idx] == '!'
            idx += 1
        elseif isgarbage
            stream[idx] == '>' ? isgarbage = false : noncanceled += 1
        elseif stream[idx] == '{'
            groups += 1
            append!(arr, groups)
        elseif stream[idx] == '<'
            isgarbage = true
        elseif stream[idx] == '}'
            groups -= 1
            total += pop!(arr)
        end
        idx += 1
    end

    return [total, noncanceled]

end


text:
function constantdiffusionproblem(stepper; nx=128, Lx=2π, kappa=1e-2, nsteps=1000, dev=CPU())
   τ = 1/kappa  # time-scale for diffusive decay
  dt = 1e-9 * τ # dynamics are resolved

  prob = Problem(nx=nx, Lx=Lx, kappa=kappa, dt=dt, stepper=stepper, dev=dev)
  g = prob.grid

  # a gaussian initial condition c(x, t=0)
  c0ampl, σ = 0.01, 0.2
  c0func(x) = @. c0ampl*exp(-x^2/(2σ^2))
  c0 = c0func.(g.x)

  # analytic solution for for 1D heat equation with constant κ
  tfinal = nsteps*dt
  σt = sqrt(2*kappa*tfinal + σ^2)
  cfinal = @. c0ampl*σ/σt * exp(-g.x^2/(2*σt^2))

  set_c!(prob, c0)
  tcomp = @elapsed stepforward!(prob, nsteps)
  updatevars!(prob)

  prob, c0, cfinal, nsteps, tcomp
end

function varyingdiffusionproblem(stepper; nx=128, Lx=2π, kappa=1e-2, nsteps=1000, dev=CPU())
   τ = 1/kappa  # time-scale for diffusive decay
  dt = 1e-9 * τ # dynamics are resolved

  kappa = kappa*ones(nx) # this is actually a constant diffusion but defining it
                         # as an array makes stepforward! call function calcN!
                         # instead of just the linear coefficients L*sol

  prob = Problem(nx=nx, Lx=Lx, kappa=kappa, dt=dt, stepper=stepper, dev=dev)
  g = prob.grid

  # a gaussian initial condition c(x, t=0)
  c0ampl, σ = 0.01, 0.2
  c0func(x) = @. c0ampl*exp(-x^2/(2σ^2))
  c0 = c0func.(g.x)

  # analytic solution for for 1D heat equation with constant κ
  tfinal = nsteps*dt
  σt = sqrt(2*kappa[1]*tfinal + σ^2)
  cfinal = @. c0ampl*σ/σt * exp(-g.x^2/(2*σt^2))

  set_c!(prob, c0)
  tcomp = @elapsed stepforward!(prob, nsteps)
  updatevars!(prob)

  prob, c0, cfinal, nsteps, tcomp
end

function constantdiffusiontest(stepper, dev::Device=CPU(); kwargs...)
  prob, c0, c1, nsteps, tcomp = constantdiffusionproblem(stepper; kwargs...)
  normmsg = "$stepper: relative error ="
  @printf("% 40s %.2e (%.3f s)\n", normmsg, norm(c1-prob.vars.c)/norm(c1), tcomp)
  isapprox(c1, prob.vars.c, rtol=nsteps*rtol_timesteppers)
end

function varyingdiffusiontest(stepper, dev::Device=CPU(); kwargs...)
  prob, c0, c1, nsteps, tcomp = varyingdiffusionproblem(stepper; kwargs...)
  normmsg = "$stepper: relative error ="
  @printf("% 40s %.2e (%.3f s)\n", normmsg, norm(c1-prob.vars.c)/norm(c1), tcomp)
  isapprox(c1, prob.vars.c, rtol=nsteps*rtol_timesteppers)
end


text:
module EGR

export MinusPlusOneVector, OutputOpts, alg, getStats, readBin, createBLOracles, L2regGradient, readData, writeBin, SGsd, Opts, VerifyGradient, VerifyRestoration, findBestStepsizeFactor, createMLOracles, StepData, writeFunction, Outputter, ResultFromOO, Problem, returnIfExists, EGRexp, createSBLOracles, writeFinal, getF, getPCC, getMCC, normalizeFeatures, trainTestRandomSeparate, createClassLabels, getStatsHIGGSSUSY, readDataHIGGSSUSY, getStatsIRIS, readDataIRIS, readDataSparse, StrainTestRandomSeparate,DSSsd,DSSexp, getSequentialFinite, getRandom, SAGinit, getSampleFunctionAt, SAGAinit, SSVRG,SSVRGsd,SSVRG,SSVRGfromEGR,SAG,SAGA,getFfinal,SAGcb,SAGAcb, DSS,getSequential,onlyAddBeta1,onlyUpdateBeta1,uLinBeta1,uQuadBeta1,uExpBeta1,USD,SGksd,SGsqrtsd,getFtrainfinal,getFtrain,UcbonlyAddBeta1,UcbonlyUpdateBeta1,UcbuLinBeta1,UcbuQuadBeta1,UcbuExpBeta1,getFoverTen,getGatWhichFisAchieved,urLinBeta1,urQuadBeta1,urExpBeta1,getSequentialR,onlyrAddBeta1,urExpBeta0,getSequentialRFinite,uExp2Beta1,urExp2Beta1

include("MinusPlusOneVector.jl")
include("DataPrep.jl")
include("LossFunctions/lossBL.jl")
include("LossFunctions/lossSBL.jl")
include("LossFunctions/lossML.jl")
include("myIO.jl")
include("CreateOracles.jl")
include("LossFunctions/L2Reg.jl")
include("Alg.jl")
include("Steps/stepSG.jl")
include("Steps/stepU.jl")
include("Steps/stepUrestore.jl")
include("Steps/stepUcb.jl")
include("Steps/stepDSS.jl")
include("Steps/stepSAGinit.jl")
include("Steps/stepSAGAinit.jl")
include("Steps/stepSSVRG.jl")
include("Steps/stepSG.jl")
include("Steps/stepSGk.jl")
include("Steps/stepSGsqrt.jl")
include("Steps/stepSAG.jl")
include("Steps/stepSAGA.jl")
include("Steps/stepSAGcb.jl")
include("Steps/stepSAGAcb.jl")
include("FindBest.jl")
include("Producers.jl")
include("RedisStuff.jl")

end


text:
using Dash, DashBootstrapComponents

app = dash(external_stylesheets=[dbc_themes.GRID]);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
type Dimension
    size::Int64
    regions
    types
end

function dim_rand_sample(dim::Dimension)
    x = []
    for i = 1:dim.size
        value = 0
        if dim.types[i] == true
            value = rand(rng, Float64) * (dim.regions[i][2] - dim.regions[i][1]) + dim.regions[i][1]
        else
            value = rand(rng, dim.regions[i][1]:dim.regions[i][2])
        end
        append!(x, value)
    end
    return x
end

function dim_limited_space(dim::Dimension)
    number = 1
    for i = 1:dim.size
        if dim.types[i] == true
            return false, 0
        else
            number = number * (dim.regions[i][2] - dim.regions[i][1] + 1)
        end
    end
    return true, number
end

function dim_print(dim::Dimension)
    zoolog("dim size $(dim.size)")
    zoolog("dim regions is: ")
    zoolog("$(dim.regions)")
    zoolog("dim types is: ")
    zoolog("$(dim.types)")
end

function is_discrete(dim::Dimension)
    for i in 1:length(dim.types)
        if dim.types[i] == true
            return false
        end
    end
    return true
end