text:
module KnowledgeGraphNote

using LightGraphs
using UUIDs
using StringDistances

struct Concept
    name::String
    link::String
    dependency::Array{String, 1}
    text::String
end

struct KnowledgeGraph
    graph::SimpleDiGraph
    nametoid::Dict{String,Int} # normalized concept name => vertex id in graph
    idtoname::Array{String} # vertex id in graph => normalized concept name
end

function print_concepts(concepts::Array{Concept, 1})
    if isnothing(concepts)
        return
    end
    for concept in concepts
        println("name: $(concept.name)")
        println("link: $(concept.link)")
        if length(concept.dependency) > 0
            println("prerequisites: $(concept.dependency)")
        end
        println(concept.text)
        println("")
    end
end

normalize_concept_name(name) = lowercase(strip(name))

get_concept_index(concepts) = Dict([(normalize_concept_name(concept.name), findfirst(x->x.name == concept.name, concepts)) for concept in concepts])

function get_digraph(concepts)
    g = SimpleDiGraph(length(concepts))
    concept_index = get_concept_index(concepts)
    for (index, concept) in enumerate(concepts)
        if length(concept.dependency) > 0
            for prerequisite in concept.dependency
                nprerequisite = normalize_concept_name(prerequisite)
                if haskey(concept_index, nprerequisite)
                    # ignore missing concepts for now
                    add_edge!(g, index, concept_index[nprerequisite])
                end
            end
        end
    end
    return (g, concept_index, [normalize_concept_name(concept.name) for concept in concepts])
end

# use an struct to hold the knowledge graph in memory so we don't have to recreate it as long as the concepts are not changed.
# it is advised to check duplicate concepts before creating knowledge graph.
function init_knowledge_graph(concepts::Array{Concept, 1})
    g, nametoid, idtoname = get_digraph(concepts)
    return KnowledgeGraph(g, nametoid, idtoname)
end

get_known_concepts(concepts) = unique(normalize_concept_name.(getfield.(concepts, :name)))

function get_duplicate_concepts(concepts::Array{Concept, 1})
    known_concepts = get_known_concepts(concepts)
    concept_occurrences = [(concept,count(c->normalize_concept_name(c.name)==concept, concepts)) for concept in known_concepts]
    duplicate_occurrences = Dict(filter(item -> item[2] > 1, concept_occurrences))
    return duplicate_occurrences
end

function get_missing_concepts(concepts::Array{Concept, 1})
    known_concepts = get_known_concepts(concepts)
    dependencies = reduce(concepts; init = Set()) do acc, concept
        acc = union(acc, Set(concept.dependency))
    end
    ret = setdiff(normalize_concept_name.(dependencies), known_concepts)
    return ret
end

function find_cycles(kg::KnowledgeGraph)
    #    savegraph("test.gml", kg.graph, "graph_name", GMLFormat())
    cycles = simplecycles(kg.graph)
    return [kg.idtoname[cycle] for cycle in cycles]
end

function generate_dot_string(graph::SimpleDiGraph, idtoname::Array{String}, name::String)
    dot = ["DiGraph $(name) {"]
    for v in vertices(graph)
        push!(dot, "$(v) [label=\"$(idtoname[v])\"];")
    end
    for e in edges(graph)
        push!(dot, "$(e.src) -> $(e.dst);")
    end
    push!(dot, "}")
    return join(dot, "\n")
end

write_string_to_file(str::String, path::String) =  open(path, "w") do io
    write(io, str)
end

function export_knowledge_graph(kg::KnowledgeGraph, path::String)
    dot = generate_dot_string(kg.graph, kg.idtoname, "KnowledgeGraph")
    write_string_to_file(dot, path)
end

function export_knowledge_graph_towards_target(kg::KnowledgeGraph, target::String, path::String)
    ntarget = normalize_concept_name(target)
    subgraph_vertices = dfs_preordering(kg.graph, kg.nametoid[ntarget])
    # println(subgraph_vertices)
    subgraph, vmap = induced_subgraph(kg.graph, subgraph_vertices)
    dot = generate_dot_string(subgraph, kg.idtoname[vmap], filter(c-> (c>='a' && c<='z'), ntarget))
    write_string_to_file(dot, path)
end

"""
please make sure there are no loops in the graph before calling this function.
return value: postordering ordered vector of vertices traversed by by dfs

using the non-recursive dfs pseudocode from https://en.wikipedia.org/wiki/Depth-first_search#Pseudocode

procedure DFS-iterative(G, v) is
    let S be a stack
    S.push(v)
    while S is not empty do
        v = S.pop()
        if v is not labeled as discovered then
            label v as discovered
            for all edges from v to w in G.adjacentEdges(v) do
                S.push(w)

this method is not recursive, because having the knowledge graph size limited by julia call stack size limit is not ok.

this method will use a function to sort a vertex's neighbors before traversing them to guarantee deterministic ordering. notice that neighbors_sort_fn needs to output the reverse of your intended traverse order because we are using a stack to hold the sorted vertex list. you can supply nothing to neighbors_sort_fn to slightly speed up the program.
"""
function dfs_postordering(graph::SimpleDiGraph, start::Int; neighbors_sort_fn = reverse ∘ sort)
    S = [start]
    n = nv(graph)
    recursive_postorder_callstack = Array{Tuple{Int, Int}}(undef, 0) # dfs post order is much easier to understand in recursive traversal. so we try to picture what would happen in recursive traversal to compute post order
    postorder = Array{Int}(undef, 0)
    # preorder = Array{Int}(undef, 0)
    discovered = BitArray(undef, n) .& 0 # a hash optimized for worst case senario

    while length(S) > 0
        v = pop!(S)
        if discovered[v] == 0
            discovered[v] = 1
            # push!(preorder, v)
            v_neighbors = neighbors(graph, v)
            if nothing != neighbors_sort_fn
                v_neighbors = neighbors_sort_fn(v_neighbors)
            end
            # println("v $(v)'s sorted neighbors: $(v_neighbors)")
            v_neighbors_discovered = 0
            v_first_undiscovered_neighbor = 0
            for w in v_neighbors
                if discovered[w] == 0
                    # only push undiscovered neighbors
                    push!(S, w)
                    if 0 == v_first_undiscovered_neighbor
                        v_first_undiscovered_neighbor = w
                    end
                else
                    v_neighbors_discovered += 1
                end
            end

            if v_neighbors_discovered == length(v_neighbors)
                # all of v's neighbors have been discovered. this is the last time we'll visit v
                push!(postorder, v)
                # now check if v is it's parent's first undiscovered neighbor
                while length(recursive_postorder_callstack) > 0 && recursive_postorder_callstack[end][2] == v
                    v_parent = recursive_postorder_callstack[end][2]
                    #all of v's parent's neighbors have been pushed to postorder.
                    while length(recursive_postorder_callstack) > 0 && recursive_postorder_callstack[end][2] == v
                        frame = pop!(recursive_postorder_callstack)
                        # println("poping $(frame) from recursive_postorder_callstack for v: $(v)")
                        v_parent = frame[1]
                        push!(postorder, v_parent)
                    end
                    v = v_parent #important! let's move on to check v's parent's parent
                end
            else
                # some of v's neighbors have not been discovered. they are pushed to S.
                # println("pushing to recursive_postorder_callstack: $((v, v_first_undiscovered_neighbor))")
                push!(recursive_postorder_callstack, (v, v_first_undiscovered_neighbor))
            end
        end
    end
    # println("preorder: $(preorder)")
    # println("postorder: $(postorder)")
    # println("recursive_postorder_callstack: $(recursive_postorder_callstack)")
    return postorder
end

"""
using the non-recursive dfs pseudocode from https://en.wikipedia.org/wiki/Depth-first_search#Pseudocode

procedure DFS-iterative(G, v) is
    let S be a stack
    S.push(v)
    while S is not empty do
        v = S.pop()
        if v is not labeled as discovered then
            label v as discovered
            for all edges from v to w in G.adjacentEdges(v) do
                S.push(w)

why dfs_preordering when i already coded dfs_postordering? because dfs_preordering is more efficient when computing the collection of vertexes in the start vertexes' subgraph.
"""
function dfs_preordering(graph::SimpleDiGraph, start::Int)
    S = [start]
    n = nv(graph)
    preorder = Array{Int}(undef, 0)
    discovered = BitArray(undef, n) .& 0 # a hash optimized for worst case senario

    while length(S) > 0
        v = pop!(S)
        if discovered[v] == 0
            discovered[v] = 1
            push!(preorder, v)
            v_neighbors = neighbors(graph, v)
            for w in v_neighbors
                if discovered[w] == 0
                    # only push undiscovered neighbors
                    push!(S, w)
                end
            end
        end
    end
    # println("preorder: $(preorder)")
    return preorder
end

"""
given a target concept, compute a sequence to learn all of its prerequisite concepts. every concept in this sequence will have its own prerequisites placed before it.

given the above definition, the result sequence is a reverse of the topological sort (https://en.wikipedia.org/wiki/Topological_sorting) of the subgraph induced by target.

please make sure the knowledge graph contains no cycles before calling this function. in other words, please make sure the input knowledge graph is a dag.

since reverse postordering of dfs travsersal of dag produces a topological sorting (https://en.wikipedia.org/wiki/Depth-first_search#Vertex_orderings), we'll just use dfs postordering to produce the result sequence.

dfs postordering ensures that every vertex is placed adjacent to its descendants. i think this is good for maintaining relevenance concepts in working memory while learning.

however dfs postordering might not be unqiue for a dag. i would use a special sorting function to determine the order in which a vertex's children is visited in dfs.

this sorting function will prioritize the children that have smaller induced subgraph. the rationale is that we'll try to make the most progress towards the target concept as early as possible, by picking the lowest hanging fruit first.

dfs traversal can provide the size of induced subgraph for a vertex. however i haven't found a way to compute induced subgraph size for all vertexes in a graph efficiently. so for now i'm using brutal force approach, doing a dfs for every single vertex concerned.
"""
function generate_learning_plan(kg::KnowledgeGraph, target::String)
    ntarget = normalize_concept_name(target)
    if !haskey(kg.nametoid, ntarget)
        println("$(target) is not in knowledge graph.")
        return []
    end
    targetid = kg.nametoid[ntarget]
    subgraph_size_dict = Dict{Int, Int}()
    subgraph_vertexes = dfs_preordering(kg.graph, targetid)
    subgraph_size_dict[targetid] = length(subgraph_vertexes)
    deleteat!(subgraph_vertexes, 1)
    for vertex in subgraph_vertexes
        subgraph_size_dict[vertex] = length(dfs_preordering(kg.graph, vertex))
    end

    # println("subgraph size:")
    # for (k, v) in subgraph_size_dict
    #     println("  $(kg.idtoname[k]): $(v)")
    # end

    sort_fn = children -> sort(children, lt=(x, y)->subgraph_size_dict[x] > subgraph_size_dict[y])
    learning_order = dfs_postordering(kg.graph, targetid, neighbors_sort_fn=sort_fn)
    return learning_order

end

"""
perhaps we can also code a function to write all notes to a file, ordered by generate_learning_plan. that requires us to identify all the root concepts (concepts which no other concepts depends on) in the knowledge graph.
but i don't think i need that. i think it's good to always have a target concept in mind when reading notes.
"""
function write_notes_for_target_concept_to_md_file(kg::KnowledgeGraph, concepts::Array{Concept, 1}, target::String, path::String)
    ntarget = normalize_concept_name(target)
    lines = ["# notes for $(target)"]
    learning_plan = generate_learning_plan(kg, ntarget)
    idtouuid = Dict{Int, UUID}()
    missing_concept_set = get_missing_concepts(concepts)
    for id in learning_plan
        concept = concepts[id]
        uuid = uuid4()
        idtouuid[id] = uuid
        push!(lines, "## $(concept.name)")
        push!(lines, "### link\n  [$(concept.link)]($(concept.link))")

        if length(concept.dependency) > 0
            push!(lines, "\n### prerequisites")
            prerequisite_lines = []
            for prerequisite_concept in concept.dependency
                nprerequisite = normalize_concept_name(prerequisite_concept)
                if nprerequisite in missing_concept_set
                    push!(prerequisite_lines, "missing: $(prerequisite_concept)")
                else
                    prerequisite_uuid = idtouuid[kg.nametoid[nprerequisite]]
                    push!(prerequisite_lines, "[$(prerequisite_concept)](#$(prerequisite_uuid))")
                end
            end
            push!(lines, join(prerequisite_lines, " | \n"))
        end
        push!(lines, "\n### definition  ")
        push!(lines, concept.text)
        push!(lines, "\n  \n----  \n")
    end

    md_text = join(lines, "\n")
    write_string_to_file(md_text, path)
end

"""
utility function that checks for duplicate, missing and cyclical concepts
return false if a problem is found. otherwise return true
"""
function check_notes(concepts::Array{Concept, 1})
    dc = get_duplicate_concepts(concepts)
    if length(dc) > 0
        println("duplicate concepts:")
        for entry in dc
            println(entry)
        end
        return false
    end

    kg = init_knowledge_graph(concepts)
    cycles = find_cycles(kg)
    if length(cycles) > 0
        println("cycle(s) found")
        for c in cycles
            println(c)
        end
        return false
    end

    mc = get_missing_concepts(concepts)
    if length(mc) > 0
        println("missing concepts:")
        for entry in mc
            println(entry)
        end
        return false
    end

    return true
end

function fuzzy_search(concepts::Array{Concept, 1}, search_term)
    idx = findall(search_term, getfield.(concepts, :name),  Levenshtein())
    return concepts[idx]
end

export print_concepts, get_missing_concepts, get_duplicate_concepts, get_known_concepts, find_cycles, init_knowledge_graph,  export_knowledge_graph_towards_target, generate_learning_plan, write_notes_for_target_concept_to_md_file, export_knowledge_graph, Concept, check_notes, fuzzy_search
end # module

max_stars_repo_name: bluepine/KnowledgeGraphNote

max_issues_repo_name: bluepine/KnowledgeGraphNote

max_forks_repo_name: bluepine/KnowledgeGraphNote


text:
using Test
using ViewerGL
GL = ViewerGL
using StaticArrays

@testset "Box.jl" begin

   # 	function Box3d()
   @testset "Box3d()" begin
      @test typeof(GL.Box3d())==ViewerGL.Box3d
      @test GL.Box3d().p1==[0.0,0.0,0.0]
      @test GL.Box3d().p2==[0.0,0.0,0.0]
      @test typeof(GL.Box3d().p1)==StaticArrays.MArray{Tuple{3},Float64,1,3}
      @test typeof(GL.Box3d().p2)==StaticArrays.MArray{Tuple{3},Float64,1,3}
      # function Box3d(p1::Point3d,p2::Point3d)
      p1=GL.Point3d(-1,-1,-1); p2=GL.Point3d(+1,+1,+1);
      @test GL.Box3d(p1::GL.Point3d,p2::GL.Point3d)==ViewerGL.Box3d(GL.Point3d(-1.0, -1.0, -1.0), GL.Point3d(1.0, 1.0, 1.0))
      @test typeof(ViewerGL.Box3d)==DataType
   end

   # function invalidBox()
   @testset "invalidBox()" begin
      @test GL.invalidBox()==ViewerGL.Box3d(GL.Point3d(Inf,Inf,Inf),GL.Point3d(-Inf,-Inf,-Inf))
      @test GL.invalidBox().p1==GL.Point3d(Inf,Inf,Inf)
      @test GL.invalidBox().p2==GL.Point3d(-Inf,-Inf,-Inf)
      @test typeof(GL.invalidBox().p1)==StaticArrays.MArray{Tuple{3},Float64,1,3}
   end

   # function addPoint(box::Box3d,point::Point3d)
   @testset "addPoint" begin
      box=GL.invalidBox(); point=GL.Point3d(0,0,0)
      @test typeof(GL.addPoint(box::GL.Box3d,point::GL.Point3d))==ViewerGL.Box3d
      @test GL.addPoint(ViewerGL.Box3d(),GL.Point3d(0,0,1))==ViewerGL.Box3d(GL.Point3d(0.0,0.0,0.0),GL.Point3d(0.0,0.0,1.0))
      @test typeof(GL.addPoint(ViewerGL.Box3d(),GL.Point3d(0,0,1)))==ViewerGL.Box3d
   end

   # function getPoints(box::Box3d)
   a = GL.addPoint(ViewerGL.Box3d(),GL.Point3d(0,0,1))
   GL.addPoint(a,GL.Point3d(1,0,1))
   GL.addPoint(a,GL.Point3d(0,1,1))
   @testset "getPoints" begin
      @test a==ViewerGL.Box3d(GL.Point3d(0.0,0.0,0.0), GL.Point3d(1.0,1.0,1.0))
      box = GL.getPoints(GL.Box3d(GL.Point3d(0,0,0),GL.Point3d(1,1,1)))
      @test typeof(box)==Array{MArray{Tuple{3},Float64,1,3},1}
      @test box==convert(Array{MArray{Tuple{3},Float64,1,3},1},[[0.0, 0.0, 0.0],
                             [1.0, 0.0, 0.0],
                             [1.0, 1.0, 0.0],
                             [0.0, 1.0, 0.0],
                             [0.0, 0.0, 1.0],
                             [1.0, 0.0, 1.0],
                             [1.0, 1.0, 1.0],
                             [0.0, 1.0, 1.0]])
   end
end

ext: jl


text:
using Downloads: download
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the shallow water equations with a discontinuous
# bottom topography function

equations = ShallowWaterEquations2D(gravity_constant=9.81)

# Note, this initial condition is used to compute errors in the analysis callback but the initialization is
# overwritten by `initial_condition_ec_discontinuous_bottom` below.
initial_condition = initial_condition_weak_blast_wave

###############################################################################
# Get the DG approximation space

volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal)
solver = DGSEM(polydeg=6, surface_flux=(flux_fjordholm_etal, flux_nonconservative_fjordholm_etal),
               volume_integral=VolumeIntegralFluxDifferencing(volume_flux))

###############################################################################
# This setup is for the curved, split form entropy conservation testing (needs periodic BCs)

# Get the unstructured quad mesh from a file (downloads the file if not available locally)
default_mesh_file = joinpath(@__DIR__, "mesh_alfven_wave_with_twist_and_flip.mesh")
isfile(default_mesh_file) || download("https://gist.githubusercontent.com/andrewwinters5000/8f8cd23df27fcd494553f2a89f3c1ba4/raw/85e3c8d976bbe57ca3d559d653087b0889535295/mesh_alfven_wave_with_twist_and_flip.mesh",
                                      default_mesh_file)
mesh_file = default_mesh_file

mesh = UnstructuredMesh2D(mesh_file, periodicity=true)

# Create the semi discretization object
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

###############################################################################
# ODE solver

tspan = (0.0, 2.0)
ode = semidiscretize(semi, tspan)

###############################################################################
# Workaround to set a discontinuous bottom topography and initial condition for debugging and testing.

# alternative version of the initial conditinon used to setup a truly discontinuous
# bottom topography function and initial condtion for this academic testcase of entropy conservation.
# The errors from the analysis callback are not important but `∑∂S/∂U ⋅ Uₜ` should be around machine roundoff
# In contrast to the usual signature of initial conditions, this one get passed the
# `element_id` explicitly. In particular, this initial conditions works as intended
# only for the specific mesh loaded above!
function initial_condition_ec_discontinuous_bottom(x, t, element_id, equations::ShallowWaterEquations2D)
  # Set up polar coordinates
  inicenter = SVector(0.7, 0.7)
  x_norm = x[1] - inicenter[1]
  y_norm = x[2] - inicenter[2]
  r = sqrt(x_norm^2 + y_norm^2)
  phi = atan(y_norm, x_norm)
  sin_phi, cos_phi = sincos(phi)

  # Set the background values
  H = 3.25
  v1 = 0.0
  v2 = 0.0
  b = 0.0

  # setup the discontinuous water height and velocities
  if element_id == 10
    H = 4.0
    v1 = 0.1882 * cos_phi
    v2 = 0.1882 * sin_phi
  end

  # Setup a discontinuous bottom topography using the element id number
  if element_id == 7
    b = 2.0 + 0.5 * sin(2.0 * pi * x[1]) + 0.5 * cos(2.0 * pi * x[2])
  end

  return prim2cons(SVector(H, v1, v2, b), equations)
end

# point to the data we want to augment
u = Trixi.wrap_array(ode.u0, semi)
# reset the initial condition
for element in eachelement(semi.solver, semi.cache)
  for j in eachnode(semi.solver), i in eachnode(semi.solver)
    x_node = Trixi.get_node_coords(semi.cache.elements.node_coordinates, equations, semi.solver, i, j, element)
    u_node = initial_condition_ec_discontinuous_bottom(x_node, first(tspan), element, equations)
    Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, j, element)
  end
end

###############################################################################
# Callbacks

summary_callback = SummaryCallback()

analysis_interval = 1000
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_solution = SaveSolutionCallback(interval=1000,
                                     save_initial_solution=true,
                                     save_final_solution=true)

stepsize_callback = StepsizeCallback(cfl=3.0)

callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution,
                        stepsize_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct ESST4B <: AVR
        Tr::Float64
        K_pr::Float64
        K_ir::Float64
        Vr_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        Ta::Float64
        K_pm::Float64
        K_im::Float64
        Vm_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        Kg::Float64
        Kp::Float64
        Ki::Float64
        VB_max::Float64
        Kc::Float64
        Xl::Float64
        θp::Float64
        V_ref::Float64
        θp_rad::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int
        states_types::Vector{StateTypes.StateType}
        internal::InfrastructureSystemsInternal
    end

In these excitation systems, voltage (and also current in compounded systems) is transformed to an appropriate level. Rectifiers, either controlled or non-controlled, provide the necessary direct current for the generator field.
Parameters of IEEE Std 421.5 Type ST4B Excitacion System. ESST4B in PSSE and PSLF

# Arguments
- `Tr::Float64`: Regulator input filter time constant in s, validation range: `(0, 0.5)`, action if invalid: `warn`
- `K_pr::Float64`: Regulator propotional gain, validation range: `(0, 75)`, action if invalid: `warn`
- `K_ir::Float64`: Regulator integral gain, validation range: `(0, 75)`, action if invalid: `warn`
- `Vr_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`: Voltage regulator limits (Vi_min, Vi_max)
- `Ta::Float64`: Voltage regulator time constant in s, validation range: `(0, 1)`, action if invalid: `warn`
- `K_pm::Float64`: Voltage regulator proportional gain output, validation range: `(0, 1.2)`, action if invalid: `warn`
- `K_im::Float64`: Voltage regulator integral gain output, validation range: `(0, 18)`, action if invalid: `warn`
- `Vm_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`: Limits for inner loop output `(Vm_min, Vm_max)`
- `Kg::Float64`: Feedback gain constant of the inner loop field regulator, validation range: `(0, 1.1)`, action if invalid: `warn`
- `Kp::Float64`: Potential circuit (voltage) gain coefficient, validation range: `(0, 10)`, action if invalid: `warn`
- `Ki::Float64`: Compound circuit (current) gain coefficient, validation range: `(0, 1.1)`, action if invalid: `warn`
- `VB_max::Float64`: Maximum available exciter voltage, validation range: `(1, 20)`, action if invalid: `warn`
- `Kc::Float64`: Rectifier loading factor proportional to commutating reactance, validation range: `(0, 1)`, action if invalid: `warn`
- `Xl::Float64`: Reactance associated with potential source, validation range: `(0, 0.5)`, action if invalid: `warn`
- `θp::Float64`: Potential circuit phase angle (degrees), validation range: `(-90, 90)`, action if invalid: `error`
- `V_ref::Float64`: Reference Voltage Set-point, validation range: `(0, nothing)`
- `θp_rad::Float64`: Potential circuit phase angle (radians)
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`: The states are:
	Vm: Sensed terminal voltage,
	Vt: Sensed Terminal Voltage,
	Vr1: Regulator Integrator,
	Vr2: Regulator Output,
	Vm: Output integrator
- `n_states::Int`: ST4B has 4 states
- `states_types::Vector{StateTypes.StateType}`: ST4B has 4 states
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct ESST4B <: AVR
    "Regulator input filter time constant in s"
    Tr::Float64
    "Regulator propotional gain"
    K_pr::Float64
    "Regulator integral gain"
    K_ir::Float64
    "Voltage regulator limits (Vi_min, Vi_max)"
    Vr_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    "Voltage regulator time constant in s"
    Ta::Float64
    "Voltage regulator proportional gain output"
    K_pm::Float64
    "Voltage regulator integral gain output"
    K_im::Float64
    "Limits for inner loop output `(Vm_min, Vm_max)`"
    Vm_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    "Feedback gain constant of the inner loop field regulator"
    Kg::Float64
    "Potential circuit (voltage) gain coefficient"
    Kp::Float64
    "Compound circuit (current) gain coefficient"
    Ki::Float64
    "Maximum available exciter voltage"
    VB_max::Float64
    "Rectifier loading factor proportional to commutating reactance"
    Kc::Float64
    "Reactance associated with potential source"
    Xl::Float64
    "Potential circuit phase angle (degrees)"
    θp::Float64
    "Reference Voltage Set-point"
    V_ref::Float64
    "Potential circuit phase angle (radians)"
    θp_rad::Float64
    ext::Dict{String, Any}
    "The states are:
	Vm: Sensed terminal voltage,
	Vt: Sensed Terminal Voltage,
	Vr1: Regulator Integrator,
	Vr2: Regulator Output,
	Vm: Output integrator"
    states::Vector{Symbol}
    "ST4B has 4 states"
    n_states::Int
    "ST4B has 4 states"
    states_types::Vector{StateTypes.StateType}
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function ESST4B(Tr, K_pr, K_ir, Vr_lim, Ta, K_pm, K_im, Vm_lim, Kg, Kp, Ki, VB_max, Kc, Xl, θp, V_ref=1.0, θp_rad=θp*π*inv(180), ext=Dict{String, Any}(), )
    ESST4B(Tr, K_pr, K_ir, Vr_lim, Ta, K_pm, K_im, Vm_lim, Kg, Kp, Ki, VB_max, Kc, Xl, θp, V_ref, θp_rad, ext, [:Vt, :Vr1, :Vr2, :Vm], 4, [StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Hybrid], InfrastructureSystemsInternal(), )
end

function ESST4B(; Tr, K_pr, K_ir, Vr_lim, Ta, K_pm, K_im, Vm_lim, Kg, Kp, Ki, VB_max, Kc, Xl, θp, V_ref=1.0, θp_rad=θp*π*inv(180), ext=Dict{String, Any}(), states=[:Vt, :Vr1, :Vr2, :Vm], n_states=4, states_types=[StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Hybrid], internal=InfrastructureSystemsInternal(), )
    ESST4B(Tr, K_pr, K_ir, Vr_lim, Ta, K_pm, K_im, Vm_lim, Kg, Kp, Ki, VB_max, Kc, Xl, θp, V_ref, θp_rad, ext, states, n_states, states_types, internal, )
end

# Constructor for demo purposes; non-functional.
function ESST4B(::Nothing)
    ESST4B(;
        Tr=0,
        K_pr=0,
        K_ir=0,
        Vr_lim=(min=0.0, max=0.0),
        Ta=0,
        K_pm=0,
        K_im=0,
        Vm_lim=(min=0.0, max=0.0),
        Kg=0,
        Kp=0,
        Ki=0,
        VB_max=0,
        Kc=0,
        Xl=0,
        θp=0,
        V_ref=0,
        θp_rad=0,
        ext=Dict{String, Any}(),
    )
end

"""Get [`ESST4B`](@ref) `Tr`."""
get_Tr(value::ESST4B) = value.Tr
"""Get [`ESST4B`](@ref) `K_pr`."""
get_K_pr(value::ESST4B) = value.K_pr
"""Get [`ESST4B`](@ref) `K_ir`."""
get_K_ir(value::ESST4B) = value.K_ir
"""Get [`ESST4B`](@ref) `Vr_lim`."""
get_Vr_lim(value::ESST4B) = value.Vr_lim
"""Get [`ESST4B`](@ref) `Ta`."""
get_Ta(value::ESST4B) = value.Ta
"""Get [`ESST4B`](@ref) `K_pm`."""
get_K_pm(value::ESST4B) = value.K_pm
"""Get [`ESST4B`](@ref) `K_im`."""
get_K_im(value::ESST4B) = value.K_im
"""Get [`ESST4B`](@ref) `Vm_lim`."""
get_Vm_lim(value::ESST4B) = value.Vm_lim
"""Get [`ESST4B`](@ref) `Kg`."""
get_Kg(value::ESST4B) = value.Kg
"""Get [`ESST4B`](@ref) `Kp`."""
get_Kp(value::ESST4B) = value.Kp
"""Get [`ESST4B`](@ref) `Ki`."""
get_Ki(value::ESST4B) = value.Ki
"""Get [`ESST4B`](@ref) `VB_max`."""
get_VB_max(value::ESST4B) = value.VB_max
"""Get [`ESST4B`](@ref) `Kc`."""
get_Kc(value::ESST4B) = value.Kc
"""Get [`ESST4B`](@ref) `Xl`."""
get_Xl(value::ESST4B) = value.Xl
"""Get [`ESST4B`](@ref) `θp`."""
get_θp(value::ESST4B) = value.θp
"""Get [`ESST4B`](@ref) `V_ref`."""
get_V_ref(value::ESST4B) = value.V_ref
"""Get [`ESST4B`](@ref) `θp_rad`."""
get_θp_rad(value::ESST4B) = value.θp_rad
"""Get [`ESST4B`](@ref) `ext`."""
get_ext(value::ESST4B) = value.ext
"""Get [`ESST4B`](@ref) `states`."""
get_states(value::ESST4B) = value.states
"""Get [`ESST4B`](@ref) `n_states`."""
get_n_states(value::ESST4B) = value.n_states
"""Get [`ESST4B`](@ref) `states_types`."""
get_states_types(value::ESST4B) = value.states_types
"""Get [`ESST4B`](@ref) `internal`."""
get_internal(value::ESST4B) = value.internal

"""Set [`ESST4B`](@ref) `Tr`."""
set_Tr!(value::ESST4B, val) = value.Tr = val
"""Set [`ESST4B`](@ref) `K_pr`."""
set_K_pr!(value::ESST4B, val) = value.K_pr = val
"""Set [`ESST4B`](@ref) `K_ir`."""
set_K_ir!(value::ESST4B, val) = value.K_ir = val
"""Set [`ESST4B`](@ref) `Vr_lim`."""
set_Vr_lim!(value::ESST4B, val) = value.Vr_lim = val
"""Set [`ESST4B`](@ref) `Ta`."""
set_Ta!(value::ESST4B, val) = value.Ta = val
"""Set [`ESST4B`](@ref) `K_pm`."""
set_K_pm!(value::ESST4B, val) = value.K_pm = val
"""Set [`ESST4B`](@ref) `K_im`."""
set_K_im!(value::ESST4B, val) = value.K_im = val
"""Set [`ESST4B`](@ref) `Vm_lim`."""
set_Vm_lim!(value::ESST4B, val) = value.Vm_lim = val
"""Set [`ESST4B`](@ref) `Kg`."""
set_Kg!(value::ESST4B, val) = value.Kg = val
"""Set [`ESST4B`](@ref) `Kp`."""
set_Kp!(value::ESST4B, val) = value.Kp = val
"""Set [`ESST4B`](@ref) `Ki`."""
set_Ki!(value::ESST4B, val) = value.Ki = val
"""Set [`ESST4B`](@ref) `VB_max`."""
set_VB_max!(value::ESST4B, val) = value.VB_max = val
"""Set [`ESST4B`](@ref) `Kc`."""
set_Kc!(value::ESST4B, val) = value.Kc = val
"""Set [`ESST4B`](@ref) `Xl`."""
set_Xl!(value::ESST4B, val) = value.Xl = val
"""Set [`ESST4B`](@ref) `θp`."""
set_θp!(value::ESST4B, val) = value.θp = val
"""Set [`ESST4B`](@ref) `V_ref`."""
set_V_ref!(value::ESST4B, val) = value.V_ref = val
"""Set [`ESST4B`](@ref) `θp_rad`."""
set_θp_rad!(value::ESST4B, val) = value.θp_rad = val
"""Set [`ESST4B`](@ref) `ext`."""
set_ext!(value::ESST4B, val) = value.ext = val
"""Set [`ESST4B`](@ref) `states_types`."""
set_states_types!(value::ESST4B, val) = value.states_types = val


text:
function instrument_rectangular(x, res)
    y = zeros(eltype(x), length(x))
    y[-res / 2 .<= x .< res / 2] .= 1 / res
    return y
end

function instrument_triangular(x, res)
    y = zeros(eltype(x), length(x))
    y[abs.(x) .<= res] .= 1 / res * ( 1 .- abs.(@view(x[abs.(x) .<= res])) ./ res)
    return y
end

instrument_gaussian(x, res) = @. 2*√(log(2)/π)/res * exp(-log(2)*(2*x/res)^2)

instrument_lorentzian(x, res) = @. res / 2. / π / (x^2 + (res / 2.)^2 )

function instrument_cosine(x, res)
    y = zeros(eltype(x), length(x))
    y[abs.(x) .<= res] .= @. cos(π*@view(x[abs.(x) .<= res])/(2*res))/res*pi/4
    return y
end

instrument_diffraction(x, res) = @. 1/res * sinc(x/res)^2

instrument_michelson(x, res) = @. 2/res * sinc(2/res*x)

const instrument_functions = Dict{Symbol, Function}(
    :rectangular => instrument_rectangular,
    :triangular => instrument_triangular,
    :gaussian => instrument_gaussian,
    :lorentzian => instrument_lorentzian,
    :cosine => instrument_cosine,
    :diffraction => instrument_diffraction,
    :michelson => instrument_michelson
)

"""
    apply_instrument_function(ν, α[, instrument_function=:rectangular, instrument_wing=10.0, instrument_resolution=0.1])

Applies an instrument function to the given input spectrum.

# Arguments
- `ν`: The wavenumber vector (NOTE: Uniform spacing is assumed but not checked!)
- `α`: The calculated absorption coefficient using [`α`](@ref)
- `instrument_function` (optional): A Symbol describing one of the instrument functions below
- `instrument_wing` (optional): The half-width of the range for calculating the instrument function in ``cm^{-1}``
- `instrument_resolution` (optional): The full-width of the instrument resolution in ``cm^{-1}``

# Output

Returns a new vector with the spectrum influenced by the instrument function

# Instrument functions

The following instrument functions ``I(x, Δ) `` are supported. Here `x`is the coordinate for evaluating the function, whose range is given by `instrument_wing`. `Δ` is the resolution parameter `instrument_resolution`. Use the stated symbol as value
for the argument `instrument_function`.

| Symbol | Equation |  Description |
| :---   |      :---: | ---: |
| `:rectangular` | `` \\begin{cases} \\frac{1}{Δ} & \\lvert x \\rvert \\leq \\frac{Δ}{2} \\\\ 0 & \\lvert x \\rvert > \\frac{Δ}{2} \\end{cases} `` |  A rectangular instrument function |
| `:triangular` | `` \\begin{cases} \\frac{1}{Δ} (1 - \\frac{\\lvert x \\rvert}{Δ}) & \\lvert x \\rvert \\leq Δ \\\\ 0 & \\lvert x \\rvert > Δ \\end{cases} ``  | A triangular instrument function |
| `:gaussian` | `` \\frac{2}{Δ} \\sqrt{\\frac{\\mathrm{ln}2}{\\pi}} \\mathrm{exp} \\left (- \\mathrm{ln}2 \\left ( \\frac{2x}{Δ}\\right)^2 \\right )`` | A Gaussian instrument function (e.g. a broadband source) |
| `:lorentzian` | `` \\frac{Δ}{2\\pi} \\frac{1}{x^2+\\left(\\frac{Δ}{2}\\right)^2} `` | A Lorentzian instrument function (e.g. a single frequency laser) |
| `:cosine` | `` \\begin{cases} \\frac{\\pi}{4Δ} \\cos \\left ( \\frac {\\pi \\lvert x \\rvert}{2Δ} \\right ) & \\lvert x \\rvert \\leq Δ \\\\ 0 & \\lvert x \\rvert > Δ \\end{cases} `` | A cosine instrument function |
| `:diffraction` | `` \\frac{1}{Δ} \\mathrm{sinc}^2 \\left(  \\frac{\\pi x}{Δ} \\right) `` | A diffraction (sinc-type) instrument function |
| `:michelson` | `` \\frac{2}{Δ} \\mathrm{sinc} \\left(  \\frac{2 \\pi x}{Δ} \\right) `` | A Michelson interferometer-type instrument function (e.g. FTIR) |

"""
function apply_instrument_function(
    ν::AbstractVector{T},
    α::AbstractVector{T},
    instrument_function::Symbol=:rectangular,
    instrument_wing::AbstractFloat=10.0,
    instrument_resolution::AbstractFloat=0.1
) where T <: AbstractFloat
    δν = ν[2] - ν[1]
    x = -instrument_wing:δν:instrument_wing
    ifn = instrument_functions[instrument_function](x, instrument_resolution)
    ifn ./= sum(ifn)
    α_inst = conv(α, ifn)
    ind_lo = length(ifn) ÷ 2 + 1
    ind_hi = length(α_inst) - ind_lo + 2
    return α_inst[ind_lo:ind_hi]
end


text:
import Base: size, show

"""
    RegularCartesianGrid{T<:AbstractFloat, R<:AbstractRange} <: Grid

A Cartesian grid with regularly spaces cells and faces so that \$Δx\$, \$Δy\$,
and \$Δz\$ are constants. Fields are stored using floating-point values of type
`T`. The cell-centered and face-centered coordinate ranges are stored as ranges
of type `A`.
"""
struct RegularCartesianGrid{T<:AbstractFloat, R<:AbstractRange} <: Grid
    dim::Int
    # Number of grid points in (x,y,z).
    Nx::Int
    Ny::Int
    Nz::Int
    # Domain size [m].
    Lx::T
    Ly::T
    Lz::T
    # Grid spacing [m].
    Δx::T
    Δy::T
    Δz::T
    # Cell face areas [m²].
    Ax::T
    Ay::T
    Az::T
    # Volume of a cell [m³].
    V::T
    # Range of coordinates at the centers of the cells.
    xC::R
    yC::R
    zC::R
    # Range of grid coordinates at the faces of the cells.
    # Note: there are Nx+1 faces in the x-dimension, Ny+1 in the y, and Nz+1 in the z.
    xF::R
    yF::R
    zF::R
end

"""
    RegularCartesianGrid(T, N, L)

Create a regular Cartesian grid with size \$N = (N_x, N_y, N_z)\$ and domain
size \$L = (L_x, L_y, L_z)\$ where fields are stored using floating-point values
of type `T`.

# Examples
```julia-repl
julia> g = RegularCartesianGrid((16, 16, 8), (2π, 2π, 2π))
```
"""
function RegularCartesianGrid(T, N, L)
    !(length(N) == 3) && throw(ArgumentError("N=$N must be a tuple of length 3."))
    !(length(L) == 3) && throw(ArgumentError("L=$L must be a tuple of length 3."))

    !all(isa.(N, Integer)) && throw(ArgumentError("N=$N should contain integers."))
    !all(isa.(L, Number))  && throw(ArgumentError("L=$L should contain numbers."))

    !all(N .>= 1) && throw(ArgumentError("N=$N must be nonzero and positive!"))
    !all(L .> 0)  && throw(ArgumentError("L=$L must be nonzero and positive!"))

    !(T in [Float32, Float64]) && throw(ArgumentError("T=$T but only Float32 and Float64 grids are supported."))

    # Count the number of dimensions with 1 grid point, i.e. the number of flat
    # dimensions, and use it to determine the dimension of the model.
    num_flat_dims = count(i->(i==1), N)
    dim = 3 - num_flat_dims
    !(1 <= dim <= 3) && throw(ArgumentError("N=$N has dimension $dim. Only 1D, 2D, and 3D grids are supported."))

    Nx, Ny, Nz = N
    Lx, Ly, Lz = L

    dim == 2 && !(Nx == 1 || Ny == 1 || Nz == 1) &&
        throw(ArgumentError("For 2D grids, Nx, Ny, or Nz must be 1."))

    dim == 3 && !(Nx != 1 && Ny != 1 && Nz != 1) &&
        throw(ArgumentError("For 3D grids, cannot have dimensions of size 1."))

    Lx = convert(T, Lx)
    Ly = convert(T, Ly)
    Lz = convert(T, Lz)

    Δx = convert(T, Lx / Nx)
    Δy = convert(T, Ly / Ny)
    Δz = convert(T, Lz / Nz)

    Ax = convert(T, Δy*Δz)
    Ay = convert(T, Δx*Δz)
    Az = convert(T, Δx*Δy)

    V = convert(T, Δx*Δy*Δz)

    xC = Δx/2:Δx:Lx
    yC = Δy/2:Δy:Ly
    zC = -Δz/2:-Δz:-Lz

    xF = 0:Δx:Lx
    yF = 0:Δy:Ly
    zF = 0:-Δz:-Lz

    # Make sure all the coordinate ranges have the same type.
    !all(typeof.([xC, yC, zC, xF, yF, zF]) .== typeof(xC)) &&
        throw(ArgumentError("At least one coordinate range type did not match."))

    RegularCartesianGrid{T, typeof(xC)}(dim, Nx, Ny, Nz, Lx, Ly, Lz, Δx, Δy, Δz, Ax, Ay,
                                        Az, V, xC, yC, zC, xF, yF, zF)
end

# Constructor aliases.
RegularCartesianGrid(N, L) = RegularCartesianGrid(Float64, N, L)

size(g::RegularCartesianGrid) = (g.Nx, g.Ny, g.Nz)
Base.eltype(g::RegularCartesianGrid{T}) where T = T

show(io::IO, g::RegularCartesianGrid) =
    print(io, "$(g.dim)-dimensional ($(typeof(g.Lx))) regular Cartesian grid\n",
              "(Nx, Ny, Nz) = ", (g.Nx, g.Ny, g.Nz), '\n',
              "(Lx, Ly, Lz) = ", (g.Lx, g.Ly, g.Lz), '\n',
              "(Δx, Δy, Δz) = ", (g.Δx, g.Δy, g.Δz))


text:
using Pollen
using Pkg
using Crayons
using ImageShow
Crayons.COLORS[:nothing] = 67
ENV["DATADEPS_ALWAYS_ACCEPT"] = "true"

using FastAI, Flux, FluxTraining
import DataAugmentation
m = FastAI
ms = [
    DataAugmentation,
    Flux,
    FluxTraining,
    m,
]

project = Project(
    Pollen.Rewriter[
        DocumentFolder(Pkg.pkgdir(m), prefix = "documents"),
        ParseCode(),
        ExecuteCode(),
        PackageDocumentation(ms),
        StaticResources(),
        DocumentGraph(),
        SearchIndex(),
        SaveAttributes((:title,), useoutputs=true),
        LoadFrontendConfig(Pkg.pkgdir(m))
    ],
)


text:
### MH holds the fields of the Metropolis-Hastings sampler
### In its most general case it accommodates an asymmetric proposal density
### For symetric proposals, the proposal correction factor equals 1, so the logproposal field is set to nothing

immutable MH <: MHSampler
  symmetric::Bool # If the proposal density is symmetric, then symmetric=true, otherwise symmetric=false
  logproposal::Union(Function, Nothing) # logpdf of asymmetric proposal density
                                        # For symmetric proposals, logproposal is set to nothing
  randproposal::Function # random sampling from proposal density

  function MH(s::Bool, l::Union(Function, Nothing), r::Function)
    if s && !isa(l, Nothing)
      error("If the symmetric field is true, then logproposal is not used in the calculations.")
    end
    new(s, l, r)
  end
end

MH(l::Function, r::Function) = MH(false, l, r) # Metropolis-Hastings sampler (asymmetric proposal)
MH(r::Function) = MH(true, nothing, r) # Metropolis sampler (symmetric proposal)

# Random-walk Metropolis, i.e. Metropolis with a normal proposal density
MH(σ::Vector{Float64}) = MH(x::Vector{Float64} -> rand(MvNormal(x, σ)))
MH(σ::Float64) = MH(x::Vector{Float64} -> rand(MvNormal(x, σ)))
MH() = MH(x::Vector{Float64} -> rand(MvNormal(x, 1.0)))

### MHHeap type holds the internal state ("local variables") of the MH sampler

type MHHeap <: MCHeap{MCBaseSample}
  instate::MCState{MCBaseSample} # Monte Carlo state used internally by the sampler
  outstate::MCState{MCBaseSample} # Monte Carlo state outputted by the sampler
  tune::MCTune
  count::Int # Current number of iterations
  ratio::Float64
end

MHHeap() =
  MHHeap(MCState(MCBaseSample(), MCBaseSample()), MCState(MCBaseSample(), MCBaseSample()), VanillaMCTune(), 0, NaN)

MHHeap(l::Int, t::MCTune=VanillaMCTune()) =
  MHHeap(MCState(MCBaseSample(l), MCBaseSample(l)), MCState(MCBaseSample(l), MCBaseSample(l)), t, 0, NaN)

### Initialize Metropolis-Hastings sampler

function initialize_heap(m::MCModel, s::MH, r::MCRunner, t::MCTuner)
  heap::MHHeap = MHHeap(m.size)

  heap.instate.current = MCBaseSample(copy(m.init))
  logtarget!(heap.instate.current, m.eval)
  @assert isfinite(heap.instate.current.logtarget) "Initial values out of model support."

  heap.tune = VanillaMCTune()

  heap.count = 1

  heap
end

function reset!(heap::MHHeap, x::Vector{Float64}, m::MCModel)
  heap.instate.current = MCBaseSample(copy(x))
  logtarget!(heap.instate.current, m.eval)
end

function initialize_task!(heap::MHHeap, m::MCModel, s::MH, r::MCRunner, t::MCTuner)
  # Hook inside Task to allow remote resetting
  task_local_storage(:reset, (x::Vector{Float64})->reset!(heap, x, m))

  while true
    iterate!(heap, m, s, r, t, produce)
  end
end

### Perform iteration for Metropolis-Hastings sampler

function iterate!(heap::MHHeap, m::MCModel, s::MH, r::MCRunner, t::MCTuner, send::Function)
  if t.verbose
    heap.tune.proposed += 1
  end

  heap.instate.successive = MCBaseSample(s.randproposal(heap.instate.current.sample))
  logtarget!(heap.instate.successive, m.eval)

  if s.symmetric
    heap.ratio = heap.instate.successive.logtarget-heap.instate.current.logtarget
  else
    heap.ratio = (heap.instate.successive.logtarget
      +s.logproposal(heap.instate.successive.sample, heap.instate.current.sample)
      -heap.instate.current.logtarget
      -s.logproposal(heap.instate.current.sample, heap.instate.successive.sample)
    )
  end
  if heap.ratio > 0 || (heap.ratio > log(rand()))
    heap.outstate = MCState(heap.instate.successive, heap.instate.current, Dict{Any, Any}("accept" => true))
    heap.instate.current = deepcopy(heap.instate.successive)

    if t.verbose
      heap.tune.accepted += 1
    end
  else
    heap.outstate = MCState(heap.instate.current, heap.instate.current, Dict{Any, Any}("accept" => false))
  end

  if t.verbose && heap.count <= r.burnin && mod(heap.count, t.period) == 0
    rate!(heap.tune)
    println("Burnin iteration $(heap.count) of $(r.burnin): ", round(100*heap.tune.rate, 2), " % acceptance rate")
  end

  heap.count += 1

  send(heap.outstate)
end


text:
function lintmacro(ex::Expr, ctx::LintContext)
    @checktoplevel(ctx, "macro")

    if !isexpr(ex.args[1], :call)
        msg(ctx, :W100, ex.args[1], "this macro syntax not understood by Lint.jl")
        return
    end

    fname = ex.args[1].args[1]
    if !isa(fname, Symbol)
        msg(ctx, :E141, ex.args[1].args[1], "invalid macro syntax")
        return
    end

    # now we pretend we're a function...
    # TODO: make this less hackish
    fname = Symbol('@', fname)
    lintfunction(Expr(:function, Expr(:call, fname, ex.args[1].args[2:end]...),
                      ex.args[2:end]...), ctx)
end

istopmacro(ex, mod, mac) = ex in (
    mac,
    GlobalRef(mod, mac),
    Expr(:(.), Symbol(string(mod)), mac))

function lintcompat(ex::Expr, ctx::LintContext)
    if VERSION < v"0.6.0-dev.2746" &&
       length(ex.args) == 2 && isexpr(ex.args[2], :abstract) &&
       length(ex.args[2].args) == 1 && isexpr(ex.args[2].args[1], :type)
        lintexpr(Compat._compat_abstract(ex.args[2].args[1]), ctx)
    elseif VERSION < v"0.6.0-dev.2746" &&
           length(ex.args) == 3 && ex.args[2] == :primitive
        lintexpr(Compat._compat_primitive(ex.args[3]), ctx)
    elseif length(ex.args) == 2
        lintexpr(Compat._compat(ex.args[2]), ctx)
    else
        msg(ctx, :E437, ex, "@compat called with wrong number of arguments")
    end
end

function lintmacrocall(ex::Expr, ctx::LintContext)
    if istopmacro(ex.args[1], Base, Symbol("@deprecate")) ||
       ex.args[1] == Symbol("@recipe")
        return
    end

    if ex.args[1] == Symbol("@lintpragma")
        lintlintpragma(ex, ctx)
        return
    end

    if istopmacro(ex.args[1], Core, Symbol("@doc")) && length(ex.args) >= 2 # see Docile.jl
        if isexpr(ex.args[2], :(->))
            lintexpr(ex.args[2].args[2], ctx) # no need to lint the doc string
            return
        elseif typeof(ex.args[2]) <: AbstractString && length(ex.args) >= 3 && isexpr(ex.args[3], :call)
            # grandfather as a docstring of a previously declared function
            return
        elseif (typeof(ex.args[2]) <: AbstractString ||
              isexpr(ex.args[2], :macrocall) && ex.args[2].args[1] == Symbol("@mstr")
             )
            if length(ex.args) >= 3
                lintexpr(ex.args[3], ctx)
            else
                msg(ctx, :W443, "did you forget an -> after @doc or make it inline?")
            end
            return
        end
        return
    end

    if ex.args[1] == Symbol("@pyimport")
        if length(ex.args) == 2 && typeof(ex.args[2]) == Symbol
            set!(ctx.current, ex.args[2], VarInfo(location(ctx)))
        elseif length(ex.args) == 4 && ex.args[3] == :as && typeof(ex.args[4]) == Symbol
            set!(ctx.current, ex.args[4], VarInfo(location(ctx)))
        end
        return
    end

    if istopmacro(ex.args[1], Compat, Symbol("@compat"))
        lintcompat(ex, ctx)
        return
    end

    if ex.args[1] == Symbol("@gensym")
        for i in 2:length(ex.args)
            if typeof(ex.args[i]) == Symbol
                set!(ctx.current, ex.args[i], VarInfo(location(ctx)))
            end
        end
        return
    end

    if ex.args[1] == Symbol("@enum")
        @checktoplevel(ctx, "@enum")
        for i in 2:length(ex.args)
            if typeof(ex.args[i]) == Symbol
                vi = VarInfo(location(ctx))
                set!(ctx.current, ex.args[i], vi)
            elseif isexpr(ex.args[i], :(=)) && typeof(ex.args[i].args[1]) == Symbol
                vi = VarInfo(location(ctx))
                set!(ctx.current, ex.args[i].args[1], vi)
            end
        end
        return
    end

    ctx.macrocallLvl = ctx.macrocallLvl + 1
    for i = 2:length(ex.args)
        lintexpr(ex.args[i], ctx)
    end
    ctx.macrocallLvl = ctx.macrocallLvl - 1
end


text:
using MatrixMath

const MeshDataPlaneQuad = rotr90flipdim2(Float32[
	-1 -1 0 0	# upleft
	1 -1 0 0	# upright
	1 1 0 0	# downright
	-1 1 0 0	# downleft
])

const MeshDataPlaneVertex = rotr90flipdim2(Float32[
	-1.0 -1.0 0.0 1.0 0.0
	1.0 -1.0 0.0 0.0 0.0
	1.0 1.0 0.0 0.0 1.0
	1.0 1.0 0.0 0.0 1.0
	-1.0 1.0 0.0 1.0 1.0
	-1.0 -1.0 0.0 1.0 1.0

	-1.0 -1.0 0.0 1.0 0.0
	1.0 -1.0 0.0 0.0 0.0
	-1.0 1.0 0.0 1.0 1.0
	1.0 -1.0 0.0 0.0 0.0
	1.0 1.0 0.0 0.0 1.0
	-1.0 1.0 0.0 1.0 1.0
])

const MeshDataCube = Float32[
	#  X, Y, Z, U, V
	# Bottom
	-1.0, -1.0, -1.0, 0.0, 0.0,
	1.0, -1.0, -1.0, 1.0, 0.0,
	-1.0, -1.0, 1.0, 0.0, 1.0,
	1.0, -1.0, -1.0, 1.0, 0.0,
	1.0, -1.0, 1.0, 1.0, 1.0,
	-1.0, -1.0, 1.0, 0.0, 1.0,

	# Top
	-1.0, 1.0, -1.0, 0.0, 0.0,
	-1.0, 1.0, 1.0, 0.0, 1.0,
	1.0, 1.0, -1.0, 1.0, 0.0,
	1.0, 1.0, -1.0, 1.0, 0.0,
	-1.0, 1.0, 1.0, 0.0, 1.0,
	1.0, 1.0, 1.0, 1.0, 1.0,

	# Front
	-1.0, -1.0, 1.0, 1.0, 0.0,
	1.0, -1.0, 1.0, 0.0, 0.0,
	-1.0, 1.0, 1.0, 1.0, 1.0,
	1.0, -1.0, 1.0, 0.0, 0.0,
	1.0, 1.0, 1.0, 0.0, 1.0,
	-1.0, 1.0, 1.0, 1.0, 1.0,

	# Back
	-1.0, -1.0, -1.0, 0.0, 0.0,
	-1.0, 1.0, -1.0, 0.0, 1.0,
	1.0, -1.0, -1.0, 1.0, 0.0,
	1.0, -1.0, -1.0, 1.0, 0.0,
	-1.0, 1.0, -1.0, 0.0, 1.0,
	1.0, 1.0, -1.0, 1.0, 1.0,

	# Left
	-1.0, -1.0, 1.0, 0.0, 1.0,
	-1.0, 1.0, -1.0, 1.0, 0.0,
	-1.0, -1.0, -1.0, 0.0, 0.0,
	-1.0, -1.0, 1.0, 0.0, 1.0,
	-1.0, 1.0, 1.0, 1.0, 1.0,
	-1.0, 1.0, -1.0, 1.0, 0.0,

	# Right
	1.0, -1.0, 1.0, 1.0, 1.0,
	1.0, -1.0, -1.0, 1.0, 0.0,
	1.0, 1.0, -1.0, 0.0, 0.0,
	1.0, -1.0, 1.0, 1.0, 1.0,
	1.0, 1.0, -1.0, 0.0, 0.0,
	1.0, 1.0, 1.0, 0.0, 1.0,
]

function createPlaneData(wdetail::Integer, hdetail::Integer)
	texture_scale = true

	W = wdetail < 1 ? 1 : wdetail
	H = hdetail < 1 ? 1 : hdetail

	MW = W + 1
	MH = H + 1

	dW = 1.f0 / W
	dH = 1.f0 / H

	s = false

	vertices=Vec3f[]
	normals=Vec3f[]
	uvs=Vec2f[]
	indicies=UInt32[]

	for h=0:H
		hs = H - h - 1
		V = h * dH

		for w=0:W
			ws = W - w - 1
			U = w * dW

			v = Vec3f(-1 + U * 2, -1 + V * 2, 0) # vec3(-cos(theta) * sin(phi),-cos(phi),0);
			uv = texture_scale ? Vec2f(V, U) : Vec2f(h, w)

			push!(vertices, v)
			push!(normals, Vec3f(0,0,1))
			push!(uvs, MatrixMath.rotate(uv - Vec2f(1,0), 3.14f0*-0.5f0))

			if h < H && w < W
				if !s
					push!(indicies, w + (hs + 1) * MW) # top left
					push!(indicies, w + hs * MW) # left
					push!(indicies, w + 1 + (hs + 1) * MW) # top right
					push!(indicies, w + 1 + hs * MW) # right
				else
					push!(indicies, ws + 1 + (hs + 1) * MW) # top right
					push!(indicies, ws + 1 + hs * MW) # right
					push!(indicies, ws + (hs + 1) * MW) # top left
					push!(indicies, ws + hs * MW) # left
				end

				#push!(indicies, w + h * MW); # left
				#push!(indicies, w + 1 + h * MW) # right
				#push!(indicies, w + (h + 1) * MW) # top left
				#push!(indicies, w + 1 + (h + 1) * MW) # top right
				#push!(indicies, w + (h + 1) * MW) # top left
				#push!(indicies, w + 1 + h * MW) # right
			end
		end

		if h < H
			if !s push!(indicies, W + hs * MW) #top right
			else push!(indicies, hs * MW) # top left
			end
			s = !s
		end
	end

	(:TRIANGLE_STRIP, vertices, normals, uvs, indicies)
end

function createCubeDataSimple()
	const vd = Vec3f[
		Vec3f(-1, -1, 0),
		Vec3f(1, -1, 0),
		Vec3f(1, 1, 0),
		Vec3f(1, 1, 0),
		Vec3f(-1, 1, 0),
		Vec3f(-1, -1, 0),

		Vec3f(-1, -1, 0),
		Vec3f(1, -1, 0),
		Vec3f(-1, 1, 0),
		Vec3f(1, -1, 0),
		Vec3f(1, 1, 0),
		Vec3f(-1, 1, 0),
	]

	len = length(vd)
	vertices=zeros(Vec3f, len)
	normals=zeros(Vec3f, 0)
	uvs=zeros(Vec2f, 0)
	indicies=zeros(UInt32, 0)

	i=0
	for v in vd; i+=1; vertices[i]=v end
	#for v in vd; i+=1; vertices[i]=v end

	(:TRIANGLES, vertices, normals, uvs, indicies)
end

function createCubeData()
		v = [
			Vec3f(-1, -1, -1), Vec3f(1, -1, -1), Vec3f(-1, 1, -1), Vec3f(1, 1, -1),
			Vec3f(-1, -1, 1), Vec3f(1, -1, 1), Vec3f(-1, 1, 1), Vec3f(1, 1, 1)
		]

		n = zeros(Vec3f, 8)
		for ni=1:length(n) n[ni] = normalize(v[ni]) end

		uv = [ Vec2f(0, 1), Vec2f(1, 1), Vec2f(0, 0), Vec2f(1, 0)	]

		const i = UInt32[
			0, 1, 2, 3, # back (0, 1, 2, 3)
			4, 5, 6, 7, # front (4, 5, 6, 7)
			0, 4, 2, 6, # left (8, 9, 10, 11)
			1, 5, 3, 7, # right (12, 13, 14, 15)
			0, 1, 4, 5, # bottom (16, 17, 18, 19)
			2, 3, 6, 7  # top (20, 21, 22, 23)
		]

		const u = UInt32[
			0, 1, 2, 3, # back (0, 1, 2, 3)
			0, 1, 2, 3, # front (4, 5, 6, 7)
			0, 1, 2, 3, # left (8, 9, 10, 11)
			0, 1, 2, 3, # right (12, 13, 14, 15)
			0, 1, 2, 3, # bottom (16, 17, 18, 19)
			0, 1, 2, 3  # top (20, 21, 22, 23)
		]

		const C = UInt32[ #36
			0 + 1, 0 + 0, 0 + 2, 0 + 2, 0 + 3, 0 + 1, # back (1, 0, 2, 2, 3, 1)
			4 + 0, 4 + 1, 4 + 3, 4 + 3, 4 + 2, 4 + 0, # front (4, 5, 7, 7, 6, 4)
			8 + 0, 5 + 4, 9 + 2, 9 + 2, 8 + 2, 8 + 0, # left (0, 4, 6, 6, 2, 0)
			12 + 1, 12 + 0, 12 + 2, 12 + 2, 12 + 3, 12 + 1, # right (5, 1, 3, 3, 7, 5)
			16 + 0, 16 + 1, 16 + 3, 16 + 3, 16 + 2, 16 + 0, # bottom (0, 1, 5, 5, 4, 0)
			20 + 1, 20 + 0, 20 + 2, 20 + 2, 20 + 3, 20 + 1 # top (3, 2, 6, 6, 7, 3)
		]

		const C2 = UInt32[ 	#35
			17, 19, 16, 18,		# bottom
			9, 11, 8, 10,			# left
			20, 0,						# skip
			0, 2, 1, 3,				# back
			14, 2,						# skip
			20, 22, 21, 23,		# top
			15, 13, 14, 12,		# right
			1, 5,							# skip
			5, 7, 4, 6,				# front
		]

		Ilen=length(i)
		Clen=length(C)

		vertices=zeros(Vec3f, Ilen)
		normals=zeros(Vec3f, Ilen)
		uvs=zeros(Vec2f, Ilen)
		indicies=zeros(UInt32, Clen)

		uvi = UInt32(0)
		for vi = 1:Ilen
			uvi+=1
			if uvi>3 uvi = 0 end
			index = i[vi] + 1

			vertices[vi] = v[index] + 1
			normals[vi] = n[index] + 1
			uvs[vi] = uv[u[vi] + 1] + 1
		end

		#uvi = UInt32(0)
		for vi=1:Clen
			#uvi+=1
			indicies[vi] = C[vi] + 1 #index
		end

		#shape->setType(DRAWS::draw_TRIANGLE_STRIP);
		(:TRIANGLE_STRIP, vertices, normals, uvs, indicies)
end

#=
createCubeData2(wdetail::Integer, hdetail::Integer)
	clearVAO();
	_uint_ index_count = 36; // 14

	if (index_count == 14) shape->setType(DRAWS::draw_TRIANGLE_STRIP);
	#else shape->setType(DRAWS::draw_TRIANGLES);

	_uint W = 1; // wdetail < 1 ? 1 : wdetail;
	_uint H = 1; //hdetail < 1 ? 1 : hdetail;
	_uint L = 1;

	_float dW = 1.f / W;
	_float dH = 1.f / H;
	_float dL = 1.f / L;

	_uint MW = W + 1;
	_uint MH = H + 1;
	_uint ML = L + 1;

	_int neg = 0;

	for (_uint d = 0; d < ML; ++d)
	{
		_float D = d * dL;
		++neg;

		for (_uint h = 0; h < MH; ++h)
		{
			_float V = h * dH;

			for (_uint w = 0; w < MW; ++w)
			{
				_float U = w * dW;

				_float x = -1 + U * 2;
				_float y = -1 + V * 2;
				_float z = -1 + D * 2;

				vec3 v = vec3(x, y, z);
				vec3 n = glm::normalize(v);
				vec2 uv = vec2(sin(U), V);

				shape->addPosition(v);
				shape->addNormal(n);
				shape->addTexcoord(uv);
			}
		}
	}

	_uint_* arr = NULL;

	if (index_count == 36)
	{
		_uint_ C[36] =
		{
			4, 5, 7, 7, 6, 4, // front side
			1, 0, 2, 2, 3, 1, // back side
			0, 4, 6, 6, 2, 0, // left side
			5, 1, 3, 3, 7, 5, // right side
			0, 1, 5, 5, 4, 0, // down side
			3, 2, 6, 6, 7, 3  // up side
		};

		arr = C;
	}
	else if (index_count == 14)
	{
		_uint_ C[14] =
		{
			1, 5, 0, 4, 6, 5, 7,
			3, 6, 2, 0, 3, 1, 5
		};

		arr = C;
	}

	_uint MWH = MW * MH;
	_uint x = 0; _uint y = 0; _uint z = 0;

	for (; z < L; ++z)
	{
		_uint b = (z + 0) * MWH; _uint f = (z + 1) * MWH;

		for (y = 0; y < H; ++y)
		{
			_uint d = (y + 0) * MW; _uint u = (y + 1) * MW;
			_uint db = d + b; _uint ub = u + b;
			_uint df = d + f; _uint uf = u + f;

			for (x = 0; x < W; ++x)
			{
				_uint l = (x + 0); _uint r = (x + 1);

				_uint ldb = l + db; _uint lub = l + ub; _uint ldf = l + df; _uint luf = l + uf;
				_uint rdb = r + db; _uint rub = r + ub; _uint rdf = r + df;  _uint ruf = r + uf;

				_uint I[8] = { ldb, rdb, lub, rub, ldf, rdf, luf, ruf };

				if (arr) for (_uint i = 0; i < index_count; ++i)
					shape->addIndex(I[arr[i]]);
			}
		}
	}
end

createTerrain(wdetail::Integer, hdetail::Integer)
	clearVAO();

	#srand(static_cast<unsigned int>(time(0))); # initialize random seed:

#	for (_uint i = 0; i < 60; ++i)
#		for (_uint j = 0; j < 60; ++j)
#			_float v1 = (_float)((sin((random(1000))*3.14f)));
#			_float v2 = (_float)((sin((random(1000))*3.14f)));
#			//_DEBUG_MESSAGE("Radom Numbers: %f, %f", v1, v2);
#			Gradient[i][j][0] = v1;
#			Gradient[i][j][1] = v2;
#		end
#	end

	generateTerrain();

	_uint W = wdetail < 1 ? 1 : wdetail;
	_uint H = hdetail < 1 ? 1 : hdetail;

	_float dW = 1.f / W;
	_float dH = 1.f / H;

	_uint MW = W + 1;
	_uint MH = H + 1;

	for (_uint h = 0; h < MH; ++h)
		_float V = h * dH;

		for (_uint w = 0; w < MW; ++w)
			_float U = w * dW;

			_float x = -1 + U * 2; // -cos(theta) * sin(phi);
			_float y = -1 + V * 2; // -cos(phi);

			vec3 v(x, y, terrain[w][h]);

			shape->addPosition(v);
			shape->addNormal(glm::normalize(v));
			shape->addTexcoord(glm::rotate(vec2(V, U), 3.14f*-.5f));

			if (h < H && w < W)
				shape->addIndex(w + h * MW);
				shape->addIndex(w + 1 + h * MW);
				shape->addIndex(w + (h + 1) * MW);
				shape->addIndex(w + 1 + (h + 1) * MW);
			end
		end
	end
end
=#

function addTriangle(list::Dict{Symbol,Any}, id::Integer, index::Integer, v::Vec3f)
	ids = list[:IDS]
	position = list[:VERTICES]
	normal = list[:NORMALS]
	texcoord = list[:UVS]

	position[ids[1]] = v
	ids[1]+=1

	n = normalize(v)

	normal[ids[2]] = n
	ids[2]+=1

	x = 1 + n[1]
	z = 1 + n[3]

	if z < 1 x = 4 + -x end
	#if ids[4] < x ids[4] = x end
	#if !(x < 0 && z < 0) x = z = 2 end

	#n2 = normalize(uid * 0.5f0)
	#if z > .5f x = mx end

	texcoord[ids[3]] = Vec2f(x * 0.25f0, 0.5f0 + n[2] * -0.5f0)
	ids[3]+=1
end

function addTriangle(list::Dict{Symbol,Any}, id::Integer, v::Array{Vec3f,1})
	for i=1:length(v) addTriangle(list, id, i, v[i]) end
end

function subdivide(list::Dict{Symbol,Any}, id::Integer, depth::Integer, v::Array{Vec3f,1})
	if depth <= 0
		addTriangle(list, id, v)
		return
	end

	v1 = v[1]
	v2 = v[2]
	v3 = v[3]

	v12 =  normalize(v1 + v2)
	v23 =  normalize(v2 + v3)
	v31 =  normalize(v3 + v1)

	#matrix 3x4 vec3f
	vs = [
		[v1, v12, v31],
		[v2, v23, v12],
		[v3, v31, v23],
		[v12, v23, v31]
	]

	#index=1+(i-1) * 3
	#vs = [v1,v2,v3,v12]

	depth = depth - 1

	for i=1:length(vs)
		subdivide(list, id, depth, vs[i])
	end
end

function computeVAO(list::Dict{Symbol,Any}, sides::Integer, depth::Integer, points::Array{Vec3f,1}, indices::Array{Vec3u,1})
	if depth>8
		error("VAO depth should not be larger than 8! Computer will freeze!")
		return
	end

	# Compute and store vertices
	numVertices = sides * 3 #/ sides *  points for one vertex
	VAOSize = numVertices * 3 * Integer(4^depth) # *(x,y,z)

	println("size: $sides*3=$numVertices*3=$(numVertices*3)*4^$depth=$VAOSize")

	list[:IDS] = [1,1,1,0]
	list[:VERTICES] = fill(zeros(Vec3f),VAOSize)
	list[:NORMALS] = fill(zeros(Vec3f),VAOSize)
	list[:UVS] = fill(zeros(Vec2f),VAOSize)
	#list[:INDICIES] = fill(UInt32(0),(numP - 1) * (numSides + 1) * 6)

	for i=1:sides
		vs = [ points[1+indices[i][1]], points[1+indices[i][2]], points[1+indices[i][3]] ]
		subdivide(list, i, depth, vs)
	end
end

function createTetrahedronSphere(depth::Integer)
	X = 0.525731112119133606f0; # sin(inc * 6.362117953f0); // 0.61803399f0,0.6662394341f0,6.362117953f0
	Z = 0.850650808352039932f0; # 1.61803399f0; // cos(0.0f0);

	list = Dict{Symbol,Any}()

	points = [
		Vec3f(X, 0, -Z), Vec3f(-X, 0, -Z), Vec3f(X, 0, Z), Vec3f(-X, 0, Z),
		Vec3f(0, -Z, -X), Vec3f(0, -Z, X), Vec3f(0, Z, -X), Vec3f(0, Z, X),
		Vec3f(-Z, -X, 0), Vec3f(Z, -X, 0), Vec3f(-Z, X, 0), Vec3f(Z, X, 0),
	]

	indices = [
		Vec3u(0, 4, 1), Vec3u(0, 9, 4), Vec3u(9, 5, 4), Vec3u(4, 5, 8), Vec3u(4, 8, 1),
		Vec3u(8, 10, 1), Vec3u(8, 3, 10), Vec3u(5, 3, 8), Vec3u(5, 2, 3), Vec3u(2, 7, 3),
		Vec3u(7, 10, 3), Vec3u(7, 6, 10), Vec3u(7, 11, 6), Vec3u(11, 0, 6), Vec3u(0, 1, 6),
		Vec3u(6, 1, 10), Vec3u(9, 0, 11), Vec3u(9, 11, 2), Vec3u(9, 2, 5), Vec3u(7, 2, 11)
	]

	computeVAO(list, 20, depth, points, indices)

	#(:TRIANGLES, list[:VERTICES], list[:NORMALS], list[:UVS], UInt32[])
	(:TRIANGLES, list[:VERTICES], [], [], [])
end

max_stars_repo_name: Gilga/JGE

max_issues_repo_name: Gilga/JGE

max_forks_repo_name: Gilga/JGE


text:
export PrecomputedAffinityMatrix, PrecomputedKnns

"""
    struct PrecomputedAffinityMatrix <: AbstractSearchContext
        dists # precomputed distances for all pairs (squared matrix)
    end

An index-like wrapper for precomputed affinity matrix.
"""
struct PrecomputedAffinityMatrix{MType<:AbstractMatrix} <: AbstractSearchContext
    dists::MType
end

function SimilaritySearch.search(p::PrecomputedAffinityMatrix, q::Integer, res::KnnResult)
    D = @view p.dists[:, q]
    @inbounds for i in eachindex(D)
        push!(res, i, D[i])
    end

    res
end

"""
    struct PrecomputedKnns <: AbstractSearchContext
        knns
        dists
    end

An index-like wrapper for precomputed all-knns (as knns and dists matrices (k, n))
"""
struct PrecomputedKnns{KnnsType<:AbstractMatrix,DistsType<:AbstractMatrix,DBType<:AbstractVector} <: AbstractSearchContext
    knns::KnnsType
    dists::DistsType
    db::DBType
end

function PrecomputedKnns(knns, dists)
    @assert size(knns) == size(dists)
    PrecomputedKnns(knns, dists, 1:size(knns, 2))
end

SimilaritySearch.getpools(::PrecomputedKnns) = SimilaritySearch.GlobalKnnResult

function SimilaritySearch.search(p::PrecomputedKnns, q::Integer, res::KnnResult)
    N = @view p.knns[:, q]
    D = @view p.dists[:, q]

    @inbounds for i in eachindex(N)
        push!(res, N[i], D[i])
    end

    res
end


text:
module ProjectMoana
export PCG, pcg_init, pcg_rand, pcg_randf

export printcol, luminosity

export HdrImage, InvalidPfmFileFormat,
       _parse_img_size, _parse_endianness, _read_float, _read_line,
       valid_coordinates, pixel_offset, get_pixel, set_pixel,
       average_luminosity, clamp_image, normalize_image, read_pfm_image, tone_mapping

export Vec, Point, Transformation, Normal, Vec2D,
       cross, squared_norm, norm, normalize, inverse, is_consistent, toVec, create_onb,
       translation, scaling, rotation_x, rotation_y, rotation_z

export BRDF, Material, Pigment, Renderer, UniformPigment, get_color, ImagePigment, CheckeredPigment, DiffuseBRDF, SpecularBRDF, scatter_ray

export Shape, Sphere, World, Plane, AAB, LightPoint, HitRecord, ray_intersection, add_shape, Ray, at

export Camera, OrthogonalCamera, PerspectiveCamera,  ImageTracer, Renderer,
       fire_ray, fire_all_rays

export  OnOff_Renderer, Flat_Renderer, PathTracer_Renderer, PointLight_Renderer, OnOff, Flat, PathTracer, PointLight

export SourceLocation, Stop, Identifier, LiteralString, LiteralNumber, Symbol, Keyword, Token, InputStream, GrammarError, read_char,
       unread_char, skip_whitespaces_and_comments, read_token, KeywordEnum, isdigit, NEW, MATERIAL, PLANE, SPHERE, DIFFUSE, SPECULAR,
       UNIFORM, CHECKERED, IMAGE, IDENTITY, TRANSLATION, ROTATION_X, ROTATION_Y, ROTATION_Z, SCALING, CAMERA, ORTHOGONAL, PERSPECTIVE,
       FLOAT, assert_is_identifier, assert_is_keyword, assert_is_number, assert_is_string, assert_is_symbol, Scene, parse_scene

using LinearAlgebra

include("PCG.jl")
include("Colors.jl")
include("HdrImages.jl")
include("Geometry.jl")
include("Surface.jl")
include("Shape.jl")
include("Cameras.jl")
include("Renderer.jl")
include("Lexer.jl")

greet(name) = println("Hello $(name)! Moana welcomes you!")
greet() = println("Hello User! Moana welcomes you!")

end # module


text:
### An alternative interface is via the function wrappers provided below
### This third interface fulfills the following purposes:
### i) It acts as a simplified interface to facilitate interaction with other packages
### ii) It provides a shorter syntax by creating the unerlying model, sampler and runner types under the hook
### iii) It fits better the mindset of some users

function ARS(f::Function, init::Vector{Float64}, nsteps::Int, burnin::Int,
  logproposal::Function, proposalscale::Float64, jumpscale::Float64;
  nchains::Int=1, thinning::Int=1)
  mcmodel::MCLikelihood = model(f, init=init)
  mcsampler::ARS = ARS(logproposal, proposalscale, jumpscale)
  mcrunner::SerialMC = SerialMC(burnin=burnin, thinning=thinning, nsteps=nsteps)
  mcsamples::Array{Float64, 3} = Array(Float64, length(mcrunner.r), length(init), nchains)

  for i = 1:nchains
    mcsamples[:, :, i] = run(mcmodel, mcsampler, mcrunner).samples
  end

  mcsamples
end

function SliceSampler(f::Function, init::Vector{Float64}, nsteps::Int, burnin::Int, widths::Vector{Float64};
  nchains::Int=1, thinning::Int=1, stepout::Bool=true)
  mcmodel::MCLikelihood = model(f, init=init)
  mcsampler::SliceSampler = SliceSampler(widths, stepout)
  mcrunner::SerialMC = SerialMC(burnin=burnin, thinning=thinning, nsteps=nsteps)
  mcsamples::Array{Float64, 3} = Array(Float64, length(mcrunner.r), length(init), nchains)

  for i = 1:nchains
    mcsamples[:, :, i] = run(mcmodel, mcsampler, mcrunner).samples
  end

  mcsamples
end

function MH(f::Function, init::Vector{Float64}, nsteps::Int, burnin::Int;
  nchains::Int=1,
  thinning::Int=1,
  logproposal::FunctionOrNothing=nothing,
  randproposal::FunctionOrNothing=(x::Vector{Float64} -> rand(IsoNormal(x, 1.))))
  mcmodel::MCLikelihood = model(f, init=init)
  mcsampler::MH
  if logproposal==nothing
    mcsampler = MH(randproposal)
  else
    mcsampler = MH(logproposal, randproposal)
  end
  mcrunner::SerialMC = SerialMC(burnin=burnin, thinning=thinning, nsteps=nsteps)
  mcsamples::Array{Float64, 3} = Array(Float64, length(mcrunner.r), length(init), nchains)

  for i = 1:nchains
    mcsamples[:, :, i] = run(mcmodel, mcsampler, mcrunner).samples
  end

  mcsamples
end

function RAM(f::Function, init::Vector{Float64}, nsteps::Int, burnin::Int;
  nchains::Int=1, thinning::Int=1, jumpscale::Float64=1., targetrate::Float64=0.234)
  mcmodel::MCLikelihood = model(f, init=init)
  mcsampler::RAM = RAM(jumpscale, targetrate)
  mcrunner::SerialMC = SerialMC(burnin=burnin, thinning=thinning, nsteps=nsteps)
  mcsamples::Array{Float64, 3} = Array(Float64, length(mcrunner.r), length(init), nchains)

  for i = 1:nchains
    mcsamples[:, :, i] = run(mcmodel, mcsampler, mcrunner).samples
  end

  mcsamples
end

function HMC(f::Function, g::Function, init::Vector{Float64}, nsteps::Int, burnin::Int;
  nchains::Int=1, thinning::Int=1, nleaps::Int=10, leapstep::Float64=0.1)
  mcmodel::MCLikelihood = model(f, grad=g, init=init)
  mcsampler::HMC = HMC(nleaps, leapstep)
  mcrunner::SerialMC = SerialMC(burnin=burnin, thinning=thinning, nsteps=nsteps)
  mcsamples::Array{Float64, 3} = Array(Float64, length(mcrunner.r), length(init), nchains)

  for i = 1:nchains
    mcsamples[:, :, i] = run(mcmodel, mcsampler, mcrunner).samples
  end

  mcsamples
end

function MALA(f::Function, g::Function, init::Vector{Float64}, nsteps::Int, burnin::Int;
  nchains::Int=1, thinning::Int=1, driftstep::Float64=1.)
  mcmodel::MCLikelihood = model(f, grad=g, init=init)
  mcsampler::MALA = MALA(driftstep=driftstep)
  mcrunner::SerialMC = SerialMC(burnin=burnin, thinning=thinning, nsteps=nsteps)
  mcsamples::Array{Float64, 3} = Array(Float64, length(mcrunner.r), length(init), nchains)

  for i = 1:nchains
    mcsamples[:, :, i] = run(mcmodel, mcsampler, mcrunner).samples
  end

  mcsamples
end

function SMMALA(f::Function, g::Function, h::Function, init::Vector{Float64}, nsteps::Int, burnin::Int;
  nchains::Int=1, thinning::Int=1, driftstep::Float64=1.)
  mcmodel::MCLikelihood = model(f, grad=g, tensor=h, init=init)
  mcsampler::SMMALA = SMMALA(driftstep=driftstep)
  mcrunner::SerialMC = SerialMC(burnin=burnin, thinning=thinning, nsteps=nsteps)
  mcsamples::Array{Float64, 3} = Array(Float64, length(mcrunner.r), length(init), nchains)

  for i = 1:nchains
    mcsamples[:, :, i] = run(mcmodel, mcsampler, mcrunner).samples
  end

  mcsamples
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module PkgTests

import Pkg

if (server = Pkg.pkg_server()) !== nothing && Sys.which("curl") !== nothing
    s = read(`curl -sLI $(server)`, String);
    @info "Pkg Server metadata:\n$s"
end

# Make sure to not start with an outdated registry
rm(joinpath(@__DIR__, "registries"); force = true, recursive = true)

include("utils.jl")
include("new.jl")
include("pkg.jl")
include("repl.jl")
include("api.jl")
include("registry.jl")
include("subdir.jl")
include("artifacts.jl")
include("binaryplatforms.jl")
include("platformengines.jl")
include("sandbox.jl")
include("resolve.jl")

# clean up locally cached registry
rm(joinpath(@__DIR__, "registries"); force = true, recursive = true)

end # module


text:
#
# 2D Position operators
#

function x_coord_matrix(p::Param)

    #Variables
    xmin        = p.xmin
    xmax        = p.xmax
    Nx          = p.xnodes

    ymin        = p.ymin
    ymax        = p.ymax
    Ny          = p.ynodes

     #Discretize space
     dx      = (xmax - xmin)/(Nx-1)
     x       = [xmin + (i-1)*dx for i in 1:Nx]

     dy      = (ymax - ymin)/(Ny-1)
     y       = [ymin + (i-1)*dy for i in 1:Ny]

     if(p.Boundaries == :periodic )
         dx      = (xmax - xmin)/(Nx)
         x       = [xmin + (i-1)*dx for i in 1:Nx]

         dy      = (ymax - ymin)/(Ny)
         y       = [ymin + (i-1)*dy for i in 1:Ny]
     end

    A = spzeros(Nx*Ny , Nx*Ny)

    index = 0

    for (ix , posx) in enumerate(x) , (iy , posy) in enumerate(y)

        index = (iy-1)*Nx+ix
        A[index , index] = posx

    end

    return A
end

function y_coord_matrix(p::Param)

    #Variables
    xmin        = p.xmin
    xmax        = p.xmax
    Nx          = p.xnodes

    ymin        = p.ymin
    ymax        = p.ymax
    Ny          = p.ynodes

    #Discretize space
    dx      = (xmax - xmin)/(Nx-1)
    x       = [xmin + (i-1)*dx for i in 1:Nx]

    dy      = (ymax - ymin)/(Ny-1)
    y       = [ymin + (i-1)*dy for i in 1:Ny]

    if(p.Boundaries == :periodic )
        dx      = (xmax - xmin)/(Nx)
        x       = [xmin + (i-1)*dx for i in 1:Nx]

        dy      = (ymax - ymin)/(Ny)
        y       = [ymin + (i-1)*dy for i in 1:Ny]
    end

    A = spzeros(Nx*Ny , Nx*Ny)

    index = 0

    for (ix , posx) in enumerate(x) , (iy , posy) in enumerate(y)

        index = (iy-1)*Nx+ix
        A[index , index] = posy

    end

    return A
end


text:
using ModelingToolkit, OrdinaryDiffEq, Test

@parameters t α β γ δ
@variables x(t) y(t)
D = Differential(t)

eqs = [D(x) ~ α * x - β * x * y,
    D(y) ~ -δ * y + γ * x * y]

sys = ODESystem(eqs)

u0 = [x => 1.0,
    y => 1.0]

p = [α => 1.5,
    β => 1.0,
    δ => 3.0,
    γ => 1.0]

tspan = (0.0, 10.0)
prob = ODEProblem(sys, u0, tspan, p)
sol = solve(prob, Tsit5())

sys2 = liouville_transform(sys)
@variables trJ

u0 = [x => 1.0,
    y => 1.0,
    trJ => 1.0]

prob = ODEProblem(sys2, u0, tspan, p, jac = true)
sol = solve(prob, Tsit5())
@test sol[end, end] ≈ 1.0742818931017244


text:
module HydrostaticFreeSurfaceModels

using KernelAbstractions: @index, @kernel, Event, MultiEvent
using KernelAbstractions.Extras.LoopInfo: @unroll

using Oceananigans.Utils: launch!

import Oceananigans: fields

#####
##### HydrostaticFreeSurfaceModel definition
#####

include("compute_w_from_continuity.jl")
include("explicit_free_surface.jl")
include("implicit_free_surface.jl")
include("implicit_free_surface_solver.jl")
include("implicit_free_surface_solver_kernels.jl")
include("rigid_lid.jl")
include("hydrostatic_free_surface_velocity_fields.jl")
include("hydrostatic_free_surface_tendency_fields.jl")
include("hydrostatic_free_surface_model.jl")
include("show_hydrostatic_free_surface_model.jl")
include("set_hydrostatic_free_surface_model.jl")

#####
##### Time-stepping HydrostaticFreeSurfaceModels
#####

"""
    fields(model::HydrostaticFreeSurfaceModel)

Returns a flattened `NamedTuple` of the fields in `model.velocities` and `model.tracers`.
"""
fields(model::HydrostaticFreeSurfaceModel) = hydrostatic_prognostic_fields(model.velocities, model.free_surface, model.tracers)

"""
    fields(model::HydrostaticFreeSurfaceModel)

Returns a flattened `NamedTuple` of the fields in `model.velocities` and `model.tracers`.
"""
hydrostatic_prognostic_fields(velocities, free_surface, tracers) = merge((u = velocities.u,
                                                                          v = velocities.v,
                                                                          η = free_surface.η),
                                                                          tracers)

include("barotropic_pressure_correction.jl")
include("hydrostatic_free_surface_advection.jl")
include("hydrostatic_free_surface_tendency_kernel_functions.jl")
include("calculate_hydrostatic_free_surface_tendencies.jl")
include("update_hydrostatic_free_surface_model_state.jl")
include("hydrostatic_free_surface_ab2_step.jl")
include("prescribed_hydrostatic_velocity_fields.jl")

#####
##### Some diagnostics
#####

include("vertical_vorticity_field.jl")

end # module


text:
"Tableau for the s-stage Gauss-Lobatto SFIRK method"
function getTableauStochasticGLRK(s::Int)

    if s==1
        name = :StochasticGLRK1
    elseif s==2
        name = :StochasticGLRK2
    elseif s==3
        name = :StochasticGLRK3
    elseif s==4
        name = :StochasticGLRK4
    elseif s==5
        name = :StochasticGLRK5
    else s==6
        name = :StochasticGLRK6
    end

    TableauSIRK(name,getCoefficientsGLRK(s),getCoefficientsGLRK(s))
end

"""
Tableau for the 2-stage stochastic symplectic DIRK method
  Tableau for the stochastic symplectic DIRK method
  Satisfies the conditions for Lagrange-d'Alembert integrators.
  Satisfies the conditions for strong convergence of order 1.0 for one Wiener process
"""
function getTableauStochasticDIRK(c::Number=0.5)

    a_drift = [[c/2  0.0    ]
                [c    (1-c)/2]]

    b_drift = [c, 1-c]

    c_drift = [c/2, (1+c)/2]

    TableauSIRK(:StochasticSymplecticDIRK, 2, a_drift, b_drift, c_drift, 2, a_drift, b_drift, c_drift)
end


text:
TensorHead(s::TensorHead) = s
TensorHead(name::AbstractString,index_types::AbstractArray{TensorIndexType},symmetry::TensorSymmetry)::TensorHead = tensor.TensorHead(name,index_types,symmetry)
TensorHead(name::AbstractString,index_types::AbstractArray{TensorIndexType})::TensorHead = tensor.TensorHead(name,index_types)
TensorHead(name::Symbol,index_types::AbstractArray{TensorIndexType},symmetry::TensorSymmetry)::TensorHead = tensor.TensorHead(name,index_types,symmetry)
TensorHead(name::Symbol,index_types::AbstractArray{TensorIndexType})::TensorHead = tensor.TensorHead(name,index_types)
(t::TensorHead)(ics::TensorIndex...) = convert(IndexedTensor,t.__pyobject__(ics...))

function Base.getproperty(A::TensorHead, k::Symbol)
    if k == :index_types
        m = getproperty(PyCall.PyObject(A),k)
        p = PyCall.PyObject(m)
        return TensorIndexType[ x for x in p]
    elseif k in fieldnames(typeof(A))
        return getfield(A,k)
    # else error?
    end
end

function Base.getproperty(A::TensMul, k::Symbol)
    if k == :coeff
        m = getproperty(PyCall.PyObject(A),k)
        return sympy_type_convert(m)
    elseif k == :nocoeff
        m = getproperty(PyCall.PyObject(A),k)
        return sympy_type_convert(m)
    elseif k in fieldnames(typeof(A))
        return getfield(A,k)
    # else error?
    end
end

scalar_exprs = Dict{SymPy.Sym,TensMul}()
scalar_names = Dict{TensMul,SymPy.Sym}()
scalar_name = "ts"
scalar_index = 1
function new_scalar(expr::TensMul)
    global scalar_exprs,scalar_names,scalar_index
    if !(expr in keys(scalar_names))
        name = scalar_name * string(scalar_index)
        symb = symbols(name,positive=true)
        scalar_exprs[symb] = expr
        scalar_names[expr] = symb
        scalar_index += 1
        return symb
    else
        return scalar_names[expr]
    end
end
function get_scalars(texpr::T)  where {T <: SymbolicObject}
    vars = [var for var in free_symbols(texpr) if startswith(string(var),scalar_name)]
    return vars
end

#function Base.:+(A::Number, B::TensScalar)
#    return TensScalar(B.var,B.expr+A,B.__pyobject__)
#end
#function Base.:-(A::Number, B::TensScalar)
#    return TensScalar(B.var,B.expr-A,B.__pyobject__)
#end
#function Base.:/(A::Number, B::TensScalar)
#    return TensScalar(B.var,A/B.expr,B.__pyobject__)
#end
#function Base.:/(A::TensScalar,B::Number)
#    return TensScalar(A.var,A.expr/B,A.__pyobject__)
#end
#function Base.:*(A::TensScalar,B::Number)
#    return TensScalar(A.var,A.expr*B,A.__pyobject__)
#end
#
#function Base.:*(A::TensScalar,B::T) where {T <: Tensor}
#    return Base.:*(B,A)
#end
#function Base.:*(B::T,A::TensScalar) where {T <:Tensor}
#    return TensMul(A.expr*B)
#end
#function Base.:^(A::TensScalar,B::Number)
#    return TensScalar(A.var,A.expr^B,A.__pyobject__)
#end
#function Base.:^(A::TensScalar,B::Integer)
#    return TensScalar(A.var,A.expr^B,A.__pyobject__)
#end
#function Base.:^(A::Number,B::TensScalar)
#    return TensScalar(B.var,A^B.expr,B.__pyobject__)
#end
#function Base.log(B::TensScalar)
#    return TensScalar(B.var,log(B.expr),B.__pyobject__)
#end
#
function Base.show(io::IO, ::MIME"text/plain", t::T) where {T <: Tensor}
    print(io, sympy.pretty(t.__pyobject__))
    sc = get_scalars(t)
    if length(sc) > 0
        println(io,"")
        for s in sc
            print(io, sympy.pretty(sympy.Eq(s,scalar_exprs[s])))
            println(io,"")
        end
    end
    return nothing
end
function Base.show(io::IO, ::MIME"text/latex", t::T) where {T <: Tensor}
    print(io, sympy.latex(t.__pyobject__, mode="equation*"))
    sc = get_scalars(t)
    if length(sc) > 0
        println(io,"")
        for s in sc
            print(io, sympy.latex(sympy.Eq(s,scalar_exprs[s]), mode="equation*"))
            println(io,"")
        end
    end
    return nothing
end
function Base.show(io::IO, ::MIME"text/plain", s::TensorSymmetry)
    if s.__pyobject__.__class__.__name__ == "TensorSymmetry"
        print(io, pycall(s.__pyobject__.__repr__,String))
    else
        print(io, s.__pyobject__)
    end
end
function Base.show(io::IO, ::MIME"text/plain", s::SymbolicObject)
    print(io, sympy.pretty(s))
    sc = get_scalars(s)
    if length(sc) > 0
        println(io,"")
        for s in sc
            print(io, sympy.pretty(sympy.Eq(s,scalar_exprs[s])))
            println(io,"")
        end
    end
end
function Base.show(io::IO, ::MIME"text/latex", s::SymbolicObject)
    print(io, sympy.latex(s, mode="equation*"))
    sc = get_scalars(s)
    if length(sc) > 0
        println(io,"")
        for s in sc
            print(io, sympy.latex(sympy.Eq(s,scalar_exprs[s]), mode="equation*"))
            println(io,"")
        end
    end
end
#function Base.show(io::IO,  t::TensScalar)
#    show(io,t.expr)
#    print(io,",")
#    print(io, jprint(sympy.Eq(t.var,t.__pyobject__)))
#    return nothing
#end

TensAdd(s::TensAdd) = s
TensMul(s::TensMul) = s
(t::TensMul)(ics::TensorIndex...) = convert(TensMul,t.__pyobject__(ics...))
(t::TensAdd)(ics::TensorIndex...) = convert(TensAdd,t.__pyobject__(ics...))
(t::IndexedTensor)(ics::TensorIndex...) = convert(IndexedTensor,t.__pyobject__(ics...))
function terms(s::T) where T <: Union{TensAdd,TensMul}
    sp = getproperty(PyCall.PyObject(s),:args)
    ret = []
    for x in sp
        push!(ret,sympy_type_convert(x))
    end
    return ret
end

function Base.:+(A::T ,B::U) where {T <: Tensor, U <: SymbolicObject}
    pyexp =  A.__pyobject__+B.__pyobject__
    return sympy_type_convert(pyexp)
end
function Base.:-(A::T ,B::U) where {T <: Tensor, U <: SymbolicObject}
    pyexp =  A.__pyobject__-B.__pyobject__
    return sympy_type_convert(pyexp)
end
## Eg. A.norm() where A = [x 1; 1 x], say
function Base.getproperty(A::TensorSymmetry, k::Symbol)
    if k == :fully_symmetric || k == :direct_product || k == :no_symmetry
        return x->pycall(getproperty(A.__pyobject__,k), TensorSymmetry, x)
    elseif k == :riemann
        return pycall(getproperty(A.__pyobject__,k), TensorSymmetry)
    elseif k in fieldnames(typeof(A))
        return getfield(A,k)
    end
end

function Base.:*(A::T ,B::U) where {T <: Tensor, U <: Tensor}
    pyexp = A.__pyobject__*B.__pyobject__
    inds = pycall(getproperty(pyexp,:get_free_indices),Array{TensorIndex})
    if length(inds) == 0
        x = new_scalar(convert(TensMul,pyexp))
        return x
    else
        return convert(TensMul,pyexp)
    end
end
function Base.:*(A::T ,B::U) where {T <: Number, U <: Tensor}
    pyexp = A*B.__pyobject__
    return convert(TensMul,pyexp)
end
function Base.:*(A::T ,B::U) where {T <: Tensor, U <: Number}
    pyexp = B*A.__pyobject__
    return convert(TensMul,pyexp)
end
function Base.:*(A::T ,B::U) where {T <: SymbolicObject, U <: Tensor}
    pyexp = A*B.__pyobject__
    return convert(TensMul,pyexp)
end
function Base.:*(A::T ,B::U) where {T <: Tensor, U <: SymbolicObject}
    pyexp = B*A.__pyobject__
    return convert(TensMul,pyexp)
end
function Base.:/(B::U,A::T) where {T <: Number, U <: Tensor}
    pyexp = B.__pyobject__/A
    return sympy_type_convert(pyexp)
end
function Base.:/(B::U,A::T) where {T <: SymbolicObject, U <: Tensor}
    pyexp = B.__pyobject__/A
    return sympy_type_convert(pyexp)
end

function (==)( a::T ,b::Real) where T <: Tensor
    return false
end
function (==)(a::T,b::U) where {T <: Tensor, U <: SymbolicObject}
    return false
end
function (==)(a::T,b::U) where {T <: Tensor, U <: Tensor}
    return a.__pyobject__ == b.__pyobject__
end
function scalarIsEqual(a::U,b::T,metric::TensorHead) where {T <: Union{SymbolicObject,Tensor}, U <: Union{SymbolicObject,Tensor}}
    sca = get_scalars(a)
    scb = get_scalars(b)
    if 2 > length(sca) && 2 > length(scb)
        exa = canon_bp(contract_metric(scalar_exprs[sca[1]],metric))
        exb = canon_bp(contract_metric(scalar_exprs[scb[1]],metric))
        return exa == exb
    elseif length(sca) > 0 || length(scb) > 0
        error("Multiple scalars, not implemented yet.")
        return
    else
        return sympy.Eq(a,b) == True
    end
end

macro heads(IT::AbstractArray{TensorIndexType},xs...)
    hs = []
    for x in xs
        push!(hs,TensorHead(x,IT))
    end
    return hs
end

function canon_bp(s::Sym)
    sc = get_scalars(s)
    subs = Dict()
    if length(sc) > 0
        for ss in sc
            subs[ss] = canon_bp(scalar_exprs[ss])
        end
    end
    pyexp = s
    for k in keys(subs)
        if typeof(subs[k]) != TensMul
            pyexp = pyexp.subs(k,subs[k])
        end
    end
    return sympy_type_convert(pyexp+0)
end

function canon_bp(t::T) where T <: Tensor
    sc = get_scalars(t)
    subs = Dict()
    if length(sc) > 0
        for ss in sc
            subs[ss] = canon_bp(scalar_exprs[ss])
        end
    end
    pyexp = tensor.canon_bp(t)
    for k in keys(subs)
        if typeof(subs[k]) != TensMul
            pyexp = pyexp.subs(k,subs[k])
        end
    end
    return sympy_type_convert(pyexp+0)
end
function contract_metric(s::Sym,metric::TensorHead)
    sc = get_scalars(s)
    subs = Dict()
    if length(sc) > 0
        for ss in sc
            subs[ss] = contract_metric(scalar_exprs[ss],metric)
        end
    end
    pyexp = tensor.contract_metric(s,metric)
    for k in keys(subs)
        if typeof(subs[k]) != TensMul
            pyexp = pyexp.subs(k,subs[k])
        end
    end
    return sympy_type_convert(pyexp+0)
end
function contract_metric(t::T,metric::TensorHead) where T <: Tensor
    sc = get_scalars(t)
    subs = Dict()
    if length(sc) > 0
        for ss in sc
            subs[ss] = contract_metric(scalar_exprs[ss],metric)
        end
    end
    pyexp = tensor.contract_metric(t,metric)
    for k in keys(subs)
        if typeof(subs[k]) != TensMul
            pyexp = pyexp.subs(k,subs[k])
        end
    end
    return sympy_type_convert(pyexp+0)
end
contract_metric(x, metric::TensorHead) = x
canon_bp(x) = x

function factor(t::T) where T <: Tensor
    return sympy_type_convert(SymPy.factor(convert(Sym,t.__pyobject__)))
end
factor(x) = x
function simplify(t::T) where T <: Tensor
    return sympy_type_convert(SymPy.simplify(convert(Sym,t.__pyobject__)))
end

#function contract_delta(t::TensAdd, delta::TensorHead)
#    return sum([ contract_delta(ti,delta) for ti in terms(t) ])
#end
#function contract_delta(t::TensMul, delta::TensorHead)
#
#end


text:
const QuasiArrayMulArray{p, q, T, V} =
    Applied{<:Any, typeof(*), <:Tuple{<:AbstractQuasiArray{T,p}, <:AbstractArray{V,q}}}

const ArrayMulQuasiArray{p, q, T, V} =
    Applied{<:Any, typeof(*), <:Tuple{<:AbstractArray{T,p}, <:AbstractQuasiArray{V,q}}}

const QuasiArrayMulQuasiArray{p, q, T, V} =
    Applied{<:Any, typeof(*), <:Tuple{<:AbstractQuasiArray{T,p}, <:AbstractQuasiArray{V,q}}}
####
# Matrix * Vector
####
const QuasiMatMulVec{T, V} = QuasiArrayMulArray{2, 1, T, V}
const QuasiMatMulMat{T, V} = QuasiArrayMulArray{2, 2, T, V}
const QuasiMatMulQuasiMat{T, V} = QuasiArrayMulQuasiArray{2, 2, T, V}

function getindex(M::Mul{<:AbstractQuasiArrayApplyStyle}, k::Number)
    A,Bs = first(M.args), tail(M.args)
    B = _mul(Bs...)
    ret = zero(eltype(M))
    for j = rowsupport(A, k) ∩ colsupport(B,1)
        ret += A[k,j] * B[j]
    end
    ret
end

getindex(M::Mul{<:AbstractQuasiArrayApplyStyle}, k::Int) =
    Base.invoke(getindex, Tuple{Mul{<:AbstractQuasiArrayApplyStyle},Number}, M, k)

function _mul_quasi_getindex(M::Mul, k::Number, j::Number)
    A,Bs = first(M.args), tail(M.args)
    B = _mul(Bs...)
    ret = zero(eltype(M))
    @inbounds for ℓ in (rowsupport(A,k) ∩ colsupport(B,j))
        ret += A[k,ℓ] * B[ℓ,j]
    end
    ret
end

getindex(M::Mul{<:AbstractQuasiArrayApplyStyle}, k::Number, j::Number) =
    _mul_quasi_getindex(M, k, j)

getindex(M::Mul{<:AbstractQuasiArrayApplyStyle}, k::Integer, j::Integer) =
    _mul_quasi_getindex(M, k, j)

function getindex(M::QuasiMatMulVec, k::AbstractArray)
    A,B = M.args
    ret = zeros(eltype(M),length(k))
    @inbounds for j in axes(A,2)
        ret .+= view(A,k,j) .* B[j]
    end
    ret
end

*(A::AbstractQuasiArray, B...) = fullmaterialize(apply(*,A,B...))
*(A::AbstractQuasiArray, B::AbstractQuasiArray, C...) = fullmaterialize(apply(*,A,B,C...))
*(A::AbstractArray, B::AbstractQuasiArray, C...) = fullmaterialize(apply(*,A,B,C...))

for op in (:pinv, :inv)
    @eval $op(A::AbstractQuasiArray) = fullmaterialize(apply($op,A))
end

axes(L::Ldiv{<:Any,<:Any,<:Any,<:AbstractQuasiMatrix}) =
    (axes(L.A, 2),axes(L.B,2))
axes(L::Ldiv{<:Any,<:Any,<:Any,<:AbstractQuasiVector}) =
    (axes(L.A, 2),)

\(A::AbstractQuasiArray, B::AbstractQuasiArray) = apply(\,A,B)
\(A::AbstractQuasiArray, B::AbstractArray) = apply(\,A,B)
\(A::AbstractArray, B::AbstractQuasiArray) = apply(\,A,B)

*(A::AbstractQuasiArray, B::Mul, C...) = apply(*,A, B.args..., C...)
*(A::Mul, B::AbstractQuasiArray, C...) = apply(*,A.args..., B, C...)

####
# MulQuasiArray
#####

const MulQuasiArray{T, N, Args<:Tuple} = ApplyQuasiArray{T, N, typeof(*), Args}

const MulQuasiVector{T, Args<:Tuple} = MulQuasiArray{T, 1, Args}
const MulQuasiMatrix{T, Args<:Tuple} = MulQuasiArray{T, 2, Args}

const Vec = MulQuasiVector

_ApplyArray(F, factors...) = ApplyQuasiArray(F, factors...)
_ApplyArray(F, factors::AbstractArray...) = ApplyArray(F, factors...)

MulQuasiOrArray = Union{MulArray,MulQuasiArray}

_factors(M::MulQuasiOrArray) = M.args
_factors(M) = (M,)

_flatten(A::MulQuasiArray, B...) = _flatten(Applied(A), B...)
flatten(A::MulQuasiArray) = ApplyQuasiArray(flatten(Applied(A)))
flatten(A::SubQuasiArray{<:Any,2,<:MulQuasiArray}) = materialize(flatten(Applied(A)))

function fullmaterialize(M::Applied{<:Any,typeof(*)})
    M_mat = materialize(flatten(M))
    typeof(M_mat) <: MulQuasiOrArray || return M_mat
    typeof(Applied(M_mat)) == typeof(M) || return(fullmaterialize(M_mat))

    ABC = M_mat.args
    length(ABC) ≤ 2 && return flatten(M_mat)

    AB = most(ABC)
    Mhead = fullmaterialize(Mul(AB...))

    typeof(_factors(Mhead)) == typeof(AB) ||
        return fullmaterialize(Mul(_factors(Mhead)..., last(ABC)))

    BC = tail(ABC)
    Mtail =  fullmaterialize(Mul(BC...))
    typeof(_factors(Mtail)) == typeof(BC) ||
        return fullmaterialize(Mul(first(ABC), _factors(Mtail)...))

    apply(*,first(ABC), Mtail.args...)
end

fullmaterialize(M::ApplyQuasiArray) = flatten(fullmaterialize(Applied(M)))
fullmaterialize(M) = flatten(M)

adjoint(A::MulQuasiArray) = ApplyQuasiArray(*, reverse(adjoint.(A.args))...)
transpose(A::MulQuasiArray) = ApplyQuasiArray(*, reverse(transpose.(A.args))...)

function similar(A::MulQuasiArray)
    B,a = A.args
    B*similar(a)
end

function similar(A::QuasiArrayMulArray)
    B,a = A.args
    applied(*, B, similar(a))
end

function copy(a::MulQuasiArray)
    @_propagate_inbounds_meta
    copymutable(a)
end

function copyto!(dest::MulQuasiArray, src::MulQuasiArray)
    d = last(dest.args)
    s = last(src.args)
    copyto!(IndexStyle(d), d, IndexStyle(s), s)
    dest
end

_mul_tail_support(j, Z) = maximum(last.(colsupport.(Ref(Z),j)))
_mul_tail_support(j, Z, Y, X...) = _mul_tail_support(OneTo(_mul_tail_support(j,Z)), Y, X...)

function _mul_getindex(k, j, A, B)
    M = min(_mul_tail_support(j,B), maximum(last.(rowsupport.(Ref(A),k))))
    A[k,1:M]*B[1:M,j]
end

function _mul_getindex(k, j, A, B, C, D...)
    N = _mul_tail_support(j, reverse(D)..., C, B)
    M = min(maximum(last.(rowsupport.(Ref(A),k))), N)
    _mul_getindex(OneTo(N), j, A[k,OneTo(M)]*B[OneTo(M),OneTo(N)], C, D...)
end

getindex(A::MulQuasiMatrix, k::AbstractVector{<:Number}, j::AbstractVector{<:Number}) =
    _mul_getindex(k, j, A.args...)

@inline quasimulapplystyle(A...) = combine_mul_styles(A...)

ApplyStyle(::typeof(*), ::Type{A}, ::Type{B}, C::Type...) where {A<:AbstractQuasiArray,B<:Union{AbstractArray,AbstractQuasiArray}} =
    quasimulapplystyle(MemoryLayout(A), MemoryLayout(B), MemoryLayout.(C)...)
ApplyStyle(::typeof(*), ::Type{A}, ::Type{B}, C::Type...) where {A<:AbstractArray,B<:AbstractQuasiArray} =
    quasimulapplystyle(MemoryLayout(A), MemoryLayout(B), MemoryLayout.(C)...)
ApplyStyle(::typeof(*), ::Type{A}, ::Type{B}, ::Type{C}, D::Type...) where {A<:AbstractQuasiArray,B<:Union{AbstractArray,AbstractQuasiArray},C<:Union{AbstractArray,AbstractQuasiArray}} =
    quasimulapplystyle(MemoryLayout(A), MemoryLayout(B), MemoryLayout(C), MemoryLayout.(D)...)
ApplyStyle(::typeof(*), ::Type{A}, ::Type{B}, ::Type{C}, D::Type...) where {A<:AbstractArray,B<:AbstractQuasiArray,C<:Union{AbstractArray,AbstractQuasiArray}} =
    quasimulapplystyle(MemoryLayout(A), MemoryLayout(B), MemoryLayout(C), MemoryLayout.(D)...)
ApplyStyle(::typeof(*), ::Type{A}, ::Type{B}, ::Type{C}, D::Type...) where {A<:AbstractArray,B<:AbstractArray,C<:AbstractQuasiArray} =
    quasimulapplystyle(MemoryLayout(A), MemoryLayout(B), MemoryLayout(C), MemoryLayout.(D)...)
ApplyStyle(::typeof(*), ::Type{A}, ::Type{B}, ::Type{C}, ::Type{D}, E::Type...) where {A<:AbstractQuasiArray,B<:Union{AbstractArray,AbstractQuasiArray},C<:Union{AbstractArray,AbstractQuasiArray},D<:Union{AbstractArray,AbstractQuasiArray}} =
    quasimulapplystyle(MemoryLayout(A), MemoryLayout(B), MemoryLayout(C), MemoryLayout(D), MemoryLayout.(E)...)
ApplyStyle(::typeof(*), ::Type{A}, ::Type{B}, ::Type{C}, ::Type{D}, E::Type...) where {A<:AbstractArray,B<:AbstractQuasiArray,C<:Union{AbstractArray,AbstractQuasiArray},D<:Union{AbstractArray,AbstractQuasiArray}} =
    quasimulapplystyle(MemoryLayout(A), MemoryLayout(B), MemoryLayout(C), MemoryLayout(D), MemoryLayout.(E)...)
ApplyStyle(::typeof(*), ::Type{A}, ::Type{B}, ::Type{C}, ::Type{D}, E::Type...) where {A<:AbstractArray,B<:AbstractArray,C<:AbstractQuasiArray,D<:Union{AbstractArray,AbstractQuasiArray}} =
    quasimulapplystyle(MemoryLayout(A), MemoryLayout(B), MemoryLayout(C), MemoryLayout(D), MemoryLayout.(E)...)
ApplyStyle(::typeof(*), ::Type{A}, ::Type{B}, ::Type{C}, ::Type{D}, E::Type...) where {A<:AbstractArray,B<:AbstractArray,C<:AbstractArray,D<:AbstractQuasiArray} =
    quasimulapplystyle(MemoryLayout(A), MemoryLayout(B), MemoryLayout(C), MemoryLayout(D), MemoryLayout.(E)...)

struct QuasiArrayLayout <: MemoryLayout end
MemoryLayout(::Type{<:AbstractQuasiArray}) = QuasiArrayLayout()
combine_mul_styles(::QuasiArrayLayout) = QuasiArrayApplyStyle()

result_mul_style(::QuasiArrayApplyStyle, ::QuasiArrayApplyStyle) = QuasiArrayApplyStyle()
result_mul_style(::QuasiArrayApplyStyle, ::IdentityMulStyle) = QuasiArrayApplyStyle()
result_mul_style(::IdentityMulStyle, ::QuasiArrayApplyStyle) = QuasiArrayApplyStyle()

####
# Matrix * Array
####

function _lmaterialize(A::MulQuasiArray, B, C...)
    As = A.args
    flatten(_ApplyArray(*, reverse(tail(reverse(As)))..., _lmaterialize(last(As), B, C...)))
end

function _rmaterialize(Z::MulQuasiArray, Y, W...)
    Zs = Z.args
    flatten(_ApplyArray(*, _rmaterialize(first(Zs), Y, W...), tail(Zs)...))
end

####
# Lazy \ ApplyArray. This applies to first arg.
#####

quasimulapplystyle(::InvLayout, _) = LdivApplyStyle()
quasimulapplystyle(::PInvLayout, _) = LdivApplyStyle()

ApplyStyle(::typeof(\), ::Type{<:AbstractQuasiArray}, ::Type{<:Applied}) = LdivApplyStyle()
\(A::AbstractQuasiArray, B::Applied) = apply(\, A, B)
function copy(L::Ldiv{LazyLayout,ApplyLayout{typeof(*)},<:AbstractQuasiMatrix})
    args = arguments(L.B)
    apply(*, L.A \  first(args),  tail(args)...)
end

# copy(A::Applied{QuasiArrayApplyStyle,typeof(*)}) = lmaterialize(A)

import LazyArrays: MulAddStyle, _αAB, scalarone
_αAB(M::Mul{MulAddStyle,<:Tuple{<:AbstractQuasiArray,<:AbstractQuasiArray}}, ::Type{T}) where T = tuple(scalarone(T), M.args...)
_αAB(M::Mul{MulAddStyle,<:Tuple{<:Number,<:AbstractQuasiArray,<:AbstractQuasiArray}}, ::Type{T}) where T = M.args