text:
mutable struct IdealActuatorCommand
    force_B::Vector{Float64}
    torque_B::Vector{Float64}
end

"""
    IdealActuator

This actuator provides exactly effects (e.g., force and torque) that are
requested from it.

"""
function IdealActuator()

    # Initialize the state based on the default (or overwritten) state,
    # constants, and true situation.
    # function init(t, constants, state, draws)
    #     return state
    # end

    # Outputs whatever effects were commanded. This runs at the minor time step.
    function effects(t, constants, state, draws, implicit, command, effects, effects_bus)
        return (BodyForce(command.force_B, [0.; 0.; 0.]),
                BodyTorque(command.torque_B))
    end

    # Create the IdealActuator as a DynamicalModel.
    DynamicalModel("ideal_actuator",
                   #init = init,
                   effects = effects,
                   timing = ModelTiming(0.01),
                   inputs = IdealActuatorCommand([0.; 0.; 0.], [0.; 0.; 0.]))

end

"""
   ReactionWheel

TODO: This should really be a 0-momentum reaction wheel array, because it the
wheels speeds don't affect the body yet. We could make a MomentumEffect that the
body dynamics could consume that would allow it to generically account for
reaction wheels cleanly.

This is a completely basic reaction wheel model. You request a torque, and it
spins the wheel with the opposite torque, providing the requested torque as
an Effect.
"""
function ReactionWheel()

    function derivatives(t, constants, speed, draws, implicit, τ_br_cmd, effects, effects_bus)
        -τ_br_cmd / constants.I_r_G[1]
    end

    function effects(t, constants, speed, draws, implicit, τ_br_cmd, effects, effects_bus)
        body, found  = find_effect(effects, BodyStateEffect())
        if !found
            error("Couldn't find BodyStateEffect.")
        end
        ω_BI_G = qrot(constants.q_GB, body.ω_BI_B)
        τ_br_G = [τ_br_cmd; 0.; 0.] - ω_BI_G × ([constants.I_r_G[1] * speed; 0.; 0.] + constants.I_r_G .* ω_BI_G)
        τ_br_B = qrot(qinv(constants.q_GB), τ_br_G)
        (BodyTorque(τ_br_B),)
    end

    function sense(t, constants, speed, draws, τ_br_cmd, effects, effects_bus)
        (speed, speed, true)
    end

    DynamicalModel("reaction_wheel",
                   constants = ReactionWheelConstants([1.; 1/sqrt(2.); 1/sqrt(2.)], [0.; 0.; 0.; 1.]),
                   state = 100. * 2pi/60., # Wheel speed (rad/s)
                   inputs = 0.,             # Requested torque (Nm)
                   outputs = 100. * 2pi/60., # Measured wheel speed (rad/s)
                   effects = effects,
                   derivatives = derivatives,
                   update = sense,
                   timing = ModelTiming(0.02))

end

mutable struct ReactionWheelConstants
    I_r_G::Vector{Float64}
    q_GB::Vector{Float64}
end

"""
   CoupledReactionWheel

This reaction wheel is modular. The body dynamics don't need to know about it at
all. It uses the differential algebraic equation framework of the simulation
environment so that we don't need to implement the coupled rotor-and-body
dynamics explicitly.

"""
function CoupledReactionWheel()

    # Get the torque on the rotor, coming from motor torque and constraints.
    function torque(τ_br_cmd, constants, state, implicit)
        τ_rb_G = [-τ_br_cmd; implicit[1]; implicit[2]]
        return τ_rb_G
    end

    # Get the derivative of the state, treating the rotor as a free 3DOF body.
    function derivatives(t, constants, state, draws, implicit, τ_br_cmd, effects, effects_bus)
        body,      = find_effect(effects, BodyStateEffect())
        ω_BI_G     = qrot(constants.q_GB, body.ω_BI_B)
        τ_rb_G     = torque(τ_br_cmd, constants, state, implicit)
        ω_RI_G_dot = constants.I_r_G .\ (τ_rb_G - ω_BI_G × (constants.I_r_G .* state.ω_RI_G))
        return CoupledReactionWheelState(ω_RI_G_dot)
    end

    # Ensure that the implicit torque keeps the rotation rates of the rotor wrt
    # the body entirely around the rotor axis (no off-axis rotation of one wrt
    # the other).
    function constraints(t, constants, state, draws, implicit, τ_br_cmd, effects, effects_bus)
        body,  = find_effect(effects, BodyStateEffect())
        ω_BI_G = qrot(constants.q_GB, body.ω_BI_B)
        return state.ω_RI_G[2:3] - ω_BI_G[2:3]
    end

    # The torque on the rotor has an equal and opposite torque on the body.
    function effects(t, constants, state, draws, implicit, τ_br_cmd, effects, effects_bus)
        τ_rb_G = torque(τ_br_cmd, constants, state, implicit)
        τ_br_B = qrot(qinv(constants.q_GB), -τ_rb_G)
        (BodyTorque(τ_br_B),)
    end

    # Produce the measurement (the rotor speed) and also reset the rotation to
    # be strictly on-axis (removes buildup of constraint torque noise)
    function update(t, constants, state, draws, τ_br_cmd, effects, effects_bus)
        body, = find_effect(effects, BodyStateEffect())
        ω_BI_G = qrot(constants.q_GB, body.ω_BI_B)
        state.ω_RI_G[2:3] = ω_BI_G[2:3] # Renormalize.
        (state, state.ω_RI_G[1] - ω_BI_G[1], true)
    end

    DynamicalModel("coupled_reaction_wheel",
                   constants = CoupledReactionWheelConstants([1.; 1/sqrt(2.); 1/sqrt(2.)], [0.; 0.; 0.; 1.]),
                   state = CoupledReactionWheelState([1000. * 2pi/60.; 0.; 0.]),
                   inputs = 0., # Requested torque on rotor from body about 1 axis of G frame (Nm)
                   outputs = 1000. * 2pi/60.,
                   implicit = [0.; 0.], # Implicit variables (constraint torque) (Nm)
                   effects = effects,
                   derivatives = derivatives,
                   constraints = constraints,
                   init = update,
                   update = update,
                   timing = ModelTiming(0.02))

end

mutable struct CoupledReactionWheelConstants
    I_r_G::Vector{Float64}
    q_GB::Vector{Float64}
end

mutable struct CoupledReactionWheelState <: ModelState
    ω_RI_G::Vector{Float64}
end


text:
module LTVModels
using LinearAlgebra, Statistics, Printf, Random, Distributions
using LinearTimeVaryingModelsBase, ControlSystemIdentification
import LinearTimeVaryingModelsBase: AbstractModel, AbstractCost, ModelAndCost, f,
dc,calculate_cost,calculate_final_cost,
predict, simulate, df,costfun, LTVStateSpaceModel,
SimpleLTVModel, rms, sse, nrmse, modelfit, aic

import ControlSystemIdentification: OutputData, InputOutputData, InputOutputStateData, AnyInput, AbstractIdData, oftype, time1, input, output, state
import ControlSystems: ninputs, noutputs, nstates

export predict, simulate, iddata, input, output, state, ninputs, noutputs, nstates

export SimpleLTVModel, KalmanModel, LTVAutoRegressive, rootspectrogram

using ControlSystems, DSP, Plots
using Test

using DiffResults, Optim
using ReverseDiff: GradientTape, GradientConfig, gradient!

"""
Struct representing a statespace model fit using a Kalman smoother
# Fields
- `At, Bt, Pt` all of type `Array{T,3}` with time in the last dimension
- `extended::Bool` indicates whether or not the model has been extended by one time step to match the length of the input data.
"""
mutable struct KalmanModel{T} <: LTVStateSpaceModel
    At::Array{T,3}
    Bt::Array{T,3}
    Pt::Array{T,3}
    extended::Bool
    ll::Float64
end
function KalmanModel(At::Array{T,3},Bt::Array{T,3},Pt::Array{T,3},extend::Bool) where T
    if extend
        At = cat(At,At[:,:,end], dims=3)
        Bt = cat(Bt,Bt[:,:,end], dims=3)
        Pt = cat(Pt,Pt[:,:,end], dims=3)
    end
    return KalmanModel{T}(At,Bt,Pt,extend,0.0)
end

KalmanModel(At,Bt,Pt,extend::Bool=false) = KalmanModel{eltype(At)}(At,Bt,Pt,extend,0.0)

mutable struct GMMModel <: AbstractModel
    M
    dynamics
    T
end

mutable struct LTVAutoRegressive{T,PT<:Union{Array{T,3},Nothing}} <: LinearTimeVaryingModelsBase.LTVModel
    na::Int
    θ::Array{T,2}
    Pt::PT
    extended::Bool
    ll::Float64
end

Base.length(m::LTVAutoRegressive) = size(m.θ,2)

include("utilities.jl")
include("peakdetection.jl")
include("statespace_fit.jl")
include("seg_bellman.jl")
include("kalmanmodel.jl")
# include("gmmmodel.jl") # Not yet working with julia v0.6
include("statespace_utils.jl")
include("wrappers.jl")
end # module

ext: jl


text:
module BenchmarkInputDataGen

export write_csv_files

using RigidBodyDynamics
using DelimitedFiles
using ForwardDiff
using RigidBodyDynamics: OneDegreeOfFreedomFixedAxis

using InteractiveUtils: subtypes
using Base.Iterators: flatten
using StaticArrays: SVector

import Random

const JointSegmentedVector{T} = RigidBodyDynamics.CustomCollections.SegmentedVector{JointID, T, Base.OneTo{JointID}, Vector{T}}

include("sample.jl")
include("csvwriters.jl")

function create_mechanism(urdf::String; floating::Bool)
    root_joint_type = floating ? QuaternionFloating{Float64}() : Fixed{Float64}()
    mechanism = parse_urdf(urdf, root_joint_type=root_joint_type)
    mechanism.gravitational_acceleration = zero(mechanism.gravitational_acceleration)
    mechanism
end

"""
    write_csv_files(urdf; num_samples, csvdir, floating=false)

Write CSV file(s) containing randomly generated benchmark data for
a mechanism parsed from URDF `urdf`.

Keyword arguments:

* `num_samples`: number of samples to generate
* `csvdir`: directory in which to store CSV files
* `floating`: whether or not to insert a floating joint.

The output CSV files have a header with column names. The CSV contains
(respectively):

1. Joint positions ('q_JOINTNAME_SUFFIX')
2. Joint velocities ('v_JOINTNAME_SUFFIX')
3. Joint accelerations ('vd_JOINTNAME_SUFFIX')
4. Joint torques ('tau_JOINTNAME_SUFFIX')
"""
function write_csv_files(urdf; num_samples, csvdir, floating=false)
    mechanism = create_mechanism(urdf, floating=floating)

    state = MechanismState(mechanism)
    result = DynamicsResult(mechanism)
    τ = fill!(similar(velocity(state)), NaN)
    samples = Sample[]
    sizehint!(samples, num_samples)

    floating_joint_velocity_range = if floating
        floating_joint = first(joints(mechanism))
        velocity_range(state, floating_joint)
    else
        1:0
    end

    Random.seed!(1)
    for _ in 1 : num_samples
        Random.rand!(state)
        Random.rand!(τ)
        τ[floating_joint_velocity_range] .= 0
        dynamics!(result, state, τ)
        sample = Sample(copy(configuration(state)), copy(velocity(state)), copy(result.v̇), copy(τ))
        push!(samples, sample)
    end

    filename = basename(urdf)
    mkpath(csvdir)
    for Writer in subtypes(CSVWriter)
        robotname, _ = splitext(filename)
        writer = Writer(mechanism, robotname)
        libname = shortlibname(writer)
        libdir = joinpath(csvdir, libname)
        mkpath(libdir)
        write_inputs_csv(writer, libdir, robotname, samples)
        write_expected_result_csvs(writer, libdir, robotname, first(samples))
    end
end

struct URDFConfig
    basename::String
    floating::Bool
end

const URDF_DIR = joinpath(@__DIR__, "..", "..", "description", "urdf")
const URDF_CONFIGS = [URDFConfig("iiwa.urdf", false), URDFConfig("atlas.urdf", true), URDFConfig("hyq.urdf", true)]
const CSV_DIR = joinpath(@__DIR__, "..", "..", "csv")

function write_csv_files(;urdfdir=URDF_DIR, csvdir=CSV_DIR, num_samples=2)
    for config in URDF_CONFIGS
        write_csv_files(joinpath(urdfdir, config.basename); num_samples=num_samples, csvdir=csvdir, floating=config.floating)
    end
end

end # module

max_stars_repo_name: CobbledSteel/rbd-benchmarks-riscv

max_issues_repo_name: CobbledSteel/rbd-benchmarks-riscv

max_forks_repo_name: CobbledSteel/rbd-benchmarks-riscv


text:
module OutputCollectors

# In this file, we define a helper class that will run subprocesses, collecting
# the stdout and stderr, timestamped such that we can merge the two streams
# intelligently after the fact, or keep them separate for proper analysis.
import Base: wait, merge

export OutputCollector, merge, collect_stdout, collect_stderr, tail, tee

# Make sure we don't have problems when `Base.have_color` is `nothing`.
function get_have_color()
    @static if isdefined(Base, :get_have_color)
        return Base.get_have_color()
    else
        # Compatibility with Julia v1.4-.  Layman check for coloured output on
        # older Julia versions.  TODO: remove this function when we drop support
        # for oldere versions of Julia.
        return something(Base.have_color, false)
    end
end

struct LineStream
    pipe::Pipe
    lines::Vector{Tuple{Float64,String}}
    task::Task
end

"""
    readuntil_many(s::IO, delims)

Given a collection of delimiter characters, read from `s` until one of those
delimiters is reached, or we reach the end of `s`.
"""
function readuntil_many(s::IO, delims)
    out = IOBuffer()
    while !eof(s)
        c = read(s, Char)
        write(out, c)
        if c in delims
            break
        end
    end
    return String(take!(out))
end

"""
    LineStream(pipe::Pipe)

Given a `Pipe` that has been initialized by `spawn()`, create an async Task to
read in lines as they come in and annotate the time the line was captured for
later replay/merging with other simultaneously captured streams.
"""
function LineStream(pipe::Pipe, event::Condition)
    # We always have to close() the input half of the stream before we can
    # read() from it.  I don't know why, and this is honestly kind of annoying
    close(pipe.in)

    lines = Tuple{Float64,String}[]
    task = @async begin
        # Read lines in until we can't anymore.
        while true
            # Push this line onto our lines, then notify() the event
            line = readuntil_many(pipe, ['\n', '\r'])
            if isempty(line) && eof(pipe)
                break
            end
            push!(lines, (time(), line))
            notify(event)
        end
    end

    # Create a second task that runs after the first just to notify()
    # This ensures that anybody that's listening to the event but gated on our
    # being alive (e.g. `tee()`) can die alongside us gracefully as well.
    @async begin
        fetch(task)
        notify(event)
    end
    return LineStream(pipe, lines, task)
end

"""
    alive(s::LineStream)

Returns `true` if the task owned by this `LineStream` is still processing
output from an underlying `Pipe`.
"""
function alive(s::LineStream)
    return !(s.task.state in [:done, :failed])
end

"""
    OutputCollector

A `run()` wrapper class that captures subprocess `stdout` and `stderr` streams
independently, resynthesizing and colorizing the streams appropriately.
"""
mutable struct OutputCollector
    cmd::Base.AbstractCmd
    P::Base.AbstractPipe
    stdout_linestream::LineStream
    stderr_linestream::LineStream
    event::Condition
    tee_stream::IO
    verbose::Bool
    tail_error::Bool

    done::Bool
    extra_tasks::Vector{Task}

    function OutputCollector(cmd, P, out_ls, err_ls, event, tee_stream,
                             verbose, tail_error)
        return new(cmd, P, out_ls, err_ls, event, tee_stream, verbose,
                   tail_error, false, Task[])
    end
end

"""
    OutputCollector(cmd::AbstractCmd; verbose::Bool = false)

Run `cmd`, and collect the output such that `stdout` and `stderr` are captured
independently, but with the time of each line recorded such that they can be
stored/analyzed independently, but replayed synchronously.
"""
function OutputCollector(cmd::Base.AbstractCmd; verbose::Bool=false,
                         tail_error::Bool=true, tee_stream::IO=stdout)
    # First, launch the command
    out_pipe = Pipe()
    err_pipe = Pipe()
    P = try
        run(pipeline(cmd, stdin=devnull, stdout=out_pipe, stderr=err_pipe); wait=false)
    catch
        @warn("Could not spawn $(cmd)")
        rethrow()
    end

    # Next, start siphoning off the first couple lines of output and error
    event = Condition()
    out_ls = LineStream(out_pipe, event)
    err_ls = LineStream(err_pipe, event)

    # Finally, wrap this up in an object so that we can merge stdout and stderr
    # back together again at the end
    self = OutputCollector(cmd, P, out_ls, err_ls, event, tee_stream,
                           verbose, tail_error)

    # If we're set as verbose, then start reading ourselves out to stdout
    if verbose
        tee(self; stream = tee_stream)
    end

    return self
end

"""
    wait(collector::OutputCollector)

Wait for the command and all line streams within an [`OutputCollector`](@ref) to
finish their respective tasks and be ready for full merging.  Return the success
of the underlying process.  Prints out the last 10 lines of the process if it
does not complete successfully unless the OutputCollector was created as
`verbose`.
"""
function wait(collector::OutputCollector)
    # If we've already done this song and dance before, then don't do it again
    if !collector.done
        wait(collector.P)
        fetch(collector.stdout_linestream.task)
        fetch(collector.stderr_linestream.task)

        # Also fetch on any extra tasks we've jimmied onto the end of this guy
        for t in collector.extra_tasks
            fetch(t)
        end

        # From this point on, we are actually done!
        collector.done = true

        # If we failed, print out the tail of the output, unless we've been
        # tee()'ing it out this whole time, but only if the user said it's okay.
        if !success(collector.P) && !collector.verbose && collector.tail_error
            our_tail = tail(collector; colored=get_have_color())
            print(collector.tee_stream, our_tail)
        end
    end

    # Shout to the world how we've done
    return success(collector.P)
end

"""
    merge(collector::OutputCollector; colored::Bool = false)

Merge the stdout and stderr streams of the [`OutputCollector`](@ref) on a
per-line basis, returning a single string containing all collected lines,
interleaved by capture time.  If `colored` is set to true, embeds terminal color
codes to print `stderr` in red.
"""
function merge(collector::OutputCollector; colored::Bool = false)
    # First, wait for things to be done.  No incomplete mergings here yet.
    wait(collector)

    # We copy here so that you can `merge()` more than once, if you want.
    stdout_lines = copy(collector.stdout_linestream.lines)
    stderr_lines = copy(collector.stderr_linestream.lines)
    output = IOBuffer()

    # Write out an stdout line, optionally with color, and pop off that line
    function write_line(lines, should_color, color)
        if should_color && colored
            print(output, color)
        end
        t, line = popfirst!(lines)
        print(output, line)
    end

    # These help us keep track of colorizing the output
    out_color = Base.text_colors[:default]
    err_color = Base.text_colors[:red]
    last_line_stderr = false

    # Merge stdout and stderr
    while !isempty(stdout_lines) && !isempty(stderr_lines)
        # Figure out if stdout's timestamp is earlier than stderr's
        if stdout_lines[1][1] < stderr_lines[1][1]
            write_line(stdout_lines,  last_line_stderr, out_color)
            last_line_stderr = false
        else
            write_line(stderr_lines, !last_line_stderr, err_color)
            last_line_stderr = true
        end
    end

    # Now drain whichever one still has data within it
    while !isempty(stdout_lines)
        write_line(stdout_lines, last_line_stderr, out_color)
        last_line_stderr = false
    end
    while !isempty(stderr_lines)
        write_line(stderr_lines, !last_line_stderr, err_color)
        last_line_stderr = true
    end

    # Clear text colors at the end, if we need to
    if last_line_stderr && colored
        print(output, Base.text_colors[:default])
    end

    # Return our ill-gotten goods
    return String(take!(output))
end

"""
    collect_stdout(collector::OutputCollector)

Returns all stdout lines collected by this collector so far.
"""
function collect_stdout(collector::OutputCollector)
    return join([l[2] for l in collector.stdout_linestream.lines], "")
end

"""
    collect_stderr(collector::OutputCollector)

Returns all stderr lines collected by this collector so far.
"""
function collect_stderr(collector::OutputCollector)
    return join([l[2] for l in collector.stderr_linestream.lines], "")
end

"""
    tail(collector::OutputCollector; len::Int = 100, colored::Bool = false)

Write out the last `len` lines, optionally writing colored lines.
"""
function tail(collector::OutputCollector; len::Int = 100,
              colored::Bool = false)
    out = merge(collector; colored=colored)

    idx = length(out)
    for line_idx in 1:len
        # We can run into UnicodeError's here
        try
            idx = findprev(isequal('\n'), out, idx-1)
            # We have to check for both `nothing` or `0` for Julia 0.6
            if idx === nothing || idx == 0
                idx = 0
                break
            end
        catch
            break
        end
    end

    return out[idx+1:end]
end

"""
    tee(c::OutputCollector; colored::Bool = false, stream::IO = stdout)

Spawn a background task to incrementally output lines from `collector` to the
standard output, optionally colored.
"""
function tee(c::OutputCollector; colored::Bool=get_have_color(),
             stream::IO=stdout)
    tee_task = @async begin
        out_idx = 1
        err_idx = 1
        out_lines = c.stdout_linestream.lines
        err_lines = c.stderr_linestream.lines

        # Helper function to print out the next line of stdout/stderr
        function print_next_line()
            timestr = Libc.strftime("[%T] ", time())
            # We know we have data, so figure out if it's for stdout, stderr
            # or both, and we need to choose which to print based on timestamp
            printstyled(stream, timestr; bold=true)
            if length(out_lines) >= out_idx
                if length(err_lines) >= err_idx
                    # If we've got input waiting from both lines, then output
                    # the one with the lowest capture time
                    if out_lines[out_idx][1] < err_lines[err_idx][1]
                        # Print the out line as it's older
                        print(stream, out_lines[out_idx][2])
                        out_idx += 1
                    else
                        # Print the err line as it's older
                        printstyled(stream, err_lines[err_idx][2]; color=:red)
                        print(stream)
                        err_idx += 1
                    end
                else
                    # Print the out line that is the only one waiting
                    print(stream, out_lines[out_idx][2])
                    out_idx += 1
                end
            else
                # Print the err line that is the only one waiting
                printstyled(stream, err_lines[err_idx][2]; color=:red)
                print(stream)
                err_idx += 1
            end
        end

        # First thing, wait for some input.  This avoids us trying to inspect
        # the liveliness of the linestreams before they've even started.
        wait(c.event)

        while alive(c.stdout_linestream) || alive(c.stderr_linestream)
            if length(out_lines) >= out_idx || length(err_lines) >= err_idx
                # If we have data to output, then do so
                print_next_line()
            else
                # Otherwise, wait for more input
                wait(c.event)
            end
        end

        # Drain the rest of stdout and stderr
        while length(out_lines) >= out_idx || length(err_lines) >= err_idx
            print_next_line()
        end
    end

    # Let the collector know that he might have to wait on this `tee()` to
    # finish its business as well.
    push!(c.extra_tasks, tee_task)

    return tee_task
end

end # module OutputCollectors


text:
using Raven
using Test

@testset "Raven" begin
 @test 2 + 2 == 5
 Raven.readEdges(open("testdata/scl.edge","r"))

end


text:
using NiLang, NiLang.AD
using LinearAlgebra

function _iqr(A)
    @assert size(A, 1) == size(A, 2)
    N = size(A, 1)
    Q = zero(A)
    R = zero(A)
    for col = 1:N
        ri = A[:,col]
        for precol = 1:col-1
            ni = Q[:,precol]'*ri
            R[precol,col] = ni
            ri = ri - ni*Q[:,precol]
        end
        R[col, col] = norm(ri)
        Q[:,col] = ri/R[col,col]
    end
    return Q, R
end

@i function iqr(Q, R, A::AbstractMatrix{T}) where T
    @anc anc_norm = zero(T)
    @anc anc_dot = zeros(T, size(A,2))
    @anc ri = zeros(T, size(A,1))
    for col = 1:size(A, 1)
        ri .+= identity.(A[:,col])
        for precol = 1:col-1
            idot(anc_dot[precol], Q[:,precol], ri)
            R[precol,col] += identity(anc_dot[precol])
            for row = 1:size(Q,1)
                ri[row] -= anc_dot[precol] * Q[row, precol]
            end
        end
        inorm2(anc_norm, ri)

        R[col, col] += anc_norm^0.5
        for row = 1:size(Q,1)
            Q[row,col] += ri[row] / R[col, col]
        end

        ~begin
            ri .+= identity.(A[:,col])
            for precol = 1:col-1
                idot(anc_dot[precol], Q[:,precol], ri)
                for row = 1:size(Q,1)
                    ri[row] -= anc_dot[precol] * Q[row, precol]
                end
            end
            inorm2(anc_norm, ri)
        end
    end
end

using Test, Random
@testset "test qr" begin
    Random.seed!(2)
    A = randn(4,4)
    Q, R = _iqr(A)
    @test Q * R ≈ A
    @test Q'*Q ≈ I
    R[abs.(R).<1e-10] .= 0.0
    @test istriu(R)

    q = zero(A)
    r = zero(A)
    @instr iqr(q, r, A)
    @test q ≈ Q
    @test r ≈ R
    @test check_inv(iqr, (q, r, A))

    @i function test1(out, q, r, A)
        iqr(q, r, A)
        out += identity(q[1,2])
    end
    @i function test2(out, q, r, A)
        iqr(q, r, A)
        out += identity(r[1,2])
    end
    @test check_grad(test1, (Loss(0.0), q, r, A); atol=0.05, verbose=true)
    @test check_grad(test2, (Loss(0.0), q, r, A); atol=0.05, verbose=true)
end


text:
mutable struct S a::Int64  end
function f(x, y::Int64) x.a = y end
f("foo", 1)

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
@testset "Tracking" begin
    gpsl1 = Tracking.GPSL1()
    carrier = Tracking.gen_carrier.(1:24000, 50Hz, 1.2, 4e6Hz)
    code = Tracking.gen_code.(Ref(gpsl1), 1:24000, 1023e3Hz, 2.0, 4e6Hz, 1)
    signal = carrier .* code * [1 1 1 1]
    correlator_outputs = zeros(SMatrix{3,4,ComplexF64})
    code_shift = Tracking.CodeShift{3}(gpsl1, 4e6Hz, 0.5)
    inits = TrackingInitials(20Hz, 1.2, 0.0Hz, 2.0)
    dopplers = Tracking.Dopplers(inits)
    phases = Tracking.Phases(inits)
    carrier_loop = Tracking.init_3rd_order_bilinear_loop_filter(18Hz)
    code_loop = Tracking.init_2nd_order_bilinear_loop_filter(1Hz)
    last_valid_correlator_outputs = copy(correlator_outputs)
    data_bits = Tracking.DataBits(gpsl1)
    cn0_state = Tracking.CN0State(20ms)
    results = @inferred Tracking._tracking(correlator_outputs, last_valid_correlator_outputs, signal, gpsl1, 4e6Hz, 30Hz, inits, dopplers, phases, code_shift, carrier_loop, code_loop, 1, x -> x[:,1], 0.5ms, 1ms, 1, 0, 0, data_bits, cn0_state, 0.0Hz)
    @test results[2].carrier_doppler ≈ 20Hz
    @test results[2].code_doppler ≈ 0Hz atol = 3e-3Hz #??
    @test results[2].code_phase ≈ 2 atol = 2e-5
end

@testset "Track L1" begin

     center_freq = 1.57542e9Hz
     code_freq = 1023e3Hz
     doppler = 10Hz
     interm_freq = 50Hz
     carrier_phase = π / 3
     sample_freq = 4e6Hz
     code_doppler = doppler * code_freq / center_freq
     data_doppler = doppler * 50Hz / center_freq
     code_phase = 2.0
     min_integration_time = 0.5ms

     run_time = 1500e-3s
     integration_time = 1e-3s
     num_integrations = convert(Int, run_time / integration_time)
     num_samples = convert(Int, run_time * sample_freq)
     integration_samples = convert(Int, integration_time * sample_freq)

     gps_l1 = Tracking.GPSL1()

     carrier = cis.(2π * (interm_freq + doppler) / sample_freq * (1:num_samples) .+ carrier_phase)
     sampled_code = Tracking.gen_code.(Ref(gps_l1), 1:num_samples, code_doppler + code_freq, code_phase, sample_freq, 1)
     signal = carrier .* sampled_code * [1 1 1 1]

     inits = TrackingInitials(0.0Hz, carrier_phase, 0.0Hz, code_phase)
     track = @inferred PhasedArrayTracking.init_tracking(Val(4), gps_l1, inits, sample_freq, interm_freq, 1, pll_bandwidth = 18.0Hz, dll_bandwidth = 1.0Hz, min_integration_time = min_integration_time, max_integration_time = integration_time)

     code_dopplers = zeros(num_integrations)
     code_phases = zeros(num_integrations)
     calculated_code_phases  = mod.((1:num_integrations) * integration_samples * (code_doppler + code_freq) / sample_freq .+ code_phase, 1023)
     carrier_dopplers = zeros(num_integrations)

     results = nothing
     for i = 1:num_integrations
         current_signal = signal[integration_samples * (i - 1) + 1:integration_samples * i, :]# .+ complex.(randn(integration_samples,2), randn(integration_samples,2)) .* 10^(5/20)
         track, results = @inferred track(current_signal, x -> x[:,1])
         @test results.num_processed_prns == 1
         @test results.num_bits == 0
         @test results.data_bits == 0
         code_phases[i] = results.code_phase
         carrier_dopplers[i] = results.carrier_doppler / Hz
         code_dopplers[i] = results.code_doppler / Hz
     end

     @test results.carrier_doppler ≈ doppler atol = 4e-4Hz
     @test mod(results.code_phase, 1023) ≈ calculated_code_phases[end] atol = 3e-6
     @test results.code_doppler ≈ code_doppler atol = 4e-6Hz
end


text:
module CircuitList

using UUIDs
import Base.append!
import Base.iterate
import Base.push!

"""
Node element in LList-type structure. Maintains forward and reverse links, data, and UUID.
"""
mutable struct Node{T}
    data::Union{Nothing, T}
    next::Union{Nothing, Node{T}}
    prev::Union{Nothing, Node{T}}
    uuid::UUID
    Node{T}() where {T} = new(nothing, nothing, nothing, UUIDs.uuid1())
    Node{T}(data) where {T} = new(data, nothing, nothing, UUIDs.uuid1())
    Node{T}(data, prev) where {T} = new(data, nothing, prev, UUIDs.uuid1())
    Node{T}(data, next, prev) where {T} = new(data, next, prev, UUIDs.uuid1())
end

"""
DLList structure for representing GateCalls. Maintains start, end and position pointer.
"""
mutable struct CList{T}
    len::Int
    head::Union{Nothing, Node{T}}
    tail::Union{Nothing, Node{T}}
    current::Union{Nothing, Node{T}}
    CList{T}() where {T} = new(0, nothing, nothing, nothing)
end

function push!(cl::CList{T}, data::T) where T
    if cl.len == 0
        n = Node{T}(data)
        cl.head = n
        cl.tail = n
        cl.current = n
        cl.current = cl.head
    else
        cl.tail.next = Node{T}(data, cl.tail)
        cl.tail.next.prev = cl.tail
        cl.tail = cl.tail.next
    end
    cl.len += 1
end

function clear!(l::CList{T}) where T
    l.head = nothing
    l.tail = nothing
    l.current = nothing
    l.len = 0
end

function append!(l1::CList{T}, l2::CList{T}) where T
    if l2.current == nothing
        return
    end
    l2c = deepcopy(l2)
    l1.tail.next = l2c.head
    l2c.head.prev = l1.tail
    l1.len += l2c.len
    clear!(l2)
    return l1
end

function copy(l::CList{T}) where T
    lc = CList{T}()
    for val in l
        push!(lc, val)
    end
    return lc
end

function iterate(list::CircuitList.CList{T}) where T
    if list.len == 0
        return nothing
    else
        list.current = list.head
        return (list.current.data, list.current.next)
    end
end

function iterate(list::CircuitList.CList{T}, n::Union{CircuitList.Node{T}, Nothing}) where T
    if list.current.next == nothing
        return nothing
    else
        res = (list.current.next.data, list.current.next)
        list.current = list.current.next
        return res
    end
end

end


text:
module SBTestDelaunay

using SpatialBoundaries
using Test

n = 21
x = rand(n)
y = rand(n)
z = rand(n)

W = wombling(x, y, z)
@test isa(W, TriangulationWomble)
@test length(W.m) == length(W.x)
@test length(W.x) == length(W.y)

end


text:
# using a closure to keep count local
function func_gen()
    count = 0
    return function(x::Vector,g::Vector)
        if length(g) > 0
            g[1] = 0 # objective function is just the value of x2
            g[2] = 1 # objective function is just the value of x2
        end
        count::Int += 1
        println("f_$count($x)")
        return x[2]
    end
end

#inequality constraint 1
function ic1(x::Vector, g::Vector)
    if length(g) > 0
        g[1] = 0
        g[2] = -1
    end
    return 1-x[2]
end

#inequality constraint 2
function ic2(x::Vector, g::Vector)
    if length(g) > 0
        g[1] = -1
        g[2] = 0
    end
    return -x[1]
end

#equality constraint 1 (empty)
function ec1(x::Vector, g::Vector)
    if length(g) > 0
        g[1] = -2*x[1]
        g[2] = 1
    end
    return x[2]-x[1]^2
end

function quadratic_ex()
    opt = NLopt.Opt(:LD_SLSQP, 2)
    NLopt.lower_bounds!(opt, [-Inf, 0.])
    NLopt.xtol_rel!(opt, 1e-4)

    NLopt.min_objective!(opt, func_gen())
    NLopt.inequality_constraint!(opt, ic1)
    NLopt.inequality_constraint!(opt, ic2)
    NLopt.equality_constraint!(opt, ec1, 0)

    (minf, minx, ret) = NLopt.optimize(opt, [5, 20])
    println("got $minf at $minx.")
end


text:
include("Graph.jl")
include("SimMetric.jl")
include("Exp.jl")

using Printf

function runMain(ags1; resultDir = "../result/")
    G = readLabeledGraph(ags1)
    println(ags1, " ", G.n, " ", G.m)
    # RoleSim
    out1 = open(resultDir*ags1*"RoleSim.txt", "w")
    @time result = topK(G, RoleSim(G))
    for i = 1 : size(result, 1)
        println(out1, i, " ", result[i])
    end
    close(out1)
    # ForestSim-Exact
    out2 = open(resultDir*ags1*"ForestSimExact.txt", "w")
    @time result = topK(G, ForestSim(G))
    for i = 1 : size(result, 1)
        println(out2, i, " ", result[i])
    end
    close(out2)
    # ForestSim-Approx
    out3 = open(resultDir*ags1*"ForestSimApprox.txt", "w")
    @time result = topKofForestSimApprox(G)
    for i = 1 : size(result, 1)
        println(out3, i, " ", result[i])
    end
    close(out3)
    # StructSim
    out4 = open(resultDir*ags1*"StructSim.txt", "w")
    @time result = topK(G, StructSim(G))
    for i = 1 : size(result, 1)
        println(out4, i, " ", result[i])
    end
    close(out4)
end

runMain("karate")
runMain("brazil-airports")
runMain("europe-airports")
runMain("expressway-starbucks")
runMain("expressway-status")
runMain("actor")


text:
# This file has simple functions that mimic (part of) the behavior of the Matlab functions with the same name
# Figure properties like for example WindowButtonDownFcn() were implemented as functions too.

include("libiup_h.jl")
include("libiup.jl")

type Handles
	figure1::Ptr{Ihandle}
	iup_canvas::Ptr{Ihandle}
	cd_canvas::Ptr{cdCanvas}	# cdCanvas is a composite type
end

### Garbage collection [prevention]
const gc_preserve = ObjectIdDict() # reference counted closures
function gc_ref(x::ANY)
	global gc_preserve
	#isbits(x) && error("can't gc-preserve an isbits object")
	gc_preserve[x] = (get(gc_preserve, x, 0)::Int)+1
	x
end
function gc_unref(x::ANY)
	global gc_preserve
	#@assert !isbits(x)
	count = get(gc_preserve, x, 0)::Int-1
	if count <= 0
		delete!(gc_preserve, x)
	end
	nothing
end

function guidata(hand::Ptr{Ihandle})
	# Get the handles stored in the Ihandle hand
	handles = unsafe_pointer_to_objref(IupGetAttribute(hand, "handles"))
	return handles
end

# -------------------------------------------------------------------------------
function guidata(hand::Ptr{Ihandle}, handles::Handles)
	# Store the handles in the Ihandle hand
	IupSetAttribute(hand, "handles", pointer_from_objref(handles))
end
#guidata(hand::Ptr{Void}, handles::Handles) = guidata(convert(Ptr{Ihandle}, hand), handles)

# -------------------------------------------------------------------------------
function setappdata(hand::Ptr{Ihandle}, name::String, val)
	# Store the pointer of the whatever in 'val' in a container names 'name'
	IupSetAttribute(hand, name, pointer_from_objref(val))
end

# -------------------------------------------------------------------------------
function getappdata(hand::Ptr{Ihandle}, name::String)
	# Do the oposite of setappdata
	att = IupGetAttribute(hand, name)
	if (att != C_NULL)
		val = unsafe_pointer_to_objref(att)
	else
		val = []
	end
	return val
end

# -------------------------------------------------------------------------------
function WindowButtonDownFcn(canvas::Ptr{Ihandle}, f::Function)
	# Sets the Button press callback function which is executed whenever a mouse button is pressed
	# The interface for the callback function is described in IUP docs under BUTTON_CB
	# int function(Ihandle* ih, int button, int pressed, int x, int y, char* status); [in C]
	# Warning: the macros that make use of the status variable are not yet implemented
	IupSetCallback(canvas, "BUTTON_CB", cfunction(f, Int, (Ptr{Ihandle}, Char, Cint, Cint, Cint, Ptr{UInt8})))
end

# -------------------------------------------------------------------------------
function WindowButtonMotionFcn(canvas::Ptr{Ihandle}, f::Function)
	# Sets the Mouse motion callback function
	# The interface for the callback function is described in IUP docs under MOTION_CB
	# int function(Ihandle *ih, int x, int y, char *status); [in C]
	# Warning: the macros that make use of the status variable are not yet implemented
	IupSetCallback(canvas, "MOTION_CB", cfunction(f, Int, (Ptr{Ihandle}, Cint, Cint, Ptr{UInt8})))
end

# -------------------------------------------------------------------------------
function KeyPressFcn(canvas::Ptr{Ihandle}, f::Function)
	# Sets the Key press callback function, which is invoked when a key is pressed or released.
	# The interface for the callback function is described in IUP docs under KEYPRESS_CB
	# int function(Ihandle *ih, int c, int press); [in C]
	IupSetCallback(canvas, "KEYPRESS_CB", cfunction(f, Int, (Ptr{Ihandle}, Char, Cint)))
end

# -------------------------------------------------------------------------------
function set_current_point(canvas::Ptr{Ihandle}, pt)
	# Stores the coordinates of the last clicked point in the canva's "CurrentPoint" attribute.
	# ATTENTION: this function MUST be called from inside the WindowButtonDownFcn() callback function
	setappdata(canvas, "CurrentPoint", pt)
end

# -------------------------------------------------------------------------------
function get_current_point(canvas::Ptr{Ihandle})
	# Get the location of last button clicked.
	# In Matlab the corresponding property depends if it refers to a Figure, case in which
	# it returns a two elements vector with [x-coordinate, y-coordinate], or an Axes object,
	# the output is a 2x3 matrix with [xfront yfront zfront; xback yback zback].
	# For the time being here we are only returning teh 1x2 Array version ... or empty is prop not set
	# MANDATORY: for this to work it is crutial that the WindowButtonDownFcn() callback function
	# is set and calls the set_current_point() function (which sets the appropriate attribute)
	cpt = getappdata(canvas, "CurrentPoint")
	return cpt
end

# -------------------------------------------------------------------------------
function set_current_character(canvas::Ptr{Ihandle}, s::Char)
	# PROBABLY, NOT FINISHED
	setappdata(canvas, "CurrentCharacter", s)
end

# -------------------------------------------------------------------------------
function get_current_character(canvas::Ptr{Ihandle})
	# Get the last key pressed. We store this info as an attribute in the canvas handle.
	# MANDATORY: for this to work it is crutial that the WindowButtonDownFcn() callback function
	# is set and calls the set_current_character() function (which sets the appropriate attribute)
	cc = getappdata(canvas, "CurrentCharacter")
	return cc
end

# -------------------------------------------------------------------------------
function isnada(x::Any)
	# Replacement for isempty that works for any type
	empty = false
	try
		isempty(x)
		empty = true
	catch
		empty = false
	end
	return empty
end


text:
"""
MOS interface struct.
"""
mutable struct Interface
    url::String
    token::String

    @doc """
        Interface(url::String, token::String)

    Create an interface to MOS at the given url and using
    the given token for authentication.
    """
    Interface(url::String, token::String) = new(url, token)

    @doc """
        Interface(url::String)

    Create an interface to MOS at the given url.
    The authentication token is assumed to be available in the
    environment variable `MOS_BACKEND_TOKEN`.
    """
    function Interface(url::String)
        if haskey(ENV, "MOS_BACKEND_TOKEN")
            token = ENV["MOS_BACKEND_TOKEN"]
        else
            token = ""
        end
        return Interface(url, token)
    end

    @doc """
        Interface()

    Create an interface to MOS
    The url is constructed from the environment variables
    `MOS_BACKEND_HOST` and `MOS_BACKEND_HOST`. The authentication token is
    assumed to be available in the environment variable `MOS_BACKEND_TOKEN`.
    """
    function Interface()
        if haskey(ENV, "MOS_BACKEND_HOST") && haskey(ENV, "MOS_BACKEND_PORT")
            host = ENV["MOS_BACKEND_HOST"]
            port = ENV["MOS_BACKEND_PORT"]
            if port == "443"
              protocol = "https"
            else
              protocol = "http"
            end
            url = "$protocol://$host:$port/api/"
        else
            url = "https://mos.fuinn.ie:443/api/"
        end
        return Interface(url)
    end

end

"""
    new_model(i::Interface, filepath::AbstractString)::Model

Create new MOS model from local annotated file.
"""
function new_model(i::Interface, filepath::AbstractString)::Model
    url = join([i.url, "model/create_from_file/"])
    open(filepath, "r") do source_file
        f = HTTP.Form(Dict("source_file" => source_file))
        h = Dict("Content-Type" => "multipart/form-data; boundary=$(f.boundary)")
        r = HTTP.post(url, add_auth!(h, i.token), body=f)
        if r.status != 200
            error("unable to create model")
        end
        return Model(i, JSON.parse(String(r.body)))
    end
end

"""
    delete_model_with_name(i::Interface, name::String)

Delete model with given name.
"""
function delete_model_with_name(i::Interface, name::String)
    url = join([i.url, "model/"])
    h = Dict()
    r = HTTP.get(url, add_auth!(h, i.token), query=["name" => name])
    if r.status != 200
        error("unable to get model")
    end
    models = JSON.parse(String(r.body))
    for m in models
        model = Model(i, m)
        id = get_id(model)
        url = join([i.url, "model/$id/"])
        h = Dict()
        r = HTTP.delete(url, add_auth!(h, i.token))
        if r.status != 204
            error("unable to delete model")
        end
    end
end

"""
    get_model(i::Interface, id::Integer)::Model

Get model with given id.
"""
function get_model(i::Interface, id::Integer)::Model
    url = join([i.url, "model/$id/"])
    h = Dict()
    r = HTTP.get(url, add_auth!(h, i.token))
    if r.status != 200
        error("unable to get model")
    end
    return Model(i, JSON.parse(String(r.body)))
end

"""
    get_model_with_name(i::Interface, name::String)::Model

Get model with given name.
"""
function get_model_with_name(i::Interface, name::String)::Model
    url = join([i.url, "model/"])
    h = Dict()
    r = HTTP.get(url, add_auth!(h, i.token), query=["name" => name])
    if r.status != 200
        error("unable to get model")
    end
    models = JSON.parse(String(r.body))
    if isempty(models)
        error("No model found with name $name")
    elseif length(models) > 1
        error("More than one model found with name $name")
    else
        m = Model(i, models[1])
        id = m.data["id"]
        return get_model(i, id)
    end
end

max_stars_repo_name: Fuinn/mos-interface-jl

max_issues_repo_name: Fuinn/mos-interface-jl

max_forks_repo_name: Fuinn/mos-interface-jl


text:
module InternalClimaParams

import CLIMAParameters
const CP = CLIMAParameters
const CPP = CP.Planet
const APS = CP.AbstractParameterSet

T_min(param_set::APS) = CPP.T_min(param_set)
T_max(param_set::APS) = CPP.T_max(param_set)
MSLP(param_set::APS) = CPP.MSLP(param_set)
R_d(param_set::APS) = CPP.R_d(param_set)
cp_d(param_set::APS) = CPP.cp_d(param_set)
kappa_d(param_set::APS) = CPP.kappa_d(param_set)
molmass_ratio(param_set::APS) = CPP.molmass_ratio(param_set)
cp_v(param_set::APS) = CPP.cp_v(param_set)
cp_l(param_set::APS) = CPP.cp_l(param_set)
cp_i(param_set::APS) = CPP.cp_i(param_set)
cv_d(param_set::APS) = CPP.cv_d(param_set)
cv_v(param_set::APS) = CPP.cv_v(param_set)
cv_l(param_set::APS) = CPP.cv_l(param_set)
cv_i(param_set::APS) = CPP.cv_i(param_set)
T_0(param_set::APS) = CPP.T_0(param_set)
e_int_v0(param_set::APS) = CPP.e_int_v0(param_set)
e_int_i0(param_set::APS) = CPP.e_int_i0(param_set)
LH_v0(param_set::APS) = CPP.LH_v0(param_set)
LH_s0(param_set::APS) = CPP.LH_s0(param_set)
LH_f0(param_set::APS) = CPP.LH_f0(param_set)
press_triple(param_set::APS) = CPP.press_triple(param_set)
R_v(param_set::APS) = CPP.R_v(param_set)
T_triple(param_set::APS) = CPP.T_triple(param_set)
T_freeze(param_set::APS) = CPP.T_freeze(param_set)
T_icenuc(param_set::APS) = CPP.T_icenuc(param_set)
pow_icenuc(param_set::APS) = CPP.pow_icenuc(param_set)
entropy_reference_temperature(param_set::APS) =
    CPP.entropy_reference_temperature(param_set)
entropy_dry_air(param_set::APS) = CPP.entropy_dry_air(param_set)
entropy_water_vapor(param_set::APS) = CPP.entropy_water_vapor(param_set)
T_surf_ref(param_set::APS) = CPP.T_surf_ref(param_set)
T_min_ref(param_set::APS) = CPP.T_min_ref(param_set)
grav(param_set::APS) = CPP.grav(param_set)

end


text:
using ScanImageTiffReader
using Test
using DotEnv
using JSON

DotEnv.config()

@testset "failures" begin
    filename = "/not/a/file"
    @test_throws Exception ScanImageTiffReader.open(filename, size)
end

if haskey(ENV, "TESTBASE")
    testbase = ENV["TESTBASE"]

    include("test_9ketamineoriginalcropped.jl")
    include("test_BigTIFF.jl")
    include("test_BigTIFFLong.jl")
    include("test_BigTIFFLong8.jl")
    include("test_BigTIFFMotorola.jl")
    include("test_BigTIFFMotorolaLongStrips.jl")
    include("test_BigTIFFSubIFD4.jl")
    include("test_BigTIFFSubIFD8.jl")
    include("test_Classic.jl")
    include("test_TR_003.jl")
    include("test_lin_00001.jl")
    include("test_linfree_00001.jl")
    include("test_linfreej_00001.jl")
    include("test_linj_00001.jl")
    include("test_oldfmt.jl")
    include("test_oldfmtj.jl")
    include("test_res_00001.jl")
    include("test_resfree_00001.jl")
    include("test_resfreej_00001.jl")
    include("test_resj_00001.jl")
    include("test_resj_2018a_00002.jl")
    include("test_single_image_si.jl")
else
    testbase = mktempdir()

    # download a test file
    dsource = "https://drive.google.com/uc?export=download&id=1Zbje2OG1QqUr9D4UGd_ag2iV6aiM_6j0"

    import HTTP
    r = HTTP.request("GET", dsource; redirect=true)
    open(joinpath(testbase, "linj_00001.tif"), "w") do io
        write(io, r.body)
    end

    include("test_linj_00001.jl")

    # cleanup
    rm(joinpath(testbase, "linj_00001.tif"))
    rm(testbase)
end


text:
"""
decimalise \\

For taking multivariate random variables and returning them as a single variable \\
in the base of the maximum value. \\
Arguments:  x, the random variable \\
            xn, number of variables in x \\
            xmax, the maximum number of possible values across the dimensions in x \\
Returns:    dec_x, Array{Int}, decimalised x \\
"""

function decimalise(x::Array{Int}, xn::Int, xmax::Int)
  if size(x)[1] != xn || maximum(x) > xmax
    error("Input dimension error!")
  end
  power_vector = xmax .^ (xn-1:-1:0)
  dec_x = size(x)[1] == 1 ? reshape(x, length(x)) : transpose(x) * power_vector
  return dec_x
end


text:
# Disjoint sets

############################################################
#
#   A forest of disjoint sets of integers
#
#   Since each element is an integer, we can use arrays
#   instead of dictionary (for efficiency)
#
#   Disjoint sets over other key types can be implemented
#   based on an IntDisjointSets through a map from the key
#   to an integer index
#
############################################################

_intdisjointsets_bounds_err_msg(T) = "the maximum number of elements in IntDisjointSets{$T} is $(typemax(T))"

"""
    IntDisjointSets{T<:Integer}(n::Integer)

A forest of disjoint sets of integers, which is a data structure
(also called a union–find data structure or merge–find set)
that tracks a set of elements partitioned
into a number of disjoint (non-overlapping) subsets.
"""
mutable struct IntDisjointSets{T<:Integer}
    parents::Vector{T}
    ranks::Vector{T}
    ngroups::T

    # creates a disjoint set comprised of n singletons
    # IntDisjointSets(n::T) where {T<:Integer} = new(collect(Base.OneTo(n)), zeros(T, n), n)
end

IntDisjointSets(n::T) where {T<:Integer} = IntDisjointSets{T}(collect(Base.OneTo(n)), zeros(T, n), n)
IntDisjointSets{T}(n::Integer) where {T<:Integer} = IntDisjointSets{T}(collect(Base.OneTo(T(n))), zeros(T, T(n)), T(n))
Base.length(s::IntDisjointSets) = length(s.parents)

"""
    num_groups(s::IntDisjointSets)

Get a number of groups.
"""
num_groups(s::IntDisjointSets) = s.ngroups
Base.eltype(::Type{IntDisjointSets{T}}) where {T<:Integer} = T

# find the root element of the subset that contains x
# path compression is implemented here
function find_root_impl!(parents::Vector{T}, x::Integer) where {T<:Integer}
    p = parents[x]
    @inbounds if parents[p] != p
        parents[x] = p = _find_root_impl!(parents, p)
    end
    return p
end

# unsafe version of the above
function _find_root_impl!(parents::Vector{T}, x::Integer) where {T<:Integer}
    @inbounds p = parents[x]
    @inbounds if parents[p] != p
        parents[x] = p = _find_root_impl!(parents, p)
    end
    return p
end

"""
    find_root!(s::IntDisjointSets{T}, x::T)

Find the root element of the subset that contains an member x.
Path compression is implemented here.
"""
find_root!(s::IntDisjointSets{T}, x::T) where {T<:Integer} = find_root_impl!(s.parents, x)

"""
    in_same_set(s::IntDisjointSets{T}, x::T, y::T)

Returns `true` if `x` and `y` belong to the same subset in `s` and `false` otherwise.
"""
in_same_set(s::IntDisjointSets{T}, x::T, y::T) where {T<:Integer} = find_root!(s, x) == find_root!(s, y)

"""
    union!(s::IntDisjointSets{T}, x::T, y::T)

Merge the subset containing x and that containing y into one
and return the root of the new set.
"""
function Base.union!(s::IntDisjointSets{T}, x::T, y::T) where {T<:Integer}
    parents = s.parents
    xroot = find_root_impl!(parents, x)
    yroot = find_root_impl!(parents, y)
    xroot != yroot ? root_union!(s, xroot, yroot) : xroot
end

"""
    root_union!(s::IntDisjointSets{T}, x::T, y::T)

Form a new set that is the union of the two sets whose root elements are
x and y and return the root of the new set.
Assume x ≠ y (unsafe).
"""
function root_union!(s::IntDisjointSets{T}, x::T, y::T) where {T<:Integer}
    parents = s.parents
    rks = s.ranks
    @inbounds xrank = rks[x]
    @inbounds yrank = rks[y]

    if xrank < yrank
        x, y = y, x
    elseif xrank == yrank
        rks[x] += one(T)
    end
    @inbounds parents[y] = x
    s.ngroups -= one(T)
    return x
end

"""
    push!(s::IntDisjointSets{T})

Make a new subset with an automatically chosen new element x.
Returns the new element. Throw an `ArgumentError` if the
capacity of the set would be exceeded.
"""
function Base.push!(s::IntDisjointSets{T}) where {T<:Integer}
    l = length(s)
    l < typemax(T) || throw(ArgumentError(_intdisjointsets_bounds_err_msg(T)))
    x = l + one(T)
    push!(s.parents, x)
    push!(s.ranks, zero(T))
    s.ngroups += one(T)
    return x
end

"""
    DisjointSets{T}(xs)

A forest of disjoint sets of arbitrary value type T.

It is a wrapper of IntDisjointSets{Int}, which uses a
dictionary to map the input value to an internal index.
"""
mutable struct DisjointSets{T} <: AbstractSet{T}
    intmap::Dict{T,Int}
    revmap::Vector{T}
    internal::IntDisjointSets{Int}

    DisjointSets{T}() where T = new{T}(Dict{T,Int}(), Vector{T}(), IntDisjointSets(0))
    function DisjointSets{T}(xs) where T    # xs must be iterable
        imap = Dict{T,Int}()
        rmap = Vector{T}()
        n = length(xs)
        sizehint!(imap, n)
        sizehint!(rmap, n)
        id = 0
        for x in xs
            imap[x] = (id += 1)
            push!(rmap,x)
        end
        new{T}(imap, rmap, IntDisjointSets(n))
    end
end

DisjointSets() = DisjointSets{Any}()
DisjointSets(xs::T...) where T = DisjointSets{T}(xs)
DisjointSets{T}(xs::T...) where T = DisjointSets{T}(xs)
DisjointSets(xs) = _DisjointSets(xs, Base.IteratorEltype(xs))
_DisjointSets(xs, ::Base.HasEltype) = DisjointSets{eltype(xs)}(xs)
function _DisjointSets(xs, ::Base.EltypeUnknown)
    T = Base.@default_eltype(xs)
    (isconcretetype(T) || T === Union{}) || return Base.grow_to!(DisjointSets{T}(), xs)
    return DisjointSets{T}(xs)
end

Base.iterate(s::DisjointSets) = iterate(s.revmap)
Base.iterate(s::DisjointSets, i) = iterate(s.revmap, i)

Base.length(s::DisjointSets) = length(s.internal)

"""
    num_groups(s::DisjointSets)

Get a number of groups.
"""
num_groups(s::DisjointSets) = num_groups(s.internal)
Base.eltype(::Type{DisjointSets{T}}) where T = T
Base.empty(s::DisjointSets{T}, ::Type{U}=T) where {T,U} = DisjointSets{U}()
function Base.sizehint!(s::DisjointSets, n::Integer)
    sizehint!(s.intmap, n)
    sizehint!(s.revmap, n)
    return s
end

"""
    find_root!{T}(s::DisjointSets{T}, x::T)

Finds the root element of the subset in `s` which has the element `x` as a member.
"""
find_root!(s::DisjointSets{T}, x::T) where {T} = s.revmap[find_root!(s.internal, s.intmap[x])]

"""
    in_same_set(s::DisjointSets{T}, x::T, y::T)

Returns `true` if `x` and `y` belong to the same subset in `s` and `false` otherwise.
"""
in_same_set(s::DisjointSets{T}, x::T, y::T) where {T} = in_same_set(s.internal, s.intmap[x], s.intmap[y])

"""
    union!(s::DisjointSets{T}, x::T, y::T)

Merge the subset containing x and that containing y into one
and return the root of the new set.
"""
Base.union!(s::DisjointSets{T}, x::T, y::T) where {T} = s.revmap[union!(s.internal, s.intmap[x], s.intmap[y])]

"""
    root_union!(s::DisjointSets{T}, x::T, y::T)

Form a new set that is the union of the two sets whose root elements are
x and y and return the root of the new set.
Assume x ≠ y (unsafe).
"""
root_union!(s::DisjointSets{T}, x::T, y::T) where {T} = s.revmap[root_union!(s.internal, s.intmap[x], s.intmap[y])]

"""
    push!(s::DisjointSets{T}, x::T)

Make a new subset with an automatically chosen new element x.
Returns the new element.
"""
function Base.push!(s::DisjointSets{T}, x::T) where T
    id = push!(s.internal)
    s.intmap[x] = id
    push!(s.revmap,x) # Note, this assumes invariant: length(s.revmap) == id
    return x
end


text:
using Flipping

test = "/home/beatriz/mainen.flipping.5ht@gmail.com/Flipping/run_task_photo/raw_data/DN1_170807a.csv"
poke = process_pokes(test)
first(poke[:,[:Streak,:Poke,:LastPoke]],10)
union(poke[:,:Protocol])
##


text:
function aft_l(ϕβ, rr)
    update_Xβ_dist!(rr, ϕβ)
    # compute likelihood
    nll = -compute_loglik!(rr)
    return nll
end

function aft_h!(grad, hes, ϕβ, rr, degreetype)
    update_Xβ_dist!(rr, ϕβ)
    # Compute auxiliary vector for derivatives.
    int_coefs = IntCoefs(rr.dist, degreetype)

    # Compute likelihoog, grad and hes
    nll = -compute_loglik!(rr)
    subtract_ders!(grad,hes, rr, int_coefs)
    return nll
end

aft{T<:Real}(S::AbstractVector{Event{T}},X::AbstractArray, pdist, degreetype; kwargs...) =
aft(EventWindow.(S),X::AbstractArray, pdist, degreetype; kwargs...)

function aft(S::AbstractVector,X::AbstractArray, pdist, degreetype; kwargs...)
    rr = AftResp(S, X, pdist)
    f1 = (ϕβ) -> aft_l(ϕβ, rr)
    h1! = (ϕβ,grad,hes) -> aft_h!(grad, hes, ϕβ, rr, degreetype)
    return newton_raphson(f1,h1!, initialize_aft(S, X, pdist, size(X,2)); kwargs...)
end

function aft(formula::Formula, data::DataFrame, pdist, degreetype = Val{50}(); kwargs...)
    M = DataFrames.ModelFrame(formula,data)
    S = convert(Array, M.df[:,1])
    model_matrix = DataFrames.ModelMatrix(M)
    X = model_matrix.m
    ϕβ, neg_ll,grad, hes =  aft(S, X, pdist, degreetype; kwargs...)
    rownms = vcat(string.(["params"],1:length(pdist.params)),coefnames(M))
    se = sqrt.(diag(pinv(hes)))
    z_score = ϕβ./se
    pvalues = 2*cdf(Normal(),-abs.(z_score))
    pdist.params[:] = ϕβ[1:length(pdist.params)]
    coefmat = CoefTable(hcat([ϕβ, se, z_score, pvalues]...),
    ["Estimate", "Std.Error", "z value", "Pr(>|z|)"], rownms, 4)
    AftModel("Accelerated Failure Time, dist = $pdist;\n", formula,
             coefmat, M, -neg_ll, -grad, hes, pdist)
end


text:
module PendulumSimulation

using Modia
@usingModiaPlot
using Test

Pendulum = Model(
   L = 0.8u"m",
   m = 1.0u"kg",
   d = 0.5u"N*m*s/rad",
   g = 9.81u"m/s^2",
   phi = Var(init = 1.57*u"rad"),
   w   = Var(init = 0u"rad/s"),
   equations = :[
          w = der(phi)
        0.0 = m*L^2*der(w) + d*w + m*g*L*sin(phi)
          r = [L*cos(phi), -L*sin(phi)]
   ]
)

pendulum = @instantiateModel(Pendulum)
simulate!(pendulum, Tsit5(), stopTime = 10.0u"s", log=true)
plot(pendulum, [("phi", "w"); "r"])

using Measurements
PendulumWithUncertainties = Pendulum | Map(L = (0.8 ± 0.2)u"m",
                                           m = (1.0 ± 0.2)u"kg",
                                           d = (0.5 ± 0.2)u"N*m*s/rad")

pendulum2 =  @instantiateModel(PendulumWithUncertainties,
                               FloatType = Measurement{Float64})

simulate!(pendulum2, Tsit5(), stopTime = 10.0u"s")
plot(pendulum2, [("phi", "w"); "r"], figure = 2)

# Linearize
using DoubleFloats
using Measurements
println("\n... Numerically linearize at stopTime = 10 with Float64 and Double64:")
(A_10, x_10) = linearize!(pendulum2, stopTime=10)

pendulum3 = SimulationModel{Measurement{Double64}}(pendulum2)
(A_10_Double64, x_10_Double64) = linearize!(pendulum3, stopTime=10)

xNames = get_xNames(pendulum2)
@show xNames
println(IOContext(stdout, :error_digits=>15), "A_10 = ", A_10, ", x_10 = ", x_10)
println(IOContext(stdout, :error_digits=>15), "A_10_Double64 = ", A_10_Double64, ", x_10_Double64 = ", x_10_Double64)

end


text:
abstract type ComponentType end
abstract type VecLike <:ComponentType end
struct Vec <:VecLike end
struct Pnt <:ComponentType end
struct Nrm <:VecLike end

name(x::Type{Vec}) = "Vector"
name(x::Type{Pnt}) = "Point"
name(x::Type{Nrm}) = "Normal"


text:
### Day 10
## Balance Bots
## Author: Thanasis Georgiou

# Shortcut to convert strings to integers
int(str) = parse(Int, str)

# Regex to find targets and numbers
targets = r"to (\w+)"
numbers = r"(\d+)"

# Dictionary of bots
bots = Dict{Int}{Array{Int}}()
function getBot(id)
  bot = get(bots, id, -1)
  if bot == -1
    bot = Array(Int, 0)
    bots[id] = bot
  end

  return bot
end

# Dictionary of bins
outputBins = Dict{Int}{Array{Int}}()
function getBin(id)
  bin = get(outputBins, id, -1)
  if bin == -1
    bin = Array(Int, 0)
    outputBins[id] = bin
  end

  return bin
end

# A command is a string plus a boolean indicated if this command is
# executed.
type Command
  string::String
  executed::Bool
end

# Read input
f = open("input.txt")
commands = map(line -> Command(chomp(line), false), readlines(f))
close(f)

# Keep executing while there are unexecuted commands
while true
  filteredCommands = filter(cmd -> cmd.executed == false, commands)
  if length(filteredCommands) == 0
    break
  end

  for cmd in filteredCommands
    # Execute command
    if cmd.string[1:5] == "value"
      # Grab useful info
      (chipId, targetBot) = map(
        x -> int(x),
        matchall(numbers, cmd.string)
      )

      # Find target bot and give chip
      bot = getBot(targetBot)
      push!(bot, chipId)

      # Mark command as executed
      cmd.executed = true
      continue
    end

    # Grab useful info
    (lowType, highType) = matchall(targets, cmd.string)
    (source, lowTarget, highTarget) = map(
      x -> int(x),
      matchall(numbers, cmd.string)
    )

    # Get source bot
    sourceBot = getBot(source)
    if (length(sourceBot) < 2)
      continue
    end

    # Answer problem 1
    if minimum(sourceBot) == 17 && maximum(sourceBot) == 61
      println("Problem 1: $source")
    end

    # Pass to output
    if lowType == "to output"
      bin = getBin(lowTarget)
      push!(bin, minimum(sourceBot))
    else
      bot = getBot(lowTarget)
      push!(bot, minimum(sourceBot))
    end

    if highType == "to output"
      bin = getBin(highTarget)
      push!(bin, maximum(sourceBot))
    else
      bot = getBot(highTarget)
      push!(bot, maximum(sourceBot))
    end

    # Mark as executed
    cmd.executed = true
  end
end

product = outputBins[0][1] * outputBins[1][1] * outputBins[2][1]
println("Problem 2: $product")


text:
invalid_imgs = [
    ("too_few_dimensions", MethodError, rand(127)),
    ("too_many_channels", AssertionError, rand(127, 257, 5)),
    ("too_many_dimensions", MethodError, rand(127, 257, 3, 1)),
]

@testset "invalid inputs" begin
for (case, exception, image) in invalid_imgs
        @testset "$(case) throws" begin
            @test_throws exception PNGFiles.save(
                joinpath(PNG_TEST_PATH, "test_img_err_$(case).png"),
                image
            )
        end
    end
end


text:
module FEMMShellT3FFModule

using LinearAlgebra: norm, Transpose, mul!, diag, eigen, I, dot, rank
using Statistics: mean
using FinEtools
import FinEtools.FESetModule: gradN!, nodesperelem, manifdim
using FinEtools.IntegDomainModule: IntegDomain, integrationdata, Jacobianvolume
import FinEtools.FEMMBaseModule: associategeometry!
import FinEtools.FEMMBaseModule: inspectintegpoints
using FinEtoolsDeforLinear.MatDeforLinearElasticModule: tangentmoduli!, update!, thermalstrain!
using FinEtools.MatrixUtilityModule: add_btdb_ut_only!, complete_lt!, locjac!, add_nnt_ut_only!, add_btsigma!
using ..FESetShellT3Module: FESetShellT3
using ..TransformerModule: QTEQTransformer, Layup2ElementAngle

const __nn = 3 # number of nodes
const __ndof = 6 # number of degrees of freedom per node

# Formulation for the transverse shear stiffness which averages the
# strain-displacement matrix.
const __TRANSV_SHEAR_FORMULATION_AVERAGE_B = 0
# Formulation for the transverse shear stiffness which averages the
# stiffness matrix.
const __TRANSV_SHEAR_FORMULATION_AVERAGE_K = 1

"""
    FEMMShellT3FF{S<:FESetT3, F<:Function} <: AbstractFEMM

Class for the finite element modeling machine of the T3 triangular Flat-Facet
shell with the Discrete Shear Gap technology and a consistent handling of the
normals.

With averaging of the transverse strain-displacement matrix or averaging of the
transverse shear stiffness matrix to provide isotropic transverse shear
response.

Also, the formulation is developed to correctly handle the coupling of twisting
moments and transverse shear (such as in the twisted beam or the Raasch hook
problems) by incorporating "nodal" normals.

Some of the programming developed consistently with the paper

[1] Cui et al, Analysis of plates and shells using an edge-based smoothed finite
element method, Comput Mech (2010) 45:141–156 DOI 10.1007/s00466-009-0429-9

The stabilization factor of the shear term of

[2] Mikko Lyly, Rolf Stenberg and Teemu Vihinen, A stable bilinear element for
the Reissner-Mindlin plate model Computer Methods in Applied Mechanics and
Engineering 110 (1993) 343-357

is incorporated. Refer to expressions (3.12) and (3.13).

The treatment of the transformation between the element and nodal coordinates
is carried out using a clean alternative to the publication

[3] Finite Elements in Analysis and Design 30 (1998) 235—242
The treatment of shell normals in ﬁnite element analysis
Richard H. MacNeal, Charles T. Wilson, Robert L. Harder, Claus C. Hoﬀ
The MacNeal-Schwendler Corporation, 815 Colorado Blvd., Los Angeles, CA 90041, USA

The formula for the element to nodal basis transformation is derived from the
expression

    [theta]_n = [A]_E^T [theta]_e + [A_3,1:3]_E^T [alpha_3]_e

by disconnecting the drilling degree of freedom from the bending degrees of
freedom in the nodal basis.

The following features are incorporated to deal with nodal normals:
- Nodal normals are averages of the normals of elements that meet at a node.
- A crease in the surface is taken into account. In that case the normal are not
  averaged across the crease. At the nodes along the crease every element uses
  the normal to its surface instead of the nodal normal.

Configuration:

These attributes of the FEMM can be set after it's been created.
- `transv_shear_formulation`: which formulation for the transverse shear stiffness?
    + `FEMMShellT3FFModule.__TRANSV_SHEAR_FORMULATION_AVERAGE_B` - averaged strains (default)
    + `FEMMShellT3FFModule.__TRANSV_SHEAR_FORMULATION_AVERAGE_K` - averaged stiffness
- `drilling_stiffness_scale`: multiplier of the generalized stiffness coefficient
- `threshold_angle`: angle in degrees. If a nodal normal subtends angle bigger
  then this threshold, the nodal normal at that note is marked as invalid.
- `mult_el_size`: multiplier of the square of the element size, used to control
  transverse shear stiffness.
"""
mutable struct FEMMShellT3FF{S<:FESetT3, F<:Function, M} <: AbstractFEMM
    integdomain::IntegDomain{S, F} # integration domain data
    mcsys::CSys # updater of the material orientation matrix
    material::M # material object.
    transv_shear_formulation::FInt
    drilling_stiffness_scale::Float64
    threshold_angle::Float64
    mult_el_size::Float64
    _associatedgeometry::Bool
    _normals::FFltMat
    _normal_valid::Vector{Bool}
    # The attributes below are buffers used in various operations.
    _loc::FFltMat
    _J0::FFltMat
    _ecoords::FFltMat
    _edisp::FFltVec
    _ecoords_e::FFltMat
    _edisp_e::FFltVec
    _dofnums::FIntMat
    _E_G::FFltMat
    _A_Es::Vector{FFltMat} # transformation nodal-element matrices
    _nvalid::Vector{Bool}
    _T::FFltMat  # element transformation matrix
    _elmat::FFltMat
    _gradN_e::FFltMat
    _Bm::FFltMat
    _Bb::FFltMat
    _Bs::FFltMat
    _DpsBmb::FFltMat
    _DtBs::FFltMat
end

"""
    FEMMShellT3FF(integdomain::IntegDomain{S, F}, mcsys::CSys, material::M) where {S<:FESetT3, F<:Function, M}

Constructor of the T3FF shell FEMM.
"""
function FEMMShellT3FF(integdomain::IntegDomain{S,F}, mcsys::CSys, material::M) where {S<:FESetT3,F<:Function,M}
    _nnmax = 0
    for j in 1:count(integdomain.fes)
        for k in eachindex(integdomain.fes.conn[j])
            _nnmax = max(_nnmax, integdomain.fes.conn[j][k])
        end
    end
    _normals = fill(0.0, _nnmax, 3)
    _normal_valid = fill(true, _nnmax)
    # Alocate the buffers
    _loc = fill(0.0, 1, 3)
    _J0 = fill(0.0, 3, 2)
    _ecoords = fill(0.0, __nn, 3)
    _edisp = fill(0.0, __nn * __ndof)
    _ecoords_e = fill(0.0, __nn, 2)
    _edisp_e = fill(0.0, __nn * __ndof)
    _dofnums = zeros(FInt, 1, __nn * __ndof)
    _E_G = fill(0.0, 3, 3)
    _A_Es = [fill(0.0, 3, 3), fill(0.0, 3, 3), fill(0.0, 3, 3)]
    _nvalid = fill(false, 3)
    _T = fill(0.0, __nn * __ndof, __nn * __ndof)
    _elmat = fill(0.0, __nn * __ndof, __nn * __ndof)
    _gradN_e = fill(0.0, __nn, 2)
    _Bm = fill(0.0, 3, __nn * __ndof)
    _Bb = fill(0.0, 3, __nn * __ndof)
    _Bs = fill(0.0, 2, __nn * __ndof)
    _DpsBmb = similar(_Bm)
    _DtBs = similar(_Bs)

    @assert delegateof(integdomain.fes) === FESetShellT3()

    return FEMMShellT3FF(integdomain, mcsys,
        material,
        # transv_shear_formulation::FInt
        # drilling_stiffness_scale::Float64
        # threshold_angle::Float64
        # mult_el_size::Float64
        __TRANSV_SHEAR_FORMULATION_AVERAGE_B, 1.0, 30.0, 5 / 12 / 1.5,
        false,
        _normals, _normal_valid,
        _loc, _J0,
        _ecoords, _edisp, _ecoords_e, _edisp_e,
        _dofnums,
        _E_G, _A_Es, _nvalid, _T,
        _elmat,
        _gradN_e,
        _Bm, _Bb, _Bs, _DpsBmb, _DtBs)
end

function isoparametric!(E_G::FFltMat, XYZ::FFltMat, J0::FFltMat, fe_label::FInt)
    return _e_g!(E_G, J0)
end

function _compute_nodal_normal!(n, mcsys::CSys, XYZ, J0::FFltMat, labl::FInt)
    updatecsmat!(mcsys, reshape(XYZ, 1, 3), J0, labl);
    n[:] .= mcsys.csmat[:, 3]
    return n
end

function FEMMShellT3FF(integdomain::IntegDomain{S, F}, material::M) where {S<:FESetT3, F<:Function, M}
    return FEMMShellT3FF(integdomain, CSys(3, 3, isoparametric!), material)
end

"""
    make(integdomain, mcsys, material)

Make a T3FF FEMM from the integration domain,  and a material.
Default isoparametric method for computing the normals is used.
"""
function make(integdomain, material)
    return FEMMShellT3FF(integdomain, material)
end

"""
    make(integdomain, mcsys, material)

Make a T3FF FEMM from the integration domain, a coordinate system to define the
orientation of the normals, and a material.
"""
function make(integdomain, mcsys, material)
    return FEMMShellT3FF(integdomain, mcsys, material)
end

function _compute_J!(J0, ecoords)
    # Compute the Jacobian matrix: the Jacobian matrix is constant within the
    # triangle
    # x, y, z = ecoords[2, :].-ecoords[1, :]
    # J0[:, 1] .= (x, y, z)
    # x, y, z = ecoords[3, :].-ecoords[1, :]
    # J0[:, 2] .= (x, y, z)
    for j in 1:3
        J0[j, 1] = ecoords[2, j] - ecoords[1, j]
        J0[j, 2] = ecoords[3, j] - ecoords[1, j]
    end
    return J0
end

function _e_g!(E_G, J0)
    # J0[:, 1] is the tangent to the coordinate curve 1
    for j in 1:3
        E_G[j,1]  = J0[j, 1]
    end
    n = sqrt(E_G[1, 1]^2 + E_G[2, 1]^2 + E_G[3, 1]^2)
    for j in 1:3
        E_G[j,1]  /= n
    end
    # J0[:, 2] is the tangent to the coordinate curve 2
    # Now compute the normal
    E_G[1, 3] = -E_G[3, 1]*J0[2, 2]+E_G[2, 1]*J0[3, 2]
    E_G[2, 3] =  E_G[3, 1]*J0[1, 2]-E_G[1, 1]*J0[3, 2]
    E_G[3, 3] = -E_G[2, 1]*J0[1, 2]+E_G[1, 1]*J0[2, 2]
    n = sqrt(E_G[1, 3]^2 + E_G[2, 3]^2 + E_G[3, 3]^2);
    for j in 1:3
        E_G[j, 3]  /= n
    end
    E_G[1, 2] = -E_G[3, 3]*E_G[2, 1]+E_G[2, 3]*E_G[3, 1]
    E_G[2, 2] =  E_G[3, 3]*E_G[1, 1]-E_G[1, 3]*E_G[3, 1]
    E_G[3, 2] = -E_G[2, 3]*E_G[1, 1]+E_G[1, 3]*E_G[2, 1]
    return E_G
end

function _ecoords_e!(ecoords_e, J0, E_G)
    ecoords_e[1, 1] = ecoords_e[1, 2] = 0.0
    ecoords_e[2, 1] = dot(view(J0, :, 1), view(E_G, :, 1))
    ecoords_e[2, 2] = dot(view(J0, :, 1), view(E_G, :, 2))
    ecoords_e[3, 1] = dot(view(J0, :, 2), view(E_G, :, 1))
    ecoords_e[3, 2] = dot(view(J0, :, 2), view(E_G, :, 2))
    return ecoords_e
end

function _centroid!(centroid, ecoords)
    centroid[1]  = (ecoords[1, 1]+ecoords[2, 1]+ecoords[3, 1])/3
    centroid[2]  = (ecoords[1, 2]+ecoords[2, 2]+ecoords[3, 2])/3
    centroid[3]  = (ecoords[1, 3]+ecoords[2, 3]+ecoords[3, 3])/3
    return centroid
end

function _shell_material_stiffness(material)
    # Compute the two material stiffness matrices: the plane stress matrix, and
    # the transverse shear matrix
    D = fill(0.0, 6, 6)
    t::FFlt, dt::FFlt, loc::FFltMat, label::FInt = 0.0, 0.0, [0.0 0.0 0.0], 0
    tangentmoduli!(material,  D,  t, dt, loc, label)
    Dps = fill(0.0, 3, 3)
    Dps[1:2, 1:2] = D[1:2, 1:2] -  (reshape(D[1:2,3], 2, 1) * reshape(D[3,1:2], 1, 2))/D[3, 3]
    ix=[1, 2, 4];
    for i = 1:3
        Dps[3,i] = Dps[i,3] = D[4, ix[i]];
    end
    Dt = fill(0.0, 2, 2)
    ix=[5, 6];
    for i = 1:2
        Dt[i,i] = D[ix[i], ix[i]];
    end
    return Dps, Dt
end

struct NodalTriadsE
    r::FFltVec
    nk_e::FFltVec
    nk::FFltVec
    f3_e::FFltVec
end

function NodalTriadsE()
    NodalTriadsE(fill(0.0, 3), fill(0.0, 3), fill(0.0, 3), [0.0, 0.0, 1.0])
end

(o::NodalTriadsE)(A_Es, nvalid, E_G, normals, normal_valid, c) = begin
    # Components of nodal cartesian ordinate systems such that the third
    # direction is the direction of the nodal normal, and the angle to rotate
    # the element normal into the nodal normal is as short as possible; these
    # components are given on the local element coordinate system basis
    # vectors (basis E). If the angle of rotation is excessively small, the
    # nodal matrix is the identity.
    #
    # The array `A_Es` holds the triads that are the nodal-to-element
    # transformation matrices. The array `nvalid` indicates for the three nodes
    # which of the normals is valid (`true`).

    for k in 1:length(A_Es)
        o.nk .= vec(view(normals, c[k], :))
        nvalid[k] = normal_valid[c[k]]
        if nvalid[k]
            # If the nodal normal is valid, pull it back into the element frame.
            mul!(o.nk_e, transpose(E_G), o.nk)
        else
            # Otherwise the element normal replaces the nodal normal.
            o.nk_e .= 0.0; o.nk_e[3] = 1.0
        end
        cross3!(o.r, o.f3_e, o.nk_e)
        if norm(o.r) > 1.0e-12
            rotmat3!(A_Es[k], o.r)
        else
            # A_Es[k] .= I(3)
            A_Es[k] .= 0.0
            A_Es[k][1, 1] = A_Es[k][2, 2] = A_Es[k][3, 3] = 1.0
        end
    end
    return A_Es, nvalid
end

struct TransfmatGToA
    Tblock::FFltMat
end

function TransfmatGToA()
    TransfmatGToA(fill(0.0, 3, 3))
end

(o::TransfmatGToA)(T, A_Es, E_G) = begin
    # Global-to-nodal transformation matrix.

    # The 3x3 blocks consist of the nodal triad expressed on the global basis.
    # The nodal basis vectors in the array `A_Es[i]` are expressed on the
    # element basis, and are then rotated with `E_G` into the global coordinate
    # system.

    # Output
    # - `T` = transformation matrix, input in the global basis, output in the
    #   nodal basis
    T .= 0.0
    for i in 1:__nn
        mul!(o.Tblock,  transpose(A_Es[i]), transpose(E_G))
        offset = (i-1)*__ndof
        r = offset+1:offset+3
        @. T[r, r] = o.Tblock
        r = offset+4:offset+6
        @. T[r, r] = o.Tblock
    end
    return T
end

function _transfmat_a_to_e!(T, A_Es, gradN_e)
    # Nodal-to-element transformation matrix.
    # - `A_Es` = matrix with the nodal triad vectors in columns, components on
    #   the element basis. Its transpose is the element triad on the nodal
    #   basis vectors.
    # - `gradN_e` = basis function gradients on the element basis.
    # Output
    # - `T` = transformation matrix, input in the nodal basis, output in the
    #   element basis: t_e = T * t_n
    # Rotation degrees of freedom: The drilling rotation of the mid surface
    # produced by the 1/2*(v,x - u,y) effect is linked to the out of plane
    # rotations.

    T .= 0.0
    for i in 1:__nn
        roffst = (i-1)*__ndof
        iA_E = A_Es[i]
        iA_E_33 = iA_E[3, 3]
        # T[r, r] .= iA_E
        for cl in 1:3
            for rw in 1:3
                T[roffst+rw, roffst+cl] = iA_E[rw, cl]
            end
        end
        for cl in 1:2
            for rw in 1:2
                T[roffst+3+rw, roffst+3+cl] = iA_E[rw, cl] - (1/iA_E_33)*iA_E[rw, 3]*iA_E[cl, 3]
            end
        end
        m1 = (1/iA_E_33)*iA_E[1, 3]
        m2 = (1/iA_E_33)*iA_E[2, 3]
        for j in 1:__nn
            coffst = (j-1)*__ndof
            for k in 1:3
                a3 = 1/2 * (iA_E[2, k] * gradN_e[j, 1] - iA_E[1, k] * gradN_e[j, 2])
                T[roffst+4, coffst+k] += m1 * a3
                T[roffst+5, coffst+k] += m2 * a3
            end
        end
    end
    return T
end

function _gradN_e_Ae!(gradN_e, ecoords_e)
    # Compute the gradient with respect to the element coordinates, three
    # gradients in two coordinates, and the area of the triangle
    # a, b = ecoords_e[2, :] .- ecoords_e[1, :]
    # c, d = ecoords_e[3, :] .- ecoords_e[1, :]
    # J = (a*d - b*c)
    # gradN_e[:] .= ((b-d)/J, d/J, -b/J, (c-a)/J, -c/J, a/J)
     a = ecoords_e[2, 1] - ecoords_e[1, 1]
     b = ecoords_e[2, 2] - ecoords_e[1, 2]
     c = ecoords_e[3, 1] - ecoords_e[1, 1]
     d = ecoords_e[3, 2] - ecoords_e[1, 2]
     J = (a*d - b*c)
     gradN_e[1, 1] = (b-d)/J
     gradN_e[2, 1] = d/J
     gradN_e[3, 1] = -b/J
     gradN_e[1, 2] = (c-a)/J
     gradN_e[2, 2] = -c/J
     gradN_e[3, 2] = a/J
    return gradN_e, J/2
end

function _add_Bsmat_o!(Bs, ecoords_e, Ae, ordering)
    # Compute the linear transverse shear strain-displacement matrix for one
    # particular ordering of the nodes. The computed entries are ADDED. The
    # matrix is NOT zeroed out initially. That is the responsibility of the
    # caller.

    # Orientation
    s, p, q = ordering
    a = ecoords_e[p, 1] - ecoords_e[s, 1]
    b = ecoords_e[p, 2] - ecoords_e[s, 2]
    c = ecoords_e[q, 1] - ecoords_e[s, 1]
    d = ecoords_e[q, 2] - ecoords_e[s, 2]
    m = (1/2/Ae) # multiplier

    # The first node in the triangle
    # Node s
    co = (s - 1) * 6 # column offset
    Bs[1, co+3] += m*(b-d);                             Bs[1, co+5] += m*(Ae)
    Bs[2, co+3] += m*(c-a); Bs[2, co+4] += m*(-Ae);
    # The other two nodes
    # Node p
    co = (p - 1) * 6 # column offset
    Bs[1, co+3] += m*(d);   Bs[1, co+4] += m*(-b*d/2);  Bs[1, co+5] += m*(a*d/2)
    Bs[2, co+3] += m*(-c);  Bs[2, co+4] += m*(b*c/2);   Bs[2, co+5] += m*(-a*c/2)
    # Node q
    co = (q - 1) * 6 # column offset
    Bs[1, co+3] += m*(-b);  Bs[1, co+4] += m*(b*d/2);   Bs[1, co+5] += m*(-b*c/2)
    Bs[2, co+3] += m*(a);   Bs[2, co+4] += m*(-a*d/2);  Bs[2, co+5] += m*(a*c/2)

    return Bs
end

function _Bsmat!(Bs, ecoords_e, Ae)
    # Compute the linear transverse shear strain-displacement matrix.
    Bs .= 0.0 # Zero out initially, then add the three contributions
    _add_Bsmat_o!(Bs, ecoords_e, Ae, (1, 2, 3))
    _add_Bsmat_o!(Bs, ecoords_e, Ae, (2, 3, 1))
    _add_Bsmat_o!(Bs, ecoords_e, Ae, (3, 1, 2))
    # Now averaging of the three contributions
    Bs .*= (1/3)
    return Bs
end

function _Bmmat!(Bm, gradN)
    # Compute the linear membrane strain-displacement matrix.
    fill!(Bm, 0.0)
    for i in 1:__nn
        Bm[1,6*(i-1)+1] = gradN[i,1];
                                             Bm[2,6*(i-1)+2] = gradN[i,2];
        Bm[3,6*(i-1)+1] = gradN[i,2];        Bm[3,6*(i-1)+2] = gradN[i,1];
    end
end

function _Bbmat!(Bb, gradN)
    # Compute the linear, displacement independent, curvature-displacement/rotation
    # matrix for a shell quadrilateral element with nfens=3 nodes. Displacements and
    # rotations are in a local coordinate system.
    fill!(Bb, 0.0)
    for i in 1:__nn
                                              Bb[1,6*(i-1)+5] = gradN[i,1];
        Bb[2,6*(i-1)+4] = -gradN[i,2];
        Bb[3,6*(i-1)+4] = -gradN[i,1];        Bb[3,6*(i-1)+5] = gradN[i,2];
    end
end

"""
    associategeometry!(self::FEMMShellT3FF,  geom::NodalField{FFlt})

Associate geometry with the FEMM.

In this case it means evaluate the nodal normals.
"""
function associategeometry!(self::FEMMShellT3FF,  geom::NodalField{FFlt})
    threshold_angle = self.threshold_angle
    J0 = self._J0
    E_G = self._E_G
    normals, normal_valid = self._normals, self._normal_valid
    nnormal = fill(0.0, 3)

    # Compute the normals at the nodes
    for el in 1:count(self.integdomain.fes)
        i, j, k = self.integdomain.fes.conn[el]
        J0[:, 1] = geom.values[j, :] - geom.values[i, :]
        J0[:, 2] = geom.values[k, :] - geom.values[i, :]
        for n in [i, j, k]
            _compute_nodal_normal!(nnormal, self.mcsys, geom.values[n, :], J0, self.integdomain.fes.label[el])
            normals[n, :] .+= nnormal
        end
    end

    # Normalize to unit length
    for j in 1:size(normals, 1)
        nn = norm(normals[j, :])
        if nn > 0.0
            normals[j, :] ./= nn
        end
    end

    # Now perform a second pass. If the nodal normal differs by a substantial
    # amount from the element normals, zero out the nodal normal to indicate
    # that the nodal normal should not be used at that vertex.
    ntolerance = 1 - sqrt(1 - sin(threshold_angle/180*pi)^2)
    for el in 1:count(self.integdomain.fes)
        i, j, k = self.integdomain.fes.conn[el]
        _compute_J!(J0, geom.values[[i, j, k], :])
        for n in [i, j, k]
            _e_g!(E_G, J0)
            nd = dot(normals[n, :], E_G[:, 3])
            if nd < 1 - ntolerance
                normal_valid[n] = false
            end
        end
    end

    # Mark the engine as being associated with geometry
    self._associatedgeometry = true
    return self
end

"""
    num_normals(self::FEMMShellT3FF)

Compute the summary of the nodal normals.
"""
function num_normals(self::FEMMShellT3FF)
    normals, normal_valid = self._normals, self._normal_valid
    total_normals = length(normal_valid)
    total_invalid_normals = length([v for v in normal_valid if v != true])
    return total_normals, total_invalid_normals
end

"""
    stiffness(self::FEMMShellT3FF, assembler::ASS, geom0::NodalField{FFlt}, u1::NodalField{TI}, Rfield1::NodalField{TI}, dchi::NodalField{TI}) where {ASS<:AbstractSysmatAssembler, TI<:Number}

Compute the material stiffness matrix.
"""
function stiffness(self::FEMMShellT3FF, assembler::ASS, geom0::NodalField{FFlt}, u1::NodalField{TI}, Rfield1::NodalField{TI}, dchi::NodalField{TI}) where {ASS<:AbstractSysmatAssembler, TI<:Number}
    @assert self._associatedgeometry == true
    fes = self.integdomain.fes
    label = self.integdomain.fes.label
    normals, normal_valid = self._normals, self._normal_valid
    centroid, J0 = self._loc, self._J0
    ecoords, edisp, dofnums = self._ecoords, self._edisp, self._dofnums
    ecoords_e, gradN_e = self._ecoords_e, self._gradN_e
    E_G, A_Es, nvalid, T = self._E_G, self._A_Es, self._nvalid, self._T
    elmat = self._elmat
    transformwith = QTEQTransformer(elmat)
    _nodal_triads_e! = NodalTriadsE()
    _transfmat_g_to_a! = TransfmatGToA()
    Bm, Bb, Bs, DpsBmb, DtBs = self._Bm, self._Bb, self._Bs, self._DpsBmb, self._DtBs
    Dps, Dt = _shell_material_stiffness(self.material)
    scf = 5/6;  # shear correction factor
    Dt .*= scf
    ipc = [1.0/3 1.0/3]
    drilling_stiffness_scale = self.drilling_stiffness_scale
    transv_shear_formulation = self.transv_shear_formulation
    mult_el_size = self.mult_el_size
    startassembly!(assembler, size(elmat, 1), size(elmat, 2), count(fes), dchi.nfreedofs, dchi.nfreedofs);
    for i in 1:count(fes) # Loop over elements
        gathervalues_asmat!(geom0, ecoords, fes.conn[i]);
        _centroid!(centroid, ecoords)
        _compute_J!(J0, ecoords)
        _e_g!(E_G, J0)
        _ecoords_e!(ecoords_e, J0, E_G)
        gradN_e, Ae = _gradN_e_Ae!(gradN_e, ecoords_e)
        t = self.integdomain.otherdimension(centroid, fes.conn[i], ipc)
        # Construct the Stiffness Matrix
        fill!(elmat,  0.0); # Initialize element matrix
        _Bmmat!(Bm, gradN_e)
        add_btdb_ut_only!(elmat, Bm, t*Ae, Dps, DpsBmb)
        _Bbmat!(Bb, gradN_e)
        add_btdb_ut_only!(elmat, Bb, (t^3)/12*Ae, Dps, DpsBmb)
        # he = sqrt(2*Ae) # we avoid taking the square root here, replacing he^2 with 2*Ae
        if transv_shear_formulation == __TRANSV_SHEAR_FORMULATION_AVERAGE_K
            for  o in [(1, 2, 3), (2, 3, 1), (3, 1, 2)]
                Bs .= 0.0; _add_Bsmat_o!(Bs, ecoords_e, Ae, o)
                add_btdb_ut_only!(elmat, Bs, (t^3/(t^2+mult_el_size*2*Ae))*Ae/3, Dt, DtBs)
            end
        else
            _Bsmat!(Bs, ecoords_e, Ae)
            add_btdb_ut_only!(elmat, Bs, (t^3/(t^2+mult_el_size*2*Ae))*Ae, Dt, DtBs)
        end
        # Complete the elementwise matrix by filling in the lower triangle
        complete_lt!(elmat)
        # Transformation from the nodal (A) to the element (E) basis
        _nodal_triads_e!(A_Es, nvalid, E_G, normals, normal_valid, fes.conn[i])
        _transfmat_a_to_e!(T, A_Es, gradN_e)
        transformwith(elmat, T)
        # Bending diagonal stiffness coefficients
        kavg = mean((elmat[4, 4], elmat[10, 10], elmat[16, 16],
            elmat[5, 5], elmat[11, 11], elmat[17, 17])) * drilling_stiffness_scale
        # Add the artificial drilling stiffness in the nodal basis, but only if
        # the nodal normal is valid.
        nvalid[1] && (elmat[6,6] += kavg)
        nvalid[2] && (elmat[12,12] += kavg)
        nvalid[3] && (elmat[18,18] += kavg)
        # Transform from global (G) into the nodal (A) basis
        _transfmat_g_to_a!(T, A_Es, E_G)
        transformwith(elmat, T)
        # Assembly
        gatherdofnums!(dchi, dofnums, fes.conn[i]);
        assemble!(assembler, elmat, dofnums, dofnums);
    end # Loop over elements
    return makematrix!(assembler);
end

function stiffness(self::FEMMShellT3FF, geom0::NodalField{FFlt}, u1::NodalField{TI}, Rfield1::NodalField{TI}, dchi::NodalField{TI}) where {TI<:Number}
    assembler = SysmatAssemblerSparseSymm();
    return stiffness(self, assembler, geom0, u1, Rfield1, dchi);
end

"""
    mass(self::FEMMShellT3FF,  assembler::A,  geom::NodalField{FFlt}, dchi::NodalField{TI}) where {A<:AbstractSysmatAssembler, TI<:Number}

Compute the diagonal (lumped) mass matrix

The mass matrix can be expected to be non-singular.
"""
function mass(self::FEMMShellT3FF,  assembler::A,  geom0::NodalField{FFlt}, dchi::NodalField{TI}) where {A<:AbstractSysmatAssembler, TI<:Number}
    @assert self._associatedgeometry == true
    fes = self.integdomain.fes
    label = self.integdomain.fes.label
    normals, normal_valid = self._normals, self._normal_valid
    centroid, J0 = self._loc, self._J0
    ecoords, edisp, dofnums = self._ecoords, self._edisp, self._dofnums
    ecoords_e, gradN_e = self._ecoords_e, self._gradN_e
    E_G, A_Es, nvalid, T = self._E_G, self._A_Es, self._nvalid, self._T
    elmat = self._elmat
    rho::FFlt = massdensity(self.material); # mass density
    Dps, Dts = _shell_material_stiffness(self.material)
    npe = nodesperelem(fes)
    ndn = ndofs(dchi)
    ipc = [1.0/3 1.0/3]
    startassembly!(assembler,  size(elmat,1),  size(elmat,2),  count(fes), dchi.nfreedofs,  dchi.nfreedofs);
    for i = 1:count(fes) # Loop over elements
        gathervalues_asmat!(geom0, ecoords, fes.conn[i]);
        _centroid!(centroid, ecoords)
        _compute_J!(J0, ecoords)
        _e_g!(E_G, J0)
        _ecoords_e!(ecoords_e, J0, E_G)
        gradN_e, Ae = _gradN_e_Ae!(gradN_e, ecoords_e)
        # Compute the translational and rotational masses corresponding to nodes
        t = self.integdomain.otherdimension(centroid, fes.conn[i], ipc)
        tmass = rho*(t*Ae)/3;
        rmass = rho*(t^3/12*Ae)/3;
        # end # Loop over quadrature points
        # Fill the elementwise matrix in the global basis.
        fill!(elmat,  0.0); # Initialize element matrix
        for k in 1:npe
            # Translation degrees of freedom
            for d in 1:3
                c = (k - 1) * __ndof + d
                elmat[c, c] += tmass
            end
            # Bending degrees of freedom
            for d in 4:6
                c = (k - 1) * __ndof + d
                elmat[c, c] += rmass
            end
        end
        # Assemble
        gatherdofnums!(dchi,  dofnums,  fes.conn[i]);# retrieve degrees of freedom
        assemble!(assembler,  elmat,  dofnums,  dofnums);# assemble symmetric matrix
    end # Loop over elements
    return makematrix!(assembler);
end

function mass(self::FEMMShellT3FF,  geom::NodalField{FFlt},  u::NodalField{TI}) where {TI<:Number}
    assembler = SysmatAssemblerSparseSymm();
    return mass(self, assembler, geom, u);
end

"""
    inspectintegpoints(self::AbstractFEMMDeforLinear,
      geom::NodalField{FFlt},  u::NodalField{TI},
      dT::NodalField{FFlt},
      felist::FIntVec,
      inspector::F,  idat, quantity=:Cauchy;
      context...) where {TI<:Number, F<:Function}

Inspect integration point quantities.

- `geom` - reference geometry field
- `u` - displacement+rotation field
- `dT` - temperature difference field
- `felist` - indexes of the finite elements that are to be inspected:
     The fes to be included are: `fes[felist]`.
- `context`    - structure: see the update!() method of the material.
- `inspector` - function with the signature
        `idat = inspector(idat, j, conn, x, out, loc);`
   where
    `idat` - a structure or an array that the inspector may
           use to maintain some state,  for instance minimum or maximum of
           stress, `j` is the element number, `conn` is the element connectivity,
           `out` is the output of the update!() method,  `loc` is the location
           of the integration point in the *reference* configuration.
### Return
The updated inspector data is returned.
"""
function inspectintegpoints(self::FEMMShellT3FF, geom0::NodalField{FFlt},  u::NodalField{TI}, dT::NodalField{FFlt}, felist::FIntVec, inspector::F, idat, quantity=:Cauchy; context...) where {TI<:Number, F<:Function}
    @assert self._associatedgeometry == true
    fes = self.integdomain.fes
    label = self.integdomain.fes.label
    normals, normal_valid = self._normals, self._normal_valid
    centroid, J0 = self._loc, self._J0
    ecoords, edisp, dofnums = self._ecoords, self._edisp, self._dofnums
    ecoords_e, gradN_e = self._ecoords_e, self._gradN_e
    E_G, A_Es, nvalid, T = self._E_G, self._A_Es, self._nvalid, self._T
    elmat = self._elmat
    _nodal_triads_e! = NodalTriadsE()
    _transfmat_g_to_a! = TransfmatGToA()
    transformwith = QTEQTransformer(elmat)
    lla = Layup2ElementAngle()
    Bm, Bb, Bs, DpsBmb, DtBs = self._Bm, self._Bb, self._Bs, self._DpsBmb, self._DtBs
    Dps, Dt = _shell_material_stiffness(self.material)
    scf = 5/6;  # shear correction factor
    Dt .*= scf
    ipc = [1.0/3 1.0/3]
    mult_el_size = self.mult_el_size
    edisp_e = deepcopy(edisp)
    edisp_n = deepcopy(edisp_e)
    out = fill(0.0, 3)
    o2_e = fill(0.0, 2, 2)
    # Sort out  the output requirements
    outputcsys = self.mcsys; # default: report the stresses in the material coord system
    for apair in pairs(context)
        sy, val = apair
        if sy == :outputcsys
            outputcsys = val
        end
    end
    BENDING_MOMENT, TRANSVERSE_SHEAR, MEMBRANE_FORCE = 1, 2, 3
    quant = BENDING_MOMENT
    if quantity == :bending || quantity == :moment || quantity == :bending_moment
        quant = BENDING_MOMENT
    end
    if quantity == :transverse_shear || quantity == :transverse || quantity == :shear
        quant = TRANSVERSE_SHEAR
    end
    if quantity == :membrane_force || quantity == :membrane
        quant = MEMBRANE_FORCE
    end
    # Loop over  all the elements and all the quadrature points within them
    for ilist = 1:length(felist) # Loop over elements
        i = felist[ilist];
        gathervalues_asmat!(geom0, ecoords, fes.conn[i]);
        gathervalues_asvec!(u, edisp, fes.conn[i]);
        _centroid!(centroid, ecoords)
        _compute_J!(J0, ecoords)
        _e_g!(E_G, J0)
        _ecoords_e!(ecoords_e, J0, E_G)
        gradN_e, Ae = _gradN_e_Ae!(gradN_e, ecoords_e)
        t = self.integdomain.otherdimension(centroid, fes.conn[i], ipc)
        # Establish nodal triads
        _nodal_triads_e!(A_Es, nvalid, E_G, normals, normal_valid, fes.conn[i])
        # Transform from global into nodal coordinates
        _transfmat_g_to_a!(T, A_Es, E_G)
        mul!(edisp_n, T, edisp)
        # Now treat the transformation from the nodal to the element triad
        _transfmat_a_to_e!(T, A_Es, gradN_e)
         # Transform the nodal vector into the elementwise coordinates
        mul!(edisp_e, T, edisp_n)
        updatecsmat!(outputcsys, centroid, J0, fes.label[i]);
        if dot(view(outputcsys.csmat, :, 3), view(E_G, :, 3)) < 0.95
            @warn "Coordinate systems mismatched?"
        end
        ocsm, ocsn = lla(E_G, outputcsys.csmat)
        o2_e[1, 1] = o2_e[2, 2]  = ocsm
        o2_e[1, 2] = ocsn; o2_e[2, 1] = -ocsn
        # Compute the Requested Quantity
        if quant == BENDING_MOMENT
            _Bbmat!(Bb, gradN_e)
            kurv = Bb * edisp_e
            mom = ((t^3)/12)*Dps * kurv
            m = [mom[1] mom[3]; mom[3] mom[2]]
            mo = o2_e' * m * o2_e
            out[:] .= mo[1, 1], mo[2, 2], mo[1, 2]
        end
        if quant == MEMBRANE_FORCE
            _Bmmat!(Bm, gradN_e)
            strn = Bm * edisp_e
            frc = (t)*Dps * strn
            f = [frc[1] frc[3]; frc[3] frc[2]]
            fo = o2_e' * f * o2_e
            out[:] .= fo[1, 1], fo[2, 2], fo[1, 2]
            # @infiltrate
        end
        if quant == TRANSVERSE_SHEAR
            _Bsmat!(Bs, ecoords_e, Ae)
            he = sqrt(2*Ae)
            shr = Bs * edisp_e
            frc = ((t^3/(t^2+mult_el_size*2*Ae)))*Dt * shr
            fo = o2_e' * frc
            out[1:2] .= fo[1], fo[2]
        end
        # Call the inspector
        idat = inspector(idat, i, fes.conn[i], ecoords, out, centroid);
        # end # Loop over quadrature points
    end # Loop over elements
    return idat; # return the updated inspector data
end

function inspectintegpoints(self::FEMMShellT3FF, geom::NodalField{FFlt},  u::NodalField{TI}, felist::FIntVec, inspector::F, idat, quantity=:Cauchy; context...) where {TI<:Number, F<:Function}
    dT = NodalField(fill(zero(FFlt), nnodes(geom), 1)) # zero difference in temperature
    return inspectintegpoints(self, geom, u, dT, felist, inspector, idat, quantity; context...);
end

end # module