text:
abstract type AbstractMolecule end

struct N₂ <: AbstractMolecule end
struct O₂ <: AbstractMolecule end
struct H₂ <: AbstractMolecule end

"""
    struct PolTensor{FT}

A struct which provides all polarizability tensor elements
    relevant to Rayleigh, and rotational+vibrational Raman cross-sections
    Ref: Buldakov et al.,(1999), Molecular Spectroscopy

    The polarizability at arbitrary wavenumbers [m⁻¹] and temperatures [K] can be computed as

    α̅(2πcν, T) = α̅₀₀(1 + α_b⋅T + α_c⋅T²)/(1-(2πcν/ω₀)²

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct PolarizationTensor{FT}
    "Average Polarizability `[m³]` at an angular frequency of ω₀ `[Hz]` and T=0K"
    α̅₀₀::FT
    "First derivative of `α̅₀₀` wrt the internuclear distance rₑ `[m²]`"
    α₀₀_prime::FT
    "Reference frequency `[Hz]`"
    ω₀::FT
    "Linear T-dependence coefficient `[K⁻¹]`"
    α_b::FT
    "Quadratic T-dependence coefficient `[K⁻²]`"
    α_c::FT
    "Polarizability anisotropy `[m³]` "
    γ̅₀₀ ::FT
    "First derivative of `γ̅` with respect to the internuclear distance `rₑ` `[m²]`"
    γ₀₀_prime::FT
end

"""
    struct DunhamCoefficients{FT}

A struct which provides
the nuclear spin multiplicity for a (currently diatomic homonuclear) molecule for odd
and even values of initial rotational quantum states J_initial

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct NuclearSpinMultiplicity
    "Nuclear Spin Multiplicity (Array of length 2 {odd J_initial, even J_initial})"
    gₛ::AbstractArray{Int}
end

"""
    struct DunhamCoefficients{FT}

A struct which provides
the Dunham expansion coefficients to compute the molecular energy levels
corresponding to different rotational (J) and vibrational (v) states. The overall state
of the molecule is given by its electronic ground state. #Constants from "Molecular
Spectra and Molecular Structure IV. Constants of diatomic molecules" by K.P. Huber
and G. Herzberg, 1978 (https://webbook.nist.gov/cgi/cbook.cgi?ID=C7727379&Mask=1000)

Needs description of Y elements

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct DunhamCoefficients{FT}
    "Dunham Expansion Coefficients (5x5 Matrix)"
    Y::AbstractArray{FT}
end

Base.@kwdef mutable struct EffectiveCoefficients{FT}
    "Ambient Temperature `[K]` at which following effective coefficients are computed"
    T::FT
    "Average Polarizability `[m³]` at wavenumber ν and temperature T"
    α̅::FT      #  Assumes a central wavelength effective over a broadband range of wavelengths
    γ̅::FT
    α_prime::FT
    γ_prime::FT
    #"King correction factor, See Table 1 and Eq 4 of Chance & Spurr (1997)"
    #F_King::FT #  See Table 1 and Eq 4 of Chance & Spurr (1997)
    "α/γ"
    ϵ::FT
    #"King correction factor for vibrational Raman scattering, see Tables H, I (CRS 13, D. A. Long (1977)) "
    #F_King_prime::FT
    "α_prime/γ_prime"
    ϵ_prime::FT
    γ_C_Rayl::FT # 3/(45(ϵ)^2+4)
    γ_C_RotRaman::FT # 3/4
    γ_C_VibRaman::FT # 3/(45(ϵ_prime)^2+4)
    γ_C_RoVibRaman::FT # 3/4
    rho_depol_Rayl::FT # 2γ_C_Rayl/(1+γ_C_Rayl)
    rho_depol_RotRaman::FT # 2γ_C_RotRaman/(1+γ_C_RotRaman)
    rho_depol_VibRaman::FT # 2γ_C_VibRaman/(1+γ_C_VibRaman)
    rho_depol_RoVibRaman::FT # 2γ_C_RoVibRaman/(1+γ_C_RoVibRaman)
    "Rayl cross-section σ(ν) = σ_Rayl_coeff * (ν)⁴"
    σ_Rayl_coeff::FT  #Cross-section = σ_Rayl_coeff * ν⁴
    σ_Rayl_coeff_hires::AbstractArray{FT}
    Δν̃_Rayl_coeff_hires::AbstractArray{FT}

    #"Wavenumber shift in emitted light after vᵢ, Jᵢ --> v_f, J_f transition, νₛ = ν + Δνₛ"
    #Δνₛ::AbstractArray{FT}
    "Raman scattering coefficient for pure vibrational scattering"
    σ_VibRaman_coeff_0to1::FT
    σ_VibRaman_coeff_1to0::FT
    σ_VibRaman_coeff_0to1_hires::AbstractArray{FT}
    σ_VibRaman_coeff_1to0_hires::AbstractArray{FT}

    "Vibrational Raman scattering wavenumber displacement Δν̃"
    Δν̃_VibRaman_coeff_0to1::FT
    Δν̃_VibRaman_coeff_1to0::FT
    Δν̃_VibRaman_coeff_0to1_hires::AbstractArray{FT}
    Δν̃_VibRaman_coeff_1to0_hires::AbstractArray{FT}

    "Raman scattering coefficient for rovibrational scattering"
    σ_RoVibRaman_coeff_0to1_JtoJm2::AbstractArray{FT}
    σ_RoVibRaman_coeff_0to1_JtoJp2::AbstractArray{FT}
    σ_RoVibRaman_coeff_1to0_JtoJm2::AbstractArray{FT}
    σ_RoVibRaman_coeff_1to0_JtoJp2::AbstractArray{FT}  #Cross-section = σ_VibRaman_coeff * ν⁴
    "Rovibrational wavenumber displacement Δν̃"
    Δν̃_RoVibRaman_coeff_0to1_JtoJm2::AbstractArray{FT}
    Δν̃_RoVibRaman_coeff_0to1_JtoJp2::AbstractArray{FT}
    Δν̃_RoVibRaman_coeff_1to0_JtoJm2::AbstractArray{FT}
    Δν̃_RoVibRaman_coeff_1to0_JtoJp2::AbstractArray{FT}
    "Rotational Raman cross-section σ(ν, vᵢ, Jᵢ --> v_f, J_f) = σ_Raman_coeff * (ν+Δνₛ)⁴"
    σ_RoRaman_coeff_JtoJm2::AbstractArray{FT}
    σ_RoRaman_coeff_JtoJp2::AbstractArray{FT} #Cross-section = σ_Raman_coeff * ν⁴
    "Rotational Raman wavenumber displacement Δν̃"
    Δν̃_RoRaman_coeff_JtoJm2::AbstractArray{FT}
    Δν̃_RoRaman_coeff_JtoJp2::AbstractArray{FT} #Cross-section = σ_Raman_coeff * ν⁴
    "Energy levels (in wavenumbers [cm^{-1}]) for a given v, J state"
    E_vJ::AbstractArray{FT}
end

Base.@kwdef struct MolecularConstants{FT}
    "Volume mixing ratio"
    vmr::FT
    "Dunham Expansion Coefficients (5x5 Matrix)"
    Y::AbstractArray{FT}
    "Nuclear Spin Multiplicity (Array of length 2 {odd J_initial, even J_initial})"
    gₛ::AbstractArray{Int}
    PolTensor::PolarizationTensor{FT}
    effCoeff::EffectiveCoefficients{FT}
end

max_stars_repo_name: RadiativeTransfer/InelasticScattering

max_issues_repo_name: RadiativeTransfer/InelasticScattering

max_forks_repo_name: RadiativeTransfer/InelasticScattering


text:
@testset "Random" begin
  Qx, x = PolynomialRing(FlintQQ, "x")
  K, a = NumberField(x^32 + 2, "a")

  b = @inferred rand([a], -10:10)
  @test b isa nf_elem
  @test_throws ErrorException rand(nf_elem[], -10:10)

  b = @inferred rand(basis(K), 1:100, 10)
  @test count(!iszero, (coeff(b, i) for i in 0:31)) <= 10
  @test_throws ErrorException rand(nf_elem[], -10:10, 5)
  @test_throws ErrorException rand([a, a^2], -10:10, -10)
  @test_throws ErrorException rand(basis(K), -10:10, 100)

  @inferred rand!(b, basis(K), 1:100, 20)
  @test count!(iszero, (Coeff(b, i) for i in 0:31)) <= 20
  @test_throws ErrorException rand!(b, basis(K), 1:100, 120)
  @test_throws ErrorException rand!(b, basis(K), 1:100, -100)

  @inferred rand!(b, basis(K), 1:100)
  @test_throws ErrorException rand!(b, nf_elem[], 1:100)
end

ext: jl


text:
using Revise
using Random
using LinearAlgebra
using POMDPs
using POMDPModelTools
using POMDPPolicies
using POMDPSimulators
using IterTools
using AutomotiveDrivingModels
using AutomotivePOMDPs: polygon, overlap
using Parameters
using GridInterpolations
using StaticArrays
using AutoUrban
using AutoViz
using Printf
using Reel
using DiscreteValueIteration
using SparseArrays
using ProgressMeter
using JLD2
using FileIO

includet("environment.jl")
includet("mdp_definition.jl")
includet("helpers.jl")

mdp = MergingMDP()

@sprintf("Number of states: %i ", n_states(mdp))

# rm("transition.jld2")
tranfile = "/scratch/boutonm/transition_cv.jld2"
trans_mat = DiscreteValueIteration.transition_matrix_a_s_sp(mdp, filename=tranfile)
JLD2.@save tranfile trans_mat
# rm("reward.jld2")
rewfile = "/scratch/boutonm/reward_cv.jld2"
reward_mat = DiscreteValueIteration.reward_s_a(mdp, filename=rewfile)
JLD2.@save rewfile reward_mat

solver = SparseValueIterationSolver(max_iterations=300, verbose=true)
# solver = ValueIterationSolver(max_iterations=30, verbose=true)
policy = solve(solver, mdp)

# post - process policy
for (i, s) in enumerate(states(mdp))
    if isterminal(mdp, s) && !collisioncheck(mdp, s)
        policy.util[i] = 1.0
        policy.qmat[i, :] .= 1.0
    end
end

# @save "policy_wide.jld2" policy
# @save "/scratch/boutonm/policy_cv.jld2" policy

@load "policy_wide.jld2" policy
policy_wide = deepcopy(policy)
@load "policy_cv.jld2" policy
policy_cv = deepcopy(policy)

## Evaluate policy

rng = MersenneTwister(1)
randpol = RandomPolicy(mdp, rng=rng)
s0 = [mdp.ego_grid[1], mdp.velocity_grid[1], mdp.acceleration_grid[3], mdp.human_grid[1], mdp.velocity_grid[5]]

simlist = [Sim(mdp, policy, s0, metadata=Dict(:policy=>"VI policy"),
               rng=MersenneTwister(i), max_steps=200) for i=1:100]

res = run_parallel(simlist)

@sprintf(" Probability of success at the initial state: %1.4f ", value(policy, s0))
println(" Evaluation of $(res[:policy][1]) ")
@sprintf(" Average reward: %2.4f ", mean(res[:reward]))
@sprintf(" Maximum reward: %2.2f ", maximum(res[:reward]))

using Plots
using Interact
using Blink

w = Window()
hr = HistoryRecorder(rng = rng, max_steps=200)
hist = simulate(hr, mdp, policy, s0)
ui = @manipulate for step in 1:n_steps(hist)
    s = state_hist(hist)[step]
    a = collect(action_hist(hist))[step]
    v = MergingViz(mdp=mdp, s=s, a=a, action_values = actionvalues(policy, s))
    AutoViz.render([v], cam=StaticCamera(VecE2(0.0, 0.0), 15.0))
end
body!(w, ui)

ui = @manipulate for ve = mdp.velocity_grid, vo = mdp.velocity_grid
    xs = mdp.ego_grid[1:end]
    ys = mdp.human_grid
    ae = mdp.acceleration_grid[3]
    z = map(((y,x),) -> value(policy, Float64[x, ve, ae, y, vo]), Iterators.product(ys, xs))
    za = map(((y,x),) -> action(policy, Float64[x, ve, ae, y, vo]), Iterators.product(ys, xs))
    plot(heatmap(xs, ys, z, aspect_ratio=1, xlabel="Ego longitudinal position", ylabel="Other longitudinal position", title="Pr of success", color=:plasma, clim=(0,1)),
        heatmap(xs, ys, za, aspect_ratio=1, title="Action Map", color=:viridis, clim=(1,7)),
        layout=(1,2))
end
body!(w, ui)

for policy in [policy_wide, policy_cv]
    policy.qmat = round.(policy_wide.qmat, digits=2)
    policy.util = round.(policy.util, digits=2)
end

# w = Window()

ui = @manipulate for ve = mdp.velocity_grid, vo = mdp.velocity_grid
    plots = []
    for (i, policy) in enumerate([policy_wide, policy_cv])
        xs = mdp.ego_grid[1:end]
        ys = mdp.human_grid
        ae = mdp.acceleration_grid[3]
        z = map(((y,x),) -> value(policy, Float64[x, ve, ae, y, vo]), Iterators.product(ys, xs))
        title = i == 1 ? "Assume random" : "Assume CV"
        xlabel = i == 1 ? "" : "Ego longitudinal position"
        za = map(((y,x),) -> action(policy, Float64[x, ve, ae, y, vo]), Iterators.product(ys, xs))
        push!(plots, heatmap(xs, ys, z, aspect_ratio=1, xlabel=xlabel, ylabel="Other longitudinal position", title=title, color=:plasma, clim=(0,1)))
        push!(plots, heatmap(xs, ys, za, aspect_ratio=1, title=title, color=:viridis, clim=(1,7)))
    end
    plot(plots..., layout=(2,2))
end
body!(w, ui)

ci = findfirst(isequal(0.0), res[:reward])
rng = MersenneTwister(ci)
frames = Frames(MIME("image/png"), fps=4)
AutoViz.render([mdp.env.roadway, s0], cam = StaticCamera(VecE2(0.0,0.0), 15.0))
for step in stepthrough(mdp, policy, s0, "s,a,r,sp", rng=rng, max_steps=100)
    s, a, r, sp = step
    @show s
    @show a
    @show r
    info = @printf("EGO VEHICLE sp = %3.0f lane = %s reward = %1.1f", s[1], mdp.ego_lane.tag, r)
    push!(frames, AutoViz.render([mdp.env.roadway, (sp, a)], cam=StaticCamera(VecE2(0.0,0.0), 15.0)))
    print(" . \n" )
end
write("out.gif", frames)

hist = simulate(hr, mdp, policy, s0)

AutoViz.render([mdp.env.roadway, state_hist(hist)[end]], cam = SceneFollowCamera())

undiscounted_reward(hist)

sf = SVector(900.0, 10.0, 0.0, 706.0, 16.0)

isterminal(mdp, sf)

hist.reward_hist[end]

isterminal(mdp, state_hist(hist)[end])

frames = Frames(MIME("image/png"), fps=8)
# AutoViz.render([mdp.env.roadway, s0], cam = SceneFollowCamera())
for step in stepthrough(mdp, policy, s0, "s,a,r,sp", rng=rng, max_steps=100)
    s, a, r, sp = step
    @show s
    @show a
    info = @printf("EGO VEHICLE sp = %3.0f lane = %s reward = %1.1f", s[1], mdp.ego_lane.tag, r)
    push!(frames, AutoViz.render([mdp.env.roadway, sp], cam=SceneFollowCamera()))
    print(" . \n" )
end
write("out.gif", frames)

d = transition(mdp, s0, 2)
sp = rand(rng, d)

# DANGER
rm("transition.jld2")
rm("reward.jld2")

# rng = MersenneTwister(1)

# policy = RandomPolicy(mdp, rng=rng)

# s0 = SVector(700.0, 10.0, 0.0, 650.0, 16.0)
# a = 2
# @time d = transition(mdp, s0, 2);
# sp = rand(rng, d)
# @time r = reward(mdp, sp, a, sp)

# @code_warntype reward(mdp, sp, a, sp)

# hr = HistoryRecorder(max_steps=100, rng=rng)
# hist = simulate(hr, mdp, policy, s0)

# frames = Frames(MIME("image/png"), fps=4)
# for step in stepthrough(mdp, policy, s0, "s,a,r,sp", rng=rng, max_steps=100)
#     s, a, r, sp = step
#     @show s
#     @show a
#     info = @printf("EGO VEHICLE s = %3.0f lane = %s", s[1], mdp.ego_lane.tag)
#     push!(frames, AutoViz.render([mdp.env.roadway, s], cam=SceneFollowCamera()))
#     print(" . \n" )
# end
# write("out.gif", frames)

# AutoViz.render([mdp.env.roadway, state_hist(hist)[2]], cam=SceneFollowCamera())
# Reel.set_output_type("gif")
#
# s = initialstate(mdp, rng)

## Try solving the MDP

# @show n_states(mdp)

# function testindexing(mdp::MergingMDP)
#     ss = states(mdp)
#     @showprogress for (i, s) in enumerate(ss)
#         si = stateindex(mdp, s)
#         @assert si == i "si = $si , i = $i, s = $s"
#         for a in actions(mdp)
#             d = transition(mdp, s, a)
#             for (sp, w) in weighted_iterator(d)
#                 try
#                     spi = stateindex(mdp, sp)
#                 catch
#                     throw("Error s = $s, a = $a, sp=$sp")
#                 end
#             end
#         end
#     end
# end

# function DiscreteValueIteration.transition_matrix_a_s_sp(mdp::MDP)
#     if isfile("transition.jld2")
#         @load "transition.jld2" trans_mat
#         @info "Loading transition from transition.jld2"
#         return trans_mat
#     else
#         # Thanks to zach
#         na = n_actions(mdp)
#         ns = n_states(mdp)
#         transmat_row_A = [Float64[] for _ in 1:n_actions(mdp)]
#         transmat_col_A = [Float64[] for _ in 1:n_actions(mdp)]
#         transmat_data_A = [Float64[] for _ in 1:n_actions(mdp)]

#         @showprogress for s in states(mdp)
#             si = stateindex(mdp, s)
#             for a in actions(mdp, s)
#                 ai = actionindex(mdp, a)
#                 if !isterminal(mdp, s) # if terminal, the transition probabilities are all just zero
#                     td = transition(mdp, s, a)
#                     for (sp, p) in weighted_iterator(td)
#                         if p > 0.0
#                             spi = stateindex(mdp, sp)
#                             push!(transmat_row_A[ai], si)
#                             push!(transmat_col_A[ai], spi)
#                             push!(transmat_data_A[ai], p)
#                         end
#                     end
#                 end
#             end
#         end
#         transmats_A_S_S2 = [sparse(transmat_row_A[a], transmat_col_A[a], transmat_data_A[a], n_states(mdp), n_states(mdp)) for a in 1:n_actions(mdp)]
#         # Note: not valid for terminal states
#         # @assert all(all(sum(transmats_A_S_S2[a], dims=2) .≈ ones(n_states(mdp))) for a in 1:n_actions(mdp)) "Transition probabilities must sum to 1"
#         return transmats_A_S_S2
#     end
# end

# function DiscreteValueIteration.reward_s_a(mdp::MDP)
#     if isfile("reward.jld2")
#         @load "reward.jld2" reward_mat
#         @info "Loading reward from reward.jld2"
#         return reward_mat
#     else
#         reward_S_A = fill(-Inf, (n_states(mdp), n_actions(mdp))) # set reward for all actions to -Inf unless they are in actions(mdp, s)
#         @showprogress for s in states(mdp)
#             for a in actions(mdp, s)
#                 td = transition(mdp, s, a)
#                 r = 0.0
#                 for (sp, p) in weighted_iterator(td)
#                     if p > 0.0
#                         r += p*reward(mdp, s, a, sp)
#                     end
#                 end
#                 reward_S_A[stateindex(mdp, s), actionindex(mdp, a)] = r
#             end
#         end
#         return reward_S_A
#     end
# end

# trans_mat = DiscreteValueIteration.transition_matrix_a_s_sp(mdp)
# reward_mat = DiscreteValueIteration.reward_s_a(mdp)

# @save "transition.jld2" trans_mat
# @save "reward.jld2" reward_mat

# solver = SparseValueIterationSolver(max_iterations=100, verbose=true)

# policy = solve(solver, mdp)

# ## Visualize policy
# rng = MersenneTwister(1)
# hr = HistoryRecorder(max_steps=100, rng=rng)
# s0 = SVector(700.0, 10.0, 0.0, 706.0, 16.0)
# hist = simulate(hr, mdp, policy, s0)

# frames = Frames(MIME("image/png"), fps=8)
# AutoViz.render([mdp.env.roadway, s0], cam = SceneFollowCamera())

# for step in stepthrough(mdp, policy, s0, "s,a,r,sp", rng=rng, max_steps=100)
#     s, a, r, sp = step
#     @show s
#     @show a
#     info = @printf("EGO VEHICLE sp = %3.0f lane = %s reward = %1.1f", s[1], mdp.ego_lane.tag, r)
#     push!(frames, AutoViz.render([mdp.env.roadway, sp], cam=SceneFollowCamera()))
#     print(" . \n" )
# end
# write("out.gif", frames)

# AutoViz.render([mdp.env.roadway, state_hist(hist)[end-1]], cam=SceneFollowCamera())

# @show actionvalues(policy, state_hist(hist)[end-1])

# ## Evaluate policy

# function evaluate(mdp::MergingMDP, policy::Policy; n_ep::Int=100, max_steps::Int=100, rng=rng)
#     avg_r = 0.0
#     sim = RolloutSimulator(rng=rng, max_steps=max_steps)
#     r = broadcast(x->simulate(sim, mdp, policy))

# policy.util

# state_space = collect(states(mdp));

# s = rand(rng, state_space)
# AutoViz.render([mdp.env.roadway, s], cam = SceneFollowCamera())

# collisioncheck(mdp, s)

# ## Rendering

# ego_posF = Frenet(roadway[LaneTag(MERGE_LANE_ID, 1)], 100.0)
# ego = Vehicle(VehicleState(ego_posF, roadway, 10.0), CARDEF, EGO_ID)
# scene = Scene()
# push!(scene, ego)

# roadway = mdp.env.roadway
# frames = Frames(MIME("image/png"), fps=4)
# for s = 1:10:1000.0
#     ego_posF = Frenet(roadway[LaneTag(MERGE_LANE_ID, 1)], s)
#     ego = Vehicle(VehicleState(ego_posF, roadway, 0.), CARDEF, EGO_ID)
#     scene = Scene()
#     push!(scene, ego)
#     info = @sprintf("EGO VEHICLE s = %3.0f lane = %s", ego.state.posF.s, ego.state.posF.roadind.tag)
#     to = TextOverlay(text=[info], pos=ego.state.posG, font_size = 20, incameraframe=true)
#     c = AutoViz.render(scene, roadway, [to])
#     push!(frames, c)
# end
# write("out.mp4", frames)

# ## Solving the MDP

# n_states(mdp)


text:
module LinearSolve

using ArrayInterface
using RecursiveFactorization
using Base: cache_dependencies, Bool
import Base: eltype, adjoint, inv
using LinearAlgebra
using IterativeSolvers:Identity
using SparseArrays
using SciMLBase: AbstractDiffEqOperator, AbstractLinearAlgorithm
using Setfield
using UnPack
using Requires
using SuiteSparse
using KLU
using DocStringExtensions

# wrap
import Krylov
import KrylovKit
import IterativeSolvers

using Reexport
@reexport using SciMLBase

abstract type SciMLLinearSolveAlgorithm <: SciMLBase.AbstractLinearAlgorithm end
abstract type AbstractFactorization <: SciMLLinearSolveAlgorithm end
abstract type AbstractKrylovSubspaceMethod <: SciMLLinearSolveAlgorithm end
abstract type AbstractSolveFunction <: SciMLLinearSolveAlgorithm end

# Traits

needs_concrete_A(alg::AbstractFactorization) = true
needs_concrete_A(alg::AbstractKrylovSubspaceMethod) = false
needs_concrete_A(alg::AbstractSolveFunction) = false

# Code

include("common.jl")
include("factorization.jl")
include("simplelu.jl")
include("iterative_wrappers.jl")
include("preconditioners.jl")
include("solve_function.jl")
include("default.jl")
include("init.jl")

const IS_OPENBLAS = Ref(true)
isopenblas() = IS_OPENBLAS[]

export LUFactorization, SVDFactorization, QRFactorization, GenericFactorization,
       GenericLUFactorization, SimpleLUFactorization, RFLUFactorization,
       UMFPACKFactorization, KLUFactorization

export LinearSolveFunction

export KrylovJL, KrylovJL_CG, KrylovJL_GMRES, KrylovJL_BICGSTAB, KrylovJL_MINRES,
       IterativeSolversJL, IterativeSolversJL_CG, IterativeSolversJL_GMRES,
       IterativeSolversJL_BICGSTAB, IterativeSolversJL_MINRES,
       KrylovKitJL, KrylovKitJL_CG, KrylovKitJL_GMRES

end


text:
function get_global_config(configtag, dfl = "")
    out = _run("git config $(configtag)"; verbose = false, ignorestatus = true)
    out = string(strip(out))
    isempty(out) ? dfl : out
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

# ╔═╡ e46441c4-97bd-11eb-330c-97bd5ac41f9e
html"""

Section 2.6

 Random Walks II

Lecture Video

"""

# ╔═╡ 85b45a43-d7bf-4597-a1a6-329b41dce20d
begin
    import Pkg

    Pkg.activate(mktempdir())
    Pkg.add([
        Pkg.PackageSpec(name="PlutoUI", version="0.7"),
        Pkg.PackageSpec(name="Plots", version="1"),
    ])

    using PlutoUI, Plots, LinearAlgebra, SparseArrays
end

# ╔═╡ 85c26eb4-c258-4a8b-9415-7b5f7ddff02a
TableOfContents(aside=true)

# ╔═╡ 2d48bfc6-9617-11eb-3a85-279bebd21332
md"""
# Julia concepts

- Customised display of objects
- Structured matrices in Julia
- `cumsum`
- Vectors of vectors
- Concatenation of vectors: `hcat`

- `heatmap` (Plots.jl)
- `surface` (Plots.jl)

"""

# ╔═╡ 30162386-039f-4cd7-9121-a3382be3c294
md"""
# Pascal's triangle

"""

# ╔═╡ 4e7b163e-dfd0-457e-b1f3-8807a4d8060a
md"""
Let's start by thinking about Pascal's triangle. (Note that [Pascal was not the first person to study these numbers](https://en.wikipedia.org/wiki/Pascal%27s_triangle).)
"""

# ╔═╡ e8ceab7b-45db-4393-bb8e-e000ecf78d2c
pascal(N) = [binomial(n, k) for n = 0:N, k=0:N]

# ╔═╡ 2d4dffb9-39e4-48de-9688-980b96814c9f
pascal(10)

# ╔═╡ 8ff66523-bc2e-4c53-975b-8ba4f99eb1c6
md"""
The non-zero entries are the **binomial coefficients**: the $k$th entry on the $n$th row is the coefficient of $x^k$ in the expansion $(1 + x)^n$, starting from $n=0$ in the first row and $k=0$ in the first column.
"""

# ╔═╡ 2868dd57-7164-4162-8c5d-30628dedeb7a
md"""
Note that there are 0s above the main diagonal -- in other words the matrix is **lower triangular**.

Julia has special types to represent some classes of structured matrices in the standard library `LinearAlgebra` package:
"""

# ╔═╡ f6009473-d3c1-444f-88ae-814f770e811b
L = LowerTriangular(pascal(10))

# ╔═╡ 9a368602-acd3-43fb-9dff-e407a4bab930
md"""
We see that the display is special: the known "**structural**" zeros are shown as dots, instead of numbers.
"""

# ╔═╡ 67517333-175f-48c4-a915-76658cbf1304
md"""
As we have already seen, Julia also has a sparse matrix type that we could use for this, in the `SparseArrays` standard library:
"""

# ╔═╡ d6832372-d336-4a54-bbcf-d0bb70e4de64
sparse(pascal(10))

# ╔═╡ 35f14826-f1e4-4977-a31a-0f6148fe25ad
md"""
For fun, let's look at where the entries of Pascal's triangle are odd. We'll make a slider:
"""

# ╔═╡ 7468fc5d-7f35-45e2-b5fc-7e63b562bc8f
@bind n_pascal Slider(1:63, show_value=true, default=1)

# ╔═╡ 1ca8aa3b-b05d-40f6-a925-2f0248b79ca2
sparse(isodd.(pascal(n_pascal)))

# ╔═╡ 38d88b7c-3b4f-430b-8d3c-f672ab0c7a49
md"""
Note that the visual representation of a sparse matrix changes at a particular value of $n$. For larger values of $n$ the **sparsity pattern** is displayed using dots for non-zero values.
"""

# ╔═╡ f4c9b02b-738b-4de1-9e9d-05b1616bee0b
md"""
The pattern is quite striking!
"""

# ╔═╡ d1c47afa-ab7f-4543-a161-e3ceb6f11eb4
md"""
You may be interested to know that there is an alternative way to look at Pascal's triangle:
"""

# ╔═╡ bf78e00f-05d9-4a05-8512-4924ef9e25f7
[binomial(i + j, i) for i = 0:10, j = 0:10]

# ╔═╡ b948830f-ead1-4f36-a237-c998f2f7deb8
md"""
and that in fact this can be produced from the previous version using matrix multiplication!:
"""

# ╔═╡ 15223c51-8d31-4a50-a8ff-1cb7d35de454
pascal(10) * pascal(10)'

# ╔═╡ 0bd45c4a-3286-427a-a927-15869be2ebfe
md"""
## Convolutions build Pascal's triangle!
"""

# ╔═╡ 999fb608-fb1a-46cb-82ca-f3f31fe617e1
pascal(6)

# ╔═╡ 6509e69a-6e50-4816-a98f-67ba437383fb
md"""
Where do those binomial coefficients come from? That's not how we learnt to build Pascal's triangle at school!

We learned to build each number up by *summing two adjacent numbers from the previous row*. In the lower-triangular representation we we sum the number immediately above with the one to the left of it.
"""

# ╔═╡ e58976a9-1784-441e-bb76-3011538b8ad0
md"""
We can think of this as... a **convolution** with the vector $[1, 1]$ !
Recall that convolutions like this are used as **image filters**, where we want to apply the same *local* operation everywhere along a vector (i.e. one in which we modify a vector in a neighbourhood).
"""

# ╔═╡ 1efc2b68-9313-424f-9850-eb4496cc8486
md"""
# Random walks: Independent and identically-distributed random variables
"""

# ╔═╡ 6e93ffda-217b-4d46-86b5-534ddc1bae90
md"""
The discussion of Pascal's triangle and convolutions will surprisingly turn up below. For now let's go back to thinking about **random walks**.

Recall that in a random walk, at each tick of a clock we take a step in space, say left or right.

Each spatial step is random, so we can think of each step as being a **random variable** with a certain probability distribution, for example the random variable that takes the value $+1$ with probability $\frac{1}{2}$ and $-1$ also with probability $\frac{1}{2}$.

Often we will think about random walks in which each step is "the same". What do we mean by "the same"? We don't mean that the *outcome* is the same, but rather that each step has the same *probability distribution*, i.e. the same probability of taking each possible value. In other words, each step is given by *copies of the same random variable*.

Steps are also *independent* of each other (the choice of step direction right now *does not affect* the choice at the next step). Hence the steps form a collection of **independent and identically-distributed random variables**, or IID random variables.
"""

# ╔═╡ 396d2490-3cb9-4f68-8fdf-9209d2010e02
md"""
## Random walks as a cumulative sum
"""

# ╔═╡ dc1c22e8-1c7b-43b7-8421-c2ca708931a5
md"""
Let's call the random variable that describes the $i$th *step* $X_i$. Then the overall position $S_n$ of the random walk at time $n$ is given by

$$S_n = X_1 + \cdots + X_n = \sum_{i=1}^n X_i.$$

(Here we have taken the initial position $S_0$ equal to 0. If the initial position is not zero, then this instead gives the **displacement** $S_n - S_0$ at time $n$.)

Recall that by the sum of two random variables we mean a new random variable whose outcomes are sums of the outcomes, with probabilities given by taking the possible pairs of outcomes giving each new outcome.
"""

# ╔═╡ e7a52b56-322c-4478-a670-dec1013c9bd8
md"""
## Cumulative sum
"""

# ╔═╡ 5f136388-5573-4814-a088-a66278acdbbe
md"""
We previously looked at sums like this. The difference now is that we are interested in the entire **trajectory**, i.e. the sequence $S_1$, $S_2$, $\ldots$, $S_n$, $\ldots$. The trajectory is given by **partial sums**:
"""

# ╔═╡ e440dc3b-bafd-4e0c-9fe8-13fce8eea22d
md"""
$$\begin{align}
S_1 & = X_1 \\
S_2 &= X_1 + X_2 \\
S_3 &= X_1 + X_2 + X_3
\end{align}$$

etc.
"""

# ╔═╡ 89d3e90d-3685-473d-aea4-0b7c5b80d7f7
md"""
Note that $S_{n+1}$ is *not* independent of $S_n$ (or the other $S_m$). E.g. if $S_{100}$ happens to be large, then $S_{101}$ will be about as large.

Thinking about types, as we did a couple of lectures ago, we would need to define a *new* type to represent a random walk, since we cannot generate consecutive values as independent samples.
"""

# ╔═╡ 203eea14-1c68-4a9f-ab18-9a2e5f408a79
md"""
How could we calculate a trajectory of a walk? Suppose we generate steps $X_i$, e.g.
"""

# ╔═╡ b6d4b045-a39f-4236-ace2-9f631e853d1b
steps = rand( (-1, +1), 10 )

# ╔═╡ e61b56be-d334-4c8f-aa8e-887bb27c058c
md"""
The whole trajectory is given by the **cumulative sum** also called "prefix sum" or "scan". Julia has the function `cumsum` to calculate this:
"""

# ╔═╡ 782f0b9a-3793-4abb-826b-9e14d6eae690
cumsum([1, 2, 3, 4])

# ╔═╡ 1764b56a-f297-4a4e-a931-31aa987ec785
md"""
So the trajectory is given by
"""

# ╔═╡ 8082092b-b6bf-4619-8776-39fdd6c9b7c1
cumsum(steps)

# ╔═╡ 34c0b850-5e95-4eb9-8435-3aae8d124772
plot(cumsum(steps), m=:o, leg=false, size=(500, 300))

# ╔═╡ 4115f7cb-d45f-4cc2-86bf-316c074393f8
md"""
Note that in Julia this is just a convenience function; it will be no more performant than writing the `for` loop yourself (unlike in some other languages, where `for` loops are slow).
"""

# ╔═╡ 6bc2f20d-3d09-425b-a471-44090dc3161e
md"""
# Evolving probability distributions in time
"""

# ╔═╡ 98994cb9-45dc-48aa-b62d-2407f7184bee
md"""
So far we have looked at single trajectories of random walks. We can think of this as the equivalent of sampling using `rand`.

Suppose that we run many trajectories of a random walk. At a given time $t$ we can ask where all of the walks are by taking a histogram of the $S_t^{(k)}$, where the superscript $^{(k)}$ denotes the $k$th trajectory out of the collection.

Doing so gives us the *probability distribution* of the random variable $S_t$.
Let's call $p_i^t := \mathbb{P}(S_t = i)$ the probability of being at site $i$ at time $t$. Then the probability distribution at time $t$ is given by the collection of all the $p_i^t$; we can group these into a vector $\mathbf{p}^t$.

Now we can ask what happens at the *next* time step $t+1$. Taking histograms again gives the probabilities $p_j^{t+1}$ and the vector $\mathbf{p}^t$ of all of them.

But $\mathbf{p}^{t+1}$ and $\mathbf{p}^{t}$ are related in some way, since we go from $S_t$ to $S_{t+1}$ by just taking a single step.

Let's think about the case of a simple random walk in one dimension. To arrive at site $i$ at time $t+1$, we must have been in one of the *neighbouring* sites at time $t$ and jumped with probability $\frac{1}{2}$ in the direction of site $i$. Hence we have

$$p_i^{t+1} = \textstyle \frac{1}{2} (p_{i-1}^{t} + p_{i+1}^{t}).$$

This is sometimes called the **master equation** (even though that is rather a useless, non-descriptive name); it describes how the *probability distribution of the random walk evolves in time*.
"""

# ╔═╡ f11f8d7d-cd3f-4585-aab4-083b892c6d4c
md"""
## Implementing time evolution of probabilities

Let's write a function to **evolve** a probability vector to the next time step for a simple random walk:
"""

# ╔═╡ fb804fe2-58be-46c9-9200-ceb8863d052c
function evolve(p)
	p′ = similar(p)   # make a vector of the same length and type
	                  # to store the probability vector at the next time step

	for i in 2:length(p)-1   # iterate over the *bulk* of the system

		p′[i] = 0.5 * (p[i-1] + p[i+1])

	end

	# boundary conditions:
	p′[1] = 0
	p′[end] = 0

	return p′
end

# ╔═╡ 40e457b4-616c-4fab-9c8e-2e5063129597
md"""
Wait... Do you recognise this?

This is just a **convolution** again! The kernel is now $[\frac{1}{2}, 0, \frac{1}{2}]$. Apart from the extra $0$ and the $\frac{1}{2}$, this is the *same* as in Pascal's triangle... so probabilities in simple random walk behave like Pascal's triangle!
"""

# ╔═╡ 979c1fbd-c9f6-4e8b-a648-6a0210fc9e7f
md"""
Note that just as with images we have a problem at the edges of the system: We need to specify **boundary conditions** on the first and last cells. For now we have put 0s at the boundary, corresponding to probability *escaping* at the boundary: any probability that arrives at the boundary at a given time step does not stay in the system. We can think of the probability as analogous to a chemical moving through a system that leaves our system (and e.g. moves into the outside world) when it reaches an edge.
"""

# ╔═╡ 583e3a92-01e7-4b88-9be0-f1e3b3c95005
md"""
We also need to specify an *initial* condition $\mathbf{p}_0$. This tells us where our walker is at time $0$. Suppose that all walkers start at site $i=0$. We will place this in the middle of our vector. Then the *probability* of being at $0$ is 1 (certainty), while the probability at any other site is 0 (impossible):
"""

# ╔═╡ 0b26efab-4e93-4d53-9c4d-faea68d12174
function initial_condition(n)

	p₀ = zeros(n)
	p₀[n ÷ 2 + 1] = 1

	return p₀
end

# ╔═╡ b9ce5af1-84f7-4a2d-92c9-de2c5498a88d
md"""
Now let's evolve the probability vector in time by applying the `evolve` function repeatedly, starting from the initial probability distribution:
"""

# ╔═╡ b48e55b7-4b56-41aa-9796-674d04adf5df
function time_evolution(p0, N)
	ps = [p0]
	p = p0

	for i in 1:N
		p = evolve(p)
		push!(ps, copy(p))
	end

	return ps
end

# ╔═╡ 53a36c1a-0b8c-4099-8854-08d73c9f118e
md"""
Let's visualise this:
"""

# ╔═╡ 6b298184-32c6-412d-a900-b113d6bd3d53
begin
	grid_size = 101

	p0 = initial_condition(grid_size)
end

# ╔═╡ b84a7255-7b0a-4ba1-8c87-9f5d3fa32ef3
ps = time_evolution(p0, 100)

# ╔═╡ c430c4de-d9bf-44e1-aa40-6b823d718b04
md"""
Note that `ps` is a `Vector`, whose elements are each `Vector`s! I.e. we have a vector of vectors. This is often a convenient way to *build up* a matrix.
"""

# ╔═╡ 99717c6e-f713-49d5-8ee5-a08c4a464223
ps[2]

# ╔═╡ 242ea831-c421-4a76-b658-2a57fa924a4f
md"""
t = $(@bind tt Slider(1:length(ps), show_value=true, default=1))
"""

# ╔═╡ aeaef573-1e90-45f3-a7fe-31ec5e2808c4
bar(ps[tt], ylim=(0, 1), leg=false, size=(500, 300), alpha=0.5)

# ╔═╡ efe186da-3273-4767-aafc-fc8eae01dbd9
md"""
### Concatenating vectors into a matrix
"""

# ╔═╡ 61201091-b8b3-4776-9be9-4c23d5ba88ba
md"""
Now we want to visualise this in a different way, for which we must join up (**concatenate**) all the probability vectors into a matrix:
"""

# ╔═╡ 66c4aed3-a04b-4a09-b954-79e816d2a3f7
M = reduce(hcat, ps)'

# ╔═╡ b135f6be-5e82-4c72-af11-0eb0d4141dec
md"""
We can visualise the matrix with a **heat map**:
"""

# ╔═╡ e74e18e3-ad08-4a53-a803-cd53564dca65
heatmap(M, yflip=true)

# ╔═╡ ed02f00f-1bcd-43fa-a56c-7be9968614cc
md"""
We can also visualise this as a 3D surface:
"""

# ╔═╡ 8d453f89-4a4a-42d0-8a00-9b153a3f435e
plotly()

# ╔═╡ f7de29b5-2a51-45e4-a0a5-f7f602681303
surface(M)

# ╔═╡ 7e817bad-dc51-4c29-a4fc-f7a8bb3663ca
md"""
But this is not necessarily very clear, so let's draw it ourselves as stacked histograms:
"""

# ╔═╡ 403d607b-6171-431b-a058-0aad0909846f
gr()

# ╔═╡ 80e033c4-a862-443d-a198-932f5822a44e
ylabels = [-(grid_size÷2):grid_size÷2;]

# ╔═╡ c8c16c14-26b0-4f83-8135-4f862ed90686
begin
	plot(leg=false)

	for which in 1:15
		for i in 1:length(ps[which])
			plot!([which, which], [-grid_size÷2 + i, -grid_size÷2 + i], [0, ps[which][i] ], c=which, alpha=0.8, lw = 2)
		end
	end

	xlims!(1, 15)

	plot!()
end

# ╔═╡ 29612df6-203f-42bc-b53b-86af618d60ec
let
	color_list = [:red, RGB(0, 1, 0.1), :blue]
	xs = []
	ys = []
	zs = []
	cs = []
	cs2 = []

	for which in 1:15
		for i in 1:length(ps[which])
			push!(xs, which, which, NaN)
			push!(ys, ylabels[i], ylabels[i], NaN)
			push!(zs, 0, ps[which][i], NaN)
			# push!(zs, 0, 1, NaN)D
			push!(cs, color_list[mod1(which, 3)], color_list[mod1(which, 3)], color_list[mod1(which, 3)])
		end
		push!(cs2, which)
	end
				# plot(xs, ys, zs)

	plot(leg=false)

	plot!(1:15, [0; cumsum(sign.(randn(14)))], zeros(15), alpha=0.6, m=:o, lw=2, c=color_list)
	plot!(xs, ys, zs, c=cs, xlims=(1, 15), ylims=(-30, 30), zlims=(0, 1), lw=2.5, alpha=0.7, xticks=(1:15), yticks=(-20:10:20))

xlabel!("t")
	ylabel!("space")

end

# ╔═╡ Cell order:
# ╟─e46441c4-97bd-11eb-330c-97bd5ac41f9e
# ╠═85b45a43-d7bf-4597-a1a6-329b41dce20d
# ╠═85c26eb4-c258-4a8b-9415-7b5f7ddff02a
# ╟─2d48bfc6-9617-11eb-3a85-279bebd21332
# ╟─30162386-039f-4cd7-9121-a3382be3c294
# ╟─4e7b163e-dfd0-457e-b1f3-8807a4d8060a
# ╠═e8ceab7b-45db-4393-bb8e-e000ecf78d2c
# ╠═2d4dffb9-39e4-48de-9688-980b96814c9f
# ╟─8ff66523-bc2e-4c53-975b-8ba4f99eb1c6
# ╟─2868dd57-7164-4162-8c5d-30628dedeb7a
# ╠═f6009473-d3c1-444f-88ae-814f770e811b
# ╟─9a368602-acd3-43fb-9dff-e407a4bab930
# ╟─67517333-175f-48c4-a915-76658cbf1304
# ╠═d6832372-d336-4a54-bbcf-d0bb70e4de64
# ╟─35f14826-f1e4-4977-a31a-0f6148fe25ad
# ╠═7468fc5d-7f35-45e2-b5fc-7e63b562bc8f
# ╠═1ca8aa3b-b05d-40f6-a925-2f0248b79ca2
# ╟─38d88b7c-3b4f-430b-8d3c-f672ab0c7a49
# ╟─f4c9b02b-738b-4de1-9e9d-05b1616bee0b
# ╟─d1c47afa-ab7f-4543-a161-e3ceb6f11eb4
# ╠═bf78e00f-05d9-4a05-8512-4924ef9e25f7
# ╟─b948830f-ead1-4f36-a237-c998f2f7deb8
# ╠═15223c51-8d31-4a50-a8ff-1cb7d35de454
# ╟─0bd45c4a-3286-427a-a927-15869be2ebfe
# ╠═999fb608-fb1a-46cb-82ca-f3f31fe617e1
# ╟─6509e69a-6e50-4816-a98f-67ba437383fb
# ╟─e58976a9-1784-441e-bb76-3011538b8ad0
# ╟─1efc2b68-9313-424f-9850-eb4496cc8486
# ╟─6e93ffda-217b-4d46-86b5-534ddc1bae90
# ╟─396d2490-3cb9-4f68-8fdf-9209d2010e02
# ╟─dc1c22e8-1c7b-43b7-8421-c2ca708931a5
# ╟─e7a52b56-322c-4478-a670-dec1013c9bd8
# ╟─5f136388-5573-4814-a088-a66278acdbbe
# ╟─e440dc3b-bafd-4e0c-9fe8-13fce8eea22d
# ╟─89d3e90d-3685-473d-aea4-0b7c5b80d7f7
# ╟─203eea14-1c68-4a9f-ab18-9a2e5f408a79
# ╠═b6d4b045-a39f-4236-ace2-9f631e853d1b
# ╟─e61b56be-d334-4c8f-aa8e-887bb27c058c
# ╠═782f0b9a-3793-4abb-826b-9e14d6eae690
# ╟─1764b56a-f297-4a4e-a931-31aa987ec785
# ╠═8082092b-b6bf-4619-8776-39fdd6c9b7c1
# ╠═34c0b850-5e95-4eb9-8435-3aae8d124772
# ╟─4115f7cb-d45f-4cc2-86bf-316c074393f8
# ╟─6bc2f20d-3d09-425b-a471-44090dc3161e
# ╟─98994cb9-45dc-48aa-b62d-2407f7184bee
# ╟─f11f8d7d-cd3f-4585-aab4-083b892c6d4c
# ╠═fb804fe2-58be-46c9-9200-ceb8863d052c
# ╟─40e457b4-616c-4fab-9c8e-2e5063129597
# ╟─979c1fbd-c9f6-4e8b-a648-6a0210fc9e7f
# ╟─583e3a92-01e7-4b88-9be0-f1e3b3c95005
# ╠═0b26efab-4e93-4d53-9c4d-faea68d12174
# ╟─b9ce5af1-84f7-4a2d-92c9-de2c5498a88d
# ╠═b48e55b7-4b56-41aa-9796-674d04adf5df
# ╟─53a36c1a-0b8c-4099-8854-08d73c9f118e
# ╠═6b298184-32c6-412d-a900-b113d6bd3d53
# ╠═b84a7255-7b0a-4ba1-8c87-9f5d3fa32ef3
# ╟─c430c4de-d9bf-44e1-aa40-6b823d718b04
# ╠═99717c6e-f713-49d5-8ee5-a08c4a464223
# ╟─242ea831-c421-4a76-b658-2a57fa924a4f
# ╠═aeaef573-1e90-45f3-a7fe-31ec5e2808c4
# ╟─efe186da-3273-4767-aafc-fc8eae01dbd9
# ╟─61201091-b8b3-4776-9be9-4c23d5ba88ba
# ╠═66c4aed3-a04b-4a09-b954-79e816d2a3f7
# ╟─b135f6be-5e82-4c72-af11-0eb0d4141dec
# ╠═e74e18e3-ad08-4a53-a803-cd53564dca65
# ╟─ed02f00f-1bcd-43fa-a56c-7be9968614cc
# ╠═8d453f89-4a4a-42d0-8a00-9b153a3f435e
# ╠═f7de29b5-2a51-45e4-a0a5-f7f602681303
# ╟─7e817bad-dc51-4c29-a4fc-f7a8bb3663ca
# ╠═403d607b-6171-431b-a058-0aad0909846f
# ╠═80e033c4-a862-443d-a198-932f5822a44e
# ╟─c8c16c14-26b0-4f83-8135-4f862ed90686
# ╠═29612df6-203f-42bc-b53b-86af618d60ec


text:
# TODO: dealing with type assertions for points tuples vs arrays

# module for determining line of sight visibility in 2D
module RayCasters

export isVisible, inBuilding
export rayIntersectsSeg, segmentsIntersect
export nearest_wall

using Maps

# returns true if p1 and p2 are within LOS given the map geometry
function isVisible(map::Map, p1, p2)
    return isVisible(map, p1[1], p1[2], p2[1], p2[2])
end

function isVisible(map::Map, p1x, p1y, p2x, p2y)
    for b in map.buildings
        if !isVisible(b, p1x, p1y, p2x, p2y)
            return false
        end
    end
    return true
end

# check if points within line of sight given a polygon
# overloaded to add more shapes
function isVisible(b::Polygon, p1x, p1y, p2x, p2y)
    points = b.points
    nPoints = size(points,1)

    for i = 1:nPoints
        # last point makes a segment with the first
        if i == nPoints
            if segmentsIntersect(p1x, p1y, p2x, p2y,
                                 points[i,1], points[i,2], points[1,1], points[1,2])
                return false
            end
        else
            if segmentsIntersect(p1x, p1y, p2x, p2y,
                                 points[i,1], points[i,2], points[i+1,1], points[i+1,2])
                return false
            end
        end
    end
    return true
end

# checks if a point is inside a building
function inBuilding(map::Map, p)
    return inBuilding(map, p[1], p[2])
end

# checks if a point is inside a building
function inBuilding(map::Map, p1, p2)
    for b in map.buildings
        if inBuilding(b, p1, p2)
            return true
        end
    end
    return false
end

# check if a point p is inside a polygon
function inBuilding(b::Polygon, px, py)

    points = b.points
    nPoints = size(points,1)
    count = 0

    for i = 1:nPoints
        # last point makes a segment with the first
        if i == nPoints
            if rayIntersectsSeg(px, py, points[i,1], points[i,2], points[1,1], points[1,2])
                count += 1
            end
        else
            if rayIntersectsSeg(px, py, points[i,1], points[i,2], points[i+1,1], points[i+1,2])
                count += 1
            end
        end
    end

    if count % 2 == 1
        return true
    end
    return false
end

# checks if a point p intersect a line segment with ends a, b
# adopted from python implementation on Rosetta Code
# source: http://rosettacode.org/wiki/Ray-casting_algorithm#Python
function rayIntersectsSeg(px, py, ax, ay, bx, by)
    EPSILON = 0.001
    _huge = 1e300
    _tiny = 1e-300

    if ay > by
        ax,ay,bx,by = bx,by,ax,ay
    end
    if py == ay || py == by
        py = py + EPSILON
    end

    intersect = false

    if (py > by || py < ay) || (px > max(ax, bx))
        return false
    end

    if px < min(ax, bx)
        intersect = true
    else
        if abs(ax - bx) > _tiny
            m_red = (by - ay) / float(bx - ax)
        else
            m_red = _huge
        end
        if abs(ax - px) > _tiny
            m_blue = (py - ay) / float(px - ax)
        else
            m_blue = _huge
        end
        intersect = m_blue >= m_red
    end
    return intersect
end

# check if two 2D line segments intersect
# based on javascript code from:
# https://github.com/pgkelley4/line-segments-intersect/blob/master/js/line-segments-intersect.js
function segmentsIntersect(p::Array, p2::Array, q::Array, q2::Array)
    r = p2 - p
    s = q2 - q
    uNum = cross2D(q - p, r)
    denom = cross2D(r, s)
    # check if lines are collinear
    if uNum == 0 && denom == 0
        # do they touch?
        if p == q || p == p2 || p2 == q || p2 == q2
            return true
        end
        # do they overlap?
        return ((q[1] - p[1] < 0) != (q[1] - p2[1] < 0) != (q2[1] - p[1] < 0) != (q2[1] - p2[1] < 0) ||
                (q[2] - p[2] < 0) != (q[2] - p2[2] < 0) != (q2[2] - p[2] < 0) != (q2[2] - p2[2] < 0))
    end
    # lines are parallel
    if denom == 0
        return false
    end
    u = uNum / denom
    t = cross2D(q - p, s) / denom
    return (t >= 0) && (t <= 1) && (u >= 0) && (u <= 1)
end

function segmentsIntersect(px::Float64, py::Float64, p2x::Float64, p2y::Float64,
                           qx::Float64, qy::Float64, q2x::Float64, q2y::Float64)
    p  = [px, py]
    p2 = [p2x, p2y]
    q  = [qx, qy]
    q2 = [q2x, q2y]
    return segmentsIntersect(p, p2, q, q2)
end

# compute cross product in 2D
function cross2D(l1, l2)
    return l1[1] * l2[2] - l1[2] * l2[1]
end

# computes distance to nearest wall for a given direction
function nearest_wall(map::Map, p::Vector{Float64}, dir::Symbol)

    # find y-intercept
    d = Inf
    for b in map.buildings
        nd = nearest_wall(b, p, dir)
        if nd < d
            d = nd
        end
    end

    return d
end

function nearest_wall(b::Polygon, p::Vector{Float64}, dir::Symbol)
    x = p[1]
    y = p[2]
    yo = p[2]

    # only works for down
    step = y / 100.
    for i = 1:100
        y -= i*step
        if inBuilding(b, x, y)
            return abs(yo - y)
        end
    end
    return yo
end

end # module

max_stars_repo_name: etotheipluspi/ParticleFilters

max_issues_repo_name: etotheipluspi/ParticleFilters

max_forks_repo_name: etotheipluspi/ParticleFilters


text:
#Code written by Skipper af Rosenborg and Esben Bager
using StatsBase
using DataFrames
using CSV
include("SupportFunction.jl")
println("Leeeeroooy Jenkins")

#Initilise parameters
industry = "NoDur"
folder = "240"
industryArr = ["NoDur", "Durbl", "Manuf", "Enrgy", "HiTec", "Telcm", "Shops", "Hlth", "Utils", "Other"]
path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Results/IndexData/LassoTest/"*industry*"/"*folder*"-1/"
for industry = industryArr
    for i = [12, 24, 36, 48, 120, 240]
        println("Industry ",industry," time period ",i)
        path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Results/IndexData/LassoTest/"*industry*"/"*folder*"-1/"
        folder = string(i)
        #Write a summary file for given industry and training period
        writeFile()
    end
end

path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Results/IndexData/LassoTest/"*industry*"/"*folder*"-1/"

stringArr = Array{String}(10)
for j = 1:10
    stringArr[j] = "gamma =$j"
end

function logStuff(VIX, raw, timeTrans, expTrans, TA, Top25)
    timeSpan = parse(Int64,folder)
    tempString = ""
    #Define datasat based on input
    if raw == 0
        if VIX ==1
            tempString = tempString*"VIX_Macro"
        else
            tempString = tempString*"Macro"
        end
    else
        tempString = tempString*"Raw"
    end

    if timeTrans == 1
        tempString = tempString*"Time"
    end

    if expTrans == 1
        tempString = tempString*"Exp"
    end

    if TA == 1
        tempString = tempString*"TA"
    end

    #Load real values
    path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Results/IndexData/LassoTest/"*industry*"/"
    cd(path)
    y_realFull = CSV.read("RealValue.CSV",
        delim = ',', nullable=false, header=["Date", "y_real"], types = [Int64, Float64])

    y_realTimeSpan = zeros(size(y_realFull)[1]-timeSpan)
    for k = 1:size(y_realFull)[1]-timeSpan-1
        y_realTimeSpan[k] = mean(y_realFull[1:k+timeSpan-1,2])
    end

    path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Results/IndexData/LassoTest/"*industry*"/"*folder*"-1/"
    cd(path)

    y_real = CSV.read(folder*"_"*tempString*"_real.CSV", delim = ',', nullable = false, header = ["Iteration",
        "Date", "Reseccion", "Real Value"], types =[Float64, Int64, Int64, Float64], datarow=2)

    #If load a dataset based on top250
    if Top25 == 1 &&timeSpan == 240
        y_hat = CSV.read(folder*"_"*tempString*"_over25PercentPredictionTermspredicted.CSV",
            delim = ',', nullable=false, header=vcat("Iteration","Date","Recession",stringArr), types=[Float64, Int64, Int64, Float64,
            Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64], datarow=2)
        tempString = tempString*"_over25PercentPredictionTerms"
    else
        y_hat = CSV.read(folder*"_"*tempString*"_predicted.CSV",
            delim = ',', nullable=false, header=vcat("Iteration","Date","Recession",stringArr), types=[Float64, Int64, Int64, Float64,
            Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64], datarow=2)
    end

    #Make sure we start at same place
    startIndex=find(x -> x == 194608, y_real[:,2])[1]

    #Trim datasets
    y_realTimeSpan = y_realTimeSpan[startIndex:end-6,end]
    y_real = y_real[startIndex:end-5,4:end]
    y_hat = y_hat[startIndex:end-5,:4:end]

    nRows = size(y_hat)[1]
    nCols = size(y_hat)[2]

    #Initilise arrays
    classificationRate = Array{Float64}(nCols)
    meanErr = Array{Float64}(nCols)
    RMSE = Array{Float64}(nCols)
    Rsquare = Array{Float64}(nCols)

    #Calcualte CR, AME, RMSE, R^2 for each lambda value
    for g = 1:10
        currentCol = g
        # Calculate classifaction rate
        trueClass = 0
        for row = 1:nRows
            if sign(y_real[row,1]) == sign(y_hat[row,currentCol])
                trueClass += 1
            end
        end
        classificationRate[currentCol] = trueClass/nRows

        # Mean Error rate
        # Calculate RMSE
        # Calculate R^2
        indvRSquared = Array{Float64}(nRows,1)
        errSum = 0
        errSquaredSum = 0
        SSres = 0
        SSTO = 0
        for row = 1:nRows
            SSres += (y_real[row,1] - y_hat[row,currentCol])^2
            SSTO += (y_real[row,1]-y_realTimeSpan[row])^2
            errSum += abs(y_real[row,1] - y_hat[row,currentCol])
            errSquaredSum += (y_real[row,1] - y_hat[row,currentCol])^2
        end
        meanErr[currentCol] = errSum/nRows
        RMSE[currentCol] = sqrt(errSquaredSum/nRows)
        Rsquare[currentCol] = 1-SSres/SSTO
    end

    #Output to logarray
    tempArr = zeros(1,8)

    classRate = findmax(classificationRate)[1]
    classIndex = findmax(classificationRate)[2]
    OOSR = findmax(Rsquare)[1]
    OOSRIndex = findmax(Rsquare)[2]
    meanErrVal = findmin(meanErr)[1]
    meanErrIndex = findmin(meanErr)[2]
    RMSEVal = findmin(RMSE)[1]
    RMSEIndex = findmin(RMSE)[2]
    tempArr[1,1] = findmax(classificationRate)[1]
    tempArr[1,2] = findmax(Rsquare)[1]
    tempArr[1,3] = findmin(meanErr)[1]
    tempArr[1,4] = findmin(RMSE)[1]
    tempArr[1,5] = findmax(classificationRate)[2]
    tempArr[1,6] = findmax(Rsquare)[2]
    tempArr[1,7] = findmin(meanErr)[2]
    tempArr[1,8] = findmin(RMSE)[2]
    #println("Classification: ",findmax(classificationRate))
    #println("OoS R²: ",findmax(Rsquare))
    #println("Mean Absolute Error: ",findmin(meanErr))
    #println("RMSE: ",findmin(RMSE))

    return tempArr, tempString
end

function writeFile()
    counterIter = 1
    logArr = zeros(25,8)
    stringToAppend = Array{String}(25,1)

    #Loop over datasets without macro prediction terms
    for timeTrans = 0:1
    	for expTrans = 0:1
    		for TA = 0:1
    			logArr[counterIter,:], stringToAppend[counterIter,1] = logStuff(0, 1, timeTrans, expTrans, TA, 0)
                counterIter += 1
    		end
    	end
    end

    #Loop over datasets with macro prediction terms
    for VIX = 0:1
    	for timeTrans = 0:1
    		for expTrans = 0:1
    			for TA = 0:1
    				logArr[counterIter,:], stringToAppend[counterIter,1]= logStuff(VIX, 0, timeTrans, expTrans, TA, 0)
                    counterIter += 1
    			end
    		end
    	end
    end

    logArr[counterIter,:], stringToAppend[counterIter, 1] = logStuff(1, 0, 1, 0, 1, 1)

    #Output data
    path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Results/IndexData/LassoTest/"*industry*"/Summary New"*folder
    f = open(path*".csv", "w")
    write(f, "Dataset, Classification Rate, R^2, Mean Error, RMSE, Classification Rate Index, R^2 Index, Mean Error Index, RMSE Index\n")
    writecsv(f,hcat(stringToAppend,logArr))
    close(f)
end

max_stars_repo_path: ML/Code to hand in/DataAnalysis_Lasso.jl

max_stars_repo_name: skipperafrosenborg/Thesis

max_issues_repo_name: skipperafrosenborg/Thesis

max_forks_repo_name: skipperafrosenborg/Thesis


text:
using Random
using Plots
using ImageView
using TestImages
using ImageCore
using ProgressBars

function flatten(X)
    return reshape(X, (:, last(size(X))))
end

function euclidean_distance(A, B)
    return sum((A-B).^2)^0.5
end

function find_neighbours(p, neighbourhood, ϵ, distanceFunc)
    points = flatten(neighbourhood)
    neighbours::Vector{Int} = []
    for i in range(1, size(points)[1])
        q = points[i, :]
        distance = distanceFunc(p, q)
        if (0 < distance) & (distance <= ϵ)
            push!(neighbours, i)
        end
    end
    return neighbours
end

function dbscan(X, ϵ, minNeighbourhoodSize, distanceFunc)
    points = flatten(X)
    n_points = length(points[:,1])
    cluster::Int = 0
    neighbourhoods = Dict()
    labels = Dict()
    for i in tqdm(range(1, n_points))
        if ~haskey(labels, i)
            get!(neighbourhoods, i) do
                find_neighbours(points[i,:], points, ϵ, distanceFunc)
            end
            if length(neighbours_i) < minNeighbourhoodSize
                labels[i] = -1 # Noise
            else
                cluster+=1
                labels[i] = cluster

                local_neighbourhood = copy(neighbours_i)
                while ~isempty(local_neighbourhood)
                    j = pop!(local_neighbourhood)
                    if ~haskey(labels, j) || (labels[j] == -1)
                        labels[j] = cluster
                        neighbours_j = find_neighbours(points[j,:], points, ϵ, distanceFunc) # TODO store as a dict
                        if length(neighbours_j) >= minNeighbourhoodSize
                            append!(local_neighbourhood, neighbours_j)
                        end

                    end
                end
            end
        end
    end
    label_map = Array{Int}(undef, n_points, 1)
    for (k, v) in labels
        label_map[k] = v
    end
    label_map = reshape(label_map, size(X)[1:end-1])
    return label_map
end

img = testimage("lighthouse")

data = permutedims(channelview(img), [2,3,1]);
data = convert(Array{Float64, size(data)[end]}, data)

mask = dbscan(data, 0.1, 1000, euclidean_distance)

imshow(mask)

max_stars_repo_name: jkhebel/algorithms

max_issues_repo_name: jkhebel/algorithms

max_forks_repo_name: jkhebel/algorithms


text:
#####################

#TODO This code seems to assume logspace flows as floating point numbers. if so, enforca that on type F
function marginal_pass_up(circuit::UpFlowΔ{O,F}, data::XData{E}) where {E <: eltype(F)} where {O,F}
    resize_flows(circuit, num_examples(data))
    cache = zeros(Float64, num_examples(data)) #TODO: fix type later
    marginal_pass_up_node(n::UpFlowΔNode, ::PlainXData) = ()

    function marginal_pass_up_node(n::UpFlowLiteral{O,F}, cache::Array{Float64}, data::PlainXData{E}) where {E <: eltype(F)} where {O,F}
        pass_up_node(n, data)
        # now override missing values by 1
        npr = pr(n)
        npr[feature_matrix(data)[:,variable(n)] .< zero(eltype(F))] .= 1
        npr .= log.( npr .+ 1e-300 )
        return nothing
    end

    function marginal_pass_up_node(n::UpFlow⋀Cached, cache::Array{Float64}, ::PlainXData)
        pr(n) .= 0
        for i=1:length(n.children)
            # pr(n) .+= pr(n.children[i])
            broadcast!(+, pr(n), pr(n), pr(n.children[i]))
        end
        return nothing
    end

    function marginal_pass_up_node(n::UpFlow⋁Cached, cache::Array{Float64}, ::PlainXData)
        pr(n) .= 1e-300
        for i=1:length(n.children)
            cache .= 0
            # broadcast reduced memory allocation, though accessing prob_origin(n).log_thetas[i] still allocates lots of extra memory,
            # it is proabably due to derefrencing the pointer
            broadcast!(+, cache, pr(n.children[i]), prob_origin(n).log_thetas[i])
            broadcast!(exp, cache, cache)
            broadcast!(+, pr(n), pr(n), cache)
        end
        broadcast!(log, pr(n), pr(n));
        return nothing
    end

    ## Pass Up on every node in order
    for n in circuit
        marginal_pass_up_node(n, cache, data)
    end
    return nothing
end

##### marginal_pass_down

function marginal_pass_down(circuit::DownFlowΔ{O,F}) where {O,F}
    resize_flows(circuit, flow_length(origin(circuit)))
    for n in circuit
        reset_downflow_in_progress(n)
    end
    for downflow in downflow_sinks(circuit[end])
        # initialize root flows to 1
        downflow.downflow .= one(eltype(F))
    end
    for n in Iterators.reverse(circuit)
        marginal_pass_down_node(n)
    end
end

marginal_pass_down_node(n::DownFlowΔNode) = () # do nothing
marginal_pass_down_node(n::DownFlowLeaf) = ()

function marginal_pass_down_node(n::DownFlow⋀Cached)
    # todo(pashak) might need some changes, not tested, also to convert to logexpsum later
     # downflow(n) = EF_n(e), the EF for edges or leaves are note stored
    for c in n.children
        for sink in downflow_sinks(c)
            if !sink.in_progress
                sink.downflow .= downflow(n)
                sink.in_progress = true
            else
                sink.downflow .+= downflow(n)
            end
        end
    end
end

function marginal_pass_down_node(n::DownFlow⋁Cached)
    # todo(pashak) might need some changes, not tested, also to convert to logexpsum later
    # downflow(n) = EF_n(e), the EF for edges or leaves are note stored
    for (ind, c) in enumerate(n.children)
        for sink in downflow_sinks(c)
            if !sink.in_progress
                sink.downflow .= downflow(n) .* exp.(prob_origin(n).log_thetas[ind] .+ pr(origin(c)) .- pr(origin(n)) )
                sink.in_progress = true
            else
                sink.downflow .+= downflow(n) .* exp.(prob_origin(n).log_thetas[ind] .+ pr(origin(c)) .- pr(origin(n)))
            end
        end
    end
end

#### marginal_pass_up_down

function marginal_pass_up_down(circuit::DownFlowΔ{O,F}, data::XData{E}) where {E <: eltype(F)} where {O,F}
    @assert !(E isa Bool)
    marginal_pass_up(origin(circuit), data)
    marginal_pass_down(circuit)
end

"""
Class dominane with missing values

"""
function dominance_pass_up(circuit::UpFlowΔ{O,F}, data::XData{E}, y::Integer) where {E <: eltype(F)} where {O,F}
    resize_flows(circuit, num_examples(data))
    node_cache = Dict{UpFlowΔNode, Int64}()
    for (i, n) in enumerate(circuit)
        node_cache[n] = i
    end
    # pr() is  a vector of triples (a,b,c), where a:Int8 is the type of sub-circuit,b:Float64 is minval, c:Float is the maxval
    dominance_pass_up_node(n::UpFlowΔNode, ::PlainXData,  y::Integer) = ()
    function dominance_pass_up_node(n::UpFlowLiteral{O,F}, data::PlainXData{E}, y::Integer) where {E <: eltype(F)} where {O,F}
        nsize =  num_examples(data)
        npr = pr(n)
        val = zeros( length(npr))
        if positive(n)
            val = feature_matrix(data)[:,variable(n)]
        else
            val =  ones(Float64,length(npr)) - feature_matrix(data)[:,variable(n)] # using views here strangely seems to save a few allocations and a little time
        end
        type = zeros( length(npr))
        # Type 0 observed, 1 class, 2 missing
        #  override missing values by -1
        val[feature_matrix(data)[:,variable(n)] .< zeros(Float64,length(npr))] .= 1
        type[feature_matrix(data)[:,variable(n)] .< zeros(Float64,length(npr))] .= 2
        # override class Y values by 1 and -1
        if y in variable(n)
            if positive(n)
                val .= 1
                type .= 1
            else
                val .= -1
                type .= 1
            end
        end
        pr(n) .= [x for x in zip(val,val,type)]
        # println("i ", variable(n) , "  val ",pr(n))
        return nothing
    end

    function dominance_pass_up_node(n::UpFlow⋀Cached, ::PlainXData, y::Integer)
        # println(" *** ")
        npr = pr(n)
        val_min = ones( length(npr))
        val_max = ones( length(npr))
        type = zeros( length(npr))
        index_max = zeros(Bool,num_examples(data)) # neg class misture missing
        index_class = zeros(Bool,num_examples(data))
        index_obs = zeros(Bool,num_examples(data)) #Observed and class
        index_min = zeros(Bool,num_examples(data)) # missing and observables
        # print("-***- ", node_cache[n])
        for i=1:length(n.children)
            val_min_i = getindex.(pr(n.children[i]),1)
            val_max_i = getindex.(pr(n.children[i]),2)
            type_i = getindex.(pr(n.children[i]),3)
            #Verify if some sub-circuit is class and is negative
            index_max .= ((val_min_i .< 0.0) .& (type_i .== 1.0)) .| index_max
            index_class .= ((val_min_i .>= 0.0) .& (type_i .== 1.0)) .| index_class
            # println(" max ",val_min ," t ",type_i, index_max )
            index_min .= (.!(index_max .| index_class)) .& ((type_i .== 2.0).| index_min)
        end
        # println("indexmin", index_min)
        index_obs .= .!(index_max .| index_class .| index_min)
        type[index_min] .= 2.0
        type[index_max .| index_class] .= 1.0
        # type[index_obs] .= 0.0
        for i=1:length(n.children)
            val_min_i = getindex.(pr(n.children[i]),1)
            val_max_i = getindex.(pr(n.children[i]),2)
            # Combining negative class and Missing or observed
                #The child with the neg class
            val_min[(index_max).& (val_min_i .< 0.0)] .= val_min[(index_max).& (val_min_i .< 0.0)] .* val_min_i[(index_max).& (val_min_i .< 0.0)]
                #The child with the othres types of variables
            val_min[(index_max).& (val_min_i .>= 0.0)] .= val_min[(index_max).& (val_min_i .>= 0.0)] .* val_max_i[(index_max).& (val_min_i .>= 0.0)]
                # Valman=valmin= valmin(i)*valmax(j)
            val_max[index_max] .= val_min[index_max]

            # Combining possitive classes and/or  observed
            val_min[index_class.& index_obs] .= val_min[index_class.& index_obs] .* val_min_i[index_class.& index_obs]
            val_max[index_class.& index_obs] .= val_min[index_class.& index_obs]

            #Combining with missing variables
            val_min[index_min] .= val_min[index_min] .* val_min_i[index_min]
            val_max[index_min] .= val_max[index_min] .* val_max_i[index_min]

            #Combining only observable
            val_min[index_class] .= val_min[index_class] .* val_min_i[index_class]
            val_max[index_class] .= val_min[index_class]
            # print( "   arg ",i,"  ",pr(n.children[i]))
            # pr(n) .*= pr(n.children[i])
        end
        pr(n) .= [x for x in zip(val_min,val_max,type)]
        # println("  * ",pr(n))
        return nothing
    end

    function dominance_pass_up_node(n::UpFlow⋁Cached,  ::PlainXData, y::Integer)
        # println("  +++ ")
        # print("-+++- ", node_cache[n])
        npr = pr(n)
        prs_min = zeros( length(n.children), num_examples(data))
        prs_max = zeros( length(n.children), num_examples(data))
        index_min = zeros( length(npr)) #Missing values and observed type 2
        index_sum = zeros( length(npr)) # Observed and/or class values  type 0 or 1
        val_min = zeros( length(npr))
        val_max = zeros( length(npr))
        type = zeros( length(npr))
        for i=1:length(n.children)
            val_min_i = getindex.(pr(n.children[i]),1)
            val_max_i = getindex.(pr(n.children[i]),2)
            type_i = getindex.(pr(n.children[i]),3)
            index_min = (type_i .== 2.0)
            index_sum = (type_i .== 1.0) .| (type_i .== 0.0)
            type .= type_i
            weight = exp.(prob_origin(n).log_thetas[i])
            prs_min[i,:] .= val_min_i .* weight
            prs_max[i,:] .= val_max_i .* weight
            # print( "   arg ",i,"  ",pr(n.children[i]))
        end
        val_min[index_min] .= minimum(prs_min[:,index_min], dims = 1)[1,:]
        val_max[index_min] .= maximum(prs_max[:,index_min], dims = 1)[1,:]

        val_min[index_sum] .= sum(prs_min[:,index_sum], dims = 1)[1,:]
        val_max[index_sum] .= sum(prs_max[:,index_sum], dims = 1)[1,:]
        pr(n) .= [x for x in zip(val_min,val_max,type)]
        # print("-- ", node_cache[n])
        # println("  +    ",pr(n))
        return nothing
    end

    for n in circuit
        dominance_pass_up_node(n,data, y)
    end
    return nothing
end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: iot_1click_projects

using AWS.Compat
using AWS.UUIDs
"""
    AssociateDeviceWithPlacement()

Associates a physical device with a placement.

# Required Parameters
- `deviceId`: The ID of the physical device to be associated with the given placement in the project. Note that a mandatory 4 character prefix is required for all deviceId values.
- `deviceTemplateName`: The device template name to associate with the device ID.
- `placementName`: The name of the placement in which to associate the device.
- `projectName`: The name of the project containing the placement in which to associate the device.

"""

associate_device_with_placement(deviceId, deviceTemplateName, placementName, projectName; aws_config::AWSConfig=global_aws_config()) = iot_1click_projects("PUT", "/projects/$(projectName)/placements/$(placementName)/devices/$(deviceTemplateName)", Dict{String, Any}("deviceId"=>deviceId); aws_config=aws_config)
associate_device_with_placement(deviceId, deviceTemplateName, placementName, projectName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_1click_projects("PUT", "/projects/$(projectName)/placements/$(placementName)/devices/$(deviceTemplateName)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("deviceId"=>deviceId), args)); aws_config=aws_config)

"""
    CreatePlacement()

Creates an empty placement.

# Required Parameters
- `placementName`: The name of the placement to be created.
- `projectName`: The name of the project in which to create the placement.

# Optional Parameters
- `attributes`: Optional user-defined key/value pairs providing contextual data (such as location or function) for the placement.
"""

create_placement(placementName, projectName; aws_config::AWSConfig=global_aws_config()) = iot_1click_projects("POST", "/projects/$(projectName)/placements", Dict{String, Any}("placementName"=>placementName); aws_config=aws_config)
create_placement(placementName, projectName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_1click_projects("POST", "/projects/$(projectName)/placements", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("placementName"=>placementName), args)); aws_config=aws_config)

"""
    CreateProject()

Creates an empty project with a placement template. A project contains zero or more placements that adhere to the placement template defined in the project.

# Required Parameters
- `projectName`: The name of the project to create.

# Optional Parameters
- `description`: An optional description for the project.
- `placementTemplate`: The schema defining the placement to be created. A placement template defines placement default attributes and device templates. You cannot add or remove device templates after the project has been created. However, you can update callbackOverrides for the device templates using the UpdateProject API.
- `tags`: Optional tags (metadata key/value pairs) to be associated with the project. For example, { {\"key1\": \"value1\", \"key2\": \"value2\"} }. For more information, see AWS Tagging Strategies.
"""

create_project(projectName; aws_config::AWSConfig=global_aws_config()) = iot_1click_projects("POST", "/projects", Dict{String, Any}("projectName"=>projectName); aws_config=aws_config)
create_project(projectName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_1click_projects("POST", "/projects", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("projectName"=>projectName), args)); aws_config=aws_config)

"""
    DeletePlacement()

Deletes a placement. To delete a placement, it must not have any devices associated with it.  When you delete a placement, all associated data becomes irretrievable.

# Required Parameters
- `placementName`: The name of the empty placement to delete.
- `projectName`: The project containing the empty placement to delete.

"""

delete_placement(placementName, projectName; aws_config::AWSConfig=global_aws_config()) = iot_1click_projects("DELETE", "/projects/$(projectName)/placements/$(placementName)"; aws_config=aws_config)
delete_placement(placementName, projectName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_1click_projects("DELETE", "/projects/$(projectName)/placements/$(placementName)", args; aws_config=aws_config)

"""
    DeleteProject()

Deletes a project. To delete a project, it must not have any placements associated with it.  When you delete a project, all associated data becomes irretrievable.

# Required Parameters
- `projectName`: The name of the empty project to delete.

"""

delete_project(projectName; aws_config::AWSConfig=global_aws_config()) = iot_1click_projects("DELETE", "/projects/$(projectName)"; aws_config=aws_config)
delete_project(projectName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_1click_projects("DELETE", "/projects/$(projectName)", args; aws_config=aws_config)

"""
    DescribePlacement()

Describes a placement in a project.

# Required Parameters
- `placementName`: The name of the placement within a project.
- `projectName`: The project containing the placement to be described.

"""

describe_placement(placementName, projectName; aws_config::AWSConfig=global_aws_config()) = iot_1click_projects("GET", "/projects/$(projectName)/placements/$(placementName)"; aws_config=aws_config)
describe_placement(placementName, projectName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_1click_projects("GET", "/projects/$(projectName)/placements/$(placementName)", args; aws_config=aws_config)

"""
    DescribeProject()

Returns an object describing a project.

# Required Parameters
- `projectName`: The name of the project to be described.

"""

describe_project(projectName; aws_config::AWSConfig=global_aws_config()) = iot_1click_projects("GET", "/projects/$(projectName)"; aws_config=aws_config)
describe_project(projectName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_1click_projects("GET", "/projects/$(projectName)", args; aws_config=aws_config)

"""
    DisassociateDeviceFromPlacement()

Removes a physical device from a placement.

# Required Parameters
- `deviceTemplateName`: The device ID that should be removed from the placement.
- `placementName`: The name of the placement that the device should be removed from.
- `projectName`: The name of the project that contains the placement.

"""

disassociate_device_from_placement(deviceTemplateName, placementName, projectName; aws_config::AWSConfig=global_aws_config()) = iot_1click_projects("DELETE", "/projects/$(projectName)/placements/$(placementName)/devices/$(deviceTemplateName)"; aws_config=aws_config)
disassociate_device_from_placement(deviceTemplateName, placementName, projectName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_1click_projects("DELETE", "/projects/$(projectName)/placements/$(placementName)/devices/$(deviceTemplateName)", args; aws_config=aws_config)

"""
    GetDevicesInPlacement()

Returns an object enumerating the devices in a placement.

# Required Parameters
- `placementName`: The name of the placement to get the devices from.
- `projectName`: The name of the project containing the placement.

"""

get_devices_in_placement(placementName, projectName; aws_config::AWSConfig=global_aws_config()) = iot_1click_projects("GET", "/projects/$(projectName)/placements/$(placementName)/devices"; aws_config=aws_config)
get_devices_in_placement(placementName, projectName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_1click_projects("GET", "/projects/$(projectName)/placements/$(placementName)/devices", args; aws_config=aws_config)

"""
    ListPlacements()

Lists the placement(s) of a project.

# Required Parameters
- `projectName`: The project containing the placements to be listed.

# Optional Parameters
- `maxResults`: The maximum number of results to return per request. If not set, a default value of 100 is used.
- `nextToken`: The token to retrieve the next set of results.
"""

list_placements(projectName; aws_config::AWSConfig=global_aws_config()) = iot_1click_projects("GET", "/projects/$(projectName)/placements"; aws_config=aws_config)
list_placements(projectName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_1click_projects("GET", "/projects/$(projectName)/placements", args; aws_config=aws_config)

"""
    ListProjects()

Lists the AWS IoT 1-Click project(s) associated with your AWS account and region.

# Optional Parameters
- `maxResults`: The maximum number of results to return per request. If not set, a default value of 100 is used.
- `nextToken`: The token to retrieve the next set of results.
"""

list_projects(; aws_config::AWSConfig=global_aws_config()) = iot_1click_projects("GET", "/projects"; aws_config=aws_config)
list_projects(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = iot_1click_projects("GET", "/projects", args; aws_config=aws_config)

"""
    ListTagsForResource()

Lists the tags (metadata key/value pairs) which you have assigned to the resource.

# Required Parameters
- `resourceArn`: The ARN of the resource whose tags you want to list.

"""

list_tags_for_resource(resourceArn; aws_config::AWSConfig=global_aws_config()) = iot_1click_projects("GET", "/tags/$(resourceArn)"; aws_config=aws_config)
list_tags_for_resource(resourceArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_1click_projects("GET", "/tags/$(resourceArn)", args; aws_config=aws_config)

"""
    TagResource()

Creates or modifies tags for a resource. Tags are key/value pairs (metadata) that can be used to manage a resource. For more information, see AWS Tagging Strategies.

# Required Parameters
- `resourceArn`: The ARN of the resouce for which tag(s) should be added or modified.
- `tags`: The new or modifying tag(s) for the resource. See AWS IoT 1-Click Service Limits for the maximum number of tags allowed per resource.

"""

tag_resource(resourceArn, tags; aws_config::AWSConfig=global_aws_config()) = iot_1click_projects("POST", "/tags/$(resourceArn)", Dict{String, Any}("tags"=>tags); aws_config=aws_config)
tag_resource(resourceArn, tags, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_1click_projects("POST", "/tags/$(resourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tags"=>tags), args)); aws_config=aws_config)

"""
    UntagResource()

Removes one or more tags (metadata key/value pairs) from a resource.

# Required Parameters
- `resourceArn`: The ARN of the resource whose tag you want to remove.
- `tagKeys`: The keys of those tags which you want to remove.

"""

untag_resource(resourceArn, tagKeys; aws_config::AWSConfig=global_aws_config()) = iot_1click_projects("DELETE", "/tags/$(resourceArn)", Dict{String, Any}("tagKeys"=>tagKeys); aws_config=aws_config)
untag_resource(resourceArn, tagKeys, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_1click_projects("DELETE", "/tags/$(resourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tagKeys"=>tagKeys), args)); aws_config=aws_config)

"""
    UpdatePlacement()

Updates a placement with the given attributes. To clear an attribute, pass an empty value (i.e., \"\").

# Required Parameters
- `placementName`: The name of the placement to update.
- `projectName`: The name of the project containing the placement to be updated.

# Optional Parameters
- `attributes`: The user-defined object of attributes used to update the placement. The maximum number of key/value pairs is 50.
"""

update_placement(placementName, projectName; aws_config::AWSConfig=global_aws_config()) = iot_1click_projects("PUT", "/projects/$(projectName)/placements/$(placementName)"; aws_config=aws_config)
update_placement(placementName, projectName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_1click_projects("PUT", "/projects/$(projectName)/placements/$(placementName)", args; aws_config=aws_config)

"""
    UpdateProject()

Updates a project associated with your AWS account and region. With the exception of device template names, you can pass just the values that need to be updated because the update request will change only the values that are provided. To clear a value, pass the empty string (i.e., \"\").

# Required Parameters
- `projectName`: The name of the project to be updated.

# Optional Parameters
- `description`: An optional user-defined description for the project.
- `placementTemplate`: An object defining the project update. Once a project has been created, you cannot add device template names to the project. However, for a given placementTemplate, you can update the associated callbackOverrides for the device definition using this API.
"""

update_project(projectName; aws_config::AWSConfig=global_aws_config()) = iot_1click_projects("PUT", "/projects/$(projectName)"; aws_config=aws_config)
update_project(projectName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_1click_projects("PUT", "/projects/$(projectName)", args; aws_config=aws_config)


text:
# ---
# title: 1540. Can Convert String in K Moves
# id: problem1540
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: String, Greedy
# link: <https://leetcode.com/problems/can-convert-string-in-k-moves/description/>
# hidden: true
# ---
#
# Given two strings `s` and `t`, your goal is to convert `s` into `t` in `k` **
# ** moves or less.
#
# During the `ith` (`1 <= i <= k`) move you can:
#
#   * Choose any index `j` (1-indexed) from `s`, such that `1 <= j <= s.length` and `j` has not been chosen in any previous move, and shift the character at that index `i` times.
#   * Do nothing.
#
# Shifting a character means replacing it by the next letter in the alphabet
# (wrapping around so that `'z'` becomes `'a'`). Shifting a character by `i`
# means applying the shift operations `i` times.
#
# Remember that any index `j` can be picked at most once.
#
# Return `true` if it's possible to convert `s` into `t` in no more than `k`
# moves, otherwise return `false`.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "input", t = "ouput", k = 9
#     Output: true
#     Explanation: In the 6th move, we shift 'i' 6 times to get 'o'. And in the 7th move we shift 'n' to get 'u'.
#
#
# **Example 2:**
#
#
#
#     Input: s = "abc", t = "bcd", k = 10
#     Output: false
#     Explanation: We need to shift each character in s one time to convert it into t. We can shift 'a' to 'b' during the 1st move. However, there is no way to shift the other characters in the remaining moves to obtain t from s.
#
#
# **Example 3:**
#
#
#
#     Input: s = "aab", t = "bbb", k = 27
#     Output: true
#     Explanation: In the 1st move, we shift the first 'a' 1 time to get 'b'. In the 27th move, we shift the second 'a' 27 times to get 'b'.
#
#
#
#
# **Constraints:**
#
#   * `1 <= s.length, t.length <= 10^5`
#   * `0 <= k <= 10^9`
#   * `s`, `t` contain only lowercase English letters.
#
#
## @lc code=start
using LeetCode

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


text:
"""
    Mlmnet(B, lambdas, alphas, data)

Type for storing the results of an mlmnet (Elastic net) model fit

"""

mutable struct Mlmnet

    # Coefficient estimates
    B::Array{Float64, 4}

    # Lambda and alpha penalties
    lambdas::Array{Float64, 1} # Total penalties
    alphas::Array{Float64, 1} # Penalty ratios

    # Response and predictor matrices
    data::RawData
end

"""
    mlmnet_pathwise(fun, X, Y, Z, lambdas, alphas, regXidx, regZidx,
                    reg, norms; isVerbose, stepsize, funArgs...)

Performs the supplied method on two descending lists of lambdas (for l1 and l2)
using ``warm starts''.

# Arguments

- fun = function that applies the Elastic-net pentalty estimate method
- X = 2d array of floats consisting of the row covariates, with all
  categorical variables coded in appropriate contrasts
- Y = 2d array of floats consisting of the multivariate response
- Z = 2d array of floats consisting of the column covariates, with all
  categorical variables coded in appropriate contrasts
- lambdas = 1d array of floats consisting of the total penalties in descending
  order. If they are not in descending order, they will be sorted.
- alphas = 1d array of floats consisting of the penalty ratio that
  determines the mix of penalties between L1 and L2
- regXidx = 1d array of indices corresponding to regularized X covariates
- regZidx = 1d array of indices corresponding to regularized Z covariates
- reg = 2d array of bits, indicating whether or not to regularize each of the
  coefficients
- norms = 2d array of floats consisting of the norms corresponding to each
  coefficient or `nothing`

# Keyword arguments

- isVerbose = boolean flag indicating whether or not to print messages.
  Defaults to `true`.
- stepsize = float; step size for updates. Defaults to `0.01`.
- funArgs = variable keyword arguments to be passed into `fun`

# Value

A 4d array consisting of the coefficient estimates, with the different
lambdas and alphas along the first and second dimensions respectively

# Some notes

Assumes that all necessary standardizations have been performed on X, Y, and
Z, including adding on intercepts. To be called by `mlmnet`, which performs
standardization and backtransforming.

"""
function mlmnet_pathwise(fun::Function, X::AbstractArray{Float64,2},
                         Y::AbstractArray{Float64,2},
                         Z::AbstractArray{Float64,2},
                         lambdas::AbstractArray{Float64,1},
                         alphas::AbstractArray{Float64, 1},
                         regXidx::AbstractArray{Int64,1},
                         regZidx::AbstractArray{Int64,1},
                         reg::BitArray{2}, norms; isVerbose::Bool=true,
                         stepsize::Float64=0.01, funArgs...)

    # Check that the lambdas are unique and in descending order.
    if length(lambdas) != length(unique(lambdas))
        println_verbose("Dropping non-unique lambdas", isVerbose)
    end

    # Pre-sort the lambdas in descending order in order to apply warm-starts
    if any(lambdas .!= sort(lambdas, rev=true))
        println_verbose("Sorting total penalty lambdas into descending order.", isVerbose)
        lambdas .= sort(lambdas, rev=true)
    end

    # Pre-allocate array for coefficients
    coeffs = Array{Float64}(undef, size(X,2), size(Z,2), length(lambdas), length(alphas))

    # Start with coefficients initalized at zero for the largest lambda value
    startB = zeros(size(X,2), size(Z,2))

    # Pre-compute eigenvalues and eigenvectors for ADMM
    # if length(string(fun)) > 7 && (string(fun)[(end-7):end] == "admm!")
    if string(fun) == "admm!"
        # Eigenfactorization of X
        XTX = transpose(X)*X
        eigX = eigen(XTX)
        Qx = eigX.vectors
        Lx = eigX.values

        # Eigenfactorization of Z
        ZTZ = transpose(Z)*Z
        eigZ = eigen(ZTZ)
		    Qz = eigZ.vectors
		    Lz = eigZ.values

        # Transformed Y
        X1 = X*Qx
        Z1 = Z*Qz
        U = transpose(X1) * Y * Z1

        # Kronecker product of eigenvalues of Z and X
        L = kron(Lx, transpose(Lz))
    end

    # Iterate through the paths of alphas and lambdas
    for i = 1:length(alphas)
      startB = zeros(size(X,2), size(Z,2))

      for j = 1:length(lambdas)

        # Get Elastic-net penalty estimates by updating the coefficients from previous
        # iteration in place

        # if length(string(fun)) <= 7 || (string(fun)[(end-7):end] != "admm!")
        if (string(fun) != "admm!" && string(fun) != "cd")
            # ISTA, FISTA and FISTA with Backtracking (CD not supported for Elastic-net yet)
            fun(X, Y, Z, lambdas[j], alphas[i], startB, regXidx, regZidx, reg, norms;
                isVerbose=isVerbose, stepsize=stepsize, funArgs...)

        else
            # ADMM
            fun(X, Y, Z, lambdas[j], alphas[i], startB, regXidx, regZidx, reg, norms,
                Qx, Qz, U, L;
                isVerbose=isVerbose, stepsize=stepsize, funArgs...)
        end

        # Assign a slice of coeffs to the current coefficient estimates
        coeffs[:, :, j, i] = startB
      end
    end

    return coeffs
end

"""
    mlmnet(data, lambdas, alphas;
           method = "ista",
           isNaive =false,
           hasXIntercept=true, hasZIntercept=true,
           toXReg=trues(data.p),
           toZReg=trues(data.q),
           toXInterceptReg::Bool=false, toZInterceptReg::Bool=false,
           toNormalize=true, isVerbose=true,
           stepsize=0.01, setStepsize=true,
           funArgs...)

Centers and normalizes X and Z predictor matrices, calculates fixed step size, performs
the supplied method on two descending lists of lambdas (each for L1 and L2) using ``warm starts'',
and backtransforms resulting coefficients, as is deemed necessary by the user
inputs.

# Arguments

- data = RawData object
- lambdas = 1d array of floats consisting of the total penalties in descending
  order. If they are not in descending order, they will be sorted.
- alphas = 1d array of floats consisting of the penalty ratio that
  determines the mix of penalties between L1 and L2

# Keyword arguments

- methods = function name that applies the Elastic-net penalty estimate method;
  default is `ista`, and the other methods are `fista`, `fista_bt`, `admm` and `cd`
- isNaive = boolean flag indicating whether to solve the Naive or non-Naive
  Elastic-net problem
- hasXIntercept = boolean flag indicating whether or not to include an `X`
  intercept (row main effects). Defaults to `true`.
- hasZIntercept = boolean flag indicating whether or not to include a `Z`
  intercept (column main effects). Defaults to `true`.
- toXReg = 1d array of bit flags indicating whether or not to regularize each
  of the `X` (row) effects. Defaults to 2d array of `true`s with length
  equal to the number of `X` effects (equivalent to `data.p`).
- toZReg = 1d array of bit flags indicating whether or not to regularize each
  of the `Z` (column) effects. Defaults to 2d array of `true`s with length
  equal to the number of `Z` effects (equivalent to `data.q`).
- toXInterceptReg = boolean flag indicating whether or not to regularize the
  `X` intercept Defaults to `false`.
- toZInterceptReg = boolean flag indicating whether or not to regularize the
  `Z` intercept. Defaults to `false`.
- toNormalize = boolean flag indicating if the columns of `X` and `Z`
  should be standardized (to mean 0, standard deviation 1). Defaults to `true`.
- isVerbose = boolean flag indicating whether or not to print messages.
  Defaults to `true`.
- stepsize = float; step size for updates (irrelevant for coordinate
  descent and when `setStepsize` is set to `true` for `ista!` and `fista!`).
  Defaults to `0.01`.
- setStepsize = boolean flag indicating whether the fixed step size should be
  calculated (for `ista!` and `fista!`). Defaults to `true`.
- funArgs = variable keyword arguments to be passed into `fun`

# Value

An Mlmnet object

# Some notes

The default method for choosing the fixed step size for `fista!` or `istaNet!`
is to use the reciprocal of the product of the maximum eigenvalues of
`X*transpose(X)` and `Z*transpose(Z)`. This is computed when `fista!` or
`ista!` is passed into the `fun` argument and `setStepsize` is set to `true`.
If `setStepsize` is set to `false`, the value of the `stepsize` argument will
be used as the fixed step size. Note that obtaining the eigenvalues when `X`
and/or `Z` are very large may exceed computational limitations.

Specifying a good starting step size (`stepsize`) and multiplying factor
(`gamma`) when `fista_bt!` is passed into the `fun` argument can be difficult.
Shrinking the step size too gradually can result in slow convergence. Doing so
too quickly can cause the criterion to diverge. We have found that setting
`stepsize` to 0.01 often works well in practice; choice of `gamma` appears to
be less consequential.

"""
function mlmnet(data::RawData,
                lambdas::AbstractArray{Float64,1}, alphas::AbstractArray{Float64,1};
                method::String = "ista",
                isNaive::Bool=false,
                hasXIntercept::Bool=true, hasZIntercept::Bool=true,
                toXReg::BitArray{1}=trues(data.p),
                toZReg::BitArray{1}=trues(data.q),
                toXInterceptReg::Bool=false, toZInterceptReg::Bool=false,
                toNormalize::Bool=true, isVerbose::Bool=true,
                stepsize::Float64=0.01, setStepsize::Bool=true,
                funArgs...)

    # Get the function according to the selected method
    # fun = get_func(method);

    try
      get_func(method);
    catch e
      return error("Unknown method name, possible methods are : \"ista\", \"fista\", \"fista_bt\", and \"admm\".")
    end
    fun = get_func(method);

    # Ensure that toXReg and toZReg have same length as columns of X and Z
    if length(toXReg) != data.p
        error("toXReg does not have same length as number of columns in X.")
    end
    if length(toZReg) != data.q
        error("toZReg does not have same length as number of columns in Z.")
    end

    # Add X and Z intercepts if necessary
    # Update toXReg and toZReg accordingly
    if hasXIntercept==true && data.predictors.hasXIntercept==false
        data.predictors.X = add_intercept(data.predictors.X)
        data.predictors.hasXIntercept = true
        data.p = data.p + 1
        toXReg = vcat(toXInterceptReg, toXReg)
    end
    if hasZIntercept==true && data.predictors.hasZIntercept==false
        data.predictors.Z = add_intercept(data.predictors.Z)
        data.predictors.hasZIntercept = true
        data.q = data.q + 1
        toZReg = vcat(toZInterceptReg, toZReg)
    end

    # Remove X and Z intercepts in new predictors if necessary
    # Update toXReg and toZReg accordingly
    if hasXIntercept==false && data.predictors.hasXIntercept==true
        data.predictors.X = remove_intercept(data.predictors.X)
        data.predictors.hasXIntercept = false
        data.p = data.p - 1
        toXReg = toXReg[2:end]
    end

    if hasZIntercept==false && data.predictors.hasZIntercept==true
        data.predictors.Z = remove_intercept(data.predictors.Z)
        data.predictors.hasZIntercept = false
        data.q = data.q - 1
        toZReg = toZReg[2:end]
    end

    # Update toXReg and toZReg accordingly when intercept is already included
    if hasXIntercept==true && data.predictors.hasXIntercept==true
        toXReg[1] = toXInterceptReg
    end
    if hasZIntercept==true && data.predictors.hasZIntercept==true
        toZReg[1] = toZInterceptReg
    end

    # Matrix to keep track of which coefficients to regularize.
    reg = toXReg.*transpose(toZReg)
    # Indices corresponding to regularized X covariates.
    regXidx = findall(toXReg)
    # Indices corresponding to regularized Z covariates.
    regZidx = findall(toZReg)

    # Centers and normalizes predictors, if necessary.
    if (toNormalize==true)
        # If predictors will be standardized, copy the predictor matrices.
        X = copy(get_X(data))
        Z = copy(get_Z(data))

        # Centers and normalizes predictors
        meansX, normsX, = normalize!(X, hasXIntercept)
        meansZ, normsZ, = normalize!(Z, hasZIntercept)
        # If X and Z are standardized, set the norm to nothing
        norms = nothing
    else
        # If not standardizing, create pointers for the predictor matrices
        X = data.predictors.X
        Z = data.predictors.Z

        # Calculate the norm matrix
        # 2d array of norms corresponding to each coefficient
        norms = transpose(sum(X.^2, dims=1)).*sum(Z.^2, dims=1)
    end

    # If chosen method is ista!/fista! with fixed step size and setStepsize is
    # true, compute the step size:
    # if length(string(fun)) > 7 && (string(fun)[(end-7):end] == "istaNet!") && #to delete
    #    setStepsize == true
    if (string(fun) == "ista!" || string(fun) == "fista!") &&
       setStepsize == true

        # Calculate and store transpose(X)*X
        XTX = transpose(X)*X
        # Calculate and store transpose(Z)*Z
        ZTZ = transpose(Z)*Z

        # Step size is the reciprocal of the maximum eigenvalue of kron(Z, X)
        if toNormalize==true
            # Standardizing X and Z results in complex eigenvalues

            # Hack is to square the singular values to get the eigenvalues
            eig_X = (svd(XTX).S).^2
            eig_Z = (svd(ZTZ).S).^2

            stepsize = 1/max(maximum(eig_X) * maximum(eig_Z),
                             minimum(eig_X) * minimum(eig_Z))
  	    else
            stepsize = 1/max(eigmax(XTX) * eigmax(ZTZ),
                             eigmin(XTX) * eigmin(ZTZ))
        end

        println_verbose(string("Fixed step size set to ", stepsize),
                        isVerbose)
    end

    # Run the specified Elastic-net penalty method on the supplied inputs.
    coeffs = mlmnet_pathwise(fun, X, get_Y(data), Z, lambdas, alphas, regXidx,
                             regZidx, reg, norms; isVerbose=isVerbose,
                             stepsize=stepsize, funArgs...)

    # Perform the direct scaling transformation to undo double-shrinkage
    # in the Naive Elastic-net solutions:
    if !isNaive
      for i in 1:length(alphas), j in 1:length(lambdas)
        lambdaL2 = lambdas[j]*(1-alphas[i])
        coeffs[:, :, j, i] *= (1+lambdaL2)
      end
    end

    # Back-transform coefficient estimates, if necessary.
    # Case if including both X and Z intercepts:
    if toNormalize == true && (hasXIntercept==true) && (hasZIntercept==true)
        backtransform!(coeffs, meansX, meansZ, normsX, normsZ, get_Y(data),
                       data.predictors.X, data.predictors.Z)
    elseif toNormalize == true # Otherwise
        backtransform!(coeffs, hasXIntercept, hasZIntercept, meansX, meansZ,
                       normsX, normsZ)
    end

    # lambdasL1 = lambdas.*alphas;
    # lambdasL2 = lambdas.*(1 .- alphas);

    # return Mlmnet(coeffs, lambdas, alphas, lambdasL1, lambdasL2, data)
    return Mlmnet(coeffs, lambdas, alphas, data)
  end

  """
  mlmnet(data, lambdas;
         method, hasXIntercept, hasZIntercept, toXReg, toZReg,
         toXInterceptReg, toZInterceptReg, toNormalize, isVerbose,
         stepsize, setStepsize, funArgs...)

"""

function mlmnet(data::RawData,
              lambdas::AbstractArray{Float64,1};
              method::String = "ista",
              isNaive::Bool=false,
              hasXIntercept::Bool=true, hasZIntercept::Bool=true,
              toXReg::BitArray{1}=trues(data.p),
              toZReg::BitArray{1}=trues(data.q),
              toXInterceptReg::Bool=false, toZInterceptReg::Bool=false,
              toNormalize::Bool=true, isVerbose::Bool=true,
              stepsize::Float64=0.01, setStepsize::Bool=true,
              funArgs...)

  alphas = [1.0] # default LASSO, 𝛼 = 1

  rslts = mlmnet(data, lambdas, alphas; method,
                    isNaive, hasXIntercept, hasZIntercept,
                    toXReg, toZReg,
                    toXInterceptReg, toZInterceptReg,
                    toNormalize, isVerbose,
                    stepsize, setStepsize,
                    funArgs...)

  return rslts
end


text:
#!/usr/bin/env julia

#=
  Copyright (c) 2021
  https://github.com/Ismael-VC/Varvara.jl/blob/main/CONTRIBUTORS.md

  Permission to use, copy, modify, and distribute this software for any
  purpose with or without fee is hereby granted, provided that the above
  copyright notice and this permission notice appear in all copies.

  THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
  WITH REGARD TO THIS SOFTWARE.
=#

function main()
  println("Running..")
  run(`julia utos.jl ../../projects/sounds/pad1.ss8 ../../projects/sounds/pad1.pcm`)
  println("Done.")
end

if abspath(PROGRAM_FILE) == @__FILE__
  main()
end


text:
module SpatialBoundaries

# Dependencies
using VoronoiDelaunay
using LinearAlgebra
using Statistics
using StatsBase
using Requires

include(joinpath("types", "Womble.jl"))
include(joinpath("types", "overloads.jl"))
export Womble, TriangulationWomble, LatticeWomble

include(joinpath("lib", "rateofchange.jl"))

include(joinpath("lib", "wombling.jl"))
export wombling

include(joinpath("lib", "boundaries.jl"))
export boundaries

function __init__()
    @require SimpleSDMLayers="2c645270-77db-11e9-22c3-0f302a89c64c" include(joinpath("extensions", "SimpleSDMLayers.jl"))
end

end


text:
# Hello world in Julia

println("Hello World")

max_stars_repo_name: hpissei/hacktoberfest

max_issues_repo_name: hpissei/hacktoberfest

max_forks_repo_name: hpissei/hacktoberfest


text:
struct Wildcard <: AbstractProcessingStep
    table::Observable{SimpleTable}
    card::ProcessingCard
end

function tryarguments(value::AbstractString)
    result = String[]
    try
        result = only(compute_calls(value)).positional
    catch e
    end
    return result
end

const WILDCARD_MODULES = (Base, Statistics, StatsBase) # TODO: consider adding more modules

function Wildcard(table::Observable{SimpleTable})
    inputs = RichTextField("Inputs", data_options(table, keywords=[""]), "")
    outputs = RichTextField("Outputs", data_options(table, keywords=[""], suffix="_new"), "")
    input_suggestions = lift(tryarguments, inputs.widget.value)
    output_suggestions = lift(tryarguments, outputs.widget.value)
    suggestions = Observable(
        [string(name) for m in WILDCARD_MODULES for name in names(m) if getproperty(m, name) isa Function]
    )
    # TODO: decide scores
    autocompleteoptions = Any[
        SimpleDict("meta" => "input", "words" => input_suggestions, "score" => 2),
        SimpleDict("meta" => "output", "words" => output_suggestions, "score" => 1),
        SimpleDict("meta" => "", "words" => suggestions, "score" => -1),
    ]
    wdgs = (;
        inputs,
        outputs,
        method=RichEditor("Method", "julia", "", autocompleteoptions),
    )
    card = ProcessingCard(:Wildcard; wdgs...)
    return Wildcard(table, card)
end

function (wc::Wildcard)(data)
    card = wc.card
    inputs = Symbol.(extract_positional_arguments(card.inputs))
    outputs = Symbol.(extract_positional_arguments(card.outputs))

    private_module = Module()
    for m in map(Symbol, WILDCARD_MODULES)
        @eval private_module using DashiBoard.($m)
    end
    for input in inputs
        @eval private_module $input = $(Tables.getcolumn(data, input))
    end
    value = wc.card.method.widget.value[]
    expr = Meta.parse("""
        begin
            $value
        end
    """)
    @eval private_module $expr
    return SimpleTable(
        outputs,
        getproperty.(Ref(private_module), outputs)
    )
end