text:
# add AprilTag sightings from this pose
function addApriltags!(fg, pssym, posetags; bnoise=0.1, rnoise=0.1, lmtype=Point2, fcttype=Pose2Pose2, DAerrors=0.0, autoinit=true )
  @show currtags = ls(fg)[2]
  for lmid in keys(posetags)
    @show lmsym = Symbol("l$lmid")
    if !(lmsym in currtags)
      @info "adding node $lmsym"
      addVariable!(fg, lmsym, lmtype)
    end
    ppbr = nothing
    if lmtype == RoME.Point2
      ppbr = Pose2Point2BearingRange(Normal(posetags[lmid][:bearing][1],bnoise),
                                     Normal(posetags[lmid][:range][1],rnoise))
    elseif lmtype == Pose2
      dx, dy = posetags[lmid][:bP2t].translation[1], posetags[lmid][:bP2t].translation[2]
      dth = convert(RotXYZ, posetags[lmid][:bP2t].linear).theta3
      # ppbr = Pose2Pose2(
      #           MvNormal([ posetags[lmid][:pos][3];
      #                     -posetags[lmid][:pos][1];
      #                      posetags[lmid][:tRYc]],
      #                    diagm([0.1;0.1;0.05].^2)) )
      ppbr = fcttype(
                MvNormal([dx;
                          dy;
                          dth],
                         Matrix(Diagonal([0.05;0.05;0.03].^2))) )
    end
    if rand() > DAerrors
      # regular single hypothesis
      addFactor!(fg, [pssym; lmsym], ppbr, autoinit=autoinit)
    else
      # artificial errors to data association occur
      info("Forcing bad data association with $lmsym")
      xx,ll = ls(fg)
      @show ll2 = setdiff(ll, [lmsym])
      @show daidx = round(Int, (length(ll2)-1)*rand()+1)
      @show rda = ll2[daidx]
      addFactor!(fg, [pssym; lmsym; rda], ppbr, multihypo=[1.0;0.5;0.5], autoinit=autoinit)
    end
  end
  nothing
end

function addnextpose!(fg, prev_psid, new_psid, pose_tag_bag; lmtype=Point2, odotype=Pose2Pose2, fcttype=Pose2Pose2, DAerrors=0.0, autoinit=true)
  prev_pssym = Symbol("x$(prev_psid)")
  new_pssym = Symbol("x$(new_psid)")
  # first pose with zero prior
  if odotype == Pose2Pose2
    addVariable!(fg, new_pssym, Pose2)
    addFactor!(fg, [prev_pssym; new_pssym], Pose2Pose2(MvNormal(zeros(3),diagm([0.4;0.1;0.4].^2))), autoinit=autoinit)
  elseif odotype == VelPose2VelPose2
    addVariable!(fg, new_pssym, DynPose2(ut=round(Int, 200_000*(new_psid))))
    addFactor!(fg, [prev_pssym; new_pssym], VelPose2VelPose2(MvNormal(zeros(3),Matrix(Diagonal([0.4;0.4;0.3].^2))),
                                                             MvNormal(zeros(2),Matrix(Diagonal([0.2;0.1].^2)))), autoinit=autoinit)
  end

  addApriltags!(fg, new_pssym, pose_tag_bag, lmtype=lmtype, fcttype=fcttype, DAerrors=DAerrors)
  new_pssym
end

function prepCamLookup(imgseq; filenamelead="camera_image")
  # camcount = readdlm(datafolder*"$(imgfolder)/cam-count.csv",',')
  camlookup = Dict{Int, String}()
  count = -1
  for i in imgseq
    count += 1
    camlookup[count] = "$(filenamelead)$(i).jpeg"
    # camlookup[camcount[i,1]] = strip(camcount[i,2])
  end
  return camlookup
end

function detectTagsViaCamLookup(camlookup, imgfolder, imgsavedir)
  # AprilTag detector
  detector = AprilTagDetector()

  # extract tags from images
  IMGS = []
  TAGS = []
  psid = 1
  for psid in 0:(length(camlookup)-1)
    img = load("$(imgfolder)/$(camlookup[psid])")
    tags = detector(img)
    push!(TAGS, deepcopy(tags))
    push!(IMGS, deepcopy(img))
    foreach(tag->drawTagBox!(IMGS[psid+1],tag, width = 5, drawReticle = false), tags)
    Images.save(imgsavedir*"/tags/img_$(psid).jpg", IMGS[psid+1])
  end
  # psid = 1
  # img = load(datafolder*"$(imgfolder)/$(camlookup[psid])")

  # free the detector memory
  freeDetector!(detector)

  return IMGS, TAGS
end

function results2csv(fg; dir=".", filename="results.csv")

if !isdir(dir*"/results")
  mkdir(dir*"/results")
end
fid = open(dir*"/results/"*filename,"w")
for sym in [ls(fg)[1]...;ls(fg)[2]...]
  p = getVertKDE(fg, sym)
  val = string(KDE.getKDEMax(p))
  println(fid, "$sym, $(val[2:(end-1)])")
end
close(fid)

end

## BUILD FACTOR GRAPH FOR SLAM SOLUTION,

function main(resultsdir::String,
              camidxs,
              tag_bagl;
              BB=20,
              N=100,
              lagLength=75,
              dofixedlag=true,
              jldfile::String="",
              failsafe::Bool=false,
              show::Bool=false  )

# Factor graph construction
fg = initfg()
prev_psid = 0

# load from previous file
if jldfile != ""
  fg, = loadjld( file=joinpath(resultsdir,jldfile) )
  xx, = ls(fg)
  prev_psid = parse(Int, string(xx[end])[2:end])
end

fg.solverParams.isfixedlag = dofixedlag
fg.solverParams.qfl = lagLength

psid = 0
pssym = Symbol("x$psid")
# first pose with zero prior
addVariable!(fg, pssym, DynPose2(ut=0))
# addFactor!(fg, [pssym], PriorPose2(MvNormal(zeros(3),diagm([0.01;0.01;0.001].^2))))
addFactor!(fg, [pssym], DynPose2VelocityPrior(MvNormal(zeros(3),Matrix(Diagonal([0.01;0.01;0.001].^2))),
                                              MvNormal(zeros(2),Matrix(Diagonal([0.1;0.05].^2)))))

addApriltags!(fg, pssym, tag_bagl[psid], lmtype=Pose2, fcttype=DynPose2Pose2)

show ? writeGraphPdf(fg) : nothing

# quick solve as sanity check
tree = batchSolve!(fg, N=N, drawpdf=true, show=show, recursive=failsafe)

# add other positions
maxlen = (length(tag_bagl)-1)

Gadfly.push_theme(:default)

@sync begin

for psid in (prev_psid+1):1:maxlen
  prev_psid
  maxlen
  @show psym = Symbol("x$psid")
  addnextpose!(fg, prev_psid, psid, tag_bagl[psid], lmtype=Pose2, odotype=VelPose2VelPose2, fcttype=DynPose2Pose2, autoinit=true)
  # writeGraphPdf(fg)

  if psid % BB == 0 || psid == maxlen
    IIF.savejld(fg, file=resultsdir*"/racecar_fg_$(psym)_presolve.jld2")
    tree = batchSolve!(fg, drawpdf=true, show=show, N=N, recursive=true)
  end
  jldfile = resultsdir*"/racecar_fg_$(psym).jld2"
  T1 = @spawn IIF.savejld(fg, file=jldfile)
  @async fetch(T1)

  ## save factor graph for later testing and evaluation
  # ensureAllInitialized!(fg)
  T2 = @spawn plotRacecarInterm(fg, resultsdir, psym)
  @async fetch(T2)

  # prepare for next iteration
  prev_psid = psid
end # for

# extract results for later use as training data
results2csv(fg, dir=resultsdir, filename="results.csv")

end #sync

return fg
end # main

#


text:
function cost(solver::Solver,X::Matrix,U::Matrix,weights::Vector,method::Symbol)
    obj = solver.obj
    f = solver.model.f
    Q = obj.cost.Q; xf = obj.xf; Qf = obj.cost.Qf; R = obj.cost.R;
    n,m = get_sizes(solver)
    m̄, = get_num_controls(solver)
    N,N_ = get_N(solver,method)
    dt = solver.dt

    if method == :hermite_simpson
        # pull out solver/objective values
        J = zeros(eltype(X),N-1)
        f1 = zeros(eltype(X),n)
        f2 = zeros(eltype(X),n)
        for k = 1:N-1
            f(f1,X[:,k], U[1:m,k])
            f(f2,X[:,k+1], U[1:m,k+1])
            x1 = X[:,k]
            x2 = X[:,k+1]

            solver.state.minimum_time ? dt = U[m̄,k] : nothing

            Xm = (x1+x2)/2 + dt/8*(f1-f2)
            Um = (U[1:m,k] + U[1:m,k+1])/2

            J[k] = dt/6*(stage_cost(obj.cost,X[:,k],U[1:m,k]) + 4*stage_cost(obj.cost,Xm,Um) + stage_cost(obj.cost,X[:,k+1],U[1:m,k+1])) # rk3 foh stage cost (integral approximation
            solver.state.minimum_time ? J[k] += solver.opts.R_minimum_time*dt : nothing
        end
        J = sum(J)
        J += stage_cost(obj.cost,X[:,N])
        return J#, XM, UM, fVal
    elseif method == :midpoint
        Xm = zeros(eltype(X),n,N-1)
        for k = 1:N-1
            Xm[:,k] = (X[:,k] + X[:,k+1])/2
        end
        J = zeros(eltype(Xm),N)
        for k = 1:N-1
            J[k] = stage_cost(obj.cost,Xm[:,k],U[:,k])
        end
        J = weights'J
        J += stage_cost(obj.cost,X[:,N])
        return J

    else
        J = zeros(eltype(X),N)
        for k = 1:N
            J[k] = stage_cost(obj.cost,X[:,k],U[:,k])
        end
        J = weights'J
        J += stage_cost(obj.cost,X[:,N])
        return J
    end
end

function cost_gradient(solver::Solver, res::DircolResults, method::Symbol)
    n,m = get_sizes(solver)
    N,N_ = get_N(solver,method)
    grad_f = zeros((n+m)N)
    cost_gradient!(solver,res.X_,res.U_,res.fVal_,res.A,res.B,res.weights,grad_f,method)
    return grad_f
end

function cost_gradient!(solver::Solver, vars::DircolVars, weights::Vector{Float64}, vals::Vector{Float64}, method::Symbol)::Nothing
    # println("simple")
    n,m = get_sizes(X,U)
    N,N_ = get_N(solver,method)
    dt = solver.dt

    obj = solver.obj
    Q = obj.cost.Q; xf = obj.xf; R = obj.cost.R; Qf = obj.cost.Qf;
    # X,U = res.X_, res.U_
    grad_f = reshape(vals, n+m, N)

    for k = 1:N
        grad_f[1:n,k] = weights[k]*Q*(X[:,k] - xf)*dt
        grad_f[n+1:end,k] = weights[k]*R*U[:,k]*dt
    end
    grad_f[1:n,N] += Qf*(X[:,N] - xf)
    return nothing
end

function collocation_constraints(solver::Solver, res::DircolResults, method::Symbol)
    n,m = get_sizes(solver)
    N,N_ = get_N(solver,method)
    g = zeros(eltype(res.X),(N-1)*n)
    collocation_constraints!(solver,res.X_,res.U_,res.fVal_,g,method)
    return g
end

function collocation_constraints(X,U,method,dt,f::Function)
    n,m,N = get_sizes(X,U)
    g = zeros(eltype(X),(N-1)*n)
    g = reshape(g,n,N-1)

    fVal = zeros(eltype(X),n,N)
    for k = 1:N
        f(view(fVal,:,k),X[:,k],U[:,k])
    end

    if method == :trapezoid
        for k = 1:N-1
            # Collocation Constraints
            g[:,k] = dt*( fVal[:,k+1] + fVal[:,k] )/2 - X[:,k+1] + X[:,k]
        end
    elseif method == :hermite_simpson_separated
        iLow = 1:2:N-1
        iMid = iLow .+ 1
        iUpp = iMid .+ 1

        midpoints =  - X[:,iMid] + (X[:,iLow]+X[:,iUpp])/2 + dt*(fVal[:,iLow] - fVal[:,iUpp])/8
        collocation = - X[:,iUpp] + X[:,iLow] + dt*(fVal[:,iLow] + 4*fVal[:,iMid] + fVal[:,iUpp])/6
        g[:,iLow] = collocation
        g[:,iMid] = midpoints

    elseif method == :hermite_simpson
        fm = zeros(eltype(X),n)
        for k = 1:N-1
            x1 = X[:,k]
            x2 = X[:,k+1]
            xm = (x1+x2)/2 + dt/8*(fVal[:,k]-fVal[:,k+1])
            um = (U[:,k] + U[:,k+1])/2
            f(fm, xm, um)
            g[:,k] = -x2 + x1 + dt*(fVal[:,k] + 4*fm + fVal[:,k+1])/6
        end

    elseif method == :midpoint
        fm = zeros(eltype(X),n)
        for k = 1:N-1
            x1 = X[:,k]
            x2 = X[:,k+1]
            xm = (x1+x2)/2
            f(fm,xm,U[:,k])
            g[:,k] = dt*fm - x2 + x1
        end
    end
    return vec(g)
end

ext: jl


text:
## Optimization Model for dtu42380a1
## Mar 11, 2020

# function get_data(args)
L = 3; # number of locations
P = 3; # number of products
T = 4; # number of periods
S = 11;  # number of materials

pi = [
    0.7	0	0	1.2	0.3 1	0	0	3	2	0
    0	0.6	0	1.2	0   0	1	0	3	2	1
    0	0	0.8	1.2	0.3 0	0	1	3	2	1
]; # pi[p,s] gives amount of material s needed to produce one unit of product p

muB = [4 5 6]; # muB[p] gives the number of seconds of preblend machine time needed to produce one unit of product p
muP = [2 2 2]; # muP[p] gives the number of seconds of packaging machine time needed to package one unit of product p

c = [4 5 4 10 3 20 20 20 15 15 30]; # c[s] gives the cost for sourcing one unit of material s

cBR = 0.07; # Blending machine regular time cost per second
cBO = 0.18; # Blending machine overtime cost per second

cPR = 0.058; # Packaging machine regular time cost per second
cPO = 0.075; # Packaging machine overtime cost  per second

h = [0.2 0.3 0.3]; # h[p] gives the inventory holding cost for product p

## I0[ip] gives the initial inventory of product p at location i
I0 = [
    2368	1834	1940
    1743	1500	2910
    4307	1799	3100
];

## D[i,p,t] gives the demand forecast for product p at location i in period t
D = zeros(L,P,T);
D[1,:,:] = [
    40539	40997	39855	38574
    20854	21562	21925	22315
    9502	8975	9347	9958
];
D[2,:,:] = [
    55948	55804	55292	53866
    38245	39364	39409	39380
    16888	16617	14462	13666
];
D[3,:,:] = [
    49559	50155	48959	46259
    27461	27718	26460	27059
    6738	10495	11565	13377
];

## ss[i,p,t] gives the safety stock needed of product p at location i in period t
ss = [
    2110	1243	814
    1650	1488	3046
    4172	1734	3014
];
# make safety stock for the rest of the periods the same
for t = 2:T
    global ss = cat(ss,ss,dims = 3);
end

## penalty for not satisfying one unit of demand
PI = 1000000;


text:
context("Not installed") do
	eval(Base, :(is_interactive = true)) # for interactive error handling

	add_format(format"NotInstalled", (), ".not_installed", [:NotInstalled])
	stdin_copy = STDIN
    stderr_copy = STDERR
    rs, wr = redirect_stdin()
	rserr, wrerr = redirect_stderr()
	ref = @async save("test.not_installed")
	println(wr, "y")
	if VERSION < v"0.4.0-dev"
		@fact_throws ErrorException wait(ref) #("unknown package NotInstalled")
	else
		@fact_throws CompositeException wait(ref) #("unknown package NotInstalled")
	end
	ref = @async save("test.not_installed")
	println(wr, "invalid") #test invalid input
	println(wr, "n") # don't install
	wait(ref)
	@fact istaskdone(ref) --> true

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

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

end


text:
# This takes over from some of the lower-level tests that were done earlier to test
# the newcommand; these are more integrated as a result.

@testset "basic" begin
    s = raw"""
        \newcommand{\abc}{hello}
        A \abc B
        """ |> fd2html
    @test s // "
A hello B
"
    s = raw"""
        \newcommand{\abc}[1]{hello#1}
        A \abc{Bob} C
        """ |> fd2html
    @test s // "
A hello Bob C
"
    s = raw"""
        \newcommand{\abc}[1]{hello!#1}
        A \abc{Bob} C
        """ |> fd2html
    @test s // "
A helloBob C
"
    s = raw"""
        \newcommand{\abc}[2]{A#1B#2C}
        AA \abc{EE}{FF} BB
        """ |> fd2html
    @test s // "
AA A EEB FFC BB
"
    s = raw"""
        \newcommand{\abc}[2]{A!#1B!#2C}
        AA \abc{EE}{FF} BB
        """ |> fd2html
    @test s // "
AA AEEBFFC BB
"
end

@testset "redef" begin
    s = raw"""
        \newcommand{\abc}{hello}
        A \abc B
        \newcommand{\abc}{bye}
        A \abc B
        """ |> fd2html
    @test s // "
A hello B
\n
A bye B
"
    s = raw"""
        \newcommand{\abc}[1]{hello#1}
        A \abc{Bob} B
        \newcommand{\abc}{bye}
        A \abc B
        """ |> fd2html
    @test s // "
A hello Bob B
\n
A bye B
"
end

@testset "nesting" begin
    s = raw"""
        \newcommand{\abc}{hello}
        \newcommand{\def}[1]{\abc#1}
        A \def{Bob} B
        """ |> fd2html
    @test s // "
A hello Bob B
"
end

@testset "errors" begin
    # malformed
    s = raw"""
        abc \newcommand abc
        """
    @test_throws F.LxObjError (s |> fd2html)
    s = raw"""
        \newcommand{abc}{hello}
        """
    @test_throws F.LxDefError (s |> fd2html)
    # ordering
    s = raw"""
        \abc
        \newcommand{\abc}{hello}
        """
    @test_throws F.LxObjError (s |> fd2html)
    # nargs not proper
    s = raw"""
        \newcommand{\abc}[d]{hello}
        """
    @test_throws F.LxDefError (s |> fd2html)
    # bad use
    s = raw"""
        \newcommand{\comb}[1]{HH#1HH}
        etc \comb then.
        """
    @test_throws F.LxObjError (s |> fd2html)
end

# ==== more integrated examples recuperated from previous more specific tests ===

@testset "integrated-1" begin
    st = raw"""
        \newcommand{\com}{blah}
        \newcommand{\comb}[ 2]{hello #1 #2 \com}
        A \comb{AA}{BB} B
        """ |> fd2html
    @test isapproxstr(st, "
A hello AA BB blah B
")
end

@testset "integrated-2" begin
    st = raw"""
        \newcommand{\com}{HH}
        \newcommand{\comb}[1]{GG#1FF}
        Blah \com and \comb{blah} etc
        ```julia
        f(x) = x^2
        ```
        etc \comb{blah} then maybe
        @@adiv inner part @@ final.
        """ |> fd2html
    @test isapproxstr(st, """

Blah HH and GG blahFF etc

f(x) = x^2

etc GG blahFF then maybe

inner part

final.

        """)
end

@testset "integrated-3" begin
    st = raw"""
        \newcommand{\E}[1]{\mathbb E\left[#1\right]}blah de blah
        ~~~
        escape b1
        ~~~
        \newcommand{\eqa}[1]{\begin{eqnarray}#1\end{eqnarray}}
        \newcommand{\R}{\mathbb R}
        Then something like
        \eqa{ \E{f(X)} \in \R &\text{if}& f:\R\maptso\R }
        and we could try to show latex:
        ```latex
        \newcommand{\brol}{\mathbb B}
        ```
        """ |> fd2html
    @test isapproxstr(st, raw"""

          blah de blah
          escape b1

          Then something like
            \[\begin{array}{rcl}
              \mathbb E\left[ f(X)\right] \in \mathbb R &\text{if}& f:\mathbb R\maptso\mathbb R
            \end{array}\]
          and we could try to show latex:

          \newcommand{\brol}{\mathbb B}

        """)
end

@testset "integrated-4" begin
    s = raw"""
        text A1 \newcommand{\com}{blah}text A2 \com and
        ~~~
        escape B1
        ~~~
        \newcommand{\comb}[ 1]{\mathrm{#1}} text C1 $\comb{b}$ text C2
        \newcommand{\comc}[ 2]{part1:#1 and part2:#2} then \comc{AA}{BB}.
        """ |> fd2html
    @test isapproxstr(s, raw"""

text A1 text A2 blah and
        escape B1
        text C1 \(\mathrm{ b}\) text C2  then part1: AA and part2: BB.

        """)
end

@testset "integrated-5" begin
    s = raw"""
        @def title = "Convex Optimisation I"
        \newcommand{\com}[1]{⭒!#1⭒}
        \com{A}

        then some
        ## blah
        end \com{B}.
        """ |> fd2html
    @test isapproxstr(s, raw"""
        ⭒A⭒

then some

blah

end ⭒B⭒.

        """)
end

@testset "integrated-6" begin
    st = raw"""
        \newcommand{\com}{HH}
        \newcommand{\comb}[1]{FF#1GG}
        A list
        * \com and \comb{blah}
        * $f$ is a function
        * a last element
        """ |> fd2html
    @test isapproxstr(st, raw"""

A list

HH and FF blahGG

\(f\) is a function

a last element

        """)

    st = raw"""
        a\newcommand{\eqa}[1]{\begin{eqnarray}#1\end{eqnarray}}b@@d .@@
        \eqa{\sin^2(x)+\cos^2(x) &=& 1}
        """ |> fd2html

    @test isapproxstr(st, raw"""

ab

.

        \[\begin{array}{rcl} \sin^2(x)+\cos^2(x) &=& 1\end{array}\]
        """)
end


text:
function make_clauses()
	outfile = open("instance.txt", "w")

	#srand(12345)

	#top = 70 #rand(1:40)
	#k = binomial(top,2)

	k = rand(20:1000)
	a = 1
	b = -1
	c = -k
	lowerbound = int((-b + sqrt(b^2 - 4*a*c))/(2a))
	upperbound = 2k
	println("Lower bound = ", lowerbound)
	println("Upper bound = ", upperbound)
	println("K = ", k)
	for i = 1:k
		x1 = rand(lowerbound:upperbound)
		x2 = rand(lowerbound:upperbound)

		if(randbool())
			x1=-x1
		end

		if(randbool())
			x2=-x2
		end

		println(outfile, join((x1, x2), " "))
	end
	close(outfile)
end

max_stars_repo_path: Julia Implemenation/makeclauses.jl

max_issues_repo_path: Julia Implemenation/makeclauses.jl

max_forks_repo_path: Julia Implemenation/makeclauses.jl


text:
module QuantumLab

lib_path = joinpath(dirname(@__FILE__),"..","deps","usr","lib")
push!(Libdl.DL_LOAD_PATH,lib_path)
try
  Libdl.dlopen("libint2-QuantumLab.so")
  global const libint2_available = true
catch e
  warn("Caught an error when trying to dlopen libint2-QuantumLab.so \n
           You can try to rebuild with Pkg.build(\"QuantumLab\"). \n
	   Alternatively, you can set the symlink by hand if you know where the
	   dynamic library is located. Until then we'll expect you can't access
	   libint2 and continue with a stub to get the code running anyways.
	   This is typical behavior under Windows and leads to a drastic performance decrease.")
  global const libint2_available = false
end

print(" + (ProgressMeter................")
using ProgressMeter
ProgressMeter.printover(STDOUT," + (Reexport...................")
using Reexport
ProgressMeter.printover(STDOUT,"")

print(" + DocumentationModule.........")
include("SubModules/DocumentationModule.jl")
@reexport using .DocumentationModule
println("Done.")

print(" + BaseModule..................")
include("SubModules/BaseModule.jl")
@reexport using .BaseModule
println("Done.")

print(" + AtomModule..................")
include("SubModules/AtomModule.jl")
@reexport using .AtomModule
println("Done.")

print(" + GeometryModule..............")
include("SubModules/GeometryModule.jl")
@reexport using .GeometryModule
println("Done.")

print(" + BasisSetModule..............")
include("SubModules/BasisSetModule.jl")
@reexport using .BasisSetModule
println("Done.")

print(" + BasisFunctionsModule........")
include("SubModules/BasisFunctionsModule.jl")
@reexport using .BasisFunctionsModule
println("Done.")

print(" + ShellModule.................")
include("SubModules/ShellModule.jl")
@reexport using .ShellModule
println("Done.")

print(" + IntegralsModule.............")
include("SubModules/IntegralsModule.jl")
@reexport using .IntegralsModule
println("Done.")

print(" + BasisModule.................")
include("SubModules/BasisModule.jl")
@reexport using .BasisModule
println("Done.")

print(" + LibInt2Module...............")
include("SubModules/LibInt2Module.jl")
@reexport using .LibInt2Module
println("Done.")

print(" + SpecialMatricesModule.......")
include("SubModules/SpecialMatricesModule.jl")
@reexport using .SpecialMatricesModule
println("Done.")

print(" + HartreeFockModule...........")
include("SubModules/HartreeFockModule.jl")
@reexport using .HartreeFockModule
println("Done.")

print(" + InitialGuessModule..........")
include("SubModules/InitialGuessModule.jl")
@reexport using .InitialGuessModule
println("Done.")

print(" + BasisSetExchangeModule......")
include("SubModules/BasisSetExchangeModule.jl")
@reexport using .BasisSetExchangeModule
println("Done.")

print(" + LaplaceModule...............")
include("SubModules/LaplaceModule.jl")
@reexport using .LaplaceModule
println("Done.")

print(" + ResolutionIdentityModule....")
include("SubModules/ResolutionIdentityModule.jl")
@reexport using .ResolutionIdentityModule
println("Done.")

end # module


text:
#module DataSet

# using AbstractFFTs
using FFTW
using Statistics

export Data1D,ind2pos,pos2ind,val,ind,cut
export FourierTransform,PhaseCorrect,BaseLineCorrect
export integral,derivative,set!
export plot,plot!
export length,shift,hard_shift
export interp
export medianBaseline

ENV["GKS_ENCODING"]="utf-8"
import Plots

import Base.iterate

#import NMR.SimplePlot

# Data structures for storing digitised data sets
# Definitions:
#
# - Index: the horizontal index of a data point
# - Position: the numerical position (integer) of the data point in
#           the internal storage array

mutable struct Data1D{Tdata,Tindex}
   dat::Array{Tdata,1}
   istart::Tindex
   istop::Tindex
end

Base.iterate(d::NMR.Data1D,state=1) =
    state > length(d) ? nothing :
      ((d.istart+state*(d.istop-d.istart)/(length(d)-1),d.dat[state]),state+1)

function ind2pos(d::Data1D,ind)
    return round(Int64,(ind-d.istart)*(length(d.dat)-1)/(d.istop-d.istart))+1
end

function pos2ind(d::Data1D,pos::Integer)
   return ind(d)[pos]
end

function ind(d::Data1D)
	return LinRange(d.istart,d.istop,length(d.dat))
end

function val(d::Data1D)
    return d.dat
end

function val(d::Data1D,ind)
    i1=(ind-d.istart)*length(d.dat)/(d.istop-d.istart)
    p1=floor(Int64,i1)+1
    ic=i1-p1+1
    return ((1-ic)*d.dat[p1]+ic*d.dat[p1+1])
end

function cut(d::Data1D,i1,i2)
    p1=ind2pos(d,i1)
    p2=ind2pos(d,i2)
    return Data1D(d.dat[p1:p2],pos2ind(d,p1),pos2ind(d,p2))
end

resolution(d::Data1D) = (d.istop-d.istart)/length(d.dat)

"""
`FourierTransform(d::Data1D)` performs an FFT assuming that the
input data set `d` is a free induction decay, and returns
a `DataSet` object with the resulting complex data.
"""
function FourierTransform(fid::Data1D;
        CTR=0.0,
        LB=0.0,
        SI=length(fid.dat),
        PPM=1,FFTSHIFT=true)

    apod=[exp(-t*LB) for t=ind(fid)];
    s=zeros(Complex{Float64},SI);
    s[1:length(apod)]=apod.*fid.dat;
    s[1]*=0.5;
    SW=length(fid.dat)/(fid.istop-fid.istart);
    spect=fft(s);
    if(FFTSHIFT)
        spect=fftshift(spect)
    end;
    return Data1D(spect,CTR-SW/2/PPM,CTR+SW/2/PPM)
end

"""
`PhaseCorrect(spect::Data1D;Ph0=0.0,Ph1=0.0,Pivot=0.0)`
returns a phase corrected `DataSet`, using the values indicated.
The parameters are given in rad, k, and rad/k, where k are the units
of the spectral domain (typically, Hz or ppm).
"""
function PhaseCorrect(spect::Data1D;Ph0=0.0,Ph1=0.0,Pivot=0.0)
    out=deepcopy(spect);
    out.dat = spect.dat .* [exp(im*(Ph0+(k-Pivot)*Ph1)) for k=ind(spect)]
    return out;
end

polyfit(x::Vector, y::Vector, deg::Int) = pinv(collect(v ^ p for v in x, p in 0:deg)) * y

function horner(x::Vector, c::Vector)
	r=ones(length(x))*c[end]
    for k=(length(c)-1):-1:1
        r.*=x
        r.+=c[k]
    end
    return r
end

@doc raw"""
## Automatic Baseline Correction
`BaseLineCorrect(spect::Data1D;regions=128,kfactor=5,wdw=32)` corrects the
base line of `spect` using the algorithm
by S. Golotvin, A. Williams, J. Magn. Reson. 146 (2000) 122-125.

### Optional Parameters:
- `regions`: the number of regions the spectrum is divided into to estimate
   the noise variance. This should be chosen large enough so that at least one
   of these regions is free from signals, and shows a flat baseline. The
   variance ``\\sigma^2_{min}`` in this region is then used to distinguish signals and noise.

- `wdw`: the window size around each data point that is used to decide whether the
  data point is part of a signal, or should be counted as base line.

- `kfactor`: a point will be counted as part of the baseline if and only if the local
  variance in its window is less than `kfactor`*``\\sigma^2_{min}``.

- `output`: defaults to `:spectrum`. `output=:points` will return an array of points,
  which indicate baseline positions (locations free from peaks) as determined
  by the algorithm. This is useful in situations where a large number of similar
  spectra need to be processed.

### Notes
- if the input data is complex, the imaginary part is ignored,
  and the result is always real.

- the first and last point of `spect` are always used as part of the baseline.
  The first and last points of the output are therefore guaranteed to be zero. This
  helps to keep the interpolation numerically stable.

## Baseline Correction With Known Locations

`BaseLineCorrect(spect::Data1D,pts::Array{Float64,1})`
uses the points in `pts` to determine the baseline. This should be a list
of positions in the spectrum that are known not to contain signal; it can be obtained
by using the option `output=:points` in the first method of `BaseLineCorrect()`.

"""
function BaseLineCorrect(spect::Data1D;regions=128,kfactor=5,lw=0.0,wdw=32,order=5,verbose=false,output=:spectrum)

    out=Data1D(real.(spect.dat),spect.istart,spect.istop);
    n=length(out.dat)
    lreg=floor(Int64,n/regions)

    # compute wdw from lw parameter if set
    if lw>0.0
        wdw = round(Int32,4*lw/resolution(spect))
    end

    # find smallest standard deviation in regions
    mindev=minimum(abs.([std(out.dat[k:k+lreg]) for k=1:lreg:(n-lreg)]))

    # find base line points
    select=[abs(std(out.dat[(k-wdw):(k+wdw)])) < kfactor*mindev for k=(wdw+1):(n-wdw)]

    pts=((wdw+1):(n-wdw))[select];
    verbose && println("Fitting $(length(pts)) Points.")
    push!(pts,1);
    push!(pts,n);

    # polynomial fitting
    coeffs=polyfit(pts/n,out.dat[pts],order)
    verbose && println("Coefficients: $coeffs")

    res=sum(abs.(out.dat[pts]-horner(pts/n,coeffs)))
    verbose && println("Residual: $res")

    out.dat .-= horner(collect(1:n)/n,coeffs)

    if output==:points
        return([pos2ind(spect,p) for p in pts])
    else
        return (out)
    end
end

function BaseLineCorrect(spect::Data1D,ind::Array{Float64,1};order=5)
    out=Data1D(real.(spect.dat),spect.istart,spect.istop);
    n=length(out.dat)
    pts=[ind2pos(out,k) for k in ind]
    coeffs=polyfit(pts/n,out.dat[pts],order)
    out.dat .-= horner(collect(1:n)/n,coeffs)
    return(out)
end

wrap(n,l)=[mod(k,l) for k in n]

function extrema(X::Array{Float64,1})
    Y=Array{Float64,1}()
    for k=2:(length(X)-1)
        if (X[k]>X[k-1] && X[k]>X[k+1]) || X[k]<X[k-1] && X[k]<X[k+1]
            push!(Y,X[k])
        end
    end
    return Y
end

@doc raw"""
    function medianBaseline(s::NMR.Data1D;wdw=length(s)>>6)

Compute baseline for the real part of `s` by the algorithm of M. S. Friedrichs,
*Journal of Biomolecular NMR*,  **5** (1995) 147  153.
The window is automatically calculated as ``\Delta\omega/2^6``, unless
overridden by giving a value to the optional key `wdw`.

"""
function medianBaseline(s::NMR.Data1D{Tdata,Tindex};wdw=length(s)>>6) where { Tdata<:Real, Tindex}
    r=s.dat
    L=length(r)
    b=zeros(L)
    c=zeros(L)
    for k=1:L
        b[k]=median(extrema(r[ wrap(k.+(-wdw:wdw),1:L)]))
    end
    gauss=exp.(-((wdw:wdw)./wdw).^2)
    gauss=gauss/sum(gauss)
    for k=1:length(r)
        c[k]=sum(gauss.*b[wrap(k.+(-wdw:wdw),1:L)])/(2*wdw+1)
    end
    return NMR.Data1D(c,s.istart,s.istop)
end

function medianBaseline(s;opts...)
    return medianBaseline(real(s))+im*medianBaseline(imag(s);opts...)
end

"""
`integral(spect::Data1D)` computes the numerical integral of `spect` by
the trapezoid rule.
"""
function integral(spect::Data1D)
    s=sum(spect.dat[2:(end-1)])+0.5*spect.dat[1]+0.5*spect.dat[end]
    return s*(spect.istop-spect.istart)/length(spect.dat)
end

"""
`derivative(spect::Data1D)` computes the derivative of `spect` and returns
it as another Data1D object
"""
function derivative(spect::Data1D)
    s=spect.dat
    inc=(spect.istop-spect.istart)/length(spect.dat)
#    d[1:(end-1)]=(spect.dat[2:end]-spect.dat[1:(end-1)] )/inc
#    d[end]=d[end-1]
    d = 1.0/12*(8*[s[2:end];0]-8*[0;s[1:(end-1)]] + [s[3:end];0;0] - [0;0;s[1:(end-2)]] )/inc
    return Data1D(d,spect.istart,spect.istop)
end

function Plot(s::Data1D;opts...)
    return Plot(ind(s),real(val(s)),style=Dict(["stroke-width"=>"1"]),Reverse=[true,false];opts...)
end

function set!(spect::Data1D,start,stop,value=0)
    rge=ind2pos(spect,start):ind2pos(spect,stop)
    spect.dat[rge]=value
end

import Base.*, Base./, Base.+,Base.-,Base.abs

function *(d::Data1D,n::Number)
    c = n * d.dat;
    return Data1D(c,d.istart,d.istop)
end

*(n::Number,d::Data1D)=d*n

function *(d1::Data1D,d2::Data1D)
	ind(d1) == ind(d2) || error("Data1D objects incompatible for multiplication")
	return Data1D(d1.dat.*d2.dat,d1.istart,d1.istop)
end

function /(d::Data1D,n::Number)
    return Data1D(d.dat./n,d.istart,d.istop)
end

function /(d1::Data1D,d2::Data1D)
	ind(d1) == ind(d2) || error("Data1D objects incompatible for division")
	return Data1D(d1.dat./d2.dat,d1.istart,d1.istop)
end

function +(d::Data1D,n::Number)
    return Data1D(n.+d.dat,d.istart,d.istop)
end

+(n::Number,d::Data1D)=d+n

function +(d1::Data1D,d2::Data1D)
	ind(d1) == ind(d2) || error("Data1D objects incompatible for addition")
	return Data1D(d1.dat.+d2.dat,d1.istart,d1.istop)
end

function -(d::Data1D,n::Number)
    return Data1D(-n.+d.dat,d.istart,d.istop)
end

function -(d1::Data1D,d2::Data1D)
	ind(d1) == ind(d2) || error("Data1D objects incompatible for subtraction")
	return Data1D(d1.dat.-d2.dat,d1.istart,d1.istop)
end

function abs(d::Data1D)
    c=Data1D(abs.(d.dat),d.istart,d.istop)
    return c
end

function plot(d::Data1D;opts...)
	return Plots.plot(real.(ind(d)),real.(val(d));opts...)
end

import Plots.RecipesBase.plot!

function plot!(d::Data1D;opts...)
	return Plots.plot!(real.(ind(d)),real.(val(d));opts...)
end

import Base.length

function length(d::Data1D)
	return length(d.dat)
end

function shift(d::Data1D,δ::Number)
	c=deepcopy(d);
	c.istart += δ
	c.istop += δ
	return c
end

"""
`function hard_shift(d::Data1D,δ::Number)`

shifts the data in `d` by `δ`. The data is re-interpolated as needed. The
resulting `Data1D` object is guaranteed to have the same index as `d`, such that
they are compatible for the purposes of addition etc. Data shifted out of the
window is lost; data points shifted into view are set to zero.
"""
function hard_shift(d::Data1D,δ::Number)
    # compute integer offset and fractional contribution from left
    inc=(d.istop-d.istart)/(length(d.dat)-1)
    n = convert(Integer,div(δ,inc))
    f = mod(δ,inc)/inc
    ndata=zero(d.dat)
    nstart=d.istart
    nstop=d.istop

    if n<0
        ndata[1:(end+n-1)]=(1-f)*d.dat[(1-n):(end-1)]+f*d.dat[(2-n):end]
    else
        ndata[n+1:end-1]=(1-f)*d.dat[1:end-n-1]+f*d.dat[2:end-n]
    end

    return Data1D(ndata,nstart,nstop)
end

import Base.real, Base.imag, Base.abs, Base.angle

function Base.real(s::NMR.Data1D)
    return(NMR.Data1D(real.(s.dat),s.istart,s.istop))
end

# function Base.abs(s::NMR.Data1D)
#     return(NMR.Data1D(abs.(s.dat),s.istart,s.istop))
# end

function Base.imag(s::NMR.Data1D)
    return(NMR.Data1D(imag.(s.dat),s.istart,s.istop))
end

function Base.angle(s::NMR.Data1D)
    return(NMR.Data1D(angle.(s.dat),s.istart,s.istop))
end

#end

max_stars_repo_name: marcel-utz/NMR

max_issues_repo_name: marcel-utz/NMR


text:
module JITTools

import Base: show, start, next, done

export load_descriptor, datadump, otool, dwarfdump, buffer,
       datapointer, datadumpr, dwarfdumpr,
# Reexport from ObjFileBase
       readmeta

immutable jit_code_entry
    next_entry::Ptr{jit_code_entry}
    prev_entry::Ptr{jit_code_entry}
    symfile_addr::Uint64
    symfile_size::Uint64
end

immutable jit_descriptor
   version::Uint32
   action_flag::Uint32
   relevant_entry::Ptr{jit_code_entry}
   first_entry::Ptr{jit_code_entry}
end

load_descriptor() = unsafe_load(cglobal(:__jit_debug_descriptor,jit_descriptor))

show(io::IO, desc::jit_descriptor) =
    print(io,"JIT Descriptor (Version ",desc.version,")")

show(io::IO, entry::jit_code_entry) =
    print(io,"JIT Entry for Symfile at 0x",hex(entry.symfile_addr)," of size ",entry.symfile_size, " bytes")

# Iterating through all the jit entries
start(desc::jit_descriptor) = desc.first_entry
function next(desc::jit_descriptor,state::Ptr{jit_code_entry})
    if state == C_NULL
        throw(BoundsError())
    end
    entry = unsafe_load(state)
    (entry, entry.next_entry)
end
done(desc::jit_descriptor, state::Ptr{jit_code_entry}) = state == C_NULL

datasize(data::jit_code_entry) = int(data.symfile_size)
datapointer(data::jit_code_entry) = convert(Ptr{Uint8},data.symfile_addr)

buffer(data...) = IOBuffer(pointer_to_array(datapointer(data...),datasize(data...),false),true,true)

datadump(io::IO,data...) = write(io,datapointer(data...),datasize(data...))

@osx_only begin

using MachO
import ObjFileBase: readmeta

readmeta(data::jit_code_entry) = readmeta(buffer(data),MachO.MachOHandle)

function datadumpr(io::IO,data::jit_code_entry)
    buf = IOBuffer()
    datadump(buf,data)
    handle = readmeta(data)
    for cmd in LoadCmds(handle)
        if eltype(cmd) <: MachO.segment_commands
            for sec in Sections(cmd)
                seek(buf,sec.offset)
                if int(datapointer(sec)) > 0x10000
                    write(buf,datapointer(sec),datasize(sec))
                end
            end
        end
    end
    write(io,takebuf_array(buf))
end

function otool(data::jit_code_entry,arg=``)
    file = tempname()
    open(file,"w") do f
        datadump(f,data)
    end
    run(`otool $arg $file`)
    Base.FS.unlink(file)
end

function otoolr(data::jit_code_entry,arg=``)
    file = tempname()
    open(file,"w") do f
        datadumpr(f,data)
    end
    run(`otool $arg $file`)
    Base.FS.unlink(file)
end

function dwarfdump(data::jit_code_entry,arg=``)
    file = tempname()
    open(file,"w") do f
        datadump(f,data)
    end
    run(`dwarfdump $arg $file`)
    Base.FS.unlink(file)
end

function dwarfdumpr(data::JITTools.jit_code_entry,arg=``)
    file = tempname()
    open(file,"w") do f
        JITTools.datadumpr(f,data)
    end
    stdout, stdin, p = readandwrite(`dwarfdump $arg $file`)
    close(stdin)
    print(readall(stdout))
    Base.FS.unlink(file)
end

function dwarfdumpr(desc::jit_descriptor, fname)
    for d in desc
        handle = readmeta(d)
        for cmd in LoadCmds(handle)
            if eltype(cmd) <: MachO.symtab_command
                for sym in Symbols(cmd)
                    if contains(bytestring(symname(cmd,sym)),fname)
                        dwarfdumpr(d)
                    end
                end
            end
        end
    end
end

datapointer(data::Union(MachO.section,MachO.section_64)) = convert(Ptr{Uint8},data.addr)
datasize(data::Union(MachO.section,MachO.section_64)) = int(data.size)

# Get unrelocated section
datapointer(entry::jit_code_entry, section::Union(MachO.section,MachO.section_64)) =
    convert(Ptr{Uint8},datapointer(entry)+section.offset)
# This is strictly speaking incorrect as the vm size may be different from the size in the file.
# For now I don't care about that though
datasize(entry::jit_code_entry,data::Union(MachO.section,MachO.section_64)) = datasize(data)

end

end # module

max_stars_repo_name: ambujone/DIDebug

max_issues_repo_name: ambujone/DIDebug

max_forks_repo_name: ambujone/DIDebug


text:
using GR, LinearAlgebra, LowRankApprox, Formatting, ImageFiltering, Images, Plots, ImageMagick, Colors, TestImages, ImageView, ImageTransformations

I=load("rgb.jpg")

G=Gray.(I)
# size(G)

A=imresize(G, ratio=1/5)
# rank(A), size(A), eltype(A)

(p,q)=size(A)
B=Array{Float64}(A)
U, S, V=psvd(B)
B010=sum(S[n]*U[1:p,n]*(V[1:q,n])' for n=1:10);
B030=sum(S[n]*U[1:p,n]*(V[1:q,n])' for n=1:30);
B050=sum(S[n]*U[1:p,n]*(V[1:q,n])' for n=1:50);
B070=sum(S[n]*U[1:p,n]*(V[1:q,n])' for n=1:70);

plot(Gray.(B),
    xaxis=false,
    xticks=false,
    yaxis=false,
    yticks=false,
    grid=false,
    title="B=B(365) 365 by 548 full rank matrix")
savefig("grayscale_svd_365.png")

plot(Gray.(B010),
    xaxis=false,
    xticks=false,
    yaxis=false,
    yticks=false,
    grid=false,
    title="B(10) rank 10 approximation")
savefig("grayscale_svd_010.png")

plot(Gray.(B030),
    xaxis=false,
    xticks=false,
    yaxis=false,
    yticks=false,
    grid=false,
    title="B(30) rank 30 approximation")
savefig("grayscale_svd_030.png")

plot(Gray.(B050),
    xaxis=false,
    xticks=false,
    yaxis=false,
    yticks=false,
    grid=false,
    title="B(50) rank 50 approximation")
savefig("grayscale_svd_050.png")

plot(Gray.(B070),
    xaxis=false,
    xticks=false,
    yaxis=false,
    yticks=false,
    grid=false,
    title="B(70) rank 70 approximation")
savefig("grayscale_svd_070.png")


text:
module SIMDops
import SIMDPirates
using  SIMDPirates: Vec, vload, vsum, vabs, vfma, vifelse, vless, vzero

# * Generic operations on SIMD packs

fma(x...) = Base.fma(x...)
fma(x::T, y::T, z::T) where {T<:NTuple} = SIMDPirates.vfma(x, y, z)

fptype(::Type{Vec{W, T}}) where {W, T} = T

# TODO: This probably belong to SIMDPirates
@inline f64(x::Number)     = Float64(x)
@inline f64(x::VecElement) = VecElement(f64(x.value))
@inline f64(x::Vec)        = broadcast(f64, x)

# * Macros rewriting mathematical operations

# ** Helper functions

replace_simd_ops(x, _, _) = x

function replace_simd_ops(sym::Symbol, ops, updates)
    # Replace e.g:  +
    #           =>  eadd
    for (op, replacement) in ops
        sym === op && return replacement
    end
    sym
end

function replace_simd_ops(expr::Expr, ops, updates)
    # Replace e.g:  lhs += rhs
    #           =>  lhs = eadd(lhs, rhs)
    for (op, replacement) in updates
        if expr.head === op
            lhs = expr.args[1]
            rhs = replace_simd_ops(expr.args[2], ops, updates)
            return :($lhs = $replacement($lhs, $rhs))
        end
    end

    newexpr = Expr(replace_simd_ops(expr.head, ops, updates))
    for arg in expr.args
        push!(newexpr.args, replace_simd_ops(arg, ops, updates))
    end
    newexpr
end

# ** Explicit/exact operations (non-fusible)

eadd(x...) = Base.:+(x...)
eadd(x::T, y::T) where {T<:NTuple} = SIMDPirates.evadd(x, y)

esub(x...) = Base.:-(x...)
esub(x::NTuple) = broadcast(esub, x)
esub(x::VecElement) = VecElement(-x.value)
esub(x::T, y::T) where {T<:NTuple} = SIMDPirates.evsub(x, y)

emul(x...) = Base.:*(x...)
emul(x::T, y::T) where {T<:NTuple} = SIMDPirates.evmul(x, y)

function explicit(expr::Expr)
    replace_simd_ops(expr,
                     (:+   => eadd,
                      :-   => esub,
                      :*   => emul,
                      :fma => fma),
                     (:+=  => eadd,
                      :-=  => esub,
                      :*=  => emul))
end

macro explicit(expr)
    explicit(expr) |> esc
end

# * Fast/fusible operations

fadd(x...) = Base.:+(x...)
fadd(x::T, y::T) where {T<:NTuple} = SIMDPirates.vadd(x, y)

fsub(x...) = Base.:-(x...)
fsub(x::NTuple) = broadcast(esub, x)
fsub(x::T, y::T) where {T<:NTuple} = SIMDPirates.vsub(x, y)

fmul(x...) = Base.:*(x...)
fmul(x::T, y::T) where {T<:NTuple} = SIMDPirates.vmul(x, y)

function fusible(expr::Expr)
    replace_simd_ops(expr,
                     (:+   => fadd,
                      :-   => fsub,
                      :*   => fmul,
                      :fma => fma),
                     (:+=  => fadd,
                      :-=  => fsub,
                      :*=  => fmul))
end

macro fusible(expr)
    fusible(expr) |> esc
end

end


text:
function set_basis(spx::SpxData, basis::Vector{Int})
    # Set basic and nonbasic indices
    spx.basis_status = Vector{Int}(undef, ncols(spx.lpdata))
    spx.basic = basis
    spx.nonbasic = Int[]
    sizehint!(spx.nonbasic, ncols(spx.lpdata) - nrows(spx.lpdata))
    for j = 1:ncols(spx.lpdata)
        if !in(j, spx.basic)
            push!(spx.nonbasic, j)
            axl = abs(spx.lpdata.v_lb[j])
            axu = abs(spx.lpdata.v_ub[j])
            if spx.lpdata.v_lb[j] > -Inf && axl <= axu
                spx.basis_status[j] = BASIS_AT_LOWER
            elseif spx.lpdata.v_ub[j] < Inf && axl > axu
                spx.basis_status[j] = BASIS_AT_UPPER
            else
                spx.basis_status[j] = BASIS_FREE
            end
        else
            spx.basis_status[j] = BASIS_BASIC
        end
    end
    set_nonbasic(spx)
end

function set_basis_status(spx::SpxData, basis_status::Vector{Int})
    # Set basic and nonbasic indices
    spx.basis_status = basis_status
    spx.basic = Int[]
    spx.nonbasic = Int[]
    sizehint!(spx.basic, nrows(spx.lpdata))
    sizehint!(spx.nonbasic, ncols(spx.lpdata) - nrows(spx.lpdata))
    for j = 1:ncols(spx.lpdata)
        if basis_status[j] == BASIS_BASIC
            push!(spx.basic, j)
        else
            push!(spx.nonbasic, j)
        end
    end
    set_nonbasic(spx)
end

function set_basis(spx::SpxData, basis::Vector{Int}, basis_status::Vector{Int})
    # Set basic and nonbasic indices
    spx.basis_status = basis_status
    spx.basic = basis
    spx.nonbasic = Int[]
    sizehint!(spx.nonbasic, ncols(spx.lpdata) - nrows(spx.lpdata))
    for j = 1:ncols(spx.lpdata)
        if !in(j, spx.basic)
            push!(spx.nonbasic, j)
        end
    end
    set_nonbasic(spx)
end

function set_basis(spx::SpxData{T}, x::T) where T
    spx.basis_status = Vector{Int}(undef, ncols(spx.lpdata))
    spx.basic = Int[]
    spx.nonbasic = Int[]
    sizehint!(spx.basic, nrows(spx.lpdata))
    sizehint!(spx.nonbasic, ncols(spx.lpdata) - nrows(spx.lpdata))
    for j = 1:ncols(spx.lpdata)
        if x[j] == spx.lpdata.v_lb[j]
            spx.basis_status[j] = BASIS_AT_LOWER
            push!(spx.nonbasic, j)
        elseif x[j] == spx.lpdata.v_ub[j]
            spx.basis_status[j] = BASIS_AT_UPPER
            push!(spx.nonbasic, j)
        elseif length(spx.basic) < nrows(spx.lpdata)
            spx.basis_status[j] = BASIS_BASIC
            push!(spx.basic, j)
        else
            spx.basis_status[j] = BASIS_FREE
            push!(spx.nonbasic, j)
        end
    end
    spx.x .= x
end

function set_nonbasic(spx::SpxData)
    for j in spx.nonbasic
        if spx.basis_status[j] == BASIS_AT_LOWER
            spx.x[j] = spx.lpdata.v_lb[j]
        elseif spx.basis_status[j] == BASIS_AT_UPPER
            spx.x[j] = spx.lpdata.v_ub[j]
        elseif spx.basis_status[j] == BASIS_FREE
            # @warn("Free nonbasic variable (x$j) is set to zero.")
            spx.x[j] = 0.0
        end
    end
end


text:
include("converters.jl")

"""
A function to run length decode an int array.
in_array: the input array of integers
"""
function runlengthdecode(in_array::Array{T,1}) where {T<:Integer}
    switch=false
    out_array=Int32[]
    this_item = 0
    for item in in_array
        if switch==false
            this_item = item
            switch=true
        else
            switch=false
            for i=1:item
                push!(out_array,Int32(this_item))
            end
        end
    end
    return out_array
end

"""
A function to run length decode an int array.
in_array: the inptut array of integers
"""
function runlengthencode(in_array::Array{T,1}) where {T<:Integer}
    if length(in_array)==0
        return Int32[]
    end
    curr_ans = in_array[1]
    out_array = Int32[curr_ans]
    counter = 1
    for in_int in in_array[2:end]
        if in_int == curr_ans
            counter+=1
        else
            push!(out_array, Int32(counter))
            push!(out_array, Int32(in_int))
            curr_ans = in_int
            counter = 1
        end
    end
    # Add the final counter
    push!(out_array, Int32(counter))
    return out_array
end

"""
A function to delta decode an int array.
in_array: the input array of integers
"""
function deltadecode(in_array::Array{T,1}) where {T<:Integer}
    if length(in_array) == 0
        return Int32[]
    end
    this_ans = in_array[1]
    out_array = Int32[this_ans]
    for i = 2:length(in_array)
        this_ans += in_array[i]
        push!(out_array,Int32(this_ans))
    end
    return out_array
end

"""
A function to delta decode an int array.
in_array: the input array to be delta encoded
"""
function deltaencode(in_array::Array{T,1}) where {T<:Integer}
    if length(in_array)==0
        return Int32[]
    end
    curr_ans = in_array[1]
    out_array = Int32[curr_ans]
    for in_int in in_array[2:end]
        push!(out_array, Int32(in_int-curr_ans))
        curr_ans = in_int
    end
    return out_array
end

function deltarecursivefloatdecode(in_array::Array{UInt8,1}, param)
    return floats(deltadecode(recursiveindexdecode(ints(in_array,2))),param)
end

function deltarecursivefloatencode(in_array::Array{Float32,1}, param)
    return bytes(recursiveindexencode(deltaencode(ints(in_array,param))))
end

function runlengthfloatdecode(in_array::Array{UInt8,1}, param)
    return floats(runlengthdecode(ints(in_array,4)),param)
end

function runlengthfloatencode(in_array::Array{Float32,1}, param)
    return bytes(runlengthencode(ints(in_array,param)))
end

function runlengthdeltaintdecode(in_array::Array{UInt8,1})
    return deltadecode(runlengthdecode(ints(in_array, 4)))
end

function runlengthdeltaintencode(in_array::Array{T,1}) where {T<:Integer}
    return bytes(runlengthencode(deltaencode(in_array)))
end

function runlengthchardecode(in_array::Array{UInt8,1})
    return chars(runlengthdecode(ints(in_array, 4)))
end

function runlengthcharencode(in_array::Array{Char,1})
    return bytes(runlengthencode(ints(in_array)))
end

function stringdecode(in_array::Array{UInt8,1})
    return decodechainlist(in_array)
end

function stringencode(in_array::Array{String,1})
    return encodechainlist(in_array)
end

#1-byte Int decode
function intdecode(in_array::Array{UInt8,1})
    return ints(in_array, 1)
end

#1-byte Int encode
function intencode(in_array::Array{Int8,1})
    return bytes(in_array)
end

#4-byte Int decode
function fourbyteintdecode(in_array::Array{UInt8,1})
    return ints(in_array, 4)
end

#4-byte Int encode
function fourbyteintencode(in_array::Array{Int32,1})
    return bytes(in_array)
end

"""
Parse the header of an input byte array and then decode using the input array,
the codec and the appropirate parameter.
input_array: the array to be decoded
"""
function decodearray(input_array::Array{UInt8,1})
    codec, length, param, input_array = parseheader(input_array)
    if codec == 10
        return deltarecursivefloatdecode(input_array, param)
    elseif codec == 9
        return runlengthfloatdecode(input_array, param)
    elseif codec == 8
        return runlengthdeltaintdecode(input_array)
    elseif codec == 6
        return runlengthchardecode(input_array)
    elseif codec == 5
        return stringdecode(input_array)
    elseif codec == 2
        return intdecode(input_array)
    elseif codec == 4
        return fourbyteintdecode(input_array)
    else
        throw(ArgumentError("Invalid codec while parsing MMTF data!"))
    end
end

"""
Encode the array using the method and then add the header to this array.
input_array: the array to be encoded
codec: the integer index of the codec to use
param: the integer parameter to use in the function
"""
function encodearray(input_array::Array{T,1}, codec, param) where {T<:Union{Real,Char,String}}
    if codec == 10
        return addheader(deltarecursivefloatencode(input_array, param), codec, length(input_array), param)
    elseif codec == 9
        return addheader(runlengthfloatencode(input_array, param), codec, length(input_array), param)
    elseif codec == 8
        return addheader(runlengthdeltaintencode(input_array), codec, length(input_array), param)
    elseif codec == 6
        return addheader(runlengthcharencode(input_array), codec, length(input_array), param)
    elseif codec == 5
        return addheader(stringencode(input_array), codec, length(input_array), param)
    elseif codec == 2
        return addheader(intencode(input_array), codec, length(input_array), param)
    elseif codec == 4
        return addheader(fourbyteintencode(input_array), codec, length(input_array), param)
    else
        throw(ArgumentError("Invalid codec while encoding to MMTF data!"))
    end
end


text:
using LinearAlgebra

⊗ = kron

function encode(s::Tuple{Array{Float64,1},Int64})
  n = sqrt(s[1][1]^2 + s[1][2]^2)
  # return (s[1]/n, s[2])
  return (s[1]/n ⊗ s[1]/n, s[2])
end

function encode(ar::Array{Tuple{Array{Float64,1},Int64}})
  return encode.(ar)
end

max_stars_repo_name: dandanua/QML

max_issues_repo_name: dandanua/QML

max_forks_repo_name: dandanua/QML


text:
using Documenter, MultivariatePolynomials

makedocs(
    sitename = "MultivariatePolynomials",

    format = Documenter.HTML(
        # See https://github.com/JuliaDocs/Documenter.jl/issues/868
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),
    # See https://github.com/jump-dev/JuMP.jl/issues/1576
    strict = true,

    pages = [
        "Introduction" => "index.md",
        "Types" => "types.md",
        "Substitution" => "substitution.md",
        "Differentiation" => "differentiation.md",
        "Division" => "division.md",
    ]
)

deploydocs(
    repo   = "github.com/JuliaAlgebra/MultivariatePolynomials.jl.git",
)


text:
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using Random
using ForwardDiff
Random.seed!(233)

function twolayer(y, x, w1, w2, b1, b2)
    f = x -> begin
        w1 = reshape(w1, 10, 10)
        w2 = reshape(w2, 10, 10)
        z = w2*tanh.(w1*x+b1)+b2
    end
    y[:] = ForwardDiff.jacobian(f, x)[:]
end
two_layer = load_op("build/libTwoLayer", "two_layer")

w1 = rand(100)
w2 = rand(100)
b1 = rand(10)
b2 = rand(10)
x = rand(10)
J = rand(100)
twolayer(J, x, w1, w2, b1, b2)

y = two_layer(constant(x), constant(w1), constant(b1), constant(w2), constant(b2))
sess = Session(); init(sess)
J0 = run(sess, y)
@show norm(J-J0)


text:
module TestCLI
using Base: decompose
import Logging
using Test
using YAML

# include source of bin file here to avoid world age issues
include("../bin/prepare_rqc_simulation_files.jl")

@testset "Test prepare rqc input cli script" begin
    # create empty temporary directory
    mktempdir() do path
        prefix = joinpath(path, "rqc_3_3_8")
        args = ["-p", prefix, "--time", "10", "-n", "15"]
        Logging.with_logger(Logging.NullLogger()) do # suppress logging
            main(args)
        end
        @test all([isfile(prefix * suffix) for suffix in [".qx", ".jld2", ".yml"]])

        params = YAML.load_file(prefix * ".yml")
        @test length(params["output"]["params"]["bitstrings"]) == 15
    end

    # create empty temporary directory
    mktempdir() do path
        prefix = joinpath(path, "rqc_3_3_8")
        N = 20
        args = ["-p", prefix, "-n", "$N", "--time", "30", "--output_method", "Rejection"]
        Logging.with_logger(Logging.NullLogger()) do # suppress logging
            main(args)
        end
        @test all([isfile(prefix * suffix) for suffix in [".qx", ".jld2", ".yml"]])

        params = YAML.load_file(prefix * ".yml")
        @test params["output"]["method"] == "Rejection"
        @test params["output"]["params"]["num_samples"] == N
    end

    # create empty temporary directory
    mktempdir() do path
        prefix = joinpath(path, "rqc_3_3_8")
        N = 20
        args = ["-p", prefix, "-n", "$N", "--time", "30"]
        Logging.with_logger(Logging.NullLogger()) do # suppress logging
            main(args)
        end
        @test all([isfile(prefix * suffix) for suffix in [".qx", ".jld2", ".yml"]])

        params = YAML.load_file(prefix * ".yml")
        @test length(params["output"]["params"]["bitstrings"]) <= N
    end
end
end


text:
"""
Vibration analysis of double-cell box structure

Reference values: 20.280 25.43 28.79 30.62 [Hz]
From: Fan, Luah, Free vibration analysis of arbitrary thin shell structures by using spline finite element, Journal Sound Vib, (179), 763-776, 1995.
"""
module dbcs_vibration_examples

using Arpack
using FinEtools
using FinEtools.MeshExportModule.VTKWrite: vtkwrite
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FEMMShellT3FFModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using VisualStructures: plot_nodes, plot_midline, render, plot_space_box, plot_midsurface, space_aspectratio, save_to_json

function _execute_dsg_model(formul, n, visualize)
    E = 207*phun("GPa")
    nu = 0.3;
    rho= 7850*phun("KG/M^3");
    L0 = 2.54 * phun("m")
    B0 = 1.27 * phun("m")
    thickness = 12.7*phun("mm");

    # Report
    @info "Mesh: $n elements per side"

    tolerance = B0/n/1000
    X = [0.0 0.0 0.0;
        B0 0.0 0.0;
        B0 B0 0.0;
        B0 2*B0 0.0;
        0.0 2*B0 0.0;
        0.0 B0 0.0;
    ]
    l2fens = FENodeSet(X);
    conn = [1 2; 6 3; 5 4; 2 3; 3 4; 5 6; 6 1]
    l2fes = FESetL2(conn);
    nLayers = 5
    ex = function (x, k)
        x[3] += k*L0/nLayers
        x
    end
    fens,fes = Q4extrudeL2(l2fens, l2fes, nLayers, ex); # Mesh
    fens, fes = Q4toT3(fens, fes)
    for r in 1:n
        fens, fes = T3refine(fens, fes)
    end
    @show count(fens)
    vtkwrite("dbcs.vtu", fens, fes)

    mater = MatDeforElastIso(DeforModelRed3D, rho, E, nu, 0.0)

    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)
    femm = formul.make(IntegDomain(fes, TriRule(1), thickness), mater)
    associate = formul.associategeometry!
    stiffness = formul.stiffness
    mass = formul.mass

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    # simple support
    l1 = selectnode(fens; box = Float64[-Inf Inf -Inf Inf 0.0 0.0], inflate = tolerance)
    for i in [1, 2, 3, 4, 5, 6]
        setebc!(dchi, l1, true, i)
    end

    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the system matrix
    associategeometry!(femm, geom0)
    K = stiffness(femm, geom0, u0, Rfield0, dchi);
    M = mass(femm, geom0, dchi);

    # Solve
    OmegaShift = 0.0*2*pi
    neigvs = 4
    d, v, nconv = eigs(K+OmegaShift*M, M; nev=neigvs, which=:SM, explicittransform=:none)
    d[:] = d .- OmegaShift;
    fs = real(sqrt.(complex(d)))/(2*pi)
    @show fs

    # Visualization
    U = v[:, 4]
    scattersysvec!(dchi, (L0/4)/maximum(abs.(U)).*U)
    update_rotation_field!(Rfield0, dchi)
    plots = cat(plot_space_box([[0 0 0]; [L0 L0 L0]]),
        #plot_nodes(fens),
        plot_midsurface(fens, fes; x = geom0.values, u = dchi.values[:, 1:3], R = Rfield0.values);
    dims = 1)
    pl = render(plots)
end

function test_convergence()
    formul = FEMMShellT3FFModule
    @info "FV12 free vibration, formulation=$(formul)"
    for n in [1, 2, 3, 4]
        _execute_dsg_model(formul, n, false)
    end
    return true
end

function allrun()
    println("#####################################################")
    println("# test_convergence ")
    test_convergence()
    return true
end # function allrun

@info "All examples may be executed with "
println("using .$(@__MODULE__); $(@__MODULE__).allrun()")

end # module
nothing


text:
using FastRefine
using LazySets
import FastRefine: forward_network, forward_affine_map, ishull

nnet = read_nnet("nnet/86442.nnet")

delta = 0.3

solver3 = DimGrid(delta)

in_hyper = Hyperrectangle(fill(1.0, 8), fill(1.0, 8))
out_hyper = Hyperrectangle(fill(0.0, 2), fill(10.0, 2))
problem = Problem(nnet, in_hyper, out_hyper)

file = open("results/group3.txt", "a")
print(file, "Test Result of Group 3:\n\n")

#solver3

time3 = 0

for i = 1:1
    timed_result =@timed solve(solver3, problem)
    print(file, "DimGrid - test " * string(i) * " - Time: " * string(timed_result.time) * " s")
    print(file, " - Output: " * string(timed_result.value) * "\n")
    global time3 += timed_result.time
end

print(file, "Average time: " * string(time3/1) * " s\n\n")

close(file)


text:
# Load modules
using Plots, LaTeXStrings
pyplot()
using ThinFilmsTools

function main()
    # Define beam
    λ = LinRange(310,1000,500) # wavelength range [nm]
    θ = [0.] # angle of incidence [degrees]
    beam = PlaneWave(λ,θ)
    # Find wavelength closest to 405.0 nm
    aux1 = Utils.findClosest(vec(beam.λ),405.0)
    # Define layers
    layers = [
        LayerTMMO(RIdb.air(beam.λ)),
        LayerTMMO(RIdb.sno2f(beam.λ); d=150., nλ0=RIdb.sno2f([beam.λ[aux1]])),
        LayerTMMO(RIdb.silicon(beam.λ)),
    ]
    # call main script
    sol = tmm_optics(beam, layers)
    return sol
end

sol = main()

# plot the R, T and A spectra
plot(Spectrum1D(),
     sol.beam.λ, [sol.Spectra.Rp, sol.Spectra.Tp, 1.0.-(sol.Spectra.Rp.+sol.Spectra.Tp)],
     label=["Reflectance" "Transmittance" "Absorbance"],
     line=([:solid :dash :dashdot]),
     ylims=(0.0,1.0),
     xlims=(sol.beam.λ[1], sol.beam.λ[end]),
)
gui()

# plot the refractive index profile
plot(RIprofile(), sol)
gui()


text:
# Copyright 2018 Erik Edin
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

using Test
using Behavior.Selection
using Behavior.Gherkin
using Behavior.Gherkin: AbstractScenario

@testset "Selection            " begin
    # These tests check that the tag selector is parsed and that the expressions are used
    # properly. It uses only tags on scenarios to check this, for simplicity. Tests in a
    # section below check that tags are inherited properly between features and scenarios.
    @testset "Tag selector" begin
        @testset "Feature has tag @foo and no scenarios; Selecting @foo returns no scenarios" begin
            # Arrange
            header = FeatureHeader("Some feature", String[], ["@foo"])
            feature = Feature(header, AbstractScenario[])

            # Act
            selector = parsetagselector("@foo")
            newfeature = select(selector, feature)

            # Assert
            @test isempty(newfeature.scenarios)
        end

        @testset "Feature has tag @foo and one scenario; selecting @foo returns the feature unchanged" begin
            # Arrange
            header = FeatureHeader("Some feature", String[], ["@foo"])
            feature = Feature(header, AbstractScenario[
                Scenario("Some scenario", String[], ScenarioStep[]),
            ])

            # Act
            selector = parsetagselector("@foo")
            newfeature = select(selector, feature)

            # Assert
            @test newfeature.scenarios == feature.scenarios
        end

        @testset "Feature has no tags and no scenarios; select @foo returns no scenarios" begin
            # Arrange
            header = FeatureHeader("Some feature", String[], String[])
            feature = Feature(header, AbstractScenario[
                Scenario("Some scenario", String[], ScenarioStep[]),
            ])

            # Act
            selector = parsetagselector("@foo")
            newfeature = select(selector, feature)

            # Assert
            @test isempty(newfeature.scenarios)
        end

        @testset "Feature has tag @bar; Select for @foo returns nothing" begin
            # Arrange
            header = FeatureHeader("Some feature", String[], ["@bar"])
            feature = Feature(header, AbstractScenario[
                Scenario("Some scenario", String[], ScenarioStep[]),
            ])

            # Act
            selector = parsetagselector("@foo")
            newfeature = select(selector, feature)

            # Assert
            @test isempty(newfeature.scenarios)
        end

        @testset "Feature has tag @foo; Select for @bar returns nothing" begin
            # Arrange
            header = FeatureHeader("Some feature", String[], ["@foo"])
            feature = Feature(header, AbstractScenario[
                Scenario("Some scenario", String[], ScenarioStep[]),
            ])

            # Act
            selector = parsetagselector("@bar")
            newfeature = select(selector, feature)

            # Assert
            @test isempty(newfeature.scenarios)
        end

        @testset "Feature has tags @bar and @foo; Select for @foo returns the feature unchanged" begin
            # Arrange
            header = FeatureHeader("Some feature", String[], ["@bar", "@foo"])
            feature = Feature(header, AbstractScenario[
                Scenario("Some scenario", String[], ScenarioStep[]),
            ])

            # Act
            selector = parsetagselector("@foo")
            newfeature = select(selector, feature)

            # Assert
            @test newfeature.scenarios == feature.scenarios
        end

        @testset "Feature has no tags; Selector is (not @foo); select returns the feature unchanged" begin
            # Arrange
            header = FeatureHeader("Some feature", String[], [])
            feature = Feature(header, AbstractScenario[
                Scenario("Some scenario", String[], ScenarioStep[]),
            ])

            # Act
            selector = parsetagselector("not @foo")
            newfeature = select(selector, feature)

            # Assert
            @test newfeature.scenarios == feature.scenarios
        end

        @testset "Feature has tag @foo; Selector is (not @foo); select returns nothing" begin
            # Arrange
            header = FeatureHeader("Some feature", String[], ["@foo"])
            feature = Feature(header, AbstractScenario[
                Scenario("Some scenario", String[], ScenarioStep[]),
            ])

            # Act
            selector = parsetagselector("not @foo")
            newfeature = select(selector, feature)

            # Assert
            @test isempty(newfeature.scenarios)
        end

        @testset "Feature has tag @foo; Empty selector; Select returns the feature unchanged" begin
            # Arrange
            header = FeatureHeader("Some feature", String[], ["@foo"])
            feature = Feature(header, AbstractScenario[
                Scenario("Some scenario", String[], ScenarioStep[]),
            ])

            # Act
            selector = parsetagselector("")
            newfeature = select(selector, feature)

            # Assert
            @test newfeature.scenarios == feature.scenarios
        end

        @testset "Feature has no tags; Empty selector; Select returns the feature unchanged" begin
            # Arrange
            header = FeatureHeader("Some feature", String[], [])
            feature = Feature(header, AbstractScenario[
                Scenario("Some scenario", String[], ScenarioStep[]),
            ])

            # Act
            selector = parsetagselector("")
            newfeature = select(selector, feature)

            # Assert
            @test newfeature.scenarios == feature.scenarios
        end

        @testset "Feature has no tags; Selector is only whitespace; Select returns the feature unchanged" begin
            # Arrange
            header = FeatureHeader("Some feature", String[], [])
            feature = Feature(header, AbstractScenario[
                Scenario("Some scenario", String[], ScenarioStep[]),
            ])

            # Act
            selector = parsetagselector("   ")
            newfeature = select(selector, feature)

            # Assert
            @test newfeature.scenarios == feature.scenarios
        end

        @testset "Feature has tag @baz; Selector is only whitespace; Select returns the feature unchanged" begin
            # Arrange
            header = FeatureHeader("Some feature", String[], [])
            feature = Feature(header, AbstractScenario[
                Scenario("Some scenario", String[], ScenarioStep[]),
            ])

            # Act
            selector = parsetagselector("   ")
            newfeature = select(selector, feature)

            # Assert
            @test newfeature.scenarios == feature.scenarios
        end
    end

    # These tests check that tags are inherited properly between features and scenarios.
    # It assumes that the parsing of tag selectors works as intended, tested above.
    @testset "Feature and Scenario selection" begin
        @testset "Feature has tag @foo and one scenario, no tags; select returns the feature with one scenario" begin
            # Arrange
            header = FeatureHeader("Some feature", String[], ["@foo"])
            feature = Feature(header, AbstractScenario[
                Scenario("Some scenario", String[], ScenarioStep[])
            ])

            # Act
            selector = parsetagselector("@foo")
            newfeature = select(selector, feature)

            # Assert
            @test newfeature.scenarios == feature.scenarios
        end

        @testset "Only one scenario has tag @bar; select returns the feature with only that scenario" begin
            # Arrange
            header = FeatureHeader("Some feature", String[], [])
            feature = Feature(header, AbstractScenario[
                Scenario("Some scenario", String["@bar"], ScenarioStep[]),
                Scenario("Some other scenario", String[], ScenarioStep[]),
            ])

            # Act
            selector = parsetagselector("@bar")
            newfeature = select(selector, feature)

            # Assert
            @test length(newfeature.scenarios) == 1
            @test newfeature.scenarios[1] == feature.scenarios[1]
        end

        @testset "One scenario has tag @ignore; selecting on (not @ignore) returns only the feature without that tag" begin
            # Arrange
            header = FeatureHeader("Some feature", String[], [])
            feature = Feature(header, AbstractScenario[
                Scenario("Some scenario", String["@ignore"], ScenarioStep[]),
                Scenario("Some other scenario", String["@other"], ScenarioStep[]),
            ])

            # Act
            selector = parsetagselector("not @ignore")
            newfeature = select(selector, feature)

            # Assert
            @test length(newfeature.scenarios) == 1
            @test newfeature.scenarios[1] == feature.scenarios[2]
        end
    end

    @testset "Or expressions" begin
        @testset "Feature has tags @foo and one scenario; selecting @foo,@bar returns the feature unchanged" begin
            # Arrange
            header = FeatureHeader("Some feature", String[], ["@foo"])
            feature = Feature(header, AbstractScenario[
                Scenario("Some scenario", String[], ScenarioStep[]),
            ])

            # Act
            selector = parsetagselector("@foo,@bar")
            newfeature = select(selector, feature)

            # Assert
            @test newfeature.scenarios == feature.scenarios
        end

        @testset "Feature has tags @foo and one scenario; selecting @bar,@foo returns the feature unchanged" begin
            # Arrange
            header = FeatureHeader("Some feature", String[], ["@foo"])
            feature = Feature(header, AbstractScenario[
                Scenario("Some scenario", String[], ScenarioStep[]),
            ])

            # Act
            selector = parsetagselector("@bar,@foo")
            newfeature = select(selector, feature)

            # Assert
            @test newfeature.scenarios == feature.scenarios
        end

        @testset "Feature has tag @foo; Selector is (not @foo,@bar); select returns nothing" begin
            # Arrange
            header = FeatureHeader("Some feature", String[], ["@foo"])
            feature = Feature(header, AbstractScenario[
                Scenario("Some scenario", String[], ScenarioStep[]),
            ])

            # Act
            selector = parsetagselector("not @foo,@bar")
            newfeature = select(selector, feature)

            # Assert
            @test isempty(newfeature.scenarios)
        end

        @testset "Feature has tag @bar; Selector is (not @foo,@bar); select returns nothing" begin
            # Arrange
            header = FeatureHeader("Some feature", String[], ["@bar"])
            feature = Feature(header, AbstractScenario[
                Scenario("Some scenario", String[], ScenarioStep[]),
            ])

            # Act
            selector = parsetagselector("not @foo,@bar")
            newfeature = select(selector, feature)

            # Assert
            @test isempty(newfeature.scenarios)
        end
    end
end