text:
## icon support

## Theme points to a directory :default...

abstract Icon
type StockIcon <: Icon
    nm::Union(Nothing, Symbol)
    theme::Union(Nothing, Symbol)
end

type FileIcon <: Icon
    file::String
end

## constructor
icon(nm::Symbol, theme::Symbol=nothing) = StockIcon(nm, theme)
icon(nm::Symbol) = StockIcon(nm, nothing)
function icon(nm::String)
    isfile(nm) || error("Not a file name: $nm")
    FileIcon(nm)
end


text:
using Test, GeometryBasics

@testset "algorithms.jl" begin
    cube = Rect(Vec3f0(-0.5), Vec3f0(1))
    cube_faces = decompose(TriangleFace{Int}, faces(cube))
    cube_vertices = decompose(Point{3,Float32}, cube)
    @test area(cube_vertices, cube_faces) == 6
    mesh = Mesh(cube_vertices, cube_faces)
    @test GeometryBasics.volume(mesh) ≈ 1
end

@testset "Cylinder" begin
    @testset "constructors" begin
        o, extr, r = Point2f0(1, 2), Point2f0(3, 4), 5.0f0
        s = Cylinder(o, extr, r)
        @test typeof(s) == Cylinder{2,Float32}
        @test typeof(s) == Cylinder2{Float32}
        @test origin(s) == o
        @test extremity(s) == extr
        @test radius(s) == r
        #@test abs(height(s)- norm([1,2]-[3,4]))<1e-5
        h = norm(o - extr)
        @test isapprox(height(s), h)
        #@test norm(direction(s) - Point{2,Float32}([2,2]./norm([1,2]-[3,4])))<1e-5
        @test isapprox(direction(s), Point2f0(2, 2) ./ h)
        v1 = rand(Point{3,Float64})
        v2 = rand(Point{3,Float64})
        R = rand()
        s = Cylinder(v1, v2, R)
        @test typeof(s) == Cylinder{3,Float64}
        @test typeof(s) == Cylinder3{Float64}
        @test origin(s) == v1
        @test extremity(s) == v2
        @test radius(s) == R
        @test height(s) == norm(v2 - v1)
        #@test norm(direction(s) - Point{3,Float64}((v2-v1)./norm(v2-v1)))<1e-10
        @test isapprox(direction(s), (v2 - v1) ./ norm(v2 .- v1))
    end

    @testset "decompose" begin

        o, extr, r = Point2f0(1, 2), Point2f0(3, 4), 5.0f0
        s = Cylinder(o, extr, r)
        positions = Point{3,Float32}[(-0.7677671, 3.767767, 0.0),
                                     (2.767767, 0.23223293, 0.0),
                                     (0.23223293, 4.767767, 0.0),
                                     (3.767767, 1.2322329, 0.0), (1.2322329, 5.767767, 0.0),
                                     (4.767767, 2.232233, 0.0)]
        @test decompose(Point3f0, Tesselation(s, (2, 3))) ≈ positions

        FT = TriangleFace{Int}
        faces = FT[(1, 2, 4), (1, 4, 3), (3, 4, 6), (3, 6, 5)]
        @test faces == decompose(FT, Tesselation(s, (2, 3)))

        v1 = Point{3,Float64}(1, 2, 3)
        v2 = Point{3,Float64}(4, 5, 6)
        R = 5.0
        s = Cylinder(v1, v2, R)
        positions = Point{3,Float64}[(4.535533905932738, -1.5355339059327373, 3.0),
                                     (7.535533905932738, 1.4644660940672627, 6.0),
                                     (3.0412414523193148, 4.041241452319315,
                                      -1.0824829046386295),
                                     (6.041241452319315, 7.041241452319315,
                                      1.9175170953613705),
                                     (-2.535533905932737, 5.535533905932738,
                                      2.9999999999999996),
                                     (0.46446609406726314, 8.535533905932738, 6.0),
                                     (-1.0412414523193152, -0.04124145231931431,
                                      7.0824829046386295),
                                     (1.9587585476806848, 2.9587585476806857,
                                      10.08248290463863), (1, 2, 3), (4, 5, 6)]

        @test decompose(Point3{Float64}, Tesselation(s, 8)) ≈ positions

        faces = TriangleFace{Int}[(3, 2, 1), (4, 2, 3), (5, 4, 3), (6, 4, 5), (7, 6, 5),
                                  (8, 6, 7), (1, 8, 7), (2, 8, 1), (3, 1, 9), (2, 4, 10),
                                  (5, 3, 9), (4, 6, 10), (7, 5, 9), (6, 8, 10), (1, 7, 9),
                                  (8, 2, 10)]
        @test faces == decompose(TriangleFace{Int}, Tesselation(s, 8))

        m = triangle_mesh(Tesselation(s, 8))

        @test GeometryBasics.faces(m) == faces
        @test GeometryBasics.coordinates(m) ≈ positions
        m = normal_mesh(s)# just test that it works without explicit resolution parameter
        @test m isa GLNormalMesh

        muv = uv_mesh(s)
        @test Rect(Point.(texturecoordinates(muv))) == FRect2D(Vec2f0(0), Vec2f0(1.0))
    end
end

@testset "HyperRectangles" begin
    a = Rect(Vec(0, 0), Vec(1, 1))
    pt_expa = Point{2,Int}[(0, 0), (1, 0), (0, 1), (1, 1)]
    @test decompose(Point{2,Int}, a) == pt_expa
    mesh = normal_mesh(a)
    @test decompose(Point2f0, mesh) == pt_expa

    b = Rect(Vec(1, 1, 1), Vec(1, 1, 1))
    pt_expb = Point{3,Int64}[[1, 1, 1], [1, 1, 2], [1, 2, 2], [1, 2, 1], [1, 1, 1],
                             [2, 1, 1], [2, 1, 2], [1, 1, 2], [1, 1, 1], [1, 2, 1],
                             [2, 2, 1], [2, 1, 1], [2, 2, 2], [1, 2, 2], [1, 1, 2],
                             [2, 1, 2], [2, 2, 2], [2, 1, 2], [2, 1, 1], [2, 2, 1],
                             [2, 2, 2], [2, 2, 1], [1, 2, 1], [1, 2, 2]]
    @test decompose(Point{3,Int}, b) == pt_expb
    mesh = normal_mesh(b)
end

NFace = NgonFace

@testset "Faces" begin
    @test convert_simplex(GLTriangleFace, QuadFace{Int}(1, 2, 3, 4)) ==
          (GLTriangleFace(1, 2, 3), GLTriangleFace(1, 3, 4))
    @test convert_simplex(NFace{3,ZeroIndex{Int}}, QuadFace{ZeroIndex{Int}}(1, 2, 3, 4)) ==
          (NFace{3,ZeroIndex{Int}}(1, 2, 3), NFace{3,ZeroIndex{Int}}(1, 3, 4))
    @test convert_simplex(NFace{3,OffsetInteger{3,Int}},
                          NFace{4,OffsetInteger{2,Int}}(1, 2, 3, 4)) ==
          (NFace{3,OffsetInteger{3,Int}}(1, 2, 3), NFace{3,OffsetInteger{3,Int}}(1, 3, 4))
    @test convert_simplex(LineFace{Int}, QuadFace{Int}(1, 2, 3, 4)) ==
          (LineFace{Int}(1, 2), LineFace{Int}(2, 3), LineFace{Int}(3, 4),
           LineFace{Int}(4, 1))

    @testset "NgonFace ambiguity" begin
        face = NgonFace((1, 2))
        @test convert_simplex(NgonFace{2,UInt32}, face) === (NgonFace{2,UInt32}((1, 2)),)
        @test convert_simplex(typeof(face), face) === (face,)
        face = NgonFace((1,))
        @test convert_simplex(NgonFace{1,UInt32}, face) === (NgonFace{1,UInt32}((1,)),)
        @test convert_simplex(typeof(face), face) === (face,)
    end
end

@testset "Normals" begin
    n64 = Vec{3,Float64}[(0.0, 0.0, -1.0), (0.0, 0.0, -1.0), (0.0, 0.0, -1.0),
                         (0.0, 0.0, -1.0), (0.0, 0.0, 1.0), (0.0, 0.0, 1.0),
                         (0.0, 0.0, 1.0), (0.0, 0.0, 1.0), (-1.0, 0.0, 0.0),
                         (-1.0, 0.0, 0.0), (-1.0, 0.0, 0.0), (-1.0, 0.0, 0.0),
                         (1.0, 0.0, 0.0), (1.0, 0.0, 0.0), (1.0, 0.0, 0.0), (1.0, 0.0, 0.0),
                         (0.0, 1.0, 0.0), (0.0, 1.0, 0.0), (0.0, 1.0, 0.0), (0.0, 1.0, 0.0),
                         (0.0, -1.0, 0.0), (0.0, -1.0, 0.0), (0.0, -1.0, 0.0),
                         (0.0, -1.0, 0.0),]
    n32 = map(Vec{3,Float32}, n64)
    r = triangle_mesh(centered(Rect3D))
    # @test normals(coordinates(r), GeometryBasics.faces(r)) == n32
    # @test normals(coordinates(r), GeometryBasics.faces(r)) == n64
end

@testset "HyperSphere" begin
    sphere = Sphere{Float32}(Point3f0(0), 1.0f0)

    points = decompose(Point, Tesselation(sphere, 3))
    point_target = Point{3,Float32}[[0.0, 0.0, 1.0], [1.0, 0.0, 6.12323e-17],
                                    [1.22465e-16, 0.0, -1.0], [-0.0, 0.0, 1.0],
                                    [-1.0, 1.22465e-16, 6.12323e-17],
                                    [-1.22465e-16, 1.49976e-32, -1.0], [0.0, -0.0, 1.0],
                                    [1.0, -2.44929e-16, 6.12323e-17],
                                    [1.22465e-16, -2.99952e-32, -1.0]]
    @test points ≈ point_target

    f = decompose(TriangleFace{Int}, Tesselation(sphere, 3))
    face_target = TriangleFace{Int}[[1, 2, 5], [1, 5, 4], [2, 3, 6], [2, 6, 5], [4, 5, 8],
                                    [4, 8, 7], [5, 6, 9], [5, 9, 8]]
    @test f == face_target
    circle = Circle(Point2f0(0), 1.0f0)
    points = decompose(Point2f0, Tesselation(circle, 20))
    @test length(points) == 20
    tess_circle = Tesselation(circle, 32)
    mesh = triangle_mesh(tess_circle)
    @test decompose(Point2f0, mesh) ≈ decompose(Point2f0, tess_circle)
end

@testset "Rectangles" begin
    rect = FRect2D(0, 7, 20, 3)
    @test (rect + 4) == FRect2D(4, 11, 20, 3)
    @test (rect + Vec(2, -2)) == FRect2D(2, 5, 20, 3)

    @test (rect - 4) == FRect2D(-4, 3, 20, 3)
    @test (rect - Vec(2, -2)) == FRect2D(-2, 9, 20, 3)

    base = Vec3f0(1, 2, 3)
    wxyz = Vec3f0(-2, 4, 2)
    rect = FRect3D(base, wxyz)
    @test (rect + 4) == FRect3D(base .+ 4, wxyz)
    @test (rect + Vec(2, -2, 3)) == FRect3D(base .+ Vec(2, -2, 3), wxyz)

    @test (rect - 4) == FRect3D(base .- 4, wxyz)
    @test (rect - Vec(2, -2, 7)) == FRect3D(base .- Vec(2, -2, 7), wxyz)

    rect = FRect2D(0, 7, 20, 3)
    @test (rect * 4) == FRect2D(0, 7 * 4, 20 * 4, 3 * 4)
    @test (rect * Vec(2, -2)) == FRect2D(0, -7 * 2, 20 * 2, -3 * 2)

    base = Vec3f0(1, 2, 3)
    wxyz = Vec3f0(-2, 4, 2)
    rect = FRect3D(base, wxyz)
    @test (rect * 4) == FRect3D(base .* 4, wxyz .* 4)
    @test (rect * Vec(2, -2, 3)) == FRect3D(base .* Vec(2, -2, 3), wxyz .* Vec(2, -2, 3))

    rect1 = Rect(Vec(0.0, 0.0), Vec(1.0, 2.0))
    rect2 = Rect(0.0, 0.0, 1.0, 2.0)
    @test rect1 isa GeometryBasics.HyperRectangle{2,Float64}
    @test rect1 == rect2

    split1, split2 = GeometryBasics.split(rect1, 2, 1)
    @test widths(split1) == widths(split2)
    @test origin(split1) == Vec(0, 0)
    @test origin(split2) == Vec(0, 1)
    @test in(split1, rect1)
    @test !in(rect1, split1)

    prim = Rect(0.0, 0.0, 1.0, 1.0)
    @test length(prim) == 2

    @test width(prim) == 1.0
    @test height(prim) == 1.0

    b1 = Rect2D(0.0, 0.0, 2.0, 2.0)
    b2 = Rect2D(0, 0, 2, 2)
    @test isequal(b1, b2)

    pt = Point(1.0, 1.0)
    b1 = Rect(0.0, 0.0, 1.0, 1.0)
    @test in(pt, b1)

    rect = Rect(0.0, 0.0, 1.0, 1.0)
    @test GeometryBasics.positive_widths(rect) isa GeometryBasics.HyperRectangle{2,Float64}

    h1 = Rect(0.0, 0.0, 1.0, 1.0)
    h2 = Rect(1.0, 1.0, 2.0, 2.0)
    @test union(h1, h2) isa GeometryBasics.HyperRectangle{2,Float64}
    @test GeometryBasics.diff(h1, h2) == h1
    @test GeometryBasics.intersect(h1, h2) isa GeometryBasics.HyperRectangle{2,Float64}

    b = Rect(0.0, 0.0, 1.0, 1.0)
    v = Vec(1, 2)
    @test update(b, v) isa GeometryBasics.HyperRectangle{2,Float64}
    v = Vec(1.0, 2.0)
    @test update(b, v) isa GeometryBasics.HyperRectangle{2,Float64}

    p = Vec(5.0, 4.0)
    rect = Rect(0.0, 0.0, 1.0, 1.0)
    @test min_dist_dim(rect, p, 1) == 4.0
    @test min_dist_dim(rect, p, 2) == 3.0
    @test max_dist_dim(rect, p, 1) == 5.0
    @test max_dist_dim(rect, p, 2) == 4.0

    rect1 = Rect(0.0, 0.0, 1.0, 1.0)
    rect2 = Rect(3.0, 1.0, 4.0, 2.0)
    @test min_dist_dim(rect1, rect2, 1) == 2.0
    @test min_dist_dim(rect1, rect2, 2) == 0.0
    @test max_dist_dim(rect1, rect2, 1) == 7.0
    @test max_dist_dim(rect1, rect2, 2) == 3.0

    @test !before(rect1, rect2)
    rect1 = Rect(0.0, 0.0, 1.0, 1.0)
    rect2 = Rect(3.0, 2.0, 4.0, 2.0)
    @test before(rect1, rect2)

    @test !meets(rect1, rect2)
    rect2 = Rect(1.0, 1.0, 4.0, 2.0)
    @test meets(rect1, rect2)

    rect1 = Rect(1.0, 1.0, 2.0, 2.0)
    rect2 = Rect(0.0, 0.0, 2.0, 1.0)
    @test !overlaps(rect1, rect2)
    rect1 = Rect(1.0, 1.0, 2.0, 2.0)
    rect2 = Rect(1.5, 1.5, 2.0, 2.0)
    @test overlaps(rect1, rect2)

    rect1 = Rect(1.0, 1.0, 2.0, 2.0)
    rect2 = Rect(0.0, 0.0, 2.0, 1.0)
    @test !GeometryBasics.starts(rect1, rect2)
    rect2 = Rect(1.0, 1.0, 1.5, 1.5)
    @test !GeometryBasics.starts(rect1, rect2)
    rect2 = Rect(1.0, 1.0, 3.0, 3.0)
    @test GeometryBasics.starts(rect1, rect2)

    rect1 = Rect(1.0, 1.0, 2.0, 2.0)
    rect2 = Rect(0.0, 0.0, 4.0, 4.0)
    @test during(rect1, rect2)
    rect1 = Rect(0.0, 0.0, 2.0, 3.0)
    rect2 = Rect(1.0, 1.0, 4.0, 2.0)
    @test !during(rect1, rect2)

    rect1 = Rect(1.0, 1.0, 2.0, 2.0)
    rect2 = Rect(0.0, 0.0, 4.0, 4.0)
    @test !finishes(rect1, rect2)
    rect1 = Rect(1.0, 0.0, 1.0, 1.0)
    rect2 = Rect(0.0, 0.0, 2.0, 1.0)
    @test !finishes(rect1, rect2)
    rect1 = Rect(1.0, 1.0, 1.0, 2.0)
    rect2 = Rect(0.0, 0.0, 2.0, 3.0)
    @test finishes(rect1, rect2)

end

ext: jl


text:
if haskey(ENV, "BUILDKITE_COMMIT") && haskey(ENV, "BUILDKITE_BRANCH")
    commit = ENV["BUILDKITE_COMMIT"]
    branch = ENV["BUILDKITE_BRANCH"]
    # Note: cluster_data_prefix is also defined in integration_tests/utils/compute_mse.jl
    cluster_data_prefix = "/central/scratch/climaci/turbulenceconvection-main"

    @info "pwd() = $(pwd())"
    @info "branch = $(branch)"
    @info "commit = $(commit)"

    using Glob
    if branch == "staging"
        commit_sha = commit[1:7]
        path = joinpath(cluster_data_prefix, commit_sha)
        mkpath(path)
        for folder_name in glob("Output.*")
            src = folder_name
            dst = joinpath(path, folder_name)
            @info "Moving $src to $dst"
            mv(src, dst; force = true)
        end
        @info "readdir(): $(readdir(path))"
    end
else
    @info "ENV keys: $(keys(ENV))"
end


text:
module JuliaInterp

using Compat, Random

include("types.jl")
include("lower.jl")
include("ast.jl")

function include(mod, path, options)
    code = read(path, String)
    tls = task_local_storage()
    prev = get(tls, :SOURCE_PATH, nothing)
    tls[:SOURCE_PATH] = path
    try
        expr::Expr = Meta.parseall(code, filename=path)
        interpret_ast(mod, expr, options)
    finally
        if prev === nothing
            delete!(tls, :SOURCE_PATH)
        else
            tls[:SOURCE_PATH] = prev
        end
    end
end

end


text:
function part1(input)
    instructions = Dict(parse_line.(eachline(input))...)
    return Int(evaluate(instructions, :a))
end

function part2(input)
    instructions = Dict(parse_line.(eachline(input))...)
    b = evaluate(instructions, :a)
    return Int(evaluate(instructions, :a, Dict(:b => b)))
end

function parse_line(line)
    operations = Dict("OR" => |,
                      "AND" => &,
                      "NOT" => ~,
                      "LSHIFT" => <<,
                      "RSHIFT" => >>)
    parts = split(line)
    output = Symbol(parts[end])
    if length(parts) > 3
        operation = operations[parts[end - 3]]
    else
        operation = identity
    end
    if length(parts) <= 4
        input = [symbol_or_number(parts[end - 2])]
    else
        input = symbol_or_number.(parts[[1, 3]])
    end
    return output => (operation, input)
end

function symbol_or_number(word)
    if '0' <= first(word) <= '9'
        return parse(UInt16, word)
    end
    return Symbol(word)
end

evaluate(instructions, value::UInt16, cache) = value

function evaluate(instructions, variable::Symbol,
                  cache = Dict{Symbol, UInt16}())
    if !haskey(cache, variable)
        operation, inputs = instructions[variable]
        values = evaluate.(Ref(instructions), inputs, Ref(cache))
        cache[variable] = operation(values...)
    end
    return cache[variable]
end


text:
# Task creation/generation

module ArbiterTasks

export ArbiterTask

import ..ArbiterGraphs: ImmutableNodeSet

immutable ArbiterTask
    name::Symbol
    func::Function
    dependencies::ImmutableNodeSet
end

"""
Create a task object

name: The name of the task.
function: The actual task function. It should take no arguments,
    and return a False-y value if it fails.
dependencies: (optional, ()) Any dependencies that this task relies
    on.
"""
function ArbiterTask(name, func::Function, dependencies::Any=())
    return ArbiterTask(symbol(name), func, ImmutableNodeSet([symbol(d) for d in dependencies]))
end

end


text:
module UnstaggeredMesh1D1P
export Mesh

struct Mesh
  x
  p
end

end


text:
immutable FirstTime{Process, Pred}
    process::Process
    pred::Pred
end

function sim{Process, Pred}(ft::FirstTime{Process, Pred}, t, k)
  s=size(vwrap(initial(ft.process)))
  gen = Generator(t, s, k)

  times=Vector{Float64}()
  mask=falses(k)

  for i in 1:gen.n
    next(ft.process, gen, i)
    t = gen.t[i]

    # todo: we can reduce k after predicate hits.
    for j in 1:gen.k
      if !mask[j] && ft.pred(t, gen.y[j])
        mask[j] = true
        push!(times, t)
      end
    end
  end

  times
end


text:
include("Path.jl")
include("Reading.jl")
include("Table.jl")
include("Plots.jl")
include("Option.jl")

using DataFrames
#=
	For New Zealand soil classification
	Webb abd Liburne Criteria for defining the soil family and soil sibling
=#
function START_PSD()
	println("=== START RUNNING PSD TEXTURE === \n")
		∑Psd, Id, N_∑Psd = reading.PSD(path.Input)

		Sand = Array{Float64}(undef, N_∑Psd)
		Silt = Array{Float64}(undef, N_∑Psd)
		Clay = Array{Float64}(undef, N_∑Psd)
		SoilClass = Array{String}(undef, N_∑Psd)

		for iSoil in 1:N_∑Psd
			Clay[iSoil] =  ∑Psd[iSoil, 1]
			Silt[iSoil] = ∑Psd[iSoil, 5] - Clay[iSoil]
			Sand[iSoil] = ∑Psd[iSoil,10] - Silt[iSoil]

			ClaySiltSand = Clay[iSoil] + Silt[iSoil] + Sand[iSoil]

			if 0.9999 < ClaySiltSand < 1.0001
				error("ID = $iSoil  Clay+Sand+Silt≠1")
			end  # if name
		end #  iSoil in 1:N_∑Psd

		if option.Classification == "NZsimple"
			N_Clayey = 0
			N_Silty = 0
			N_Loamy = 0
			N_Sandy = 0
			N_Error = 0

		elseif option.Classification == "NZ"
			N_Sand = 0
			N_LoamySand = 0
			N_SandyLoam = 0
			N_LoamySilt = 0
			N_Silt = 0
			N_SandyClayLoam = 0
			N_ClayLoam = 0
			N_SiltLoam = 0
			N_SiltyClay = 0
			N_LoamyClay = 0
			N_Clay = 0
			N_Error = 0
		end

		for iSoil in 1:N_∑Psd
			if option.Classification == "NZsimple"
				if Clay[iSoil] >= 0.35
					SoilClass[iSoil] = "Clayey"
					N_Clayey += 1

				elseif Clay[iSoil] < 0.35 && Silt[iSoil] >= 0.40
					SoilClass[iSoil] = "Silty"
					N_Silty += 1

				elseif Clay[iSoil] < 0.35 && Silt[iSoil] < 0.40  && Clay[iSoil] >= 0.08
					SoilClass[iSoil] = "Loamy"
					N_Loamy += 1

				elseif Clay[iSoil] < 0.35 && Silt[iSoil] < 0.40  && Clay[iSoil] < 0.08
					SoilClass[iSoil] = "Sandy"
					N_Sandy += 1

				else
					SoilClass[iSoil] = "Error"
					N_Error += 1
				end

			elseif option.Classification == "NZ" # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
				if Clay[iSoil] <= 0.08 && Sand[iSoil] >= 0.8
					SoilClass[iSoil] = "Sand"
					N_Sand += 1

				elseif Clay[iSoil] <= 0.08 && Silt[iSoil] < 0.40 &&  Sand[iSoil] < 0.8
					SoilClass[iSoil] = "LoamySand"
					N_LoamySand += 1

				elseif Clay[iSoil] > 0.08 && Silt[iSoil] < 0.40
					SoilClass[iSoil] = "SandyLoam"
					N_SandyLoam += 1

				elseif 0.4 < Silt[iSoil] < 0.82
					SoilClass[iSoil] = "LoamySilt"
					N_LoamySilt += 1

				elseif 0.82 < Silt[iSoil]
					SoilClass[iSoil] = "Silt"
					N_Silt += 1

				elseif Silt[iSoil] <= 0.15
					SoilClass[iSoil] = "SandyClayLoam"
					N_SandyClayLoam += 1

				elseif 0.15 < Silt[iSoil] < 0.4
					SoilClass[iSoil] = "ClayLoam"
					N_ClayLoam += 1

				elseif 0.4 <= Silt[iSoil]
					SoilClass[iSoil] = "SiltLoam"
					N_SiltLoam += 1

				elseif Clay[iSoil] < 0.60 && 0.30 < Silt[iSoil]
					SoilClass[iSoil] = "SiltyClay"
					N_SiltyClay += 1

				elseif Clay[iSoil] < 0.60 && Silt[iSoil] < 0.3
					SoilClass[iSoil] = "LoamyClay"
					N_LoamyClay += 1

				elseif Clay[iSoil] >= 0.60
					SoilClass[iSoil] = "Clay"
					N_Clay += 1

				else
					N_Error += 1
				end

			end
		end # for iSoil in 1:N_∑Psd

		# OUTPUTS
			table.Output(path.Output, ∑Psd, Id, N_∑Psd, Sand, Silt, Clay, SoilClass)
			emptyDF = DataFrame(N_Clayey=Int64[], N_Silty=Int64[], N_Loamy=Int64[], N_Sandy=Int64[], N_Error=Int64[])

			if option.Classification == "NZsimple"
				Statistics_N = DataFrame(N_Clayey=N_Clayey, N_Silty=N_Silty, N_Loamy=N_Loamy, N_Sandy=N_Sandy, N_Error=N_Error)
				println(Statistics_N, "\n")

				emptyDF = DataFrame(Sand=Float64[], Silt=Float64[], Clay=Float64[])
				Statistics_SandSiltClay = DataFrame(Sand=Sand, Silt=Silt, Clay=Clay)
				println(describe(Statistics_SandSiltClay))

			elseif option.Classification == "NZ"
				Statistics_N = DataFrame(N_Sand = N_Sand, N_LoamySand = N_LoamySand, N_SandyLoam = N_SandyLoam, N_LoamySilt = N_LoamySilt, N_Silt = N_Silt, N_SandyClayLoam = N_SandyClayLoam, N_ClayLoam = N_ClayLoam,  N_SiltLoam = N_SiltLoam, N_SiltyClay = N_SiltyClay, N_LoamyClay = N_LoamyClay, N_Clay = N_Clay, N_Error=N_Error)
				println(Statistics_N, "\n")
			end

			# Work in progress
			plots.TERNARY(path.Plots, Sand, Silt, Clay)

		println("\n === END RUNNING PSD TEXTURE ===")

end # function START_PSD()

START_PSD()

max_stars_repo_name: joseph-pollacco/TextureTriangleNZ

max_issues_repo_name: joseph-pollacco/TextureTriangleNZ

max_forks_repo_name: joseph-pollacco/TextureTriangleNZ


text:
f1=open("words_alpha.txt")
words=Set{String}()
for w in readlines(f1)
    push!(words,w)
end
println("# of words: $(length(words))");
f2=open("p059_cipher.txt")
cipher=read(f2,String) # if `String` is not supplied, it will be considered as a json-like array
# cipher=map(x->parse(Int32,x),convert(Array{String},split(cipher,",")))
cipher=map(split(cipher,",")) do x
    parse(Int32,x)
end
t=zeros(Int32,3)
ans=(0,"")
for t[1] in 'a':'z'
    for t[2] in 'a':'z'
        for t[3] in 'a':'z'
            dec=convert(Array{Char},map(x->cipher[x]⊻t[x%3+1],1:length(cipher)))|>String
            cnt_word=0
            tot_spl=0
            sp=split(dec," ")
            for u in sp
                cnt_word+=in(u,words)
                tot_spl+=1
            end
            global ans=max(ans,(cnt_word,dec)) # https://docs.julialang.org/en/v1/manual/variables-and-scoping/, stupid design imo
        end
    end
end
println("$ans")
println("sum of ascii: $(sum(map(x->Int(x),collect(ans[2]))))")


text:
# call each time on activating this package
using Revise

##  for Travis - must be a better way ...
if ! ( "src/" in LOAD_PATH )
    push!( LOAD_PATH, "src/")
end

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
module MyExample2

include("extra_file.jl")

end


text:
module TestInitialValues

include("preamble.jl")
using InitialValues: Init, hasinitialvalue

@testset begin
    @test hasinitialvalue(push!!)
    @test hasinitialvalue(append!!)
    @test push!!(Init(push!!), 1) == [1]
    @test append!!(Init(append!!), [1]) == [1]
    @test append!!([1], Init(append!!)) == [1]
end

end  # module


text:
using DashBootstrapComponents

nav = dbc_nav(
    [
        dbc_navitem(dbc_navlink("Active", active=true, href="#")),
        dbc_navitem(dbc_navlink("A link", href="#")),
        dbc_navitem(dbc_navlink("Another link", href="#")),
        dbc_navitem(dbc_navlink("Disabled", disabled=true, href="#")),
    ],
    vertical="md",
);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
# function generate_random_pose{T<:Number}(n::Number, range::UnitRange{T})
#    vec = [Pose2D(rand(range),rand(range),deg2rad(rand(0:360))) for p in 1:n]
#    return vec
# end
#
function is_rotation(m::AbstractMatrix{T})

end

max_stars_repo_name: elsuizo/Robot

max_issues_repo_name: elsuizo/Robot

max_forks_repo_name: elsuizo/Robot


text:
# AIML Interpreter

using EzXML, DataFrames

# Build graph DataFrame for quick access
# Need to learn more about tree in memory
function parseaiml(fileName)
    doc = readxml(fileName)
    aiml = root(doc)

    for item in eachelement(aiml)
        for phrase in eachelement(item)
            ph = strip(nodecontent(phrase))
            println(ph)
        end
    end
end

# Preprocess user input
function parseprompt(prompt)
    matchPrompt(uppercase(prompt))
    # Define default sentence structures to fit? KZ
end

# Find pattern matching processed user input and choose proper handler
function matchprompt(prompt)
end

parseaiml("basic.aiml")

max_stars_repo_name: seantfish/chatbot

max_issues_repo_name: seantfish/chatbot

max_forks_repo_name: seantfish/chatbot


text:
using StaticArrays
import ACE
import ACE: evaluate, evaluate_d, evaluate_dd, valtype, gradtype,
            write_dict, read_dict,
            DState,
            rrule_evaluate!, frule_evaluate!

using LinearAlgebra: I, norm

# TODO:
#   - retire GetNorm
#   - polish GetVal, and expand to GetVali
#   - polish the frule and rrule implementations

# ------------------ Some different ways to produce an argument

abstract type StateTransform end
abstract type StaticGet <: StateTransform end

ACE.evaluate(fval::StaticGet, X) = getval(X, fval)
ACE.evaluate_d(fval::StaticGet, X) = getval_d(X, fval)

valtype(fval::StaticGet, X) = typeof(evaluate(fval, X))
gradtype(fval::StaticGet, X) = typeof(evaluate_d(fval, X))

struct GetVal{VSYM} <: StaticGet end

getval(X, ::GetVal{VSYM}) where {VSYM} = getproperty(X, VSYM)

_one(x::Number) = one(x)
_one(x::SVector{3, T}) where {T}  = SMatrix{3, 3, T}(I)

getval_d(X, ::GetVal{VSYM}) where {VSYM} =
      DState( NamedTuple{(VSYM,)}( (_one(getproperty(X, VSYM)),) ) )

get_symbols(::GetVal{VSYM}) where {VSYM} = (VSYM,)

function rrule_evaluate!(dB, dP, ::GetVal{VSYM}, X) where {VSYM}
   x = getproperty(X, VSYM)
   TDX = eltype(dB)
   for n = 1:length(dB)
      dB[n] = TDX(  DState( NamedTuple{(VSYM,)}( ( dP[n], ) ) ) )
   end
   return dB
end

grad_type_dP(TDP, ::GetVal{VSYM}, X) where {VSYM} =
      typeof(DState( NamedTuple{(VSYM,)}( (zero(TDP),) ) ))

function rrule_evaluate(dP, ::GetVal{VSYM}, X) where {VSYM}
   x = getproperty(X, VSYM)
   return [ DState( NamedTuple{(VSYM,)}( ( dP[n], ) ) )
            for n = 1:length(dP) ]
end

# TODO - this is incomplete for now
# struct GetVali{VSYM, IND} <: StaticGet end
# getval(X, ::GetVali{VSYM, IND}) where {VSYM, IND} = getproperty(X, VSYM)[IND]
# getval_d(X, ::GetVali{VSYM, IND}) where {VSYM, IND} = __e(getproperty(X, VSYM), Val{IND}())

struct GetNorm{VSYM} <: StaticGet end

getval(X, ::GetNorm{VSYM}) where {VSYM} = norm(getproperty(X, VSYM))

function getval_d(X, ::GetNorm{VSYM}) where {VSYM}
   x = getproperty(X, VSYM)
   return DState( NamedTuple{(VSYM,)}( (x/norm(x),) ) )
end

function evaluate_dd(::GetNorm{VSYM}, X) where {VSYM}
   𝐫 = getproperty(X,VSYM)
   r = norm(𝐫)
   𝐫̂ = 𝐫 / r
   ddx = (I - 𝐫̂ * 𝐫̂') / r
   return DState( NamedTuple{(VSYM,)}( (ddx,) ) )
end

get_symbols(::GetNorm{VSYM}) where {VSYM} = (VSYM,)

write_dict(fval::StaticGet) = Dict("__id__" => "ACE_StaticGet",
                                   "expr" => string(typeof(fval)) )

read_dict(::Val{:ACE_StaticGet}, D::Dict) = eval( Meta.parse(D["expr"]) )()

function rrule_evaluate!(dB, dP, ::GetNorm{VSYM}, X) where {VSYM}
   x = getproperty(X, VSYM)
   dx = x/norm(x)
   TDX = eltype(dB)
   for n = 1:length(dB)
      dB[n] = TDX(  DState( NamedTuple{(VSYM,)}( ( dx * dP[n], ) ) ) )
   end
   return dB
end

grad_type_dP(TDP, ::GetNorm{VSYM}, X) where {VSYM} =
      typeof(DState( NamedTuple{(VSYM,)}( (zero(SVector{3,TDP}),) ) ))


text:
function init_pseudo_observable_mappings!(m::myModel1010depo)

    pseudo_names = if subspec(m) in ["ss2", "ss4", "ss22"]
	    [:y_t, :y_f_t, :OutputGap,
	     :π_t, :LongRunInflation, :MarginalCost,
	     :Wages, :FlexibleWages, :z_t, :Hours, :FlexibleHours,
	     :RealNaturalRate, :ExAnteRealRate, :NominalNaturalRate,
	     :NominalFFR, :ExpectedAvgNominalNaturalRate, :NominalRateGap,
	     :ExpectedAvg10YearRealRate,    :ExpectedAvg10YearRealNaturalRate,
	     :ExpectedAvg10YearNominalRate, :ExpectedAvg10YearNominalNaturalRate,:ExpectedAvg10YearRateGap,
	     :ExpectedAvg5YearNominalRate,  :ExpectedAvg5YearNominalNaturalRate,
	     :ExpectedAvg5YearRealRate,     :ExpectedAvg5YearRealNaturalRate, :ExpectedAvg5YearRateGap,
	     :ExpectedAvg20YearNominalRate, :ExpectedAvg20YearNominalNaturalRate,
	     :ExpectedAvg20YearRealRate,    :ExpectedAvg20YearRealNaturalRate,   :ExpectedAvg20YearRateGap,
	     :Forward5YearRateGap,
	     :Forward10YearRateGap,
	     :Forward20YearRealRate, :Forward20YearRealNaturalRate,
	     :Forward30YearRealRate, :Forward30YearRealNaturalRate,
		 :RealDepositRate, :NominalRate, :NominalDepRate, :gdp, :gdi, :tfp, :InflationPCE,
		 :GDPdeflator, :BBBspread, :AAAspread, :LongRate, :LongRealRate, :LongInflation]
    else
        [:y_t, :y_f_t, :OutputGap,
         :Hours,
         :RealNaturalRate, :ExAnteRealRate,
         :NominalNaturalRate,
         :NominalFFR,
         :ExpectedAvg5YearRealRate, :ExpectedAvg5YearRealNaturalRate,
         :ExpectedAvg10YearRealRate,:ExpectedAvg10YearRealNaturalRate,
         :ExpectedAvg20YearRealRate, :ExpectedAvg20YearRealNaturalRate,
         :Forward5YearRealRate,  :Forward5YearRealNaturalRate,
         :Forward10YearRealRate, :Forward10YearRealNaturalRate,
         :Forward20YearRealRate, :Forward20YearRealNaturalRate,
         :Forward30YearRealRate, :Forward30YearRealNaturalRate,
		 :RealDepositRate, :NominalRate, :NominalDepRate, :gdp, :gdi, :tfp, :InflationPCE,
		 :GDPdeflator, :BBBspread, :AAAspread, :LongRate, :LongRealRate, :LongInflation]
    end

    if subspec(m) in ["ss22"]
        to_add = [:g_t, :b_liq_t, :b_safe_t, :μ_t, :z_t, :λ_f_t, :λ_w_t, :rm_t, :σ_ω_t, :μ_e_t,
                  :γ_t, :π_star_t, :lr_t, :tfp_t, :e_gdpdef_t, :e_corepce_t, :e_gdp_t, :e_gdi_t,
                  :e_BBB_t, :e_AAA_t]
        if subspec(m) in ["ss22"]
            to_add = setdiff(to_add, [:z_t])
        end
        push!(pseudo_names, to_add...)
    end

    if subspec(m) in ["ss13", "ss18"]
        push!(pseudo_names, :LiquidityConvenienceYield, :SafetyConvenienceYield,
              :Forward20YearLiquidityConvenienceYield, :Forward20YearSafetyConvenienceYield,
              :SDF, :Forward20YearSDF)
    end

	# Create PseudoObservable objects
	pseudo = OrderedDict{Symbol,PseudoObservable}()
	for k in pseudo_names
	    pseudo[k] = PseudoObservable(k)
	end

    # Fill in names and reverse transforms
	pseudo[:y_t].name     = "Output Growth"
	pseudo[:y_t].longname = "Output Growth Per Capita"

################################################################################
# The following variables are added to compare the model estimate of the variables
# estimated with measurement error with their observed counterparts
	pseudo[:gdp].name     	   = "GDP growth"
	pseudo[:gdp].longname      = "Model-implied GDP Growth"
	pseudo[:gdp].rev_transform = loggrowthtopct_annualized_percapita

	pseudo[:gdi].name 		   = "GDI growth"
	pseudo[:gdi].longname	   = "Model-implied GDI Growth"
	pseudo[:gdi].rev_transform = loggrowthtopct_annualized_percapita

	pseudo[:tfp].name          = "Model-implied TFP"
	pseudo[:tfp].longname      = "Model-implied Fernald's TFP"
	pseudo[:tfp].rev_transform = quartertoannual

	pseudo[:InflationPCE].name          = "PCE inflation"
	pseudo[:InflationPCE].longname      = "Core PCE Inflation"
	pseudo[:InflationPCE].rev_transform = loggrowthtopct_annualized

	pseudo[:GDPdeflator].name          = "GDP deflator"
	pseudo[:GDPdeflator].longname      = "GDP Deflator Inflation"
	pseudo[:GDPdeflator].rev_transform = loggrowthtopct_annualized

	pseudo[:BBBspread].name          = "Baa spread"
	pseudo[:BBBspread].longname      = "Baa - 20-year Treasury Spread"
	pseudo[:BBBspread].rev_transform = quartertoannual

	pseudo[:AAAspread].name 		 = "Aaa spread"
	pseudo[:AAAspread].longname 	 = "Aaa - 20-year Treasury Spread"
	pseudo[:AAAspread].rev_transform = quartertoannual

	pseudo[:LongRate].name 	 	    = "10y nominal yield"
	pseudo[:LongRate].longname 		= "10-year Nominal Bond Yield"
	pseudo[:LongRate].rev_transform = quartertoannual

	pseudo[:LongRealRate].name 	 	    = "10y real yield"
	pseudo[:LongRealRate].longname 		= "10-year Real Bond Yield"
	pseudo[:LongRealRate].rev_transform = quartertoannual

	pseudo[:LongInflation].name 	 	    = "..."
	pseudo[:LongInflation].longname 		= "..."
	pseudo[:LongInflation].rev_transform    = loggrowthtopct_annualized
################################################################################

	pseudo[:y_f_t].name = "Flexible Output Growth"
	pseudo[:y_f_t].longname = "Output that would prevail in a flexible-price economy."

	pseudo[:OutputGap].name = "Output Gap"
	pseudo[:OutputGap].longname = "Output Gap"

    if subspec(m) in ["ss2", "ss4", "ss22"]
	    pseudo[:π_t].name = "Inflation"
	    pseudo[:π_t].longname = "Inflation"
	    pseudo[:π_t].rev_transform = quartertoannual

	    pseudo[:LongRunInflation].name = "Long Run Inflation"
	    pseudo[:LongRunInflation].longname = "Long Run Inflation"
	    pseudo[:LongRunInflation].rev_transform = quartertoannual

        pseudo[:MarginalCost].name = "Marginal Cost"
        pseudo[:MarginalCost].longname = "Marginal Cost"

	    pseudo[:Wages].name = "Wages"
	    pseudo[:Wages].longname = "Wages"

	    pseudo[:FlexibleWages].name = "Flexible Wages"
	    pseudo[:FlexibleWages].longname = "Wages that would prevail in a flexible-wage economy"

	    pseudo[:z_t].name     = "z_t"
	    pseudo[:z_t].longname = "z_t"
    end

	pseudo[:Hours].name = "Hours"
	pseudo[:Hours].longname = "Hours"

    if subspec(m) in ["ss2", "ss4", "ss22"]
	    pseudo[:FlexibleHours].name     = "Flexible Hours"
	    pseudo[:FlexibleHours].longname = "Flexible Hours"
    end

	pseudo[:RealNaturalRate].name = "Real Natural Rate"
	pseudo[:RealNaturalRate].longname = "The real interest rate that would prevail in a flexible-price economy."
	pseudo[:RealNaturalRate].rev_transform = quartertoannual

	pseudo[:ExAnteRealRate].name = "Ex Ante Real Rate"
	pseudo[:ExAnteRealRate].longname = "Ex Ante Real Rate"
	pseudo[:ExAnteRealRate].rev_transform = quartertoannual

	pseudo[:NominalRate].name = "Nominal Riks-Free Rate"
	pseudo[:NominalRate].longname = "Nominal Riks-Free Rate"
	pseudo[:NominalRate].rev_transform = quartertoannual

	pseudo[:RealDepositRate].name = "Ex Ante Real Deposit Rate"
	pseudo[:RealDepositRate].longname = "Ex Ante Real Interest Rate on Mutual Funds' Deposits"
	pseudo[:RealDepositRate].rev_transform = quartertoannual

	pseudo[:NominalDepRate].name = "Nominal Deposit Rate"
	pseudo[:NominalDepRate].longname = "Nominal Interest Rate on Deposits"
	pseudo[:NominalDepRate].rev_transform = quartertoannual

	pseudo[:NominalFFR].name     = "Nominal FFR"
	pseudo[:NominalFFR].longname = "Nominal FFR at an annual rate"
	pseudo[:NominalFFR].rev_transform = quartertoannual

    if subspec(m) in ["ss2", "ss4", "ss22"]
        pseudo[:ExpectedAvgNominalNaturalRate].name     = "Expected Average Nominal Natural Rate"
        pseudo[:ExpectedAvgNominalNaturalRate].longname = "Natural Rate + Expected Inflation"
        pseudo[:ExpectedAvgNominalNaturalRate].rev_transform = quartertoannual

        pseudo[:NominalRateGap].name     = "Nominal Rate Gap"
        pseudo[:NominalRateGap].longname = "Nominal FFR - Nominal Natural Rate"
        pseudo[:NominalRateGap].rev_transform = quartertoannual
    end

    pseudo[:ExpectedAvg10YearRealRate].name     = "Expected Average 10-Year Real Rate"
    pseudo[:ExpectedAvg10YearRealRate].longname = "Expected Average 10-Year Real Interest Rate"
    pseudo[:ExpectedAvg10YearRealRate].rev_transform = quartertoannual

    pseudo[:ExpectedAvg10YearRealNaturalRate].name     = "Expected Average 10-Year Real Natural Rate"
    pseudo[:ExpectedAvg10YearRealNaturalRate].longname = "Expected Average 10-Year Real Natural Rate of Interest"
    pseudo[:ExpectedAvg10YearRealNaturalRate].rev_transform = quartertoannual

    if subspec(m) in ["ss2", "ss4", "ss22"]
        pseudo[:ExpectedAvg10YearNominalRate].name     = "Expected Average 10-Year Nominal Rate"
        pseudo[:ExpectedAvg10YearNominalRate].longname = "Expected Average 10-Year Nominal Interest Rate"
        pseudo[:ExpectedAvg10YearNominalRate].rev_transform = quartertoannual

        pseudo[:ExpectedAvg10YearNominalNaturalRate].name     = "Expected Average 10-Year Nominal Natural Rate"
        pseudo[:ExpectedAvg10YearNominalNaturalRate].longname = "Expected Average 10-Year Nominal Natural Rate of Interest"
        pseudo[:ExpectedAvg10YearNominalNaturalRate].rev_transform = quartertoannual

        pseudo[:ExpectedAvg10YearRateGap].name     = "Expected Average 10-Year Rate Gap"
        pseudo[:ExpectedAvg10YearRateGap].longname = "Expected Average 10-Year Rate Gap"
        pseudo[:ExpectedAvg10YearRateGap].rev_transform = quartertoannual
    end

    pseudo[:ExpectedAvg5YearRealRate].name     = "Expected Average 5-Year Real Rate"
    pseudo[:ExpectedAvg5YearRealRate].longname = "Expected Average 5-Year Real Interest Rate"
    pseudo[:ExpectedAvg5YearRealRate].rev_transform = quartertoannual

    pseudo[:ExpectedAvg5YearRealNaturalRate].name     = "Expected Average 5-Year Real Natural Rate"
    pseudo[:ExpectedAvg5YearRealNaturalRate].longname = "Expected Average 5-Year Real Natural Rate of Interest"
    pseudo[:ExpectedAvg5YearRealNaturalRate].rev_transform = quartertoannual

    if subspec(m) in ["ss2", "ss4", "ss22"]
        pseudo[:ExpectedAvg5YearNominalRate].name     = "Expected Average 5-Year Rate"
        pseudo[:ExpectedAvg5YearNominalRate].longname = "Expected Average 5-Year Interest Rate"
        pseudo[:ExpectedAvg5YearNominalRate].rev_transform = quartertoannual

        pseudo[:ExpectedAvg5YearNominalNaturalRate].name     = "Expected Average 5-Year Natural Rate"
        pseudo[:ExpectedAvg5YearNominalNaturalRate].longname = "Expected Average 5-Year Natural Rate of Interest"
        pseudo[:ExpectedAvg5YearNominalNaturalRate].rev_transform = quartertoannual

        pseudo[:ExpectedAvg5YearRateGap].name     = "Expected Average 5-Year Rate Gap"
        pseudo[:ExpectedAvg5YearRateGap].longname = "Expected Average 5-Year Rate Gap"
        pseudo[:ExpectedAvg5YearRateGap].rev_transform = quartertoannual
    end

    pseudo[:ExpectedAvg20YearRealRate].name     = "Expected Average 20-Year Real Rate"
    pseudo[:ExpectedAvg20YearRealRate].longname = "Expected Average 20-Year Real Interest Rate"
    pseudo[:ExpectedAvg20YearRealRate].rev_transform = quartertoannual

    pseudo[:ExpectedAvg20YearRealNaturalRate].name     = "Expected Average 20-Year Real Natural Rate"
    pseudo[:ExpectedAvg20YearRealNaturalRate].longname = "Expected Average 20-Year Real Natural Rate of Interest"
    pseudo[:ExpectedAvg20YearRealNaturalRate].rev_transform = quartertoannual

    if subspec(m) in ["ss2", "ss4", "ss22"]
        pseudo[:ExpectedAvg20YearNominalRate].name     = "Expected Average 20-Year Nominal Rate"
        pseudo[:ExpectedAvg20YearNominalRate].longname = "Expected Average 20-Year Nominal Interest Rate"
        pseudo[:ExpectedAvg20YearNominalRate].rev_transform = quartertoannual

        pseudo[:ExpectedAvg20YearNominalNaturalRate].name     = "Expected Average 20-Year Nominal Natural Rate"
        pseudo[:ExpectedAvg20YearNominalNaturalRate].longname = "Expected Average 20-Year Nominal Natural Rate of Interest"
        pseudo[:ExpectedAvg20YearNominalNaturalRate].rev_transform = quartertoannual

        pseudo[:ExpectedAvg20YearRateGap].name     = "Expected Average 20-Year Rate Gap"
        pseudo[:ExpectedAvg20YearRateGap].longname = "Expected Average 20-Year Rate Gap"
        pseudo[:ExpectedAvg20YearRateGap].rev_transform = quartertoannual

        pseudo[:Forward5YearRateGap].name     = "5-Year Forward Rate Gap"
        pseudo[:Forward5YearRateGap].longname = "5-Year Forward Rate Gap"
        pseudo[:Forward5YearRateGap].rev_transform = quartertoannual

        pseudo[:Forward10YearRateGap].name     = "10-Year Forward Rate Gap"
        pseudo[:Forward10YearRateGap].longname = "10-Year Forward Rate Gap"
        pseudo[:Forward10YearRateGap].rev_transform = quartertoannual
    end

    pseudo[:Forward20YearRealRate].name     = "Forward 20-Year Real Rate"
    pseudo[:Forward20YearRealRate].longname = "Forward 20-Year Real Rate of Interest (not the average, computed by projecting TTT foreward 80 periods)"
    pseudo[:Forward20YearRealRate].rev_transform = quartertoannual

    pseudo[:Forward20YearRealNaturalRate].name     = "Forward 20-Year Real Natural Rate"
    pseudo[:Forward20YearRealNaturalRate].longname = "Forward 20-Year Real Natural Rate of Interest (not the average, computed by projecting TTT foreward 80 periods)"
    pseudo[:Forward20YearRealNaturalRate].rev_transform = quartertoannual

    pseudo[:Forward30YearRealRate].name     = "Forward 30-Year Real Rate"
    pseudo[:Forward30YearRealRate].longname = "Forward 30-Year Real Rate of Interest (not the average, computed by projecting TTT foreward 80 periods)"
    pseudo[:Forward30YearRealRate].rev_transform = quartertoannual

    pseudo[:Forward30YearRealNaturalRate].name     = "Forward 30-Year Real Natural Rate"
    pseudo[:Forward30YearRealNaturalRate].longname = "Forward 30-Year Real Natural Rate of Interest (not the average, computed by projecting TTT foreward 80 periods)"
    pseudo[:Forward30YearRealNaturalRate].rev_transform = quartertoannual

    if !(subspec(m) in ["ss2", "ss4", "ss22"])
        pseudo[:NominalNaturalRate].name     = "Nominal Natural Rate"
        pseudo[:NominalNaturalRate].longname = "Natural Rate + Expected Inflation"
        pseudo[:NominalNaturalRate].rev_transform = quartertoannual

        pseudo[:Forward5YearRealRate].name     = "Forward 5-Year Real Rate"
        pseudo[:Forward5YearRealRate].longname = "Forward 5-Year Real Rate of Interest (not the average, computed by projecting TTT foreward 80 periods)"
        pseudo[:Forward5YearRealRate].rev_transform = quartertoannual

        pseudo[:Forward5YearRealNaturalRate].name     = "Forward 5-Year Real Natural Rate"
        pseudo[:Forward5YearRealNaturalRate].longname = "Forward 5-Year Real Natural Rate of Interest (not the average, computed by projecting TTT foreward 80 periods)"
        pseudo[:Forward5YearRealNaturalRate].rev_transform = quartertoannual

        pseudo[:Forward10YearRealRate].name     = "Forward 10-Year Real Rate"
        pseudo[:Forward10YearRealRate].longname = "Forward 10-Year Real Rate of Interest (not the average, computed by projecting TTT foreward 80 periods)"
        pseudo[:Forward10YearRealRate].rev_transform = quartertoannual

        pseudo[:Forward10YearRealNaturalRate].name     = "Forward 10-Year Real Natural Rate"
        pseudo[:Forward10YearRealNaturalRate].longname = "Forward 10-Year Real Natural Rate of Interest (not the average, computed by projecting TTT foreward 80 periods)"
        pseudo[:Forward10YearRealNaturalRate].rev_transform = quartertoannual
    end

    if subspec(m) in ["ss13", "ss18"]
        pseudo[:LiquidityConvenienceYield].name = "Liquidity convenience yield"
        pseudo[:LiquidityConvenienceYield].longname = "Liquidity convenience yield"
        pseudo[:LiquidityConvenienceYield].rev_transform = quartertoannual

        pseudo[:SafetyConvenienceYield].name = "Safety convenience yield"
        pseudo[:SafetyConvenienceYield].longname = "Safety convenience yield"
        pseudo[:SafetyConvenienceYield].rev_transform = quartertoannual

        pseudo[:SDF].name = "Stochastic discount factor"
        pseudo[:SDF].longname = "Stochastic discount factor"
        pseudo[:SDF].rev_transform = quartertoannual

        pseudo[:Forward20YearLiquidityConvenienceYield].name = "20-year forward liquidity convenience yield"
        pseudo[:Forward20YearLiquidityConvenienceYield].longname = "20-year forward liquidity convenience yield"
        pseudo[:Forward20YearLiquidityConvenienceYield].rev_transform = quartertoannual

        pseudo[:Forward20YearSafetyConvenienceYield].name = "20-year forward safety convenience yield"
        pseudo[:Forward20YearSafetyConvenienceYield].longname = "20-year forward safety convenience yield"
        pseudo[:Forward20YearSafetyConvenienceYield].rev_transform = quartertoannual

        pseudo[:Forward20YearSDF].name = "20-year forward stochastic discount factor"
        pseudo[:Forward20YearSDF].longname = "20-year forward stochastic discount factor"
        pseudo[:Forward20YearSDF].rev_transform = quartertoannual
    end

    # Exogenous processes
    if subspec(m) in ["ss22"]
        for i in to_add
            pseudo[i].name = DSGE.detexify(i)
            pseudo[i].longname = DSGE.detexify(i)
        end
    end

    # Add to model object
    m.pseudo_observable_mappings = pseudo
end


text:
module Ganja

using Blink, Requires, AbstractTensors #, Grassmann
import Blink: JSString
import AbstractTensors: value

export loadganja!, loadganja, loadexample!, loadexample

function readfile(filename)
    open(joinpath(filename),"r") do f
        read(f,String)
    end
end

readexample(filename) = readfile(joinpath(@__DIR__,"..","test",filename))

#loadganja!(w) = loadjs!(w,joinpath(@__DIR__,"..","deps","ganja.js"))
loadganja!(w) = loadjs!(w,"https://unpkg.com/ganja.js")
loadganja() = (w=Window();loadganja!(w);w)

loadexample!(w,filename) = js(w,JSString(readexample(filename)))
loadexample(filename) = (w=loadganja();loadexample(w,filename);w)

function __init__()
    @require Grassmann="4df31cd9-4c27-5bea-88d0-e6a7146666d8" begin
        export addelement!
        addelement!(w,t) = js(w,JSString("addElement($(string(value(Grassmann.MultiVector(t)))))"))
    end
end

end # module


text:
using Test
using ForestBiometrics

@test sdi(200,12.3) == 278.82242465294496


text:
import DrWatson
import PowerSystems
import Dates
import JuMP
import Feather
import DataFrames
import Gurobi
import MathOptInterface
import JSON

const MOI = MathOptInterface
const MOIU = MathOptInterface.Utilities
const PSY = PowerSystems
optimizer = JuMP.with_optimizer(Gurobi.Optimizer, OutputFlag=1)

############################ Load Processed System Data ####################################
uc_system = PSY.System(joinpath(DrWatson.datadir(), "uc_system.json"))
ed_system = PSY.System(joinpath(DrWatson.datadir(), "ed_system.json"))
suc_system = PSY.System(joinpath(DrWatson.datadir(), "suc_system.json"))

include(joinpath(DrWatson.srcdir(), "experiment_utilities.jl"))
################################# Load Functions ###########################################
include(joinpath(DrWatson.srcdir(), "computation/SUC_model.jl"))
include(joinpath(DrWatson.srcdir(), "computation/ED_model.jl"))
include(joinpath(DrWatson.srcdir(), "computation/UC_model.jl"))
include(joinpath(DrWatson.srcdir(), "computation/utility_functions.jl"))
include(joinpath(DrWatson.srcdir(), "computation/trial.jl"))

################################# Results Functions ########################################
include(joinpath(DrWatson.srcdir(), "results/get_results.jl"))
include(joinpath(DrWatson.srcdir(), "results/export_results.jl"))
include(joinpath(DrWatson.srcdir(), "results/collect_results.jl"))
include(joinpath(DrWatson.srcdir(), "results/metrics.jl"))

raw_ouput_folder = prepare_workspace!(joinpath(DrWatson.datadir(), "exp_raw"))
results_date = remove_char("$(round(Dates.now(),Dates.Minute))", ":")
results_path = joinpath(raw_ouput_folder,results_date)

for ix in 1:10
    trial(ix, "UC", range(1+(ix-1)*30, length = 30), optimizer, results_path)
end

for ix in 1:10
    trial(ix, "SUC", range(1+(ix-1)*30, length = 30), optimizer, results_path)
end

fuel_cost_suc = Vector()
load_shedding_suc = Vector()
curtailment_suc = Vector()
computation_times = Vector()
for ix in 1:10
    push!(fuel_cost_suc, total_fuel_cost(results_date, ix, range(1+(ix-1)*30, length = 30), ed_system))
    push!(load_shedding_suc, total_load_not_supplied(results_date, ix,range(1+(ix-1)*30, length = 30), ed_system))
    push!(curtailment_suc, total_curtailment(results_date, ix, range(1+(ix-1)*30, length = 30), ed_system))
    push!(computation_times, total_computing_time(results_date, ix, range(1+(ix-1)*30, length = 30), "SUC"))
end

fuel_cost_uc = Vector()
load_shedding_uc = Vector()
curtailment_uc = Vector()
computation_times = Vector()
for ix in 1:10
    push!(fuel_cost_uc, total_fuel_cost(results_date, ix, range(1+(ix-1)*30, length = 30), ed_system))
    push!(load_shedding_uc, total_load_not_supplied(results_date, ix,range(1+(ix-1)*30, length = 30), ed_system))
    push!(curtailment_uc, total_curtailment(results_date, ix, range(1+(ix-1)*30, length = 30), ed_system))
    push!(computation_times, total_computing_time(results_date, ix, range(1+(ix-1)*30, length = 30), "UC"))
end


text:
### This support BLAS style multiplication
#           α * A * B + β C
# but avoids the broadcast machinery

# Lazy representation of α*A*B + β*C
struct MulAdd{StyleA, StyleB, StyleC, T, AA, BB, CC}
    α::T
    A::AA
    B::BB
    β::T
    C::CC
end

@inline MulAdd{StyleA,StyleB,StyleC}(α::T, A::AA, B::BB, β::T, C::CC) where {StyleA,StyleB,StyleC,T,AA,BB,CC} =
    MulAdd{StyleA,StyleB,StyleC,T,AA,BB,CC}(α,A,B,β,C)

@inline function MulAdd{StyleA,StyleB,StyleC}(αT, A, B, βV, C) where {StyleA,StyleB,StyleC}
    α,β = promote(αT,βV)
    MulAdd{StyleA,StyleB,StyleC}(α, A, B, β, C)
end

@inline MulAdd(α, A::AA, B::BB, β, C::CC) where {AA,BB,CC} =
    MulAdd{typeof(MemoryLayout(AA)), typeof(MemoryLayout(BB)), typeof(MemoryLayout(CC))}(α, A, B, β, C)

@inline _mul_eltype(A) = A
@inline _mul_eltype(A, B) = Base.promote_op(*, A, B)
@inline _mul_eltype(A, B, C, D...) = _mul_eltype(Base.promote_op(*, A, B), C, D...)

@inline eltype(::MulAdd{StyleA,StyleB,StyleC,T,AA,BB,CC}) where {StyleA,StyleB,StyleC,T,AA,BB,CC} =
     promote_type(_mul_eltype(T, eltype(AA), eltype(BB)), _mul_eltype(T, eltype(CC)))

size(M::MulAdd, p::Int) = size(M)[p]
axes(M::MulAdd, p::Int) = axes(M)[p]
length(M::MulAdd) = prod(size(M))
size(M::MulAdd) = map(length,axes(M))
axes(M::MulAdd) = axes(M.C)

similar(M::MulAdd, ::Type{T}, axes) where {T,N} = similar(Array{T}, axes)
similar(M::MulAdd, ::Type{T}) where T = similar(M, T, axes(M))
similar(M::MulAdd) = similar(M, eltype(M))

check_mul_axes(A) = nothing
_check_mul_axes(::Number, ::Number) = nothing
_check_mul_axes(::Number, _) = nothing
_check_mul_axes(_, ::Number) = nothing
_check_mul_axes(A, B) = axes(A,2) == axes(B,1) || throw(DimensionMismatch("Second axis of A, $(axes(A,2)), and first axis of B, $(axes(B,1)) must match"))
function check_mul_axes(A, B, C...)
    _check_mul_axes(A, B)
    check_mul_axes(B, C...)
end

# we need to special case AbstractQ as it allows non-compatiple multiplication
function check_mul_axes(A::AbstractQ, B, C...)
    axes(A.factors, 1) == axes(B, 1) || axes(A.factors, 2) == axes(B, 1) ||
        throw(DimensionMismatch("First axis of B, $(axes(B,1)) must match either axes of A, $(axes(A))"))
    check_mul_axes(B, C...)
end

function instantiate(M::MulAdd)
    @boundscheck check_mul_axes(M.α, M.A, M.B)
    @boundscheck check_mul_axes(M.β, M.C)
    @boundscheck axes(M.A,1) == axes(M.C,1) || throw(DimensionMismatch("First axis of A, $(axes(M.A,1)), and first axis of C, $(axes(M.C,1)) must match"))
    @boundscheck axes(M.B,2) == axes(M.C,2) || throw(DimensionMismatch("Second axis of B, $(axes(M.B,2)), and second axis of C, $(axes(M.C,2)) must match"))
    M
end

const ArrayMulArrayAdd{StyleA,StyleB,StyleC} = MulAdd{StyleA,StyleB,StyleC,<:Any,<:AbstractArray,<:AbstractArray,<:AbstractArray}
const MatMulVecAdd{StyleA,StyleB,StyleC} = MulAdd{StyleA,StyleB,StyleC,<:Any,<:AbstractMatrix,<:AbstractVector,<:AbstractVector}
const MatMulMatAdd{StyleA,StyleB,StyleC} = MulAdd{StyleA,StyleB,StyleC,<:Any,<:AbstractMatrix,<:AbstractMatrix,<:AbstractMatrix}
const VecMulMatAdd{StyleA,StyleB,StyleC} = MulAdd{StyleA,StyleB,StyleC,<:Any,<:AbstractVector,<:AbstractMatrix,<:AbstractMatrix}

broadcastable(M::MulAdd) = M

const BlasMatMulVecAdd{StyleA,StyleB,StyleC,T<:BlasFloat} = MulAdd{StyleA,StyleB,StyleC,T,<:AbstractMatrix{T},<:AbstractVector{T},<:AbstractVector{T}}
const BlasMatMulMatAdd{StyleA,StyleB,StyleC,T<:BlasFloat} = MulAdd{StyleA,StyleB,StyleC,T,<:AbstractMatrix{T},<:AbstractMatrix{T},<:AbstractMatrix{T}}
const BlasVecMulMatAdd{StyleA,StyleB,StyleC,T<:BlasFloat} = MulAdd{StyleA,StyleB,StyleC,T,<:AbstractVector{T},<:AbstractMatrix{T},<:AbstractMatrix{T}}

muladd!(α, A, B, β, C) = materialize!(MulAdd(α, A, B, β, C))
materialize(M::MulAdd) = copy(instantiate(M))
copy(M::MulAdd) = copyto!(similar(M), M)

@inline function copyto!(dest::AbstractArray{T}, M::MulAdd) where T
    M.C === dest  || copyto!(dest, M.C)
    muladd!(M.α, M.A, M.B, M.β, dest)
end

@inline function copyto!(dest::AbstractArray{T}, M::MulAdd{<:Any,<:Any,ZerosLayout}) where T
    α,A,B,β,C = M.α, M.A, M.B, M.β, M.C
    if !isbitstype(T) # instantiate
        dest .= β .* view(A,:,1) .* Ref(B[1])  # get shape right
    end
    muladd!(α, A, B, β, dest)
end

# Modified from LinearAlgebra._generic_matmatmul!
function tile_size(T, S, R)
    tile_size = 0
    if isbitstype(R) && isbitstype(T) && isbitstype(S)
        tile_size = floor(Int, sqrt(tilebufsize / max(sizeof(R), sizeof(S), sizeof(T))))
    end
    tile_size
end

function tiled_blasmul!(tile_size, α, A::AbstractMatrix{T}, B::AbstractMatrix{S}, β, C::AbstractMatrix{R}) where {S,T,R}
    mA, nA = size(A)
    mB, nB = size(B)
    nA == mB || throw(DimensionMismatch("Dimensions must match"))
    size(C) == (mA, nB) || throw(DimensionMismatch("Dimensions must match"))

    @inbounds begin
        sz = (tile_size, tile_size)
        # FIXME: This code is completely invalid!!!
        Atile = unsafe_wrap(Array, convert(Ptr{T}, pointer(Abuf[Threads.threadid()])), sz)
        Btile = unsafe_wrap(Array, convert(Ptr{S}, pointer(Bbuf[Threads.threadid()])), sz)

        z1 = zero(A[1, 1]*B[1, 1] + A[1, 1]*B[1, 1])
        z = convert(promote_type(typeof(z1), R), z1)

        if mA < tile_size && nA < tile_size && nB < tile_size
            copy_transpose!(Atile, 1:nA, 1:mA, 'N', A, 1:mA, 1:nA)
            copyto!(Btile, 1:mB, 1:nB, 'N', B, 1:mB, 1:nB)
            for j = 1:nB
                boff = (j-1)*tile_size
                for i = 1:mA
                    aoff = (i-1)*tile_size
                    s = z
                    for k = 1:nA
                        s += Atile[aoff+k] * Btile[boff+k]
                    end
                    C[i,j] = α*s + β*C[i,j]
                end
            end
        else
            # FIXME: This code is completely invalid!!!
            Ctile = unsafe_wrap(Array, convert(Ptr{R}, pointer(Cbuf[Threads.threadid()])), sz)
            for jb = 1:tile_size:nB
                jlim = min(jb+tile_size-1,nB)
                jlen = jlim-jb+1
                for ib = 1:tile_size:mA
                    ilim = min(ib+tile_size-1,mA)
                    ilen = ilim-ib+1
                    copyto!(Ctile, 1:ilen, 1:jlen, C, ib:ilim, jb:jlim)
                    lmul!(β,Ctile)
                    for kb = 1:tile_size:nA
                        klim = min(kb+tile_size-1,mB)
                        klen = klim-kb+1
                        copy_transpose!(Atile, 1:klen, 1:ilen, 'N', A, ib:ilim, kb:klim)
                        copyto!(Btile, 1:klen, 1:jlen, 'N', B, kb:klim, jb:jlim)
                        for j=1:jlen
                            bcoff = (j-1)*tile_size
                            for i = 1:ilen
                                aoff = (i-1)*tile_size
                                s = z
                                for k = 1:klen
                                    s += Atile[aoff+k] * Btile[bcoff+k]
                                end
                                Ctile[bcoff+i] += α*s
                            end
                        end
                    end
                    copyto!(C, ib:ilim, jb:jlim, Ctile, 1:ilen, 1:jlen)
                end
            end
        end
    end

    C
end

function default_blasmul!(α, A::AbstractMatrix, B::AbstractMatrix, β, C::AbstractMatrix)
    mA, nA = size(A)
    mB, nB = size(B)
    nA == mB || throw(DimensionMismatch("Dimensions must match"))
    size(C) == (mA, nB) || throw(DimensionMismatch("Dimensions must match"))

    @inbounds for k in colsupport(A), j in rowsupport(B)
        z2 = zero(A[k, 1]*B[1, j] + A[k, 1]*B[1, j])
        Ctmp = convert(promote_type(eltype(C), typeof(z2)), z2)
        @simd for ν = rowsupport(A,k) ∩ colsupport(B,j)
            Ctmp = muladd(A[k, ν],B[ν, j],Ctmp)
        end
        C[k,j] = muladd(α,Ctmp, β*C[k,j])
    end
    C
end

function default_blasmul!(α, A::AbstractMatrix, B::AbstractVector, β, C::AbstractVector)
    mA, nA = size(A)
    mB = length(B)
    nA == mB || throw(DimensionMismatch("Dimensions must match"))
    length(C) == mA || throw(DimensionMismatch("Dimensions must match"))

    lmul!(β, C)
    (nA == 0 || mB == 0)  && return C

    z = zero(A[1]*B[1] + A[1]*B[1])
    Astride = size(A, 1) # use size, not stride, since its not pointer arithmetic

    @inbounds for k in colsupport(B,1)
        aoffs = (k-1)*Astride
        b = B[k]
        for i = 1:mA
            C[i] += α * A[aoffs + i] * b
        end
    end

    C
end

function materialize!(M::MatMulMatAdd)
    α, A, B, β, C = M.α, M.A, M.B, M.β, M.C
    if C ≡ B
        B = copy(B)
    end
    ts = tile_size(eltype(A), eltype(B), eltype(C))
    if iszero(β) # false is a "strong" zero to wipe out NaNs
        if ts == 0 || !(axes(A) isa NTuple{2,OneTo{Int}}) || !(axes(B) isa NTuple{2,OneTo{Int}}) || !(axes(C) isa NTuple{2,OneTo{Int}})
            default_blasmul!(α, A, B, false, C)
        else
            tiled_blasmul!(ts, α, A, B, false, C)
        end
    else
        if ts == 0 || !(axes(A) isa NTuple{2,OneTo{Int}}) || !(axes(B) isa NTuple{2,OneTo{Int}}) || !(axes(C) isa NTuple{2,OneTo{Int}})
            default_blasmul!(α, A, B, β, C)
        else
            tiled_blasmul!(ts, α, A, B, β, C)
        end
    end
end

function materialize!(M::MatMulVecAdd)
    α, A, B, β, C = M.α, M.A, M.B, M.β, M.C
    if C ≡ B
        B = copy(B)
    end
    default_blasmul!(α, A, B, iszero(β) ? false : β, C)
end

# make copy to make sure always works
@inline function _gemv!(tA, α, A, x, β, y)
    if x ≡ y
        BLAS.gemv!(tA, α, A, copy(x), β, y)
    else
        BLAS.gemv!(tA, α, A, x, β, y)
    end
end

# make copy to make sure always works
@inline function _gemm!(tA, tB, α, A, B, β, C)
    if B ≡ C
        BLAS.gemm!(tA, tB, α, A, copy(B), β, C)
    else
        BLAS.gemm!(tA, tB, α, A, B, β, C)
    end
end

@inline materialize!(M::BlasMatMulVecAdd{<:AbstractColumnMajor,<:AbstractStridedLayout,<:AbstractStridedLayout}) =
    _gemv!('N', M.α, M.A, M.B, M.β, M.C)
@inline materialize!(M::BlasMatMulVecAdd{<:AbstractRowMajor,<:AbstractStridedLayout,<:AbstractStridedLayout}) =
    _gemv!('T', M.α, transpose(M.A), M.B, M.β, M.C)
@inline materialize!(M::BlasMatMulVecAdd{<:ConjLayout{<:AbstractRowMajor},<:AbstractStridedLayout,<:AbstractStridedLayout,<:BlasComplex}) =
    _gemv!('C', M.α, M.A', M.B, M.β, M.C)

@inline materialize!(M::BlasVecMulMatAdd{<:AbstractColumnMajor,<:AbstractColumnMajor,<:AbstractColumnMajor}) =
    _gemm!('N', 'N', M.α, M.A, M.B, M.β, M.C)
@inline materialize!(M::BlasVecMulMatAdd{<:AbstractColumnMajor,<:AbstractRowMajor,<:AbstractColumnMajor}) =
    _gemm!('N', 'T', M.α, M.A, transpose(M.B), M.β, M.C)
@inline materialize!(M::BlasVecMulMatAdd{<:AbstractColumnMajor,<:ConjLayout{<:AbstractRowMajor},<:AbstractColumnMajor,<:BlasComplex}) =
    _gemm!('N', 'C', M.α, M.A, M.B', M.β, M.C)

@inline materialize!(M::BlasMatMulMatAdd{<:AbstractColumnMajor,<:AbstractColumnMajor,<:AbstractColumnMajor}) =
    _gemm!('N', 'N', M.α, M.A, M.B, M.β, M.C)
@inline materialize!(M::BlasMatMulMatAdd{<:AbstractColumnMajor,<:AbstractRowMajor,<:AbstractColumnMajor}) =
    _gemm!('N', 'T', M.α, M.A, transpose(M.B), M.β, M.C)
@inline materialize!(M::BlasMatMulMatAdd{<:AbstractColumnMajor,<:ConjLayout{<:AbstractRowMajor},<:AbstractColumnMajor,<:BlasComplex}) =
    _gemm!('N', 'C', M.α, M.A, M.B', M.β, M.C)

@inline materialize!(M::BlasMatMulMatAdd{<:AbstractRowMajor,<:AbstractColumnMajor,<:AbstractColumnMajor}) =
    _gemm!('T', 'N', M.α, transpose(M.A), M.B, M.β, M.C)
@inline materialize!(M::BlasMatMulMatAdd{<:ConjLayout{<:AbstractRowMajor},<:AbstractColumnMajor,<:AbstractColumnMajor,<:BlasComplex}) =
    _gemm!('C', 'N', M.α, M.A', M.B, M.β, M.C)

@inline materialize!(M::BlasMatMulMatAdd{<:AbstractRowMajor,<:AbstractRowMajor,<:AbstractColumnMajor}) =
    _gemm!('T', 'T', M.α, transpose(M.A), transpose(M.B), M.β, M.C)
@inline materialize!(M::BlasMatMulMatAdd{<:AbstractRowMajor,<:ConjLayout{<:AbstractRowMajor},<:AbstractColumnMajor,<:BlasComplex}) =
    _gemm!('T', 'C', M.α, transpose(M.A), M.B', M.β, M.C)

@inline materialize!(M::BlasMatMulMatAdd{<:ConjLayout{<:AbstractRowMajor},<:AbstractRowMajor,<:AbstractColumnMajor,<:BlasComplex}) =
    _gemm!('C', 'T', M.α, M.A', M.B', M.β, M.C)
@inline materialize!(M::BlasMatMulMatAdd{<:ConjLayout{<:AbstractRowMajor},<:ConjLayout{<:AbstractRowMajor},<:AbstractColumnMajor,<:BlasComplex}) =
    _gemm!('C', 'C', M.α, M.A', M.B', M.β, M.C)

@inline materialize!(M::BlasMatMulMatAdd{<:AbstractColumnMajor,<:AbstractColumnMajor,<:AbstractRowMajor}) =
    _gemm!('T', 'T', M.α, M.B, M.A, M.β, transpose(M.C))
@inline materialize!(M::BlasMatMulMatAdd{<:AbstractColumnMajor,<:AbstractColumnMajor,<:ConjLayout{<:AbstractRowMajor},<:BlasComplex}) =
    _gemm!('C', 'C', M.α, M.B, M.A, M.β, M.C')

@inline materialize!(M::BlasMatMulMatAdd{<:AbstractColumnMajor,<:AbstractRowMajor,<:AbstractRowMajor}) =
    _gemm!('N', 'T', M.α, transpose(M.B), M.A, M.β, transpose(M.C))
@inline materialize!(M::BlasMatMulMatAdd{<:AbstractColumnMajor,<:AbstractRowMajor,<:ConjLayout{<:AbstractRowMajor},<:BlasComplex}) =
    _gemm!('N', 'T', M.α, transpose(M.B), M.A, M.β, M.C')
@inline materialize!(M::BlasMatMulMatAdd{<:AbstractColumnMajor,<:ConjLayout{<:AbstractRowMajor},<:ConjLayout{<:AbstractRowMajor},<:BlasComplex}) =
    _gemm!('N', 'C', M.α, M.B', M.A, M.β, M.C')

@inline materialize!(M::BlasMatMulMatAdd{<:AbstractRowMajor,<:AbstractColumnMajor,<:AbstractRowMajor}) =
    _gemm!('T', 'N', M.α, M.B, transpose(M.A), M.β, transpose(M.C))
@inline materialize!(M::BlasMatMulMatAdd{<:ConjLayout{<:AbstractRowMajor},<:AbstractColumnMajor,<:ConjLayout{<:AbstractRowMajor},<:BlasComplex}) =
    _gemm!('C', 'N', M.α, M.B, M.A', M.β, M.C')

@inline materialize!(M::BlasMatMulMatAdd{<:AbstractRowMajor,<:AbstractRowMajor,<:AbstractRowMajor}) =
    _gemm!('N', 'N', M.α, transpose(M.B), transpose(M.A), M.β, transpose(M.C))
@inline materialize!(M::BlasMatMulMatAdd{<:ConjLayout{<:AbstractRowMajor},<:ConjLayout{<:AbstractRowMajor},<:ConjLayout{<:AbstractRowMajor},<:BlasComplex}) =
    _gemm!('N', 'N', M.α, M.B', M.A', M.β, M.C')

###
# Symmetric
###

# make copy to make sure always works
@inline function _symv!(tA, α, A, x, β, y)
    if x ≡ y
        BLAS.symv!(tA, α, A, copy(x), β, y)
    else
        BLAS.symv!(tA, α, A, x, β, y)
    end
end

@inline function _hemv!(tA, α, A, x, β, y)
    if x ≡ y
        BLAS.hemv!(tA, α, A, copy(x), β, y)
    else
        BLAS.hemv!(tA, α, A, x, β, y)
    end
end

materialize!(M::BlasMatMulVecAdd{<:SymmetricLayout{<:AbstractColumnMajor},<:AbstractStridedLayout,<:AbstractStridedLayout}) =
    _symv!(symmetricuplo(M.A), M.α, symmetricdata(M.A), M.B, M.β, M.C)

materialize!(M::BlasMatMulVecAdd{<:SymmetricLayout{<:AbstractRowMajor},<:AbstractStridedLayout,<:AbstractStridedLayout}) =
    _symv!(symmetricuplo(M.A) == 'L' ? 'U' : 'L', M.α, transpose(symmetricdata(M.A)), M.B, M.β, M.C)

materialize!(M::BlasMatMulVecAdd{<:HermitianLayout{<:AbstractColumnMajor},<:AbstractStridedLayout,<:AbstractStridedLayout,<:BlasComplex}) =
    _hemv!(symmetricuplo(M.A), M.α, hermitiandata(M.A), M.B, M.β, M.C)

materialize!(M::BlasMatMulVecAdd{<:HermitianLayout{<:AbstractRowMajor},<:AbstractStridedLayout,<:AbstractStridedLayout,<:BlasComplex}) =
    _hemv!(symmetricuplo(M.A) == 'L' ? 'U' : 'L', M.α, hermitiandata(M.A)', M.B, M.β, M.C)

####
# Diagonal
####

# Diagonal multiplication never changes structure
similar(M::MulAdd{<:DiagonalLayout,<:DiagonalLayout}, ::Type{T}, axes) where T = similar(M.B, T, axes)
similar(M::MulAdd{<:DiagonalLayout}, ::Type{T}, axes) where T = similar(M.B, T, axes)
similar(M::MulAdd{<:Any,<:DiagonalLayout}, ::Type{T}, axes) where T = similar(M.A, T, axes)
# equivalent to rescaling
function materialize!(M::MulAdd{<:DiagonalLayout{<:AbstractFillLayout}})
    M.C .= (M.α * getindex_value(M.A.diag)) .* M.B .+ M.β .* M.C
    M.C
end

function materialize!(M::MulAdd{<:Any,<:DiagonalLayout{<:AbstractFillLayout}})
    M.C .= M.α .* M.A .* getindex_value(M.B.diag) .+ M.β .* M.C
    M.C
end

copy(M::MulAdd{<:DiagonalLayout{<:AbstractFillLayout}}) = (M.α * getindex_value(M.A.diag)) .* M.B .+ M.β .* M.C
copy(M::MulAdd{<:DiagonalLayout{<:AbstractFillLayout},<:Any,ZerosLayout}) = (M.α * getindex_value(M.A.diag)) .* M.B
copy(M::MulAdd{<:AbstractFillLayout,<:AbstractFillLayout,<:AbstractFillLayout}) = M.α*M.A*M.B + M.β*M.C
copy(M::MulAdd{<:Any,<:DiagonalLayout{<:AbstractFillLayout}}) = (M.α * getindex_value(M.B.diag)) .* M.A .+ M.β .* M.C
copy(M::MulAdd{<:Any,<:DiagonalLayout{<:AbstractFillLayout},ZerosLayout}) = (M.α * getindex_value(M.B.diag)) .* M.A

BroadcastStyle(::Type{<:MulAdd}) = ApplyBroadcastStyle()

scalarone(::Type{T}) where T = one(T)
scalarone(::Type{<:AbstractArray{T}}) where T = scalarone(T)
scalarzero(::Type{T}) where T = zero(T)
scalarzero(::Type{<:AbstractArray{T}}) where T = scalarzero(T)

fillzeros(::Type{T}, ax) where T = Zeros{T}(ax)

function mul!(dest::AbstractArray{W}, A::AbstractArray{T}, b::AbstractArray{V}) where {T,V,W}
    TVW = promote_type(W, _mul_eltype(T,V))
    muladd!(scalarone(TVW), A, b, scalarzero(TVW), dest)
end

function MulAdd(A::AbstractArray{T}, B::AbstractVector{V}) where {T,V}
    TV = _mul_eltype(eltype(A), eltype(B))
    MulAdd(scalarone(TV), A, B, scalarzero(TV), fillzeros(TV,(axes(A,1))))
end

function MulAdd(A::AbstractArray{T}, B::AbstractMatrix{V}) where {T,V}
    TV = _mul_eltype(eltype(A), eltype(B))
    MulAdd(scalarone(TV), A, B, scalarzero(TV), fillzeros(TV,(axes(A,1),axes(B,2))))
end

mul(A::AbstractArray, B::AbstractArray) = materialize(MulAdd(A,B))

macro layoutmul(Typ)
    ret = quote
        LinearAlgebra.mul!(dest::AbstractVector, A::$Typ, b::AbstractVector) =
            ArrayLayouts.mul!(dest,A,b)

        LinearAlgebra.mul!(dest::AbstractMatrix, A::$Typ, b::AbstractMatrix) =
            ArrayLayouts.mul!(dest,A,b)
        LinearAlgebra.mul!(dest::AbstractMatrix, A::$Typ, b::$Typ) =
            ArrayLayouts.mul!(dest,A,b)

        Base.:*(A::$Typ, B::$Typ) = ArrayLayouts.mul(A,B)
        Base.:*(A::$Typ, B::AbstractMatrix) = ArrayLayouts.mul(A,B)
        Base.:*(A::$Typ, B::AbstractVector) = ArrayLayouts.mul(A,B)
        Base.:*(A::AbstractMatrix, B::$Typ) = ArrayLayouts.mul(A,B)
        Base.:*(A::LinearAlgebra.AdjointAbsVec, B::$Typ) = ArrayLayouts.mul(A,B)
        Base.:*(A::LinearAlgebra.TransposeAbsVec, B::$Typ) = ArrayLayouts.mul(A,B)

        Base.:*(A::LinearAlgebra.AbstractQ, B::$Typ) = ArrayLayouts.lmul(A,B)
        Base.:*(A::$Typ, B::LinearAlgebra.AbstractQ) = ArrayLayouts.rmul(A,B)
    end
    for Struc in (:AbstractTriangular, :Diagonal)
        ret = quote
            $ret

            Base.:*(A::LinearAlgebra.$Struc, B::$Typ) = ArrayLayouts.mul(A,B)
            Base.:*(A::$Typ, B::LinearAlgebra.$Struc) = ArrayLayouts.mul(A,B)
        end
    end
    for Mod in (:Adjoint, :Transpose, :Symmetric, :Hermitian)
        ret = quote
            $ret

            LinearAlgebra.mul!(dest::AbstractMatrix, A::$Typ, b::$Mod{<:Any,<:AbstractMatrix}) =
                ArrayLayouts.mul!(dest,A,b)

            LinearAlgebra.mul!(dest::AbstractVector, A::$Mod{<:Any,<:$Typ}, b::AbstractVector) =
                ArrayLayouts.mul!(dest,A,b)

            Base.:*(A::$Mod{<:Any,<:$Typ}, B::$Mod{<:Any,<:$Typ}) = ArrayLayouts.mul(A,B)
            Base.:*(A::$Mod{<:Any,<:$Typ}, B::AbstractMatrix) = ArrayLayouts.mul(A,B)
            Base.:*(A::AbstractMatrix, B::$Mod{<:Any,<:$Typ}) = ArrayLayouts.mul(A,B)
            Base.:*(A::$Mod{<:Any,<:$Typ}, B::AbstractVector) = ArrayLayouts.mul(A,B)

            Base.:*(A::$Mod{<:Any,<:$Typ}, B::$Typ) = ArrayLayouts.mul(A,B)
            Base.:*(A::$Typ, B::$Mod{<:Any,<:$Typ}) = ArrayLayouts.mul(A,B)

            Base.:*(A::$Mod{<:Any,<:$Typ}, B::Diagonal) = ArrayLayouts.mul(A,B)
            Base.:*(A::Diagonal, B::$Mod{<:Any,<:$Typ}) = ArrayLayouts.mul(A,B)

            Base.:*(A::LinearAlgebra.AbstractTriangular, B::$Mod{<:Any,<:$Typ}) = ArrayLayouts.mul(A,B)
            Base.:*(A::$Mod{<:Any,<:$Typ}, B::LinearAlgebra.AbstractTriangular) = ArrayLayouts.mul(A,B)
        end
    end

    esc(ret)
end


text:
module Ifpack2
  using CxxWrap, MPI
  import ..libjltrilinos
  import ..Teuchos
  import ..Tpetra

  @wrapmodule(libjltrilinos, :register_ifpack2)

  function default_parameters(ifpack_prec = "ILUT")
    pl = Teuchos.ParameterList("Ifpack2")
    ifpack2 = Teuchos.ParameterList("Ifpack2")
    pl["Preconditioner Type"] = ifpack_prec
    pl[ifpack_prec] = Teuchos.ParameterList(ifpack_prec)
  end

end


text:
module HDGElasticity

using LinearAlgebra
using SparseArrays
using StaticArrays
using IntervalArithmetic
using Roots
using CartesianMesh
using PolynomialBasis
using ImplicitDomainQuadrature
import ImplicitDomainQuadrature: extend

IDQ = ImplicitDomainQuadrature

include("affine_map.jl")
include("utils.jl")
include("variational_form_utils.jl")
include("dg_mesh.jl")
include("fit_interface_hybrid_element.jl")
include("function_space.jl")
include("isotropic_elasticity.jl")
include("local_operator.jl")
include("local_hybrid_coupling.jl")
include("hybrid_coupling.jl")
include("hybrid_local_coupling.jl")
include("local_solver.jl")
include("assembly.jl")

end # module


text:
# Model

function init_model()
    demog = init_demographics();
    endow = init_test_endowments(demog.lastWorkAge);
    kGrid = init_test_capital_grid();
    budget = init_test_budget(demog.lastWorkAge);
    util = UtilityCRRA(2.0);
    discFactor = 0.97;
    m = Model(demog, endow, kGrid, budget, util, discFactor);
    @assert validate_model(m);
    return m
end

function validate_model(m :: Model)
    isValid = validate_demographics(m.demog) &&
        validate_endowments(m.endow) &&
        validate_capital_grid(m.kGrid) &&
        validate_budget(m.budget);
    return isValid
end

## ----------  Helpers

n_efficiencies(m :: Model, t :: Integer) =
    isretired(m, t)  ?  1  :  n_efficiencies(m.endow);
betar(m :: Model) = m.discFactor * m.budget.intRate;
euler_dev(m :: Model, ctV) = UtilityFunctions890.euler_dev(m.util, ctV, betar(m));

function kprime_range(m :: Model, t :: Integer, k, eff)
    kPrimeMin = kprime_min(m.kGrid, t);
    kPrimeMax = min(kprime_max(m.budget, t, k, eff), kprime_max(m.kGrid, t));
    return kPrimeMin, kPrimeMax
end

# A little trick -- the same as writing
# `lifespan(m :: Model) = lifespan(m.demog)` etc.
Lazy.@forward Model.demog (
    lifespan, workspan, isretired
);

Lazy.@forward Model.kGrid (
    make_k_grid, k_min, k_max, kgrid_min, kgrid_max
);

Lazy.@forward Model.budget (
    consumption, wage
);

# -------------


text:
@testset "559.maximum-depth-of-n-ary-tree.jl" begin
    @test max_depth_559([[3, 2, 4], Int[], [5, 6], Int[], Int[], Int[]]) == 3
    @test max_depth_559([[2, 3, 4, 5], Int[], [6, 7], [8], [9, 10], Int[], [11], [12], [13], Int[], [14], Int[], Int[], Int[]]) == 5
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 311bed50-ee23-4546-9140-a8b67e76d755
# setup package environment over Tethys packages
begin
    import Pkg
    Pkg.activate(Base.current_project())

    using Tethys

	using Plots, StatsPlots
end

# ╔═╡ 89a4e2a6-d7eb-11ea-3b22-fb24e94f2aa9
using BenchmarkTools

# ╔═╡ da1bb8e9-92b0-4c1d-a48d-fa1046f441c7
gr(format=:png, size=(800,600)) # Force raster / PNG, as otherwise millions of data points get sent to web browser ^_^

# ╔═╡ 141e94f8-7362-11eb-3606-6360ade0fefc
Ĥ=Tethys.FrohlichHamiltonian(α=5.0, μ=-6)

# ╔═╡ 9a8c77a0-a593-11ea-09e0-5177021a936e
begin
	MAX_ORDER=2
	diag1=Tethys.Diag(0,1.0,MAX_ORDER,[0.1,0.2,1.0 , 0.3,0.4,1.0], [0.0,0.0,0.0,0.0,0.0,0.0])
	GF1=Tethys.BareExpansion(diag1, Ĥ)
end

# ╔═╡ 6c3748a8-a5a7-11ea-2490-d5e6dce272a8
diag1

# ╔═╡ ddcd5f98-a593-11ea-1b0f-6f6067269fc4
GF2=Tethys.Monte!(diag1, Ĥ)

# ╔═╡ 69b81526-a5a7-11ea-272f-21121dbbfbb9
diag1

# ╔═╡ ed0ad5da-a593-11ea-0d19-1fa974170c9d
begin
	diag=Tethys.Diag(0,1.0,MAX_ORDER,[0.1,0.2,1.0 , 0.3,0.4,1.0], [0,0,0,0,0,0])
	MC=[]
	for i in 1:1_000_000
		GF=Tethys.Monte!(diag, Ĥ)
		append!(MC,GF)
	end
end

# ╔═╡ 6bef50b2-a594-11ea-3ab0-dd19b06da250
MC

# ╔═╡ 6f29880e-a594-11ea-3a1e-b18693818f43
plot(MC)

# ╔═╡ e8d22018-a593-11ea-2d09-c921536e9287
density(MC)

# ╔═╡ d08f9624-a599-11ea-3af6-5f6dc915480f
diag

# ╔═╡ d16e91d0-a599-11ea-23f1-1963e0d078bb
Tethys.BareExpansion(diag, Ĥ)

# ╔═╡ abd6ba19-16c4-4844-98bc-b1c7ee5716bd
# Benchmarking below here...

# ╔═╡ a93ba57a-d7eb-11ea-26c8-1180ccb530d5
@benchmark Tethys.Monte!(diag1, H)

# ╔═╡ 4bd42778-d7ec-11ea-095b-8f80fca95014
@benchmark Tethys.BareExpansion(diag, H)

# ╔═╡ c2e29b42-7362-11eb-10e1-e1653cdd521e
@benchmark(exp(2.0))

# ╔═╡ Cell order:
# ╠═311bed50-ee23-4546-9140-a8b67e76d755
# ╠═da1bb8e9-92b0-4c1d-a48d-fa1046f441c7
# ╠═141e94f8-7362-11eb-3606-6360ade0fefc
# ╠═9a8c77a0-a593-11ea-09e0-5177021a936e
# ╠═6c3748a8-a5a7-11ea-2490-d5e6dce272a8
# ╠═ddcd5f98-a593-11ea-1b0f-6f6067269fc4
# ╠═69b81526-a5a7-11ea-272f-21121dbbfbb9
# ╠═ed0ad5da-a593-11ea-0d19-1fa974170c9d
# ╠═6bef50b2-a594-11ea-3ab0-dd19b06da250
# ╠═6f29880e-a594-11ea-3a1e-b18693818f43
# ╠═e8d22018-a593-11ea-2d09-c921536e9287
# ╠═d08f9624-a599-11ea-3af6-5f6dc915480f
# ╠═d16e91d0-a599-11ea-23f1-1963e0d078bb
# ╠═abd6ba19-16c4-4844-98bc-b1c7ee5716bd
# ╠═89a4e2a6-d7eb-11ea-3b22-fb24e94f2aa9
# ╠═a93ba57a-d7eb-11ea-26c8-1180ccb530d5
# ╠═4bd42778-d7ec-11ea-095b-8f80fca95014
# ╠═c2e29b42-7362-11eb-10e1-e1653cdd521e


text:
using TopOpt, Zygote, FiniteDifferences, LinearAlgebra, Test, Random, SparseArrays
const FDM = FiniteDifferences
using Ferrite: ndofs_per_cell, getncells
using TopOpt: ndofs
using Arpack
include("..//truss_topopt_problems//utils.jl")

gm_ins_dir = joinpath(@__DIR__, "..", "truss_topopt_problems", "instances", "ground_meshes");

@testset "AssembleK" begin
    nels = (2, 2)
    problem = PointLoadCantilever(Val{:Quadratic}, nels, (1.0, 1.0), 1.0, 0.3, 1.0)
    ak = AssembleK(problem)
    dh = problem.ch.dh
    total_ndof = ndofs(dh)
    T = eltype(problem.E)
    einfo = ElementFEAInfo(problem)
    k = size(einfo.Kes[1], 1)
    N = length(einfo.Kes)
    for _ in 1:3
        v = rand(T, total_ndof)
        f = Kx -> sum(ak(Kx) * v)
        Kes = [rand(T, k, k) for _ in 1:N]
        Kes .= transpose.(Kes) .+ Kes
        val1, grad1 = NonconvexCore.value_gradient(f, Kes)
        val2, grad2 = f(Kes), Zygote.gradient(f, Kes)[1]
        grad3 = FDM.grad(central_fdm(5, 1), f, Kes)[1]
        @test val1 == val2
        @test norm(grad1 - grad2) == 0
        map(1:length(grad2)) do i
            g1 = grad2[i]
            _g2 = grad3[i]
            g2 = (_g2' + _g2) / 2
            @test norm(g1 - g2) <= 1e-4
        end
    end
end

@testset "ElementK" begin
    nels = (2, 2)
    problem = PointLoadCantilever(Val{:Quadratic}, nels, (1.0, 1.0), 1.0, 0.3, 1.0)
    solver = FEASolver(Direct, problem; xmin=0.01, penalty=TopOpt.PowerPenalty(1.0))

    ek = ElementK(solver)
    dh = problem.ch.dh
    T = eltype(problem.E)
    N = getncells(dh.grid)
    k = ndofs_per_cell(dh)

    # * check stiffness matrix consistency
    Kes_1 = ek(ones(T, prod(nels)))
    for (ci, (k1, k0)) in enumerate(zip(Kes_1, solver.elementinfo.Kes))
        @test k1 ≈ k0
    end

    for _ in 1:3
        vs = [rand(T, k, k) for i in 1:N]
        f = x -> begin
            Kes = ek(x)
            sum([sum(Kes[i] * vs[i]) for i in 1:length(x)])
        end

        x = clamp.(rand(prod(nels)), 0.1, 1.0)

        val1, grad1 = NonconvexCore.value_gradient(f, x)
        val2, grad2 = f(x), Zygote.gradient(f, x)[1]
        grad3 = FDM.grad(central_fdm(5, 1), f, x)[1]
        @test val1 == val2
        @test norm(grad1 - grad2) == 0
        @test norm(grad1 - grad3) <= 1e-5
    end
end

@testset "TrussElementKσ $problem_dim" for problem_dim in ["2d", "3d"]
    file_name = "tim_$(problem_dim).json"
    problem_file = joinpath(gm_ins_dir, file_name)

    mats = TrussFEAMaterial(1.0, 0.3)
    crossecs = TrussFEACrossSec(800.0)

    node_points, elements, _, _, fixities, load_cases = load_truss_json(problem_file)
    ndim, nnodes, ncells = length(node_points[1]), length(node_points), length(elements)
    loads = load_cases["0"]

    problem = TrussProblem(
        Val{:Linear}, node_points, elements, loads, fixities, mats, crossecs
    )
    solver = FEASolver(Direct, problem)
    solver()
    u = solver.u

    esigk = TrussElementKσ(problem, solver)
    nels = length(solver.vars)
    dh = problem.ch.dh
    T = eltype(u)
    N = getncells(dh.grid)
    k = ndofs_per_cell(dh)

    # * check geometric stiffness matrix consistency
    Kσs_0 = get_truss_Kσs(problem, u, solver.elementinfo.cellvalues)

    for _ in 1:3
        vs = [rand(T, k, k) for i in 1:N]
        f = x -> begin
            Keσs = esigk(u, x)
            sum([sum(Keσs[i] * vs[i]) for i in 1:length(x)])
        end

        x = clamp.(rand(nels), 0.1, 1.0)

        Kσs_1 = esigk(u, x)
        for (ci, (k1, k0)) in enumerate(zip(Kσs_1, Kσs_0))
            @test k1 ≈ k0 * x[ci]
        end

        val1, grad1 = NonconvexCore.value_gradient(f, x)
        val2, grad2 = f(x), Zygote.gradient(f, x)[1]
        grad3 = FDM.grad(central_fdm(5, 1), f, x)[1]
        @test val1 == val2
        @test norm(grad1 - grad2) == 0
        @test norm(grad1 - grad3) <= 1e-5
    end
end

@testset "apply_boundary" begin
    nels = (2, 2)
    problem = PointLoadCantilever(Val{:Quadratic}, nels, (1.0, 1.0), 1.0, 0.3, 1.0)
    ch = problem.ch
    dh = problem.ch.dh
    T = eltype(problem.E)
    total_ndof = ndofs(dh)

    for _ in 1:3
        v = rand(T, total_ndof)
        K = sprand(Float64, total_ndof, total_ndof, 0.75)
        K = K + K'

        function f1(x)
            M = K * sum(x)
            M = apply_boundary_with_zerodiag!(M, ch)
            return sum(M * v)
        end

        function f2(x)
            M = K * sum(x)
            M = apply_boundary_with_meandiag!(M, ch)
            return sum(M * v)
        end

        x = rand(total_ndof)
        for f in [f2] #f1,
            val1, grad1 = NonconvexCore.value_gradient(f, x)
            val2, grad2 = f(x), Zygote.gradient(f, x)[1]
            grad3 = FDM.grad(central_fdm(5, 1), f, x)[1]
            # @test val1 == val2
            # @test norm(grad1 - grad2) == 0
            @test norm(grad1 - grad3) <= 1e-5
        end
    end
end

@testset "Buckling SDP constraint $problem_dim" for problem_dim in ["2d", "3d"]
    file_name = "tim_$(problem_dim).json"
    problem_file = joinpath(gm_ins_dir, file_name)

    mats = TrussFEAMaterial(1.0, 0.3)
    crossecs = TrussFEACrossSec(800.0)

    node_points, elements, _, _, fixities, load_cases = load_truss_json(problem_file)
    ndim, nnodes, ncells = length(node_points[1]), length(node_points), length(elements)
    loads = load_cases["0"]

    problem = TrussProblem(
        Val{:Linear}, node_points, elements, loads, fixities, mats, crossecs
    )

    xmin = 0.0001 # minimum density
    p = 4.0 # penalty
    V = 0.5 # maximum volume fraction
    x0 = fill(1.0, ncells) # initial design

    solver = FEASolver(Direct, problem)
    ch = problem.ch
    dh = problem.ch.dh
    total_ndof = ndofs(dh)

    comp = TopOpt.Compliance(problem, solver)
    dp = TopOpt.Displacement(solver)
    assemble_k = TopOpt.AssembleK(problem)
    element_k = ElementK(solver)
    truss_element_kσ = TrussElementKσ(problem, solver)

    # * comliance minimization objective
    obj = comp
    c = 1.0 # buckling load multiplier

    function buckling_matrix_constr(x)
        # * Array(K + c*Kσ) ⋟ 0, PSD
        # * solve for the displacement
        u = dp(x)

        # * x -> Kes, construct all the element stiffness matrices
        # a list of small matrices for each element (cell)
        Kes = element_k(x)

        # * Kes -> K (global linear stiffness matrix)
        K = assemble_k(Kes)
        K = apply_boundary_with_meandiag!(K, ch)

        # * u_e, x_e -> Ksigma_e
        Kσs = truss_element_kσ(u, x)

        # * Kσs -> Kσ
        Kσ = assemble_k(Kσs)
        Kσ = apply_boundary_with_zerodiag!(Kσ, ch)

        return Array(K + c * Kσ)
    end

    # * check initial design stability
    @test isfinite(logdet(cholesky(buckling_matrix_constr(x0))))

    for _ in 1:3
        v = rand(eltype(x0), total_ndof)
        f = x -> sum(buckling_matrix_constr(x) * v)

        x = clamp.(rand(ncells), 0.1, 1.0)

        solver.vars = x
        solver()
        K, G = buckling(problem, solver.globalinfo, solver.elementinfo, x; u=solver.u)
        @test K + G ≈ buckling_matrix_constr(x)

        val1, grad1 = NonconvexCore.value_gradient(f, x)
        val2, grad2 = f(x), Zygote.gradient(f, x)[1]
        grad3 = FDM.grad(central_fdm(5, 1), f, x)[1]

        @test val1 == val2
        @test norm(grad1 - grad2) == 0
        @test norm(grad1 - grad3) <= 1e-5
    end
end