text:
# Define small function that returns true if a value should be replaced, false otherwise
@inline _isreplaceable_(x)::Bool = DataArrays.isna(x)||isnan(x)
@inline _isreplaceable_(x::Nullable)::Bool = isnull(x)||DataArrays.isna(x.value)||isnan(x.value)

###########################
# FunctionCell Interface  #
###########################
"""
	filterg([f],[g],opts)

Generic data filtering, returns an untrained filter. When data is piped into it, a train filter is generated.
`f` is function to obtain targets by calling `Learn.Base(data,f)`, defaults to `identity`. `g` is a function
that must return a `::Bool` and specifies if `true` that a given variable value of an observation is
to be filtered (e.g. if `g=isnan`, `NaN` values of the variable will be replaced with the results of the
filtering the variable). `opts` specifies the filtering:

* `"mean"` replaces values by the mean of the variable
* `"c-mean"` replaces values by the class-wise mean of the variable
* `"median"` replaces values by the median of the variable
* `"c-median"` replaces values by the median of the class
* `"majority"` replaces values by the most frequent value of the variable
* `"c-majority"` replaces values by the most frequent value of the variable in the same class
* `::Function` a 2 argument function which is applied on each variable value and its associated label
* `::Dict` a dictionary indicating the variable (key) and one of the above options (values)

# Examples
```
julia> a = datacell([1 NaN 1 2 2; NaN 0 1 2 NaN; NaN NaN 0 0 1])
DataCell, 5 obs, 3 vars, 0 target(s)/obs

julia> wu = filterg(isnan, Dict(1=>"mean", 2=>(a,b)->123, 3=>"median"))
Generic filter (mixed), varying I/O dimensions, untrained

julia> wt = a |> wu
Generic filter (mixed), 3->3, trained

julia> +wt
Model
`- data: Dict{Int64,Tuple{Function,Any}}
`- properties: Int/Int/MLLabelUtils.LabelEnc.NativeLabels{Any,0}/Dict{String,Any}

julia> +a
3×5 Array{Float64,2}:
   1.0  NaN    1.0  2.0    2.0
 NaN      0.0  1.0  2.0  NaN
 NaN    NaN    0.0  0.0    1.0

julia> +(a |> wt)
3×5 Array{Float64,2}:
   1.0  1.5  1.0  2.0    2.0
 123.0  0.0  1.0  2.0  123.0
   0.0  0.0  0.0  0.0    1.0

julia> a = datacell([1 2 3 4 5; 6 7 8 9 0.0],["a", "a", "a", "b", "b"])
DataCell, 5 obs, 2 vars, 1 target(s)/obs, 2 distinct values: "b"(2),"a"(3)

julia> wu = filterg(x->mod(x,5)>2, "mean")
Generic filter (mean), varying I/O dimensions, untrained

julia> wt = a |> wu
Generic filter (mean), 2->2, trained

julia> +a
2×5 Array{Float64,2}:
 1.0  2.0  3.0  4.0  5.0
 6.0  7.0  8.0  9.0  0.0

julia> +(a |> wt)
2×5 Array{Float64,2}:
 1.0  2.0  2.66667  2.66667  5.0
 6.0  7.0  4.33333  4.33333  0.0
```
"""
filterg(f::Function, g::Function, opts::T where T<:AbstractString) =
	FunctionCell(filterg, (f,g,opts), ModelProperties(), "Generic filter ("*opts*")")

filterg(f::Function, g::Function, opts::T where T<:Function) =
	FunctionCell(filterg, (f,g,opts), ModelProperties(), "Generic filter (function)")

filterg(f::Function, g::Function, opts::T where T<:Dict) =
	FunctionCell(filterg, (f,g,opts), ModelProperties(), "Generic filter (mixed)")

filterg(g::Function, opts) = filterg(identity, g, opts)

filterg(opts) = filterg(identity, _isreplaceable_, opts)

############################
# DataCell/Array Interface #
############################
"""
	filterg(data,[f],[g],opts)

Generates a trained filter.
"""
# Training

filterg(x::T where T<:CellData, opts) = filterg(strip(x), identity, _isreplaceable_, opts)
filterg(x::T where T<:CellData, g::Function, opts) = filterg(strip(x), identity, g, opts)
filterg(x::T where T<:CellData, f::Function, g::Function, opts) = filterg(strip(x), f, g, opts)

filterg(x::Tuple{T,S} where T<:AbstractArray where S<:AbstractArray, opts) = filterg(x[1], x[2], identity, _isreplaceable_, opts)
filterg(x::Tuple{T,S} where T<:AbstractArray where S<:AbstractArray, g::Function, opts) = filterg(x[1], x[2], identity, g, opts)
filterg(x::Tuple{T,S} where T<:AbstractArray where S<:AbstractArray, f::Function, g::Function, opts) = filterg(x[1], x[2], f, g, opts)

filterg(x::T where T<:AbstractArray, opts) = filterg(x, Void[], identity, _isreplaceable_, opts)
filterg(x::T where T<:AbstractArray, g::Function, opts) = filterg(x, Void[], identity, g, opts)
filterg(x::T where T<:AbstractArray, f::Function, g::Function, opts) = filterg(x, Void[], f, g, opts)

filterg(x::T where T<:AbstractArray, y::S where S<:AbstractArray, f::Function, g::Function, opts) = begin

	# Get dictionary or construct a proper one from original;
	# Inputs: filtering options and total number of variables
	_sdict_(x::T, m) where T<:String = Dict(i=>x for i in 1:m)
	_sdict_(x::T, m) where T<:Function = Dict(i=>x for i in 1:m)
	_sdict_(x::T, m) where T<:Dict = begin
		out = Dict{Int,Any}()
		for (k,v) in x
			if(k isa Int)
				@assert (k>=1)&&(k<= m) "[filterg] Index $(k) for a variable out of bounds."
				push!(out, k=>v)
			elseif (k isa UnitRange{Int} || k isa Vector{Int})
				@assert (minimum(k)>=1)&&(maximum(k)<= m) "[filterg] Index $(k) for a variable out of bounds."
				push!(out, (ki=>v for ki in k)...)
			else
				error("[filterg] Unsupported variable index syntax.")
			end
		end
		return out
	end
	dopts = _sdict_(opts, nvars(x))

	# Determine filterg name (for information only)
	_filtername_(x) = begin
		if isempty(x) return "Generic filter ()"
		elseif (length(unique(x)) == 1)
			if (collect(x)[1] isa Function) return "Generic filter (function)"
			else return "Generic filter ("*collect(x)[1]*")"
			end
		else return "Generic filter (mixed)"
		end
	end
	filtername = _filtername_(values(dopts))

	# Get targets (nothing for unlabeled data cells and arrays)
	_targets_(f,y) = targets(f,y)
	_targets_(f,::AbstractArray{Void}) = nothing
	labels = _targets_(f,y)
	@assert !any(isna.(labels))
	@assert !any(isnan.(labels))
	@assert (labels isa Vector)||(labels isa Void) "[filterg] Labels need to be a Vector or ::Void."

	# Get priors
	priors = countmapn(labels)

	# Create model properties

	modelprops = ModelProperties(nvars(x), nvars(x), nothing,
			Dict("priors" => priors,				# Priors
			"variables_to_process" => collect(keys(dopts)),		# Which columns to process (e.g. can be any sort of vector
			"targets_function" => f					# Function used to obtain targets (needed at runtime)
			)
	)

	# Define filtering functions (v - train data vector, t - train targets array, x - runtime data, y - runtime label)
	_mean_(v::T, ::Void) where T<:AbstractVector = begin m = mean(v); return (x,y)->m; end
	_mean_(v::T, t::S) where T<:AbstractVector where S<:AbstractArray = _mean_(v,nothing)

	_median_(v::T, ::Void) where T<:AbstractVector = begin m = median(v); (x,y)->m; end
	_median_(v::T, t::S) where T<:AbstractVector where S<:AbstractArray = _median_(v,nothing)

	_majority_(v::T, ::Void) where T<:AbstractVector = begin
			u = unique(v)
			cm = countapp(v,u)
			_, idx = findmax(cm)
			m = u[idx]
			return (x,y)->m
		end
	_majority_(v::T, t::S) where T<:AbstractVector where S<:AbstractArray = _majority_(v,nothing)

	_cmean_(v::T, ::Void) where T<:AbstractVector = _mean_(v,nothing)
	_cmean_(v::T, t::S) where T<:AbstractVector where S<:AbstractVector = begin
		d = Dict(ti=>mean(v[t.==ti]) for ti in unique(t))
		(x,y) -> d[y]
	end

	_cmedian_(v::T, ::Void) where T<:AbstractVector = _median_(v,nothing)
	_cmedian_(v::T, t::S) where T<:AbstractVector where S<:AbstractVector = begin
		d = Dict(ti=>median(v[t.==ti]) for ti in unique(t))
		(x,y) -> d[y]
	end

	_cmajority_(v::T, ::Void) where T<:AbstractVector = _majority_(v, nothing)
	_cmajority_(v::T, t::S) where T<:AbstractVector where S<:AbstractArray = begin
		#(x,y) -> Dict(k=>_majority_(v[t.==k], nothing) for k in unique(t))[y]
		d = Dict{eltype(t), eltype(v)}()
		@inbounds for ti in unique(t)
			vt = view(v,t.==ti)
			u = unique(vt)
			cm = countapp(vt,u)
			_, idx = findmax(cm)
			m = u[idx]
			push!(d, ti=>m)
		end
		return (x,y)->d[y]
	end

	# Create dictionary for each processed variable Dict(column=>(isreplaceable, filter function))
	modeldata = Dict{Int, Tuple{Function,Function}}()
	for (idx, v) in dopts
		fmask = .!g.(_variable_(x,idx))					# Mask for samples that are not to be replaced
		fvariable = _variable_(x, idx)[fmask]				# Variable values for the mask
		flabels = labels[fmask]						# Labels for the mask

		if v=="mean"
			w = (g, _mean_(fvariable, flabels))
		elseif v=="c-mean"
			w = (g, _cmean_(fvariable, flabels))
		elseif v=="median"
			w = (g, _median_(fvariable, flabels))
		elseif v=="c-median"
			w = (g, _cmedian_(fvariable, flabels))
		elseif v=="majority"
			w = (g, _majority_(fvariable, flabels))
		elseif v=="c-majority"
			w = (g, _cmajority_(fvariable, flabels))
		elseif v isa Function
			w = (g, v)
		else # do not filter if option not recognized
			warn("[filterg] Option $(v) not recognized, will not process variable $(idx)")
			w = (g, (x,y)->x)
		end
		push!(modeldata, idx => w )
	end

	# Return model cell
	return FunctionCell(filterg, Model(modeldata, modelprops), filtername)
end

# Execution
filterg(x::T where T<:CellData, model::Model) = datacell(filterg(strip(x), model))

filterg(x::Tuple{T,S} where T<:AbstractArray where S<:AbstractArray, model::Model) = (filterg(x[1], x[2], model),x[2])

filterg(x::T where T<:AbstractArray, y::S where S<:AbstractArray, model::Model) = begin

	f = model.properties.other["targets_function"]
	xc = deepcopy(x)				# Create local copy of input
	labels = targets(f, y)				# Get labels

	# Iterate over variable index, replaceable and filter functions
	@inbounds for (idx, (isreplaceable, ff)) in model.data
		variable = _variable_(xc,idx) 		# Assign temp variable
		rmask = isreplaceable.(variable)

		# Process variable vector
		variable[rmask] = ff.(variable[rmask], labels[rmask])
	end
	return xc
end

filterg(x::T where T<:AbstractArray, model::Model) = begin

	f = model.properties.other["targets_function"]
	xc = deepcopy(x)				# Create local copy of input

	# Iterate over variable index, replaceable and filter functions
	@inbounds for (idx, (isreplaceable, ff)) in model.data
		variable = _variable_(xc,idx) 		# Assign temp variable
		rmask = isreplaceable.(variable)

		# Process variable vector
		variable[rmask] = ff.(variable[rmask], nothing)
	end
	return xc
end


text:
struct ZZIdl <: NumFieldOrdIdl
  gen::fmpz

  function ZZIdl(x::fmpz)
    if x < 0
      return new(abs(x))
    else
      return new(x)
    end
  end
end

struct ZZFracIdl <: NumFieldOrdFracIdl
  gen::fmpq

  function ZZFracIdl(x::fmpq)
    if x < 0
      return new(abs(x))
    else
      return new(x)
    end
  end
end

order(::ZZIdl) = FlintZZ

order(::ZZFracIdl) = FlintZZ

# constructors
*(::FlintIntegerRing, x::IntegerUnion) = ideal(ZZ, x)

*(x::IntegerUnion, ::FlintIntegerRing) = ideal(ZZ, x)

ideal(::FlintIntegerRing, x::fmpz) = ZZIdl(x)

ideal(::FlintIntegerRing, x::Integer) = ZZIdl(fmpz(x))

ideal(::FlintIntegerRing, x::AbstractVector{fmpz}) = ZZIdl(gcd(x))

ideal(::FlintIntegerRing, x::AbstractVector{<:Integer}) = ZZIdl(fmpz(gcd(x)))

fractional_ideal(::FlintIntegerRing, x::fmpq) = ZZFracIdl(x)

fractional_ideal(::FlintIntegerRing, x::RingElement) = ZZFracIdl(fmpq(x))

*(x::fmpq, ::FlintIntegerRing) = ZZFracIdl(x)

*(::FlintIntegerRing, x::fmpq) = ZZFracIdl(x)

#

norm(x::ZZIdl) = gen(x)

norm(x::ZZFracIdl) = gen(x)

minimum(x::ZZIdl) = gen(x)

minimum(x::ZZFracIdl) = gen(x)

prime_decomposition(O::NfOrd, p::ZZIdl) = prime_decomposition(O, gen(p))

uniformizer(x::ZZIdl) = gen(x)

# printing
function Base.show(io::IO, x::ZZIdl)
  print(io, "$(x.gen)ZZ")
end

function Base.show(io::IO, x::ZZFracIdl)
  print(io, "Fractional ideal $(x.gen)ZZ")
end

# comparison
function ==(I::ZZIdl, J::ZZIdl)
  return I.gen == J.gen
end

# access
gen(I::ZZIdl) = I.gen

gens(I::ZZIdl) = fmpz[I.gen]

gen(I::ZZFracIdl) = I.gen
gens(I::ZZFracIdl) = fmpq[I.gen]

#TODO

# arithmetic

function +(I::ZZIdl, J::ZZIdl)
  g = gcd(I.gen, J.gen)
  return ZZIdl(g)
end

function *(s::fmpz, J::ZZIdl)
  return ZZIdl(s*J.gen)
end

function (J::ZZIdl, s::fmpz)
  return ZZIdl(s*J.gen)
end

# TODO

gcd(I::ZZIdl, J::ZZIdl) = ZZIdl(gcd(I.gen, J.gen))
gcd(I::ZZIdl, n::T) where T <: Union{fmpz, Int} = ZZIdl(gcd(I.gen, n))
gcd(n::T, I::ZZIdl) where T <: Union{fmpz, Int} = ZZIdl(gcd(I.gen, n))

isone(I::ZZIdl) = isone(I.gen)

maximal_order(::FlintRationalField) = ZZ

ideal_type(::FlintIntegerRing) = ZZIdl
order_type(::FlintRationalField) = FlintIntegerRing
ideal_type(FlintIntegerRing) = ZZIdl
order_type(FlintRationalField) = FlintIntegerRing
place_type(::FlintRationalField) = PosInf
place_type(::Type{FlintRationalField}) = PosInf

fractional_ideal_type(::FlintRationalField) = ZZFracIdl

elem_in_nf(x::fmpz) = FlintQQ(x)

nf(::FlintIntegerRing) = FlintQQ
number_field(::FlintIntegerRing) = FlintQQ

# Infinite places

isreal(::PosInf) = true
iscomplex(::PosInf) = false

infinite_places(::FlintRationalField) = [inf]
infinite_place(::FlintRationalField) = inf

function infinite_place(::FlintRationalField, i::Int)
  i !=1 && error("Index must be 1")
  return inf
end

real_places(::FlintRationalField) = PosInf[inf]

complex_places(::FlintRationalField) = PosInf[]

function sign(x::Union{fmpq,fmpz},p::PosInf)
  return sign(x)
end

function signs(a::Union{fmpq,fmpz}, l::Vector{PosInf})
  return Dict((inf, sign(a)))
end

function ispositive(x::Union{fmpq,fmpz},p::PosInf)
  return x > 0
end

function istotally_positive(x::Union{fmpq,fmpz},p::PosInf)
  return x > 0
end

function isnegative(x::Union{fmpq,fmpz},p::PosInf)
  return x < 0
end

number_field(::PosInf) = QQ

uniformizer(::PosInf) = QQ(-1)

infinite_places_uniformizers(::FlintRationalField) = fmpq[QQ(1)]

function hilbert_symbol(a,b, p::ZZIdl)
  return hilbert_symbol(a,b, gen(p))
end

islocal_norm(K, x, p::ZZIdl) = islocal_norm(K, x, gen(p))

function quadratic_defect(q::fmpq, p::ZZIdl)
  return quadratic_defect(q, gen(p))
end

################################################################################
#
#  Support
#
################################################################################

function support(a::fmpq, R::FlintIntegerRing)
  return ZZIdl[p*R for (p, _) in factor(a, R)]
end

################################################################################
#
#  CRT
#
################################################################################

function crt(a::Vector{fmpz}, b::Vector{ZZIdl})
  return crt(a, fmpz[gen(x) for x in b])
end

ext: jl


text:
using MyExample
using Test

@testset "MyExample.jl" begin
    # Write your own tests here.
    #2x + 3y
    @test my_f(2, 1) ==  7
    @test my_f(2, 3) ==  13
    @test my_f(2, 3) == 11
end


text:
"""
Dynamic path planning problem of hovercraft trying to minimize thrust usage
to hit all the waypoints at the alotted target times.
"""

using InfiniteOpt, JuMP, Ipopt, Plots

# Set problem information
max_time = 60
time_points = Vector(0:max_time)
xw = [1 4 6 1; 1 3 0 1] # waypoints
times = [0; 25; 50; 60]

# Initialize the model
m = InfiniteModel(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))

# Set the parameters and variables
@infinite_parameter(m, t in [0, max_time])
@infinite_variable(m, x[1:2](t), start = 1) # position
@infinite_variable(m, v[1:2](t), start = 0) # velocity
@infinite_variable(m, u[1:2](t), start = 0) # thruster input

# Set the initial conditions
@BDconstraint(m, initial_velocity[i = 1:2](t == 0), v[i] == 0)

# Manually implement euler scheme for motion equations
@point_variable(m, x[i](time_points[j]), xp[i = 1:2, j = 1:length(time_points)])
@point_variable(m, v[i](time_points[j]), vp[i = 1:2, j = 1:length(time_points)])
@point_variable(m, u[i](time_points[j]), up[i = 1:2, j = 1:length(time_points)])
@constraint(m, [i = 1:2, j = 1:length(time_points)-1], xp[i, j+1] == xp[i, j] + vp[i, j])
@constraint(m, [i = 1:2, j = 1:length(time_points)-1], vp[i, j+1] == vp[i, j] + up[i, j])

# Hit all the waypoints
@BDconstraint(m, [i = 1:2, j = 1:length(times)](t == times[j]), x[i] == xw[i, j])

# Specify the objective
@objective(m, Min, support_sum(u[1]^2 + u[2]^2, t))

# Optimize the model
optimize!(m)

# Get the results
if has_values(m)
    x_opt = value.(x)
    v_opt = value.(v)
    u_opt = value.(u)
    opt_supports = supports.(x)
    obj_opt = objective_value(m)
end

# Plot the results
scatter(xw[1,:], xw[2,:], label = "Waypoints")
plot!(x_opt[1,:], x_opt[2,:], label = "Trajectory")
xlabel!("x1")
ylabel!("x2")


text:
"""
    ZN{N} <: DAS
ZN symmetry singleton type where N is an Int.
"""
struct ZN{N} <: DAS end
Base.show(io::IO, ::ZN{N}) where N = print(io,"Z",N)

"""
    ZNCharges{N} <: DASCharges
singleton-type for collections of `ZNCharge{N}` objects. Only the type parameter
is provided - `ZNCharges{N}` doesn't have a field.
Upon iteration yields `ZNCharge{N}(i)` for `i = 0,...,N-1`

# Example
```julia-repl
julia>a = ZNCharges{2}()
julia>foreach(println, a)
ZNCharge{2}(0)
ZNCharge{2}(1)
```
"""
struct ZNCharges{N} <: DASCharges end
Base.show(io::IO, s::ZNCharges{N}) where N = print(io,"Z",N,"Charges(",0:(N-1),")")

⊕(::ZNCharges{N}, ::ZNCharges{N}) where N = ZNCharges{N}()
Base.inv(::ZNCharges{N}) where N = ZNCharges{N}()
Base.length(::ZNCharges{N}) where N = N
Base.getindex(s::ZNCharges{N}, i::Int) where N = ZNCharge{N}(i-1)
chargestype(::ZN{N}) where N = ZNCharges{N}
Base.eltype(::Union{ZNCharges{N},Type{ZNCharges{N}}}) where N = ZNCharge{N}
Base.intersect(::ZNCharges{N}, ::ZNCharges{N}) where N = ZNCharges{N}()
Base.intersect(::ZNCharges, ::ZNCharges) = ArgumentError("cannot intersect different Z charges")

"""
    ZNCharge{N} <: DASCharge
holds the charge of a ZN symmetry as an integer.
The integer is taken `mod N` s.t. the charge is always between `0` and `N-1`

# Example
```julia-repl
julia>a = ZNCharge{2}(1);
julia>a ⊕ a
ZNCharge{2}(0)

julia>a = ZNCharge{2}(3) == ZNCharge{2}(1)
true
```
"""
struct ZNCharge{N} <: DASCharge
    ch::Int
    ZNCharge{N}(a) where {N} = new(mod(a,N))
end
Base.show(io::IO, s::ZNCharge{N}) where N = print(io,"Z",N,"Charge(",s.ch,")")

⊕(a::ZNCharge{N}, b::ZNCharge{N}) where N = ZNCharge{N}(a.ch + b.ch)
function chargeindex(ch::ZNCharge{N}, chs::ZNCharges{N}) where N
    ch in chs || throw(ArgumentError(string(ch, " not in ", chs)))
    (ch.ch + 1)
end
chargetype(::ZN{N}) where N = ZNCharge{N}
Base.zero(::Type{ZNCharge{N}}) where N = ZNCharge{N}(zero(Int))
Base.:+(chs::ZNCharges{N}, ch::ZNCharge{N}) where N = ZNCharges{N}()

const Z2 = ZN{2}
const Z2Charges = ZNCharges{2}
const Z2Charge = ZNCharge{2}


text:
# Unit tests for integrations/dftk_integration.jl

@testset "integrations/dftk_integration.jl" begin
    pspkey = list_psp(:Ar, functional = "lda")[1].identifier
    particles = [Atom(:Ar, [i & 1 == 0 ? 7 : 21, i & 2 == 0 ? 7 : 21, i & 4 == 0 ? 7 : 21]u"bohr"; pseudopotential = pspkey) for i ∈ 0:7]
    box = [[28.0, 0, 0], [0, 28.0, 0], [0, 0, 28.0]]u"bohr"
    boundary_conditions = [Periodic(), Periodic(), Periodic()]
    system1 = FlexibleSystem(particles, box, boundary_conditions)

    potential1 = DFTKPotential(5, [1, 1, 1]; damping = 0.7, tol = 1e-4)
    potential2 = DFTKPotential(136.05693123044495u"eV", [1, 1, 1]; damping = 0.7, tol = 1e-4)

    @test potential1.Ecut ≈ potential2.Ecut

    @test Atomistic.calculate_scf(system1, potential1) == potential1.previous_scfres[]
    @test InteratomicPotentials.potential_energy(system1, potential1) isa Float64
    @test InteratomicPotentials.force(system1, potential1) isa AbstractVector{<:StaticVector{3,<:Float64}}

    system2 = DynamicSystem(system1, 0u"s")
    system3 = DynamicSystem(system1, 1u"ħ_au / hartree")

    @test InteratomicPotentials.potential_energy(system2, potential1) == potential1.potential_energy_cache[0]
    @test InteratomicPotentials.force(system3, potential1) isa AbstractVector{<:StaticVector{3,<:Float64}}
    @test keys(potential1.potential_energy_cache) == Set([0.0, 1.0])
    pe2 = potential1.potential_energy_cache[0]
    pe3 = potential1.potential_energy_cache[1]

    # workaround way to validate that cached results are used
    modified_particles = [Atom(:Ar, [i & 1 == 0 ? 5 : 15, i & 2 == 0 ? 5 : 15, i & 4 == 0 ? 5 : 15]u"bohr"; pseudopotential = pspkey) for i ∈ 0:7]
    system4 = FlexibleSystem(modified_particles, box, boundary_conditions)
    system5 = DynamicSystem(system4, 0u"s")
    system6 = DynamicSystem(system4, 1u"ħ_au / hartree")

    @test InteratomicPotentials.potential_energy(system5, potential1) == pe2
    @test InteratomicPotentials.potential_energy(system6, potential1) == pe3
    @test keys(potential1.potential_energy_cache) == Set([0.0, 1.0])
end


text:
@checked struct PassRegistry
    ref::API.LLVMPassRegistryRef
end

Base.unsafe_convert(::Type{API.LLVMPassRegistryRef}, pr::PassRegistry) = pr.ref

export GlobalPassRegistry

GlobalPassRegistry() = PassRegistry(API.LLVMGetGlobalPassRegistry())


text:
using FortranFiles
include("../setting/constants.jl")
include("../types/types.jl")

# * mesh
# nspec for my simulation setting is 4480
function sem_mesh_read(basedir::String, iproc::Int64)
    # * init
    mesh_data = sem_mesh_data()
    f = sem_io_open_file_for_read(basedir, iproc, "solver_data")

    nspec = Int64(read(f, Int32))
    nglob = Int64(read(f, Int32))

    sem_mesh_init!(mesh_data, nspec, nglob)
    mesh_data.nspec = Int64(nspec)
    mesh_data.nglob = Int64(nglob)

    dummy = zeros(Float32, nglob)
    dummy4 = zeros(Float32, NGLLX, NGLLY, NGLLZ, nspec)

    read(f, dummy)
    mesh_data.xyz_glob[1,:] .= Float64.(dummy)
    read(f, dummy)
    mesh_data.xyz_glob[2,:] .= Float64.(dummy)
    read(f, dummy)
    mesh_data.xyz_glob[3,:] .= Float64.(dummy)

    dummy_int1 = similar(mesh_data.ibool, Int32)
    read(f, dummy_int1)
    mesh_data.ibool = Int64.(dummy_int1)

    dummy_int2 = similar(mesh_data.idoubling, Int32)
    read(f, dummy_int2)
    mesh_data.idoubling = Int64.(dummy_int2)

    read(f, mesh_data.ispec_is_tiso)

    read(f, dummy4)
    mesh_data.xix .= Float64.(dummy4)
    read(f, dummy4)
    mesh_data.xiy .= Float64.(dummy4)
    read(f, dummy4)
    mesh_data.xiz .= Float64.(dummy4)
    read(f, dummy4)
    mesh_data.etax .= Float64.(dummy4)
    read(f, dummy4)
    mesh_data.etay .= Float64.(dummy4)
    read(f, dummy4)
    mesh_data.etaz .= Float64.(dummy4)
    read(f, dummy4)
    mesh_data.gammax .= Float64.(dummy4)
    read(f, dummy4)
    mesh_data.gammay .= Float64.(dummy4)
    read(f, dummy4)
    mesh_data.gammaz .= Float64.(dummy4)

    close(f)

    # separate crustal mesh layers for REGIONAL_MOHO_MESH
    # 3-layer crust: 10(third layer), 11, 12(shallowest layer)
    if REGIONAL_MOHO_MESH
        num = 0
        for ispec = 1:nspec
            if mesh_data.idoubling[ispec] == IFLAG_CRUST
                id = num - fld(num, 3) * 3
                mesh_data.idoubling[ispec] = 10 * IFLAG_CRUST + id
                num = num + 1
            end
        end
    end

    # separate mesh layers across 410-km
    # 40: above 410, 41: below 410
    for ispec = 1:nspec
        if mesh_data.idoubling[ispec] == IFLAG_670_220
            iglob = mesh_data.ibool[MIDX, MIDY, MIDZ, ispec]
            # element center coordinate
            xyz_center = mesh_data.xyz_glob[:,iglob]
            depth = (1.0 - sqrt(sum(mesh_data.xyz_glob[:,iglob].^2))) * R_EARTH_KM
            # this is dangerous due to 410 undulation
            # depth < 410 ? mesh_data.idoubling(ispec) = 10 * IFLAG_670_220 : mesh_data.idoubling(ispec) = 10 * IFLAG_670_220 + 1
            if depth < 410
                mesh_data.idoubling[ispec] = 10 * IFLAG_670_220
            else
                mesh_data.idoubling[ispec] = 10 * IFLAG_670_220 + 1
            end
        end
    end
    return mesh_data
end

function sem_io_open_file_for_read(basedir::String, iproc::Int64, tag::String)
    filename = "$(basedir)/proc$(lpad(string(iproc), 6, '0'))_reg1_$(tag).bin"
    f = FortranFile(filename)
    return f
end

function sem_io_open_file_for_write(basedir::String, iproc::Int64, tag::String)
    filename = "$(basedir)/proc$(lpad(string(iproc), 6, '0'))_reg1_$(tag).bin"
    f = FortranFile(filename, "w")
    return f
end

function sem_mesh_init!(mesh_data::sem_mesh_data, nspec::Int64, nglob::Int64)
    mesh_data.xyz_glob = zeros(Float64, 3, nglob)
    mesh_data.ibool = zeros(Int64, NGLLX, NGLLY, NGLLZ, nspec)
    mesh_data.idoubling = zeros(Int64, nspec)
    mesh_data.ispec_is_tiso = zeros(Bool, nspec)
    mesh_data.xix = zeros(Float64, NGLLX, NGLLY, NGLLZ, nspec)
    mesh_data.xiy = zeros(Float64, NGLLX, NGLLY, NGLLZ, nspec)
    mesh_data.xiz = zeros(Float64, NGLLX, NGLLY, NGLLZ, nspec)
    mesh_data.etax = zeros(Float64, NGLLX, NGLLY, NGLLZ, nspec)
    mesh_data.etay = zeros(Float64, NGLLX, NGLLY, NGLLZ, nspec)
    mesh_data.etaz = zeros(Float64, NGLLX, NGLLY, NGLLZ, nspec)
    mesh_data.gammax = zeros(Float64, NGLLX, NGLLY, NGLLZ, nspec)
    mesh_data.gammay = zeros(Float64, NGLLX, NGLLY, NGLLZ, nspec)
    mesh_data.gammaz = zeros(Float64, NGLLX, NGLLY, NGLLZ, nspec)
end

# * data
function sem_io_read_gll_file_1!(basedir::String, iproc::Int64, model_name::String, model_gll::Array{Float64,4})
    nspec = size(model_gll)[4]
    f = sem_io_open_file_for_read(basedir, iproc, model_name)

    dummy = zeros(Float32, NGLLX, NGLLY, NGLLZ, nspec)
    read(f, dummy)
    close(f)
    model_gll .= Float64.(dummy)
end

function sem_io_read_gll_file_n!(basedir::String, iproc::Int64, model_names::Vector{String}, nmodel::Int64, model_gll::Array{Float64,5})
    dummy = similar(model_gll[1,:,:,:,:])
    for imodel = 1:nmodel
        sem_io_read_gll_file_1!(basedir, iproc, model_names[imodel], dummy)
        model_gll[imodel,:,:,:,:] = dummy
    end
end

function sem_io_write_gll_file_n(basedir::String, iproc::Int64, model_names::Vector{String}, nmodel::Int64, model_gll::Array{Float64,5})
    for imodel = 1:nmodel
        sem_io_write_gll_file_1(basedir, iproc, model_names[imodel], model_gll[imodel,:,:,:,:])
    end
end

function sem_io_write_gll_file_1(basedir::String, iproc::Int64, model_name::String, model_gll::Array{Float64,4})
    f = sem_io_open_file_for_write(basedir, iproc, model_name)
    dummy = Float32.(model_gll)
    write(f, dummy)
    close(f)
end

max_stars_repo_name: ziyixi/SeisScripts

max_issues_repo_name: ziyixi/SeisScripts

max_forks_repo_name: ziyixi/SeisScripts


text:
module GCP
include("util.jl")
include("credential.jl")
include("session.jl")
include("bigquery.jl")
using .Bigquery
mutable struct Project
    id::String
    credential::Credential
    session::Session
    authorize
    bigquery::Bigquery.Service
    function Project(credential, scopes::AbstractVector{<:AbstractString})
        p = new()
        begin
            if credential isa String && isfile(credential)
                p.credential = Credential(credential)
            else
                p.credential = credential
            end
            p.session = Session(p.credential, scopes)
            p.id = (p.credential).project_id
            p.authorize = function authorize(; kws...)
                    auth = GCP.authorize(p.session; kws...)
                    headers = Dict{String, String}("Authorization" => "$(auth[:token_type]) $(auth[:access_token])")
                    headers
                end
        end
        p.bigquery = Bigquery.Service(p)
        p
    end
end
end

max_stars_repo_name: ProfLeao/mais

max_issues_repo_name: ProfLeao/mais

max_forks_repo_name: ProfLeao/mais


text:
using Test
module SingleTrait
    using BinaryTraits, Test
    struct Bird end

    @test check(SingleTrait, Bird).result == true # everything ok without traits defined

    @trait Fly
    @assign Bird with CanFly

    function test()
        @testset "Single Trait" begin
            @test istrait(FlyTrait) == true
            @test istrait(Int) == false
            @test supertype(FlyTrait) === Any
            @test supertype(CanFly) <: FlyTrait
            @test supertype(CannotFly) <: FlyTrait
            @test flytrait(Bird()) == CanFly()
        end
    end
end

module MultipleTraits
    using BinaryTraits, Test
    struct Duck end
    struct Dog end
    @trait Fly
    @trait Swim
    @assign Duck with CanFly, CanSwim
    @assign Dog with CanSwim
    function test()
        @testset "Multiple Traits" begin
            @test flytrait(Dog()) == CannotFly()
            @test swimtrait(Dog()) == CanSwim()
            @test flytrait(Duck()) == CanFly()
            @test swimtrait(Duck()) == CanSwim()
        end
    end
end

module TraitSuperType
    using BinaryTraits, Test
    struct Duck end
    abstract type Mobility end
    @trait Fly as Mobility
    function test()
        @testset "Super Type" begin
            @test supertype(FlyTrait) <: Mobility
        end
    end
end

module CustomPrefixes
    using BinaryTraits, Test
    @trait Iterable prefix Is,Not
    @assign AbstractArray with IsIterable
    next(x) = next(iterabletrait(x), x)
    next(::IsIterable, x) = iterate(x)
    next(::NotIterable, x) = :toobad
    function test()
        @testset "Custom Prefix" begin
            @test next([1,2,3]) !== nothing
            @test next(:hello) === :toobad
        end
    end
end

module CompositeTraits
    using BinaryTraits, Test

    # possible traits
    @trait Move
    @trait CarryPassenger
    @trait FourWheels prefix Has,No
    @trait Engine prefix Has,No

    # assignments
    struct Acura end
    @assign Acura with CanMove, CanCarryPassenger, HasFourWheels, HasEngine
    struct Tricycle end
    @assign Tricycle with CanMove, CanCarryPassenger

    # composite
    @trait Car prefix Is,Not with CanMove, CanCarryPassenger, HasFourWheels, HasEngine

    function test()
        @testset "Composite" begin
            @test cartrait(Acura()) == IsCar()
            @test cartrait(Tricycle()) == NotCar()
        end
    end
end

module SyntaxErrors
    using BinaryTraits, Test
    using BinaryTraits: SyntaxError

    # This macro expands to testing code that returns true when syntax error
    # is detected properly
    macro testme(ex)
        esc(quote
                try
                    @macroexpand($ex)
                    false
                catch e
                    # @info "expansion error" e
                    e isa LoadError
                end
            end)
    end

    @trait Eat
    @trait Drink
    struct Dog end

    function test()
        @testset "Syntax" begin
            @test @testme @trait Fly as                # missing type
            @test @testme @trait Fly as 1              # 1 is not a type
            @test @testme @trait Fly prefix            # missing prefix tuple
            @test @testme @trait Fly prefix 1,2        # wrong type
            @test @testme @trait Fly prefix Is         # needs two symbols
            @test @testme @trait Fly prefix Is,Not,Ha  # needs two symbols
            @test @testme @trait Fly with Eat          # must have at least 2 sub-traits
            @test @testme @trait Fly with Eat,Drink,1  # 1 is not a symbol

            @test @testme @assign Dog with 1           # 1 is not a symbol
            # invalid argument
            @test @testme @implement CanCreep by creep2(_, []) # invalid argument
        end
    end
end

module Interfaces
    using BinaryTraits, Test
    using BinaryTraits: SyntaxError, extract_type

    const SUPPORT_KWARGS = VERSION >= v"1.2"
    const mod = @__MODULE__

    struct Bird end
    # Fly trait requires multiple contracts with variety of func signatures
    @trait Fly
    @assign Bird with CanFly
    bird_check = @check(Bird)
    @test bird_check.result == true # everything ok without interface contracts

    @implement CanFly by liftoff(_)
    @implement CanFly by speed(_, resistence::Float64)::Float64
    @implement CanFly by flyto(::Float64, ::Float64, _)::String  # fly to (x,y)

    # Pretty trait requires a single contract
    @trait Pretty prefix Is,Not
    @implement IsPretty by look_at_the_mirror_daily(_)::Bool

    # Bird satisfies all contracts from FlyTrait by concrete type
    liftoff(::Bird) = "hi ho!"
    speed(::Bird, resistence::Float64) = 100 - resistence
    flyto(x::Float64, y::Float64, ::Bird) = "Arrvied at ($x, $y)"

    # Duck satisfies partial contracts from FlyTrait by concrete type
    struct Duck end
    @assign Duck with CanFly
    liftoff(::Duck) = "hi ho!"

    # Chicken does not satisfy any contract from FlyTrait
    struct Chicken end
    @assign Chicken with CanFly

    # Flamingo exhibits both Fly and Pretty traits
    struct Flamingo end
    @assign Flamingo with CanFly, IsPretty  # composite trait
    liftoff(::Flamingo) = "wee!"
    speed(::Flamingo, resistence::Float64) = 150 - resistence
    flyto(x::Float64, y::Float64, ::Flamingo) = "Arrvied at ($x, $y)"
    look_at_the_mirror_daily(::Flamingo) = true

    # Test composite traits - total underlying 4 contracts required for this!
    @trait FlyPretty prefix Is,Not with CanFly, IsPretty

    # Crane partially satisfies FlyTrait and fully satisfies Pretty trait
    struct Crane end
    @assign Crane with IsFlyPretty
    speed(::Crane, resistence::Float64) = 150 - resistence
    flyto(x::Float64, y::Float64, ::Crane) = "Arrvied at ($x, $y)"
    look_at_the_mirror_daily(::Crane) = true

    struct Penguin end
    @trait Dive
    @assign Penguin with CanDive
    @implement CanDive by dive1(_, ::Integer)           # no argument name
    @implement CanDive by dive2(_, ::Vector{<:Integer}) # parameterized type
    if SUPPORT_KWARGS
        @implement CanDive by dive31(_, x::Real;)            # keyword arguments
        @implement CanDive by dive32(_, x::Real; kw::Real)   # keyword arguments
        @implement CanDive by dive33(_, y; kw1::Real, kw2)   # keyword arguments
        @implement CanDive by dive34(_, x; kw1)
    end
    @implement CanDive by dive4(_, ::Base.Bottom)
    @implement CanDive by dive5(_, x)
    @implement CanDive by dive6(_, ::Number)

    dive1(::Penguin, ::Real) = 1                # Real >: Integer
    dive2(::Penguin, ::Vector) = 2              # Vector >: Vector{<:Integer}
    dive31(::Penguin, ::Number) = 31            # no kw argument
    dive32(::Penguin, ::Number; kw::Complex) = 32 # kw argument type is ignored!
    dive33(::Penguin, ::Int; kw...) = 33        # any number of kw arguments
    dive34(::Penguin, ::Float64) = 34           # keyword argument missing
    dive4(::Penguin, ::Integer) = 4             # Integer >: Base.Bottom
    dive5(::Penguin, ::Int) = 5                 # Int >: Bottom
    dive6(::Penguin, ::Int) = 6                 # not Int >: Number

    # Issue #30 - propagate interfaces according to type inheritance
    abstract type Animal end
    struct Rabbit <: Animal end
    @trait Eat
    @assign Animal with CanEat
    @implement CanEat by eat(_)
    eat(::Animal) = 1

    # no contract requirements (code coverage)
    struct Kiwi end

    # weird argument types in contract specs
    @trait Creep
    @implement CanCreep by creep1(_, a::Int=5) # argument assignment

    struct Snake end
    @assign Snake with CanCreep
    creep1(::Snake, ::Integer) = 1

    function test()
        @testset "Interface validation" begin

            bird_check = @check(Bird)
            @test bird_check.result == true
            @test bird_check.implemented |> length == 3
            @test bird_check.misses |> length == 0

            chicken_check = @check(Chicken)
            @test chicken_check.result == false
            @test chicken_check.implemented |> length == 0
            @test chicken_check.misses |> length == 3

            duck_check = @check(Duck)
            @test duck_check.result == false
            @test duck_check.implemented |> length == 1
            @test duck_check.misses |> length == 2

            flamingo_check = @check(Flamingo)
            @test flamingo_check.result == true
            @test flamingo_check.implemented |> length == 4
            @test flamingo_check.misses |> length == 0

            crane_check = @check(Crane)
            @test crane_check.result == false
            @test crane_check.implemented |> length == 3
            @test crane_check.misses |> length == 1

            penguin_check = @check(Penguin)
            @test penguin_check.result == false
            @test penguin_check.implemented |> length == (SUPPORT_KWARGS ? 7 : 4)
            @test penguin_check.misses |> length == (SUPPORT_KWARGS ? 2 : 1)

            rabbit_check = @check(Rabbit)
            @test rabbit_check.result == true

            # test `show` function
            buf = IOBuffer()
            contains(s) = x -> occursin(s, x)

            show(buf, flamingo_check)
            @test buf |> take! |> String |> contains("has implemented")

            show(buf, bird_check)
            @test buf |> take! |> String |> contains("has implemented")

            show(buf, duck_check)
            @test buf |> take! |> String |> contains("is missing")

            # Bird is assigned with 1 FlyTrait and that requires 3 contracts
            @test required_contracts(mod, Bird) |> length == 3

            # Crane requires 4 contracts because it has both Fly and Pretty traits
            @test required_contracts(mod, Crane) |> length == 4

            # Penguin
            if SUPPORT_KWARGS
                show(buf, penguin_check.misses)
                @test buf |> take! |> String |> contains("dive34")
            end
            show(buf, penguin_check.misses)
            @test buf |> take! |> String |> contains("dive6")

            # has no interface requirements
            check_kiwi = @check(Kiwi)
            @test check_kiwi.result
            show(buf, check_kiwi)
            @test buf |> take! |> String |> contains("has no interface contract")

            # strange argument types are both accepted
            check_snake = @check(Snake)
            @test check_snake.result
            @test check_snake.implemented |> length == 1

            # code coverage
            @test_throws SyntaxError extract_type(:([]), :(CanFly), nothing)

        end
    end
end

module Verbose
    using BinaryTraits, Test, Logging

    # Capture output and make sure that `expected` appears in the log
    # when expression `ex` is evaluated
    macro testme(expected, ex)
        esc(quote
            buf = IOBuffer()
            with_logger(ConsoleLogger(buf)) do
                @macroexpand($ex)
            end
            s = String(take!(buf))
            @test occursin($expected, s)
        end)
    end

    @trait Scratch
    struct Cat end

    # Testing verbose mode
    function test()
        BinaryTraits.set_verbose(true)
        @testset "Verbose" begin
            @testme "struct CanScratch" @trait Scratch
            @testme "scratchtrait(::Cat) = CanScratch()" @assign Cat with CanScratch
        end
    end
end

module CrossModule
    using Test, Logging, BinaryTraits
    module X
        using BinaryTraits
        @trait RowTable prefix Is,Not
        @implement IsRowTable by row(_, ::Integer)
        __init__() = inittraits(@__MODULE__)
    end

    module Y
        using Test
        using BinaryTraits
        using ..X
        struct AwesomeTable end
        @assign AwesomeTable with X.IsRowTable
        r = @check(AwesomeTable)
        @test r.implemented |> length == 0
        X.row(::AwesomeTable, ::Number) = 1
        __init__() = inittraits(@__MODULE__)
    end

    function test()
        @testset "cross-module implementation" begin
            r = @check(Y.AwesomeTable)
            @test r.implemented |> length == 1
        end
    end
end

@testset "BinaryTraits Tests" begin
    import .SingleTrait;        SingleTrait.test()
    import .MultipleTraits;     MultipleTraits.test()
    import .TraitSuperType;     TraitSuperType.test()
    import .CustomPrefixes;     CustomPrefixes.test()
    import .CompositeTraits;    CompositeTraits.test()
    import .SyntaxErrors;       SyntaxErrors.test()
    import .Interfaces;         Interfaces.test()
    import .Verbose;            Verbose.test()
    import .CrossModule;        CrossModule.test()
end


text:
facts("Assign ranks within tolerance") do

  context("Ranks correctly if none are within tolerance of each other") do
    res = BlackBoxOptim.Utils.assign_ranks_within_tolerance([1.0, 2.0, 3.0])
    @fact length(res) => 3
    @fact res[1] => (1, 1.0, 1.0)
    @fact res[2] => (2, 2.0, 2.0)
    @fact res[3] => (3, 3.0, 3.0)

    res = BlackBoxOptim.Utils.assign_ranks_within_tolerance([3.0])
    @fact length(res) => 1
    @fact res[1] => (1, 3.0, 3.0)
  end

  context("Ranks correctly if some are within tolerance of each other") do
    res = BlackBoxOptim.Utils.assign_ranks_within_tolerance([1.0, 1.0+1e-6, 200.0]; tolerance = 1e-5)
    @fact length(res) => 3
    @fact res[1] => (1, 1.0, 1.0)
    @fact res[2] => (1, 1.0+1e-6, 1.0+1e-6)
    @fact res[3] => (3, 200.0, 200.0)
  end

  context("Ranks correctly if all are within tolerance of each other") do
    res = BlackBoxOptim.Utils.assign_ranks_within_tolerance([1.0, 1.0+1e-6, 1.0-1e-6]; tolerance = 1e-5)
    @fact length(res) => 3
    @fact res[1] => (1, 1.0-1e-6, 1.0-1e-6)
    @fact res[2] => (1, 1.0, 1.0)
    @fact res[3] => (1, 1.0+1e-6, 1.0+1e-6)
  end

  context("Ranks in reverse if none are within tolerance of each other") do
    res = BlackBoxOptim.Utils.assign_ranks_within_tolerance([1.0, 20.0, 13.0]; rev = true)
    @fact length(res) => 3
    @fact res[1] => (1, 20.0, 20.0)
    @fact res[2] => (2, 13.0, 13.0)
    @fact res[3] => (3, 1.0, 1.0)
  end

  context("Ranks correctly for a complex example testing many aspects") do
    values = shuffle([(-11.0, :a), (1.0, :b), (1.0+1e-6, :c), (1.0-1e-6, :d),
      (2000, :e), (3000, :f)])
    res = BlackBoxOptim.Utils.assign_ranks_within_tolerance(values;
      tolerance = 1e-5, by = (p) -> p[1], rev = true)
    @fact length(res) => length(values)
    @fact res[1] => (1, (3000, :f), 3000)
    @fact res[2] => (2, (2000, :e), 2000)
    @fact res[3] => (3, (1.0+1e-6, :c), 1.0+1e-6)
    @fact res[4] => (3, (1.0, :b), 1.0)
    @fact res[5] => (3, (1.0-1e-6, :d), 1.0-1e-6)
    @fact res[6] => (6, (-11.0, :a), -11.0)
  end

end


text:
function DelaunayArguments(arg0::AbsoluteDate, arg1::jdouble, arg10::jdouble, arg11::jdouble, arg12::jdouble, arg13::jdouble, arg2::jdouble, arg3::jdouble, arg4::jdouble, arg5::jdouble, arg6::jdouble, arg7::jdouble, arg8::jdouble, arg9::jdouble)
    return DelaunayArguments((AbsoluteDate, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg10, arg11, arg12, arg13, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9)
end


text:
using ITensors
using Test

@testset "Test argsdict function" begin
  empty!(ARGS)
  push!(
    ARGS,
    "x",
    "n = 1",
    "nf :: Float64 = 2",
    "ni :: AutoType = 2",
    "ns :: String = 2",
    "nc :: ComplexF64 = 3",
    "x = 2e-1, 2e-3, 0.1",
    "2::AutoType",
    "N = 1e-3",
    "vc = ComplexF64[1 + 2im, 3]",
    "3",
    "--autotype",
    "vf = [1.0, 3]",
    "y = true",
    "1+2im",
    "s = \"use_qr\"",
    "--stringtype",
    "y",
  )
  args = argsdict()
  empty!(ARGS)

  @test args["_arg1"] == "x"
  @test args["nf"] == 2.0
  @test args["nc"] == 3.0 + 0.0im
  @test args["ns"] == "2"
  @test args["ni"] == 2
  @test args["y"] == true
  @test args["N"] == 0.001
  @test args["x"] == (0.2, 0.002, 0.1)
  @test args["_arg2"] == 2
  @test args["vc"] == Complex{Float64}[1.0 + 2.0im, 3.0 + 0.0im]
  @test args["_arg3"] == "3"
  @test args["vf"] == [1.0, 3.0]
  @test args["n"] == 1
  @test args["_arg4"] == 1 + 2im
  @test args["s"] == "use_qr"
  @test args["_arg5"] == "y"

  push!(
    ARGS,
    "x",
    "n -> 1",
    "nf :: Float64 -> 2",
    "ni :: AutoType -> 2",
    "ns :: String -> 2",
    "nc :: ComplexF64 -> 3",
    "x -> 2e-1, 2e-3, 0.1",
    "2",
    "N -> 1e-3",
    "vc -> ComplexF64[1 + 2im, 3]",
    "3",
    "--autotype",
    "vf -> [1.0, 3]",
    "y -> true",
    "1+2im",
    "s -> \"use_qr\"",
    "--stringtype",
    "y",
  )
  args = argsdict(;
    first_arg=2,
    delim="->",
    as_symbols=true,
    default_named_type=String,
    default_positional_type=ITensors.AutoType,
    prefix="test",
  )
  empty!(ARGS)

  @test args[:nf] == 2.0
  @test args[:ni] == 2
  @test args[:ns] == "2"
  @test args[:nc] == 3.0 + 0.0im
  @test args[:y] == true
  @test args[:N] == "1e-3"
  @test args[:x] == "2e-1, 2e-3, 0.1"
  @test args[:test1] == 2
  @test args[:vc] == "ComplexF64[1 + 2im, 3]"
  @test args[:test2] == 3
  @test args[:vf] == [1.0, 3.0]
  @test args[:n] == "1"
  @test args[:test3] == 1 + 2im
  @test args[:s] == "use_qr"
  @test args[:test4] == "y"

  #
  # Check for some syntax errors
  #

  push!(ARGS, "x y=2")
  @test_throws ErrorException argsdict()
  empty!(ARGS)

  push!(ARGS, "x=y=2")
  @test_throws ErrorException argsdict()
  empty!(ARGS)

  push!(ARGS, "x::MyType = 2")
  @test_throws UndefVarError argsdict()
  empty!(ARGS)

  push!(ARGS, "x = y")
  @test_throws UndefVarError argsdict()
  empty!(ARGS)
end


text:
""" proposal(samples; transform=identity)
returns an MvNormal distribution and the log determinant of the jacobian matrix
"""

## MvNormal proposal
function proposal(samples::AbstractMatrix)
    μ = mean(samples, dims=2)
    Σ = cov(samples')
    return MvNormal(vec(μ), Σ)
end

## TODO : WRAP-3 proposal


text:
"""
**Represents an occurrence in the GBIF format**

This is currently a subset of all the fields. This `struct` is *not* mutable --
this ensures that the objects returned from the GBIF database are never modified
by the user.

The `taxon` field is a `GBIFTaxon` object, and can therefore be manipulated as
any other `GBIFTaxon`.
"""
struct GBIFRecord
    key::Int64
    datasetKey::AbstractString
    dataset::Union{Missing, AbstractString}
    publishingOrgKey::Union{Missing, AbstractString}
    publishingCountry::Union{Missing, AbstractString}
    institutionCode::Union{Missing, AbstractString}
    protocol::Union{Missing, AbstractString}
    countryCode::Union{Missing, AbstractString}
    country::Union{Missing, AbstractString}
    basisOfRecord::Symbol
    individualCount::Union{Missing, Integer}
    latitude::Union{Missing, AbstractFloat}
    longitude::Union{Missing, AbstractFloat}
    precision::Union{Missing, AbstractFloat}
    uncertainty::Union{Missing, AbstractFloat}
    geodetic::Union{Missing, AbstractString}
    date::Union{Missing, DateTime}
    identified::Union{Missing, DateTime}
    issues::Vector{Symbol}
    taxonKey::Union{Missing, Integer}
    rank::Union{Missing, AbstractString}
    taxon::GBIFTaxon
    generic::Union{Missing, AbstractString}
    epithet::Union{Missing, AbstractString}
    vernacular::Union{Missing, AbstractString}
    scientific::Union{Missing, AbstractString}
    observer::Union{Missing, AbstractString}
    license::Union{Missing, AbstractString}
end

"""
**Internal function to format dates in records**
"""
function format_date(o, d)
    t = get(o, d, missing)
    if t === nothing || ismissing(t)
        return missing
    else
        try
            return DateTime(t[1:19])
        catch
            return missing
        end
    end
end

"""
**Generates an occurrence from the JSON response of GBIF**
"""
function GBIFRecord(o::Dict{String, Any})
    # Prepare the taxon object
    levels = ["kingdom", "phylum", "class", "order", "family", "genus", "species"]
    r = Dict{Any,Any}()
    for l in levels
        if haskey(o, l)
            r[l] = Pair(o[l], o[l*"Key"])
        else
            r[l] = missing
        end
    end

    # The name of the taxon is the generic name + specific epithet, but if the
    # generic name is asbent, this defaults to the scientificName.
    if !ismissing(get(o, "genericName", missing))
        this_name = o["genericName"]
        if !ismissing(get(o, "specificEpithet", missing))
            this_name *= " "*get(o, "specificEpithet", missing)
        end
    else
        this_name = o["scientificName"]
    end

    this_record_taxon = GBIFTaxon(
        this_name,
        o["scientificName"],
        :ACCEPTED,
        :EXACT,
        r["kingdom"],
        r["phylum"],
        r["class"],
        r["order"],
        r["family"],
        r["genus"],
        r["species"],
        100,
        false
    )

    formatted_record =  GBIFRecord(
        o["key"],
        o["datasetKey"],
        get(o, "datasetName", missing),
        get(o, "publishingOrgKey", missing),
        get(o, "publishingCountry", missing),
        get(o, "institutionCode", missing),
        get(o, "protocol", missing),
        get(o, "countryCode", missing),
        get(o, "country", missing),
        Symbol(get(o, "basisOfRecord", "UNKNOWN")),
        get(o, "individualCount", missing),
        get(o, "decimalLatitude", missing),
        get(o, "decimalLongitude", missing),
        get(o, "precision", missing),
        get(o, "coordinateUncertaintyInMeters", missing),
        get(o, "geodeticDatum", missing),
        format_date(o, "eventDate"),
        format_date(o, "dateIdentified"),
        Symbol.(o["issues"]),
        get(o, "taxonKey", missing),
        get(o, "taxonRank", missing),
        this_record_taxon,
        get(o, "genericName", missing),
        get(o, "specificEpithet", missing),
        get(o, "vernacularName", missing),
        get(o, "scientificName", missing),
        get(o, "recordedBy", missing),
        get(o, "license", missing)
    )
    return formatted_record
end

"""
**List of occurrences and metadata**

This type has actually very few information: the `query` field stores the query
parameters. This type is mutable and fully iterable.

The `occurrences` field is pre-allocated, meaning that it will contain `#undef`
elements up to the total number of hits on GBIF. When iterating, this is taken
care of automatically, but this needs to be accounted for if writing code that
accesses this field directly.
"""
mutable struct GBIFRecords
    query::Union{Vector{Pair}, Nothing}
    occurrences::Vector{GBIFRecord}
end


text:
function test15()
  filename = "test15.asc"
  # exectuablepath = null string will not run LTspice.exe.  Test parsing only.
  sim = LTspiceSimulation(filename,executablepath="")
  @test LTspice.does_circuitfilearray_file_match(sim)
  show(IOBuffer(),sim)
  @test measurementnames(sim) == ("σ", "a", "ψπππππ")
  @test parameternames(sim) == ("m", "μ", "MEG", "Ω", "θ", "Δ", "Φ", "ψ")
end
test15()


text:
module SamplingDataFrame

    # Write your package code here.
    using DataFrames,StatsBase

    export TestDataSetsNumbering,CreateTestDataSets

    #テストデータと検証用データを作成するためのサンプリングを行なう
    include("TestDataSetsNumbering.jl")
end


text:
module BenchFilterMap

using BenchmarkTools
using Restacker: restack

@noinline function filter_map_norestack!(ys, xs)
    @inbounds for i in eachindex(ys, xs)
        x = xs[i]
        if -0.5 < x < 0.5
            ys[i] = 2x
        end
    end
end

@noinline function filter_map_restack!(ys, xs)
    ys = restack(ys)
    @inbounds for i in eachindex(ys, xs)
        x = xs[i]
        if -0.5 < x < 0.5
            ys[i] = 2x
        end
    end
end

suite = BenchmarkGroup()

xs = view(randn(10_000), :)
ys = view(randn(10_000), :)
suite["norestack"] = @benchmarkable filter_map_norestack!($ys, $xs)
suite["restack"] = @benchmarkable filter_map_restack!($ys, $xs)

end  # module
BenchFilterMap.suite


text:
### A Pluto.jl notebook ###
# v0.17.7

using Markdown
using InteractiveUtils

# ╔═╡ fffc6cac-c044-4cd9-a3b7-084341878faf
using HypertextLiteral

# ╔═╡ dd3e39ac-8b2c-11ec-1774-7543d2573dc0
md"# 5. The Mass Action Law

## 5.1 Definition of the Law of Mass Action

Now that we know how to obtain the stoichiometric matrices and the order of any given reaction, the next step is to write the differential equations that capture the dynamics, i.e., the changes in time of the amount of all the species involved.  To illustrate this, we will use a generic scheme of interaction in vector notation:

$$\begin{align*}
X_1 + X_2 \overset{k}{\longrightarrow} X_3  \tag{58}\\
\end{align*}$$
The stoichometrix matrix is: "

# ╔═╡ 7c7a696a-8716-4739-9a40-44a34c621485
md" The matrix A is [1,1,0]"

# ╔═╡ 30a79e11-2165-4908-94a3-d020cc6d2e16
md"The order of the reaction is $(sum([1,1,0]))"

# ╔═╡ e11754b1-6993-403b-a0ec-e0aaf86a5dd3
md"The Units of k are s^-1 x k^$(1-(sum([1,1,0])))"

# ╔═╡ ab0aa754-78e2-47b5-a398-a7e1ab8d6c58
md"We will write down the differential equations that account for the dynamics of each species.

```math
\begin{align}
\frac{\mathrm{d} X_1}{\mathrm{d} t}  &=&  \tag{59}\\
\frac{\mathrm{d} X_2}{\mathrm{d} t} &=& \tag{60}\\
\frac{\mathrm{d} X_3}{\mathrm{d} t} &=&  \tag{61}
\end{align}
```

The __Mass Action Law__ describes how a chemical reaction takes place under ideal and constant conditions. In brief, it is based on the fact that the speed of any reaction is proportional to the concentrations of the reactants involved, which has been discussed previously in this course (the proportionalinty constant is the kinetic rate constant). From here, it follows that the rate of change of any reactant is also proportional to the concetrations of the reactants involved. This velocity depends on the affinity of the reaction and the concentrations of each of the reactants (in any scheme of interactions, the ones at the left of the arrow). The affinity of the reactions is defined by the rate constant `k` of the reaction. The `k` is the proportionality constant that correlates the speed fo the reaction with the concentrations of the species involved, which in this case are $X_1$ and $X_2$:

```math
\begin{align}
\frac{\mathrm{d} X_1}{\mathrm{d} t}  &= - k X_1 X_2 \tag{62}\\
\frac{\mathrm{d} X_2}{\mathrm{d} t} &=- k X_1 X_2 \tag{63}\\
\frac{\mathrm{d} X_3}{\mathrm{d} t} &=  k X_1 X_2\tag{64}
\end{align}
```

Lets see if the units are correct:

```math
\begin{align}
\frac{[M]}{[s]}  &= - \frac{1}{[s][M]}[M] [M] \tag{65}\\
\frac{[M]}{[s]}  &= - \frac{1}{[s][M]}[M] [M] \tag{66}\\
\frac{[M]}{[s]}  &=  \frac{1}{[s][M]}[M] [M] \tag{67}\\
\end{align}
```

For the case where more than one molecule of reactant is involved, we will use the reaction of formation of water.
"

# ╔═╡ d2e26cb1-f136-4b96-8ef3-d3ee30f61b3d
md" The matrix A is [2 1 0]"

# ╔═╡ 10f98ec7-fcad-4587-90f8-62f7b1b74f4d
md" The matrix B is [0 0 2]"

# ╔═╡ fee636f1-4944-46a4-bcbe-32b7273268c7
md"The order of the reaction is $(sum([2 1 0]))"

# ╔═╡ 56a68cf6-6e64-4875-b614-18987857bc06
md"The Units of k are s^-1 x k^$(1-(sum([2 1 0])))"

# ╔═╡ 05faae5d-fdb6-4f1d-a77f-0a1337e1666d
StMat=([0 0 2]-[2 1 0])'

# ╔═╡ 0d524884-afbf-4594-80ee-9fdfcdabbea0
md"we can decompose the equation as

$$\begin{align*}
X_1 +  X_1 + X_2 \overset{k}{\longrightarrow}  X_3 + X_3 \tag{68}\\
\end{align*}$$

We will write down the differential equations that account for the dynamics of each species. The speed of each reaction depends on the concentration of the species involved (what is before the arrow). In this case is two units of $X_1$ and one unit of $X_2$. In addition, since the reaction consumes 2 units of $X_1$ for each unit of $X_2$, the rate of consumption of $X_1$ is twice the rate of consumption $X_2$. Another more automaticic way of doing it is to use the matrix of stochoimetric coefficients $(B-A)'$. Therefore:

```math
\begin{align}
\frac{\mathrm{d} X_1}{\mathrm{d} t}  &= - 2k X_1^2 X_2 \tag{69}\\
\frac{\mathrm{d} X_2}{\mathrm{d} t} &=- k X_1^2 X_2 \tag{70}\\
\frac{\mathrm{d} X_3}{\mathrm{d} t} &=2  k X_1^2 X_2 \tag{71}
\end{align}
```

Lets check if units are correct:

```math
\begin{align}
\frac{[M]}{[s]}  &= - \frac{1}{[s][M]^2} [M]^2 [M] \tag{72}\\
\frac{[M]}{[s]}  &= - \frac{1}{[s][M]^2} [M]^2 [M] \tag{73}\\
\frac{[M]}{[s]}  &= \frac{1}{[s][M]^2} [M]^2 [M] \tag{74}
\end{align}
```

"

# ╔═╡ 61a1a25f-21b1-40a6-aafa-17fddf4a848c
md"## 5.2 General formulation of the Mass Action Law

When we have to write the differenetial equations of more complex systems of interactions, such as Eq 21, we can use a general method based omn the following equation:
```math
\begin{align}
\frac{\mathrm{d} X}{\mathrm{d} t}&= (B-A)^T \cdot K \cdot X^A \tag{75}\\
\end{align}
```

where  `A` and `B` are the matrices with the stoichiometric coefficients, `X` is the state vector $X=[X_1,X_2,X_3]^T$, and `K` is a matrix in the form:

```math
K=\begin{pmatrix}
 k_1 & 0   & 0  & \dots & 0 \\
 0 &  k_2  & 0  & \dots & 0\\
 0 & 0  & k_3 &\dots  &0 \\
 \vdots &  \vdots   &  \vdots  &  \vdots  &  \vdots  \tag{76}\\
 0 & 0  & 0  &  \dots & k_r
\end{pmatrix}
```

$X^A$ is a column vector, calculated as:

```math
X^A=\begin{pmatrix}
X_1^{A_{11}} \cdot X_2^{A_{12}} \cdot X_3^{A_{13}} \cdot  ... \cdot  X_s^{A_{1s}} \\
X_1^{A_{21}} \cdot X_2^{A_{22}} \cdot X_3^{A_{23}} \cdot  ... \cdot  X_s^{A_{2s}} \\
 \vdots  \,\,\,\,\,\,\,\,\,\,\,  \vdots  \,\,\,\,\,\,\,\,\,\,\,\,    \vdots \,\,\,\,\,\,\,\,\,   \vdots  \,\,\,\,\,\,\,\,   \vdots  \\
X_1^{A_{r1}} \cdot X_2^{A_{r2}} \cdot X_3^{A_{r3}} \cdot  ... \cdot  X_s^{A_{rs}} \\
\end{pmatrix} =\begin{pmatrix}
 \prod_{i=1}^{s} X_{i}^{A_{1i}}\\
\prod_{i=1}^{s} X_{i}^{A_{2i}}\\
   \vdots  \tag{77} \\
\prod_{i=1}^{s} X_{i}^{A_{ri}}\\
\end{pmatrix}
```

"

# ╔═╡ 2b58c174-c501-47ed-ade6-0ee9b88b76eb
md" ### 5.3 Example of differential form of Mass Action Law

As an example we will use the previous reaction for the formation of salt.

```math
NaCO_3 + CaCl_2 \overset{k_1}{\underset{k_2}{\longleftrightarrow}} CaCO_3 + 2 \cdot NaCl \tag{78}\\
```

with two reactions $r=2$ and four species $s=4$, with the stoichiometric matrices `A` and `B` listed in Eq. 30"

# ╔═╡ 6c491ea3-10b9-46b3-a3ff-2c27d029efe0
begin
	A=[1 1 0 0; 0 0 1 2]
	B=[0 0 1 2; 1 1 0 0]
	stoichiometric_matrix= (B-A)'
	println("The stoichiometric matrix is $stoichiometric_matrix ")
end

# ╔═╡ eb79f46c-5b6d-4b36-a6e0-f92fe6938b46
md"in our case

```math
K=\begin{pmatrix}
 k_1 & 0   \tag{79}\\
 0 &  k_2  \\
\end{pmatrix}
```

and

```math
X^A=\begin{pmatrix}
X_1\cdot X_2  \tag{80}\\
X_3\cdot X_4^2
\end{pmatrix}
```

therefore

```math
\begin{align}
 \begin{bmatrix}
\frac{\mathrm{d} X_1}{\mathrm{d} t}\\ \frac{\mathrm{d} X_2}{\mathrm{d} t} \\ \frac{\mathrm{d} X_3}{\mathrm{d} t} \\ \frac{\mathrm{d} X_4}{\mathrm{d} t}
\end{bmatrix}&= \begin{bmatrix}
 -1   &  1 \\
    -1  &   1\\
     1  &  -1\\
     2  &  -2
\end{bmatrix}\begin{bmatrix}
k_1 & 0\\
 0& k_2
\end{bmatrix}\begin{bmatrix}
X_1 \cdot X_2\\
X_3 \cdot X_4^2
\end{bmatrix} \tag{81}
\end{align}
```

so calculating

```math
\begin{align}
\begin{bmatrix}
\frac{\mathrm{d} X_1}{\mathrm{d} t}\\ \frac{\mathrm{d} X_2}{\mathrm{d} t} \\ \frac{\mathrm{d} X_3}{\mathrm{d} t} \\ \frac{\mathrm{d} X_4}{\mathrm{d} t}
\end{bmatrix}&= \begin{bmatrix}
 -1   &  1 \\
    -1  &   1\\
     1  &  -1\\
     2  &  -2
\end{bmatrix}\begin{bmatrix}
k_1 \cdot  X_1 \cdot X_2\\
k_2 \cdot X_3 \cdot X_4^2
\end{bmatrix} \tag{82}
\end{align}
```

In other words, the differential equations are obtained by multiplying the matrix $(B-A)'$ by a vector of speeds of the reactions involved. This way, the final set of equations is:

```math
\begin{align}
            \frac{ \mathrm{d}X_1 }{\mathrm{d}t} &= - k_1 \cdot X_1 \cdot X_2 + k_2 \cdot X_3 \cdot X_4^2  \tag{83}\\
            \frac{ \mathrm{d}X_2 }{\mathrm{d}t} &= - k_1 \cdot X_1 \cdot X_2 + k_2 \cdot X_3 \cdot X_4^2   \tag{84}  \\
             \frac{ \mathrm{d} X_3 }{\mathrm{d}t} &=  k_1 \cdot X_1 \cdot X_2 - k_2 \cdot X_3 \cdot X_4^2  \tag{85}  \\
              \frac{ \mathrm{d} X_4 }{\mathrm{d}t} &= 2 k_1 \cdot X_1  \cdot X_2 - 2 k_2 \cdot X_3 \cdot X_4^2   \tag{86}  \\
\end{align}
```
"

# ╔═╡ 3f13dd7b-5430-41b2-a6c2-e185f924ba6b
@htl("""

<div class='blue-background'>
Hello!
</div>

<script>
// more about selecting elements later!
currentScript.previousElementSibling.innerText = "Computer Task 4: Differential Equations '"

</script>

<style>
.blue-background {
	padding: .5em;
	background: lightblue;
	color: black;
}
</style>

""")

# ╔═╡ c5a8887c-6875-481d-97bc-cb7da256d2fb
md" 1. Find the order of the reactions, the units of the kinetic rate constants, and write the differential equations that govern the dynamics of the following system:

```math
X_1  \overset{k_1}{\underset{k_2}{\longleftrightarrow}} X_2
```

2. Find the order of the reactions, the units of the kinetic rate constants, and write the differential equations that govern the dynamics of the following system:

```math
X_1 + X_2 \overset{k_1}{\underset{k_2}{\longleftrightarrow}} 2 X_1
```

3.  Find the order of the reactions, the units of the kinetic rate constants, and write the differential equations that govern the dynamics of the following system:

```math
X_1\overset{k_1}{\longrightarrow}X_2\overset{k_2}{\longrightarrow}X_3
```
"

# ╔═╡ fba4e0b3-52a5-4d1a-b7da-6740f2c2a836
md"## Solution Task 1

Decomposing the reversible reaction into two irreversible reactions
```math
\begin{align*}
X_1  \overset{k_1}{\underset{}{\longrightarrow}} X_2 \\
X_2  \overset{k_2}{\underset{}{\longrightarrow}} X_1
\end{align*}
```

let's put the stoichiometric values

```math
\begin{align*}
1 \cdot  X_1 + 0 \cdot X_2 \overset{k_1}{\underset{}{\longrightarrow}} 0 \cdot X_1+ 1 \cdot X_2 \\
0 \cdot X_1+ 1 \cdot X_2 \overset{k_2}{\underset{}{\longrightarrow}} 1 \cdot  X_1 + 0 \cdot X_2
\end{align*}
```

the stoichiometric matrices are:

"

# ╔═╡ 0731c83b-3b56-4b80-9654-1760dcbfca9b
AA=[1 0; 0 1];

# ╔═╡ 51ac8bbb-314d-4e82-ab3c-4d65be6f5e74
BB=[0 1; 1 0];

# ╔═╡ 7155b893-ec48-4213-b763-fe4d604efdaa
(BB-AA)'

# ╔═╡ e6a154b3-1f4f-4759-b0ab-62429b0c5200
function Calculate_Keq_stoichiometric_generic(A,B)

 for i=1:size(A)[1]
    println("The order of the reaction $i is ", sum(A[i,:]))
    println("The Units of k$i are s^-1 . k^$(1-(sum(A[i,:])))")
 end
    stoichiometric_matrix= (B-A)'
    println("The stoichiometric matrix is $stoichiometric_matrix ")

end

# ╔═╡ 7206af1d-1d79-4c18-82c2-7979d3baa2d5
Calculate_Keq_stoichiometric_generic(AA,BB)

# ╔═╡ b23e2b04-83a3-46f3-86e4-f49d8d42314d
md"The $K$ matrix and the $X^A$ vector required for the equation for the general form of the Mass action are:

```math
K=\begin{pmatrix}
 k_1 & 0\\
 0 & k_2
\end{pmatrix} \tag{26}
```

```math
X^A=\begin{pmatrix}
X_1^1\cdot X_2^0\\
X_1^0\cdot X_2^1
\end{pmatrix} = \begin{pmatrix}
X_1\\
 X_2
\end{pmatrix}  \tag{27}
```

so

```math
\begin{align}
 \begin{bmatrix}
\frac{\mathrm{d} X_1}{\mathrm{d} t}\\ \frac{\mathrm{d} X_2}{\mathrm{d} t} \end{bmatrix}= \begin{bmatrix}
-1 & 1\\
 1 & -1   \end{bmatrix} \begin{pmatrix}
 k_1 & 0\\
 0 & k_2
\end{pmatrix}\begin{pmatrix}
X_1\\
 X_2
\end{pmatrix}
\end{align}
```

```math
\begin{align}
\begin{bmatrix}
\frac{\mathrm{d} X_1}{\mathrm{d} t}\\ \frac{\mathrm{d} X_2}{\mathrm{d} t} \end{bmatrix}=  \begin{bmatrix}
-1 & 1\\
 1 & -1   \end{bmatrix}\begin{pmatrix}
k_1 \cdot X_1\\
  k_2 \cdot X_2
\end{pmatrix} =   \begin{bmatrix}
- k_1 \cdot X_1 +  k_2 \cdot X_2\\
 k_1 \cdot X_1 -k_2 \cdot X_2   \end{bmatrix}
\end{align} \tag{28}
```

The differential equations of the system can be written using the differential form of Mass Action Law as:

```math
\begin{align}
\frac{dX_1}{dt} = -k_1 \cdot X_1 +k_2 \cdot X_2 \tag{29}\\
\frac{dX_2}{dt} = k_1 \cdot X_1 - k_2 \cdot X_2 \tag{30}\\
\end{align}
```

	"

# ╔═╡ 0ca58de1-28bd-4df1-88ca-efa8d17250f6
md"## Solution Task 2

Decomposing the reversible reaction into two irreversible reactions
```math
\begin{align}
X_1  + X_2 \overset{k_1}{\underset{}{\longrightarrow}} 2 X_1 \\
 2 X_1  \overset{k_2}{\underset{}{\longrightarrow}}  X_1  + X_2
\end{align}
```

lest put the stoichiometric values

```math
\begin{align}
1 \cdot  X_1 + 1 \cdot X_2 \overset{k_1}{\underset{}{\longrightarrow}} 2 \cdot X_1+ 0 \cdot X_2 \\
2 \cdot X_1+ 0 \cdot X_2 \overset{k_2}{\underset{}{\longrightarrow}} 1 \cdot  X_1 + 1 \cdot X_2
\end{align}
```

the stoichiometric matrices are:"

# ╔═╡ 52a42411-35b9-472a-bd7e-99e7e6395015
AAA=[1 1; 2 0];

# ╔═╡ 637ae77a-a61c-4742-b737-aa4985ce6a32
BBB=[2 0; 1 1];

# ╔═╡ 45f11cc5-59bd-44a6-977b-78e1bea3d9c5
(BBB-AAA)'

# ╔═╡ a21e1203-57a8-4677-968b-c0a90cec9403
Calculate_Keq_stoichiometric_generic(AAA,BBB)

# ╔═╡ a796442b-85fa-45af-8c90-9ce9a8d4fcd2
md"The $K$ matrix and the $X^A$ vector required for the equation for the general form of the Mass action are:

```math
K=\begin{pmatrix}
 k_1 & 0\\
 0 & k_2
\end{pmatrix} \tag{26}
```

```math
X^A=\begin{pmatrix}
X_1^1\cdot X_2^1\\
X_1^2\cdot X_2^0
\end{pmatrix} = \begin{pmatrix}
X_1 \cdot X_2\\
 X_1^2
\end{pmatrix}  \tag{27}
```

so

```math
\begin{align}
 \begin{bmatrix}
\frac{\mathrm{d} X_1}{\mathrm{d} t}\\ \frac{\mathrm{d} X_2}{\mathrm{d} t} \end{bmatrix}= \begin{bmatrix}
1 & -1\\
 -1 & 1   \end{bmatrix} \begin{pmatrix}
 k_1 & 0\\
 0 & k_2
\end{pmatrix}\begin{pmatrix}
X_1 \cdot X_2\\
 X_1^2
\end{pmatrix}
\end{align}
```

```math
\begin{align}
\begin{bmatrix}
\frac{\mathrm{d} X_1}{\mathrm{d} t}\\ \frac{\mathrm{d} X_2}{\mathrm{d} t} \end{bmatrix}=  \begin{bmatrix}
1 & -1\\
 -1 & 1   \end{bmatrix}\begin{pmatrix}
k_1 \cdot X_1 \cdot X_2\\
  k_2 \cdot X_1^2
\end{pmatrix} =   \begin{bmatrix}
 k_1 \cdot X_1 \cdot X_2 -  k_2 \cdot X_1^2\\
- k_1 \cdot X_1 \cdot X_2 + k_2 \cdot X_1^2   \end{bmatrix}
\end{align} \tag{28}
```

The differential equations of the system can be written using the differential form of Mass Action Law as:

```math
\begin{align}
\frac{dX_1}{dt} =  k_1 \cdot X_1 \cdot X_2 -  k_2 \cdot X_1^2 \tag{29}\\
\frac{dX_2}{dt} = - k_1 \cdot X_1 \cdot X_2 + k_2 \cdot X_1^2 \tag{30}\\
\end{align}
```
"

# ╔═╡ 3ac754bf-e6de-417f-9093-35a73f036a1e
md"## Solution Task 3

Decomposing the reversible reaction into two irreversible reactions
```math
\begin{align}
X_1   \overset{k_1}{\underset{}{\longrightarrow}}  X_2 \\
X_2  \overset{k_2}{\underset{}{\longrightarrow}}  X_3
\end{align}
```

lest put the stoichiometric values

```math
\begin{align}
1 \cdot  X_1 + 0 \cdot X_2 + 0 \cdot X_3 \overset{k_1}{\underset{}{\longrightarrow}} 0 \cdot X_1+ 1 \cdot X_2 + 0 \cdot X_3\\
 0 \cdot X_1+ 1 \cdot X_2 + 0 \cdot X_3 \overset{k_2}{\underset{}{\longrightarrow}}  0 \cdot X_1+ 0 \cdot X_2 + 1 \cdot X_3
\end{align}
```

the stoichiometric matrices are:"

# ╔═╡ 6c027717-d272-4a3d-824c-54ea64f4fd92
AAAA=[1 0 0 ; 0 1 0];

# ╔═╡ d3950d79-5918-4ed7-bdad-dd04f22d3a30
BBBB=[0 1 0; 0 0 1];

# ╔═╡ 3895430e-9726-4c54-b0d2-9d7901f66f01
Calculate_Keq_stoichiometric_generic(AAAA,BBBB)

# ╔═╡ 3b18856d-27b8-4e1b-bdef-927f3e487e49
md"The $K$ matrix and the $X^A$ vector required for the equation for the general form of the Mass action are:

```math
K=\begin{pmatrix}
 k_1 & 0\\
 0 & k_2
\end{pmatrix} \tag{26}
```

```math
X^A=\begin{pmatrix}
X_1^1\cdot X_2^0 \cdot X_3^0\\
X_1^0\cdot X_2^1 \cdot X_3^0
\end{pmatrix} = \begin{pmatrix}
X_1 \\
 X_2
\end{pmatrix}  \tag{27}
```

so

```math
\begin{align}
 \begin{bmatrix}
\frac{\mathrm{d} X_1}{\mathrm{d} t}\\ \frac{\mathrm{d} X_2}{\mathrm{d} t} \\ \frac{\mathrm{d} X_3}{\mathrm{d} t} \end{bmatrix}= \begin{bmatrix}
-1 & 0\\
 1 & -1  \\
 0 & 1 \end{bmatrix} \begin{pmatrix}
 k_1 & 0\\
 0 & k_2
\end{pmatrix}\begin{pmatrix}
X_1 \\
 X_2
\end{pmatrix}
\end{align}
```

```math
\begin{align}
\begin{bmatrix}
\frac{\mathrm{d} X_1}{\mathrm{d} t}\\ \frac{\mathrm{d} X_2}{\mathrm{d} t} \\ \frac{\mathrm{d} X_3}{\mathrm{d} t} \end{bmatrix}=  \begin{bmatrix}
-1 & 0\\
 1 & -1  \\
 0 & 1 \end{bmatrix} \begin{pmatrix}
k_1 X_1 \\
 k_2 X_2
\end{pmatrix}   =   \begin{bmatrix}
 - k_1 \cdot X_1 \\
 k_1 \cdot X_1 - k_2 \cdot X_2  \\
 k_2 \cdot X_2 \end{bmatrix}
\end{align} \tag{28}
```

The differential equations of the system can be written using the differential form of Mass Action Law as:

```math
\begin{align}
\frac{dX_1}{dt} =  -k_1 \cdot X_1 \\
\frac{dX_2}{dt} =   k_1 \cdot X_1 \cdot X_2 - k_2 \cdot X_2 \tag{30}\\
\frac{dX_3}{dt} =   k_2 \cdot X_2
\end{align}
```
"

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
HypertextLiteral = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"

[compat]
HypertextLiteral = "~0.9.3"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.1"
manifest_format = "2.0"

[[deps.HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"
"""

# ╔═╡ Cell order:
# ╠═fffc6cac-c044-4cd9-a3b7-084341878faf
# ╟─dd3e39ac-8b2c-11ec-1774-7543d2573dc0
# ╠═7c7a696a-8716-4739-9a40-44a34c621485
# ╠═30a79e11-2165-4908-94a3-d020cc6d2e16
# ╠═e11754b1-6993-403b-a0ec-e0aaf86a5dd3
# ╟─ab0aa754-78e2-47b5-a398-a7e1ab8d6c58
# ╠═d2e26cb1-f136-4b96-8ef3-d3ee30f61b3d
# ╠═10f98ec7-fcad-4587-90f8-62f7b1b74f4d
# ╠═fee636f1-4944-46a4-bcbe-32b7273268c7
# ╠═56a68cf6-6e64-4875-b614-18987857bc06
# ╠═05faae5d-fdb6-4f1d-a77f-0a1337e1666d
# ╟─0d524884-afbf-4594-80ee-9fdfcdabbea0
# ╟─61a1a25f-21b1-40a6-aafa-17fddf4a848c
# ╟─2b58c174-c501-47ed-ade6-0ee9b88b76eb
# ╠═6c491ea3-10b9-46b3-a3ff-2c27d029efe0
# ╟─eb79f46c-5b6d-4b36-a6e0-f92fe6938b46
# ╟─3f13dd7b-5430-41b2-a6c2-e185f924ba6b
# ╟─c5a8887c-6875-481d-97bc-cb7da256d2fb
# ╟─fba4e0b3-52a5-4d1a-b7da-6740f2c2a836
# ╠═0731c83b-3b56-4b80-9654-1760dcbfca9b
# ╠═51ac8bbb-314d-4e82-ab3c-4d65be6f5e74
# ╠═7155b893-ec48-4213-b763-fe4d604efdaa
# ╠═e6a154b3-1f4f-4759-b0ab-62429b0c5200
# ╠═7206af1d-1d79-4c18-82c2-7979d3baa2d5
# ╟─b23e2b04-83a3-46f3-86e4-f49d8d42314d
# ╟─0ca58de1-28bd-4df1-88ca-efa8d17250f6
# ╠═52a42411-35b9-472a-bd7e-99e7e6395015
# ╠═637ae77a-a61c-4742-b737-aa4985ce6a32
# ╠═45f11cc5-59bd-44a6-977b-78e1bea3d9c5
# ╠═a21e1203-57a8-4677-968b-c0a90cec9403
# ╟─a796442b-85fa-45af-8c90-9ce9a8d4fcd2
# ╟─3ac754bf-e6de-417f-9093-35a73f036a1e
# ╠═6c027717-d272-4a3d-824c-54ea64f4fd92
# ╠═d3950d79-5918-4ed7-bdad-dd04f22d3a30
# ╠═3895430e-9726-4c54-b0d2-9d7901f66f01
# ╟─3b18856d-27b8-4e1b-bdef-927f3e487e49
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_path: 7 General Form of Mass Action.jl

max_stars_repo_name: davidgmiguez/PlutoNotebooks

max_issues_repo_path: 7 General Form of Mass Action.jl

max_issues_repo_name: davidgmiguez/PlutoNotebooks

max_forks_repo_path: 7 General Form of Mass Action.jl

max_forks_repo_name: davidgmiguez/PlutoNotebooks


text:
include("../optimizers/linesearch.jl")
include("../utils/functions.jl")
include("../utils/plot.jl")

using Gadfly

# general function to plot the output of the algorithms
function plot_values(fx1, fx2, fx3, fx4, xlab, ylab, title)
    nsamps = length(fx1)
    nsamps2 = length(fx2)
    nsamps3 = length(fx3)
    nsamps4 = length(fx4)

    Gadfly.plot(layer(x=1:nsamps, y=fx1, Geom.line,
            Theme(default_color=color("blue"))),
        layer(x=1:nsamps2, y=fx2, Geom.line,
            Theme(default_color=color("red"))),
        layer(x=1:nsamps3, y=fx3, Geom.line,
            Theme(default_color=color("orange"))),
        layer(x=1:nsamps4, y=fx4, Geom.line,
            Theme(default_color=color("purple"))),
        Guide.xlabel(xlab), Guide.ylabel(ylab),
        Guide.title(title),
        Guide.manual_color_key("Legend",
            ["Newton", "steepest", "Newton-CG", "BFGS"],
            ["blue", "red", "orange", "purple"]),
        Scale.x_log10, Scale.y_log10)
end

@time xvals = line_search(fenton, [3.;4.],
    fenton_g, fenton_h, "newton", 1000);
@time svals = line_search(fenton, [3.;4.],
    fenton_g, fenton_h, "steepest", 1000);
@time cvals = line_search(fenton, [3.;4.],
    fenton_g, fenton_h, "newton_CG", 1000);
@time qvals = line_search(fenton, [3.;4.],
    fenton_g, fenton_h, "BFGS", 1000);

nsamps = length(xvals)
nsamps2 = length(svals)
nsamps3 = length(cvals)
nsamps4 = length(qvals)

fx = [fenton(xvals[i]) for i in 1:nsamps]
fx2 = [fenton(svals[i]) for i in 1:nsamps2]
fx3 = [fenton(cvals[i]) for i in 1:nsamps3]
fx4 = [fenton(qvals[i]) for i in 1:nsamps4]

plot_values(fx, fx2, fx3, fx4, "iteration", "f(x)", "Value of function")

grads = [norm(fenton_g(xvals[i]), 2) for i in 1:nsamps]
grads2 = [norm(fenton_g(svals[i]), 2) for i in 1:nsamps2]
grads3 = [norm(fenton_g(cvals[i]), 2) for i in 1:nsamps3]
grads4 = [norm(fenton_g(qvals[i]), 2) for i in 1:nsamps4]

r1 = compute_grad_ratio(grads)
r2 = compute_grad_ratio(grads2)
r3 = compute_grad_ratio(grads3)
r4 = compute_grad_ratio(grads4)

plot_values(r1, r2, r3, r4, "iteration",
    "gradient norm ratios", "gradient norm ratios")

plot_values(grads, grads2, grads3, grads4, "iteration",
    "gradient norms", "gradient norms")

@time xvals = line_search(rosenbrock, ones(100)*3, rosenbrock_g, rosenbrock_h,
    "newton", 2000, 1e-8);
@time svals = line_search(rosenbrock, ones(100)*3, rosenbrock_g, rosenbrock_h,
    "steepest", 5000, 1e-8);
@time cvals = line_search(rosenbrock, ones(100)*3, rosenbrock_g, rosenbrock_h,
    "newton_CG", 2000, 1e-8);
@time qvals = line_search(rosenbrock, ones(100)*3, rosenbrock_g, rosenbrock_h,
    "BFGS", 2000, 1e-8);

func = rosenbrock
func_g = rosenbrock_g

nsamps = length(xvals)
nsamps2 = length(svals)
nsamps3 = length(cvals)
nsamps4 = length(qvals)

fx = [func(xvals[i]) for i in 1:nsamps]
fx2 = [func(svals[i]) for i in 1:nsamps2]
fx3 = [func(cvals[i]) for i in 1:nsamps3]
fx4 = [func(qvals[i]) for i in 1:nsamps4]

plot_values(fx, fx2, fx3, fx4, "iteration", "f(x)", "Value of function")

nsamps = length(xvals)

grads = [norm(func_g(xvals[i]), 2) for i in 1:nsamps]
grads2 = [norm(func_g(svals[i]), 2) for i in 1:nsamps2]
grads3 = [norm(func_g(cvals[i]), 2) for i in 1:nsamps3]
grads4 = [norm(func_g(qvals[i]), 2) for i in 1:nsamps4]

r1 = compute_grad_ratio(grads)
r2 = compute_grad_ratio(grads2)
r3 = compute_grad_ratio(grads3)
r4 = compute_grad_ratio(grads4)

plot_values(r1, r2, r3, r4, "iteration", "gradient norm ratios",
    "gradient norm ratios")

plot_values(grads, grads2, grads3, grads4, "iteration", "gradient norms",
    "gradient norms")

@time xvals = line_search(cute, ones(100)*10, cute_g, cute_h, "newton",
    2000, 1e-8);
@time svals = line_search(cute, ones(100)*10, cute_g, cute_h, "steepest",
    2000, 1e-8);
@time cvals = line_search(cute, ones(100)*10, cute_g, cute_h, "newton_CG",
    2000, 1e-8);
@time qvals = line_search(cute, ones(100)*10, cute_g, cute_h, "BFGS",
    2000, 1e-8);

func = cute
func_g = cute_g

nsamps = length(xvals)
nsamps2 = length(svals)
nsamps3 = length(cvals)
nsamps4 = length(qvals)

fx = [func(xvals[i]) for i in 1:nsamps]
fx2 = [func(svals[i]) for i in 1:nsamps2]
fx3 = [func(cvals[i]) for i in 1:nsamps3]
fx4 = [func(qvals[i]) for i in 1:nsamps4]

plot_values(fx, fx2, fx3, fx4, "iteration", "f(x)", "Value of function")

grads = [norm(func_g(xvals[i]), 2) for i in 1:nsamps]
grads2 = [norm(func_g(svals[i]), 2) for i in 1:nsamps2]
grads3 = [norm(func(cvals[i]), 2) for i in 1:nsamps3]
grads4 = [norm(func_g(qvals[i]), 2) for i in 1:nsamps4]

r1 = compute_grad_ratio(grads)
r2 = compute_grad_ratio(grads2)
r3 = compute_grad_ratio(grads3)
r4 = compute_grad_ratio(grads4)

plot_values(r1, r2, r3, r4, "iteration", "gradient norm ratio",
    "gradient norm ratios")

plot_values(grads, grads2, grads3, grads4, "iteration", "gradient norm",
    "gradient norms")


text:
@testset "009" begin
    @test tdpe.e009(3+4+5) == 3*4*5
end

max_stars_repo_name: mithodin/ttpe

max_issues_repo_name: mithodin/ttpe

max_forks_repo_name: mithodin/ttpe


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMPlates.jl/blob/master/LICENSE

using Test, FEMPlates

@testset "FEMPlates.jl" begin
    @testset "test_mindlin_stiffness_matrix" begin
        include("test_mindlin_stiffness_matrix.jl")
    end
    @testset "test_mindlin_force_vector" begin
        include("test_mindlin_force_vector.jl")
    end
    @testset "test_dkt" begin
        include("test_dkt.jl")
    end
end


text:
"""
    Current assumption is that the pairs are passed in so the first qubit is on top of or
    or to the left of the second qubit

    Designed to be called by hinton plots and plot the qubit to qubit interactions.
"""
function outlineQubits(qbitPairs,ax)
    colors = ["darkslategray","slategray"]
    for x in qbitPairs
        firstP = translate_14Q_toIndex(x[1])
        secondP = translate_14Q_toIndex(x[2])
        if firstP[1]<secondP[1] # To the left
            rect = plt.Rectangle([firstP[1]+0.55,firstP[2]-0.45], 1.9, 0.9, facecolor="slategray",edgecolor="darkslategray",linewidth=2)
        else # To the bottom
            rect = plt.Rectangle([firstP[1]+0.55,firstP[2]-0.45], 0.9, 1.9, facecolor="slategray",edgecolor="darkslategray",linewidth=2)
        end
        ax.add_patch(rect)
        #rect = []
        #if x[1] < 6
        #    rect = plt[:Rectangle]([x[1]+0.55,0.55], 1.9, 0.9, facecolor="slategray",edgecolor="darkslategray",linewidth=2)
        #else
        #    rect = plt[:Rectangle]([(13-x[1])+1.55,1.55], 1.9, 0.9, facecolor="slategray",edgecolor="darkslategray",linewidth=2)
        #end
        #ax[:add_patch](rect)
    end
end

"""
Adapted from demo algorithm found at
https://matplotlib.org/gallery/specialty_plots/hinton_demo.html

No default on ax or max_weight - pass it in (allows different plots to keep it constant)

Modified to accept keyword, which allows for a fixed highlight of negative values.

## Arguments:
-   `matrix: Array{Float64,1}` of values to be plotted

-   `max_weight: Float` The value required to 'fill a box'

-   `ax: PyCall.PyObject`, the graphical axis (get from gca())

## Optional Arguments:

-   `addScale: Bool` If true will draw the scale of the box (from max_weight) to the rhs of the plot.
-   `highlightNegative: Bool` If true will print negative values as a qubit, value given by qubit and font by fontsize.
-   `fontsize: Float64` default is 12.
-   `qubit: Int` qubit to highlight i.e. the qubit wrt which all the other information is printed.

"""
function hinton(matrix,max_weight, ax;
                        errorsH=[],
                        errorsL=[],
                        qubitPairs=[],
                        highlightNegative = false,
                        fontsize = 12,
                        qubit=0,
                        addScale=false,
                        xoffset=0,
                        yoffset=0,
                        joined=[])
    #"""Draw Hinton diagram for visualizing a weight matrix."""
    if errorsH == []
        errorsH = deepcopy(matrix)
    end
    if errorsL == []
        errorsL = deepcopy(matrix)
    end
    matrix=transpose(matrix)
    errorsH=transpose(errorsH)
    errorsL = transpose(errorsL)
    ax.patch.set_facecolor("gray")
    ax.set_aspect("equal", "box")
    ax.xaxis.set_major_locator(plt.NullLocator())
    ax.yaxis.set_major_locator(plt.NullLocator())
    for x in 1:size(matrix)[1]
        for y in 1:size(matrix)[2]
            color="lightgray"
            size = 1.0
            if matrix[x,y]!= -3
                rect = plt.Rectangle([x - size / 2, y - size / 2], size, size,
                             facecolor="gray", edgecolor="lightgray")
                ax.add_patch(rect)
            end
        end
    end
    if addScale
        ax.add_patch(rect)
        rect = plt.Rectangle([size(matrix)[1]+xoffset+0.5,yoffset+1.3],0.1,0.1,
                facecolor = "gray")

        ax.add_patch(rect)

        plt.annotate(s="", xy=(size(matrix)[1]-1+0.7,1.5), xytext=(size(matrix)[1]-1+0.7,0.5), arrowprops=Dict([(:arrowstyle,"<->")]))
        plt.annotate( string(max_weight),
                    xy=(size(matrix)[1]-1+0.8,0.8),
                    xytext=(size(matrix)[1]-1+0.8,0.8),
                    va = "top", ha="left")

    end

    joined=[]
    for x in 1:size(matrix)[1]
        for y in 1:size(matrix)[2]
            w= matrix[x,y]
            if w == -1
                 push!(joined,(x,y))
            end
        end
    end

    outlineQubits(qubitPairs,ax)
    for x in 1:size(matrix)[1]
        for y in 1:size(matrix)[2]
            w= matrix[x,y]

            barsH = errorsH[x,y]
            barsL = errorsL[x,y]

            size = min(sqrt(abs(w) / max_weight),1)
            largestError = 0
            smallestError = 0
            if w < 0
                largestError = min(sqrt(abs(barsL)/max_weight),1)
                smallestError = min(sqrt(abs(barsH)/max_weight),1)
                if barsH > 0
                    smallestError = 0
                end
            else
                largestError = min(sqrt(abs(barsH)/max_weight),1)
                smallestError = min(sqrt(abs(barsL)/max_weight),1)
                if barsL < 0
                    smallestError = 0
                end
            end
            if largestError < size
                    #print("CHECK ($x,$y): largestError = $largestError, size = $size\n")
                largestError = size
            end
            if smallestError > size
                    #print("CHECK ($x,$y): smallestError = $smallestError, size = $size\n")
                smallestError = size
            end
            if w >= 0

                color = "lightgray"
                rect = plt.Rectangle([x - largestError / 2, y - largestError / 2], largestError, largestError,
                                 facecolor=color, edgecolor= color)
                ax.add_patch(rect)

                rect = plt.Rectangle([x - size / 2, y - size / 2], size, size,
                            facecolor=color, edgecolor= "black",linestyle=":")
                ax.add_patch(rect)
                color = "white"
                rect = plt.Rectangle([x - smallestError / 2, y - smallestError / 2], smallestError, smallestError,
                                            facecolor="white", edgecolor= "white")
                ax.add_patch(rect)

            else
                if highlightNegative || w == -3 || w == -1
                    if w == -3
            #            rect = plt[:Rectangle]([x - size / 2, y - size / 2], size, size,
            #            facecolor="gray", edgecolor="gray")
            #            ax[:add_patch](rect)

                    elseif w == -2
                        color = "red"
                        plt.text(x,y,"X",color="red",fontsize=fontsize,
                                horizontalalignment="center",
                                verticalalignment="center"
                                )
                    else

                        if length(joined) == 2 # need to generalise if we are putting in more than one 'pair'
                            (a,b),(c,d) = joined
                            if (b==d) # join on same plane
                                rect = plt.Rectangle([a ,b - 0.35],1,0.62,facecolor="blue",edgecolor="black")
                                ax.add_patch(rect)
                            else
                                rect = plt.Rectangle([a-0.33,b],0.62,1,facecolor="blue",edgecolor="black")
                                ax.add_patch(rect)
                            end
                        end

                        color = "blue"
                        size=0.7
                        circ = plt.Circle([x - 0.01  , y - 0.03  ], size/2,
                             facecolor=color, edgecolor= "black")
                        ax.add_patch(circ)
                        q=99
                        if y==1
                            q = x-1
                        else
                            q = (15-x)
                        end
                        plt.text(x,y,"Q"*string(q),color="white",fontsize=fontsize,
                                horizontalalignment="center",
                                verticalalignment="center"
                                )
                        push!(joined,(x,y))
                    end
                else # just a plain old negative number
                    color = "darkslategrey"
                    rect = plt.Rectangle([x - largestError / 2, y - largestError / 2], largestError, largestError,
                                 facecolor=color, edgecolor= color)
                    ax.add_patch(rect)

                    rect = plt.Rectangle([x - size / 2, y - size / 2], size, size,
                            facecolor=color, edgecolor= "white",linestyle=":")
                    ax.add_patch(rect)

                    color = "black"
                    rect = plt.Rectangle([x - smallestError / 2, y - smallestError / 2], smallestError, smallestError,
                                            facecolor="black", edgecolor= "black")
                    ax.add_patch(rect)
                end
            end

        end
    end
    ax.autoscale_view()
    ax.invert_yaxis()

end

"""
    translateLocationd(data,qubit)

Rather specific IBMQX16 function. Basically remaps the qubits into the same order as they appear in the IBM 16 Qubit machine.
"""
function translateLocation(data,i)
    a = Array{Float64}(undef,2,8)
    a[2,1]=data[1]
    a[1,1]=data[2]
    a[1,2:8]= data[3:9]
    a[2,2:8]=reverse(data[10:16])

    return a
end

"""
    translateLocationd(data,qubit)

Rather specific IBMQX14 function. Basically remaps the qubits into the same order as they appear in the IBM 14 Qubit machine.
"""
function translate_14Q_Location(data,i)
    a = Array{Float64}(undef,2,8)
    a[1,1:7]= data[1:7]
    a[2,2:8]=reverse(data[8:14])
    a[2,1]=-3
    a[1,8]=-3
    return a
end

"""
Specific to a particular device. Translates the qubit to an x,y grid reference.
"""
function translate_14Q_toIndex(i)
   if i < 7
        return (i,1)
   end
   return (14-i,2)
end

"""
    drawIBMHinton(titleText)

Another specific function (specific to Melbourne). Basically draws hinton plots for each of the qubits, given the mutualInformation Ps.
If you specify a save file it will save it.

"""
function drawIBM_14Q_Hinton(mutualPs,titleText;saveFile = "",titleFontSize=24)
    addScale = false
    fig = figure("Slightly larger",figsize=(18,6))
    scale = round(maximum(maximum.(mutualPs)/5),RoundUp,digits=3)
    fig[:suptitle](titleText*" Full box = $scale",fontsize=titleFontSize)
    ax = gca()
    ax[:set_facecolor]("gray")
    for (ix,i) in enumerate(1:4)
      subplot(4,4,ix)
      ax=gca()
     hinton(translate_14Q_Location(mutualPs[i],i),scale,ax,highlightNegative=true,qubit=i-1)
    end

    for (ix,i) in enumerate(5:7)
        subplot(4,4,4+ix)
        ax=gca()
        hinton(translate_14Q_Location(mutualPs[i]',i),scale,ax,highlightNegative=true,qubit=i-1)
    end
    for (ix,i) in enumerate(14:-1:12)
        subplot(4,4,9+ix)
        ax=gca()
        hinton(translate_14Q_Location(mutualPs[i]',i),scale,ax,highlightNegative=true,qubit=i-1,
            fontsize=i > 10 ? 9 : 12)
    end
    for (ix,i) in enumerate(11:-1:8)
        subplot(4,4,12+ix)
        ax=gca()
        hinton(translate_14Q_Location(mutualPs[i]',i),scale,ax,highlightNegative=true,qubit=i-1,
            fontsize=i >10 ? 9 : 12)
    end
    fig[:subplots_adjust](top=1.3)
    plt[:tight_layout]()
    if saveFile != ""
        savefig(saveFile)
    end
end

"""
    drawIBMHinton(titleText)

Another specific function (specific to IBMQX16). Basically draws hinton plots for each of the qubits, given the mutualInformation Ps.
If you specify a save file it will save it.

"""
function drawIBMHinton(mutualPs,titleText;saveFile = "",titleFontSize=24)
    addScale = false
    fig = figure("Slightly larger",figsize=(18,6))
    scale = round(maximum(maximum.(mutualPs)),RoundUp,digits=3)
    fig[:suptitle](titleText*" Full box = $scale",fontsize=titleFontSize)
    ax = gca()
    ax[:set_facecolor]("gray")
    for (ix,i) in enumerate(2:4)
      subplot(4,4,ix)
      ax=gca()
     hinton(translateLocation(mutualPs[i]',i),scale,ax,highlightNegative=true,qubit=i-1)
    end

    for (ix,i) in enumerate(5:5)
     subplot(4,4,4)
      ax=gca()
      hinton(translateLocation(mutualPs[i]',i),scale,ax,highlightNegative=true,qubit=i-1,addScale=addScale)
    end
    for (ix,i) in enumerate(6:9)
        subplot(4,4,4+ix)
        ax=gca()
        hinton(translateLocation(mutualPs[i]',i),scale,ax,highlightNegative=true,qubit=i-1)
    end
    for (ix,i) in enumerate(vcat([1],16:-1:14))
        subplot(4,4,8+ix)
        ax=gca()
        hinton(translateLocation(mutualPs[i]',i),scale,ax,highlightNegative=true,qubit=i-1,
            fontsize=i > 10 ? 9 : 12)
    end
    for (ix,i) in enumerate(13:-1:10)
        subplot(4,4,12+ix)
        ax=gca()
        hinton(translateLocation(mutualPs[i]',i),scale,ax,highlightNegative=true,qubit=i-1,
            fontsize=i >10 ? 9 : 12)
    end
    fig[:subplots_adjust](top=1.3)
    plt[:tight_layout]()
    if saveFile != ""
        savefig(saveFile)
    end
end

"""
    Draw Hinton diagram for visualizing a weight matrix.

    Adapted from demo algorithm found at
    https://matplotlib.org/gallery/specialty_plots/hinton_demo.html

    No default on ax or max_weight - pass it in (allows different plots to keep it constant)

    We appear to have suffered from parameter creep! There are a large number of adjustable parameters now, all with semi-sensible defaults

## Arguments:
-   `matrix: Array{Float64,1}` of values to be plotted
-   `max_weight: Float` The value required to 'fill a box'
-   `ax`, the graphical axis (get from gca())
-   `highCorr/lowCorr : Array{Float64,1}`- possible high/low values. Used for drawing error boxes.
-   `highlightNegative`: If true will print negative values as a qubit, value given by qubit and font by fontsize.
-   `addScale`: If true will draw the scale of the box (from max_weight) to the rhs of the plot.
-    `fontsize = 12`
-    `qubit=0`
-    `addScale=false`
-    `addAxis=false`
-    `addOneAxis=false`
-    `adjust = 0.5`
-    `noQ=false`
-    `ind=[]`
-    `showQubits=false`
-    `qlabels=[]`
-    `stagger=true`
-    `voffset=0`
-    `useBars=false`
-    `thickness=0.01`
-    `altColor = "gray"`

"""
function covhinton(matrix, max_weight, ax;
    highCorr= [],
    lowCorr=[],
    highlightNegative = false,
    fontsize = 12,
    qubit=0,
    addScale=false,
    addAxis=false,
    addOneAxis=false,
    adjust = 0.5,
    noQ=false,
    ind=[],
    showQubits=false,
    qlabels=[],
    stagger=true,
    voffset=0,
    useBars=false,
    thickness=0.01,
    altColor = "gray")
    matrix=transpose(matrix)
    if highCorr ==[]
        highCorr = matrix
    else
        highCorr=transpose(highCorr)
    end
    if lowCorr==[]
        lowCorr = matrix
    else
        lowCorr = transpose(lowCorr)
    end
    maxX = size(matrix)[1]
    maxY = size(matrix)[2]
    xoffset=10
    yoffset=10
    ax.axis("off")

    ax.patch.set_facecolor("lightgray")
    ax.set_aspect("equal", "box")
    ax.xaxis.set_major_locator(plt.NullLocator())
    ax.yaxis.set_major_locator(plt.NullLocator())
    #xlim([0,size(matrix[1])+2])
    xlim([0,maxX+3])
    ylim([0,maxY+3])
    rect = plt.Rectangle([1.45,1.45],maxX+0.1,maxY+0.1,facecolor="black",edgecolor="black")
    ax.add_patch(rect)

    for x in 1:size(matrix)[1]
        for y in 1:size(matrix)[2]
            if x > y
                color = altColor
            else
                color="gray"
            end
            size = 1.0
            rect = plt.Rectangle([x+1 - size / 2, y+1 - size / 2], size, size,
               facecolor=color, edgecolor="lightgray")
            ax.add_patch(rect)
        end
    end
    if addAxis
        if ind==[]
            ind = collect(0:13)
        end
        qstring = "Q"
        if noQ
            qstring=""
        end
        if qlabels == []
            qlabels=["$qstring$(ind[x])" for x =1:maxX]
        end
        for x in 1:maxX
            plt.annotate( "$(qlabels[x])",
                xy=(2,0.2),
                xytext=(x+adjust,0.4+(stagger ? x%2 : 0)*0.4+voffset),
                va = "top", ha="left",fontsize=fontsize)
            plt.annotate( "$(qlabels[x])",
                xy=(1.5,0.2),
                xytext=(adjust-0.2,x+adjust+0.2),
                va = "top", ha="left",fontsize=fontsize)
            if !(addOneAxis)
                plt.annotate( "$(qlabels[x])",
                    xy=(2,0.2),
                    xytext=(maxX+1.7+adjust,x+adjust+0.2),
                    va = "top", ha="left",fontsize=fontsize)
                plt.annotate( "$(qlabels[x])",
                    xy=(2,0.2),
                    xytext=(x+adjust,maxY+2+(stagger ? x%2 : 0)*0.4+voffset),
                    va = "top", ha="left",fontsize=fontsize)
            end
        end

    end
    if addScale
        ax.add_patch(rect)
        rect = plt.Rectangle([size(matrix)[1]+xoffset+0.5,yoffset+1.3],0.1,0.1,
            facecolor = "blue")

        ax.add_patch(rect)

        plt.annotate(s="", xy=(size(matrix)[1]+1.7,2.5), xytext=(size(matrix)[1]+1.7,1.5), arrowprops=Dict([(:arrowstyle,"<->")]))
        plt.annotate( string(max_weight),
            xy=(size(matrix)[1]+1.8,1.8),
            xytext=(size(matrix)[1]+1.8,1.8), fontsize=fontsize,
            va = "top", ha="left")

    end
    for x in 1:size(matrix)[1]
        for y in 1:size(matrix)[2]
            w= matrix[x,y]
            if isapprox(w,1) && showQubits # assume its the qubit
                 color = "blue"
                size=0.7
                circ = plt.Circle([x+1 - 0.01  , y+1 - 0.03  ], size/2,
                   facecolor=color, edgecolor= "black")
                ax.add_patch(circ)
                q=x-1
                plt.text(x+1,y+1,"Q"*string(q),color="white",fontsize=fontsize,
                    horizontalalignment="center",
                    verticalalignment="center"
                   )

            else
                size = min(sqrt(abs(w) / max_weight),1)
                maxw = highCorr[x,y]
                minw = lowCorr[x,y]
                largestError = 0
                smallestError = 0
                if w < 0
                    largestError = min(sqrt(abs(minw)/max_weight),1)
                    smallestError = min(sqrt(abs(maxw)/max_weight),1)
                    if maxw > 0
                        if useBars
                            smallestError = -1*smallestError
                        else
                            smallestError = 0
                        end
                    end
                else
                    largestError = min(sqrt(abs(maxw)/max_weight),1)
                    smallestError = min(sqrt(abs(minw)/max_weight),1)
                    if minw < 0
                        if useBars
                            smallestError = -1*smallestError
                        else
                            smallestError = 0
                        end
                    end
                end
                if largestError < size
                    #print("CHECK ($x,$y): largestError = $largestError, size = $size\n")
                    largestError = size
                end
                if smallestError > size
                    #print("CHECK ($x,$y): smallestError = $smallestError, size = $size\n")
                    smallestError = size
                end

                if w >= 0
                    if !useBars
                        # large square light gray
                        # actual square light gray, black border
                        # small sqaure (if !=0) white

                        color = "lightgray"
                        rect = plt.Rectangle([x+1 - largestError / 2, y+1 - largestError / 2], largestError, largestError,
                                 facecolor=color, edgecolor= color)
                        ax.add_patch(rect)

                        rect = plt.Rectangle([x+1 - size / 2, y+1 - size / 2], size, size,
                        facecolor=color, edgecolor= "black",linestyle=":")
                        ax.add_patch(rect)
                        color = "white"
                        rect = plt.Rectangle([x+1 - smallestError / 2, y+1 - smallestError / 2], smallestError, smallestError,
                                            facecolor="white", edgecolor= "white")
                        ax.add_patch(rect)
                    else
                        color = "white"
                        rect = plt.Rectangle([x+1 - size / 2, y+1 - size / 2], size, size,
                        facecolor=color, edgecolor= color)
                        ax.add_patch(rect)
                        # plot the error bar
                        if x!=y # not on the diagonals
                            color = "red"
                            moresize = 0
                            if smallestError < 0
                               moresize = size/2
                            end
                            rect = plt.Rectangle([x+1 + size/3 , y+1- smallestError/2],  thickness,-moresize-largestError/2+smallestError/2,
                                facecolor=color, edgecolor=color)
                            ax.add_patch(rect)
                #print("($x,$y): Size is $(size/2), going from $(-smallestError/2) to $(-moresize-largestError/2+smallestError/2)\n")
                        end
                    end

                else
                    if highlightNegative
                        color = "blue"
                        size=0.7
                        circ = plt.Circle([x+1 - 0.01  , y+1 - 0.03  ], size/2,
                        facecolor=color, edgecolor= "black")
                        ax.add_patch(circ)
                        if x==y
                            plt.text(x+1,y+1,"Q"*string(y-1),color="white",fontsize=fontsize,
                             horizontalalignment="center",
                             verticalalignment="center"
                            )
                        end
                        rect = plt.Rectangle([x+0.5- 0.01, y+0.5 - 0.01], 0.98,0.98,
                                 facecolor="slategrey", edgecolor= "darkslategrey")
                            ax.add_patch(rect)

                    else
                         if !useBars
                        # large square dim gray
                        # actual square dim gray, white border
                        # small sqaure (if !=0) white

                            color = "darkslategrey"
                            rect = plt.Rectangle([x+1 - largestError / 2, y+1 - largestError / 2], largestError, largestError,
                                 facecolor=color, edgecolor= color)
                            ax.add_patch(rect)

                            rect = plt.Rectangle([x+1 - size / 2, y+1 - size / 2], size, size,
                            facecolor=color, edgecolor= "white",linestyle=":")
                            ax.add_patch(rect)
                            color = "black"
                            rect = plt.Rectangle([x+1 - smallestError / 2, y+1 - smallestError / 2], smallestError, smallestError,
                                            facecolor="black", edgecolor= "black")
                            ax.add_patch(rect)
                        else
                            color = "black"
                            rect = plt.Rectangle([x+1 - size / 2, y+1 - size / 2], size, size,
                            facecolor=color, edgecolor= color)
                            ax.add_patch(rect)
                            # plot the error bar
                            color = "red"
                            moresize = 0
                            if smallestError < 0
                               moresize = size/2
                            end
                            rect = plt.Rectangle([x+1 + size/3 , y+1- smallestError/2],  thickness,-moresize-largestError/2+smallestError/2,
                                facecolor=color, edgecolor=color)
                            ax.add_patch(rect)
                        end
                     end

                end
            end
        end #for y
    end # for x
    ax.autoscale_view()
    ax.invert_yaxis()

end