text:
import Dolo
import Distributions

# test IID-Normal
mu = [0.0, 0.0]
Sigma = [1.0 0.0; 0.0 2.0]
d = Dolo.MvNormal(mu, Sigma)
sim = Dolo.simulate(d, 1, 200, [0.1, -0.01])
sim = Dolo.simulate(d, 1, 200, [0.1, -0.01]; stochastic=false)
dp = Dolo.discretize(d)

# test VAR
var = Dolo.VAR1([0.0,0.0],[0.99 0.0; 0.0 0.08],Matrix(I,2,2)*0.01)

@time sim = Dolo.simulate(var, 500, 500)

Tvalues = linspace(1, sim[Axis{:T}][end], sim[Axis{:T}][end])
import PyPlot
plt = PyPlot;
fig = plt.figure("Markov Chaine")
PyPlot.plot(Tvalues, sim[1,1,:])

sim = Dolo.simulate(var, 5, 200, [0.8, 0.8])
sim = Dolo.simulate(var, 1, 200, [0.8, 0.8]; stochastic=false)

dp = Dolo.discretize(var)

irf = Dolo.response(var, [0.1, 0.3])
irf = Dolo.response(var, [0.8, 0.8], [0.1, 0.3])
irf = Dolo.response(var, [0.8, 0.8], 2)
irf = Dolo.response(var)

import temp

# markov chain
values = [ 0   0;
           0.1 0;
           0.1 0.1 ]

transitions = [ 0.1 0.4 0.5;
                0.1 0.8 0.1;
                0.2 0.0 0.8 ]

mc = DiscreteMarkovProcess(transitions, values)
dproc_mc = discretize(mc)
dproc_mc.nodes
dproc_mc.integration_weights
dproc_mc.integration_weights
@assert n_nodes(dproc_mc) == 3
dproc_mc

# normal law
Sigma = [0.05 0.01; 0.01 0.1]
mvn = MvNormal(Sigma)
dproc_mvn = discretize(mvn)
dproc_mvn.integration_weights
dproc_mvn.integration_nodes
dproc_mvn
@assert n_nodes(dproc_mvn) == 0

###############################################################################
M = rand(2)
Sigma = [0.5 0.5; 0.5 0.5]
R = [0.9 0.1; 0 .8];
n_states= [5,10]
n_states::Array{Int,1}
n = 2
n_integration=[7,4]
# 4/ Create a VAR1 type (in the module). Add new functions that operate on this type. Add them to the test file.

VAR_info = VAR1(M,R,Sigma)
using Base.product
D_Var = discretize(VAR_info, n_states, n_integration)
D_Var.nodes
D_Var.integration_nodes
D_Var.integration_weights[2]
[sum(D_Var.integration_weights[j]) for j in 1:kron([n_integration[j] for j in 1:n]...)]

#

M = rand(1)
Sigma = [0.5]'
R = [0.9]'
n_states= [5]
n_states::Array{Int,1}
n = 2
n_integration=[7]
# 4/ Create a VAR1 type (in the module). Add new functions that operate on this type. Add them to the test file.

VAR_info = VAR1(M,R,Sigma)
using Base.product
D_Var = discretize(VAR_info, n_states, n_integration)
D_Var.nodes
D_Var.integration_nodes
D_Var.integration_weights[2]
[sum(D_Var.integration_weights[j]) for j in 1:kron([n_integration[j] for j in 1:n]...)]


text:
# Optimized Product Quantization. Adapted from Mohammad Norouzi's code.
export train_opq, quantize_opq

"""
    quantize_opq(X, R, C, V=false) -> B

Given data and PQ/OPQ codeboks, quantize.

# Arguments
- `X::Matrix{T}`: `d`-by-`n` data to quantize
- `R::Matrix{T}`: `d`-by-`d` rotation to apply to the data before quantizing
- `C::Vector{Matrix{T}}`: `m`-long vector with `d/m`-by-`h` matrix entries. Each matrix is a (O)PQ codebook.
- `V::Bool`: Whether to print progress

# Returns
- `B::Matrix{Int16}`: `m`-by-`n` matrix with the codes that approximate `X`
"""
function quantize_opq(
  X::Matrix{T},          # d-by-n matrix of data points to quantize
  R::Matrix{T},          # d-by-d matrix. Learned rotation for X
  C::Vector{Matrix{T}},  # m-long array. Each entry is a d-by-h codebooks
  V::Bool=false) where T <: AbstractFloat

  # Apply rotation and quantize as in PQ
  return quantize_pq(R'*X, C, V)
end

"""
    train_opq(X, m, h, niter, init, V=false) -> C, B, R, error

Trains an optimized product quantizer.

# Arguments
- `X::Matrix{T}`: `d`-by-`n` data to quantize
- `m::Integer`: Number of codebooks
- `h::Integer`: Number of entries in each codebook (typically 256)
- `niter::Integer`: Number of iterations to use
- `init::String`: Method used to intiialize `R`, either `"natural"` (identity) or `"random"`.
- `V::Bool`: Whether to print progress

# Returns
- `B::Matrix{Int16}`: `m`-by-`n` matrix with the codes
- `C::Vector{Matrix{T}}`: `m`-long vector with `d`-by-`h` matrix entries. Each matrix is a codebook of size approximately `d/m`-by-`h`.
- `R::Matrix{T}`: `d`-by-`d` learned rotation for the data
- `obj::Vector{T}`: The quantization error each iteration
"""
function train_opq(
  X::Matrix{T},      # d-by-n matrix of data points to train on.
  m::Integer,        # number of codebooks
  h::Integer,        # number of entries per codebook
  niter::Integer,    # number of optimization iterations
  init::String,      # how to initialize the optimization
  V::Bool=false) where T <: AbstractFloat

  if V; @printf("Training an optimized product quantizer\n"); end

  d, n = size(X)

  C = Vector{Matrix{T}}(undef, m) # codebooks

  obj = zeros(Float32, niter+1 )

  # Number of bits in the final codes.
  nbits = log2(h) * m
  CB    = zeros(T, size(X))

  if init == "natural" # Initialize R with identity
    R = Matrix{T}(1.0I, d, d)
  elseif init == "random"
    R, _, _ = svd( randn( T, d, d ))
  else
    error("Intialization $init unknown")
  end

  RX = R' * X # Rotate the data

  subdims = splitarray( 1:d, m )

  # Initialization sampling RX
  for i = 1:m
    perm = sample(1:n, h, replace=false)
    C[i] = RX[ subdims[i], perm ]
  end

  # Variables needed for methods in Clustering.kmeans
  costs      = zeros(Float32, n)
  counts     = zeros(Int, h)
  cweights   = zeros(Float64, h)
  to_update  = zeros(Bool, h)
  to_update2 = ones(Bool, h)
  unused     = Int[]

  # Initialize the codes -- B
  B    = Vector{Vector{Int}}(undef, m) # codes
  for i = 1:m; B[i] = zeros(Int, n); end # Allocate codes

  for i=1:m
    dmat = Distances.pairwise(Distances.SqEuclidean(), C[i], RX[subdims[i],:])
    Clustering.update_assignments!(dmat, true, B[i], costs, counts, to_update, unused)
    CB[subdims[i],:] .= C[i][:, B[i]]
  end

  for iter=0:niter
    if V; start_time = time_ns(); end # Take time if asked to

    # Compute objective function
    obj[iter+1] = sum( (R*CB - X).^2 ) ./ n
    if V; @printf("%3d %e... ", iter, obj[iter+1]); end

    # update R
    U, S, VV = svd(X * CB', full=false)
    R = U * VV'

    # update R*X
    RX = R' * X

    for i=1:m
      # update C
      Clustering.update_centers!(RX[subdims[i],:], nothing, B[i], to_update2, C[i], cweights)

      # update B
      dmat = Distances.pairwise(Distances.SqEuclidean(), C[i], RX[subdims[i],:])
      Clustering.update_assignments!(dmat, false, B[i], costs, counts, to_update, unused)

      # update D*B
      CB[subdims[i], :] .= C[i][:, B[i]]
    end # for i=1:m

    if V; @printf("done in %.2f secs\n", (time_ns() - start_time)/1e9); end
  end # for iter=0:niter

  B = hcat(B...)
  B = convert(Matrix{Int16}, B)
  C = convert(Vector{Matrix{T}}, C)

  return C, collect(B'), R, obj
end

function experiment_opq(
  Xt::Matrix{T}, # d-by-n. Data to learn codebooks from
  Xb::Matrix{T}, # d-by-n. Base set
  Xq::Matrix{T}, # d-by-n. Queries
  gt::Vector{UInt32}, # ground truth
  m::Integer,    # number of codebooks
  h::Integer,    # number of entries per codebook
  init::String="natural", # initialization method for the rotation
  niter::Integer=25, # Number of k-means iterations for training
  knn::Integer=1000,
  V::Bool=false) where T <: AbstractFloat # whether to print progress

  # === Train ===
  C, B, R, train_error = train_opq(Xt, m, h, niter, init, V)
  if V; @printf("Error in training is %e\n", train_error[end]); end

  # === Encode the base set ===
  B_base     = quantize_opq(Xb, R, C, V)
  base_error = qerror_opq(Xb, B_base, C, R)
  if V; @printf("Error in base is %e\n", base_error); end

  # === Compute recall ===
  if V; println("Querying m=$m ... "); end
  b = Int(log2(h) * m)
  @time dists, idx = linscan_opq(B_base, Xq, C, b, R, knn)
  if V; println("done"); end

  recall = eval_recall(gt, idx, knn)
  return C, B, R, train_error, B_base, recall
end

function experiment_opq_query_base(
  Xt::Matrix{T}, # d-by-n. Data to learn codebooks from
  Xq::Matrix{T}, # d-by-n. Queries
  gt::Vector{UInt32}, # ground truth
  m::Integer,    # number of codebooks
  h::Integer,    # number of entries per codebook
  init::String="natural", # initialization method for the rotation
  niter::Integer=25, # Number of k-means iterations for training
  knn::Integer=1000,
  V::Bool=false) where T <: AbstractFloat # whether to print progress

  # === Train ===
  C, B, R, train_error = train_opq(Xt, m, h, niter, init, V)
  if V; @printf("Error in training is %e\n", train_error[end]); end

  # === Compute recall ===
  if V; println("Querying m=$m ... "); end
  b = Int(log2(h) * m)
  @time dists, idx = linscan_opq(B, Xq, C, b, R, knn)
  if V; println("done"); end

  recall = eval_recall(gt, idx, knn)
  return C, B, R, train_error, recall
end

ext: jl


text:
export row, col, reshape
export split, concat, cat1, cat2, cat3, cat4, cat5
export dropdim1, dropdim2, dropdim3, dropdim4, dropdim5, squeeze
export subtoind, indtosub
export stack, flatten, unstack, vflatten, unflatten
export riffle
export matrix, unmatrix
export lines, unlines, unzip
export findsub, find_
export randsample, flip, flip!, flipdims

#######################################
## row, col, reshape

row(a::Number) = row([a])
row(a) = reshape(a,1,length(a))
row(a::Tuple) = reshape([a...],1,length(a))
row(a...) = row(vcat(a...))

col(a::Number) = col([a])
col(a) = reshape(a,length(a),1)
col(a::Tuple) = reshape([a...],length(a),1)
col(a...) = col(vcat(a...))

import Base.reshape
reshape(a::AbstractArray, siz::Array{T,2}) where {T} = reshape(a, siz...)
function reshape(a::AbstractArray)
    if !(ndims(a)==1 || len(a)==1)
        error("FuncionalData.reshape(a): a was expected to be a vector or len(a)==1, actual size was: $(size(a))")
    end
    r2 = sqrt(length(a))
    r3 = cbrt(length(a))
    if round(r2)==r2
        return Base.reshape(a, round(Int,r2), round(Int,r2))
    elseif round(r3)==r3
        return Base.reshape(a, round(Int,r3), round(Int,r3), round(Int,r3))
    else
        error("cannot reshape array of size $(size(a)) to a square or cube")
    end
end

#######################################
##  split, cat

import Base.split
split(a::AbstractArray,x) = split(a, y->y==x)

function split(a::AbstractArray,f::Function)
    start = 1
    r = Any[]
    for i = 1:length(a)
        if f(a[i])
            push!(r, a[start:i-1])
            start = i + 1
        elseif i == length(a)
            push!(r, a[start:i])
        end
    end
    return r
end

concat(a) = concat(a...)
concat(a...) = @p flatten [reject(collect(a),x->(try return isempty(x) catch; return false end))...]
cat1(a) = cat(a..., dims=1)
cat2(a) = cat(a..., dims=2)
cat3(a) = cat(a..., dims=3)
cat4(a) = cat(a..., dims=4)
cat5(a) = cat(a..., dims=5)
cat1(a, b, c...) = cat(a, b, c..., dims=1)
cat2(a, b, c...) = cat(a, b, c..., dims=2)
cat3(a, b, c...) = cat(a, b, c..., dims=3)
cat4(a, b, c...) = cat(a, b, c..., dims=4)
cat5(a, b, c...) = cat(a, b, c..., dims=5)

dropdim1(a) = dropdims(a, dims = 1)
dropdim2(a) = dropdims(a, dims = 2)
dropdim3(a) = dropdims(a, dims = 3)
dropdim4(a) = dropdims(a, dims = 4)
dropdim5(a) = dropdims(a, dims = 5)

function squeeze(a)
    t = tuple(findall(size(a).==1)...)
    t == () ? a : dropdims(a, dims = t)
end

#######################################
## subtoind, indtosub

function subtoind(subs, a)
    strides_ = strides(a)
    r = oneel(subs)
    for i in 1:length(strides_)
        r += (subs[i]-1) * strides_[i]
    end
    return r
end

indtosub(i::Int, a) = indtosub([i], a)
function indtosub(inds::AbstractArray, a::Union{Array,BitArray})
    if length(a) <= 5
        s = tuple(a...)
    else
        s = size(a)
    end
    @p map CartesianIndices(s)[inds] (x->col([convert(Tuple,x)...])) | flatten
end

#######################################
## flatten, stack, unstack, unflatten

function stack(a::Array{T,1}) where {T}
    r = arraylike(fst(a), len(a))
    for i = 1:len(a)
        setat!(r, i, at(a,i))
    end
    return r
end
stack(a::DenseArray{T}) where {T<:Number} = a

const StringLike  = Union{Char, AbstractString}
tostring(a) = string(a)
tostring(a::AbstractString) = a
flatten(a::StringLike) = a
flatten(a::Array{T,1}) where {T<:StringLike} = join(a)

function flatten(a::Array{T,2}) where {T<:StringLike}
    @assert size(a,1)==1
    join(vec(a))
end

flatten(a::AbstractArray{T,N}) where {T<:Number,N} = a
function flatten(a::Array{T,1}) where {T}
    if isempty(a)
        return similar(a)
    end
    if isa(a[1], StringLike)
        return join(map(a,tostring))
    end
    isa(a[1], AbstractArray) || return a
    if ndims(fst(a)) == 1
        return vcat(a...)
    end
    firstnonempty = @p map a not*isempty | findfirst
    r = arraylike(fst(at(a,firstnonempty)), sum([len(x) for x in a]))
    ind = 0
    for i = 1:len(a)
        for j = 1:len(a[i])
            setat!(r, ind+j, at(a[i],j))
        end
        ind += len(a[i])
    end
    return r
end

function flatten(a::Array{T,2}) where {T}
    if isempty(a)
        return arraylike(a)
    end
    if isa(a[1], StringLike)
        @assert size(a,1) == 1
        return flatten(vec(a))
    end
    isa(a[1], AbstractArray) || return a
    ms = Base.map(x->size(x,1), a)
    ns = Base.map(x->size(x,2), a)
    msum = 0
    nsum = 0
    for mi = 1:size(a,1)
        for ni = 2:size(a,2)
            if ms[mi,ni] != ms[mi,1]
                error("flatten: ms[mi,ni] was $(ms[mi,ni]) but ms[mi,1] was $(ms[mi,1])")
            end
        end
        msum += ms[mi,1]
    end
    for ni = 1:size(a,2)
        for mi = 2:size(a,1)
            if ns[mi,ni] != ns[1,ni]
                error("flatten: ns[mi,ni] was $(ns[mi,ni]) but ns[1,ni] was $(ns[1,ni])")
            end
        end
        nsum += ns[1,ni]
    end

    typ = eltype(a[1])
    for i = 1:length(a)
        typ = promote_type(typ, eltype(a[i]))
    end

    r = Array{typ,2}(undef, msum, nsum)

    ncum = 0
    for n = 1:size(a,2)
        mcum = 0
        for m = 1:size(a,1)
            r[mcum .+ (1:ms[m,n]), ncum .+ (1:ns[m,n])] .= a[m,n]
            mcum += ms[m,n]
        end
        ncum += ns[1,n]
    end
    r
end

flatten(a::Array{T,3}) where {T} = @p map a flatten | stack

vflatten(a) = @p permutedims a | flatten | permutedims

unstack(a) = Any[at(a,i) for i in 1:len(a)]

function unflatten(a,template)
    lens = @p mapvec template len
    ends = cumsum(lens)
    starts = @p droplast [1; ends.+1]
    @p mapvec2 starts ends (i,j)->part(a,i:j)
end

#######################################
##  riffle

function riffle(a::AbstractString,x::Union{Char,AbstractString})
    join(a,x)
end

function riffle(a::Array{Any,1},x)
    if len(a)==1 return a end

    dims = size(a)
    #dims[end] = dims[end]*2-1
    dims = tuple(dims[1:end-1]..., dims[end]*2-1)
    r = similar(a,dims...)

    for i = 1:len(a)
        r[i*2-1] = a[i]
    end
    for i = 2:2:len(r)-1
        r[i] = x
    end
    return r
end

function riffle(a,x)
    if len(a)==1 return a end

    dims = size(a)
    #dims[end] = dims[end]*2-1
    dims = tuple(dims[1:end-1]..., dims[end]*2-1)
    # r = similar(a,dims...)
    r = Array{Any}(undef, dims...)
    # println("dims: $dims   r: $r")
    for i = 1:len(a)
        setat!(r, i*2-1, at(a,i))
    end
    for i = 2:2:len(r)-1
        setat!(r, i, x)
    end
    return r
end

#######################################
##  matrix, unmatrix

matrix(a::AbstractArray{T}) where {T<:Number} = @p reshape a div(length(a),len(a)) len(a)
matrix(a) = @p map a col | flatten
function unmatrix(a, example)
    r = @p map a (x->reshape(x,siz(fst(example))))
    if eltype(example)==Any
        r = unstack(r)
    end
    r
end

#######################################
##  lines, unlines

function lines(a)
    r = split(a,r"\r\n|\n")
    if !isempty(r) && last(r)==""
        droplast(r)
    else
        r
    end
end
unlines(a) = join(a,'\n')

#######################################
##  randperm, randsample

import Random.randperm
_randperm = randperm
randperm(a::Number) = _randperm(a)
randperm(a) = part(a, randperm(len(a)))

randsample(a, n = 1) = part(a, rand(1:len(a), round(Int,n)))

flip!(a::Array) = a[:] = flip(a)
flip(a::Pair) = Pair(snd(a),fst(a))
flip(a::Dict) = @p vec a | map flip | Dict
flip(a) = part(a, len(a):-1:1)
function flipdims(a,d1,d2)
    dims = collect(1:ndims(a))
    dims[d1] = d2
    dims[d2] = d1
    permutedims(a, dims)
end

#####################################################
##   unzip

function unzip(a)
    forth(a) = at(a,4)
    fifth(a) = at(a,5)
    if len(fst(a)) == 2
        map(a, fst), map(a, snd)
    elseif len(fst(a)) == 3
        map(a, fst), map(a, snd), map(a, third)
    elseif len(fst(a)) == 4
        map(a, fst), map(a, snd), map(a, third), map(a, forth)
    elseif len(fst(a)) == 5
        map(a, fst), map(a, snd), map(a, third), map(a, forth), map(a, fifth)
    else
        error("FunctionalData.unzip: cannot unzip of len(fst(a))==$(len(fst(a)))")
    end
end

#####################################################
##   find

findsub(a::Array) = findsub(a.!=0)
findsub(a::BitArray) = indtosub(findall(vec(a)), a)

find_(a) = (LinearIndices(a))[findall(a)]


text:
using Base.Test
using LyapunovExponents
using LyapunovExponents: stable_le_error, FixedPointConvDetail,
    NonNegativeAutoCovConvDetail, NonPeriodicConvDetail,
    report, dimension, arnold_cat_map, van_der_pol

@time @testset "stable_le_error" begin
    with_negative(xs) = [xs, -xs]

    for xs in with_negative(repeat(1:3, outer=10))
        err, detail = stable_le_error(xs)
        @test detail.period == 3
    end

    for xs in with_negative(repeat([1], outer=10))
        err, detail = stable_le_error(xs)
        @test detail isa FixedPointConvDetail
        @test detail.var == 0
        @test err == 0
    end

    for xs in with_negative(collect(1:20))
        err, detail = stable_le_error(xs, cutoff=5)
        @test detail isa NonNegativeAutoCovConvDetail
        @test err == Inf
    end

    for xs in with_negative(collect(1:20))
        err, detail = stable_le_error(xs)
        @test detail isa NonPeriodicConvDetail
        @test err == Inf
    end
end

@time @testset "Terminate stable system $(ex.name)" for ex in [
        arnold_cat_map(M = [0.5 0.1; 1 0.1],
                       t_attr = 100000).example,
        van_der_pol(params = [:a => 3],
                    t_attr = 100000).example,
        van_der_pol(params = [:a => 0],
                    t_attr = 100000).example,
        van_der_pol(params = [:a => 0.1],
                    t_attr = 100000).example,
        ]
    for dim_lyap in 1:dimension(ex)
        println("$(ex.name) dim_lyap=$dim_lyap")
        @time sol = solve(ex; dim_lyap=dim_lyap)
        report(sol)
        @test lyapunov_exponents(sol)[1] < 0
        @test sol.converged
    end
end


text:
# This file is part of the TaylorSeries.jl Julia package, MIT license
#
# Parameters for HomogeneousPolynomial and TaylorN

@doc """
    ParamsTaylorN

DataType holding the current parameters for `TaylorN` and
`HomogeneousPolynomial`.

**Fields:**

- `order          :: Int`  Order (degree) of the polynomials
- `num_vars       :: Int`  Number of variables
- `variable_names :: Array{String,1}` Name of the variables

These parameters can be changed using `set_params_TaylorN(order, numVars)`.
""" ->
type ParamsTaylorN
    order          :: Int
    num_vars       :: Int
    variable_names :: Array{String,1}
end

@compat const _params_TaylorN_ = ParamsTaylorN(6, 2, String["x₁", "x₂"])

## Utilities to get the maximum order and number of variables
get_order() = _params_TaylorN_.order
get_numvars() = _params_TaylorN_.num_vars
get_variable_names() = _params_TaylorN_.variable_names

set_variable_names{T<:AbstractString}(names::Vector{T}) = _params_TaylorN_.variable_names = names

get_variables() = [TaylorN(i) for i in 1:get_numvars()]

@doc """
    set_variables(R, names; [order=6])
    set_variables(names; [order=6])
    set_variables(R, names; [order=6, numvars=-1])
    set_variables(names; [order=6, numvars=-1])

Return a `TaylorN{R}` vector with `numvars` polynomials, each representing an
independent variable, using `names` as the `String` for the output.

If `numvars` is not specified, it is inferred from the length of `names`.
If `length(names)==1` and `numvars>1`, it uses this name with subscripts for
the different variables. When changing the `order` or `numvars`, the
hash_tables are reset.
""" ->
function set_variables{T<:AbstractString}(R::Type, names::Vector{T}; order=6)
    order >= 1 || error("Order must be at least 1")

    num_vars = length(names)
    num_vars >= 1 || error("Number of variables must be at least 1")

    _params_TaylorN_.variable_names = names

    if !(order == get_order() && num_vars == get_numvars())
        # if these are unchanged, no need to regenerate tables

        _params_TaylorN_.order = order
        _params_TaylorN_.num_vars = num_vars

        resize!(coeff_table,order+1)
        resize!(index_table,order+1)
        resize!(size_table,order+1)
        resize!(pos_table,order+1)

        coeff_table[:], index_table[:], size_table[:], pos_table[:] = generate_tables(num_vars, order)
        gc();
    end

    # return a list of the new variables
    TaylorN{R}[TaylorN(R,i) for i in 1:get_numvars()]
end

set_variables{T}(names::Vector{T}; order=6) = set_variables(Float64, names, order=order)

function set_variables{T<:AbstractString}(R::Type, names::T; order=6, numvars=-1)
    variable_names = split(names)

    if length(variable_names) == 1 && numvars >= 1
        name = variable_names[1]
        variable_names = [string(name, subscriptify(i)) for i in 1:numvars]
    end

    set_variables(R, variable_names, order=order)
end

set_variables{T<:AbstractString}(names::T; order=6, numvars=-1) =
    set_variables(Float64, names, order=order, numvars=numvars)

@doc """
    show_params_TaylorN()

Display the current parameters for `TaylorN` and `HomogeneousPolynomial` types.
""" ->
function show_params_TaylorN()
    info( """
    Parameters for `TaylorN` and `HomogeneousPolynomial`:
    Maximum order       = $(get_order())
    Number of variables = $(get_numvars())
    Variable names      = $(get_variable_names())
    """)
    nothing
end


text:
"""
    intersect(chartA, chartB) -> [chart1, chart2, ...]

Compute the intersection of two charts of equal dimension.

Compute an array of charts such that the disjoint union of those charts produces the intersection of the two charts provided as inputs. In particular the sum of integrals over the returned charts will equal the integral over the intersection of the two given charts.
"""
function intersection(p1::Simplex{U,1,C,2,T}, p2::Simplex{U,1,C,2,T}) where {U,C,T}

    tol = sqrt(eps(T))
    PT = eltype(p1.vertices)

    W = [u for u in p2.vertices]
    clipconvex!(W, p1.vertices[1],  p1.tangents[1])
    clipconvex!(W, p1.vertices[2], -p1.tangents[1])

    # For consistency an array needs to be return. In higher
    # dimensions the intersection could be the union of
    # multiple simplices.
    [simplex(W[1], W[2])]
end

function clipconvex!(W, v, m)
    m2 = dot(m,m)
    for i in 1:length(W)
        s = dot(W[i]-v,m)
        s = min(0, s)
        W[i] = v + s / m2 * m
    end
end

"""
    intersection(triangleA, triangle B)

The output inherits the orientation from the first operand
"""
function intersection(p1::Simplex{U,2,C,3}, p2::Simplex{U,2,C,3}) where {U,C}
  pq = sutherlandhodgman(p1.vertices, p2.vertices)
  [ simplex(pq[1], pq[i], pq[i+1]) for i in 2:length(pq)-1 ]
end

function intersection(p1::Simplex{U,3,C,4}, p2::Simplex{U,3,C,4}) where {U,C}
  @assert overlap(p1,p2)
  volume(p1) <= volume(p2) ? [p1] : [p2]
end

"""
    intersectline(a,b,p,q)

Computes the intersection of the lines (in a 2D space) defined
by points [a,b] and [p,q]
"""
function intersectlines(a,b,p,q)

    P = typeof(a)

    d1 = det(@SMatrix [a[1] a[2]; b[1] b[2]])
    d2 = det(@SMatrix [p[1] p[2]; q[1] q[2]])

    id = one(eltype(a))
    d1x = det(@SMatrix [a[1] id; b[1] id])
    d1y = det(@SMatrix [a[2] id; b[2] id])

    d2x = det(@SMatrix [p[1] id; q[1] id])
    d2y = det(@SMatrix [p[2] id; q[2] id])

    den = det(@SMatrix [d1x d1y; d2x d2y])
    @assert !isinf(1/den)

    P(
        det(@SMatrix [d1 d1x; d2 d2x]) / den,
        det(@SMatrix [d1 d1y; d2 d2y]) / den,
    )

end

"""
    inside(p,a,b)

Determines is p is on the interior side of the segment b of the boundary,
assuming that the boundary is oriented counter-clockwise.
"""
function leftof(p,a,b)

    tol = sqrt(eps(eltype(a)))
    ap = @SVector [ p[1]-a[1], p[2]-a[2] ]
    ab = @SVector [ b[1]-a[1], b[2]-a[2] ]
    ap[1]*ab[2] - ap[2]*ab[1] <= tol ? true : false

end

#export sutherlandhodgman2d

"""
    sutherlandhodgman2d(subject,clipper)

Computes the intersection of the coplanar triangles
subject and clipper.
"""
function sutherlandhodgman2d(subject,clipper)

    PT = eltype(subject)

    clipped = Array(copy(subject))
    sizehint!(clipped, 8)

    input = Array(copy(clipped))
    sizehint!(input, 8)

    b = last(clipper)
    for a in clipper

        resize!(input, length(clipped))
        copyto!(input, clipped)
        resize!(clipped, 0)

        isempty(input) || (q = last(input))

        for p in input
            if leftof(p, b, a)
                if !leftof(q, b, a)
                    ist = intersectlines(b,a,q,p)
                    push!(clipped, ist)
                end

                push!(clipped, p)

            elseif leftof(q, b, a)

                ist = intersectlines(b,a,q,p)
                push!(clipped, ist)
            end

            q = p
        end
        b = a
    end

    return clipped
end

#export sutherlandhodgman

"""
    sutherlandhodgman(subject, clipper)

Compute the intersection of two coplanar triangles, potentially
embedded in a higher dimensional space.
"""
function sutherlandhodgman(subject, clipper)

    triangle = simplex(clipper, Val{2})
    subject2d = [carttobary(triangle,p) for p in subject]
    for p in subject2d for x in p @assert !isinf(x) end end
    clipper2d = [carttobary(triangle,q) for q in clipper]
    for p in clipper2d for x in p @assert !isinf(x) end end
    clipped2d = sutherlandhodgman2d(subject2d, clipper2d)
    clipped = [barytocart(triangle,q) for q in clipped2d ]

end


text:
immutable BInt{ R, T } <: Integer
        n::T
        function BInt( x::T , _ )
                new( x )
        end
end

# get BInt as a non-BInt integer in the simplest way:
int{R,T}( b::BInt{R,T} ) = b.n

Base.big{R,T}( b::BInt{R,T} ) = Base.big( b.n )

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

immutable UncheckedBInt{R,T}
end

@generated function Base.call{R,T,T2}( ::Type{ UncheckedBInt{R,T} }, x::T2 )
        # check that the given type of T is suitable
        # BigInts are always suitable, but don't have typemin,typemax
        if T == BigInt
                return :( BInt{R,T}( T(x), 0 ) )
        else
                if !( range_within_range( R, T ) )
                        throw( InexactError() );
                end
                return :( BInt{R,T}( x % T, 0 ) )
        end
end

@generated function Base.call{R,T2}( ::Type{ UncheckedBInt{R} }, x::T2 )
        T = default_int_type(R)
        if T == BigInt
                :( BInt{R,$T}( $T(x), 0 ) )
        else
                :( BInt{R,$T}( x % $T, 0 ) )
        end
end

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

# this is for constructing constants - where the range is
# just a single number (or a tuple in the case of very large numbers)
@generated function Base.call{R}( ::Type{ BInt{R} } )
#        v = from_type_constant(R)
        T = default_int_type(R)
        v = T( from_type_constant(R) )
        :( BInt{R,$T}( $v, 0 ) )
end

# construct BInt from any Integer type
# includes run-time check on bounds

@generated function Base.call{R,Tx<:Integer}( ::Type{ BInt{R} }, x::Tx )

        T = default_int_type( R )

        Tworking = default_int_type( range_union( R, Tx ) )

        l,u = get_range_bounds( R, Tworking )

        :(
                if ( $l <= $Tworking(x) );
                        if ( $Tworking(x) <= $u );
                                return BInt{R,$T}( x % $T, 0 );
                        end;
                end;
                throw( InexactError() );
        )
end

# making a BInt from a BInt only supported
# from narrow to an at-least-as-wide range
@generated function Base.call{R,Rx,Tx<:Integer}( ::Type{ BInt{R} }, x::BInt{Rx,Tx} )
        if !range_within_range( Rx, R )
                throw( InexactError() );
        end
        # know that R is wider than Rx
        # so no further need to check anything
        T = default_int_type( R )
        :(
                UncheckedBInt{R,$T}( x );
        )
end

# @code_llvm BInt{1:10 }( 5 )
# @code_llvm BInt{1:10 }( Int128(5) )
# @code_llvm BInt{1:10 }( Int16(5) )

@generated function Base.call{R,T,Tx<:Integer}( ::Type{ BInt{R,T} }, x::Tx )

        if !range_within_range( get_range_bounds( R )..., typemin(T), typemax(T) )
                throw( InexactError() );
        end

        Tworking = default_int_type( range_union( R, T ) )
        l,u = get_range_bounds( R, Tworking )

        :(
                if ( $Tworking(x) >= $l ) ;
                        if ( $Tworking(x) <= $u ) ;
                                return UncheckedBInt{R,T}( x )
                        end;
                end;
                throw( InexactError() );
        )
end
@generated function Base.call{R,T,Rx,Tx<:Integer}( ::Type{ BInt{R,T} }, x::BInt{Rx,Tx} )
        if !range_within_range( Rx, R )
                throw( InexactError() );
        end
        :(      #UncheckedBInt does the check that T is suitable for R
                UncheckedBInt{R,T}( x );
        )
end

# @code_llvm BInt{1:10,Int32 }( Int32(5) )
# @code_llvm BInt{1:10,Int32 }( Int16(5) )
# @code_llvm BInt{1:10,Int32 }( 5 )

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

# The % operator:
# The converts any Integer (including BInts ) to a BInt by truncating MSBS
# Should be very quick because no run-time range checks are required

@inline function sign_extend_from( n::Integer, b )
        sgn = ( n >> b ) & 1 == 1
        msbs = sgn ? ~zero(n) << b : ~( ~zero(n)<<b )

        return sgn ? n | msbs : n & msbs
end

import Base.%

# conversion from BInt to any Integer:
function %{Tres <: Integer,R,T}( n::BInt{R,T}, ::Type{Tres} )
        return n.n % Tres
end

@generated function %{Rn,Tn,R}( n::BInt{Rn,Tn}, ::Type{BInt{R}} )
        Td = default_int_type(R)
        l,u = get_range_bounds(R)
#        ln,un = get_range_bounds(Tn)

        if range_within_range( Rn, R )
                # the target range R includes all of the source range Rn
                # just convert - no modulo required
                return :( UncheckedBInt{R,$Td}( $Td(n.n) ) )
        end
        if is_unsigned_binary_range( l, u )
                nbits = number_of_bits_unsigned(u)
                mask = Td( ( big(1)<<nbits ) - 1 )
                return :( UncheckedBInt{R,$Td}( $Td(n.n) & $mask ) )
        elseif is_signed_2sComplement_range( l, u )
                b = number_of_bits_signed(u) -1
                msbs1 = ~zero(Tn) << b;
                msbs0 = ~( ~zero(Tn)<<b );
                :(
                        # nn = sign_extend_from( n, $nbits );
                        # ideally would exploit llvm sext instruction...
                        sgn = ( n.n >> $b ) & 1 == 1;
                        nn = sgn ? n.n | $msbs1 : n.n & $msbs0;
                        UncheckedBInt{R,$Td}( nn )
                        )
        else
                error( "Type range not supported for % operator" )
        end
end

@generated function %{Tn<:Integer, R}( n::Tn, ::Type{ BInt{R} } )
        Td = default_int_type(R)
        l,u = get_range_bounds(R)
#        ln,un = get_range_bounds(Tn)

        if range_within_range( Tn, R )
                # just convert - no modulo required
                return :( UncheckedBInt{$R,$Td}( $Td(n) ) )
        end
        if is_unsigned_binary_range( l, u )
                nbits = number_of_bits_unsigned(u)
                mask = Td( ( big(1)<<nbits ) - 1 )
                return :( UncheckedBInt{$R,$Td}( $Td(n) & $mask ) )
        elseif is_signed_2sComplement_range( l, u )
                b = number_of_bits_signed(u) -1
                msbs1 = ~zero(n) << b;
                msbs0 = ~( ~zero(n)<<b );
                return :(
                        # nn = sign_extend_from( n, $nbits );
                        # ideally would exploit llvm sext instruction...
                        sgn = ( n >> $b ) & 1 == 1;
                        nn = sgn ? n | $msbs1 : n & $msbs0;
                        UncheckedBInt{$R,$Td}( nn )
                        )
        else
                error( "Type range not supported for % operator" )
        end
end
@generated function %{Tn <: Integer,R,T}( n::Tn, ::Type{ BInt{R,T} } )
        l,u = get_range_bounds(R)

        if is_unsigned_binary_range( l, u )
                nbits = number_of_bits_unsigned(u)
                mask = T( ( big(1)<<nbits ) - 1 )
                return :( UncheckedBInt{$R,$T}( $T(n) & $mask ) )
        elseif is_signed_2sComplement_range( l, u )
                nbits = number_of_bits_signed(u) -1
                return :(
                        # nn = sign_extend_from( n, $nbits );
                        b = $nbits;
                        sgn = ( n >> b ) & 1 == 1;
                        msbs = sgn ? ~zero(n) << b : ~( ~zero(n)<<b );

                        nn = sgn ? n | msbs : n & msbs;
                        UncheckedBInt{$R,$T}( nn )
                        )
        else
                error( "Type range not supported for % operator\nRHS BInt needs power of two range" )
        end
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiExtensionsV1beta1IngressRule <: SwaggerModel
    host::Any # spec type: Union{ Nothing, String } # spec name: host
    http::Any # spec type: Union{ Nothing, IoK8sApiExtensionsV1beta1HTTPIngressRuleValue } # spec name: http

    function IoK8sApiExtensionsV1beta1IngressRule(;host=nothing, http=nothing)
        o = new()
        validate_property(IoK8sApiExtensionsV1beta1IngressRule, Symbol("host"), host)
        setfield!(o, Symbol("host"), host)
        validate_property(IoK8sApiExtensionsV1beta1IngressRule, Symbol("http"), http)
        setfield!(o, Symbol("http"), http)
        o
    end
end # type IoK8sApiExtensionsV1beta1IngressRule

const _property_map_IoK8sApiExtensionsV1beta1IngressRule = Dict{Symbol,Symbol}(Symbol("host")=>Symbol("host"), Symbol("http")=>Symbol("http"))
const _property_types_IoK8sApiExtensionsV1beta1IngressRule = Dict{Symbol,String}(Symbol("host")=>"String", Symbol("http")=>"IoK8sApiExtensionsV1beta1HTTPIngressRuleValue")
Base.propertynames(::Type{ IoK8sApiExtensionsV1beta1IngressRule }) = collect(keys(_property_map_IoK8sApiExtensionsV1beta1IngressRule))
Swagger.property_type(::Type{ IoK8sApiExtensionsV1beta1IngressRule }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiExtensionsV1beta1IngressRule[name]))}
Swagger.field_name(::Type{ IoK8sApiExtensionsV1beta1IngressRule }, property_name::Symbol) =  _property_map_IoK8sApiExtensionsV1beta1IngressRule[property_name]

function check_required(o::IoK8sApiExtensionsV1beta1IngressRule)
    true
end

function validate_property(::Type{ IoK8sApiExtensionsV1beta1IngressRule }, name::Symbol, val)
end


text:
const VIS = Union{VI, Vector{VI}}

function newblock(m::SOItoMOIBridge, n)
    push!(m.blockdims, n)
    m.nblocks += 1
end

isfree(m, v::VI) = v.value in m.free
function unfree(m, v)
    @assert isfree(m, v)
    delete!(m.free, v.value)
end

function _constraintvariable!(m::SOItoMOIBridge{T}, vs::VIS, s::ZS) where T
    blk = newblock(m, -_length(vs))
    for (i, v) in _enumerate(vs)
        setvarmap!(m, v, (blk, i, i, one(T), _getconstant(m, s)))
        unfree(m, v)
    end
    blk
end
vscaling(::Type{<:NS}) = 1.
vscaling(::Type{<:PS}) = -1.
_length(vi::VI) = 1
_length(vi::Vector{VI}) = length(vi)
_enumerate(vi::VI) = enumerate((vi,))
_enumerate(vi::Vector{VI}) = enumerate(vi)
function _constraintvariable!(m::SOItoMOIBridge, vs::VIS, s::S) where S<:Union{NS, PS}
    blk = newblock(m, -_length(vs))
    cst = _getconstant(m, s)
    m.blkconstant[blk] = cst
    for (i, v) in _enumerate(vs)
        setvarmap!(m, v, (blk, i, i, vscaling(S), cst))
        unfree(m, v)
    end
    blk
end
function getmatdim(k::Integer)
    # n*(n+1)/2 = k
    # n^2+n-2k = 0
    # (-1 + sqrt(1 + 8k))/2
    n = div(isqrt(1 + 8k) - 1, 2)
    if n * (n+1) != 2*k
        error("sd dim not consistent")
    end
    n
end
function _constraintvariable!(m::SOItoMOIBridge{T}, vs::VIS, ::DS) where T
    d = getmatdim(length(vs))
    k = 0
    blk = newblock(m, d)
    for i in 1:d
        for j in 1:i
            k += 1
            setvarmap!(m, vs[k], (blk, i, j, i == j ? one(T) : one(T)/2, zero(T)))
            unfree(m, vs[k])
        end
    end
    blk
end
_var(f::SVF) = f.variable
_var(f::VVF) = f.variables
function _throw_error_if_unfree(m, vi::MOI.VariableIndex)
    if !isfree(m, vi)
        error("A variable cannot be constrained by multiple ",
              "`MOI.SingleVariable` or `MOI.VectorOfVariables` constraints.")
    end
end
function _throw_error_if_unfree(m, vis::MOI.Vector)
    for vi in vis
        _throw_error_if_unfree(m, vi)
    end
end
function MOIU.allocate_constraint(m::SOItoMOIBridge{T}, f::VF, s::SupportedSets) where T
    vis = _var(f)
    _throw_error_if_unfree(m, vis)
    blk = _constraintvariable!(m, vis, s)
    if isa(s, ZS)
        ci = _allocate_constraint(m, f, s)
        m.zeroblock[ci] = blk
        return ci
    else
        return CI{typeof(f), typeof(s)}(-blk)
    end
end

_getconstant(m::SOItoMOIBridge, s::MOI.AbstractScalarSet) = MOIU.getconstant(s)
_getconstant(m::SOItoMOIBridge{T}, s::MOI.AbstractSet) where T = zero(T)

_var(f::SVF, j) = f.variable
_var(f::VVF, j) = f.variables[j]
function MOIU.load_constraint(m::SOItoMOIBridge, ci::CI, f::VF, s::SupportedSets)
    if ci.value >= 0 # i.e. s is ZS or _var(f) wasn't free at allocate_constraint
        setconstant!(m, ci, s)
        cs = m.constrmap[ci]
        @assert !isempty(cs)
        for k in 1:length(cs)
            vm = varmap(m, _var(f, k))
            # For free variables, the length of vm is 2, clearly not the case here
            @assert length(vm) == 1
            (blk, i, j, coef, shift) = first(vm)
            c = cs[k]
            setconstraintcoefficient!(m.sdoptimizer, coef, c, blk, i, j)
            setconstraintconstant!(m.sdoptimizer,  _getconstant(m, s) - coef * shift, c)
        end
    end
end

function loadfreevariables!(m::SOItoMOIBridge{T}) where T
    for vi in m.free
        blk = newblock(m, -2)
        # x free transformed into x = y - z with y, z >= 0
        setvarmap!(m, VI(vi), [(blk, 1, 1, one(T), zero(T)), (blk, 2, 2, -one(T), zero(T))])
    end
end


text:
# Simple snake game using Makie.jl, made by Nand Vinchhi for GCI'19
using Makie

score = 0
println("  ______   __    _       __       _   _   _______")
println("/ ______| |  \\  | |     /  \\     | | / / |  _____|")
println("| |____   | \\ \\ | |    / /\\ \\    | |/ /  | |____")
println("|_____  \\ | |\\ \\| |   / /__\\ \\   | | \\   |  ____|")
println(" _____| | | | \\ \\ |  / /----\\ \\  | |\\ \\  | |_____")
println("|_______/ |_|  \\__| /_/      \\_\\ |_| \\_\\ |_______|")
println("")
println("use the arrow keys to move and collect the red fruit. The snake gets faster and faster as the game progresses!")
scene = Scene(resolution = (400, 400), raw = true, camera = campixel!)
display(scene)

k = [[0, 20], [0, 20], [0, 20], [0, 20]]
dir = lift(scene.events.keyboardbuttons) do but

    if ispressed(but, Keyboard.left)
        global k[1] = [-20, 0]

    elseif ispressed(but, Keyboard.up)
        global k[1] = [0, 20]

    elseif ispressed(but, Keyboard.right)
        global k[1] = [20, 0]

    elseif ispressed(but, Keyboard.down)
        global k[1] = [0, -20]

    end

end

y = [70, 50, 30, 10]
x = [10, 10, 10, 10]
fruit_x = 90
fruit_y = 110
game_over = false
speed = 0.3
scene = scatter!([fruit_x], [fruit_y], color = :red, marker = :rect, markersize = 20)

function iter()

    global scene = scatter!(x, y, color = :white, marker = :rect, markersize = 25)

    for i1 = 0:(length(x) - 2)
        global n1 = length(x) - i1
        if k[n1 - 1] == [(x[n1 - 1] - x[n1]), (y[n1 - 1] - y[n1])]
            global k[n1] = k[n1 - 1]
        end
    end

    for f = 2:length(x)
        if x[f] == x[1] && y[f] == y[1]
            global game_over = true
        end
    end

    if x[1] < 10 || x[1] > 390 || y[1] < 10 || y[1] > 390
        global game_over = true
    end

    for i = 1:length(x)
        global x[i] += k[i][1]
        global y[i] += k[i][2]
    end

    if x[1] == fruit_x && y[1] == fruit_y
        global scene = scatter!([fruit_x], [fruit_y], color = :white, marker = :rect, markersize = 25)
        global fruit_x = rand([10, 30, 50, 70, 90, 110, 130, 150, 170, 190, 210, 230, 250, 270, 290, 310, 330, 350, 370, 390])
        global fruit_y = rand([10, 30, 50, 70, 90, 110, 130, 150, 170, 190, 210, 230, 250, 270, 290, 310, 330, 350, 370, 390])
        global scene = scatter!([fruit_x], [fruit_y], color = :red, marker = :rect, markersize = 20)
        push!(x, x[end] - k[end][1])
        push!(y, y[end] - k[end][2])
        push!(k, k[end])
        global speed -= 0.01
        global score += 1
    end
    global scene = scatter!(x, y, color = :green, marker = :rect, markersize = 20)

    sleep(speed)

end

while game_over == false
    iter()
end

for i = 1:100
    println("")
end

println("Game Over. Better Luck next time!")
println("Your score was " * string(score) * ".")
for i3 = 1:3
    scene = scatter!(x, y, color = :white, marker = :rect, markersize = 25)
    sleep(0.5)
    scene = scatter!(x, y, color = :green, marker = :rect, markersize = 20)
    sleep(0.5)
end


text:
using Example2018

@show Polynomial([1,2,3])
@show evaluate(Polynomial([1,2,3]), 2)
@show Polynomial([1,2,3]) + Polynomial([4,5])
@show 3 * Polynomial([1,2,4,42])
@show deriv(Polynomial([1,1,1,1,1]))
@show deriv(Polynomial([0,0,1,0,0]))


text:
##
using CSV
datadir = "/Users/dgleich/Downloads/spamdir/"
using SparseArrays
function writeSMAT(filename::AbstractString, A::SparseMatrixCSC{T,Int}; values::Bool=true) where T
    open(filename, "w") do outfile
        write(outfile, join((size(A,1), size(A,2), nnz(A)), " "), "\n")

        rows = rowvals(A)
        vals = nonzeros(A)
        m, n = size(A)
        for j = 1:n
           for nzi in nzrange(A, j)
              row = rows[nzi]
              val = vals[nzi]
              if values
                write(outfile, join((row-1, j-1, val), " "), "\n")
              else
                write(outfile, join((row-1, j-1, 1), " "), "\n")
              end
           end
        end
    end
end

## 2000 - harvested area
data2000h = CSV.read("$datadir/spam_h.csv")
data = data2000h
data = data[data.stat_code .== "USA",:]
# filter out hawaii & alaska...
data=data[data.x .>= -140,:]
crops = names(data)[10:29]
X = Matrix(data[:,10:29])
##
#=
using LinearAlgebra
f = norm.(eachrow(X)) .> 0.1)
data = data[f,:]
X = X[f,:]
=#
##
include("nn.jl")
Y = copy(X')
@time idxs,dsts = nnsimple(Y,25)
##
# G = nngraph(idxs[1:45,:],dsts[1:45,:])
G = nngraph(idxs,dsts)

##
using Pkg
if haskey(Pkg.installed(), "MatrixNetworks")
  using MatrixNetworks
  @show is_connected(G)
end

##
using DelimitedFiles
writeSMAT("cropsim-2000-harvestarea.smat",G)
writedlm("cropsim-2000-harvestarea.xy", zip(data.x,data.y))
##
writedlm("cropsim-2000-harvestarea-metadata.csv",
    zip(data.alloc_key, data.hc_seq5m, sum.(eachrow(X)), norm.(eachrow(X)), maximum.(eachrow(X))), ',')

## 2000 - yield

data2000y = CSV.read("$datadir/spam_y.csv")
data = data2000y
data = data[data.stat_code .== "USA",:]
# filter out hawaii & alaska...
data=data[data.x .>= -140,:]
crops = names(data)[10:29]
X = Matrix(data[:,10:29])

##
include("nn.jl")
Y = copy(X')
@time idxs,dsts = nnsimple(Y,25)
G = nngraph(idxs,dsts)

##
using DelimitedFiles
writeSMAT("cropsim-2000-yield.smat",G)
writedlm("cropsim-2000-yield.xy", zip(data.x,data.y))
##
# hc_seq5m is the same as cell5m in later datasets
writedlm("cropsim-2000-yield-metadata.csv",
    zip(data.alloc_key, data.hc_seq5m, sum.(eachrow(X)), norm.(eachrow(X)), maximum.(eachrow(X))), ',')

## 2005 - harvested area
# crops = names(data)[10:51]
data2005h = CSV.read("$datadir/spam2005V3r2_global_H_TA.csv")
data = data2005h
## This data doesn't have x,y so we add them.
"""
CELL5M is a 5-minute arc resolution view of the -180,-90
lat-long system. This corresponds to a cell-size of
0.08333333333. The following code converts between cell5m
values and latlong.
"""
function cell5m_to_latlong(cell::Int)
  nrows = 4320
  ncols = 2160
  col, row = divrem(cell, nrows)
  x = row*0.0833333333333333333-(180-0.0833333333333333333/2)
  y = (90-0.0833333333333333333/2) - col*0.0833333333333333333
  return x,y
end
function add_xy_coords!(df)
  cell5m = df.cell5m
  x = zeros(length(cell5m))
  y = copy(x)
  for i=1:length(cell5m)
    x[i],y[i] = cell5m_to_latlong(cell5m[i])
  end
  df.x = x
  df.y = y
end
add_xy_coords!(data)
data = data[data.iso3 .== "USA",:]
# filter out hawaii & alaska...
data=data[data.x .>= -140,:]
crops = names(data)[8:49]
X = Matrix(data[:,8:49])

##
include("nn.jl")
Y = copy(X')
@time idxs,dsts = nnsimple(Y,25)
G = nngraph(idxs,dsts)
##
using DelimitedFiles
writeSMAT("cropsim-2005-harvestarea.smat",G)
writedlm("cropsim-2005-harvestarea.xy", zip(data.x,data.y))
##
using LinearAlgebra
writedlm("cropsim-2005-harvestarea-metadata.csv",
    zip(data.alloc_key, data.cell5m, sum.(eachrow(X)), norm.(eachrow(X)), maximum.(eachrow(X))), ',')

## 2005 - yield
data2005y = CSV.read("$datadir/spam2005V3r2_global_Y_TA.csv")
data = data2005y
add_xy_coords!(data)
data = data[data.iso3 .== "USA",:]
# filter out hawaii & alaska...
data=data[data.x .>= -140,:]
crops = names(data)[8:49]
X = Matrix(data[:,8:49])

##
include("nn.jl")
Y = copy(X')
@time idxs,dsts = nnsimple(Y,25)
G = nngraph(idxs,dsts)
##
using DelimitedFiles
using LinearAlgebra

writeSMAT("cropsim-2005-yield.smat",G)
writedlm("cropsim-2005-yield.xy", zip(data.x,data.y))
writedlm("cropsim-2005-yield-metadata.csv",
    zip(data.alloc_key, data.cell5m, sum.(eachrow(X)), norm.(eachrow(X)), maximum.(eachrow(X))), ',')

## 2010 - harvested area
# crops = names(data)[10:51]
data2010h = CSV.read("$datadir/spam2010V1r1_global_H_TA.csv")
##
data = data2010h
data = data[data.iso3 .== "USA",:]
# filter out hawaii & alaska...
data=data[data.x .>= -140,:]
crops = names(data)[10:51]
X = Matrix(data[:,10:51])

##
include("nn.jl")
Y = copy(X')
@time idxs,dsts = nnsimple(Y,25)
G = nngraph(idxs,dsts)
##
using DelimitedFiles
using LinearAlgebra

writeSMAT("cropsim-2010-harvestarea.smat",G)
writedlm("cropsim-2010-harvestarea.xy", zip(data.x,data.y))
writedlm("cropsim-2010-harvestarea-metadata.csv",
    zip(data.alloc_key, data.cell5m, sum.(eachrow(X)), norm.(eachrow(X)), maximum.(eachrow(X))), ',')

## 2010 - yield
data2010y = CSV.read("$datadir/spam2010V1r1_global_Y_TA.csv")
##
data = data2010y
data = data[data.iso3 .== "USA",:]
# filter out hawaii & alaska...
data=data[data.x .>= -140,:]
crops = names(data)[10:51]
X = Matrix(data[:,10:51])

##
include("nn.jl")
Y = copy(X')
@time idxs,dsts = nnsimple(Y,25)
G = nngraph(idxs,dsts)
##
using DelimitedFiles
using LinearAlgebra

writeSMAT("cropsim-2010-yield.smat",G)
writedlm("cropsim-2010-yield.xy", zip(data.x,data.y))
writedlm("cropsim-2010-yield-metadata.csv",
    zip(data.alloc_key, data.cell5m, sum.(eachrow(X)), norm.(eachrow(X)), maximum.(eachrow(X))), ',')

max_stars_repo_name: empet/more-graphs

max_issues_repo_name: empet/more-graphs

max_forks_repo_name: empet/more-graphs


text:
module Sources

using Adapt
using DocStringExtensions

import ..Parameters as RP

export AbstractSourceLW, SourceLWNoScat, SourceLW2Str, SourceSW2Str, source_func_longwave, source_func_shortwave

"""
    AbstractSourceLW{FT,FTA1D,FTA2D}

Abstract longwave source for no-scattering and two stream longwave solvers.
"""
abstract type AbstractSourceLW{FT, FTA1D, FTA2D} end

"""
    SourceLWNoScat{FT,FTA1D,FTA2D} <: AbstractSourceLW{FT,FTA1D,FTA2D}

Longwave sources: computed at layer center, layer edges,
and at the surface for no scattering calculations

# Fields

$(DocStringExtensions.FIELDS)
"""
struct SourceLWNoScat{
    FT <: AbstractFloat,
    FTA1D <: AbstractArray{FT, 1},
    FTA2D <: AbstractArray{FT, 2},
    PS <: RP.ARP,
} <: AbstractSourceLW{FT, FTA1D, FTA2D}
    "Parameter set"
    param_set::PS
    "Planck source at layer average temperature `[W/m2]` `(nlay, ncol)`"
    lay_source::FTA2D
    "Planck source at layer edge in increasing ilay direction `[W/m2]` `(nlay, ncol)`, includes spectral weighting that accounts for state-dependent frequency to g-space mapping"
    lev_source_inc::FTA2D
    "Planck source at layer edge in decreasing ilay direction `[W/m2]` `(nlay, ncol)`, includes spectral weighting that accounts for state-dependent frequency to g-space mapping"
    lev_source_dec::FTA2D
    "Surface source `[W/m2]` `(ncol)`"
    sfc_source::FTA1D
    "temporary storage array used in no scattering GPU calculations `(nlay, ncol)`"
    src_up::FTA2D
    "temporary storage array used in no scattering GPU calculations `(nlay, ncol)`"
    src_dn::FTA2D
end
Adapt.@adapt_structure SourceLWNoScat

function SourceLWNoScat(param_set::RP.ARP, ::Type{FT}, ::Type{DA}, nlay::Int, ncol::Int) where {FT <: AbstractFloat, DA}
    FTA1D = DA{FT, 1}
    FTA2D = DA{FT, 2}

    lay_source = FTA2D(undef, nlay, ncol) # lay_source
    lev_source_inc = FTA2D(undef, nlay, ncol) # lev_source_inc
    lev_source_dec = FTA2D(undef, nlay, ncol) # lev_source_dec
    sfc_source = FTA1D(undef, ncol)       # sfc_source
    src_up = FTA2D(undef, nlay, ncol) # src_up
    src_dn = FTA2D(undef, nlay, ncol) # src_dn

    return SourceLWNoScat{eltype(lay_source), typeof(sfc_source), typeof(lay_source), typeof(param_set)}(
        param_set,
        lay_source,
        lev_source_inc,
        lev_source_dec,
        sfc_source,
        src_up,
        src_dn,
    )
end

"""
    SourceLW2Str{FT,FTA1D,FTA2D} <: AbstractSourceLW{FT,FTA1D,FTA2D}

Longwave sources: computed at layer center, layer edges,
and at the surface for 2-stream calculations

# Fields

$(DocStringExtensions.FIELDS)
"""
struct SourceLW2Str{FT <: AbstractFloat, FTA1D <: AbstractArray{FT, 1}, FTA2D <: AbstractArray{FT, 2}, PS <: RP.ARP} <:
       AbstractSourceLW{FT, FTA1D, FTA2D}
    "Parameter set"
    param_set::PS
    "Planck source at layer average temperature `[W/m2]` `(nlay, ncol)`"
    lay_source::FTA2D
    "Planck source at layer edge in increasing ilay direction `[W/m2]` `(nlay, ncol)`, includes spectral weighting that accounts for state-dependent frequency to g-space mapping"
    lev_source_inc::FTA2D
    "Planck source at layer edge in decreasing ilay direction `[W/m2]` `(nlay, ncol)`, includes spectral weighting that accounts for state-dependent frequency to g-space mapping"
    lev_source_dec::FTA2D
    "Surface source `[W/m2]` `(ncol)`"
    sfc_source::FTA1D
    "level source `[W/m2]` `(nlay+1, ncol)`, used in 2 stream calculations"
    lev_source::FTA2D
    "temporary storage array, used in 2 stream calculations `(nlay, ncol)`"
    src_up::FTA2D
    "temporary storage array, used in 2 stream calculations `(nlay, ncol)`"
    src_dn::FTA2D
    "temporary storage array, used in 2 stream calculations `(nlay, ncol)`"
    Rdif::FTA2D
    "temporary storage array, used in 2 stream calculations `(nlay, ncol)`"
    Tdif::FTA2D
    "temporary storage array, used in 2 stream calculations `(nlay + 1, ncol)`"
    albedo::FTA2D
    "temporary storage array, used in 2 stream calculations `(nlay + 1, ncol)`"
    src::FTA2D
end
Adapt.@adapt_structure SourceLW2Str

function SourceLW2Str(param_set::RP.ARP, ::Type{FT}, ::Type{DA}, nlay::Int, ncol::Int) where {FT <: AbstractFloat, DA}
    FTA1D = DA{FT, 1}
    FTA2D = DA{FT, 2}

    lay_source = FTA2D(undef, nlay, ncol) # lay_source
    lev_source_inc = FTA2D(undef, nlay, ncol) # lev_source_inc
    lev_source_dec = FTA2D(undef, nlay, ncol) # lev_source_dec
    sfc_source = FTA1D(undef, ncol) # sfc_source
    lev_source = FTA2D(undef, nlay + 1, ncol) # lev_source
    src_up = FTA2D(undef, nlay, ncol) # src_up
    src_dn = FTA2D(undef, nlay, ncol) # src_dn
    Rdif = FTA2D(undef, nlay, ncol) # Rdif
    Tdif = FTA2D(undef, nlay, ncol) # Tdif
    albedo = FTA2D(undef, nlay + 1, ncol) # albedo
    src = FTA2D(undef, nlay + 1, ncol) # src

    return SourceLW2Str{eltype(lay_source), typeof(sfc_source), typeof(lay_source), typeof(param_set)}(
        param_set,
        lay_source,
        lev_source_inc,
        lev_source_dec,
        sfc_source,
        lev_source,
        src_up,
        src_dn,
        Rdif,
        Tdif,
        albedo,
        src,
    )
end

"""
    source_func_longwave(
        param_set,
        ::Type{FT},
        ncol::Int,
        nlay::Int,
        OPC::Symbol,
        ::Type{DA},
    ) where {FT<:AbstractFloat,DA}

Initializes the longwave source for one scalar and two stream simulations.
"""
function source_func_longwave(
    param_set::RP.ARP,
    ::Type{FT},
    ncol::Int,
    nlay::Int,
    OPC::Symbol,
    ::Type{DA},
) where {FT <: AbstractFloat, I <: Int, DA}
    lay_source = DA{FT, 2}(undef, nlay, ncol)
    lev_source_inc = DA{FT, 2}(undef, nlay, ncol)
    lev_source_dec = DA{FT, 2}(undef, nlay, ncol)
    sfc_source = DA{FT, 1}(undef, ncol)
    src_up = DA{FT, 2}(undef, nlay, ncol)
    src_dn = DA{FT, 2}(undef, nlay, ncol)

    if OPC === :OneScalar
        return SourceLWNoScat{eltype(lay_source), typeof(sfc_source), typeof(src_up), typeof(param_set)}(
            param_set,
            lay_source,
            lev_source_inc,
            lev_source_dec,
            sfc_source,
            src_up,
            src_dn,
        )
    else
        PS = typeof(param_set)
        lev_source = DA{FT, 2}(undef, nlay + 1, ncol)
        Rdif = DA{FT, 2}(undef, nlay, ncol)
        Tdif = DA{FT, 2}(undef, nlay, ncol)
        albedo = DA{FT, 2}(undef, nlay + 1, ncol)
        src = DA{FT, 2}(undef, nlay + 1, ncol)
        return SourceLW2Str{FT, typeof(sfc_source), typeof(src_up), PS}(
            param_set,
            lay_source,
            lev_source_inc,
            lev_source_dec,
            sfc_source,
            lev_source,
            src_up,
            src_dn,
            Rdif,
            Tdif,
            albedo,
            src,
        )
    end
end

"""
    SourceSW2Str{FT,FTA1D,FTA2D}

Shortwave sources: computed at layer center, layer edges,
and at the surface for 2-stream calculations

# Fields

$(DocStringExtensions.FIELDS)
"""
struct SourceSW2Str{FT <: AbstractFloat, FTA1D <: AbstractArray{FT, 1}, FTA2D <: AbstractArray{FT, 2}}
    "surface source `(ncol)`"
    sfc_source::FTA1D
    "albedo `(nlay + 1, ncol)`"
    albedo::FTA2D
    "temporary storage array, used in 2 stream calculations `(nlay, ncol)`"
    src_up::FTA2D
    "temporary storage array, used in 2 stream calculations `(nlay, ncol)`"
    src_dn::FTA2D
    "temporary storage array, used in 2 stream calculations `(nlay, ncol)`"
    Rdif::FTA2D
    "temporary storage array, used in 2 stream calculations `(nlay, ncol)`"
    Tdif::FTA2D
    "temporary storage array, used in 2 stream calculations `(nlay, ncol)`"
    Rdir::FTA2D
    "temporary storage array, used in 2 stream calculations `(nlay, ncol)`"
    Tdir::FTA2D
    "temporary storage array, used in 2 stream calculations `(nlay, ncol)`"
    Tnoscat::FTA2D
    "temporary storage array, used in 2 stream calculations `(nlay + 1, ncol)`"
    src::FTA2D
end
Adapt.@adapt_structure SourceSW2Str

function SourceSW2Str(::Type{FT}, ::Type{DA}, nlay::I, ncol::I) where {FT <: AbstractFloat, I <: Int, DA}
    FTA1D = DA{FT, 1}
    FTA2D = DA{FT, 2}

    sfc_source = FTA1D(undef, ncol)
    albedo = FTA2D(undef, nlay + 1, ncol)
    src_up = FTA2D(undef, nlay, ncol)
    src_dn = FTA2D(undef, nlay, ncol)
    Rdif = FTA2D(undef, nlay, ncol)
    Tdif = FTA2D(undef, nlay, ncol)
    Rdir = FTA2D(undef, nlay, ncol)
    Tdir = FTA2D(undef, nlay, ncol)
    Tnoscat = FTA2D(undef, nlay, ncol)
    src = FTA2D(undef, nlay + 1, ncol)

    return SourceSW2Str{eltype(sfc_source), typeof(sfc_source), typeof(albedo)}(
        sfc_source,
        albedo,
        src_up,
        src_dn,
        Rdif,
        Tdif,
        Rdir,
        Tdir,
        Tnoscat,
        src,
    )
end

"""
    source_func_shortwave(
        ::Type{FT},
        ncol::I,
        nlay::I,
        opc::Symbol,
        ::Type{DA},
    ) where {FT<:AbstractFloat,I<:Int,DA}

Initializes the shortwave source for one scalar and two stream simulations.
"""
function source_func_shortwave(
    ::Type{FT},
    ncol::I,
    nlay::I,
    opc::Symbol,
    ::Type{DA},
) where {FT <: AbstractFloat, I <: Int, DA}
    if opc == :OneScalar
        return nothing
    else
        sfc_source = DA{FT}(undef, ncol)
        albedo = DA{FT}(undef, nlay + 1, ncol)
        src_up = DA{FT}(undef, nlay, ncol)
        src_dn = DA{FT}(undef, nlay, ncol)
        Rdif = DA{FT}(undef, nlay, ncol)
        Tdif = DA{FT}(undef, nlay, ncol)
        Rdir = DA{FT}(undef, nlay, ncol)
        Tdir = DA{FT}(undef, nlay, ncol)
        Tnoscat = DA{FT}(undef, nlay, ncol)
        src = DA{FT}(undef, nlay + 1, ncol)
        return SourceSW2Str{FT, typeof(sfc_source), typeof(src_up)}(
            sfc_source,
            albedo,
            src_up,
            src_dn,
            Rdif,
            Tdif,
            Rdir,
            Tdir,
            Tnoscat,
            src,
        )
    end
end

end

max_stars_repo_name: climate-machine/JRRTMGP

max_issues_repo_name: climate-machine/JRRTMGP

max_forks_repo_name: climate-machine/JRRTMGP


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

# See https://github.com/JuliaLang/Pkg.jl/issues/2942
# Once this Pkg issue is resolved, this must be removed
uuid = Base.UUID("a83860b7-747b-57cf-bf1f-3e79990d037f")
delete!(Pkg.Types.get_last_stdlibs(v"1.6.3"), uuid)

julia_versions = [v"1.6.3", v"1.7.0", v"1.8.0", v"1.9.0"]

name = "jlqml"
version = v"0.3"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/barche/jlqml.git", "4d21e2fcf0bcae29c45e76c8ac9d1d45893a77e6"),
]

# Bash recipe for building across all platforms
script = raw"""
if test -f "$prefix/lib/cmake/Qt5Gui/Qt5GuiConfigExtras.cmake"; then
    sed -i 's/_qt5gui_find_extra_libs.*AGL.framework.*//' $prefix/lib/cmake/Qt5Gui/Qt5GuiConfigExtras.cmake
fi

# Override compiler ID to silence the horrible "No features found" cmake error
if [[ $target == *"apple-darwin"* ]]; then
  macos_extra_flags="-DCMAKE_CXX_COMPILER_ID=AppleClang -DCMAKE_CXX_COMPILER_VERSION=10.0.0 -DCMAKE_CXX_STANDARD_COMPUTED_DEFAULT=11"
fi

mkdir build
cd build
cmake \
    -DCMAKE_BUILD_TYPE=Release \
    -DCMAKE_FIND_ROOT_PATH=${prefix} \
    -DCMAKE_INSTALL_PREFIX=$prefix \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DJulia_PREFIX=${prefix} \
    -DQt5Core_DIR=$prefix/lib/cmake/Qt5Core \
    -DQt5Quick_DIR=$prefix/lib/cmake/Qt5Quick \
    -DQt5Svg_DIR=$prefix/lib/cmake/Qt5Svg \
    -DQt5Widgets_DIR=$prefix/lib/cmake/Qt5Widgets \
    $macos_extra_flags \
    ../jlqml/
VERBOSE=ON cmake --build . --config Release --target install -- -j${nproc}
install_license $WORKSPACE/srcdir/jlqml*/LICENSE.md
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
include("../../L/libjulia/common.jl")
platforms = vcat(libjulia_platforms.(julia_versions)...)
platforms = expand_cxxstring_abis(platforms)

filter!(p -> libc(p) != "musl", platforms) # Qt_jll is currently not available for muslc
# Qt5Declarative_jll is not available for these architectures:
filter!(p -> arch(p) != "armv6l", platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libjlqml", :libjlqml),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("libcxxwrap_julia_jll"),
    Dependency("Qt5Declarative_jll"),
    Dependency("Qt5Svg_jll"),
    BuildDependency("Libglvnd_jll"),
    BuildDependency("libjulia_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
    preferred_gcc_version = v"9",
    julia_compat = "1.6")

max_stars_repo_name: MichelJuillard/Yggdrasil

max_issues_repo_name: MichelJuillard/Yggdrasil

max_forks_repo_name: MichelJuillard/Yggdrasil


text:
"""
    Coulomb(nl_only)

The Coulomb electrostatic interaction.
"""
struct Coulomb{C, T} <: GeneralInteraction
    cutoff::C
    nl_only::Bool
    coulomb_const::T
end

Coulomb(cutoff, nl_only) = Coulomb(
    cutoff,
    nl_only,
    138.935458 / 70.0 # Treat ϵr as 70 for now
)

Coulomb(nl_only=false) = Coulomb(
    NoCutoff(),
    nl_only,
    138.935458 / 70.0 # Treat ϵr as 70 for now
)

@inline @inbounds function force(inter::Coulomb{C},
                                    coord_i,
                                    coord_j,
                                    atom_i,
                                    atom_j,
                                    box_size) where C
    dr = vector(coord_i, coord_j, box_size)
    r2 = sum(abs2, dr)

    cutoff = inter.cutoff
    coulomb_const = inter.coulomb_const
    qi, qj = atom_i.charge, atom_j.charge

    params = (coulomb_const, qi, qj)

    if cutoff_points(C) == 0
        f = force_nocutoff(inter, r2, inv(r2), params)
    elseif cutoff_points(C) == 1
        sqdist_cutoff = cutoff.sqdist_cutoff
        r2 > sqdist_cutoff && return zero(coord_i)

        f = force_cutoff(cutoff, r2, inter, params)
    elseif cutoff_points(C) == 2
        sqdist_cutoff = cutoff.sqdist_cutoff
        activation_dist = cutoff.activation_dist

        r2 > sqdist_cutoff && return zero(coord_i)

        if r2 < activation_dist
            f = force_nocutoff(inter, r2, inv(r2), params)
        else
            f = force_cutoff(cutoff, r2, inter, params)
        end
    end

    return f * dr
end

@fastmath function force_nocutoff(::Coulomb, r2, invr2, (coulomb_const, qi, qj))
    (coulomb_const * qi * qj) / √(r2 ^ 3)
end

@inline @inbounds function potential_energy(inter::Coulomb{C},
                                    s::Simulation,
                                    i::Integer,
                                    j::Integer) where C
    U = eltype(s.coords[i]) # this is not Unitful compatible
    i == j && return zero(U)

    dr = vector(s.coords[i], s.coords[j], s.box_size)
    r2 = sum(abs2, dr)

    cutoff = inter.cutoff
    coulomb_const = inter.coulomb_const
    qi = s.atoms[i].charge
    qj = s.atoms[j].charge
    params = (coulomb_const, qi, qj)

    if cutoff_points(C) == 0
        potential(inter, r2, inv(r2), params)
    elseif cutoff_points(C) == 1
        sqdist_cutoff = cutoff.sqdist_cutoff * σ2
        r2 > sqdist_cutoff && return zero(U)

        potential_cutoff(cutoff, r2, inter, params)
    elseif cutoff_points(C) == 2
        r2 > cutoff.sqdist_cutoff && return zero(U)

        if r2 < activation_dist
            potential(inter, r2, inv(r2), params)
        else
            potential_cutoff(cutoff, r2, inter, params)
        end
    end
end

@fastmath function potential(::Coulomb, r2, invr2, (coulomb_const, qi, qj))
    (coulomb_const * qi * qj) * √invr2
end


text:
# SYSTEM INFO:
julia> versioninfo()
Julia Version 1.1.0
Commit 80516ca202 (2019-01-21 21:24 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i7-6600U CPU @ 2.60GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, skylake)

# colored bar chart:
bar(rand(10); fillcolor = 1:10)

# to update package AND Toml
] free <package>
] update <package>

# CREATING CODE FROM A STRING
"""
The built-in function Meta.parse takes a string and
transforms it into an expression. This expression
can be evaluated in Julia with the function Core.eval.
For example:
"""
julia> expr = Meta.parse("1+2*3")
:(1 + 2 * 3)
julia> eval(expr)
7
julia> expr = Meta.parse("sqrt(π)")
:(sqrt(π))
julia> eval(expr)
1.7724538509055159

# MULTICORE PARALLEL example, 1=head, 0=tail
# How many cores do we have>
Sys.CPU_THREADS

using Distributed
heads = @distributed (+) for i=1:200000000
  Int(floor(rand()+0.5))
end
print(heads,"\n")

# Calling a c lib function
#using Distributed already called
#t1=ccall((:clock, "libc"), UInt32,())
@time heads = @distributed (+) for i=1:200000000
  Int(floor(rand()+0.5))
end
print(heads,"\n")
#t2=ccall((:clock, "libc"), UInt32,())
#print((t2-t1)/1000000.," seconds\n")
# Don't need the C clock function since Julia has a
# @time macro and other timing operations.

# UNICODE characters in strings and formulas
s="abπ∑ef\n"
print(s)
ß=2π/3

function Σ(x)  # yes, function names can use UTF-8 chars
    s=0
    for i=1:length(x)
        s+=x[i]
    end
    return s
end

x=1:100
println(Σ(x))

x=[1.0 2; 3 4]
println(Σ(x))
# the base function sum() does nearly the same thing so
# could have written Σ()=sum()

# With LaTeXStrings you can use
L"y = 0.19^{+1.10}_{-1.05}"

# FWIW, you can definitely use latexstrings in IJulia with no TeX installation
# I think it uses MathJax (https://www.mathjax.org/).

# plotting unicode with GR (not working?)
ENV["GKS_ENCODING"] = "utf-8"

# OS UTILITIES & COMMANDS using run() and read()
pwd()
run(`echo hello`)  # note (`) is a backtick`
files=String(read(`ls`))
files=split(files,"\n")
println(files)

# COMPlEX floating-point numbers
x=2.1+3.2im

# SWAP TWO NUMBERS: Don't need swap macro
a,b=b,a

# checking approx. equality function:
isapprox(3.0, 3.01, rtol=0.1)

# COMPARISON OPERATORS
a = 2
b = 3
c = 3
@show(a < c && b < c)  # then and operator
@show(a < c || b < c)  # the or operator
@show(a != b) ; # not equal

# DICE ROLLS
rand(1:5) # Create a DISCRETE RANDOM DRAW from 1-5

# argmin function:
A = [3 5 2 1 6]
min_loc = argmin(A)
A[min_loc]

## FILTERING: A few ways to filter a variable:
a = [1,2,3,5,1,6]

#Higher-order filtering

z = filter(x -> x > 4, a)
@show(z)

# This calls the filter function with the predicate x -> x > 4 (see Anonymous functions in the Julia manual).
#  Broadcasting and indexing

z2 = a[Bool[a[i] > 4 for i = 1:length(a)]]
@show(z2)

# This performs a broadcasting comparision between the
# elements of a and 4, then uses the resulting array of
# booleans to index a.

# It can be written more compactly using a broadcasting
# operator:

z3 = a[a .> 4]
@show(z3)

## METAPROGRAMMING machinery exposed to the user in Julia
# eval(), Expr(), parse(), Symbol(), typeof(), Regex()

# always check what Julia has to offer before rolling
# your own functions, macros and extensions.

# GETTING USER RESPONSE TO TEXT
println("Who are you?")
s=readline()
println("Hello $s and Hello World!")

# Works much better in a function
# as a function:
function whoru()
	println("Who are you?")
	s=readline()
	println("Hello $s and Hello World!")
end

# USE OF DICTIONARY example
# Bad name dictionary can be as big as one wants
# and severity of baddness can be designated by the number
badnames=Dict("Dope" => 1, "Knothead" => 1,"Stupid" => 1)
# The bot wants a name
println("Please give me a name for I have none.")
name=readline()
# Is it a bad name or a good name
n=get(badnames,name,0)
if n==0
    println("$name is nice, I will take it." )
else
    println("I think that's not a very nice name.")
    name="Sky"  #could use a random name selector here
    println("I chose the name $name for myself.")
    println("$name is a strong name and I like it.")
end

# PLOTTING A BAR GRAPH of a dictionary
c = Dict(1 => 7, 2 => 5, 3 => 2, 4 => 14, 5 => 22)
@show(cx = collect(keys(c)));
@show(cy = collect(values(c)));
plot(cx, cy, st=:bar, bar_width=0.3, label="",size=(300,200), color=1:5)

# with string keys
cs = Dict("awful" => 7, "bad" => 5, "so so" => 2, "good" => 14, "great" => 22)
@show(cx = collect(keys(cs)));
@show(cy = collect(values(cs)));
plot(cx, cy, st=:bar, bar_width=0.3, label="",size=(300,200),color = 1:5)

# SORTPERM to order an array
# define dictionary c with numbers as above
c = Dict(1 => 7, 2 => 5, 3 => 2, 4 => 14, 5 => 22)
cx = collect(keys(c)));
cy = collect(values(c));
@show(zz = [cx cy])
@show(z = zz[sortperm(zz[:, 1]), :])  # sort the array zz by column 1
# BAR PLOT of the dictionary
plot(z[:,1], z[:,2], st=:bar, bar_width=0.3, label="", color = 1:5, title= "strings",
			xticks = (z[:,1],["awful", "bad", "so so", "good", "great"]), size = (300,200))

# STRING TO FLOAT CONVERSION
st=readline();       #Arg STDIN assumed if left out
t=parse(Float32,st)  #string to Float32 conversion

# SPLITTING/FORMATTING PRINT statement
depthu = 10; tu = 2
println("depth $depthu meters \nsound time $tu seconds")

# GLM and DATAFRAMES
# GLM Generalized Linear Model
#Pkg.add("Gadfly")
#using Gadfly
#Pkg.add("DataFrames")
#Pkg.add("GLM")

# Gadfly is similar to ggplot - based on the "grammar of graphics"
using Gadfly
using DataFrames, GLM

## COMPREHENSIONS
[f(x) for i in 1:n]

# Some alternatives for bigger, multiline "comprehension":
# You can use map
map(1:n) do x
    blah
    blah
    blah
    f(x)
end

# You can use begin end expression
[   begin
               y = x + 1
               y^2
           end
           for x = 1 : 5 ]
5-element Array{Int64,1}:
  4
  9
 16
 25
 36

data=DataFrame(x=[0.9,2.0,3.0],y=[2.0,4.0,6.5])
ols=glm(@formula(y ~ x), data, Normal(), IdentityLink())
c=coef(ols)
#  y=c[1]+c[2]*x
#  y=-0.0370091 + 2.13746*x
plot(data,x=:x,y=:y,Geom.point,intercept=[c[1]],slope=[c[2]],Geom.abline(color="red",style=:dash))

using GLM, DataFrames
data = DataFrame(x1=[0.9,2.0,3.0,3.9],x2=[0.291,0.210,0.105,0.021],y=[2.0,4.0,6.5,7.9])
reg = glm(@formula(y ~ x1 + x2), data, Normal(), IdentityLink())
@show(coef(reg))
stderror(reg)
predict(reg)


text:
# Multiple dispatch
#https://www.juliabox.com/notebook/notebooks/tutorials/intro-to-julia/10.%20Multiple%20dispatch.ipynb

using InteractiveUtils
versioninfo()

#=
In this script we'll explore multiple dispatch, which is as key
feature of Julia.

Multiple dispatch makes software generic and fast!

Starting with the familiar
=#

f(x) = x^2

#=
We could declare functions in Julia without giving Julia any
information about the types of the input arguments that
function will receive.
And then Julia will determine on its own which input argument
types make sense and which do not.
=#

@show f(10)
@show f([1 2; 3 4])

#=
However, we also have the option to tell Julia explicitly what types
our input arguments are allowed to have.
=#

foo(x::String, y::String) = println("My inputs x=$x and y=$y are both strings")
foo("Hello","Hi")
foo(x::Int, y::Int) = println("My inputs x=$x and y=$y are both integers")
foo(3, 4) #Now function `foo` works on integers!
foo("Hello","World") # But look!, `foo` also still works when `x` and `y` are strings!

@show methods(foo)
#=
So, we now call `foo` function on integers or strings. When you call
`foo` on a particular set of arguments, Julia will infer the types of the inputs and dispatch
the appropriate method. This is multiple dispatch
=#

println(@which foo(3,4))

#=
And we can continue to add other methods to our generic function `foo`.
Let's add one that takes the abstract type `Number`, which includes subtypes such as
`Int`, `Float64`, and other objects you would think of as numbers:
=#

foo(x::Number, y::Number) = println("My inputs x=$x and y=$y are both numbers!")
foo(3,4)
foo(3.0,4.0)

#=
We can also add a fallback, duck-typed method for `foo` that takes inputs of any type:
=#

foo(x,y)=println("I accept inputs of any type!!!",typeof(x),typeof(y))

v = rand(3)
foo(v,v)

#= output
Julia Version 0.7.0
Commit a4cb80f3ed (2018-08-08 06:46 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin14.5.0)
  CPU: Intel(R) Core(TM) M-5Y51 CPU @ 1.10GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.0 (ORCJIT, broadwell)
f(10) = 100
f([1 2; 3 4]) = [7 10; 15 22]
My inputs x=Hello and y=Hi are both strings
My inputs x=3 and y=4 are both integers
My inputs x=Hello and y=World are both strings
methods(foo) = # 2 methods for generic function "foo":
[1] foo(x::Int64, y::Int64) in Main at /Users/terasakisatoshi/work/juliaLang/officialTutorial/multiple_dispatch.jl:37
[2] foo(x::String, y::String) in Main at /Users/terasakisatoshi/work/juliaLang/officialTutorial/multiple_dispatch.jl:35
foo(x::Int64, y::Int64) in Main at /Users/terasakisatoshi/work/juliaLang/officialTutorial/multiple_dispatch.jl:37
My inputs x=3 and y=4 are both integers
My inputs x=3.0 and y=4.0 are both numbers!
I accept inputs of any type!!!Array{Float64,1}Array{Float64,1}
=#

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using DashBootstrapComponents

items = [
    dbc_dropdownmenuitem("First"),
    dbc_dropdownmenuitem(divider=true),
    dbc_dropdownmenuitem("Second"),
];

dropdown = dbc_row(
    [
        dbc_col(
            dbc_dropdownmenu(
                label="Dropdown (default)",
                children=items,
                direction="down",
            ),
            width="auto",
        ),
        dbc_col(
            dbc_dropdownmenu(label="Dropleft", children=items, direction="left"),
            width="auto",
        ),
        dbc_col(
            dbc_dropdownmenu(label="Dropright", children=items, direction="right"),
            width="auto",
        ),
        dbc_col(
            dbc_dropdownmenu(label="Dropup", children=items, direction="up"),
            width="auto",
        ),
    ],
    justify="between",
);

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:
@testset "Table" begin

@testset "Constructors and macros" begin
    @test @inferred(Table{(:A, :B)}(([1,2,3], [2.0,4.0,6.0]))) == Table{(:A, :B), Tuple{Vector{Int},Vector{Float64}}}(([1,2,3],[2.0,4.0,6.0]))
    @test @inferred(Table{(:A, :B), Tuple{Vector{Int}, Vector{Float64}}}()) == Table{(:A, :B), Tuple{Vector{Int},Vector{Float64}}}((Int[],Float64[]))

    @test @inferred(Table{(:A, :B)}(([1,2,3],[2.0,4.0,6.0]), Val{false})) == Table{(:A, :B), Tuple{Vector{Int},Vector{Float64}}}(([1,2,3],[2.0,4.0,6.0]), Val{true})

    @test @Table(A=[1,2,3], B=[2.0,4.0,6.0])  == Table{(:A, :B), Tuple{Vector{Int},Vector{Float64}}}(([1,2,3],[2.0,4.0,6.0]))
    @test (@Table A=[1,2,3] B=[2.0,4.0,6.0]) == Table{(:A, :B), Tuple{Vector{Int},Vector{Float64}}}(([1,2,3],[2.0,4.0,6.0]))
    @test @Table(A::Vector{Int}=[1,2,3], B::Vector{Float64}=[2.0,4.0,6.0]) == Table{(:A, :B), Tuple{Vector{Int},Vector{Float64}}}(([1,2,3],[2.0,4.0,6.0]))
    @test (@Table A::Vector{Int}=[1,2,3] B::Vector{Float64}=[2.0,4.0,6.0]) == Table{(:A, :B), Tuple{Vector{Int},Vector{Float64}}}(([1,2,3],[2.0,4.0,6.0]))
end

@testset "Introspection" begin
    t = @Table(A=[1,2,3],B=[2.0,4.0,6.0])

    @test (show(t); println(); true)

    @test names(t) == (:A,:B)
    @test names(typeof(t)) == (:A,:B)
    @test names(Table{(:A,:B)}) == (:A,:B)
    @test eltypes(t) == Tuple{Int,Float64}
    @test eltypes(typeof(t)) == Tuple{Int,Float64}
    @test storagetypes(t) == Tuple{Vector{Int},Vector{Float64}}
    @test storagetypes(typeof(t)) == Tuple{Vector{Int},Vector{Float64}}

    @test @inferred(rename(t, Val{(:A_new, :B_new)})) == @Table(A_new=[1,2,3],B_new=[2.0,4.0,6.0])
    @test @inferred(rename(t, Val{:A}, Val{:A_new})) == @Table(A_new=[1,2,3],B=[2.0,4.0,6.0])

    @test length(t) == 3
    @test ncol(t) == 2
    @test nrow(t) == 3
    @test size(t) == (3,)
    @test size(t,1) == 3
    @test ndims(t) == 1
    @test isempty(t) == false
    @test endof(t) == 3

    @test @Table(A=[1,2,3],B=[2.0,4.0,6.0]) == @Table(B=[2.0,4.0,6.0],A=[1,2,3])
    @test @Table(A=[1,2,3],B=[2.0,4.0,6.0]) != @Table(B=[2.0,4.0,6.0],A_new=[1,2,3])
    @test @Table(A=[1,2,3],B=[2.0,4.0,6.0]) != @Table(B=[2.0,4.0,6.0],A=[1,2,4])

    @test copy(t) == t
end

@testset "Accessing and iterating rows" begin
    # push/pop/append with data.

    @test (t=@Table(A=[1,2,3],B=[2.0,4.0,6.0]); empty!(t); isempty(t) == true )

    @test first(@Table(A=[1,2,3],B=[2.0,4.0,6.0])) == @Row(A=1,B=2.0)
    @test next(@Table(A=[1,2,3],B=[2.0,4.0,6.0]),1) == (@Row(A=1,B=2.0),2)
    @test last(@Table(A=[1,2,3],B=[2.0,4.0,6.0])) == @Row(A=3,B=6.0)

    @test (t=@Table(A=[1,2,3],B=[2.0,4.0,6.0]); @inferred(getindex(t, 1)) == @Row(A=1,B=2.0) )
    @test (t=@Table(A=[1,2,3],B=[2.0,4.0,6.0]); @inferred(getindex(t, :)) == @Table(A=[1,2,3],B=[2.0,4.0,6.0]) )
    @test (t=@Table(A=[1,2,3],B=[2.0,4.0,6.0]); @inferred(getindex(t, 3:-1:1)) == @Table(A=[3,2,1],B=[6.0,4.0,2.0]) )

    @test (t=@Table(A=[1,2,3],B=[2.0,4.0,6.0]); t[1] = @Row(A=4,B=0.0); t == @Table(A=[4,2,3],B=[0.0,4.0,6.0]) )
    @test (t=@Table(A=[1,2,3],B=[2.0,4.0,6.0]); t[1] = @Row(B=0.0,A=4); t == @Table(A=[4,2,3],B=[0.0,4.0,6.0]) )
    @test (t=@Table(A=[1,2,3],B=[2.0,4.0,6.0]); t[1:2] = t[2:3]; t == @Table(A=[2,3,3],B=[4.0,6.0,6.0]) )

    # unsafe get/setindex

    @test (t=@Table(A=[1,2],B=[2.0,4.0]); push!(t,@Row(A=3,B=6.0)); t == @Table(A=[1,2,3],B=[2.0,4.0,6.0]) )
    @test (t=@Table(A=[1,2,3],B=[2.0,4.0,6.0]); pop!(t) == @Row(A=3,B=6.0) && t == @Table(A=[1,2],B=[2.0,4.0]) )
    @test (t=@Table(A=[1,2,3],B=[2.0,4.0,6.0]); append!(t,t); t == @Table(A=[1,2,3,1,2,3],B=[2.0,4.0,6.0,2.0,4.0,6.0]) )
    @test (t=@Table(A=[1,2,3],B=[2.0,4.0,6.0]); TypedTables.head(t) == t )
    @test (t=@Table(A=[1,2,3],B=[2.0,4.0,6.0]); TypedTables.tail(t) == t )
end

@testset "Accessing columns and cells" begin
    @test (t=@Table(A=[1,2,3],B=[2.0,4.0,6.0]); @inferred(getindex(t, Val{:A})) == [1,2,3] )
    @test (t=@Table(A=[1,2,3],B=[2.0,4.0,6.0]); @inferred(getindex(t, Val{(:A,)})) == @Table(A=[1,2,3]) )

    # TODO Simultaneous indexing row and column
    # Single Row
    @test (t=@Table(A=[1,2,3],B=[2.0,4.0,6.0]); @inferred(getindex(t, 1, Val{:A})) == 1 )
    @test (t=@Table(A=[1,2,3],B=[2.0,4.0,6.0]); @inferred(getindex(t, 1, Val{(:A,)})) == @Row(A=1) )

    # Multiple rows
    @test (t=@Table(A=[1,2,3],B=[2.0,4.0,6.0]); @inferred(getindex(t, :, Val{:A})) == [1,2,3] )
    @test (t=@Table(A=[1,2,3],B=[2.0,4.0,6.0]); @inferred(getindex(t, :, Val{(:A,)})) == @Table(A=[1,2,3]) )

    @test (t=@Table(A=[1,2,3],B=[2.0,4.0,6.0]); @inferred(getindex(t, 1, :)) == @Row(A=1,B=2.0) )
    @test (t=@Table(A=[1,2,3],B=[2.0,4.0,6.0]); @inferred(getindex(t, :, :)) == @Table(A=[1,2,3],B=[2.0,4.0,6.0]) )

end

@testset "Concatenation" begin
    @test (r=@Row(A=1,B=2.0); t = @inferred(vcat(r)); t == @Table(A=[1],B=[2.0]) )
    @test (r=@Row(A=1,B=2.0); t = @inferred(vcat(r,r)); t == @Table(A=[1,1],B=[2.0,2.0]) )
    @test (r=@Row(A=1,B=2.0); t = (vcat(r,r,r)); t == @Table(A=[1,1,1],B=[2.0,2.0,2.0]) )

    @test (r=@Row(A=1,B=2.0); t = vcat(r); t2 = (vcat(r,t)); t2 == @Table(A=[1,1],B=[2.0,2.0]) )
    @test (r=@Row(A=1,B=2.0); t = vcat(r); t2 = (vcat(t,r)); t2 == @Table(A=[1,1],B=[2.0,2.0]) )
    @test (r=@Row(A=1,B=2.0); t = vcat(r); t2 = @inferred(vcat(t,t)); t2 == @Table(A=[1,1],B=[2.0,2.0]) )
    @test (r=@Row(A=1,B=2.0); t = vcat(r); t2 = @inferred(vcat(t,t,t)); t2 == @Table(A=[1,1,1],B=[2.0,2.0,2.0]) )

    @test (c=@Column(A=[1,2,3]); t = @inferred(hcat(c)); t == @Table(A=[1,2,3]) )
    @test (c1=@Column(A=[1,2,3]); c2=@Column(B=[2.0,4.0,6.0]); t = @inferred(hcat(c1,c2)); t == @Table(A=[1,2,3], B=[2.0,4.0,6.0]) )
    @test (c1=@Column(A=[1,2,3]); c2=@Column(B=[2.0,4.0,6.0]); c3=@Column(C=[true,false,false]); t = @inferred(hcat(c1,c2,c3)); t == @Table(A=[1,2,3],B=[2.0,4.0,6.0],C=[true,false,false]) )

    @test (t=@Table(A=[1,2,3]); t2 = @inferred(hcat(t)); t2 == @Table(A=[1,2,3]) )
    @test (c1=@Column(A=[1,2,3]); c2=@Column(B=[2.0,4.0,6.0]); t = @inferred(hcat(hcat(c1),c2)); t == @Table(A=[1,2,3], B=[2.0,4.0,6.0]) )
    @test (c1=@Column(A=[1,2,3]); c2=@Column(B=[2.0,4.0,6.0]); t = @inferred(hcat(c1,hcat(c2))); t == @Table(A=[1,2,3], B=[2.0,4.0,6.0]) )
    @test (c1=@Column(A=[1,2,3]); c2=@Column(B=[2.0,4.0,6.0]); t = @inferred(hcat(hcat(c1),hcat(c2))); t == @Table(A=[1,2,3], B=[2.0,4.0,6.0]) )
    @test (c1=@Column(A=[1,2,3]); c2=@Column(B=[2.0,4.0,6.0]); c3=@Column(C=[true,false,false]); t = @inferred(hcat(hcat(c1),c2,c3)); t == @Table(A=[1,2,3],B=[2.0,4.0,6.0],C=[true,false,false]) )
    @test (c1=@Column(A=[1,2,3]); c2=@Column(B=[2.0,4.0,6.0]); c3=@Column(C=[true,false,false]); t = @inferred(hcat(c1,hcat(c2),c3)); t == @Table(A=[1,2,3],B=[2.0,4.0,6.0],C=[true,false,false]) )
    @test (c1=@Column(A=[1,2,3]); c2=@Column(B=[2.0,4.0,6.0]); c3=@Column(C=[true,false,false]); t = @inferred(hcat(hcat(c1),hcat(c2),c3)); t == @Table(A=[1,2,3],B=[2.0,4.0,6.0],C=[true,false,false]) )
end

end


text:
using HorizonSideRobots
include("../basics.jl")

md"""
        ДАНО: Робот - в произвольной клетке ограниченного прямоугольного поля (без внутренних перегородок)

        РЕЗУЛЬТАТ: Робот - в исходном положении, в клетке с роботом стоит маркер, и все остальные клетки поля промаркированы в шахматном порядке
"""

 """Решение задачи 7"""
 function solve_7!(r)
    (M, N), (I, J) = get_field_size_and_position!(r)
    visited = Set()
    cells = Set()

    for i in 1:M
        for j in 1:N
            if (I % 2 != i % 2 && J % 2 != j % 2) || (I % 2 == i % 2 && J % 2 == j % 2)
                push!(cells, (i, j))
            end
        end
    end

    mark_some_cells_with_dfs!(r, visited, (I, J), false, cells)

    show(r)
end


text:
# Adaptive Gauss-Kronrod quadrature routines (arbitrary precision),
# originally written and contributed to Julia by Steven G. Johnson, 2013.
#
# This file was formerly a part of Julia. License is MIT: http://julialang.org/license

"""
The `QuadGK` module implements 1d numerical integration by an adaptive Gauss–Kronrod algorithm.

The package provides three functions: `quadgk`, `gauss`, and `kronrod`.

* `quadgk` performs adaptive integration over arbitrary intervals
* `gauss` computes Gaussian quadrature points and weights for integrating over the interval [-1, 1]
* `kronrod` computes Kronrod points, weights, and embedded Gaussian quadrature weights for integrating over [-1, 1].

Typical usage looks like:
```
using QuadGK
integral, err = quadgk(x -> exp(-x^2), 0, 1, rtol=1e-8)
```
which computes the integral of exp(–x²) from x=0 to x=1 to a relative tolerance of 10⁻⁸,
and returns the approximate `integral = 0.746824132812427` and error estimate
`err = 7.887024366937112e-13`.
"""
module QuadGK

export quadgk, gauss, kronrod

using DataStructures, LinearAlgebra
import Base.Order.Reverse

include("gausskronrod.jl")
include("evalrule.jl")
include("adapt.jl")
include("weightedgauss.jl")

end # module QuadGK


text:
using EzXML
using Revise

includet("..\\src\\purecompbase.jl")
includet("..\\src\\propsbase.jl")
includet("..\\src\\readDB.jl")

fileName = String("data\\chemsep1.xml")

dbRoot = getDBRoot(fileName)

idx = 1
tmp = EzXML.findfirst("//compounds/compound[" * string(idx) *"]", dbRoot)
props = nodes(tmp)


text:
#=
test_earley:
- Julia version: 1.1.0
- Author: bramb
- Date: 2019-03-08
=#
using Test
using HyperGraphTools.Earley

@testset "Earley 1" begin
    g = Grammar()

    N = Rule("N", [Production("time"), Production("flight"), Production("banana"), Production("flies"), Production("boy"), Production("telescope")])
    D = Rule("D", [Production("the"), Production("a"), Production("an")])
    V = Rule("V", [Production("book"), Production("eat"), Production("sleep"), Production("saw")])
    P = Rule("P", [Production("with"), Production("in"), Production("on"), Production("at"), Production("through")])

    PP = Rule("PP")
    NP = Rule("NP", [Production([D, N]), Production("john"), Production("houston")])
    push!(NP,Production([NP, PP]))
    push!(PP,Production([P, NP]))

    VP = Rule("VP", Production([V, NP]))
    push!(VP,Production([VP, PP]))
    S = Rule("S", [Production([NP, VP]), Production([VP])])
#     @show(S)

    push!(g,S)
    push!(g,N)
    push!(g,D)
    push!(g,V)
    push!(g,P)
    push!(g,PP)
    push!(g,NP)
    push!(g,VP)
#     @show(g)

    grammar_string = str(g)
    expected = """
    S -> NP VP | VP
    PP -> P NP
    NP -> D N | "john" | "houston" | NP PP
    P -> "with" | "in" | "on" | "at" | "through"
    N -> "time" | "flight" | "banana" | "flies" | "boy" | "telescope"
    V -> "book" | "eat" | "sleep" | "saw"
    VP -> V NP | VP PP
    D -> "the" | "a" | "an"
    """

    @test grammar_string == expected
    println(grammar_string)

    tokens = ["book", "the", "flight", "through", "houston"]
    state = earley_parse(tokens,g)
    println(str(state))

    tokens = ["john", "saw", "the", "boy", "with", "the", "telescope"]
    state = earley_parse(tokens,g)
    println(str(state))
end

@testset "Earley 2" begin
    g = Grammar()

    SYM = Rule("SYM", Production("a"))
    OP = Rule("OP", Production("+"))
    EXPR = Rule("EXPR", Production([SYM]))
    push!(EXPR,Production([EXPR, OP, EXPR]))
    S = Rule("S", [Production([EXPR])])

    push!(g,S)
    push!(g,SYM)
    push!(g,OP)
    push!(g,EXPR)

    grammar_string = str(g)
    expected = """
    EXPR -> SYM | EXPR OP EXPR
    S -> EXPR
    SYM -> "a"
    OP -> "+"
    """

    @test grammar_string == expected
    println(grammar_string)

    tokens = ["a", "+", "a", "+", "a", "+", "a"]
    state = earley_parse(tokens,g)
    println(str(state))
end

max_stars_repo_name: brambg/HyperGraphTools

max_issues_repo_name: brambg/HyperGraphTools

max_forks_repo_name: brambg/HyperGraphTools


text:
using URIs: URI

const GOOGLE_STORAGE_API = "https://storage.googleapis.com"
const GOOGLE_STORAGE_REST_API = GOOGLE_STORAGE_API * "/storage/v1"
const GOOGLE_STORAGE_CREDENTIALS = Dict{String,String}()

"""
    Zarr.gcs_credentials(user_project,access_token,token_type)

Set the user project, access token and and token type for the Google Cloud
Store.
"""
function gcs_credentials(user_project,access_token,token_type)
  GOOGLE_STORAGE_CREDENTIALS["user_project"] = user_project
  GOOGLE_STORAGE_CREDENTIALS["access_token"] = access_token
  GOOGLE_STORAGE_CREDENTIALS["token_type"] = token_type
  nothing
end

"""
    Zarr.gcs_credentials(; metadata_url = "http://metadata.google.internal/computeMetadata/v1/")

Set (or renew) the user project, access token and and token type for the Google
Cloud Store from the Meatadata server (assuming the function is executed from
the Google Cloud).
For some data sets, the error message "Bucket is requester pays bucket but no
user project provided" is returned if the credentials are not provided.
"""
function gcs_credentials(;metadata_url = "http://metadata.google.internal/computeMetadata/v1/")
  headers = Dict("Metadata-Flavor" => "Google")

  url = joinpath(metadata_url,"project","project-id")
  user_project = String(HTTP.get(url, headers=headers).body)

  url = joinpath(metadata_url,"instance","service-accounts","default","token")
  auth = JSON.parse(String(HTTP.get(url,headers).body));

  gcs_credentials(user_project,auth["access_token"],auth["token_type"])
end

function _gcs_request_headers()
  headers = Dict{String,String}()
  if haskey(GOOGLE_STORAGE_CREDENTIALS,"user_project")
    headers["x-goog-user-project"] = GOOGLE_STORAGE_CREDENTIALS["user_project"]
  end

  if haskey(GOOGLE_STORAGE_CREDENTIALS,"token_type") &&
    haskey(GOOGLE_STORAGE_CREDENTIALS,"access_token")

    headers["Authorization"] = string(
      GOOGLE_STORAGE_CREDENTIALS["token_type"]," ",
      GOOGLE_STORAGE_CREDENTIALS["access_token"])
  end

  return headers
end

struct GCStore <: AbstractStore
  bucket::String

  function GCStore(url::String)
    uri = URI(url)

    if uri.scheme == "gs"
      bucket = uri.host
    else
      parts = split(uri.path,'/',limit=3)
      bucket = parts[2]
    end
    @debug "GCS bucket: $bucket"
    new(bucket)
  end
end

function Base.getindex(s::GCStore, k::String)
  url = string(GOOGLE_STORAGE_API,"/",s.bucket,"/",k)
  headers = _gcs_request_headers()
  r = HTTP.request("GET",url,headers,status_exception = false)
  if r.status >= 300
    if r.status == 404
      @debug "get: $url: not found"
      nothing
    else
      error("Error connecting to $url :", String(r.body))
    end
  else
    @debug "get: $url"
    r.body
  end
end

function cloud_list_objects(s::GCStore,p)
  prefix = (isempty(p) || endswith(p,"/")) ? p : string(p,"/")

  url = string(GOOGLE_STORAGE_REST_API, "/b/", s.bucket, "/o")

  @debug "call: $url"
  headers = _gcs_request_headers()
  params = Dict("prefix" => prefix, "delimiter" => "/")
  r = JSON.parse(String(HTTP.get(url,headers,
                                 query = params).body))

  return r
end

function storagesize(s::GCStore,p)
  r = cloud_list_objects(s,p)
  items = r["items"]
  datafiles = filter(entry -> !any(filename -> endswith(entry["name"], filename), [".zattrs",".zarray",".zgroup"]), items)
  if isempty(datafiles)
    0
  else
    sum(datafiles) do f
      parse(Int, f["size"])
    end
  end
end

function subkeys(s::GCStore, p)
  r = cloud_list_objects(s, p)
  keys = map(item -> String(split(item["name"],'/')[end]),  r["items"])
  return keys
end

function subdirs(s::GCStore, p)
  r = cloud_list_objects(s,p)
  dirs = map(prefix -> String(split(prefix,'/')[end-1]), r["prefixes"])
  return dirs
end

pushfirst!(storageregexlist,r"^https://storage.googleapis.com"=>GCStore)
pushfirst!(storageregexlist,r"^http://storage.googleapis.com"=>GCStore)
push!(storageregexlist,r"^gs://"=>GCStore)

function storefromstring(::Type{<:GCStore}, url,_)
  uri = URI(url)
  if uri.scheme == "gs"
    p = lstrip(uri.path,'/')
  else
    parts = split(uri.path,'/',limit=2, keepempty=false)
    p = (length(parts) == 2 ? parts[2] : "")
  end

  @debug "path: $p"
  return GCStore(url),p
end


text:
abstract type OrdinaryDiffEqCache <: DiffEqBase.DECache end
abstract type OrdinaryDiffEqConstantCache <: OrdinaryDiffEqCache end
abstract type OrdinaryDiffEqMutableCache <: OrdinaryDiffEqCache end
struct ODEEmptyCache <: OrdinaryDiffEqConstantCache end
struct ODEChunkCache{CS} <: OrdinaryDiffEqConstantCache end

@cache mutable struct CompositeCache{T,F} <: OrdinaryDiffEqCache
  caches::T
  choice_function::F
  current::Int
end

function alg_cache(alg::CompositeAlgorithm,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false})
  caches = map(alg.algs) do x
    alg_cache(x,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,Val(false))
  end
  CompositeCache(caches,alg.choice_function,1)
end

function alg_cache(alg::CompositeAlgorithm,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true})
  caches = map(alg.algs) do x
    alg_cache(x,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,Val(true))
  end
  CompositeCache(caches,alg.choice_function,1)
end

alg_cache(alg::OrdinaryDiffEqAlgorithm,prob,callback::F) where {F} = ODEEmptyCache()

@cache struct FunctionMapCache{uType,rateType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::rateType
end

function alg_cache(alg::FunctionMap,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true})
  FunctionMapCache(u,uprev,FunctionMap_scale_by_time(alg) ? rate_prototype : similar(u))
end

struct FunctionMapConstantCache <: OrdinaryDiffEqConstantCache end

alg_cache(alg::FunctionMap,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false}) = FunctionMapConstantCache()

@cache struct ExplicitRKCache{uType,rateType,uNoUnitsType,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  utilde::rateType
  atmp::uNoUnitsType
  fsalfirst::rateType
  fsallast::rateType
  kk::Vector{rateType}
  tab::TabType
end

function alg_cache(alg::ExplicitRK,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true})
  kk = Vector{typeof(rate_prototype)}(undef, 0)
  for i = 1:alg.tableau.stages
    push!(kk,zero(rate_prototype))
  end
  fsalfirst = kk[1]
  if isfsal(alg.tableau)
    fsallast = kk[end]
  else
    fsallast = zero(rate_prototype)
  end
  utilde = zero(rate_prototype)
  tmp = similar(u)
  atmp = similar(u,uEltypeNoUnits)
  tab = ExplicitRKConstantCache(alg.tableau,rate_prototype)
  ExplicitRKCache(u,uprev,tmp,utilde,atmp,fsalfirst,fsallast,kk,tab)
end

struct ExplicitRKConstantCache{MType,VType,KType} <: OrdinaryDiffEqConstantCache
  A::MType
  c::VType
  α::VType
  αEEst::VType
  stages::Int
  kk::KType
end

function ExplicitRKConstantCache(tableau,rate_prototype)
  @unpack A,c,α,αEEst,stages = tableau
  A = A' # Transpose A to column major looping
  kk = Array{typeof(rate_prototype)}(undef, stages) # Not ks since that's for integrator.opts.dense
  ExplicitRKConstantCache(A,c,α,αEEst,stages,kk)
end

alg_cache(alg::ExplicitRK,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false}) = ExplicitRKConstantCache(alg.tableau,rate_prototype)

get_chunksize(cache::DiffEqBase.DECache) = error("This cache does not have a chunksize.")
get_chunksize(cache::ODEChunkCache{CS}) where {CS} = CS


text:
# __BEGIN_LICENSE__
#
# ThreeDeconv.jl
#
# Copyright (c) 2018, Stanford University
#
# All rights reserved.
#
# Redistribution and use in source and binary forms for academic and other
# non-commercial purposes with or without modification, are permitted provided
# that the following conditions are met:
#
# * Redistributions of source code, including modified source code, must retain
#   the above copyright notice, this list of conditions and the following
#   disclaimer.
#
# * Redistributions in binary form or a modified form of the source code must
#   reproduce the above copyright notice, this list of conditions and the
#   following disclaimer in the documentation and/or other materials provided with
#   the distribution.
#
# * Neither the name of The Leland Stanford Junior University, any of its
#   trademarks, the names of its employees, nor contributors to the source code
#   may be used to endorse or promote products derived from this software without
#   specific prior written permission.
#
# * Where a modified version of the source code is redistributed publicly in
#   source or binary forms, the modified source code must be published in a freely
#   accessible manner, or otherwise redistributed at no charge to anyone
#   requesting a copy of the modified source code, subject to the same terms as
#   this agreement.
#
# THIS SOFTWARE IS PROVIDED BY THE TRUSTEES OF THE LELAND STANFORD JUNIOR
# UNIVERSITY "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE LELAND STANFORD JUNIOR
# UNIVERSITY OR ITS TRUSTEES BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
# BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
# IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
# OF SUCH DAMAGE.
#
# __END_LICENSE__

using Printf

include("foi.jl")

"""
noise estimation
"""
function noise_estimation(img::AbstractArray;
                          maxnum_pairs::Int=300,
                          verbose::Bool=false)
    @assert sizeof(img) <= 10^9 * 8 # Prevent input image size more than 1GB

    offset = -minimum(img)
    subimg = img .+ offset
    scale = 1 / (maximum(subimg) * 1.1)
    scaled_img = Float32.(subimg.*scale)

    scaled_a, scaled_b = foi_noiseestimation(scaled_img, maxnum_pairs=maxnum_pairs, verbose=verbose)

    # Rescale back the estimated parameters.
    a = scaled_a / scale
    b = scaled_b / scale^2
    @printf "Estimated noise parameters are a = %4.3f, b = %4.3f.\n" a b

    γ_est = a
    σsq_est = a * offset + b

    if σsq_est < .0
        @printf "Noise estimation might have failed. σ^2 = %4.3f.\n
                 Proceeding by setting σ^2 = .0 to assure the variance's
                 nonnegativity.\n" σsq_est
        σsq_est = .0 + 1e-10
    end
    σ_est = sqrt(σsq_est)
    println("Estimated parameters for Poisson-Gaussian noise model")
    println("γ: $γ_est, σ: $σ_est")
    return γ_est, σ_est
end


text:
using LCPSim
using Base.Test
using RigidBodyDynamics
using RigidBodyDynamics: Bounds
using StaticArrays: SVector
using Gurobi: GurobiSolver
using Rotations: RotY
using MeshCat
using MeshCatMechanisms
using CoordinateTransformations: Translation

function load_model(urdf, θ)
    mechanism = parse_urdf(Float64, urdf)
    env = LCPSim.parse_contacts(mechanism, urdf, 1.0, :xz)
    planar_joint = findjoint(mechanism, "plane_to_box")
    planar_joint.position_bounds .= Bounds(-10, 10)
    planar_joint.velocity_bounds .= Bounds(-1000, 1000)
    state = MechanismState(mechanism)
    set_configuration!(state, planar_joint, [-1, -1.5, θ])
    mechanism, env, state
end

@testset "bricks on inclined planes" begin
    vis = Visualizer()
    if !haskey(ENV, "CI")
        open(vis)
        wait(vis)
    end

    # A brick on an incline which is slightly too shallow, so the brick does not
    # slide at all
    urdf = joinpath(@__DIR__, "urdf", "shallow_incline.urdf")
    mechanism, env, x1 = load_model(urdf, -π/4 - 0.1)
    mv1 = MechanismVisualizer(mechanism, URDFVisuals(urdf), vis[:stick][:model])
    Δt = 0.02
    N = 100
    controller = LCPSim.passive_controller()

    q0 = copy(configuration(x1))
    results_stick = LCPSim.simulate(x1, controller, env, Δt, N, GurobiSolver(Gurobi.Env(), OutputFlag=0))
    @testset "sticking" begin
        @test length(results_stick) == N
        for i in 50:length(results_stick)
            @test norm(configuration(results_stick[i].state) .- configuration(results_stick[i - 1].state)) <= 1e-5
            @test norm(velocity(results_stick[i].state)) <= 1e-5
        end
    end

    # A slightly steeper incline causes the brick to begin to slide
    urdf = joinpath(@__DIR__, "urdf", "steep_incline.urdf")
    mechanism, env, x2 = load_model(urdf, -π/4 + 0.1)
    mv2 = MechanismVisualizer(mechanism, URDFVisuals(urdf), vis[:slide][:model])
    settransform!(vis[:slide], Translation(0, 1, 0))
    q0 = copy(configuration(x2))
    results_slide = LCPSim.simulate(x2, controller, env, Δt, N, GurobiSolver(Gurobi.Env(), OutputFlag=0))
    @testset "sliding" begin
        @test length(results_slide) == N
        for i in 50:length(results_slide)
            @test norm(configuration(results_slide[i].state) .- configuration(results_slide[i - 1].state)) >= 1e-2
            @test norm(velocity(results_slide[i].state)) >= 0.5
        end
    end

    for i in 1:length(results_stick)
        set_configuration!(mv1, configuration(results_stick[i].state))
        set_configuration!(mv2, configuration(results_slide[i].state))
        sleep(Δt)
    end

end