text:
module TestCSets
using Test

using Catlab, Catlab.Theories, Catlab.Graphs, Catlab.CategoricalAlgebra,
  Catlab.CategoricalAlgebra.FinSets
using Catlab.Graphs.BasicGraphs: TheoryGraph
using Catlab.Present

@present TheoryDDS(FreeSchema) begin
  X::Ob
  Φ::Hom(X,X)
end
const DDS = CSetType(TheoryDDS, index=[:Φ])

# FinSets interop
#################

g = Graph(6)
add_edges!(g, 2:4, 3:5)
f = FinFunction(g, :V)
@test collect(f) == 1:6
@test is_indexed(f)
f = FinFunction(g, :src)
@test codom(f) == FinSet(6)
@test collect(f) == 2:4
@test is_indexed(f)

f = FinDomFunction(g, :E)
@test collect(f) == 1:3
@test is_indexed(f)
f = FinDomFunction(g, :tgt)
@test codom(f) == TypeSet(Int)
@test collect(f) == 3:5
@test is_indexed(f)

g = path_graph(WeightedGraph{Float64}, 3, E=(weight=[0.5, 1.5],))
f = FinDomFunction(g, :weight)
@test codom(f) == TypeSet(Float64)
@test collect(f) == [0.5, 1.5]

# C-set morphisms
#################

# Constructors and accessors.
g, h = path_graph(Graph, 4), cycle_graph(Graph, 2)
α = CSetTransformation((V=[1,2,1,2], E=[1,2,1]), g, h)
@test components(α) == (V=α[:V], E=α[:E])
@test α[:V] isa FinFunction{Int} && α[:E] isa FinFunction{Int}
@test α[:V](3) == 1
@test α[:E](2) == 2

α′ = CSetTransformation(g, h, V=[1,2,1,2], E=[1,2,1])
@test components(α′) == components(α)
α′′ = CSetTransformation(g, h, V=FinFunction([1,2,1,2]), E=FinFunction([1,2,1]))
@test components(α′′) == components(α)

# Naturality.
@test is_natural(α)
β = CSetTransformation((V=[1,2,1,2], E=[1,1,1]), g, h)
@test !is_natural(β)
β = CSetTransformation((V=[2,1], E=[2,1]), h, h)
@test is_natural(β)

# Category of C-sets.
@test dom(α) === g
@test codom(α) === h
@test compose(α,β) == CSetTransformation((V=α[:V]⋅β[:V], E=α[:E]⋅β[:E]), g, h)
@test force(compose(id(g), α)) == α
@test force(compose(α, id(h))) == α

# Limits
#-------

# Terminal object in Graph: the self-loop.
term = Graph(1)
add_edge!(term, 1, 1)
lim = terminal(Graph)
@test ob(lim) == term
@test force(delete(lim, g)) ==
  CSetTransformation((V=fill(1,4), E=fill(1,3)), g, term)

# Products in Graph: unitality.
lim = product(g, term)
@test ob(lim) == g
@test force(proj1(lim)) == force(id(g))
@test force(proj2(lim)) ==
  CSetTransformation((V=fill(1,4), E=fill(1,3)), g, term)

# Product in Graph: two directed intervals (Reyes et al 2004, p. 48).
I = Graph(2)
add_edge!(I, 1, 2)
prod = ob(product(I, I))
@test nv(prod) == 4
@test ne(prod) == 1
@test src(prod) != tgt(prod)

# Product in Graph: deleting edges by multiplying by an isolated vertex.
g = path_graph(Graph, 4)
g0 = ob(product(g, Graph(1)))
@test nv(g0) == nv(g)
@test ne(g0) == 0

# Product in Graph: copying edges by multiplying by the double self-loop.
loop2 = Graph(1)
add_edges!(loop2, [1,1], [1,1])
lim = product(g, loop2)
g2 = ob(lim)
@test nv(g2) == nv(g)
@test ne(g2) == 2*ne(g)
@test src(g2) == repeat(src(g), 2)
@test tgt(g2) == repeat(tgt(g), 2)
α = CSetTransformation((V=[2,3], E=[2]), I, g)
β = CSetTransformation((V=[1,1], E=[2]), I, loop2)
γ = pair(lim, α, β)
@test force(γ⋅proj1(lim)) == α
@test force(γ⋅proj2(lim)) == β

# Equalizer in Graph from (Reyes et al 2004, p. 50).
g, h = cycle_graph(Graph, 2), Graph(2)
add_edges!(h, [1,2,2], [2,1,1])
ϕ = CSetTransformation((V=[1,2], E=[1,2]), g, h)
ψ = CSetTransformation((V=[1,2], E=[1,3]), g, h)
@test is_natural(ϕ) && is_natural(ψ)
eq = equalizer(ϕ, ψ)
@test ob(eq) == I
@test force(incl(eq)[:V]) == FinFunction([1,2])
@test force(incl(eq)[:E]) == FinFunction([1], 2)

# Pullback in Graph from (Reyes et al 2004, p. 53).
# This test exposes an error in the text: using their notation, there should be
# a second edge between the vertices (1,3) and (1,4).
g0, g1, g2 = Graph(2), Graph(3), Graph(2)
add_edges!(g0, [1,1,2], [1,2,2])
add_edges!(g1, [1,2,3], [2,3,3])
add_edges!(g2, [1,2,2], [1,2,2])
ϕ = CSetTransformation((V=[1,2,2], E=[2,3,3]), g1, g0)
ψ = CSetTransformation((V=[1,2], E=[1,3,3]), g2, g0)
@test is_natural(ϕ) && is_natural(ψ)
lim = pullback(ϕ, ψ)
@test nv(ob(lim)) == 3
@test sort!(collect(zip(src(ob(lim)), tgt(ob(lim))))) ==
  [(2,3), (2,3), (3,3), (3,3)]
@test is_natural(proj1(lim)) && is_natural(proj2(lim))

# Same pullback using generic limit interface.
diagram = FreeDiagram([g1, g2, g0], [(ϕ,1,3), (ψ,2,3)])
π1, π2, _ = lim′ = limit(diagram)
@test ob(lim′) == ob(lim)
@test force(π1) == force(proj1(lim))
@test force(π2) == force(proj2(lim))

# Colimits
#---------

# Initial object in graph: the empty graph.
colim = initial(Graph)
@test ob(colim) == Graph()
@test create(colim, g) == CSetTransformation((V=Int[], E=Int[]), Graph(), g)

# Coproducts in Graph: unitality.
g = path_graph(Graph, 4)
colim = coproduct(g, Graph())
@test ob(colim) == g
@test force(coproj1(colim)) == force(id(g))
@test force(coproj2(colim)) ==
  CSetTransformation((V=Int[], E=Int[]), Graph(), g)

# Coproduct in Graph.
h = cycle_graph(Graph, 2)
colim = coproduct(g, h)
coprod = ob(colim)
@test nv(coprod) == 6
@test src(coprod) == [1,2,3,5,6]
@test tgt(coprod) == [2,3,4,6,5]
α = CSetTransformation((V=[1,2,1,2], E=[1,2,1]), g, h)
β = id(h)
γ = copair(colim, α, β)
@test force(coproj1(colim)⋅γ) == α
@test force(coproj2(colim)⋅γ) == force(β)

# Coequalizer in Graph: collapsing a segment to a loop.
g = Graph(2)
add_edge!(g, 1, 2)
α = CSetTransformation((V=[1], E=Int[]), Graph(1), g)
β = CSetTransformation((V=[2], E=Int[]), Graph(1), g)
@test is_natural(α) && is_natural(β)
coeq = coequalizer(α, β)
@test ob(coeq) == ob(terminal(Graph))
@test force(proj(coeq)[:V]) == FinFunction([1,1])
@test force(proj(coeq)[:E]) == FinFunction([1])

# Pushout in Graph from (Reyes et al 2004, p. 59).
α = CSetTransformation((V=[2], E=Int[]), Graph(1), g)
β = CSetTransformation((V=[1], E=Int[]), Graph(1), ob(terminal(Graph)))
@test is_natural(α) && is_natural(β)
colim = pushout(α, β)
@test nv(ob(colim)) == 2
@test src(ob(colim)) == [1,2]
@test tgt(ob(colim)) == [2,2]
@test is_natural(coproj1(colim)) && is_natural(coproj2(colim))

# Same pushout using generic colimit interface.
#diagram = FreeDiagram([Graph(1), g, ob(terminal(Graph))], [(α,1,2), (β,1,3)])
# ^ This diagram gives colimit that is isomorphic but not equal!
diagram = FreeDiagram([g, ob(terminal(Graph)), Graph(1)], [(α,3,1), (β,3,2)])
ι1, ι2, _ = colim′ = colimit(diagram)
@test ob(colim′) == ob(colim)
@test force(ι1) == force(coproj1(colim))
@test force(ι2) == force(coproj2(colim))

# Attributed C-set morphisms
############################

# Constructors and accessors.
g = path_graph(WeightedGraph{Float64}, 2, E=(weight=2.0,))
h = path_graph(WeightedGraph{Float64}, 4, E=(weight=[1.,2.,3.],))
α = ACSetTransformation((V=[2,3], E=[2]), g, h)
@test length(components(α)) == 2

# Naturality.
@test is_natural(α)
β = ACSetTransformation((V=[2,1], E=[2]), g, h)
@test !is_natural(β) # Preserves weight but not graph homomorphism
β = ACSetTransformation((V=[1,2], E=[1]), g, h)
@test !is_natural(β) # Graph homomorphism but does not preserve weight

# Colimits
#---------

@present TheoryVELabeledGraph <: TheoryGraph begin
  Label::Data
  vlabel::Attr(V,Label)
  elabel::Attr(E,Label)
end
const VELabeledGraph = ACSetType(TheoryVELabeledGraph, index=[:src,:tgt])

# Initial labeled graph.
@test ob(initial(VELabeledGraph{Symbol})) == VELabeledGraph{Symbol}()

# Coproduct of labeled graphs.
g = VELabeledGraph{Symbol}()
add_vertices!(g, 2, vlabel=[:u,:v])
add_edge!(g, 1, 2, elabel=:e)
h = VELabeledGraph{Symbol}()
add_vertex!(h, vlabel=:u)
add_edge!(h, 1, 1, elabel=:f)
coprod = ob(coproduct(g, h))
@test subpart(coprod, :vlabel) == [:u, :v, :u]
@test subpart(coprod, :elabel) == [:e, :f]

# Pushout of labeled graph.
g0 = VELabeledGraph{Symbol}()
add_vertex!(g0, vlabel=:u)
α = ACSetTransformation((V=[1], E=Int[]), g0, g)
β = ACSetTransformation((V=[1], E=Int[]), g0, h)
@test is_natural(α) && is_natural(β)
colim = pushout(α, β)
@test src(ob(colim)) == [1,1]
@test tgt(ob(colim)) == [2,1]
@test subpart(ob(colim), :vlabel) == [:u, :v]
@test subpart(ob(colim), :elabel) == [:e, :f]

α′ = ACSetTransformation(V=[2], E=Int[], g0, g)
@test !is_natural(α′) # Vertex labels don't match.
@test_throws ErrorException pushout(α′, β)

# Finding C-set morphisms
#########################

# Graphs
#-------

g, h = path_graph(Graph, 3), path_graph(Graph, 4)
@test homomorphisms(g, h) == [CSetTransformation((V=[1,2,3], E=[1,2]), g, h),
                              CSetTransformation((V=[2,3,4], E=[2,3]), g, h)]
@test !is_isomorphic(g, h)

I = ob(terminal(Graph))
@test homomorphism(g, I) == CSetTransformation((V=[1,1,1], E=[1,1]), g, I)
@test !is_homomorphic(g, I, monic=true)
@test !is_homomorphic(I, h)

# Graph homomorphism starting from partial assignment, e.g. vertex assignment.
α = CSetTransformation((V=[2,3,4], E=[2,3]), g, h)
@test homomorphisms(g, h, initial=(V=[2,3,4],)) == [α]
@test homomorphisms(g, h, initial=(V=Dict(1 => 2, 3 => 4),)) == [α]
@test homomorphisms(g, h, initial=(E=Dict(1 => 2),)) == [α]
# Inconsistent initial assignment.
@test !is_homomorphic(g, h, initial=(V=Dict(1 => 1), E=Dict(1 => 3)))
# Consistent initial assignment but no extension to complete assignment.
@test !is_homomorphic(g, h, initial=(V=Dict(1 => 2, 3 => 3),))

# Monic and iso on a componentwise basis.
g1, g2 = path_graph(Graph, 3), path_graph(Graph, 2)
add_edges!(g1, [1,2,3,2], [1,2,3,3])  # loops on each node and one double arrow
add_edge!(g2, 1, 2)  # double arrow
@test length(homomorphisms(g2, g1)) == 8 # each vertex + 1->2, and four for 2->3
@test length(homomorphisms(g2, g1, monic=[:V])) == 5 # remove vertex solutions
@test length(homomorphisms(g2, g1, monic=[:E])) == 2 # two for 2->3
@test length(homomorphisms(g2, g1, iso=[:E])) == 0

# Symmetic graphs
#-----------------

g, h = path_graph(SymmetricGraph, 4), path_graph(SymmetricGraph, 4)
αs = homomorphisms(g, h)
@test all(is_natural(α) for α in αs)
@test length(αs) == 16
αs = isomorphisms(g, h)
@test length(αs) == 2
@test map(α -> collect(α[:V]), αs) == [[1,2,3,4], [4,3,2,1]]
g = path_graph(SymmetricGraph, 3)
@test length(homomorphisms(g, h, monic=true)) == 4

# Graph colorability via symmetric graph homomorphism.
# The 5-cycle has chromatic number 3 but the 6-cycle has chromatic number 2.
K₂, K₃ = complete_graph(SymmetricGraph, 2), complete_graph(SymmetricGraph, 3)
C₅, C₆ = cycle_graph(SymmetricGraph, 5), cycle_graph(SymmetricGraph, 6)
@test !is_homomorphic(C₅, K₂)
@test is_homomorphic(C₅, K₃)
@test is_natural(homomorphism(C₅, K₃))
@test is_homomorphic(C₆, K₂)
@test is_natural(homomorphism(C₆, K₂))

# Labeled graphs
#---------------

@present TheoryLabeledGraph <: TheoryGraph begin
  Label::Data
  label::Attr(V,Label)
end
const LabeledGraph = ACSetType(TheoryLabeledGraph, index=[:src,:tgt])

g = cycle_graph(LabeledGraph{Symbol}, 4, V=(label=[:a,:b,:c,:d],))
h = cycle_graph(LabeledGraph{Symbol}, 4, V=(label=[:c,:d,:a,:b],))
@test homomorphism(g, h) == ACSetTransformation((V=[3,4,1,2], E=[3,4,1,2]), g, h)
h = cycle_graph(LabeledGraph{Symbol}, 4, V=(label=[:a,:b,:d,:c],))
@test !is_homomorphic(g, h)

# Sub-C-sets
############

# Construct sub-C-sets.
X = path_graph(Graph, 4)
A = Subobject(X, V=[2,3,4], E=[2,3])
@test Subobject(X, V=[false,true,true,true], E=[false,true,true]) == A
α = hom(A)
@test is_natural(α)
@test dom(α) == path_graph(Graph, 3)
@test codom(α) == X

# Lattice of sub-C-sets.
X = Graph(6)
add_edges!(X, [1,2,3,4,4], [3,3,4,5,6])
A, B = Subobject(X, V=1:4, E=1:3), Subobject(X, V=3:6, E=3:5)
@test A ∧ B == Subobject(X, V=3:4, E=3:3)
@test A ∨ B == Subobject(X, V=1:6, E=1:5)
@test ⊤(X) |> force == Subobject(X, V=1:6, E=1:5)
@test ⊥(X) |> force == Subobject(X, V=1:0, E=1:0)

# Bi-Heyting algebra of sub-C-sets.
#
# Implication in Graph (Reyes et al 2004, Sec 9.1, p. 139).
I = Graph(1)
Y = path_graph(Graph, 3) ⊕ path_graph(Graph, 2) ⊕ path_graph(Graph, 2)
add_vertex!(Y)
add_edge!(Y, 2, 8)
Z = cycle_graph(Graph, 1) ⊕ cycle_graph(Graph, 1)
ιY, ιZ = colim = pushout(CSetTransformation(I, Y, V=[3]),
                         CSetTransformation(I, Z, V=[1]))
B_implies_C, B = Subobject(ιY), Subobject(ιZ)
C = Subobject(ob(colim), V=2:5, E=2:3)
@test (B ⟹ C) == B_implies_C

# Subtraction in Graph (Reyes et al 2004, Sec 9.1, p. 144).
X = ob(colim)
C = Subobject(X, V=2:5, E=[2,3,ne(X)-1])
@test (B \ C) == Subobject(X, V=[nv(X)], E=[ne(X)])

# Negation in Graph (Reyes et al 2004, Sec 9.1, p. 139-140).
X = cycle_graph(Graph, 1) ⊕ path_graph(Graph, 2) ⊕ cycle_graph(Graph, 4)
add_vertex!(X)
add_edge!(X, 4, 8)
A = Subobject(X, V=[2,3,4,5,8], E=[3,7])
neg_A = Subobject(X, V=[1,6,7], E=[1,5])
@test is_natural(hom(A)) && is_natural(hom(neg_A))
@test ¬A == neg_A
@test ¬neg_A == Subobject(X, V=[2,3,4,5,8], E=[2,3,7])

# Non in Graph (Reyes et al 2004, Sec 9.1, p. 144).
X = path_graph(Graph, 5) ⊕ path_graph(Graph, 2) ⊕ cycle_graph(Graph, 1)
A = Subobject(X, V=[1,4,5], E=[4])
non_A = Subobject(X, V=setdiff(vertices(X), 5), E=setdiff(edges(X), 4))
@test ~A == non_A
@test ~non_A == Subobject(X, V=[4,5], E=[4])

# Negation and non in DDS.
S₁ = DDS(); add_parts!(S₁, :X, 5, Φ=[3,3,4,5,5])
S₂ = DDS(); add_parts!(S₂, :X, 3, Φ=[2,3,3])
ι₁, ι₂ = colim = coproduct(S₁, S₂)
S = ob(colim)
A = Subobject(S, X=[3,4,5])
@test ¬A == Subobject(ι₂)
@test ¬Subobject(ι₂) == Subobject(ι₁)
@test ~A == ⊤(S) |> force

# Serialization
###############

function roundtrip_json_acset(x::T) where T <: AbstractACSet
  mktempdir() do dir
    path = joinpath(dir, "acset.json")
    write_json_acset(x, path)
    read_json_acset(T, path)
  end
end

g = star_graph(Graph, 5)
@test roundtrip_json_acset(g) == g

g = path_graph(WeightedGraph{Float64}, 3, E=(weight=[0.5, 1.5],))
@test roundtrip_json_acset(g) == g

g = VELabeledGraph{Symbol}()
add_vertices!(g, 2, vlabel=[:u,:v])
add_edge!(g, 1, 2, elabel=:e)
@test roundtrip_json_acset(g) == g

@present TheoryLabeledDDS <: TheoryDDS begin
  Label::Data
  label::Attr(X, Label)
end
const LabeledDDS = ACSetType(TheoryLabeledDDS, index=[:Φ, :label])

ldds = LabeledDDS{Int}()
add_parts!(ldds, :X, 4, Φ=[2,3,4,1], label=[100, 101, 102, 103])
@test roundtrip_json_acset(ldds) == ldds

end


text:
# see
# https://stackoverflow.com/questions/43066957/adding-global-title-to-plots-jl-subplots
"""
`supertitle(titletext; kw...)`
Create plot axes for super title.

# Arguments
	`titletext`
	- title string, e.g. "My title"
	- Input argument of `text`, i.e. `annotate!(..., text(titletext...))`.
	- For example, one can define `titletext = L"F=ma", 48)`, which results in `annotate!([(..., text(L"F=ma", 48))])`

# How to use:
	p1 = plot(...);
	(title, layout) = supertitle("distribution of step sizes"; size = (200,100), height=0.1);
	plot(title,p1,layout=layout);

# Options:
    size: specifying the size of the superlegend axis
    - if larger size is specified,
        - the font size will be smaller
        - the upper white edge will be larger

    height: determine the ratio of height of the upper axes to the whole. Default is 0.1 (axes for super title has 0.1, and the main plot has 0.9, for `plot(p1,title,...))`.

    ylim: you can set the distance between legend strings to the subplot below and to the edge of the very top. Default is [0,0]

"""
function supertitle(titletext;size=(200,100),height=0.1, ylim=[0,0])
	y = zeros(3);
	# line, scatter, heatmap, etc all just pass through to the plot function. You can imagine it looks like:
	# vline(args...; kw...) = plot(args...; kw..., seriestype = :vline)
	if isa(titletext, String) || isa(titletext, LaTeXString)
		titletext = (titletext,);
	end
	title = Plots.scatter(y, marker=0,markeralpha=0,
		ylim = ylim,
		axis=false, leg=false,grid = false,ticks = nothing,size=size); # annotations=(2, y[2], Plots.text(titletext)),
	annotate!(title, [(2.0, y[2], Plots.text(titletext...))])
	return title, Plots.grid(2,1,heights = [height,1-height])
end

ext: jl


max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
using Test
using ModelConstructors
using Distributions, DSGE, InteractiveUtils
using Nullables, Printf, Random

my_tests = [
            "parameters",
            "distributions_ext",
            "settings",
            "statistics",
            "paths",
            ]

for test in my_tests
    test_file = string("$test.jl")
    @printf " * %s\n" test_file
    include(test_file)
end


text:
using LightGraphs
using GraphIO

function createcnf(g, ncolors)
    colors = 1:ncolors

    getcolor(v, c) = (v - 1) * ncolors + c

    clauses = String[]
    # Vertex clauses.
    for v in vertices(g)
        vars = [-getcolor(v, c) for c in colors]
        for c in colors
            vars[c] *= -1
            push!(clauses, join(map(string, vars), ' ') * " 0")
            vars[c] *= -1
        end
    end

    for e in edges(g)
        v1 = src(e)
        v2 = dst(e)
        for c in colors
            push!(clauses, "$(-getcolor(v1, c)) $(-getcolor(v2, c)) 0")
        end
    end

    nvars = nv(g) * ncolors
    header = "p cnf $nvars $(length(clauses))\n"

    header * join(clauses, '\n') * "\n"
end

function fromfile(ncolors, infile, outfile, fmt=NETFormat())
    g = loadgraph(infile, "g", fmt)
    formula = createcnf(g, ncolors)
    write(outfile, formula)
end

max_stars_repo_name: mtsch/LVR-SAT-Solver

max_issues_repo_name: mtsch/LVR-SAT-Solver

max_forks_repo_name: mtsch/LVR-SAT-Solver


text:
# Rose Teague 06/07/2018
# AdiabaticPropagationTest.jl
#
# Runs the Adiabatic propagation of an eigenstate
push!(LOAD_PATH,"../src/")

using TheDancer
using Base.Test

N, Edisorder, Jdisorder, modelJ, B, dipolestrength, r = init!(3, 0.0, 0.0, 0.001, 298, 20)

"""
Test full adiabatic propagation for a simple 3-site system. With small transfer
integral (J=0.001eV) the electron density should end up all at the central
site.
"""
function main()
    S,E,H,psi,dipoles=prepare_model()
    density = abs.(psi.^2)

    for i in 1:200
        S,H,psi,density,dipoles = AdiabaticPropagation(S,dipoles,E,false)
    end

    ε = 1E-6
    @test  1 ≈ density[2] atol = ε
end

main()


text:
module Tody

export @pseudo_str

import Base: pop!, push!, parse, run

const Nil = Nothing
const nil = nothing

using Base.Iterators
import Base: ==, show, parse, push!, peek, isempty, length, last, insert!

## Expressions ###############################################################

abstract type Expression end

const Ex = Expression

parts(ex::Ex) = ex.parts
parts!(ex::Ex, p) = (empty!(ex.parts); append!(ex.parts, p); ex.parts)

isempty(ex::Ex) = isempty(parts(ex))
length(ex::Ex) = length(parts(ex))
push!(ex::Ex, v) = push!(parts(ex), v)
insert!(ex::Ex, i, v) = insert!(parts(ex), i, v)
head(ex::Ex) = parts(ex)[1]
args(ex::Ex) = @view parts(ex)[2:end]
argn(ex::Ex) = length(args(ex))

==(x::Ex, y::Ex) = head(x) == head(y) && args(x) == args(y)

show(io::IO, x::Ex) = pprint(io, x)

## Function applications #####################################################

struct Call <: Expression
    parts    :: AbstractVector
end
Call() = Call([])
Call(head, args) = Call([head, args...])

## Incremental evaluations ###################################################

## Environments

struct Env
    vars    :: Dict{Symbol, Any}
    par     :: Union{Env, Nil}
end
Env(par=nil) = Env(Dict(), par)

locate(e::Env, name)      = haskey(e.vars, name) ? e : locate(e.par, name)
locate( ::Nil, name)      = error("unknown variable $name")

trylocate(e::Env, name)   = haskey(e.vars, name) ? e : trylocate(e.par, name)
trylocate( ::Nil, name)   = nil

defined(e::Env, name)     = haskey(e.vars, name) ? true : defined(e.par, name)
defined( ::Nil, name)     = false

lookup(e::Env, name)      = locate(e, name).vars[name]

define(e::Env, name, val) = e.vars[name] = val
delete(e::Env, name)      = delete!(e.vars, name)

# Try to rebind existing variable. Define variable if not found:
function assign(e::Env, name, val)
    env = trylocate(e, name)
    env ≡ nil && (env = e)
    define(env, name, val)
end

parent(e::Env) = e.par
parent!(e::Env, par) = e.par = par

## An expression with some of its components evaluated

mutable struct Eval <: Expression
    src     :: Call
    tgt     :: Call
    env     :: Env
    par     :: Union{Eval, Nil}
    idx     :: Int
    frame   :: Bool
end

Eval(ex::Call; env=Env(), par=nil, frame=false) =
    Eval(ex, Call(), env, par, 0, frame)

lookup(e::Eval, name)      = lookup(e.env, name)
define(e::Eval, name, val) = define(e.env, name, val)
assign(e::Eval, name, val) = assign(e.env, name, val)

defined(e::Eval, name)     = defined(e.env, name)

head(e::Eval) = isempty(e.tgt) ? head(e.src) : head(e.tgt)

args(e::Eval) = flatten((args(e.tgt), drop(args(e.src), argn(e.tgt))))
argn(e::Eval) = argn(e.src)

parent(e::Eval)  = e.par
parent!(e::Eval, par) = e.par = par
isframe(e::Eval) = e.frame

done(e::Eval) =      e.idx == length(e.src)
peek(e::Eval) =      (@assert !done(e); parts(e.src)[e.idx + 1])
last(e::Eval) =      last(parts(e.tgt))
skip!(e::Eval) =     (e.idx += 1; e)
push!(e::Eval, v) =  (@assert !done(e); e.idx += 1; push!(e.tgt, v))
rewind!(e::Eval) =   (e.idx = 1; parts!(e.tgt, [head(e.tgt)]); e)

function rewind!(e::Eval, n::Int)
    @assert 1 ≤ n < e.idx
    # Might still trigger indexing error, if some args have been skipped.
    e.idx -= n
    parts!(e.tgt, parts(e.tgt)[1:end-n])
    return e
end

return!(::Nil) = nil

function return!(ev::Eval, v)
    (par = parent(ev)) ≡ nil && return v
    push!(par, v)
    return par
end

## One step in a stepwise evaluation

step!(ev::Eval) = done(ev) ? apply!(ev) : eval!(ev)
eval!(ev::Eval) = eval!(ev, head(ev), peek(ev))
apply!(ev::Eval) = apply!(ev, head(ev), args(ev))

# By default, the head and all arguments are evaluated eagerly. Note that when
# hd ≡ ex, this dispatches on the pre-evaluated head (e.g., a symbol), and for
# the arguments, it dispatches on the evaluated head (e.g., a function).
eval!(ev::Eval, hd, ex)     = eval!(ev, ex)
eval!(ev::Eval, ex)         = noeval!(ev, ex)
eval!(ev::Eval, ex::Symbol) = noeval!(ev, lookup(ev, ex))
eval!(ev::Eval, ex::Call)   = Eval(ex, par=ev, env=ev.env)

apply!(ev::Eval, func::Function, args) = return!(ev, func(args...))

noeval!(ev::Eval, ex)       = (push!(ev, ex); ev)
noapply!(ev::Eval)          = return!(ev, ev)

run(ex::Ex) = Eval(ex, env=BUILTINS())

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

include("util/expr.jl")
include("util/orders.jl")
include("parser/Parsers.jl")
include("core/prelude.jl")

const GRAMMAR  = Parsers.program
const OPPOLICY = Parsers.Ops()
BUILTINS() = prelude(Env())

prelude(OPPOLICY)

function parse(str::AbstractString)

    compacted = replace(str, r"\s*"s => "")
    compacted = compacted == "{}" ? "" : compacted
    isempty(compacted) && return Call([:do])
    res = parse(GRAMMAR, OPPOLICY, str)
    res = res isa Ex ? res : Call([:do, res])

    return res
end

macro pseudo_str(s)
    :($(parse(s)))
end

end


text:
module TestPipelineStatic

using Test
using MLJBase

t  = MLJBase.WrappedFunction(f=log)
f, = fit(t, 0)
@test transform(t, f, 5) ≈ log(5)

end
true


text:
a = open("c:/scratch/python-test.parquet", "w")
seek(a, 0)
write(a,"par1")

close(a)

x = Array{Int32, 1}(undef, 100)


text:
using Plots
using LinearAlgebra
#using Polynomials
using SparseArrays
using Arpack
using ForwardDiff
using BenchmarkTools
using Formatting
using FileIO
mutable struct L
    L1::Float64
    L2::Float64
    L3::Float64
end
mutable struct XYs
    Xs::Array
    Ys::Array
end

println("opening file... \n")
a = open("../programing/efelkismos-hard.unv","r")
line = readlines(a)
line = line[21:end]
nodes = Array{Float64}[]
#read nodes
println("reading nodes...\n")
global i=1
while true
    global i
    l = line[i]
    l = split(l)
    arr = [parse(Float64,k) for k in l]
    if length(arr)==1 & (arr[1]== -1.)
        break
    end
    i = i+2
    push!(nodes,arr)
end
line = line[i+2:end]
#read elements
println("reading elements....\n")
elem = Array{Int64}[]
edge = Array{Int64}[]

i = 1
while true

    global i
    l = line[i]
    l = split(l)
    arr = [parse(Float64,k) for k in l]
    if length(arr)==1 & (arr[1]==-1)
        break
    elseif arr[2]==22
        i = i+3
        l = line[i-1]
        l = split(l)
        arr = [parse(Float64,k) for k in l]
        push!(edge,arr)
    elseif arr[2]==42
        l = line[i+1]
        l = split(l)
        arr = [parse(Float64,k) for k in l]
        push!(elem,arr)
        i = i+2
    end

end
#read boundaries
println("reading boundaries...\n")

function splitToInt(l)
    l = split(l)
    l = [parse(Int64,k) for k in l]
    return l
end
i = i+3 #record 1 of 2467 format
n = splitToInt(line[i])[end] #shows how many elems are to read
line = line[i+2:end]
bc = Array{Int64}[]
global n = div(n,2)+mod(n,2)

while true
global i = 1
global line,n

println("splitting boundaries \n")
while true
    global i,line,n
    if (i>n)
        break
    end
    if tryparse(Float64, split(line[i])[1]) == nothing
        break
    end
    arr = splitToInt(line[i])
    push!(bc,[arr[2]])
    if length(arr)>5
        push!(bc,[arr[6]])
    end
    i = i+1
end
    push!(bc,[-1])
    arr = splitToInt(line[i])
    if arr[1]==-1
        break
    end
    n = arr[end]
    n = div(n,2)+mod(n,2)
    line = line[i+2:end]
end

#=
using RecursiveArrayTools
elem = VectorOfArray(elem); elem = convert(Array,elem)
edge = VectorOfArray(edge); edge = convert(Array,edge)
=#
println("making the arrays \n")
using Flux:batch
nodes = Transpose(batch(nodes))
edge = Transpose(batch(edge))
elem = Transpose(batch(elem))
elem = elem[:,[1 3 5 2 4 6]][:,1,:]
bc = Transpose(batch(bc))

xs = 10*[-1 0 0 -0.5 0 -0.5]
ys = 10*[0 0 1 0 0.5 0.5]
l = [1 0 0]
#defining functions
function tria6(l,xs,ys)
    (l1,l2,l3) = l
    N1 = (2*l1-1)*l1; N2 = (2*l2-1)*l2; N3 = (2*l3-1)*l3;
    N4 = 4*l1*l2; N5 = 4*l2*l3; N6 = 4*l3*l1;
    x = dot(xs,[N1 N2 N3 N4 N5 N6])
    y = dot(ys,[N1 N2 N3 N4 N5 N6])
    return x,y
end
function N(l)
    (l1,l2,l3) = l;
    N1 = (2*l1-1)*l1; N2 = (2*l2-1)*l2; N3 = (2*l3-1)*l3;
    N4 = 4*l1*l2; N5 = 4*l2*l3; N6 = 4*l3*l1;

    N = [N1 0 N2 0 N3 0 N4 0 N5 0 N6 0;
         0 N1 0 N2 0 N3 0 N4 0 N5 0 N6];
    return N
end
function J(l,xs,ys)
    (l1,l2,l3) = l
    Nl1 = [4*l1-1 0 0 4*l2 0 4*l3]
    Nl2 = [0 4*l2-1 0 4*l1 4*l3 0]
    Nl3 = [0 0 4*l3-1 0 4*l2 4*l1]
    Nl = [Nl1; Nl2; Nl3]
    J = Nl*[xs' ys']
end
function Nl(l)
    (l1,l2,l3) = l
    Nl1 = [4*l1-1 0 0 4*l2 0 4*l3]
    Nl2 = [0 4*l2-1 0 4*l1 4*l3 0]
    Nl3 = [0 0 4*l3-1 0 4*l2 4*l1]
    Nl = [Nl1; Nl2; Nl3]
end
function Nxy(l,xs,ys)
    (l1,l2,l3) = l
    j = J(l,xs,ys)
    nl = Nl(l)
    nxy = j\nl
    return nxy
end

function me(xs,ys)
    #neeeds multiplication by rho*h
    me  = zeros(12,12)
    a1 = 0.0597158717; b1 = 0.4701420641; a2 = 0.7974269853; b2 = 0.1012865073
    w1 = 0.2250000000; w2 =0.1323941527; w3=0.1323941527; w4=0.1323941527
    w5 =0.1259391805; w6 = 0.1259391805; w7 = 0.1259391805;
    g1 = [1 1 1]./3; g2 = [a1 b1 b1]; g3 = [b1 a1 b1]; g4 = [b1 b1 a1];
    g5 = [a2 b2 b2]; g6 = [b2 a2 b2]; g7 = [b2 b2 a2]
    w = [w1 w2 w3 w4 w5 w6 w7]
    g = [[g1] [g2] [g3] [g4] [g5] [g6] [g7]]

    for i=1:7
        n = N(g[i])
        me = me + w[i]*n'*n
    end
    return me
end
function Ematrix(e,v,cond)
    if cond=="stress"
        E = (e/(1-v^2))*[1 v 0;
                        v 1 0;
                        0 0 (1-v)/2]
    elseif cond=="strain"
        E = (e*(1-v)/((1+v)*(1-2*v)))*[1 v/(v-1) 0;
                                       v/(v-1) 1 0;
                                        0 0 (1-2*v)/(2-2*v)]
    else
        println("motherfucker no planestress or planestrain condition")
        exit()
    end
    return E
end
function ke(E,xs,ys)
    #neeeds multiplication by h
    me  = zeros(12,12)
    a1 = 0.0597158717; b1 = 0.4701420641; a2 = 0.7974269853; b2 = 0.1012865073
    w1 = 0.2250000000; w2 =0.1323941527; w3=0.1323941527; w4=0.1323941527
    w5 =0.1259391805; w6 = 0.1259391805; w7 = 0.1259391805;
    g1 = [1 1 1]./3; g2 = [a1 b1 b1]; g3 = [b1 a1 b1]; g4 = [b1 b1 a1];
    g5 = [a2 b2 b2]; g6 = [b2 a2 b2]; g7 = [b2 b2 a2]
    w = [w1 w2 w3 w4 w5 w6 w7]
    g = [[g1] [g2] [g3] [g4] [g5] [g6] [g7]]

    ke = zeros(12,12)
    for i=1:7
        nxy = Nxy(g[i],xs,ys)
        nx = nxy[1,:]; ny = nxy[2,:];
        nx = [1 0 0 0 0 0;0 0 0 0 0 0;0 1 0 0 0 0;0 0 0 0 0 0;0 0 1 0 0 0;0 0 0 0 0 0;0 0 0 1 0 0;0 0 0 0 0 0;0 0 0 0 1 0;0 0 0 0 0 0;0 0 0 0 0 1;0 0 0 0 0 0]*nx
        ny = [0 0 0 0 0 0;1 0 0 0 0 0;0 0 0 0 0 0;0 1 0 0 0 0;0 0 0 0 0 0;0 0 1 0 0 0;0 0 0 0 0 0;0 0 0 1 0 0;0 0 0 0 0 0;0 0 0 0 1 0;0 0 0 0 0 0;0 0 0 0 0 1]*ny;
        b = [1 0;0 0;0 1]*[nx';ny'] + [0 0;0 1;1 0]*[nx';ny']

        ke = ke + w[i]*b'*E*b
    end
    return ke
end
function BC(f,xedge,yedge,map)
    #f is vertical array
    #needs to multiplied by h
    #(fx,fy)
    if map[1]==0
        l1 = zeros(7,1)
        l2 = [-0.949108  -0.741531  -0.405845  0.405845  0.741531  0.949108  0.0]
        w = [0.129485  0.279705  0.38183  0.38183  0.279705  0.129485  0.417959]
        #transform to 0~1
        l2 = 0.5 .*l2 .+ 0.5
        l3 = 1 .- l2
        w = 0.5 .* w
    elseif map[2]==0
        l2 = zeros(7,1)
        l3 = [-0.949108  -0.741531  -0.405845  0.405845  0.741531  0.949108  0.0]
        w = [0.129485  0.279705  0.38183  0.38183  0.279705  0.129485  0.417959]
        #transform to 0~1
        l3 = 0.5 .*l3 .+ 0.5
        l1 = 1 .- l3
        w = 0.5 .* w
    elseif map[3]==0
        l3 = zeros(7,1)
        l1 = [-0.949108  -0.741531  -0.405845  0.405845  0.741531  0.949108  0.0]
        w = [0.129485  0.279705  0.38183  0.38183  0.279705  0.129485  0.417959]
        #transform to 0~1
        l1 = 0.5 .*l1 .+ 0.5
        l2 = 1 .- l1
        w = 0.5 .* w
    else
        println("this is not a valid edge map")
        exit()
    end
    bc = zeros(12,1)
    for i=1:7
        l = (l1[i],l2[i],l3[i])
        (x,y) = tria6(l,xedge,yedge)
        bc = bc + N(l)'*f([x,y])
    end
    return bc
end
function ε(l,xs,ys)
    nxy = Nxy(l,xs,ys)
    nx = nxy[1,:]; ny = nxy[2,:];
    nx = [1 0 0 0 0 0;0 0 0 0 0 0;0 1 0 0 0 0;0 0 0 0 0 0;0 0 1 0 0 0;0 0 0 0 0 0;0 0 0 1 0 0;0 0 0 0 0 0;0 0 0 0 1 0;0 0 0 0 0 0;0 0 0 0 0 1;0 0 0 0 0 0]*nx
    ny = [0 0 0 0 0 0;1 0 0 0 0 0;0 0 0 0 0 0;0 1 0 0 0 0;0 0 0 0 0 0;0 0 1 0 0 0;0 0 0 0 0 0;0 0 0 1 0 0;0 0 0 0 0 0;0 0 0 0 1 0;0 0 0 0 0 0;0 0 0 0 0 1]*ny;
    nxy = [nx';ny']
    xy = [collect(zip(xs,ys))[i][k] for i=1:6 for k=1:2]
    b = [1 0;0 0;0 1]*[nx';ny'] + [0 0;0 1;1 0]*[nx';ny']
    epsilon = b*xy
    return epsilon
end
n = length(elem[:,1])
dofs = 2*length(nodes[:,1])
global E = Ematrix(220e9,0.33,"strain")
global K = sparse([],[],Float64[],dofs,dofs);
global M = sparse([],[],Float64[],dofs,dofs);

println("calculating stiffness/mass matrix...\n")
@time begin
for i =1:n
    i1 = elem[i,:]'
    xs = nodes[i1,1]
    ys = nodes[i1,2]
    i1 = [2 .*i1 (i1 .* 2 .-1)]
    kre = 0.1 .* ke(E,xs,ys)
    mre = 0.1 .* 7e3 .* me(xs,ys)
    K[i1,i1] = K[i1,i1][1,:,1,:] + kre
    M[i1,i1] = M[i1,i1][1,:,1,:] + mre
end
end
#start boundary conditions
println("calculating right-hand-side...\n")

F = zeros(dofs,1)
#find indexes of edge that are to be bounded
nbc = count(x->x==-1,bc)+1
ix = findall(x->x==-1,bc)
ix = [ix[i][1] for i=1:nbc-1]
prepend!(ix,0)
ix = [collect(ix[i-1]+1:ix[i]-1) for i=2:nbc]
#ix now contains all indexes, seperated so you can now iterate ix[i], for each i you get the list of edges that are on boundarycondition group i
global bc,edge,elem,nodes,i,F
i=1
#iterate through groups
f(x) = [[0 0]', [0 1e6]']
while true
    if i>nbc-1
        break
    end
    global bc,edge,elem,nodes,rows,xedge,yedge,mapx,i,F
    edix = bc[ix[i]] #contains edge id
    ed = edge[edix,:]#contains edge nodes id
    rows = [i for k=1:length(ix[i]) for i=1:length(elem[:,1]) if issubset(ed[k,1],elem[i,:]) & issubset(ed[k,2],elem[i,:]) & issubset(ed[k,3],elem[i,:])] #contains superset id element
    xedge = nodes[elem[rows,:],:][:,:,1]; yedge = nodes[elem[rows,:],:][:,:,2] #contains x,y of edges
    #zero the non edge indexes of element
    superelem = elem[rows,:]
    mapx = [sum(dims=2,superelem[i,:] .== ed[i,:]') for i=1:length(superelem[:,1])]
    for k=1:length(mapx)
        global bc,edge,elem,nodes,rows,xedge,yedge,mapx,i,F
        g(x) = f(x)[i]
        a = BC(g,xedge[k,:],yedge[k,:],mapx[k])
        id = [2 .*elem[rows[k],:] (2 .* elem[rows[k],:].-1)]
        F[id] = a
    end
    i = i+1
end
#gib traction nodes where (u,v) are known
#reminder, ix contains boundary condition nodes INDEXES
println("solving for U....\n")
freedofs = zeros(dofs,1)
fixedofs = edge[bc[ix[1]],:]
fixedofs = elem[[i for k=1:length(ix[1]) for i=1:length(elem[:,1]) if issubset(fixedofs[k,1],elem[i,:]) & issubset(fixedofs[k,2],elem[i,:]) & issubset(fixedofs[k,3],elem[i,:])],:][:] #contains superset id element

freedofs[[2 .*fixedofs' (2 .* fixedofs' .-1)]] .= 1
freedofs = [i for i=1:dofs if freedofs[i]==0]
fixedofs = [2 .* fixedofs' (2 .* fixedofs' .-1)]'
U = zeros(dofs,1)
rhs = F
rhs = rhs[freedofs] - K[freedofs,fixedofs][:,:,1]*U[fixedofs]
Kfree = K[freedofs,freedofs]
Mfree = M[freedofs,freedofs]

Ufree = U[freedofs]
U[freedofs] = Kfree\rhs

n = length(elem[:,1])
for i=1:n
    i1 = elem[i,:]'
    xs = nodes[i1,1]
    ys = nodes[i1,2]
    i1 = [2 .*i1 (i1 .* 2 .-1)]
end

println("calculating element average stress....\n")
#calculate element stress
stress = zeros(n,3)
for i=1:length(elem[:,1])
    l =  [1 0 0;
          0 1 0;
          0 0 1;
          0.5 0.5 0;
          0 0.5 0.5;
          0.5 0 0.5]

    i1 = elem[i,:]'
    xs = nodes[i1,1]
    ys = nodes[i1,2]
    mean = zeros(3,1)
    for k =1:6
        epsilon = ε(l[k,:],xs,ys)
        sigma = E*epsilon
        mean = mean + sigma/6;
     end
     stress[i,:] = mean
 end

#println(stress')
println("DONE! \n")

max_stars_repo_name: billyon/first

max_issues_repo_name: billyon/first

max_forks_repo_name: billyon/first


text:
# A set of simple mesh makers for testing
export simple_line_mesh, simple_quad_mesh, simple_hex_mesh

#=
# Builds a 1D interval mesh.
# nx = number of vertices
# bn = number of boundary regions
# interval = limits of the square domain
=#
function simple_line_mesh(nx, bn, interval)
    # mesh only
    Nv = nx;                    # number of vertex nodes
    xv = zeros(1,Nv);           # coordinates of vertices
    ind = Array(1:Nv);          # indices are in order
    nel = nx-1;                 # number of elements
    el = zeros(Int, 2, nel);    # element vertex maps
    etypes = ones(Int, nel);    # element types (gmsh number)
    numvert = 2*ones(Int, nel); # number of vertices per element
    invind = invert_index(ind); # inverse index ind[i] = j -> invind[j] = i
    f2n = zeros(Int, 1,Nv);     # face2node mapping
    f2e = zeros(Int, 2,Nv);     # face2element mapping
    e2f = zeros(Int, 2,nel);    # element2face mapping
    normals = ones(1,Nv);      # normals of faces
    bdryID = zeros(Int, Nv);    # BID of faces (0=interior)

    scale = interval[2]-interval[1];
    h = scale/(nx-1);               # uniformly divided

    # vertex nodes are initially ordered lexicographically (that's a 17 letter word!)
    for j=1:nx
        xv[1,j] = interval[1] + (j-1)*h;
    end

    # Elements are ordered the same way
    for ei=1:nel
        el[1,ei] = ei;
        el[2,ei] = ei+1;
        e2f[1,ei] = ei;
        e2f[2,ei] = ei+1;
    end

    # face2node, face2element, normals
    for fi=1:Nv
        f2n[1,fi] = fi;
        f2e[1,fi] = fi-1;
        f2e[2,fi] = fi;
        # normals[fi] = 1; # already set to 1
        # bdryID[fi] = 0; # already set to 0
    end
    f2e[2,Nv] = 0;
    bdryID[1] = 1;
    f2e[:,1] = [1,0];
    normals[1] = -1;
    if bn == 2
        bdryID[Nv] = 2; # for two BIDs
    else
        bdryID[Nv] = 1; # for one BID
    end

    mesh = MeshData(Nv, xv, ind, nel, el, etypes, numvert, invind, f2n, f2e, e2f, normals, bdryID); # MeshData struct

    return mesh;
end

#=
# Builds a 2D triangle mesh
# nx = number of vertices in each direction
# bn = number of boundary regions
# interval = limits of the square domain
=#
function simple_tri_mesh(nxy, bn, interval)
    if length(nxy) == 2
        nx = nxy[1];
        ny = nxy[2];
    else
        nx = nxy;
        ny = nx;
    end
    if length(interval) == 2
        interval = [interval[1], interval[2], interval[1], interval[2]];
    end

    # mesh
    Nv = nx*ny;                 # number of vertex nodes
    xv = zeros(2,Nv);           # coordinates of vertices
    ind = Array(1:Nv);          # indices are in order
    nel = (nx-1)*(ny-1)*2;      # number of elements
    el = zeros(Int, 3, nel);    # element vertex maps
    etypes = fill(2, nel);      # element types (gmsh number)
    numvert = fill(3, nel);     # number of vertices per element
    invind = invert_index(ind); # inverse index ind[i] = j -> invind[j] = i
    Nf = nx*(ny-1) + ny*(nx-1) + (nx-1)*(ny-1); # number of faces
    f2n = zeros(Int, 2,Nf);     # face2node mapping
    f2e = zeros(Int, 2,Nf);     # face2element mapping
    e2f = zeros(Int, 3,nel);    # element2face mapping
    normals = ones(2,Nf);       # normals of faces
    bdryID = zeros(Int, Nf);    # BID of faces (0=interior)

    bids = [1]; # everywhere
    allbids = [1,1,1,1];

    scalex = interval[2]-interval[1];
    hx = scalex/(nx-1);         # uniformly divided
    scaley = interval[4]-interval[3];
    hy = scaley/(ny-1);         # uniformly divided
    hd = sqrt(hy*hy + hx*hx);   # diagonal length

    # vertex nodes are ordered lexicographically
    for j=1:ny
        for i=1:nx
            k = i + (j-1)*nx;
            xv[1, k] = interval[1] + (i-1)*hx;
            xv[2, k] = interval[3] + (j-1)*hy;
        end
    end

    # Elements are ordered like
    #  _____ _____
    # |\    |\    |
    # | \ 2 | \ 4 |
    # |  \  |  \  |
    # | 1 \ | 3 \ |
    # |____\|____\|
    for j=1:(ny-1)
        for i=1:(nx-1)
            ei1 = i*2-1 + (j-1)*2*(nx-1); # element index
            ei2 = i*2   + (j-1)*2*(nx-1);

            el[1,ei1] = i + (j-1)*nx;
            el[2,ei1] = i + (j-1)*nx + 1;
            el[3,ei1] = i + (j)*nx;

            el[3,ei2] = i + (j-1)*nx + 1;
            el[1,ei2] = i + (j)*nx + 1;
            el[2,ei2] = i + (j)*nx;

            # face2node, face2element, element2face, normals
            f1 = 3*(ei-1)+1 + (j-1); # left face index
            f2 = f1+1; # bottom
            f3 = f1+2; # center
            f4 = f1+3; # right
            f5 = j==(ny-1) ? Nf-(nx-1)+i : j*3*(nx-1) + j + i*3 - 1; # top

            f2n[1,f1] = el[3,ei1];
            f2n[2,f1] = el[1,ei1];
            f2n[1,f2] = el[1,ei1];
            f2n[2,f2] = el[2,ei1];
            f2n[1,f3] = el[2,ei1];
            f2n[2,f3] = el[3,ei1];
            f2n[1,f4] = el[3,ei2];
            f2n[2,f4] = el[1,ei2];
            f2n[1,f5] = el[1,ei2];
            f2n[2,f5] = el[2,ei2];

            f2e[2,f1] = ei1;
            f2e[2,f2] = ei1;
            f2e[1,f3] = ei1;
            f2e[2,f3] = ei2;
            f2e[1,f4] = ei2;
            f2e[1,f4] = ei2;

            e2f[1,ei1] = f2;
            e2f[2,ei1] = f3;
            e2f[3,ei1] = f1;
            e2f[1,ei2] = f5;
            e2f[2,ei2] = f3;
            e2f[3,ei2] = f4;

            normals[:,f1] = [1,0];
            normals[:,f2] = [0,1];
            normals[:,f3] = [hy/hd, hx/hd];
            normals[:,f4] = [1,0];
            normals[:,f5] = [0,1];
        end
    end

    # boundaries
    for j=1:(ny-1)
        eleft = (j-1)*2*(nx-1) + 1;
        eright = j*2*(nx-1);

        fleft = e2f[3,eleft];
        fright = e2f[3,eright];

        bdryID[fleft] = 1; # always 1
        bdryID[fright] = allbids[2];

        # need to change normals, f2e for left side
        normals[:,fleft] = [-1,0];
        f2e[:,fleft] = [eleft, 0];
    end
    for i=1:(nx-1)
        ebottom = i*2-1;
        etop = nel - 2*(nx-1) + 2*i;

        fbottom = e2f[1,ebottom];
        ftop = e2f[1,etop];

        bdryID[fbottom] = allbids[3];
        bdryID[ftop] = allbids[4];

        # need to change normals, f2e for bottom side
        normals[:,fbottom] = [0,-1];
        f2e[:,fbottom] = [ebottom, 0];
    end

    mesh = MeshData(Nv, xv, ind, nel, el, etypes, numvert, invind, f2n, f2e, e2f, normals, bdryID); # MeshData struct

    return mesh;
end

#=
# Builds a 2D quad mesh
# nx = number of vertices
# bn = number of boundary regions
# interval = limits of the square domain
=#
function simple_quad_mesh(nxy, bn, interval)
    if length(nxy) == 2
        nx = nxy[1];
        ny = nxy[2];
    else
        nx = nxy;
        ny = nx;
    end
    if length(interval) == 2
        interval = [interval[1], interval[2], interval[1], interval[2]];
    end

    # mesh
    Nv = nx*ny;                 # number of vertex nodes
    xv = zeros(2,Nv);           # coordinates of vertices
    ind = Array(1:Nv);          # indices are in order
    nel = (nx-1)*(ny-1);        # number of elements
    el = zeros(Int, 4, nel);    # element vertex maps
    etypes = 3*ones(Int, nel);  # element types (gmsh number)
    numvert = 4*ones(Int, nel); # number of vertices per element
    invind = invert_index(ind); # inverse index ind[i] = j -> invind[j] = i
    Nf = nx*(ny-1) + ny*(nx-1); # number of faces
    f2n = zeros(Int, 2,Nf);     # face2node mapping
    f2e = zeros(Int, 2,Nf);     # face2element mapping
    e2f = zeros(Int, 4,nel);    # element2face mapping
    normals = ones(2,Nf);       # normals of faces
    bdryID = zeros(Int, Nf);    # BID of faces (0=interior)

    bids = [1]; # everywhere
    allbids = [1,1,1,1];

    scalex = interval[2]-interval[1];
    hx = scalex/(nx-1);         # uniformly divided
    scaley = interval[4]-interval[3];
    hy = scaley/(ny-1);         # uniformly divided

    # vertex nodes are ordered lexicographically
    for j=1:ny
        for i=1:nx
            k = i + (j-1)*nx;
            xv[1, k] = interval[1] + (i-1)*hx;
            xv[2, k] = interval[3] + (j-1)*hy;
        end
    end

    # Elements are ordered the same way
    for j=1:(ny-1)
        for i=1:(nx-1)
            ei = i + (j-1)*(nx-1); # element index

            el[1,ei] = i + (j-1)*nx;
            el[2,ei] = i + (j-1)*nx + 1;
            el[3,ei] = i + (j)*nx + 1;
            el[4,ei] = i + (j)*nx;

            # face2node, face2element, element2face, normals
            f1 = 2*(ei-1)+1 + (j-1); # left face index
            f2 = f1+1; # bottom
            f3 = f1+2; # right
            f4 = j==(ny-1) ? Nf-(nx-1)+i : j*2*(nx-1) + j + i*2; # top

            f2n[1,f1] = el[1,ei];
            f2n[2,f1] = el[4,ei];
            f2n[1,f2] = el[1,ei];
            f2n[2,f2] = el[2,ei];
            f2n[1,f3] = el[2,ei];
            f2n[2,f3] = el[3,ei];
            f2n[1,f4] = el[4,ei];
            f2n[2,f4] = el[3,ei];

            f2e[2,f1] = ei;
            f2e[2,f2] = ei;
            f2e[1,f3] = ei;
            f2e[1,f4] = ei;

            e2f[1,ei] = f1;
            e2f[2,ei] = f2;
            e2f[3,ei] = f3;
            e2f[4,ei] = f4;

            normals[:,f1] = [1,0];
            normals[:,f2] = [0,1];
            normals[:,f3] = [1,0];
            normals[:,f4] = [0,1];
        end
    end

    # boundaries
    for j=1:(ny-1)
        eleft = (j-1)*(nx-1) + 1;
        eright = j*(nx-1);

        fleft = e2f[1,eleft];
        fright = e2f[3,eright];

        bdryID[fleft] = 1; # always 1
        bdryID[fright] = allbids[2];

        # need to change normals, f2e for left side
        normals[:,fleft] = [-1,0];
        f2e[:,fleft] = [eleft, 0];
    end
    for i=1:(nx-1)
        ebottom = i;
        etop = nel - (nx-1) + i;

        fbottom = e2f[2,ebottom];
        ftop = e2f[4,etop];

        bdryID[fbottom] = allbids[3];
        bdryID[ftop] = allbids[4];

        # need to change normals, f2e for bottom side
        normals[:,fbottom] = [0,-1];
        f2e[:,fbottom] = [ebottom, 0];
    end

    mesh = MeshData(Nv, xv, ind, nel, el, etypes, numvert, invind, f2n, f2e, e2f, normals, bdryID); # MeshData struct

    return mesh;
end

#=
# Builds a 3D hex mesh
# nx = number of vertices
# bn = number of boundary regions
# interval = limits of the square domain
=#
function simple_hex_mesh(nxyz, bn, interval)
    if length(nxyz) == 3
        nx = nxyz[1];
        ny = nxyz[2];
        nz = nxyz[3]
    else
        nx = nxyz;
        ny = nx;
        nz = nx;
    end
    if length(interval) == 2
        interval = [interval[1], interval[2], interval[1], interval[2], interval[1], interval[2]];
    end

    # mesh
    Nv = nx*ny*nz;              # number of vertex nodes
    xv = zeros(3,Nv);           # coordinates of vertices
    ind = Array(1:Nv);          # indices are in order
    nel = (nx-1)*(ny-1)*(nz-1); # number of elements
    el = zeros(Int, 8, nel);    # element vertex maps
    etypes = 5*ones(Int, nel);  # element types (gmsh number)
    numvert = 8*ones(Int, nel); # number of vertices per element
    invind = invert_index(ind); # inverse index ind[i] = j -> invind[j] = i
    Nf = nz*(nx-1)*(ny-1) + ny*(nx-1)*(nz-1) + nx*(nz-1)*(ny-1); # number of faces
    f2n = zeros(Int, 4,Nf);     # face2node mapping
    f2e = zeros(Int, 2,Nf);     # face2element mapping
    e2f = zeros(Int, 6,nel);    # element2face mapping
    normals = ones(3,Nf);       # normals of faces
    bdryID = zeros(Int, Nf);    # BID of faces (0=interior)

    bids = [1]; # everywhere
    allbids = [1,1,1,1,1,1];

    scalex = interval[2]-interval[1];
    hx = scalex/(nx-1);         # uniformly divided
    scaley = interval[4]-interval[3];
    hy = scaley/(ny-1);         # uniformly divided
    scalez = interval[6]-interval[5];
    hz = scalez/(nz-1);         # uniformly divided

    # vertex nodes are ordered lexicographically
    for k=1:nz
        for j=1:ny
            for i=1:nx
                ni = i + (j-1)*nx + (k-1)*nx*ny;
                xv[1, ni] = interval[1] + (i-1)*hx;
                xv[2, ni] = interval[3] + (j-1)*hy;
                xv[3, ni] = interval[5] + (k-1)*hz;
            end
        end
    end

    # Elements are ordered lexicographically
    # Elemental vertices are ordered according to GMSH order
    nextface = 1;
    for k=1:(nz-1)
        for j=1:(ny-1)
            for i=1:(nx-1)
                ei = i + (j-1)*(nx-1) + (k-1)*(ny-1)*(nx-1); # element index

                el[1,ei] = i + (j-1)*nx + (k-1)*nx*ny;
                el[2,ei] = i + (j-1)*nx + (k-1)*nx*ny + 1;
                el[3,ei] = i + (j)*nx + (k-1)*nx*ny + 1; #GMSH
                el[4,ei] = i + (j)*nx + (k-1)*nx*ny; #GMSH
                el[5,ei] = i + (j-1)*nx + (k)*nx*ny;
                el[6,ei] = i + (j-1)*nx + (k)*nx*ny + 1;
                el[7,ei] = i + (j)*nx + (k)*nx*ny + 1; #GMSH
                el[8,ei] = i + (j)*nx + (k)*nx*ny; #GMSH

                # face2node, face2element, element2face, normals
                if i==1
                    f1 = nextface;
                    nextface = nextface+1;
                else
                    f1 = e2f[2,ei-1];# left face index
                end
                f2 = nextface; # right
                nextface = nextface+1;

                if j==1
                    f3 = nextface;
                    nextface = nextface+1;
                else
                    f3 = e2f[4,ei-(nx-1)]; # bottom
                end
                f4 = nextface; # top
                nextface = nextface+1;

                if k==1
                    f5 = nextface;
                    nextface = nextface+1;
                else
                    f5 = e2f[6,ei-((nx-1)*(ny-1))]; # front
                end
                f6 = nextface; # back
                nextface = nextface+1;

                f2n[:,f1] = el[[1,4,5,8],ei]; # modified to match GMSH
                f2n[:,f2] = el[[2,3,6,7],ei];
                f2n[:,f3] = el[[1,2,5,6],ei];
                f2n[:,f4] = el[[4,3,7,8],ei];
                f2n[:,f5] = el[[1,2,3,4],ei];
                f2n[:,f6] = el[[5,6,7,8],ei];

                f2e[2,f1] = ei;
                f2e[1,f2] = ei;
                f2e[2,f3] = ei;
                f2e[1,f4] = ei;
                f2e[2,f5] = ei;
                f2e[1,f6] = ei;

                e2f[1,ei] = f1;
                e2f[2,ei] = f2;
                e2f[3,ei] = f3;
                e2f[4,ei] = f4;
                e2f[5,ei] = f5;
                e2f[6,ei] = f6;

                normals[:,f1] = [1,0,0];
                normals[:,f2] = [1,0,0];
                normals[:,f3] = [0,1,0];
                normals[:,f4] = [0,1,0];
                normals[:,f5] = [0,0,1];
                normals[:,f6] = [0,0,1];
            end
        end
    end

    # boundaries
    for k=1:(nz-1)
        for j=1:(ny-1)
            eleft = (k-1)*(nx-1)*(ny-1) + (j-1)*(nx-1) + 1;
            eright = (k-1)*(nx-1)*(ny-1) + j*(nx-1);

            fleft = e2f[1,eleft];
            fright = e2f[2,eright];

            bdryID[fleft] = 1; # always 1
            bdryID[fright] = allbids[2];

            # need to change normals, f2e for left side
            normals[:,fleft] = [-1,0,0];
            f2e[:,fleft] = [eleft, 0];
        end
    end

    for k=1:(nz-1)
        for i=1:(nx-1)
            ebottom = (k-1)*(nx-1)*(ny-1) + i;
            etop = (k-1)*(nx-1)*(ny-1) + (nx-1)*(ny-2) + i;

            fbottom = e2f[3,ebottom];
            ftop = e2f[4,etop];

            bdryID[fbottom] = allbids[3];
            bdryID[ftop] = allbids[4];

            # need to change normals, f2e for bottom side
            normals[:,fbottom] = [0,-1, 0];
            f2e[:,fbottom] = [ebottom, 0];
        end
    end

    for j=1:(ny-1)
        for i=1:(nx-1)
            efront = (j-1)*(nx-1) + i;
            eback = (nz-2)*(nx-1)*(ny-1) + efront;

            ffront = e2f[5,efront];
            fback = e2f[6,eback];

            bdryID[ffront] = allbids[5];
            bdryID[fback] = allbids[6];

            # need to change normals, f2e for bottom side
            normals[:,ffront] = [0,0,-1];
            f2e[:,ffront] = [efront, 0];
        end
    end

    mesh = MeshData(Nv, xv, ind, nel, el, etypes, numvert, invind, f2n, f2e, e2f, normals, bdryID); # MeshData struct

    return mesh
end

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
using BAT, ValueShapes, IntervalSets, Distributions, Plots, EmpiricalDistributions
using StatsBase, ArraysOfArrays, LinearAlgebra, LaTeXStrings, QuadGK, PrettyTables, HypothesisTests, Statistics
ENV["JULIA_DEBUG"] = "BAT"

function setup_benchmark()
    if(!(("plots1D" in readdir()) && ("plots2D" in readdir()) && ("results" in readdir())))
        mkpath("benchmark_results")
        cd("benchmark_results")
        mkpath("plots1D")
        mkpath("plots2D")
        mkpath("results")
    end

    include("utils.jl")
    #include("functions_1D.jl")
    include("functions_2D.jl")
    #include("run_benchmark_1D.jl")
    include("run_benchmark_2D.jl")
    include("run_benchmark_ND.jl")
    end

function do_benchmarks(;algorithm=MetropolisHastings(), n_steps=10^5, n_chains=8)
    #run_1D_benchmark(algorithm=algorithm, n_steps=n_steps, n_chains=n_chains)
    run_2D_benchmark(algorithm=algorithm, n_steps=n_steps, n_chains=n_chains)
    run_ND_benchmark(n_dim=2:2:20,algorithm=MetropolisHastings(), n_steps=2*10^5, n_chains=4)
    run_ks_ahmc_vs_mh(n_dim=20:5:35)
end

setup_benchmark()
do_benchmarks()


text:
using DrWatson
@quickactivate "OptimalEstimationProject"
using MATLAB
using IndirectTrajOpt
using StaticArrays
using LinearAlgebra

# Data info
tf      = 33.1

# Init parameters
ps = initCR3BPIndirectParams("Low Thrust 10 CR3BP")
μ = ps.crp.μ
LU = ps.crp.LU
TU = ps.crp.TU

# Trajectory time span
tspan = (0.0, tf)
tspanPlot = (0.0, tspan[2] * 86400 / TU)

# Initial conditions
ics = [-0.019488511458668, -0.016033479812051, 0.0,
        8.918881923678198, -4.081793688818725, 0.0,
        1.0]

# Target final conditions
fcs = [(1.0 - μ) + 6060.483/LU, 19452.284/LU, -34982.968/LU, 0.082677*(TU/LU), 0.006820*(TU/LU), -0.368434*(TU/LU), 0.0]
xf = @SVector [fcs[1], fcs[2], fcs[3], fcs[4], fcs[5], fcs[6]]

# Integrate target trajectory
solTarg = integrate(xf, (0.0, 1.6), "Low Thrust 10 CR3BP", CR3BP(), FullSolutionHistoryNoControl())
rxTarg = [u[1] for u in solTarg.u]
ryTarg = [u[2] for u in solTarg.u]
rzTarg = [u[3] for u in solTarg.u]

#for i in 1:6
#    fcs[i] = solTarg(0.16)[i]
#end

# Grab requirements from DataFrame
#λ0 = [3.0872240069228414, 7.585751775356703, -0.8828499338288641, -0.022292151740247245, 0.01139245771460113, -0.00015851392516777712, 0.13268783730938283]
#ps.ϵ = 0.16662720666508826
#λ0 = [5.963713007087376, 12.99948134789672, 1.6584130970338573, -0.04062357611978672, 0.018558730103857064, -0.0004993971529679529, 0.1220053310292062]
λ0 = [5.963713007227692, 12.999481348176722, 1.6584130970663207, -0.04062357612068792, 0.018558730104232632, -0.0004993971529784666, 0.12200533102974148]
ps.ϵ = 0.0
#λ0 = GetHomotopySolutionVector(fss)[1]

prob = IndirectOptimizationProblem("Low Thrust 10 CR3BP", ics, fcs, tspan)
fss = FSSSolver(λ0, ics, fcs, prob.BVPFunc, prob.BVPWithSTMFunc;
                homotopy = true, homotopyParamVec = [0.0])
#solve!(fss; factor = 3.0, ftol = 1e-10)
#λ0 = GetHomotopySolutionVector(fss)[1]

# Initialize final integration inputs
y0 = @SVector [ics[1], ics[2], ics[3], ics[4], ics[5], ics[6], ics[7],
            λ0[1], λ0[2], λ0[3], λ0[4], λ0[5], λ0[6], λ0[7]]

# Integrate
sol = integrate(y0, tspanPlot, ps, CR3BP(), FullSolutionHistory(), MEMF())

# Compute thrusting and coasting arcs
cSc = ps.sp.c * ps.crp.TU / (ps.crp.LU * 1000.0)
rxt = Vector{Float64}(undef, length(sol)); fill!(rxt, NaN)
ryt = Vector{Float64}(undef, length(sol)); fill!(ryt, NaN)
rzt = Vector{Float64}(undef, length(sol)); fill!(rzt, NaN)
rxc = Vector{Float64}(undef, length(sol)); fill!(rxc, NaN)
ryc = Vector{Float64}(undef, length(sol)); fill!(ryc, NaN)
rzc = Vector{Float64}(undef, length(sol)); fill!(rzc, NaN)
Sold = 0.0
for i in 1:length(sol)
    λv = norm(@view(sol.u[i][11:13]))
    S = IndirectTrajOpt.computeS(sol.u[i], λv, cSc)
    #println(S)
    i == 1 ? Sold = S : ()
    if S > ps.ϵ || Sold > ps.ϵ
        rxc[i] = sol.u[i][1]
        ryc[i] = sol.u[i][2]
        rzc[i] = sol.u[i][3]
    end
    if S <= ps.ϵ || Sold <= ps.ϵ
        rxt[i] = sol.u[i][1]
        ryt[i] = sol.u[i][2]
        rzt[i] = sol.u[i][3]
    end
    global Sold = S
end

mat"""
figure
plot3(-$μ*$LU, 0, 0, 'xk')
hold on
plot3((1-$μ)*$LU, 0, 0, 'xk')
plot3($rxc*$LU,$ryc*$LU,$rzc*$LU, 'b')
plot3($rxt*$LU,$ryt*$LU,$rzt*$LU, 'r')
plot3($rxTarg*$LU, $ryTarg*$LU, $rzTarg*$LU, '--k')
axis equal
grid on
xlabel("x, km")
ylabel("y, km")
zlabel("z, km")
"""

max_stars_repo_name: GrantHecht/OptimalEstimationProject


text:
include(joinpath("inc", "algorithm_help_msg.jl"));

using ArgParse;
using Plots;

s = ArgParseSettings();
@add_arg_table! s begin
  "--E0", "-z"
    help = "magnitude of electric field"
    arg_type = Float64
    default = 0.0
  "--dk", "-D"
    help = "kappa2 - kappa1"
    arg_type = Float64
    default = 1.0
  "--kT", "-k"
    help = "dimensionless temperature"
    arg_type = Float64
    default = 1.0
  "--gz", "-g"
    help = "dimensionless stretch in the z-direction (direction of E-field)"
    arg_type = Float64
    default = 0.0
  "--gx", "-f"
    help = "dimensionless stretch in the x-direction"
    arg_type = Float64
    default = 0.0
  "--num-monomers", "-N"
    help = "number of monomers"
    arg_type = Int
    default = 100
  "--algorithm", "-A"
    help = "algorithm (ISRES | PRAXIS | SBPLX | SLSQP | LBFGS | AUGLAG | Newton | Anneal). Run with --print-algo-help for more detailed information"
    default = "try Newton then PRAXIS,SBPLX,Newton"
  "--print-algo-help"
    help = "Get detailed information regarding algorithms themselves, and how to run solvers in stages"
    action = :store_true
  "--print-every-percent-compl"
    help = "Print diagnostic info every percentage complete of maximum iterations (should be between 0 and 1)"
    arg_type = Float64
    default = 0.1
  "--verbose", "-v"
    help = "print diagnostic information"
    action = :store_true
  "--reltol-int", "-r"
    help = "relative tolerance for numerical integration"
    arg_type = Float64
    default = 1e-6
  "--maxevals-int", "-e"
    help = "maximum evaluations for numerical integration"
    arg_type = Int
    default = 2500
  "--reltol-opt", "-R"
    help = "relative tolerance for optimization"
    arg_type = Float64
    default = 0.0
  "--abstol-opt", "-B"
    help = "absolute tolerance for optimization"
    arg_type = Float64
    default = 1e-8
  "--maxevals-opt", "-E"
    help = "maximum evaluations for numerical integration"
    arg_type = String
    default = "try 1000 then 10000,5000,50"
  "--init-guess", "-X"
    help = "initial guess"
    arg_type = String
  "--anneal-factor", "-m"
    help = "annealing factor, gamma i.e. T_new = T * (1 - gamma)"
    arg_type = Float64
    default = 0.005
  "--anneal-init-T", "-t"
    help = "initial temperature"
    arg_type = Float64
    default = 25.0
  "--anneal-final-T", "-T"
    help = "final temperature"
    arg_type = Float64
    default = 0.0
  "--anneal-relax-steps", "-S"
    help = "number of steps before annealing and step size adaptation"
    arg_type = Int
    default = 100
  "--anneal-init-delta", "-s"
    help = "initial step size for annealing"
    arg_type = Float64
    default = 0.5
  "--anneal-adapt-factor", "-P"
    help = "step size adapation factor (should be greater than 1.0)"
    arg_type = Float64
    default = 1.2
  "--anneal-acc-thrshld", "-a"
    help = "step size adapation acceptance threshold (should be in [0, 1))"
    arg_type = Float64
    default = 0.7
  "--plot-dist", "-p"
    help = "plot monomer distribution"
    action = :store_true
end

pa = parse_args(s);
if pa["print-algo-help"]
  println(_ALGO_HELP_MSG_);
  exit();
end

const dk = pa["dk"];
const E0 = pa["E0"];
const kT = pa["kT"];
const ω0 = dk * E0*E0 / (kT);
const _u0_ = (dk > 0) ? -dk*(E0^2) : 0.0;
const γz = pa["gz"];
const γx = pa["gx"];
const N = pa["num-monomers"];
const reltol_opt = pa["reltol-opt"];
const abstol_opt = pa["abstol-opt"];

const verbose = pa["verbose"];
const print_every_percent_compl = pa["print-every-percent-compl"];
const RELTOL_INT = pa["reltol-int"];
const MAXEVALS_INT = pa["maxevals-int"];
const ANNEAL_FACTOR = pa["anneal-factor"];
const ANNEAL_INIT_T = pa["anneal-init-T"];
const ANNEAL_FINAL_T = pa["anneal-final-T"];
const ANNEAL_RELAX_STEPS = pa["anneal-relax-steps"];
const ANNEAL_INIT_DELTA = pa["anneal-init-delta"];
const ANNEAL_ADAPT_FACTOR = pa["anneal-adapt-factor"];
const ANNEAL_ACC_THRSHLD = pa["anneal-acc-thrshld"];

include(joinpath("inc", "solvers.jl"));
include(joinpath("inc", "free_energy.jl"));
include(joinpath("inc", "dipole.jl"));
include(joinpath("inc", "output.jl"));

_RHS_ = Float64[N; γz*N; γx*N];
ω(ϕ::Real, θ::Real) = (ω0 * (cos(θ))^2);

main_algos = parse_algostr(pa["algorithm"]);
main_max_iters = parse_max_iters(pa["maxevals-opt"]);
 DEFAULT_INIT_GUESS = (
  if pa["init-guess"] != nothing
    eval(parse(pa["init-guess"]));
  #=elseif ω == 0
    initial_guess(γx, γz);
  elseif abs(ω) < 1
    ((1-abs(ω))*initial_guess(γx, γz) +
     (abs(ω))*initial_guess(ω, (γx, γz)));
  else
    initial_guess(ω, (γx, γz));=#
  else
    initial_guess(γx, γz);
  end
);
(xs, rs, converged, iterations) = run_solver(main_algos, main_max_iters;
                                             x0s = DEFAULT_INIT_GUESS,
                                             reltol = reltol_opt,
                                             abstol = abstol_opt);
if verbose
  println((converged) ? "Solution converged." : "Solution did NOT converge.");
  println("Iterations: ", iterations);
end
residuals_output(rs);
standard_solution_output(xs);
if pa["plot-dist"]
  ρ_polar(θ) = (pquadrature(ϕ -> ρ(ϕ, θ, xs[1], xs[2], xs[3]), 0.0, 2*π))[1];
  p = plot(range(0, 2*π; length=1000), ρ_polar; proj=:polar);
  display(p);
  println("Press <Enter> to continue...");
  readline();
end

max_stars_repo_name: grasingerm/polymer-stats

max_issues_repo_name: grasingerm/polymer-stats

max_forks_repo_name: grasingerm/polymer-stats


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg
using Base.BinaryPlatforms

include("../../fancy_toys.jl")

name = "libclangex"
version = v"0.1.6"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/Gnimuc/libclangex.git", "527f22f606f43682c89ef76746d5d9652b2db1dc")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd libclangex/
mkdir build && cd build
cmake .. -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
     -DLLVM_DIR=$prefix \
     -DCLANG_DIR=$prefix \
     -DLLVM_ASSERT_BUILD=false \
     -DCMAKE_BUILD_TYPE=Release \
     -DCMAKE_EXPORT_COMPILE_COMMANDS=true
make -j${nproc}
make install
install_license ../COPYRIGHT ../LICENSE-APACHE ../LICENSE-MIT
"""

function configure(julia_version, llvm_version)
    # These are the platforms we will build for by default, unless further
    # platforms are passed in on the command line
    platforms = expand_cxxstring_abis(supported_platforms(; experimental=true))

    foreach(platforms) do p
        BinaryPlatforms.add_tag!(p.tags, "julia_version", string(julia_version))
    end

    # The products that we will ensure are always built
    products = Product[
        LibraryProduct("libclangex", :libclangex),
    ]
    # ver = "$(llvm_version.major).$(llvm_version.minor).$(llvm_version.patch)"
    dependencies = [
        RuntimeDependency("Clang_jll"),
        BuildDependency(get_addable_spec("LLVM_full_jll", llvm_version))
        # Dependency("LLVM_full_jll", llvm_version, compat=ver),
        # ^ is given through julia_version tag
    ]

    return platforms, products, dependencies
end

# TODO: Don't require build-id on LLVM version
supported = (
    # (v"1.6", v"11.0.1+3"),
    (v"1.7", v"12.0.1+3"),
    (v"1.8", v"12.0.1+3"),
)

for (julia_version, llvm_version) in supported
    platforms, products, dependencies = configure(julia_version, llvm_version)

    any(should_build_platform.(triplet.(platforms))) || continue

    # Build the tarballs.
    build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
                   preferred_gcc_version=v"8", julia_compat="1.7")
end

max_stars_repo_name: staticfloat/BinaryBuilderRepos

max_issues_repo_name: staticfloat/BinaryBuilderRepos

max_forks_repo_name: staticfloat/BinaryBuilderRepos


text:
test_path = mktempdir()
@testset "ThermalGen data misspecification" begin
    # See https://discourse.julialang.org/t/how-to-use-test-warn/15557/5 about testing for warning throwing
    info_message = "The data doesn't include devices of type ThermalStandard, consider changing the device models"
    device_model = DeviceModel(ThermalStandard, ThermalStandardUnitCommitment)
    c_sys5_re_only = PSB.build_system(PSITestSystems, "c_sys5_re_only")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5_re_only)
    @test_logs (:info, info_message) match_mode = :any mock_construct_device!(
        model,
        device_model,
    )
end

################################### Unit Commitment tests ##################################
@testset "Thermal UC With DC - PF" begin
    bin_variable_keys = [
        PSI.VariableKey(OnVariable, PSY.ThermalStandard),
        PSI.VariableKey(StartVariable, PSY.ThermalStandard),
        PSI.VariableKey(StopVariable, PSY.ThermalStandard),
    ]

    uc_constraint_keys = [
        PSI.ConstraintKey(RampConstraint, PSY.ThermalStandard, "up"),
        PSI.ConstraintKey(RampConstraint, PSY.ThermalStandard, "dn"),
        PSI.ConstraintKey(DurationConstraint, PSY.ThermalStandard, "up"),
        PSI.ConstraintKey(DurationConstraint, PSY.ThermalStandard, "dn"),
    ]

    aux_variables_keys = [
        PSI.AuxVarKey(PSI.TimeDurationOff, ThermalStandard),
        PSI.AuxVarKey(PSI.TimeDurationOn, ThermalStandard),
    ]
    device_model = DeviceModel(ThermalStandard, ThermalStandardUnitCommitment)

    c_sys5_uc = PSB.build_system(PSITestSystems, "c_sys5_uc")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5_uc)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 480, 0, 480, 120, 120, true)
    psi_constraint_test(model, uc_constraint_keys)
    psi_checkbinvar_test(model, bin_variable_keys)
    psi_checkobjfun_test(model, GAEVF)
    psi_aux_variable_test(model, aux_variables_keys)

    c_sys14 = PSB.build_system(PSITestSystems, "c_sys14")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys14)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 480, 0, 240, 120, 120, true)
    psi_checkbinvar_test(model, bin_variable_keys)
    psi_checkobjfun_test(model, GQEVF)
end

@testset "Thermal UC With AC - PF" begin
    bin_variable_keys = [
        PSI.VariableKey(OnVariable, PSY.ThermalStandard),
        PSI.VariableKey(StartVariable, PSY.ThermalStandard),
        PSI.VariableKey(StopVariable, PSY.ThermalStandard),
    ]
    uc_constraint_keys = [
        PSI.ConstraintKey(RampConstraint, PSY.ThermalStandard, "up"),
        PSI.ConstraintKey(RampConstraint, PSY.ThermalStandard, "dn"),
        PSI.ConstraintKey(DurationConstraint, PSY.ThermalStandard, "up"),
        PSI.ConstraintKey(DurationConstraint, PSY.ThermalStandard, "dn"),
    ]

    aux_variables_keys = [
        PSI.AuxVarKey(PSI.TimeDurationOff, ThermalStandard),
        PSI.AuxVarKey(PSI.TimeDurationOn, ThermalStandard),
    ]

    device_model = DeviceModel(ThermalStandard, ThermalStandardUnitCommitment)

    c_sys5_uc = PSB.build_system(PSITestSystems, "c_sys5_uc")
    model = DecisionModel(MockOperationProblem, ACPPowerModel, c_sys5_uc)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 600, 0, 600, 240, 120, true)
    psi_constraint_test(model, uc_constraint_keys)
    psi_checkbinvar_test(model, bin_variable_keys)
    psi_checkobjfun_test(model, GAEVF)
    psi_aux_variable_test(model, aux_variables_keys)

    c_sys14 = PSB.build_system(PSITestSystems, "c_sys14")
    model = DecisionModel(MockOperationProblem, ACPPowerModel, c_sys14;)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 600, 0, 360, 240, 120, true)
    psi_checkbinvar_test(model, bin_variable_keys)
    psi_checkobjfun_test(model, GQEVF)
end

@testset "Thermal MultiStart UC With DC - PF" begin
    bin_variable_keys = [
        PSI.VariableKey(OnVariable, PSY.ThermalMultiStart),
        PSI.VariableKey(StartVariable, PSY.ThermalMultiStart),
        PSI.VariableKey(StopVariable, PSY.ThermalMultiStart),
    ]
    uc_constraint_keys = [
        PSI.ConstraintKey(RampConstraint, PSY.ThermalMultiStart, "up"),
        PSI.ConstraintKey(RampConstraint, PSY.ThermalMultiStart, "dn"),
        PSI.ConstraintKey(DurationConstraint, PSY.ThermalMultiStart, "up"),
        PSI.ConstraintKey(DurationConstraint, PSY.ThermalMultiStart, "dn"),
    ]
    device_model = DeviceModel(ThermalMultiStart, ThermalStandardUnitCommitment)

    c_sys5_uc = PSB.build_system(PSITestSystems, "c_sys5_pglib")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5_uc;)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 384, 0, 240, 48, 144, true)
    psi_constraint_test(model, uc_constraint_keys)
    psi_checkbinvar_test(model, bin_variable_keys)
    psi_checkobjfun_test(model, GAEVF)
end

@testset "Thermal MultiStart UC With AC - PF" begin
    bin_variable_keys = [
        PSI.VariableKey(OnVariable, PSY.ThermalMultiStart),
        PSI.VariableKey(StartVariable, PSY.ThermalMultiStart),
        PSI.VariableKey(StopVariable, PSY.ThermalMultiStart),
    ]
    uc_constraint_keys = [
        PSI.ConstraintKey(RampConstraint, PSY.ThermalMultiStart, "up"),
        PSI.ConstraintKey(RampConstraint, PSY.ThermalMultiStart, "dn"),
        PSI.ConstraintKey(DurationConstraint, PSY.ThermalMultiStart, "up"),
        PSI.ConstraintKey(DurationConstraint, PSY.ThermalMultiStart, "dn"),
    ]
    device_model = DeviceModel(ThermalMultiStart, ThermalStandardUnitCommitment)

    c_sys5_uc = PSB.build_system(PSITestSystems, "c_sys5_pglib")
    model = DecisionModel(MockOperationProblem, ACPPowerModel, c_sys5_uc;)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 432, 0, 288, 96, 144, true)
    psi_constraint_test(model, uc_constraint_keys)
    psi_checkbinvar_test(model, bin_variable_keys)
    psi_checkobjfun_test(model, GAEVF)
end

################################### Basic Unit Commitment tests ############################
@testset "Thermal Basic UC With DC - PF" begin
    bin_variable_keys = [
        PSI.VariableKey(OnVariable, PSY.ThermalStandard),
        PSI.VariableKey(StartVariable, PSY.ThermalStandard),
        PSI.VariableKey(StopVariable, PSY.ThermalStandard),
    ]
    device_model = DeviceModel(ThermalStandard, ThermalBasicUnitCommitment)

    c_sys5_uc = PSB.build_system(PSITestSystems, "c_sys5_uc")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5_uc)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 480, 0, 240, 120, 120, true)
    psi_checkbinvar_test(model, bin_variable_keys)
    psi_checkobjfun_test(model, GAEVF)

    c_sys14 = PSB.build_system(PSITestSystems, "c_sys14")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys14;)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 480, 0, 240, 120, 120, true)
    psi_checkbinvar_test(model, bin_variable_keys)
    psi_checkobjfun_test(model, GQEVF)
end

@testset "Thermal Basic UC With AC - PF" begin
    bin_variable_keys = [
        PSI.VariableKey(OnVariable, PSY.ThermalStandard),
        PSI.VariableKey(StartVariable, PSY.ThermalStandard),
        PSI.VariableKey(StopVariable, PSY.ThermalStandard),
    ]
    device_model = DeviceModel(ThermalStandard, ThermalBasicUnitCommitment)

    c_sys5_uc = PSB.build_system(PSITestSystems, "c_sys5_uc")
    model = DecisionModel(MockOperationProblem, ACPPowerModel, c_sys5_uc)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 600, 0, 360, 240, 120, true)
    psi_checkbinvar_test(model, bin_variable_keys)
    psi_checkobjfun_test(model, GAEVF)

    c_sys14 = PSB.build_system(PSITestSystems, "c_sys14")
    model = DecisionModel(MockOperationProblem, ACPPowerModel, c_sys14;)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 600, 0, 360, 240, 120, true)
    psi_checkbinvar_test(model, bin_variable_keys)
    psi_checkobjfun_test(model, GQEVF)
end

@testset "Thermal MultiStart Basic UC With DC - PF" begin
    bin_variable_keys = [
        PSI.VariableKey(OnVariable, PSY.ThermalMultiStart),
        PSI.VariableKey(StartVariable, PSY.ThermalMultiStart),
        PSI.VariableKey(StopVariable, PSY.ThermalMultiStart),
    ]
    device_model = DeviceModel(ThermalMultiStart, ThermalBasicUnitCommitment)

    c_sys5_uc = PSB.build_system(PSITestSystems, "c_sys5_pglib")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5_uc;)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 384, 0, 96, 48, 144, true)
    psi_checkbinvar_test(model, bin_variable_keys)
    psi_checkobjfun_test(model, GAEVF)
end

@testset "Thermal MultiStart Basic UC With AC - PF" begin
    bin_variable_keys = [
        PSI.VariableKey(OnVariable, PSY.ThermalMultiStart),
        PSI.VariableKey(StartVariable, PSY.ThermalMultiStart),
        PSI.VariableKey(StopVariable, PSY.ThermalMultiStart),
    ]
    device_model = DeviceModel(ThermalMultiStart, ThermalBasicUnitCommitment)

    c_sys5_uc = PSB.build_system(PSITestSystems, "c_sys5_pglib")
    model = DecisionModel(MockOperationProblem, ACPPowerModel, c_sys5_uc;)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 432, 0, 144, 96, 144, true)
    psi_checkbinvar_test(model, bin_variable_keys)
    psi_checkobjfun_test(model, GAEVF)
end

################################### Basic Dispatch tests ###################################
@testset "ThermalStandard with ThermalBasicDispatch With DC - PF" begin
    device_model = DeviceModel(ThermalStandard, ThermalBasicDispatch)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 120, 0, 120, 120, 0, false)
    psi_checkobjfun_test(model, GAEVF)

    c_sys14 = PSB.build_system(PSITestSystems, "c_sys14")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys14)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 120, 0, 120, 120, 0, false)
    psi_checkobjfun_test(model, GQEVF)
end

@testset "ThermalStandard  with ThermalBasicDispatch With AC - PF" begin
    device_model = DeviceModel(ThermalStandard, ThermalBasicDispatch)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")

    model = DecisionModel(MockOperationProblem, ACPPowerModel, c_sys5)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 240, 0, 240, 240, 0, false)
    psi_checkobjfun_test(model, GAEVF)

    c_sys14 = PSB.build_system(PSITestSystems, "c_sys14")
    model = DecisionModel(MockOperationProblem, ACPPowerModel, c_sys14;)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 240, 0, 240, 240, 0, false)
    psi_checkobjfun_test(model, GQEVF)
end

# This Formulation is currently broken
@testset "ThermalMultiStart with ThermalBasicDispatch With DC - PF" begin
    device_model = DeviceModel(ThermalMultiStart, ThermalBasicDispatch)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5_pglib")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 240, 0, 48, 48, 96, false)
    psi_checkobjfun_test(model, GAEVF)
end

@testset "ThermalMultiStart with ThermalBasicDispatch With AC - PF" begin
    device_model = DeviceModel(ThermalMultiStart, ThermalBasicDispatch)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5_pglib")
    model = DecisionModel(MockOperationProblem, ACPPowerModel, c_sys5)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 288, 0, 96, 96, 96, false)
    psi_checkobjfun_test(model, GAEVF)
end

################################### No Minimum Dispatch tests ##############################
@testset "Thermal Dispatch NoMin With DC - PF" begin
    device_model = DeviceModel(ThermalStandard, ThermalDispatchNoMin)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 120, 0, 120, 120, 0, false)
    key = PSI.ConstraintKey(ActivePowerVariableLimitsConstraint, ThermalStandard, "lb")
    moi_lbvalue_test(model, key, 0.0)
    psi_checkobjfun_test(model, GAEVF)

    c_sys14 = PSB.build_system(PSITestSystems, "c_sys14")

    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys14)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 120, 0, 120, 120, 0, false)
    key = PSI.ConstraintKey(ActivePowerVariableLimitsConstraint, ThermalStandard, "lb")
    moi_lbvalue_test(model, key, 0.0)
    psi_checkobjfun_test(model, GQEVF)
end

@testset "Thermal Dispatch NoMin With AC - PF" begin
    device_model = DeviceModel(ThermalStandard, ThermalDispatchNoMin)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    model = DecisionModel(MockOperationProblem, ACPPowerModel, c_sys5)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 240, 0, 240, 240, 0, false)
    key = PSI.ConstraintKey(ActivePowerVariableLimitsConstraint, ThermalStandard, "lb")
    moi_lbvalue_test(model, key, 0.0)
    psi_checkobjfun_test(model, GAEVF)

    c_sys14 = PSB.build_system(PSITestSystems, "c_sys14")

    model = DecisionModel(MockOperationProblem, ACPPowerModel, c_sys14;)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 240, 0, 240, 240, 0, false)
    key = PSI.ConstraintKey(ActivePowerVariableLimitsConstraint, ThermalStandard, "lb")
    moi_lbvalue_test(model, key, 0.0)
    psi_checkobjfun_test(model, GQEVF)
end

@testset "Thermal Dispatch NoMin With DC - PF" begin
    device_model = DeviceModel(ThermalMultiStart, ThermalDispatchNoMin)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5_pglib")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5)
    @test_throws IS.ConflictingInputsError mock_construct_device!(model, device_model)
end

@testset "ThermalMultiStart Dispatch NoMin With AC - PF" begin
    device_model = DeviceModel(ThermalMultiStart, ThermalDispatchNoMin)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5_pglib")
    model = DecisionModel(MockOperationProblem, ACPPowerModel, c_sys5;)
    @test_throws IS.ConflictingInputsError mock_construct_device!(model, device_model)
end

@testset "Operation Model ThermalDispatchNoMin - and PWL Non Convex" begin
    c_sys5_pwl_ed_nonconvex = PSB.build_system(PSITestSystems, "c_sys5_pwl_ed_nonconvex")
    template = get_thermal_dispatch_template_network()
    set_device_model!(template, DeviceModel(ThermalStandard, ThermalDispatchNoMin))
    model = DecisionModel(
        MockOperationProblem,
        CopperPlatePowerModel,
        c_sys5_pwl_ed_nonconvex;
        export_pwl_vars=true,
        initialize_model=false,
    )
    @test_throws IS.InvalidValue mock_construct_device!(
        model,
        DeviceModel(ThermalStandard, ThermalDispatchNoMin),
    )
end

################################## Ramp Limited Testing ##################################
@testset "ThermalStandard with ThermalStandardDispatch With DC - PF" begin
    constraint_keys = [
        PSI.ConstraintKey(RampConstraint, PSY.ThermalStandard, "up"),
        PSI.ConstraintKey(RampConstraint, PSY.ThermalStandard, "dn"),
    ]
    device_model = DeviceModel(ThermalStandard, ThermalStandardDispatch)
    c_sys5_uc = PSB.build_system(PSITestSystems, "c_sys5_uc")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5_uc;)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 120, 0, 168, 168, 0, false)
    psi_constraint_test(model, constraint_keys)
    psi_checkobjfun_test(model, GAEVF)

    c_sys14 = PSB.build_system(PSITestSystems, "c_sys14")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys14;)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 120, 0, 120, 120, 0, false)
    psi_checkobjfun_test(model, GQEVF)
end

@testset "ThermalStandard with ThermalStandardDispatch With AC - PF" begin
    constraint_keys = [
        PSI.ConstraintKey(RampConstraint, PSY.ThermalStandard, "up"),
        PSI.ConstraintKey(RampConstraint, PSY.ThermalStandard, "dn"),
    ]
    device_model = DeviceModel(ThermalStandard, ThermalStandardDispatch)
    c_sys5_uc = PSB.build_system(PSITestSystems, "c_sys5_uc")
    model = DecisionModel(MockOperationProblem, ACPPowerModel, c_sys5_uc;)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 240, 0, 288, 288, 0, false)
    psi_constraint_test(model, constraint_keys)
    psi_checkobjfun_test(model, GAEVF)

    c_sys14 = PSB.build_system(PSITestSystems, "c_sys14")
    model = DecisionModel(MockOperationProblem, ACPPowerModel, c_sys14;)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 240, 0, 240, 240, 0, false)
    psi_checkobjfun_test(model, GQEVF)
end

@testset "ThermalMultiStart with ThermalStandardDispatch With DC - PF" begin
    constraint_keys = [
        PSI.ConstraintKey(RampConstraint, PSY.ThermalMultiStart, "up"),
        PSI.ConstraintKey(RampConstraint, PSY.ThermalMultiStart, "dn"),
    ]
    device_model = DeviceModel(ThermalMultiStart, ThermalStandardDispatch)
    c_sys5_uc = PSB.build_system(PSITestSystems, "c_sys5_pglib")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5_uc;)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 240, 0, 96, 96, 96, false)
    psi_constraint_test(model, constraint_keys)
    psi_checkobjfun_test(model, GAEVF)
end

@testset "ThermalMultiStart with ThermalStandardDispatch With AC - PF" begin
    constraint_keys = [
        PSI.ConstraintKey(RampConstraint, PSY.ThermalMultiStart, "up"),
        PSI.ConstraintKey(RampConstraint, PSY.ThermalMultiStart, "dn"),
    ]
    device_model = DeviceModel(ThermalMultiStart, ThermalStandardDispatch)
    c_sys5_uc = PSB.build_system(PSITestSystems, "c_sys5_pglib")
    model = DecisionModel(MockOperationProblem, ACPPowerModel, c_sys5_uc;)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 288, 0, 144, 144, 96, false)
    psi_constraint_test(model, constraint_keys)
    psi_checkobjfun_test(model, GAEVF)
end

################################### ThermalMultiStart Testing ##############################

@testset "Thermal MultiStart with MultiStart UC and DC - PF" begin
    constraint_keys = [
        PSI.ConstraintKey(ActiveRangeICConstraint, PSY.ThermalMultiStart),
        PSI.ConstraintKey(StartTypeConstraint, PSY.ThermalMultiStart),
        PSI.ConstraintKey(MustRunConstraint, PSY.ThermalMultiStart),
        PSI.ConstraintKey(
            StartupTimeLimitTemperatureConstraint,
            PSY.ThermalMultiStart,
            "warm",
        ),
        PSI.ConstraintKey(
            StartupTimeLimitTemperatureConstraint,
            PSY.ThermalMultiStart,
            "hot",
        ),
        PSI.ConstraintKey(StartupInitialConditionConstraint, PSY.ThermalMultiStart, "lb"),
        PSI.ConstraintKey(StartupInitialConditionConstraint, PSY.ThermalMultiStart, "ub"),
    ]
    device_model = DeviceModel(PSY.ThermalMultiStart, PSI.ThermalMultiStartUnitCommitment)
    no_less_than = Dict(true => 334, false => 282)
    c_sys5_pglib = PSB.build_system(PSITestSystems, "c_sys5_pglib")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5_pglib;)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 528, 0, no_less_than[false], 108, 192, true)
    psi_constraint_test(model, constraint_keys)
    psi_checkobjfun_test(model, GAEVF)
end

@testset "Thermal MultiStart with MultiStart UC and AC - PF" begin
    constraint_keys = [
        PSI.ConstraintKey(ActiveRangeICConstraint, PSY.ThermalMultiStart),
        PSI.ConstraintKey(StartTypeConstraint, PSY.ThermalMultiStart),
        PSI.ConstraintKey(MustRunConstraint, PSY.ThermalMultiStart),
        PSI.ConstraintKey(
            StartupTimeLimitTemperatureConstraint,
            PSY.ThermalMultiStart,
            "warm",
        ),
        PSI.ConstraintKey(
            StartupTimeLimitTemperatureConstraint,
            PSY.ThermalMultiStart,
            "hot",
        ),
        PSI.ConstraintKey(StartupInitialConditionConstraint, PSY.ThermalMultiStart, "lb"),
        PSI.ConstraintKey(StartupInitialConditionConstraint, PSY.ThermalMultiStart, "ub"),
    ]
    device_model = DeviceModel(PSY.ThermalMultiStart, PSI.ThermalMultiStartUnitCommitment)
    no_less_than = Dict(true => 382, false => 330)
    c_sys5_pglib = PSB.build_system(PSITestSystems, "c_sys5_pglib")
    model = DecisionModel(MockOperationProblem, ACPPowerModel, c_sys5_pglib;)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 576, 0, no_less_than[false], 156, 192, true)
    psi_constraint_test(model, constraint_keys)
    psi_checkobjfun_test(model, GAEVF)
end

################################ Thermal Compact UC Testing ################################
@testset "Thermal Standard with Compact UC and DC - PF" begin
    device_model = DeviceModel(PSY.ThermalStandard, PSI.ThermalCompactUnitCommitment)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 480, 0, 480, 120, 120, true)
    psi_checkobjfun_test(model, GAEVF)
end

@testset "Thermal MultiStart with Compact UC and DC - PF" begin
    device_model = DeviceModel(PSY.ThermalMultiStart, PSI.ThermalCompactUnitCommitment)
    c_sys5_pglib = PSB.build_system(PSITestSystems, "c_sys5_pglib")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5_pglib;)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 384, 0, 240, 48, 144, true)
    psi_checkobjfun_test(model, GAEVF)
end

@testset "Thermal Standard with Compact UC and AC - PF" begin
    device_model = DeviceModel(PSY.ThermalStandard, PSI.ThermalCompactUnitCommitment)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    model = DecisionModel(MockOperationProblem, ACPPowerModel, c_sys5)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 600, 0, 600, 240, 120, true)
    psi_checkobjfun_test(model, GAEVF)
end

@testset "Thermal MultiStart with Compact UC and AC - PF" begin
    device_model = DeviceModel(PSY.ThermalMultiStart, PSI.ThermalCompactUnitCommitment)
    c_sys5_pglib = PSB.build_system(PSITestSystems, "c_sys5_pglib")
    model = DecisionModel(MockOperationProblem, ACPPowerModel, c_sys5_pglib;)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 432, 0, 288, 96, 144, true)
    psi_checkobjfun_test(model, GAEVF)
end

################################ Thermal Basic Compact UC Testing ################################
@testset "Thermal Standard with Compact UC and DC - PF" begin
    device_model = DeviceModel(PSY.ThermalStandard, PSI.ThermalBasicCompactUnitCommitment)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 480, 0, 240, 120, 120, true)
    psi_checkobjfun_test(model, GAEVF)
end

@testset "Thermal MultiStart with Compact UC and DC - PF" begin
    device_model = DeviceModel(PSY.ThermalMultiStart, PSI.ThermalBasicCompactUnitCommitment)
    c_sys5_pglib = PSB.build_system(PSITestSystems, "c_sys5_pglib")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5_pglib;)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 384, 0, 96, 48, 144, true)
    psi_checkobjfun_test(model, GAEVF)
end

@testset "Thermal Standard with Compact UC and AC - PF" begin
    device_model = DeviceModel(PSY.ThermalStandard, PSI.ThermalBasicCompactUnitCommitment)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    model = DecisionModel(MockOperationProblem, ACPPowerModel, c_sys5)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 600, 0, 360, 240, 120, true)
    psi_checkobjfun_test(model, GAEVF)
end

@testset "Thermal MultiStart with Compact UC and AC - PF" begin
    device_model = DeviceModel(PSY.ThermalMultiStart, PSI.ThermalBasicCompactUnitCommitment)
    c_sys5_pglib = PSB.build_system(PSITestSystems, "c_sys5_pglib")
    model = DecisionModel(MockOperationProblem, ACPPowerModel, c_sys5_pglib;)
    mock_construct_device!(model, device_model)
    moi_tests(model, false, 432, 0, 144, 96, 144, true)
    psi_checkobjfun_test(model, GAEVF)
end

############################ Thermal Compact Dispatch Testing ##############################
@testset "Thermal Standard with Compact Dispatch and DC - PF" begin
    device_model = DeviceModel(PSY.ThermalStandard, PSI.ThermalCompactDispatch)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5)
    mock_construct_device!(model, device_model; built_for_recurrent_solves=true)
    moi_tests(model, true, 120, 0, 144, 144, 0, false)
    psi_checkobjfun_test(model, GAEVF)
end

@testset "Thermal MultiStart with Compact Dispatch and DC - PF" begin
    device_model = DeviceModel(PSY.ThermalMultiStart, PSI.ThermalCompactDispatch)
    c_sys5_pglib = PSB.build_system(PSITestSystems, "c_sys5_pglib")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5_pglib)
    mock_construct_device!(model, device_model; built_for_recurrent_solves=true)
    moi_tests(model, true, 240, 0, 96, 96, 96, false)
    psi_checkobjfun_test(model, GAEVF)
end

@testset "Thermal Standard with Compact Dispatch and AC - PF" begin
    device_model = DeviceModel(PSY.ThermalStandard, PSI.ThermalCompactDispatch)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    model = DecisionModel(MockOperationProblem, ACPPowerModel, c_sys5)
    mock_construct_device!(model, device_model; built_for_recurrent_solves=true)
    moi_tests(model, true, 240, 0, 264, 264, 0, false)
    psi_checkobjfun_test(model, GAEVF)
end

@testset "Thermal MultiStart with Compact Dispatch and AC - PF" begin
    device_model = DeviceModel(PSY.ThermalMultiStart, PSI.ThermalCompactDispatch)
    c_sys5_pglib = PSB.build_system(PSITestSystems, "c_sys5_pglib")
    model = DecisionModel(MockOperationProblem, ACPPowerModel, c_sys5_pglib)
    mock_construct_device!(model, device_model; built_for_recurrent_solves=true)
    moi_tests(model, true, 288, 0, 144, 144, 96, false)
    psi_checkobjfun_test(model, GAEVF)
end

############################# Model validation tests #######################################
@testset "Solving ED with CopperPlate for testing Ramping Constraints" begin
    ramp_test_sys = PSB.build_system(PSITestSystems, "c_ramp_test")
    template = ProblemTemplate(CopperPlatePowerModel)
    set_device_model!(template, ThermalStandard, ThermalStandardDispatch)
    set_device_model!(template, PowerLoad, StaticPowerLoad)
    ED = DecisionModel(
        EconomicDispatchProblem,
        template,
        ramp_test_sys;
        optimizer=HiGHS_optimizer,
        initialize_model=false,
    )
    @test build!(ED; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
    moi_tests(ED, false, 10, 0, 15, 15, 5, false)
    psi_checksolve_test(ED, [MOI.OPTIMAL], 11191.00)
end

# Testing Duration Constraints
@testset "Solving UC with CopperPlate for testing Duration Constraints" begin
    template = get_thermal_standard_uc_template()
    UC = DecisionModel(
        UnitCommitmentProblem,
        template,
        PSB.build_system(PSITestSystems, "c_duration_test");
        optimizer=HiGHS_optimizer,
        initialize_model=false,
    )
    @test build!(UC; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
    moi_tests(UC, false, 56, 0, 56, 14, 21, true)
    psi_checksolve_test(UC, [MOI.OPTIMAL], 8223.50)
end

## PWL linear Cost implementation test
@testset "Solving UC with CopperPlate testing Convex PWL" begin
    template = get_thermal_standard_uc_template()
    UC = DecisionModel(
        UnitCommitmentProblem,
        template,
        PSB.build_system(PSITestSystems, "c_linear_pwl_test");
        optimizer=HiGHS_optimizer,
        initialize_model=false,
    )
    @test build!(UC; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
    moi_tests(UC, false, 32, 0, 8, 4, 14, true)
    psi_checksolve_test(UC, [MOI.OPTIMAL], 9336.736919354838)
end

@testset "Solving UC with CopperPlate testing PWL-SOS2 implementation" begin
    template = get_thermal_standard_uc_template()
    UC = DecisionModel(
        UnitCommitmentProblem,
        template,
        PSB.build_system(PSITestSystems, "c_sos_pwl_test");
        optimizer=cbc_optimizer,
        initialize_model=false,
    )
    @test build!(UC; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
    moi_tests(UC, false, 32, 0, 8, 4, 14, true)
    # Cbc can have reliability issues with SoS. The objective function target in the this
    # test was calculated with CPLEX do not change if Cbc gets a bad result
    psi_checksolve_test(UC, [MOI.OPTIMAL], 8500.0, 10.0)
end

@testset "UC with MarketBid Cost in ThermalGenerators" begin
    template = get_thermal_standard_uc_template()
    set_device_model!(
        template,
        DeviceModel(ThermalMultiStart, ThermalMultiStartUnitCommitment),
    )
    UC = DecisionModel(
        UnitCommitmentProblem,
        template,
        PSB.build_system(PSITestSystems, "c_market_bid_cost");
        optimizer=cbc_optimizer,
        initialize_model=false,
    )
    @test build!(UC; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
    moi_tests(UC, false, 38, 0, 16, 8, 16, true)
end

@testset "Solving UC Models with Linear Networks" begin
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    c_sys5_dc = PSB.build_system(PSITestSystems, "c_sys5_dc")
    systems = [c_sys5, c_sys5_dc]
    networks = [DCPPowerModel, NFAPowerModel, StandardPTDFModel, CopperPlatePowerModel]
    commitment_models = [ThermalStandardUnitCommitment, ThermalCompactUnitCommitment]
    PTDF_ref = IdDict{System, PTDF}(c_sys5 => PTDF(c_sys5), c_sys5_dc => PTDF(c_sys5_dc))

    for net in networks, sys in systems, model in commitment_models
        template =
            get_thermal_dispatch_template_network(NetworkModel(net, PTDF=PTDF_ref[sys]))
        set_device_model!(template, ThermalStandard, model)
        UC = DecisionModel(template, sys; optimizer=HiGHS_optimizer)
        @test build!(UC; output_dir=mktempdir(cleanup=true)) == PSI.BuildStatus.BUILT
        psi_checksolve_test(UC, [MOI.OPTIMAL, MOI.LOCALLY_SOLVED], 340000, 100000)
    end
end

@testset "Test Feedforwards to ThermalStandard with ThermalStandardDispatch" begin
    device_model = DeviceModel(ThermalStandard, ThermalStandardDispatch)
    ff_sc = SemiContinuousFeedforward(
        component_type=ThermalStandard,
        source=OnVariable,
        affected_values=[ActivePowerVariable],
    )

    ff_ub = UpperBoundFeedforward(
        component_type=ThermalStandard,
        source=ActivePowerVariable,
        affected_values=[ActivePowerVariable],
    )

    PSI.attach_feedforward!(device_model, ff_sc)
    PSI.attach_feedforward!(device_model, ff_ub)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5)
    mock_construct_device!(model, device_model; built_for_recurrent_solves=true)
    moi_tests(model, true, 120, 0, 264, 144, 0, false)
end

@testset "Test Feedforwards to ThermalStandard with ThermalBasicDispatch" begin
    device_model = DeviceModel(ThermalStandard, ThermalBasicDispatch)
    ff_sc = SemiContinuousFeedforward(
        component_type=ThermalStandard,
        source=OnVariable,
        affected_values=[ActivePowerVariable],
    )

    ff_ub = UpperBoundFeedforward(
        component_type=ThermalStandard,
        source=ActivePowerVariable,
        affected_values=[ActivePowerVariable],
    )

    PSI.attach_feedforward!(device_model, ff_sc)
    PSI.attach_feedforward!(device_model, ff_ub)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5)
    mock_construct_device!(model, device_model; built_for_recurrent_solves=true)
    moi_tests(model, true, 120, 0, 240, 120, 0, false)
end

@testset "Test Feedforwards to ThermalStandard with ThermalCompactDispatch" begin
    device_model = DeviceModel(PSY.ThermalStandard, PSI.ThermalCompactDispatch)
    ff_sc = SemiContinuousFeedforward(
        component_type=ThermalStandard,
        source=OnVariable,
        affected_values=[PowerAboveMinimumVariable],
    )

    ff_ub = UpperBoundFeedforward(
        component_type=ThermalStandard,
        source=PSI.PowerAboveMinimumVariable,
        affected_values=[PSI.PowerAboveMinimumVariable],
    )

    PSI.attach_feedforward!(device_model, ff_sc)
    PSI.attach_feedforward!(device_model, ff_ub)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5)
    mock_construct_device!(model, device_model; built_for_recurrent_solves=true)
    moi_tests(model, true, 120, 0, 264, 144, 0, false)
end

@testset "Test Feedforwards to ThermalMultiStart with ThermalStandardDispatch" begin
    device_model = DeviceModel(ThermalMultiStart, ThermalStandardDispatch)
    ff_sc = SemiContinuousFeedforward(
        component_type=ThermalMultiStart,
        source=OnVariable,
        affected_values=[ActivePowerVariable],
    )

    ff_ub = UpperBoundFeedforward(
        component_type=ThermalMultiStart,
        source=ActivePowerVariable,
        affected_values=[ActivePowerVariable],
    )

    PSI.attach_feedforward!(device_model, ff_sc)
    PSI.attach_feedforward!(device_model, ff_ub)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5_pglib")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5)
    mock_construct_device!(model, device_model; built_for_recurrent_solves=true)
    moi_tests(model, true, 240, 0, 144, 96, 96, false)
end

@testset "Test Feedforwards to ThermalMultiStart with ThermalBasicDispatch" begin
    device_model = DeviceModel(ThermalMultiStart, ThermalBasicDispatch)
    ff_sc = SemiContinuousFeedforward(
        component_type=ThermalMultiStart,
        source=OnVariable,
        affected_values=[ActivePowerVariable],
    )

    ff_ub = UpperBoundFeedforward(
        component_type=ThermalMultiStart,
        source=ActivePowerVariable,
        affected_values=[ActivePowerVariable],
    )

    PSI.attach_feedforward!(device_model, ff_sc)
    PSI.attach_feedforward!(device_model, ff_ub)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5_pglib")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5)
    mock_construct_device!(model, device_model; built_for_recurrent_solves=true)
    moi_tests(model, true, 240, 0, 96, 48, 96, false)
end

@testset "Test Feedforwards to ThermalMultiStart with ThermalCompactDispatch" begin
    device_model = DeviceModel(PSY.ThermalMultiStart, PSI.ThermalCompactDispatch)
    ff_sc = SemiContinuousFeedforward(
        component_type=ThermalMultiStart,
        source=OnVariable,
        affected_values=[PSI.PowerAboveMinimumVariable],
    )

    ff_ub = UpperBoundFeedforward(
        component_type=ThermalMultiStart,
        source=PSI.PowerAboveMinimumVariable,
        affected_values=[PSI.PowerAboveMinimumVariable],
    )

    PSI.attach_feedforward!(device_model, ff_sc)
    PSI.attach_feedforward!(device_model, ff_ub)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5_pglib")
    model = DecisionModel(MockOperationProblem, DCPPowerModel, c_sys5)
    mock_construct_device!(model, device_model; built_for_recurrent_solves=true)
    moi_tests(model, true, 240, 0, 144, 96, 96, false)
end


text:
mutable struct ProjMPO_MPS
  PH::ProjMPO
  pm::Vector{ProjMPS}
  weight::Float64
end

function ProjMPO_MPS(H::MPO, mpsv::Vector{MPS}; weight=1.0)
  return ProjMPO_MPS(ProjMPO(H), [ProjMPS(m) for m in mpsv], weight)
end

ProjMPO_MPS(H::MPO, Ms::MPS...; weight=1.0) = ProjMPO_MPS(H, [Ms...], weight)

nsite(P::ProjMPO_MPS) = nsite(P.PH)

Base.length(P::ProjMPO_MPS) = length(P.PH)

function product(P::ProjMPO_MPS, v::ITensor)::ITensor
  Pv = product(P.PH, v)
  for p in P.pm
    Pv += P.weight * product(p, v)
  end
  return Pv
end

function Base.eltype(P::ProjMPO_MPS)
  elT = eltype(P.PH)
  for p in P.pm
    elT = promote_type(elT, eltype(p))
  end
  return elT
end

(P::ProjMPO_MPS)(v::ITensor) = product(P, v)

Base.size(P::ProjMPO_MPS) = size(P.H)

function position!(P::ProjMPO_MPS, psi::MPS, pos::Int)
  position!(P.PH, psi, pos)
  for p in P.pm
    position!(p, psi, pos)
  end
end

noiseterm(P::ProjMPO_MPS, phi::ITensor, dir::String) = noiseterm(P.PH, phi, dir)