text:
isgpu(x) = false
function __init__()
    @require CuArrays="3a865a2d-5b23-5a0f-bc46-62713ec82fae" begin
        gpu_or_cpu(x::CuArrays.CuArray) = CuArrays.CuArray
        gpu_or_cpu(x::TrackedArray{<:Any,<:Any,<:CuArrays.CuArray}) = CuArrays.CuArray
        gpu_or_cpu(x::Transpose{<:Any,<:CuArrays.CuArray}) = CuArrays.CuArray
        gpu_or_cpu(x::Adjoint{<:Any,<:CuArrays.CuArray}) = CuArrays.CuArray
        gpu_or_cpu(x::Adjoint{<:Any,TrackedArray{<:Any,<:Any,<:CuArrays.CuArray}}) = CuArrays.CuArray
        gpu_or_cpu(x::Transpose{<:Any,TrackedArray{<:Any,<:Any,<:CuArrays.CuArray}}) = CuArrays.CuArray
        isgpu(::CuArrays.CuArray) = true
        isgpu(::TrackedArray{<:Any,<:Any,<:CuArrays.CuArray}) = true
        isgpu(::Transpose{<:Any,<:CuArrays.CuArray}) = true
        isgpu(::Adjoint{<:Any,<:CuArrays.CuArray}) = true
        isgpu(::Adjoint{<:Any,TrackedArray{<:Any,<:Any,<:CuArrays.CuArray}}) = true
        isgpu(::Transpose{<:Any,TrackedArray{<:Any,<:Any,<:CuArrays.CuArray}}) = true
    end

    @require NLopt="76087f3c-5699-56af-9a33-bf431cd00edd" begin
        function sciml_train(loss, θ, opt::NLopt.Opt, data = DEFAULT_DATA; maxeval=100)
            local x, cur, state
            cur,state = iterate(data)

            function nlopt_grad!(θ,grad)
              _x,lambda = Flux.Zygote.pullback(θ) do θ
                x = loss(θ)
                first(x)
              end

              if length(grad) > 0
                grad .= first(lambda(1))
                @printf("fval:%2.2e  norm:%2.2e\n", _x[1], norm(θ))
              end

              return _x
            end

            NLopt.min_objective!(opt, nlopt_grad!)
            NLopt.maxeval!(opt, maxeval)

            t0= time()
            (minf,minx,ret) = NLopt.optimize(opt, θ)
            _time = time()

            Optim.MultivariateOptimizationResults(opt,
                                                    θ,# initial_x,
                                                    minx, #pick_best_x(f_incr_pick, state),
                                                    minf, # pick_best_f(f_incr_pick, state, d),
                                                    maxeval, #iteration,
                                                    maxeval >= maxeval, #iteration == options.iterations,
                                                    false, # x_converged,
                                                    0.0,#T(options.x_tol),
                                                    0.0,#T(options.x_tol),
                                                    NaN,# x_abschange(state),
                                                    NaN,# x_abschange(state),
                                                    false,# f_converged,
                                                    0.0,#T(options.f_tol),
                                                    0.0,#T(options.f_tol),
                                                    NaN,#f_abschange(d, state),
                                                    NaN,#f_abschange(d, state),
                                                    false,#g_converged,
                                                    0.0,#T(options.g_tol),
                                                    NaN,#g_residual(d),
                                                    false, #f_increased,
                                                    nothing,
                                                    maxeval,
                                                    maxeval,
                                                    0,
                                                    ret,
                                                    NaN,
                                                    _time-t0,)
        end

        function sciml_train(loss, θ, opt::NLopt.Opt, lower_bounds, upper_bounds, data = DEFAULT_DATA; maxeval=100, nstart=1)
          local x, cur, state
          cur,state = iterate(data)

          function nlopt_grad!(θ,grad)
            _x,lambda = Flux.Zygote.pullback(θ) do θ
              x = loss(θ)
              first(x)
            end

            if length(grad) > 0
              grad .= first(lambda(1))
            end

            return _x
          end

          NLopt.min_objective!(opt, nlopt_grad!)
          NLopt.lower_bounds!(opt, lower_bounds)
          NLopt.upper_bounds!(opt, upper_bounds)
          NLopt.maxeval!(opt, maxeval)

          if nstart > 1
            localopt = opt
            opt = NLopt.Opt(:G_MLSL_LDS, length(lower_bounds))
            NLopt.lower_bounds!(opt, lower_bounds)
            NLopt.upper_bounds!(opt, upper_bounds)
            NLopt.local_optimizer!(opt, localopt)
            NLopt.maxeval!(opt, nstart * maxeval)
          end

          t0 = time()
          NLopt.optimize(opt, θ)
          _time = time()

          Optim.MultivariateOptimizationResults(opt,
                                                    θ,# initial_x,
                                                    minx, #pick_best_x(f_incr_pick, state),
                                                    minf, # pick_best_f(f_incr_pick, state, d),
                                                    maxeval, #iteration,
                                                    maxeval >= maxeval, #iteration == options.iterations,
                                                    false, # x_converged,
                                                    0.0,#T(options.x_tol),
                                                    0.0,#T(options.x_tol),
                                                    NaN,# x_abschange(state),
                                                    NaN,# x_abschange(state),
                                                    false,# f_converged,
                                                    0.0,#T(options.f_tol),
                                                    0.0,#T(options.f_tol),
                                                    NaN,#f_abschange(d, state),
                                                    NaN,#f_abschange(d, state),
                                                    false,#g_converged,
                                                    0.0,#T(options.g_tol),
                                                    NaN,#g_residual(d),
                                                    false, #f_increased,
                                                    nothing,
                                                    maxeval,
                                                    maxeval,
                                                    0,
                                                    ret,
                                                    NaN,
                                                    _time-t0,)
        end
    end

    @require MultistartOptimization = "3933049c-43be-478e-a8bb-6e0f7fd53575" begin
        function sciml_train(loss, _θ, opt::MultistartOptimization.TikTak, data = DEFAULT_DATA;lower_bounds, upper_bounds, local_method,
                              maxiters = get_maxiters(data), kwargs...)
          local x, cur, state
          cur,state = iterate(data)

          _loss = function (θ)
            x = loss(θ,cur...)
            first(x)
          end

          t0 = time()

          P = MultistartOptimization.MinimizationProblem(_loss, lower_bounds, upper_bounds)
          multistart_method = opt
          local_method = MultistartOptimization.NLoptLocalMethod(local_method)
          p = MultistartOptimization.multistart_minimization(multistart_method, local_method, P)

          t1 = time()

          Optim.MultivariateOptimizationResults(opt,
                                                [NaN],# initial_x,
                                                p.location, #pick_best_x(f_incr_pick, state),
                                                p.value, # pick_best_f(f_incr_pick, state, d),
                                                0, #iteration,
                                                false, #iteration == options.iterations,
                                                false, # x_converged,
                                                0.0,#T(options.x_tol),
                                                0.0,#T(options.x_tol),
                                                NaN,# x_abschange(state),
                                                NaN,# x_abschange(state),
                                                false,# f_converged,
                                                0.0,#T(options.f_tol),
                                                0.0,#T(options.f_tol),
                                                NaN,#f_abschange(d, state),
                                                NaN,#f_abschange(d, state),
                                                false,#g_converged,
                                                0.0,#T(options.g_tol),
                                                NaN,#g_residual(d),
                                                false, #f_increased,
                                                nothing,
                                                maxiters,
                                                maxiters,
                                                0,
                                                true,
                                                NaN,
                                                t1 - t0)
        end
    end
    @require QuadDIRECT = "dae52e8d-d666-5120-a592-9e15c33b8d7a" begin
        export QuadDirect
        struct QuadDirect
        end
        function sciml_train(loss, _θ, opt::QuadDirect, data = DEFAULT_DATA;lower_bounds, upper_bounds, splits,
                              maxiters = get_maxiters(data), kwargs...)
          local x, cur, state
          cur,state = iterate(data)

          _loss = function (θ)
            x = loss(θ,cur...)
            first(x)
          end

          t0 = time()

          root, x0 = QuadDIRECT.analyze(_loss, splits, lower_bounds, upper_bounds; kwargs...)

          t1 = time()

          Optim.MultivariateOptimizationResults(opt,
                                                [NaN],# initial_x,
                                                QuadDIRECT.position(root), #pick_best_x(f_incr_pick, state),
                                                QuadDIRECT.minimum(root), # pick_best_f(f_incr_pick, state, d),
                                                0, #iteration,
                                                false, #iteration == options.iterations,
                                                false, # x_converged,
                                                0.0,#T(options.x_tol),
                                                0.0,#T(options.x_tol),
                                                NaN,# x_abschange(state),
                                                NaN,# x_abschange(state),
                                                false,# f_converged,
                                                0.0,#T(options.f_tol),
                                                0.0,#T(options.f_tol),
                                                NaN,#f_abschange(d, state),
                                                NaN,#f_abschange(d, state),
                                                false,#g_converged,
                                                0.0,#T(options.g_tol),
                                                NaN,#g_residual(d),
                                                false, #f_increased,
                                                nothing,
                                                maxiters,
                                                maxiters,
                                                0,
                                                true,
                                                NaN,
                                                t1 - t0)
        end
    end
end


text:
using GZip
const file0 = "GFP ChIP-SEQ-antibody1.fastq.gz"
const file1 = "NeuroD2 ChIP-SEQ-antibody1.fastq.gz"

# The files are too big, we'll just write an iterator that reads from
# the file as needed.

import Base: start, done, next

type GulayseData
    nbatch::Int
    flip::Bool
    f0::GZipStream
    f1::GZipStream
    x::Array{Float32,4}
    y::Array{Float32,2}
    GulayseData(n=128,f=true)=new(n,f)
end

function start(d::GulayseData)
    isdefined(d,:f0) && close(f0)
    isdefined(d,:f1) && close(f1)
    d.f0 = GZip.open(file0)
    d.f1 = GZip.open(file1)
    d.x = zeros(Float32, 1, 51, 4, d.nbatch)
    d.y = zeros(Float32, 2, d.nbatch)
    return nothing
end

function done(d::GulayseData, s)
    eof(d.f0) || eof(d.f1)
end

function next(d::GulayseData, s)
    d.x[:] = d.y[:] = 0
    for n=1:d.nbatch      # how do we handle train vs test?
        d.flip ? d.y[1,n]=1 : d.y[2,n]=1 # do not have a loss function for 1-D output
        xread(d.flip ? d.f1 : d.f0, d.x, n) # alternate reading from positive and negative examples
        d.flip = !d.flip
    end
    return ((d.x, d.y), nothing)
end

function xread(fp::GZipStream, x::Array{Float32,4}, n::Int, s::ASCIIString="")
    for i=1:4
        eof(fp) && return
        tmp = readline(fp)
        i==2 && (s = tmp)
    end
    length(s)==52 || error("length(s)==$(length(s))")
    for i=1:51
        s[i] == 'A' ? x[1,i,1,n] = 1 :
        s[i] == 'C' ? x[1,i,2,n] = 1 :
        s[i] == 'G' ? x[1,i,3,n] = 1 :
        s[i] == 'T' ? x[1,i,4,n] = 1 :
        s[i] == 'N' ? nothing :
        error("s[$i]==$(s[i])")
    end
end

ext: jl

max_stars_repo_name: JuliaTagBot/BindingPrediction

max_issues_repo_name: JuliaTagBot/BindingPrediction

max_forks_repo_name: JuliaTagBot/BindingPrediction


text:
# --------------------------------------------------------------------------
# ACE.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# All rights reserved.
# --------------------------------------------------------------------------

using LinearAlgebra: norm, dot

import JuLIP: evaluate!, evaluate_d!, alloc_temp
import JuLIP.MLIPs: alloc_B, alloc_dB, IPBasis

struct OneOrthogonal{T,B} <: IPBasis
   coeffs::Matrix{T}
   nested::B
   # ----------------- used only for construction ...
   #                   but useful to have since it defines the notion or orth.
   tdf::Vector{T}
   ww::Vector{T}
end

Base.length(P::OneOrthogonal) = length(P.nested)

ACE.rand_radial(J::OneOrthogonal) = ACE.rand_radial(J.nested)

"""
    OneOrthogonal(P, tdf = P.tdf, ww = P.ww)

Orthogonal basis `Q` such that `span(Q) == span(P)` and `Q[:,1:end-1]` are
orthogonal to `1`, where `P` is an orthogonal basis with respect to the inner
product defined by `(tdf,ww)`.
"""
function OneOrthogonal(P, tdf=P.tdf, ww = P.ww)
   T = eltype(alloc_B(P))
   dotw = (f1, f2) -> dot(f1,ww.*f2)
   normw = f->sqrt(dotw(f1,f2))

   tmp = alloc_temp(P)
   p = Matrix{T}(undef, length(tdf), length(P))
   for i = 1:length(tdf)
      evaluate!(@view(p[i,:]), tmp, P, tdf[i])
   end

   rotate! = (qe,p1,p2) -> begin
      local c = (dot(p1,ww), dot(p2,ww))
      c = c ./ norm(c)
      qe .= c[1]*p1 .+ c[2]*p2
      return c
   end

   qe = Vector{T}(undef, length(tdf))
   c = Matrix{T}(undef, 2, length(P)-1)
   c[:,1] .= rotate!(qe, p[:,1], p[:,2])
   for j = 2:length(P)-1
      c[:,j] .= rotate!(qe, p[:,j+1],qe)
   end

   return OneOrthogonal(c, P, tdf, ww)
end

alloc_B( J::OneOrthogonal{T}) where {T} = zeros(T, length(J))
alloc_dB(J::OneOrthogonal{T}) where {T} = zeros(T, length(J))

alloc_B( J::OneOrthogonal{T}, x::TX) where {T, TX} = zeros(TX, length(J))
alloc_dB(J::OneOrthogonal{T}, x::TX) where {T, TX} = zeros(TX, length(J))

function evaluate!(q, tmp, Q::OneOrthogonal, t)
   P = Q.nested
   c = Q.coeffs
   evaluate!(q, tmp, P, t)
   q[1],qe = c[2,1]*q[1] - c[1,1]*q[2], c[1,1]*q[1] + c[2,1]*q[2]
   for j = 2:length(Q)-1
       q[j],qe = c[2,j]*q[j+1] - c[1,j]*qe, c[1,j]*q[j+1] + c[2,j]*qe
   end
   q[end] = qe
   return q
end

function evaluate_d!(q, dq, tmp, Q::OneOrthogonal, t)
   P = Q.nested
   c = Q.coeffs
   evaluate_d!(q, dq, tmp, P, t)
    q[1], qe = c[2,1]* q[1] - c[1,1]* q[2], c[1,1]* q[1] + c[2,1]* q[2]
   dq[1],dqe = c[2,1]*dq[1] - c[1,1]*dq[2], c[1,1]*dq[1] + c[2,1]*dq[2]
   for j = 2:length(Q)-1
        q[j], qe = c[2,j]* q[j+1] - c[1,j]* qe, c[1,j]* q[j+1] + c[2,j]* qe
       dq[j],dqe = c[2,j]*dq[j+1] - c[1,j]*dqe, c[1,j]*dq[j+1] + c[2,j]*dqe
   end
    q[end] = qe
   dq[end] = dqe
   return dq
end

# ------------------------------------------------------------------
#   a filtering mechanism to sort out the non-orthogonal functions

function filter_oneorth(ν, P::OneOrthogonal)
   degP = length(P)-1
   # accept all ν with ki < degP, i.e., don't allow the last function
   return maximum(ν.kk) < degP
end


text:
function calculateTransfiniteNodes(c1::Function, c2::Function, c3::Function, c4::Function, P12::Vector{Float64},
    P14::Vector{Float64}, P34::Vector{Float64}, P32::Vector{Float64}, us::Vector{Float64}, vs::Vector{Float64})

    S = Array{Union{Nothing, Tuple}}(nothing, size(us,1), size(vs,1))

    for (i, u) in enumerate(us)
        for (j, v) in enumerate(vs)
            tmp = (1-v)*c1(u) + v*c3(u) + (1-u)*c2(v) + u*c4(v) - ( (1-u)*(1-v)*P12 + u*v*P34 + u*(1-v)*P14 + (1-u)*v*P32 )
            S[i,j] = tuple(tmp...)
        end
    end

    return S

end


text:
using PrettyPlots

# Set The Problem Up
main_dir=pwd()
results_name = "testing_clp76"; # Name Folder to Store Comparison Results in
results_dir = string(main_dir,"/results/",results_name)
warn("you may have the wrong senario!")
senario = "s1.jl"  #TODO save this with the actual test data!

# define vehicle parameters; has constraints on vehicle model for optimization
pa = Vpara();   # initialize parameter set
@unpack_Vpara pa # unpack parameters

include(string(main_dir,"/initialize.jl"))

# Choose Results to Compare
k = 4 # number of results to compare
path = ["","","",""] # should match k
path[1] = "/home/febbo/Documents/workspace/OCP/results/testing_clp76eulerIPOPT";
path[2] = "/home/febbo/Documents/workspace/OCP/results/testing_clp76trap2IPOPT";
path[3] = "/home/febbo/Documents/workspace/OCP/results/testing_clp76eulerKNITRO";
path[4] = "/home/febbo/Documents/workspace/OCP/results/testing_clp76KNITRO";

# Plot Settings
label_string = ["Euler & IPOPT", "Trap. & IPOPT","Euler & KNITRO", "Trap. & KNITRO"];

if isdir(results_dir)
		rm(results_dir; recursive=true)
		print("\n The old results have all been deleted! \n \n")
end

dfs=data_sets(k,"STATES")  # put the data from each test into a DataFrame
dfs_opt=data_sets(k,"OPT_info")
mkdir(results_dir)
cd(results_dir)

sa=saplot(dfs,label_string,false,pa)
sr=srplot(dfs,label_string,false,pa)
yaw=yawplot(dfs,label_string,false,pa)
yawr=yawrplot(dfs,label_string,false,pa)

longv=longvplot(dfs,label_string,false,pa)
ax=axplot(dfs,label_string,false,pa)
jx=jxplot(dfs,label_string,false,pa)

latv=latvplot(dfs,label_string,false,pa)

pp=pplot(dfs,1,label_string,false,false,obs_data,s_data,pa)

lt=ltplot(dfs,label_string,false,pa)
vt=vtplot(dfs,label_string,false,pa)

# make a figure with all states and controls
plot(longv, ax, jx, sa, sr, yaw, yawr, latv, pp, size=(1800,1200))
savefig("main.png")

# make a figure with the tire forces
plot(lt, vt, size=(1800,1200),layout = grid(2,1) )
savefig("tire_f.png")

# make a simulation with all vehicles
#panim_fun(label_string)

optplot(dfs_opt,label_string,false)

tplot(dfs_opt,dfs,label_string,false)

cd(main_dir)


text:
using Test, QuadGK, GrapheneQFT, Cubature

function random_atom()
      u = rand(-20:20)
      v = rand(-20:20)
      return GrapheneQFT.GrapheneCoord(
            u,
            v,
            rand() < 0.5 ? GrapheneQFT.A : GrapheneQFT.B,
      )
end

r1 = rand(-20:20, 10)
r2 = rand(-20:20, 10)
@test map((x, y) -> graphene_A(x, y), r1, r2) ==
      map((x, y) -> GrapheneQFT.GrapheneCoord(x, y, GrapheneQFT.A), r1, r2)
@test map((x, y) -> graphene_B(x, y), r1, r2) ==
      map((x, y) -> GrapheneQFT.GrapheneCoord(x, y, GrapheneQFT.B), r1, r2)

a1 = graphene_A(0, 17)
a2 = graphene_A(1, 1)
a3 = graphene_B(3, -2)
a4 = graphene_A(2, 7)
a5 = graphene_A(0, 1)

ϵ1 = 1.7
ϵ2 = -0.2

V1 = 0.4
V2 = 2.0
V3 = -7.1

imp1 = ImpurityState(ϵ1, [(V1, a4), (V2, a3)])
imp2 = ImpurityState(ϵ2, [(V3, a5)])

c1 = 3.0
c2 = 5.0
c3 = -9.0

p1 = (a1, a1, c1)
p2 = (a2, a1, c2)
p3 = (a3, a1, c3)

my_system = mkGrapheneSystem(0.0, 0.0, [imp1, imp2], [p1, p2, p3])

@test my_system.imps == [ϵ1, ϵ2]
@test my_system.scattering_atoms == [a5, a2, a4, a1, a3]

@test my_system.Δ == [
      0 0 0 0 0
      0 0 0 c2 0
      0 0 0 0 0
      0 c2 0 c1 c3
      0 0 0 c3 0
]
@test my_system.V == [
      0 V3
      0 0
      V1 0
      0 0
      V2 0
]

rand_num = ((rand() - 1 / 2) + 1im * (rand() - 1 / 2)) * 20
scattering_pairs =
      [
            (x, y) for x in my_system.scattering_atoms,
            y in my_system.scattering_atoms
      ] |> vec

@test G_R(rand_num, scattering_pairs, my_system) |> real ≈
      G_R(conj.(rand_num), scattering_pairs, my_system) |> real

@test G_R(rand_num, scattering_pairs, my_system) |> imag ≈
      -G_R(conj.(rand_num), scattering_pairs, my_system) |> imag

@test abs(
      hcubature(
            r -> 2 * π * Ψ_pz(r[1], r[2]) .^ 2 * r[1]^2 * sin(r[2]),
            [0, 0],
            [40, π],
      )[1] - 1,
) < 1e-6

R_rand = 40 * rand()
τ_rand = π * rand()
exact_coulomb = coulomb_potential_pz(R_rand, τ_rand)[1]
@test abs(coulomb_potential_pz_interp(R_rand, τ_rand) - exact_coulomb) /
      exact_coulomb < 1e-2


text:
# Compile-time constants
# logic derived from
# https://github.com/open-mpi/ompi/blob/master/ompi/include/mpi.h.in

# Implementation limits:
const MPI_MAX_DATAREP_STRING = Cint(128)
const MPI_MAX_ERROR_STRING = Cint(256)
const MPI_MAX_INFO_KEY = Cint(36)
const MPI_MAX_INFO_VAL = Cint(256)
const MPI_MAX_LIBRARY_VERSION_STRING = Cint(256)
const MPI_MAX_OBJECT_NAME = Cint(64)
const MPI_MAX_PORT_NAME = Cint(1024)
const MPI_MAX_PROCESSOR_NAME = Cint(256)

# Types

# Various (signed) integer types:
const MPI_Aint = Int
const MPI_Fint = Int32
const MPI_Count = Int
const MPI_Offset = Int

# Status:
struct MPI_Status
    MPI_SOURCE::Cint
    MPI_TAG::Cint
    MPI_ERROR::Cint
    _private0::Cint
    _private1::Csize_t
end

# Opaque handles:
const MPI_Comm = Ptr{Cvoid}
const MPI_Datatype = Ptr{Cvoid}
const MPI_Errhandler = Ptr{Cvoid}
const MPI_File = Ptr{Cvoid}
const MPI_Group = Ptr{Cvoid}
const MPI_Info = Ptr{Cvoid}
const MPI_Message = Ptr{Cvoid}
const MPI_Op = Ptr{Cvoid}
const MPI_Request = Ptr{Cvoid}
const MPI_Win = Ptr{Cvoid}

# Function pointers:
const MPI_Comm_copy_attr_function = Ptr{Cvoid}
const MPI_Comm_delete_attr_function = Ptr{Cvoid}
const MPI_Comm_errhandler_function = Ptr{Cvoid}
const MPI_Comm_errhandler_fn = MPI_Comm_errhandler_function
const MPI_Copy_function = MPI_Comm_copy_attr_function
const MPI_Datarep_conversion_function = Ptr{Cvoid}
const MPI_Datarep_extent_function = Ptr{Cvoid}
const MPI_Delete_function = Ptr{Cvoid}
const MPI_File_errhandler_function = Ptr{Cvoid}
const MPI_File_errhandler_fn = Ptr{Cvoid}
const MPI_Grequest_cancel_function = Ptr{Cvoid}
const MPI_Grequest_free_function = Ptr{Cvoid}
const MPI_Grequest_query_function = Ptr{Cvoid}
const MPI_Type_copy_attr_function = Ptr{Cvoid}
const MPI_Type_delete_attr_function = Ptr{Cvoid}
const MPI_User_function = Ptr{Cvoid}
const MPI_Win_copy_attr_function = Ptr{Cvoid}
const MPI_Win_delete_attr_function = Ptr{Cvoid}
const MPI_Win_errhandler_function = Ptr{Cvoid}
const MPI_Win_errhandler_fn = MPI_Win_errhandler_function

# Run-time constants

@const_ref MPI_ANY_SOURCE Cint -1
@const_ref MPI_ANY_TAG    Cint -1
@const_ref MPI_PROC_NULL  Cint -2
@const_ref MPI_ROOT       Cint -4

@const_ref MPI_CART       Cint 1
@const_ref MPI_GRAPH      Cint 2
@const_ref MPI_DIST_GRAPH Cint 3

@const_ref MPI_IDENT      Cint 0
@const_ref MPI_CONGRUENT  Cint 1
@const_ref MPI_SIMILAR    Cint 2
@const_ref MPI_UNEQUAL    Cint 3

@const_ref MPI_KEYVAL_INVALID     Cint -1
@const_ref MPI_UNDEFINED          Cint -32766

@const_ref MPI_TAG_UB             Cint 0
@const_ref MPI_HOST               Cint 1
@const_ref MPI_IO                 Cint 2
@const_ref MPI_WTIME_IS_GLOBAL    Cint 3
@const_ref MPI_APPNUM             Cint 4
@const_ref MPI_LASTUSEDCODE       Cint 5
@const_ref MPI_UNIVERSE_SIZE      Cint 6

@const_ref MPI_WIN_BASE           Cint 7
@const_ref MPI_WIN_SIZE           Cint 8
@const_ref MPI_WIN_DISP_UNIT      Cint 9
@const_ref MPI_WIN_CREATE_FLAVOR  Cint 10
@const_ref MPI_WIN_MODEL          Cint 11

@const_ref MPI_COMBINER_NAMED            Cint 0
@const_ref MPI_COMBINER_DUP              Cint 1
@const_ref MPI_COMBINER_CONTIGUOUS       Cint 2
@const_ref MPI_COMBINER_VECTOR           Cint 3
@const_ref MPI_COMBINER_HVECTOR_INTEGER  Cint 4
@const_ref MPI_COMBINER_HVECTOR          Cint 5
@const_ref MPI_COMBINER_INDEXED          Cint 6
@const_ref MPI_COMBINER_HINDEXED_INTEGER Cint 7
@const_ref MPI_COMBINER_HINDEXED         Cint 8
@const_ref MPI_COMBINER_INDEXED_BLOCK    Cint 9
@const_ref MPI_COMBINER_STRUCT_INTEGER   Cint 10
@const_ref MPI_COMBINER_STRUCT           Cint 11
@const_ref MPI_COMBINER_SUBARRAY         Cint 12
@const_ref MPI_COMBINER_DARRAY           Cint 13
@const_ref MPI_COMBINER_F90_REAL         Cint 14
@const_ref MPI_COMBINER_F90_COMPLEX      Cint 15
@const_ref MPI_COMBINER_F90_INTEGER      Cint 16
@const_ref MPI_COMBINER_RESIZED          Cint 17
@const_ref MPI_COMBINER_HINDEXED_BLOCK   Cint 18

@const_ref MPI_COMM_TYPE_SHARED Cint 0

@const_ref MPI_DISTRIBUTE_BLOCK  Cint 0
@const_ref MPI_DISTRIBUTE_CYCLIC Cint 1
@const_ref MPI_DISTRIBUTE_NONE   Cint 2

@const_ref MPI_SUCCESS                   Cint 0
@const_ref MPI_ERR_BUFFER                Cint 1
@const_ref MPI_ERR_COUNT                 Cint 2
@const_ref MPI_ERR_TYPE                  Cint 3
@const_ref MPI_ERR_TAG                   Cint 4
@const_ref MPI_ERR_COMM                  Cint 5
@const_ref MPI_ERR_RANK                  Cint 6
@const_ref MPI_ERR_REQUEST               Cint 7
@const_ref MPI_ERR_ROOT                  Cint 8
@const_ref MPI_ERR_GROUP                 Cint 9
@const_ref MPI_ERR_OP                    Cint 10
@const_ref MPI_ERR_TOPOLOGY              Cint 11
@const_ref MPI_ERR_DIMS                  Cint 12
@const_ref MPI_ERR_ARG                   Cint 13
@const_ref MPI_ERR_UNKNOWN               Cint 14
@const_ref MPI_ERR_TRUNCATE              Cint 15
@const_ref MPI_ERR_OTHER                 Cint 16
@const_ref MPI_ERR_INTERN                Cint 17
@const_ref MPI_ERR_IN_STATUS             Cint 18
@const_ref MPI_ERR_PENDING               Cint 19
@const_ref MPI_ERR_ACCESS                Cint 20
@const_ref MPI_ERR_AMODE                 Cint 21
@const_ref MPI_ERR_ASSERT                Cint 22
@const_ref MPI_ERR_BAD_FILE              Cint 23
@const_ref MPI_ERR_BASE                  Cint 24
@const_ref MPI_ERR_CONVERSION            Cint 25
@const_ref MPI_ERR_DISP                  Cint 26
@const_ref MPI_ERR_DUP_DATAREP           Cint 27
@const_ref MPI_ERR_FILE_EXISTS           Cint 28
@const_ref MPI_ERR_FILE_IN_USE           Cint 29
@const_ref MPI_ERR_FILE                  Cint 30
@const_ref MPI_ERR_INFO_KEY              Cint 31
@const_ref MPI_ERR_INFO_NOKEY            Cint 32
@const_ref MPI_ERR_INFO_VALUE            Cint 33
@const_ref MPI_ERR_INFO                  Cint 34
@const_ref MPI_ERR_IO                    Cint 35
@const_ref MPI_ERR_KEYVAL                Cint 36
@const_ref MPI_ERR_LOCKTYPE              Cint 37
@const_ref MPI_ERR_NAME                  Cint 38
@const_ref MPI_ERR_NO_MEM                Cint 39
@const_ref MPI_ERR_NOT_SAME              Cint 40
@const_ref MPI_ERR_NO_SPACE              Cint 41
@const_ref MPI_ERR_NO_SUCH_FILE          Cint 42
@const_ref MPI_ERR_PORT                  Cint 43
@const_ref MPI_ERR_QUOTA                 Cint 44
@const_ref MPI_ERR_READ_ONLY             Cint 45
@const_ref MPI_ERR_RMA_CONFLICT          Cint 46
@const_ref MPI_ERR_RMA_SYNC              Cint 47
@const_ref MPI_ERR_SERVICE               Cint 48
@const_ref MPI_ERR_SIZE                  Cint 49
@const_ref MPI_ERR_SPAWN                 Cint 50
@const_ref MPI_ERR_UNSUPPORTED_DATAREP   Cint 51
@const_ref MPI_ERR_UNSUPPORTED_OPERATION Cint 52
@const_ref MPI_ERR_WIN                   Cint 53
@const_ref MPI_T_ERR_MEMORY              Cint 54
@const_ref MPI_T_ERR_NOT_INITIALIZED     Cint 55
@const_ref MPI_T_ERR_CANNOT_INIT         Cint 56
@const_ref MPI_T_ERR_INVALID_INDEX       Cint 57
@const_ref MPI_T_ERR_INVALID_ITEM        Cint 58
@const_ref MPI_T_ERR_INVALID_HANDLE      Cint 59
@const_ref MPI_T_ERR_OUT_OF_HANDLES      Cint 60
@const_ref MPI_T_ERR_OUT_OF_SESSIONS     Cint 61
@const_ref MPI_T_ERR_INVALID_SESSION     Cint 62
@const_ref MPI_T_ERR_CVAR_SET_NOT_NOW    Cint 63
@const_ref MPI_T_ERR_CVAR_SET_NEVER      Cint 64
@const_ref MPI_T_ERR_PVAR_NO_STARTSTOP   Cint 65
@const_ref MPI_T_ERR_PVAR_NO_WRITE       Cint 66
@const_ref MPI_T_ERR_PVAR_NO_ATOMIC      Cint 67
@const_ref MPI_ERR_RMA_RANGE             Cint 68
@const_ref MPI_ERR_RMA_ATTACH            Cint 69
@const_ref MPI_ERR_RMA_FLAVOR            Cint 70
@const_ref MPI_ERR_RMA_SHARED            Cint 71
@const_ref MPI_T_ERR_INVALID             Cint 72
@const_ref MPI_T_ERR_INVALID_NAME        Cint 73
@const_ref MPI_ERR_PROC_ABORTED          Cint 74
@const_ref MPI_ERR_PROC_FAILED           Cint 75
@const_ref MPI_ERR_PROC_FAILED_PENDING   Cint 76
@const_ref MPI_ERR_REVOKED               Cint 77

@const_ref MPI_LOCK_EXCLUSIVE  Cint 1
@const_ref MPI_LOCK_SHARED     Cint 2

@const_ref MPI_MODE_CREATE           Cint   1
@const_ref MPI_MODE_RDONLY           Cint   2
@const_ref MPI_MODE_WRONLY           Cint   4
@const_ref MPI_MODE_RDWR             Cint   8
@const_ref MPI_MODE_DELETE_ON_CLOSE  Cint  16
@const_ref MPI_MODE_UNIQUE_OPEN      Cint  32
@const_ref MPI_MODE_EXCL             Cint  64
@const_ref MPI_MODE_APPEND           Cint 128
@const_ref MPI_MODE_SEQUENTIAL       Cint 256

@const_ref MPI_MODE_NOCHECK          Cint  1
@const_ref MPI_MODE_NOPRECEDE        Cint  2
@const_ref MPI_MODE_NOPUT            Cint  4
@const_ref MPI_MODE_NOSTORE          Cint  8
@const_ref MPI_MODE_NOSUCCEED        Cint 16

@const_ref MPI_ORDER_C       Cint 0
@const_ref MPI_ORDER_FORTRAN Cint 1

@const_ref MPI_SEEK_SET Cint 600
@const_ref MPI_SEEK_CUR Cint 602
@const_ref MPI_SEEK_END Cint 604

@const_ref MPI_THREAD_SINGLE Cint 0
@const_ref MPI_THREAD_FUNNELED Cint 1
@const_ref MPI_THREAD_SERIALIZED Cint 2
@const_ref MPI_THREAD_MULTIPLE Cint 3

@const_ref MPI_TYPECLASS_INTEGER Cint 1
@const_ref MPI_TYPECLASS_REAL    Cint 2
@const_ref MPI_TYPECLASS_COMPLEX Cint 3

@const_ref MPI_ARGV_NULL      Ptr{Cvoid} C_NULL
@const_ref MPI_ARGVS_NULL     Ptr{Cvoid} C_NULL

@const_ref MPI_UNWEIGHTED     Ptr{Cvoid} 2
@const_ref MPI_WEIGHTS_EMPTY  Ptr{Cvoid} 3
@const_ref MPI_BOTTOM         Ptr{Cvoid} 0
@const_ref MPI_IN_PLACE       Ptr{Cvoid} 1

@const_ref MPI_COMM_NULL  MPI_Comm cglobal((:ompi_mpi_comm_null, libmpi))
@const_ref MPI_COMM_SELF  MPI_Comm cglobal((:ompi_mpi_comm_self, libmpi))
@const_ref MPI_COMM_WORLD MPI_Comm cglobal((:ompi_mpi_comm_world, libmpi))

@const_ref MPI_COMM_DUP_FN           MPI_Comm_copy_attr_function   cglobal((:OMPI_C_MPI_COMM_DUP_FN, libmpi))
@const_ref MPI_COMM_NULL_COPY_FN     MPI_Comm_copy_attr_function   cglobal((:OMPI_C_MPI_COMM_NULL_COPY_FN, libmpi))
@const_ref MPI_COMM_NULL_DELETE_FN   MPI_Comm_delete_attr_function cglobal((:OMPI_C_MPI_COMM_NULL_DELETE_FN, libmpi))

@const_ref MPI_DATATYPE_NULL MPI_Datatype cglobal((:ompi_mpi_datatype_null, libmpi))

@const_ref MPI_BYTE MPI_Datatype cglobal((:ompi_mpi_byte, libmpi))
@const_ref MPI_PACKED MPI_Datatype cglobal((:ompi_mpi_packed, libmpi))
@const_ref MPI_CHAR MPI_Datatype cglobal((:ompi_mpi_char, libmpi))
@const_ref MPI_SHORT MPI_Datatype cglobal((:ompi_mpi_short, libmpi))
@const_ref MPI_INT MPI_Datatype cglobal((:ompi_mpi_int, libmpi))
@const_ref MPI_LONG MPI_Datatype cglobal((:ompi_mpi_long, libmpi))
@const_ref MPI_FLOAT MPI_Datatype cglobal((:ompi_mpi_float, libmpi))
@const_ref MPI_DOUBLE MPI_Datatype cglobal((:ompi_mpi_double, libmpi))
# @const_ref MPI_LONG_DOUBLE MPI_Datatype cglobal((:ompi_mpi_long_double, libmpi))
@const_ref MPI_UNSIGNED_CHAR MPI_Datatype cglobal((:ompi_mpi_unsigned_char, libmpi))
@const_ref MPI_SIGNED_CHAR MPI_Datatype cglobal((:ompi_mpi_signed_char, libmpi))
@const_ref MPI_UNSIGNED_SHORT MPI_Datatype cglobal((:ompi_mpi_unsigned_short, libmpi))
@const_ref MPI_UNSIGNED_LONG MPI_Datatype cglobal((:ompi_mpi_unsigned_long, libmpi))
@const_ref MPI_UNSIGNED MPI_Datatype cglobal((:ompi_mpi_unsigned, libmpi))
@const_ref MPI_FLOAT_INT MPI_Datatype cglobal((:ompi_mpi_float_int, libmpi))
@const_ref MPI_DOUBLE_INT MPI_Datatype cglobal((:ompi_mpi_double_int, libmpi))
@const_ref MPI_LONG_DOUBLE_INT MPI_Datatype cglobal((:ompi_mpi_longdbl_int, libmpi))
@const_ref MPI_LONG_INT MPI_Datatype cglobal((:ompi_mpi_long_int, libmpi))
@const_ref MPI_SHORT_INT MPI_Datatype cglobal((:ompi_mpi_short_int, libmpi))
@const_ref MPI_2INT MPI_Datatype cglobal((:ompi_mpi_2int, libmpi))
@const_ref MPI_UB MPI_Datatype cglobal((:ompi_mpi_ub, libmpi))
@const_ref MPI_LB MPI_Datatype cglobal((:ompi_mpi_lb, libmpi))
@const_ref MPI_WCHAR MPI_Datatype cglobal((:ompi_mpi_wchar, libmpi))
@const_ref MPI_LONG_LONG_INT MPI_Datatype cglobal((:ompi_mpi_long_long_int, libmpi))
const MPI_LONG_LONG = MPI_LONG_LONG_INT
@const_ref MPI_UNSIGNED_LONG_LONG MPI_Datatype cglobal((:ompi_mpi_unsigned_long_long, libmpi))
@const_ref MPI_2COMPLEX MPI_Datatype cglobal((:ompi_mpi_2cplex, libmpi))
@const_ref MPI_2DOUBLE_COMPLEX MPI_Datatype cglobal((:ompi_mpi_2dblcplex, libmpi))

@const_ref MPI_INT8_T MPI_Datatype cglobal((:ompi_mpi_int8_t, libmpi))
@const_ref MPI_UINT8_T MPI_Datatype cglobal((:ompi_mpi_uint8_t, libmpi))
@const_ref MPI_INT16_T MPI_Datatype cglobal((:ompi_mpi_int16_t, libmpi))
@const_ref MPI_UINT16_T MPI_Datatype cglobal((:ompi_mpi_uint16_t, libmpi))
@const_ref MPI_INT32_T MPI_Datatype cglobal((:ompi_mpi_int32_t, libmpi))
@const_ref MPI_UINT32_T MPI_Datatype cglobal((:ompi_mpi_uint32_t, libmpi))
@const_ref MPI_INT64_T MPI_Datatype cglobal((:ompi_mpi_int64_t, libmpi))
@const_ref MPI_UINT64_T MPI_Datatype cglobal((:ompi_mpi_uint64_t, libmpi))
@const_ref MPI_AINT MPI_Datatype cglobal((:ompi_mpi_aint, libmpi))
@const_ref MPI_OFFSET MPI_Datatype cglobal((:ompi_mpi_offset, libmpi))
@const_ref MPI_C_BOOL MPI_Datatype cglobal((:ompi_mpi_c_bool, libmpi))
@const_ref MPI_C_FLOAT_COMPLEX MPI_Datatype cglobal((:ompi_mpi_c_float_complex, libmpi))
const MPI_C_COMPLEX = MPI_C_FLOAT_COMPLEX
@const_ref MPI_C_DOUBLE_COMPLEX MPI_Datatype cglobal((:ompi_mpi_c_double_complex, libmpi))
# @const_ref MPI_C_LONG_DOUBLE_COMPLEX MPI_Datatype cglobal((:ompi_mpi_c_long_double_complex, libmpi))
@const_ref MPI_COUNT MPI_Datatype cglobal((:ompi_mpi_count, libmpi))

@const_ref MPI_ERRHANDLER_NULL MPI_Errhandler  cglobal((:ompi_mpi_errhandler_null, libmpi))
@const_ref MPI_ERRORS_ARE_FATAL MPI_Errhandler cglobal((:ompi_mpi_errors_are_fatal, libmpi))
@const_ref MPI_ERRORS_RETURN MPI_Errhandler    cglobal((:ompi_mpi_errors_return, libmpi))
# @const_ref MPI_ERRORS_ABORT MPI_Errhandler     cglobal((:ompi_mpi_errors_abort, libmpi))

@const_ref MPI_FILE_NULL MPI_File cglobal((:ompi_mpi_file_null, libmpi))

@const_ref MPI_GROUP_EMPTY MPI_Group cglobal((:ompi_mpi_group_empty, libmpi))
@const_ref MPI_GROUP_NULL MPI_Group cglobal((:ompi_mpi_group_null, libmpi))

@const_ref MPI_INFO_ENV MPI_Info cglobal((:ompi_mpi_info_env, libmpi))
@const_ref MPI_INFO_NULL MPI_Info cglobal((:ompi_mpi_info_null, libmpi))

@const_ref MPI_MESSAGE_NO_PROC MPI_Message cglobal((:ompi_message_no_proc, libmpi))
@const_ref MPI_MESSAGE_NULL MPI_Message cglobal((:ompi_message_null, libmpi))

@const_ref MPI_DISPLACEMENT_CURRENT MPI_Offset -54278278

@const_ref MPI_OP_NULL MPI_Op cglobal((:ompi_mpi_op_null, libmpi))
@const_ref MPI_MAX     MPI_Op cglobal((:ompi_mpi_op_max, libmpi))
@const_ref MPI_MIN     MPI_Op cglobal((:ompi_mpi_op_min, libmpi))
@const_ref MPI_SUM     MPI_Op cglobal((:ompi_mpi_op_sum, libmpi))
@const_ref MPI_PROD    MPI_Op cglobal((:ompi_mpi_op_prod, libmpi))
@const_ref MPI_LAND    MPI_Op cglobal((:ompi_mpi_op_land, libmpi))
@const_ref MPI_BAND    MPI_Op cglobal((:ompi_mpi_op_band, libmpi))
@const_ref MPI_LOR     MPI_Op cglobal((:ompi_mpi_op_lor, libmpi))
@const_ref MPI_BOR     MPI_Op cglobal((:ompi_mpi_op_bor, libmpi))
@const_ref MPI_LXOR    MPI_Op cglobal((:ompi_mpi_op_lxor, libmpi))
@const_ref MPI_BXOR    MPI_Op cglobal((:ompi_mpi_op_bxor, libmpi))
@const_ref MPI_MINLOC  MPI_Op cglobal((:ompi_mpi_op_minloc, libmpi))
@const_ref MPI_MAXLOC  MPI_Op cglobal((:ompi_mpi_op_maxloc, libmpi))
@const_ref MPI_REPLACE MPI_Op cglobal((:ompi_mpi_op_replace, libmpi))
@const_ref MPI_NO_OP   MPI_Op cglobal((:ompi_mpi_op_no_op, libmpi))

@const_ref MPI_REQUEST_NULL MPI_Request cglobal((:ompi_request_null, libmpi))

@const_ref MPI_STATUS_IGNORE    Ptr{MPI_Status} C_NULL
@const_ref MPI_STATUSES_IGNORE  Ptr{MPI_Status} C_NULL

@const_ref MPI_TYPE_DUP_FN          MPI_Comm_copy_attr_function   cglobal((:OMPI_C_MPI_TYPE_DUP_FN, libmpi), MPI_Comm_copy_attr_function)
@const_ref MPI_TYPE_NULL_COPY_FN    MPI_Type_copy_attr_function   cglobal((:OMPI_C_MPI_TYPE_NULL_COPY_FN, libmpi), MPI_Comm_copy_attr_function)
@const_ref MPI_TYPE_NULL_DELETE_FN  MPI_Type_delete_attr_function cglobal((:OMPI_C_MPI_TYPE_NULL_DELETE_FN, libmpi), MPI_Type_delete_attr_function)

@const_ref MPI_WIN_NULL MPI_Win cglobal((:ompi_mpi_win_null, libmpi))

@const_ref MPI_WIN_DUP_FN           MPI_Win_copy_attr_function   cglobal((:OMPI_C_MPI_WIN_DUP_FN, libmpi), MPI_Win_copy_attr_function)
@const_ref MPI_WIN_NULL_COPY_FN     MPI_Win_copy_attr_function   cglobal((:OMPI_C_MPI_WIN_NULL_COPY_FN, libmpi), MPI_Win_copy_attr_function)
@const_ref MPI_WIN_NULL_DELETE_FN   MPI_Win_delete_attr_function cglobal((:OMPI_C_MPI_WIN_NULL_DELETE_FN, libmpi), MPI_Win_delete_attr_function)


text:
module ConductanceInjection
export ConductanceStates, ConductanceParams, gen_syn_current

struct ConductanceStates{T <: AbstractFloat}
    g_ampa::T
    g_nmda::T
    g_gaba_a::T
    g_gaba_b::T
end

struct ConductanceParams{T <: AbstractFloat}
    tau_ampa::T
    tau_nmda::T
    tau_gaba_a::T
    tau_gaba_b::T
end

function gen_syn_current(dt,
                         delta_exct,
                         delta_inhbt,
                         states::ConductanceStates,
                         params::ConductanceParams,
                         state_v::AbstractFloat,
                         updater)

    states = ConductanceStates(agent.states.g_ampa * (1 - dt / agent.params.tau_ampa) + delta_exct,
                               agent.states.g_nmda * (1 - dt / agent.params.tau_nmda) + delta_exct,
                               agent.states.g_gaba_a * (1 - dt / agent.params.tau_gaba_a) + delta_inhbt,
                               agent.states.g_gaba_b * (1 - dt / agent.params.tau_gaba_b) + delta_inhbt)
    updater(states)
    return - (stts.g_ampa * state_v +
              stts.g_nmda * state_v * ((state_v + 80) / 60)^2 / (1 + ((state_v + 80) / 60)^2) +
              stts.g_gaba_a * (state_v + 70) + stts.g_gaba_b(state_v + 90))

end

end # Module end


text:
using Printf

@enum State empty tree fire

function evolution(nepoch::Int=100, init::Matrix{State}=fill(tree, 30, 50))
    # Single evolution
    function evolve!(forest::Matrix{State}; f::Float64=0.12, p::Float64=0.5)
        dir = [-1 -1; -1 0; -1 1; 0 -1; 0 1; 1 -1; 1 0; 1 1]
        # A tree will burn if at least one neighbor is burning
        for i in 1:size(forest, 1), j in 1:size(forest, 2)
            for k in 1:size(dir, 1)
                if checkbounds(Bool, forest, i + dir[k, 1], j + dir[k, 2]) &&
                    get(forest, i + dir[k, 1], j + dir[k, 2]) == fire
                    forest[i, j] = fire
                    break
                end
            end
        end
        for i in LinearIndices(forest)
            # A burning cell turns into an empty cell
            if forest[i] == fire forest[i] = empty end
            # A tree ignites with probability f even if no neighbor is burning
            if forest[i] == tree && rand() < f forest[i] = fire end
            # An empty space fills with a tree with probability p
            if forest[i] == empty && rand() < p forest[i] = tree end
        end
    end

    # Print functions
    function printforest(f::Matrix{State})
        for i in 1:size(f, 1)
            for j in 1:size(f, 2)
                print(f[i, j] == empty ? ' ' : f[i, j] == tree ? '🌲' : '🔥')
            end
            println()
        end
    end
    function printstats(f::Matrix{State})
        tot = length(f)
        nt  = count(x -> x in (tree, fire), f)
        nb  = count(x -> x == fire, f)
        @printf("\n%6i cell(s), %6i tree(s), %6i currently burning (%6.2f%%, %6.2f%%)\n",
                tot, nt, nb, nt / tot * 100, nb / nt * 100)
    end

    # Main
    printforest(init)
    printstats(init)
    for i in 1:nepoch
        # println("\33[2J")
        evolve!(init)
        # printforest(init)
        # printstats(init)
        # sleep(1)
    end
    printforest(init)
    printstats(init)
end

evolution()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
@testset "radial basis function" begin
    xc = rand(10)
    yc = rand(10)
    e = rand(10)
    c = rand(10)
    d = rand(3)
    r = RBF2D(xc, yc; c=c, eps=e, d=d, kind = 0)

    x = rand(5)
    y = rand(5)
    o = r(x, y)

    O = zeros(5)
    for i = 1:5
        for j = 1:10
            dist = sqrt((x[i]-xc[j])^2 + (y[i]-yc[j])^2)
            O[i] += c[j] * exp(-(e[j]*dist)^2)
        end
        O[i] += d[1] + d[2] * x[i] + d[3] * y[i]
    end
    init(sess)
    @test norm(run(sess, o)-O)<1e-5
end

@testset "interp1" begin
    Random.seed!(233)
    x = sort(rand(10))
    y = @. x^2 + 1.0
    z = [x[1]; x[2]; rand(5) * (x[end]-x[1]) .+ x[1]; x[end]]
    u = interp1(x,y,z)
    @test norm(run(sess, u)-[1.026422850882909
                    1.044414684090653
                    1.312604319732756
                    1.810845361128137
                    1.280789421523103
                    1.600084940795178
                    1.930560200260898
                    1.972130181835701]) < 1e-10
end


text:
function UniformIntegerDistribution(arg0::jint, arg1::jint)
    return UniformIntegerDistribution((jint, jint), arg0, arg1)
end

function cumulative_probability(obj::UniformIntegerDistribution, arg0::jint)
    return jcall(obj, "cumulativeProbability", jdouble, (jint,), arg0)
end

function get_numerical_mean(obj::UniformIntegerDistribution)
    return jcall(obj, "getNumericalMean", jdouble, ())
end

function get_numerical_variance(obj::UniformIntegerDistribution)
    return jcall(obj, "getNumericalVariance", jdouble, ())
end

function get_support_lower_bound(obj::UniformIntegerDistribution)
    return jcall(obj, "getSupportLowerBound", jint, ())
end

function get_support_upper_bound(obj::UniformIntegerDistribution)
    return jcall(obj, "getSupportUpperBound", jint, ())
end

function is_support_connected(obj::UniformIntegerDistribution)
    return jcall(obj, "isSupportConnected", jboolean, ())
end

function probability(obj::UniformIntegerDistribution, arg0::jint)
    return jcall(obj, "probability", jdouble, (jint,), arg0)
end


text:
abstract type QuiqboxContainer <: Any end

abstract type ParameterizedFunction <: Function end

abstract type QuiqboxVariableBox <: QuiqboxContainer end

abstract type QuiqboxParameter{ParameterT, ContainerT} <: QuiqboxVariableBox end
abstract type QuiqboxDataBox <: QuiqboxVariableBox end

abstract type ImmutableParameter{ContainerT, DataT} <: QuiqboxParameter{ImmutableParameter, ContainerT} end
abstract type SemiMutableParameter{ContainerT, DataT} <: QuiqboxParameter{SemiMutableParameter, ContainerT} end
abstract type MutableParameter{ContainerT, DataT} <: QuiqboxParameter{MutableParameter, ContainerT} end

abstract type ImmutableDataBox <: QuiqboxDataBox end
abstract type MutableDataBox <: QuiqboxDataBox end

abstract type AbstractHartreeFockFinalValue <: ImmutableDataBox end
abstract type AbstractBasisSetData <: ImmutableDataBox end
abstract type MolecularDataBox <: ImmutableDataBox end

abstract type HartreeFockintermediateData <: MutableDataBox end

abstract type DifferentiableParameter{ContainerT, DataT} <: MutableParameter{ContainerT, DataT} end

abstract type HartreeFockFinalValue{T} <: AbstractHartreeFockFinalValue end

abstract type BasisSetData{N} <: AbstractBasisSetData end

abstract type MolecularCoefficients <: MolecularDataBox end

abstract type QuiqboxBasis <: QuiqboxContainer end

abstract type SpatialOrbital <: QuiqboxBasis end
abstract type SpinOrbital <: QuiqboxBasis end

abstract type StructSpatialBasis <: SpatialOrbital end

abstract type AbstractMolOrbital <: SpinOrbital end

abstract type NucleusCenteredBasis <: StructSpatialBasis end
abstract type FloatingBasis <: StructSpatialBasis end

abstract type AbstractGaussFunc <: NucleusCenteredBasis end

abstract type AbstractGTBasisFuncs <: FloatingBasis end

abstract type MolecularHartreeFockCoefficient{NucleiC, ElectronC} <: MolecularCoefficients end

abstract type GTBasisFuncs{OrbitalN}  <: AbstractGTBasisFuncs end

abstract type CompositeGTBasisFuncs{NofLinearlyCombinedBasis, NofOrbital}  <: GTBasisFuncs{NofOrbital} end

abstract type FloatingGTBasisFuncs{Subshell, GaussFuncN, OrbitalN} <: CompositeGTBasisFuncs{1, OrbitalN} end


text:
using LinearAlgebra
using CompScienceMeshes
using SauterSchwab3D

const pI = point(1,5,3)
const pII = point(2,5,3)
const pIII = point(7,1,0)
const pIV = point(5,1,-3)
const pV = point(0,0,0)

const P = simplex(pI,pII,pIII)
const Q = simplex(pI,pIV,pV)

const sing = SauterSchwab3D.singularity_detection(P,Q)

Accuracy2 = 30
cv_ref = CommonVertex4D(sing,SauterSchwab3D._legendre(Accuracy2,0.0,1.0))

function integrand(x,y)
      return ((x-pI)'*(y-pV))*exp(-im*1*norm(x-y))/(4pi*norm(x-y))
end

function INTEGRAND(u,v)
   n1 = neighborhood(P,u)
   n2 = neighborhood(Q,v)
   x = cartesian(n1)
   y = cartesian(n2)
   output = integrand(x,y)*jacobian(n1)*jacobian(n2)
   return(output)
end

print("Ref: ")
ref = sauterschwab_parameterized(INTEGRAND, cv_ref)
println(ref)
println()

res_tp =[]
res_sp = []
res_gm = []
n1 = []
n2 = []
n3 = []
#Grauss tensor product
for i in 2:1:25
   Accuracy = i
   cv = CommonVertex4D(sing,SauterSchwab3D._legendre(Accuracy,0.0,1.0))

   int_tp = sauterschwab_parameterized(INTEGRAND, cv)

   num_pts = length(cv.qps)^4
   push!(n1,num_pts)

   push!(res_tp,int_tp)
end

#ShunnHam simplex product
for i in 2:1:7
   Accuracy = i
   cv_s = CommonVertex4D_S(sing,SauterSchwab3D._shunnham2D(Accuracy))

   int_sp = sauterschwab_parameterized(INTEGRAND, cv_s)

   num_pts = length(cv_s.qps)^2
   push!(n2,num_pts)
   push!(res_sp,int_sp)
end

#GrundmannMoeller simplex product
for i in 2:1:12
   Accuracy = i
   cv_gm = CommonVertex4D_S(sing,SauterSchwab3D._grundmannMoeller2D(2*Accuracy-1))

   int_gm = sauterschwab_parameterized(INTEGRAND, cv_gm)

   num_pts = length(cv_gm.qps)^2
   push!(n3,num_pts)
   push!(res_gm,int_gm)
end

err_tp = norm.(res_tp.-ref)/norm(ref)
err_sp = norm.(res_sp.-ref)/norm(ref)
err_gm = norm.(res_gm.-ref)/norm(ref)

println(err_tp)
println(err_sp)
println(err_gm)

using Plots
plot( yaxis=:log, xaxis=:log, fontfamily="Times")
plot!(n1,err_tp, label="Gauss Tensor-Product",markershape=:circle)
plot!(n2,err_sp, label="Simplex Tensor-Product",markershape=:rect)
#plot!(n3,err_gm, label="Simplex-Product GM",markershape=:x)
plot!(xlims=(5,10^4),ylims=(1e-6,1))
plot!(xlabel="Quad. points/Func. evals.", ylabel="Rel. Error.", title="Common Vertex 4D",legend=:bottomleft)

using BenchmarkTools

ref = -1.0761994366251377 + 1.0886007687110666im

cv = CommonVertex4D(sing,SauterSchwab3D._legendre(6,0.0,1.0))

int_tp = sauterschwab_parameterized(INTEGRAND, cv)

num_pts = length(cv.qps)^4
println(num_pts)
err_tp = norm.(int_tp.-ref)/norm(ref)
println(err_tp)
@time for i in 1:100 sauterschwab_parameterized(INTEGRAND, cv) end

cv_s = CommonVertex4D_S(sing,SauterSchwab3D._shunnham2D(6))

int_sp = sauterschwab_parameterized(INTEGRAND, cv_s)

num_pts = length(cv_s.qps)^2
println(num_pts)
err_sp = norm.(int_sp.-ref)/norm(ref)
println(err_sp)
 @time for i in 1:100 sauterschwab_parameterized(INTEGRAND, cv_s) end

#=

cv_gm = CommonVertex4D_S(SauterSchwab3D._grundmannMoeller2D(2*6-1))

int_gm = sauterschwab_parameterized(INTEGRAND, cv_gm)

num_pts = length(cv_gm.qps)^2
println(num_pts)
err_gm = norm.(int_gm.-ref)/norm(ref)
println(err_gm)
=#


text:
# --------------------------------------------------------------------------- #
# Lessons, aka job descriptions for training processes

const Region      = NTuple{4, Int}
const ImageOrFile = Union{Image, String}
const LabelOrFile = Union{Label, String}
const Selection   = Tuple{ImageOrFile, LabelOrFile, Region}

# Remark: Images and labels can either be given as data or as path.
# The selection region is *only* applied for data that is loaded via
# a path. So if an image or a label is given as data, the selection
# region is therefore ignored!
# In other words: If data is given, it is expected that the data
# was created with crop(data, region...)

mutable struct Lesson

  # Model
  model
  imgtype :: DataType
  kwargs  :: Vector{Any}

  # Algorithm
  optimizer :: String
  lr        :: Float64

  # Image and label data
  folder     :: String
  selections :: Array{Selection}

  # Data augmentation
  imageop :: ImageOp

  # Parameters
  epochs       :: Int
  batchsize    :: Int
  patchsize    :: Int
  patchmode    :: Int
  kernelsize   :: Int
  kernelheight :: Int

end

# --------------------------------------------------------------------------- #
# Lesson constructor

function Lesson(modelc;      # model or constructor function for a model
                imgtype      = RGBImage,
                kwargs       = Any[],
                folder       = "",
                optimizer    = "adam", # allowed: adam, rmsprop, nesterov
                lr           = 1e-4,
                selections   = [],
                imageop      = Id(),
                epochs       = 10,
                batchsize    = 1,
                patchsize    = 256,
                patchmode    = 0,
                kernelsize   = 7,
                kernelheight = 100)

    return Lesson(modelc, imgtype, kwargs, optimizer,
                  lr, folder, selections, imageop, epochs,
                  batchsize, patchsize, patchmode,
                  kernelsize, kernelheight)
end

# --------------------------------------------------------------------------- #
# Training with lessons

function resolvesel(s, folder)
  img, lbl, region = s
  if typeof(img) == String
    img = crop(imgload(joinpath(folder, img)), region...)
  else
    # If data is given, expect that it was created by calling
    # `crop` on the full image
    @assert imgsize(img) == region[4:-1:3]
  end
  if typeof(lbl) == String
    lbl = crop(lblload(joinpath(folder, img)), region...)
  else
    # If data is given, expect that it was created by calling
    # `crop` on the full label
    @assert length(lbl) == length(crop(lbl, 1, 1, region[3:4]...))
  end
  return (img, lbl)
end

function train(lesson :: Lesson; kwargs...)
  if typeof(lesson.model) <: Model
    model = deepcopy(lesson.model)
  else
    model = lesson.model(lesson.imgtype; lesson.kwargs...)
  end

  if lesson.optimizer == "adam"
    opt = Knet.Adam
  elseif lesson.optimizer == "rmsprop"
    opt = Knet.Rmsprop
  elseif lesson.optimizer == "nesterov"
    opt = Knet.Nesterov
  else
    warn("Optimizer $optimizer not valid. Using fallback optimizer adam")
    opt = Knet.Adam
  end

  lmgs = resolvesel.(lesson.selections, lesson.folder)

  train!(model, lmgs;
         epochs     = lesson.epochs,
         batchsize  = lesson.batchsize,
         patchsize  = lesson.patchsize,
         patchmode  = lesson.patchmode,
         kernelsize = lesson.kernelsize,
         peakheight = lesson.kernelheight,
         imageop    = lesson.imageop,
         lr         = lesson.lr,
         opt        = opt,
         shuffle    = true,
         kwargs...)

  return model
end

max_stars_repo_name: tscode/DCellC

max_issues_repo_name: tscode/DCellC


text:
using JuliaChem.JCModules

using HDF5

function compute_formation_energy(mol::Molecule, basis::Basis,
  total_energy::Float64)

  #== generate list of unique atoms in system ==#
  atom_list = Dict{String, Float64}([])
  h5open(joinpath(@__DIR__, "../../../records/eatom.h5"),"r") do eatom
    for atom in mol
      if !haskey(atom_list, atom.symbol)
        atom_list[atom.symbol] = read(eatom[atom.symbol][basis.model])
      end
    end
  end

  #== compute energy sum of constituent atoms ==#
  constituent_energy_sum = 0.0
  for atom in mol
    constituent_energy_sum += atom_list[atom.symbol]
  end

  #== compute formation energy ==#
  E_form = total_energy - constituent_energy_sum

  #== we are done ==#
  return E_form
end


text:
"""
$(TYPEDEF)

XY Euclidean manifold variable node softtype.
"""
struct Point2 <: IncrementalInference.InferenceVariable end
getDimension(::Point2) = 2
getManifolds(::Point2) = (:Euclid, :Euclid)


text:
### Hastings & Powell 1991
## Reproducing figure 5

#### Figure 5 (Poincaré section & map) ####
#### b1 = 3.0

u0 = [0.76, 0.16, 9.9]  #initial conditions for x,y,z
tspan = (0.0, 10000.0) #timespan
p_b3 = [5.0, 0.1, 3.0, 2.0, 0.4, 0.01] # parameter values (a1, a2, b1, b2, c1, c2)

# Solving the system of differential equations
prob_b3 = ODEProblem(parameterized_model!, u0, tspan, p_b3)
sol_b3 = solve(prob_b3, alg = RK4(), reltol=1e-14)

sol_x_b3 = sol_b3[1, :] # timeseries for x
sol_y_b3 = sol_b3[2, :] # timeseries for y
sol_z_b3 = sol_b3[3, :] # timeseries for z

constant_b3 = 9.0 # poincaré section (z value)
ϵ_b3 = 0.05 # margin of error for z value

# Find timesteps where points are in poincaré section
ind_x_b3 = findall(i -> 0.9 <= i <= 1.0, sol_x_b3) # x between 0.9 and 1 (high)
ind_y_b3 = findall(i -> 0.0 <= i <= 0.1, sol_y_b3) # y between 0 and 0.1 (low)
ind_z_b3 = findall(i -> constant_b3 - ϵ_b3 <= i <= constant_b3 + ϵ_b3, sol_z_b3) # z around constant

ind_xyz_b3 = zeros(Int64, 10000) # initialize the vector of concordant indices
# for loop : concordant indices
for j in 1:length(ind_x_b3)
    if any(i -> i == ind_x_b3[j], ind_y_b3) # concordant x and y indices
        if any(i -> i == ind_x_b3[j], ind_z_b3) # concordant x, y and z indices
            ind_xyz_b3[j] = ind_x_b3[j] # if the index if concordant, keep it
        end
    end
end

# for loop : non adjacent indices
for j in 1:(length(ind_xyz_b3) - 1)
    if ind_xyz_b3[j] == ind_xyz_b3[j + 1] - 1 # if 2 consecutive concordant indices are also
                                      # consecutive in the solution, don't keep them
        ind_xyz_b3[j] = 0
    end
end
ind_xyz_b3 = ind_xyz_b3[findall(i -> i > 0, ind_xyz_b3)] # remove extra values (0 values)

# Figure 5A - Poincaré section
xmin = 0.95 ; xmax = .983 ; ymin = 0.015 ; ymax = .04
fig5a = plot(sol_x_b3[ind_xyz_b3], sol_y_b3[ind_xyz_b3],
             legend = false,
             seriestype = :scatter,
             xlabel = "x(n)", xlim = (xmin, xmax), xticks = xmin:0.01:xmax, xminorticks = 2,
             ylabel = "y(n)", ylim = (ymin, ymax), yticks = ymin:0.01:ymax, yminorticks = 2,
             xgrid = :none, ygrid = :none,
             title = "A", titleloc = :left,
             c = :black, ms = 3)

# Figure 5B - Poincaré map
ind_n_b3 = ind_xyz_b3[1:length(ind_xyz_b3)-1] # indices for x(n)
ind_nplus1_b3 = ind_xyz_b3[2:length(ind_xyz_b3)] # indices for x(n+1)

xmin = 0.95 ; xmax = 0.98
fig5b = plot(sol_x_b3[ind_n_b3], sol_x_b3[ind_nplus1_b3],
             legend = false,
             seriestype = :scatter,
             xlabel = "x(n)", xlim = (xmin, xmax), xticks = xmin:0.01:xmax, xminorticks = 2,
             ylabel = "x(n+1)", ylim=(xmin, xmax), yticks = xmin:0.01:xmax, yminorticks = 2,
             xgrid = :none, ygrid = :none,
             title = "B", titleloc = :left,
             c = :black, ms = 3)
plot!(xmin:0.01:xmax, xmin:0.01:xmax, c = :black)

#### FIGURES 5C and 5D ####
#### b1 = 6.0

#Values
p_b6 = [5.0, 0.1, 6.0, 2.0, 0.4, 0.01] # parameter values (a1, a2, b1, b2, c1, c2)

# Solving the system of differential equations
prob_b6 = ODEProblem(parameterized_model!, u0, tspan, p_b6)
sol_b6 = solve(prob_b6, alg = RK4(), reltol=1e-14) # RK4 algorithm

sol_x_b6 = sol_b6[1, :] # timeseries for x
sol_y_b6 = sol_b6[2, :] # timeseries for y
sol_z_b6 = sol_b6[3, :] # timeseries for z

constant_b6 = 3.0 # poincaré section (z value)
ϵ_b6 = 0.05 # margin of error for z value

# Find timesteps where points are in poincaré section
ind_x_b6 = findall(i -> 0.93 <= i <= 1.0, sol_x_b6) # x between 0.9 and 1 (high)
ind_y_b6 = findall(i -> 0.0 <= i <= 0.085, sol_y_b6) # y between 0 and 0.1 (low)
ind_z_b6 = findall(i -> constant_b6 - ϵ_b6 <= i <= constant_b6 + ϵ_b6, sol_z_b6) # z around constant

ind_xyz_b6 = zeros(Int64, 15000) # initialize the vector of concordant indices
# for loop : concordant indices
for j in 1:length(ind_x_b6)
    if any(i -> i == ind_x_b6[j], ind_y_b6) # concordant x and y indices
        if any(i -> i == ind_x_b6[j], ind_z_b6) # concordant x, y and z indices
            ind_xyz_b6[j] = ind_x_b6[j] # if the index if concordant, keep it
        end
    end
end

# for loop : non adjacent indices
for j in 1:(length(ind_xyz_b6)-1)
    if ind_xyz_b6[j] == ind_xyz_b6[j + 1] - 1 # if 2 consecutive concordant indices are also
                                      # consecutive in the solution, don't keep them
        ind_xyz_b6[j] = 0
    end
end
ind_xyz_b6 = ind_xyz_b6[findall(i -> i > 0, ind_xyz_b6)] # remove extra values (0 values)

# Figure 5C - Poincaré section
xmin = 0.93 ; xmax = 1.003 ; ymin = -0.003 ; ymax = 0.09
fig5c = plot(sol_x_b6[ind_xyz_b6], sol_y_b6[ind_xyz_b6],
             legend = false,
             seriestype = :scatter,
             xlabel = "x(n)", xlim = (xmin, xmax), xticks = xmin:0.02:xmax, xminorticks = 2,
             ylabel = "y(n)", ylim = (ymin, ymax), yticks = 0.0:0.04:ymax, yminorticks = 2,
             xgrid = :none, ygrid = :none,
             title = "C", titleloc = :left,
             c = :black, ms = 3)

# Figure 5D - Poincaré map
ind_n_b6 = ind_xyz_b6[1:length(ind_xyz_b6)-1] # indices for x(n)
ind_nplus1_b6 = ind_xyz_b6[2:length(ind_xyz_b6)] # indices for x(n+1)

xmin = 0.93 ; xmax = 1.003
fig5d = plot(sol_x_b6[ind_n_b6], sol_x_b6[ind_nplus1_b6],
             legend = false,
             seriestype = :scatter,
             xlabel = "x(n)", xlim = (xmin, xmax), xticks = xmin:0.02:xmax, xminorticks = 2,
             ylabel = "x(n+1)", ylim = (xmin, xmax), yticks = xmin:0.02:xmax, yminorticks = 2,
             xgrid = :none, ygrid = :none,
             title = "D", titleloc = :left,
             c = :black, ms = 3)
plot!(xmin:0.01:xmax, xmin:0.01:xmax, color = :black)

# Export complete fig5
fig5 = plot(fig5a, fig5b, fig5c, fig5d, layout = 4,
            titlefontsize = 11, tickfontsize = 8, guidefontsize = 10,
            framestyle = :box
            )
savefig(fig5, joinpath("..", "article", "figures", "fig5.pdf"))


text:
@testset "shapes" begin
    inputsize = 20
    poollength = 10
    datalength = 100
    x = rand(Float32, inputsize, poollength, 1, datalength)
    m = TPALSTM(inputsize, 10, poollength, 2, 32, 1)
    @test size(m(x)) == (1, datalength)
    if Flux.CUDA.functional()
        @test size(gpu(m)(gpu(x))) == (1, datalength)
    end
    FluxArchitectures.initialize_bias!(m)
    @test all([all(m.lstm.chain.chain[i].cell.b .== 1) for i in 1:length(m.lstm.chain.chain)])
end

@testset "misc" begin
    @test repr(TPALSTM(10, 10, 10, 2, 32, 1)) == "TPALSTM(10, 10, 10, 2)"
    @test repr(TPALSTM(10, 10, 10, 2, 20, 4)) == "TPALSTM(10, 10, 10, 2, 20, 4)"
end

max_stars_repo_name: sdobber/FluxArchitectures

max_issues_repo_name: sdobber/FluxArchitectures

max_forks_repo_name: sdobber/FluxArchitectures


text:
"""
    erlang_rng(k, λ, shape=1; seed=nothing)

Generate a `shape` element array of random variables from a Erlang_{`k`}(`λ`) distribution. Optionally you can set a specific seed.

# Notes

The pdf of an Erlang_{k}(λ) distribution is given as:

``f(x, k, λ) = \\frac{λ^k e^{-λx} x^{k-1}}{(k-1)!} \\quad x ≥ 0``

# Examples

```julia-repl
julia> erlang_rng(5, .5)
1-element Vector{Float64}:
 10.803989701023117

julia> erlang_rng(3, 1, (2,2))
2×2×1 Array{Float64, 3}:
[:, :, 1] =
 2.19956  4.18505
 5.46892  2.5633
```

# References

D. Goldsman, P. Goldsman. A first course in probability and statistics. 2021.

L. Martino, D. Luengo. Extremely efficient generation of Gamma random variables for α ≥ 1. 2013.
"""
function erlang_rng(k::Int, λ::Real, shape::Union{Int, Tuple{Vararg{Int}}}=1; seed::Union{Int, Nothing}=nothing)
    U = get_std_uniform((shape..., k), seed=seed)
    k_dim = ndims(U)  # k is final dimension
    X = (-1/λ) .* log.(prod(U, dims=k_dim))  # Here (-1/λ) represents mean
    return X
end


text:
import Printf
import JSON
import Healpix
import LinearAlgebra: dot, ×
import Base: show

using RecipesBase

export MINUTES_PER_DAY,
       DAYS_PER_YEAR,
       rpm2angfreq,
       angfreq2rpm,
       period2rpm,
       rpm2period,
       ScanningStrategy,
       update_scanning_strategy,
       load_scanning_strategy,
       to_dict,
       save,
       time2pointing,
       time2pointing!,
       genpointings,
       genpointings!,
       CORE_SCANNING_STRATEGY,
       PICO_SCANNING_STRATEGY,
       SegmentedTimeSpan

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

rpm2angfreq(rpm) = 2π * rpm / 60
angfreq2rpm(ω) = ω / 2π * 60

@doc raw"""
    rpm2angfreq(rpm)
    angfreq2rpm(ω)

Convert rotations per minute into angular frequency 2πν (in Hertz), and vice
versa.
"""
rpm2angfreq, angfreq2rpm

period2rpm(p) = 60 / p
rpm2period(rpm) = 60 / rpm

@doc raw"""
    period2rpm(p)
    rpm2perriod(rpm)

Convert a period (time span) in seconds into a number of rotations per minute,
and vice versa.
"""
period2rpm, rpm2period

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

const MINUTES_PER_DAY = 60 * 24
const DAYS_PER_YEAR = 365.25

@doc raw"""
The structure `ScanningStrategy` encodes the information needed to build a
set of pointing directions. It contains the following fields:

- `omega_spin_hz`: angular speed of the rotation around the spin axis (equal to 2πν)
- `omega_prec_hz`: angular speed of the rotation around the precession axis (equal to 2πν)
- `omega_year_hz`: angular speed of the rotation around the Elliptic axis (equal to 2πν)
- `omega_hwp_hz`: angular speed of the rotation of the half-wave plate (equal to 2πν)
- `spinsunang_rad`: angle between the spin axis and the Sun-Earth direction
- `borespinang_rad`: angle between the boresight direction and the spin axis
- `q1`, `q3`: quaternions used to generate the pointings

Each field has its measure unit appended to the name. For instance, field
`spinsunang_rad` must be expressed in radians.

"""
struct ScanningStrategy
    omega_spin_hz::Float64
    omega_prec_hz::Float64
    omega_year_hz::Float64
    omega_hwp_hz::Float64
    spinsunang_rad::Float64
    borespinang_rad::Float64
    # First quaternion used in the rotation
    q1::Quaternion
    # Third quaternion used in the rotation
    q3::Quaternion

    ScanningStrategy(; spin_rpm = 0,
        prec_rpm = 0,
        yearly_rpm = 1  / (MINUTES_PER_DAY * DAYS_PER_YEAR),
        hwp_rpm = 0,
        spinsunang_rad = deg2rad(45.0),
        borespinang_rad = deg2rad(50.0),
    ) = new(rpm2angfreq(spin_rpm),
        rpm2angfreq(prec_rpm),
        rpm2angfreq(yearly_rpm),
        rpm2angfreq(hwp_rpm),
        spinsunang_rad,
        borespinang_rad,
        compute_q1(borespinang_rad),
        compute_q3(spinsunang_rad))

    ScanningStrategy(io::IO) = load_scanning_strategy(io)
    ScanningStrategy(filename::AbstractString) = load_scanning_strategy(filename)
end

compute_q1(borespinang_rad) = qrotation_y(borespinang_rad)
compute_q3(spinsunang_rad) = qrotation_y(π / 2 - spinsunang_rad)

@doc raw"""
    update_scanning_strategy(sstr::ScanningStrategy)

Update the internal fields of a `ScanningStrategy` object. If you change any of the
fields in a `ScanningStrategy` object after it has been created using the constructors,
call this function
before using one of the functions `time2pointing`, `time2pointing!`, `genpointings`,
and `genpointings!`, as they rely on a number of internal parameters that need to be
updated.

```julia
sstr = ScanningStrategy()
# ... use sstr ...

sstr.borespinang_rad *= 0.5
update_scanning_strategy(sstr)
```

"""
function update_scanning_strategy(sstr::ScanningStrategy)
    sstr.q1 = compute_q1(sstr.borespinang_rad)
    sstr.q3 = compute_q3(sstr.spinsunang_rad)
end

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

# «Exploring cosmic origins with CORE: Survey requirements and mission
# design», Delabrouille et al. (2018), Sect. 5.3
# https://doi.org/10.1088/1475-7516/2018/04/014

@doc raw"""
    CORE_SCANNING_STRATEGY::ScanningStrategy

A [`ScanningStrategy`](@ref) object describing the proposed scanning
strategy of the CORE spacecraft ([Delabrouille et al.,
2018](https://doi.org/10.1088/1475-7516/2018/04/014)).

"""
const CORE_SCANNING_STRATEGY = ScanningStrategy(
    spin_rpm = 1 / 2,
    prec_rpm = 1 / (4 * 24 * 60),      # Precession period is 4 days
    spinsunang_rad = deg2rad(30),
    borespinang_rad = deg2rad(65),
)

# «PICO: Probe of Inflation and Cosmic Origins», Hanany et al. (2019), Sect. 4.1.2
# https://ui.adsabs.harvard.edu/abs/2019arXiv190210541H/abstract

@doc raw"""
    PICO_SCANNING_STRATEGY::ScanningStrategy

A [`ScanningStrategy`](@ref) object describing the proposed scanning
strategy of the PICO spacecraft ([Hanany et al.,
2019](https://ui.adsabs.harvard.edu/abs/2019arXiv190210541H/abstract)).

"""
const PICO_SCANNING_STRATEGY = ScanningStrategy(
    spin_rpm = 1,
    prec_rpm = 1 / (10 * 60),      # Precession period is 10 hr
    spinsunang_rad = deg2rad(26),
    borespinang_rad = deg2rad(69),
)

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

function load_scanning_strategy(io::IO)
    data = JSON.Parser.parse(io)

    sstr_data = data["scanning_strategy"]

    ScanningStrategy(spin_rpm = sstr_data["spin_rpm"],
        prec_rpm = sstr_data["prec_rpm"],
        yearly_rpm = sstr_data["yearly_rpm"],
        hwp_rpm = sstr_data["hwp_rpm"],
        spinsunang_rad = sstr_data["spinsunang_rad"],
        borespinang_rad = sstr_data["borespinang_rad"])
end

function load_scanning_strategy(filename::AbstractString)
    open(filename) do inpf
        load_scanning_strategy(inpf)
    end
end

@doc raw"""
    load_scanning_strategy(io::IO) -> ScanningStrategy
    load_scanning_strategy(filename::AbstractString) -> ScanningStrategy

Create a `ScanningStrategy` object from the definition found in the JSON file
`io`, or from the JSON file with path `filename`. See also
[`load_scanning_strategy`](@ref).

"""
load_scanning_strategy

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

@doc raw"""
    to_dict(sstr::ScanningStrategy) -> Dict{String, Any}

Convert a scanning strategy into a dictionary suitable to be serialized using
JSON or any other structured format. See also [`save`](@ref).

"""
function to_dict(sstr::ScanningStrategy)
    Dict("scanning_strategy" => Dict("spin_rpm" => angfreq2rpm(sstr.omega_spin_hz),
        "prec_rpm" => angfreq2rpm(sstr.omega_prec_hz),
        "yearly_rpm" => angfreq2rpm(sstr.omega_year_hz),
        "hwp_rpm" => angfreq2rpm(sstr.omega_hwp_hz),
        "spinsunang_rad" => sstr.spinsunang_rad,
        "borespinang_rad" => sstr.borespinang_rad))
end

function save(io::IO, sstr::ScanningStrategy)
    print(io, JSON.json(to_dict(sstr), 4))
end

function save(filename::AbstractString, sstr::ScanningStrategy)
    open(filename, "w") do outf
        save(outf, sstr)
    end
end

@doc raw"""
    save(io::IO, sstr::ScanningStrategy)
    save(filename::AbstractString, sstr::ScanningStrategy)

Write a definition of the scanning strategy in a self-contained JSON file.
You can reload this definition using one of the constructors of
[`ScanningStrategy`](@ref).

"""
save

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

function show(io::IO, sstr::ScanningStrategy)
    Printf.@printf(io, """Scanning strategy:
    Spin angular velocity........................................ %g rot/s
    Precession angular velocity.................................. %g rot/s
    Yearly angular velocity around the Sun....................... %g rot/s
    Half-wave plate angular velocity............................. %g rot/s
    Angle between the spin axis and the Sun-Earth direction (α).. %f°
    Angle between the boresight direction and the spin axis (β).. %f°
""",
        sstr.omega_spin_hz,
        sstr.omega_prec_hz,
        sstr.omega_year_hz,
        sstr.omega_hwp_hz,
        rad2deg(sstr.spinsunang_rad),
        rad2deg(sstr.borespinang_rad))
end

function drawline(x1, y1, x2, y2; closepath = true)
    points = []
    append!(points, [(x1, y1), (x2, y2)])

    closepath && push!(points, (missing, missing))

    points
end

function drawarrow(x1, y1, x2, y2, head_length)
    points = []
    append!(points, [(x1, y1), (x2, y2)])

    head_aperture = π / 6
    angle = atan(y2 - y1, x2 - x1)
    sin1, cos1 = sincos(angle + π - head_aperture)
    sin2, cos2 = sincos(angle + π + head_aperture)
    append!(points, [(x2 + head_length * cos1, y2 + head_length * sin1),
                     (x2, y2),
                     (x2 + head_length * cos2, y2 + head_length * sin2)])

    push!(points, (missing, missing))

    points
end

function drawarc(center_x, center_y, radius, angle1, angle2; closepath = true)
    points = []
    append!(points, [(center_x + radius * cos(θ), center_y + radius * sin(θ))
                     for θ in range(angle1, angle2;
                                    length=round(Int, abs(angle2 - angle1) / 0.05))])

    closepath && push!(points, (missing, missing))

    points
end

function drawdoublearc(center_x, center_y, radius1, radius2, angle1, angle2)
    points = []

    append!(points, drawarc(center_x, center_y, radius1, angle1, angle2))
    append!(points, drawarc(center_x, center_y, radius2, angle1, angle2))

    points
end

function drawrotatedellipse(center_x, center_y, radius_x, radius_y, tilt_angle;
                            n = 40, closepath = true)

    sintilt, costilt = sincos(tilt_angle)
    rotmatr = [costilt -sintilt; sintilt costilt]

    points = []
    append!(points, [Float64[center_x, center_y] + rotmatr * Float64[radius_x * cos(θ), radius_y * sin(θ)]
                     for θ in range(0, 2π, length=n)])

    closepath && push!(points, (missing, missing))

    points
end

@recipe function plot_sstr(sstr::ScanningStrategy)
    α = sstr.spinsunang_rad
    β = sstr.borespinang_rad

    paths = []
    label_array = []

    arrow_head_size = 0.05

    #################################################################################
    # Sun-Earth axis

    append!(paths, drawarrow(-0.5, 0.0, 1.0, 0.0, arrow_head_size))

    #################################################################################
    # Spin-Sun angle

    append!(paths, drawarrow(0.0, 0.0, cos(α), sin(α), arrow_head_size))
    append!(paths, drawarc(0.0, 0.0, 0.10, 0.0, α))

    push!(label_array, (0.14 * cos(α/2), 0.14 * sin(α/2), "\\alpha"))

    #################################################################################
    # Boresight-Spin angle

    pt1 = Float64[cos(α + β), sin(α + β)]
    pt2 = Float64[cos(α - β), sin(α - β)]
    append!(paths, drawarrow(0.0, 0.0, pt1[1], pt1[2], arrow_head_size))
    append!(paths, drawarrow(0.0, 0.0, pt2[1], pt2[2], arrow_head_size))

    append!(paths, drawdoublearc(0.0, 0.0, 0.20, 0.22, α, α + β))
    push!(label_array, (0.26 * cos(α + β/2), 0.26 * sin(α + β/2), "\\beta"))

    append!(paths, drawdoublearc(0.0, 0.0, 0.21, 0.23, α - β, α))
    push!(label_array, (0.27 * cos(α - β/2), 0.27 * sin(α - β/2), "\\beta"))

    # Orbit of the spin axis around the Sun-Earth axis (vertical ellipse)
    @series begin
        orbits = []

        # Orbit of the spin axis around the Sun-Earth axis (vertical ellipse)
        append!(orbits, drawrotatedellipse(cos(α), 0.0, 0.1, sin(α), 0.0))

        seriestype := :path
        linecolor := :gray
        line := (0.5, :dash)

        [pt[1] for pt in orbits], [pt[2] for pt in orbits]
    end

    # Orbit of the boresight axis around the spin axis (tilted ellipse)
    @series begin
        orbits = []

        append!(orbits, drawrotatedellipse((pt1[1] + pt2[1]) / 2,
                                           (pt1[2] + pt2[2]) / 2,
                                           sqrt((pt2[1] - pt1[1])^2 + (pt2[2] - pt1[2])^2) / 2,
                                           0.1,
                                           π/2 + α))

        seriestype := :path
        linecolor := :gray
        line := (0.5, :dash)

        [pt[1] for pt in orbits], [pt[2] for pt in orbits]
    end

    seriestype := :path
    framestyle := :none
    axis := nothing
    legend := false
    linecolor := :black
    aspect_ratio := 1.0
    annotation := label_array

    [pt[1] for pt in paths], [pt[2] for pt in paths]
end

@doc raw"""
    plot(sstr::ScanningStrategy)

Plot a diagram representing the scanning strategy. You must have
loaded the [Plots](https://github.com/JuliaPlots/Plots.jl) package
before running this command.

"""
plot_scanning_strategy

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

function time2pointing!(sstr::ScanningStrategy, time_s, beam_dir, polangle_rad, resultvec)
    curpolang = mod2pi(polangle_rad + 4 * sstr.omega_hwp_hz * time_s)
    # The polarization vector lies on the XY plane; if polangle_rad=0 then
    # the vector points along the X direction at t=0.
    poldir = StaticArrays.SVector(cos(curpolang), sin(curpolang), 0.0)

    q2 = qrotation_z(sstr.omega_spin_hz * time_s)
    q4 = qrotation_x(sstr.omega_prec_hz * time_s)
    q5 = qrotation_z(sstr.omega_year_hz * time_s)

    qtot = q5 * (q4 * (sstr.q3 * (q2 * sstr.q1)))
    rot = rotationmatrix_normalized(qtot)
    # Direction in the sky of the beam main axis
    resultvec[4:6] = rot * beam_dir
    # Direction in the sky of the beam polarization axis
    poldir = rot * poldir

    # The North for a vector v is just -dv/dθ, as θ is the
    # colatitude and moves along the meridian
    (θ, ϕ) = Healpix.vec2ang(resultvec[4:6]...)
    northdir = StaticArrays.SVector(-cos(θ) * cos(ϕ), -cos(θ) * sin(ϕ), sin(θ))

    cosψ = clamp(dot(northdir, poldir), -1, 1)
    crosspr = northdir × poldir
    sinψ = clamp(sqrt(dot(crosspr, crosspr)), -1, 1)
    resultvec[3] = atan(cosψ, sinψ)

    resultvec[1], resultvec[2] = θ, ϕ
end

function time2pointing(sstr::ScanningStrategy, time_s, beam_dir, polangle_rad)
    resultvec = Float64[0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
    time2pointing!(sstr, time_s, beam_dir, polangle_rad, resultvec)
    resultvec
end

@doc raw"""
    time2pointing!(sstr::ScanningStrategy, time_s, beam_dir, polangle_rad, resultvec)
    time2pointing(sstr::ScanningStrategy, time_s, beam_dir, polangle_rad)

Calculate the pointing direction of a beam along the direction `beam_dir`, with
a detector sensitive to the polarization along the angle `polangle_rad`. The
result is saved in `resultvec` for `time2pointing!`, and it is the return value
of `time2pointing`; it is a 6-element array containing the following fields:

1. The colatitude (in radians) of the point in the sky
2. The longitude (in radians) of the point in the sky
3. The polarization angle, in the reference frame of the sky
4. The X component of the normalized pointing vector
5. The Y component
6. The Z component

Fields #4, #5, #6 are redundant, as they can be derived from the colatitude
(field #1) and longitude (field #2). They are returned as the code already
computes them.

The vector `beam_dir` and the angle `polangle_rad` must be expressed in the
reference system of the focal plane. If `polangle_rad == 0`, the detector
measures polarization along the x axis of the focal plane. The normal direction
to the focal plane is along the z axis; thus, the boresight director is such
that `beam_dir = [0., 0., 1.]`.

"""
time2pointing!, time2pointing

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

function genpointings!(sstr::ScanningStrategy, timerange_s, beam_dir, polangle_rad, result)
    @assert size(result)[1] == length(timerange_s)

    @inbounds for (idx, t) in enumerate(timerange_s)
        time2pointing!(sstr, t, beam_dir, polangle_rad, view(result, idx, :))
    end
end

function genpointings(sstr::ScanningStrategy, timerange_s, beam_dir, polangle_rad)
    result = Array{Float64}(undef, length(timerange_s), 6)
    genpointings!(sstr, timerange_s, beam_dir, polangle_rad, result)

    result
end

@doc raw"""
    genpointings!(sstr::ScanningStrategy, timerange_s, beam_dir, polangle_rad, result)
    genpointings(sstr::ScanningStrategy, timerange_s, beam_dir, polangle_rad)

Generate a set of pointing directions and angles for a given orientation
`beam_dir` (a 3-element vector) of the boresight beam, assuming the scanning
strategy in `sstr`. The pointing directions are calculated over all the elements
of the list `timerange_s`. The angle `polangle_rad` is the value of the
polarization angle at time ``t = 0``.

The two functions only differ in the way the result is returned to the caller.
Function `genpointings` returns a N×6 matrix containing the following fields:

1. The colatitude (in radians)
2. The longitude (in radians)
3. The polarization angle (in radians)
4. The X component of the one-length pointing vector
5. The Y component
6. The Z component

Function `genpointings!` works like `genpointings`, but it accept a
pre-allocated matrix as input (the `result` parameter) and will save the result
in it. The matrix must have two dimensions with size `(N, 6)` at least.

Both functions are simple iterators wrapping [`time2pointing!`](@ref) and
[`time2pointing`](@ref).

"""
genpointings!, genpointings

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

@doc raw"""
    struct SegmentedTimeSpan

An immutable structure representing a long time span, split into units
of equal length called *segments*. This structure is typically used to
make consecutive calls to [`genpointings!`](@ref) and
[`genpointings`](@ref).

The fields are the following:

- `start_time` is the time of the first sample in the time span
- `sampling_time` is the integration time for one sample
- `segment_duration` is the duration of one segment
- `num_of_segments` is the number of segments in the time span

The following example defines a time span of 1 year as the composition
of multiple spans, each lasting one day:

```julia
sts = SegmentedTimeSpan(
    start_time = 0.0,
    sampling_time = 0.1,
    segment_duration = 24 * 3600,
    num_of_segments = 365,
)
```

Once a `SegmentedTimeSpan` has been constructed, it behaves like an
array, whose elements are time ranges of type `StepRangeLen`. You can
iterate over it, to run functions like [`genpointings`](@ref):

```julia
for cur_time_span in sts
    pointings = genpointings(PICO_SCANNING_STRATEGY, cur_time_span, ...)

    # "pointings" contain the pointings for the time span
end
```

"""
struct SegmentedTimeSpan
    start_time::Float64
    sampling_time::Float64
    segment_duration::Float64
    num_of_segments::Int

    SegmentedTimeSpan(;
        start_time=0.0,
        sampling_time=1.0,
        segment_duration=3600.0,
        num_of_segments=24,
    ) = new(start_time, sampling_time, segment_duration, num_of_segments)
end

last_time_in_segment(start_time, sts::SegmentedTimeSpan) = start_time + sts.segment_duration - sts.sampling_time

function Base.iterate(iter::SegmentedTimeSpan)
    iter.num_of_segments > 0 || return nothing
    (iter.start_time : iter.sampling_time : last_time_in_segment(iter.start_time, iter),
     iter.start_time + iter.segment_duration)
end

function Base.iterate(iter::SegmentedTimeSpan, start_time)
    start_time >= iter.start_time + iter.segment_duration * iter.num_of_segments && return nothing

    (start_time : iter.sampling_time : last_time_in_segment(start_time, iter),
     start_time + iter.segment_duration)
end

function Base.getindex(sts::SegmentedTimeSpan, idx::Number)
    1 <= idx <= sts.num_of_segments || throw(BoundsError(sts, idx))

    start_time = sts.start_time + sts.segment_duration * (idx - 1)
    start_time : sts.sampling_time : last_time_in_segment(start_time, sts)
end

Base.getindex(sts::SegmentedTimeSpan, I) = [sts[x] for x in I]
Base.length(sts::SegmentedTimeSpan) = sts.num_of_segments
Base.firstindex(sts::SegmentedTimeSpan) = 1
Base.lastindex(sts::SegmentedTimeSpan) = length(sts)
Base.eachindex(sts::SegmentedTimeSpan) = 1:length(sts)


text:
@testset "bbox" begin

    line = LineString([[102.0, -10.0], [103.0, 1.0], [104.0, 0.0], [130.0, 4.0]])
    mp = MultiPolygon([[[[102.0, 2.0], [103.0, 2.0], [103.0, 3.0], [102.0, 3.0], [102.0, 2.0]]],
           [[[100.0, 0.0], [101.0, 0.0], [101.0, 1.0], [100.0, 1.0], [100.0, 0.0]], [[100.2, 0.2], [100.8, 0.2], [100.8, 0.8], [100.2, 0.8], [100.2, 0.2]]]])

    p1 = Point([102.0, 0.5])
    poly = Polygon([[[101.0, 0.0], [101.0, 1.0], [100.0, 1.0], [100.0, 0.0], [101.0, 0.0]]])
    ml = MultiLineString( [[[100.0, 0.0], [101.0, 1.0]],
    [[102.0, 2.0], [103.0, 3.0]]])

    fc = FeatureCollection([Feature(p1), Feature(mp), Feature(poly), Feature(line), Feature(ml)])

    @test bbox(fc) == [100, -10, 130, 4]
    @test bbox(p1) == [102, 0.5, 102, 0.5]
    @test bbox(line) == [102, -10, 130, 4]
    @test bbox(poly) == [100, 0, 101, 1]
    @test bbox(ml) == [100, 0, 103, 3]
    @test bbox(mp) == [100, 0, 103, 3]

    @test bbox(Feature(p1)) == bbox(p1)

    bbpoly = [0, 0, 10, 10]

    coords = bbox_polygon(bbpoly).coordinates
    @test length(coords[1]) === 5
    @test coords[1][1][1] == coords[1][end][1]
    @test coords[1][1][2] == coords[1][1][2]

    @test_throws ErrorException bbox_polygon([1, 2, 3, 4, 5, 6])

end