text:
"""
    ExpMethodHigham2005(A::AbstractMatrix);
    ExpMethodHigham2005(b::Bool=true);

Computes the matrix exponential using the algorithm Higham, N. J. (2005). "The scaling and squaring method for the matrix exponential revisited." SIAM J. Matrix Anal. Appl.Vol. 26, No. 4, pp. 1179–1193" based on generated code. If a matrix is specified, balancing is determined automatically.
"""
struct ExpMethodHigham2005
    do_balancing::Bool
end
ExpMethodHigham2005(A::AbstractMatrix)=ExpMethodHigham2005(A isa StridedMatrix)
ExpMethodHigham2005()=ExpMethodHigham2005(true)

function alloc_mem(A,::ExpMethodHigham2005)
    T = eltype(A)
    scale = T <: LinearAlgebra.BlasFloat ? similar(A, real(T), size(A, 1)) : nothing
    return [similar(A) for i=1:5], scale
end

# Import the generated code
for i=1:13
    include("exp_generated/exp_$i.jl")
end

function getmem(cache,k) # Called from generated code
    return cache[k-1];
end
function ldiv_for_generated!(C,A,B) # C=A\B. Called from generated code
    F = lu!(A) # This allocation is unavoidable, due to the interface of LinearAlgebra
    ldiv!(F,B); # Result stored in B
    if (pointer_from_objref(C) != pointer_from_objref(B)) # Aliasing allowed
        copyto!(C,B)
    end
    return C
end

const RHO_V = (0.015, 0.25, 0.95, 2.1, 5.4, 10.8, 21.6, 43.2, 86.4, 172.8, 345.6, 691.2)

# From LinearAlgebra
const LIBLAPACK = VERSION >= v"1.7" ? BLAS.libblastrampoline : LAPACK.liblapack
using LinearAlgebra: BlasInt, checksquare
for (gebal, gebak, elty, relty) in
    ((:dgebal_, :dgebak_, :Float64, :Float64),
     (:sgebal_, :sgebak_, :Float32, :Float32),
     (:zgebal_, :zgebak_, :ComplexF64, :Float64),
     (:cgebal_, :cgebak_, :ComplexF32, :Float32))
    @eval begin
        #     SUBROUTINE DGEBAL( JOB, N, A, LDA, ILO, IHI, SCALE, INFO )
        #*     .. Scalar Arguments ..
        #      CHARACTER          JOB
        #      INTEGER            IHI, ILP, INFO, LDA, N
        #     .. Array Arguments ..
        #      DOUBLE PRECISION   A( LDA, * ), SCALE( * )
        function gebal_noalloc!(job::AbstractChar, A::AbstractMatrix{$elty}, scale)
            BLAS.chkstride1(A)
            n = checksquare(A)
            LAPACK.chkfinite(A) # balancing routines don't support NaNs and Infs
            ihi = Ref{BlasInt}()
            ilo = Ref{BlasInt}()
            info = Ref{BlasInt}()
            ccall((BLAS.@blasfunc($gebal), LIBLAPACK), Cvoid,
                  (Ref{UInt8}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
                   Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$relty}, Ptr{BlasInt}, Clong),
                  job, n, A, max(1,stride(A,2)), ilo, ihi, scale, info, 1)
            LAPACK.chklapackerror(info[])
            ilo[], ihi[], scale
        end
    end
end

# Inplace add of a UniformScaling object (support julia 1.6.2)
@inline function inplace_add!(A,B::UniformScaling) # Called from generated code
    s = B.λ
    @inbounds for i in diagind(A)
        A[i] += s
    end
end
function exponential!(A,method::ExpMethodHigham2005,_cache=alloc_mem(A,method))
    cache, _scale = _cache
    n = LinearAlgebra.checksquare(A)
    nA = opnorm(A,1);

    # Maybe to balancing
    if method.do_balancing
        if A isa StridedMatrix{<:LinearAlgebra.BLAS.BlasFloat}
            ilo, ihi, scale = gebal_noalloc!('B', A, _scale)    # modifies A and _scale
        else
            A, bal = GenericSchur.balance!(A)
            ilo, ihi, scale = bal.ilo, bal.ihi, bal.D
        end
    end

    # Make the call to the appropriate exp_gen! function
    X = Base.Cartesian.@nif 13 d -> begin
        nA < RHO_V[d]
    end d -> begin # if condition
        exp_gen!(cache, A, Val(d))
    end d -> begin # fallback (d == 13)
        exp_gen!(cache, A, Val(d))
    end

    # Undo the balancing
    if method.do_balancing
        for j = ilo:ihi
            scj = scale[j]
            for i = 1:n
                X[j,i] *= scj
            end
            for i = 1:n
                X[i,j] /= scj
            end
        end

        if A isa StridedMatrix{<:LinearAlgebra.BLAS.BlasFloat}
            if ilo > 1       # apply lower permutations in reverse order
                for j in (ilo-1):-1:1; LinearAlgebra.rcswap!(j, Int(scale[j]), X) end
            end
            if ihi < n       # apply upper permutations in forward order
                for j in (ihi+1):n;    LinearAlgebra.rcswap!(j, Int(scale[j]), X) end
            end
        else
            if ilo > 1       # apply lower permutations in reverse order
                for j in (ilo-1):-1:1; LinearAlgebra.rcswap!(j, bal.prow[j], X) end
            end
            if ihi < n       # apply upper permutations in forward order
                for j in (ihi+1):n;    LinearAlgebra.rcswap!(j, bal.pcol[j-ihi], X) end
            end
        end
    end

    return X
end

max_stars_repo_name: JuliaDiffEq/ExponentialUtilities

max_issues_repo_name: JuliaDiffEq/ExponentialUtilities

max_forks_repo_name: JuliaDiffEq/ExponentialUtilities


text:
using Lorenz96
using LogFixPoint16s
using PyPlot
using JLD2, FileIO
using ColorSchemes
using Statistics, StatsBase
cmap = ColorMap(ColorSchemes.oslo.colors)

# LogFixPoint16s.set_nfrac(10)

N = 300
n = 200
Δt = 0.1
F = 8.0

t = collect(0:N)*Δt
nvec = collect(1:n)

αs = [1.0,1.3,2.0]
# Xf1 = L96(Float32;N,n,Δt,α=αs[1])
# Xf2 = L96(Float32;N,n,Δt,α=αs[2])
# Xf3 = L96(Float32;N,n,Δt,α=αs[3])

X  = fill(F,n)
X[100] = F+0.1

Xa1 = L96(LogFixPoint16,X=X/αs[1];N,n,Δt,α=αs[1])
Xa2 = L96(LogFixPoint16,X=X/αs[2];N,n,Δt,α=αs[2])
Xa3 = L96(LogFixPoint16,X=X/αs[3];N,n,Δt,α=αs[3])

## ENSEMBLE

# parameters
Nens = 100
N2 = 10
α2 = collect(1:1e-3:2.0)
nαs = length(α2)
E = 5.13    # climatological error

f16,lfx16 = load("N10n200.jld2","f16","lfx16")
bf16,blfx16 = load("BF_N10n200.jld2","f16","lfx16")

# monte carlo sampling for mean
m = Nens
f16m = fill(0.0,nαs,m)
lfx16m = fill(0.0,nαs,m)
bf16m = fill(0.0,nαs,m)
blfx16m = fill(0.0,nαs,m)

for i in 1:m
    # calculate median from a random half of the ensemble members
    f16m[:,i] = mean(f16[:,rand(1:Nens,Nens÷2)],dims=2)
    lfx16m[:,i] = mean(lfx16[:,rand(1:Nens,Nens÷2)],dims=2)
    bf16m[:,i] = mean(bf16[:,rand(1:Nens,Nens÷2)],dims=2)
    blfx16m[:,i] = mean(blfx16[:,rand(1:Nens,Nens÷2)],dims=2)
end

# confidence interval from those
f16lo = [percentile(f16m[i,:],2.5) for i in 1:nαs]
f16hi = [percentile(f16m[i,:],97.5) for i in 1:nαs]

lfx16lo = [percentile(lfx16m[i,:],2.5) for i in 1:nαs]
lfx16hi = [percentile(lfx16m[i,:],97.5) for i in 1:nαs]

bf16lo = [percentile(bf16m[i,:],2.5) for i in 1:nαs]
bf16hi = [percentile(bf16m[i,:],97.5) for i in 1:nαs]

blfx16lo = [percentile(blfx16m[i,:],2.5) for i in 1:nαs]
blfx16hi = [percentile(blfx16m[i,:],97.5) for i in 1:nαs]

## PLOT

ioff()
fig,axs = subplots(2,3,sharey=true,figsize=(10,6))
tight_layout(rect=[0.03,0.04,0.95,0.96],w_pad=0.01)
pos = axs[1,end].get_position()
pos1 = axs[1,1].get_position()
cax = fig.add_axes([pos.x1+0.02,pos.y0,0.015,pos.y1-pos.y0])
bax = fig.add_axes([pos1.x0,0.1,pos.x1-pos1.x0,0.95*(pos.y1-pos.y0)])
axs[2,1].remove()
axs[2,2].remove()
axs[2,3].remove()

vmin,vmax = -8,8
axs[1,1].pcolormesh(t,nvec,Xa1;vmin,vmax,cmap)
axs[1,2].pcolormesh(t,nvec,Xa2;vmin,vmax,cmap)
q = axs[1,3].pcolormesh(t,nvec,Xa3;vmin,vmax,cmap)
cbar = colorbar(q;cax)
cbar.set_ticks([-8,-4,0,4,8])

axs[1,1].set_ylabel("X")
axs[1,1].set_yticks([])

l = ["a","b","c"]
for (ia,ax) in enumerate(axs[1,:])
    ax.set_title("LogFixPoint16, α=$(αs[ia])",loc="left")
    ax.set_title(l[ia],loc="right",fontweight="bold")
    ax.set_xticks([0,10,20,30])
    ax.set_xticklabels([0,10,20,30])
end

axs[1,2].set_xlabel("time [mtu]")
axs[1,3].set_xlabel("time [mtu]")

bax.semilogy(α2,mean(f16,dims=2)/E,"k",lw=2.5,label="Float16")
bax.plot(α2,mean(lfx16,dims=2)/E,"C1",lw=2.5,label="LogFixPoint16")
bax.plot(α2,mean(bf16,dims=2)/E,"grey",lw=2.5,label="BFloat16")
bax.plot(α2,mean(blfx16,dims=2)/E,"C3",lw=2.5,label="BLogFixPoint16")

alfa=0.5
bax.fill_between(α2,f16lo/E,f16hi/E,alpha=alfa,color="k",label="95% confidence")
bax.fill_between(α2,lfx16lo/E,lfx16hi/E,alpha=alfa,color="C1",label="95% confidence")
bax.fill_between(α2,bf16lo/E,bf16hi/E,alpha=alfa,color="grey",label="95% confidence")
bax.fill_between(α2,blfx16lo/E,blfx16hi/E,alpha=alfa,color="C3",label="95% confidence")

bax.legend(loc=1,ncol=2,fontsize=9)
bax.set_ylim(5e-4,6e-2)

bax.set_ylabel("RMSE (normalised)")
bax.set_xlabel("Friction parameter α")
bax.set_xlim(1.0,2.0)

bax.set_title("Error to Float64 after 1 mtu",loc="left")
bax.set_title("d",loc="right",fontweight="bold")

savefig("lorenz96_lfx16.png",dpi=200)
close(fig)

ext: jl


text: any_of(cols::AbstractVector{T}) where T = df -> intersect(T.(names(df)), cols)


text:
# using ModelExploration
include("../src/Slice.jl")
include("../src/ModelExploration.jl")

import AlgebraicPetri
using Test
using Catlab.CategoricalAlgebra
using Catlab.Graphs
using Catlab.Graphs.SimpleGraphs: TheoryGraph

const Petri = AlgebraicPetri.PetriNet

"""
Loose ACSet Transformations

Make a wrapper for Slice and implement limits for that.
"""

# Petri net example
###################
function dezero(s::StructACSet)::StructACSet
  for h in keys(s.homs)
    set_subpart!(s, h, replace(x -> x == 0 ? 1 : x, s[h]))
  end
  return s
end

function graph_to_petri(G::Graph)::Slice
  idE = id(TheoryGraph[:E])
  obs = Dict(:S => :V, :T=>:E, :I => :E, :O=>:E)
  morphs = Dict(:is => :src, :it => idE, :os => :src, :ot => idE)
  return Slice(migrate(PetriNet, G, obs, morphs), A, [:S=>[1=>1]])
end

# Slices over A are PetriNets that have only unary and binary transitions
UnaryT = dezero(@acset Petri begin S = 1; T = 1; I = 1; O = 1 end)
BinaryT = dezero(@acset Petri begin S = 1; T = 1; I = 2; O = 2 end)
S1 = @acset Petri begin S=1 end
A = apex(pushout(homomorphism(S1, UnaryT), homomorphism(S1, BinaryT)))

"""
Create linear Petri Net
Could we do this with functor from Graph?
"""
multi_city(n::Int)::Slice = graph_to_petri(complete_graph(Graph, n))

MultiCity = Linear(multi_city)

young_old, rich_poor = multi_city(2), multi_city(2)
young_mid_old = multi_city(3)
young_old_rich_poor = multi_city(4)

AgeRefine = Finite([young_old, young_mid_old, rich_poor, young_old_rich_poor],
                    [[2,4],    Int[],         [4],       Int[]])

exp = Explore(PetriNet, [Slice(A)], SIR, [Free(), MultiCity, AgeRefine])

data = DataSet(path)

exp.run(data)


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct FunctionDef_AttrEntry <: ProtoType
    key::AbstractString
    value::AttrValue
    FunctionDef_AttrEntry(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct FunctionDef_AttrEntry (mapentry)

mutable struct FunctionDef_RetEntry <: ProtoType
    key::AbstractString
    value::AbstractString
    FunctionDef_RetEntry(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct FunctionDef_RetEntry (mapentry)

mutable struct FunctionDef_ControlRetEntry <: ProtoType
    key::AbstractString
    value::AbstractString
    FunctionDef_ControlRetEntry(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct FunctionDef_ControlRetEntry (mapentry)

mutable struct FunctionDef <: ProtoType
    signature::OpDef
    attr::Base.Dict{AbstractString,AttrValue} # map entry
    node_def::Base.Vector{NodeDef}
    ret::Base.Dict{AbstractString,AbstractString} # map entry
    control_ret::Base.Dict{AbstractString,AbstractString} # map entry
    FunctionDef(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct FunctionDef
const __fnum_FunctionDef = Int[1,5,3,4,6]
meta(t::Type{FunctionDef}) = meta(t, ProtoBuf.DEF_REQ, __fnum_FunctionDef, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct GradientDef <: ProtoType
    function_name::AbstractString
    gradient_func::AbstractString
    GradientDef(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GradientDef

mutable struct FunctionDefLibrary <: ProtoType
    _function::Base.Vector{FunctionDef}
    gradient::Base.Vector{GradientDef}
    FunctionDefLibrary(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct FunctionDefLibrary

export FunctionDefLibrary, FunctionDef_AttrEntry, FunctionDef_RetEntry, FunctionDef_ControlRetEntry, FunctionDef, GradientDef
# mapentries: "FunctionDef_AttrEntry" => ("AbstractString", "AttrValue"), "FunctionDef_RetEntry" => ("AbstractString", "AbstractString"), "FunctionDef_ControlRetEntry" => ("AbstractString", "AbstractString")


text:
## top-level submission file

"""
    f(a, b) = c
A function that adds the inputs `a` and `b`
"""
function f(a, b)
    # write your code here. Remember to return the sum of a and b!
end


text:
using FluxRC, Test

# Kou's benchmark results
"""
r
[-0.81684757,  0.63369515, -0.81684757, -0.10810302, -0.78379396, -0.10810302]

s
[0.63369515, -0.81684757, -0.81684757, -0.78379396, -0.10810302, -0.10810302]

rf
[-0.77459667, 0., 0.77459667, 0.77459667, 0., -0.77459667, -1., -1., -1.]

sf
[-1., -1., -1., -0.77459667, 0., 0.77459667, 0.77459667, 0., -0.77459667]

V
[[ 0.70710678,  1.45054272,  1.39329302,  0.        ,  0.        ,
    -0.04593312],
[ 0.70710678, -0.72527136,  0.43019448,  1.25620684, -0.83406781,
    1.03084378],
[ 0.70710678, -0.72527136,  0.43019448, -1.25620684,  0.83406781,
    1.03084378],
[ 0.70710678, -0.67569095,  0.30868284,  0.        , -0.        ,
    -1.08925616],
[ 0.70710678,  0.33784547, -0.70898933, -0.58516552, -0.88132995,
    0.04853591],
[ 0.70710678,  0.33784547, -0.70898933,  0.58516552,  0.88132995,
    0.04853591]]

correction_coeffs
[[ 0.3928371   0.62853936  0.3928371   0.55555556  0.88888889  0.55555556
  0.3928371   0.62853936  0.3928371 ]
[-0.55555556 -0.88888889 -0.55555556 -0.52003383  0.62853936  1.30570803
  0.923275    0.44444444 -0.36771945]
[ 0.68041382  1.08866211  0.68041382  0.21689446 -0.76980036  1.70760644
  1.20746009 -0.54433105  0.15336754]
[-0.74535599  0.          0.74535599  1.20746009  1.08866211  0.15336754
 -0.10844723 -0.76980036 -0.85380322]
[ 0.91287093 -0.         -0.91287093 -0.64549722  2.          0.64549722
 -0.45643546 -1.41421356  0.45643546]
[ 0.60858062 -1.21716124  0.60858062  1.6939963   0.86066297  0.02732963
  0.01932497  0.60858062  1.19783627]]

phifj_solution
[[ 0.39198253  0.72780648  0.39198253 -0.13710719  0.42817214  4.66476318
  3.29848568  0.30276343 -0.09694943]
[-0.09694943  0.30276343  3.29848568  4.66476318  0.42817214 -0.13710719
  0.39198253  0.72780648  0.39198253]
[ 3.29848568  0.30276343 -0.09694943  0.55434701  1.02927379  0.55434701
 -0.09694943  0.30276343  3.29848568]
[ 0.20029351  2.7069103   0.20029351 -1.03402506 -0.97126559  0.00792184
  0.00560158 -0.68678848 -0.73116613]
[-0.73116613 -0.68678848  0.00560158  0.00792184 -0.97126559 -1.03402506
  0.20029351  2.7069103   0.20029351]
[ 0.00560158 -0.68678848 -0.73116613  0.2832578   3.82814926  0.2832578
 -0.73116613 -0.68678848  0.00560158]]
"""

py_V = [
    0.70710678 1.45054272 1.39329302 0.0 0.0 -0.04593312
    0.70710678 -0.72527136 0.43019448 1.25620684 -0.83406781 1.03084378
    0.70710678 -0.72527136 0.43019448 -1.25620684 0.83406781 1.03084378
    0.70710678 -0.67569095 0.30868284 0.0 -0.0 -1.08925616
    0.70710678 0.33784547 -0.70898933 -0.58516552 -0.88132995 0.04853591
    0.70710678 0.33784547 -0.70898933 0.58516552 0.88132995 0.04853591
]

py_Vf = [
    0.70710678 -1.0 1.22474487 -1.34164079 1.64316767 1.09544512
    0.70710678 -1.0 1.22474487 0.0 -0.0 -1.36930639
    0.70710678 -1.0 1.22474487 1.34164079 -1.64316767 1.09544512
    0.70710678 -0.661895 0.27606157 1.5368458 -0.82158384 2.15610529
    0.70710678 0.5 -0.61237244 0.8660254 1.59099026 0.6846532
    0.70710678 1.661895 2.17342817 0.19520501 0.82158384 0.03478494
    0.70710678 1.661895 2.17342817 -0.19520501 -0.82158384 0.03478494
    0.70710678 0.5 -0.61237244 -0.8660254 -1.59099026 0.6846532
    0.70710678 -0.661895 0.27606157 -1.5368458 0.82158384 2.15610529
]

phifj_ref = [
    0.39198253 0.72780648 0.39198253 -0.13710719 0.42817214 4.66476318 3.29848568 0.30276343 -0.09694943
    -0.09694943 0.30276343 3.29848568 4.66476318 0.42817214 -0.13710719 0.39198253 0.72780648 0.39198253
    3.29848568 0.30276343 -0.09694943 0.55434701 1.02927379 0.55434701 -0.09694943 0.30276343 3.29848568
    0.20029351 2.7069103 0.20029351 -1.03402506 -0.97126559 0.00792184 0.00560158 -0.68678848 -0.73116613
    -0.73116613 -0.68678848 0.00560158 0.00792184 -0.97126559 -1.03402506 0.20029351 2.7069103 0.20029351
    0.00560158 -0.68678848 -0.73116613 0.2832578 3.82814926 0.2832578 -0.73116613 -0.68678848 0.00560158
]

# my workflow
N = deg = 2
Np = (N + 1) * (N + 2) ÷ 2

pl, wl = tri_quadrature(N)

V = vandermonde_matrix(N, pl[:, 1], pl[:, 2])
Vr, Vs = ∂vandermonde_matrix(N, pl[:, 1], pl[:, 2])
∂l = ∂lagrange(V, Vr, Vs)

ϕ = correction_field(N, V)

pf, wf = triface_quadrature(N)
ψf = zeros(3, N + 1, Np)
for i = 1:3
    ψf[i, :, :] .= vandermonde_matrix(N, pf[i, :, 1], pf[i, :, 2])
end

# Vandermonde -> solution points
V1 = hcat(V[1, :], V[3, :], V[5, :], V[2, :], V[4, :], V[6, :]) |> permutedims
@test V1 ≈ py_V

# Vandermonde -> flux points
Vf = zeros(9, 6)
for i = 1:3, j = 1:3
    Vf[(i-1)*3+j, :] .= ψf[i, j, :]
end
@test Vf ≈ py_Vf

# coefficients
ϕ1 = zeros(6, 9)
for i = 1:3, j = 1:3
    ϕ1[:, (i-1)*3+j] .= ϕ[i, j, :]
end
ϕ2 = hcat(ϕ1[1, :], ϕ1[3, :], ϕ1[5, :], ϕ1[2, :], ϕ1[4, :], ϕ1[6, :]) |> permutedims
@test ϕ2 ≈ phifj_ref

# reproduce from python script
cd(@__DIR__)
include("python_phi.jl")

py_sigma, py_fi, py_figrad =
    py"get_phifj_solution_grad_tri"(N, N + 1, 6, 3 * (N + 1), N + 1, py_V, py_Vf)
@test ϕ2 ≈ py_fi


text:
function updateyaml(mech, p)
    yaml = YAML.load_file(mech)
    n_species = length(yaml["phases"][1]["species"])
    n_reactions = length(yaml["reactions"])
    species_names = yaml["phases"][1]["species"]
    elements = yaml["phases"][1]["elements"]
    n_elements = length(elements)

    _p = reshape(p, npr, 3)

    for i = 1:npr
        p_vec = _p[i, :]
        fA = exp(p_vec[1])
        fb = p_vec[2]
        fEa = p_vec[3]
        reaction = yaml["reactions"][i]
        if haskey(reaction, "type")
            if (reaction["type"] == "falloff")
                yaml["reactions"][i]["low-P-rate-constant"]["A"] *= fA
                yaml["reactions"][i]["low-P-rate-constant"]["b"] += fb
                yaml["reactions"][i]["low-P-rate-constant"]["Ea"] += fEa
                yaml["reactions"][i]["high-P-rate-constant"]["A"] *= fA
                yaml["reactions"][i]["high-P-rate-constant"]["b"] += fb
                yaml["reactions"][i]["high-P-rate-constant"]["Ea"] += fEa
            else
                yaml["reactions"][i]["rate-constant"]["A"] *= fA
                yaml["reactions"][i]["rate-constant"]["b"] += fb
                yaml["reactions"][i]["rate-constant"]["Ea"] += fEa
            end
        else
            yaml["reactions"][i]["rate-constant"]["A"] *= fA
            yaml["reactions"][i]["rate-constant"]["b"] += fb
            yaml["reactions"][i]["rate-constant"]["Ea"] += fEa * 1000.0
        end
    end

    YAML.write_file(mech[1:end - 5] * "_op.yaml", yaml)
end

pyplot()

epoch = 58
@load string(ckpt_path, "/model_$(epoch).bson") p opt l_loss l_loss_val l_grad l_pnorm iter
# @load string(ckpt_path, "/mymodel.bson") p opt l_loss l_loss_val l_grad l_pnorm iter

# @load string("./results/JetA/checkpoint/model_103.bson") p opt l_loss l_loss_val l_grad l_pnorm iter
# p[33:48] .= p[33:48]/1000.0

writedlm(string(fig_path, "loss.txt"), hcat([l_loss, l_loss_val, l_grad, l_pnorm]...))

updateyaml(mech, p)

plt_loss = plot(l_loss, yscale=:log10, label="train", lw=2);
plot!(plt_loss, l_loss_val, yscale=:log10, label="validation", lw=2, ls=:dash);
plt_grad = plot(l_grad, yscale=:log10, label="grad_norm");
plt_pnorm = plot(l_pnorm .+ 1.e-8, yscale=:identity, label="p_norm");
xlabel!(plt_loss, "Epoch");
xlabel!(plt_grad, "Epoch");
xlabel!(plt_pnorm, "Epoch");
ylabel!(plt_loss, "Loss");
ylabel!(plt_grad, "Grad Norm");
ylabel!(plt_pnorm, "p Norm");
plt = plot(plt_loss, legend=:best);
plot!(
        plt,
        xtickfontsize = 11,
        ytickfontsize = 11,
        xguidefontsize = 12,
        yguidefontsize = 12,
        size = (350, 350),
    )
png(plt, string(fig_path, "/loss"));

plt = plot(plt_pnorm, legend=:best);
plot!(
        plt,
        xtickfontsize = 11,
        ytickfontsize = 11,
        xguidefontsize = 12,
        yguidefontsize = 12,
        size = (350, 350),
    )
png(plt, string(fig_path, "/pnorm"));

p_ga = readdlm(string(ckpt_path, "/../p_ga.txt"), '\t')
p_ga = reshape(p_ga, npr*3, 1)

l_phi = [0.5, 1.0, 1.3] .* 1.08831169
l_plt = []
for i in 1:3
    expdata = l_exp[i]
    l_idt = zeros(size(expdata)[1], 5)
    for j in 1:size(expdata)[1]
        local T0 = 1000.0 / expdata[j, 1]
        l_idt[j, 1] = expdata[j, 1]
        l_idt[j, 2] = expdata[j, 2]
        local P = 20.0 / 1.013 * one_atm
        local phi = l_phi[i]
        l_idt[j, 3] = f_idt(T0, P, phi, p.*0; dT=800) * 1000.0
        l_idt[j, 4] = f_idt(T0, P, phi, p; dT=800) * 1000.0
        l_idt[j, 5] = f_idt(T0, P, phi, p_ga; dT=800) * 1000.0
    end
    # plt = Plots.scatter(l_idt[:, 1], l_idt[:, 2], yscale=:log10, label="exp")
    # Plots.plot!(plt, l_idt[:, 1], l_idt[:, 3], lw=1, yscale=:log10, label="base")
    # Plots.plot!(plt, l_idt[:, 1], l_idt[:, 4], lw=3, yscale=:log10, label="train")
    # # plot!(plt, [minimum(l_idt), maximum(l_idt)], [minimum(l_idt), maximum(l_idt)], label="y=x")
    # xlabel!(plt, "1000 / T [1/K]")
    # ylabel!(plt, "IDT [ms]")
    # xlims!(plt, (0.7, 1.7))
    # title!(plt, "20 bar phi=$(l_phi[i])")
    # plot!(plt, legend=:topleft, framestyle = :box, foreground_color_legend = nothing)
    # push!(l_plt, plt)

    writedlm(string(fig_path, "phi=$(l_phi[i])"), l_idt)
end

plt = plot(l_plt..., layout=(3, 1), size=(400, 1200))
png(plt, string(fig_path, "/regression_T"))
gr()

max_stars_repo_name: DENG-MIT/DeepReduction

max_issues_repo_name: DENG-MIT/DeepReduction

max_forks_repo_name: DENG-MIT/DeepReduction


text:
using Plots
using Statistics

function plot_mean_endmembers(endmember_library; output_name = "")
    plots = []
    for (_u, u) in enumerate(endmember_library.class_valid_keys)
        mean_spectra = mean(endmember_library.spectra[endmember_library.classes .== u,:],dims=1)[:]
        plot!(endmember_library.wavelengths, mean_spectra, label=u, xlim=[300,3200])
    end
    xlabel!("Wavelength [nm]")
    ylabel!("Reflectance")
    xticks!([500, 1000, 1500, 2000, 2500, 3000])
    plot!(dpi=200)
    if output_name != ""
        savefig(output_name)
    end
    plot!()
end

function plot_endmembers(endmember_library::SpectralLibrary; output_name::String = "")

    for (_u, u) in enumerate(endmember_library.class_valid_keys)
        if _u == 1
            plot(endmember_library.wavelengths, endmember_library.spectra[endmember_library.classes .== u,:]', lab="", xlim=[300,3200], color=palette(:tab10)[Mod(_u)],dpi=200)
        else
            plot!(endmember_library.wavelengths, endmember_library.spectra[endmember_library.classes .== u,:]', lab="",xlim=[300,3200], color=palette(:tab10)[Mod(_u)])
        end
    end
    xlabel!("Wavelenth [nm]")
    ylabel!("Reflectance")
    xticks!([500, 1000, 1500, 2000, 2500, 3000])
    for (_u, u) in enumerate(endmember_library.class_valid_keys)
        plot!([1:2],[0,0.3], color=palette(:tab10)[Mod(_u)], labels=u)
    end
    if output_name != ""
        savefig(output_name)
    end
    plot!()
end

function plot_endmembers_individually(endmember_library::SpectralLibrary; output_name::String = "")
    plots = []
    spectra = endmember_library.spectra
    classes = endmember_library.classes
    for (_u, u) in enumerate(endmember_library.class_valid_keys)
        sp = spectra[classes .== u,:]
        sp[broadcast(isnan,sp)] .= 0
        toplot = spectra[classes .== u,:]

        if size(toplot)[2] == sum(endmember_library.good_bands)
            push!(plots, plot(endmember_library.wavelengths[endmember_library.good_bands], toplot', title=u, xlabel="Wavelength [nm]", ylabel="Reflectance"))
        else
            push!(plots, plot(endmember_library.wavelengths, toplot', title=u, xlabel="Wavelength [nm]", ylabel="Reflectance"))
        end
        xticks!([500, 1000, 1500, 2000, 2500])
    end
    plot(plots...,size=(1000,600),dpi=200)
    plot!(legend=:outertopright)
    if output_name != ""
        savefig(output_name)
    end
    plot!()
end

max_stars_repo_name: fotxoa-geo/SpectralUnmixing

max_issues_repo_name: fotxoa-geo/SpectralUnmixing

max_forks_repo_name: fotxoa-geo/SpectralUnmixing


text:
# module Bukdu

export render

struct UnknownModuleError <: Exception
    msg::String
end

"""
    render(::Type{Text}, data)::Render
"""
function render(::Type{Text}, data)::Render
    Render("text/plain; charset=utf-8", unsafe_wrap(Vector{UInt8}, string(data)))
end

"""
    render(::Type{HTML}, data)::Render
"""
function render(::Type{HTML}, data)::Render
    Render("text/html; charset=utf-8", unsafe_wrap(Vector{UInt8}, string(data)))
end

"""
    render(::Type{JSON}, data)::Render
"""
function render(::Type{JSON}, data)::Render
    Render("application/json; charset=utf-8", unsafe_wrap(Vector{UInt8}, JSON2.write(data)))
end

"""
    render(::Type{JavaScript}, data)::Render
"""
function render(::Type{JavaScript}, data)::Render
    Render("application/javascript; charset=utf-8", unsafe_wrap(Vector{UInt8}, string(data)))
end

"""
    render(::Type{EventStream})::EventStream
"""
function render(::Type{EventStream})::EventStream
    EventStream("text/event-stream; charset=utf-8", [])
end

# application/wasm

function render(m::Module, data)::AbstractRender # throw UnknownModuleError
    if nameof(m) == :JSON
        render(JSON, data)
    elseif nameof(m) == :HTML5
        render(HTML, data)
    else
        throw(UnknownModuleError(string(m)))
    end
end

# module Bukdu

max_stars_repo_name: wookay/TestBukdu

max_issues_repo_name: wookay/TestBukdu

max_forks_repo_name: wookay/TestBukdu


text:
module Agent

using ...dlgo

export AbstractAgent
export RandomAgent
export is_point_an_eye, select_move

export MCTSTree, MCTSNode, MCTSAgent, NodeInfo
export is_leaf, expand_leaf!, update_tree!, search_probabilities, thompson_sampling,
        num_visits, num_wins, info

include("utils.jl")
include("randomplay.jl")
include("mctstree.jl")

end


text:
module EMLtank_basic
	using DanaTypes
	using DotPlusInheritance
	using Reexport
	@reexport using ...streams.EMLstreams
	import EMLtypes.length
	@reexport using ...reactors.EMLvol_tank
	import EMLtypes.length
	include("tank_basic/tank_basic.jl")
	include("tank_basic/tank_vap.jl")
	include("tank_basic/tank_liq.jl")
	include("tank_basic/tank_liqvap.jl")
end


text:
s = BenchmarkGroup()
suite["traversals"] = s

g = DiGraph(1000,5000, seed=1)
s["minimum cut"] = @benchmarkable minimum_cut($g)


text:
using Econometrics, DelimitedFiles, StatsPlots, MCMCChains

chain = readdlm("chain")
chn = Chains(chain,["μ", "κ","α","σ","ρ","λ0","λ1","τ"])
display(chn)
plot(chn)
savefig("chain.png")

p = npdensity(chain[:,1])
plot!(p, legend=false)
savefig(p, "mu.png")

p = npdensity(chain[:,2])
plot!(legend=false)
savefig(p, "kappa.png")

p = npdensity(chain[:,3])
plot!(p, legend=false)
savefig(p, "alpha.png")

p = npdensity(chain[:,4])
plot!(p, legend=false)
savefig(p, "sigma.png")

p = npdensity(chain[:,5])
plot!(p, legend=false)
savefig(p, "rho.png")

p = npdensity(chain[:,6])
plot!(p, legend=false)
savefig(p, "lambda0.png")

p = npdensity(chain[:,7])
plot!(p, legend=false)
savefig(p, "lambda1.png")

p = npdensity(chain[:,8])
plot!(p, legend=false)
savefig(p, "tau.png")

max_stars_repo_name: mcreel/ILSNM

max_issues_repo_name: mcreel/ILSNM

max_forks_repo_name: mcreel/ILSNM


text:
using ECOSProblems
using LinearAlgebra
using Test

function dual_inf(x, z, bl, bu)
    dInf = 0.0; jInf = 0
    for j = 1:n
	      if bl[j] < bu[j]
	          dj = z[j]
	          if dj > 0
		            dj =  dj * min(x[j] - bl[j], 1.0)
	          elseif dj < 0
		            dj = -dj * min(bu[j] - x[j], 1.0)
            end
	          if dInf < dj
		            dInf = dj
		            jInf = j
            end
	      end
    end
    return dInf, jInf
end

@testset "ECOSProblems.jl" begin
    atol = 1e-6
    n = 3; m = 100
    A = randn(m, n)
    x0 = rand(n)
    b = A*x0
    u = rand(n)
    x, r, rNorm = lsbox(A, b, u, verbose=false)
    @test isapprox(norm(r), rNorm, atol=atol)
    @test isapprox(norm(b-A*x), norm(r), atol=atol)
    dInf, jInf = dual_inf(x, -(A'r), zeros(n), u)
    @test dInf ≤ 1e-3
    # println("dual Inf = $dInf, var = $jInf, $(g[jInf])   $(u[jInf]-x[jInf]) ")
end


text:
using BrainFlow

# specify logging library to use
BrainFlow.enable_dev_logger(BrainFlow.BOARD_CONTROLLER)

params = BrainFlowInputParams()
board_shim = BrainFlow.BoardShim(BrainFlow.SYNTHETIC_BOARD, params)

BrainFlow.prepare_session(board_shim)
BrainFlow.start_stream(board_shim)
sleep(5)
BrainFlow.stop_stream(board_shim)
data = BrainFlow.get_current_board_data(32, board_shim)
BrainFlow.release_session(board_shim)

BrainFlow.write_file(data, "test.csv", "w")
restored_data = BrainFlow.read_file("test.csv")

println("Original Data")
println(data)
println("Restored Data")
println(restored_data)


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using Test

Test.@testset "nested_sampling" begin
    # Python ultranest package doesn't seem to be available via Conda on
    # 32-bit systems:
    if Int == Int64
        include("test_ultranest.jl")
    end
end


text:
@testset "Geometric operations" begin
  @testset "unique" begin
    X = [i*j for i in 1:2, j in 1:1_000_000]
    z = rand(1_000_000)
    d = georef((z=[z;z],), [X X])
    u = uniquecoords(d)
    p = [centroid(u, i) for i in 1:nelements(u)]
    U = reduce(hcat, coordinates.(p))
    @test nelements(u) == 1_000_000
    @test Set(eachcol(U)) == Set(eachcol(X))

    X = rand(3,100)
    z = rand(100)
    n = [string(i) for i in 1:100]
    Xd = hcat(X, X[:,1:10])
    zd = vcat(z, z[1:10])
    nd = vcat(n, n[1:10])
    sdata = georef((z=zd, n=nd), PointSet(Xd))
    ndata = uniquecoords(sdata)
    @test nelements(ndata) == 100
  end

  @testset "groupby" begin
    setify(lists) = Set(Set.(lists))

    d = georef((z=[1,2,3],x=[4,5,6]), rand(2,3))
    g = groupby(d, :z)
    @test all(nelements.(g) .== 1)
    rows = [[1 4], [2 5], [3 6]]
    for i in 1:3
      @test Tables.matrix(values(g[i])) ∈ rows
    end

    z = vec([1 1 1; 2 2 2; 3 3 3])
    sdata = georef((z=z,), CartesianGrid(3,3))
    p = groupby(sdata, :z)
    @test setify(subsets(p)) == setify([[1,4,7],[2,5,8],[3,6,9]])

    # groupby with missing values
    z = vec([missing 1 1; 2 missing 2; 3 3 missing])
    sdata = georef((z=z,), CartesianGrid(3,3))
    p = groupby(sdata, :z)
    @test setify(subsets(p)) == setify([[4,7],[2,8],[3,6],[1,5,9]])
  end

  @testset "filter" begin
    𝒟 = georef((a=[1,2,3], b=[1,1,missing]))
    𝒫 = filter(s -> !ismissing(s.b), 𝒟)
    @test 𝒫[:a] == [1,2]
    @test 𝒫[:b] == [1,1]

    𝒟 = georef((a=[1,2,3],b=[3,2,1]))
    𝒫ₐ = filter(s -> s.a > 1, 𝒟)
    𝒫ᵦ = filter(s -> s.b > 1, 𝒟)
    𝒫ₐᵦ = filter(s -> s.a > 1 && s.b > 1, 𝒟)
    @test nelements(𝒫ₐ) == 2
    @test nelements(𝒫ᵦ) == 2
    @test nelements(𝒫ₐᵦ) == 1
    @test 𝒫ₐ[:a] == [2,3]
    @test 𝒫ₐ[:b] == [2,1]
    @test 𝒫ᵦ[:a] == [1,2]
    @test 𝒫ᵦ[:b] == [3,2]
    @test 𝒫ₐᵦ[:a] == [2]
    @test 𝒫ₐᵦ[:b] == [2]
  end
end


text:
module MinionsREPL

function __init__()
    if ccall(:jl_generating_output, Cint, ()) == 0
        include(joinpath(@__DIR__, "repl.jl"))
    end
end

end


text:
for jy in ("j","y"), nu in (0,1)
    jynu = Expr(:quote, symbol(string(jy,nu)))
    jynuf = Expr(:quote, symbol(string(jy,nu,"f")))
    bjynu = symbol(string("bessel",jy,nu))
    if jy == "y"
        @eval begin
            $bjynu(x::Float64) = nan_dom_err(ccall(($jynu,libm),  Float64, (Float64,), x), x)
            $bjynu(x::Float32) = nan_dom_err(ccall(($jynuf,libm), Float32, (Float32,), x), x)
        end
    else
        @eval begin
            $bjynu(x::Float64) = ccall(($jynu,libm),  Float64, (Float64,), x)
            $bjynu(x::Float32) = ccall(($jynuf,libm), Float32, (Float32,), x)
        end
    end
    @eval begin
        $bjynu(x::Real) = $bjynu(float(x))
        $bjynu(x::Complex) = $(symbol(string("bessel",jy)))($nu,x)
        @vectorize_1arg Number $bjynu
    end
end

type AmosException <: Exception
    info::Int32
end

let
    const ai::Array{Float64,1} = Array(Float64,2)
    const ae::Array{Int32,1} = Array(Int32,2)
    global _airy, _biry
    function _airy(z::Complex128, id::Int32, kode::Int32)
        ccall((:zairy_,openspecfun), Void,
              (Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32},
               Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}),
              &real(z), &imag(z),
              &id, &kode,
              pointer(ai,1), pointer(ai,2),
              pointer(ae,1), pointer(ae,2))
        if ae[2] == 0 || ae[2] == 3
            return complex(ai[1],ai[2])
        else
            throw(AmosException(ae[2]))
        end
    end
    function _biry(z::Complex128, id::Int32, kode::Int32)
        ccall((:zbiry_,openspecfun), Void,
              (Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32},
               Ptr{Float64}, Ptr{Float64}, Ptr{Int32}),
              &real(z), &imag(z),
              &id, &kode,
              pointer(ai,1), pointer(ai,2),
              pointer(ae,1))
        if ae[1] == 0 || ae[1] == 3  # ignore underflow
            return complex(ai[1],ai[2])
        else
            throw(AmosException(ae[2]))
        end
    end
end

function airy(k::Int, z::Complex128)
    id = int32(k==1 || k==3)
    if k == 0 || k == 1
        return _airy(z, id, int32(1))
    elseif k == 2 || k == 3
        return _biry(z, id, int32(1))
    else
        error("invalid argument")
    end
end

airy(z) = airy(0,z)
@vectorize_1arg Number airy
airyprime(z) = airy(1,z)
@vectorize_1arg Number airyprime
airyai(z) = airy(0,z)
@vectorize_1arg Number airyai
airyaiprime(z) = airy(1,z)
@vectorize_1arg Number airyaiprime
airybi(z) = airy(2,z)
@vectorize_1arg Number airybi
airybiprime(z) = airy(3,z)
@vectorize_1arg Number airybiprime

airy(k::Number, x::FloatingPoint) = oftype(x, real(airy(k, complex(x))))
airy(k::Number, x::Real) = airy(k, float(x))
airy(k::Number, z::Complex64) = complex64(airy(k, complex128(z)))
airy(k::Number, z::Complex) = airy(convert(Int,k), complex128(z))
@vectorize_2arg Number airy

function airyx(k::Int, z::Complex128)
    id = int32(k==1 || k==3)
    if k == 0 || k == 1
        return _airy(z, id, int32(2))
    elseif k == 2 || k == 3
        return _biry(z, id, int32(2))
    else
        error("invalid argument")
    end
end

airyx(z) = airyx(0,z)
@vectorize_1arg Number airyx

airyx(k::Number, x::FloatingPoint) = oftype(x, real(airyx(k, complex(x))))
airyx(k::Number, x::Real) = airyx(k, float(x))
airyx(k::Number, z::Complex64) = complex64(airyx(k, complex128(z)))
airyx(k::Number, z::Complex) = airyx(convert(Int,k), complex128(z))
@vectorize_2arg Number airyx

const cy = Array(Float64,2)
const ae = Array(Int32,2)
const wrk = Array(Float64,2)

function _besselh(nu::Float64, k::Int32, z::Complex128, kode::Int32)
    ccall((:zbesh_,openspecfun), Void,
          (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32},
           Ptr{Int32}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}),
          &real(z), &imag(z), &nu, &kode, &k, &1,
          pointer(cy,1), pointer(cy,2),
          pointer(ae,1), pointer(ae,2))
    if ae[2] == 0 || ae[2] == 3
        return complex(cy[1],cy[2])
    else
        throw(AmosException(ae[2]))
    end
end

function _besseli(nu::Float64, z::Complex128, kode::Int32)
    ccall((:zbesi_,openspecfun), Void,
          (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32},
           Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}),
          &real(z), &imag(z), &nu, &kode, &1,
          pointer(cy,1), pointer(cy,2),
          pointer(ae,1), pointer(ae,2))
    if ae[2] == 0 || ae[2] == 3
        return complex(cy[1],cy[2])
    else
        throw(AmosException(ae[2]))
    end
end

function _besselj(nu::Float64, z::Complex128, kode::Int32)
    ccall((:zbesj_,openspecfun), Void,
          (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32},
           Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}),
          &real(z), &imag(z), &nu, &kode, &1,
          pointer(cy,1), pointer(cy,2),
          pointer(ae,1), pointer(ae,2))
    if ae[2] == 0 || ae[2] == 3
        return complex(cy[1],cy[2])
    else
        throw(AmosException(ae[2]))
    end
end

function _besselk(nu::Float64, z::Complex128, kode::Int32)
    ccall((:zbesk_,openspecfun), Void,
          (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32},
           Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}),
          &real(z), &imag(z), &nu, &kode, &1,
          pointer(cy,1), pointer(cy,2),
          pointer(ae,1), pointer(ae,2))
    if ae[2] == 0 || ae[2] == 3
        return complex(cy[1],cy[2])
    else
        throw(AmosException(ae[2]))
    end
end

function _bessely(nu::Float64, z::Complex128, kode::Int32)
    ccall((:zbesy_,openspecfun), Void,
          (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32},
           Ptr{Int32}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32},
           Ptr{Float64}, Ptr{Float64}, Ptr{Int32}),
          &real(z), &imag(z), &nu, &kode, &1,
          pointer(cy,1), pointer(cy,2),
          pointer(ae,1), pointer(wrk,1),
          pointer(wrk,2), pointer(ae,2))
    if ae[2] == 0 || ae[2] == 3
        return complex(cy[1],cy[2])
    else
        throw(AmosException(ae[2]))
    end
end

function besselh(nu::Float64, k::Integer, z::Complex128)
    if nu < 0
        s = (k == 1) ? 1 : -1
        return _besselh(-nu,int32(k),z,int32(1)) * complex(cospi(nu),-s*sinpi(nu))
    end
    return _besselh(nu,int32(k),z,int32(1))
end

function besselhx(nu::Float64, k::Integer, z::Complex128)
    if nu < 0
        s = (k == 1) ? 1 : -1
        return _besselh(-nu,int32(k),z,int32(2)) * complex(cospi(nu),-s*sinpi(nu))
    end
    return _besselh(nu,int32(k),z,int32(2))
end

function besseli(nu::Float64, z::Complex128)
    if nu < 0
        return _besseli(-nu,z,int32(1)) - 2_besselk(-nu,z,int32(1))*sinpi(nu)/pi
    else
        return _besseli(nu,z,int32(1))
    end
end

function besselix(nu::Float64, z::Complex128)
    if nu < 0
        return _besseli(-nu,z,int32(2)) - 2_besselk(-nu,z,int32(2))*exp(-abs(real(z))-z)*sinpi(nu)/pi
    else
        return _besseli(nu,z,int32(2))
    end
end

function besselj(nu::Float64, z::Complex128)
    if nu < 0
        return _besselj(-nu,z,int32(1))*cospi(nu) + _bessely(-nu,z,int32(1))*sinpi(nu)
    else
        return _besselj(nu,z,int32(1))
    end
end

besselj(nu::Integer, x::FloatingPoint) = typemin(Int32) <= nu <= typemax(Int32) ?
    oftype(x, ccall((:jn, libm), Float64, (Cint, Float64), nu, x)) :
    besselj(float64(nu), x)

besselj(nu::Integer, x::Float32) = typemin(Int32) <= nu <= typemax(Int32) ?
    ccall((:jnf, libm), Float32, (Cint, Float32), nu, x) :
    besselj(float64(nu), x)

function besseljx(nu::Float64, z::Complex128)
    if nu < 0
        return _besselj(-nu,z,int32(2))*cospi(nu) + _bessely(-nu,z,int32(2))*sinpi(nu)
    else
        return _besselj(nu,z,int32(2))
    end
end

besselk(nu::Float64, z::Complex128) = _besselk(abs(nu), z, int32(1))

besselkx(nu::Float64, z::Complex128) = _besselk(abs(nu), z, int32(2))

function bessely(nu::Float64, z::Complex128)
    if nu < 0
        return _bessely(-nu,z,int32(1))*cospi(nu) - _besselj(-nu,z,int32(1))*sinpi(nu)
    else
        return _bessely(nu,z,int32(1))
    end
end

function besselyx(nu::Float64, z::Complex128)
    if nu < 0
        return _bessely(-nu,z,int32(2))*cospi(nu) - _besselj(-nu,z,int32(2))*sinpi(nu)
    else
        return _bessely(nu,z,int32(2))
    end
end

besselh(nu, z) = besselh(nu, 1, z)
besselh(nu::Real, k::Integer, z::Complex64) = complex64(besselh(float64(nu), k, complex128(z)))
besselh(nu::Real, k::Integer, z::Complex) = besselh(float64(nu), k, complex128(z))
besselh(nu::Real, k::Integer, x::Real) = besselh(float64(nu), k, complex128(x))
@vectorize_2arg Number besselh

hankelh1(nu, z) = besselh(nu, 1, z)
@vectorize_2arg Number hankelh1

hankelh2(nu, z) = besselh(nu, 2, z)
@vectorize_2arg Number hankelh2

besselhx(nu::Real, k::Integer, z::Complex64) = complex64(besselhx(float64(nu), k, complex128(z)))
besselhx(nu::Real, k::Integer, z::Complex) = besselhx(float64(nu), k, complex128(z))
besselhx(nu::Real, k::Integer, x::Real) = besselhx(float64(nu), k, complex128(x))

hankelh1x(nu, z) = besselhx(nu, 1, z)
@vectorize_2arg Number hankelh1x

hankelh2x(nu, z) = besselhx(nu, 2, z)
@vectorize_2arg Number hankelh2x

function besseli(nu::Real, x::FloatingPoint)
    if x < 0 && !isinteger(nu)
        throw(DomainError())
    end
    oftype(x, real(besseli(float64(nu), complex128(x))))
end

function besselix(nu::Real, x::FloatingPoint)
    if x < 0 && !isinteger(nu)
        throw(DomainError())
    end
    oftype(x, real(besselix(float64(nu), complex128(x))))
end

function besselj(nu::FloatingPoint, x::FloatingPoint)
    if isinteger(nu)
        if typemin(Int32) <= nu <= typemax(Int32)
            return besselj(int(nu), x)
        end
    elseif x < 0
        throw(DomainError())
    end
    oftype(x, real(besselj(float64(nu), complex128(x))))
end

function besseljx(nu::Real, x::FloatingPoint)
    if x < 0 && !isinteger(nu)
        throw(DomainError())
    end
    oftype(x, real(besseljx(float64(nu), complex128(x))))
end

function besselk(nu::Real, x::FloatingPoint)
    if x < 0
        throw(DomainError())
    end
    if x == 0
        return oftype(x, Inf)
    end
    oftype(x, real(besselk(float64(nu), complex128(x))))
end

function besselkx(nu::Real, x::FloatingPoint)
    if x < 0
        throw(DomainError())
    end
    if x == 0
        return oftype(x, Inf)
    end
    oftype(x, real(besselkx(float64(nu), complex128(x))))
end

function bessely(nu::Real, x::FloatingPoint)
    if x < 0
        throw(DomainError())
    end
    if isinteger(nu) && typemin(Int32) <= nu <= typemax(Int32)
        return bessely(int(nu), x)
    end
    oftype(x, real(bessely(float64(nu), complex128(x))))
end
function bessely(nu::Integer, x::FloatingPoint)
    if x < 0
        throw(DomainError())
    end
    return oftype(x, ccall((:yn, libm), Float64, (Cint, Float64), nu, x))
end
function bessely(nu::Integer, x::Float32)
    if x < 0
        throw(DomainError())
    end
    return ccall((:ynf, libm), Float32, (Cint, Float32), nu, x)
end

function besselyx(nu::Real, x::FloatingPoint)
    if x < 0
        throw(DomainError())
    end
    oftype(x, real(besselyx(float64(nu), complex128(x))))
end

for f in ("i", "ix", "j", "jx", "k", "kx", "y", "yx")
    bfn = symbol(string("bessel", f))
    @eval begin
        $bfn(nu::Real, z::Complex64) = complex64($bfn(float64(nu), complex128(z)))
        $bfn(nu::Real, z::Complex) = $bfn(float64(nu), complex128(z))
        $bfn(nu::Real, x::Integer) = $bfn(nu, float64(x))
        @vectorize_2arg Number $bfn
    end
end

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia