text:
#-------------------------------------------------------------------------------
# Test suite for 'lowder' function
#-------------------------------------------------------------------------------

@testset "lowder" begin

    @testset "linear models" begin

        # --------------------------------------------------------------------------------
        # Test 01 - Two linear functions (sol.index == 1)
        # --------------------------------------------------------------------------------

        function f1(x)
            return 2.0 * x[1] - 5.0 * x[2]
        end

        function f2(x)
            return - 1.5 * x[1] + x[2]
        end

        fmin_list = [f1, f2]
        a = [0.0, 0.0]
        b = [5.0, 5.0]

        x = [1.0, 2.0]

        sol = LOWDER.lowder(fmin_list, x, a, b; δ = 1.0, Δ = 1.5, m = 3)

        @test( sol.status == :success )
        @test( sol.index == 1 )
        @test( sol.solution == [0.0, 5.0] )

        # --------------------------------------------------------------------------------
        # Test 02 - Two linear functions (sol.index == 2)
        # --------------------------------------------------------------------------------

        x = [4.0, 1.0]

        sol = LOWDER.lowder(fmin_list, x, a, b; δ = 1.0, Δ = 1.5, m = 3)

        @test( sol.status == :success )
        @test( sol.index == 2 )
        @test( sol.solution == [5.0, 0.0] )

        # --------------------------------------------------------------------------------
        # Test 03 - A linear function and a quadratic function (sol.index == 1)
        # --------------------------------------------------------------------------------

        function f3(x)
            return - 0.5 * x[1] + 10 * x[2]
        end

        function f4(x)
            return x[1] ^ 2.0 + x[2] ^ 2.0
        end

        fmin_list = [f3, f4]
        a = [0.0, 0.0]
        b = [5.0, 5.0]

        x = [4.0, 1.0]

        sol = LOWDER.lowder(fmin_list, x, a, b; δ = 1.0, Δ = 1.2, m = 3)

        @test( sol.status == :success )
        @test( sol.index == 1 )
        @test( sol.solution == [5.0, 0.0] )

        # --------------------------------------------------------------------------------
        # Test 04 - A linear function and a quadratic function (sol.index == 2)
        # --------------------------------------------------------------------------------

        x = [1.0, 1.0]

        sol = LOWDER.lowder(fmin_list, x, a, b; δ = 1.0, Δ = 1.2, m = 3)

        @test( sol.status == :success )
        @test( sol.index == 2 )
        @test( sol.solution == [0.0, 0.0] )

        # --------------------------------------------------------------------------------
        # Test 05 - Two quadratic functions (sol.index == 1)
        # --------------------------------------------------------------------------------

        function f5(x)
            return - x[2] ^ 2.0 + x[1]
        end

        function f6(x)
            return 0.1 * x[1] ^ 2.0 - x[1] * x[2]
        end

        fmin_list = [f5, f6]
        a = [0.0, 0.0]
        b = [5.0, 5.0]

        x = [1.0, 3.0]

        sol = LOWDER.lowder(fmin_list, x, a, b; δ = 0.5, Δ = 1.0, m = 3)

        @test( sol.status == :success )
        @test( sol.index == 1 )
        @test( sol.solution == [0.0, 5.0] )

        # --------------------------------------------------------------------------------
        # Test 06 - Two quadratic functions (sol.index == 2)
        # --------------------------------------------------------------------------------

        x = [5.0, 1.0]

        sol = LOWDER.lowder(fmin_list, x, a, b; δ = 0.5, Δ = 1.0, m = 3)

        @test( sol.status == :success )
        @test( sol.index == 2 )
        @test( sol.solution == [5.0, 5.0] )

        # --------------------------------------------------------------------------------
        # Test 07 - Rosenbrock & Freudenstein and Roth functions
        #           (full_active_set and η = 0.0)
        # --------------------------------------------------------------------------------

        function f7(x)
            return ( 10.0 * ( x[2] - x[1] ^ 2.0 ) ) ^ 2.0 + ( 1.0 - x[1] ) ^ 2.0
        end

        function f8(x)
            return ( -13.0 + x[1] + ( ( 5.0 - x[2] ) * x[2] - 2.0 ) * x[2] ) ^ 2.0 + ( -29.0 + x[1] + ( ( x[2] + 1.0 ) * x[2] - 14.0 ) * x[2] ) ^ 2.0
        end

        fmin_list = [f7, f8]
        a = [0.0, 0.0]
        b = [5.0, 5.0]

        x = [5.0, 0.0]

        sol = LOWDER.lowder(fmin_list, x, a, b; δ = 1.0, Δ = 1.2, m = 3, η = 0.0)

        @test( sol.status == :success )
        @test( sol.index == 2 )
        @test( sol.solution == [5.0, 0.0] )

        # --------------------------------------------------------------------------------
        # Test 08 - Rosenbrock & Freudenstein and Roth functions (full_active_set)
        # --------------------------------------------------------------------------------

        x = [5.0, 0.0]

        sol = LOWDER.lowder(fmin_list, x, a, b; δ = 1.0, Δ = 1.2, m = 3)

        @test( sol.status == :success )
        @test( sol.index == 2 )
        @test( sol.solution == [5.0, 0.0] )

        # --------------------------------------------------------------------------------
        # Test 09 - Rosenbrock & Freudenstein and Roth functions
        # --------------------------------------------------------------------------------

        x = [5.0, 3.0]

        sol = LOWDER.lowder(fmin_list, x, a, b; δ = 1.0, Δ = 2.0, m = 3)

        @test( sol.status == :stalled )
        @test( sol.index == 2 )
        @test( sol.solution == [5.0, 4.0] )

        # --------------------------------------------------------------------------------
        # Test 10 - Rosenbrock & Freudenstein and Roth functions
        # --------------------------------------------------------------------------------

        x = [5.0, 5.0]

        sol = LOWDER.lowder(fmin_list, x, a, b; δ = 1.0, Δ = 2.0, m = 3)

        @test( sol.status == :stalled )
        @test( sol.index == 2 )
        @test( sol.solution == [5.0, 4.0] )

        # --------------------------------------------------------------------------------
        # Test 11 - Linear and Quadratic functions
        # --------------------------------------------------------------------------------

        function f9(x)
            return x[1] ^ 2.0 + x[2] ^ 2.0
        end

        function f10(x)
            return 0.5 * x[1] + 0.1 * x[2] + 1.0
        end

        fmin_list = [f9, f10]
        a = [- 1.0, - 2.0]
        b = [1.5, 2.0]

        x = [1.5, 1.0]

        sol = LOWDER.lowder(fmin_list, x, a, b; δ = 0.5, Δ = 0.5, m = 3)

        @test( sol.status == :stalled )
        @test( sol.index == 1 )
        @test( isapprox( sol.solution, [0.0, 0.0], atol=1e-4 ) )

    end

end

max_stars_repo_name: aschwertner/LOWDER

max_issues_repo_name: aschwertner/LOWDER

max_forks_repo_name: aschwertner/LOWDER


text:
using Winston

p = FramedPlot(
        aspect_ratio=1,
        xrange=(0,100),
        yrange=(0,100))

n = 21
x = linspace(0, 100, n)
yA = 40 + 10randn(n)
yB = x + 5randn(n)

a = Points(x, yA, kind="circle")
setattr(a, label="a points")

b = Points(x, yB)
setattr(b, label="b points")
style(b, kind="filled circle")

s = Slope(1, (0,0), kind="dotted")
setattr(s, label="slope")

l = Legend(.1, .9, {a,b,s})

add(p, s, a, b, l)

file(p, "example2.png")

ext: jl


text:
# Julia 0.6 macros to aid in vectorization: @view, @views, @__dot__ (@.),
# backported from Julia 0.6.

# prior to julia#20247, the replace_ref_end! macro had hygiene bugs
if VERSION < v"0.6.0-dev.2406"
    function trailingsize(A, n)
        s = 1
        for i=n:ndims(A)
            s *= size(A,i)
        end
        return s
    end
    replace_ref_end!(ex) = replace_ref_end_!(ex, nothing)[1]
    # replace_ref_end_!(ex,withex) returns (new ex, whether withex was used)
    function replace_ref_end_!(ex, withex)
        used_withex = false
        if isa(ex,Symbol) && ex == :end
            withex === nothing && error("Invalid use of end")
            return withex, true
        elseif isa(ex,Expr)
            if ex.head == :ref
                ex.args[1], used_withex = replace_ref_end_!(ex.args[1],withex)
                S = isa(ex.args[1],Symbol) ? ex.args[1]::Symbol : gensym(:S) # temp var to cache ex.args[1] if needed
                used_S = false # whether we actually need S
                # new :ref, so redefine withex
                nargs = length(ex.args)-1
                if nargs == 0
                    return ex, used_withex
                elseif nargs == 1
                    # replace with endof(S)
                    ex.args[2], used_S = replace_ref_end_!(ex.args[2],:($endof($S)))
                else
                    n = 1
                    J = endof(ex.args)
                    for j = 2:J-1
                        exj, used = replace_ref_end_!(ex.args[j],:($size($S,$n)))
                        used_S |= used
                        ex.args[j] = exj
                        if isa(exj,Expr) && exj.head == :...
                            # splatted object
                            exjs = exj.args[1]
                            n = :($n + length($exjs))
                        elseif isa(n, Expr)
                            # previous expression splatted
                            n = :($n + 1)
                        else
                            # an integer
                            n += 1
                        end
                    end
                    ex.args[J], used = replace_ref_end_!(ex.args[J],:($trailingsize($S,$n)))
                    used_S |= used
                end
                if used_S && S !== ex.args[1]
                    S0 = ex.args[1]
                    ex.args[1] = S
                    ex = Expr(:let, ex, :($S = $S0))
                end
            else
                # recursive search
                for i = eachindex(ex.args)
                    ex.args[i], used = replace_ref_end_!(ex.args[i],withex)
                    used_withex |= used
                end
            end
        end
        ex, used_withex
    end
end

# convert x[...] on lhs of .= to a view in broadcast! call
if VERSION < v"0.5.0-dev+5575" #17510
    dotview(args...) = getindex(args...)
    dotview(A::AbstractArray, args...) = view(A, args...)
    dotview{T<:AbstractArray}(A::AbstractArray{T}, args::Integer...) = getindex(A, args...)
    todotview(x) = x
    function todotview(ex::Expr)
        if ex.head == :ref
            ex = replace_ref_end!(ex)
            if Meta.isexpr(ex, :ref)
                ex = Expr(:call, dotview, ex.args...)
            else # ex replaced by let ...; foo[...]; end
                assert(Meta.isexpr(ex, :let) && Meta.isexpr(ex.args[1], :ref))
                ex.args[1] = Expr(:call, dotview, ex.args[1].args...)
            end
        end
    end
end

if !isdefined(Base, Symbol("@view"))
    macro view(ex)
        if Meta.isexpr(ex, :ref)
            ex = replace_ref_end!(ex)
            if Meta.isexpr(ex, :ref)
                ex = Expr(:call, view, ex.args...)
            else # ex replaced by let ...; foo[...]; end
                assert(Meta.isexpr(ex, :let) && Meta.isexpr(ex.args[1], :ref))
                ex.args[1] = Expr(:call, view, ex.args[1].args...)
            end
            Expr(:&&, true, esc(ex))
        else
            throw(ArgumentError("Invalid use of @view macro: argument must be a reference expression A[...]."))
        end
    end
    export @view
end

if !isdefined(Base, Symbol("@views"))
    maybeview(A, args...) = getindex(A, args...)
    maybeview(A::AbstractArray, args...) = view(A, args...)
    maybeview(A::AbstractArray, args::Number...) = getindex(A, args...)
    maybeview(A) = getindex(A)
    maybeview(A::AbstractArray) = getindex(A)

    _views(x) = x
    function _views(ex::Expr)
        if ex.head in (:(=), :(.=))
            # don't use view for ref on the lhs of an assignment,
            # but still use views for the args of the ref:
            lhs = ex.args[1]
            Expr(ex.head, Meta.isexpr(lhs, :ref) ?
                          Expr(:ref, map(_views, lhs.args)...) : _views(lhs),
                 _views(ex.args[2]))
        elseif VERSION < v"0.5.0-dev+5575" && isexpr(ex, :comparison) && ex.args[2] == :.= #17510
            # as above, but in Julia 0.4 a .= produces a comparison expression
            lhs_ = ex.args[1]
            lhs = Meta.isexpr(lhs_, :ref) ? Expr(:ref, map(_views, lhs_.args)...) : _views(lhs_)
            if length(ex.args) == 3
                Expr(:.=, lhs, _views(ex.args[3]))
            else
                Expr(:.=, lhs, _views(Expr(:comparison, ex.args[3:end]...)))
            end
        elseif ex.head == :ref
            Expr(:call, maybeview, map(_views, ex.args)...)
        else
            h = string(ex.head)
            # don't use view on the lhs of an op-assignment a[i...] += ...
            if last(h) == '=' && Meta.isexpr(ex.args[1], :ref)
                lhs = ex.args[1]

                # temp vars to avoid recomputing a and i,
                # which will be assigned in a let block:
                a = gensym(:a)
                i = [gensym(:i) for k = 1:length(lhs.args)-1]

                # for splatted indices like a[i, j...], we need to
                # splat the corresponding temp var.
                I = similar(i, Any)
                for k = 1:length(i)
                    if Meta.isexpr(lhs.args[k+1], :...)
                        I[k] = Expr(:..., i[k])
                        lhs.args[k+1] = lhs.args[k+1].args[1] # unsplat
                    else
                        I[k] = i[k]
                    end
                end

                Expr(:let,
                     Expr(first(h) == '.' ? :(.=) : :(=), :($a[$(I...)]),
                          Expr(:call, Symbol(h[1:end-1]),
                               :($maybeview($a, $(I...))),
                               map(_views, ex.args[2:end])...)),
                     :($a = $(_views(lhs.args[1]))),
                     [:($(i[k]) = $(_views(lhs.args[k+1]))) for k=1:length(i)]...)
            else
                Expr(ex.head, map(_views, ex.args)...)
            end
        end
    end

    macro views(x)
        esc(_views(replace_ref_end!(x)))
    end
    export @views
end

# we can't define @. because that doesn't parse in Julia < 0.6, but
# we can define @__dot__, which is what @. is sugar for:
if !isdefined(Base, Symbol("@__dot__"))
    dottable(x) = false # avoid dotting spliced objects (e.g. view calls inserted by @view)
    dottable(x::Symbol) = !Base.isoperator(x) || first(string(x)) != '.' || x == :.. # don't add dots to dot operators
    dottable(x::Expr) = x.head != :$
    undot(x) = x
    function undot(x::Expr)
        if x.head == :.=
            Expr(:(=), x.args...)
        elseif x.head == :block # occurs in for x=..., y=...
            Expr(:block, map(undot, x.args)...)
        else
            x
        end
    end
    __dot__(x) = x
    function __dot__(x::Expr)
        dotargs = map(__dot__, x.args)
        if x.head == :call && dottable(x.args[1])
            Expr(:., dotargs[1], Expr(:tuple, dotargs[2:end]...))
        elseif x.head == :$
            x.args[1]
        elseif x.head == :let # don't add dots to "let x=... assignments
            Expr(:let, dotargs[1], map(undot, dotargs[2:end])...)
        elseif x.head == :for # don't add dots to for x=... assignments
            Expr(:for, undot(dotargs[1]), dotargs[2])
        elseif (x.head == :(=) || x.head == :function || x.head == :macro) &&
               Meta.isexpr(x.args[1], :call) # function or macro definition
            Expr(x.head, x.args[1], dotargs[2])
        else
            head = string(x.head)
            if last(head) == '=' && first(head) != '.'
                Expr(Symbol('.',head), dotargs...)
            else
                Expr(x.head, dotargs...)
            end
        end
    end
    macro __dot__(x)
        esc(__dot__(x))
    end
    macro dotcompat(x)
        esc(_compat(__dot__(x)))
    end
    export @__dot__, @dotcompat
else
    # in 0.6, use the __dot__ function from Base.Broadcast
    macro dotcompat(x)
        esc(_compat(Base.Broadcast.__dot__(x)))
    end
    export @dotcompat
end


text:
using PolarimetryCalibration, Plots, LaTeXStrings

σ = 0.01
PAsrc = 36.0
PArange = 0:0.2:180
const (I,Q,U,V)=(1,2,3,4)

stokes = [(müller_rec()*
           müller_sky(PAaz*π/180)*
           lin_pol_source(pa = PAsrc*π/180))[i] + σ*randn()
          for PAaz = PArange, i in 1:4]
plot(PArange, stokes[:,V], label=L"ideal", xlabel = "PA", ylabel = "Stokes Q")

stokes = [(müller_rec(ΔG = 0.1)*
           müller_sky(PAaz*π/180)*
           lin_pol_source(pa = PAsrc*π/180))[i] + σ*randn()
          for PAaz = PArange, i in 1:4]
plot!(PArange, stokes[:,V], label=L"\Delta G = 0.1")

stokes = [(müller_rec(ψ = 52*π/180)*
           müller_sky(PAaz*π/180)*
           lin_pol_source(pa = PAsrc*π/180))[i] + σ*randn()
          for PAaz = PArange, i in 1:4]
plot!(PArange, stokes[:,V], label=L"$\psi = 52$ $^\ocirc$")

max_stars_repo_name: KajWiik/PolarimetryCalibration

max_issues_repo_name: KajWiik/PolarimetryCalibration

max_forks_repo_name: KajWiik/PolarimetryCalibration


text:
# test components

@testset "--- Checking Generic Mapping object ---" begin

  # Create Test mesh for tests
  nnodes = [3,3,3]  # Number of nodes of the FE grid that need to be mapped
  nodes_xyz = zeros(nnodes[1], nnodes[2], nnodes[3], 3)
  origin = [1,1,1]
  incz = 0.0
  for k = 1:nnodes[3]
    incy = 0.0
    for j = 1:nnodes[2]
      incx = 0.0
      for i = 1:nnodes[1]
        nodes_xyz[i,j,k,1] = origin[1] + incx
        nodes_xyz[i,j,k,2] = origin[2] + incy
        nodes_xyz[i,j,k,3] = origin[3] + incz
        incx += 1
      end
      incy += 1
    end
    incz += 1
  end

  # Create Mapping Object
  ndim = 3
  order = [2,2,2]  # Order of B-splines in the 3 directions
  nControlPts = [3,3,3]
  map = Mapping(ndim, order, nControlPts, nnodes)

  @test ( map.ndim )== 3
  @test ( map.nctl )== [3,3,3]
  @test ( map.order )== [2,2,2]
  @test ( map.numnodes )== [3,3,3]

  @testset "--- Checking Knot calculations ---" begin

    calcKnot(map)
    for i = 1:map.ndim
      @test ( map.edge_knot[i][1] )== 0.0
      @test ( map.edge_knot[i][2] )== 0.0
      @test isapprox( map.edge_knot[i][3], 0.5) atol=1e-15
      @test ( map.edge_knot[i][4] )== 1.0
      @test ( map.edge_knot[i][5] )== 1.0
    end
  end  # End context("--- Checking Knot calculations ---")

end  # End facts("--- Checking Mapping object ---")

# MPI Declarations
comm = MPI.COMM_WORLD
comm_world = MPI.MPI_COMM_WORLD
comm_self = MPI.COMM_SELF
my_rank = MPI.Comm_rank(comm)
comm_size = MPI.Comm_size(comm)

# Pumi Specific Tests
@testset "--- Checking FFD Types and Functions For Full Serial DG Pumi Meshes ---" begin

  # 2D mesh
  opts = Dict{ASCIIString, Any}()
  opts["order"] = 1
  opts["dimensions"] = 2
  opts["use_DG"] = true
  opts["Tsbp"] = Float64
  opts["Tmsh"] = Complex128
  opts["operator_type"] = "SBPOmega"
  opts["dmg_name"] = "../src/mesh_files/2D_Airfoil.dmg"
  opts["smb_name"] = "../src/mesh_files/2D_Airfoil.smb"
  opts["numBC"] = 2

  # For 2DAirfoil
  opts["BC1"] = [8,11,14,17]
  opts["BC1_name"] = "FarField"
  opts["BC2"] = [5]
  opts["BC2_name"] = "Airfoil"

  opts["coloring_distance"] = 2 # 0 For CG Mesh 2 for DG Mesh
  opts["jac_type"] = 2
  opts["jac_method"] = 2
  opts["run_type"] = 5

  # Create PumiMesh and SBP objects
  dofpernode = 1
  ref_verts = [-1. 1 -1; -1 -1 1]
  Tsbp = opts["Tsbp"]
  Tmsh = opts["Tmsh"]
  sbp = getTriSBPOmega(degree=opts["order"], Tsbp=opts["Tsbp"])
  sbpface = TriFace{Tsbp}(opts["order"], sbp.cub, ref_verts.')
  topo = 0
  shape_type = 2
  mesh = PumiMeshDG2(Tmsh, sbp, opts, sbpface, dofpernode=dofpernode,
                     shape_type=shape_type)
#  mesh = PumiMeshDG2{Tmsh}(opts["dmg_name"], opts["smb_name"], opts["order"],
#                           sbp, opts, sbpface; dofpernode=dofpernode,
#                           coloring_distance=opts["coloring_distance"],
#                           shape_type=shape_type)

  @testset "--- Checking Linear Mapping For DG Mesh ---" begin

    # Free Form deformation parameters
    ndim = 2
    order = [4,4,2]  # Order of B-splines in the 3 directions
    nControlPts = [4,4,2]

    # Create Mapping object
    map = PumiMapping{Tmsh}(ndim, order, nControlPts, mesh, bc_nums=[1, 2])

    # Create knot vector
    calcKnot(map)

    # Create Bounding box
    offset = [0., 0., 0.5] # No offset in the X & Y direction
    box = PumiBoundingBox{Tmsh}(map, mesh, sbp, offset)
#=
    # Control points
    controlPoint(map, box)

    # Populate map.xi
    calcParametricMappingLinear(map, box, mesh, )

    @test ( map.ndim )== 2
    @test ( map.nctl )== [4,4,2]
    @test ( map.order )== [4,4,2]
    @test ( map.edge_knot[1] )== [0.0,0.0,0.0,0.0,1.0,1.0,1.0,1.0]
    @test ( map.edge_knot[2] )== [0.0,0.0,0.0,0.0,1.0,1.0,1.0,1.0]
    @test ( map.edge_knot[3] )== [0.0,0.0,1.0,1.0]
    @test ( size(map.aj) )== (3,4,3)
    @test ( size(map.dl) )== (3,3)
    @test ( size(map.dr) )== (3,3)
    @test ( size(map.work) )== (12,4,4,2)
    @test ( size(map.cp_xyz) )== (3,4,4,2)
    @test ( size(map.xi) )== (3,3,1498)
=#
#=
    # Check Control Point Coordinates
    outname = string("./testvalues/control_points_full_linear_mapping.dat")
    orig_control_pts = readdlm(outname)
    for i = 1:length(map.cp_xyz)
      err = norm(map.cp_xyz[i] - orig_control_pts[i], 2)
      @test  err  < 1e-14
    end
=#
#=
    # Check map.xi
    outname = string("./testvalues/xi_full_linear_mapping_DG.dat")
    orig_xi_vals = readdlm(outname)
    for i = 1:10:length(map.xi)
      err = norm(orig_xi_vals[i] - map.xi[i], 2)
      @test  err  < 1e-14
    end
=#
    # Bounding Box Base.Test
    @test ( box.ndim )== 3
    @test ( box.offset )== [0.0,0.0,0.5]
    @test ( box.origin )== [-0.1,-0.1,-0.5]
    @test ( box.unitVector )== [1.0 0.0 0.0
                              0.0 1.0 0.0
                              0.0 0.0 1.0]
    @test ( box.geom_bounds )== [-0.1 -0.1 -0.0
                                0.2 0.1 0.0]
    @test ( box.box_bounds )== [-0.1 -0.1 -0.5
                               0.2 0.1 0.5]
    @test isapprox( box.lengths[1], 0.3) atol= 1e-14
    @test isapprox( box.lengths[2], 0.2) atol= 1e-14
    @test isapprox( box.lengths[3], 1.0) atol= 1e-14
  end # End context("--- Checking Linear mapping For Entire DG Mesh ---")
end  # end facts block
#=
  # Free Form deformation parameters
  ndim = 2
  order = [4,4,2]  # Order of B-splines in the 3 directions
  nControlPts = [4,4,2]

  # Create Mapping object
  map = PumiMapping{Tmsh}(ndim, order, nControlPts, mesh)

  # Create knot vector
  calcKnot(map)

  # Create Bounding box
  offset = [0., 0., 0.5] # No offset in the X & Y direction
  box = PumiBoundingBox{Tmsh}(map, mesh, sbp, offset)

  # Control points
  controlPoint(map, box)

  # Populate map.xi
  calcParametricMappingNonlinear(map, box, mesh)

  @testset "--- Checking Nonlinear Mapping For Entire DG Mesh ---" begin
#=
    # Check if the values are correct
    outname = string("./testvalues/xi_full_linear_mapping_DG.dat")
    orig_xi_vals = readdlm(outname)
    for i = 1:10:length(map.xi)
      err = norm(orig_xi_vals[i] - map.xi[i], 2)
      @test  err  < 1e-14
    end
=#
  end # End context("--- Checking Nonlinear mapping For Entire DG Mesh ---")

  @testset "--- Checking Control Point Manipulation ---" begin
#=
    # Rigid body rotation
    theta = -20*pi/180  # Rotate wall coordinates by 10 degrees
    rotMat = [cos(theta) -sin(theta) 0
              sin(theta) cos(theta)  0
              0          0           1] # Rotation matrix
    # Rotate the control points
    for k = 1:map.nctl[3]
      for j = 1:map.nctl[2]
        for i = 1:map.nctl[1]
          map.cp_xyz[:,i,j,k] = rotMat*map.cp_xyz[:,i,j,k]
        end
      end
    end

    # Rigid body translation
    map.cp_xyz[1,:,:,:] += 0.2
    map.cp_xyz[2,:,:,:] += 0.3

    vertices = evalVolume(map, mesh)
    commitToPumi(map, mesh, sbp, vertices, opts)

    outname = string("./testvalues/volume_coords_full_DG_mesh_2D_airfoil.dat")
    test_vert_coords = readdlm(outname)
    @test ( length(test_vert_coords) )== length(mesh.vert_coords)
    for i = 1:20:length(mesh.vert_coords)
      err = norm(test_vert_coords[i] - mesh.vert_coords[i], 2)
      @test  err  < 1e-14
    end

    writeVisFiles(mesh, "translation_plus_rotation_DG")
=#
  end

end # End facts("--- Checking FFD Types and Functions For Serial DG Pumi Meshes ---")

@testset "--- Checking Specific Geometry Faces in Pumi DG Mesh Embedded in FFD ---" begin

  comm = MPI.COMM_WORLD
  comm_world = MPI.MPI_COMM_WORLD
  comm_self = MPI.COMM_SELF
  my_rank = MPI.Comm_rank(comm)
  comm_size = MPI.Comm_size(comm)

  opts = Dict{ASCIIString, Any}()
  opts["order"] = 1
  opts["dimensions"] = 2
  opts["use_DG"] = true
  opts["Tsbp"] = Float64
  opts["Tmsh"] = Complex128
  opts["operator_type"] = "SBPOmega"
  opts["dmg_name"] = "../src/mesh_files/2D_Airfoil.dmg"
  opts["smb_name"] = "../src/mesh_files/2D_Airfoil.smb"
  opts["numBC"] = 2

  # For 2DAirfoil
  opts["BC1"] = [8,11,14,17]
  opts["BC1_name"] = "FarField"
  opts["BC2"] = [5]
  opts["BC2_name"] = "Airfoil"

  opts["coloring_distance"] = 2 # 0 For CG Mesh 2 for DG Mesh
  opts["jac_type"] = 2
  opts["jac_method"] = 2
  opts["run_type"] = 5

  # Create PumiMesh and SBP objects
  dofpernode = 1
  ref_verts = [-1. 1 -1; -1 -1 1]
  Tsbp = opts["Tsbp"]
  Tmsh = opts["Tmsh"]
  sbp = getTriSBPOmega(degree=opts["order"], Tsbp=opts["Tsbp"])
  sbpface = TriFace{Tsbp}(opts["order"], sbp.cub, ref_verts.')
  topo = 0
  shape_type = 2
  mesh = PumiMeshDG2(Tmsh, sbp, opts, sbpface, dofpernode=dofpernode,
                     shape_type=shape_type)

#  mesh = PumiMeshDG2{Tmsh}(opts["dmg_name"], opts["smb_name"], opts["order"],
#                           sbp, opts, sbpface; dofpernode=dofpernode,
#                           coloring_distance=opts["coloring_distance"],
#                           shape_type=shape_type)

  @testset "--- Checking Linear Mapping for 2D DG Mesh ---" begin

    # geometry faces to be embedded in FFD Box
#    geom_faces = opts["BC2"]
    bc_nums =[2]

    # Free Form deformation parameters
    ndim = 2
    order = [4,4,2]  # Order of B-splines in the 3 directions
    nControlPts = [4,4,2]

    # Create Mapping object
    map = PumiMapping{Tmsh}(ndim, order, nControlPts, mesh, full_geom=false, bc_nums=bc_nums)

    # Create knot vector
    calcKnot(map)

    # Create Bounding box
    offset = [0., 0., 0.5] # No offset in the X & Y direction
    box = PumiBoundingBox{Tmsh}(map, mesh, sbp, offset)

    # Control points
    controlPoint(map, box)
#=
    calcParametricMappingLinear(map, box, mesh, bc_nums)

    @test ( map.ndim )== 2
    @test ( map.full_geom )== false
    @test ( map.nctl )== [4,4,2]
    @test ( map.order )== [4,4,2]
    @test ( map.edge_knot[1] )== [0.0,0.0,0.0,0.0,1.0,1.0,1.0,1.0]
    @test ( map.edge_knot[2] )== [0.0,0.0,0.0,0.0,1.0,1.0,1.0,1.0]
    @test ( map.edge_knot[3] )== [0.0,0.0,1.0,1.0]
    @test ( size(map.aj) )== (3,4,3)
    @test ( size(map.dl) )== (3,3)
    @test ( size(map.dr) )== (3,3)
    @test ( size(map.work) )== (12,4,4,2)
    @test ( size(map.cp_xyz) )== (3,4,4,2)
    @test ( size(map.xi) )== (1,)
    @test ( size(map.xi[1]) )== (3,2,102) # Essentially tests defineMapXi
=#
#=
    outname = string("./testvalues/xi_values_2D_airfoil_face5.dat")
    ctr = 1
    test_xi_values = readdlm(outname)
    for itr = 1:length(geom_faces)
      geom_face_number = geom_faces[itr]
      # get the boundary array associated with the geometric edge
      itr2 = 0
      for itr2 = 1:mesh.numBC
        if findfirst(mesh.bndry_geo_nums[itr2],geom_face_number) > 0
          break
        end
      end
      start_index = mesh.bndry_offsets[itr2]
      end_index = mesh.bndry_offsets[itr2+1]
      idx_range = start_index:end_index
      bndry_facenums =sview(mesh.bndryfaces, start_index:(end_index - 1))
      nfaces = length(bndry_facenums)
      for i = 1:nfaces
        bndry_i = bndry_facenums[i]
        # get the local index of the vertices
        vtx_arr = mesh.topo.face_verts[:,bndry_i.face]
        for j = 1:length(vtx_arr)
          err1 = norm(map.xi[itr][1,j,i] - test_xi_values[ctr],2)
          err2 = norm(map.xi[itr][2,j,i] - test_xi_values[ctr+1],2)
          @test  err1  < 1e-14
          @test  err2  < 1e-14
          ctr += 2
        end  # End for j = 1:length(vtx_arr)
      end    # End for i = 1:nfaces
    end      # End for itr = 1:length(geomfaces)
=#

  end # End context("--- Checking Linear Mapping for DG Mesh ---")

  # geometry faces to be embedded in FFD Box
#  geom_faces = opts["BC2"]
  bc_nums =[2]

  # Free Form deformation parameters
  ndim = 2
  order = [4,4,2]  # Order of B-splines in the 3 directions
  nControlPts = [4,4,2]

  # Create Mapping object
  map = PumiMapping{Tmsh}(ndim, order, nControlPts, mesh, full_geom=false, bc_nums=bc_nums)

  # Create knot vector
  calcKnot(map)

  # Create Bounding box
  offset = [0., 0., 0.5] # No offset in the X & Y direction
  box = PumiBoundingBox{Tmsh}(map, mesh, sbp, offset)

  # Control points
  controlPoint(map, box)

  @testset "--- Checking Nonlinear Mapping for DG Mesh ---" begin
#=
    # Populate map.xi
    calcParametricMappingNonlinear(map, box, mesh, geom_faces)
    @test ( size(map.xi) )== (1,)
    @test ( size(map.xi[1]) )== (3,2,102) # Essentially tests defineMapXi

    outname = string("./testvalues/xi_values_2D_airfoil_face5.dat")
    ctr = 1
    test_xi_values = readdlm(outname)
    for itr = 1:length(geom_faces)
      geom_face_number = geom_faces[itr]
      # get the boundary array associated with the geometric edge
      itr2 = 0
      for itr2 = 1:mesh.numBC
        if findfirst(mesh.bndry_geo_nums[itr2],geom_face_number) > 0
          break
        end
      end
      start_index = mesh.bndry_offsets[itr2]
      end_index = mesh.bndry_offsets[itr2+1]
      idx_range = start_index:end_index
      bndry_facenums =sview(mesh.bndryfaces, start_index:(end_index - 1))
      nfaces = length(bndry_facenums)
      for i = 1:nfaces
        bndry_i = bndry_facenums[i]
        # get the local index of the vertices
        vtx_arr = mesh.topo.face_verts[:,bndry_i.face]
        for j = 1:length(vtx_arr)
          err1 = norm(map.xi[itr][1,j,i] - test_xi_values[ctr],2)
          err2 = norm(map.xi[itr][2,j,i] - test_xi_values[ctr+1],2)
          @test  err1  < 1e-14
          @test  err2  < 1e-14
          ctr += 2
        end  # End for j = 1:length(vtx_arr)
      end    # End for i = 1:nfaces
    end      # End for itr = 1:length(geomfaces)
=#
  end # End context("--- Checking Nonlinear Mapping for DG Mesh ---")

  @testset "--- Checking Surface evaluation for 2D DG Mesh ---" begin

#=
#    commitToPumi(map, mesh, sbp, vertices, opts)
#    writeVisFiles(mesh, "mesh_warped")

    outname = string("./testvalues/modified_coordinates_2D_airfoil_face5.dat")
    test_surface_coords = readdlm(outname)
    ctr = 1
    for itr = 1:length(map.geom_faces)
      geom_face_number = map.geom_faces[itr]
      # get the boundary array associated with the geometric edge
      itr2 = 0
      for itr2 = 1:mesh.numBC
        if findfirst(mesh.bndry_geo_nums[itr2],geom_face_number) > 0
          break
        end
      end
      start_index = mesh.bndry_offsets[itr2]
      end_index = mesh.bndry_offsets[itr2+1]
      idx_range = start_index:(end_index-1)
      bndry_facenums =sview(mesh.bndryfaces, idx_range) # faces on geometric edge i
      nfaces = length(bndry_facenums)
      for i = 1:nfaces
        bndry_i = bndry_facenums[i]
        # get the local index of the vertices on the boundary face (local face number)
        vtx_arr = mesh.topo.face_verts[:,bndry_i.face]
        for j = 1:length(vtx_arr)
          err1 = norm(mesh.vert_coords[1,vtx_arr[j],bndry_i.element] -
                 test_surface_coords[ctr], 2)
          err2 = norm(mesh.vert_coords[2,vtx_arr[j],bndry_i.element] -
                 test_surface_coords[ctr+1], 2)
          @test  err1  < 1e-14
          @test  err2  < 1e-14
          ctr += 2
        end  # End for j = 1:length(vtx_arr)
      end    # End for i = 1:nfaces
    end  # End for itr = 1:length(map.geom_faces)
=#
  end # End context("--- Checking Surface evaluation for 2D DG Mesh ---")
#=
  @testset "--- Checking evaldXdControlPointProduct for 2D DG Mesh ---" begin

    fill!(map.work, 0.0)
    pert = 1e-6

    # Create seed vector
    # - Get original wall coordinates
    orig_wallCoords = evalSurface(map, mesh)[1]  # only do the first geometric entity
#    orig_wallCoords = FFD.getUniqueWallCoordsArray(mesh, geom_faces)
    nwall_faces = FFD.getnWallFaces(mesh, geom_faces)
    Xs_bar = randn(3, size(orig_wallCoords,2))
    Xs_bar[3,:] = 0.0 # To accurately simulate a 2D mesh
    cp_xyz_bar = zeros(map.cp_xyz)
    evaldXdControlPointProduct(map, mesh, vec(Xs_bar))
    for i = 1:size(map.work, 4)
      for j = 1:size(map.work, 3)
        for k = 1:size(map.work, 2)
            cp_xyz_bar[1:3,k,j,i] = map.work[1:3, k,j,i]
        end
      end
    end

    # Check against finite difference
    cp_jacobian = zeros(length(orig_wallCoords), length(map.cp_xyz))
    for i = 1:length(map.cp_xyz)
      map.cp_xyz[i] += pert
      vertices = evalSurface(map, mesh)
      commitToPumi(map, mesh, sbp, vertices, opts)
      new_wallCoords = FFD.getUniqueWallCoordsArray(mesh, geom_faces)
      cp_jacobian[:,i] = (vec(new_wallCoords) - vec(orig_wallCoords))/pert
      map.cp_xyz[i] -= pert
    end # End for i = 1:length(map.cp_xyz)
    prod_val = transpose(cp_jacobian)*vec(Xs_bar)

    error = vec(cp_xyz_bar) - prod_val
    for i = 1:length(error)
      @test isapprox( error[i], 0.0) atol=1e-9
    end

  end # End context("--- Checking evaldXdControlPointProduct for 2D DG Mesh ---")
=#
  @testset "---- Checking Transposed Jacobian-vector Product -----" begin
#=
    delta_B = zeros(eltype(map.work), 3, map.nctl[1], map.nctl[2], map.nctl[3])
    Xs_orig = evalSurface(map, mesh)
    println("typeof(Xs_orig) = ", typeof(Xs_orig))
    println("size(Xs_orig) = ", size(Xs_orig))
    println("size(Xs_orig[1]) = ", size(Xs_orig[1]))

    delta_S = Array{Array{Complex128, 3}}(length(map.geom_faces))

    h = 1e-6

    for itr=1:length(map.geom_faces)
      println("itr = ", itr)
      delta_S[itr] = zeros(Complex128, size(Xs_orig[itr]))

      # compute finite differences delta_s^T jac * delta_ba
      # and compare against the AD version
      # here we take all combinations of delta_S and delta_x = 1 at a single
      # entry
      for i=4:length(map.cp_xyz)
        # in 2D, skip z components
        if (i % 3) == 0
          continue
        end
        println("i = ", i)
        map.cp_xyz[i] += h
        verts_new = evalSurface(map, mesh)
        map.cp_xyz[i] -= h

        for j=1:length(verts_new[itr])
#          if (j % 3) == 0
#            continue
#          end
          println("j = ", j)
          fill!(delta_B, 0.0)
          delta_S[itr][j] = 1
          FFD.evaldXdControlPointTransposeProduct(map, mesh, delta_S, delta_B)

          delta_xj = (verts_new[itr][j] - Xs_orig[itr][j])/h
          println("delta_xj = ", delta_xj)
          println("delta_B = ", delta_B[i])

          @test isapprox( delta_B[i], delta_xj) atol=1e-6
          delta_S[itr][j] = 0
        end
      end
    end
=#

      end  # end context

end # End facts("--- Checking Specific Geometry Faces in Pumi DG Mesh Embedded in FFD ---")

#=
@testset "--- Checking Functions Specific to CG Pumi Meshes in Serial ---" begin

  # MPI Declarations
  comm = MPI.COMM_WORLD
  comm_world = MPI.MPI_COMM_WORLD
  comm_self = MPI.COMM_SELF
  my_rank = MPI.Comm_rank(comm)
  comm_size = MPI.Comm_size(comm)

  opts = Dict{ASCIIString, Any}()
  opts["order"] = 1
  opts["dimensions"] = 2
  opts["use_DG"] = true
  opts["Tsbp"] = Float64
  opts["Tmsh"] = Complex128
  opts["operator_type"] = "SBPGamma"
  opts["dmg_name"] = "../src/mesh_files/2D_Airfoil.dmg"
  opts["smb_name"] = "../src/mesh_files/2D_Airfoil.smb"
  opts["numBC"] = 2

  # For 2DAirfoil
  opts["BC1"] = [8,11,14,17]
  opts["BC1_name"] = "FarField"
  opts["BC2"] = [5]
  opts["BC2_name"] = "Airfoil"

  opts["coloring_distance"] = 0 # For CG Mesh 2 for DG Mesh
  opts["jac_type"] = 2
  opts["jac_method"] = 2
  opts["run_type"] = 5
  opts["use_edge_res"] = false
  opts["write_edge_vertnums"] = false
  opts["write_face_vertnums"] = false

  # Create PumiMesh and SBP objects
  dofpernode = 1
  ref_verts = [-1. 1 -1; -1 -1 1]
  Tsbp = opts["Tsbp"]
  Tmsh = opts["Tmsh"]

  shape_type = 1
  sbp = getTriSBPGamma(degree=opts["order"], Tsbp=Tsbp)
  sbpface = TriFace{Float64}(opts["order"], sbp.cub, sbp.vtx)
  mesh = PumiMesh2{Tmsh}(opts["dmg_name"], opts["smb_name"], opts["order"], sbp, opts, sbpface;
                             dofpernode=dofpernode,
                             coloring_distance=opts["coloring_distance"],
                             shape_type=shape_type)

  @testset "--- Checking Linear Mapping For Entire CG Mesh ---" begin

    # Free Form deformation parameters
    ndim = 2
    order = [4,4,2]  # Order of B-splines in the 3 directions
    nControlPts = [4,4,2]

    # Create Mapping object
    map = PumiMapping{Tmsh}(ndim, order, nControlPts, mesh)

    # Create knot vector
    calcKnot(map) # No need to test this since doesnt depend on mesh type

    # Create Bounding box
    # Doesn't depend on mesh type if entire geometry is embedded in FFD box
    offset = [0., 0., 0.5] # No offset in the X & Y direction
    box = PumiBoundingBox{Tmsh}(map, mesh, sbp, offset)

    # Control points
    controlPoint(map, box) # No Need to test this since it doesnt depend on mesh type

    # Populate map.xi
    calcParametricMappingLinear(map, box, mesh)

    # Check map.xi
    outname = string("./testvalues/xi_full_linear_mapping_DG.dat")
    orig_xi_vals = readdlm(outname)
    for i = 1:10:length(map.xi)
      err = norm(orig_xi_vals[i] - map.xi[i], 2)
      @test  err  < 1e-14
    end

  end # End context("--- Checking Linear mapping For Entire DG Mesh ---")

  @testset "--- Checking Nonlinear Mapping For Entire CG Mesh ---" begin

    # Free Form deformation parameters
    ndim = 2
    order = [4,4,2]  # Order of B-splines in the 3 directions
    nControlPts = [4,4,2]

    # Create Mapping object
    map = PumiMapping{Tmsh}(ndim, order, nControlPts, mesh)

    # Create knot vector
    calcKnot(map) # No need to test this since doesnt depend on mesh type

    # Create Bounding box
    # Doesn't depend on mesh type if entire geometry is embedded in FFD box
    offset = [0., 0., 0.5] # No offset in the X & Y direction
    box = PumiBoundingBox{Tmsh}(map, mesh, sbp, offset)

    # Control points
    controlPoint(map, box) # No Need to test this since it doesnt depend on mesh type

    # Populate map.xi
    calcParametricMappingNonlinear(map, box, mesh)

    # Check map.xi
    outname = string("./testvalues/xi_full_linear_mapping_DG.dat")
    orig_xi_vals = readdlm(outname)
    for i = 1:10:length(map.xi)
      err = norm(orig_xi_vals[i] - map.xi[i], 2)
      @test  err  < 1e-14
    end

  end # End context("--- Checking Linear mapping For Entire DG Mesh ---")

end # End facts("--- Checking Functions Specific to CG Pumi Meshes in Serial ---")

@testset "--- Checking Specific Geometry Faces in Pumi CG Mesh Embedded in FFD ---" begin

  comm = MPI.COMM_WORLD
  comm_world = MPI.MPI_COMM_WORLD
  comm_self = MPI.COMM_SELF
  my_rank = MPI.Comm_rank(comm)
  comm_size = MPI.Comm_size(comm)

  opts = PdePumiInterface.get_defaults()

  # 2D mesh
  opts["order"] = 1
  opts["dimensions"] = 2
  opts["use_DG"] = false
  opts["operator_type"] = "SBPGamma"
  opts["dmg_name"] = "../src/mesh_files/2D_Airfoil.dmg"
  opts["smb_name"] = "../src/mesh_files/2D_Airfoil.smb"
  opts["numBC"] = 2

  # For 2DAirfoil
  opts["BC1"] = [8,11,14,17]
  opts["BC1_name"] = "FarField"
  opts["BC2"] = [5]
  opts["BC2_name"] = "Airfoil"

  opts["coloring_distance"] = 0 # 0 For CG Mesh 2 for DG Mesh
  opts["jac_type"] = 2
  opts["jac_method"] = 2
  opts["run_type"] = 5

  # Create PumiMesh and SBP objects
  sbp, mesh, pmesh, Tsol, Tres, Tmsh, mesh_time = createMeshAndOperator(opts, 1)

  @testset "--- Checking Linear Mapping for 2D CG Mesh ---" begin

    # geometry faces to be embedded in FFD Box
    geom_faces = opts["BC2"]

    # Free Form deformation parameters
    ndim = 2
    order = [4,4,2]  # Order of B-splines in the 3 directions
    nControlPts = [4,4,2]

    # Create Mapping object
    map = PumiMapping{Tmsh}(ndim, order, nControlPts, mesh, full_geom=false, geom_faces=geom_faces)

    # Create knot vector
    calcKnot(map)

    # Create Bounding box
    offset = [0., 0., 0.5] # No offset in the X & Y direction
    box = PumiBoundingBox{Tmsh}(map, mesh, sbp, offset)

    # Control points
    controlPoint(map, box)

    calcParametricMappingLinear(map, box, mesh, geom_faces)

    @test ( map.ndim )== 2
    @test ( map.full_geom )== false
    @test ( map.nctl )== [4,4,2]
    @test ( map.order )== [4,4,2]
    @test ( map.edge_knot[1] )== [0.0,0.0,0.0,0.0,1.0,1.0,1.0,1.0]
    @test ( map.edge_knot[2] )== [0.0,0.0,0.0,0.0,1.0,1.0,1.0,1.0]
    @test ( map.edge_knot[3] )== [0.0,0.0,1.0,1.0]
    @test ( size(map.aj) )== (3,4,3)
    @test ( size(map.dl) )== (3,3)
    @test ( size(map.dr) )== (3,3)
    @test ( size(map.work) )== (12,4,4,2)
    @test ( size(map.cp_xyz) )== (3,4,4,2)

    outname = string("./testvalues/xi_values_2D_airfoil_face5.dat")
    test_xi_values = readdlm(outname)
    ctr = 1
    for itr = 1:length(geom_faces)
      geom_face_number = geom_faces[itr]
      # get the boundary array associated with the geometric edge
      itr2 = 0
      for itr2 = 1:mesh.numBC
        if findfirst(mesh.bndry_geo_nums[itr2],geom_face_number) > 0
          break
        end
      end
      start_index = mesh.bndry_offsets[itr2]
      end_index = mesh.bndry_offsets[itr2+1]
      idx_range = start_index:end_index
      bndry_facenums =sview(mesh.bndryfaces, start_index:(end_index - 1))
      nfaces = length(bndry_facenums)
      for i = 1:nfaces
        bndry_i = bndry_facenums[i]
        # get the local index of the vertices
        vtx_arr = mesh.topo.face_verts[:,bndry_i.face]
        for j = 1:length(vtx_arr)
          err1 = norm(map.xi[itr][1,j,i] - test_xi_values[ctr],2)
          err2 = norm(map.xi[itr][2,j,i] - test_xi_values[ctr+1],2)
          @test  err1  < 1e-14
          @test  err2  < 1e-14
          ctr += 2
        end  # End for j = 1:length(vtx_arr)
      end    # End for i = 1:nfaces
    end      # End for itr = 1:length(geomfaces)

  end # End context("--- Checking Linear Mapping for DG Mesh ---")

  # geometry faces to be embedded in FFD Box
  geom_faces = opts["BC2"]

  # Free Form deformation parameters
  ndim = 2
  order = [4,4,2]  # Order of B-splines in the 3 directions
  nControlPts = [4,4,2]

  # Create Mapping object
  map = PumiMapping{Tmsh}(ndim, order, nControlPts, mesh, full_geom=false, geom_faces=geom_faces)

  # Create knot vector
  calcKnot(map)

  # Create Bounding box
  offset = [0., 0., 0.5] # No offset in the X & Y direction
  box = PumiBoundingBox{Tmsh}(map, mesh, sbp, offset)

  # Control points
  controlPoint(map, box)

  @testset "--- Checking Nonlinear Mapping for CG Mesh ---" begin

    # Populate map.xi
    calcParametricMappingNonlinear(map, box, mesh, geom_faces)

    # Check for errors in map.xi
    outname = string("./testvalues/xi_values_2D_airfoil_face5.dat")
    test_xi_values = readdlm(outname)
    ctr = 1
    for itr = 1:length(geom_faces)
      geom_face_number = geom_faces[itr]
      # get the boundary array associated with the geometric edge
      itr2 = 0
      for itr2 = 1:mesh.numBC
        if findfirst(mesh.bndry_geo_nums[itr2],geom_face_number) > 0
          break
        end
      end
      start_index = mesh.bndry_offsets[itr2]
      end_index = mesh.bndry_offsets[itr2+1]
      idx_range = start_index:end_index
      bndry_facenums =sview(mesh.bndryfaces, start_index:(end_index - 1))
      nfaces = length(bndry_facenums)
      for i = 1:nfaces
        bndry_i = bndry_facenums[i]
        # get the local index of the vertices
        vtx_arr = mesh.topo.face_verts[:,bndry_i.face]
        for j = 1:length(vtx_arr)
          err1 = norm(map.xi[itr][1,j,i] - test_xi_values[ctr],2)
          err2 = norm(map.xi[itr][2,j,i] - test_xi_values[ctr+1],2)
          @test  err1  < 1e-14
          @test  err2  < 1e-14
          ctr += 2
        end  # End for j = 1:length(vtx_arr)
      end    # End for i = 1:nfaces
    end      # End for itr = 1:length(geomfaces)

  end # End context("--- Checking Nonlinear Mapping for DG Mesh ---")

  @testset "--- Checking Surface evaluation for 2D CG Mesh ---" begin

    evalSurface(map, mesh, sbp)
    outname = string("./testvalues/modified_coordinates_2D_airfoil_face5.dat")
    #=
    test_surface_coords = readdlm(outname)
    ctr = 1
    for itr = 1:length(map.geom_faces)
      geom_face_number = map.geom_faces[itr]
      # get the boundary array associated with the geometric edge
      itr2 = 0
      for itr2 = 1:mesh.numBC
        if findfirst(mesh.bndry_geo_nums[itr2],geom_face_number) > 0
          break
        end
      end
      start_index = mesh.bndry_offsets[itr2]
      end_index = mesh.bndry_offsets[itr2+1]
      idx_range = start_index:(end_index-1)
      bndry_facenums =sview(mesh.bndryfaces, idx_range) # faces on geometric edge i
      nfaces = length(bndry_facenums)
      for i = 1:nfaces
        bndry_i = bndry_facenums[i]
        # get the local index of the vertices on the boundary face (local face number)
        vtx_arr = mesh.topo.face_verts[:,bndry_i.face]
        for j = 1:length(vtx_arr)
          err1 = norm(mesh.coords[1,vtx_arr[j],bndry_i.element] -
                 test_surface_coords[ctr], 2)
          err2 = norm(mesh.coords[2,vtx_arr[j],bndry_i.element] -
                 test_surface_coords[ctr+1], 2)
          @test  err1  < 1e-14
          @test  err2  < 1e-14
          ctr += 2
        end  # End for j = 1:length(vtx_arr)
      end    # End for i = 1:nfaces
    end  # End for itr = 1:length(map.geom_faces)
    =#

  end # End context("--- Checking Nonlinear Mapping for DG Mesh ---")
end # End facts("--- Checking Specific Geometry Faces in Pumi CG Mesh Embedded in FFD ---")
=#
=#


text:
module Bioinformatics195

export normalizeDNA,
       composition,
       gc_content,
       complement,
       reverse_complement,
       parse_fasta

# # uncomment the following line if you intend to use BioSequences types
#using BioSequences

"""
    normalizeDNA(::AbstractString)

Ensures that a sequence only contains valid bases
(or `'N'` for unknown bases).
Returns a String.
"""
function normalizeDNA(seq)
    seq = uppercase(string(seq))
    for base in seq
        # note: `N` indicates an unknown base
        occursin(base, "AGCTN") || error("invalid base $base")
    end
    return seq # change to `return LongDNASeq(seq)` if you want to try to use BioSequences types
end

function composition(seq)
    seq = normalizeDNA(seq)

    bases = Dict()

    for base in seq
        if haskey(bases, base)
            bases[base] += 1
        else
            bases[base] = 1
        end
    end
    return bases
end

function gc_content(seq)
    bases = composition(seq)

    return (get(bases,'G', 0) + get(bases, 'C', 0)) / length(seq)
end

function complementbase(base)
    base = uppercase(string(base))
    comp = Dict("A"=>'T',
                "T"=>'A',
                "G"=>'C',
                "C"=>'G',
                "N"=>'N')
    return comp[uppercase(base)]
end

function complement(seq)
    seq = normalizeDNA(seq)
    result = ""

    for base in seq
        result = result * complementbase(base)
    end
    return result
end

function reverse_complement(seq)
    seq = normalizeDNA(seq)
    result = ""
    index = length(seq)

    while( index > 0)
        result = result * complementbase(seq[index])
        index = index - 1
    end
    return result
end

function parse_fasta(path)
    heading = Vector()
    seqs = Vector()
    sequence = ""

    for line in eachline(path)
        if startswith(line, '>')
            if length(sequence) != 0
                push!(seqs, sequence)
            end
            sequence = ""
            push!(heading, line[2:end])
        else
            line = normalizeDNA(line)
            sequence = sequence * line
        end
    end
    if length(sequence) != 0
        push!(seqs, sequence)
    end

    return heading, seqs
end

end # module Bioinformatics195


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

# This file contains analysis of constant values.
# So far, we assume there is no alias between different symbols.
# TODO: when alias is considered, we need add a post-processing that
# is a contamination process, to remove potential variants, due to aliases.
# But there is no need to change the basic analysis.

@doc """ Find symbols/GenSyms whose values are constant in the region. """
function find_constant_values(
    region   :: Region,
    liveness :: Liveness,
    cfg      :: CFG
)
    constants = Set{Sym}()
    blocks    = cfg.basic_blocks

    if isa(region, LoopRegion)
        bb_idxs = region.loop.members
    else
        bb_idxs = keys(blocks)
    end

    # Add all the variables that are read in the region
    for bb_idx in bb_idxs
        bb = blocks[bb_idx]
        for stmt in bb.statements
            use = LivenessAnalysis.use(stmt, liveness)
            union!(constants, use)
        end
    end

    # Remove all the variables that are written in the region
    for bb_idx in bb_idxs
        bb = blocks[bb_idx]
        for stmt in bb.statements
            def = LivenessAnalysis.def(stmt, liveness)
            setdiff!(constants, def)
        end
    end

    dprintln(1, 0, "\nConstants discovered:")
    dprintln(1, 1, "", constants)

    constants
end

@doc """
Find symbols/GenSyms that are statically defined exactly once in the region.
They are not constants, but finding them is similar to finding constants.
"""
function find_single_defs(
    region   :: Region,
    liveness :: Liveness,
    cfg      :: CFG
)
    single_defs = Set{Sym}()
    blocks      = cfg.basic_blocks

    # ISSUE: What is exactly the definition and usage of single-def?
    # Should the initial values of single_defs include the symbols live into the
    # region, even if they are not defined in the region?
    if isa(region, LoopRegion)
        bb_idxs     = region.loop.members
#        head_bb     = blocks[region.loop.head]
#        single_defs = LivenessAnalysis.live_in(head_bb, liveness)
    else
        bb_idxs     = keys(blocks)
#        entry_bb    = blocks[-1]
#        single_defs = LivenessAnalysis.live_in(entry_bb, liveness)
    end

    for bb_idx in bb_idxs
        bb = blocks[bb_idx]
        for stmt in bb.statements
            def            = LivenessAnalysis.def(stmt, liveness)
            total          = union(single_defs, def)
            double_defined = intersect(single_defs, def)
            single_defs    = setdiff(total, double_defined)
        end
    end

    dprintln(1, 0, "\nSingle-defs discovered:")
    dprintln(1, 1, "", single_defs)

    single_defs
end

max_stars_repo_name: IntelLabs/Sparso

max_issues_repo_name: IntelLabs/Sparso

max_forks_repo_name: IntelLabs/Sparso


text:
module FEMMPointMassModule

using LinearAlgebra: norm, Transpose, mul!
using FinEtools
using FinEtools.IntegDomainModule: IntegDomain
import FinEtoolsDeforLinear.MatDeforElastIsoModule: MatDeforElastIso
using ..FESetCorotBeamModule: FESetL2CorotBeam, local_frame_and_def!, local_mass!, local_stiffness!, natural_forces!, local_geometric_stiffness!, local_forces!, MASS_TYPE_CONSISTENT_WITH_ROTATION_INERTIA

"""
    FEMMPointMass{S<:AbstractFESet}

Class for linear deformation finite element modeling machine.
"""
mutable struct FEMMPointMass{S<:AbstractFESet} <: AbstractFEMM
    integdomain::IntegDomain{S} # integration domain data
    _massmatrix::FFltMat
end

"""
    mass(self::FEMMCorotBeam,  assembler::A,
      geom::NodalField{FFlt},
      u::NodalField{T}) where {A<:AbstractSysmatAssembler, T<:Number}

Compute the consistent mass matrix

This is a general routine for the abstract linear-deformation  FEMM.
"""
function mass(self::FEMMPointMass, assembler::ASS, geom0::NodalField{FFlt}, u1::NodalField{T}, Rfield1::NodalField{T}, dchi::NodalField{TI}) where {ASS<:AbstractSysmatAssembler, T<:Number, TI<:Number}
    fes = self.integdomain.fes
    dofnums = zeros(FInt, 1, 6);
    startassembly!(assembler, 3, 3, count(fes), dchi.nfreedofs, dchi.nfreedofs);
    for i = 1:count(fes) # Loop over elements
        gatherdofnums!(dchi, dofnums, fes.conn[i]); # degrees of freedom
        assemble!(assembler, self._massmatrix, dofnums[1:3], dofnums[1:3]);
    end # Loop over elements
    return makematrix!(assembler);
end

function mass(self::FEMMPointMass, geom0::NodalField{FFlt}, u1::NodalField{T}, Rfield1::NodalField{T}, dchi::NodalField{TI}) where {T<:Number, TI<:Number}
    assembler = SysmatAssemblerSparseSymm();
    return mass(self, assembler, geom0, u1, Rfield1, dchi);
end

end # module


text:
mutable struct Bard1 <:AbstractMPCCModel
  meta :: MPCCModelMeta
  counters :: MPCCCounters
end

function Bard1(; x0 = zeros(5))

  nvar = length(x0)
  lvar = -Inf*ones(nvar)
  uvar =  Inf*ones(nvar)
  ncon = 1
  lcon = [0.]
  ucon = [0.]
  y0   = [0.]
  ncc  = 3
  lccG = zeros(3)
  lccH = zeros(3)

  nnzh = nvar * (nvar + 1) / 2
  nnzj = nvar * ncon

  meta = MPCCModelMeta(nvar, x0 = x0, y0=y0, ncc = ncc, lccG = lccG, lccH = lccH,
                             lvar = lvar, uvar = uvar,
                             ncon = ncon, lcon = lcon, ucon = ucon,
                             nnzj=nnzj, nnzh=nnzh,
                             minimize=true, islp=false, name="Bard1")

  return Bard1(meta, MPCCCounters())
end

function MPCC.obj(nlp :: Bard1, x :: AbstractVector)
  #increment!(nlp, :neval_obj)
  return (x[1]-5)^2+(2*x[2]+1)^2
end

function MPCC.grad!(nlp :: Bard1, x :: AbstractVector, gx :: AbstractVector)
  #increment!(nlp, :neval_grad)
  gx .= vcat(2*(x[1] - 5), 4 * (2 * x[2] + 1), zeros(3))
  return gx
end
function MPCC.cons_nl!(nlp :: Bard1, x :: AbstractVector, c :: AbstractVector)
  #increment!(nlp, :neval_cons)
  c[1:nlp.meta.ncon] .= [2*(x[2]-1)-1.5*x[1]+x[3]-0.5*x[4]+x[5]]
  return c
end

function MPCC.consG!(nlp :: Bard1, x :: AbstractVector, c :: AbstractVector)
  #increment!(nlp, :neval_consG)
  c[1:nlp.meta.ncc] .= vcat(3x[1]-x[2]-3, -x[1]+0.5x[2]+4, -x[1]-x[2]+7 )
  return c
end

function MPCC.consH!(nlp :: Bard1, x :: AbstractVector, c :: AbstractVector)
  #increment!(nlp, :neval_consH)
  c[1:nlp.meta.ncc] .= vcat(x[3], x[4], x[5])
  return c
end

function MPCC.jac_nl_structure!(nlp :: Bard1, rows :: AbstractVector{<: Integer}, cols :: AbstractVector{<: Integer})
  m, n = nlp.meta.ncon, nlp.meta.nvar
  I = ((i,j) for i = 1:m, j = 1:n)
  rows[1 : nlp.meta.nnzj] .= getindex.(I, 1)[:]
  cols[1 : nlp.meta.nnzj] .= getindex.(I, 2)[:]
  return rows, cols
end

function MPCC.jac_nl_coord!(nlp :: Bard1, x :: AbstractVector, vals :: AbstractVector)
  #increment!(nlp, :neval_jac)
  Jx = [-1.5 2. 1. -0.5 1.]
  vals[1 : nlp.meta.nnzj] .= Jx[:]
  return vals
end

function MPCC.jacG_structure!(nlp :: Bard1, rows :: AbstractVector{<: Integer}, cols :: AbstractVector{<: Integer})
  m, n = nlp.meta.ncc, nlp.meta.nvar
  I = ((i,j) for i = 1:m, j = 1:n)
  rows[1 : n*m] .= getindex.(I, 1)[:]
  cols[1 : n*m] .= getindex.(I, 2)[:]
  return rows, cols
end

function MPCC.jacG_coord!(nlp :: Bard1, x :: AbstractVector, vals :: AbstractVector)
  m, n = nlp.meta.ncc, nlp.meta.nvar
  #increment!(nlp, :neval_jacG)
  Jx = vcat([3 -1 0 0 0], [-1 0.5 0 0 0], [-1 -1 0 0 0])
  vals[1 : n*m] .= Jx[:]
  return vals
end

function MPCC.jacH_structure!(nlp :: Bard1, rows :: AbstractVector{<: Integer}, cols :: AbstractVector{<: Integer})
  m, n = nlp.meta.ncc, nlp.meta.nvar
  I = ((i,j) for i = 1:m, j = 1:n)
  rows[1 : n*m] .= getindex.(I, 1)[:]
  cols[1 : n*m] .= getindex.(I, 2)[:]
  return rows, cols
end

function MPCC.jacH_coord!(nlp :: Bard1, x :: AbstractVector, vals :: AbstractVector)
  m, n = nlp.meta.ncc, nlp.meta.nvar
  #increment!(nlp, :neval_jacH)
  Jx = zeros(3, 5)
  Jx[1,1] = 1.
  Jx[2,2] = 1.
  Jx[3,3] = 1.
  vals[1 : n*m] .= Jx[:]
  return vals
end

function MPCC.jnlprod!(nlp :: Bard1, x :: AbstractVector, v :: AbstractVector, Jv :: AbstractVector)
  #increment!(nlp, :neval_jprod)
  Jv[1:nlp.meta.ncon] .= [-1.5 2. 1. -0.5 1.] * v
  return Jv
end

function MPCC.jnltprod!(nlp :: Bard1, x :: AbstractVector, v :: AbstractVector, Jtv :: AbstractVector)
  #increment!(nlp, :neval_jtprod)
  Jtv[1:nlp.meta.nvar] .= [-1.5 2. 1. -0.5 1.]' * v
  return Jtv
end

function MPCC.jGprod!(nlp :: Bard1, x :: AbstractVector, v :: AbstractVector, Jv :: AbstractVector)
  #increment!(nlp, :neval_jGprod)
  Jv[1:nlp.meta.ncc] .=  vcat([3 -1 0 0 0], [-1 0.5 0 0 0], [-1 -1 0 0 0]) * v
  return Jv
end

function MPCC.jGtprod!(nlp :: Bard1, x :: AbstractVector, v :: AbstractVector, Jtv :: AbstractVector)
  #increment!(nlp, :neval_jGtprod)
  Jtv[1:nlp.meta.nvar] .=  vcat([3 -1 0 0 0], [-1 0.5 0 0 0], [-1 -1 0 0 0])' * v
  return Jtv
end

function MPCC.jHprod!(nlp :: Bard1, x :: AbstractVector, v :: AbstractVector, Jv :: AbstractVector)
  #increment!(nlp, :neval_jHprod)
  Jx = zeros(3, 5)
  Jx[1,1] = 1.
  Jx[2,2] = 1.
  Jx[3,3] = 1.
  Jv[1:nlp.meta.ncc] .= Jx * v
  return Jv
end

function MPCC.jHtprod!(nlp :: Bard1, x :: AbstractVector, v :: AbstractVector, Jtv :: AbstractVector)
  #increment!(nlp, :neval_jHtprod)
  Jx = zeros(3, 5)
  Jx[1,1] = 1.
  Jx[2,2] = 1.
  Jx[3,3] = 1.
  Jtv[1:nlp.meta.nvar] .= Jx' * v
  return Jtv
end

#function hess(nlp :: Bard1, x :: AbstractVector; obj_weight :: Real = one(eltype(x)))
#  #increment!(nlp, :neval_hess)
#  ℓ(x) = obj_weight * nlp.f(x)
#  Hx = ForwardDiff.hessian(ℓ, x)
#  return tril(Hx)
#end

#function hess(nlp :: Bard1, x :: AbstractVector, y :: AbstractVector; obj_weight :: Real = one(eltype(x)))
#  #increment!(nlp, :neval_hess)
#  ℓ(x) = obj_weight * nlp.f(x) + dot(nlp.c(x), y)
#  Hx = ForwardDiff.hessian(ℓ, x)
#  return tril(Hx)
#end

function MPCC.hess_structure!(nlp :: Bard1, rows :: AbstractVector{<: Integer}, cols :: AbstractVector{<: Integer})
  n = nlp.meta.nvar
  I = ((i,j) for i = 1:n, j = 1:n if i ≥ j)
  rows[1 : nlp.meta.nnzh] .= getindex.(I, 1)
  cols[1 : nlp.meta.nnzh] .= getindex.(I, 2)
  return rows, cols
end

function MPCC.hess_coord!(nlp :: Bard1, x :: AbstractVector, vals :: AbstractVector; obj_weight :: Real = one(eltype(x)))
  #increment!(nlp, :neval_hess)
  ℓ(x) = obj_weight * nlp.f(x)
  Hx = zeros(5,5)
  k = 1
  for j = 1 : nlp.meta.nvar
    for i = j : nlp.meta.nvar
      vals[k] = Hx[i, j]
      k += 1
    end
  end
  return vals
end

function MPCC.hess_coord!(nlp :: Bard1, x :: AbstractVector, y :: AbstractVector, vals :: AbstractVector; obj_weight :: Real = one(eltype(x)))
  #increment!(nlp, :neval_hess)
  ncon, ncc = nlp.meta.ncon, nlp.meta.ncc
  # ℓ(x) = obj_weight * f(x) + λ₁ c(x) - λ₂ G(x) - λ₃ H(x)
  #ℓ(x) = obj_weight * nlp.f(x) + dot(nlp.c(x), y[1:ncon]) - dot(nlp.G(x), y[ncon+1:ncon+ncc]) - dot(nlp.H(x), y[ncon+ncc+1:ncon+2*ncc])
  Hx = zeros(5,5)
  k = 1
  for j = 1 : nlp.meta.nvar
    for i = j : nlp.meta.nvar
      vals[k] = Hx[i, j]
      k += 1
    end
  end
  return vals
end

function MPCC.hessG(nlp :: Bard1, x :: AbstractVector, y :: AbstractVector)
  #increment!(nlp, :neval_hessG)
  #ℓ(x) = dot(nlp.G(x), y)
  Hx = zeros(5,5)
  return tril(Hx)
end

function MPCC.hessG_structure!(nlp :: Bard1, rows :: AbstractVector{<: Integer}, cols :: AbstractVector{<: Integer})
  n = nlp.meta.nvar
  I = ((i,j) for i = 1:n, j = 1:n if i ≥ j)
  rows[1 : nlp.meta.nnzh] .= getindex.(I, 1)
  cols[1 : nlp.meta.nnzh] .= getindex.(I, 2)
  return rows, cols
end

function MPCC.hessG_coord!(nlp :: Bard1, x :: AbstractVector, y :: AbstractVector, vals :: AbstractVector)
  #increment!(nlp, :neval_hessG)
  #ℓ(x) = dot(nlp.G(x), y)
  Hx = zeros(5,5)
  k = 1
  for j = 1 : nlp.meta.nvar
    for i = j : nlp.meta.nvar
      vals[k] = Hx[i, j]
      k += 1
    end
  end
  return vals
end

function MPCC.hessH(nlp :: Bard1, x :: AbstractVector, y :: AbstractVector)
  #increment!(nlp, :neval_hessH)
  #ℓ(x) = dot(nlp.H(x), y)
  Hx = zeros(5,5)
  return tril(Hx)
end

function MPCC.hessH_structure!(nlp :: Bard1, rows :: AbstractVector{<: Integer}, cols :: AbstractVector{<: Integer})
  n = nlp.meta.nvar
  I = ((i,j) for i = 1:n, j = 1:n if i ≥ j)
  rows[1 : nlp.meta.nnzh] .= getindex.(I, 1)
  cols[1 : nlp.meta.nnzh] .= getindex.(I, 2)
  return rows, cols
end

function MPCC.hessH_coord!(nlp :: Bard1, x :: AbstractVector, y :: AbstractVector, vals :: AbstractVector)
  #increment!(nlp, :neval_hess)
  #ℓ(x) = dot(nlp.H(x), y)
  Hx = zeros(5,5)
  k = 1
  for j = 1 : nlp.meta.nvar
    for i = j : nlp.meta.nvar
      vals[k] = Hx[i, j]
      k += 1
    end
  end
  return vals
end

function MPCC.hprod!(nlp :: Bard1, x :: AbstractVector, v :: AbstractVector, Hv :: AbstractVector; obj_weight :: Real = one(eltype(x)))
  #increment!(nlp, :neval_hprod)
  #ℓ(x) = obj_weight * nlp.f(x)
  Hv .= zeros(5)
  return Hv
end

function MPCC.hprod!(nlp :: Bard1, x :: AbstractVector, y :: AbstractVector, v :: AbstractVector, Hv :: AbstractVector; obj_weight :: Real = one(eltype(x)))
  #increment!(nlp, :neval_hprod)
  ncon, ncc = nlp.meta.ncon, nlp.meta.ncc
  #ℓ(x) = obj_weight * nlp.f(x) + dot(nlp.c(x), y[1:ncon]) - dot(nlp.G(x), y[ncon+1:ncon+ncc]) - dot(nlp.H(x), y[ncon+ncc+1:ncon+2*ncc])
  Hv .= zeros(5)
  return Hv
end

function MPCC.hGprod!(nlp :: Bard1, x :: AbstractVector, y :: AbstractVector, v :: AbstractVector, Hv :: AbstractVector)
  #increment!(nlp, :neval_hGprod)
  #ℓ(x) = dot(nlp.G(x), y)
  Hv .= zeros(5)
  return Hv
end

function MPCC.hHprod!(nlp :: Bard1, x :: AbstractVector, y :: AbstractVector, v :: AbstractVector, Hv :: AbstractVector)
  #increment!(nlp, :neval_hHprod)
  #ℓ(x) = dot(nlp.H(x), y)
  Hv .= zeros(5) #ForwardDiff.derivative(t -> ForwardDiff.gradient(ℓ, x + t * v), 0)
  return Hv
end


text:
# This example uses the TMAR Filter from
#
#    @article{doi:10.1175/MWR-D-16-0220.1,
#      author = {Light, Devin and Durran, Dale},
#      title = {Preserving Nonnegativity in Discontinuous Galerkin
#               Approximations to Scalar Transport via Truncation and Mass
#               Aware Rescaling (TMAR)},
#      journal = {Monthly Weather Review},
#      volume = {144},
#      number = {12},
#      pages = {4771-4786},
#      year = {2016},
#      doi = {10.1175/MWR-D-16-0220.1},
#    }
#
# to reproduce the example in section 4b.  It is a shear swirling
# flow deformation of a transported quantity from LeVeque 1996.  The exact
# solution at the final time is the same as the initial condition.
#
using MPI
using CLIMA
using CLIMA.Mesh.Topologies
using CLIMA.Mesh.Grids
using CLIMA.Mesh.Filters
using CLIMA.DGBalanceLawDiscretizations
using CLIMA.MPIStateArrays
using CLIMA.ODESolvers
using CLIMA.GenericCallbacks
using CLIMA.VTK
using LinearAlgebra
using Logging
using Dates
using Printf
using StaticArrays

const finaltime = 5

function velocity(x, y, t)
  sx, cx = sinpi(x), cospi(x)
  sy, cy = sinpi(y), cospi(y)
  ct = cospi(t/finaltime)

  u =  2sx^2*sy*cy*ct
  v = -2sy^2*sx*cx*ct
  (u, v)
end

cosbell(τ, q) = τ ≤ 1 ? ((1 + cospi(τ))/2)^q : zero(τ)

function initialcondition!(Q, x, y, z, _)
  FT = eltype(Q)
  x0, y0 = FT(1//4), FT(1//4)
  τ = 4hypot(x-x0, y-y0)
  @inbounds Q[1] = cosbell(τ, 3)
end

const num_aux_states = 4
const _a_x, _a_y, _a_u, _a_v = 1:num_aux_states

function aux_init!(aux, x, y, z)
  u, v = velocity(x, y, zero(x))

  @inbounds aux[_a_x], aux[_a_y], aux[_a_u], aux[_a_v] = x, y, u, v
end

function advectionflux!(F, state, _, aux, _)
  @inbounds begin
    U, q = @SVector([aux[_a_u], aux[_a_v], 0]), state[1]
    F[:, 1] .= U*q
  end
end

function upwindflux!(fs, nM, stateM, _, auxM, stateP, viscP, auxP, t)
  @inbounds begin
    UM, qM, qP = @SVector([auxM[_a_u], auxM[_a_v], 0]), stateM[1], stateP[1]
    un = dot(nM, UM)

    fs[1] = un ≥ 0 ? un * qM : un * qP
  end
end

function upwindboundaryflux!(fs, nM, stateM, viscM, auxM, stateP, viscP, auxP,
                             bctype, t)
  stateP .= 0
  upwindflux!(fs, nM, stateM, viscM, auxM, stateP, viscP, auxP, t)
end

function preodefun!(disc, Q, t)
  # Update the velocity
  DGBalanceLawDiscretizations.dof_iteration!(disc.auxstate, disc, Q) do R, _, _, aux
      @inbounds R[_a_u], R[_a_v] = velocity(aux[_a_x], aux[_a_y], t)
  end

  Filters.apply!(Q, 1, disc.grid, TMARFilter())
end

function setupDG(mpicomm, dim, Ne, polynomialorder, FT=Float64, ArrayType=Array)
  brickrange = (range(FT(0); length=Ne+1, stop=1),
                range(FT(0); length=Ne+1, stop=1),
                range(FT(0); length=Ne+1, stop=1))

  topology = BrickTopology(mpicomm, brickrange[1:dim])

  grid = DiscontinuousSpectralElementGrid(topology; polynomialorder =
                                          polynomialorder, FloatType = FT,
                                          DeviceArray = ArrayType,)

  spatialdiscretization = DGBalanceLaw(grid = grid, length_state_vector = 1,
                                       flux! = advectionflux!,
                                       numerical_flux! = upwindflux!,
                                       numerical_boundary_flux! =
                                       upwindboundaryflux!,
                                       auxiliary_state_length = num_aux_states,
                                       auxiliary_state_initialization! = aux_init!,
                                       preodefun! = preodefun!)
end

function run()
  CLIMA.init()
  DeviceArrayType = CLIMA.array_type()

  mpicomm = MPI.COMM_WORLD
  mpi_logger = ConsoleLogger(MPI.Comm_rank(mpicomm) == 0 ? stderr : devnull)
  rank = MPI.Comm_rank(mpicomm)

  dim = 2
  Ne = 20
  polynomialorder = 4
  FT = Float64

  spatialdiscretization = setupDG(mpicomm, dim, Ne, polynomialorder, FT,
                                  DeviceArrayType)
  Q = MPIStateArray(spatialdiscretization, initialcondition!)

  maxvelosity = 2
  elementsize = 1 / Ne
  dx = elementsize / polynomialorder^2
  CFL = 1

  dt = CFL * dx / maxvelosity
  @info @sprintf "dt = %1.2e" dt
  sork = SSPRK33ShuOsher(spatialdiscretization, Q; dt = dt, t0 = 0)

  initialsumQ = weightedsum(Q)

  vtk_step = 0
  mkpath("vtk")
  function vtkoutput()
    Filters.apply!(Q, 1, spatialdiscretization.grid, TMARFilter())

    filename = @sprintf("vtk/q_rank%04d_step%04d", rank, vtk_step)
    writevtk(filename, Q, spatialdiscretization, ("q",))

    minQ = MPI.Reduce([minimum(Q.realdata)], MPI.MIN, 0, Q.mpicomm)
    maxQ = MPI.Reduce([maximum(Q.realdata)], MPI.MAX, 0, Q.mpicomm)
    sumQ = weightedsum(Q)

    with_logger(mpi_logger) do
      @info @sprintf("""step = %d
                           min Q = %25.16e
                           max Q = %25.16e
                     sum error Q = %25.16e
                     """, vtk_step, minQ[1], maxQ[1], (initialsumQ -
                                                       sumQ)/initialsumQ)
    end

    vtk_step += 1
    nothing
  end

  cb_vtk = GenericCallbacks.EveryXSimulationSteps(vtkoutput, 40)

  vtkoutput()

  # We integrate so that the final solution is equal to the initial solution
  Qe = copy(Q)

  solve!(Q, sork; timeend = finaltime, callbacks = (cb_vtk, ))

  vtkoutput()

  minQ = MPI.Reduce([minimum(Q.realdata)], MPI.MIN, 0, Q.mpicomm)
  maxQ = MPI.Reduce([maximum(Q.realdata)], MPI.MAX, 0, Q.mpicomm)
  finalsumQ = weightedsum(Q)
  sumerror = (initialsumQ - finalsumQ) / initialsumQ
  error = euclidean_distance(Q, Qe)
  with_logger(mpi_logger) do
    @info @sprintf("""Run with
                   dim              = %d
                   Ne               = %d
                   polynomial order = %d
                   min              = %e
                   max              = %e
                   L2 error         = %e
                   sum error        = %e
                   """, dim, Ne, polynomialorder, minQ[1], maxQ[1], error,
                   sumerror)
  end
end

run()

max_stars_repo_name: favba/CLIMA

max_issues_repo_name: favba/CLIMA

max_forks_repo_name: favba/CLIMA


text:
# Warp Vote (B.13)

for mode in (:all, :any, :uni)
    fname = Symbol("vote_$mode")
    fname_sync = Symbol("vote_$(mode)_sync")
    @eval export $fname, $fname_sync

    intrinsic = "llvm.nvvm.vote.$mode"
    @eval begin
        # FIXME: ccall($intrinsic, llvmcall, $rettyp, (Bool,), pred)
        #        doesn't use i1 for Bool
        @inline $fname(pred) =
            Base.llvmcall($("""
                declare i1 @$intrinsic(i1)

                define i8 @entry(i8) #0 {
                    %predicate = icmp eq i8 %0, 1
                    %llvmbool = call i1 @$intrinsic(i1 %predicate)
                    %jlbool = zext i1 %llvmbool to i8
                    ret i8 %jlbool
                }

                attributes #0 = { alwaysinline }""", "entry"),
            Bool, Tuple{Bool}, pred)
    end

    # warp-synchronous
    intrinsic = "llvm.nvvm.vote.$mode.sync"
    @eval begin
        @inline $fname_sync(mask, pred) =
            Base.llvmcall($("""
                declare i1 @$intrinsic(i32, i1)

                define i8 @entry(i32 %mask, i8) #0 {
                    %predicate = icmp eq i8 %0, 1
                    %llvmbool = call i1 @$intrinsic(i32 %mask, i1 %predicate)
                    %jlbool = zext i1 %llvmbool to i8
                    ret i8 %jlbool
                }

                attributes #0 = { alwaysinline }""", "entry"),
            Bool, Tuple{UInt32, Bool}, mask, pred)
    end
end

# ballot returns an integer, so we need to repeat the above
for mode in (:ballot, )
    fname = Symbol("vote_$mode")
    fname_sync = Symbol("vote_$(mode)_sync")
    @eval export $fname, $fname_sync

    intrinsic = "llvm.nvvm.vote.$mode"
    @eval begin
        @inline $fname(pred) =
            Base.llvmcall($("""
                declare i32 @$intrinsic(i1)

                define i32 @entry(i8) #0 {
                    %predicate = icmp eq i8 %0, 1
                    %ret = call i32 @$intrinsic(i1 %predicate)
                    ret i32 %ret
                }

                attributes #0 = { alwaysinline }""", "entry"),
            UInt32, Tuple{Bool}, pred)
    end

    # warp-synchronous
    intrinsic = "llvm.nvvm.vote.$mode.sync"
    @eval begin
        @inline $fname_sync(mask, pred) =
            Base.llvmcall($("""
                declare i32 @$intrinsic(i32, i1)

                define i32 @entry(i32 %mask, i8) #0 {
                    %predicate = icmp eq i8 %0, 1
                    %ret = call i32 @$intrinsic(i32 %mask, i1 %predicate)
                    ret i32 %ret
                }

                attributes #0 = { alwaysinline }""", "entry"),
            UInt32, Tuple{UInt32, Bool}, mask, pred)
    end
end

"""
    vote_all(predicate::Bool)
    vote_all_sync(mask::UInt32, predicate::Bool)

Evaluate `predicate` for all active threads of the warp and return whether `predicate`
is true for all of them.
"""
vote_all
@doc (@doc vote_all) vote_all_sync

"""
    vote_any(predicate::Bool)
    vote_any_sync(mask::UInt32, predicate::Bool)

Evaluate `predicate` for all active threads of the warp and return whether `predicate`
is true for any of them.
"""
vote_any
@doc (@doc vote_any) vote_any_sync

"""
    vote_uni(predicate::Bool)
    vote_uni_sync(mask::UInt32, predicate::Bool)

Evaluate `predicate` for all active threads of the warp and return whether `predicate` is
the same for any of them.
"""
vote_uni
@doc (@doc vote_uni) vote_uni_sync

"""
    vote_ballot(predicate::Bool)
    vote_ballot_sync(mask::UInt32, predicate::Bool)

Evaluate `predicate` for all active threads of the warp and return an integer whose Nth bit
is set if and only if `predicate` is true for the Nth thread of the warp and the Nth thread
is active.
"""
vote_ballot
@doc (@doc vote_ballot) vote_ballot_sync


text:
#=
list of predefined signomials and domains from various applications

prefixes:
- CS16 refers to
"Relative Entropy Relaxations for Signomial Optimization"
(2016) by Chandrasekaran & Shah
- MCW19 refers to
"Signomial and Polynomial Optimization via Relative Entropy and Partial Dualization"
(2019) by Murray, Chandrasekaran, & Wierman

c are coefficients
A are powers
x is a feasible point in domain, if known
obj_ub is an objective upper bound (value at feasible point), if known
=#

using LinearAlgebra
import Distributions
import Random
using SparseArrays

signomialmin_data = Dict{Symbol, NamedTuple}(
    :motzkin2 => (
        # f = 1 - 3*x1^2*x2^2 + x1^2*x2^4 + x1^4*x2^2
        fc = [1, -3, 1, 1],
        fA = [0 0; 2 2; 2 4; 4 2],
        gc = [],
        gA = [],
        obj_ub = 0.0,
        ),
    :motzkin3 => (
        # f = x3^6 - 3*x1^2*x2^2*x3^2 + x1^2*x2^4 + x1^4*x2^2
        fc = [0, 1, -3, 1, 1],
        fA = [0 0 0; 0 0 6; 2 2 2; 2 4 0; 4 2 0],
        gc = [],
        gA = [],
        obj_ub = 0.0,
        ),
    :CS16ex8_13 => (
        fc = [0, 10, 10, 10, -14.6794, -7.8601, 8.7838],
        fA = [0 0 0; 10.2 0 0; 0 9.8 0; 0 0 8.2; 1.5089 1.0981 1.3419;
            1.0857 1.9069 1.6192; 1.0459 0.0492 1.6245],
        gc = [],
        gA = [],
        obj_ub = -0.9748,
        ),
    :CS16ex8_14 => (
        fc = [0, 10, 10, 10, 7.5907, -10.9888, -13.9164],
        fA = [0 0 0; 10.2 0 0; 0 9.8 0; 0 0 8.2; 1.9864 0.2010 1.0855;
            2.8242 1.9355 2.0503; 0.1828 2.7772 1.9001],
        gc = [],
        gA = [],
        obj_ub = -0.739,
        ),
    :CS16ex18 => (
        fc = [0, 10, 10, 10, -14.6794, -7.8601, 8.7838],
        fA = [0 0 0; 10.2070 0.0082 -0.0039; -0.0081 9.8024 -0.0097;
            0.0070 -0.0156 8.1923; 1.5296 1.0927 1.3441;
            1.0750 1.9108 1.6339; 1.0513 0.0571 1.6188],
        gc = [],
        gA = [],
        obj_ub = -0.9441,
        ),
    :CS16ex12 => (
        fc = [0, 10, 10, 10, -14.6794, -7.8601, 8.7838],
        fA = [0 0 0; 10.2 0 0; 0 9.8 0; 0 0 8.2; 1.5089 1.0981 1.3419;
            1.0857 1.9069 1.6192; 1.0459 0.0492 1.6245],
        gc = [[1, -8, -8, -8, -6.4]],
        gA = [[0 0 0; 10.2 0 0; 0 9.8 0; 0 0 8.2; 1.0857 1.9069 1.6192]],
        obj_ub = -0.6144,
        ),
    :CS16ex13 => (
        fc = [0, 10, 10, 10, -14.6794, -7.8601, 8.7838],
        fA = [0 0 0; 10.2 0 0; 0 9.8 0; 0 0 8.2; 1.5089 1.0981 1.3419;
            1.0857 1.9069 1.6192; 1.0459 0.0492 1.6245],
        gc = [[-8, -8, -8, 0.7410, -0.4492, 1.4240]],
        gA = [[10.2 0 0; 0 9.8 0; 0 0 8.2; 1.5089 1.0981 1.3419;
            1.0857 1.9069 1.6192; 1.0459 0.0492 1.6245]],
        obj_ub = -0.7372,
        ),
    :MCW19ex1_mod => (
        fc = [0, 0.5, -1, 0],
        fA = [0 0 0; 1 -1 0; 1 0 0; 0 -1 0],
        gc = [[100, -1, -1, -0.05], [-70, 1], [-1, 1], [-0.5, 1], [150, -1],
            [30, -1], [21, -1]],
        gA = [[0 0 0; 0 1 -1; 0 1 0; 1 0 1], [0 0 0; 1 0 0], [0 0 0; 0 1 0],
            [0 0 0; 0 0 1], [0 0 0; 1 0 0], [0 0 0; 0 1 0], [0 0 0; 0 0 1]],
        obj_ub = -147.5,
        ),
    :MCW19ex8 => (
        fc = [0, 0.05, 0.05, 0.05, 1],
        fA = sparse([2, 3, 4, 5], [1, 2, 3, 9], [1, 1, 1, 1], 5, 10),
        gc = [[1, 0.5, -1], [1, 0.5, -1], [1, 0.5, -1], [1, -0.25, -0.5],
            [1, -0.79681], [1, -0.79681], [1, -0.79681]],
        gA = [
            sparse([2, 2, 2, 3, 3], [1, 4, 7, 10, 7], [1, 1, -1, 1, -1], 3, 10),
            sparse([2, 2, 2, 3, 3], [2, 5, 8, 7, 8], [1, 1, -1, 1, -1], 3, 10),
            sparse([2, 2, 2, 3, 3], [3, 6, 9, 8, 9], [1, 1, -1, 1, -1], 3, 10),
            sparse([2, 3, 3], [10, 9, 10], [-1, 1, -1], 3, 10),
            sparse([2, 2], [4, 7], [1, -1], 2, 10),
            sparse([2, 2], [5, 8], [1, -1], 2, 10),
            sparse([2, 2], [6, 9], [1, -1], 2, 10),
            ],
        obj_ub = 0.2056534,
        ),
    )

function signomialmin_random(
    m::Int,
    n::Int;
    num_samples::Int = 100,
    neg_c_frac::Real = 0.0,
    sparsity::Real = 0.3,
    rseed::Int = 1,
    )
    Random.seed!(rseed)

    # generate random objective signomial
    fc = rand(m)
    for k in 1:m
        if rand() < neg_c_frac
            fc[k] *= -1
        end
    end
    fA = vcat(zeros(1, n), Matrix(sprandn(m - 1, n, sparsity)))
    for k in 2:size(fA, 1)
        if iszero(norm(fA[k, :]))
            fA[k, :] = randn(n)
        end
    end

    # bounded domain set (in exp space) is intersection of positives and ball
    gc = [vcat(1, fill(-1, n))]
    gA = [sparse(2:(n + 1), 1:n, fill(2, n))]

    # sample points to get an objective upper bound
    eval_signomial(c::Vector, A::AbstractMatrix, x::Vector) =
        sum(c_k * exp(dot(A_k, x)) for (c_k, A_k) in zip(c, eachrow(A)))
    obj_ub = Inf
    for i in 1:num_samples
        x = abs.(randn(n))
        r = rand(Distributions.Exponential(0.5))
        x /= sqrt(sum(abs2, x) + r)
        x = log.(x)
        @assert eval_signomial(gc[1], gA[1], x) >= 0
        obj_ub = min(obj_ub, eval_signomial(fc, fA, x))
    end

    return (fc, fA, gc, gA, obj_ub)
end


text:
"""
$(SIGNATURES)

Python linspace function
"""
linspace(start, stop, n::T) where {T<:Integer} =
    collect(range(start, stop = stop, length = n))

"""
$(SIGNATURES)

Heaviside step function
"""
heaviside(x::T) where {T<:Real} = ifelse(x >= 0, one(x), zero(x))

"""
$(SIGNATURES)

Fortran sign function
"""
fortsign(x, y) = abs(x) * sign(y)

"""
$(SIGNATURES)

Split matrix into row vectors
"""
function mat_split(m::AM{T}) where {T<:Number}

    if length(m[1, :]) == 2
        nx = eltype(m).([1.0 0.0])
        ny = eltype(m).([0.0 1.0])

        return nx * m, ny * m
    elseif length(m[1, :]) == 3
        nx = eltype(m).([1.0 0.0 0.0])
        ny = eltype(m).([0.0 1.0 0.0])
        nz = eltype(m).([0.0 0.0 1.0])

        return nx * m, ny * m, nz * m
    end

end

"""
$(SIGNATURES)

Gauss Legendre integral for fast spectral method

* @args: number of quadrature points N, integral range [a, b]
* @args: quadrature points x & weights w

"""
function lgwt(N::Integer, a::Real, b::Real)
    x = zeros(N)
    w = zeros(N)

    N1 = N
    N2 = N + 1

    y = zeros(N1)
    y0 = zeros(N1)
    Lp = zeros(N1)
    L = zeros(N1, N2)

    # initial guess
    for i = 1:N1
        y[i] =
            cos((2.0 * (i - 1.0) + 1.0) * 4.0 * atan(1.0) / (2.0 * (N - 1.0) + 2.0)) +
            0.27 / N1 *
            sin(4.0 * atan(1.0) * (-1.0 + i * 2.0 / (N1 - 1.0)) * (N - 1.0) / N2)
        y0[i] = 2.0
    end

    # compute the zeros of the N+1 legendre Polynomial
    # using the recursion relation and the Newton method
    while maximum(abs.(y .- y0)) > 0.0000000000001
        L[:, 1] .= 1.0
        L[:, 2] .= y
        for k = 2:N1
            @. L[:, k+1] = ((2.0 * k - 1.0) * y * L[:, k] - (k - 1) * L[:, k-1]) / k
        end
        @. Lp = N2 * (L[:, N1] - y * L[:, N2]) / (1.0 - y^2)
        @. y0 = y
        @. y = y0 - L[:, N2] / Lp
    end

    # linear map from [-1 1] to [a,b]
    @. x = (a * (1.0 - y) + b * (1.0 + y)) / 2.0
    @. w = N2^2 * (b - a) / ((1.0 - y^2) * Lp^2) / N1^2

    return x, w
end

"""
$(SIGNATURES)

Extract subarray except the last column
"""
function extract_last(a::AA{T}, idx::Integer; mode = :view::Symbol) where {T}
    if mode == :copy

        if ndims(a) == 2
            sw = a[:, idx]
        elseif ndims(a) == 3
            sw = a[:, :, idx]
        elseif ndims(a) == 4
            sw = a[:, :, :, idx]
        elseif ndims(a) == 5
            sw = a[:, :, :, :, idx]
        end

    elseif mode == :view

        if ndims(a) == 2
            sw = @view a[:, idx]
        elseif ndims(a) == 3
            sw = @view a[:, :, idx]
        elseif ndims(a) == 4
            sw = @view a[:, :, :, idx]
        elseif ndims(a) == 5
            sw = @view a[:, :, :, :, idx]
        end

    else

        throw("Error in extraction mode setup")

    end

    return sw
end

"""
$(SIGNATURES)

Calculate rate of convergence
"""
convergence_order(e1, e2, r = 2) = log(e1 / e2) / log(r)

"""
$(SIGNATURES)

Calculate L1 error
"""
L1_error(u::T, ue::T, Δx) where {T<:AA} = sum(abs.(u .- ue) .* Δx)

"""
$(SIGNATURES)

Calculate L2 error
"""
L2_error(u::T, ue::T, Δx) where {T<:AA} = sqrt(sum((abs.(u .- ue) .* Δx) .^ 2))

"""
$(SIGNATURES)

Calculate L∞ error
"""
L∞_error(u::T, ue::T, Δx) where {T<:AA} = maximum(abs.(u .- ue) .* Δx)


text:
# Automatically generated using Clang.jl

const IVISCOPE_MAJOR_VERSION = Int32(4)
const IVISCOPE_MINOR_VERSION = Int32(0)
const IVISCOPE_SPECIFICATION_MAJOR_VERSION = Int32(3)
const IVISCOPE_SPECIFICATION_MINOR_VERSION = Int32(0)
const IVISCOPE_DRIVER_VENDOR = "National Instruments"
const IVISCOPE_DRIVER_DESCRIPTION = "IviScope Class Driver"
const IVISCOPE_ATTR_CACHE = IVI_ATTR_CACHE
const IVISCOPE_ATTR_RANGE_CHECK = IVI_ATTR_RANGE_CHECK
const IVISCOPE_ATTR_QUERY_INSTRUMENT_STATUS = IVI_ATTR_QUERY_INSTRUMENT_STATUS
const IVISCOPE_ATTR_RECORD_COERCIONS = IVI_ATTR_RECORD_COERCIONS
const IVISCOPE_ATTR_SIMULATE = IVI_ATTR_SIMULATE
const IVISCOPE_ATTR_INTERCHANGE_CHECK = IVI_ATTR_INTERCHANGE_CHECK
const IVISCOPE_ATTR_SPY = IVI_ATTR_SPY
const IVISCOPE_ATTR_USE_SPECIFIC_SIMULATION = IVI_ATTR_USE_SPECIFIC_SIMULATION
const IVISCOPE_ATTR_GROUP_CAPABILITIES = IVI_ATTR_GROUP_CAPABILITIES
const IVISCOPE_ATTR_FUNCTION_CAPABILITIES = IVI_ATTR_FUNCTION_CAPABILITIES
const IVISCOPE_ATTR_CLASS_DRIVER_PREFIX = IVI_ATTR_CLASS_DRIVER_PREFIX
const IVISCOPE_ATTR_CLASS_DRIVER_VENDOR = IVI_ATTR_CLASS_DRIVER_VENDOR
const IVISCOPE_ATTR_CLASS_DRIVER_DESCRIPTION = IVI_ATTR_CLASS_DRIVER_DESCRIPTION
const IVISCOPE_ATTR_CLASS_DRIVER_CLASS_SPEC_MAJOR_VERSION = IVI_ATTR_CLASS_DRIVER_CLASS_SPEC_MAJOR_VERSION
const IVISCOPE_ATTR_CLASS_DRIVER_CLASS_SPEC_MINOR_VERSION = IVI_ATTR_CLASS_DRIVER_CLASS_SPEC_MINOR_VERSION
const IVISCOPE_ATTR_SPECIFIC_DRIVER_PREFIX = IVI_ATTR_SPECIFIC_DRIVER_PREFIX
const IVISCOPE_ATTR_SPECIFIC_DRIVER_LOCATOR = IVI_ATTR_SPECIFIC_DRIVER_LOCATOR
const IVISCOPE_ATTR_IO_RESOURCE_DESCRIPTOR = IVI_ATTR_IO_RESOURCE_DESCRIPTOR
const IVISCOPE_ATTR_LOGICAL_NAME = IVI_ATTR_LOGICAL_NAME
const IVISCOPE_ATTR_SPECIFIC_DRIVER_VENDOR = IVI_ATTR_SPECIFIC_DRIVER_VENDOR
const IVISCOPE_ATTR_SPECIFIC_DRIVER_DESCRIPTION = IVI_ATTR_SPECIFIC_DRIVER_DESCRIPTION
const IVISCOPE_ATTR_SPECIFIC_DRIVER_CLASS_SPEC_MAJOR_VERSION = IVI_ATTR_SPECIFIC_DRIVER_CLASS_SPEC_MAJOR_VERSION
const IVISCOPE_ATTR_SPECIFIC_DRIVER_CLASS_SPEC_MINOR_VERSION = IVI_ATTR_SPECIFIC_DRIVER_CLASS_SPEC_MINOR_VERSION
const IVISCOPE_ATTR_INSTRUMENT_FIRMWARE_REVISION = IVI_ATTR_INSTRUMENT_FIRMWARE_REVISION
const IVISCOPE_ATTR_INSTRUMENT_MANUFACTURER = IVI_ATTR_INSTRUMENT_MANUFACTURER
const IVISCOPE_ATTR_INSTRUMENT_MODEL = IVI_ATTR_INSTRUMENT_MODEL
const IVISCOPE_ATTR_SUPPORTED_INSTRUMENT_MODELS = IVI_ATTR_SUPPORTED_INSTRUMENT_MODELS
const IVISCOPE_ATTR_CLASS_DRIVER_REVISION = IVI_ATTR_CLASS_DRIVER_REVISION
const IVISCOPE_ATTR_SPECIFIC_DRIVER_REVISION = IVI_ATTR_SPECIFIC_DRIVER_REVISION
const IVISCOPE_ATTR_DRIVER_SETUP = IVI_ATTR_DRIVER_SETUP
const IVISCOPE_ATTR_CHANNEL_COUNT = IVI_ATTR_CHANNEL_COUNT
const IVISCOPE_ATTR_VERTICAL_RANGE = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(1)
const IVISCOPE_ATTR_VERTICAL_OFFSET = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(2)
const IVISCOPE_ATTR_VERTICAL_COUPLING = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(3)
const IVISCOPE_ATTR_PROBE_ATTENUATION = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(4)
const IVISCOPE_ATTR_CHANNEL_ENABLED = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(5)
const IVISCOPE_ATTR_MAX_INPUT_FREQUENCY = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(6)
const IVISCOPE_ATTR_INPUT_IMPEDANCE = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(103)
const IVISCOPE_ATTR_ACQUISITION_TYPE = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(101)
const IVISCOPE_ATTR_ACQUISITION_START_TIME = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(109)
const IVISCOPE_ATTR_HORZ_TIME_PER_RECORD = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(7)
const IVISCOPE_ATTR_HORZ_RECORD_LENGTH = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(8)
const IVISCOPE_ATTR_HORZ_MIN_NUM_PTS = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(9)
const IVISCOPE_ATTR_HORZ_SAMPLE_RATE = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(10)
const IVISCOPE_ATTR_TRIGGER_TYPE = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(12)
const IVISCOPE_ATTR_TRIGGER_SOURCE = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(13)
const IVISCOPE_ATTR_TRIGGER_COUPLING = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(14)
const IVISCOPE_ATTR_TRIGGER_HOLDOFF = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(16)
const IVISCOPE_ATTR_TRIGGER_LEVEL = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(17)
const IVISCOPE_ATTR_TRIGGER_SLOPE = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(18)
const IVISCOPE_ATTR_TV_TRIGGER_SIGNAL_FORMAT = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(201)
const IVISCOPE_ATTR_TV_TRIGGER_EVENT = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(205)
const IVISCOPE_ATTR_TV_TRIGGER_LINE_NUMBER = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(206)
const IVISCOPE_ATTR_TV_TRIGGER_POLARITY = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(204)
const IVISCOPE_ATTR_RUNT_HIGH_THRESHOLD = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(301)
const IVISCOPE_ATTR_RUNT_LOW_THRESHOLD = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(302)
const IVISCOPE_ATTR_RUNT_POLARITY = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(303)
const IVISCOPE_ATTR_GLITCH_WIDTH = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(401)
const IVISCOPE_ATTR_GLITCH_POLARITY = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(402)
const IVISCOPE_ATTR_GLITCH_CONDITION = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(403)
const IVISCOPE_ATTR_WIDTH_LOW_THRESHOLD = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(501)
const IVISCOPE_ATTR_WIDTH_HIGH_THRESHOLD = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(502)
const IVISCOPE_ATTR_WIDTH_POLARITY = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(503)
const IVISCOPE_ATTR_WIDTH_CONDITION = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(504)
const IVISCOPE_ATTR_AC_LINE_TRIGGER_SLOPE = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(701)
const IVISCOPE_ATTR_NUM_ENVELOPES = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(105)
const IVISCOPE_ATTR_MEAS_HIGH_REF = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(607)
const IVISCOPE_ATTR_MEAS_LOW_REF = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(608)
const IVISCOPE_ATTR_MEAS_MID_REF = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(609)
const IVISCOPE_ATTR_TRIGGER_MODIFIER = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(102)
const IVISCOPE_ATTR_NUM_AVERAGES = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(104)
const IVISCOPE_ATTR_SAMPLE_MODE = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(106)
const IVISCOPE_ATTR_INITIATE_CONTINUOUS = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(107)
const IVISCOPE_ATTR_PROBE_SENSE_VALUE = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(108)
const IVISCOPE_ATTR_INTERPOLATION = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(19)
const IVISCOPE_VAL_MAX_TIME_INFINITE = IVI_VAL_MAX_TIME_INFINITE
const IVISCOPE_VAL_MAX_TIME_IMMEDIATE = IVI_VAL_MAX_TIME_IMMEDIATE
const IVISCOPE_VAL_ACQ_COMPLETE = Int32(1)
const IVISCOPE_VAL_ACQ_IN_PROGRESS = Int32(0)
const IVISCOPE_VAL_ACQ_STATUS_UNKNOWN = -(Int32(1))
const IVISCOPE_VAL_RISE_TIME = Int32(0)
const IVISCOPE_VAL_FALL_TIME = Int32(1)
const IVISCOPE_VAL_FREQUENCY = Int32(2)
const IVISCOPE_VAL_PERIOD = Int32(3)
const IVISCOPE_VAL_VOLTAGE_RMS = Int32(4)
const IVISCOPE_VAL_VOLTAGE_PEAK_TO_PEAK = Int32(5)
const IVISCOPE_VAL_VOLTAGE_MAX = Int32(6)
const IVISCOPE_VAL_VOLTAGE_MIN = Int32(7)
const IVISCOPE_VAL_VOLTAGE_HIGH = Int32(8)
const IVISCOPE_VAL_VOLTAGE_LOW = Int32(9)
const IVISCOPE_VAL_VOLTAGE_AVERAGE = Int32(10)
const IVISCOPE_VAL_WIDTH_NEG = Int32(11)
const IVISCOPE_VAL_WIDTH_POS = Int32(12)
const IVISCOPE_VAL_DUTY_CYCLE_NEG = Int32(13)
const IVISCOPE_VAL_DUTY_CYCLE_POS = Int32(14)
const IVISCOPE_VAL_AMPLITUDE = Int32(15)
const IVISCOPE_VAL_VOLTAGE_CYCLE_RMS = Int32(16)
const IVISCOPE_VAL_VOLTAGE_CYCLE_AVERAGE = Int32(17)
const IVISCOPE_VAL_OVERSHOOT = Int32(18)
const IVISCOPE_VAL_PRESHOOT = Int32(19)
const IVISCOPE_VAL_MEASUREMENT_FUNCTION_CLASS_EXT_BASE = Int32(100)
const IVISCOPE_VAL_MEASUREMENT_FUNCTION_SPECIFIC_EXT_BASE = Int32(1000)
const IVISCOPE_VAL_AC = Int32(0)
const IVISCOPE_VAL_DC = Int32(1)
const IVISCOPE_VAL_GND = Int32(2)
const IVISCOPE_VAL_VERTICAL_COUPLING_CLASS_EXT_BASE = Int32(100)
const IVISCOPE_VAL_VERTICAL_COUPLING_SPECIFIC_EXT_BASE = Int32(1000)
const IVISCOPE_VAL_EDGE_TRIGGER = Int32(1)
const IVISCOPE_VAL_WIDTH_TRIGGER = Int32(2)
const IVISCOPE_VAL_RUNT_TRIGGER = Int32(3)
const IVISCOPE_VAL_GLITCH_TRIGGER = Int32(4)
const IVISCOPE_VAL_TV_TRIGGER = Int32(5)
const IVISCOPE_VAL_IMMEDIATE_TRIGGER = Int32(6)
const IVISCOPE_VAL_AC_LINE_TRIGGER = Int32(7)
const IVISCOPE_VAL_TRIGGER_TYPE_CLASS_EXT_BASE = Int32(200)
const IVISCOPE_VAL_TRIGGER_TYPE_SPECIFIC_EXT_BASE = Int32(1000)
const IVISCOPE_VAL_POSITIVE = Int32(1)
const IVISCOPE_VAL_NEGATIVE = Int32(0)
const IVISCOPE_VAL_EXTERNAL = "VAL_EXTERNAL"
const IVISCOPE_VAL_TTL0 = "VAL_TTL0"
const IVISCOPE_VAL_TTL1 = "VAL_TTL1"
const IVISCOPE_VAL_TTL2 = "VAL_TTL2"
const IVISCOPE_VAL_TTL3 = "VAL_TTL3"
const IVISCOPE_VAL_TTL4 = "VAL_TTL4"
const IVISCOPE_VAL_TTL5 = "VAL_TTL5"
const IVISCOPE_VAL_TTL6 = "VAL_TTL6"
const IVISCOPE_VAL_TTL7 = "VAL_TTL7"
const IVISCOPE_VAL_ECL0 = "VAL_ECL0"
const IVISCOPE_VAL_ECL1 = "VAL_ECL1"
const IVISCOPE_VAL_PXI_STAR = "VAL_PXI_STAR"
const IVISCOPE_VAL_RTSI_0 = "VAL_RTSI_0"
const IVISCOPE_VAL_RTSI_1 = "VAL_RTSI_1"
const IVISCOPE_VAL_RTSI_2 = "VAL_RTSI_2"
const IVISCOPE_VAL_RTSI_3 = "VAL_RTSI_3"
const IVISCOPE_VAL_RTSI_4 = "VAL_RTSI_4"
const IVISCOPE_VAL_RTSI_5 = "VAL_RTSI_5"
const IVISCOPE_VAL_RTSI_6 = "VAL_RTSI_6"
const IVISCOPE_VAL_PROBE_SENSE_ON = -(Int32(1))
const IVISCOPE_VAL_PROBE_ATTENUATION_CLASS_EXT_BASE = -(Int32(100))
const IVISCOPE_VAL_PROBE_ATTENUATION_SPECIFIC_EXT_BASE = -(Int32(1000))
const IVISCOPE_VAL_HF_REJECT = Int32(3)
const IVISCOPE_VAL_LF_REJECT = Int32(4)
const IVISCOPE_VAL_NOISE_REJECT = Int32(5)
const IVISCOPE_VAL_TRIGGER_COUPLING_CLASS_EXT_BASE = Int32(100)
const IVISCOPE_VAL_TRIGGER_COUPLING_SPECIFIC_EXT_BASE = Int32(1000)
const IVISCOPE_VAL_NO_INTERPOLATION = Int32(1)
const IVISCOPE_VAL_SINE_X = Int32(2)
const IVISCOPE_VAL_LINEAR = Int32(3)
const IVISCOPE_VAL_INTERPOLATION_CLASS_EXT_BASE = Int32(100)
const IVISCOPE_VAL_INTERPOLATION_SPECIFIC_EXT_BASE = Int32(1000)
const IVISCOPE_VAL_NTSC = Int32(1)
const IVISCOPE_VAL_PAL = Int32(2)
const IVISCOPE_VAL_SECAM = Int32(3)
const IVISCOPE_VAL_TV_SIGNAL_FORMAT_CLASS_EXT_BASE = Int32(100)
const IVISCOPE_VAL_TV_SIGNAL_FORMAT_SPECIFIC_EXT_BASE = Int32(1000)
const IVISCOPE_VAL_TV_EVENT_FIELD1 = Int32(1)
const IVISCOPE_VAL_TV_EVENT_FIELD2 = Int32(2)
const IVISCOPE_VAL_TV_EVENT_ANY_FIELD = Int32(3)
const IVISCOPE_VAL_TV_EVENT_ANY_LINE = Int32(4)
const IVISCOPE_VAL_TV_EVENT_LINE_NUMBER = Int32(5)
const IVISCOPE_VAL_TV_TRIGGER_EVENT_CLASS_EXT_BASE = Int32(100)
const IVISCOPE_VAL_TV_TRIGGER_EVENT_SPECIFIC_EXT_BASE = Int32(1000)
const IVISCOPE_VAL_TV_POSITIVE = Int32(1)
const IVISCOPE_VAL_TV_NEGATIVE = Int32(2)
const IVISCOPE_VAL_TV_TRIGGER_POLARITY_CLASS_EXT_BASE = Int32(100)
const IVISCOPE_VAL_TV_TRIGGER_POLARITY_SPECIFIC_EXT_BASE = Int32(1000)
const IVISCOPE_VAL_RUNT_POSITIVE = Int32(1)
const IVISCOPE_VAL_RUNT_NEGATIVE = Int32(2)
const IVISCOPE_VAL_RUNT_EITHER = Int32(3)
const IVISCOPE_VAL_GLITCH_POSITIVE = Int32(1)
const IVISCOPE_VAL_GLITCH_NEGATIVE = Int32(2)
const IVISCOPE_VAL_GLITCH_EITHER = Int32(3)
const IVISCOPE_VAL_GLITCH_LESS_THAN = Int32(1)
const IVISCOPE_VAL_GLITCH_GREATER_THAN = Int32(2)
const IVISCOPE_VAL_WIDTH_POSITIVE = Int32(1)
const IVISCOPE_VAL_WIDTH_NEGATIVE = Int32(2)
const IVISCOPE_VAL_WIDTH_EITHER = Int32(3)
const IVISCOPE_VAL_WIDTH_WITHIN = Int32(1)
const IVISCOPE_VAL_WIDTH_OUTSIDE = Int32(2)
const IVISCOPE_VAL_AC_LINE_POSITIVE = Int32(1)
const IVISCOPE_VAL_AC_LINE_NEGATIVE = Int32(2)
const IVISCOPE_VAL_AC_LINE_EITHER = Int32(3)
const IVISCOPE_VAL_NORMAL = Int32(0)
const IVISCOPE_VAL_PEAK_DETECT = Int32(1)
const IVISCOPE_VAL_HI_RES = Int32(2)
const IVISCOPE_VAL_ENVELOPE = Int32(3)
const IVISCOPE_VAL_AVERAGE = Int32(4)
const IVISCOPE_VAL_ACQUISITION_TYPE_CLASS_EXT_BASE = Int32(100)
const IVISCOPE_VAL_ACQUISITION_TYPE_SPECIFIC_EXT_BASE = Int32(1000)
const IVISCOPE_VAL_NO_TRIGGER_MOD = Int32(1)
const IVISCOPE_VAL_AUTO = Int32(2)
const IVISCOPE_VAL_AUTO_LEVEL = Int32(3)
const IVISCOPE_VAL_TRIGGER_MOD_CLASS_EXT_BASE = Int32(100)
const IVISCOPE_VAL_TRIGGER_MOD_SPECIFIC_EXT_BASE = Int32(1000)
const IVISCOPE_VAL_REAL_TIME = Int32(0)
const IVISCOPE_VAL_EQUIVALENT_TIME = Int32(1)
const IVISCOPE_WARN_INVALID_WFM_ELEMENT = IVI_CLASS_WARN_BASE + Int32(0x0001)
const IVISCOPE_ERROR_CHANNEL_NOT_ENABLED = IVI_CLASS_ERROR_BASE + Int32(0x0001)
const IVISCOPE_ERROR_UNABLE_TO_PERFORM_MEASUREMENT = IVI_CLASS_ERROR_BASE + Int32(0x0002)
const IVISCOPE_ERROR_MAX_TIME_EXCEEDED = IVI_CLASS_ERROR_BASE + Int32(0x0003)
const IVISCOPE_ERROR_INVALID_ACQ_TYPE = IVI_CLASS_ERROR_BASE + Int32(0x0004)
const IVISCOPE_WARNMSG_INVALID_WFM_ELEMENT = "One of the elements in the waveform array is invalid."
const IVISCOPE_ERRMSG_CHANNEL_NOT_ENABLED = "Specified channel is not enabled."
const IVISCOPE_ERRMSG_UNABLE_TO_PERFORM_MEASUREMENT = "Unable to perform desired measurement operation."
const IVISCOPE_ERRMSG_MAX_TIME_EXCEEDED = "Maximum time exceeded before the operation completed."
const IVISCOPE_ERRMSG_INVALID_ACQ_TYPE = "Invalid acquisition type."

# Skipping MacroDefinition: IVISCOPE_ERROR_CODES_AND_MSGS { IVISCOPE_ERROR_CHANNEL_NOT_ENABLED , IVISCOPE_ERRMSG_CHANNEL_NOT_ENABLED } , { IVISCOPE_ERROR_UNABLE_TO_PERFORM_MEASUREMENT , IVISCOPE_ERRMSG_UNABLE_TO_PERFORM_MEASUREMENT } , { IVISCOPE_ERROR_MAX_TIME_EXCEEDED , IVISCOPE_ERRMSG_MAX_TIME_EXCEEDED } , { IVISCOPE_ERROR_INVALID_ACQ_TYPE , IVISCOPE_ERRMSG_INVALID_ACQ_TYPE } , { IVISCOPE_WARN_INVALID_WFM_ELEMENT , IVISCOPE_WARNMSG_INVALID_WFM_ELEMENT }


text:
# # Near-surface backscattering
#
# Near-surface backscattering is a method of accurately calculating the backscattering from an infinite halfspace. For just the code see [backscattering.jl](backscattering.jl)
# First, let us see why it is difficult to approximate the scattering from a halfspace filled with particles. That is, let us find out how many particles are required before the backscattering converges.
#
# ## Generate a large material filled with particles.
#
using MultipleScattering
using Plots
pyplot(linewidth=2)

host_medium = Acoustic(1.0, 1.0, 2)

radius = 0.8
volfrac = 0.10
max_width = 70.

particle_medium = Acoustic(0.2, 0.1, 2)
particle_shape = Circle(radius)

bottomleft = [0.,-max_width]
topright = [max_width,max_width]

shape = Box([bottomleft,topright])
particles = random_particles(particle_medium, particle_shape; region_shape = shape, volume_fraction = volfrac, seed = 1)

# We send an incoming harmonic plane wave and receive the backscattering at `x`:
x = [-10.,0.]
source =  plane_source(host_medium; position = x,
        direction = [1.0,0.],
        amplitude = 1.0)

plot(particles)
scatter!([x[1]],[x[2]], lab="")
annotate!([(x[1], x[2] -max_width/10., "Receiver")])
plot!(shape, linecolor = :red)

# ## Calculate backscattering for different quantity of particles
# We will shave off particles on the right of this group of particles (above), and then calculate the backscattered waves for a range of angular frequencies `ωs`.
ωs = collect(0.01:0.01:1.)
widths = 10.:5.:max_width
num_particles = zeros(length(widths))

#This part below my take a while!
results = map(eachindex(widths)) do i
    bottomleft = [0.,-widths[i]]
    topright = [widths[i],widths[i]]
    shape = Box([bottomleft, topright])

    ps = filter(p -> p ⊆ shape, particles) # select particles that are inside shape
    num_particles[i] = Int(length(ps))

    simulation = FrequencySimulation(ps, source)
    run(simulation, x, ωs)
end

#save("results.jld", "$(typeof(results))",results)
#save("num_particles.jld", "$(typeof(num_particles))",num_particles)

#To load results uncomment
    # results = first(values(load("results.jld")))
    # num_particles = first(values(load("num_particles.jld")))

backscattered_waves = field.(results)

M = length(backscattered_waves)
bM = backscattered_waves[M] # backscattering from largest material
differences = [norm(b - bM) for b in backscattered_waves[1:(M-1)]]./norm(bM)

plot_converge = plot(num_particles[1:(M-1)], differences,
    xlabel = "number of particles", yguide ="error %",
    label="frequency convergence"
)
savefig("freq_convergence.png")

time_simulations = frequency_to_time.(results)
receiver = results[1].x[1]
times = 2*(widths .- receiver[1]) # time it takes for an incident plane wave to reach the furthest particles and then return to the receiver

plot()
for i in [1,3,6,9,12,13]
    plot!(time_simulations[i],label="$(num_particles[i]) particles"
        , xlims=(0,maximum(times)+10.), ylims=(-0.2,0.1)
        , xticks = [0.; 30.; times]
    )
end
gui()
savefig("time_response.png")

time_vec = 0.:pi:34.2
time_results = frequency_to_time.(results; t_vec = time_vec, impulse = GaussianImpulse(maximum(ωs)))

backscattered_waves = field.(time_results)
bM = backscattered_waves[M] # backscattering from largest material
differences = [norm(b - bM) for b in backscattered_waves[1:(M-1)]]./norm(bM)
plot(plot_converge)
plot!(num_particles[1:(M-1)], differences, xlabel = "number of particles", yguide ="error %", label="time convergence")
savefig("compare_convergence.png")

## Using near surface backscattering
shape = TimeOfFlight(receiver,80.0)
scatter([receiver[1]],[receiver[2]]);
annotate!([(receiver[1], receiver[2] -max_width/10., "Receiver")])
plot!(particles);
plot!(shape, linecolor=:red)

savefig("time_of_flight_shape.png")

times = 40.:15.:80.
near_surface_simulations = map(times) do t
    shape = TimeOfFlight(receiver,t) # choose a material with particles only in the near surface region
    ps = filter(p -> p ⊆ shape, particles) # select particles that are inside shape
    run(FrequencySimulation(ps, source), x, ωs) # calculate backscattering
end
save("near_surface_simulations.jld","Array{FrequencySimulation{Float64},1}",near_surface_simulations)

time_near_simulations = frequency_to_time.(near_surface_simulations; impulse = GaussianImpulse(maximum(ωs)))

plot()
for i in 1:length(times)
    plot!(time_near_simulations[i],label="time of flight $(times[i])"
        , xlims=(0,maximum(times)+10.), ylims=(-0.6,0.3)
        , xticks = [0.; times], title=""
    )
end
gui()

savefig("time_of_flight_response.png")