text:
"""
`SolutionSDE`: Solution of a stochastic differential equation

Contains all fields necessary to store the solution of an SDE.

### Fields
* `conv`: type of the solution: "strong" or "weak"
* `nd`: dimension of the dynamical variable ``q``
* `nm`: dimension of the Wiener process
* `nt`: number of time steps to store
* `ns`: number of sample paths
* `ni`: number of initial conditions
* `t`:  time steps
* `q`:  solution `q[nd, nt+1, ns, ni]` with `q[:,0,:,:]` the initial conditions
* `W`:  Wiener process driving the stochastic process q
* `K`:  integer parameter defining the truncation of the increments of the Wiener process (for strong solutions),
*       A = √(2 K Δt |log Δt|) due to Milstein & Tretyakov; if K=0 no truncation
* `ntime`: number of time steps to compute
* `nsave`: save every nsave'th time step

"""
mutable struct SolutionSDE{dType, tType, NQ, NW} <: StochasticSolution{dType, tType, NQ, NW}
    conv::String
    nd::Int
    nm::Int
    nt::Int
    ns::Int
    ni::Int
    t::TimeSeries{tType}
    q::SStochasticDataSeries{dType,NQ}
    W::WienerProcess{dType,tType,NW}
    K::Int
    ntime::Int
    nsave::Int
    counter::Int
end

function SolutionSDE(equation::SDE{DT,TT,VT,BT}, Δt::TT, ntime::Int, nsave::Int=1; K::Int=0, conv::String="strong") where {DT,TT,VT,BT}
    nd = equation.d
    nm = equation.m
    ns = equation.ns
    ni = equation.n
    nt = div(ntime, nsave)

    if nd==ns==ni==1
        NQ = 1
    elseif ns==ni==1
        NQ = 2
    elseif ns==1 || ni==1
        NQ = 3
    else
        NQ = 4
    end

    @assert conv=="strong" || (conv=="weak" && K==0)
    @assert DT <: Number
    @assert TT <: Real
    @assert nd > 0
    @assert ns > 0
    @assert ni > 0
    @assert nsave > 0
    @assert ntime == 0 || ntime ≥ nsave
    @assert mod(ntime, nsave) == 0

    t = TimeSeries{TT}(nt, Δt, nsave)

    q = SStochasticDataSeries(DT, nd, nt, ns, ni)

    # Holds the Wiener process data for ALL computed time steps
    # Wiener process increments are automatically generated here
    W = WienerProcess(DT, nm, ntime, ns, Δt, conv)
    NW = ndims(W.ΔW)

    s = SolutionSDE{DT,TT,NQ,NW}(conv, nd, nm, nt, ns, ni, t, q, W, K, ntime, nsave, 0)
    set_initial_conditions!(s, equation)
    return s
end

function SolutionSDE(equation::SDE{DT,TT,VT,BT}, Δt::TT, dW::Array{DT, NW}, dZ::Array{DT, NW}, ntime::Int, nsave::Int=1; K::Int=0, conv::String="strong") where {DT,TT,VT,BT,NW}
    nd = equation.d
    nm = equation.m
    ns = equation.ns
    ni = equation.n
    nt = div(ntime, nsave)

    if nd==ns==ni==1
        NQ = 1
    elseif ns==ni==1
        NQ = 2
    elseif ns==1 || ni==1
        NQ = 3
    else
        NQ = 4
    end

    @assert size(dW) == size(dZ)

    if NW==1
        @assert ns==nm==1
        @assert ntime==length(dW)
    elseif NW==2
        @assert nm==size(dW,1)
        @assert ntime==size(dW,2)
        @assert ns==1
    elseif NW==3
        @assert nm==size(dW,1)
        @assert ntime==size(dW,2)
        @assert ns==size(dW,3)
    end

    @assert conv=="strong" || (conv=="weak" && K==0)
    @assert DT <: Number
    @assert TT <: Real
    @assert nd > 0
    @assert ns > 0
    @assert ni > 0
    @assert nsave > 0
    @assert ntime == 0 || ntime ≥ nsave
    @assert mod(ntime, nsave) == 0

    t = TimeSeries{TT}(nt, Δt, nsave)

    q = SStochasticDataSeries(DT, nd, nt, ns, ni)

    # Holds the Wiener process data for ALL computed time steps
    # Wiener process increments are prescribed by the arrays ΔW and ΔZ
    W = WienerProcess(Δt, dW, dZ, conv)

    s = SolutionSDE{DT,TT,NQ,NW}(conv, nd, nm, nt, ns, ni, t, q, W, K, ntime, nsave, 0)
    set_initial_conditions!(s, equation)
    return s
end

function SolutionSDE(t::TimeSeries{TT}, q::SStochasticDataSeries{DT,NQ}, W::WienerProcess{DT,TT,NW}; K::Int=0) where {DT,TT,NQ,NW}
    # extract parameters
    conv = W.conv
    nd = q.nd
    ns = q.ns
    ni = q.ni
    nt = t.n
    nm = W.nd
    ntime = W.nt
    nsave = t.step

    @assert conv=="strong" || (conv=="weak" && K==0)
    @assert ntime==nt*nsave
    @assert q.nt == nt
    @assert W.ns == q.ns

    # create solution
    SolutionSDE{DT,TT,NQ,NW}(conv, nd, nm, nt, ns, ni, t, q, W, K, ntime, nsave, 0)
end

# If the Wiener process W data are not available, creates a one-element zero array instead
# For instance used when reading a file with no Wiener process data saved
function SolutionSDE(t::TimeSeries{TT}, q::SStochasticDataSeries{DT,NQ}; K::Int=0, conv::String="strong") where {DT,TT,NQ}
    # extract parameters
    nd = q.nd
    ns = q.ns
    ni = q.ni
    nt = t.n
    nsave = t.step
    ntime = nt*nsave

    @assert conv=="strong" || (conv=="weak" && K==0)
    @assert q.nt == nt

    W = WienerProcess(t.Δt, [0.0], [0.0], conv)

    nm = W.nd
    NW = ndims(W.ΔW)

    # create solution
    SolutionSDE{DT,TT,NQ,NW}(conv, nd, nm, nt, ns, ni, t, q, W, K, ntime, nsave, 0)
end

function SolutionSDE(file::String)
    # open HDF5 file
    get_config(:verbosity) > 1 ? @info("Reading HDF5 file ", file) : nothing
    h5 = h5open(file, "r")

    # read attributes
    nsave = read(attrs(h5)["nsave"])
    ni    = read(attrs(h5)["ni"])

    # reading data arrays
    t = TimeSeries(read(h5["t"]), nsave)

    if exists(attrs(h5),"conv")
        conv = read(attrs(h5)["conv"])
    else
        conv = "strong"
    end

    W_exists = exists(h5, "ΔW") && exists(h5, "ΔZ")

    if W_exists == true
        W = WienerProcess(t.Δt, read(h5["ΔW"]), read(h5["ΔZ"]), conv)
    end

    if exists(attrs(h5),"K")
        K = read(attrs(h5)["K"])
    else
        K=0
    end

    q_array = read(h5["q"])

    close(h5)

    if ndims(q_array)==3 && ni>1
        q = SStochasticDataSeries(q_array,IC=true)
    else
        q = SStochasticDataSeries(q_array)
    end

    # create solution
    if W_exists == true
        return SolutionSDE(t, q, W, K=K)
    else
        return SolutionSDE(t, q, K=K, conv=conv)
    end

end

time(sol::SolutionSDE)  = sol.t.t
ntime(sol::SolutionSDE) = sol.ntime
nsave(sol::SolutionSDE) = sol.nsave

function set_initial_conditions!(sol::SolutionSDE{DT,TT}, equ::SDE{DT,TT}) where {DT,TT}
    set_initial_conditions!(sol, equ.t₀, equ.q₀)
end

function set_initial_conditions!(sol::SolutionSDE{DT,TT}, t₀::TT, q₀::Union{Array{DT}, Array{TwicePrecision{DT}}}) where {DT,TT}
    # Sets the initial conditions sol.q[0] with the data from q₀
    # q₀ may be 1D (nd elements - single deterministic initial condition),
    # 2D (nd x ns or nd x ni matrix - single random or multiple deterministic initial condition),
    # or 3D (nd x ns x ni matrix - multiple random initial condition)
    set_data!(sol.q, q₀, 0)
    compute_timeseries!(sol.t, t₀)
end

# copies the m-th initial condition for the k-th sample path from sol.q to q
function get_initial_conditions!(sol::SolutionSDE{DT,TT}, q::Union{Vector{DT}, Vector{TwicePrecision{DT}}}, k, m) where {DT,TT}

    @assert k ≤ sol.ns
    @assert m ≤ sol.ni

    N = ndims(sol.q)

    if N==1
        # 1D space, 1 sample path and 1 initial condition, k==m==1
        q[1] = get_data!(sol.q, 0)
    elseif N==2
        # Multidimensional space, 1 sample path and 1 initial condition, k==m==1
        get_data!(sol.q, q, 0)
    elseif N==3
        # Multidimensional space, with either 1 sample path or 1 initial condition
        if sol.ns==1
            # 1 sample path, k==1, reading the m-th initial condition
            get_data!(sol.q, q, 0, m)
        else
            #1 initial condition, m==1, reading the k-th sample path
            get_data!(sol.q, q, 0, k)
        end
    elseif N==4
        # Multidimensional space, multiple sample paths and initial conditions
        get_data!(sol.q, q, 0, k, m)
    end

end

function CommonFunctions.set_solution!(sol::SolutionSDE{DT,TT,NQ,NW}, q::Union{Vector{DT}, Vector{TwicePrecision{DT}}}, n, k, m) where {DT,TT,NQ,NW}

    if mod(n, sol.nsave) == 0

        if NQ ∈ (1,2)
            # Single sample path and a single initial condition, k==m==1
            @assert k==m==1
            set_data!(sol.q, q, div(n, sol.nsave))
        elseif NQ==3
            if sol.ni==1
                #Single initial condition, multiple sample paths, m==1
                @assert m==1
                set_data!(sol.q, q, div(n, sol.nsave), k)
            else
                #Single sample path, multiple initial conditions, k==1
                @assert k==1
                set_data!(sol.q, q, div(n, sol.nsave), m)
            end
        elseif NQ==4
            # Multiple sample paths and initial conditions
            set_data!(sol.q, q, div(n, sol.nsave), k, m)
        end

        sol.counter += 1
    end
end

function reset!(sol::SolutionSDE)
    reset!(sol.q)
    compute_timeseries!(sol.t, sol.t[end])
    generate_wienerprocess!(sol.W)
    sol.counter = 0
end

"""
Creates HDF5 file and initialises datasets for SDE solution object.
  It is implemented as one function for all NQ and NW cases, rather than several
  separate cases as was done for SolutionODE.
  nt - the total number of time steps to store
  ntime - the total number of timesteps to be computed
"""
function create_hdf5(solution::SolutionSDE{DT,TT,NQ,NW}, file::AbstractString, nt::Int=solution.nt, ntime::Int=solution.ntime; save_W=true) where {DT,TT,NQ,NW}
    @assert nt ≥ 1
    @assert ntime ≥ 1

    # create HDF5 file
    h5 = createHDF5(solution, file)

    # Adding the attributes specific to SolutionSDE that were not added above
    save_attributes(solution, h5)
    attrs(h5)["nt"] = nt
    attrs(h5)["ntime"] = ntime

    # create dataset
    # nt and ntime can be used to set the expected total number of timesteps to be saved,
    # so that the size of the array does not need to be adapted dynamically.
    # Right now, it has to be set as dynamical size adaptation is not yet
    # working. The default value is the size of the solution structure.
    if NQ==1
        q = d_create(h5, "q", datatype(DT), dataspace((nt+1,)), "chunk", (1,))
        # copy initial conditions
        q[1] = solution.q[0]
    elseif NQ==2
        q = d_create(h5, "q", datatype(DT), dataspace(solution.nd, nt+1), "chunk", (solution.nd,1))
        # copy initial conditions
        q[1:solution.nd, 1] = solution.q[1:solution.nd, 0]
    elseif NQ==3
        if solution.ns>1
            q = d_create(h5, "q", datatype(DT), dataspace(solution.nd, nt+1, solution.ns), "chunk", (solution.nd,1,1))
        else
            q = d_create(h5, "q", datatype(DT), dataspace(solution.nd, nt+1, solution.ni), "chunk", (solution.nd,1,1))
        end
        # copy initial conditions
        q[:, 1, :] = solution.q.d[:, 1, :]
    else
        q = d_create(h5, "q", datatype(DT), dataspace(solution.nd, nt+1, solution.ns, solution.ni), "chunk", (solution.nd,1,1,1))
        # copy initial conditions
        q[:, 1, :, :] = solution.q.d[:, 1, :, :]
    end

    if save_W==true
        # creating datasets to store the Wiener process increments
        if NW==1
            dW = d_create(h5, "ΔW", datatype(DT), dataspace((ntime,)), "chunk", (1,))
            dZ = d_create(h5, "ΔZ", datatype(DT), dataspace((ntime,)), "chunk", (1,))
        elseif NW==2
            dW = d_create(h5, "ΔW", datatype(DT), dataspace(solution.nm, ntime), "chunk", (solution.nm,1))
            dZ = d_create(h5, "ΔZ", datatype(DT), dataspace(solution.nm, ntime), "chunk", (solution.nm,1))
        elseif NW==3
            dW = d_create(h5, "ΔW", datatype(DT), dataspace(solution.nm, ntime, solution.ns), "chunk", (solution.nm,1,1))
            dZ = d_create(h5, "ΔZ", datatype(DT), dataspace(solution.nm, ntime, solution.ns), "chunk", (solution.nm,1,1))
        end
    end

    # Creating a dataset for storing the time series
    t = d_create(h5, "t", datatype(TT), dataspace((nt+1,)), "chunk", (1,))
    t[1] = solution.t[0]

    return h5
end

"""
Append solution to HDF5 file.
  offset - start writing q at the position offset+2
  offset2- start writing ΔW, ΔZ at the position offset2+1
"""
function CommonFunctions.write_to_hdf5(solution::SolutionSDE{DT,TT,NQ,NW}, h5::HDF5File, offset=0, offset2=offset) where {DT,TT,NQ,NW}
    # set convenience variables and compute ranges
    d   = solution.nd
    m   = solution.nm
    n   = solution.nt
    s   = solution.ns
    i   = solution.ni
    ntime = solution.ntime
    j1  = offset+2
    j2  = offset+1+n
    jw1 = offset2+1
    jw2 = offset2+ntime

    # # extend dataset if necessary
    # if size(x, 2) < j2
    #     set_dims!(x, (d, j2))
    # end

    # saving the time time series
    h5["t"][j1:j2] = solution.t[1:n]

    # copy data from solution to HDF5 dataset
    if NQ==1
        h5["q"][j1:j2] = solution.q[1:n]
    elseif NQ==2
        h5["q"][:, j1:j2] = solution.q[:, 1:n]
    elseif NQ==3
        h5["q"][:, j1:j2, :] = solution.q[:, 1:n, :]
    else
        h5["q"][:, j1:j2, :, :] = solution.q[:, 1:n, :, :]
    end

    if exists(h5, "ΔW") && exists(h5, "ΔZ")
        # copy the Wiener process increments from solution to HDF5 dataset
        if NW==1
            h5["ΔW"][jw1:jw2] = solution.W.ΔW.d[1:ntime]
            h5["ΔZ"][jw1:jw2] = solution.W.ΔZ.d[1:ntime]
        elseif NW==2
            h5["ΔW"][:,jw1:jw2] = solution.W.ΔW.d[:,1:ntime]
            h5["ΔZ"][:,jw1:jw2] = solution.W.ΔZ.d[:,1:ntime]
        elseif NW==3
            h5["ΔW"][:,jw1:jw2,:] = solution.W.ΔW.d[:,1:ntime,:]
            h5["ΔZ"][:,jw1:jw2,:] = solution.W.ΔZ.d[:,1:ntime,:]
        end
    end

    return nothing
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    BlockWeighting(sides)
    BlockWeighting(side₁, side₂, ...)

A weighting method that assigns weights to points in spatial object
based on blocks of given `sides`. The number `n` of points inside a
block determines the weights `1/n` of these points.
"""
struct BlockWeighting{T,N} <: WeightingMethod
  sides::SVector{N,T}
end

BlockWeighting(sides::NTuple{N,T}) where {N,T} =
  BlockWeighting{T,N}(sides)

BlockWeighting(sides::Vararg) = BlockWeighting(sides)

weight(object, method::BlockWeighting) =
  weight(domain(object), method)

function weight(object::AbstractDomain, method::BlockWeighting)
  p = partition(object, BlockPartition(method.sides))

  weights = Vector{Float64}(undef, nelms(object))
  for s in subsets(p)
    n = length(s)
    weights[s] .= 1/n
  end

  SpatialWeights(object, weights)
end

ext: jl


text:
@testset "ranking" begin
    e = Set([1, 2, 4])
    q = Set([1])
    S = [1 0 0.5; 0 1 0.5; 0.5 0.5 1]
    rl = [3, 2]
    rs = [3.0, 2]
    r = Ranking(e, q, S, rl, rs)

    @testset "constructors" begin
        @test isa(r, Ranking)
        @test r.e == e && r.q == q && r.S == S && r.rl == rl && r.rs == rs
    end
end


text:
##
function get_pdf_samples(samples, trans_samples, log_posterior, prop_dist, lb, ub, names)
    n = size(samples, 2)
    n == size(trans_samples,2) || throw(DimensionMismatch("samples dimensions don't match"))

    g = zeros(n)
    p = zeros(n)

    for i in 1:n
        g_sample = @views(trans_samples[:,i])
        p_sample = @views(samples[:,i])
        logdetJ = sum(logjacobian_contribution.(transform.(p_sample, lb, ub), lb, ub))
        g[i] = logpdf(prop_dist, g_sample)
        p[i] = isnothing(names) ? log_posterior(p_sample) + logdetJ : log_posterior(NamedTuple{Tuple(names)}(p_sample)) + logdetJ
    end

    return p, g, p-g
end

## Iterative algo
function iterative_algorithm(l₁, l₂, n₁, n₂; tol, maxiter)

    lstar = median(l₁)
    r = 0.0
    s₁ = n₁ / (n₁ + n₂)
    s₂ = n₂ / (n₁ + n₂)

    logml = -Inf
    ϵ = 0.0
    for i = 1:maxiter
        numterm = 0.0
        denomterm = 0.0
        for (l₁ⱼ, l₂ⱼ) in zip(l₁, l₂)
            el1 = clamp(exp(l₁ⱼ - lstar), 0.0, 1e100) # Avoiding numerical errors
            el2 = clamp(exp(l₂ⱼ - lstar), 1e-100, 1e100)
            numterm += el2 / (s₁ * el2 + s₂ * r)
            denomterm += 1 / (s₁ * el1 + s₂ * r)

        end
        rnew = (1/n₂) * numterm / (denomterm / n₁)

        logmlnew = log(rnew) + lstar
        ϵ = abs((logmlnew - logml) / logmlnew)
        if ϵ < tol
            return logmlnew, i
        end
        logml = logmlnew
        r = rnew
    end
    @warn "Maximum number of iterations ($maxiter) reached before convergence under the tolerance level $tol"
    return logml, maxiter
end


text:
#################### Diag

struct Diag{Tv <: AbstractFloat, Ti <: Integer} <: AbstractArray{Tv,1}
    index::Vector{Ti}
    val::Vector{Tv}
end

function spdiag(A::Matrix{Tv}) where {Tv}
    m, n = size(A)
    index = Vector{Int}(undef, min(m,n))
    z = 1
    for i = 1:min(m,n)
        index[i] = z
        z += m+1
    end
    val = reshape(A, length(A))
    Diag(index, val)
end

function spdiag(A::SparseArrays.SparseMatrixCSC{Tv,Ti}) where {Tv,Ti}
    m, n = size(A)
    index = zeros(Ti, min(m,n))
    val = reshape(A.nzval, length(A.nzval))
    for j = 1:min(m,n)
        for z = A.colptr[j]:A.colptr[j+1]-1
            i = A.rowval[z]
            if i == j
                index[i] = z
                break
            elseif i > j
                break
            end
        end
    end
    Diag(index, val)
end

function spdiag(A::SparseCSR{Tv,Ti}) where {Tv,Ti}
    m, n = size(A)
    index = zeros(Ti, min(m,n))
    val = reshape(A.val, length(A.val))
    for i = 1:min(m,n)
        for z = A.rowptr[i]:A.rowptr[i+1]-1
            j = A.colind[z]
            if i == j
                index[i] = z
                break
            elseif i < j
                break
            end
        end
    end
    Diag(index, val)
end

function spdiag(A::SparseCSC{Tv,Ti}) where {Tv,Ti}
    m, n = size(A)
    index = zeros(Ti, min(m,n))
    val = reshape(A.val, length(A.val))
    for j = 1:min(m,n)
        for z = A.colptr[j]:A.colptr[j+1]-1
            i = A.rowind[z]
            if i == j
                index[i] = z
                break
            elseif i > j
                break
            end
        end
    end
    Diag(index, val)
end

function spdiag(A::SparseCOO{Tv,Ti}) where {Tv,Ti}
    m, n = size(A)
    index = zeros(Ti, min(m,n))
    val = reshape(A.val, length(A.val))
    for z = 1:SparseArrays.nnz(A)
        i = A.rowind[z]
        j = A.colind[z]
        if i == j
            index[i] = z
        end
    end
    Diag(index, val)
end

function Base.length(A::Diag{Tv,Ti}) where {Tv,Ti}
    return length(A.index)
end

function Base.size(A::Diag{Tv,Ti}) where {Tv,Ti}
    return (length(A.index),)
end

function Base.getindex(A::Diag{Tv,Ti}, i::Ti) where {Tv,Ti}
    z = A.index[i]
    if z == 0
        return Tv(0)
    else
        return A.val[z]
    end
end

function Base.setindex!(A::Diag{Tv,Ti}, value::Tv, i::Ti) where {Tv,Ti}
    z = A.index[i]
    if z != 0
        A.val[z] = value
    end
end


text:
#!/usr/bin/env julia

using Pkg
haskey(Pkg.installed(), "Documenter") || Pkg.add("Documenter")

using Documenter
using PrettyPrinting

# Highlight indented code blocks as Julia code.
using Markdown
Markdown.Code(code) = Markdown.Code("julia", code)
Documenter.Utilities.Markdown2._convert_inline(s::Markdown.Code) =
    Documenter.Utilities.Markdown2.CodeSpan(s.code)

makedocs(
    sitename = "PrettyPrinting.jl",
    pages = [
        "Home" => "index.md",
    ],
    modules = [PrettyPrinting])

deploydocs(
    repo = "github.com/rbt-lang/PrettyPrinting.jl.git",
)


text:
prettify_expr(expr) = expr
prettify_expr(f::Function) = nameof(f)
prettify_expr(expr::Expr) = Expr(expr.head, prettify_expr.(expr.args)...)

@latexrecipe function f(eqs::Vector{Equation})
    # Set default option values.
    env --> :align
    cdot --> false

    # Convert both the left- and right-hand sides to expressions of basic types
    # that latexify can deal with

    rhs = getfield.(eqs, :rhs)
    rhs = prettify_expr.(_toexpr(rhs; canonicalize=false))
    rhs = [postwalk(x -> x isa Expr && length(x.args) == 1 ? x.args[1] : x, eq) for eq in rhs]
    rhs = [postwalk(x -> x isa Expr && x.args[1] == :_derivative && length(x.args[2].args) == 2 ? :($(Symbol(:d, x.args[2]))/($(Symbol(:d, x.args[2].args[2])))) : x, eq) for eq in rhs]
    rhs = [postwalk(x -> x isa Expr && x.args[1] == :_derivative ? "\\frac{d\\left($(Latexify.latexraw(x.args[2]))\\right)}{d$(Latexify.latexraw(x.args[3]))}" : x, eq) for eq in rhs]

    lhs = getfield.(eqs, :lhs)
    lhs = prettify_expr.(_toexpr(lhs; canonicalize=false))
    lhs = [postwalk(x -> x isa Expr && length(x.args) == 1 ? x.args[1] : x, eq) for eq in lhs]
    lhs = [postwalk(x -> x isa Expr && x.args[1] == :_derivative && length(x.args[2].args) == 2 ? :($(Symbol(:d, x.args[2]))/($(Symbol(:d, x.args[2].args[2])))) : x, eq) for eq in lhs]
    lhs = [postwalk(x -> x isa Expr && x.args[1] == :_derivative ? "\\frac{d\\left($(Latexify.latexraw(x.args[2]))\\right)}{d$(Latexify.latexraw(x.args[3]))}" : x, eq) for eq in lhs]

    return lhs, rhs
end

Base.show(io::IO, ::MIME"text/latex", x::Num) = print(io, latexify(x))
Base.show(io::IO, ::MIME"text/latex", x::Symbolic) = print(io, latexify(x))
Base.show(io::IO, ::MIME"text/latex", x::Vector{Equation}) = print(io, latexify(x))
Base.show(io::IO, ::MIME"text/latex", x::AbstractArray{Num}) = print(io, latexify(x))

# `_toexpr` is only used for latexify
function _toexpr(O; canonicalize=true)
    if canonicalize
        canonical, O = canonicalexpr(O)
        canonical && return O
    else
        !istree(O) && return O
    end

    op = operation(O)
    args = arguments(O)
    if op isa Differential
        ex = _toexpr(args[1]; canonicalize=canonicalize)
        wrt = _toexpr(op.x; canonicalize=canonicalize)
        return :(_derivative($ex, $wrt))
    elseif op isa Sym
        isempty(args) && return nameof(op)
        return Expr(:call, _toexpr(op; canonicalize=canonicalize), _toexpr(args; canonicalize=canonicalize)...)
    end
    return Expr(:call, op, _toexpr(args; canonicalize=canonicalize)...)
end
_toexpr(s::Sym; kw...) = nameof(s)

function canonicalexpr(O)
    !istree(O) && return true, O
    op = operation(O)
    args = arguments(O)
    if op === (^)
        if length(args) == 2 && args[2] isa Number && args[2] < 0
            ex = _toexpr(args[1])
            if args[2] == -1
                expr = Expr(:call, inv, ex)
            else
                expr = Expr(:call, ^, Expr(:call, inv, ex), -args[2])
            end
            return true, expr
        end
    end
    return false, O
end
for fun in [:_toexpr]
    @eval begin
        function $fun(eq::Equation; kw...)
            Expr(:(=), $fun(eq.lhs; kw...), $fun(eq.rhs; kw...))
        end

        $fun(eqs::AbstractArray; kw...) = map(eq->$fun(eq; kw...), eqs)
        $fun(x::Integer; kw...) = x
        $fun(x::AbstractFloat; kw...) = x
    end
end
_toexpr(x::Num; kw...) = _toexpr(value(x); kw...)


text:
using PiecewiseOrthogonalPolynomials, ClassicalOrthogonalPolynomials, BlockArrays, Test, FillArrays, LinearAlgebra

@testset "transform" begin
    for r in (range(-1, 1; length=2), range(-1, 1; length=4), range(0, 1; length=4)), T in (Chebyshev(), Legendre())
        P = PiecewisePolynomial(T, r)
        x = axes(P,1)
        @test P[:,Block.(Base.OneTo(3))] \ x ≈ (P\ x)[Block.(1:3)]
        @test (P/P\x)[0.1] ≈ 0.1
        @test (P/P\exp.(x))[0.1] ≈ exp(0.1)
    end

    r = range(-1, 1; length=10_000)
    P = PiecewisePolynomial(Chebyshev(), r); Pₙ = P[:,Block.(1:1000)]; x = axes(P,1)
    @time u = Pₙ / Pₙ \ cos.(10_000x.^2);
    @test u[0.1] ≈ cos(10_000*0.1^2)

    P = ContinuousPolynomial{0}(r)
    x = axes(P,1)
    @test P[:,Block.(Base.OneTo(3))] \ x ≈ (P\ x)[Block.(1:3)]
    @test (P/P\x)[0.1] ≈ 0.1
    @test (P/P\exp.(x))[0.1] ≈ exp(0.1)
end

@testset "lowering" begin
    c = [randn(20); zeros(∞)]
    for r in (range(-1, 1; length=2), range(0, 1; length=2), range(0, 1; length=4))
        P = ContinuousPolynomial{0}(r)
        C = ContinuousPolynomial{1}(r)
        R = P \ C
        for x in range(first(r), last(r); length=100)
            @test (C*c)[x] ≈ (P*(R*c))[x]
        end
    end
end

@testset "mass" begin
    for r in (range(-1,1; length=2), range(0,1; length=4), range(0, 1; length=4))
        P = ContinuousPolynomial{0}(r)
        C = ContinuousPolynomial{1}(r)
        JR = Block.(1:10)
        KR = Block.(1:11)
        R = P\C
        @test R[KR,JR]'*((P'P)[KR,KR]*R[KR,JR]) ≈ (C'C)[JR,JR]
        @test (P'C)[JR,JR] ≈ (C'P)[JR,JR]'

    end
end

@testset "derivative" begin
    r = range(0,1; length=4)
    C = ContinuousPolynomial{1}(r)
    P = ContinuousPolynomial{0}(r)
    x = axes(C,1)
    D = Derivative(x)
    P\D*C
end

@testset "multiplication" begin
    r = range(-1,1; length=4)
    P = ContinuousPolynomial{0}(r)
    C = ContinuousPolynomial{1}(r)
    x = axes(P,1)
    a = P / P \ broadcast(x -> abs(x) ≤ 1/3 ? 1.0 : 0.5, x)
    @test (P \ (a .* P))[Block.(1:2), Block.(1:2)] ≈ Diagonal([0.5,1,0.5,0.5,1,0.5])

    c = [randn(4); Zeros(∞)]
    @test ((a .* C) * c)[0.1] ≈ (C*c)[0.1]
    @test ((a .* C) * c)[0.4] ≈ (C*c)[0.4]/2
end

@testset "weak Laplacian" begin
    r = range(0,1; length=5)
    C = ContinuousPolynomial{1}(r)
    x = axes(C,1)
    D = Derivative(x)
    Δ = -(D*C)'*(D*C)
    M = C'C
    L = Δ + M
    KR = Block.(1:10)
    @time L[KR,KR]
end


text:
import ParticleSwarmOptimizer
const pso = ParticleSwarmOptimizer
using Test

@testset "ParticleSwarmOptimizer" begin
    include("TestFunctions.jl")
    include("helper.jl")
    include("neighbourhood.jl")
    include("pso.jl")
end


text:
# Placeholder entry point script file

using Logging

function main()
   @info "Inside main"
   println("Hello world")
end

main()

max_stars_repo_name: NREL-SIIP/SIIP-Tutorial

max_issues_repo_name: NREL-SIIP/SIIP-Tutorial

max_forks_repo_name: NREL-SIIP/SIIP-Tutorial


text:
using PyTensorFlow
using PyCall
using LinearAlgebra
using PyPlot
using Random
Random.seed!(233)
include("fwi_util.jl")
include("fwi_util_op.jl")
np = pyimport("numpy")

# argsparse.jl
# ENV["CUDA_VISIBLE_DEVICES"] = 1
# ENV["PARAMDIR"] = "Src/params/"
# config = tf.ConfigProto(device_count = Dict("GPU"=>0))

nz = 200
nx = 200
dz = 20
dx = 20
nSteps = 2001
dt = 0.0025
f0 = 4.5
filter_para = [0, 0.1, 100.0, 200.0]
nPml = 32
isAc = true
nPad = 0
# x_src = collect(5:10:nx-2nPml-5)
# z_src = 2ones(Int64, size(x_src))
# x_rec = collect(5:100-nPml)
# z_rec = 2ones(Int64, size(x_rec))

x_src = [100-nPml]
z_src = [100-nPml]

z = (5:10:nz-2nPml-5)|>collect
x = (5:10:nx-2nPml-5)|>collect
x_rec, z_rec = np.meshgrid(x, z)
x_rec = x_rec[:]
z_rec = z_rec[:]

# x_rec = collect(5:1:nx-2nPml-5)
# z_rec = 60ones(Int64, size(x_rec))

para_fname = "./para_file.json"
survey_fname = "./survey_file.json"
data_dir_name = "./Data"
paraGen(nz, nx, dz, dx, nSteps, dt, f0, nPml, nPad, filter_para, isAc, para_fname, survey_fname, data_dir_name)
surveyGen(z_src, x_src, z_rec, x_rec, survey_fname)

cp = 3000ones(nz, nx)
# cp = (1. .+ 0.1*rand(nz, nx)) .* 3000.
cs = zeros(nz, nx)
den = 1000.0 .* ones(nz, nx)

tf_cp = constant(cp)
tf_cs = constant(cs)
tf_den = constant(den)

src = Matrix{Float64}(undef, 1, 2001)
src[1,:] = Float64.(reinterpret(Float32, read("./Src/params/ricker_10Hz.bin")))
tf_stf = constant(repeat(src, outer=length(z_src)))
tf_para_fname = tf.strings.join([para_fname])
tf_gpu_id0 = constant(0, dtype=Int32)
tf_gpu_id1 = constant(1, dtype=Int32)
tf_shot_ids0 = constant(collect(Int32, 0:length(x_src)-1), dtype=Int32)
tf_shot_ids1 = constant(collect(Int32, 13:25), dtype=Int32)

res1 = fwi_obs_op(tf_cp, tf_cs, tf_den, tf_stf, tf_gpu_id0, tf_shot_ids0, tf_para_fname)
# res2 = fwi_obs_op(tf_cp2, tf_cs2, tf_den2, tf_stf, tf_gpu_id1, tf_shot_ids0, tf_para_fname)

sess=Session();init(sess);
@time run(sess, res1)
error("")
# error("")

# function obj()
#     res = 0.0
#     for i = 1:29
#         gpu_id = mod(i, 2)
#         res += fwi_op(tf_cp, tf_cs, tf_den, tf_stf, constant(gpu_id, dtype=Int32), constant([i], dtype=Int32), tf_para_fname)
#     end
#     return res
# end
# J = obj()

# J1 = fwi_op(tf_cp, tf_cs, tf_den, tf_stf, tf_gpu_id0, tf_shot_ids0, tf_para_fname)
# J2 = fwi_op(tf_cp, tf_cs, tf_den, tf_stf, tf_gpu_id1, tf_shot_ids0, tf_para_fname)
# J = J1 + J2
# # config = tf.ConfigProto()
# # config.allow_growth
# # config.intra_op_parallelism_threads = 2
# # config.inter_op_parallelism_threads = 2
# sess=Session();init(sess);
# # @time run(sess, J)
# gg = gradients(J1, tf_cp)
# grad_cp = run(sess, gg)
# imshow(grad_cp);colorbar();

# error("")

# gradient check -- v
function scalar_function(m)
    return fwi_op(m, tf_cs, tf_den, tf_stf, tf_gpu_id0, tf_shot_ids0, tf_para_fname)
end

# open("./Data/Shot1.bin","w") do f
#     write(f, zeros(nz*nx,1))
# end

m_ = constant(3100ones(nz, nx))
# m_ = constant(cp)
# v_ = 100. .* (1. .+ rand(Float32, 384, 134))

v0 = zeros(nz, nx)
# PLEASE!!!!!!!!!!!!!! Don't perturb in the CPML region!!!!!!!!!!!!!!!!!!!!!!!
v0[nPml+5:nz-nPml-5, nPml+5:nx-nPml-5] .= 1.0
v_ = constant(Float64.((1. .+ 0.1*rand(nz, nx)) .* 500 .* v0))
y_ = scalar_function(m_)
dy_ = gradients(y_, m_)
ms_ = Array{Any}(undef, 5)
ys_ = Array{Any}(undef, 5)
s_ = Array{Any}(undef, 5)
w_ = Array{Any}(undef, 5)
gs_ =  @. 1 / 10^(1:5)

for i = 1:5
    g_ = gs_[i]
    ms_[i] = m_ + g_ * v_
    ys_[i] = scalar_function(ms_[i])
    s_[i] = ys_[i] - y_
    w_[i] = s_[i] - g_*sum(v_.*dy_)
end

sess = Session()
init(sess)
sval_ = run(sess, s_)
wval_ = run(sess, w_)
# error("")
sval_ = [x[1] for x in sval_]
wval_ = [x[1] for x in wval_]
close("all")
loglog(gs_, abs.(sval_), "*-", label="finite difference")
loglog(gs_, abs.(wval_), "+-", label="automatic differentiation")
loglog(gs_, gs_.^2 * 0.5*abs(wval_[1])/gs_[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
loglog(gs_, gs_ * 0.5*abs(sval_[1])/gs_[1], "--",label="\$\\mathcal{O}(\\gamma)\$")

plt.gca().invert_xaxis()
legend()
xlabel("\$\\gamma\$")
ylabel("Error")


text:
module InverseProblems

export
    InverseProblem,
    forward_map,
    forward_run!,
    observation_map,
    observation_map_variance_across_time,
    ConcatenatedOutputMap

using OffsetArrays, Statistics, OrderedCollections
using Suppressor: @suppress

using ..Transformations: transform_field_time_series
using ..Parameters: new_closure_ensemble, transform_to_constrained, build_parameters_named_tuple

using ..Observations:
    AbstractObservation,
    SyntheticObservations,
    initialize_forward_run!,
    FieldTimeSeriesCollector,
    observation_times,
    forward_map_names

using Oceananigans: run!, fields, FieldTimeSeries, CPU
using Oceananigans.Architectures: architecture
using Oceananigans.OutputReaders: InMemory
using Oceananigans.Fields: interior, location
using Oceananigans.Grids: Flat, Bounded,
                          Face, Center,
                          RectilinearGrid, offset_data,
                          topology, halo_size,
                          interior_parent_indices

using Oceananigans.Models.HydrostaticFreeSurfaceModels: SingleColumnGrid, YZSliceGrid, ColumnEnsembleSize

import ..Transformations: normalize!

#####
##### Output maps (maps from simulation output to observation space)
#####

output_map_type(fp) = output_map_str(fp)

"""
    struct ConcatenatedOutputMap

Forward map transformation of simulation output to the concatenated
vectors of the simulation output.
"""
struct ConcatenatedOutputMap end

output_map_str(::ConcatenatedOutputMap) = "ConcatenatedOutputMap"

#####
##### InverseProblems
#####

struct InverseProblem{F, O, S, T, P, I}
    observations :: O
    simulation :: S
    time_series_collector :: T
    free_parameters :: P
    output_map :: F
    initialize_simulation :: I
end

nothingfunction(simulation) = nothing

"""
    InverseProblem(observations,
                   simulation,
                   free_parameters;
                   output_map = ConcatenatedOutputMap(),
                   time_series_collector = nothing,
                   initialize_simulation = nothingfunction)

Return an `InverseProblem`.
"""
function InverseProblem(observations,
                        simulation,
                        free_parameters;
                        output_map = ConcatenatedOutputMap(),
                        time_series_collector = nothing,
                        initialize_simulation = nothingfunction)

    if isnothing(time_series_collector) # attempt to construct automagically
        simulation_fields = fields(simulation.model)
        collected_fields = NamedTuple(name => simulation_fields[name] for name in forward_map_names(observations))
        time_series_collector = FieldTimeSeriesCollector(collected_fields, observation_times(observations))
    end

    return InverseProblem(observations, simulation, time_series_collector, free_parameters, output_map, initialize_simulation)
end

Base.summary(ip::InverseProblem) =
    string("InverseProblem{", summary(ip.output_map), "} with free parameters ", ip.free_parameters.names)

function Base.show(io::IO, ip::InverseProblem)
    sim_str = "Simulation on $(summary(ip.simulation.model.grid)) with Δt=$(ip.simulation.Δt)"
    out_map_str = summary(ip.output_map)

    print(io, summary(ip), '\n',
        "├── observations: $(summary(ip.observations))", '\n',
        "├── simulation: $sim_str", '\n',
        "├── free_parameters: $(summary(ip.free_parameters))", '\n',
        "└── output map: $out_map_str")

    return nothing
end

"""
    expand_parameters(ip, θ::Vector)

Convert parameters `θ` to `Vector{<:NamedTuple}`, where the elements
correspond to `ip.free_parameters`.

`θ` may represent an ensemble of parameter sets via:

* `θ::Vector{<:Vector}` (caution: parameters must be ordered correctly!)
* `θ::Matrix` (caution: parameters must be ordered correctly!)
* `θ::Vector{<:NamedTuple}`

or a single parameter set if `θ::Vector{<:Number}`.

If `length(θ)` is less the the number of ensemble members in `ip.simulation`, the
last parameter set is copied to fill the parameter set ensemble.
"""
function expand_parameters(ip, θ::Vector)
    Nfewer = Nensemble(ip) - length(θ)
    Nfewer < 0 && throw(ArgumentError("There are $(-Nfewer) more parameter sets than ensemble members!"))

    θ = [build_parameters_named_tuple(ip.free_parameters, θi) for θi in θ]

    # Fill out parameter set ensemble
    Nfewer > 0 && append!(θ, [θ[end] for _ = 1:Nfewer])

    return θ
end

# Expand single parameter set
expand_parameters(ip, θ::Union{NamedTuple, Vector{<:Number}}) = expand_parameters(ip, [θ])

# Convert matrix to vector of vectors
expand_parameters(ip, θ::Matrix) = expand_parameters(ip, [θ[:, k] for k = 1:size(θ, 2)])

#####
##### Forward map evaluation given vector-of-vector (one parameter vector for each ensemble member)
#####

const OneDimensionalEnsembleGrid = RectilinearGrid{<:Any, Flat, Flat, Bounded}
const TwoDimensionalEnsembleGrid = RectilinearGrid{<:Any, Flat, Bounded, Bounded}

Nobservations(grid::OneDimensionalEnsembleGrid) = grid.Ny
Nobservations(grid::TwoDimensionalEnsembleGrid) = 1

Nensemble(grid::Union{OneDimensionalEnsembleGrid, TwoDimensionalEnsembleGrid}) = grid.Nx
Nensemble(ip::InverseProblem) = Nensemble(ip.simulation.model.grid)

"""
    observation_map(ip::InverseProblem)

Transform and return `ip.observations` appropriate for `ip.output_map`.
"""
observation_map(ip::InverseProblem) = observation_map(ip.output_map, ip.observations)

observation_map(map::ConcatenatedOutputMap, observations) = transform_time_series(map, observations)

"""
    forward_run!(ip, parameters)

Initialize `ip.simulation` with `parameters` and run it forward. Output is stored
in `ip.time_series_collector`.
"""
function forward_run!(ip::InverseProblem, parameters; suppress=false)
    observations = ip.observations
    simulation = ip.simulation
    closures = simulation.model.closure

    # Ensure there are enough parameters for ensemble members in the simulation
    θ = expand_parameters(ip, parameters)

    # Set closure parameters
    simulation.model.closure = new_closure_ensemble(closures, θ, architecture(simulation.model.grid))

    initialize_forward_run!(simulation, observations, ip.time_series_collector, ip.initialize_simulation)

    if suppress
        @suppress run!(simulation)
    else
        run!(simulation)
    end

    return nothing
end

"""
    forward_map(ip, parameters)

Run `ip.simulation` forward with `parameters` and return the data,
transformed into an array format expected by `EnsembleKalmanProcesses.jl`.
"""
function forward_map(ip, parameters; suppress=true)

    # Run the simulation forward and populate the time series collector
    # with model data.
    forward_run!(ip, parameters; suppress)

    # Verify that data was collected properly
    all(ip.time_series_collector.times .≈ ip.time_series_collector.collection_times) ||
        error("FieldTimeSeriesCollector.collection_times does not match FieldTimeSeriesCollector.times. \n" *
              "Field time series data may not have been properly collected")

    # Transform the model data according to `ip.output_map` into
    # the array format expected by EnsembleKalmanProcesses.jl
    # The result has `size(output) = (output_size, ensemble_capacity)`,
    # where `output_size` is determined by both the `output_map` and the
    # data collection dictated by `ip.observations`.
    output = transform_forward_map_output(ip.output_map, ip.observations, ip.time_series_collector)

    # (Nobservations, Nensemble)
    return output
end

(ip::InverseProblem)(θ) = forward_map(ip, θ)

"""
    inverting_forward_map(ip::InverseProblem, X)

Transform unconstrained parameters `X` into constrained,
physical-space parameters `θ` and execute `forward_map(ip, θ)`.
"""
function inverting_forward_map(ip::InverseProblem, X)
    θ = transform_to_constrained(ip.free_parameters.priors, X)
    return forward_map(ip, θ)
end

#####
##### ConcatenatedOutputMap
#####

"""
    transform_time_series(::ConcatenatedOutputMap, observation::SyntheticObservations)

Transforms, normalizes, and concatenates data for field time series in `observation`.
"""
function transform_time_series(::ConcatenatedOutputMap, observation::SyntheticObservations)
    data_vector = []

    for field_name in forward_map_names(observation)
        # Transform time series data observation-specified `transformation`
        field_time_series = observation.field_time_serieses[field_name]
        transformation = observation.transformation[field_name]
        transformed_datum = transform_field_time_series(transformation, field_time_series)

        # Build out array
        push!(data_vector, transformed_datum)
    end

    # Concatenate!
    concatenated_data = hcat(data_vector...)

    return Matrix(transpose(concatenated_data))
end

"""
    transform_time_series(map, batched_observations::Vector)

Concatenate the output of `transform_time_series` of each observation
in `batched_observations`.
"""
transform_time_series(map, batched_observations::Vector) =
    vcat(Tuple(transform_time_series(map, obs) for obs in batched_observations)...)

const BatchedOrSingletonObservations = Union{SyntheticObservations,
                                             Vector{<:SyntheticObservations}}

function transform_forward_map_output(map::ConcatenatedOutputMap,
                                      observations::BatchedOrSingletonObservations,
                                      time_series_collector)

    # transposed_output isa Vector{SyntheticObservations} where SyntheticObservations is Nx by Nz by Nt
    transposed_forward_map_output = transpose_model_output(time_series_collector, observations)

    return transform_time_series(map, transposed_forward_map_output)
end

vectorize(observation) = [observation]
vectorize(observations::Vector) = observations

# Dispatch transpose_model_output based on collector grid
transpose_model_output(time_series_collector, observations) =
    transpose_model_output(time_series_collector.grid, time_series_collector, observations)

transpose_model_output(collector_grid::YZSliceGrid, time_series_collector, observations) =
    SyntheticObservations(time_series_collector.field_time_serieses,
                          observations.forward_map_names,
                          collector_grid,
                          time_series_collector.times,
                          nothing,
                          nothing,
                          observations.transformation)

"""
    transpose_model_output(collector_grid, time_series_collector, observations)

Transpose a `NamedTuple` of 4D `FieldTimeSeries` model output collected by `time_series_collector`
into a Vector of `SyntheticObservations` for each member of the observation batch.

Return a 1-vector in the case of singleton observations.
"""
function transpose_model_output(collector_grid::SingleColumnGrid, time_series_collector, observations)
    observations = vectorize(observations)
    times = time_series_collector.times

    transposed_output = []

    Nensemble = collector_grid.Nx
    Nbatch =  collector_grid.Ny
    Nz = collector_grid.Nz
    Hz = collector_grid.Hz
    Nt = length(times)

    grid = drop_y_dimension(collector_grid)

    for j = 1:Nbatch
        observation = observations[j]
        time_serieses = OrderedDict{Any, Any}()

        for name in forward_map_names(observation)
            loc = LX, LY, LZ = location(observation.field_time_serieses[name])
            topo = topology(grid)

            field_time_series = time_series_collector.field_time_serieses[name]

            indices = field_time_series.indices
            raw_data = parent(field_time_series.data)
            data = OffsetArray(view(raw_data, :, j:j, :, :), 0, 0, -Hz, 0)

            time_series = FieldTimeSeries{LX, LY, LZ, InMemory}(data, grid, nothing, times, indices)
            time_serieses[name] = time_series
        end

        # Convert to NamedTuple
        time_serieses = NamedTuple(name => time_series for (name, time_series) in time_serieses)

        batch_output = SyntheticObservations(time_serieses,
                                             observation.forward_map_names,
                                             grid,
                                             times,
                                             nothing,
                                             nothing,
                                             observation.transformation)

        push!(transposed_output, batch_output)
    end

    return transposed_output
end

function drop_y_dimension(grid::SingleColumnGrid)
    new_size = ColumnEnsembleSize(Nz=grid.Nz, ensemble=(grid.Nx, 1), Hz=grid.Hz)
    new_halo_size = ColumnEnsembleSize(Nz=1, Hz=grid.Hz)
    z_domain = (grid.zᵃᵃᶠ[1], grid.zᵃᵃᶠ[grid.Nz])
    new_grid = RectilinearGrid(size=new_size, halo=new_halo_size, z=z_domain, topology=(Flat, Flat, Bounded))
    return new_grid
end

#####
##### ConcatenatedVectorNormMap
#####

"""
    ConcatenatedVectorNormMap()

Forward map transformation of simulation output to a scalar by
taking a naive `norm` of the difference between concatenated vectors of the
observations and simulation output.
"""
struct ConcatenatedVectorNormMap end

output_map_str(::ConcatenatedVectorNormMap) = "ConcatenatedVectorNormMap"
observation_map(::ConcatenatedVectorNormMap, observations) = reshape([0], 1, 1)

function transform_forward_map_output(::ConcatenatedVectorNormMap, obs, time_series_collector)
    # Collected concatenated output and observations
    G = transform_forward_map_output(ConcatenatedOutputMap(), obs, time_series_collector)
    y = observation_map(ConcatenatedOutputMap(), obs)

    # Compute vector norm across ensemble members. result should be
    # (1, Nensemble)
    return mapslices(Gᵏ -> norm(Gᵏ - y), G, dims=1)
end

observation_map(map::ConcatenatedVectorNormMap, observations) = hcat(0.0)

#####
##### Utils
#####

"""
    observation_map_variance_across_time(map::ConcatenatedOutputMap, observation::SyntheticObservations)

Return an array of size `(Nensemble, Ny * Nz * Nfields, Ny * Nz * Nfields)` that stores the covariance of
each element of the observation map measured across time, for each ensemble member, where `Nensemble` is
the ensemble size, `Ny` is either the number of grid elements in `y` or the batch size, `Nz` is the number
of grid elements in the vertical, and `Nfields` is the number of fields in `observation`.
"""
function observation_map_variance_across_time(map::ConcatenatedOutputMap, observation::SyntheticObservations)
    # These aren't right because every field can have a different transformation, so...
    Nx, Ny, Nz = size(observation.grid)
    Nt = length(first(observation.transformation).time)

    Nfields = length(forward_map_names(observation))

    y = transform_time_series(map, observation)
    @assert length(y) == Nx * Ny * Nz * Nt * Nfields # otherwise we're headed for trouble...

    y = transpose(y) # (Nx, Ny*Nz*Nt*Nfields)

    # Transpose `Nfields` dimension
    permuted_y = permutedims(y, [1, 2, 4, 3])
    reshaped_permuted_y = reshape(permuted_y, Nx, Ny * Nz * Nfields, Nt)

    # Compute `var`iance across time
    dataset = [reshape(var(reshaped_permuted_y[:, :, 1:n], dims = 3), Nx, Ny * Nz, Nfields) for n = 1:Nt]
    concatenated_dataset = cat(dataset..., dims = 2)
    replace!(concatenated_dataset, NaN => 0) # variance for first time step is zero

    return reshape(concatenated_dataset, Nx, Ny * Nz * Nt * Nfields)
end

observation_map_variance_across_time(map::ConcatenatedOutputMap, observations::Vector) =
    hcat(Tuple(observation_map_variance_across_time(map, observation) for observation in observations)...)

observation_map_variance_across_time(ip::InverseProblem) = observation_map_variance_across_time(ip.output_map, ip.observations)

end # module


text:
function poly(scene::makie, points, attributes::Dict)
    last(points) == first(points) && pop!(points)
    polys = GeometryTypes.split_intersections(points)
    result = []
    for poly in polys
        mesh = GLNormalMesh(poly) # make polygon
        if !isempty(GeometryTypes.faces(mesh)) # check if polygonation has any faces
            push!(result, GLVisualize.visualize(mesh, Style(:default), kw_args))
        else
            warn("Couldn't draw the polygon: $points")
        end
    end
    result
end

function shape(d, kw_args)
    points = Plots.extract_points(d)
    result = []
    for rng in iter_segments(d[:x], d[:y])
        ps = points[rng]
        meshes = poly(ps, kw_args)
        append!(result, meshes)
    end
    result
end


text:
using Plots, Statistics, JLD, LaTeXStrings, StatsBase
function movement(first_pos, step_num)
    directions=[[0,1],[0,-1],[1,0],[-1,0]]
    walking=rand(directions, step_num)
    return first_pos .+ cumsum(hcat(walking),dims=1)[step_num]
end
function Rg(data)
    COM=(mean(data, dims=1))
    return sqrt(mean((data[:,1] .- COM[1]) .^2 + (data[:,2] .- COM[2]) .^2))
end

step_num_list=[i for i in 100:100:1000]
first_pos=[0,0]
run_num=100000
AllData=[]
MeanSquaredDistance=[]
RgData=[]
for step_num in step_num_list
    final_destinations=zeros(Int, run_num, 2)
    squared_distance=[]
    for i in 1:run_num
        final_destinations[i,:]=movement(first_pos, step_num)
        push!(squared_distance, (final_destinations[i,1]^2 + final_destinations[i,2]^2))
    end
    push!(SavedData, final_destinations)
    push!(MeanSquaredDistance, mean(squared_distance))
    push!(RgData, Rg(final_destinations))
end
#Radius of Gyration plot:
scatter(log.(RgData), log.(step_num_list), legend=false, xlabel="log(Time)", ylabel="log(Rg)",
    title="Logarithm of Radious of Gyration over the Logarithm of Time", dpi=400, titlefontsize=9)
savefig("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter5\\4.5\\Figs\\Rg_Time.png")
# mean(r^2) plot:
scatter(MeanSquaredDistance, step_num_list, ylabel=L"<r^2>", xlabel="Time",
    title="Mean Squared Distance over Time", legend=false, dpi=400)
savefig("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter5\\4.5\\Figs\\MeanSquare(r)_Time.png")

save("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter5\\4.5\\2DWalker.jld",
    "data", AllData, "MeanSquaredDistanceData",MeanSquaredDistance, "RgData",RgData )


text:
# Towers.jl by Alex Van Mechelen
# Creation Date: 10/04/2020
# Made for Julia 1.2.0

global const AUTO_START_GAME = true # Change startup behavior when included in REPL

# Colors used in printstyled()
global const COL_PLAYVALS = :yellow
global const COL_GAMEVALS = :normal
global const COL_PENCIL = :light_cyan
global const COL_GRID = :light_black
global const COL_ERROR = :light_red
global const COL_HIGHLIGHT = :blue
global const COL_QUESTION = :light_green
global const COL_EXPL = :light_blue
global const COL_GUESS = :light_magenta

# Input keys (followed by [ENTER] press (readline()))
global const KEY_UP = ('u','e')
global const KEY_DOWN = ('d',)
global const KEY_LEFT = ('l','s')
global const KEY_RIGHT = ('r','f')
global const KEY_CTRL = ('_','-')
global const KEY_PENCIL = ('p',)
global const KEY_TOGGLE_ERRORS = ('o',)

# Char to display to indicate highLighted line if no clue in that line
global const HIGH_CHAR = "●"

# Maximum dimension for generated 'PLAY' Games
global const MAX_GEN_DIM = 6

# Make higher for more difficult 'easy' levels (I would recommend values from 10 to 13), higher values mean more computation time
global const EASY_DIFF = 10

# Minimum screen dimensions
global const MIN_SCREEN_WIDTH = 83
global const MIN_SCREEN_HEIGHT = 30

# Color of highlighted box in game
global Col_Box = COL_HIGHLIGHT # Not const (Pencil-Mark mode)

# Show outBoardErrors or not
global ShowOutBoardErrors = true # Change this to change the default value

# Represents a game board
mutable struct Board
	boardStr::String
	dim::Int
	clues::Array{Int64,2}
	grid::Array{Int64,2}
	gameVals::Array{Bool,2}
	guessVals::Array{Bool,2}
	highLight::Array{Int64,1}
	highLightLine::Int
	pencilMarks::Array{Array{Int64,1},2}

	# Used to create a new Board from a string ( = Game-ID in Simon Tatham's Portable Puzzle Collection)
	function Board(boardStr::String)
		# Check if valid boardStr
		try
			(dim,code) = split(boardStr,":")
			dim = parse(Int,dim)
			@assert(0<dim<10,"Grid dimension must be between 0 and 10.")
			onlyOutClues = false
			outBoardClues = code #If only outside clues given
			inBoardClues = code #Predefine inBoardClues, otherwise not defined error
			try
				(outBoardClues,inBoardClues) = split(code,",")
			catch #Only outside board clues given (no ",")
				onlyOutClues = true
			end
			arrOutBoardClues = split(outBoardClues,"/")

			# Outside Board Clues
			@assert(Int((length(arrOutBoardClues))/4)==dim,"Invalid Board String. (First Part)")
			intArr = Array{Int,1}(undef,length(arrOutBoardClues))
			for i in 1:length(arrOutBoardClues)
				if isempty(arrOutBoardClues[i])
					intArr[i] = 0
				else
					intArr[i] = parse(Int,arrOutBoardClues[i])
				end
			end
			clues = reshape(intArr,(dim,4))

			# Inside Board Clues
			grid = zeros(Int,(dim,dim))
			gameVals = zeros(Bool,(dim,dim))
			if !onlyOutClues
				j = 1
				for i in inBoardClues
					(subArrIndex,inSubArrIndex) = divrem(j,dim)
					num = Int(collect(i)[1])-Int('a')+1 #Map alphabet to numbers
					if inSubArrIndex == 0
						subArrIndex -= 1
						inSubArrIndex = dim
					end
					if num < -38 #If this character is a number
						num = parse(Int,i)
						grid[subArrIndex+1,inSubArrIndex] = num
						gameVals[subArrIndex+1,inSubArrIndex] = true
						j += 1
					else #Character is not a number
						if num != -1 #Character is not an underscore
							j += num
						end
					end
				end
				@assert(sqrt(j-1)==dim,"Invalid Board String. (Second Part)")
			end

			# Highlight top left corner when initialized
			highLight = [1,1]

			# Do not highLight any line when initialized
			highLightLine = 0

			# No values guessed by solver when initialized
			guessVals = zeros(Bool,(dim,dim))

			return true,new(boardStr,dim,clues,grid,gameVals,guessVals,highLight,highLightLine,initPencilMarks(dim))

		catch #Invalid boardStr -> succes = false
			return false,nothing
		end
	end
end

# Create Array to store Pencil-Marks
function initPencilMarks(dim::Int)
	pm = Array{Array{Int,1},2}(undef,dim,dim)
	for i in 1:length(pm)
		pm[i] = zeros(Int,dim)
	end
	return pm
end

function Base.show(io::IO,b::Board)
	Calibrate()
	drawGrid(io,b)
end

function Base.write(io::IO,b::Board)
	Calibrate()
	drawGrid(io,b)
end

function loadingGame()
	clearScreen()
	printstyled("Loading Game...\n",color=COL_HIGHLIGHT)
end

# To show user input keys without ''
function Base.show(io::IO,t::Tuple)
	write(io,"(")
	for i in 1:length(t)
		write(io,t[i])
		i == length(t) || write(io,",")
	end
	write(io,")")
end

# Show Pencil-Marks if present in highlighted cell
function showPencilMarks(io::IO,pms::Array{Int64,1})
	printstyled("Pencil-Marks on highlighted cell: ",color=COL_PENCIL)
	for i in 1:length(pms)
		if pms[i] != 0
			printstyled(io,pms[i],color=COL_QUESTION)
			i != count(pms.!=0) && printstyled(io,",",color=COL_PENCIL)
		end
	end
end

# Load a new game from a gameStr
function loadGame(gameStr::String)
	succes,newBoard = Board(gameStr)
	if succes
		return true,newBoard
	else
		return false,nothing
	end
end

using Random

shuffleRows(mat) = mat[shuffle(1:end), :]
shuffleCols(mat) = mat[:, shuffle(1:end)]

function addAtDiagonal(mat::Array)
    n = size(mat)[1] + 1
    newMat = similar(mat, size(mat) .+ 1)
    for j in 1:n, i in 1:n
        newMat[i, j] = (i == n && j < n) ? mat[1, j] : (i == j) ? n :
            (i < j) ? mat[i, j - 1] : mat[i, j]
    end
    return newMat
end

function makeLatinSquare(N::Int)
    mat = [2 1 ; 1 2]
    for i in 3:N
        mat = addAtDiagonal(mat)
    end
    return shuffleCols(shuffleRows(mat))
end

# Adds last part to BoardStr
function generateLastPartOfBoardStrFromGrid(boardStr,grid::Array{Int64,2})
	# If no clues inside the grid, the boardStr is already okay
	grid == zeros(Int64,size(grid)) && return true,boardStr

	boardStr *= ","

	countEmpty = 0
	grid = transpose(grid)

	for i in 1:length(grid)
		if countEmpty > Int('y') # Impossible to convert this board to a boardStr
			return false,""
		elseif grid[i] == 0
			countEmpty += 1

			# If last cell reached
			if i == length(grid)
				boardStr *= string(Char(countEmpty+Int('a')-1))
			end
		else
			# Alphabet separator
			if i != 1
				boardStr *= countEmpty == 0 ? "_" : string(Char(countEmpty+Int('a')-1))
			end

			# Numerical value of the clue
			boardStr *= string(grid[i])

			countEmpty = 0 # Reset empty-cell count
		end
	end
	return true,boardStr
end

# NEW VERSION OF GENERATEGAME (GENERATION BY ADDITION) -> Still to improve...
function generateGame(dim::Int,easyMode::Bool)
	succes = false
	numGuessed = 174 # Initialize (non-zero)
	numTries = 3 # Tries to find a good clue to remove
	maxTries = dim*EASY_DIFF

	# Make a Latin Square of the given dimension
	latSqr = makeLatinSquare(dim)

	# Initialize clues
	clues = zeros(Int,(dim,4))

	# Fill in clue values based on Latin Square
	for i in 1:dim
		# Top / Bottom / Left / Right
		clues[i,1] = towersVisible(latSqr[:,i])
		clues[i,2] = towersVisible(reverse(latSqr[:,i]))
		clues[i,3] = towersVisible(latSqr[i,:])
		clues[i,4] = towersVisible(reverse(latSqr[i,:]))
	end

	# Generate BoardStr
	# First part of boardStr is dim
	boardStr = string(dim)*":"

	# Clues (outside the grid) separated by "/"
	for clueIndex in 1:length(clues)
		clues[clueIndex] != 0 && (boardStr *= string(clues[clueIndex]))
		clueIndex != length(clues) && (boardStr *= "/")
	end

	# Clues (inside the grid) separated by "_" (because all next to one another)
	boardStrInsideOnly = deepcopy(boardStr)
	boardStr = generateLastPartOfBoardStrFromGrid(boardStr,latSqr)[2]

	# Initialize game, solvedGame
	game = Board(boardStr)[2]
	solvedGame = deepcopy(game)

	if easyMode
		counter = 0
		while counter < maxTries && (succes || numTries > 0) && solvedGame.guessVals == zeros(Bool,(dim,dim))
			counter += 1
			# Reset BoardStr to its only-inside-version
			boardStr = deepcopy(boardStrInsideOnly)

			# Make a copy of the old game grid
			oldGameGrid = deepcopy(game.grid)

			# Remove a random value inside the grid
			randIndex = rand(1:length(latSqr))
			game.grid[randIndex] = 0

			# Generate a new boardStr
			succes,boardStr = generateLastPartOfBoardStrFromGrid(boardStr,game.grid)
			if !succes
				game.grid = deepcopy(oldGameGrid)
				numTries -= 1
				succes = false
				continue
			end

			# Generate the new game with this boardStr
			lastGame = deepcopy(game)
			succes,game = Board(boardStr)
			succes || error("Invalid BoardStr generated: $(boardStr)")

			# Try to solve the new game
			lastSolvedGame = deepcopy(solvedGame)
			succes,solvedGame = mainSolver(Board(game.boardStr)[2]) # Standard is solve invisibly

			# If unsolvable or guesses had to be made, retry by removing an other value
			if !succes || solvedGame.guessVals != zeros(Bool,(dim,dim))
				game.grid = deepcopy(oldGameGrid)
				numTries -= 1
				succes = false

				# Return last playable easy game
				game = deepcopy(lastGame)
				boardStr = game.boardStr
				solvedGame = deepcopy(lastSolvedGame)
			end
		end
	else
		for i in 1:dim^2
			while true
				# Remove a random value inside the grid
				randIndex = rand(1:length(latSqr)+4*dim)
				if randIndex > length(latSqr)
					randIndex -= length(latSqr)

					# Index of the clue in a line of outBoardClues
					lineIndex = rem(randIndex,dim)

					# Correction for when randIndex is a multiple of dim
					lineIndex == 0 && (lineIndex = dim)

					# Index of the direction from which you're looking at the grid
					dirIndex = div(randIndex,dim)+1

					# Correction for when randIndex = 4*dim
					dirIndex > 4 && (dirIndex = 4)

					# Try to remove it, if not, choose another clue to remove
					if clues[lineIndex,dirIndex] != 0
						clues[lineIndex,dirIndex] = 0
					else
						break
					end
				else
					# Try to remove it, if not, choose another clue to remove
					if game.grid[randIndex] != 0
						game.grid[randIndex] = 0
					else
						break
					end
				end
			end
		end

		# Generate BoardStr
		# First part of boardStr is dim
		boardStr = string(dim)*":"

		# Clues (outside the grid) separated by "/"
		for clueIndex in 1:length(clues)
			clues[clueIndex] != 0 && (boardStr *= string(clues[clueIndex]))
			clueIndex != length(clues) && (boardStr *= "/")
		end

		# Prepare the game for playing
		succes,boardStr = generateLastPartOfBoardStrFromGrid(boardStr,game.grid)
		succes,game = Board(boardStr)
		succes || error("Invalid BoardStr generated: $(boardStr)")
		succes,solvedGame = mainSolver(Board(game.boardStr)[2]) # Standard is solve invisibly
	end

	return game,solvedGame
end

# OLD VERSION OF GENERATEGAME (GENERATION BY ADDITION)
function OLDgenerateGame(dim::Int,easyMode::Bool)
	succes = false
	numGuessed = 174 # Initialize (non-zero)
	firstTime = true
	boardStr = ""

	# Initialize game, solvedGame
	game = Board("1:///")[2]
	solvedGame = deepcopy(game)

	while !succes || (numGuessed != 0 && easyMode)
		if firstTime
			firstTime = false

			# Make a Latin Square of the given dimension
			latSqr = makeLatinSquare(dim)

			# Initialize clues
			clues = zeros(Int,(dim,4))

			# Fill in clue values based on Latin Square
			for i in 1:dim
				# Top / Bottom / Left / Right
				clues[i,1] = towersVisible(latSqr[:,i])
				clues[i,2] = towersVisible(reverse(latSqr[:,i]))
				clues[i,3] = towersVisible(latSqr[i,:])
				clues[i,4] = towersVisible(reverse(latSqr[i,:]))
			end

			# Generate BoardStr
			# First part of boardStr is dim
			boardStr = string(dim)*":"

			# Clues (outside the grid) separated by "/"
			for clueIndex in 1:length(clues)
				clues[clueIndex] != 0 && (boardStr *= string(clues[clueIndex]))
				clueIndex != length(clues) && (boardStr *= "/")
			end
		end

		succes,game = Board(boardStr)
		succes || error("Invalid BoardStr generated: $(boardStr)")

		succes,solvedGame = mainSolver(Board(game.boardStr)[2]) # Standard is solve invisibly

		# Temporary copy of solvedGame guessVals
		copySolvedGuessVals = deepcopy(solvedGame.guessVals)

		# For easy games, test if number of guesses is zero, else add clue inside grid
		guessedValueIndices = findall(reshape(transpose(copySolvedGuessVals),length(copySolvedGuessVals)))
		numGuessed = length(guessedValueIndices)

		# On difficulty [Hard], generate new game if no guesses were needed to solve puzzle
		if !easyMode && numGuessed == 0
			succes = false
			firstTime = true # Generate new boardStr
		end

		easyMode || continue # Go to next while-iteration if hard mode
		numGuessed == 0 && continue # Go to next while-iteration if no guesses were made while solving

		# Pick a random index to fill in an extra in-board clue
		fillInIndex = rand(guessedValueIndices)

		if numGuessed != 0 # No need to include && easyMode, because checked earlier
			solvedGame.grid = transpose(solvedGame.grid)
			toFillIn = solvedGame.grid[fillInIndex]
			solvedGame.grid = transpose(solvedGame.grid)

			# Make a backup of the boardStr
			copyBoardStr = deepcopy(boardStr)

			# Remove in-grid clues from boardStr and generate new coded in-grid clues
			boardStr = split(boardStr,",")[1]
			boardStr *= ","

			# Add new clue to the game
			game.grid = transpose(game.grid)
			game.grid[fillInIndex] = toFillIn
			game.grid = transpose(game.grid)

			# Make a temporary copy of the grid
			gridCopy = deepcopy(game.grid)

			# Reshape the game grid to one long vector
			reshGame = reshape(gridCopy,length(gridCopy))

			# Calculate new end of gameStr
			countEmpty = 0

			game.grid = transpose(game.grid)
			for i in 1:length(game.grid)
				if countEmpty > Int('y') # Impossible to convert this board to a boardStr
					boardStr = deepcopy(copyBoardStr)
					break # Pick another cell
				elseif game.grid[i] == 0
					countEmpty += 1

					# If last cell reached
					if i == length(game.grid)
						boardStr *= string(Char(countEmpty+Int('a')-1))
					end
				else
					# Alphabet separator
					boardStr *= countEmpty == 0 ? "_" : string(Char(countEmpty+Int('a')-1))

					# Numerical value of the clue
					boardStr *= string(game.grid[i])

					countEmpty = 0 # Reset empty-cell count
				end
			end
		end
	end

	return game,solvedGame
end

function convHlLine(game::Board,hlL::Int)
	dim = game.dim

	# Convert hlL index to line/col index
	convHlL = (hlL-1)%dim+1

	# Check what line/col to highLight
	rowLHl = (hlL in 1:dim)
	rowRHl = (hlL in dim+1:2*dim)
	colTHl = (hlL in 2*dim+1:3*dim)
	colBHl = (hlL in 3*dim+1:4*dim)

	return convHlL,rowLHl,rowRHl,colTHl,colBHl
end

# Print full game screen (including emptyLines, errorMessages)
function drawGrid(io::IO,b::Board)
	inErrors,outErrors,gameCorrect = getErrors(b)

	# Shorter names for frequently used variables
	hl = b.highLight
	dim = b.dim
	pms = b.pencilMarks

	convHlL,rowLHl,rowRHl,colTHl,colBHl = convHlLine(b,b.highLightLine)

    (!(0 < hl[1] <= dim) || !(0 < hl[2] <= dim)) && errorMessage("Cannot highlight that cell")
    indHl = 0 # Standard -> no highLight in this line

	# Above grid
	printstyled("Boardstring: $(b.boardStr)",color=COL_HIGHLIGHT)
	emptyLines(Int(height/2-dim-3))
	hlIndex = colTHl ? convHlL : 0
	drawClueSeg(io,b.dim,b.clues[:,1],outErrors[:,1],hlIndex)

	# Grid itself
    for i in 1:dim+1
        chars = ("├","┼","┤") # Standard -> 'linking' chars

        i-1 != hl[1] && (indHl = 0) # Only reset to zero if highLight occurred two horizontal lines ago
		i == hl[1] && (indHl = hl[2]) # If highLight in this row -> Set indHl to column value of hl

		# Top or Bottom of grid -> 'partial-linking' chars
        i == 1 && (chars = ("┌","┬","┐"))
        i == dim+1 && (chars = ("└","┴","┘"))

        drawLineSeg(io,dim,chars,indHl)

        i-1 == hl[1] && (indHl = 0) # If preceding line had highLight, now DO reset it to zero to avoid highlighting the grid VALUE below
		# After bottom grid edge, don't try to draw in-grid values
		doHL = (convHlL == i)
        i != dim+1 && drawValSeg(io,dim,b.grid[i,:],b.gameVals[i,:],b.guessVals[i,:],pms[i,:],inErrors[i,:],[outErrors[i,3],outErrors[i,4]],zeroToSpace(b.clues[i,3],rowLHl&&i==convHlL),zeroToSpace(b.clues[i,4],rowRHl&&i==convHlL),doHL && rowLHl,doHL && rowRHl,indHl)
    end

	# Below grid
	hlIndex = colBHl ? convHlL : 0
	drawClueSeg(io,b.dim,b.clues[:,2],outErrors[:,2],hlIndex)
	emptyLines(Int(height/2-dim)-4)
	count(pms[hl...].!=0) != 0 && showPencilMarks(io,pms[hl...])
	emptyLines(2)
	gameCorrect && gameWon()
end

# Print segment of clues (top & bottom)
function drawClueSeg(io::IO,dim::Int,clues::Array{Int,1},outErrors::Array{Bool,1},hlIndex::Int=0)
	spaces(Int(width/2-2*dim)-1)
	colVals = COL_GAMEVALS
	for i in 1:length(clues)
		ShowOutBoardErrors && outErrors[i] == 1 && (colVals = COL_ERROR) # If this clue value has an error-label
		i == hlIndex && (colVals = COL_EXPL)
		printstyled(io," "^3,zeroToSpace(clues[i],i==hlIndex),color = colVals)
		colVals = COL_GAMEVALS
	end
	print("\n")
end

# Print horizontal segment of game-grid (only lines, no numbers)
function drawLineSeg(io::IO,dim::Int,chars::Tuple,ind::Int=0)
    ind > dim && errorMessage("BoundsError: Index exceeds dimension")
	spaces(Int(width/2-2*dim))

    first = chars[1]; middle = chars[2]; final = chars[3]
    line = "─"^3; highlight = "═"^3

	printstyled(io,first,color = COL_GRID)
	if ind != 0
		printstyled(io,(line*middle)^(ind-1),color = COL_GRID) # Print standard lines until highLight position
		printstyled(io,highlight,color = Col_Box)
		if ind != dim # Print rest of grid if not last cell highlighted
			printstyled(io,middle,color = COL_GRID)
			printstyled(io,(line*middle)^(dim-ind-1),color = COL_GRID)
        end
    else # No highLights on this line segment
		printstyled(io,(line*middle)^(dim-1),color = COL_GRID)
    end
	ind != dim && printstyled(io,line,color = COL_GRID) # Draw last normal line segment if that cell shouldn't be highlighted
    printstyled(io,final*"\n",color = COL_GRID)
end

# Print horizontal segment of game-grid (lines & numbers (in & out game-grid))
function drawValSeg(io::IO,dim::Int,values::Array{Int64,1},gameVals::Array{Bool,1},guessVals::Array{Bool,1},pms::Array{Array{Int64,1},1},inErrors::Array{Bool,1},outErrors::Array{Bool,1},clueFirst::String,clueLast::String,rowLHl::Bool,rowRHl::Bool,ind::Int=0)
    ind > dim && errorMessage("BoundsError: Index exceeds dimension")
    length(values) != dim && errorMessage("Dimensions do not match")
	length(outErrors) > 2 && error("outErrors Array should not have more than two elements in drawValSeg.")

	spaces(Int(width/2-2*dim)-2)

    line = "│"; highlight = "║"
	useLine = line # Which of the above to use

	cluesVisible = 3
	clues = Array{Int,1}(undef,cluesVisible)

	colGrid = COL_GRID; colVals = COL_GAMEVALS

	# Print preceding (left) clue
	ShowOutBoardErrors && outErrors[1] == 1 && (colVals = COL_ERROR)
	rowLHl && (colVals = COL_EXPL)
	printstyled(io,clueFirst*" ",color = colVals)

	colVals = COL_GAMEVALS # Reset after potential change above

	# Grid segment with values (or spaces)
    for i in 1:dim
		val = " "*string(values[i])*" "

        i-1 != ind && (useLine = line; colGrid = COL_GRID) # Reset if previous cell not highlighted
        i == ind && (useLine = highlight; colGrid = Col_Box)

        printstyled(io,useLine,color = colGrid)

		gameVals[i] == 1 || (colVals = COL_PLAYVALS) # Draw values not entered by user in other color
		inErrors[i] == 1 && (colVals = COL_ERROR)
		guessVals[i] == 1 && (colVals = COL_GUESS)

		try
			clues = pms[i][1:cluesVisible] #[1:cluesVisible] To make only first $cluesVisible clues visible, otherwise grid had to be much larger
		catch # Grid dimension is lower than $cluesVisible -> add zeros to the end
			clues = [pms[i]...,zeros(Int64,cluesVisible-length(pms[i]))...]
		end

		values[i] == 0 && (val = zeroToSpace(clues); colVals = COL_PENCIL) # No grid value -> draw Pencil-Marks
        printstyled(io,val,color = colVals)
		colVals = COL_GAMEVALS # Reset for next iteration
    end

    dim-1 == ind && (useLine = line; colGrid = COL_GRID)
    printstyled(io,useLine,color = colGrid)

	# Print posterior (right) clue
	ShowOutBoardErrors && outErrors[2] == 1 && (colVals = COL_ERROR)
	rowRHl && (colVals = COL_EXPL)
	printstyled(io," "*clueLast*"\n",color = colVals)
end

# Returns string of values given an Array, zeros become spaces
function zeroToSpace(arrVals::Array{Int,1})
	res = ""
	for elt in arrVals
		res *= zeroToSpace(elt)
	end
	return res
end

# Returns string given integer, zeros become spaces
function zeroToSpace(intClue::Int,highLighted::Bool=false)
	highLighted && intClue == 0 && return HIGH_CHAR
	return intClue == 0 ? " " : string(intClue)
end

# Get all the in & out errors of a Board
function getErrors(b::Board)
	inErrors = getInBoardErrors(b)
	outErrors = getOutBoardErrors(b) # Check correct num of towers visible
	gameCorrect = false
	isempty(findall(inErrors .== true)) && isempty(findall(outErrors .== true)) && isempty(findall(b.grid .== 0)) && (gameCorrect = true)
	return inErrors,outErrors,gameCorrect
end

# Returns Array of Booleans (1 = error (more than one occurrence of a number in same row/column)), same size as Board.grid
function getInBoardErrors(b::Board)
	boolErrorGrid = zeros(Bool,size(b.grid))
	for i in 1:b.dim # For all the rows
		boolArr = getDoubles(b.grid[i,:])
		boolErrorGrid[i,:] = boolArr
	end
	for i in 1:b.dim # For all the columns
		boolArr = getDoubles(b.grid[:,i])
		boolErrorGrid[:,i] = (boolErrorGrid[:,i] .| boolArr) # Bitwise or to keep earlier detected errors
	end
	return boolErrorGrid
end

function getDoubles(a::Array{Int64,1},boolArr::Array{Bool,1}=zeros(Bool,size(a)))
    size(a) != size(boolArr) && error("Both arrays must have the same size.")
    checked = () # Remember already checked doubles to speed up the process
    for i in 1:length(a)
        (i in checked || a[i] == 0) && continue # Already checked or 0 -> no need to check (again)
        pair = findall(a .== a[i])
        if length(pair)>1
            checked = (checked...,pair...)
            for i in pair
                boolArr[i] = true
            end
        end
    end
    return boolArr
end

# Returns Array of Booleans (1 = error (based on number of towers visible)), same size as Board.clues
function getOutBoardErrors(b::Board)
	boolErrorClues = zeros(Bool,size(b.clues))
	boolErrorClues[:,1] = checkOutBoardErrors(b.clues[:,1],b.grid)
	boolErrorClues[:,2] = checkOutBoardErrors(b.clues[:,2],b.grid,true)
	boolErrorClues[:,3] = checkOutBoardErrors(b.clues[:,3],transpose(b.grid))
	boolErrorClues[:,4] = checkOutBoardErrors(b.clues[:,4],transpose(b.grid),true)
	return boolErrorClues
end

# For one of four directions to look at the Board.grid
function checkOutBoardErrors(clues::Array{Int,1},gridVals,rev::Bool=false)
	boolArr = zeros(Bool,size(clues))
	for i in 1:length(clues)
		clues[i] == 0 && continue # No clue given -> continue
		if rev
			boolArr[i] = checkClueError(clues[i],reverse(gridVals[:,i]))
		else
			boolArr[i] = checkClueError(clues[i],gridVals[:,i])
		end
	end
	return boolArr
end

# For individual outBoardClue
function checkClueError(clue::Int,seg::Array{Int64,1})
	return towersVisible(seg)!=clue
end

function towersVisible(seg::Array{Int64,1})
	maxHeight = length(seg)
	numVis = 0
	maxHVis = 0
	for i in 1:maxHeight
		if seg[i] > maxHVis
			maxHVis = seg[i]
			numVis += 1
		end
		seg[i] == maxHeight && return numVis
	end
	return numVis
end

function gameWon()
	printstyled("Congrats! You won the game!\n",color = COL_QUESTION)
	readline()
	return nothing
end

function Calibrate()
	scrHeight,scrWidth = displaysize(stdout) # Determine screen size

	# Screen size too small to display main screen.
	if scrHeight < MIN_SCREEN_HEIGHT || scrWidth < MIN_SCREEN_WIDTH
		clearScreen()
		printstyled("Please enlarge your screen.",color = COL_ERROR)
		while scrHeight < MIN_SCREEN_HEIGHT || scrWidth < MIN_SCREEN_WIDTH
			scrHeight,scrWidth = displaysize(stdout)
		end
		clearScreen()
	end

	global height = Int(round(scrHeight/2)*2) #Make screen size even (easily divisible by two)
	global width = Int(round(scrWidth/2)*2)

	return nothing
end

function emptyLines(n::Int)
	print("\n"^n)
end

function spaces(n::Int)
	print(" "^n)
end

# Print with preceding spaces and specific chars highlighted, used in Main Menu
function spacePrint(n::Int,str::String,highLights::Tuple=(),newLine::Bool=true)
	spaces(n)
	for character in str
		if character in highLights
			printstyled(character,color=COL_HIGHLIGHT)
		else
			print(character)
		end
	end
	newLine && print("\n")
end

function clearScreen()
	try
		emptyLines(height)
	catch # Global height not yet defined
		emptyLines(100)
	end
end

function errorMessage(msg::String)
	printstyled(msg*".\n",color=COL_ERROR)
end

function errorMessage(msg::String,game::Board)
	print(game)
	errorMessage(msg)
end

function outBoardErrorToggleMessage(ShowOutBoardErrors::Bool)
	if ShowOutBoardErrors
		printstyled("Succesfully enabled ShowOutBoardErrors.\n",color = COL_QUESTION)
	else
		printstyled("Succesfully disabled ShowOutBoardErrors.\n",color = COL_QUESTION)
	end
end

function moveSelection!(dir::Char,game::Board)
	succes = true

	dim = game.dim
	highLight = game.highLight

	# Check direction
	if dir in KEY_UP && highLight[1]>1
		highLight[1] -= 1
	elseif dir in KEY_UP
		errorMessage("Cannot move Up",game)
		succes = false
	elseif dir in KEY_DOWN && highLight[1]<dim
		highLight[1] += 1
	elseif dir in KEY_DOWN
		errorMessage("Cannot move Down",game)
		succes = false
	elseif dir in KEY_LEFT && highLight[2]>1
		highLight[2] -= 1
	elseif dir in KEY_LEFT
		errorMessage("Cannot move Left",game)
		succes = false
	elseif dir in KEY_RIGHT && highLight[2]<dim
		highLight[2] += 1
	elseif dir in KEY_RIGHT
		errorMessage("Cannot move Right",game)
		succes = false
	else
		errorMessage("Invalid character",game)
		succes = false
	end

	succes && (game.highLight = highLight)
	return succes,game
end

# Place a number (not Pencil-Mark) in a cell
function playMove!(num::UInt8,game::Board)
	succes = true
	try
		if game.gameVals[game.highLight...] == 1
			errorMessage("You can't change game clues",game)
			return false,game
		end
		game.grid[game.highLight...] = num
		game.pencilMarks[game.highLight...] = zeros(Int,game.dim)
	catch
		errorMessage("Cannot do that",game)
		succes = false
	end
	return succes,game
end

function getInput(game::Board,solvedGame::Board=game,noSolvedGameGiven::Bool=true)
	succes = false
	input = lowercase(readline())
	Calibrate() # Screen size could have changed while prompting user
	if (length(input) > 1 || isempty(input))
		errorMessage("Please only enter single characters or numbers",game)
	else
		try # Check if number
			input = parse(UInt8,input)
			if input > game.dim
			errorMessage("No digits higher than $(game.dim) allowed",game)
			else
			succes,game = playMove!(input,game)
			end
		catch # Input is non-number
			input = collect(input)[1] # Convert to Char
			if input in KEY_TOGGLE_ERRORS
				global ShowOutBoardErrors = !ShowOutBoardErrors
				print(game)
				outBoardErrorToggleMessage(ShowOutBoardErrors)
			elseif input in KEY_CTRL # CTRL mode
				print(game)
				printstyled("You have entered CTRL mode: Please enter command.\n",color = COL_QUESTION)
				try # Comment for full details of error
					line = readline()
					Calibrate() # Screen size could have changed while prompting user
					lowercase(line) == "exit" && return true,game # Exit command given
					if lowercase(line) in ("solve","solvestep")
						solveStepByStep = lowercase(line)=="solvestep"
						solvable = true
						if noSolvedGameGiven || solveStepByStep
							print(game)
							solveStepByStep || printstyled("Solving Game...\n",color = COL_HIGHLIGHT)
							solvable,solvedGame = mainSolver(Board(game.boardStr)[2],solveStepByStep) # Standard is solve invisibly
							emptyLines(1) # To make it a full screen
							if !solvable
								print(game)
								printstyled("Cannot solve that game.\n",color = COL_ERROR)
							end
						end
						if solvable
							game = deepcopy(solvedGame)
							succes = true
						end
					else
						eval(Meta.parse(line)) # Try to evaluate the command
						succes = true
					end
				catch ex # Comment for full details of error
					print(game)
					printstyled("Something went wrong trying to run that command: $(ex).\n",color = COL_ERROR)
				end # Comment for full details of error
			elseif input in KEY_PENCIL # Pencil-Mark mode
				if game.grid[game.highLight...] == 0
					global Col_Box = COL_PENCIL
					print(game)
					global Col_Box = COL_HIGHLIGHT
					printstyled("You have entered Pencil-Mark mode: Please enter a number between zero and $(game.dim+1) (both excluded).\n",color = COL_PENCIL)
					try
						input = parse(UInt8,readline()) # Results in an error if not a number
						Calibrate() # Screen size could have changed while prompting user
						if input > game.dim || input == 0
							pencilNumError(game)
						elseif input in game.pencilMarks[game.highLight...] # If already Pencil-Mark, remove it
							index = findall(game.pencilMarks[game.highLight...].==input)[1]
							game.pencilMarks[game.highLight...][index] = 0
							sort!(game.pencilMarks[game.highLight...],rev=true)
							succes = true
						elseif 0 in game.pencilMarks[game.highLight...] # If spaces left to place Pencil-Mark, do so
							index = findall(game.pencilMarks[game.highLight...].==0)[1]
							game.pencilMarks[game.highLight...][index] = input
							sort!(game.pencilMarks[game.highLight...],rev=true)
							succes = true
						else
							errorMessage("This cell is already full of pencil-marks, you can remove these by clearing the cell [zero]",game)
						end
					catch
						pencilNumError(game)
					end
				else
					errorMessage("You cannot enter Pencil-Mark mode when the cell already contains a number",game)
				end

			else
				succes,game = moveSelection!(input,game)
			end
		end
		succes && print(game,"\n") #\n makes it as long as if there was an error message
	end
	return false,game # No exit command given
end

function pencilNumError(game::Board)
	errorMessage("Only numbers between zero and $(game.dim+1) allowed (both excluded)",game)
end

# In game -> keep asking for input until gameWon()
function loopGetInput(game::Board,solvedGame::Board=game)
	noSolvedGameGiven = (solvedGame == game)
	print(game,"\n")
	while !getErrors(game)[3] # While there are still errors
		exitGame,game = getInput(game,solvedGame,noSolvedGameGiven)
		exitGame && return true # Exit command given
	end
	return false # No exit command given
end

function printLogo(precedingSpaces::Int=Int(floor((width-55)/2))) #Make logo centered, .... #TODO
	toHighLight = ('╗','╚','═','╔','╝','║')
	spacePrint(precedingSpaces,"████████╗░█████╗░░██╗░░░░░░░██╗███████╗██████╗░░██████╗",toHighLight)
	spacePrint(precedingSpaces,"╚══██╔══╝██╔══██╗░██║░░██╗░░██║██╔════╝██╔══██╗██╔════╝",toHighLight)
	spacePrint(precedingSpaces,"░░░██║░░░██║░░██║░╚██╗████╗██╔╝█████╗░░██████╔╝╚█████╗░",toHighLight)
	spacePrint(precedingSpaces,"░░░██║░░░██║░░██║░░████╔═████║░██╔══╝░░██╔══██╗░╚═══██╗",toHighLight)
	spacePrint(precedingSpaces,"░░░██║░░░╚█████╔╝░░╚██╔╝░╚██╔╝░███████╗██║░░██║██████╔╝",toHighLight)
	spacePrint(precedingSpaces,"░░░╚═╝░░░░╚════╝░░░░╚═╝░░░╚═╝░░╚══════╝╚═╝░░╚═╝╚═════╝░",toHighLight)
	spacePrint(precedingSpaces,"                               by Alex Van Mechelen    ",toHighLight)
end

function printPlayButton(highLight::Bool=false,precedingSpaces::Int=Int(floor((width-15)/2)))
	toHighLight = ()
	highLight && (toHighLight = ('█','▀','▄'))
	spacePrint(precedingSpaces,"█▀█ █░░ ▄▀█░█▄█",toHighLight)
	spacePrint(precedingSpaces,"█▀▀░█▄▄░█▀█ ░█░",toHighLight)
end

function printLoadButton(highLight::Bool=false,precedingSpaces::Int=Int(floor((width-15)/2)))
	toHighLight = ()
	highLight && (toHighLight = ('█','▀','▄'))
	spacePrint(precedingSpaces,"█░░ █▀█░▄▀█ █▀▄",toHighLight)
	spacePrint(precedingSpaces,"█▄▄░█▄█ █▀█░█▄▀",toHighLight)
end

function printExitButton(highLight::Bool=false,precedingSpaces::Int=Int(floor((width-13)/2)))
	toHighLight = ()
	highLight && (toHighLight = ('█','▀','▄'))
	spacePrint(precedingSpaces,"█▀▀░▀▄▀ █░▀█▀",toHighLight)
	spacePrint(precedingSpaces,"██▄ █░█░█ ░█░",toHighLight)
end

function printMainScreen(selectedButton::Int=1,subError::Bool=false,doubleLine::Bool=false)
	emptyLines(1)
	printStr = "^^^ Scroll up to see history ^^^"
	spaces(width-length(printStr)-1)
	printstyled(printStr,color=COL_HIGHLIGHT)
	emptyLines(Int(round(height-20)/2)-2)
	printLogo()
	emptyLines(4)
	printPlayButton(selectedButton == 1)
	emptyLines(2)
	printLoadButton(selectedButton == 2)
	emptyLines(2)
	printExitButton(selectedButton == 3)
	emptyLines(Int(round((height-22)/2)-Int(subError)-2*Int(doubleLine)))
end

# Main function, calls itself until user asks for exit
function mainScreen(selectedButton::Int=1,pastError::Bool=false)
	Calibrate()
	pastError || printMainScreen(selectedButton,selectedButton==1,selectedButton==1)

	# Print UI instructions if PLAY Button selected (and no past error)
	if selectedButton == 1 && !pastError
		printstyled("Controls: [Up] $(KEY_UP), [Down] $(KEY_DOWN), [Left] $(KEY_LEFT), [Right] $(KEY_RIGHT) followed by [Enter].\nAlso: [Pencil-Marks] $(KEY_PENCIL), [Toggle Error Display] $(KEY_TOGGLE_ERRORS) and [Ctrl Mode] $(KEY_CTRL).\nIn CTRL-Mode: [Solve Game] (solve), [Solve Game Step By Step] (solvestep) and [Exit Game] (exit).\n",color=COL_HIGHLIGHT)
	end

	pastError = false
	numButtons = 3
	try # Leave commented for debugging
		input = readline() #Gives error if unexpected character
		Calibrate() # Screen size could have changed while prompting user
		if isempty(input)
			if selectedButton == 1
				printMainScreen(selectedButton,true)
				printstyled("Please enter a Game-size [",color = COL_QUESTION)
				for i in 3:MAX_GEN_DIM
					printstyled(i,color = COL_QUESTION)
					if i != MAX_GEN_DIM
						printstyled(",",color = COL_QUESTION)
					else
						printstyled("].\n",color = COL_QUESTION)
					end
				end
				try # Leave commented for debugging
					inputSize = parse(Int,readline()) # Results in an error if not a number
					Calibrate() # Screen size could have changed while prompting user
					if 3 <= inputSize <= MAX_GEN_DIM
						easyMode = true # Standard difficulty is [Easy]
						if inputSize > 4
							printMainScreen(selectedButton,true)
							printstyled("Please enter a Game-difficulty [Easy (e) Hard (h)].\n",color = COL_QUESTION)
							inputDiff = lowercase(readline())
							Calibrate() # Screen size could have changed while prompting user

							if typeof(inputDiff) == String && inputDiff in ("e","h")
								easyMode = (inputDiff == "e")
							else
								mainError("Invalid Game-difficulty",selectedButton)
							end
						end

						loadingGame()
						game,solvedGame = generateGame(inputSize,easyMode)

						exitGame = loopGetInput(game,solvedGame)
						exitGame && return # Exit mainScreen (main function) if exit command given
						selectedButton = 1 # Reset to first button after playing.
					else
						mainError("Only Game-sizes between 3 and $(MAX_GEN_DIM) allowed",selectedButton)
					end
				catch # Leave commented for debugging
					mainError("Invalid Game-size",selectedButton)
				end # Leave commented for debugging
			elseif selectedButton == 2
				printMainScreen(selectedButton,true)
				printstyled("Please enter a Game-ID. [From: https://www.chiark.greenend.org.uk/~sgtatham/puzzles/js/towers.html]\n",color = COL_QUESTION)
				succes,game = loadGame(readline())
				Calibrate() # Screen size could have changed while prompting user
				if succes
					exitGame = loopGetInput(game)
					selectedButton = 1 # Reset to first button after playing.
					exitGame && return # Exit mainScreen (main function) if exit command given
				else
					mainError("Invalid Game-ID",selectedButton)
				end
			elseif selectedButton == 3
				exitScreen()
				print("\n")
				closeGame(t::Timer) = exit() # Close REPL
				wait(Timer(closeGame,3))
			else
				mainError("Something went wrong, please restart the game",selectedButton)
			end
		else
			length(input) > 1 && mainKeysError(selectedButton)
			input = collect(input)[1] # Convert to char
			if input in KEY_DOWN
				selectedButton = Int(selectedButton%numButtons+1) # Loop forwards through Menu buttons
			elseif input in KEY_UP
				selectedButton = Int((selectedButton-numButtons-1)%numButtons+numButtons) # Loop backwards through Menu buttons
			elseif input in KEY_TOGGLE_ERRORS
				global ShowOutBoardErrors = !ShowOutBoardErrors
				printMainScreen(selectedButton,true)
				outBoardErrorToggleMessage(ShowOutBoardErrors)
				mainScreen(selectedButton,true)
			elseif input in KEY_CTRL
				printMainScreen(selectedButton,true)
				printstyled("You have entered CTRL mode: Please enter command.\n",color = COL_QUESTION)
				try
					line = readline()
					Calibrate() # Screen size could have changed while prompting user
					lowercase(line) == "exit" && return
					if lowercase(line) == "solve"
						printMainScreen(selectedButton,true)
						printstyled("You cannot use this command in Main Menu.\n",color = COL_ERROR)
						pastError = true
					else
						eval(Meta.parse(line))
					end
				catch ex
					printMainScreen(selectedButton,true)
					printstyled("Something went wrong trying to run that command: $(ex).\n",color = COL_ERROR)
					pastError = true
				end
			else
				mainKeysError(selectedButton)
			end
		end
		mainScreen(selectedButton,pastError)
	catch # Leave commented for debugging
		mainKeysError(selectedButton)
	end # Leave commented for debugging
end

function mainKeysError(selectedButton::Int)
	mainError("Please only use [Up] $(KEY_UP) or [Down] $(KEY_DOWN) keys or [Enter] in Main Menu",selectedButton)
end

function mainError(str::String,selectedButton::Int=1)
	printMainScreen(selectedButton,true)
	printstyled(str*".\n",color = COL_ERROR)
	mainScreen(selectedButton,true)
end

function exitScreen()
	printMainScreen(1,true) # Selectedbutton is first button
	printstyled("Thank you for playing Towers by Alex Van Mechelen.",color = COL_HIGHLIGHT)
end

Base.active_repl.interface.modes[1].prompt="Towers> " #Changed REPL prompt text

include("Solver.jl")

if AUTO_START_GAME
	emptyLines(10000) # Make it impossible to scroll down
	Calibrate()
	loadingGame()
	mainScreen() # Automatically start on main screen when included in REPL (if AUTO_START_GAME = true)
	exitScreen()
else
	#This message is displayed when this file is included in the REPL and AUTO_START_GAME is turned off
	printstyled("Succesfully included Towers.jl by Alex Van Mechelen.",color = COL_HIGHLIGHT)
end

max_stars_repo_path: Towers Game/Towers.jl

max_issues_repo_path: Towers Game/Towers.jl

max_forks_repo_path: Towers Game/Towers.jl


text:
using Distributions, JLD, CSV, DataFrames,Pkg, DataStructures

austin_stations_pre = CSV.File("raw/austin-stations.csv") |> DataFrame

# Yeesian has his coords in 3-stations.csv in this string format:
# "GeoInterface.Point([-77.00514742400335,38.83089499410552])"
formatted_coordinates = [string("GeoInterface.Point([", austin_stations_pre["LONGITUDE"][i], ", ", austin_stations_pre["LATITUDE"][i], "])") for i in 1:size(austin_stations_pre["LONGITUDE"], 1)]

austin_stations_post = DataFrame(geometry = formatted_coordinates,
                                 OBJECTID_1 = austin_stations_pre["OBJECTID_1"]
                                 )

CSV.write("formatted/austin_stations_formatted.csv", austin_stations_post)

max_stars_repo_name: michaelhilborn/AmbulanceDeployement

max_issues_repo_name: joshua-ong/AmbulanceDeployment

max_forks_repo_name: joshua-ong/AmbulanceDeployment


text:
using NCDatasets, Revise, Statistics, PyPlot, DelimitedFiles
using tele
using Dates, DSP, DataFrames

data_dir = "/home/brynn/Documents/JP/860/tsa_data/"
oc2k_p = data_dir * "Ocean2kLR2015/input_data/Ocean2kLR2015sst.csv"
cesm_p1 = data_dir * "tos_Omon_CCSM4_past1000_r1i1p1_085001-134912.nc"
cesm_p2 = data_dir * "tos_Omon_CCSM4_past1000_r1i1p1_135001-185012.nc"

#READ OCEAN2K Data
mat = readdlm(oc2k_p, ',')
c1 = 109 #age column index Shevenell TS - got this from excel sheetlen
c2 = 110
age = mat[begin:20, c1]
d_age = [DateTime(x) for x in age]
val = mat[begin:20, c2]
savetype = ".pdf"

#READ NETCDFS TO BOUNDING BOX
#Shevenell = -64.87, -64.2
bbox = [71, 55, 290, 350]

(latc, lonc, timec1, datac1) = read_netcdf(cesm_p1, "tos", bbox)
(latc, lonc, timec, datac) = read_netcdf(cesm_p2, "tos", bbox)
#concatenate two sets of data
datac = cat(datac1, datac; dims = 3)
timec = cat(timec1, timec; dims = 1)
datac = datac .- 273.15 #convert to celsius
anomc = anomaly(datac, timec)

#Plot data - mostly just to make sure we're in the right spot
contplot(latc,lonc,temp_mean(datac), "images/oc2k"*savetype)

#FILTER CHARACTERISTICS
responsetype = Lowpass(0.1*10^(-7),fs = 3.80265176 * 10^(-7))
designmethod = Butterworth(4)

#Compute spatial average of CESM data
timec = reinterpret(DateTime, timec) #convert from DatetimeNoLeap to DateTime
tavgc_raw = spatial_mean(datac)

#Convert to annual average
tavgc_raw_a, years = yearly_mean(timec, tavgc_raw)
years_d = [DateTime(x) for x in years]

#Plot annual average CESM and Shevenell dataset
figure()
plot(years_d, tavgc_raw_a, alpha = 0.5)
plot(d_age, val, ".-")
legend(["CESM Annual Average", "Shevenell et al"])
xlabel("Time (years)")
ylabel("SST [°C]")
title("Alkenone Reconstruction of SST vs CESM Annual Average")
savefig("images/oc2k_cesm"*savetype)

#compute statistics
df = DataFrame(min = [], max = [], mean = [], median = [], std = [])
push!(df, character(tavgc_raw))
push!(df, character(val))

#Make boxplot of statistics
figure(figsize = (6, 8))
boxplot([val, tavgc_raw], notch = true, labels =  ["Shevenell et al", "CESM (850 - 1850)"])
ylabel("SST [°C]")
grid(true)
savefig("images/boxes_shev"*savetype)
close()

#Extreme value analysis
quant = quantile(tavgc_raw, [0.25, 0.5, 0.75])
thresh = (quant[3] - quant[1])* 1.5 + quant[3]
ind = findall(>(thresh), tavgc_raw)

extreme1 = temp_mean(datac[:,:, ind])
extreme2 = temp_mean(anomc[:,:, ind])

#Plot extreme value anomalies
contplot(latc,lonc,extreme2, "images/ext_anom_lm"*savetype)


text:
using RemoteREPL
using Test
using Sockets
using RemoteREPL: repl_prompt_text, DEFAULT_PORT

ENV["JULIA_DEBUG"] = "RemoteREPL"

@testset "Protocol header handshake" begin
    # Happy case
    io = IOBuffer()
    RemoteREPL.send_header(io)
    seek(io, 0)
    @test RemoteREPL.verify_header(io)

    # Broken magic number
    io = IOBuffer()
    write(io, "BrokenMagic", RemoteREPL.PROTOCOL_VERSION)
    seek(io, 0)
    @test_throws ErrorException RemoteREPL.verify_header(io)
    # Version mismatch
    io = IOBuffer()
    write(io, RemoteREPL.PROTOCOL_VERSION, typemax(UInt32))
    seek(io, 0)
    @test_throws ErrorException RemoteREPL.verify_header(io)

    # ser_version=10 in julia 1.5
    for (local_ver, remote_ver) in [(10,13), (13,10)]
        io = IOBuffer()
        seek(io, 0)
        RemoteREPL.send_header(io, remote_ver)
        @test_throws ErrorException RemoteREPL.verify_header(io, local_ver)
    end
end

@testset "Prompt text" begin
    function fake_conn(host, port; is_open=true)
        io = IOBuffer()
        is_open || close(io)
        RemoteREPL.Connection(host, port, nothing, nothing, nothing, nothing, io, :Main)
    end
    @test repl_prompt_text(fake_conn(Sockets.localhost, DEFAULT_PORT)) == "julia@localhost> "
    @test repl_prompt_text(fake_conn("localhost",       DEFAULT_PORT)) == "julia@localhost> "
    @test repl_prompt_text(fake_conn(ip"192.168.1.1",   DEFAULT_PORT)) == "julia@192.168.1.1> "

    @test repl_prompt_text(fake_conn("ABC", DEFAULT_PORT)) == "julia@ABC> "
    @test repl_prompt_text(fake_conn("ABC", 12345))        == "julia@ABC:12345> "
    @test repl_prompt_text(fake_conn("ABC", DEFAULT_PORT, is_open=false)) == "julia@ABC [disconnected]> "
end

# Connect to a non-default loopback address to test SSH integration
test_interface = ip"127.111.111.111"

# Detect whether an ssh server is running on the default port 22
use_ssh = if "use_ssh=true" in ARGS
    true
elseif "use_ssh=false" in ARGS
    false
else
    # Autodetct
    try
        socket = Sockets.connect(test_interface, 22)
        # https://tools.ietf.org/html/rfc4253#section-4.2
        id_string = String(readavailable(socket))
        startswith(id_string, "SSH-")
    catch
        false
    end
end

if !use_ssh
    test_interface = Sockets.localhost
end
@info use_ssh ? "Running tests with SSH tunnel" : "Testing without SSH tunnel - localhost only"

# Use non-default port to avoid clashes with concurrent interactive use or testing.
test_port = RemoteREPL.find_free_port(Sockets.localhost)
server_proc = run(`$(Base.julia_cmd()) -e "using Sockets; using RemoteREPL; serve_repl($test_port)"`, wait=false)

try

@testset "RemoteREPL.jl" begin
    local conn = nothing
    max_tries = 4
    for i=1:max_tries
        try
            conn = RemoteREPL.Connection(host=test_interface, port=test_port,
                                         tunnel=use_ssh ? :ssh : :none,
                                         ssh_opts=`-o StrictHostKeyChecking=no`)
            break
        catch exc
            if i == max_tries
                rethrown()
            end
            # Server not yet started - continue waiting
            sleep(2)
        end
    end
    @assert isopen(conn)

    # Some basic tests of the transport and server side and partial client side.
    #
    # More full testing of the client code would requires some tricky mocking
    # of the REPL environment.
    function runcommand(cmdstr)
        result = RemoteREPL.run_remote_repl_command(conn, IOBuffer(), cmdstr)
        # Unwrap Text for testing purposes
        return result isa Text ? result.content : result
    end

    @test runcommand("asdf = 42") == "42"
    @test runcommand("Main.asdf") == "42"
    @test !isdefined(Main, :asdf) # asdf not defined locally

    # Output Limiting
    @test '⋮' in runcommand("ones(1000)")

    # Error formatting
    @test occursin(r"DivideError.*Stacktrace"s, runcommand("1÷0"))

    # Logging
    @test occursin(r"Info:.*xxx"s, runcommand("""@info "hi" xxx=[1,2]"""))

    # Semicolon suppresses output
    @test isnothing(runcommand("asdf;"))

    # Help mode
    @test occursin("helpmodetest documentation!",
        begin
            runcommand("function helpmodetest end")
            runcommand("@doc \"helpmodetest documentation!\" helpmodetest")
            runcommand("?helpmodetest")
        end)

    # Completions
    runcommand("complete_me_1 = 1")
    runcommand("complete_me_2 = 2")
    completion_msg = RemoteREPL.send_and_receive(conn,
                        (:repl_completion, ("complete_m", "complete_m")))
    @test completion_msg[1] == :completion_result
    @test completion_msg[2] == (["complete_me_1", "complete_me_2"], "complete_m", true)

    # Evaluation in other modules
    runcommand("""module TestMod
                  struct SomeStruct
                  end
                  var_in_test_mod = 123
               end""")
    @test runcommand("%module TestMod") == "Evaluating commands in module Main.TestMod"
    @test runcommand("var_in_test_mod") == "123"
    # Test that show() on the remote side uses the eval module as the context
    # module in the show IOContext
    @test runcommand("SomeStruct") == "SomeStruct"
    completion_msg = RemoteREPL.send_and_receive(conn,
                        (:repl_completion, ("complete_m", "complete_m")))
    @test completion_msg[1] == :completion_result
    @test completion_msg[2] == ([], "complete_m", true)

    # Test the @remote macro
    Main.eval(:(clientside_var = 0:41))
    @test runcommand("serverside_var = 1 .+ @remote clientside_var") == "1:42"
    @test Main.clientside_var == 0:41
    @test @remote(conn, serverside_var) == 1:42
    @test_throws RemoteREPL.RemoteException @remote(conn, error("hi"))

    # Special case handling of stdout
    @test runcommand("println(@remote(stdout), \"hi\")") == "hi\n"

    # Issue #26: world age errors showing output
    @test runcommand("""
        struct A
        end

        Base.show(io::IO, a::A) = print(io, "Surprise")

        A()
        """) == "Surprise"
    @test occursin("Surprise", runcommand("""
        struct MyExc <: Exception
        end

        Base.showerror(io::IO, e::MyExc) = print(io, "Surprise")

        throw(MyExc())
        """))

    # Execute a single command on a separate connection
    @test (RemoteREPL.remote_eval(test_interface, test_port, "asdf")::Text).content == "42"

    # Test interrupts
    # 1. The straightforward reliable way, based on directly sending an
    # interrupt message to the server
    @sync begin
        # Interrupt after 0.5 seconds
        @async begin
            sleep(0.5)
            RemoteREPL.send_interrupt(conn)
        end
        # Remote command which attempts to sleep for 10 seconds and returns the
        # actual time spent sleeping
        @test @remote(conn, (@timed try sleep(10) ; catch ; end).time) < 1
    end
    # 2. Approximately simulate the way that a user-driven ^C gets delivered to
    # the task waiting on the socket
    @sync begin
        task = current_task()
        @async begin
            sleep(0.5)
            schedule(task, InterruptException(), error=true)
        end
        # Remote command which attempts to sleep for 10 seconds and returns the
        # actual time spent sleeping
        @test @remote(conn, (@timed try sleep(10) ; catch ; end).time) < 1
    end
    # 3. Test that we disconnect from the server when trying to interrupt
    #    non-interruptable code
    @sync begin
        task = @async try
            # Execute some non-interruptable code on the server.
            #
            # Normally this would be a compute-heavy task which doesn't yield,
            # but a non-interruptable sleep which catches InterruptException is
            # a more controlled way to simulate this.
            @remote(conn, begin
                t0 = time()
                sleep_time = 5
                while true
                    dt = sleep_time + t0 - time()
                    if dt < 0
                        break
                    end
                    try
                        sleep(dt)
                    catch exc
                        if exc isa InterruptException
                            continue
                        else
                            rethrow()
                        end
                    end
                end
            end)
        catch exc
            exc
        end
        # Three attempts to interrupt will result in the connection being
        # closed to allow the client to recover.
        sleep(0.5)
        schedule(task, InterruptException(), error=true)
        @test !istaskdone(task)
        sleep(0.5)
        schedule(task, InterruptException(), error=true)
        @test !istaskdone(task)
        sleep(0.5)
        schedule(task, InterruptException(), error=true)
        @test fetch(task) == RemoteREPL.RemoteException("ERROR: Failed to interrupt server, connection closed!")
    end
end

finally
    kill(server_proc)
end