text:
module GaussianDistributions

# Gaussian
using LinearAlgebra, Random, Statistics
using Distributions
using LinearAlgebra: norm_sqr

import Random: rand, GLOBAL_RNG
import Statistics: mean, cov, var, std
import Distributions: pdf, logpdf, sqmahal, cdf, quantile
import LinearAlgebra: cholesky
import Base: size, iterate, length

export PSD, Gaussian

"""
    PSD{T}

Simple wrapper for the lower triangular Cholesky root of a positive (semi-)definite element `σ`.
"""
struct PSD{T}
    σ::T
    PSD(σ::T) where {T} = istril(σ) ? new{T}(σ) : throw(ArgumentError("Argument not lower triangular"))
end
cholesky(P::PSD) = (U=P.σ',)   # using a named tuple for now. FIXME? Also: TESTME

"""
Sum of the log of the diagonal elements. Second argument `d` is used
to handle `UniformScaling` and other linear operators whose dimensions
are determined by the dimension of the argument they work on.
"""
sumlogdiag(Σ::Float64, _) = log(Σ)
sumlogdiag(Σ, d) = sum(log.(diag(Σ)))
sumlogdiag(J::UniformScaling, d) = log(J.λ)*d

_logdet(Σ::PSD, d) = 2*sumlogdiag(Σ.σ, d)

_logdet(Σ, d) = logdet(Σ)
_logdet(J::UniformScaling, d) = log(J.λ) * d

_symmetric(Σ) = Symmetric(Σ)
_symmetric(J::UniformScaling) = J

"""
    Gaussian(μ, Σ) -> P

Gaussian distribution with mean `μ` and covariance `Σ`. Defines `rand(P)` and `(log-)pdf(P, x)`.
Designed to work with `Number`s, `UniformScaling`s, `StaticArrays` and `PSD`-matrices.

Implementation details: On `Σ` the functions `logdet`, `whiten` and `unwhiten`
(or `chol` as fallback for the latter two) are called.
"""
struct Gaussian{T,S}
    μ::T
    Σ::S
    Gaussian{T,S}(μ, Σ) where {T,S} = new(μ, Σ)
    Gaussian(μ::T, Σ::S) where {T,S} = new{T,S}(μ, Σ)
end

iterate(P::Gaussian) = P.μ, true
iterate(P::Gaussian, s) = s ? (P.Σ, false) : nothing
length(P::Gaussian) = 2

pair(u) = u[1], u[2]
pair(p::Gaussian) = p.μ, p.Σ

Base.:(==)(g1::Gaussian, g2::Gaussian) = g1.μ == g2.μ && g1.Σ == g2.Σ
Base.isapprox(g1::Gaussian, g2::Gaussian; kwargs...) =
    isapprox(g1.μ, g2.μ; kwargs...) && isapprox(g1.Σ, g2.Σ; kwargs...)
Gaussian() = Gaussian(0.0, 1.0)
mean(P::Gaussian) = P.μ
cov(P::Gaussian) = P.Σ
var(P::Gaussian{<:Number}) = P.Σ
std(P::Gaussian{<:Number}) = sqrt(var(P))
Base.convert(::Type{Gaussian{T, S}}, g::Gaussian) where {T, S} =
    Gaussian(convert(T, g.μ), convert(S, g.Σ))

dim(P::Gaussian) = length(P.μ)
whiten(Σ::PSD, z) = Σ.σ\z
whiten(Σ, z) = cholesky(Σ).U'\z
whiten(Σ::Number, z) = sqrt(Σ)\z
whiten(Σ::UniformScaling, z) = sqrt(Σ.λ)\z

unwhiten(Σ::PSD, z) = Σ.σ*z
unwhiten(Σ, z) = cholesky(Σ).U'*z
unwhiten(Σ::Number, z) = sqrt(Σ)*z
unwhiten(Σ::UniformScaling, z) = sqrt(Σ.λ)*z

sqmahal(P::Gaussian, x) = norm_sqr(whiten(P.Σ, x - P.μ))

rand(P::Gaussian) = rand(GLOBAL_RNG, P)
rand(RNG::AbstractRNG, P::Gaussian) = P.μ + unwhiten(P.Σ, randn(RNG, typeof(P.μ)))
rand(RNG::AbstractRNG, P::Gaussian{Vector{T}}) where T =
    P.μ + unwhiten(P.Σ, randn(RNG, T, length(P.μ)))
rand(RNG::AbstractRNG, P::Gaussian{<:Number}) = P.μ + sqrt(P.Σ)*randn(RNG, typeof(one(P.μ)))

logpdf(P::Gaussian, x) = -(sqmahal(P,x) + _logdet(P.Σ, dim(P)) + dim(P)*log(2pi))/2
pdf(P::Gaussian, x) = exp(logpdf(P::Gaussian, x))
cdf(P::Gaussian{Number}, x) = Distributions.normcdf(P.μ, sqrt(P.Σ), x)

Base.:+(g::Gaussian, vec) = Gaussian(g.μ + vec, g.Σ)
Base.:+(vec, g::Gaussian) = g + vec
Base.:-(g::Gaussian, vec) = g + (-vec)
Base.:*(M, g::Gaussian) = Gaussian(M * g.μ, M * g.Σ * M')
⊕(g1::Gaussian, g2::Gaussian) = Gaussian(g1.μ + g2.μ, g1.Σ + g2.Σ)
⊕(vec, g::Gaussian) = vec + g
⊕(g::Gaussian, vec) = g + vec
const independent_sum = ⊕

function rand_scalar(RNG::AbstractRNG, P::Gaussian{T}, dims) where {T}
    X = zeros(T, dims)
    for i in 1:length(X)
        X[i] = rand(RNG, P)
    end
    X
end

function rand_vector(RNG::AbstractRNG, P::Gaussian{Vector{T}}, dims::Union{Integer, NTuple}) where {T}
    X = zeros(T, dim(P), dims...)
    for i in 1:prod(dims)
        X[:, i] = rand(RNG, P)
    end
    X
end
rand(RNG::AbstractRNG, P::Gaussian, dim::Integer) = rand_scalar(RNG, P, dim)
rand(RNG::AbstractRNG, P::Gaussian, dims::Tuple{Vararg{Int64,N}} where N) = rand_scalar(RNG, P, dims)

rand(RNG::AbstractRNG, P::Gaussian{Vector{T}}, dim::Integer) where {T} = rand_vector(RNG, P, dim)
rand(RNG::AbstractRNG, P::Gaussian{Vector{T}}, dims::Tuple{Vararg{Int64,N}} where N) where {T} = rand_vector(RNG, P, dims)
rand(P::Gaussian, dims::Tuple{Vararg{Int64,N}} where N) = rand(GLOBAL_RNG, P, dims)
rand(P::Gaussian, dim::Integer) = rand(GLOBAL_RNG, P, dim)

"""
    logpdfnormal(x, Σ)

Logarithm of the probability density function of centered Gaussian with covariance Σ
"""
function logpdfnormal(x, Σ)

    S = cholesky(_symmetric(Σ)).U'

    d = length(x)
     -((norm(S\x))^2 + 2sumlogdiag(S,d) + d*log(2pi))/2
end
function logpdfnormal(x::Float64, Σ)
     -(x^2/Σ + log(Σ) + log(2pi))/2
end

"""
    conditional(P::Gaussian, A, B, xB)

Conditional distribution of `X[i for i in A]` given
`X[i for i in B] == xB` if ``X ~ P``.
"""
function conditional(P::Gaussian, A, B, xB)
    Z = P.Σ[A,B]*inv(P.Σ[B,B])
    Gaussian(P.μ[A] + Z*(xB - P.μ[B]), P.Σ[A,A] - Z*P.Σ[B,A])
end

"""
    correct(u, v, H) = u, yres, S

Joseph form correction step of a Kalman filter with `u = Gaussian(x, P)` state
and `v = Gaussian(y, R)` the observation with uncertainty `R`.
`H` is the observation operator. Returns corrected/conditional
distribution u, the residual and the innovation covariance.

See https://en.wikipedia.org/wiki/Kalman_filter#Update.
"""
function correct(u::T, v, H) where {T}
    x, Ppred = pair(u)
    y, R = pair(v)
    yres = y - H*x # innovation residual

    S = (H*Ppred*H' + R) # innovation covariance

    K = Ppred*H'/S # Kalman gain
    x = x + K*yres
    P = (I - K*H)*Ppred*(I - K*H)' + K*R*K' #  Joseph form

    T(x, P), yres, S
end

include("bivariate.jl")

end # module


text:
# Generate figures on the Toronto street network

# Load necessary packages
using Plots, SparseArrays, JLD2, LinearAlgebra, MultiscaleGraphSignalTransforms
include("auxilaries.jl")

# Set up the resolution and display size
gr(dpi=200, size=(800,600))

# Load the Toronto street network (vehicular volume counts as its graph signal)
JLD2.@load "Toronto_new.jld2"
tmp1 = toronto["G"]
G=GraphSig(tmp1["W"],xy=tmp1["xy"],f=tmp1["f"],name =tmp1["name"],
    plotspecs = tmp1["plotspecs"])

# Generate Fig. 10a
# GraphSig_Plot(G, linewidth = 1., markersize = 4.,
#    markercolor = :viridis, markerstrokealpha = 0., notitle = true)
gplot(G.W, G.xy; width=1);
plt = scatter_gplot!(G.xy; marker = G.f, plotOrder = :s2l,
    ms = G.f*8/maximum(G.f), mswidth=1)
display(plt)
savefig("toronto_vv.pdf")
savefig("toronto_vv.png")

# Assign correct edge weights via 1/Euclidean distance
G = Adj2InvEuc(G)
# Perform the hierarchical graph partitioning using the Fiedler vectors
GP = partition_tree_fiedler(G) # Lrw by default
# Compute the expansion coefficients of the full GHWT c2f dictionary
dmatrix = ghwt_analysis!(G, GP=GP)

# Extract the coefficients corresponding to the Haar basis
BS_haar = bs_haar(GP)
dvec_haar = dmatrix2dvec(dmatrix, GP, BS_haar)

# Extract the coefficients corresponding to the Walsh basis
BS_walsh = bs_walsh(GP)
dvec_walsh = dmatrix2dvec(dmatrix, GP, BS_walsh)

# Compute the c2f GHWT best basis and the coefficients
dvec_c2f, BS_c2f = ghwt_c2f_bestbasis(dmatrix, GP)

# Compute the f2c GHWT best basis and the coefficients
dvec_f2c, BS_f2c = ghwt_f2c_bestbasis(dmatrix, GP)

# Compute the eGHWT best basis and the coefficients
dvec_eghwt, BS_eghwt = eghwt_bestbasis(dmatrix, GP)

################################################################################
###################### Generate approximaation errors ##########################
################################################################################
# Generate Fig. 10b up to 50% of the coefficients retained.
DVEC = [ dvec_haar[:], dvec_walsh[:], dvec_c2f[:], dvec_f2c[:], dvec_eghwt[:] ]
T = [ "Graph Haar","Graph Walsh","GHWT_c2f", "GHWT_f2c", "eGHWT" ]
L = [ (:dashdot, :orange), (:dashdot, :blue), (:solid, :red), (:solid, :green),
      (:solid, :black) ]
approx_error2(DVEC, T, L, 0.5)
savefig("toronto_vv_approx_error.pdf")
savefig("toronto_vv_approx_error.png")

################################################################################
######################### Synthesis by top vectors #############################
################################################################################
color_limit_residual = (0., 0.15)

function top_vectors_residual(p::Int64, dvec::Array{Float64,2},
                                BS::BasisSpec, GP::GraphPart, G::GraphSig)
    sorted_dvec = sort(abs.(dvec[:]), rev = true)
    dvecT = copy(dvec)
    dvecT[abs.(dvec) .< sorted_dvec[p]].= 0
    (recon, GS)= ghwt_synthesis(dvecT, GP, BS, G)
    GS.name = "Square of the residual"
    #GS.f =  (G.f - GS.f).^2
    #print((maximum(GS.f), minimum(GS.f)))
    #GraphSig_Plot(GS)
    GS.f = (G.f - recon).^2 ./(G.f.^2)
    GraphSig_Plot(GS, linewidth = 1., markersize = 4., markercolor = :viridis,
        markerstrokealpha =0., clim = color_limit_residual)
    display(current())
end

################################################################################
############## Visuallization of Top basis vectors #############################
################################################################################

### function to plot the basis vectors of user specified range
# Recommended number of basis vectors are squared numbers, e.g., 9, 16, 25.
function top_vectors_plot(dvec::Array{Float64, 2}, BS::BasisSpec, GP::GraphPart,
                             G::GraphSig; istart = 2, iend = 17, ms = 2)
    sorted_ind = sortperm(dvec[:].^2; rev = true)
    ibs = iend-istart+1
    n1 = Int64(sqrt(ibs))
    plot(layout = Plots.grid(n1, n1))
    idsp = 1
    for ib in istart:iend
        dvecT = fill(0., size(dvec))
        dvecT[sorted_ind[ib]] = 1
        bv = ghwt_synthesis(dvecT, GP, BS)
        gplot!(G.W, G.xy; width=0.25, subplot=idsp, color=:gray)
        scatter_gplot!(G.xy; marker = bv, plotOrder = :propabs, ms, subplot=idsp)
        plt = plot!(cbar = false, clims = (-0.001,0.001), axis=([], false), subplot=idsp)
        idsp += 1
    end
    display(current())
end # of top_vectors_plot

### Fig. 11a: Haar
top_vectors_plot(dvec_haar, BS_haar, GP, G, iend=10)
savefig("toronto_vv_haar09.png")
savefig("toronto_vv_haar09.pdf")

### Walsh
top_vectors_plot(dvec_walsh, BS_walsh, GP, G, iend=10)
savefig("toronto_vv_walsh09.png")
savefig("toronto_vv_walsh09.pdf")

### Fig. 11b: GHWT_c2f (= Walsh in this case)
top_vectors_plot(dvec_c2f, BS_c2f, GP, G, iend=10)
savefig("toronto_vv_c2f09.png")
savefig("toronto_vv_c2f09.pdf")

### Fig. 11c: GHWT_f2c
top_vectors_plot(dvec_f2c, BS_f2c, GP, G, iend=10)
savefig("toronto_vv_f2c09.png")
savefig("toronto_vv_f2c09.pdf")

### Fig. 11d: eGHWT
top_vectors_plot(dvec_eghwt, BS_eghwt, GP, G, iend=10)
savefig("toronto_vv_eghwt09.png")
savefig("toronto_vv_eghwt09.pdf")

ext: jl


text:
# This file defines the reflection and transmission coefficients of a multilayer

abstract type Polarization end
struct te <: Polarization end
struct tm <: Polarization end

abstract type Structure end

struct Bulk{T <:OptProp,U <:OptProp} <: Structure
    ep1 :: T
    ep2 :: U
end
# Outer constructor for the Bulk type when one medium is vacuum
Bulk(ep2) = Bulk(Cst(1.0+im*0.0),ep2)

struct Layer <: Structure
    material  :: OptProp
    thickness :: Float64 # in meter. For semi-infinite media thickness = 0
end

Layer(material)  = Layer(material,0.0)

struct Film <: Structure
    film :: Layer
    substrate :: Layer
end

# new Vector type
const MultiLayer = Vector{Layer}

# Union types
const LayerOrMultiLayer = Union{Layer,MultiLayer}

# Accessor functions
substrate(structure :: Bulk) = Layer(structure.ep2)
substrate(structure :: MultiLayer) = structure[end]
substrate(structure :: Film) = structure.substrate

top_medium(structure :: Bulk) = Layer(structure.ep1)
top_medium(structure :: MultiLayer) = structure[1]
top_medium(structure :: Film) = Layer(Cst())

OpticalProperties.permittivity(layer::Layer,w) = permittivity(layer.material,w)

# Normal functions

function compute_kz(kx,eps,w)
    sqrt(eps*(w/c0)^2 - kx*kx)
end

function compute_kx(theta_0,eps,w)
    sqrt(eps)*sin(theta_0)*w/c0
end

function compute_ang(kx,eps,w)
    asin(real(kx)/real(sqrt(eps))/(w/c0))
end

function outer_media(b,kx,w)
    eps1 = permittivity(top_medium(b),w)
    k1z  = compute_kz(kx,eps1,w)

    eps2 = permittivity(substrate(b),w)
    k2z  = compute_kz(kx,eps2,w)
    return eps1,k1z,eps2,k2z
end

# Fresnel coefficient of a semi-infinite medium : perpendicular
function rt(pol :: te, eps1,eps2, k0z,k2z,w)
  # k0z = incident wavevector from medium 1
    r = (k0z - k2z)/(k0z + k2z)
    t = (2.0+im*0.0)*k0z/(k0z + k2z)
    return r,t
end

# Fresnel coefficient of a semi-infinite medium : parallel
function rt(pol :: tm, eps1,eps2, k0z,k2z,w)
  # k0z = incident wavevector from medium 1
    r = (eps2*k0z-eps1*k2z)/(eps2*k0z+eps1*k2z)
    t = (2.0+im*0.0)*k0z/(eps2*k0z+eps1*k2z)*sqrt(eps1)*sqrt(eps2)
    return r,t
end

function power_t(pol::te,eps1,eps2, k0z,k2z,t)
    real(k2z)/real(k0z)*abs(t)^2
end

function power_t(pol::tm,eps1,eps2, k0z,k2z,t)
    real(k2z/eps2)/real(k0z/eps1)*abs(eps2/eps1)*abs(t)^2
end

# Fresnel coefficient of a semi-infinite medium : as a function of incidence angle and frequency
function rt(structure :: Bulk, pol :: Polarization, kx ,w)
    eps1 = permittivity(structure.ep1,w)
    eps2 = permittivity(structure.ep2,w)

    k0z   = compute_kz(kx,eps1,w)
    k2z   = compute_kz(kx,eps2,w)

    r,t=rt(pol,eps1,eps2,k0z,k2z,w)
    return r, t
end

# Fresnel coefficient of a multilayered semi-infinite medium : as a function of incidence angle and frequency

function rt(structure :: MultiLayer, pol :: Polarization, kx ,w)
      S   =[1.0+0.0*im  0.0+0.0*im ;
              0.0+0.0*im  1.0+0.0*im ]
    if real(kx)<w/c0 || real(kx)>w/c0
        scattering_matrix!(S,structure,pol,kx,w)
        return S[2,1], S[1,1]
    else
        return 1.0+0.0*im, 0.0+0.0*im
  end
end

function scattering_matrix!(S,structure :: MultiLayer, pol :: Polarization, kx, w)
        for i=2:length(structure)
            eps1 :: Complex{Float64} = permittivity(structure[i-1].material,w)
            eps2 :: Complex{Float64} = permittivity(structure[i].material,w)
            k0z   = compute_kz(kx,eps1,w)
            k2z   = compute_kz(kx,eps2,w)
            r,t = rt(pol, eps1,eps2, k0z,k2z,w)

             S[1,1] =(S[1,1]*t*exp(im*structure[i-1].thickness*k0z))/(1.0+0.0*im - S[1,2]*r*exp(2.0*im*structure[i-1].thickness*k0z))
             S[1,2] =(S[1,2]*exp(2.0*im*structure[i-1].thickness*k0z)-r)/(1.0+0.0*im - S[1,2]*r*exp(2.0*im*structure[i-1].thickness*k0z))
             S[2,1] =(S[1,1]*S[2,2]*r*exp(im*structure[i-1].thickness*k0z))/t+S[2,1]
             S[2,2] =(S[2,2]*exp(im*structure[i-1].thickness*k0z)*(r*S[1,2] + 1.0+0.0*im))/t
        end
    return nothing
end

function rt(structure :: Film , pol :: Polarization , kx, w)
    filmmaterial      = structure.film.material
    substratematerial = structure.substrate.material
    th                = structure.film.thickness
    eps2 = permittivity(Cst(),w)
    eps1 = permittivity(filmmaterial,w)
    eps0 = permittivity(substratematerial,w)

    k0z = compute_kz(kx,eps0,w)
    k1z = compute_kz(kx,eps1,w)
    k2z = compute_kz(kx,eps2,w)

    r10 ,t10 = rt(pol, eps1,eps0, k1z,k0z,w)
    r21 ,t21 = rt(pol, eps2,eps1, k2z,k1z,w)

    r = (r21 + r10*exp(2*im*k1z*th))/(1 + r21*r10*exp(2*im*k1z*th))
    t = (t21*t10*exp(im*k1z*th))/(1 + r21*r10*exp(2*im*k1z*th))

    return r,t
end

function power_rt(structure, pol :: Polarization, kx ,w)

    eps1,k1z,eps2,k2z = outer_media(structure,kx,w)

    r,t=rt(structure,pol,kx,w)
    R = abs(r)^2
    T = power_t(pol,eps1,eps2, k1z,k2z,t)

    return R, T
end


text:
#=
using Revise, Test, Random, LinearAlgebra, RCall, Distributions
using Cytof5
=#
Random.seed!(10)
printDebug = true

@testset "Compile MCMC.metropolis." begin
  lfc(x) = x
  MCMC.metropolis(1.0, lfc, 1.0)
  @test true

  lfc(x::Vector{Float64}) = sum(x)
  MCMC.metropolis(randn(4), lfc, Matrix{Float64}(LinearAlgebra.I, 4, 4))
  @test true
end

@testset "Compile MCMC.TuningParam." begin
  tval_init = 1.23
  t = MCMC.TuningParam(tval_init)
  MCMC.update_tuning_param_default(t, true)
  MCMC.update_tuning_param_default(t, false)
  MCMC.update_tuning_param_default(t, true)
  @test t.value == tval_init
  #@test MCMC.acceptanceRate(t) == 2/3
end

mutable struct State
  x::Int
  y::Float64
  z::Vector{Float64}
end

@testset "Compile MCMC.gibbs." begin
  function update(s::State, i::Int, out::Any)
    s.x += 1
    s.y -= 1
    s.z[1] += 1
  end

  s = State(0, 0, [0,0])
  out, lastState = MCMC.gibbs(s, update, monitors=[[:x,:y], [:z]], thins=[1,2],
                              printFreq=printDebug ? 0 : -1)
  @test true
end

@testset "Test logit / sigmoid" begin
  @test MCMC.logit(.6) ≈ log(.6 / .4)
  @test MCMC.sigmoid(3.) ≈ 1 / (1 + exp(-3.))
end

mutable struct Param
  mu::Vector{Float64}
  sig2::Float64
  cs_sig2::Float64
end

@testset "MCMC.metropolisAdaptive" begin
  #muTrue = [1., 3.]
  #J = length(muTrue)
  J = 3
  muTrue = randn(J) * 10
  sig2True = .5
  nj = 300
  y = [ rand(Normal(m, sqrt(sig2True)), nj) for m in muTrue ]
  n = nj * J

  muPriorMean = 0.
  muPriorSd = 10.

  sig2Prior_a = 3.
  sig2Prior_b = 2.

  tunerMu = [ deepcopy(MCMC.TuningParam(1.0)) for j in 1:J ]
  tunerSig2 = MCMC.TuningParam(1.0)

  function update(s::Param, i::Int, out::Any)

    function update_muj(j::Int)
      function lfc(muj::Float64)
        ll = sum(logpdf.(Normal(muj, sqrt(s.sig2)), y[j]))
        lp = logpdf(Normal(muPriorMean, muPriorSd), muj)
        return ll + lp
      end

      s.mu[j] = MCMC.metropolisAdaptive(s.mu[j], lfc, tunerMu[j])
      return
    end

    function update_mu()
      for j in 1:J
        update_muj(j)
      end
      return
    end

    function update_sig2()

      # If true, test `metLogAdaptive`
      # otherwise, test `metLogitAdaptive`
      test_metlogadaptive = true

      function log_prob(sig2::Float64)
        ll = sum( sum(logpdf.(Normal(s.mu[j], sqrt(sig2)), y[j])) for j in 1:J )
        if test_metlogadaptive
          lp = logpdf(InverseGamma(sig2Prior_a, sig2Prior_b), sig2)
        else
          lp = logpdf(Uniform(0.0, 5.0), sig2)
        end

        return ll + lp
      end

      if test_metlogadaptive
        s.sig2 = MCMC.metLogAdaptive(s.sig2, log_prob, tunerSig2)
      else
        s.sig2 = MCMC.metLogitAdaptive(s.sig2, log_prob, tunerSig2, a=0.0, b=5.0)
      end

      return
    end

    # Update
    update_mu()
    update_sig2()
    s.cs_sig2 = tunerSig2.value + 0

    return
  end

  init = Param(fill(0., J), sig2True, tunerSig2.value)
  @time out, state = MCMC.gibbs(deepcopy(init), update, nmcmc=2000, nburn=20000,
                                printFreq=printDebug ? 0 : -1)

  if printDebug
    # Print acceptance rates
    for j in 1:J
      println("μ$j Acceptance rate: $(MCMC.acceptanceRate(tunerMu[j]))")
      println("μ$j Tuner value:     $(tunerMu[j].value)")
    end
    println("sig2 Acceptance rate:  $(MCMC.acceptanceRate(tunerSig2))")
    println("sig2 Tuner value:      $(tunerSig2.value)")

    # Print Means
    println("muTrue: $(muTrue)")
    muPost = hcat([o[:mu] for o in out[1]]...)
    sig2Post = hcat([o[:sig2] for o in out[1]]...)
    cs_sig2 = hcat([o[:cs_sig2] for o in out[1]]...)
    R"""
    library(rcommon)
    pdf('result/amcmc.pdf')
    plotPosts($(muPost[1:min(4,J), :]'), cnames=paste('truth:', round($muTrue, 3)))
    plotPost($(sig2Post'), main=paste('truth:', $sig2True))
    plot(1:length($cs_sig2), $cs_sig2, type='l')
    dev.off()
    """
  end
end

@testset "Test logsumexp" begin
  x = randn(10)
  @test MCMC.logsumexp(x) == log(sum(exp.(x)))
  @test MCMC.logsumexp(x[1:3]) == MCMC.logsumexp(x[1], x[2], x[3])
end


text:
function compile(net;debug=false,gpu=false,eltype=Float64)
    # (c) David Barber, University College London 2015

    # It's good to ensure that we only compuile on the CPU since then we don't need to write inplace versions of the GPU functions (we only need the inplace versions of the GPU derivatives. This makes coding a bit easier)

    if debug
        println("Compiling into a tree and storing messages:")
    end
    node=net.node
    N=length(node)
#    NN=net.node[end].index
    if isempty(net.FunctionNode)
        net.FunctionNode=N
    end

    # get the computation tree:
    G=zeros(Bool,N,N)
    returnderivative=zeros(Bool,N)
        for i in net.validnodes
            G[node[i].index,node[i].parents]=1
#            if !(node[i].input)
                if HasParents(node[i])
#                node[i].takederivative=any(map( (x)-> any(x.takederivative), node[node[i].parents]) ) # set to true for those parents that require derivative
                node[i].takederivative=any( (x)-> any(x.takederivative), node[node[i].parents] ) # set to true for those parents that require derivative
            end
    end
    for i in net.validnodes
        if !(node[i]==nothing)
            node[i].children=setdiff(find(G[:,i]),node[i].index)
        end
    end

    # forward pass: assume the nodes are defined in ancestral order
    if debug;  println("Forward Pass compilation:");  end

    #TODO: need to fix this for conversion.....
    for n in net.validnodes
        if node[n].isconst==true
            #if typeof(node[n].constval)==Float32
            if isaScalar(node[n].constval)
                net.value[n]=cArray(gpu,node[n].constval*ones(1,1))
            else
                net.value[n]=cArray(gpu,node[n].constval)
            end
        end
    end

    ForwardPassList!(net) # call this as ForwardPassList(net;ExcludeNodes=[1 2 3 ....]) if there are nodes 1 2 3 ... that are not needed to be calculated in the forward pass

    value=net.value
    gradient=net.gradient
    auxvalue=similar(value)
    net.auxvalue=auxvalue
    net.value=value
    net.gradient=gradient
    net.gpu=gpu

    ADforward!(net;debug=debug,AllocateMemory=true)
    if debug;  println("Done forward pass compilation:");  end

    # backward compilation:
    if debug;  println("Backward Pass compilation:");  end

    gradient=Array(Any,N)

    for i in net.validnodes
        gradient[i]=cArray(gpu,Float64,size(value[i]))
    end
    fill!(gradient[net.FunctionNode],1.0)

    Nend=net.FunctionNode
    if Nend!=N
        anc=ancestors(node,Nend)
    else
        anc=flipdim(net.validnodes,1); anc=anc[2:end]
    end

    relevantchildren=Array(Any,N)
    ancUnionNend=union(anc,Nend)
    for n in anc
       relevantchildren[n]=intersect(node[n].children,ancUnionNend)
    end

    parentIDX=Dict()
    for c in net.validnodes
        for n in unique(node[c].parents)
            inds=findin(node[c].parents,n)
            parentIDX[c,n]=inds
        end
    end

    net.ancestors=anc
    net.relevantchildren=relevantchildren
    net.gradient=gradient
    net.parentIDX=parentIDX

    if debug
        println("Compiled into a DAG with $(length(node)) nodes")
    end

    return net

end

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
include("types.jl")

_DIR = "stochastic_process"
include(joinpath(_DIR, "bounded_diffusion_domain.jl"))
include(joinpath(_DIR, "guid_prop_bridge.jl"))
include(joinpath(_DIR, "path_to_wiener.jl"))

_DIR = "solvers"
include(joinpath(_DIR, "vern7.jl"))
include(joinpath(_DIR, "tsit5.jl"))
include(joinpath(_DIR, "rk4.jl"))
include(joinpath(_DIR, "ralston3.jl"))
include(joinpath(_DIR, "euler_maruyama_dom_restr.jl"))

_DIR = "transition_kernels"
include(joinpath(_DIR, "random_walk.jl"))

_DIR = "mcmc_extras"
include(joinpath(_DIR, "adaptation.jl"))
include(joinpath(_DIR, "blocking_schedule.jl"))
include(joinpath(_DIR, "first_passage_times.jl"))
include(joinpath(_DIR, "starting_pt.jl"))

_DIR = "examples"
include(joinpath(_DIR, "fitzhugh_nagumo.jl"))
include(joinpath(_DIR, "radial_ornstein_uhlenbeck.jl"))
include(joinpath(_DIR, "lorenz_system.jl"))
include(joinpath(_DIR, "lorenz_system_const_vola.jl"))

_DIR = "mcmc"
include(joinpath(_DIR, "priors.jl"))
include(joinpath(_DIR, "setup.jl"))
include(joinpath(_DIR, "workspace.jl"))
include(joinpath(_DIR, "conjugate_updt.jl"))
include(joinpath(_DIR, "mcmc.jl"))
include(joinpath(_DIR, "repeated.jl"))


text:
# using Revise
using IMPES
using JPhreeqc
# using Plots, LinearAlgebra, IMPES
# using AnalyticalEOR

#* Define Relative Perms
kr = RelPerms(
        swr= 0.2,
        sor= 0.2,
        krw0= 0.3,
        kro0= 1.0,
        nw= 2.0,
        no= 3.0,
            )

#* Define reservoir
res = Reservoir(
    L = 300.,
    W = 200.0 * 0.3048,
    h = 50.0 * 0.3048,
    k = 100.0 * 9.869233e-16,
    phi = 0.2,
    kr = kr,
    )

#* Define fluids
fluid = Fluids(
            μw=1e-3,
            μo=5e-3,
)

#* Define grid
grid = Grid(
            res=res,
            Nx = 100,
            )

#* Define BC
bc = BoundaryConditions(
                        [:neumann, :dirichlet],
                        [426.5 * 0.3048 ^ 3 / 3600 / 24, 1000 * 6894.76]
                        )

tmax = 4000. * 3600 * 24 # days
dt = 1. * 3600 * 24

#* Initialize simulation and postprocessing
sim = Simulation(
                1000*ones(grid.Nx)* 6894.76,
                0.2*ones(grid.Nx),
                tmax,
                dt
)

post = FlowResults(sim, grid, tmax, dt)

ncomps = 2
cinj = ones(ncomps) * 1.
ci = zeros(ncomps, grid.Nx)
ntimes = ceil(Int64, sim.tmax / sim.dt)
c = zeros(ncomps, grid.Nx, ntimes+1)
c[:,:, 1] = ci
ctrans = CompTransport(c, ncomps, cinj)

while sim.t < tmax
    IMPES!(sim, res, fluid, grid, bc, post)
    IMPEC!(ctrans, sim, res, grid, post)
end

using Plots

i = 1000
plot(grid.x, post.s[:,i], seriestype=:scatter)
plot(grid.x, ctrans.c[1,:,i], seriestype=:scatter)
plot!(grid.x, ctrans.c[2,:,i], seriestype=:scatter)

max_stars_repo_path: examples/Ex1 - BL1D with PHREEQC.jl

max_issues_repo_path: examples/Ex1 - BL1D with PHREEQC.jl

max_forks_repo_path: examples/Ex1 - BL1D with PHREEQC.jl


text:
""" Returns the convolution of 2 arrays, arr1 and arr2
Pre: arr1 and arr2 must not be empty.
"""
function convolve(arr1, arr2)
    if length(arr1) > length(arr2)
        larger_arr = arr1
        smaller_arr = arr2
    else
        larger_arr = arr2
        smaller_arr = arr1
    end

    number_of_columns = length(arr1) + length(arr2) - 1
    conv_matrix = zeros(length(larger_arr), number_of_columns)
    for i in 1:length(larger_arr)
        mulitplied_row = smaller_arr * larger_arr[i]
        # Insert the multiplied row into the convolution matrix
        conv_matrix[CartesianIndex.(i, i:length(mulitplied_row)+i-1)] .= mulitplied_row
    end

    convolution = sum(conv_matrix, dims=1)  # Find sum of columns
    return convolution
end

""" Returns the product of 2 matrices.
Pre: The number of columns in matrix1 must equal the number of rows in matrix2 (throws
DimensionMismatch if not).
Post: The resulting matrix will have the same number of rows as matrix1 and the same number of
columns as matrix2.
"""
function matrixMultiply(matrix1, matrix2)
    # Handle the case when 1 matrix is empty and the other has a length of 1, or when both matrices
    # are emtpy.
    if length(matrix1) == 0 && length(matrix2) == 1 ||
            length(matrix1) == 1 && length(matrix2) == 0 ||
            length(matrix1) == 0 && length(matrix2) == 0
        return reshape([], 0, 0)
    elseif size(matrix1, 2) != size(matrix2, 1)
        throw(DimensionMismatch("The number of columns in matrix1 must equal the number of rows in matrix2"))
    end

    resultMatrix = zeros(size(matrix1, 1), size(matrix2, 2))
    for row in 1:size(matrix1, 1)
        for column in 1:size(matrix2, 2)
            row1 = matrix1[row, :]
            column2 = matrix2[:, column]
            resultMatrix[row, column] = sum(dotProduct(row1, column2))
        end
    end
    return resultMatrix
end

#  Returns the dot product of array1 and array2
function dotProduct(array1, array2)
    return array1 .* array2
end

""""""
function DFT(array)

end

""""""
function IDFT(array)

end


text:
using PyPlot

"""
    plotUsage(sv, colors=["b","g","r","y"])

Plot the element usage for a multi-element DAC.
The colors are for `sv` = 1,-1,i,-i.
"""
function plotUsage(sv, colors=["b","g","r","y"])
    T = size(sv, 2)
    M = size(sv, 1)

    # Plot the grid
    x = repeat(0:T, inner=2)
    T2 = ceil(Int, (T+1)/2)
    y = [zeros(1, T2); M*ones(2, T2); zeros(1,T2)]
    y = y[1:2*(T+1)]
    plot(x, y, color="k")
    M2 = ceil(Int, (M+1)/2)
    x = [zeros(1, M2); T*ones(2, M2); zeros(1, M2)]
    x = x[1:2*(M+1)]
    y = repeat(0:M, inner=2)
    plot(x, y, color="k")
    axis("image")

    for t = 1:T
        for i = 1:M
            if sv[i,t] == 1
                fill([t-1, t-1, t, t], [i-1, i, i, i-1], color=colors[1])
            elseif sv[i,t] == -1
                fill([t-1, t-1, t, t], [i-1, i, i, i-1], color=colors[2])
            elseif sv[i,t] == 1im
                fill([t-1, t-1, t, t], [i-1, i, i, i-1], color=colors[3])
            elseif sv[i,t] == -1im
                fill([t-1, t-1, t, t], [i-1, i, i, i-1], color=colors[4])
            end
        end
    end
end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
import Base: *

*(lhs::Bool, rhs::Sigmoid{N, ES, mode}) where {N, ES, mode} = reinterpret(Sigmoid{N, ES, mode}, @s(rhs) * lhs)
*(lhs::Sigmoid{N, ES, mode}, rhs::Bool) where {N, ES, mode} = reinterpret(Sigmoid{N, ES, mode}, @s(lhs) * rhs)

#############################################################################3##
## return types for valid division.

const multiplication_types = Dict((:guess, :guess) => :guess,
                                  (:inward_exact,  :inward_exact)  => :inward_exact,
                                  (:inward_ulp,    :inward_exact)  => :inward_ulp,
                                  (:inward_exact,  :inward_ulp)    => :inward_ulp,
                                  (:inward_ulp,    :inward_ulp)    => :inward_ulp,
                                  (:outward_exact, :outward_exact) => :outward_exact,
                                  (:outward_ulp,   :outward_exact) => :outward_ulp,
                                  (:outward_exact, :outward_ulp)   => :outward_ulp,
                                  (:outward_ulp,   :outward_ulp)   => :outward_ulp)

const nanerror_code = :(throw(NaNError(*, Any[lhs, rhs])))

const multiplication_inf_zero = Dict((:guess,         :guess)         => nanerror_code,
                                     (:inward_exact,  :inward_exact)  => nanerror_code,
                                     (:inward_ulp,    :inward_exact)  => :(zero(Sigmoid{N,ES,:inward_exact})),
                                     (:inward_exact,  :inward_ulp)    => :(Sigmoid{N,ES,:inward_exact}(Inf)),
                                     (:inward_ulp,    :inward_ulp)    => :(zero(Sigmoid{N,ES,:inward_ulp})),
                                     (:outward_exact, :outward_exact) => nanerror_code,
                                     (:outward_ulp,   :outward_exact) => :(zero(Sigmoid{N,ES,:outward_exact})),
                                     (:outward_exact, :outward_ulp)   => :(Sigmoid{N,ES,:outward_exact}(Inf)),
                                     (:outward_ulp,   :outward_ulp)   => :(Sigmoid{N,ES,:outward_ulp}(Inf)))

const multiplication_left_inf = Dict((:guess, :guess)                 => :(Sigmoid{N,ES,:guess}(Inf)),
                                     (:inward_exact,  :inward_exact)  => :(Sigmoid{N,ES,:inward_exact}(Inf)),
                                     (:inward_ulp,    :inward_exact)  => :(Sigmoid{N,ES,:inward_ulp}(Inf)),
                                     (:inward_exact,  :inward_ulp)    => :(Sigmoid{N,ES,:inward_exact}(Inf)),
                                     (:inward_ulp,    :inward_ulp)    => :(Sigmoid{N,ES,:inward_ulp}(Inf)),
                                     (:outward_exact, :outward_exact) => :(Sigmoid{N,ES,:outward_exact}(Inf)),
                                     (:outward_ulp,   :outward_exact) => :(Sigmoid{N,ES,:outward_ulp}(Inf)),
                                     (:outward_exact, :outward_ulp)   => :(Sigmoid{N,ES,:outward_exact}(Inf)),
                                     (:outward_ulp,   :outward_ulp)   => :(Sigmoid{N,ES,:outward_ulp}(Inf)))

const multiplication_left_zero = Dict((:guess,         :guess)         => :(zero(Stype)),
                                      (:inward_exact,  :inward_exact)  => :(zero(Sigmoid{N,ES,:inward_exact})),
                                      (:inward_ulp,    :inward_exact)  => :(zero(Sigmoid{N,ES,:inward_ulp})),
                                      (:inward_exact,  :inward_ulp)    => :(zero(Sigmoid{N,ES,:inward_exact})),
                                      (:inward_ulp,    :inward_ulp)    => :(zero(Sigmoid{N,ES,:inward_ulp})),
                                      (:outward_exact, :outward_exact) => :(zero(Sigmoid{N,ES,:outward_exact})),
                                      (:outward_ulp,   :outward_exact) => :(zero(Sigmoid{N,ES,:outward_ulp})),
                                      (:outward_exact, :outward_ulp)   => :(zero(Sigmoid{N,ES,:outward_exact})),
                                      (:outward_ulp,   :outward_ulp)   => :(zero(Sigmoid{N,ES,:outward_ulp})))

@generated function *(lhs::Sigmoid{N, ES, lhs_mode}, rhs::Sigmoid{N, ES, rhs_mode}) where {N, ES, lhs_mode, rhs_mode}

    #dealing with modes for multiplication
    if !haskey(multiplication_types, (lhs_mode, rhs_mode))

        #inward_exact and outward_exact could be annihilators which makes some processes
        #try to send "incorrect types" to the multiplication algorithm.  This segment
        #of metacode deals with these situations, and creates stub multiplication
        #algorithms which handle these cases.

        rhs_infcheck,rhs_zercheck = (rhs_mode == :inward_exact) || (rhs_mode == :outward_exact) ? (:(!isfinite(rhs)),:(Base.iszero(rhs))) : (:(false), :(false))

        lhs_test = if (lhs_mode == :inward_exact) || (lhs_mode == :outward_exact)
            quote
                if Base.iszero(lhs)
                    $rhs_infcheck && $nanerror_code
                    return zero(Sigmoid{N,ES,:inward_exact})
                end
                if !isfinite(lhs)
                    $rhs_zercheck && $nanerror_code
                    return Sigmoid{N,ES,:outward_exact}(Inf)
                end
            end
        else
            :()
        end

        rhs_test = if (rhs_mode == :inward_exact) || (rhs_mode == :outward_exact)
            quote
                Base.iszero(rhs) && return zero(Sigmoid{N,ES,:inward_exact})
                !isfinite(rhs) && return Sigmoid{N,ES,:outward_exact}(Inf)
            end
        else
            :()
        end

        return quote
            $lhs_test
            $rhs_test
            throw(ArgumentError("incompatible types"))
        end
    end

    #throw(ArgumentError("incompatible types passed to multiplication function! ($lhs_mode, $rhs_mode)"))
    mode = multiplication_types[(lhs_mode, rhs_mode)]
    infzero_code = multiplication_inf_zero[(lhs_mode,rhs_mode)]
    zeroinf_code = multiplication_inf_zero[(rhs_mode,lhs_mode)]
    leftinf_code = multiplication_left_inf[(lhs_mode,rhs_mode)]
    rightinf_code = multiplication_left_inf[(rhs_mode,lhs_mode)]
    leftzero_code = multiplication_left_zero[(lhs_mode,rhs_mode)]
    rightzero_code = multiplication_left_zero[(rhs_mode,lhs_mode)]

    S = Sigmoid{N,ES,mode}

    quote
        Stype = $S
        #multiplying infinities is infinite.
        if !isfinite(lhs)
            (reinterpret((@UInt), rhs) == zero(@UInt)) && return $infzero_code
            return $leftinf_code
        end
        if !isfinite(rhs)
            (reinterpret((@UInt), lhs) == zero(@UInt)) && return $zeroinf_code
            return $rightinf_code
        end
        #mulitplying zeros is zero
        (reinterpret((@UInt), lhs) == zero(@UInt)) && return $leftzero_code
        (reinterpret((@UInt), rhs) == zero(@UInt)) && return $rightzero_code
        return mul_algorithm(lhs, rhs)
    end
end

@generated function mul_algorithm(lhs::Sigmoid{N, ES, lhs_mode}, rhs::Sigmoid{N, ES, rhs_mode}) where {N, ES, lhs_mode, rhs_mode}

    #dealing with modes for multiplication
    haskey(multiplication_types, (lhs_mode, rhs_mode)) || throw(ArgumentError("incompatible types passed to multiplication function!"))
    mode = multiplication_types[(lhs_mode, rhs_mode)]

    S = Sigmoid{N,ES,mode}
    quotemode = QuoteNode(mode)

    quote
        #generate the lhs and rhs subcomponents.
        mode = $quotemode
        @breakdown lhs
        @breakdown rhs
        #sign is the xor of both signs.
        mul_sgn = lhs_sgn ⊻ rhs_sgn
        #the multiplicative exponent is the sum of the two exponents.
        mul_exp = lhs_exp + rhs_exp
        lhs_frc = (lhs_frc >> 1) | (@signbit)
        rhs_frc = (rhs_frc >> 1) | (@signbit)
        #then calculate the fraction.
        mul_frc = demote(promote(lhs_frc) * promote(rhs_frc))
        shift = leading_zeros(mul_frc)
        mul_frc <<= shift + 1
        mul_exp -= (shift - 1)
        __round(build_numeric($S, mul_sgn, mul_exp, mul_frc))
    end
end

################################################################################
## division

const division_types = Dict((:guess, :guess) => :guess,
                            (:exact, :exact) => :exact,  #these three are supported for inverse calculation.
                            (:exact, :upper) => :lower,
                            (:exact, :lower) => :upper,
                            (:inward_exact,  :outward_exact)  => :inward_exact,
                            (:inward_ulp,    :outward_exact)  => :inward_ulp,
                            (:inward_exact,  :outward_ulp)    => :inward_ulp,
                            (:inward_ulp,    :outward_ulp)    => :inward_ulp,
                            (:outward_exact, :inward_exact) => :outward_exact,
                            (:outward_ulp,   :inward_exact) => :outward_ulp,
                            (:outward_exact, :inward_ulp)   => :outward_ulp,
                            (:outward_ulp,   :inward_ulp)   => :outward_ulp)

const division_left_inf = Dict((:guess, :guess)                 => :(Sigmoid{N,ES,:guess}(Inf)),
                               (:exact, :exact)                 => :(Sigmoid{N,ES,:exact}(Inf)),  #these three are supported for inverse calculation.
                               (:exact, :upper)                 => :(Sigmoid{N,ES,:exact}(Inf)),
                               (:exact, :lower)                 => :(Sigmoid{N,ES,:exact}(Inf)),
                               (:inward_exact,  :outward_exact) => :(Sigmoid{N,ES,:inward_exact}(Inf)),
                               (:inward_ulp,    :outward_exact) => :(Sigmoid{N,ES,:inward_ulp}(Inf)),
                               (:inward_exact,  :outward_ulp)   => :(Sigmoid{N,ES,:inward_exact}(Inf)),
                               (:inward_ulp,    :outward_ulp)   => :(Sigmoid{N,ES,:inward_ulp}(Inf)),
                               (:outward_exact, :inward_exact)  => :(Sigmoid{N,ES,:outward_exact}(Inf)),
                               (:outward_ulp,   :inward_exact)  => :(nothing),
                               (:outward_exact, :inward_ulp)    => :(Sigmoid{N,ES,:outward_exact}(Inf)),
                               (:outward_ulp,   :inward_ulp)    => :(nothing))

const division_right_inf = Dict((:guess, :guess)                 => :(zero(Sigmoid{N,ES,:guess})),
                                (:exact, :exact)                 => :(zero(Sigmoid{N,ES,:exact})),  #these three are supported for inverse calculation.
                                (:exact, :upper)                 => :(zero(Sigmoid{N,ES,:lower})),
                                (:exact, :lower)                 => :(zero(Sigmoid{N,ES,:upper})),
                                (:upper, :exact)                 => :(zero(Sigmoid{N,ES,:exact})),
                                (:lower, :exact)                 => :(zero(Sigmoid{N,ES,:exact})),
                                (:inward_exact,  :outward_exact) => :(zero(Sigmoid{N,ES,:inward_exact})),
                                (:inward_ulp,    :outward_exact) => :(zero(Sigmoid{N,ES,:inward_exact})),
                                (:inward_exact,  :outward_ulp)   => :(nothing),
                                (:inward_ulp,    :outward_ulp)   => :(nothing),
                                (:outward_exact, :inward_exact)  => :(zero(Sigmoid{N,ES,:outward_exact})),
                                (:outward_ulp,   :inward_exact)  => :(zero(Sigmoid{N,ES,:outward_exact})),
                                (:outward_exact, :inward_ulp)    => :(zero(Sigmoid{N,ES,:outward_ulp})),
                                (:outward_ulp,   :inward_ulp)    => :(zero(Sigmoid{N,ES,:outward_ulp})))

const division_zero_inf =  Dict((:guess, :guess)                 => :(zero(Sigmoid{N,ES,:guess})),
                                (:exact, :exact)                 => :(zero(Sigmoid{N,ES,:exact})),  #these three are supported for inverse calculation.
                                (:exact, :upper)                 => :(zero(Sigmoid{N,ES,:lower})),
                                (:exact, :lower)                 => :(zero(Sigmoid{N,ES,:upper})),
                                (:upper, :exact)                 => :(zero(Sigmoid{N,ES,:exact})),
                                (:lower, :exact)                 => :(zero(Sigmoid{N,ES,:exact})),
                                (:inward_exact,  :outward_exact) => :(zero(Sigmoid{N,ES,:inward_exact})),
                                (:inward_ulp,    :outward_exact) => :(nothing),
                                (:inward_exact,  :outward_ulp)   => :(zero(Sigmoid{N,ES,:inward_exact})),
                                (:inward_ulp,    :outward_ulp)   => :(nothing),
                                (:outward_exact, :inward_exact)  => :(zero(Sigmoid{N,ES,:inward_exact})),
                                (:outward_ulp,   :inward_exact)  => :(zero(Sigmoid{N,ES,:outward_exact})),
                                (:outward_exact, :inward_ulp)    => :(zero(Sigmoid{N,ES,:outward_exact})),
                                (:outward_ulp,   :inward_ulp)    => :(zero(Sigmoid{N,ES,:outward_ulp})))

const division_both_inf = Dict((:guess, :guess)                 => :(throw(NaNError(/,[lhs, rhs]))),
                               (:exact, :exact)                 => :(throw(NaNError(/,[lhs, rhs]))),  #these three are supported for inverse calculation.
                               (:exact, :upper)                 => :(Sigmoid{N,ES,:exact}(Inf)),
                               (:exact, :lower)                 => :(Sigmoid{N,ES,:exact}(Inf)),
                               (:upper, :exact)                 => :(zero(Sigmoid{N,ES,:exact})),
                               (:lower, :exact)                 => :(zero(Sigmoid{N,ES,:exact})),
                               (:inward_exact,  :outward_exact) => :(throw(NaNError(/,[lhs, rhs]))),
                               (:inward_ulp,    :outward_exact) => :(zero(Sigmoid{N,ES,:inward_exact})),
                               (:inward_exact,  :outward_ulp)   => :(Sigmoid{N,ES,:inward_exact}(Inf)),
                               (:inward_ulp,    :outward_ulp)   => :(nothing),
                               (:outward_exact, :inward_exact)  => :(throw(NaNError(/,[lhs, rhs]))),
                               (:outward_ulp,   :inward_exact)  => :(zero(Sigmoid{N,ES,:outward_exact})),
                               (:outward_exact, :inward_ulp)    => :(Sigmoid{N,ES,:outward_exact}(Inf)),
                               (:outward_ulp,   :inward_ulp)    => :(nothing))

const division_both_zero = Dict((:guess, :guess)                 => :(throw(NaNError(/,[lhs, rhs]))),
                                (:exact, :exact)                 => :(throw(NaNError(/,[lhs, rhs]))),  #these three are supported for inverse calculation.
                                (:upper, :exact)                 => :(Sigmoid{N,ES,:exact}(Inf)),
                                (:lower, :exact)                 => :(Sigmoid{N,ES,:exact}(Inf)),
                                (:exact, :upper)                 => :(zero(Sigmoid{N,ES,:exact})),
                                (:exact, :lower)                 => :(zero(Sigmoid{N,ES,:exact})),
                                (:inward_exact,  :outward_exact) => :(throw(NaNError(/,[lhs, rhs]))),
                                (:inward_ulp,    :outward_exact) => :(Sigmoid{N,ES,:inward_exact}(Inf)),
                                (:inward_exact,  :outward_ulp)   => :(zero(Sigmoid{N,ES,:inward_exact})),
                                (:inward_ulp,    :outward_ulp)   => :(nothing),
                                (:outward_exact, :inward_exact)  => :(throw(NaNError(/,[lhs, rhs]))),
                                (:outward_ulp,   :inward_exact)  => :(Sigmoid{N,ES,:outward_exact}(Inf)),
                                (:outward_exact, :inward_ulp)    => :(zero(Sigmoid{N,ES,:outward_exact})),
                                (:outward_ulp,   :inward_ulp)    => :(nothing))

const division_left_zero = Dict((:guess, :guess)                 => :(zero(Sigmoid{N,ES,:guess})),
                                (:exact, :exact)                 => :(zero(Sigmoid{N,ES,:exact})),  #these three are supported for inverse calculation.
                                (:exact, :upper)                 => :(zero(Sigmoid{N,ES,:exact})),
                                (:exact, :lower)                 => :(zero(Sigmoid{N,ES,:exact})),
                                (:upper, :exact)                 => :(zero(Sigmoid{N,ES,:upper})),
                                (:lower, :exact)                 => :(zero(Sigmoid{N,ES,:lower})),
                                (:inward_exact,  :outward_exact) => :(zero(Sigmoid{N,ES,:inward_exact})),
                                (:inward_ulp,    :outward_exact) => :(nothing),
                                (:inward_exact,  :outward_ulp)   => :(zero(Sigmoid{N,ES,:inward_exact})),
                                (:inward_ulp,    :outward_ulp)   => :(nothing),
                                (:outward_exact, :inward_exact)  => :(zero(Sigmoid{N,ES,:outward_exact})),
                                (:outward_ulp,   :inward_exact)  => :(zero(Sigmoid{N,ES,:outward_ulp})),
                                (:outward_exact, :inward_ulp)    => :(zero(Sigmoid{N,ES,:outward_exact})),
                                (:outward_ulp,   :inward_ulp)    => :(zero(Sigmoid{N,ES,:outward_ulp})))

const division_right_zero = Dict((:guess, :guess)                 => :(Sigmoid{N,ES,:guess}(Inf)),
                                 (:exact, :exact)                 => :(Sigmoid{N,ES,:exact}(Inf)),  #these three are supported for inverse calculation.
                                 (:exact, :upper)                 => :(Sigmoid{N,ES,:lower}(Inf)),
                                 (:exact, :lower)                 => :(Sigmoid{N,ES,:upper}(Inf)),
                                 (:upper, :exact)                 => :(Sigmoid{N,ES,:exact}(Inf)),
                                 (:lower, :exact)                 => :(Sigmoid{N,ES,:exact}(Inf)),
                                 (:inward_exact,  :outward_exact) => :(Sigmoid{N,ES,:inward_exact}(Inf)),
                                 (:inward_ulp,    :outward_exact) => :(Sigmoid{N,ES,:inward_exact}(Inf)),
                                 (:inward_exact,  :outward_ulp)   => :(Sigmoid{N,ES,:inward_ulp}(Inf)),
                                 (:inward_ulp,    :outward_ulp)   => :(Sigmoid{N,ES,:inward_ulp}(Inf)),
                                 (:outward_exact, :inward_exact)  => :(Sigmoid{N,ES,:outward_exact}(Inf)),
                                 (:outward_ulp,   :inward_exact)  => :(Sigmoid{N,ES,:outward_exact}(Inf)),
                                 (:outward_exact, :inward_ulp)    => :(nothing),
                                 (:outward_ulp,   :inward_ulp)    => :(nothing))

@generated function Base.:/(lhs::Sigmoid{N,ES,lhs_mode}, rhs::Sigmoid{N,ES,rhs_mode}) where {N, ES, lhs_mode, rhs_mode}
     if !haskey(division_types, (lhs_mode, rhs_mode))
        #inward_exact and outward_exact could be annihilators which makes some processes
        #try to send "incorrect types" to the division algorithm.  This segment
        #of metacode deals with these situations, and creates stub division
        #algorithms which handle these cases.

        rhs_infcheck,rhs_zercheck = (rhs_mode == :inward_exact) || (rhs_mode == :outward_exact) ? (:(!isfinite(rhs)),:(Base.iszero(rhs))) : (:(false), :(false))

        lhs_test = if (lhs_mode == :inward_exact) || (lhs_mode == :outward_exact)
            quote
                if Base.iszero(lhs)
                    $rhs_zercheck && $nanerror_code
                    return zero(Sigmoid{N,ES,:inward_exact})
                end
                if !isfinite(lhs)
                    $rhs_infcheck && $nanerror_code
                    return Sigmoid{N,ES,:outward_exact}(Inf)
                end
            end
        else
            :()
        end

        rhs_test = if (rhs_mode == :inward_exact) || (rhs_mode == :outward_exact)
            quote
                Base.iszero(rhs) && return Sigmoid{N,ES,:outward_exact}(Inf)
                !isfinite(rhs) && return zero(Sigmoid{N,ES,:inward_exact})
            end
        else
            :()
        end

        return quote
            $lhs_test
            $rhs_test
            throw(ArgumentError("incompatible types"))
        end
    end

    mode = division_types[(lhs_mode, rhs_mode)]
    quotemode = QuoteNode(mode)

    #specialized code for certain results
    division_left_inf_code  = division_left_inf[(lhs_mode, rhs_mode)]
    division_right_inf_code = division_right_inf[(lhs_mode, rhs_mode)]
    division_both_inf_code  = division_both_inf[(lhs_mode, rhs_mode)]
    division_both_zero_code = division_both_zero[(lhs_mode, rhs_mode)]
    division_left_zero_code = division_left_zero[(lhs_mode, rhs_mode)]
    division_right_zero_code = division_right_zero[(lhs_mode, rhs_mode)]
    division_zero_inf_code  = division_zero_inf[(lhs_mode, rhs_mode)]

    #calculate the number of rounds we should apply the goldschmidt method.
    rounds = Int(ceil(log(2,N))) + 1
    top_bit = promote(one(@UInt) << (__BITS - 1))
    bot_bit = (one(@UInt) << (__BITS - N - 1))

    quote
        mode = $quotemode
        #ZERO AND INFINITY EXCEPTION CASES.

        #dividing infinities or by zero is infinite.
        if !isfinite(lhs)
            isfinite(rhs) || return $division_both_inf_code
            return $division_left_inf_code
        end

        if rhs == zero(Sigmoid{N, ES, rhs_mode})
            (lhs == zero(Sigmoid{N, ES, lhs_mode})) && return $division_both_zero_code
            return $division_right_zero_code
        end

        #dividing zeros or by infinity is zero
        if !isfinite(rhs)
            Base.iszero(lhs) || return $division_right_inf_code
            return $division_zero_inf_code
        end

        lhs == zero(Sigmoid{N, ES, lhs_mode}) && return $division_left_zero_code

    cq_mask = promote(-one(@UInt))

    #generate the lhs and rhs subcomponents.  Unlike multiplication, however,
    #we want there to 'always be a hidden bit', so we should use the "numeric" method.
    @breakdown lhs numeric
    @breakdown rhs numeric

    #sign is the xor of both signs.
    div_sgn = lhs_sgn ⊻ rhs_sgn

    #do something different if rhs_frc is zero (aka power of two)
    if rhs_frc == 0
      div_exp = lhs_exp - rhs_exp
      return __round(build_numeric(Sigmoid{N, ES, mode}, div_sgn, div_exp, lhs_frc))
    end

    #the multiplicative exponent is the product of the two exponents.
    div_exp = lhs_exp - rhs_exp - 1

    #calculate the number of zeros in the solution.
    lhs_zeros = trailing_zeros(lhs_frc) - (__BITS - N)
    rhs_zeros = trailing_zeros(rhs_frc) - (__BITS - N)

    cumulative_quotient = promote(lhs_frc)
    cumulative_zpower   = promote(-rhs_frc) >> 1
    power_gain = 0

    #then calculate the fraction, using binomial goldschmidt.
    # binomial goldschmidt algorithm:  x ∈ [1, 2), y ∈ [0.5, 1)
    #   define z ≡ 1 - y ⇒ y == 1 - z.
    #   Note (1 - z)(1 + z)(1 + z^2)(1 + z^4) == (1 - z^2n) → 1
    for rd = 1:($rounds - 1)
      #update the quotient
      cumulative_quotient += ((cumulative_quotient * cumulative_zpower) >> __BITS) + cumulative_zpower
      cumulative_zpower = (cumulative_zpower ^ 2) >> __BITS
      shift = __BITS - leading_zeros(cumulative_quotient)
      if (shift > 0)
        cumulative_quotient &= cq_mask
        cumulative_quotient >>= 1
        power_gain += 1
        (shift == 2) && (cumulative_quotient |= $top_bit)
      end
    end
    #update the cumulative quotient one last time.
    cumulative_quotient += ((cumulative_quotient * cumulative_zpower) >> __BITS) + cumulative_zpower
    shift = __BITS - leading_zeros(cumulative_quotient)
    if (shift > 0)
      cumulative_quotient &= cq_mask
      cumulative_quotient >>= 1
      power_gain += 1
      (shift == 2) && (cumulative_quotient |= $top_bit)
    end

    div_frc = demoteright(cumulative_quotient)

    result_ones = trailing_ones(div_frc >> (__BITS - N))

    if (lhs_zeros + N + (result_ones != N) == result_ones + rhs_zeros)
      #increment the lowest bit
      div_frc += $bot_bit
      #mask out all the other ones that were trailing.
      div_frc &= -$bot_bit

      #bump it if we triggered a carry.
      power_gain += (div_frc == 0)
    end

    __round(build_numeric(Sigmoid{N, ES, mode}, div_sgn, div_exp + power_gain, div_frc))
  end
end

Base.inv(x::Sigmoid{N,ES,mode}) where {N, ES, mode} = one(Sigmoid{N,ES,mode}) / x
Base.inv(x::Sigmoid{N,ES,:lower}) where {N, ES} = one(Sigmoid{N,ES,:exact}) / x
Base.inv(x::Sigmoid{N,ES,:upper}) where {N, ES} = one(Sigmoid{N,ES,:exact}) / x

max_issues_repo_name: MohHizzani/SigmoidNumbers

max_forks_repo_name: MohHizzani/SigmoidNumbers


text:
import jInv.ForwardShare
import jInv.Mesh
import Base.Test
import jInv.Utils

@everywhere begin
  using jInv.ForwardShare
  using jInv.Mesh
  using Test
  using jInv.Utils
  mutable struct TestProb <: ForwardProbType
  	Mesh::AbstractMesh
  end
end

Minv = getRegularMesh([0 1 0 2 0 3],[12 14 16])
Mfor = getRegularMesh([0 1 0 2 0 3],[6 7 8])
pFor = TestProb(Mfor)
pFors = Array{RemoteChannel}(undef,2)
workerList  = workers()
nw          = nworkers()
pFors[1]    = initRemoteChannel(identity,workerList[1%nw+1],pFor)
pFors[2]    = initRemoteChannel(identity,workerList[2%nw+1],pFor)

# prepare single mesh2mesh
M2Mc = prepareMesh2Mesh(pFor,Minv,true) # compact
M2Mf = prepareMesh2Mesh(pFor,Minv,false) # not compact
M2Mrc = prepareMesh2Mesh(pFors,Minv,true) # compact
M2Mrf = prepareMesh2Mesh(pFors,Minv,false) # compact

# compare remote versions with local ones
@test all(fetch(M2Mrc[1]).==M2Mc)
@test all(fetch(M2Mrf[1]).==M2Mf)

# test interpGlobal2Local
x   = randn(Minv.nc)
xf1 = interpGlobalToLocal(x,M2Mc)
xf2 = interpGlobalToLocal(x,M2Mf)
@test norm(xf1-xf2)/norm(xf2) < 1e-12

# test interpLocal2Global
x = randn(Mfor.nc)
xf1 = interpLocalToGlobal(x,M2Mc)
xf2 = interpLocalToGlobal(x,M2Mf)
@test norm(xf1-xf2)/norm(xf2) < 1e-12


text:
export l1Reg

"""
l1Reg(m,mref,M)

l1-regularizer: R = sum(abs(m-mref))
"""
function l1Reg(m::Vector{Float64},mref::Vector{Float64},M::AbstractMesh)
	dm   = m .- mref
    dR   = ones(size(m))
    Rc   = sum(dm)
    return Rc,dR,1e-4*speye(length(m))
end


text:
############################################
## default convention: clockwise rotation ##
############################################

## type declarations
##------------------

immutable Rot0 <: CWRotation
    description::String
    degrees::Int

    function Rot0(descr::String, degr::Int)
        if degr != 0
            error("Wrong degrees for zero rotation!")
        end
        new(descr, degr)
    end

    Rot0() = Rot0("No rotation", 0)
end

immutable Rot90 <: CWRotation
    description::String
    degrees::Int

    function Rot90(descr::String, degr::Int)
        if (degr != 90)
            error("Wrong degrees for 90 degrees rotation!")
        end
        new(descr, degr)
    end

    Rot90() = Rot90("Clockwise rotation of 90 degrees", 90)
end

immutable Rot180 <: CWRotation
    description::String
    degrees::Int

    function Rot180(descr::String, degr::Int)
        if (degr != 180)
            error("Wrong degrees for 180 degrees rotation!")
        end
        new(descr, degr)
    end

    Rot180() = Rot180("Clockwise rotation of 180 degrees", 180)
end

immutable Rot270 <: CWRotation
    description::String
    degrees::Int

    function Rot270(descr::String, degr::Int)
        if (degr != 270)
            error("Wrong degrees for 270 degrees rotation!")
        end
        new(descr, degr)
    end

    Rot270() = Rot270("Clockwise rotation of 270 degrees", 270)
end

#######################################################
## resolving pdf and cdf arguments for modifications ##
#######################################################

## clockwise rotation, array of points
##------------------------------------

function resolve_args(rot::Rot0,
                      u1::FloatVec, u2::FloatVec)
    return (u1, u2)
end

function resolve_args(rot::Rot90,
                      u1::FloatVec, u2::FloatVec)
    return (1-u2, u1)
end

function resolve_args(rot::Rot180,
                     u1::FloatVec, u2::FloatVec)
    return (1-u1, 1-u2)
end

function resolve_args(rot::Rot270,
                     u1::FloatVec, u2::FloatVec)
    return (u2, 1-u1)
end

##########################
## function evaluations ##
##########################

## pdf
##----

function pdf(cop::PairCop, rot::CWRotation,
             u1::FloatVec, u2::FloatVec)
    v1, v2 = resolve_args(rot, u1, u2)
    return pdf(cop, v1, v2)
end

## cdf
##----

function cdf(cop::PairCop, rot::Rot0,
             u1::FloatVec, u2::FloatVec)
    return cdf(cop, u1, u2)
end

function cdf(cop::PairCop, rot::Rot90,
             u1::FloatVec, u2::FloatVec)
    return u1 - cdf(cop, 1-u2, u1)
end

function cdf(cop::PairCop, rot::Rot180,
             u1::FloatVec, u2::FloatVec)
    return u1 + u2 + cdf(cop, 1-u1, 1-u2) - 1
end

function cdf(cop::PairCop, rot::Rot270,
             u1::FloatVec, u2::FloatVec)
    return u2 - cdf(cop, u2, 1-u1)
end

## h-functions
##------------

function hfun(cop::PairCop, rot::Rot0,
              u1::FloatVec, u2::FloatVec)
    return hfun(cop, u1, u2)
end

function hfun(cop::PairCop, rot::Rot90,
              u1::FloatVec, u2::FloatVec)
    return vfun(cop, u1, 1-u2)
end

function hfun(cop::PairCop, rot::Rot180,
              u1::FloatVec, u2::FloatVec)
    return 1 - hfun(cop, 1-u1, 1-u2)
end

function hfun(cop::PairCop, rot::Rot270,
              u1::FloatVec, u2::FloatVec)
    return 1 - vfun(cop, 1-u1, u2)
end

## v-functions
##------------

function vfun(cop::PairCop, rot::Rot0,
              u2::FloatVec, u1::FloatVec)
    return vfun(cop, u2, u1)
end

function vfun(cop::PairCop, rot::Rot90,
              u2::FloatVec, u1::FloatVec)
    return 1 - hfun(cop, 1-u2, u1)
end

function vfun(cop::PairCop, rot::Rot180,
              u2::FloatVec, u1::FloatVec)
    return 1 - vfun(cop, 1-u2, 1-u1)
end

function vfun(cop::PairCop, rot::Rot270,
              u2::FloatVec, u1::FloatVec)
    return hfun(cop, u2, 1-u1)
end

## inverse h-functions
##--------------------

function hinv(cop::PairCop, rot::Rot0,
              u1::FloatVec, u2::FloatVec)
    return hinv(cop, u1, u2)
end

function hinv(cop::PairCop, rot::Rot90,
              u1::FloatVec, u2::FloatVec)
    return vinv(cop, u1, 1-u2)
end

function hinv(cop::PairCop, rot::Rot180,
              u1::FloatVec, u2::FloatVec)
    return 1 - hinv(cop, 1-u1, 1-u2)
end

function hinv(cop::PairCop, rot::Rot270,
              u1::FloatVec, u2::FloatVec)
    return 1 - vinv(cop, 1-u1, u2)
end

## inverse v-functions
##--------------------

function vinv(cop::PairCop, rot::Rot0,
              u2::FloatVec, u1::FloatVec)
    return vinv(cop, u2, u1)
end

function vinv(cop::PairCop, rot::Rot90,
              u2::FloatVec, u1::FloatVec)
    return 1 - hinv(cop, 1-u2, u1)
end

function vinv(cop::PairCop, rot::Rot180,
              u2::FloatVec, u1::FloatVec)
    return 1 - vinv(cop, 1-u2, 1-u1)
end

function vinv(cop::PairCop, rot::Rot270,
              u2::FloatVec, u1::FloatVec)
    return hinv(cop, u2, 1-u1)
end


text:
using FVM_PDEsolver
using Documenter

DocMeta.setdocmeta!(FVM_PDEsolver, :DocTestSetup, :(using FVM_PDEsolver); recursive=true)

makedocs(;
    modules=[FVM_PDEsolver],
    authors="Amir Farzin",
    repo="https://github.com/amirfarzin/FVM_PDEsolver.jl/blob/{commit}{path}#{line}",
    sitename="FVM_PDEsolver.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)


text:
import Pkg
Pkg.instantiate() # downloads all dependencies for the current project

Pkg.add("Distributions")
Pkg.add("Random")

using StatsBase, StatsPlots, Distributions, Random

# Example data
data = (randn(10_000), randn(10_000))

# Plot StatsPlots 2D histogram
histogram2d(data)

# Fit a histogram with StatsBase
h = fit(Histogram, data)
x = searchsortedfirst(h.edges[1], 0.1)      # returns 10
y = searchsortedfirst(h.edges[2], 0.1)      # returns 11
h.weights[x, y]                             # returns 243

# Or as a function
function get_freq(h, xval, yval)
    x = searchsortedfirst(h.edges[1], xval)
    y = searchsortedfirst(h.edges[2], yval)
    h.weights[x, y]
end

get_freq(h, 1.4, 0.6)                       # returns 32


text:
include("../matplotlib_settings.jl")

using GeometricExamples

include("guiding_center_4d_settings_firk.jl")

tableau_list = get_tableau_list_firk()

equ = problem.guiding_center_4d_ode()

include("guiding_center_4d.jl")

run_tableau_list(tableau_list, equ, Δt, ntime)

generate_html(tableau_list, figures, prefix, fig_suff, fig_dir, out_dir, html_file;
              title=coll_title, params=get_params(Δt, ntime, problem))


text:
abstract type Task end

for pre in ("AI", "AO", "DI", "DO", "CI", "CO")
  @eval mutable struct $(Symbol(pre*"Task")) <: Task
      th::TaskHandle
  end
  @eval $(Symbol(pre*"Task"))() = $(Symbol(pre*"Task"))(task())
  @eval $(Symbol(pre*"Task"))(s::String) = $(Symbol(pre*"Task"))(task(s))
end

function task(name::String)
    t = TaskHandle[0]
    catch_error( DAQmxCreateTask(Ref(codeunits(name),1), pointer(t)) )
    t[1]
end
task() = task("")

for (cfunction, jfunction) in (
        (DAQmxStartTask, :start),
        (DAQmxStopTask,  :stop),
        (DAQmxClearTask, :clear))

    @eval function $jfunction(t::Task)
        catch_error( $cfunction(t.th) )
        nothing
    end

    @eval @doc $(string("`", jfunction, """(task)`

    """,jfunction," the specified NIDAQ task")) $jfunction
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module ConversionTests

using Test
using Dates

@testset "conversion to/from UNIX" begin
    # Test conversion to and from unix
    @test Dates.unix2datetime(Dates.datetime2unix(DateTime(2000, 1, 1))) == DateTime(2000, 1, 1)
    @test Dates.value(Dates.DateTime(1970)) == Dates.UNIXEPOCH

    # Tests from here: https://en.wikipedia.org/wiki/Unix_time
    @test string(Dates.unix2datetime(1095379198.75)) == string("2004-09-16T23:59:58.75")
    @test string(Dates.unix2datetime(1095379199.00)) == string("2004-09-16T23:59:59")
    @test string(Dates.unix2datetime(1095379199.25)) == string("2004-09-16T23:59:59.25")
    @test string(Dates.unix2datetime(1095379199.50)) == string("2004-09-16T23:59:59.5")
    @test string(Dates.unix2datetime(1095379199.75)) == string("2004-09-16T23:59:59.75")
    @test string(Dates.unix2datetime(1095379200.00)) == string("2004-09-17T00:00:00")
    @test string(Dates.unix2datetime(1095379200.25)) == string("2004-09-17T00:00:00.25")
    @test string(Dates.unix2datetime(1095379200.50)) == string("2004-09-17T00:00:00.5")
    @test string(Dates.unix2datetime(1095379200.75)) == string("2004-09-17T00:00:00.75")
    @test string(Dates.unix2datetime(1095379201.00)) == string("2004-09-17T00:00:01")
    @test string(Dates.unix2datetime(1095379201.25)) == string("2004-09-17T00:00:01.25")
    @test string(Dates.unix2datetime(915148798.75)) == string("1998-12-31T23:59:58.75")
    @test string(Dates.unix2datetime(915148799.00)) == string("1998-12-31T23:59:59")
    @test string(Dates.unix2datetime(915148799.25)) == string("1998-12-31T23:59:59.25")
    @test string(Dates.unix2datetime(915148799.50)) == string("1998-12-31T23:59:59.5")
    @test string(Dates.unix2datetime(915148799.75)) == string("1998-12-31T23:59:59.75")
    @test string(Dates.unix2datetime(915148800.00)) == string("1999-01-01T00:00:00")
    @test string(Dates.unix2datetime(915148800.25)) == string("1999-01-01T00:00:00.25")
    @test string(Dates.unix2datetime(915148800.50)) == string("1999-01-01T00:00:00.5")
    @test string(Dates.unix2datetime(915148800.75)) == string("1999-01-01T00:00:00.75")
    @test string(Dates.unix2datetime(915148801.00)) == string("1999-01-01T00:00:01")
    @test string(Dates.unix2datetime(915148801.25)) == string("1999-01-01T00:00:01.25")
end

@testset "conversion to/from Rata Die" begin
    @test Date(Dates.rata2datetime(734869)) == Dates.Date(2013, 1, 1)
    @test Dates.datetime2rata(Dates.rata2datetime(734869)) == 734869
end

@testset "conversion to/from Julian calendar" begin
    # Tests from here: http://mysite.verizon.net/aesir_research/date/back.htm#JDN
    @test Dates.julian2datetime(1721119.5) == Dates.DateTime(0, 3, 1)
    @test Dates.julian2datetime(1721424.5) == Dates.DateTime(0, 12, 31)
    @test Dates.julian2datetime(1721425.5) == Dates.DateTime(1, 1, 1)
    @test Dates.julian2datetime(2299149.5) == Dates.DateTime(1582, 10, 4)
    @test Dates.julian2datetime(2415020.5) == Dates.DateTime(1900, 1, 1)
    @test Dates.julian2datetime(2415385.5) == Dates.DateTime(1901, 1, 1)
    @test Dates.julian2datetime(2440587.5) == Dates.DateTime(1970, 1, 1)
    @test Dates.julian2datetime(2444239.5) == Dates.DateTime(1980, 1, 1)
    @test Dates.julian2datetime(2452695.625) == Dates.DateTime(2003, 2, 25, 3)
    @test Dates.datetime2julian(Dates.DateTime(2013, 12, 3, 21)) == 2456630.375
end
@testset "now, today, and UTC" begin
    @test typeof(Dates.now()) <: Dates.DateTime
    @test typeof(Dates.today()) <: Dates.Date
    @test typeof(Dates.now(Dates.UTC)) <: Dates.DateTime

    if Sys.isapple()
        withenv("TZ" => "UTC") do
            @test abs(Dates.now() - now(Dates.UTC)) < Dates.Second(1)
        end
    end
    @test abs(Dates.now() - now(Dates.UTC)) < Dates.Hour(16)
end
@testset "Issue #9171, #9169" begin
    let t = Dates.Period[Dates.Week(2), Dates.Day(14), Dates.Hour(14 * 24), Dates.Minute(14 * 24 * 60), Dates.Second(14 * 24 * 60 * 60), Dates.Millisecond(14 * 24 * 60 * 60 * 1000)]
        for i = 1:length(t)
            Pi = typeof(t[i])
            for j = 1:length(t)
                @test t[i] == t[j]
            end
            for j = i+1:length(t)
                Pj = typeof(t[j])
                tj1 = t[j] + Pj(1)
                @test t[i] < tj1
                @test_throws InexactError Pi(tj1)
                @test_throws InexactError Pj(Pi(typemax(Int64)))
                @test_throws InexactError Pj(Pi(typemin(Int64)))
            end
        end
    end
    @test Dates.Year(3) == Dates.Month(36)
    #@test_throws MethodError Int(Dates.Month(36)) # re-enable when deprecation removed
    @test Dates.Year(3) < Dates.Month(37)
    @test_throws InexactError convert(Dates.Year, Dates.Month(37))
    @test_throws InexactError Dates.Month(Dates.Year(typemax(Int64)))
end

@testset "Conversions to/from numbers" begin
    # Ensure that conversion of 32-bit integers work
    let dt = DateTime(1915, 1, 1, 12)
        unix = Int32(Dates.datetime2unix(dt))
        julian = Int32(Dates.datetime2julian(dt))

        @test Dates.unix2datetime(unix) == dt
        @test Dates.julian2datetime(julian) == dt
    end

    a = Dates.DateTime(2000)
    b = Dates.Date(2000)
    @test Dates.value(b) == 730120
    @test Dates.value(a) == 63082368000000
    @test convert(Dates.DateTime, Dates.Millisecond(63082368000000)) == a
    @test convert(Dates.Millisecond, a) == Dates.Millisecond(63082368000000)
    @test Dates.DateTime(Dates.UTM(63082368000000)) == a
    @test Dates.DateTime(Dates.UTM(63082368000000.0)) == a
    @test convert(Dates.Date, Dates.Day(730120)) == b
    @test convert(Dates.Day, b) == Dates.Day(730120)
    @test Dates.Date(Dates.UTD(730120)) == b
    @test Dates.Date(Dates.UTD(730120.0)) == b
    @test Dates.Date(Dates.UTD(Int32(730120))) == b

    dt = Dates.DateTime(2000, 1, 1, 23, 59, 59, 50)
    t = Dates.Time(dt)
    @test Dates.hour(t) == 23
    @test Dates.minute(t) == 59
    @test Dates.second(t) == 59
    @test Dates.millisecond(t) == 50
    @test Dates.microsecond(t) == 0
    @test Dates.nanosecond(t) == 0
end

end

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
# WARNING this file is only loaded if DataFrames.jl is also active
# This all happens thanks to the Requires.jl package

import Base: getindex
import Base: setindex!
import SimpleSDMLayers: clip, latitudes, longitudes

"""
    Base.getindex(layer::T, df::DataFrames.DataFrame; latitude = :latitude, longitude = :longitude) where {T <: SimpleSDMLayer}

Returns the values of a layer at all occurrences in a `DataFrame`.
"""
function Base.getindex(layer::T, df::DataFrames.DataFrame; latitude = :latitude, longitude = :longitude) where {T <: SimpleSDMLayer}
    lats = df[:, latitude]
    lons = df[:, longitude]
    return [layer[lon, lat] for (lon, lat) in zip(lons, lats)]
end

"""
    clip(layer::T, df::DataFrames.DataFrame; latitude = :latitude, longitude = :longitude) where {T <: SimpleSDMLayer}

Returns a clipped version (with a 10% margin) around all occurences in a
`DataFrame`.
"""
function SimpleSDMLayers.clip(layer::T, df::DataFrames.DataFrame; latitude = :latitude, longitude = :longitude) where {T <: SimpleSDMLayer}
   occ_latitudes = filter(!ismissing, df[:, latitude])
   occ_longitudes = filter(!ismissing, df[:, longitude])

   lat_min = minimum(occ_latitudes)
   lat_max = maximum(occ_latitudes)
   lon_min = minimum(occ_longitudes)
   lon_max = maximum(occ_longitudes)

   lat_Δ = abs(lat_max - lat_min)
   lon_Δ = abs(lon_max - lon_min)

   scaling = 0.1
   lon_s = scaling*lon_Δ
   lat_s = scaling*lat_Δ

   lat_max = min(layer.top, lat_max+lat_s)
   lat_min = max(layer.bottom, lat_min-lat_s)
   lon_max = min(layer.right, lon_max+lon_s)
   lon_min = max(layer.left, lon_min-lon_s)

   return layer[left=lon_min, right=lon_max, bottom=lat_min, top=lat_max]
end

"""
    DataFrames.DataFrame(layer::T) where {T <: SimpleSDMLayer}

Returns a DataFrame from a `SimpleSDMLayer` element, with columns for latitudes,
longitudes and grid values.
"""
function DataFrames.DataFrame(layer::T; kw...) where {T <: SimpleSDMLayer}
    lats = repeat(latitudes(layer), outer = size(layer, 2))
    lons = repeat(longitudes(layer), inner = size(layer, 1))
    return DataFrames.DataFrame(latitude = lats, longitude = lons, values = vec(layer.grid); kw...)
end

"""
    DataFrames.DataFrame(layers::Array{SimpleSDMLayer})

Returns a single DataFrame from an `Array` of compatible`SimpleSDMLayer`
elements, with every layer as a column, as well as columns for latitudes and longitudes.
"""
function DataFrames.DataFrame(layers::Array{T}; kw...) where {T <: SimpleSDMLayer}
    l1 = layers[1]
    all(x -> SimpleSDMLayers._layers_are_compatible(x, l1), layers)

    lats = repeat(latitudes(l1), outer = size(l1, 2))
    lons = repeat(longitudes(l1), inner = size(l1, 1))
    values = mapreduce(x -> vec(x.grid), hcat, layers)

    df = DataFrames.DataFrame(values; kw...)
    DataFrames.insertcols!(df, 1, :latitude => lats)
    DataFrames.insertcols!(df, 1, :longitude => lons)
    return df
end

for ty in (:SimpleSDMResponse, :SimpleSDMPredictor)
    eval(
        quote
            """
                $($ty)(df::DataFrame, col::Symbol, layer::T; latitude::Symbol = :latitude, longitude::Symbol = :longitude) where {T <: SimpleSDMLayer}

            Returns a `$($ty)` from a `DataFrame`.
            """
            function SimpleSDMLayers.$ty(df::DataFrames.DataFrame, col::Symbol, layer::SimpleSDMLayer; latitude::Symbol = :latitude, longitude::Symbol = :longitude)
                lats = df[:, latitude]
                lons = df[:, longitude]

                uniquelats = unique(lats)
                uniquelons = unique(lons)

                grid = Array{Any}(nothing, size(layer))

                if uniquelats == latitudes(layer) && uniquelons == longitudes(layer)
                    grid[:] = df[:, col]
                else
                    lats_idx = [SimpleSDMLayers._match_latitude(layer, lat) for lat in lats]
                    lons_idx = [SimpleSDMLayers._match_longitude(layer, lon) for lon in lons]
                    for (lat, lon, value) in zip(lats_idx, lons_idx, df[:, col])
                        grid[lat, lon] = value
                    end
                end

                internal_types = unique(typeof.(grid))
                return SimpleSDMLayers.$ty(Array{Union{internal_types...}}(grid), layer)
            end
        end,
    )
end


text:
function transpose_strings(input::AbstractArray)
    # This saves us from having to make a special case for maximum below
    input == [] && return []

    # Construct a matrix from the strings where empty spaces are replaced with missing
    zenith = maximum(length,input)
    result = Array{Union{Missing,Char},2}(undef,length(input),zenith)
    for (i,s) in enumerate(input)
        for j in 1:zenith
            @inbounds result[i,j] = get(s,j,missing)
        end
    end

    # For every column, remove trailing missings, replace the rest with spaces and convert into a string
    [string_lpadmissing!(r) for r in eachcol(result)]
end

# This modifies its input and returns something else which may be surprsing
# but as we're using eachcol we get a copy anyway
function string_lpadmissing!(maybechars)
    # replace the missings up to the last character with spaces
    last = findlast(!ismissing,maybechars)
    last = isnothing(last) ? 0 : last
    for i in 1:last-1 # we know last isn't missing
        if ismissing(maybechars[i]); maybechars[i] = ' '; end
    end

    # leave out the missings after the last character
    String(Array{Char,1}(maybechars[1:last]))
end

# For testing:
strs = [rand("qwertyuiop",rand(900:1000,1)|>first)|>join for _ in 1:1_000];

# The tests could be harder: they don't test for preserving trailing spaces in the input, e.g.
# ```
# julia> transpose_strings(["abcedf","abc","def "])
# 6-element Array{String,1}:
#  "aad"
#  "bbe"
#  "ccf"
#  "e  "
#  "d"
#  "f"
# ```
#
# From the task description:
# In general, all characters from the input should also be present in the transposed output. That means that if a column in the input text contains only spaces on its bottom-most row(s), the corresponding output row should contain the spaces in its right-most column(s).
#
# (without this constraint a ~two liner with hcat, lpad, and rstrip suffices)


text:
using Crayons: Crayon
using Match: @match
using PyCall: PyError

INFO_STYLE = Crayon(foreground = :blue)
INPUT_STYLE  = Crayon(foreground = :light_yellow)
FILE_INPUT_STYLE = Crayon(foreground = :magenta)
OUTPUT_STYLE = Crayon(foreground = :white, bold=false)
ERROR_STYLE  = Crayon(foreground = :red)

function help()
    println(INFO_STYLE, "Input uses EQN format.")
    println(OUTPUT_STYLE, """
        \t▪ Use curly brackets
        \t▪ Use 'over' for division and 'sup' for power""")
    println(INFO_STYLE, "Math functions:")
    println(OUTPUT_STYLE, """
        \tsqrt   log     log2    log10   exp
        \tsin    cos    tan    cot    sec    csc
        \tsinh   cosh   tanh   coth   sech   csch
        \tasin   acos   atan   acot   asec   acsc
        \tasinh  acosh  atanh  acoth  asech  acsch""")
    println(INFO_STYLE, "Other functions:")
    println(OUTPUT_STYLE, """
        \t▪ compute(x): evaluate the symbol x or an expression.
        \t▪ compute!(x): same as above, but also sets the value of x.
        \t▪ file(filename): execute commands from specified file.""")
    println(INFO_STYLE, "\nCall exit() to close REPL.")
end

function process_line(line)
    nline = strip(line)
    if nline == "?"
        help()
    else
        ex = parseEqn(nline)
        repl_dispatch(ex)
    end
end

function file(filename)
    lines = filter(!isempty, [strip(ln) for ln in readlines(filename)])
    for line ∈ lines
        println(FILE_INPUT_STYLE, "> $(line)")
        process_line(line)
    end
end

function repl()
    print(INPUT_STYLE, "> ")
    try
        input = readline()
        process_line(input)
    catch e
        println(ERROR_STYLE,
            if isa(e, Meta.ParseError) || isa(e, PyError)
                "Parsing error: $(e.msg)"
            elseif isa(e, MethodError)
                "Signature error in function '$(e.f)'\nPassed args: " * join(["$(a)::$(typeof(a))" for a in e.args], ", ")
            elseif isa(e, SystemError)
                "System error $(e.errnum): $(e.prefix)"
            elseif isa(e, DomainError)
                "Domain error with $(e.val):\n$(e.msg)"
            elseif isa(e, ErrorException)
                "Error: $(e.msg)"
            elseif isa(e, StackOverflowError)
                "Error: stack overflow"
            else
                rethrow()
            end)
    end
end

repl_dispatch(::Nothing) = nothing
repl_dispatch(x) = println(OUTPUT_STYLE, x)

function repl_dispatch(x::Symbol)
    value = Context.get(x)
    println(OUTPUT_STYLE, string(value))
end

function repl_dispatch(ex::Assignment)
    l, r = ex.lvalue, ex.rvalue
    if isa(l, Expr)
        Context.eval(:($l = $r))
        println(OUTPUT_STYLE, "$l = $r")
    else
        Context.set(l, r)
        val = Context.get(l)
        println(OUTPUT_STYLE, "$l = $(string(val))")
    end
end

function repl_dispatch(int::Integral)
    result = Context.eval(int)
    println(OUTPUT_STYLE, string(result))
end

function repl_dispatch(ex::Expr)
    @assert ex.head == :call

    result = Context.eval(Context.eagerCalc(ex))
    if result ≠ nothing
        println(OUTPUT_STYLE, string(result))
    end
end


text:
using GeoStats
using StaticArrays

@testset "Interpolating rasters" begin
    @testset "Regular Grid interpolation" begin
        ga = GeoArray(Array{Union{Missing, Float64}}(rand(10, 10, 1)), GeoArrays.geotransform_to_affine(SVector(0.,1.,0.,0.,0.,1.)), "")
        ga.A[2,2,1] = missing
        @test count(ismissing, ga) == 1
        GeoArrays.interpolate!(ga, Kriging())
        @test count(ismissing, ga) == 0
    end
    @testset "Regular Grid interpolation with negative spacing" begin
        ga = GeoArray(Array{Union{Missing, Float64}}(rand(10, 10, 1)), GeoArrays.geotransform_to_affine(SVector(0.,1.,0.,0.,0.,-1.)), "")
        ga.A[2,2,1] = missing
        @test count(ismissing, ga) == 1
        GeoArrays.interpolate!(ga, Kriging())
        @test count(ismissing, ga) == 0
    end
    @testset "Irregular Grid interpolation (not yet supported)" begin
        ga = GeoArray(Array{Union{Missing, Float64}}(rand(10, 10, 1)), GeoArrays.geotransform_to_affine(SVector(0.,1.,1.,0.,0.,1.)), "")
        @test_throws ErrorException GeoArrays.interpolate!(ga, Kriging())
    end
end


text:
function border(x,y,rib)
    rib1, rib2 = if Plots.istuple(rib)
        first(rib), last(rib)
    else
        rib, rib
    end
    yline = vcat(y-rib1,(y+rib2)[end:-1:1])
    xline = vcat(x,x[end:-1:1])
    return xline, yline
end

@recipe function f(::Type{Val{:shadederror}},plt::Plot; shade = 0.)

    # set up the subplots
    x,y = d[:x],d[:y]
    xline, yline = border(x,y,shade)

    # line plot
    @series begin
        primary := true
        x := x
        y := y
        seriestype := :path
        ()
    end

    @series begin
        # shaded error bar
        primary := false
        fillrange := 0
        fillalpha --> 0.5
        linewidth := 0
        x := xline
        y := yline
        seriestype := :path
        ()
    end
end

@shorthands shadederror


text:
using Catalyst, ModelingToolkit, OrdinaryDiffEq, Test, LinearAlgebra, DiffEqJump

@parameters k α
@variables t, A(t), B(t), C(t), D(t)
rxs = [Reaction(t*k, [A], [B], [2*α^2], [k+α*C])
       Reaction(1.0, [A,B], [C,D], [α,2], [k,α])
]
@named rs = ReactionSystem(rxs, t)
@test issetequal(states(rs), [A,B,C,D])
@test issetequal(parameters(rs), [k,α])
osys = convert(ODESystem, rs)

g = (k+α*C)
rs2 = @reaction_network rs begin
    t*k, 2*α^2*A --> $g*B
    1.0, α*A + 2*B --> k*C + α*D
end k α
@test rs2 == rs

rxs2 = [(@reaction t*k, 2*α^2*A --> $g*B),
        (@reaction 1.0, α*A + 2*B --> k*C + α*D)]
rs3  = ReactionSystem(rxs2, t; name=:rs)
@test rs3 == rs

u0map = [A => 3.0, B => 2.0, C => 3.0, D => 1.5]
pmap  = (k => 2.5, α => 2)
tspan = (0.0,5.0)
oprob = ODEProblem(osys, u0map, tspan, pmap)
# this is a hack because of https://github.com/SciML/ModelingToolkit.jl/issues/1475
oprob = remake(oprob, p=Tuple(pv[2] for pv in pmap))
du1 = zeros(size(oprob.u0))
oprob.f(du1,oprob.u0,oprob.p,1.5)

function oderhs(du,u,p,t)
    k = p[1]; α = p[2]
    A = u[1]; B = u[2]; C = u[3]; D = u[4]
    n = 2*α^2
    rl = t*k/factorial(n) * A^n
    rl2 = A^α*B^2 / (2 * factorial(α))
    du[1] = -n*rl - α*rl2
    du[2] = (k+α*C)*rl - 2*rl2
    du[3] = k*rl2
    du[4] = α*rl2
end
u0 = [uv[2] for uv in u0map]
p  = Tuple(pv[2] for pv in pmap)
oprob2 = ODEProblem(oderhs, u0, tspan, p)
du2 = copy(du1)
oprob2.f(du2,oprob2.u0,oprob2.p,1.5)
@test norm(du1 .- du2) < 100*eps()

# test without rate law scalings
osys = convert(ODESystem, rs, combinatoric_ratelaws=false)
oprob = ODEProblem(osys, u0map, tspan, pmap)
function oderhs(du,u,p,t)
    k = p[1]; α = p[2]
    A = u[1]; B = u[2]; C = u[3]; D = u[4]
    n = 2*α^2
    rl = t*k * A^n
    rl2 = A^α*B^2
    du[1] = -n*rl - α*rl2
    du[2] = (k+α*C)*rl - 2*rl2
    du[3] = k*rl2
    du[4] = α*rl2
end
oprob2 = ODEProblem(oderhs, [uv[2] for uv in u0map], tspan, oprob.p)
du1 .= 0; du2 .=0
oprob.f(du1,oprob.u0,oprob.p,1.5)
oprob2.f(du2,oprob2.u0,oprob2.p,1.5)
@test norm(du1 .- du2) < 100*eps()

# SDESystem test
ssys = convert(SDESystem, rs)
sf = SDEFunction{false}(ssys, states(ssys), parameters(ssys))
G  = sf.g(u0,p,1.0)
function sdenoise(u,p,t)
    k = p[1]; α = p[2]
    A = u[1]; B = u[2]; C = u[3]; D = u[4]
    n = 2*α^2
    rl = sqrt(t*k/factorial(n) * A^n)
    rl2 = sqrt(A^α*B^2 / (2 * factorial(α)))
    G = [-n*rl (-α*rl2);
         (k+α*C)*rl (-2*rl2);
         0.0 k*rl2;
         0.0 α*rl2]
end
G2 = sdenoise(u0,p,1.0)
@test norm(G-G2) < 100*eps()

# SDESystem test with no combinatoric rate laws
ssys = convert(SDESystem, rs, combinatoric_ratelaws=false)
sf = SDEFunction{false}(ssys, states(ssys), parameters(ssys))
G  = sf.g(u0,p,1.0)
function sdenoise(u,p,t)
    k = p[1]; α = p[2]
    A = u[1]; B = u[2]; C = u[3]; D = u[4]
    n = 2*α^2
    rl = sqrt(t*k*A^n)
    rl2 = sqrt(A^α*B^2)
    G = [-n*rl (-α*rl2);
         (k+α*C)*rl (-2*rl2);
         0.0 k*rl2;
         0.0 α*rl2]
end
G2 = sdenoise(u0,p,1.0)
@test norm(G-G2) < 100*eps()

# JumpSystem test
js  = convert(JumpSystem, rs)
u0map = [A => 3, B => 2, C => 3, D => 5]
u0    = [uv[2] for uv in u0map]
pmap  = (k => 5, α => 2)
p     = Tuple(pv[2] for pv in pmap)
statetoid = Dict(state => i for (i,state) in enumerate(states(js)))
function r1(u,p,t)
    α = p[2]; A = u[1]
    t*k * binomial(A, 2*α^2)
end
function affect1!(integrator)
    A = integrator.u[1]
    B = integrator.u[2]
    C = integrator.u[3]
    k = p[1]; α = p[2]
    integrator.u[1] -= 2*α^2
    integrator.u[2] += (k+α*C)
    nothing
end
ttt = 1.5
j1 = VariableRateJump(r1,affect1!)
vrj = ModelingToolkit.assemble_vrj(js, equations(js)[2], statetoid)
@test isapprox(vrj.rate(u0,p,ttt), r1(u0,p,ttt))
fake_integrator1 = (u=copy(u0),p=p,t=ttt); fake_integrator2 = deepcopy(fake_integrator1);
vrj.affect!(fake_integrator1);
affect1!(fake_integrator2);
@test fake_integrator1 == fake_integrator2

function r2(u,p,t)
    α = p[2]; A = u[1]; B = u[2]
    binomial(A,α) * B * (B-1) / 2
end
function affect2!(integrator)
    A = integrator.u[1]
    B = integrator.u[2]
    C = integrator.u[3]
    D = integrator.u[4]
    k = p[1]; α = p[2]
    integrator.u[1] -= α
    integrator.u[2] -= 2
    integrator.u[3] += k
    integrator.u[4] += α
    nothing
end
j2 = ConstantRateJump(r2,affect2!)
crj = ModelingToolkit.assemble_crj(js, equations(js)[1], statetoid)
@test isapprox(crj.rate(u0,p,ttt), r2(u0,p,ttt))
fake_integrator1 = (u=copy(u0),p=p,t=ttt); fake_integrator2 = deepcopy(fake_integrator1);
crj.affect!(fake_integrator1);
affect2!(fake_integrator2);
@test fake_integrator1 == fake_integrator2

# a few simple solving tests via the SIR Model
@parameters α β γ k
@variables t, S(t), I(t), R(t)
rxs = [Reaction(α, [S,I], [I], [1,1], [2]),
       Reaction(β, [I], [R], [1], [1])]
@named sir_ref = ReactionSystem(rxs, t)

rxs2 = [Reaction(α, [S,I], [I], [γ,1], [k]),
        Reaction(β, [I], [R], [γ], [γ])]
@named sir = ReactionSystem(rxs2, t)

@test issetequal(states(sir_ref), states(sir))

# ODEs
p1    = (α => .1/1000, β => .01)
p2    = (α => .1/1000, β => .01, γ => 1, k => 2)
tspan = (0.0,250.0)
u0    = [S => 999.0, I => 1.0, R => 0.0]
oprob = ODEProblem(sir_ref, u0, tspan, p1)
sol   = solve(oprob, Tsit5())
# here we hack around https://github.com/SciML/ModelingToolkit.jl/issues/1475
p2dict = Dict(p2)
pvs    = Tuple(p2dict[s] for s in parameters(sir))
@test all(isequal.(parameters(sir), parameters(convert(ODESystem, sir))))
oprob2 = ODEProblem(sir, u0, tspan, pvs)
sol2   = solve(oprob2, Tsit5())
@test norm(sol - sol2(sol.t)) < 1e-10

# jumps
Nsims = 10000
u0 = [S => 999, I => 1, R => 0]
jsys  = convert(JumpSystem, sir_ref)
dprob = DiscreteProblem(jsys, u0, tspan, p1)
jprob = JumpProblem(jsys, dprob, Direct(), save_positions=(false,false))
function getmean(jprob, tf)
    m = zeros(3)
    for i=1:Nsims
        sol = solve(jprob, SSAStepper())
        m .+= sol(tspan[2])
    end
    m /= Nsims
    m
end
m1 = getmean(jprob, tspan[2])
jsys2 = convert(JumpSystem, sir)
p2dict = Dict(p2)
pvs    = Tuple(p2dict[s] for s in parameters(sir))
@test all(isequal.(parameters(sir), parameters(jsys2)))
dprob2 = DiscreteProblem(jsys2, u0, tspan, pvs)
jprob2 = JumpProblem(jsys2, dprob2, Direct(), save_positions=(false,false))
m2 = getmean(jprob2, tspan[2])
@test maximum(abs.(m1[2:3] .- m2[2:3]) ./ m1[2:3]) < .05


text:
using WassersteinDictionaries
using Test

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

max_stars_repo_name: JuliaRCM/WassersteinDictionaries

max_issues_repo_name: JuliaRCM/WassersteinDictionaries

max_forks_repo_name: JuliaRCM/WassersteinDictionaries


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "ADIOS2"
version = v"2.7.1"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/ornladios/ADIOS2/archive/refs/tags/v2.7.1.tar.gz", "c8e237fd51f49d8a62a0660db12b72ea5067512aa7970f3fcf80b70e3f87ca3e"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd ADIOS2-2.7.1
# PR <https://github.com/ornladios/ADIOS2/pull/2712>
atomic_patch -p1 ${WORKSPACE}/srcdir/patches/ndims.patch
atomic_patch -p1 ${WORKSPACE}/srcdir/patches/shlwapi.patch
mkdir build
cd build
if [[ "$target" == *-apple-* ]]; then
    # Set up a wrapper script for the assembler. GCC's assembly
    # output isn't accepted by the LLVM assembler.
    as=$(which as)
    mv "$as" "$as.old"
    export AS="$as.old"
    ln -s "$WORKSPACE/srcdir/scripts/as.llvm" "$as"
fi
mpiopts=
winopts=
if [[ "$target" == x86_64-w64-mingw32 ]]; then
    # cmake's auto-detection doesn't work on Windows.
    # The SST and Table ADIOS2 components don't build on Windows
    # (reported in <https://github.com/ornladios/ADIOS2/issues/2705>)
    mpiopts="-DMPI_HOME=$prefix -DMPI_GUESS_LIBRARY_NAME=MSMPI -DMPI_C_LIBRARIES=msmpi64 -DMPI_CXX_LIBRARIES=msmpi64"
    winopts="-DADIOS2_USE_SST=OFF -DADIOS2_USE_Table=OFF"
elif [[ "$target" == *-mingw* ]]; then
    mpiopts="-DMPI_HOME=$prefix -DMPI_GUESS_LIBRARY_NAME=MSMPI"
    winopts="-DADIOS2_USE_SST=OFF -DADIOS2_USE_Table=OFF"
fi
# We disable Fortran because this would require building many more
# versions of the library
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN%.*}_gcc.cmake -DCMAKE_BUILD_TYPE=Release -DADIOS2_USE_Fortran=OFF -DADIOS2_BUILD_EXAMPLES=OFF -DBUILD_TESTING=OFF ${mpiopts} ${winopts} ..
make -j${nproc}
make -j${nproc} install
install_license ../Copyright.txt ../LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("aarch64", "linux"; libc="glibc"),
    Platform("aarch64", "linux"; libc="musl"),
    Platform("powerpc64le", "linux"; libc="glibc"),
    Platform("x86_64", "linux"; libc="glibc"),
    Platform("x86_64", "linux"; libc="musl"),
    Platform("x86_64", "macos"),
    Platform("x86_64", "windows"),

    # These platforms fail:

    # [22:03:24] /workspace/srcdir/ADIOS2-2.7.1/source/adios2/engine/ssc/SscReader.cpp:420:71: error: narrowing conversion of ‘18446744073709551613ull’ from ‘long long unsigned int’ to ‘unsigned int’ inside { } [-Wnarrowing]
    # [22:03:24]                  m_IO.DefineVariable<T>(b.name, {adios2::LocalValueDim});       \
    # (32-bit architectures are officially not supported.)
    #FAIL Platform("armv7l", "linux"; libc="glibc"),
    #FAIL Platform("i686", "linux"; libc="glibc"),
    #TODO Platform("i686", "linux"; libc="musl"),
    #TODO Platform("i686", "windows"),

    # [22:40:43] /workspace/srcdir/ADIOS2-2.7.1/source/adios2/toolkit/profiling/taustubs/tautimer.cpp:208:31: error: ‘__NR_gettid’ was not declared in this scope
    # [22:40:43]      mytid = (uint64_t)syscall(__NR_gettid);
    # (Likely the respective syscall does not exist on FreeBSD;
    # reported as <https://github.com/ornladios/ADIOS2/issues/2705>.)
    #FAIL Platform("x86_64", "freebsd"),
]
# Apparently, macOS doesn't use different C++ string APIs
platforms = expand_cxxstring_abis(platforms; skip=Sys.isapple)

# The products that we will ensure are always built
products = [
    LibraryProduct("libadios2_c", :libadios2_c),
    LibraryProduct("libadios2_c_mpi", :libadios2_c_mpi),
    LibraryProduct("libadios2_core", :libadios2_core),
    LibraryProduct("libadios2_core_mpi", :libadios2_core_mpi),
    LibraryProduct("libadios2_cxx11", :libadios2_cxx11),
    LibraryProduct("libadios2_cxx11_mpi", :libadios2_cxx11_mpi),
    LibraryProduct("libadios2_taustubs", :libadios2_taustubs),

    # Missing on Windows:
    # LibraryProduct("libadios2_atl", :libadios2_atl),
    # LibraryProduct("libadios2_dill", :libadios2_dill),
    # LibraryProduct("libadios2_evpath", :libadios2_evpath),
    # LibraryProduct("libadios2_ffs", :libadios2_ffs),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="Blosc_jll")),
    # We don't want to use Bzip2 because this would lock us into Julia ≥1.6
    # Dependency(PackageSpec(name="Bzip2_jll")),
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae")),
    Dependency(PackageSpec(name="MPICH_jll")),
    Dependency(PackageSpec(name="MicrosoftMPI_jll")),
    Dependency(PackageSpec(name="ZeroMQ_jll")),
    Dependency(PackageSpec(name="libpng_jll")),
    Dependency(PackageSpec(name="zfp_jll")),
    # We cannot use HDF5 because we need a HDF5 configuration with MPI support
    # Dependency(PackageSpec(name="HDF5_jll")),
]

# Build the tarballs, and possibly a `build.jl` as well.
# GCC 4 is too old for Windows; it doesn't have <regex.h>
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"5")

max_stars_repo_name: Roger-luo/Yggdrasil

max_issues_repo_name: Roger-luo/Yggdrasil

max_forks_repo_name: Roger-luo/Yggdrasil