text:
using PowerSimulationsDynamics
using Sundials

"""
Case 23:
This case study a 15-state droop grid forming inverter against an infinite bus located at bus 1, with inverter located at bus 2.
The perturbation increase the reference power (analogy for mechanical power) from 0.5 to 0.7.
"""

##################################################
############### LOAD DATA ########################
##################################################

include(joinpath(dirname(@__FILE__), "data_tests/test23.jl"))

##################################################
############### SOLVE PROBLEM ####################
##################################################

#time span
tspan = (0.0, 4.0);

#Define Fault using Callbacks
Pref_change = ControlReferenceChange(1.0, case_inv, PSID.P_ref_index, 0.7)

path = (joinpath(pwd(), "test-23"))
!isdir(path) && mkdir(path)
try
    #Define Simulation Problem
    sim = Simulation!(
        ImplicitModel,
        omib_sys, # system
        path,
        tspan,
        Pref_change,
    )

    small_sig = small_signal_analysis(sim)

    #Solve problem in equilibrium
    execute!(sim, Sundials.IDA())

    #Obtain data for angles
    series = get_state_series(sim, ("generator-102-1", :θ_oc))

    print_device_states(sim)

    diff = [0.0]
    res = get_init_values_for_comparison(sim)
    for (k, v) in test23_x0_init
        diff[1] += LinearAlgebra.norm(res[k] - v)
    end
    @test (diff[1] < 1e-3)
    @test sim.solution.retcode == :Success
    @test small_sig.stable
finally
    @info("removing test files")
    rm(path, force = true, recursive = true)
end


text:
module ego_dataset

using NPZ
using Random: randperm

export EgoDataset, compile_data

struct EgoDataset
    T::Int
    num_frames::Int
    episodes::Array
    episode_map::Dict
    idx_map::Dict

    function EgoDataset(data_dir::AbstractString, T = 10)
        num_frames = 0
        episodes = []
        episode_map = Dict{Int,Int}()
        idx_map = Dict{Int,Int}()

        for fullpath in readdir(data_dir, join = true)
            len = length(readdir(fullpath)) ÷ 4 # folder contains loc,rot,spd,act for each frame
            @show len
            if len < T
                continue
            end

            locs = zeros(len, 3) # location 3-element vector
            rots = zeros(len, 3) # rotation 3-element vector
            spds = zeros(len, 1) # speed 1-element
            acts = zeros(len, 3) # action 3-element vector

            for file in readdir(fullpath, join = true)
                base = basename(file)
                idx = parse(Int, base[end-8:end-4]) + 1
                if startswith(base, "act")
                    acts[idx, :] = npzread(file)
                elseif startswith(base, "loc")
                    locs[idx, :] = npzread(file)
                elseif startswith(base, "rot")
                    rots[idx, :] = npzread(file)
                elseif startswith(base, "spd")
                    spds[idx, 1] = npzread(file)
                end
            end

            offset = num_frames
            for i = 1:len-T+1
                num_frames += 1
                # each index for the current batch maps to the current episode in episodes list.
                # length(episodes)+1 because we will add the episode just in a sec.
                episode_map[offset+i] = length(episodes) + 1
                idx_map[offset+i] = i
            end
            push!(episodes, (locs, rots, spds, acts))
        end
        new(T, num_frames, episodes, episode_map, idx_map)
    end
end

function Base.length(d::EgoDataset)
    return d.num_frames
end

function Base.getindex(d::EgoDataset, idx::Int)
    episode_idx = d.episode_map[idx]
    index = d.idx_map[idx]

    episode = d.episodes[episode_idx]
    locs = episode[1][index:index+d.T-1, :]
    rots = episode[2][index:index+d.T-1, :]
    spds = episode[3][index:index+d.T-1, :]
    acts = episode[4][index:index+d.T-1, :]
    return locs, rots, spds, acts
end

mutable struct EgoDatasetLoader
    ego_dataset::EgoDataset
    batchsize::Int
    shuffle::Bool
    imax
    indices

    function EgoDatasetLoader(ego_dataset::EgoDataset, batchsize; shuffle=false)
        imax = length(ego_dataset) - (length(ego_dataset) % batchsize)
        new(ego_dataset, batchsize, shuffle, imax, 1:length(ego_dataset))
    end
end

function Base.iterate(d::EgoDatasetLoader, i=1)
    len = length(d.ego_dataset)
    if i > d.imax
        return nothing
    end

    if d.shuffle && i==1
        d.indices = randperm(len)
    end
    nexti = min(i+d.batchsize-1, len)
    ids = d.indices[i:nexti]
    @show ids

    locs_batch = zeros(d.ego_dataset.T, 3, d.batchsize)
    rots_batch = zeros(d.ego_dataset.T, 3, d.batchsize)
    spds_batch = zeros(d.ego_dataset.T, 1, d.batchsize)
    acts_batch = zeros(d.ego_dataset.T, 3, d.batchsize)

    for idx = ids
        locs, rots, spds, acts = d.ego_dataset[idx]
        batch_idx = (idx-1)%16 + 1
        locs_batch[:,:,batch_idx] = locs
        rots_batch[:,:,batch_idx] = rots
        spds_batch[:,:,batch_idx] = spds
        acts_batch[:,:,batch_idx] = acts
    end

    return ((locs_batch, rots_batch, spds_batch, acts_batch), nexti+1)
end

function compile_data(args)
    ego_data = EgoDataset(args.data_dir, args.traj_len)
    return EgoDatasetLoader(ego_data, args.batchsize; shuffle=args.shuffle)
end

ext: jl


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*--------------------------------------------------------------------
#* Author: Tiago Os��io
#* $Id$
#*-------------------------------------------------------------------
type PID
	PID()=begin
		[
			:(PID_intTerm = 0),
			:(diff(PID_dFilt) = 0/"s"),
			:(diff(PID_inputFilt) = 0/"s"),
			:(diff(PID_setPointFilt) = 0/"s"),
		],
		[
			"","","","",
		],
		new(
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Type of PID Incremental",
				:Valid=>["Ideal","Parallel","Series","Ideal_AWBT","Parallel_AWBT","Series_AWBT","Ideal_AW","Parallel_AW","Series_AW"],
				:Default=>"Ideal"
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Controller action",
				:Valid=>["Direct","Reverse"],
				:Default=>"Reverse"
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Controller mode",
				:Valid=>["Automatic","Manual"],
				:Default=>"Automatic"
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Controller mode",
				:Valid=>["Clipped","Unclipped"],
				:Default=>"Clipped"
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Derivative term filter constant",
				:Default=>1
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Proportional term setPoint change filter"
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Previous scaled bias",
				:Default=>0.5
			)),
			time_sec(Dict{Symbol,Any}(
				:Brief=>"Derivative time constant"
			)),
			time_sec(Dict{Symbol,Any}(
				:Brief=>"Integral time constant"
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Controller gain",
				:Default=>0.5
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Derivative term SP change filter"
			)),
			time_sec(Dict{Symbol,Any}(
				:Brief=>"Input filter time constant"
			)),
			time_sec(Dict{Symbol,Any}(
				:Brief=>"Input filter time constant"
			)),
			control_signal(Dict{Symbol,Any}(
				:Default=>0
			)),
			control_signal(Dict{Symbol,Any}(
				:Default=>1000
			)),
			control_signal(Dict{Symbol,Any}(
				:Default=>0
			)),
			control_signal(Dict{Symbol,Any}(
				:Default=>1
			)),
			control_signal(Dict{Symbol,Any}(
				:Protected=>true,
				:PosX=>0,
				:PosY=>0.5
			)),
			control_signal(Dict{Symbol,Any}(
				:Protected=>true,
				:PosX=>0.54,
				:PosY=>1
			)),
			control_signal(),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Derivative term",
				:Hidden=>true ,
				:Default=>0
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Derivative term filtered",
				:Hidden=>true ,
				:Default=>0.5
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Error definition for proportional term",
				:Hidden=>true
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Error definition for derivative term",
				:Hidden=>true
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Error definition for integral term",
				:Hidden=>true
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Filtered input",
				:Hidden=>true
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Integral term",
				:Hidden=>true ,
				:Default=>0
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Sum of proportional, integral and derivative terms",
				:Hidden=>true
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Variable outp scaled between -1 and 1",
				:Hidden=>true
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Proportional term",
				:Default=>0 ,
				:Hidden=>true
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Filtered setPoint",
				:Default=>0,
				:Hidden=>true
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Integral term multiplier used in anti-reset windup",
				:Hidden=>true
			)),
			DanaReal(Dict{Symbol,Any}(
				:Protected=>true,
				:Hidden=>true
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Previous scaled input signal",
				:Default=>0.5,
				:Hidden=>true
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Scaled output signal",
				:Default=>0.5,
				:Hidden=>true
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Scaled setPoint",
				:Default=>0.5,
				:Hidden=>true
			)),
			[
				:(PID_input*(MaxInput - MinInput) = Input - MinInput),
				:(Output = PID_output*(MaxOutput - MinOutput) + MinOutput),
				:(PID_setPoint*(MaxInput - MinInput) = SetPoint - MinInput),
				:((tau + 1e-3*"s")*diff(PID_inputFilt)= PID_input - PID_inputFilt),
				:(tau*diff(PID_inputFilt)= PID_input - PID_inputFilt),
				:((tauSet + 1e-3*"s")*diff(PID_setPointFilt)= PID_setPoint - PID_setPointFilt),
				:(tauSet*diff(PID_setPointFilt)= PID_setPoint - PID_setPointFilt),
				:(PID_error = PID_inputFilt*(beta-1.0)),
				:(PID_errorD= PID_inputFilt*(gamma-1.0)),
				:(PID_errorI= 0),
				:(PID_error = beta*PID_setPointFilt - PID_inputFilt),
				:(PID_errorD = gamma*PID_setPointFilt - PID_inputFilt),
				:(PID_errorI = PID_setPointFilt-PID_inputFilt),
				:(PID_propTerm=PID_error),
				:(alpha*(derivTime + 1e-3*"s")*diff(PID_dFilt) = PID_errorD - PID_dFilt),
				:(alpha*(derivTime)*diff(PID_dFilt) = PID_errorD - PID_dFilt),
				:(PID_derivTerm = derivTime*diff(PID_dFilt)),
				:(PID_outps=2*PID_outp-1),
				:(PID_output=(sign(PID_outps)+1)/2),
				:(PID_output=PID_outp),
				:(PID_output=PID_outp),
				:(PID_action = -1.0),
				:(PID_action = 1.0),
				:(intTime*diff(PID_intTerm) = PID_AWFactor*PID_errorI),
				:(PID_outp = bias + PID_action*gain*(PID_propTerm + PID_intTerm + PID_derivTerm)),
				:(PID_AWFactor=-tanh(sign(PID_outps)*PID_outps*100-102)),
				:(PID_AWFactor=1),
				:(intTime*diff(PID_intTerm) = PID_AWFactor*PID_errorI),
				:(PID_outp = bias + PID_action*(gain*PID_propTerm + PID_intTerm + PID_derivTerm)),
				:(PID_AWFactor=-tanh(sign(PID_outps)*PID_outps*100-102)),
				:(PID_AWFactor=1),
				:(intTime*diff(PID_intTerm) = PID_AWFactor*PID_errorI),
				:(PID_outp = bias + PID_action*(gain*(PID_propTerm + PID_intTerm)*(1 + PID_derivTerm))),
				:(PID_AWFactor=-tanh(sign(PID_outps)*PID_outps*100-102)),
				:(PID_AWFactor=1),
				:(intTime*diff(PID_intTerm) = PID_errorI),
				:(PID_outp = bias + PID_action*gain*(PID_propTerm + PID_intTerm + PID_derivTerm)),
				:(PID_AWFactor=1),
				:(intTime*diff(PID_intTerm) = PID_errorI),
				:(PID_outp = bias + PID_action*(gain*PID_propTerm + PID_intTerm + PID_derivTerm)),
				:(PID_AWFactor=1),
				:(intTime*diff(PID_intTerm) = PID_errorI),
				:(PID_outp = bias + PID_action*(gain*(PID_propTerm + PID_intTerm)*(1 + PID_derivTerm))),
				:(PID_AWFactor=1),
				:(PID_action*gain*(intTime*diff(PID_intTerm)-PID_errorI) = PID_output-PID_outp),
				:(PID_outp = bias + PID_action*gain*(PID_propTerm + PID_intTerm + PID_derivTerm)),
				:(PID_AWFactor=1),
				:(PID_action*gain*(intTime*diff(PID_intTerm)-PID_errorI) = PID_output-PID_outp),
				:(PID_outp = bias + PID_action*(gain*PID_propTerm + PID_intTerm + PID_derivTerm)),
				:(PID_AWFactor=1),
				:(PID_action*gain*(intTime*diff(PID_intTerm)-PID_errorI) = PID_output-PID_outp),
				:(PID_outp = bias + PID_action*(gain*(PID_propTerm + PID_intTerm)*(1 + PID_derivTerm))),
				:(PID_AWFactor=1),
			],
			[
				"Input ","Output ","Set Point ","Input first order filter","Input first order filter","setPoint first order filter","setPoint first order filter","Error definition for proportional term","Error definition for derivative term","Error definition for integral term","Error definition for proportional term","Error definition for derivative term","Error definition for integral term","Calculate proportional term","Derivative term filter","Derivative term filter","Calculate derivative term","Scale outp","Calculate clipped output when saturated","Calculate clipped output when not saturated","Calculate unclipped output","","","Calculate integral term with anti-windup","Sum of proportional, integral and derivative terms","Calculate AWFactor","Calculate AWFactor","Calculate integral term with anti-windup","Sum of proportional, integral and derivative terms","Calculate AWFactor","Calculate AWFactor","Calculate integral term with anti-windup","Sum of proportional, integral and derivative terms","Calculate AWFactor","Calculate AWFactor","Calculate integral term","Sum of proportional, integral and derivative terms","Calculate AWFactor - Not in use in this mode","Calculate integral term","Sum of proportional, integral and derivative terms","Calculate AWFactor - Not in use in this mode","Calculate integral term","Sum of proportional, integral and derivative terms","Calculate AWFactor - Not in use in this mode","Calculate integral term with anti-windup and bumpless transfer","Sum of proportional, integral and derivative terms","Calculate AWFactor - Not in use in this mode","Calculate integral term with anti-windup and bumpless transfer","Sum of proportional, integral and derivative terms","Calculate AWFactor - Not in use in this mode","Calculate integral term with anti-windup and bumpless transfer","Sum of proportional, integral and derivative terms","Calculate AWFactor - Not in use in this mode",
			],
			[:PID_Select,:Action,:Mode,:Clip,:alpha,:beta,:bias,:derivTime,:intTime,:gain,:gamma,:tau,:tauSet,:MinInput,:MaxInput,:MinOutput,:MaxOutput,],
			[:Input,:Output,:SetPoint,:PID_derivTerm,:PID_dFilt,:PID_error,:PID_errorD,:PID_errorI,:PID_inputFilt,:PID_intTerm,:PID_outp,:PID_outps,:PID_propTerm,:PID_setPointFilt,:PID_AWFactor,:PID_action,:PID_input,:PID_output,:PID_setPoint,]
		)
	end
	PID_Select::DanaSwitcher
	Action::DanaSwitcher
	Mode::DanaSwitcher
	Clip::DanaSwitcher
	alpha::positive
	beta::positive
	bias::control_signal
	derivTime::time_sec
	intTime::time_sec
	gain::positive
	gamma::positive
	tau::time_sec
	tauSet::time_sec
	MinInput::control_signal
	MaxInput::control_signal
	MinOutput::control_signal
	MaxOutput::control_signal
	Input::control_signal
	Output::control_signal
	SetPoint::control_signal
	PID_derivTerm::control_signal
	PID_dFilt::control_signal
	PID_error::control_signal
	PID_errorD::control_signal
	PID_errorI::control_signal
	PID_inputFilt::control_signal
	PID_intTerm::control_signal
	PID_outp::control_signal
	PID_outps::control_signal
	PID_propTerm::control_signal
	PID_setPointFilt::control_signal
	PID_AWFactor::DanaReal
	PID_action::DanaReal
	PID_input::control_signal
	PID_output::control_signal
	PID_setPoint::control_signal
	initials::Array{Expr,1}
	initialNames::Array{String,1}
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export PID
function initial(in::PID)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
end
function setEquationFlow(in::PID)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	if (tau < 1e-3*"s")
		addEquation(4)
	else
		addEquation(5)
	end
	if (tauSet < 1e-3*"s")
		addEquation(6)
	else
		addEquation(7)
	end
	let switch=Mode
		if switch=="Manual"
			addEquation(8)
			addEquation(9)
			addEquation(10)
		elseif switch=="Automatic"
			addEquation(11)
			addEquation(12)
			addEquation(13)
		end
	end
	addEquation(14)
	if (derivTime < 1e-3*"s")
		addEquation(15)
	else
		addEquation(16)
	end
	addEquation(17)
	addEquation(18)
	let switch=Clip
		if switch=="Clipped"
			if abs(PID_outps)>1
				addEquation(19)
			else
				addEquation(20)
			end
		elseif switch=="Unclipped"
			addEquation(21)
		end
	end
	let switch=Action
		if switch=="Direct"
			addEquation(22)
		elseif switch=="Reverse"
			addEquation(23)
		end
	end
	let switch=PID_Select
		if switch=="Ideal_AW"
			addEquation(24)
			addEquation(25)
			if abs(PID_outps)>1 && (PID_action*sign(PID_outps)*PID_errorI)>0
				addEquation(26)
			else
				addEquation(27)
			end
		elseif switch=="Parallel_AW"
			addEquation(28)
			addEquation(29)
			if abs(PID_outps)>1 && (PID_action*sign(PID_outps)*PID_errorI)>0
				addEquation(30)
			else
				addEquation(31)
			end
		elseif switch=="Series_AW"
			addEquation(32)
			addEquation(33)
			if abs(PID_outps)>1 && (PID_action*sign(PID_outps)*PID_errorI)>0
				addEquation(34)
			else
				addEquation(35)
			end
		elseif switch=="Ideal"
			addEquation(36)
			addEquation(37)
			addEquation(38)
		elseif switch=="Parallel"
			addEquation(39)
			addEquation(40)
			addEquation(41)
		elseif switch=="Series"
			addEquation(42)
			addEquation(43)
			addEquation(44)
		elseif switch=="Ideal_AWBT"
			addEquation(45)
			addEquation(46)
			addEquation(47)
		elseif switch=="Parallel_AWBT"
			addEquation(48)
			addEquation(49)
			addEquation(50)
		elseif switch=="Series_AWBT"
			addEquation(51)
			addEquation(52)
			addEquation(53)
		end
	end
end
function atributes(in::PID,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/PID"
	drive!(fields,_)
	return fields
end
PID(_::Dict{Symbol,Any})=begin
	newModel=PID()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
struct Model{model} end
Model(model::Symbol) = Model{model}()
Model(model::String) = Model{Symbol(model)}()

macro Model_str(s)
  return :(Model{$(Expr(:quote, Symbol(s)))})
end

struct Observable{obs} end
Observable(obs::Symbol) = Observable{obs}()
Observable(obs::String) = Observable{Symbol(obs)}()

macro Observable_str(s)
  return :(Observable{$(Expr(:quote, Symbol(s)))})
end

∫(f, a, b) = quadgk(f, a, b)[1]
ITensorInfiniteMPS.nrange(model::Model) = 2; #required to keep everything compatible with the current implementation for 2 band models

# Create an infinite sum of Hamiltonian terms
function InfiniteSum{T}(model::Model, s::Vector; kwargs...) where {T}
  return InfiniteSum{T}(model, infsiteinds(s); kwargs...)
end

function InfiniteSum{MPO}(model::Model, s::CelledVector; kwargs...)
  N = length(s)
  mpos = [MPO(model, s, n; kwargs...) for n in 1:N] #slightly improved version. Note: the current implementation does not really allow for staggered potentials for example
  return InfiniteSum{MPO}(mpos, translator(s))
end

function InfiniteSum{ITensor}(model::Model, s::CelledVector; kwargs...)
  N = length(s)
  itensors = [ITensor(model, s, n; kwargs...) for n in 1:N] #slightly improved version. Note: the current implementation does not really allow for staggered potentials for example
  return InfiniteSum{ITensor}(itensors, translator(s))
end

# MPO building version
function ITensors.MPO(model::Model, s::CelledVector, n::Int64; kwargs...)
  n1, n2 = 1, 2
  opsum = OpSum(model, n1, n2; kwargs...)
  return MPO(opsum, [s[x] for x in n:(n + nrange(model) - 1)]) #modification to allow for more than two sites per term in the Hamiltonians
end

# Version accepting IndexSet
function ITensors.ITensor(model::Model, s::CelledVector, n::Int64; kwargs...)
  return prod(MPO(model, s, n; kwargs...)) #modification to allow for more than two sites per term in the Hamiltonians
end

# Helper function to make an MPO
import ITensors: op
op(::OpName"Zero", ::SiteType, s::Index) = ITensor(s', dag(s))

function InfiniteMPOMatrix(model::Model, s::CelledVector; kwargs...)
  return InfiniteMPOMatrix(model, s, translator(s); kwargs...)
end

function InfiniteMPOMatrix(model::Model, s::CelledVector, translator::Function; kwargs...)
  N = length(s)
  temp_H = InfiniteSum{MPO}(model, s; kwargs...)
  range_H = nrange(temp_H)[1]
  ls = CelledVector([Index(1, "Link,c=1,n=$n") for n in 1:N], translator)
  mpos = [Matrix{ITensor}(undef, 1, 1) for i in 1:N]
  for j in 1:N
    Hmat = fill(op("Zero", s[j]), range_H + 1, range_H + 1)
    identity = op("Id", s[j])
    Hmat[1, 1] = identity
    Hmat[end, end] = identity
    for n in 0:(range_H - 1)
      idx = findfirst(x -> x == j, findsites(temp_H[j - n]; ncell=N))
      if isnothing(idx)
        Hmat[range_H + 1 - n, range_H - n] = identity
      else
        Hmat[range_H + 1 - n, range_H - n] = temp_H[j - n][idx]#replacetags(linkinds, temp_H[j - n][idx], "Link, l=$n", tags(ls[j-1]))
      end
    end
    mpos[j] = Hmat
    #mpos[j] += dense(Hmat) * setelt(ls[j-1] => total_dim) * setelt(ls[j] => total_dim)
  end
  #return mpos
  return InfiniteMPOMatrix(mpos, translator)
end


text:
export ContinuousSystem,
       DiscreteSystem,
       add_dimension,
       next_set

# name alises
const LCS = LinearContinuousSystem
const LDS = LinearDiscreteSystem
const CLCS = ConstrainedLinearContinuousSystem
const CLDS = ConstrainedLinearDiscreteSystem
const CLCCS = ConstrainedLinearControlContinuousSystem
const CLCDS = ConstrainedLinearControlDiscreteSystem
const ACS = AffineContinuousSystem
const ADS = AffineDiscreteSystem
const CACCS = ConstrainedAffineControlContinuousSystem
const CACDS = ConstrainedAffineControlDiscreteSystem
const CACS = ConstrainedAffineContinuousSystem
const CADS = ConstrainedAffineDiscreteSystem
const BBCS = BlackBoxContinuousSystem
const CBBCS = ConstrainedBlackBoxContinuousSystem
const CBBCCS = ConstrainedBlackBoxControlContinuousSystem

export LCS, LDS, CLCS, CLDS, CLCCS, CLCDS, CACCS, CACDS, CACS, CADS, IVP, BBCS,
       CBBCS, CBBCCS

import Base: *
import LazySets.constrained_dimensions

*(M::AbstractMatrix, input::ConstantInput) =  ConstantInput(M * input.U)

# convenience functions
next_set(inputs::ConstantInput) = collect(nextinput(inputs, 1))[1]
next_set(inputs::AbstractInput, state::Int64) = collect(nextinput(inputs, state))[1]

"""
    add_dimension(A::AbstractMatrix, m=1)

Adds an extra zero row and column to a matrix.

### Input

- `A` -- matrix
- `m` -- (optional, default: `1`) the number of extra dimensions

### Examples

```jldoctest add_dimension_test
julia> A = [0.4 0.25; 0.46 -0.67]
2×2 Array{Float64,2}:
 0.4    0.25
 0.46  -0.67
julia> add_dimension(A)
3×3 Array{Float64,2}:
 0.4    0.25  0.0
 0.46  -0.67  0.0
 0.0    0.0   0.0
```

```jldoctest add_dimension_test
julia> add_dimension(A, 2)
4×4 Array{Float64,2}:
 0.4    0.25  0.0  0.0
 0.46  -0.67  0.0  0.0
 0.0    0.0   0.0  0.0
 0.0    0.0   0.0  0.0
```
"""
function add_dimension(A::AbstractMatrix, m=1)
    n = size(A, 1)
    return vcat(hcat(A, zeros(n, m)), zeros(m, n+m))
end

"""
    add_dimension(X::LazySet, m=1)::LazySet

Adds an extra dimension to a LazySet through a Cartesian product.

### Input

- `X` -- a lazy set
- `m` -- (optional, default: `1`) the number of extra dimensions

### Examples

```jldoctest add_dimension_set
julia> X = BallInf(ones(9), 0.5);

julia> dim(X)
9

julia> Xext = add_dimension(X);

julia> dim(Xext)
10

julia> X = ZeroSet(4);

julia> dim(add_dimension(X))
5

julia> typeof(X)
ZeroSet{Float64}
```

More than one dimension can be added passing the second argument:

```jldoctest add_dimension_set
julia> Xext = add_dimension(BallInf(zeros(10), 0.1), 4);

julia> dim(Xext)
14
```

### Notes

In the special case that the given set is a zero set, instead of cartesian product
a new zero set with extended dimensions is returned.
"""
function add_dimension(X::LazySet, m=1)::LazySet
    return X * ZeroSet(m)
end

function add_dimension(X::ZeroSet, m=1)::ZeroSet
    return ZeroSet(dim(X)+m)
end

"""
    add_dimension(cs, m=1)

Adds an extra dimension to a continuous system.

### Input

- `cs` -- continuous system
- `m` -- (optional, default: `1`) the number of extra dimensions

### Examples

```jldoctest add_dimension_cont_sys
julia> using MathematicalSystems, SparseArrays

julia> A = sprandn(3, 3, 0.5);

julia> X0 = BallInf(zeros(3), 1.0);

julia> s = InitialValueProblem(LinearContinuousSystem(A), X0);

julia> sext = add_dimension(s);

julia> statedim(sext)
4
```

If there is an input set, it is also extended:

```jldoctest add_dimension_cont_sys
julia> U = ConstantInput(Ball2(ones(3), 0.1));

julia> s = InitialValueProblem(ConstrainedLinearControlContinuousSystem(A, Matrix(1.0I, size(A)), nothing, U), X0);

julia> sext = add_dimension(s);

julia> statedim(sext)
4

julia> dim(next_set(inputset(sext)))
4
```

Extending a system with a varying input set:

If there is an input set, it is also extended:

```jldoctest add_dimension_cont_sys
julia> U = VaryingInput([Ball2(ones(3), 0.1 * i) for i in 1:3]);

julia> s = InitialValueProblem(ConstrainedLinearControlContinuousSystem(A, Matrix(1.0I, size(A)), nothing, U), X0);

julia> sext = add_dimension(s);

julia> statedim(sext)
4

julia> dim(next_set(inputset(sext), 1))
4
```

Extending a varing input set with more than one extra dimension:

```jldoctest add_dimension_cont_sys
julia> sext = add_dimension(s, 7);

julia> statedim(sext)
10

julia> dim(next_set(inputset(sext), 1))
10
```
"""
function add_dimension(cs, m=1)
    Aext = add_dimension(cs.s.A, m)
    X0ext = add_dimension(cs.x0, m)
    if hasmethod(inputset, Tuple{typeof(cs.s)})
        Uext = map(x -> add_dimension(x, m), inputset(cs))
        s = ConstrainedLinearControlContinuousSystem(Aext, Matrix(1.0I, size(Aext)), nothing, Uext)
    else
        s = LinearContinuousSystem(Aext)
    end
    return InitialValueProblem(s, X0ext)
end

"""
    constrained_dimensions(HS::HybridSystem)::Dict{Int,Vector{Int}}

For each location, compute all dimensions that are constrained in the invariant
or the guard of any outgoing transition.

### Input

- `HS`  -- hybrid system

### Output

A dictionary mapping the index of each location ``ℓ`` to the dimension indices
that are constrained in ``ℓ``.
"""
function constrained_dimensions(HS::HybridSystem)::Dict{Int,Vector{Int}}
    result = Dict{Int,Vector{Int}}()
    sizehint!(result, nstates(HS))
    for mode in states(HS)
        vars = Vector{Int}()
        append!(vars, constrained_dimensions(stateset(HS, mode)))
        for transition in out_transitions(HS, mode)
            append!(vars, constrained_dimensions(stateset(HS, transition)))
        end
        result[mode] = unique(vars)
    end

    return result
end


text:
struct SavingPathCallback <: DiffusionSpecificCallback
end

#NOTE this needs to be coordinated with DiffusionVis.jl
struct PathPlottingCallback <: DiffusionSpecificCallback
end

struct AuxLawAdaptationCallback <: DiffusionSpecificCallback
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#    Tests of formatters.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "Formatters" begin
    expected = """
┌────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│      1 │  false │      1 │      1 │
│      0 │   true │      0 │      0 │
│      3 │  false │      3 │      3 │
│      0 │   true │      0 │      0 │
│      5 │  false │      5 │      5 │
│      0 │   true │      0 │      0 │
└────────┴────────┴────────┴────────┘
"""

    # Single formatter.
    formatter = (data,i,j) -> begin
        if j != 2
            return isodd(i) ? i : 0
        else
            return data
        end
    end
    result = pretty_table(String, data; formatters = formatter)
    @test result == expected

    # Two formatters.
    expected = """
┌────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│      1 │  false │      1 │      1 │
│     -1 │   true │     -1 │     -1 │
│      3 │  false │      3 │      3 │
│     -1 │   true │     -1 │     -1 │
│      5 │  false │      5 │      5 │
│     -1 │   true │     -1 │     -1 │
└────────┴────────┴────────┴────────┘
"""

    f1 = (data,i,j) -> begin
        if j != 2
            return isodd(i) ? i : 0
        else
            return data
        end
    end

    f2 = (data,i,j) -> begin
        j != 2 && return data == 0 ? -1 : data
        return data
    end

    result = pretty_table(String, data; formatters = (f1, f2))
    @test result == expected
end

@testset "Pre-defined formatters" begin

    # ft_round
    # ==========================================================================
    expected = """
┌────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│    1.0 │    0.0 │    1.0 │    1.0 │
│    2.0 │    1.0 │    2.0 │    2.0 │
│    3.0 │    0.0 │    3.0 │    3.0 │
│    4.0 │    1.0 │    4.0 │    4.0 │
│    5.0 │    0.0 │    5.0 │    5.0 │
│    6.0 │    1.0 │    6.0 │    6.0 │
└────────┴────────┴────────┴────────┘
"""
    result = pretty_table(String, data; formatters = ft_round(1))
    @test result == expected

    expected = """
┌────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│    1.0 │  false │    1.0 │      1 │
│    2.0 │   true │    2.0 │      2 │
│    3.0 │  false │    3.0 │      3 │
│    4.0 │   true │    4.0 │      4 │
│    5.0 │  false │    5.0 │      5 │
│    6.0 │   true │    6.0 │      6 │
└────────┴────────┴────────┴────────┘
"""
    result = pretty_table(String, data; formatters = ft_round(1, [3, 1]))
    @test result == expected

    # Check if `ft_round` correctly avoid unsupported types.

    vec = ["Test", :symbol, 'a', π, exp(1), log(19)]

    expected = """
┌────────┐
│ Col. 1 │
├────────┤
│   Test │
│ symbol │
│      a │
│   3.14 │
│   2.72 │
│   2.94 │
└────────┘
"""

    result = pretty_table(String, vec; formatters = ft_round(2))
    @test result == expected

    expected = """
┌─────────┐
│  Col. 1 │
├─────────┤
│  "Test" │
│ :symbol │
│     'a' │
│    3.14 │
│    2.72 │
│    2.94 │
└─────────┘
"""

    result = pretty_table(String, vec; formatters = ft_round(2), renderer = :show)
    @test result == expected

    # ft_printf
    # ==========================================================================

    expected = """
┌──────────┬──────────┬──────────┬──────────┐
│   Col. 1 │   Col. 2 │   Col. 3 │   Col. 4 │
├──────────┼──────────┼──────────┼──────────┤
│    1.000 │    0.000 │    1.000 │    1.000 │
│    2.000 │    1.000 │    2.000 │    2.000 │
│    3.000 │    0.000 │    3.000 │    3.000 │
│    4.000 │    1.000 │    4.000 │    4.000 │
│    5.000 │    0.000 │    5.000 │    5.000 │
│    6.000 │    1.000 │    6.000 │    6.000 │
└──────────┴──────────┴──────────┴──────────┘
"""
    result = pretty_table(String, data; formatters = ft_printf("%8.3f"))
    @test result == expected

    expected = """
┌──────────┬────────┬────────┬──────────┐
│   Col. 1 │ Col. 2 │ Col. 3 │   Col. 4 │
├──────────┼────────┼────────┼──────────┤
│    1.000 │  false │    1.0 │    1.000 │
│    2.000 │   true │    2.0 │    2.000 │
│    3.000 │  false │    3.0 │    3.000 │
│    4.000 │   true │    4.0 │    4.000 │
│    5.000 │  false │    5.0 │    5.000 │
│    6.000 │   true │    6.0 │    6.000 │
└──────────┴────────┴────────┴──────────┘
"""
    result = pretty_table(String, data; formatters = ft_printf("%8.3f",[1,4]))
    @test result == expected

    expected = """
┌──────────┬────────┬────────┬──────────┐
│   Col. 1 │ Col. 2 │ Col. 3 │   Col. 4 │
├──────────┼────────┼────────┼──────────┤
│     1.00 │  false │    1.0 │   1.0000 │
│     2.00 │   true │    2.0 │   2.0000 │
│     3.00 │  false │    3.0 │   3.0000 │
│     4.00 │   true │    4.0 │   4.0000 │
│     5.00 │  false │    5.0 │   5.0000 │
│     6.00 │   true │    6.0 │   6.0000 │
└──────────┴────────┴────────┴──────────┘
"""
    result = pretty_table(
        String,
        data;
        formatters = ft_printf(["%8.2f", "%8.4f"], [1, 4])
    )
    @test result == expected
end


text:
"""
    FEMMHeatDiffModule

Module for operations on interiors of domains to construct system matrices and
system vectors for linear heat conduction/diffusion.
"""
module FEMMHeatDiffModule

using FinEtools.FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import FinEtools.FENodeSetModule: FENodeSet
import FinEtools.FESetModule: FESet, nodesperelem, manifdim, gradN!
import FinEtools.MatHeatDiffModule: MatHeatDiff
import FinEtools.IntegDomainModule: IntegDomain, integrationdata, Jacobianvolume
import FinEtools.CSysModule: CSys, updatecsmat!
import FinEtools.FieldModule: ndofs, gatherdofnums!, gatherfixedvalues_asvec!, gathervalues_asvec!, gathervalues_asmat!
import FinEtools.NodalFieldModule: NodalField
import FinEtools.ElementalFieldModule: ElementalField
import FinEtools.AssemblyModule: SysvecAssemblerBase, SysmatAssemblerBase, SysmatAssemblerSparseSymm, startassembly!, assemble!, makematrix!, makevector!, SysvecAssembler
import FinEtools.ForceIntensityModule: ForceIntensity
import FinEtools.FEMMBaseModule: FEMMAbstractBase, inspectintegpoints
import FinEtools.MatrixUtilityModule: add_gkgt_ut_only!, complete_lt!, locjac!
import LinearAlgebra: mul!, Transpose
At_mul_B!(C, A, B) = mul!(C, Transpose(A), B)
A_mul_B!(C, A, B) = mul!(C, A, B)
import LinearAlgebra: norm, dot

# Type for heat diffusion finite element modeling machine.
mutable struct FEMMHeatDiff{S<:FESet, F<:Function, M<:MatHeatDiff} <: FEMMAbstractBase
    integdomain::IntegDomain{S, F} # geometry data
    mcsys::CSys # updater of the material orientation matrix
    material::M # material object
end

"""
    FEMMHeatDiff(integdomain::IntegDomain{S, F}, material::M) where {S<:FESet, F<:Function, M<:MatHeatDiff}

Construct with the default orientation matrix (identity).
"""
function FEMMHeatDiff(integdomain::IntegDomain{S, F}, material::M) where {S<:FESet, F<:Function, M<:MatHeatDiff}
    return FEMMHeatDiff(integdomain, CSys(manifdim(integdomain.fes)), material)
end

function  buffers1(self::FEMMHeatDiff, geom::NodalField{FFlt}, temp::NodalField{FFlt})
    # Constants
    fes = self.integdomain.fes
    nfes = count(fes); # number of finite elements in the set
    ndn = ndofs(temp); # number of degrees of freedom per node
    nne = nodesperelem(fes); # number of nodes for element
    sdim = ndofs(geom);   # number of space dimensions
    mdim = manifdim(fes); # manifold dimension of the element
    Kedim = ndn*nne;      # dimension of the element matrix
    elmat = fill(zero(FFlt), Kedim, Kedim); # buffer
    elvec = fill(zero(FFlt), Kedim); # buffer
    elvecfix = fill(zero(FFlt), Kedim); # buffer
    dofnums = fill(zero(FInt), Kedim); # buffer
    loc = fill(zero(FFlt), 1, sdim); # buffer
    J = fill(zero(FFlt), sdim, mdim); # buffer
    RmTJ = fill(zero(FFlt), mdim, mdim); # buffer
    gradN = fill(zero(FFlt), nne, mdim); # buffer
    kappa_bar = fill(zero(FFlt), mdim, mdim); # buffer
    kappa_bargradNT = fill(zero(FFlt), mdim, nne); # buffer
    return dofnums, loc, J, RmTJ, gradN, kappa_bar, kappa_bargradNT, elmat, elvec, elvecfix
end

"""
    conductivity(self::FEMMHeatDiff,
      assembler::A, geom::NodalField{FFlt},
      temp::NodalField{FFlt}) where {A<:SysmatAssemblerBase}

Compute the conductivity matrix.
"""
function conductivity(self::FEMMHeatDiff,  assembler::A, geom::NodalField{FFlt},  temp::NodalField{FFlt}) where {A<:SysmatAssemblerBase}
    fes = self.integdomain.fes
    npts,  Ns,  gradNparams,  w,  pc = integrationdata(self.integdomain);
    # Prepare assembler and buffers
    dofnums, loc, J, RmTJ, gradN, kappa_bar, kappa_bargradNT, elmat = buffers1(self, geom, temp)
    # Thermal conductivity matrix is in local  material coordinates.
    kappa_bar = self.material.tangentmoduli!(self.material, kappa_bar)
    startassembly!(assembler, size(elmat,1), size(elmat,2), count(fes), temp.nfreedofs, temp.nfreedofs);
    for i = 1:count(fes) # Loop over elements
        fill!(elmat,  0.0); # Initialize element matrix
        for j=1:npts # Loop over quadrature points
            locjac!(loc, J, geom.values, fes.conn[i], Ns[j], gradNparams[j])
            Jac = Jacobianvolume(self.integdomain, J, loc, fes.conn[i], Ns[j]);
            updatecsmat!(self.mcsys, loc, J, fes.label[i]);
            At_mul_B!(RmTJ,  self.mcsys.csmat,  J); # local Jacobian matrix
            gradN!(fes, gradN, gradNparams[j], RmTJ);
            # Add the product gradN*kappa_bar*gradNT*(Jac*w[j])
            add_gkgt_ut_only!(elmat, gradN, (Jac*w[j]), kappa_bar, kappa_bargradNT)
        end # Loop over quadrature points
        complete_lt!(elmat)
        gatherdofnums!(temp, 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 conductivity(self::FEMMHeatDiff, geom::NodalField{FFlt},  temp::NodalField{FFlt})
    assembler = SysmatAssemblerSparseSymm();
    return conductivity(self, assembler, geom, temp);
end

"""
    nzebcloadsconductivity(self::FEMMHeatDiff,
      assembler::A,  geom::NodalField{FFlt},
      temp::NodalField{FFlt}) where {A<:SysvecAssemblerBase}

Compute load vector for nonzero EBC of prescribed temperature.
"""
function nzebcloadsconductivity(self::FEMMHeatDiff, assembler::A,  geom::NodalField{FFlt},  temp::NodalField{FFlt}) where {A<:SysvecAssemblerBase}
    fes = self.integdomain.fes
    npts,  Ns,  gradNparams,  w,  pc = integrationdata(self.integdomain);
    # Prepare assembler and buffers
    dofnums, loc, J, RmTJ, gradN, kappa_bar, kappa_bargradNT, elmat, elvec, elvecfix = buffers1(self, geom, temp)
    # Thermal conductivity matrix is in local  material coordinates.
    kappa_bar = self.material.tangentmoduli!(self.material, kappa_bar)
    startassembly!(assembler,  temp.nfreedofs);
    # Now loop over all finite elements in the set
    for i = 1:count(fes) # Loop over elements
        gatherfixedvalues_asvec!(temp, elvecfix, fes.conn[i]);# retrieve element coordinates
        if norm(elvecfix, Inf) != 0. # Is the load nonzero?
            fill!(elmat,  0.0);
            for j=1:npts # Loop over quadrature points
                locjac!(loc, J, geom.values, fes.conn[i], Ns[j], gradNparams[j])
                Jac = Jacobianvolume(self.integdomain, J, loc, fes.conn[i], Ns[j]);
                updatecsmat!(self.mcsys, loc, J, fes.label[i]);
                At_mul_B!(RmTJ,  self.mcsys.csmat,  J); # local Jacobian matrix
                gradN!(fes, gradN, gradNparams[j], RmTJ);
                # Add the product gradN*kappa_bar*gradNT*(Jac*w[j])
                add_gkgt_ut_only!(elmat, gradN, (Jac*w[j]), kappa_bar, kappa_bargradNT)
            end # Loop over quadrature points
            complete_lt!(elmat)
            A_mul_B!(elvec, elmat, elvecfix) # compute  the load vector
            gatherdofnums!(temp, dofnums, fes.conn[i]); # retrieve degrees of freedom
            assemble!(assembler,  -elvec,  dofnums); # assemble element load vector
        end
    end
    return makevector!(assembler);
end

function nzebcloadsconductivity(self::FEMMHeatDiff,  geom::NodalField{FFlt},   temp::NodalField{FFlt})
    assembler = SysvecAssembler()
    return  nzebcloadsconductivity(self, assembler, geom, temp);
end

"""
    energy(self::FEMMHeatDiff, geom::NodalField{FFlt},  temp::NodalField{FFlt})

Compute the "energy" integral over the interior domain.
"""
function energy(self::FEMMHeatDiff, geom::NodalField{FFlt},  temp::NodalField{FFlt})
    fes = self.integdomain.fes
    npts,  Ns,  gradNparams,  w,  pc = integrationdata(self.integdomain);
    # Prepare assembler and buffers
    dofnums, loc, J, RmTJ, gradN, kappa_bar, kappa_bargradNT, elmat, elvec, elvecfix = buffers1(self, geom, temp)
    # Thermal conductivity matrix is in local  material coordinates.
    kappa_bar = self.material.tangentmoduli!(self.material, kappa_bar)
    gradT = fill(0.0, 1, size(gradN, 2))
    fluxT = deepcopy(gradT)
    energy = 0.0
    # Now loop over all finite elements in the set
    for i = 1:count(fes) # Loop over elements
        gathervalues_asvec!(temp, elvec, fes.conn[i]);# retrieve element coordinates
        for j=1:npts # Loop over quadrature points
            locjac!(loc, J, geom.values, fes.conn[i], Ns[j], gradNparams[j])
            Jac = Jacobianvolume(self.integdomain, J, loc, fes.conn[i], Ns[j]);
            updatecsmat!(self.mcsys, loc, J, fes.label[i]);
            At_mul_B!(RmTJ,  self.mcsys.csmat,  J); # local Jacobian matrix
            gradN!(fes, gradN, gradNparams[j], RmTJ);
            At_mul_B!(gradT, elvec, gradN)
            A_mul_B!(fluxT, gradT, kappa_bar)
            energy += dot(vec(gradT), vec(fluxT)) * (Jac*w[j])
        end # Loop over quadrature points
    end
    return energy;
end

"""
    inspectintegpoints(self::FEMMHeatDiff, geom::NodalField{FFlt}, u::NodalField{T}, temp::NodalField{FFlt}, felist::FIntVec, inspector::F, idat, quantity=:heatflux; context...) where {F<:Function}

Inspect integration point quantities.

# Inputs
- `geom` - reference geometry field
- `u` - displacement field (ignored)
- `temp` - temperature 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::FEMMHeatDiff, geom::NodalField{FFlt}, u::NodalField{T}, temp::NodalField{FFlt}, felist::FIntVec, inspector::F, idat, quantity=:heatflux; context...) where {T<:Number, F<:Function}
    fes = self.integdomain.fes
    npts,  Ns,  gradNparams,  w,  pc = integrationdata(self.integdomain);
    dofnums, loc, J, RmTJ, gradN, kappa_bar, kappa_bargradNT, elmat, elvec, elvecfix = buffers1(self, geom, temp)
    # Thermal conductivity matrix is in local  material coordinates.
    kappa_bar =  self.material.thermal_conductivity;
    # Sort out  the output requirements
    outputcsys = self.mcsys; # default: report the vector quantities in the material coord system
    for apair in pairs(context)
        sy, val = apair
        if sy == :outputcsys
            outputcsys = val
        end
    end
    t= 0.0
    dt = 0.0
    Te = fill(zero(FFlt), nodesperelem(fes)) # nodal temperatures -- buffer
    nne = nodesperelem(fes); # number of nodes for element
    sdim = ndofs(geom);            # number of space dimensions
    xe = fill(zero(FFlt), nne, sdim); # array of node coordinates -- buffer
    qpgradT = fill(zero(FFlt), 1, sdim); # Temperature gradient -- buffer
    qpflux = fill(zero(FFlt), sdim); # thermal strain -- buffer
    out1 = fill(zero(FFlt), sdim); # output -- buffer
    out =  fill(zero(FFlt), sdim);# output -- buffer
    # 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!(geom, xe, fes.conn[i]);# retrieve element coords
        gathervalues_asvec!(temp, Te, fes.conn[i]);# retrieve element temperatures
        for j = 1:npts # Loop over quadrature points
            locjac!(loc, J, geom.values, fes.conn[i], Ns[j], gradNparams[j])
            Jac = Jacobianvolume(self.integdomain, J, loc, fes.conn[i], Ns[j]);
            updatecsmat!(self.mcsys, loc, J, fes.label[i]);
            At_mul_B!(RmTJ,  self.mcsys.csmat,  J); # local Jacobian matrix
            gradN!(fes, gradN, gradNparams[j], RmTJ);
            # Quadrature point quantities
            A_mul_B!(qpgradT, reshape(Te, 1, :), gradN); # temperature gradient in material coordinates
            out = self.material.update!(self.material, qpflux, out, vec(qpgradT), 0.0, 0.0, loc, fes.label[i], quantity)
            if (quantity == :heatflux)   # Transform heat flux vector,  if that is "out"
                A_mul_B!(out1, transpose(self.mcsys.csmat), out);# To global coord sys
                A_mul_B!(out, outputcsys.csmat, out1);# To output coord sys
            end
            # Call the inspector
            idat = inspector(idat, i, fes.conn[i], xe, out, loc);
        end # Loop over quadrature points
    end # Loop over elements
    return idat; # return the updated inspector data
end

end


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# McCormick.jl
# A McCormick operator library in Julia
# See https://github.com/PSORLab/McCormick.jl
#############################################################################
# src/forward_operators/multiplication.jl
# Defines multiplication operator.
#############################################################################

# Differentiable multiplication kernel definition
@inline function sigu(y::Float64)
	x = y/MC_DIFF_MU1T
	(0.0 <= x) ? x^(1.0/MC_DIFF_MU) : -abs(x)^(1.0/MC_DIFF_MU)
end
@inline abs_pu(x::Float64) = abs(x)^MC_DIFF_MU1T
@inline psi_pow(x::Float64) = MC_DIFF_MU1T*x*abs(x)^MC_DIFF_MU1N

@inline function gCxAcv(alpha::Interval{Float64}, beta::Interval{Float64},
	            lambda::Interval{Float64}, nu::Interval{Float64}, x1::MC, x2::MC)

	# gCxA pre-terms
	alplo = alpha.lo
	alphi = alpha.hi
	betlo = beta.lo
	bethi = beta.hi
	LmdDel = lambda.hi - lambda.lo
	NuDel = nu.hi - nu.lo
	LmdSum = lambda.lo + lambda.hi
	NuSum = nu.lo + nu.hi

	xslo = lambda.lo + LmdDel*(((nu.hi - betlo)/NuDel) - sigu(NuSum/NuDel))
	xshi = lambda.lo + LmdDel*(((nu.hi - bethi)/NuDel) - sigu(NuSum/NuDel))
	yslo = nu.lo + NuDel*(((lambda.hi - alplo)/LmdDel) - sigu(LmdSum/LmdDel))
	yshi = nu.lo + NuDel*(((lambda.hi - alphi)/LmdDel) - sigu(LmdSum/LmdDel))

	term1, _ = mid3(alplo, alphi, xslo)
	term2, _ = mid3(alplo, alphi, xshi)
	term3, _ = mid3(betlo, bethi, yslo)
	term4, _ = mid3(betlo, bethi, yshi)

	# calculates the convex relaxation
	tempGxA1a = LmdDel*NuDel*abs_pu((betlo - nu.lo)/NuDel - (lambda.hi - term1)/LmdDel)
	tempGxA2a = LmdDel*NuDel*abs_pu((bethi - nu.lo)/NuDel - (lambda.hi - term2)/LmdDel)
	tempGxA3a = LmdDel*NuDel*abs_pu((term3 - nu.lo)/NuDel - (lambda.hi - alplo)/LmdDel)
	tempGxA4a = LmdDel*NuDel*abs_pu((term4 - nu.lo)/NuDel - (lambda.hi - alphi)/LmdDel)

	NuDotLmd = lambda.lo*nu.lo + lambda.hi*nu.hi
	tempGxA1 = 0.5*(term1*NuSum + betlo*LmdSum - NuDotLmd + tempGxA1a)
	tempGxA2 = 0.5*(term2*NuSum + bethi*LmdSum - NuDotLmd + tempGxA2a)
	tempGxA3 = 0.5*(alplo*NuSum + term3*LmdSum - NuDotLmd + tempGxA3a)
	tempGxA4 = 0.5*(alphi*NuSum + term4*LmdSum - NuDotLmd + tempGxA4a)

	# gets minima which is cv/cc/Intv depending on arguments
	cv = min(tempGxA1, tempGxA2, tempGxA3, tempGxA4)

	if (cv == tempGxA1)
		psi_eval = psi_pow((betlo - nu.lo)/NuDel - (lambda.hi - term1)/LmdDel)
		psi_mx1 = 0.5*(NuSum + NuDel*psi_eval)
		psi_my1 = 0.5*(LmdSum + LmdDel*psi_eval)
	elseif (cv == tempGxA2)
		psi_eval = psi_pow((bethi - nu.lo)/NuDel - (lambda.hi - term2)/LmdDel)
		psi_mx1 = 0.5*(NuSum + NuDel*psi_eval)
		psi_my1 = 0.5*(LmdSum + LmdDel*psi_eval)
	elseif (cv == tempGxA3)
		psi_eval = psi_pow((term3 - nu.lo)/NuDel - (lambda.hi - alplo)/LmdDel)
		psi_mx1 = 0.5*(NuSum + NuDel*psi_eval)
		psi_my1 = 0.5*(LmdSum + LmdDel*psi_eval)
	else
		psi_eval = psi_pow((term4 - nu.lo)/NuDel - (lambda.hi - alphi)/LmdDel)
		psi_mx1 = 0.5*(NuSum + NuDel*psi_eval)
		psi_my1 = 0.5*(LmdSum + LmdDel*psi_eval)
	end

    cv_grad = max(0.0, psi_mx1)*x1.cv_grad + min(0.0, psi_mx1)*x1.cc_grad +
		      max(0.0, psi_my1)*x2.cv_grad + min(0.0, psi_my1)*x2.cc_grad

    return cv, cv_grad
end

@inline function gCxAcc(alpha::Interval{Float64},beta::Interval{Float64},lambda::Interval{Float64},
	            nu::Interval{Float64},x1::MC,x2::MC)

	# gCxA pre-terms
	alplo = alpha.lo
	alphi = alpha.hi
	betlo = beta.lo
	bethi = beta.hi
	LmdDel = lambda.hi - lambda.lo
	NuDel = nu.hi - nu.lo
	LmdSum = lambda.lo + lambda.hi
	NuSum = nu.lo + nu.hi

	xslo = lambda.lo + LmdDel*(((nu.hi - betlo)/NuDel) - sigu(NuSum/NuDel))
	xshi = lambda.lo + LmdDel*(((nu.hi - bethi)/NuDel) - sigu(NuSum/NuDel))
	yslo = nu.lo + NuDel*(((lambda.hi - alplo)/LmdDel) - sigu(LmdSum/LmdDel))
	yshi = nu.lo + NuDel*(((lambda.hi - alphi)/LmdDel) - sigu(LmdSum/LmdDel))

	# calculates term 1
	term1, _ = mid3(alplo, alphi, xslo)
	term2, _ = mid3(alplo, alphi, xshi)
	term3, _ = mid3(betlo, bethi, yslo)
	term4, _ = mid3(betlo, bethi, yshi)

	# calculates the convex relaxation
	tempGxA1a = LmdDel*NuDel*abs_pu((betlo - nu.lo)/NuDel - (lambda.hi - term1)/LmdDel)
	tempGxA2a = LmdDel*NuDel*abs_pu((bethi - nu.lo)/NuDel - (lambda.hi - term2)/LmdDel)
	tempGxA3a = LmdDel*NuDel*abs_pu((term3 - nu.lo)/NuDel - (lambda.hi - alplo)/LmdDel)
	tempGxA4a = LmdDel*NuDel*abs_pu((term4 - nu.lo)/NuDel - (lambda.hi - alphi)/LmdDel)

	NuDotLmd = lambda.lo*nu.lo+lambda.hi*nu.hi
	tempGxA1 = 0.5*(term1*NuSum + betlo*LmdSum - NuDotLmd + tempGxA1a)
	tempGxA2 = 0.5*(term2*NuSum + bethi*LmdSum - NuDotLmd + tempGxA2a)
	tempGxA3 = 0.5*(alplo*NuSum + term3*LmdSum - NuDotLmd + tempGxA3a)
	tempGxA4 = 0.5*(alphi*NuSum + term4*LmdSum - NuDotLmd + tempGxA4a)

	# gets minima which is cv/cc/Intv depending on arguments
	ncc = min(tempGxA1, tempGxA2, tempGxA3, tempGxA4)

   if (ncc == tempGxA1)
	   psi_eval = psi_pow((betlo - nu.lo)/NuDel - (lambda.hi - term1)/LmdDel)
	   psi_mx1 = 0.5*(NuSum + NuDel*psi_eval)
	   psi_my1 = 0.5*(LmdSum + LmdDel*psi_eval)
	elseif (ncc == tempGxA2)
		psi_eval = psi_pow((bethi - nu.lo)/NuDel - (lambda.hi - term2)/LmdDel)
		psi_mx1 = 0.5*(NuSum + NuDel*psi_eval)
		psi_my1 = 0.5*(LmdSum + LmdDel*psi_eval)
	elseif (ncc == tempGxA3)
		psi_eval = psi_pow((term3 - nu.lo)/NuDel - (lambda.hi - alplo)/LmdDel)
		psi_mx1 = 0.5*(NuSum + NuDel*psi_eval)
		psi_my1 = 0.5*(LmdSum + LmdDel*psi_eval)
	else
		psi_eval = psi_pow((term4 - nu.lo)/NuDel - (lambda.hi - alphi)/LmdDel)
		psi_mx1 = 0.5*(NuSum + NuDel*psi_eval)
		psi_my1 = 0.5*(LmdSum + LmdDel*psi_eval)
	end
	ncc_grad = max(0.0,psi_mx1)*x1.cc_grad + min(0.0,psi_mx1)*x1.cv_grad -
		       max(0.0,psi_my1)*x2.cc_grad - min(0.0,psi_my1)*x2.cv_grad

    return ncc, ncc_grad
end

@inline function gCxAIntv(alpha::Interval{Float64},beta::Interval{Float64},lambda::Interval{Float64},
	              nu::Interval{Float64},x1::MC,x2::MC)

	# gCxA pre-terms
	alplo = alpha.lo
	alphi = alpha.hi
	betlo = beta.lo
	bethi = beta.hi
	LmdDel = lambda.hi-lambda.lo
	NuDel = nu.hi-nu.lo
	LmdSum = lambda.lo+lambda.hi
	NuSum = nu.lo+nu.hi

	xslo = lambda.lo + LmdDel*(((nu.hi - betlo)/NuDel) - sigu(NuSum/NuDel))
	xshi = lambda.lo + LmdDel*(((nu.hi - bethi)/NuDel) - sigu(NuSum/NuDel))
	yslo = nu.lo + NuDel*(((lambda.hi - alplo)/LmdDel) - sigu(LmdSum/LmdDel))
	yshi = nu.lo + NuDel*(((lambda.hi - alphi)/LmdDel) - sigu(LmdSum/LmdDel))

	# calculates terms
	term1, _ = mid3(alplo, alphi, xslo)
	term2, _ = mid3(alplo, alphi, xshi)
	term3, _ = mid3(betlo, bethi, yslo)
	term4, _ = mid3(betlo, bethi, yshi)

	# calculates the convex relaxation
	tempGxA1a = LmdDel*NuDel*abs_pu((betlo - nu.lo)/NuDel - (lambda.hi - term1)/LmdDel)
	tempGxA2a = LmdDel*NuDel*abs_pu((bethi - nu.lo)/NuDel - (lambda.hi - term2)/LmdDel)
	tempGxA3a = LmdDel*NuDel*abs_pu((term3 - nu.lo)/NuDel - (lambda.hi - alplo)/LmdDel)
	tempGxA4a = LmdDel*NuDel*abs_pu((term4 - nu.lo)/NuDel - (lambda.hi - alphi)/LmdDel)

	NuLmd = lambda.lo*nu.lo+lambda.hi*nu.hi
	tempGxA1 = 0.5*(term1*NuSum + betlo*LmdSum - NuLmd + tempGxA1a)
	tempGxA2 = 0.5*(term2*NuSum + bethi*LmdSum - NuLmd + tempGxA2a)
	tempGxA3 = 0.5*(alplo*NuSum + term3*LmdSum - NuLmd + tempGxA3a)
	tempGxA4 = 0.5*(alphi*NuSum + term4*LmdSum - NuLmd + tempGxA4a)

	# gets minima which is cv/cc/Intv depending on arguments
	return min(tempGxA1, tempGxA2, tempGxA3, tempGxA4)
end

@inline function multiply_MV(x1::MC{N,Diff}, x2::MC{N,Diff}, z::Interval{Float64}) where {N}

	x1cv = x1.cv
	x2cv = x2.cv
	x1cc = x1.cc
	x2cc = x2.cc
	x1lo = x1.Intv.lo
	x2lo = x2.Intv.lo
	x1hi = x1.Intv.hi
	x2hi = x2.Intv.hi
	x1h_l = x1hi - x1lo
	x1h_v = x1hi - x1cv
	x2h_l = x2hi - x2lo
	x2v_l = x2cv - x2lo
	cv = 0.0
	cc = 0.0

	alpha0 = Interval{Float64}(x1cv, x1cc)
	beta0 = Interval{Float64}(x2cv, x2cc)

	if (0.0 <= x1lo) && (0.0 <= x2lo)

		btemp1 = max(0.0, x2v_l/x2h_l - x1h_v/x1h_l)
		btemp2 = x2lo + MC_DIFF_MU1T*x2h_l*btemp1^MC_DIFF_MU
		btemp3 = x1lo + MC_DIFF_MU1T*x1h_l*btemp1^MC_DIFF_MU

		cv = x1cv*x2lo + x1lo*x2v_l + x1h_l*x2h_l*btemp1^MC_DIFF_MU1T
		cv_grad = btemp2*x1.cv_grad + btemp3*x2.cv_grad

	elseif (x1hi <= 0.0) && (x2hi <= 0.0)

		btemp1 = max(0.0,(x2hi - x2cc)/x2h_l - (x1cc-x1lo)/x1h_l)
		btemp2 = -x2hi + MC_DIFF_MU1T*(x2h_l)*max(0.0, btemp1)^MC_DIFF_MU
		btemp3 = -x1hi + MC_DIFF_MU1T*(x1h_l)*max(0.0, btemp1)^MC_DIFF_MU

		cv = x1cc*x2hi + x1hi*(x2cc - x2hi) + (x1h_l)*(x2h_l)*btemp1^MC_DIFF_MU1T
		cv_grad = (-btemp2)*x1.cc_grad - btemp3*x2.cc_grad

	else
		cv,cv_grad = gCxAcv(alpha0, beta0, x1.Intv, x2.Intv, x1, x2)
	end

	if (x1hi <= 0.0) && (0.0 <= x2lo)

		btemp1 = max(0.0, x2v_l/x2h_l - (x1cc-x1lo)/x1h_l)
		btemp2 = x2lo + MC_DIFF_MU1T*x2h_l*btemp1^MC_DIFF_MU
		btemp3 = -x1hi + MC_DIFF_MU1T*x1h_l*btemp1^MC_DIFF_MU
		cc = x1cc*x2lo - x1hi*(x2lo - x2cv) - x1h_l*x2h_l*btemp1^MC_DIFF_MU1T
		cc_grad = btemp2*x1.cc_grad - btemp3*x2.cv_grad

	elseif (0.0 <= x1.Intv.lo) && (x2.Intv.hi <= 0.0)

		btemp1 = max(0.0, (x2hi-x2cc)/x2h_l - x1h_v/x1h_l)
		btemp2 = -x2hi + MC_DIFF_MU1T*(x2h_l)*btemp1^MC_DIFF_MU
		btemp3 = x1lo + MC_DIFF_MU1T*x1h_l*btemp1^MC_DIFF_MU
		cc = x1cv*x2hi + x1lo*(x2cc - x2hi) - x1h_l*x2h_l*btemp1^MC_DIFF_MU1T
		cc_grad = -btemp2*x1.cv_grad + btemp3*x2.cc_grad

	else
		cct::Float64, cc_gradt = gCxAcc(-alpha0, beta0, -x1.Intv, x2.Intv, x1, x2)
		cc = -cct
		cc_grad = cc_gradt
	end

	cnst = x2.cnst ? x1.cnst : (x1.cnst ? x2.cnst : x1.cnst || x2.cnst)

	return MC{N,Diff}(cv, cc, z, cv_grad, cc_grad, cnst)
end

@inline function mult_kernel(x1::MC{N,Diff}, x2::MC{N,Diff}, y::Interval{Float64}) where N
	degen1 = (x1.Intv.hi - x1.Intv.lo) == 0.0
	degen2 = (x2.Intv.hi - x2.Intv.lo) == 0.0
    if degen1 && !degen2
		out = x1.Intv.hi*x2
	elseif !degen1 && degen2
		out = x2.Intv.hi*x1
	elseif degen1 && degen2
		out = MC{N,Diff}(x1.Intv.hi*x2.Intv.hi)
    else
    	out = multiply_MV(x1, x2, y)
		check_relaxation_error!(out)
	end
	return out
end

# Nonsmooth multiplication kernel definition
@inline function mul1_u1pos_u2pos(x1::MC{N,T}, x2::MC{N,T}, z::Interval{Float64}, cnst::Bool) where {N,T<:RelaxTag}
    xLc = z.lo
    xUc = z.hi
  	cv1 = x2.Intv.hi*x1.cv + x1.Intv.hi*x2.cv - x1.Intv.hi*x2.Intv.hi
  	cv2 = x2.Intv.lo*x1.cv + x1.Intv.lo*x2.cv - x1.Intv.lo*x2.Intv.lo
  	if (cv1 > cv2)
    	cv = cv1
    	cv_grad = x2.Intv.hi*x1.cv_grad
  	else
    	cv = cv2
    	cv_grad = x2.Intv.lo*x1.cv_grad
  	end

  	cc1 = x2.Intv.lo*x1.cc + x1.Intv.hi*x2.cc - x1.Intv.hi*x2.Intv.lo
  	cc2 = x2.Intv.hi*x1.cc + x1.Intv.lo*x2.cc - x1.Intv.lo*x2.Intv.hi
  	if (cc1 < cc2)
    	cc = cc1
    	cc_grad = x2.Intv.lo*x1.cc_grad
  	else
    	cc = cc2
    	cc_grad = x2.Intv.hi*x1.cc_grad
  	end
  	cv, cc, cv_grad, cc_grad = cut(xLc, xUc, cv, cc, cv_grad, cc_grad)
  	return MC{N,T}(cv, cc, z, cv_grad, cc_grad, cnst)
end
@inline function mul1_u1pos_u2mix(x1::MC{N,T}, x2::MC{N,T}, z::Interval{Float64}, cnst::Bool) where {N,T<:RelaxTag}
    xLc = z.lo
    xUc = z.hi
    cv1 = x2.Intv.hi*x1.cv + x1.Intv.hi*x2.cv - x1.Intv.hi*x2.Intv.hi
    cv2 = x2.Intv.lo*x1.cc + x1.Intv.lo*x2.cv - x1.Intv.lo*x2.Intv.lo
    if (cv1 > cv2)
        cv = cv1
        cv_grad = x2.Intv.hi*x1.cv_grad
    else
        cv = cv2
        cv_grad = x2.Intv.lo*x1.cc_grad
    end

    cc1 = x2.Intv.lo*x1.cv + x1.Intv.hi*x2.cc - x1.Intv.hi*x2.Intv.lo
    cc2 = x2.Intv.hi*x1.cc + x1.Intv.lo*x2.cc - x1.Intv.lo*x2.Intv.hi
    if (cc1 < cc2)
        cc = cc1
        cc_grad = x2.Intv.lo*x1.cv_grad
    else
        cc = cc2
        cc_grad = x2.Intv.hi*x1.cc_grad
    end
    cv,cc,cv_grad,cc_grad = cut(xLc,xUc,cv,cc,cv_grad,cc_grad)
    return MC{N,T}(cv, cc, z, cv_grad, cc_grad, cnst)
end
@inline function mul1_u1mix_u2mix(x1::MC{N,T}, x2::MC{N,T}, z::Interval{Float64}, cnst::Bool) where {N,T<:RelaxTag}
    xLc = z.lo
    xUc = z.hi
    cv1 = x2.Intv.hi*x1.cv + x1.Intv.hi*x2.cv - x1.Intv.hi*x2.Intv.hi
    cv2 = x2.Intv.lo*x1.cc + x1.Intv.lo*x2.cc - x1.Intv.lo*x2.Intv.lo
    cv = cv1
    if (cv1 > cv2)
        cv = cv1
        cv_grad = x2.Intv.hi*x1.cv_grad
    else
        cv = cv2
        cv_grad = x2.Intv.lo*x1.cc_grad
    end
    cc1 = x2.Intv.lo*x1.cv + x1.Intv.hi*x2.cc - x1.Intv.hi*x2.Intv.lo
    cc2 = x2.Intv.hi*x1.cc + x1.Intv.lo*x2.cv - x1.Intv.lo*x2.Intv.hi
    cc = cc1
    if (cc1 < cc2)
        cc = cc1
        cc_grad = x2.Intv.lo*x1.cv_grad
    else
        cc = cc2
        cc_grad = x2.Intv.hi*x1.cc_grad
    end
    cv,cc,cv_grad,cc_grad = cut(xLc,xUc,cv,cc,cv_grad,cc_grad)
    return MC{N,T}(cv, cc, z, cv_grad, cc_grad, cnst)
end
@inline function mul2_u1pos_u2pos(x1::MC{N,T}, x2::MC{N,T}, z::Interval{Float64}) where {N,T<:RelaxTag}
	xLc = z.lo
	xUc = z.hi
	cnst = x2.cnst ? x1.cnst : (x1.cnst ? x2.cnst : x1.cnst || x2.cnst)
	cv1 = x2.Intv.hi*x1.cv + x1.Intv.hi*x2.cv - x1.Intv.hi*x2.Intv.hi
	cv2 = x2.Intv.lo*x1.cv + x1.Intv.lo*x2.cv - x1.Intv.lo*x2.Intv.lo
	if (cv1 > cv2)
		cv = cv1
		cv_grad = x2.Intv.hi*x1.cv_grad + x1.Intv.hi*x2.cv_grad
	else
		cv = cv2
		cv_grad = x2.Intv.lo*x1.cv_grad + x1.Intv.lo*x2.cv_grad
	end
	cc1 = x2.Intv.lo*x1.cc + x1.Intv.hi*x2.cc - x1.Intv.hi*x2.Intv.lo
	cc2 = x2.Intv.hi*x1.cc + x1.Intv.lo*x2.cc - x1.Intv.lo*x2.Intv.hi
	if (cc1 < cc2)
		cc = cc1
		cc_grad = x2.Intv.lo*x1.cc_grad + x1.Intv.hi*x2.cc_grad
	else
		cc = cc2
		cc_grad = x2.Intv.hi*x1.cc_grad + x1.Intv.lo*x2.cc_grad
	end
	cv, cc, cv_grad, cc_grad = cut(xLc, xUc, cv, cc, cv_grad, cc_grad)
	return MC{N,T}(cv, cc, z, cv_grad, cc_grad, cnst)
end
@inline function mul2_u1pos_u2mix(x1::MC{N,T}, x2::MC{N,T}, z::Interval{Float64}, cnst::Bool) where {N,T<:RelaxTag}
    xLc = z.lo
    xUc = z.hi
    cv1 = x2.Intv.hi*x1.cv + x1.Intv.hi*x2.cv - x1.Intv.hi*x2.Intv.hi
    cv2 = x2.Intv.lo*x1.cc + x1.Intv.lo*x2.cv - x1.Intv.lo*x2.Intv.lo
    if (cv1 > cv2)
        cv = cv1
        cv_grad = x1.Intv.hi*x2.cv_grad
    else
        cv = cv2
        cv_grad = x1.Intv.lo*x2.cc_grad
    end
    cc1 = x2.Intv.lo*x1.cv + x1.Intv.hi*x2.cc - x1.Intv.hi*x2.Intv.lo
    cc2 = x2.Intv.hi*x1.cc + x1.Intv.lo*x2.cc - x1.Intv.lo*x2.Intv.hi
    if (cc1 < cc2)
        cc = cc1
        cc_grad = x1.Intv.hi*x2.cc_grad
    else
        cc = cc2
        cc_grad = x1.Intv.lo*x2.cc_grad
    end
    cv, cc, cv_grad, cc_grad = cut(xLc, xUc, cv, cc, cv_grad, cc_grad)
    return MC{N,T}(cv, cc, z, cv_grad, cc_grad, cnst)
end
@inline function mul2_u1mix_u2mix(x1::MC{N,T}, x2::MC{N,T}, z::Interval{Float64}) where {N,T<:RelaxTag}
	xLc = z.lo
	xUc = z.hi
  	cnst = x2.cnst ? x1.cnst : (x1.cnst ? x2.cnst : x1.cnst || x2.cnst)
	cv1 = x2.Intv.hi*x1.cv + x1.Intv.hi*x2.cv - x1.Intv.hi*x2.Intv.hi
	cv2 = x2.Intv.lo*x1.cc + x1.Intv.lo*x2.cc - x1.Intv.lo*x2.Intv.lo
	if (cv1 > cv2)
		cv = cv1
		cv_grad = x2.Intv.hi*x1.cv_grad + x1.Intv.hi*x2.cv_grad
	else
		cv = cv2
		cv_grad = x2.Intv.lo*x1.cc_grad + x1.Intv.lo*x2.cc_grad
	end
	cc1 = x2.Intv.lo*x1.cv + x1.Intv.hi*x2.cc - x1.Intv.hi*x2.Intv.lo
	cc2 = x2.Intv.hi*x1.cc + x1.Intv.lo*x2.cv - x1.Intv.lo*x2.Intv.hi
	if (cc1 < cc2)
		cc = cc1
		cc_grad = x2.Intv.lo*x1.cv_grad + x1.Intv.hi*x2.cc_grad
	else
		cc = cc2
		cc_grad = x2.Intv.hi*x1.cc_grad + x1.Intv.lo*x2.cv_grad
	end
	cv, cc, cv_grad, cc_grad = cut(xLc, xUc, cv, cc, cv_grad, cc_grad)
	return MC{N,T}(cv, cc, z, cv_grad, cc_grad, cnst)
end
@inline function mul3_u1pos_u2mix(x1::MC{N,T}, x2::MC{N,T}, z::Interval{Float64}) where {N,T<:RelaxTag}
	xLc = z.lo
	xUc = z.hi
    cnst = x2.cnst ? x1.cnst : (x1.cnst ? x2.cnst : x1.cnst || x2.cnst)
	cv1 = x2.Intv.hi*x1.cv + x1.Intv.hi*x2.cv - x1.Intv.hi*x2.Intv.hi
	cv2 = x2.Intv.lo*x1.cc + x1.Intv.lo*x2.cv - x1.Intv.lo*x2.Intv.lo
	if cv1 > cv2
		cv = cv1
		cv_grad = x2.Intv.hi*x1.cv_grad + x1.Intv.hi*x2.cv_grad
	else
		cv = cv2
		cv_grad = x2.Intv.lo*x1.cc_grad + x1.Intv.lo*x2.cv_grad
	end
	cc1 = x2.Intv.lo*x1.cv + x1.Intv.hi*x2.cc - x1.Intv.hi*x2.Intv.lo
	cc2 = x2.Intv.hi*x1.cc + x1.Intv.lo*x2.cc - x1.Intv.lo*x2.Intv.hi
	if cc1 < cc2
		cc = cc1
		cc_grad = x2.Intv.lo*x1.cv_grad + x1.Intv.hi*x2.cc_grad
	else
		cc = cc2
		cc_grad = x2.Intv.hi*x1.cc_grad + x1.Intv.lo*x2.cc_grad
	end
	cv, cc, cv_grad, cc_grad = cut(xLc, xUc, cv, cc, cv_grad, cc_grad)
	return MC{N,T}(cv, cc, z, cv_grad, cc_grad, cnst)
end

@inline function multiply_STD_NS(x1::MC, x2::MC, y::Interval{Float64})
	if x2.Intv.lo >= 0.0
    	x2.cnst && (return mul1_u1pos_u2pos(x1, x2, y, x1.cnst))
    	x1.cnst && (return mul1_u1pos_u2pos(x2, x1, y, x2.cnst))
  		return mul2_u1pos_u2pos(x1, x2, y)
	elseif x2.Intv.hi <= 0.0
		return -mult_kernel(x1, -x2, -y)
	else
    	x2.cnst && (return mul1_u1pos_u2mix(x1, x2, y, x1.cnst))
    	x1.cnst && (return mul2_u1pos_u2mix(x1, x2, y, x2.cnst))
	end
	mul3_u1pos_u2mix(x1, x2, y)
end

@inline function mult_kernel(x1::MC{N,NS}, x2::MC{N,NS}, y::Interval{Float64}) where N
	isone(x1) && (return x2)
	isone(x2) && (return x1)
	if x1.Intv.lo >= 0.0
		return multiply_STD_NS(x1, x2, y)
	elseif x1.Intv.hi <= 0.0
		(x2.Intv.lo >= 0.0) && (return -mult_kernel(-x1, x2, -y))
	    (x2.Intv.hi <= 0.0) && (return mult_kernel(-x1, -x2, y))
		return -mult_kernel(x2, -x1, -y)
	elseif x2.Intv.lo >= 0.0
		return mult_kernel(x2, x1, y)
	elseif x2.Intv.hi <= 0.0
		return -mult_kernel(-x2, x1, -y)
	else
    	x2.cnst && (return mul1_u1mix_u2mix(x1, x2, y, x1.cnst))
    	x1.cnst && (return mul1_u1mix_u2mix(x2, x1, y, x2.cnst))
	end
	return mul2_u1mix_u2mix(x1, x2, y)
end

# Nonsmooth multivariate multiplication kernel definition
@inline mul_MV_ns1cv(x1::Float64, x2::Float64, MC1::MC, MC2::MC) = MC2.Intv.hi*x1 + MC1.Intv.hi*x2 - MC2.Intv.hi*MC1.Intv.hi
@inline mul_MV_ns2cv(x1::Float64, x2::Float64, MC1::MC, MC2::MC) = MC2.Intv.lo*x1 + MC1.Intv.lo*x2 - MC2.Intv.lo*MC1.Intv.lo
@inline mul_MV_ns3cv(x1::Float64, x2::Float64, MC1::MC, MC2::MC) = max(mul_MV_ns1cv(x1,x2,MC1,MC2), mul_MV_ns2cv(x1,x2,MC1,MC2))

@inline mul_MV_ns1cc(x1::Float64, x2::Float64, MC1::MC, MC2::MC) = MC2.Intv.lo*x1 + MC1.Intv.hi*x2 - MC2.Intv.lo*MC1.Intv.hi
@inline mul_MV_ns2cc(x1::Float64, x2::Float64, MC1::MC, MC2::MC) = MC2.Intv.hi*x1 + MC1.Intv.lo*x2 - MC2.Intv.hi*MC1.Intv.lo
@inline mul_MV_ns3cc(x1::Float64, x2::Float64, MC1::MC, MC2::MC) = min(mul_MV_ns1cc(x1,x2,MC1,MC2), mul_MV_ns2cc(x1,x2,MC1,MC2))

@inline tol_MC(x::Float64, y::Float64) = abs(x-y) <= (MC_MV_TOL + MC_MV_TOL*0.5*abs(x+y))

@inline function multiply_MV_NS(x1::MC{N,MV}, x2::MC{N,MV}, zIntv::Interval{Float64}, cnst::Bool) where N

	flag = true

	# convex calculation
    k = -diam(x2)/diam(x1)
    z = (x1.Intv.hi*x2.Intv.hi - x1.Intv.lo*x2.Intv.lo)/diam(x1)

 	x1vta = mid3v(x1.cv, x1.cc, (x2.cv - z)/k)
 	x1vtb = mid3v(x1.cv, x1.cc, (x2.cc - z)/k)
 	x2vta = mid3v(x2.cv, x2.cc, k*x1.cv + z)
 	x2vtb = mid3v(x2.cv ,x2.cc, k*x1.cc + z)

    x1vt = (x1.cv, x1.cc, x1vta, x1vtb, x1.cv, x1.cc)
    x2vt = (x2vta, x2vtb, x2.cv, x2.cc, x2.cv, x2.cc)

	vt = mul_MV_ns3cv.(x1vt, x2vt, x1, x2)
    cv, cvi = findmin(vt)

 	if tol_MC(mul_MV_ns1cv(x1vt[cvi], x2vt[cvi], x1, x2), mul_MV_ns2cv(x1vt[cvi], x2vt[cvi], x1, x2))

		alph1 = 0.0
		alph2 = 1.0

		MC1thin = tol_MC(x1.cv, x1.cc)
 		if ~MC1thin && (x1vt[cvi] > x1.cv) && ~tol_MC(x1vt[cvi], x1.cv)
 			alph2 = min(alph2, -x2.Intv.lo/diam(x2))
 		end

 		if ~MC1thin && (x1vt[cvi] < x1.cc) && ~tol_MC(x1vt[cvi], x1.cc)
			alph1 = max(alph1, -x2.Intv.lo/diam(x2))
		end

		MC2thin = tol_MC(x2.cv, x2.cc)
		if ~MC2thin && (x2vt[cvi] > x2.cv) && ~tol_MC(x2vt[cvi], x2.cv)
			alph2 = min(alph2, -x1.Intv.lo/diam(x1))
		end
		if ~MC2thin && (x2vt[cvi] < x2.cc) && ~tol_MC(x2vt[cvi], x2.cc)
			alph1 = max(alph1, -x1.Intv.lo/diam(x1))
		end

		alphthin = tol_MC(alph1, alph2)
		(~alphthin && (alph1 > alph2)) && (return false, x1)
		myalph = (alph1 + alph2)/2.0
	elseif mul_MV_ns1cv(x1vt[cvi], x2vt[cvi], x1,x2) > mul_MV_ns2cv(x1vt[cvi], x2vt[cvi], x1, x2)
		myalph = 1.0
	else
		myalph = 0.0
	end

	sigma1 = x2.Intv.lo + myalph*diam(x2)
	sigma2 = x1.Intv.lo + myalph*diam(x1)
	cv_grad1 = (sigma1 > 0.0) ? x1.cv_grad : x1.cc_grad
	cv_grad2 = (sigma2 > 0.0) ? x2.cv_grad : x2.cc_grad

	x1.cnst && (cv_grad = sigma2*cv_grad2)
	x2.cnst && (cv_grad = sigma1*cv_grad1)
	(~x1.cnst && ~x1.cnst) && (cv_grad = sigma1*cv_grad1 + sigma2*cv_grad2)
	(x1.cnst && x1.cnst) && (cv_grad = zero(SVector{N,Float64}))

	 # concave calculation
	 k = diam(x2)/diam(x1)
	 z = (x1.Intv.hi*x2.Intv.lo - x1.Intv.lo*x2.Intv.hi)/diam(x1)

	 x1vta = mid3v(x1.cv, x1.cc, (x2.cv - z)/k)
	 x1vtb = mid3v(x1.cv, x1.cc, (x2.cc - z)/k)
	 x2vta = mid3v(x2.cv, x2.cc, k*x1.cv + z)
	 x2vtb = mid3v(x2.cv, x2.cc, k*x1.cc + z)

	 x1vt = (x1.cv, x1.cc, x1vta, x1vtb, x1.cv, x1.cc)
	 x2vt = (x2vta, x2vtb, x2.cv, x2.cc, x2.cc, x2.cv)
	 vt = mul_MV_ns3cc.(x1vt, x2vt, x1, x2)
	 cc, cci = findmax(vt)

	 if tol_MC(mul_MV_ns1cc(x1vt[cci], x2vt[cci], x1, x2), mul_MV_ns2cc(x1vt[cci], x2vt[cci], x1, x2))

		alph1 = 0.0
	 	alph2 = 1.0

		 MC1thin = tol_MC(x1.cv, x1.cc)
		 if ~MC1thin && (x1vt[cci] > x1.cv) && ~tol_MC(x1vt[cci], x1.cv)
		 	alph1 = max(alph1, -x2.Intv.lo/diam(x2))
		 end
		 if ~MC1thin && (x1vt[cci] < x1.cc) && ~tol_MC(x1vt[cci], x1.cc)
		 	alph2 = min(alph2, -x2.Intv.lo/diam(x2))
		 end

		 MC2thin = tol_MC(x2.cv, x2.cc)
		 if ~MC2thin && (x2vt[cci] > x2.cv) && ~tol_MC(x2vt[cci], x2.cv)
		 	alph2 = min(alph2, x1.Intv.hi/diam(x1))
		 end
		 if ~MC2thin && (x2vt[cci] < x2.cc) && ~tol_MC(x2vt[cci], x2.cc)
		 	alph1 = max(alph1, x1.Intv.hi/diam(x1))
		 end

        alphthin = tol_MC(alph1, alph2)
		(~alphthin && (alph1 > alph2)) && (return false, x1)
	 	myalph = (alph1 + alph2)/2.0
	elseif mul_MV_ns1cc(x1vt[cci], x2vt[cci], x1, x2) > mul_MV_ns2cc(x1vt[cci], x2vt[cci], x1, x2)
	 	myalph = 1.0
	else
		myalph = 0.0
	end

	sigma1 = x2.Intv.lo + myalph*diam(x2)
	sigma2 = x1.Intv.hi - myalph*diam(x1)
	cc_grad1 = (sigma1 > 0.0) ? x1.cc_grad :  x1.cv_grad
	cc_grad2 = (sigma2 > 0.0) ? x2.cc_grad :  x2.cv_grad

	x1.cnst && (cc_grad = sigma2*cc_grad2)
	x2.cnst && (cc_grad = sigma1*cc_grad1)
	(~x1.cnst && ~x1.cnst) && (cc_grad = sigma1*cc_grad1 + sigma2*cc_grad2)
	(x1.cnst && x1.cnst) && (cc_grad = zero(SVector{N,Float64}))

    return flag, MC{N,MV}(cv, cc, zIntv, cv_grad, cc_grad, cnst)
end

@inline function mult_kernel(x1::MC{N,MV}, x2::MC{N,MV}, y::Interval{Float64}) where N
	degen1 = (x1.Intv.hi - x1.Intv.lo) <= MC_DEGEN_TOL
	degen2 = (x2.Intv.hi - x2.Intv.lo) <= MC_DEGEN_TOL
	if degen1 || degen2
		multiply_STD_NS(x1, x2, y)
	end
	flag, x3 = multiply_MV_NS(x1, x2, y, x1.cnst && x2.cnst)
	~flag && (return nan(MC{N,MV}))
	return x3
end

#
@inline function *(x1::MC{N,T}, x2::MC{N,T}) where {N, T<:Union{NS,MV}}
	mult_kernel(x1, x2, x1.Intv*x2.Intv)
end
@inline function *(x1::MC{N,Diff}, x2::MC{N,Diff}) where N
	degen1 = ((x1.Intv.hi - x1.Intv.lo) == 0.0)
	degen2 = ((x2.Intv.hi - x2.Intv.lo) == 0.0)
	if ~(degen1 || degen2)
		if (min(x1.Intv.lo, x2.Intv.lo) < 0.0 < max(x1.Intv.hi, x2.Intv.hi))
			lo_Intv_calc::Float64 = gCxAIntv(x1.Intv, x2.Intv, x1.Intv, x2.Intv, x1, x2)
			hi_Intv_calc::Float64 = -gCxAIntv(-x1.Intv, x2.Intv, -x1.Intv, x2.Intv, x1, x2)
			z = Interval{Float64}(lo_Intv_calc, hi_Intv_calc)
		else
			z = x1.Intv*x2.Intv
		end
	else
		z = x1.Intv*x2.Intv
	end
	return mult_kernel(x1, x2, z)
end

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
@testset "thin_subpixel_edges" begin
    @info "Test: thin_subpixel_edges"

    @testset "API" begin
        img_ref = Gray{N0f8}.(testimage("lake_gray.tif"))
        img = copy(img_ref)

        f = SubpixelNonmaximaSuppression()
        g₁, g₂ = imgradients(img, KernelFactors.scharr)
        mag = hypot.(g₁, g₂)

        g₁_cpy = copy(g₁)
        g₂_cpy = copy(g₂)
        mag_cpy = copy(mag)

        thinned_edges₀ = zeros(eltype(mag), axes(mag))
        offsets₀ = zeros(SVector{2,Float64}, axes(mag))
        thin_subpixel_edges!(thinned_edges₀, offsets₀, mag, g₁, g₂, f)
        thinned_edges₁, offsets₁ = thin_subpixel_edges(mag, g₁, g₂, f)

        @test mag == mag_cpy # mag unchanged
        @test g₁ == g₁_cpy # g₁ unchanged
        @test g₂ == g₂_cpy # g₂ unchanged

        @test all(map(x-> .!all(isnan.(x)), offsets₀))
        @test all(map(x-> .!all(isnan.(x)), offsets₁))

        @test all(offsets₀ .== offsets₁)
        @test thinned_edges₀ == thinned_edges₁
        @test eltype(thinned_edges₀) == eltype(mag)
        @test eltype(thinned_edges₁) == eltype(mag)

        for T in [Float16, Float32, Float64]
            thinned_edges₂ = zeros(T, axes(mag))
            offsets₂ = zeros(SVector{2,T}, axes(mag))
            thin_subpixel_edges!(thinned_edges₂, offsets₂, mag, g₁, g₂, f)
            @test eltype(thinned_edges₂) == T
            @test eltype(offsets₂) == SVector{2,T}

            thinned_edges₃, offsets₃ = thin_subpixel_edges(T, SVector{2, T}, mag, g₁, g₂, f)
            @test eltype(thinned_edges₃) == T
            @test eltype(offsets₃) == SVector{2, T}
        end
    end

    @testset "Numerical" begin
        img = Gray{N0f8}.(load("algorithms/References/circle.png"))
        nms_ref = Gray{N0f8}.(load("algorithms/References/circle_nms.png"))

        f = SubpixelNonmaximaSuppression()

        g₁, g₂ = imgradients(img, KernelFactors.scharr)
        mag = hypot.(g₁, g₂)

        thinned_edges₀, subpixel_offsets₀ = thin_subpixel_edges(mag, g₁, g₂, f)

        for i in CartesianIndices(nms_ref)
            nms_ref[i] == 0.0 ? (@test thinned_edges₀[i] == 0) : (@test thinned_edges₀[i] == mag[i])
            thinned_edges₀[i] == 0.0 ? (@test subpixel_offsets₀[i] == SVector(0.0, 0.0)) : (@test subpixel_offsets₀[i] != SVector(0.0, 0.0))
        end
    end
end


text:
using (LinearAlgebra)
include("mpii.annotation.reader.jl")
include("variant.data.jl")

#=
function headSize = util_get_head_size(rect)
SC_BIAS = 0.6; % 0.8*0.75
headSize = SC_BIAS*norm([rect.x2 rect.y2] - [rect.x1 rect.y1]);
end
=#

function get_head_size(corner_one, corner_two)
    # TODO: This can be parameterized
    SC_BIAS = 0.6
    head_size = SC_BIAS * norm(corner_one .- corner_two)
    return head_size
end

#=
function dist = getDistPCKh(pred,gt,refDist)

assert(size(pred,1) == size(gt,1) && size(pred,2) == size(gt,2) && size(pred,3) == size(gt,3));
assert(size(refDist,1) == size(gt,3));

dist = nan(1,size(pred,2),size(pred,3));

for imgidx = 1:size(pred,3)

    % distance to gt joints
    dist(1,:,imgidx) = sqrt(sum((pred(:,:,imgidx) - gt(:,:,imgidx)).^2,1))./refDist(imgidx);

end
=#

function get_distance_in_PCKh(pred, gt, head_size)
    pckh = sqrt(sum((pred .- gt) .^ 2)) / head_size
    return pckh
end

# TODO: impl PCKH
# NOW it should not be that hard

# Data Items needs to be scaled beforehand
function create_acc_dict()
    person_typed_h_ranged_joint_by_joint_accuracy_dict = Dict()
    person_types = ["single", "multi", "mixed"]
    h_ranges = [hr for hr = 0:0.1:1]
    person_typed_h_ranged_joint_by_joint_accuracy_dict = Dict()

    for person_type in person_types
        h_ranged_joint_by_joint_accuracy_dict = Dict()
        for h in h_ranges
            joint_by_joint_accuracy_dict = Dict()
            # (accuracate_preds / total_number_of_gt)
            for j = 1:global_num_joints
                joint_by_joint_accuracy_dict[j] = (0, 0)
            end
            h_ranged_joint_by_joint_accuracy_dict[h] = joint_by_joint_accuracy_dict
        end
        person_typed_h_ranged_joint_by_joint_accuracy_dict[person_type] =
            h_ranged_joint_by_joint_accuracy_dict
    end
    return person_typed_h_ranged_joint_by_joint_accuracy_dict
end

function original_pckh(
    output,
    data_items;
    h_range = PCKh_range,
    consider_threshold = output_consider_threshold,
    offmat = nothing,
    person_typed_h_ranged_joint_by_joint_accuracy_dict = Dict(),
)

    num_elements = size(output)[4]
    num_joints = global_num_joints
    correct_counter = 0
    non_existant_joint_counter = 0

    for idx = 1:num_elements
        data_item =  isa(data_items, DataItem) ? data_items : data_items[idx]
        visible_joint_ids = data_item.joints[:, 1]

        head_corner_one = (data_item.annorect["x1"], data_item.annorect["y1"])
        head_corner_two = (data_item.annorect["x2"], data_item.annorect["y2"])
        head_size = get_head_size(head_corner_one, head_corner_two)

        for j = 1:num_joints
            joint_sc = output[:, :, j, idx]
            joint_sc[joint_sc.<consider_threshold] .= 0
            max_coord = argmax(joint_sc)

            if !(j in visible_joint_ids)
                non_existant_joint_counter += 1
            else
                data_joint_idx = findlast(x -> x[1] == j, data_item.joints)[1]
                joint_gold = data_item.joints[data_joint_idx, :][2:3]

                pred = (max_coord[2], max_coord[1])

                pred_f8 = pred .* preprocess_stride .+ 0.5 * preprocess_stride

                if offmat != nothing
                    dx = offmat[:, :, 2*j-1, idx][max_coord] * global_locref_stdev
                    dy = offmat[:, :, 2*j, idx][max_coord] * global_locref_stdev
                    pred_f8 = (pred_f8[1] + dx, pred_f8[2] + dy)
                end

                dist_PCKh = get_distance_in_PCKh(pred_f8, joint_gold, head_size)

                if dist_PCKh <= h_range
                    correct_counter += 1
                end

                person_types = ["single", "multi", "mixed"]
                h_ranges = [hr for hr = 0:0.1:1]
                for person_type in person_types
                    if (data_item.is_single && person_type == "single") ||
                       (!data_item.is_single && person_type == "multi") ||
                       (person_type == "mixed")

                        h_ranged_joint_by_joint_accuracy_dict =
                            person_typed_h_ranged_joint_by_joint_accuracy_dict[person_type]
                        for h in h_ranges
                            current_values_correct =
                                h_ranged_joint_by_joint_accuracy_dict[h][j][1]
                            current_values_total =
                                h_ranged_joint_by_joint_accuracy_dict[h][j][2]
                            current_values_total += 1
                            if dist_PCKh <= h
                                current_values_correct += 1
                            end
                            h_ranged_joint_by_joint_accuracy_dict[h][j] =
                                (current_values_correct, current_values_total)
                        end
                    end
                end
            end
        end
    end

    acc =
        100.0 * correct_counter / ((num_elements * num_joints) - non_existant_joint_counter)
    return (acc, person_typed_h_ranged_joint_by_joint_accuracy_dict)
end

function modelized_PCKh_sigm(
    model,
    data,
    data_items;
    h_range = PCKh_range,
    consider_threshold = output_consider_threshold,
)
    results = []
    acc_dist_dict = create_acc_dict()

    if isa(data, VariantData)
        prev_shuffle_value = data.shuffle
        data.shuffle = false
        for (i, (x, ygold)) in enumerate(data)
            ygold = Array{Float32}(ygold)
            output = Array{Float32}(model(x))
            part_detection_scores = output[:, :, 1:global_num_joints, :]
            part_detection_probs = Array{Float32}(sigm.(part_detection_scores))
            channel_dim = size(output)[end-1]
            offmat = nothing
            if channel_dim >= 3 * global_num_joints
                offmat = output[:, :, global_num_joints+1:global_num_joints*3, :]
            end
            pck, acc_dist_dict = original_pckh(
                part_detection_probs,
                data_items[i];
                h_range = h_range,
                consider_threshold = consider_threshold,
                offmat = offmat,
                person_typed_h_ranged_joint_by_joint_accuracy_dict = acc_dist_dict,
            )
            push!(results, pck)
        end
        data.shuffle = prev_shuffle_value
    else # means KnetData
        original_x_data = Knet.atype()(reshape(data.x, data.xsize))
        batchsize = data.batchsize
        total_element_count = size(original_x_data)[4]

        batch_coefficient = floor(total_element_count / batchsize) |> Int
        total_element_count = batch_coefficient * batchsize

        for i = 1:batchsize:total_element_count
            datapart = original_x_data[:, :, :, i:i+batchsize-1]
            dataitem_part = data_items[i:i+batchsize-1]

            output = model(datapart)
            part_detection_scores = output[:, :, 1:global_num_joints, :]
            part_detection_probs = Array{Float32}(sigm.(part_detection_scores))

            channel_dim = size(output)[end-1]
            offmat = nothing
            if channel_dim >= 3 * global_num_joints
                offmat = output[:, :, global_num_joints+1:global_num_joints*3, :]
            end

            pck, acc_dist_dict = original_pckh(
                part_detection_probs,
                dataitem_part;
                h_range = h_range,
                consider_threshold = consider_threshold,
                offmat = offmat,
                person_typed_h_ranged_joint_by_joint_accuracy_dict = acc_dist_dict,
            )

            push!(results, pck)
        end
    end

    final = sum(results) / length(results)
    final /= 100
    return (final, acc_dist_dict)
end

# Naive Accuracy

# https://www.programmersought.com/article/7644537351/#:~:text=1.,PCK%20%2D%20Percentage%20of%20Correct%20Keypoints&text=calculates%20the%20percentage%20of%20detections,used%20as%20a%20normalized%20reference.

function naive_pck(
    output,
    y_gold;
    consider_threshold = output_consider_threshold,
    wrong_scs_cache = nothing,
)
    num_elements = size(y_gold)[4]
    num_joints = global_num_joints
    correct_counter = 0
    for idx = 1:num_elements
        for j = 1:num_joints
            joint_sc = output[:, :, j, idx]
            joint_sc_initial = deepcopy(joint_sc)
            joint_sc[joint_sc.<consider_threshold] .= 0
            max_coord = argmax(joint_sc)
            gt_sc = y_gold[:, :, j, idx]
            if gt_sc[max_coord] == 1
                correct_counter += 1
            elseif joint_sc[max_coord] == 0 && maximum(gt_sc) == 0
                correct_counter += 1
            else
                if wrong_scs_cache != nothing
                    push!(wrong_scs_cache[1], gt_sc)
                    push!(wrong_scs_cache[2], joint_sc)
                    push!(wrong_scs_cache[3], joint_sc_initial)
                end
            end

        end
    end
    acc = 100.0 * correct_counter / (num_elements * num_joints)
    return acc
end

function modelized_naive_pck_sigm(model, data)
    results = []
    for (x, ygold) in data
        ygold = Array{Float32}(ygold)
        ygold = ygold[:, :, 1:global_num_joints, :]

        output = model(x)
        output = output[:, :, 1:global_num_joints, :]

        output = Array{Float32}(sigm.(output))
        pck = naive_pck(output, ygold)
        push!(results, pck)
    end
    final = sum(results) / length(results)
    return final / 100
end

function modelized_naive_pck_sigm_with_wrong_cache(model, data, wrong_scs_cache)
    results = []
    for (x, ygold) in data
        ygold = Array{Float32}(ygold)
        output = Array{Float32}(sigm.(model(x)))
        pck = naive_pck(output, ygold; wrong_scs_cache = wrong_scs_cache)
        push!(results, pck)
    end
    final = sum(results) / length(results)
    return final / 100
end


text:
function errorcheck(error::Integer)
    error == 0 && return nothing
    str = " "^Int(LJM.MAX_NAME_SIZE)
    LJM.ErrorToString(error, str)
    error > 0 && @warn unsafe_string(pointer(str))
    error < 0 && throw(unsafe_string(pointer(str)))
    return nothing
end


text:
# An instruction builder represents a point within a basic block and is the exclusive means
# of building instructions using the C interface.

export Builder,
       position!,
       debuglocation, debuglocation!

@checked struct Builder
    ref::API.LLVMBuilderRef
end
reftype(::Type{Builder}) = API.LLVMBuilderRef

Builder() = Builder(API.LLVMCreateBuilder())
Builder(ctx::Context) = Builder(API.LLVMCreateBuilderInContext(ref(ctx)))

dispose(builder::Builder) = API.LLVMDisposeBuilder(ref(builder))

function Builder(f::Core.Function, args...)
    builder = Builder(args...)
    try
        f(builder)
    finally
        dispose(builder)
    end
end

Base.position(builder::Builder) = BasicBlock(API.LLVMGetInsertBlock(ref(builder)))
position!(builder::Builder, inst::Instruction) =
    API.LLVMPositionBuilderBefore(ref(builder), ref(inst))
position!(builder::Builder, bb::BasicBlock) =
    API.LLVMPositionBuilderAtEnd(ref(builder), blockref(bb))
position!(builder::Builder) = API.LLVMClearInsertionPosition(ref(builder))

Base.insert!(builder::Builder, inst::Instruction) =
    API.LLVMInsertIntoBuilder(ref(builder), ref(inst))
Base.insert!(builder::Builder, inst::Instruction, name::String) =
    API.LLVMInsertIntoBuilderWithName(ref(builder), ref(inst), name)

debuglocation(builder::Builder) =
    MetadataAsValue(API.LLVMGetCurrentDebugLocation(ref(builder)))
debuglocation!(builder::Builder) =
    API.LLVMSetCurrentDebugLocation(ref(builder), convert(API.LLVMValueRef, C_NULL))
debuglocation!(builder::Builder, loc::MetadataAsValue) =
    API.LLVMSetCurrentDebugLocation(ref(builder), ref(loc))
debuglocation!(builder::Builder, inst::Instruction) =
    API.LLVMSetInstDebugLocation(ref(builder), ref(inst))

## build methods

# TODO/IDEAS:
# - dynamic dispatch based on `llvmtype` (eg. disambiguating `add!` and `fadd!`)
# - auto `ref(...)` using `@generated deref` for API calling
#   (cfr. `cconvert` for `ccall`)

# NOTE: the return values for these operations are, according to the C API, always a Value.
#       however, the C++ API learns us that we can be more strict.

export ret!, br!, switch!, indirectbr!, invoke!, resume!, unreachable!,

       binop!, add!, nswadd!, nuwadd!, fadd!, sub!, nswsub!, nuwsub!, fsub!, mul!, nswmul!,
       nuwmul!, fmul!, udiv!, sdiv!, exactsdiv!, fdiv!, urem!, srem!, frem!, neg!, nswneg!,
       nuwneg!, fneg!,

       shl!, lshr!, ashr!, and!, or!, xor!, not!,

       extract_element!, insert_element!, shuffle_vector!,

       extract_value!, insert_value!,

       alloca!, array_alloca!, malloc!, array_malloc!, free!, load!, store!, fence!,
       atomic_rmw!, atomic_cmpxchg!, gep!, inbounds_gep!, struct_gep!,

       trunc!, zext!, sext!, fptoui!, fptosi!, uitofp!, sitofp!, fptrunc!, fpext!,
       ptrtoint!, inttoptr!, bitcast!, addrspacecast!, zextorbitcast!, sextorbitcast!,
       truncorbitcast!, cast!, pointercast!, intcast!, fpcast!,

       icmp!, fcmp!, phi!, select!, call!, va_arg!, landingpad!,

       globalstring!, globalstring_ptr!, isnull!, isnotnull!, ptrdiff!

# terminator instructions

ret!(builder::Builder) =
    Instruction(API.LLVMBuildRetVoid(ref(builder)))

ret!(builder::Builder, V::Value) =
    Instruction(API.LLVMBuildRet(ref(builder), ref(V)))

ret!(builder::Builder, RetVals::Vector{T}) where {T<:Value} =
    Instruction(API.LLVMBuildAggregateRet(ref(builder), ref.(RetVals), length(RetVals)))

br!(builder::Builder, Dest::BasicBlock) =
    Instruction(API.LLVMBuildBr(ref(builder), blockref(Dest)))

br!(builder::Builder, If::Value, Then::BasicBlock, Else::BasicBlock) =
    Instruction(API.LLVMBuildCondBr(ref(builder), ref(If), blockref(Then), blockref(Else)))

switch!(builder::Builder, V::Value, Else::BasicBlock, NumCases::Integer=10) =
    Instruction(API.LLVMBuildSwitch(ref(builder), ref(V), blockref(Else), NumCases))

indirectbr!(builder::Builder, Addr::Value, NumDests::Integer=10) =
    Instruction(API.LLVMBuildIndirectBr(ref(builder), ref(Addr), NumDests))

invoke!(builder::Builder, Fn::Value, Args::Vector{T}, Then::BasicBlock, Catch::BasicBlock, Name::String="") where {T<:Value} =
    Instruction(API.LLVMBuildInvoke(ref(builder), ref(Fn), ref.(Args), length(Args), blockref(Then), blockref(Catch), Name))

resume!(builder::Builder, Exn::Value) =
    Instruction(API.LLVMBuildResume(ref(builder), ref(Exn)))

unreachable!(builder::Builder) =
    Instruction(API.LLVMBuildUnreachable(ref(builder)))

# binary operations

binop!(builder::Builder, Op::API.LLVMOpcode, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildBinOp(ref(builder), Op, ref(LHS), ref(RHS), Name))

add!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildAdd(ref(builder), ref(LHS), ref(RHS), Name))

nswadd!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildNSWAdd(ref(builder), ref(LHS), ref(RHS), Name))

nuwadd!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildNUWAdd(ref(builder), ref(LHS), ref(RHS), Name))

fadd!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildFAdd(ref(builder), ref(LHS), ref(RHS), Name))

sub!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildSub(ref(builder), ref(LHS), ref(RHS), Name))

nswsub!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildNSWSub(ref(builder), ref(LHS), ref(RHS), Name))

nuwsub!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildNUWSub(ref(builder), ref(LHS), ref(RHS), Name))

fsub!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildFSub(ref(builder), ref(LHS), ref(RHS), Name))

mul!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildMul(ref(builder), ref(LHS), ref(RHS), Name))

nswmul!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildNSWMul(ref(builder), ref(LHS), ref(RHS), Name))

nuwmul!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildNUWMul(ref(builder), ref(LHS), ref(RHS), Name))

fmul!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildFMul(ref(builder), ref(LHS), ref(RHS), Name))

udiv!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildUDiv(ref(builder), ref(LHS), ref(RHS), Name))

sdiv!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildSDiv(ref(builder), ref(LHS), ref(RHS), Name))

exactsdiv!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildExactSDiv(ref(builder), ref(LHS), ref(RHS), Name))

fdiv!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildFDiv(ref(builder), ref(LHS), ref(RHS), Name))

urem!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildURem(ref(builder), ref(LHS), ref(RHS), Name))

srem!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildSRem(ref(builder), ref(LHS), ref(RHS), Name))

frem!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildFRem(ref(builder), ref(LHS), ref(RHS), Name))

# bitwise binary operations

shl!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildShl(ref(builder), ref(LHS), ref(RHS), Name))

lshr!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildLShr(ref(builder), ref(LHS), ref(RHS), Name))

ashr!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildAShr(ref(builder), ref(LHS), ref(RHS), Name))

and!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildAnd(ref(builder), ref(LHS), ref(RHS), Name))

or!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildOr(ref(builder), ref(LHS), ref(RHS), Name))

xor!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildXor(ref(builder), ref(LHS), ref(RHS), Name))

# vector operations

extract_element!(builder::Builder, VecVal::Value, Index::Value, Name::String="") =
    Value(API.LLVMBuildExtractElement(ref(builder), ref(VecVal), ref(Index), Name))

insert_element!(builder::Builder, VecVal::Value, EltVal::Value, Index::Value, Name::String="") =
    Value(API.LLVMBuildInsertElement(ref(builder), ref(VecVal), ref(EltVal), ref(Index), Name))

shuffle_vector!(builder::Builder, V1::Value, V2::Value, Mask::Value, Name::String="") =
    Value(API.LLVMBuildShuffleVector(ref(builder), ref(V1), ref(V2), ref(Mask), Name))

# aggregate operations

extract_value!(builder::Builder, AggVal::Value, Index, Name::String="") =
    Value(API.LLVMBuildExtractValue(ref(builder), ref(AggVal), Index, Name))

insert_value!(builder::Builder, AggVal::Value, EltVal::Value, Index, Name::String="") =
    Value(API.LLVMBuildInsertValue(ref(builder), ref(AggVal), ref(EltVal), Index, Name))

# memory access and addressing operations

alloca!(builder::Builder, Ty::LLVMType, Name::String="") =
    Instruction(API.LLVMBuildAlloca(ref(builder), ref(Ty), Name))

array_alloca!(builder::Builder, Ty::LLVMType, Val::Value, Name::String="") =
    Instruction(API.LLVMBuildArrayAlloca(ref(builder), ref(Ty), ref(Val), Name))

malloc!(builder::Builder, Ty::LLVMType, Name::String="") =
    Instruction(API.LLVMBuildMalloc(ref(builder), ref(Ty), Name))

array_malloc!(builder::Builder, Ty::LLVMType, Val::Value, Name::String="") =
    Instruction(API.LLVMBuildArrayMalloc(ref(builder), ref(Ty), ref(Val), Name))

free!(builder::Builder, PointerVal::Value) =
    Instruction(API.LLVMBuildFree(ref(builder), ref(PointerVal)))

load!(builder::Builder, PointerVal::Value, Name::String="") =
    Instruction(API.LLVMBuildLoad(ref(builder), ref(PointerVal), Name))

store!(builder::Builder, Val::Value, Ptr::Value) =
    Instruction(API.LLVMBuildStore(ref(builder), ref(Val), ref(Ptr)))

fence!(builder::Builder, ordering::API.LLVMAtomicOrdering, singleThread::Core.Bool=false, Name::String="") =
    Instruction(API.LLVMBuildFence(ref(builder), ordering, convert(Bool, singleThread), Name))

atomic_rmw!(builder::Builder, op::API.LLVMAtomicRMWBinOp, PTR::Value, Val::Value, ordering::API.LLVMAtomicOrdering, singleThread::Core.Bool) =
    Instruction(API.LLVMBuildAtomicRMW(ref(builder), op, ref(PTR), ref(Val), ordering, convert(Bool, singleThread)))

atomic_cmpxchg!(builder::Builder, Ptr::Value, Cmp::Value, New::Value, SuccessOrdering::API.LLVMAtomicOrdering, FailureOrdering::API.LLVMAtomicOrdering, SingleThread::Core.Bool) =
    Instruction(API.LLVMBuildAtomicCmpXchg(ref(builder), ref(Ptr), ref(Cmp), ref(New), SuccessOrdering,FailureOrdering, convert(Bool, SingleThread)))

gep!(builder::Builder, Pointer::Value, Indices::Vector{T}, Name::String="") where {T<:Value} =
    Value(API.LLVMBuildGEP(ref(builder), ref(Pointer), ref.(Indices), length(Indices), Name))

inbounds_gep!(builder::Builder, Pointer::Value, Indices::Vector{T}, Name::String="") where {T<:Value} =
    Value(API.LLVMBuildInBoundsGEP(ref(builder), ref(Pointer), ref.(Indices), length(Indices), Name))

struct_gep!(builder::Builder, Pointer::Value, Idx, Name::String="") =
    Value(API.LLVMBuildStructGEP(ref(builder), ref(Pointer), Idx, Name))

# conversion operations

trunc!(builder::Builder, Val::Value, DestTy::LLVMType, Name::String="") =
    Value(API.LLVMBuildTrunc(ref(builder), ref(Val), ref(DestTy), Name))

zext!(builder::Builder, Val::Value, DestTy::LLVMType, Name::String="") =
    Value(API.LLVMBuildZExt(ref(builder), ref(Val), ref(DestTy), Name))

sext!(builder::Builder, Val::Value, DestTy::LLVMType, Name::String="") =
    Value(API.LLVMBuildSExt(ref(builder), ref(Val), ref(DestTy), Name))

fptoui!(builder::Builder, Val::Value, DestTy::LLVMType, Name::String="") =
    Value(API.LLVMBuildFPToUI(ref(builder), ref(Val), ref(DestTy), Name))

fptosi!(builder::Builder, Val::Value, DestTy::LLVMType, Name::String="") =
    Value(API.LLVMBuildFPToSI(ref(builder), ref(Val), ref(DestTy), Name))

uitofp!(builder::Builder, Val::Value, DestTy::LLVMType, Name::String="") =
    Value(API.LLVMBuildUIToFP(ref(builder), ref(Val), ref(DestTy), Name))

sitofp!(builder::Builder, Val::Value, DestTy::LLVMType, Name::String="") =
    Value(API.LLVMBuildSIToFP(ref(builder), ref(Val), ref(DestTy), Name))

fptrunc!(builder::Builder, Val::Value, DestTy::LLVMType, Name::String="") =
    Value(API.LLVMBuildFPTrunc(ref(builder), ref(Val), ref(DestTy), Name))

fpext!(builder::Builder, Val::Value, DestTy::LLVMType, Name::String="") =
    Value(API.LLVMBuildFPExt(ref(builder), ref(Val), ref(DestTy), Name))

ptrtoint!(builder::Builder, Val::Value, DestTy::LLVMType, Name::String="") =
    Value(API.LLVMBuildPtrToInt(ref(builder), ref(Val), ref(DestTy), Name))

inttoptr!(builder::Builder, Val::Value, DestTy::LLVMType, Name::String="") =
    Value(API.LLVMBuildIntToPtr(ref(builder), ref(Val), ref(DestTy), Name))

bitcast!(builder::Builder, Val::Value, DestTy::LLVMType, Name::String="") =
    Value(API.LLVMBuildBitCast(ref(builder), ref(Val), ref(DestTy), Name))

addrspacecast!(builder::Builder, Val::Value, DestTy::LLVMType, Name::String="") =
    Value(API.LLVMBuildAddrSpaceCast(ref(builder), ref(Val), ref(DestTy), Name))

zextorbitcast!(builder::Builder, Val::Value, DestTy::LLVMType, Name::String="") =
    Value(API.LLVMBuildZExtOrBitCast(ref(builder), ref(Val), ref(DestTy), Name))

sextorbitcast!(builder::Builder, Val::Value, DestTy::LLVMType, Name::String="") =
    Value(API.LLVMBuildSExtOrBitCast(ref(builder), ref(Val), ref(DestTy), Name))

truncorbitcast!(builder::Builder, Val::Value, DestTy::LLVMType, Name::String="") =
    Value(API.LLVMBuildTruncOrBitCast(ref(builder), ref(Val), ref(DestTy), Name))

cast!(builder::Builder, Op::API.LLVMOpcode, Val::Value, DestTy::LLVMType, Name::String="") =
    Value(API.LLVMBuildCast(ref(builder), Op, ref(Val), ref(DestTy), Name))

pointercast!(builder::Builder, Val::Value, DestTy::LLVMType, Name::String="") =
    Value(API.LLVMBuildPointerCast(ref(builder), ref(Val), ref(DestTy), Name))

intcast!(builder::Builder, Val::Value, DestTy::LLVMType, Name::String="") =
    Value(API.LLVMBuildIntCast(ref(builder), ref(Val), ref(DestTy), Name))

fpcast!(builder::Builder, Val::Value, DestTy::LLVMType, Name::String="") =
    Value(API.LLVMBuildFPCast(ref(builder), ref(Val), ref(DestTy), Name))

# other operations

icmp!(builder::Builder, Op::API.LLVMIntPredicate, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildICmp(ref(builder), Op, ref(LHS), ref(RHS), Name))

fcmp!(builder::Builder, Op::API.LLVMRealPredicate, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildFCmp(ref(builder), Op, ref(LHS), ref(RHS), Name))

phi!(builder::Builder, Ty::LLVMType, Name::String="") =
    Instruction(API.LLVMBuildPhi(ref(builder), ref(Ty), Name))

select!(builder::Builder, If::Value, Then::Value, Else::Value, Name::String="") =
    Value(API.LLVMBuildSelect(ref(builder), ref(If), ref(Then), ref(Else), Name))

call!(builder::Builder, Fn::Value, Args::Vector{T}=Value[], Name::String="") where {T<:Value} =
    Instruction(API.LLVMBuildCall(ref(builder), ref(Fn), ref.(Args), length(Args), Name))

va_arg!(builder::Builder, List::Value, Ty::LLVMType, Name::String="") =
    Instruction(API.LLVMBuildVAArg(ref(builder), ref(List), ref(Ty), Name))

landingpad!(builder::Builder, Ty::LLVMType, PersFn::Value, NumClauses::Integer, Name::String="") =
    Instruction(API.LLVMBuildLandingPad(ref(builder), ref(Ty), ref(PersFn), NumClauses, Name))

neg!(builder::Builder, V::Value, Name::String="") =
    Value(API.LLVMBuildNeg(ref(builder), ref(V), Name))

nswneg!(builder::Builder, V::Value, Name::String="") =
    Value(API.LLVMBuildNSWNeg(ref(builder), ref(V), Name))

nuwneg!(builder::Builder, V::Value, Name::String="") =
    Value(API.LLVMBuildNUWNeg(ref(builder), ref(V), Name))

fneg!(builder::Builder, V::Value, Name::String="") =
    Value(API.LLVMBuildFNeg(ref(builder), ref(V), Name))

not!(builder::Builder, V::Value, Name::String="") =
    Value(API.LLVMBuildNot(ref(builder), ref(V), Name))

# other build methods

globalstring!(builder::Builder, Str::String, Name::String="") =
    Value(API.LLVMBuildGlobalString(ref(builder), Str, Name))

globalstring_ptr!(builder::Builder, Str::String, Name::String="") =
    Value(API.LLVMBuildGlobalStringPtr(ref(builder), Str, Name))

isnull!(builder::Builder, Val::Value, Name::String="") =
    Value(API.LLVMBuildIsNull(ref(builder), ref(Val), Name))

isnotnull!(builder::Builder, Val::Value, Name::String="") =
    Value(API.LLVMBuildIsNotNull(ref(builder), ref(Val), Name))

ptrdiff!(builder::Builder, LHS::Value, RHS::Value, Name::String="") =
    Value(API.LLVMBuildPtrDiff(ref(builder), ref(LHS), ref(RHS), Name))


text:
module Args

import Random

struct Bencharg
    gen::Function
    minsize::NTuple{N,Int} where N
    tag

    function Bencharg(gen, minsize::NTuple{N,Int}, tag) where N
        @assert hasmethod(gen, typeof(minsize)) "gen must have a method accepting arguments of type Int"
        new(gen, minsize, tag)
    end
end

Bencharg(gen::Function, minsize::Int=1, tag=tag(gen(minsize))) = Bencharg(gen, (minsize,), tag)

show(io::IO, a::Bencharg) = print(io, "$(name(a)) with size ≥ $(a.minsize)")

generate(o, n) = o

function generate(a::Bencharg, n)
    Random.seed!(42)
    a.gen( n... )
end

generate(t::Tuple, n) = map(a -> generate(a, n), t)

smallest(o) = o

smallest(a::Bencharg) = generate(a, a.minsize)

smallest(t::Tuple) = map(smallest, t)

function tag(o)
    t = typeof(o)
    if hasmethod(size, (t,))
        return (type=t, size=size(o))
    end
    return t
end

tag(f::Function) = f

tag(a::Bencharg) = a.tag

tag(t::Tuple) = map(tag, t)

end # module


text:
function quadratic1(a, b, c)
    t(s) = (-b + s*sqrt(Δ))/(2a)
    a == 0 && error("a must be different than zero")
    Δ = Complex(b^2-4*a*c)
    t(1), t(-1)
end

function quadratic2(a, b, c)
    Δ = Complex(b^2-4*a*c)
    t(s) = (-b + s*sqrt(Δ))/(2a)
    a == 0 && error("a must be different than zero")
    t(1), t(-1)
end

max_stars_repo_path: Chapter06/08. Ensuring type stability of your code/quad.jl