text:
#  Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at https://mozilla.org/MPL/2.0/.
#############################################################################
# JuMP
# An algebraic modelling language for Julia
# See https://github.com/jump-dev/JuMP.jl
#############################################################################
# Hock-Schittkowski Nonlinear Test Suite
# These tests include JuMP implementations of the models described in
#  W. Hock, K. Schittkowski, Test Examples for Nonlinear Programming
#  Codes, Lecture Notes in Economics and Mathematical Systems,
#  Springer, No, 187, 1981
# More information, including original model descriptions, is available at
# http://www.klaus-schittkowski.de/downloads.htm
#############################################################################

# These tests are not executed as part of runtests because they depend on an
# NLP solver.
# TODO: Move these tests out of JuMP and into MOI so that JuMP can be tested
# separately from the solver, as it is for everything except NLP.

using Ipopt, JuMP
using Test
using MathOptInterface
const MOI = MathOptInterface

@testset "NLP solver tests" begin
    @testset "HS071" begin
        # hs071
        # Polynomial objective and constraints
        # min x1 * x4 * (x1 + x2 + x3) + x3
        # st  x1 * x2 * x3 * x4 >= 25
        #     x1^2 + x2^2 + x3^2 + x4^2 = 40
        #     1 <= x1, x2, x3, x4 <= 5
        # Start at (1,5,5,1)
        # End at (1.000..., 4.743..., 3.821..., 1.379...)
        m = Model(Ipopt.Optimizer)
        initval = [1, 5, 5, 1]
        @variable(m, 1 <= x[i = 1:4] <= 5, start = initval[i])
        @NLobjective(m, Min, x[1] * x[4] * (x[1] + x[2] + x[3]) + x[3])
        @NLconstraint(m, x[1] * x[2] * x[3] * x[4] >= 25)
        @NLconstraint(m, sum(x[i]^2 for i in 1:4) == 40)

        set_silent(m)
        JuMP.optimize!(m)

        @test JuMP.has_values(m)
        @test JuMP.termination_status(m) == MOI.LOCALLY_SOLVED
        @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT

        @test JuMP.value.(x) ≈ [1.000000, 4.742999, 3.821150, 1.379408] atol =
            1e-3
    end

    @testset "HS071 (no macros)" begin
        # hs071
        # Polynomial objective and constraints
        # min x1 * x4 * (x1 + x2 + x3) + x3
        # st  x1 * x2 * x3 * x4 >= 25
        #     x1^2 + x2^2 + x3^2 + x4^2 = 40
        #     1 <= x1, x2, x3, x4 <= 5
        # Start at (1,5,5,1)
        # End at (1.000..., 4.743..., 3.821..., 1.379...)
        m = Model(Ipopt.Optimizer)
        initval = [1, 5, 5, 1]
        @variable(m, 1 <= x[i = 1:4] <= 5, start = initval[i])
        JuMP.set_NL_objective(
            m,
            MOI.MIN_SENSE,
            :($(x[1]) * $(x[4]) * ($(x[1]) + $(x[2]) + $(x[3])) + $(x[3])),
        )
        JuMP.add_NL_constraint(
            m,
            :($(x[1]) * $(x[2]) * $(x[3]) * $(x[4]) >= 25),
        )
        JuMP.add_NL_constraint(
            m,
            :($(x[1])^2 + $(x[2])^2 + $(x[3])^2 + $(x[4])^2 == 40),
        )
        bad_expr = :(x[1]^2 + x[2]^2 + x[3]^2 + x[4]^2 == 40)
        @test_throws ErrorException JuMP.add_NL_constraint(m, bad_expr)

        set_silent(m)
        JuMP.optimize!(m)

        @test JuMP.has_values(m)
        @test JuMP.termination_status(m) == MOI.LOCALLY_SOLVED
        @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT

        @test JuMP.value.(x) ≈ [1.000000, 4.742999, 3.821150, 1.379408] atol =
            1e-3
    end

    @testset "HS071 (epigraph)" begin
        # hs071, with epigraph formulation
        # Linear objective, nonlinear constraints
        # min t
        # st  t >= x1 * x4 * (x1 + x2 + x3) + x3
        #     ...
        model = Model(Ipopt.Optimizer)
        start = [1.0, 5.0, 5.0, 1.0]
        @variable(model, 1 <= x[i = 1:4] <= 5, start = start[i])
        @variable(model, t, start = 100)
        @objective(model, Min, t)
        @NLconstraint(model, t >= x[1] * x[4] * (x[1] + x[2] + x[3]) + x[3])
        @NLconstraint(model, x[1] * x[2] * x[3] * x[4] >= 25)
        @NLconstraint(model, sum(x[i]^2 for i in 1:4) == 40)

        set_silent(model)
        JuMP.optimize!(model)
        @test JuMP.termination_status(model) == MOI.LOCALLY_SOLVED
        @test JuMP.primal_status(model) == MOI.FEASIBLE_POINT
        @test JuMP.value.(x) ≈ [1.000000, 4.742999, 3.821150, 1.379408] atol =
            1e-3
    end

    @testset "HS109" begin
        a = 50.176
        b1 = 0.25
        b = sin(b1)
        c = cos(b1)

        L = [0.0, 0.0, -0.55, -0.55, 196, 196, 196, -400, -400]
        U = [Inf, Inf, 0.55, 0.55, 252, 252, 252, 800, 800]

        m = Model(Ipopt.Optimizer)
        @variable(m, L[i] <= x[i = 1:9] <= U[i], start = 0.0)

        @NLobjective(
            m,
            Min,
            3 * x[1] + 1e-6 * x[1]^3 + 2 * x[2] + 0.522074e-6 * x[2]^3
        )

        @NLconstraint(m, x[4] - x[3] + 0.55 >= 0)
        @NLconstraint(m, x[3] - x[4] + 0.55 >= 0)
        @NLconstraint(m, 2250000 - x[1]^2 - x[8]^2 >= 0)
        @NLconstraint(m, 2250000 - x[2]^2 - x[9]^2 >= 0)
        @NLconstraint(
            m,
            x[5] * x[6] * sin(-x[3] - 0.25) +
            x[5] * x[7] * sin(-x[4] - 0.25) +
            2 * b * x[5]^2 - a * x[1] + 400 * a == 0
        )
        @NLconstraint(
            m,
            x[5] * x[6] * sin(x[3] - 0.25) +
            x[6] * x[7] * sin(x[3] - x[4] - 0.25) +
            2 * b * x[6]^2 - a * x[2] + 400 * a == 0
        )
        @NLconstraint(
            m,
            x[5] * x[7] * sin(x[4] - 0.25) +
            x[6] * x[7] * sin(x[4] - x[3] - 0.25) +
            2 * b * x[7]^2 +
            881.779 * a == 0
        )
        @NLconstraint(
            m,
            a * x[8] +
            x[5] * x[6] * cos(-x[3] - 0.25) +
            x[5] * x[7] * cos(-x[4] - 0.25) - 200 * a - 2 * c * x[5]^2 +
            0.7533e-3 * a * x[5]^2 == 0
        )
        @NLconstraint(
            m,
            a * x[9] +
            x[5] * x[6] * cos(x[3] - 0.25) +
            x[6] * x[7] * cos(x[3] - x[4] - 0.25) - 2 * c * x[6]^2 +
            0.7533e-3 * a * x[6]^2 - 200 * a == 0
        )
        @NLconstraint(
            m,
            x[5] * x[7] * cos(x[4] - 0.25) +
            x[6] * x[7] * cos(x[4] - x[3] - 0.25) - 2 * c * x[7]^2 +
            22.938 * a +
            0.7533e-3 * a * x[7]^2 == 0
        )

        set_silent(m)
        JuMP.optimize!(m)

        @test JuMP.has_values(m)
        @test JuMP.termination_status(m) == MOI.LOCALLY_SOLVED
        @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT

        @test JuMP.objective_value(m) ≈ 5326.851310161077 atol = 1e-5
    end

    @testset "HS110" begin
        m = Model(Ipopt.Optimizer)
        @variable(m, -2.001 <= x[1:10] <= 9.999, start = 9)

        @NLobjective(
            m,
            Min,
            sum(log(x[j] - 2)^2 + log(10 - x[j])^2 for j in 1:10) -
            prod(x[i] for i in 1:10)^0.2
        )

        set_silent(m)
        JuMP.optimize!(m)

        @test JuMP.has_values(m)
        # Ipopt returns AlmostLOCALLY_SOLVED and NearlyFEASIBLE_POINT on this instance.
        # @test JuMP.termination_status(m) == MOI.LOCALLY_SOLVED
        # @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT

        @test JuMP.objective_value(m) ≈ -45.77846971 atol = 1e-5
    end

    @testset "HS111" begin
        c = [
            -6.089,
            -17.164,
            -34.054,
            -5.914,
            -24.721,
            -14.986,
            -24.100,
            -10.708,
            -26.662,
            -22.179,
        ]

        m = Model(Ipopt.Optimizer)
        @variable(m, -100 <= x[1:10] <= 100, start = -2.3)

        @NLobjective(
            m,
            Min,
            sum(
                exp(x[j]) * (c[j] + x[j] - log(sum(exp(x[k]) for k in 1:10)))
                for j in 1:10
            )
        )

        @NLconstraint(
            m,
            exp(x[1]) +
            2 * exp(x[2]) +
            2 * exp(x[3]) +
            exp(x[6]) +
            exp(x[10]) == 2
        )
        @NLconstraint(m, exp(x[4]) + 2 * exp(x[5]) + exp(x[6]) + exp(x[7]) == 1)
        @NLconstraint(
            m,
            exp(x[3]) + exp(x[7]) + exp(x[8]) + 2 * exp(x[9]) + exp(x[10]) == 1
        )

        set_silent(m)
        JuMP.optimize!(m)

        @test JuMP.has_values(m)
        @test JuMP.termination_status(m) == MOI.LOCALLY_SOLVED
        @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT

        @test JuMP.objective_value(m) ≈ -47.76109026 atol = 1e-5
    end

    @testset "HS112" begin
        c = [
            -6.089,
            -17.164,
            -34.054,
            -5.914,
            -24.721,
            -14.986,
            -24.100,
            -10.708,
            -26.662,
            -22.179,
        ]

        m = Model(Ipopt.Optimizer)
        @variable(m, x[1:10] >= 1e-6, start = 0.1)

        @NLobjective(
            m,
            Min,
            sum(
                x[j] * (c[j] + log(x[j] / sum(x[k] for k in 1:10))) for
                j in 1:10
            )
        )

        @NLconstraint(m, x[1] + 2 * x[2] + 2 * x[3] + x[6] + x[10] == 2)
        @NLconstraint(m, x[4] + 2 * x[5] + x[6] + x[7] == 1)
        @NLconstraint(m, x[3] + x[7] + x[8] + 2 * x[9] + x[10] == 1)

        set_silent(m)
        JuMP.optimize!(m)

        @test JuMP.has_values(m)
        @test JuMP.termination_status(m) == MOI.LOCALLY_SOLVED
        @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT

        @test JuMP.objective_value(m) ≈ -47.76109026 atol = 1e-5
    end

    @testset "HS114" begin
        n = 10
        a = 0.99
        b = 0.9

        lower = [0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 85, 90, 3, 1.2, 145]
        upper = [2000, 16000, 120, 5000, 2000, 93, 95, 12, 4, 162]
        start = [1745, 12000, 110, 3048, 1974, 89.2, 92.8, 8, 3.6, 145]

        m = Model(Ipopt.Optimizer)
        @variable(m, lower[i] <= x[i = 1:n] <= upper[i], start = start[i])

        @NLobjective(
            m,
            Min,
            5.04 * x[1] + 0.035 * x[2] + 10 * x[3] + 3.36 * x[5] -
            0.063 * x[4] * x[7]
        )

        @NLconstraint(m, 35.82 - 0.222 * x[10] - b * x[9] >= 0)
        @NLconstraint(m, -133 + 3 * x[7] - a * x[10] >= 0)
        @NLconstraint(
            m,
            -(35.82 - 0.222 * x[10] - b * x[9]) + x[9] * (1 / b - b) >= 0
        )
        @NLconstraint(
            m,
            -(-133 + 3 * x[7] - a * x[10]) + (1 / a - a) * x[10] >= 0
        )
        @NLconstraint(
            m,
            1.12 * x[1] + 0.13167 * x[1] * x[8] - 0.00667 * x[1] * x[8]^2 -
            a * x[4] >= 0
        )
        @NLconstraint(
            m,
            57.425 + 1.098 * x[8] - 0.038 * x[8]^2 + 0.325 * x[6] - a * x[7] >=
            0
        )
        @NLconstraint(
            m,
            -(
                1.12 * x[1] + 0.13167 * x[1] * x[8] - 0.00667 * x[1] * x[8]^2 -
                a * x[4]
            ) + (1 / a - a) * x[4] >= 0
        )
        @NLconstraint(
            m,
            -(
                57.425 + 1.098 * x[8] - 0.038 * x[8]^2 + 0.325 * x[6] -
                a * x[7]
            ) + (1 / a - a) * x[7] >= 0
        )
        @NLconstraint(m, 1.22 * x[4] - x[1] - x[5] == 0)
        @NLconstraint(m, 98000 * x[3] / (x[4] * x[9] + 1000 * x[3]) - x[6] == 0)
        @NLconstraint(m, (x[2] + x[5]) / x[1] - x[8] == 0)

        set_silent(m)
        JuMP.optimize!(m)

        @test JuMP.has_values(m)
        @test JuMP.termination_status(m) == MOI.LOCALLY_SOLVED
        @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT

        @test JuMP.objective_value(m) ≈ -1768.80696 atol = 1e-3
    end

    @testset "HS116" begin
        N = 13
        a = 0.002
        b = 1.262626
        c = 1.231059
        d = 0.03475
        e = 0.975
        f = 0.00975

        lower = [
            0.1,
            0.1,
            0.1,
            0.0001,
            0.1,
            0.1,
            0.1,
            0.1,
            500,
            0.1,
            1.0,
            0.0001,
            0.0001,
            0.0,
            0.0,
            0.0,
        ]
        upper = [
            1.0,
            1.0,
            1.0,
            0.1,
            0.9,
            0.9,
            1000,
            1000,
            1000,
            500,
            150,
            150,
            150,
            Inf,
            Inf,
            Inf,
        ]
        start =
            [0.5 2 0.8 3 0.9 4 0.1 5 0.14 6 0.5 7 489 8 80 9 650 0.5 2 0.8 3 0.9 4 0.1 5 0.14 6 0.5 7 489 8 80 9 650]

        m = Model(Ipopt.Optimizer)
        @variable(m, lower[i] <= x[i = 1:N] <= upper[i], start = start[i])
        @NLobjective(m, Min, x[11] + x[12] + x[13])

        @NLconstraints m begin
            x[3] - x[2] >= 0
            x[2] - x[1] >= 0
            1 - a * x[7] + a * x[8] >= 0
            x[11] + x[12] + x[13] >= 50
            x[13] - b * x[10] + c * x[3] * x[10] >= 0
            x[5] - d * x[2] - e * x[2] * x[5] + f * x[2]^2 >= 0
            x[6] - d * x[3] - e * x[3] * x[6] + f * x[3]^2 >= 0
            x[4] - d * x[1] - e * x[1] * x[4] + f * x[1]^2 >= 0
            x[12] - b * x[9] + c * x[2] * x[9] >= 0
            x[11] - b * x[8] + c * x[1] * x[8] >= 0
            x[5] * x[7] - x[1] * x[8] - x[4] * x[7] + x[4] * x[8] >= 0
            1 - a * (x[2] * x[9] + x[5] * x[8] - x[1] * x[8] - x[6] * x[9]) -
            x[5] - x[6] >= 0
            x[2] * x[9] - x[3] * x[10] - x[6] * x[9] - 500 * x[2] +
            500 * x[6] +
            x[2] * x[10] >= 0
            x[2] - 0.9 - a * (x[2] * x[10] - x[3] * x[10]) >= 0
            x[11] + x[12] + x[13] <= 250
        end

        set_silent(m)
        JuMP.optimize!(m)

        @test JuMP.has_values(m)
        @test JuMP.termination_status(m) == MOI.LOCALLY_SOLVED
        @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT

        # This test occasionally fails, for unknown reasons.
        @test JuMP.objective_value(m) ≈ 97.588409 atol = 1e-3
    end

    @testset "HS118" begin
        m = Model(Ipopt.Optimizer)

        L = zeros(15)
        L[1] = 8.0
        L[2] = 43.0
        L[3] = 3.0

        U = zeros(15)
        U[1] = 21.0
        U[2] = 57.0
        U[3] = 16.0
        for k in 1:4
            U[3*k+1] = 90.0
            U[3*k+2] = 120.0
            U[3*k+3] = 60.0
        end

        start = [
            20.0,
            55.0,
            15.0,
            20.0,
            60.0,
            20.0,
            20.0,
            60.0,
            20.0,
            20.0,
            60.0,
            20.0,
            20.0,
            60.0,
            20.0,
        ]

        @variable(m, L[i] <= x[i = 1:15] <= U[i])

        # Initial solution (could also use 'start' keyword in @variable)
        JuMP.set_start_value.(x, start)

        @NLobjective(
            m,
            Min,
            sum(
                2.3 * x[3*k+1] +
                0.0001 * x[3*k+1]^2 +
                1.7 * x[3*k+2] +
                0.0001 * x[3*k+2]^2 +
                2.2 * x[3*k+3] +
                0.00015 * x[3*k+3]^2 for k in 0:4
            )
        )

        # constr1
        for j in 1:4
            @constraint(m, x[3*j+1] - x[3*j-2] + 7 <= 13)
            @constraint(m, x[3*j+1] - x[3*j-2] + 7 >= 0)
        end

        # constr2
        for j in 1:4
            @constraint(m, x[3*j+2] - x[3*j-1] + 7 <= 14)
            @constraint(m, x[3*j+2] - x[3*j-1] + 7 >= 0)
        end

        # constr3
        for j in 1:4
            @constraint(m, x[3*j+3] - x[3*j] + 7 <= 13)
            @constraint(m, x[3*j+3] - x[3*j] + 7 >= 0)
        end

        @constraint(m, x[1] + x[2] + x[3] >= 60)
        @constraint(m, x[4] + x[5] + x[6] >= 50)
        @constraint(m, x[7] + x[8] + x[9] >= 70)
        @constraint(m, x[10] + x[11] + x[12] >= 85)
        @constraint(m, x[13] + x[14] + x[15] >= 100)

        set_silent(m)
        JuMP.optimize!(m)

        @test JuMP.has_values(m)
        @test JuMP.termination_status(m) == MOI.LOCALLY_SOLVED
        @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT

        @test JuMP.value.(x[1:4]) ≈ [8.0, 49.0, 3.0, 1.0] atol = 1e-4
        @test JuMP.objective_value(m) ≈ 664.82045 atol = 1e-5
    end

    @testset "Two-sided constraints" begin
        m = Model(Ipopt.Optimizer)
        @variable(m, x)
        @NLobjective(m, Max, x)
        l = -1
        u = 1
        @NLconstraint(m, l <= x <= u)
        JuMP.optimize!(m)

        @test JuMP.has_values(m)
        @test JuMP.termination_status(m) == MOI.LOCALLY_SOLVED
        @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
        @test JuMP.objective_value(m) ≈ u atol = 1e-6

        @NLobjective(m, Min, x)
        set_silent(m)
        JuMP.optimize!(m)

        @test JuMP.has_values(m)
        @test JuMP.termination_status(m) == MOI.LOCALLY_SOLVED
        @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
        @test JuMP.objective_value(m) ≈ l atol = 1e-6
    end

    @testset "Two-sided constraints (no macros)" begin
        m = Model(Ipopt.Optimizer)
        @variable(m, x)
        JuMP.set_NL_objective(m, MOI.MAX_SENSE, x)
        l = -1
        u = 1
        JuMP.add_NL_constraint(m, :($l <= $x <= $u))
        set_silent(m)
        JuMP.optimize!(m)

        @test JuMP.has_values(m)
        @test JuMP.termination_status(m) == MOI.LOCALLY_SOLVED
        @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
        @test JuMP.objective_value(m) ≈ u atol = 1e-6

        JuMP.set_NL_objective(m, MOI.MIN_SENSE, x)
        JuMP.optimize!(m)

        @test JuMP.has_values(m)
        @test JuMP.termination_status(m) == MOI.LOCALLY_SOLVED
        @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
        @test JuMP.objective_value(m) ≈ l atol = 1e-6
    end

    @testset "Duals" begin
        m = Model(Ipopt.Optimizer)
        @variable(m, x >= 0)
        @variable(m, y <= 5)
        @variable(m, 2 <= z <= 4)
        @variable(m, 0 <= r[i = 3:6] <= i)
        @NLobjective(m, Min, -((x + y) / 2.0 + 3.0) / 3.0 - z - r[3])
        @constraint(m, cons1, x + y >= 2)
        @constraint(m, cons2, sum(r[i] for i in 3:5) <= (2 - x) / 2.0)
        @NLconstraint(m, cons3, 7.0 * y <= z + r[6] / 1.9)

        function test_result()
            @test JuMP.has_values(m)
            @test JuMP.termination_status(m) == MOI.LOCALLY_SOLVED
            @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT

            @test JuMP.value(x) ≈ 0.9774436 atol = 1e-6
            @test JuMP.value(y) ≈ 1.0225563 atol = 1e-6
            @test JuMP.value(z) ≈ 4.0 atol = 1e-6
            @test JuMP.value(r[3]) ≈ 0.5112781 atol = 1e-6
            @test JuMP.value(r[4]) ≈ 0.0 atol = 1e-6
            @test JuMP.value(r[5]) ≈ 0.0 atol = 1e-6
            @test JuMP.value(r[6]) ≈ 6.0 atol = 1e-6
            @test JuMP.dual_status(m) == MOI.FEASIBLE_POINT
            # Reduced costs
            @test JuMP.dual(JuMP.LowerBoundRef(x)) ≈ 0.0 atol = 1e-6
            @test JuMP.dual(JuMP.UpperBoundRef(y)) ≈ 0.0 atol = 1e-6
            @test JuMP.dual(JuMP.UpperBoundRef(z)) ≈ -1.0714286 atol = 1e-6
            @test JuMP.dual(JuMP.LowerBoundRef(r[3])) ≈ 0.0 atol = 1e-6
            @test JuMP.dual(JuMP.UpperBoundRef(r[3])) ≈ 0.0 atol = 1e-6
            @test JuMP.dual(JuMP.LowerBoundRef(r[4])) ≈ 1.0 atol = 1e-6
            @test JuMP.dual(JuMP.UpperBoundRef(r[4])) ≈ 0.0 atol = 1e-6
            @test JuMP.dual(JuMP.LowerBoundRef(r[5])) ≈ 1.0 atol = 1e-6
            @test JuMP.dual(JuMP.UpperBoundRef(r[5])) ≈ 0.0 atol = 1e-6
            @test JuMP.dual(JuMP.UpperBoundRef(r[6])) ≈ -0.03759398 atol = 1e-6
            @test JuMP.dual(JuMP.LowerBoundRef(r[6])) ≈ 0.0 atol = 1e-6

            # Constraint duals
            @test JuMP.dual(cons1) ≈ 0.333333 atol = 1e-6
            @test JuMP.dual(cons2) ≈ -1.0 atol = 1e-6
            @test JuMP.dual(cons3) ≈ -0.0714286 atol = 1e-6
        end

        set_silent(m)
        JuMP.optimize!(m)
        test_result()
        @test JuMP.objective_value(m) ≈ -5.8446115 atol = 1e-6

        # Same objective with sense/sign flipped.
        @NLobjective(m, Max, ((x + y) / 2.0 + 3.0) / 3.0 + z + r[3])

        JuMP.optimize!(m)
        test_result()
        @test JuMP.objective_value(m) ≈ 5.8446115 atol = 1e-6
    end

    @testset "Quadratic inequality constraints, linear objective" begin
        m = Model(Ipopt.Optimizer)
        @variable(m, -2 <= x <= 2)
        @variable(m, -2 <= y <= 2)
        @objective(m, Min, x - y)
        @constraint(m, x + x^2 + x * y + y^2 <= 1)
        set_silent(m)
        JuMP.optimize!(m)

        @test JuMP.has_values(m)
        @test JuMP.termination_status(m) == MOI.LOCALLY_SOLVED
        @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
        @test JuMP.objective_value(m) ≈ -1 - 4 / sqrt(3) atol = 1e-6
        @test JuMP.value(x) + JuMP.value(y) ≈ -1 / 3 atol = 1e-3
    end

    @testset "Quadratic inequality constraints, NL objective" begin
        m = Model(Ipopt.Optimizer)
        @variable(m, -2 <= x <= 2)
        @variable(m, -2 <= y <= 2)
        @NLobjective(m, Min, x - y)
        @constraint(m, x + x^2 + x * y + y^2 <= 1)
        set_silent(m)
        JuMP.optimize!(m)

        @test JuMP.has_values(m)
        @test JuMP.termination_status(m) == MOI.LOCALLY_SOLVED
        @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
        @test JuMP.objective_value(m) ≈ -1 - 4 / sqrt(3) atol = 1e-6
        @test JuMP.value(x) + JuMP.value(y) ≈ -1 / 3 atol = 1e-3
    end

    @testset "Quadratic equality constraints" begin
        m = Model(Ipopt.Optimizer)
        @variable(m, 0 <= x[1:2] <= 1)
        @constraint(m, x[1]^2 + x[2]^2 == 1 / 2)
        @NLobjective(m, Max, x[1] - x[2])
        set_silent(m)
        JuMP.optimize!(m)

        @test JuMP.has_values(m)
        @test JuMP.termination_status(m) == MOI.LOCALLY_SOLVED
        @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
        @test JuMP.objective_value(m) ≈ sqrt(1 / 2) atol = 1e-6
        @test JuMP.value.(x) ≈ [sqrt(1 / 2), 0] atol = 1e-6
    end

    @testset "Fixed variables" begin
        m = Model(Ipopt.Optimizer)
        set_silent(m)
        @variable(m, x == 0)
        @variable(m, y ≥ 0)
        @objective(m, Min, y)
        @NLconstraint(m, y ≥ x^2)
        for α in 1:4
            JuMP.fix(x, α)
            JuMP.optimize!(m)
            @test JuMP.value(y) ≈ α^2 atol = 1e-6
        end
    end

    @testset "ifelse" begin
        model = Model(Ipopt.Optimizer)
        @variable(model, x, start = 2)
        # The minimizer is at smooth point, so solvers should be okay.
        @NLobjective(model, Min, ifelse(x <= 1, x^2, x))
        set_silent(model)
        JuMP.optimize!(model)

        @test JuMP.termination_status(model) == MOI.LOCALLY_SOLVED
        @test JuMP.primal_status(model) == MOI.FEASIBLE_POINT
        @test JuMP.value(x) ≈ 0.0 atol = 1e-5
    end

    @testset "infeasible problem" begin
        # (Attempt to) solve an infeasible problem
        model = Model(Ipopt.Optimizer)
        n = 10
        @variable(model, 0 <= x[i = 1:n] <= 1)
        @NLobjective(model, Max, x[n])
        for i in 1:n-1
            @NLconstraint(model, x[i+1] - x[i] == 0.15)
        end
        set_silent(model)
        JuMP.optimize!(model)

        @test JuMP.termination_status(model) == MOI.LOCALLY_INFEASIBLE
    end

    @testset "unbounded problem" begin
        model = Model(Ipopt.Optimizer)
        @variable(model, x >= 0)
        @NLobjective(model, Max, x)
        @NLconstraint(model, x >= 5)
        set_silent(model)
        JuMP.optimize!(model)
        @test JuMP.termination_status(model) == MOI.NORM_LIMIT
    end

    @testset "Derivatives of x^4, x < 0" begin
        model = Model(Ipopt.Optimizer)
        @variable(model, x >= -1, start = -0.5)
        @NLobjective(model, Min, x^4)
        set_silent(model)
        JuMP.optimize!(model)

        @test JuMP.termination_status(model) == MOI.LOCALLY_SOLVED
        @test JuMP.primal_status(model) == MOI.FEASIBLE_POINT
        @test JuMP.value(x) ≈ 0.0 atol = 1e-2
    end

    # This test seems to be checking for a bug in expression handling.
    @testset "Entropy maximization" begin
        model = Model(Ipopt.Optimizer)
        @variable(model, x[1:4] >= 0, start = 1)
        @variable(model, z[1:4], start = 0)
        @NLexpression(model, entropy[i = 1:4], -x[i] * log(x[i]))
        @NLobjective(
            model,
            Max,
            sum(z[i] for i in 1:2) + sum(z[i] / 2 for i in 3:4)
        )
        @NLconstraint(model, z_constr1[i = 1], z[i] <= entropy[i])
        @NLconstraint(model, z_constr1_dup[i = 2], z[i] <= entropy[i]) # duplicate expressions
        @NLconstraint(model, z_constr2[i = 3:4], z[i] <= 2 * entropy[i])
        @constraint(model, sum(x) == 1)
        set_silent(model)
        JuMP.optimize!(model)

        @test JuMP.termination_status(model) == MOI.LOCALLY_SOLVED
        @test JuMP.primal_status(model) == MOI.FEASIBLE_POINT
        @test JuMP.value.(x) ≈ [1 / 4, 1 / 4, 1 / 4, 1 / 4] atol = 1e-3
        @test JuMP.value(entropy[1]) ≈ -(1 / 4) * log(1 / 4) atol = 1e-4
    end

    @testset "Changing objectives" begin
        model = Model(Ipopt.Optimizer)
        @variable(model, x >= 0)
        @variable(model, y >= 0)
        @objective(model, Max, x + y)
        @NLconstraint(model, x + 2y <= 1)
        JuMP.optimize!(model)
        @test JuMP.value(x) ≈ 1.0 atol = 1e-4
        @test JuMP.value(y) ≈ 0.0 atol = 1e-4
        @test JuMP.objective_value(model) ≈ 1.0 atol = 1e-4

        @objective(model, Max, 2x + y)
        set_silent(model)
        JuMP.optimize!(model)
        @test JuMP.value(x) ≈ 1.0 atol = 1e-4
        @test JuMP.value(y) ≈ 0.0 atol = 1e-4
        @test JuMP.objective_value(model) ≈ 2.0 atol = 1e-4
    end

    @testset "Setting NLobjective then objective" begin
        model = Model(Ipopt.Optimizer)
        @variable(model, x >= 0)
        @variable(model, y >= 0)
        @NLobjective(model, Max, x + y)
        @NLconstraint(model, x + 2y <= 1)
        JuMP.optimize!(model)
        @test JuMP.value(x) ≈ 1.0 atol = 1e-4
        @test JuMP.value(y) ≈ 0.0 atol = 1e-4
        @test JuMP.objective_value(model) ≈ 1.0 atol = 1e-4

        @objective(model, Max, 2x + y)
        set_silent(model)
        JuMP.optimize!(model)
        @test JuMP.value(x) ≈ 1.0 atol = 1e-4
        @test JuMP.value(y) ≈ 0.0 atol = 1e-4
        @test JuMP.objective_value(model) ≈ 2.0 atol = 1e-4
    end

    my_square(x) = x^2
    function my_f(x, y)
        return (x - 1)^2 + (y - 2)^2
    end

    @testset "User-defined functions" begin
        model = Model(Ipopt.Optimizer)
        JuMP.register(model, :my_f, 2, my_f, autodiff = true)
        JuMP.register(model, :my_square, 1, my_square, autodiff = true)

        @variable(model, x[1:2] >= 0.5)
        @NLobjective(model, Min, my_f(x[1], my_square(x[2])))

        set_silent(model)
        JuMP.optimize!(model)
        @test JuMP.termination_status(model) == MOI.LOCALLY_SOLVED
        @test JuMP.primal_status(model) == MOI.FEASIBLE_POINT
        @test JuMP.value.(x) ≈ [1, sqrt(2.0)]
    end

    @testset "Univariate user-defined functions" begin
        model = Model(Ipopt.Optimizer)
        JuMP.register(model, :my_square, 1, my_square, autodiff = true)

        # Test just univariate functions because this is a path where hessians
        # are enabled.
        @variable(model, x[1:2] >= 0.5)
        @NLobjective(model, Min, my_square(x[1] - 0.4) + my_square(x[2] - 2.0))
        set_silent(model)
        JuMP.optimize!(model)
        @test JuMP.termination_status(model) == MOI.LOCALLY_SOLVED
        @test JuMP.primal_status(model) == MOI.FEASIBLE_POINT
        @test JuMP.value.(x) ≈ [0.5, 2.0] atol = 1e-4
    end

    @testset "Issue #927" begin
        model = Model(Ipopt.Optimizer)
        JuMP.register(model, :my_f, 2, my_f, autodiff = true)
        @variable(model, x)
        @NLobjective(model, Min, my_f(x, x))
        set_silent(model)
        JuMP.optimize!(model)
        @test JuMP.termination_status(model) == MOI.LOCALLY_SOLVED
        @test JuMP.primal_status(model) == MOI.FEASIBLE_POINT
        @test JuMP.value.(x) ≈ 1.5 atol = 1e-4
    end
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
const USE_GPU = false  # Use GPU? If this is set false, then no GPU needs to be available
using ParallelStencil
using ParallelStencil.FiniteDifferences2D
@static if USE_GPU
    @init_parallel_stencil(CUDA, Float64, 2)
else
    @init_parallel_stencil(Threads, Float64, 2)
end
using Plots, Printf, Statistics

@parallel function compute_V!(Vx::Data.Array, Vy::Data.Array, P::Data.Array, dt::Data.Number, ρ::Data.Number, dx::Data.Number, dy::Data.Number)
    @inn(Vx) = @inn(Vx) - dt/ρ*@d_xi(P)/dx
    @inn(Vy) = @inn(Vy) - dt/ρ*@d_yi(P)/dy
    return
end

@parallel function compute_P!(P::Data.Array, Vx::Data.Array, Vy::Data.Array, dt::Data.Number, k::Data.Number, dx::Data.Number, dy::Data.Number)
    @all(P) = @all(P) - dt*k*(@d_xa(Vx)/dx + @d_ya(Vy)/dy)
    return
end

##################################################
@views function acoustic2D()
    # Physics
    lx, ly    = 40.0, 40.0  # domain extends
    k         = 1.0         # bulk modulus
    ρ         = 1.0         # density
    t         = 0.0         # physical time
    # Numerics
    nx, ny    = 255, 255    # numerical grid resolution; should be a mulitple of 32-1 for optimal GPU perf
    nt        = 1000        # number of timesteps
    nout      = 10          # plotting frequency
    # Derived numerics
    dx, dy    = lx/(nx-1), ly/(ny-1) # cell sizes
    # Array allocations
    P         = @zeros(nx  ,ny  )
    Vx        = @zeros(nx+1,ny  )
    Vy        = @zeros(nx  ,ny+1)
    # Initial conditions
    P        .= Data.Array([exp(-((ix-1)*dx-0.5*lx)^2 -((iy-1)*dy-0.5*ly)^2) for ix=1:size(P,1), iy=1:size(P,2)])
    dt        = min(dx,dy)/sqrt(k/ρ)/4.1
    # Preparation of visualisation
    ENV["GKSwstype"]="nul"; if isdir("viz2D_out")==false mkdir("viz2D_out") end; loadpath = "./viz2D_out/"; anim = Animation(loadpath,String[])
    println("Animation directory: $(anim.dir)")
    X, Y      = -lx/2:dx:lx/2, -ly/2:dy:ly/2
    # Time loop
    for it = 1:nt
        if (it==11)  global wtime0 = Base.time()  end
        @parallel compute_V!(Vx, Vy, P, dt, ρ, dx, dy)
        @parallel compute_P!(P, Vx, Vy, dt, k, dx, dy)
        t = t + dt
        # Visualisation
        if mod(it,nout)==0
            heatmap(X, Y, Array(P)', aspect_ratio=1, xlims=(X[1],X[end]), ylims=(Y[1],Y[end]), c=:viridis, title="Pressure"); frame(anim)
        end
    end
    # Performance
    wtime    = Base.time()-wtime0
    A_eff    = (3*2)/1e9*nx*ny*sizeof(Data.Number)  # Effective main memory access per iteration [GB] (Lower bound of required memory access: H and dHdτ have to be read and written (dHdτ for damping): 4 whole-array memaccess; B has to be read: 1 whole-array memaccess)
    wtime_it = wtime/(nt-10)                        # Execution time per iteration [s]
    T_eff    = A_eff/wtime_it                       # Effective memory throughput [GB/s]
    @printf("Total steps=%d, time=%1.3e sec (@ T_eff = %1.2f GB/s) \n", nt, wtime, round(T_eff, sigdigits=2))
    gif(anim, "acoustic2D.gif", fps = 15)
    return
end

acoustic2D()

ext: jl


text:
using Test
using FlightMechanics

@testset "ISA1978" begin
    # Test sea level
    T, p, rho, a = atmosphere_isa(0.0)
    @test isapprox(T, 288.15)
    @test isapprox(p, 101325)
    @test isapprox(rho, 1.225)
    @test isapprox(a, 340.29398, rtol=1e-7)

    # Test 0-11 Km
    h = [0.0, 50.0, 550.0, 6500.0, 10000.0, 11000.0]
    # TODO: have a look again at unzip options
    # I would expect a tuple of Arrays, however an Array of tuples is returned.
    # https://discourse.julialang.org/t/broadcast-map-return-type-for-a-function-returning-a-tuple/574/6
    rv = atmosphere_isa.(h)
    rv_arr = reduce(hcat, rv)
    T = rv_arr[1, :]
    p = rv_arr[2, :]
    rho = rv_arr[3, :]
    a = rv_arr[4, :]
    @test isapprox(T, [288.150, 287.825, 284.575, 245.900, 223.150, 216.650])
    @test isapprox(rho, [1.2250, 1.2191, 1.1616, 0.62384, 0.41271, 0.36392], rtol=1e-4)
    @test isapprox(a, [340.29, 340.10, 338.18, 314.36, 299.46, 295.07], rtol=1e-5)

    # Test 11-20 Km
    h = [12000, 14200, 17500, 20000, ]  # m
    rv = atmosphere_isa.(h)
    rv_arr = reduce(hcat, rv)
    T = rv_arr[1, :]
    p = rv_arr[2, :]
    rho = rv_arr[3, :]
    a = rv_arr[4, :]
    @test isapprox(T, [216.650, 216.650, 216.650, 216.650])
    @test isapprox(rho, [0.31083, 0.21971, 0.13058, 0.088035], rtol=1e-4)
    @test isapprox(a, [295.07, 295.07, 295.07, 295.07], rtol=1e-5)

    # Test 20-32 Km
    h = [22100, 24000, 28800, 32000]  # m
    rv = atmosphere_isa.(h)
    rv_arr = reduce(hcat, rv)
    T = rv_arr[1, :]
    p = rv_arr[2, :]
    rho = rv_arr[3, :]
    a = rv_arr[4, :]
    @test isapprox(T, [218.750, 220.650, 225.450, 228.650])
    @test isapprox(rho, [0.062711, 0.046267, 0.021708, 0.013225], rtol=1e-4)
    @test isapprox(a, [296.50, 297.78, 301.00, 303.13], rtol=1e-5)

    # Test 32-47 Km
    h = [32200, 36000, 42000, 47000]  # m
    rv = atmosphere_isa.(h)
    rv_arr = reduce(hcat, rv)
    T = rv_arr[1, :]
    p = rv_arr[2, :]
    rho = rv_arr[3, :]
    a = rv_arr[4, :]
    @test isapprox(T, [229.210, 239.850, 256.650, 270.650])
    @test isapprox(rho, [0.012805, 0.0070344, 0.0028780, 0.0014275], rtol=1e-4)
    @test isapprox(a, [303.50, 310.47, 321.16, 329.80], rtol=1e-5)

    # Test 47-51 Km
    h = [47200, 49000, 51000]  # m
    rv = atmosphere_isa.(h)
    rv_arr = reduce(hcat, rv)
    T = rv_arr[1, :]
    p = rv_arr[2, :]
    rho = rv_arr[3, :]
    a = rv_arr[4, :]
    @test isapprox(T, [270.650, 270.650, 270.650])
    @test isapprox(rho, [0.0013919, 0.0011090, 0.00086160], rtol=1e-4)
    @test isapprox(a, [329.80, 329.80, 329.80], rtol=1e-5)

    # Test 51-71 Km
    h = [51500, 60000, 71000]  # m
    rv = atmosphere_isa.(h)
    rv_arr = reduce(hcat, rv)
    T = rv_arr[1, :]
    p = rv_arr[2, :]
    rho = rv_arr[3, :]
    a = rv_arr[4, :]
    @test isapprox(T, [269.250, 245.450, 214.650])
    @test isapprox(rho, [0.00081298, 2.8832e-4, 6.4211e-5], rtol=1e-4)
    @test isapprox(a, [328.94, 314.07, 293.70], rtol=1e-4)

    # Test 71-84 Km
    h = [52000, 60000, 84500]  # m
    rv = atmosphere_isa.(h)
    rv_arr = reduce(hcat, rv)
    T = rv_arr[1, :]
    p = rv_arr[2, :]
    rho = rv_arr[3, :]
    a = rv_arr[4, :]
    @test isapprox(T, [267.850, 245.450, 187.650])
    @test isapprox(rho, [7.6687e-4, 2.8832e-4, 7.3914e-6], rtol=1e-4)
    @test isapprox(a, [328.09, 314.07, 274.61], rtol=1e-5)

    @test_throws DomainError atmosphere_isa(84500.1)
    @test_throws DomainError atmosphere_isa(-0.1)
end

@testset "AtmosphereF16" begin
    import FlightMechanics: atmosphere_f16
    # Test sea level
    T, p, rho, a = atmosphere_f16(0.0)

    # Passing for ISA1978
    @test_broken isapprox(T, 288.15)
    @test_broken isapprox(p, 101325)
    @test_broken isapprox(rho, 1.225)
    @test_broken isapprox(a, 340.29398, rtol=1e-7)

    @test isapprox(T, 288.15, atol=0.5)
    @test isapprox(p, 101325, atol=30)
    @test isapprox(rho, 1.225, atol=1e-4)
    @test isapprox(a, 340.29398, atol=0.1)

    # Test 0-11 Km
    h = [0.0, 50.0, 550.0, 6500.0, 10000.0, 11000.0]
    # TODO: have a look again at unzip options
    # I would expect a tuple of Arrays, however an Array of tuples is returned.
    # https://discourse.julialang.org/t/broadcast-map-return-type-for-a-function-returning-a-tuple/574/6
    rv = atmosphere_f16.(h)
    rv_arr = reduce(hcat, rv)
    T = rv_arr[1, :]
    p = rv_arr[2, :]
    rho = rv_arr[3, :]
    a = rv_arr[4, :]

    # Passing for ISA1978
    @test_broken isapprox(T, [288.150, 287.825, 284.575, 245.900, 223.150, 216.650])
    @test_broken isapprox(rho, [1.2250, 1.2191, 1.1616, 0.62384, 0.41271, 0.36392], rtol=1e-4)
    @test_broken isapprox(a, [340.29, 340.10, 338.18, 314.36, 299.46, 295.07], rtol=1e-5)

    @test isapprox(T, [288.150, 287.825, 284.575, 245.900, 223.150, 216.650], rtol=0.1)
    @test isapprox(rho, [1.2250, 1.2191, 1.1616, 0.62384, 0.41271, 0.36392], rtol=1e-2)
    @test isapprox(a, [340.29, 340.10, 338.18, 314.36, 299.46, 295.07], rtol=1e-2)

    # Test 11-20 Km
    h = [12000, 14200, 17500, 20000, ]  # m
    rv = atmosphere_f16.(h)
    rv_arr = reduce(hcat, rv)
    T = rv_arr[1, :]
    p = rv_arr[2, :]
    rho = rv_arr[3, :]
    a = rv_arr[4, :]

    # Passing for ISA1978
    @test_broken isapprox(T, [216.650, 216.650, 216.650, 216.650])
    @test_broken isapprox(rho, [0.31083, 0.21971, 0.13058, 0.088035], rtol=1e-4)
    @test_broken isapprox(a, [295.07, 295.07, 295.07, 295.07], rtol=1e-5)

    @test isapprox(T, [216.650, 216.650, 216.650, 216.650], rtol=0.1)
    @test isapprox(rho, [0.31083, 0.21971, 0.13058, 0.088035], rtol=6*1e-2)
    @test isapprox(a, [295.07, 295.07, 295.07, 295.07], rtol=1e-2)
end


text:
function buildN(A::fPEPS,
                L::Environments,
                R::Environments,
                IEnvs,
                row::Int,
                col::Int,
                ϕ::ITensor)::ITensor
    Ny, Nx   = size(A)
    N        = spinI(firstind(A[row, col], "Site"); is_gpu=is_gpu(A))
    workingN = N
    if row > 1
        workingN *= IEnvs[:below][row - 1]
    end
    workingN *= L.I[row]
    workingN *= ϕ
    workingN *= R.I[row]
    if row < Ny
        workingN *= IEnvs[:above][end - row]
    end
    return workingN
end

function buildHIedge(A::fPEPS,
                     E::Environments,
                     row::Int,
                     col::Int,
                     side::Symbol,
                     ϕ::ITensor )
    Ny, Nx = size(A)
    is_cu  = is_gpu(A)
    IH_a   = is_cu ? cuITensor(1.0) : ITensor(1.0)
    IH_b   = is_cu ? cuITensor(1.0) : ITensor(1.0)
    IH     = is_cu ? cuITensor(1.0) : ITensor(1.0)
    next_col = side == :left ? 2 : Nx - 1
    @inbounds for work_row in 1:row-1
        IH_b *= E.H[work_row] * A[work_row, col]
        IH_b *= dag(prime(A[work_row, col], "Link"))
    end
    @inbounds for work_row in row+1:Ny
        IH_a *= E.H[work_row] * A[work_row, col]
        IH_a *= dag(prime(A[work_row, col], "Link"))
    end
    IH *= E.H[row]
    IH *= IH_b
    op  = spinI(firstind(A[row, col], "Site"); is_gpu=is_cu)
    op  = is_cu ? cuITensor(op) : op
    IH *= ϕ
    IH *= op
    IH *= IH_a
    AAinds = inds(prime(ϕ))
    @assert hasinds(inds(IH), AAinds)
    @assert hasinds(AAinds, inds(IH))
    return (IH,)
end

function buildHIs(A::fPEPS,
                  L::Environments,
                  R::Environments,
                  row::Int,
                  col::Int,
                  ϕ::ITensor)
    Ny, Nx = size(A)
    is_cu  = is_gpu(A)
    @timeit "buildHIedge" begin
        if col == 1
            return buildHIedge(A, R, row, col, :left, ϕ)
        end
        if col == Nx
            return buildHIedge(A, L, row, col, :right, ϕ)
        end
    end
    HLI_a = is_cu ? cuITensor(1.0) : ITensor(1.0)
    HLI_b = is_cu ? cuITensor(1.0) : ITensor(1.0)
    IHR_a = is_cu ? cuITensor(1.0) : ITensor(1.0)
    IHR_b = is_cu ? cuITensor(1.0) : ITensor(1.0)
    HLI   = is_cu ? cuITensor(1.0) : ITensor(1.0)
    IHR   = is_cu ? cuITensor(1.0) : ITensor(1.0)
    @timeit "build HI lower" begin
        @inbounds for work_row in 1:row-1
            tL     = A[work_row, col]*L.H[work_row]
            HLI_b *= tL
            HLI_b *= R.I[work_row]
            HLI_b *= dag(prime(A[work_row, col], "Link"))
            tR     = A[work_row, col]*R.H[work_row]
            IHR_b *= tR
            IHR_b *= L.I[work_row]
            IHR_b *= dag(prime(A[work_row, col], "Link"))
        end
    end
    @timeit "build HI upper" begin
        @inbounds for work_row in reverse(row+1:Ny)
            tL     = A[work_row, col]*L.H[work_row]
            HLI_a *= tL
            HLI_a *= R.I[work_row]
            HLI_a *= dag(prime(A[work_row, col], "Link"))
            tR     = A[work_row, col]*R.H[work_row]
            IHR_a *= tR
            IHR_a *= L.I[work_row]
            IHR_a *= dag(prime(A[work_row, col], "Link"))
        end
    end
    @timeit "join up HIs" begin
        HLI  *= L.H[row]
        HLI  *= HLI_b
        HLI  *= ϕ
        HLI  *= R.I[row]
        IHR  *= L.I[row]
        IHR  *= IHR_b
        IHR  *= ϕ
        IHR  *= R.H[row]
        HLI *= HLI_a
        IHR *= IHR_a
        op   = spinI(firstind(A[row, col], "Site"); is_gpu=is_cu)
        HLI *= op
        IHR *= op
    end
    AAinds = inds(prime(ϕ))
    @assert hasinds(inds(IHR), AAinds)
    @assert hasinds(inds(HLI), AAinds)
    @assert hasinds(AAinds, inds(IHR))
    @assert hasinds(AAinds, inds(HLI))
    return (HLI, IHR)
end

function verticalTerms(A::fPEPS,
                       L::Environments,
                       R::Environments,
                       AI,
                       AV,
                       H,
                       row::Int,
                       col::Int,
                       ϕ::ITensor)::Vector{ITensor}
    Ny, Nx = size(A)
    is_cu  = is_gpu(A)
    vTerms = ITensor[]
    AAinds = inds(prime(ϕ))
    dummy  = is_cu ? cuITensor(1.0) : ITensor(1.0)
    @inbounds for opcode in 1:length(H)
        thisVert = dummy
        op_row_a = H[opcode].sites[1][1]
        op_row_b = H[opcode].sites[2][1]
        if op_row_b < row || row < op_row_a
            local V, I
            if op_row_a > row
                V = AV[:above][opcode][end - row]
                I = row > 1 ? AI[:below][row - 1] : dummy
            elseif op_row_b < row
                V = AV[:below][opcode][row - op_row_b]
                I = row < Ny ? AI[:above][end - row] : dummy
            end
            thisVert *= V
            thisVert *= L.I[row]
            thisVert *= ϕ
            thisVert *= R.I[row]
            thisVert *= I
            thisVert *= spinI(firstind(A[row, col], "Site"); is_gpu=is_cu)
        elseif row == op_row_a
            low_row  = op_row_a - 1
            high_row = op_row_b
            AIL = low_row > 0 ? AI[:below][low_row] : dummy
            AIH = high_row < Ny ? AI[:above][end - high_row] : dummy
            sA   = firstind(A[op_row_a, col], "Site")
            op_a = replaceind!(copy(H[opcode].ops[1]), H[opcode].site_ind, sA)
            op_a = replaceind!(op_a, H[opcode].site_ind', sA')
            sB   = firstind(A[op_row_b, col], "Site")
            op_b = replaceind!(copy(H[opcode].ops[2]), H[opcode].site_ind, sB)
            op_b = replaceind!(op_b, H[opcode].site_ind', sB')
            thisVert  = AIH
            thisVert *= L.I[op_row_b]
            thisVert *= A[op_row_b, col] * op_b
            thisVert *= R.I[op_row_b]
            thisVert *= dag(A[op_row_b, col])'
            thisVert *= L.I[op_row_a]
            thisVert *= ϕ
            thisVert *= R.I[op_row_a]
            thisVert *= AIL
            thisVert *= op_a
        elseif row == op_row_b
            low_row  = op_row_a - 1
            high_row = op_row_b
            AIL = low_row > 0 ? AI[:below][low_row] : dummy
            AIH = high_row < Ny ? AI[:above][end - high_row] : dummy
            sA = firstind(A[op_row_a, col], "Site")
            op_a = replaceind!(copy(H[opcode].ops[1]), H[opcode].site_ind, sA)
            op_a = replaceind!(op_a, H[opcode].site_ind', sA')
            sB = firstind(A[op_row_b, col], "Site")
            op_b = replaceind!(copy(H[opcode].ops[2]), H[opcode].site_ind, sB)
            op_b = replaceind!(op_b, H[opcode].site_ind', sB')
            thisVert  = AIL
            thisVert *= L.I[op_row_a]
            thisVert *= A[op_row_a, col] * op_a
            thisVert *= R.I[op_row_a]
            thisVert *= dag(A[op_row_a, col])'
            thisVert *= L.I[op_row_b]
            thisVert *= ϕ
            thisVert *= R.I[op_row_b]
            thisVert *= AIH
            thisVert *= op_b
        end
        @assert hasinds(inds(thisVert), AAinds) "inds of thisVert and AAinds differ!\n$(inds(thisVert))\n$AAinds\n"
        @assert hasinds(AAinds, inds(thisVert)) "inds of thisVert and AAinds differ!\n$(inds(thisVert))\n$AAinds\n"
        if hasinds(AAinds, inds(thisVert)) && hasinds(inds(thisVert), AAinds)
            push!(vTerms, thisVert)
        end
    end
    return vTerms
end

function diagonalTerms(A::fPEPS,
                       L::Environments,
                       R::Environments,
                       AI,
                       AD,
                       H,
                       row::Int,
                       col::Int,
                       ϕ::ITensor)::Vector{ITensor}
    Ny, Nx = size(A)
    is_cu  = is_gpu(A)
    dTerms = Vector{ITensor}(undef, length(H))
    AAinds = inds(prime(ϕ))
    dummy  = is_cu ? cuITensor(1.0) : ITensor(1.0)
    for opcode in 1:length(H)
        thisDiag = dummy
        op_row_a = H[opcode].sites[1][1]
        op_row_b = H[opcode].sites[2][1]
        op_col_a = H[opcode].sites[1][2]
        op_col_b = H[opcode].sites[2][2]
        side     = op_col_a < col
        max_op_row = max(op_row_a, op_row_b)
        min_op_row = min(op_row_a, op_row_b)
        if max_op_row < row || row < min_op_row
            local D, I
            if min_op_row > row
                D = AD[:above][opcode][end - row]
                I = row > 1 ? AD[:below][opcode][row - 1] : dummy
            elseif max_op_row < row
                D = AD[:below][opcode][row - 1]
                I = row < Ny ? AD[:above][opcode][end - row] : dummy
            end
            thisDiag *= D
            thisDiag *= side ? L.DiagInProgress[row, opcode] : L.I[row]
            thisDiag *= ϕ
            thisDiag *= side ? R.I[row] : R.DiagInProgress[row, opcode]
            thisDiag *= I
            thisDiag *= spinI(firstind(A[row, col], "Site"); is_gpu=is_cu)
            @assert hasinds(inds(thisDiag), AAinds) "inds of thisDiag and AAinds differ!\n$(inds(thisDiag))\n$AAinds\n"
            @assert hasinds(AAinds, inds(thisDiag)) "inds of thisDiag and AAinds differ!\n$(inds(thisDiag))\n$AAinds\n"
        else
            low_row  = min(op_row_a, op_row_b) - 1
            high_row = max(op_row_a, op_row_b)
            AIL      = low_row  >= 1 ? AD[:below][opcode][low_row]        : dummy
            AIH      = high_row < Ny ? AD[:above][opcode][end - high_row] : dummy
            sA       = firstind(A[op_row_a, col], "Site")
            sB       = firstind(A[op_row_b, col], "Site")
            local op_a, op_b
            if row == op_row_a
                if col == op_col_a
                    op_a     = replaceind!(copy(H[opcode].ops[1]), H[opcode].site_ind, sA)
                    op_a     = replaceind!(op_a, H[opcode].site_ind', sA')
                    op_b     = spinI(sB; is_gpu=is_cu)
                elseif col == op_col_b
                    op_b     = replaceind!(copy(H[opcode].ops[1]), H[opcode].site_ind, sB)
                    op_b     = replaceind!(op_b, H[opcode].site_ind', sB')
                    op_a     = spinI(sA; is_gpu=is_cu)
                end
                thisDiag  = op_row_b > op_row_a ? AIH : AIL
                thisDiag *= side ? L.DiagInProgress[op_row_b, opcode] : L.I[op_row_b]
                thisDiag *= A[op_row_b, col] * op_b
                thisDiag *= side ? R.I[op_row_b] : R.DiagInProgress[op_row_b, opcode]
                thisDiag *= dag(A[op_row_b, col])'
                thisDiag *= side ? L.DiagInProgress[op_row_a, opcode] : L.I[op_row_a]
                thisDiag *= ϕ
                thisDiag *= side ? R.I[op_row_a] : R.DiagInProgress[op_row_a, opcode]
                thisDiag *= op_row_b > op_row_a ? AIL : AIH
                thisDiag *= op_a
            elseif row == op_row_b
                if col == op_col_b
                    op_a     = spinI(sA; is_gpu=is_cu)
                    op_b     = replaceind!(copy(H[opcode].ops[2]), H[opcode].site_ind, sB)
                    op_b     = replaceind!(op_b, H[opcode].site_ind', sB')
                elseif col == op_col_a
                    op_a     = replaceind!(copy(H[opcode].ops[1]), H[opcode].site_ind, sA)
                    op_a     = replaceind!(op_a, H[opcode].site_ind', sA')
                    op_b     = spinI(sB; is_gpu=is_cu)
                end
                thisDiag  = op_row_b > op_row_a ? AIL : AIH
                thisDiag *= side ? L.DiagInProgress[op_row_a, opcode] : L.I[op_row_a]
                thisDiag *= A[op_row_a, col] * op_a
                thisDiag *= side ? R.I[op_row_a] : R.DiagInProgress[op_row_a, opcode]
                thisDiag *= dag(A[op_row_a, col])'
                thisDiag *= side ? L.DiagInProgress[op_row_b, opcode] : L.I[op_row_b]
                thisDiag *= ϕ
                thisDiag *= side ? R.I[op_row_b] : R.DiagInProgress[op_row_b, opcode]
                thisDiag *= op_row_b > op_row_a ? AIH : AIL
                thisDiag *= op_b
            end
            @assert hasinds(inds(thisDiag), AAinds) "inds of thisDiag and AAinds differ!\n$(inds(thisDiag))\n$AAinds\n"
            @assert hasinds(AAinds, inds(thisDiag)) "inds of thisDiag and AAinds differ!\n$(inds(thisDiag))\n$AAinds\n"
        end
        if hasinds(AAinds, inds(thisDiag)) && hasinds(inds(thisDiag), AAinds)
            dTerms[opcode] = thisDiag
        end
    end
    return dTerms
end

function fieldTerms(A::fPEPS,
                    L::Environments,
                    R::Environments,
                    AI,
                    AF,
                    H,
                    row::Int,
                    col::Int,
                    ϕ::ITensor)::Vector{ITensor}
    Ny, Nx = size(A)
    is_cu  = is_gpu(A)
    fTerms = Vector{ITensor}(undef, length(H))
    dummy  = is_cu ? cuITensor(1.0) : ITensor(1.0)
    AAinds = inds(prime(ϕ))
    @inbounds for opcode in 1:length(H)
        thisField = copy(dummy)
        op_row    = H[opcode].sites[1][1]
        if op_row != row
            local F, I
            if op_row > row
                F = AF[:above][opcode][end - row]
                I = row > 1 ? AI[:below][row - 1] : dummy
            else
                F = AF[:below][opcode][row - 1]
                I = row < Ny ? AI[:above][end - row] : dummy
            end
            thisField *= F
            thisField *= L.I[row]
            thisField *= ϕ
            thisField *= R.I[row]
            thisField *= I
            thisField *= spinI(firstind(A[row, col], "Site"); is_gpu=is_cu)
        else
            low_row  = op_row - 1
            high_row = op_row
            AIL = low_row > 0   ? AI[:below][low_row]        : dummy
            AIH = high_row < Ny ? AI[:above][end - high_row] : dummy
            thisField  = AIL
            thisField *= L.I[row]
            thisField *= ϕ
            thisField *= R.I[row]
            thisField *= AIH
            sA = firstind(A[row, col], "Site")
            op = copy(H[opcode].ops[1])
            op = replaceind!(op, H[opcode].site_ind, sA)
            op = replaceind!(op, H[opcode].site_ind', sA')
            thisField *= op
        end
        @assert hasinds(inds(thisField), AAinds)
        @assert hasinds(AAinds, inds(thisField))
        fTerms[opcode] = thisField
    end
    return fTerms
end

function connectLeftTerms(A::fPEPS,
                          L::Environments,
                          R::Environments,
                          AI, AL, H,
                          row::Int, col::Int,
                          ϕ::ITensor)::Vector{ITensor}
    Ny, Nx = size(A)
    is_cu  = is_gpu(A)
    lTerms = Vector{ITensor}(undef, length(H))
    AAinds = inds(prime(ϕ))
    dummy  = is_cu ? cuITensor(1.0) : ITensor(1.0)
    @inbounds for opcode in 1:length(H)
        op_row_b = H[opcode].sites[2][1]
        op_b = copy(H[opcode].ops[2])
        as   = firstind(A[op_row_b, col], "Site")
        op_b = replaceind!(op_b, H[opcode].site_ind, as)
        op_b = replaceind!(op_b, H[opcode].site_ind', as')
        thisHori = is_cu ? cuITensor(1.0) : ITensor(1.0)
        if op_row_b != row
            local ancL, I
            if op_row_b > row
                ancL = AL[:above][opcode][end - row]
                I = row > 1 ? AL[:below][opcode][row - 1] : dummy
            else
                ancL = AL[:below][opcode][row - 1]
                I = row < Ny ? AL[:above][opcode][end - row] : dummy
            end
            thisHori  = ancL
            thisHori *= L.InProgress[row, opcode]
            thisHori *= ϕ
           thisHori *= R.I[row]
            thisHori *= I
            thisHori *= spinI(firstind(A[row, col], "Site"); is_gpu=is_cu)
        else
            low_row = (op_row_b <= row) ? op_row_b - 1 : row - 1;
            high_row = (op_row_b >= row) ? op_row_b + 1 : row + 1;
            if low_row >= 1
                thisHori *= AL[:below][opcode][low_row]
            end
            thisHori *= R.I[row]
            thisHori *= ϕ
            thisHori *= L.InProgress[row, opcode]
            thisHori *= op_b
            if high_row <= Ny
                thisHori *= AL[:above][opcode][end - high_row + 1]
            end
        end
        @assert hasinds(inds(thisHori), AAinds)
        @assert hasinds(AAinds, inds(thisHori))
        lTerms[opcode] = thisHori
    end
    return lTerms
end

function connectRightTerms(A::fPEPS,
                           L::Environments,
                           R::Environments,
                           AI, AR, H,
                           row::Int, col::Int,
                           ϕ::ITensor)::Vector{ITensor}
    Ny, Nx = size(A)
    is_cu  = is_gpu(A)
    rTerms = Vector{ITensor}(undef, length(H))
    AAinds = inds(prime(ϕ))
    dummy  = is_cu ? cuITensor(1.0) : ITensor(1.0)
    @inbounds for opcode in 1:length(H)
        op_row_a = H[opcode].sites[1][1]
        op_a = copy(H[opcode].ops[1])
        as   = firstind(A[op_row_a, col], "Site")
        op_a = replaceind!(op_a, H[opcode].site_ind, as)
        op_a = replaceind!(op_a, H[opcode].site_ind', as')
        thisHori = is_cu ? cuITensor(1.0) : ITensor(1.0)
        if op_row_a != row
            local ancR, I
            if op_row_a > row
                ancR = AR[:above][opcode][end - row]
                I = row > 1 ? AR[:below][opcode][row - 1] : dummy
            else
                ancR = AR[:below][opcode][row - 1]
                I = row < Ny ? AR[:above][opcode][end - row] : dummy
            end
            thisHori = ancR
            thisHori *= R.InProgress[row, opcode]
            thisHori *= ϕ
            thisHori *= L.I[row]
            thisHori *= I
            thisHori *= spinI(firstind(A[row, col], "Site"); is_gpu=is_cu)
        else
            low_row = (op_row_a <= row) ? op_row_a - 1 : row - 1;
            high_row = (op_row_a >= row) ? op_row_a + 1 : row + 1;
            if low_row >= 1
                thisHori *= AR[:below][opcode][low_row]
            end
            thisHori *= L.I[row]
            thisHori *= R.InProgress[row, opcode]
            thisHori *= ϕ
            thisHori *= op_a
            if high_row <= Ny
                thisHori *= AR[:above][opcode][end - high_row + 1]
            end
        end
        @assert hasinds(inds(thisHori), AAinds)
        @assert hasinds(AAinds, inds(thisHori))
        rTerms[opcode] = thisHori
    end
    return rTerms
end

function buildLocalH(A::fPEPS,
                     L::Environments, R::Environments,
                     AncEnvs, H,
                     row::Int, col::Int,
                     ϕ::ITensor;
                     verbose::Bool=false)
    field_H_terms = getDirectional(vcat(H[:, col]...), Field)
    vert_H_terms  = getDirectional(vcat(H[:, col]...), Vertical)
    term_count    = 1 + length(field_H_terms) + length(vert_H_terms)
    Ny, Nx        = size(A)
    @timeit "build Ns" begin
        N   = buildN(A, L, R, AncEnvs[:I], row, col, ϕ)
    end
    den = scalar(collect(N * dag(ϕ)'))
    if verbose
        println("--- NORM ---")
        println(den)
    end
    local left_H_terms, right_H_terms, diag_left_H_terms, diag_right_H_terms
    if col > 1
        left_H_terms = getDirectional(vcat(H[:, col - 1]...), Horizontal)
        term_count += length(left_H_terms)
        diag_left_H_terms = getDirectional(vcat(H[:, col - 1]...), Diag)
        term_count += length(diag_left_H_terms)
    end
    if col < Nx
        right_H_terms = getDirectional(vcat(H[:, col]...), Horizontal)
        term_count += length(right_H_terms)
        diag_right_H_terms = getDirectional(vcat(H[:, col]...), Diag)
        term_count += length(diag_right_H_terms)
    end
    if 1 < col < Nx
        term_count += 1
    end
    Hs = Vector{ITensor}(undef, term_count)
    term_counter = 1
    @debug "\t\tBuilding I*H and H*I row $row col $col"
    @timeit "build HIs" begin
        HIs = buildHIs(A, L, R, row, col, ϕ)
        Hs[term_counter] = HIs[1]
        if length(HIs) == 2
            Hs[term_counter+1] = HIs[2]
        end
        term_counter += length(HIs)
        if verbose
            println("--- HI TERMS ---")
            for HI in HIs
                println(scalar(HI * dag(ϕ)'))
            end
        end
    end
    @debug "\t\tBuilding vertical H terms row $row col $col"
    @timeit "build vertical terms" begin
        vTs = verticalTerms(A, L, R, AncEnvs[:I], AncEnvs[:V], vert_H_terms, row, col, ϕ)
        Hs[term_counter:term_counter+length(vTs) - 1] = vTs
        term_counter += length(vTs)
        if verbose
            println( "--- vT TERMS col $col row $row ---")
            for (vv, vT) in enumerate(vTs)
                println(vert_H_terms[vv])
                println(scalar(vT * dag(ϕ)')/den)
            end
        end
    end
    @debug "\t\tBuilding field H terms row $row col $col"
    @timeit "build field terms" begin
        fTs = fieldTerms(A, L, R, AncEnvs[:I], AncEnvs[:F], field_H_terms, row, col, ϕ)
        Hs[term_counter:term_counter+length(fTs) - 1] = fTs[:]
        term_counter += length(fTs)
        if verbose
            println( "--- fT TERMS col $col row $row ---")
            for (ff, fT) in enumerate(fTs)
                println(field_H_terms[ff])
                println(scalar(fT * dag(ϕ)')/den)
            end
        end
    end
    if col > 1
        @debug "\t\tBuilding left H terms row $row col $col"
        @timeit "build left terms" begin
            lTs = connectLeftTerms(A, L, R, AncEnvs[:I], AncEnvs[:L], left_H_terms, row, col, ϕ)
            Hs[term_counter:term_counter+length(lTs) - 1] = lTs[:]
            term_counter += length(lTs)
            if verbose
                println( "--- lT TERMS col $col row $row ---")
                for (ll, lT) in enumerate(lTs)
                    println(left_H_terms[ll])
                    println(scalar(lT * dag(ϕ)')/den)
                end
            end
        end
        @debug "\t\tBuilt left terms"
        @debug "\t\tBuilding left diag H terms row $row col $col"
        @timeit "build diag left terms" begin
            lTs = diagonalTerms(A, L, R, AncEnvs[:I], AncEnvs[:DL], diag_left_H_terms, row, col, ϕ)
            Hs[term_counter:term_counter+length(lTs) - 1] = lTs[:]
            term_counter += length(lTs)
            if verbose
                println( "--- diag lT TERMS col $col row $row ---")
                for (ll, lT) in enumerate(lTs)
                    println(left_H_terms[ll])
                    println(scalar(lT * dag(ϕ)')/den)
                end
            end
        end
        @debug "\t\tBuilt diag left terms"
    end
    if col < Nx
        @debug "\t\tBuilding right H terms row $row col $col"
        @timeit "build right terms" begin
            rTs = connectRightTerms(A, L, R, AncEnvs[:I], AncEnvs[:R], right_H_terms, row, col, ϕ)
            Hs[term_counter:term_counter+length(rTs) - 1] = rTs[:]
            term_counter += length(rTs)
            if verbose
                println( "--- rT TERMS col $col row $row ---")
                for (rr, rT) in enumerate(rTs)
                    println(right_H_terms[rr])
                    println(scalar(rT * dag(ϕ)')/den)
                end
            end
        end
        @debug "\t\tBuilding right diag H terms row $row col $col"
        @timeit "build diag right terms" begin
            rTs = diagonalTerms(A, L, R, AncEnvs[:I], AncEnvs[:DR], diag_right_H_terms, row, col, ϕ)
            Hs[term_counter:term_counter+length(rTs) - 1] = rTs[:]
            term_counter += length(rTs)
            if verbose
                println( "--- diag rT TERMS col $col row $row ---")
                for (rr, rT) in enumerate(rTs)
                    println(right_H_terms[rr])
                    println(scalar(rT * dag(ϕ)')/den)
                end
            end
        end
        @debug "\t\tBuilt right terms"
    end
    return Hs, N
end

function buildAncs(A::fPEPS, L::Environments, R::Environments, H, col::Int)
    Ny, Nx = size(A)
    is_cu  = is_gpu(A)
    dummy  = is_cu ? cuITensor(1.0) : ITensor(1.0)
    @timeit "makeAncIs" begin
        Ia = makeAncillaryIs(A, L, R, col)
        Ib = Vector{ITensor}(undef, Ny)
        Is = (above=Ia, below=Ib)
    end
    @timeit "makeAncVs" begin
        vH = getDirectional(vcat(H[:, col]...), Vertical)
        Va = makeAncillaryVs(A, L, R, vH, col)
        Vb = [Vector{ITensor}() for ii in 1:length(Va)]
        Vs = (above=Va, below=Vb)
    end
    @timeit "makeAncFs" begin
        fH = getDirectional(vcat(H[:, col]...), Field)
        Fa = makeAncillaryFs(A, L, R, fH, col)
        Fb = [Vector{ITensor}() for ii in 1:length(Fa)]
        Fs = (above=Fa, below=Fb)
    end
    Ls = (above=Vector{ITensor}(), below=Vector{ITensor}())
    Rs = (above=Vector{ITensor}(), below=Vector{ITensor}())
    DLs = (above=Vector{ITensor}(), below=Vector{ITensor}())
    DRs = (above=Vector{ITensor}(), below=Vector{ITensor}())
    if col > 1
        @timeit "makeAncLs" begin
            lH = getDirectional(vcat(H[:, col-1]...), Horizontal)
            La = makeAncillarySide(A, L, R, lH, col, :left)
            Lb = [Vector{ITensor}() for ii in  1:length(La)]
            Ls = (above=La, below=Lb)
        end
        @timeit "makeAncDLs" begin
            lH = getDirectional(vcat(H[:, col-1]...), Diag)
            La = makeAncillarySideDiag(A, L, R, lH, col, :left)
            Lb = [Vector{ITensor}() for ii in  1:length(La)]
            DLs = (above=La, below=Lb)
        end
    end
    if col < Nx
        @timeit "makeAncRs" begin
            rH = getDirectional(vcat(H[:, col]...), Horizontal)
            Ra = makeAncillarySide(A, R, L, rH, col, :right)
            Rb = [Vector{ITensor}() for ii in  1:length(Ra)]
            Rs = (above=Ra, below=Rb)
        end
        @timeit "makeAncRs" begin
            rH = getDirectional(vcat(H[:, col]...), Diag)
            Ra = makeAncillarySideDiag(A, R, L, rH, col, :right)
            Rb = [Vector{ITensor}() for ii in  1:length(Ra)]
            DRs = (above=Ra, below=Rb)
        end
    end
    Ancs = (I=Is, V=Vs, F=Fs, L=Ls, R=Rs, DL=DLs, DR=DRs)
    return Ancs
end

function updateAncs(A::fPEPS,
                    L::Environments, R::Environments,
                    AncEnvs, H,
                    row::Int, col::Int)
    Ny, Nx = size(A)
    is_cu  = is_gpu(A)

    Is, Vs, Fs, Ls, Rs, DLs, DRs = AncEnvs
    @debug "\tUpdating ancillary identity terms for col $col row $row"
    Ib = updateAncillaryIs(A, Is[:below], L, R, row, col)
    Is = (above=Is[:above], below=Ib)

    @debug "\tUpdating ancillary vertical terms for col $col row $row"
    vH = getDirectional(vcat(H[:, col]...), Vertical)
    Vb = updateAncillaryVs(A, Vs[:below], Ib, L, R, vH, row, col)
    Vs = (above=Vs[:above], below=Vb)

    @debug "\tUpdating ancillary field terms for col $col row $row"
    fH = getDirectional(vcat(H[:, col]...), Field)
    Fb = updateAncillaryFs(A, Fs[:below], Ib, L, R, fH, row, col)
    Fs = (above=Fs[:above], below=Fb)

    if col > 1
        @debug "\tUpdating ancillary left terms for col $col row $row"
        lH = getDirectional(vcat(H[:, col-1]...), Horizontal)
        Lb = updateAncillarySide(A, Ls[:below], Ib, L, R, lH, row, col, :left)
        Ls = (above=Ls[:above], below=Lb)
        @debug "\tUpdating ancillary diag left terms for col $col row $row"
        lH = getDirectional(vcat(H[:, col-1]...), Diag)
        DLb = updateAncillarySideDiag(A, DLs[:below], Ib, L, R, lH, row, col, :left)
        DLs = (above=DLs[:above], below=DLb)
    end
    if col < Nx
        @debug "\tUpdating ancillary right terms for col $col row $row"
        rH = getDirectional(vcat(H[:, col]...), Horizontal)
        Rb = updateAncillarySide(A, Rs[:below], Ib, R, L, rH, row, col, :right)
        Rs = (above=Rs[:above], below=Rb)
        @debug "\tUpdating ancillary diag right terms for col $col row $row"
        rH = getDirectional(vcat(H[:, col]...), Diag)
        DRb = updateAncillarySideDiag(A, DRs[:below], Ib, R, L, rH, row, col, :right)
        DRs = (above=DRs[:above], below=DRb)
    end
    Ancs = (I=Is, V=Vs, F=Fs, L=Ls, R=Rs, DL=DLs, DR=DRs)
    return Ancs
end

struct ITensorMap
  A::fPEPS
  H::Matrix{Vector{Operator}}
  L::Environments
  R::Environments
  AncEnvs::NamedTuple
  row::Int
  col::Int
end
Base.eltype(M::ITensorMap)  = eltype(M.A)
Base.size(M::ITensorMap)    = ITensors.dim(M.A[M.row, M.col])
function (M::ITensorMap)(v::ITensor)
    Hs, N  = buildLocalH(M.A, M.L, M.R, M.AncEnvs, M.H, M.row, M.col, v)
    localH = sum(Hs)
    return noprime(localH)
end

function optimizeLocalH(A::fPEPS,
                        L::Environments, R::Environments,
                        AncEnvs, H,
                        row::Int, col::Int;
                        kwargs...)
    Ny, Nx        = size(A)
    is_cu         = is_gpu(A)
    field_H_terms = getDirectional(vcat(H[:, col]...), Field)
    vert_H_terms  = getDirectional(vcat(H[:, col]...), Vertical)
    @debug "\tBuilding H for col $col row $row"
    @timeit "build H" begin
        Hs, N = buildLocalH(A, L, R, AncEnvs, H, row, col, A[row, col])
    end
    initial_N = real(scalar(collect(N * dag(A[row, col])')))
    @timeit "sum H terms" begin
        localH = sum(Hs)
    end
    initial_E = real(scalar(collect(deepcopy(localH) * dag(A[row, col])')))
    @info "Initial energy at row $row col $col : $(initial_E/(initial_N*Nx*Ny)) and norm : $initial_N"
    @debug "\tBeginning davidson for col $col row $row"
    mapper   = ITensorMap(A, H, L, R, AncEnvs, row, col)
    λ, new_A = davidson(mapper, A[row, col]; maxiter=1, kwargs...)
    new_E    = λ #real(scalar(collect(new_A * localH * dag(new_A)')))
    N        = buildN(A, L, R, AncEnvs[:I], row, col, new_A)
    new_N    = real(scalar(collect(N * dag(new_A)')))
    @info "Optimized energy at row $row col $col : $(new_E/(new_N*Nx*Ny)) and norm : $new_N"
    #println("Optimized energy at row $row col $col : $(new_E/(new_N*Nx*Ny)) and norm : $new_N")
    @timeit "restore intraColumnGauge" begin
        if row < Ny
            @debug "\tRestoring intraColumnGauge for col $col row $row"
            cmb_is     = IndexSet(firstind(A[row, col], "Site"))
            if col > 1
                cmb_is = IndexSet(cmb_is..., commonind(A[row, col], A[row, col - 1]))
            end
            if col < Nx
                cmb_is = IndexSet(cmb_is..., commonind(A[row, col], A[row, col + 1]))
            end
            cmb     = combiner(cmb_is, tags="CMB")
            ci      = combinedind(cmb)
            Lis     = IndexSet(ci) #cmb_is
            if row > 1
                Lis = IndexSet(Lis..., commonind(A[row, col], A[row - 1, col]))
            end
            old_ci  = commonind(A[row, col], A[row+1, col])
            svdA    = new_A*cmb
            Ris     = uniqueinds(inds(svdA), Lis)
            U, S, V = svd(svdA, Ris; kwargs...)
            new_ci  = commonind(V, S)
            replaceind!(V, new_ci, old_ci)
            A[row, col]    = V*cmb
            newU = S*U*A[row+1, col]
            replaceind!(newU, new_ci, old_ci)
            A[row+1, col] = newU
            if row < Ny - 1
                nI    = spinI(firstind(A[row+1, col], "Site"); is_gpu=is_cu)
                newAA = AncEnvs[:I][:above][end - row - 1]
                newAA *= L.I[row+1]
                newAA *= A[row+1, col] * nI
                newAA *= R.I[row+1]
                newAA *= dag(A[row+1, col])'
                AncEnvs[:I][:above][end - row] = newAA
            end
        else
            A[row, col] = initial_N > 0 ? new_A/√new_N : new_A
        end
    end
    return A, AncEnvs
end

function sweepColumn(A::fPEPS,
                     L::Environments, R::Environments,
                     H,
                     col::Int;
                     kwargs...)
    Ny, Nx = size(A)
    @timeit "intraColumnGauge" begin
        A = intraColumnGauge(A, col; kwargs...)
    end
    if col == div(Nx,2)
        L_s = buildLs(A, H; kwargs...)
        R_s = buildRs(A, H; kwargs...)
        EAncEnvs = buildAncs(A, L_s[col - 1], R_s[col + 1], H, col)
        N, E = measureEnergy(A, L_s[col - 1], R_s[col + 1], EAncEnvs, H, 1, col)
        println("Energy at MID: ", E/(Nx*Ny))
        println("Nx: ", Nx)
    end
    @debug "Beginning buildAncs for col $col"
    two_site_vertical::Bool   = get(kwargs, :two_site_vertical, false)
    if two_site_vertical
        AncEnvs = buildAncs(A, L, R, H, col)
        @inbounds for row in 1:Ny-1
            if row > 1
                @timeit "updateAncs" begin
                    AncEnvs = updateAncs(A, L, R, AncEnvs, H, row-1, col)
                end
            end
            @debug "Beginning optimizing H for col $col"
            A, AncEnvs = optimizeLocalHTwoSiteVertical(A, L, R, AncEnvs, H, row, col; kwargs...)
        end
    else
        AncEnvs = buildAncs(A, L, R, H, col)
        @inbounds for row in 1:Ny
            if row > 1
                @timeit "updateAncs" begin
                    AncEnvs = updateAncs(A, L, R, AncEnvs, H, row-1, col)
                end
            end
            @debug "Beginning optimizing H for col $col"
            A, AncEnvs = optimizeLocalH(A, L, R, AncEnvs, H, row, col; kwargs...)
        end
    end
    return A
end


text:
module Err
export PatternUnsolvedException, InternalException, SyntaxError
struct PatternUnsolvedException <: Exception
    pattern :: Expr
end

struct InternalException <: Exception
    msg :: String
end

struct SyntaxError <: Exception
    msg :: String
end

struct DataTypeUsageError <: Exception
    msg :: String
end

end


text:
module TestExamples2DMHD

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "2d")

@testset "MHD" begin
  @testset "elixir_mhd_alfven_wave.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"),
      l2   = [0.000112085405989522, 5.908049046758828e-6, 5.908049046772597e-6, 8.493815951190774e-6, 1.297306491315482e-6, 1.2316941466595997e-6, 1.2316941466702522e-6, 1.8368487700950752e-6, 4.4903164608485606e-7],
      linf = [0.0002693844074639351, 1.6289386154083596e-5, 1.6289386155082797e-5, 2.747802438275715e-5, 5.344948833307939e-6, 8.137608476288527e-6, 8.137608476843639e-6, 1.2121292898431557e-5, 2.306551665486275e-6])
  end

  @testset "elixir_mhd_alfven_wave_mortar.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave_mortar.jl"),
      l2   = [3.843336695097301e-6, 1.465357374916213e-6, 1.416088654606427e-6, 1.6806652628357514e-6, 1.4256305721007091e-6, 1.4507398918063128e-6, 1.397023085105777e-6, 1.704347968034822e-6, 8.210831113360175e-7],
      linf = [3.4691159990107856e-5, 1.4347253854782305e-5, 1.39581811973849e-5, 1.4252304932890758e-5, 8.244903569987194e-6, 1.0371383968976744e-5, 1.1133207181823757e-5, 1.600372637142189e-5, 5.572707635298199e-6],
      tspan = (0.0, 1.0))
  end

  @testset "elixir_mhd_ec.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ec.jl"),
      l2   = [0.03628262478842946, 0.043011555956537974, 0.042997357872762064, 0.02574553086994775, 0.1620610267679157, 0.01745438427233919, 0.017455233781000938, 0.026879482381500206, 0.00020380087569640704],
      linf = [0.2351786266785303, 0.3156206778985494, 0.30941696929809714, 0.21167563519254826, 0.9688251298546073, 0.09076254289154928, 0.0916058976949814, 0.15698032974768517, 0.006131914796905391])
  end

  @testset "elixir_mhd_orszag_tang.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_orszag_tang.jl"),
      l2   = [0.21671818699492046, 0.26364533451882755, 0.31388339884931576, 0.0, 0.5122729051056656, 0.2291828683723837, 0.34308184804821307, 0.0, 0.003210081660131638],
      linf = [1.2648775120640154, 0.673836824735164, 0.8585376516993447, 0.0, 2.8148181376188246, 0.657323372510843, 0.9595541188847632, 0.0, 0.0525390951527749],
      tspan = (0.0, 0.09))
  end

  @testset "elixir_mhd_orszag_tang.jl with flux_hll" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_orszag_tang.jl"),
      l2   = [0.10793703703273708, 0.20177434827232335, 0.2296817363329907, 0.0, 0.2994119440979309, 0.15675567943351448, 0.24281245338083307, 0.0, 0.003548624701792127],
      linf = [0.560160833798496, 0.5098933680011135, 0.6566913038761877, 0.0, 0.9905839416293134, 0.39936698379939284, 0.6734754381366532, 0.0, 0.12739518123975463],
      tspan = (0.0, 0.06), surface_flux = flux_hll)
  end
end

end # module


text:
abstract SmoothingKernel

type MomentFunction
  gᵢ::Function      ## Moment Function
  sᵢ::Function      ## Smoothed moment function
  sn::Function      ## (m×1) ∑pᵢ sᵢ
  ∂∑pᵢsᵢ::Function  ## (k×m)
  ∂sᵢλ::Function    ## (n×k)
  ∂²sᵢλ::Function   ## (kxk)
  kern::SmoothingKernel
  nobs::Int64
  nmom::Int64
  npar::Int64
end

type MinimumDivergenceProblem <: MathProgSolverInterface.AbstractNLPEvaluator
  momf::MomentFunction
  div::Divergence
  nobs::Int64
  nmom::Int64
  npar::Int64
  gele::Int64
  hele::Int64
  lmbd::Array{Float64, 1}
  wght::Array{Float64, 1}
end

typealias MDP MinimumDivergenceProblem

type MinimumDivergenceProblem <: MathProgSolverInterface.AbstractNLPEvaluator
  momf::MomentFunction
  div::Divergence
  nobs::Int64
  nmom::Int64
  npar::Int64
  gele::Int64
  hele::Int64
  lmbd::Array{Float64, 1}
  wght::Array{Float64, 1}
end


text:
push!(LOAD_PATH,"../src/")

using Documenter

using Genie, Genie.App, Genie.AppServer, Genie.Assets
using Genie.Cache, Genie.Commands, Genie.Configuration, Genie.Cookies
using Genie.Deploy, Genie.Encryption, Genie.Exceptions
using Genie.FileTemplates, Genie.Flash, Genie.Generator
using Genie.Headers, Genie.HTTPUtils, Genie.Input, Genie.Plugins
using Genie.Renderer, Genie.Requests, Genie.Responses, Genie.Router
using Genie.Sessions, Genie.Toolbox, Genie.Util, Genie.WebChannels

push!(LOAD_PATH,  "../../src",
                  "../../src/cache_adapters",
                  "../../src/session_adapters",
                  "../../src/renderers")

makedocs(
    sitename = "Genie - The Highly Productive Julia Web Framework",
    format = Documenter.HTML(prettyurls = false),
    pages = [
        "Home" => "index.md",
        "Guides" => [
          "Working with Genie Apps" => "guides/Working_With_Genie_Apps.md",
          "Using Genie in an interactive environment" => "guides/Interactive_environment.md",
          "Developing an API backend" => "guides/Simple_API_backend.md",
          "Using Genie Plugins" => "guides/Genie_Plugins.md",
          "Working With Genie Apps: Intermediate Topics [WIP]" => "guides/Working_With_Genie_Apps_Intermediary_Topics.md"
        ],
        "Tutorials" => [
          "Welcome to Genie"  => "tutorials/1--Overview.md",
          "Installing Genie"  => "tutorials/2--Installing_Genie.md",
          "Getting started"   => "tutorials/3--Getting_Started.md",
          "Creating a web service" => "tutorials/4--Developing_Web_Services.md",
          "Developing MVC web applications" => "tutorials/4-1--Developing_MVC_Web_Apps.md",
          "Handling URI/query params" => "tutorials/5--Handling_Query_Params.md",
          "Working with forms and POST payloads" => "tutorials/6--Working_with_POST_Payloads.md",
          "Using JSON payloads" => "tutorials/7--Using_JSON_Payloads.md",
          "Uploading files" => "tutorials/8--Handling_File_Uploads.md",
          "Adding your libraries into Genie" => "tutorials/9--Publishing_Your_Julia_Code_Online_With_Genie_Apps.md",
          "Loading and starting Genie apps" => "tutorials/10--Loading_Genie_Apps.md",
          "Managing Genie app's dependencies" => "tutorials/11--Managing_External_Packages.md",
          "Advanced routing" => "tutorials/12--Advanced_Routing_Techniques.md",
          "Auto-loading configuration code with initializers" => "tutorials/13--Initializers.md",
          "The secrets file" => "tutorials/14--The_Secrets_File.md",
          "Auto-loading user libraries" => "tutorials/15--The_Lib_Folder.md",
          "Using Genie with Docker" => "tutorials/16--Using_Genie_With_Docker.md",
          "Working with WebSockets" => "tutorials/17--Working_with_Web_Sockets.md",
          "Force compiling route handlers" => "tutorials/80--Force_Compiling_Routes.md",
          "Deploying to Heroku with Buildpacks" => "tutorials/90--Deploying_With_Heroku_Buildpacks.md",
          "Deploying to Heroku with Docker" => "tutorials/91--Deploying_Genie_Docker_Apps_on_Heroku.md",
          "Deploying to server with Nginx" => "tutorials/92--Deploying_Genie_Server_Apps_with_Nginx.md"
        ],
        "API" => [
          "App" => "API/app.md",
          "AppServer" => "API/appserver.md",
          "Assets" => "API/assets.md",
          "Cache" => "API/cache.md",
          "Commands" => "API/commands.md",
          "Configuration" => "API/configuration.md",
          "Cookies" => "API/cookies.md",
          "Deploy" => [
            "Docker" => "API/deploy_docker.md",
            "Heroku" => "API/deploy_heroku.md"
          ],
          "Encryption" => "API/encryption.md",
          "Exceptions" => "API/exceptions.md",
          "FileTemplates" => "API/filetemplates.md",
          "Flash" => "API/flash.md",
          "Generator" => "API/generator.md",
          "Genie" => "API/genie.md",
          "Headers" => "API/headers.md",
          "HttpUtils" => "API/httputils.md",
          "Input" => "API/input.md",
          "Plugins" => "API/plugins.md",
          "Renderer" => "API/renderer.md",
          "HTML Renderer" => "API/renderer_html.md",
          "JS Renderer" => "API/renderer_js.md",
          "JSON Renderer" => "API/renderer_json.md",
          "Requests" => "API/requests.md",
          "Responses" => "API/responses.md",
          "Router" => "API/router.md",
          "Sessions" => "API/sessions.md",
          "Toolbox" => "API/toolbox.md",
          "Util" => "API/util.md",
          "WebChannels" => "API/webchannels.md"
        ]
    ],
)

deploydocs(
  repo = "github.com/GenieFramework/Genie.jl.git",
)


text:
using MPI

MPI.Initialized() || MPI.Init()

     comm = MPI.COMM_WORLD
mpi_ranks = MPI.Comm_size(comm)

@assert mpi_ranks == 4

using Statistics
using Oceananigans
using Oceananigans.Distributed

     ranks = (2, 2, 1)
      topo = (Periodic, Periodic, Bounded)
 full_grid = RegularRectilinearGrid(topology=topo, size=(128, 128, 1), extent=(4π, 4π, 1), halo=(3, 3, 3))
      arch = MultiCPU(grid=full_grid, ranks=ranks)
local_rank = MPI.Comm_rank(MPI.COMM_WORLD)

model = DistributedShallowWaterModel(
                  architecture = arch,
                          grid = full_grid,
                   timestepper = :RungeKutta3,
                     advection = UpwindBiasedFifthOrder(),
    gravitational_acceleration = 1.0
)

set!(model, h=model.grid.Lz)

uh₀ = rand(size(model.grid)...);
uh₀ .-= mean(uh₀);
set!(model, uh=uh₀, vh=uh₀)

progress(sim) = @info "Iteration: $(sim.model.clock.iteration), time: $(sim.model.clock.time)"
simulation = Simulation(model, Δt=0.001, stop_time=100.0, iteration_interval=1, progress=progress)

uh, vh, h = model.solution
outputs = (ζ=ComputedField(∂x(vh/h) - ∂y(uh/h)),)
filepath = "mpi_shallow_water_turbulence_rank$(local_rank).nc"
simulation.output_writers[:fields] =
    NetCDFOutputWriter(model, outputs, filepath=filepath, schedule=TimeInterval(1.0), mode="c")

MPI.Barrier(MPI.COMM_WORLD)

run!(simulation)

using Printf
using NCDatasets
using CairoMakie

nranks = prod(ranks)

if local_rank == 0

    ds = [NCDataset("mpi_shallow_water_turbulence_rank$r.nc") for r in 0:nranks-1]

    frame = Node(1)
    plot_title = @lift @sprintf("Oceananigans.jl + MPI: 2D turbulence t = %.2f", ds[1]["time"][$frame])
    ζ = [@lift ds[r]["ζ"][:, :, 1, $frame] for r in 1:nranks]

    fig = Figure(resolution=(1600, 1200))

    for rx in 1:ranks[1], ry in 1:ranks[2]
        ax = fig[rx, ry] = Axis(fig)
        r = (ry-1)*ranks[2] + rx
        hm = CairoMakie.heatmap!(ax, ds[r]["xF"], ds[r]["yF"], ζ[r], colormap=:balance, colorrange=(-2, 2))
    end

    record(fig, "mpi_shallow_water_turbulence.mp4", 1:length(ds[1]["time"])-1, framerate=30) do n
        @info "Animating MPI turbulence frame $n/$(length(ds[1]["time"]))..."
        frame[] = n
    end

    [close(d) for d in ds]
end


text:
import SemanticModels.ModelTools: model
import Base: show

"""    SimpleProblem <: AbstractProblem

represents a generic scientific model, where there are blocks of code that get run in order some of which
are function definitions. This type assumes a minimal structure on code. Namely,

1. There are imports/usingf expressions at the top of the expression
2. Code is broken into chunks with either function defintions or begin/end pairs
3. There is a single entrypoint function that "runs the model" This defaults to main().
"""
mutable struct SimpleProblem <: AbstractProblem
    expr::Expr
    imports::Any
    blocks::Vector{Expr}
    functions::Vector{Expr}
    entry::Expr
end

entrypoint(m::AbstractProblem) = m.entry
function entryname(m::AbstractProblem)
    c = entrypoint(m)
    iscall(c) || error("Entrypoint is not a :call")
    if iscall(c)
        return c.args[1]
    elseif isa(c, Symbol)
        return c
    end
end

function body(expr)
    if expr.head == :block
        return body(expr.args[end])
    elseif expr.head == :module
        return expr.args[end]
    end
    return :nothing
end

function model(::Type{SimpleProblem}, expr::Expr, entrypoint=:(main()))
    if expr.head == :block
        expr = expr.args[end]
    end
    b = body(expr)
    statements = b.args
    imports = filter(or(isusing,isimport), statements)
    blocks = filter(or(isblock, isfunc), statements)
    funcs = filter(isfunc, statements)
    return SimpleProblem(expr, imports, blocks, funcs, entrypoint)
end

function show(io::IO, m::SimpleProblem)
    fns = funcname.(m.functions)
    write(io, "SimpleProblem(\n")
    write(io, "  imports=$(repr(m.imports)),\n")
    write(io, "  blocks=$(repr(m.blocks)),\n")
    write(io, "  functions=$(repr(fns))\n")
    write(io, "  entry=$(repr(m.entry))")
end
# -


text:
function NonLinearEvolution(arg0::jdouble, arg1::RealVector, arg2::RealMatrix, arg3::RealMatrix, arg4::RealMatrix)
    return NonLinearEvolution((jdouble, RealVector, RealMatrix, RealMatrix, RealMatrix), arg0, arg1, arg2, arg3, arg4)
end

function get_current_state(obj::NonLinearEvolution)
    return jcall(obj, "getCurrentState", RealVector, ())
end

function get_current_time(obj::NonLinearEvolution)
    return jcall(obj, "getCurrentTime", jdouble, ())
end

function get_measurement_jacobian(obj::NonLinearEvolution)
    return jcall(obj, "getMeasurementJacobian", RealMatrix, ())
end

function get_process_noise_matrix(obj::NonLinearEvolution)
    return jcall(obj, "getProcessNoiseMatrix", RealMatrix, ())
end

function get_state_transition_matrix(obj::NonLinearEvolution)
    return jcall(obj, "getStateTransitionMatrix", RealMatrix, ())
end


text:
module GridPaths

using Vec
using Parameters

export
    GridTransition,
        WEST,EAST,NORTH,SOUTH,WAIT,
    get_translation,
    get_heading_angle,
    get_heading_vector,
    get_angular_offset,
    GridWaypoint,
    GridWorldPath,
    construct_grid_world_path

@enum GridTransition begin
    EAST    = 0
    NORTH   = 1
    WEST    = 2
    SOUTH   = 3
    WAIT    = 4
end
"""
    `get_translation`

    Returns the direction associated with a particular `GridTransition`
"""
function get_translation(a::GridTransition)
    if a == WEST
        return VecE2(-1.0, 0.0)
    elseif a == EAST
        return VecE2(1.0, 0.0)
    elseif a == NORTH
        return VecE2(0.0,1.0)
    elseif a == SOUTH
        return VecE2(0.0,-1.0)
    elseif a == WAIT
        return VecE2(0.0,0.0)
    end
    return VecE2(0.0,0.0)
end
"""
    `get_heading_angle`

    Returns the heading angle corresponding to a `GridTransition`s
"""
function get_heading_angle(a::GridTransition)
    if a == WAIT
        return 0.0
    elseif a == EAST
        return 0.0
    elseif a == NORTH
        return π/2
    elseif a == WEST
        return π
    elseif a == SOUTH
        return 3π/2
    end
    return NaN
end
"""
    `get_heading_vector`

    Returns the heading vector corresponding to a `GridTransition`s
"""
function get_heading_vector(a::GridTransition)
    if a == WAIT
        return VecE2(0.0,0.0)
    elseif a == EAST
        return VecE2(1.0,0.0)
    elseif a == NORTH
        return VecE2(0.0,1.0)
    elseif a == WEST
        return VecE2(-1.0,0.0)
    elseif a == SOUTH
        return VecE2(0.0,-1.0)
    end
    return VecE2(NaN,NaN)
end
"""
    `get_angular_offset`

    Returns the angular offset between two `GridTransition`s
"""
function get_angular_offset(θ1,θ2)
    Δθ = θ2 - θ1
    if Δθ > π
        Δθ = Δθ - 2π
    elseif Δθ < -π
        Δθ = Δθ + 2π
    end
    return Δθ
end

struct GridWaypoint
    pt::VecE2
    t::Float64
    transition::GridTransition
end
GridWaypoint(x::Float64,y::Float64,t::Float64,a::GridTransition) = GridWaypoint(VecE2(x,y),t,a)

"""
    `GridWorldPath`
"""
@with_kw struct GridWorldPath
    start_pt::VecE2                 = VecE2()
    start_time::Float64             = 0.0
    waypoints::Vector{GridWaypoint} = Vector{GridWaypoint}()
    # environment parameters
    cellwidth::Float64              = 1.0
end
function construct_grid_world_path(
    start_pt::VecE2,
    start_time::Float64,
    cellwidth::Float64,
    transition_time::Float64,
    transitions::Vector{GridTransition})
    waypoints = Vector{GridWaypoint}()
    pt = start_pt
    t = start_time
    for a in transitions
        next_pt = pt + cellwidth * get_translation(a)
        next_t = t + transition_time
        push!(waypoints, GridWaypoint(next_pt, next_t, a))
        pt = next_pt
        t = next_t
    end
    GridWorldPath(
        start_pt=start_pt,
        start_time=start_time,
        waypoints=waypoints,
        cellwidth=cellwidth)
end
function construct_grid_world_path(
    start_pt::VecE2,
    start_time::Float64,
    transitions::Vector{GridTransition},
    cellwidth::Float64,
    transition_time::Float64=1.0)
    construct_grid_world_path(start_pt,start_time,cellwidth,transition_time,transitions)
end
function construct_grid_world_path(
    x0::Float64,
    y0::Float64,
    t0::Float64,
    cellwidth::Float64,
    transition_time::Float64,
    transitions::Vector{GridTransition}
    )
    construct_grid_world_path(
        VecE2(x0,y0),
        t0,
        cellwidth,
        transition_time,
        transitions
        )
end
function construct_grid_world_path(
    x0::Float64,
    y0::Float64,
    t0::Float64,
    cellwidth::Float64,
    transition_time::Float64,
    transitions::Vector{Int}
    )
    construct_grid_world_path(
        x0, y0, t0, cellwidth, transition_time,
        map(i->GridTransition(i),transitions)
        )
end

end # End of module Paths


text:
@testset "build" begin
    mktempdir() do dir
        files = map(1:2) do i
            without_extension = "file$i"
            file = "$(without_extension).jl"
            content = pluto_notebook_content("""
                x = begin
                    sleep(3)
                    x = 3000 + $i
                end
                """)
            path = joinpath(dir, file)
            write(path, content)
            return file
        end
        build_notebooks(BuildOptions(dir, use_distributed=false))

        html_file = joinpath(dir, "file1.html")
        @test contains(read(html_file, String), "3001")

        html_file = joinpath(dir, "file2.html")
        @test contains(read(html_file, String), "3002")
    end
end

@testset "is_pluto_file" begin
    cd(mktempdir()) do
        nb_text = """
            ### A Pluto.jl notebook ###
            # v0.14.0
            """
        write("true.jl", nb_text)
        @test PlutoStaticHTML._is_pluto_file("true.jl")

        jl_text = """
            module Foo
            end # module
            """
        write("false.jl", jl_text)
        @test !PlutoStaticHTML._is_pluto_file("false.jl")
    end
end

@testset "invalid_notebook" begin
    try
        mktempdir() do dir
            content = pluto_notebook_content("@assert false")
            file = "file.jl"
            path = joinpath(dir, file)
            write(path, content)
            build_notebooks(BuildOptions(dir), [file], use_distributed=false)
        end
        error("Test should have failed")
    catch AssertionError
        @test true # Success.
    end
end

@testset "extract_previous_output" begin
    # function extract_previous_output(html::AbstractString)::String
    html = """
        lorem
        $(PlutoStaticHTML.BEGIN_IDENTIFIER)
        ipsum
        $(PlutoStaticHTML.END_IDENTIFIER)
        dolar
        """
    actual = PlutoStaticHTML.extract_previous_output(html)
    expected = """
        $(PlutoStaticHTML.BEGIN_IDENTIFIER)
        ipsum
        $(PlutoStaticHTML.END_IDENTIFIER)
        """
    @test strip(actual) == strip(expected)
end


text:
##############################################################################
##
## FixedEffect
##
## The categoricalarray may have pools that are never referred. Note that the pool does not appear in FixedEffect anyway.
##
##############################################################################

struct FixedEffect{R <: AbstractVector{<:Integer}, I <: AbstractVector{<: Real}}
	refs::R                 # refs must be between 0 and n
	interaction::I          # the continuous interaction
	n::Int                  # Number of potential values (= maximum(refs))
	function FixedEffect{R, I}(refs, interaction, n) where {R <: AbstractVector{<:Integer}, I <: AbstractVector{<: Real}}
		length(refs) == length(interaction) || error("refs and interaction don't have the same length")
		new(refs, interaction, n)
	end
end

function FixedEffect(args...; interaction::AbstractVector = uweights(length(args[1])))
	g = group(args...)
	FixedEffect{typeof(g.refs), typeof(interaction)}(g.refs, interaction, g.n)
end

function Base.show(io::IO, fe::FixedEffect)
	println(io, "Refs:        ", length(fe.refs), "-element ", typeof(fe.refs))
	println(io, "             ", Int.(fe.refs[1:min(5, length(fe.refs))]), "...")
	println(io, "Interaction: ", length(fe.interaction), "-element ", typeof(fe.interaction))
	println(io, "             ", fe.interaction[1:min(5, length(fe.interaction))], "...")
end

Base.length(fe::FixedEffect) = length(fe.refs)
Base.eltype(fe::FixedEffect) = eltype(I)

Base.getindex(fe::FixedEffect, esample::Colon) = fe
function Base.getindex(fe::FixedEffect{R, I}, esample::AbstractVector) where {R, I}
	if fe.interaction isa UnitWeights
		FixedEffect{R, I}(fe.refs[esample], uweights(eltype(I), sum(esample)),  fe.n)
	else
		FixedEffect{R, I}(fe.refs[esample], fe.interaction[esample],  fe.n)
	end
end

##############################################################################
##
## group combines multiple refs
## Missings have a ref of 0
##
##############################################################################

mutable struct GroupedArray{N} <: AbstractArray{UInt32, N}
	refs::Array{UInt32, N}   # refs must be between 0 and n. 0 means missing
	n::Int                   # Number of potential values (= maximum(refs))
end
group(xs::GroupedArray) = xs

function group(xs::AbstractArray)
	refs = Array{UInt32}(undef, size(xs))
	invpool = Dict{eltype(xs), UInt32}()
	n = 0
	has_missing = false
	@inbounds for i in eachindex(xs)
		x = xs[i]
		if x === missing
			refs[i] = 0
			has_missing = true
		else
			lbl = get(invpool, x, 0)
			if lbl !== 0
				refs[i] = lbl
			else
				n += 1
				refs[i] = n
				invpool[x] = n
			end
		end
	end
	return GroupedArray{ndims(xs)}(refs, n)
end

function group(args...)
	g1 = deepcopy(group(args[1]))
	for j = 2:length(args)
		gj = group(args[j])
		length(g1.refs) == length(gj.refs) || throw(DimensionError())
		combine!(g1, gj)
	end
	factorize!(g1)
end

function combine!(g1::GroupedArray, g2::GroupedArray)
	@inbounds for i in eachindex(g1.refs, g2.refs)
		# if previous one is missing or this one is missing, set to missing
		g1.refs[i] = (g1.refs[i] == 0 || g2.refs[i] == 0) ? 0 : g1.refs[i] + (g2.refs[i] - 1) * g1.n
	end
	g1.n = g1.n * g2.n
	return g1
end

function factorize!(x::GroupedArray{N}) where {N}
	refs = x.refs
	uu = sort!(unique(refs))
	has_missing = uu[1] == 0
	ngroups = length(uu) - has_missing
	dict = Dict{UInt32, UInt32}(zip(uu, UInt32(1-has_missing):UInt32(ngroups)))
	for i in eachindex(refs)
		refs[i] = dict[refs[i]]
	end
	GroupedArray{N}(refs, ngroups)
end

##############################################################################
##
## Find connected components
##
##############################################################################
# Return a vector of sets that contains the indices of each unique value
function refsrev(fe::FixedEffect)
	out = Vector{Int}[Int[] for _ in 1:fe.n]
	for i in eachindex(fe.refs)
		push!(out[fe.refs[i]], i)
	end
	return out
end

# Returns a vector of all components
# A component is a vector that, for each fixed effect,
# contains all the refs that are included in the component.
function components(fes::AbstractVector{<:FixedEffect})
	refs_vec = Vector{UInt32}[fe.refs for fe in fes]
	refsrev_vec = Vector{Vector{Int}}[refsrev(fe) for fe in fes]
	visited = falses(length(refs_vec[1]))
	out = Vector{Set{Int}}[]
	for i in eachindex(visited)
		if !visited[i]
			# obs not visited yet, so create new component
			component_vec = Set{UInt32}[Set{UInt32}() for _ in 1:length(refsrev_vec)]
			# visit all obs in the same components
			tovisit = Set{Int}(i)
			while !isempty(tovisit)
				for (component, refs, refsrev) in zip(component_vec, refs_vec, refsrev_vec)
					ref = refs[i]
					# if group is not in component yet
					if ref ∉ component
						# add group to the component
						push!(component, ref)
						# visit other observations in same group
						union!(tovisit, refsrev[ref])
					end
				end
				# mark obs as visited
				i = pop!(tovisit)
				visited[i] = true
			end
			push!(out, component_vec)
		end
	end
	return out
end

##############################################################################
##
## normalize! a vector of fixedeffect coefficients using connected components
##
##############################################################################

function normalize!(fecoefs::AbstractVector{<: Vector{<: Real}}, fes::AbstractVector{<:FixedEffect}; kwargs...)
	# The solution is generally not unique. Find connected components and scale accordingly
	idx = findall(fe -> isa(fe.interaction, UnitWeights), fes)
	length(idx) >= 2 && rescale!(view(fecoefs, idx), view(fes, idx))
	return fecoefs
end

function rescale!(fecoefs::AbstractVector{<: Vector{<: Real}}, fes::AbstractVector{<:FixedEffect})
	for component_vec in components(fes)
		m = 0.0
		# demean all fixed effects except the first
		for j in length(fecoefs):(-1):2
			fecoef, component = fecoefs[j], component_vec[j]
			mj = 0.0
			for k in component
				mj += fecoef[k]
			end
			mj = mj / length(component)
			for k in component
				fecoef[k] -= mj
			end
			m += mj
		end
		# rescale the first fixed effects
		fecoef, component = fecoefs[1], component_vec[1]
		for k in component
			fecoef[k] += m
		end
	end
end

function full(fecoefs::AbstractVector{<: Vector{<: Real}}, fes::AbstractVector{<:FixedEffect})
	[fecoef[fe.refs] for (fecoef, fe) in zip(fecoefs, fes)]
end