text:
################################################################################
# Bfield module
################################################################################
module Bfield
using Unitful
using UnitfulAtomic
using UnitfulAstro
using LinearAlgebra
using Measurements
using IterTools
using Measurements
using LinearAlgebra

function constant(x::Array{T};
    direction::Array = [0.0, 0.0, 6.0e-10]) where {T <: Number}
    out = vec(direction)
end #constx

"""
    dipole(x::Number,y::Number,z::Number; m = [0.,0.,1.1e59])

Returns dipole field given a position and magnetic moment.
Assumes units are in SI.
"""
function dipole(x::Array{T}; m::Array = [0.,0.,1.1e59]) where {T <: Number}

    # println("pos = ", pos)
    rMag = norm(x)
    # println("pos = ", rMag)
    if rMag == 0
        rMag = 1.0e-6
    end
    rHat = x/rMag
    # println("rHat = ", rHat)
    μ0 = 1u"μ0"
    μ0 = μ0.val
    num = μ0*(3.0*rHat*(m⋅rHat)-m)
    den = 4*π*rMag^3
    out = num/den

end #dipole

g = 1e-4u"T"
μg = 1e-6*g

janssonDict = Dict(
    :bi    =>[measurement(0.1  , 1.8 )*μg,
              measurement(3.0  , 0.6 )*μg,
              measurement(-0.9 , 0.8 )*μg,
              measurement(-0.8 , 0.3 )*μg,
              measurement(-2.0 , 0.1 )*μg,
              measurement(-4.2 , 0.5 )*μg,
              measurement(0.   , 1.8 )*μg,
              measurement(2.7  , 1.8 )*μg
             ],
    :bring => measurement(0.1  , 0.1 )*μg,
    :hdisk => measurement(0.4  , 0.03)*u"kpc",
    :wdisk => measurement(0.27 , 0.08)*u"kpc",
    :bn    => measurement(1.4  , 0.1 )*μg,
    :bs    => measurement(-1.1 , 0.1 )*μg,
    :rn    => measurement(9.22 , 0.08)*u"kpc",
    :rs    => measurement(16.7 , 0.0 )*u"kpc",
    :wh    => measurement(0.2  , .12 )*u"kpc",
    :z0    => measurement(5.3  , 1.6 )*u"kpc",
    :bx    => measurement(4.6  , 0.3 )*μg,
    :Θ0x   => measurement(49.  , 1.  )*u"°",
    :rcx   => measurement(4.8  , 0.2 )*u"kpc",
    :rx    => measurement(2.9  , 0.1 )*u"kpc",
    :γ     => measurement(2.92 , 0.14),
    :rnx   => [5.1,6.3,7.1,8.3,9.8,11.4,12.7,15.5]*u"kpc",
    :rinner  => 1.0u"kpc",
    :rdisklims => [3.0,5.0,20.0]*u"kpc", #[molecular ring, start arms, stop arms]
    :θarms => 11.5u"°"
    )

export janssonDict

"""
    spiral(ϕ,rnx,b)

Returns ``r(ϕ) = rnx*exp(b*ϕ) ``
"""
spiral(ϕ,rnx,b) = rnx*exp(b*ϕ)
export spiral

"""
    galspiral(ϕ, params::Dict)

# Returns shifted logarithmic spiral for jansson2012
"""
function galspiral(ϕ, rnx, θarms = janssonDict[:θarms])
    b = tan(θarms)
    return spiral(ϕ-π,rnx,b)
end
export galspiral

logistic(z,h,w) = 1.0/(1.0 + exp(-2.0*(uconvert(Unitful.NoUnits,(abs(z)-h)/w))))
export logistic

function jansson2012(x::Array{T}; params::Dict = janssonDict) where {T <: Quantity}
    @assert typeof(x[1]) <: Unitful.Length "x-coordinate requires units of distance"
    # define μgauss in this way so it plays nicely with various notebooks (Pluto, Jupyter...)
    g = 1e-4u"T"
    μg = 1e-6*g

    # initialize field output as zero
    out = zeros(3)*μg

    # ensure that coordinates are in kpc
    x = map(d->uconvert(u"kpc",d),x)
    ρ = norm(x)

    # return zero field if inside the unfitted center sphere
    if ρ < params[:rinner]
        return out
    end

    # calculate cylindrical coordinates
    r = sqrt(x[1]^2+x[2]^2)

    # return zero field if outside the galaxy
    if r > params[:rdisklims][end]
        return out
    end

    # define coordinate array with proper units store cylindrical coordinates
    ϕ = x[1].val == 0.0 ? 0.0u"rad" : atan(x[2],x[1])*u"rad"
    z = x[3]

    # generate cylindrical unit vectors
    rhat = [cos(ϕ),sin(ϕ),0.0]
    ϕhat = [-sin(ϕ),cos(ϕ),0.0]
    zhat = [0.0,0.0,1.0]

    # DISK COMPONENT
    if r > params[:rdisklims][2] # arm region
        # determine spiral arm in which point resides

        #   TODO: spiral search could be faster generate spiral transition
        #   radii at given ϕ
        """
            Returns index of spiral in which point resides
        """
        function whichspiral(r,ϕ,rnx = params[:rnx],θarms = params[:θarms])
            for ϕ0 in [-2*π,0.0,2*π,4*π]
                for (i,r0) in enumerate(rnx)
                    rSpiral = galspiral(ϕ+ϕ0,r0,θarms)
                    if r < rSpiral
                        return i
                    end
                end
            end
            @assert true "There is a bug in whichspiral"
        end

        # used to get field magnitude later
        iSpiral = whichspiral(r,ϕ)
        # disk field direction
        bhat = sin(params[:θarms])*rhat + cos(params[:θarms])*ϕhat
        # the spiral arm disk field scales as 1/r
        bdisk = Measurements.value(params[:bi][iSpiral])*uconvert(Unitful.NoUnits,params[:rdisklims][2]/r)*bhat
        # use logistic function to smoothly transition to halo field
        bdisk *= 1.0 - logistic(z,Measurements.value(params[:hdisk]),Measurements.value(params[:wdisk]))
        # superpose disk field onto total field
        out += bdisk

    elseif r > params[:rdisklims][1] # molecular ring
        # molecular ring is purely toroidal
        bdisk = Measurements.value(params[:bring])*ϕhat
        # use logistic function transition
        bdisk *= 1.0 - logistic(z,Measurements.value(params[:hdisk]),Measurements.value(params[:wdisk]))
        # superpose disk field onto total field
        out += bdisk
    end

    # TOROIDAL HALO COMPONENT
    bhmag = exp(-1*uconvert(Unitful.NoUnits,abs(z)/Measurements.value(params[:z0])))
    bhmag *= logistic(z,Measurements.value(params[:hdisk]),Measurements.value(params[:wdisk]))
    if z.val >= 0.0
        bhmag *= Measurements.value(params[:bn])*(1.0
            - logistic(z,Measurements.value(params[:rn]),Measurements.value(params[:wh])))
    else
        bhmag *= Measurements.value(params[:bs])*(1.0
            - logistic(z,Measurements.value(params[:rs]),Measurements.value(params[:wh])))
    end
    out += bhmag*ϕhat

    # X FIELD COMPONENT

    # get values from parameters for easier reading
    # TODO: substitute these back in after some readability adjustments
    Θ0x = Measurements.value(params[:Θ0x])
    rx = Measurements.value(params[:rx])
    bx0 = Measurements.value(params[:bx])
    rcx = Measurements.value(params[:rcx])

    rp = r - abs(z)/tan(Θ0x)
    rp = rp.val < 0.0 ? 0.0*u"kpc" : rp
    if rp > rcx
        bx = bx0*exp(uconvert(Unitful.NoUnits,-1*rp/rx))
        bmag = bx*uconvert(Unitful.NoUnits,rp/r)
        bhat = z.val > 0.0 ? [cos(ϕ)*cos(Θ0x),sin(ϕ)*cos(Θ0x),sin(Θ0x)] : [cos(ϕ)*cos(Θ0x),sin(ϕ)*cos(Θ0x),-1*sin(Θ0x)]
        out += bmag*bhat
    else
        rp = r*rcx/(rcx+abs(z)/tan(Θ0x))
        bx = bx0*exp(uconvert(Unitful.NoUnits,-1*rp/rx))
        bmag = bx*uconvert(Unitful.NoUnits,rp/r)^2
        Θx = r == rp ? 0.0u"rad" : atan(abs(z)/(r-rp))
        bhat = z.val > 0.0 ? [cos(ϕ)*cos(Θx),sin(ϕ)*cos(Θx),sin(Θx)] : [cos(ϕ)*cos(Θx),sin(ϕ)*cos(Θx),-1*sin(Θx)]
        out += bmag*bhat
    end

    return out
end #jansson2012

export jansson2012
########################### B-field wrappers ##################################
function galacticDipole(x::Array{T}) where {T <: Number}
    dipole(x)
end #galacticDipole

function galacticConstantZ(x::Array{T}) where {T <: Number}
    constant(x)
end #galacticConstantX

function jansson2012m(x::Array{T}) where {T <: Number}
    xin = uconvert.(u"kpc",x*u"m")
    return ustrip.(u"T",jansson2012(xin))
end

function jansson2012kpc(x::Array{T}) where {T <: Number}
    xin = x*u"kpc"
    return ustrip.(u"T",jansson2012(xin))
end
######################## end B-field wrappers ##################################

"""
`Bfield.model`: dictionary containing commonly used magnetic fields.
"""
model = Dict(
    :galacticDipole         => galacticDipole,
    :galacticConstantZ      => galacticConstantZ,
    :jansson2012m           => jansson2012m,
    :jansson2012kpc         => jansson2012kpc
)

end #Bfield

export Bfield


text:
# ---
# title: 1410. HTML Entity Parser
# id: problem1410
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: String, Stack
# link: <https://leetcode.com/problems/html-entity-parser/description/>
# hidden: true
# ---
#
# **HTML entity parser** is the parser that takes HTML code as input and replace
# all the entities of the special characters by the characters itself.
#
# The special characters and their entities for HTML are:
#
#   * **Quotation Mark:**  the entity is `"` and symbol character is `"`.
#   * **Single Quote  Mark:** the entity is `'` and symbol character is `'`.
#   * **Ampersand:**  the entity is `&` and symbol character is `&`.
#   * **Greater Than Sign:**  the entity is `>` and symbol character is `>`.
#   * **Less Than Sign:**  the entity is `<` and symbol character is `<`.
#   * **Slash:**  the entity is `⁄` and symbol character is `/`.
#
# Given the input `text` string to the HTML parser, you have to implement the
# entity parser.
#
# Return _the text_ after replacing the entities by the special characters.
#
#
#
# **Example 1:**
#
#
#
#     Input: text = "& is an HTML entity but &ambassador; is not."
#     Output: "& is an HTML entity but &ambassador; is not."
#     Explanation: The parser will replace the & entity by &
#
#
# **Example 2:**
#
#
#
#     Input: text = "and I quote: "...""
#     Output: "and I quote: \"...\""
#
#
# **Example 3:**
#
#
#
#     Input: text = "Stay home! Practice on Leetcode :)"
#     Output: "Stay home! Practice on Leetcode :)"
#
#
# **Example 4:**
#
#
#
#     Input: text = "x > y && x < y is always false"
#     Output: "x > y && x < y is always false"
#
#
# **Example 5:**
#
#
#
#     Input: text = "leetcode.com⁄problemset⁄all"
#     Output: "leetcode.com/problemset/all"
#
#
#
#
# **Constraints:**
#
#   * `1 <= text.length <= 10^5`
#   * The string may contain any possible characters out of all the 256 ASCII characters.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end

ext: jl


text:
# using .Bijectors

struct Transform{N,T} <: Bijector{N}
    t::T

    function Transform(t::T) where {T <:  TransformVariables.TransformTuple}
        new{t.dimension, T}(t)
    end
end

function (b::Transform{N,T})(x::NamedTuple) where {N,T <: TransformVariables.AbstractTransform}
    inverse(b.t, x)
end

function (b::Inverse{<: Transform{N,T}})(y::AbstractVector) where {N,T <:  TransformVariables.AbstractTransform}
    b.orig.t(y)
end

Bijectors.bijector(d::ModelClosure) = Transform(xform(d))

function Bijectors.logabsdetjac(b::Transform, x::NamedTuple)
    -transform_and_logjac(b.t, b.t(x))[2]
end


text:
@testset "991.broken-calculator.jl" begin
    @test broken_calc(2, 3) == 2
    @test broken_calc(3, 10) == 3
    @test broken_calc(5, 8) == 2
    @test broken_calc(1024, 1) == 1023
end


text:
module TestMultiReturningFunction4A

using Modia
using Modia.StaticArrays
@usingModiaPlot

mutable struct Mbs
    derw1
    derw2
end

function twoCoupledInertias(mbs::Mbs,J1,J2,tau0)
    tau1 = tau0 - J1*mbs.derw1
    tau2 = tau1 - J2*mbs.derw2
    return (tau1,tau2)
end

ThreeCoupledInertias = Model(
    J1 = 1.1,
    J2 = 1.2,
    J3 = 1.3,
    phi1 = Var(init=0.0),
    phi2 = Var(init=0.0),
    phi3 = Var(init=0.0),
    w1   = Var(init=0.0),
    w2   = Var(init=0.0),
    w3   = Var(init=0.0),
    tau1 = Var(start=0.0),
    tau2 = Var(start=0.0),
    equations = :[
        tau0 = sin(time/u"s")
        w1 = der(phi1)
        w2 = der(phi2)
        w3 = der(phi3)
        phi2 = phi1
        phi3 = phi2
        mbs = Mbs(der(w1), der(w2))
        (tau1,tau2) = twoCoupledInertias(mbs,J1,J2,tau0)
        J3*der(w3) = tau2
    ]
)

threeCoupledInertias = @instantiateModel(ThreeCoupledInertias, unitless=true, log=false, logDetails=false, logCode=true, logStateSelection=false)
simulate!(threeCoupledInertias, stopTime = 2.0, log=true, requiredFinalStates=[0.3933746028781301, 0.3029735305821084])

plot(threeCoupledInertias, [("phi1", "phi2", "phi3"),
                            ("w1","w2","w3"),
                            ("tau0", "tau1", "tau2")])

end


text:
const test_3d_inputfile = "input_vals_3d.jl"  # input file used by all tests

"""
  Test converting back and forth between conservative and entropy variables
  in 3d
"""
function test_3d_conversion(mesh, sbp, eqn, opts)

# create entropy variable param type
  params_e = EulerEquationMod.ParamType{3, :entropy, Float64, Float64, Float64}(mesh, sbp, opts, 1)

  @testset "----- Testing Conversion -----" begin

    # test unsafe conversion kernels
    q = [1., 2, 3, 4, 15]
    q2 = zeros(q)
    EulerEquationMod.convertToEntropy_(eqn.params, q, q2)
    q3 = zeros(q)
    EulerEquationMod.convertToConservative_(eqn.params, q2, q3)
    @test isapprox( q3, q) atol=1e-13

    vIR = copy(q)
    EulerEquationMod.convertToIR(eqn.params, vIR, vIR)
    diff = vIR - q2./eqn.params.gamma_1
    @test isapprox( norm(diff), 0.0) atol=1e-12

    EulerEquationMod.convertToConservativeFromIR_(eqn.params, vIR, vIR)
    diff = vIR - q
    @test isapprox( norm(diff), 0.0) atol=1e-12

    # test in-place conversion
    q2 = copy(q)
    EulerEquationMod.convertToEntropy_(eqn.params, q2, q2)
    EulerEquationMod.convertToConservative_(eqn.params, q2, q2)
    @test isapprox( q2, q) atol=1e-13

    # test safe interface
    EulerEquationMod.convertToEntropy(eqn.params, q, q2)
    EulerEquationMod.convertToEntropy_(eqn.params, q, q3)
    @test isapprox( q2, q3) atol=1e-13

    EulerEquationMod.convertToConservative(eqn.params, q, q2)
    @test isapprox( q2, q) atol=1e-13

    EulerEquationMod.convertToEntropy_(eqn.params, q, q2)
    EulerEquationMod.convertToEntropy(params_e, q2, q3)
    @test isapprox( q3, q2) atol=1e-13

    EulerEquationMod.convertToConservative(params_e, q2, q3)
    @test isapprox( q3, q) atol=1e-13
  end  # end facts block

  return nothing
end  # end function

#test_3d_conversion(mesh, sbp, eqn, opts)
add_func2!(EulerTests, test_3d_conversion,  test_3d_inputfile, [TAG_ENTROPYVARS, TAG_SHORTTEST])

function test_3d_eigensystem(mesh, sbp, eqn, opts)

  @testset "----- Testing Eigensystem -----" begin

    params = eqn.params
    q = [1., 2, 3, 4, 15]
    numDofPerNode = length(q)
    qc = convert(Array{Complex128}, q)
    aux_vars = Array{Complex128}(1)
    qg = deepcopy(q)
    dxidx = mesh.dxidx[:, :, 1, 1]  # arbitrary
    F = zeros(Complex128, numDofPerNode)

    Ax = zeros(numDofPerNode, numDofPerNode)
    Ax3 = zeros(Ax)
    Ay = zeros(Ax)
    Ay3 = zeros(Ay)
    Az = zeros(Ax)
    Az3 = zeros(Ax)

    EulerEquationMod.calcA(params, q, Ax3)
    EulerEquationMod.calcB(params, q, Ay3)
    EulerEquationMod.calcC(params, q, Az3)
    h = 1e-20
    pert = Complex128(0, h)
    for i=1:numDofPerNode
      qc[i] += pert
      dir = [1.0, 0.0, 0]
      p = EulerEquationMod.calcPressure(params, qc)
      aux_vars[1] = p
      EulerEquationMod.calcEulerFlux(params, qc, aux_vars, dir, F)
      Ax[:, i] = imag(F)/h

      dir = [0.0, 1.0, 0.0]
      EulerEquationMod.calcEulerFlux(params, qc, aux_vars, dir, F)
      Ay[:, i] = imag(F)/h

      dir = [0.0, 0.0, 1.0]
      EulerEquationMod.calcEulerFlux(params, qc, aux_vars, dir, F)
      Az[:, i] = imag(F)/h

      qc[i] -= pert
    end

    # now compute Ax and and Ay from their eigensystem and compare
    Yx = zeros(Ax)
    Yy = zeros(Ax)
    Yz = zeros(Ax)
    Lambdax = zeros(numDofPerNode)
    Lambday = zeros(numDofPerNode)
    Lambdaz = zeros(numDofPerNode)

    EulerEquationMod.calcEvecsx(params, q, Yx)
    EulerEquationMod.calcEvecsy(params, q, Yy)
    EulerEquationMod.calcEvecsz(params, q, Yz)
    EulerEquationMod.calcEvalsx(params, q, Lambdax)
    EulerEquationMod.calcEvalsy(params, q, Lambday)
    EulerEquationMod.calcEvalsz(params, q, Lambdaz)

    Ax2 = Yx*diagm(Lambdax)*inv(Yx)
    Ay2 = Yy*diagm(Lambday)*inv(Yy)
    Az2 = Yz*diagm(Lambdaz)*inv(Yz)

    @test isapprox( Ax2, Ax) atol=1e-12
    @test isapprox( Ay2, Ay) atol=1e-12
    @test isapprox( Az2, Az) atol=1e-12
    @test isapprox( Ax3, Ax) atol=1e-12
    @test isapprox( Ay3, Ay) atol=1e-12
    @test isapprox( Az3, Az) atol=1e-12

    # check A0 = Y*S2*Y.'
    A0 = zeros(Ax)
    Sx = zeros(numDofPerNode)
    Sy = zeros(numDofPerNode)
    Sz = zeros(numDofPerNode)

    EulerEquationMod.getIRA0(params, q, A0)
    EulerEquationMod.calcEScalingx(params, q, Sx)
    EulerEquationMod.calcEScalingy(params, q, Sy)
    EulerEquationMod.calcEScalingz(params, q, Sz)
    A02 = Yx*diagm(Sx)*Yx.'
    A03 = Yy*diagm(Sy)*Yy.'
    A04 = Yz*diagm(Sz)*Yz.'

    @test isapprox( A0, A02) atol=1e-12
    @test isapprox( A0, A03) atol=1e-12
    @test isapprox( A0, A04) atol=1e-12

  end  # end facts block

  return nothing
end

add_func2!(EulerTests, test_3d_eigensystem,  test_3d_inputfile, [TAG_ENTROPYVARS, TAG_SHORTTEST])

"""
  Test calculation of Euler flux in 3D, for both conservative and entropy
  variables.
"""
function test_3d_flux(mesh, sbp, eqn, opts)
  params_e = EulerEquationMod.ParamType{3, :entropy, Float64, Float64, Float64}(mesh, sbp, opts, 1)

  @testset "----- Testing Flux Calculation -----" begin
    q = [1., 2, 3, 4, 15]
    F1 = [2, 4.2, 6, 8, 30.4]
    F2 = [3, 6, 9.2, 12, 45.6]
    F3 = [4, 8, 12, 16.2, 60.8]
    q2 = zeros(q)
    EulerEquationMod.convertToEntropy(eqn.params, q, q2)

    p = EulerEquationMod.calcPressure(eqn.params, q)
    p2 = EulerEquationMod.calcPressure(params_e, q2)
    @test isapprox( p, 0.2) atol=1e-12
    @test isapprox( p2, p) atol=1e-10
    aux_vars = [p]

    F = zeros(Float64, 5)
    Fe = zeros(F)
    nrm = [1., 0, 0]
    EulerEquationMod.calcEulerFlux(eqn.params, q, aux_vars, nrm, F)
    EulerEquationMod.calcEulerFlux(params_e, q2, aux_vars, nrm, Fe)
    @test isapprox( F, F1) atol=1e-12
    @test isapprox( Fe, F) atol=1e-10

    nrm = [0., 1, 0]
    EulerEquationMod.calcEulerFlux(eqn.params, q, aux_vars, nrm, F)
    EulerEquationMod.calcEulerFlux(params_e, q2, aux_vars, nrm, Fe)
    @test isapprox( F, F2) atol=1e-12
    @test isapprox( Fe, F) atol=1e-10

    nrm = [0., 0, 1]
    EulerEquationMod.calcEulerFlux(eqn.params, q, aux_vars, nrm, F)
    EulerEquationMod.calcEulerFlux(params_e, q2, aux_vars, nrm, Fe)
    @test isapprox( F, F3) atol=1e-12
    @test isapprox( Fe, F) atol=1e-10
  end  # end facts block

  return nothing
end  # end function

#test_3d_flux(mesh, sbp, eqn, opts)
add_func2!(EulerTests, test_3d_flux,  test_3d_inputfile, [TAG_ENTROPYVARS, TAG_FLUX, TAG_SHORTTEST])

"""
  Test some auxiliary calculation functions in 3D
"""
function test_3d_misc(mesh, sbp, eqn, opts)
  params_e = EulerEquationMod.ParamType{3, :entropy, Float64, Float64, Float64}(mesh, sbp, opts, 1)

  @testset "----- Testing Misc. Functions -----" begin

    q = [1., 2, 3, 4, 15]
    q2 = zeros(q)
    EulerEquationMod.convertToEntropy(eqn.params, q, q2)
    a = EulerEquationMod.calcSpeedofSound(eqn.params, q)
    @test isapprox( a, sqrt(0.28)) atol=1e-13
    a2 = EulerEquationMod.calcSpeedofSound(params_e, q2)
    @test isapprox( a2, a) atol=1e-13

    s = EulerEquationMod.calcEntropy(eqn.params, q)
    @test isapprox( s, log(0.4*0.5)) atol=1e-12
    s2 = EulerEquationMod.calcEntropy(params_e, q2)
    @test isapprox( s2, s) atol=1e-12

  end  # end facts block

  return nothing
end  # end function

#test_3d_misc(mesh, sbp, eqn, opts)
add_func2!(EulerTests, test_3d_misc,  test_3d_inputfile, [TAG_ENTROPYVARS, TAG_MISC, TAG_SHORTTEST])

function test_3d_secondary_quantities()

  fname = "input_vals_3d.jl"
  opts = read_input(fname)

  for p=1:4
    opts["order"] = p
    fname = "input_vals_3d_p$p"
    make_input(opts, fname)
    mesh, sbp, eqn, opts = solvePDE(fname*".jl")

    coords = mesh.coords[:, :, 1]
    dxidx = mesh.dxidx[:, :, :, 1]
    jac = mesh.jac[:, 1]
    q = zeros(mesh.numDofPerNode, mesh.numNodesPerElement)
    vorticity = zeros(3, mesh.numNodesPerElement)

    ufunc = (x, y, z) -> 2*x^p + y^p + z^p + 1
    vfunc = (x, y, z) -> x^p + 2*y^p + z^p + 1
    wfunc = (x, y, z) -> x^p + y^p + 2*z^p + 1

    dudx = (x, y, z) -> 2*p*x^(p-1)
    dudy = (x, y, z) -> 1*p*y^(p-1)
    dudz = (x, y, z) -> 1*p*z^(p-1)

    dvdx = (x, y, z) -> 1*p*x^(p-1)
    dvdy = (x, y, z) -> 2*p*y^(p-1)
    dvdz = (x, y, z) -> 1*p*z^(p-1)

    dwdx = (x, y, z) -> 1*p*x^(p-1)
    dwdy = (x, y, z) -> 1*p*y^(p-1)
    dwdz = (x, y, z) -> 2*p*z^(p-1)

    for itr=0:1
      # populate q
      for i=1:mesh.numNodesPerElement
        x = coords[1, i] + itr
        y = coords[2, i] + itr
        z = coords[3, i] + itr

        q[1, i] = 1
        q[2, i] = ufunc(x, y, z)
        q[3, i] = vfunc(x, y, z)
        q[4, i] = wfunc(x, y, z)
        q[5, i] = 10
      end
      EulerEquationMod.calcVorticity(eqn.params, sbp, q, dxidx, jac, vorticity)

      for i=1:mesh.numNodesPerElement
        x = coords[1, i] + itr
        y = coords[2, i] + itr
        z = coords[3, i] + itr

        vortx = dwdy(x, y, z) - dvdz(x, y, z)
        vorty = -dwdx(x, y, z) + dudz(x, y, z)
        vortz = dvdx(x, y, z) - dudy(x, y, z)

        @test isapprox( vortx, vorticity[1, i]) atol=1e-12
        @test isapprox( vorty, vorticity[2, i]) atol=1e-12
        @test isapprox( vortz, vorticity[3, i]) atol=1e-12
      end
    end  # end loop itr
  end  # end loop over p

  # the enstrophy integral is O(2p + 1), so for integration to be exact, make
  # solution 2*p - 1
  for pprime = 3:4
    p = pprime - 2
#    p = 1
    opts["order"] = pprime
    fname = "input_vals_3d_p$p"
    make_input(opts, fname)
    mesh, sbp, eqn, opts = solvePDE(fname*".jl")

    ufunc = (x, y, z) -> 2*x^p + y^p + z^p + 1
    vfunc = (x, y, z) -> x^p + 2*y^p + z^p + 1
    wfunc = (x, y, z) -> x^p + y^p + 2*z^p + 1

    q = eqn.q
    for el = 1:mesh.numEl
      for i=1:mesh.numNodesPerElement
        x = mesh.coords[1, i, el]
        y = mesh.coords[2, i, el]
        z = mesh.coords[3, i, el]

        q[1, i, el] = 1
        q[2, i, el] = ufunc(x, y, z)
        q[3, i, el] = vfunc(x, y, z)
        q[4, i, el] = wfunc(x, y, z)
        q[5, i, el] = 10
      end
    end

    # ----- test enstrophy ------
    enstrophy_exact = (24*3^p - (24*p^2)/(2*p - 1) - 12*3^(2*p) + (8*3^(2*p)*p^2)/(2*p - 1) - 12)/(2*mesh.volume)

    enstrophy_numerical = EulerEquationMod.calcEnstrophy(mesh, sbp, eqn, opts, q)

    @test isapprox( enstrophy_exact, enstrophy_numerical) atol=1e-12

    #----- test Kinetic Energy -----
    q_vec = reshape(q, mesh.numDof)

    if p == 1
      ke_exact = 1992
    elseif p== 2
      ke_exact = 132712/15
    end

    ke_exact /= 2*mesh.volume

    ke_numerical = EulerEquationMod.calcKineticEnergy(mesh, sbp, eqn, opts, q_vec)

    @test isapprox( ke_exact, ke_numerical) atol=1e-12

    # ----- test kinetic energy dissipation rate ------
    t = 0.5
    ufunc = (x, y, z) -> 2*x^p + y^p + z^p + 1 + t^p
    vfunc = (x, y, z) -> x^p + 2*y^p + z^p + 1 + t^p
    wfunc = (x, y, z) -> x^p + y^p + 2*z^p + 1 + t^p
    rhofunc = (x, y, z) -> t^p

    dudt = (x, y, z) -> p*t^(p-1)
    dvdt = (x, y, z) -> p*t^(p-1)
    dwdt = (x, y, z) -> p*t^(p-1)
    drhodt = (x, y, z) -> p*t^(p-1)

    q = eqn.q
    res = zeros(q)
    for el = 1:mesh.numEl
      for i=1:mesh.numNodesPerElement
        x = mesh.coords[1, i, el]
        y = mesh.coords[2, i, el]
        z = mesh.coords[3, i, el]

        rho = rhofunc(x, y, z)
        q[1, i, el] = rho
        q[2, i, el] = rho*ufunc(x, y, z)
        q[3, i, el] = rho*vfunc(x, y, z)
        q[4, i, el] = rho*wfunc(x, y, z)
        q[5, i, el] = 10

        res[1, i, el] = drhodt(x, y, z)
        res[2, i, el] = rhofunc(x,  y, z)*dudt(x, y, z) + ufunc(x, y, z)*drhodt(x, y, z)
        res[3, i, el] = rhofunc(x,  y, z)*dvdt(x, y, z) + vfunc(x, y, z)*drhodt(x, y, z)

        res[4, i, el] = rhofunc(x,  y, z)*dwdt(x, y, z) + wfunc(x, y, z)*drhodt(x, y, z)

        res[5, i, el] = 0

      end
    end

    q_vec = reshape(q, mesh.numDof)
    res_vec = reshape(res, mesh.numDof)

    if p == 1
      dkedt_exact = 24*t + 216
    else
      dkedt_exact = 16*t*(3*t^2 + 55)
    end
    dkedt_exact *= rhofunc(1, 1, 1)/mesh.volume

    dkedt_numerical = EulerEquationMod.calcKineticEnergydt(mesh, sbp, eqn, opts, q_vec, res_vec)

    @test isapprox( dkedt_exact, dkedt_numerical) atol=1e-12

  end  # end loop over p

  return nothing

end

add_func1!(EulerTests, test_3d_secondary_quantities, [TAG_MISC, TAG_SHORTTEST])

"""
  Test calculation of A0
"""
function test_3d_matrices(mesh, sbp, eqn, opts)
  params_e = EulerEquationMod.ParamType{3, :entropy, Float64, Float64, Float64}(mesh, sbp, opts, 1)

  @testset "----- Testing Coefficient Matrices calculations -----" begin
    q = [1., 2, 3, 4, 15]
    q2 = zeros(q)
    EulerEquationMod.convertToEntropy(eqn.params, q, q2)

    A0 = [-4 -8 -12 -16 -60;
           -8 -16.8 -24 -32 -121.6;
           -12 -24 -36.8 -48 -182.4;
           -16 -32 -48 -64.8 -243.2;
           -60 -121.6 -182.4 -243.2 -923.6]
    fac = -0.625
    A0 = A0*fac

    A02 = zeros(A0)
    EulerEquationMod.calcA0(params_e, q2, A02)
    for i=1:length(A0)
      @test isapprox( A0[i], A02[i]) atol=1e-10
    end

    A0inv = inv(A0)
    A02inv = zeros(A02)
    EulerEquationMod.calcA0Inv(params_e, q2, A02inv)
    for i=1:length(A0)
      @test isapprox( A0inv[i], A02inv[i]) atol=1e-8
    end

    # test IRA0 against A0 and a jacobian computed with complex step
    A03 = A0 * eqn.params.gamma_1
    A04_test = zeros(A02)
    A04_code = zeros(A02)
    EulerEquationMod.getIRA0(eqn.params, q, A04_code)
    q3 = zeros(Complex128, length(q))
    q4 = zeros(Complex128, length(q))

    q3[:] = q
    EulerEquationMod.convertToEntropy_(eqn.params, q3, q4)
    for i=1:length(q)
      h = 1e-20
      q4[i] += complex(0, h)
      EulerEquationMod.convertToConservative_(eqn.params, q4, q3)
      A04_test[:, i] = imag(q3)/h
      q4[i] -= complex(0, h)
    end
    scale!(A04_test, eqn.params.gamma_1)

    for i=1:length(A04_test)
      @test isapprox( A04_code[i], A04_test[i]) atol=1e-13
      @test isapprox( A04_code[i], A03[i]) atol=1e-10
    end

    test_IRA0inv(eqn.params, q3)

  end  # end facts block

  return nothing
end  # end function

#test_3d_matrices(mesh, sbp, eqn, opts)
add_func2!(EulerTests, test_3d_matrices,  test_3d_inputfile, [TAG_SHORTTEST, TAG_ENTROPYVARS])

"""
  Test Roe solver in 3D
"""
function test_3d_bc(mesh, sbp, eqn, opts)
  @testset "----- Testing BC Solvers -----" begin

    params = eqn.params
    q = [1., 2, 3, 4, 15]
    F = zeros(q)
    F2 = zeros(q)
    nrm = [1., 1, 1]
    dir = mesh.dxidx[:, :, 1, 1].'*nrm
    p = EulerEquationMod.calcPressure(eqn.params, q)
    aux_vars = [p]
    dxidx = mesh.dxidx[:, :, 1, 1]
    nrm_xy = zeros(mesh.dim)

    calcBCNormal(params, dxidx, nrm, nrm_xy)
    EulerEquationMod.RoeSolver(params, q, q, aux_vars, nrm_xy, F2)
    EulerEquationMod.calcEulerFlux(params, q, aux_vars, dir, F)

    @test isapprox( F2, F) atol=1e-13

    func1 = EulerEquationMod.noPenetrationESBC(mesh, eqn)
    # make velocity parallel to the boundary
    nrm = mesh.nrm_bndry[:, 1, 1]
    tngt = mesh.coords_bndry[:, 2, 1] - mesh.coords_bndry[:, 1, 1]
    tngt = tngt./norm(tngt)
    vval = 3
    q = [1., vval*tngt[1], vval*tngt[2], vval*tngt[3], 15]
    aux_vars = eqn.aux_vars[:, 1, 1]
    coords = mesh.coords_bndry[:, 1, 1]

    func1(params, q, aux_vars, coords, nrm, F2)
    EulerEquationMod.calcEulerFlux(params, q, aux_vars, nrm, F)

    @test isapprox( F2, F) atol=1e-13

    func1 = EulerEquationMod.noPenetrationBC(mesh, eqn)
    func1(params, q, aux_vars, coords, nrm, F2)

    @test isapprox( F2, F) atol=1e-13

  end  # end facts block

  return nothing
end  # end functiona

function test_3d_functional(mesh, sbp, eqn, opts)
  @testset "----- Testing 3D functional -----" begin

    q_bndry_orig = eqn.q_bndry
    q_bndry2  = zeros(eqn.q_bndry)
    eqn.q_bndry = q_bndry2
    s1 = sview(eqn.q_bndry, 1, :, :)
    fill!(s1, 1.0)
    s2 = sview(eqn.q_bndry, 2, :, :)
    fill!(s2, 0.355)
    s3 = sview(eqn.q_bndry, 3, :, :)
    fill!(s3, 0.355)
    s4 = sview(eqn.q_bndry, 4, :, :)
    fill!(s4, 0.355)
    s5 = sview(eqn.q_bndry, 5, :, :)
    fill!(s5, 3.0)

    moment_about = [2., 2, 2]
    mom = EulerEquationMod.calcMomentContribution!(mesh, eqn, [1], moment_about)
    println("moment = ", mom)
    for i=1:length(mom)
      @test isapprox( mom[i], 0.0) atol=1e-13
    end

    # test reverse mode
    moment_about = [0.0, 0, 0]
    nout = 3
    nin = mesh.dim*mesh.numNodesPerFace*mesh.numBoundaryFaces
    jac = zeros(nout, nin)

    mom_0 = EulerEquationMod.calcMomentContribution!(mesh, eqn, [1], moment_about)
    pert = 1e-6
    for i=1:nin
      mesh.nrm_bndry[i] += pert
      mom_i = EulerEquationMod.calcMomentContribution!(mesh, eqn, [1], moment_about)

      for j=1:nout
        jac[j, i] = (mom_i[j] - mom_0[j])/pert
      end

      mesh.nrm_bndry[i] -= pert
    end

    # reverse mode
    jac2 = zeros(jac)
    mom_bar = zeros(mom_0)
    for i=1:nout
      mom_bar[i] = 1
      fill!(mesh.nrm_bndry_bar, 0.0)
      EulerEquationMod.calcMomentContribution_revm!(mesh, eqn, [1], moment_about, mom_bar)

      for j=1:nin
        jac2[i, j] = mesh.nrm_bndry_bar[j]
      end

      mom_bar[i] = 0
    end

    #=
    for i=1:nin
      @test isapprox( norm(jac[:, i] - jac2[:, i])/size(jac, 1), 0.0) atol=1e-5
    end
    =#
    @test isapprox( norm(jac - jac2)/sqrt(length(jac)), 0.0) atol=1e-5

    # test other method
    println("checking second method")
    face_range = mesh.bndry_offsets[1]:(mesh.bndry_offsets[2] - 1)
    bndryfaces = mesh.bndryfaces[face_range]
    coords_faces = mesh.coords_bndry[:, :, face_range]
    nrm = mesh.nrm_bndry[:, :, face_range]

    dforce = EulerEquationMod.computeDForce(mesh, eqn, bndryfaces, nrm)
    mom_0 = EulerEquationMod.calcMomentContribution!(mesh.sbpface, coords_faces, dforce, moment_about)
    nin = length(dforce)
    nout = 3
    jac = zeros(nout, nin)
    for i=1:nin
      dforce[i] += pert
      mom_i = EulerEquationMod.calcMomentContribution!(mesh.sbpface, coords_faces, dforce, moment_about)

      for j=1:nout
        jac[j, i] = (mom_i[j] - mom_0[j])/pert
      end

      dforce[i] -= pert
    end

    jac2 = zeros(jac)
    dforce_bar = zeros(dforce)
    coords_bndry_bar = zeros(coords_faces)
    mom_bar = zeros(mom_0)
    for i=1:nout
      mom_bar[i] = 1
      fill!(dforce_bar, 0.0)
      fill!(coords_bndry_bar, 0.0)

      EulerEquationMod.calcMomentContribution_rev!(mesh.sbpface, coords_faces, coords_bndry_bar, dforce, dforce_bar, moment_about, mom_bar)
      for j=1:nin
        jac2[i, j] = dforce_bar[j]
      end

      mom_bar[i] = 0
    end

    #=
    for i=1:size(jac, 2)
      @test isapprox( norm(jac[:, i] - jac2[:, i]), 0.0) atol=1e-5
    end
    =#
    @test isapprox( norm(jac - jac2)/sqrt(length(jac)), 0.0) atol=1e-5

    eqn.q_bndry = q_bndry_orig
  end

  return nothing
end

add_func2!(EulerTests, test_3d_bc,  test_3d_inputfile, [TAG_BC, TAG_SHORTTEST])
add_func2!(EulerTests, test_3d_functional,  test_3d_inputfile, [TAG_REVERSEMODE, TAG_SHORTTEST])


text:
# module Lupin

function vecjoin(elements::Array{E, N}, delim::D)::Union{Vector{typejoin(E, D)}, Vector{Union{E, D}}} where {E, N, D}
    first = true
    T = typejoin(E, D)
    if T === Any
        T = Union{E, D}
    end
    result = Vector{T}()
    for el in elements
        if first
            first = false
        else
            push!(result, delim)
        end
        push!(result, el)
    end
    result
end

function vecjoin(r::AbstractRange{E}, delim::D)::Union{Vector{typejoin(E, D)}, Vector{Union{E, D}}} where {E, D}
    vecjoin(collect(r), delim)
end

# module Lupin


text:
include("tabular_cfr.jl")
include("outcome_sampling_mccfr.jl")
include("external_sampling_mccfr.jl")


text:
#=
title: VRearth( manipulatable virtual earth map )
version: 0.1.0
date: 2020/06/21
author: okimebarun
url: https://github.com/okimebarun/
url: https://qiita.com/oki_mebarun/
=#

module VRMapEarth

greet() = print("Hello World!")

############################################################################
using Interact
using Plots
using Images, ImageView, ImageMagick
using Blink

#=
 calculation part
=#

# parameters
mutable struct ParamDict
    p::Float32  # z-axis of eye point
    r::Float32  # sphere radius
    v0::Int     # mid-point of v
    w0::Int     # mid-point of w
    m0::Int     # mid-point of m
    n0::Int     # mid-point of n
    z1::Float32 # z-axis of img-plane
    z2::Float32 # z-axis of background-plane (dont use now)
    msh::Int    # map size for longtitude
    msv::Int    # map size for latitude
end;

# mapping functions

# get (a,b) from (v,w)
function imcoef(d::ParamDict, v::Int, w::Int)
    t=d.p - d.z1
    a=Float32(v - d.v0)/t
    b=Float32(w - d.w0)/t
    (a,b)
end

# get (flag, x,u,z) from (a,b)
function spxyz(d::ParamDict, a::Float32, b::Float32)
    ab1 = a^2+b^2+1
    r1 = d.p^2 -(d.p^2 - d.r^2)*ab1
    if r1 < 1e-10
        (false, 0.0, 0.0, 0.0)
    else
        t = (d.p - sqrt(r1)) / ab1
        x = a*t
        y = b*t
        z = -t + d.p
        (true, x, y, z)
    end
end

# get (long, lat ,coef) from (x,y,z) and (long0, lat0)
function polerax(d::ParamDict, x::Float32, y::Float32, z::Float32,
        long0::Int, lat0::Int)

    phi = pi*(Float32(lat0)/180.0+1.0)
    ax = 0
    ay = cos(phi)
    az = sin(phi)
    polerv = [ax, ay,  az] # earth's axis
    polerh = [ 0, az, -ay] # orthogonal to polerv

    # get latitude
    r = d.r
    c1 = [x,y,z]'*polerv/r
    lat = acos(c1*0.9999)

    # get longtitude
    # inline function
    cross3d(a,b)=[a[2]b[3]-a[3]b[2], a[3]b[1]-a[1]b[3], a[1]b[2]-a[2]b[1]]
    b  = [x,y,z]/r - polerv*c1 # projection on polerv
    b  = b/sqrt(b'*b) # unitilize
    cv = cross3d(b, polerh) # outer product
    c2 = cv'*polerv
    c3 = b' *polerh
    long = acos(c3*0.9999) * sign(c2)+ pi*Float32(long0)/180.0

    # convert to 360 / 180 degree
    long2 = Float32(360.0*(long/(2pi) + 2.0) % 360.0) # Mod on float
    lat2  = Float32(180.0*( lat/(pi)  + 2.0) % 180.0)

    # reflection coefficient
    coef = [x, y, z]'*[0.0, 0.0, 1.0]/d.r
    (long2, lat2, coef)
end

# get (mx,my) from (long, lat)
function mapxy(d::ParamDict, long::Float32, lat::Float32)
    mx = round(Int, long/360.0*Float32(d.msh-1))+1
    my = round(Int, lat /180.0*Float32(d.msv-1))+1
    (mx, my)
end

# get (m, n) from (v,w)
function backax(d::ParamDict, v::Int, w::Int)
    m = round(Int, ((v-d.v0)*d.m0)/d.v0 + d.m0)+1
    n = round(Int, ((w-d.w0)*d.n0)/d.w0 + d.n0)+1
    (m, n)
end

# get (flag, my, mx, coef) from (w, v, long0, lat0)
function convToMap(d::ParamDict, w::Int, v::Int, long0::Int, lat0::Int)
    (a,b) = imcoef(d, v, w)
    (bon, x, y, z) = spxyz(d, a, b)
    if bon
        (long, lat, coef) = polerax(d, x,y,z, long0, lat0)
        (mx, my) = mapxy(d, long, lat)
        (true, my, mx, coef)
    else
        (m, n) = backax(d, v, w)
        (false, n,  m, 0.0)
    end
end

#=
 GUI part
=#

# main
function real_body()
    # load images
    img_galaxy = load("img_galaxy.jpg")
    img_earth  = load("img_earth.jpg")

    # setup parameter
    d1 = ParamDict(5000.0, 300.0, 400, 300, 500, 300, 510.0, -600.0, 1000, 800)
    d1.msh = size(img_earth)[2]
    d1.msv = size(img_earth)[1]
    d1.m0 = round(Int,size(img_galaxy)[2]/2-1)
    d1.n0 = round(Int,size(img_galaxy)[1]/2-1)

    # view plane
    vs=1:800
    ws=1:600

    # functions
    # convert to image
    function imconv1(p)
        # shade effect
        rgb0 = RGB{N0f8}(0.1,0.1,0.2)
        # select pixel from two images depending on the first flat
        p[1] ? RGB{N0f8}( img_earth[p[2],p[3]]*p[4] + rgb0*(1.0-p[4]) ) : img_galaxy[p[2],p[3]]
    end

    # visualize the result
    function view1(long::Int, lat::Int, z1::Float32, d1=d1, ws=ws, vs=vs)
        d1.z1 = z1
        imap1a(w,v)=convToMap(d1, w, v, long, lat)
        # Julia like mapping!
        m   = imap1a.( ws, vs' )
        img = imconv1.(m)
        colorview(RGB, img)
    end

    # GUI objects
    long = slider(0:10:360,  label="longitude");
    lat  = slider(-90:5:90,  label="latitude ");
    zoom = slider(-5000:100:4500,label="zoom in/out");
    slider01 = vbox(long,lat,zoom);

    # interacting function
    img_plot01=map(
        (long, lat, zoom) ->view1(long, lat, Float32(zoom)),
        map(observe, [long, lat, zoom])...,
    );

    # UI
    ui = vbox(slider01, img_plot01);

    # show window via Blink
    w = Window(Dict(:title=>"VR earth"));
    size(w,818,800);
    body!(w, ui);

    while( active(w))
        sleep(1)
    end
end
#

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

#
function julia_main()
    try
        real_main()
    catch
        Base.invokelatest(Base.display_error, Base.catch_stack())
        return 1
    end
    return 0
end

#
function real_main()
    @show ARGS
    @show Base.PROGRAM_FILE

    # do something
    real_body()
	#

    return
end

if abspath(PROGRAM_FILE) == @__FILE__
    real_main()
end
#

end # module


text:
function digamma(::Type{Gamma}, arg0::jdouble)
    return jcall(Gamma, "digamma", jdouble, (jdouble,), arg0)
end

function gamma(::Type{Gamma}, arg0::jdouble)
    return jcall(Gamma, "gamma", jdouble, (jdouble,), arg0)
end

function inv_gamma1pm1(::Type{Gamma}, arg0::jdouble)
    return jcall(Gamma, "invGamma1pm1", jdouble, (jdouble,), arg0)
end

function lanczos(::Type{Gamma}, arg0::jdouble)
    return jcall(Gamma, "lanczos", jdouble, (jdouble,), arg0)
end

function log_gamma(::Type{Gamma}, arg0::jdouble)
    return jcall(Gamma, "logGamma", jdouble, (jdouble,), arg0)
end

function log_gamma1p(::Type{Gamma}, arg0::jdouble)
    return jcall(Gamma, "logGamma1p", jdouble, (jdouble,), arg0)
end

function regularized_gamma_p(::Type{Gamma}, arg0::jdouble, arg1::jdouble)
    return jcall(Gamma, "regularizedGammaP", jdouble, (jdouble, jdouble), arg0, arg1)
end

function regularized_gamma_p(::Type{Gamma}, arg0::jdouble, arg1::jdouble, arg2::jdouble, arg3::jint)
    return jcall(Gamma, "regularizedGammaP", jdouble, (jdouble, jdouble, jdouble, jint), arg0, arg1, arg2, arg3)
end

function regularized_gamma_q(::Type{Gamma}, arg0::jdouble, arg1::jdouble)
    return jcall(Gamma, "regularizedGammaQ", jdouble, (jdouble, jdouble), arg0, arg1)
end

function regularized_gamma_q(::Type{Gamma}, arg0::jdouble, arg1::jdouble, arg2::jdouble, arg3::jint)
    return jcall(Gamma, "regularizedGammaQ", jdouble, (jdouble, jdouble, jdouble, jint), arg0, arg1, arg2, arg3)
end

function trigamma(::Type{Gamma}, arg0::jdouble)
    return jcall(Gamma, "trigamma", jdouble, (jdouble,), arg0)
end


text:
"Check image headers to determine image format."
module ImgHdr

export imgtype

"Read content of the file and data as `Array{UInt8}`."
function _read(file)
    try
        open(file, "r") do s
            return read(s,String)
        end
    catch
        throw(MethodError("file is not present"))
    end
end

const tests = Function[]

"""GIF ('87 and '89 variants)"""
function isa_gif(f)
    h=_read(f)
    if h[1:6] == "GIF87a" || h[1:6] == "GIF89a"
        return true
    end
    return false
end

push!(tests,isa_gif)

"JPEG data in JFIF or Exif format"
function isa_jpeg(f)
    h=_read(f)
    if h[7:10] == "JFIF" || h[7:10] =="Exif"
        return true
    end
    return false
end

push!(tests,isa_jpeg)

"SGI image library"
function isa_rgb(f)
    h=_read(f)
    if startswith(h,"\001\332")
        return true
    end
    return false
end

push!(tests,isa_rgb)

"TIFF (can be in Motorola or Intel byte order)"
function isa_tiff(f)
    h=_read(f)
    if h[1:2] == "MM" || h[1:2] == "II"
        return true
    end
    return false
end

push!(tests,isa_tiff)

function isa_webp(f)
    h=_read(f)
    if startswith(h, "RIFF") && h[9:12] == "WEBP"
        return true
    end
    return false
end

push!(tests,isa_webp)

function isa_png(f)
    h=_read(f)
    if startswith(h,"\x89PNG\r\n\x1a\n")
        return true
    end
    return false
end

push!(tests,isa_png)

function isa_exr(f)
    h=_read(f)
    if startswith(h,"\x76\x2f\x31\x01")
        return true
    end
    return false
end

push!(tests,isa_exr)

"X bitmap (X10 or X11)"
function isa_xbm(f)
    h=_read(f)
    if startswith(h,"#define ")
        return true
    end
    return false
end
push!(tests,isa_xbm)

"PBM (portable bitmap)"
function isa_pbm(f)
    h=_read(f)
    if length(h) >= 3 && h[1] == 'P' && occursin(h[2],"14") && occursin(h[3],"\t\n\r")
        return true
    end
    return false
end

push!(tests,isa_pbm)

"Sun raster file"
function isa_rast(f)
    h=_read(f)
    if startswith(h,"\x59\xA6\x6A\x95")
        return return true
    end
    return false
end

push!(tests,isa_rast)

function isa_bmp(f)
    h=_read(f)
    if startswith(h,"BM")
        return true
    end
    return false
end

push!(tests,isa_bmp)

"PGM (portable graymap)"
function isa_pgm(f)
    h=_read(f)
    if length(h) >= 3 && h[1] == 'P' && occursin(h[2],"25") && occursin(h[3],"\t\n\r")
        return true
    end
    return false
end

push!(tests,isa_pgm)

"PPM (portable pixmap)"
function isa_ppm(f)
    h=_read(f)
    if length(h) >= 3 && h[1] == 'P' && occursin(h[2],"36") && occursin(h[3],"\t\n\r")
        return true
    end
    return false
end

push!(tests,isa_ppm)

# println(length(tests))

"""
       imgtype("path/to/image_file")
       Return a string with image type
       # Examples
       ```
       julia> using ImgHdr

       julia> imgtype("path/to/file.png")
       "png"
       ```
"""
function imgtype(filename::String)
    for f in tests
        p=f(filename)
        if p
            return split(string(f), '_')[2]
        end
    end
    return nothing  # just make it explicit
end

end


text:
"""
    SciTail

SciTail is a NLI dataset created from multiple-choice science exams and web sentences.

For details, see the [2018 paper](http://ai2-website.s3.amazonaws.com/publications/scitail-aaai-2018_cameraready.pdf) "SciTail: A Textual Entailment Dataset from Science Question Answering"
by Tushar Khot, Asish Sabharwal, and Peter Clark.

```julia
SciTail.train_tsv()
SciTail.train_jsonl()
SciTail.dev_tsv()
SciTail.dev_jsonl()
SciTail.test_tsv()
SciTail.test_jsonl()
```
"""
module SciTail

import ...register_data
using DataDeps: @datadep_str

scitail_base() = joinpath(datadep"SciTailV1.1", "SciTailV1.1")
scitail_jsonl(filename) = joinpath(scitail_base(), "snli_format", filename)
scitail_tsv(filename)   = joinpath(scitail_base(), "tsv_format", filename)

train_tsv()   = scitail_tsv("scitail_1.0_train.tsv")
train_jsonl() = scitail_jsonl("scitail_1.0_train.txt")
dev_tsv()     = scitail_tsv("scitail_1.0_dev.tsv")
dev_jsonl()   = scitail_jsonl("scitail_1.0_dev.txt")
test_tsv()    = scitail_tsv("scitail_1.0_test.tsv")
test_jsonl()  = scitail_jsonl("scitail_1.0_test.txt")

function __init__()
    register_data(
        "SciTailV1.1",
        """
        The SciTail dataset is an entailment dataset created from multiple-choice science exams and web sentences.

        See:
        http://data.allenai.org/scitail/
        """,
        "http://data.allenai.org.s3.amazonaws.com/downloads/SciTailV1.1.zip",
        "3fccd37350a94ca280b75998568df85fc2fc62843a3198d644fcbf858e6943d5")
end

end # module


text:
@testset "rendezvous.jl" begin

    # @testset "experiment" begin

        @process lock() begin
            while true
                @accept caller lock(i::Integer) begin
                    println("locked for p1($i)")
                end
                @accept caller unlock()
            end
        end

        the_lock = Nothing

        @process p1(i::Integer) begin
            println("$(current_time()): Process p1($i) started")
            @send the_lock lock(i)
            println("$(current_time()): Process p1($i) acquired lock")
            work(rand(Uniform(1.0, 10.0)))
            println("$(current_time()): Process p1($i) released lock")
            @send the_lock unlock()
            println("$(current_time()): Process p1($i) ended")
        end

        @simulation begin
            # current_trace!(true)
            global the_lock = @schedule now lock()
            for i = 1:10
                @schedule at rand(0:10) p1(i)
            end
            start_simulation()
        end
    # end
    # @testset "experiment" begin

    #     function Hello end

    #     @process MessageCenter() begin
    #         while true
    #             @accept caller Hello(s::String) begin
    #                 println(s)
    #             end
    #         end
    #     end

    #     message_center = Nothing

    #     @process p(i::Integer) begin
    #         @send message_center Hello("Hello from $i !")
    #     end

    #     @simulation begin
    #         global message_center = @schedule now MessageCenter()
    #         for i = 1:10
    #             @schedule at rand(1.0:10.0) p(i)
    #         end
    #         start_simulation()
    #     end
    # end
end


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

name = "YAJL"
version = v"2.1.0"

# Collection of sources required to build YAJL
sources = [
    GitSource("https://github.com/lloyd/yajl.git",
              "a0ecdde0c042b9256170f2f8890dd9451a4240aa"),
    DirectorySource("bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd yajl

if [[ "$target" == *-w64-* ]]; then
  atomic_patch -p1 "$WORKSPACE/srcdir/patches/01-fix-windows-compiler.patch"
fi

mkdir build
cd build

cmake \
  -DCMAKE_INSTALL_PREFIX="$prefix" \
  -DCMAKE_TOOLCHAIN_FILE="$CMAKE_TARGET_TOOLCHAIN" \
  -DCMAKE_BUILD_TYPE="Release" \
  ..

make "-j$nproc"
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libyajl", :libyajl),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil