text:
const VECTOR_TYPE = true
const SCALAR_TYPE = false

struct Argument{T,U,V}
    name::String
    flags::Vector{String}
    nargs::Union{Int,Symbol}
    default_flag::Optional{String}
    action::Optional{Symbol}
    constant::Optional{T}
    default::U
    choices::Optional{AbstractVector{T}}
    required::Bool
    help::Optional{String}
    metavar::Optional{String}
    dest::Symbol
end

function Argument(name_or_flags::Union{Symbol,String}...;
    type::DataType = Nothing,
    nargs::Union{Int,Char,String,Nothing} = nothing,
    action::Union{Symbol,String,Nothing} = nothing,
    constant::R = nothing,
    default::Union{Vector{S},S} = nothing,
    choices::Union{AbstractVector{T},T} = nothing,
    required::Bool = false,
    help::Optional{String} = nothing,
    metavar::Optional{String} = nothing,
    dest::Union{Symbol,String,Nothing} = nothing) where {R,S,T}

    name, default_flag, flags = parse_name_flags(name_or_flags)
    if dest === nothing
        dest = Symbol(name)
    end

    action = to_symbol(action)
    dest = to_symbol(dest)
    nargs = to_symbol(nargs)

    validate_action(action)

    if type === Nothing
        if action in [:store_true, :store_false]
            type = Bool
        elseif action === :count
            type = Int
        else
            type = coalesce_promote_types(R, S, T)
        end
    end

    if default === nothing
        if action in [:store_true, :store_false]
            default = action === :store_true ? false : true
        elseif action === :count
            default = 0
        end
    end

    nargs = get_nargs(nargs, type, action)
    validate_nargs(nargs)

    if isempty(flags)
        # This is a positional argument, so check if it is required
        required = nargs !== :? && nargs !== :*
    end

    dest_is_vector = need_vector(action, nargs, default)

    validate_args(action, nargs, constant, default, type, choices, dest_is_vector)

    Argument{type,typeof(default),dest_is_vector}(
        name,
        flags,
        nargs,
        default_flag,
        action,
        constant,
        default,
        choices,
        required,
        help,
        metavar,
        dest)
end

to_symbol(::Nothing) = nothing
to_symbol(n::Number) = n
to_symbol(s) = Symbol(s)

function need_vector(action, nargs, default)
    return action in [:append, :append_const] || multiple_values(nargs) || default isa AbstractVector
end

multiple_values(nargs::Integer) = nargs > 1
multiple_values(nargs::Symbol) = nargs in [:+, :*]

function parse_name_flags(name_or_flags::Tuple{Vararg{Union{Symbol,String}}})
    if length(name_or_flags) === 1 && Base.isidentifier(name_or_flags[1])
        name = String(name_or_flags[1])
        return name, nothing, String[]
    end

    flag_names, flags = extract_flag_names(name_or_flags)
    default_flag = flags[1]
    name = longest(flag_names)

    return name, default_flag, flags
end

function extract_flag_names(input_flags)
    flags = String[]
    flag_names = String[]

    for flag in input_flags
        if !(flag isa AbstractString)
            throw(ArgumentError("Invalid flag: `$flag`"))
        end

        dash_count, rest = remove_dashes(flag)

        if !(1 <= dash_count <= 2) || length(rest) === 0
            throw(ArgumentError("Invalid flag: `$flag`"))
        end

        push!(flag_names, replace(rest, '-' => '_'))
        push!(flags, flag)
    end

    return flag_names, flags
end

function coalesce_promote_types(types::DataType...)
    length(types) === 0 && return String

    ret_type = Nothing

    for type in types
        if type !== Nothing
            ret_type = ret_type === Nothing ? type : promote_type(ret_type, type)
        end
    end

    ret_type === Nothing && return String

    return ret_type
end

function validate_args(action, nargs, constant, default, type, choices, dest_is_vector)
    @nospecialize

    if action in [:store_true, :store_false, :store_const, :append_const, :count]
        # Test that nargs and action are consistent
        if nargs !== nothing && nargs !== 0
            throw(ArgumentError("nargs=`$nargs` cannot be specified with action=`$action`"))
        end

        # Test that action and choices are consistent
        if choices !== nothing
            throw(ArgumentError("You cannot specify `choices` for action=`$action`"))
        end
    end

    # Test that nargs and constant are consistent
    if action in [:store_const, :append_const] && constant === nothing
        throw(ArgumentError("Please provide a value for `constant` to use with $action"))
    end

    # Test that default is a vector if necessary
    if default !== nothing && dest_is_vector && !(default isa AbstractVector)
        throw(ArgumentError("Default value must be a vector"))
    end

    # Test that default and choices are consistent
    if choices !== nothing && default !== nothing
        if dest_is_vector
            # We made sure that default was a vector above
            for value in default
                test_value_in_choices(value, choices)
            end
        else
            test_value_in_choices(default, choices)
        end
    end

end

function test_value_in_choices(value, choices)
    if !(value in choices)
        throw(ArgumentError("Default value `$value` must be one of $choices"))
    end
end

validate_action(action) = nothing
function validate_action(action::Symbol)
    if !(action in [:store_true, :store_false, :store_const, :append, :append_const, :count, :help])
        throw(ArgumentError("Invalid action: `$(String(action))`"))
    end
end

validate_nargs(nargs) = nothing
function validate_nargs(nargs::Symbol)
    if !(nargs in [:?, :+, :*])
        throw(ArgumentError("Invalid specifier `$nargs` for `nargs`"))
    end
end

get_nargs(nargs::Nothing, ::Type{Bool}, action) = 0
get_nargs(nargs::Nothing, ::Type{<:Vector}, action) = :*
get_nargs(nargs::Nothing, _, action) = action in [:store_true, :store_false, :store_const, :append_const, :count, :help] ? 0 : 1
get_nargs(nargs::Union{Int,Symbol}, _, _) = nargs

function arg_name(arg::Argument, to_uppercase = true)
    name = to_uppercase ? uppercase(arg.name) : arg.name
    return something(arg.metavar, name)
end


text:
println("these tests will print warnings: ")

module PathError
    import FileIO: File, @format_str
    save(file::File{format"PATHERROR"}, data) = nothing
    load(file::File{format"PATHERROR"}) = nothing
end
add_format(format"PATHERROR", (), ".patherror", [PathError])

@testset "Path errors" begin
    # handling a nonexistent parent directory, during save
    temp_dir = joinpath(mktempdir(), "dir_that_did_not_exist")
    @assert !isdir(temp_dir) "Testing error. This dir shouldn't exist"
    fn = joinpath(temp_dir, "file.patherror")
    save(fn, "test content")
    @test isdir(temp_dir)

    # handling a filepath that's an existing directory, during save
    @test_throws ArgumentError save(format"PATHERROR", mktempdir(), "test content")

    # handling a nonexistent filepath, during load
    @test_throws ArgumentError load(joinpath(mktempdir(), "dummy.patherror"))
end

@testset "Not installed" begin
    # The logging part of this test only works when running under Pkg.test (TODO: fix this)
    @test_logs (:warn, r"supply `pkg` as a Module or `name=>uuid`") @test_throws(ArgumentError, add_format(format"NotInstalled", (), ".not_installed", [:NotInstalled]))
    # Give it a fake UUID
    add_format(format"NotInstalled", (), ".not_installed", [:NotInstalled=>UUID("79e393ae-7a7b-11eb-1530-bf4d98024096")])
    @test_throws CapturedException save("test.not_installed", nothing)

    # Core.eval(Base, :(is_interactive = true)) # for interactive error handling
    # add_format(format"NotInstalled", (), ".not_installed", [:NotInstalled])
    # stdin_copy = stdin
    # stderr_copy = stderr
    # rs, wr = redirect_stdin()
    # rserr, wrerr = redirect_stderr()
    # ref = @async save("test.not_installed", nothing)
    # println(wr, "y")
    # @test_throws ArgumentError fetch(ref) #("unknown package NotInstalled")
    # ref = @async save("test.not_installed", nothing)
    # println(wr, "invalid") #test invalid input
    # println(wr, "n") # don't install
    # fetch(ref)
    # @test istaskdone(ref)

    # close(rs);close(wr);close(rserr);close(wrerr)
    # redirect_stdin(stdin_copy)
    # redirect_stderr(stderr_copy)

    # Core.eval(Base, :(is_interactive = false)) # for interactive error handling

end

# Missing load/save functions
module BrokenIO
    using FileIO
end
add_format(format"BROKEN", (), ".brok", [BrokenIO])

@testset "Absent implementation" begin
    stderr_copy = stderr
    rserr, wrerr = redirect_stderr()
    ex = try load(Stream{format"BROKEN"}(stdin)) catch e e end
    @test isa(ex, CapturedException)
    @test isa(ex.ex, FileIO.LoaderError)
    @test isa(ex.ex.ex, FileIO.SpecError)
    ex = try save(Stream{format"BROKEN"}(stdout), nothing) catch e e end
    @test isa(ex, CapturedException)
    @test isa(ex.ex, FileIO.WriterError)
    @test isa(ex.ex.ex, FileIO.SpecError)
    redirect_stderr(stderr_copy)
    close(rserr);close(wrerr)
end

module BrokenIO2
    using FileIO
    function save(f::File{format"BROKEN2"}, data)
        error("whoops")
    end
end
add_format(format"BROKEN2", (), ".brok2", [BrokenIO2])

@testset "Broken implementation/limit" begin
    io = IOBuffer()
    stderr_copy = stderr
    rserr, wrerr = redirect_stderr()
    ex = try save(Stream{format"BROKEN2"}(stdout), [1:1000;]) catch e e end
    @test isa(ex, CapturedException)
    @test isa(ex.ex, MethodError)
    showerror(io, ex)
    str = String(take!(io))
    @test startswith(str, "MethodError: no method matching save(::Stream{DataFormat{:BROKEN2}")
    exw = FileIO.WriterError("BrokenIO2", "a message", ex.ex)
    showerror(io, exw)
    str = String(take!(io))
    @test startswith(str, "BrokenIO2 writer error: a message\n  due to MethodError(")
    @test !occursin("499", str)   # printing should be with `:limit`
    @test  occursin("…", str)     # printing should be with `:limit`
    ex = try save(File{format"BROKEN2"}(tempname()), [1:1000;]) catch e e end
    @test isa(ex, CapturedException)
    @test ex.ex == ErrorException("whoops")
    showerror(io, ex)
    str = String(take!(io))
    @test !occursin("499", str)   # printing should be with `:limit`
    redirect_stderr(stderr_copy)
    close(rserr);close(wrerr)
end

module MultiError1
    import FileIO: @format_str, File
    load(file::File{format"MultiError"}) = error("1")
end
module MultiError2
    import FileIO: @format_str, File, magic
    load(file::File{format"MultiError"}) = error("2")
end
@testset "multiple errors" begin
    println("this test will print warnings: ")
    add_format(
        format"MultiError",
        (),
        ".multierr",
        [MultiError1],
        [MultiError2]
    )
    tmpfile = joinpath(mktempdir(), "test.multierr")
    open(tmpfile, "w") do io
        println(io, "dummy content")
    end
    ex = try load(tmpfile) catch e e end
    @test isa(ex, CapturedException)
    @test isa(ex.ex, ErrorException)
    rm(dirname(tmpfile), recursive=true)
end

@testset "stacktrace trimming" begin
    sfs = try
        error("whoops")
    catch err
        stacktrace(catch_backtrace())
    end
    n = length(sfs)
    FileIO.trim!(sfs)
    @test length(sfs) < 5 < n
end

module BustedME
using FileIO
load(file, args...) = foo()
add_format(format"BUSTEDME", "BUSTEDME", ".bus", [BustedME])
end
@testset "error accuracy" begin
    tmpfile = joinpath(mktempdir(), "test.bus")
    open(tmpfile, "w") do io
        write(io, magic(format"BUSTEDME"))
        println(io, "dummy content")
    end
    ex = try load(tmpfile) catch e e end
    @test isa(ex, CapturedException)
    @test isa(ex.ex, UndefVarError)
    rm(dirname(tmpfile), recursive=true)
end

ext: jl


text:
export plot, plot_rectangle, plot_circle, plot_ellipse, plot_bounds, plot_graph, plot_tree, plot_path, plot_line_segments, save_plot

rectangle_corners(lo, hi) = ([lo[1],hi[1],hi[1],lo[1],lo[1]], [lo[2],lo[2],hi[2],hi[2],lo[2]])

function plot_rectangle(lo, hi; xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf, kwargs...)
    X, Y = rectangle_corners(lo, hi)
    plt.fill(clamp(X, xmin, xmax),
             clamp(Y, ymin, ymax),
             edgecolor="black", zorder=0; kwargs...)
end

function plot_polygon(pts; xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf, kwargs...)
    XY = statevec2mat(pts)'
    plt.fill(clamp(XY[:,1], xmin, xmax),
             clamp(XY[:,2], ymin, ymax),
             edgecolor="black", zorder=0; kwargs...)
end

function plot_circle(c, r; xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf, kwargs...)
    plt.fill(clamp(c[1] + r*cos(linspace(0, 2pi, 40)), xmin, xmax),
             clamp(c[2] + r*sin(linspace(0, 2pi, 40)), ymin, ymax),
             edgecolor="black", zorder=0; kwargs...)
end

function plot_ellipse(c, a, b, t; xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf, kwargs...)
    XY = [a*cos(linspace(0, 2pi, 40)) b*sin(linspace(0, 2pi, 40))]*[cos(t) sin(t); -sin(t) cos(t)]
    plt.fill(clamp(c[1] + XY[:,1], xmin, xmax),
             clamp(c[2] + XY[:,2], ymin, ymax),
             edgecolor="black", zorder=0; kwargs...)
end

function plot_ellipse(c, Sigma; xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf, kwargs...)
    XY = [cos(linspace(0, 2pi, 40)) sin(linspace(0, 2pi, 40))]*chol(full(Sigma))
    plt.fill(clamp(c[1] + XY[:,1], xmin, xmax),
             clamp(c[2] + XY[:,2], ymin, ymax),
             edgecolor="black", zorder=0; kwargs...)
end

function plot_bounds(lo = zeros(2), hi = ones(2))
    plt.plot(rectangle_corners(lo, hi)..., color="black", linewidth=1.0, linestyle="-")
    plt.axis("equal")
end

function plot_graph(V::Matrix, F; kwargs...)  # learn how to just pass kwargs
    plt.scatter(V[1,:], V[2,:], zorder=1; kwargs...)
    X = vcat([V[1,idx_list] for idx_list in findn(triu(F))]..., fill(nothing, 1, sum(triu(F))))[:]
    Y = vcat([V[2,idx_list] for idx_list in findn(triu(F))]..., fill(nothing, 1, sum(triu(F))))[:]
    plt.plot(X, Y, linewidth=.5, linestyle="-", zorder=1; kwargs...)
end
plot_graph(V::Vector, F; kwargs...) = plot_graph(statevec2mat(V), F; kwargs...)

function plot_tree(V::Matrix, A; kwargs...)
    plt.scatter(V[1,:], V[2,:], zorder=1; kwargs...)
    X = vcat(V[1,find(A)], V[1,A[find(A)]], fill(nothing, 1, countnz(A)))[:]
    Y = vcat(V[2,find(A)], V[2,A[find(A)]], fill(nothing, 1, countnz(A)))[:]
    plt.plot(X, Y, linewidth=.5, linestyle="-", zorder=1; kwargs...)
end
plot_tree(V::Vector, A; kwargs...) = plot_tree(statevec2mat(V), A; kwargs...)

function plot_line_segments(P1::Vector, P2::Vector; kwargs...)
    X = [[v[1] for v in P1]'; [v[1] for v in P2]'; fill(nothing, 1, length(P1))][:]
    Y = [[v[2] for v in P1]'; [v[2] for v in P2]'; fill(nothing, 1, length(P1))][:]
    plt.plot(X, Y; kwargs...)
end

function plot_path(V::Matrix, idx_list = 1:size(V,2); kwargs...)
    plt.plot(V[1,idx_list], V[2,idx_list], linewidth=1.0, linestyle="-", zorder=2; kwargs...)
end
plot_path(V::Vector, idx_list = 1:length(V); kwargs...) = plot_path(statevec2mat(V), idx_list; kwargs...)

function save_plot(fname, title_string = "")
    plt.axis("off")
    plt.axes()[:get_xaxis]()[:set_visible](false)
    plt.axes()[:get_yaxis]()[:set_visible](false)
    title_string != "" && plt.title(title_string, fontsize=20)
    plt.savefig(fname, bbox_inches="tight")
end

# using Compose

# function composable(C::Circle, opts = (fill("red"), fillopacity(1.)))
#     compose(context(), circle(C.c..., C.r), opts...)
# end

# function composable(P::Polygon, opts = (fill("red"), fillopacity(1.)))
#     compose(context(), polygon([(p[1],p[2]) for p in P.points]), opts...)
# end

# function composable(C::Compound2D, opts = (fill("red"), fillopacity(1.)))
#     compose(context(), [composable(P, opts) for P in C.parts]...)
# end

# function composable(L::Line, opts = (stroke("black"),))
#     compose(context(), line([(L.v[1], L.v[2]), (L.w[1], L.w[2])]), opts...)
# end

# function composable(SS::StateSpace, opts = (linewidth(2px), fill(nothing), stroke("black")); flipy=true)
#     w, h = SS.hi - SS.lo
#     compose(context(-SS.lo[1]/w, -SS.lo[2]/h, 1/w, 1/h),
#             compose(context(), rectangle(SS.lo[1], flipy ? SS.hi[2] : SS.lo[2], w, h), opts...))
# end

# flip(C::Context) = compose(context(mirror=Mirror()), C)

# function plot(P::MPProblem)
#     p = compose(composable(P.SS), composable(P.CC))
#     if P.status == :solved
#         println("IMPLEMENT ME!")
#     end
#     flip(p)
# end


text:
function interpolate(obj::TimeInterpolable, arg0::AbsoluteDate, arg1::Collection)
    return jcall(obj, "interpolate", TimeInterpolable, (AbsoluteDate, Collection), arg0, arg1)
end

function interpolate(obj::TimeInterpolable, arg0::AbsoluteDate, arg1::Stream)
    return jcall(obj, "interpolate", TimeInterpolable, (AbsoluteDate, Stream), arg0, arg1)
end


text:
using Base.Test
using EGR

println("TestEGRlikeSG")

createOracleOutputLevel = 1
numEquivalentPasses = 1
algOutputLevel = 0
maxOutputNum=20

myREfunction(problem, opts, sd, expIndices,n) = false
myWriteFunction(problem, sd, opts, k, gnum, origWant, fromOutputsFunction) = false

for (gradientOracle, numVars, numTP, csDataType, LossFunctionString, myOutputter, L2reg, thisDataName, thisProblem) in Oracles

	print(" $thisDataName $LossFunctionString L2reg = $L2reg")

	maxG  = Int(round(numEquivalentPasses*numTP))

	expIndices=unique(round(Int, logspace(0,log10(maxG+1),maxOutputNum)))-1

	myOutputOpts =  OutputOpts(myOutputter,expIndices)

	myOpts(stepSizePower) = Opts(zeros(numVars,1); stepSizePower=stepSizePower, maxG=maxG, outputLevel=algOutputLevel)

	s(k,I) = 0
	u(k,I) = 1
	beta(k) = 1

	a = Task(() -> getSequentialFinite(numTP, gradientOracle))
	b = Task(() -> getSequentialFinite(numTP, gradientOracle))
	(outString, results_k, results_gnum,results_fromOutputsFunction,xFromEGR) = alg(thisProblem(a)	, myOpts(0), USD(s,u,beta,numVars,  "EGR.SGlike.a=1",1,true), myOutputOpts, myWriteFunction)

	(outString, results_k, results_gnum,results_fromOutputsFunction,xFromSG)  = alg(thisProblem(b)	, myOpts(0), SGsd( "SG.a=1")                                               , myOutputOpts, myWriteFunction)

	relError = norm(xFromEGR-xFromSG)/norm(xFromSG)

	println(" relError between EGR and SG = $relError")
	if relError>1e-13
		error("SG and EGR do not coincide")
	end

end

println("TestEGRlikeSG successful!")
println()


text:
using SafeTestsets

@safetestset "Siegler Model" begin
    cd(@__DIR__)
    using Pkg
    Pkg.activate("../")
    using ACTRModels, Test, Distributions, Random, StatsPlots
    include("../models/Siegler/Siegler_Model.jl")
    Random.seed!(187302)
    δ = 16.0
    τ = -.3
    s = .5
    parms = (mmp = true,noise = true,mmp_fun = sim_fun,ter = 2.05)
    stimuli = [(num1 = 1,num2 = 1),(num1 = 1,num2 = 2),(num1 = 1,num2 = 3),(num1 = 2,num2 = 2),
        (num1 = 2,num2 = 3),(num1 = 3,num2 = 3)]
    temp = map(x -> simulate(stimuli, parms; δ, τ, s), 1:500)
    data = vcat(vcat(temp...)...)

    x = range(δ * .8, δ * 1.2, length=50)
    y = map(x -> computeLL(parms, data; δ=x, τ, s), x)
    mxv,mxi = findmax(y)
    δ′ = x[mxi]
    @test δ′ ≈ δ atol = .5

    x = range(τ * .8, τ * 1.2, length=50)
    y = map(x -> computeLL(parms, data; δ, τ=x, s), x)
    mxv,mxi = findmax(y)
    τ′ = x[mxi]
    @test τ′ ≈ τ atol = .1

    x = range(s * .8, s * 1.2, length=50)
    y = map(x -> computeLL(parms, data; δ, τ, s=x), x)
    mxv,mxi = findmax(y)
    s′ = x[mxi]
    @test s′ ≈ s atol = .03
end


text:
# This file has implementations of Galerkin Least-Squares (GLS) method and
# Streamwise Upwind Petrov-Galerkin method in it.

@doc """
### EulerEquationMod.GLS

Add Galerkin Least-Squares stabilization to the weak residual. This
implementation is only for steady problems and conservative variables

**Inputs**

*  `mesh`: AbstractMesh type
*  `sbp` : Summation-by-parts operator
*  `eqn` : Equation object used elsewhere

**Outputs**

* None

"""->
function GLS(mesh::AbstractMesh{Tmsh}, sbp::AbstractOperator,
       eqn::EulerData{Tsol, Tres, Tdim}) where {Tmsh, Tsol, Tres, Tdim}

  #  println("Entered GLS")
  parametricFluxJacobian(mesh, sbp, eqn) # Calculate the euler flux jacobian
  tau = zeros(Tsol, mesh.numDofPerNode, mesh.numDofPerNode,
              mesh.numNodesPerElement, mesh.numEl) # Stabilization term
  calcTau(mesh, sbp, eqn, tau, eqn.params.order)
  # println("Tau = \n", tau)

  gls_res = zeros(eqn.res) # input array for volumeintegrate!

  # Get shape function derivatives
  Hinv = 1./sbp.w
  shapefuncderiv = zeros(Tsol, sbp.numnodes, sbp.numnodes, Tdim)
  for k = 1:Tdim
    for i = 1:sbp.numnodes
      for j = 1:sbp.numnodes
        shapefuncderiv[i,j,k] = Hinv[i]*sbp.Q[i,j,k]
      end
    end
  end

  # Populate gls_res
  for i = 1:mesh.numEl
    endof = mesh.numDofPerNode*mesh.numNodesPerElement # dofs in an element
    ndof = mesh.numDofPerNode
    u = zeros(Tsol,endof) # element level q_vec
    # Populate u
    for j = 1:mesh.numNodesPerElement
      for k = 1:mesh.numDofPerNode
        u[(j-1)*ndof+k] = eqn.q[k,j,i]
      end
    end
    # println("eqn.q = \n", round(eqn.q,2))
    # println("u = \n", round(u,2))
    Axidxi = zeros(Tsol, endof, endof)
    Axi = sview(eqn.Axi,:,:,:,i) # Get flux jacobians for all nodes in an element
    Aeta = sview(eqn.Aeta,:,:,:,i)
    calcAxidxi(Axidxi, shapefuncderiv, Axi, Aeta, ndof, mesh.numNodesPerElement)
    # println("Axidxi = \n", round(Axidxi,2))

    intArr = zeros(Axidxi) # intermediate array for storing tau*Axidxi
    for j = 1:mesh.numNodesPerElement
      m = (j-1)*ndof + 1
      intArr[m:m+ndof-1,:] = sbp.w[j]*mesh.jac[j,i]*tau[:,:,j,i]*Axidxi[m:m+ndof-1,:]
    end

    intvec = Axidxi.'*intArr*u
    for j = 1:mesh.numNodesPerElement
      for k = 1:mesh.numDofPerNode
        gls_res[k,j,i] = intvec[(j-1)*ndof+k]
      end
    end
  end   # end for i = 1:mesh.numEl

  # volumeintegrate!(sbp, gls_res, supg_res)

  # println(supg_res)
  for i = 1:mesh.numEl
    for j = 1:mesh.numNodesPerElement
      for k = 1:mesh.numDofPerNode
        eqn.res[k,j,i] -= gls_res[k,j,i]  # -ve sign because of right hand side.
      end                                 # Fluxes are computed explicitly and
    end                                   # not from weak redidual.
  end

  return nothing
end

@doc """
### EulerEquationMod.parametricFluxJacobian

Claculates the flux jacobian in the parametric space for all the nodes in the
mesh. It uses conservative variables

**Inputs**

*  `mesh` : Abstract mesh object
*  `sbp`  : Summation-by-parts operator
*  `eqn`  : Equation object

**Outputs**

*  None
"""->

function parametricFluxJacobian(mesh::AbstractMesh{Tmsh},
        sbp::AbstractOperator,
        eqn::EulerData{Tsol, Tres, Tdim}) where {Tmsh, Tsol, Tres, Tdim}

# Calculates the Flux Jacobian in the parametric space

  fill!(eqn.Axi, 0.0)   # Reset Axi & Aeta to zero before calculating flux-
  fill!(eqn.Aeta, 0.0)  # jacobaian

  for i = 1:mesh.numEl
    for j = 1:mesh.numNodesPerElement
      q = sview(eqn.q,:, j, i)
      dxidx = sview(mesh.dxidx,:,:,j,i)
      # Calculate the flux jacobian get Ax & Ay
      Ax= zeros(Tsol, 4, 4) # Flux jacobian in the x direction
      Ay = zeros(Ax)        # Flux jacobian in the y direction
      calcFluxJacobian(eqn, q, Ax, Ay)
      # get the fluz jacobian in xi & eta direction
      eqn.Axi[:,:,j,i] = Ax*dxidx[1,1] + Ay*dxidx[1,2]
      eqn.Aeta[:,:,j,i] = Ax*dxidx[2,1] + Ay*dxidx[2,2]
    end  # end for j = 1:mesh.numNodesPerElement
  end    # end for i = 1:mesh.numEL

  return nothing
end

@doc """
### EulerEquationMod.calcFluxJacobian

It computes the Euler flux Jacobian at the nodal level in te physical space
(X & Y axes). It uses conservative variables.

**Inputs**

*  `eqn`  : Equation object
*  `q`    : Conservative variables at the node
*  `Ax`   : Flux Jacobian at a node in the X-direction
*  `Ay`   : Flux Jacobian at a node in the Y-direction

**Outputs**

*  None
"""->
function calcFluxJacobian(eqn::EulerData{Tsol, Tres, Tdim},
        q::AbstractArray{Tsol,1}, Ax::AbstractArray{Tsol,2},
        Ay::AbstractArray{Tsol,2}) where {Tsol, Tres, Tdim}

  gamma_1 = eqn.params.gamma_1
  gamma = eqn.params.gamma
  R = eqn.params.R     # Gas constant
  cv = eqn.params.cv   # Specific heat at constant volume
  u = q[2]/q[1] # Get velocity in the x-direction
  v = q[3]/q[1] # Get velocity in the x-direction
  intvar = (R/cv)*(q[4]/q[1] - 0.5*(u*u + v*v)) # intermediate variable

  # Populating Ax
  Ax[1,1] = 0
  Ax[1,2] = 1
  Ax[1,3] = 0
  Ax[1,4] = 0
  Ax[2,1] = -u*u + 0.5*R*(u*u + v*v)/cv
  Ax[2,2] = 2*u - R*u/cv
  Ax[2,3] = -R*v/cv
  Ax[2,4] = R/cv
  Ax[3,1] = -u*v
  Ax[3,2] = v
  Ax[3,3] = u
  Ax[3,4] = 0
  Ax[4,1] = -q[2]*q[4]/(q[1]*q[1]) - u*intvar + 0.5*u*R*(u*u + v*v)/cv
  Ax[4,2] = q[4]/q[1] + intvar - R*u*u/cv
  Ax[4,3] = -R*u*v/cv
  Ax[4,4] = u + R*u/cv

  # Populating Ay
  Ay[1,1] = 0
  Ay[1,2] = 0
  Ay[1,3] = 1
  Ay[1,4] = 0
  Ay[2,1] = -v*u
  Ay[2,2] = v
  Ay[2,3] = u
  Ay[2,4] = 0
  Ay[3,1] = -v*v + 0.5*R*(u*u + v*v)/cv
  Ay[3,2] = -R*u/cv
  Ay[3,3] = 2*v - R*v/cv
  Ay[3,4] = R/cv
  Ay[4,1] = -q[3]*q[4]/(q[1]*q[1]) - v*intvar + v*(R/cv)*0.5*(u*u + v*v)
  Ay[4,2] = -R*v*u/cv
  Ay[4,3] = q[4]/q[1] + intvar - R*v*v/cv
  Ay[4,4] = v + R*v/cv

  return nothing
end # end function calcFluxJacobian

@doc """
### EulerEquationMod.calcTau

Calculates the stabilization term tau for GLS

**Inputs**

*  `mesh` : Abstract mesh type
*  `sbp`  : Summation-by-parts operator
*  `eqn`  : Euler equation object
*  `tau`  : Stabilization term

**Outputs**

* None
"""->
function calcTau(mesh::AbstractMesh{Tmsh},
sbp::AbstractOperator, eqn::EulerData{Tsol, Tres, Tdim},
tau::AbstractArray{Tsol,4}, order::Integer) where {Tmsh, Tsol, Tres, Tdim}

  # Using Glasby's Method Equation 26 & 27

  # Get shape function derivatives
  Hinv = 1./sbp.w
  shapefuncderiv = zeros(Tsol, sbp.numnodes, sbp.numnodes, Tdim)
  for k = 1:Tdim
    for i = 1:sbp.numnodes
      for j = 1:sbp.numnodes
        shapefuncderiv[i,j,k] = Hinv[i]*sbp.Q[i,j,k]
      end
    end
  end

  for i = 1:mesh.numEl
    for j = 1:mesh.numNodesPerElement
      q = sview(eqn.q,:,j,i)
      dxidx = sview(mesh.dxidx,:,:,j,i) # get mapping jacobian at node j
      for k = 1:mesh.numNodesPerElement
        for l = 1:Tdim
          jac_vector = Tsol[dxidx[l,1];dxidx[l,2]] # get a jacobian vector for eigen value factorization
          T = zeros(Tsol,4,4)
          Tinv = zeros(T)
          Lambda = eye(T)
          calcEigenFactorization(eqn, q, jac_vector, T, Tinv, Lambda)
          tau[:,:,j,i] += T*abs(shapefuncderiv[j,k,l]*Lambda)*Tinv
        end # end for l = 1:Tdim
      end   # for k = 1:mesh.numNodesPerElement
      tau[:,:,j,i] = inv(tau[:,:,j,i])
    end     # end for j = 1:mesh.numNodesPerElement
  end       # end for i = 1:mesh.numEl

  #= Using Masayuki's thesis and Barth's chapter. Using latter's convention

  for i = 1:mesh.numEl
    # Get the normal vectors at the nodes
    for j = 1:mesh.numNodesPerElement
      q = sview(eqn.q,:,j,i)
      dxidx = sview(mesh.dxidx,:,:,j,i)
      jac = sview(mesh.jac,j,i)
      Ax= zeros(Tsol, 4, 4) # Flux jacobian in the x direction
      Ay = zeros(Ax)        # Flux jacobian in the y direction
      calcFluxJacobian(eqn, q, Ax, Ay)
      n = zeros(Tsol,2,2)
      n[:,1] = jac.*[dxidx[1,1];dxidx[1,2]]
      n[:,2] = jac.*[dxidx[2,1];dxidx[2,2]]
      n[:,1] = n[:,1]/sqrt(n[1,1]*n[1,1] + n[2,1]*n[2,1])
      n[:,2] = n[:,2]/sqrt(n[1,2]*n[1,2] + n[2,2]*n[2,2])
      for k = 1:Tdim
        Bn = zeros(Tsol, mesh.numDofPerNode, mesh.numDofPerNode)
        Bn = (n[1,k]*Ax + n[2,k]*Ay)*sqrt(n[1,k]*n[1,k] + n[2,k]*n[2,k])
        E = eigfact(Bn)
        # println("E.values = \n", E.values)
        # println("size of E.values = ", size(E.values))
        absLambda = eye(Bn)
        for l = 1:mesh.numDofPerNode
          absLambda[l,l] = abs(E.values[l])*absLambda[l,l]
        end
        tau[:,:,j,i] += E.vectors*absLambda*inv(E.vectors)
      end
      tau[:,:,j,i] = inv(tau[:,:,j,i])
    end # end for j = 1:mesh.numNodesPerElement
  end   # end for i = 1:mesh.numEl

  =#

  #=  ORIGINAL NOT WORKING TAU

  for i = 1:mesh.numEl
    for j = 1:mesh.numNodesPerElement
      q = sview(eqn.q,:,j,i)
      T = (q[4] - 0.5*(q[2]*q[2] + q[3]*q[3])/q[1])*(1/(q[1]*eqn.params.cv))
      c = sqrt(eqn.params.gamma*eqn.params.R*T)  # Speed of sound
      ux = q[2]/q[1]
      uy = q[3]/q[1]
      h = circumcircleDiameter(mesh.coords[:,:,i])
      # println("h = ", h)
      tau[j,i] = h/(c + sqrt(ux*ux + uy*uy))
    end
  end  =#

  return nothing
end

@doc """
### EulerEquationMod.calcAxidxi

Calculates Axi*Dxi + Aeta*Deta at the element level

**Inputs**

*  `Axidxi` : Product of flux jacobian with shape function derivative
              Axi*Dxi + Aeta*Deta
*  `shapefuncderiv`: Shape function derivative (Dxi and Deta above)
*  `Axi` : Flux jacobian in the xi direction
*  `Aeta` : Flux jacobian in the eta direction
*  `ndof` : Number of degrees of freedom per node
*  `nnpe` : Number of nodes per element

**Outputs**

*  None
"""->
function calcAxidxi(Axidxi::AbstractArray{Tsol, 2},
                    shapefuncderiv::AbstractArray{Tsol,3},
                    Axi::AbstractArray{Tsol,3},
                    Aeta::AbstractArray{Tsol,3}, ndof, nnpe) where Tsol

  for i = 1:nnpe
    for j = 1:nnpe
      m = (i-1)*ndof+1
      n = (j-1)*ndof+1
      nAxi = sview(Axi,:,:,i) # Flux Jacobian at a node in xi direction
      nAeta = sview(Aeta,:,:,i)
      Axidxi[m:(m+ndof-1),n:(n+ndof-1)] = nAxi*shapefuncderiv[i,j,1] +
                                          nAeta*shapefuncderiv[i,j,2]
    end # end for j = 1:nnpe
  end   # end for i = 1:nnpe

  return nothing
end

@doc """
### EulerEquationMod.calcEigenFactorization

Computes the eigen value facotrization of the flux jacobian in either xi or eta
direction. This is a nodal level function. It is only for 2D conservative
variables.

Reference: ``Efficient Solution Methods for the Navier–Stokes Equations``, T.H,
           Pulliam, NASA Ames Research Center

**Inputs**

*  `eqn`  : Euler equation object
*  `q`    : Conservative variables at a node
*  `dxidx`: mapping jacobian vector. It can be either [dξ/dx dξ/dy] or
            [dη/dx dη/dy] depending on what Axi or Aeta being factorized.
*  `T`    : Matrix of eigen vectors.
*  `Tinv` : Inverse of T
*  `Lambda` : Diagonal matrix of eigen values

**Outputs**

*  None

"""->
function calcEigenFactorization(eqn::EulerData{Tsol, Tres, Tdim},
                          q::AbstractArray{Tsol,1}, dxidx::AbstractArray{Tsol,1},
                          T::AbstractArray{Tsol,2}, Tinv::AbstractArray{Tsol,2},
                          Lambda::AbstractArray{Tsol,2}) where {Tsol, Tres, Tdim}

  kappax = dxidx[1]
  kappay = dxidx[2]
  kappa = sqrt(kappax*kappax + kappay*kappay)
  u = q[2]/q[1] # velocity in the x-direction
  v = q[3]/q[1] # velocity in the y-direction
  Temp = (q[4] - 0.5*(q[2]*q[2] + q[3]*q[3])/q[1])*(1/(q[1]*eqn.params.cv)) # Temperature
  c = sqrt(eqn.params.gamma*eqn.params.R*Temp)  # Speed of sound
  phi2 = 0.5*eqn.params.gamma_1*(u*u + v*v) # phi^2
  beta = 0.5/(c*c)
  thetacap = (kappax*u + kappay*v)/kappa
  U = kappax*u + kappay*v

  # Matrix of eigen vectors
  T[1,1] = 1.0
  T[2,1] = u
  T[3,1] = v
  T[4,1] = phi2/eqn.params.gamma_1
  T[1,2] = 0.0
  T[2,2] = kappay/kappa
  T[3,2] = -kappax/kappa
  T[4,2] = (kappay*u - kappax*v)/kappa
  T[1,3] = 1.0
  T[2,3] = u + kappax*c/kappa
  T[3,3] = v + kappay*c/kappa
  T[4,3] = (phi2 + c*c)/eqn.params.gamma_1 + c*thetacap
  T[1,4] = 1.0
  T[2,4] = u - kappax*c/kappa
  T[3,4] = v - kappay*c/kappa
  T[4,4] = (phi2 + c*c)/eqn.params.gamma_1 - c*thetacap

  # Inverse of matrix of eigen vectors
  Tinv[1,1] = 1 - phi2/(c*c)
  Tinv[2,1] = -(kappay*u - kappax*v)/kappa
  Tinv[3,1] = beta*(phi2 - c*thetacap)
  Tinv[4,1] = beta*(phi2 + c*thetacap)
  Tinv[1,2] = eqn.params.gamma_1*u/(c*c)
  Tinv[2,2] = kappay/kappa
  Tinv[3,2] = beta*(kappax*c/kappa - eqn.params.gamma_1*u)
  Tinv[4,2] = -beta*(kappax*c/kappa + eqn.params.gamma_1*u)
  Tinv[1,3] = eqn.params.gamma_1*v/(c*c)
  Tinv[2,3] = -kappax/kappa
  Tinv[3,3] = beta*(kappay*c/kappa - eqn.params.gamma_1*v)
  Tinv[4,3] = -beta*(kappay*c/kappa + eqn.params.gamma_1*v)
  Tinv[1,4] = -eqn.params.gamma_1/(c*c)
  Tinv[2,4] = 0.0
  Tinv[3,4] = beta*eqn.params.gamma_1
  Tinv[4,4] = beta*eqn.params.gamma_1

  # Get the eigen values
  Lambda[1,1] = U
  Lambda[2,2] = U
  Lambda[3,3] = U + c*kappa
  Lambda[4,4] = U - c*kappa

  return nothing
end

@doc """
### EulerEquationMod.calcElementArea

Calculates the element area of a 2D triangular element

**Inputs**

*  `coords` : Vertex coordinates of the element

**Outputs**

*  `area` : Area of the triangular element
"""->
function calcElementArea(coords::AbstractArray{Tmsh, 2}) where Tmsh
  # Calculates the element area using coordinates
  # 2D function for linear mapping

  A = coords[:,1]
  B = coords[:,2]
  C = coords[:,3]
  area = 0.5*(A[1]*(B[2] - C[2]) + B[1]*(C[2] - A[2]) + C[1]*(A[2] - B[2]))

  return area
end # end function calcElementArea

@doc """
### EulerEquationMod.circumcircleDiameter

Calculates the circumcircle diameter of a triangular element

**Inputs**

*  `coords` : Vertex coordinates of the element

**Outputs**

*  `dia` : Diameter of the circumcircle
"""->
function circumcircleDiameter(coords::AbstractArray{Tmsh, 2}) where Tmsh
  # Calculates the circumcircle diameter of the triangular element
  # 2D function with linear mapping
  # Reference: http://geometryatlas.com/entries/109

  x0 = coords[1,1]
  y0 = coords[2,1]
  x1 = coords[1,2]
  y1 = coords[2,2]
  x2 = coords[1,3]
  y2 = coords[2,3]

  a = sqrt((x0-x1)^2 + (y0-y1)^2)
  b = sqrt((x0-x2)^2 + (y0-y2)^2)
  c = sqrt((x1-x2)^2 + (y1-y2)^2)

  dia = a*b*c/sqrt((a+b+c)*(b+c-a)*(c+a-b)*(a+b-c))

  return dia
end  # end function circumcircleDiameter

@doc """
### EulerEquationMod.SUPG

Add Streamwise Upwind Petrov-Galerkin stabilization to the weak residual. This
implementation is only for steady problems

**Inputs**

*  `mesh`: AbstractMesh type
*  `sbp` : Summation-by-parts operator
*  `eqn` : Equation object used elsewhere

**Outputs**

* None

"""->

# SUPG implementation
function SUPG(mesh::AbstractMesh{Tmsh}, sbp::AbstractOperator,
      eqn::EulerData{Tsol, Tres, Tdim}) where {Tmsh, Tsol, Tres, Tdim}

  FluxJacobian(mesh, sbp, eqn) # Calculate the euler flux jacobian
  tau = zeros(Tsol, mesh.numNodesPerElement, mesh.numEl) # Stabilization term
  calcStabilizationTerm(mesh, sbp, eqn, tau)
  # println("tau = \n", tau)

  #=
  # Calculate strong residual
  strong_res = zeros(eqn.res)
  for i = 1:Tdim
    flux_parametric_i = sview(eqn.flux_parametric,:,:,:,i)
    differentiate!(sbp, i, flux_parametric_i, strong_res)
  end
  =#
  supg_res = zeros(eqn.res)
  intvec = zeros(Tsol, mesh.numDofPerNode, mesh.numNodesPerElement, mesh.numEl,
                 Tdim) # intermediate vector for calculating integral

  # SUPG works at the element level so we need to do do a loop over elements
  # get the strong residual from the weak residual. since it also includes the
  # boundary conditions

  for i = 1:mesh.numEl
    for j = 1:mesh.numNodesPerElement
      strong_res = zeros(Tsol, mesh.numDofPerNode)
      JHinverse = mesh.jac[j,i]/sbp.w[j]
      for k = 1:mesh.numDofPerNode
        strong_res[k] = JHinverse*eqn.res[k,j,i]
      end
      Axi = sview(eqn.Axi,:,:,j,i)
      Aeta = sview(eqn.Aeta,:,:,j,i)
      intvec[:,j,i,1] = (tau[j,i]*Axi).'*strong_res # [:,j,i]
      intvec[:,j,i,2] = (tau[j,i]*Aeta).'*strong_res # [:,j,i]
    end # end for j = 1:mesh.numNodesPerElement
  end   # end for i = 1:mesh.numEl

  # calculate the SUPG residual
  for i = 1:Tdim
    intvec_i = sview(intvec,:,:,:,i)
    weakdifferentiate!(sbp, i, intvec_i,supg_res, trans=true)
  end

  #=
  supg_res_vec = zeros(eqn.res_vec)
  for i=1:mesh.numEl  # loop over elements
    for j=1:mesh.numNodesPerElement
      for k=1:4  # loop over dofs on the node
        dofnum_k = mesh.dofs[k, j, i]
        supg_res_vec[dofnum_k] += supg_res[k,j,i]
      end
    end
  end
  =#
  # Add tthe SUPG residual to the weak residual
  for i = 1:mesh.numEl
    for j = 1:mesh.numNodesPerElement
      for k = 1:mesh.numDofPerNode
        eqn.res[k,j,i] += supg_res[k,j,i] # because the negative sign is
      end                                 # already incorporated in the weak
    end                                   # residual.
  end

  # innerprod_supg = eqn.q_vec.'*supg_res_vec
  # println("innerprod_SUPG = ", innerprod_supg)

  #  println("eqn.res = \n", eqn.res)
  #  println("eqn.q = \n", eqn.q)
  return nothing
end # end function SUPG

@doc """
### EulerEquationMod.calcStabilizationTerm

Calculates the stabilization term tau for all the nodes in the mesh

**Inputs**

*  `mesh` : Abstract mesh type
*  `sbp`  : Summation-by-parts operator
*  `eqn`  : Equation object within EulerEquationMod
*  `tau`  : Stabilization term being calculated

**Outputs**

*  None
"""->
# Stabilization Term 3
function calcStabilizationTerm(mesh::AbstractMesh{Tmsh},
       sbp::AbstractOperator, eqn::EulerData{Tsol, Tres, Tdim},
       tau::AbstractArray{Tsol,2}) where {Tmsh, Tsol, Tres, Tdim}

  # Reference: http://enu.kz/repository/2010/AIAA-2010-1183.pdf, eqn 15

  # Get shape function derivatives
  Hinv = 1./sbp.w
  # println(Hinv)

  shapefuncderiv = zeros(sbp.numnodes, sbp.numnodes, Tdim)
  for k = 1:Tdim
    for i = 1:sbp.numnodes
      for j = 1:sbp.numnodes
        shapefuncderiv[j,i,k] = Hinv[i]*sbp.Q[j,i,k]
      end
    end
  end
  # println(eqn.Axi)

  for i = 1:mesh.numEl
    for j = 1:mesh.numNodesPerElement
      q = sview(eqn.q,:,j,i)
      T = (q[4] - 0.5*(q[2]*q[2] + q[3]*q[3])/q[1])*(1/(q[1]*eqn.params.cv))
      c = sqrt(eqn.params.gamma*eqn.params.R*T)  # Speed of sound
      ux = q[2]/q[1]
      uy = q[3]/q[1]
      h_supg = 0.0
      for k = 1:sbp.numnodes
        h_supg += abs(ux*shapefuncderiv[j,k,1] + uy*shapefuncderiv[j,k,2])
      end
      h_supg = 2/h_supg
      tau[j,i] = 0.5*h_supg/(c + sqrt(ux*ux + uy*uy))
    end
  end

  return nothing
end # end calcStabilizationTerm

#=
# Stabilization Term 1
function calcStabilizationTerm(mesh::AbstractMesh{Tmsh},
                               sbp::AbstractOperator, eqn::EulerData{Tsol, Tres, Tdim},
                               tau::AbstractArray{Tsol,2}) where {Tmsh, Tsol, Tres, Tdim}

  # q in the parametric space. Since everything happens in this space
  q_param = zeros(Tsol, mesh.numDofPerNode, mesh.numNodesPerElement, mesh.numEl)

  for i = 1:mesh.numEl
    for j = 1:mesh.numNodesPerElement
      q = sview(eqn.q,:,j,i)
      dxidx = sview(mesh.dxidx,:,:,j,i)
      rhoe = q[4] -0.5*(q[2]*q[2] + q[3]*q[3])/(q[1]*q[1])
      uxi = (q[2]*dxidx[1,1] + q[3]*dxidx[1,2])/q[1]
      ueta = (q[2]*dxidx[2,1] + q[3]*dxidx[2,2])/q[1]
      q_param[1,j,i] = q[1]
      q_param[2,j,i] = q[1]*uxi
      q_param[3,j,i] = q[1]*ueta
      q_param[4,j,i] = rhoe
    end # end for j = 1:mesh.numNodesPerElement
  end   # end for i = 1:mesh.numEl

  beta = zeros(Tsol,2, mesh.numNodesPerElement, mesh.numEl)

  for k = 1:Tdim
    res = zeros(q_param)
    differentiate!(sbp, k, q_param, res)
    for i = 1:mesh.numEl
      for j = 1:mesh.numNodesPerElement
        for l = 1:mesh.numDofPerNode
          beta[k,j,i] += q_param[l,j,i]*res[l,j,i]
        end # end for l = 1:mesh.numDofPerNode
      end   # end for j = 1:mesh.numNodesPerElement
    end     # end for i = 1:mesh.numEl
  end       # end for k = 1:Tdim

  for i = 1:mesh.numEl
    elem_area = calcElementArea(mesh.coords[:,:,i])
    h = sqrt(2*elem_area)
    for j = 1:mesh.numNodesPerElement
      beta[:,j,i] = beta[:,j,i]/norm(beta[:,j,i],2)
      q = sview(eqn.q,:,j,i)
      T = (q[4] - 0.5*(q[2]*q[2] + q[3]*q[3])/q[1])*(1/(q[1]*eqn.params.cv))
      c = sqrt(eqn.params.gamma*eqn.params.R*T)  # Speed of sound
      uxi = zeros(Tsol,2) # Array of velocities in the xi & eta direction
      uxi[1] = q_param[2,j,i]/q_param[1,j,i]
      uxi[2] = q_param[3,j,i]/q_param[1,j,i]
      # Advective stabilization term
      tau_a = 0.5*h/(c + norm(uxi.'*beta[:,j,i], 1))
      # tau[j,i] = 0.005*max(0.0, tau_a)
      tau[j,i] = max(0.0, tau_a)
    end # end for j = 1:mesh.numNodesPerElement
  end   # for i = 1:mesh.numEl

  return nothing
end # end calcStabilizationTerm
=#

#=
# Stabilization term 2
function calcStabilizationTerm(mesh::AbstractMesh{Tmsh},
                               sbp::AbstractOperator, eqn::EulerData{Tsol, Tres, Tdim},
                               tau::AbstractArray{Tsol,2}) where {Tmsh, Tsol, Tres, Tdim}

  # Reference for stabilization: http://enu.kz/repository/2010/AIAA-2010-1183.pdf

  for i = 1:mesh.numEl
    for j = 1:mesh.numNodesPerElement
      C = 2 # for 2D quad with ref. coord [-1,1]
            # for 2D tri with ref coord [0,1], C = 1
      q = sview(eqn.q,:,j,i)
      T = (q[4] - 0.5*(q[2]*q[2] + q[3]*q[3])/q[1])*(1/(q[1]*eqn.params.cv))
      c = sqrt(eqn.params.gamma*eqn.params.R*T)  # Speed of sound
      ux = q[2]/q[1]
      uy = q[3]/q[1]
      dxidx = sview(mesh.dxidx,:,:, j, i)
      if abs(ux) < 1e-13 || abs(uy) < 1e-13
        tau[j,i] = 0.0
      else
        g_ij = dxidx[1,1]*dxidx[1,2] + dxidx[2,1]*dxidx[2,2]
        if abs(g_ij) < 1e-13
          tau[j,i] = 0.0
        else
          elem_area = calcElementArea(mesh.coords[:,:,i])
          # h_supg = sqrt(2*elem_area)
          h_supg = C*sqrt((ux*ux + uy*uy)/(ux*g_ij*uy))
          tau[j,i] = 0.5*h_supg/(c + sqrt(ux*ux + uy*uy))
        end # end if
      end   # end if
    end # end for j = 1:mesh.numNodesPerElement
  end   # end for i = 1:mesh.numEl

  return nothing
end # end calcStabilizationTerm
=#
#=
# Stabilization term 4
function calcStabilizationTerm(mesh::AbstractMesh{Tmsh},
                               sbp::AbstractOperator, eqn::EulerData{Tsol, Tres, Tdim},
                               tau::AbstractArray{Tsol,2}) where {Tmsh, Tsol, Tres, Tdim}

  # Reference: Three-Dimensional Stabilized Finite Elements for Compressible
  #            Navier–Stokes, T. Taylor Erwin, AIAA Journal Vol 51, No. 6,
  #            June 2013

  # Get shape function derivatives
  Hinv = 1./sbp.w
  # println(Hinv)

  shapefuncderiv = zeros(sbp.numnodes, sbp.numnodes, Tdim)
  for k = 1:Tdim
    for i = 1:sbp.numnodes
      for j = 1:sbp.numnodes
        shapefuncderiv[j,i,k] = Hinv[i]*sbp.Q[j,i,k]
      end
    end
  end

  for i = 1:mesh.numEl
    for j = 1:mesh.numNodesPerElement
      Axi = sview(eqn.Axi,:,:,j,i)
      Aeta = sview(eqn.Aeta,:,:,j,i)
      invtau = zeros(mesh.numDofPerNode, mesh.numDofPerNode)
      for k = 1:sbp.numnodes
        invtau += shapefuncderiv[j,k,1]*Axi + shapefuncderiv[j,k,2]*Aeta
        #println("invtau = ", invtau)
      end
      T = inv(invtau)
      println("T = \n", T)
    end
  end

  return nothing
end # end calcStabilizationTerm
=#

#------------------------------------------------------------------------------
# Debugging code
# calculate the boundary integral using the actual euler flux

function getPhysBCFluxes(mesh, sbp, eqn, opts, bndryfluxPhysical)
  # Calculate the physical BC flux

  #println("mesh.bndry_funcs = ", mesh.bndry_funcs)

  functor_i = BCDict["isentropicVortexBC_physical"]
  for i=1:mesh.numBC
  #  println("computing flux for boundary condition ", i)
    # functor_i = mesh.bndry_funcs[i]
    start_index = mesh.bndry_offsets[i]
    end_index = mesh.bndry_offsets[i+1]
    bndry_facenums_i = sview(mesh.bndryfaces, start_index:(end_index - 1))
    bndryflux_i = sview(bndryfluxPhysical, :, :, start_index:(end_index - 1))
    #functor_i(q, flux_parametric, aux_vars, x, dxidx, nrm, bndryflux_i, eqn.params)

    calcBoundaryFlux(mesh, sbp, eqn, functor_i, bndry_facenums_i, bndryflux_i)
  end

  return nothing
end

function residualComparison(mesh, sbp, eqn, opts)

 # Get strong residual by differentiating
  differentiation_strong_res = zeros(eqn.res)

  for i = 1:Tdim
    flux_parametric_i = sview(eqn.flux_parametric,:,:,:,i)
    differentiate!(sbp, i, flux_parametric_i, differentiation_strong_res)
  end

  # Get strong residual from the weak form
  strong_res_from_weak = zeros(eqn.res)
  for i = 1:mesh.numEl
    for j = 1:mesh.numNodesPerElement
      JHinverse = mesh.jac[j,i]/sbp.w[j]
      for k = 1:mesh.numDofPerNode
        strong_res_from_weak[k,j,i] = JHinverse*eqn.res[k,j,i]
      end # end for k = 1:mesh.numDofPerNode
    end # end for i = 1:mesh.numEl
  end   # end for j = 1:mesh.numNodesPerElement

  println("differentiation_strong_res = \n", differentiation_strong_res)
  println("strong_res_from_weak = \n", strong_res_from_weak)
  # println("\n eqn.res = \n", eqn.res)
  Error = zeros(eqn.res)
  for i = 1:mesh.numEl
    Error[:,:,i] = differentiation_strong_res[:,:,i] + strong_res_from_weak[:,:,i]
  end
  println("\nError = \n", Error)

  return nothing
end # end function residualComparison

#----------------------------------------------------------------------------


text:
module CURAND

export create_generator,
       destroy_generator,
       # generator types
       CURAND_RNG_TEST,
       CURAND_RNG_PSEUDO_DEFAULT,
       CURAND_RNG_PSEUDO_XORWOW,
       CURAND_RNG_PSEUDO_MRG32K3A,
       CURAND_RNG_PSEUDO_MTGP32,
       CURAND_RNG_PSEUDO_MT19937,
       CURAND_RNG_PSEUDO_PHILOX4_32_10,
       CURAND_RNG_QUASI_DEFAULT,
       CURAND_RNG_QUASI_SOBOL32,
       CURAND_RNG_QUASI_SCRAMBLED_SOBOL32,
       CURAND_RNG_QUASI_SOBOL64,
       CURAND_RNG_QUASI_SCRAMBLED_SOBOL64,
       # high-level API
       curand,
       curandn,
       curand_logn,
       curand_poisson,
       # wrappers
       get_version,
       set_pseudo_random_generator_seed,
       set_generator_offset,
       set_generator_ordering,
       set_quasi_random_generator_dimensions,
       generate,
       generate_long_long,
       generate_uniform,
       generate_uniform_double,
       generate_normal,
       generate_normal_double,
       generate_log_normal,
       generate_log_normal_double,
       create_poisson_distribtion,
       destroy_distribution,
       generate_poisson,
       generate_seeds

include("core.jl")

end # module


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

using Base: @deprecate, depwarn

# BEGIN 0.7 deprecations

# PR #22475
import Base: cat
@deprecate cat(::Type{Val{N}}, A::_SparseConcatGroup...) where {N} cat(Val(N), A...)
@deprecate cat(::Type{Val{N}}, A::_DenseConcatGroup...) where {N} cat(Val(N), A...)

# deprecate remaining vectorized methods over SparseVectors (zero-preserving)
for op in (:floor, :ceil, :trunc, :round,
        :log1p, :expm1,  :sinpi,
        :sin,   :tan,    :sind,   :tand,
        :asin,  :atan,   :asind,  :atand,
        :sinh,  :tanh,   :asinh,  :atanh)
    @eval import Base.Math: $op
    @eval @deprecate ($op)(x::AbstractSparseVector{<:Number,<:Integer}) ($op).(x)
end
# deprecate remaining vectorized methods over SparseVectors (not-zero-preserving)
for op in (:exp, :exp2, :exp10, :log, :log2, :log10,
           :cos, :cosd, :acos, :cosh, :cospi,
           :csc, :cscd, :acot, :csch, :acsch,
           :cot, :cotd, :acosd, :coth,
           :sec, :secd, :acotd, :sech, :asech)
    @eval import Base.Math: $op
    @eval @deprecate ($op)(x::AbstractSparseVector{<:Number,<:Integer}) ($op).(x)
end

# PR 23341
import LinearAlgebra: diagm
@deprecate diagm(A::SparseMatrixCSC) sparse(Diagonal(sparsevec(A)))

# PR #23757
@deprecate spdiagm(x::AbstractVector) sparse(Diagonal(x))
function spdiagm(x::AbstractVector, d::Number)
    depwarn(string("`spdiagm(x::AbstractVector, d::Number)` is deprecated, use ",
        "`spdiagm(d => x)` instead, which now returns a square matrix. To preserve the old ",
        "behaviour, use `sparse(SparseArrays.spdiagm_internal(d => x)...)`"), :spdiagm)
    I, J, V = spdiagm_internal(d => x)
    return sparse(I, J, V)
end
function spdiagm(x, d)
    depwarn(string("`spdiagm((x1, x2, ...), (d1, d2, ...))` is deprecated, use ",
        "`spdiagm(d1 => x1, d2 => x2, ...)` instead, which now returns a square matrix. ",
        "To preserve the old behaviour, use ",
        "`sparse(SparseArrays.spdiagm_internal(d1 => x1, d2 => x2, ...)...)`"), :spdiagm)
    I, J, V = spdiagm_internal((d[i] => x[i] for i in 1:length(x))...)
    return sparse(I, J, V)
end
function spdiagm(x, d, m::Integer, n::Integer)
    depwarn(string("`spdiagm((x1, x2, ...), (d1, d2, ...), m, n)` is deprecated, use ",
        "`spdiagm(d1 => x1, d2 => x2, ...)` instead, which now returns a square matrix. ",
        "To specify a non-square matrix and preserve the old behaviour, use ",
        "`I, J, V = SparseArrays.spdiagm_internal(d1 => x1, d2 => x2, ...); sparse(I, J, V, m, n)`"), :spdiagm)
    I, J, V = spdiagm_internal((d[i] => x[i] for i in 1:length(x))...)
    return sparse(I, J, V, m, n)
end

@deprecate sparse(s::UniformScaling, m::Integer) sparse(s, m, m)

# PR #25037
@deprecate spones(A::SparseMatrixCSC) LinearAlgebra.fillstored!(copy(A), 1)
@deprecate spones(A::SparseVector)    LinearAlgebra.fillstored!(copy(A), 1)
export spones

# full for sparse arrays
import Base: full
function full(S::Union{SparseVector,SparseMatrixCSC})
    (arrtypestr, desttypestr) =
        isa(S, SparseVector)    ? ("SparseVector",    "Vector") :
        isa(S, SparseMatrixCSC) ? ("SparseMatrixCSC", "Matrix") :
            error("should not be reachable!")
    depwarn(string(
        "`full(S::$(arrtypestr))` (and `full` in general) has been deprecated. ",
        "To replace `full(S::$(arrtypestr))`, consider `$(desttypestr)(S)` or, ",
        "if that option is too narrow, `Array(S)`."), :full)
    return Array(S)
end

# issue #22849
import Base: reinterpret
@deprecate reinterpret(::Type{T}, a::SparseMatrixCSC{S}, dims::NTuple{N,Int}) where {T, S, N} reinterpret(T, reshape(a, dims))

# deprecate speye
export speye
function speye(n::Integer)
    depwarn(string("`speye(n::Integer)` has been deprecated in favor of `I`, `sparse`, and ",
                    "`SparseMatrixCSC` constructor methods. For a direct replacement, consider ",
                    "`sparse(1.0I, n, n)`, `SparseMatrixCSC(1.0I, n, n)`, or `SparseMatrixCSC{Float64}(I, n, n)`. ",
                    "If `Float64` element type is not necessary, consider the shorter `sparse(I, n, n)` ",
                    "or `SparseMatrixCSC(I, n, n)` (with default `eltype(I)` of `Bool`)."), :speye)
    return sparse(1.0I, n, n)
end
function speye(m::Integer, n::Integer)
    depwarn(string("`speye(m::Integer, n::Integer)` has been deprecated in favor of `I`, ",
                    "`sparse`, and `SparseMatrixCSC` constructor methods. For a direct ",
                    "replacement, consider `sparse(1.0I, m, n)`, `SparseMatrixCSC(1.0I, m, n)`, ",
                    "or `SparseMatrixCSC{Float64}(I, m, n)`. If `Float64` element type is not ",
                    " necessary, consider the shorter `sparse(I, m, n)` or `SparseMatrixCSC(I, m, n)` ",
                    "(with default `eltype(I)` of `Bool`)."), :speye)
    return sparse(1.0I, m, n)
end
function speye(::Type{T}, n::Integer) where T
    depwarn(string("`speye(T, n::Integer)` has been deprecated in favor of `I`, `sparse`, and ",
                    "`SparseMatrixCSC` constructor methods. For a direct replacement, consider ",
                    "`sparse(T(1)I, n, n)` if `T` is concrete or `SparseMatrixCSC{T}(I, n, n)` ",
                    "if `T` is either concrete or abstract. If element type `T` is not necessary, ",
                    "consider the shorter `sparse(I, n, n)` or `SparseMatrixCSC(I, n, n)` ",
                    "(with default `eltype(I)` of `Bool`)."), :speye)
    return SparseMatrixCSC{T}(I, n, n)
end
function speye(::Type{T}, m::Integer, n::Integer) where T
    depwarn(string("`speye(T, m::Integer, n::Integer)` has been deprecated in favor of `I`, ",
                    "`sparse`, and `SparseMatrixCSC` constructor methods. For a direct ",
                    "replacement, consider `sparse(T(1)I, m, n)` if `T` is concrete or ",
                    "`SparseMatrixCSC{T}(I, m, n)` if `T` is either concrete or abstract. ",
                    "If element type `T` is not necessary, consider the shorter ",
                    "`sparse(I, m, n)` or `SparseMatrixCSC(I, m, n)` (with default `eltype(I)` ",
                    "of `Bool`)."), :speye)
    return SparseMatrixCSC{T}(I, m, n)
end
function speye(S::SparseMatrixCSC{T}) where T
    depwarn(string("`speye(S::SparseMatrixCSC{T})` has been deprecated in favor of `I`, ",
                    "`sparse`, and `SparseMatrixCSC` constructor methods. For a direct ",
                    "replacement, consider `sparse(T(1)I, size(S)...)` if `T` is concrete or ",
                    "`SparseMatrixCSC{eltype(S)}(I, size(S))` if `T` is either concrete or abstract. ",
                    "If preserving element type `T` is not necessary, consider the shorter ",
                    "`sparse(I, size(S)...)` or `SparseMatrixCSC(I, size(S))` (with default ",
                    "`eltype(I)` of `Bool`)."), :speye)
    return SparseMatrixCSC{T}(I, size(S)...)
end

# former imports into SparseArrays
import Base: Ac_mul_B, At_mul_B
import Base: A_mul_Bc, A_mul_Bt, Ac_mul_Bc, At_mul_Bt
import Base: At_ldiv_B, Ac_ldiv_B
import LinearAlgebra: A_mul_B!, Ac_mul_B!, At_mul_B!, A_ldiv_B!
import LinearAlgebra: At_ldiv_B!, Ac_ldiv_B!, A_rdiv_B!, A_rdiv_Bc!, mul!, ldiv!, rdiv!

# A[ct]_(mul|ldiv|rdiv)_B[ct][!] methods from base/sparse/linalg.jl, to deprecate
using LinearAlgebra: Adjoint, Transpose
@deprecate Ac_ldiv_B(A::SparseMatrixCSC, B::RowVector)  (\)(adjoint(A), B)
@deprecate At_ldiv_B(A::SparseMatrixCSC, B::RowVector)  (\)(transpose(A), B)
@deprecate Ac_ldiv_B(A::SparseMatrixCSC, B::AbstractVecOrMat)   (\)(adjoint(A), B)
@deprecate At_ldiv_B(A::SparseMatrixCSC, B::AbstractVecOrMat)   (\)(transpose(A), B)
@deprecate A_rdiv_Bc!(A::SparseMatrixCSC{T}, D::Diagonal{T}) where {T}  rdiv!(A, adjoint(D))
@deprecate A_rdiv_Bt!(A::SparseMatrixCSC{T}, D::Diagonal{T}) where {T}  rdiv!(A, transpose(D))
@deprecate A_rdiv_B!(A::SparseMatrixCSC{T}, D::Diagonal{T}) where {T}   rdiv!(A, D)
@deprecate A_ldiv_B!(L::LowerTriangular{T,<:SparseMatrixCSCUnion{T}}, B::StridedVecOrMat) where {T}     ldiv!(L, B)
@deprecate A_ldiv_B!(U::UpperTriangular{T,<:SparseMatrixCSCUnion{T}}, B::StridedVecOrMat) where {T}     ldiv!(U, B)
@deprecate A_mul_Bt(A::SparseMatrixCSC{Tv,Ti}, B::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti}     (*)(A, transpose(B))
@deprecate A_mul_Bc(A::SparseMatrixCSC{Tv,Ti}, B::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti}     (*)(A, adjoint(B))
@deprecate At_mul_B(A::SparseMatrixCSC{Tv,Ti}, B::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti}     (*)(transpose(A), B)
@deprecate Ac_mul_B(A::SparseMatrixCSC{Tv,Ti}, B::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti}     (*)(adjoint(A), B)
@deprecate At_mul_Bt(A::SparseMatrixCSC{Tv,Ti}, B::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti}    (*)(transpose(A), transpose(B))
@deprecate Ac_mul_Bc(A::SparseMatrixCSC{Tv,Ti}, B::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti}    (*)(adjoint(A), adjoint(B))
@deprecate A_mul_B!(C::StridedVecOrMat, A::SparseMatrixCSC, B::StridedVecOrMat)     mul!(C, A, B)
@deprecate Ac_mul_B!(C::StridedVecOrMat, A::SparseMatrixCSC, B::StridedVecOrMat)    mul!(C, adjoint(A), B)
@deprecate At_mul_B!(C::StridedVecOrMat, A::SparseMatrixCSC, B::StridedVecOrMat)    mul!(C, transpose(A), B)
@deprecate A_mul_B!(α::Number, A::SparseMatrixCSC, B::StridedVecOrMat, β::Number, C::StridedVecOrMat)   mul!(C, A, B, α, β)
@deprecate A_mul_B(A::SparseMatrixCSC{TA,S}, x::StridedVector{Tx}) where {TA,S,Tx}  (*)(A, x)
@deprecate A_mul_B(A::SparseMatrixCSC{TA,S}, B::StridedMatrix{Tx}) where {TA,S,Tx}  (*)(A, B)
@deprecate Ac_mul_B!(α::Number, A::SparseMatrixCSC, B::StridedVecOrMat, β::Number, C::StridedVecOrMat)  mul!(C, adjoint(A), B, α, β)
@deprecate Ac_mul_B(A::SparseMatrixCSC{TA,S}, x::StridedVector{Tx}) where {TA,S,Tx}     (*)(adjoint(A), x)
@deprecate Ac_mul_B(A::SparseMatrixCSC{TA,S}, B::StridedMatrix{Tx}) where {TA,S,Tx}     (*)(adjoint(A), B)
@deprecate At_mul_B!(α::Number, A::SparseMatrixCSC, B::StridedVecOrMat, β::Number, C::StridedVecOrMat)  mul!(C, transpose(A), B, α, β)
@deprecate At_mul_B(A::SparseMatrixCSC{TA,S}, x::StridedVector{Tx}) where {TA,S,Tx}     (*)(transpose(A), x)
@deprecate At_mul_B(A::SparseMatrixCSC{TA,S}, B::StridedMatrix{Tx}) where {TA,S,Tx}     (*)(transpose(A), B)
@deprecate A_mul_Bt(A::SparseMatrixCSC{TvA,TiA}, B::SparseMatrixCSC{TvB,TiB}) where {TvA,TiA,TvB,TiB}   (*)(A, transpose(B))
@deprecate A_mul_Bc(A::SparseMatrixCSC{TvA,TiA}, B::SparseMatrixCSC{TvB,TiB}) where {TvA,TiA,TvB,TiB}   (*)(A, adjoint(B))
@deprecate At_mul_B(A::SparseMatrixCSC{TvA,TiA}, B::SparseMatrixCSC{TvB,TiB}) where {TvA,TiA,TvB,TiB}   (*)(transpose(A), B)
@deprecate Ac_mul_B(A::SparseMatrixCSC{TvA,TiA}, B::SparseMatrixCSC{TvB,TiB}) where {TvA,TiA,TvB,TiB}   (*)(adjoint(A),B)
@deprecate At_mul_Bt(A::SparseMatrixCSC{TvA,TiA}, B::SparseMatrixCSC{TvB,TiB}) where {TvA,TiA,TvB,TiB}  (*)(transpose(A), transpose(B))
@deprecate Ac_mul_Bc(A::SparseMatrixCSC{TvA,TiA}, B::SparseMatrixCSC{TvB,TiB}) where {TvA,TiA,TvB,TiB}  (*)(adjoint(A), adjoint(B))

# A[ct]_(mul|ldiv|rdiv)_B[ct][!] methods from base/sparse/sparsevector.jl, to deprecate
for isunittri in (true, false), islowertri in (true, false)
    unitstr = isunittri ? "Unit" : ""
    halfstr = islowertri ? "Lower" : "Upper"
    tritype = :(LinearAlgebra.$(Symbol(unitstr, halfstr, "Triangular")))
    @eval #=Base.SparseArrays=# begin
        using LinearAlgebra: Adjoint, Transpose
        @deprecate At_ldiv_B(A::$tritype{TA,<:AbstractMatrix}, b::SparseVector{Tb}) where {TA<:Number,Tb<:Number}   (\)(transpose(A), b)
        @deprecate At_ldiv_B(A::$tritype{TA,<:StridedMatrix}, b::SparseVector{Tb}) where {TA<:Number,Tb<:Number}    (\)(transpose(A), b)
        @deprecate At_ldiv_B(A::$tritype, b::SparseVector)  (\)(transpose(A), b)
        @deprecate Ac_ldiv_B(A::$tritype{TA,<:AbstractMatrix}, b::SparseVector{Tb}) where {TA<:Number,Tb<:Number}   (\)(adjoint(A), b)
        @deprecate Ac_ldiv_B(A::$tritype{TA,<:StridedMatrix}, b::SparseVector{Tb}) where {TA<:Number,Tb<:Number}    (\)(adjoint(A), b)
        @deprecate Ac_ldiv_B(A::$tritype, b::SparseVector)  (\)(adjoint(A), b)
        @deprecate A_ldiv_B!(A::$tritype{<:Any,<:StridedMatrix}, b::SparseVector)   ldiv!(A, b)
        @deprecate At_ldiv_B!(A::$tritype{<:Any,<:StridedMatrix}, b::SparseVector)  ldiv!(transpose(A), b)
        @deprecate Ac_ldiv_B!(A::$tritype{<:Any,<:StridedMatrix}, b::SparseVector)  ldiv!(adjoint(A), b)
    end
end

using LinearAlgebra: Adjoint, Transpose
@deprecate Ac_mul_B(A::SparseMatrixCSC, x::AbstractSparseVector)    (*)(adjoint(A), x)
@deprecate At_mul_B(A::SparseMatrixCSC, x::AbstractSparseVector)    (*)(transpose(A), x)
@deprecate Ac_mul_B!(α::Number, A::SparseMatrixCSC, x::AbstractSparseVector, β::Number, y::StridedVector)   mul!(y, adjoint(A), x, α, β)
@deprecate Ac_mul_B!(y::StridedVector{Ty}, A::SparseMatrixCSC, x::AbstractSparseVector{Tx}) where {Tx,Ty}   mul!(y, adjoint(A), x)
@deprecate At_mul_B!(α::Number, A::SparseMatrixCSC, x::AbstractSparseVector, β::Number, y::StridedVector)   mul!(y, transpose(A), x, α, β)
@deprecate At_mul_B!(y::StridedVector{Ty}, A::SparseMatrixCSC, x::AbstractSparseVector{Tx}) where {Tx,Ty}   mul!(y, transpose(A), x)
@deprecate A_mul_B!(α::Number, A::SparseMatrixCSC, x::AbstractSparseVector, β::Number, y::StridedVector)    mul!(y, A, x, α, β)
@deprecate A_mul_B!(y::StridedVector{Ty}, A::SparseMatrixCSC, x::AbstractSparseVector{Tx}) where {Tx,Ty}    mul!(y, A, x)
@deprecate At_mul_B!(α::Number, A::StridedMatrix, x::AbstractSparseVector, β::Number, y::StridedVector)     mul!(y, transpose(A), x, α, β)
@deprecate At_mul_B!(y::StridedVector{Ty}, A::StridedMatrix, x::AbstractSparseVector{Tx}) where {Tx,Ty}     mul!(y, transpose(A), x)
@deprecate At_mul_B(A::StridedMatrix{Ta}, x::AbstractSparseVector{Tx}) where {Ta,Tx}    (*)(transpose(A), x)
@deprecate A_mul_B!(α::Number, A::StridedMatrix, x::AbstractSparseVector, β::Number, y::StridedVector)  mul!(y, A, x, α, β)
@deprecate A_mul_B!(y::StridedVector{Ty}, A::StridedMatrix, x::AbstractSparseVector{Tx}) where {Tx,Ty}  mul!(y, A, x)

# methods involving RowVector from base/sparse/linalg.jl, to deprecate
\(::SparseMatrixCSC, ::RowVector) = throw(DimensionMismatch("Cannot left-divide matrix by transposed vector"))
\(::Adjoint{<:Any,<:SparseMatrixCSC}, ::RowVector) = throw(DimensionMismatch("Cannot left-divide matrix by transposed vector"))
\(::Transpose{<:Any,<:SparseMatrixCSC}, ::RowVector) = throw(DimensionMismatch("Cannot left-divide matrix by transposed vector"))

# methods involving RowVector from base/sparse/higherorderfns.jl, to deprecate
@eval SparseArrays.HigherOrderFns begin
    BroadcastStyle(::Type{<:RowVector{T,<:Vector}}) where T = Broadcast.MatrixStyle()
end

import Base: asyncmap
@deprecate asyncmap(f, s::AbstractSparseArray...; kwargs...) sparse(asyncmap(f, map(Array, s)...; kwargs...))

# PR 26347: implicit scalar broadcasting within setindex!
@deprecate setindex!(A::SparseMatrixCSC, x::Number, i::Integer, J::AbstractVector{<:Integer}) (A[i, J] .= x; A)
@deprecate setindex!(A::SparseMatrixCSC, x::Number, I::AbstractVector{<:Integer}, j::Integer) (A[I, j] .= x; A)
@deprecate setindex!(A::SparseMatrixCSC, x, ::Colon)          fill!(A, x)
@deprecate setindex!(A::SparseMatrixCSC, x, ::Colon, ::Colon) fill!(A, x)
@deprecate setindex!(A::SparseMatrixCSC, x, ::Colon, j::Union{Integer, AbstractVector}) (A[:, j] .= x; A)
@deprecate setindex!(A::SparseMatrixCSC, x, i::Union{Integer, AbstractVector}, ::Colon) (A[i, :] .= x; A)
@deprecate setindex!(A::SparseMatrixCSC, x::Number, I::AbstractVector{<:Integer}, J::AbstractVector{<:Integer}) (A[I, J] .= x; A)

#25395 keywords unlocked
@deprecate dropzeros(x, trim)     dropzeros(x, trim = trim)
@deprecate dropzeros!(x, trim)    dropzeros!(x, trim = trim)
@deprecate droptol!(A, tol, trim) droptol!(A, tol, trim = trim)

Base.@deprecate_binding blkdiag blockdiag

@deprecate complex(x::AbstractSparseVector{<:Real}, y::AbstractSparseVector{<:Real}) complex.(x, y)
@deprecate complex(x::AbstractVector{<:Real}, y::AbstractSparseVector{<:Real}) complex.(x, y)
@deprecate complex(x::AbstractSparseVector{<:Real}, y::AbstractVector{<:Real}) complex.(x, y)

@deprecate diff(a::SparseMatrixCSC, dim::Integer) diff(a, dims=dim)

# END 0.7 deprecations

# BEGIN 1.0 deprecations

# END 1.0 deprecations


text:
using Test
using OmegaCore
using Distributions
using OmegaDistributions
using OmegaTest
using OmegaCore.Interventions

function test_merge_1()
    xx = 1 ~ Normal(0, 1)
    y(ω) = xx(ω) + 10
    yi = intervene(y, xx => (ω::Ω -> 200.0))
    yi2 = intervene(yi, xx => (ω::Ω -> 300.0))
    @test randsample(yi2) == 210
  # yi2(ω)
    @test isinferred(randsample, yi2)
end

function test_merge_2()
    xx = 1 ~ Normal(0, 1)
    y(ω) = xx(ω) + 10
    xr = 2 ~ Normal(30, 1)
    yi = intervene(y, xx => xr)
    yi2 = intervene(yi, xr => (ω::Ω -> 300.0))
    @test randsample(yi2) == 310
    @test isinferred(randsample, yi2)
end

function test_merge_3()
    xx = 1 ~ Normal(0, 1)
    y(ω) = xx(ω) + 10
    yi = intervene(y, (xx => (ω::Ω -> 200.0), xx => (ω::Ω -> 300.0)))
    yi3 = intervene(yi, xx => (ω::Ω -> 400.0))
    @test randsample(yi3) == 210
    @test isinferred(randsample, yi3)
end

function test_merge_4()
    xx = 1 ~ Normal(0, 1)
    y(ω) = xx(ω) + 10
    xr = 2 ~ Normal(30, 1)
    xrr = xr .* xr
    yi = intervene(y, xx => xrr)
    yi2 = intervene(yi, xr => (ω::Ω -> 30.0))
    @test randsample(yi2) == 910.0
    @test isinferred(randsample, yi2)
end

function test_merge_more_than_5_interventions()
    xx = 1 ~ Normal(0, 1)
    y(ω) = xx(ω) + 10
    yi = intervene(y, xx => (ω::Ω -> 200.0))
    yi2 = intervene(yi, xx => (ω::Ω -> 300.0))
    yi3 = intervene(yi2, xx => (ω::Ω -> 400.0))
    yi4 = intervene(yi3, xx => (ω::Ω -> 500.0))
    yi5 = intervene(yi4, xx => (ω::Ω -> 600.0))
    yi6 = intervene(yi5, xx => (ω::Ω -> 700.0))
    @test randsample(yi6) == 210
  # isinferred fails
end

function minimal_examplea()
    xx = 1 ~ Normal(0, 1)
    y(ω) = xx(ω) + 10
    yi = intervene(y, xx => (ω::Ω -> 200.0))
    yi2 = intervene(yi, xx => (ω::Ω -> 300.0))
    @test randsample(yi2) == 210
end

function minimal_exampleb()
    xx = 1 ~ Normal(0, 1)
    y(ω) = xx(ω) + 10
    xr = 2 ~ Normal(30, 1)
    yi = intervene(y, xx => xr)
    yi2 = intervene(yi, xr => (ω::Ω -> 300.0))
    @test randsample(yi2) == 310
end

function test_model()
  # Normally distributed random variable with id 1
    x = 1 ~ Normal(0, 1)

    # Normally distributed random variable with id 2 and x as mean
    function y(ω)
        xt = x(ω)
        u = 2 ~ Uniform(xt, xt + 1)
        u(ω)
    end
    x_ = 0.1
    y_ = 0.3

    # model -- tuple-valued random variable ω::Ω -> (x(ω), y(ω)), becase we want joint pdf
    m(ω) = (x(ω), y(ω))
    (x, y, m)
end

function samplemodel()
    x, y, z = test_model()
    ω = defω()
    y(ω)
end

function test_intervention()
    x, y, m = test_model()
    yⁱ = y |ᵈ (x => (ω::Ω -> 100.0))
    @test 100.0 <= randsample(yⁱ) <= 101.0
    @test isinferred(randsample, yⁱ)
end

function test_intervene_diff_parents()
    x = 1 ~ Normal(0, 1)
    y = 2 ~ Normal.(x, 1)
    x2 = 3 ~ Normal(0, 1)
    yi = y |ᵈ (x => ω::Ω -> 100.0)
    yi2 = y |ᵈ (x2 => ω::Ω -> 100.0)
    yi_, yi2_ = randsample((yi, yi2))
    @test yi_ != yi2_
end

function test_two_interventions()
    x = 1 ~ Normal(0, 1)
    y = 2 ~ Uniform(10.0, 20.0)
    z = 3 ~ Normal.(x, y)
    zi = z |ᵈ (x => (ω::Ω -> 100.0), y => (ω::Ω -> 0.1))
    @test 99 <= randsample(zi) <= 101
    @test isinferred(randsample, zi)
end

function test_three_interventions()
    x = 1 ~ Normal(0, 1)
    y = 2 ~ Uniform(10.0, 20.0)
    c = 3 ~ Uniform(2.0, 3.0)
    z = 4 ~ Normal.(x .* c, y)
    zi = z |ᵈ (x => (ω::Ω -> 100.0), y => (ω::Ω -> 0.1), c => (ω::Ω -> 1.0))
    @test 99 <= randsample(zi) <= 101
end

# function test_intervention_logpdf()
#     x, y, m = test_model()
#   # Log density of model wrt ω
#     l = logpdf(m, ω)

#   # Check it is what it should be
#     @test l == logpdf(Normal(0, 1), x_) + logpdf(Normal(x_, 1), y_)

#   # Intervened model
#     v_ = 100.0

#   # y had x been v_
#     yⁱ = y | had(x => v_)

#   # new model with Intervened variables
#     mⁱ = rt(x, yⁱ)

#   # log pdf of interved model on same ω
#     lⁱ = logpdf(mⁱ, ω)
#     logpdf(x, ω)

#     @test lⁱ == logpdf(Normal(0, 1), x_) + logpdf(Normal(v_, 1), y_)
#     @test lⁱ < l
# end

function test_self_intervene()
    p = 0.7
    q = 0.3
    E = ~ 1 ~ Bernoulli(p)     # Execution order
    C = ~ 2 ~ Uniform(0, 1)    # Calmness
    N = C .< q                  # Nerves
    A = E .| N                 # A shoots
    B = E                      # B shoots on order
    D = A .| B                 # Prisoner Dies
    cf = (D |ᵈ (A => 0)) |ᶜ D
    r = 0.8
    na4 = D |ᵈ (A => false, B => ifelse.(3 ~ Bernoulli(r), false, B))
    ω = defω()
    na4(ω)
end

function test_many_intervene_tuple()
  x = 1 ~ Normal(0, 1)
  σ = 2 ~ Uniform(0, 0.1)
  z =3 ~  Normal(x, σ)
  @test isinferred(randsample, intervene(z, (x, σ), (100.0, 20.0)))
end

@testset "intervene" begin
    test_intervention()
    test_intervene_diff_parents()
    test_two_interventions()
    test_three_interventions()
    test_merge_1()
    test_merge_2()
    test_merge_3()
    test_merge_4()
    # test_changed_rettype_merge()
    test_merge_more_than_5_interventions()
    test_self_intervene()
    minimal_examplea()
    minimal_exampleb()
end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
using LMCLUS
using Distributions
using Test
using Random

@testset "Thresholding" begin

	Random.seed!(283739285);
	N = 10000
	bins = 100
	generate_sample(N, P1, μ1, σ1, P2, μ2, σ2) = vcat(rand(Normal(μ1,σ1), round(Int, N*P1)),rand(Normal(μ2,σ2), round(Int, N*P2)))

	# Following tests' parameters are taken from origin work:
	# J. Kittler & J. Illingworth: "Minimum Error Thresholding" Pattern Recognition, Vol 19, nr 1. 1986, pp. 41-47.
	# See fig.1 and fih.2

	# Test 1
	D = generate_sample(N, 0.5, 50, 15, 0.5, 150, 15)
	mindint, maxdist = extrema(D)
	res = separation(LMCLUS.Kittler, D, bins=bins)
	@test threshold(res) ≈ 100.0 atol=2.0
	@test res.mindist == mindint
	@test res.maxdist == maxdist
	@test res.bins == bins

	# Test 2
	res = separation(LMCLUS.Kittler, generate_sample(N, 0.25, 38, 9, 0.75, 121, 44), bins=bins)
	@test threshold(res) ≈ 63.0 atol=2.0

	# Test 3
	res = separation(LMCLUS.Kittler, generate_sample(N, 0.45, 47, 13, 0.55, 144, 25), bins=bins)
	@test threshold(res) ≈ 82.0 atol=2.0

	# Test 4
	res = separation(LMCLUS.Kittler, generate_sample(N, 0.5, 50, 4, 0.5, 150, 30), bins=bins)
	@test threshold(res) ≈ 65.0 atol=2.

	# Try unimodal histogram
	D = rand(Normal(1, 10), N)
	mindint, maxdist = extrema(D)
	res = separation(LMCLUS.Kittler, D, bins=bins)
	@test criteria(res) ≈ 0.0 atol=1e-2
	@test threshold(res) == maxdist # threshold is set to maximal element
	@test res.mindist == mindint
	@test res.maxdist == maxdist
	@test res.bins == bins

	## Otsu ##
	# Test 1
	res = separation(LMCLUS.Otsu, generate_sample(N, 0.5, 50, 15, 0.5, 150, 15), bins=bins)
	@test threshold(res) ≈ 100.0 atol=5.0

	# Test 2
	res = separation(LMCLUS.Otsu, generate_sample(N, 0.25, 38, 9, 0.75, 121, 44), bins=bins)
	@test threshold(res) ≈ 108.0 atol=2.0

	# Test 3
	res = separation(LMCLUS.Otsu, generate_sample(N, 0.45, 47, 13, 0.55, 144, 25), bins=bins)
	@test threshold(res) ≈ 112.0 atol=2.0

	# Test 4
	res = separation(LMCLUS.Otsu, generate_sample(N, 0.5, 50, 4, 0.5, 150, 30), bins=bins)
	@test threshold(res) ≈ 121.0 atol=2.

	# Try unimodal histogram
	res = separation(LMCLUS.Otsu, rand(Normal(1, 10), N), bins=bins)
	@test threshold(res) ≈ 1.0 atol=1.

	# Extreme case
	N2 = 10
	D = [rand(N2); 10.0 .+ rand(N2)]
	res = separation(LMCLUS.Kittler, D, bins=7)
	@test threshold(res) ≈ 1.0 atol=1.0

	res = separation(LMCLUS.Otsu, D, bins=7)
	@test_broken threshold(res) ≈ 1.0 atol=1.0

end


text:
# Data structure and methods for degrees of freedom
# Alex Tam, 12/10/2020

"Mutable data structure of degrees of freedom"
mutable struct State{T}
    an::Vector{Vector{Vector{T}}} # Dimensionless positions (x, y) of actin filament nodes
    mp::Vector{Vector{T}} # Relative positions (m1, m2) of myosin motors along attached filaments
end

"Convert the state structure to a vector for use in the optimiser"
function build_dof(s::State{Float64})
    # Type specification is required here
    an = Vector{Float64}();  mp = Vector{Float64}(); # Pre-allocate
    # Actin filament node positions
    for i = 1:length(s.an) # Loop over filaments
        for j = 1:length(s.an[i]) # Loop over nodes
            push!(an, s.an[i][j][1]); # Store node x-position
            push!(an, s.an[i][j][2]); # Store node y-position
        end
    end
    # Myosin motor relative positions
    for i = 1:length(s.mp) # Loop over motors
        push!(mp, s.mp[i][1]); # Store motor position on 1st filament
        push!(mp, s.mp[i][2]); # Store motor position on 2nd filament
    end
    return vcat(an, mp)
end

"Convert the state vector to a data structure for easier access"
function build_state(x::Vector{T}, af, mm) where {T}
    an = Vector{Vector{Vector{T}}}(undef, length(af)); mp = Vector{Vector{T}}(undef, length(mm)); # Pre-allocate
    # Actin filament node positions
    dof_count = 0; # Initialise DOF counter
    for i = 1:length(af) # Loop over filaments
        nodes = Vector{Vector{T}}(undef, length(af[i].segments)+1); # Pre-allocate all nodes on an current filament
        for j = 1:(length(af[i].segments)+1) # Loop over nodes
            nodes[j] = Vector{T}(undef,2); # Pre-allocate current node position
            nodes[j][1] = x[dof_count+1]; nodes[j][2] = x[dof_count+2]; # Obtain node positions
            dof_count += 2; # Increment counter
        end
        an[i] = nodes; # Construct vector of vectors of vectors
    end
    # Myosin motor relative positions
    for i = 1:length(mm) # Loop over motors
        mp[i] = Vector{T}(undef, 2); # Pre-allocate current motor position
        mp[i][1] = x[dof_count + (i-1)*2 + 1]; mp[i][2] = x[dof_count + (i-1)*2 + 2]; # Obtain motor positions
    end
    return State{T}(an, mp)
end


text:
using Distributions, SpecialFunctions
import Distributions: mean
"""
TSOROriginal: Original algorithm of Fearnhead & Liu 2007
            for the Categorical task (section 2.3.2)
Update functions can be found in tbayesfilter.jl
The type of approximation can be found in approximations.jl
"""
struct TSOR{D, A} <: BayesFilter
    ns::Int
    N::Int
    posteriors::Array{D, 1}
    prior::D
    probabilities::Array{Float64, 1}
    stochasticity::Float64
    changeprobability::Float64
    approximation::A
    Ps1::Array{Float64, 1}
    seed::Any
    rng::MersenneTwister
end
function TSOR(;ns=5, stochasticity = .01,
                    changeprobability = 0.01, N = 20, seed = 3)
    prior = MyDirichlet(ones(ns) * stochasticity)
    approximation = SORApprox(N)
    rng = MersenneTwister(seed)
    TSOR(ns, N, [deepcopy(prior)], prior, [1.], stochasticity, changeprobability,
        approximation, Array{Float64, 1}(undef, ns), seed, rng)
end


text:
# __precompile__(true)
module Glimpse
using Reexport
@reexport using ColorTypes

@reexport using Overseer
@reexport using GeometryBasics
using GeometryBasics.StaticArrays
using Overseer: update
using Packing
using FreeTypeAbstraction
using Serialization
using SignedDistanceFields
using Parameters

using ModernGL
using Quaternions

using GLAbstraction
const GLA = GLAbstraction

using LinearAlgebra
using GLFW
using Observables
using Parameters

using CImGui
using CImGui.GLFWBackend
using CImGui.OpenGLBackend
using ThreadPools
const Gui = CImGui

using TimerOutputs
const to = TimerOutput()
include("utils.jl")
include("extensions.jl")
include("types.jl")
include("entities.jl")
include("maths/matrices.jl")
include("maths/vecmath.jl")
include("callbacks.jl")
include("shader.jl")
include("geometries.jl")
include("marching_cubes.jl")

export RGBAf0
export RGBf0

#package exports, types & enums
export Diorama
export expose, center_camera!
# For now only one context allowed, could change later I guess

const GLFW_context = Ref{GLFW.Window}()
#package exports, default geometries
end # module


text:
# infrastructure.jl

module Infrastructure

include("./db.jl")
include("./doc.jl") # database functions
include("./counter.jl") # sequence number invoice

using ..AppliAR

using CSV

import ..AppliAR: Domain, API
using .Domain
using .API

import AppliSales: Order
import AppliGeneralLedger: JournalEntry
using Dates, Sockets

export process, read_bank_statements, retrieve_unpaid_invoices, retrieve_paid_invoices, connect, disconnect

#CONST start_invoice_nbr = 1000
const START_INVOICE_NBR = 1000
#const file_invoice_nbr = "./invoicenbr.txt"
const FILE_INVOICE_NBR = "./invoicenbr.txt"
#const file_unpaid_invoices = "./test_invoicing.txt"
const FILE_UNPAID_INVOICES = "./test_invoicing.txt"
#const file_paid_invoices = "./test_invoicing_paid.txt"
const FILE_PAID_INVOICES = "./test_invoicing_paid.txt"

const SALES = 8000
const BANK = 1150
const AR = 1300
const VAT = 4000

read_bank_statements(path::String) = begin
    # read the CSV file containing bank statements
    df = CSV.read(path) # returns a DataFrame

    # return an array with BankStatement's
    # row[1] is the first value of row, row[2] the second value, etc.
    return [BankStatement(row[1], row[2], row[3], row[4]) for row in eachrow(df)]
end # read_bank_statements

process(orders::Array{Order, 1}; path=FILE_UNPAID_INVOICES) = begin
    # get last invoice number
    #try
    #    read_from_file(FILE_INVOICE_NBR)
    #catch e
    #    add_to_file(FILE_INVOICE_NBR, [START_INVOICE_NBR])
    #end

    #invnbr = last(read_from_file(FILE_INVOICE_NBR))

    # run in kubernetes
    #key = haskey(ENV, "HOSTNAME") ? split(ENV["HOSTNAME"], "-")[1] * split(ENV["HOSTNAME"], "-")[3] : "A"

    # create invoices
    #invoices = [create(order, key * "-" * string(invnbr += 1)) for order in orders]
    #invoices = [create(order, key * "-" * string(retrieve_invoice_nbr("ABC", ip"192.168.2.40", 30014))) for order in orders]

    #invoices = [create(order, key * "-" * string(retrieve_invoice_nbr(ENV["CNTNAME"], IPv4(ENV["CNTIP"]), parse(Int64, ENV["CNTPORT"])))) for order in orders]
    cntname =  readline("/etc/name/cntname")
    cntip = readline("/etc/name/cntip")
    cntport = readline("/etc/name/cntport")
    #invoices = [create(order, string(retrieve_invoice_nbr(ENV["CNTNAME"], IPv4(ENV["CNTIP"]), parse(Int64, ENV["CNTPORT"])))) for order in orders]
    invoices = [create(order, string(retrieve_invoice_nbr(cntname, IPv4(cntip), parse(Int64, cntport)))) for order in orders]

    # save invoice number
    #add_to_file(FILE_INVOICE_NBR, [invnbr])

    # archive invoices
    add_to_file(path, invoices)

    # create journal entries from invoices
    return entries = [conv2entry(inv, AR, SALES) for inv in invoices]

end # process orders

process(invoices::Array{UnpaidInvoice, 1}, stms::Array{BankStatement, 1}; path=FILE_PAID_INVOICES) = begin

    # create array with potential paid invoices based on received bank statements
    paid_invoices = PaidInvoice[]
    for unpaid_invoice in invoices
      for s in stms # get potential paid invoices
        if occursin(id(unpaid_invoice), descr(s)) # description contains invoice number
          push!(paid_invoices, create(unpaid_invoice, s))
        end
      end
    end

    # archive paid invoices
    add_to_file(path, paid_invoices)

    # return array with JournalEntry's
    return entries = [conv2entry(inv, BANK, AR) for inv in paid_invoices]
end # process invoices

retrieve_unpaid_invoices(;path=FILE_UNPAID_INVOICES)::Array{UnpaidInvoice, 1} = begin

    # retrieve unpaid invoices
    unpaid_invoices = UnpaidInvoice[invoice for invoice in read_from_file(path)]

    # return the array with UnpaidInvoice's that
    return unpaid_invoices
end # retrieve_unpaid_invoices

retrieve_paid_invoices(;path=FILE_PAID_INVOICES)::Array{PaidInvoice, 1} = begin

    # retrieve unpaid invoices as dataframe
    paid_invoices = PaidInvoice[invoice for invoice in read_from_file(path)]

    # return the array with UnpaidInvoice's that
    return paid_invoices
end # retrieve_unpaid_invoices

end # module

max_stars_repo_name: rbontekoe/AppliAR

max_issues_repo_name: rbontekoe/AppliAR

max_forks_repo_name: rbontekoe/AppliAR


text:
mutable struct SMS_EMOA <: AbstractParameters
    N::Int
    η_cr::Float64
    p_cr::Float64
    η_m::Float64
    p_m::Float64
    n_samples::Int
end

include("update-population.jl")

"""
    SMS_EMOA(;
        N = 100,
        η_cr = 20,
        p_cr = 0.9,
        η_m = 20,
        p_m = 1.0 / D,
        n_samples = 10_000,
        information = Information(),
        options = Options(),
    )

Parameters for the metaheuristic SMS-EMOA.

Parameters:

- `N` Population size.
- `η_cr`  η for the crossover.
- `p_cr` Crossover probability.
- `η_m`  η for the mutation operator.
- `p_m` Mutation probability (1/D for D-dimensional problem by default).
- `n_samples` number of samples to approximate hypervolume in many-objective (M > 2).

To use SMS_EMOA, the output from the objective function should be a 3-touple
`(f::Vector, g::Vector, h::Vector)`, where `f` contains the objective functions,
`g` and `h` are inequality, equality constraints respectively.

A feasible solution is such that `g_i(x) ≤ 0 and h_j(x) = 0`.

```julia
using Metaheuristics

# Dimension
D = 2

# Objective function
f(x) = ( x, [sum(x.^2) - 1], [0.0] )

# bounds
bounds = [-1 -1;
           1  1.0
        ]

# define the parameters (use `SMS_EMOA()` for using default parameters)
sms_emoa = SMS_EMOA(N = 100, p_cr = 0.85)

# optimize
status = optimize(f, bounds, sms_emoa)

# show results
display(status)
```

"""
function SMS_EMOA(;
    N = 100,
    η_cr = 20,
    p_cr = 0.9,
    η_m = 20,
    p_m = -1,
    n_samples = 10_000,
    information = Information(),
    options = Options(),
)

    parameters = SMS_EMOA(N, promote( Float64(η_cr), p_cr, η_m, p_m)..., n_samples)
    Algorithm(
        parameters,
        information = information,
        options = options,
    )

end

function update_state!(
    status::State,
    parameters::SMS_EMOA,
    problem::AbstractProblem,
    information::Information,
    options::Options,
    args...;
    kargs...
    )

    I = randperm(parameters.N)
    J = randperm(parameters.N)

    for i = 1:parameters.N

        pa = tournament_selection(status.population, I[i])
        pb = tournament_selection(status.population, J[i])

        # crossover
        _, c = SBX_crossover(get_position(pa),
                             get_position(pb),
                             problem.bounds,
                             parameters.η_cr,
                             parameters.p_cr)

        # mutation
        polynomial_mutation!(c, problem.bounds, parameters.η_m, parameters.p_m)

        # repair solutions if necesary
        reset_to_violated_bounds!(c, problem.bounds)

        # evaluate offspring
        offspring = create_solution(c, problem)

        # Reduce: save offspring in population and reduce population according to
        # front contribution
        update_population!(status.population, offspring, parameters.n_samples)
    end

end

function initialize!(
    status,
    parameters::SMS_EMOA,
    problem::AbstractProblem,
    information::Information,
    options::Options,
    args...;
    kargs...
)
    D = size(problem.bounds, 2)

    if parameters.p_m < 0.0
        parameters.p_m = 1.0 / D
    end

    if options.iterations == 0
        options.iterations = 500
    end

    if options.f_calls_limit == 0
        options.f_calls_limit = options.iterations * parameters.N + 1
    end

    status = gen_initial_state(problem,parameters,information,options,status)

    status

end

function final_stage!(
    status::State,
    parameters::SMS_EMOA,
    problem::AbstractProblem,
    information::Information,
    options::Options,
    args...;
    kargs...
    )
    status.final_time = time()

    #status.best_sol = get_pareto_front(status.population, is_better)

end


text:
__precompile__(true)

module DynClust
using StatsBase
using Distributions
using ComputationalResources
using Distributed
using ProgressMeter
using Images
using OffsetArrays

function __init__()
    # Enable `using` to load additional modules in this folder
    push!(LOAD_PATH, dirname(@__FILE__))
    # Now check for any resources that your package supports
    if haveresource(ArrayFireLibs)
        # User has indicated support for the ArrayFire libraries, so load your relevant code
        @eval using DynClustAF
    end
    # Put additional resource checks here
    # Don't forget to clean up!
    pop!(LOAD_PATH)
end

include("utilityFunctions.jl")
include("multiTest.jl")
include("denoising.jl")
include("clustering.jl")

export  runDenoising,
        getDenoisingResults,
        runClustering,
        getClusteringResults

end # module


text:
using Main.CaesarM: rot13

while true

    println("\nText to shift [leave blank to break]: ")
    human_input = readline()
    human_input == "" ? break :
    println("\nAlphabet [leave blank for default]: ")
    human_alphabet = readline()
    human_alphabet == "" ? human_alphabet = "abcdefghijklmnopqrstuvwxyz" :
    println()

    println("Here's our decryption:")
    println(rot13(human_input, human_alphabet))

end # while

max_stars_repo_name: Squalm/Cipher-Tools

max_issues_repo_name: Squalm/Cipher-Tools

max_forks_repo_name: Squalm/Cipher-Tools


text:
using InvertibleNetworks, Test, LinearAlgebra, Random

Random.seed!(11)

##############################################################################################################################
# 4D Tensor
X = randn(Float32, 28, 28, 2, 4)

# Squeeze and unsqueeze
@test isapprox(norm(X - unsqueeze(squeeze(X; pattern="column"); pattern="column")), 0f0; atol=1f-6)
@test isapprox(norm(X - unsqueeze(squeeze(X; pattern="patch"); pattern="patch")), 0f0; atol=1f-6)
@test isapprox(norm(X - unsqueeze(squeeze(X; pattern="checkerboard"); pattern="checkerboard")), 0f0; atol=1f-6)

#test norm preservation of orthogonal transforms
@test norm(X) ≈ norm(Haar_squeeze(X))

# Wavelet transform invertibility
Y = wavelet_squeeze(X)
X_ = wavelet_unsqueeze(Y)
@test isapprox(norm(X - X_)/norm(X), 0f0; atol=1f-5)
@test X ≈ invHaar_unsqueeze(Haar_squeeze(X))

# Wavelet transform adjoint test
X = randn(Float32, 28, 28, 2, 4)
Y = randn(Float32, 14, 14, 8, 4)
a = dot(Y, wavelet_squeeze(X))
b = dot(X, wavelet_unsqueeze(Y))
@test isapprox(a/b - 1f0, 0f0; atol=1f-5)
a = dot(Y, Haar_squeeze(X))
b = dot(X, invHaar_unsqueeze(Y))
@test isapprox(a/b - 1f0, 0f0; atol=1f-5)

# Split and cat
@test isapprox(norm(X - tensor_cat(tensor_split(X))), 0f0; atol=1f-5)

# Unsqueeze works for any size as long as channel size divisible by 4 (for 4D Tensor)
x_size  = 27
y_size  = 27
ch_size = 8
batch_size = 4
X = randn(Float32, x_size, y_size, ch_size, batch_size)
@test isequal(size(unsqueeze(X)), (x_size*2, y_size*2, div(ch_size,4), batch_size))

##############################################################################################################################
# 5D Tensor
X = randn(Float32, 16, 16, 16, 2, 4)

# Squeeze and unsqueeze
@test isapprox(norm(X - unsqueeze(squeeze(X; pattern="column"); pattern="column")), 0f0; atol=1f-6)
@test isapprox(norm(X - unsqueeze(squeeze(X; pattern="patch"); pattern="patch")), 0f0; atol=1f-6)

#test norm preservation of orthogonal transforms
@test norm(X) ≈ norm(Haar_squeeze(X))

# Wavelet transform invertibility
Y = wavelet_squeeze(X)
X_ = wavelet_unsqueeze(Y)
@test isapprox(norm(X - X_)/norm(X), 0f0; atol=1f-5)
@test X ≈ invHaar_unsqueeze(Haar_squeeze(X))

# Wavelet transform adjoint test
X = randn(Float32, 16, 16, 16, 2, 4)
Y = randn(Float32, 8, 8, 8, 16, 4)
a = dot(Y, wavelet_squeeze(X))
b = dot(X, wavelet_unsqueeze(Y))
@test isapprox(a/b - 1f0, 0f0; atol=1f-5)
a = dot(Y, Haar_squeeze(X))
b = dot(X, invHaar_unsqueeze(Y))
@test isapprox(a/b - 1f0, 0f0; atol=1f-5)

# Split and cat
@test isapprox(norm(X - tensor_cat(tensor_split(X))), 0f0; atol=1f-5)

# Unsqueeze works for any size as long as channel size divisible by 8 (for 5D Tensor)
x_size  = 27
y_size  = 27
z_size  = 27
ch_size = 8
batch_size = 4
X = randn(Float32, x_size, y_size, z_size, ch_size, batch_size)
@test isequal(size(unsqueeze(X)), (x_size*2, y_size*2, z_size*2, div(ch_size,8), batch_size))


text:
#Helper function for a regularizer that's actually a constraint
function inf_or_zero(r::Regularizer, a::AbstractMatrix)
  for i in 1:size(a,2)
    if evaluate(r, a[:,i]) === Inf
      return Inf
    else
      continue
    end
  end
  return 0
end

## Takes a prox step for each column, mutable method
function colwise_prox!(r::Regularizer, u::AbstractMatrix, alpha::Number)
  for i in 1:size(u,2)
    prox!(r, view(u, :, i), alpha)
  end
  u
end

## Immutable wrapper around colwise_prox!
function colwise_prox(r::Regularizer, u::AbstractMatrix, alpha::Number)
  v = copy(u)
  colwise_prox!(r, v, alpha)
  v
end

## Unit one-sparse constraint
##
function evaluate(r::UnitOneSparseConstraint, a::AbstractMatrix)
  inf_or_zero(r, a)
end

function prox!(r::UnitOneSparseConstraint, u::Matrix, alpha::Number)
  colwise_prox!(r, u, alpha)
end

function prox(r::UnitOneSparseConstraint, u::AbstractMatrix, alpha::Number)
  colwise_prox(r, u, alpha)
end

## K-sparse constraint
##
function evaluate(r::KSparseConstraint, a::AbstractMatrix)
  inf_or_zero(r, a)
end

function prox!(r::KSparseConstraint, u::Matrix, alpha::Number)
  colwise_prox!(r, u, alpha)
end

function prox(r::KSparseConstraint, u::AbstractMatrix, alpha::Number)
  colwise_prox(r, u, alpha)
end

## One-sparse constraint
##
function evaluate(r::OneSparseConstraint, a::AbstractMatrix)
  inf_or_zero(r, a)
end

function prox!(r::OneSparseConstraint, u::Matrix, alpha::Number)
  colwise_prox!(r, u, alpha)
end

function prox(r::OneSparseConstraint, u::AbstractMatrix, alpha::Number)
  colwise_prox(r, u, alpha)
end


text:
module InternalBallistics
end


text:
# These tests are designed to ensure that Immerse.jl works. If you
# need to edit these tests to make them pass, that's fine, but please
# submit the corresponding fix to Immerse.

using Compose, Base.Test
import Cairo

### The Immerse backend
srf = Cairo.CairoImageSurface(10, 10, Cairo.FORMAT_RGB24)
ctx = compose(compose(context(), rectangle(0.0w, 0.0h, 0.8w, 0.7h, :rect)), fill("tomato"))
be = Compose.ImmerseBackend(srf)
draw(be, ctx)
@test be.coords[:rect][2] == Compose.UnitBox()

### Finding tagged objects
const ContainersWithChildren = Union{Context,Compose.Table}
const Iterables = Union{ContainersWithChildren, AbstractArray}

iterable(ctx::ContainersWithChildren) = ctx.children
iterable(a::AbstractArray) = a

function find_tagged(root)
    handles = Dict{Symbol,Context}()
    find_tagged!(handles, root)
end

function find_tagged!(handles, obj::Iterables)
    for item in iterable(obj)
        if has_tag(item)
            handles[item.tag] = obj
        else
            find_tagged!(handles, item)
        end
    end
    handles
end

function find_tagged!(handles, obj::Context)
    for item in obj.form_children
        if has_tag(item)
            handles[item.tag] = obj
        end
    end
    for item in obj.container_children
        find_tagged!(handles, item)
    end
    handles
end

find_tagged!(handles, obj) = obj

has_tag(form::Compose.Form, tag) = form.tag == tag
has_tag(form::Compose.Form)      = form.tag != Compose.empty_tag

has_tag(obj, tag) = false
has_tag(obj)      = false

@test find_tagged(ctx)[:rect] == ctx

### Coordinate computations
function absolute_to_data(x, y, transform, unit_box, parent_box)
    xt, yt = invert_transform(transform, x, y)
    (unit_box.x0 + unit_box.width *(xt-parent_box.x0[1])/Measures.width(parent_box),
     unit_box.y0 + unit_box.height*(yt-parent_box.x0[2])/Measures.height(parent_box))
end

invert_transform(::Compose.IdentityTransform, x, y) = x, y

function invert_transform(t::Compose.MatrixTransform, x, y)
    @assert t.M[3,1] == t.M[3,2] == 0
    xyt = t.M\[x, y, 1.0]
    xyt[1], xyt[2]
end

box, units, transform = be.coords[:rect]
xd, yd = absolute_to_data(0.5mm, 0.5mm, transform, units, box)
@test isapprox(xd,0.1417; atol=0.0001)


text:
module LightLattices

using StaticArrays, RecursiveArrayTools, DocStringExtensions

export AbstractNodeCollection

"""
$(TYPEDEF)

Abstract collection of nodes in `D`-dimensional space. Type `T` is used to represent coordinates.
"""
abstract type AbstractNodeCollection{D, T} end

"""
$(TYPEDSIGNATURES)

Returns the coordinate of node with index `i1` relative to coordinate of the node with index `i2`.
"""
Base.@propagate_inbounds function relative_coordinate(collection::AbstractNodeCollection, i1, i2) end

include("cells.jl")
include("lattices.jl")
include("utils.jl")

export AbstractNodeCollection
export AbstractCell, TrivialCell, HomogeneousCell, InhomogeneousCell
export AbstractLattice, RegularLattice
export switch_coord_type
export relative_coordinate
export num_of_groups, group_size

end


text:
# strategy:
# * @model macro converts formula exprs in original callinto Formula calls and
#   plugs into ModelBuilder{M} where M is the original call
# * If provided data, calls ModelFrame{M} constructor with builder and data.
#   This summarizes and then sets schema on all the formulas in args/kw,

mutable struct ModelBuilder{M}
    args::Tuple
    kw::Dict{Symbol,Any}
    formula_idx::Vector{Int}
    formula_kws::Vector{Symbol}
end

ModelBuilder(head, args...; kw...) = ModelBuilder{head}(args, Dict(kw))
function ModelBuilder{M}(args::Tuple, kw::Dict) where M
    formula_idx = find(x -> x isa Formula, args)
    formula_kws = Symbol[k for (k,v) in kw if isa(v, Formula)]
    ModelBuilder{M}(args, kw, formula_idx, formula_kws)
end

_replace_formula!(x) = x
function _replace_formula!(ex::Expr)
    if Meta.isexpr(ex, :kw)
        ex.args[2] = _replace_formula!(ex.args[2])
        return ex
    end

    if is_formula(ex)
        return _formula!(ex)
    else
        return ex
    end
end

"""
    @model f(args...; kw...)

Capture formula arguments to a model constructor (or any other) function, and
transform to a `ModelBuilder` constructor with `Formula`s in their place.

"""
macro model(ex)
    @argcheck check_call(ex)
    # strategy is to replace each argument that's a formula with an expr for a
    # formula, then plug into a call to ModelBuilder
    builder = Expr(:call, :(StreamModels.ModelBuilder), ex.args[1],
                   _replace_formula!.(ex.args[2:end])...)
    # prevent macro hygine from transforming captured variables to refer to
    # `StatsModels`
    return esc(builder)
end

"""
    build(builder::ModelBuilder{M}, source) where M
    build(builder::ModelBuilder{M}, source, schema::Data.Schema) where M

Construct a model of type M from a model builder and a data source.  This
extracts the Formula [keyword] arguments, sets the schema on them, performs any
necessary summarization, and instantiates them as numeric arrays, before calling
the constructor.

"""
function build(builder::ModelBuilder{M}, source) where M
    sch = summarize(source,
                    builder.args[builder.formula_idx]...,
                    builder.kw[builder.formula_kw]...)
    for idx in builder.formula_idx
        set_schema!(builder.args[idx], sch)
    end
    for kw in builder.formula_kws
        set_schema!(builder.kw[kw], sch)
    end

    ## call model matrix and plug into constructor
end

# """
#     predict(mf::ModelFrame)

# Predicted values for the model (needs to be fitted first; just delegates to the
# stored model

# """
# StatsBase.predict(mf::ModelFrame) = predict(mf.model)

# """
#     predict(mf::ModelFrame{M,S}, source::S) where {M,S}

# Generate predictions for new data.  The source is converted into numeric
# arguments and plugged in in the same order they occur in the builder, without
# the non-formula arguments.

# Packages can customize this behavior by specializing on the model type.

# """
# function StatsBase.predict(mf::ModelFrame{M,S}, source::S) where {M,S} end

# mutable struct ModelFrame{M,S}
#     model::M
#     source::S
#     schema::Data.Schema
# end


text:
module VisualSearchACTR
    using Reexport, Distributions, StatsBase, ArgCheck
    using Crayons, ACTRModels, ConcreteStructs
    @reexport using Random, ACTRModels
    import ACTRModels: AbstractParms, AbstractACTR
    export ACTRV, Experiment, VisualObject, Feature, Data, Fixation, Parm
    export populate_features, populate_visicon, initialize_trial!
    export feature_visibility!, compute_angular_distance, compute_acuity_threshold
    export compute_angular_size, update_visibility!, run_trial!, run_condition!
    export orient!, compute_activations!, update_decay!, draw_object!
    export conjunctive_ratio, conjunctive_set, feature_set, generate_stimuli
    export fixation_probs, fixation_prob, import_gui
    include("structs.jl")
    include("model.jl")
    include("task.jl")
    include("trace.jl")
end