text:
using SciMLBase, EllipsisNotation

"""
    NODEDataloader{T,N} <: AbstractNODEDataloader{T,N}

Struct containing batched data for sequence learning of ODEs. Can be indexed and interated over. Each batch returns a tuple (t, data(t)).

# Initialized with

    NODEData(sol::SciMLBase.AbstractTimeseriesSolution, N_length::Integer; dt=nothing, valid_set=nothing)

* `sol`: DE solution
* `N_length`: length of each batch
* `dt`: time increment that `sol` is interpolated at. If `nothing` then the `sol.t` is used as the time steps of the data
* `valid_set` if valid_set ∈ [0,1] splits the data into a train and valid set with `valid_set` of the share of the data belonging to the valid_set.

    NODEDataloader(data::AbstractArray{T,N}, t::AbstractArray{T,1}, N_length::Integer)

* `data`:: Data that is already in N_dim_1 x ... x N_t format
* `t`: time axis
* `N_length`: length of each batch

"""
struct NODEDataloader{T,U,N} <: AbstractNODEDataloader{T,U,N}
    data::AbstractArray{T,N}
    t::AbstractArray{U,1}
    N::Integer
    N_length::Integer
end

function NODEDataloader(sol::Union{SciMLBase.AbstractTimeseriesSolution, SciMLBase.AbstractDiffEqArray}, N_length::Integer; dt=nothing, valid_set=nothing)

    if isnothing(dt)
        data = DeviceArray(sol)
        t = sol.t
    else
        t = sol.t[1]:eltype(sol)(dt):sol.t[end]
        data = DeviceArray(sol(t))
    end

    if isnothing(valid_set)
        N_t = length(t)
        N = N_t - N_length

        return NODEDataloader(DeviceArray(data), t, N, N_length)
    else
        @assert 0<valid_set<1

        N_t = length(t)
        N_t_valid = Int(floor(valid_set*N_t))
        N_t_train = N_t - N_t_valid

        return NODEDataloader(DeviceArray(data[..,1:N_t_train]), t[1:N_t_train], N_t_train - N_length, N_length), NODEDataloader(DeviceArray(data[..,N_t_train+1:N_t]), t[N_t_train+1:N_t], N_t_valid - N_length, N_length)
    end
end

NODEDataloader(data::AbstractArray{T,N}, t::AbstractArray{U,1}, N_length::Integer) where {T,U,N} = NODEDataloader(DeviceArray(data), t, length(t) - N_length, N_length)
NODEDataloader(data::NODEDataloader, N_length::Integer) = NODEDataloader(data.data, data.t, N_length)

function Base.getindex(iter::NODEDataloader{T,U,N}, i::Integer) where {T,U,N}
    @assert 0 < i <= iter.N
    return (iter.t[i:i+iter.N_length-1] ,iter.data[..,i:i+iter.N_length-1])
end

function Base.iterate(iter::AbstractNODEDataloader, state=1)
    if state>iter.N
        return nothing
    else
        return (iter[state], state+1)
    end
end

Base.length(iter::AbstractNODEDataloader) = iter.N
Base.eltype(iter::AbstractNODEDataloader) = eltype(iter.data)

Base.firstindex(iter::AbstractNODEDataloader) = 1
Base.lastindex(iter::AbstractNODEDataloader) = iter.N
Base.show(io::IO,seq::NODEDataloader{T,U,N}) where {T,U,N} = print(io, "NODEData{",T,",",N,"} with ",seq.N," batches with length ",seq.N_length)

"""
    get_trajectory(data::NODEDataloader, N)

Returns a (t, x(t)) tuple of length `N`.
"""
function get_trajectory(data::NODEDataloader, N)
    ind = 1:N
    (data.t[ind], data.data[..,ind])
end

cpu(data::NODEDataloader) = NODEDataloader(Array(data.data), data.t, data.N, data.N_length)
gpu(data::NODEDataloader) = NODEDataloader(DeviceArray(data.data), data.t, data.N, data.N_length)


text:
__precompile__()

module CardinalDicts

export CardinalDict, clearindex!, keymax, isfull

include("for_indexing.jl")
include("CardinalDict.jl")
include("CardinalDict_api.jl")

end # module

ext: jl


text:
using Conda
#===============================================================================
NOTE: PyDAQmx is a conda package I built from a PyPI package,
following the directions at https://conda.io/docs/build_tutorials/pkgs.html.
===============================================================================#
Conda.add_channel("haberdashPI"); Conda.add("PyDAQmx")


text:
using HTTP
using PyCall
using NPZ
# pikcleを扱うためのPythonのモジュールの読み込み
pickle = pyimport("pickle")

url_base = "https://raw.githubusercontent.com/tomsercu/lstm/master/data/"
key_file = Dict(
    "train" => "ptb.train.txt",
    "test" => "ptb.test.txt",
    "valid" => "ptb.valid.txt"
)
save_file = Dict(
    "train" => "ptb.train.npy",
    "test" => "ptb.test.npy",
    "valid" => "ptb.valid.npy"
)
vocab_file = "ptb.vocab.pkl"

dataset_dir = dirname(@__FILE__)

function _download(file_name)
    file_path = dataset_dir * "/" * file_name
    if isfile(file_path)
        return
    end

    println("Downloading " * file_name * " ...")

    try
        HTTP.download(url_base * file_name, file_path)
    catch e
        println(e)
        return
    end

    println("Done!")
end

function load_vocab()
    vocab_path = dataset_dir * "/" * vocab_file

    # 既にツリーバンクがダウンロードされている場合
    if isfile(vocab_path)
        @pywith pybuiltin("open")(vocab_path, "rb") as f begin
            word_to_id, id_to_word = pickle.load(f)
            return word_to_id, id_to_word
        end
    end

    word_to_id = Dict()
    id_to_word = Dict()
    data_type = "train"
    file_name = key_file[data_type]
    file_path = dataset_dir * "/" * file_name

    _download(file_name)

    f = open(file_path)
    words = split(strip(replace(read(f, String), "\n" => "<eos>")))
    close(f)

    for (i, word) = enumerate(words)
        if false == get(word_to_id, word, false)
            tmp_id = length(word_to_id) + 1
            word_to_id[word] = tmp_id
            id_to_word[tmp_id] = word
        end
    end

    f = open(vocab_path, "w")
    pickle.dump((word_to_id, id_to_word), f)
    close(f)

    return word_to_id, id_to_word
end

function load_data(;data_type="train")
    #=
    パラメータ:
    data_type データの種類（"train" or "test" or "valid (val)"）

    返り値:
    なし =#
    if data_type == "val" data_type = "valid" end
    save_path = dataset_dir * "/" * save_file[data_type]

    word_to_id, id_to_word = load_vocab()

    if isfile(save_path)
        corpus = npzread(save_path)
        return corpus, word_to_id, id_to_word
    end

    file_name = key_file[data_type]
    file_path = dataset_dir * "/" * file_name
    _download(file_name)

    f = open(file_path)
    words = split(strip(replace(read(f, String), "\n" => "<eos>")))
    close(f)
    corpus = [word_to_id[w] for w = words]

    npzwrite(save_path, corpus)
    return corpus, word_to_id, id_to_word
end

for data_type in ("train", "val", "test")
    load_data(data_type=data_type)
end


text:
using StaticArrays
using UniqueVectors
using Bravais
using ExactDiag

using JLD
using Combinatorics

using Test

using SparseArrays
using LinearAlgebra
using Arpack

debug = false

@inferred ExactDiag.exp_2πiη(0//1)
@test ExactDiag.exp_2πiη(0//1) == 1
@test ExactDiag.exp_2πiη(1//2) == -1
@test ExactDiag.exp_2πiη(-1//2) == -1
@test ExactDiag.exp_2πiη(-1//1) == 1

@test ExactDiag.exp_2πiη(-1//1) ≈ exp(-2π * im)
@test ExactDiag.exp_2πiη(1//2) ≈ exp(π * im)
@test ExactDiag.exp_2πiη(1//4) ≈ exp(π/2 * im)
@test ExactDiag.exp_2πiη(1//3) ≈ exp(2π/3 * im)
@test ExactDiag.exp_2πiη(1//3) ≈ exp(2π/3 * im)

@testset "Spin half         " begin include("spin_half.jl") end
@testset "Hubbard           " begin include("hubbard.jl") end
@testset "Abelian           " begin include("abelian.jl") end
@testset "Abelian spin-flip " begin include("abelian-spinflip.jl") end
@testset "Abelian reflection" begin include("abelian-reflection.jl") end
@testset "Time evolution    " begin include("time-evolution.jl") end


text:
import Oceananigans.Grids: required_halo_size

struct IsotropicDiffusivity{TD, N, K} <: AbstractIsotropicDiffusivity{TD}
    ν :: N
    κ :: K

    function IsotropicDiffusivity{TD}(ν::N, κ::K) where {TD, N, K}
        return new{TD, N, K}(ν, κ)
    end
end

"""
    IsotropicDiffusivity([FT=Float64;] ν=0, κ=0, time_discretization = ExplicitTimeDiscretization())

Returns parameters for an isotropic diffusivity model with viscosity `ν`
and thermal diffusivities `κ` for each tracer field in `tracers`
`ν` and the fields of `κ` may be constants, arrays, fields, or functions of `(x, y, z, t)`.

`κ` may be a `NamedTuple` with fields corresponding
to each tracer, or a single number to be a applied to all tracers.
"""
function IsotropicDiffusivity(FT=Float64; ν=0, κ=0, time_discretization::TD = ExplicitTimeDiscretization()) where TD
    if ν isa Number && κ isa Number
        κ = convert_diffusivity(FT, κ)
        return IsotropicDiffusivity{TD}(FT(ν), κ)
    else
        return IsotropicDiffusivity{TD}(ν, κ)
    end
end

required_halo_size(closure::IsotropicDiffusivity) = 1

function with_tracers(tracers, closure::IsotropicDiffusivity{TD}) where TD
    κ = tracer_diffusivities(tracers, closure.κ)
    return IsotropicDiffusivity{TD}(closure.ν, κ)
end

calculate_diffusivities!(K, arch, grid, closure::IsotropicDiffusivity, args...) = nothing

@inline diffusivity(closure::IsotropicDiffusivity, ::Val{tracer_index}, args...) where tracer_index = closure.κ[tracer_index]
@inline viscosity(closure::IsotropicDiffusivity, args...) = closure.ν

Base.show(io::IO, closure::IsotropicDiffusivity) = print(io, "IsotropicDiffusivity: ν=$(closure.ν), κ=$(closure.κ)")


text:
begin
    using FileIO
    import Makie
    using ColorTypes
    using RegionTrees
    using AdaptiveDistanceFields
    using DistMesh
    using Zygote
    using Roots
    using GeometryBasics
    using LinearAlgebra
    using StaticArrays
    #using FractalGeometry
    using ForwardDiff
    using DifferentialEquations
    using SparseArrays
    using Plots
    import GeometryTypes
    using NearestNeighbors
end

begin
    include("geom_example.jl")
    include("simplicialsurface.jl")
end

begin
    testp = Point3{Float64}(1., -1., 0.64)
    const tau = 2pi
    const τ = tau
    const TAU = τ

    ori = Point3{Float64}(-1.1, -1.1, -1.1)
    wid = Point3{Float64}(2.2, 2.2, 2.2)
    orisv = SVector(ori)
    widsv = SVector(wid)
    #ori = GeometryBasics.Point{3, Float64}(-1.5, -1.5, -1.5)
    #wid = GeometryBasics.Point{3, Float64}(3.0, 3.0, 3.0)

    #result = distmesh(mandeldist, HUniform(), 0.04, DistMeshSetup(); origin=ori, widths=wid, maxiters=100)
    #result = distmesh(spheredist, HUniform(), 0.1, DistMeshSetup(); origin=ori, widths=wid, maxiters=100)

    outerbox(v) = boxdist(v, Point3{Float64}(0.8, 0.8, 0.8))

    boxdistcur(v) = boxdist(v, Point3{Float64}(0.2, 0.2, 0.2))
    mandeldist = p->-mandelbulbdistfunc(1)(p)
end

begin
    curve = Array{Tuple{Point3{Float64}, Float64}, 1}(undef, 0)
    push!(curve, (Point3{Float64}(-0.4, -0.4, -0.4), 1.0))
    push!(curve, (Point3{Float64}(0.4, -0.4, -0.4), 0.7))
    push!(curve, (Point3{Float64}(0.4, 0.4, -0.4), 0.45))
    push!(curve, (Point3{Float64}(0.4, 0.4, 0.4), 0.3))
    multibox = curveofshapes(curve, boxdistcur)
end

begin
    linecurveextrudedist(v) = curvedistgen(linecurve, -10., 10.)(v) - 0.4

    linecurveextrude_intersect_box_dist(v) = max(outerbox(v), linecurveextrudedist(v))
end

#=
function centraldiff(f::Function,p::VT) where VT
    println("Using our central diff")
    @show f, VT
    deps = sqrt(eps(eltype(VT)))
    dx = (f(p.+VT(deps,0,0)) - f(p.-VT(deps,0,0)))
    dy = (f(p.+VT(0,deps,0)) - f(p.-VT(0,deps,0)))
    dz = (f(p.+VT(0,0,deps)) - f(p.-VT(0,0,deps)))
    grad = VT(dx,dy,dz)./(2deps) #normalize?

end
=#

begin
    helixcurve(t) = Point3(0.3*cos(tau * t), 0.3*sin(tau * t), t)
    #helixcurve(t) = Point3{Float64}(0., 0., 1) .* t
    #helixcurvedist = curvedistgen(helixcurve, -10., 10.)
    #disttocurve(t) = normsq(helixcurve(t) - testp)
    #Ddisttocurve = D(disttocurve)
    helixcurveextrudedist(v) = curvedistgen(helixcurve, -10., 10.)(v) - 0.2
    helixcurveextrude_intersect_box_dist(v) = max(outerbox(v), helixcurveextrudedist(v))
end

# v expensive!!!
#helixadf = AdaptiveDistanceField(helixcurveextrude_intersect_box_dist, orisv, widsv)

begin
    weirdcurve(t) = Point3(0.5*cos(tau * t), 0.5*sin(tau * t), 0.4*cos(2tau*t))
    weirdcurveextrudedist(v) = curvedistgen(weirdcurve, -0.51, 0.51)(v) - 0.125
    weirdcurveextrude_intersect_box_dist(v) = max(outerbox(v), weirdcurveextrudedist(v))
end

begin
    arc_length_sampling_freq=0.01
    linecurvedist, γsampled_points = curvedistgenkdtree(weirdcurve, -0.50, 0.50; arc_length_sampling_freq=arc_length_sampling_freq)
    linecurveextrudedist(v) = linecurvedist(v) - 0.125
    linecurveextrude_intersect_box_dist(v) = max(outerbox(v), linecurveextrudedist(v))
end

begin
    curves = [linecurve]
    validtimes = [(-1., 1.)]
    extrusions = [-0.125]
    uniondistfuncs = [curvedistgenkdtree(curve, tmin, tmax; arc_length_sampling_freq=arc_length_sampling_freq)[1] for (curve, (tmin, tmax)) in zip(curves, validtimes)]
    uniondist(p) = minimum([uniondistfunc(p) - extrusion for (uniondistfunc, extrusion) in zip(uniondistfuncs, extrusions)])

end

begin
    curdist = uniondist
    curh0 = 0.1

    fh(v) = 0.05 + (0.5 * curdist(v))
    #fh = HUniform()
    result = distmesh(curdist, HUniform(), curh0, DistMeshSetup(); origin=ori, widths=wid, maxiters=100)
    resultouterbox = distmesh(outerbox, HUniform(), 0.5, DistMeshSetup(); origin=ori, widths=wid, maxiters=100)
    #resultouterbox = distmesh(outerbox, HUniform(), 5curh0, DistMeshSetup(); origin=ori, widths=wid, maxiters=2)

    #Makie.scatter(result.points)

    boundarytris, alltris = boundary_triangles(result.tetrahedra)
    boundarytrisouterbox, alltrisouterbox = boundary_triangles(resultouterbox.tetrahedra)

    boundarymesh = GeometryBasics.Mesh(result.points, boundarytris)
    boundarymeshouterbox = GeometryBasics.Mesh(resultouterbox.points, boundarytrisouterbox)
    #vertexsets = connectedvertices(faces(boundarytris))
    #removedmesh = removesmallcomponents(boundarymesh, vertexsets; sizetokeep=2)
    #allmesh = GeometryBasics.Mesh(result.points, alltris)
end
begin

    show_axis = true

    scene = Makie.Scene()
    Makie.mesh!(scene, boundarymesh, color=[norm(v[1:2]) for v in coordinates(boundarymesh)], show_axis=show_axis)

    Makie.wireframe!(scene, boundarymesh, show_axis=show_axis)
    Makie.wireframe!(scene, boundarymeshouterbox, show_axis=show_axis, color=RGB(0.8, 0.3, 0.8))
end

begin
    #scene = Makie.Scene()
    stride = 1
    Makie.scatter!(scene, γsampled_points[1:stride:end], show_axis=show_axis, color=[norm(point) for point in γsampled_points[1:stride:end]])
end

#Makie.wireframe!(scene, allmesh)
#Makie.wireframe!(scene, removedmesh)
#Makie.mesh!(scene, removedmesh, color=[norm(v[1:2]) for v in coordinates(removedmesh)])
#=
begin
aleaf = nothing
outv = nothing
for leaf in allleaves(adf.root)
    aleaf = leaf
    outv = RegionTrees.center(leaf)
    #outv = hcat(collect(RegionTrees.vertices(leaf.boundary))...)

end
end

adfpoints = map(RegionTrees.center, allleaves(helixadf.root))
=#

begin
    #savename = "meshes/generated/superhighreshelix1.stl"
    #save(savename, boundarymesh)
end
begin
    mesh = load(savename)
    scene = Makie.Scene()
    Makie.mesh!(scene, mesh, color=[norm(v[1:2]) for v in coordinates(mesh)])
    Makie.wireframe!(scene, mesh)
end

max_stars_repo_name: zoemcc/geometry

max_issues_repo_name: zoemcc/geometry

max_forks_repo_name: zoemcc/geometry


text:
using CubicHermiteSpline
using Test

@testset "basis function" begin
    h00, h10, h01, h11 = CubicHermiteSpline.basis(0)
    @test h00 == 1
    @test h10 == 0
    @test h01 == 0
    @test h11 == 0

    h00, h10, h01, h11 = CubicHermiteSpline.basis(1)
    @test h00 == 0
    @test h10 == 0
    @test h01 == 1
    @test h11 == 0

    t = 0.5
    h00, h10, h01, h11 = CubicHermiteSpline.basis(t)
    @test h00 ≈ (1 + 2t) * (t - 1)^2
    @test h10 ≈ t * (t - 1)^2
    @test h01 ≈ t^2 * (3 - 2t)
    @test h11 ≈ t^2 * (t - 1)
end

@testset "basis_derivative function" begin
    h00, h10, h01, h11 = CubicHermiteSpline.basis_derivative(0)
    @test h00 == 0
    @test h10 == 1
    @test h01 == 0
    @test h11 == 0

    h00, h10, h01, h11 = CubicHermiteSpline.basis_derivative(1)
    @test h00 == 0
    @test h10 == 0
    @test h01 == 0
    @test h11 == 1

    t = 0.5
    h00, h10, h01, h11 = CubicHermiteSpline.basis_derivative(t)
    @test h00 ≈ 6t * (t - 1)
    @test h10 ≈ (3t - 1) * (t - 1)
    @test h01 ≈ -6t * (t - 1)
    @test h11 ≈ t * (3t - 2)
end

@testset "Interval enclosing the data point" begin
    x = [0, 3, 5, 8, 10, 100]
    idx, x1, x2 = CubicHermiteSpline.findinterval(5, x)
    @test idx == 3
    @test x1 == 5
    @test x2 === nothing
    idx, x1, x2 = CubicHermiteSpline.findinterval(11, x)
    @test idx == 5
    @test x1 == 10
    @test x2 == 100

    x = [0.0, 0.3, 0.5, 0.8, 10.5, 100.1]
    idx, x1, x2 = CubicHermiteSpline.findinterval(0.5, x)
    @test idx == 3
    @test x1 == 0.5
    @test x2 === nothing
    idx, x1, x2 = CubicHermiteSpline.findinterval(0.99, x)
    @test idx == 4
    @test x1 == 0.8
    @test x2 == 10.5
end

@testset "Interpolation" begin
    # Cubic Hermite spline interpolation should exactly interpolate the 3rd order polynomial.
    f(x) = x^3 - 3x^2 + 2x - 5
    df(x) = 3x^2 - 6x + 2
    x = range(0, 2.5, step=0.5)
    y = f.(x)
    gradient = df.(x)
    spl = CubicHermiteSplineInterpolation(x, y, gradient)

    # Interpolate at the input data point, simply return input data
    xi = 0.5
    @test spl(xi) ≈ f(xi)
    @test interp(spl, xi) ≈ f(xi)
    # Interpolate at arbitrary position
    xi = 1.2
    @test spl(xi) ≈ f(xi)
    @test interp(spl, xi) ≈ f(xi)

    # Test array interpolation
    xi = [0.5, 1.2]
    yi = spl(xi)
    @test yi[1] ≈ f(xi[1])
    @test yi[2] ≈ f(xi[2])
    yi = interp(spl, xi)
    @test yi[1] ≈ f(xi[1])
    @test yi[2] ≈ f(xi[2])
end

@testset "Interpolation grad" begin
    # Cubic Hermite spline interpolation should exactly interpolate the 3rd order polynomial.
    f(x) = x^3 - 3x^2 + 2x - 5
    df(x) = 3x^2 - 6x + 2
    # x = range(0, 2.5, step=0.5)
    x = [0, 0.3, 0.8, 1.5, 2.7, 3]
    y = f.(x)
    gradient = df.(x)
    spl = CubicHermiteSplineInterpolation(x, y, gradient)

    # Interpolate at the input data point, simply return input data
    xi = 0.5
    @test spl(xi; grad=true) ≈ df(xi)
    @test grad(spl, xi) ≈ df(xi)
    # Interpolate at arbitrary position
    xi = 2.1
    @test spl(xi; grad=true) ≈ df(xi)
    @test grad(spl, xi) ≈ df(xi)

    # Test array interpolation
    xi = [0.5, 2.1]
    yi = spl(xi; grad=true)
    @test yi[1] ≈ df(xi[1])
    @test yi[2] ≈ df(xi[2])
    yi = grad(spl, xi)
    @test yi[1] ≈ df(xi[1])
    @test yi[2] ≈ df(xi[2])
end


text:
### A Pluto.jl notebook ###
# v0.17.2

# using Markdown
# using InteractiveUtils

# ╔═╡ f9eff180-3c4d-49c2-9f4e-e3f425a96966
#=╠═╡ notebook_exclusive
begin
	import Pkg
	Pkg.activate("..")
end
  ╠═╡ notebook_exclusive =#

# ╔═╡ 941fa912-61d4-4847-8bb2-5e83eb31ca34
begin
	using Revise
	using PlutoDevMacros
end

# ╔═╡ f73c9b70-4c68-4d3b-be3e-067a118681e8
#=╠═╡ notebook_exclusive
using BenchmarkTools
  ╠═╡ notebook_exclusive =#

# ╔═╡ d9b0ef51-60cc-4a30-9f8b-f6b6b2dbeec5
@macroexpand @plutoinclude "test1.jl" "all"

# ╔═╡ e2ac5a43-d683-486b-a2a1-f26f48313c47
@plutoinclude "test1.jl" "all"

# ╔═╡ f3287e4b-c92a-4eef-8b79-78346a55803b
asd(x::Float64) = "FLOAT"

# ╔═╡ 4d8e68cb-910f-4f6c-89cd-2bf2da7b5a70
struct TestStruct2 end

# ╔═╡ 53ec02f7-ca14-41d9-b3b7-757c36f526f7
asd(::TestStruct2) = "TESTSTRUCT2"

# ╔═╡ f5b8ccdd-ad66-4ae4-b3e9-82452464168b
#=╠═╡ notebook_exclusive
asd(3)
  ╠═╡ notebook_exclusive =#

# ╔═╡ a33fdee5-037c-41de-aeec-bc4383fe4826
#=╠═╡ notebook_exclusive
@benchmark asd(TestStruct2())
  ╠═╡ notebook_exclusive =#

# ╔═╡ a646c602-529e-45c1-a170-3650318c1c2d
#=╠═╡ notebook_exclusive
@benchmark asd(TestStruct1())
  ╠═╡ notebook_exclusive =#

# ╔═╡ Cell order:
# ╠═f9eff180-3c4d-49c2-9f4e-e3f425a96966
# ╠═941fa912-61d4-4847-8bb2-5e83eb31ca34
# ╠═f73c9b70-4c68-4d3b-be3e-067a118681e8
# ╠═d9b0ef51-60cc-4a30-9f8b-f6b6b2dbeec5
# ╠═e2ac5a43-d683-486b-a2a1-f26f48313c47
# ╠═f5b8ccdd-ad66-4ae4-b3e9-82452464168b
# ╠═f3287e4b-c92a-4eef-8b79-78346a55803b
# ╠═4d8e68cb-910f-4f6c-89cd-2bf2da7b5a70
# ╠═53ec02f7-ca14-41d9-b3b7-757c36f526f7
# ╠═a33fdee5-037c-41de-aeec-bc4383fe4826
# ╠═a646c602-529e-45c1-a170-3650318c1c2d


text:
using Distributed
using JLD2

@time @everywhere include("../BioClim/src/required.jl")

## Get & prepare data
@time @everywhere begin
    # Load data from CSV files
    df = CSV.read("/Users/kiristern/Documents/GitHub/SDM/data/pred_prey/Puma_concolor.csv", header=true)
    # Prepare data (select columns, arrange values)
    df = prepare_gbif_data(df)
    # Separate species
    taxa_occ = [df[df.species .== u,:] for u in unique(df.species)]
    # # Define coordinates range
    # lon_range = (-165.0, -30.0)
    # lat_range = (75.0, -60.0)
end

## Load predictions
@load "data/predictions.jld2" predictions
pred = predictions[1]
prey1 = predictions[2]
prey2 = predictions[3]

## expand
p1 = pred
p2 = prey1
p3 = prey2

# Get coordinates of predictions
#pred
lons_p1 = longitudes(p1)
lats_p1 = latitudes(p1)
#prey1
lons_p2 = longitudes(p2)
lats_p2 = latitudes(p2)
#prey2
lons_p3 = longitudes(p3)
lats_p3 = latitudes(p3)

# Get position of p2 on p1
m_lon = findmin(abs.(minimum(lons_p2) .- lons_p1))[2]
M_lon = findmin(abs.(maximum(lons_p2) .- lons_p1))[2]
m_lat = findmin(abs.(minimum(lats_p2) .- lats_p1))[2]
M_lat = findmin(abs.(maximum(lats_p2) .- lats_p1))[2]
p1[(m_lat:M_lat), (m_lon:M_lon)]
# Check if size matches
size(p2)
size(p3)

# Create new layer
new1 = fill(NaN, size(p1))
new2 = fill(NaN, size(p1))
new1[(m_lat:M_lat), (m_lon:M_lon)] .= p1[(m_lat:M_lat), (m_lon:M_lon)]
new2[(m_lat:M_lat), (m_lon:M_lon)] .= p2.grid
newlayer1 = SDMLayer(new1, p1.left, p1.right, p1.bottom, p1.top)
newlayer2 = SDMLayer(new2, p1.left, p1.right, p1.bottom, p1.top)

# Plot result
plotp1 = plotSDM(p1)
plotp2 = plotSDM(p2)
plotn1 = plotSDM(newlayer1)
plotn2 = plotSDM(newlayer2)

# Combine heatmaps
sdm_plot = plotSDM(newlayer1)
heatmap!(
    sdm_plot,
    longitudes(newlayer2), latitudes(newlayer2), # layer range
    newlayer2.grid, # evenness values
    aspectratio=92.60/60.75, # aspect ratio
    c=:viridis, # ~color palette
    clim=(0.0, maximum(filter(!isnan, newlayer2.grid))) # colorbar limits
)
plot!(
    sdm_plot,
    xlims=(p2.left, p2.right),
    ylims=(p2.bottom, p2.top),
    aspectratio=92.60/60.75
)
# Compare to p2
plotp2

max_stars_repo_name: kiristern/SDM

max_issues_repo_name: kiristern/SDM

max_forks_repo_name: kiristern/SDM


text:
#==========================================================
This routine initializes an optimization problem to recover
a |0⟩ to |2⟩ swap gate on a single qudit with 3 energy
levels (and 1 guard state). The  drift Hamiltonian in the
rotating frame is
        H0 = 2π*diagm([0 0 -2.2538e-1 -7.0425e-1]).
Here the control Hamiltonian includes the usual symmetric
and anti-symmetric terms
     H_{sym} = p(t)(a + a^†),    H_{asym} = q(t)(a - a^†),
where a is the annihilation operator for the qudit.
The problem parameters for this example are:
                ω_a =  2π × 4.09947    Grad/s,
                ξ_a =  2π × 2.2538e-01 Grad/s.
We use Bsplines with carrier waves with frequencies
0, ξ_a Grad/s.
==========================================================#
using LinearAlgebra
#using Plots
#pyplot()
#using FFTW
#using DelimitedFiles
using Printf
#using Ipopt
using Random

Base.show(io::IO, f::Float64) = @printf(io, "%20.13e", f)

import Juqbox

verbose = false
Nosc=1 # Number of oscillators

N = 3 # 4 # Number of essential energy levels
Nguard = 1 # 0 # Number of guard/forbidden energy levels
Ntot = N + Nguard # Total number of energy levels

samplerate = 32 # for output files
casename = "swap02" # naming output files"

T = 150.0 # 150.0 # 200.0 # 100.0 # Duration of gate

# frequencies (in GHz, will be multiplied by 2*pi to get angular frequencies in the Hamiltonian matrix)
freq_alice=[0, 4.09947, 3.87409, 3.6206] # GHz
fa = freq_alice[2]
@assert(Ntot <= 4) # we don't know any higher frequencies

utarget = zeros(ComplexF64,Ntot,N)

# 0>  to  |2>  swap gate
if N >= 3
    utarget[1,1] = 0 #1/sqrt(2)
    utarget[2,1] = 0
    utarget[3,1] = 1 #1/sqrt(2)
#
    utarget[1,2] = 0
    utarget[2,2] = 1
    utarget[3,2] = 0
    #
    utarget[1,3] = 1 #1/sqrt(2)
    utarget[2,3] = 0
    utarget[3,3] = 0 #-1/sqrt(2)
#
end

if N==4
    utarget[4,4] = 1
end

omega1 = Juqbox.setup_rotmatrices([N], [Nguard], [freq_alice[2]])

# Compute Ra*utarget
rot1 = Diagonal(exp.(im*omega1*T))

# target in the rotating frame
vtarget = rot1*utarget

startFromScratch = true # false
startFile="pcof-60-t100.dat"

bsname = "pcof-60-t100-new.dat"
use_bcarrier = true

# frequencies (in GHz, will be multiplied by 2*pi to get angular frequencies in the Hamiltonian matrix)
fa = 4.10595
xa = 2* 0.1099

number = Diagonal(collect(0:Ntot-1))

H0 = -0.5*(2*pi)*xa* (number*number - number) # xa is in GHz

# lowering matrix
amat = Bidiagonal(zeros(Ntot),sqrt.(collect(1:Ntot-1)),:U) # standard loweing matrix
adag = Array(transpose(amat));
Hsym_ops=[Array(amat + adag)]
Hanti_ops=[Array(amat - adag)]
H0 = Array(H0)

# END NEW SETUP

Ncoupled = length(Hsym_ops) # Number of paired Hamiltonians
Nfreq= 2 # number of carrier frequencies 3 gives a cleaner sol than 2

# setup carrier frequencies
om = zeros(Ncoupled,Nfreq)
# seg=1 is the real, seg=2 is the imaginary. Use the same frequencies for both
if Nfreq == 3
    om[1:Ncoupled,2] .= H0[3,3]
    om[1:Ncoupled,3] .= H0[4,4] - H0[3,3]
elseif Nfreq == 2
    om[1:Ncoupled,2] .= H0[3,3]
end
#println("Carrier frequencies [GHz]: ", om[1,:]./(2*pi))

#max amplitude (in angular frequency) 2*pi*GHz
# NOTE: optimize-once.jl uses maxpar to generate filenames
maxpar =2*pi*0.0132/Nfreq/2 # attempt to limit lab frame amplitude to 0.0132 = 0.5/6/2/pi

# Estimate time step
K1 =  H0 + maxpar.*( amat +  amat') + 1im* maxpar.*(amat - amat')
lamb = eigvals(Array(K1))
maxeig = maximum(abs.(lamb))

Pmin = 80
samplerate1 = maxeig*Pmin/(2*pi)
nsteps = ceil(Int64, T*samplerate1)
#println("Duration T = ", T, " number of time-steps = ", nsteps)

# Initial conditions
Ident = Matrix{Float64}(I, Ntot, Ntot)
U0 = Ident[1:Ntot,1:N]

# try the new Juqbox2 module
params = Juqbox.objparams([N], [Nguard], T, nsteps, Uinit=U0, Utarget=vtarget, Cfreq=om, Rfreq=[freq_alice[2]],
                          Hconst=H0, Hsym_ops=Hsym_ops, Hanti_ops=Hanti_ops)

# Quiet mode for testing
params.quiet = true

# set random number seed to make algorithm deterministic
Random.seed!(2456)

# initial parameter guess
if startFromScratch
  D1 = 10 # Number of B-spline coefficients per frequency, sin/cos and real/imag
  nCoeff = 2*Ncoupled*Nfreq*D1 # factor '2' is for sin/cos
  #pcof0  = zeros(nCoeff)
  pcof0 = (rand(nCoeff) .- 0.5).*maxpar*0.1
#  println("*** Starting from random pcof with amplitude ", maxpar*0.1)
else
  # use if you want to have initial coefficients read from file
  pcof0 = vec(readdlm(startFile))
    #  println("*** Starting from B-spline coefficients in file: ", startFile)
    nCoeff = length(pcof0)
    D1 = div(nCoeff, 2*Ncoupled*Nfreq)  # number of B-spline coeff per control function
end

# min and max coefficient values
useBarrier = true
minCoeff = -maxpar*ones(nCoeff);
maxCoeff = maxpar*ones(nCoeff);

if verbose
    println("*** Settings ***")
    println("Resonant frequencies [GHz] = ", freq_alice[1:Ntot])
    println("Total number of states, Ntot = ", Ntot, " Total number of guard states, Nguard = ", Nguard)
    if use_bcarrier
        println("Using B-spline basis functions with carrier wave, # freq = ", Nfreq)
    else
        println("Using regular B-spline basis functions")
    end
    println("Number of coefficients per spline = ", D1, " Total number of parameters = ", nCoeff)
    println("Max parameter amplitudes: maxpar = ", maxpar)
    println("Tikhonov coefficients: tik0 = ", params.tik0)
end

# optional arguments to setup_ipopt_problem()
maxIter = 50
lbfgsMax = 250

# Estimate number of terms in Neumann series for time stepping (Default 3)
tol = eps(1.0); # machine precision
Juqbox.estimate_Neumann!(tol, params, [maxpar])

# Allocate all working arrays
wa = Juqbox.Working_Arrays(params, nCoeff)
prob = Juqbox.setup_ipopt_problem(params, wa, nCoeff, minCoeff, maxCoeff, maxIter=maxIter, lbfgsMax=lbfgsMax)

# uncomment to run the gradient checker for the initial pcof
# addOption( prob, "derivative_test", "first-order"); # for testing the gradient

#println("Initial coefficient vector stored in 'pcof0'")


text:
######################################################################
# max.jl
######################################################################
using ErgodicControl

# Ergodic Manager
d = Domain([1,1], 100)
means = [[.2,.5], [.8,.5]]
Sigmas = [.015*eye(2), .015*eye(2)]
weights = [1., 2.]
phi = gaussian(d, means, Sigmas, weights)
em = ErgodicManagerR2(d, phi, 5)

# Trajectory Manager
x0 = [0.49,0.01]
x0 = [0.49,0.01, 0, 0]
N = 30
h = 0.5
ci = ConstantInitializer([0.000, 0.000])
tm = TrajectoryManager(x0, h, N, ci)
dynamics!(tm, DoubleIntegrator(2,h))

# max trajectory
mi = 1000
tm.R = 50eye(2)
xd, ud = max_trajectory(em, tm, means, Sigmas, weights, max_iters=mi)


text:
"""
Generates various Genie files.
"""
module Generator

import SHA, Dates, Pkg, Logging, UUIDs
import Inflector
import Genie

const JULIA_PATH = joinpath(Sys.BINDIR, "julia")

function validname(name::String)
  filter(! isempty, [x.match for x in collect(eachmatch(r"[0-9a-zA-Z_]*", name))]) |> join
end

"""
    newcontroller(resource_name::String) :: Nothing

Generates a new Genie controller file and persists it to the resources folder.
"""
function newcontroller(resource_name::String; path::Union{String,Nothing} = nothing, pluralize::Bool = true) :: Nothing
  resource_name = validname(resource_name)

  Inflector.is_singular(resource_name) && pluralize && (resource_name = Inflector.to_plural(resource_name))
  resource_name = uppercasefirst(resource_name)

  resource_path = path === nothing ? setup_resource_path(resource_name, path = ".") : (ispath(path) ? path : mkpath(path))
  cfn = controller_file_name(resource_name)
  write_resource_file(resource_path, cfn, resource_name, :controller, pluralize = pluralize) &&
    @info "New controller created at $(abspath(joinpath(resource_path, cfn)))"

  nothing
end

"""
    newresource(resource_name::String, config::Settings) :: Nothing

Generates all the files associated with a new resource and persists them to the resources folder.
"""
function newresource(resource_name::String; path::String = ".", pluralize::Bool = true) :: Nothing
  resource_name = validname(resource_name)

  Inflector.is_singular(resource_name) && pluralize &&
    (resource_name = Inflector.to_plural(resource_name))

  resource_path = setup_resource_path(resource_name, path = path)
  for (resource_file, resource_type) in [(controller_file_name(resource_name), :controller)]
    write_resource_file(resource_path, resource_file, resource_name, resource_type, pluralize = pluralize) &&
      @info "New $resource_file created at $(abspath(joinpath(resource_path, resource_file)))"
  end

  views_path = joinpath(resource_path, "views")
  isdir(views_path) || mkpath(views_path)

  nothing
end

"""
    setup_resource_path(resource_name::String) :: String

Computes and creates the directories structure needed to persist a new resource.
"""
function setup_resource_path(resource_name::String; path::String = ".") :: String
  isdir(Genie.config.path_app) || Genie.Generator.mvc_support(path)

  resource_path = joinpath(path, Genie.config.path_resources, lowercase(resource_name))

  if ! isdir(resource_path)
    mkpath(resource_path)
    push!(LOAD_PATH, resource_path)
  end

  resource_path
end

"""
    write_resource_file(resource_path::String, file_name::String, resource_name::String) :: Bool

Generates all resouce files and persists them to disk.
"""
function write_resource_file(resource_path::String, file_name::String, resource_name::String, resource_type::Symbol; pluralize::Bool = true) :: Bool
  resource_name = (pluralize ? (Inflector.to_plural(resource_name)) : resource_name) |> Inflector.from_underscores

  try
    if resource_type == :controller
      resource_does_not_exist(resource_path, file_name) || return true
      open(joinpath(resource_path, file_name), "w") do f
        write(f, Genie.FileTemplates.newcontroller(resource_name))
      end
    end
  catch ex
    @error ex
  end

  try
    if resource_type == :test
      resource_does_not_exist(resource_path, file_name) || return true
      open(joinpath(resource_path, file_name), "w") do f
        name = pluralize ? (Inflector.to_singular(resource_name)) : resource_name
        write(f, Genie.FileTemplates.newtest(resource_name,  name))
      end
    end
  catch ex
    @error ex
  end

  try
    Genie.load_resources()
  catch ex
    @error ex
  end

  true
end

function binfolderpath(path::String) :: String
  bin_folder_path = joinpath(path, Genie.config.path_bin)
  isdir(bin_folder_path) || mkpath(bin_folder_path)

  bin_folder_path
end

"""
    setup_windows_bin_files(path::String = ".") :: Nothing

Creates the bin/server and bin/repl binaries for Windows
"""
function setup_windows_bin_files(path::String = ".") :: Nothing
  bin_folder_path = binfolderpath(path)

  open(joinpath(bin_folder_path, "repl.bat"), "w") do f
    write(f, "\"$JULIA_PATH\" --color=yes --depwarn=no --project=@. -q -i -- \"%~dp0..\\$(Genie.BOOTSTRAP_FILE_NAME)\" %*")
  end

  open(joinpath(bin_folder_path, "server.bat"), "w") do f
    write(f, "\"$JULIA_PATH\" --color=yes --depwarn=no --project=@. -q -i -- \"%~dp0..\\$(Genie.BOOTSTRAP_FILE_NAME)\" s %*")
  end

  open(joinpath(bin_folder_path, "runtask.bat"), "w") do f
    write(f, "\"$JULIA_PATH\" --color=yes --depwarn=no --project=@. -q -- \"%~dp0..\\$(Genie.BOOTSTRAP_FILE_NAME)\" -r %*")
  end

  nothing
end

"""
    setup_nix_bin_files(path::String = ".") :: Nothing

Creates the bin/server and bin/repl binaries for *nix systems
"""
function setup_nix_bin_files(path::String = ".") :: Nothing
  bin_folder_path = binfolderpath(path)

  open(joinpath(bin_folder_path, "repl"), "w") do f
    write(f, "#!/bin/sh\n" * raw"julia --color=yes --depwarn=no --project=@. -q -L $(dirname $0)/../bootstrap.jl -- \"$@\"")
  end

  open(joinpath(bin_folder_path, "server"), "w") do f
    write(f, "#!/bin/sh\n" * raw"julia --color=yes --depwarn=no --project=@. -q -i -- $(dirname $0)/../bootstrap.jl s \"$@\"")
  end

  open(joinpath(bin_folder_path, "runtask"), "w") do f
    write(f, "#!/bin/sh\n" * raw"julia --color=yes --depwarn=no --project=@. -q -- $(dirname $0)/../bootstrap.jl -r \"$@\"")
  end

  chmod(joinpath(bin_folder_path, "server"), 0o700)
  chmod(joinpath(bin_folder_path, "repl"), 0o700)
  chmod(joinpath(bin_folder_path, "runtask"), 0o700)

  nothing
end

"""
    resource_does_not_exist(resource_path::String, file_name::String) :: Bool

Returns `true` if the indicated resources does not exists - false otherwise.
"""
function resource_does_not_exist(resource_path::String, file_name::String) :: Bool
  if isfile(joinpath(resource_path, file_name))
    @warn "File already exists, $(joinpath(resource_path, file_name)) - skipping"

    return false
  end

  true
end

"""
    controller_file_name(resource_name::Union{String,Symbol})

Computes the controller file name based on the resource name.
"""
function controller_file_name(resource_name::Union{String,Symbol}) :: String
  uppercasefirst(string(resource_name)) * Genie.GENIE_CONTROLLER_FILE_POSTFIX
end

"""
    secret_token() :: String

Generates a random secret token to be used for configuring the call to `Genie.secret_token!`.
"""
function secret_token() :: String
  SHA.sha256("$(randn()) $(Dates.now())") |> bytes2hex
end

"""
    write_secrets_file(app_path=".")

Generates a valid `config/secrets.jl` file with a random secret token.
"""
function write_secrets_file(app_path::String = ".") :: Nothing
  secrets_path = joinpath(app_path, Genie.config.path_config)
  ispath(secrets_path) || mkpath(secrets_path)

  open(joinpath(secrets_path, Genie.SECRETS_FILE_NAME), "w") do f
    write(f, """Genie.secret_token!("$(secret_token())") """)
  end

  nothing
end

"""
    migrate_secrets_file(app_path=".")

Replace the previous content of `config/secrets.jl` with the new syntax but the same token.

Specifically, the previous syntax of the file is the following, which is now deprecated:

    const SECRET_TOKEN = "token"

This syntax is replaced by the new syntax:

    Genie.secret_token!("token")
"""
function migrate_secrets_file(app_path::String = ".") :: Nothing
  secrets_path = joinpath(app_path, Genie.config.path_config, Genie.SECRETS_FILE_NAME)
  if isfile(secrets_path)
    match_deprecated = match(r"SECRET_TOKEN\s*=\s*\"(.*)\"", readline(secrets_path))
    if match_deprecated != nothing # does the file use the deprecated syntax?
      open(secrets_path, "w") do f
        write(f, """Genie.secret_token!("$(match_deprecated.captures[1])") """)
      end # replace the previous content of the file
      @info "Successfully migrated $(secrets_path) to a valid syntax."
    else
      error("No migration possible; $(secrets_path) is not using a migrate-able syntax.")
    end
  else
    error("No migration possible; $(secrets_path) is not a file.")
  end
  return nothing
end

"""
    fullstack_app(app_path::String = ".") :: Nothing

Writes the files necessary to create a full stack Genie app.
"""
function fullstack_app(app_name::String = ".", app_path::String = ".") :: Nothing
  cp(joinpath(@__DIR__, "..", Genie.NEW_APP_PATH), app_path)

  scaffold(app_name, app_path)

  nothing
end

"""
    new(app_name::String, app_path::String = abspath(app_name), autostart::Bool = true) :: Nothing

Creates a minimal Genie app.
"""
function new(app_name::String, app_path::String = "", autostart::Bool = true) :: Nothing
  app_name = validname(app_name)
  app_path = abspath(app_name)

  scaffold(app_name, app_path)

  post_create(app_name, app_path, autostart = autostart)

  nothing
end

"""
    scaffold(app_path::String = ".") :: Nothing

Writes the file necessary to scaffold a minimal Genie app.
"""
function scaffold(app_name::String, app_path::String = "") :: Nothing
  app_name = validname(app_name)
  app_path = abspath(app_name)

  isdir(app_path) || mkpath(app_path)

  for f in [Genie.config.path_src, Genie.GENIE_FILE_NAME, Genie.ROUTES_FILE_NAME,
            ".gitattributes", ".gitignore"]
    try
      cp(joinpath(@__DIR__, "..", Genie.NEW_APP_PATH, f), joinpath(app_path, f))
    catch ex
    end
  end

  write_app_custom_files(app_name, app_path)

  nothing
end

"""
    microstack_app(app_path::String = ".") :: Nothing

Writes the file necessary to create a microstack app.
"""
function microstack_app(app_name::String = ".", app_path::String = ".") :: Nothing
  isdir(app_path) || mkpath(app_path)

  for f in [Genie.config.path_bin, Genie.config.path_config, Genie.config.server_document_root]
    cp(joinpath(@__DIR__, "..", Genie.NEW_APP_PATH, f), joinpath(app_path, f))
  end

  scaffold(app_name, app_path)

  nothing
end

"""
    mvc_support(app_path::String = ".") :: Nothing

Writes the files used for rendering resources using the MVC stack and the Genie templating system.
"""
function mvc_support(app_path::String = ".") :: Nothing
  cp(joinpath(@__DIR__, "..", Genie.NEW_APP_PATH, Genie.config.path_app), joinpath(app_path, Genie.config.path_app))

  nothing
end

"""
    db_support(app_path::String = ".") :: Nothing

Writes files used for interacting with the SearchLight ORM.
"""
function db_support(app_path::String = ".", include_env::Bool = true, add_dependencies::Bool = true;
                    testmode::Bool = false) :: Nothing
  cp(joinpath(@__DIR__, "..", Genie.NEW_APP_PATH, Genie.config.path_db), joinpath(app_path, Genie.config.path_db), force = true)

  initializers_dir = joinpath(app_path, Genie.config.path_initializers)
  initializer_path = joinpath(initializers_dir, Genie.SEARCHLIGHT_INITIALIZER_FILE_NAME)
  source_path = joinpath(@__DIR__, "..", Genie.NEW_APP_PATH, Genie.config.path_initializers, Genie.SEARCHLIGHT_INITIALIZER_FILE_NAME) |> normpath

  if !isfile(initializer_path)
    ispath(initializers_dir) || mkpath(initializers_dir)
    include_env && cp(joinpath(@__DIR__, "..", Genie.NEW_APP_PATH, Genie.config.path_env), joinpath(app_path, Genie.config.path_env), force = true)
    cp(source_path, initializer_path)
  end

  add_dependencies && install_db_dependencies(testmode = testmode)

  nothing
end

"""
    write_app_custom_files(path::String, app_path::String) :: Nothing

Writes the Genie app main module file.
"""
function write_app_custom_files(path::String, app_path::String) :: Nothing
  moduleinfo = Genie.FileTemplates.appmodule(path)

  open(joinpath(app_path, Genie.config.path_src, moduleinfo[1] * ".jl"), "w") do f
    write(f, moduleinfo[2])
  end

  open(joinpath(app_path, Genie.BOOTSTRAP_FILE_NAME), "w") do f
    write(f,
    """
    pwd() == joinpath(@__DIR__, "bin") && cd(@__DIR__) # allow starting app from bin/ dir

    using $(moduleinfo[1])
    $(moduleinfo[1]).main()
    """)
  end

  isdir(joinpath(app_path, "test")) || mkpath(joinpath(app_path, "test"))
  open(joinpath(app_path, "test", "runtests.jl"), "w") do f
    write(f,
      """
      using $(moduleinfo[1]), Test
      # implement your tests here
      @test 1 == 1
      """)
  end

  nothing
end

"""
    install_app_dependencies(app_path::String = ".") :: Nothing

Installs the application's dependencies using Julia's Pkg
"""
function install_app_dependencies(app_path::String = "."; testmode::Bool = false, dbsupport::Bool = false) :: Nothing
  @info "Installing app dependencies"
  Pkg.activate(".")

  pkgs = ["Dates", "Logging", "LoggingExtras", "MbedTLS"]

  testmode ? Pkg.develop("Genie") : push!(pkgs, "Genie")
  testmode ? Pkg.develop("Inflector") : push!(pkgs, "Inflector")

  Pkg.add(pkgs)

  dbsupport && install_db_dependencies(testmode = testmode)

  @info "Installing dependencies for unit tests"

  Pkg.activate("test")

  Pkg.add("Test")

  Pkg.activate(".") # return to the main project

  nothing
end

"""
    generate_project(name)

Generate the `Project.toml` with a name and a uuid.
If this file already exists, generate `Project_sample.toml` as a reference instead.
"""
function generate_project(name::String) :: Nothing
  name = Genie.FileTemplates.appmodule(name)[1] # convert to camel case

  mktempdir() do tmpdir
    tmp = joinpath(tmpdir, name, "Project.toml")

    pkgproject(Pkg.API.Context(), name, tmpdir) # generate tmp

    # Pkg.project(Pkg.stdout_f(), name, tmpdir) # generate tmp

    if !isfile("Project.toml")
      mv(tmp, "Project.toml") # move tmp here
      @info "Project.toml has been generated"
    else
      mv(tmp, "Project_sample.toml"; force = true)
      @warn "$(abspath("."))/Project.toml already exists and will not be replaced. " *
        "Make sure that it specifies a name and a uuid, using Project_sample.toml as a reference."
    end
  end # remove tmpdir on completion

  nothing
end

function pkgproject(ctx::Pkg.API.Context, pkg::String, dir::String) :: Nothing
  name = email = nothing

  gitname = Pkg.LibGit2.getconfig("user.name", "")
  isempty(gitname) || (name = gitname)

  gitmail = Pkg.LibGit2.getconfig("user.email", "")
  isempty(gitmail) || (email = gitmail)

  if name === nothing
      for env in ["GIT_AUTHOR_NAME", "GIT_COMMITTER_NAME", "USER", "USERNAME", "NAME"]
          name = get(ENV, env, nothing)
          name !== nothing && break
      end
  end

  name === nothing && (name = "Unknown")

  if email === nothing
      for env in ["GIT_AUTHOR_EMAIL", "GIT_COMMITTER_EMAIL", "EMAIL"];
          email = get(ENV, env, nothing)
          email !== nothing && break
      end
  end

  authors = ["$name " * (email === nothing ? "" : "<$email>")]

  uuid = UUIDs.uuid4()

  pkggenfile(ctx, pkg, dir, "Project.toml") do io
      toml = Dict{String,Any}("authors" => authors,
                              "name" => pkg,
                              "uuid" => string(uuid),
                              "version" => "0.1.0",
                              )
      Pkg.TOML.print(io, toml, sorted=true, by=key -> (Pkg.Types.project_key_order(key), key))
  end
end

function pkggenfile(f::Function, ctx::Pkg.API.Context, pkg::String, dir::String, file::String) :: Nothing
  path = joinpath(dir, pkg, file)
  println(ctx.io, "    $(Base.contractuser(path))")
  mkpath(dirname(path))
  open(f, path, "w")
end

function install_db_dependencies(; testmode::Bool = false) :: Nothing
  try
    Pkg.add("SearchLight")
    testmode || install_searchlight_dependencies()
  catch ex
    @error ex
  end

  nothing
end

function install_searchlight_dependencies() :: Nothing # TODO: move this to SearchLight post install
  backends = ["SQLite", "MySQL", "PostgreSQL"]

  println("Please choose the DB backend you want to use: ")
  for i in 1:length(backends)
    println("$i. $(backends[i])")
  end
  println("Input $(join([1:length(backends)...], ", ", " or ")) and press ENTER to confirm")
  println()

  choice = try
    parse(Int, readline())
  catch _
    0
  end

  (choice in [1, 2, 3]) || return install_searchlight_dependencies()

  Pkg.add("SearchLight$(backends[choice])")

  nothing
end

"""
    autostart_app(path::String = "."; autostart::Bool = true) :: Nothing

If `autostart` is `true`, the newly generated Genie app will be automatically started.
"""
function autostart_app(path::String = "."; autostart::Bool = true) :: Nothing
  if autostart
    @info "Starting your brand new Genie app - hang tight!"
    Genie.loadapp(abspath(path), autostart = autostart)
  else
    @info("Your new Genie app is ready!
        Run \njulia> Genie.loadapp() \nto load the app's environment
        and then \njulia> up() \nto start the web server on port 8000.")
  end

  nothing
end

"""
    remove_searchlight_initializer(app_path::String = ".") :: Nothing

Removes the SearchLight initializer file if it's unused
"""
function remove_searchlight_initializer(app_path::String = ".") :: Nothing
  rm(joinpath(app_path, Genie.config.path_initializers, Genie.SEARCHLIGHT_INITIALIZER_FILE_NAME), force = true)

  nothing
end

"""
    newapp(app_name::String; autostart::Bool = true, fullstack::Bool = false, dbsupport::Bool = false, mvcsupport::Bool = false) :: Nothing

Scaffolds a new Genie app, setting up the file structure indicated by the various arguments.

# Arguments
- `path::String`: the name of the app and the path where to bootstrap it. Spaces not allowed
- `autostart::Bool`: automatically start the app once the file structure is created
- `fullstack::Bool`: the type of app to be bootstrapped. The fullstack app includes MVC structure, DB connection code, and asset pipeline files.
- `dbsupport::Bool`: bootstrap the files needed for DB connection setup via the SearchLight ORM
- `mvcsupport::Bool`: adds the files used for HTML+Julia view templates rendering and working with resources

# Examples
```julia-repl
julia> Genie.newapp("MyGenieApp")
2019-08-06 16:54:15:INFO:Main: Done! New app created at MyGenieApp
2019-08-06 16:54:15:DEBUG:Main: Changing active directory to MyGenieApp
2019-08-06 16:54:15:DEBUG:Main: Installing app dependencies
 Resolving package versions...
  Updating `~/Dropbox/Projects/GenieTests/MyGenieApp/Project.toml`
  [c43c736e] + Genie v0.10.1
  Updating `~/Dropbox/Projects/GenieTests/MyGenieApp/Manifest.toml`

2019-08-06 16:54:27:INFO:Main: Starting your brand new Genie app - hang tight!
 _____         _
|   __|___ ___|_|___
|  |  | -_|   | | -_|
|_____|___|_|_|_|___|

┌ Info:
│ Starting Genie in >> DEV << mode
└
[ Info: Logging to file at MyGenieApp/log/dev.log
[ Info: Ready!
2019-08-06 16:54:32:DEBUG:Main: Web Server starting at http://127.0.0.1:8000
2019-08-06 16:54:32:DEBUG:Main: Web Server running at http://127.0.0.1:8000
```
"""
function newapp(app_name::String; autostart::Bool = true, fullstack::Bool = false,
                dbsupport::Bool = false, mvcsupport::Bool = false, testmode::Bool = false) :: Nothing
  app_name = validname(app_name)
  app_path = abspath(app_name)

  fullstack ? fullstack_app(app_name, app_path) : microstack_app(app_name, app_path)

  (dbsupport || fullstack) ? db_support(app_path, testmode = testmode) : remove_searchlight_initializer(app_path)

  mvcsupport && (fullstack || mvc_support(app_path))

  write_secrets_file(app_path)

  try
    setup_windows_bin_files(app_path)
  catch ex
    @error ex
  end

  try
    setup_nix_bin_files(app_path)
  catch ex
    @error ex
  end

  post_create(app_name, app_path; autostart = autostart, testmode = testmode, dbsupport = (dbsupport || fullstack))

  nothing
end

function post_create(app_name::String, app_path::String; autostart::Bool = true, testmode::Bool = false, dbsupport::Bool = false)
  @info "Done! New app created at $app_path"

  @info "Changing active directory to $app_path"
  cd(app_path)

  generate_project(app_name)

  install_app_dependencies(app_path, testmode = testmode, dbsupport = dbsupport)

  set_files_mod()

  autostart_app(app_path, autostart = autostart)

  nothing
end

function set_files_mod() :: Nothing
  for f in ["Manifest.toml", "Project.toml", "routes.jl"]
    try
      chmod(f, 0o700)
    catch _
      @warn "Can't change mod for $f. File might be read-only."
    end
  end

  nothing
end

"""
    newapp_webservice(path::String = "."; autostart::Bool = true, dbsupport::Bool = false) :: Nothing

Template for scaffolding a new Genie app suitable for nimble web services.

# Arguments
- `path::String`: the name of the app and the path where to bootstrap it
- `autostart::Bool`: automatically start the app once the file structure is created
- `dbsupport::Bool`: bootstrap the files needed for DB connection setup via the SearchLight ORM
"""
function newapp_webservice(path::String = "."; autostart::Bool = true, dbsupport::Bool = false) :: Nothing
  newapp(path, autostart = autostart, fullstack = false, dbsupport = dbsupport, mvcsupport = false)
end

"""
    newapp_mvc(path::String = "."; autostart::Bool = true) :: Nothing

Template for scaffolding a new Genie app suitable for MVC web applications (includes MVC structure and DB support).

# Arguments
- `path::String`: the name of the app and the path where to bootstrap it
- `autostart::Bool`: automatically start the app once the file structure is created
"""
function newapp_mvc(path::String = "."; autostart::Bool = true) :: Nothing
  newapp(path, autostart = autostart, fullstack = false, dbsupport = true, mvcsupport = true)
end

"""
    newapp_fullstack(path::String = "."; autostart::Bool = true) :: Nothing

Template for scaffolding a new Genie app suitable for full stack web applications (includes MVC structure, DB support, and frontend asset pipeline).

# Arguments
- `path::String`: the name of the app and the path where to bootstrap it
- `autostart::Bool`: automatically start the app once the file structure is created
"""
function newapp_fullstack(path::String = "."; autostart::Bool = true) :: Nothing
  newapp(path, autostart = autostart, fullstack = true, dbsupport = true, mvcsupport = true)
end

end


text:
import HoTRG_lemon.SpinModule:
	initialize!

type Classical3dSquareLattice{T} <: Classical3dLattice{T}
	spinmodel::ClassicalSpinModel{T}

	tensorT::Array{T,6}
	tensorTtilde::Array{T,6}
 	function Classical3dSquareLattice{T}(spinmodel::ClassicalSpinModel{T})
		this = new{T}()
 		this.spinmodel = spinmodel
 		initialize!(this)
 		return this
 	end
end

Classical3dSquareLattice{T}(spinmodel::ClassicalSpinModel{T}) = Classical3dSquareLattice{T}(spinmodel)

#---
# initialize functions for constructor
function initialize!{T}(lattice::Classical3dSquareLattice{T})
	lattice.tensorT = initializeTensorT3D(lattice.spinmodel)
	lattice.tensorTtilde = initializeTensorT3D(lattice.spinmodel; tilde = true)
end

function initializeTensorT3D{T}(spinmodel::ClassicalSpinModel{T}; tilde::Bool = false)
	if tilde
		measureOperator = getMeasureOperator(spinmodel)
	else
		measureOperator = ones(T,getStates(spinmodel))
	end
	tensorW = getFactorW(spinmodel)
	leng = size(tensorW,1)==size(tensorW,2) ? size(tensorW,1) : throw(DomainError())
	tensorT = zeros(T,leng,leng,leng,leng,leng,leng)
	for a = 1:leng, l = 1:leng, r = 1:leng, f=1:leng, b=1:leng, u=1:leng, v=1:leng
		tensorT[l,r,f,b,u,v] += measureOperator[a] *
		tensorW[a,l] * tensorW[a,r] * tensorW[a,f] * tensorW[a,b] *
		tensorW[a,u] * tensorW[a,v]
	end
	return tensorT
end

#---
# functions for Classical3dSquareLattice
function getCoarserate(lattice::Classical3dSquareLattice)
	return 2
end

#---
# functions for test
function testClassical3dSquareLattice()
	spinmodel = testIsingModel()
	return buildLattice("classical_3d_square", spinmodel)
end


text:
using Arpack

function circle_line_intersections(a1::T, a2::T, b::T, r::T) where T
    """
    Calculates the points intersecting points between the 2d plane defined
    as a1*x + a2*y = b
    and the circle x^2 + y^2 = r^2
    Returns x1, y1, x2, y2 such that [x1; y1] and [x2; y2] are the two intersecting points.
    If no intersection(s) exist then NaN values are returned.
    Reference: http://mathworld.wolfram.com/Circle-LineIntersection.html
    """
    # Normalize, i.e. make [a1; a2] unit norm
    norm_a = sqrt(a1^2 + a2^2)
    if norm_a < 1e-11  # Very badly scaled problem
        @warn "Badly scaled circe-line"
        return T(NaN), T(NaN), T(NaN), T(NaN)
    end
    a1 /= norm_a; a2 /= norm_a; b /= norm_a

    # First calculate a point (x, y) on the line
    # and dx, dy such that (x + dx, y + dy) lies also on the line.
    if abs(a1) > abs(a2) # Avoid dividing with something small
        x = (b - a2)/a1
        y = one(T)

        dx = -a2/a1
        dy = one(T)
    else
        x = one(T)
        y = (b - a1)/a2

        dx = one(T)
        dy = -a1/a2
    end

    # Make (dx, dy) unit norm
    norm_d = sqrt(dx^2 + dy^2)
    dx /= norm_d; dy /= norm_d

    D = x*(y + dy) - (x + dx)*y
    #= Allow for slightly infeasible solutions
    if abs(r)/abs(D) <= one(T)
        @show abs(r)/abs(D) - 1
        if abs(r)/abs(D) >= one(T) - 1e-11
            D = r
        end
    end
    =#
    if r >= abs(D)
        c = sqrt(r^2 - D^2)

        x1 = D*dy; y1 = -D*dx
        x2 = x1;  y2 = y1;

        δ1 = dx*c
        if dy < 0; δ1 = -δ1 end
        δ2 = abs(dy)*c

        x1 += δ1; y1 += δ2
        x2 -= δ1; y2 -= δ2

        return x1, y1, x2, y2
    else
        return T(NaN), T(NaN), T(NaN), T(NaN)
    end
end

function trs_robust(P::AbstractArray{T}, q::AbstractVector{T}, r::T; tol=0.0, v0=zeros((0,)), kwargs...) where T
    n = length(q)
    if n < 15
        X, info =  trs_boundary_small(Symmetric(P), q, r; kwargs...)
    else
        try
            X, info = trs_boundary(Symmetric(P), q, r; tol=tol, v0=v0, kwargs...)
        catch e
            @warn "Error:", e, " - switching to a direct solver"
            X, info = trs_boundary_small(Symmetric(P), q, r; kwargs...)
        end
    end
    # @show info, norm(P*X[:,1] + info.λ[1]*C*X[:, 1] + q)
    return X, info
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:2] <= 1)

                     @variable(m, -6.103907844236108 <= q <= 8.378675476959124)

                     add_NL_constraint(m, :(softplus(-0.8297607126794757 + -0.5064858067335005*softplus(-0.6188655598844357 + -0.38627203797486676*$(x[1]) + 0.4878933200630442*$(x[2])) + -0.8619658571743893*softplus(-0.7914291135758744 + 0.8682104207735972*$(x[1]) + 0.3684602934879275*$(x[2])) + 0.8415580445654567*softplus(-0.26405078449033637 + -0.6452849815457378*$(x[1]) + 0.9586530591287103*$(x[2]))) + softplus(-0.0809505831099564 + 0.6185835236357913*softplus(-0.6188655598844357 + -0.38627203797486676*$(x[1]) + 0.4878933200630442*$(x[2])) + -0.5288114683363174*softplus(-0.7914291135758744 + 0.8682104207735972*$(x[1]) + 0.3684602934879275*$(x[2])) + -0.42070009030061817*softplus(-0.26405078449033637 + -0.6452849815457378*$(x[1]) + 0.9586530591287103*$(x[2]))) + softplus(0.16011867292674786 + -0.14898607305663658*softplus(-0.6188655598844357 + -0.38627203797486676*$(x[1]) + 0.4878933200630442*$(x[2])) + -0.8220978033005912*softplus(-0.7914291135758744 + 0.8682104207735972*$(x[1]) + 0.3684602934879275*$(x[2])) + -0.8518875556726564*softplus(-0.26405078449033637 + -0.6452849815457378*$(x[1]) + 0.9586530591287103*$(x[2]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using ReservoirComputing

train_data = ones(Int, 2, 10)
rule = 30
generations = 3
expansion_size = 10
permutations = 2

reca = RECA_discrete(train_data, rule, generations, expansion_size, permutations)

W_out = ESNtrain(reca, 0.001, train_data = ones(Float64, 2, 10))
output = RECAdirect_predict_discrete(reca, W_out, ones(Int, 2, 10))

@test isequal(size(output), size(ones(Int, 2, 10)))


text:
push!(LOAD_PATH, joinpath(@__DIR__, "..", ".."))

import ClimaCore.Geometry, LinearAlgebra, UnPack
import ClimaCore:
    Fields,
    Domains,
    Topologies,
    Meshes,
    DataLayouts,
    Operators,
    Geometry,
    Spaces

using OrdinaryDiffEq: ODEProblem, solve, SSPRK33

using Logging: global_logger
using TerminalLoggers: TerminalLogger
global_logger(TerminalLogger())

const FT = Float64

a = FT(0.0)
b = FT(4pi)
n = 128
α = FT(0.1)

domain = Domains.IntervalDomain(
    Geometry.ZPoint{FT}(a),
    Geometry.ZPoint{FT}(b),
    boundary_tags = (:left, :right),
)
mesh = Meshes.IntervalMesh(domain, nelems = n)

cs = Spaces.CenterFiniteDifferenceSpace(mesh)
fs = Spaces.FaceFiniteDifferenceSpace(cs)

V = Geometry.Cartesian3Vector.(ones(FT, fs))
θ = sin.(Fields.coordinate_field(cs).z)

# Solve advection Equation: ∂θ/dt = -∂(vθ)

# upwinding
function tendency1!(dθ, θ, _, t)
    UB = Operators.UpwindBiasedProductC2F(
        left = Operators.SetValue(sin(a - t)),
        right = Operators.SetValue(sin(b - t)),
    )
    ∂ = Operators.DivergenceF2C()

    return @. dθ = -∂(UB(V, θ))
end

# use the advection operator
function tendency2!(dθ, θ, _, t)
    A = Operators.AdvectionC2C(
        left = Operators.SetValue(sin(-t)),
        right = Operators.Extrapolate(),
    )
    return @. dθ = -A(V, θ)
end

@show tendency2!(similar(θ), θ, nothing, 0.0)

# Solve the ODE operator
Δt = 0.01
prob = ODEProblem(tendency2!, θ, (0.0, 10.0))
sol = solve(
    prob,
    SSPRK33(),
    dt = Δt,
    saveat = 10 * Δt,
    progress = true,
    progress_message = (dt, u, p, t) -> t,
);

ENV["GKSwstype"] = "nul"
import Plots
Plots.GRBackend()

dirname = "advect"
path = joinpath(@__DIR__, "output", dirname)
mkpath(path)

anim = Plots.@animate for u in sol.u
    Plots.plot(u, xlim = (-1, 1))
end
Plots.mp4(anim, joinpath(path, "advect.mp4"), fps = 10)
Plots.png(
    Plots.plot(sol.u[end], xlim = (-1, 1)),
    joinpath(path, "advect_end.png"),
)

function linkfig(figpath, alt = "")
    # buildkite-agent upload figpath
    # link figure in logs if we are running on CI
    if get(ENV, "BUILDKITE", "") == "true"
        artifact_url = "artifact://$figpath"
        print("\033]1338;url='$(artifact_url)';alt='$(alt)'\a\n")
    end
end

linkfig("output/$(dirname)/advect_end.png", "Advect End Simulation")