text:
# Author: Ritchie Lee, ritchie.lee@sv.cmu.edu
# Date: 11/21/2014

#ACASX implementation: based on CCAS interface to libCAS

module ACASX_CCAS_Impl

export
    addObserver,

    initialize,
    update,

    ACASX_CCAS,
    ACASXInput,
    ACASXOutput

import Compat.ASCIIString

using AbstractCollisionAvoidanceSystemImpl
using AbstractCollisionAvoidanceSystemInterfaces
using CommonInterfaces

using AbstractCASCoordImpl

using ACASXCommonImpl
using CASCoordination
using CASInterface
using CCAS
using RLESUtils, Observers

import CommonInterfaces.initialize
import CommonInterfaces.update

type ACASX_CCAS <: AbstractCollisionAvoidanceSystem
  my_id::Int64 #aircraft number
  max_intruders::Int64
  quant::Int64
  constants::Constants
  casShared::CASShared
  version::ASCIIString
  equipage::EQUIPAGE
  input::Input
  output::Output
  coord::AbstractCASCoord

  function ACASX_CCAS(aircraft_id::Int64, libcas::AbstractString, config_file::AbstractString,
                 num_aircraft::Int, coord::AbstractCASCoord, equipage::EQUIPAGE=EQUIPAGE_TCAS)
    cas = new()
    cas.my_id = aircraft_id
    cas.max_intruders = num_aircraft - 1
    cas.quant = equipage == EQUIPAGE_TCAS ? 25 : 100
    cas.constants = Constants(cas.quant, config_file, cas.max_intruders)
    cas.casShared = CASShared(libcas, cas.constants)
    cas.version = version(cas.casShared)
    cas.equipage = equipage

    @assert cas.max_intruders == max_intruders(cas.casShared) #Will fail if library did not open properly

    cas.coord = coord
    cas.input = Input(cas.max_intruders)
    cas.output = Output(cas.max_intruders)
    setRecord(cas.coord, cas.my_id,
              ACASXCoordRecord(cas.my_id, equipage, cas.quant, cas.max_intruders))
    reset(cas.casShared)
    return cas
  end
end

function update(cas::ACASX_CCAS, input::ACASXInput)
  ACASXCommonImpl.update_from_coord!(input, cas.coord, cas.my_id)
  if cas.equipage == EQUIPAGE_TCAS
    update!(cas.casShared, input, cas.output) #cas.output is modified in place
    ACASXCommonImpl.update_to_coord!(cas.coord, cas.my_id, cas.output)
  end

  cas.input = input #keep a record
  return cas.output
end

function initialize(cas::ACASX_CCAS)
  ACASXCommonImpl.initialize(cas)
  reset(cas.casShared)
end

end #module


text:
# Abstract type hierarchy
abstract type Asset end

abstract type Property <: Asset end
abstract type Investment <: Asset end
abstract type Cash <: Asset end

abstract type House <: Property end
abstract type Apartment <: Property end

abstract type FixedIncome <: Investment end
abstract type Equity <: Investment end

# simple functions on abstract types
describe(a::Asset) = "Something valuable"
describe(e::Investment) = "Financial investment"
describe(e::Property) = "Physical property"

# function that just raise an error
"""
    location(p::Property)

Returns the location of the property as a tuple of (latitude, longitude).
"""
location(p::Property) = error("Location is not defined in the concrete type")

# an empty function
"""
    location(p::Property)

Returns the location of the property as a tuple of (latitude, longitude).
"""
function location(p::Property) end

# interaction function that works at the abstract type level
function walking_disance(p1::Property, p2::Property)
    loc1 = location(p1)
    loc2 = location(p2)
    return abs(loc1.x - loc2.x) + abs(loc1.y - loc2.y)
end

ext: jl


text:
"""
# Notes
- Currently, only iip (isinplace) method is supported.
- Default solver: Tsit5() in OrdinaryDiffEq.jl
"""
function sim(state0, dyn, p=nothing;
        solver=Tsit5(),
        kwargs...
    )
    FSimBase.sim(state0, dyn, p;
                 solver=solver,
                 kwargs...,
                )
end


text:
module TestConstants

using Test
using Canopy.Constants

"""Read CODATA 2014 physical constants table."""
function read_codata(filename)
    codata_constants = Dict{String, Float64}()
    codata_raw = open(readlines, filename)
    for line in codata_raw[eachindex(codata_raw)[10:end]]
        if length(line) < 85
            break
        end
        key = strip(line[1:59])
        val = parse(Float64, replace(line[61:85], r"(\s+|\.\.\.)" => ""))
        push!(codata_constants, key => val)
    end
    return codata_constants
end

codata_constants = read_codata((@__DIR__) * "/assets/codata/codata2014.txt")

validation_pairs = [
    (c, codata_constants["speed of light in vacuum"]),
    (h, codata_constants["Planck constant"]),
    (N_A, codata_constants["Avogadro constant"]),
    (k_B, codata_constants["Boltzmann constant"]),
    (sigma_SB, codata_constants["Stefan-Boltzmann constant"]),
    (c_1, codata_constants["first radiation constant"]),
    (c_1L, codata_constants["first radiation constant for spectral radiance"]),
    (c_2, codata_constants["second radiation constant"]),
    (R, codata_constants["molar gas constant"]),
    (g, codata_constants["standard acceleration of gravity"]),
    (atm, codata_constants["standard atmosphere"]),
    (T_0, 273.15),
    (kappa, 0.4),
    (M_w, 18.015e-3),
    (M_d, 28.965e-3)]

@time @testset "Canopy.Constants" begin
    for (val1, val2) in validation_pairs
        @test isapprox(val1, val2, rtol=1e-3)
    end
end

end  # module


text:
function thingorfunc(code, cursor, mod = Main; name = getqualifiedname(code, cursor))
  name == "" && (name = lastcall(scopes(code, cursor)))
  name == nothing ? name : getthing(mod, name, nothing)
end

function doc(code, cursor, mod = Main)
  docs = AbstractString[]
  name = getqualifiedname(code, cursor)

  thing = thingorfunc(code, cursor, mod; name = name)
  thing == nothing || push!(docs, sprint(Base.help, thing))

  texcmds = texcommands(name, code, cursor)
  texcmds == nothing || push!(docs, texcmds)

  help = join(docs, "\n\n")
  help in ("No help information found.\n", "") ? nothing : help
end

function texcommands(name, code, cursor)
  chars = collect(name)
  if isempty(chars) # fallback if getqualifiedname failed
    line = collect(lines(code)[cursor.line])
    c = cursor.column
    chars = line[max(c - 1, 1):min(c, length(line))]
  end

  syms = AbstractString[]
  for char in chars
    cmd = get(reverse_latex_commands, char, "")
    isempty(cmd) || push!(syms, "  $(char)\u00a0 $(cmd)")
  end
  isempty(syms) && return

  "LaTeX command$(length(syms) > 1 ? "s" : ""):\n\n$(join(unique(syms), "\n"))"
end

function methodsorwith(code, cursor, mod = Main)
  thing = thingorfunc(code, cursor, mod)
  thing == nothing && return
  return thing == Module ? methodswith(Module) :
         (isa(thing, Function) && isgeneric(thing)) || isleaftype(thing) ? methods(thing) :
         eval(Main, :(methodswith($(typeof(thing)), true))) # I have no idea why I thought this was necessary
end


text:
# For custom types, test get/setproperty for each "kind"
# of wrapped field:
# 1) isbits types like Cint/SArray (e.g. jlModel.nq or jlData.solver) -> pass through to inner mj type (e.g. jlData.d::mjData)
# 2) Ptr to "statically" sized Array (e.g. d.qpos) -> Base.Array
# 3) Ptr to dynamically sized UnsafeArray (e.g. d.contact) -> Dynamically sized UnsafeArrays.UnsafeArray

@testset "Sugar" begin

    @testset "@set!!" begin
        jm = jlModel(TESTMODELXML)
        jd = jlData(jm)
        mj_resetData(jm, jd)
        @test jm.opt.timestep == 0.001
        qpos0 = copy(jd.qpos)
        mj_step(jm, jd)
        @set!! jm.opt.timestep = -0.001
        @test jm.opt.timestep == -0.001
        mj_step(jm, jd)
        @test isapprox(qpos0, jd.qpos; atol = 1e-3)
        @test jd.time == 0.0
    end

    @testset "Basic" begin

        @testset "jlModel" begin
            _m, _d, _pm, _pd = makemd()
            _jm = @inferred jlModel(_pm)
            f = x -> x.nq
            @test @inferred(f(_jm)) isa Cint

            f = x -> x.body_mass
            @test @inferred(f(_jm)) isa Array{mjtNum,1}

            jm1 = @inferred(jlModel(TESTMODELXML))
            jm2 = @inferred(jlModel(TESTMODELMJB))
            @testset "$name" for name in propertynames(jm1)
                f1 = getproperty(jm1, name)
                if typeof(f1) <: Ptr
                    true
                else
                    f2 = getproperty(jm2, name)
                    ret = isapprox(f1, f2) # TODO
                    #ret=if Sys.iswindows() && typeof(f1) <: AbstractArray
                    #    # For some reason on Windows, subsequent allocations yield oh-so-slightly different values
                    #    return isapprox(f1, f2)
                    #else
                    #    return f1 == f2
                    #end
                    #if !ret
                    #    @warn name
                    #    @info f1 f2
                    #end
                end
            end

            let m1 = jlModel(TESTMODELXML)
                m2 = copy(m1)
                m3 = deepcopy(m1)
                m4 = jlModel(m1)
                @test all(propertynames(jlModel)) do n
                    x1 = getproperty(m1, n)
                    x2 = getproperty(m2, n)
                    x3 = getproperty(m3, n)
                    x4 = getproperty(m4, n)
                    if x1 isa Number || x2 isa AbstractArray
                        return x1 == x2 == x3 == x4
                    else
                        return true
                    end
                end
            end

        end

        @testset "jlData" begin
            _m, _d, _pm, _pd = makemd()
            _jd = @inferred jlData(_pm, _pd)
            f = x -> x.ncon
            @test @inferred(f(_jd)) isa Cint

            f = x -> x.qpos
            @test @inferred(f(_jd)) isa Array{mjtNum,1}

            f = x -> x.contact
            @test @inferred(f(_jd)) isa UnsafeArray{mjContact,1}
            @test length(_jd.contact) == _d.ncon
            @test length(getfield(_jd, :contact)) == _m.nconmax

            f = x -> x.solver
            @test @inferred(f(_jd)) isa Array{mjSolverStat,1}
            @test length(_jd.solver) == length(_d.solver) == mjNSOLVER
        end

        @testset "Simulation" begin
            m, d, _, _ = jlmakemd()
            jm = jlModel(TESTMODELXML)
            jd = jlData(jm)
            mj_resetData(m, d)
            mj_resetData(jm, jd)

            @test m.opt.timestep == jm.opt.timestep == 0.001

            qpos0 = copy(jd.qpos)

            mj_step(m, d)
            mj_step(jm, jd)
            qpos1 = copy(jd.qpos)

            @test jd.qpos != qpos0
            @test unsafe_wrap(Array, d.qpos, size(qpos0)) == jd.qpos
            @test d.time == jd.time == 0.001

            # backwards step
            @set!! m.opt.timestep = -0.001
            @set!! jm.m.opt.timestep = -0.001

            mj_step(m, d)
            mj_step(jm, jd)

            @test isapprox(qpos0, unsafe_wrap(Array, d.qpos, size(qpos0)); atol = 1e-3)
            @test isapprox(qpos0, jd.qpos; atol = 1e-3)

            @test d.time == 0.0
            @test jd.time == 0.0

            mj_resetData(m, d)
            mj_resetData(jm, jd)

            @test qpos0 == unsafe_wrap(Array, d.qpos, size(qpos0)) == jd.qpos
            @test d.time == jd.time == 0.0

            @set!! m.opt.timestep = 0.001
            @set!! jm.m.opt.timestep = 0.001

            mj_step1(m, d)
            mj_step2(m, d)

            mj_step1(jm, jd)
            mj_step2(jm, jd)

            @test jd.qpos == unsafe_wrap(Array, d.qpos, size(qpos0)) == qpos1
        end
    end

    # Test one of each of our custom axes
    @testset "namify" begin
        m = jlModel(TESTMODELXML)
        d = jlData(m)
        mn, dn = namify(m, d)

        ans = [
            0 1 0
            1 0 0
            0 0 -1
        ]
        @test isapprox(dn.xmat[:, :, :box1_body], ans)
        @test isapprox(dn.xquat[:w, :box1_body], 0.0)
        @test isapprox(dn.xquat[:x, :box1_body], sqrt(2)/2)
        @test isapprox(dn.xquat[:y, :box1_body], sqrt(2)/2)
        @test isapprox(dn.xquat[:z, :box1_body], 0.0)

        @test mn.actuator_forcerange[:min, :box1_motor] == -1.0
        @test mn.actuator_forcerange[:max, :box1_motor] == 1.0

        @test mn.mat_rgba[:] == Float32[0.1, 0.2, 0.3, 0.4]

        @test all(propertynames(mn)) do name
            fn = getproperty(mn, name)
            f = getproperty(m, name)
            length(fn) == length(f)
        end

        @test all(propertynames(dn)) do name
            fn = getproperty(dn, name)
            f = getproperty(d, name)
            length(fn) == length(f)
        end
    end

    @testset "jlmujoco" begin
        jm = jlModel(TESTMODELXML)
        @test jl_enabled(jm, MJCore.mjENBL_ENERGY)
        @test jl_disabled(jm, MJCore.mjDSBL_CLAMPCTRL)
    end

end


text:
"""
    SRM0

A SRM0 neuron described by

``v(t) = \\eta(t - t^f) \\Theta(t - t^f) + I``

  where ``\\Theta`` is the Heaviside function and
  ``t^f`` is the last output spike time.

*Note:* The SRM0 model normally includes an EPSP term which
  is modeled by [`Synapse.EPSP`](@ref)

For more details see:
  [Spiking Neuron Models: Single Neurons, Populations, Plasticity]
  (https://icwww.epfl.ch/~gerstner/SPNM/node27.html#SECTION02323400000000000000)

Fields:
- `voltage::VT`: membrane potential
- `current::VT`: injected (unprocessed) current
- `lastspike::VT`: last time this neuron spiked
- `η::F`: refractory response function
"""
mutable struct SRM0{VT<:Real, F<:Function} <: AbstractCell
    voltage::VT
    current::VT
    lastspike::VT
    η::F
end

Base.show(io::IO, ::MIME"text/plain", neuron::SRM0) =
    print(io, """SRM0:
                     voltage:   $(neuron.voltage)
                     current:   $(neuron.current)
                     lastspike: $(neuron.lastspike)
                     η:         $(neuron.η)""")
Base.show(io::IO, neuron::SRM0) =
    print(io, "SRM0(voltage: $(neuron.voltage), current = $(neuron.current))")

"""
    SRM0{T}(η) where T<:Real
    SRM0(η)

Create a SRM0 neuron with refractory response function `η`.
"""
SRM0{T}(η::F) where {T<:Real, F<:Function} = SRM0{T, F}(0, 0, -Inf, η)
SRM0(η::Function) = SRM0{Float32}(η)

"""
    SRM0(η₀, τᵣ, v_th)

Create a SRM0 neuron with refractory response function:

``-\\eta_0 \\exp\\left(-\\frac{\\Delta}{\\tau_r}\\right)``
"""
function SRM0{T}(η₀::Real, τᵣ::Real) where T<:Real
    η = Δ -> -η₀ * exp(-Δ / τᵣ)
    SRM0{T}(η)
end
SRM0(η₀::Real, τᵣ::Real) = SRM0{Float32}(η₀, τᵣ)

isactive(neuron::SRM0, t::Integer; dt::Real = 1.0) = (neuron.current > 0)
getvoltage(neuron::SRM0) = neuron.voltage

"""
    excite!(neuron::SRM0, current)
    excite!(neurons::AbstractArray{<:SRM0}, current)

Excite an SRM0 `neuron` with external `current`.
"""
excite!(neuron::SRM0, current) = (neuron.current += current)
excite!(neurons::T, current) where T<:AbstractArray{<:SRM0} = (neurons.current .+= current)

"""
    spike!(neuron::SRM0, t::Integer; dt::Real = 1.0)
    spike!(neurons::AbstractArray{<:SRM0}, spikes; dt::Real = 1.0)

Record a output spike from the threshold function with the `neuron` body.
Sets `neuron.lastspike`.
"""
spike!(neuron::SRM0, t::Integer; dt::Real = 1.0) = (t > 0) && (neuron.lastspike = dt * t)
function spike!(neurons::T, spikes; dt::Real = 1.0) where T<:AbstractArray{<:SRM0}
    map!((x, s) -> (s > 0) ? dt * s : x, neurons.lastspike, neurons.lastspike, spikes)
end

"""
    srm0(t::Real, I, V; lastspike, eta)
    srm0!(t::Real, I::AbstractArray{<:Real}, V::AbstractArray{<:Real}; lastspike::AbstractArray{<:Real}, eta)
    srm0!(t::Real, I::CuVector{<:Real}, V::CuVector{<:Real}; lastspike::CuVector{<:Real}, eta)

Evaluate a SRM0 neuron.
Use `CuVector` instead of `Vector` to evaluate on GPU.

# Fields
- `t`: current time in seconds
- `I`: external current
- `V`: current membrane potential
- `lastspike`: time of last output spike in seconds
- `eta`: post-synaptic response function
"""
srm0(t::Real, I, V; lastspike, eta) = eta(t - lastspike) + I
function srm0!(t::Real, I::AbstractArray{<:Real}, V::AbstractArray{<:Real}; lastspike::AbstractArray{<:Real}, eta)
    V .= map.(eta, (t .- lastspike))
    V .+= I
end
srm0!(t::Real, I::CuVector{<:Real}, V::CuVector{<:Real}; lastspike::CuVector{<:Real}, eta) =
    V .= map.(eta, t .- lastspike) .+ I

"""
    evaluate!(neuron::SRM0, t::Integer; dt::Real = 1.0)
    (neuron::SRM0)(t::Integer; dt::Real = 1.0)
    evaluate!(neurons::AbstractArray{<:SRM0}, t::Integer; dt::Real = 1.0)

Evaluate the neuron model at time `t`. Return resulting membrane potential.
"""
function evaluate!(neuron::SRM0, t::Integer; dt::Real = 1.0)
    neuron.voltage = srm0(t * dt, neuron.current, neuron.voltage; lastspike = neuron.lastspike, eta = neuron.η)
    neuron.current = 0

    return neuron.voltage
end
(neuron::SRM0)(t::Integer; dt::Real = 1.0) = evaluate!(neuron, t; dt = dt)
function evaluate!(neurons::T, t::Integer; dt::Real = 1.0) where T<:AbstractArray{<:SRM0}
    srm0!(t * dt, neurons.current, neurons.voltage; lastspike = neurons.lastspike, eta = neurons.η)
    neurons.current .= 0

    return neurons.voltage
end

"""
    reset!(neuron::SRM0)
    reset!(neurons::AbstractArray{<:SRM0})

Reset `neuron`.
"""
function reset!(neuron::SRM0)
    neuron.voltage = 0
    neuron.lastspike = -Inf
    neuron.current = 0
end
function reset!(neurons::T) where T<:AbstractArray{<:SRM0}
    neurons.voltage .= 0
    neurons.lastspike .= -Inf
    neurons.current .= 0
end


text:
# This file is a part of ChangesOfVariables.jl, licensed under the MIT License (MIT).

"""
    ChangesOfVariables

Lightweight package that defines functionality to calculate volume element
changes for functions that perform a change of variables (like coordinate
transformations).
"""
module ChangesOfVariables

using ChainRulesCore
using LinearAlgebra
using Test

include("with_ladj.jl")
include("test.jl")

end # module


text:
@testset "multilevel_PARSDMM" begin

Random.seed!(123)

#test serial multilevel PARSDMM
    TF = Float32
    options = PARSDMM_options()
    default_PARSDMM_options(options,TF)
    options.parallel     = false
    options.evol_rel_tol = 100*eps(TF)
    options.maxit        = 5000

    n_levels          = 2
    coarsening_factor = 1.87

    x         = randn(TF,100,201)
    comp_grid = compgrid((1.0, 1.0),(100, 201))
    x         = vec(x)
    m         = deepcopy(x)

    constraint = Vector{SetIntersectionProjection.set_definitions}()

    #total variation
    (TD_OP, AtA_diag, dense, TD_n)=get_TD_operator(comp_grid,"TV",options.FL)
    m_min     = 0.0
    m_max     = 0.5f0*norm(TD_OP*x,1)
    set_type  = "l1"
    TD_OP     = "TV"
    app_mode  = ("matrix","")
    custom_TD_OP = ([],false)
    push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

    #(P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid,options.FL)
    #(TD_OP,AtA,l,y)        = PARSDMM_precompute_distribute(TD_OP,set_Prop,comp_grid,options)
    (TD_OP_levels,AtA_levels,P_sub_levels,set_Prop_levels,comp_grid_levels)=setup_multi_level_PARSDMM(m,n_levels,coarsening_factor,comp_grid,constraint,options)

    (x1,log_PARSDMM) = PARSDMM(m,AtA,TD_OP,set_Prop,P_sub,comp_grid,options);

    result = deepcopy(x1)
    for i=1:length(TD_OP)-1
      copy!(result,x1)
      @test norm(P_sub[i](TD_OP[i]*result)-(TD_OP[i]*result))/norm((TD_OP[i]*result)) <= 1.5*options.feas_tol
    end

#test parallel multilevel PARSDMM
        options=PARSDMM_options()
        default_PARSDMM_options(options,Float64)
        options.parallel = true

            (TD_OP_levels,AtA_levels,P_sub_levels,set_Prop_levels,comp_grid_levels,constraint_level)=setup_multi_level_PARSDMM(m,n_levels,coarsening_factor,comp_grid,constraint,options)

            (x2,log_PARSDMM) = PARSDMM_multi_level(m,TD_OP_levels,AtA_levels,P_sub_levels,set_Prop_levels,comp_grid_levels,options);
            result=Vector{typeof(x[1])}(length(x2))
            for i=1:length(TD_OP_levels[1])-1
              copy!(result,x2)
              @test norm(P_sub_levels[1][i](TD_OP_levels[1][i]*result)-(TD_OP_levels[1][i]*result))/norm((TD_OP_levels[1][i]*result)) <= options.feas_tol
            end
@test isapprox(x1,x2,rtol=options.feas_tol)

#run serial PARSDMM and compare results
default_PARSDMM_options(options,Float64)
options.parallel    = false

(P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid,options.FL);
(TD_OP,AtA,l,y) = PARSDMM_precompute_distribute(TD_OP,set_Prop,comp_grid,options)

(x3,log_PARSDMM) = PARSDMM(m,AtA,TD_OP,set_Prop,P_sub,comp_grid,options);
result=Vector{typeof(x[1])}(length(x3))
for i=1:length(TD_OP)-1
  copy!(result,x3)
  @test norm(P_sub[i](TD_OP[i]*result)-(TD_OP[i]*result))/norm((TD_OP[i]*result)) <= options.feas_tol
end

@test isapprox(x1,x3,rtol=options.feas_tol)

end

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
############################################################################################
### Types
############################################################################################
@doc raw"""
Model implied covariance and means via RAM notation.

# Constructor

    RAM(;
        specification,
        meanstructure = false,
        gradient = true,
        kwargs...)

# Arguments
- `specification`: either a `RAMMatrices` or `ParameterTable` object
- `meanstructure::Bool`: does the model have a meanstructure?
- `gradient::Bool`: is gradient-based optimization used

# Extended help

## Implementation
Subtype of `SemImply`.

## RAM notation

The model implied covariance matrix is computed as
```math
    \Sigma = F(I-A)^{-1}S(I-A)^{-T}F^T
```
and for models with a meanstructure, the model implied means are computed as
```math
    \mu = F(I-A)^{-1}M
```

## Interfaces
- `identifier(::RAM) `-> Dict containing the parameter labels and their position
- `n_par(::RAM)` -> Number of parameters

- `Σ(::RAM)` -> model implied covariance matrix
- `μ(::RAM)` -> model implied mean vector

RAM matrices for the current parameter values:
- `A(::RAM)`
- `S(::RAM)`
- `F(::RAM)`
- `M(::RAM)`

Jacobians of RAM matrices w.r.t to the parameter vector `θ`
- `∇A(::RAM)` -> ``∂vec(A)/∂θᵀ``
- `∇S(::RAM)` -> ``∂vec(S)/∂θᵀ``
- `∇M(::RAM)` = ``∂M/∂θᵀ``

Vector of indices of each parameter in the respective RAM matrix:
- `A_indices(::RAM)`
- `S_indices(::RAM)`
- `M_indices(::RAM)`

Additional interfaces
- `F⨉I_A⁻¹(::RAM)` -> ``F(I-A)^{-1}``
- `F⨉I_A⁻¹S(::RAM)` -> ``F(I-A)^{-1}S``
- `I_A(::RAM)` -> ``I-A``
- `has_meanstructure(::RAM)` -> `Val{Bool}` does the model have a meanstructure?

Only available in gradient! calls:
- `I_A⁻¹(::RAM)` -> ``(I-A)^{-1}``
"""
mutable struct RAM{A1, A2, A3, A4, A5, A6, V, V2, I1, I2, I3, M1, M2, M3, M4, S1, S2, S3, B, D} <: SemImply
    Σ::A1
    A::A2
    S::A3
    F::A4
    μ::A5
    M::A6

    n_par::V
    ram_matrices::V2
    has_meanstructure::B

    A_indices::I1
    S_indices::I2
    M_indices::I3

    F⨉I_A⁻¹::M1
    F⨉I_A⁻¹S::M2
    I_A::M3
    I_A⁻¹::M4

    ∇A::S1
    ∇S::S2
    ∇M::S3

    identifier::D
end

using StructuralEquationModels

############################################################################################
### Constructors
############################################################################################

function RAM(;
        specification,
        #vech = false,
        gradient = true,
        meanstructure = false,
        kwargs...)

    ram_matrices = RAMMatrices(specification)
    identifier = StructuralEquationModels.identifier(ram_matrices)

    # get dimensions of the model
    n_par = length(ram_matrices.parameters)
    n_var, n_nod = ram_matrices.size_F
    parameters = ram_matrices.parameters
    F = zeros(ram_matrices.size_F); F[CartesianIndex.(1:n_var, ram_matrices.F_ind)] .= 1.0

    # get indices
    A_indices = copy(ram_matrices.A_ind)
    S_indices = copy(ram_matrices.S_ind)
    !isnothing(ram_matrices.M_ind) ? M_indices = copy(ram_matrices.M_ind) : M_indices = nothing

    #preallocate arrays
    A_pre = zeros(n_nod, n_nod)
    S_pre = zeros(n_nod, n_nod)
    !isnothing(M_indices) ? M_pre = zeros(n_nod) : M_pre = nothing

    set_RAMConstants!(A_pre, S_pre, M_pre, ram_matrices.constants)

    A_pre = check_acyclic(A_pre, n_par, A_indices)

    # pre-allocate some matrices
    Σ = zeros(n_var, n_var)
    F⨉I_A⁻¹ = zeros(n_var, n_nod)
    F⨉I_A⁻¹S = zeros(n_var, n_nod)
    I_A = similar(A_pre)

    if gradient
        ∇A = get_matrix_derivative(A_indices, parameters, n_nod^2)
        ∇S = get_matrix_derivative(S_indices, parameters, n_nod^2)
    else
        ∇A = nothing
        ∇S = nothing
    end

    # μ
    if meanstructure

        has_meanstructure = Val(true)

        if gradient
            ∇M = get_matrix_derivative(M_indices, parameters, n_nod)
        else
            ∇M = nothing
        end

        μ = zeros(n_var)

    else
        has_meanstructure = Val(false)
        M_indices = nothing
        M_pre = nothing
        μ = nothing
        ∇M = nothing
    end

    return RAM(
        Σ,
        A_pre,
        S_pre,
        F,
        μ,
        M_pre,

        n_par,
        ram_matrices,
        has_meanstructure,

        A_indices,
        S_indices,
        M_indices,

        F⨉I_A⁻¹,
        F⨉I_A⁻¹S,
        I_A,
        copy(I_A),

        ∇A,
        ∇S,
        ∇M,

        identifier
    )
end

############################################################################################
### methods
############################################################################################

# dispatch on meanstructure
objective!(imply::RAM, par, model::AbstractSemSingle) =
    objective!(imply, par, model, imply.has_meanstructure)
gradient!(imply::RAM, par, model::AbstractSemSingle) =
    gradient!(imply, par, model, imply.has_meanstructure)

# objective and gradient
function objective!(imply::RAM, parameters, model, has_meanstructure::Val{T}) where T

    fill_A_S_M(
        imply.A,
        imply.S,
        imply.M,
        imply.A_indices,
        imply.S_indices,
        imply.M_indices,
        parameters)

    imply.I_A .= I - imply.A

    copyto!(imply.F⨉I_A⁻¹, imply.F)
    rdiv!(imply.F⨉I_A⁻¹, factorize(imply.I_A))

    Σ_RAM!(
        imply.Σ,
        imply.F⨉I_A⁻¹,
        imply.S,
        imply.F⨉I_A⁻¹S)

    if T
        μ_RAM!(imply.μ, imply.F⨉I_A⁻¹, imply.M)
    end

end

function gradient!(imply::RAM, parameters, model::AbstractSemSingle, has_meanstructure::Val{T}) where T

    fill_A_S_M(
        imply.A,
        imply.S,
        imply.M,
        imply.A_indices,
        imply.S_indices,
        imply.M_indices,
        parameters)

    imply.I_A .= I - imply.A
    copyto!(imply.I_A⁻¹, imply.I_A)

    imply.I_A⁻¹ .= LinearAlgebra.inv!(factorize(imply.I_A⁻¹))
    imply.F⨉I_A⁻¹ .= imply.F*imply.I_A⁻¹

    Σ_RAM!(
        imply.Σ,
        imply.F⨉I_A⁻¹,
        imply.S,
        imply.F⨉I_A⁻¹S)

    if T
        μ_RAM!(imply.μ, imply.F⨉I_A⁻¹, imply.M)
    end

end

hessian!(imply::RAM, par, model::AbstractSemSingle, has_meanstructure) =
    gradient!(imply, par, model, has_meanstructure)
objective_gradient!(imply::RAM, par, model::AbstractSemSingle, has_meanstructure) =
    gradient!(imply, par, model, has_meanstructure)
objective_hessian!(imply::RAM, par, model::AbstractSemSingle, has_meanstructure) =
    gradient!(imply, par, model, has_meanstructure)
gradient_hessian!(imply::RAM, par, model::AbstractSemSingle, has_meanstructure) =
    gradient!(imply, par, model, has_meanstructure)
objective_gradient_hessian!(imply::RAM, par, model::AbstractSemSingle, has_meanstructure) =
    gradient!(imply, par, model, has_meanstructure)

############################################################################################
### Recommended methods
############################################################################################

identifier(imply::RAM) = imply.identifier
n_par(imply::RAM) = imply.n_par

function update_observed(imply::RAM, observed::SemObserved; kwargs...)
    if n_man(observed) == size(imply.Σ, 1)
        return imply
    else
        return RAM(;observed = observed, kwargs...)
    end
end

############################################################################################
### additional methods
############################################################################################

Σ(imply::RAM) = imply.Σ
μ(imply::RAM) = imply.μ

A(imply::RAM) = imply.A
S(imply::RAM) = imply.S
F(imply::RAM) = imply.F
M(imply::RAM) = imply.M

∇A(imply::RAM) = imply.∇A
∇S(imply::RAM) = imply.∇S
∇M(imply::RAM) = imply.∇M

A_indices(imply::RAM) = imply.A_indices
S_indices(imply::RAM) = imply.S_indices
M_indices(imply::RAM) = imply.M_indices

F⨉I_A⁻¹(imply::RAM) = imply.F⨉I_A⁻¹
F⨉I_A⁻¹S(imply::RAM) = imply.F⨉I_A⁻¹S
I_A(imply::RAM) = imply.I_A
I_A⁻¹(imply::RAM) = imply.I_A⁻¹ # only for gradient available!

has_meanstructure(imply::RAM) = imply.has_meanstructure

ram_matrices(imply::RAM) = imply.ram_matrices

############################################################################################
### additional functions
############################################################################################

function Σ_RAM!(Σ, F⨉I_A⁻¹, S, pre2)
    mul!(pre2, F⨉I_A⁻¹, S)
    mul!(Σ, pre2, F⨉I_A⁻¹')
end

function μ_RAM!(μ, F⨉I_A⁻¹, M)
    mul!(μ, F⨉I_A⁻¹, M)
end

function check_acyclic(A_pre, n_par, A_indices)
    # fill copy of A-matrix with random parameters
    A_rand = copy(A_pre)
    randpar = rand(n_par)

    fill_matrix(
        A_rand,
        A_indices,
        randpar)

    # check if the model is acyclic
    acyclic = isone(det(I-A_rand))

    # check if A is lower or upper triangular
    if iszero(A_rand[.!tril(ones(Bool, size(A_pre)...))])
        A_pre = LowerTriangular(A_pre)
    elseif iszero(A_rand[.!tril(ones(Bool, size(A_pre)...))'])
        A_pre = UpperTriangular(A_pre)
    elseif acyclic
        @info "Your model is acyclic, specifying the A Matrix as either Upper or Lower Triangular can have great performance benefits.\n" maxlog=1
    end

    return A_pre
end


text:
function test_craig()
  craig_tol = 1.0e-6

  function test_craig(A, b; λ=0.0)
    (x, y, stats) = craig(A, b, λ=λ)
    r = b - A * x
    # if λ > 0
    #   s = r / sqrt(λ);
    #   r = r - sqrt(λ) * s;
    # end
    resid = norm(r) / norm(b)
    @printf("craig: residual: %7.1e\n", resid)
    return (x, y, stats, resid)
  end

  # Underdetermined consistent.
  A = rand(10, 25); b = A * ones(25)
  (x, y, stats, resid) = test_craig(A, b)
  @test(norm(x - A' * y) <= craig_tol * norm(x))
  @test(resid <= craig_tol)
  @test(stats.solved)
  (xI, xmin, xmin_norm) = check_min_norm(A, b, x)
  @test(norm(xI - xmin) <= cond(A) * craig_tol * xmin_norm)

  # Underdetermined inconsistent.
  A = ones(10, 25); A[1, :] .= 0; b = ones(10)
  (x, y, stats, resid) = test_craig(A, b)
  # @test(norm(x - A' * y) <= craig_tol * norm(x))
  show(stats)
  @test(stats.inconsistent)

  # Square consistent.
  A = rand(10, 10); b = A * ones(10)
  (x, y, stats, resid) = test_craig(A, b)
  @test(norm(x - A' * y) <= craig_tol * norm(x))
  @test(resid <= craig_tol)
  @test(stats.solved)
  (xI, xmin, xmin_norm) = check_min_norm(A, b, x)
  @test(norm(xI - xmin) <= cond(A) * craig_tol * xmin_norm)

  # Square inconsistent.
  A = ones(10, 10); A[1, :] .= 0; b = ones(10)
  (x, y, stats, resid) = test_craig(A, b)
  # @test(norm(x - A' * y) <= craig_tol * norm(x))
  @test(stats.inconsistent)

  # Overdetermined consistent.
  A = rand(25, 10); b = A * ones(10)
  (x, y, stats, resid) = test_craig(A, b)
  @test(norm(x - A' * y) <= craig_tol * norm(x))
  @test(resid <= craig_tol)
  @test(stats.solved)
  (xI, xmin, xmin_norm) = check_min_norm(A, b, x)
  @test(norm(xI - xmin) <= cond(A) * craig_tol * xmin_norm)

  # Overdetermined inconsistent.
  A = ones(5, 3); A[1, :] .= 0; b = ones(5)
  (x, y, stats, resid) = test_craig(A, b)
  # @test(norm(x - A' * y) <= craig_tol * norm(x))
  @test(stats.inconsistent)

  # With regularization, all systems are underdetermined and consistent.
  # (x, y, stats, resid) = test_craig(A, b, λ=1.0e-3)
  # @test(norm(x - A' * y) <= craig_tol * norm(x))
  # @test(resid <= craig_tol)
  # @test(stats.solved)
  # (xI, xmin, xmin_norm) = check_min_norm(A, b, x, λ=1.0e-3)
  # @test(norm(xI - xmin) <= cond(A) * craig_tol * xmin_norm)

  # Code coverage.
  (x, y, stats) = craig(sparse(A), b, λ=1.0e-3)
  show(stats)

  # Test b == 0
  (x, y, stats) = craig(A, zeros(size(A,1)), λ=1.0e-3)
  @test x == zeros(size(A,2))
  @test y == zeros(size(A,1))
  @test stats.status == "x = 0 is a zero-residual solution"

  # Test integer values
  A = [I; rand(1:10, 2, 3)]
  b = A * ones(Int, 3)
  (x, y, stats) = craig(A, b)
  @test stats.solved
end

test_craig()


text:
using GridapGmsh
using GridapGmsh: gmsh, GmshDiscreteModel

function add_pin!(gmsh, o, r, t, lc)
    factory = gmsh.model.geo

    # inner and outer circle
    l1 = add_circle!(gmsh, o, r, lc)
    l2 = add_circle!(gmsh, o, r + t, lc)

    s1 = factory.addPlaneSurface([l1])
    s2 = factory.addPlaneSurface([l2, l1]) # l1 is a hole

    return s1, s2
end

function add_circle!(gmsh, o, r, lc)
    factory = gmsh.model.geo
    ox, oy = o

    p1 = factory.addPoint(ox, oy, 0, lc)  # center - origin
    p2 = factory.addPoint(ox + r, oy, 0, lc)  # right
    p3 = factory.addPoint(ox, oy + r, 0, lc)  # up
    p4 = factory.addPoint(ox - r, oy, 0, lc)  # left
    p5 = factory.addPoint(ox, oy - r, 0, lc)  # down

    c1 = factory.addCircleArc(p2, p1, p3)
    c2 = factory.addCircleArc(p3, p1, p4)
    c3 = factory.addCircleArc(p4, p1, p5)
    c4 = factory.addCircleArc(p5, p1, p2)

    l1 = factory.addCurveLoop([c1, c2, c3, c4])

    return l1
end

function add_square!(gmsh, s)
    factory = gmsh.model.geo

    p1 = factory.addPoint(0, 0, 0, lc)
    p2 = factory.addPoint(s, 0, 0, lc)
    p3 = factory.addPoint(s, s, 0, lc)
    p4 = factory.addPoint(0, s, 0, lc)

    l1 = factory.addLine(p1, p2)
    l2 = factory.addLine(p2, p3)
    l3 = factory.addLine(p3, p4)
    l4 = factory.addLine(p4, p1)

    cl1 = factory.addCurveLoop([l1, l2, l3, l4])

    return cl1
end

N = 17

# in cm
p = 1.26
r = 0.54
lc = 0.5

gmsh.initialize()
gmsh.model.add("c5g7")

# tags
guidetube = Int32[]
UO₂ = Int32[]
MOX_43 = Int32[]
MOX_7 = Int32[]
MOX_87 = Int32[]
fissionchamber = Int32[]

GT_pos = [( 3, 6), ( 3, 9), (3, 12),
          ( 4, 4), ( 4, 14),
          ( 6, 3), ( 6, 6), ( 6, 9), ( 6, 12), ( 6, 15),
          ( 9, 3), ( 9, 6),          ( 9, 12), ( 9, 15),
          (12, 3), (12, 6), (12, 9), (12, 12), (12, 15),
          (14, 4), (14, 14),
          (15, 6), (15, 9), (15, 12)]

FC_pos = [(9, 9)]

# UO₂ fuel assemblies
for n in 1:2
    # corner
    xc = (n - 1) * N * p
    yc = (n - 1) * N * p

    for i in 1:N, j in 1:N
        xo = xc + p / 2 + (i - 1) * p
        yo = yc + p / 2 + (j - 1) * p

        tag = add_circle!(gmsh, (xo, yo), r, lc)
        gmsh.model.geo.addPlaneSurface([tag])

        pos = (i, j)
        if pos in GT_pos
            push!(guidetube, tag)
        elseif pos in FC_pos
            push!(fissionchamber, tag)
        else
            push!(UO₂, tag)
        end
    end
end

MOX_43_pos = Tuple{Int64,Int64}[]
for i in 1:N, j in 1:N
    if i == 1 || i == N || j == 1 || j == N
        push!(MOX_43_pos, (i, j))
    end
end

MOX_87_pos = Tuple{Int64,Int64}[]
for i in 4:N-3, j in 4:N-3
    pos = (i, j)
    if pos in GT_pos
        continue
    elseif pos in FC_pos
        continue
    elseif pos in ((5, 4), (13, 4), (4, 5), (14, 5), (4, 13), (14, 13), (5, 14), (13, 14))
        continue
    else
        push!(MOX_87_pos, pos)
    end
end

MOX_7_pos = Tuple{Int64,Int64}[]
for i in 1:N, j in 1:N
    pos = (i, j)
    if pos in GT_pos
        continue
    elseif pos in FC_pos
        continue
    elseif pos in MOX_43_pos
        continue
    elseif pos in MOX_87_pos
        continue
    # elseif pos in ((5, 4), (13, 4), (4, 5), (14, 5), (4, 13), (14, 13), (5, 14), (13, 14))
    #     push!(MOX_7_pos, pos)
    else
        push!(MOX_7_pos, pos)
    end
end

# MOX fuel assemblies
for n in 1:2
    # corner
    xc = isone(n) ? N * p : 0
    yc = isone(n) ? 0 : N * p

    for i in 1:N, j in 1:N
        xo = xc + p / 2 + (i - 1) * p
        yo = yc + p / 2 + (j - 1) * p

        tag = add_circle!(gmsh, (xo, yo), r, lc)
        gmsh.model.geo.addPlaneSurface([tag])

        pos = (i, j)
        if pos in GT_pos
            push!(guidetube, tag)
        elseif pos in MOX_43_pos
            push!(MOX_43, tag)
        elseif pos in MOX_7_pos
            push!(MOX_7, tag)
        elseif pos in MOX_87_pos
            push!(MOX_87, tag)
        elseif pos in FC_pos
            push!(fissionchamber, tag)
        end
    end
end

s = 2 * N * p
lc = 0.6
factory = gmsh.model.geo

p1 = factory.addPoint(0, 0, 0, lc)
p2 = factory.addPoint(s, 0, 0, lc)
p3 = factory.addPoint(s, s, 0, lc)
p4 = factory.addPoint(0, s, 0, lc)

ml1 = factory.addLine(p1, p2)
ml2 = factory.addLine(p2, p3)
ml3 = factory.addLine(p3, p4)
ml4 = factory.addLine(p4, p1)

cl1 = factory.addCurveLoop([ml1, ml2, ml3, ml4])

moderator = factory.addPlaneSurface(vcat(cl1, guidetube, fissionchamber, UO₂, MOX_43, MOX_7, MOX_87))

s = 2 * N * p
δs = 1 / 3 * N * p
lc = 1.
factory = gmsh.model.geo

p1 = factory.addPoint(s, 0, 0, lc)
p2 = factory.addPoint(s, s, 0, lc)
p3 = factory.addPoint(0, s, 0, lc)
p4 = factory.addPoint(0, s + δs, 0, lc)
p5 = factory.addPoint(s + δs, s + δs, 0, lc)
p6 = factory.addPoint(s + δs, 0, 0, lc)

r1l1 = factory.addLine(p1, p2)
r1l2 = factory.addLine(p2, p3)
r1l3 = factory.addLine(p3, p4)
r1l4 = factory.addLine(p4, p5)
r1l5 = factory.addLine(p5, p6)
r1l6 = factory.addLine(p6, p1)

cl1 = factory.addCurveLoop([r1l1, r1l2, r1l3, r1l4, r1l5, r1l6])

ref1 = factory.addPlaneSurface([cl1])

s += δs
δs = 1 / 3 * N * p
lc = 1.5
factory = gmsh.model.geo

p1 = factory.addPoint(s, 0, 0, lc)
p2 = factory.addPoint(s, s, 0, lc)
p3 = factory.addPoint(0, s, 0, lc)
p4 = factory.addPoint(0, s + δs, 0, lc)
p5 = factory.addPoint(s + δs, s + δs, 0, lc)
p6 = factory.addPoint(s + δs, 0, 0, lc)

r2l1 = factory.addLine(p1, p2)
r2l2 = factory.addLine(p2, p3)
r2l3 = factory.addLine(p3, p4)
r2l4 = factory.addLine(p4, p5)
r2l5 = factory.addLine(p5, p6)
r2l6 = factory.addLine(p6, p1)

cl1 = factory.addCurveLoop([r2l1, r2l2, r2l3, r2l4, r2l5, r2l6])

ref2 = factory.addPlaneSurface([cl1])

s += δs
δs = 1 / 3 * N * p
lc = 2.
factory = gmsh.model.geo

p1 = factory.addPoint(s, 0, 0, lc)
p2 = factory.addPoint(s, s, 0, lc)
p3 = factory.addPoint(0, s, 0, lc)
p4 = factory.addPoint(0, s + δs, 0, lc)
p5 = factory.addPoint(s + δs, s + δs, 0, lc)
p6 = factory.addPoint(s + δs, 0, 0, lc)

r3l1 = factory.addLine(p1, p2)
r3l2 = factory.addLine(p2, p3)
r3l3 = factory.addLine(p3, p4)
r3l4 = factory.addLine(p4, p5)
r3l5 = factory.addLine(p5, p6)
r3l6 = factory.addLine(p6, p1)

cl1 = factory.addCurveLoop([r3l1, r3l2, r3l3, r3l4, r3l5, r3l6])

ref3 = factory.addPlaneSurface([cl1])

pg1 = gmsh.model.geo.addPhysicalGroup(2, guidetube)
pg2 = gmsh.model.geo.addPhysicalGroup(2, fissionchamber)
pg3 = gmsh.model.geo.addPhysicalGroup(2, UO₂)
pg4 = gmsh.model.geo.addPhysicalGroup(2, MOX_43)
pg5 = gmsh.model.geo.addPhysicalGroup(2, MOX_7)
pg6 = gmsh.model.geo.addPhysicalGroup(2, MOX_87)
pg7 = gmsh.model.geo.addPhysicalGroup(2, [moderator, ref1, ref2, ref3])
pg8 = gmsh.model.geo.addPhysicalGroup(1, [ml1, r1l6, r2l6, r3l6])
pg9 = gmsh.model.geo.addPhysicalGroup(1, [ml4, r1l3, r2l3, r3l3])
pg10 = gmsh.model.geo.addPhysicalGroup(1, [r3l4])
pg11 = gmsh.model.geo.addPhysicalGroup(1, [r3l5])

gmsh.model.setPhysicalName(2, pg1, "guide-tube")
gmsh.model.setPhysicalName(2, pg2, "fission-chamber")
gmsh.model.setPhysicalName(2, pg3, "UO2")
gmsh.model.setPhysicalName(2, pg4, "MOX_43")
gmsh.model.setPhysicalName(2, pg5, "MOX_7")
gmsh.model.setPhysicalName(2, pg6, "MOX_87")
gmsh.model.setPhysicalName(2, pg7, "water")
gmsh.model.setPhysicalName(1, pg8, "bottom")
gmsh.model.setPhysicalName(1, pg9, "right")
gmsh.model.setPhysicalName(1, pg10, "top")
gmsh.model.setPhysicalName(1, pg11, "left")

# removemos puntos de la geometria duplicados (origenes por ejemplo)
gmsh.model.geo.removeAllDuplicates()

gmsh.model.geo.synchronize()

gmsh.model.mesh.generate(2)

gmsh.write("c5g7-2.msh")

if !("-nopopup" in ARGS)
    gmsh.fltk.run()
end

gmsh.finalize()

# move to json file format
using Gridap
mshfile = joinpath(@__DIR__,"../c5g7-2.msh")
model = GmshDiscreteModel(mshfile; renumber=true)
Gridap.Io.to_json_file(model, "c5g7-2.json")


text:
arg_dict = Dict{String,Any}(
"physics" => "Advection",
"run_type" => 20,
"order" => 1,
"dimensions" => 2,
"smb_name" => "SRCMESHES/square_x0-2_y0-2_20x20_tri.smb",
"IC_name" => "ICsinwave",
"numBC" => 1,
"BC1" => [0, 1, 2, 3],
# "BC1" => [0],
"BC1_name" => "sinwaveBC",
"delta_t" => 0.1,
"t_max" => 1.0,
# "res_abstol" => -1.0,
# "res_reltol" => -1.0,
"use_DG" => true,
"operator_type" => "SBPGamma",
# "Flux_name" => "RoeFlux",
"Flux_name" => "LFFlux",
"output_freq" => 100,
"real_time" => true,
"vortex_x0" => 5.0,
"use_itermax" => true,
"itermax" => 1000,
"output_freq" => 1,
"jac_method" => 1,        # 1: FD
"jac_type" => 2,          # 1: dense matrix, 2: Julia sparse
"do_postproc" => true,
"exact_soln_func" => "ICsinwave",
"use_Minv" => true,
)


text:
using ApplicationBuilderAppUtils

# Create a temporary .html file, and open it to share the greetings.
Base.@ccallable function julia_main(ARGS::Vector{String})::Cint
    tmpdir = mktempdir()
    filename = joinpath(tmpdir, "hello.html")
    open(filename, "w") do io
        println(io, "Hello, World!")
        println(io, "
    -- Love, $PROGRAM_FILE")
        println(io, """

Current working directory: $(pwd())""")
        ApplicationBuilderAppUtils.cd_to_bundle_resources()
        println(io, """

After cd to resources directory: $(pwd())""")
    end
    run(`open file://$filename`)
    return 0
end


text:
using Base.Test

using DataReaders

@testset "google_quotes" begin
    dr = DataReader("google-quotes")

    symb = DataSymbol("MSFT")

    data = get(dr, symb)
end


text:
_inner_type(::SimpleSDMResponse{T}) where {T <: Any} = T
_inner_type(::SimpleSDMPredictor{T}) where {T <: Any} = T

"""
    mask!(l1::T1, l2::T2) where {T1 <: SimpleSDMLayer, T2 <: SimpleSDMLayer}

Changes the second layer so that the positions for which the first layer is zero
(of the appropriate type) or `nothing` are set to `nothing`. This is mostly
useful in cases where you have a `Bool` layer.
"""
function mask!(l1::T1, l2::T2) where {T1 <: SimpleSDMLayer, T2 <: SimpleSDMLayer}
    _itype = _inner_type(l1)
    dropfunc = (x) -> isnothing(x) || (x == zero(_itype))
    todrop = findall(dropfunc, l1.grid)
    l2.grid[todrop] .= nothing
end

"""
    mask(l1::T1, l2::T2) where {T1 <: SimpleSDMLayer, T2 <: SimpleSDMLayer}

Returns a copy of the second layer in which the positions for which the first
layer is zero (of the appropriate type) or `nothing` are set to `nothing`. This
is mostly useful in cases where you have a `Bool` layer.
"""
function mask(l1::T1, l2::T2) where {T1 <: SimpleSDMLayer, T2 <: SimpleSDMLayer}
    l3 = copy(l2)
    mask!(l1, l3)
    return l3
end

_format_polygon(xy) = Polygon([Point(c[1], c[2]) for c in xy[1]])

function _points_from_polygon(poly)
    pl = Point[]
    for i in 1:length(poly.exterior)
        push!(pl, poly.exterior[i][1])
    end
    push!(pl, pl[1])
    return pl
end

function _mask_from_polygon(polygon::Polygon, layer::T) where {T <: SimpleSDMLayer}
    msk = similar(layer, Bool)
    pts = _points_from_polygon(polygon)
    keep = filter(p -> inpolygon(Array(p), pts)!=0, keys(layer))
    if length(keep) > 0
        msk[keep] = fill(true, length(keep))
    end
    return msk
end

function mask(multipolygon::P, layer::T) where {T <: SimpleSDMLayer, P <: Vector{<:GeometryBasics.Polygon}}
    msks = [_mask_from_polygon(polygon, layer) for polygon in multipolygon]
    return mask(maximum(msks), layer)
end

function mask(polygon::Polygon, layer::T) where {T <: SimpleSDMLayer}
    msk = _mask_from_polygon(polygon, layer)
    return mask(msk, layer)
end

function mask(circle::Circle, layer::T) where {T <: SimpleSDMLayer}
    msk = similar(layer, Bool)
    keep = filter(k -> k in circle, keys(layer))
    if length(keep) > 0
        msk[keep] = fill(true, length(keep))
    end
    return mask(msk, layer)
end


text:
abstract type TriangularHarmonicPlan{T} end

function *(P::TriangularHarmonicPlan, X::AbstractMatrix)
    mul!(zero(X), P, X)
end

function \(P::TriangularHarmonicPlan, X::AbstractMatrix)
    At_mul_B!(zero(X), P, X)
end

include("trifunctions.jl")
include("slowplan.jl")

function plan_tri2cheb(A::AbstractMatrix, α, β, γ; opts...)
    M, N = size(A)
    # if M ≤ 1023
        SlowTriangularHarmonicPlan(A, α, β, γ)
    # else
    #     ThinTriangularHarmonicPlan(A, α, β, γ; opts...)
    # end
end

tri2cheb(A::AbstractMatrix, α, β, γ; opts...) = plan_tri2cheb(A, α, β, γ; opts...)*A
cheb2tri(A::AbstractMatrix, α, β, γ; opts...) = plan_tri2cheb(A, α, β, γ; opts...)\A


text:
"""
	delete_arms(config::REDCap.Config, arms::Array)

Delete Arms from project. Removing all arms reverts the project into a non-longitudinal project.

#### NOTE: This only works for longitudinal projects.

#### Parameters:
* `config` - Struct containing url and api-key
* `arms` - Array of arm names to delete

#### Returns:
Number of succesfully deleted arms
"""
function delete_arms(config::REDCap.Config, arms::Array)
	return api_pusher("delete", "arm", config, arms=arms)
end

"""
	delete_events(config::REDCap.Config, events::Array)

Delete Events from project. Removing all but one event reverts the project into a non-longitudinal project.

#### NOTE: This only works for longitudinal projects.

#### Parameters:
* `config` - Struct containing url and api-key
* `events` - Array of event names to delete

#### Returns:
Number of successfully deleted events
"""
function delete_events(config::REDCap.Config, events::Array)
	return api_pusher("delete", "event", config, events=events)
end

"""
	delete_file(config::REDCap.Config, record::String, field::String, event::String; repeat_instance::Integer=1, returnFormat::String="json")

Delete document attached to record.

#### Parameters:
* `config` - Struct containing url and api-key
* `record` - Name of record containing file
* `field` - Name of field containing file
* `event` - Name of event containing file
* `repeat_instance` - Number of repeated instances (long project)
* `returnFormat` - Error message format

#### Returns:
Nothing/error
"""
function delete_file(config::REDCap.Config, record::String, field::String, event::String; repeat_instance::Integer=1, returnFormat::String="json")
	return api_pusher("delete", "file", config, record=record, field=field, event=event, repeat_instance=repeat_instance, returnFormat=returnFormat)
end

"""
	delete_records(config::REDCap.Config, records::Array; arm::Integer=0)

Delete one or more records from project.

#### Parameters:
* `config` - Struct containing the url and api-key
* `records` - Array of record names to delete
* `arm` - Number of arm containing records

#### Returns:
Number of records successfully deleted
"""
function delete_records(config::REDCap.Config, records::Array; arm::Integer=-1)
	if arm != -1	#REDCap treats the request differently if arm is sent as a field
		return api_pusher("delete", "record", config, records=records, arm=arm)
	else
		return api_pusher("delete", "record", config, records=records)
	end
end


text:
import ..UncertainDatasets:
    UVAL_COLLECTION_TYPES

import ..SamplingConstraints:
    SamplingConstraint

import ..UVAL_COLLECTION_TYPES

"""
    resample_elwise(uvd::UVAL_COLLECTION_TYPES)

Resample each element in `uvals` once. The i-th entry in the returned
vector is a `1`-element vector containing one unique draw of `uvals[i]`.
"""
function resample_elwise(uvd::UVAL_COLLECTION_TYPES)
    [resample(uvd[i], 1) for i = 1:length(uvd)]
end

"""
    resample_elwise(uvd::UVAL_COLLECTION_TYPES, n::Int)

Resample each element in `uvals` `n` times. The i-th entry in the returned
vector is a `n`-element vector consisting of `n` unique draws of `uvals[i]`.
"""
function resample_elwise(uvd::UVAL_COLLECTION_TYPES, n::Int)
    [resample(uvd[i], n) for i = 1:length(uvd)]
end

"""
    resample_elwise(uvd::UVAL_COLLECTION_TYPES,
        constraint::Union{SamplingConstraint, Vector{SamplingConstraint}},
        n::Int)

Resample each element in `uvals` `n` times. The i-th entry in the returned
vector is a `n`-element vector consisting of `n` unique draws of `uvals[i]`, drawn
after first truncating the support of `uvals[i]` according to the provided `constraint`(s).
"""
resample_elwise(uvd::UVAL_COLLECTION_TYPES,
    constraint::Union{SamplingConstraint, Vector{SamplingConstraint}},
    n::Int)

function resample_elwise(uvd::UVAL_COLLECTION_TYPES,
        constraint::SamplingConstraint, n::Int)

    [resample(uvd[i], constraint, n) for i = 1:length(uvd)]
end

function resample_elwise(uvd::UVAL_COLLECTION_TYPES,
        constraint::Vector{<:SamplingConstraint},
        n::Int)

    Lc, Luv = length(constraint), length(uvd)

    if Lc != Luv
        error("""The number of constraints must match the number of uncertain values in the dataset. Got $Lc constraint but needed $Luv.""")
    end

    [resample(uvd[i], constraint[i], n) for i = 1:length(uvd)]
end


text:
function calc_forces(surf::TwoDSurf, vels::Vector{Float64})

    # First term in eqn (2.30) Ramesh et al. in coefficient form
    cnc = 2*pi*((surf.kinem.u + vels[1])*cos(surf.kinem.alpha)/surf.uref + (surf.kinem.hdot - vels[2])*sin(surf.kinem.alpha)/surf.uref)*(surf.a0[1] + surf.aterm[1]/2.)

    # Second term in eqn (2.30) Ramesh et al. in coefficient form
    cnnc = 2*pi*(3*surf.c*surf.a0dot[1]/(4*surf.uref) + surf.c*surf.adot[1]/(4*surf.uref) + surf.c*surf.adot[2]/(8*surf.uref))

    # Suction force given in eqn (2.31) Ramesh et al.
    cs = 2*pi*surf.a0[1]*surf.a0[1]

    #The components of normal force and moment from induced velocities are calulcated in dimensional units and nondimensionalized later
    nonl=0
    nonl_m=0
    for ib = 1:surf.ndiv-1
        nonl = nonl + (surf.uind[ib]*cos(surf.kinem.alpha) - surf.wind[ib]*sin(surf.kinem.alpha))*surf.bv[ib].s
        nonl_m = nonl_m + (surf.uind[ib]*cos(surf.kinem.alpha) - surf.wind[ib]*sin(surf.kinem.alpha))*surf.x[ib]*surf.bv[ib].s
    end
    nonl = nonl*2. /(surf.uref*surf.uref*surf.c)
    nonl_m = nonl_m*2. /(surf.uref*surf.uref*surf.c*surf.c)

    # Normal force coefficient
    cn = cnc + cnnc + nonl

    # Lift and drag coefficients
    cl = cn*cos(surf.kinem.alpha) + cs*sin(surf.kinem.alpha)
    cd = cn*sin(surf.kinem.alpha)-cs*cos(surf.kinem.alpha)

    #Pitching moment is clockwise or nose up positive
    cm = cn*surf.pvt - 2*pi*(((surf.kinem.u + vels[1])*cos(surf.kinem.alpha)/surf.uref + (surf.kinem.hdot - vels[2])*sin(surf.kinem.alpha)/surf.uref)*(surf.a0[1]/4. + surf.aterm[1]/4. - surf.aterm[2]/8.) + (surf.c/surf.uref)*(7. *surf.a0dot[1]/16. + 3. *surf.adot[1]/16. + surf.adot[2]/16. - surf.adot[3]/64.)) - nonl_m
    return cl, cd, cm
end

function writeStamp(dirname::String, t::Float64, surf::TwoDSurf, curfield::TwoDFlowField)
    dirvec = readdir()
    if dirname in dirvec
        rm(dirname, recursive=true)
    end
    mkdir(dirname)
    cd(dirname)

    f = open("timeKinem", "w")
    Serialization.serialize(f, ["#time \t", "alpha (deg) \t", "h/c \t", "u/uref \t"])
    DelimitedFiles.writedlm(f, [t, surf.kinem.alpha, surf.kinem.h, surf.kinem.u]')
    close(f)

    f = open("FourierCoeffs", "w")
    Serialization.serialize(f, ["#Fourier coeffs (0-n) \t", "d/dt (Fourier coeffs) \n"])
    matfour = zeros(surf.naterm+1, 2)
    matfour[:,1] = [surf.a0[1];surf.aterm[:]]
    matfour[:,2] = [surf.a0dot[1];surf.adot[:]]
    DelimitedFiles.writedlm(f, matfour)
    close(f)

    # f = open("forces", "w")
    # write(f, ["#cl \t", "cd \t", "cm \t", "Gamma \t", "cn \t", "cs \t", "cnc \t",
    #           "cnnc \t", "nonl \t", "cm_n \t", "cm_pvt \t", "nonl_m \n"])
    # cl, cd, cm, gamma, cn, cs, cnc, cncc, nonl, cm_n, cm_pvt, nonl_m = calc_forces_more(surf)
    # writedlm(f, [cl, cd, cm, gamma, cn, cs, cnc, cnnc, nonl, cm_n, cm_pvt, nonl_m])
    # close(f)

    f = open("tev", "w")
    Serialization.serialize(f, ["#strength \t", "x-position \t", "z-position \n"])
    tevmat = zeros(length(curfield.tev), 3)
    for i = 1:length(curfield.tev)
        tevmat[i,:] = [curfield.tev[i].s curfield.tev[i].x curfield.tev[i].z]
    end
    DelimitedFiles.writedlm(f, tevmat)
    close(f)

    f = open("lev", "w")
    Serialization.serialize(f, ["#strength \t", "x-position \t", "z-position \n"])
    levmat = zeros(length(curfield.lev), 3)
    for i = 1:length(curfield.lev)
        levmat[i,:] = [curfield.lev[i].s curfield.lev[i].x curfield.lev[i].z]
    end
    DelimitedFiles.writedlm(f, levmat)
    close(f)

    f = open("boundv", "w")
    Serialization.serialize(f, ["#strength \t", "x-position \t", "z-position \n"])
    bvmat = zeros(length(surf.bv), 3)
    for i = 1:length(surf.bv)
        bvmat[i,:] = [surf.bv[i].s surf.bv[i].x surf.bv[i].z]
    end
    DelimitedFiles.writedlm(f, bvmat)
    close(f)

    f = open("delcp_edgevel", "w")
    Serialization.serialize(f, ["#x \t", "delcp \t", "delcp_inner \t", "delcp_outer \t", "qu \t", "ql \n"])
    mat = zeros(length(surf.x), 6)
    mat[:,1] = surf.x[:]
    delcp, delcp_in, delcp_out = calc_delcp(surf, [curfield.u[1]; curfield.w[1]])
    qu, ql = calc_edgeVel(surf, [curfield.u[1]; curfield.w[1]])
    mat[:,2] = delcp[:]
    mat[:,3] = delcp_in[:]
    mat[:,4] = delcp_out[:]
    mat[:,5] = qu[:]
    mat[:,6] = ql[:]
    DelimitedFiles.writedlm(f, mat)
    close(f)

    cd("..")
end

function writeStamp(dirname::String, t::Float64, surf::Vector{TwoDSurf}, curfield::TwoDFlowField)
    dirvec = readdir()
    if dirname in dirvec
        rm(dirname, recursive=true)
    end
    mkdir(dirname)
    cd(dirname)

    nsurf = length(surf)
    for i = 1:nsurf
        f = open("timeKinem-$i", "w")
        Serialization.Serialization.serialize(f, ["#time \t", "alpha (deg) \t", "h/c \t", "u/uref \t"])
        DelimitedFiles.writedlm(f, [t, surf[i].kinem.alpha, surf[i].kinem.h, surf[i].kinem.u]')
        close(f)

        f = open("FourierCoeffs-$i", "w")
        Serialization.serialize(f, ["#Fourier coeffs (0-n) \t", "d/dt (Fourier coeffs) \n"])
        matfour = zeros(surf[i].naterm+1, 2)
        matfour[:,1] = [surf[i].a0[1];surf[i].aterm[:]]
        matfour[:,2] = [surf[i].a0dot[1];surf[i].adot[:]]
        DelimitedFiles.writedlm(f, matfour)
        close(f)

        # f = open("forces", "w")
        # write(f, ["#cl \t", "cd \t", "cm \t", "Gamma \t", "cn \t", "cs \t", "cnc \t",
        #           "cnnc \t", "nonl \t", "cm_n \t", "cm_pvt \t", "nonl_m \n"])
        # cl, cd, cm, gamma, cn, cs, cnc, cncc, nonl, cm_n, cm_pvt, nonl_m = calc_forces_more(surf)
        # writedlm(f, [cl, cd, cm, gamma, cn, cs, cnc, cnnc, nonl, cm_n, cm_pvt, nonl_m])
        # close(f)
    end

    f = open("tev", "w")
    Serialization.serialize(f, ["#strength \t", "x-position \t", "z-position \n"])
    tevmat = zeros(length(curfield.tev), 3)
    for i = 1:length(curfield.tev)
        tevmat[i,:] = [curfield.tev[i].s curfield.tev[i].x curfield.tev[i].z]
    end
    DelimitedFiles.writedlm(f, tevmat)
    close(f)

    f = open("lev", "w")
    Serialization.serialize(f, ["#strength \t", "x-position \t", "z-position \n"])
    levmat = zeros(length(curfield.lev), 3)
    levcount = 0
    for i = 1:length(curfield.lev)
        if curfield.lev[i].vc != 0.
            levcount += 1
            levmat[levcount,:] = [curfield.lev[i].s curfield.lev[i].x curfield.lev[i].z]
        end
    end

    levmat = levmat[1:levcount,:]
    DelimitedFiles.writedlm(f, levmat)
    close(f)

    for is = 1:nsurf
        f = open("boundv-$is", "w")
        Serialization.serialize(f, ["#strength \t", "x-position \t", "z-position \n"])
        bvmat = zeros(length(surf[is].bv), 3)
        for i = 1:length(surf[is].bv)
            bvmat[i,:] = [surf[is].bv[i].s surf[is].bv[i].x surf[is].bv[i].z]
        end
        DelimitedFiles.writedlm(f, bvmat)
        close(f)
    end

    cd("..")
end

function calc_delcp(surf::TwoDSurf, vels::Vector{Float64})

    p_in = zeros(surf.ndiv)
    p_out = zeros(surf.ndiv)
    gam = zeros(surf.ndiv)
    gamint = zeros(surf.ndiv)
    p_com = zeros(surf.ndiv)
    gammod = zeros(surf.ndiv)

    for i = 1:surf.ndiv
        udash = surf.kinem.u*cos(surf.kinem.alpha) + surf.kinem.hdot*sin(surf.kinem.alpha) + surf.uind[i]*cos(surf.kinem.alpha) - surf.wind[i]*sin(surf.kinem.alpha)

        p_in[i] = 8*surf.uref*sqrt(surf.c)*surf.a0[1]*sqrt(surf.x[i])*udash/(surf.rho*surf.c + 2*surf.x[i]) + 4*surf.uref*sqrt(surf.c)*surf.a0dot[1]*sqrt(surf.x[i])
    end
    for i = 2:surf.ndiv
        udash = surf.kinem.u*cos(surf.kinem.alpha) + surf.kinem.hdot*sin(surf.kinem.alpha) + surf.uind[i]*cos(surf.kinem.alpha) - surf.wind[i]*sin(surf.kinem.alpha)
        gam[i] = surf.a0[1]*cot(surf.theta[i]/2)
        gamint[i] = surf.a0dot[1]*(surf.theta[i] + sin(surf.theta[i]))
        for n = 1:surf.naterm
            gam[i] += surf.aterm[n]*sin(n*surf.theta[i])
            gammod[i] += surf.aterm[n]*sin(n*surf.theta[i])
            if n == 1
                gamint[i] += surf.adot[1]*(surf.theta[i]/2 - sin(2*surf.theta[i])/4)
            else
                gamint[i] += surf.adot[n]/2*(sin((n-1)*surf.theta[i])/(n-1) - sin((n+1)*surf.theta[i])/(n+1))
            end
        end
        gam[i] = gam[i]*surf.uref
        gammod[i] = gammod[i]*surf.uref
        gamint[i] = gamint[i]*surf.uref*surf.c

        p_out[i] = 2*udash*gam[i] + gamint[i]
        p_com[i] = 2*udash*surf.uref*surf.a0[1]*(cos(surf.theta[i]/2)-1)/sin(surf.theta[i]/2) + 2*udash*gammod[i] + gamint[i] + 8*surf.uref*udash*surf.c*surf.a0[1]*sin(surf.theta[i]/2)/(surf.rho*surf.c + 2*surf.c*(sin(surf.theta[i]/2))^2)
    end

    return p_com, p_in, p_out
end

function calc_edgeVel(surf::TwoDSurf, vels::Vector{Float64})

    gammod = zeros(surf.ndiv)
    q_com_u = zeros(surf.ndiv)
    q_com_l = zeros(surf.ndiv)

    q_com_u[1] = surf.uref*surf.a0[1]/sqrt(0.5*surf.rho)
    q_com_l[1] = surf.uref*surf.a0[1]/sqrt(0.5*surf.rho)

    for i = 2:surf.ndiv
        udash = (surf.kinem.u + vels[1])*cos(surf.kinem.alpha) + (surf.kinem.hdot - vels[2])*sin(surf.kinem.alpha) + surf.uind[i]*cos(surf.kinem.alpha) - surf.wind[i]*sin(surf.kinem.alpha)

        for n = 1:surf.naterm
            gammod[i] += surf.aterm[n]*sin(n*surf.theta[i])
        end
        gammod[i] = gammod[i]*surf.uref

        q_com_u[i] = surf.uref*surf.a0[1]*(cos(surf.theta[i]/2) - 1)/sin(surf.theta[i]/2) + surf.uref*gammod[i] + (sqrt(surf.c)*surf.uref*surf.a0[1] + sqrt(surf.x[i])*udash)/sqrt(surf.x[i] + surf.rho*surf.c/2)
        q_com_l[i] = -surf.uref*surf.a0[1]*(cos(surf.theta[i]/2) - 1)/sin(surf.theta[i]/2) - surf.uref*gammod[i] + (-sqrt(surf.c)*surf.uref*surf.a0[1] + sqrt(surf.x[i])*udash)/sqrt(surf.x[i] + surf.rho*surf.c/2)
    end

    return q_com_u, q_com_l
end

function getEndCycle(mat, k)

    T = pi/k

    end_cycle = mat[end,1]

    #Find number of cycles
    ncyc = 0
    for i = 1:1000
        t_l = real(i)*T
        if t_l > end_cycle
            break
        end
        ncyc = ncyc + 1
    end

    start_t = real(ncyc-1)*T
    end_t = real(ncyc)*T
    start_ind = argmin(abs.(mat[:,1] .- start_t))
    end_ind = argmin(abs.(mat[:,1] .- end_t))

    nlast = end_ind - start_ind + 1

    newmat = zeros(nlast, 8)
    newmat[:,1] = (mat[start_ind:end_ind,1] .- start_t)/T
    for i = 2:7
        newmat[:,i] = mat[start_ind:end_ind,i]
    end

    f = open("resultsSummaryEndCycle", "w")
    Serialization.serialize(f, ["#t/T \t", "alpha (deg) \t", "h/c \t", "u/uref \t", "A0 \t", "Cl \t", "Cd \t", "Cm \n"])
    DelimitedFiles.writedlm(f, newmat)
    close(f)

    return newmat
end

max_stars_repo_name: hjabird/UNSflow

max_issues_repo_name: hjabird/UNSflow

max_forks_repo_name: hjabird/UNSflow


text:
using Base: splat
using Combinatorics: combinations, permutations
using LinearAlgebra: I, ×

function parse_input(io::IO)
    scanners = Vector{Vector{Int}}[]
    while !eof(io)
        line = readline(io)
        m = match(r"^--- scanner ([0-9]+) ---$", line)
        @assert m !== nothing
        scanner = parse(Int, m.captures[1])
        beacons = Vector{Int}[]
        line = readline(io)
        while !isempty(line)
            push!(beacons, parse.(Int, split(line, ',')))
            line = readline(io)
        end
        @assert scanner == length(scanners)
        push!(scanners, beacons)
    end
    return scanners
end

const PERMUTATIONS = [Matrix{Int}(I(3)[:, xyz]) for xyz in permutations(1:3)]

const TRANSFORMATIONS = let vs = Vector.(eachcol([I(3) -I(3)]))
    [[x y (x × y)] for x in vs for y in vs if !iszero(x × y)]
end

function distances(beacons::AbstractVector)::Vector
    return map(combinations(beacons, 2)) do (a, b)
        return abs.(a .- b)
    end
end

function overlap(scanner1::AbstractVector,
                 scanner2::AbstractVector;
                 min_overlap::Integer = 66
                )::Bool
    dists1 = distances(scanner1)
    dists2 = distances(scanner2)
    olap = map(PERMUTATIONS) do P
        intersect(dists1, [P * d for d in dists2]) |> length
    end |> maximum
    return olap >= min_overlap
end

function overlaps(scanners::AbstractVector)::Matrix
    return [overlap(s1, s2) for s1 in scanners, s2 in scanners]
end

struct Transform
    T
    d
end

(t::Transform)(v) = t.T * v + t.d

function find_transform(scanner1::AbstractVector, scanner2::AbstractVector)
    ((olap, d), T) = map(TRANSFORMATIONS) do T
        sc2 = [T * beacon for beacon in scanner2]
        return map(scanner1) do b1
            displacements = [b1 .- b2 for b2 in sc2]
            (olap, d) = map(displacements) do d
                return length(intersect(scanner1, [b2 .+ d for b2 in sc2]))
            end |> findmax
            return (olap, displacements[d])
        end |> maximum
    end |> findmax
    @assert olap >= 12
    return Transform(TRANSFORMATIONS[T], d)
end

function find_transforms(scanners::AbstractVector)
    olaps = overlaps(scanners)
    transforms = Vector{Any}(nothing, length(scanners))
    transforms[1] = identity
    unaligned = collect(2:length(scanners))
    while !isempty(unaligned)
        aligned = setdiff(keys(scanners), unaligned)
        (a, u) = findfirst(olaps[aligned, unaligned]).I
        (i1, i2) = (aligned[a], unaligned[u])
        T = find_transform(scanners[i1], scanners[i2])
        transforms[i2] = transforms[i1] ∘ T
        setdiff!(unaligned, i2)
    end
    return transforms
end

function q1(scanners::AbstractVector, transforms::AbstractVector)
    return union(map((T, sc) -> T.(sc), transforms, scanners)...) |> length
end

function q2(transforms::AbstractVector)
    scanners = map(T -> T([0, 0, 0]), transforms)
    return map(combinations(scanners, 2)) do (sc1, sc2)
        return abs.(sc1 .- sc2) |> sum
    end |> maximum
end

open("day19-test.in", "r") do io
    scanners = parse_input(io)
    transforms = find_transforms(scanners)

    @assert q1(scanners, transforms) == 79
    @assert q2(transforms) == 3621
end

open("day19.in", "r") do io
    scanners = parse_input(io)
    transforms = find_transforms(scanners)

    println("Q1: ", q1(scanners, transforms))
    println("Q2: ", q2(transforms))
end

max_stars_repo_name: dhanak/competitive-coding

max_issues_repo_name: dhanak/competitive-coding

max_forks_repo_name: dhanak/competitive-coding