text:
module ft2d
using SpecialFunctions
export rectft,circft,ellipft,recip2real,recipvec2real,real2recip

function rectft(frax,fray,dnx,dny)
    #rectangle transforms to sinc
    return frax*fray*sinc.(frax*dnx).*sinc.(fray*dny)
end

function ellipft(frax,fray,dnx,dny)
    #argument of the bessel function
    radix=.5*sqrt.((dnx*frax).^2+(dny*fray).^2)
    result=besselj.(1,2*pi*radix)./radix
    #asymptotic value for going towards 0
    result[radix.==0].=pi
    #scale
    return result*frax*fray/4
end

function circft(fra,dnx,dny)
    return ellipft(fra,fra,dnx,dny)
end

function real2recip(dnx,dny,f,x,y)
    F=zeros(size(dnx))*1im
    for i=1:size(dnx,1)
        for j=1:size(dnx,2)
            F[i,j]=size(x,1)^-1*size(x,2)^-1*sum(f.*exp.(-1im*dnx[i,j]*x*2*pi-1im*dny[i,j]*y*2*pi))
        end
    end
    return F
end
function recip2real(dnx,dny,F,x,y)
    f=zeros(size(x))*1im
    for i=1:size(x,1)
        for j=1:size(y,2)
            f[i,j]=size(dnx,1)^-1*sum(F.*exp.(1im*x[i,j]*dnx*2*pi+1im*y[i,j]*dny*2*pi))
        end
    end
    return f
end

function recipvec2real(nx,ny,F,x,y)
    f=zeros(size(x))*1im
    for i=1:size(x,1)
        for j=1:size(y,2)
            f[i,j]=sum(F.*exp.(1im*x[i,j]*nx*2*pi+1im*y[i,j]*ny*2*pi))
        end
    end
    return f
end
end


text:
# 4.0 Deprecations
@deprecate space_neighbors nearby_ids
@deprecate node_neighbors nearby_positions
@deprecate get_node_contents ids_in_position
@deprecate get_node_agents agents_in_position
@deprecate pick_empty random_empty
@deprecate find_empty_nodes empty_positions
@deprecate has_empty_nodes has_empty_positions
@deprecate nodes positions

# 4.2 Deprecations
@deprecate fastest Schedulers.fastest
@deprecate by_id Schedulers.by_id
@deprecate as_added Schedulers.by_id
@deprecate random_activation Schedulers.randomly
@deprecate partial_activation Schedulers.partially
@deprecate property_activation Schedulers.by_property
@deprecate by_type Schedulers.by_type

@deprecate osm_random_road_position OSM.random_road_position
@deprecate osm_plan_route OSM.plan_route
@deprecate osm_map_coordinates OSM.map_coordinates
@deprecate osm_road_length OSM.road_length
@deprecate osm_random_route! OSM.random_route!
@deprecate osm_latlon OSM.latlon
@deprecate osm_intersection OSM.intersection
@deprecate osm_road OSM.road
@deprecate move_agent!(a, model::ABM{<:OpenStreetMapSpace}, distance) move_along_route!(a, model, distance)

ext: jl


text:
# Patches should allow using imported bindings in the body of the patch
@testset "imported binding in body" begin
    @test_throws UndefVarError Minute
    @test isdefined(Dates, :Minute)
    using Dates: Minute, Hour

    myminute(x::Integer) = Minute(x)

    # Patches should work when referencing bindings imported in the file where the patch
    # is created.
    patch = @patch myminute(x::Integer) = Minute(Hour(x))
    apply(patch) do
        @test (@mock myminute(5)) == Minute(300)
    end
end

# Patches should allow using <module>.<binding> syntax in the signature
@testset "qualified binding in signature" begin
    @test_throws UndefVarError AbstractCmd
    @test isdefined(Base, :AbstractCmd)

    patch = @patch read(cmd::Base.AbstractCmd, ::Type{String}) = "bar"
    apply(patch) do
        @test (@mock read(`foo`, String)) == "bar"
    end
end

# Patches should allow using imported bindings syntax in the signature
@testset "imported binding in signature" begin
    using Base: AbstractCmd

    patch = @patch read(cmd::AbstractCmd, ::Type{String}) = "bar"
    apply(patch) do
        @test (@mock read(`foo`, String)) == "bar"
    end
end


text:
module GeneralizedCRT

export crt
using Base.Threads
import Base.Threads.@spawn

const THRESHOLD1 = 7  # split longer input arrays to use binary instead of sequential algo
const THRESHOLD2 = 15 # split longer input arrays to different cpu threads

"""
    crt(a, b, p, q)

Given `0 <= a < p` and `0 <= b < q` with `a == b mod(gcd(p, q)`
return unique solution `x ≡ a mod p` and `x ≡ b mod q`
with `0 <= x < lcm(p, q)`.

Return `x, lcm(p, q)`.
"""
function crt(a::A, b::B, p::C, q::D) where {A,B,C,D}
    crt(promote(a, b, p, q)...)
end

function crt(a::T, b::T, p::T, q::T) where T
    c, u = gcdx(p, q)
    q1 = div(q, c)
    u = mod(u, q1)
    iszero(mod(b - a, c)) || throw(ArgumentError("remainders must have same rest when diveded by ($a, $b, gcd($p, $q))"))
    bac = div(b - a, c)
    x = mod(bac * u, q1) * p + a
    lcm = Base.checked_mul(p, q1)
    return x, lcm
end

function crt(ap::AbstractVector{<:Tuple{<:Integer,<:Integer}})
    crt([x[1] for x in ap], [x[2] for x in ap])
end
crt(t::Tuple{S,T}...) where {S<:Integer,T<:Integer} = crt(collect(t))

"""
    crt(a::Vector, m::Vector; nthreads, splitstreads, splitbinary)

Given `0 <= a[i] < p[i] for all i` with `a[i] == a[j] mod(gcd(p[i], p[j]) for i != j`,
return unique solution `x ≡ a[i] mod p[i] for all i`
with `0 <= x < lcm(p, q)`.

Optional arguments:
Use `nthreads=nthreads()` threads if `length(a) > splitthreads=($THRESHOLD2)`.
Use binary algorithm is `length(a) > splitbinary=($THRESHOLD1)`.

Return `x, lcm(p, q)`.
"""
function crt(a::AbstractVector{T}, m::AbstractVector{S}; nt=nthreads()-1, splitthreads=THRESHOLD2, splitbinary=THRESHOLD1) where {S,T}
    # split input array and execute in parallel threads
    n = length(a)
    n == length(m) || throw(ArgumentError("vectors of same size required"))
    if nt <= 1 || n <= splitthreads
        crts(a, m; splitbinary)
    else
        nt = max(nt, nthreads() - 1)
        R = promote_type(S, T)
        ni = max(n ÷ nt, 2)
        nt = (n + ni - 1) ÷ ni
        tasks = Vector{Any}(undef, nt)
        n2 = 0
        nth = nt
        for i = 1:nt
            n1 = n2 + 1
            n2 = n2 + (n - n2) ÷ nth
            nth -= 1
            av = view(a, n1:n2)
            mv = view(m, n1:n2)
            tasks[i] = @spawn crts($av, $mv; splitbinary=$splitbinary)
        end
        a0 = Vector{R}(undef, nt)
        m0 = Vector{R}(undef, nt)
        for i = 1:nt
            res = try
                fetch(tasks[i])
            catch ex
                ex isa TaskFailedException ? throw(ex.task.result) : rethrow()
            end
            a0[i] = first(res)
            m0[i] = last(res)
        end
        x, l = crt(a0, m0)
        x, l
    end
end

# split input array in halves a recombine output
function crts(a::AbstractVector{T}, m::AbstractVector{S}; splitbinary) where {S,T}
    n = length(a)
    n == length(m) || throw(ArgumentError("vectors of same size required"))
    if n <= max(splitbinary, 3)
        crtm(a, m)
    else
        n2 = (n + 1) ÷ 2
        xI, lcmI = crts(view(a,1:n2), view(m, 1:n2); splitbinary)
        xJ, lcmJ = crts(view(a, n2+1:n), view(m, n2+1:n); splitbinary)
        crt(xI, xJ, lcmI, lcmJ)
    end
end

# process input array sequentially
function crtm(a::AbstractVector{T}, m::AbstractVector{S}) where {S,T}
    n = length(a)
    # println("crtm($n, $(length(m)))")
    n == length(m) || throw(ArgumentError("vectors of same size required"))
    R = promote_type(S, T)
    n == 0 && return zero(R), one(R)
    xI, lcmI = promote(a[1], m[1])
    xI = _mod(xI, lcmI)
    for i = 2:n
        xI, lcmI = crt(xI, _mod(a[i], m[i]), lcmI, m[i])
    end
    xI, lcmI
end

_mod(a, b) = 0 <=a < b ? oftype(b, a) : mod(a, b)

end # module


text:
module Dashboards
import HTTP, JSON2
using MacroTools
include("ComponentPackages.jl")
include("ComponentMetas.jl")
include("Components.jl")
include("Front.jl")

import .ComponentPackages
import .Front
using .ComponentMetas
using .Components

export Dash, Component, Front, @use, <|, @callid_str, CallbackId, callback!,
 link_type!, make_handler, PreventUpdate, no_update

ComponentPackages.@reg_components()

@doc """
    module Dashboards

Julia backend for [Plotly Dash](https://github.com/plotly/dash)

# Examples
```julia
import HTTP
using Dashboards
app = Dash("Test", external_stylesheets=["https://codepen.io/chriddyp/pen/bWLwgP.css"]) do
    html_div() do
        dcc_input(id="graphTitle", value="Let's Dance!", type = "text"),
        html_div(id="outputID"),
        dcc_graph(id="graph",
            figure = (
                data = [(x = [1,2,3], y = [3,2,8], type="bar")],
                layout = Dict(:title => "Graph")
            )
        )

    end
end
callback!(app, callid"{graphTitle.type} graphTitle.value => outputID.children") do type, value
    "You've entered: '\$(value)' into a '\$(type)' input control"
end
callback!(app, callid"graphTitle.value => graph.figure") do value
    (
        data = [
            (x = [1,2,3], y = abs.(randn(3)), type="bar"),
            (x = [1,2,3], y = abs.(randn(3)), type="scatter", mode = "lines+markers", line = (width = 4,))
        ],
        layout = (title = value,)
    )
end
handle = make_handler(app, debug = true)
HTTP.serve(handle, HTTP.Sockets.localhost, 8080)
```

# Available components
$(ComponentPackages.component_doc_list())
""" Dashboards

const IdProp = Tuple{Symbol, Symbol}

struct PreventUpdate <: Exception

end

struct CallbackId
    state ::Vector{IdProp}
    input ::Vector{IdProp}
    output ::Vector{IdProp}
end

Base.convert(::Type{Vector{IdProp}}, v::IdProp) = [v]

CallbackId(;input ::Union{Vector{IdProp}, IdProp},  output ::Union{Vector{IdProp}, IdProp}, state ::Union{Vector{IdProp}, IdProp} = Vector{IdProp}()) = CallbackId(state, input, output)

struct Callback
    func ::Function
    id ::CallbackId
end

struct NoUpdate
end

no_update() = NoUpdate()

"""
    struct Dash <: Any

Representation of Dash application
"""
struct Dash
    name ::String
    layout ::Component
    callbacks ::Dict{Symbol, Callback}
    external_stylesheets ::Vector{String}
    url_base_pathname ::String
    assets_folder ::String
    callable_components ::Dict{Symbol, Component}
    type_links ::Dict{Symbol, Dict{Symbol, Type}}
    function Dash(name::String, layout::Component; external_stylesheets ::Vector{String} = Vector{String}(), url_base_pathname="/", assets_folder::String = "assets")
        new(name, layout, Dict{Symbol, Callback}(), external_stylesheets, url_base_pathname, assets_folder, Components.collect_with_ids(layout), Dict{Symbol, Dict{Symbol, Type}}())
    end
end

"""
    Dash(layout_maker::Function, name::String; external_stylesheets ::Vector{String} = Vector{String}(), url_base_pathname::String="/")::Dash

Construct a Dash app using callback for layout creation

# Arguments
- `layout_maker::Function` - function for layout creation. Must has signature ()::Component
- `name::String` - Dashboard name
- `external_stylesheets::Vector{String} = Vector{String}()` - vector of external css urls
- `url_base_pathname::String="/"` - base url path for dashboard, default "/"
- `assets_folder::String` - a path, relative to the current working directory,
for extra files to be used in the browser. Default `"assets"`

# Examples
```jldoctest
julia> app = Dash("Test") do
    html_div() do
        html_h1("Test Dashboard")
    end
end
```
"""
function Dash(layout_maker::Function, name::String;  external_stylesheets ::Vector{String} = Vector{String}(), url_base_pathname="/", assets_folder::String = "assets")
    Dash(name, layout_maker(), external_stylesheets=external_stylesheets, url_base_pathname=url_base_pathname, assets_folder = assets_folder)
end

function parse_props(s)
    function make_prop(part)
        m = match(r"^(?<id>[A-Za-z]+[\w\-\:\.]*)\.(?<prop>[A-Za-z]+[\w\-\:\.]*)$", strip(part))
        if isnothing(m)
            error("expected <id>.<property>[,<id>.<property>...] in $(part)")
        end
        return (Symbol(m[:id]), Symbol(m[:prop]))
    end

    props_parts = split(s, ",", keepempty = false)

    return map(props_parts) do part
        return make_prop(part)
    end
end

"""
    @callid_str"

Macro for crating Dash CallbackId.
Parse string in form "[{State1[, ...]}] Input1[, ...] => Output1[, ...]"

#Examples
```julia
    id1 = callid"{inputDiv.children} input.value => output1.value, output2.value"
```
"""
macro callid_str(s)
    rex = r"(\{(?<state>.*)\})?(?<input>.*)=>(?<output>.*)"ms
    m = match(rex, s)
    if isnothing(m)
        error("expected {state} input => output")
    end
    input = parse_props(strip(m[:input]))
    output = parse_props(strip(m[:output]))
    state = isnothing(m[:state]) ? Vector{IdProp}() : parse_props(strip(m[:state]))
    return CallbackId(state, input, output)
end

idprop_string(idprop::IdProp) = "$(idprop[1]).$(idprop[2])"

function check_idprop(app::Dash, id::IdProp)
    if !haskey(app.callable_components, id[1])
        error("The layout havn't component with id `$(id[1])]`")
    end
    if !is_prop_available(app.callable_components[id[1]], id[2])
        error("The component with id `$(id[1])` havn't property `$(id[2])``")
    end
end

function output_string(id::CallbackId)
    if length(id.output) == 1
        return idprop_string(id.output[1])
    end
    return ".." *
    join(map(idprop_string, id.output), "...") *
    ".."
end

"""
    callback!(func::Function, app::Dash, id::CallbackId)

Create a callback that updates the output by calling function `func`.

#Examples
```julia
app = Dash("Test") do
    html_div() do
        dcc_input(id="graphTitle", value="Let's Dance!", type = "text"),
        html_div(id="outputID"),
        html_div(id="outputID2")

    end
end
callback!(app, CallbackId(
    state = [(:graphTitle, :type)],
    input = [(:graphTitle, :value)],
    output = [(:outputID, :children), (:outputID2, :children)]
    )
    ) do stateType, inputValue
    return (stateType * "..." * inputValue, inputValue)
end
```
You can use macro `callid` string macro for make CallbackId :
```julia
callback!(app, callid"{graphTitle.type} graphTitle.value => outputID.children, outputID2.children") do stateType, inputValue

    return (stateType * "..." * inputValue, inputValue)
end
```
"""
function callback!(func::Function, app::Dash, id::CallbackId)
    for out in id.output
        if any(x->out in x.id.output, values(app.callbacks))
            error("output \"$(out)\" already registered")
        end
    end

    foreach(x->check_idprop(app,x), id.state)
    foreach(x->check_idprop(app,x), id.input)
    foreach(x->check_idprop(app,x), id.output)

    out_symbol = Symbol(output_string(id))

    push!(app.callbacks, out_symbol => Callback(func, id))
end

function callback_argument_type(app::Dash, id::AbstractString, prop::AbstractString)::Type
    id_sym = Symbol(id)
    prop_sym = Symbol(prop)
    if haskey(app.type_links, id_sym) && haskey(app.type_links[id_sym], prop_sym)
        return app.type_links[id_sym][prop_sym]
    end
    return Any
end

function push_link_type!(app::Dash, idprop::IdProp, t::Type)
    check_idprop(app, idprop)
    if !haskey(app.type_links, idprop[1])
        push!(app.type_links, idprop[1] => Dict{Symbol, Type}())
    end
    app.type_links[idprop[1]][idprop[2]] = t
end

function link_type!(app::Dash, idprop::AbstractString, t::Type)
    m = match(r"^((?<id>[A-Za-z]+[\w\-\:\.]*)|(?<any_id>\*))\.(?<prop>[A-Za-z]+[\w\-\:\.]*)$", strip(idprop))
    if isnothing(m)
        error("expected <id>.<property> or *.<property>")
    end

    if !isnothing(m[:id])
        push_link_type!(app, (Symbol(m[:id]), Symbol(m[:prop])), t)
    elseif !isnothing(m[:any_id])
        prop = Symbol(m[:prop])
        for (id, comp) in pairs(app.callable_components)
            if is_prop_available(comp, prop)
                push_link_type!(app, (id, prop), t)
            end
        end
    else
        @assert false "not reachable"
    end

end

function index_page(app::Dash; debug = false)
    metas = ["""<meta http-equiv="X-UA-Compatible" content="IE=edge">""", """<meta charset="UTF-8">"""]
    title = app.name
    css = join(map(app.external_stylesheets) do s
        """<link rel="stylesheet" href="$(s)">"""
    end, "\n"
    )

    app_entry = """

                Loading...

    """
    config = (
        url_base_pathname = "$(app.url_base_pathname)",
        requests_pathname_prefix = "$(app.url_base_pathname)",
        ui = true,
        props_check = debug,
        show_undo_redo = false
    )

    scripts = ComponentPackages.components_js_include(app.url_base_pathname, debug = debug)

    """<!DOCTYPE html>

            $(join(metas, ""))
            <title>$(title)</title>
            $(css)

            $(app_entry)

            $scripts

    """
end

function dependencies_json(app::Dash)
    id_prop_named(p::IdProp) = (id = p[1], property = p[2])
    result = map(values(app.callbacks)) do dep
        (inputs = id_prop_named.(dep.id.input),
        state = id_prop_named.(dep.id.state),
        output = output_string(dep.id)
        )
    end
    return JSON2.write(result)
end

function process_callback(app::Dash, body::String)
    params = JSON2.read(body)
    output = Symbol(params[:output])
    if !haskey(app.callbacks, output)
        return []
    end
    convert_values(inputs) = map(inputs) do x
        type = callback_argument_type(app, x.id, x.property)
        return haskey(x, :value) ? Front.from_dash(type, x.value) : nothing
    end
    args = []
    if haskey(params, :state)
        append!(args, convert_values(params.state))
        #append!(args, map(x->haskey(x, :value) ? x[:value] : nothing, params.state))
    end
    if haskey(params, :inputs)
        append!(args, convert_values(params.inputs))
        #append!(args, map(x->haskey(x, :value) ? x[:value] : nothing, params.inputs))
    end

    res = app.callbacks[output].func(args...)
    if length(app.callbacks[output].id.output) == 1
        if !(res isa NoUpdate)
            return Dict(
                :response => Dict(
                    :props => Dict(
                        Symbol(app.callbacks[output].id.output[1][2]) => Front.to_dash(res)
                    )
                )
            )
        else
            return Dict()
        end
    end
    response = Dict{Symbol, Any}()
    for (ind, out) in enumerate(app.callbacks[output].id.output)
        if !(res[ind] isa NoUpdate)
            push!(response,
            Symbol(out[1]) => Dict(
                Symbol(out[2]) => Front.to_dash(res[ind])
            )
            )
        end
    end
    return Dict(:response=>response, :multi=>true)

end

function process_assets(app::Dash, path)
    assets_path = "$(app.url_base_pathname)assets/"
    filename = joinpath(app.assets_folder, replace(path, assets_path=>""))
    try
        return HTTP.Response(200, [], body = read(filename))
    catch
        return HTTP.Response(404)
    end
end

"""
    make_handler(app::Dash; debug = false)

Make handler for routing Dash application in HTTP package

#Arguments
- `app::Dash` - Dash application
- `debug::Bool = false` - Enable/disable all the dev tools

#Examples
```jldoctest
julia> app = Dash("Test") do
    html_div() do
        html_h1("Test Dashboard")
    end
end
julia> handler = make_handler(app)
julia> HTTP.serve(handler, HTTP.Sockets.localhost, 8080)
```

"""
function make_handler(app::Dash; debug::Bool = false)
    function (req::HTTP.Request)
        uri = HTTP.URI(req.target)
        ComponentPackages.@register_js_sources(uri.path, app.url_base_pathname)
        if uri.path == "$(app.url_base_pathname)"
            return HTTP.Response(200, index_page(app, debug = debug))
        end
        if uri.path == "$(app.url_base_pathname)_dash-layout"
            return HTTP.Response(200, ["Content-Type" => "application/json"], body = JSON2.write(app.layout))
        end
        if uri.path == "$(app.url_base_pathname)_dash-dependencies"
            return HTTP.Response(200, ["Content-Type" => "application/json"], body = dependencies_json(app))
        end
        if startswith(uri.path, "$(app.url_base_pathname)assets/")
            return process_assets(app, uri.path)
        end
        if uri.path == "$(app.url_base_pathname)_dash-update-component" && req.method == "POST"
            try
                return HTTP.Response(200, ["Content-Type" => "application/json"],
                    body = JSON2.write(
                        process_callback(app, String(req.body))
                    )
                )
            catch e
                if isa(e,PreventUpdate)
                    return HTTP.Response(200)
                else
                    throw(e)
                end
            end
        end
        return HTTP.Response(404)
    end
end

end # module


text:
using Jags
using Base.Test

old = pwd()
ProjDir = Pkg.dir("Jags", "Examples", "Line1")
cd(ProjDir)

inits1 = [
  Dict("alpha" => 0,"beta" => 0,"tau" => 1),
  Dict("alpha" => 1,"beta" => 2,"tau" => 1),
  Dict("alpha" => 3,"beta" => 3,"tau" => 2),
  Dict("alpha" => 5,"beta" => 2,"tau" => 5)
]

function test(init::Array{Dict{String, Int64},1})
  res = map((x)->convert(Dict{String, Any}, x), init)
end

function test(init::Array{Dict{String, Float64},1})
  res = map((x)->convert(Dict{String, Any}, x), init)
end

function test(init::Array{Dict{String, Number},1})
  res = map((x)->convert(Dict{String, Any}, x), init)
end

function test(init::Array{Dict{String, Any},1})
  res = map((x)->convert(Dict{String, Any}, x), init)
end

line = "
model {
  for (i in 1:n) {
        mu[i] <- alpha + beta*(x[i] - x.bar);
        y[i]   ~ dnorm(mu[i],tau);
  }
  x.bar   <- mean(x[]);
  alpha    ~ dnorm(0.0,1.0E-4);
  beta     ~ dnorm(0.0,1.0E-4);
  tau      ~ dgamma(1.0E-3,1.0E-3);
  sigma   <- 1.0/sqrt(tau);
}
"

data = Dict{String, Any}()
data["x"] = [1, 2, 3, 4, 5]
data["y"] = [1, 3, 3, 3, 5]
data["n"] = 5

inits = test(inits1)

monitors = Dict{String, Bool}(
  "alpha" => true,
  "beta" => true,
  "tau" => true,
  "sigma" => true,
)

jagsmodel = Jagsmodel(name="line1", model=line, monitor=monitors,
  ncommands=1, nchains=4, adapt=1000, update=10000, thin=1,
  deviance=true, dic=true, popt=true, pdir=ProjDir);

println("\nJagsmodel that will be used:")
jagsmodel |> display
println("Input observed data dictionary:")
data |> display
println("\nInput initial values dictionary:")
inits |> display

cd(ProjDir)
isdir("tmp") &&
  rm("tmp", recursive=true);

cd(old)


text:
for Typ in (:Lmul, :Rmul)
    @eval begin
        struct $Typ{StyleA, StyleB, TypeA, TypeB}
            A::TypeA
            B::TypeB
        end

        $Typ(A::TypeA, B::TypeB) where {TypeA,TypeB} = $Typ{typeof(MemoryLayout(TypeA)),typeof(MemoryLayout(TypeB)),TypeA,TypeB}(A,B)

        BroadcastStyle(::Type{<:$Typ}) = ApplyBroadcastStyle()
        broadcastable(M::$Typ) = M

        eltype(::$Typ{<:Any,<:Any,A,B}) where {A,B} = promote_type(eltype(A), eltype(B))
        size(M::$Typ, p::Int) = size(M)[p]
        axes(M::$Typ, p::Int) = axes(M)[p]
        length(M::$Typ) = prod(size(M))
        size(M::$Typ) = map(length,axes(M))
        axes(M::$Typ) = (axes(M.A,1),axes(M.B,2))

        similar(M::$Typ, ::Type{T}, axes) where {T,N} = similar(Array{T}, axes)
        similar(M::$Typ, ::Type{T}) where T = similar(M, T, axes(M))
        similar(M::$Typ) = similar(M, eltype(M))
    end
end

const MatLmulVec{StyleA,StyleB} = Lmul{StyleA,StyleB,<:AbstractMatrix,<:AbstractVector}
const MatLmulMat{StyleA,StyleB} = Lmul{StyleA,StyleB,<:AbstractMatrix,<:AbstractMatrix}

const BlasMatLmulVec{StyleA,StyleB,T<:BlasFloat} = Lmul{StyleA,StyleB,<:AbstractMatrix{T},<:AbstractVector{T}}
const BlasMatLmulMat{StyleA,StyleB,T<:BlasFloat} = Lmul{StyleA,StyleB,<:AbstractMatrix{T},<:AbstractMatrix{T}}

const MatRmulMat{StyleA,StyleB} = Rmul{StyleA,StyleB,<:AbstractMatrix,<:AbstractMatrix}
const BlasMatRmulMat{StyleA,StyleB,T<:BlasFloat} = Rmul{StyleA,StyleB,<:AbstractMatrix{T},<:AbstractMatrix{T}}

####
# LMul materialize
####

axes(M::MatLmulVec) = (axes(M.A,1),)

lmul(A, B) = materialize(Lmul(A, B))

materialize(L::Lmul) = copy(instantiate(L))

copy(M::Lmul) = materialize!(Lmul(M.A,copy(M.B)))
copy(M::Rmul) = materialize!(Rmul(copy(M.A),M.B))

@inline function copyto!(dest::AbstractArray, M::Lmul)
    M.B ≡ dest || copyto!(dest, M.B)
    materialize!(Lmul(M.A,dest))
end

@inline function copyto!(dest::AbstractArray, M::Rmul)
    M.A ≡ dest || copyto!(dest, M.A)
    materialize!(Rmul(dest,M.B))
end

materialize!(M::Lmul) = lmul!(M.A,M.B)
materialize!(M::Rmul) = rmul!(M.A,M.B)

macro lazylmul(Typ)
    esc(quote
        LinearAlgebra.lmul!(A::$Typ, x::AbstractVector) = ArrayLayouts.materialize!(ArrayLayouts.Lmul(A,x))
        LinearAlgebra.lmul!(A::$Typ, x::AbstractMatrix) = ArrayLayouts.materialize!(ArrayLayouts.Lmul(A,x))
        LinearAlgebra.lmul!(A::$Typ, x::StridedVector) = ArrayLayouts.materialize!(ArrayLayouts.Lmul(A,x))
        LinearAlgebra.lmul!(A::$Typ, x::StridedMatrix) = ArrayLayouts.materialize!(ArrayLayouts.Lmul(A,x))
    end)
end


text:
module MatlabDecomp

using MATLAB: @mat_str
using TensorDecompositions: CANDECOMP

import SparseArrays

export
matlab_load_sptensor,
matlab_nncp_loaded_spt,
@mat_str,
put_adj_mat_mat

# Performance
# Runs NPM iterations at about 0.5 to 1 per minute

## using an attached matlab session
#
# This way is a bit clunky and slow compared to just using files.

mat"""
addpath ../matlab
addpath ../matlab/tensor_toolbox
"""

function matlab_load_sptensor(adj_mats)
    sparsefloat(am1) =
        convert(SparseArrays.SparseMatrixCSC{Float64,Int64}, am1)
    adj_mats = map(sparsefloat, adj_mats)
    mat"spt = sparse_matrix_list_to_sptensor($adj_mats)"
end

"Convert dict from matlab to CANDECOMP struct"
function _matlab_dict_to_cp(D, r)
    if r == 1
        lmbda = Array{Float64}(undef, r)
        lmbda .= D["lambda"]
        # Make it a 10x1 array, not a 10x0 array
        factors = map(permutedims ∘ permutedims, D["u"]|>Tuple)
        CANDECOMP(factors, lmbda)
    else
        CANDECOMP(D["u"]|>Tuple, D["lambda"])
    end
end

"""
    ncp(r, method="apg")

'apg' is the 'apg-tf' method from Xu and Yin 2013.

"""
function matlab_ncp_loaded_spt(r, method="apg")
    if method == "apg"
        mat"$D = ncp_apg(spt, $r, {});"
    end
    _matlab_dict_to_cp(D, r)
end

function relerror_loaded(D)
    D = Dict(("lambda" => D.lambdas, "u" => Any[D.factors...]))
    mat"relerror(spt, $D)"
end

matlab_nncp_loaded_spt(r) = matlab_ncp_loaded_spt(r, "apg")

## Trade .mat files

import MAT

function put_adj_mat_mat(platform, adj_mats=missing)
    filename = "../data/processed/$platform-adj-mats.mat"
    if isfile(filename)
        @info "Assuming pre-existing file is fine"
    else
        if ismissing(adj_mats)
            _, adj_mats = get_depversions_and_adj_mats(platform)
        end
        MAT.matwrite(filename, Dict("adj_mats"=>adj_mats,))
    end

    println("""
        Feed matlab something like this:

load('$filename')
X = sparse_matrix_list_to_sptensor(adj_mats);
r = 4;
D = ncp(X, r, {});
save(sprintf('data/processed/$platform-D%d.mat', r))
""")
end

matlab_load_npm_decomps(r) = begin
    platform = "NPM"
    D = MAT.matread("../data/processed/$platform-D$r.mat")["D"]
    # This might not work for r = 1; but who cares
    D = CANDECOMP(D["u"]|>Tuple, D["lambda"][:])
end

end


text:
function make_sh()
    open("submit.sh", "w") do io
        println(io, "#!/bin/bash")
        println(io, "#")
        println(io, "#SBATCH --job-name=prs_sim")
        println(io, "#")
        println(io, "#SBATCH --time=48:00:00")
        println(io, "#SBATCH --cpus-per-task=1")
        println(io, "#SBATCH --mem-per-cpu=20G")
        println(io, "#SBATCH --partition=owners")
        println(io, "")
        println(io, "#save job info on joblog:")
        println(io, "echo \"Job \$JOB_ID started on:   \" `hostname -s`")
        println(io, "echo \"Job \$JOB_ID started on:   \" `date `")
        println(io, "")
        println(io, "# load the job environment:")
        println(io, "module load julia")
        println(io, "module load python/3.9.0")
        # println(io, "module load R")
        # println(io, "module load gcc")
        # println(io, "module load cmake")
        println(io, "")
        println(io, "# run code")
        println(io, "echo 'running PRS genotype simulations'")
        println(io, "conda activate prs_sim")
        println(io, "python simulate_genotypes.py --sim 3 --threads 16")
        println(io, "")
        println(io, "#echo job info on joblog:")
        println(io, "echo \"Job \$JOB_ID ended on:   \" `hostname -s`")
        println(io, "echo \"Job \$JOB_ID ended on:   \" `date `")
        println(io, "#echo \" \"")
    end
end
make_sh()


text:
module AbstractOperations

export ∂x, ∂y, ∂z, @at, @unary, @binary, @multiary
export Δx, Δy, Δz, Ax, Ay, Az, volume
export Average, Integral, KernelFunctionOperation
export UnaryOperation, Derivative, BinaryOperation, MultiaryOperation, ConditionalOperation

using Base: @propagate_inbounds

import Adapt
using CUDA

using Oceananigans
using Oceananigans.Architectures
using Oceananigans.Grids
using Oceananigans.Operators
using Oceananigans.BoundaryConditions
using Oceananigans.Fields

using Oceananigans.Operators: interpolation_operator
using Oceananigans.Architectures: device
using Oceananigans: AbstractModel

import Oceananigans.Architectures: architecture
import Oceananigans.BoundaryConditions: fill_halo_regions!
import Oceananigans.Fields: compute_at!

#####
##### Basic functionality
#####

abstract type AbstractOperation{LX, LY, LZ, G, T} <: AbstractField{LX, LY, LZ, G, T, 3} end

const AF = AbstractField # used in unary_operations.jl, binary_operations.jl, etc

# We have no halos to fill
@inline fill_halo_regions!(::AbstractOperation, args...; kwargs...) = nothing

architecture(a::AbstractOperation) = architecture(a.grid)

# AbstractOperation macros add their associated functions to this list
const operators = Set()

"""
    at(loc, abstract_operation)

Return `abstract_operation` relocated to `loc`ation.
"""
at(loc, f) = f # fallback

include("grid_validation.jl")
include("grid_metrics.jl")
include("metric_field_reductions.jl")
include("unary_operations.jl")
include("binary_operations.jl")
include("multiary_operations.jl")
include("derivatives.jl")
include("kernel_function_operation.jl")
include("conditional_operations.jl")
include("computed_field.jl")
include("at.jl")
include("broadcasting_abstract_operations.jl")
include("show_abstract_operations.jl")

# Make some operators!

# Some operators:
import Base: sqrt, sin, cos, exp, tanh, -, +, /, ^, *

@unary sqrt sin cos exp tanh
@unary -
@unary +

@binary +
@binary -
@binary /
@binary ^

@multiary +

# For unknown reasons, the operator definition macros @binary and @multiary fail to work
# properly for :*. We thus manually define :* for fields.
import Base: *

eval(define_binary_operator(:*))
push!(operators, :*)
push!(binary_operators, :*)

eval(define_multiary_operator(:*))
push!(operators, :*)
push!(multiary_operators, :*)

end # module