text:
function buildNTwoSiteVertical(A::fPEPS,
                               L::Environments,
                               R::Environments,
                               IEnvs,
                               row::Int,
                               col::Int,
                               ϕ::NamedTuple)::ITensor
    Ny, Nx   = size(A)
    up_row   = row + 1
    Nbelow   = L.I[row]
    if row > 1
        Nbelow *= IEnvs[:below][row - 1]
    end
    Nbelow *= ϕ[:Ab]
    Nbelow *= R.I[row]
    Nbelow *= dag(ϕ[:Ab]')
    Nabove = L.I[up_row]
    if up_row < Ny
        Nabove *= IEnvs[:above][end - up_row]
    end
    Nabove *= ϕ[:Aa]
    Nabove *= R.I[up_row]
    Nabove *= dag(ϕ[:Aa]')
    workingN  = ϕ[:sab]
    workingN *= spinI(firstind(A[up_row, col], "Site"); is_gpu=is_gpu(A))
    workingN *= spinI(firstind(A[row, col], "Site"); is_gpu=is_gpu(A))
    workingN *= Nbelow
    workingN *= Nabove
    AAinds    = inds(prime(ϕ[:sab]))
    @assert hasinds(inds(workingN), AAinds) "row $row\nIndsN: $(inds(workingN))\nIndsAA: $AAinds"
    @assert hasinds(AAinds, inds(workingN)) "row $row\nIndsN: $(inds(workingN))\nIndsAA: $AAinds"
    return workingN
end

function buildHIedgeTwoSiteVertical(A::fPEPS,
                                    E::Environments,
                                    row::Int,
                                    col::Int,
                                    side::Symbol,
                                    ϕ::NamedTuple)
    Ny, Nx = size(A)
    is_cu  = is_gpu(A)
    up_row = row + 1
    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 up_row+1:Ny
        IH_a *= E.H[work_row] * A[work_row, col]
        IH_a *= dag(prime(A[work_row, col], "Link"))
    end
    IH_b *= E.H[row]
    IH_a *= E.H[up_row]
    IH_b *= ϕ[:Ab]
    IH_b *= dag(ϕ[:Ab]')
    IH_a *= ϕ[:Aa]
    IH_a *= dag(ϕ[:Aa]')

    IH = ϕ[:sab]
    IH *= spinI(firstind(A[row, col], "Site"); is_gpu=is_cu)
    IH *= spinI(firstind(A[up_row, col], "Site"); is_gpu=is_cu)
    IH *= IH_b
    IH *= IH_a
    AAinds    = inds(prime(ϕ[:sab]))
    @assert hasinds(inds(IH), AAinds)
    @assert hasinds(AAinds, inds(IH))
    return (IH,)
end

function buildHIsTwoSiteVertical(A::fPEPS,
                                 L::Environments,
                                 R::Environments,
                                 row::Int,
                                 col::Int,
                                 ϕ::NamedTuple)
    Ny, Nx = size(A)
    is_cu  = is_gpu(A)
    up_row = row + 1
    @timeit "buildHIedge" begin
        if col == 1
            return buildHIedgeTwoSiteVertical(A, R, row, col, :left, ϕ)
        end
        if col == Nx
            return buildHIedgeTwoSiteVertical(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(up_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_b *= L.H[row]
        HLI_a *= L.H[up_row]
        HLI_b *= ϕ[:Ab]
        HLI_a *= ϕ[:Aa]
        HLI_b *= R.I[row]
        HLI_a *= R.I[up_row]
        HLI_b *= dag(ϕ[:Ab]')
        HLI_a *= dag(ϕ[:Aa]')
        IHR_b *= R.H[row]
        IHR_a *= R.H[up_row]
        IHR_b *= ϕ[:Ab]
        IHR_a *= ϕ[:Aa]
        IHR_b *= L.I[row]
        IHR_a *= L.I[up_row]
        IHR_b *= dag(ϕ[:Ab]')
        IHR_a *= dag(ϕ[:Aa]')
        HLI *= ϕ[:sab]
        IHR *= ϕ[:sab]
        HLI *= spinI(firstind(A[row, col], "Site"); is_gpu=is_cu)
        HLI *= spinI(firstind(A[up_row, col], "Site"); is_gpu=is_cu)
        IHR *= spinI(firstind(A[row, col], "Site"); is_gpu=is_cu)
        IHR *= spinI(firstind(A[up_row, col], "Site"); is_gpu=is_cu)
        HLI *= HLI_a
        IHR *= IHR_a
        HLI *= HLI_b
        IHR *= IHR_b
    end
    AAinds    = inds(prime(ϕ[:sab]))
    @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 verticalTermsTwoSiteVertical(A::fPEPS,
                                      L::Environments,
                                      R::Environments,
                                      AI,
                                      AV,
                                      H,
                                      row::Int,
                                      col::Int,
                                      ϕ::NamedTuple)::Vector{ITensor}
    Ny, Nx = size(A)
    is_cu  = is_gpu(A)
    vTerms = Vector{ITensor}(undef, length(H))
    AAinds = inds(prime(ϕ[:sab]))
    dummy  = is_cu ? cuITensor(1.0) : ITensor(1.0)
    up_row = row + 1
    @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 || up_row < op_row_a
            local Va, Vb
            if op_row_a > up_row
                Va = AV[:above][opcode][end - up_row]
                Vb = row > 1 ? AI[:below][row - 1] : dummy
            elseif op_row_b < row
                Vb = AV[:below][opcode][row - op_row_b]
                Va = up_row < Ny ? AI[:above][end - up_row] : dummy
            end
            Vb *= L.I[row]
            Vb *= ϕ[:Ab]
            Vb *= R.I[row]
            Vb *= dag(ϕ[:Ab]')
            Va *= L.I[up_row]
            Va *= ϕ[:Aa]
            Va *= R.I[up_row]
            Va *= dag(ϕ[:Aa]')
            thisVert = ϕ[:sab]
            thisVert *= spinI(firstind(A[row, col], "Site"); is_gpu=is_cu)
            thisVert *= spinI(firstind(A[up_row, col], "Site"); is_gpu=is_cu)
            thisVert *= Vb
            thisVert *= Va
        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')
            AIH *= L.I[op_row_b]
            AIH *= ϕ[:Aa]
            AIH *= R.I[op_row_b]
            AIH *= dag(ϕ[:Aa]')
            AIL *= L.I[op_row_a]
            AIL *= ϕ[:Ab]
            AIL *= R.I[op_row_a]
            AIL *= dag(ϕ[:Ab]')
            thisVert = ϕ[:sab]
            thisVert *= op_a
            thisVert *= op_b
            thisVert *= AIL
            thisVert *= AIH
        elseif row == op_row_b
            low_row  = op_row_a - 1
            high_row = op_row_b + 1
            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')
            AIL *= L.I[op_row_a]
            AIL *= A[op_row_a, col] * op_a
            AIL *= R.I[op_row_a]
            AIL *= dag(A[op_row_a, col])'
            AIL *= L.I[op_row_b]
            AIL *= ϕ[:Ab]
            AIL *= R.I[op_row_b]
            AIL *= dag(ϕ[:Ab]')
            AIH *= L.I[up_row]
            AIH *= ϕ[:Aa]
            AIH *= R.I[up_row]
            AIH *= dag(ϕ[:Aa]')
            thisVert *= ϕ[:sab]
            thisVert *= op_b
            thisVert *= spinI(firstind(A[up_row, col], "Site"); is_gpu=is_cu)
            thisVert *= AIL
            thisVert *= AIH
        elseif up_row == op_row_a
            low_row  = op_row_a - 2
            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')
            AIH *= L.I[op_row_b]
            AIH *= A[op_row_b, col] * op_b
            AIH *= R.I[op_row_b]
            AIH *= dag(A[op_row_b, col])'
            AIH *= L.I[op_row_a]
            AIH *= ϕ[:Aa]
            AIH *= R.I[op_row_a]
            AIH *= dag(ϕ[:Aa]')
            AIL *= L.I[row]
            AIL *= ϕ[:Ab]
            AIL *= R.I[row]
            AIL *= dag(ϕ[:Ab]')
            thisVert  = ϕ[:sab]
            thisVert *= op_a
            thisVert *= spinI(firstind(A[row, col], "Site"); is_gpu=is_cu)
            thisVert *= AIL
            thisVert *= AIH
        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)
            vTerms[opcode] = thisVert
        end
    end
    return vTerms
end

function fieldTermsTwoSiteVertical(A::fPEPS,
                                   L::Environments,
                                   R::Environments,
                                   AI,
                                   AF,
                                   H,
                                   row::Int,
                                   col::Int,
                                   ϕ::NamedTuple)::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(ϕ[:sab]))
    up_row = row + 1
    @inbounds for opcode in 1:length(H)
        thisField = dummy
        op_row    = H[opcode].sites[1][1]
        if op_row != row && up_row != op_row
            local Fa, Fb
            if op_row > up_row
                Fa = AF[:above][opcode][end - up_row]
                Fb = row > 1 ? AI[:below][row - 1] : dummy
            else
                Fb = AF[:below][opcode][row - 1]
                Fa = up_row < Ny ? AI[:above][end - up_row] : dummy
            end
            Fb *= L.I[row]
            Fb *= ϕ[:Ab]
            Fb *= R.I[row]
            Fb *= dag(ϕ[:Ab]')
            Fa *= L.I[up_row]
            Fa *= ϕ[:Aa]
            Fa *= R.I[up_row]
            Fa *= dag(ϕ[:Aa]')
            thisField = ϕ[:sab]
            thisField *= spinI(firstind(A[row, col], "Site"); is_gpu=is_cu)
            thisField *= spinI(firstind(A[up_row, col], "Site"); is_gpu=is_cu)
            thisField *= Fa
            thisField *= Fb
        else
            low_row  = row - 1
            high_row = up_row
            AIL = low_row > 0   ? AI[:below][low_row]        : dummy
            AIH = high_row < Ny ? AI[:above][end - high_row] : dummy
            thisField  = AIL
            AIL *= L.I[row]
            AIH *= L.I[up_row]
            AIL *= ϕ[:Ab]
            AIH *= ϕ[:Aa]
            AIL *= R.I[row]
            AIH *= R.I[up_row]
            AIL *= dag(ϕ[:Ab]')
            AIH *= dag(ϕ[:Aa]')
            thisField = ϕ[:sab]
            sA = firstind(A[row, col], "Site")
            sB = firstind(A[up_row, col], "Site")
            if op_row == row
                op = copy(H[opcode].ops[1])
                op = replaceind!(op, H[opcode].site_ind, sA)
                op = replaceind!(op, H[opcode].site_ind', sA')
                thisField *= op
                thisField *= spinI(sB; is_gpu=is_cu)
            elseif op_row == up_row
                op = copy(H[opcode].ops[1])
                op = replaceind!(op, H[opcode].site_ind, sB)
                op = replaceind!(op, H[opcode].site_ind', sB')
                thisField *= op
                thisField *= spinI(sA; is_gpu=is_cu)
            end
            thisField *= AIH
            thisField *= AIL
        end
        @assert hasinds(inds(thisField), AAinds) "row $row\ntf: $(inds(thisField))\naa: $AAinds"
        @assert hasinds(AAinds, inds(thisField)) "row $row\ntf: $(inds(thisField))\naa: $AAinds"
        fTerms[opcode] = thisField
    end
    return fTerms
end

function connectLeftTermsTwoSiteVertical(A::fPEPS,
                                         L::Environments,
                                         R::Environments,
                                         AI, AL, H,
                                         row::Int, col::Int,
                                         ϕ::NamedTuple)::Vector{ITensor}
    Ny, Nx = size(A)
    up_row = row + 1
    is_cu  = is_gpu(A)
    lTerms = Vector{ITensor}(undef, length(H))
    AAinds = inds(prime(ϕ[:sab]))
    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 && op_row_b != up_row
            if op_row_b > up_row
                La = AL[:above][opcode][end - up_row]
                Lb = row > 1 ? AL[:below][opcode][row - 1] : dummy
            else
                Lb = AL[:below][opcode][row - 1]
                La = up_row < Ny ? AL[:above][opcode][end - up_row] : dummy
            end
            Lb *= L.InProgress[row, opcode]
            La *= L.InProgress[up_row, opcode]
            Lb *= ϕ[:Ab]
            La *= ϕ[:Aa]
            Lb *= R.I[row]
            La *= R.I[up_row]
            Lb *= dag(ϕ[:Ab]')
            La *= dag(ϕ[:Aa]')
            thisHori  = ϕ[:sab]
            thisHori *= spinI(firstind(A[row, col], "Site"); is_gpu=is_cu)
            thisHori *= spinI(firstind(A[up_row, col], "Site"); is_gpu=is_cu)
            thisHori *= La
            thisHori *= Lb
        else
            low_row  = (op_row_b <= row)    ? op_row_b - 1 : row - 1
            high_row = (op_row_b >= up_row) ? op_row_b + 1 : up_row + 1
            THB  = low_row >= 1   ? AL[:below][opcode][low_row]            : dummy
            THB *= R.I[row]
            THB *= ϕ[:Ab]
            THB *= L.InProgress[row, opcode]
            THB *= dag(ϕ[:Ab]')
            thisHori = ϕ[:sab]
            if row == op_row_b
                thisHori *= op_b
                thisHori *= spinI(firstind(A[up_row, col], "Site"); is_gpu=is_cu)
            elseif up_row == op_row_b
                thisHori *= op_b
                thisHori *= spinI(firstind(A[row, col], "Site"); is_gpu=is_cu)
            end
            thisHori *= THB
            THA  = high_row <= Ny ? AL[:above][opcode][end - high_row + 1] : dummy
            THA *= R.I[up_row]
            THA *= ϕ[:Aa]
            THA *= L.InProgress[up_row, opcode]
            THA *= dag(ϕ[:Aa]')
            thisHori *= THA
        end
        @assert hasinds(inds(thisHori), AAinds) "row $row\ntH: $(inds(thisHori))\nAA: $AAinds"
        @assert hasinds(AAinds, inds(thisHori)) "row $row\ntH: $(inds(thisHori))\nAA: $AAinds"
        lTerms[opcode] = thisHori
    end
    return lTerms
end

function connectRightTermsTwoSiteVertical(A::fPEPS,
                                          L::Environments,
                                          R::Environments,
                                          AI, AR, H,
                                          row::Int, col::Int,
                                          ϕ::NamedTuple)::Vector{ITensor}
    Ny, Nx = size(A)
    is_cu  = is_gpu(A)
    up_row = row + 1
    rTerms = Vector{ITensor}(undef, length(H))
    AAinds = inds(prime(ϕ[:sab]))
    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 && op_row_a != up_row
            if op_row_a > up_row
                Ra = AR[:above][opcode][end - up_row]
                Rb = row > 1 ? AR[:below][opcode][row - 1] : dummy
            else
                Rb = AR[:below][opcode][row - 1]
                Ra = up_row < Ny ? AR[:above][opcode][end - up_row] : dummy
            end
            Rb *= R.InProgress[row, opcode]
            Ra *= R.InProgress[up_row, opcode]
            Rb *= ϕ[:Ab]
            Ra *= ϕ[:Aa]
            Rb *= L.I[row]
            Ra *= L.I[up_row]
            Rb *= dag(ϕ[:Ab]')
            Ra *= dag(ϕ[:Aa]')
            thisHori *= ϕ[:sab]
            thisHori *= spinI(firstind(A[row, col], "Site"); is_gpu=is_cu)
            thisHori *= spinI(firstind(A[up_row, col], "Site"); is_gpu=is_cu)
            thisHori *= Ra
            thisHori *= Rb
        else
            low_row  = (op_row_a <= row)    ? op_row_a - 1 : row - 1;
            high_row = (op_row_a >= up_row) ? op_row_a + 1 : up_row + 1;
            THB  = low_row >= 1   ? AR[:below][opcode][low_row]            : dummy
            THA  = high_row <= Ny ? AR[:above][opcode][end - high_row + 1] : dummy
            THB *= L.I[row]
            THA *= L.I[up_row]
            THB *= ϕ[:Ab]
            THA *= ϕ[:Aa]
            THB *= R.InProgress[row, opcode]
            THA *= R.InProgress[up_row, opcode]
            THB *= dag(ϕ[:Ab]')
            THA *= dag(ϕ[:Aa]')
            thisHori = ϕ[:sab]
            thisHori *= op_a
            if op_row_a == row
                thisHori *= spinI(firstind(A[up_row, col], "Site"); is_gpu=is_cu)
            elseif op_row_a == up_row
                thisHori *= spinI(firstind(A[row, col], "Site"); is_gpu=is_cu)
            end
            thisHori *= THB
            thisHori *= THA
        end
        @assert hasinds(inds(thisHori), AAinds)
        @assert hasinds(AAinds, inds(thisHori))
        rTerms[opcode] = thisHori
    end
    return rTerms
end

function buildLocalHTwoSiteVertical(A::fPEPS,
                                    L::Environments, R::Environments,
                                    AncEnvs, H,
                                    row::Int, col::Int,
                                    ϕ::NamedTuple;
                                    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   = buildNTwoSiteVertical(A, L, R, AncEnvs[:I], row, col, ϕ)
    end
    expϕ = ϕ[:Ab] * ϕ[:sab] * ϕ[:Aa]
    den = scalar(collect(N * dag(ϕ[:sab])'))
    local left_H_terms, right_H_terms
    if col > 1
        left_H_terms = getDirectional(vcat(H[:, col - 1]...), Horizontal)
        term_count  += length(left_H_terms)
    end
    if col < Nx
        right_H_terms = getDirectional(vcat(H[:, col]...), Horizontal)
        term_count   += length(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 = buildHIsTwoSiteVertical(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(ϕ[:sab])'))
            end
        end
    end
    @debug "\t\tBuilding vertical H terms row $row col $col"
    @timeit "build vertical terms" begin
        vTs = verticalTermsTwoSiteVertical(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 ---")
            for vT in vTs
                println(scalar(vT * dag(ϕ[:sab])'))
            end
        end
    end
    @debug "\t\tBuilding field H terms row $row col $col"
    @timeit "build field terms" begin
        fTs = fieldTermsTwoSiteVertical(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 ---")
            for fT in fTs
                println(scalar(fT * dag(ϕ[:sab])'))
            end
        end
    end
    if col > 1
        @debug "\t\tBuilding left H terms row $row col $col"
        @timeit "build left terms" begin
            lTs = connectLeftTermsTwoSiteVertical(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 ---")
                for lT in lTs
                    println(scalar(lT * dag(ϕ[:sab])'))
                end
            end
        end
        @debug "\t\tBuilt left terms"
    end
    if col < Nx
        @debug "\t\tBuilding right H terms row $row col $col"
        @timeit "build right terms" begin
            rTs = connectRightTermsTwoSiteVertical(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 ---")
                for rT in rTs
                    println(scalar(rT * dag(ϕ[:sab])'))
                end
            end
        end
        @debug "\t\tBuilt right terms"
    end
    return Hs, N
end

struct ITensorMapTwoSiteVertical
  A::fPEPS
  H::Matrix{Vector{Operator}}
  L::Environments
  R::Environments
  AncEnvs::NamedTuple
  row::Int
  col::Int
  ϕa::ITensor
  ϕb::ITensor
end
Base.eltype(M::ITensorMapTwoSiteVertical) = eltype(M.A)
function Base.size(M::ITensorMapTwoSiteVertical)
    vci  = dim(commonindex(M.A[M.row, M.col], M.A[M.row+1, M.col]))
    saci = dim(firstind(M.A[M.row, M.col], "Site"))
    sbci = dim(firstind(M.A[M.row+1, M.col], "Site"))
    return 2*vci*saci*sbci
end

function (M::ITensorMapTwoSiteVertical)(v)
    ϕ = (Aa=M.ϕa, Ab=M.ϕb, sab=v)
    Hs, N  = buildLocalHTwoSiteVertical(M.A, M.L, M.R, M.AncEnvs, M.H, M.row, M.col, ϕ)
    localH = sum(Hs)
    return noprime(localH)
end

function optimizeLocalHTwoSiteVertical(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)
    Rainds        = IndexSet(firstind(A[row, col], "Site"), commonindex(A[row, col], A[row+1, col]))
    Lainds        = uniqueinds(inds(A[row, col]), IndexSet(firstind(A[row, col], "Site"), commonindex(A[row, col], A[row+1, col])))
    Qa, Ra, qra   = qr(A[row, col], Lainds; tags="sa,Link,r$row,c$col")
    Rbinds        = IndexSet(firstind(A[row+1, col], "Site"), commonindex(A[row, col], A[row+1, col]))
    Lbinds        = uniqueinds(inds(A[row+1, col]), IndexSet(firstind(A[row+1, col], "Site"), commonindex(A[row, col], A[row+1, col])))
    Qb, Rb, qrb   = qr(A[row+1, col], Lbinds; tags="sb,Link,r$(row+1),c$col")
    ϕa            = Qa
    ϕb            = Qb
    ϕs            = Ra * Rb
    ϕ             = (Ab=ϕa, Aa=ϕb, sab=ϕs)
    @debug "\tBuilding H for col $col row $row"
    @timeit "build H" begin
        Hs, N = buildLocalHTwoSiteVertical(A, L, R, AncEnvs, H, row, col, ϕ)
    end
    expϕ = ϕ[:Ab] * ϕ[:sab] * ϕ[:Aa]
    initial_N = real(scalar(collect(N * dag(ϕ[:sab]'))))
    @timeit "sum H terms" begin
        localH = sum(Hs)
    end
    initial_E = real(scalar(collect(deepcopy(localH) * dag(ϕ[:sab])')))
    @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   = ITensorMapTwoSiteVertical(A, H, L, R, AncEnvs, row, col, ϕ[:Aa], ϕ[:Ab])
    λ, new_ϕ = davidson(mapper, ϕ[:sab]; maxiter=2, kwargs...)
    new_E    = λ

    # split the result
    Lis      = IndexSet(firstind(A[row, col], "Site"), commonindex(ϕa, new_ϕ))
    U, S, Vt = svd(new_ϕ, Lis; tags=tags(commonindex(A[row, col], A[row+1, col])), kwargs...)
    A[row, col]   = ϕa * U
    A[row+1, col] = ϕb * S * Vt
    #=N        = buildNTwoSite(A, L, R, AncEnvs[:I], row, col, new_ϕ)
    new_N    = real(scalar(collect(N * dag(new_ϕ)')))
    @info "Optimized energy at row $row col $col : $(new_E/(new_N*Nx*Ny)) and norm : $new_N"=#
    @timeit "restore intraColumnGauge" begin
        if row <= Ny - 1
            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    = A[row, col]*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 = copy(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
        end
    end
    return A, AncEnvs
end


text:
using Query
using DataFrames

df = DataFrame(name=["John", "Sally", "Kirk"], age=[23., 42., 59.], children=[3,2,2])

x = @from i in df begin
    @group i.name by i.children
    @collect
end

println(x)

x = @from i in df begin
    @group i by i.children
    @collect
end

println(x)

ext: jl


text:
function readinput(filename)
    return filter(!isempty, split(read(filename, String),"\n"))
end

function gentiles(lines)
    tiles = Dict()
    tileid = 0
    for line in lines
        if startswith(line, "Tile ")
            tileid = line[6:end-1]
        else
            push!(tiles, tileid => append!((get!(tiles, tileid, [])),[line]))
        end
    end
    return tiles
end

function getedges(tile)
    prefix :: Array{Any,1} = [tile[1], join(x[end] for x in tile), tile[end], join(x[1] for x in tile)]
    append!(prefix, [tile])
    return prefix
end

function fliptile(tile)
    return collect(reverse(x) for x in tile)
end

function rotate(tile, num)
    newtile = tile
    for _ in 1:num
        reversetile = reverse(newtile)
        newtile = []
        for i in 1:length(reversetile[1])
            slice = ""
            for r in reversetile
                slice *= r[i]
            end
            push!(newtile, slice)
        end
    end
    return newtile
end

function genedges(tiles)
    edges = Dict()
    for (tileid, tile) in tiles
        e = []
        append!(e, [getedges(tile)])
        append!(e, [getedges(rotate(tile, 1))])
        append!(e, [getedges(rotate(tile, 2))])
        append!(e, [getedges(rotate(tile, 3))])
        append!(e, [getedges(fliptile(tile))])
        append!(e, [getedges(fliptile(rotate(tile, 1)))])
        append!(e, [getedges(fliptile(rotate(tile, 2)))])
        append!(e, [getedges(fliptile(rotate(tile, 3)))])
        push!(edges, tileid => e)
    end
    return edges
end

function geteligibles(curr, eligibles)
    return collect([(a,b) for (a,b) in eligibles if !(a in [c[1] for c in curr])])
end

function multcorner(nextedges, sidesize, totalsize)
    a = parse(Int,nextedges[1][1])
    b = parse(Int,nextedges[sidesize][1])
    c = parse(Int,nextedges[totalsize][1])
    d = parse(Int,nextedges[totalsize - sidesize + 1][1])
    return a * b * c * d
end

function mountimage(lines)
    alltiles = gentiles(lines)
    alledges = genedges(alltiles)
    sidesize = Int(sqrt(length(alltiles)))
    totalsize = length(alltiles)
    northedges, eastedges, southedges, westedges = Dict(), Dict(), Dict(), Dict()

    for (tileid, edges) in alledges
        for edge in edges
            e = pop!(northedges, edge[1], [])
            push!(e, (tileid, edge))
            push!(northedges, edge[1] => e)

            e = pop!(eastedges, edge[2], [])
            push!(e, (tileid, edge))
            push!(eastedges, edge[2] => e)

            e = pop!(southedges, edge[3], [])
            push!(e, (tileid, edge))
            push!(southedges, edge[3] => e)

            e = pop!(westedges, edge[4], [])
            push!(e, (tileid, edge))
            push!(westedges, edge[4] => e)
        end
    end

    for (itileid, iedges) in alledges
        for iedge in iedges
            eligibles = geteligibles([(itileid, iedge)], get(westedges, iedge[2], []))
            stack = [([(itileid, iedge)], eligibles)]

            while !isempty(stack)
                (currentedges, eligibles) = pop!(stack)

                if isempty(eligibles)
                    continue
                end

                currenteligible = pop!(eligibles)
                push!(stack, (currentedges, eligibles))
                nextedges = deepcopy(currentedges)
                push!(nextedges, currenteligible)

                if length(nextedges) == totalsize
                    return (multcorner(nextedges, sidesize, totalsize), nextedges)
                end

                nexteligibles = []

                if length(nextedges) < sidesize
                    nexteligibles = geteligibles(nextedges, get(westedges, currenteligible[2][2], []))
                else
                    above = nextedges[length(nextedges)-sidesize+1]
                    nexteligibles = geteligibles(nextedges, get(northedges, above[2][3], []))
                    if length(nextedges) % sidesize > 1
                        prev = nextedges[end]
                        nexteligibles = collect([n for n in nexteligibles if n[2][4] == prev[2][2]])
                    end
                end

                if !isempty(nexteligibles)
                    push!(stack, (nextedges, nexteligibles))
                end
            end
        end
    end
    return nothing
end

function combine(tiles)
    tilessize = length(tiles)
    side = Int(sqrt(tilessize))
    lines = []
    for i in collect(1:side:tilessize)
        tls = tiles[i:i+side-1]
        for idx in 1:length(tls[1])
            push!(lines, join(x[idx] for x in tls))
        end
    end
    return lines
end

function readmonster(filename)
    lines = readlines(filename)
    head, offset = -1, -1
    mid, low = [], []
    for (i, line) in enumerate(lines)
        for (j, char) in enumerate(line)
            if i == 1 && char == '#'
                head = j
            elseif i == 2 && char == '#'
                append!(mid, j - head)
            elseif i == 3 && char == '#'
                append!(low, j - head)
            elseif i != 1 && j > head && j - head > offset
                offset = j - head
            end
        end
    end
    return (head, offset, mid, low)
end

function findmonsters(image, monster="d20.monster")
    side = length(image[1])
    (head, offset, mid, low) = readmonster(monster)
    count = 0
    for i in 19:side
        for j in 1:length(image)-1
            if image[j][i] == '#'
                if i < side-offset && all([image[j+1][i + k] == '#' for k in mid])
                    if  i < side-offset && all([image[j+2][i + k] == '#' for k in low])
                        count += 1
                    end
                end
            end
        end
    end
    return count
end

function monsterhunt(image; showimage=false)
    bordertiles = [i[2][5] for i in image]
    nobordertiles = [[i[2:end-1] for i in j[2:end-1]] for j in bordertiles]

    bordertiles = combine(bordertiles)
    nobordertiles = combine(nobordertiles)
    if showimage
        println("------BORDER------")
        println(join(bordertiles, '\n'))
        println("-----NOBORDER-----")
        println(join(rotate(bordertiles, 1),'\n'))
        println("------------------")
    end
    total = findmonsters(nobordertiles)
    total += findmonsters(rotate(nobordertiles, 2))
    total += findmonsters(rotate(nobordertiles, 3))
    total += findmonsters(rotate(nobordertiles, 4))
    total += findmonsters(fliptile(nobordertiles))
    total += findmonsters(fliptile(rotate(nobordertiles, 2)))
    total += findmonsters(fliptile(rotate(nobordertiles, 3)))
    total += findmonsters(fliptile(rotate(nobordertiles, 4)))

    return (count("#",join(nobordertiles)) - total * 15)
end

lines = readinput("d20.in")
P1,image = mountimage(lines)
println("P1: ", P1)
P2 = monsterhunt(image)
println("P2: ", P2)

max_stars_repo_name: eduellery/adventofcode

max_issues_repo_name: eduellery/adventofcode

max_forks_repo_name: eduellery/adventofcode


text:
##### Beginning of file

module PredictMDFull # Begin module PredictMDFull

__precompile__(true)

include(joinpath("registry_url_list.jl"))
include(joinpath("package_directory.jl"))
include(joinpath("version.jl"))

include(joinpath("package_list.jl"))

include(joinpath("import_required_packages.jl"))

include(joinpath("import_all.jl"))

include(joinpath("welcome.jl"))
include(joinpath("init.jl"))

end # End module PredictMDFull

##### End of file


text:
using Pope
using Pope.LJH

@testset "ljhutil" begin
    noipath = "ReferenceMicrocalFiles/ljh/20150707_C_chan13.noi"
    a,b,c=LJH.dir_base_ext(noipath)
    @test a=="ReferenceMicrocalFiles/ljh"
    @test b=="20150707_C"
    @test c==".noi"

    ljhpath = "ReferenceMicrocalFiles/ljh/20150707_D_chan13.ljh"
    a,b,c=LJH.dir_base_ext(ljhpath)
    @test a=="ReferenceMicrocalFiles/ljh"
    @test b=="20150707_D"
    @test c==".ljh"

    @test LJH.channel(ljhpath)==13
    @test LJH.channel(noipath)==13

    @test LJH.fnames(ljhpath,[1,2,3,4])==["ReferenceMicrocalFiles/ljh/20150707_D_chan1.ljh", "ReferenceMicrocalFiles/ljh/20150707_D_chan2.ljh", "ReferenceMicrocalFiles/ljh/20150707_D_chan3.ljh", "ReferenceMicrocalFiles/ljh/20150707_D_chan4.ljh"]

    localdirljhpath = "20150707_D_chan13.ljh"
    a,b,c = LJH.dir_base_ext(localdirljhpath)
    @test a=="."
    @test b=="20150707_D"
    @test c==".ljh"

    dir = "artifacts"
    isdir(dir) || mkdir(dir)
    a,b,c = LJH.dir_base_ext("artifacts")
    @test a == "artifacts"
    @test b == "artifacts"

    a,b,c = LJH.dir_base_ext("artifacts/")
    @test a == "artifacts"
    @test b == "artifacts"

dir = "artifacts/ljh_test1"
isdir(dir) || mkdir(dir)
fnames = LJH.fnames(joinpath(dir,"ljh_test1"),1:2:480)
for fname in fnames
    touch(fname)
end
ljhdict = LJH.allchannels(first(fnames))
ljhdict2 = LJH.allchannels(dir)
@test ljhdict == ljhdict2
@test collect(keys(ljhdict))==collect(1:2:480)
@test collect(values(ljhdict))==fnames

dir_noi = "artifacts/ljh_test2"
isdir(dir_noi) || mkdir(dir_noi)
fnames_noi = LJH.fnames(joinpath(dir_noi,"ljh_test2.noi"),1:2:480)
for fname in fnames_noi
    touch(fname)
end
ljhdict_noi = LJH.allchannels(first(fnames_noi))
@test collect(keys(ljhdict_noi))==collect(1:2:480)
@test collect(values(ljhdict_noi))==fnames_noi
@test ljhdict_noi != ljhdict

outputname(x) = Pope.outputname(x,"annotation")
@test outputname("abc_chan1.ljh") == outputname("abc_chan1.noi") == outputname("abc")
@test Pope.outputname("abc_chan1.ljh","model") == "./abc_model.hdf5"
@test Pope.outputname("abc_chan1.ljh","model","pdf") == "./abc_model.pdf"
@test Pope.outputname("abc_chan1.ljh","model",".pdf") == "./abc_model.pdf"
end


text:
function get_optimum end
function get_optimum_last end

############################################################
## Helpers for pre-allocating objects
############################################################

function proto_policy(endo, exo, value, agg_state, params, hh::Household)
  mc = MarkovChain(exo)

  state = merge(endo.grid[1], exo.grid[1])
  𝔼V = extrapolated_𝔼V(endo, BSpline(Linear()), value, exo, 1, 𝔼(hh))
  @unpack pol, pol_full = get_optimum(state, agg_state, 𝔼V, params, endo, hh::Household)

  (proto_pol=pol, proto_pol_full=pol_full)
end

############################################################
## Iterating the bellman operator etc
############################################################

function iterate_bellman!(value_new, value_old, policy, policies_full, endo, exo, converged, agg_state, params, hh::Household)
  mc = MarkovChain(exo)
  n = length(value_new)
  prog = Progress(n, desc="Iterating", offset=1, dt=1)
  for (i_exo, exo_state) in enumerate(mc.state_values)
    # Create interpolated expected value function
    itp_scheme = BSpline(Cubic(Line(OnGrid())))
    𝔼V = extrapolated_𝔼V(endo, itp_scheme, value_old, exo, i_exo, 𝔼(hh))

    for (i_endo, endo_state) in enumerate(endo.grid)
      states = merge(endo_state, exo_state)
      ProgressMeter.next!(prog)
      @unpack pol, pol_full, val, conv = get_optimum(states, agg_state, 𝔼V, params, endo, hh)

      policy[i_endo, i_exo]    = pol
      policies_full[i_endo, i_exo] = pol_full
      value_new[i_endo, i_exo] = val
      converged[i_endo, i_exo] = conv
    end
  end
end

############################################################
## Finite Horizon (OLG)
############################################################

function solve_bellman(endo, exo, aggregate_state, params, hh::Household, t_grid)
  containers = initialize_values_policies(endo, exo, aggregate_state, params, hh, t_grid)
  @unpack value, policy, policies_full, converged = containers

  solve_bellman_T!(value, policy, policies_full, endo, exo, converged, aggregate_state, params, hh, t_grid)

  (val = value, policy = policy, policies_full=StructArray(policies_full), converged=converged)
end

function initialize_values_policies(endo, exo, aggregate_state, params, hh::Household, t_grid)
  container_size = (length(endo), length(exo), length(t_grid))

  value = zeros(container_size)

  @unpack proto_pol, proto_pol_full = proto_policy(endo, exo, value[:,:,1], aggregate_state, params, hh)

  policy = fill(proto_pol, container_size)
  policies_full = fill(proto_pol_full, container_size)
  converged = trues(container_size)

  (value=value, policy=policy, policies_full=policies_full, converged=converged)
end

function solve_bellman_T!(value, policy, policies_full, endo, exo, converged, aggregate_state, params, hh::Household, t_grid)

  # Solving the last period (special-cased for bequests, etc)
  iterate_bellman_last!(@view(value[:,:,end]),
                   @view(policy[:,:,end]),
                   @view(policies_full[:,:,end]),
                   endo, exo,
                   @view(converged[:,:,end]), aggregate_state, params, hh)

  # Solving all other periods
  @showprogress "Bellman(T): VFI" for i_t in reverse(1:length(t_grid)-1)
    iterate_bellman!(@view(value[:,:,i_t]),
                     @view(value[:,:,i_t+1]),
                     @view(policy[:,:,i_t]),
                     @view(policies_full[:,:,i_t]),
                     endo, exo,
                     @view(converged[:,:,i_t]), aggregate_state, params, hh)
  end

  number_conv = sum(converged)

  length(converged) == number_conv || @warn "Bellman didn't converge at $(round((1-number_conv / length(converged)) * 100, digits=4))% ($(length(converged) - number_conv) states)"
end

# For OLG: solve last period
function iterate_bellman_last!(value, policy, policies_full, endo, exo, converged, agg_state, params, hh::Household)
  mc = MarkovChain(exo)
  n = length(value)
  prog = Progress(n, desc="Iterating", offset=1, dt=1)
  for (i_exo, exo_state) in enumerate(mc.state_values)
    # Create interpolated expected value function

    for (i_endo, endo_state) in enumerate(endo.grid)
      states = merge(endo_state, exo_state)
      ProgressMeter.next!(prog)
      @unpack pol, pol_full, val, conv = get_optimum_last(states, agg_state, params, endo, hh)

      policy[i_endo, i_exo]    = pol
      policies_full[i_endo, i_exo] = pol_full
      value[i_endo, i_exo] = val
      converged[i_endo, i_exo] = conv
    end
  end
end

############################################################
## Infinite horizon
############################################################

function solve_bellman(endo, exo, aggregate_state, params, hh::Household; maxiter=200, rtol=eps()^0.4)
  containers = initialize_values_policies(endo, exo, aggregate_state, params, hh)
  @unpack value_old, value_new, policy, policies_full, converged = containers

  solve_bellman!(value_old, value_new, policy, policies_full, endo, exo, converged, aggregate_state, params, hh::Household; maxiter=maxiter, rtol=rtol)

  (val = value_new, policy = policy, policies_full=StructArray(policies_full), converged=converged)
end

function initialize_values_policies(endo, exo, aggregate_state, params, hh::Household)
  container_size = (length(endo), length(exo))

  value_old = zeros(container_size)
  value_new = zeros(container_size)

  @unpack proto_pol, proto_pol_full = proto_policy(endo, exo, value_new, aggregate_state, params, hh)

  policy = fill(proto_pol, container_size)
  policies_full = fill(proto_pol_full, container_size)
  converged = trues(container_size)

  (value_old=value_old, value_new=value_new, policy=policy, policies_full=policies_full, converged=converged)
end

function solve_bellman!(value_old, value_new, policy, policies_full, endo, exo, converged, aggregate_state, params, hh::Household; maxiter, rtol)

  prog = ProgressThresh(rtol, "Bellman: VFI")
  for i in 1:maxiter
    iterate_bellman!(value_new, value_old, policy, policies_full, endo, exo, converged, aggregate_state, params, hh)
    diff = norm(value_old - value_new)

    adj_fact = max(norm(value_old), norm(value_new))

    ProgressMeter.update!(prog, diff / adj_fact)
    value_old .= value_new

    if diff < rtol * adj_fact
      break
    end
    if i == maxiter
      print("\n"^2)
      @warn "reached $maxiter, diff= $diff"
    end
  end
  number_conv = sum(converged)

  length(converged) == number_conv || @warn "Bellman didn't converge at $(round((1-number_conv / length(converged)) * 100, digits=4))% ($(length(converged) - number_conv) states)"

end

############################################################
## Special cases for coupled value functions
############################################################

function solve_bellman(endo, exo, aggregate_state, params, chh::CoupledHouseholds; maxiter=200, rtol=eps()^0.4)
  containers = initialize_values_policies(endo, exo, aggregate_state, params, chh)
  @unpack W_old, W_new, V, policy, policies_full, policy_hh, converged = containers

  solve_bellman!(W_old, W_new, V, policy, policies_full, policy_hh, endo, exo, converged, aggregate_state, params, chh; maxiter=maxiter, rtol=rtol)

  (val = W_new, policy = policy, policy_hh=policy_hh, policies_full=StructArray.(policies_full), converged=converged)
end

function initialize_values_policies(endo, exo, aggregate_state, params, chh::CoupledHouseholds)
  container_size = (length(endo), length(exo))

  hh_tup = households(chh)
  n = length(chh)

  V         = [zeros(container_size) for i in 1:n]
  W_old     = [zeros(container_size) for i in 1:n]
  W_new     = [zeros(container_size) for i in 1:n]
  converged = [trues(container_size) for i in 1:n]

  policy_hh = [zeros(Int, container_size) for i in 1:n]

  # @unpack proto_pol, proto_pol_full
  proto = map(1:n) do i
    proto_policy(endo, exo, V[1], aggregate_state, params[i], hh_tup[i])
  end

  proto_pol = [p.proto_pol for p in proto]
  proto_pol_full = [p.proto_pol_full for p in proto]

  policy = fill.(proto_pol, Ref(container_size))
  policies_full = fill.(proto_pol_full, Ref(container_size))

  (W_old=W_old, W_new=W_new, V=V, policy=policy, policies_full=policies_full, policy_hh=policy_hh, converged=converged)
end

function solve_bellman!(W_old::Vector, W_new::Vector, V::Vector, policy::Vector, policies_full::Vector, policies_hh::Vector, endo, exo, converged::Vector, aggregate_state, params::Vector, chh::CoupledHouseholds; maxiter, rtol)

  prog = ProgressThresh(rtol, "Solving Bellman equation")

  for i in 1:maxiter
    iterate_bellman!.(V, W_old, policy, policies_full, Ref(endo), Ref(exo), converged, Ref(aggregate_state), params, households(chh))

    update_coupled_values!(W_new, V, policies_hh)

    diff = norm.(W_old .- W_new)

    adj_fact = max.(norm.(W_old), norm.(W_new))

    relative_error = maximum(diff ./ adj_fact)

    ProgressMeter.update!(prog, relative_error)

    for i in 1:length(W_old)
      W_old[i] .= W_new[i]
    end

    if relative_error < rtol
      break
    end

    if i == maxiter
      print("\n"^2)
      @warn "reached $maxiter, diff= $relative_error"
    end
  end

  all(all.(converged)) || @warn "optimization didn't converge at $((1 .- mean.(converged)) .* 100)%"

end

#function update_coupled_values! end # needs to be provided

function update_coupled_values!(W, V, policies_hh)
  W_own, W_rent = W
  V_own, V_rent = V
  policies_own, policies_rent = policies_hh

  W_own  .= max.(V_own, V_rent)
  policies_own .= V_own .> V_rent

  W_rent .= max.(V_own, V_rent)
  policies_rent .= V_own .> V_rent

  nothing
end


text:
using DataDeps

register(DataDep(
  "ecoli",
  "http://archive.ics.uci.edu/ml/datasets/Ecoli",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/ecoli/ecoli.data",
  "008bd8fbb1d8b34040c3c8c4e987cac2a7ebf116e008b140cc6441be5261ba1d",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=9,
        feature_indices=2:8,
        id_indices=1,
        delim=' ',
        ignorerepeated=true,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
print("threshold: execute functions ...\n")

#threshold
let
    x = randn(200)*2
    for th in (BiggestTH(), HardTH(), SoftTH(), SemiSoftTH(), SteinTH())
        threshold(x, th, 2)
    end
    for th in (PosTH(), NegTH())
        threshold(x, th)
    end
end

# denoising
let
    vs = VisuShrink(10)
    n = 2^8
    x0 = testfunction(n, "Doppler")
    x = x0 + 0.05*randn(n)
    y = denoise(x, TI=true)
    y = denoise(x, TI=true, nspin = 8)
    y = denoise(x, TI=false)
    y = denoise(x, nothing)
    y = denoise(randn(32,32), TI=true)
end

# best basis
let
    wt = wavelet(WT.db4)
    x = sin.(4*range(0, stop=2*pi-eps(), length=1024))
    tree = bestbasistree(x, wt)
    xtb = wpt(x, wt, tree)
    @test iwpt(xtb, wt, tree) ≈ x

    x = sin.(4*range(0, stop=2*pi-eps(), length=5*64)) # non-dyadic
    tree = bestbasistree(x, wt)
    xtb = wpt(x, wt, tree)
    @test iwpt(xtb, wt, tree) ≈ x
end

#matching pursuit
let
    N = 128
    M = 64
    y=randn(N)
    Arr = randn(M,N)
    func(a::AbstractVector) = Arr*a
    funct(a::AbstractVector) = Arr'*a
    x = func(y)
    matchingpursuit(x, func, funct, 0.1)
end


text:
module Optimization

using Logging
using KrylovMethods

using DeformableRegistration
include("helpers/objectiveFunctionCreation.jl")

export checkStoppingCriteria, ArmijoLineSearch, optimizeGaussNewton

function optimizeGaussNewton(Jfunc::Function,  # objective Function
                             y::Array{Float64,1}, yInitial::Array{Float64,1}, options::regOptions)

    # Initilization of JRef: reference value of objective Function
    JRef = Jfunc(yInitial)[1]
    JOld = NaN

    output = (Logging.LogLevel == Logging.DEBUG) |  (Logging.LogLevel == Logging.INFO)

    for iter = 1:options.maxIterGaussNewton

        # get derivatives of objective function
        J, dJ, d2J = Jfunc(y,doDerivative=true,doHessian=true)

        # solve system of linear equations d2J*dy=-dJ
          # (dy: change of the variables / search direction)
        # number of parameters < 10 => use backslash operator
        # else (large systems) => use conjugate gradient (cg) method

        cgIterations=0

        if(length(y)<10 && !(typeof(d2J) <: Function))
            dy=d2J\-dJ
        else
            dy,flag,resvec,cgIterations = KrylovMethods.cg(d2J,-dJ,maxIter=options.maxIterCG, tol=1e-5)[1:4]
        end

		    # check descent direction
		    if( (dJ'*dy)[1] > 0)
            dy = -dy
            warn("Changing sign of computed descent direction. This is a suspicious move.")
		    end

        # armijo line search (LS) method
        stepLength,LSiter,LSfailed = ArmijoLineSearch(Jfunc,J,dJ,y,dy)
        if(LSfailed)
            break
        end

        # output
        if(cgIterations==0)
          s = @sprintf("%3d: J %8.4e     LSiter: %2d    J/Jref: %1.2f \n",iter, J[1], LSiter, J[1]/JRef[1])
          info(s)
        else
          s = @sprintf("%3d: J %8.4e     LSiter: %2d     CGiter: %3d     J/Jref: %1.2f \n",iter, J[1], LSiter,cgIterations, J[1]/JRef[1])
          info(s)
        end

        # update parameter y
        y = y + stepLength.*dy

        # stopping criteria
        if (iter>1)
            if checkStoppingCriteria(J, JOld, JRef, dJ, y, stepLength*dy) #printActiveStoppingCirteria=output
                break
            end
        end

        # updating JOld
        JOld = J[1]

    end

    return y

end

function ArmijoLineSearch(Jfunc::Function,         # objective function
                          Jc::Float64,             # Jc = Jfunc(y)
                          dJ::Array{Float64,1},    # gradient of the objective funtion
                          y::Array{Float64,1},     # old variables
                          dy::Array{Float64,1};    # change of the variables / search direction
                          tolLS::Float64 = 1e-4)   # tolerance

    stepLength = 1.0; LSiter = 1; LSfailed = false;

    while( (Jfunc(y+stepLength*dy)[1]) >= (Jc + (stepLength.*tolLS.*dJ'*dy)[1]) )

        if(stepLength < 1e-5)
            LSfailed = true
            break
        end

        stepLength = 0.5 * stepLength
        LSiter = LSiter + 1

    end

    if(LSfailed)
        s = @sprintf("   0. Line search failed after %2d iterations.",LSiter)
        Logging.debug(s)
    end

    return stepLength,LSiter,LSfailed

end

function checkStoppingCriteria(J,JOld,JRef,    # value of the current objective function J(y+dy), the old J(y) and the reference J(yInital)
                               dJ,             # gradient of the objective funtion
                               y,              # old variables
                               dy;             # change of the variables / search direction
                               tolJ = 1e-3,    # tolerance: change of the objective function
                               tolY = 1e-2,    # tolerance: change of the variables
                               tolG = 1e-2,    # tolerance: change of the gradient
                               tolQ = 1e-5     # tolerance: change of quotient
                               )

    # Initialize STOP-Array with false
    STOP = Array{Bool}(5); STOP[:] = false

    # Either STOP[1:3] or STOP[4]
    STOP[1] = abs(JOld-J) <= tolJ  * (1+abs(JOld))
    STOP[2] = norm(dy)    <= tolY  * (1+norm(y))
    STOP[3] = norm(dJ)    <= tolG  * (1+abs(JOld))
    STOP[4] = norm(dJ)    <= 1e6 * eps()
    STOP[5] = abs((JOld-J)/(JRef-J)) <= tolQ

    if(all(STOP[1:3]))
        s = @sprintf("STOPPING CRITERIA:\n")
        Logging.info(s)
        s = @sprintf("   1. |JOld-J| = %5e <= %5e \n",abs(JOld-J), tolJ*(1+abs(JOld)) )
        Logging.info(s)
        s = @sprintf("   2. ||dy|| = %5e <= %5e \n",norm(dy),    tolY*(1+norm(y))  )
        Logging.info(s)
        s = @sprintf("   3. ||dJ|| = %5e <= %5e \n",norm(dJ),    tolG*(1+abs(JOld)) )
        Logging.info(s)
    end

    if(STOP[4])
        s = @sprintf("STOPPING CRITERION:\n")
        Logging.info(s)
        s = @sprintf("   4. ||dJ|| = %5e <= %5e \n",norm(dJ), 1e6*eps())
        Logging.info(s)
    end

    if(STOP[5])
        s = @sprintf("STOPPING CRITERION:\n")
        Logging.info(s)
        s = @sprintf("   5. |(JOld-J)/(JRef-J)|= %5e <= %5e \n",abs((JOld-J)/(JRef-J)), tolQ)
        Logging.info(s)
    end

    return all(STOP[1:3]) | STOP[4] | STOP[5]

end

end


text:
"""
    MESS{T_samplespace<:SampleSpace,TR<:Real} <: PointAlternative

    MESS{T_samplespace<:SampleSpace}(
      samplespace::T_samplespace,
      p0, p1,
      alpha, beta,
      pess;
      minstoppingforfutility::Real   = 0.0,
      minconditionalpower::Real      = 0.0,
      MONOTONECONDITIONALPOWER::Bool = true
    )

This type represents a set of parameters for finding optimal two-stage designs
minimizing the expected sample size on a point in the parameter space subject to
type one and two error rate constraints (Minimal Expected Sample Size).

    > [1] Simon R. Optimal two-stage designs for phase II clinical trials. `Controlled Clinical Trials` 1989; 10, 1-10.

    > [2] Kunzmann K and Kieser M. Optimal adaptive two-stage designs for single-arm trials with binary endpoint. `arxive.org` 2016; arXiv:1605.00249.

# Details

Please note that due to technical reasons the power constraint is actually imlemented in a relaxed way by
adding lambda*(1-beta - power)_+ to the objective function. In practice, lambda should be large enough
to ensure a power of 1-beta, but depending on the setting, a larger lambda might be required.

# Parameters

| Parameter    | Description |
| -----------: | :---------- |
| samplespace  | a sample space object |
| p0           | upper boundary of the null hypothesis |
| p1           | point alternative to power on |
| alpha        | maximal tolerable type one error rate |
| beta         | maximal tolerable type two error rate on p1 |
| pess         | response rate under which to minimize expected sample size |
| minstoppingforfutility | minimal probability for stopping for futility under p0 |
| minconditionalpower | minimal conditional power upon continuation to stage two |
| MONOTONECONDITIONALPOWER | if true, the conditional power must be monotonously increasing, this constraint is only relevant if nmax is set very restrictively |
"""
mutable struct MESS{TR<:Real} <: PointAlternative

  label::String
  samplespace::SampleSpace
  p0::TR
  p1::TR
  mtoer::TR
  beta::TR
  pess::TR
  minconditionalpower::TR
  MONOTONECONDITIONALPOWER::Bool
  monotoneconditionalpowerpivots
  UNIMODAL::Bool
  minstoppingforfutility::TR

  # penalty parameter for relaxing power
  lambda::TR

  function MESS{TR}(
    label::String,
    samplespace::SampleSpace,
    p0::TR, p1::TR,
    alpha::TR, beta::TR,
    pess::TR,
    minconditionalpower::TR,
    MONOTONECONDITIONALPOWER::Bool,
    monotoneconditionalpowerpivots,
    UNIMODAL::Bool,
    minstoppingforfutility::TR,
    lambda::TR
  ) where {TR<:Real}

    @checkprob p0 p1 alpha beta pess minconditionalpower minstoppingforfutility
    p0 >= p1 ?
      error("p0 must be smaller than p1") : nothing
    new(
      label, samplespace, p0, p1, alpha, beta, pess, minconditionalpower,
      MONOTONECONDITIONALPOWER, monotoneconditionalpowerpivots, UNIMODAL, minstoppingforfutility, lambda
    )

  end # inner constructor

end # MESS

function MESS(
  samplespace::SampleSpace,
  p0::Real, p1::Real;
  alpha::Real = 0.05, beta::Real = 0.2,
  pess::Real = p1,
  minstoppingforfutility::Real = 0.0,
  minconditionalpower::Real = 0.0,
  MONOTONECONDITIONALPOWER::Bool = true,
  monotoneconditionalpowerpivots = [p1],
  UNIMODAL::Bool = true,
  label::String = "",
  lambda = 10.0^9
)

  p0, p1, alpha, beta, pess, minstoppingforfutility, minconditionalpower, lambda =
    promote(p0, p1, alpha, beta, pess, minstoppingforfutility, minconditionalpower, lambda)

  MESS{typeof(p0)}(
    label, samplespace, p0, p1, alpha, beta, pess, minconditionalpower,
    MONOTONECONDITIONALPOWER, monotoneconditionalpowerpivots, UNIMODAL, minstoppingforfutility, lambda
  )

end # MESS

maxsamplesize(params::MESS) = maxsamplesize(params.samplespace)

isgroupsequential(params::MESS) = isgroupsequential(params.ss)

hasmonotoneconditionalpower(params::MESS) = params.MONOTONECONDITIONALPOWER

isunimodal(params::MESS) = params.UNIMODAL

minconditionalpower(params::MESS) = params.minconditionalpower

Base.show(io::IO, params::MESS) = print("MESS")

function score(design::Design, params::MESS, p::T) where {T<:Real}

  @checkprob p
  n = SampleSize(design, p)
  return mean(n)

end # score

function score(design::Design, params::MESS)

  n = SampleSize(design, params.pess)
  return mean(n)

end # score

function completemodel(ipm::IPModel, params::MESS, n1::T) where {T<:Integer}

  ss = samplespace(params)
  !possible(n1, ss) ? warn("completemodel(EB): n1 and sample space incompatible") : nothing
  # extract ip model
  m             = ipm.m
  y             = ipm.y
  nvals         = ipm.nvals
  cvals         = ipm.cvals
  cvalsfinite   = ipm.cvalsfinite
  cvalsinfinite = [-Inf; Inf]
  # extract other parameters
  nmax          = maxsamplesize(ss, n1)
  p0            = null(params)
  p1            = params.p1
  # add type one error rate constraint
  JuMP.@constraint(m,
      sum(dbinom(x1, n1, p0)*_cpr(x1, n1, n, c, p0)*y[x1, n, c] for
          x1 in 0:n1,
          n  in nvals,
          c  in cvals
      ) <= mtoer(params)
  )
  JuMP.@expression(m, power,
    sum(dbinom(x1, n1, p1)*_cpr(x1, n1, n, c, p1)*y[x1, n, c] for
      x1 in 0:n1,
      n  in nvals,
      c  in cvals
    )
  )
  JuMP.@constraint(m, power >= 1 - params.beta)
  # add conditional type two error rate constraint (power)
  for x1 in 0:n1
      # ensure monotonicity if required
      if (x1 >= 1) & hasmonotoneconditionalpower(params)
        for p in params.monotoneconditionalpowerpivots
          JuMP.@constraint(m,
              sum(_cpr(x1, n1, n, c, p)*y[x1, n, c] - _cpr(x1 - 1, n1, n, c, p)*y[x1 - 1, n, c] for
                  n in nvals, c in cvals
              ) >= 0
          )
        end
      end
      JuMP.@constraint(m,
          sum(_cpr(x1, n1, n, c, p1)*y[x1, n, c] for
              n  in nvals,
              c  in cvalsfinite
          ) + sum(y[x1, n, c] for
              n  in nvals,
              c  in cvalsinfinite
          ) >= minconditionalpower(params)
      )
  end
  # add constraint for minimal stopping-for-futility probability
  JuMP.@constraint(m,
      sum(dbinom(x1, n1, p0)*y[x1, n1, Inf] for
          x1 in 0:n1
      ) >= params.minstoppingforfutility
  )
  # add optimality criterion
  JuMP.@objective(m, Min,
      sum(dbinom(x1, n1, params.pess)*n*y[x1, n, c] for
          x1 in 0:n1,
          n  in nvals,
          c  in cvals
      )
  )
  return true

end # completemodel

function _isfeasible(design::Design, params::MESS)

    all(power.(design, linspace(0, null(params))) .<= mtoer(params) + .001) ? nothing : throw(InexactError())
    (power(design, params.p1) >= 1 - params.beta - .001) ? nothing : throw(InexactError())
    return true

end # _isfeasible

max_stars_repo_name: kkmann/binaryTwoStageDesigns

max_issues_repo_name: kkmann/binaryTwoStageDesigns

max_forks_repo_name: kkmann/binaryTwoStageDesigns


text:
module Horde

include("GVF/GVF.jl")
export GVF

end # module


text:
include("sphmagic.jl")

function computeSimulationStep(step, x, ρ, u, v, p, L, H)
    # First Leap-Frog integration
    if step > 1
        # Keep track of previous energy and velocity
        previous_u = copy(u)
        previous_v = copy(v)

        # Compute half-step integration
        for i in FIXED_SIDES_PARTICLES_COUNT:PARTICLES_COUNT-FIXED_SIDES_PARTICLES_COUNT
            u[i] += TIME_DELTA/2 * H[i]
            v[i] += TIME_DELTA/2 * L[i]
        end
    end

    (x, ρ, u, v, p, L, H) = doSPHmagic(x, ρ, u, v, p, L, H)

    # Second Leap-Frog integration
    for i in FIXED_SIDES_PARTICLES_COUNT:PARTICLES_COUNT-FIXED_SIDES_PARTICLES_COUNT
        if step > 1
            # Compute half-step integration with previous value
            u[i] = previous_u[i] + TIME_DELTA * H[i]
            v[i] = previous_v[i] + TIME_DELTA * L[i]
        else
            # Compute first half-step integration
            u[i] += TIME_DELTA * H[i]
            v[i] += TIME_DELTA * L[i]
        end
        # Compute regular integration
        x[i] += TIME_DELTA * v[i]
    end

    (x, ρ, u, v, p, L, H)
end


text:
module Penrose

export ζ, ϕ, Qζ, hkite, hdart,  penrose, PenroseElem, forced_penrose, conj, galois_gen

using Luxor
using StructEquality

using ...CoreDefs
using ...NumFields

@NumFields.simple_number_field_concrete Qζ [-1, 1, -1, 1] ζ
promote_rule(::Type{Qζ}, ::Type{Rational{Int64}}) = Qζ

const ϕ = ζ + ζ^9
const L = Qζ

function conj(x :: Qζ)
    return embed_field(a -> a, ζ^9, x)
end

function galois_gen(x :: Qζ)
    return embed_field(a -> a, ζ^3, x)
end

function Base.://(x :: Qζ, y :: Qζ)
    y1 = galois_gen(y)
    y2 = galois_gen(y1)
    y3 = galois_gen(y2)
    alg_conj = y1*y2*y3
    norm = y*alg_conj
    return x*alg_conj*norm.denom//norm.coeffs[1]
end

@def_structequal struct PenroseElem <: DGroupElem
    rot :: Int
    refl :: Bool
    z :: Qζ
end

function Base.:*(g :: PenroseElem, h :: PenroseElem)
    if !g.refl
        return PenroseElem(
            mod(g.rot + h.rot, 10),
            h.refl,
            g.z + ζ^(10-g.rot) * h.z
        )
    else
        return PenroseElem(
            mod(g.rot - h.rot, 10),
            !h.refl,
            g.z + ζ^(10-g.rot) * conj(h.z)
        )
    end
end
function Base.:*(g :: PenroseElem, x :: Qζ)
    if !g.refl
        return g.z + ζ^(10-g.rot) * x
    else
        return g.z + ζ^(10-g.rot) * conj(x)
    end
end
function CoreDefs.dilate(λ :: Qζ, g :: PenroseElem)
    return PenroseElem(g.rot, g.refl, λ*g.z)
end
function CoreDefs.id(::Type{PenroseElem})
    return PenroseElem(0,0,L(0))
end
function Base.inv(p::PenroseElem)
    return PenroseElem(0,p.refl,L(0))*PenroseElem(10-p.rot,0,L(0))*PenroseElem(0,0,-p.z)
end

@enum PenrosePTile begin
    Hkite=1
    Hdart=2
end
function hkite(r,s,z)
    return PenroseElem(r,s,z) => Hkite
end
function hkite()
    return hkite(0,0,L(0))
end
function hdart(r,s,z)
    return PenroseElem(r,s,z) => Hdart
end
function hdart()
    return hdart(0,0,L(0))
end

function embed_nf(x)
    ζ_float = complex(cos(2*pi/10), sin(2*pi/10))
    return embed_field(Complex{Float64}, ζ_float, x)
end
function embed_nf_p(x)
    z = embed_nf(x)
    return Point(real(z), imag(z))
end
function CoreDefs.embed_aff(g :: PenroseElem)
    m = [cos(2*pi/10) sin(2*pi/10); -sin(2*pi/10) cos(2*pi/10)]^g.rot

    z = embed_nf(g.z)
    return [m[1,1], m[2,1], (-1)^g.refl*m[1,2], (-1)^g.refl*m[2,2], real(z), imag(z)]
end

function CoreDefs.draw(ptile::PenrosePTile, action)
    if ptile == Hkite
        Luxor.poly([
            embed_nf_p(L(1)),
            embed_nf_p(L(ζ^4)),
            embed_nf_p(L(ζ^6)),
        ], close = true, action)
    else
        Luxor.poly([
            embed_nf_p(L(1-ζ^0//ϕ)),
            embed_nf_p(L(ζ^2-ζ^0//ϕ)),
            embed_nf_p(L(ζ^8-ζ^0//ϕ))
        ], close = true, action)
    end
end

function in_interval(x :: Qζ)
    x_float = embed_nf(x)
    return 0 ≤ real(x_float) && real(x_float) ≤ 1 && imag(x_float) ≈ 0
end
function in_border(x, ptile :: PenrosePTile)
    if ptile == Hkite
        return any([
            in_interval( (x - 1)//(ζ^4 - 1)),
            in_interval( (x - ζ^4)//(ζ^6 - ζ^4)),
            in_interval( (x - ζ^6)//(1 - ζ^6)),
        ])
    else
        return any([ # TODO fix!
            in_interval( (x - 1+ζ^0//ϕ)//(ζ^2 - 1)),
            in_interval( (x - ζ^2+ζ^0//ϕ)//(ζ^8 - ζ^2)),
            in_interval( (x - ζ^8+ζ^0//ϕ)//(1 - ζ^8)),
        ])
    end
end

function penrose()
    # This is a function only so that Revise works
    pen_subst = Dict([
        (Hkite, [
            hkite(7, 0, ϕ*ζ^4-ζ^3),
            hkite(6, 1, ζ^4//ϕ),
            hdart(7, 0, ζ^4//ϕ+1+(1-ζ^0//ϕ)*ζ^8),
            ]),
        (Hdart, [
            hdart(4,0, ϕ*ζ^8-ζ^0+ζ+ζ^2+(ζ^8-ζ^0//ϕ)*ζ^7),
            hkite(3,1,ϕ*ζ^8-ζ^0+ζ^2)
        ])
    ])
    return SubSystem(pen_subst, ϕ) :: CoreDefs.SubSystem{PenroseElem, Qζ, PenrosePTile}
end

function color(x)
    l = x[2] == Hdart ? 0 : 1
    return 1 + l + 2*x[1].refl
end

function in_bounds(tile, n, window)
    phi = (1 + sqrt(5))/2
    center = embed_nf(tile[1].z)
    return abs(real(center))*phi^(n+1) < window.w + 3/2*phi^(n+1) && abs(imag(center))*phi^(n+1) < window.h + 3/2*phi^(n+1)
end

function force(tiling)
    fkite = [
        hkite(0,0,L(0))
        hkite(9, 1, ζ^8//ϕ)
    ]

    fdart = unique([
        hdart(0,0,L(0));
        hdart(7,1,L(0));
        PenroseElem(9, 1, ζ^2)*fkite;
        PenroseElem(8, 0, ζ)*fkite
    ])

    result = convert(Vector{Pair{PenroseElem, PenrosePTile}}, [])

    for tile in tiling
        if tile[2] == Hdart
            append!(result, tile[1]*fdart)
        else
            append!(result, tile[1]*fkite)
        end
    end
    return result
end

function forced_penrose()
    forced_subst = Dict([
        (label, force(penrose().sub[label]) ) for label in [Hkite, Hdart]
    ])

    return SubSystem(forced_subst, ϕ) :: CoreDefs.SubSystem{PenroseElem, Qζ, PenrosePTile}
end

@NumFields.simple_number_field_concrete Qψ [1,-1] ψ

function embed_psi(x)
    psi_float = (sqrt(5) - 1)/2
    return embed_field(rational_to_float, psi_float, x)
end
# ψ = 1/ϕ, except that ψ is an element of the field Q(ψ) and ϕ is an element of the field Q(ζ_10)

function frequency(patch :: AbstractVector, depth)
    return frequency(Dict(patch), depth)
end
function frequency(patch :: AbstractDict, depth)
    freq = Qψ(0)
    harmonic = Dict([
        (Hkite, (ψ^(2*depth - 1))),
        (Hdart, (ψ^(2*depth)))
    ])

    center_ptile = patch[id(PenroseElem)]

    for label in [Hkite, Hdart]
        domain = substitute(penrose(), [(PenroseElem(0,0,L(0)) => label)], depth-1)
        forced_domain = Dict(substitute(forced_penrose(), ([PenroseElem(0,0,L(0)) => label]), depth-1))
        for tile in domain
            if tile[2] == center_ptile
                translated_patch = tile[1] * patch
                if translated_patch ⊆ forced_domain
                    freq += harmonic[label]
                end
            end
        end
    end
    return freq
end

end


text:
# ---
# title: 537. Complex Number Multiplication
# id: problem537
# author: Indigo
# date: 2021-06-23
# difficulty: Medium
# categories: Math, String
# link: <https://leetcode.com/problems/complex-number-multiplication/description/>
# hidden: true
# ---
#
# Given two strings representing two [complex
# numbers](https://en.wikipedia.org/wiki/Complex_number).
#
# You need to return a string representing their multiplication. Note i2 = -1
# according to the definition.
#
# **Example 1:**
#
#
#
#     Input: "1+1i", "1+1i"
#     Output: "0+2i"
#     Explanation: (1 + i) * (1 + i) = 1 + i2 + 2 * i = 2i, and you need convert it to the form of 0+2i.
#
#
# **Example 2:**
#
#
#
#     Input: "1+-1i", "1+-1i"
#     Output: "0+-2i"
#     Explanation: (1 - i) * (1 - i) = 1 + i2 - 2 * i = -2i, and you need convert it to the form of 0+-2i.
#
#
# **Note:**
#
#   1. The input strings will not have extra blank.
#   2. The input strings will be given in the form of **a+bi** , where the integer **a** and **b** will both belong to the range of [-100, 100]. And **the output should be also in this form**.
#
#
## @lc code=start
using LeetCode

function complex_number_multiply(num1::String, num2::String)
    Ct = Complex{Int}
    mult = parse(Ct, num1) * parse(Ct, num2)
    return "$(mult.re)+$(mult.im)i"
end
## @lc code=end


text:
using Bokeh
autoopen(true)

plotfile("glyphs.html")

x = linspace(0, 2pi, 10)
y = ones(x) .+ 0.1 .* sin.(x)
glyphs = "rs|bo|g*|kd|m+|cx|y^|rv|.g|oxb|o+b|sxg|s+g"
ys = Vector[i*y for i in -1:-1:-length(split(glyphs, '|'))]
plot(ys, glyphs)


text:
#=

=#

"""
    Struct to define the properties of a strongly connected component(SCC) object.
    No arguments is needed to construct this object.

    Properties:
        number_nodes:
            Integer value for the size of the SCC.
        have_input:
            Boolean value to signal whether the SCC has input from other SCC.
        nodes:
            Array of integers regarding the indexes of the nodes belonging
            to the fiber.
        type:
            Integer value ranging from 0 to 2. Each value correspond to a different
            type of SCC. This classify is related to the inputs that the SCC receives
            and it is used for a proper initialization of the fibration partitioning
            algorithm.
"""
mutable struct StrongComponent
    number_nodes::Int
    have_input::Bool
    nodes::Array{Int}
    type::Int
    function StrongComponent()
        number_nodes = 0
        have_input = false
        nodes = Int[]
        type = -1
        new(number_nodes, have_input, nodes, type)
    end
end

"""
    Add 'node' to 'strong' object. If 'node' is an array
    of nodes, then all nodes are inserted.
"""
function insert_nodes(node::Int, strong::StrongComponent)
    strong.number_nodes += 1
    append!(strong.nodes, [node])
end

function insert_nodes(node::Array{Int}, strong::StrongComponent)
    strong.number_nodes += length(node)
    append!(strong.nodes, node)
end

function get_nodes(strong::StrongComponent)
    return strong.nodes
end

function get_input_bool(strong::StrongComponent)
    return strong.have_input
end

"""
    Check if the given SCC receives or not input from another
    components in the 'graph'.

    the field 'have_input' of 'strong' is modified to 'true'
    if the component receives external information. Otherwise,
    'have_input' maintains its default ('false').
"""
function check_input(strong::StrongComponent, graph::Graph)
    from_out = false
    for u in strong.nodes
        input_nodes = get_in_neighbors(u, graph)
        for w in input_nodes
            if w in strong.nodes
                from_out = false
            else
                from_out = true
            end
            if from_out
                strong.have_input = true
                break
            end
        end
        if strong.have_input
            break
        end
    end
    return
end

"""
    This function should be called after 'check_input'
"""
function classify_strong(strong::StrongComponent, graph::Graph)
    if strong.have_input
        strong.type = 0
    else
        """
            If it doesn't receive any external input, then we
            must check if it is an isolated self-loop node.
        """
        if length(strong.nodes)==1
            in_neighbors = get_in_neighbors(strong.nodes[1], graph)
            if length(in_neighbors)==0
                strong.type = 1
            else
                strong.type = 2 # Isolated self-loop node.
            end
        else
            strong.type = 1 # SCC does not have external input.
        end
    end
    return
end


text:
# Example for basic 2D modeling:
# The receiver positions and the source wavelets are the same for each of the four experiments.
# Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: January 2017
#

using JUDI, PyPlot, JLD, SegyIO, LinearAlgebra

# Load Sigsbee model
if !isfile("sigsbee2A_model.jld")
    run(`wget ftp://slim.gatech.edu/data/SoftwareRelease/Imaging.jl/CompressiveLSRTM/sigsbee2A_model.jld`)
end

M = load("sigsbee2A_model.jld")
m = M["m"][1200:2200, 150:800]
m0 = M["m0"][1200:2200, 150:800]
dm = M["dm"][1200:2200, 150:800]
n = size(m)

# Setup info and model structure
model = Model(n, M["d"], M["o"], m)
model0 = Model(n, M["d"], M["o"], m0)
dm = vec(dm)

## Set up receiver geometry
nsrc = 1    # number of sources
nxrec = 1200
xrec = range(700f0, stop=6700f0, length=nxrec)
yrec = 0f0
zrec = range(50f0, stop=50f0, length=nxrec)

# receiver sampling and recording time
timeR = 8000f0   # receiver recording time [ms]
dtR = 4f0    # receiver sampling interval

# Set up receiver structure
recGeometry = Geometry(xrec, yrec, zrec; dt=dtR, t=timeR, nsrc=nsrc)

## Set up source geometry (cell array with source locations for each shot)
xsrc = convertToCell([6800f0])
ysrc = convertToCell([0f0])
zsrc = convertToCell([20f0])

# source sampling and number of time steps
timeS = 8000f0
dtS = 2f0

# Set up source structure
srcGeometry = Geometry(xsrc, ysrc, zsrc; dt=dtS, t=timeS)

# setup wavelet
f0 = 0.015
wavelet = ricker_wavelet(timeS,dtS,f0)
q = judiVector(srcGeometry,wavelet)

# Set up info structure for linear operators
ntComp = get_computational_nt(srcGeometry,recGeometry,model)
info = Info(prod(model.n),nsrc,ntComp)

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

opt = Options(isic=false,
              save_data_to_disk=false,
              optimal_checkpointing=false
              )

# Setup operators
Pr = judiProjection(info, recGeometry)
F = judiModeling(info, model; options=opt)
F0 = judiModeling(info, model0; options=opt)
Ps = judiProjection(info, srcGeometry)
J = judiJacobian(Pr*F0*Ps', q)

# Set frequencies
q_dist = generate_distribution(q)
J.options.frequencies = Array{Any}(undef, nsrc)
J.options.frequencies[1] = select_frequencies(q_dist; fmin=0.003, fmax=0.04, nf=20)

# Velocity contrast
d_lin1 = J*dm
rtm1 = J'*d_lin1

# Impedance
J.options.isic = true
d_lin2 = J*dm
rtm2 = J'*d_lin2

############################################## Plots ###############################################

# Normalize images
dm1 = rtm1 / norm(rtm1, Inf)
dm2 = rtm2 / norm(rtm2, Inf)

# Plot migration velocity
fig=figure(figsize=(6.66, 5))
s1 = subplot(2, 2, 1)
im1 = imshow(adjoint(sqrt.(1f0 ./ m0)), vmin=1.5, vmax=4.511,cmap="gray", extent=(9.144, 16.764, 6.096, 1.143))
s1[:tick_params]("both", labelsize=8)
xlabel("Lateral position [km]", size=8);
ylabel("Depth [km]", size=8)
t1 = title("a)", fontweight="bold", size=10, loc="left")
t1["set_position"]([-.14, 1])

# Plot true image
s2 = subplot(2, 2, 2)
im2 = imshow(adjoint(reshape(dm, n)), vmin=-6e-2, vmax=6e-2,cmap="gray", extent=(9.144, 16.764, 6.096, 1.143))
s2[:tick_params]("both", labelsize=8)
xlabel("Lateral position [km]", size=8);
ylabel("Depth [km]", size=8)
t2 = title("b)", fontweight="bold", size=10, loc="left")
t2["set_position"]([-.14, 1])

# RTM image w/ zero-lag cross-correlation imaging condition
s3 = subplot(2, 2, 3)
im3 = imshow(adjoint(reshape(dm1,n)), vmin=-1e-1, vmax=1e-1,cmap="gray", extent=(9.144, 16.764, 6.096, 1.143))
s3[:tick_params]("both", labelsize=8)
xlabel("Lateral position [km]", size=8);
ylabel("Depth [km]", size=8)
t3 = title("c)", fontweight="bold", size=10, loc="left")
t3["set_position"]([-.14, 1])

# RTM image w/ inverse-scattering/impedance imaging condition
s4 = subplot(2, 2, 4)
im4 = imshow(adjoint(reshape(dm2, n)), vmin=-4e-2, vmax=4e-2,cmap="gray", extent=(9.144, 16.764, 6.096, 1.143))
s4[:tick_params]("both", labelsize=8)
xlabel("Lateral position [km]", size=8);
ylabel("Depth [km]", size=8)
t4 = title("d)", fontweight="bold", size=10, loc="left")
t4["set_position"]([-.14, 1])

tight_layout()

savefig("figure1", dpi=300, format="png")
run(`mogrify -trim figure1`)

max_stars_repo_name: nogueirapeterson/JUDI

max_issues_repo_name: nogueirapeterson/JUDI

max_forks_repo_name: nogueirapeterson/JUDI


text:
using DataFrames
using DataFramesMeta
using CSV
using Iterators

solver_order = ["pomcpow" => "POMCPOW",
                "pft" => "PFT-DPW",
                "qmdp" => "QMDP",
                "pomcpdpw" => "POMCP-DPW",
                "despot" => "DESPOT",
                "d_pomcp" => "POMCP\\textsuperscript{D}",
                "d_despot" => "DESPOT\\textsuperscript{D}"]

cardinality = Dict("lasertag" => "(D, D, D)",
                   "lightdark" => "(D, D, C)",
                   "subhunt" => "(D, D, C)",
                   "vdpbarrier" => "(C, C, C)",
                   "vdptag" => "(C, C, C)",
                   "multilane" => "(C, D, C)"
                  )

problem_order = [["lasertag", "lightdark", "subhunt"], ["vdpbarrier", "multilane"]]
# problem_order = ["multilane"]

filenames = Dict("simpleld"=>"/home/zach/.julia/v0.6/ContinuousPOMDPTreeSearchExperiments/icaps_2018/data/simpleld_Monday_5_Mar_21_39.csv",
                 "subhunt"=>"/home/zach/.julia/v0.6/ContinuousPOMDPTreeSearchExperiments/icaps_2018/data/subhunt_Monday_5_Mar_22_08.csv",
                 "lasertag"=>"/home/zach/.julia/v0.6/ContinuousPOMDPTreeSearchExperiments/icaps_2018/data/lasertag_Monday_5_Mar_20_17.csv",
                 "vdpbarrier"=>"/home/zach/.julia/v0.6/ContinuousPOMDPTreeSearchExperiments/icaps_2018/data/vdpbarrier_Monday_5_Mar_23_01.csv",
                 "multilane"=>"/home/zach/.julia/v0.6/ContinuousPOMDPTreeSearchExperiments/data/multilane_Monday_30_Apr_23_41.csv"
                )
filenames["lightdark"] = filenames["simpleld"]

# filenames = Dict("lasertag" => "$(Pkg.dir("ContinuousPOMDPTreeSearchExperiments"))/icaps_2018/data/lasertag_Monday_26_Feb_18_46.csv",
#                  "lightdark" => "$(Pkg.dir("ContinuousPOMDPTreeSearchExperiments"))/icaps_2018/data/simpleld_Monday_26_Feb_20_13.csv",
#                  "subhunt" => "$(Pkg.dir("ContinuousPOMDPTreeSearchExperiments"))/icaps_2018/data/subhunt_Monday_26_Feb_20_44.csv",
#                  "vdpbarrier" => "$(Pkg.dir("ContinuousPOMDPTreeSearchExperiments"))/icaps_2018/data/bdpbarrier_Monday_26_Feb_21_42.csv")

data = Dict(
    "lightdark" => Dict(
        "limits" => (-20.0, 80.0),
        "name" => "Light Dark"
    ),

    "subhunt" => Dict(
        "limits" => (0.0, 80.0),
        "name" => "Sub Hunt"
    ),

    "lasertag" => Dict(
        "limits" => (-20.0, -8.0),
        "name" => "Laser Tag"
    ),

    "vdptag" => Dict(
        "limits" => (-20.0, 40.0),
        "name" => "VDP Tag"
    ),

    "vdpbarrier" => Dict(
        "limits" => (0.0, 31.0),
        "name" => "VDP Tag"
    ),

    "multilane" => Dict(
        "limits" => (10.0, 40.0),
        "name" => "Multilane",
    )

)

for p in chain(problem_order...)
    df = CSV.read(filenames[p])
    d = data[p]
    for s in unique(df[:solver])
        rs = @where(df, :solver.==s)[:reward]
        m = mean(rs)
        sem = std(rs)/sqrt(length(rs))
        if s == "ar_despot"
            s = "despot"
        elseif s == "despot_01"
            s = "despot"
        elseif p == "lasertag" && s == "pomcp"
            s = "d_pomcp"
        elseif p == "multilane" && s == "pftdpw"
            s = "pft"
        elseif p == "multilane" && s == "qmdp"
            s = "pomcpdpw"
        end
        if p == "multilane"
            m *= 100.0
            sem *= 100.0
        end
        d[s] = (m, sem)
    end
end

hbufs = []
tbufs = []

for po_chunk in problem_order
    hbuf = IOBuffer()

    for k in po_chunk
        d = data[k]
        print(hbuf, "& & $(d["name"]) \\makebox[0pt][l]{$(cardinality[k])} & ")
    end
    print(hbuf, "\\\\")

    tbuf = IOBuffer()
    for (k, n) in solver_order
        print(tbuf, n*" ")
        for p in po_chunk
            d = data[p]
            if haskey(d, k)
                m, sem = d[k]
                lo, hi = d["limits"]
                frac = (m-lo)/(hi-lo)
                @printf(tbuf, "& & \\result{%.1f}{%.1f}{%d}{%.2f} ",
                        m, sem, round(Int, 100*frac), frac
                       )
            else
                print(tbuf,  "& & \\noresult{} ")
            end
        end
        print(tbuf, "\\\\\n")
    end

    push!(hbufs, hbuf)
    push!(tbufs, tbuf)
end

tablebuf = IOBuffer()

columns = "l"*"Xrl"^maximum(length(chunk) for chunk in problem_order)

println(tablebuf, "\\begin{tabularx}{\\linewidth}{$columns}")
println(tablebuf, "\\toprule")

for i in 1:length(hbufs)
    println(tablebuf, String(hbufs[i]))
    println(tablebuf, "\\midrule")
    println(tablebuf, String(tbufs[i]))

    if i < length(hbufs)
        println(tablebuf, "\\midrule")
    else
        println(tablebuf, "\\bottomrule")
        println(tablebuf, "\\end{tabularx}")
    end
end

open("/home/zach/Devel/thesis/results_table.tex", "w") do f
    println(f, String(tablebuf))
end


text:
module FastGaussQuadrature

using LinearAlgebra
using SpecialFunctions
using StaticArrays

export gausslegendre
export gausschebyshev
export gausslaguerre
export gausshermite
export gaussjacobi
export gausslobatto
export gaussradau
export approx_besselroots
export besselroots

import SpecialFunctions: besselj, airyai, airyaiprime

include("constants.jl")
include("gausslegendre.jl")
include("gausschebyshev.jl")
include("gausslaguerre.jl")
include("gausshermite.jl")
include("gaussjacobi.jl")
include("gausslobatto.jl")
include("gaussradau.jl")
include("besselroots.jl")

end


text:
module RobotZoo

using StaticArrays
using Parameters
using RobotDynamics
using Rotations
using LinearAlgebra

import RobotDynamics: dynamics, forces, moments, wrenches, inertia, inertia_inv, orientation
import RobotDynamics: state_dim, control_dim

include("acrobot.jl")
include("car.jl")
include("cartpole.jl")
include("double_integrator.jl")
include("pendulum.jl")
include("quadrotor.jl")
include("planar_quad.jl")
include("yak_plane.jl")
include("satellite.jl")
include("freebody.jl")
include("rocket.jl")

include("LinearModels.jl")

end # module


text:
""" multiple strings replacement, for multiple on Julia version """
function replace_multi(text, pairs...) ::String
    VERSION >= v"1.7" && return replace(text, pairs...)
    VERSION < v"1.7" && begin
        for pair in pairs
            text = replace(text, pair)
        end
    end
    return text
end

# ---------------------------------------------------------------------------- #
#                                     REGEX                                    #
# ---------------------------------------------------------------------------- #
# ---------------------------------- markup ---------------------------------- #

"""
This regex uses lookahead and lookbehind to exclude [[
at the beginning of a tag, with this:
    (?<!\\[)\\[(?!\\[)
"""
const OPEN_TAG_REGEX = r"(?<!\[)\[(?!\[)[a-zA-Z _0-9. ,()#]*\]"
const CLOSE_TAG_REGEX = r"(?<!\[)\[(?!\[)\/[a-zA-Z _0-9. ,()#]+[^/\[]\]"
const GENERIC_CLOSER_REGEX = r"(?<!\[)\[(?!\[)\/\]"

"""
    remove_markup(input_text::AbstractString)::AbstractString

Remove all markup tags from a string of text.
"""
remove_markup(input_text)::String = replace_multi(input_text,
                                        OPEN_TAG_REGEX => "",
                                        GENERIC_CLOSER_REGEX => "",
                                        CLOSE_TAG_REGEX => ""
                                            )

"""
    has_markup(text::String)

Returns `true` if `text` includes a `MarkupTag`
"""
has_markup(text::String)::Bool = occursin(OPEN_TAG_REGEX, text)

# ----------------------------------- ansi ----------------------------------- #
const ANSI_REGEXE = r"\e\[[0-9;]*m"

"""
    remove_ansi(input_text::AbstractString)::AbstractString

Remove all ANSI tags from a string of text
"""
remove_ansi(input_text)::String = replace(input_text, ANSI_REGEXE => "", )

"""
    has_markup(text::String)

Returns `true` if `text` includes a `MarkupTag`
"""
has_ansi(text)::Bool = occursin(ANSI_REGEXE, text)

"""
    get_last_ANSI_code(text)::String

Get the last ANSI code in a sting, returns "" if no ANSI code found.
"""
function get_last_ANSI_code(text)::String
    has_ansi(text) || return ""

    # get the last matching regex
    rmatch = collect((eachmatch(ANSI_REGEXE, text)))[end]
    return rmatch.match
end

"""
    get_ANSI_codes(text)::String

Returns a string with all ANSI codes in the input.
"""
function get_ANSI_codes(text)::String
    has_ansi(text) || return ""
    matches = collect((eachmatch(ANSI_REGEXE, text)))
    return *(map(m->m.match, matches)...)
end

"""
    replace_ansi(input_text)

Replace ANSI tags with ¦.

The number of '¦' matches the length of the ANSI tags.
Used when we want to hide ANSI tags but keep the string length intact.
"""
function replace_ansi(input_text)
    while occursin(rx, input_text)
        mtch = match(ANSI_REGEXE, input_text)
        input_text = replace_text(input_text, mtch.offset-1, mtch.offset+length(mtch.match)-1, '¦')
    end
    return input_text
end

# --------------------------- clean text / text len -------------------------- #
"""
    cleantext(str::AbstractString)

Remove all style information from a string.
"""
cleantext(str)::String = (remove_ansi ∘ remove_markup)(str)

"""
    textlen(x::AbstractString)

Get length of text after all style information is removed.
"""
textlen(x::String)::Int = (textwidth ∘ remove_markup ∘ remove_ansi)(x)
textlen(x::SubString)::Int = (textwidth ∘ remove_markup ∘ remove_ansi)(x)

# --------------------------------- brackets --------------------------------- #
const brackets_regexes = [
    r"(?<!\[)\[(?!\[)",
    r"(?<!\])\](?!\])",
]

"""
    remove_ansi(str)::String

Replace each squared bracket with a double copy of itself
"""
escape_brackets(text)::String = replace_multi(text,
        brackets_regexes[1]=>"[[",
        brackets_regexes[2]=>"]]",
)

const remove_brackets_regexes = [
    r"\[\[",
    r"\]\]",
]

"""
    unescape_brackets(text)::String

Replece every double squared parenthesis with a single copy of itself
"""
unescape_brackets(text)::String = replace_multi(text,
    remove_brackets_regexes[1]=>"[",
    remove_brackets_regexes[2]=>"]",
)

# ---------------------------------------------------------------------------- #
#                                      I/O                                     #
# ---------------------------------------------------------------------------- #
"""
    read_file_lines(path::String, start::Int, stop::Int)

Read a file and select only lines in range `start` -> `stop`.

Returns a vector of tuples with the line number and line content.
"""
function read_file_lines(path::AbstractString, start::Int, stop::Int)
    start = start < 1 ? 1 : start
    stop = stop >= countlines(path) ? countlines(path) : stop
    lines = readlines(path; keep = true)
    return collect(enumerate(lines))[start:stop]
end

# ---------------------------------------------------------------------------- #
#                                     MISC                                     #
# ---------------------------------------------------------------------------- #
"""
    tview(text, start::Int, stop::Int)

Get a view object with appropriate indices
"""
tview(text, start::Int, stop::Int) = view(text,thisind(text, start):thisind(text, stop))
tview(text, start::Int, stop::Int, simple::Symbol) = view(text, start:stop)

"""
    replace_text(text::AbstractString, start::Int, stop::Int, replace::AbstractString)

Replace a section of a `text` between `start` and `stop` with `replace`.
"""
function replace_text(text, start::Int, stop::Int, replace::String)::String
    if start == 0
        return  replace * text[stop+1:end]
    end

    start = isvalid(text, start) ? start : max(prevind(text, start), 1)
    if start == 1
        return text[1] * replace * text[stop+1:end]
    elseif stop == length(text)
        return text[1:start] * replace
    else
        return text[1:start] * replace * text[stop+1:end]
    end
end

"""
    replace_text(text::AbstractString, start::Int, stop::Int, char::Char='_')

Replace a section of a `text`  between `start` and `stop` with another string composed of repeats of a given character `char`.
"""
function replace_text(text, start::Int, stop::Int, char::Char = '_')::String
    replacement = char^(stop - start)
    return replace_text(text, start, stop, replacement)
end

"""
    nospaces(text::AbstractString)

Remove all spaces from a string.
"""
nospaces(text::AbstractString) = replace(text, " " => "")

"""
    remove_brackets(text::AbstractString)

Remove all () brackets from a string.
"""
remove_brackets(text)::String = replace_multi(text, "(" => "", ")" => "")

"""
    unspace_commas(text::AbstractString)

Remove spaces after commas.
"""
unspace_commas(text)::String = replace_multi(text, ", " => ",", ". " => ".")

"""
Split a string into a vector of Chars.
"""
chars(text::AbstractString)::Vector{Char} = collect(text)

"""
    join_lines(lines)

Merge a vector of strings in a single string.
"""
join_lines(lines::Vector{String})::String = join(lines, "\n")
join_lines(lines::Vector)::String = join(lines, "\n")
join_lines(lines...) = join(lines, "\n")

"""
    split_lines(text::AbstractString)

Split a string into its composing lines.
"""
split_lines(text::String)::Vector{String} = split(text, "\n")
split_lines(text::SubString)::Vector{String} = String.(split(text, "\n"))

"""
    split_lines(renderable)

Split a renderable's text.
"""
function split_lines(renderable)
    string(typeof(renderable)) == "Segment" && return split_lines(renderable.text)
    return [s.text for s in renderable.segments]

end

"""
    do_by_line(fn::Function, text::String)

Apply `fn` to each line in the `text`.

The function `fn` should accept a single `::String` argument.
"""
function do_by_line(fn::Function, text)::String
    out = ""
    for (last, line) in loop_last(split_lines(text))
        out *= fn(line) * (last ? "" : "\n")
    end
    return out
end

do_by_line(fn::Function, text::Vector)::String = join_lines(fn.(text))

# ------------------------------- reshape text ------------------------------- #
"""
    fillin(text::String)::String

Ensure that each line in a multi-line text has the same width.
"""
function fillin(text)::String
    lines = split_lines(text)
    length(lines) == 1 && return text

    w = max(map(textlen, lines)...)
    return join_lines(map(
        (ln) -> ln * " "^(w - textlen(ln)),
        lines
    ))
end

"""
    truncate(text::AbstractString, width::Int)

Shorten a string of text to a target width
"""
function truncate(text::AbstractString, width::Int)
    textlen(text) <= width && return text
    return text[1:prevind(text, width - 2)] * "..."
end

# ---------------------------------------------------------------------------- #
#                                 RESHAPE TEXT                                 #
# ---------------------------------------------------------------------------- #

"""
Store information about the location and style of an ansi tag
"""
struct AnsiTag
    start::Int
    stop::Int
    style::SubString
end

"""
    excise_style(text)::Tuple{String, Vector{AnsiTag}, Bool}

Cut ANSI style information out of a string and store ANSI tags locations.
"""
function excise_style(text)::Tuple{String, Vector{AnsiTag}, Bool}
    hasstyle = has_markup(text) || has_ansi(text)
    text = apply_style(text)

    tags::Vector{AnsiTag} = []
    while occursin(ANSI_REGEXE, text)
        mtch = match(ANSI_REGEXE, text)
        _start = thisind(text, max(mtch.offset-1, 1))
        _stop = thisind(text, mtch.offset + length(mtch.match))

        push!(tags, AnsiTag(_start, _stop, mtch.match))
        text = text[1:_start] * text[_stop:end]
    end
    return text, tags, hasstyle
end

"""
    reinject_style(text, tags::Vector{AnsiTag}, cuts::Vector{Int})

Re-insert previously excised ANSI style tags.
If the cleaned text was reshaped, `cuts` stores information about
where new lines were added so that the ANSI style can be added
in the right place.
"""
function reinject_style(text, tags::Vector{AnsiTag}, cuts::Vector{Int}, issimple::Bool)
    for tag in reverse(tags)
        pads = findall(cuts .< tag.start)

        if issimple
            pads = isnothing(pads) ? 0 : max(0, length(pads)-1)

            cut1 = min(tag.start + pads, length(text))
            cut1 = text[cut1] == '\n' ? cut1+1 : cut1
            cut2 = cut1+1

        else
            pads = isnothing(pads) ? 0 : max(0, length(pads))

            cut1 = thisind(text, tag.start+pads)
            cut2 = thisind(text, cut1 + ncodeunits(text[cut1]))
        end
        cut1 > ncodeunits(text) && continue
        # cut2 = min(cut2, ncodeunits(text))

        text = text[1:cut1] * tag.style * text[cut2:end]

    end
    return text
end

"""
    correct_newline_style(text)

Given a multi-line text with ANSI style, add an ANSI close
tag at the end of each line and restore the style in the next.

This is important for reshaped text that needs to be put in
layouts with other text and you don't want the ANSI style
to bleed into other renderables.
"""
function correct_newline_style(text)
    ansi = ""
    out = ""
    for (last, line) in loop_last(split(text, "\n"))
        out *= ansi * line * (last ? "\e[0m" : "\e[0m\n")
        ansi = last ? "" : get_ANSI_codes(line)
    end
    out
end

"""
    get_text_info(text)

Extract a bunch of information from a string of text.
Info includes with at each char, number of codeunits at each char,
which characters are spaces and wether the text is 'simple' (i.e.
each char has ncodeunits=1) or not.
"""
function get_text_info(text)
    ctext = collect(text)
    widths = cumsum(textwidth.(ctext))
    nunits = cumsum(ncodeunits.(ctext))

    issimple = length(text) == ncodeunits(text)
    if issimple
        isspace = ' ' .== ctext
        nchar = identity
    else
        charidxs = collect(eachindex(text))  # codeunits idx of start of each char
        nchar(unit) = or(findfirst(charidxs .== prevind(text, unit)), 1)  # n chars at codeunit
        isspace = zeros(Bool, length(text))
        isspace[nchar.(only.(findall(" ", text)))] .= 1  # only.(findall(" ", text) can become findall(' ', text) for Julia > 1.6
    end

    return widths, isspace, issimple, nunits, nchar
end

"""
    reshape_text(text, width::Int)

Reshape a string to have max width when printed out.

Cut the string into multiple lines so that each line
has at most `width` when printed to terminal: ignoring
ANSI style but taking into account characters with
widths > 1. Text is preferentially split at a space
when possible to improve readability.

Future developer, this function is highly optimized to
reshape text as accurately as possible and as fast as possible.
It took several days to get it to be able to handle
style information correctly, handle multi codeunit
characters, handle characters with width >1, handle
situations in which the string can't be split up
at a space etc... It's a very delicate balance to get
everything right and changing anything can cause
the whole thing to break, so do so at your peril!
"""
function reshape_text(text, width::Int)
    if occursin("\n", text)
        return do_by_line(l -> reshape_text(l, width), text)
    end

    textlen(text) <= width && return text

    # Remove style information to simplify reshaping
    text, tags, hasstyle = excise_style(text)

    # extract style information
    original_widths, isspace, issimple, nunits, nchar = get_text_info(text)
    widths = view(original_widths, :) .- 1

    # infer places where to cut the text
    cuts = [1]
    N = length(original_widths)
    while cuts[end] < ncodeunits(text)
        # start with the worst case: we cut mid-word
        lastcut = cuts[end]
        candidate = issimple ? findlast(widths .<= width) : findlast(widths .< width)

        if isspace[candidate] == 1
            # we got lucky
            cut = candidate
        else
            # try to get a space in the current line and cut there if possible
            lastspace = findlast(view(isspace, lastcut:candidate))
            cut = isnothing(lastspace) ? candidate : lastspace + lastcut
        end

        if cut <= N
            # adjust widths so that we ignore already cut text
            push!(cuts, nunits[cut])
            widths = original_widths .- original_widths[cut]
        else
            # done
            break
        end
    end

    # create output text by slicing up the input text
    out = ""
    applied_cuts::Vector{Int} = []
    for (last, (pre, post)) in loop_last(zip(cuts[1:end-1], cuts[2:end]))
        post - pre <= 1 && continue
        Δ = last ? 0 : 1

        if issimple
            # simple means all characters have 1 codeunit length, easy
            !isspace[pre] && push!(applied_cuts, pre)
            newline = tview(text, pre, post-Δ, :simple)
        else
            # some chars have >1 codeunit length, we need to be careful
            _pre = thisind(text, pre)
            _post = thisind(text, post)
            newline = tview(text, _pre, _post, :simple)
            !isspace[nchar(pre)] && push!(applied_cuts, _pre)
        end

        # append to output text
        out *= lstrip(newline) * (last ? "" : "\n")
    end

    # if necessary re-insert style information in the reshaped text
    if hasstyle
        return correct_newline_style(reinject_style(out, tags, applied_cuts, issimple))
    else
        return out
    end
end