text:
function compress_mpo!(env, mpstarget, mpo, maxsweeps)

    mps = MPO_to_MPS(mpo)

    compress_var!(mpstarget, mps, env, maxsweeps)

    MPS_to_MPO(mpstarget)

end

"""
function to compress mps A using maxsweeps left and right of the
variational algorithm with Ain as initial guess. Ain must be in left
canonical form (normalisation not necessary)
"""

function compress_var(Ain,A,maxsweeps)

    na = length(Ain)
    FL = Array{Array{Complex{Float64}, 2}, 1}(na)
    FR = Array{Array{Complex{Float64}, 2}, 1}(na)
    A0 = Array{Array{Complex{Float64}, 3}, 1}(na)

    #initializing left environment

    FL[1] = eye(1)
    for jj = 1:(na - 1)
        FL[jj + 1] = update_lenv(Ain[jj], A[jj], FL[jj])
    end
    FR[na] = eye(1)

    for sweep = 1:maxsweeps

        #sweep left

        for jj = na:-1:2
            AR = prod_LR(A[jj], FR[jj])
            A0[jj] = prod_LR(FL[jj], AR)
            A0[jj] = rightorth!(A0[jj])
            FR[jj - 1] = update_renv(A0[jj], AR)
        end
        AR = prod_LR(A[1], FR[1])
        A0[1] = prod_LR(FL[1], AR)

        #sweep right

        for jj = 1:(na - 1)
            AL = prod_LR(FL[jj], A[jj])
            A0[jj] = prod_LR(AL, FR[jj])
            A0[jj] = leftorth!(A0[jj])
            FL[jj + 1] = update_lenv(A0[jj], AL)
        end
        AL = prod_LR(FL[na], A[na])
        A0[na] = prod_LR(AL, FR[na])

    end

    return A0

end

function compress_var!(A0, A, env, maxsweeps)

    na = length(A0)

    #initializing left environment

    env[1][1] = 1.0
    for jj = 1:(na - 1)
        update_lenv_co!(env[jj + 1], A0[jj], A[jj], env[jj])
    end
    env[na + 1][1] = 1.0

    for sweep = 1:maxsweeps

        #sweep left

        for jj = na:-1:2
            AR = prod_LR(A[jj], env[jj + 1])
            prod_LR_co!(A0[jj], env[jj], AR)
            A0[jj] = rightorth!(A0[jj])
            update_renv!(env[jj], A0[jj], AR)
        end
        AR = prod_LR(A[1], env[2])
        prod_LR!(A0[1], env[1], AR)

        #sweep right

        for jj = 1:(na - 1)
            AL = prod_LR_co(env[jj], A[jj])
            prod_LR!(A0[jj], AL, env[jj + 1])
            A0[jj] = leftorth!(A0[jj])
            update_lenv_co!(env[jj + 1], A0[jj], AL)
        end
        AL = prod_LR_co(env[na], A[na])
        prod_LR!(A0[na], AL, env[na + 1])

    end

end
"""
function to compress mps A using maxsweeps left and right of the two
site variational algorithm with Ain as initial guess. Ain must be in left
canonical form (normalisation not necessary)
"""

function compress_var_two_site(Ain,A,maxsweeps,svd_tol,dims,dlim)

    na = length(Ain)
    FL = Array{Array{Complex{Float64},2},1}(na)
    FR = Array{Array{Complex{Float64},2},1}(na)
    A0 = Array{Array{Complex{Float64},3},1}(na)

    #initializing left environment

    FL[1] = eye(1)
    for jj = 1:na-1
        FL[jj+1] = update_lenv(Ain[jj],A[jj],FL[jj])
    end
    FR[na] = eye(1);

    for sweep = 1:maxsweeps

        U,S,Vd = eye(1),eye(1),eye(1)
        #sweep left

        for jj = na:-1:2;

            av = A_prod_LR_two_site(FL[jj-1],A[jj-1],A[jj],FR[jj])
            U,S,Vd = svdecon(reshape(av,dims[jj-1]*size(A[jj-1],2),
                size(A[jj],2)*dims[jj+1]))
            S = S[S.>0]
            stot = sum(S.^2)
            dims[jj] = min(sum(cumsum(S.^2/stot) .< 1-svd_tol)+1,dlim)
            A0[jj] = reshape(Vd[:,1:dims[jj]]',dims[jj],size(A[jj],2),dims[jj+1])
            FR[jj-1] = update_renv(A0[jj],A[jj],FR[jj])

        end
        #A0[1] = reshape(U[:,1:Dims[2]]*diagm(S[1:Dims[2]]),Dims[1],size(A[1],2),Dims[2])

        #sweep right

        for jj = 1:na-1;

            av = A_prod_LR_two_site(FL[jj],A[jj],A[jj+1],FR[jj+1])
            U,S,Vd = svdecon(reshape(av,dims[jj]*size(A[jj],2),
                size(A[jj+1],2)*dims[jj+2]))
            S = S[S.>0]
            stot = sum(S.^2)
            dims[jj+1] = min(sum(cumsum(S.^2/stot) .< 1-svd_tol)+1,dlim)
            A0[jj] = reshape(U[:,1:dims[jj+1]],dims[jj],size(A[jj],2),dims[jj+1])
            FL[jj+1] = update_lenv(A0[jj],A[jj],FL[jj])

        end
        A0[na] = reshape(diagm(S[1:dims[na]])*(Vd[:,1:dims[na]]'),dims[na],size(A[na],2),dims[na+1])

    end

    return A0

end

"""
function to compress mps resulting from the product of H and A using
maxsweeps left and right of the variational algorithm with Ain as
initial guess. Ain must be in left canonical form (normalisation not necessary)
"""

function compress_var_apply_H(Ain, A, H, maxsweeps)

    na = length(Ain)
    FL = Array{Array{Complex{Float64}, 3}, 1}(na)
    FR = Array{Array{Complex{Float64}, 3}, 1}(na)
    A0 = Array{Array{Complex{Float64}, 3}, 1}(na)

    #initializing left environment

    FL[1] = ones(1, 1, 1)
    for jj = 1:(na - 1)
        FL[jj + 1] = update_le(Ain[jj], H[jj], A[jj], FL[jj])
    end
    FR[na] = ones(1, 1, 1)

    for sweep = 1:maxsweeps

        #sweep left

        for jj = na:-1:2;
            AOR = prod_AOR(A[jj], H[jj], FR[jj])
            A0[jj] = prod_LR(FL[jj], AOR)
            A0[jj] = rightorth!(A0[jj])
            FR[jj - 1] = update_re(A0[jj], AOR)
        end

        #sweep right

        for jj = 1:(na - 1)
            AOL = prod_AOL(A[jj], H[jj], FL[jj])
            A0[jj] = prod_LR(AOL, FR[jj])
            A0[jj] = leftorth!(A0[jj])
            FL[jj + 1] = update_le(A0[jj], AOL)
        end

        AOL = prod_AOL(A[na], H[na], FL[na])
        A0[na] = prod_LR(AOL, FR[na])

    end

    return A0

end

function compress_var_apply_H!(A0, envop, Ain, A, H, maxsweeps)

    na = length(Ain)

    #initializing left environment

    envop[1][1] = 1.0
    for jj = 1:(na - 1)
        update_le!(envop[jj + 1], Ain[jj], H[jj], A[jj], envop[jj])
    end
    envop[na + 1][1] = 1.0

    for sweep = 1:maxsweeps

        #sweep left

        for jj = na:-1:2

            AOR = prod_AOR(A[jj], H[jj], envop[jj + 1])
            prod_LR!(A0[jj], envop[jj], AOR)
            A0[jj] = rightorth!(A0[jj])
            update_re!(envop[jj], A0[jj], AOR)

        end

        #sweep right

        for jj = 1:(na - 1)

            AOL = prod_AOL(A[jj], H[jj], envop[jj])
            prod_LR!(A0[jj], AOL, envop[jj + 1])
            A0[jj] = leftorth!(A0[jj])
            update_le!(envop[jj + 1], A0[jj], AOL)

        end

        AOL = prod_AOL(A[na], H[na], envop[na])
        prod_LR!(A0[na], AOL, envop[na + 1])

    end

end

function compress_var_apply_H_right(Ain, A, H, maxsweeps)

    na = length(Ain)
    FR = Array{Array{Complex{Float64}, 3}, 1}(na)
    A0 = Array{Array{Complex{Float64}, 3}, 1}(na)

    #initializing right environment

    FR[na] = ones(1, 1, 1)
    for jj = na:-1:2;
        FR[jj - 1] = update_re(Ain[jj], H[jj], A[jj], FR[jj])
    end
    FL = complex(ones(1, 1, 1))

    #sweep right

    for jj = 1:(na - 1)
        AOL = prod_AOL(A[jj], H[jj], FL)
        A0[jj] = prod_LR(AOL, FR[jj])
        A0[jj] = leftorth!(A0[jj])
        FL = update_le(A0[jj], AOL)
    end

    AOL = prod_AOL(A[na], H[na], FL)
    A0[na] = prod_LR(AOL, FR[na])

    return A0

end

function compress_var_apply_H_left(Ain, A, H)

    na = length(Ain)
    FL = Array{Array{Complex{Float64}, 3}, 1}(na)
    A0 = Array{Array{Complex{Float64}, 3}, 1}(na)

    #initializing left environment

    FL[1] = ones(1, 1, 1)
    for jj = 1:(na - 1)
        FL[jj + 1] = update_le(Ain[jj], H[jj], A[jj], FL[jj])
    end
    FR = complex(ones(1, 1, 1))

    #sweep left

    for jj = na:-1:2;
        AOR = prod_AOR(A[jj], H[jj], FR)
        A0[jj] = prod_LR(FL[jj], AOR)
        A0[jj] = rightorth!(A0[jj])
        FR = update_re(A0[jj], AOR)
    end

    AOR = prod_AOR(A[1], H[1], FR)
    A0[1] = prod_LR(FL[1], AOR)

    return A0

end

function compress_var_apply_H_left!(A0, envop, Ain, A, H)

    na = length(Ain)

    #initializing left environment

    envop[1][1] = 1.0
    for jj = 1:(na - 1)
        update_le!(envop[jj + 1], Ain[jj], H[jj], A[jj], envop[jj])
    end
    envop[na + 1][1] = 1.0

    #sweep left

    for jj = na:-1:2

        AOR = prod_AOR(A[jj], H[jj], envop[jj + 1])
        prod_LR!(A0[jj], envop[jj], AOR)
        A0[jj] = rightorth!(A0[jj])
        update_re!(envop[jj], A0[jj], AOR)

    end

    AOR = prod_AOR(A[1], H[1], envop[2])
    prod_LR!(A0[1], envop[1], AOR)

end

"""
function to compress mps resulting from the product of H and A using
maxsweeps left and right of the variational two site algorithm with Ain as
initial guess. Ain must be in left canonical form (normalisation not necessary)
"""

function compress_var_apply_H_two_site(Ain,A,H,maxsweeps,svd_tol,dims,dlim)

    na = length(Ain)
    FL = Array{Array{Complex{Float64},3},1}(na)
    FR = Array{Array{Complex{Float64},3},1}(na)
    A0 = Array{Array{Complex{Float64},3},1}(na)

    #initializing left environment

    FL[1] = ones(1,1,1)
    for jj = 1:na-2
        FL[jj+1] = update_le(Ain[jj],H[jj],A[jj],FL[jj])
    end
    FR[na] = ones(1,1,1);

    for sweep = 1:maxsweeps

        U,S,Vd = eye(1),eye(1),eye(1)
        #sweep left

        for jj = na:-1:2;

            av = A_prod_LOR_two_site(FL[jj-1],A[jj-1],H[jj-1],A[jj],H[jj],FR[jj])
            U,S,Vd = svdecon(reshape(av,dims[jj-1]*size(A[jj-1],2),
                size(A[jj],2)*dims[jj+1]))
            S = S[S.>0]
            stot = sum(S.^2)
            dims[jj] = min(sum(cumsum(S.^2/stot) .< 1-svd_tol)+1,dlim)
            A0[jj] = reshape(Vd[:,1:dims[jj]]',dims[jj],size(A[jj],2),dims[jj+1])
            FR[jj-1] = update_re(A0[jj],H[jj],A[jj],FR[jj])

        end

        #A0[1] = reshape(U[:,1:Dims[2]]*diagm(S[1:Dims[2]]),Dims[1],size(A[1],2),Dims[2])

        #sweep right

        for jj = 1:na-1;

            av = A_prod_LOR_two_site(FL[jj],A[jj],H[jj],A[jj+1],H[jj+1],FR[jj+1])
            U,S,Vd = svdecon(reshape(av,dims[jj]*size(A[jj],2),
                size(A[jj+1],2)*dims[jj+2]))
            S = S[S.>0]
            stot = sum(S.^2)
            dims[jj+1] = min(sum(cumsum(S.^2/stot) .< 1-svd_tol)+1,dlim)
            A0[jj] = reshape(U[:,1:dims[jj+1]],dims[jj],size(A[jj],2),dims[jj+1])
            FL[jj+1] = update_le(A0[jj],H[jj],A[jj],FL[jj])

        end

        A0[na] = reshape(diagm(S[1:dims[na]])*(Vd[:,1:dims[na]]'),dims[na],size(A[na],2),dims[na+1])

    end

    return A0

end

function compress_var_apply_H_two_site_left(Ain,A,H,maxsweeps,svd_tol,dims,dlim)

    na = length(Ain)
    FL = Array{Array{Complex{Float64},3},1}(na)
    FR = Array{Array{Complex{Float64},3},1}(na)
    A0 = Array{Array{Complex{Float64},3},1}(na)

    #initializing left environment

    FL[1] = ones(1,1,1)
    for jj = 1:na-2
        FL[jj+1] = update_le(Ain[jj],H[jj],A[jj],FL[jj])
    end
    FR[na] = ones(1,1,1);

    U,S,Vd = eye(1),eye(1),eye(1)
    #sweep left

    for jj = na:-1:2;

        av = A_prod_LOR_two_site(FL[jj-1],A[jj-1],H[jj-1],A[jj],H[jj],FR[jj])
        U,S,Vd = svdecon(reshape(av,dims[jj-1]*size(A[jj-1],2),
            size(A[jj],2)*dims[jj+1]))
        S = S[S.>0]
        stot = sum(S.^2)
        dims[jj] = min(sum(cumsum(S.^2/stot) .< 1-svd_tol)+1,dlim)
        A0[jj] = reshape(Vd[:,1:dims[jj]]',dims[jj],size(A[jj],2),dims[jj+1])
        FR[jj-1] = update_re(A0[jj],H[jj],A[jj],FR[jj])

    end

    A0[1] = reshape(U[:,1:dims[2]]*diagm(S[1:dims[2]]),dims[1],size(A[1],2),dims[2])

    return A0

end

function compress_var_apply_H_two_site_right(Ain,A,H,maxsweeps,svd_tol,dims,dlim)

    na = length(Ain)
    FL = Array{Array{Complex{Float64},3},1}(na)
    FR = Array{Array{Complex{Float64},3},1}(na)
    A0 = Array{Array{Complex{Float64},3},1}(na)

    #initializing right environment

    FR[na] = ones(1,1,1);
    for jj = na:-1:2;
        FR[jj-1] = update_re(Ain[jj],H[jj],A[jj],FR[jj])
    end
    FL[1] = ones(1,1,1);

    U,S,Vd = eye(1),eye(1),eye(1)

    #sweep right

    for jj = 1:na-1;

        av = A_prod_LOR_two_site(FL[jj],A[jj],H[jj],A[jj+1],H[jj+1],FR[jj+1])
        U,S,Vd = svdecon(reshape(av,dims[jj]*size(A[jj],2),
            size(A[jj+1],2)*dims[jj+2]))
        S = S[S.>0]
        stot = sum(S.^2)
        dims[jj+1] = min(sum(cumsum(S.^2/stot) .< 1-svd_tol)+1,dlim)
        A0[jj] = reshape(U[:,1:dims[jj+1]],dims[jj],size(A[jj],2),dims[jj+1])
        FL[jj+1] = update_le(A0[jj],H[jj],A[jj],FL[jj])

    end

    A0[na] = reshape(diagm(S[1:dims[na]])*(Vd[:,1:dims[na]]'),dims[na],size(A[na],2),dims[na+1])

    return A0

end

"""
function to compress linear combination of mps stored in array A, which
have weights given by coef, using a
variational algorithm with Ain as initial guess. Ain must be in left
canonical form (normalisation not necessary)
"""

function compress_sum_var(Ain, coef, A, maxsweeps)

    na = length(Ain)
    numv = length(coef)
    FL = Array{Array{Array{Complex{Float64}, 2}, 1}, 1}(numv)
    FR = Array{Array{Array{Complex{Float64}, 2}, 1}, 1}(numv)
    A0 = Array{Array{Complex{Float64}, 3}, 1}(na)
    AR = Array{Array{Complex{Float64}, 3}, 1}(numv)

    #initializing left environment

    for ll = 1:numv
        #A[ll][na] = coef[ll]*Av[ll][na]
        FL[ll] = Array{Array{Complex{Float64}, 2}, 1}(na)
        FL[ll][1] = eye(1)
        for jj = 1:(na - 1)
            FL[ll][jj + 1] = update_lenv(Ain[jj], A[ll][jj], FL[ll][jj])
        end
        FR[ll] = Array{Array{Complex{Float64}, 2}, 1}(na)
        FR[ll][na] = coef[ll]*eye(1)
    end

    for sweep = 1:maxsweeps

        #sweep left

        for jj = na:-1:2

            AR[1] = prod_LR(A[1][jj], FR[1][jj])
            A0[jj] = prod_LR(FL[1][jj], AR[1])
            for ll = 2:numv
                AR[ll] = prod_LR(A[ll][jj], FR[ll][jj])
                A0[jj] += prod_LR(FL[ll][jj], AR[ll])
            end
            A0[jj] = rightorth!(A0[jj])
            for ll = 1:numv
                FR[ll][jj - 1] = update_renv(A0[jj], AR[ll])
            end

        end

        #sweep right

        for jj = 1:(na - 1)

            AR[1] = prod_LR(FL[1][jj], A[1][jj])
            A0[jj] = prod_LR(AR[1], FR[1][jj])
            for ll = 2:numv
                AR[ll] = prod_LR(FL[ll][jj], A[ll][jj])
                A0[jj] += prod_LR(AR[ll], FR[ll][jj])
            end
            A0[jj] = leftorth!(A0[jj])
            for ll = 1:numv
                FL[ll][jj + 1] = update_lenv(A0[jj], AR[ll])
            end

        end

    end

    AR[1] = prod_LR(FL[1][na], A[1][na])
    A0[na] = prod_LR(AR[1], FR[1][na])
    for ll = 2:numv
        AR[ll] = prod_LR(FL[ll][na], A[ll][na])
        A0[na] += prod_LR(AR[ll], FR[ll][na])
    end

    return A0

end

"""
function to compress linear combination of mps stored array A, which
have weights given by coef, using the two site
variational algorithm with Ain as initial guess. Ain must be in left
canonical form (normalisation not necessary)
"""

function compress_sum_var_two_site(Ain,coef,A,maxsweeps,svd_tol,dims,dlim)

    na = length(Ain)
    numv = length(coef)
    FL = Array{Array{Array{Complex{Float64},2},1},1}(numv)
    FR = Array{Array{Array{Complex{Float64},2},1},1}(numv)
    A0 = Array{Array{Complex{Float64},3},1}(na)

    #initializing left environment

    for ll = 1:numv

        #A[ll][na] = coef[ll]*Av[ll][na]
        FL[ll] = Array{Array{Complex{Float64},2},1}(na)
        FL[ll][1] = eye(1)
        for jj = 1:na-2
            FL[ll][jj+1] = update_lenv(Ain[jj],A[ll][jj],FL[ll][jj])
        end
        FR[ll] = Array{Array{Complex{Float64},2},1}(na)
        FR[ll][na] = coef[ll]*eye(1);

    end

    for sweep = 1:maxsweeps

        U,S,Vd = eye(1),eye(1),eye(1)
        #sweep left

        for jj = na:-1:2

            av = A_prod_LR_two_site(FL[1][jj-1],A[1][jj-1],A[1][jj],FR[1][jj])
            for ll = 2:numv
                av += A_prod_LR_two_site(FL[ll][jj-1],A[ll][jj-1],A[ll][jj],FR[ll][jj])
            end
            U,S,Vd = svdecon(reshape(av,dims[jj-1]*size(A[1][jj-1],2),
                size(A[1][jj],2)*dims[jj+1]))
            S = S[S.>0]
            stot = sum(S.^2)
            dims[jj] = min(sum(cumsum(S.^2/stot) .< 1-svd_tol)+1,dlim)
            A0[jj] = reshape(Vd[:,1:dims[jj]]',dims[jj],size(A[1][jj],2),dims[jj+1])
            for ll = 1:numv
                FR[ll][jj-1] = update_renv(A0[jj],A[ll][jj],FR[ll][jj])
            end

        end

        #A0[1] = reshape(U[:,1:Dims[2]]*diagm(S[1:Dims[2]]),Dims[1],size(A[1][1],2),Dims[2])

        #sweep right

        for jj = 1:na-1;

            av = A_prod_LR_two_site(FL[1][jj],A[1][jj],A[1][jj+1],FR[1][jj+1])
            for ll = 2:numv
                av += A_prod_LR_two_site(FL[ll][jj],A[ll][jj],A[ll][jj+1],FR[ll][jj+1])
            end
            U,S,Vd = svdecon(reshape(av,dims[jj]*size(A[1][jj],2),
                size(A[1][jj+1],2)*dims[jj+2]))
            S = S[S.>0]
            stot = sum(S.^2)
            dims[jj+1] = min(sum(cumsum(S.^2/stot) .< 1-svd_tol)+1,dlim)
            A0[jj] = reshape(U[:,1:dims[jj+1]],dims[jj],size(A[1][jj],2),dims[jj+1])
            for ll = 1:numv
                FL[ll][jj+1] = update_lenv(A0[jj],A[ll][jj],FL[ll][jj])
            end
        end

        A0[na] = reshape(diagm(S[1:dims[na]])*(Vd[:,1:dims[na]]'),dims[na],size(A[1][na],2),dims[na+1])

    end

    return A0

end

"""
function to compress linear combination of mps stored in cell array A, which
have weights given by coef and where H is applied to the last MPS.
The compression uses the variational algorithm with A0 as initial guess.
A0 must be in left canonical form (normalisation not necessary)
"""

function compress_sum_var_apply_H(Ain, coef, A, H, maxsweeps)

    na = length(Ain)
    numv = length(coef)
    FL = Array{Array{Array{Complex{Float64}, 2}, 1}, 1}(numv - 1)
    FR = Array{Array{Array{Complex{Float64}, 2}, 1}, 1}(numv - 1)
    FLH = Array{Array{Complex{Float64}, 3}, 1}(na)
    FRH = Array{Array{Complex{Float64} ,3}, 1}(na)
    A0 = Array{Array{Complex{Float64}, 3}, 1}(na)
    AR = Array{Array{Complex{Float64}, 3}, 1}(numv - 1)

    #initializing left environment

    for ll = 1:(numv - 1)
        #A[ll][na] = coef[ll]*Av[ll][na]
        FL[ll] = Array{Array{Complex{Float64}, 2}, 1}(na)
        FL[ll][1] = eye(1)
        for jj = 1:(na - 1)
            FL[ll][jj + 1] = update_lenv(Ain[jj], A[ll][jj], FL[ll][jj])
        end
        FR[ll] = Array{Array{Complex{Float64}, 2}, 1}(na)
        FR[ll][na] = coef[ll]*eye(1)
    end

    FLH[1] = ones(1, 1, 1)
    for jj = 1:(na - 1)
        FLH[jj + 1] = update_le(Ain[jj], H[jj], A[numv][jj], FLH[jj])
    end
    FRH[na] = coef[numv]*ones(1, 1, 1)

    for sweep = 1:maxsweeps

        #sweep left

        for jj = na:-1:2;

            AR[1] = prod_LR(A[1][jj], FR[1][jj])
            A0[jj] = prod_LR(FL[1][jj], AR[1])
            for ll = 2:(numv - 1)
                AR[ll] = prod_LR(A[ll][jj], FR[ll][jj])
                A0[jj] += prod_LR(FL[ll][jj], AR[ll])
            end
            AOR = prod_AOR(A[numv][jj], H[jj], FRH[jj])
            A0[jj] += prod_LR(FLH[jj], AOR)
            A0[jj] = rightorth!(A0[jj])
            for ll = 1:(numv - 1)
                FR[ll][jj - 1] = update_renv(A0[jj], AR[ll])
            end
            FRH[jj - 1] = update_re(A0[jj], AOR)
        end

        #sweep right

        for jj = 1:(na - 1)

            AR[1] = prod_LR(FL[1][jj], A[1][jj])
            A0[jj] = prod_LR(AR[1], FR[1][jj])
            for ll = 2:(numv - 1)
                AR[ll] = prod_LR(FL[ll][jj], A[ll][jj])
                A0[jj] += prod_LR(AR[ll], FR[ll][jj])
            end
            AOL = prod_AOL(A[numv][jj], H[jj], FLH[jj])
            A0[jj] += prod_LR(AOL, FRH[jj])
            A0[jj] = leftorth!(A0[jj])
            for ll = 1:(numv - 1)
                FL[ll][jj + 1] = update_lenv(A0[jj], AR[ll])
            end
            FLH[jj + 1] = update_le(A0[jj], AOL)
        end

    end

    AR[1] = prod_LR(FL[1][na], A[1][na])
    A0[na] = prod_LR(AR[1], FR[1][na])
    for ll = 2:numv-1
        AR[ll] = prod_LR(FL[ll][na], A[ll][na])
        A0[na] += prod_LR(AR[ll], FR[ll][na])
    end
    AOL = prod_AOL(A[numv][na], H[na], FLH[na])
    A0[na] += prod_LR(AOL, FRH[na])

    return A0

end

function compress_sum_var_apply_H_left(Ain, coef, A, H)

    na = length(Ain)
    numv = length(coef)
    FL = Array{Array{Array{Complex{Float64}, 2}, 1}, 1}(numv - 1)
    FR = Array{Array{Complex{Float64}, 2}, 1}(numv - 1)
    FLH = Array{Array{Complex{Float64}, 3}, 1}(na)
    A0 = Array{Array{Complex{Float64}, 3}, 1}(na)
    AR = Array{Array{Complex{Float64}, 3}, 1}(numv - 1)

    #initializing left environment

    for ll = 1:(numv - 1)
        #A[ll][na] = coef[ll]*Av[ll][na]
        FL[ll] = Array{Array{Complex{Float64}, 2}, 1}(na)
        FL[ll][1] = eye(1)
        for jj = 1:(na - 1)
            FL[ll][jj + 1] = update_lenv(Ain[jj], A[ll][jj], FL[ll][jj])
        end
        FR[ll] = coef[ll]*eye(1)

    end

    FLH[1] = ones(1, 1, 1)
    for jj = 1:(na - 1)
        FLH[jj + 1] = update_le(Ain[jj], H[jj], A[numv][jj], FLH[jj])
    end

    FRH = coef[numv]*complex(ones(1, 1, 1))

    #sweep left

    for jj = na:-1:2;

        AR[1] = prod_LR(A[1][jj], FR[1])
        A0[jj] = prod_LR(FL[1][jj], AR[1])
        for ll = 2:(numv - 1)
            AR[ll] = prod_LR(A[ll][jj], FR[ll])
            A0[jj] += prod_LR(FL[ll][jj], AR[ll])
        end
        AOR = prod_AOR(A[numv][jj], H[jj], FRH)
        A0[jj] += prod_LR(FLH[jj], AOR)
        A0[jj] = rightorth!(A0[jj])
        for ll = 1:(numv - 1)
            FR[ll] = update_renv(A0[jj], AR[ll])
        end
        FRH = update_re(A0[jj], AOR)
    end

    AR[1] = prod_LR(A[1][1], FR[1])
    A0[1] = prod_LR(FL[1][1], AR[1])
    for ll = 2:numv-1
        AR[ll] = prod_LR(A[ll][1], FR[ll])
        A0[1] += prod_LR(FL[ll][1], AR[ll])
    end
    AOR = prod_AOR(A[numv][1], H[1], FRH)
    A0[1] += prod_LR(FLH[1], AOR)

    return A0

end

function compress_sum_var_apply_H!(A0, envop, coef, A, H, maxsweeps)

    na = length(A0)
    numv = length(coef)
    env = Array{Array{Array{Complex{Float64}, 2}, 1}, 1}(numv - 1)
    AR = Array{Array{Complex{Float64}, 3}, 1}(numv - 1)

    #initializing left environment

    for ll = 1:(numv - 1)
        #A[ll][na] = coef[ll]*Av[ll][na]
        env[ll] = Array{Array{Complex{Float64}, 2}, 1}(na + 1)
        env[ll][1][1] = 1.0
        for jj = 1:(na - 1)
            env[ll][jj + 1] = update_lenv(A0[jj], A[ll][jj], env[ll][jj])
        end
        env[ll][na + 1][1] = coef[ll]
    end

    envop[1][1] = 1.0
    for jj = 1:(na - 1)

        update_le!(envop[jj + 1], A0[jj], H[jj], A[numv][jj], envop[jj])

    end
    envop[na + 1][1] = coef[numv]

    for sweep = 1:maxsweeps

        #sweep left

        for jj = na:-1:2;

            AR[1] = prod_LR(A[1][jj], env[1][jj + 1])
            prod_LR!(A0[jj], env[1][jj], AR[1])
            for ll = 2:(numv - 1)
                AR[ll] = prod_LR(A[ll][jj], env[ll][jj + 1])
                A0[jj] += prod_LR(env[ll][jj], AR[ll])
            end
            AOR = prod_AOR(A[numv][jj], H[jj], envop[jj + 1])
            A0[jj] += prod_LR(envop[jj], AOR)
            A0[jj] = rightorth!(A0[jj])
            for ll = 1:(numv - 1)
                update_renv!(env[ll][jj], A0[jj], AR[ll])
            end
            update_re!(envop[jj], A0[jj], AOR)

        end

        #sweep right

        for jj = 1:(na - 1)

            AR[1] = prod_LR(env[1][jj], A[1][jj])
            prod_LR!(A0[jj], AR[1], env[1][jj + 1])
            for ll = 2:(numv - 1)
                AR[ll] = prod_LR(env[ll][jj], A[ll][jj])
                A0[jj] += prod_LR(AR[ll], env[ll][jj + 1])
            end
            AOL = prod_AOL(A[numv][jj], H[jj], envop[jj])
            A0[jj] += prod_LR(AOL, envop[jj + 1])
            A0[jj] = leftorth!(A0[jj])
            for ll = 1:(numv - 1)
                update_lenv!(env[ll][jj + 1], A0[jj], AR[ll])
            end
            update_le!(envop[jj + 1], A0[jj], AOL)

        end

    end

    AR[1] = prod_LR(env[1][na], A[1][na])
    prod_LR!(A0[na], AR[1], env[1][na + 1])
    for ll = 2:numv-1
        AR[ll] = prod_LR(env[ll][na], A[ll][na])
        A0[na] += prod_LR(AR[ll], env[ll][na + 1])
    end
    AOL = prod_AOL(A[numv][na], H[na], envop[na])
    A0[na] += prod_LR(AOL, envop[na + 1])

    nothing

end

function compress_sum_var_apply_H_left!(A0, env, envop, Ain, coef, A, H)

    na = length(Ain)
    numv = length(coef)
    AR = similar(A0, numv - 1)

    #initializing left environment

    for ll = 1:(numv - 1)

        env[ll][1][1] = 1.0
        for jj = 1:(na - 1)
            update_lenv!(env[ll][jj + 1], Ain[jj], A[ll][jj], env[ll][jj])
        end
        env[ll][na + 1][1] = coef[ll]

    end

    envop[1][1] = 1.0
    for jj = 1:(na - 1)

        update_le!(envop[jj + 1], Ain[jj], H[jj], A[numv][jj], envop[jj])

    end
    envop[na + 1][1] = coef[numv]

    #sweep left

    for jj = na:-1:2;

        AR[1] = prod_LR(A[1][jj], env[1][jj + 1])
        prod_LR!(A0[jj], env[1][jj], AR[1])
        for ll = 2:(numv - 1)
            AR[ll] = prod_LR(A[ll][jj], env[ll][jj + 1])
            A0[jj] += prod_LR(env[ll][jj], AR[ll])
        end
        AOR = prod_AOR(A[numv][jj], H[jj], envop[jj + 1])
        A0[jj] += prod_LR(envop[jj], AOR)
        A0[jj] = rightorth!(A0[jj])
        for ll = 1:(numv - 1)
            update_renv!(env[ll][jj], A0[jj], AR[ll])
        end
        update_re!(envop[jj], A0[jj], AOR)
    end

    AR[1] = prod_LR(A[1][1], env[1][2])
    prod_LR!(A0[1], env[1][1], AR[1])
    for ll = 2:numv-1
        AR[ll] = prod_LR(A[ll][1], env[ll][2])
        A0[1] += prod_LR(env[ll][1], AR[ll])
    end
    AOR = prod_AOR(A[numv][1], H[1], envop[2])
    A0[1] += prod_LR(envop[1], AOR)
    nothing

end

function compress_sum_var_apply_H_right(Ain, coef, A, H)

    na = length(Ain)
    numv = length(coef)
    FL = Array{Array{Complex{Float64}, 2}, 1}(numv - 1)
    FR = Array{Array{Array{Complex{Float64}, 2}, 1}, 1}(numv - 1)
    FRH = Array{Array{Complex{Float64}, 3}, 1}(na)
    A0 = Array{Array{Complex{Float64}, 3}, 1}(na)
    AL = Array{Array{Complex{Float64}, 3}, 1}(numv - 1)

    #initializing right environment

    for ll = 1:(numv - 1)
        FR[ll] = Array{Array{Complex{Float64}, 2}, 1}(na)
        FR[ll][na] = eye(1)
        for jj = na:-1:2
            FR[ll][jj - 1] = update_renv(Ain[jj], A[ll][jj], FR[ll][jj])
        end
        FL[ll] = coef[ll]*eye(1)
    end

    FRH[na] = ones(1, 1, 1)
    for jj = na:-1:2
        FRH[jj - 1] = update_re(Ain[jj], H[jj], A[numv][jj], FRH[jj])
    end

    FLH = coef[numv]*complex(ones(1, 1, 1))

    #sweep right

    for jj = 1:(na - 1)

        AL[1] = prod_LR(FL[1], A[1][jj])
        A0[jj] = prod_LR(AL[1], FR[1][jj])
        for ll = 2:(numv - 1)
            AL[ll] = prod_LR(FL[ll], A[ll][jj])
            A0[jj] += prod_LR(AL[ll], FR[ll][jj])
        end
        AOL = prod_AOL(A[numv][jj], H[jj], FLH)
        A0[jj] += prod_LR(AOL, FRH[jj])
        A0[jj] = leftorth!(A0[jj])
        for ll = 1:(numv - 1)
            FL[ll] = update_lenv(A0[jj], AL[ll])
        end
        FLH = update_le(A0[jj], AOL)

    end

    AL[1] = prod_LR(FL[1], A[1][na])
    A0[na] = prod_LR(AL[1], FR[1][na])
    for ll = 2:(numv - 1)
        AL[ll] = prod_LR(FL[ll], A[ll][na])
        A0[na] += prod_LR(AL[ll], FR[ll][na])
    end
    AOL = prod_AOL(A[numv][na], H[na], FLH)
    A0[na] += prod_LR(AOL, FRH[na])

    return A0

end

function compress_sum_var_apply_H_right!(A0, env, envop, Ain, coef, A, H)

    na = length(Ain)
    numv = length(coef)
    AL = similar(A0, numv - 1)

    #initializing right environment

    for ll = 1:(numv - 1)

        env[ll][na + 1][1] = 1.0
        for jj = na:-1:2
            update_renv!(env[ll][jj], Ain[jj], A[ll][jj], env[ll][jj + 1])
        end
        env[ll][1][1] = coef[ll]

    end

    envop[na + 1][1] = 1.0
    for jj = na:-1:2

        update_re!(envop[jj], Ain[jj], H[jj], A[numv][jj], envop[jj + 1])

    end
    envop[1][1] = coef[numv]

    #sweep right

    for jj = 1:(na - 1)

        AL[1] = prod_LR(env[1][jj], A[1][jj])
        prod_LR!(A0[jj], AL[1], env[1][jj + 1])
        for ll = 2:(numv - 1)
            AL[ll] = prod_LR(env[ll][jj], A[ll][jj])
            A0[jj] += prod_LR(AL[ll], env[ll][jj + 1])
        end
        AOL = prod_AOL(A[numv][jj], H[jj], envop[jj])
        A0[jj] += prod_LR(AOL, envop[jj + 1])
        A0[jj] = leftorth!(A0[jj])
        for ll = 1:(numv - 1)
            update_lenv!(env[ll][jj + 1], A0[jj], AL[ll])
        end
        update_le!(envop[jj + 1], A0[jj], AOL)

    end

    AL[1] = prod_LR(env[1][na], A[1][na])
    prod_LR!(A0[na], AL[1], env[1][na + 1])
    for ll = 2:(numv - 1)
        AL[ll] = prod_LR(env[ll][na], A[ll][na])
        A0[na] += prod_LR(AL[ll], env[ll][na + 1])
    end
    AOL = prod_AOL(A[numv][na], H[na], envop[na])
    A0[na] += prod_LR(AOL, envop[na + 1])
    nothing

end

"""
function to compress linear combination of mps stored in array A, which
have weights given by coef and where H is applied to the last MPS.
The compression uses the two site variational algorithm with A0 as initial guess.
A0 must be in left canonical form (normalisation not necessary)
"""

function compress_sum_var_apply_H_two_site(Ain,coef,A,H,maxsweeps,svd_tol,dims,dlim)

    na = length(Ain)
    numv = length(coef)
    FL = Array{Array{Array{Complex{Float64},2},1},1}(numv-1)
    FR = Array{Array{Array{Complex{Float64},2},1},1}(numv-1)
    FLH = Array{Array{Complex{Float64},3},1}(na)
    FRH = Array{Array{Complex{Float64},3},1}(na)
    A0 = Array{Array{Complex{Float64},3},1}(na)

    #initializing left environment

    for ll = 1:numv-1
        #A[ll][na] = coef[ll]*Av[ll][na]
        FL[ll] = Array{Array{Complex{Float64},2},1}(na)
        FL[ll][1] = eye(1)
        for jj = 1:na-2
            FL[ll][jj+1] = update_lenv(Ain[jj],A[ll][jj],FL[ll][jj])
        end
        FR[ll] = Array{Array{Complex{Float64},2},1}(na)
        FR[ll][na] = coef[ll]*eye(1);
    end

    FLH[1] = ones(1,1,1)
    for jj = 1:na-1
        FLH[jj+1] = update_le(Ain[jj],H[jj],A[numv][jj],FLH[jj])
    end
    FRH[na] = coef[numv]*ones(1,1,1)

    for sweep = 1:maxsweeps

        U,S,Vd = eye(1),eye(1),eye(1)
        #sweep left

        for jj = na:-1:2;

            av = A_prod_LR_two_site(FL[1][jj-1],A[1][jj-1],A[1][jj],FR[1][jj])
            for ll = 2:numv-1
                av += A_prod_LR_two_site(FL[ll][jj-1],A[ll][jj-1],A[ll][jj],FR[ll][jj])
            end
            av += A_prod_LOR_two_site(FLH[jj-1],A[numv][jj-1],H[jj-1],A[numv][jj],H[jj],FRH[jj])
            U,S,Vd = svdecon(reshape(av,dims[jj-1]*size(A[1][jj-1],2),
                size(A[1][jj],2)*dims[jj+1]))
            S = S[S.>0]
            stot = sum(S.^2)
            dims[jj] = min(sum(cumsum(S.^2/stot) .< 1-svd_tol)+1,dlim)
            A0[jj] = reshape(Vd[:,1:dims[jj]]',dims[jj],size(A[1][jj],2),dims[jj+1])
            for ll = 1:numv-1
                FR[ll][jj-1] = update_renv(A0[jj],A[ll][jj],FR[ll][jj])
            end
            FRH[jj-1] = update_re(A0[jj],H[jj],A[numv][jj],FRH[jj])

        end

        #A0[1] = reshape(U[:,1:Dims[2]]*diagm(S[1:Dims[2]]),Dims[1],size(A[1][1],2),Dims[2])

        #sweep right

        for jj = 1:na-1;

            av = A_prod_LR_two_site(FL[1][jj],A[1][jj],A[1][jj+1],FR[1][jj+1])
            for ll = 2:numv-1
                av += A_prod_LR_two_site(FL[ll][jj],A[ll][jj],A[ll][jj+1],FR[ll][jj+1])
            end
            av += A_prod_LOR_two_site(FLH[jj],A[numv][jj],H[jj],A[numv][jj+1],H[jj+1],FRH[jj+1])
            U,S,Vd = svdecon(reshape(av,dims[jj]*size(A[1][jj],2),
                size(A[1][jj+1],2)*dims[jj+2]))
            S = S[S.>0]
            stot = sum(S.^2)
            dims[jj+1] = min(sum(cumsum(S.^2/stot) .< 1-svd_tol)+1,dlim)
            A0[jj] = reshape(U[:,1:dims[jj+1]],dims[jj],size(A[1][jj],2),dims[jj+1])
            for ll = 1:numv-1
                FL[ll][jj+1] = update_lenv(A0[jj],A[ll][jj],FL[ll][jj])
            end
            FLH[jj+1] = update_le(A0[jj],H[jj],A[numv][jj],FLH[jj])

        end

        A0[na] = reshape(diagm(S[1:dims[na]])*(Vd[:,1:dims[na]]'),dims[na],size(A[1][na],2),dims[na+1])

    end

    return A0

end

function compress_sum_var_apply_H_two_site_left(Ain,coef,A,H,maxsweeps,svd_tol,dims,dlim)

    na = length(Ain)
    numv = length(coef)
    FL = Array{Array{Array{Complex{Float64},2},1},1}(numv-1)
    FR = Array{Array{Array{Complex{Float64},2},1},1}(numv-1)
    FLH = Array{Array{Complex{Float64},3},1}(na)
    FRH = Array{Array{Complex{Float64},3},1}(na)
    A0 = Array{Array{Complex{Float64},3},1}(na)

    #initializing left environment

    for ll = 1:numv-1
        #A[ll][na] = coef[ll]*Av[ll][na]
        FL[ll] = Array{Array{Complex{Float64},2},1}(na)
        FL[ll][1] = eye(1)
        for jj = 1:na-2
            FL[ll][jj+1] = update_lenv(Ain[jj],A[ll][jj],FL[ll][jj])
        end
        FR[ll] = Array{Array{Complex{Float64},2},1}(na)
        FR[ll][na] = coef[ll]*eye(1);
    end

    FLH[1] = ones(1,1,1)
    for jj = 1:na-2
        FLH[jj+1] = update_le(Ain[jj],H[jj],A[numv][jj],FLH[jj])
    end
    FRH[na] = coef[numv]*ones(1,1,1)

    U,S,Vd = eye(1),eye(1),eye(1)
    #sweep left

    for jj = na:-1:2;

        av = A_prod_LR_two_site(FL[1][jj-1],A[1][jj-1],A[1][jj],FR[1][jj])
        for ll = 2:numv-1
            av += A_prod_LR_two_site(FL[ll][jj-1],A[ll][jj-1],A[ll][jj],FR[ll][jj])
        end
        av += A_prod_LOR_two_site(FLH[jj-1],A[numv][jj-1],H[jj-1],A[numv][jj],H[jj],FRH[jj])
        U,S,Vd = svdecon(reshape(av,dims[jj-1]*size(A[1][jj-1],2),
            size(A[1][jj],2)*dims[jj+1]))
        S = S[S.>0]
        stot = sum(S.^2)
        dims[jj] = min(sum(cumsum(S.^2/stot) .< 1-svd_tol)+1,dlim)
        A0[jj] = reshape(Vd[:,1:dims[jj]]',dims[jj],size(A[1][jj],2),dims[jj+1])
        for ll = 1:numv-1
            FR[ll][jj-1] = update_renv(A0[jj],A[ll][jj],FR[ll][jj])
        end
        FRH[jj-1] = update_re(A0[jj],H[jj],A[numv][jj],FRH[jj])

    end

    A0[1] = reshape(U[:,1:dims[2]]*diagm(S[1:dims[2]]),dims[1],size(A[1][1],2),dims[2])

    return A0

end

function compress_sum_var_apply_H_two_site_right(Ain,coef,A,H,maxsweeps,svd_tol,dims,dlim)

    na = length(Ain)
    numv = length(coef)
    FL = Array{Array{Array{Complex{Float64},2},1},1}(numv-1)
    FR = Array{Array{Array{Complex{Float64},2},1},1}(numv-1)
    FLH = Array{Array{Complex{Float64},3},1}(na)
    FRH = Array{Array{Complex{Float64},3},1}(na)
    A0 = Array{Array{Complex{Float64},3},1}(na)

    #initializing right environment

    for ll = 1:numv-1
        FR[ll] = Array{Array{Complex{Float64},2},1}(na)
        FR[ll][na] = eye(1)
        for jj = na:-1:2
            FR[ll][jj-1] = update_renv(Ain[jj],A[ll][jj],FR[ll][jj])
        end
        FL[ll] = Array{Array{Complex{Float64},2},1}(na)
        FL[ll][1] = coef[ll]*eye(1)
    end

    FLH[1] = coef[numv]*ones(1,1,1)
    FRH[na] = ones(1,1,1)
    for jj = na:-1:2;
        FRH[jj-1] = update_re(Ain[jj],H[jj],A[numv][jj],FRH[jj])
    end

    U,S,Vd = eye(1),eye(1),eye(1)
    #sweep right

    for jj = 1:na-1;

        av = A_prod_LR_two_site(FL[1][jj],A[1][jj],A[1][jj+1],FR[1][jj+1])
        for ll = 2:numv-1
            av += A_prod_LR_two_site(FL[ll][jj],A[ll][jj],A[ll][jj+1],FR[ll][jj+1])
        end
        av += A_prod_LOR_two_site(FLH[jj],A[numv][jj],H[jj],A[numv][jj+1],H[jj+1],FRH[jj+1])
        U,S,Vd = svdecon(reshape(av,dims[jj]*size(A[1][jj],2),
            size(A[1][jj+1],2)*dims[jj+2]))
        S = S[S.>0]
        stot = sum(S.^2)
        dims[jj+1] = min(sum(cumsum(S.^2/stot) .< 1-svd_tol)+1,dlim)
        A0[jj] = reshape(U[:,1:dims[jj+1]],dims[jj],size(A[1][jj],2),dims[jj+1])
        for ll = 1:numv-1
            FL[ll][jj+1] = update_lenv(A0[jj],A[ll][jj],FL[ll][jj])
        end
        FLH[jj+1] = update_le(A0[jj],H[jj],A[numv][jj],FLH[jj])

    end

    A0[na] = reshape(diagm(S[1:dims[na]])*(Vd[:,1:dims[na]]'),dims[na],size(A[1][na],2),dims[na+1])

    return A0

end

"""
Compresses a left-canonical MPS to a maximum bond dimension D_max using a
relative tolerance tol in the trimming, giving a right-canonical MPS.
"""
function compressionSVD!(mps, dmax, tol = 1e-15)

    n = length(mps)
    lambda = Array{Array}(n - 1)
    err_trim = Array{Float64}(n)
    err_glob = Float64
    ent_ent = Array{Float64}(n)
    dims = Array{Int64}(n - 1)
    dnew = 1
    r = 1

    for i = n:-1:1

        dl, d, dr = size(mps[i])
        mps_site = reshape(mps[i], (dl*d, dr))*r
        #u, s, v = svdecon(reshape(mps_site, (dl, dnew*d)))
        u, s, v = svd(reshape(mps_site, (dl, dnew*d)))
        s = s[s .> 0]
        stot = sum(abs2.(s))
        b = dnew
        dnew = minimum([sum(cumsum(abs2.(s)/stot) .< (1 - tol)) + 1, dmax])
        err_trim[i] = sum(abs2.(s[(dnew + 1):end]))/stot
        ent_ent[i] = -2*sum(abs2.(s).*log.(s))
        mps[i] = reshape(v[:, 1:dnew]', (dnew, d, b))
        r = u[:, 1:dnew]*diagm(s[1:dnew])

    end

    err_glob = 1 - abs2(r[1])

    return err_glob, err_trim, ent_ent

end


text:
using LinearAlgebra

# Fourier Types {Parameter type, full or half flag, Differentiable Flag}
mutable struct FourierBasis{T,TF} <: AbstractFuncApprox where {T<:Real,TF<:Bool}
    C::Array{T, 2}
    range_low::Array{T, 1}
    range_width::Array{T, 1}
    function FourierBasis(ranges::Array{T, 2}, dorder::Int, iorder::Int, full::Bool=false) where {T<:Real}
        tp = eltype(ranges)
        new{tp,full}(make_Cmat(T, size(ranges)[1], dorder, iorder), ranges[:, 1], ranges[:, 2] - ranges[:, 1])
    end
    FourierBasis(C::Array{T, 2}, range_low::Vector{T}, range_width::Vector{T}, full::Bool) where {T} = new{T,full}(deepcopy(C), deepcopy(range_low), deepcopy(range_width))
end

function increment_counter!(counter::Vector{Int}, maxDigit::Int)
    for i in length(counter):-1:1
        counter[i] += 1
        if (counter[i] > maxDigit)
            counter[i] = 0
        else
            break
        end
    end
end

function make_Cmat(T::Type, num_inputs::Int, dorder::Int, iorder::Int)::Array{T, 2}
    iTerms = iorder * num_inputs
    dTerms = (dorder+1)^num_inputs
    oTerms = min(iorder, dorder) * num_inputs
    num_feats = iTerms + dTerms - oTerms
    C = zeros(T, (num_feats, num_inputs))
    counter = zeros(Int, num_inputs)
    termCount::Int = 1

    while termCount <= dTerms
        for i in 1:num_inputs
            C[termCount, i] = counter[i]
        end
        increment_counter!(counter, dorder)
        termCount += 1
    end
    for i in 1:num_inputs
        for j in (dorder+1):(iorder)
            C[termCount, i] = j
            termCount += 1
        end
    end
    C .*= π

    return C  # transpose(C)
end

function processinputs!(feats, ϕ::FourierBasis{T}, x::Array{T,1}) where {T<:Real}
    scaled = @. (x - ϕ.range_low) / ϕ.range_width
    # feats .= ϕ.C * scaled
    @. feats = 0.0
    for i in 1:length(x)
        feats .+= view(ϕ.C, :, i) .* scaled[i]
    end
    nothing
end

function call(ϕ::FourierBasis{T,false}, x::Array{T, 1})::Array{T,1} where {T<:Real}
    feats = zeros(T, get_num_outputs(ϕ))
    processinputs!(feats, ϕ, x)
    feats .= cos.(feats)
    return feats
end

function call!(out, ϕ::FourierBasis{T,false}, x::Array{T, 1}) where {T<:Real}
    # feats = zeros(T, get_num_outputs(ϕ))
    # processinputs!(feats, ϕ, x)
    # out .= cos.(feats)
    processinputs!(out, ϕ, x)
    @. out = cos(out)
end

function call(ϕ::FourierBasis{T,true}, x::Array{T, 1})::Array{T,1} where {T<:Real}
    feats = zeros(T, get_num_outputs(ϕ))
    call!(feats, ϕ, x)
    # numhalf= size(ϕ.C)[1]
    # tmp = zeros(T, numhalf)
    # processinputs!(tmp, ϕ, x)
    # feats[numhalf+1:end] .= sin.(tmp)
    # feats[1:numhalf] .= cos.(tmp)
    return feats
end

function call!(out, ϕ::FourierBasis{T,true}, x::Array{T, 1}) where {T<:Real}
    numhalf= size(ϕ.C)[1]
    # tmp = zeros(T, numhalf)
    # processinputs!(tmp, ϕ, x)
    # out[numhalf+1:end] .= sin.(tmp)
    # out[1:numhalf] .= cos.(tmp)
    processinputs!(view(out, 1:numhalf), ϕ, x)
    out[numhalf+1:end] .= sin.(view(out, 1:numhalf))
    out[1:numhalf] .= cos.(view(out, 1:numhalf))
end

function get_params(ϕ::FourierBasis)
    return ϕ.C
end

function copy_params!(params::Array{T}, ϕ::FourierBasis{T}) where {T<:Real}
    p = reshape(params, size(ϕ.C))
    copyto!(p, ϕ.C)
end

function copy_params(ϕ::FourierBasis{T}) where {T<:Real}
    params = deepcopy(ϕ.C)
    return params
end

function set_params!(ϕ::FourierBasis{T}, θ::Array{T}) where {T<:Real}
    ϕ.C .= reshape(θ,size(ϕ.C))
end

function add_to_params!(ϕ::FourierBasis{T}, θ::Array{T}) where {T<:Real}
    ϕ.C .+= reshape(θ, size(ϕ.C))
end

function get_num_params(ϕ::FourierBasis)::Int64
    return length(ϕ.C)
end

function get_num_outputs(ϕ::FourierBasis{T,false})::Int where {T}
    return size(ϕ.C)[1]
end

function get_num_outputs(ϕ::FourierBasis{T,true})::Int where {T}
    return 2. * size(ϕ.C)[1]
end

function get_output_range(ϕ::FourierBasis{T}) where {T}
    ranges = zeros(T, (get_num_outputs(ϕ), 2))
    ranges[:, 1] .= -1.
    ranges[:, 2] .= 1.
    #ranges[1, 1]  = 1.  # first term is always bias
    return ranges
end

function clone(ϕ::FourierBasis{T, TF})::FourierBasis{T, TF} where {T<:Real,TF}
    ϕ2 = FourierBasis(ϕ.C, ϕ.range_low, ϕ.range_width, TF)
    return ϕ2
end

ext: jl

max_stars_repo_name: ScottJordan/EvaluationOfRLAlgs

max_issues_repo_name: ScottJordan/EvaluationOfRLAlgs

max_forks_repo_name: ScottJordan/EvaluationOfRLAlgs


text:
# TODO: maybe more fancy file format and correctness checking should be done

abstract type Abstract_HG_format end
struct HGF_Format <: Abstract_HG_format end
struct JSON_Format <: Abstract_HG_format end

"""
    hg_save(io::IO, h::Hypergraph, format::HGF_Format)

Saves a hypergraph `h` to an output stream `io` in `hgf` format.

"""
function hg_save(io::IO, h::Hypergraph, format::HGF_Format)
    println(io, length(h.v2he), " ", length(h.he2v))
    for he in h.he2v
        skeys = sort(collect(keys(he)))
        println(io, join(["$k=$(he[k])" for k in skeys], ' '))
    end
end

"""
    hg_save(io::IO, h::Hypergraph, format::JSON_Format)

Saves a hypergraph `h` to an output stream `io` in `json` format.

If `h` has `Composite Types` either for vertex metadata or hyperedges metadata,
the user has to explicit tell the JSON3 package about it, for instance using:

`JSON3.StructType(::Type{MyType}) = JSON3.Struct()`.

See the (JSON3.jl documentation)[https://github.com/quinnj/JSON3.jl] for more details.

The `json` in output contains the following information (keys):

* `n` : number of vertices
* `k` : number of hyperedges
* `m` : a matrix representation of `h` where rows are vertices and columns are hyperedges
* `v2he` : mapping vertices to hyperedges
* `v_meta` : vertices metadata
* `he_meta` : hyperedges metadata

"""
function hg_save(io::IO, h::Hypergraph, format::JSON_Format)
    json_hg = Dict{Symbol, Any}()

    json_hg[:n] = nhv(h)
    json_hg[:k] = nhe(h)

    json_hg[:m] = JSON3.write(Matrix(h))
    json_hg[:v2he] = JSON3.write(h. v2he)

    json_hg[:v_meta] = JSON3.write(h.v_meta)
    json_hg[:he_meta] = JSON3.write(h.he_meta)

    JSON3.write(io, json_hg)
end

"""
    hg_save(
        fname::AbstractString, h::Hypergraph;
        format::Abstract_HG_format=HGF_Format()
    )

Saves a hypergraph `h` to a file `fname` in the specified `format`.
The default saving format is `hgf`.

"""
hg_save(
    fname::AbstractString, h::Hypergraph;
    format::Abstract_HG_format=HGF_Format()) =
    open(io -> hg_save(io, h, format), fname, "w")

"""
    hg_load(
        io::IO,
        format::HGF_Format();
        T::Type{U} = Bool,
        D::Type{<:AbstractDict{Int, U}} = Dict{Int,U},
        V=Nothing,
        E=Nothing
    ) where {U <: Real}

Loads a hypergraph from a stream `io` from `hgf` format.

**Arguments**

* `T` : type of weight values stored in the hypergraph's adjacency matrix
* `D` : dictionary for storing values the default is `Dict{Int, T}`
* `V` : type of values stored in the vertices of the hypergraph
* `E` : type of values stored in the edges of the hypergraph

Skips a single initial comment.

"""
function hg_load(
    io::IO,
    format::HGF_Format;
    T::Type{U} = Bool,
    D::Type{<:AbstractDict{Int, U}} = Dict{Int, T},
    V=Nothing,
    E=Nothing
) where {U <: Real}
    line = readline(io)

    if startswith(line, "\"\"\"")
      singleline = true
        while(
            !( (!singleline && endswith(line, "\"\"\"")) ||
            (singleline && endswith(line, "\"\"\"") && length(line)>5)
            ) &&
            !eof(io)
            )
                line = readline(io)
                singleline = false
        end
        if eof(io)
            throw(ArgumentError("malformed input"))
        end
       line = readline(io)
    end

    l = split(line)
    length(l) == 2 || throw(ArgumentError("expected two integers"))
    n, k = parse.(Int, l)
    h = Hypergraph{T, V, E, D}(n, k)

    for i in 1:k
        lastv = 0
        for pos in split(readline(io))
            entry = split(pos, '=')
            length(entry) == 2 || throw(ArgumentError("expected vertex=weight"))
            v = parse(Int, entry[1])
            w = parse(T, entry[2])

            if v > lastv
                lastv = v
            else
                throw(ArgumentError("vertices in hyperedge must be sorted"))
            end
            h[v, i] = w
        end
    end
    # we ignore lines beyond k+1 in the file
    h
end

"""
    hg_load(
        io::IO,
        format::JSON_Format;
        T::Type{U} = Bool,
        D::Type{<:AbstractDict{Int, U}} = Dict{Int,U},
        V = Nothing,
        E = Nothing
    ) where {U <: Real}

Loads a hypergraph from a stream `io` from `json` format.

**Arguments**

* `T` : type of weight values stored in the hypergraph's adjacency matrix
* `D` : dictionary for storing values the default is `Dict{Int, T}`
* `V` : type of values stored in the vertices of the hypergraph
* `E` : type of values stored in the edges of the hypergraph

"""
function hg_load(
        io::IO,
        format::JSON_Format;
        T::Type{U} = Bool,
        D::Type{<:AbstractDict{Int, U}} = Dict{Int, T},
        V = Nothing,
        E = Nothing
    ) where {U <: Real}
    json_hg = JSON3.read(readline(io))

    m = reshape(JSON3.read(json_hg.m, Array{Union{T, Nothing}}), json_hg.n, json_hg.k)

    v_meta = Vector{Union{Nothing,V}}(nothing, size(m, 1))
    he_meta = Vector{Union{Nothing,E}}(nothing, size(m, 2))

    if V != Nothing
        v_meta = JSON3.read(json_hg.v_meta, Array{Union{V, Nothing}})
    end

    if E != Nothing
        he_meta = JSON3.read(json_hg.he_meta, Array{Union{E, Nothing}})
    end

    h = Hypergraph{T, V, E, D}(m; v_meta=v_meta, he_meta=he_meta)

    h
end

"""
    hg_load(
        fname::AbstractString;
        format::Abstract_HG_format = HGF_Format(),
        T::Type{U} = Bool,
        D::Type{<:AbstractDict{Int, U}} = Dict{Int,U},
        V = Nothing,
        E = Nothing) where {U <: Real}
    )

Loads a hypergraph from a file `fname`.
The default saving format is `hgf`.

**Arguments**

* `T` : type of weight values stored in the hypergraph's adjacency matrix
* `D` : dictionary for storing values the default is `Dict{Int, T}`
* `V` : type of values stored in the vertices of the hypergraph
* `E` : type of values stored in the edges of the hypergraph

"""
hg_load(
    fname::AbstractString;
    format::Abstract_HG_format = HGF_Format(),
    T::Type{U} = Bool,
    D::Type{<:AbstractDict{Int, U}} = Dict{Int, T},
    V = Nothing,
    E = Nothing) where {U <: Real} =
    open(io -> hg_load(io, format; T=T, D=D, V=V, E=E), fname, "r")


text:
using BoxLeastSquares
using Distributions
using StableRNGs
using Test
using Unitful

rng = StableRNG(129325)

function make_data(N=500)
    t = rand(rng, Uniform(0, 10), N)
    y = ones(N)
    dy = rand(rng, Uniform(0.005, 0.010), N)
    period = 2.0
    t0 = 0.5
    duration = 0.16
    depth = 0.2
    mask = @. abs((t - t0 + 0.5 * period) % period - 0.5 * period) < 0.5 * duration
    y[mask] .-= depth
    y .+= dy .* randn(rng, N)
    return t, y, dy, (;period, t0, duration, depth)
end

@testset "BoxLeastSquares.jl" begin
    include("bls.jl")
    include("plotting.jl")
    include("printing.jl")
end


text:
# extensionsequences.jl

#################################################
# Abstract extension type and generic interface
#################################################

# In general, extension sequences specialize on the type of the underlying
# data. This means operations with extension sequences are typically fast.
# Since a point to an array is stored, memory is allocated every time an
# extension sequence is constructed. This may be avoided in future versions
# of Julia.

"""
Any subtype of ExtensionSequence embeds an indexable vectorlike object with finite length
and extends it into a bi-infinite sequence. The extension can be indexed with any
integer.

Indexing with an index in the allowed range of the embedded vector simply returns
the corresponding element. Indexing with indices outside this range result in
a computation. For example, for periodic extensions, the indices are mapped to
the allowed range using periodicity.

Usually the indices of the embedded vector `a` range from `1` to `length(a)`, but
that need not be the case. The extension respects the range returned by
`linearindices(a)`. One could use an OffsetArray with indices `-2:2`.

In current versions of Julia, construction of an ExtensionSequence allocates
memory, because the types store a reference to the embedded vector.
"""
abstract ExtensionSequence{T} <: Sequence{T}

# We assume that the embedded vector is in the field 'a'.
# A convention we follow is that variable 'k' refers to indices of the sequence,
# and variable 'i' to indices of the embedded vector.

"The embedded vector of the extension sequence."
subvector(s::ExtensionSequence) = s.a

"The length of the embedded vector of the extension sequence."
sublength(s::ExtensionSequence) = length(linearindices(s.a))

"The first linear index of the embedded vector."
first_subindex(s::ExtensionSequence) = first(linearindices(s.a))

"The last linear index of the embedded vector."
last_subindex(s::ExtensionSequence) = last(linearindices(s.a))

# The indexing logic can refer to the `element` function defined for all extension types
# Concrete subtypes of ExtensionSequence have to implement `extensiontype`.
getindex(s::ExtensionSequence, k) = element(extensiontype(s), subvector(s), k)

# By implementing `conj` and `reverse`, we get `transpose` and `ctranspose` for free
# Concrete types should implement a `similar` function that returns a sequence of
# the same type but with a different embedded vector.
for op in (:conj,)
    @eval $op(s::ExtensionSequence) = similar(s, $op(subvector(s)))
end

reverse(s::ExtensionSequence) = _reverse(s, subvector(s))
_reverse(s, a) = similar(s, reverse_indices(a))

#######################
# Periodic extensions
#######################

"""
A PeriodicSequence extends a vector `a` periodically to a bi-infinite sequence.

The indices of `a` map to `a`. Indices outside this range are mapped
by periodization modulo the length of `a`.
"""
immutable PeriodicSequence{A,T} <: ExtensionSequence{T}
    a       :: A
end

PeriodicSequence(a::AbstractVector) = PeriodicSequence{typeof(a),eltype(a)}(a)

similar(p::PeriodicSequence, a) = PeriodicSequence(a)

extensiontype(s::PeriodicSequence) = PeriodicExtension()

#######################
# CompactSequence
#######################

"""
CompactSequence extends a vector `a` with zeros to a bi-infinite sequence.

The indices of `a` map to `a`. Indices outside this range correspond
to zero values.
"""
immutable CompactSequence{A,T} <: ExtensionSequence{T}
    a :: A
end

typealias ZeroPaddingSequence CompactSequence

CompactSequence{A}(a::A) = CompactSequence{A,eltype(A)}(a)

similar(s::CompactSequence, a) = CompactSequence(a)

extensiontype(s::CompactSequence) = ZeroPadding()

iscompact(::CompactSequence) = true

nzrange(s::CompactSequence) = linearindices(subvector(s))

#######################
# Constant padding
#######################

"""
A ConstantPaddingSequence extends a vector `a` with a given `constant` to a bi-infinite sequence.

The indices of `a` map to `a`. Indices outside this range correspond to the constant
value.
"""
immutable ConstantPaddingSequence{A,T} <: ExtensionSequence{T}
    a           ::  A
    constant    ::  T
end

ConstantPaddingSequence(a::AbstractVector, constant) = ConstantPaddingSequence{typeof(a),eltype(a)}(a, constant)

constant(s::ConstantPaddingSequence) = s.constant

similar(s::ConstantPaddingSequence, a) = ConstantPaddingSequence(a, constant(s))

extensiontype(s::ConstantPaddingSequence) = ConstantPadding(s.constant)

#######################
# Symmetric extensions
#######################

"""
A SymmetricSequence extends a vector 'a' symmetrically to a bi-infinite sequence.

The indices of `a` map to `a`. Indices outside this range are mapped by symmetrizing.

The symmetry around each of the endpoints can be whole-point (the endpoint is not repeated)
or half-point (the endpoint is repeated). The symmetry can also be even (symmetric)
or odd (anti-symmetric).

Definition:

immutable SymmetricSequence{SYM_LEFT,SYM_RIGHT,PAR_LEFT,PAR_RIGHT,A,T} <: ExtensionSequence{T}

Type parameters indicate the type of symmetry and its parity:
- SYM_LEFT:  either :wp (whole point) or :hp (half point) near left endpoint
- SYM_RIGHT: either :wp or :hp for right endpoint
- PAR_LEFT: either :odd or :even symmetry near left endpoint
- PAR_RIGHT: also :odd or :even
- A: the type of the embedded vector
- T: the eltype of the sequence
"""
immutable SymmetricSequence{SYM_LEFT,SYM_RIGHT,PAR_LEFT,PAR_RIGHT,A,T} <: ExtensionSequence{T}
    a :: A
end

SymmetricSequence(a::AbstractVector) = symmetric_extension_wholepoint_even(a)

similar{SL,SR,PL,PR}(s::SymmetricSequence{SL,SR,PL,PR}, a) =
    SymmetricSequence{SL,SR,PL,PR,typeof(a),eltype(a)}(a)

# Provide four of the sixteen combinations for convenience. Construct the others by explicitly
# specifying the type parameters symbols.
symmetric_extension_wholepoint_even(a) =
    SymmetricSequence{:wp,:wp,:even,:even,typeof(a),eltype(a)}(a)

symmetric_extension_halfpoint_even(a) =
    SymmetricSequence{:hp,:hp,:even,:even,typeof(a),eltype(a)}(a)

symmetric_extension_wholepoint_odd(a) =
    SymmetricSequence{:wp,:wp,:odd,:odd,typeof(a),eltype(a)}(a)

symmetric_extension_halfpoint_odd(a) =
    SymmetricSequence{:hp,:hp,:odd,:odd,typeof(a),eltype(a)}(a)

# We have to be careful when reversing a symmetric extension: the endpoints change places
function reverse{SL,SR,PL,PR}(s::SymmetricSequence{SL,SR,PL,PR})
    b = reverse_indices(subvector(s))
    SymmetricSequence{SR,SL,PR,PL,typeof(b),eltype(b)}(b)
end

for op in (:left_parity, :right_parity, :left_symmetry, :right_symmetry)
    @eval $op(s::SymmetricSequence) = $op(extensiontype(s))
end

extensiontype{SL,SR,PL,PR}(s::SymmetricSequence{SL,SR,PL,PR}) =
    SymmetricExtension{SL,SR,PL,PR}()


text:
using OneHot
using Base: tail, front

X = OneHotArray(rand(1:4, 5), 4)

for i = 1:size(X,2)
    @test OneHot.decode(X[:,i]) == X.c[i]
end

A = randn(3,4)
@test A[:, X.c] == @inferred A * X
@test A * X == A * Array(X)

B = randn(3,3)
@test_throws DimensionMismatch B * X

@test X[:,1:3] isa OneHotArray{2}
@inferred X[:,1:3]

for I in eachindex(X)
    i = Tuple(I)
    @test X[I] == (first(i) == X.c[tail(i)...])
end

@testset "convert to/from BitArray" begin
    X = OneHotArray(rand(1:4, 10, 5), 4)
    @test Array(X) == BitArray(X) == X
    @test OneHotArray(BitArray(X)) == X
    @test OneHotArray(X) == X
end


text:
using VDPTag2
using Plots
using POMDPToolbox
using Reel
using ProgressMeter
using ParticleFilters
using StaticArrays

frames = Frames(MIME("image/png"), fps=2)

# pomdp = VDPTagPOMDP()
pomdp = VDPTagPOMDP(mdp=VDPTagMDP(barriers=CardinalBarriers(0.2, 2.8)))
policy = ManageUncertainty(pomdp, 0.01)
# policy = ToNextML(mdp(pomdp))

rng = MersenneTwister(5)

hr = HistoryRecorder(max_steps=30, rng=rng)
filter = SIRParticleFilter(pomdp, 200, rng=rng)
hist = simulate(hr, pomdp, policy, filter)

gr()
@showprogress "Creating gif..." for i in 1:n_steps(hist)
    push!(frames, plot(pomdp, view(hist, 1:i)))
end

filename = string(tempname(), "_vdprun.gif")
write(filename, frames)
println(filename)
run(`setsid xdg-open $filename`)


text:
using Makie, Dates, Observables, Printf, Images

mss = [Millisecond(25), Millisecond(125)]
t = Millisecond(Second(5))
fs = Int.(t./mss)
fs[end] += 1
Δt = vcat((fill(ms, f) for (ms, f) in zip(mss, fs))...)
Δt[1] = Millisecond(0)
ts = cumsum(Δt)
cd(tempdir())
dar = 2
w = length(ts)
h = dar*w
frame = Observable(1)
point = lift(frame) do fn
    [(fn, dar*fn)]
end
msg = lift(frame, point) do fn, p
    ti = @sprintf "%03i" fn
    t = @sprintf "%-012s" Time(0) + ts[fn]
    x, y = p[]
    """
Frame #: $ti
Time: $t
x: $x
y: $y"""
end
limits = FRect3D(Vec3f0(1, 2, 0), Vec3f0(w - 1, 2w - 2, 0))
scene = Scene(resolution = (w, h), limits = limits,  scale_plot = false, show_axis = false)#, padding=(0,0), backgroundcolor=:red)
# poly!(scene, limits, color = :gray)
text!(scene, msg, position = (2, w), align = (:left,  :bottom))#, color = :white)
scatter!(scene, point, markersize = 10, marker = :+)#, color = :white)
for i in 1:w
    frame[] = i
    AbstractPlotting.update_limits!(scene)
    AbstractPlotting.update!(scene)
    fname = @sprintf "%03i.png" i
    save(fname, scene)
end
run(`mogrify -resize $(h)x$h\! "*".png`)
run(`ffmpeg -y -i %03d.png -pix_fmt yuv420p cfr.mp4`)
f1, f2 = fs
fps1, fps2 = Millisecond(Second(1)) .÷ mss
run(`mp4fpsmod -r $f1:$fps1  -r $f2:$fps2 -o vfr.mp4 cfr.mp4`)

run(`ffmpeg -y -i vfr.mp4 -c copy -bsf:v h264_metadata=sample_aspect_ratio=0.5 testvideo.mp4`)

#=using VideoIO
f = VideoIO.openvideo("vfr2.mp4")
pixelaspectratio = VideoIO.aspect_ratio(f)

f = VideoIO.openvideo("cfr.mp4")
ts2 = Millisecond[]
while !eof(f)
    read(f)
    # l = gettime(f)
    # push!(ts2, Millisecond(round(Int, 1000l)))
    push!(ts2, Millisecond(1))
end
ts2=#

# here we create `n` frames in the form of individual images
using FileIO, Colors, Printf
n = 6
levels = UInt8.(range(0, stop = 255, length = n))
path = mktempdir()
img = Matrix{UInt8}(undef, 10, 10)
for (i, l) in enumerate(levels)
    name = @sprintf "%03i.jpg" i
    fill!(img, l)
    save(joinpath(path, name), img)
end

# here we convert these `n` images into a video
x = joinpath(path, "%03d.jpg")
y = joinpath(path, "out.mp4")
fps = 25
run(`ffmpeg -i $x -r $fps -pix_fmt gray $y`)

duration = n/fps
times = range(0, step = 1/fps, length = n + 1)
@assert duration == last(times)

# now let's see how many frames VideoIO thinks are in this video
# and what images "made it"
using VideoIO, Statistics
f = VideoIO.openvideo(y)
fn = 0
pix = Float64[]
times = Float64[]
while !eof(f)
    global fn
    img = read(f)
    push!(pix, mean(Gray.(img)))
    push!(times, gettime(f))
    fn += 1
end
@show fn # I get `n - 2` here
@show pix # I get [0.0, 0.2, 0.4, 0.6] So the last 2 are missing
@show sum(times)

duration = VideoIO.get_duration(f.avin.io)
run(`ffmpeg -i $y -f null -`)


text:
abstract type Backdrop

end

@valueequals @fieldproxy mutable struct Parallax <: Backdrop
    data::Dict{String, Any}

    Parallax(data::Dict{String, Any}) = new(data)
    Parallax(; kwargs...) = new(Dict(string(a) => b for (a, b) in kwargs))
end

@valueequals @fieldproxy mutable struct Effect{T} <: Backdrop
    name::String
    data::Dict{String, Any}

    Effect{T}(name::String, data::Dict{String, Any}) where T = new{T}(name, data)
    Effect{T}(name::String; kwargs...) where T = new{T}(name, Dict(string(a) => b for (a, b) in kwargs))
end

Effect(name::String, data::Dict{String, Any}) = Effect{Symbol(name)}(name, data)
Effect(name::String; kwargs...) = Effect{Symbol(name)}(name, Dict(string(a) => b for (a, b) in kwargs))

@fieldproxy mutable struct Apply <: Backdrop
    data::Dict{String, Any}
    parallax::Array{Union{Effect, Parallax}, 1}

    Apply(data::Dict{String, Any}, parallax::Array{Union{Effect, Parallax}, 1}) = new(data, parallax)
    Apply(parallax::Array{Union{Effect, Parallax}, 1}; kwargs...) = new(Dict(string(a) => b for (a, b) in kwargs), parallax)
end

@mapdef Effect "snowFg" SnowFg(only::String="*", exclude::String="")
@mapdef Effect "snowBg" SnowBg(only::String="*", exclude::String="")
@mapdef Effect "heatwave" Heatwave(only::String="*", exclude::String="")
@mapdef Effect "dreamstars" DreamStars(only::String="*", exclude::String="")
@mapdef Effect "stars" Stars(only::String="*", exclude::String="")
@mapdef Effect "mirrorfg" MirrorFg(only::String="*", exclude::String="")
@mapdef Effect "reflectionfg" ReflectionFg(only::String="*", exclude::String="")
@mapdef Effect "corestarsfg" CoreStarFg(only::String="*", exclude::String="")
@mapdef Effect "godrays" GodRays(only::String="*", exclude::String="")
@mapdef Effect "bossStarField" BossStarField(only::String="*", exclude::String="")
@mapdef Effect "windsnow" Wind(only::String="*", exclude::String="", color::String="")
@mapdef Effect "planets" Planets(only::String="*", exclude::String="", count::Number=32, size::String="small", color::String="", scrollx::Number=1.0, scrolly::Number=1.0)
@mapdef Effect "starfield" Starfield(only::String="*", exclude::String="", color::String="", scrollx::Number=1.0, scrolly::Number=1.0, speed::Number=1.0)
@mapdef Effect "petals" Petals(only::String="*", exclude::String="")
@mapdef Effect "northernlights" NorthernLights(only::String="*", exclude::String="")
@mapdef Effect "tentacles" TentacleEffect(only::String="*", exclude::String="", side::String="Right", color::String="", offset::Number=0)
@mapdef Effect "stardust" StarDust(only::String="*", exclude::String="")
@mapdef Effect "rain" Rain(only::String="*", exclude::String="", color::String="161933") # Color added by Everest
@mapdef Effect "blackhole" BlackHole(only::String="*", exclude::String="")

mutable struct Style
    foregrounds::Array{Backdrop, 1}
    backgrounds::Array{Backdrop, 1}

    Style(fg::Array{Backdrop, 1}, bg::Array{Backdrop, 1}) = new(fg, bg)
    Style() = new(Backdrop[], Backdrop[])
end

function expandApply(style::Backdrop)
    if isa(style, Maple.Apply)
        res = Backdrop[]

        for p in style.parallax
            if isa(p, Maple.Parallax)
                push!(res, Maple.Parallax(merge(style.data, p.data)))

            elseif isa(p, Maple.Effect)
                push!(res, Maple.Effect(p.name, merge(style.data, p.data)))
            end
        end

        return res

    else
        return Backdrop[style]
    end
end

# Expand Apply types into parallax/effect
# Easier to use for both coding and user perspective
function expandStylegroundApplies(backdrops::Array{Backdrop, 1})
    res = Backdrop[]

    for backdrop in backdrops
        append!(res, expandApply(backdrop))
    end

    return res
end

function expandStylegroundApplies!(style::Style)
    style.foregrounds = expandStylegroundApplies(style.foregrounds)
    style.backgrounds = expandStylegroundApplies(style.backgrounds)
end

function Base.Dict(s::Style)
    return Dict{String, Any}[
        Dict{String, Any}(
            "__name" => "Foregrounds",
            "__children" => Dict.(s.foregrounds)
        ),

        Dict{String, Any}(
            "__name" => "Backgrounds",
            "__children" => Dict.(s.backgrounds)
        )
    ]
end

function Base.Dict(p::Parallax)
    res = copy(p.data)

    res["__name"] = "parallax"

    return res
end

function Base.Dict(a::Apply)
    res = copy(a.data)

    res["__name"] = "apply"
    res["__children"] = Dict.(a.parallax)

    return res
end

function Base.Dict(e::Effect)
    res = copy(e.data)

    res["__name"] = e.name

    return res
end

max_stars_repo_name: CelestialCartographers/Maple

max_issues_repo_name: CelestialCartographers/Maple

max_forks_repo_name: CelestialCartographers/Maple


text:
export had

# Intervention Context
Cassette.@context HadCtx

struct Intervention{X, V}
  x::X
  v::V
end

"`x` had `i` been the case"
function had(x, i::Intervention)
  let ctx = HadCtx(metadata = i)
    ω -> Cassette.overdub(ctx, x, ω)
  end
end

function Cassette.overdub(ctx::HadCtx{Intervention{X, V}}, x::X, ω::Ω) where {X, V}
  @show x
  @show ω
  @show ctx.metadata.v
  if x == ctx.metadata.x
    @show :true
    ctx.metadata.v
  else
    @show :false
    Cassette.recurse(ctx, x, ω)
  end
end

# Sugar

"Examples Syntax for `x | had(3 => x)`"
had(xtov::Pair) = Intervention(xtov.first, xtov.second)
Base.:|(x, i::Intervention) = had(x, i)


text:
mutable struct CachedCovariate{T<:AbstractFloat} <: AbstractCovariate{T}
    cache::Dict{Tuple{Integer, Integer}, AbstractVector{T}}
    basecovariate::AbstractCovariate{T}
    lockobj::Threads.TatasLock
end

Base.length(var::CachedCovariate) = length(var.basecovariate)

getname(var::CachedCovariate) = getname(var.basecovariate)

function CachedCovariate(basecovariate::AbstractCovariate{T}) where {T<:AbstractFloat}
    CachedCovariate{T}(Dict{Tuple{Integer, Integer}, AbstractVector{T}}(), basecovariate, Threads.TatasLock())
end

function cache(basecovariate::AbstractCovariate{T}) where {T<:AbstractFloat}
    CachedCovariate{T}(Dict{Tuple{Integer, Integer}, AbstractVector{T}}(), basecovariate, Threads.TatasLock())
end

function slice(covariate::CachedCovariate{T}, fromobs::Integer, toobs::Integer, slicelength::Integer) where {T<:AbstractFloat}
    basecov = covariate.basecovariate
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    if isa(basecov, CachedCovariate) || isa(basecov, Covariate)
        slice(basecov, fromobs, toobs, slicelength)
    else
        lockobj = covariate.lockobj
        lock(lockobj)
        try
            if !((fromobs, toobs) in keys(covariate.cache))
                v, _ = tryread(slice(basecov, fromobs, toobs, toobs - fromobs + 1))
                covariate.cache[(fromobs, toobs)] = get(v)
            end
        finally
            unlock(lockobj)
        end
        slice(covariate.cache[(fromobs, toobs)], 1, toobs - fromobs + 1, slicelength)
    end
end

function Base.map(covariate::CachedCovariate, dataframe::AbstractDataFrame)
    CachedCovariate(map(covariate.basecovariate, dataframe))
end

max_stars_repo_name: vikashdwd/JuML

max_issues_repo_name: vikashdwd/JuML

max_forks_repo_name: vikashdwd/JuML


text:
# Invertible network obtained from stringing invertible networks together
# Author: Gabrio Rizzuti, grizzuti3@gatech.edu
# Date: September 2020

export ComposedInvertibleNetwork, Composition
import Base.length, Base.∘

struct ComposedInvertibleNetwork <: InvertibleNetwork
    layers::Array{T, 1} where {T <: Union{NeuralNetLayer, InvertibleNetwork}}
    logdet_array::Array{Bool, 1}
    logdet::Bool
    npars::Array{Int64, 1}
end

@Flux.functor ComposedInvertibleNetwork

## Constructors

function Composition(layer...)

    # Initializing output
    depth = length(layer)
    net_array = Array{Union{NeuralNetLayer, InvertibleNetwork}, 1}(undef, depth)
    logdet_array = Array{Bool, 1}(undef, depth)
    logdet = false
    npars = Array{Int64, 1}(undef, depth)

    # Loop over layers
    for i = 1:depth

        # Selecting layer (last-to-first)
        layer_i = layer[depth-i+1]

        # Setting layer/logdets/# of parameters
        net_array[i] = layer_i
        if hasproperty(layer_i, :logdet) && layer_i.logdet
            logdet_array[i] = true
            logdet = true
        else
            logdet_array[i] = false
        end
        npars[i] = length(get_params(layer_i))

    end

    return ComposedInvertibleNetwork(net_array, logdet_array, logdet, npars)

end

## Composition utilities

function ∘(net1::ComposedInvertibleNetwork, net2::ComposedInvertibleNetwork)
    return Composition(cat(net1.layers[end:-1:1], net2.layers[end:-1:1]; dims=1)...)
end

function ∘(net1::Union{NeuralNetLayer, InvertibleNetwork}, net2::Union{NeuralNetLayer, InvertibleNetwork})
    return Composition(cat(net1, net2; dims=1)...)
end

function ∘(net1::Union{NeuralNetLayer, InvertibleNetwork}, net2::ComposedInvertibleNetwork)
    return Composition(cat(net1, net2.layers[end:-1:1]; dims=1)...)
end

function ∘(net1::ComposedInvertibleNetwork, net2::Union{NeuralNetLayer, InvertibleNetwork})
    return Composition(cat(net1.layers[end:-1:1], net2; dims=1)...)
end

function length(N::ComposedInvertibleNetwork)
    return length(N.layers)
end

## Forward/inverse/backward

function forward(X, N::ComposedInvertibleNetwork)
    N.logdet && (logdet = 0f0)
    for i = 1:length(N)
        if N.logdet_array[i]
            X, logdet_ = N.layers[i].forward(X)
            logdet += logdet_
        else
            X = N.layers[i].forward(X)
        end
    end
    N.logdet ? (return X, logdet) : (return X)
end

function inverse(Y, N::ComposedInvertibleNetwork)
    for i = length(N):-1:1
        Y = N.layers[i].inverse(Y)
    end
    return Y
end

function backward(ΔY, Y, N::ComposedInvertibleNetwork; set_grad::Bool = true)
    if ~set_grad
        Δθ = Array{Parameter, 1}(undef, 0)
        N.logdet && (∇logdet = Array{Parameter, 1}(undef, 0))
    end
    for i = length(N):-1:1
        if set_grad
            ΔY, Y = N.layers[i].backward(ΔY, Y)
        else
            if N.logdet_array[i]
                ΔY, Δθi, Y, ∇logdet_i = N.layers[i].backward(ΔY, Y; set_grad=set_grad)
                ∇logdet = cat(∇logdet_i, ∇logdet; dims=1)
            else
                ΔY, Δθi, Y = N.layers[i].backward(ΔY, Y; set_grad=set_grad)
            end
            Δθ = cat(Δθi, Δθ; dims=1)
        end
    end
    if set_grad
        return ΔY, Y
    else
        N.logdet ? (return ΔY, Δθ, Y, ∇logdet) : (return ΔY, Δθ, Y)
    end
end

## Jacobian-related utilities

function jacobian(ΔX, Δθ::Array{Parameter, 1}, X, N::ComposedInvertibleNetwork)
    N.logdet && (l = 0f0; GNΔθ = Array{Parameter, 1}(undef, 0))
    idx_pars = 0
    for i = 1:length(N)
        npars_i = N.npars[i]
        Δθi = Δθ[idx_pars+1:idx_pars+npars_i]
        if N.logdet_array[i]
            ΔX, X, li, GNΔθi = N.layers[i].jacobian(ΔX, Δθi, X)
            l += li
            GNΔθ = cat(GNΔθ, GNΔθi; dims=1)
        else
            ΔX, X = N.layers[i].jacobian(ΔX, Δθi, X)
        end
        idx_pars += npars_i
    end
    N.logdet ? (return ΔX, X, l, GNΔθ) : (return ΔX, X)
end

function adjointJacobian(ΔY, Y, N::ComposedInvertibleNetwork)
    return backward(ΔY, Y, N; set_grad = false)
end

## Other utilities

function clear_grad!(N::ComposedInvertibleNetwork)
    for i = 1:length(N)
        clear_grad!(N.layers[i])
    end
end

function get_params(N::ComposedInvertibleNetwork)
    p = Array{Parameter, 1}(undef, 0)
    for i = 1:length(N)
        p = cat(p, get_params(N.layers[i]); dims=1)
    end
    return p
end


text:
julia> split(s, r"==|!=|=")
 5-element Array{SubString{String},1}:
  "a"
  ""
  "b"
  ""
  "c"

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Smoothing
using Test

@test Smoothing.binomial([1.0,2.0,3.0,4.0,5.0], 1) == [1.0,2.0,3.0,4.0,5.0]
@test Smoothing.binomial([1.0,2.0,3.0,4.0,5.0], 3) == [1.0,2.0,3.0,4.0,5.0]
println("Passed Smoothing of counting numbers tests (1 and 3 pass)")
@test Smoothing.binomial([0.0,2.0,4.0,6.0,4.0, 2.0, 0.0], 1) == [0.0, 2.0, 4.0, 5.0, 4.0, 2.0, 0.0]
println("Passed Smoothing of manually created triangle pulse")
println("3 Tests passed!")


text:
include("AlgAssRelOrd/Order.jl")
include("AlgAssRelOrd/Elem.jl")
include("AlgAssRelOrd/Ideal.jl")