text:
module PartitionedPoisson

using Gridap
using Gridap.Geometry
using Gridap.FESpaces
using LinearAlgebra
import PartitionedArrays
const PArrays = PartitionedArrays
using MPI
using FileIO

export poisson

u_2d(x) = x[1]+x[2]
u_3d(x) = x[1]+x[2]+x[3]

function poisson(;
  mode::Symbol=:seq,
  nc::Tuple=(10,10,10),
  np::Tuple=(2,2,2),
  nr::Int=1,
  title::String="seq",
  verbose::Bool=false)

  if ! (mode in (:seq,:mpi))
    throw(ArgumentError("mode=$mode is not a valid value for the kw-argument."))
  end

  if mode == :mpi && ! MPI.Initialized()
    MPI.Init()
  end

  backend = mode == :seq ? PArrays.sequential : PArrays.mpi
  parts = PArrays.get_part_ids(backend,np)
  for ir in 1:nr
    str_r = lpad(ir,ceil(Int,log10(nr)),'0')
    title_r = "$(title)_ir$(str_r)"
    _poisson(parts,nc,title_r,ir,verbose)
  end
end

function _poisson(parts,nc,title,ir,verbose)

  domain_2d = (0,1,0,1)
  domain_3d = (0,1,0,1,0,1)
  np = size(parts)

  domain = length(nc) == 3 ? domain_3d : domain_2d
  u = length(nc) == 3 ? u_3d : u_2d
  order = 1

  t = PArrays.PTimer(parts,verbose=verbose)
  PArrays.tic!(t)

  # Partition of the Cartesian ids
  # with ghost layer
  cell_gcis = PArrays.PCartesianIndices(parts,nc,PArrays.with_ghost)
  PArrays.toc!(t,"cell_gcis")

  # Partitioned range of cells
  # with ghost layer
  cell_range = PArrays.PRange(parts,nc,PArrays.with_ghost)
  neighbors = cell_range.exchanger.parts_snd
  ngcells = length(cell_range)
  PArrays.toc!(t,"cell_range")

  # Local discrete models
  model = PArrays.map_parts(cell_gcis) do gcis
    cmin = first(gcis)
    cmax = last(gcis)
    desc = CartesianDescriptor(domain,nc)
    CartesianDiscreteModel(desc,cmin,cmax)
  end
  PArrays.toc!(t,"model")

  # Local FE spaces
  U, V = PArrays.map_parts(model) do model
    reffe = ReferenceFE(lagrangian,Float64,order)
    V = TestFESpace(model,reffe;dirichlet_tags="boundary")
    U = TrialFESpace(u,V)
    U, V
  end
  PArrays.toc!(t,"U, V")

  # Cell-wise local dofs
  cell_to_ldofs, nldofs = PArrays.map_parts(V) do V
    get_cell_dof_ids(V), num_free_dofs(V)
  end
  PArrays.toc!(t,"cell_to_ldofs, nldofs")

  # Find and count number owned dofs
  ldof_to_part, nodofs = PArrays.map_parts(
    cell_range.partition,cell_to_ldofs,nldofs) do partition,cell_to_ldofs,nldofs

    ldof_to_part = fill(Int32(0),nldofs)
    cache = array_cache(cell_to_ldofs)
    for cell in 1:length(cell_to_ldofs)
      owner = partition.lid_to_part[cell]
      ldofs = getindex!(cache,cell_to_ldofs,cell)
      for ldof in ldofs
        if ldof>0
          #TODO this simple approach concentrates dofs
          # in the last part and creates inbalances
          ldof_to_part[ldof] = max(owner,ldof_to_part[ldof])
        end
      end
    end
    nodofs = count(p->p==partition.part,ldof_to_part)
    ldof_to_part, nodofs
  end
  PArrays.toc!(t,"ldof_to_part, nodofs")

  # Find the global range of owned dofs
  first_gdof, ngdofsplus1 = PArrays.xscan(+,reduce,nodofs,init=1)
  ngdofs = ngdofsplus1 - 1
  PArrays.toc!(t,"first_gdof, ngdofs")

  # Distribute gdofs to owned ones
  ldof_to_gdof = PArrays.map_parts(
    parts,first_gdof,ldof_to_part) do part,first_gdof,ldof_to_part

    offset = first_gdof-1
    ldof_to_gdof = Vector{Int}(undef,length(ldof_to_part))
    odof = 0
    gdof = 0
    for (ldof,owner) in enumerate(ldof_to_part)
      if owner == part
        odof += 1
        ldof_to_gdof[ldof] = odof
      else
        ldof_to_gdof[ldof] = gdof
      end
    end
    for (ldof,owner) in enumerate(ldof_to_part)
      if owner == part
        ldof_to_gdof[ldof] += offset
      end
    end
    ldof_to_gdof
  end
  PArrays.toc!(t,"ldof_to_gdof (owned)")

  # Create cell-wise global dofs
  cell_to_gdofs = PArrays.map_parts(
    parts,
    ldof_to_gdof,cell_to_ldofs,cell_range.partition) do part,
    ldof_to_gdof,cell_to_ldofs,partition

    cache = array_cache(cell_to_ldofs)
    ncells = length(cell_to_ldofs)
    ptrs = Vector{Int32}(undef,ncells+1)
    for cell in 1:ncells
      ldofs = getindex!(cache,cell_to_ldofs,cell)
      ptrs[cell+1] = length(ldofs)
    end
    PArrays.length_to_ptrs!(ptrs)
    ndata = ptrs[end]-1
    data = Vector{Int}(undef,ndata)
    gdof = 0
    for cell in partition.oid_to_lid
      ldofs = getindex!(cache,cell_to_ldofs,cell)
      p = ptrs[cell]-1
      for (i,ldof) in enumerate(ldofs)
        if ldof > 0 && ldof_to_gdof[ldof] != gdof
          data[i+p] = ldof_to_gdof[ldof]
        end
      end
    end
    PArrays.Table(data,ptrs)
  end
  PArrays.toc!(t,"cell_to_gdofs (owned)")

  # Exchange the global dofs
  PArrays.exchange!(cell_to_gdofs,cell_range.exchanger)
  PArrays.toc!(t,"cell_to_gdofs (ghost)")

  # Distribute global dof ids also to ghost
  PArrays.map_parts(
    parts,
    cell_to_ldofs,cell_to_gdofs,ldof_to_gdof,ldof_to_part,cell_range.partition) do part,
    cell_to_ldofs,cell_to_gdofs,ldof_to_gdof,ldof_to_part,partition

    gdof = 0
    cache = array_cache(cell_to_ldofs)
    for cell in partition.hid_to_lid
      ldofs = getindex!(cache,cell_to_ldofs,cell)
      p = cell_to_gdofs.ptrs[cell]-1
      for (i,ldof) in enumerate(ldofs)
        if ldof > 0 && ldof_to_part[ldof] == partition.lid_to_part[cell]
          ldof_to_gdof[ldof] = cell_to_gdofs.data[i+p]
        end
      end
    end
  end
  PArrays.toc!(t,"ldof_to_gdof (ghost)")

  # Setup dof partition
  dof_partition = PArrays.map_parts(parts,ldof_to_gdof,ldof_to_part) do part,ldof_to_gdof,ldof_to_part
    PArrays.IndexSet(part,ldof_to_gdof,ldof_to_part)
  end
  PArrays.toc!(t,"dof_partition")

  # Setup dof exchanger
  dof_exchanger = PArrays.Exchanger(dof_partition,neighbors)
  PArrays.toc!(t,"dof_exchanger")

  # Setup dof range
  dofs = PArrays.PRange(ngdofs,dof_partition,dof_exchanger)
  PArrays.toc!(t,"dofs")

  # Setup Integration (only for owned cells)
  Ω, dΩ = PArrays.map_parts(cell_range.partition,model) do partition, model
    Ω = Triangulation(model,partition.oid_to_lid)
    dΩ = Measure(Ω,2*order)
    Ω, dΩ
  end
  PArrays.toc!(t,"Ω, dΩ")

  # Integrate the coo vectors
  I,J,C,vec = PArrays.map_parts(Ω,dΩ,U,V,ldof_to_gdof) do Ω,dΩ,U,V,ldof_to_gdof
    dv = get_cell_shapefuns(V)
    du = get_cell_shapefuns_trial(U)
    cellmat = ∫( ∇(du)⋅∇(dv) )dΩ
    cellvec = 0
    uhd = zero(U)
    matvecdata = collect_cell_matrix_and_vector(cellmat,cellvec,uhd)
    assem = SparseMatrixAssembler(U,V)
    ncoo = count_matrix_and_vector_nnz_coo(assem,matvecdata)
    I = zeros(Int,ncoo)
    J = zeros(Int,ncoo)
    C = zeros(Float64,ncoo)
    vec = zeros(Float64,num_free_dofs(V))
    fill_matrix_and_vector_coo_numeric!(I,J,C,vec,assem,matvecdata)
    I,J,C,vec
  end
  PArrays.toc!(t,"I,J,C,vec")

  PArrays.map_parts(I,J,ldof_to_gdof) do I,J,ldof_to_gdof
    for i in 1:length(I)
      I[i] = ldof_to_gdof[I[i]]
      J[i] = ldof_to_gdof[J[i]]
    end
  end
  PArrays.toc!(t,"I,J (global)")

  # Find the ghost rows
  hrow_to_hdof = PArrays.touched_hids(dofs,I)
  PArrays.toc!(t,"hrow_to_hdof")

  hrow_to_gid, hrow_to_part = PArrays.map_parts(
    hrow_to_hdof, dof_partition) do hrow_to_hdof, dof_partition
    hrow_to_ldof = view(dof_partition.hid_to_lid,hrow_to_hdof)
    hrow_to_gdof = dof_partition.lid_to_gid[hrow_to_ldof]
    hrow_to_part = dof_partition.lid_to_part[hrow_to_ldof]
    hrow_to_gdof, hrow_to_part
  end
  PArrays.toc!(t,"hrow_to_gid, hrow_to_part")

  # Create the range for rows
  rows = PArrays.PRange(
    parts,ngdofs,nodofs,first_gdof,hrow_to_gid,hrow_to_part,neighbors)
  PArrays.toc!(t,"rows")

  # Move values to the owner part
  # since we have integrated only over owned cells
  PArrays.assemble!(I,J,C,rows)
  PArrays.toc!(t,"I,J,C (assemble!)")

  # Find the ghost cols
  hcol_to_hdof = PArrays.touched_hids(dofs,J)
  PArrays.toc!(t,"hcol_to_hdof")

  hcol_to_gid, hcol_to_part = PArrays.map_parts(
    hcol_to_hdof, dof_partition) do hcol_to_hdof, dof_partition
    hcol_to_ldof = view(dof_partition.hid_to_lid,hcol_to_hdof)
    hcol_to_gdof = dof_partition.lid_to_gid[hcol_to_ldof]
    hcol_to_part = dof_partition.lid_to_part[hcol_to_ldof]
    hcol_to_gdof, hcol_to_part
  end
  PArrays.toc!(t,"hcol_to_gid, hcol_to_part")

  # Create the range for cols
  cols = PArrays.PRange(
    parts,ngdofs,nodofs,first_gdof,hcol_to_gid,hcol_to_part,neighbors)
  PArrays.toc!(t,"cols")

  # Create the sparse matrix
  A_exchanger = PArrays.empty_exchanger(parts)
  A = PArrays.PSparseMatrix(I,J,C,rows,cols,A_exchanger,ids=:global)
  PArrays.toc!(t,"A")

  # Rhs aligned with the FESpace
  dof_values = PArrays.PVector(vec,dofs)
  PArrays.toc!(t,"dof_values")

  # Allocate rhs aligned with the matrix
  b = PArrays.PVector(0.0,rows)
  PArrays.toc!(t,"b (allocate)")

  # Fill rhs
  # TODO A more elegant way of doing this?
  PArrays.map_parts(
    b.values,dof_values.values,rows.partition,dofs.partition,first_gdof) do b1, b2, p1, p2, first_gdof
    offset = first_gdof - 1
    for i in p1.oid_to_lid
      gdof = p1.lid_to_gid[i]
      odof = gdof-offset
      ldof = p2.oid_to_lid[odof]
      b1[i] = b2[ldof]
    end
    for i in p1.hid_to_lid
      gdof = p1.lid_to_gid[i]
      ldof = p2.gid_to_lid[gdof]
      b1[i] = b2[ldof]
    end
  end
  PArrays.toc!(t,"b (fill)")

  # Import and add remote contributions
  PArrays.assemble!(b)
  PArrays.toc!(t,"b (assemble!)")

  # Interpolate exact solution
  # (this is aligned with the FESpace)
  PArrays.map_parts(dof_values.values,U) do dof_values,U
    interpolate!(u,dof_values,U)
  end

  # Allocate solution
  # aligned with the matrix
  x = PArrays.PVector(0.0,cols)
  PArrays.toc!(t,"x (allocate)")

  #### Fill
  #### only needed to fill owned values
  #### since A*x will do the exchange
  x .= dof_values
  PArrays.toc!(t,"x (fill)")

  r = similar(b)
  PArrays.toc!(t,"r (allocate)")

  mul!(r,A,x)
  PArrays.toc!(t,"r (A*x)")

  r .= r .- b
  PArrays.toc!(t,"r (-b)")

  errnorm = norm(r)
  PArrays.toc!(t,"norm(r)")

  display(t)

  PArrays.map_main(t.data) do data
    out = Dict{String,Any}()
    merge!(out,data)
    out["errnorm"] = errnorm
    out["ngdofs"] = ngdofs
    out["ngcells"] = ngcells
    out["nc"] = nc
    out["np"] = np
    out["ir"] = ir
    save("$title.bson",out)
  end

  nothing
end

end # module


text:
using CART
using Test

tests = ["basic"] # the test file names are stored as strings...

println("Running tests:")

for t in tests
    include("$(t).jl") # ... so that they can be evaluated in a loop
end

ext: jl


text:
using Rock
using Base.Test
w = addprocs(3)
@everywhere module KVStore
import Rock
@enum Op Put Append
type Command <: Rock.ExternalCriticalCommand
    op::Op
    k::ASCIIString
    v::ASCIIString
    Command(op, k,v) = new(op, k,v)
end
typealias Snapshot Dict{ASCIIString, ASCIIString}
function transition(px, t)
    if t.op == Put
        SQLite.query(Rock.db(px), "insert into snapshot values (?,?)",[t.k,t.v])
    elseif t.op == Append
        SQLite.query(Rock.db(px),"insert into snapshot values (?, ?) on duplicate key update snapshot set V = V||? where K = ?", [t.k, t.v, t.v, t.k])
    end
end
function mk(port, peers::Array{Tuple{ASCIIString, Int64},1}, path::ASCIIString)
    ctx = Rock.Context(port, peers, transition, "$path/sqlite.db")
    Rock.node(ctx,  port-8080)
end
function command(h, p, c)
    Rock.command(Rock.client("localhost", 8081), c)
end
function append(h,p,k,v)
    command(h,p, Command(Append, k, v))
end

function put(h,p,k,v)
    command(h,p, Command(Put, k, v))
end
end

import KVStore
peers = [("localhost", 8081), ("localhost", 8082), ("localhost", 8083) ]::Array{Tuple{ASCIIString, Int64},1}
for (i, worker) = enumerate(w)
    @spawnat worker (mktempdir(".") do f
                     KVStore.mk(8080+i, peers, f)
                     end)
end
wait(Timer(1))
KVStore.put("localhost",8081, "Test", "Put")
# write your own tests here
@test 1 == 1

max_stars_repo_name: JeremyRubin/Rock


text:
using MLDatasets
using MLDataUtils
using Plots
using Images
using Interact

include("vae2.jl")

function reconstruct(m::VAE, x::AbstractVector)
    x = reshape(x, length(x), 1)
    x_rec = decode(m, m(x))
    return x_rec
end

function show_pic(x)
    a = reshape(x, 28, 28)'
    return plot(Gray.(a))
end

function show_recon(m, X, device; n=5)
    subplots = []
    cpu = CPU()
    for i in rand(1:size(X, 2), n)
        x = X[:, i]
        x_ = reconstruct(m, device(x))
        p = show_pic(cpu(x))
        p_ = show_pic(cpu(x_))
        push!(subplots, p, p_)
    end
    plot(subplots..., layout=(n, 2))
end

function interpolate_latent_var(m, x, z_idx, device)
    vals = collect(-2:0.5:2)
    z = m(x)
    xs_ = []
    cpu = CPU()
    for v in vals
        z[z_idx, :] = v
        x_ = decode(m, device(z)) |> cpu
        push!(xs_, x_)
    end
    subplots = [show_pic(x_) for x_ in xs_]
    # plot(subplots..., layout=length(subplots))
    return subplots
end

function show_latent_vars(m, x, z_idxs, device)
    groups = [interpolate_latent_var(m, x, z_idx, device) for z_idx in z_idxs]
    subplots = vcat(groups...)
    n_cols = length(groups[1])
    plot(subplots..., layout=(length(z_idxs), n_cols))
end

function show_samples(m, n, device)
    z_len = length(m.enc2mu.b)  # assuming Linear
    z = randn(z_len, n)
    z = device(z)
    x_ = decode(m, z) |> CPU()
    subplots = [show_pic(x_[:, i]) for i in 1:size(x_, 2)]
    plot(subplots..., layout=length(subplots))
end

function main()
    device = best_available_device()
    m = VAE(
        Sequential(
            Linear(784 => 400),
            x -> relu.(x)),
        Linear(400 => 20),
        Linear(400 => 20),
        Sequential(
            Linear(20 => 400),
            x -> relu.(x),
            Linear(400 => 784),
            x -> logistic.(x));
        beta=5)
    m = m |> device

    X, _ = MNIST.traindata()
    X = convert(Matrix{Float64}, reshape(X, 784, 60000))
    @time m = fit!(m, X, device=device, opt=Adam(; lr=1e-3), n_epochs=10)

    show_recon(m, X, device, n=5)
    show_latent_vars(m, device(X[:, 2]), 1:4, device)
    show_samples(m, 10, device)
end


text:
"""
	gmtgravmag3d(cmd0::String=""; kwargs...)

Compute the gravity/magnetic anomaly of a 3-D body by the method of Okabe.

Full option list at [`gmtgravmag3d`]($(GMTdoc)gmtgravmag3d.html)

```julia
	G = gmtgravmag3d(M=(shape=:prism, params=(1,1,1,5)), I=1.0, R="-15/15/-15/15", H="10/60/10/-10/40");
	imshow(G)
```
"""
function gmtgravmag3d(cmd0::String="", arg1=nothing; kwargs...)

	length(kwargs) == 0 && return monolitic("gmtgravmag3d", cmd0, arg1)

	d = init_module(false, kwargs...)[1]		# Also checks if the user wants ONLY the HELP mode

	cmd = parse_common_opts(d, "", [:I :R :V_params :bi :f])[1]
	cmd = parse_these_opts(cmd, d, [[:C :density], [:E :thickness], [:F :track], [:G :grid :outgrid],
	                                [:L :observation_level], [:S :radius], [:Z :reference_level]])
	cmd = add_opt(d, cmd, 'H', [:H :mag_params], (field_dec="", field_dip="", mag="", mag_dec="", mag_dip=""))

	arg2 = nothing;
	if ((val = find_in_dict(d, [:Tv :index], true, "GMTdataset | Mx3 array | String")[1]) !== nothing && val != "")  opt = " -Tv"
	elseif ((val = find_in_dict(d, [:Tr :raw_triang], true, "GMTdataset | Mx3 array | String")[1]) !== nothing && val != "") opt = " -Tr"
	elseif ((val = find_in_dict(d, [:Ts :stl :STL], true, "String (file name)")[1]) !== nothing && val != "")  opt = " -Ts"
	elseif ((val = find_in_dict(d, [:M :body], false, "String | NamedTuple")[1]) !== nothing && val != "")
		cmd = add_opt(d, cmd, 'M', [:M :body], (shape="+s", params=","))
		opt = " -Ts"
	elseif (show_kwargs[1])  opt = ""
	else   error("Missing one of 'index', 'raw_triang' or 'str' data")
	end
	if (opt != " -Ts")		# The STL format can only be requested via file
		cmd *= opt
		if (isa(val, Array{<:Real}) || isa(val, GMTdataset) || isa(val, Vector{<:GMTdataset}))
			(arg1 === nothing) ? arg1 = val : arg2 = val		# Find the free slot
		else
			cmd *= arg2str(val)
		end
		if (opt == " -Tv")
			(!occursin(" -I", cmd)) && error("For grid output MUST specify grid increment ('I' or 'inc')")
		end
	end
	(!occursin(" -F", cmd) && !occursin(" -G", cmd)) && (cmd *= " -G")

	(opt != " -Tv") && return finish_PS_module(d, "gmtgravmag3d " * cmd, "", true, false, false, arg1, arg2)
	common_grd(d, cmd0, cmd, "gmtgravmag3d ", arg1, arg2)[2]		# Finish build cmd and run it
end


text:
"Square the sum of the first `n` positive integers"
function square_of_sum(n)
    (n*(n +1 ) ÷ 2)^2
end

"Sum the squares of the first `n` positive integers"
function sum_of_squares(n)
    (n*(n + 1)*(2n + 1)) ÷ 6

end

"Subtract the sum of squares from square of the sum of the first `n` positive ints"
function difference(n)
    square_of_sum(n)-sum_of_squares(n)
end

max_stars_repo_name: Rigante-pl/Exercism-Julia

max_issues_repo_name: Rigante-pl/Exercism-Julia

max_forks_repo_name: Rigante-pl/Exercism-Julia


text:
using Statistics: mean
using Base.Iterators: cycle
using IterTools: ncycle
using Knet: Knet, AutoGrad, Data, param, param0, mat, RNN, dropout, value, nll,
        adam, minibatch, progress!, progress, converge, cpucopy, logsoftmax, softmax

struct Embed; w; end

Embed(embed::Int, vocab::Int)=Embed(param(embed,vocab))

(e::Embed)(x) = e.w[:,x]

struct Linear; w; b; end

Linear(input::Int, output::Int)=Linear(param(output,input), param0(output))

(l::Linear)(x) = l.w * mat(x,dims=1) .+ l.b  # (H,B,T)->(H,B*T)->(V,B*T)

struct LM
        embed::Embed
        rnn::RNN
        linear::Linear
end

function (m::LM)(x)
        x = m.embed(x)
        x = m.rnn(x)
        x = m.linear(x)
        x
end

(m::LM)(x,y) = nll(m(x),y)
(m::LM)(d::Data) = mean(m(x,y) for (x,y) in d)

# The h=0,c=0 options to RNN enable a persistent state between iterations
LM(embedsz::Int, vocabsz::Int, hiddensz::Int, outputsz::Int; o...) =
    LM(Embed(embedsz, vocabsz), RNN(embedsz,hiddensz;h=0,c=0,o...), Linear(hiddensz, outputsz))

LM(name::String) = Knet.load(KongYiji.dir("lm.jld2"), name)

function prepare(xs::Vector{<:Integer}, batchsz, seqlen)
        nbatch = length(xs) ÷ batchsz
        x = reshape(xs[1:nbatch * batchsz], nbatch, batchsz)' # reshape full data to (B,N) with contiguous rows
        minibatch(x[:, 1:nbatch - 1], x[:, 2:nbatch], seqlen) # split into (B,T) blocks
end

function train!(m::LM, seq::Vector{<:Integer}, epochs, batchsz, seqlen)
        data = prepare(seq, batchsz, seqlen)
        losses = [begin
                        @show "GC"
                        m.rnn.h = 0
                        m.rnn.c = 0
                        Knet.gc()
                        loss
                  end
                  for loss in every(progress(adam(m, ncycle(data, epochs))), length(data))
        ]
        #Knet.gc()
        return losses
end

every(itr, n) = (x for (i,x) in enumerate(itr) if i % n == 0);

hiddens(m::LM) = (m.rnn.h, m.rnn.c)

function reset!(m::LM, hc)
        m.rnn.h = hc[1]
        m.rnn.c = hc[2]
end

function LM(corpus; epochs=100, batchsz=256, seqlen=100, o...)
        seqs = posidsents(corpus)
        np = length(postags(corpus))
        embedsz::Int = np * 2
        vocabsz::Int = np
        hiddensz::Int = np #todo
        outputsz::Int = np

        nn = LM(embedsz, vocabsz, hiddensz, outputsz; o...)
        seq = collect(Iterators.flatten(seqs))
        train!(nn, seq, epochs, batchsz, seqlen)
        nn = nn |> cpucopy
        reset!(nn, (0, 0))
        return nn
end

function generate!(lm::LM, start::Int; words=nothing, len=100)

        r = String[]
        getword(i) = isnothing(words) ? string(i) : words[i]
        x = start
        for _ in 1:len
                push!(r, getword(x))
                py = softmax(lm(x))
                rd, acc = rand(), 0.
                for (i, p) in enumerate(py)
                        acc += p
                        if rd <= acc; x = i; break; end
                end
        end
        r
end


text:
"""
    newdims = Dims_create(nnodes::Integer, dims)

A convenience function for selecting a balanced Cartesian grid of a total of `nnodes`
nodes, for example to use with [`MPI.Cart_create`](@ref).

`dims` is an array or tuple of integers specifying the number of nodes in each dimension.
The function returns an array `newdims` of the same length, such that if `newdims[i] =
dims[i]` if `dims[i]` is non-zero, and `prod(newdims) == nnodes`, and values `newdims` are
as close to each other as possible.

`nnodes` should be divisible by the product of the non-zero entries of `dims`.

# External links
$(_doc_external("MPI_Dims_create"))
"""
function Dims_create(nnodes::Integer, dims)
    dims = Cint[dim for dim in dims]
    @mpichk ccall((:MPI_Dims_create, libmpi), Cint,
                  (Cint, Cint, Ptr{Cint}), nnodes, length(dims), dims)
    return dims
end

"""
    comm_cart = Cart_create(comm::Comm, dims; periodic=map(_->false, dims), reorder=false)

Create new MPI communicator with Cartesian topology information attached.

`dims` is an array or tuple of integers specifying the number of MPI processes in each
coordinate direction, and `periodic` is an array or tuple of `Bool`s indicating the
periodicity of each coordinate. `prod(dims)` must be less than or equal to the size of
`comm`; if it is smaller than some processes are returned a null communicator.

If `reorder == false` then the rank of each process in the new group is identical to its
rank in the old group, otherwise the function may reorder the processes.

See also [`MPI.Dims_create`](@ref).

# External links
$(_doc_external("MPI_Cart_create"))
"""
function Cart_create(comm::Comm, dims; periodic = map(_->false, dims), reorder=false)
    if !(dims isa Array{Cint})
        dims = Cint[dim for dim in dims]
    end
    if !(periodic isa Array{Cint})
        periodic = Cint[flag for flag in periodic]
    end
    length(dims) == length(periodic) || error("dims and periodic arguments are not the same length")
    comm_cart = Comm()
    # int MPI_Cart_create(MPI_Comm comm_old, int ndims, const int dims[],
    #                     const int periods[], int reorder, MPI_Comm *comm_cart)
    @mpichk ccall((:MPI_Cart_create, libmpi), Cint,
                  (MPI_Comm, Cint, Ptr{Cint}, Ptr{Cint}, Cint, Ptr{MPI_Comm}),
                  comm, length(dims), dims, periodic, reorder, comm_cart)
    if comm_cart != COMM_NULL
        finalizer(free, comm_cart)
    end
    comm_cart
end

"""
    rank = Cart_rank(comm::Comm, coords)

Determine process rank in communicator `comm` with Cartesian structure.  The `coords`
array specifies the 0-based Cartesian coordinates of the process. This is the inverse of [`MPI.Cart_coords`](@ref)

# External links
$(_doc_external("MPI_Cart_rank"))
"""
function Cart_rank(comm::Comm, coords)
    if !(coords isa Vector{Cint})
        coords = Cint[coord for coord in coords]
    end
    rank = Ref{Cint}()
    # int MPI_Cart_rank(MPI_Comm comm, const int coords[], int *rank)
    @mpichk ccall((:MPI_Cart_rank, libmpi), Cint,
                  (MPI_Comm, Ptr{Cint}, Ptr{Cint}),
                  comm, coords, rank)
    Int(rank[])
end

"""
    dims, periods, coords = Cart_get(comm::Comm)

Obtain information on the Cartesian topology of dimension `N` underlying the
communicator `comm`. This is specified by two `Cint` arrays of `N` elements
for the number of processes and periodicity properties along each Cartesian dimension.
A third `Cint` array is returned, containing the Cartesian coordinates of the calling process.

# External links
$(_doc_external("MPI_Cart_get"))
"""
function Cart_get(comm::Comm)
    maxdims = Cartdim_get(comm)
    # preallocate with nontrivial values
    dims    = Cint[-1 for i = 1:maxdims]
    periods = Cint[-1 for i = 1:maxdims]
    coords  = Cint[-1 for i = 1:maxdims]
    # int MPI_Cart_get(MPI_Comm comm, int maxdims, int dims[], int periods[], int coords[])
    @mpichk ccall((:MPI_Cart_get, libmpi), Cint,
                  (MPI_Comm, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}),
                  comm, maxdims, dims, periods, coords)
    return dims, periods, coords
end

"""
    ndims = Cartdim_get(comm::Comm)

Return number of dimensions of the Cartesian topology associated with the communicator `comm`.

# External links
$(_doc_external("MPI_Cartdim_get"))
"""
function Cartdim_get(comm::Comm)
    dims = Ref{Cint}()
    # int MPI_Cartdim_get(MPI_Comm comm, int *ndims)
    @mpichk ccall((:MPI_Cartdim_get, libmpi), Cint,
                  (MPI_Comm, Ptr{Cint}),
                  comm, dims)
    return Int(dims[])
end

"""
    coords = Cart_coords(comm::Comm, rank::Integer=Comm_rank(comm))

Determine coordinates of a process with rank `rank` in the Cartesian communicator
`comm`. If no `rank` is provided, it returns the coordinates of the current process.

Returns an integer array of the 0-based coordinates. The inverse of [`Cart_rank`](@ref).

# External links
$(_doc_external("MPI_Cart_coords"))
"""
function Cart_coords(comm::Comm, rank::Integer=Comm_rank(comm))
    maxdims = Cartdim_get(comm)
    coords = Vector{Cint}(undef, maxdims)
    # int MPI_Cart_coords(MPI_Comm comm, int rank, int maxdims, int coords[])
    @mpichk ccall((:MPI_Cart_coords, libmpi), Cint,
                  (MPI_Comm, Cint, Cint, Ptr{Cint}),
                  comm, rank, maxdims, coords)
    return Int.(coords)
end

"""
    rank_source, rank_dest = Cart_shift(comm::Comm, direction::Integer, disp::Integer)

Return the source and destination ranks associated to a shift along a given
direction.

# External links
$(_doc_external("MPI_Cart_shift"))
"""
function Cart_shift(comm::Comm, direction::Integer, disp::Integer)
    rank_source = Ref{Cint}()
    rank_dest   = Ref{Cint}()
    # int MPI_Cart_shift(MPI_Comm comm, int direction, int disp,
    #                    int *rank_source, int *rank_dest)
    @mpichk ccall((:MPI_Cart_shift, libmpi), Cint,
                  (MPI_Comm, Cint, Cint, Ptr{Cint}, Ptr{Cint}),
                  comm, direction, disp, rank_source, rank_dest)
    Int(rank_source[]), Int(rank_dest[])
end

"""
    comm_sub = Cart_sub(comm::Comm, remain_dims)

Create lower-dimensional Cartesian communicator from existent Cartesian
topology.

`remain_dims` should be a boolean vector specifying the dimensions that should
be kept in the generated subgrid.

# External links
$(_doc_external("MPI_Cart_sub"))
"""
function Cart_sub(comm::Comm, remain_dims)
    if !(remain_dims isa Array{Cint})
        remain_dims = Cint[dim for dim in remain_dims]
    end
    comm_sub = Comm()
    # int MPI_Cart_sub(MPI_Comm comm, const int remain_dims[], MPI_Comm *comm_new)
    @mpichk ccall((:MPI_Cart_sub, libmpi), Cint,
                  (MPI_Comm, Ptr{Cint}, Ptr{MPI_Comm}),
                  comm, remain_dims, comm_sub)
    if comm_sub != COMM_NULL
        finalizer(free, comm_sub)
    end
    comm_sub
end


text:
struct Monomial{N}
    θ::SVector{N,Int}
end
Monomial(p::Int) = Monomial(SVector(p))
Monomial(p::Tuple) = Monomial(SVector(p))

dimension(m::Monomial{N}) where {N} = N

Base.exponent(m::Monomial) = m.θ

function (m::Monomial)(x::SVector)
    @assert length(x) == dimension(m)
    x .^ exponent(m) |> prod
end
(m::Monomial)(x)  = m(SVector(x))

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
### derivatives-type.jl
#
# Copyright (C) 2016, 2017 Mosè Giordano.
#
# Maintainer: Mosè Giordano <mose AT gnu DOT org>
# Keywords: uncertainty, error propagation, physics
#
# This file is a part of Measurements.jl.
#
# License is MIT "Expat".
#
### Commentary:
#
# This file containes the definition of Derivatives type.  This is borrowed from
# Base.ImmutableDict type (see base/dict.jl in Julia source code), introduced in
# Julia 0.5.  `get' and `getindex' have methods different from those of the
# standard Base.ImmutabileDict, otherwise Derivatives and Base.ImmutableDict
# would be identical.
#
### Code:

struct Derivatives{T} <: AbstractDict{Tuple{T,T,UInt64},T}
    parent::Derivatives{T}
    key::Tuple{T,T,UInt64}
    value::T
    Derivatives{T}() where {T} = new() # represents an empty dictionary
    Derivatives{T}(parent::Derivatives, key, value) where {T} = new(parent, key, value)
end

Derivatives(t::Derivatives{T}, KV::Pair) where {T} = Derivatives{T}(t, KV[1], KV[2])

Base.convert(::Type{Derivatives{T}}, ders::Derivatives{T}) where {T} = ders
# XXX: this function is very slow, but shouldn't be called often.
function Base.convert(::Type{Derivatives{T}}, ders::Derivatives{S}) where {T,S}
    new = Derivatives{T}()
    for der in ders
        new = Derivatives{T}(new, Tuple{T,T,T}(der[1]), T(der[2]))
    end
    return new
end

function Base.getindex(dict::Derivatives, key)
    while isdefined(dict, :parent)
        # We want to compare the three elements of the key tuples.  In order to allow for
        # calculations involving NaNs we need to use `isequal` for comparing the first two
        # elements, then we can use `==` to compare the tag.  Since `==` is faster than
        # `isequal`, we first compare the last element and then the other two.  This is the
        # only difference with the plain Base.ImmutableDict type.  The same in `get'.
        dict.key[3] == key[3] && isequal(dict.key[1], key[1]) &&
            isequal(dict.key[2], key[2]) && return dict.value
        dict = dict.parent
    end
    throw(KeyError(key))
end
function Base.get(dict::Derivatives, key, default)
    while isdefined(dict, :parent)
        # See comment in `getindex'.
        dict.key[3] == key[3] && isequal(dict.key[1], key[1]) &&
            isequal(dict.key[2], key[2]) && return dict.value
        dict = dict.parent
    end
    return default
end

# this actually defines reverse iteration (e.g. it should not be used for merge/copy/filter type operations)
function Base.iterate(derivs::Derivatives{T}, t=derivs) where T
    isdefined(t, :parent) || return nothing
    (Pair{Tuple{T,T,UInt64},T}(t.key, t.value), t.parent)
end

Base.length(t::Derivatives) = count(x->true, t)


text:
function get_humanoid(;
    timestep=0.01,
    gravity=[0.0; 0.0; -9.81],
    friction_coefficient=0.8,
    spring=0.0,
    damper=0.0,
	contact_feet=true,
    contact_body=false,
    T=Float64)

    path = joinpath(@__DIR__, "../deps/humanoid.urdf")
    mech = Mechanism(path, true, T,
        gravity=gravity,
        timestep=timestep,
        spring=spring,
        damper=damper)

    if contact_feet
        origin = Origin{T}()
        bodies = mech.bodies
        eqs = mech.joints

        # Foot contact
        left_foot = get_body(mech, :left_foot)

		aa = -0.43000 * [-0.44721, 0.00000, 0.89442]
		ql = axis_angle_to_quaternion(aa)
        qll = ql * RotX(-1.57080)*RotY(1.47585)*RotZ(-1.47585) # Quaternion(RotXYZ(roll=-1.57080, pitch=1.47585, yaw=-1.47585)) # roll pitch yaw
        qlr = ql * RotX(+1.57080)*RotY(1.47585)*RotZ(+1.47585) # Quaternion(RotXYZ(roll=+1.57080, pitch=1.47585, yaw=+1.47585)) # roll pitch yaw

        pfll = vector_rotate([ 0.5 * left_foot.shape.shape[1].rh[2] + 0.03500; -0.03; 0.0], qll)
        pbll = vector_rotate([-0.5 * left_foot.shape.shape[1].rh[2] + 0.03500; -0.03; 0.0], qll)
        pflr = vector_rotate([ 0.5 * left_foot.shape.shape[1].rh[2] + 0.03500; +0.01; 0.0], qlr)
		pblr = vector_rotate([-0.5 * left_foot.shape.shape[1].rh[2] + 0.03500; +0.01; 0.0], qlr)
        p = [0.0,0.054,0.]
        o = left_foot.shape.shape[1].rh[1]
        contacts = [
					p,
                    # pfll,
                    # pbll,
                    # pflr,
                    # pblr,
                   ]
        contact_radius = [
                    o,
                    # o,
                    # o,
                    # o,
                  ]
        n = length(contacts)
        normal = [[0.0; 0.0; 1.0] for i = 1:n]
        friction_coefficients = friction_coefficient * ones(T, n)

        contacts_left = contact_constraint(left_foot, normal,
            friction_coefficient=friction_coefficients,
            contact_origins=contacts,
            contact_radius=contact_radius)

        right_foot = get_body(mech, :right_foot)

        pfr = [0.5 * right_foot.shape.shape[1].rh[2]; 0.0; 0.0]
        ofr = right_foot.shape.shape[1].rh[1]
        pbr = [-0.5 * right_foot.shape.shape[1].rh[2]; 0.0; 0.0]
        obr = right_foot.shape.shape[1].rh[1]

        contacts = [
                    pfr,
                    pbr,
                   ]

        contact_radius = [
                            ofr,
                            obr,
                         ]

        n = length(contacts)
        normal = [[0.0; 0.0; 1.0] for i = 1:n]
        friction_coefficients = friction_coefficient * ones(T, n)

        contacts_right = contact_constraint(right_foot, normal,
            friction_coefficient=friction_coefficients,
            contact_origins=contacts,
            contact_radius=contact_radius)

        set_minimal_coordinates!(mech, get_joint(mech, :floating_base), [0.0; 0.0; 1.2; 0.1; 0.0; 0.0])
        # mech = Mechanism(origin, bodies, eqs, [contacts_left; contacts_right],
            # gravity=gravity,
            # timestep=timestep,
            # spring=spring,
            # damper=damper)
        mech = Mechanism(origin, bodies, eqs, [contacts_left; ],
            gravity=gravity,
            timestep=timestep,
            spring=spring,
            damper=damper)
    end

    return mech
end

function initialize_humanoid!(mechanism::Mechanism{T};
    body_position=[0.0, 0.0, 1.5],
    body_orientation=[0.1, 0.0, 0.0]) where T
    set_minimal_coordinates!(mechanism,
        get_joint(mechanism, :floating_base),
        [body_position; body_orientation])
    zero_velocity!(mechanism)
end


text:
using CompoundPeriods, Dates
using Test

cperiod = Day(5) + Hour(17) + Minute(35)
rperiod = reverse(cperiod)

@test rperiod.cperiod == cperiod

cperiod = canonical(Minute(-3600) + Second(22))
@test (cperiod - Day(cperiod)) > Hour(12)

@test TimeUnits(Minute(10)+Microsecond(100)) == Microsecond(600000100)

@test Month(Year(3)+Month(2)) == Month(2)
@test hour(Day(5)+Hour(3)+Second(1)) == 3

@test isempty(Hour(0) + Minute(0))
@test canonical(Hour(0)+Minute(0)) == Nanosecond(0)

@test Year(2)+Month(12)+Hour(7) > Year(1)+Month(24)+Minute(5)
@test Year(2)+Month(12)+Hour(7) == Year(1)+Month(24)+Hour(7)


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

name = "MariaDB_Connector_ODBC"
version = v"3.1.7"

# Collection of sources required to build MariaDB_Connector_ODBC
sources = [
    ArchiveSource("https://downloads.mariadb.com/Connectors/odbc/connector-odbc-3.1.7/mariadb-connector-odbc-3.1.7-ga-src.tar.gz",
                  "699c575e169d770ccfae1c1e776aa7725d849046476bf6579d292c89e8c8593e"),
    FileSource("https://downloads.mariadb.com/Connectors/odbc/connector-odbc-3.1.7/mariadb-connector-odbc-3.1.7-win64.msi",
               "ef1ad796cc4aba67aa6606b35239302c579b7894604239d1caaa1f2d10623bd2"; filename = "x86_64-w64-mingw32.msi"),
    FileSource("https://downloads.mariadb.com/Connectors/odbc/connector-odbc-3.1.7/mariadb-connector-odbc-3.1.7-win32.msi",
               "a60b74dc3af0b450d892731c2973037beef6fc972e0813f3db9f28f083a63402"; filename = "i686-w64-mingw32.msi"),
    ## Keep the patches just in case some day we decide to build for Windows
    ## from source
    # DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
install_license $WORKSPACE/srcdir/mariadb-connector*/COPYING

if [[ "${target}" == *-mingw* ]]; then
    # For Windows just use the prebuilt library
    mkdir -p "${libdir}"
    cd $WORKSPACE/srcdir
    apk add p7zip
    7z x "${target}.msi"
    cp mariadb_odbc_dll "${libdir}/mariadb_odbc.dll"
    exit
fi

cd $WORKSPACE/srcdir/mariadb-connector*/

# Skip building of macOS package
sed -i 's/ADD_SUBDIRECTORY(osxinstall)/# ADD_SUBDIRECTORY(osxinstall)/' CMakeLists.txt

## Keep this for reference in case we decide to build for Windows from source
# if [[ "${target}" == *-mingw* ]]; then
#     # Disable dialog plugin
#     sed -i 's/-DPLUGINS_LIB_DIR.*:dialog.*//' CMakeLists.txt
#     atomic_patch -p1 "${WORKSPACE}/srcdir/patches/0001-fix-case-headers-windows.patch"
#     # The following is not really a patch: Microsoft Visual C++ thinks that it's
#     # a good idea to put some stray characters in the first bits of the file,
#     # making it hard to digest for MinGW.  We can't diff the file as the
#     # original one appears to be binary, because of the stray character.
#     cp "${WORKSPACE}/srcdir/patches/resource.h" dsn/
# fi

mkdir build && cd build

export CFLAGS="-I${includedir}/mariadb"
# Find the MariaDB lib and link to the right OpenSSL and Zlib libraries
export LDFLAGS="-L${libdir}/mariadb ${libdir}/libssl.${dlext} ${libdir}/libcrypto.${dlext} -lz"
cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    -DWITH_OPENSSL=ON \
    -DWITH_SSL=OPENSSL \
    -DOPENSSL_FOUND=ON \
    -DOPENSSL_CRYPTO_LIBRARY=${libdir}/libcrypto.${dlext} \
    -DOPENSSL_SSL_LIBRARY=${libdir}/libssl.${dlext} \
    -DICONV_LIBRARIES=${libdir}/libiconv.${dlext} \
    -DICONV_INCLUDE_DIR=${includedir} \
    ..
make -j${nproc}
make install
"""

platforms = [
    Platform("i686", "linux"; libc="glibc"),
    Platform("x86_64", "linux"; libc="glibc"),
    # Platform("aarch64", "linux"; libc="glibc"),
    # Platform("armv7l", "linux"; libc="glibc"),
    # Platform("powerpc64le", "linux"; libc="glibc"),
    Platform("i686", "linux"; libc="musl"),
    Platform("x86_64", "linux"; libc="musl"),
    # Platform("aarch64", "linux"; libc="musl"),
    # Platform("armv7l", "linux"; libc="musl"),
    Platform("x86_64", "macos"),
    # Platform("x86_64", "freebsd"),
    Platform("i686", "windows"),
    Platform("x86_64", "windows"),
]

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

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency("MariaDB_Connector_C_jll"),
    Dependency("iODBC_jll"),
    Dependency("Libiconv_jll"),
    Dependency("unixODBC_jll"),
    Dependency("OpenSSL_jll"),
    Dependency("Zlib_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
function dataviewer(t; throttle = 0.1, nbins = 30, nbins_range = 1:100)
    (t isa AbstractObservable) || (t = Observable{Any}(t))

    columns_and_names = @map create_columns_from_iterabletable(getiterator(&t))

    names = @map (&columns_and_names)[2]

    dict = @map Dict((key, convert_missing.(val)) for (val, key)  in zip((&columns_and_names)...))
    x =  @nodeps dropdown(names, placeholder = "First axis", multiple = true)
    y =  @nodeps dropdown(names, placeholder = "Second axis", multiple = true)
    y_toggle = @nodeps togglecontent(y, value = false, label = "Second axis")
    plot_type = @nodeps dropdown(
        OrderedDict(
            "line"         => Plots.plot,
            "scatter"      => Plots.scatter,
            "bar"          => StatPlots.groupedbar,
            "boxplot"      => StatPlots.boxplot,
            "corrplot"     => StatPlots.corrplot,
            "cornerplot"   => StatPlots.cornerplot,
            "density"      => StatPlots.density,
            "histogram"    => StatPlots.histogram,
            "marginalhist" => StatPlots.marginalhist,
            "violin"       => StatPlots.violin
        ),
        placeholder = "Plot type")

    # Add bins if the plot allows it
    display_nbins = @map (&plot_type) in [corrplot, cornerplot, histogram, marginalhist] ? "block" : "none"
    nbins =  (@nodeps slider(nbins_range, extra_obs = ["display" => display_nbins], value = nbins, label = "number of bins"))
    nbins.scope.dom = Widgets.div(nbins.scope.dom, attributes = Dict("data-bind" => "style: {display: display}"))
    nbins_throttle = Observables.throttle(throttle, nbins)

    by = @nodeps dropdown(names, multiple = true, placeholder="Group by")
    by_toggle = @nodeps togglecontent(by, value = false, label = "Split data")
    plt = @nodeps button("plot")
    output = @map begin
        &plt
        if (plt[] == 0)
            plot()
        else
            x_cols = hcat(getindex.((&dict,), x[])...)
            has_y = y_toggle[] && !isempty(y[])
            has_by = by_toggle[] && !isempty(by[])
            y_cols = has_y ? [hcat(getindex.((&dict,), y[])...)] : []
            by_tup = Tuple(getindex(&dict, b) for b in by[])
            by_kwarg = has_by ? [(:group, NamedTuple{Tuple(by[])}(by_tup))] : []
            label = length(x[]) > 1 ? [(:label, x[])] :
                    (y_toggle[] && length(y[]) > 1) ? [(:label, y[])] : []
            densityplot1D = plot_type[] in [density, histogram]
            xlabel = (length(x[]) == 1 && (densityplot1D || has_y)) ? [(:xlabel, x[][1])] : []
            ylabel = (has_y && length(y[]) == 1) ? [(:ylabel, y[][1])] :
                     (!has_y && !densityplot1D && length(x[]) == 1) ? [(:ylabel, x[][1])] : []
            plot_type[](x_cols, y_cols...; nbins = &nbins_throttle, by_kwarg..., label..., xlabel..., ylabel...)
        end
    end
    wdg = Widget{:dataviewer}(["x" => x, "y" => y, "y_toggle" => y_toggle, "by" => by, "by_toggle" => by_toggle,
        "plot_type" => plot_type, "plot_button" => plt, "nbins" => nbins], output = output)
    @layout! wdg Widgets.div(
        Widgets.div(
            :x,
            :y_toggle,
            :plot_type,
            :by_toggle,
            :plot_button
        ),
        Widgets.div(style = Dict("width" => "3em")),
        Widgets.div(
            Widgets.observe(_),
            :nbins
        ),
        style = Dict("display" => "flex", "direction" => "row")
    )
end


text:
# This file was generated, do not modify it. # hide
ridge_model = RidgeRegressor(lambda=0.1)
ridge = machine(ridge_model, W, z)
ẑ = predict(ridge, W)


text:
"""
    PlanetaryConstants

Module containing physical constants and parameters characterizing the planet.
"""
module PlanetaryConstants

using ..UniversalConstants
using DocStringExtensions

export PlanetConstants,
       AbstractPlanetConstants,
       EarthConstants
       # MarsConstants

"""
    PlanetConstants{FT<:AbstractFloat}

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct PlanetConstants{FT<:AbstractFloat}
  "Von Karman constant (1)"
  k_Karman::FT
  "Molecular weight dry air (kg/mol)"
  molmass_dryair::FT
  "Gas constant dry air (J/kg/K)"
  R_d::FT = FT(gas_constant(FT)/molmass_dryair)
  "Adiabatic exponent dry air"
  kappa_d::FT
  "Isobaric specific heat dry air"
  cp_d::FT = FT(R_d/kappa_d)
  "Isochoric specific heat dry air"
  cv_d::FT = FT(cp_d - R_d)
  "Density of water (kg/m^3)"
  ρ_liq::FT
  "Density of liquid water (kg/m^3)"
  ρ_cloud_liq::FT = FT(ρ_liq)
  "Density of ice (kg/m^3)"
  ρ_ice::FT
  "Density of water ice (kg/m^3)"
  ρ_cloud_ice::FT = FT(ρ_ice)
  "Molecular weight (kg/mol)"
  molmass_water::FT
  "Molar mass ratio dry air/water"
  molmass_ratio::FT = FT(molmass_dryair/molmass_water)
  "Gas constant water vapor (J/kg/K)"
  R_v::FT = FT(gas_constant(FT)/molmass_water)
  "Isobaric specific heat vapor (J/kg/K)"
  cp_v::FT
  "Isobaric specific heat liquid (J/kg/K)"
  cp_l::FT
  "Isobaric specific heat ice (J/kg/K)"
  cp_i::FT
  "Isochoric specific heat vapor (J/kg/K)"
  cv_v::FT = FT(cp_v - R_v)
  "Isochoric specific heat liquid (J/kg/K)"
  cv_l::FT = FT(cp_l)
  "Isochoric specific heat ice (J/kg/K)"
  cv_i::FT = FT(cp_i)
  "Freezing point temperature (K)"
  T_freeze::FT
  "Minimum temperature guess in saturation adjustment (K)"
  T_min::FT
  "Maximum temperature guess in saturation adjustment (K)"
  T_max::FT
  "Homogeneous nucleation temperature (K)"
  T_icenuc::FT
  "Triple point temperature (K)"
  T_triple::FT
  "Reference temperature (K)"
  T_0::FT = FT(T_triple)
  "Latent heat vaporization at T_0 (J/kg)"
  LH_v0::FT
  "Latent heat sublimation at T_0 (J/kg)"
  LH_s0::FT
  "Latent heat of fusion at T_0 (J/kg)"
  LH_f0::FT = FT(LH_s0 - LH_v0)
  "Specific internal energy of vapor at T_0 (J/kg)"
  e_int_v0::FT = FT(LH_v0 - R_v*T_0)
  "Specific internal energy of ice at T_0 (J/kg)"
  e_int_i0::FT = FT(LH_f0)
  "Triple point vapor pressure (Pa)"
  press_triple::FT
  "Reference density sea water (kg/m^3)"
  ρ_ocean::FT
  "Specific heat sea water (J/kg/K)"
  cp_ocean::FT
  "Mean planetary radius (m)"
  planet_radius::FT
  "Length of day (s)"
  day::FT
  "Ang. velocity planetary rotation (1/s)"
  Omega::FT
  "Gravitational acceleration (m/s^2)"
  grav::FT
  "Length of anomalistic year (s)"
  year_anom::FT = FT(365.26*day)
  "Length of semimajor orbital axis (m)"
  orbit_semimaj::FT
  "Total solar irradiance (W/m^2)"
  TSI::FT
  "Mean sea level pressure (Pa)"
  MSLP::FT
end

abstract type AbstractPlanetConstants end

struct EarthConstants <: AbstractPlanetConstants end
include("EarthConstants.jl")

# struct MarsConstants <: AbstractPlanetConstants end
# include("MarsConstants.jl")

end

max_stars_repo_name: yairchn/CLIMA

max_issues_repo_name: yairchn/CLIMA

max_forks_repo_name: yairchn/CLIMA


text:
#---------------------------------------------------------------------#
#This function computes the Initial & Exact Solutions.
#Written by F.X. Giraldo on April 19, 2019
#           Department of Applied Mathematics
#           Naval Postgraduate School
#           Monterey; CA 93943-5216
#---------------------------------------------------------------------#
function exact_solution(coord,Npoin,time,case,DFloat)

#Set some constants
xmin=0
xmax=+2
xc=0.5*(xmax+xmin)
xl=xmax-xmin
c=2.0

#Initialize
qe=zeros(DFloat,Npoin,2)
qb=zeros(DFloat,Npoin,1)
gravity=DFloat(1.0)
hmean=DFloat(1.0)

#Construct Solution
for I=1:Npoin
    x=coord[I]
    h=0.5*cos(c*pi*x)*cos(c*π*time)
    u=0.5*sin(c*pi*x)*sin(c*π*time)
    hb=hmean
    qe[I,1]=h
    qe[I,2]=u
    qb[I]=hb
end #I
return (qe,qb,gravity)
end


text:
"""
A Continuable is a function which takes a single argument function a -> b

We could have made most of the core functionality work for general functions a... -> b
however the semantics become too complicated very soon.

For instance for product. How should the product look like? say we have two continuables and
name the continuation's arguments a... and  b... respectively.
How to call the follow up continuation then?
cont(a..., b...) or cont(a, b)?
Probably something like if isa(a and b, SingletonTuple) then cont(a...,b...)
else cont(a,b). However this seems to complicate things unnecessarily complex.

Another point is for example the interaction with iterables which always deliver tuples.

See also this github issue https://github.com/JuliaLang/julia/issues/6614 for tuple deconstructing which
would give a similar familarity.

Many documentation strings are taken and adapted from
https://github.com/JuliaCollections/Iterators.jl/blob/master/src/Iterators.jl.
"""
# TODO copy documentation strings
module Continuables
export
  cont, @cont, AbstractContinuable, Continuable, innerfunctype,
  @Ref, stoppable, stop,
  emptycontinuable, singleton, repeated, iterated,
  aschannel, ascontinuable, i2c, @i2c,
  reduce, reduce!, zip, product, chain, flatten, cycle, foreach, map, all, any, sum, prod,
  take, takewhile, drop, dropwhile, partition, groupbyreduce, groupby,
  nth

using Compat
using ExprParsers
using DataTypesBasic
using OrderedCollections
import Base.Iterators: cycle, flatten, take, drop, partition, product

include("utils.jl")
include("itertools.jl")

## Continuable Core -------------------------------------------------------------------------

"""
`cont` is reserved function parameter name
"""
cont(args...; kwargs...) = error("`cont` is reserved to be used within `Continuables.@cont`")

"""
    AbstractContinuable

Abstract type which all continuable helper functions use for dispatch.

The interface for a continuable just consists of
```
Base.foreach(cont, continuable::YourContinuable)
```

For julia 1.0 you further need to provide
```
(continuable::YourContinuable)(cont) = foreach(cont, continuable)
```
which is provided automatically for more recent julia versions.
"""
abstract type AbstractContinuable end

"""
    foreach(func, continuable)

Runs the continuable with `func` as the continuation. This is the core interface of a continuable.

It is especially handy when using `do` syntax.
"""
function Base.foreach(cont, ::AbstractContinuable)
  error("You need to specialize `Base.foreach` for your custom Continuable type.")
end

"""
  Continuable(func)

Assumes func to have a single argument, the continuation function (usually named `cont`).

Example
-------

```
Continuable(function (cont)
  for i in 1:10
    cont(i)
  end
end)
```
"""
struct Continuable{Func} <: AbstractContinuable
  f::Func
  # making sure typevar `Func` has always the correct meaning so that we can dispatch on it
  # otherwise one could construct `Continuable{Any}(func)` which would break dispatching on the typevariable.
  Continuable(f) = new{typeof(f)}(f)
end
Base.foreach(cont, c::Continuable) = c.f(cont)
(c::Continuable)(cont) = foreach(cont, c)

"""
    innerfunctype(continuable)

Returns the generic type of `continuable.f`.
This can be used to specialize functions on specific Continuables, like `Base.IteratorSize` and else.

# Examples

```julia
mycont(x) = @cont cont(x)
Base.length(::Continuable{<:innerfunctype(mycont(:examplesignature))}) = :justatest
length(mycont(1)) == :justatest
length(mycont("a")) == :justatest
# other continuables are not affected
anothercont(x) = @cont cont(x)
length(anothercont(42)) == 1
```
"""
function innerfunctype(continuable::Continuable)
  typeof(continuable.f).name.wrapper
end

include("./syntax.jl")  # parts of syntax.jl require the Continuable definition, hence here

## conversions ----------------------------------------------

"""
    aschannel(continuable) -> Channel

Convert the continuable into a channel. Performance is identical compared to when you would have build a Channel
directly.
"""
function aschannel(continuable::AbstractContinuable, size=0; elemtype=Any, taskref=nothing, spawn=false)
  Channel{elemtype}(size, taskref=taskref, spawn=spawn) do channel
    continuable() do x
      put!(channel, x)
    end
  end
end

"""
    ascontinuable(iterable)

Converts an iterable to a Continuable. There should not be any performance loss.
"""
ascontinuable(iterable) = @cont foreach(cont, iterable)
"""
    i2c(iterable)

Alias for [`ascontinuable`](@ref). "i2c" is meant as abbreviation for "iterable to continuable".
Also comes in macro form as [`@i2c`](@ref).
"""
const i2c = ascontinuable
"""
    @i2c iterable

Alias for [`ascontinuable`](@ref). "i2c" is meant as abbreviation for "iterable to continuable".
Also comes in function form as [`i2c`](@ref).
"""
macro i2c(expr)
  esc(:(ascontinuable($expr)))
end

## factories -----------------------------------------------

"""
    emptycontinuable

The continuable with no elements.
"""
const emptycontinuable = @cont ()  # mind the space!!

"""
    singleton(value)

Construct a Continuable containing only the one given value.
"""
singleton(value) = @cont cont(value)

"""
    repeated(() -> 2[, n])
    repeated([n]) do
      # ...
      "returnvalue"
    end

Constructs a Continuable which repeatedly yields the returnvalue from calling the given function again and again.
Until infinity or if `n` is given, exactly `n` times.
"""
@cont function repeated(f)
  while true
    cont(f())
  end
end

@cont function repeated(f, n::Integer)
  for _ in 1:n
    cont(f())
  end
end

"""
    iterated((x) -> x*x, startvalue)
    iterated(startvalue) do x
      # ...
      x*x
    end

Constructs an infinite Continuable which
"""
@cont @Ref function iterated(f, x)
  a = Ref(x)
  cont(a)
  while true
    a = f(a)
    cont(a)
  end
end

## core helpers ----------------------------------------------------------

# Continuable cannot implement Base.iterate efficiently

try
  # only in more recent julia, we can specialize call syntax for AbstractTypes
  (c::AbstractContinuable)(cont) = foreach(cont, c)
catch
end

Base.IteratorSize(::AbstractContinuable) = Base.SizeUnknown()
@Ref function Base.length(continuable::AbstractContinuable)
  i = Ref(0)
  continuable() do _
    i += 1
  end
  i
end

Base.IteratorEltype(::AbstractContinuable) = Base.EltypeUnknown()
Base.eltype(::AbstractContinuable) = Any

"""
    collect(continuable[, n]) -> Vector

Constructs Vector out of the given Continuable. If also given the length `n` explicitly,
a Vector of respective is preallocated. IMPORTANTLY `n` needs to be the true length of the continuable.
Smaller `n` will result in error.
"""
function Base.collect(c::AbstractContinuable)
  everything = Vector(undef, 0)
  reduce!(push!, c, init = everything)
end

@Ref function Base.collect(c::AbstractContinuable, n)
  a = Vector(undef, n)
  # unfortunately the nested call of enumerate results in slower code, hence we have a manual index here
  # this is so drastically that for small `n` a preallocate version with enumerate would be slower than the non-preallocate version
  i = Ref(1)
  c() do x
    a[i] = x
    i += 1
  end
  a
end

"""
    enumerate(continuable)

Constructs new Continuable with elements `(i, x)` for each `x` in the continuable, where `i` starts at `1` and
increments by `1` for each element.
"""
@cont @Ref function Base.enumerate(continuable::AbstractContinuable)
  i = Ref(1)
  continuable() do x
    cont((i, x))
    i += 1
  end
end

"""
    map(func, continuable)

Constructs new Continuable where the given `func` was applied to each element.
"""
Base.map(func, continuable::AbstractContinuable) = @cont continuable(x -> cont(func(x)))

"""
    filter(predicate, continuable)

Constructs new Continuable where only elements `x` with `predicate(x) == true` are kept.
"""
@cont function Base.filter(bool, continuable::AbstractContinuable)
  continuable() do x
    if bool(x)
      cont(x)
    end
  end
end

"""
    reduce(operator, continuable; [init])

Like Base.reduce this will apply `operator` iteratively to combine all elements into one accumulated result.
"""
Base.reduce(op, continuable::AbstractContinuable; init = nothing) = foldl_continuable(op, continuable, init)

"""
    reduce(operator, continuable; [init])

Like Base.foldl this will apply `operator` iteratively to combine all elements into one accumulated result.
The order is guaranteed to be left to right.
"""
Base.foldl(op, continuable::AbstractContinuable; init = nothing) = foldl_continuable(op, continuable, init)

struct EmptyStart end
@Ref function foldl_continuable(op, continuable, init::Nothing)
  acc = Ref{Any}(EmptyStart())
  lifted_op(acc::EmptyStart, x) = x
  lifted_op(acc, x) = op(acc, x)
  continuable() do x
    acc = lifted_op(acc, x)
  end
  acc
end
@Ref function foldl_continuable(op, continuable, init)
  acc = Ref(init)
  continuable() do x
    acc = op(acc, x)
  end
  acc
end

"""
    reduce!(op!, continuable; [init])

Mutating version of Base.reduce

If no `init` is given
    `op!` is assumed to mutate a hidden state (equivalent to mere continuation)
else
    `init` is the explicit state and will be passed to `op!` as first argument (the accumulator)
"""
function reduce!(op!, continuable::AbstractContinuable; init = nothing)
  if isnothing(init)
    continuable(op!)
  else
    reduce!(op!, continuable, init)
  end
end
function reduce!(op!, continuable::AbstractContinuable, acc)
  continuable() do x
    op!(acc, x)
  end
  acc
end

"""
    sum(continuable)

sums up all elements
"""
Base.sum(c::AbstractContinuable) = reduce(+, c)

"""
    sum(continuable)

multiplies up all elements
"""
Base.prod(c::AbstractContinuable) = reduce(*, c)

"""
    all([func, ]continuable; [lazy])

Checks whether all elements in the continuable are true.
If a function `func` is given, it is first applied to the elements before comparing for truth.

If `lazy=true` (default) the Continuable will only be evaluated until the first `false` value.
Elseif `lazy=false` all elements of the Continuable will be combined.
"""
@Ref function Base.all(continuable::AbstractContinuable; lazy=true)
  if lazy
    stoppable(continuable, true) do b
      if !b
        stop(false)
      end
    end
  else # non-lazy
    b = Ref(true)
    continuable() do x
      b &= x
    end
    b
  end
end
Base.all(f, continuable::AbstractContinuable; kwargs...) = all(map(f, continuable); kwargs...)

"""
    any([func, ]continuable; [lazy])

Checks whether at least one element in the continuable is true.
If a function `func` is given, it is first applied to the elements before comparing for truth.

If `lazy=true` (default) the Continuable will only be evaluated until the first `true` value.
Elseif `lazy=false` all elements of the Continuable will be combined.
"""
@Ref function Base.any(continuable::AbstractContinuable; lazy=true)
  if lazy
    stoppable(continuable, false) do b
      if b
        stop(true)
      end
    end

  else # non-lazy
    b = Ref(false)
    continuable() do x
      b |= x
    end
    b
  end
end
Base.any(f, continuable::AbstractContinuable; kwargs...) = any(map(f, continuable); kwargs...)

## zip ----------------------------
# zip is the only method which seems to be unimplementable with continuations
# hence we have to go to tasks or arrays

"""
    azip(continuables...)

Zipping continuables via intermediate array representation

CAUTION: loads everything into memory
"""
azip(cs::AbstractContinuable...) = @cont begin
  # not possible with continuations... bring it to memory and apply normal zip
  array_cs = collect.(cs)
  for t in zip(array_cs...)
    cont(t)
  end
end

"""
    chzip(continuables...)

Zipping continuables via Channel
"""
chzip(cs::AbstractContinuable...) = @cont begin
  # or use aschannel and iterate
  channel_cs = aschannel.(cs)
  for t in zip(channel_cs...)
    cont(t)
  end
end

"""
    zip(continuables...; [lazy])

Constructs new Continuable with elements from the given continuables zipped up.
I.e. will yield for each position in the original continuables a tuple `(x, y, ...)`
where `x`, `y`, ... are the elements from `continuables` at the same position respectively.

If `lazy=true` (default), it will use Channels to do the zipping.
Elseif `lazy=false`, it will use Arrays instead.

!!! warning CAUTION
    `zip` on Continuables is not performant, but will fallback to either Channels (`lazy=true`, default) which are
    very slow, or Arrays (`lazy=false`) which will load everything into Memory.
"""
function Base.zip(cs::AbstractContinuable...; lazy=true)
  if lazy
    chzip(cs...)
  else
    azip(cs...)
  end
end

"""
    cycle(continuable[, n])

Constructs new Continuable which loops through the given continuable.
If `n` is given, it will loop `n` times, otherwise endlessly.
"""
cycle(continuable::AbstractContinuable) = @cont while true
  continuable(cont)
end

cycle(continuable::AbstractContinuable, n::Integer) = @cont for _ in 1:n
  continuable(cont)
end

## combine continuables --------------------------------------------

# IMPORTANT we cannot overload the empty product as it would conflict with iterables

"""
    product(continuables...)

Construct a new Continuable which yields all combinations of the given continuables, analog
to how Iterators.product work for iterables.

Mind that `product()` will still return an empty iterator instead of an empty Continuable.
Use [`emptycontinuable`](@ref) instead if you need an empty Continuable.
"""
product(c1::AbstractContinuable) = c1

# note this function in fact returns a continuable, however it is written as highlevel as that no explicit "f(...) = cont -> begin ... end" is needed
product(c1::AbstractContinuable, c2::AbstractContinuable) = @cont begin
  c1() do x
    c2() do y
      cont((x,y))
    end
  end
end

# this method is underscored because we assume the first continuation to deliver tuples and not values
_product(c1::AbstractContinuable, c2::AbstractContinuable) = @cont begin
  c1() do t
    c2() do x
      cont(tuple(t..., x))
    end
  end
end

@Ref function product(c1::AbstractContinuable, c2::AbstractContinuable, cs::Vararg{<:AbstractContinuable})
  acc = Ref{Any}(product(c1, c2))  # make first into singleton tuple to start recursion
  for continuable in cs
    acc = _product(acc, continuable)
  end
  acc
end

"""
    flatten(continuable_of_continuables)

Constructs new Continuable by concatinating all continuables in the given `continuable_of_continuables`.
Analog to Iterators.flatten.

For iterables of continuable use `Continuables.chain(iterable_of_continuables...)` instead.
"""
@cont function flatten(continuable::AbstractContinuable)
  continuable() do subcontinuable
    subcontinuable(cont)
  end
end

"""
    chain(continuables...)
    chain(iterables...) = flatten(iterables)

When given Continuables it will construct a new continuable by concatinating all given continuables.
When given anything else it will default to use `Iterator.flatten`.
"""
chain(iterables::Vararg) = flatten(iterables)
chain(cs::Vararg{<:AbstractContinuable}) = @cont begin
  for continuable in cs
    continuable(cont)
  end
end

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

"""
    take(continuable, n)
    take(n, continuable)

Construct a new Continuable which only yields the first `n` elements.
`n` can be larger as the total length, no problem.
"""
@cont @Ref function take(continuable::AbstractContinuable, n::Integer)
  i = Ref(0)
  stoppable(continuable) do x
    i += 1
    if i > n
      stop()
    end
    cont(x)
  end
end
take(n::Integer, continuable::AbstractContinuable) = take(continuable, n)

"""
    takewhile(predicate, continuable)
    takewhile(predicate, iterable)

If given a Continuable, it constructs a new Continuable yielding elements until `predicate(element)` returns `false`.

Also implements a respective functionality for iterables for convenience.
"""
takewhile(bool, iterable) = TakeWhile(bool,  iterable)
@cont function takewhile(bool, continuable::AbstractContinuable)
  stoppable(continuable) do x
    if !bool(x)
      stop()
    end
    cont(x)
  end
end

"""
    drop(continuable, n)
    drop(n, continuable)

Construct a new Continuable which yields all elements but the first `n`.
`n` can be larger as the total length, no problem.
"""
@cont @Ref function drop(continuable::AbstractContinuable, n::Integer)
  i = Ref(0)
  continuable() do x
    i += 1
    if i > n
      cont(x)
    end
  end
end
drop(n::Integer, continuable::AbstractContinuable) = drop(continuable, n)

"""
    dropwhile(predicate, continuable)
    dropwhile(predicate, iterable)

If given a Continuable, it constructs a new Continuable yielding elements until `predicate(element)` returns `true`.

Also implements a respective functionality for iterables for convenience.
"""
dropwhile(bool, iterable) = DropWhile(bool, iterable)
@cont @Ref function dropwhile(bool, continuable::AbstractContinuable)
  dropping = Ref(true)
  continuable() do x
    if dropping
      dropping &= bool(x)
      # without nested "if" statement we would have to use two separate if statements at the top (instead of using if else)
      !dropping && cont(x)
    else
      cont(x)
    end
  end
end

"""
    partition(continuable, n[, step])

Constructs new Continuable which yields whole subsections of the given continuable, gathered as Vectors.
`n` is the length of a subsection. The very last subsection might be of length `n` or smaller respectively, collecting
the remaining elements.

If `step` is given, the second subsection is exactly `step`-number of elements apart from the previous subsection,
and hence overlapping if `n > step`.
Further, importantly, if `step` is given, there is no rest, but each subsection will be guaranteed to have the same
length. This semantics is copied from [IterTools.jl](https://juliacollections.github.io/IterTools.jl/latest/#partition(xs,-n,-[step])-1)

# Examples
```jldoctest
julia> using Continuables

julia> partition(i2c(1:10), 3) |> collect
4-element Array{Any,1}:
 Any[1, 2, 3]
 Any[4, 5, 6]
 Any[7, 8, 9]
 Any[10]
julia> partition(i2c(1:10), 5, 2) |> collect
3-element Array{Any,1}:
 Any[1, 2, 3, 4, 5]
 Any[3, 4, 5, 6, 7]
 Any[5, 6, 7, 8, 9]
julia> partition(i2c(1:10), 3, 3) |> collect
3-element Array{Any,1}:
 Any[1, 2, 3]
 Any[4, 5, 6]
 Any[7, 8, 9]
```
"""
@cont @Ref function partition(continuable::AbstractContinuable, n::Integer)
  i = Ref(1)
  part = Ref(Vector(undef, n))
  continuable() do x
    part[i] = x
    i += 1
    if i > n
      cont(part)
      part = Vector(undef, n)
      i = 1
    end
  end
  # final bit # TODO is this wanted? with additional step parameter I think this is mostly unwanted
  if i > 1
    # following the implementation for iterable, we cut the length to the defined part
    cont(_takewhile_isassigned(part))
  end
end

@cont @Ref function partition(continuable::AbstractContinuable, n::Integer, step::Integer)
  i = Ref(0)
  n_overlap = n - step
  part = Ref(Vector(undef, n))
  continuable() do x
    i += 1
    if i > 0  # if i is negative we simply skip these
      part[i] = x
    end
    if i == n
      cont(part)
      if n_overlap > 0
        overlap = part[1+step:n]
        part = Vector(undef, n)
        part[1:n_overlap] = overlap
      else
        # we need to recreate new part because of references
        part = Vector(undef, n)
      end
      i = n_overlap
    end
  end
end

function _takewhile_isassigned(vec::Vector)
  n = length(vec)
  for i in 1:n
    if !isassigned(vec, i)
      return vec[1:(i-1)]
    end
  end
  return vec
end

"""
    groupbyreduce(by, continuable, op2[, op1])
    groupbyreduce(by, iterable, op2[, op1])

Group elements and returns OrderedDict of keys (constructed by `by`) and values (aggregated with `op2`/`op1`)
If given anything else then a continuable, we interpret it as an iterable and provide the same functionality.

# Parameters

by: function of element to return the key for the grouping/dict
continuable: will get grouped
op2: f(accumulator, element) = new_accumulator
op1: f(element) = initial_accumulator

# Examples
```jldoctest
julia> using Continuables

julia> groupbyreduce(x -> x % 4, @i2c(1:10), (x, y) -> x + y)
OrderedCollections.OrderedDict{Any,Any} with 4 entries:
  1 => 15
  2 => 18
  3 => 10
  0 => 12
julia> groupbyreduce(x -> x % 4, @i2c(1:10), (x, y) -> x + y, x -> x+5)
OrderedCollections.OrderedDict{Any,Any} with 4 entries:
  1 => 20
  2 => 23
  3 => 15
  0 => 17
```
"""
function groupbyreduce(by, continuable::AbstractContinuable, op2, op1=identity)
  d = OrderedDict()
  continuable() do x
    key = by(x)
    if key in keys(d)
      d[key] = op2(d[key], x)
    else
      d[key] = op1(x)
    end
  end
  d
end
# adding iterable versions for the general case (tests showed that these are actually compiling to the iterable version in terms of code and speed, awesome!)
groupbyreduce(by, iterable, op2, op1=identity) = groupbyreduce(by, ascontinuable(iterable), op2, op1)

"""
    groupby(f, continuable)
    groupby(f, iterable)

Wrapper around the more general [`groupbyreduce`](@ref) which combines elements to a Vector.
If you happen to aggregate your resulting grouped Vectors, think about using `groupbyreduce` directly, as
this can massively speed up aggregations.

Note that the interface is different from `IterTools.groupby`, as we directly return an OrderedDict
(instead of a iterable of values).

# Examples

```jldoctest
julia> using Continuables

julia> groupby(x -> x % 4, @i2c 1:10)
OrderedCollections.OrderedDict{Any,Any} with 4 entries:
  1 => [1, 5, 9]
  2 => [2, 6, 10]
  3 => [3, 7]
  0 => [4, 8]
```
"""
groupby(f, continuable::AbstractContinuable) = groupbyreduce(f, continuable, push!, x -> [x])
groupby(f, iterable) = groupby(f, ascontinuable(iterable))

## subsets & peekiter -------------------------------------------------------

# subsets seem to be implemented for arrays in the first place (and not iterables in general)
# hence better use IterTools.subsets directly

# peekiter is the only method of Iterators.jl missing. However it in fact makes no sense for continuables
# as they are functions and don't get consumed

## extract values from continuables  ----------------------------------------

"""
    nth(continuable, n)
    nth(n, continuable)

Extracts the `n`th element from the given continuable.

# Examples
```jldoctest
julia> using Continuables

julia> nth(i2c(4:10), 3)
6
julia> nth(1, i2c(4:10))
4
```
"""
@Ref function nth(continuable::AbstractContinuable, n::Integer)
  i = Ref(0)
  ret = stoppable(continuable) do x
    i += 1
    if i==n
      # CAUTION: we cannot use return here as usual because this is a subfunction. Return works here more like continue
      stop(x)
    end
  end
  if ret === nothing
    error("given continuable has length $i < $n")
  end
  ret
end
nth(n::Integer, continuable::AbstractContinuable) = nth(continuable, n)

end  # module


text:
@doc doc"""
    MisclassLoss <: SupervisedLoss

Misclassification loss that assigns `1` for misclassified
examples and `0` otherwise. It is a generalization of
`ZeroOneLoss` for more than two classes.
"""
struct MisclassLoss <: SupervisedLoss end

agreement(target, output) = target == output

value(::MisclassLoss, agreement::Bool) = agreement ? 0 : 1
deriv(::MisclassLoss, agreement::Bool) = 0
deriv2(::MisclassLoss, agreement::Bool) = 0
value_deriv(::MisclassLoss, agreement::Bool) = agreement ? (0, 0) : (1, 0)

value(loss::MisclassLoss, target::Number, output::Number) = value(loss, agreement(target, output))
deriv(loss::MisclassLoss, target::Number, output::Number) = deriv(loss, agreement(target, output))
deriv2(loss::MisclassLoss, target::Number, output::Number) = deriv2(loss, agreement(target, output))

isminimizable(::MisclassLoss) = false
isdifferentiable(::MisclassLoss) = false
isdifferentiable(::MisclassLoss, at) = at != 0
istwicedifferentiable(::MisclassLoss) = false
istwicedifferentiable(::MisclassLoss, at) = at != 0
isnemitski(::MisclassLoss) = false
islipschitzcont(::MisclassLoss) = false
isconvex(::MisclassLoss) = false
isclasscalibrated(::MisclassLoss) = false
isclipable(::MisclassLoss) = false

# ===============================================================

@doc doc"""
    PoissonLoss <: SupervisedLoss

Loss under a Poisson noise distribution (KL-divergence)

``L(target, output) = exp(output) - target*output``
"""
struct PoissonLoss <: SupervisedLoss end

value(loss::PoissonLoss, target::Number, output::Number) = exp(output) - target*output
deriv(loss::PoissonLoss, target::Number, output::Number) = exp(output) - target
deriv2(loss::PoissonLoss, target::Number, output::Number) = exp(output)
function value_deriv(loss::PoissonLoss, target::Number, output::Number)
    exp_o = exp(output)
    return (exp_o-(target*output), exp_o-target)
end

isdifferentiable(::PoissonLoss) = true
isdifferentiable(::PoissonLoss, y, t) = true
istwicedifferentiable(::PoissonLoss) = true
istwicedifferentiable(::PoissonLoss, y, t) = true
islipschitzcont(::PoissonLoss) = false
isconvex(::PoissonLoss) = true
# TODO: isstrictlyconvex(::PoissonLoss) = ?
isstronglyconvex(::PoissonLoss) = false

# ===============================================================

@doc doc"""
    CrossEntropyLoss <: SupervisedLoss

Cross-entropy loss also known as log loss and logistic loss is defined as:

``L(target, output) = - target*ln(output) - (1-target)*ln(1-output)``
"""
struct CrossEntropyLoss <: SupervisedLoss end
const LogitProbLoss = CrossEntropyLoss

function value(loss::CrossEntropyLoss, target::Number, output::Number)
    target >= 0 && target <=1 || error("target must be in [0,1]")
    output >= 0 && output <=1 || error("output must be in [0,1]")
    if target == 0
        -log(1 - output)
    elseif target == 1
        -log(output)
    else
        -(target * log(output) + (1-target) * log(1-output))
    end
end
deriv(loss::CrossEntropyLoss, target::Number, output::Number) = (1-target) / (1-output) - target / output
deriv2(loss::CrossEntropyLoss, target::Number, output::Number) = (1-target) / (1-output)^2 + target / output^2
value_deriv(loss::CrossEntropyLoss, target::Number, output::Number) = (value(loss,target,output), deriv(loss,target,output))

isdifferentiable(::CrossEntropyLoss) = true
isdifferentiable(::CrossEntropyLoss, y, t) = t != 0 && t != 1
istwicedifferentiable(::CrossEntropyLoss) = true
istwicedifferentiable(::CrossEntropyLoss, y, t) = t != 0 && t != 1
isconvex(::CrossEntropyLoss) = true

# ===============================================================
# L(target, output) = sign(agreement) < 0 ? 1 : 0


text:
#compare to precomputed values
k0 = 0.3
kin = 0.5k0
tols = [1e-9;
        1e-7;
        1e-5]
N_res = [250;
         100;
         313]

@testset "minimum N" begin
    shape_funs =   [N -> squircle(1, N);
                    N -> rounded_star(0.08, 0.02, 3, N);
                    N -> rounded_star(25, 5, 5, N)]

    for (i,f) in enumerate(shape_funs)
        N, err = minimumN(k0, kin, f; tol = tols[i], N_points = 20_000,
                    N_start = 250, N_min = 100, N_max = 400)
        @test N == N_res[i] && err <= tols[i]
    end
end

P_res = [13;
         9;
         34]
dists = [2.0;
    	 2.0;
         1.0]
@testset "minimum P" begin
    shapes = [squircle(1, N_res[1]);
                rounded_star(0.08, 0.02, 3, N_res[2]);
                rounded_star(25, 5, 5, N_res[3])]

    for i in eachindex(shapes)
        P, err = minimumP(k0, kin, shapes[i]; tol = tols[i], N_points = 20_000,
                    P_min = 1, P_max = 60, dist = dists[i])
        @test P == P_res[i] && err <= tols[i]
    end
end


text:
@time @safetestset "Prerequisites" begin
    @time @safetestset "Tools" begin include("Tools.jl") end
    @time @safetestset "Combinatorics" begin include("Combinatorics.jl") end
    @time @safetestset "Traits" begin include("Traits.jl") end
    @time @safetestset "CompositeStructures" begin include("CompositeStructures.jl") end
    @time @safetestset "SimpleTrees" begin include("SimpleTrees.jl") end
    @time @safetestset "NamedVectors" begin include("NamedVectors.jl") end
    @time @safetestset "VectorSpaces" begin include("VectorSpaces.jl") end
end


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Tests related to the NRLMSISE-00 Atmospheric Model.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] https://ccmc.gsfc.nasa.gov/pub/modelweb/atmospheric/msis/nrlmsise00/nrlmsise00_output.txt
#   [2] https://ccmc.gsfc.nasa.gov/modelweb/models/nrlmsise00.php
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

# File: ./src/earth/atmospheric_models/nrlmsise00
# ===============================================

# Functions: gtd7 and gtd7d
# -------------------------

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Tests available in the text file `output.txt` at [1].
#
#  DAY           172          81         172         172         172
#  UT         29000.      29000.      75000.      29000.      29000.
#  ALT          400.        400.       1000.        100.        400.
#  LAT           60.         60.         60.         60.          0.
#  LONG         -70.        -70.        -70.        -70.        -70.
#  LST           16.         16.         16.         16.         16.
#  F107A        150.        150.        150.        150.        150.
#  F107         150.        150.        150.        150.        150.
#  AP             4.          4.          4.          4.          4.
#
#  TINF      1250.54     1166.75     1239.89     1027.32     1212.40
#  TG        1241.42     1161.71     1239.89      206.89     1208.14
#  HE      6.665E+05   3.407E+06   1.124E+05   5.412E+07   1.851E+06
#  O       1.139E+08   1.586E+08   6.934E+04   1.919E+11   1.477E+08
#  N2      1.998E+07   1.391E+07   4.247E+01   6.116E+12   1.579E+07
#  O2      4.023E+05   3.263E+05   1.323E-01   1.225E+12   2.634E+05
#  AR      3.557E+03   1.560E+03   2.619E-05   6.023E+10   1.589E+03
#  H       3.475E+04   4.854E+04   2.017E+04   1.060E+07   5.816E+04
#  N       4.096E+06   4.381E+06   5.741E+03   2.616E+05   5.479E+06
#  ANM O   2.667E+04   6.957E+03   2.374E+04   0.000E+00   1.264E+03
#  RHO     4.075E-15   5.002E-15   2.757E-18   3.584E-10   4.810E-15
#
#
#  DAY           172         172         172         172         172
#  UT         29000.      29000.      29000.      29000.      29000.
#  ALT          400.        400.        400.        400.        400.
#  LAT           60.         60.         60.         60.         60.
#  LONG           0.        -70.        -70.        -70.        -70.
#  LST           16.          4.         16.         16.         16.
#  F107A        150.        150.         70.        150.        150.
#  F107         150.        150.        150.        180.        150.
#  AP             4.          4.          4.          4.         40.
#
#  TINF      1220.15     1116.39     1031.25     1306.05     1361.87
#  TG        1212.71     1113.00     1024.85     1293.37     1347.39
#  HE      8.673E+05   5.776E+05   3.740E+05   6.748E+05   5.529E+05
#  O       1.279E+08   6.979E+07   4.783E+07   1.245E+08   1.198E+08
#  N2      1.823E+07   1.237E+07   5.240E+06   2.369E+07   3.496E+07
#  O2      2.922E+05   2.493E+05   1.760E+05   4.912E+05   9.340E+05
#  AR      2.403E+03   1.406E+03   5.502E+02   4.579E+03   1.096E+04
#  H       3.686E+04   5.292E+04   8.897E+04   3.245E+04   2.686E+04
#  N       3.897E+06   1.070E+06   1.980E+06   5.371E+06   4.890E+06
#  ANM O   2.667E+04   2.667E+04   9.122E+03   2.667E+04   2.805E+04
#  RHO     4.356E-15   2.471E-15   1.572E-15   4.564E-15   4.975E-15
#
#
#  DAY           172         172         172         172         172
#  UT         29000.      29000.      29000.      29000.      29000.
#  ALT            0.         10.         30.         50.         70.
#  LAT           60.         60.         60.         60.         60.
#  LONG         -70.        -70.        -70.        -70.        -70.
#  LST           16.         16.         16.         16.         16.
#  F107A        150.        150.        150.        150.        150.
#  F107         150.        150.        150.        150.        150.
#  AP             4.          4.          4.          4.          4.
#
#  TINF      1027.32     1027.32     1027.32     1027.32     1027.32
#  TG         281.46      227.42      237.44      279.56      219.07
#  HE      1.375E+14   4.427E+13   2.128E+12   1.412E+11   1.255E+10
#  O       0.000E+00   0.000E+00   0.000E+00   0.000E+00   0.000E+00
#  N2      2.050E+19   6.598E+18   3.171E+17   2.104E+16   1.875E+15
#  O2      5.499E+18   1.770E+18   8.506E+16   5.645E+15   4.923E+14
#  AR      2.452E+17   7.892E+16   3.793E+15   2.517E+14   2.240E+13
#  H       0.000E+00   0.000E+00   0.000E+00   0.000E+00   0.000E+00
#  N       0.000E+00   0.000E+00   0.000E+00   0.000E+00   0.000E+00
#  ANM O   0.000E+00   0.000E+00   0.000E+00   0.000E+00   0.000E+00
#  RHO     1.261E-03   4.059E-04   1.951E-05   1.295E-06   1.148E-07
#
################################################################################

@testset "Function gtd7d" begin
    # Scenario 01
    # ===========

    test_list = ["nrlmsise00_test01.txt",
                 "nrlmsise00_test02.txt",
                 "nrlmsise00_test03.txt"]

    for file in test_list
        raw = readdlm(file)

        # Loop for each test case.
        for i = 2:size(raw,2)
            # Parse the inputs.
            doy       = raw[1,i]
            sec       = raw[2,i]
            alt       = raw[3,i]
            g_lat     = raw[4,i]
            g_long    = raw[5,i]
            lst       = raw[6,i]
            f107A     = raw[7,i]
            f107      = raw[8,i]
            ap        = raw[9,i]

            # Parse the outputs.
            tinf      = raw[10,i]
            tg        = raw[11,i]
            den_He    = raw[12,i]
            den_O     = raw[13,i]
            den_N2    = raw[14,i]
            den_O2    = raw[15,i]
            den_Ar    = raw[16,i]
            den_H     = raw[17,i]
            den_N     = raw[18,i]
            den_aO    = raw[19,i]
            den_Total = raw[20,i]

            # Configure and run NRLMSISE-00.
            nrlmsise00d = conf_nrlmsise00(0,
                                          doy,
                                          sec,
                                          alt,
                                          g_lat,
                                          g_long,
                                          lst,
                                          f107A,
                                          f107,
                                          ap)

            out = gtd7(nrlmsise00d)

            # Test the results.
            @test out.T_exo     ≈ tinf      rtol = 1e-3 atol = 1e-9
            @test out.T_alt     ≈ tg        rtol = 1e-3 atol = 1e-9
            @test out.den_He    ≈ den_He    rtol = 1e-3 atol = 1e-9
            @test out.den_O     ≈ den_O     rtol = 1e-3 atol = 1e-9
            @test out.den_N2    ≈ den_N2    rtol = 1e-3 atol = 1e-9
            @test out.den_O2    ≈ den_O2    rtol = 1e-3 atol = 1e-9
            @test out.den_Ar    ≈ den_Ar    rtol = 1e-3 atol = 1e-9
            @test out.den_H     ≈ den_H     rtol = 1e-3 atol = 1e-9
            @test out.den_N     ≈ den_N     rtol = 1e-3 atol = 1e-9
            @test out.den_aO    ≈ den_aO    rtol = 1e-3 atol = 1e-9

            @test out.den_Total ≈ den_Total rtol = 1e-3

            # Compare the results with `gtd7d`, in which the only modification
            # is that the density of anomalous oxygen is summed in `den_Total`.

            outd = gtd7d(nrlmsise00d)

            @test outd.T_exo     ≈ tinf      rtol = 1e-3 atol = 1e-9
            @test outd.T_alt     ≈ tg        rtol = 1e-3 atol = 1e-9
            @test outd.den_He    ≈ den_He    rtol = 1e-3 atol = 1e-9
            @test outd.den_O     ≈ den_O     rtol = 1e-3 atol = 1e-9
            @test outd.den_N2    ≈ den_N2    rtol = 1e-3 atol = 1e-9
            @test outd.den_O2    ≈ den_O2    rtol = 1e-3 atol = 1e-9
            @test outd.den_Ar    ≈ den_Ar    rtol = 1e-3 atol = 1e-9
            @test outd.den_H     ≈ den_H     rtol = 1e-3 atol = 1e-9
            @test outd.den_N     ≈ den_N     rtol = 1e-3 atol = 1e-9
            @test outd.den_aO    ≈ den_aO    rtol = 1e-3 atol = 1e-9

            @test outd.den_Total ≈ den_Total + 1.66e-24*16den_aO rtol = 1e-3
        end
    end
end

# Function: nrlmsise00
# ---------------------

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Simulation using only the daily AP.
#
# Data obtained from the online version of NRLMSISE-00 [2]:
#
# Input parameters
#
# year= 1986, month= 6, day= 19, hour=21.50,
# Time_type = Universal
# Coordinate_type = Geographic
# latitude= -16.00, longitude= 312.00, height= 100.00
# Prof. parameters: start= 0.00 stop= 1000.00 step= 100.00
#
# Optional parametes: F10.7(daily) =121.; F10.7(3-month avg) =80.; ap(daily) = 7.
#
#    Selected parameters are:
# 1 Height, km
# 2 O, cm-3
# 3 N2, cm-3
# 4 O2, cm-3
# 5 Mass_density, g/cm-3
# 6 Temperature_neutral, K
# 7 Temperature_exospheric, K
# 8 He, cm-3
# 9 Ar, cm-3
# 10 H, cm-3
# 11 N, cm-3
# 12 Anomalous_Oxygen, cm-3
#
#       1          2          3          4         5      6     7          8          9         10         11         12
#     0.0  0.000E+00  1.918E+19  5.145E+18 1.180E-03  297.7  1027  1.287E+14  2.294E+17  0.000E+00  0.000E+00  0.000E+00
#   100.0  4.244E+11  9.498E+12  2.240E+12 5.783E-10  165.9  1027  1.061E+08  9.883E+10  2.209E+07  3.670E+05  0.000E+00
#   200.0  2.636E+09  2.248E+09  1.590E+08 1.838E-13  829.3   909  5.491E+06  1.829E+06  2.365E+05  3.125E+07  1.802E-09
#   300.0  3.321E+08  6.393E+07  2.881E+06 1.212E-14  900.6   909  3.173E+06  1.156E+04  1.717E+05  6.708E+06  4.938E+00
#   400.0  5.088E+07  2.414E+06  6.838E+04 1.510E-15  907.7   909  1.979E+06  1.074E+02  1.518E+05  1.274E+06  1.237E+03
#   500.0  8.340E+06  1.020E+05  1.839E+03 2.410E-16  908.5   909  1.259E+06  1.170E+00  1.355E+05  2.608E+05  4.047E+03
#   600.0  1.442E+06  4.729E+03  5.500E+01 4.543E-17  908.6   909  8.119E+05  1.455E-02  1.214E+05  5.617E+04  4.167E+03
#   700.0  2.622E+05  2.394E+02  1.817E+00 1.097E-17  908.6   909  5.301E+05  2.050E-04  1.092E+05  1.264E+04  3.173E+03
#   800.0  4.999E+04  1.317E+01  6.608E-02 3.887E-18  908.6   909  3.503E+05  3.254E-06  9.843E+04  2.964E+03  2.246E+03
#   900.0  9.980E+03  7.852E-01  2.633E-03 1.984E-18  908.6   909  2.342E+05  5.794E-08  8.900E+04  7.237E+02  1.571E+03
#  1000.0  2.082E+03  5.055E-02  1.146E-04 1.244E-18  908.6   909  1.582E+05  1.151E-09  8.069E+04  1.836E+02  1.103E+03
#
# Scenario 02
# ===========
#
# Simulation using the AP vector.
#
# Data obtained from the online version of NRLMSISE-00 [2]:
#
# Input parameters
# year= 2016, month= 6, day= 1, hour=11.00,
# Time_type = Universal
# Coordinate_type = Geographic
# latitude= -23.00, longitude= 315.00, height= 100.00
# Prof. parameters: start= 0.00 stop= 1000.00 step= 100.00
#
# Optional parametes: F10.7(daily) =not specified; F10.7(3-month avg) =not specified; ap(daily) = not specified
#
#    Selected parameters are:
# 1 Height, km
# 2 O, cm-3
# 3 N2, cm-3
# 4 O2, cm-3
# 5 Mass_density, g/cm-3
# 6 Temperature_neutral, K
# 7 Temperature_exospheric, K
# 8 He, cm-3
# 9 Ar, cm-3
# 10 H, cm-3
# 11 N, cm-3
# 12 Anomalous_Oxygen, cm-3
# 13 F10_7_daily
# 14 F10_7_3_month_avg
# 15 ap_daily
# 16 ap_00_03_hr_prior
# 17 ap_03_06_hr_prior
# 18 ap_06_09_hr_prior
# 19 ap_09_12_hr_prior
# 20 ap_12_33_hr_prior
# 21 ap_33_59_hr_prior
#
#       1          2          3          4         5      6     7          8          9         10         11         12      13     14     15      16     17     18      19     20     21
#     0.0  0.000E+00  1.919E+19  5.149E+18 1.181E-03  296.0  1027  1.288E+14  2.296E+17  0.000E+00  0.000E+00  0.000E+00   89.0   89.4    4.9    3.0    4.0    4.0    6.0    9.9   10.6
#   100.0  6.009E+11  9.392E+12  2.213E+12 5.764E-10  165.2  1027  1.339E+08  9.580E+10  2.654E+07  2.737E+05  0.000E+00   89.0   89.4    4.9    3.0    4.0    4.0    6.0    9.9   10.6
#   200.0  3.882E+09  1.978E+09  1.341E+08 2.028E-13  687.5   706  2.048E+07  9.719E+05  5.172E+05  1.632E+07  1.399E-09   89.0   89.4    4.9    3.0    4.0    4.0    6.0    9.9   10.6
#   300.0  3.140E+08  2.481E+07  9.340E+05 9.666E-15  705.5   706  1.082E+07  1.879E+03  3.860E+05  2.205E+06  3.876E+00   89.0   89.4    4.9    3.0    4.0    4.0    6.0    9.9   10.6
#   400.0  2.855E+07  3.737E+05  7.743E+03 8.221E-16  706.0   706  5.938E+06  4.688E+00  3.317E+05  2.633E+05  9.738E+02   89.0   89.4    4.9    3.0    4.0    4.0    6.0    9.9   10.6
#   500.0  2.787E+06  6.372E+03  7.382E+01 9.764E-17  706.0   706  3.319E+06  1.397E-02  2.868E+05  3.424E+04  3.186E+03   89.0   89.4    4.9    3.0    4.0    4.0    6.0    9.9   10.6
#   600.0  2.910E+05  1.222E+02  8.047E-01 2.079E-17  706.0   706  1.887E+06  4.919E-05  2.490E+05  4.742E+03  3.281E+03   89.0   89.4    4.9    3.0    4.0    4.0    6.0    9.9   10.6
#   700.0  3.240E+04  2.622E+00  9.973E-03 8.474E-18  706.0   706  1.090E+06  2.034E-07  2.171E+05  6.946E+02  2.498E+03   89.0   89.4    4.9    3.0    4.0    4.0    6.0    9.9   10.6
#   800.0  3.835E+03  6.264E-02  1.398E-04 4.665E-18  706.0   706  6.393E+05  9.809E-10  1.900E+05  1.074E+02  1.768E+03   89.0   89.4    4.9    3.0    4.0    4.0    6.0    9.9   10.6
#   900.0  4.816E+02  1.659E-03  2.204E-06 2.817E-18  706.0   706  3.806E+05  5.481E-12  1.669E+05  1.747E+01  1.236E+03   89.0   89.4    4.9    3.0    4.0    4.0    6.0    9.9   10.6
#  1000.0  6.400E+01  4.853E-05  3.891E-08 1.772E-18  706.0   706  2.298E+05  3.528E-14  1.471E+05  2.988E+00  8.675E+02   89.0   89.4    4.9    3.0    4.0    4.0    6.0    9.9   10.6
#
#
################################################################################

@testset "Function nrlmsise00" begin
    ## Scenario 01
    ## ===========

    # Test outputs.
    test_out = [
       0.0  0.000E+00  1.918E+19  5.145E+18 1.180E-03  297.7  1027  1.287E+14 2.294E+17  0.000E+00  0.000E+00  0.000E+00;
     100.0  4.244E+11  9.498E+12  2.240E+12 5.783E-10  165.9  1027  1.061E+08 9.883E+10  2.209E+07  3.670E+05  0.000E+00;
     200.0  2.636E+09  2.248E+09  1.590E+08 1.838E-13  829.3   909  5.491E+06 1.829E+06  2.365E+05  3.125E+07  1.802E-09;
     300.0  3.321E+08  6.393E+07  2.881E+06 1.212E-14  900.6   909  3.173E+06 1.156E+04  1.717E+05  6.708E+06  4.938E+00;
     400.0  5.088E+07  2.414E+06  6.838E+04 1.510E-15  907.7   909  1.979E+06 1.074E+02  1.518E+05  1.274E+06  1.237E+03;
     500.0  8.340E+06  1.020E+05  1.839E+03 2.410E-16  908.5   909  1.259E+06 1.170E+00  1.355E+05  2.608E+05  4.047E+03;
     600.0  1.442E+06  4.729E+03  5.500E+01 4.543E-17  908.6   909  8.119E+05 1.455E-02  1.214E+05  5.617E+04  4.167E+03;
     700.0  2.622E+05  2.394E+02  1.817E+00 1.097E-17  908.6   909  5.301E+05 2.050E-04  1.092E+05  1.264E+04  3.173E+03;
     800.0  4.999E+04  1.317E+01  6.608E-02 3.887E-18  908.6   909  3.503E+05 3.254E-06  9.843E+04  2.964E+03  2.246E+03;
     900.0  9.980E+03  7.852E-01  2.633E-03 1.984E-18  908.6   909  2.342E+05 5.794E-08  8.900E+04  7.237E+02  1.571E+03;
    1000.0  2.082E+03  5.055E-02  1.146E-04 1.244E-18  908.6   909  1.582E+05 1.151E-09  8.069E+04  1.836E+02  1.103E+03;
    ]

    # Constant input parameters.
    year   = 1986
    month  = 6
    day    = 19
    hour   = 21
    minute = 30
    second = 00
    g_lat  = -16*pi/180
    g_long = 312*pi/180
    f107   = 121
    f107A  = 80
    ap     = 7
    JD     = date_to_jd(year, month, day, hour, minute, second)

    for i = 1:size(test_out,1)
        # Run the NRLMSISE-00 model wih the input parameters.
        out = nrlmsise00(JD,
                         test_out[i,1]*1000,
                         g_lat,
                         g_long,
                         f107A,
                         f107,
                         ap;
                         output_si = false,
                         dversion = false)

        @test out.den_O     ≈ test_out[i, 2]  rtol = 1e-3 atol = 1e-9
        @test out.den_N2    ≈ test_out[i, 3]  rtol = 1e-3 atol = 1e-9
        @test out.den_O2    ≈ test_out[i, 4]  rtol = 1e-3 atol = 1e-9
        @test out.T_alt     ≈ test_out[i, 6]  rtol = 1e-1 atol = 1e-9
        @test out.T_exo     ≈ test_out[i, 7]  rtol = 1e-1 atol = 1e-9
        @test out.den_He    ≈ test_out[i, 8]  rtol = 1e-3 atol = 1e-9
        @test out.den_Ar    ≈ test_out[i, 9]  rtol = 1e-3 atol = 1e-9
        @test out.den_H     ≈ test_out[i,10]  rtol = 1e-3 atol = 1e-9
        @test out.den_N     ≈ test_out[i,11]  rtol = 1e-3 atol = 1e-9
        @test out.den_aO    ≈ test_out[i,12]  rtol = 1e-3 atol = 1e-9
        @test out.den_Total ≈ test_out[i, 5]  rtol = 1e-3

        # Test the version that calls `gtd7d` instead of `gtd7`.
        outd = nrlmsise00(JD,
                          test_out[i,1]*1000,
                          g_lat,
                          g_long,
                          f107A,
                          f107,
                          ap;
                          output_si = false,
                          dversion = true)

        @test outd.den_O     ≈ test_out[i, 2]  rtol = 1e-3 atol = 1e-9
        @test outd.den_N2    ≈ test_out[i, 3]  rtol = 1e-3 atol = 1e-9
        @test outd.den_O2    ≈ test_out[i, 4]  rtol = 1e-3 atol = 1e-9
        @test outd.T_alt     ≈ test_out[i, 6]  rtol = 1e-1 atol = 1e-9
        @test outd.T_exo     ≈ test_out[i, 7]  rtol = 1e-1 atol = 1e-9
        @test outd.den_He    ≈ test_out[i, 8]  rtol = 1e-3 atol = 1e-9
        @test outd.den_Ar    ≈ test_out[i, 9]  rtol = 1e-3 atol = 1e-9
        @test outd.den_H     ≈ test_out[i,10]  rtol = 1e-3 atol = 1e-9
        @test outd.den_N     ≈ test_out[i,11]  rtol = 1e-3 atol = 1e-9
        @test outd.den_aO    ≈ test_out[i,12]  rtol = 1e-3 atol = 1e-9

        @test outd.den_Total ≈ test_out[i, 5] + 1.66e-24*16outd.den_aO rtol = 1e-3
    end

    ## Scenario 02
    ## ===========

    # Test outputs.
    test_out = [
       0.0  0.000E+00  1.919E+19  5.149E+18 1.181E-03  296.0  1027  1.288E+14 2.296E+17  0.000E+00  0.000E+00  0.000E+00   89.0   89.4    4.9    3.0 4.0    4.0    6.0    9.9   10.6;
     100.0  6.009E+11  9.392E+12  2.213E+12 5.764E-10  165.2  1027  1.339E+08 9.580E+10  2.654E+07  2.737E+05  0.000E+00   89.0   89.4    4.9    3.0 4.0    4.0    6.0    9.9   10.6;
     200.0  3.882E+09  1.978E+09  1.341E+08 2.028E-13  687.5   706  2.048E+07 9.719E+05  5.172E+05  1.632E+07  1.399E-09   89.0   89.4    4.9    3.0 4.0    4.0    6.0    9.9   10.6;
     300.0  3.140E+08  2.481E+07  9.340E+05 9.666E-15  705.5   706  1.082E+07 1.879E+03  3.860E+05  2.205E+06  3.876E+00   89.0   89.4    4.9    3.0 4.0    4.0    6.0    9.9   10.6;
     400.0  2.855E+07  3.737E+05  7.743E+03 8.221E-16  706.0   706  5.938E+06 4.688E+00  3.317E+05  2.633E+05  9.738E+02   89.0   89.4    4.9    3.0 4.0    4.0    6.0    9.9   10.6;
     500.0  2.787E+06  6.372E+03  7.382E+01 9.764E-17  706.0   706  3.319E+06 1.397E-02  2.868E+05  3.424E+04  3.186E+03   89.0   89.4    4.9    3.0 4.0    4.0    6.0    9.9   10.6;
     600.0  2.910E+05  1.222E+02  8.047E-01 2.079E-17  706.0   706  1.887E+06 4.919E-05  2.490E+05  4.742E+03  3.281E+03   89.0   89.4    4.9    3.0 4.0    4.0    6.0    9.9   10.6;
     700.0  3.240E+04  2.622E+00  9.973E-03 8.474E-18  706.0   706  1.090E+06 2.034E-07  2.171E+05  6.946E+02  2.498E+03   89.0   89.4    4.9    3.0 4.0    4.0    6.0    9.9   10.6;
     800.0  3.835E+03  6.264E-02  1.398E-04 4.665E-18  706.0   706  6.393E+05 9.809E-10  1.900E+05  1.074E+02  1.768E+03   89.0   89.4    4.9    3.0 4.0    4.0    6.0    9.9   10.6;
     900.0  4.816E+02  1.659E-03  2.204E-06 2.817E-18  706.0   706  3.806E+05 5.481E-12  1.669E+05  1.747E+01  1.236E+03   89.0   89.4    4.9    3.0 4.0    4.0    6.0    9.9   10.6;
    1000.0  6.400E+01  4.853E-05  3.891E-08 1.772E-18  706.0   706  2.298E+05 3.528E-14  1.471E+05  2.988E+00  8.675E+02   89.0   89.4    4.9    3.0 4.0    4.0    6.0    9.9   10.6;
    ]

    # Constant input parameters.
    year   = 2016
    month  = 6
    day    = 1
    hour   = 11
    minute = 00
    second = 00
    g_lat  = -23*pi/180
    g_long = -45*pi/180
    JD     = date_to_jd(year, month, day, hour, minute, second)

    # TODO: The tolerances here are much higher than in the previous test with
    # the daily AP only. This must be analyzed. However, it was observed that in
    # the NRLMSISE-00 version, the last component of the AP array is the
    # "Average of eight 3 hour AP indices from 36 to 57 hours prior to current
    # time." Whereas the label from the online version is "ap_33_59_hr_prior",
    # which seems to be different. On the other hand, the MSISE90 describe the
    # last vector of AP array as "Average of eight 3 hr AP indicies [sic] from
    # 36 to 59 hrs prior to current time." Hence, it seems that the online
    # version is using the old algorithm related to the AP array. We must change
    # the algorithm in our implementation to the old one and see if the result
    # is more similar to the online version.
    #
    # The information can be obtained at:
    #
    #   https://ccmc.gsfc.nasa.gov/modelweb/
    #   https://ccmc.gsfc.nasa.gov/pub/modelweb/atmospheric/msis/msise90/
    #

    for i = 1:size(test_out,1)
        f107  = test_out[i, 13]
        f107A = test_out[i, 14]
        ap_a  = test_out[i, 15:21]

        # Run the NRLMSISE-00 model wih the input parameters.
        out = nrlmsise00(JD,
                         test_out[i,1]*1000,
                         g_lat,
                         g_long,
                         f107A,
                         f107,
                         ap_a;
                         output_si = false,
                         dversion = false)

        @test out.den_O     ≈ test_out[i, 2]  rtol = 5e-2 atol = 1e-9
        @test out.den_N2    ≈ test_out[i, 3]  rtol = 5e-2 atol = 1e-9
        @test out.den_O2    ≈ test_out[i, 4]  rtol = 5e-2 atol = 1e-9
        @test out.T_alt     ≈ test_out[i, 6]  rtol = 1e-1 atol = 1e-9
        @test out.T_exo     ≈ test_out[i, 7]  rtol = 1e-1 atol = 1e-9
        @test out.den_He    ≈ test_out[i, 8]  rtol = 5e-2 atol = 1e-9
        @test out.den_Ar    ≈ test_out[i, 9]  rtol = 5e-2 atol = 1e-9
        @test out.den_H     ≈ test_out[i,10]  rtol = 5e-2 atol = 1e-9
        @test out.den_N     ≈ test_out[i,11]  rtol = 5e-2 atol = 1e-9
        @test out.den_aO    ≈ test_out[i,12]  rtol = 9e-2 atol = 1e-9
        @test out.den_Total ≈ test_out[i, 5]  rtol = 5e-2

        # Test the version that calls `gtd7d` instead of `gtd7`.
        outd = nrlmsise00(JD,
                          test_out[i,1]*1000,
                          g_lat,
                          g_long,
                          f107A,
                          f107,
                          ap_a;
                          output_si = false,
                          dversion = true)

        @test outd.den_O     ≈ test_out[i, 2]  rtol = 5e-2 atol = 1e-9
        @test outd.den_N2    ≈ test_out[i, 3]  rtol = 5e-2 atol = 1e-9
        @test outd.den_O2    ≈ test_out[i, 4]  rtol = 5e-2 atol = 1e-9
        @test outd.T_alt     ≈ test_out[i, 6]  rtol = 1e-1 atol = 1e-9
        @test outd.T_exo     ≈ test_out[i, 7]  rtol = 1e-1 atol = 1e-9
        @test outd.den_He    ≈ test_out[i, 8]  rtol = 5e-2 atol = 1e-9
        @test outd.den_Ar    ≈ test_out[i, 9]  rtol = 5e-2 atol = 1e-9
        @test outd.den_H     ≈ test_out[i,10]  rtol = 5e-2 atol = 1e-9
        @test outd.den_N     ≈ test_out[i,11]  rtol = 5e-2 atol = 1e-9
        @test outd.den_aO    ≈ test_out[i,12]  rtol = 9e-2 atol = 1e-9

        @test outd.den_Total ≈ test_out[i, 5] + 1.66e-24*16outd.den_aO rtol = 5e-2
    end
end


text:
using FSimZoo
using Test

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


text:
using BayesianDataFusion
using Base.Test

X = sprand(50, 100, 0.1)
Xr = sparse_csr(X)

@test size(X) == size(Xr)
@test size(X,1) == size(Xr,1)
@test size(X,2) == size(Xr,2)
@test eltype(X) == eltype(Xr)

y1 = rand(100)
y2 = rand(50)

u1 = X * y1
u2 = X' * y2

v1 = Xr * y1
v2a = Ac_mul_B(Xr, y2)
v2b = At_mul_B(Xr, y2)

@test_approx_eq u1 v1
@test_approx_eq u2 v2a
@test_approx_eq u2 v2b

@test_approx_eq full(At_mul_B(X,X)) full(At_mul_B(Xr,Xr))

rows = [1, 2, 2, 4]
cols = [2, 1, 3, 3]
vals = [0.1, 0.2, 0.15, 0.3]
z    = rand(3)
w1   = sparse(rows, cols, vals) * z
w2   = sparse_csr(rows, cols, vals) * z
@test_approx_eq w1 w2