text:
# ----------------------------------------------------------------------------
# This file was autogenerated by /dev/parse-SeeK.jl.
# It contains data included in the SeeK-path package (github.com/giovannipizzi/seekpath).
# The SeeK-path package is licensed under the MIT license, a copy of which is
# included below.
# ----------------------------------------------------------------------------

const pathsd_3d = Dict(
   :hR1 => [[:Γ, :T, :H₂], [:H₀, :L, :Γ, :S₀], [:S₂, :F, :Γ]],
   :tI2 => [[:Γ, :X, :P, :N, :Γ, :M, :S], [:S₀, :Γ], [:X, :R], [:G, :M]],
   :tI1 => [[:Γ, :X, :M, :Γ, :Z], [:Z₀, :M], [:X, :P, :N, :Γ]],
   :cP1 => [[:Γ, :X, :M, :Γ, :R, :X], [:R, :M, :X₁]],
   :oC2 => [[:Γ, :Y, :F₀], [:Δ₀, :Γ, :Z, :B₀], [:G₀, :T, :Y], [:Γ, :S, :R, :Z, :T]],
   :oI3 => [[:Γ, :X, :F₀], [:Σ₀, :Γ, :Λ₀], [:G₀, :X], [:Γ, :R, :W, :S, :Γ, :T, :W]],
   :oP1 => [[:Γ, :X, :S, :Y, :Γ, :Z, :U, :R, :T, :Z], [:X, :U], [:Y, :T], [:S, :R]],
   :oA1 => [[:Γ, :Y, :C₀], [:Σ₀, :Γ, :Z, :A₀], [:E₀, :T, :Y], [:Γ, :S, :R, :Z, :T]],
   :mC1 => [[:Γ, :C], [:C₂, :Y₂, :Γ, :M₂, :D], [:D₂, :A, :Γ], [:L₂, :Γ, :V₂]],
   :oI2 => [[:Γ, :X, :U₂], [:Y₀, :Γ, :Λ₀], [:G₂, :X], [:Γ, :R, :W, :S, :Γ, :T, :W]],
   :hP1 => [[:Γ, :M, :K, :Γ, :A, :L, :H, :A], [:L, :M], [:H, :K, :H₂]],
   :tP1 => [[:Γ, :X, :M, :Γ, :Z, :R, :A, :Z], [:X, :R], [:M, :A]],
   :cI1 => [[:Γ, :H, :N, :Γ, :P, :H], [:P, :N]],
   :mC3 => [[:Γ, :A, :I₂], [:I, :M₂, :Γ, :Y], [:L₂, :Γ, :V₂]],
   :hP2 => [[:Γ, :M, :K, :Γ, :A, :L, :H, :A], [:L, :M], [:H, :K]],
   :hR2 => [[:Γ, :L, :T, :P₀], [:P₂, :Γ, :F]],
   :cF1 => [[:Γ, :X, :U], [:K, :Γ, :L, :W, :X, :W₂]],
   :oA2 => [[:Γ, :Y, :F₀], [:Δ₀, :Γ, :Z, :B₀], [:G₀, :T, :Y], [:Γ, :S, :R, :Z, :T]],
   :cP2 => [[:Γ, :X, :M, :Γ, :R, :X], [:R, :M]],
   :aP2 => [[:Γ, :X], [:Y, :Γ, :Z], [:R, :Γ, :T], [:U, :Γ, :V]],
   :oF1 => [[:Γ, :Y, :T, :Z, :Γ, :Σ₀], [:U₀, :T], [:Y, :C₀], [:A₀, :Z], [:Γ, :L]],
   :oC1 => [[:Γ, :Y, :C₀], [:Σ₀, :Γ, :Z, :A₀], [:E₀, :T, :Y], [:Γ, :S, :R, :Z, :T]],
   :aP3 => [[:Γ, :X], [:Y, :Γ, :Z], [:R₂, :Γ, :T₂], [:U₂, :Γ, :V₂]],
   :mC2 => [[:Γ, :Y, :M, :A, :Γ], [:L₂, :Γ, :V₂]],
   :oI1 => [[:Γ, :X, :F₂], [:Σ₀, :Γ, :Y₀], [:U₀, :X], [:Γ, :R, :W, :S, :Γ, :T, :W]],
   :mP1 => [[:Γ, :Z, :D, :B, :Γ, :A, :E, :Z, :C₂, :Y₂, :Γ]],
   :cF2 => [[:Γ, :X, :U], [:K, :Γ, :L, :W, :X]],
   :oF2 => [[:Γ, :T, :Z, :Y, :Γ, :Λ₀], [:Q₀, :Z], [:T, :G₀], [:H₀, :Y], [:Γ, :L]],
   :oF3 => [[:Γ, :Y, :C₀], [:A₀, :Z, :B₀], [:D₀, :T, :G₀], [:H₀, :Y], [:T, :Γ, :Z], [:Γ, :L]],
)

const pointsd_3d = Dict(
   :hR1 => [:Γ => :([0, 0, 0]), :T => :([1 / 2, 1 / 2, 1 / 2]), :L => :([1 / 2, 0, 0]), :F => :([1 / 2, 0, 1 / 2]), :S₀ => :([N, -N, 0]), :S₂ => :([1 - N, 0, N]), :H₀ => :([1 / 2, -1 + Y, 1 - Y]), :H₂ => :([Y, 1 - Y, 1 / 2])],
   :tI2 => [:Γ => :([0, 0, 0]), :M => :([1 / 2, 1 / 2, -1 / 2]), :X => :([0, 0, 1 / 2]), :P => :([1 / 4, 1 / 4, 1 / 4]), :N => :([0, 1 / 2, 0]), :S₀ => :([-H, H, H]), :S => :([H, 1 - H, -H]), :R => :([-Z, Z, 1 / 2]), :G => :([1 / 2, 1 / 2, -Z])],
   :tI1 => [:Γ => :([0, 0, 0]), :M => :([-1 / 2, 1 / 2, 1 / 2]), :X => :([0, 0, 1 / 2]), :P => :([1 / 4, 1 / 4, 1 / 4]), :Z => :([H, H, -H]), :Z₀ => :([-H, 1 - H, H]), :N => :([0, 1 / 2, 0])],
   :cP1 => [:Γ => :([0, 0, 0]), :R => :([1 / 2, 1 / 2, 1 / 2]), :M => :([1 / 2, 1 / 2, 0]), :X => :([0, 1 / 2, 0]), :X₁ => :([1 / 2, 0, 0])],
   :oC2 => [:Γ => :([0, 0, 0]), :Y => :([1 / 2, 1 / 2, 0]), :T => :([1 / 2, 1 / 2, 1 / 2]), :Z => :([0, 0, 1 / 2]), :S => :([0, 1 / 2, 0]), :R => :([0, 1 / 2, 1 / 2]), :Δ₀ => :([-X, X, 0]), :F₀ => :([X, 1 - X, 0]), :B₀ => :([-X, X, 1 / 2]), :G₀ => :([X, 1 - X, 1 / 2])],
   :oI3 => [:Γ => :([0, 0, 0]), :X => :([1 / 2, -1 / 2, 1 / 2]), :S => :([1 / 2, 0, 0]), :R => :([0, 1 / 2, 0]), :T => :([0, 0, 1 / 2]), :W => :([1 / 4, 1 / 4, 1 / 4]), :Σ₀ => :([-Y, Y, Y]), :F₀ => :([Y, -Y, 1 - Y]), :Λ₀ => :([Z, Z, -Z]), :G₀ => :([1 - Z, -Z, Z])],
   :oP1 => [:Γ => :([0, 0, 0]), :X => :([1 / 2, 0, 0]), :Z => :([0, 0, 1 / 2]), :U => :([1 / 2, 0, 1 / 2]), :Y => :([0, 1 / 2, 0]), :S => :([1 / 2, 1 / 2, 0]), :T => :([0, 1 / 2, 1 / 2]), :R => :([1 / 2, 1 / 2, 1 / 2])],
   :oA1 => [:Γ => :([0, 0, 0]), :Y => :([-1 / 2, 1 / 2, 0]), :T => :([-1 / 2, 1 / 2, 1 / 2]), :Z => :([0, 0, 1 / 2]), :S => :([0, 1 / 2, 0]), :R => :([0, 1 / 2, 1 / 2]), :Σ₀ => :([X, X, 0]), :C₀ => :([-X, 1 - X, 0]), :A₀ => :([X, X, 1 / 2]), :E₀ => :([-X, 1 - X, 1 / 2])],
   :mC1 => [:Γ => :([0, 0, 0]), :Y₂ => :([-1 / 2, 1 / 2, 0]), :A => :([0, 0, 1 / 2]), :M₂ => :([-1 / 2, 1 / 2, 1 / 2]), :V₂ => :([0, 1 / 2, 0]), :L₂ => :([0, 1 / 2, 1 / 2]), :C => :([1 - S, 1 - S, 0]), :C₂ => :([-1 + S, S, 0]), :D => :([-1 + P, P, 1 / 2]), :D₂ => :([1 - P, 1 - P, 1 / 2])],
   :oI2 => [:Γ => :([0, 0, 0]), :X => :([-1 / 2, 1 / 2, 1 / 2]), :S => :([1 / 2, 0, 0]), :R => :([0, 1 / 2, 0]), :T => :([0, 0, 1 / 2]), :W => :([1 / 4, 1 / 4, 1 / 4]), :Y₀ => :([Z, -Z, Z]), :U₂ => :([-Z, Z, 1 - Z]), :Λ₀ => :([H, H, -H]), :G₂ => :([-H, 1 - H, H])],
   :hP1 => [:Γ => :([0, 0, 0]), :A => :([0, 0, 1 / 2]), :K => :([1 / 3, 1 / 3, 0]), :H => :([1 / 3, 1 / 3, 1 / 2]), :H₂ => :([1 / 3, 1 / 3, -1 / 2]), :M => :([1 / 2, 0, 0]), :L => :([1 / 2, 0, 1 / 2])],
   :tP1 => [:Γ => :([0, 0, 0]), :Z => :([0, 0, 1 / 2]), :M => :([1 / 2, 1 / 2, 0]), :A => :([1 / 2, 1 / 2, 1 / 2]), :R => :([0, 1 / 2, 1 / 2]), :X => :([0, 1 / 2, 0])],
   :cI1 => [:Γ => :([0, 0, 0]), :H => :([1 / 2, -1 / 2, 1 / 2]), :P => :([1 / 4, 1 / 4, 1 / 4]), :N => :([0, 0, 1 / 2])],
   :mC3 => [:Γ => :([0, 0, 0]), :Y => :([1 / 2, 1 / 2, 0]), :A => :([0, 0, 1 / 2]), :M₂ => :([-1 / 2, 1 / 2, 1 / 2]), :V₂ => :([0, 1 / 2, 0]), :L₂ => :([0, 1 / 2, 1 / 2]), :I => :([-1 + R, R, 1 / 2]), :I₂ => :([1 - R, 1 - R, 1 / 2])],
   :hP2 => [:Γ => :([0, 0, 0]), :A => :([0, 0, 1 / 2]), :K => :([1 / 3, 1 / 3, 0]), :H => :([1 / 3, 1 / 3, 1 / 2]), :M => :([1 / 2, 0, 0]), :L => :([1 / 2, 0, 1 / 2])],
   :hR2 => [:Γ => :([0, 0, 0]), :T => :([1 / 2, -1 / 2, 1 / 2]), :P₀ => :([H, -1 + H, H]), :P₂ => :([H, H, H]), :L => :([1 / 2, 0, 0]), :F => :([1 / 2, -1 / 2, 0])],
   :cF1 => [:Γ => :([0, 0, 0]), :X => :([1 / 2, 0, 1 / 2]), :L => :([1 / 2, 1 / 2, 1 / 2]), :W => :([1 / 2, 1 / 4, 3 / 4]), :W₂ => :([3 / 4, 1 / 4, 1 / 2]), :K => :([3 / 8, 3 / 8, 3 / 4]), :U => :([5 / 8, 1 / 4, 5 / 8])],
   :oA2 => [:Γ => :([0, 0, 0]), :Y => :([1 / 2, 1 / 2, 0]), :T => :([1 / 2, 1 / 2, 1 / 2]), :Z => :([0, 0, 1 / 2]), :S => :([0, 1 / 2, 0]), :R => :([0, 1 / 2, 1 / 2]), :Δ₀ => :([-X, X, 0]), :F₀ => :([X, 1 - X, 0]), :B₀ => :([-X, X, 1 / 2]), :G₀ => :([X, 1 - X, 1 / 2])],
   :cP2 => [:Γ => :([0, 0, 0]), :R => :([1 / 2, 1 / 2, 1 / 2]), :M => :([1 / 2, 1 / 2, 0]), :X => :([0, 1 / 2, 0])],
   :aP2 => [:Γ => :([0, 0, 0]), :Z => :([0, 0, 1 / 2]), :Y => :([0, 1 / 2, 0]), :X => :([1 / 2, 0, 0]), :V => :([1 / 2, 1 / 2, 0]), :U => :([1 / 2, 0, 1 / 2]), :T => :([0, 1 / 2, 1 / 2]), :R => :([1 / 2, 1 / 2, 1 / 2])],
   :oF1 => [:Γ => :([0, 0, 0]), :T => :([1, 1 / 2, 1 / 2]), :Z => :([1 / 2, 1 / 2, 0]), :Y => :([1 / 2, 0, 1 / 2]), :Σ₀ => :([0, H, H]), :U₀ => :([1, 1 - H, 1 - H]), :A₀ => :([1 / 2, 1 / 2 + J, J]), :C₀ => :([1 / 2, 1 / 2 - J, 1 - J]), :L => :([1 / 2, 1 / 2, 1 / 2])],
   :oC1 => [:Γ => :([0, 0, 0]), :Y => :([-1 / 2, 1 / 2, 0]), :T => :([-1 / 2, 1 / 2, 1 / 2]), :Z => :([0, 0, 1 / 2]), :S => :([0, 1 / 2, 0]), :R => :([0, 1 / 2, 1 / 2]), :Σ₀ => :([X, X, 0]), :C₀ => :([-X, 1 - X, 0]), :A₀ => :([X, X, 1 / 2]), :E₀ => :([-X, 1 - X, 1 / 2])],
   :aP3 => [:Γ => :([0, 0, 0]), :Z => :([0, 0, 1 / 2]), :Y => :([0, 1 / 2, 0]), :X => :([1 / 2, 0, 0]), :V₂ => :([1 / 2, -1 / 2, 0]), :U₂ => :([-1 / 2, 0, 1 / 2]), :T₂ => :([0, -1 / 2, 1 / 2]), :R₂ => :([-1 / 2, -1 / 2, 1 / 2])],
   :mC2 => [:Γ => :([0, 0, 0]), :Y => :([1 / 2, 1 / 2, 0]), :A => :([0, 0, 1 / 2]), :M => :([1 / 2, 1 / 2, 1 / 2]), :V₂ => :([0, 1 / 2, 0]), :L₂ => :([0, 1 / 2, 1 / 2])],
   :oI1 => [:Γ => :([0, 0, 0]), :X => :([1 / 2, 1 / 2, -1 / 2]), :S => :([1 / 2, 0, 0]), :R => :([0, 1 / 2, 0]), :T => :([0, 0, 1 / 2]), :W => :([1 / 4, 1 / 4, 1 / 4]), :Σ₀ => :([-Z, Z, Z]), :F₂ => :([Z, 1 - Z, -Z]), :Y₀ => :([H, -H, H]), :U₀ => :([1 - H, H, -H])],
   :mP1 => [:Γ => :([0, 0, 0]), :Z => :([0, 1 / 2, 0]), :B => :([0, 0, 1 / 2]), :Y₂ => :([-1 / 2, 0, 0]), :C₂ => :([-1 / 2, 1 / 2, 0]), :D => :([0, 1 / 2, 1 / 2]), :A => :([-1 / 2, 0, 1 / 2]), :E => :([-1 / 2, 1 / 2, 1 / 2])],
   :cF2 => [:Γ => :([0, 0, 0]), :X => :([1 / 2, 0, 1 / 2]), :L => :([1 / 2, 1 / 2, 1 / 2]), :W => :([1 / 2, 1 / 4, 3 / 4]), :K => :([3 / 8, 3 / 8, 3 / 4]), :U => :([5 / 8, 1 / 4, 5 / 8])],
   :oF2 => [:Γ => :([0, 0, 0]), :T => :([0, 1 / 2, 1 / 2]), :Z => :([1 / 2, 1 / 2, 1]), :Y => :([1 / 2, 0, 1 / 2]), :Λ₀ => :([K, K, 0]), :Q₀ => :([1 - K, 1 - K, 1]), :G₀ => :([1 / 2 - J, 1 - J, 1 / 2]), :H₀ => :([1 / 2 + J, J, 1 / 2]), :L => :([1 / 2, 1 / 2, 1 / 2])],
   :oF3 => [:Γ => :([0, 0, 0]), :T => :([0, 1 / 2, 1 / 2]), :Z => :([1 / 2, 1 / 2, 0]), :Y => :([1 / 2, 0, 1 / 2]), :A₀ => :([1 / 2, 1 / 2 + H, H]), :C₀ => :([1 / 2, 1 / 2 - H, 1 - H]), :B₀ => :([1 / 2 + K, 1 / 2, K]), :D₀ => :([1 / 2 - K, 1 / 2, 1 - K]), :G₀ => :([P, 1 / 2 + P, 1 / 2]), :H₀ => :([1 - P, 1 / 2 - P, 1 / 2]), :L => :([1 / 2, 1 / 2, 1 / 2])],
)

const paramsd_3d = Dict(
   :hR1 => [:D => :((((a * a) / 4) / c) / c), :Y => :(5 / 6 - 2D), :N => :(1 / 3 + D)],
   :tI2 => [:H => :((1 + ((a * a) / c) / c) / 4), :Z => :((((a * a) / 2) / c) / c)],
   :tI1 => [:H => :((1 + ((c * c) / a) / a) / 4)],
   :oC2 => [:X => :((1 + ((b * b) / a) / a) / 4)],
   :oI3 => [:Z => :((1 + ((c * c) / b) / b) / 4), :Y => :((1 + ((a * a) / b) / b) / 4), :D => :((((a * a - c * c) / 4) / b) / b), :M => :((((c * c + a * a) / 4) / b) / b)],
   :oA1 => [:X => :((1 + ((b * b) / c) / c) / 4)],
   :mC1 => [:Z => :((((2 + (a / c) * cosβ) / 4) / sinβ) / sinβ), :H => :(1 / 2 - (2 * Z * c * cosβ) / a), :S => :(3 / 4 - (((((b * b) / 4) / a) / a) / sinβ) / sinβ), :P => :(S - ((3 / 4 - S) * a * cosβ) / c)],
   :oI2 => [:Z => :((1 + ((b * b) / a) / a) / 4), :H => :((1 + ((c * c) / a) / a) / 4), :D => :((((c * c - b * b) / 4) / a) / a), :N => :((((b * b + c * c) / 4) / a) / a)],
   :mC3 => [:Z => :((((a * a) / b) / b + ((1 + (a / c) * cosβ) / sinβ) / sinβ) / 4), :R => :(1 - ((Z * b * b) / a) / a), :E => :(1 / 2 - (2 * Z * c * cosβ) / a), :F => :(E / 2 + (((a * a) / 4) / b) / b + (((a * c * cosβ) / 2) / b) / b), :U => :(2F - Z), :W => :((((c / 2) / a) / cosβ) * ((1 - 4U) + ((a * a * sinβ * sinβ) / b) / b)), :D => :((-1 / 4 + W / 2) - (Z * c * cosβ) / a)],
   :hR2 => [:Z => :(1 / 6 - (((c * c) / 9) / a) / a), :H => :(1 / 2 - 2Z), :N => :(1 / 2 + Z)],
   :oA2 => [:X => :((1 + ((c * c) / b) / b) / 4)],
   :oF1 => [:J => :(((1 + ((a * a) / b) / b) - ((a * a) / c) / c) / 4), :H => :((1 + ((a * a) / b) / b + ((a * a) / c) / c) / 4)],
   :oC1 => [:X => :((1 + ((a * a) / b) / b) / 4)],
   :mC2 => [:Z => :((((a * a) / b) / b + ((1 + (a / c) * cosβ) / sinβ) / sinβ) / 4), :M => :((1 + ((a * a) / b) / b) / 4), :D => :((((-a * c * cosβ) / 2) / b) / b), :X => :(1 / 2 - (2 * Z * c * cosβ) / a), :P => :((1 + Z) - 2M), :S => :(X - 2D)],
   :oI1 => [:Z => :((1 + ((a * a) / c) / c) / 4), :H => :((1 + ((b * b) / c) / c) / 4), :D => :((((b * b - a * a) / 4) / c) / c), :N => :((((a * a + b * b) / 4) / c) / c)],
   :mP1 => [:Y => :((((1 + (a / c) * cosβ) / 2) / sinβ) / sinβ), :N => :(1 / 2 + (Y * c * cosβ) / a)],
   :oF2 => [:J => :(((1 + ((c * c) / a) / a) - ((c * c) / b) / b) / 4), :K => :((1 + ((c * c) / a) / a + ((c * c) / b) / b) / 4)],
   :oF3 => [:H => :(((1 + ((a * a) / b) / b) - ((a * a) / c) / c) / 4), :K => :(((1 + ((b * b) / a) / a) - ((b * b) / c) / c) / 4), :P => :(((1 + ((c * c) / b) / b) - ((c * c) / a) / a) / 4)],
)

# ------------------ COPY OF GIOVANNIPIZZI/SEEKPATH LICENSE ------------------
# The MIT License (MIT)
#
# Copyright (c), 2016-2018, Giovanni Pizzi, ECOLE POLYTECHNIQUE FEDERALE DE
# LAUSANNE (Theory and Simulation of Materials (THEOS) and National Centre for
# Computational Design and Discovery of Novel Materials (NCCR MARVEL)),
# Switzerland. All rights reserved.
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
#
# ---
#
# For a list of open-source software included in this repository, see the file
# open_source_linceses.txt.
# ----------------------------------------------------------------------------


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_script/substitute.jl
# Utilities that are used to substitute expressions in the computational tape.
#############################################################################

"""
    Template_Node

A structure which holds a symbol indicating whether the node is an operator,
a number, or an expression `type`, a value which identifies the function or
symbol `value`, potentially a numeric value `num_value`, and a check that can
be run to verify the node is correct `check`.
"""
struct Template_Node <: Any
    type::Symbol # op, num, expr
    value::Symbol
    num_value::Float64
    check::Function
end

"""
    Template_Graph

Holds a list of Template_Nodes, set of directed edges, lengths, an adjacency
matrix and the number of children.
"""
struct Template_Graph <: Any
    nd::Vector{Template_Node}
    dag::Vector{Pair{Int,Int}}
    ndlen::Int
    daglen::Int
    adj::SparseMatrixCSC{Bool,Int}
    num_children::Vector{Int}
end
function Template_Graph(nd::Dict{Int,Template_Node}, dag::Vector{Pair{Int,Int}})
    ndlen = length(nd)
    daglen = length(dag)
    ndlist = [nd[i] for i in 1:ndlen]
    adj = spzeros(Bool, ndlen, ndlen)
    for i in 1:daglen
        @inbounds p = dag[i]
        @inbounds x = p[1]
        @inbounds y = p[2]
        @inbounds adj[x, y] = true
    end
    num_children = Int[length(nzrange(adj, i)) for i in 1:ndlen]
    Template_Graph(ndlist, dag, ndlen, daglen, adj, num_children)
end

always_true(x) = true
function Template_Node(type::Symbol, value::Symbol; check::Function = always_true)
    @assert (type == :op || type == :expr || type == :num)
    Template_Node(type, value, NaN, check)
end
function Template_Node(type::Symbol, num_value::Float64; check::Function = always_true)
    @assert (type == :num)
    Template_Node(type, :null, num_value, check)
end

const DAG_PATTERNS = Template_Graph[]
const DAG_SUBSTITUTIONS = Template_Graph[]
const DAG_SPDICT = Dict{Int,Int}[]
const DAG_LENGTHS = Int[0,0]

#=
conventions for substition, the expression to be checked always appears at key 1
in the Template_Graph and operations are ordered from low value to high value left to right
so if 1 is a -, and 4 => 1, 3 => 1 then the expression is 4 - 3
=#
"""
    register_substitution!

Specifies that the `src::Template_Graph` should be subsituted out for the
`trg::Template_Graph`.

Conventions for substition, the expression to be checked always appears at key 1
in the Template_Graph and operations are ordered from low value to high value left to right
so if 1 is a -, and 4 => 1, 3 => 1 then the expression is 4 - 3
"""
function register_substitution!(src::Template_Graph, trg::Template_Graph)
    push!(DAG_PATTERNS, src)
    push!(DAG_SUBSTITUTIONS, trg)
    DAG_LENGTHS[1] += 1
    DAG_LENGTHS[2] += 1
    d = Dict{Int,Int}()
    for (i, src_nd) in enumerate(src.nd)
        src_nd_type = src_nd.type
        if (src_nd_type === :expr || src_nd_type === :num)
            for (j, trg_nd) in enumerate(trg.nd)
                if trg_nd.type === src_nd_type
                    if src_nd.value === trg_nd.value
                        @inbounds d[j] = i
                        break
                    end
                end
            end
        end
    end
    push!(DAG_SPDICT, d)
end

function matching_info(x::Template_Node, y::NodeData,
                       const_values::Vector{Float64}, parameter_values::Vector{Float64})
    match_flag = false
    if (x.type == :op)
        if (y.nodetype == CALL)
            if haskey(operator_to_id, x.value)
                @inbounds indx_id = operator_to_id[x.value]
                if indx_id == y.index
                    match_flag = true
                end
            else
                match_flag = false
            end
        elseif (y.nodetype == CALLUNIVAR)
            if haskey(univariate_operator_to_id, x.value)
                @inbounds indx_id = univariate_operator_to_id[x.value]
                if indx_id == y.index
                    match_flag = true
                end
            else
                match_flag = false
            end
        end
    elseif (x.type == :num)
        if (y.nodetype == VALUE)
            @inbounds num_value = const_values[y.index]
        elseif (y.nodetype == PARAMETER)
            @inbounds num_value = parameter_values[y.index]
        else
            match_flag = false
        end
        if (y.nodetype == VALUE) || (y.nodetype == PARAMETER)
            if x.check(num_value)
                match_flag = true
            else
                match_flag = false
            end
        end
    elseif (x.type == :expr)
        match_flag = true
    end
    return match_flag
end

function is_match(pattern::Template_Graph, indx::Int, nd::Vector{NodeData}, dag_adj::SparseMatrixCSC{Bool,Int},
                  const_values::Vector{Float64}, parameter_values::Vector{Float64})
    match_flag = true
    match_dict = Dict{Int,Int}()
    pattern_adj = pattern.adj
    pat_children_arr = rowvals(pattern_adj)
    dag_children_arr = rowvals(dag_adj)

    # do a breadth first search with paired template, nd data,
    # if any pair of children fail then
    pindx_initial = 1
    queue = Tuple{Int,Int}[(pindx_initial, indx)]
    while (~isempty(queue) && (match_flag == true))
        (num_pat, num_dag) = popfirst!(queue)
        @inbounds patt_nd = pattern.nd[num_pat]
        @inbounds dag_nd = nd[num_dag]
        match_info_flag = matching_info(patt_nd, dag_nd, const_values, parameter_values)
        if match_info_flag
            if patt_nd.type === :expr
                @inbounds match_dict[num_pat] = num_dag
            elseif (patt_nd.type === :num) && (patt_nd.num_value === NaN)
                @inbounds match_dict[num_pat] = num_dag
            end
            @inbounds pat_children_idx = nzrange(pattern_adj, num_pat)
            pat_length = length(pat_children_idx)
            if pat_length > 0
                @inbounds dag_children_idx = nzrange(dag_adj, num_dag)
                if pat_length == length(dag_children_idx)
                    for i in 1:pat_length
                        @inbounds dchild = dag_children_idx[i]
                        @inbounds pchild = pat_children_idx[i]
                        @inbounds didx = dag_children_arr[dchild]
                        @inbounds pidx = pat_children_arr[pchild]
                        push!(queue, (pidx, didx))
                    end
                else
                    match_flag = false
                    break
                end

            end
        else
            match_flag = false
            break
        end
    end
    return match_flag, match_dict
end

function find_match(indx::Int, nd::Vector{NodeData}, adj::SparseMatrixCSC{Bool,Int},
                    const_values::Vector{Float64}, parameter_values::Vector{Float64})

    #println("started find match:")
    flag = false
    pattern_number = -1
    match_dict = Dict{Int,Int}()
    @inbounds sub_len = DAG_LENGTHS[1]
    for i in 1:sub_len
        #println("checking pattern i = $i")
        @inbounds pattern = DAG_PATTERNS[i]
        inner_flag, match_dict = is_match(pattern, indx, nd, adj,
                                          const_values, parameter_values)
        #println("inner_flag = $(inner_flag)")
        if inner_flag
            flag = true
            pattern_number = i
            break
        end
    end
    #println("flag = $flag")
    return flag, pattern_number, match_dict
end

#=
Takes a template node and makes the appropriate JuMP node, takes the parent index,
number of child for a pattern element, constant storage vector and it's length
=#
function op_node_to_dag!(x::Template_Node, parent::Int, child_len::Int)
    if child_len > 1
        @inbounds op = operator_to_id[x.value]
        node = NodeData(CALL, op, parent)
    else
        @inbounds op = univariate_operator_to_id[x.value]
        node = NodeData(CALLUNIVAR, op, parent)
    end
    return node
end

function bfs_expr_add!(new_nds::Vector{NodeData}, node_count::Int, num_prt::Int,
                       parent_dict::Dict{Int,Int}, match_dict::Dict{Int,Int},
                       expr_loc::Int, nd::Vector{NodeData}, adj::SparseMatrixCSC{Bool,Int},
                       children_arr::Vector{Int})
    queue = Tuple{Int,Int}[(expr_loc, num_prt)]
    inner_node_count = node_count
    while ~isempty(queue)
        (node_num, prior_prt) = popfirst!(queue) # pop node
        @inbounds active_node = nd[node_num] # store node
        new_node = NodeData(active_node.nodetype, active_node.index, prior_prt)
        inner_node_count += 1 # update node count
        @inbounds parent_dict[num_prt] = inner_node_count
        push!(new_nds, new_node)
        if (active_node.nodetype !== SUBEXPRESSION &&
            active_node.nodetype !== MOIVARIABLE &&
            active_node.nodetype !== VARIABLE &&
            active_node.nodetype !== VALUE)
            @inbounds children_idx = nzrange(adj, node_num)
            if (length(children_idx) > 0)
                for child in children_idx
                    @inbounds idx = children_arr[child]
                    push!(queue, (idx, inner_node_count))
                end
            end
        end
    end
    inner_node_count
end

# we assume a tree structure, so if we don't load child nodes,
# then they are effectively deleted
function substitute!(match_num::Int, node_num::Int, prior_prt::Int, nd::Vector{NodeData},
                     const_list::Vector{Float64}, const_len::Int, node_count::Int,
                     parent_dict::Dict{Int,Int}, match_dict::Dict{Int,Int},
                     queue::Vector{Tuple{Int,Int}}, new_nds::Vector{NodeData},
                     adj::SparseMatrixCSC{Bool,Int}, children_arr::Vector{Int})
    @inbounds subs_template = DAG_SUBSTITUTIONS[match_num]
    @inbounds subs_patt_dict = DAG_SPDICT[match_num]
    subs_adj = subs_template.adj
    subs_children_arr = rowvals(subs_adj)
    queue = Tuple{Int,Int,Int}[(prior_prt, 1, -1)] # node_num, prior_prt, dag_num, dag_prt
    inner_node_count = node_count
    while ~isempty(queue)
        (num_prt, num_sub, _) = popfirst!(queue)
        @inbounds active_node = subs_template.nd[num_sub]
        active_type = active_node.type
        if active_type === :op
            @inbounds active_node_children = subs_template.num_children[num_sub]
            node = op_node_to_dag!(active_node, num_prt, active_node_children)
            push!(new_nds, node)
            inner_node_count += 1
            @inbounds parent_dict[num_prt] = inner_node_count
            @inbounds subs_children_idx = nzrange(subs_adj, num_sub)
            if ~isempty(subs_children_idx)
                for child in subs_children_idx
                    @inbounds sub_child = subs_children_arr[child]
                    push!(queue, (inner_node_count, sub_child, num_sub))
                end
            end
        elseif active_type === :num
            if (active_node.num_value !== NaN)
                push!(const_list, active_node.num_value)
                const_len += 1
                inner_node_count += 1
                @inbounds parent_dict[num_prt] = inner_node_count
                node = NodeData(VALUE, const_len, num_prt)
                push!(new_nds, node)
            else
                @inbounds pindx = subs_patt_dict[num_sub]
                @inbounds expr_loc = match_dict[pindx]
                inner_node_count += 1
                @inbounds parent_dict[num_prt] = inner_node_count
                @inbounds node = NodeData(VALUE, nd[expr_loc].index, num_prt)
                push!(new_nds, node)
            end
        elseif active_type === :expr # Need to only use
            @inbounds pindx = subs_patt_dict[num_sub]
            @inbounds expr_loc = match_dict[pindx]
            temp_node_count = bfs_expr_add!(new_nds, inner_node_count, num_prt,
                                            parent_dict, match_dict, expr_loc,
                                            nd, adj, children_arr)
            inner_node_count = temp_node_count
        end
    end
    return inner_node_count, const_len
end

# searchs through expression breadth first search that short-cirucits
"""
    flatten_expression!

Flattens (usually) the dag by making all registered substitutions for the
expression `expr::_NonlinearExprData`. Performs a depth-first search through
the expression adding the terminal node to the stack, then checking to determine
if it matches a registered substitution pattern. If it doesn't not then node is
added to the new expression graph representation and it's children are added to
the queue. If an expression (node) is identified as a pattern then it is
substituted and any children expression nodes are then checked for patterns until
the depth first search is exhausted.
"""
function flatten_expression!(expr::_NonlinearExprData, parameter_values::Vector{Float64})
    nd = expr.nd
    adj = adjmat(nd)
    children_arr = rowvals(adj)
    node_count = 0
    const_list = expr.const_values
    const_len = length(const_list)
    parent_dict = Dict{Int,Int}(-1 => -1)
    queue = Tuple{Int,Int}[(1,-1)]
    new_nds = NodeData[]
    while ~isempty(queue)
        (node_num, prior_prt) = popfirst!(queue)
        @inbounds active_node = nd[node_num]
        if (active_node.nodetype !== SUBEXPRESSION &&
            active_node.nodetype !== MOIVARIABLE &&
            active_node.nodetype !== VARIABLE &&
            active_node.nodetype !== VALUE)
            is_match, match_num, match_dict = find_match(node_num, nd, adj, const_list, parameter_values)
            if ~is_match
                @inbounds parent_num = parent_dict[prior_prt]
                push!(new_nds, NodeData(active_node.nodetype, active_node.index, parent_num))
                node_count += 1
                @inbounds parent_dict[node_num] = node_count
                @inbounds children_idx = nzrange(adj, node_num) # ADD CHILDREN
                if (length(children_idx) > 0)
                    for child in children_idx
                        @inbounds idx = children_arr[child]
                        push!(queue, (idx, node_num))
                    end
                end
             else
                node_count, const_len = substitute!(match_num, node_num, prior_prt,
                                                    nd, const_list, const_len,
                                                    node_count, parent_dict, match_dict,
                                                    queue, new_nds, adj, children_arr)
            end
        else
            new_parent = parent_dict[prior_prt]
            push!(new_nds, NodeData(active_node.nodetype, active_node.index, new_parent))
            #push!(new_nds, NodeData(active_node.nodetype, active_node.index, prior_prt))
            node_count += 1
            @inbounds parent_dict[node_num] = node_count
        end
    end
    #println("parent_dict: $(parent_dict)")
    expr.nd = new_nds
    return
end

"""
    dag_flattening!

Flattens (usually) the dag by making all registered substitutions for every
nonlinear term in the Optimizer.
"""
function dag_flattening!(x::T) where T <: AbstractOptimizer

    if isa(x._nlp_data.evaluator, NLPEvaluator)
        nlp_data = x._nlp_data.evaluator.m.nlp_data
        params = nlp_data.nlparamvalues
        if ~isnothing(nlp_data.nlobj)
            flatten_expression!(nlp_data.nlobj, params)
        end

        for i = 1:length(nlp_data.nlconstr)
            flatten_expression!(nlp_data.nlconstr[i].terms, params)
        end

        for i = 1:length(nlp_data.nlexpr)
            flatten_expression!(nlp_data.nlexpr[i], params)
        end
    end

    return nothing
end

ext: jl

max_stars_repo_name: MatthewStuber/EAGO

max_issues_repo_name: MatthewStuber/EAGO

max_forks_repo_name: MatthewStuber/EAGO


text:
#Semantic Segmentation on CityScape dataset
#with the convolutional neural network know as PSPNET.
#This script also combines various
# packages from the Julia ecosystem  with Flux.

using Flux
using NNlib
using Parameters: @with_kw
using Flux.Data: DataLoader
using Flux.Optimise: Optimiser, WeightDecay
using CuArrays
import BSON: @save, @load
using CUDAapi
using ProgressMeter
using FileIO
using Images
import DrWatson: savename, struct2dict
using Statistics, Random

#We create our model which has three parts
@with_kw mutable struct PSPNet
    #first part is making the residual blocks,In total there are three residual blocks where each block contains three convolutional layers in series and one skip connection
    number_of_classes = 3
    Filters = (32,32,64)
    #Start of first residual block
    conv1 = Conv((1,1),number_of_classes=>Filters[1], pad = 0, dilation = 1)
    batch1 = BatchNorm(Filters[1])
    conv2 = Conv((3,3),Filters[1]=>Filters[2], pad = 1, dilation = 1)
    batch2 = BatchNorm(Filters[2])
    conv3 = Conv((1,1),Filters[2]=>Filters[3], pad = 0, dilation = 1)
    batch3 = BatchNorm(Filters[3])

    conv_skip = Conv((3,3), number_of_classes=>Filters[3];
                    stride = 1, pad = 1, dilation = 1)
    batch_skip = BatchNorm(Filters[3])

    #End of first residual block

    #Start of second residual block
    Filters_2 = (64,64,128)
    conv4 = Conv((1,1), 64=>Filters_2[1], pad = 0, dilation = 1)
    batch4 = BatchNorm(Filters_2[1])
    conv5 = Conv((3,3),Filters_2[1]=>Filters_2[2], pad = 1, dilation = 1)
    batch5 = BatchNorm(Filters_2[2])
    conv6 = Conv((1,1),Filters_2[2] => Filters_2[3], pad = 0, dilation = 1)
    batch6 = BatchNorm(Filters_2[3])

    conv_skip_2 = Conv((3,3), 64=>Filters_2[3];
                    stride = 1, pad = 1, dilation = 1)
    batch_skip_2 = BatchNorm(Filters_2[3])

    #End of second residual block

    #Start of third residual block
    Filters_3 = (128,128,256)
    conv7 = Conv((1,1), 128=>Filters_3[1], pad = 0, dilation = 1)
    batch7 = BatchNorm(Filters_3[1])
    conv8 = Conv((3,3),Filters_3[1]=>Filters_3[2], pad = 1, dilation = 1)
    batch8 = BatchNorm(Filters_3[2])
    conv9 = Conv((1,1),Filters_3[2] => Filters_3[3], pad = 0, dilation = 1)
    batch9 = BatchNorm(Filters_3[3])

    conv_skip_3 = Conv((3,3), 128=>Filters_3[3];
        stride = 1, pad = 1, dilation = 1)
    batch_skip_3 = BatchNorm(Filters_3[3])

    #End of third residual block

    #Now we start making the pyramid feature map
    mp1 = MeanPool((1,1); pad = 0, stride = 1)
    ct1 = Conv((1,1), 256 =>64 ,stride = 1, pad = 0, dilation = 1)

    mp2 = MeanPool((2,2); pad = 0, stride = 2)
    ct2 = Conv((1,1),  256=>64, stride = 1, pad = 0, dilation = 1)

    mp3 = MeanPool((4,4); pad = 0, stride = 4)
    ct3 = Conv((1,1),  256=>64 ,stride = 1, pad = 0, dilation = 1)

    mp4 = MeanPool((8,8); pad = 0, stride = 8)
    ct4 = Conv((1,1),  256=>64, stride = 1, pad = 0, dilation = 1)

    #End of feature pyramid

    #Now we make the last convolutional layer which has input from the pyramid feature map plus the fetaures extracted from residual blocks
    #so the output from residual blocks is (256,256,256,1) and the output from the pyramid feature map is four feature maps
    #with shape (256,256,64,1).
    #The last convolutional layer takes input which is concatenation of output from residual blocks and Pyramid feature map
    conv_last = Conv((3,3),512=>number_of_classes, pad = 1, dilation = 1)
    batch_last = BatchNorm(number_of_classes)

end

mutable struct Model
   x_output
   x_skip_output
   x_output_2
   x_skip_output_2
   x_output_3
   x_skip_output_3
   pyramid_feature_1
   pyramid_feature_2
   pyramid_feature_3
   pyramid_feature_4
   last_conv
end

function pspnet(m::PSPNet)

    x_output = Chain(m.conv1, m.batch1, x -> leakyrelu.(x,0.01), m.conv2, m.batch2, x->leakyrelu.(x,0.01), m.conv3, m.batch3, x->leakyrelu.(x,0.01))
    x_skip_output = Chain(m.conv_skip,m.batch_skip)
    x_output_2 = Chain(m.conv4, m.batch4, x->leakyrelu.(x,0.01), m.conv5, m.batch5, x->leakyrelu.(x,0.01), m.conv6, m.batch6, x->leakyrelu.(x,0.01))
    x_skip_output_2 = Chain(m.conv_skip_2,m.batch_skip_2)
    x_output_3 = Chain(m.conv7, m.batch7, x->leakyrelu.(x,0.01), m.conv8, m.batch8, x->leakyrelu.(x,0.01), m.conv9, m.batch9, x->leakyrelu.(x,0.01))
    x_skip_output_3 = Chain(m.conv_skip_3, m.batch_skip_3)
    pyramid_feature_1 = Chain(m.mp1, m.ct1, x->upsample(x,(1,1,1,1)))
    pyramid_feature_2 = Chain(m.mp2, m.ct2, x->upsample(x,(2,2,1,1)))
    pyramid_feature_3 = Chain(m.mp3, m.ct3, x->upsample(x,(4,4,1,1)))
    pyramid_feature_4 = Chain(m.mp4, m.ct4, x->upsample(x,(8,8,1,1)))
    last_conv = Chain(m.conv_last, m.batch_last)
    return Model(x_output, x_skip_output, x_output_2, x_skip_output_2, x_output_3, x_skip_output_3,
                pyramid_feature_1, pyramid_feature_2, pyramid_feature_3, pyramid_feature_4, last_conv)

end

@with_kw mutable struct Args
    batchsize = 10
    λ=0
    η = 0.01
    ρ = 0.99
    seed=0
    epochs = 10
    cuda = true
    checktime = 5        # Save the model every `checktime` epochs. Set to 0 for no checkpoints.
    savepath = nothing
    train_img_path = joinpath(homedir(), "Downloads", "CityScape_data","TrainImg")  #you can change the path as per your choice
    train_mask_path = joinpath(homedir(), "Downloads", "CityScape_data","TrainMask")
    test_img_path = joinpath(homedir(), "Downloads", "CityScape_data","TestImg")
    test_mask_path = joinpath(homedir(), "Downloads", "CityScape_data","TestMask")
end

function upsample(x,ratio)
  (h, w, c, n) = size(x)
  y = similar(x, (1, ratio[1], 1, ratio[2], 1, 1))
  fill!(y, 1)
  z = reshape(x, (h, 1, w, 1, c, n))  .* y
  reshape(permutedims(z, (2,1,4,3,5,6)), size(x) .* ratio)
end

function forward(feature_extractor::Model,input::Array)
    sum = feature_extractor.x_output(input)+feature_extractor.x_skip_output(input)
    sum2 = feature_extractor.x_output_2(sum)+feature_extractor.x_skip_output_2(sum)
    sum3 = feature_extractor.x_output_3(sum2)+feature_extractor.x_skip_output_3(sum2)
    pf1 = feature_extractor.pyramid_feature_1(sum3)
    pf2 = feature_extractor.pyramid_feature_2(sum3)
    pf3 = feature_extractor.pyramid_feature_3(sum3)
    pf4 = feature_extractor.pyramid_feature_4(sum3)
    matrices = (sum3,pf1,pf2,pf3,pf4)
    concat_matrices = cat(matrices...;dims=(3))
    output = feature_extractor.last_conv(concat_matrices)
    return output
end

rgb_to_array(img)=begin
    H, W = height(img), width(img)
    a = zeros(H, W, 3)
    for w in 1:W
        for h in 1:H
            rgb = img[h,w]
            a[h,w,1] = rgb.r
            a[h,w,2] = rgb.g
            a[h,w,3] = rgb.b
        end
    end
    a
end

function get_data(Train_img_path::String, Train_Label_path::String, Test_img_path::String, Test_Label_path::String)
    Train_data = readdir(Train_img_path, join = true, sort = true)
    Train_Label = readdir(Train_Label_path, join = true, sort = true)
    Test_data = readdir(Test_img_path, join = true, sort = true)
    Test_Label = readdir(Test_Label_path, join = true, sort = true)
    Xtrain = zeros(256,256,3,length(Train_data))
    Ytrain = zeros(256,256,3,length(Train_Label))
    Xtest = zeros(256,256,3,length(Test_data))
    Ytest = zeros(256,256,3,length(Test_Label))
    i,j = 1,1
    for (img_path,mask_path) in zip(Train_data,Train_Label)

        img = load(img_path)
        mask = load(mask_path)
        img = rgb_to_array(img)
        mask = rgb_to_array(mask)
        Xtrain[:,:,:,i] = img
        Ytrain[:,:,:,i] = mask
        i += 1
    end
    for (test_img_path, test_mask_path) in zip(Test_data,Test_Label)
        test_img = load(test_img_path)
        test_mask = load(test_mask_path)
        img = rgb_to_array(test_img)
        mask = rgb_to_array(test_mask)
        Xtest[:,:,:,j] = img
        Ytest[:,:,:,j] = mask
        j += 1
    end

    train_loader = DataLoader(Xtrain,Ytrain,batchsize = args.batchsize,shuffle = false)
    test_loader = DataLoader(Xtest,Ytest,batchsize = args.batchsize,shuffle = false)

    return train_loader, test_loader
end

function train(psp_layers::PSPNet, Args::Args)
    args = Args()
    args.seed > 0 && Random.seed!(args.seed)
    use_cuda = args.cuda && CUDAapi.has_cuda_gpu()
    if use_cuda
        device = gpu
        @info "Training on GPU"
    else
        device = cpu
        @info "Training on CPU"
    end

    if args.savepath == nothing
     experiment_folder = savename("PSPNet", args, scientific=3,
                 accesses = [:batchsize, :η, :λ]) # construct path from these fields
     args.savepath = joinpath("runs", experiment_folder)
    end

    ## DATA
    train_loader, test_loader = get_data(args.train_img_path , args.train_mask_path , args.test_img_path , args.test_mask_path)
    @info "Dataset CityScape: $(train_loader.nobs) train and $(test_loader.nobs) test examples"

    ## MODEL AND OPTIMIZER
    layers = psp_layers
    model = pspnet(layers)|> device
    ps = Flux.params(model.x_output , model.x_skip_output , model.x_output_2 , model.x_skip_output_2 , model.x_output_3 , model.x_skip_output_3 ,
                    model.pyramid_feature_1 , model.pyramid_feature_2 , model.pyramid_feature_3 , model.pyramid_feature_4 ,
                    model.last_conv)
    loss(x , y) = Flux.mse(x , y)
    opt = ADAM(args.η)
    if args.λ > 0
        opt = Optimiser(opt , WeightDecay(args.λ))
    end

    ## TRAINING
    @info "Start Training"
    report(0)
    for epoch in 1:args.epochs
        p = ProgressMeter.Progress(length(train_loader))

        for (x, y) in train_loader
            x, y = x |> device, y |> device
            gs = Flux.gradient(ps) do
                ŷ = forward(model,x)
                loss(ŷ, y)
            end
            Flux.Optimise.update!(opt, ps, gs)
            ProgressMeter.next!(p)   # comment out for no progress bar
        end

        if args.checktime > 0 && epoch % args.checktime == 0
            !ispath(args.savepath) && mkpath(args.savepath)
            modelpath = joinpath(args.savepath, "model.bson")
            let model=cpu(model), args=struct2dict(args)
                BSON.@save modelpath model epoch args
            end
            @info "Model saved in \"$(modelpath)\""
        end
    end
end

max_stars_repo_name: ocissor/model-zoo

max_issues_repo_name: ocissor/model-zoo

max_forks_repo_name: ocissor/model-zoo


text:
# Copyright (c) Microsoft Corporation.
# Licensed under the MIT license.

module LstmTests

dir = @__DIR__
lstm_test_implementations = Tuple{String, Float64}[
    ("$dir/../../../../src/julia/modules/Zygote/ZygoteLSTM.jl", 1e-8)
]

map!(tup -> (abspath(tup[1]), tup[2]), lstm_test_implementations, lstm_test_implementations)

using Test
include("../../../../src/julia/shared/load.jl")
include("../../../../src/julia/runner/load.jl")
include("load.jl")
import ADPerfTest
import TestLoader
using LSTMData
using TestUtils

@testset "LSTM Module Tests" begin

@testset "LSTM Module Test ($(basename(module_path)))" for (module_path, tolerance) in lstm_test_implementations
    module_dir, module_filename = splitdir(module_path)
    module_name, module_ext = splitext(module_filename)
    # Assert that the path to the module is correct
    ext_r = @test module_ext == ".jl"
    ispath_r = @test (ispath(module_path) && !isdir(module_path))
    if isa(ext_r, Test.Pass) && isa(ispath_r, Test.Pass)
        need_modify_load_path = !(module_dir ∈ LOAD_PATH)
        if need_modify_load_path
            push!(LOAD_PATH, module_dir)
        end
        try
            test = TestLoader.get_lstm_test(module_name)
            # Module loads
            @test isa(test, ADPerfTest.Test{LSTMInput, LSTMOutput})
            input = load_lstm_input("$dir/../../../../data/lstm/lstm_l2_c1024.txt")
            output = empty_lstm_output()
            test.prepare!(test.context, input)
            correct_loss = 6.66665179558852183e-01
            correct_grad = [
                3.14387024996101187e-06, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -3.67840807336070475e-05, -1.54848612683180522e-04, 5.48717267869773910e-06, -4.04527413533285253e-06, -8.77353908686236070e-07, -1.81973718092835457e-06, -6.61404620766666418e-07, 1.25432314145999742e-04, 5.26459549995064566e-06, 3.42833197850334681e-06, 2.79076684010150125e-05, 4.45603055235501376e-05, 4.04185405141412220e-05, 1.13487138666022078e-05, -7.36119414163935681e-05, -3.92318596883106641e-04, 6.03607595779890245e-06, -2.50749024433336479e-05, -4.83955404587018941e-06, -4.94175793354897013e-06, -1.49848451794559008e-06, 4.01789457601843220e-06, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -6.16501771215585166e-05, -2.14471058586458561e-04, 5.56661202545761896e-06, -5.22021372686648383e-06, -1.95401365542668346e-06, -4.68776054707821880e-06, -2.67879854428802963e-06, 8.49896476968742872e-04, 2.14480696808520274e-04, 1.27240956373907973e-04, 1.03380278076147588e-04, 3.26983387274766885e-05, 1.47321495437521366e-04, 2.03767485364521309e-05, -1.71027241883012373e-04, -9.97102906947628905e-04, 1.05573445256271496e-04, -4.10813160791978987e-05, -4.66363145389210376e-06, -5.61458724856520127e-06, -1.72719535440704236e-06, 7.41562400488615767e-04, 1.35905892236421730e-04, 9.18319427682989299e-05, 1.14039861417365726e-04, 1.31037490514354042e-04, 1.89958003915723621e-04, 6.66667089111214168e-05, -2.99753093683376040e-04, -1.39428488543169235e-03, 7.81997486731591678e-05, -7.11003463557882681e-05, -1.56689414959509608e-05, -1.59132673852719420e-05, -5.88701856095873723e-06, 3.65214776221949554e-04, 1.08440867176087962e-04, 6.75251197117822909e-05, 5.07904844905876160e-05, 6.94964603070496247e-05, 9.71513382958139477e-05, 2.40897410354102865e-05, -1.31891358419510490e-04, -7.36772740642002493e-04, 4.48757608883305795e-05, -4.51743229379305454e-05, -7.63715831354149397e-06, -8.36900345795258848e-06, -2.41858043904271068e-06, 8.00930878140701909e-04, 8.05738754618121506e-05, 7.04205639535801850e-05, 1.34615500556859374e-04, 2.07383290113708586e-04, 2.94203848388590759e-04, 2.08296283404346797e-04, -5.00619564307732703e-04, -1.90384053275097134e-03, 6.83889421204583172e-05, -8.85415990620916452e-05, -4.01278522469193820e-05, -3.19293895829058307e-05, -2.48575295654621768e-05, 2.47822472805507175e-03, 5.07502578529807780e-03, 2.53859920727141910e-03, 1.88350319998056696e-04, 5.10157230919509725e-05, 3.54335963584264193e-04, 4.32702650224108475e-05, -3.06537980830570744e-04, -1.87256106103047193e-03, 7.86527299400905083e-04, -7.40280147056494078e-05, -7.36222276800657044e-06, -9.51938152613299431e-06, -2.78936845972344488e-06, 5.38571851736894072e-04, 4.12509633131798620e-05, 6.47788903296268593e-05, 3.76860255482900047e-04, 1.00145332877384172e-03, 2.29657803047558681e-04, 4.44627297218581029e-04, -9.73952035222110210e-04, -4.51211367568013846e-03, 1.41925554887538418e-04, -1.66675086094364628e-05, -2.28915871203058926e-05, -6.90861176975165640e-05, -2.81513790989548166e-04, 2.53531909024735675e-04, 2.91530880231834222e-04, 2.87839568382086426e-04, 6.61417346851010207e-05, 6.11690024226054017e-04, 1.08165251111585849e-04, 9.94833840089197921e-05, -5.39732825546851164e-04, -1.95313692267145384e-03, 2.87262803065020395e-04, -6.42028361615228470e-06, -6.81225636233755266e-06, -3.58438396937832646e-05, -1.41607977576401395e-04, 9.90063334626659942e-04, 6.71457692678222891e-05, 9.64171460882922153e-05, 2.83657185457934097e-04, 1.41431021522050158e-03, 7.42601545855336369e-04, 1.63135992873617838e-04, -1.80905573085369175e-03, -5.78938583485935632e-03, 1.08933179139588087e-04, -2.82053334001221257e-04, -1.58448932394055300e-04, -1.62053197372598410e-04, -4.44062656277905885e-04, 1.23784405389406552e-03, 3.70797204635395672e-04, 8.95872992208362269e-04, 9.54650103551505725e-04, 7.94492683033071450e-04, 1.33045599071493839e-04, 1.90972131022078425e-03, -1.27723370300947075e-03, -7.98761497515965146e-03, 1.31804902759735014e-03, -1.18749502523366873e-05, -1.17265333639628955e-05, -5.72852084344175575e-05, -1.76267731147033947e-04, 1.57111882911252715e-03, 9.63872765804848195e-04, 1.29692708705586861e-03, 6.85018812295237929e-04, 1.56084730076406770e-03, 5.51507203835322228e-04, 9.44664300904442429e-04, -1.74595527620669054e-03, -8.47272040239691569e-03, 1.06146398161773513e-03, -3.01166484114288631e-05, -3.61150963335128990e-05, -1.17040277318999298e-04, -4.53714566010890040e-04, 4.83883598142451641e-04, 5.36532886336024261e-04, 7.10491025781043260e-04, 4.03370948422198858e-04, 9.21362107223384799e-04, 1.83544588672272390e-04, 5.84842906913008364e-04, -9.12226698645542925e-04, -3.71928760349100754e-03, 5.24341732641706525e-04, -8.78104574502159225e-06, -9.40356235117202346e-06, -6.16753351894278713e-05, -2.59729129852274867e-04, 2.88643409216300618e-03, 1.56939967321338843e-03, 1.92928469338899343e-03, 5.15538874537588353e-04, 2.20418224227594915e-03, 1.78306235870001779e-03, 3.46590677977152962e-04, -3.24323971790235277e-03, -1.08700028019174302e-02, 8.13601709951329772e-04, -5.09026696834261072e-04, -2.50224828500900109e-04, -2.74229202243154427e-04, -7.15862244954266883e-04, 2.36221822129175068e-03, 6.82455302403955932e-04, 2.21184971904590886e-03, 5.81458374769466737e-03, 1.19699412544996883e-03, 2.26008787507052335e-04, 1.12207195875873196e-02, -2.15943584477549123e-03, -1.52105422913329149e-02, 2.40692059621100754e-03, -1.62908509623434489e-05, -1.62307954324598644e-05, -9.93149651364052892e-05, -3.23349485361925845e-04, 8.95056656284422921e-06, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.73683965955792591e-04, -3.21979564391077671e-04, 7.71925547683623476e-06, -2.98807224102069401e-05, -1.11738854533001905e-05, -2.92505411923477579e-06, -2.04187281292674588e-06, 1.14042748169088035e-02, 9.93947389065250562e-03, 4.53876797372907911e-03, 1.46753585074672422e-03, 1.59578501610211192e-02, 7.13215996147348839e-03, 2.52120644391976836e-03, -1.79990239825079217e-02, -2.62659614088943671e-02, 2.44293992080160194e-03, -3.67074714937877696e-03, -1.31444481103537205e-02, -4.54999892021307108e-03, -1.09270983790279372e-02, 2.17348662792538061e-02, 1.82857375496154884e-02, 1.11887441899222798e-02, 8.97672919990902263e-03, 2.39853893219980456e-02, 1.20800375241284126e-02, 1.48638788575309582e-02, -3.02324690101714345e-02, -5.00116775406580991e-02, 4.46620942840666828e-03, -5.23602406001136191e-03, -1.79316285235904226e-02, -7.73167106472800580e-03, -2.00225690895980656e-02
            ]
            test.calculate_objective!(test.context, 1)
            test.output!(output, test.context)
            # Objective is calculated correctly
            @test correct_loss ≈ output.objective atol=tolerance
            test.calculate_objective!(test.context, 3)
            test.output!(output, test.context)
            # Objective is calculated correctly for times = 3
            @test correct_loss ≈ output.objective atol=tolerance
            test.calculate_jacobian!(test.context, 1)
            test.output!(output, test.context)
            # Gradient is calculated correctly
            @test correct_grad ≈ output.gradient atol=tolerance
            test.calculate_jacobian!(test.context, 3)
            test.output!(output, test.context)
            # Gradient is calculated correctly for times = 3
            @test correct_grad ≈ output.gradient atol=tolerance
            # Objective runs multiple times
            @test can_objective_run_multiple_times!(test.context, test.calculate_objective!)
            # Jacobian runs multiple times
            @test can_objective_run_multiple_times!(test.context, test.calculate_jacobian!)
        finally
            pop!(LOAD_PATH)
        end
    end
end

end

end

max_stars_repo_name: dsyme/ADBench

max_issues_repo_name: dsyme/ADBench

max_forks_repo_name: dsyme/ADBench


text:
@testset "999.available-captures-for-rook.jl" begin
    board = [
        '.' '.' '.' '.' '.' '.' '.' '.'
        '.' '.' '.' 'p' '.' '.' '.' '.'
        '.' '.' '.' 'R' '.' '.' '.' 'p'
        '.' '.' '.' '.' '.' '.' '.' '.'
        '.' '.' '.' '.' '.' '.' '.' '.'
        '.' '.' '.' 'p' '.' '.' '.' '.'
        '.' '.' '.' '.' '.' '.' '.' '.'
        '.' '.' '.' '.' '.' '.' '.' '.'
    ]
    @test num_rook_captures(board) == 3

    board = [
        '.' '.' '.' '.' '.' '.' '.' '.'
        '.' 'p' 'p' 'p' 'p' 'p' '.' '.'
        '.' 'p' 'p' 'B' 'p' 'p' '.' '.'
        '.' 'p' 'B' 'R' 'B' 'p' '.' '.'
        '.' 'p' 'p' 'B' 'p' 'p' '.' '.'
        '.' 'p' 'p' 'p' 'p' 'p' '.' '.'
        '.' '.' '.' '.' '.' '.' '.' '.'
        '.' '.' '.' '.' '.' '.' '.' '.'
    ]
    @test num_rook_captures(board) == 0

    board = [
        '.' '.' '.' '.' '.' '.' '.' '.'
        '.' '.' '.' 'p' '.' '.' '.' '.'
        '.' '.' '.' 'p' '.' '.' '.' '.'
        'p' 'p' '.' 'R' '.' 'p' 'B' '.'
        '.' '.' '.' '.' '.' '.' '.' '.'
        '.' '.' '.' 'B' '.' '.' '.' '.'
        '.' '.' '.' 'p' '.' '.' '.' '.'
        '.' '.' '.' '.' '.' '.' '.' '.'
    ]
    @test num_rook_captures(board) == 3
end


text:
srand(5678)

p = 5
n = 8
k = 3

for T in (Float64, Float32)
    Wg = max(rand(T, p, k) .- 0.3, 0)
    Hg = max(rand(T, k, n) .- 0.3, 0)
    X = Wg * Hg

    for alg in (:multmse, :multdiv, :projals, :alspgrad)
        for init in (:random, :nndsvd, :nndsvda, :nndsvdar)
            ret = NMF.nnmf(X, k, alg=alg, init=init)
        end
    end
end


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

# mapping formula from a point (xi, eta) in reference space [-1,1]^2 to a point (x,y)
# in physical coordinate space for a quadrilateral element with straight sides
#     Alg. 95 from the blue book of Kopriva
function straight_side_quad_map(xi, eta, corner_points)

  x = (0.25 * (  corner_points[1,1] * (1.0 - xi) * (1.0 - eta)
               + corner_points[2,1] * (1.0 + xi) * (1.0 - eta)
               + corner_points[3,1] * (1.0 + xi) * (1.0 + eta)
               + corner_points[4,1] * (1.0 - xi) * (1.0 + eta)) )

  y = (0.25 * (  corner_points[1,2] * (1.0 - xi) * (1.0 - eta)
               + corner_points[2,2] * (1.0 + xi) * (1.0 - eta)
               + corner_points[3,2] * (1.0 + xi) * (1.0 + eta)
               + corner_points[4,2] * (1.0 - xi) * (1.0 + eta)) )

  return x, y
end

# Compute the metric terms for the straight sided quadrilateral mapping
#     Alg. 100 from the blue book of Kopriva
function straight_side_quad_map_metrics(xi, eta, corner_points)

  X_xi  = ( 0.25 * (  (1.0 - eta) * (corner_points[2,1] - corner_points[1,1])
                    + (1.0 + eta) * (corner_points[3,1] - corner_points[4,1])) )

  X_eta = ( 0.25 * (  (1.0 - xi) * (corner_points[4,1] - corner_points[1,1])
                    + (1.0 + xi) * (corner_points[3,1] - corner_points[2,1])) )

  Y_xi  = ( 0.25 * (  (1.0 - eta) * (corner_points[2,2] - corner_points[1,2])
                    + (1.0 + eta) * (corner_points[3,2] - corner_points[4,2])) )

  Y_eta = ( 0.25 * (  (1.0 - xi) * (corner_points[4,2] - corner_points[1,2])
                    + (1.0 + xi) * (corner_points[3,2] - corner_points[2,2])) )

  return X_xi, X_eta, Y_xi, Y_eta
end

# construct the (x,y) node coordinates in the volume of a straight sided element
function calc_node_coordinates!(node_coordinates::AbstractArray{<:Any, 4}, element, nodes, corners)

  for j in eachindex(nodes), i in eachindex(nodes)
    node_coordinates[:, i ,j ,element] .= straight_side_quad_map(nodes[i], nodes[j], corners)
  end

  return node_coordinates
end

# construct the metric terms for a straight sided element
function calc_metric_terms!(jacobian_matrix, element, nodes, corners)

  # storage format:
  #   jacobian_matrix[1,1,:,:,:] <- X_xi
  #   jacobian_matrix[1,2,:,:,:] <- X_eta
  #   jacobian_matrix[2,1,:,:,:] <- Y_xi
  #   jacobian_matrix[2,2,:,:,:] <- Y_eta
  for j in eachindex(nodes), i in eachindex(nodes)
    (jacobian_matrix[1, 1, i, j, element],
     jacobian_matrix[1, 2, i, j, element],
     jacobian_matrix[2, 1, i, j, element],
     jacobian_matrix[2, 2, i, j, element]) = straight_side_quad_map_metrics(nodes[i], nodes[j],
                                                                            corners)
  end

  return jacobian_matrix
end

# construct the normal direction vectors (but not actually normalized) for a straight sided element
# normalization occurs on the fly during the surface flux computation
function calc_normal_directions!(normal_directions, element, nodes, corners)

  # normal directions on the boundary for the left (local side 4) and right (local side 2)
  for j in eachindex(nodes)
    # side 2
    X_xi, X_eta, Y_xi, Y_eta = straight_side_quad_map_metrics(1.0, nodes[j], corners)
    Jtemp = X_xi * Y_eta - X_eta * Y_xi
    normal_directions[1, j, 2, element] = sign(Jtemp) * ( Y_eta )
    normal_directions[2, j, 2, element] = sign(Jtemp) * (-X_eta )

    # side 4
    X_xi, X_eta, Y_xi, Y_eta = straight_side_quad_map_metrics(-1.0, nodes[j], corners)
    Jtemp =  X_xi * Y_eta - X_eta * Y_xi
    normal_directions[1, j, 4, element] = -sign(Jtemp) * ( Y_eta )
    normal_directions[2, j, 4, element] = -sign(Jtemp) * (-X_eta )
  end

  # normal directions on the boundary for the top (local side 3) and bottom (local side 1)
  for i in eachindex(nodes)
    # side 1
    X_xi, X_eta, Y_xi, Y_eta = straight_side_quad_map_metrics(nodes[i], -1.0, corners)
    Jtemp =  X_xi * Y_eta - X_eta * Y_xi
    normal_directions[1, i, 1, element] = -sign(Jtemp) * (-Y_xi )
    normal_directions[2, i, 1, element] = -sign(Jtemp) * ( X_xi )

    # side 3
    X_xi, X_eta, Y_xi, Y_eta = straight_side_quad_map_metrics(nodes[i], 1.0, corners)
    Jtemp = X_xi * Y_eta - X_eta * Y_xi
    normal_directions[1, i, 3, element] = sign(Jtemp) * (-Y_xi )
    normal_directions[2, i, 3, element] = sign(Jtemp) * ( X_xi )
  end

  return normal_directions
end

end # @muladd


text:
using SpecialFunctions
function  poisson(theta, y, x)
    lambda = exp.(x*theta)
    n = size(y,1)
    logdensity = y .* (x*theta) .- lambda .- loggamma.(y .+ 1.)
end

max_stars_repo_name: Hiroakiyusheng/Econometrics

max_issues_repo_name: Hiroakiyusheng/Econometrics

max_forks_repo_name: Hiroakiyusheng/Econometrics


text:
struct BoolBridge{T, B1<:MOIBC.SetMapBridge{T}, B2<:MOIBC.SetMapBridge{T}, F1, F2, S1, S2} <: MOIBC.SetMapBridge{T, S2, S1, F2, F1}
    con_idx::CI
end

function MOI.supports_constraint(
    bridge::Type{<:BoolBridge{T, B1, B2}},
    fct::Type{<:MOI.AbstractFunction},
    set::Type{<:CS.BoolSet{F1,F2,F1dim,F2dim,S1,S2}}
) where {T, B1, B2, F1, F2, F1dim, F2dim, S1<:BoolSet, S2<:BoolSet}
    return MOI.supports_constraint(bridge, F1, S1) && MOI.supports_constraint(bridge, F2, S2)
end

function MOI.supports_constraint(
    bridge::Type{<:BoolBridge{T, B1, B2}},
    fct::Type{<:MOI.AbstractFunction},
    set::Type{<:CS.BoolSet{F1,F2,F1dim,F2dim,S1,S2}}
) where {T, B1, B2, F1, F2, F1dim, F2dim, S1<:BoolSet, S2}
    direct_support, inner_bridge = how_supported(B1, B2, F2, S2)
    if direct_support || inner_bridge !== nothing
        return MOI.supports_constraint(bridge, F1, S1)
    end
end

function MOI.supports_constraint(
    bridge::Type{<:BoolBridge{T, B1, B2}},
    fct::Type{<:MOI.AbstractFunction},
    set::Type{<:CS.BoolSet{F1,F2,F1dim,F2dim,S1,S2}}
) where {T, B1, B2, F1, F2, F1dim, F2dim, S1, S2<:BoolSet}
    direct_support, inner_bridge = how_supported(B1, B2, F1, S1)
    if direct_support || inner_bridge !== nothing
        return MOI.supports_constraint(bridge, F2, S2)
    end
end

function MOI.supports_constraint(
    bridge::Type{<:BoolBridge{T, B1, B2}},
    fct::Type{<:MOI.AbstractFunction},
    set::Type{<:CS.BoolSet{F1,F2,F1dim,F2dim,S1,S2}}
) where {T, B1, B2, F1, F2, F1dim, F2dim, S1, S2}
    direct_support, inner_bridge = how_supported(B1, B2, F1, S1)
    if !direct_support && inner_bridge === nothing
        return false
    end
    direct_support, inner_bridge = how_supported(B1, B2, F2, S2)
    if !direct_support && inner_bridge === nothing
        return false
    end
    return true
end

function added_constraint_types(
    bridge::Type{<:BoolBridge{T, B1, B2}},
    ::Type{<:CS.BoolSet{F1,F2,F1dim,F2dim,S1,S2}}
) where {T,B1,B2,F1, F2, F1dim, F2dim, S1, S2}
    set = get_bridged_and_set(bridge,F1,F2,F1dim,F2dim,S1,S2)
end

function MOIBC.concrete_bridge_type(
    ::Type{<:BoolBridge{T,B1,B2}},
    G::Type{<:MOI.AbstractFunction},
    ::Type{<:CS.BoolSet{F1,F2,F1dim,F2dim,S1,S2}}
) where {T,B1,B2,F1, F2, F1dim, F2dim, S1, S2}
    return BoolBridge{T,B1,B2,F1,F2,S1,S2}
end

function MOIBC.bridge_constraint(B::Type{<:BoolBridge{T, B1, B2, F1, F2, S1, S2}}, model, fct, set::CS.BoolSet) where {T, B1, B2, F1, F2, S1, S2}
    new_fct = map_function(B, fct, set)
    new_set = map_set(B, set)
    return BoolBridge{T,B1,B2,F1,F2,S1,S2}(MOI.add_constraint(model, new_fct, new_set))
end

function unpack_constraint_types(ct1, ct2)
    new_F1 = ct1[1][1]
    new_F2 = ct2[1][1]
    new_S1 = ct1[1][2]
    new_S2 = ct2[1][2]
    return new_F1, new_F2, new_S1, new_S2
end

function how_supported(
    B1::Type{<:MOIBC.SetMapBridge{T}},
    B2::Type{<:MOIBC.SetMapBridge{T}},
    F, S
) where {T}
    supported_by_B1 = MOI.supports_constraint(B1, F, S)
    supported_by_B2 = MOI.supports_constraint(B2, F, S)
    if !supported_by_B1 && !supported_by_B2
        return true, nothing
    elseif supported_by_B1
        return false, B1
    elseif supported_by_B2
        return false, B2
    else
        return false, nothing
    end
end

function get_constraint_types(
    B1::Type{<:MOIBC.SetMapBridge{T}},
    B2::Type{<:MOIBC.SetMapBridge{T}},
    F, S
) where {T}
    direct_support, inner_bridge = how_supported(B1, B2, F, S)
    if direct_support
        return [(F,S)]
    else
        return MOIB.added_constraint_types(inner_bridge, F, S)
    end
end

function get_mapped_fct(
    B1::Type{<:MOIBC.SetMapBridge{T}},
    B2::Type{<:MOIBC.SetMapBridge{T}},
    F, S,
    fct,
) where {T}
    direct_support, inner_bridge = how_supported(B1, B2, F, S)
    if direct_support
        return fct
    else
        return MOIBC.map_function(inner_bridge, fct)
    end
end

function get_mapped_set(
    B1::Type{<:MOIBC.SetMapBridge{T}},
    B2::Type{<:MOIBC.SetMapBridge{T}},
    F, S,
    set,
) where {T}
    direct_support, inner_bridge = how_supported(B1, B2, F, S)
    if direct_support
        return set
    else
        return MOIBC.map_set(inner_bridge, set)
    end
end

function get_bridged_and_set(
    bridge::Type{<:BoolBridge{T, B1, B2}},
    F1::Type{<:MOI.AbstractFunction},
    F2::Type{<:MOI.AbstractFunction},
    F1dim,
    F2dim,
    S1::Type{<:BoolSet{F11,F12,F11dim,F12dim,S11,S12}},
    S2::Type{<:BoolSet{F21,F22,F21dim,F22dim,S21,S22}}
) where {T, B1, B2, F11, F12,F11dim,F12dim, S11, S12, F21, F22, F21dim,F22dim, S21, S22}
    lhs_constraint_types = get_bridged_and_set(bridge, F11, F12, F11dim, F12dim, S11, S12)
    rhs_constraint_types = get_bridged_and_set(bridge, F21, F22, F21dim, F22dim, S21, S22)
    new_F1, new_F2, new_S1, new_S2 = unpack_constraint_types(lhs_constraint_types, rhs_constraint_types)
    return [(MOI.VectorAffineFunction{T},BoolSet{new_F1,new_F2,F1dim,F2dim,new_S1,new_S2})]
end

function get_bridged_and_set(
    bridge::Type{<:BoolBridge{T, B1, B2}},
    F1::Type{<:MOI.AbstractFunction},
    F2::Type{<:MOI.AbstractFunction},
    F1dim,
    F2dim,
    S1::Type{<:BoolSet{F11,F12,F11dim,F12dim,S11,S12}},
    S2
) where {T, B1, B2, F11, F12, F11dim, F12dim, S11, S12}
    lhs_constraint_types = get_bridged_and_set(bridge, F11, F12, F11dim, F12dim, S11, S12)
    rhs_constraint_types = get_constraint_types(B1, B2, F2, S2)
    new_F1, new_F2, new_S1, new_S2 = unpack_constraint_types(lhs_constraint_types, rhs_constraint_types)
    return [(MOI.VectorAffineFunction{T},BoolSet{new_F1,new_F2,F1dim,F2dim,new_S1,new_S2})]
end

function get_bridged_and_set(
    bridge::Type{<:BoolBridge{T, B1, B2}},
    F1::Type{<:MOI.AbstractFunction},
    F2::Type{<:MOI.AbstractFunction},
    F1dim,
    F2dim,
    S1,
    S2::Type{<:BoolSet{F21,F22,F21dim,F22dim,S21,S22}}
) where {T, B1, B2, F21, F22, F21dim, F22dim, S21, S22}
    lhs_constraint_types = get_constraint_types(B1, B2, F1, S1)
    rhs_constraint_types = get_bridged_and_set(bridge, F21, F22, F21dim, F22dim, S21, S22)
    new_F1, new_F2, new_S1, new_S2 = unpack_constraint_types(lhs_constraint_types, rhs_constraint_types)
    return [(MOI.VectorAffineFunction{T},BoolSet{new_F1,new_F2,F1dim,F2dim,new_S1,new_S2})]
end

function get_bridged_and_set(
    bridge::Type{<:BoolBridge{T, B1, B2}},
    F1::Type{<:MOI.AbstractFunction},
    F2::Type{<:MOI.AbstractFunction},
    F1dim,
    F2dim,
    S1,
    S2,
) where {T, B1, B2}
    lhs_constraint_types = get_constraint_types(B1, B2, F1, S1)
    rhs_constraint_types = get_constraint_types(B1, B2, F2, S2)
    new_F1, new_F2, new_S1, new_S2 = unpack_constraint_types(lhs_constraint_types, rhs_constraint_types)
    return [(MOI.VectorAffineFunction{T}, BoolSet{new_F1,new_F2,F1dim,F2dim,new_S1,new_S2})]
end

function map_function(
    bridge::Type{<:BoolBridge{T, B1, B2}},
    fct,
    set::CS.BoolSet{F1,F2,F1dim,F2dim,S1,S2}
) where {T,B1,B2,F1,F2,F1dim,F2dim,S1<:BoolSet,S2<:BoolSet}
    f = MOIU.eachscalar(fct)
    lhs_fct = map_function(bridge, f[1:get_value(F1dim)], set.lhs_set)
    rhs_fct = map_function(bridge, f[end-get_value(F2dim)+1:end], set.rhs_set)
    return MOIU.operate(vcat, T, lhs_fct, rhs_fct)
end

function map_function(
    bridge::Type{<:BoolBridge{T, B1, B2}},
    fct,
    set::CS.BoolSet{F1,F2,F1dim,F2dim,S1,S2}
) where {T,B1,B2,F1,F2,F1dim,F2dim,S1<:BoolSet,S2}
    f = MOIU.eachscalar(fct)
    lhs_fct = map_function(bridge, f[1:get_value(F1dim)], set.lhs_set)
    rhs_fct = get_mapped_fct(B1, B2, F2, S2, f[end-get_value(F2dim)+1:end])
    return MOIU.operate(vcat, T, lhs_fct, rhs_fct)
end

function map_function(
    bridge::Type{<:BoolBridge{T, B1, B2}},
    fct,
    set::CS.BoolSet{F1,F2,F1dim,F2dim,S1,S2}
) where {T,B1,B2,F1,F2,F1dim,F2dim,S1,S2<:BoolSet}
    f = MOIU.eachscalar(fct)
    lhs_fct = get_mapped_fct(B1, B2, F1, S1, f[1:get_value(F1dim)])
    rhs_fct = map_function(bridge, f[end-get_value(F2dim)+1:end], set.rhs_set)
    return MOIU.operate(vcat, T, lhs_fct, rhs_fct)
end

function map_function(
    bridge::Type{<:BoolBridge{T, B1, B2}},
    fct,
    set::CS.BoolSet{F1,F2,F1dim,F2dim,S1,S2}
) where {T,B1,B2,F1,F2,F1dim,F2dim,S1,S2}
    f = MOIU.eachscalar(fct)
    lhs_fct = get_mapped_fct(B1, B2, F1, S1, f[1:get_value(F1dim)])
    rhs_fct = get_mapped_fct(B1, B2, F2, S2, f[end-get_value(F2dim)+1:end])
    return MOIU.operate(vcat, T, lhs_fct, rhs_fct)
end

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

function map_set(
    bridge::Type{<:BoolBridge{T, B1, B2}},
    set::CS.BoolSet{F1,F2,F1dim,F2dim,S1,S2}
) where {T,B1,B2,F1,F2,F1dim,F2dim,S1<:BoolSet,S2<:BoolSet}
    lhs_set = map_set(bridge, set.lhs_set)
    rhs_set = map_set(bridge, set.rhs_set)
    return typeof_without_parmas(set){F1,F2}(lhs_set, rhs_set)
end

function map_set(
    bridge::Type{<:BoolBridge{T, B1, B2}},
    set::CS.BoolSet{F1,F2,F1dim,F2dim,S1,S2}
) where {T,B1,B2,F1,F2,F1dim,F2dim,S1<:BoolSet,S2}
    lhs_set = map_set(bridge, set.lhs_set)
    rhs_set = get_mapped_set(B1, B2, F2, S2, set.rhs_set)
    return typeof_without_parmas(set){F1,F2}(lhs_set, rhs_set)
end

function map_set(
    bridge::Type{<:BoolBridge{T, B1, B2}},
    set::CS.BoolSet{F1,F2,F1dim,F2dim,S1,S2}
) where {T,B1,B2,F1,F2,F1dim,F2dim,S1,S2<:BoolSet}
    lhs_set = get_mapped_set(B1, B2, F1, S1, set.lhs_set)
    rhs_set = map_set(bridge, set.rhs_set)
    return typeof_without_parmas(set){F1,F2}(lhs_set, rhs_set)
end

function map_set(
    bridge::Type{<:BoolBridge{T, B1, B2}},
    set::CS.BoolSet{F1,F2,F1dim,F2dim,S1,S2}
) where {T,B1,B2,F1,F2,F1dim,F2dim,S1,S2}
    lhs_set = get_mapped_set(B1, B2, F1, S1, set.lhs_set)
    rhs_set = get_mapped_set(B1, B2, F2, S2, set.rhs_set)

    return typeof_without_parmas(set){F1,F2}(lhs_set, rhs_set)
end


text:
const td = mktempdir()
flush_td() = (isdir(td) && rm(td; recursive=true); mkdir(td))
F.FOLDER_PATH[] = td

fd2html_td(e)  = fd2html(e; dir=td)
fd2html_tdv(e) = F.fd2html_v(e; dir=td)

F.def_GLOBAL_VARS!()
F.def_GLOBAL_LXDEFS!()

@testset "Paths" begin
    P = F.set_paths!()

    @test F.PATHS[:folder]   == td
    @test F.PATHS[:src]      == joinpath(td, "src")
    @test F.PATHS[:src_css]  == joinpath(td, "src", "_css")
    @test F.PATHS[:src_html] == joinpath(td, "src", "_html_parts")
    @test F.PATHS[:libs]     == joinpath(td, "libs")
    @test F.PATHS[:pub]      == joinpath(td, "pub")
    @test F.PATHS[:css]      == joinpath(td, "css")

    @test P == F.PATHS

    mkdir(F.PATHS[:src])
    mkdir(F.PATHS[:src_pages])
    mkdir(F.PATHS[:libs])
    mkdir(F.PATHS[:src_css])
    mkdir(F.PATHS[:src_html])
    mkdir(F.PATHS[:assets])
end

# copying _libs/katex in the F.PATHS[:libs] so that it can be used in testing
# the js_prerender_math
cp(joinpath(dirname(dirname(pathof(Franklin))), "test", "_libs", "katex"), joinpath(F.PATHS[:libs], "katex"))

@testset "Set vars" begin
    d = F.PageVars(
    	"a" => 0.5 => (Real,),
    	"b" => "hello" => (String, Nothing))
    F.set_vars!(d, ["a"=>"5", "b"=>"nothing"])

    @test d["a"].first == 5
    @test d["b"].first === nothing

    @test_logs (:warn, "Page var 'a' (type(s): (Real,)) can't be set to value 'blah' (type: String). Assignment ignored.") F.set_vars!(d, ["a"=>"\"blah\""])

    @test_throws F.PageVariableError F.set_vars!(d, ["a"=> "sqrt(-1)"])

    # assigning new variables

    F.set_vars!(d, ["blah"=>"1"])
    @test d["blah"].first == 1
end

@testset "Def+coms" begin # see #78
    st = raw"""
        @def title = "blah" <!-- comment -->
        @def hasmath = false
        etc
        """
    m = F.convert_md(st)
    @test F.locvar("title") == "blah"
    @test F.locvar("hasmath") == false
end


text:
"
    Take the union of 2 algorithms using &
"
struct UnionAlg{A,B} <: Algorithm
    alg1 :: A
    alg2 :: B
end

Base.:&(alg1::Algorithm,alg2::Algorithm) = UnionAlg(alg1,alg2)

function changebonds(state,alg::UnionAlg)
    state = changebonds(state,alg.alg1);
    state = changebonds(state,alg.alg2);
    return state
end

function changebonds(state,H,alg::UnionAlg,envs=environments(state,H))
    (state,envs) = changebonds(state,H,alg.alg1,envs)
    (state,envs) = changebonds(state,H,alg.alg2,envs)
    return (state,envs)
end

function find_groundstate(state,opperator,alg::UnionAlg,envs=environments(state,H))
    (state,envs) = find_groundstate(state,opperator,alg.alg1,envs);
    (state,envs,delta) = find_groundstate(state,opperator,alg.alg2,envs);
end


text:
using Documenter, Literate
using AugmentedGaussianProcesses

# Create the notebooks and examples markdown using Literate

EXAMPLES = joinpath(@__DIR__, "examples")
MD_OUTPUT = joinpath(@__DIR__, "src", "examples")
NB_OUTPUT = joinpath(@__DIR__, "src", "examples", "notebooks")

ispath(MD_OUTPUT) && rm(MD_OUTPUT; recursive=true)
ispath(NB_OUTPUT) && rm(NB_OUTPUT; recursive=true)

# add links to binder and nbviewer below the first heading of level 1
function preprocess(content)
    sub = s"""
\0
#
# [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/examples/notebooks/@__NAME__.ipynb)
"""
    return replace(content, r"^# # .*$"m => sub; count=1)
end

for file in readdir(EXAMPLES; join=true)
    endswith(file, ".jl") || continue
    Literate.markdown(file, MD_OUTPUT; documenter=true, preprocess=preprocess)
    Literate.notebook(file, NB_OUTPUT; documenter=true)
end

# Make the docs

makedocs(;
    modules=[AugmentedGaussianProcesses],
    format=Documenter.Writers.HTMLWriter.HTML(;
        assets=["assets/icon.ico"], analytics="UA-129106538-2"
    ),
    sitename="AugmentedGaussianProcesses",
    authors="Théo Galy-Fajou",
    pages=[
        "Home" => "index.md",
        "Background" => "background.md",
        "User Guide" => "userguide.md",
        "Kernels" => "kernel.md",
        "Examples" =>
            joinpath.("examples", filter(x -> endswith(x, ".md"), readdir(MD_OUTPUT))),
        "Julia GP Packages" => "comparison.md",
        "API" => "api.md",
    ],
)

# Deploy the docs

deploydocs(;
    deps=Deps.pip("mkdocs", "python-markdown-math"),
    repo="github.com/theogf/AugmentedGaussianProcesses.jl.git",
    target="build",
    push_preview=true,
)


text:
module ConferenceCall

using MacroTools: splitdef, combinedef, @capture

export @confcall, @confcall_fast

include("utils.jl")

val_value(::Type{Val{T}}) where T = T

available_vals(fn) = (m.sig.types[2] for m in methods(fn).ms)
sorted_available_vals(fn) = sort(collect(available_vals(fn)), by=val_value)

""" See Cassette#6 for why this is not Revise-able. """
@generated function call_all_methods_tuple_fast(fn, args...)
    quote
        tuple($([:(fn($v(), args...))
                 for v in sorted_available_vals(fn.instance)
                 if hasmethod(fn.instance, Tuple{v, args...})]...))
    end
end

call_all_methods_vector(fn, args...) =
    [fn(v(), args...) for v in sorted_available_vals(fn)
     if applicable(fn, v(), args...)]

impl_name(fname::Symbol) = Symbol(fname, "_impl")
function impl_name(fname::Expr)
    @assert @capture(fname, mod_.fname2_)
    return :($mod.$(impl_name(fname2)))
end

""" `@confcall_fast function foo end` """
macro confcall_fast(fn_def)
    @assert(@capture(fn_def, function fname_ end),
            "Use `@confcall` for function definitions")
    impl = impl_name(fname)
    esc(quote
        function $impl end
        $fname(args...) = $ConferenceCall.call_all_methods_tuple_fast($impl, args...)
        end)
end

macro confcall(fn_def)
    if @capture(fn_def, function fname_ end)
        impl = impl_name(fname)
        esc(quote
            function $impl end
            $fname(args...) = $ConferenceCall.call_all_methods_vector($impl, args...)
        end)
    else
        di = splitdef(fn_def)
        esc(:($ConferenceCall.@confcall $(hash(fn_def)) $fn_def))
    end
end

macro confcall(key, fn_def)
    di = splitdef(fn_def)
    di2 = copy(di)
    di2[:args] = tuple(:(::Val{$key}), di[:args]...)
    di2[:name] = impl_name(di[:name])
    return esc(combinedef(di2))
end

end # module


text:
using Test

@testset "src/ConvexPolytope.jl" begin

using BellComm: ConvexPolytope

@testset "facet_constraints()" begin
    dir = "./test/integration/files/"

    @testset "cleanup = false" begin
        test_dir = mkpath(dir * "test_porta_tmp")

        @test_throws SystemError readdir(test_dir * "porta_tmp")

        poi_file = test_dir * "/facet_constraints_test.poi"
        ieq_file = test_dir * "/facet_constraints_test.poi.ieq"

        @test !(isfile(poi_file))
        @test !(isfile(ieq_file))

        vertices = [[0;0;0],[1;0;0],[0;1;0],[1;1;1]]
        constraints = ConvexPolytope.facet_constraints(vertices, filename="facet_constraints_test", dir=test_dir, cleanup=false)
        @test constraints == ([[0 0 0 -1],[0 -1 0 1],[0 0 -1 1],[-1 1 1 -1]],[],[])

        @test isfile(poi_file)
        @test isfile(ieq_file)

        @test length(readdir(test_dir)) == 2

        rm(test_dir, recursive=true)
        @test_throws SystemError readdir(test_dir)
    end

    @testset "cleanup = true" begin

        vertices = [[0;0;0],[1;0;0],[0;1;0],[1;1;1]]
        constraints = ConvexPolytope.facet_constraints(vertices, dir=dir)
        @test constraints == (
            [[0 0 0 -1],[0 -1 0 1],[0 0 -1 1],[-1 1 1 -1]],
            [],
            []
        )

        @test_throws SystemError readdir(dir * "porta_tmp")
    end

    @testset "default filename used" begin
        test_dir = mkpath(dir * "test_porta_tmp")

        @test_throws SystemError readdir(test_dir * "porta_tmp")

        poi_file = test_dir * "/facet_constraints_tmp.poi"
        ieq_file = test_dir * "/facet_constraints_tmp.poi.ieq"

        @test !(isfile(poi_file))
        @test !(isfile(ieq_file))

        vertices = [[0;0;0],[1;0;0],[0;1;0],[1;1;1]]
        constraints = ConvexPolytope.facet_constraints(vertices, dir=test_dir, cleanup=false)
        @test constraints == (
            [[0 0 0 -1],[0 -1 0 1],[0 0 -1 1],[-1 1 1 -1]],
            [],
            []
        )

        @test isfile(poi_file)
        @test isfile(ieq_file)
        @test length(readdir(test_dir)) == 2

        rm(test_dir, recursive=true)
        @test_throws SystemError readdir(test_dir)
    end

    @testset "all defaults used" begin

        vertices = [[0;0;0],[1;0;0],[0;1;0],[1;1;1]]
        constraints = ConvexPolytope.facet_constraints(vertices)
        @test constraints == (
            [[0 0 0 -1],[0 -1 0 1],[0 0 -1 1],[-1 1 1 -1]],
            [],
            []
        )

        @test_throws SystemError readdir("./porta_tmp")
    end
end

end


text:
using AbstractPlotting
using AbstractPlotting.MakieLayout
using GLMakie

using LinearAlgebra

function bezier()
descr = """
Bezier Curves: B(t) = ∑(binomial(n,i) * tⁱ * (1-t)⁽ⁿ⁻¹⁾ * Pᵢ)
"""

## From http://juliaplots.org/MakieReferenceImages/gallery//edit_polygon/index.html:
function add_move!(scene, points, pplot)
    idx = Ref(0); dragstart = Ref(false); startpos = Base.RefValue(Point2f0(0))
    on(events(scene).mousedrag) do drag
        if ispressed(scene, Mouse.left)
            if drag == Mouse.down
                plot, _idx = mouse_selection(scene)
                if plot == pplot
                    idx[] = _idx; dragstart[] = true
                    startpos[] = to_world(scene, Point2f0(scene.events.mouseposition[]))
                end
            elseif drag == Mouse.pressed && dragstart[] && checkbounds(Bool, points[], idx[])
                pos = to_world(scene, Point2f0(scene.events.mouseposition[]))

                # very wierd, but we work with components
                # not vector
                z = zero(eltype(pos))
                x,y = pos

                ptidx = idx[]

                x = clamp(x, -1, 1)
                y = clamp(y, -1, 1)

                points[][idx[]] = [x,y]
                points[] = points[]
            end
        else
            dragstart[] = false
        end
        return
    end
end

upperpoints = Point2f0[(0,0), (5, 0), (5, 5), (0,5)]
lowerpoints = (Point2f0[(0,0), (5, 0), (5, -5), (0,-5)])

points = Node(Point2f0[(1, 4), (3, 0), (4,-4.0)])

# where we lay our scene:
scene, layout = layoutscene()
layout.halign = :left
layout.valign = :top

p = layout[1, 1:2] = LScene(scene)
rowsize!(layout, 1, Auto(1))
colsize!(layout, 1, Auto(1))
colsize!(layout, 2, Auto(1))

npts = layout[2,1:2] = LSlider(scene, range=3:12, startvalue=4)
layout[3,1:2] = LText(scene, chomp(descr))

# points = Node(Point2f0[(-1/2, -1/2),
#                        (-1/2, 1/2),
#                        (1/2, 1/2),
#                        (1/2, -1/2)])

#npts = 6
#ts = range(3pi/2, -pi/2, length=npts+2)
#points = Node(Point2f0[(cos(t),sin(t)) for t in ts[2:end-1]])

points = lift(npts.value) do val
    ts = range(3pi/2, -pi/2, length=val+2)
    Point2f0[(cos(t),sin(t)) for t in ts[2:end-1]]
end

bcurve = lift(points) do pts
    n = length(pts) - 1
    B = t -> begin
        tot = 0.0
        for (i′, Pᵢ) in enumerate(pts)
            i = i′ - 1
            tot += binomial(n, i) * t^i * (1-t)^(n-i) * Pᵢ
        end
        tot
    end
    ts = range(0, 1, length=200)
    Point2f0[B(t) for t in ts]
end

lines!(p.scene, bcurve, strokecolor=:black, strokewidth=15)
lines!(p.scene, points, strokecolor=:gray90, strokewidth=5, linestyle=:dash)
scatter!(p.scene, points)
xlims!(p.scene, (-1,1))
ylims!(p.scene, (-1,1))
add_move!(p.scene, points, p.scene[end])

scene

end


text:
# ******************************************************************************************
# Notices:
#
# Copyright © 2022 United States Government as represented by the Administrator of the
# National Aeronautics and Space Administration.  All Rights Reserved.
#
# Disclaimers
#
# No Warranty: THE SUBJECT SOFTWARE IS PROVIDED "AS IS" WITHOUT ANY WARRANTY OF ANY KIND,
# EITHER EXPRESSED, IMPLIED, OR STATUTORY, INCLUDING, BUT NOT LIMITED TO, ANY WARRANTY
# THAT THE SUBJECT SOFTWARE WILL CONFORM TO SPECIFICATIONS, ANY IMPLIED WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, OR FREEDOM FROM INFRINGEMENT, ANY
# WARRANTY THAT THE SUBJECT SOFTWARE WILL BE ERROR FREE, OR ANY WARRANTY THAT DOCUMENTATION,
# IF PROVIDED, WILL CONFORM TO THE SUBJECT SOFTWARE. THIS AGREEMENT DOES NOT, IN ANY MANNER,
# CONSTITUTE AN ENDORSEMENT BY GOVERNMENT AGENCY OR ANY PRIOR RECIPIENT OF ANY RESULTS,
# RESULTING DESIGNS, HARDWARE, SOFTWARE PRODUCTS OR ANY OTHER APPLICATIONS RESULTING FROM
# USE OF THE SUBJECT SOFTWARE.  FURTHER, GOVERNMENT AGENCY DISCLAIMS ALL WARRANTIES AND
# LIABILITIES REGARDING THIRD-PARTY SOFTWARE, IF PRESENT IN THE ORIGINAL SOFTWARE, AND
# DISTRIBUTES IT "AS IS."
#
# Waiver and Indemnity:  RECIPIENT AGREES TO WAIVE ANY AND ALL CLAIMS AGAINST THE UNITED
# STATES GOVERNMENT, ITS CONTRACTORS AND SUBCONTRACTORS, AS WELL AS ANY PRIOR RECIPIENT.
# IF RECIPIENT'S USE OF THE SUBJECT SOFTWARE RESULTS IN ANY LIABILITIES, DEMANDS, DAMAGES,
# EXPENSES OR LOSSES ARISING FROM SUCH USE, INCLUDING ANY DAMAGES FROM PRODUCTS BASED ON,
# OR RESULTING FROM, RECIPIENT'S USE OF THE SUBJECT SOFTWARE, RECIPIENT SHALL INDEMNIFY AND
# HOLD HARMLESS THE UNITED STATES GOVERNMENT, ITS CONTRACTORS AND SUBCONTRACTORS, AS WELL
# AS ANY PRIOR RECIPIENT, TO THE EXTENT PERMITTED BY LAW.  RECIPIENT'S SOLE REMEDY FOR ANY
# SUCH MATTER SHALL BE THE IMMEDIATE, UNILATERAL TERMINATION OF THIS AGREEMENT.
# ******************************************************************************************

"""
    MCTS <: LocalSolver

Monte Carlo tree search algorithm.
"""
Base.@kwdef mutable struct MCTS <: LocalSolver
    num_iterations::Int64      = 1000
    top_k::Int64               = 10
    k::Float64                 = 1.0
    α::Float64                 = 0.7
    c::Float64                 = 1.0
    tree::Union{Node, Nothing} = nothing
end

"""
    MCTSResult <: LocalResult

Monte Carlo tree search result.
"""
mutable struct MCTSResult <: LocalResult
    path::Vector
end

"""
Roll out random policy from given state, returning total return and final state.
"""
function rollout(mdp::CommonRLInterface.AbstractEnv, s::Node)
    d = terminated(mdp)
    a = rand(actions(mdp))
    r = act!(mdp, a)
    d && return r, s

    s′ = add(s, a; forward=false)
    rf, sf = rollout(mdp, s′)
    return r + rf, sf
end

"""
Perform one full simulation of MCTS.
"""
function simulate(mcts::MCTS, mdp::CommonRLInterface.AbstractEnv, s::Node=mcts.tree)
    s.n += 1

    # expansion
    s.n == 1 && return rollout(mdp, s)

    if length(s.transitions) < mcts.k * s.n ^ mcts.α
        # progressive widening
        a = rand(actions(mdp))
        s′ = add(s, a; forward=false)
    else
        # maximize upper confidence bound
        a, s′ = top_transition(s)
    end

    # transition
    d = terminated(mdp)
    r = act!(mdp, a)
    d && return r, s′

    # update value estimate
    rf, sf = simulate(mcts, mdp, s′)
    q = r + rf
    s′.q = ((s′.n - 1) * s′.q + q) / (s′.n)

    # update upper confidence bound
    s.transitions[a => s′] = s′.q + mcts.c * sqrt(log(s.n) / s′.n)
    return q, sf
end

function Solvers.solve(mcts::MCTS, env_fn::Function)
    mdp = env_fn()
    A = typeof(rand(actions(mdp)))
    mcts.tree = Node{A}()
    best_paths = PriorityQueue()

    @showprogress for _ in 1:mcts.num_iterations
        reset!(mdp)
        r, s = simulate(mcts, mdp)
        path = trace(s)
        enqueue!(best_paths, MCTSResult(path), r, mcts.top_k)
    end

    return best_paths
end


text:
using Pkg
Pkg.add("Images")
Pkg.add("ImageMetadata")
Pkg.add("ImageView")
Pkg.add("TestImages")
Pkg.add("ImageMagick")
Pkg.update()

# validation
using Images, ImageMetadata, TestImages, ImageView, ImageMagick

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
export interpolate_y_l

"""
Interpolate l (vectors of Lagrangian multipliers in ADMM-based algorithms)
and y (auxiliary variable vectors) to finer grids
"""
function interpolate_y_l(
                        l                ::Vector{Vector{TF}},
                        y                ::Vector{Vector{TF}},
                        set_Prop_levels  ::Vector{Any},
                        comp_grid_levels ::Vector{Any},
                        dim3             ::Bool,
                        i                ::Integer
                        ) where {TF<:Real}

  for j=1:length(l) #loop over the various constraint sets

    if set_Prop_levels[i].tag[j][2] == "TV" || set_Prop_levels[i].tag[j][2]=="D2D" || set_Prop_levels[i].tag[j][2]=="D3D"#this tag contains the transform-domain operator type as a string
      if dim3 #use 3D
        p1e = (comp_grid_levels[i+1].n[1]-1)*comp_grid_levels[i+1].n[2]*comp_grid_levels[i+1].n[3]
        p2e = (comp_grid_levels[i+1].n[1]-1)*comp_grid_levels[i+1].n[2]*comp_grid_levels[i+1].n[3] + comp_grid_levels[i+1].n[1]*(comp_grid_levels[i+1].n[2]-1)*comp_grid_levels[i+1].n[3]

        y_part_1=y[j][1:p1e]
        y_part_2=y[j][p1e+1:p2e]
        y_part_3=y[j][p2e+1:end]

        l_part_1=l[j][1:p1e]
        l_part_2=l[j][p1e+1:p2e]
        l_part_3=l[j][p2e+1:end]

        itp_l_p1   = interpolate(reshape(l_part_1,comp_grid_levels[i+1].n[1]-1,comp_grid_levels[i+1].n[2],comp_grid_levels[i+1].n[3]), BSpline(Interpolations.Constant()))
        itp_y_p1   = interpolate(reshape(y_part_1,comp_grid_levels[i+1].n[1]-1,comp_grid_levels[i+1].n[2],comp_grid_levels[i+1].n[3]), BSpline(Interpolations.Constant()))

        itp_l_p2   = interpolate(reshape(l_part_2,comp_grid_levels[i+1].n[1],comp_grid_levels[i+1].n[2]-1,comp_grid_levels[i+1].n[3]), BSpline(Interpolations.Constant()))
        itp_y_p2   = interpolate(reshape(y_part_2,comp_grid_levels[i+1].n[1],comp_grid_levels[i+1].n[2]-1,comp_grid_levels[i+1].n[3]), BSpline(Interpolations.Constant()))

        itp_l_p3   = interpolate(reshape(l_part_3,comp_grid_levels[i+1].n[1],comp_grid_levels[i+1].n[2],comp_grid_levels[i+1].n[3]-1), BSpline(Interpolations.Constant()))
        itp_y_p3   = interpolate(reshape(y_part_3,comp_grid_levels[i+1].n[1],comp_grid_levels[i+1].n[2],comp_grid_levels[i+1].n[3]-1), BSpline(Interpolations.Constant()))

        l_1_fine = itp_l_p1(range(1,stop=comp_grid_levels[i+1].n[1]-1,length=comp_grid_levels[i].n[1]-1), range(1,stop=comp_grid_levels[i+1].n[2],length=comp_grid_levels[i].n[2]), range(1,stop=comp_grid_levels[i+1].n[3],length=comp_grid_levels[i].n[3]))
        y_1_fine = itp_y_p1(range(1,stop=comp_grid_levels[i+1].n[1]-1,length=comp_grid_levels[i].n[1]-1), range(1,stop=comp_grid_levels[i+1].n[2],length=comp_grid_levels[i].n[2]), range(1,stop=comp_grid_levels[i+1].n[3],length=comp_grid_levels[i].n[3]))

        l_2_fine = itp_l_p2(range(1,stop=comp_grid_levels[i+1].n[1],length=comp_grid_levels[i].n[1]), range(1,stop=comp_grid_levels[i+1].n[2]-1,length=comp_grid_levels[i].n[2]-1), range(1,stop=comp_grid_levels[i+1].n[3],length=comp_grid_levels[i].n[3]))
        y_2_fine = itp_y_p2(range(1,stop=comp_grid_levels[i+1].n[1],length=comp_grid_levels[i].n[1]), range(1,stop=comp_grid_levels[i+1].n[2]-1,length=comp_grid_levels[i].n[2]-1), range(1,stop=comp_grid_levels[i+1].n[3],length=comp_grid_levels[i].n[3]))

        l_3_fine = itp_l_p3(range(1,stop=comp_grid_levels[i+1].n[1],length=comp_grid_levels[i].n[1]), range(1,stop=comp_grid_levels[i+1].n[2],length=comp_grid_levels[i].n[2]), range(1,stop=comp_grid_levels[i+1].n[3]-1,length=comp_grid_levels[i].n[3]-1))
        y_3_fine = itp_y_p3(range(1,stop=comp_grid_levels[i+1].n[1],length=comp_grid_levels[i].n[1]), range(1,stop=comp_grid_levels[i+1].n[2],length=comp_grid_levels[i].n[2]), range(1,stop=comp_grid_levels[i+1].n[3]-1,length=comp_grid_levels[i].n[3]-1))

        l[j]=[vec(l_1_fine);vec(l_2_fine);vec(l_3_fine)]
        y[j]=[vec(y_1_fine);vec(y_2_fine);vec(y_3_fine)]
      else #use 2D version
        y_part_1=y[j][1:(comp_grid_levels[i+1].n[1]-1)*comp_grid_levels[i+1].n[2]]
        y_part_2=y[j][1+(comp_grid_levels[i+1].n[1]-1)*comp_grid_levels[i+1].n[2]:end]

        l_part_1=l[j][1:(comp_grid_levels[i+1].n[1]-1)*comp_grid_levels[i+1].n[2]]
        l_part_2=l[j][1+(comp_grid_levels[i+1].n[1]-1)*comp_grid_levels[i+1].n[2]:end]

        itp_l_p1   = interpolate(reshape(l_part_1,comp_grid_levels[i+1].n[1]-1,comp_grid_levels[i+1].n[2]), BSpline(Interpolations.Constant()))
        itp_y_p1   = interpolate(reshape(y_part_1,comp_grid_levels[i+1].n[1]-1,comp_grid_levels[i+1].n[2]), BSpline(Interpolations.Constant()))

        itp_l_p2   = interpolate(reshape(l_part_2,comp_grid_levels[i+1].n[1],comp_grid_levels[i+1].n[2]-1), BSpline(Interpolations.Constant()))
        itp_y_p2   = interpolate(reshape(y_part_2,comp_grid_levels[i+1].n[1],comp_grid_levels[i+1].n[2]-1), BSpline(Interpolations.Constant()))

        l_1_fine = itp_l_p1(range(1,stop=comp_grid_levels[i+1].n[1]-1,length=comp_grid_levels[i].n[1]-1), range(1,stop=comp_grid_levels[i+1].n[2],length=comp_grid_levels[i].n[2]))
        y_1_fine = itp_y_p1(range(1,stop=comp_grid_levels[i+1].n[1]-1,length=comp_grid_levels[i].n[1]-1), range(1,stop=comp_grid_levels[i+1].n[2],length=comp_grid_levels[i].n[2]))

        l_2_fine = itp_l_p2(range(1,stop=comp_grid_levels[i+1].n[1],length=comp_grid_levels[i].n[1]), range(1,stop=comp_grid_levels[i+1].n[2]-1,length=comp_grid_levels[i].n[2]-1))
        y_2_fine = itp_y_p2(range(1,stop=comp_grid_levels[i+1].n[1],length=comp_grid_levels[i].n[1]), range(1,stop=comp_grid_levels[i+1].n[2]-1,length=comp_grid_levels[i].n[2]-1))

        l[j]=[vec(l_1_fine);vec(l_2_fine)]
        y[j]=[vec(y_1_fine);vec(y_2_fine)]
      end
    else #for TD-operator is identity,D_z,D_x , etc

      #compute how many gridpoints more/less the transform-domain 'image' has
      #in each dimension compared to the model grid for x and m
      s = ( comp_grid_levels[i].n .- set_Prop_levels[i].TD_n[j] )
      #y and l have same dimensions
      itp_l   = interpolate(reshape(l[j],set_Prop_levels[i+1].TD_n[j]), BSpline(Interpolations.Constant()))
      itp_y   = interpolate(reshape(y[j],set_Prop_levels[i+1].TD_n[j]), BSpline(Interpolations.Constant()))

      if dim3
        l_fine = itp_l(range(1,stop=set_Prop_levels[i+1].TD_n[j][1],length=comp_grid_levels[i].n[1]-s[1]), range(1,stop=set_Prop_levels[i+1].TD_n[j][2],length=comp_grid_levels[i].n[2]-s[2]), range(1,stop=set_Prop_levels[i+1].TD_n[j][3],length=comp_grid_levels[i].n[3]-s[3]))
        y_fine = itp_y(range(1,stop=set_Prop_levels[i+1].TD_n[j][1],length=comp_grid_levels[i].n[1]-s[1]), range(1,stop=set_Prop_levels[i+1].TD_n[j][2],length=comp_grid_levels[i].n[2]-s[2]), range(1,stop=set_Prop_levels[i+1].TD_n[j][3],length=comp_grid_levels[i].n[3]-s[3]))
      else
        l_fine = itp_l(range(1,stop=set_Prop_levels[i+1].TD_n[j][1],length=comp_grid_levels[i].n[1]-s[1]), range(1,stop=set_Prop_levels[i+1].TD_n[j][2],length=comp_grid_levels[i].n[2]-s[2]))
        y_fine = itp_y(range(1,stop=set_Prop_levels[i+1].TD_n[j][1],length=comp_grid_levels[i].n[1]-s[1]), range(1,stop=set_Prop_levels[i+1].TD_n[j][2],length=comp_grid_levels[i].n[2]-s[2]))
      end
      l[j]=vec(l_fine)
      y[j]=vec(y_fine)
    end #end if flow for TD-operator type

  end # end loop over constraint sets

return l, y
end

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
using DifferentialEquations, DiffEqParamEstim, Plots, Optim

# Monte Carlo Problem Set Up for solving set of ODEs with different initial conditions

# Set up Lotka-Volterra system
function pf_func(du,u,p,t)
  du[1] = p[1] * u[1] - p[2] * u[1]*u[2]
  du[2] = -3 * u[2] + u[1]*u[2]
end
p = [1.5,1.0]
prob = ODEProblem(pf_func,[1.0,1.0],(0.0,10.0),p)

# Setting up to solve the problem N times (for the N different initial conditions)
N = 10;
initial_conditions = [[1.0,1.0], [1.0,1.5], [1.5,1.0], [1.5,1.5], [0.5,1.0], [1.0,0.5], [0.5,0.5], [2.0,1.0], [1.0,2.0], [2.0,2.0]]
function prob_func(prob,i,repeat)
  ODEProblem(prob.f,initial_conditions[i],prob.tspan,prob.p)
end
monte_prob = MonteCarloProblem(prob,prob_func=prob_func)

# Check above does what we want
sim = solve(monte_prob,Tsit5(),num_monte=N)
plot(sim)

# Generate a dataset from these runs
data_times = 0.0:0.1:10.0
sim = solve(monte_prob,Tsit5(),num_monte=N,saveat=data_times)
data = Array(sim)

# Building a loss function
losses = [L2Loss(data_times,data[:,:,i]) for i in 1:N]

loss(sim) = sum(losses[i](sim[i]) for i in 1:N)

prob = ODEProblem(pf_func,[1.0,1.0],(0.0,10.0),[1.2,0.8])
function prob_func(prob,i,repeat)
  ODEProblem(prob.f,initial_conditions[i],prob.tspan,prob.p)
end
monte_prob = MonteCarloProblem(prob,prob_func=prob_func)
sim = solve(monte_prob,Tsit5(),num_monte=N,saveat=data_times)
loss(sim)

obj = build_loss_objective(monte_prob,Tsit5(),loss,num_monte=N,
                           saveat=data_times)

lower = zeros(2)
upper = fill(2.0,2)
result = optimize(obj, lower, upper, [1.3,0.9], Fminbox(BFGS()))

result

obj = build_loss_objective(monte_prob,Tsit5(),loss,num_monte=N,
                           abstol=1e-8,reltol=1e-8,
                           saveat=data_times)
result = optimize(obj, lower, upper, [1.3,0.9], Fminbox(BFGS()))

result

using DiffEqTutorials
DiffEqTutorials.tutorial_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
struct Scenario{R,G1,G2,G3,I,T,P,IP<:AbstractProblem{R,G1,G2,G3,I,T,P}
} <: AbstractScenario

    name::String
    probability::Float64
    parent::Union{Scenario{R,G1,G2,G3,I,T,P,IP},Nothing}
    children::Vector{Scenario{R,G1,G2,G3,I,T,P,IP}}
    investmentproblem::IP

    investments::Investments{R,G1,G2,G3}
    operations::Operations{R,G1,G2,G3,I,T,P}
    markets::Markets{R,T,P}

    periodweights::Vector{Float64}

    function Scenario{}(
        name::String,
        probability::Float64,
        parentscenario::Union{Scenario{R,G1,G2,G3,I,T,P,IP},Nothing},
        childscenarios::Vector{Scenario{R,G1,G2,G3,I,T,P,IP}},
        investmentproblem::IP,
        investments::Investments{R,G1,G2,G3},
        operations::Operations{R,G1,G2,G3,I,T,P},
        markets::Markets{R,T,P},
        periodweights::Vector{Float64}
) where {R,G1,G2,G3,I,T,P,IP<:AbstractProblem{R,G1,G2,G3,I,T,P}}

        @assert length(periodweights) == P
        @assert 0 < probability <= 1

        s = new{R,G1,G2,G3,I,T,P,IP}(
            name, probability, parentscenario, childscenarios,
            investmentproblem, investments, operations, markets, periodweights)

        isnothing(parentscenario) || push!(parentscenario.children, s)

        return s

    end

end

function Scenario(
    parent::Scenario{R,G1,G2,G3,I,T,P,IP},
    name::String, p::Float64,
    investments::Investments{R,G1,G2,G3},
    operations::Operations{R,G1,G2,G3,I,T,P},
    markets::Markets{R,T,P},
    weights::Vector{Float64}
) where {R,G1,G2,G3,I,T,P,IP}

    s = Scenario(
        name, p, parent, Scenario{R,G1,G2,G3,I,T,P,IP}[],
        parent.investmentproblem, investments, operations, markets, weights)

    return setup!(s)

end

function Scenario(
    invprob::IP,
    name::String,
    investments::Investments{R,G1,G2,G3},
    operations::Operations{R,G1,G2,G3,I,T,P},
    markets::Markets{R,T,P},
    weights::Vector{Float64}
) where {R,G1,G2,G3,I,T,P,IP<:AbstractProblem{R,G1,G2,G3,I,T,P}}

    s = Scenario(name, 1.0, nothing, Scenario{R,G1,G2,G3,I,T,P,IP}[], invprob,
                 investments, operations, markets, weights)

    return setup!(s)

end

function loadscenario(
    regions::Dict{String,Int}, techs::Technologies,
    periods::Dict{String,Int}, n_timesteps::Int,
    scenariofolder::String)

    resourcepath = joinpath(scenariofolder, "resources")
    investments = loadinvestments(techs, regions, resourcepath)
    operations = loadoperations(
        techs, regions, periods, n_timesteps, resourcepath)
    markets = loadmarkets(
        regions, periods, n_timesteps, joinpath(scenariofolder, "markets"))

    return investments, operations, markets

end

function setup!(s::Scenario)

    ip = s.investmentproblem
    m = ip.model

    invs = s.investments
    ops = s.operations
    markets = s.markets
    weights = s.periodweights

    # Investments
    setup!(s, :thermal, m, ip.initialconditions.thermal)
    setup!(s, :variable, m, ip.initialconditions.variable)
    setup!(s, :storage, m, ip.initialconditions.storage)

    # Operations
    setup!(ops.thermal,  ip.technologies.thermal,   m, invs.thermal, weights, s)
    setup!(ops.variable, ip.technologies.variable,  m, invs.variable, weights, s)
    setup!(ops.storage,     ip.technologies.storage,   m, invs.storage, weights, s)
    setup!(ops.transmission, ip.technologies.interface, m)

    # Markets
    setup!(markets.capacity,     m, ops, s)
    setup!(markets.rec,          m, ops, weights, s)
    setup!(markets.energy,       m, ops, weights, s)
    setup!(markets.raisereserve, m, ops, weights, s)
    setup!(markets.lowerreserve, m, ops, weights, s)

    return s

end

function maturing(
    s::Scenario, r::Int, g::Int,
    invtype::Symbol, leadtime::Symbol, action::Symbol)

    stepsback = 0
    count = ExpressionRef()
    historical = s

    while !isnothing(historical)

        invs = getfield(historical.investments, invtype)
        if getfield(invs, leadtime)[r,g] == stepsback
            count += getfield(invs, action)[r,g]
        end

        historical = historical.parent
        stepsback += 1

    end

    return count

end

function welfare(s::Scenario)

    discountrate = s.investmentproblem.discountrate
    oneoffwelfare = welfare(s.investments)
    recurringwelfare = welfare(s.operations) + welfare(s.markets)

    return if length(s.children) > 0
        # Current + discounted future welfare
        oneoffwelfare + recurringwelfare + discountrate *
        sum(cs.probability * welfare(cs) for cs in s.children)
    else
        # Repeat recurring welfare forever
        oneoffwelfare + recurringwelfare / (1 - discountrate)
    end

end

function discount(s::Scenario)

    haschildren = length(s.children) > 0
    discountrate = s.investmentproblem.discountrate

    discountings = 0
    while s.parent != nothing
        discountings += 1
        s = s.parent
    end

    oneoffrate = discountrate ^ discountings
    recurringrate = haschildren ? oneoffrate : oneoffrate / (1 - discountrate)

    return oneoffrate, recurringrate

end


text:
module Utils

function is_prime(n::Int)
    if n == 2
        return true
    else
        for i = 2:Int(ceil(sqrt(n)))
            if n%i==0
                return false
            end
        end
        return true
    end
end

function next_prime(n::Int)
    i = n+1
    while !is_prime(i)
        i += 1
    end
    return i
end

# Prime number iterator.
struct Primes end
Base.start(::Primes) = 2
Base.next(P::Primes, state) = (state, next_prime(state))
Base.done(P::Primes, state) = false

function factorize(n::Int)
    # Factorizes n = p1^(e1) * p2^(e2) * ... * pk^(ek)
    # into array of tuples [(p1,e1),(p2,e2),...,(pk,ek)]
    # We'll do this by iteratively dividing n by ever
    # larger primes (and powers thereof) until n becomes 1,
    # at which point we'll know factorization is complete.
    factors = []

    for prime in Primes()
        if n % prime == 0 # We've found a factor
            # Find the power to which the prime
            # factor is raised by progressively dividing
            # n by the factor until it no longer divides
            # evenly.
            power = 0
            while n % prime == 0
                n /= prime
                power += 1
            end
            push!(factors, (prime, power))
        end
        if n == 1
            break # We've fully factored n
        end
    end
    return factors
end

end