text:
"""
$SIGNATURES

Clear the environment dictionaries.
"""
function clear_dicts()
    recursive_empty!.((GLOBAL_LXDEFS, GLOBAL_VARS, LOCAL_VARS, ALL_PAGE_VARS))
end

function recursive_empty!(d::AbstractDict)
    for (k, v) in d
        v isa AbstractDict && recursive_empty!(d[k])
    end
    empty!(d)
end

"""
$(SIGNATURES)

Runs Franklin in the current directory.

Keyword arguments:

* `clear=false`:            whether to remove any existing output directory
* `verb=false`:             whether to display messages
* `port=8000`:              the port to use for the local server (should pick a number
                            between 8000 and 9000)
* `single=false`:           whether to run a single pass or run continuously
* `nomess=false`:           suppresses all messages (internal use).
* `is_final_pass=false`:    whether we're in a "final pass" (if so, links are
                            fixed in case of a project website, see
                            [`convert_and_write`](@ref).
* `prerender=false`:        whether to pre-render javascript (KaTeX and highlight.js)
* `no_fail_prerender=true`: whether, in a prerendering phase, ignore errors and
                            try to produce an output
* `eval_all=false`:         whether to force re-evaluation of all code blocks
* `silent=false`:           switch this on to suppress all output (including eval
                            statements).
* `cleanup=true`:           whether to clear environment dictionaries, see
                            [`cleanup`](@ref).
* `on_write(pg, fd_vars)`:  callback function after the page is rendered,
                            passing as arguments the rendered page and the page
                            variables
* `host="127.0.0.1"`:       the host to use for the local server
* `show_warnings=true`:     whether to show franklin  warnings
* `launch=!single`:         whether to launch the browser when serving
"""
function serve(; clear::Bool             = false,
                 verb::Bool              = false,
                 port::Int               = 8000,
                 single::Bool            = false,
                 prerender::Bool         = false,
                 nomess::Bool            = false,
                 is_final_pass::Bool     = false,
                 no_fail_prerender::Bool = true,
                 eval_all::Bool          = false,
                 silent::Bool            = false,
                 cleanup::Bool           = true,
                 on_write::Function      = (_, _) -> nothing,
                 log::Bool               = false,
                 host::String            = "127.0.0.1",
                 show_warnings::Bool     = true,
                 fail_on_warning::Bool   = false,
                 launch::Bool            = !single,
                 no_set_paths::Bool      = false,
                 join_to_prepath::String = ""
                 )::Union{Nothing,Int}

    LOGGING[] = log
    # set the global path
    FOLDER_PATH[] = pwd()
    # silent mode?
    silent && (FD_ENV[:SILENT_MODE] = true; verb = false)

    if silent || !show_warnings
        FD_ENV[:SHOW_WARNINGS] = false
    end

    FD_ENV[:FAIL_ON_WARNING] = fail_on_warning

    # in case of optim, there may be a prepath given which should be
    # kept
    prepath = get(GLOBAL_VARS, "prepath", nothing)
    def_GLOBAL_VARS!()
    isnothing(prepath) || set_var!(GLOBAL_VARS, "prepath", prepath.first)

    # Set context out of kwargs vars
    FD_ENV[:CLEAR]      = clear
    FD_ENV[:VERB]       = verb
    FD_ENV[:FINAL_PASS] = is_final_pass
    FD_ENV[:PRERENDER]  = prerender
    FD_ENV[:NO_FAIL_PRERENDER] = no_fail_prerender
    FD_ENV[:ON_WRITE]   = on_write

    # check if there's a config file, if there is, check the variable
    # definitions looking at the ones that would affect overall structure etc.
    process_config()

    if !all(isdir, (joinpath(FOLDER_PATH[], "_layout"),
                    joinpath(FOLDER_PATH[], "_css")))
        throw(ArgumentError(
            "The current directory doesn't  have a `_layout` or `_css` " *
            "folder; change directory to a valid Franklin folder."))
    end

    # check if a Project.toml file is available, if so activate the folder
    flag_env = false
    if isfile(joinpath(FOLDER_PATH[], "Project.toml"))
        Pkg.activate(".")
        flag_env = true
    end

    # construct the set of files to watch
    watched_files = fd_setup(no_set_paths)

    # set a verbosity var that we'll use in the rest of the function
    nomess && (verb = false)

    # do a first full pass
    nomess || println("→ Initial full pass...")
    start = time()
    FD_ENV[:FORCE_REEVAL] = eval_all
    sig = fd_fullpass(watched_files, join_to_prepath)
    FD_ENV[:FORCE_REEVAL] = false
    sig < 0 && return sig
    fmsg = rpad("✔ full pass...", 40)
    verb && (println(""); print(fmsg); print_final(fmsg, start); println(""))

    # start the continuous loop
    if !single
        FD_ENV[:CLEAR] = false
        nomess || println("→ Starting the server...")
        coreloopfun = (cntr, fw) -> fd_loop(cntr, fw, watched_files)
        # start the liveserver in the current directory
        live_server_dir = "__site"
        LiveServer.setverbose(verb)
        LiveServer.serve(
            port=port,
            coreloopfun=coreloopfun,
            dir=live_server_dir,
            host=host,
            launch_browser=launch
        )
    end
    flag_env &&
        rprint("→ Use Pkg.activate() to go back to your main environment.")

    cleanup && clear_dicts()
    return nothing
end

"""
$(SIGNATURES)

Sets up the collection of watched files by doing an initial scan of the input
directory. It also sets the paths variables and prepares the output directory.

See also [`serve`](@ref).
"""
function fd_setup(no_set_paths::Bool=false)::NamedTuple
    # . setting up:
    # -- reading and storing the path variables
    # -- setting up the output directory (see `clear`)
    no_set_paths || set_paths!()
    prepare_output_dir()

    # . recovering the list of files in the input dir we care about
    # -- these are stored in dictionaries, the key is the full path and the
    # value is the time of last change (useful for continuous monitoring)
    md_pages         = TrackedFiles()
    html_pages       = TrackedFiles()
    other_files      = TrackedFiles()
    infra_files      = TrackedFiles()
    literate_scripts = TrackedFiles()
    # named tuples of all the watched files (order matters)
    watched_files = (other    = other_files,
                     infra    = infra_files,
                     md       = md_pages,
                     html     = html_pages,
                     literate = literate_scripts)
    # fill the dictionaries
    scan_input_dir!(watched_files...)
    return watched_files
end

"""
$(SIGNATURES)

A single full pass of Franklin looking at all watched files and processing them
as appropriate.

**Keyword arguments**

* `clear=false`:     whether to remove any existing output directory
* `verb=false`:      whether to display messages
* `is_final_pass=false` : whether it's the final pass before deployment
* `prerender=false`: whether to prerender katex and code blocks
* `no_fail_prerender=true`: whether to skip if a prerendering goes wrong in which
                            case don't prerender

See also [`fd_loop`](@ref), [`serve`](@ref) and [`publish`](@ref).
"""
function fd_fullpass(watched_files::NamedTuple, join_to_prepath::String="")::Int
    # keep track of context (some things either will or won't be done on
    # the full pass, e.g. see tag generation)
    FD_ENV[:FULL_PASS]  = true

    # reset global page variables and latex definitions
    # NOTE: need to keep track of pre-path if specified, see optimize
    prepath = get(GLOBAL_VARS, "prepath", "")
    def_GLOBAL_VARS!()
    def_GLOBAL_LXDEFS!()
    empty!.((RSS_ITEMS, SITEMAP_DICT))
    # reinsert prepath if specified
    isempty(prepath) || (GLOBAL_VARS["prepath"] = prepath)

    # process configuration file (see also `process_mddefs!`)
    # note the order (utils the config) is important, see also #774
    process_utils()
    process_config()

    if !isempty(join_to_prepath)
        set_var!(GLOBAL_VARS, "prepath",
                 joinpath(globvar("prepath"), join_to_prepath))
    end

    # form page segments
    root    = path(:folder)
    layout  = path(:layout)
    head    = read(joinpath(layout, "head.html"),      String)
    pg_foot = read(joinpath(layout, "page_foot.html"), String)
    foot    = read(joinpath(layout, "foot.html"),      String)

    # look for an index file to process
    hasindexmd   = isfile(joinpath(root, "index.md"))
    hasindexhtml = isfile(joinpath(root, "index.html"))

    if !(hasindexmd || hasindexhtml)
        FD_ENV[:SOURCE] = "full pass"
        print_warning("""
            No 'index.md' or 'index.html' found in the root directory. There
            should be one. Ignoring.
            """)
    end

    # if RSS is to be generated, check the template and the variables
    # the function will not do much if generate_rss is set to false
    prepare_for_rss()

    # go over all pages note that the html files are processed AFTER the
    # markdown files and so if you both have an `index.md` and an `index.html`
    # with otherwise the same path, it's the latter that will be considered.
    s = 0
    for (case, dict) ∈ pairs(watched_files), (fpair, t) ∈ dict
        # help keep track of what we're doing for debugging
        (:fd_fullpass, joinpath(fpair...)) |> logger
        # process
        a = process_file(case, fpair, head, pg_foot, foot, t)
        # in case of failure of prerendering, if no_fail_prerender, we force
        # prerender=false
        if a < 0 && FD_ENV[:PRERENDER] && FD_ENV[:NO_FAIL_PRERENDER]
            FD_ENV[:PRERENDER] = false
            process_file(case, fpair, head, pg_foot, foot, t)
            FD_ENV[:PRERENDER]  = true
        end
        s += a
    end

    # Finalize
    # > generate RSS if appropriate
    globvar("generate_rss")::Bool && rss_generator()
    # > generate tags if appropriate
    generate_tag_pages()
    # > generate sitemap if appropriate
    globvar("generate_sitemap")::Bool && sitemap_generator()
    # > generate robots if appropriate
    globvar("generate_robots")::Bool && robots_generator()

    # done
    FD_ENV[:FULL_PASS] = false

    # re-evaluate delayed pages
    if !isempty(DELAYED)
        cp_DELAYED = copy(DELAYED)
        empty!(DELAYED) # so that functions are effectively applied
        for page in cp_DELAYED
            case = Symbol(strip(splitext(page)[2], '.'))  # md / html
            fpair = path(:folder) => page
            a = process_file(case, fpair, head, pg_foot, foot)
            if a < 0 && FD_ENV[:PRERENDER] && FD_ENV[:NO_FAIL_PRERENDER]
                FD_ENV[:PRERENDER] = false
                process_file(case, fpair, head, pg_foot, foot)
                FD_ENV[:PRERENDER] = true
            end
            s += a
        end
    end

    # return -1 if any page failed to build, 0 otherwise
    return ifelse(s < 0, -1, 0)
end

"""
$(SIGNATURES)

This is the function that is continuously run, checks if files have been
modified and if so, processes them. Every 30 cycles, it checks whether any file
was added or deleted and consequently updates the `watched_files`.
"""
function fd_loop(
            cycle_counter::Int,
            ::LiveServer.FileWatcher,
            watched_files::NamedTuple
            )::Nothing
    verb = FD_ENV[:VERB]
    # every 30 cycles (3 seconds), scan directory to check for new or deleted
    # files and update dicts accordingly
    if mod(cycle_counter, 30) == 0
        # 1) check if some files have been deleted; note that we don't do
        # anything, we just remove the file reference from the corresponding
        # dictionary.
        for d ∈ watched_files, (fpair, _) ∈ d
            fpath = joinpath(fpair...)
            if !isfile(fpath)
                delete!(d, fpair)
                rp = splitext(get_rpath(fpath))[1]
                haskey(ALL_PAGE_VARS, rp) && delete!(ALL_PAGE_VARS, rp)
            end
        end
        # 2) scan the input folder, if new files have been added then this will
        # update the dictionaries
        scan_input_dir!(watched_files..., verb; in_loop=true)
    else
        layout = path(:layout)
        # do a pass over the files, check if one has changed and if so trigger
        # the appropriate file processing mechanism
        for (case, dict) ∈ pairs(watched_files), (fpair, t) ∈ dict
            # check if there was a modification to the file
            fpath = joinpath(fpair...)
            cur_t = mtime(fpath)
            cur_t <= t && continue
            # if there was then the file has been modified and should be
            # re-processed + copied
            fmsg = rpad("→ file $(fpath[length(FOLDER_PATH[])+1:end]) was " *
                        "modified ", 30)
            verb && print(fmsg)
            dict[fpair] = cur_t

            # if it's an infra_file trigger a fullpass as potentially
            # the whole website depends upon it
            # if it's a CSS file then don't however, that's directly read by browser
            if haskey(watched_files[:infra], fpair) && splitext(fpair.second)[2] != ".css"
                verb && println("→ full pass...")
                start = time()
                fd_fullpass(watched_files)
                if verb
                    print_final(rpad("✔ full pass...", 15), start)
                    println("")
                end

            # if it's a literate_file, check that it's included on some page
            # and if it is, then trigger build of that (or those) page(s)
            elseif haskey(watched_files[:literate], fpair)
                fmsg = fmsg * rpad("→ updating... ", 15)
                verb && print("\r" * fmsg)
                start = time()
                #
                literate_path = splitext(unixify(joinpath(fpair...)))[1]
                # retrieve head, foot etc
                head   = read(joinpath(layout, "head.html"),      String)
                pgfoot = read(joinpath(layout, "page_foot.html"), String)
                foot   = read(joinpath(layout, "foot.html"),      String)
                # process all md files that have `\literate` with something
                # that matches
                for mdfpair in keys(watched_files.md)
                    # read the content and look for `{...}` that may refer
                    # to that script.
                    content = read(joinpath(mdfpair...), String)
                    for m in eachmatch(r"\{(.*?)(\.jl)?\}", content)
                        if endswith(literate_path, m.captures[1])
                            process_file(:md, mdfpair, head, pgfoot, foot,
                                         cur_t)
                            # no need to look at further matches on that page
                            # since we've already triggered the page build
                            break
                        end
                    end
                end
                verb && print_final(fmsg, start)

            # for any other file, just call `process_file`
            else
                fmsg = fmsg * rpad("→ updating... ", 15)
                verb && print("\r" * fmsg)
                start = time()
                # retrieve head, foot etc
                head   = read(joinpath(layout, "head.html"),      String)
                pgfoot = read(joinpath(layout, "page_foot.html"), String)
                foot   = read(joinpath(layout, "foot.html"),      String)
                # then process
                process_file(case, fpair, head, pgfoot, foot, cur_t)
                verb && print_final(fmsg, start)
            end
        end
    end
    return nothing
end


text:
using ShallowWaters
using Test

@testset "No Forcing" begin
    Prog = RunModel(Ndays=1,Fx0=0)
    @test all(Prog.u .== 0.0f0)
    @test all(Prog.v .== 0.0f0)
    @test all(Prog.η .== 0.0f0)
end

@testset "Boundary Conditions" begin
    Prog = RunModel(Ndays=1,bc="periodic")
    @test all(abs.(Prog.η) .< 10)    # sea surface height shouldn't exceed +-10m

    Prog = RunModel(Ndays=1,bc="nonperiodic")
    @test all(abs.(Prog.η) .< 10)

    Prog = RunModel(Ndays=1,α=0)     # free-slip
    @test all(abs.(Prog.η) .< 10)
end

@testset "Continuity Forcing" begin
    Prog = RunModel(Ndays=1,surface_relax=true)
    @test all(abs.(Prog.η) .< 10)    # sea surface height shouldn't exceed +-10m
    @test all(Prog.u .!= 0.0f0)

    Prog = RunModel(Ndays=1,surface_forcing=true)
    @test all(abs.(Prog.η) .< 10)
    @test all(Prog.u .!= 0.0f0)
end

@testset "Mixed Precision" begin
    Prog = RunModel(Float32,Tprog=Float64,Ndays=1)
    @test all(abs.(Prog.η) .< 10)    # sea surface height shouldn't exceed +-10m
    @test all(Prog.u .!= 0.0f0)

    Prog = RunModel(Float16,Tprog=Float32,Ndays=1)
    @test all(abs.(Prog.η) .< 10)    # sea surface height shouldn't exceed +-10m
    @test all(Prog.u .!= 0.0f0)
end

ext: jl


text:
using Catalyst
using PWS
using StaticArrays
using Test

sn = @reaction_network begin
    κ, ∅ --> L
    λ, L --> ∅
end κ λ

xn = @reaction_network begin
    ρ, L --> L + X
    μ, X --> ∅
end ρ μ

u0 = SA[10, 20]
dtimes = 0:0.5:10.0
ps = [5.0, 1.0]
px = [3.0, 0.1]

system = PWS.SimpleSystem(sn, xn, u0, ps, px, dtimes)

algorithms = [DirectMCEstimate(128), SMCEstimate(128), TIEstimate(0, 4, 128)]
for algorithm in algorithms
    result = mutual_information(system, algorithm, num_responses = 10)
    for v in result[!, :MutualInformation]
        @test v[1] == 0
        @test length(v) == length(system.dtimes)
    end

    initial_condition = PWS.empirical_dist(rand(50, 50), 0:49, 0:49)

    system2 = PWS.SimpleSystem(sn, xn, initial_condition, ps, px, dtimes)
    result2 = mutual_information(system2, algorithm, num_responses = 10)
    for v in result2[!, :MutualInformation]
        @test v[1] != 0
        @test length(v) == length(system.dtimes)
    end
end


text:
module Hydroelastics

using RigidBodyDynamics
using LinearAlgebra
using Statistics
using StaticArrays
using GeometryBasics, Meshes
using CoordinateTransformations, MeshCat
using Polyhedra
using SpatialIndexing

const SI = SpatialIndexing

include("object.jl")
include("shapes.jl")
include("rendering.jl")
include("forces.jl")
include("simulation.jl")

end

max_stars_repo_name: ekzhang/hydroelastics

max_issues_repo_name: ekzhang/hydroelastics

max_forks_repo_name: ekzhang/hydroelastics


text:
### Utilities ###

struct ForcingParams{GF}
    field :: GF
    t0 :: Float64
    σt :: Float64
end

_process_forcing(::Nothing,s,grid) = nothing

function _process_forcing(params::ForcingParams,s,grid)
  gf = GeneratedField(s,params.field,grid)
  pulse = PulseField(gf,params.t0,params.σt)
  return [pulse]
end

function _process_forcing(params::Vector{<:ForcingParams},s,grid)
  pulses = ModulatedField[]
  for p in params
    gf = GeneratedField(s,p.field,grid)
    push!(pulses,PulseField(gf,p.t0,p.σt))
  end
  return pulses
end


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

# -----------------------------------------------------------------------------------------------
# GEOMETRY
# -----------------------------------------------------------------------------------------------
using ...OpticSim
using StaticArrays
using LinearAlgebra

#region Vec3

"""
`Vec3{T}` provides an immutable vector of fixed length 3 and type `T`.

`Vec3` defines a series of convenience constructors, so you can just type e.g. `Vec3(1, 2, 3)` or `Vec3([1.0, 2.0, 3.0])`.
It also supports comprehensions, and the `zeros()`, `ones()`, `fill()`, `rand()` and `randn()` functions, such as `Vec3(rand(3))`.
"""
Vec3{T} = SVector{3,T}
export Vec3

# empty constructor - initialized with zeros
Vec3(::Type{T} = Float64) where {T<:Real} = zeros(Vec3{T})

"""
returns the unit vector `[1, 0, 0]`
"""
unitX3(::Type{T} = Float64) where {T<:Real} = Vec3{T}(one(T), zero(T), zero(T))
"""
returns the unit vector `[0, 1, 0]`
"""
unitY3(::Type{T} = Float64) where {T<:Real} = Vec3{T}(zero(T), one(T), zero(T))
"""
returns the unit vector `[0, 0, 1]`
"""
unitZ3(::Type{T} = Float64) where {T<:Real} = Vec3{T}(zero(T), zero(T), one(T))

export unitX3, unitY3, unitZ3

#endregion Vec3

#region Vec4

"""
`Vec4{T}` provides an immutable vector of fixed length 4 and type `T`.

`Vec4` defines a series of convenience constructors, so you can just type e.g. `Vec3(1, 2, 3, 4)` or `Vec3([1.0, 2.0, 3.0, 4.0])`.
It also supports comprehensions, and the `zeros()`, `ones()`, `fill()`, `rand()` and `randn()` functions, such as `Vec4(rand(4))`.
"""
Vec4{T} = SVector{4,T}

# empty constructor - initialized with zeros
Vec4(::Type{T} = Float64) where {T<:Real} = zeros(Vec4{T})
export Vec4

"""
    Vec4(v::SVector{3, T}) where {T<:Real} -> Vec4{T}

Accept `SVector` and create a `Vec4` type [v[1], v[2], v[3], 1]
"""
function Vec4(v::SVector{3, T}) where {T<:Real}
    return Vec4{T}(v[1], v[2], v[3], one(T))
end

"""
    Vec4(m::SMatrix{3,N,T} where{N,T<:Real} -> SMatrix{3,N,T})
Input is matrix of 3d points, each column is one point. Returns matrix of 3d points with 1 appended in the last row.
"""
function Vec4(m::SMatrix{3,N,T}) where {N,T<:Real}
    return vcat(m,ones(SMatrix{1,N,T}))::SMatrix{4,N,T}
end
"""
returns the unit vector `[1, 0, 0, 0]`
"""
unitX4(::Type{T} = Float64) where {T<:Real} = Vec4{T}(one(T), zero(T), zero(T), zero(T))
"""
returns the unit vector `[0, 1, 0, 0]`
"""
unitY4(::Type{T} = Float64) where {T<:Real} = Vec4{T}(zero(T), one(T), zero(T), zero(T))
"""
returns the unit vector `[0, 0, 1, 0]`
"""
unitZ4(::Type{T} = Float64) where {T<:Real} = Vec4{T}(zero(T), zero(T), one(T), zero(T))
"""
returns the unit vector `[0, 0, 0, 1]`
"""
unitW4(::Type{T} = Float64) where {T<:Real} = Vec4{T}(zero(T), zero(T), zero(T), one(T))

export unitX4, unitY4, unitZ4, unitW4
#endregion Vec4

#region Transform

# utility function that given a vector, return 2 orthogonal vectors to that one
function get_orthogonal_vectors(direction::Vec3{T}) where {T<:Real}
    axis1 = normalize(direction)

    dp = dot(unitX3(T), axis1)
    # check the special case where the input vector is parallel to the X axis
    if	( dp >= one(T) - eps(T))
        # axis1 = unitX(T);
        axis2 = unitY3(T);
        axis3 = unitZ3(T);
        return (axis2, axis3)
    elseif  ( dp <= -one(T) + eps(T))
        # axis1 = -unitX(T);
        axis2 = -unitY3(T);
        axis3 = -unitZ3(T);
        return (axis2, axis3)
    end

    axis3 = normalize(cross(axis1, unitX3(T)))
    axis2 = normalize(cross(axis3, axis1))
    return (axis2, axis3)
end
export get_orthogonal_vectors

#---------------------------------------
# 3D Transform / Local Frame
#---------------------------------------
"""
    Transform{S<:Real}

Transform encapsulating rotation, translation and scale in 3D space. Translation happens **after** rotation.

```julia
Transform{S}(θ::T, ϕ::T, ψ::T, x::T, y::T, z::T)
Transform(rotation::SMatrix{3,3,S}, translation::SVector{3,S})
Transform(rotation::AbstractArray{S,2}, translation::AbstractArray{S,1})
```
`θ`, `ϕ` and `ψ` in first constructor are in **radians**.
"""
struct Transform{T}
    matrix::SMatrix{4,4,T,16}

    """ This is a private internal function. In general don't want to allow users to populate Transform matrices with arbitrary elements. Not to be called by code outside of the Transform module. Don't use Transform{Float64}(...) for example. Instead use Transform(..)"""
    function Transform{T}(a11::T,a21::T,a31::T,a41::T,a12::T,a22::T,a32::T,a42::T,a13::T,a23::T,a33::T,a43::T,a14::T,a24::T,a34::T,a44::T) where{T<:Real}
        return new{T}(SMatrix{4,4,T,16}(a11,a21,a31,a41,a12,a22,a32,a42,a13,a23,a33,a43,a14,a24,a34,a44))
    end

    Transform{T}(mat::SMatrix{4,4,T,16}) where{T<:Real} = new{T}(mat)
end
export Transform

#functions to make Transform compatible with base matrix API
matrix(a::Transform) = a.matrix

# Base.length(a::Transform) = length(matrix(a))
Base.getindex(a::Transform, indices::Vararg{Int,N}) where{N} = getindex(matrix(a),indices...)
# Base.iterate(a::Transform) = iterate(matrix(a))
# Base.iterate(a::Transform{Float64}, b::Tuple{StaticArrays.SOneTo{16}, Int64}) = iterate(matrix(a),b)

Base.collect(a::Transform) = collect(matrix(a))

Base.:*(transa::Transform{T},transb::Transform{T}) where{T<:Real} = Transform{T}(matrix(transa)*matrix(transb))

function Base.:*(transform::Transform{T}, v::SVector{3,S}) where {T<:Real,S<:Number}
    t = matrix(transform)

    res = t * Vec4(v)
    if (t[4,4] == one(T))
        return SVector(res[1], res[2], res[3])
    else
        return SVector(res[1]/res[4], res[2]/res[4], res[3]/res[4])
    end
end

#Transform is not necessarily constrained to be a rigid body transformation so use general invers.
Base.inv(a::Transform{T}) where{T<:Real}= Transform{T}(inv(matrix(a)))

""" The t and m matrices are allowed to be of different element type. This allows transforming a Unitful matrix for example:
```
id = identitytransform()
m = fill(1mm,3,4)
id*m #returns a matrix filled with Unitful quantities. If both matrices had to be the same type this would not work
```
"""
function Base.:*(transform::Transform{T}, m::SMatrix{3,N,S}) where{N,T<:Real,S<:Number}
    res = MMatrix{3,N,T}(undef)
    t = matrix(transform)

    for outcol in 1:N
        for row in 1:3
            sum = T(0)
            for incol in 1:3
                sum += t[row,incol]*m[incol,outcol]
            end
            #implicit 1 w coordinate value
            sum += t[row,4]
            res[row,outcol] = sum
        end
        if t[4,4] != 1
            res[:,outcol] /= t[4,4]
        end
    end
    return SMatrix{3,N,S}(res)
end

""" The t and m matrices are allowed to be of different element type. This allows transforming a Unitful matrix for example:
```
id = identitytransform()
m = fill(1mm,3,4)
id*m #returns a matrix filled with Unitful quantities. If both matrices had to be the same type this would not work
```
"""
Base.:*(transform::Transform{T},m::SMatrix{4,N,S}) where{N,T<:Real,S<:Number} = matrix(transform)*m

Base.transpose(a::Transform{T}) where{T<:Real} = Transform{T}(a')

# END of functions for compatibility with base matrix API

# for compatability ith the "old" RigidBodyTransform
"""
identitytransform([S::Type]) -> Transform{S}

Returns the [`Transform`](@ref) of type `S` (default `Float64`) representing the identity transform.
"""
identitytransform(::Type{T} = Float64) where {T<:Real} = Transform{T}(
    one(T), zero(T), zero(T), zero(T),
    zero(T), one(T), zero(T), zero(T),
    zero(T), zero(T), one(T), zero(T),
    zero(T), zero(T), zero(T), one(T)
)
export identitytransform

"""
    Transform([S::Type]) -> Transform{S}

Returns the [`Transform`](@ref) of type `S` (default `Float64`) representing the identity transform.
"""
function Transform(::Type{T} = Float64) where {T<:Real}
    return identitytransform(T)
end

"""
    Transform(colx::Vec3{T}, coly::Vec3{T},colz::Vec3{T}, colw::Vec3{T}, ::Type{T} = Float64) where {T<:Real}

Costruct a transform from the input columns.
"""
function Transform(colx::Vec3{T}, coly::Vec3{T}, colz::Vec3{T}, colw::Vec3{T} = zero(Vec3{T})) where {T<:Real}
    return Transform{T}(vcat(hcat(colx,coly,colz,colw),SMatrix{1,4,T}(zero(T),zero(T),zero(T),one(T)) ))
end

"""
    Transform(colx::Vec3{T}, coly::Vec3{T},colz::Vec3{T}, colw::Vec3{T}, ::Type{T} = Float64) where {T<:Real}

Costruct a transform from the input columns.
"""
function Transform(colx::Vec4{T}, coly::Vec4{T}, colz::Vec4{T}, colw::Vec4{T}) where {T<:Real}
    return Transform{T}(hcat(colx,coly,colz,colw))
end

"""
    Transform(origin, forward) -> Transform{S}

Returns the [`Transform`](@ref) of type `S` (default `Float64`) representing the local frame with origin and forward direction. the other 2 axes are computed automaticlly.
"""
function Transform(origin::Vec3{T}, forward::Vec3{T} = unitZ3()) where {T<:Real}
    forward = normalize(forward)
    right, up = get_orthogonal_vectors(forward)
    return Transform(right, up, forward, origin)
end

function Transform(θ::T, ϕ::T, ψ::T, x::T, y::T, z::T) where {T<:Number}
    return Transform(rotmat(T, θ, ϕ, ψ), Vec3{T}(x, y, z))
end

"""
    Transform(rotation::SMatrix{3,3,T}, translation::SVector{3,T}) where {T<:Real} -> Transform{S}

Returns the [`Transform`](@ref) of type `S` (default `Float64`) created by a rotation matrix and translation vector.
"""
function Transform(rotation::SMatrix{3,3,T}, translation::SVector{3,T}) where {T<:Real}
    return Transform{T}(
        rotation[1,1], rotation[2,1], rotation[3,1], zero(T),
        rotation[1,2], rotation[2,2], rotation[3,2], zero(T),
        rotation[1,3], rotation[2,3], rotation[3,3], zero(T),
        translation[1], translation[2], translation[3], one(T))
end

"""
    Transform(rotation::AbstractArray{T,2}, translation::AbstractArray{T,1}) where {T<:Real} -> Transform{S}

Returns the [`Transform`](@ref) of type `S` (default `Float64`) created by a rotation matrix (3x3) and translation vector of length 3.
"""
function Transform(rotation::AbstractArray{T,2}, translation::AbstractArray{T,1}) where {T<:Real}
    @assert size(rotation)[1] == size(rotation)[2] == length(translation) == 3
    return Transform{T}(
        rotation[1,1], rotation[2,1], rotation[3,1], zero(T),
        rotation[1,2], rotation[2,2], rotation[3,2], zero(T),
        rotation[1,3], rotation[2,3], rotation[3,3], zero(T),
        translation[1], translation[2], translation[3], one(T))
end

# define some utility functions

"""
    right(t::Transform{<:Real}) -> Vec3

Assuming t is a 3D rigid transform representing a local left-handed coordinate system, this function will return the first column, representing the "X" axis.
"""
right(t::Transform{<:Real}) = normalize(Vec3(t[1,1], t[2,1], t[3,1]))

"""
    up(t::Transform{<:Real}) -> Vec3

Assuming t is a 3D rigid transform representing a local left-handed coordinate system, this function will return the second column, representing the "Y" axis.
"""
up(t::Transform{<:Real}) = normalize(Vec3(t[1,2], t[2,2], t[3,2]))

"""
    forward(t::Transform{<:Real}) -> Vec3

Assuming t is a 3D rigid transform representing a local left-handed coordinate system, this function will return the third column, representing the "Z" axis.
"""
forward(t::Transform{<:Real}) = normalize(Vec3(t[1,3], t[2,3], t[3,3]))

"""
    origin(t::Transform{<:Real}) -> Vec3

Assuming t is a 3D rigid transform representing a local left-handed coordinate system, this function will return the fourth column, containing the translation part of the transform in 3D.
"""
OpticSim.origin(t::Transform{<:Real}) = Vec3(t[1,4], t[2,4], t[3,4])

export right, up, forward, origin

"""
    rotationX(angle::T) where {T<:Real} -> Transform

Builds a rotation matrix for a rotation around the x-axis.
Parameters:
    The counter-clockwise `angle` in radians.
"""
function rotationX(angle::T) where {T<:Real}
    c = cos(angle);
    s = sin(angle);

    row1 = unitX4()
    row2 = Vec4(zero(T), c, s, zero(T))
    row3 = Vec4(zero(T), -s, c, zero(T))
    row4 = unitW4()

    # transposing because the constructors treat these vectors as columns instead of rows
    return transpose(Transform(row1, row2, row3, row4))
end
export rotationX

"""
    rotationY(angle::T) where {T<:Real} -> Transform

Builds a rotation matrix for a rotation around the y-axis.
Parameters:
    The counter-clockwise `angle` in radians.
"""
function rotationY(angle::T) where {T<:Real}
    c = cos(angle);
    s = sin(angle);

    row1 = Vec4(c, zero(T), -s, zero(T))
    row2 = unitY4()
    row3 = Vec4(s, zero(T), c, zero(T))
    row4 = unitW4()

    # transposing because the constructors treat these vectors as columns instead of rows
    return transpose(Transform(row1, row2, row3, row4))
end
export rotationY

"""
    rotationZ(angle::T) where {T<:Real} -> Transform

Builds a rotation matrix for a rotation around the z-axis.
Parameters:
    The counter-clockwise `angle` in radians.
"""
function rotationZ(angle::T) where {T<:Real}
    c = cos(angle);
    s = sin(angle);

    row1 = Vec4(c, s, zero(T), zero(T))
    row2 = Vec4(-s, c, zero(T), zero(T))
    row3 = unitZ4()
    row4 = unitW4()

    # transposing because the constructors treat these vectors as columns instead of rows
    return transpose(Transform(row1, row2, row3, row4))
end
export rotationZ

"""
    rotation(t::Transform{T}) where {T<:Real} -> SMatrix{3,3,T}

returns the rotation part of the transform `t` - a 3x3 matrix.
"""
function rotation(t::Transform{T}) where {T<:Real}
    rot = SMatrix{3, 3, T}(
        t[1, 1], t[2, 1], t[3, 1],
        t[1, 2], t[2, 2], t[3, 2],
        t[1, 3], t[2, 3], t[3, 3])
    return Transform(rot, zero(Vec3{T}))
end

"""
    rotate(a::Transform{T}, vector::Union{Vec3{T}, SVector{3,T}}) where {T<:Real} -> Vec3{T}

apply the rotation part of the transform `a` to the vector `vector` - this operation is usually used to rotate direction vectors.
"""
rotate(a::Transform{T}, vector::Union{Vec3{T}, SVector{3,T}}) where {T<:Real} = rotation(a) * vector

"""
    rotation([S::Type], θ::T, ϕ::T, ψ::T) -> Transform{S}

Returns the [`Transform`](@ref) of type `S` (default `Float64`) representing the rotation by `θ`, `ϕ` and `ψ` around the *x*, *y* and *z* axes respectively **in radians**.
"""
rotation(θ::T, ϕ::T, ψ::T) where {T<:Number} = rotation(Float64, θ, ϕ, ψ)
rotation(::Type{S}, θ::T, ϕ::T, ψ::T) where {T<:Number,S<:Real} = Transform(rotmat(S, θ, ϕ, ψ), zeros(SVector{3,S}))
"""
    rotationd([S::Type], θ::T, ϕ::T, ψ::T) -> Transform{S}

Returns the [`Transform`](@ref) of type `S` (default `Float64`) representing the rotation by `θ`, `ϕ` and `ψ` around the *x*, *y* and *z* axes respectively **in degrees**.
"""
rotationd(θ::T, ϕ::T, ψ::T) where {T<:Number} = rotationd(Float64, θ, ϕ, ψ)
rotationd(::Type{S}, θ::T, ϕ::T, ψ::T) where {T<:Number,S<:Real} = Transform(rotmatd(S, θ, ϕ, ψ), zeros(SVector{3,S}))
export translation, rotation, rotationd
"""
    translation(x::T, y::T, z::T) where {T<:Real}

Creates a translation transform
"""
translation(::Type{S}, x::T, y::T, z::T) where {T<:Number,S<:Real} = convert(Transform{S},translation(x, y, z))

function translation(x::T, y::T, z::T) where {T<:Real}
    col1 = unitX4(T)
    col2 = unitY4(T)
    col3 = unitZ4(T)
    col4 = Vec4(x, y, z, one(T))

    return Transform(col1, col2, col3, col4)
end

"""
    translation(x::T, y::T, z::T) where {T<:Real}

Creates a translation transform
"""
function translation(t::Vec3{T}) where {T<:Real}
    return translation(t[1], t[2], t[3])
end
export translation

"""
    scale(x::T, y::T, z::T) where {T<:Real}

Creates a scaling transform
"""
function scale(x::T, y::T, z::T) where {T<:Real}
    col1 = unitX4(T) * x
    col2 = unitY4(T) * y
    col3 = unitZ4(T) * z
    col4 = unitW4(T)

    return Transform(col1, col2, col3, col4)
end

"""
    scale(s::T) where {T<:Real}

Creates a uniform scaling transform
"""
function scale(s::T) where {T<:Real}
    return scale(s, s, s)
end

"""
    scale(t::Vec3{T}) where {T<:Real}

Creates a scaling transform
"""
function scale(t::Vec3{T}) where {T<:Real}
    return scale(t[1], [2], [3])
end
export scale

"""
    local2world(t::Transform{T}) where {T<:Real}

return the transform matrix that takes a point in the local coordinate system to the global one
"""
function local2world(t::Transform{T}) where {T<:Real}
    return t
end
export local2world

"""
    world2local(t::Transform{T}) where {T<:Real}

return the transform matrix that takes a point in the global coordinate system to the local one
"""
function world2local(t::Transform{T}) where {T<:Real}
    return inv(t)
end
export world2local

"""
    decomposeRTS(tr::Transform{T}) where {T<:Real}

return a touple containing the rotation matrix, the translation vector and the scale vecto represnting the transform.
"""
function decomposeRTS(tr::Transform{T}) where {T<:Real}
    t = Vec3(tr[1,4], tr[2,4], tr[3,4])
    sx = norm(Vec3(tr[1,1], tr[2,1], tr[3,1]))
    sy = norm(Vec3(tr[1,2], tr[2,2], tr[3,2]))
    sz = norm(Vec3(tr[1,3], tr[2,3], tr[3,3]))
    s = Vec3(sx, sy, sz)
    rot = SMatrix{4, 4, T}(tr[1,1]/sx, tr[2, 1]/sx, tr[3,1]/sx, 0, tr[1,2]/sy, tr[2, 2]/sy, tr[3,2]/sy, 0, tr[1,3]/sz, tr[2, 3]/sz, tr[3,3]/sz, 0, 0, 0, 0, 1)

    return rot, t, s
end
export decomposeRTS

"""
    rotmatbetween([S::Type], a::SVector{3,T}, b::SVector{3,T}) -> SMatrix{3,3,S}

Returns the rotation matrix of type `S` (default `Float64`) representing the rotation between vetors `a` and `b`, i.e. rotation(a,b) * a = b.
"""
rotmatbetween(a::SVector{3,T}, b::SVector{3,T}) where {T<:Real} = rotmatbetween(Float64, a, b)
function rotmatbetween(::Type{S}, a::SVector{3,T}, b::SVector{3,T}) where {T<:Real,S<:Real}
    # TODO: Brian, is there a hidden assumption that a and b are normalized?
    v = cross(a, b)
    c = dot(a, b)
    V = SMatrix{3,3,S,9}(0, v[3], -v[2], -v[3], 0, v[1], v[2], -v[1], 0)
    R = I + V + V^2 * one(T) / (one(T) + c)
    return SMatrix{3,3,S,9}(R)
end

"""
    rotmatd([S::Type], θ::T, ϕ::T, ψ::T) -> SMatrix{3,3,S}

Returns the rotation matrix of type `S` (default `Float64`) representing the rotation by `θ`, `ϕ` and `ψ` around the *x*, *y* and *z* axes respectively **in degrees**.
"""
rotmatd(θ::T, ϕ::T, ψ::T) where {T<:Number} = rotmat(Float64, deg2rad(θ), deg2rad(ϕ), deg2rad(ψ))
rotmatd(::Type{S}, θ::T, ϕ::T, ψ::T) where {T<:Number,S<:Real} = rotmat(S, deg2rad(θ), deg2rad(ϕ), deg2rad(ψ))
export rotmatd
"""
    rotmat([S::Type], θ::T, ϕ::T, ψ::T) -> SMatrix{3,3,S}

Returns the rotation matrix of type `S` (default `Float64`) representing the rotation by `θ`, `ϕ` and `ψ` around the *x*, *y* and *z* axes respectively **in radians**.
"""
rotmat(θ::T, ϕ::T, ψ::T) where {T<:Number} = rotmat(Float64, θ, ϕ, ψ)
function rotmat(::Type{S}, θ::T, ϕ::T, ψ::T) where {T<:Number,S<:Real}
    sinψ = sin(ψ)
    sinϕ = sin(ϕ)
    sinθ = sin(θ)
    cosψ = cos(ψ)
    cosϕ = cos(ϕ)
    cosθ = cos(θ)
    return SMatrix{3,3,S,9}(cosψ * cosϕ, sinψ * cosϕ, -sinϕ, cosψ * sinϕ * sinθ - sinψ * cosθ, sinψ * sinϕ * sinθ + cosψ * cosθ, cosϕ * sinθ, cosψ * sinϕ * cosθ + sinψ * sinθ, sinψ * sinϕ * cosθ - cosψ * sinθ, cosϕ * cosθ)
end
export rotmat

#endregion Transform


text:
using AutomotiveDrivingModels
using AutoRisk
using HDF5
using NGSIM

function build_feature_extractor()
    subexts = [
        CoreFeatureExtractor(),
        TemporalFeatureExtractor(),
        WellBehavedFeatureExtractor(),
        CarLidarFeatureExtractor(20, carlidar_max_range = 50.),
        ForeForeFeatureExtractor()
    ]
    ext = MultiFeatureExtractor(subexts)
    return ext
end

function extract_features(
        ext,
        trajdata,
        roadway,
        timestep_delta,
        record_length,
        offset,
        prime,
        maxframes)
    n_features = length(ext)
    max_n_objects = maximum(n_objects_in_frame(trajdata, i) for i in 1 : nframes(trajdata))
    scene = Scene(max_n_objects)
    rec = SceneRecord(record_length, 0.1, max_n_objects)
    features = Dict{Int, Array{Float64}}()
    ctr = 0
    n_frames = nframes(trajdata)

    for frame in (offset - prime : offset - 1)
        # prime the rec
        AutomotiveDrivingModels.update!(rec, get!(scene, trajdata, frame))
    end

    veh_features = pull_features!(ext, rec, roadway, 1)

    for frame in offset : (n_frames - offset)
        ctr += 1
        if maxframes != nothing && ctr >= maxframes
            break
        end

        print("\rframe $(frame) / $(n_frames - offset)")

        # update the rec
        AutomotiveDrivingModels.update!(rec, get!(scene, trajdata, frame))

        # every timestep_delta step, extract features
        if frame % timestep_delta == 0

            for (vidx, veh) in enumerate(scene)
                # extract features
                veh_features = pull_features!(ext, rec, roadway, vidx)

                # add entry to features if vehicle not yet encountered
                if !in(veh.id, keys(features))
                    features[veh.id] = zeros(n_features, 0)
                end

                # stack onto existing features
                features[veh.id] = cat(2, features[veh.id],
                    reshape(veh_features, (n_features, 1)))
            end
        end
    end
    return features
end

function write_features(features, output_filepath, ext)
    n_features = length(ext)

    # compute max length across samples
    maxlen = 0
    for (traj_idx, feature_dict) in features
        for (veh_id, veh_features) in feature_dict
            maxlen = max(maxlen, size(veh_features, 2))
        end
    end
    println("max length across samples: $(maxlen)")

    # write trajectory features
    h5file = h5open(output_filepath, "w")
    for (traj_idx, feature_dict) in features

        feature_array = zeros(n_features, maxlen, length(feature_dict))
        for (idx, (veh_id, veh_features)) in enumerate(feature_dict)
            feature_array[:, 1:size(veh_features, 2), idx] = reshape(veh_features, (n_features, size(veh_features, 2), 1))
        end
        h5file["$(traj_idx)"] = feature_array

    end

    # write feature names
    attrs(h5file)["feature_names"] = feature_names(ext)
    close(h5file)
end

function extract_ngsim_features(
        timestep_delta = 1, # timesteps between feature extractions
        record_length = 10, # number of frames for record to track in the past
        offset = 500, # from ends of the trajectories
        prime = 10,
        maxframes = nothing; # nothing for no max
        output_filename = "ngsim.h5",
        n_expert_files = 1) # number of time periods for which to extract.

    ext = build_feature_extractor()
    features = Dict{Int, Dict{Int, Array{Float64}}}()

    tic()
    # extract
    for traj_idx in 1:n_expert_files

        # setup
        trajdata = load_trajdata(traj_idx)
        roadway = get_corresponding_roadway(traj_idx)
        features[traj_idx] = extract_features(
            ext,
            trajdata,
            roadway,
            timestep_delta,
            record_length,
            offset,
            prime,
            maxframes
        )
    end
    toc()

    output_filepath = joinpath("../data/trajectories/", output_filename)
    println("output filepath: $(output_filepath)")
    write_features(features, output_filepath, ext)

end

function extract_simple_features(
        filepath,
        output_filepath,
        timestep_delta = 1, # timesteps between feature extractions
        record_length = 20, # number of frames for record to track in the past
        offset = 3, # from ends of the trajectories
        prime = 2,
        maxframes = nothing) # nothing for no max

    ext = build_feature_extractor()
    features = Dict{Int, Dict{Int, Array{Float64}}}()

    tic()
    # setup
    trajdata = load_trajdata(filepath)
    roadway = get_corresponding_roadway(1)
    features[1] = extract_features(
        ext,
        trajdata,
        roadway,
        timestep_delta,
        record_length,
        offset,
        prime,
        maxframes
    )
    toc()
    write_features(features, output_filepath, ext)
end

# NGSIM
extract_ngsim_features(output_filename="ngsim_all.h5", n_expert_files=6)

# DEBUG
# trajdata_filepath = "/Users/wulfebw/.julia/v0.5/NGSIM/data/2_simple.txt"
# output_filepath = "../data/trajectories/2_simple.h5"
# extract_simple_features(trajdata_filepath, output_filepath)


text:
# ---
# title: 456. 132 Pattern
# id: problem456
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Stack
# link: <https://leetcode.com/problems/132-pattern/description/>
# hidden: true
# ---
#
# Given an array of `n` integers `nums`, a **132 pattern** is a subsequence of
# three integers `nums[i]`, `nums[j]` and `nums[k]` such that `i < j < k` and
# `nums[i] < nums[k] < nums[j]`.
#
# Return _`true` if there is a **132 pattern** in `nums`, otherwise, return
# `false`._
#
# **Follow up:** The `O(n^2)` is trivial, could you come up with the `O(n logn)`
# or the `O(n)` solution?
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [1,2,3,4]
#     Output: false
#     Explanation: There is no 132 pattern in the sequence.
#
#
# **Example 2:**
#
#
#
#     Input: nums = [3,1,4,2]
#     Output: true
#     Explanation: There is a 132 pattern in the sequence: [1, 4, 2].
#
#
# **Example 3:**
#
#
#
#     Input: nums = [-1,3,2,0]
#     Output: true
#     Explanation: There are three 132 patterns in the sequence: [-1, 3, 2], [-1, 3, 0] and [-1, 2, 0].
#
#
#
#
# **Constraints:**
#
#   * `n == nums.length`
#   * `1 <= n <= 104`
#   * `-109 <= nums[i] <= 109`
#
#
## @lc code=start
using LeetCode

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


text:
module Utils

export @get, MappedVector, centraldiff

"""
A macro for extracting fields from an object.  For example, instead of a statement
like

    (obj.a + obj.b)^obj.c

it is more readable if we first locally bind the variables

    a = obj.a
    b = obj.b
    c = obj.c

    (a + b)^c

Using the `@get` marco, this becomes

    @get obj (a, b, c)

    (a + b)^c

We can also locally assign different names to the binding

    @get obj (a, b, c) (α, β, γ)
    (α + β)^γ
"""
macro get(object, fields...)
    if length(fields) == 1
        try
            @assert typeof(fields[1]) == Expr
            @assert fields[1].head == :tuple
            @assert all([typeof(arg) == Symbol for arg in fields[1].args])
        catch
            throw(ArgumentError("second argument must be a tuple of field names"))
        end
        esc(Expr(:block, [:($sym = $object.$sym) for sym in fields[1].args]...))
    elseif length(fields) == 2
        try
            @assert typeof(fields[1]) == Expr
            @assert typeof(fields[2]) == Expr
            @assert (fields[1].head == :tuple && fields[2].head == :tuple)
            @assert all([typeof(arg) == Symbol for arg in fields[1].args])
            @assert all([typeof(arg) == Symbol for arg in fields[2].args])
        catch
            throw(ArgumentError("second and third argument must be tuples of field names"))
        end

        nargs = length(fields[1].args)
        @assert nargs == length(fields[2].args) "field name tuples must have the same length"
        esc(Expr(:block, [:($(fields[2].args[i]) = $object.$(fields[1].args[i])) for i in 1:nargs]...))
    else
        throw(ArgumentError("Usage: @get <object> (field names...) [(reference names...)]"))
    end
end

# Based on Tim Holy's JuliaCon 2016 Keynote
"""
A wrapper around an array that applies a function to any index element

# Fields
- `f`: the function to apply
- `data`: the underlying array
- `offset`: an optional offset (0 by default)

# Constructor

    MappedVector(f, data::AbstractVector{T}, T₀ = typeof(f(one(T))), offset = 0) where {T}

where `T₀` should be the type of the mapped values.

# Examples:

```julia
julia> x = [-π, 0.0, π];

julia> y = MappedVector(cos, x, Float64, 1) # apply `cos` on demand and start indexing from 0
Array{Float64,1} → Base.#cos

julia> y[0]
-1.0
```
"""
struct MappedVector{T,A <: AbstractVector,F} <: AbstractVector{T}
    f::F
    data::A
    offset::Int
end

Base.size(A::MappedVector) = size(A.data)
Base.@propagate_inbounds Base.getindex(A::MappedVector, i::Int) = A.f(A.data[i + A.offset])

function MappedVector(f, data::AbstractVector{T}, offset = 0) where {T}
    T₀ = typeof(f(one(T)))
    A  = typeof(data)
    F  = typeof(f)

    MappedVector{T₀, A, F}(f, data, offset)
end

function Base.show(io::IO, M::MIME"text/plain", m::MappedVector{T, A, F}) where {T, A, F}
    print(io, "$A → $F ($(1-m.offset):$(length(m.data)-m.offset))")
end
Base.show(io::IO, m::MappedVector) = Base.show(io, MIME("text/plain"), m)

# Base.gradient stopped working, so this is the fix
function centraldiff(F::AbstractVector, h::Vector)
    n = length(F)
    T = typeof(oneunit(eltype(F))/oneunit(eltype(h)))
    g = similar(F, T)
    if n == 1
        g[1] = zero(T)
    elseif n > 1
        g[1] = (F[2] - F[1]) / (h[2] - h[1])
        g[n] = (F[n] - F[n-1]) / (h[end] - h[end-1])
        if n > 2
            h = h[3:n] - h[1:n-2]
            g[2:n-1] = (F[3:n] - F[1:n-2]) ./ h
        end
    end
    g
end

centraldiff(F::AbstractVector) = centraldiff(F, [1:length(F);])

centraldiff(F::AbstractVector, h::Real) = centraldiff(F, [h*(1:length(F));])

include("autodiff/forwarddiff.jl")

end


text:
#=

JuliaBolt is ported from Neobolt: Neo4j Bolt connector for Python
https://github.com/neo4j-drivers/neobolt

Original notices:

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

=#

#####################################
#             Structure             #
#####################################

mutable struct Structure
    tag::UInt8
    fields::Array{Any}

    Structure(tag, fields) = new(tag, fields)
    Structure(tag) = new(tag, [])
end

function Base.getindex(structure::Structure, i::Int)
    return Base.getindex(structure.fields, i)
end

function Base.setindex(structure::Structure, v, i::Int)
    return Base.setindex(structure.fields, v, i)
end

function Base.firstindex(structure::Structure)
    return Base.firstindex(structure.fields)
end

function Base.lastindex(structure::Structure)
    return Base.lastindex(structure.fields)
end

# Structure ----------------------- #

Base.:(==)(x::Structure, y::Structure) = x.tag == y.tag && x.fields == y.fields

# Structure ----------------------- #
# # # # # # # # # # # # # # # # # # #

#####################################
#              strpack              #
#####################################

function struct_pack(value::UInt8)
    return [hton(value)]
end

# strpack ------------------------- #

function struct_pack(value::Union{UInt16, UInt32, Int16, Int32, Int64, Float64})
    return copy(reinterpret(UInt8, [hton(value)]))
end

# strpack ------------------------- #

function struct_unpack(t::DataType, v::AbstractArray{UInt8})
    _value = reinterpret(t, v)
    return ntoh(_value[1])
end

# strpack ------------------------- #
# # # # # # # # # # # # # # # # # # #

include("packer.jl")
include("unpacker.jl")

max_stars_repo_name: virtualgraham/JuliaBolt

max_issues_repo_name: virtualgraham/JuliaBolt

max_forks_repo_name: virtualgraham/JuliaBolt


text:
abstract type PlasticityBounds end

struct PlasticityBoundsNone <: PlasticityBounds  end

struct PlasticityBoundsNonnegative <: PlasticityBounds end

struct PlasticityBoundsLowHigh <: PlasticityBounds
  low::Float64
  high::Float64
end

function (plb::PlasticityBoundsNonnegative)(w::R,Δw::R) where R
  ret = w+Δw
  return min(0.0,ret)
end
function (plb::PlasticityBoundsLowHigh)(w::R,Δw::R) where R
  ret = w+Δw
  return min(plb.high,max(plb.low,ret))
end

# no plasticity, nothing to do
@inline function plasticity_update!(::Real,::Real,
     ::PSSpikingType,::Connection,::PSSpikingType,::NoPlasticity)
  return nothing
end

# Pairwise plasticity
struct PairSTDP <: PlasticityRule
  τplus::Float64
  τminus::Float64
  Aplus::Float64
  Aminus::Float64
  o::Vector{Float64} # pOst
  r::Vector{Float64} # pRe
  bounds::PlasticityBounds
  function PairSTDP(τplus,τminus,Aplus,Aminus,n_post,n_pre;
       plasticity_bounds=PlasticityBoundsNonnegative())
    new(τplus,τminus,Aplus,Aminus,zeros(n_post),zeros(n_pre),plasticity_bounds)
  end
end

function reset!(pl::PairSTDP)
  fill!(pl.r,0.0)
  fill!(pl.o,0.0)
  return nothing
end

function plasticity_update!(::Real,dt::Real,
     pspost::PSSpikingType,conn::Connection,pspre::PSSpikingType,
     plast::PairSTDP)
  # elements of sparse matrix that I need
  _colptr = SparseArrays.getcolptr(conn.weights) # column indexing
	row_idxs = rowvals(conn.weights) # postsynaptic neurons
	weightsnz = nonzeros(conn.weights) # direct access to weights
  idx_pre_spike = findall(pspre.isfiring)
  idx_post_spike = findall(pspost.isfiring)
  # update synapses
  # presynpatic spike go along w column
  for j_pre in idx_pre_spike
		_posts_nz = nzrange(conn.weights,j_pre) # indexes of corresponding pre in nz space
		@inbounds for _pnz in _posts_nz
      ipost = row_idxs[_pnz]
      Δw = -plast.o[ipost]*plast.Aminus
      weightsnz[_pnz] = max(0.0,weightsnz[_pnz]-Δw)
    end
  end
  # postsynaptic spike: go along w row
  # innefficient ... need to search i element for each column
  for i_post in idx_post_spike
    for j_pre in (1:pspre.n)
      _start = _colptr[j_pre]
      _end = _colptr[j_pre+1]-1
      _pnz = searchsortedfirst(row_idxs,i_post,_start,_end,Base.Order.Forward)
      if (_pnz<=_end) && (row_idxs[_pnz] == i_post) # must check!
        Δw = plast.r[j_pre]*plast.Aplus
        weightsnz[_pnz]+= Δw
      end
    end
  end
  # update the plasticity trace variables
  for j_pre in idx_pre_spike
    plast.r[j_pre]+=1.0
  end
  for i_post in idx_post_spike
    plast.o[i_post]+=1.0
  end
  @. plast.r -= plast.r*dt/plast.τplus
  @. plast.o -= plast.o*dt/plast.τminus
  return nothing
end

### Triplets rule
struct PlasticityTriplets <: PlasticityRule
  τplus::Float64
  τminus::Float64
  τx::Float64
  τy::Float64
  A2plus::Float64
  A3plus::Float64
  A2minus::Float64
  A3minus::Float64
  o1::Vector{Float64} # pOst
  o2::Vector{Float64} # pOst
  r1::Vector{Float64} # pRe
  r2::Vector{Float64} # pRe
  bounds::PlasticityBounds
  function PlasticityTriplets(τplus,τminus,τx,τy,
      A2plus,A3plus,A2minus,A3minus,n_post,n_pre;
        plasticity_bounds=PlasticityBoundsNonnegative())
    new(τplus,τminus,τx,τy,A2plus,A3plus,A2minus,A3minus,
        ntuple(_->zeros(n_post),2)...,
        ntuple(_->zeros(n_pre),2)...,plasticity_bounds)
  end
end

function reset!(pl::PlasticityTriplets)
  fill!(pl.r1,0.0)
  fill!(pl.r2,0.0)
  fill!(pl.o1,0.0)
  fill!(pl.o2,0.0)
  return nothing
end

function plasticity_update!(t_now::Real,dt::Real,
     pspost::PSSpikingType,conn::Connection,pspre::PSSpikingType,
     plast::PlasticityTriplets)
  # elements of sparse matrix that I need
  _colptr = SparseArrays.getcolptr(conn.weights) # column indexing
	row_idxs = rowvals(conn.weights) # postsynaptic neurons
	weightsnz = nonzeros(conn.weights) # direct access to weights
  idx_pre_spike = findall(pspre.isfiring)
  idx_post_spike = findall(pspost.isfiring)
  # update the plasticity traces 1
  for j_pre in idx_pre_spike
    plast.r1[j_pre]+=1.0
  end
  for i_post in idx_post_spike
    plast.o1[i_post]+=1.0
  end
  # update synapses
  # presynpatic spike go along w column
  for j_pre in idx_pre_spike
		_posts_nz = nzrange(conn.weights,j_pre) # indexes of corresponding pre in nz space
		@inbounds for _pnz in _posts_nz
      i_post = row_idxs[_pnz]
      Δw = -plast.o1[i_post]*(plast.A2minus+plast.A3minus*plast.r2[j_pre])
      weightsnz[_pnz] = plast.bounds(weightsnz[_pnz],Δw)
    end
  end
  # postsynaptic spike: go along w row
  # innefficient ... need to search i element for each column
  for i_post in idx_post_spike
    for j_pre in (1:pspre.n)
      _start = _colptr[j_pre]
      _end = _colptr[j_pre+1]-1
      _pnz = searchsortedfirst(row_idxs,i_post,_start,_end,Base.Order.Forward)
      if (_pnz<=_end) && (row_idxs[_pnz] == i_post) # must check!
        Δw = plast.r1[j_pre]*(plast.A2plus+plast.A3plus*plast.o2[i_post])
        weightsnz[_pnz] = plast.bounds(weightsnz[_pnz],Δw)
      end
    end
  end
  # update the plasticity traces 2
  for j_pre in idx_pre_spike
    plast.r2[j_pre]+=1.0
  end
  for i_post in idx_post_spike
    plast.o2[i_post]+=1.0
  end
  # and timestep all
  @. plast.r1 -= plast.r1*dt/plast.τplus
  @. plast.r2 -= plast.r2*dt/plast.τx
  @. plast.o1 -= plast.o1*dt/plast.τminus
  @. plast.o2 -= plast.o2*dt/plast.τy
  return nothing
end

# Vogels et al 2011 plasticity rule

struct PlasticityInhibitoryVogels <: PlasticityRule
  τ::Float64
  η::Float64
  α::Float64
  o::Vector{Float64} # pOst
  r::Vector{Float64} # pRe
  bounds::PlasticityBounds
  function PlasticityInhibitoryVogels(τ,η,n_post,n_pre;
      r_target=5.0,
      plasticity_bounds::PlasticityBounds=PlasticityBoundsNonnegative())
    α = 2*r_target*τ
    new(τ,η,α,zeros(n_post),zeros(n_pre),plasticity_bounds)
  end
end
function reset!(pl::PlasticityInhibitoryVogels)
  fill!(pl.r,0.0)
  fill!(pl.o,0.0)
  return nothing
end

function plasticity_update!(::Real,dt::Real,
     pspost::PSSpikingType,conn::Connection,pspre::PSSpikingType,
     plast::PlasticityInhibitoryVogels)
  # elements of sparse matrix that I need
  _colptr = SparseArrays.getcolptr(conn.weights) # column indexing
	row_idxs = rowvals(conn.weights) # postsynaptic neurons
	weightsnz = nonzeros(conn.weights) # direct access to weights
  idx_pre_spike = findall(pspre.isfiring)
  idx_post_spike = findall(pspost.isfiring)
  # update traces first  (mmmh)
  for j_pre in idx_pre_spike
    plast.r[j_pre]+=1.0
  end
  for i_post in idx_post_spike
    plast.o[i_post]+=1.0
  end
  # update synapses
  # presynpatic spike go along w column
  for j_pre in idx_pre_spike
		_posts_nz = nzrange(conn.weights,j_pre) # indexes of corresponding pre in nz space
		@inbounds for _pnz in _posts_nz
      ipost = row_idxs[_pnz]
      Δw = plast.η*(plast.o[ipost]-plast.α)
      #@show plast.o[ipost]
      weightsnz[_pnz] = plast.bounds(weightsnz[_pnz],Δw)
    end
  end
  # postsynaptic spike: go along w row
    # innefficient ... need to search i element for each column
  for i_post in idx_post_spike
    for j_pre in (1:pspre.n)
      _start = _colptr[j_pre]
      _end = _colptr[j_pre+1]-1
      _pnz = searchsortedfirst(row_idxs,i_post,_start,_end,Base.Order.Forward)
      if (_pnz<=_end) && (row_idxs[_pnz] == i_post) # must check!
        Δw = plast.η*plast.r[j_pre]
        #@show Δw
        weightsnz[_pnz] = plast.bounds(weightsnz[_pnz],Δw)
      end
    end
  end
  @. plast.r -= plast.r*dt/plast.τ
  @. plast.o -= plast.o*dt/plast.τ
  return nothing
end


text:
include("../src/PowerFlowNetworks.jl")
using .PowerFlowNetworks
using SQLite
using ArgParse

ArgParse.parse_item(::Type{Vector{String}}, x::AbstractString) = split(x, ",")

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "--dbpath", "-p"
            help = "Path of the sqlite database file."
            default = "data/PowerFlowNetworks.sqlite"
        "--extension_alg"
            help = "Algorithm used for the chordal extension."
            default = "cholesky"
        "--cliques_path"
            help = "Directory where to store the cliques."
            default = "data/cliques"
        "--cliquetrees_path"
            help = "Directory where to store the cliquetrees."
            default = "data/cliquetrees"
        "--graphs_path"
            help = "Where to store the graphs."
            default = "data/graphs/"
        "--preprocess_path"
            help = "JSON file containing the preprocess option."
            default = "data/configs/preprocess_default.json"
        "--min_nv"
            help = "Minimum number of vertices a network had to have to be processed."
            arg_type = Int
            default = typemin(Int)
        "--max_nv"
            help = "Maximum number of vertices a network had to have to be processed."
            arg_type = Int
            default = typemax(Int)
    end
    return parse_args(s)
end

function main()
    parsed_args = parse_commandline()
    db = SQLite.DB(parsed_args["dbpath"])
    generate_decompositions!(db;
                             cliques_path=parsed_args["cliques_path"],
                             cliquetrees_path=parsed_args["cliquetrees_path"],
                             graphs_path=parsed_args["graphs_path"],
                             extension_alg=parsed_args["extension_alg"],
                             preprocess_path=parsed_args["preprocess_path"],
                             min_nv=parsed_args["min_nv"],
                             max_nv=parsed_args["max_nv"])
end

main()

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks