text:
struct OutputFlags
    write_volt_maps::Bool
    write_cur_maps::Bool
    write_cum_cur_map_only::Bool
    write_max_cur_maps::Bool
    set_null_currents_to_nodata::Bool
    set_null_voltages_to_nodata::Bool
    compress_grids::Bool
    log_transform_maps::Bool
end

function compute_3col(resistances::Matrix{T}) where {T}
    fp = deleteat!(resistances[:,1], 1)
    l = length(fp)
    r3col = zeros(T, div(l * (l-1), 2), 3)
    k = 1
    for i = 1:l
        for j = i+1:l
            r3col[k,1] = fp[i]
            r3col[k,2] = fp[j]
            r3col[k,3] = resistances[j+1,i+1]
            k += 1
        end
    end
    r3col
end

function write_cur_maps(name, output, component_data, finitegrounds, flags, cfg)

    # Get desired data
    G = component_data.matrix
    voltages = flags.is_advanced ? output : output.voltages
    cc = component_data.cc
    nodemap = component_data.local_nodemap
    hbmeta = component_data.hbmeta
    cellmap = component_data.cellmap

    # Flags
    log_transform = flags.outputflags.log_transform_maps
    set_null_currents_to_nodata =
        flags.outputflags.set_null_currents_to_nodata
    write_max_cur_maps = flags.outputflags.write_max_cur_maps
    write_cum_cur_map_only = flags.outputflags.write_cum_cur_map_only

    node_currents, branch_currents = _create_current_maps(G, voltages, finitegrounds, cfg, nodemap = nodemap, hbmeta = hbmeta)

    if !flags.is_raster

        # Branch currents
        branch_currents_array = _convert_to_3col(branch_currents, cc)

        # TODO: implement cumulative maps for netowrk mode

        # Accumulate branch currents
        # cum_branch_curr[mycsid()] .+= branch_currents_array[:,3]

        # Node currents
        node_currents_array = _append_name_to_node_currents(node_currents, cc)

        # Accumulate node currents
        # cum_node_curr[mycsid()] .+= node_currents_array[:,2]

        # !write_cum_cur_map_only &&
        write_currents(node_currents_array, branch_currents_array, name, cfg)
    else

        cmap = node_currents
        cum_curr = output.cum.cum_curr
        max_curr = output.cum.max_curr

        # Process the current map
        process_grid!(cmap, cellmap, hbmeta, log_transform = log_transform,
                            set_null_to_nodata = set_null_currents_to_nodata)

        # Accumulate by default
        cum_curr[mycsid()] .+= cmap

        # Max current if user asks for it
        if write_max_cur_maps
            max_curr[mycsid()] .= max.(max_curr[mycsid()], cmap)
        end

        # Write current maps
        !write_cum_cur_map_only &&
                        write_grid(cmap, name, cfg, hbmeta)

   end
end

function write_currents(node_curr_arr, branch_curr_arr, name, cfg)
   pref = split(cfg["output_file"], ".out")[1]
   writedlm("$(pref)_node_currents$(name).txt", node_curr_arr, '\t')
   writedlm("$(pref)_branch_currents$(name).txt", branch_curr_arr, '\t')
end

_append_name_to_node_currents(node_currents, cc) = [cc node_currents]

function _convert_to_3col(branch_currents, cc)

    l = length(branch_currents.nzval)
    graph = zeros(l, 3)
    graph[:,3] = branch_currents.nzval

    # Inspired by show method for sparse matrices in Julia Base
    b = branch_currents

    k = 1
    for i = 1:size(b, 1)
        for j in nzrange(b, i)
            row = b.rowval[j]
            graph[k,1] = cc[row]
            graph[k,2] = cc[i]
            k += 1
        end
    end

    graph
end

function _create_current_maps(G, voltages, finitegrounds, cfg; nodemap = Matrix{Float64}(), hbmeta = RasterMeta())

    node_currents = get_node_currents(G, voltages, finitegrounds)

    if cfg["data_type"] == "network"

        branch_currents = _get_branch_currents(G, voltages, true)
        branch_currents = abs.(branch_currents)
        return node_currents, branch_currents

    else

        current_map = zeros(eltype(G), hbmeta.nrows, hbmeta.ncols)
        for j = 1:size(nodemap, 2)
            for i = 1:size(nodemap, 1)
                idx = nodemap[i,j]
                if idx == 0
                    continue
                else
                    current_map[i,j] = node_currents[idx]
                end
            end
        end

        return current_map, spzeros(0,0)
    end
end

function get_node_currents(G, voltages, finitegrounds)

    node_currents_pos = _get_node_currents_posneg(G, voltages, finitegrounds, true)
    node_currents_neg = _get_node_currents_posneg(G, voltages, finitegrounds, false)
    node_currents = map((x,y) -> x > y ? x : y, node_currents_pos, node_currents_neg)

end

function _get_node_currents_posneg(G::SparseMatrixCSC{T,V},
                            voltages, finitegrounds, pos) where {T,V}

    branch_currents = _get_branch_currents(G, voltages, pos)
    branch_currents = branch_currents - branch_currents'
    dropnonzeros!(branch_currents)

	if finitegrounds[1]!= -9999
        finiteground_currents = finitegrounds .* voltages
        if pos
            map!(x -> x < 0 ? -x : 0, finiteground_currents, finiteground_currents)
        else
            map!(x -> x > 0 ? x : 0, finiteground_currents, finiteground_currents)
        end
        n = size(G, 1)
        branch_currents = branch_currents + spdiagm(0 => finiteground_currents)
    end

    s = vec(sum(branch_currents, dims=1))

    s
end

function dropnonzeros!(G)
    for i = 1:size(G, 1)
        for j in nzrange(G, i)
            row = G.rowval[j]
            val = G.nzval[j]
            if val < 0
                G.nzval[j] = 0
            end
        end
    end
    dropzeros!(G)
end

function _get_branch_currents(G::SparseMatrixCSC{T,V}, voltages, pos) where {T,V}

    branch_currents = _get_branch_currents_posneg(G, voltages, pos)

    # Make sparse matrix with branch_currents as right upper triangle
    N = size(G, 1)
    n = size(branch_currents, 1)
    I = zeros(V, n)
    J = zeros(V, n)
    k = 1
    for i = 1:N
        for j in nzrange(G, i)
            row = G.rowval[j]
            if i > row
                I[k] = row
                J[k] = i
                k += 1
            end
        end
    end
    @assert n + 1 == k

    B = sparse(I, J, branch_currents, N, N)

    B
	# branch_currents
end

function _get_branch_currents_posneg(G::SparseMatrixCSC{T,V},
                                v::Vector{T}, pos) where {T,V}

    n = count_upper(G)
    b = zeros(T, n)
    k = 1
    if pos
        for i = 1:size(G, 1)
            for j in nzrange(G, i)
                row = G.rowval[j]
                val = G.nzval[j]
                if i > row
                    b[k] = abs(val) * (v[row] - v[i])
                    k += 1
                end
            end
        end
    else
        for i = 1:size(G, 1)
            for j in nzrange(G, i)
                row = G.rowval[j]
                val = G.nzval[j]
                if i > row
                    b[k] = abs(val) * (v[i] - v[row])
                    k += 1
                end
            end
        end
    end
    # @show n, k
    @assert n + 1 == k
    maxcur = maximum(b)
    # map!(x -> abs(x / maxcur) < 1e-8 ? 0 : x, b, b)
    for i = 1:size(b, 1)
        if abs(b[i]/ maxcur) < 1e-8
            b[i] = 0
        end
    end

    b
end

function count_upper(G)
    n = 0
    for i = 1:size(G, 1)
        for j in nzrange(G,i)
            row = G.rowval[j]
            if i > row
                n += 1
            end
        end
    end
    n
end

function process_grid!(cmap, cellmap, hbmeta; log_transform = false,
                                set_null_to_nodata = false)
    if log_transform
        map!(x -> x > 0 ? log10(x) : float(hbmeta.nodata), cmap, cmap)
    end

    if set_null_to_nodata
        for i in eachindex(cmap)
            if cellmap[i] == 0
                cmap[i] = hbmeta.nodata
            end
        end
    end

end

function write_grid(cmap, name, cfg, hbmeta;
                        voltage = false, cum = false,
                        max = false)

    str = "curmap"
    if cum
        str = "cum_$(str)"
    elseif max
        str = "max_$(str)"
    elseif voltage
        str = "voltmap"
    end

    pref = split(cfg["output_file"], ".out")[1]
    filename = "$(pref)_$(str)$(name)"

    cfg["write_as_tif"] in TRUELIST ? (file_format = "tif") :
            (file_format = "asc")

    write_raster(filename,
                 cmap,
                 hbmeta.wkt,
                 hbmeta.transform,
                 file_format)
end

function write_grid(cmap, name, cfg, hbmeta, cellmap;
                        voltage = false, cum = false, max = false,
                        log_transform = false, set_null_to_nodata = false)

    pref = split(cfg["output_file"], ".out")[1]

    if log_transform
        map!(x -> x > 0 ? log10(x) : float(hbmeta.nodata), cmap, cmap)
    end

    if set_null_to_nodata
        for i in eachindex(cmap)
            if cellmap[i] == 0
                cmap[i] = hbmeta.nodata
            end
        end
    end

    str = "curmap"
    if cum
        str = "cum_$(str)"
    elseif max
        str = "max_$(str)"
    elseif voltage
        str = "voltmap"
    end

    filename = "$(pref)_$(str)$(name)"

    cfg["write_as_tif"] in TRUELIST ? (file_format = "tif") :
            (file_format = "asc")

    write_raster(filename,
                 cmap,
                 hbmeta.wkt,
                 hbmeta.transform,
                 file_format)
end

function write_volt_maps(name, output, component_data, flags, cfg)

    voltages = flags.is_advanced ? output : output.voltages

    if !flags.is_raster

        cc = component_data.cc
        write_voltages(cfg["output_file"], name, voltages, cc)

    else

        # Desired data
        # voltages = output.voltages
        cc = component_data.cc
        hbmeta = component_data.hbmeta
        nodemap = component_data.local_nodemap
        set_null_voltages_to_nodata = flags.outputflags.set_null_voltages_to_nodata

        vm = _create_voltage_map(voltages, nodemap, hbmeta)
        write_grid(vm, name, cfg, hbmeta, component_data.cellmap, voltage = true,
                        set_null_to_nodata = set_null_voltages_to_nodata)
    end
end

function write_voltages(output, name, voltages::Vector{T}, cc) where {T}

    volt_arr = zeros(T, size(voltages, 1), 2)
    volt_arr[:,1] = cc
    volt_arr[:,2] = voltages

    pref = split(output, ".out")[1]
    writedlm("$(pref)_voltages$(name).txt", volt_arr)

end

function _create_voltage_map(voltages::Vector{T}, nodemap, hbmeta) where {T}
    voltmap = zeros(T, hbmeta.nrows, hbmeta.ncols)
    for j = 1:size(nodemap, 2)
        for i = 1:size(nodemap, 1)
            idx = nodemap[i,j]
            if idx == 0
                continue
            else
                voltmap[i,j] = voltages[idx]
            end
        end
    end
    voltmap
end

alloc_map(hbmeta) = zeros(hbmeta.nrows, hbmeta.ncols)

function accum_voltages!(base, newvolt, nodemap, hbmeta)
    voltmap = _create_voltage_map(newvolt, nodemap, hbmeta)
    for i in eachindex(base)
        base[i] += voltmap[i]
    end
end

function accum_currents!(base, newcurr, cfg, G, voltages, finitegrounds, nodemap, hbmeta)
    node_currents, branch_currents = _create_current_maps(G, voltages, finitegrounds, cfg,
                                                            nodemap = nodemap, hbmeta = hbmeta)

    for i in eachindex(base)
        base[i] += node_currents[i]
    end
end

function save_resistances(r, cfg)
    pref = split(cfg["output_file"], ".out")[1]
    filename = "$(pref)_resistances.out"
    filename_3col = "$(pref)_resistances_3columns.out"
    rcol = compute_3col(r)
    open(filename, "w") do f
        writedlm(f, r, ' ')
    end
    open(filename_3col, "w") do f
        writedlm(f, rcol, ' ')
    end
end

function write_cum_maps(cum, cellmap::Matrix{T}, cfg, hbmeta, write_max, write_cum) where T

    if write_cum || cfg["write_cur_maps"] in TRUELIST
        cum_curr = zeros(T, size(cellmap)...)
        for i = 1:nprocs()
            cum_curr .+= cum.cum_curr[i]
        end
        postprocess_cum_curmap!(cum_curr)
        write_grid(cum_curr, "", cfg, hbmeta, cum = true)
    end

    if write_max
        max_curr = fill(T(-9999), size(cellmap)...)
        for i = 1:nprocs()
            max_curr .= max.(cum.max_curr[i], max_curr)
        end
        postprocess_cum_curmap!(max_curr)
        write_grid(max_curr, "", cfg, hbmeta, max = true)
    end

end

# Write a single band raster, either in .tif or .asc format,
# inspired by GeoArrays.write()
function write_raster(fn_prefix::String,
                      array::Matrix{T} where T <: Number,
                      wkt::String,
                      transform,
                      file_format::String)
    # transponse array back to columns by rows
    array_t = permutedims(array, [2, 1])

    width, height = size(array_t)

    # Define extension and driver based in file_format
    file_format == "tif" ? (ext = ".tif"; driver = "GTiff") :
            (ext = ".asc"; driver = "AAIGrid")

    file_format == "tif" ? (options = ["COMPRESS=LZW"]) :
                           (options = [])

    # Append file extention to filename
    fn = string(fn_prefix, ext)

    # Create raster in memory *NEEDED* because no create driver for .asc
    ArchGDAL.create(fn_prefix,
                    driver = ArchGDAL.getdriver("MEM"),
                    width = width,
                    height = height,
                    nbands = 1,
                    dtype = eltype(array_t),
                    options = options) do dataset
        band = ArchGDAL.getband(dataset, 1)
        # Write data to band
        ArchGDAL.write!(band, array_t)

        # Write nodata and projection info
        ArchGDAL.setnodatavalue!(band, -9999.0)
        ArchGDAL.setgeotransform!(dataset, transform)
        ArchGDAL.setproj!(dataset, wkt)

        # Copy memory object to disk (necessary because ArchGDAL.create
        # does not support creation of ASCII rasters)
        ArchGDAL.write(dataset, fn,
                       driver = ArchGDAL.getdriver(driver),
                       options = options)
    end

end


text:
##############################################################################
##
## Least Squares using LSMR
##
##############################################################################

##############################################################################
##
## FixedEffectVector : vector x in A'Ax = A'b
##
## We define these methods used in lsmr! (duck typing):
## copy!, fill!, scale!, axpy!, norm
##
##############################################################################

type FixedEffectVector
    _::Vector{Vector{Float64}}
end

function FixedEffectVector(fes::Vector{FixedEffect})
    out = Vector{Float64}[]
    for fe in fes
        push!(out, similar(fe.scale))
    end
    return FixedEffectVector(out)
end

eltype(fem::FixedEffectVector) = Float64

length(fev::FixedEffectVector) = reduce(+, map(length, fev._))

function copy!(fev2::FixedEffectVector, fev1::FixedEffectVector)
    for i in 1:length(fev1._)
        copy!(fev2._[i], fev1._[i])
    end
    return fev2
end

function fill!(fev::FixedEffectVector, x)
    for i in 1:length(fev._)
        fill!(fev._[i], x)
    end
end

function scale!(fev::FixedEffectVector, α::Number)
    for i in 1:length(fev._)
        scale!(fev._[i], α)
    end
    return fev
end

function axpy!(α::Number, fev1::FixedEffectVector, fev2::FixedEffectVector)
    for i in 1:length(fev1._)
        axpy!(α, fev1._[i], fev2._[i])
    end
    return fev2
end

function norm(fev::FixedEffectVector)
    out = zero(Float64)
    for i in 1:length(fev._)
        out += sumabs2(fev._[i])
    end
    return sqrt(out)
end

##############################################################################
##
## FixedEffectMatrix
##
## A is the model matrix of categorical variables
## normalized by diag(1/a1, ..., 1/aN) (Jacobi preconditoner)
##
## We define these methods used in lsmr! (duck typing):
## A_mul_B!(α, A, b, β, c) updates c -> α Ab + βc
## Ac_mul_B!(α, A, b, β, c) updates c -> α A'b + βc
##
##############################################################################

type FixedEffectMatrix
    _::Vector{FixedEffect}
    m::Int
    n::Int
    cache::Vector{Vector{Float64}}
end

function FixedEffectMatrix(fev::Vector{FixedEffect})
    m = length(fev[1].refs)
    n = reduce(+, map(x -> length(x.scale),  fev))
    caches = Vector{Float64}[]
    for i in 1:length(fev)
        push!(caches, cache(fev[i]))
    end
    return FixedEffectMatrix(fev, m, n, caches)
end

function cache(fe::FixedEffect)
    out = zeros(Float64, length(fe.refs))
    @inbounds @simd for i in 1:length(out)
        out[i] = fe.scale[fe.refs[i]] * fe.interaction[i] * fe.sqrtw[i]
    end
    return out
end

eltype(fem::FixedEffectMatrix) = Float64

size(fem::FixedEffectMatrix, dim::Integer) = (dim == 1) ? fem.m :
                                            (dim == 2) ? fem.n : 1

# Define x -> A * x
function A_mul_B_helper!(α::Number, fe::FixedEffect,
                        x::Vector{Float64}, y::AbstractVector{Float64}, cache::Vector{Float64})
    @inbounds for (i, j) in zip(1:length(y), eachindex(y))
        y[j] += α * x[fe.refs[i]] * cache[i]
    end
end
function A_mul_B!(α::Number, fem::FixedEffectMatrix, fev::FixedEffectVector,
                β::Number, y::AbstractVector{Float64})
    safe_scale!(y, β)
    for i in 1:length(fev._)
        A_mul_B_helper!(α, fem._[i], fev._[i], y, fem.cache[i])
    end
    return y
end

# Define x -> A' * x
function Ac_mul_B_helper!(α::Number, fe::FixedEffect,
                        y::AbstractVector{Float64}, x::Vector{Float64}, cache::Vector{Float64})
    @inbounds for (i, j) in zip(1:length(y), eachindex(y))
        x[fe.refs[i]] += α * y[j] * cache[i]
    end
end
function Ac_mul_B!(α::Number, fem::FixedEffectMatrix,
                y::AbstractVector{Float64}, β::Number, fev::FixedEffectVector)
   safe_scale!(fev, β)
    for i in 1:length(fev._)
       # @code_warntype Ac_mul_B_helper!(α, fem._[i], y, fev._[i])
        Ac_mul_B_helper!(α, fem._[i], y, fev._[i], fem.cache[i])
    end
    return fev
end

function safe_scale!(x, β)
    if β != 1
        β == 0 ? fill!(x, zero(eltype(x))) : scale!(x, β)
    end
end

##############################################################################
##
## FixedEffectProblem is a wrapper around a FixedEffectMatrix
## with some storage arrays used when solving (A'A)X = A'y
##
##############################################################################

type LSMRFixedEffectProblem <: FixedEffectProblem
    m::FixedEffectMatrix
    x::FixedEffectVector
    v::FixedEffectVector
    h::FixedEffectVector
    hbar::FixedEffectVector
    u::Vector{Float64}
end

function FixedEffectProblem(fes::Vector{FixedEffect}, ::Type{Val{:lsmr}})
    m = FixedEffectMatrix(fes)
    x = FixedEffectVector(fes)
    v = FixedEffectVector(fes)
    h = FixedEffectVector(fes)
    hbar = FixedEffectVector(fes)
    u = Array(Float64, size(m, 1))
    return LSMRFixedEffectProblem(m, x, v, h, hbar, u)
end

get_fes(fep::LSMRFixedEffectProblem) = fep.m._

function solve!(fep::LSMRFixedEffectProblem, r::AbstractVector{Float64};
    tol::Real = 1e-8, maxiter::Integer = 100_000)
    fill!(fep.x, zero(Float64))
    copy!(fep.u, r)
    x, ch = lsmr!(fep.x, fep.m, fep.u, fep.v, fep.h, fep.hbar;
        atol = tol, btol = tol, conlim = 1e8, maxiter = maxiter)
    return div(ch.mvps, 2), ch.isconverged
end

function solve_residuals!(fep::LSMRFixedEffectProblem, r::AbstractVector{Float64}; kwargs...)
    iterations, converged = solve!(fep, r; kwargs...)
    A_mul_B!(-1.0, fep.m, fep.x, 1.0, r)
    return r, iterations, converged
end

function solve_coefficients!(fep::LSMRFixedEffectProblem, r::AbstractVector{Float64}; kwargs...)
    iterations, converged = solve!(fep, r; kwargs...)
    for i in 1:length(fep.x._)
        broadcast!(*, fep.x._[i], fep.x._[i], fep.m._[i].scale)
    end
    return fep.x._, iterations, converged
end

ext: jl


text:
# Autogenerated wrapper script for beefalo_jll for armv7l-linux-gnueabihf
export beefalo

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("beefalo")
JLLWrappers.@declare_executable_product(beefalo)
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_executable_product(
        beefalo,
        "bin/beefalo",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# # Steady-state flow around a cylinder in 2D using immersed boundaries

using Statistics
using Plots
using JLD2
using Printf

using Oceananigans
using Oceananigans.Advection
using Oceananigans.Grids
using Oceananigans.Fields
using Oceananigans.AbstractOperations
using Oceananigans.OutputWriters

# setting boundary condition topology
topology=(Periodic, Bounded, Bounded)

# setting up 2D grid
grid = RegularCartesianGrid(topology=topology, size=(350, 350, 1), x=(20, 40), y=(10, 30), z=(0, 1))
#grid = RegularCartesianGrid(topology=topology, size=(1500, 1500, 1), x=(0, 60), y=(0, 60), z=(0, 1))

# reynolds number
Re = 40

#cylinder with center at (30,20)
const R = 1 # radius
inside_cylinder(x, y, z) = ((x - 30)^2 + (y - 20)^2) <= R # immersed solid

# boundary conditions: inflow and outflow in y
v_bcs = VVelocityBoundaryConditions(grid,
                                    north = OpenBoundaryCondition(1),
                                    south = OpenBoundaryCondition(1))

# setting up incompressible model with immersed boundary
model = IncompressibleModel(timestepper = :RungeKutta3,
                              advection = UpwindBiasedFifthOrder(),
                                   grid = grid,
                               buoyancy = nothing,
                                tracers = nothing,
                                closure = IsotropicDiffusivity(ν=1/Re),
                    boundary_conditions = (v=v_bcs,),
                      immersed_boundary = inside_cylinder)

# initial condition
# setting velocitiy to zero inside the cylinder and 1 everywhere else
v₀(x, y, z) = ifelse(inside_cylinder(x,y,z),0.,1.)
set!(model, v=v₀)

progress(sim) = @info @sprintf("Iteration: % 4d, time: %.2f, max(v): %.5f, min(v): %.5f",
sim.model.clock.iteration,
sim.model.clock.time,
maximum(sim.model.velocities.v.data),
minimum(sim.model.velocities.v.data))

simulation = Simulation(model, Δt=5.7e-3, stop_time=5, iteration_interval=10, progress=progress)

# ## Output

simulation.output_writers[:fields] = JLD2OutputWriter(model,
                                                      merge(model.velocities, model.pressures),
                                                      schedule = TimeInterval(0.5),
                                                      prefix = "flow_around_cylinder",
                                                      force = true)

# run it
run!(simulation)

# Analyze Results
file = jldopen(simulation.output_writers[:fields].filepath)

iterations = parse.(Int, keys(file["timeseries/t"]))

u, v, w = model.velocities

xv, yv, zv = nodes(v)

# defining a function to mark boundary
function circle_shape(h, k, r)
            θ = LinRange(0,2*π,500)
            h.+r*sin.(θ),k.+r*cos.(θ)
end

@info "Making a neat movie of velocity..."

anim = @animate for (i, iteration) in enumerate(iterations)

    @info "Plotting frame $i from iteration $iteration..."

    t = file["timeseries/t/$iteration"]

    v_slice = file["timeseries/v/$iteration"][:, :, 1]
    @info maximum(v_slice) minimum(v_slice)
    v_max = maximum(abs, v_slice)
    v_lim = 0.8 * v_max

    v_levels = vcat([-v_max], range(-v_lim, stop=v_lim, length=50), [v_max])

    v_plot = contourf(xv, yv, v_slice';
                      linewidth = 0,
                          color = :balance,
                    aspectratio = 1,
                          title = @sprintf("v(x, y, t = %.1f) around a cylinder", t),
                         xlabel = "x",
                         ylabel = "y",
                         levels = v_levels,
                          xlims = (grid.xF[1], grid.xF[grid.Nx]),
                          ylims = (grid.yF[1], grid.yF[grid.Ny]),
                          clims = (-v_lim, v_lim))
    plot!(circle_shape(30,20,1),seriestype=[:shape,],linecolor=:black,
    legend=false,fillalpha=0, aspect_ratio=1)
end

gif(anim, "flow_around_cyl_velocity.gif", fps = 8) # hide


text:
for (A,F) in ((:(+), :mag_add), (:(-), :mag_sub), (:(*), :mag_mul), (:(/), :mag_div))
    @eval begin
        function ($A)(x::Mag, y::Mag)
            z = Mag()
            ccall(@libarb($F), Cvoid, (Ref{Mag}, Ref{Mag}, Ref{Mag}), z, x, y)
            return z
         end
    end
end

for (A,F) in ((:(+), :arf_add), (:(-), :arf_sub), (:(*), :arf_mul_rnd_any), (:(/), :arf_div))
    @eval begin
        function ($A)(x::ArbFloat{P}, y::ArbFloat{P}, roundingmode::RoundingMode=RoundNearest) where {P}
            z = ArbFloat{P}()
            rounding = match_rounding_mode(roundingmode)
            rounddir = ccall(@libarb($F), Cint, (Ref{ArbFloat}, Ref{ArbFloat}, Ref{ArbFloat}, Clong, Cint), z, x, y, P, rounding)
            return z
         end
    end
end

for (A,F) in ((:(+), :arb_add), (:(-), :arb_sub), (:(*), :arb_mul), (:(/), :arb_div), (:(^), :arb_pow))
    @eval begin
        function ($A)(x::ArbReal{P}, y::ArbReal{P}) where {P}
            z = ArbReal{P}()
            ccall(@libarb($F), Cvoid, (Ref{ArbReal}, Ref{ArbReal}, Ref{ArbReal}, Clong), z, x, y, P)
            return z
         end
    end
end

function (^)(x::ArbFloat{P}, y::ArbFloat{P}) where {P}
    rx = ArbReal{P}(x)
    ry = ArbReal{P}(y)
    r = rx^ry
    return ArbFloat{P}(r)
end

for (A,F) in ((:(+), :acb_add), (:(-), :acb_sub), (:(*), :acb_mul), (:(/), :acb_div), (:(^), :acb_pow))
    @eval begin
        function ($A)(x::ArbComplex{P}, y::ArbComplex{P}) where {P}
            z = ArbComplex{P}()
            ccall(@libarb($F), Cvoid, (Ref{ArbComplex}, Ref{ArbComplex}, Ref{ArbComplex}, Clong), z, x, y, P)
            return z
         end
    end
end

(+)(x::ArbComplex{P}, y::ArbReal{P}) where {P} = (+)(promote(x, y)...,)
(-)(x::ArbComplex{P}, y::ArbReal{P}) where {P} = (-)(promote(x, y)...,)
(*)(x::ArbComplex{P}, y::ArbReal{P}) where {P} = (*)(promote(x, y)...,)
(/)(x::ArbComplex{P}, y::ArbReal{P}) where {P} = (/)(promote(x, y)...,)
(^)(x::ArbComplex{P}, y::ArbReal{P}) where {P} = (^)(promote(x, y)...,)
(+)(x::ArbReal{P}, y::ArbComplex{P}) where {P} = (+)(promote(x, y)...,)
(-)(x::ArbReal{P}, y::ArbComplex{P}) where {P} = (-)(promote(x, y)...,)
(*)(x::ArbReal{P}, y::ArbComplex{P}) where {P} = (*)(promote(x, y)...,)
(/)(x::ArbReal{P}, y::ArbComplex{P}) where {P} = (/)(promote(x, y)...,)
(^)(x::ArbReal{P}, y::ArbComplex{P}) where {P} = (^)(promote(x, y)...,)

(+)(x::ArbComplex{P}, y::ArbFloat{P}) where {P} = (+)(promote(x, y)...,)
(-)(x::ArbComplex{P}, y::ArbFloat{P}) where {P} = (-)(promote(x, y)...,)
(*)(x::ArbComplex{P}, y::ArbFloat{P}) where {P} = (*)(promote(x, y)...,)
(/)(x::ArbComplex{P}, y::ArbFloat{P}) where {P} = (/)(promote(x, y)...,)
(^)(x::ArbComplex{P}, y::ArbFloat{P}) where {P} = (^)(promote(x, y)...,)
(+)(x::ArbFloat{P}, y::ArbComplex{P}) where {P} = (+)(promote(x, y)...,)
(-)(x::ArbFloat{P}, y::ArbComplex{P}) where {P} = (-)(promote(x, y)...,)
(*)(x::ArbFloat{P}, y::ArbComplex{P}) where {P} = (*)(promote(x, y)...,)
(/)(x::ArbFloat{P}, y::ArbComplex{P}) where {P} = (/)(promote(x, y)...,)
(^)(x::ArbFloat{P}, y::ArbComplex{P}) where {P} = (^)(promote(x, y)...,)

(+)(x::ArbReal{P}, y::ArbFloat{P}) where {P} = (+)(promote(x, y)...,)
(-)(x::ArbReal{P}, y::ArbFloat{P}) where {P} = (-)(promote(x, y)...,)
(*)(x::ArbReal{P}, y::ArbFloat{P}) where {P} = (*)(promote(x, y)...,)
(/)(x::ArbReal{P}, y::ArbFloat{P}) where {P} = (/)(promote(x, y)...,)
(^)(x::ArbReal{P}, y::ArbFloat{P}) where {P} = (^)(promote(x, y)...,)
(+)(x::ArbFloat{P}, y::ArbReal{P}) where {P} = (+)(promote(x, y)...,)
(-)(x::ArbFloat{P}, y::ArbReal{P}) where {P} = (-)(promote(x, y)...,)
(*)(x::ArbFloat{P}, y::ArbReal{P}) where {P} = (*)(promote(x, y)...,)
(/)(x::ArbFloat{P}, y::ArbReal{P}) where {P} = (/)(promote(x, y)...,)
(^)(x::ArbFloat{P}, y::ArbReal{P}) where {P} = (^)(promote(x, y)...,)

(+)(x::Mag, y::ArbFloat{P}) where {P} = (+)(promote(x, y)...,)
(-)(x::Mag, y::ArbFloat{P}) where {P} = (-)(promote(x, y)...,)
(*)(x::Mag, y::ArbFloat{P}) where {P} = (*)(promote(x, y)...,)
(/)(x::Mag, y::ArbFloat{P}) where {P} = (/)(promote(x, y)...,)
(^)(x::Mag, y::ArbFloat{P}) where {P} = (^)(promote(x, y)...,)
(+)(x::ArbFloat{P}, y::Mag) where {P} = (+)(promote(x, y)...,)
(-)(x::ArbFloat{P}, y::Mag) where {P} = (-)(promote(x, y)...,)
(*)(x::ArbFloat{P}, y::Mag) where {P} = (*)(promote(x, y)...,)
(/)(x::ArbFloat{P}, y::Mag) where {P} = (/)(promote(x, y)...,)
(^)(x::ArbFloat{P}, y::Mag) where {P} = (^)(promote(x, y)...,)

(+)(x::Mag, y::ArbReal{P}) where {P} = (+)(promote(x, y)...,)
(-)(x::Mag, y::ArbReal{P}) where {P} = (-)(promote(x, y)...,)
(*)(x::Mag, y::ArbReal{P}) where {P} = (*)(promote(x, y)...,)
(/)(x::Mag, y::ArbReal{P}) where {P} = (/)(promote(x, y)...,)
(^)(x::Mag, y::ArbReal{P}) where {P} = (^)(promote(x, y)...,)
(+)(x::ArbReal{P}, y::Mag) where {P} = (+)(promote(x, y)...,)
(-)(x::ArbReal{P}, y::Mag) where {P} = (-)(promote(x, y)...,)
(*)(x::ArbReal{P}, y::Mag) where {P} = (*)(promote(x, y)...,)
(/)(x::ArbReal{P}, y::Mag) where {P} = (/)(promote(x, y)...,)
(^)(x::ArbReal{P}, y::Mag) where {P} = (^)(promote(x, y)...,)

for F in (:(+), :(-), :(*), :(/), :(^))
  @eval begin
    ($F)(x::ArbFloat{P}, y::Integer) where {P} = ($F)(promote(x, y)...,)
    ($F)(x::ArbReal{P}, y::Integer) where {P} = ($F)(promote(x, y)...,)
    ($F)(x::ArbComplex{P}, y::Integer) where {P} = ($F)(promote(x, y)...,)
    ($F)(x::ArbFloat{P}, y::AbstractFloat) where {P} = ($F)(promote(x, y)...,)
    ($F)(x::ArbReal{P}, y::AbstractFloat) where {P} = ($F)(promote(x, y)...,)
    ($F)(x::ArbComplex{P}, y::AbstractFloat) where {P} = ($F)(promote(x, y)...,)
    ($F)(x::ArbComplex{P}, y::Complex) where {P} = ($F)(promote(x, y)...,)
    ($F)(x::Integer, y::ArbFloat{P}) where {P} = ($F)(promote(x, y)...,)
    ($F)(x::Integer, y::ArbReal{P}) where {P} = ($F)(promote(x, y)...,)
    ($F)(x::Integer, y::ArbComplex{P}) where {P} = ($F)(promote(x, y)...,)
    ($F)(x::IEEEFloat, y::ArbFloat{P}) where {P} = ($F)(promote(x, y)...,)
    ($F)(x::IEEEFloat, y::ArbReal{P}) where {P} = ($F)(promote(x, y)...,)
    ($F)(x::IEEEFloat, y::ArbComplex{P}) where {P} = ($F)(promote(x, y)...,)
    ($F)(x::Complex, y::ArbComplex{P}) where {P} = ($F)(promote(x, y)...,)
  end
end

# interact with Bool

@inline ArbFloat{P}(b::Bool) where {P} = b ? one(ArbFloat{P}) : zero(ArbFloat{P})
@inline Base.Bool(x::ArbFloat{P}) where {P} = iszero(x) ? false : (isone(x) ? true : throw(InexactError("")))
@inline ArbReal{P}(b::Bool) where {P} = b ? one(ArbReal{P}) : zero(ArbReal{P})
@inline Base.Bool(x::ArbReal{P}) where {P} = iszero(x) ? false : (isone(x) ? true : throw(InexactError("")))
@inline ArbComplex{P}(b::Bool) where {P} = b ? one(ArbComplex{P}) : zero(ArbComplex{P})
@inline Base.Bool(x::ArbComplex{P}) where {P} = iszero(x) ? false : (isone(x) ? true : throw(InexactError("")))

(+)(x::ArbFloat{P}, b::Bool) where {P} = b ? one(ArbFloat{P}) + x : zero(ArbFloat{P})
(+)(b::Bool, x::ArbFloat{P}) where {P} = b ? x + one(ArbFloat{P}) : zero(ArbFloat{P})

(+)(x::ArbReal{P}, b::Bool) where {P} = b ? one(ArbReal{P}) + x : zero(ArbReal{P})
(+)(b::Bool, x::ArbReal{P}) where {P} = b ? x + one(ArbReal){P} : zero(ArbReal{P})

(+)(x::ArbComplex{P}, b::Bool) where {P} = b ? one(ArbComplex{P}) + x : zero(ArbComplex{P})
(+)(b::Bool, x::ArbComplex{P}) where {P} = b ? x + one(ArbComplex){P} : zero(ArbComplex{P})

(-)(x::ArbFloat{P}, b::Bool) where {P} = b ? one(ArbFloat{P}) - x : zero(ArbFloat{P})
(-)(b::Bool, x::ArbFloat{P}) where {P} = b ? x - one(ArbFloat{P}) : zero(ArbFloat{P})

(-)(x::ArbReal{P}, b::Bool) where {P} = b ? one(ArbReal{P}) - x : zero(ArbReal{P})
(-)(b::Bool, x::ArbReal{P}) where {P} = b ? x - one(ArbReal){P} : zero(ArbReal{P})

(-)(x::ArbComplex{P}, b::Bool) where {P} = b ? one(ArbComplex{P}) - x : zero(ArbComplex{P})
(-)(b::Bool, x::ArbComplex{P}) where {P} = b ? x - one(ArbComplex){P} : zero(ArbComplex{P})

(*)(x::ArbFloat{P}, b::Bool) where {P} = b ? x : zero(ArbFloat{P})
(*)(b::Bool, x::ArbFloat{P}) where {P} = b ? x : zero(ArbFloat{P})

(*)(x::ArbReal{P}, b::Bool) where {P} = b ? x : zero(ArbReal{P})
(*)(b::Bool, x::ArbReal{P}) where {P} = b ? x : zero(ArbReal{P})

(*)(x::ArbComplex{P}, b::Bool) where {P} = b ? x : zero(ArbComplex{P})
(*)(b::Bool, x::ArbComplex{P}) where {P} = b ? x : zero(ArbComplex{P})

(\)(x::ArbFloat{P}, b::Bool) where {P} = b ? x : inf(ArbFloat{P})
(\)(b::Bool, x::ArbFloat{P}) where {P} = b ? inv(x) : zero(ArbFloat{P})

(\)(x::ArbReal{P}, b::Bool) where {P} = b ? x : inf(ArbReal{P})
(\)(b::Bool, x::ArbReal{P}) where {P} = b ? inv(x) : zero(ArbReal{P})

(\)(x::ArbComplex{P}, b::Bool) where {P} = b ? x : inf(ArbComplex{P})
(\)(b::Bool, x::ArbComplex{P}) where {P} = b ? inv(x) : zero(ArbComplex{P})


text:
using DiscO
using Random
data = db()

meta(data, "name", "(1+1) EA")
meta(data, "description", "This is for test.")

# Size of the bitstring.
n = 100

# Instance contains data relating to the problem, e.g. Graph
# Here, we don't have any specific information regarding the problem just the size
instance = Dict(["size"=> n])

# We want to optimize OneMax funciton, OneMax := sum(candidate)
fitness = onemax

# Using bitrand we create a string with n bits chosen
# unifromly random from 0 and 1
start = bitrand(n)
candidate = (point=start, fitness=fitness(instance, start))

# The algorithm must be stop when the fitness becomes n
stop(instance, best, iter) = (best.fitness == instance["size"] || iter == 100000)

# We use static bit mutation with probability 1/n.
mutation = staticBitMutation(1/n)

# Run algorithm (1+1) EA with standard mutation probability
optimum = ea1p1(instance, candidate, mutation, fitness, stop, data)

max_stars_repo_name: ahrajabi/DiscO

max_issues_repo_name: ahrajabi/DiscO

max_forks_repo_name: ahrajabi/DiscO


text:
# Copyright 2018 IBM Corp.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

module TestDoctrine
using Test

using Catlab, Catlab.WiringDiagrams, Catlab.Graphics
import Catlab.Graphics: Graphviz
using SemanticFlowGraphs.Doctrine

# Monocl category
#################

A, B, C = Ob(Monocl, :A, :B, :C)
A0, B0, C0, A1, B1, C1 = Ob(Monocl, :A0, :B0, :C0, :A1, :B1, :C1)
I = munit(Monocl.Ob)
f = Hom(:f, A, B)
f0, f1 = Hom(:f0, A0, B0), Hom(:f1, A1, B1)

# Subobjects
subA = SubOb(A0, A)
subB = SubOb(B0, B)
@test dom(subA) == A0
@test codom(subA) == A

sub = compose(SubOb(A0, A), SubOb(A, A1))
@test dom(sub) == A0
@test codom(sub) == A1
@test dom(subob_id(A)) == A
@test codom(subob_id(A)) == A
@test_throws SyntaxDomainError compose(subA, subB)

sub = otimes(subA, subB)
@test dom(sub) == otimes(A0,B0)
@test codom(sub) == otimes(A,B)

# Submorphisms
subf = SubHom(f0, f, subA, subB)
@test dom(subf) == f0
@test codom(subf) == f
@test subob_dom(subf) == subA
@test subob_codom(subf) == subB

subf1 = SubHom(f, f1, SubOb(A, A1), SubOb(B, B1))
sub = compose(subf, subf1)
@test dom(sub) == f0
@test codom(sub) == f1
@test subob_dom(sub) == compose(SubOb(A0,A), SubOb(A,A1))
@test subob_codom(sub) == compose(SubOb(B0,B), SubOb(B,B1))
@test dom(subhom_id(f)) == f
@test codom(subhom_id(f)) == f
@test subob_dom(subhom_id(f)) == subob_id(A)
@test subob_codom(subhom_id(f)) == subob_id(B)

g, g0 = Hom(:g, B, C), Hom(:g0, B0, C0)
subg = SubHom(g0, g, SubOb(B0, B), SubOb(C0, C))
sub = compose2(subf, subg)
@test dom(sub) == compose(f0,g0)
@test codom(sub) == compose(f,g)
@test subob_dom(sub) == SubOb(A0, A)
@test subob_codom(sub) == SubOb(C0, C)

# Explicit coercions
@test dom(coerce(subA)) == A0
@test codom(coerce(subA)) == A
@test coerce(subob_id(A)) == id(A)
@test compose(coerce(subob_id(A)), f) == f
@test compose(f, coerce(subob_id(B))) == f

# Constructors
@test dom(construct(A)) == I
@test codom(construct(A)) == A
@test dom(construct(f)) == B
@test codom(construct(f)) == A

# Pairs
f, g = Hom(:f, A, B), Hom(:g, A, C)
@test dom(pair(f,g)) == A
@test codom(pair(f,g)) == otimes(B,C)
@test pair(f,g) == compose(mcopy(A), otimes(f,g))
@test dom(pair(A0,f,g)) == A0
@test codom(pair(A0,f,g)) == otimes(B,C)

# Monocl wiring diagram
#######################

A0, A, B0, B = Ob(Monocl, :A0, :A, :B0, :B)
f = Hom(:f, A, B)
g = Hom(:g, B, A)

diagram = to_wiring_diagram(f)
@test boxes(diagram) == [ Box(f) ]
@test input_ports(diagram) == [ A ]
@test output_ports(diagram) == [ B ]

# Coercion
sub = SubOb(A0, A)
diagram = to_wiring_diagram(compose(coerce(sub), f))
@test boxes(diagram) == [ Box(f) ]
@test input_ports(diagram) == [ A0 ]
@test output_ports(diagram) == [ B ]

# Graphviz support.
diagram = to_wiring_diagram(compose(coerce(SubOb(A0,A)), construct(g)))
@test isa(to_graphviz(diagram), Graphviz.Graph)

end

max_stars_repo_name: cclauss/semanticflowgraph

max_issues_repo_name: cclauss/semanticflowgraph

max_forks_repo_name: cclauss/semanticflowgraph


text:
# ## Loading modules
# using Distributed
# using DataFrames
# using CSV
#
# addprocs(6)
#
# ## Load modules in parallel
# @everywhere include("RP_Cournot/modules/MyModules.jl")
#
# @everywhere begin
#     using .MyModules
#     using Random
#     using SharedArrays
#     ## Setting seed
#     Random.seed!(66636)
# end
#
# ## Settin up directories
#
# dir = @__DIR__ ;
# datadir=dir*"/data/oildata.csv" ;
# resultsdir=dir*"/results/"
#
# ## Loading data
#
# pt, qt = load_data_select(datadir)

## Initializing array for results

mchain=SharedArray{Float64,3}(zeros(Integer(size(pt,1)/6),size(qt,2),Param.nsim))
decade=SharedArray{Float64,2}(zeros(4,2))
acc_r_t=SharedArray{Float64,1}(zeros(5))
TS, pvalue = nothing, nothing
## Testing

# Markov Chain
@time MyModules.shared_fill_chain!(mchain, pt, qt, fill_chain_sem_algo!)
# Joint test
@time TS, pvalue, acc_r_t[1] = MyModules.RP_Cournot_test_par(mchain)
# Testing by decade
@time begin
    dec = [0, 16, 36, 56, 72]
    for i=1:4
        m = mchain[dec[i]+1:dec[i+1],:,:]
        decade[i,1], decade[i,2], acc_r_t[i+1] = MyModules.RP_Cournot_test_par(m)
    end
end

## Saving results to csv

joint_sem_algo=DataFrame(Time="1980-2008",
                 TS=TS,
                 pvalue=pvalue,
                 stars=mystars.(pvalue))
CSV.write(resultsdir*"joint_ak_$set.csv", joint_sem_algo)

decade_sem_algo=DataFrame(Time=["1973-1980", "1981-1990","1991-2000","2001-2008"],
                 TS=decade[:,1],
                 pvalue=decade[:,2],
                 stars=mystars.(decade[:,2]))
CSV.write(resultsdir*"decade_ak_$set.csv", decade_sem_algo)

acc_rates_algo=DataFrame(Ar=acc_r_t)
CSV.write(resultsdir*"accp_r_ak_$set.csv",acc_rates_algo)

max_stars_repo_name: hans-mtz/SCR

max_issues_repo_name: hans-mtz/SCR

max_forks_repo_name: hans-mtz/SCR


text:
using AlignedSpans
using Documenter

makedocs(; modules=[AlignedSpans],
         sitename="AlignedSpans",
         authors="Beacon Biosignals, Inc.",
         pages=["Introduction" => "index.md",
                "API Documentation" => "API.md"],
         strict=true)

deploydocs(; repo="github.com/beacon-biosignals/AlignedSpans.jl.git",
           devbranch="main",
           push_preview=true)


text:
struct BoxList{P<:AbstractBoxPartition,L<:AbstractVector}
    partition::P
    keylist::L
end

Base.length(list::BoxList) = length(list.keylist)
Base.getindex(list::BoxList, x::AbstractVector) = BoxList(list.partition, list.keylist[x])
Base.getindex(list::BoxList, i::Int64) = list.keylist[i]

BoxList(set::BoxSet) = BoxList(set.partition, collect(set.set))

BoxSet(list::BoxList) = BoxSet(list.partition, Set(list.keylist))

struct TransferOperator{L<:BoxList,W}
    vertices::L
    edges::Dict{Tuple{Int,Int},W}
end

function Base.show(io::IO, T::TransferOperator)
    n = length(T.vertices)
    m = length(T.edges)
    print(io, "TransferOperator on $(n) boxes with $(m) edges")
end

function invert_vector(x::AbstractVector{T}) where T
    dict = Dict{T,Int}()
    sizehint!(dict, length(x))

    for i in 1:length(x)
        dict[x[i]] = i
    end

    return dict
end

# TODO: this code is generally incorrect. only valid for BoxPartition and special choices of points
function TransferOperator(g::SampledBoxMap, boxset::BoxSet)
    Q = g.domain
    n = length(g.domain_points(Q.center, Q.radius))
    P = boxset.partition
    edges = [ Dict{Tuple{Int64,Int64},Float64}() for k = 1:nthreads() ]
    boxlist = BoxList(boxset)
    key_to_index = invert_vector(boxlist.keylist)

    @threads for i = 1:length(boxlist)
        box = key_to_box(P, boxlist[i])
        c, r = box.center, box.radius
        points = g.domain_points(c, r)
        for p in points
            fp = g.map(@muladd p .* r .+ c)
            hit = point_to_key(P, fp)
            if !isnothing(hit)
                if hit in boxset.set
                    j = key_to_index[hit]
                    e = (i,j)
                    edges[threadid()][e] = get(edges[threadid()], e, 0) + 1.0/n
                end
            end
        end
    end
    edges = merge(edges...)

    return TransferOperator(boxlist, edges)
end

function strongly_connected_components(gstar::TransferOperator)
    graph = SimpleDiGraph(
        [Edge(edge[1], edge[2]) for edge in keys(gstar.edges)]
    )

    sccs = LightGraphs.strongly_connected_components(graph)

    connected_vertices = Int[]

    for scc in sccs
        if length(scc) > 1 || has_edge(graph, scc[1], scc[1])
            append!(connected_vertices, scc)
        end
    end

    return BoxSet(gstar.vertices[connected_vertices])
end

function matrix(gstar::TransferOperator{L,W}) where {L,W}
    num_edges = length(gstar.edges)

    I = Vector{Int}()
    sizehint!(I, num_edges)

    J = Vector{Int}()
    sizehint!(J, num_edges)

    V = Vector{W}()
    sizehint!(V, num_edges)

    for (edge, weight) in gstar.edges
        push!(I, edge[1])
        push!(J, edge[2])
        push!(V, weight)
    end

    n = length(gstar.vertices)
    return sparse(I, J, V, n, n)
end

function Arpack.eigs(gstar::TransferOperator{BoxList{P,L}}; nev::Int=1) where {P,L}
    G = matrix(gstar)

    λ, ϕ, nconv = eigs(G'; nev=nev)

    ϕ_funs = BoxFun{P,keytype(P),ComplexF64}[]

    for i in 1:nev
        ϕi = ϕ[:,i]
        dict = Dict{keytype(P),ComplexF64}()
        sizehint!(dict, length(ϕi))

        for j in 1:length(ϕi)
            dict[gstar.vertices.keylist[j]] = ϕi[j]
        end

        ϕ_fun = BoxFun{P,keytype(P),ComplexF64}(gstar.vertices.partition, dict)
        normalize!(ϕ_fun)
        push!(ϕ_funs, ϕ_fun)
    end

    return λ, ϕ_funs, nconv
end


text:
### A Pluto.jl notebook ###
# v0.12.11

using Markdown
using InteractiveUtils

# ╔═╡ 9e6dfb4a-f7c4-11ea-0371-03ff747c9f61
using Pkg, DrWatson

# ╔═╡ 9e6e3ad6-f7c4-11ea-0980-81a920a19f14
begin
  @quickactivate "StatisticalRethinkingStan"
  using DynamicHMC
  using LogDensityProblems, TransformVariables
  using DynamicHMC
  using StatisticalRethinking
end

# ╔═╡ f385d68a-f7c3-11ea-3fc2-ad7453c86414
md"## Intro-dhmc-2d.jl"

# ╔═╡ 9e6ec7aa-f7c4-11ea-37b6-c30c9bc9b70b
md"### snippet 9.3"

# ╔═╡ 9e7cb554-f7c4-11ea-3a22-5d7dca89f1ce
md"##### Construct the logdensity problem."

# ╔═╡ 9e7d4954-f7c4-11ea-37ec-a53f2a326de6
Base.@kwdef mutable struct clip_9_3_1_model{
  TY <: AbstractVector, TX <: AbstractVector}
    "Observations."
    y::TY
    "Covariate"
    x::TX
end

# ╔═╡ 9e8a56c8-f7c4-11ea-1b9c-f1874c10c04b
md"##### Write a function to return properly dimensioned transformation."

# ╔═╡ 9e8af36e-f7c4-11ea-0bc2-c770176975d7
function make_transformation(model::clip_9_3_1_model)
  as((muy = asℝ, mux = asℝ))
end

# ╔═╡ 9e97b33e-f7c4-11ea-17f0-cf849e4a1164
md"##### Instantiate the model with data and inits."

# ╔═╡ 9e9ff8c8-f7c4-11ea-1e71-c7f189d3f5ed
Random.seed!(1234591);

# ╔═╡ 9eaa03c4-f7c4-11ea-2b34-89e2f92810cc
begin
	N = 100
	x = rand(Normal(0, 1), N)
	y = rand(Normal(0, 1), N)
	model1 = clip_9_3_1_model(;y=y, x=x)
end

# ╔═╡ 9eac4c1c-f7c4-11ea-1fd4-4573a58933a7
θ = (muy = 0.0, mux=0.0)

# ╔═╡ 9ebc38fa-f7c4-11ea-2456-933dc993b6ab
md"##### Make the type callable with the parameters *as a single argument*."

# ╔═╡ 9ebe4348-f7c4-11ea-1ac0-8dd4036acfef
function (model:: clip_9_3_1_model)(θ)
    @unpack y, x, = model    # extract the data
    @unpack muy, mux = θ     # works on the named tuple too
    ll = 0.0
    ll += loglikelihood(Normal(mux, 1), x)
    ll += loglikelihood(Normal(muy, 1), y)
    ll += logpdf(Normal(0, 1), mux)
    ll += logpdf(Normal(0, 1), muy)
    ll
end

# ╔═╡ 9ec7d1e4-f7c4-11ea-0c8a-e36dbaa21dd7
model1(θ)

# ╔═╡ 9ed0a3c6-f7c4-11ea-1857-1bfa81bb165b
md"##### Wrap the problem with a transformation, then use Flux for the gradient."

# ╔═╡ 9edac1f8-f7c4-11ea-2964-b56a3c28fb3c
begin
	P = TransformedLogDensity(make_transformation(model1), model1)
	∇P = ADgradient(:ForwardDiff, P);
end

# ╔═╡ 9ee47d1a-f7c4-11ea-0353-157e1b36e02e
md"##### Tune and sample."

# ╔═╡ 9eecc290-f7c4-11ea-101e-0d4097815b8f
results = mcmc_with_warmup(Random.GLOBAL_RNG, ∇P, 1000);

# ╔═╡ 9ef6da32-f7c4-11ea-3f73-a5617190ffd2
md"##### We use the transformation to obtain the posterior from the chain."

# ╔═╡ 9effe6ae-f7c4-11ea-19cd-6d79b43b4278
posterior = P.transformation.(results.chain)

# ╔═╡ e3ae2806-2b90-11eb-070f-39ddd6208ccd
typeof(posterior)

# ╔═╡ 9f087fbc-f7c4-11ea-3264-9b235c4757c8
md"##### Extract the posterior means,"

# ╔═╡ 9f145896-f7c4-11ea-1104-b3a151572004
begin
	[mean(first, posterior), mean(last, posterior)]
	a3d = Array{Float64, 3}(undef, 1000, 2, 1)
	for j in 1:1
	  for i in 1:1000
		a3d[i, 1, j] = values(posterior[i].muy)
		a3d[i, 2, j] = values(posterior[i].mux)
	  end
	end
end

# ╔═╡ 9f1eba20-f7c4-11ea-278e-8b535e3f4fe5
begin
	pnames = ["muy", "mux"]
	sections = Dict(:parameters =>pnames,)
	chn = MCMCChains.Chains(a3d, pnames, sections, start=1)
	CHNS(chn)
end

# ╔═╡ 9f282556-f7c4-11ea-2b1d-6d19a869b9bf
# Draw 200 samples:

function draw_n_samples(model, grad;
  epsilon = 0.03, # Step size
  L = 11, # No of leapfrog steps
  n_samples = 1000, # No of samples
  q = [-0.1, 0.2]) # Initial position

  samples = zeros(n_samples, 2)
  for i in 1:n_samples
    q, ptraj, qtraj, accept,
      dH = StatisticalRethinking.HMC(model, grad, 0.03, 11, q)
    samples[i, :] = q
  end

  samples
end

# ╔═╡ 9f33c5ca-f7c4-11ea-394c-59b681e11b40
begin
	samples = draw_n_samples(model1, ∇P; n_samples=200);
	mean(samples, dims=1)
end

# ╔═╡ 9f49119e-f7c4-11ea-1d55-dfcf0d3f164a
md"## End of intro-dhmc-2d.jl"

# ╔═╡ Cell order:
# ╟─f385d68a-f7c3-11ea-3fc2-ad7453c86414
# ╠═9e6dfb4a-f7c4-11ea-0371-03ff747c9f61
# ╠═9e6e3ad6-f7c4-11ea-0980-81a920a19f14
# ╟─9e6ec7aa-f7c4-11ea-37b6-c30c9bc9b70b
# ╟─9e7cb554-f7c4-11ea-3a22-5d7dca89f1ce
# ╠═9e7d4954-f7c4-11ea-37ec-a53f2a326de6
# ╟─9e8a56c8-f7c4-11ea-1b9c-f1874c10c04b
# ╠═9e8af36e-f7c4-11ea-0bc2-c770176975d7
# ╟─9e97b33e-f7c4-11ea-17f0-cf849e4a1164
# ╠═9e9ff8c8-f7c4-11ea-1e71-c7f189d3f5ed
# ╠═9eaa03c4-f7c4-11ea-2b34-89e2f92810cc
# ╠═9eac4c1c-f7c4-11ea-1fd4-4573a58933a7
# ╠═9ebc38fa-f7c4-11ea-2456-933dc993b6ab
# ╠═9ebe4348-f7c4-11ea-1ac0-8dd4036acfef
# ╠═9ec7d1e4-f7c4-11ea-0c8a-e36dbaa21dd7
# ╟─9ed0a3c6-f7c4-11ea-1857-1bfa81bb165b
# ╠═9edac1f8-f7c4-11ea-2964-b56a3c28fb3c
# ╟─9ee47d1a-f7c4-11ea-0353-157e1b36e02e
# ╠═9eecc290-f7c4-11ea-101e-0d4097815b8f
# ╠═9ef6da32-f7c4-11ea-3f73-a5617190ffd2
# ╠═9effe6ae-f7c4-11ea-19cd-6d79b43b4278
# ╠═e3ae2806-2b90-11eb-070f-39ddd6208ccd
# ╠═9f087fbc-f7c4-11ea-3264-9b235c4757c8
# ╠═9f145896-f7c4-11ea-1104-b3a151572004
# ╠═9f1eba20-f7c4-11ea-278e-8b535e3f4fe5
# ╠═9f282556-f7c4-11ea-2b1d-6d19a869b9bf
# ╠═9f33c5ca-f7c4-11ea-394c-59b681e11b40
# ╟─9f49119e-f7c4-11ea-1d55-dfcf0d3f164a


text:
export macau_blocked

function macau_blocked(data::RelationData;
              num_latent::Int = 10,
              burnin          = 500,
              psamples        = 200,
              verbose::Bool   = true)
end


text:
module CausalityTools
    using Reexport
    import DynamicalSystemsBase: trajectory, DiscreteDynamicalSystem, ContinuousDynamicalSystem
    import DelayEmbeddings: Dataset
    export trajectory, DiscreteDynamicalSystem, ContinuousDynamicalSystem, Dataset

    import TransferEntropy
    import TransferEntropy: transferentropy, mutualinfo, conditional_mutualinfo, Hilbert
    export transferentropy, mutualinfo, conditional_mutualinfo, Hilbert
    @reexport using Entropies
    @reexport using TransferEntropy
    @reexport using TimeseriesSurrogates

    include("JointDistanceDistribution/JointDistanceDistribution.jl")
    include("CrossMappings/CrossMappings.jl")
    include("SMeasure/smeasure.jl")
    include("PredictiveAsymmetry/PredictiveAsymmetry.jl")
    include("example_systems/ExampleSystems.jl")

    using Requires
    function __init__()
        #@require UncertainData="dcd9ba68-c27b-5cea-ae21-829cd07325bf" begin
        #   include("integrations/uncertaindata.jl")
        #end

        #@require Simplices="d5428e67-3037-59ba-9ab1-57a04f0a3b6a" begin
        # #   import PerronFrobenius: SimplexExact, SimplexPoint
        #    export SimplexExact, SimplexPoint
        #end
    end
end


text:
### A Pluto.jl notebook ###
# v0.14.0

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 0db6ee04-81b7-11eb-330c-11b578b72c90
html"""

Section 1.9

 Taking Advantage of Structure

Lecture Video

"""

# ╔═╡ ca1a1072-81b6-11eb-1fee-e7df687cc314
PlutoUI.TableOfContents(aside = true)

# ╔═╡ b0ba5b8c-f5d1-11ea-1304-3f0e47f935fe
md"# Examples of structure"

# ╔═╡ 864e1180-f693-11ea-080e-a7d5aabc9ca5
begin
	import Pkg
    Pkg.activate(mktempdir())
    Pkg.add([
		Pkg.PackageSpec(name="ImageIO", version="0.5"),
		Pkg.PackageSpec(name="ImageShow", version="0.2"),
		Pkg.PackageSpec(name="FileIO", version="1.6"),
		Pkg.PackageSpec(name="PNGFiles", version="0.3.6"),
		Pkg.PackageSpec(name="Colors", version="0.12"),
		Pkg.PackageSpec(name="ColorVectorSpace", version="0.8"),
        Pkg.PackageSpec(name="ColorSchemes", version="3"),

        Pkg.PackageSpec(name="PlutoUI", version="0.7"),
    ])
	using Colors, ColorVectorSpace, ImageShow, FileIO
	using ImageShow.ImageCore
	using ColorSchemes

	using InteractiveUtils, PlutoUI
	using LinearAlgebra, SparseArrays, Statistics
end

# ╔═╡ 261c4df2-f5d2-11ea-2c72-7d4b09c46098
md"""
# One-hot vectors
This example comes from machine learning.
"""

# ╔═╡ 3cada3a0-81cc-11eb-04c8-bde26d36a84e
md"""
A **one-hot** vector has a single "hot" element, i.e. a single 1 in a sea of zeros. For example:
"""

# ╔═╡ fe2028ba-f6dc-11ea-0228-938a81a91ace
myonehatvector = [0, 1, 0, 0, 0, 0]

# ╔═╡ 8d2c6910-f5d4-11ea-1928-1baf09815687
md"""How much "information" (numbers) do you need to represent a one-hot vector? Is it $n$ or is it two?
"""

# ╔═╡ 54649792-81cc-11eb-1038-9161a4037acf
md"""
(There are also "1-cold" vectors:)
"""

# ╔═╡ 0a902426-f6dd-11ea-0ae4-fb0c47863fe7
# also one "cold"
1 .- myonehatvector

# ╔═╡ 4794e860-81b7-11eb-2c91-8561c20f308a
md"""
## Julia: `structs` (creating a new type in Julia)
"""

# ╔═╡ 67827da8-81cc-11eb-300e-278104d2d958
md"""
We can create our own new types in Julia. Let's create a new type to represent one-hot vectors. It will be a subtype of `AbstractVector`, meaning that it *behaves like* a vector.
"""

# ╔═╡ 4624cd26-f5d3-11ea-1cf8-7d6555eae5fa
struct OneHot <: AbstractVector{Int}
	n::Int
	k::Int
end

# ╔═╡ 9bdabef8-81cc-11eb-14a1-67a9a7d968c0
md"""
We need to specify how long the vector is:
"""

# ╔═╡ 397ac764-f5fe-11ea-20cc-8d7cab19d410
Base.size(x::OneHot) = (x.n, )

# ╔═╡ a22dcd2c-81cc-11eb-1252-13ace134192d
md"""
and how to extract the $i$th component:
"""

# ╔═╡ 82c7046c-f5d3-11ea-04e2-ef7c0f4db5da
Base.getindex(x::OneHot, i::Int) = Int(x.k == i)

# ╔═╡ b024c318-81cc-11eb-018c-e1f7830ff51b
md"""
(Note that `x.k == i` returns a Boolean value, `true` or `false`, which we are then converting to an `Int`.)
"""

# ╔═╡ 93bfe3ac-f756-11ea-20fb-8f7d586b42f3
myonehotvector = OneHot(6,2)

# ╔═╡ 175039aa-f758-11ea-251a-5db57d7c4b32
myonehotvector[3]

# ╔═╡ c2a4b0a2-81cc-11eb-37a7-db601a6ddfdf
myonehotvector[2]

# ╔═╡ c5ed7d3e-81cc-11eb-3386-15b72db8155d
md"""
This behaves as if it were the original vector, but we are storing only 2 integers.
This is an example of taking advantage of structure.
"""

# ╔═╡ e2e354a8-81b7-11eb-311a-35151063c2a7
md"""
## Julia: dump and Dump
"""

# ╔═╡ dc5a96ba-81cc-11eb-3189-25920df48afa
md"""
`dump` shows the internal data stored inside a given object:
"""

# ╔═╡ af0d3c22-f756-11ea-37d6-11b630d2314a
with_terminal() do
	dump(myonehotvector)
end

# ╔═╡ 06e3a842-26b8-4417-9cf5-8a083ccdb264
md"""
Since `dump` writes to a terminal, you can't use it directly inside Pluto, which is why we wrap it inside a `with_terminal` block. You can also use the function `Dump` from `PlutoUI`:
"""

# ╔═╡ 91172a3e-81b7-11eb-0953-9f5e0207f863
Dump(myonehotvector)

# ╔═╡ 4bbf3f58-f788-11ea-0d24-6b0fb070829e
myonehotvector

# ╔═╡ fe70d104-81b7-11eb-14d0-eb5237d8ea6c
md"""
### Visualizing a one-hot vector
"""

# ╔═╡ ef8f44b2-f5fc-11ea-1e4d-bd873cd39d6c
md"""
n=$(@bind nn Slider(1:20, show_value=true))
"""

# ╔═╡ fd9211c0-f5fc-11ea-1745-7f2dae88af9e
md"""
k=$(@bind kk Slider(1:nn, default=1, show_value=true))
"""

# ╔═╡ f1154df8-f693-11ea-3b16-f32835fcc470
x = OneHot(nn, kk)

# ╔═╡ 982590d4-f5ff-11ea-3802-73292c75ad6c
show_image(x)

# ╔═╡ 81c35324-f5d4-11ea-2338-9f982d38732c
md"# Diagonal matrices"

# ╔═╡ 2cfda0dc-f5d5-11ea-16c4-b5a33b90e37f
md"Another example is diagonal matrices. Here's how you might see them in high school:"

# ╔═╡ 150432d4-f5d5-11ea-32b2-19a2a91d9637
denseD = [5 0 0
	    0 6 0
	    0 0 -10]

# ╔═╡ 44215aa4-f695-11ea-260e-b564c6fbcd4a
md"Julia has a better way of representing them:"

# ╔═╡ 21328d1c-f5d5-11ea-288e-4171ad35326d
D = Diagonal(denseD)

# ╔═╡ 75761cc0-81cd-11eb-1186-7d47debd68ca
md"""
It even displays nicely, with dots instead of zeros.

We can also create a diagonal matrix from the values on the diagonal:
"""

# ╔═╡ 6bd8a886-f758-11ea-2587-870a3fa9d710
Diagonal([5, 6, -10])

# ╔═╡ 4c533ac6-f695-11ea-3724-b955eaaeee49
md"How much information is stored for each representation? We can use Julia's `dump` function to find out:"

# ╔═╡ 466901ea-f5d5-11ea-1db5-abf82c96eabf
Dump(denseD)

# ╔═╡ b38c4aae-f5d5-11ea-39b6-7b0c7d529019
Dump(D)

# ╔═╡ 93e04ed8-81cd-11eb-214a-a761ef8c406f
md"""
We see that `Diagonal` stores only the diagonal entries, not the zeros!
"""

# ╔═╡ e90c55fc-f5d5-11ea-10f1-470ff772985d
md"""We should always look for *structure* where it exists!"""

# ╔═╡ 19775c3c-f5d6-11ea-15c2-89618e654a1e
md"# Sparse matrices"

# ╔═╡ 653792a8-f695-11ea-1ae0-43761c502583
md"A *sparse* matrix is a matrix that has many zeros, and is hence worth storing in a *sparse* representation:"

# ╔═╡ 79c94d2a-f75a-11ea-031d-09d70d229e15
denseM = [0 0 9; 0 0 0;12 0 4]

# ╔═╡ 10bc5d50-81b9-11eb-2ac7-354a6c6c826b
md"""
The above displays a sparse matrix in so-called `(i, j, value)` form.  We could
store sparse matrices in this way:
"""

# ╔═╡ 77d6a952-81ba-11eb-24e3-cb6510a59455
M = sparse(denseM)

# ╔═╡ 1f3ba55a-81b9-11eb-001f-593b9d8639ca
md"""
Although it looks like it's stored like this, in fact the actual storage format is different. In the Julia `SparseArrays.jl` package, the storage format
is **compressed sparse column** (CSC) format,
which is generally considered favorable for arithmetic, matrix-vector
	products and column slicing.  Of course, for specific matrices, other
		formats might be better.

* `nzval` contains the nonzero matrix entries
* `rowval` is the "i" or row entry for the corresponding value in nzval
  * length(rowval) == length(nzval)
* `colptr[j]` points into nzval and tells you the first nonzero in or after column j
* The last entry of colptr points beyond the end of nzval to indicate no more columns.
  * length(colptr) == number of columns + 1
"""

# ╔═╡ 3d4a702e-f75a-11ea-031c-333d591fc442
Dump(sparse(M))

# ╔═╡ 80ff4010-81bb-11eb-374e-215a57defb0b
md"""
 An example where CSC may not be a great choice is the following. The reason is that  `colptr` must have an entry in each column:
"""

# ╔═╡ 5de72b7c-f5d6-11ea-1b6f-35b830b5fb34
M2 = sparse([1, 2, 10^6], [4, 9, 10^6], [7, 8, 9])

# ╔═╡ 8b60629e-f5d6-11ea-27c8-d934460d3a57
with_terminal() do
	dump(M2)
end

# ╔═╡ 2fd7e52e-f5d7-11ea-3b5a-1f338e2451e0
M3 = [1 0 2 0 10; 0 3 4 0 9; 0 0 0 5 8; 0 0 0 0 7]

# ╔═╡ 2e87d4fe-81bc-11eb-0d16-b988bcedcc73
M4 = M3 .* 0

# ╔═╡ cde79f38-f5d6-11ea-3297-0b5b240f7b9e
Dump(sparse(M4))

# ╔═╡ aa09c008-f5d8-11ea-1bdc-b51ee6eb2478
sparse(M4)

# ╔═╡ 62a6ec62-f5d9-11ea-071e-ed33c5dea0cd
md"""# Random vectors
"""

# ╔═╡ 7f63daf6-f695-11ea-0b80-8702a83103a4
md"How much structure is there in a *random* vector?"

# ╔═╡ 67274c3c-f5d9-11ea-3475-c9d228e3bd5a
v = rand(1:9, 1_000_000)

# ╔═╡ 765c6552-f5d9-11ea-29d3-bfe7b4b04612
md"""You might guess that there is "no structure". But you can actually think of randomness itself as a structure.

For example, take the mean and standard deviation -- some would say that's the structure.

"""

# ╔═╡ 126fb3ea-f5da-11ea-2f7d-0b3259a296ce
mean(v), std(v), 5, sqrt(10 * 2/3)

# ╔═╡ ed0b2358-81ce-11eb-3339-93abcc06fd91
md"""
If you repeat the calculation, to 3 or 4 digits the mean and standard deviation don't change, and are approximately equal to the theoretical values on the right.
"""

# ╔═╡ 24ce92fa-81cf-11eb-30f0-b1e357d79d53
md"""
We can also count how many times each digit occurs in the data set:
"""

# ╔═╡ 2d4500e0-81cf-11eb-1699-d310074fddf5
[sum(v .== i) for i in 1:9]

# ╔═╡ 3546ff30-81cf-11eb-3afc-05c5db61366f
md"""
We see that each number occurs roughly the same number of times.
"""

# ╔═╡ 9b9e2c2a-f5da-11ea-369b-b513b196515b
md"Statisticians (and professors who've just graded exams) might say that under certain circumstances the mean and the variance give you the necessary structure, and the rest can be thrown away."

# ╔═╡ e68b98ea-f5da-11ea-1a9d-db45e4f80241
m = sum(v) / length(v)  # mean

# ╔═╡ f20ccac4-f5da-11ea-0e69-413b5e49f423
σ² = sum( (v .- m) .^ 2 ) / (length(v) - 1)

# ╔═╡ 12a2e96c-f5db-11ea-1c3e-494ae7446886
σ = sqrt(σ²)

# ╔═╡ 22487ce2-f5db-11ea-32e9-6f70ab2c0353
std(v)

# ╔═╡ 389ae62e-f5db-11ea-1557-c3adbbee0e5c
md"Sometimes the summary statistics are all you want. (But sometimes not.)"

# ╔═╡ 0c2b6408-f5d9-11ea-2b7f-7fece2eecc1f
md"# Multiplication tables"

# ╔═╡ 5d767290-f5dd-11ea-2189-81198fd216ce
outer(v, w) = [x * y for x ∈ v, y ∈ w]  # just a multiplication table

# ╔═╡ 587790ce-f6de-11ea-12d9-fde2a17ae314
outer(1:10, 1:10)

# ╔═╡ a39e8256-f6de-11ea-3170-c923b56609da
md"Did you memorize this in third grade?"

# ╔═╡ 8c84edd0-f6de-11ea-2180-61c6b81aac3b
@bind k Slider(1:14, show_value=true)

# ╔═╡ 22b73baa-f6df-11ea-197f-bbb4bd1a7ef5
outer(1:k, 1:k)

# ╔═╡ b2332814-f6e6-11ea-1c7d-556c7d4687f1
outer([2,4,6], [10,100,1000])

# ╔═╡ 9ab7a72e-81cf-11eb-2b78-073ff51cae58
md"""
A multiplication table is clearly a structure, but it's not sparse -- there are no zeros. Nonetheless you need much less information to reconstruct the matrix.
"""

# ╔═╡ fd8dd108-f6df-11ea-2f7c-3d99d054ac15
md"In the context of 1:k times 1:k, just one number k is needed."

# ╔═╡ 165788b2-f601-11ea-3e69-cdbbb6558e54
md"If you look at the following matrix? Does it have any structure? It's certainly more hidden."

# ╔═╡ 22941bb8-f601-11ea-1d6e-0d955297bc2e
outer( rand(3), rand(4) )  # but it's just a multiplication table

# ╔═╡ c33bf00e-81cf-11eb-1e1a-e5a879a45093
md"""
You might guess by visualizing the matrix that it is a multiplication table:
"""

# ╔═╡ 2f75df7e-f601-11ea-2fc2-aff4f335af33
show_image( outer( rand(10), rand(10) ))

# ╔═╡ 7ff664f0-f74b-11ea-0d2d-b53f19e4f4bf
md"We can factor out a multiplication table, if it's there:"

# ╔═╡ a0611eaa-81bc-11eb-1d23-c12ab14138b1
md"""
### Julia: Exceptions are thrown (generated) using `error`

An exception is anything that can interrupt a program, e.g. invalid input data.
"""

# ╔═╡ a4728944-f74b-11ea-03c3-9123908c1f8e
function factor( mult_table )
	v = mult_table[:, 1]
	w = mult_table[1, :]

	if v[1] ≠ 0
		w /= v[1]
	end

	# Good code has a check:
	if outer(v, w) ≈ mult_table
	   return v, w
	else
		error("Input is not a multiplication table")
	end
end

# ╔═╡ 05c35402-f752-11ea-2708-59cf5ef74fb4
factor( outer([1, 2, 3], [2, 2, 2] ) )

# ╔═╡ 8c11b19e-81bc-11eb-184b-bf6ffefe29de
md"""
A random 2x2 matrix is not a multiplication table. Most matrices are not given by multiplication tables.
"""

# ╔═╡ 8baaa994-f752-11ea-18d9-b3d0a6b9f7d9
factor( rand(2,2) )

# ╔═╡ d92bece4-f754-11ea-0242-99f198bb5b7b
md" Let's add two (or more) multiplication tables:"

# ╔═╡ e740999c-f754-11ea-2089-4b7a9aec6030
A = sum( outer(rand(3),rand(3)) for i=1:2 )

# ╔═╡ 0a79a7b4-f755-11ea-1b2d-21173567b257
md"Is it possible, given the matrix, to find the structure? E.g. to show that a matrix is a sum of outer products (multiplication table)."

# ╔═╡ 5adb98c2-f6e0-11ea-1fde-53b0fd6639c3
md"The answer is yes: The **Singular-Value Decomposition** (SVD) from algebra can find the structure!"

# ╔═╡ 487d6f9c-81d0-11eb-3bb0-336a4beb9b38
md"""
Let's take the SVD and calculate the sum of two outer products:
"""

# ╔═╡ 5a493052-f601-11ea-2f5f-f940412905f2
begin
	U, Σ, V = svd(A)

    outer( U[:, 1], V[:, 1] * Σ[1] ) + outer( U[:, 2], V[:, 2] * Σ[2] )
end

# ╔═╡ 55b76aee-81d0-11eb-0bcc-413f5bd14360
md"""
We see that we reconstruct the original matrix!"
"""

# ╔═╡ 709bf30a-f755-11ea-2e82-bd511e598c77
B = rand(3,3)

# ╔═╡ 782532b0-f755-11ea-1385-cd1a28c4b9d5
begin
	UU, ΣΣ, VV = svd( B )
    outer( UU[:,1], VV[:,1] * ΣΣ[1] ) + outer( UU[:,2], VV[:,2] * ΣΣ[2] )
end

# ╔═╡ 5bc4ab0a-f755-11ea-0fad-4987ad9fc02f
md"and it can approximate too!"

# ╔═╡ a5d637ea-f5de-11ea-3b70-877e876bc9c9
flag = outer([1,1,1,2,2,2,1,1,1], [1,1,1,1,1,1,1,1,1])

# ╔═╡ 21bbb60a-f5df-11ea-2c1b-dd716a657df8
cs = distinguishable_colors(100)

# ╔═╡ 2668e100-f5df-11ea-12b0-073a578a5edb
cs[flag]

# ╔═╡ e8d727f2-f5de-11ea-1456-f72602e81e0d
cs[flag + flag']

# ╔═╡ f5fcdeea-f75c-11ea-1fc3-731f0ef1ad14
outer([1,1,1,2,2,2,1,1,1], [1,1,1,1,1,1,1,1,1]) + outer([1,1,1,1,1,1,1,1,1], [1,1,1,2,2,2,1,1,1])

# ╔═╡ 0373fbf6-f75d-11ea-2a9e-cbb714d69cf4
cs[outer([1,1,1,2,2,2,1,1,1], [1,1,1,1,1,1,1,1,1]) + outer([1,1,1,1,1,1,1,1,1], [1,1,1,2,2,2,1,1,1])]

# ╔═╡ ebd72fb8-f5e0-11ea-0630-573337dff753
md"""
# Singular Value Decomposition (SVD): A tool to find structure
"""

# ╔═╡ b6478e1a-f5f6-11ea-3b92-6d4f067285f4
tree_url = "https://user-images.githubusercontent.com/6933510/110924885-d7f1b200-8322-11eb-9df7-7abf29c8db7d.png"

# ╔═╡ f2c11f88-f5f8-11ea-3e02-c1d4fa22031e
image = load(download(tree_url))

# ╔═╡ 29062f7a-f5f9-11ea-2682-1374e7694e32
picture = Float64.(channelview(image));

# ╔═╡ 5471fd30-f6e2-11ea-2cd7-7bd48c42db99
size(picture)

# ╔═╡ 6156fd1e-f5f9-11ea-06a9-211c7ab813a4
pr, pg, pb = eachslice(picture, dims=1)

# ╔═╡ a9766e68-f5f9-11ea-0019-6f9d02050521
[RGB.(pr, 0, 0) RGB.(0, pg, 0) RGB.(0, 0, pb)]

# ╔═╡ 0c0ee362-f5f9-11ea-0f75-2d2810c88d65
begin
	Ur, Σr, Vr = svd(pr)
	Ug, Σg, Vg = svd(pg)
	Ub, Σb, Vb = svd(pb)
end;

# ╔═╡ b95ce51a-f632-11ea-3a64-f7c218b9b3c9
@bind n Slider(1:200, show_value=true)

# ╔═╡ 7ba6e6a6-f5fa-11ea-2bcd-616d5a3c898b
RGB.(sum(outer(Ur[:,i], Vr[:,i]) .* Σr[i] for i in 1:n),
	 sum(outer(Ug[:,i], Vg[:,i]) .* Σg[i] for i in 1:n),
	 sum(outer(Ub[:,i], Vb[:,i]) .* Σb[i] for i in 1:n))

# ╔═╡ 8df84fcc-f5d5-11ea-312f-bf2a3b3ce2ce
md"# Appendix"

# ╔═╡ 0edd7cca-834f-11eb-0232-ff0850027f76
md"## Syntax Learned"

# ╔═╡ 69be8194-81b7-11eb-0452-0bc8b9f22286
md"""
Syntax to be learned:

* A `struct` is a great way to embody structure.
* `dump` and `Dump`: to see what's inside a data structure.
* `Diagonal`, `sparse`
* `error` (throws an exception)
* `svd` (Singular Value Decomposition)
"""

# ╔═╡ 1c462f68-834f-11eb-1447-85848814769b
[ dump, Diagonal, error, svd]

# ╔═╡ 5813e1b2-f5ff-11ea-2849-a1def74fc065
begin
	show_image(M) = get.([ColorSchemes.rainbow], M ./ maximum(M))
	show_image(x::AbstractVector) = show_image(x')
end

# ╔═╡ Cell order:
# ╟─0db6ee04-81b7-11eb-330c-11b578b72c90
# ╟─ca1a1072-81b6-11eb-1fee-e7df687cc314
# ╟─b0ba5b8c-f5d1-11ea-1304-3f0e47f935fe
# ╠═864e1180-f693-11ea-080e-a7d5aabc9ca5
# ╟─261c4df2-f5d2-11ea-2c72-7d4b09c46098
# ╟─3cada3a0-81cc-11eb-04c8-bde26d36a84e
# ╠═fe2028ba-f6dc-11ea-0228-938a81a91ace
# ╟─8d2c6910-f5d4-11ea-1928-1baf09815687
# ╟─54649792-81cc-11eb-1038-9161a4037acf
# ╠═0a902426-f6dd-11ea-0ae4-fb0c47863fe7
# ╟─4794e860-81b7-11eb-2c91-8561c20f308a
# ╟─67827da8-81cc-11eb-300e-278104d2d958
# ╠═4624cd26-f5d3-11ea-1cf8-7d6555eae5fa
# ╟─9bdabef8-81cc-11eb-14a1-67a9a7d968c0
# ╠═397ac764-f5fe-11ea-20cc-8d7cab19d410
# ╟─a22dcd2c-81cc-11eb-1252-13ace134192d
# ╠═82c7046c-f5d3-11ea-04e2-ef7c0f4db5da
# ╟─b024c318-81cc-11eb-018c-e1f7830ff51b
# ╠═93bfe3ac-f756-11ea-20fb-8f7d586b42f3
# ╠═175039aa-f758-11ea-251a-5db57d7c4b32
# ╠═c2a4b0a2-81cc-11eb-37a7-db601a6ddfdf
# ╟─c5ed7d3e-81cc-11eb-3386-15b72db8155d
# ╟─e2e354a8-81b7-11eb-311a-35151063c2a7
# ╟─dc5a96ba-81cc-11eb-3189-25920df48afa
# ╠═af0d3c22-f756-11ea-37d6-11b630d2314a
# ╟─06e3a842-26b8-4417-9cf5-8a083ccdb264
# ╠═91172a3e-81b7-11eb-0953-9f5e0207f863
# ╠═4bbf3f58-f788-11ea-0d24-6b0fb070829e
# ╟─fe70d104-81b7-11eb-14d0-eb5237d8ea6c
# ╟─ef8f44b2-f5fc-11ea-1e4d-bd873cd39d6c
# ╟─fd9211c0-f5fc-11ea-1745-7f2dae88af9e
# ╠═f1154df8-f693-11ea-3b16-f32835fcc470
# ╠═982590d4-f5ff-11ea-3802-73292c75ad6c
# ╟─81c35324-f5d4-11ea-2338-9f982d38732c
# ╟─2cfda0dc-f5d5-11ea-16c4-b5a33b90e37f
# ╠═150432d4-f5d5-11ea-32b2-19a2a91d9637
# ╟─44215aa4-f695-11ea-260e-b564c6fbcd4a
# ╠═21328d1c-f5d5-11ea-288e-4171ad35326d
# ╟─75761cc0-81cd-11eb-1186-7d47debd68ca
# ╟─6bd8a886-f758-11ea-2587-870a3fa9d710
# ╟─4c533ac6-f695-11ea-3724-b955eaaeee49
# ╠═466901ea-f5d5-11ea-1db5-abf82c96eabf
# ╠═b38c4aae-f5d5-11ea-39b6-7b0c7d529019
# ╟─93e04ed8-81cd-11eb-214a-a761ef8c406f
# ╟─e90c55fc-f5d5-11ea-10f1-470ff772985d
# ╟─19775c3c-f5d6-11ea-15c2-89618e654a1e
# ╟─653792a8-f695-11ea-1ae0-43761c502583
# ╠═79c94d2a-f75a-11ea-031d-09d70d229e15
# ╟─10bc5d50-81b9-11eb-2ac7-354a6c6c826b
# ╠═77d6a952-81ba-11eb-24e3-cb6510a59455
# ╟─1f3ba55a-81b9-11eb-001f-593b9d8639ca
# ╠═3d4a702e-f75a-11ea-031c-333d591fc442
# ╟─80ff4010-81bb-11eb-374e-215a57defb0b
# ╠═5de72b7c-f5d6-11ea-1b6f-35b830b5fb34
# ╠═8b60629e-f5d6-11ea-27c8-d934460d3a57
# ╠═2fd7e52e-f5d7-11ea-3b5a-1f338e2451e0
# ╠═2e87d4fe-81bc-11eb-0d16-b988bcedcc73
# ╠═cde79f38-f5d6-11ea-3297-0b5b240f7b9e
# ╠═aa09c008-f5d8-11ea-1bdc-b51ee6eb2478
# ╟─62a6ec62-f5d9-11ea-071e-ed33c5dea0cd
# ╟─7f63daf6-f695-11ea-0b80-8702a83103a4
# ╠═67274c3c-f5d9-11ea-3475-c9d228e3bd5a
# ╟─765c6552-f5d9-11ea-29d3-bfe7b4b04612
# ╠═126fb3ea-f5da-11ea-2f7d-0b3259a296ce
# ╟─ed0b2358-81ce-11eb-3339-93abcc06fd91
# ╟─24ce92fa-81cf-11eb-30f0-b1e357d79d53
# ╠═2d4500e0-81cf-11eb-1699-d310074fddf5
# ╟─3546ff30-81cf-11eb-3afc-05c5db61366f
# ╟─9b9e2c2a-f5da-11ea-369b-b513b196515b
# ╠═e68b98ea-f5da-11ea-1a9d-db45e4f80241
# ╠═f20ccac4-f5da-11ea-0e69-413b5e49f423
# ╠═12a2e96c-f5db-11ea-1c3e-494ae7446886
# ╠═22487ce2-f5db-11ea-32e9-6f70ab2c0353
# ╟─389ae62e-f5db-11ea-1557-c3adbbee0e5c
# ╟─0c2b6408-f5d9-11ea-2b7f-7fece2eecc1f
# ╠═5d767290-f5dd-11ea-2189-81198fd216ce
# ╠═587790ce-f6de-11ea-12d9-fde2a17ae314
# ╟─a39e8256-f6de-11ea-3170-c923b56609da
# ╠═8c84edd0-f6de-11ea-2180-61c6b81aac3b
# ╠═22b73baa-f6df-11ea-197f-bbb4bd1a7ef5
# ╠═b2332814-f6e6-11ea-1c7d-556c7d4687f1
# ╟─9ab7a72e-81cf-11eb-2b78-073ff51cae58
# ╟─fd8dd108-f6df-11ea-2f7c-3d99d054ac15
# ╟─165788b2-f601-11ea-3e69-cdbbb6558e54
# ╠═22941bb8-f601-11ea-1d6e-0d955297bc2e
# ╟─c33bf00e-81cf-11eb-1e1a-e5a879a45093
# ╠═2f75df7e-f601-11ea-2fc2-aff4f335af33
# ╟─7ff664f0-f74b-11ea-0d2d-b53f19e4f4bf
# ╟─a0611eaa-81bc-11eb-1d23-c12ab14138b1
# ╠═a4728944-f74b-11ea-03c3-9123908c1f8e
# ╠═05c35402-f752-11ea-2708-59cf5ef74fb4
# ╟─8c11b19e-81bc-11eb-184b-bf6ffefe29de
# ╠═8baaa994-f752-11ea-18d9-b3d0a6b9f7d9
# ╟─d92bece4-f754-11ea-0242-99f198bb5b7b
# ╠═e740999c-f754-11ea-2089-4b7a9aec6030
# ╟─0a79a7b4-f755-11ea-1b2d-21173567b257
# ╟─5adb98c2-f6e0-11ea-1fde-53b0fd6639c3
# ╟─487d6f9c-81d0-11eb-3bb0-336a4beb9b38
# ╟─5a493052-f601-11ea-2f5f-f940412905f2
# ╟─55b76aee-81d0-11eb-0bcc-413f5bd14360
# ╠═709bf30a-f755-11ea-2e82-bd511e598c77
# ╠═782532b0-f755-11ea-1385-cd1a28c4b9d5
# ╟─5bc4ab0a-f755-11ea-0fad-4987ad9fc02f
# ╠═a5d637ea-f5de-11ea-3b70-877e876bc9c9
# ╠═21bbb60a-f5df-11ea-2c1b-dd716a657df8
# ╠═2668e100-f5df-11ea-12b0-073a578a5edb
# ╠═e8d727f2-f5de-11ea-1456-f72602e81e0d
# ╠═f5fcdeea-f75c-11ea-1fc3-731f0ef1ad14
# ╠═0373fbf6-f75d-11ea-2a9e-cbb714d69cf4
# ╟─ebd72fb8-f5e0-11ea-0630-573337dff753
# ╠═b6478e1a-f5f6-11ea-3b92-6d4f067285f4
# ╠═f2c11f88-f5f8-11ea-3e02-c1d4fa22031e
# ╠═29062f7a-f5f9-11ea-2682-1374e7694e32
# ╠═5471fd30-f6e2-11ea-2cd7-7bd48c42db99
# ╠═6156fd1e-f5f9-11ea-06a9-211c7ab813a4
# ╠═a9766e68-f5f9-11ea-0019-6f9d02050521
# ╠═0c0ee362-f5f9-11ea-0f75-2d2810c88d65
# ╠═b95ce51a-f632-11ea-3a64-f7c218b9b3c9
# ╠═7ba6e6a6-f5fa-11ea-2bcd-616d5a3c898b
# ╟─8df84fcc-f5d5-11ea-312f-bf2a3b3ce2ce
# ╟─0edd7cca-834f-11eb-0232-ff0850027f76
# ╟─69be8194-81b7-11eb-0452-0bc8b9f22286
# ╠═1c462f68-834f-11eb-1447-85848814769b
# ╟─5813e1b2-f5ff-11ea-2849-a1def74fc065


text:
using Test
using LinearAlgebra
using CartesianMesh
using PolynomialBasis
using ImplicitDomainQuadrature
using HDGElasticity

function allapprox(v1, v2)
    return all(v1 .≈ v2)
end

function allapprox(v1, v2, atol)
    @assert length(v1) == length(v2)
    return all([isapprox(v1[i], v2[i], atol = atol) for i = 1:length(v1)])
end

function circle_distance_function(coords, xc, r)
    dist = [r - norm(coords[:, i] - xc) for i = 1:size(coords)[2]]
end

function plane_distance_function(coords,n,x0)
    return [n'*(coords[:,idx]-x0) for idx in 1:size(coords)[2]]
end

function bc_displacement(coords; alpha = 0.1, beta = 0.1)
    disp = copy(coords)
    disp[1, :] .*= alpha
    disp[2, :] .*= beta
    return disp
end

polyorder = 2
numqp = 6
Dhalf = sqrt(HDGElasticity.plane_strain_voigt_hooke_matrix_2d(1.0, 2.0))
stabilization = 0.1
normals = HDGElasticity.reference_normals()
levelset = InterpolatingPolynomial(1, 2, polyorder)
mesh = UniformMesh([0.0, 0.0], [1.0, 1.0], [1, 1])
circle_center = [2.5, 0.5]
radius = 2.0
levelsetcoeffs = HDGElasticity.levelset_coefficients(
    x -> circle_distance_function(x, circle_center, radius),
    mesh,
    levelset.basis,
)
# levelsetcoeffs = HDGElasticity.levelset_coefficients(
#     x -> plane_distance_function(x, [1.,0.], [0.5,0.]),
#     mesh,
#     levelset.basis,
# )

dgmesh = HDGElasticity.DGMesh(mesh, levelsetcoeffs, levelset)
ufs = HDGElasticity.UniformFunctionSpace(
    dgmesh,
    polyorder,
    numqp,
    levelsetcoeffs,
    levelset,
)
coords = HDGElasticity.nodal_coordinates(mesh,ufs.vbasis)

function curve_normals(imap,points)
    npts = length(points)
    normals = zeros(2,npts)
    for (idx,x) in enumerate(points)
        t = gradient(imap,x)
        n = [-t[2],t[1]]
        normals[:,idx] .= n/norm(n)
    end
    return normals
end

cellmap = HDGElasticity.CellMap(dgmesh.domain[1])

ufs.icoeffs[1][2] = -1.0
ufs.icoeffs[1][6] = +1.0
update!(ufs.imap, ufs.icoeffs[1])

update!(ufs.imap,ufs.icoeffs[1])
refnormals = curve_normals(ufs.imap,ufs.iquad.points)
invjac = HDGElasticity.inverse_jacobian(cellmap)
transformed_normals = diagm(invjac)*refnormals
magnitude = mapslices(x->norm(x),transformed_normals,dims=1)
inormals = -hcat([transformed_normals[:,i]/magnitude[i] for i = 1:size(magnitude)[2]]...)

LL = HDGElasticity.LLop(ufs.vbasis,ufs.vquads[1,1],cellmap)
LU = HDGElasticity.LUop(ufs.vbasis,ufs.vquads[1,1],Dhalf,cellmap)
fLH = HDGElasticity.LHop(ufs.vbasis,ufs.sbasis,
    ufs.fquads[1,1],dgmesh.facemaps,normals,Dhalf,cellmap)
LH = hcat(fLH...)
iLH = HDGElasticity.LHop_on_interface(ufs.vbasis,ufs.sbasis,ufs.iquad,
        ufs.imap,inormals,Dhalf,cellmap)
tLH = hcat(LH,iLH)
UU = HDGElasticity.UUop(ufs.vbasis,ufs.fquads[1,1],dgmesh.facemaps,
    ufs.iquad,inormals,ufs.imap,stabilization,cellmap)
fUH = HDGElasticity.UHop(ufs.vbasis,ufs.sbasis,ufs.fquads[1,1],
    dgmesh.facemaps,stabilization,cellmap)
UH = hcat(fUH...)
iUH = HDGElasticity.UHop_on_interface(ufs.vbasis,ufs.sbasis,ufs.iquad,
    ufs.imap,inormals,stabilization,cellmap)
tUH = hcat(UH,iUH)

HIrefcoords = reshape(ufs.icoeffs[1], 2, :)
HIcoords = hcat([cellmap(HIrefcoords[:, i]) for i = 1:size(HIrefcoords)[2]]...)
HIdisp = vec(bc_displacement(HIcoords))

H1c = [
    0.0 0.5 1.0
    0.0 0.0 0.0
]
H2c = [
    1.0 1.0 1.0
    0.0 0.5 1.0
]
H3c = [
    0.0 0.5 1.0
    1.0 1.0 1.0
]
H4c = [
    0.0 0.0 0.0
    0.0 0.5 1.0
]

Hcoords = [H1c, H2c, H3c, H4c]
isactiveface = [length(fq) > 0 ? true : false for fq in ufs.fquads[1, 1]]
faceids = findall(isactiveface)
Hfacedisp = vcat(vec.(bc_displacement.(Hcoords[faceids]))...)
Hdisp = vcat(Hfacedisp, HIdisp)
Udisp = vec(bc_displacement(coords))

lop = HDGElasticity.local_operator(
    ufs.vbasis,
    ufs.vquads[1, 1],
    ufs.fquads[1, 1],
    dgmesh.facemaps,
    ufs.iquad,
    ufs.imap,
    -inormals,
    Dhalf,
    stabilization,
    cellmap,
)
facelhops = HDGElasticity.local_hybrid_operator(
    ufs.vbasis,
    ufs.sbasis,
    ufs.fquads[1, 1],
    dgmesh.facemaps,
    normals,
    Dhalf,
    stabilization,
    cellmap,
)
lhop = hcat(facelhops...)
ilhop = HDGElasticity.local_hybrid_operator_on_interface(
    ufs.vbasis,
    ufs.sbasis,
    ufs.iquad,
    ufs.imap,
    -inormals,
    Dhalf,
    stabilization,
    cellmap,
)

HIrefcoords = reshape(ufs.icoeffs[1], 2, :)
HIcoords = hcat([cellmap(HIrefcoords[:, i]) for i = 1:size(HIrefcoords)[2]]...)
HIdisp = vec(bc_displacement(HIcoords))

H1c = [
    0.0 0.5 1.0
    0.0 0.0 0.0
]
H2c = [
    1.0 1.0 1.0
    0.0 0.5 1.0
]
H3c = [
    0.0 0.5 1.0
    1.0 1.0 1.0
]
H4c = [
    0.0 0.0 0.0
    0.0 0.5 1.0
]

lochyb = hcat(lhop, ilhop)

rhs = lochyb * Hdisp

sol = lop \ rhs
L = -Dhalf * reshape(sol[1:27], 3, :)
U = reshape(sol[28:45], 2, :)

# testL = zeros(3, 9)
# testL[1:2, :] .= 0.6
# @test allapprox(L, testL, 1e-10)


text:
module SimpleHypergraphs

using LightGraphs

export Hypergraph, getvertices, gethyperedges, hg_load, hg_save
export add_vertex!, add_hyperedge!
export set_vertex_meta!, get_vertex_meta
export set_hyperedge_meta!, get_hyperedge_meta
export BipartiteView, shortest_path
export TwoSectionView

include("hypergraph.jl")
include("bipartite.jl")
include("io.jl")
include("twosection.jl")

end # module


text:
using DataStructures

"""
    counter is used to make conditional distribution, which is used by score functions to
    calculate conditonal frequency distribution
"""
function counter2(data, min::Integer, max::Integer)
    data = everygram(data, min_len=min, max_len=max)
    data = split.(data)
    temp_lm = DefaultDict{SubString{String}, Accumulator{String,Int64}}(counter(SubString{String}))
    for i in 1:length(data)
        history,word = data[i][1:end-1], data[i][end]
        temp_lm[join(history, " ")][word] += 1
    end
    return temp_lm
end


text:
# ---
# title: 385. Mini Parser
# id: problem385
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: String, Stack
# link: <https://leetcode.com/problems/mini-parser/description/>
# hidden: true
# ---
#
# Given a nested list of integers represented as a string, implement a parser to
# deserialize it.
#
# Each element is either an integer, or a list -- whose elements may also be
# integers or other lists.
#
# **Note:** You may assume that the string is well-formed:
#
#   * String is non-empty.
#   * String does not contain white spaces.
#   * String contains only digits `0-9`, `[`, `-` `,`, `]`.
#
#
#
# **Example 1:**
#
#
#
#     Given s = "324",
#
#     You should return a NestedInteger object which contains a single integer 324.
#
#
#
#
# **Example 2:**
#
#
#
#     Given s = "[123,[456,[789]]]",
#
#     Return a NestedInteger object containing a nested list with 2 elements:
#
#     1. An integer containing value 123.
#     2. A nested list containing two elements:
#         i.  An integer containing value 456.
#         ii. A nested list with one element:
#              a. An integer containing value 789.
#
#
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
import GPUArrays: allowscalar, @allowscalar

function _getindex(xs::CuArray{T}, i::Integer) where T
  buf = Mem.view(buffer(xs), (i-1)*sizeof(T))
  return Mem.download(T, buf)[1]
end

function _setindex!(xs::CuArray{T}, v::T, i::Integer) where T
  buf = Mem.view(buffer(xs), (i-1)*sizeof(T))
  Mem.upload!(buf, T[v])
end

## logical indexing

Base.getindex(xs::CuArray, bools::AbstractArray{Bool}) = getindex(xs, CuArray(bools))

function Base.getindex(xs::CuArray{T}, bools::CuArray{Bool}) where {T}
  bools = reshape(bools, prod(size(bools)))
  indices = cumsum(bools)  # unique indices for elements that are true

  n = _getindex(indices, length(indices))  # number that are true
  ys = CuArray{T}(undef, n)

  if n > 0
    num_threads = min(n, 256)
    num_blocks = ceil(Int, length(indices) / num_threads)

    function kernel(ys::CuDeviceArray{T}, xs::CuDeviceArray{T}, bools, indices)
        i = threadIdx().x + (blockIdx().x - 1) * blockDim().x

        if i <= length(xs) && bools[i]
            b = indices[i]   # new position
            ys[b] = xs[i]

        end

        return
    end

    @cuda blocks=num_blocks threads=num_threads kernel(ys, xs, bools, indices)
  end

  return ys
end


text:
using Test

STARTING_SUBJECT_NUM = 7
MAX_VAL = 20201227 # ~20M

function transform_subject_number(subject_number, loop_size)
    value = 1
    for i in 1:loop_size
        # Set the value to itself multiplied by the subject number.
        value = value * subject_number
        # Set the value to the remainder after dividing the value by 20201227.
        value = mod(value, MAX_VAL)
    end
    return value
end

function search_for_loop_size(target)
    value = 1
    for i in 1:MAX_VAL
        # Set the value to itself multiplied by the subject number.
        value = value * STARTING_SUBJECT_NUM
        # Set the value to the remainder after dividing the value by 20201227.
        value = mod(value, MAX_VAL)
        if value == target
            return i
        end
    end
    throw("failed to find loop size")
end

function run(public_keys)
    # Discover loop size
    loop_sizes = map(search_for_loop_size, public_keys)

    @show public_keys
    @show loop_sizes

    # The card transforms the subject number of the door's public key according
    # to the card's loop size. The result is the encryption key.

    # The door transforms the subject number of the card's public key according
    # to the door's loop size. The result is the same encryption key as the
    # card calculated.
    encryption_keys = [
        transform_subject_number(public_keys[2], loop_sizes[1])
        transform_subject_number(public_keys[1], loop_sizes[2])
    ]
    @assert encryption_keys[1] == encryption_keys[2]

    result = encryption_keys[1]
    println("Result = $(result)")
    return result
end

@test run([5764801,17807724]) == 14897079

run([1965712, 19072108])


text:
__precompile__(true)
module AutoEnvs

using AutoRisk
using AutoViz
using HDF5
using JLD
using NGSIM
using PyPlot

import AutoViz: render
import Base: reset, step

# module
include("make.jl")
include("env.jl")
include("debug_envs.jl")
include("ngsim_utils.jl")
include("ngsim_env.jl")
include("vectorized_ngsim_env.jl")
include("multiagent_ngsim_env.jl")
include("multiagent_ngsim_env_videomaker.jl")
end


text:
#-------------------------------------------------------------------------------
# Set of useful functions related to LOVO functions.
#-------------------------------------------------------------------------------

"""

    fmin_eval!(func_list::Array{Function, 1}, r:: Int64,
                y::Vector{Float64}, imin_set::Vector{Bool})

Computes the value of the objective function fmin(y) and an index belonging to
the set I_{min}(y).

    - 'func_list': list containing the functions that determine the objective
    function fmin.

    - 'r': number of functions that make up the objective function fmin.

    - 'y': n-dimensional vector.

The function modifies the argument:

    - 'imin_set': boolean vector with the indexes belonging to the I_{min}(y) set.

Returns the function value 'fmin_y' and the index 'imin_y'.

"""
function fmin_eval!(
                    func_list::Array{Function, 1},
                    r::Int64,
                    y::Vector{Float64},
                    imin_set::Vector{Bool}
                    )

    fmin_y = func_list[1](y)
    imin_y = 1
    imin_set[1] = true

    for i = 2:r

        tmp = func_list[i](y)

        if tmp < fmin_y

            fmin_y = tmp
            imin_y = i
            imin_set[1:(i - 1)] .= false
            imin_set[i] = true

        elseif tmp == fmin_y

            imin_set[i] = true

        else

            imin_set[i] = false

        end

    end

    return fmin_y, imin_y

end

"""

    fmin_partial_eval!(func_list::Array{Function, 1}, r:: Int64, idx::Int64,
                        fi_y::Float64, y::Vector{Float64}, imin_set::Vector{Bool})

Computes the value of the objective function fmin(y) and an index belonging to
the set Imin(y), using the precalculated value of f_idx(y).

    - 'func_list': list containing the functions that determine the objective
    function fmin.

    - 'r': number of functions that make up the objective function fmin.

    - 'idx': index corresponding to the available function value.

    - 'fi_y': available function value.

    - 'y': n-dimensional vector.

The function modifies the argument:

    - 'imin_set': boolean vector with the indexes belonging to the I_{min}(y) set.

Returns the function value 'fmin_y' and the index 'imin_y'.

"""
function fmin_partial_eval!(
                    func_list::Array{Function, 1},
                    r::Int64,
                    idx::Int64,
                    fi_y::Float64,
                    y::Vector{Float64},
                    imin_set::Vector{Bool}
                    )

    fmin_y = fi_y
    imin_y = idx
    imin_set[idx] = true

    for i = 1:r

        if i == idx

            tmp = fi_y

        else

            tmp = func_list[i](y)

        end

        if tmp < fmin_y

            fmin_y = tmp
            imin_y = i
            imin_set[1:(i - 1)] .= false
            imin_set[i] = true

        elseif tmp == fmin_y

            imin_set[i] = true

        else

            imin_set[i] = false

        end

    end

    return fmin_y, imin_y

end

"""

    fi_eval(func_list::Array{Function, 1}, idx:: Int64,
                y::Vector{Float64})

Computes the value of the function f_index(y).

    - 'func_list': list containing the functions that determine the objective
    function fmin.

    - 'idx': index of the funtion in 'func_list'.

    - 'y': n-dimensional vector.

Returns the function value.

"""
function fi_eval(
                    func_list::Array{Function, 1},
                    idx:: Int64,
                    y::Vector{Float64}
                    )

    return func_list[idx](y)

end

max_stars_repo_name: aschwertner/LOWDER

max_issues_repo_name: aschwertner/LOWDER

max_forks_repo_name: aschwertner/LOWDER


text:
if basename(pwd()) == "aoc"
    cd("2018/9")
end

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
# Configures plotting backend to use for demos
#-------------------------------------------------------------------------------

using CMDimData
using CMDimData.EasyPlot
CMDimData.@includepkg EasyPlotInspect
EasyPlot.defaults.guibuilder = EasyPlot.getbuilder(:gui, :InspectDR)

# last line


text:
using RandomNumbers
using Test

@testset "Check types" begin
    @test typeof(randfloat()) == Float64
    @test typeof(randfloat(Float32)) == Float32
    @test typeof(randfloat(Float64)) == Float64
    @test typeof(randfloat(Float16)) == Float16

    # arrays
    @test typeof(randfloat(Float16,2)) <: Vector{Float16}
    @test typeof(randfloat(Float32,2)) <: Vector{Float32}
    @test typeof(randfloat(Float64,2)) <: Vector{Float64}

    @test typeof(randfloat(Float16,2,3)) <: Matrix{Float16}
    @test typeof(randfloat(Float32,2,3)) <: Matrix{Float32}
    @test typeof(randfloat(Float64,2,3)) <: Matrix{Float64}
end

@testset "Non-default RNG" begin
    import RandomNumbers.Xorshifts.Xorshift64
    rng = Xorshift64()
    @test typeof(randfloat(rng)) == Float64
    @test typeof(randfloat(rng,Float16)) == Float16
    @test typeof(randfloat(rng,Float32)) == Float32
    @test typeof(randfloat(rng,Float64)) == Float64

    # arrays
    @test typeof(randfloat(rng,Float16,2)) <: Vector{Float16}
    @test typeof(randfloat(rng,Float32,2)) <: Vector{Float32}
    @test typeof(randfloat(rng,Float64,2)) <: Vector{Float64}

    @test typeof(randfloat(rng,Float16,2,3)) <: Matrix{Float16}
    @test typeof(randfloat(rng,Float32,2,3)) <: Matrix{Float32}
    @test typeof(randfloat(rng,Float64,2,3)) <: Matrix{Float64}

    import RandomNumbers.MersenneTwisters.MT19937
    rng = MT19937()
    @test typeof(randfloat(rng)) == Float64
    @test typeof(randfloat(rng,Float16)) == Float16
    @test typeof(randfloat(rng,Float32)) == Float32
    @test typeof(randfloat(rng,Float64)) == Float64

    # arrays
    @test typeof(randfloat(rng,Float16,2)) <: Vector{Float16}
    @test typeof(randfloat(rng,Float32,2)) <: Vector{Float32}
    @test typeof(randfloat(rng,Float64,2)) <: Vector{Float64}

    @test typeof(randfloat(rng,Float16,2,3)) <: Matrix{Float16}
    @test typeof(randfloat(rng,Float32,2,3)) <: Matrix{Float32}
    @test typeof(randfloat(rng,Float64,2,3)) <: Matrix{Float64}
end

@testset "Distribution is uniform [0,1) Float64" begin
    N = 100_000_000

    x = randfloat(N)
    @test maximum(x) > 0.999999
    @test minimum(x) < 1e-7

    exponent_mask = 0x7ff0000000000000
    exponent_bias = 1023

    # histogram for exponents
    H = fill(0,64)
    for xi in x     # 0.5 is in H[1], 0.25 is in H[2] etc
        e = reinterpret(UInt64,xi) & exponent_mask
        H[exponent_bias-Int(e >> 52)] += 1
    end

    # test that the most frequent exponents occur at 50%, 25%, 12.5% etc.
    for i in 1:10
        @test isapprox(H[i]/N,2.0^-i,atol=5e-4)
    end
end

@testset "Distribution is uniform [0,1) Float32" begin
    N = 100_000_000

    x = randfloat(Float32,N)
    @test maximum(x) > prevfloat(1f0,5)
    @test minimum(x) < 1e-7

    exponent_mask = 0x7f800000
    exponent_bias = 127

    # histogram for exponents
    H = fill(0,64)
    for xi in x     # 0.5f0 is in H[1], 0.25f0 is in H[2] etc
        e = reinterpret(UInt32,xi) & exponent_mask
        H[exponent_bias-Int(e >> 23)] += 1
    end

    # test that the most frequent exponents occur at 50%, 25%, 12.5% etc.
    for i in 1:10
        @test isapprox(H[i]/N,2.0^-i,atol=5e-4)
    end
end

@testset "Distribution is uniform [0,1) Float16" begin
    N = 100_000_000

    x = randfloat(Float16,N)
    @test maximum(x) == prevfloat(one(Float16))
    @test minimum(x) == zero(Float16)               # zero can be produced

    exponent_mask = 0x7c00
    exponent_bias = 15

    # histogram for exponents
    H = fill(0,16)
    for xi in x     # 0.5f0 is in H[1], 0.25f0 is in H[2] etc
        e = reinterpret(UInt16,xi) & exponent_mask
        H[exponent_bias-Int(e >> 10)] += 1
    end

    # test that the most frequent exponents occur at 50%, 25%, 12.5% etc.
    for i in 1:10
        @test isapprox(H[i]/N,2.0^-i,atol=5e-4)
    end
end


text:
# Sample usage

using Revise
using TypeDBClient: dbconnect, open, read, write, match, insert, commit, create_database

dbconnect("127.0.0.1") do client
    create_database(client, "my-typedb")
    open(client, "my-typedb") do session

        write(session) do transaction
            insert(transaction, raw"insert $_ isa person;")
            commit(transaction)
        end
        read(session) do transaction
            answers = match(transaction, raw"match $p isa person;")
        end
    end
end


text:
using CollectiveSpins, QuantumOptics
using Test

@testset "reducedspin" begin

# Test creation
N = 4
b_full = SpinBasis(1//2)^N
b_red1 = CollectiveSpins.ReducedSpinBasis(N, 1)
b_red2 = CollectiveSpins.ReducedSpinBasis(N, 2)

# Operators
σp = [embed(b_full, i, sigmap(SpinBasis(1//2))) for i=1:N]

# Get all single excitation and two excitation states
GS_full = tensor([spindown(SpinBasis(1//2)) for i=1:N]...)
single_exc = Ket[GS_full]
two_exc = Ket[]
for i=1:N
    psi_ = σp[i]*GS_full
    push!(single_exc, normalize(psi_))
    for j=i+1:N
        push!(two_exc, normalize(σp[j]*psi_))
    end
end

# Subspace and projectors
b_sub1 = SubspaceBasis(b_full, single_exc)
b_sub2 = SubspaceBasis(b_full, [single_exc; two_exc])
@test length(b_red1) == length(b_sub1)
@test length(b_red2) == length(b_sub2)

P1 = projector(b_sub1, b_full)
P2 = projector(b_sub2, b_full)

psi1 = P1*single_exc[2]
psi2 = CollectiveSpins.reducedspinstate(b_red1, [1])
@test psi1 != psi2
@test psi1.data == psi2.data

psi1 = P2*two_exc[end] # Corresponds to state where N-1 and N are excited
psi2 = CollectiveSpins.reducedspinstate(b_red2, [N-1,N])
@test psi1 != psi2
@test psi1.data == psi2.data

# Test matrix elements on reduced bases
sz(i) = embed(b_full, i, sigmaz(SpinBasis(1//2)))
sp(i) = embed(b_full, i, sigmap(SpinBasis(1//2)))
for i=1:N
    @test CollectiveSpins.reducedsigmap(b_red1, i).data == sparse(P1*sp(i)*P1').data
    @test CollectiveSpins.reducedsigmam(b_red1, i).data == sparse(P1*sp(i)'*P1').data
    @test CollectiveSpins.reducedsigmap(b_red2, i).data == sparse(P2*sp(i)*P2').data
    @test CollectiveSpins.reducedsigmam(b_red2, i).data == sparse(P2*sp(i)'*P2').data
    @test CollectiveSpins.reducedsigmaz(b_red1, i).data == sparse(P1*sz(i)*P1').data
    @test CollectiveSpins.reducedsigmaz(b_red2, i).data == sparse(P2*sz(i)'*P2').data
end

i=2
j=4
tmp = CollectiveSpins.reducedsigmap(b_red2,i)*CollectiveSpins.reducedsigmap(b_red2,j)*CollectiveSpins.reducedspinstate(b_red2,[])
ind = findfirst(!iszero, tmp.data)
@test b_red2.indexMapper[findfirst(y->y[2]==ind, b_red2.indexMapper)][1] == [i,j]

# Test transitions
@test CollectiveSpins.reducedspintransition(b_red2, [i,j], [i-1,j-1]) == CollectiveSpins.reducedspintransition(b_red2, [i,j], [j-1,i-1])
@test_throws BoundsError CollectiveSpins.reducedspintransition(b_red2, [i+1,j+1], [1,1])

# Test a full example
N = 8     # Number of spins
M = 3     # Maximum number of excitations
a = 0.3   # spin-spin distance
eta = 0.5 # Drive amplitude
geometry = CollectiveSpins.geometry.chain(a, N)

e_list = [[0,0,1] for i=1:N]   # Quantization axis
S = CollectiveSpins.SpinCollection(geometry, e_list)

b = CollectiveSpins.ReducedSpinBasis(N, M)

# Drive
sx(j) = CollectiveSpins.reducedsigmax(b, j)
H_drive = eta*sum(sx(j) for j=1:N)

# Dipole-Dipole
H_dipole = CollectiveSpins.reducedspin.Hamiltonian(S, M)
Gammas, Jumps = CollectiveSpins.reducedspin.JumpOperators(S, M)

# Time Evolution
T = [0:0.05:3.;]
psi0 = CollectiveSpins.reducedspinstate(b, []) # start in the ground state
tout, rho = QuantumOptics.timeevolution.master(T, psi0, H_drive + H_dipole, Jumps; rates=Gammas)

# Test sigmaz
N = 4
M = 4
b_red = CollectiveSpins.ReducedSpinBasis(N,M)
SZ = sum(CollectiveSpins.reducedsigmaz(b_red,j) for j=1:N)

@test sum(SZ.data)==0
sz_diag = [SZ.data[i,i] for i=1:2^N]
for m=0:N
    inds = findall(isequal(-N+2m),sz_diag)
    @test length(inds) == binomial(N,m)
    state = CollectiveSpins.reducedspinstate(b_red, [1:m;])
    @test expect(SZ,state) == -N+2m
end

# Test with lower bound
b_red1 = CollectiveSpins.ReducedSpinBasis(N, 1, 1)
b_red2 = CollectiveSpins.ReducedSpinBasis(N, 2, 1)

# Subspace and projectors
tmp = copy(single_exc)
single_exc = tmp[2:end] # Remove ground state
b_sub1 = SubspaceBasis(b_full, single_exc)
b_sub2 = SubspaceBasis(b_full, [single_exc; two_exc])
@test length(b_red1) == length(b_sub1)
@test length(b_red2) == length(b_sub2)

P1 = projector(b_sub1, b_full)
P2 = projector(b_sub2, b_full)

psi1 = P1*single_exc[1]
psi2 = CollectiveSpins.reducedspinstate(b_red1, [1])
@test psi1 != psi2
@test psi1.data == psi2.data

psi1 = P2*two_exc[end] # Corresponds to state where N-1 and N are excited
psi2 = CollectiveSpins.reducedspinstate(b_red2, [N-1,N])
@test psi1 != psi2
@test psi1.data == psi2.data

# Test matrix elements on reduced bases
sz(i) = embed(b_full, i, sigmaz(SpinBasis(1//2)))
sp(i) = embed(b_full, i, sigmap(SpinBasis(1//2)))
for i=1:N
    @test CollectiveSpins.reducedsigmap(b_red2, i).data == sparse(P2*sp(i)*P2').data
    @test CollectiveSpins.reducedsigmam(b_red2, i).data == sparse(P2*sp(i)'*P2').data
    @test CollectiveSpins.reducedsigmaz(b_red1, i).data == sparse(P1*sz(i)*P1').data
    @test CollectiveSpins.reducedsigmaz(b_red2, i).data == sparse(P2*sz(i)'*P2').data
end

SZ1 = sum(CollectiveSpins.reducedsigmaz(b_red1, i) for i=1:N)
@test all([SZ1.data[i,i] for i=1:N] .== -N+2)
SZ2 = sum(CollectiveSpins.reducedsigmaz(b_red2, i) for i=1:N)
@test all([SZ2.data[i,i] for i=1:length(b_red2)] .== [[-2.0 for i=1:N]; zeros(binomial(N,2))])

end #testset