text:
# helper test functions to check for SDD properties holding

function validate(sdd::Sdd)
    for n in sdd
      validate(n)
    end
    #TODO make one of these for structured decomposability
    @assert is_decomposable(sdd)
    is_unique(sdd, 5)
 end

function validate(n::SddNode)
   validate(GateType(n), n)
end

function validate(::⋁Gate, n::SddNode)
   size = num_children(n)
   primes = compile(false)
   for i = 1:size
      element = children(n)[i]
      # has alternating layers
      @test GateType(element) isa ⋀Gate
      for j = i+1:size
         other_element = children(n)[j]
         # is deterministic
         @test is_false(prime(element) & prime(other_element))
         # is compressed
         @test sub(element) !== sub(other_element)
      end
      primes = primes | prime(element)
   end
   # is exhaustive
   @test primes === compile(true)
   # cannot be trimmed to the sub
   @test size >= 1
   # cannot be trimmed to the prime
   if size == 2
      e1 = children(n)[1]
      e2 = children(n)[2]
      has_false_sub = (is_false(sub(e1)) || is_false(sub(e2)))
      has_true_sub = (is_true(sub(e1)) || is_true(sub(e2)))
      @test !(has_false_sub && has_true_sub)
   end
   @test NodeType(vtree(n)) isa Inner
end

function validate(::⋀Gate, n::SddNode)
   @test num_children(n) == 2
   @test !(GateType(prime(n)) isa ⋀Gate)
   @test !(GateType(sub(n)) isa ⋀Gate)
   # has no false prime
   @test !is_false(prime(n))
   @test NodeType(vtree(n)) isa Inner
   @test GateType(prime(n)) isa ConstantGate || descends_left_from(prime(n), n)
   @test GateType(sub(n)) isa ConstantGate || descends_right_from(sub(n), n)
end

function validate(::LiteralGate, l::SddNode)
   @test variable(l) == first(variables(vtree(l)))
end

function validate(::ConstantGate, ::SddNode)
   # nothing to check?
end

function is_unique(sdd::Sdd, k::Int)
   signatures = prob_equiv_signature(sdd, k)
   decision_nodes_by_signature = groupby(n -> signatures[n], ⋁_nodes(sdd))
   for (signature, nodes) in decision_nodes_by_signature
      if length(nodes) > 1
         println("Equivalent Nodes:")
         for node in nodes
            println("  - Node: $node Pr: $(sat_prob(node2dag(node)))")
         end
      end
      @test length(nodes) == 1
   end
end


text:
module Logistic

using LogicCircuits
using ...ProbabilisticCircuits

include("queries.jl")
include("parameters.jl")

# TODO structure learning

end

ext: jl


text:
function gpdps_potts!(u,v,data,params;
                      uref=0*u,
                      vref=0*v)
    σ,τ,ω,maxit     = params
    α,γ,f,isotropic = data

    fᵧ = t -> 2t^2/(2t^2 + γ)
    if isotropic
        Fᵧ = (u₁,u₂) -> sum(@. fᵧ(sqrt(u₁^2+u₂^2)))
        κ₁ = (u₁,u₂,v₁,v₂) -> 2(1-u₁*v₁-u₂*v₂)
        κ₂ = (u₁,u₂,v₁,v₂) -> 2(1-u₁*v₁-u₂*v₂)
    else
        Fᵧ = (u₁,u₂) -> sum(@. fᵧ(u₁)+fᵧ(u₂))
        κ₁ = (u₁,u₂,v₁,v₂) -> 2(1-u₁*v₁)
        κ₂ = (u₁,u₂,v₁,v₂) -> 2(1-u₂*v₂)
    end

    # preallocate
    uold = similar(u)
    ξ    = similar(u)
    v₁,v₂ = view(v,:,:,1), view(v,:,:,2)
    u₁,u₂ = similar(u), similar(u)
    z₁,z₂ = similar(u), similar(u)
    Kᵤ,Kᵥ = similar(u), similar(v)
    for it = 1:maxit
        # Kx = ∇ᵀ(κ(∇u⋅v)⋅v)
        u₁,u₂ = ∇!(u₁,u₂,u)
        @. z₁ = κ₁(u₁,u₂,v₁,v₂)*v₁
        @. z₂ = κ₂(u₁,u₂,v₁,v₂)*v₂
        Kᵤ = ∇ᵀ!(Kᵤ,z₁,z₂)

        # print output
        if mod(it,1000) == 0
            J   = norm(u.-f)^2/(2α) + Fᵧ(u₁,u₂)
            err = norm(u.-uref)^2+norm(v.-vref)^2
            @printf("it %6d: J = %1.6e, err = %1.6e\n", it, J, err)
            imshow(u', colormap=2)
        end

        # x = prox_τG(x-τKx)
        @. uold = u
        @. u = (u-τ*Kᵤ+τ/α*f)/(1+τ/α)
        @. ξ = u + ω*(u-uold)

        # Ky = κ(∇ξ⋅v)⋅∇ξ
        u₁,u₂ = ∇!(u₁,u₂,ξ)
        @. Kᵥ[:,:,1] = κ₁(u₁,u₂,v₁,v₂)*u₁
        @. Kᵥ[:,:,2] = κ₂(u₁,u₂,v₁,v₂)*u₂

        # y = prox_σF*(y+σKy)
        @. v = (v+σ*Kᵥ)/(1+γ*σ)
    end
end

function ∇!(u₁,u₂,u)
    @. @views begin
        u₁[1:(end-1), :] = u[2:end, :] - u[1:(end-1), :]
        u₁[end, :] = 0

        u₂[:, 1:(end-1)] = u[:, 2:end] - u[:, 1:(end-1)]
        u₂[:, end] = 0
    end
    return u₁, u₂
end

function ∇ᵀ!(v,v₁,v₂)
    @. @views begin
        v[2:(end-1), :] = v₁[1:(end-2), :] - v₁[2:(end-1), :]
        v[1, :] = -v₁[1, :]
        v[end, :] = v₁[end-1, :]

        v[:, 2:(end-1)] += v₂[:, 1:(end-2)] - v₂[:, 2:(end-1)]
        v[:, 1] += -v₂[:, 1]
        v[:, end] += v₂[:, end-1]
    end
    return v
end

function estimate_steplengths(f,α,γ,isotropic;expected_max_jump,γ_for_my)
    dyn_range = maximum(abs.(f))-minimum(abs.(f))
    expected_max_jump = 1.0
    γ_for_my = 10

    L = √8
    R_K = 2*L

    if isotropic
        mx = √2*expected_max_jump*dyn_range
    else
        mx = expected_max_jump*dyn_range
    end

    fmy = t -> 2t/(2t^2+γ_for_my)
    # fmy has its maximum in the interval [0, mx] either at mx or at √(γ_for_my/2)
    my = max(fmy(mx), fmy(√(γ_for_my/2)))

    # This estimate ignores the difficult-to-estimate ρ_y term, essentially
    # assuming ρ_y=0. To offset this, we set δ≪1.
    δ = 0.1
    μ = (1+δ)/2
    γ̃ = γ/2
    α̃ = α*2
    ξx = 1/α-1/α̃
    @assert(ξx>2*L*my^2)
    λx = 2*L^2*my^4/(ξx-2*L*my^2)
    λy = mx^2
    τ₀ = δ/λx
    τ₁ = 2*γ̃*α̃/(λy+√(λy^2+4*γ̃*α̃*(R_K^2/(1-μ)+2*λy/α̃)))
    τ  = min(τ₀, τ₁)
    σ = τ/(γ̃*α̃)
    ω = 1/(1+τ/α̃)

    return τ,σ,ω
end


text:
"""
    Abstract error type for errors in configuration files
"""
abstract type ConfigurationError <: Exception end

"""
    Specialized error type for errors in scanner configuration files
"""
struct ScannerConfigurationError <: ConfigurationError
  message::String
end

"""
    Specialized error type for errors in protocol configuration files
"""
struct ProtocolConfigurationError <: ConfigurationError
  message::String
end

"""
    Specialized error type for errors in sequence configuration files
"""
struct SequenceConfigurationError <: ConfigurationError
  message::String
end

"""
    Base.show(io::IO, ex::ConfigurationError)
Custom printing of `ConfigurationError` subtypes
"""
Base.show(io::IO, ex::ConfigurationError) = print(io, "$(typeof(ex)): $(ex.message)")

"""
    Abstract exception type for exceptions in devices
"""
abstract type AbstractDeviceException <: Exception end

"""
    General exception type for exceptions in devices
"""
struct DeviceException <: AbstractDeviceException
  message::String
  device::Union{Device, Nothing}
end

"""
    Base.show(io::IO, ex::DeviceException)
Custom printing of `DeviceException`
"""
function Base.show(io::IO, ex::DeviceException)
  if isnothing(ex.device)
    print(io, "$(typeof(ex)): $(ex.message)")
  else
    print(io, "$(typeof(ex)) in device with ID `$(ex.device.deviceID)`: $(ex.message)")
  end
end

abstract type ProtocolException end
struct CancelException <: ProtocolException end
struct IllegalStateException <: ProtocolException
  message::AbstractString
end


text:
using FITSIO

include("rbf_compress.jl")

const TILE_SIZE = 256
#const XCLUST = min(Int(round(TILE_SIZE / 8)), 32)#/16
#const YCLUST = min(Int(round(TILE_SIZE / 8)), 32)#/16
#const NCLUST = XCLUST * YCLUST

dir = "/home/chris/ダウンロード"
file = "ALMA01030862.fits"
#file = "ALMA01157085.fits"
#file = "Cygnus_sp46_vs-150_ve100_dv0.50_CN_Tmb.fits.gz"

fitspath = "/home/chris/NAO/NRO/SF/orion_12co_all_SF7.5arcsec_dV1.0kms.fits"
#fitspath = string(dir, "/", file)

if size(ARGS)[1] == 0
    println("usage: julia compress.jl <FITS filename>")
    #exit()
else
    fitspath = ARGS[1]
end

println("FITS file to compress: $(fitspath)")

f = FITS(fitspath)
N = ndims(f[1])
println("ndims: ", N, ", size: ", size(f[1]))

header = read_header(f[1])

width = 0
height = 0
depth = 1

if haskey(header, "NAXIS1")
    width = header["NAXIS1"]
end

if haskey(header, "NAXIS2")
    height = header["NAXIS2"]
end

if haskey(header, "NAXIS3")
    depth = header["NAXIS3"]
end

println("width: $(width), height: $(height), depth: $(depth)")

if depth < 1
    println("depth must be >= 1")
    close(f)
    exit()
end

if width < TILE_SIZE && height < TILE_SIZE
    println("the width and/or the height must be >= $(TILE_SIZE)")
    close(f)
    exit()
end

device, ctx, queue = cl.create_compute_context()
println(device)

#for frame = 1:depth
for frame = 1:1
    if N == 4
        data = read(f[1], :, :, frame, :);
    elseif N == 3
        data = read(f[1], :, :, frame);
    elseif N == 2
        data = read(f[1], :, :);
    end

    #println("HDU $(frame): ", size(data))
    sub = view(data, :, :, 1, 1)
    println("frame : ", frame, "\tdims: ", size(sub))

    ncols = Int(ceil(width / TILE_SIZE))
    nrows = Int(ceil(height / TILE_SIZE))

    println("nrows: $(nrows),  ncols: $(ncols)")

    for row in nrows:nrows
        for col in ncols:ncols
            x₁ = (col - 1) * TILE_SIZE
            x₂ = min(width, x₁ + TILE_SIZE)
            y₁ = (row - 1) * TILE_SIZE
            y₂ = min(height, y₁ + TILE_SIZE)
            println("processing row $(row) column $(col) :> x₁=$(x₁) x₂=$(x₂) y₁=$(y₁) y₂=$(y₂)")
            tile = view(sub, (x₁+1):x₂, (y₁+1):y₂)
            (N, padding, c1, c2, p0, p1, p2, w) = rbf_compress_tile(tile, device, ctx, queue)
        end
    end
end

close(f)


text:
struct QuadOverLinAtom <: AbstractExpr
    head::Symbol
    id_hash::UInt64
    children::Tuple{AbstractExpr,AbstractExpr}
    size::Tuple{Int,Int}

    function QuadOverLinAtom(x::AbstractExpr, y::AbstractExpr)
        if x.size[2] != 1 && y.size != (1, 1)
            error("quad over lin arguments must be a vector and a scalar")
        end
        children = (x, y)
        return new(:qol, hash(children), children, (1, 1))
    end
end

function sign(q::QuadOverLinAtom)
    return Positive()
end

function monotonicity(q::QuadOverLinAtom)
    return (sign(q.children[1]) * Nondecreasing(), Nonincreasing())
end

function curvature(q::QuadOverLinAtom)
    return ConvexVexity()
end

function evaluate(q::QuadOverLinAtom)
    x = evaluate(q.children[1])
    return x' * x / evaluate(q.children[2])
end

function conic_form!(q::QuadOverLinAtom, unique_conic_forms::UniqueConicForms)
    if !has_conic_form(unique_conic_forms, q)
        t = Variable()
        qol_objective = conic_form!(t, unique_conic_forms)
        x, y = q.children
        conic_form!(SOCConstraint(y + t, y - t, 2 * x), unique_conic_forms)
        conic_form!(y >= 0, unique_conic_forms)
        cache_conic_form!(unique_conic_forms, q, qol_objective)
    end
    return get_conic_form(unique_conic_forms, q)
end

quadoverlin(x::AbstractExpr, y::AbstractExpr) = QuadOverLinAtom(x, y)


text:
using DelimitedFiles
function to_phase(wfc_file, wfci_file, xDim, yDim)
    wfc = reshape(readdlm(wfc_file) + im*readdlm(wfci_file), xDim, yDim)
    hval = maximum(abs2.(wfc))
    return angle.(wfc).*(abs2.(wfc)/hval)
end

function to_density(wfc_file, wfci_file, xDim, yDim)
    wfc = reshape(readdlm(wfc_file) + im*readdlm(wfci_file), xDim, yDim)
    return abs2.(wfc)
end

wfc_file = "data_L10/wfc_L10"
wfci_file = "data_L10/wfci_L10"
xDim = 512
yDim = 512

writedlm("phase_L10.dat", to_phase(wfc_file, wfci_file, xDim, yDim))
writedlm("density_L10.dat", to_density(wfc_file, wfci_file, xDim, yDim))

max_stars_repo_name: leios/thesis

max_issues_repo_name: leios/thesis

max_forks_repo_name: leios/thesis


text:
# @todo change "Observation" to "Image"
struct Observation
    channels::Array{Float32, 4}
    function Observation(channel1, channel2, channel3)
        @assert size(channel1) == size(channel2) == size(channel3)
        channels = cat(channel1, channel2, channel3; dims=3)
        new(Flux.batch([channels]))
    end
    function Observation(images::Vector{Observation})
        channels = map(x->x.channels[:,:,:,1], images)
        new(Flux.batch(channels))
    end
end
Base.length(image::Observation) = size(image)[end]
Base.size(image::Observation) = size(image.channels)


text:
"""
    Dataset(::PyObject)

Loose wrapper around `xarray.Dataset`, mostly used for dispatch.

To create a `Dataset`, use [`convert_to_dataset`](@ref).
"""
struct Dataset
    o::PyObject

    function Dataset(o::PyObject)
        pyisinstance(
            o,
            xarray.Dataset,
        ) || throw(ArgumentError("$o is not an `xarray.Dataset`."))
        return new(o)
    end
end

@inline Dataset(data::Dataset) = data

@inline PyObject(data::Dataset) = getfield(data, :o)

Base.convert(::Type{Dataset}, o::PyObject) = Dataset(o)

Base.hash(data::Dataset) = hash(PyObject(data))

Base.propertynames(data::Dataset) = propertynames(PyObject(data))

function Base.getproperty(data::Dataset, name::Symbol)
    o = PyObject(data)
    name === :o && return o
    return getproperty(o, name)
end

function Base.show(io::IO, data::Dataset)
    out = pycall(pybuiltin("str"), String, data)
    out = replace(out, "<xarray.Dataset>" => "Dataset")
    print(io, out)
end

"""
    convert_to_dataset(obj; group = :posterior, kwargs...) -> Dataset

Convert a supported object to a `Dataset`. In most cases, this function calls
[`convert_to_inference_data`](@ref) and returns the corresponding `group`.
"""
function convert_to_dataset(obj; group = :posterior, kwargs...)
    group = Symbol(group)
    idata = convert_to_inference_data(obj; group = group, kwargs...)
    dataset = getproperty(idata, group)
    return dataset
end

convert_to_dataset(data::Dataset; kwargs...) = data

"""
    dict_to_dataset(data::Dict{String,Array}; kwargs...) -> Dataset

Convert a dictionary with data and keys as variable names to a [`Dataset`](@ref).

# Keywords
- `attrs::Dict{String,Any}`: Json serializable metadata to attach to the
    dataset, in addition to defaults.
- `library::String`: Name of library used for performing inference. Will be
    attached to the attrs metadata.
- `coords::Dict{String,Array}`: Coordinates for the dataset
- `dims::Dict{String,Vector{String}}`: Dimensions of each variable. The keys
    are variable names, values are vectors of coordinates.

# Examples

```@example
using ArviZ
ArviZ.dict_to_dataset(Dict("x" => randn(4, 100), "y" => randn(4, 100)))
```
"""
function dict_to_dataset(data; library = nothing, attrs = nothing, kwargs...)
    if library !== nothing
        ldict = Dict("inference_library" => string(library))
        attrs = (attrs === nothing ? ldict : merge(attrs, ldict))
    end
    return arviz.dict_to_dataset(data; attrs = attrs, kwargs...)
end

"""
    dataset_to_dict(ds::Dataset) -> Tuple{Dict{String,Array},NamedTuple}

Convert a `Dataset` to a dictionary of `Array`s. The function also
returns keyword arguments to [`dict_to_dataset`](@ref).
"""
function dataset_to_dict(ds::Dataset)
    ds_dict = ds.to_dict()
    data_vars = ds_dict["data_vars"]
    attrs = ds_dict["attrs"]

    coords = ds_dict["coords"]
    delete!(coords, "chain")
    delete!(coords, "draw")
    coords = Dict(k => v["data"] for (k, v) in coords)

    data = Dict{String,Array}()
    dims = Dict{String,Vector{String}}()
    for (k, v) in data_vars
        data[k] = v["data"]
        dim = v["dims"][3:end]
        if !isempty(dim)
            dims[k] = [dim...]
        end
    end

    return data, (attrs = attrs, coords = coords, dims = dims)
end


text:
# Autogenerated wrapper script for libnabo_jll for aarch64-apple-darwin
export libnabo

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("libnabo")
JLLWrappers.@declare_library_product(libnabo, "@rpath/libnabo.1.dylib")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libnabo,
        "lib/libnabo.1.0.7.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module OpenTelemetryExporterPrometheus

export PrometheusExporter

using OpenTelemetrySDK
using HTTP
using Sockets

mutable struct PrometheusExporter <: AbstractExporter
    server::Sockets.TCPServer
    provider::Union{MeterProvider,Nothing}
    taskref::Ref{Task}
    function PrometheusExporter(; host = "127.0.0.1", port = 9966, kw...)
        server = Sockets.listen(Sockets.InetAddr(parse(IPAddr, host), port))
        exporter = new(server, nothing, Ref{Task}())
        exporter.taskref[] = @async HTTP.listen(host, port; server = server, kw...) do http
            HTTP.setstatus(http, 200)
            HTTP.setheader(http, "Content-Type" => "text/plain")

            if isnothing(exporter.provider)
                write(http, "MeterProvider is not set yet!!!")
            else
                for ins in exporter.provider.async_instruments
                    ins()
                end
                text_based_format(http, exporter.provider)
            end
            return
        end
        exporter
    end
end

function (r::MetricReader{<:MeterProvider,<:PrometheusExporter})()
    if isnothing(r.exporter.provider)
        r.exporter.provider = r.provider
    else
        @info "The prometheus exporter is already properly set!"
    end
end

# TODO: support exemplars
function text_based_format(io, provider::MeterProvider)
    for m in values(provider.metrics)
        write(io, "# HELP $(m.name) $(m.description)\n")
        write(io, "# TYPE $(m.name) $(prometheus_type(m.aggregation))\n")
        for (attrs, point) in m
            if point.value isa OpenTelemetrySDK.HistogramValue
                val = point.value
                for (i, c) in enumerate(Iterators.accumulate(+, val.counts))
                    if i == length(val.counts)
                        write(io, "$(m.name)_bucket{le=\"+Inf\"} $c\n")
                        write(io, "$(m.name)_count $c\n")
                    else
                        write(io, "$(m.name)_bucket{le=\"$(point.boundaries[i])\"} $c\n")
                    end
                end
                # ???
                if !isnothing(point.sum)
                    write(io, "$(m.name)_sum $(point.sum)\n")
                end
            else
                write(io, "$(m.name){")
                # TODO: escape
                join(io, ("$k=\"$v\"" for (k, v) in pairs(attrs)), ",")
                write(io, "} $(point.value) $(point.time_unix_nano ÷ 10^6) \n")
            end
        end
        write(io, "\n")
    end
end

prometheus_type(::SumAgg) = "counter"
prometheus_type(::LastValueAgg) = "gauge"
prometheus_type(::HistogramAgg) = "histogram"

end # module


text:
using Test
import UUIDs
import Pluto: PlutoRunner, Notebook, WorkspaceManager, Cell, ServerSession, ClientSession, update_run!
using Pluto.WorkspaceManager: poll

@testset "Logging" begin
    🍭 = ServerSession()
    🍭.options.evaluation.workspace_use_distributed = true

    fakeclient = ClientSession(:fake, nothing)
    🍭.connected_clients[fakeclient.id] = fakeclient

    @testset "Logging respects maxlog" begin
        @testset "Single log" begin
            notebook = Notebook(Cell.([
                """
                for i in 1:10
                    @info "logging" i maxlog=2
                end
                """,
            ]))

            update_run!(🍭, notebook, notebook.cells)
            @test notebook.cells[begin] |> noerror

            @test poll(5, 1/60) do
                length(notebook.cells[begin].logs) == 2
            end

            # Check that maxlog doesn't occur in the message
            @test all(notebook.cells[begin].logs) do log
                all(log["kwargs"]) do kwarg
                    kwarg[1] != "maxlog"
                end
            end

            WorkspaceManager.unmake_workspace((🍭, notebook))
        end

        @testset "Multiple log" begin
            notebook = Notebook(Cell.([
                """
                for i in 1:10
                    @info "logging" i maxlog=2
                    @info "logging more" maxlog = 4
                    @info "even more logging"
                end
                """,
            ]))

            update_run!(🍭, notebook, notebook.cells)
            @test notebook.cells[begin] |> noerror

            # Wait until all 16 logs are in
            @test poll(5, 1/60) do
                length(notebook.cells[begin].logs) == 16
            end

            # Get the ids of the three logs and their counts. We are
            # assuming that the logs are ordered same as in the loop.
            ids = unique(getindex.(notebook.cells[begin].logs, "id"))
            counts = [count(log -> log["id"] == id, notebook.cells[begin].logs) for id in ids]
            @test counts == [2, 4, 10]

            # Check that maxlog doesn't occur in the messages
            @test all(notebook.cells[begin].logs) do log
                all(log["kwargs"]) do kwarg
                    kwarg[1] != "maxlog"
                end
            end

            WorkspaceManager.unmake_workspace((🍭, notebook))
        end
    end
end


text:
include("polygon.jl")
include("special.jl")
include("sector.jl")


text:
using Compat, ClassicalLaminateTheory

forcesAndMoments =
 @compat Dict(
    :Nx => 300,
    :Ny => 100,
    :Nxy => 200,
    :Mx => 5,
    :My => 3,
    :Mxy => 1)

materialProperties =@compat Dict(
  :E1 => 20.0*10^6,
  :E2 => 1.5*10^6,
  :G12 => 1.0*10^6,
  :ν12 => .34,
  :alpha1 => 0.2e-6,
  :alpha2 => 20e-6,
  :beta1 => 0.2e-4,
  :beta2 => 20e-4,
  :thickness => 0.0075)

laminateProperties =@compat Dict(
  :nplies => 5,
  :repeats => 1,
  :symmetric => false,
  :orientation => [0, 45, 90, 45, -45],
  :deltaTemp => -270.0,
  :deltaMoisture => 1
)

createLaminate!(laminateProperties, materialProperties, forcesAndMoments)

fm = [
157.552,
-42.4476,
210.226,
4.84661,
3.15339,
0.846607
]

@assert(round(laminateProperties[:fTotal], 3) == round(fm, 3))
println("Tests for CLT2 passed!")


text:
using DataFrames, SQLite
export ContinuousSpace, index!

#######################################################################################
# Continuous space structure
#######################################################################################
struct ContinuousSpace{F, E} <: AbstractSpace
  D::Int
  update_vel!::F
  periodic::Bool
  extend::E
  metric::String
  db::SQLite.DB
  insertq::SQLite.Stmt
  searchq::SQLite.Stmt
  deleteq::SQLite.Stmt
  updateq::SQLite.Stmt
  searchqNoId::SQLite.Stmt
end

const COORDS = 'a':'z' # letters representing coordinates in database

"""
    ContinuousSpace(D::Int [, update_vel!]; periodic::Bool = true, extend = nothing, metric = "cityblock")
Create a `ContinuousSpace` of dimensionality `D`.
In this case, your agent positions (field `pos`) should be of type `NTuple{D, F}`
where `F <: AbstractFloat`.
In addition, the agent type must have a third field `vel::NTuple{D, F}` representing
the agent's velocity.

The optional argument `update_vel!` is a **function**, `update_vel!(agent, model)` that updates
the agent's velocities **before** the agent has been moved, see [`move_agent!`](@ref).
You can of course change the agents velocities
during the agent interaction, the `update_vel!` functionality targets arbitrary forces.
By default no update is done this way.

## Keywords

* `periodic = true` : whether continuous space is periodic or not
* `extend::NTuple{D} = ones` : Extend of space. The `d` dimension starts at 0
  and ends at `extend[d]`. If `periodic = true`, this is also when
  periodicity occurs. If `periodic ≠ true`, `extend` is only used at plotting.
"""
function ContinuousSpace(D::Int, update_vel! = defvel;
  periodic = true, extend = Tuple(1.0 for i in 1:D), metric = "cityblock")

  # TODO: implement using different metrics in space_neighbors
  @assert metric ∈ ("cityblock", "euclidean")

  # TODO: allow extend to be useful even without periodicity: agents bounce of walls then
  # (improve to do this `move_agent!`)

  db, q, q2, q3, q4, q5 = prepare_database(D)
  ContinuousSpace(D, update_vel!, periodic, extend, metric, db, q, q2, q3, q4, q5)
end

# Deprecate Space constructor
@deprecate Space(D::Int, update_vel!::Function) ContinuousSpace(D::Int, update_vel!::Function)

function prepare_database(D)
  db = SQLite.DB()
  dimexpression = join("$x REAL, " for x in COORDS[1:D])
  stmt = "CREATE TABLE tab ("*dimexpression*"id INTEGER PRIMARY KEY)"
  DBInterface.execute(db, stmt)
  insertedxpression = join("$x, " for x in COORDS[1:D])
  qmarks = join("?, " for _ in 1:D)
  insertstmt = "INSERT INTO tab ($(insertedxpression)id) VALUES ($(qmarks)?)"
  q = DBInterface.prepare(db, insertstmt)
  searchexpr = join("$x BETWEEN ? AND ? AND " for x in COORDS[1:D])
  searchq = "SELECT id FROM tab WHERE $(searchexpr)id != ?"
  q2 = DBInterface.prepare(db, searchq)
  deleteq = "DELETE FROM tab WHERE id = ?"
  q3 = DBInterface.prepare(db, deleteq)
  updateexpr = strip(join("$x = ?," for x in COORDS[1:D]), ',')
  updateq = "UPDATE tab SET $updateexpr WHERE id = ?"
  q4 = DBInterface.prepare(db, updateq)
  searchexpr2 = join("$x BETWEEN ? AND ? AND " for x in COORDS[1:D])
  searchq2 = "SELECT id FROM tab WHERE $(searchexpr2)"[1:end-4]
  q5 = DBInterface.prepare(db, searchq2)
  return db, q, q2, q3, q4, q5
end

defvel(a, m) = nothing

# TODO: re-check this at the end.
function Base.show(io::IO, abm::ContinuousSpace)
    s = "$(abm.D)-dimensional $(abm.periodic ? "periodic " : "")ContinuousSpace"
    abm.update_vel! ≠ defvel && (s *= " with velocity updates")
    print(io, s)
end

#######################################################################################
# SQLite database functions
#######################################################################################
"Collect IDs from an SQLite.Query where IDs are stored in `colname`"
function collect_ids(q::SQLite.Query, colname=:id)
  output = Union{Int, Missing}[]
  for row in q
    push!(output, row[colname])
  end
  return output
end

# TODO: index! needs a better name
"""
    index!(model)
Index the database underlying the `ContinuousSpace` of the model.

This can drastically improve performance for retrieving data, but adding new
data can become slower because after each addition, index needs to be called again.

Lack of index won't be noticed for small databases. Only use it when you have
many agents and not many additions of agents.
"""
function index!(model)
  D = model.space.D
  expr = join("$x," for x in COORDS[1:D])
  DBInterface.execute(model.space.db, "CREATE INDEX pos ON tab ($(expr)id)")
end

#######################################################################################
# Extention of Agents.jl Model-Space interaction API
#######################################################################################
# central, low level function that is always called by all others!
function add_agent_pos!(agent::A, model::ABM{A, <: ContinuousSpace}) where {A<:AbstractAgent}
  DBInterface.execute(model.space.insertq, (agent.pos..., agent.id))
  model[agent.id] = agent
  return agent
end

function randompos(space::ContinuousSpace)
  pos = Tuple(rand(space.D))
  space.extend ≠ nothing && (pos = pos .* space.extend)
  return pos
end

function add_agent!(agent::A, model::ABM{A, <: ContinuousSpace}) where {A<:AbstractAgent}
  agent.pos = randompos(model.space)
  add_agent_pos!(agent, model)
end

function add_agent!(agent::A, pos, model::ABM{A, <: ContinuousSpace}) where {A<:AbstractAgent}
  agent.pos = pos
  add_agent_pos!(agent, model)
end

# versions that create the agents
function add_agent!(model::ABM{A, <: ContinuousSpace}, args...) where {A<:AbstractAgent}
  add_agent!(randompos(model.space), model, args...)
end

function add_agent!(pos::Tuple, model::ABM{A, <: ContinuousSpace}, args...) where {A<:AbstractAgent}
  id = nextid(model)
  agent = A(id, pos, args...)
  add_agent_pos!(agent, model)
end

"""
    move_agent!(agent::A, model::ABM{A, ContinuousSpace}, dt = 1.0)
Propagate the agent forwards one step according to its velocity,
_after_ updating the agent's velocity (see [`ContinuousSpace`](@ref)).

For this continuous space version of `move_agent!`, the "evolution algorithm"
is a trivial Euler scheme with `dt` the step size, i.e. the agent position is updated
as `agent.pos += agent.vel * dt`.
"""
function move_agent!(agent::A, model::ABM{A, S, F, P}, dt = 1.0) where {A<:AbstractAgent, S <: ContinuousSpace, F, P}
  model.space.update_vel!(agent, model)
  agent.pos = agent.pos .+ dt .* agent.vel
  if model.space.periodic
    agent.pos = mod.(agent.pos, model.space.extend)
  end
  DBInterface.execute(model.space.updateq, (agent.pos..., agent.id))
  return agent.pos
end

function kill_agent!(agent::AbstractAgent, model::ABM{A, S}) where {A, S<:ContinuousSpace}
  DBInterface.execute(model.space.deleteq, (agent.id,))
  delete!(model.agents, agent.id)
  return model
end

function genocide!(model::ABM{A, S}, n::Int) where {A, S<:ContinuousSpace}
  ids = strip(join("$id," for id in keys(model.agents) if id > n), ',')
  DBInterface.execute(model.space.db, "DELETE FROM tab WHERE id IN ($ids)")
  for id in keys(model.agents)
    id > n && (delete!(model.agents, id))
  end
  return model
end

function genocide!(model::ABM{A, S}, f::Function) where {A, S<:ContinuousSpace}
  ids = strip(join("$(agent.id)," for agent in values(model.agents) if f(agent)), ',')
  DBInterface.execute(model.space.db, "DELETE FROM tab WHERE id IN ($ids)")
  for agent in values(model.agents)
    f(agent) && (delete!(model.agents, agent.id))
  end
  return model
end

function genocide!(model::ABM{A, S}) where {A, S<:ContinuousSpace}
  DBInterface.execute(model.space.db, "DELETE FROM tab")
  for agent in model.agents
    delete!(model.agents, agent.id)
  end
end

# TODO: at the moment this function doesn't check the metric and uses only cityblock
# we can easily adjust to arbitrary metric by doing a final check
# filter!(...) where the filtering function checks distances w.r.t. `r`.
"""
    space_neighbors(pos::Tuple, model::ABM, r::Real)
Return IDs of all agents within radius `r` from a particular position `pos` for any space.
"""
function space_neighbors(pos::Tuple, model, r::Real)
  left = pos .- r
  right = pos .+ r
  res = interlace(left, right)
  collect_ids(DBInterface.execute(model.space.searchqNoId, res))
end

"""
    space_neighbors(agent::AbstractAgent, model::ABM, r::Real)
Return neighbours of a particular agent, within radius `r` for any space.
"""
function space_neighbors(agent::A, model::ABM{A, <:ContinuousSpace}, r::Real) where {A<:AbstractAgent}
  left = agent.pos .- r
  right = agent.pos .+ r
  res = interlace(left, right)
  collect_ids(DBInterface.execute(model.space.searchq, (res...,agent.id)))
end

@generated function interlace(left::NTuple{D}, right::NTuple{D}) where {D}
  a = [[:(left[$i]), :(right[$i])] for i=1:D]
  b = vcat(a...)
  quote
    tuple($(b...))
  end
end

#######################################################################################
# Continuous space exclusive
#######################################################################################
export nearest_neighbor, elastic_collision!, interacting_pairs

"""
    nearest_neighbor(agent, model, r) → nearest
Return the agent that has the closest distance to given `agent`, according to the
space's metric.
Valid only in continuous space.
"""
function nearest_neighbor(agent, model, r)
  n = space_neighbors(agent, model, r)
  length(n) == 0 && return nothing
  d, j = Inf, 0
  for i in 1:length(n)
    @inbounds dnew = sqrt(sum(abs2.(agent.pos .- model[n[i]].pos)))
    _, j = findmin(d)
    if dnew < d
      d, j = dnew, i
    end
  end
  return model[n[j]]
end

using LinearAlgebra

"""
    elastic_collision!(a, b, f = nothing)
Resolve a (hypothetical) elastic collision between the two agents `a, b`.
They are assumed to be disks of equal size touching tangentially.
Their velocities (field `vel`) are adjusted for an elastic collision happening between them.
This function works only for two dimensions.
Notice that collision only happens if both disks face each other, to avoid
collision-after-collision.

If `f` is a `Symbol`, then the agent property `f`, e.g. `:mass`, is taken as a mass
to weight the two agents for the collision. By default no weighting happens.

One of the two agents can have infinite "mass", and then acts as an immovable object
that specularly reflects the other agent. In this case of course momentum is not
conserved, but kinetic energy is still conserved.
"""
function elastic_collision!(a, b, f = nothing)
  # Do elastic collision according to
  # https://en.wikipedia.org/wiki/Elastic_collision#Two-dimensional_collision_with_two_moving_objects
  v1, v2, x1, x2 = a.vel, b.vel, a.pos, b.pos
  length(v1) != 2 && error("This function works only for two dimensions.")
  r1 = x1 .- x2; r2 = x2 .- x1
  m1, m2 = f == nothing ? (1.0, 1.0) : (getfield(a, f), getfield(b, f))
  # mass weights
  m1 == m2 == Inf && return false
  if m1 == Inf
    @assert v1 == (0, 0) "An agent with ∞ mass cannot have nonzero velocity"
    dot(r1, v2) ≤ 0 && return false
    v1 = ntuple(x -> zero(eltype(v1)), length(v1))
    f1, f2 = 0.0, 2.0
  elseif m2 == Inf
    @assert v2 == (0, 0) "An agent with ∞ mass cannot have nonzero velocity"
    dot(r2, v1) ≤ 0 && return false
    v2 = ntuple(x -> zero(eltype(v1)), length(v1))
    f1, f2 = 2.0, 0.0
  else
    # Check if disks face each other, to avoid double collisions
    !(dot(r2, v1) > 0 && dot(r2, v1) > 0) && return false
    f1 = (2m2/(m1+m2))
    f2 = (2m1/(m1+m2))
  end
  ken = norm(v1)^2 + norm(v2)^2
  dx = a.pos .- b.pos
  dv = a.vel .- b.vel
  n = norm(dx)^2
  n == 0 && return false # do nothing if they are at the same position
  a.vel = v1 .- f1 .* ( dot(v1 .- v2, r1) / n ) .* (r1)
  b.vel = v2 .- f2 .* ( dot(v2 .- v1, r2) / n ) .* (r2)
  return true
end

"""
    interacting_pairs(model, r)
Return an iterator that yields pairs of agents `(a1, a2)` that are closest
neighbors to each other, within some interaction radius `r`.

This function is usefully combined with `model_step!`, when one wants to perform
some pairwise interaction across all pairs of closest agents once
(and does not want to trigger the event twice, both with `a1` and with `a2`, which
is unavoidable when using `agent_step!`).

Internally uses [`nearest_neighbor`](@ref).
"""
function interacting_pairs(model, r)
  pairs = Tuple{Int, Int}[]
  #TODO: This can be optimized further I assume
  for id in keys(model.agents)
    # Skip already checked agents
    any(isequal(id), p[2] for p in pairs) && continue
    a1 = model[id]
    a2 = nearest_neighbor(a1, model, r)
    a2 ≠ nothing && push!(pairs, (id, a2.id))
  end
  return PairIterator(pairs, model.agents)
end

struct PairIterator{A}
  pairs::Vector{Tuple{Int, Int}}
  agents::Dict{Int, A}
end

Base.length(iter::PairIterator) = length(iter.pairs)
function Base.iterate(iter::PairIterator, i = 1)
  i > length(iter) && return nothing
  p = iter.pairs[i]
  id1, id2 = p
  return (iter.agents[id1], iter.agents[id2]), i+1
end


text:
using FiniteVolumeMethod

grid = TreeGrid([0., 0.], periodic = true)


text:
const required_fields = [
    "byr",  # Birth Year
    "ecl",  # Eye Color
    "eyr",  # Expiration Year
    "hcl",  # Hair Color
    "hgt",  # Height
    "iyr",  # Issue Year
    "pid"   # Passport ID
]

function new_passport(fields)

    p = Dict{String,Any}()
    for field in vcat(split.(fields)...)
       key, value = split(field,":")
       p[key] = value
    end
    return p

end

function read_passports(filename)

   passports = []
   open(filename) do f

       fields  = []
       for line in eachline(f)

           if isempty(line)
               push!(passports, new_passport(fields))
               fields = []
           else
               push!(fields, line)
           end
       end

       # end of the file, last passport
       push!(passports, new_passport(fields))

   end

   return passports

end

function valid1( p )

    valid_fields = collect(filter(k->(k!="cid"), keys(p)))

    if length(valid_fields) == 7
        return all( sort(valid_fields) .== required_fields)
    else
        return false
    end

end

const eye_colors = ["amb","blu","brn","gry","grn","hzl","oth"]

function valid2( p ) # not finish...

    if valid1(p)

        1920 <= parse(Int, p["byr"]) <= 2002 || (return false)
        2010 <= parse(Int, p["iyr"]) <= 2020 || (return false)
        2020 <= parse(Int, p["eyr"]) <= 2030 || (return false)
        hgt = p["hgt"]
        if endswith(hgt, "cm")
            150 <= parse(Int,hgt[1:end-2]) <= 193 || (return false)
        elseif endswith(hgt, "in")
            59 <= parse(Int,hgt[1:end-2]) <= 76 || (return false)
        else
            return false
        end
        startswith(p["hcl"], "#") || (return false)
        tryparse(Int, p["hcl"][2:end], base=16) == nothing && (return false)
        p["ecl"] in eye_colors || (return false)
        length(p["pid"]) == 9 || (return false)
        tryparse(Int, p["pid"]) == nothing && (return false)

        return true

    else
        return false
    end

end

passports = read_passports( "input04.txt")
println(count(valid1.(passports)))
println(count(valid2.(passports)))

max_stars_repo_name: pnavaro/AdventOfCode

max_issues_repo_name: pnavaro/AdventOfCode

max_forks_repo_name: pnavaro/AdventOfCode


text:
module SpikeSynchrony

include("SPIKEdistance.jl")
include("vanRossum.jl")

end # module


text:
module test_h3_api_inspection

using Test
using H3.API # h3GetResolution h3GetBaseCell stringToH3 h3ToString h3IsValid h3IsResClassIII h3IsPentagon

@test h3GetResolution(0x085283473fffffff) == 5

@test h3GetBaseCell(0x085283473fffffff) == 20

@test stringToH3("85283473fffffff") == 0x085283473fffffff

@test h3ToString(0x085283473fffffff) == "85283473fffffff"

@test h3IsValid(0x85283473fffffff)
@test !h3IsValid(0x5004295803a8800)

@test h3IsResClassIII(0x085283473fffffff)

@test !h3IsPentagon(0x085283473fffffff)

end # module test_h3_api_inspection


text:
@testset "AlgAssRelOrdIdl" begin
  Qx, x = FlintQQ["x"]
  f = x^2 - 10x - 8
  K, a = number_field(f, "a")
  KG = group_algebra(K, GrpAbFinGen([ 2 ]))

  @testset "Arithmetic" begin
    O = any_order(KG)
    I = 2*O
    J = 4*O

    @test I + J == I
    @test I*J == 8*O
    @test intersect(I, J) == J
    @test I^2 == J
    @test I^fmpz(2) == J

    @test norm(I) == 4*base_ring(O)
  end

  @testset "Locally free basis" begin
    Ky, y = K["y"]
    OK = maximal_order(K)
    g = y^2 + 1
    L, b = number_field(g, "b") # Gal(L/K) == C_2
    OL = maximal_order(L)

    LtoKG, KGtoL = Hecke._find_isomorphism(L, KG)
    basisOL = Vector{elem_type(KG)}()
    for i = 1:degree(L)
      for b in basis(pseudo_basis(OL, copy = false)[i][2])
        push!(basisOL, LtoKG(b*pseudo_basis(OL, copy = false)[i][1]))
      end
    end
    d = lcm([ denominator(b) for b in basisOL ])
    OKG = Order(KG, basis(KG))
    I = Hecke.ideal_from_lattice_gens(KG, OKG, [ d*b for b in basisOL ])

    p = prime_decomposition(OK, 3)[1][1]
    g = Hecke.locally_free_basis(I, p)
    OKGg = OKG*g
    mat_I = Hecke.coprime_bases(OKG, I, p)[4]
    mat_OKGg = Hecke.coprime_bases(OKG, OKGg, p)[4]
    t = det(mat_I*inv(mat_OKGg))
    @test valuation(t, p) == 0

    p = prime_decomposition(OK, 2)[1][1]
    @test_throws ErrorException  Hecke.locally_free_basis(I, p)
  end

end


text:
# Run this in the REPL, not Juno!!!
using PkgBenchmark
bresult = benchmarkpkg("DynamicalSystemsBase"; retune = true)


text:
module cookstress_examples
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsDeforLinear.AlgoDeforLinearModule
using FinEtools.MeshExportModule
import LinearAlgebra: cholesky

function cookstress()
    println("Cook membrane problem,  plane stress."        )
    t0 = time()

    E = 1.0;
    nu = 1.0/3;
    width = 48.0; height = 44.0; thickness  = 1.0;
    free_height  = 16.0;
    Mid_edge  = [48.0, 52.0];# Location of tracked  deflection
    magn = 1.0/free_height;# Magnitude of applied load
    convutip = 23.97;
    n = 32;#*int(round(sqrt(170.)/2.)); # number of elements per side
    tolerance = minimum([width, height])/n/1000.;#Geometrical tolerance

    fens, fes = T3block(width, height,  n,  n)

    # Reshape into a trapezoidal panel
    for i=1:count(fens)
        fens.xyz[i, 2]=fens.xyz[i, 2]+(fens.xyz[i, 1]/width)*(height -fens.xyz[i, 2]/height*(height-free_height));
    end

    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz, 1), 2)) # displacement field

    l1 = selectnode(fens; box=[0, 0, -Inf,  Inf],  inflate = tolerance)
    setebc!(u, l1, 1, val=0.0)
    setebc!(u, l1, 2, val=0.0)
    applyebc!(u)
    numberdofs!(u)

    @time boundaryfes =  meshboundary(fes);
    Toplist = selectelem(fens, boundaryfes,  box= [width,  width,  -Inf,  Inf ],  inflate=  tolerance);
    el1femm =  FEMMBase(IntegDomain(subset(boundaryfes, Toplist),  GaussRule(1, 2)))
    fi = ForceIntensity([0.0, +magn]);
    F2 = distribloads(el1femm,  geom,  u,  fi,  2);

    MR = DeforModelRed2DStress
    material = MatDeforElastIso(MR,  0.0, E, nu, 0.0)

    femm = FEMMDeforLinear(MR, IntegDomain(fes,  TriRule(1)),  material)

    K = stiffness(femm,  geom,  u)
    K = cholesky(K)
    U=  K\(F2)
    scattersysvec!(u, U[:])

    nl = selectnode(fens,  box=[Mid_edge[1], Mid_edge[1], Mid_edge[2], Mid_edge[2]], inflate=tolerance);
    theutip = zeros(FFlt, 1, 2)
    gathervalues_asmat!(u, theutip, nl);
    println("$(time()-t0) [s];  displacement =$(theutip[2]) as compared to converged $convutip")

    File =  "a.vtk"
    vtkexportmesh(File,  fes.conn,  geom.values+u.values,
    FinEtools.MeshExportModule.T3; vectors=[("u", u.values)])

    true

end # cookstress

function cookstress_algo()
    E = 1.0;
    nu = 1.0/3;
    width = 48.0; height = 44.0; thickness  = 1.0;
    free_height  = 16.0;
    Mid_edge  = [48.0, 52.0];# Location of tracked  deflection
    magn = 1.0/free_height;# Magnitude of applied load
    convutip = 23.97;
    n = 20;#*int(round(sqrt(170.)/2.)); # number of elements per side
    tolerance = minimum([width, height])/n/1000.;#Geometrical tolerance

    fens,fes = T3block(width, height, n, n)

    # Reshape into a trapezoidal panel
    for i = 1:count(fens)
        fens.xyz[i,2] = fens.xyz[i,2]+(fens.xyz[i,1]/width)*(height -fens.xyz[i,2]/height*(height-free_height));
    end

    # Clamped edge of the membrane
    l1 = selectnode(fens; box=[0.,0.,-Inf, Inf], inflate = tolerance)
    ess1 = FDataDict("displacement"=>  0.0, "component"=> 1, "node_list"=>l1)
    ess2 = FDataDict("displacement"=>  0.0, "component"=> 2, "node_list"=>l1)

    # Traction on the opposite edge
    boundaryfes =  meshboundary(fes);
    Toplist  = selectelem(fens, boundaryfes, box= [width, width, -Inf, Inf ], inflate=  tolerance);
    el1femm = FEMMBase(IntegDomain(subset(boundaryfes, Toplist), GaussRule(1, 2)))
    flux1 = FDataDict("traction_vector"=>[0.0,+magn], "femm"=>el1femm)

    # Make the region
    MR = DeforModelRed2DStress
    material = MatDeforElastIso(MR,  0.0, E, nu, 0.0)
    region1 = FDataDict("femm"=>FEMMDeforLinear(MR,
    IntegDomain(fes, TriRule(1)), material))

    modeldata = FDataDict("fens"=>fens,    "regions"=>[region1],    "essential_bcs"=>[ess1, ess2],    "traction_bcs"=>[flux1])

    # Call the solver
    modeldata = AlgoDeforLinearModule.linearstatics(modeldata)

    u = modeldata["u"]
    geom = modeldata["geom"]

    # Extract the solution
    nl = selectnode(fens, box=[Mid_edge[1],Mid_edge[1],Mid_edge[2],Mid_edge[2]],
    inflate=tolerance);
    theutip = u.values[nl,:]
    println("displacement =$(theutip[2]) as compared to converged $convutip")

    modeldata["postprocessing"] = FDataDict("file"=>"cookstress", "quantity"=>:Cauchy, "component"=>:xy)
    modeldata = AlgoDeforLinearModule.exportstress(modeldata)

    true

end # cookstress_algo

function cookstress_algo_export()
    E = 1.0;
    nu = 1.0/3;
    width = 48.0; height = 44.0; thickness  = 1.0;
    free_height  = 16.0;
    Mid_edge  = [48.0, 52.0];# Location of tracked  deflection
    magn = 1.0/(free_height*thickness);# Density of applied load
    convutip = 23.97;
    n = 30;# number of elements per side
    tolerance = minimum([width, height])/n/1000.;#Geometrical tolerance

    fens,fes = T3block(width, height, n, n)

    # Reshape into a trapezoidal panel
    for i = 1:count(fens)
        fens.xyz[i,2] = fens.xyz[i,2]+(fens.xyz[i,1]/width)*(height -fens.xyz[i,2]/height*(height-free_height));
    end

    # Clamped edge of the membrane
    l1 = selectnode(fens; box=[0.,0.,-Inf, Inf], inflate = tolerance)
    ess1 = FDataDict("displacement"=>  0.0, "component"=> 1, "node_list"=>l1)
    ess2 = FDataDict("displacement"=>  0.0, "component"=> 2, "node_list"=>l1)

    # Traction on the opposite edge
    boundaryfes =  meshboundary(fes);
    Toplist  = selectelem(fens, boundaryfes, box= [width, width, -Inf, Inf ], inflate=  tolerance);
    el1femm = FEMMBase(IntegDomain(subset(boundaryfes, Toplist), GaussRule(1, 2), thickness))
    flux1 = FDataDict("traction_vector"=>[0.0,+magn], "femm"=>el1femm )

    # Make the region
    MR = DeforModelRed2DStress
    material = MatDeforElastIso(MR,  0.0, E, nu, 0.0)
    region1 = FDataDict("femm"=>FEMMDeforLinear(MR,
    IntegDomain(fes, TriRule(1), thickness), material))

    modeldata = FDataDict("fens"=>fens, "regions"=>[region1], "essential_bcs"=>[ess1, ess2], "traction_bcs"=>[flux1])

    # Call the solver
    modeldata = AlgoDeforLinearModule.linearstatics(modeldata)

    u = modeldata["u"]
    geom = modeldata["geom"]

    # Extract the solution
    nl = selectnode(fens, box=[Mid_edge[1],Mid_edge[1],Mid_edge[2],Mid_edge[2]],
    inflate=tolerance);
    theutip = u.values[nl,:]
    println("displacement =$(theutip[2]) as compared to converged $convutip")

    modeldata["postprocessing"] = FDataDict("file"=>"cookstress-ew",
    "quantity"=>:Cauchy, "component"=>:xy)
    modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)
    fld = modeldata["postprocessing"]["exported"][1]["field"]
    println("range of Cauchy_xy = $((minimum(fld.values), maximum(fld.values)))")
    File = modeldata["postprocessing"]["exported"][1]["file"]
    @async run(`"paraview.exe" $File`)

    modeldata["postprocessing"] = FDataDict("file"=>"cookstress-ew-vm",
    "quantity"=>:vm, "component"=>1)
    modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)
    fld = modeldata["postprocessing"]["exported"][1]["field"]
    println("range of vm = $((minimum(fld.values), maximum(fld.values)))")
    File = modeldata["postprocessing"]["exported"][1]["file"]
    @async run(`"paraview.exe" $File`)

    modeldata["postprocessing"] = FDataDict("file"=>"cookstress-ew-pressure",
    "quantity"=>:pressure, "component"=>1)
    modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)
    fld = modeldata["postprocessing"]["exported"][1]["field"]
    println("range of pressure = $((minimum(fld.values), maximum(fld.values)))")
    File = modeldata["postprocessing"]["exported"][1]["file"]
    @async run(`"paraview.exe" $File`)

    modeldata["postprocessing"] = FDataDict("file"=>"cookstress-ew-princ1",
    "quantity"=>:princCauchy, "component"=>1)
    modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)
    fld = modeldata["postprocessing"]["exported"][1]["field"]
    println("range of princCauchy Max = $((minimum(fld.values), maximum(fld.values)))")
    File = modeldata["postprocessing"]["exported"][1]["file"]
    @async run(`"paraview.exe" $File`)

    modeldata["postprocessing"] = FDataDict("file"=>"cookstress-ew-princ2",
    "quantity"=>:princCauchy, "component"=>2)
    modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)
    fld = modeldata["postprocessing"]["exported"][1]["field"]
    println("range of princCauchy Min = $((minimum(fld.values), maximum(fld.values)))")
    File = modeldata["postprocessing"]["exported"][1]["file"]
    @async run(`"paraview.exe" $File`)

    AE = AbaqusExporter("Cookstress_algo_stress");
    HEADING(AE, "Cook trapezoidal panel, plane stress");
    COMMENT(AE, "Converged free mid-edge displacement = 23.97");
    PART(AE, "part1");
    END_PART(AE);
    ASSEMBLY(AE, "ASSEM1");
    INSTANCE(AE, "INSTNC1", "PART1");
    NODE(AE, fens.xyz);
    COMMENT(AE, "We are assuming three node triangles in plane-stress");
    COMMENT(AE, "CPS3 are pretty poor-accuracy elements, but here we don't care about it.");
    @assert nodesperelem(modeldata["regions"][1]["femm"].integdomain.fes) == 3
    ELEMENT(AE, "CPS3", "AllElements", connasarray(modeldata["regions"][1]["femm"].integdomain.fes))
    NSET_NSET(AE, "clamped", modeldata["essential_bcs"][1]["node_list"])
    ORIENTATION(AE, "GlobalOrientation", vec([1. 0 0]), vec([0 1. 0]));
    SOLID_SECTION(AE, "elasticity", "GlobalOrientation", "AllElements", thickness);
    END_INSTANCE(AE);
    END_ASSEMBLY(AE);
    MATERIAL(AE, "elasticity")
    ELASTIC(AE, E, nu)
    STEP_PERTURBATION_STATIC(AE)
    BOUNDARY(AE, "ASSEM1.INSTNC1.clamped", 1)
    BOUNDARY(AE, "ASSEM1.INSTNC1.clamped", 2)
    bfes = modeldata["traction_bcs"][1]["femm"].integdomain.fes
    COMMENT(AE, "Concentrated loads: we are assuming that the elements on the boundary");
    COMMENT(AE, "have two nodes each and also that they are the same length.");
    COMMENT(AE, "Then the concentrated loads below will be correctly lumped.");
    nl = connectednodes(bfes)
    F = zeros(count(modeldata["fens"]))
    for ix = 1:count(bfes)
        for jx = 1:2
            F[bfes.conn[ix][jx]] += 1.0/n/2/thickness
        end
    end
    for ixxxx = 1:length(F)
        if F[ixxxx] != 0.0
            CLOAD(AE, "ASSEM1.INSTNC1.$(ixxxx)", 2, F[ixxxx])
        end
    end
    END_STEP(AE)
    close(AE)

    true

end # cookstress_algo_export

function cookstress_algo_export_ortho()
    E = 1.0;
    nu = 1.0/3;
    width = 48.0; height = 44.0; thickness  = 1.0;
    free_height  = 16.0;
    Mid_edge  = [48.0, 52.0];# Location of tracked  deflection
    magn = 1.0/(free_height*thickness);# Density of applied load
    convutip = 23.97;
    n = 30;# number of elements per side
    tolerance = minimum([width, height])/n/1000.;#Geometrical tolerance

    fens,fes = T3block(width, height, n, n)

    # Reshape into a trapezoidal panel
    for i = 1:count(fens)
        fens.xyz[i,2] = fens.xyz[i,2]+(fens.xyz[i,1]/width)*(height -fens.xyz[i,2]/height*(height-free_height));
    end

    # Clamped edge of the membrane
    l1 = selectnode(fens; box=[0.,0.,-Inf, Inf], inflate = tolerance)
    ess1 = FDataDict("displacement"=>  0.0, "component"=> 1, "node_list"=>l1)
    ess2 = FDataDict("displacement"=>  0.0, "component"=> 2, "node_list"=>l1)

    # Traction on the opposite edge
    boundaryfes =  meshboundary(fes);
    Toplist  = selectelem(fens, boundaryfes, box= [width, width, -Inf, Inf ], inflate=  tolerance);
    el1femm = FEMMBase(IntegDomain(subset(boundaryfes, Toplist), GaussRule(1, 2), thickness))
    flux1 = FDataDict("traction_vector"=>[0.0,+magn],
    "femm"=>el1femm
    )

    # Make the region
    MR = DeforModelRed2DStress
    material = MatDeforElastOrtho(MR,  0.0, E, nu, 0.0)
    region1 = FDataDict("femm"=>FEMMDeforLinear(MR,
    IntegDomain(fes, TriRule(1), thickness), material))

    modeldata = FDataDict("fens"=>fens,
    "regions"=>[region1],
    "essential_bcs"=>[ess1, ess2],
    "traction_bcs"=>[flux1]
    )

    # Call the solver
    modeldata = AlgoDeforLinearModule.linearstatics(modeldata)

    u = modeldata["u"]
    geom = modeldata["geom"]

    # Extract the solution
    nl = selectnode(fens, box=[Mid_edge[1],Mid_edge[1],Mid_edge[2],Mid_edge[2]],
    inflate=tolerance);
    theutip = u.values[nl,:]
    println("displacement =$(theutip[2]) as compared to converged $convutip")

    modeldata["postprocessing"] = FDataDict("file"=>"cookstress-ew",
    "quantity"=>:Cauchy, "component"=>:xy)
    modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)
    fld = modeldata["postprocessing"]["exported"][1]["field"]
    println("range of Cauchy_xy = $((minimum(fld.values), maximum(fld.values)))")
    File = modeldata["postprocessing"]["exported"][1]["file"]
    @async run(`"paraview.exe" $File`)

    modeldata["postprocessing"] = FDataDict("file"=>"cookstress-ew-vm",
    "quantity"=>:vm, "component"=>1)
    modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)
    fld = modeldata["postprocessing"]["exported"][1]["field"]
    println("range of vm = $((minimum(fld.values), maximum(fld.values)))")
    File = modeldata["postprocessing"]["exported"][1]["file"]
    @async run(`"paraview.exe" $File`)

    modeldata["postprocessing"] = FDataDict("file"=>"cookstress-ew-pressure",
    "quantity"=>:pressure, "component"=>1)
    modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)
    fld = modeldata["postprocessing"]["exported"][1]["field"]
    println("range of pressure = $((minimum(fld.values), maximum(fld.values)))")
    File = modeldata["postprocessing"]["exported"][1]["file"]
    @async run(`"paraview.exe" $File`)

    modeldata["postprocessing"] = FDataDict("file"=>"cookstress-ew-princ1",
    "quantity"=>:princCauchy, "component"=>1)
    modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)
    fld = modeldata["postprocessing"]["exported"][1]["field"]
    println("range of princCauchy Max = $((minimum(fld.values), maximum(fld.values)))")
    File = modeldata["postprocessing"]["exported"][1]["file"]
    @async run(`"paraview.exe" $File`)

    modeldata["postprocessing"] = FDataDict("file"=>"cookstress-ew-princ2",
    "quantity"=>:princCauchy, "component"=>2)
    modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)
    fld = modeldata["postprocessing"]["exported"][1]["field"]
    println("range of princCauchy Min = $((minimum(fld.values), maximum(fld.values)))")
    File = modeldata["postprocessing"]["exported"][1]["file"]
    @async run(`"paraview.exe" $File`)

    AE = AbaqusExporter("Cookstress_algo_stress");
    HEADING(AE, "Cook trapezoidal panel, plane stress");
    COMMENT(AE, "Converged free mid-edge displacement = 23.97");
    PART(AE, "part1");
    END_PART(AE);
    ASSEMBLY(AE, "ASSEM1");
    INSTANCE(AE, "INSTNC1", "PART1");
    NODE(AE, fens.xyz);
    COMMENT(AE, "We are assuming three node triangles in plane-stress");
    COMMENT(AE, "CPS3 are pretty poor-accuracy elements, but here we don't care about it.");
    @assert nodesperelem(modeldata["regions"][1]["femm"].integdomain.fes) == 3
    ELEMENT(AE, "CPS3", "AllElements", connasarray(modeldata["regions"][1]["femm"].integdomain.fes))
    NSET_NSET(AE, "clamped", modeldata["essential_bcs"][1]["node_list"])
    ORIENTATION(AE, "GlobalOrientation", vec([1. 0 0]), vec([0 1. 0]));
    SOLID_SECTION(AE, "elasticity", "GlobalOrientation", "AllElements", thickness);
    END_INSTANCE(AE);
    END_ASSEMBLY(AE);
    MATERIAL(AE, "elasticity")
    ELASTIC(AE, E, nu)
    STEP_PERTURBATION_STATIC(AE)
    BOUNDARY(AE, "ASSEM1.INSTNC1.clamped", 1)
    BOUNDARY(AE, "ASSEM1.INSTNC1.clamped", 2)
    bfes = modeldata["traction_bcs"][1]["femm"].integdomain.fes
    COMMENT(AE, "Concentrated loads: we are assuming that the elements on the boundary");
    COMMENT(AE, "have two nodes each and also that they are the same length.");
    COMMENT(AE, "Then the concentrated loads below will be correctly lumped.");
    nl = connectednodes(bfes)
    F = zeros(count(modeldata["fens"]))
    for ix = 1:count(bfes)
        for jx = 1:2
            F[bfes.conn[ix][jx]] += 1.0/n/2/thickness
        end
    end
    for ixxxx = 1:length(F)
        if F[ixxxx] != 0.0
            CLOAD(AE, "ASSEM1.INSTNC1.$(ixxxx)", 2, F[ixxxx])
        end
    end
    END_STEP(AE)
    close(AE)

    true

end # cookstress_algo_export_ortho

function cookstress_t6_algo()
    println("Cook plane stress, with quadratic triangles.")
    E = 1.0;
    nu = 1.0/3;
    width = 48.0; height = 44.0; thickness  = 1.0;
    free_height  = 16.0;
    Mid_edge  = [48.0, 52.0];# Location of tracked  deflection
    magn = 1.0/free_height;# Magnitude of applied load
    convutip = 23.97;
    n = 10;#*int(round(sqrt(170.)/2.)); # number of elements per side
    tolerance = minimum([width, height])/n/1000.;#Geometrical tolerance

    fens,fes = T6block(width, height, n, n)

    # Reshape into a trapezoidal panel
    for i = 1:count(fens)
        fens.xyz[i,2] = fens.xyz[i,2]+(fens.xyz[i,1]/width)*(height -fens.xyz[i,2]/height*(height-free_height));
    end

    # Clamped edge of the membrane
    l1 = selectnode(fens; box=[0.,0.,-Inf, Inf], inflate = tolerance)
    ess1 = FDataDict("displacement"=>  0.0, "component"=> 1, "node_list"=>l1)
    ess2 = FDataDict("displacement"=>  0.0, "component"=> 2, "node_list"=>l1)

    # Traction on the opposite edge
    boundaryfes =  meshboundary(fes);
    Toplist  = selectelem(fens, boundaryfes, box= [width, width, -Inf, Inf ], inflate=  tolerance);
    el1femm = FEMMBase(IntegDomain(subset(boundaryfes, Toplist), GaussRule(1, 3)))
    flux1 = FDataDict("traction_vector"=>[0.0,+magn],
    "femm"=>el1femm
    )

    # Make the region
    MR = DeforModelRed2DStress
    material = MatDeforElastIso(MR,  0.0, E, nu, 0.0)
    region1 = FDataDict("femm"=>FEMMDeforLinear(MR,
    IntegDomain(fes, TriRule(3)), material))

    modeldata = FDataDict("fens"=>fens,
    "regions"=>[region1],
    "essential_bcs"=>[ess1, ess2],
    "traction_bcs"=>[flux1]
    )

    # Call the solver
    modeldata = AlgoDeforLinearModule.linearstatics(modeldata)

    u = modeldata["u"]
    geom = modeldata["geom"]

    # Extract the solution
    nl = selectnode(fens, box=[Mid_edge[1],Mid_edge[1],Mid_edge[2],Mid_edge[2]],
    inflate=tolerance);
    theutip = u.values[nl,:]
    println("displacement =$(theutip[2]) as compared to converged $convutip")

    modeldata["postprocessing"] = FDataDict("file"=>"cookstress",
    "quantity"=>:Cauchy, "component"=>:xy)
    modeldata = AlgoDeforLinearModule.exportstress(modeldata)
    File = modeldata["postprocessing"]["exported"][1]["file"]
    @async run(`"paraview.exe" $File`)
    fld = modeldata["postprocessing"]["exported"][1]["field"]
    println("$(minimum(fld.values)) $(maximum(fld.values))")
    true

end # cookstress_t6_algo

function cookstress_t6_ortho2iso_algo()
    println("Cook plane stress, with quadratic triangles. With orthotropic  material model.")
    E = 1.0;
    nu = 1.0/3;
    E1 = E2 = E3 = E;
    nu12 = nu13 = nu23 = nu;
    G12 = G13 = G23 = E/2.0/(1+nu);
    width = 48.0; height = 44.0; thickness  = 1.0;
    free_height  = 16.0;
    Mid_edge  = [48.0, 52.0];# Location of tracked  deflection
    magn = 1.0/free_height;# Magnitude of applied load
    convutip = 23.97;
    n = 10;#*int(round(sqrt(170.)/2.)); # number of elements per side
    tolerance = minimum([width, height])/n/1000.;#Geometrical tolerance

    fens,fes = T6block(width, height, n, n)

    # Reshape into a trapezoidal panel
    for i = 1:count(fens)
        fens.xyz[i,2] = fens.xyz[i,2]+(fens.xyz[i,1]/width)*(height -fens.xyz[i,2]/height*(height-free_height));
    end

    # Clamped edge of the membrane
    l1 = selectnode(fens; box=[0.,0.,-Inf, Inf], inflate = tolerance)
    ess1 = FDataDict("displacement"=>  0.0, "component"=> 1, "node_list"=>l1)
    ess2 = FDataDict("displacement"=>  0.0, "component"=> 2, "node_list"=>l1)

    # Traction on the opposite edge
    boundaryfes =  meshboundary(fes);
    Toplist  = selectelem(fens, boundaryfes, box= [width, width, -Inf, Inf ], inflate=  tolerance);
    el1femm = FEMMBase(IntegDomain(subset(boundaryfes, Toplist), GaussRule(1, 3)))
    flux1 = FDataDict("traction_vector"=>[0.0,+magn], "femm"=>el1femm)

    # Make the region
    MR = DeforModelRed2DStress
    # This material model is orthotropic,  but the input parameters correspond to an
    # isotropiic material  model..
    material=MatDeforElastOrtho(MR, E1,E2,E3,nu12,nu13,nu23,G12,G13,G23)
    region1 = FDataDict("femm"=>FEMMDeforLinear(MR, IntegDomain(fes, TriRule(3)), material))

    modeldata = FDataDict("fens"=>fens, "regions"=>[region1], "essential_bcs"=>[ess1, ess2], "traction_bcs"=>[flux1])

    # Call the solver
    modeldata = AlgoDeforLinearModule.linearstatics(modeldata)

    u = modeldata["u"]
    geom = modeldata["geom"]

    # Extract the solution
    nl = selectnode(fens, box=[Mid_edge[1],Mid_edge[1],Mid_edge[2],Mid_edge[2]],    inflate=tolerance);
    theutip = u.values[nl,:]
    println("displacement =$(theutip[2]) as compared to converged $convutip")

    modeldata["postprocessing"] = FDataDict("file"=>"cookstress",
    "quantity"=>:Cauchy, "component"=>:xy)
    modeldata = AlgoDeforLinearModule.exportstress(modeldata)
    File = modeldata["postprocessing"]["exported"][1]["file"]
    @async run(`"paraview.exe" $File`)
    fld = modeldata["postprocessing"]["exported"][1]["field"]
    println("$(minimum(fld.values)) $(maximum(fld.values))")
    true

end # cookstress_t6_ortho2iso_algo

function allrun()
    println("#####################################################")
    println("# cookstress ")
    cookstress()
    println("#####################################################")
    println("# cookstress_algo ")
    cookstress_algo()
    println("#####################################################")
    println("# cookstress_algo_export ")
    cookstress_algo_export()
    println("#####################################################")
    println("# cookstress_algo_export_ortho ")
    cookstress_algo_export_ortho()
    println("#####################################################")
    println("# cookstress_t6_algo ")
    cookstress_t6_algo()
    println("#####################################################")
    println("# cookstress_t6_ortho2iso_algo ")
    cookstress_t6_ortho2iso_algo()
    return true
end # function allrun

end # module cookstress_examples


text:
int_rules_4_3_1_1 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 4.3.1.1*(a+b*tan)^n =#
    @apply_utils Antiderivative((~(b') * tan(~(c') + ~(d') * ~x)) ^ ~n, ~x) => (~b * (~b * tan(~c + ~d * ~x)) ^ (~n - 1)) / (~d * (~n - 1)) - (~b) ^ 2 * Antiderivative((~b * tan(~c + ~d * ~x)) ^ (~n - 2), ~x) <-- FreeQ([~b, ~c, ~d], ~x) && GtQ(~n, 1)
    @apply_utils Antiderivative((~(b') * tan(~(c') + ~(d') * ~x)) ^ ~n, ~x) => (~b * tan(~c + ~d * ~x)) ^ (~n + 1) / (~b * ~d * (~n + 1)) - (1 / (~b) ^ 2) * Antiderivative((~b * tan(~c + ~d * ~x)) ^ (~n + 2), ~x) <-- FreeQ([~b, ~c, ~d], ~x) && LtQ(~n, -1)
    @apply_utils Antiderivative(tan(~(c') + ~(d') * ~x), ~x) => -(log(RemoveContent(cos(~c + ~d * ~x), ~x))) / ~d <-- FreeQ([~c, ~d], ~x)
    #= Antiderivative(1/tan((~c')+(~d')*(~x)),~x) :=  log(RemoveContent(sin(c+d*x),x))/d <-- FreeQ([c,d],x) =#
    @apply_utils Antiderivative((~(b') * tan(~(c') + ~(d') * ~x)) ^ ~n, ~x) => (~b / ~d) * Subst(Antiderivative((~x) ^ ~n / ((~b) ^ 2 + (~x) ^ 2), ~x), ~x, ~b * tan(~c + ~d * ~x)) <-- FreeQ([~b, ~c, ~d, ~n], ~x) && Not(IntegerQ(~n))
    @apply_utils Antiderivative((~a + ~(b') * tan(~(c') + ~(d') * ~x)) ^ 2, ~x) => ((~a) ^ 2 - (~b) ^ 2) * ~x + ((~b) ^ 2 * tan(~c + ~d * ~x)) / ~d + 2 * ~a * ~b * Antiderivative(tan(~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x)
    #= Antiderivative(((~a)+(~b')*tan((~c')+(~d')*(~x)))^(~n),~x) :=  Antiderivative(ExpandIntegrand((a+b*tan(c+d*x))^n,x),x) <-- FreeQ([a,b,c,d],x) &&  IGtQ(n,0) =#
    @apply_utils Antiderivative((~a + ~(b') * tan(~(c') + ~(d') * ~x)) ^ ~n, ~x) => (~b * (~a + ~b * tan(~c + ~d * ~x)) ^ (~n - 1)) / (~d * (~n - 1)) + 2 * ~a * Antiderivative((~a + ~b * tan(~c + ~d * ~x)) ^ (~n - 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (EqQ((~a) ^ 2 + (~b) ^ 2, 0) && GtQ(~n, 1))
    @apply_utils Antiderivative((~a + ~(b') * tan(~(c') + ~(d') * ~x)) ^ ~n, ~x) => (~a * (~a + ~b * tan(~c + ~d * ~x)) ^ ~n) / (2 * ~b * ~d * ~n) + (1 / (2 * ~a)) * Antiderivative((~a + ~b * tan(~c + ~d * ~x)) ^ (~n + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (EqQ((~a) ^ 2 + (~b) ^ 2, 0) && LtQ(~n, 0))
    @apply_utils Antiderivative(sqrt(~a + ~(b') * tan(~(c') + ~(d') * ~x)), ~x) => ((-2 * ~b) / ~d) * Subst(Antiderivative(1 / (2 * ~a - (~x) ^ 2), ~x), ~x, sqrt(~a + ~b * tan(~c + ~d * ~x))) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && EqQ((~a) ^ 2 + (~b) ^ 2, 0)
    @apply_utils Antiderivative((~a + ~(b') * tan(~(c') + ~(d') * ~x)) ^ ~n, ~x) => (-(~b) / ~d) * Subst(Antiderivative((~a + ~x) ^ (~n - 1) / (~a - ~x), ~x), ~x, ~b * tan(~c + ~d * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && EqQ((~a) ^ 2 + (~b) ^ 2, 0)
    @apply_utils Antiderivative((~a + ~(b') * tan(~(c') + ~(d') * ~x)) ^ ~n, ~x) => (~b * (~a + ~b * tan(~c + ~d * ~x)) ^ (~n - 1)) / (~d * (~n - 1)) + Antiderivative((((~a) ^ 2 - (~b) ^ 2) + 2 * ~a * ~b * tan(~c + ~d * ~x)) * (~a + ~b * tan(~c + ~d * ~x)) ^ (~n - 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && GtQ(~n, 1))
    @apply_utils Antiderivative((~a + ~(b') * tan(~(c') + ~(d') * ~x)) ^ ~n, ~x) => (~b * (~a + ~b * tan(~c + ~d * ~x)) ^ (~n + 1)) / (~d * (~n + 1) * ((~a) ^ 2 + (~b) ^ 2)) + (1 / ((~a) ^ 2 + (~b) ^ 2)) * Antiderivative((~a - ~b * tan(~c + ~d * ~x)) * (~a + ~b * tan(~c + ~d * ~x)) ^ (~n + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && LtQ(~n, -1))
    @apply_utils Antiderivative(1 / (~a + ~(b') * tan(~(c') + ~(d') * ~x)), ~x) => (~a * ~x) / ((~a) ^ 2 + (~b) ^ 2) + (~b / ((~a) ^ 2 + (~b) ^ 2)) * Antiderivative((~b - ~a * tan(~c + ~d * ~x)) / (~a + ~b * tan(~c + ~d * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && NeQ((~a) ^ 2 + (~b) ^ 2, 0)
    @apply_utils Antiderivative((~a + ~(b') * tan(~(c') + ~(d') * ~x)) ^ ~n, ~x) => (~b / ~d) * Subst(Antiderivative((~a + ~x) ^ ~n / ((~b) ^ 2 + (~x) ^ 2), ~x), ~x, ~b * tan(~c + ~d * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && NeQ((~a) ^ 2 + (~b) ^ 2, 0)
end

max_stars_repo_path: src/rules/4 Trig functions/4.3 Tangent/4.3.1.1 (a+b tan)^n.jl


text:
"""
    invgeotransform!(gt_in::Vector{Float64}, gt_out::Vector{Float64})

Invert Geotransform.

This function will invert a standard 3x2 set of GeoTransform coefficients. This
converts the equation from being pixel to geo to being geo to pixel.

### Parameters
* `gt_in`       Input geotransform (six doubles - unaltered).
* `gt_out`      Output geotransform (six doubles - updated).

### Returns
`gt_out`
"""
function invgeotransform!(
    gt_in::Vector{Float64},
    gt_out::Vector{Float64},
)::Vector{Float64}
    result = Bool(GDAL.gdalinvgeotransform(pointer(gt_in), pointer(gt_out)))
    result || error("Geotransform coefficients is uninvertable")
    return gt_out
end

invgeotransform(gt_in::Vector{Float64})::Vector{Float64} =
    invgeotransform!(gt_in, Array{Float64}(undef, 6))

"""
    applygeotransform(geotransform::Vector{Float64}, pixel::Float64,
        line::Float64)

Apply GeoTransform to x/y coordinate.

Applies the following computation, converting a (pixel,line) coordinate into a
georeferenced `(geo_x,geo_y)` location.
```
    geo_x = geotransform[1] + pixel*geotransform[2] + line*geotransform[3]
    geo_y = geotransform[4] + pixel*geotransform[5] + line*geotransform[6]
```
### Parameters
* `geotransform`  Six coefficient GeoTransform to apply.
* `pixel`           input pixel position.
* `line`            input line position.
"""
function applygeotransform(
    geotransform::Vector{Float64},
    pixel::Float64,
    line::Float64,
)::Vector{Float64}
    geo_xy = Vector{Float64}(undef, 2)
    geo_x = pointer(geo_xy)
    geo_y = geo_x + sizeof(Float64)
    GDAL.gdalapplygeotransform(pointer(geotransform), pixel, line, geo_x, geo_y)
    return geo_xy
end

"""
    composegeotransform!(gt1::Vector{Float64}, gt2::Vector{Float64},
        gtout::Vector{Float64})

Compose two geotransforms.

The resulting geotransform is the equivalent to `padfGT1` and then `padfGT2`
being applied to a point.

### Parameters
* `gt1`     the first geotransform, six values.
* `gt2`     the second geotransform, six values.
* `gtout`   the output geotransform, six values.
"""
function composegeotransform!(
    gt1::Vector{Float64},
    gt2::Vector{Float64},
    gtout::Vector{Float64},
)::Vector{Float64}
    GDAL.gdalcomposegeotransforms(pointer(gt1), pointer(gt2), pointer(gtout))
    return gtout
end

function composegeotransform(
    gt1::Vector{Float64},
    gt2::Vector{Float64},
)::Vector{Float64}
    return composegeotransform!(gt1, gt2, Vector{Float64}(undef, 6))
end


text:
#=
jinc.jl
=#

export jinc

using SpecialFunctions: besselj1 #grab bessel function
using Plots: plot
using Test: @inferred

"""
    jinc(x)
return `jinc(x) = J1(pi*x)/(2x)`, where `J1` is a Bessel function of the first kind.
units of `x` are typically cycles/m
return type is `promote_type(typeof(x), Float32)`
"""
function jinc(x::Real)
    T = promote_type(typeof(x), Float32)
    if (x == 0) return convert(T, pi/4) end
    y = abs(x)
    convert(T, besselj1(pi*y) / (2*y))
end

"""
    jinc(:test) jinc(:plot)
self test, plot
"""
function jinc(test::Symbol)
    if test === :plot
        r = LinRange(-10,10,201)
        return plot(r, jinc.(r))
    end
    !(test === :test) && throw("non-test, symbolic input to jinc")
    jinx = x -> jinc.(x)
    r = LinRange(-10,10,201)
    @inferred jinx(r)
    for T in (Int, Float16, Float32, Float64)
        @inferred jinc(T(0))
        @inferred jinc(T(1))
    end
    jinc(:plot)
    true
end


text:
module Stacks

export NNTopo, @nntopo_str, @nntopo
export Stack, show_stackfunc, stack

include("./code.jl")
include("./topology.jl")
include("./stack.jl")

end


text:
x0 = 1.0
P0 = 1.0

Φ = 0.5
b = 0.0
Q = 2.0

yshadow = 0.0
H = 1.0
R = 1.0

E = LinearEvolution(Φ, Gaussian(b, Q))
Obs = LinearObservationModel(H, R)
M = LinearStateSpaceModel(E, Obs)

O = LinearObservation(LinearEvolution(Φ, Gaussian(b, Q)), H, R)


text:
# JuMP Variable interface
# ------------------------------------------------------------------------------

# main interface

JuMP.is_fixed(p::ParameterRef) = true
JuMP.FixRef(p::ParameterRef) =
    error("Parameters do not have have explicit constraints, hence no constraint reference.")
JuMP.unfix(p::ParameterRef) = error("Parameters cannot be unfixed.")

function JuMP.fix(p::ParameterRef, val::Real)
    data = _getparamdata(p)::ParameterData
    data.sync = false
    data.future_values[index(data, p)] = val
    return nothing
end

"""
    fix(p::ParameterRef, val::Real)::Nothing

Sets the parameter `p` to the new value `val`.
"""
function fix(p::ParameterRef, val::Real)
    params = _getparamdata(p)::ParameterData
    params.sync = false
    params.future_values[index(p)] = val
    return nothing
end

function JuMP.value(p::ParameterRef)
    data = _getparamdata(p)::ParameterData
    data.future_values[index(data, p)]
end

# interface continues

JuMP.owner_model(p::ParameterRef) = p.model

struct ParameterNotOwned <: Exception
    parameter::ParameterRef
end

function JuMP.check_belongs_to_model(p::ParameterRef, model::AbstractModel)
    if owner_model(p) !== model
        throw(ParameterNotOwned(p))
    end
end

Base.iszero(::ParameterRef) = false
Base.copy(p::ParameterRef) = ParameterRef(p.ind, p.model)
# Base.broadcastable(v::VariableRef) = Ref(v) # NEEDED???

"""
    delete(model::Model, param::ParameterRef)

Delete the parameter `param` from the model `model`.

Note.
After the first deletion you might experience performance reduction.
Therefore, only use thid command if there is no other way around.
"""
function JuMP.delete(model::Model, param::ParameterRef)
    error("Parameters can be deleted currently.")
    if model !== owner_model(param)
        error("The variable reference you are trying to delete does not " *
              "belong to the model.")
    end
    # create dictionary map
    # turn flag has_deleted
    # delete names ?
end

"""
    is_valid(model::Model, parameter::ParameterRef)

Return `true` if `parameter` refers to a valid parameter in `model`.
"""
function JuMP.is_valid(model::Model, parameter::ParameterRef)
    return model === owner_model(parameter)
end

# The default hash is slow. It's important for the performance of AffExpr to
# define our own.
# https://github.com/JuliaOpt/MathOptInterface.jl/issues/234#issuecomment-366868878
function Base.hash(p::ParameterRef, h::UInt)
    return hash(objectid(owner_model(p)), hash(p.ind, h))
end
function Base.isequal(p1::ParameterRef, p2::ParameterRef)
    return owner_model(p1) === owner_model(p2) && p1.ind == p2.ind
end

function JuMP.name(p::ParameterRef)
    dict = _getparamdata(p).names
    if haskey(dict, p)
        return dict[p]
    else
        return ""
    end
end

function JuMP.set_name(p::ParameterRef, s::String)
    dict = _getparamdata(p).names
    dict[p] = s
end

function parameter_by_name(model::Model, name::String)
    # can be improved with a lazy rev_names map
    dict = _getparamdata(model).names
    for (par, n) in dict
        if n == name
            return par
        end
    end
    return nothing
end

JuMP.has_lower_bound(p::ParameterRef) = false
JuMP.LowerBoundRef(p::ParameterRef) =
    error("Parameters do not have bounds.")
JuMP.set_lower_bound(p::ParameterRef, lower::Number) =
    error("Parameters do not have bounds.")
JuMP.delete_lower_bound(p::ParameterRef) =
    error("Parameters do not have bounds.")
JuMP.lower_bound(p::ParameterRef) =
    error("Parameters do not have bounds.")

JuMP.has_upper_bound(p::ParameterRef) = false
JuMP.UpperBoundRef(p::ParameterRef) =
    error("Parameters do not have bounds.")
JuMP.set_upper_bound(p::ParameterRef, lower::Number) =
    error("Parameters do not have bounds.")
JuMP.delete_upper_bound(p::ParameterRef) =
    error("Parameters do not have bounds.")
JuMP.upper_bound(p::ParameterRef) =
    error("Parameters do not have bounds.")

JuMP.is_integer(p::ParameterRef) = false
JuMP.set_integer(p::ParameterRef) =
    error("Parameters do not have integrality constraints.")
JuMP.unset_integer(p::ParameterRef) =
    error("Parameters do not have integrality constraints.")
JuMP.IntegerRef(p::ParameterRef) =
    error("Parameters do not have integrality constraints.")

JuMP.is_binary(p::ParameterRef) = false
JuMP.set_binary(p::ParameterRef) =
    error("Parameters do not have binary constraints.")
JuMP.unset_binary(p::ParameterRef) =
    error("Parameters do not have binary constraints.")
JuMP.BinaryRef(p::ParameterRef) =
    error("Parameters do not have binary constraints.")

JuMP.start_value(p::ParameterRef) =
    error("Parameters do not have start values.")
JuMP.set_start_value(p::ParameterRef, value::Number) =
    error("Parameters do not have start values.")


text:
using ThreeBodyProblem
using DifferentialEquations
using Plots

gr() # This command sets GR as our plotting backend

μ₁ = 398600 # {km³/s²} gravitational parameter of Earth
μ₂ = 4902   # {km³/s²} gravitational parameter of the Moon
d = 384400  # {km} average distance between Earth and the Moon

p = [μ₁, μ₂, d]

L1 = computeL1(p) # calculate the Lagrange point L1 for the given system parameters p

Lpts = computeLpts(p)
L4 = Lpts[4]

R₁,R₂ = computeR1R2(p) # {km} distances of Primary and Secondary bodies from origin
ωₛ = sqrt((μ₁ + μ₂)/d^3) # {rad/s} rotation rate of system

N = 100
X = range(-1.5*d,1.5*d,length=N)
Y = range(-d,d,length=N)

gr()
f(x,y) = begin
    -(x^2 + y^2)*ωₛ^2/2 - μ₁/sqrt((x+R₁)^2 + y^2) - μ₂/sqrt((x-R₂)^2 + y^2)
    # rv = [x; y; 0; 0; 0; 0]
    # findUeff(rv,p)
end
contour(X,Y,f,levels=200,fill=true)

Rₑ = 6378.0 # {km} radius of the Earth
Rₘ = 1738.0 # {km} radius of the Moon

h = 200.0           # {km} altitude of parking orbit
vᵢ = 10.92367104    # {km/s} injection velocity in rotating frame
ϕ = 47.70061087     # {°} injection angle, measured from +y

r₀ = [-R₁ - (Rₑ + 200)*cosd(ϕ); -(Rₑ + 200)*sind(ϕ); 0]
v₀ = vᵢ*[sind(ϕ); -cosd(ϕ); 0];
rv₀ = [r₀;v₀] # {km; km/s} our initial state

tspan = (0.,86400*6.) # {sec} 6 day span

prob = ODEProblem(CR3BPdynamics!,rv₀,tspan,p) # CR3BPdynamics! is our in-place dynamics function
sol = solve(prob,reltol=1e-6)

plot(sol,vars=(1,2),title="Lunar return trajectory in rotating frame")
plot!(circle(Rₑ,[-R₁;0]),color="blue",label="Earth")
plot!(circle(Rₘ,[R₂;0]),color="gray",label="Moon",aspect_ratio=:equal)

rv_inert = [rot2inert(sol.u[i],sol.t[i],p) for i = 1:length(sol)]
x = [rv_inert[i][1] for i = 1:length(rv_inert)]
y = [rv_inert[i][2] for i = 1:length(rv_inert)]
plot(x,y,legend=:topleft)


text:
module BasicPatterns
using MatchCore.AbstractPatterns
using MatchCore.AbstractPatterns.RedyFlavoured

export P_bind, P_tuple, P_type_of, P_vector, P_capture, P_vector3, P_slow_view, P_fast_view
export P_svec, P_svec3
export SimpleCachablePre, see_captured_vars, see_captured_vars!

const EQ = 0b001
const GT = 0b100
const LT = 0b010

@nospecialize
OptionalLn = Union{LineNumberNode,Nothing}

function see_captured_vars(inner::Any, in_scope::ChainDict{Symbol, Symbol})
    bind = Expr(:block)
    for_chaindict(in_scope) do k, v
        push!(bind.args, :($k = $v))
    end
    isempty(bind.args) ? inner : Expr(:let, bind, inner)
end

function see_captured_vars!(inner::Any, in_scope::ChainDict{Symbol, Symbol})
    bind = Expr(:block)
    for_chaindict(in_scope) do k, v
        if k !== v
            assign = :($k = $v)
            push!(bind.args, assign)
        end
    end
    isempty(bind.args) ? inner : Expr(:let, bind, inner)
end

struct SimpleCachablePre <: APP
    f :: Function
end
(f::SimpleCachablePre)(target) = f.f(target)

function sequence_index(viewed, i::Integer, ::Any, ::Any)
    :($viewed[$i])
end

function self_index(viewed, i::Integer, ::Any, ::Any)
    @assert i === 1
    viewed
end

function length_eq_check(seq, n::Int)
    if n === 0
        :(isempty($seq))
    else
        :(length($seq) === $n)
    end
end

function mk_type_object(i::Int, ::Type{T}) where {T}
    if isabstracttype(T)
        TypeVar(Symbol(:var, i), T)
    else
        T
    end
end

"""match by type
"""
function P_type_of(t, prepr::AbstractString="isa $t")
    recog_type() = t
    comp = PComp(prepr, recog_type)
    decons(comp, [])
end

"""bind a symbol
"""
function P_bind(n::Symbol, expr::Any; see_capture=false)
    function bind_effect!(target, scope::ChainDict{Symbol,Symbol}, ln::OptionalLn)
        expr′ = see_capture ? see_captured_vars(expr, scope) : expr
        n′ = scope[n] = gensym(n)
        :($n′ = $expr′)
    end
    effect(bind_effect!)
end

"""bind a symbol
"""
function P_capture(n::Symbol)
    function capture_effect!(target, scope::ChainDict{Symbol,Symbol}, ln::OptionalLn)
        if target isa Symbol
            scope[n] = target
            return nothing
        end
        n′ = scope[n] = gensym(n)
        :($(n′) = $target)
    end
    effect(capture_effect!)
end

"""deconstruct a tuple
"""
function P_tuple(fields::AbstractArray, prepr::AbstractString="Tuple")
    function type_of_tuple(xs...)

        ts = [mk_type_object(i, xs[i]) for i in eachindex(xs)]
        foldl(ts, init=Tuple{ts...}) do last, t
            t isa TypeVar ?
                UnionAll(t, last) :
                last
        end
    end
    comp = PComp(prepr, type_of_tuple)

    decons(comp, sequence_index, fields)
end

function type_of_vector(types...)
    AbstractArray
end
"""deconstruct a vector
"""
function P_vector(fields::AbstractArray, prepr::AbstractString="1DVector")

    n_fields = length(fields)

    function pred(target)
        length_eq_check(target, n_fields)
    end
    comp = PComp(prepr, type_of_vector; guard1=NoncachablePre(pred))
    decons(comp, sequence_index, fields)
end

"""deconstruct a vector
"""
function P_svec(fields::AbstractArray, prepr::AbstractString="svec")
    function type_of_svec(_...)
        Core.SimpleVector
    end
    n_fields = length(fields)
    function pred(target)
        :(ndims($target) === 1 && length($target) === $n_fields)
    end
    comp = PComp(prepr, type_of_svec; guard1=NoncachablePre(pred))
    decons(comp, sequence_index, fields)
end

"""deconstruct a vector in this way: [a, b, c, pack..., d, e]
"""
function P_vector3(init::AbstractArray, pack::Function, tail::AbstractArray, prepr::AbstractString = "1DVector Pack")
    n1 = length(init)
    n2 = length(tail)
    min_len = length(init) + length(tail)
    function extract(arr, i::Int, ::Any, ::Any)
        ex = if i <= n1
            :($arr[$i])
        elseif i === n1 + 1
            n2 === 0 ?
            :(  $SubArray($arr, ($(n1+1):length($arr),     )) ) :
            :(  $SubArray($arr, ($(n1+1):length($arr)-$n2, )) )
        else
            incr = i - n1 - 1
            j = n2 - incr
            ex = j == 0 ? :($arr[end]) : :($arr[end-$j])
            :($ex)
        end
    end
    function pred(target)
        :(ndims($target) === 1 && length($target) >= $min_len)
    end
    comp = PComp(prepr, type_of_vector; guard1=NoncachablePre(pred))
    decons(comp, extract, [init; pack; tail])
end

"""deconstruct a vector in this way: [a, b, c, pack..., d, e]
"""
function P_svec3(init::AbstractArray, pack::Function, tail::AbstractArray, prepr::AbstractString = "svec Pack")
    n1 = length(init)
    n2 = length(tail)
    min_len = length(init) + length(tail)
    function type_of_svec(types...)
        Core.SimpleVector
    end
    function extract(arr, i::Int, ::Any, ::Any)
        if i <= n1
            :($arr[$i])
        elseif i === n1 + 1
            n2 === 0 ?
                :($arr[$(n1+1):end]) :
                :($arr[$(n1+1):end-$n2])
        else
            incr = i - n1 - 1
            j = n2 - incr
            j == 0 ? :($arr[end]) : :($arr[end-$j])
        end
    end
    function pred(target)
        :(length($target) >= $min_len)
    end
    comp = PComp(prepr, type_of_svec; guard1=NoncachablePre(pred))
    decons(comp, extract, [init; pack; tail])
end

"""untyped view pattern
"""
function P_slow_view(trans, p::Function, prepr::AbstractString="ViewBy($trans)")
    function type_of_slow_view(args...)
        Any
    end

    comp = PComp(
        prepr, type_of_slow_view;
        view=SimpleCachablePre(trans)
    )
    decons(comp, self_index, [p])
end

"""typed view pattern
"""
function P_fast_view(tcons, trans, p::Function, prepr::AbstractString="ViewBy($trans, typecons=$tcons)")

    comp = PComp(
        prepr, tcons;
        view=SimpleCachablePre(trans)
    )
    decons(comp, self_index, [p])
end

@specialize
end