text:
# In this tutorial, we will learn
# - How to use `DrWatson.jl` to accelerate and reproduce our Gridap simulation workflows

# [DrWatson.jl](https://github.com/JuliaDynamics/DrWatson.jl) is a Julia package that helps managing a typical scientific workflow thorought all its phases, see a summary [here](https://juliadynamics.github.io/DrWatson.jl/stable/workflow/).

# All its functionalities can be accessed with (non-invasive) simple function calls.

# In order to illustrate how to benefit from `DrWatson.jl` in `Gridap.jl` simulations, we refactor here the convergence test from the `Code validation` tutorial.

# Instead of implementing a helper function to carry out the convergence test, we will generate them using `DrWatson.jl` functions.

# ## 1. Activate your project

# The first step is to activate our project using `quickactivate`. This does not only activate the project, it also sets the relative paths within the project, so you can safely use the functions `projectdir()` and its derivatives `datadir()`, `plotsdir()`, `srcdir()`, etc. Beware of this [warning](https://juliadynamics.github.io/DrWatson.jl/dev/project/#DrWatson.quickactivate), you must activate the project before using other packages.

using DrWatson
@quickactivate "Tutorials"

# Although this tutorial is already in a Project (and git repo), we could also start our scientific project from scratch with `DrWatson.jl`, using function `initialize_project`. This function initiates, on the working directory, (1) a git repo with a folder structure enriched for scientific workflows, e.g. folders `data`, `plots`, `papers`, etc., and (2)`Project.toml` and `Manifest.toml` files. More details [here](https://juliadynamics.github.io/DrWatson.jl/dev/workflow/#.-Setup-the-project-1).

# Once the project is activated, we ensure that all packages we use have the versions dictated by our activated project.

using Gridap
import Gridap: ∇

# ## 2. Prepare the simulations

# We consider the Poisson equation in the unit square $\Omega\doteq (0,1)^2$ as a model problem,
#
# ```math
# \left\lbrace
# \begin{aligned}
# -\Delta u = f  \ \text{in} \ \Omega\\
# u = g \ \text{on}\ \partial\Omega.\\
# \end{aligned}
# \right.
# ```
#
# We are going to perform a convergence test with the manufactured solution $u(x) = x_1^3 + x_2^3$.

# To this end, we want to solve our computational model for many combinations of mesh size and order of FE approximation (*parameters*) and extract the L2- and H1-norm errors (*output data*).

# We first group all parameters and parameter values in a single ditionary

params = Dict(
  "cells_per_axis" => [8,16,32,64],
  "fe_order" => [1,2]
)

# and then we use DrWatson's `dict_list` to expand all the parameters into a vector of dictionaries. Each dictionary contains the parameter-value combinations corresponding to a single simulation case.

dicts = dict_list(params)

# **Warning!** Be careful when combining parameters of different value type. You may end up with dictionaries that do not have a concrete type and experience a significant type-inference overhead when running the simulations.

# We wrap next in a function a run of our computational model for a single pair `(cells_per_axis,fe_order)`. The function returns the L2- and H1-error norms.

# We define the manufactured function, as usual

p = 3
u(x) = x[1]^p+x[2]^p
∇u(x) = VectorValue(p*x[1]^(p-1),p*x[2]^(p-1))
f(x) = -p*(p-1)*(x[1]^(p-2)+x[2]^(p-2))
∇(::typeof(u)) = ∇u

# And the function that runs a single case of our parametric space reads

function run(n::Int,k::Int)

  domain = (0,1,0,1)
  partition = (n,n)
  model = CartesianDiscreteModel(domain,partition)

  reffe = ReferenceFE(lagrangian,Float64,k)
  V0 = TestFESpace(model,reffe,conformity=:H1,dirichlet_tags="boundary")
  U = TrialFESpace(V0,u)

  degree = 2*p
  Ω = Triangulation(model)
  dΩ = Measure(Ω,degree)

  a(u,v) = ∫( ∇(u)⊙∇(v) ) * dΩ
  b(v) = ∫( v*f ) * dΩ

  op = AffineFEOperator(a,b,U,V0)

  uh = solve(op)

  e = u - uh

  el2 = sqrt(sum( ∫( e*e )*dΩ ))
  eh1 = sqrt(sum( ∫( e*e + ∇(e)⋅∇(e) )*dΩ ))

  (el2, eh1)

end

# In order to communicate with `DrWatson.jl` helper functions, we need to add an extra layer on top of `run`, such that the input and output are dictionaries.

# Note the use of functions [@unpack](https://juliadynamics.github.io/DrWatson.jl/dev/name/#UnPack.@unpack) and [@dict](https://juliadynamics.github.io/DrWatson.jl/dev/name/#DrWatson.@dict) to decompose and compose the dictionaries. You can check in `DrWatson.jl`'s documentation further functions to manipulate dictionaries.

function run(case::Dict)
  @unpack cells_per_axis, fe_order = case
  el2, eh1 = run(cells_per_axis,fe_order)
  h = 1.0/cells_per_axis
  results = @strdict el2 eh1 h
  merge(case,results)
end

# ## 3. Run and save

# While running the simulations, we need to save the results. `DrWatson.jl` frees you from the burden of generating the filenames for each case. For this purpose, it provides the functions [savename](https://juliadynamics.github.io/DrWatson.jl/stable/name/#DrWatson.savename), [@tagsave](https://juliadynamics.github.io/DrWatson.jl/stable/save/#DrWatson.@tagsave) or [produce_or_load](https://juliadynamics.github.io/DrWatson.jl/stable/save/#DrWatson.produce_or_load), among others.

# Among them, we recommend using [produce_or_load](https://juliadynamics.github.io/DrWatson.jl/stable/save/#DrWatson.produce_or_load). The special feature of this function is that it checks whether the file containing the output data of the case already exists. If that happens, then the function loads the file, instead of running the case. In this way, we avoid repeating simulations that have already been run.

# Thus, in order to run all simulation cases, it suffices to map all cases in `dicts` to the `produce_or_load` function:

function run_or_load(case::Dict)
  produce_or_load(
    projectdir("assets","validation_DrWatson"),
    case,
    run,
    prefix="res",
    tag=true,
    verbose=true
  )
  return true
end

map(run_or_load,dicts)

# Note that the results of each case are stored in a binary database file in the `projectdir("assets","validation_DrWatson")` folder. Each result file stores the output dictionary that returns from `run(case)`.

# We also observe that we set `tag=true` in `produce_or_load`. This option is *key to preserve reproducibility*. It adds to the output dictionary the field `:gitcommit`, thus allowing us to trace the status of the code, at which we obtained those results. Furthermore, if the git repo is dirty, one more field `:gitpatch` is added, storing the difference string.

# In some situations, you will prefer to repeat all simulations and track their evolution as you change the code. To this end, check out [safesave](https://juliadynamics.github.io/DrWatson.jl/dev/save/#DrWatson.safesave).

# ## 4. Listing the simulations

# Results stored by `DrWatson.jl` in databases are handled with the [DataFrames.jl](https://dataframes.juliadata.org/stable/) package, a powerful Julia package to manipulate tabular data.

using DataFrames

# To collect all simulation results, it suffices to use the `collect_results!` function from `DrWatson.jl` from the folder where the results are stored.

df = collect_results(projectdir("assets","validation_DrWatson"))

# We order next the database by (ascending) mesh size and we extract the arrays of mesh sizes and errors

sort!(df,:h)
hs = df[(df.fe_order .== 1),:h]
el2s1 = df[(df.fe_order .== 1),:el2]
eh1s1 = df[(df.fe_order .== 1),:eh1]
el2s2 = df[(df.fe_order .== 2),:el2]
eh1s2 = df[(df.fe_order .== 2),:eh1]

# ## 5. Generate the plot

# With the generated data, we do the classical convergence plot and interpret it in the same way as in the validation tutorial.

using Plots

plot(hs,[el2s1 eh1s1 el2s2 eh1s2],
    xaxis=:log, yaxis=:log,
    label=["L2 k=1" "H1 k=1" "L2 k=2" "H1 k=2"],
    shape=:auto,
    xlabel="h",ylabel="error norm")

#src savefig("conv.png")

#md # If you run the code in a notebook, you will see a figure like this one:
#md # ![](../assets/validation_DrWatson/conv.png)
#
# Congrats, another tutorial done!

# ---

# If you use DrWatson.jl in a scientific project that leads to a publication, please do not forget to cite the paper associated with it:
# ```
# @article{Datseris2020,
#  doi = {10.21105/joss.02673},
#  url = {https://doi.org/10.21105/joss.02673},
#  year = {2020},
#  publisher = {The Open Journal},
#  volume = {5},
#  number = {54},
#  pages = {2673},
#  author = {George Datseris and Jonas Isensee and Sebastian Pech and Tamás Gál},
#  title = {DrWatson: the perfect sidekick for your scientific inquiries},
#  journal = {Journal of Open Source Software}
# }
# ```


text:
# Compare line drawing performance of Matplotlib vs. GR
#
# These are the results on a MacBook Pro 2,6 GHz Intel Core i5::
#
# MPL:   10 fps  GR: 1400 fps  => speedup:  140.0

import PyPlot
import GR
ENV["GKSwstype"] = "gksqt"

function mpl()
    x = collect(0:0.01:2*pi)

    PyPlot.plot(x, sin.(x))
    t = time_ns()
    for i = 1:100
        PyPlot.cla()
        PyPlot.plot(x .+ i / 10.0, sin.(x .+ i / 10.0))
        PyPlot.pause(0.0001)
    end

    fps = round(Int64, 100 / (1e-9 * (time_ns() - t)))
end

function gr()
    x = collect(0:0.01:2*pi)

    GR.plot(x, sin.(x))
    t = time_ns()
    for i = 1:100
        GR.plot(x .+ i / 10.0, sin.(x .+ i / 10.0))
    end

    fps = round(Int64, 100 / (1e-9 * (time_ns() - t)))
end

using Printf

for i in 1:10
    fps_mpl = mpl()
    fps_gr = gr()
    speedup = float(fps_gr) / fps_mpl

    @printf("MPL: %4d fps  GR: %4d fps  => speedup: %6.1f\r",
            fps_mpl, fps_gr, speedup)
end

@printf("\n")

ext: jl


text:
abstract type SubstitutionModel end

const SM = SubstitutionModel

"""
`NucleicAcidSubstitutionModel` is an abstract type that contains all models
describing a substitution process impacting biological sequences of `DNA` or
`RNA` with continous time Markov models.
"""
abstract type NucleicAcidSubstitutionModel <: SubstitutionModel end

const NASM = NucleicAcidSubstitutionModel

const Qmatrix = SMatrix{4, 4, Float64}

const Pmatrix = SMatrix{4, 4, Float64}

function Base.convert(::Type{T}, θ::F...; safe::Bool=true) where {T <: NASM, F <: Float64}
  return T(θ..., safe=safe)
end

function Base.convert(::Type{T}, θ_vec::A; safe::Bool=true) where {T <: NASM, A <: AbstractArray}
  return T(θ_vec, safe=safe)
end

function Base.convert(::Type{T}, θ_vec::A, π_vec::A; safe::Bool=true) where {T <: NASM, A <: AbstractArray}
  return T(θ_vec, π_vec, safe=safe)
end


text:
using DiffEqSensitivity, OrdinaryDiffEq, ForwardDiff, Calculus
using Test
function fb(du,u,p,t)
  du[1] = dx = p[1]*u[1] - p[2]*u[1]*u[2]
  du[2] = dy = -t*p[3]*u[2] + t*u[1]*u[2]
end
function jac(J,u,p,t)
  (x, y, a, b, c) = (u[1], u[2], p[1], p[2], p[3])
  J[1,1] = a + y * b * -1
  J[2,1] = t * y
  J[1,2] = b * x * -1
  J[2,2] = t * c * -1 + t * x
end
function paramjac(pJ,u,p,t)
  (x, y, a, b, c) = (u[1], u[2], p[1], p[2], p[3])
  pJ[1,1] = x
  pJ[2,1] = 0.0
  pJ[1,2] = - x * y
  pJ[2,2] = 0.0
  pJ[1,3] = 0.0
  pJ[2,3] = - t * y
end

f = ODEFunction(fb,jac=jac,paramjac=paramjac)
p = [1.5,1.0,3.0]
prob = ODEForwardSensitivityProblem(f,[1.0;1.0],(0.0,10.0),p)
probInpl = ODEForwardSensitivityProblem(fb,[1.0;1.0],(0.0,10.0),p)
probnoad = ODEForwardSensitivityProblem(fb,[1.0;1.0],(0.0,10.0),p,
                                        ForwardSensitivity(autodiff=false))
probnoadjacvec = ODEForwardSensitivityProblem(fb,[1.0;1.0],(0.0,10.0),p,
                                        ForwardSensitivity(autodiff=false,autojacvec=true))
probnoad2 = ODEForwardSensitivityProblem(f,[1.0;1.0],(0.0,10.0),p,
                                        ForwardSensitivity(autodiff=false))
probvecmat = ODEForwardSensitivityProblem(fb,[1.0;1.0],(0.0,10.0),p,
                                          ForwardSensitivity(autojacvec=false,autojacmat=true))
sol = solve(prob,Tsit5(),abstol=1e-14,reltol=1e-14)
@test_broken solve(probInpl,KenCarp4(),abstol=1e-14,reltol=1e-14).retcode == :Success
solInpl = solve(probInpl,KenCarp4(autodiff=false),abstol=1e-14,reltol=1e-14)
solInpl2 = solve(probInpl,Rodas4(autodiff=false),abstol=1e-10,reltol=1e-10)
solnoad = solve(probnoad,KenCarp4(autodiff=false),abstol=1e-14,reltol=1e-14)
solnoadjacvec = solve(probnoadjacvec,KenCarp4(autodiff=false),abstol=1e-14,reltol=1e-14)
solnoad2 = solve(probnoad,KenCarp4(autodiff=false),abstol=1e-14,reltol=1e-14)
solvecmat = solve(probvecmat,Tsit5(),abstol=1e-14,reltol=1e-14)

x = sol[1:sol.prob.f.numindvar,:]

@test sol(5.0) ≈ solnoad(5.0)
@test sol(5.0) ≈ solnoad2(5.0)
@test sol(5.0) ≈ solnoadjacvec(5.0) atol=1e-6 rtol=1e-6
@test sol(5.0) ≈ solInpl(5.0)
@test isapprox(solInpl(5.0), solInpl2(5.0),rtol=1e-5)
@test sol(5.0) ≈ solvecmat(5.0)

# Get the sensitivities

da = sol[sol.prob.f.numindvar+1:sol.prob.f.numindvar*2,:]
db = sol[sol.prob.f.numindvar*2+1:sol.prob.f.numindvar*3,:]
dc = sol[sol.prob.f.numindvar*3+1:sol.prob.f.numindvar*4,:]

sense_res1 = [da[:,end] db[:,end] dc[:,end]]

prob = ODEForwardSensitivityProblem(f.f,[1.0;1.0],(0.0,10.0),p,
                                    ForwardSensitivity(autojacvec=true))
sol = solve(prob,Tsit5(),abstol=1e-14,reltol=1e-14,saveat=0.01)
x = sol[1:sol.prob.f.numindvar,:]

# Get the sensitivities

res = sol[1:sol.prob.f.numindvar,:]
da = sol[sol.prob.f.numindvar+1:sol.prob.f.numindvar*2,:]
db = sol[sol.prob.f.numindvar*2+1:sol.prob.f.numindvar*3,:]
dc = sol[sol.prob.f.numindvar*3+1:sol.prob.f.numindvar*4,:]

sense_res2 = [da[:,end] db[:,end] dc[:,end]]

function test_f(p)
  prob = ODEProblem(f,eltype(p).([1.0,1.0]),(0.0,10.0),p)
  solve(prob,Tsit5(),abstol=1e-14,reltol=1e-14,save_everystep=false)[end]
end

p = [1.5,1.0,3.0]
fd_res = ForwardDiff.jacobian(test_f,p)
calc_res = Calculus.finite_difference_jacobian(test_f,p)

@test sense_res1 ≈ sense_res2 ≈ fd_res
@test sense_res1 ≈ sense_res2 ≈ calc_res

## Check extraction

xall, dpall = extract_local_sensitivities(sol)
@test xall == res
@test dpall[1] == da

_,dpall_matrix = extract_local_sensitivities(sol,Val(true))
@test mapreduce(x->x[:, 2], hcat, dpall) == dpall_matrix[2]

x, dp = extract_local_sensitivities(sol,length(sol.t))
sense_res2 = reduce(hcat,dp)
@test sense_res1 == sense_res2

@test extract_local_sensitivities(sol,sol.t[3]) == extract_local_sensitivities(sol,3)

tmp = similar(sol[1])
@test extract_local_sensitivities(tmp,sol,sol.t[3]) == extract_local_sensitivities(sol,3)

# asmatrix=true
@test extract_local_sensitivities(sol, length(sol), true) == (x, sense_res2)
@test extract_local_sensitivities(sol, sol.t[end], true) == (x, sense_res2)
@test extract_local_sensitivities(tmp, sol, sol.t[end], true) == (x, sense_res2)

# Return type inferred
@inferred extract_local_sensitivities(sol, 1)
@inferred extract_local_sensitivities(sol, 1, Val(true))
@inferred extract_local_sensitivities(sol, sol.t[3])
@inferred extract_local_sensitivities(sol, sol.t[3], Val(true))
@inferred extract_local_sensitivities(tmp, sol, sol.t[3])
@inferred extract_local_sensitivities(tmp, sol, sol.t[3], Val(true))

### ForwardDiff version

prob = ODEForwardSensitivityProblem(f.f,[1.0;1.0],(0.0,10.0),p,
                                    ForwardDiffSensitivity())
sol = solve(prob,Tsit5(),abstol=1e-14,reltol=1e-14,saveat=0.01)

xall, dpall = extract_local_sensitivities(sol)
@test xall ≈ res
@test dpall[1] ≈ da atol=1e-9

_,dpall_matrix = extract_local_sensitivities(sol,Val(true))
@test mapreduce(x->x[:, 2], hcat, dpall) == dpall_matrix[2]

x, dp = extract_local_sensitivities(sol,length(sol.t))
sense_res2 = reduce(hcat,dp)
@test fd_res == sense_res2

@test extract_local_sensitivities(sol,sol.t[3]) == extract_local_sensitivities(sol,3)

tmp = similar(sol[1])
@test extract_local_sensitivities(tmp,sol,sol.t[3]) == extract_local_sensitivities(sol,3)

# asmatrix=true
@test extract_local_sensitivities(sol, length(sol), true) == (x, sense_res2)
@test extract_local_sensitivities(sol, sol.t[end], true) == (x, sense_res2)
@test extract_local_sensitivities(tmp, sol, sol.t[end], true) == (x, sense_res2)

# Return type inferred
@inferred extract_local_sensitivities(sol, 1)
@inferred extract_local_sensitivities(sol, 1, Val(true))
@inferred extract_local_sensitivities(sol, sol.t[3])
@inferred extract_local_sensitivities(sol, sol.t[3], Val(true))
@inferred extract_local_sensitivities(tmp, sol, sol.t[3])
@inferred extract_local_sensitivities(tmp, sol, sol.t[3], Val(true))

# Test mass matrix
function rober_MM(du, u, p, t)
    y₁, y₂, y₃ = u
    k₁, k₂, k₃ = p
    du[1] = -k₁ * y₁ + k₃ * y₂ * y₃
    du[2] = k₁ * y₁ - k₂ * y₂^2 - k₃ * y₂ * y₃
    du[3] = y₁ + y₂ + y₃ - 1
    nothing
end
function rober_no_MM(du, u, p, t)
    y₁, y₂, y₃ = u
    k₁, k₂, k₃ = p
    du[1] = -k₁ * y₁ + k₃ * y₂ * y₃
    du[2] = k₁ * y₁ - k₂ * y₂^2 - k₃ * y₂ * y₃
    du[3] = k₂*y₂^2
    nothing
end

M = [1.0 0 0; 0 1.0 0; 0 0 0]
p = [0.04, 3e7, 1e4]
u0 = [1.0, 0.0, 0.0]
tspan = (0.0, 12.0)

f_MM= ODEFunction(rober_MM, mass_matrix = M)
f_no_MM= ODEFunction(rober_no_MM)

prob_MM_ForwardSensitivity = ODEForwardSensitivityProblem(f_MM, u0, tspan, p, ForwardSensitivity())
sol_MM_ForwardSensitivity = solve(prob_MM_ForwardSensitivity , Rodas4(autodiff = false), reltol = 1e-14, abstol = 1e-14)

prob_MM_ForwardDiffSensitivity = ODEForwardSensitivityProblem(f_MM, u0, tspan, p, ForwardDiffSensitivity())
sol_MM_ForwardDiffSensitivity = solve(prob_MM_ForwardDiffSensitivity, Rodas4(autodiff = false), reltol = 1e-14, abstol = 1e-14)

prob_no_MM = ODEForwardSensitivityProblem(f_no_MM, u0, tspan, p, ForwardSensitivity())
sol_no_MM= solve(prob_no_MM, Rodas4(autodiff = false), reltol = 1e-14, abstol = 1e-14)

sen_MM_ForwardSensitivity = extract_local_sensitivities(sol_MM_ForwardSensitivity,10.0,true)
sen_MM_ForwardDiffSensitivity = extract_local_sensitivities(sol_MM_ForwardDiffSensitivity,10.0,true)
sen_no_MM = extract_local_sensitivities(sol_no_MM,10.0,true)

@test sen_MM_ForwardSensitivity[2] ≈ sen_MM_ForwardDiffSensitivity[2] atol=1e-10 rtol=1e-10
@test sen_MM_ForwardSensitivity[2] ≈ sen_no_MM[2] atol=1e-10 rtol=1e-10
# Test Float32

function f32(du,u,p,t)
  du[1] = dx = p[1]*u[1] - p[2]*u[1]*u[2]
  du[2] = dy = -p[3]*u[2] + u[1]*u[2]
end
p = [1.5f0,1.0f0,3.0f0]
prob = ODEForwardSensitivityProblem(f32,[1.0f0;1.0f0],(0.0f0,10.0f0),p)
sol = solve(prob,Tsit5())


text:
function create_model(setup)
    coeffs = get_coefficients(setup)
    mesh, = create_geometry(setup)

    Model(; setup.name, mesh, coeffs...)
end


text:
"""
    uniform_bsplines_eval_basis( spline_degree, normalized_offset, bspl )

# UNIFORM B-SPLINE FUNCTIONS

## Evaluate all non vanishing uniform B-Splines in unit cell.

Returns an array with the values of the b-splines of the
requested degree, evaluated at a given cell offset. The cell size is
normalized between 0 and 1, thus the offset given must be a number
between 0 and 1.

Output:

```math
bspl(1:d+1)= B_d(-(d+1)/2+d+x),...,B_d(-(d+1)/2+x)
```

with ``d``=`spline_degree` and ``x``=`normalized_offset`
where ``B_d=B_{d-1}*B_0`` and ``B_0=1_[-1/2,1/2]`` and `*` is convolution
the following FORTRAN code can be used for comparison with
[deboor](http://pages.cs.wisc.edu/~deboor/)

```fortran
do i=-d,d+1
    t(i+d+1)=real(i,8)
end do
call bsplvb(t,d+1,1,normalized_offset,d+1,out)
```

We also have the property (from the symmetry of the B-spline)
```math
out[1:d+1]= B_d(-(d+1)/2+xx),...,B_d(-(d+1)/2+d+xx),...,
```
where ``xx=1-`` `normalized_offset`

"""
function uniform_bsplines_eval_basis(spline_degree::Int, normalized_offset::Float64)
    @assert spline_degree >= 0
    @assert normalized_offset >= 0.0
    @assert normalized_offset <= 1.0

    bspl = zeros(Float64, spline_degree + 1)

    bspl[1] = 1.0
    @inbounds for j in 1:spline_degree
        xx = -normalized_offset
        j_real = Float64(j)
        inv_j = 1.0 / j_real
        saved = 0.0
        for r in 0:(j - 1)
            xx = xx + 1.0
            temp = bspl[r + 1] * inv_j
            bspl[r + 1] = saved + xx * temp
            saved = (j_real - xx) * temp
        end
        bspl[j + 1] = saved
    end

    return bspl
end

function uniform_bsplines_eval_basis!(
    bspl::Vector{Float64}, spline_degree::Int, normalized_offset::Float64
)
    bspl[1] = 1.0
    @inbounds for j in 1:spline_degree
        xx = -normalized_offset
        j_real = Float64(j)
        inv_j = 1.0 / j_real
        saved = 0.0
        for r in 0:(j - 1)
            xx = xx + 1.0
            temp = bspl[r + 1] * inv_j
            bspl[r + 1] = saved + xx * temp
            saved = (j_real - xx) * temp
        end
        bspl[j + 1] = saved
    end
end

function uniform_bsplines_eval_basis!(
    bspl::Array{Float64,2},
    spline_degree::Int,
    normalized_offset_x::Float64,
    normalized_offset_y::Float64,
)
    bspl[1, 1] = 1.0
    bspl[1, 2] = 1.0
    @inbounds for j in 1:spline_degree
        xx = -normalized_offset_x
        yy = -normalized_offset_y
        j_real = Float64(j)
        inv_j = 1.0 / j_real
        saved_x = 0.0
        saved_y = 0.0
        for r in 0:(j - 1)
            xx = xx + 1
            yy = yy + 1
            temp_x = bspl[r + 1, 1] * inv_j
            temp_y = bspl[r + 1, 2] * inv_j
            bspl[r + 1, 1] = saved_x + xx * temp_x
            bspl[r + 1, 2] = saved_y + yy * temp_y
            saved_x = (j_real - xx) * temp_x
            saved_y = (j_real - yy) * temp_y
        end
        bspl[j + 1, 1] = saved_x
        bspl[j + 1, 2] = saved_y
    end
end

export eval_uniform_periodic_spline_curve

"""
    eval_uniform_periodic_spline_curve( degree, scoef )

Evaluate uniform periodic spline curve defined by coefficients scoef at
knots (which are the grid points)

"""
function eval_uniform_periodic_spline_curve(degree::Int, scoef::Vector{Float64})

    # get bspline values at knots
    bspl = uniform_bsplines_eval_basis(degree, 0.0)
    n = length(scoef)
    sval = similar(scoef)
    for i in 1:n
        val = 0.0
        for j in 1:degree
            imj = mod1(i - j + n, n)
            val = val + bspl[j] * scoef[imj]
        end
        sval[i] = val
    end

    return sval
end


text:
using GuaranteedEllipseFitting
using Test, StaticArrays, LinearAlgebra

@testset "GuaranteedEllipseFitting.jl" begin
    include("data_points.jl")
    include("test_sampling.jl")
    include("test_observations.jl")
    include("test_conversion.jl")
    include("test_normalise_data_context.jl")
    include("test_direct_ellipse_fit.jl")
    include("guaranteed_ellipse_fit.jl")
    include("test_covariance.jl")
end


text:
module Textual

using ..FastAI
using ..FastAI:
    # blocks
    Block, WrapperBlock, AbstractBlock, OneHotTensor, OneHotTensorMulti, Label,
    LabelMulti, wrapped, Continuous, getencodings, getblocks, encodetarget, encodeinput,
    # encodings
    Encoding, StatefulEncoding, OneHot,
    # visualization
    ShowText,
    # other
    Context, Training, Validation, FASTAI_METHOD_REGISTRY, registerlearningtask!

import Requires: @require

using InlineTest
using Random

include("recipes.jl")
include("blocks/text.jl")
include("transform.jl")

function __init__()
    _registerrecipes()
end

export Paragraph
end


text:
"""
# PersistenceDiagrams.jl

Types and functions for working with persistence diagrams.

See [the docs](https://mtsch.github.io/PersistenceDiagrams.jl/dev/) for documentation.
"""
module PersistenceDiagrams

export PersistenceDiagram,
    PersistenceInterval,
    birth,
    death,
    persistence,
    midlife,
    representative,
    birth_simplex,
    death_simplex,
    dim,
    threshold,
    Bottleneck,
    Wasserstein,
    weight,
    matching,
    PersistenceImage,
    PersistenceCurve,
    BettiCurve,
    Life,
    Midlife,
    LifeEntropy,
    MidlifeEntropy,
    PDThresholding,
    Landscape,
    Landscapes,
    Silhuette,
    barcode,
    PersistenceImageVectorizer,
    PersistenceCurveVectorizer,
    PersistenceLandscapeVectorizer

using Compat
using Hungarian
using PersistenceDiagramsBase
using RecipesBase
using Statistics
using Tables
import ScientificTypes

include("matching.jl")

include("persistencecurves.jl")
include("persistenceimages.jl")

include("plotsrecipes.jl")

include("scitypes.jl")
include("mlj.jl")

end


text:
struct U0 <: Exception end
struct U1 <: Exception end

function foo()
    for i in 1:2
        try
            bar()
        catch err
            if isa(err, U0) println("catched U0")
            else rethrow(err) end
        end
    end
end

function bar()
    baz()
end

function baz()
    if isdefined(:_called) && _called
        throw(U1())
    else
        global _called = true
        throw(U0())
    end
end

foo()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module Solve

using JuLIP, QMMM2, LinearAlgebra
using JuLIPMaterials, StaticArrays
using JuLIP.MLIPs

function autocell12(atref, atnew)
   C = Matrix(cell(atref))
   x, y, _ =  xyz(atnew)
   @assert minimum(x) >= 0
   @assert minimum(y) >= 0
   C[1,1] = maximum(x)
   C[2,2] = maximum(y)
   return C
end

"""
`function cluster_sequence(atmax, RDOM, RBUF,

* `get_data(atmax, "xcore")` must return the position of the defect core

"""
function cluster_sequence(atmax, RDOM, RBUF,
                          dims = findall(.!(pbc(atmax))))
   @assert sort(RDOM) == RDOM
   X_max = positions(atmax)
   x0 = get_data(atmax, "xcore")::JVecF

   Iinmax = Vector{Vector{Int}}(undef, length(RDOM))
   Iinnext = Vector{Vector{Int}}(undef, length(RDOM))
   AT = Vector{Atoms{Float64}}(undef, length(RDOM))

   at_next = atmax
   Iinmax_next = 1:length(atmax)

   for n = length(RDOM):-1:1
      X_next = positions(at_next)
      r_next = [ norm(x[dims] - x0[dims]) for x in X_next ]
      Iinnext[n] =  findall( r_next .<= RDOM[n] + RBUF[n] )
      Iinmax[n] = Iinmax_next[Iinnext[n]]
      Iinmax_next = Iinmax[n]
      at_next = Atoms(:Si, X_next[Iinnext[n]])
      empty!(at_next.Z)
      append!(at_next.Z, fill(atomic_number(:Si), length(at_next)) )
      r = r_next[Iinnext[n]]
      set_free!(at_next, findall(r .<=  RDOM[n]))
      set_pbc!(at_next, pbc(atmax))
      set_cell!(at_next, autocell12(atmax, at_next))
      AT[n] = at_next
   end

   # test this construction
   for n = 1:length(RDOM)-1
      at = AT[n]
      X = positions(at)
      at_next = AT[n+1]
      X_next = positions(at_next)
      @assert X == X_max[Iinmax[n]]
      @assert X == X_next[Iinnext[n]]
   end
   X = positions(AT[end])
   @assert X ==  X_max[Iinmax[end]]

   return AT, Iinmax, Iinnext
end

function solve_sequence(atmax, AT, Iinmax, Iinnext,
                        calc = calculator(atmax),
                        preconcalc = calc)
   UU = Vector{Vector{JVecF}}(undef, length(AT))
   EE = Vector{Float64}(undef, length(AT))

   # solve the first one
   @show length(AT[1])
   set_calculator!(AT[1], calc)
   X0 = positions(AT[1])
   E0 = energy(AT[1])
   minimise!(AT[1]; verbose=2, precond=FF(AT[1]))
   UU[1] = positions(AT[1]) .- X0
   EE[1] = energy(AT[1]) - E0

   # now solve sequentially
   for n = 2:length(AT)
      @show length(AT[n])
      set_calculator!(AT[n], calc)
      E0 = energy(AT[n])
      X0 = positions(AT[n])
      X = copy(X0)
      X[Iinnext[n-1]] += UU[n-1]
      set_positions!(AT[n], X)
      minimise!(AT[n]; verbose=2, precond=FF(AT[n], preconcalc))
      UU[n] = positions(AT[n]) - X0
      EE[n] = energy(AT[n]) - E0
   end

   # and the last one - the reference calculation
   @show length(atmax)
   set_calculator!(atmax, calc)
   X0max = positions(atmax)
   E0max = energy(atmax)
   X = copy(X0max)
   X[Iinmax[end]] += UU[end]
   set_positions!(atmax, X)
   minimise!(atmax; verbose=2, precond=FF(atmax, preconcalc))
   Umax = positions(atmax) - X0max
   Emax = energy(atmax) - E0max

   return UU, EE, Umax, Emax
end

function errors(atmax, AT, Iinmax, UU, Umax; rcut=3.5)
   errs2 = zeros(length(AT))
   errsinf = zeros(length(AT))
   for n = 1:length(AT)
      at = AT[n]
      Uerr = copy(Umax)
      Uerr[Iinmax[n]] -= UU[n]
      ee = JuLIPMaterials.strains(Uerr, atmax, rcut = rcut)
      errs2[n] = norm(ee, 2)
      errsinf[n] = norm(ee, Inf)
   end
   return errs2, errsinf
end

end


text:
#####
##### Utility Functions
#####

# Construct a `shape` vector from `x`.
# Since Julia is column major while C++ is row major, we need to be careful about the order
# we pass things.
#
# In general, the C++ code will not do any reversing schenanigans.
# Thus, it's up to the Julia interfacing code to make sure the appropriate conversions
# between column major and row major are made.
#
# nGraph Ordered Types shadowed
_reversed(x) = [Int64(i) for i in reverse(x)]
shape(x) = _reversed(x)
strides(x) = _reversed(x)

# It is convenient to construct Shape, Strides etc. from Tuples, Numbers, Vectors etc.
# Here, is the machinery that does the dispatch.
_expand(N, x) = fill(Int(x), N)

# Numbers
shape(N, x::Number) = _expand(N, x)
strides(N, x::Number) = _expand(N, x)

# Tuples and Vectors
maybecollect(x::Vector) = x
maybecollect(x) = collect(x)

function ordered(N, x)
    # Determine how many repetitions we need
    repitions, remainder = divrem(length(x), N)
    if !iszero(remainder)
        error("The length of `x` must be divisible by $N")
    end

    return _reversed(repeat(maybecollect(x), repitions))
end
shape(N, x::Union{Tuple,Vector}) = ordered(N, x)
strides(N, x::Union{Tuple,Vector}) = ordered(N, x)

#####
##### Element Type Maps
#####

# This should be kept inline with `ngraph/type/element_type.hpp`
@enum NGElements::Int32 begin
    ng_undefined=0
    ng_dynamic
    ng_boolean
    ng_bf16
    ng_f16
    ng_f32
    ng_f64
    ng_i8
    ng_i16
    ng_i32
    ng_i64
    ng_u1
    ng_u8
    ng_u16
    ng_u32
    ng_u64
end

# This is kind of a gross way of mapping Julia types to ngraph types.
# TODO: Think of a better way of doing this.
const TYPEMAPS = Dict(
   ng_boolean => Bool,
   ng_f32     => Float32,
   ng_f64     => Float64,
   ng_i8      => Int8,
   ng_i16     => Int16,
   ng_i32     => Int32,
   ng_i64     => Int64,
   ng_u8      => UInt8,
   ng_u16     => UInt16,
   ng_u32     => UInt32,
   ng_u64     => UInt64,
)

const Element = Lib.CxxWrap.CxxWrapCore.ConstCxxPtr{Lib.Element}

# Mapping from Julia => nGraph
Element(::Type{T}) where {T} = error("No translation defined for $T")
for (S,T) in TYPEMAPS
    @eval Element(::Type{$T}) = Lib.ngraph_type($(Int32(S)))
end

back(x) = TYPEMAPS[NGElements(@ngraphcall Type_t(x[]))]
ngraph_convert(::Type{T}) where {T <: Number} = Element(T)[]

#####
##### Nodes
#####

const NodeCppType = Lib.CxxWrap.StdLib.SharedPtrAllocated{nGraph.Lib.Node}
Base.ndims(x::NodeCppType) = length(@ngraphcall get_output_shape(x, 0))
Base.eltype(x::NodeCppType) = back(@ngraphcall get_output_element_type(x, 0))

# Subtype AbstractArray to get all of the nice fallbacks.
struct Node{T,N} <: AbstractArray{T,N}
    obj::NodeCppType
end
ngraph_convert(x::Node) = x.obj

Base.show(io::IO, x::Node{T,N}) where {T,N} = println(io, "Node{$T,$N} - $(name(x))")
Base.display(x::Node) = show(stdout, x)

Node{T}(x::NodeCppType) where {T} = Node{T,ndims(x)}(x)
Node(x::NodeCppType) = Node{eltype(x),ndims(x)}(x)

Node(x::AbstractArray{T,N}) where {T,N} = Node{T,N}(x)
Node{T,N}(x::AbstractArray{T,N}) where {T,N} = parameter(T, size(x))

Node(x::T) where {T <: Number} = Node{T}(x)
Node{T}(x::T) where {T <: Number} = Node{T,0}(fill(x))

# Array style arguments
Base.ndims(x::Node) = ndims(x.obj)

function Base.size(x::Node)
    nd = ndims(x)
    dims = @ngraphcall get_output_shape(x, 0)
    # Reversing is done on the ngraph side.
    return Tuple(Int.(reverse(dims)))
end

Base.size(x::Node, i::Integer) = size(x)[i]
Base.length(x) = prod(size(x))
Base.eltype(x::Node{T}) where {T} = T
Base.IndexStyle(::Node) = Base.IndexLinear()

name(x::Node) = String(@ngraphcall get_name(x))
description(x::Node) = String(@ngraphcall description(x))

# So these can be used as keys in a Dict
Base.:(==)(x::T, y::T) where {T <: Node} = name(x) == name(y)
Base.hash(x::Node, h::UInt = UInt(0x209348)) = hash(name(x), h)

function Base.getindex(::Node, i::Int)
    error("Yeah, yeah, I know \"Node\" is an AbstractArray ... but please don't index into it.")
end

#####
##### NGFunction
#####

mutable struct NGFunction
    ptr::Lib.CxxWrap.StdLib.SharedPtrAllocated{nGraph.Lib.NGFunction}
end
ngraph_convert(x::NGFunction) = x.ptr

function NGFunction(parameters::Vector, results::Vector)
    ptr = make_function(parameters, results)
    return NGFunction(ptr)
end

function make_function(inputs::Vector, outputs::Vector)
    return @ngraphcall make_function(outputs, inputs)
end

name(f::NGFunction) = String(@ngraphcall get_name(f))
poolsize(f::NGFunction) = Int(@ngraphcall get_temporary_pool_size(f))

#####
##### Backend
#####

struct Backend{T}
    ptr::Lib.CxxWrap.StdLib.SharedPtrAllocated{nGraph.Lib.Backend}
end
ngraph_convert(x::Backend) = x.ptr

#Backend(str::String = "CPU") = Backend(Lib.create(str))
Backend{T}() where {T} = Backend{T}(@ngraphcall create(string(T)))
version(x::Backend) = String(Lib.get_version(ngraph_convert(x)))

#####
##### Tensor
#####

struct Tensor
    ptr::Lib.CxxWrap.StdLib.SharedPtrAllocated{nGraph.Lib.RuntimeTensor}
    backend::Backend

    # Keep track of the base array to avoid GC
    # TODO: Rethink the whole - "casting to pointer" thing.
    parent::AbstractArray

    function Tensor(backend::Backend, v::Array{T}) where {T}
        # This is kind of scary - we ... just have to make sure that the parent array
        # doesn't get moved (i.e. resized ... )
        vptr = Base.unsafe_convert(Ptr{Cvoid}, pointer(v))
        ptr = Lib.create_tensor(
            backend.ptr,
            Element(T)[],
            shape(size(v)),
            vptr
        )

        return new(ptr, backend, v)
    end
end
ngraph_convert(x::Tensor) = x.ptr

Base.eltype(x::Tensor) = back(Lib.get_element_type(ngraph_convert(x)))
Base.sizeof(x::Tensor) = convert(Int, Lib.get_size_in_bytes(ngraph_convert(x)))

Tensor(backend, x::Tensor) = x

# Capture scalars in a 0-dimensional array
function Tensor(backend, x::T) where {T}
    A = Array{T,0}(undef)
    A[] = x
    return Tensor(backend, A)
end

# If we're trying to create a tensor view from a node
# create an undefined array and return a view of that.
function Tensor(backend::Backend, v::Node{T,N}) where {T,N}
    A = Array{T}(undef, size(v))
    return Tensor(backend, A)
end

function Base.size(t::Tensor)
    shape = Lib.get_shape(ngraph_convert(t))
    return Tuple(reverse(shape))
end

Base.parent(t::Tensor) = t.parent

function Base.show(io::IO, TV::Tensor)
    printstyled(io, "Tensor View\n"; color = :yellow)
    println(io, parent(TV))
end

#####
##### Executable
#####

mutable struct Executable
    ptr::Lib.CxxWrap.StdLib.SharedPtrAllocated{nGraph.Lib.Executable}
    ngraph_function::NGFunction
    backend::Backend

    function Executable(ptr, ngraph_function::NGFunction, backend::Backend)
        ex = new(ptr, ngraph_function, backend)

        # On cleanup - remove the compiled function.
        finalizer(ex) do x
            Lib.remove_compiled_function(x.backend.ptr, x.ptr)
        end
        return ex
    end
end
ngraph_convert(ex::Executable) = ex.ptr

function (ex::Executable)(
        inputs::Vector{Tensor},
        outputs::Vector{Tensor}
    )
    return @ngraphcall call(ex, outputs, inputs)
end

function compile(backend::Backend, inputs::Vector, outputs::Vector; kw...)
    # Gather up the raw shared pointers
    fn = NGFunction(inputs, outputs)
    compile(backend::Backend, fn; kw...)
end

function compile(
        backend::Backend,
        ngraph_function::NGFunction;
        emit_timing::Bool = false,
    )

    @time pointer = Lib.compile(backend.ptr, ngraph_function.ptr, emit_timing)

    # Indicate that the compiler has been invoked.
    return Executable(pointer, ngraph_function, backend)
end

# Conversions for vectors of Nodes and Tensors
ngraph_convert(x::Vector{<:Node}) = Lib.CxxWrap.CxxRef.(ngraph_convert.(x))
ngraph_convert(x::Vector{<:Tensor}) = Lib.CxxWrap.CxxRef.(ngraph_convert.(x))

#####
##### Extract performance data
#####

# TODO: Fix This
# function get_performance(ex::Executable)
#     # Construct a PerfCounterTranslator
#     translator = Lib.PerfCounterTranslator(getpointer(ex))
#
#     # Create a dictionary of timing results. Iterate through the CounterTranslator to
#     # construct the dict.
#     times = Dict{String,Int}()
#     for i in 1:Lib._length(translator)
#         name, time = Lib._getindex(translator, i-1)
#         times[name] = convert(Int, time)
#     end
#     return times
# end


text:
using EcologicalNetworks
using EcologicalNetworksPlots
using Plots
using CSV

##################################
##################################
df = CSV.read("WoodEtal_Append1_v2.csv")

# Separate based on WebScale
scale = [df[df.WebScale .== u,:] for u in unique(df.WebScale)]

# select for edge interaction columns only
edgelist_A = scale[3][:, [:PredTSN, :PreyTSN]]

# remove the "san" from the dataframe itself
edgelist_A.PredTSN = replace.(edgelist_A.PredTSN, "san" .=> "")

### Do the same for PreyTSN
# remove the "san" from the dataframe itself
edgelist_A.PreyTSN = replace.(edgelist_A.PreyTSN, "san" .=> "")
edgelist_A

# convert "string" to "integer"
# edgelist_A.PredTSN = parse.(Int64, edgelist_A.PredTSN)
# edgelist_A.PreyTSN = parse.(Int64, edgelist_A.PreyTSN)

#############################################
##Change PredTSN to newID starting from 1
# Sort by ID
sort!(edgelist_A, :PredTSN)
# Get unique IDs
unique_ids = unique(edgelist_A.PredTSN)
#create new IDs
newid_pred = collect(1:length(unique_ids))
# Create Dict where id => newid
id_dict = Dict(unique_ids[i] => newid_pred[i] for i in 1:length(unique_ids))
# Add new ids in dataframe
edgelist_A.newPredID = [id_dict[id] for id in edgelist_A.PredTSN]

### Do same for prey
## Change PreyTSN to new ID starting from end of newPredID
# Sort by ID
sort!(edgelist_A, :PreyTSN)
# Get unique IDs
unique_ids_prey = unique(edgelist_A.PreyTSN)
# Create new IDs
newid_prey = collect(length(id_dict)+1:length(id_dict)+length(unique_ids_prey))
# Create Dict
id_dict2 = Dict(unique_ids_prey[i] => newid_prey[i] for i in 1:length(unique_ids_prey))
# Add new ids in dataframe
edgelist_A.newPreyID = [id_dict2[id] for id in edgelist_A.PreyTSN]

A = edgelist_A[:, [:newPredID, :newPreyID]]

#sort
sort!(A, :newPredID)

#get unique values
unique!(A)

#############################################
##Create Adjacency matrix
#create matrix of zeros with n = colPred, m = colPrey
M = zeros(Int16, length(unique(A.newPredID)), length(unique(A.newPreyID)))

# Fill adjacency matrix with ones where there are interactions
for i in 1:length(A.newPredID)
    M[A[i, 1], A[i, 2] - length(unique(A.newPredID))] = 1
end
M

# Convert ones and zeros in boolean
M_bool = convert(Array{Bool}, M_wcol .== 1)

# Convert matrix in Bipartite network
B = BipartiteNetwork(M_bool)

########################################
#Plot

I = initial(RandomInitialLayout, B)
for step in 1:4000
  position!(ForceDirectedLayout(true, false, 2.5), I, B)
end
plot(I, B, title="Réseau d'intéractions de l'Archipel Sanak", lc=:lightblue)
scatter!(I, B, mc=:lightgrey)
png("figures/web_plot.png")

#species with a degree larger than 3 in red:
I3 = initial(RandomInitialLayout, B)
for step in 1:4000
  position!(ForceDirectedLayout(2.5), I3, B)
end
plot(I3, B, aspectratio=1)
scatter!(I3, B)
N = convert(AbstractUnipartiteNetwork, convert(BinaryNetwork, B))
core3 = collect(keys(filter(p -> p.second ≥ 3, degree(N))))
plot!(I3, N[core3], lc=:red)
scatter!(I3, N[core3], mc=:red)

max_stars_repo_name: kiristern/reseaux

max_issues_repo_name: kiristern/reseaux

max_forks_repo_name: kiristern/reseaux


text:
function Transform(arg0::AbsoluteDate, arg1::AngularCoordinates)
    return Transform((AbsoluteDate, AngularCoordinates), arg0, arg1)
end

function Transform(arg0::AbsoluteDate, arg1::PVCoordinates)
    return Transform((AbsoluteDate, PVCoordinates), arg0, arg1)
end

function Transform(arg0::AbsoluteDate, arg1::Rotation)
    return Transform((AbsoluteDate, Rotation), arg0, arg1)
end

function Transform(arg0::AbsoluteDate, arg1::Rotation, arg2::Vector3D)
    return Transform((AbsoluteDate, Rotation, Vector3D), arg0, arg1, arg2)
end

function Transform(arg0::AbsoluteDate, arg1::Rotation, arg2::Vector3D, arg3::Vector3D)
    return Transform((AbsoluteDate, Rotation, Vector3D, Vector3D), arg0, arg1, arg2, arg3)
end

function Transform(arg0::AbsoluteDate, arg1::Transform, arg2::Transform)
    return Transform((AbsoluteDate, Transform, Transform), arg0, arg1, arg2)
end

function Transform(arg0::AbsoluteDate, arg1::Vector3D)
    return Transform((AbsoluteDate, Vector3D), arg0, arg1)
end

function Transform(arg0::AbsoluteDate, arg1::Vector3D, arg2::Vector3D)
    return Transform((AbsoluteDate, Vector3D, Vector3D), arg0, arg1, arg2)
end

function Transform(arg0::AbsoluteDate, arg1::Vector3D, arg2::Vector3D, arg3::Vector3D)
    return Transform((AbsoluteDate, Vector3D, Vector3D, Vector3D), arg0, arg1, arg2, arg3)
end

function freeze(obj::Transform)
    return jcall(obj, "freeze", Transform, ())
end

function get_acceleration(obj::Transform)
    return jcall(obj, "getAcceleration", Vector3D, ())
end

function get_angular(obj::Transform)
    return jcall(obj, "getAngular", AngularCoordinates, ())
end

function get_cartesian(obj::Transform)
    return jcall(obj, "getCartesian", PVCoordinates, ())
end

function get_date(obj::Transform)
    return jcall(obj, "getDate", AbsoluteDate, ())
end

function get_inverse(obj::Transform)
    return jcall(obj, "getInverse", Transform, ())
end

function get_jacobian(obj::Transform, arg0::CartesianDerivativesFilter, arg1::Vector{Vector{jdouble}})
    return jcall(obj, "getJacobian", void, (CartesianDerivativesFilter, Vector{Vector{jdouble}}), arg0, arg1)
end

function get_rotation(obj::Transform)
    return jcall(obj, "getRotation", Rotation, ())
end

function get_rotation_acceleration(obj::Transform)
    return jcall(obj, "getRotationAcceleration", Vector3D, ())
end

function get_rotation_rate(obj::Transform)
    return jcall(obj, "getRotationRate", Vector3D, ())
end

function get_translation(obj::Transform)
    return jcall(obj, "getTranslation", Vector3D, ())
end

function get_velocity(obj::Transform)
    return jcall(obj, "getVelocity", Vector3D, ())
end

function interpolate(::Type{Transform}, arg0::AbsoluteDate, arg1::CartesianDerivativesFilter, arg2::AngularDerivativesFilter, arg3::Collection)
    return jcall(Transform, "interpolate", Transform, (AbsoluteDate, CartesianDerivativesFilter, AngularDerivativesFilter, Collection), arg0, arg1, arg2, arg3)
end

function interpolate(obj::TimeInterpolable, arg0::AbsoluteDate, arg1::Collection)
    return jcall(obj, "interpolate", TimeInterpolable, (AbsoluteDate, Collection), arg0, arg1)
end

function interpolate(obj::Transform, arg0::AbsoluteDate, arg1::Stream)
    return jcall(obj, "interpolate", Transform, (AbsoluteDate, Stream), arg0, arg1)
end

function shifted_by(obj::Transform, arg0::jdouble)
    return jcall(obj, "shiftedBy", Transform, (jdouble,), arg0)
end

function transform_line(obj::Transform, arg0::Line)
    return jcall(obj, "transformLine", Line, (Line,), arg0)
end

function transform_position(obj::Transform, arg0::FieldVector3D)
    return jcall(obj, "transformPosition", FieldVector3D, (FieldVector3D,), arg0)
end

function transform_position(obj::Transform, arg0::Vector3D)
    return jcall(obj, "transformPosition", Vector3D, (Vector3D,), arg0)
end

function transform_pv_coordinates(obj::Transform, arg0::FieldPVCoordinates)
    return jcall(obj, "transformPVCoordinates", FieldPVCoordinates, (FieldPVCoordinates,), arg0)
end

function transform_pv_coordinates(obj::Transform, arg0::PVCoordinates)
    return jcall(obj, "transformPVCoordinates", PVCoordinates, (PVCoordinates,), arg0)
end

function transform_pv_coordinates(obj::Transform, arg0::TimeStampedFieldPVCoordinates)
    return jcall(obj, "transformPVCoordinates", TimeStampedFieldPVCoordinates, (TimeStampedFieldPVCoordinates,), arg0)
end

function transform_pv_coordinates(obj::Transform, arg0::TimeStampedPVCoordinates)
    return jcall(obj, "transformPVCoordinates", TimeStampedPVCoordinates, (TimeStampedPVCoordinates,), arg0)
end

function transform_vector(obj::Transform, arg0::FieldVector3D)
    return jcall(obj, "transformVector", FieldVector3D, (FieldVector3D,), arg0)
end

function transform_vector(obj::Transform, arg0::Vector3D)
    return jcall(obj, "transformVector", Vector3D, (Vector3D,), arg0)
end


text:
export proxProj!, normProj

"""
    proxProj!(x::Vector{T}, λ::Float64; projFunc=x->x, kargs...)

applies the projection given by `projFunc`.
"""
function proxProj!(x::Vector{T}, λ::Float64; projFunc=x->x, kargs...) where T
  x[:] = projFunc(x)
end

"""
    normProj(x::Vector{T}, λ::Float64=0.0; projFunc=x->x, kargs...) where T

evaluate indicator function of set to be projected onto.
"""
function normProj(x::Vector{T}, λ::Float64=0.0; projFunc=x->x, kargs...) where T
  y = copy(x)
  y[:] = proxProj!(y,λ,projFunc=projFunc)
  if y != x
    return Inf
  end
  return 0.
end


text:
mutable struct DE <: AbstractParameters
    N::Int
    F::Float64
    CR::Float64
    CR_min::Float64
    CR_max::Float64
    F_min::Float64
    F_max::Float64
    strategy::Symbol
end

"""
    DE(;
        N  = 0,
        F  = 1.0,
        CR = 0.9,
        strategy = :rand1,
        information = Information(),
        options = Options()
    )

Parameters for Differential Evolution (DE) algorithm: step-size `F`,`CR` controlls the binomial
crossover, `N` is the population size. The parameter `trategy` is related to the variation
operator (`:rand1`, `:rand2`, `:best1`, `:best2`, `:randToBest1`).

# Example

```jldoctest
julia> f(x) = sum(x.^2)
f (generic function with 1 method)

julia> optimize(f, [-1 -1 -1; 1 1 1.0], DE())
+=========== RESULT ==========+
  iteration: 1000
    minimum: 0
  minimizer: [0.0, 0.0, 0.0]
    f calls: 30000
 total time: 0.0437 s
+============================+

julia> optimize(f, [-1 -1 -1; 1 1 1.0], DE(N=50, F=1.5, CR=0.8))
+=========== RESULT ==========+
  iteration: 600
    minimum: 8.68798e-25
  minimizer: [3.2777877981303293e-13, 3.7650459509488005e-13, -7.871487597385812e-13]
    f calls: 30000
 total time: 0.0319 s
+============================+
```

"""
function DE(;
    N::Int = 0,
    F = 1.0,
    CR = 0.9,
    CR_min = CR,
    CR_max = CR,
    F_min = F,
    F_max = F,
    strategy::Symbol = :rand1,
    information = Information(),
    options = Options(),
)

    parameters =
        DE(N, promote(F, CR, CR_min, CR_max, F_min, F_max)..., strategy)

    Algorithm(
        parameters,
        information = information,
        options = options,
    )

end

function update_state!(
        status,
        parameters::DE,
        problem::AbstractProblem,
        information::Information,
        options::Options,
        args...;
        kargs...
)
    population = status.population

    F = parameters.F
    CR = parameters.CR

    # stepsize
    if parameters.F_min < parameters.F_max
        F = parameters.F_min + (F_max - parameters.F_min) * rand()
    end

    if parameters.CR_min < parameters.CR_max
        CR =
            parameters.CR_min + (parameters.CR_max - parameters.CR_min) * rand()
    end

    new_vectors = reproduction(status, parameters, problem)

    # evaluate solutions
    new_solutions = create_solutions(new_vectors, problem,ε=options.h_tol)
    append!(status.population, new_solutions)

    environmental_selection!(status.population, parameters)
    status.best_sol = get_best(status.population)
end

function environmental_selection(population, parameters::DE)
    @assert length(population) == 2*parameters.N

    new_solutions = population[parameters.N+1:end]
    population = population[1:parameters.N]

    survivals = Int[]

    # select survival
    for (i, h) in enumerate(new_solutions)
        if is_better(h, population[i])
            push!(survivals, parameters.N + i)
        else
            push!(survivals, i)
        end
    end
    return survivals
end

function environmental_selection!(population, parameters::DE)
    mask = environmental_selection(population, parameters)
    ignored = ones(Bool, length(population))
    ignored[mask] .= false
    deleteat!(population, ignored)

    return
end

function initialize!(
        status,
        parameters::DE,
        problem::AbstractProblem,
        information::Information,
        options::Options,
        args...;
        kargs...
)
    D = size(problem.bounds, 2)

    if parameters.N <= 5
        parameters.N = 10 * D
    end

    if parameters.CR < 0 || parameters.CR > 1
        parameters.CR = 0.5
        options.debug &&
            @warn("CR should be from interval [0,1]; set to default value 0.5")
    end

    if options.f_calls_limit == 0
        options.f_calls_limit = 10000D
        options.debug &&
            @warn( "f_calls_limit increased to $(options.f_calls_limit)")
    end

    if options.iterations == 0
        options.iterations = div(options.f_calls_limit, parameters.N) + 1
    end

    return gen_initial_state(problem,parameters,information,options,status)
end

function final_stage!(
        status,
        parameters::DE,
        problem::AbstractProblem,
        information::Information,
        options::Options,
        args...;
        kargs...
    )
    status.final_time = time()
end


text:
## Standard Listeners
## ================

printfl = Listener(:print, x->print(x...),true,false)
printlnfl = Listener(:print, x->println(x...),true,false)


text: # YAde Yadi Yada

max_stars_repo_name: OculusMode/Julification

max_issues_repo_name: OculusMode/Julification

max_forks_repo_name: OculusMode/Julification


text:
# Introduced in 0.6
# NOTE: Deprecated Imputor docstrings use julia-repl rather than jldoctest since depwarn
# output isn't consistent across installs.
"""
    Fill(; value=mean)

Fills in the missing data with a specific value.
The current implementation is univariate, so each variable in a table or matrix will
be handled independently.

!!! Use Impute.Replace for constants or Impute.Substitue for median/mode substitution.

# Keyword Arguments
* `value::Any`: A scalar or a function that returns a scalar if
  passed the data with missing data removed (e.g, `mean`)

# Example
```julia-repl
julia> using Impute: Fill, impute

julia> M = [1.0 2.0 missing missing 5.0; 1.1 2.2 3.3 missing 5.5]
2×5 Matrix{Union{Missing, Float64}}:
 1.0  2.0   missing  missing  5.0
 1.1  2.2  3.3       missing  5.5

julia> impute(M, Fill(); dims=:rows)
2×5 Matrix{Union{Missing, Float64}}:
 1.0  2.0  2.66667  2.66667  5.0
 1.1  2.2  3.3      3.025    5.5
```
"""
struct Fill{T} <: Imputor
    value::T

    function Fill(value::T) where T
        Base.depwarn(
            "Impute.Fill is deprecated in favour of Impute.Replace for constants and " *
            "Impute.Substitute for calculating summary statistics over non-missing data.",
            :Fill
        )
        return new{T}(value)
    end
end

# TODO: Switch to using Base.@kwdef on 1.1
Fill(; value=mean) = Fill(value)

function _impute!(data::AbstractArray{Union{T, Missing}}, imp::Fill) where T
    fill_val = if isa(imp.value, Function)
        available = skipmissing(data)

        if isempty(available)
            @debug "Cannot apply fill function $(imp.value) as all values are missing"
            return data
        else
            imp.value(available)
        end
    else
        imp.value
    end

    for i in eachindex(data)
        if ismissing(data[i])
            data[i] = fill_val
        end
    end

    return data
end

# A couple utility methods to avoid messing up var and obs dimensions
# NOTE: We aren't deprecating these as they were always internal function that weren't
# intended for public use.
obsdim(dims::Int) = dims
vardim(dims::Int) = dims == 1 ? 2 : 1  # dims is obsdims, so we want the other one.

function obswise(data::AbstractMatrix; dims=1)
    return (selectdim(data, obsdim(dims), i) for i in axes(data, obsdim(dims)))
end

function varwise(data::AbstractMatrix; dims=2)
    return (selectdim(data, vardim(dims), i) for i in axes(data, vardim(dims)))
end

function filterobs(f::Function, data::AbstractMatrix; dims=1)
    mask = [f(x) for x in obswise(data; dims=dims)]
    return dims == 1 ? data[mask, :] : data[:, mask]
end

function filtervars(f::Function, data::AbstractMatrix; dims=2)
    mask = [f(x) for x in varwise(data; dims=dims)]
    return dims == 1 ? data[:, mask] : data[mask, :]
end

"""
    DropObs()

Removes missing observations from the `AbstractArray` or `Tables.table`
provided.

!!! Use `Impute.filter` instead

# Example
```julia-repl
julia> using Impute: DropObs, impute

julia> M = [1.0 1.1; 2.0 2.2; missing 3.3; missing missing; 5.0 5.5]
5×2 Matrix{Union{Missing, Float64}}:
 1.0       1.1
 2.0       2.2
  missing  3.3
  missing   missing
 5.0       5.5

julia> impute(M, DropObs())
3×2 Matrix{Union{Missing, Float64}}:
 1.0  1.1
 2.0  2.2
 5.0  5.5
```
"""
struct DropObs <: Imputor
    function DropObs()
        Base.depwarn(
            "Impute.DropObs is deprecated in favour of the more general Impute.Filter.",
            :DropObs
        )
        return new()
    end
end

# Special case impute! for vectors because we know filter! will work
impute!(data::Vector, imp::DropObs) = Base.filter!(!ismissing, data)

function impute!(data::Vector{<:NamedTuple}, imp::DropObs)
    return Base.filter!(r -> all(!ismissing, propertyvalues(r)), data)
end

impute(data::AbstractVector, imp::DropObs) = Base.filter(!ismissing, data)

function impute(data::Vector{<:NamedTuple}, imp::DropObs)
    return Base.filter(r -> all(!ismissing, propertyvalues(r)), data)
end

function impute(data::AbstractMatrix{Union{T, Missing}}, imp::DropObs; dims=1) where T
    return filterobs(obs -> all(!ismissing, obs), data; dims=dims)
end

function impute(table, imp::DropObs)
    istable(table) || throw(MethodError(impute, (table, imp)))
    rows = Tables.rows(table)

    # Unfortunately, we'll need to construct a new table
    # since Tables.rows is just an iterator
    filtered = Iterators.filter(rows) do r
        all(!ismissing, propertyvalues(r))
    end

    table = materializer(table)(filtered)
    return table
end

"""
    DropVars()

Finds variables with too many missing values in a `AbstractMatrix` or
`Tables.table` and removes them from the input data.

!!! Use `Impute.filter` instead

# Examples
```julia-repl
julia> using Impute: DropVars, impute

julia> M = [1.0 2.0 missing missing 5.0; 1.1 2.2 3.3 missing 5.5]
2×5 Matrix{Union{Missing, Float64}}:
 1.0  2.0   missing  missing  5.0
 1.1  2.2  3.3       missing  5.5

julia> impute(M, DropVars())
2×3 Matrix{Union{Missing, Float64}}:
 1.0  2.0  5.0
 1.1  2.2  5.5
```
"""
struct DropVars <: Imputor
    function DropVars()
        Base.depwarn(
            "Impute.DropVars is deprecated in favour of the more general Impute.Filter.",
            :DropVars
        )
        return new()
    end
end

function impute!(data::Vector{<:NamedTuple}, imp::DropVars)
    return materializer(data)(impute(Tables.columns(data), imp))
end

function impute(data::AbstractMatrix{Union{T, Missing}}, imp::DropVars; dims=1) where T
    return filtervars(data; dims=dims) do vars
        all(!ismissing, vars)
    end
end

function impute(table, imp::DropVars)
    istable(table) || throw(MethodError(impute, (table, imp)))
    cols = Tables.columns(table)

    cnames = Iterators.filter(propertynames(cols)) do cname
        all(!ismissing, getproperty(cols, cname))
    end

    selected = TableOperations.select(table, cnames...)
    table = materializer(table)(selected)
    return table
end

function impute!(data::AbstractMatrix{Union{T, Missing}}, imp::Union{DropObs, DropVars}) where T
    data = impute(data, imp)
    return data
end

impute!(data, imp::Union{DropObs, DropVars}) = impute(data, imp)

@deprecate impute(data, C::Chain) C(data) false
@deprecate impute!(data, C::Chain) C(data) false


text:
inpt = """
jio a, +18
inc a
tpl a
inc a
tpl a
tpl a
tpl a
inc a
tpl a
inc a
tpl a
inc a
inc a
tpl a
tpl a
tpl a
inc a
jmp +22
tpl a
inc a
tpl a
inc a
inc a
tpl a
inc a
tpl a
inc a
inc a
tpl a
tpl a
inc a
inc a
tpl a
inc a
inc a
tpl a
inc a
inc a
tpl a
jio a, +8
inc b
jie a, +4
tpl a
inc a
jmp +2
hlf a
jmp -7 """

a = 1
b = 0
instructs = split(inpt, '\n')
i = 1
while true
	items = split(instructs[i])
	if items[1] == "inc"
		eval(parse(items[2]*" += 1"))
	elseif items[1] == "tpl"
		eval(parse(items[2]*" *= 3"))
	elseif items[1] == "hlf"
		eval(parse(items[2]*" = div("*items[2]*", 2)"))
	elseif items[1] == "jmp"
		j = parse(Int, items[2])
		if i+j > length(instructs)
			break
		end
		i += j
		continue
	elseif items[1] == "jie"
		j = parse(Int, items[3])
		if eval(parse(" iseven("*items[2][1:end-1]*")"))
			if i+j > length(instructs)
				break
			end
			i += j
			continue
		end
	elseif items[1] == "jio"
		j = parse(Int, items[3])
		if eval(parse(items[2][1:end-1]*" == 1"))
			if i+j > length(instructs)
				break
			end
			i += j
			continue
		end
	end
	i += 1
end
println("a = $a\tb = $b")

max_stars_repo_name: jstrube/adventofcode

max_issues_repo_name: jstrube/adventofcode

max_forks_repo_name: jstrube/adventofcode


text:
#############################################################################
# types.jl
# Defines SubmodFunc and Sets.
# SubmodFunc subtypes of Convex.AbstractExpr, and is subtyped by atoms which
# are combinatorial functions.
# Each type which subtypes SubmodFunc must have:
#
## head::Symbol                  -- a symbol such as :vecnorm, :+ etc
## children::(AbstractExpr,)     -- the expressions on which the current expression
##                               -- is operated
## id_hash::UInt64               -- identifier hash, can be a hash of children
##                                  or a unique identifier of the object
## size::(Int, Int)              -- size of the resulting expression.
## setvariables::SetVariable      -- the set variables of this combinatorial function.
#
#
# In addition, each atom must implement the following functions:
## sign: Returns the sign of the result of the expression;
## monotonicity: The monotonicity of the arguments with respect to the function,
##      i.e if the argument is nondecreasing, will the function be nonincreasing
##      or nondecreasing? eg. negate(x) will have Nonincreasing monotonicity;
## exprmodularity: Returns the modularity of the result of the expression;
## evaluate: Evaluates the value of the expression, assuming the problem has been
##           solved.
#
# Sets is subtyped by combinatorial sets.
#############################################################################

export Values, ValOrNothing
export SubmodFunc, CombiSet, AllCombiSet, ContiSet, SCOPEModel

# Type of values
const Values = Union{Number, AbstractArray}
const ValOrNothing = Union{Values, Nothing}

### Combinatorial functions
abstract type SubmodFunc <: AbstractExpr end

### Combinatorial sets
abstract type CombiSet <: AbstractExpr end
### All combinatorial Sets
const AllCombiSet = Union{AbstractArray, CombiSet}

### Sets on continuous variables
abstract type ContiSet <: AbstractExpr end

### Problem models
abstract type SCOPEModel end

# only works for expressions with one variable
function evaluate(f::AbstractExpr, w::AbstractArray)
  var = get_v(f)
  if isa(var[1], Variable)
    var[1].value = w
  else
    if w != []
      if isa(w, AbstractArray{Int}) == false
        error("The elements assigned must be integers.")
      elseif maximum(w) > maximum(var[1].baseset) || minimum(w) < minimum(var[1].baseset)
        error("The elements assigned exceeds the cardinality of the base set of the variable.")
      else
        var[1].elements = w
      end
    else
      var[1].elements = []
    end
  end
  evaluate(f)
end

function evaluate(f::AbstractExpr, w::Values...)
  var = get_v(f)
  for i = 1:length(w)
    if isa(var[1], Variable)
      var[1].value = w
    else
      if isa(w, AbstractArray{Int}) == false
        error("The elements assigned must be integers.")
      elseif maximum(w) > maximum(var[1].baseset) || minimum(w) < minimum(var[1].baseset)
        error("The elements assigned exceeds the cardinality of the base set of the variable.")
      else
        var[1].elements = w
      end
    end
  end
  evaluate(f)
end


text:
using Pkg

Pkg.update()
Pkg.add("DataFrames")
Pkg.add("StatsBase")
Pkg.add("Combinatorics")
Pkg.add("Distances")
Pkg.add("Clustering")
Pkg.add("FFTW")
Pkg.add("Query")
Pkg.add("HTTP")
Pkg.add("JSON")
Pkg.add("Glob")

Pkg.add("Conda")
using Conda
Conda.add("scikit-learn")
Conda.add("pyclustering", channel="conda-forge")
Conda.update()

Pkg.add("PyCall")

pkg"activate . "

cd(@__DIR__)
ROOT_PATH = pwd()
push!(LOAD_PATH, ROOT_PATH, "src")

pkg"instantiate"
pkg"precompile"


text:
StructTypes.StructType(::Type{Coord}) = StructTypes.ArrayType()

StructTypes.names(::Type{Facet}) = (
    (:normal, :normal),
    (:loop, :loop)
)
StructTypes.StructType(::Type{Facet}) = StructTypes.Struct()

StructTypes.names(::Type{Design}) = (
    (:facets, :facets),
)
StructTypes.StructType(::Type{Design}) = StructTypes.Struct()


text:
"""
    function decodePoint(point::Int32, scale::Float32, boundingBoxMin::Float32) -> Float32

decodes a point of a bin file (https://github.com/PropellerAero/potree-propeller-private/blob/master/docs/file_format.md )
"""
function decodePoint(point::Int32, scale::Float64, boundingBoxMin::Float64)::Float64
	return (point * scale) + boundingBoxMin
end

"""
binToPointsUgly(fname::String) -> nothing

First version of decoder doesn't return nothing.
It just shows the decoded points and Binary code for rgba
"""
function binToPointsUgly(fname::String)
	read(fname, data)

	#divido in pacchetti di 4 * 4 byte
	rawdata = reshape(data, (16, div(length(data), 16)))

	for i in 1:size(rawdata,2)
		x = parse(Int32, bitstring(UInt8(rawdata[4,i]))*bitstring(UInt8(rawdata[3,i]))*bitstring(UInt8(rawdata[2,i]))*bitstring(UInt8(rawdata[1,i])); base=2)
		y = parse(Int32, bitstring(UInt8(rawdata[8,i]))*bitstring(UInt8(rawdata[7,i]))*bitstring(UInt8(rawdata[6,i]))*bitstring(UInt8(rawdata[5,i])); base=2)
		z = parse(Int32, bitstring(UInt8(rawdata[12,i]))*bitstring(UInt8(rawdata[11,i]))*bitstring(UInt8(rawdata[10,i]))*bitstring(UInt8(rawdata[9,i])); base=2)
		r = parse(Int32, bitstring(UInt8(rawdata[13,i])))
		g = parse(Int32, bitstring(UInt8(rawdata[14,i])))
		b = parse(Int32, bitstring(UInt8(rawdata[15,i])))
		a = parse(Int32, bitstring(UInt8(rawdata[16,i])))
		@show decodeX(x), decodeY(y), decodeZ(z)
		@show r,g,b,a
	end
end

"""
bin2points(fname::String) -> Array{Float64,2}(undef, 3, 0)

Returns every point of the passed bin file.
"""
function bin2points(fname::String, cloud_metadata::CloudMetadata)::Array{Float64,2}
	data = read(fname)
	scale = cloud_metadata.scale
	bb = cloud_metadata.boundingBox

	#divido in pacchetti di 4 * 4 byte
	rawdata = reshape(data, (16, div(length(data), 16)))

	allPoints = Array{Float64,2}(undef, 3, 0)

	for i in 1:size(rawdata,2)
		x = parse(Int32, bitstring(UInt8(rawdata[4,i]))*bitstring(UInt8(rawdata[3,i]))*bitstring(UInt8(rawdata[2,i]))*bitstring(UInt8(rawdata[1,i])); base=2)
		y = parse(Int32, bitstring(UInt8(rawdata[8,i]))*bitstring(UInt8(rawdata[7,i]))*bitstring(UInt8(rawdata[6,i]))*bitstring(UInt8(rawdata[5,i])); base=2)
		z = parse(Int32, bitstring(UInt8(rawdata[12,i]))*bitstring(UInt8(rawdata[11,i]))*bitstring(UInt8(rawdata[10,i]))*bitstring(UInt8(rawdata[9,i])); base=2)

		allPoints = hcat(allPoints, vcat(decodePoint(x, scale, bb.x_min),decodePoint(y, scale, bb.y_min),decodePoint(z, scale, bb.z_min)))

	end
	return allPoints
end

"""
bin2rgb(fname::String) -> {LasIO.N0f16,2}(undef, 3, 0)

(BUGGED RETURN) TODO: parse rgba to LasIO.N0f16
"Should" return rgb(a) values of every point of the bin file.
"""
function bin2rgb(fname::String)
	data = read(fname)

	allRGB = Array{LasIO.N0f16,2}(undef, 3, 0)

	#divido in pacchetti di 4 * 4 byte
	rawdata = reshape(data, (16, div(length(data), 16)))

	for i in 1:size(rawdata,2)
		r = parse(Int32, bitstring(UInt8(rawdata[13,i])))
		g = parse(Int32, bitstring(UInt8(rawdata[14,i])))
		b = parse(Int32, bitstring(UInt8(rawdata[15,i])))
		a = parse(Int32, bitstring(UInt8(rawdata[16,i])))
		allRGB = hcat(allRGB, vcat(r,g,b))
	end
	return allRGB
end

"""
bin2pointcloud(source::String) ->  Array{Float64,2}(undef, 3, 0)

TODO: check bin2rgb + create the PC
Will return a PointCloud
"""
function bin2pointcloud(source::String)::Array{Float64,2}

	cloud_metadata = Potree.CloudMetadata(source)
	trie = Potree.potree2trie(source)
	all_files = Potree.get_all_values(trie)

	Vtot = Array{Float64,2}(undef, 3, 0)
	rgbtot = Array{LasIO.N0f16,2}(undef, 3, 0)

	l = Threads.ReentrantLock()

	Threads.@threads for fname in all_files
		partialV = bin2points(fname, cloud_metadata)
		partialRGB = bin2rgb(fname)
		Threads.lock(l)
		Vtot = hcat(Vtot, partialV)
		rgbtot = hcat(rgbtot, partialRGB)
		Threads.unlock(l)
	end
	return Vtot
end

"""
bin2pointcloudNoMultithreading(source::String) ->  Array{Float64,2}(undef, 3, 0)

TODO: check bin2rgb + create the PC
Will return a PointCloud
"""
function bin2pointcloudNoMultithreading(source::String)::Array{Float64,2}

	cloud_metadata = Potree.CloudMetadata(source)
	trie = Potree.potree2trie(source)
	all_files = Potree.get_all_values(trie)

	Vtot = Array{Float64,2}(undef, 3, 0)
	rgbtot = Array{LasIO.N0f16,2}(undef, 3, 0)

	for fname in all_files
		partialV = bin2points(fname, cloud_metadata)
		partialRGB = bin2rgb(fname)
		Vtot = hcat(Vtot, partialV)
		rgbtot = hcat(rgbtot, partialRGB)
	end
	return Vtot
end

"""
bin2parsing(fname::String) -> Array{Float64,2}(undef, 3, 0)

Returns every point of the passed bin file.
"""
function binParsing(fname::String, cloud_metadata::CloudMetadata)
	data = read(fname)
	scale = cloud_metadata.scale
	bb = cloud_metadata.boundingBox

	allPoints = Array{Float64,2}(undef, 3, 0)
	allRGB = Array{LasIO.N0f16,2}(undef, 3, 0)

	io = open(fname, "r");

	while !eof(io)
		x = Int32(read(io, UInt32))
		y = Int32(read(io, UInt32))
		z = Int32(read(io, UInt32))
		r = read(io, UInt8)
		g = read(io, UInt8)
		b = read(io, UInt8)
		a = read(io, UInt8)
		allPoints = hcat(allPoints, vcat(decodePoint(x, scale, bb.x_min),decodePoint(y, scale, bb.y_min),decodePoint(z, scale, bb.z_min)))
		# allRGB = hcat(allRGB, vcat(r,g,b))
	end

	# return allPoints, allRGB
	return allPoints
end

function bin2pointcloudNEW(source::String)::PointCloud

	cloud_metadata = Potree.CloudMetadata(source)
	trie = Potree.potree2trie(source)
	all_files = Potree.get_all_values(trie)

	Vtot = Array{Float64,2}(undef, 3, 0)
	rgbtot = Array{LasIO.N0f16,2}(undef, 3, 0)

	for fname in all_files
		# partialV, partialRGB = binParsing(fname, cloud_metadata)
		partialV = binParsing(fname, cloud_metadata)
		Vtot = hcat(Vtot, partialV)
		# rgbtot = hcat(rgbtot, partialRGB)
	end
	return PointCloud(Vtot)
end


text:
function create_tree2(
#modified the default create_tree function to include the number of vertical slices as an argument
#which was by default set to 5
            z_root::FT,
            z_trunk::FT,
            z_canopy::FT,
            soil_bounds::Array{FT,1},
            air_bounds::Array{FT,1},
			N_subunit::Int,
            root_dist_par::FT
) where {FT<:AbstractFloat}
    # determine how many layers in roots and canopy
    #=
	_n_root  = 0;
	_r_index = Int[]
	for i in eachindex(soil_bounds)
        _z = soil_bounds[i];
        if _z > z_root
            _n_root += 1;
            push!(_r_index, i);
        else
            break
        end
    end
	=#
	#_n_root = 4;
	#_r_index = [4,5,6,7]

	#distribute roots in lower half of soil column
	#_r_index = [i for i in eachindex(soil_bounds[1:(end-1)]) if abs(soil_bounds[i]) >= abs(soil_bounds[end]/2)];
	_r_index = eachindex(soil_bounds[1:(end-1)]);

	_n_root = length(_r_index);

    _n_canopy = 0;
    _c_index  = Int[];
    for i in eachindex(air_bounds)
        _z0 = air_bounds[i];
        _z1 = air_bounds[i+1];
        if _z0 <= z_trunk < z_canopy <= _z1
            _n_canopy += 1;
            push!(_c_index, i);
            break
        elseif _z0 < z_trunk < _z1 < z_canopy
            _n_canopy += 1;
            push!(_c_index, i);
        elseif z_trunk <= _z0 < _z1 < z_canopy
            _n_canopy += 1;
            push!(_c_index, i);
        elseif z_trunk <= _z0 <= z_canopy < _z1
            _n_canopy += 1;
            push!(_c_index, i-1);
            break
        elseif z_canopy <= _z0 < _z1
            break
        end
    end

    # create evenly distributed root system for now
    _roots = RootHydraulics{FT}[];

	soil_thick = abs.(soil_bounds[2:end] - soil_bounds[1:(end-1)]);
	soil_mid = abs.(soil_bounds[2:end] + soil_bounds[1:(end-1)])/2;
	#root_dist_len = exp.(-1*abs(2/z_root)*soil_mid);
	root_dist_len = exp.(-root_dist_par*soil_mid);
	root_dist_layer = root_dist_len .* soil_thick / sum(root_dist_len .* soil_thick);

	for i in _r_index
        _Δh = abs(soil_bounds[i+1] + soil_bounds[i]) / 2;
		dZ = abs(soil_bounds[i+1]) - abs(soil_bounds[i]);
		zfrac = dZ / abs(z_root);
		#=
		if i <= 4
			zfrac = 0;
		else
			zfrac = 1/4;
		end
		=#
        _rt = RootHydraulics{FT}(N=N_subunit,
                    #area=1*zfrac,
					area=FT(root_dist_layer[i]),
                    k_max=25,#*zfrac,
                    k_rhiz=5e14,#*zfrac,
                    Δh=_Δh);
        push!(_roots, _rt);
    end

	#=
	for i in _r_index
        _Δh = abs(soil_bounds[i+1] + soil_bounds[i]) / 2;
        _rt = RootHydraulics{FT}(N=N_subunit,
                    area=1/_n_root,
                    k_max=25/_n_root,
                    k_rhiz=5e14/_n_root,
                    Δh=_Δh);
        push!(_roots, _rt);
    end
	=#

    # create Trunk
    _k_trunk = z_canopy / z_trunk * 50;
    _trunk   = StemHydraulics{FT}(N=N_subunit,
                    k_max=_k_trunk,
                    Δh=z_trunk);

    # create Branches
    _k_branch = z_canopy / (z_canopy - z_trunk) * 50;
    _branch   = StemHydraulics{FT}[];
    for i in _c_index
        _Δh = (air_bounds[i] + max(air_bounds[i+1], z_trunk)) / 2 - z_trunk;
        _st = StemHydraulics{FT}(N=N_subunit,
                    area=1/_n_canopy,
                    k_max=_k_branch/_n_canopy,
                    Δh=_Δh);
        push!(_branch, _st);
    end

    # create leaves
    _leaves = [LeafHydraulics{FT}(N=N_subunit,area=1500/_n_canopy) for i in 1:_n_canopy];

    # return plant
    return TreeLikeOrganism{FT}(_n_root,
                                _n_canopy,
                                _roots,
                                _trunk,
                                _branch,
                                _leaves,
                                _r_index,
                                _c_index,
                                zeros(FT,_n_root),
                                zeros(FT,_n_root),
                                zeros(FT,_n_root))
end


text:
#include("./mujoco_common.jl")
#include("./mj_common.jl")

using MuJoCo

############################################################# mujoco startup

val = mj.activate("/home/klowrey/.mujoco/mjkey.txt")
println("Acivation: $val")

modelfile = joinpath(dirname(@__FILE__), "humanoid.xml")
pm = mj.loadXML(modelfile, C_NULL)

if pm == nothing
   println("Bad load")
   quit()
end

pd = mj.makeData(pm)

m = unsafe_load(pm)
d = unsafe_load(pd)

#for i=1:100
#   mj.step(pm, pd)
#end

function data(pd, field::Symbol, value)
   d = unsafe_load(pd) # memory copy, without pointers
   setfield!(d, field, value)
end

function data(pd, field::Symbol)
   d = unsafe_load(pd) # memory copy, without pointers
   return getfield(d, field)
end

# function for struct infos
#structinfo(T) = [(fieldoffset(T,i), fieldname(T,i), fieldtype(T,i)) for i = 1:nfields(T)]
structinfo(T) = Dict(fieldname(T,i)=>(fieldoffset(T,i),  fieldtype(T,i)) for i = 1:nfields(T))

dinfo = structinfo(mj.Data)
minfo = structinfo(mj.Model)
oinfo = structinfo(mj.Option)

optionoffset = minfo[:opt][1] # bit offset
viscosityoffset = oinfo[:viscosity][1] # bit offset

# sets viscosity in model.options
#unsafe_store!(convert(Ptr{mjtNum}, (pm + optionoffset + viscosityoffset)), 1234.5678)

mtypes = [fieldtype(mj.Model, i) for i=1:nfields(mj.Model)]
dtypes = [fieldtype(mj.Data, i) for i=1:nfields(mj.Data)]

# BROKEN::
function update_model(m::mjModel, field::Symbol, val::mjtNum)
   setfield!(m, field, val)
   return m
end
function update_model(m::mjModel, field::Symbol, val::Integer)
   setfield!(m, field, convert(Cint, val))
   return m
end

# convert to macro?
#function update_model(pm::Ptr{mjModel}, field::Symbol, val::Integer)
#   f_off, f_type = minfo[field] # get offset & type
#   #@assert f_type == typeof(val)
#   # pointer math in julia is bytes, not of type offsets
#   unsafe_store!(convert(Ptr{f_type}, (pm + f_off)), convert(Cint, val))
#end

# pointer math in julia is bytes, not of type size
function update_model(pm::Ptr{mjModel}, offset::Integer, val::Integer)
   unsafe_store!(convert(Ptr{Cint}, (pm + offset)), convert(Cint, val))
end
function update_model(pm::Ptr{mjtNum}, offset::Integer, val::mjtNum)
   unsafe_store!(convert(Ptr{mjtNum}, (pm + offset)), val)
end

function update_model(pm::Ptr{mjModel}, field::Symbol, val::Union{Integer, mjtNum})
   f_off, f_type = minfo[field] # get offset & type
   @assert f_type == typeof(val) || f_type == typeof(convert(Cint, val))

   update_model(pm, f_off, convert(f_type, val)) # so that we can write ints
end

function mutatemodel!(pm)
   m = unsafe_load(pm)
   update_model(m, :nq, 222)
   unsafe_store!(pm, m)
end

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

function bad_map_data(pd::Ptr{mj.Data})
   # FUCKING NOPE
   c_data = unsafe_load(pd)
   mydata = unsafe_load( Ptr{jlData}(Libc.malloc(sizeof(jlData))) ) # make our data
   mydata.d = pd

   ptr_fields = intersect( fieldnames(jlData), fieldnames(mj.Data) )

   for f in ptr_fields
      #qpos = unsafe_wrap(Array, d.qpos, m.nq)
      mydata.f = unsafe_wrap(Array, c_data.f, f_size[f])
   end
end

function map_data(pm::Ptr{mjModel}, pd::Ptr{mj.Data})
   c_data = unsafe_load(pd)
   #mydata = unsafe_load( Ptr{jlData}(Libc.malloc(sizeof(jlData))) ) # make our data
   #mydata.d = pd
   args = Array{Any}(1)
   args[1] = pd

   ptr_fields = intersect( fieldnames(jlData), fieldnames(mjData) )

   f_size = mj.getdatasize(unsafe_load(pm), unsafe_load(pd))

   for f in ptr_fields
      push!(args, unsafe_wrap(Array, getfield(c_data, f), f_size[f]) )
   end
   return jlData(args...)
end


text:
# Tests for the PArray type.
# Author: Noah C. Benson <n@nben.net>

@testset "PArray" begin
    numops = 100
    @testset "1D" begin
        a = Real[]
        p = Air.PVector{Real}()
        ops = [:push, :pushfirst, :pop, :popfirst, :set, :get]
        for ii in 1:numops
            q = rand(ops)
            if q == :push
                x = rand(Float64)
                push!(a, x)
                p = Air.push(p, x)
            elseif q == :pop && length(a) > 0
                x = pop!(a)
                @test x == p[end]
                p = Air.pop(p)
            elseif q == :pushfirst
                x = rand(Float64)
                pushfirst!(a, x)
                p = Air.pushfirst(p, x)
            elseif q == :popfirst && length(a) > 0
                x = popfirst!(a)
                @test x == p[1]
                p = Air.popfirst(p)
            elseif q == :set && length(a) > 0
                k = rand(1:length(a))
                v = rand(Float64)
                a[k] = v
                p = Air.setindex(p, v, k)
            elseif q == :get && length(a) > 0
                k = rand(1:length(a))
                @test a[k] == p[k]
            end
            @test a == p
            @test size(a) == size(p)
        end
    end
    @testset "2D" begin
        # These should be more fleshed out, but for now, we can do just a
        # few simple tests.
        a = Array(reshape(1:200, (20,10)))
        p = PArray(a)
        @test p == a
        @test p[5,8] == a[5,8]
        @test p[:,4] == a[:,4]
        @test p[9,:] == a[9,:]
        p1 = Air.setindex(p, -5, 15, 8)
        @test p1 != a
        @test p1[:,4] == a[:,4]
        @test p1[9,:] == a[9,:]
        @test p1[:,8] != a[:,8]
        @test p1[15,:] != a[15,:]
    end
    @testset "3D" begin
        a = Array(reshape(1:1000, (20,10,5)))
        p = PArray(a)
        @test p == a
        @test p[5,8,2] == a[5,8,2]
        @test p[:,4,1] == a[:,4,1]
        @test p[9,:,3] == a[9,:,3]
        @test p[:,2,4] == a[:,2,4]
        p1 = Air.setindex(p, -5, 15, 8, 5)
        @test p1 != a
        @test p1[9,:,:] == a[9,:,:]
        @test p1[:,3,:] == a[:,3,:]
        @test p1[:,:,4] == a[:,:,4]
        @test p1[15,:,:] != a[15,:,:]
        @test p1[:,8,:] != a[:,8,:]
        @test p1[:,:,5] != a[:,:,5]
    end
end


text:
module TimeSteppers

export
    AdamsBashforthTimeStepper,
    time_step!,
    compute_w_from_continuity!

using Oceananigans: device

using GPUifyLoops: @launch, @loop, @unroll

import Oceananigans: TimeStepper

using Oceananigans: AbstractGrid, Model, Tendencies, tracernames,
                    @hascuda, CPU, GPU, launch_config, datatuples, datatuple,
                    @loop_xyz,

                    buoyancy_perturbation,
                    x_f_cross_U, y_f_cross_U, z_f_cross_U,

                    fill_halo_regions!, apply_z_bcs!, apply_y_bcs!,

                    time_to_run

using Oceananigans.SurfaceWaves: x_curl_Uˢ_cross_U, y_curl_Uˢ_cross_U, z_curl_Uˢ_cross_U,
                                 ∂t_uˢ, ∂t_vˢ, ∂t_wˢ

@hascuda using CUDAnative, CUDAdrv, CuArrays

using Oceananigans.Operators
using Oceananigans.Solvers
using Oceananigans.Diagnostics
using Oceananigans.OutputWriters

using Oceananigans.Solvers: solve_poisson_3d!, PoissonBCs, PPN, PNN
using Oceananigans.Diagnostics: run_diagnostic
using Oceananigans.OutputWriters: write_output

using ..TurbulenceClosures: calculate_diffusivities!, ∂ⱼ_2ν_Σ₁ⱼ, ∂ⱼ_2ν_Σ₂ⱼ, ∂ⱼ_2ν_Σ₃ⱼ, ∇_κ_∇c

"""
    TimeStepper(name, args...)

Returns a timestepper with name `name`, instantiated with `args...`.

Example
=======

julia> stepper = TimeStepper(:AdamsBashforth, Float64, CPU(), grid, tracernames)
"""
function TimeStepper(name::Symbol, args...)
    fullname = Symbol(name, :TimeStepper)
    return eval(Expr(:call, fullname, args...))
end

# Fallback
TimeStepper(stepper, args...) = stepper

"""Returns the arguments passed to boundary conditions functions."""
boundary_condition_function_arguments(model) =
    (model.clock.time, model.clock.iteration, datatuple(model.velocities),
     datatuple(model.tracers), model.parameters)

#####
##### Time-stepping kernels/function that are independent of the TimeStepper
#####

"""
    time_step!(model; Nt, Δt, kwargs...)

Step forward `model` `Nt` time steps with step size `Δt`.

The kwargs are passed to the `time_step!` function specific to `model.timestepper`.
"""
time_step!(model; Nt, Δt, kwargs...) = time_step!(model, Nt, Δt; kwargs...)

function time_step!(model, Nt, Δt; kwargs...)

    if model.clock.iteration == 0
        [ run_diagnostic(model, diag) for diag in values(model.diagnostics) ]
        [ write_output(model, out)    for out  in values(model.output_writers) ]
    end

    for n in 1:Nt
        time_step!(model, Δt; kwargs...)

        [ time_to_run(model.clock, diag) && run_diagnostic(model, diag) for diag in values(model.diagnostics) ]
        [ time_to_run(model.clock, out) && write_output(model, out) for out in values(model.output_writers) ]
    end

    return nothing
end

"""
    time_step_precomputations!(diffusivities, pressures, velocities, tracers, model)

Perform precomputations necessary for an explicit timestep or substep.
"""
function time_step_precomputations!(diffusivities, pressures, velocities, tracers, model)

    fill_halo_regions!(merge(velocities, tracers), model.boundary_conditions.solution, model.architecture,
                       model.grid, boundary_condition_function_arguments(model)...)

    calculate_diffusivities!(diffusivities, model.architecture, model.grid, model.closure, model.buoyancy,
                             velocities, tracers)

    # Diffusivities share bcs with pressure:
    fill_halo_regions!(diffusivities, model.boundary_conditions.pressure, model.architecture, model.grid)

    @launch(device(model.architecture), config=launch_config(model.grid, :xy),
            update_hydrostatic_pressure!(pressures.pHY′, model.grid, model.buoyancy, tracers))

    fill_halo_regions!(pressures.pHY′, model.boundary_conditions.pressure, model.architecture, model.grid)

    return nothing
end

"""
    calculate_tendencies!(diffusivities, pressures, velocities, tracers, model)

Calculate the interior and boundary contributions to tendency terms without the
contribution from non-hydrostatic pressure.
"""
function calculate_tendencies!(tendencies, velocities, tracers, pressures, diffusivities, model)

    calculate_interior_source_terms!(tendencies, model.architecture, model.grid, model.coriolis, model.buoyancy,
                                     model.surface_waves, model.closure, velocities, tracers, pressures.pHY′,
                                     diffusivities, model.forcing, model.parameters, model.clock.time)

    calculate_boundary_source_terms!(tendencies, model.boundary_conditions.solution, model.architecture,
                                     model.grid, boundary_condition_function_arguments(model)...)

    return nothing
end

"""
    calculate_pressure_correction!(nonhydrostatic_pressure, Δt, tendencies, velocities, model)

Calculate the (nonhydrostatic) pressure correction associated `tendencies`, `velocities`, and step size `Δt`.
"""
function calculate_pressure_correction!(nonhydrostatic_pressure, Δt, tendencies, velocities, model)
    velocity_tendencies = (u=tendencies.u, v=tendencies.v, w=tendencies.w)

    velocity_tendency_boundary_conditions = (u=model.boundary_conditions.tendency.u,
                                             v=model.boundary_conditions.tendency.v,
                                             w=model.boundary_conditions.tendency.w)

    fill_halo_regions!(velocity_tendencies, velocity_tendency_boundary_conditions, model.architecture, model.grid)

    @launch(device(model.architecture), config=launch_config(model.grid, :xyz),
            calculate_poisson_right_hand_side!(model.poisson_solver.storage, model.architecture, model.grid,
                                               model.poisson_solver.bcs, velocities, tendencies, Δt))

    solve_for_pressure!(nonhydrostatic_pressure, model.architecture, model.grid, model.poisson_solver,
                        model.poisson_solver.storage)

    fill_halo_regions!(nonhydrostatic_pressure, model.boundary_conditions.pressure, model.architecture, model.grid)

    return nothing
end

calculate_pressure_correction!(::Nothing, args...) = nothing

"""
    calculate_explicit_substep!(tendencies, velocities, tracers, pressures, diffusivities, model)

Calculate the initial and explicit substep of the two-step fractional step method with pressure correction.
"""
function calculate_explicit_substep!(tendencies, velocities, tracers, pressures, diffusivities, model)
    time_step_precomputations!(diffusivities, pressures, velocities, tracers, model)
    calculate_tendencies!(tendencies, velocities, tracers, pressures, diffusivities, model)
    return nothing
end

"""
    complete_pressure_correction_step!(velocities, Δt, tracers, pressures, tendencies, model)

After calculating the pressure correction, complete the pressure correction step by updating
the velocity and tracer fields.
"""
function complete_pressure_correction_step!(velocities, Δt, tracers, pressures, tendencies, model)
    update_solution!(velocities, tracers, model.architecture, model.grid, Δt, tendencies, pressures.pNHS)

    velocity_boundary_conditions = (u=model.boundary_conditions.solution.u,
                                    v=model.boundary_conditions.solution.v,
                                    w=model.boundary_conditions.solution.w)

    # Recompute vertical velocity w from continuity equation to ensure incompressibility
    fill_halo_regions!(velocities, velocity_boundary_conditions, model.architecture, model.grid,
                       boundary_condition_function_arguments(model)...)

    compute_w_from_continuity!(model)

    return nothing
end

include("kernels.jl")
include("adams_bashforth.jl")

end # module


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

mutable struct PolicyAssignment <: SwaggerModel
    properties::Any # spec type: Union{ Nothing, PolicyAssignmentProperties } # spec name: properties
    id::Any # spec type: Union{ Nothing, String } # spec name: id
    type::Any # spec type: Union{ Nothing, String } # spec name: type
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    sku::Any # spec type: Union{ Nothing, PolicySku } # spec name: sku
    location::Any # spec type: Union{ Nothing, String } # spec name: location
    identity::Any # spec type: Union{ Nothing, Identity } # spec name: identity

    function PolicyAssignment(;properties=nothing, id=nothing, type=nothing, name=nothing, sku=nothing, location=nothing, identity=nothing)
        o = new()
        validate_property(PolicyAssignment, Symbol("properties"), properties)
        setfield!(o, Symbol("properties"), properties)
        validate_property(PolicyAssignment, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        validate_property(PolicyAssignment, Symbol("type"), type)
        setfield!(o, Symbol("type"), type)
        validate_property(PolicyAssignment, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(PolicyAssignment, Symbol("sku"), sku)
        setfield!(o, Symbol("sku"), sku)
        validate_property(PolicyAssignment, Symbol("location"), location)
        setfield!(o, Symbol("location"), location)
        validate_property(PolicyAssignment, Symbol("identity"), identity)
        setfield!(o, Symbol("identity"), identity)
        o
    end
end # type PolicyAssignment

const _property_map_PolicyAssignment = Dict{Symbol,Symbol}(Symbol("properties")=>Symbol("properties"), Symbol("id")=>Symbol("id"), Symbol("type")=>Symbol("type"), Symbol("name")=>Symbol("name"), Symbol("sku")=>Symbol("sku"), Symbol("location")=>Symbol("location"), Symbol("identity")=>Symbol("identity"))
const _property_types_PolicyAssignment = Dict{Symbol,String}(Symbol("properties")=>"PolicyAssignmentProperties", Symbol("id")=>"String", Symbol("type")=>"String", Symbol("name")=>"String", Symbol("sku")=>"PolicySku", Symbol("location")=>"String", Symbol("identity")=>"Identity")
Base.propertynames(::Type{ PolicyAssignment }) = collect(keys(_property_map_PolicyAssignment))
Swagger.property_type(::Type{ PolicyAssignment }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_PolicyAssignment[name]))}
Swagger.field_name(::Type{ PolicyAssignment }, property_name::Symbol) =  _property_map_PolicyAssignment[property_name]

function check_required(o::PolicyAssignment)
    true
end

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


text:
@testset "error handling" begin
    n = 5
    A = diagm(-1 => fill(-1.0+1.0im,n-1)) + triu(rand(n,n))
    H = GenericSchur.HessenbergArg(A,zeros(n-1))
    @test_throws ArgumentError GenericSchur.gschur!(H)

    A = triu(rand(n,n),-1)
    H = GenericSchur.HessenbergArg(A,zeros(n-1))
    Z = rand(n-1,n-1)
    @test_throws DimensionMismatch GenericSchur.gschur!(H,Z)

    @test_throws DimensionMismatch GenericSchur.gschur!(rand(5,4))
    @test_throws DimensionMismatch GenericSchur.gschur!(rand(ComplexF64,5,4))

    # WARNING: These depend on absence of surprises in stdlib. Revise as needed.
    A = rand(Int, 4,4)
    @test_throws MethodError schur!(A)
    A = rand(Int, 4,4) + im * rand(Int, 4,4)
    @test_throws MethodError schur!(A)
    A = rand(1:10, 4,4) .// rand(1:10, 4,4)
    @test_throws MethodError schur!(A)
    @test_throws MethodError hessenberg!(A)
end


text:
abstract type Assoc{I, T}; end

indices(a::Assoc) = error("`indices` not defined for $(typeof(a))")
getindex(::Assoc{I}, ::I) where {I} = error("`getindex` not defined for $(typeof(a))")
function setindex!(::Assoc{I, T}, ::T, ::I) where {I, T}
    error("`setindex` not defined for $(typeof(a))")
end

struct IndexDirect <: IndexStyle; end
struct IndexToken <: IndexStyle; end
IndexStyle(::Assoc) = IndexDirect()

# Implement this interface:
#    - `indices(::Assoc)`
#    - `getindex(::Assoc{I}, ::I) where I`
#    - `setindex!(::Assoc{I, T}, ::T, ::I) where I` (optional)
#    - `Base.IndexStyle(::Assoc)` -> `IndexDirect` or `IndexToken`
#      (if `IndexToken`, implement `gettokenindex` and `gettokenvalue`)
#
# Get for free:
#    - `indextype`, `eltype`
#    - `pairs` - iterates index=>value `Pair`s (also indexible)
#    - `getindex(::Assoc{I,T}, ::Assoc{K2, I})` returning an Assoc{K2,T}
#    - similarly for `setindex!`
#    - length, iteration, etc...
#    - similar
#
# TODO:
#    - view (same indexing behavior as `getindex`)
#    - map (values, preserves indices)
#    - filter (preserves a subset of indices)
#    - SAC.jl functions ?

indextype(::Assoc{I}) where {I} = I
indextype(::Type{<:Assoc{I}}) where {I} = I
eltype(::Assoc{<:Any, T}) where {T} = T
eltype(::Type{<:Assoc{<:Any, T}}) where {T} = T

"""
    Pairs(a::Assoc)

Returns a new `Assoc` with the same indices as `a`, where each value is now a `Pair` of
the index and value.
"""
struct Pairs{I, T, A <: Assoc{I, T}} <: Assoc{I, Pair{I, T}}
    a::A
end
Pairs(a::Assoc{I,T}) where {I, T} = Pairs{I, T, typeof(a)}(a)

indices(p::Pairs) = indices(p.a)
IndexStyle(p::Pairs) = IndexStyle(p.a)
@propagate_inbounds getindex(p::Pairs{I}, i::I) where {I} = (i => p.a[i])
@propagate_inbounds getindex(p::Pairs{<:CartesianIndex}, i::Int) = (i => p.a[i])
tokens(p::Pairs) = tokens(p.a)
gettokenindex(p::Pairs, t) = gettokenindex(p.a, t)
gettokenvalue(p::Pairs, t) = (gettokenindex(p.a, t) => gettokenvalue(p.a, t))

pairs(a::Assoc) = Pairs(a)

length(a::Assoc) = length(indices(a))

start(a::Assoc) = _start(IndexStyle(a), a)
_start(::IndexDirect, a::Assoc) = start(indices(a))
_start(::IndexToken, a::Assoc) = start(tokens(a)) #?

@propagate_inbounds next(a::Assoc, i) = _next(IndexStyle(a), a, i)
@propagate_inbounds function _next(::IndexDirect, a::Assoc, i)
    (j, i2) = next(indices(a), i)
    return (a[j], i2)
end
@propagate_inbounds function _next(::IndexToken, a::Assoc, i)
    (j, i2) = next(tokens(a), i)
    return (gettokenvalue(a, j), i2)
end

done(a::Assoc, i) = _done(IndexStyle(a), a, i)
_done(::IndexDirect, a::Assoc, i) = done(indices(a), i)
_done(::IndexToken, a::Assoc, i) = done(tokens(a), i)

@propagate_inbounds function getindex(a::Assoc, inds::Assoc)
    out = similar(a, indices(inds))
    for (i, j) in pairs(inds)
        out[i] = a[j]
    end
    return out
end

@propagate_inbounds function setindex!(a::Assoc{I,T}, v, i::I) where {I,T}
    a[i] = convert(T, v)
    return v
end

@propagate_inbounds function setindex!(a::Assoc{I,T}, v::T, inds::Associative) where {I,T}
    for i in inds
        a[i] = v
    end
    return v
end

@propagate_inbounds function setindex!(a::Assoc{I,T}, v, inds::Associative) where {I,T}
    a[inds] = convert(T, v)
end

@propagate_inbounds function setindex!(a::Assoc{I,T}, v::Associative, inds::Associative) where {I,T}
    for (i, j) in pais(inds)
        a[j] = v[i]
    end
    return v
end

# Similar works a bit like for `Array`...
similar(a::Assoc) = similar(a, eltype(a), indices(a))
similar(a::Assoc, ::Type{T}) where {T} = similar(a, T, indices(a))
similar(a::Assoc, i) = similar(a, eltype(a), i)
similar(a::Assoc, dims::Union{Integer, AbstractUnitRange}...) = similar(a, eltype(a), dims...)

# `empty` is much like `similar`, but with no defined indices (expects `push` or `setindex`
# or whatever to create indices)
empty(a::Assoc) = empty(a, indextype(a), eltype(a))
empty(a::Assoc, ::Type{T}) where {T} = empty(a, indextype(a), T)

# A couple of basic methods for now:
function Base.:(==)(a::Assoc, b::Assoc)
    if !issetequal(a, b)
        return false
    end

    for i in indices(a)
        if a[i] != b[i]
            return false
        end
    end
    return true
end

function Base.map(f, a::Assoc)
    out = similar(a, Base.promote_op(f, eltype(a)))
    for i in indices(a)
        @inbounds out[i] = f(a[i])
    end
end


text:
# Autogenerated wrapper script for Gnuastro_jll for x86_64-linux-gnu
export astarithmetic, astconvertt, astconvolve, astcosmiccal, astcrop, astfits, astmatch, astmkcatalog, astmknoise, astmkprof, astnoisechisel, astscript_sort_by_night, astsegment, aststatistics, asttable, astwarp, libgnuastro

using CFITSIO_jll
using GSL_jll
using JpegTurbo_jll
using Libtiff_jll
using LibGit2_jll
using WCS_jll
JLLWrappers.@generate_wrapper_header("Gnuastro")
JLLWrappers.@declare_executable_product(astarithmetic)
JLLWrappers.@declare_executable_product(astconvertt)
JLLWrappers.@declare_executable_product(astconvolve)
JLLWrappers.@declare_executable_product(astcosmiccal)
JLLWrappers.@declare_executable_product(astcrop)
JLLWrappers.@declare_executable_product(astfits)
JLLWrappers.@declare_executable_product(astmatch)
JLLWrappers.@declare_executable_product(astmkcatalog)
JLLWrappers.@declare_executable_product(astmknoise)
JLLWrappers.@declare_executable_product(astmkprof)
JLLWrappers.@declare_executable_product(astnoisechisel)
JLLWrappers.@declare_executable_product(astscript_sort_by_night)
JLLWrappers.@declare_executable_product(astsegment)
JLLWrappers.@declare_executable_product(aststatistics)
JLLWrappers.@declare_executable_product(asttable)
JLLWrappers.@declare_executable_product(astwarp)
JLLWrappers.@declare_library_product(libgnuastro, "libgnuastro.so.12")
function __init__()
    JLLWrappers.@generate_init_header(CFITSIO_jll, GSL_jll, JpegTurbo_jll, Libtiff_jll, LibGit2_jll, WCS_jll)
    JLLWrappers.@init_executable_product(
        astarithmetic,
        "bin/astarithmetic",
    )

    JLLWrappers.@init_executable_product(
        astconvertt,
        "bin/astconvertt",
    )

    JLLWrappers.@init_executable_product(
        astconvolve,
        "bin/astconvolve",
    )

    JLLWrappers.@init_executable_product(
        astcosmiccal,
        "bin/astcosmiccal",
    )

    JLLWrappers.@init_executable_product(
        astcrop,
        "bin/astcrop",
    )

    JLLWrappers.@init_executable_product(
        astfits,
        "bin/astfits",
    )

    JLLWrappers.@init_executable_product(
        astmatch,
        "bin/astmatch",
    )

    JLLWrappers.@init_executable_product(
        astmkcatalog,
        "bin/astmkcatalog",
    )

    JLLWrappers.@init_executable_product(
        astmknoise,
        "bin/astmknoise",
    )

    JLLWrappers.@init_executable_product(
        astmkprof,
        "bin/astmkprof",
    )

    JLLWrappers.@init_executable_product(
        astnoisechisel,
        "bin/astnoisechisel",
    )

    JLLWrappers.@init_executable_product(
        astscript_sort_by_night,
        "bin/astscript-sort-by-night",
    )

    JLLWrappers.@init_executable_product(
        astsegment,
        "bin/astsegment",
    )

    JLLWrappers.@init_executable_product(
        aststatistics,
        "bin/aststatistics",
    )

    JLLWrappers.@init_executable_product(
        asttable,
        "bin/asttable",
    )

    JLLWrappers.@init_executable_product(
        astwarp,
        "bin/astwarp",
    )

    JLLWrappers.@init_library_product(
        libgnuastro,
        "lib/libgnuastro.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# glcPointsLocal(N::Int) = [0.5(1-cos((i-1)pi/(N-1))) for i in 1:N]
# glcPointsGlobal(left::Number,right::Number,localPnts::Vector) = (right-left).*localPnts.-left

function glcPointsGlobal(left::Number,right::Number,N::Int)
    glcPointsLocal= [0.5(1-cos((i-1)pi/(N-1))) for i in 1:N]
    (right-left).*glcPointsLocal.-left
end

function numstring(x::Vector)::Float64
    pop!(cumprod(x))
end

function CofDQM(x::Vector,order::Int)
    #x: vector of sampling points;
    #order: largest derivative order;
    #return: matrix with dimension of (x,x,m)

    n = length(x)
    C = fill(0.0,(n,n,order))
    (n>order) || error("need more sampling points")

    # 1st order

    for i in 1:n, j in 1:n
        if i != j
            C[i,j,1] = numstring([((m!=i && m!=j) ? x[i]-x[m] : 1) for  m in 1:n])/numstring([((m!=j) ? x[j]-x[m] : 1) for  m in 1:n])

            # m_range_denominater = filter!(xx->xx!=i,[m for m in 1:n])
            # m_range_numerator = filter!(xx->xx!=j,m_range_denominater)
            #C[i,j,1] = numstring(filter(xx->xx!=0,[x[i]-x[m] for m in m_range_numerator]))/numstring(filter(xx->xx!=0,[x[j]-x[m] for m in m_range_denominater]))
            # C[i,j,1] = pop!(cumprod(filter(xx->xx!=0,[x[i]-x[m] for m in m_range_numerator])))/pop!(cumprod(filter(xx->xx!=0,[x[j]-x[m] for m in 1:n])))
        else

            # C[i,i,1] = sum(filter(xx->xx!=Inf,[1/(x[i]-x[m]) for m in 1:n]))
            C[i,i,1] = sum([((m!=i) ? 1/(x[i]-x[m]) : 0) for m in 1:n]) #这个慢了4个微秒

        end
    end

    # higher order
    for i in 1:n, j in 1:n, k in 2:order
        if i != j
            C[i,j,k] = k*(C[i,i,k-1]*C[i,j,1] - C[i,j,k-1]/(x[i]-x[j]))
        end
    end
    for i in 1:n, j in 1:n, k in 2:order
        if i == j
            C[i,j,k] = -sum([((m!=i) ? C[i,m,k] : 0) for m in 1:n])
        end
    end
    C
end

`
demoooo: test 1st-order approximation accuracy
`
pnts = glcPointsGlobal(-10,1,100)
# pnts = collect(0:2:10)

Coef = CofDQM(pnts,2)
f4(x::Number) = 4x^3+3x-2
f4d(x) = 12x^2+3

f4d.(pnts)-Coef[:,:,1]*f4.(pnts)

max_stars_repo_name: frydaydeep/QDEM

max_issues_repo_name: frydaydeep/QDEM

max_forks_repo_name: frydaydeep/QDEM


text:
# %%

using Binance.BinanceF,Dates, DataFrames, Plots
# %%

hr24 = Binance.get24HR()
hr24ETHBTC = Binance.get24HR("ETHBTC")

market = Binance.getMarket()
market_BNBBTC = Binance.getMarket("BNBBTC")

function getBinanceKlineDataframe(symbol; startDateTime = nothing, endDateTime = nothing, interval="1m")
    klines = Binance.getKlines(symbol; startDateTime = startDateTime, endDateTime = endDateTime, interval = interval)
    result = hcat(map(z -> map(x -> typeof(x) == String ? parse(Float64, x) : x, z), klines)...)';

    if size(result,2) == 0
        return nothing
    end

    symbolColumnData = map(x -> symbol, collect(1:size(result, 1)));
    df = DataFrame([symbolColumnData, Dates.unix2datetime.(result[:,1]/1000) ,result[:,2],result[:,3],result[:,4],result[:,5],result[:,6],result[:,8],Dates.unix2datetime.(result[:,7] / 1000),result[:,9],result[:,10],result[:,11]], [:symbol,:startDate,:open,:high,:low,:close,:volume,:quoteAVolume, :endDate, :trades, :tbBaseAVolume,:tbQuoteAVolume]);
end

dfKlines = getBinanceKlineDataframe("ETHBTC");

plot(dfKlines[:close];label="ETHBTC interval = '1m'")