text:
"""
```
gensysct(Γ0, Γ1, c, Ψ, Π; ϵ, div)
gensysct(F::LinearAlgebra.GeneralizedSchur, c, Ψ, Π; ϵ, div)
```

Generate state-space solution to canonical-form DSGE model.

System given as
```
Γ0*Dy(t) = Γ1*y(t) + c + Ψ*z(t) + Π*η(t),
```
with z an exogenous variable process and η being endogenously white noise expectational errors.

Returned system is
```
Dy(t) = G1*y(t) + C + impact * z(t)
```

Returned values are
```
G1, C, impact, qt', a, b, z, eu
```

This code is based on a [routine](http://sims.princeton.edu/yftp/gensys/) originally copyright Chris Sims.

Here for completeness, but should not be used for most HACT/HANK models.

Also returned is the qz decomposition, qt'az' = Γ0, qt'bz' = Γ1, with a and b
upper triangular and the system ordered so that all zeros on the diagonal of b are in
the lower right corner, all cases where the real part of bii/aii is greater than or
equal to div appear in the next block above the zeros, and the remaining bii/aii's
all have bii/aii<div .  These elements can be used to construct the full backward and
forward solution.  See the paper \"Solving Linear Rational Expectations Models\",
http://eco-072399b.princeton.edu/yftp/gensys .  Note that if one simply wants the backwar
and forward projection of y on eΨlon, ignoring existence and uniqueness questions, the
projection can be computed by Fourier methods.

If `div` is omitted from argument list, a `div`>1 is calculated.

### Return codes

* `eu[1]==1` for existence
* `eu[2]==1` for uniqueness
* `eu[1]==-1` for existence for white noise η
* `eu==[-2,-2]` for coincident zeros.
"""
function gensysct(Γ0::Matrix{Float64}, Γ1::Matrix{Float64}, c::Array{Float64},
                  Ψ::Matrix{Float64}, Π::Matrix{Float64};
                  check_existence::Bool = true, check_uniqueness::Bool = true,
                  ϵ::Float64 = sqrt(eps()) * 10,
                  div::Float64 = -1.0, complex_decomposition::Bool = false)

    F = complex_decomposition ? schurfact(complex(Γ0), complex(Γ1)) : schurfact(real(Γ0), real(Γ1))
    gensysct(F, c, Ψ, Π; check_existence = check_existence, check_uniqueness = check_uniqueness,
             ϵ = ϵ, div = div)
end

# In place version. Changes Γ0, Γ1
function gensysct!(Γ0::Matrix{Float64}, Γ1::Matrix{Float64}, c::Array{Float64},
                   Ψ::Matrix{Float64}, Π::Matrix{Float64};
                   check_existence::Bool = true, check_uniqueness::Bool = true,
                   ϵ::Float64 = sqrt(eps()) * 10,
                   div::Float64 = -1.0, complex_decomposition::Bool = false)

    F = complex_decomposition ? schur!(complex(Γ0), complex(Γ1)) : schur!(real(Γ0), real(Γ1))
    gensysct(F, c, Ψ, Π; check_existence = check_existence, check_uniqueness = check_uniqueness,
             ϵ = ϵ, div = div)
end

# Uses basic schur decomposition
function gensysct(Γ1::Matrix{Float64}, c::Array{Float64}, Ψ::Matrix{Float64}, Π::Matrix{Float64};
                  check_existence::Bool = true, check_uniqueness::Bool = true,
                  ϵ::Float64 = sqrt(eps()) * 10,
                  div::Float64 = -1.0, complex_decomposition::Bool = false)
    # Assumes that Γ0 is the identity
    F = complex_decomposition ? schurfact(complex(Γ1)) : schurfact(real(Γ1))
    gensysct(F, c, Ψ, Π; check_existence = check_existence, check_uniqueness = check_uniqueness,
             ϵ = ϵ, div = div)
end

# In place version of above gensysct
function gensysct!(Γ1::Matrix{Float64}, c::Array{Float64}, Ψ::Matrix{Float64}, Π::Matrix{Float64};
                  check_existence::Bool = true, check_uniqueness::Bool = true, ϵ::Float64 = sqrt(eps()) * 10,
                  div::Float64 = -1.0, complex_decomposition::Bool = false)
    # Assumes that Γ0 is the identity
    F = complex_decomposition ? schur!(complex(Γ1)) : schur!(real(Γ1))
    gensysct(F, c, Ψ, Π; check_existence = check_existence, check_uniqueness = check_uniqueness,
             ϵ = ϵ, div = div)
end

function gensysct(F::LinearAlgebra.GeneralizedSchur, c::Array{Float64}, Ψ::Matrix{Float64}, Π::Matrix{Float64};
                  check_existence::Bool = true, check_uniqueness::Bool = true, ϵ::Float64 = sqrt(eps()) * 10,
                  div::Float64 = -1.0)
    div < 0.0 && (div = new_divct(F; ϵ = ϵ))
    eu = [0, 0]
    a, b = F.S, F.T
    n = size(a, 1)

    for i in 1:n
        if (abs(a[i, i]) < ϵ) && (abs(b[i, i]) < ϵ)
            info("Coincident zeros.  Indeterminacy and/or nonexistence.")
            eu = [-2, -2]
            G1 = Array{Float64, 2}() ;  C = Array{Float64, 1}() ; impact = Array{Float64, 2}()
            a, b, qt, z = FS.S, FS.T, FS.Q, FS.Z
            return G1, C, impact, qt', a, b, z, eu
        end
    end
    movelast = fill(false, n)
    for i in 1:n
        movelast[i] = (real(b[i, i] / a[i, i]) > div) || (abs(a[i, i]) < ϵ)
    end
    nunstab = sum(movelast)
    FS = ordschur!(F, .!movelast)
    a, b, qt, z = FS.S, FS.T, FS.Q, FS.Z

    qt1 = qt[:, 1:(n - nunstab)]
    qt2 = qt[:, (n - nunstab + 1):n]
    a2 = a[(n - nunstab + 1):n, (n - nunstab + 1):n]
    b2 = b[(n - nunstab + 1):n, (n - nunstab + 1):n]
    etawt = qt2' * Π # Ac_mul_B(qt2, Π)
    ~, ueta, deta, veta = decomposition_svdct!(etawt, ϵ = ϵ)
    zwt = qt2' * Ψ # Ac_mul_B(qt2, Ψ)
    bigev, uz, dz, vz = decomposition_svdct!(zwt, ϵ = ϵ)

    # check existence
    if isempty(bigev) || !check_existence
        exist = true
        existx = true
    else
        exist = norm(uz- (ueta * adjoint(ueta)) * uz, 2) < ϵ * n

        zwtx0 = b2 \ zwt
        zwtx = zwtx0
        M = b2 \ a2
        M = scale!(M, 1 / norm(M))
        for i in 2:nunstab
            zwtx = hcat(M * zwtx, zwtx0)
        end
        zwtx = b2 * zwtx
        bigev, ux, dx, vx = decomposition_svdct!(zwtx, ϵ = ϵ)
        existx = norm(ux - ueta * ueta' * ux, 2) < ϵ * n
    end

    etawt1 = qt1' * Π # Ac_mul_B(qt1, Π)
    bigev, ueta1, deta1, veta1 = decomposition_svdct!(etawt1, ϵ = ϵ)

    if nunstab == 0
       eu[1] = 1
    else
        if exist
            eu[1] = -1
        end
    end

    # check uniqueness
    if isempty(veta1) || !check_uniqueness
        eu[2] = 1
    else
        eu[2] = Int64(norm(veta1- (veta * adjoint(veta)) * veta1, 2) < ϵ * n)
    end

    tmat = hcat(eye(n - nunstab), -ueta1 * deta1 * conj(veta1)' * veta * (deta \ ueta'))
    G0 =  vcat(tmat * a, hcat(zeros(nunstab, n - nunstab), eye(nunstab)))
    G1 =  vcat(tmat * b, zeros(nunstab, n))

    G0I = inv(G0)
    G1 = G0I * G1
    usix = (n - nunstab + 1):n
    Busix = b[usix,usix]
    Ausix = a[usix,usix]
    C = G0I * vcat(tmat * qt' * c, (Ausix - Busix) \ (qt2' * c))
    impact = G0I * vcat(tmat * qt' * Ψ, zeros(nunstab, size(Ψ, 2)))
    # C = G0I * vcat(tmat * Ac_mul_B(qt, c), (Ausix - Busix) \ Ac_mul_B(qt2, c))
    # impact = G0I * vcat(tmat * Ac_mul_B(qt, Ψ), zeros(nunstab, size(Ψ, 2)))

    #G1 = z * A_mul_Bc(G1, z)
    G1 = z * G1 * z'
    G1 = real(G1)
    C = real(z * C) .* ones(size(z, 1), 1)
    impact = real(z * impact)

    return G1::Matrix{Float64}, C::Matrix{Float64}, impact::Matrix{Float64}, qt'::Matrix, a::Matrix, b::Matrix, z::Matrix, eu::Vector{Int64}
end

# Port of SeHyoun's schur_solver
## Inputs
# div: number of unstable roots known a priori
## Outputs
# eu[1] = 1 solution exists
#         0 no stable solution
#         -1 solution exists for white noise η
# eu[2] = 1 solution is unique
#         0 solution is not unique
function gensysct(F::LinearAlgebra.Schur, c::Array{Float64}, Ψ::Matrix{Float64}, Π::Matrix{Float64};
                  check_existence::Bool = true, check_uniqueness::Bool = true, ϵ::Float64 = sqrt(eps()) * 10,
                  div::Float64 = -1.0)
    eu = [0, 0]
    U, T = F.Z, F.T        # U is unitary matrix, T is Schur matrix
    n = size(U, 1)
    g_eigs = real(F.values) # Get real eigenvalues
    stable_eigs = g_eigs .<= 0 # stable eigenvalues
    nunstab = n - sum(stable_eigs)  # Unstable eigenvalues
    if div > -1
        sort_order = sortperm(g_eigs, rev = true) # g_eigs[sort_order] -> descending sorted g_eigs
        locs = ones(n, 1)
        locs[sort_order[1:div]] = 0 # Set to zero the eigenvalues which are unstable
        FS = ordschur!(F, locs)     # Change F in place, FS creates reference to same object, so no copying
        if nunstab > div
            @warn "<gensysct>: There are more than div number of positive eigenvalues with smallest values:\n"
            println(g_eigs[sort_order][div + 1:nunstab])
        elseif nunstab < div
            @warn "<gensysct>: There are less than div number of positive eigenvalues:\n"
            println(g_eigs[sort_order][nunstab + 1:div])
        end
        nunstab = div
    else
        FS = ordschur!(F, stable_eigs)
        U, T = FS.Z, FS.T
    end

    U1 = U[:, 1:(n - nunstab)]'
    U2 = U[:, (n - nunstab + 1):n]'
    etawt = U2 * Π # Ac_mul_B(U2, Π)
    ~, ueta, deta, veta = decomposition_svdct!(etawt, ϵ = ϵ) # Wrapper function here, see end of script, this checks for bigev and returns ueta, etc., after selecting only the relevant columns

    # check existence
    if check_existence
        zwt = U2 * Ψ   # Ac_mul_B(U2, Ψ)
        bigev, uz, dz, vz = decomposition_svdct!(zwt, ϵ = ϵ)
        if isempty(bigev)
            eu[1] = 1
        else
            eu[1] = norm(uz- (ueta * adjoint(ueta)) * uz, 2) < ϵ * n
        end

        if (eu[1] == 0 && (div == -1))
            @warn "<gensysct>: Solution does not exist"
        end
        impact = real(-Π * veta * (deta \ ueta') * uz * dz * vz' + Ψ)
    else
        eu[1] = 1
        impact = real(-Π * veta * (deta \ ueta') * U2 + Ψ)
    end

    # Check uniqueness
    if check_uniqueness
        etawt1 = U1 * Π # Ac_mul_B(U1, Π)
        ~, ~, deta1, veta1 = decomposition_svdct!(etawt1)
        if isempty(veta1)
            eu[2] = 1
        else
            eu[2] = Int64(norm(veta1 - (veta * adjoint(veta)) * veta1, 2) < ϵ * n)
        end
    end
    I, J, V = SparseArrays.spdiagm_internal(0 => [ones(n - nunstab); zeros(nunstab)])
    diag_m = sparse(I, J, V, n, n)
    G1 = real(U * T * diag_m * U')#spdiagm([ones(n - nunstab); zeros(nunstab)], 0, n, n) * U')
    F = U1[:, 1:nunstab]' * inv(U1[:, nunstab + 1:end]')
    impact = [F * Ψ[nunstab + 1:end, :]; Ψ[nunstab + 1:end, :]]
    C = real(U * c) .* ones(size(U, 1), 1)

    return G1::Matrix{Float64}, C::Matrix{Float64}, impact::Matrix{Float64}, U::Matrix, T::Matrix::Matrix, eu::Vector{Int64}
end

function new_divct(F::LinearAlgebra.GeneralizedSchur; ϵ::Float64 = sqrt(eps()) * 10)
    a, b = F.S, F.T
    n = size(a, 1)
    div = 0.001
    for i in 1:n
        if abs(a[i, i]) > ϵ
            divhat = real(b[i, i] / a[i, i])
            if (ϵ < divhat) && (divhat < div)
                div = 0.5 * divhat
            end
        end
    end
    return div
end

function decomposition_svdct!(A; ϵ::Float64 = sqrt(eps()) * 10)
    Asvd = svd!(A)
    bigev = findall(Asvd.S .> ϵ)
    Au = Asvd.U[:, bigev]
    Ad = diagm(0 => Asvd.S[bigev])
    Av = Asvd.V[:, bigev]
    return bigev, Au, Ad, Av
end

# Additional cases when inputs are sparse matrices
function gensysct(Γ0::SparseMatrixCSC{Float64,Int64}, Γ1::SparseMatrixCSC{Float64,Int64},
                  c::SparseMatrixCSC{Float64,Int64}, Ψ::SparseMatrixCSC{Float64,Int64},
                  Π::SparseMatrixCSC{Float64,Int64};
                  check_existence::Bool = true, check_uniqueness::Bool = true, ϵ::Float64 = sqrt(eps()) * 10,
                  div::Float64 = -1.0, complex_decomposition::Bool = false)
    F = complex_decomposition ? schurfact(complex(Matrix{Float64}(Γ0)), complex(Matrix{Float64}(Γ1))) : schurfact(real(Matrix{Float64}(Γ0)), real(Matrix{Float64}(Γ1)))
    gensysct(F, Matrix{Float64}(c), Matrix{Float64}(Ψ), Matrix{Float64}(Π); check_existence = check_existence,
             check_uniqueness = check_uniqueness, ϵ = ϵ, div = div)
end

function gensysct!(Γ0::SparseMatrixCSC{Float64,Int64}, Γ1::SparseMatrixCSC{Float64,Int64},
                   c::SparseMatrixCSC{Float64,Int64}, Ψ::SparseMatrixCSC{Float64,Int64},
                   Π::SparseMatrixCSC{Float64,Int64};
                  check_existence::Bool = true, check_uniqueness::Bool = true, ϵ::Float64 = sqrt(eps()) * 10,
                  div::Float64 = -1.0, complex_decomposition::Bool = false)
    F = complex_decomposition ? schur!(complex(Matrix{Float64}(Γ0)), complex(Matrix{Float64}(Γ1))) : schur!(real(Matrix{Float64}(Γ0)), real(Matrix{Float64}(Γ1)))
    gensysct(F, Matrix{Float64}(c), Matrix{Float64}(Ψ), Matrix{Float64}(Π); check_existence = check_existence,
             check_uniqueness = check_uniqueness, ϵ = ϵ, div = div)
end

function gensysct(Γ1::SparseMatrixCSC{Float64,Int64}, c::SparseMatrixCSC{Float64,Int64},
                  Ψ::SparseMatrixCSC{Float64,Int64}, Π::SparseMatrixCSC{Float64,Int64};
                  check_existence::Bool = true, check_uniqueness::Bool = true, ϵ::Float64 = sqrt(eps()) * 10,
                  div::Float64 = -1.0, complex_decomposition::Bool = false)
    # Assumes that Γ0 is the identity
    F = complex_decomposition ? schurfact(complex(Matrix{Float64}(Γ1))) : schurfact(real(Matrix{Float64}(Γ1)))
    gensysct(F, Matrix{Float64}(c), Matrix{Float64}(Ψ), Matrix{Float64}(Π); check_existence = check_existence,
             check_uniqueness = check_uniqueness, ϵ = ϵ, div = div)
end

function gensysct!(Γ1::SparseMatrixCSC{Float64,Int64}, c::SparseMatrixCSC{Float64,Int64},
                   Ψ::SparseMatrixCSC{Float64,Int64}, Π::SparseMatrixCSC{Float64,Int64};
                  check_existence::Bool = true, check_uniqueness::Bool = true, ϵ::Float64 = sqrt(eps()) * 10,
                  div::Float64 = -1.0, complex_decomposition::Bool = false)
    # Assumes that Γ0 is the identity
    F = complex_decomposition ? schur!(complex(Matrix{Float64}(Γ1))) : schur!(real(Matrix{Float64}(Γ1)))
    gensysct(F, Matrix{Float64}(c), Matrix{Float64}(Ψ), Matrix{Float64}(Π); check_existence = check_existence,
             check_uniqueness = check_uniqueness, ϵ = ϵ, div = div)
end


text:
# # Tutorial 04 Perturbation Theory
#
# ## Introduction
#
# This tutorial demonstrates how perturbation theory is utilized using the
# WavesAndEigenvalues package. You may ask: 'What is perturbation theory?'
# Well, perturbation theory deals with the approximation of a solution of a
# mathematical problem based on a *known* solution of a problem similar to the
# problem of interest. (Puhh, that was a long sentence...) The problem is said
# to be perturbed from a baseline solution.
# You can find a comprehensive presentation of the internal algorithms in
# [1,2,3].
## #jl
# ## Model set-up and baseline-solution
# The model is the same Rijke tube configuration as in Tutorial 01:
using WavesAndEigenvalues.Helmholtz
mesh=Mesh("Rijke_mm.msh",scale=0.001) #load mesh
dscrp=Dict() #initialize model descriptor
dscrp["Interior"]=(:interior, ()) #define resonant cavity
dscrp["Outlet"]=(:admittance, (:Y,1E15)) #specify outlet BC
γ=1.4 #ratio of specific heats
ρ=1.225 #density at the reference location upstream to the flame in kg/m^3
Tu=300.0    #K unburnt gas temperature
Tb=1200.0    #K burnt gas temperature
P0=101325.0 # ambient pressure in Pa
A=pi*0.025^2 # cross sectional area of the tube
Q02U0=P0*(Tb/Tu-1)*A*γ/(γ-1) #the ratio of mean heat release to mean velocity Q02U0
x_ref=[0.0; 0.0; -0.00101] #reference point
n_ref=[0.0; 0.0; 1.00] #directional unit vector of reference velocity
n=0.01 #interaction index
τ=0.001 #time delay
dscrp["Flame"]=(:flame,(γ,ρ,Q02U0,x_ref,n_ref,:n,:τ,n,τ)) #flame dynamics
R=287.05 # J/(kg*K) specific gas constant (air)
speedofsound(x,y,z) = z<0. ? sqrt(γ*R*Tu) : sqrt(γ*R*Tb)
c=generate_field(mesh,speedofsound)
L=discretize(mesh,dscrp,c)

# To obtain a baseline solution we solve it using the housholder iteration with
# a very small stopping tollerance of `tol=1E-11`.
sol,nn,flag=householder(L,340*2*pi,maxiter=20,tol=1E-11)

# this small tolerance is necessary because as the name suggests the base-line
# solution is the basis to our subsequent steps. If it is inaccurate we will
# definetly also encounter inaccurate approximations to other configurations
# than the base-line set-up.
## #jl
# ## Taylor series
#
# Probably, you are somewhat familiar with the concept of Taylor series
# which is a basic example of perturbation theory. There is some function
# f from which the value f(x0) is known together with some derivatives at the
# same point, the first N say. Then, we can approximate f(x0+Δ) as:
#
# f(x0+Δ)≈∑_n=0^N f_n(x0)Δ^n

# Let's consider the time delay `τ`. Our problem depends exponentially on this
# value.  We can utilize a fast perturbation algorithm to compute the first
# 20 Taylor-series coefficients of the eigenfrequency `ω` w.r.t. `τ` by just
# typing
perturb_fast!(sol,L,:τ,20)
# The output shows you how long it takes to compute the coefficients.
# Obviously, it takes longer and longer for higher order coefficients.
# Note, that there is also an algorithm `perturb!(sol,L,:τ,20)` which does
# exactly the same as the fast algorithm but is slower, when it comes to high
# orders (N>5).

# Both algorithms populate a field in the solution object holding the Taylor
# coefficients.
sol.eigval_pert

# We can use these values to form the Taylor-series approximation and for
# convenience we can just do so by calling to the solution object itself.
# For instance let's assume we would like to approximate the value of the
# eigenfrequency when "τ==0.00125" based on the 20th order series expansion.
Δ=0.00001
ω_approx=sol(:τ,τ+Δ,20)

# Let's compare this value to the true solution:
L.params[:τ]=τ+Δ #change parameter
sol_exact,nn,flag=householder(L,340*2*pi,maxiter=20,tol=1E-11) #solve again
ω_exact=sol_exact.params[:ω]
println(" exact=$(ω_exact/2/pi)  vs  approx=$(ω_approx/2/pi))")
# Clearly, the approximation matches the first 4 digits after the point!

# We can also compute the approximation at any lower order than 20.
# For instance, the first-order approximation is:
ω_approx=sol(:τ,τ+Δ,1)
println(" first-order approx=$(ω_approx/2/pi)")

# Note that the accuracy is less than at twentieth order.
#
# Also note that we cannot compute the perturbation at higher order than 20
# because we have only computed the Taylor series coefficients up to 20th order.
# Of course we could, prepare higher order coefficients, let's say up to 30th
# order by
perturb_fast!(sol,L,:τ,30)
# and then
ω_approx=sol(:τ,τ+Δ,30)
println(" 30th-order approx=$(ω_approx/2/pi)")

# Indeed, `30` is a special limit because per default the WavesAndEigenvalues
# package installs the perturbation algorithm on your machine up to 30th order.
# This is because higher orders would consume significantly more memory in
# the installation directory and also the computation of higher orders may take
# some time.  However, if you really want to use higher orders. You can rebuild
# your package by first setting a special environment variable to your desired
# order -- say 40 -- by `ENV["JULIA_WAE_PERT_ORDER"]=40` and then run
# `Pkg.build("WavesAndEigenvalues")`
# Note, that if you don't make the environment variable permanent,
# these steps will be necessary once every time you got any update to the
# package from julias package manager.

#TODO:speed

## #jl
# Ok, we've seen how we can compute Taylor-series coefficients and how to
# evaluate them as a Taylor series. But how do we choose parameters like
# the baseline set-up or the eprturbation order to get a reasonable estimate?
# Well there is a lot you can do but, unfornately, there are a lot of
# misunderstandings when it comes to the quality perturbative approximations.
#
# Take a second and try to answer the following questions:
# 1. What is the range of Δ in which you can expect reliable results from
# the approximation, i.e., the difference from the true solution is small?
# 2. How does the quality of your approximation improve if you increase the
# number of known derivatives?
# 3. What else can you do to improve your solution?
#
# Regarding the first question, many people misbelieve that the approximation is
# good as long as Δ (the shift from the expansion point) is small. This is right
# and wrong at the same time. Indeed, to classify Δ as small, you need to
# compare it agains some value. For Taylor series this is the radius of
# convergence. Convergence radii can be quite huge and sometimes super small.
# Also keep in mind that the nuemrical value of Δ in most engineering
# applications is meaningless if it is not linked to some unit of measure.
# (You know the joke: What is larger, 1 km or a million mm?)

# So how do we get the radius of convergence? It's as simple like
r=conv_radius(sol,:τ)
#The return value here is an array that holds N-1 values, where N is the order
# of Taylor-series coefficients that have been computed for `:τ`. This is
# because the convergenvce radius is estimated based on the ratio of two
# consecutive Taylor-series coefficients. Usually, the last entry of r should
# give you the best approximate.
println("Best estimate available: r=$(r[end])")
# However, there are cases where the estimation procedure for the convergence
# radius is not appropriate. That is why you can have a look on the other entries
# in r to judge whether there is smooth convergence.
## #jl
# Let's see what's happening if we evaluate the Taylor series beyond it's radius
# of convergence.
ω_approx=sol(:τ,τ+r[end]+0.001,20)
L.params[:τ]=τ+r[end]+0.001
sol_exact,nn,flag=householder(L,340*2*pi,maxiter=20,tol=1E-11) #solve again
ω_exact=sol_exact.params[:ω]
println(" exact=$(ω_exact/2/pi)  vs  approx=$(ω_approx/2/pi))")
# Outch, the approximation is completely off.
# Remember question 2? Maybe the estimate gets better if we increase the
# perturbation order ? At least the last time we've seen an improvement.
# But this time...
ω_approx=sol(:τ,τ+r[end]+0.001,30)
println(" exact=$(ω_exact/2/pi)  vs  approx=$(ω_approx/2/pi))")
# things get *way* worse! This is exaclty the reason why some clever person
# named r the radius of *convergence*. Beyond that radius we cannot make the
# Taylor-series converge without shifting the expansion point. You might think:
# 'Well, then let's shift the expansion point!' While this is definitely a valid
# approach, there is something better you can do...

# ## Series accelartion by Padé approximation
#
# The Taylor-series cannot converge beyond its radius of convergence. So why
# not trying someting else than a Taylor-series?. The truncated Taylor series is
# a polynomial approximation to our unknown relation ω=ω(τ). An alternative
# (if not to say a generalazation) of this approach is to consider rational
# polynomial approximations.  So isntead of
# f(x0+Δ)≈∑_n=0^N f_n(x0)Δ^n
# we try something like
# f(x0+Δ)≈[ ∑_l=0^L a_l(x0)Δ^l ] / [ 1 + ∑_m=0^M b_m(x0)Δ^m ]
# In order to get the same asymptotic behaviour close to our expansion point, we
# demand that the truncated Taylor series expansion of our Padé approximant is
# identical to the approximation of our unknown function. Here comes the clou:
# we already know these values because we computed the Taylor-series
# coefficients. Only little algebra is needed to convert the Taylor coefficients
# to Padé coefficients and all of this is build in to the solution type. All you
# need to know is that the number of Padé coefficients is related to the number
# of Taylor coefficients by the simple formula L+M=N. Let's give it a try and
# compute the Padé approximant for L=M=10 (a so-called diagonal Padé
# approximant). This is just achieve by an extra argument for the call to our
# solution object.
ω_approx=sol(:τ,τ+r[end]+0.001,10,10)
println(" exact=$(ω_exact/2/pi)  vs  approx=$(ω_approx/2/pi))")
# Wow! There is now only a diffrence of about 1hz between the true solution and
# it's estimate. I'would say that's fine for many egineering applications.
# What's your opinion?
## #jl
# ## Conclusion
# You learned how to use perturbation theory. The main computational costs for
# this approach is the computation of Taylor-series coefficients. Once you have
# them everything comes down to ta few evaluations of scalar polynomials.
# Especially, when you like to rapidly evalaute your model for a bunch of values
# in a given parameter range. This approach should speed up your computations.
#
# The next tutorial will familiarize you with the use of higher order
# finite elements.

ext: jl


text:
# This file was generated, do not modify it. # hide
by(iris, :Species, MPL = :PetalLength => mean, SPL = :PetalLength => std)


text:
# This file calculates the solution to Laplace's equation via the adaptive spectral method.
# Δu = 0,
# u|Γ = u_0,
# u^i = log|z-z_s|,
# u = u^i + u^s.
# The normal derivative ∂u/∂n of the solution is calculated on the constant-charge boundary Γ.
# The reflected solution is calculated by convolving ∂u/∂n with the fundamental solution.
# Then, the total solution is obtained by summing the incident and the reflected solutions.

using ApproxFun, SingularIntegralEquations

z_s = 2.0
ui(x,y) = logabs(complex(x,y)-z_s)
g1(x,y) = 1/2

# Set the domains.
N = 10
r = 1e-1
cr = exp(im*2π*(0:N-1)/N)
crl = (1-2im*r)cr
crr = (1+2im*r)cr
dom = ∪(Interval,crl,crr) # All infinitesimal plates
#dom = ∪(Circle,cr,ones(length(cr))r) # All wires
#dom = ∪(Interval,crl[1:2:end],crr[1:2:end]) ∪ ∪(Circle,cr[2:2:end],ones(length(cr[2:2:end]))r) # Interlaced wires and plates

sp = Space(dom)
cwsp = CauchyWeight(sp⊗sp,0)
uiΓ,⨍ = Fun(t->ui(real(t),imag(t))+0im,sp),DefiniteLineIntegral(dom)

@time G = GreensFun(g1,cwsp;method=:Cholesky)

@time φ0,∂u∂n=[0 ⨍;1 ⨍[G]]\Any[0.,uiΓ]

println("The length of ∂u∂n is: ",length(∂u∂n))

us(x,y) = -logkernel(∂u∂n,complex(x,y))/2
ut(x,y) = ui(x,y) + us(x,y)
println("This is the approximate gradient: ",((ut(1e-5,0.)-ut(-1e-5,0.))/2e-5))


text:
#This script calculates the squeezing for different quadratures.
using LinearAlgebra
using Plots; pyplot()
using NPZ
using LaTeXStrings

include("MagnusExpansion.jl")
include("boson.jl")
include("HarmonicCorrectionBoson.jl")

"""
Calculate the uncertainties on the x(θ) and y(θ) quadratures for
a range imax values of θ in the interval [θi, θi+θrange].
"""
function quadrature_scan(moments, imax; θi=-pi/4, θrange=pi/2)
	θ_vec, Δx2_vec, Δy2_vec = zeros(imax), zeros(imax), zeros(imax)
	for i = 1: imax
		θ = θi + (i - 1)*θrange/(imax - 1)
		θ_vec[i] = θ
		Δx2_vec[i] = (moments[3] + cos(2*θ)*moments[1] + sin(2*θ)*moments[2])/2
		Δy2_vec[i] = (moments[3] - cos(2*θ)*moments[1] - sin(2*θ)*moments[2])/2
	end
	return θ_vec, Δx2_vec, Δy2_vec
end
"""
Apply the function quadrature_scan a couple of times, always around the value of θ
where maximum squeezing has been found.
"""
function fine_quadrature_scan(moments, imax)
	θ_vec, Δx2_vec, Δy2_vec = quadrature_scan(moments, imax)
	for i = 1: 4
		k = argmin(Δy2_vec)
		θi = θ_vec[k - 1]
		θrange = θ_vec[k + 1] - θ_vec[k - 1]
		output = quadrature_scan(moments, imax; θi=θi, θrange=θrange)
		θ_vec .= output[1]
		Δx2_vec .= output[2]
		Δy2_vec .= output[3]
	end
	return θ_vec, Δx2_vec, Δy2_vec
end

#array with different gate times tf
tf0, tf1 = 2.0, 15.0
nt = 200
tf_vec = zeros(nt)
nt = size(tf_vec)[1]
dB_vec = zeros(nt, 3)
moments = zeros(nt, 3)
coeffs_array = zeros(nt, 3)

imax = 10
θopt = zeros(nt,2)

for j = 1: nt
	println(j)
	tf = tf0*(tf1/tf0)^((j - 1)/(nt - 1))
	tf_vec[j] = tf
	κ = 0.0
	# solution in the absence of counter-rotating terms
	cr_terms = false
	moments[j, :] = solve_2moments(tf, κ, cr_terms, zeros(3))
	Δx2 = (moments[j, 3] + moments[j, 1])/2
	Δy2 = (moments[j, 3] - moments[j, 1])/2
	dB_vec[j, 1] = -10.0*log10(2*Δy2)

	# solution in the presence of counter-rotating terms
	cr_terms = true
	moments[j, :] = solve_2moments(tf, κ, cr_terms, zeros(3))
	Δx2 = (moments[j, 3] + moments[j, 1])/2
	Δy2 = (moments[j, 3] - moments[j, 1])/2
	dB_vec[j, 2] = -10.0*log10(2*Δy2)

	θ_vec, Δx2_vec, Δy2_vec = fine_quadrature_scan(moments[j, :], imax)
	k = argmin(Δy2_vec)
	θopt[j, 1] = θ_vec[k]

	# solution in the presence of counter-rotating terms and correction.
	cr_terms = true
	M = correction_matrix(tf)
	order = 6
	coeffs = get_coeffs(M, order, tf)

	coeffs_array[j, :] = coeffs
	moments[j, :] = solve_2moments(tf, κ, cr_terms, coeffs)
	Δx2 = (moments[j, 3] + moments[j, 1])/2
	Δy2 = (moments[j, 3] - moments[j, 1])/2
	dB_vec[j, 3] = -10.0*log10(2*Δy2)

	θ_vec, Δx2_vec, Δy2_vec = fine_quadrature_scan(moments[j, :], imax)
	k = argmin(Δy2_vec)
	θopt[j, 2] = θ_vec[k]
end
# save figures
plot(tf_vec, dB_vec,
	 linewidth = 2,
	 labels = ["No correction" "6th order\nMagnus" "RWA"])
plot!(tickfont = font(12, "Serif"),
	  xlabel = "Pulse width "*L"\omega_\mathrm{a} t_\mathrm{f}",
	  ylabel = "Squeezing (dB)",
	  guidefont = font(14, "Serif"),
	  legendfont = font(12, "Serif"),
	  background_color_legend = false,
	  margin = 5Plots.mm)
savefig("squeezing.pdf")

y = hcat(θopt, zeros(size(θopt)[1]))
plot(tf_vec, y,
	 linewidth = 2,
	 labels = ["No correction" "6th order\nMagnus" "RWA"])
plot!(tickfont = font(12, "Serif"),
	  xlabel = "Pulse width "*L"\omega_\mathrm{a} t_\mathrm{f}",
	  ylabel = L"\varphi",
	  guidefont = font(14, "Serif"),
	  legendfont = font(12, "Serif"),
	  background_color_legend = false,
	  margin = 5Plots.mm)
savefig("squeezing_angle.pdf")

y = hcat(1.0./tf_vec, abs.(coeffs_array))
plot(tf_vec, y,
	 linewidth = 2,
	 labels = [L"1/t_\mathrm{f}" L"|c_{x,1}|" L"|c_{y,1}|" L"|\Delta|"])
plot!(yaxis = :log,
	  tickfont = font(12, "Serif"),
	  xlabel = "Pulse width "*L"\omega_\mathrm{a} t_\mathrm{f}",
	  ylabel = "coeff."*L"\times \omega_\mathrm{a}^{-1}",
	  guidefont = font(14, "Serif"),
	  legendfont = font(12, "Serif"),
	  background_color_legend = false,
	  margin = 5Plots.mm)
savefig("coeffs.pdf")

# save data files
output = zeros(nt,4)
output[:, 1] = tf_vec
output[:, 2: end] = dB_vec
npzwrite("dB.npy", output)
output = zeros(nt,3)
output[:, 1] = tf_vec
output[:, 2: end] = θopt
npzwrite("angle.npy", output)
output = zeros(nt,4)
output[:, 1] = tf_vec
output[:, 2: end] = coeffs_array
npzwrite("coeffs.npy", output)

max_stars_repo_name: thalesfr/MagnusCorrection

max_issues_repo_name: thalesfr/MagnusCorrection

max_forks_repo_name: thalesfr/MagnusCorrection


text:
using Pkg
using Pkg: TOML
using Pkg: Operations, Types
using UUIDs

#=
genfile & create_project_from_require have been taken from the PR
https://github.com/JuliaLang/PkgDev.jl/pull/144
which was created by https://github.com/KristofferC

THIS IS JUST A TEMPORARY SOLUTION FOR PACKAGES WITHOUT A TOML AND WILL GET MOVED OUT!
=#

function packages_from_require(reqfile::String)
    ctx = Pkg.Types.Context()
    pkgs = Types.PackageSpec[]
    compatibility = Pair{String, String}[]
    for r in Pkg.Pkg2.Reqs.read(reqfile)
        r isa Pkg.Pkg2.Reqs.Requirement || continue
        r.package == "julia" && continue
        push!(pkgs, Types.PackageSpec(r.package))
        intervals = r.versions.intervals
        if length(intervals) != 1
            @warn "Project.toml creator cannot handle multiple requirements for $(r.package), ignoring"
        else
            l = intervals[1].lower
            h = intervals[1].upper
            if l != v"0.0.0-"
                # no upper bound
                if h == typemax(VersionNumber)
                    push!(compatibility, r.package => string(">=", VersionNumber(l.major, l.minor, l.patch)))
                else # assume semver
                    push!(compatibility, r.package => string(">=", VersionNumber(l.major, l.minor, l.patch), ", ",
                                                             "<", VersionNumber(h.major, h.minor, h.patch)))
                end
            end
        end
    end
    Operations.registry_resolve!(ctx.env, pkgs)
    Operations.ensure_resolved(ctx.env, pkgs)
    pkgs
end
function create_project_from_require(pkgname::String, path::String, toml_path::String)
    ctx = Pkg.Types.Context()
    # Package data
    path = abspath(path)
    mainpkg = Types.PackageSpec(pkgname)
    Pkg.Operations.registry_resolve!(ctx.env, [mainpkg])
    if !Operations.has_uuid(mainpkg)
        uuid = UUIDs.uuid1()
        @info "Unregistered package $pkgname, giving it a new UUID: $uuid"
        mainpkg.version = v"0.1.0"
    else
        uuid = mainpkg.uuid
        @info "Registered package $pkgname, using already given UUID: $(mainpkg.uuid)"
        Pkg.Operations.set_maximum_version_registry!(ctx.env, mainpkg)
        v = mainpkg.version
        # Remove the build
        mainpkg.version = VersionNumber(v.major, v.minor, v.patch)
    end
    # Dependency data
    dep_pkgs = Types.PackageSpec[]
    test_pkgs = Types.PackageSpec[]
    compatibility = Pair{String, String}[]

    reqfiles = [joinpath(path, "REQUIRE"), joinpath(path, "test", "REQUIRE")]
    for (reqfile, pkgs) in zip(reqfiles, [dep_pkgs, test_pkgs])
        if isfile(reqfile)
            append!(pkgs, packages_from_require(reqfile))
        end
    end

    stdlib_deps = Pkg.Operations.find_stdlib_deps(ctx, path)
    for (stdlib_uuid, stdlib) in stdlib_deps
        pkg = Types.PackageSpec(stdlib, stdlib_uuid)
        if stdlib == "Test"
            push!(test_pkgs, pkg)
        else
            push!(dep_pkgs, pkg)
        end
    end

    # Write project

    project = Dict(
        "name" => pkgname,
        "uuid" => string(uuid),
        "version" => string(mainpkg.version),
        "deps" => Dict(pkg.name => string(pkg.uuid) for pkg in dep_pkgs)
    )

    if !isempty(compatibility)
        project["compat"] =
            Dict(name => ver for (name, ver) in compatibility)
    end

    if !isempty(test_pkgs)
        project["extras"] = Dict(pkg.name => string(pkg.uuid) for pkg in test_pkgs)
        project["targets"] = Dict("test" => [pkg.name for pkg in test_pkgs])
    end

    open(toml_path, "w") do io
        Pkg.TOML.print(io, project, sorted=true, by=key -> (Types.project_key_order(key), key))
    end
end

function package_toml(package::Symbol)
    pstr = string(package)
    # could use eval using here?! Not sure what is actually better
    pkg_module = Base.require(Module(), package)
    pkg_root = normpath(joinpath(dirname(pathof(pkg_module)), ".."))
    toml = joinpath(pkg_root, "Project.toml")
    # Check for snoopfile and fall back to runtests.jl
    # if it can't be found
    snoopfile = get_snoopfile(pkg_root)
    # We will create a new toml, based that will include all test dependencies etc
    # We're also using the precompile toml as a temp toml for packages not having a toml
    precompile_toml = package_folder(pstr, "Project.toml")
    isdir(dirname(precompile_toml)) || mkpath(dirname(precompile_toml))
    test_deps = Dict()
    if !isfile(toml)
        create_project_from_require(pstr, pkg_root, precompile_toml)
    else
        testreq = joinpath(pkg_root, "test", "REQUIRE")
        if isfile(testreq)
            pkgs = packages_from_require(testreq)
            test_deps = Dict(pkg.name => string(pkg.uuid) for pkg in pkgs)
        end
        cp(toml, precompile_toml, force = true)
        chmod(precompile_toml, 0o644)
    end
    # remove any old manifest
    # if isfile(package_folder(pstr, "Manifest.toml"))
    #     rm(package_folder(pstr, "Manifest.toml"))
    # end
    # add ourselves as dependencies and ensure we have a manifest
    # run_julia("""
    # using Pkg
    # Pkg.instantiate()
    # pkg"add PackageCompiler Pkg"
    # """, project = precompile_toml)

    toml = TOML.parsefile(precompile_toml)

    deps = merge(get(toml, "deps", Dict()), test_deps)
    # Add the package itself
    deps[toml["name"]] = toml["uuid"]
    # Add the packages we need
    test_deps = get(toml, "extras", Dict())
    compile_toml = Dict()
    compile_toml["name"] = string(package, "Precompile")
    compile_toml["deps"] = merge(test_deps, deps)
    if haskey(toml, "compat")
        compile_toml["compat"] = toml["compat"]
    end
    write_toml(precompile_toml, compile_toml)
    precompile_toml, snoopfile
end

function write_toml(path, dict)
    open(path, "w") do io
        TOML.print(
            io, dict,
            sorted = true, by = key-> (Types.project_key_order(key), key)
        )
    end
end

function get_snoopfile(pkg_root)
    snoopfileroot = joinpath(pkg_root,"snoopfile.jl")
    snoopfilesnoopdir = joinpath(pkg_root, "snoop", "snoopfile.jl")
    snoopfilesrc = joinpath(pkg_root, "src", "snoopfile.jl")
    if isfile(snoopfileroot)
        snoopfile = snoopfileroot
    elseif isfile(snoopfilesnoopdir)
        snoopfile = snoopfilesnoopdir
    elseif isfile(snoopfilesrc)
        snoopfile = snoopfilesrc
    else
        snoopfile = joinpath(pkg_root, "test", "runtests.jl")
    end
    return snoopfile
end


text:
using Test, EllZeroTrendFiltering
using Interpolations

Random.seed!(31415)

# Auxilliary function for evaluating the approximation error for a specific
# piecewise linear approximation (I, Y), given the transition costs ℓ
function compute_cost(ℓ, I, Y)
    cost = 0.0;
    for k=1:length(I)-1
     cost += ℓ[I[k], I[k+1]](Y[k], Y[k+1])
    end
    return cost
end

## Test 1
# Compare the approximation error for specific piecewise linear approximations
# (I, Y) when the costs are computed  by
# (1) using ℓ computed by compute_discrete_transition_costs and
#     then evaluated using the above auxilliary function
# (2) simple linear interpolation

g1 = sin.(range(0, stop=π, length=50))
I1 = [1,10,30,40,50]
Y1 = g1[I1]
t_ss1 = [1, 10, 21, 33, 50]

g2 = g1
I2 = [1,15,30,36,50]
Y2 = randn(5)
t_ss2 = t_ss1

g3 = randn(20)
I3 = [1; 3:2:19; 20]
Y3 = randn(length(I3))
t_ss3 = [1,5,10,15,20]

g4 = g3
I4 = 1:length(g3)
Y4 = g4
t_ss4 = t_ss3

for (g, I, Y, t_subsampled) in ((g1, I1, Y1, t_ss1), (g2, I2, Y2, t_ss2), (g3, I3, Y3, t_ss3))#, (g4, I4, Y4))
    ℓ = compute_discrete_transition_costs(g)
    cost1 = compute_cost(ℓ,I,Y)

    y = interpolate((I,), Y, Gridded(Linear()))(1:length(g))
    cost2 = sum((y[1:end-1]-g[1:end-1]).^2) # Note: cost at i=N should not be included

    @test cost1 ≈ cost2

    # Test transition cost with t grid
    ℓ_subsampled = compute_discrete_transition_costs(g, t_subsampled)

    for i in 1:length(t_subsampled)-1
        for ip in (i+1):length(t_subsampled)
            t_i = t_subsampled[i]
            t_ip = t_subsampled[ip]
            @test ℓ_subsampled[i, ip] == ℓ[t_i, t_ip]
        end
    end
end

## Test 2
# Check the transition costs for a very simple example
ℓ = compute_discrete_transition_costs([1.0, 2.0, 3.0])

# (y - 1)^2
@test ℓ[1,2].P == [1.0 0.0; 0.0 0.0]
@test ℓ[1,2].q == [-2.0, 0]
@test ℓ[1,2].r == 1

# (y - 1)^2 + (y/2 + y'/2 - 2)^2
@test ℓ[1,3].P == [1.25 0.25; 0.25 0.25]
@test ℓ[1,3].q == [-4.0, -2.0]
@test ℓ[1,3].r == 5.0


text:
"""
    horner(c,x)

Evaluate a polynomial whose coefficients are given in descending order in `c`,
at the point `x`, using Horner's rule.
"""
function horner(c,x)

    n = length(c)
    y = c[1]
    for k = 2:n
        y = x*y + c[k]
    end

    return y

end


text:
module KissMCMC

using ProgressMeter, StaticArrays

# std-lib:
using Statistics

export metropolis, emcee, squash_walkers, make_theta0s

include("samplers.jl")
include("analysis.jl")

end # module


text:
# Automatically generated using Clang.jl

const HID_API_VERSION_MAJOR = 0
const HID_API_VERSION_MINOR = 10
const HID_API_VERSION_PATCH = 0

struct hid_api_version
    major::Cint
    minor::Cint
    patch::Cint
end

const hid_device_ = Cvoid
const hid_device = hid_device_

struct hid_device_info
    path::Cstring
    vendor_id::UInt16
    product_id::UInt16
    serial_number::Ptr{Cwchar_t}
    release_number::UInt16
    manufacturer_string::Ptr{Cwchar_t}
    product_string::Ptr{Cwchar_t}
    usage_page::UInt16
    usage::UInt16
    interface_number::Cint
    next::Ptr{hid_device_info}
end


text:
module HydrostaticBoussinesq

export HydrostaticBoussinesqModel, AbstractHydrostaticBoussinesqProblem

using StaticArrays
using LinearAlgebra: dot, Diagonal
using CLIMAParameters.Planet: grav

using ..VariableTemplates
using ..MPIStateArrays
using ..Mesh.Filters: apply!
using ..Mesh.Grids: VerticalDirection
using ..DGMethods:
    BalanceLaw,
    LocalGeometry,
    DGModel,
    indefinite_stack_integral!,
    reverse_indefinite_stack_integral!,
    nodal_update_auxiliary_state!,
    copy_stack_field_down!
using ..DGMethods.NumericalFluxes: RusanovNumericalFlux

import ..DGMethods.NumericalFluxes: update_penalty!
import ..DGMethods:

    vars_state_conservative,
    init_state_conservative!,
    vars_state_auxiliary,
    init_state_auxiliary!,
    vars_state_gradient,
    compute_gradient_argument!,
    vars_state_gradient_flux,
    compute_gradient_flux!,
    vars_integrals,
    integral_load_auxiliary_state!,
    integral_set_auxiliary_state!,
    vars_reverse_integrals,
    reverse_integral_load_auxiliary_state!,
    reverse_integral_set_auxiliary_state!,
    flux_first_order!,
    flux_second_order!,
    source!,
    wavespeed,
    update_auxiliary_state!,
    update_auxiliary_state_gradient!

×(a::SVector, b::SVector) = StaticArrays.cross(a, b)
⋅(a::SVector, b::SVector) = StaticArrays.dot(a, b)

abstract type AbstractHydrostaticBoussinesqProblem end

"""
    HydrostaticBoussinesqModel <: BalanceLaw

A `BalanceLaw` for ocean modeling.

write out the equations here

ρₒ = reference density of sea water
cʰ = maximum horizontal wave speed
cᶻ = maximum vertical wave speed
αᵀ = thermal expansitivity coefficient
νʰ = horizontal viscosity
νᶻ = vertical viscosity
κʰ = horizontal diffusivity
κᶻ = vertical diffusivity
fₒ = first coriolis parameter (constant term)
β  = second coriolis parameter (linear term)

# Usage

    HydrostaticBoussinesqModel(problem)

"""
struct HydrostaticBoussinesqModel{PS, P, T} <: BalanceLaw
    param_set::PS
    problem::P
    ρₒ::T
    cʰ::T
    cᶻ::T
    αᵀ::T
    νʰ::T
    νᶻ::T
    κʰ::T
    κᶻ::T
    fₒ::T
    β::T
    function HydrostaticBoussinesqModel{FT}(
        param_set::PS,
        problem;
        ρₒ = FT(1000),  # kg / m^3
        cʰ = FT(0),     # m/s
        cᶻ = FT(0),     # m/s
        αᵀ = FT(2e-4),  # (m/s)^2 / K
        νʰ = FT(5e3),   # m^2 / s
        νᶻ = FT(5e-3),  # m^2 / s
        κʰ = FT(1e3),   # m^2 / s
        κᶻ = FT(1e-4),  # m^2 / s
        fₒ = FT(1e-4),  # Hz
        β = FT(1e-11), # Hz / m
    ) where {FT <: AbstractFloat, PS}
        return new{PS, typeof(problem), FT}(
            param_set,
            problem,
            ρₒ,
            cʰ,
            cᶻ,
            αᵀ,
            νʰ,
            νᶻ,
            κʰ,
            κᶻ,
            fₒ,
            β,
        )
    end
end
HBModel = HydrostaticBoussinesqModel

"""
    vars_state_conservative(::HBModel)

prognostic variables evolved forward in time

u = (u,v) = (zonal velocity, meridional velocity)
η = sea surface height
θ = temperature
"""
function vars_state_conservative(m::HBModel, T)
    @vars begin
        u::SVector{2, T}
        η::T # real a 2-D variable TODO: should be 2D
        θ::T
    end
end

"""
    init_state_conservative!(::HBModel)

sets the initial value for state variables
dispatches to ocean_init_state! which is defined in a problem file such as SimpleBoxProblem.jl
"""
function ocean_init_state! end
function init_state_conservative!(m::HBModel, Q::Vars, A::Vars, coords, t)
    return ocean_init_state!(m.problem, Q, A, coords, t)
end

"""
    vars_state_auxiliary(::HBModel)
helper variables for computation

second half is because there is no dedicated integral kernels
these variables are used to compute vertical integrals
w = vertical velocity
wz0 = w at z = 0
pkin = bulk hydrostatic pressure contribution

first half of these are fields that are used for computation
y = north-south coordinate

"""
function vars_state_auxiliary(m::HBModel, T)
    @vars begin
        y::T     # y-coordinate of the box
        w::T     # ∫(-∇⋅u)
        pkin::T  # ∫(-αᵀθ)
        wz0::T   # w at z=0
    end
end

"""
    init_state_auxiliary!(::HBModel)

sets the initial value for auxiliary variables (those that aren't related to vertical integrals)
dispatches to ocean_init_aux! which is defined in a problem file such as SimpleBoxProblem.jl
"""
function ocean_init_aux! end
function init_state_auxiliary!(m::HBModel, A::Vars, geom::LocalGeometry)
    return ocean_init_aux!(m, m.problem, A, geom)
end

"""
    vars_state_gradient(::HBModel)

variables that you want to take a gradient of
these are just copies in our model
"""
function vars_state_gradient(m::HBModel, T)
    @vars begin
        ∇u::SVector{2, T}
        ∇θ::T
    end
end

"""
    compute_gradient_argument!(::HBModel)

copy u and θ to var_gradient
this computation is done pointwise at each nodal point

# arguments:
- `m`: model in this case HBModel
- `G`: array of gradient variables
- `Q`: array of state variables
- `A`: array of aux variables
- `t`: time, not used
"""
@inline function compute_gradient_argument!(m::HBModel, G::Vars, Q::Vars, A, t)
    G.∇u = Q.u
    G.∇θ = Q.θ

    return nothing
end

"""
    vars_state_gradient_flux(::HBModel)

the output of the gradient computations
multiplies ∇u by viscosity tensor and ∇θ by the diffusivity tensor
"""
function vars_state_gradient_flux(m::HBModel, T)
    @vars begin
        ν∇u::SMatrix{3, 2, T, 6}
        κ∇θ::SVector{3, T}
    end
end

"""
    compute_gradient_flux!(::HBModel)

copy ∇u and ∇θ to var_diffusive
this computation is done pointwise at each nodal point

# arguments:
- `m`: model in this case HBModel
- `D`: array of diffusive variables
- `G`: array of gradient variables
- `Q`: array of state variables
- `A`: array of aux variables
- `t`: time, not used
"""
@inline function compute_gradient_flux!(
    m::HBModel,
    D::Vars,
    G::Grad,
    Q::Vars,
    A::Vars,
    t,
)
    ν = viscosity_tensor(m)
    D.ν∇u = -ν * G.∇u

    κ = diffusivity_tensor(m, G.∇θ[3])
    D.κ∇θ = -κ * G.∇θ

    return nothing
end

"""
    viscosity_tensor(::HBModel)

uniform viscosity with different values for horizontal and vertical directions

# Arguments
- `m`: model object to dispatch on and get viscosity parameters
"""
@inline viscosity_tensor(m::HBModel) = Diagonal(@SVector [m.νʰ, m.νʰ, m.νᶻ])

"""
    diffusivity_tensor(::HBModel)

uniform diffusivity in the horizontal direction
applies convective adjustment in the vertical, bump by 1000 if ∂θ∂z < 0

# Arguments
- `m`: model object to dispatch on and get diffusivity parameters
- `∂θ∂z`: value of the derivative of temperature in the z-direction
"""
@inline function diffusivity_tensor(m::HBModel, ∂θ∂z)
    ∂θ∂z < 0 ? κ = (@SVector [m.κʰ, m.κʰ, 1000 * m.κᶻ]) : κ =
        (@SVector [m.κʰ, m.κʰ, m.κᶻ])

    return Diagonal(κ)
end

"""
    vars_integral(::HBModel)

location to store integrands for bottom up integrals
∇hu = the horizontal divegence of u, e.g. dw/dz
"""
function vars_integrals(m::HBModel, T)
    @vars begin
        ∇hu::T
        αᵀθ::T
    end
end

"""
    integral_load_auxiliary_state!(::HBModel)

copy w to var_integral
this computation is done pointwise at each nodal point

arguments:
m -> model in this case HBModel
I -> array of integrand variables
Q -> array of state variables
A -> array of aux variables
"""
@inline function integral_load_auxiliary_state!(
    m::HBModel,
    I::Vars,
    Q::Vars,
    A::Vars,
)
    I.∇hu = A.w # borrow the w value from A...
    I.αᵀθ = -m.αᵀ * Q.θ # integral will be reversed below

    return nothing
end

"""
    integral_set_auxiliary_state!(::HBModel)

copy integral results back out to aux
this computation is done pointwise at each nodal point

arguments:
m -> model in this case HBModel
A -> array of aux variables
I -> array of integrand variables
"""
@inline function integral_set_auxiliary_state!(m::HBModel, A::Vars, I::Vars)
    A.w = I.∇hu
    A.pkin = I.αᵀθ

    return nothing
end

"""
    vars_reverse_integral(::HBModel)

location to store integrands for top down integrals
αᵀθ = density perturbation
"""
function vars_reverse_integrals(m::HBModel, T)
    @vars begin
        αᵀθ::T
    end
end

"""
    reverse_integral_load_auxiliary_state!(::HBModel)

copy αᵀθ to var_reverse_integral
this computation is done pointwise at each nodal point

arguments:
m -> model in this case HBModel
I -> array of integrand variables
A -> array of aux variables
"""
@inline function reverse_integral_load_auxiliary_state!(
    m::HBModel,
    I::Vars,
    Q::Vars,
    A::Vars,
)
    I.αᵀθ = A.pkin

    return nothing
end

"""
    reverse_integral_set_auxiliary_state!(::HBModel)

copy reverse integral results back out to aux
this computation is done pointwise at each nodal point

arguments:
m -> model in this case HBModel
A -> array of aux variables
I -> array of integrand variables
"""
@inline function reverse_integral_set_auxiliary_state!(
    m::HBModel,
    A::Vars,
    I::Vars,
)
    A.pkin = I.αᵀθ

    return nothing
end

"""
    flux_first_order!(::HBModel)

calculates the hyperbolic flux contribution to state variables
this computation is done pointwise at each nodal point

# arguments:
m -> model in this case HBModel
F -> array of fluxes for each state variable
Q -> array of state variables
A -> array of aux variables
t -> time, not used

# computations
∂ᵗu = ∇⋅(g*η + g∫αᵀθdz + v⋅u)
∂ᵗθ = ∇⋅(vθ) where v = (u,v,w)
"""
@inline function flux_first_order!(
    m::HBModel,
    F::Grad,
    Q::Vars,
    A::Vars,
    t::Real,
)
    FT = eltype(Q)
    _grav::FT = grav(m.param_set)
    @inbounds begin
        u = Q.u # Horizontal components of velocity
        η = Q.η
        θ = Q.θ
        w = A.w   # vertical velocity
        pkin = A.pkin

        v = @SVector [u[1], u[2], w]
        Iʰ = @SMatrix [
            1 -0
            -0 1
            -0 -0
        ]

        # ∇h • (g η)
        F.u += _grav * η * Iʰ

        # ∇h • (- ∫(αᵀ θ))
        F.u += _grav * pkin * Iʰ

        # ∇h • (v ⊗ u)
        # F.u += v * u'

        # ∇ • (u θ)
        F.θ += v * θ
    end

    return nothing
end

"""
    flux_second_order!(::HBModel)

calculates the parabolic flux contribution to state variables
this computation is done pointwise at each nodal point

# arguments:
- `m`: model in this case HBModel
- `F`: array of fluxes for each state variable
- `Q`: array of state variables
- `D`: array of diff variables
- `A`: array of aux variables
- `t`: time, not used

# computations
∂ᵗu = -∇⋅(ν∇u)
∂ᵗθ = -∇⋅(κ∇θ)
"""
@inline function flux_second_order!(
    m::HBModel,
    F::Grad,
    Q::Vars,
    D::Vars,
    HD::Vars,
    A::Vars,
    t::Real,
)
    F.u += D.ν∇u
    F.θ += D.κ∇θ

    return nothing
end

"""
    source!(::HBModel)
    calculates the source term contribution to state variables
    this computation is done pointwise at each nodal point

    arguments:
    m -> model in this case HBModel
    F -> array of fluxes for each state variable
    Q -> array of state variables
    A -> array of aux variables
    t -> time, not used

    computations
    ∂ᵗu = -f×u
    ∂ᵗη = w|(z=0)
"""
@inline function source!(
    m::HBModel,
    S::Vars,
    Q::Vars,
    D::Vars,
    A::Vars,
    t::Real,
    direction,
)
    @inbounds begin
        u, v = Q.u # Horizontal components of velocity
        wz0 = A.wz0

        # f × u
        f = coriolis_force(m, A.y)
        S.u -= @SVector [-f * v, f * u]

        S.η += wz0
    end

    return nothing
end

"""
    coriolis_force(::HBModel)

northern hemisphere coriolis

# Arguments
- `m`: model object to dispatch on and get coriolis parameters
- `y`: y-coordinate in the box
"""
@inline coriolis_force(m::HBModel, y) = m.fₒ + m.β * y

"""
    wavespeed(::HBModel)

calculates the wavespeed for rusanov flux
"""
@inline wavespeed(m::HBModel, n⁻, _...) = abs(SVector(m.cʰ, m.cʰ, m.cᶻ)' * n⁻)

"""
    update_penalty(::HBModel)
    set Δη = 0 when computing numerical fluxes
"""
# We want not have jump penalties on η (since not a flux variable)
function update_penalty!(
    ::RusanovNumericalFlux,
    ::HBModel,
    n⁻,
    λ,
    ΔQ::Vars,
    Q⁻,
    A⁻,
    Q⁺,
    A⁺,
    t,
)
    ΔQ.η = -0

    return nothing
end

"""
    update_auxiliary_state!(::HBModel)

    applies the vertical filter to the zonal and meridional velocities to preserve numerical incompressibility
    applies an exponential filter to θ to anti-alias the non-linear advective term

    doesn't actually touch the aux variables any more, but we need a better filter interface than this anyways
"""
function update_auxiliary_state!(
    dg::DGModel,
    m::HBModel,
    Q::MPIStateArray,
    t::Real,
    elems::UnitRange,
)
    FT = eltype(Q)
    MD = dg.modeldata

    # `update_aux!` gets called twice, once for the real elements and once for
    # the ghost elements.  Only apply the filters to the real elems.
    if elems == dg.grid.topology.realelems
        # required to ensure that after integration velocity field is divergence free
        vert_filter = MD.vert_filter
        apply!(Q, (:u,), dg.grid, vert_filter, direction = VerticalDirection())

        exp_filter = MD.exp_filter
        apply!(Q, (:θ,), dg.grid, exp_filter, direction = VerticalDirection())
    end

    return true
end

"""
    update_auxiliary_state_gradient!(::HBModel)

    ∇hu to w for integration
    performs integration for w and pkin (should be moved to its own integral kernels)
    copies down w and wz0 because we don't have 2D structures

    now for actual update aux stuff
    implements convective adjustment by bumping the vertical diffusivity up by a factor of 1000 if dθdz < 0
"""
function update_auxiliary_state_gradient!(
    dg::DGModel,
    m::HBModel,
    Q::MPIStateArray,
    t::Real,
    elems::UnitRange,
)
    FT = eltype(Q)
    A = dg.state_auxiliary

    # store ∇ʰu as integrand for w
    function f!(m::HBModel, Q, A, D, t)
        @inbounds begin
            ν = viscosity_tensor(m)
            ∇u = ν \ D.ν∇u # minus sign included in gradient flux
            A.w = (∇u[1, 1] + ∇u[2, 2])
        end

        return nothing
    end
    nodal_update_auxiliary_state!(f!, dg, m, Q, t, elems; diffusive = true)

    # compute integrals for w and pkin
    indefinite_stack_integral!(dg, m, Q, A, t, elems) # bottom -> top
    reverse_indefinite_stack_integral!(dg, m, Q, A, t, elems) # top -> bottom

    # project w(z=0) down the stack
    # [1] to convert from range to integer
    # copy_stack_field_down! doesn't like ranges
    # eventually replace this with a reshape and broadcast
    index_w = varsindex(vars_state_auxiliary(m, FT), :w)[1]
    index_wz0 = varsindex(vars_state_auxiliary(m, FT), :wz0)[1]
    copy_stack_field_down!(dg, m, A, index_w, index_wz0, elems)

    return true
end

include("SimpleBoxProblem.jl")
include("LinearHBModel.jl")
include("BoundaryConditions.jl")
include("Courant.jl")

end

max_stars_repo_name: leios/CLIMA

max_issues_repo_name: leios/CLIMA

max_forks_repo_name: leios/CLIMA


text:
## dummy model to pack to WeibullPH models together

struct DoubleWeibullData <: AbstractData
    x::WeibullPHData
    y::WeibullPHData
end

nevents(data::DoubleWeibullData) = sum(data.x.status) + sum(data.y.status)

DoubleWeibullData() = DoubleWeibullData(WeibullPHData(), WeibullPHData())

mutable struct DoubleWeibullParam <: AbstractParam
    x::WeibullPHParam
    y::WeibullPHParam
end

struct DoubleWeibullFixedParam <: AbstractFixedParam
    x::WeibullPHFixedParam
    y::WeibullPHFixedParam
end

struct DoubleWeibullPrior <: AbstractPrior
    x::WeibullPHPrior
    y::WeibullPHPrior
end

DoubleWeibullPrior() = DoubleWeibullPrior(WeibullPHPrior(), WeibullPHPrior())

function DoubleWeibullModel()
    m1 = WeibullPHModel()
    m2 = WeibullPHModel()
    Model(DoubleWeibullParam(m1.param, m2.param), DoubleWeibullFixedParam(m1.fparam, m2.fparam), generate_data)
end

DoubleWeibullModel(param, fparam, gen_data=generate_data) = Model(param, fparam, gen_data)

struct DoubleWeibullFit <: AbstractFit
    x::WeibullPHFit
    y::WeibullPHFit
end

slice(i::Int, x::DoubleWeibullFit) = Param(slice(x.x), slice(x.y))

function sample_prior(B::Int, prior::DoubleWeibullPrior, fpar)
    par1 = sample_prior(B, prior.x, fpar.x)
    par2 = sample_prior(B, prior.y, fpar.y)

    par = Vector{Param}(undef, B)
    for b in 1:B
        par[b] = Param(par1[b], par2[b])
    end

    par
end

## return prior distribution (before sampling) or posterior sample (after sampling)
## of treatment effect parameter
function get_outcome(x::DoubleWeibullFit)
    if size(x.beta, 1) == 0
        [x.prior.x.beta, x.prior.y.beta]
    else
        hcat(x.x.beta, x.y.beta)
    end
end

generate_data(par, fpar) = predict_data(par, fpar, DoubleWeibullData(), 0.0)

function predict_data(par::DoubleWeibullParam, fpar::DoubleWeibullFixedParam, data0::DoubleWeibullData, caltime::Float64, nevents=NaN)
    DoubleWeibullData(predict_data(par.x, fpar.x, data0.x, caltime),
                        predict_data(par.y, fpar.y, data0.y, caltime))
end

function select_data(t::Float64, data::DoubleWeibullData; entry_stop=t)
    DoubleWeibullData(select_data(t, data.x; entry_stop=entry_stop),
         select_data(t, data.y; entry_stop=entry_stop))
end

function laplace_approx(data::DoubleWeibullData, param::DoubleWeibullParam, prior::DoubleWeibullPrior)
    vcat(laplace_approx(data.x, param.x, prior.x),
         laplace_approx(data.y, param.y, prior.y))
end

function get_approx(param::DoubleWeibullParam, new_data::DoubleWeibullData, prior::DoubleWeibullPrior)
    laplace_approx(new_data, param, prior)
end

function fit(M::Int, data::DoubleWeibullData; prior::DoubleWeibullPrior=DoubleWeibullPrior(), warmup::Int=Int(floor(M/2)))
    DoubleWeibullFit(fit(M, data.x; prior=prior.x, warmup=warmup), fit(M, data.y, prior=prior.y, warmup=warmup))
end

summary(x::DoubleWeibullFit) = summary(x.x), summary(x.y)

function doubleweibull_test_laplace_approx()
    m = DoubleWeibullModel()
    data = m.generate_data(m.param, m.fparam)

    la = laplace_approx(data, m.param, DoubleWeibullPrior())
    x = fit(10000, data; warmup=5000)

    print(summary(x))
    x, la
end


text:
"""
    TypeH()

Lorem ipsum dolor sit amet.
"""
struct TypeH

end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using Test

Test.@testset "io" begin
    include("test_hdf5.jl")
    include("test_bat_io.jl")
end


text:
using DungAnalyse
using Test


text:
function bug1_w1(w1, plt1)
    dt = 0
    start = (x = 0, y = 0)
    goal = (x = 10, y = 10)
    x_start = [0];
    y_start = [0];
    x_goal = [10];
    y_goal = [10];
    xx = [];yy = [];

    points = []
    jj = 0
    while (true)
        step = start.x + 0.01
        X, Y, obs = start_toward_goal(w1, start, goal, step)
        plot!(plt1, X, Y, legend=false, color=:black, linewidth=4)
        if (X[end] == goal.x && Y[end] == goal.y)
            println("Completed")
            println("Total distance is :", dt)

            return
        end

        path = make_path(w1, obs, X[end], Y[end])
        xx, yy, dt = follow_path(path, (x = x_goal[1], y = y_goal[1]))
        start = (x = xx, y = yy)

        plot!(plt1, [xx], [yy], legend=false, marker=([:dot :d]))

        jj += 1;
    end
    # display(plt1)
end

function get_commun(w1, obs, pointx, pointy)
    first = []
    second = []
    comon_edge = ()
    mid_way_dege = ()
    for arr in w1
        for ii = 1:size(arr, 1)
            first = arr[ii]
            if ii == size(arr, 1)
                second = arr[1]
            else
                second = arr[ii + 1]
            end
            if first.x == pointx && first.y == pointy && arr != obs
                comon_edge = first
                break;
            end
            if second.x == pointx && second.y == pointy && arr != obs
                comon_edge = second
                break;
            end
            if first.x < second.x && pointx in collect(first.x:0.01:second.x) && first.y == pointy
                mid_way_dege = (x = pointx, y = pointy)
                break;
            end
            if first.y < second.y && pointy in collect(first.y:0.01:second.y) && first.x == pointx
                mid_way_dege = (x = pointx, y = pointy)
                break;
            end

            if first.x > second.x && pointx in collect(second.x:0.01:first.x) && first.y == pointy
                mid_way_dege = (x = pointx, y = pointy)
                break;
            end
            if first.y > second.y && pointy in collect(second.y:0.01:first.y) && first.x == pointx
                mid_way_dege = (x = pointx, y = pointy)
                break;
            end
        end
        if !isempty(comon_edge) && !isempty(mid_way_dege)
            return comon_edge, mid_way_dege, arr
        end
    end

    return (), (), ()
end

function make_array(cords::Array)
    x = [];
    y = [];
    for tup in cords
        push!(x, tup.x)
        push!(y, tup.y)
    end
    push!(x, getindex(cords, 1).x)
    push!(y, getindex(cords, 1).y)
    return x, y
end

function path_eqn(start_point::NamedTuple, end_point::NamedTuple)
    slope = (end_point.y - start_point.y) / (end_point.x - start_point.x);

    y_cept = end_point.y - slope * end_point.x;
    return slope, y_cept;
end

function check_obstacle(edges::Array, postx, posty)
    slope = 0;y_cept = 0
    edge1 = ()
    edge2 = ()
    for i = 1:size(edges, 1)
        if i == size(edges, 1)
            edge1 = edges[i]
            edge2 = edges[1]
            slope, y_cept = path_eqn(edge1, edge2)
        else
            edge1 = edges[i]
            edge2 = edges[i + 1]
            slope, y_cept = path_eqn(edge1, edge2)
        end
        if slope == Inf
            eq_val = posty
        end
        if y_cept == Inf || slope == 0
            eq_val = edge1.y

        end
        if slope != Inf || slope != 0
            eqn = Polynomial([y_cept,slope]);
            eq_val = round(eqn(postx), digits=2)
        end
        if validate_point(edge1, edge2, (x = postx, y = posty))
            return true
        end
    end
    return false
end

function validate_point(edge1, edge2, point)
    order = []
    if (edge1.x == edge2.x && edge2.x == point.x)
        if (edge1.y < edge2.y && point.y in collect(edge1.y:0.01:edge2.y) )
            return true
        end
        if (edge1.y > edge2.y && point.y in collect(edge2.y:0.01:edge1.y) )
            return true
        end
    end
    if (edge1.y == edge2.y && edge2.y == point.y)
        if (edge1.x < edge2.x && point.x in collect(edge1.x:0.01:edge2.x) )
            return true
        end
        if (edge1.x > edge2.x && point.x in collect(edge2.x:0.01:edge1.x) )
            return true
        end
    end

    return false
end

function start_toward_goal(w1, start, goal, step)
    obs = false
    slope, y_cep = path_eqn(start, goal)
    eqn = Polynomial([y_cep,slope])
    plot!(eqn, step, goal.x)
    X = [];
    Y = [];
    while (step <= goal.x)
        eq_val = round(eqn(step), digits=2)
        push!(X, step);
        push!(Y, eq_val)
        for arr in w1
            obs = check_obstacle(arr, step, eq_val)
            if (obs)
                return X, Y, arr;
            end
        end
        step = round(step + 0.01, digits=2)
    end
    return [start.x, goal.x], [start.y, goal.y], []
end

function make_path(w1, obs, obx, oby)
    path = []
    for i = 1:size(obs, 1)
        comon_edge, mid_way_dege, array = get_commun(w1, obs, obs[i].x, obs[i].y)
        if (isempty(comon_edge) && isempty(mid_way_dege))
            push!(path, obs[i])
        elseif (!isempty(comon_edge) && !isempty(mid_way_dege))
            path = [(x = obx, y = oby),(x = 3, y = 4), (x = 4, y = 4), (x = 4, y = 12),
            (x = 12, y = 12), (x = 12, y = 6),(x = 6, y = 6),
            (x = 6, y = 5),(x = 12, y = 5),(x = 13, y = 5), (x = 13, y = 13),
            (x = 12, y = 13),(x = 3, y = 13),(x = 3, y = 12),(x = obx, y = oby)]
        end

    end
    return path
end

function follow_path(path, goal)
    X = [];
    Y = [];
    dist_total = 0
    start = path[1];
    pointx = typemax(Int32)
    pointy = typemax(Int32)
    rpointx = typemax(Int32)
    rpointy = typemax(Int32)
    shortest_dis = typemax(Int32)
    old_shotest_dis = 0
    next_point = (x = typemax(Int32), y = typemax(Int32))
    Xp = []
    Yp = []
    retx = []
    rety = []
    for i = 1:size(path, 1)

        if (i + 1 ) > size(path, 1)
            i = 0
        end
        next_point = path[i + 1];
        if start.x != next_point.x
            if start.x < next_point.x
                X = collect(LinRange(start.x, next_point.x, 100))
            else
                X = collect(LinRange(next_point.x, start.x, 100))
            end
            if start.y == next_point.y
                Y = next_point.y * ones(size(X, 1))
            end
        end
        if start.y != next_point.y
            if start.y < next_point.y
                Y = collect(LinRange(start.y, next_point.y, 100))
            else
                Y = collect(LinRange(next_point.y, start.y, 100))
            end
            if start.x == next_point.x
                X = next_point.x * ones(size(Y, 1))
            end
        end
        push!(Xp, start.x)
        push!(Yp, start.y)
        push!(retx, X)
        push!(rety, Y)

        dis, pointx, pointy = shortest_distance(goal, (x = X, y = Y))
        if dis <= shortest_dis
            shortest_dis = dis
            rpointx = round(pointx, digits=1)
            rpointy = round(pointy, digits=1)
        end
        dist_total += round(sqrt((start.x - next_point.x)^2 + (start.y - next_point.y)^2))
        start = next_point

    end
    push!(Xp, Xp[1]), push!(Yp, Yp[1])
    plot!(Xp, Yp, legend=false, color=:black, linewidth=4)
    return rpointx, rpointy, dist_total
end

# bug2()

function shortest_distance(goal, points)
    shortest_dis = typemax(Int32);
    pointx = typemax(Int32)
    pointy = typemax(Int32)
    dis = 0;
    for i = 1:size(points.x, 1)
        dis =  round(sqrt((goal.x - points.x[i])^2 + (goal.y - points.y[i])^2), digits=2)
        if dis < shortest_dis
            shortest_dis = dis;
            pointx = round(points.x[i], digits=1)
            pointy = round(points.y[i], digits=1)

        end
    end
    return shortest_dis, pointx, pointy;
end

function bug1_w2(w2)
    goal = (x = 35, y = 0)
    start = (x = 0, y = 0)
    path = [(x = 30, y = 0), (x = 29, y = 0), (x = 29, y = 5), (x = 20, y = 5),
    (x = 20, y = 0), (x = 19, y = 0), (x = 19, y = 5), (x = 10, y = 5),
    (x = 10, y = 0), (x = 9, y = 0), (x = 9, y = 5), (x = -5, y = 5),
    (x = -5, y = -5), (x = 4, y = -5), (x = 4, y = 1), (x = 5, y = 1),
    (x = 5, y = -5), (x = 14, y = -5), (x = 14, y = 1), (x = 15, y = 1),
     (x = 15, y = -5), (x = 24, y = -5), (x = 24, y = 1), (x = 25, y = 1),
     (x = 25, y = -6),(x = -6, y = -6), (x = -6, y = 6), (x = 30, y = 6) ,
     (x = 30, y = 0)]

    dist = follower(path, goal)
     #  (x=30, y=0), (x=35, y=0)]
    println("Total Distance is : ", dist)

    plot!([30, goal.x], [ goal.y, goal.y], legend=false, color=:black, arrow=:arrow, linewidth=4)
    plot!([start.x, 4], [ start.y, 0], legend=false, color=:black, arrow=:arrow, linewidth=4)
end

function follower(path, goal)

    start = path[1];
    next_point = ()
    dist_total = 0
    for i = 1:size(path, 1) - 1

        if (i + 1 ) > size(path, 1)
            i = 0
        end
        next_point = path[i + 1];
        # println(start, next_point)
        plot!([start.x,  next_point.x], [  start.y, next_point.y], legend=false, arrow=:arrow, color=:black, linewidth=4)
        dist_total += round(sqrt((start.x - next_point.x)^2 + (start.y - next_point.y)^2))
        start = next_point;
    end
    # plot!([ goal.x, path[1].x], [ goal.y, path[1].y], legend=false, color=:blue, linewidth=1)
    return dist_total
end
# bug1()


text:
module SimpleNetworks

using DataFrames: DataFrameRow

import StaticTrafficAssignment:
        AbstractNetwork, AbstractLink, AbstractZone, AbstractBush,
        CostFunction, dijkstra, topologicalorder, orderednodes,
        has_link, add_link!, rem_link!, add_zone!, rem_zone!,
        idx, numnodes, numzones, numlinks, numsources, numsinks,
        sources, sinks, links, zones,
        outneighbors, inneighbors, outdegree, indegree, fstar, bstar,
        upn, dwn, props,
        id, issource, issink, throughflowallowed, src, net

export SimpleNetwork, SimpleLink, SimpleZone, SimpleBush

include("simplelink.jl")
include("simplezone.jl")
include("simplenetwork.jl")
include("simplebush.jl")

end  # module SimpleNetworks