text:
# ----------------------------------------------------
# --- BF-OFDM structure
# ----------------------------------------------------
"""
---
BF-OFDM structure
# --- Syntax
- nFBMC		: Number of carriers (PPN size)
- nOFDM		: Number of subbcarriers per carrier (OFDM precoder size)
- K			: Overlapping factor of the PPN
- GI		: CP size of the precoder
- δ			: Rate factor (compression factor)
- allocatedSubcarriers	: Vector of allocated subcarriers
- filterName  : Type of filter used (name)
- filterTaps  : FIR coefficient values
- BT		: Gaussian parameters (if gaussian filter is used)
- filterStopBand : DolphChebyshev attenuation factor (if DC is used)
- fS		: Frqeuency spreading coefficients (if FS filter is used)
- nFFT		: Equivalent OFDM FFT size
- nCP		: Equivalent OFDM CP size
# ---
# v 1.0 - Robin Gerzaguet.
"""
struct StrucBFOFDM<: Waveform
	nFBMC::Int;
	nOFDM::Int;
	K::Int;
	GI::Int;
	δ::Float64;
	allocatedSubcarriers::Array{Int};
	filterName::String;
	filterTaps::Array{Float64};
	BT::Float64;
	filterStopBand::Float64;
	fS::Array{Float64};
	nFFT::Int;
	nCP::Int;
end

"""
---
BF-OFDM initialisation
# --- Syntax
bfofdm = initBFOFDM(nFBMC::Int,nOFDM::Int,K::Int,GI::Int,δ::Float64,allocatedSubcarriers::Array{Int},filterType::String;BT=-1,filterStopBand=-1,fS=[],nFFT=-1,nCP=-1)
# --- Input parameters
- nFBMC		: Number of carriers (PPN size)
- nOFDM		: Number of subbcarriers per carrier (OFDM precoder size)
- K			: Overlapping factor of the PPN
- GI		: CP size of the precoder
- δ			: Rate factor (compression factor)
- allocatedSubcarriers	: Vector of allocated subcarriers
- filterName  : Type of filter used (name)
- filterTaps  : FIR coefficient values
- BT		: Gaussian parameters (if gaussian filter is used)
- filterStopBand : DolphChebyshev attenuation factor (if DC is used)
- fS		: Frqeuency spreading coefficients (if FS filter is used)
- nFFT		: Equivalent OFDM FFT size
- nCP		: Equivalent OFDM CP size
# --- Output parameters
bfofdm		: BF-OFDM structure
# ---
# v 1.0 - Robin Gerzaguet.
"""
function initBFOFDM(nFBMC,nOFDM,K,GI,δ,allocatedSubcarriers,filterType;BT=-1,filterStopBand=-1,fS=[],nFFT=-1,nCP=-1)
	# --- Create direclty the filter taps + filter name.
	(filterTaps,fS) = getBFOFDMFilter(nFBMC,K,filterType,GI=GI,BT=BT,filterStopBand=filterStopBand,fS=fS);
	# --- Create the BF-OFDM structure
	return StrucBFOFDM(nFBMC,nOFDM,K,GI,δ,allocatedSubcarriers,filterType,filterTaps,BT,filterStopBand,fS,nFFT,nCP)
end

# ----------------------------------------------------
# --- Main BF-OFDM signal generation method
# ----------------------------------------------------
"""
---
Creates a Block Filtered - OFDM (BF-OFDM) signal parametrized by its numerlogy and its waveform parameter. Generate a time domain signal based on the input matrix qamMat. The input matrix is a T/F matrix of size (nRe x nbSymb) with nRe the number of allocated subcarrier and nbSymb the number of symbols.

The waveform is parametrized by the PPN size (number of carriers) nFBMC, the size of the precoding stage nOFDM, the CP size of the precoding stage GI and the compression factor parameter δ. See [1], [2] for waveform description and [3,5] for importance and design of compression rate (set to 0.5 in former works on BF-OFDM).

BF-OFDM is characterized by its pulse shape filter. Pulse shape can be automatically tuned to optimal (in terms of  innterference management) [4]. In such a case, filterType should be a string (can be "gaussian_opt", "dc_opt" and "fs_opt"). For a specific pulse shape filter, an array of the nFBMC*K filter taps should be given (and generated for example  with getBFOFDMFilter).
# --- Syntax
sigId	  =bfofdmSigGen(qamMat,nFBMC,nOFDM,K,GI,δ,allocatedSubcarriers,filterType;typeTx="PPN")
# --- Input parameters
- qamMat		: T/F symbols to transmit (QAM symbols) [Array{Complex{Float64},nRe,nbSymb] with nbSymb number of BF-OFDM symbols and nRe number of allocated subcarriers.
- nFBMC		: PPN size
- nOFDM		: OFDM precoding size
- GI			: CP size of precoding stage
- K			: Overlapping factor
- δ			: Compression rate
- allocatedSubcarriers : Vector or allocated subcarriers Array{Int,nRe} with maximum(nRe) < nOFDM*nFBMC*δ
- filterType	: Pulse shape type. Can be a string (see getBFOFDMFilter) or a vector of size K x nFBMC with filter taps.
- typeTx		: Filterbank architecture. Can be "PPN" (polyphase network) or "FS" (frequency spreading) if filterType is of type "fs" or "fs_opt".
# --- Output parameters
- sigId		: Time domainn BF-OFDM signal. [Array{Complex{Float64},nChip}.
# ---
# References
- [1] Gerzaguet, R and Demmer, D and Doré, J-B and Le Ruyet, D. and Kténas, D, "Block-Filtered OFDM: A new Promising Waveform for Multi-service Scenarios", 2017
- [2] Demmer, D and Gerzaguet, R and Doré, J-B and Le Ruyet, D. and Kténas, D, "Block-filtered OFDM: A novel waveform for future wireless technologies", 2017.
- [3] Demmer, D and Rostom, Z and Gerzaguet, R and Doré, J-B and Le Ruyet, D. "Study of OFDM Precoded Filter-Bank Waveforms", 2018
- [4] Demmer, D and Gerzaguet, R and Doré, J-B and Le Ruyet, D. and Kténas, D, "Filter Design for 5G BF-OFDM Waveform", 2017.
- [5] Demmer, D.; Zakaria, R.; Gerzaguet, R.; Doré, J. & Le Ruyet, D. Study of OFDM Precoded Filter-Bank Waveforms, IEEE Transactions on Wireless Communications, 2019.
# v 1.0 - Robin Gerzaguet.
"""
 function bfofdmSigGen(qamMat::Array{Complex{Float64}},nFBMC::Int,nOFDM::Int,K::Int,GI::Int,δ::Float64,allocatedSubcarriers::Array{Int},filterType;typeTx="PPN")
	 # ----------------------------------------------------
	 # --- Overall parameters
	 # ----------------------------------------------------
	 # Set 4 core for FFT computation
	 #FFTW.set_num_threads(4)
	 # --- Define compression alloc
	 nOPB				= Int(nOFDM * δ );
	 nbDataSubcarriers	= length(allocatedSubcarriers);
	 nbSymb				= size(qamMat,2);

	 # ----------------------------------------------------
	 # --- Subcarrier and carrier management
	 # ----------------------------------------------------
	 # We have the subcarrier indexes, as with CP-OFDM
	 # --> We should deduce which carriers is enable (i.e with PPN stage holds data)
	 # --> We also should deduce how the oversampled grid (nOFDM*nFBMC) bears data.
	 # --- Getting PPN alloc.
	 fbmcCarriers	= getBFOFDM_carrierFromSubcarriers(nFBMC, nOFDM, δ , allocatedSubcarriers);
	 nbDataCarriers	= length(fbmcCarriers);
	 # --- Getting oversampling grid index
	 # Indead in BF-OFDM the precoding stage are IFFT blocks of size nOFDM on which nOFDM*δ are allocated.
	 # BF-OFDM can be seen as a global grid of nOFDM*nFBMC, with nRe allocated subcarriers.
	 # Purpose is to get the nRe vector that points which elemnts are allocated on the nOFDM*nFBMC grid
	 allSubcarriers	= getBFOFDM_oversampledGridSubcarriers(nFBMC,nOFDM,δ,fbmcCarriers);
	 # In subcarrier mode, we can have only part of a carrier that is allocated
	 # In this case, size(qamMap,2) is not equal to nOFDM*δ*nbDataCarriers
	 # The given matrix have always lower number of subcarriers, so zero padding approach here
	 # We should remap with 0, and extract only submatrixe of appropriate size
	 if size(qamMat,1) != nOFDM*δ*nbDataCarriers
		 qTmp		= zeros(Complex{Float64},Int(nOFDM*δ*nFBMC),nbSymb);
		 # --- Get all the index associadted to fbmcCarriers alloc.
		 carrierAllmapped = getBFOFDM_subCarrierFromFBMCCarriers(nFBMC, nOFDM, δ, fbmcCarriers)
		 # --- Mapp in the total T/F matrix
		 qTmp[ allocatedSubcarriers,: ]	= qamMat;
		 # --- Getting only full carrier allocated part
		 qamMat = qTmp[ carrierAllmapped,: ];
	 end
	 # At this stage qamMat should be [Complex{Float64},nbDataCarriers*nOFDM*δ,nbSymb]

	 # ----------------------------------------------------
	 # --- Getting filter  and Pre-distortion coefficients
	 # ----------------------------------------------------
	 # --- Time impulse response
	 if isa(filterType,String)
		 # --- Filter is defined as optimal, we launch getBFOFDMFilter without additional parameters
		 # This can only be done with optimal filters, but this is checked in BFOFDM_filters.
		 # Lazy approach here, launch directly without proper check.
		 pulseShape,	  = getBFOFDMFilter(nFBMC,K,filterType,GI=GI);
		 pulseShape	  = 1/(K*nFBMC) .*pulseShape;
	 else
		 # --- Filter is given as filter taps
		 # We check the length
		 if length(filterType) != K * nFBMC
			 error("Unapproriate filter for BF-OFDM: It seems that  filter taps are given, but the length of the filter does not match the PPN expected size (getting %d instead of %d)",length(filterType),K*nFBMC);
		 end
		 # Map to pulseShape
		 pulseShape	  = 1/(K*nFBMC) .* filterType;
	 end
	 # --- Frequency response
	 # Frequency response is necessary to compute the pre-distortion coefficient (based on inverse of frequency response)
	 # Calculate the frequency response of the filter
	 # FFT is done on nOFDM*n FBMC *δ , should be doing zero padding
	 # We assume that nOB > K as otherwise Near Perfect Reconstruction (NPR) cannnot be achieved.
	 pZP	= [pulseShape;zeros(Complex{Float64},nFBMC*nOPB - K * nFBMC)];
	 G		= fft(pZP,1);
	 # --- Calculate pre-distortion coefficient
	 # Get frequency location where subcarrier are effectively mapped and calculate inverse. See [2].
	 zoneInt	= [G[ end-nOPB÷2+1:end ] ;G[ 1:nOPB÷2 ]];
	 Gcorr		= 1 ./abs.(zoneInt).^2 .*conj.(zoneInt);
	 # Same PD coefficients for each FBMC carrier
	 GCorrFreq  = repeat(Gcorr,nbDataCarriers);
	 # --- Apply Predistortion
	 qamMap		= mapslices( x-> x.* GCorrFreq,qamMat,dims=1);

	#----------------------------------------------------
	## --- OFDM preprocessing
	#----------------------------------------------------
	#TODO Adding window for BF-OFDM
	winLenght   = 0;
	# --- OFDM modulation + GI insertion
	output_ifft_full			  = zeros(Complex{Float64},nbDataCarriers,(GI+nOFDM)*nbSymb+winLenght);
	# --- Init a signal with Extended size for windowing
    sizeSymb         = nOFDM+GI;
	sizeExt			 = sizeSymb + winLenght;
	sigExtend		 = zeros(Complex{Float64},nbDataCarriers,sizeExt);
	filterTx         = ones(Int, sizeExt);
	wLD2             = winLenght÷2;
	apodWindow       = transpose(repeat(filterTx,1,nbDataCarriers));
	for idx_j=1:nbSymb
		# --- Getting data of poreprocessing blocks
		input_ifft                                = reshape(qamMap[ :,idx_j ],nOPB ,nbDataCarriers);
		# Mapping carriers depending on delta value
		input_fft_padded                          = zeros(Complex{Float64},nOFDM,nbDataCarriers);
		# --- Map data to appropriate location for each carrier
		input_fft_padded[ allSubcarriers ]        = input_ifft;
		# IFFT preprocessing
		output_ifft			                      = transpose(ifft(input_fft_padded,1));

		#FIXME Speed up with CP insertion with matrix approach.
		# --- Adding CP
		sigExtend[ :,1:GI+wLD2 ]	              = output_ifft[ :,1+end-(GI+wLD2):end ];
		# --- Adding signal
		sigExtend[ :,GI+wLD2 .+ collect(1:nOFDM) ]= output_ifft;
		# --- Adding CS
		sigExtend[ :,end+1-wLD2:end ]             = output_ifft[ 1:wLD2 ];

		# --- Apodisation
		sigApod                                   = sigExtend .* apodWindow;
		# --- Folding in output vector
		output_ifft_full[ :,1 + (idx_j-1)*sizeSymb: idx_j*sizeSymb+ winLenght ] = output_ifft_full[ :,1 + (idx_j-1)*sizeSymb : idx_j*sizeSymb+ winLenght ] .+ sigApod;
	end
	# --- Map preprocessing block to synthesis block
	synthesisMat					= zeros(Complex{Float64},nFBMC, (nOFDM+GI)*nbSymb);
	synthesisMat[fbmcCarriers,:]  	= output_ifft_full;

	# ----------------------------------------------------
	## --- Filter bank stage
	# ----------------------------------------------------
	# --- Init vector
	#TODO Frequency spreading architecture
	# --- Windowing is applied for each symbol
	pKron	= repeat(nFBMC*pulseShape,1,(nOFDM+GI)*nbSymb+winLenght);
	G       =  [];
	# --- Repeat and copy in polyphase filter
	x		= repeat(ifft(synthesisMat,1),K,1);
	# --- Apply filter
	x		= x .* pKron;
	# --- Overlap and Add
	nbFBMCSymbTotal     = (nOFDM+GI) * nbSymb + winLenght;
	y = zeros(Complex{Float64},Int((nbFBMCSymbTotal-1)*nFBMC*δ +nFBMC*K));
	for k=1: nbFBMCSymbTotal
		#FIXME Heavy line in term of computation -> switch to C ?
		y[ (k-1)*Int(nFBMC*δ)  .+ collect(1:( nFBMC*K)) ] = y[ (k-1)*Int(nFBMC*δ)  .+ collect(1:( nFBMC*K)) ] .+ (x[:,k]);
	end

	# ----------------------------------------------------
	# --- Output
	# ----------------------------------------------------
	return y;
 end

 # ----------------------------------------------------
 # ---  Multiple dispatch for structure managment
 # ----------------------------------------------------
 function bfofdmSigGen(qamMat,bfofdm::StrucBFOFDM)
	 return bfofdmSigGen(qamMat,bfofdm.nFBMC,bfofdm.nOFDM,bfofdm.K,bfofdm.GI,bfofdm.δ,bfofdm.allocatedSubcarriers,bfofdm.filterTaps);
 end


text:
# Make an enum to indicate "Source" or "Sink". Move convenient that using symbols
"""
Enum indicating a direction. Values: `Source`, `Sink`.
"""
@enum Direction Source Sink

################################################################################
#                                  PORT TYPES                                  #
################################################################################

"Classification of port types."
@enum PortClass Input Output

struct Port
    name::String

    "The class of this port. Must be a [`PortClass`](@ref)'"
    class::PortClass
    metadata::Dict{String,Any}

    Port(name, class::PortClass; metadata = emptymeta()) = new(name, class, metadata)
end

const _port_compat = Dict(Source => (Input,), Sink => (Output,))

const _port_inverses = Dict(Input => Output, Output => Input)

"""
$(SIGNATURES)

Return `true` if `port` is the correct class for the given [`Direction`](@ref)
"""
checkclass(port::Port, direction::Direction) = port.class in _port_compat[direction]

"""
$(SIGNATURES)

Return a version of `port` with the class inverted.
"""
invert(port::Port) = Port(port.name, _port_inverses[port.class]; metadata = port.metadata)

############
# Port Doc #
############
@doc """
Port type for modeling input/output ports of a `Component`.

API
---
* [`checkclass`](@ref)
* [`invert`](@ref)
""" Port

################################################################################
#                                  LINK TYPE                                   #
################################################################################

struct Link
    name::String
    sources::Vector{Path{Port}}
    dests::Vector{Path{Port}}
    metadata::Dict{String,Any}

    function Link(name, srcs::T, dsts::T, metadata) where {T<:Vector{Path{Port}}}
        return new(name, srcs, dsts, Dict{String,Any}(metadata))
    end
end

"""
$(SIGNATURES)

Return [`Vector{Path{Port}}`](@ref Path) of sources for `link`.
"""
sources(link::Link) = link.sources

"""
$(SIGNATURES)

Return [`Vector{Path{Port}}`](@ref Path) of destinations for `link`.
"""
dests(link::Link) = link.dests

############
# Link Doc #
############
@doc """
    struct Link{P <: AbstractComponentPath}

Link data type for describing which ports are connected. Can have multiple
sources and multiple sinks.

API
---
* [`sources`](@ref)
* [`dests`](@ref)
""" Link

################################################################################
#                               COMPONENT TYPES                                #
################################################################################

# Master abstract type from which all component types will subtype
abstract type AbstractComponent end

"Return an iterator for the children within a component."
children(c::AbstractComponent) = values(c.children)
childnames(component::AbstractComponent) = keys(component.children)
"Return an iterator for links within the component."
links(c::AbstractComponent) = values(c.links)

#-------------------------------------------------------------------------------
# Component
#-------------------------------------------------------------------------------
"""
Basic building block of architecture models. Can be used to construct
hierarchical models.
"""
struct Component <: AbstractComponent
    name::String
    primitive::String
    "Sub-components of this component. Indexed by instance name."
    children::Dict{String,Component}
    "Ports instantiated directly by this component. Indexed by instance name."
    ports::Dict{String,Port}
    "Links instantiated directly by this component. Indexed by instance name."
    links::Dict{String,Link}
    """
    Record of the `Link` (by name) attached to a `Port`, keyed by `Path{Port}`.
    Length of each `Path{Port}` must be 1 or 2, to reference ports either
    instantiated by this component directly, or by one of this component's
    immediate children.j
    """
    portlink::Dict{Path{Port},Link}

    """
    `Dict{String,Any}` for holding any extra data needed by the user.
    """
    metadata::Dict{String,Any}
end

function Component(name; primitive::String = "", metadata = Dict{String,Any}())
    children = Dict{String,Component}()
    ports = Dict{String,Port}()
    links = Dict{String,Link}()
    portlink = Dict{Path{Port},String}()

    # Return the newly constructed type.
    return Component(name, primitive, children, ports, links, portlink, metadata)
end

# Promote types for paths
path_promote(::Type{Component}, ::Type{T}) where {T<:Union{Port,Link}} = T
path_demote(::Type{T}) where {T<:Union{Component,Port,Link}} = Component

# String macros for constructing port and link paths.
macro component_str(s)
    return :(Path{Component}($s))
end
macro link_str(s)
    return :(Path{Link}($s))
end
macro port_str(s)
    return :(Path{Port}($s))
end

function relative_port(component::AbstractComponent, portpath::Path{Port})
    # If the port is defined in the component, just return the port itself
    if length(portpath) == 1
        return component[portpath]
        # If the port is defined one level down in the component hierarchy,
        # extract the port from the level and "invert" it so the directionality
        # of the port is relative to the component "c"
    elseif length(portpath) == 2
        return invert(component[portpath])
    else
        error("Invalid relative port path $portpath for component $(component.name)")
    end
end

#-------------------------------------------------------------------------------
ports(c::Component) = values(c.ports)
ports(c::Component, classes) = Iterators.filter(x -> x.class in classes, values(c.ports))

portpaths(component) = [Path{Port}(name) for name in keys(component.ports)]
function portpaths(component, classes)
    return [Path{Port}(k) for (k, v) in component.ports if v.class in classes]
end
connected_ports(a::AbstractComponent) = collect(keys(a.portlink))

@doc """
    ports(component, [classes])

Return an iterator for all the ports of the given component. Ports of children
are not given. If `classes` are provided, only ports matching the specified
classes will be returned.
""" ports

@doc """
    portpaths(component, [classes])::Vector{Path{Port}}

Return `Path`s to all ports immediately instantiated in `component`.
""" portpaths

#------------------------------------------------------------------------------
# TopLevel
#-------------------------------------------------------------------------------
"""
Top level component for an architecture mode. Main difference is between a
`TopLevel` and a `Component` is that children of a `TopLevel` are accessed
via address instead of instance name. A `TopLevel` also does not have any
ports of its own.

Parameter `D` is the dimensionality of the `TopLevel`.
"""
struct TopLevel{D} <: AbstractComponent
    name::String

    "Direct children of the `TopLevel`, indexed by instance name."
    children::Dict{String,Component}

    "Translation from child instance name to the `Address` that child occupies."
    child_to_address::Dict{String,Address{D}}

    "Translation from `Address` to the `Component` at that address."
    address_to_child::Dict{Address{D},String}

    "`Link`s instantiated directly by the `TopLevel`."
    links::Dict{String,Link}

    "Record of which `Link`s are attached to which `Port`s. Indexed by `Path{Port}."
    portlink::Dict{Path{Port},Link}
    metadata::Dict{String,Any}

    # --Constructor
    function TopLevel{D}(name, metadata = Dict{String,Any}()) where {D}
        # Create a bunch of empty items.
        links = Dict{String,Link}()
        portlink = Dict{Path{Port},Link}()
        children = Dict{String,Component}()
        child_to_address = Dict{String,Address{D}}()
        address_to_child = Dict{Address{D},String}()

        return new{D}(
            name, children, child_to_address, address_to_child, links, portlink, metadata
        )
    end
end

################################################################################
# Convenience methods.
################################################################################

isgloballink(p::Path{Link}) = length(p) == 1
isgloballink(p::Path) = false

isglobalport(p::Path{Port}) = length(p) == 2
isglobalport(p::Path) = false

Base.string(::Type{Component}) = "Component"
Base.string(::Type{Port}) = "Port"
Base.string(::Type{Link}) = "Link"

"""
$(SIGNATURES)

Return an iterator of all addresses with subcomponents in `toplevel`.
"""
addresses(toplevel::TopLevel) = keys(toplevel.address_to_child)
portpaths(toplevel::TopLevel) = Vector{Path{Port}}()

"""
    getaddress(item)

Single argument version: Return an address encapsulated in `item`.

    getaddress(item, index)

Get an address from `item` referenced by `index`.

Method List
-----------
$(METHODLIST)
"""
function getaddress end

"""
    hasaddress(toplevel, path)

Return `true` if the item referenced by `path` has an `Address`.
"""
function hasaddress end

"""
    isaddress(toplevel, address)

Return `true` if `address` exists in `toplevel`.
"""
function isaddress end

getaddress(toplevel::TopLevel, path::Path) = getaddress(toplevel, first(path))
getaddress(toplevel::TopLevel, str::String) = toplevel.child_to_address[str]
hasaddress(toplevel::TopLevel, path::Path) = hasaddress(toplevel, first(path))
hasaddress(toplevel::TopLevel, str::String) = haskey(toplevel.child_to_address, str)
isaddress(toplevel, address) = haskey(toplevel.address_to_child, address)

getname(toplevel::TopLevel, a::Address) = toplevel.address_to_child[a]
function Base.size(t::TopLevel{D}) where {D}
    return dim_max(addresses(t)) .- dim_min(addresses(t)) .+ Tuple(1 for _ in 1:D)
end

"""
    fullsize(toplevel, ruleset)

"""
function fullsize(toplevel::TopLevel{D}, ruleset::RuleSet) where {D}
    # Step 1: Compute the span of addresses.
    address_iter = addresses(toplevel)
    span = dim_max(address_iter) .- dim_min(address_iter) .+ Tuple(1 for _ in 1:D)

    # Run through all addresses of this toplevel. If there is more than one
    # mappable component in a given address, we need to add a dimension to
    # account for that.
    max_mappables = maximum(
        (length ∘ mappables)(toplevel, ruleset, addr) for addr in address_iter
    )
    if max_mappables > 1
        return (max_mappables, span...)
    else
        return span
    end
end

#-------------------------------------------------------------------------------
# Various overloadings of the method "getindex"
#-------------------------------------------------------------------------------

function Base.getindex(c::AbstractComponent, p::Path{T}) where {T<:Union{Port,Link}}
    length(p) == 0 && error("Paths to Ports and Links must have non-zero length")

    c = descend(c, p.steps, length(p) - 1)
    return gettarget(c, T, last(p))
end

Base.getindex(c::AbstractComponent, p::Path{Component}) = descend(c, p.steps, length(p))
Base.getindex(c::AbstractComponent, s::Address) = c.children[c.address_to_child[s]]
function descend(c::AbstractComponent, steps::Vector{String}, n::Integer)
    for i in 1:n
        c = c.children[steps[i]]
    end
    return c
end

gettarget(c::Component, ::Type{Port}, target) = c.ports[target]
gettarget(c::AbstractComponent, ::Type{Link}, target) = c.links[target]

@doc """
    getindex(component, path::Path{T})::T where T <: Union{Port,Link,Component}

Return the architecture type referenced by `path`. Error horribly if `path` does
not exist.

    getindex(toplevel, address)::Component

Return the top level component of `toplevel` at `address`.
""" getindex

################################################################################

"""
    walk_children(component::AbstractComponent, [address]) :: Vector{Path{Component}}

Return relative paths to all the children of `component`. If `address` is given
return relative paths to all components at `address`.
"""
function walk_children(component::AbstractComponent)
    # Recurse on all children.
    paths = Vector{Path{Component}}()
    for (name, child) in component.children
        # Recurse on all children - append each child's instance name in front
        # of the path.
        child_paths = catpath.(name, walk_children(child))
        append!(paths, child_paths)
    end
    # Append an empty path for this component.
    push!(paths, Path{Component}())
    return paths
end

function walk_children(toplevel::TopLevel, address::Address)
    # Walk the component at this address
    paths = walk_children(toplevel[address])
    component_name = toplevel.address_to_child[address]
    # Append the first part of the component path to each of the sub paths.
    return catpath.(component_name, paths)
end

function connected_components(tl::TopLevel{D}) where {D}
    # Construct the associative for the connected components.
    # Use a set to automatically deal with duplicates.
    cc = Dict{Address{D},Set{Address{D}}}()
    # Iterate through all links - record adjacency information
    for link in links(tl)
        for source_port in sources(link), sink_port in dests(link)
            src_address = getaddress(tl, source_port)
            snk_address = getaddress(tl, sink_port)

            push_to_dict(cc, src_address, snk_address)
        end
    end
    # Default unseen addresses to an empty set of addresses.
    for address in addresses(tl)
        if !haskey(cc, address)
            cc[address] = Set{Address{D}}()
        end
    end
    return cc
end

################################################################################
# METHODS FOR NAVIGATING THE HIERARCHY
################################################################################
function search_metadata(c::AbstractComponent, key, value, f::Function = (==))
    return haskey(c.metadata, key) ? f(value, c.metadata[key]) : false
end
search_metadata(c::AbstractComponent, key) = haskey(c.metadata, key)

function search_metadata!(c::AbstractComponent, key, value, f::Function = (==))
    # check top component
    search_metadata(c, key, value, f) && return true
    # recursively call search_metadata! on all subcomponents
    for child in values(c.children)
        search_metadata!(child, key, value, f) && return true
    end
    return false
end

function get_metadata!(c::AbstractComponent, key)
    if haskey(c.metadata, key)
        return c.metadata[key]
    end

    for child in values(c.children)
        val = get_metadata!(child, key)
        if val !== nothing
            return val
        end
    end

    return nothing
end

################################################################################
# ASSERTION METHODS.
################################################################################

function assert_no_children(c::AbstractComponent)
    passed = true
    if length(c.children) != 0
        passed = false
        @error "Cmponent $(c.name) is not expected to have any children."
    end
    return passed
end

assert_no_intrarouting(c::AbstractComponent) = length(c.links) == 0
isfree(c::AbstractComponent, p::Path{Port}) = !haskey(c.portlink, p)

################################################################################
# Documentation for TopLevel
################################################################################

@doc """
    connected_components(toplevel::TopLevel{A,D})::Dict{Address{D}, Set{Address{D}}

Return `d` where key `k` is a valid address of `tl` and where `d[k]` is the set
of valid addresses of `tl` whose components are the destinations of links
originating at address `k`.
""" connected_components

@doc """
    search_metadata(c::AbstractComponent, key, value, f::Function = ==)

Search the metadata of field of `c` for `key`. If `c.metadata[key]` does not
exist, return `false`. Otherwise, return `f(value, c.metadata[key])`.
""" search_metadata

@doc """
    search_metadata!(c::AbstractComponent, key, value, f::Function = ==)

Call `search_metadata` on each subcomponent of `c`. Return `true` if function
call return `true` for any subcomponent.
""" search_metadata!

ext: jl


text:
using Knet: @diff, Param, value, grad, params
using Knet: sigm, tanh, softmax

const in_size        = 10
const hiddens        = (15)
const out_size       = 10

const memory_size    = 5
const turing_hiddens = (8) # TODO

mutable struct Layer

    wri::Param
    wrs::Param
    wrm::Param

    wwi::Param
    wws::Param
    wwm::Param

    wki::Param
    wks::Param
    wkm::Param

    wfi::Param
    wfs::Param
    wfm::Param

    wii::Param
    wis::Param
    wim::Param

    wsi::Param
    wss::Param
    wsm::Param

    w::Param

end

Layer(in_size,layer_size) =
begin

    # read
    wri = Param(randn(in_size, memory_size))
    wrs = Param(randn(layer_size, memory_size))
    wrm = Param(randn(memory_size, memory_size))

    # write
    wwi = Param(randn(in_size, memory_size))
    wws = Param(randn(layer_size, memory_size))
    wwm = Param(randn(memory_size, memory_size))

    # keep
    wki = Param(randn(in_size, layer_size))
    wks = Param(randn(layer_size, layer_size))
    wkm = Param(randn(memory_size, layer_size))

    # forget
    wfi = Param(randn(in_size, layer_size))
    wfs = Param(randn(layer_size, layer_size))
    wfm = Param(randn(memory_size, layer_size))

    # intermediate
    wii = Param(randn(in_size, layer_size))
    wis = Param(randn(layer_size, layer_size))
    wim = Param(randn(memory_size, layer_size))

    # show
    wsi = Param(randn(in_size, layer_size))
    wss = Param(randn(layer_size, layer_size))
    wsm = Param(randn(memory_size, layer_size))

    # intermediate memory
    w = Param(randn(layer_size, memory_size))
    # TODO : hidden params wll be added here.

Layer(wri,wrs,wrm,wwi,wws,wwm,wki,wks,wkm,wfi,wfs,wfm,wii,wis,wim,wsi,wss,wsm,w)
end

(layer::Layer)(in, state, memory) =
begin

    read   = softmax(in * layer.wri + state * layer.wrs + memory * layer.wrm)
    write  = softmax(in * layer.wwi + state * layer.wws + memory * layer.wwm)
    attn_m = read .* memory

    keep   = sigm.(in * layer.wki + state * layer.wks + attn_m * layer.wkm)
    forget = sigm.(in * layer.wfi + state * layer.wfs + attn_m * layer.wfm)
    interm = tanh.(in * layer.wii + state * layer.wis + attn_m * layer.wim)
    show   = sigm.(in * layer.wsi + state * layer.wss + attn_m * layer.wsm)

    interm_m = tanh.(interm * layer.w)
    state    = forget .* state + keep .* interm
    out      = show .* tanh.(state)

    memory += write .* interm_m

(out, state, memory)
end

make() =
begin
    model = []
    hm_layers = length(hiddens)+1
    for i in 1:hm_layers
        if     i == 1         i, o = in_size, hiddens[1]
        elseif i == hm_layers i, o = hiddens[end], out_size
        else                  i, o = hiddens[i-1], hiddens[i]
        end
        push!(model, Layer(i,o))
    end

model
end

zero_states = [zeros(1,l_size) for l_size in (hiddens..., out_size)]
zero_memory = zeros(1,memory_size)

prop(model, data; states=zero_states, memory=zero_memory) =
begin
    response = []
    for timestep in data
        neu_states = []
        for (layer, state) in zip(model, states)
            timestep, state, memory = layer(timestep, state, memory)
            push!(neu_states, state)
        end
        states = neu_states
        push!(response, timestep)
    end

(response, states, memory)
end

loss(seq1, seq2) = sum([sum((s1 - s2).^2) for (s1,s2) in zip(seq1, seq2)])


text:
function crc32(crc::Int, str::String)
    table = zeros(UInt32, 256)

    for i in 0:255
        tmp = i
        for j in 0:7
            if tmp & 1 == 1
                tmp >>= 1
                tmp ⊻= 0xedb88320
            else
                tmp >>= 1
            end
        end

        table[i + 1] = tmp
    end

    crc ⊻= 0xffffffff

    for i in UInt32.(collect(str))
        crc = (crc >> 8) ⊻ table[(crc & 0xff) ⊻ i + 1]
    end

    crc ⊻ 0xffffffff
end

str = "The quick brown fox jumps over the lazy dog"
crc = crc32(0, str)
assert(crc == 0x414fa339)
println("Message: ", str)
println("Checksum: ", hex(crc))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using PlaywithZygoteFlux
using Test
using Zygote

# @testset "PlaywithZygoteFlux.jl" begin
#     # Write your tests here.
# end

f(x::Int) = 0
# f(10)
f'(30)

dx,= gradient(x->2,1)

dx, =gradient(x->x-x+2,1)

gradient(x -> 2*sum(x), [1 2])

# 1+0.1*dx
something(dx,0)

B =
dx, = gradient(A -> reduce(*,A), rand(2,2))
typeof(dx)


text:
using FITSIO
using Base.Test
using Compat

# -----------------------------------------------------------------------------
# Images

# Create a FITS instance and loop over supported types.
fname = tempname() * ".fits"
f = FITS(fname, "w")
for T in [UInt8, Int8, UInt16, Int16, UInt32, Int32, Int64,
          Float32, Float64]
    indata = reshape(T[1:100;], 5, 20)

    # Test writing the data to a new extension
    write(f, indata)

    # test reading the full array
    outdata = read(f[end])
    @test indata == outdata
    @test eltype(indata) == eltype(outdata)

    # test reading subsets of the array
    @test read(f[end], :, :) == indata
    @test read(f[end], 4, 1:10) == indata[4, 1:10]  # 2-d array
    @test read(f[end], :, 4) == indata[:, 4]  # 1-d array
    @test read(f[end], 2, 3) == indata[2, 3]  # scalar
    @test read(f[end], :, 1:2:10) == indata[:, 1:2:10]
    @test read(f[end], 1:3, :) == indata[1:3, :]

    # test expected errors
    @test_throws DimensionMismatch read(f[end], :)
    @test_throws DimensionMismatch read(f[end], :, :, 1)
    @test_throws BoundsError read(f[end], 1:6, :)
    @test_throws BoundsError read(f[end], 1, 0)

end

# test iteration
for hdu in f
    @test size(hdu) == (5, 20)
end
close(f)
if isfile(fname)
    rm(fname)
end

# copy_section()
fname1 = tempname() * ".fits"
f1 = FITS(fname1, "w")
indata = reshape(Float32[1:400;], 20, 20)
write(f1, indata)

fname2 = tempname() * ".fits"
f2 = FITS(fname2, "w")
copy_section(f1[1], f2, 1:10, 1:10)
copy_section(f1[1], f2, 1:10, 1:2:20)
outdata = read(f2[1])
@test outdata == indata[1:10, 1:10]
outdata = read(f2[2])
@test outdata == indata[1:10, 1:2:20]
close(f1)
close(f2)
if isfile(fname1)
    rm(fname1)
end
if isfile(fname2)
    rm(fname2)
end

# -----------------------------------------------------------------------------
# Tables

fname = tempname() * ".fits"
f = FITS(fname, "w")

## Binary table
indata = Dict{Compat.ASCIIString, Array}()
for (i, T) in enumerate([UInt8, Int8, UInt16, Int16, UInt32, Int32, Int64,
                         Float32, Float64, Complex64, Complex128])
    indata["col$i"] = T[1:20;]
end
i = length(indata) + 1
indata["col$i"] = [randstring(10) for j=1:20]  # ASCIIString column
i += 1
indata["col$i"] = ones(Bool, 20)  # Bool column
i += 1
indata["col$i"] = reshape([1:40;], (2, 20))  # vector Int64 column
i += 1
indata["col$i"] = [randstring(5) for j=1:2, k=1:20]  # vector ASCIIString col
indata["vcol1"] = [randstring(j) for j=1:20]  # variable length column
indata["vcol2"] = [collect(1.:j) for j=1.:20.] # variable length

# test writing
write(f, indata; varcols=["vcol1", "vcol2"])

# test reading
for (colname, incol) in indata
    outcol = read(f[2], colname)  # table is in extension 2 (1 = primary hdr)
    @test outcol == incol
    @test eltype(outcol) == eltype(incol)
end

## ASCII tables

indata = Dict{Compat.ASCIIString, Array}()
for (i, T) in enumerate([Int16, Int32, Float32, Float64])
    indata["col$i"] = T[1:20;]
end
i = length(indata) + 1
indata["col$i"] = [randstring(10) for j=1:20]

write(f, indata; hdutype=ASCIITableHDU)

# For ASCII tables, the types don't round trip so we need to define the
# expected output type for each input type.
expected_type = @compat(Dict(Int16=>Int32, Int32=>Int32,
                     Float32=>Float64, Float64=>Float64,
                     Compat.ASCIIString=>Compat.ASCIIString))
for (colname, incol) in indata
    outcol = read(f[3], colname)  # table is in extension 3
    @test outcol == incol
    @test eltype(outcol) == expected_type[eltype(incol)]
end
close(f)
isfile(fname) && rm(fname)

# -----------------------------------------------------------------------------
# FITSHeader

fname = tempname() * ".fits"
f = FITS(fname, "w")

# test that show() works on an empty file and that the beginning and end
# arre what we expect.
io = IOBuffer()
show(io, f)
s = takebuf_string(io)
@test s[1:6] == "File: "
@test s[end-7:end] == "No HDUs."

inhdr = FITSHeader(["FLTKEY", "INTKEY", "BOOLKEY", "STRKEY", "COMMENT",
                    "HISTORY"],
                   [1.0, 1, true, "string value", nothing, nothing],
                   ["floating point keyword",
                    "",
                    "boolean keyword",
                    "string value",
                    "this is a comment",
                    "this is a history"])

inhdr["INTKEY"] = 2  # test setting by key
inhdr[1] = 2.0  # test settting by index
set_comment!(inhdr, "INTKEY", "integer keyword") # test setting a comment

indata = reshape(Float32[1:100;], 5, 20)
write(f, indata; header=inhdr)

# Write a second block.
inhdr2 = deepcopy(inhdr)
inhdr2["INTKEY"] = 3 # Set it to a different value.
write(f, indata; header=inhdr2)

outhdr = read_header(f[1])
@test outhdr["FLTKEY"] === 2.0
@test outhdr["INTKEY"] === 2
@test outhdr["BOOLKEY"] === true
@test outhdr["STRKEY"] == "string value"
@test get_comment(outhdr, 13) == "this is a comment"
@test get_comment(outhdr, 14) == "this is a history"
@test length(outhdr) == 14
@test haskey(outhdr, "FLTKEY")

# Read entire header as a single string
s = read_header(f[1], Compat.ASCIIString)
@test s[1:9] == "SIMPLE  ="  # all headers should start with this.
@test length(s) == (9 + length(inhdr)) * 80  # 9 lines = 8 default + "END"

# Test to check that read_header gets the right block even after reading another.
s_reread = read_header(f[1])
s_reread = read_header(f[2])
s_reread = read_header(f[1], Compat.ASCIIString)
@assert s == s_reread

# Read single keywords
@test read_key(f[1], 9) == ("FLTKEY", 2.0, "floating point keyword")
@test read_key(f[1], "FLTKEY") == (2.0, "floating point keyword")

# Test that show() works and that the beginning of output is what we expect.
io = IOBuffer()
show(io, f)
s = takebuf_string(io)
@test s[1:6] == "File: "

# Clean up from last test.
close(f)
if isfile(fname)
    rm(fname)
end

# -----------------------------------------------------------------------------
# parsing non-standard keyword records

# `create_test_file` : Create a simple FITS file for testing, with the
# given header string added after the required keywords. The length of
# `header` must be a multiple of 80.  The purpose of creating such
# files is to test the parsing of non-standard FITS keyword records
# (non-standard files can't be created with cfitsio).

function create_test_file(fname::AbstractString, header::Compat.ASCIIString)
    if length(header) % 80 != 0
        error("length of header must be multiple of 80")
    end

    f = open(fname, "w")

    stdhdr = "SIMPLE  =                    T / file does conform to FITS standard             BITPIX  =                  -64 / number of bits per data pixel                  NAXIS   =                    2 / number of data axes                            NAXIS1  =                   10 / length of data axis 1                          NAXIS2  =                   10 / length of data axis 2                          EXTEND  =                    T / FITS dataset may contain extensions            "
    endline = "END                                                                             "
    data = fill(0., (10, 10))  # 10x10 array of big-endian Float64 zeros

    # write header
    write(f, stdhdr)
    write(f, header)
    write(f, endline)

    # add padding
    block_position = (length(stdhdr) + length(header) + length(endline)) % 2880
    padding = (block_position == 0) ? 0 : 2880 - block_position
    write(f, " "^padding)

    # write data
    write(f, data)

    # add padding
    block_position = sizeof(data) % 2880
    padding = (block_position == 0) ? 0 : 2880 - block_position
    write(f, fill(0x00, (padding,)))

    close(f)
end

# Create a test file with a few non-standard keyword records
fname = tempname() * ".fits"
header = "        Warning: CROTA2 is inaccurate due to considerable skew                  SKEW    =  1.9511305786508E+00,  1.9234924037208E+00 /Measure of skew           "
create_test_file(fname, header)

# check that we can read the header (and data).
f = FITS(fname)
hdr = read_header(f[1])
data = read(f[1])
close(f)

# test that we can write it back out.
fname2 = tempname() * ".fits"
f2 = FITS(fname2, "w")
write(f2, data; header=hdr)
close(f2)

# clean up.
if isfile(fname)
    rm(fname)
end
if isfile(fname2)
    rm(fname2)
end

# -----------------------------------------------------------------------------
# Miscellaneous

# test that this function works and returns the right type.
@test typeof(FITSIO.libcfitsio_version()) === VersionNumber

println("All tests passed.")


text:
using PatModules
@auxmodule! utils () begin
    export do_nothing

    function do_nothing()
    end
end


text:
# MATLAB

# specific to Julia, that is, related to the output file format (language).

#=
F: read the file containing constants
I:
O:
BA
=#
function m_include_constants_from_literature(src_file_name, pad_string)
  # create src_buffer -
  src_buffer = ""
  # path to distrubtion -
  path_to_src_file = src_file_name
  open(path_to_src_file,"r") do src_file
    for line in eachline(src_file)
      src_buffer *= pad_string*"$line"
    end
  end
  return src_buffer
end

#=
F: generate kinetic buffer
I:
O: return monod affinity constant symbol, mRNA species array,
  protein species array, W_string_array, disassociation_const_string_array.
BA
=#
function m_build_kinetics_buffer(all_species_dict::Dict{String, Int},
  all_rnx_list::Array, all_txtl_dict::Dict, sys2user::Dict)
  kinetics = "function [rnx_rate_vector, TX_rate_vector, TL_rate_vector] = Kinetics(X, data_dictionary)"
            #  "\n\n\t% load all species dictionary" *
            #  "\n\tall_species_dict = data_dictionary(\'all_species_dict\');  % String-->Int"
  # for signaling
  kinetics *= "\n\n\t% generate kinetics for signaling, assume saturation kinetics"
  totalRnxNo = length(all_rnx_list)
  kinetics *= "\n\t% load Monod affinity constants" *
              "\n\tMonodK = data_dictionary(\'MonodAffinityConstantDict\');  % String-->Float"
  MonodAffinityConstant_String_Array = Array{String,1}()  # collect monod affinity constant symbol
  kinetics *= "\n\t% load reaction kinetic constants" *
              "\n\tkcat = data_dictionary(\'kcat_signaling\');  % in order of rnx" *
              "\n\n\t% write reaction rate equations" *
              "\n\trnx_rate_vector = zeros($totalRnxNo, 1);"
  for (index, rnx) in enumerate(all_rnx_list)  # go thru every rnx
    tmp_rnx_rate = "kcat($index)"
    if isdefined(rnx, :catalysts)  # enzyme
      for token in rnx.catalysts
        tmp_rnx_rate *= "*X($(all_species_dict[token.oriBioName]))"
      end
    end
    if isdefined(rnx,:reactants)  # substrate
      for token in rnx.reactants
        tmp_rnx_rate *= "*X($(all_species_dict[token.oriBioName]))/(MonodK(\'MonodK~rnx$(index)~$(token.oriBioName)\')+X($(all_species_dict[token.oriBioName])))"
        push!(MonodAffinityConstant_String_Array, "MonodK~rnx$(index)~$(token.oriBioName)")
      end
    end
    kinetics *= "\n\trnx_rate_vector($index) = "*tmp_rnx_rate * ";"
  end

  # for TXTL
  W_string_array = Array{String,1}()  # W: weight of protein action on transcription
  disassociation_const_string_array = Array{String,1}()  # disassociation constant in tranfer function
  kinetics *= "\n\n\t% generate kinetics & control terms for TXTL;" *
              "\n\t% load data from data dictionary;" *
              "\n\tW_value_dict = data_dictionary(\'W_value_dict\');  % String-->Float" *
              "\n\tcoop = data_dictionary(\'cooperativity\');" *
              "\n\tdisassociation_constant_dict = data_dictionary(\'transferFunctionDisassociationConstantDict\');" *
              "\n\tbackground_control_term = data_dictionary(\'backgroundGeneExpressionControlTermDict\');" *
              "\n\tkcatTX = data_dictionary(\'kcatTranscription\');" *
              "\n\tRNAPconc = data_dictionary(\'RNAPConcentration\');" *
              "\n\tgeneConc = data_dictionary(\'avgGeneConcentration\');" *
              "\n\tsaturationTX = data_dictionary(\'transcriptionSaturationConstant\');" *
              "\n\tkcatTL = data_dictionary(\'kcatTranslation\');" *
              "\n\tRIBOconc = data_dictionary(\'RIBOConcentration\');" *
              "\n\tsaturationTL = data_dictionary(\'translationSaturationConstant\');"

  totalTXTLNo = length(all_txtl_dict)
  kinetics *= "\n\n\t% write control terms, transcription rate equations, and translation rate equations" *
             "\n\tTX_control_term = containers.Map();  % control term" *
             "\n\tTX_rate_vector = containers.Map();  % transcription rate" *
             "\n\tTL_rate_vector = containers.Map();  % translation rate"
  for (key, txtl) in all_txtl_dict # go thru every txtl
    tmp_protein_string = replace(key, sys2user["MRNA"] => sys2user["PROTEIN"], count=1)
    kinetics *= "\n\t% $key and $tmp_protein_string"
    up_factors_array = Array{String,1}()  # collection of upregulation factors name: targetedmRNA_factor(s)
    if !isempty(txtl.activationProtein)  # upregulation
      for token_array in txtl.activationProtein
        tmp_up = ""
        tmp_up_name = key*"_"
        for token in token_array
          tmp_W = "W~$(token.oriBioName)~$(key)"
          tmp_up *= "W_value_dict(\'$(tmp_W)\') * X($(all_species_dict[token.oriBioName]))^coop/((disassociation_constant_dict(\'KD~$(token.oriBioName)\'))^coop+X($(all_species_dict[token.oriBioName]))^coop)*"
          push!(W_string_array, tmp_W)
          push!(disassociation_const_string_array, "KD~$(token.oriBioName)")
          tmp_up_name *= token.oriBioName*"_"
        end
        tmp_up = chop(tmp_up)
        tmp_up_name = chop(tmp_up_name)
        push!(up_factors_array, tmp_up_name)
        kinetics *= "\n\t$(tmp_up_name) = $(tmp_up);"
      end
    end
    down_factors_array = Array{String,1}()  # collection of downregulation factors name: targetedmRNA_factor(s)
    if !isempty(txtl.inhibitionProtein)  # downregulation
      for token_array in txtl.inhibitionProtein
        tmp_down = ""
        tmp_down_name = key*"_"
        for token in token_array
          tmp_W = "W~$(token.oriBioName)~$(key)"
          tmp_down *= "W_value_dict(\'$(tmp_W)\') * X($(all_species_dict[token.oriBioName]))^coop/((disassociation_constant_dict(\'KD~$(token.oriBioName)\'))^coop+X($(all_species_dict[token.oriBioName]))^coop)*"
          push!(W_string_array, tmp_W)
          push!(disassociation_const_string_array, "KD~$(token.oriBioName)")
          tmp_down_name *= token.oriBioName*"_"
        end
        tmp_down = chop(tmp_down)
        tmp_down_name = chop(tmp_down_name)
        push!(down_factors_array, tmp_down_name)
        kinetics *= "\n\t$(tmp_down_name) = $(tmp_down);"
      end
    end
    # combine to form up and down, then merge up and down to form control term
    if !isempty(up_factors_array) && !isempty(down_factors_array)
      up_action = ""
      for up in up_factors_array
        up_action *= up*" +"
      end
      kinetics *= "\n\tup_action_on_$key = $(chop(up_action));"
      down_action = ""
      for down in down_factors_array
        down_action *= down*" +"
      end
      kinetics *= "\n\tdown_action_on_$key = $(chop(down_action));" *
                  "\n\tTX_control_term(\'$key\') = (background_control_term(\'$key\') + up_action_on_$key)/(1 + background_control_term(\'$key\') + up_action_on_$key + down_action_on_$key);"
    elseif !isempty(up_factors_array)
      up_action = ""
      for up in up_factors_array
        up_action *= up*" +"
      end
      kinetics *= "\n\tup_action_on_$key = $(chop(up_action));" *
                  "\n\tTX_control_term(\'$key\') = (background_control_term(\'$key\') + up_action_on_$key)/(1 + background_control_term(\'$key\') + up_action_on_$key);"
    elseif !isempty(down_factors_array)
      down_action = ""
      for down in down_factors_array
        down_action *= down*" +"
      end
      kinetics *= "\n\tdown_action_on_$key = $(chop(down_action));" *
                  "\n\tTX_control_term(\'$key\') = background_control_term(\'$key\')/(1 + background_control_term(\'$key\') + down_action_on_$key);"
    else
      # seems impossible, a TXTL dict is constructed b/c something acts on the key
    end
    kinetics *= "\n\tTX_rate_vector(\'$key\') = TX_control_term(\'$key\') * kcatTX * RNAPconc * geneConc / (saturationTX + geneConc);" *
                "\n\tTL_rate_vector(\'$tmp_protein_string\') = kcatTL * RIBOconc * X($(all_species_dict[key])) / (saturationTL + X($(all_species_dict[key])));"
  end

  kinetics *= "\n\nend"
  return (kinetics, MonodAffinityConstant_String_Array, W_string_array, disassociation_const_string_array)
end

function m_build_data_dictionary_buffer(host_type::AbstractString, all_species_array::Array,
  all_species2index_dict::Dict,
  rnx_species_array::Array, all_rnx_list::Array, all_txtl_dict::Dict,
  Monod_affinity_constant_array::Array, W_string_array::Array,
  disassociation_const_array::Array, mRNA_species_array::Array, protein_species_array::Array)
  # data dictionary
  buffer = "function dataDictionary = DataDictionary()"

  # stoichiometry
  buffer *= "\n\n\t% load stoichiometry"
  buffer *= "\n\tstoichiometric_matrix = load(\'stoichiometry.dat\');"
  No_species = length(all_species_array)
  # initial concentration
  buffer *= "\n\n\t% initial species concentration" *
    "\n\tinitial_condition = [\n\t\t"
  buffer *= join(["0.001, ...  % $i $(all_species_array[i])" for i = 1:(No_species - 1)], "\n\t\t")
  buffer *= "\n\t\t0.001  ...  % $No_species $(all_species_array[No_species])"
  buffer *= "\n\t];"
  # all species dictionary
  buffer *= "\n\n\t% all species dictionary"
  buffer *= "\n\tall_species_dict = containers.Map();  % String --> Int"
  #all_species_dict = Dict{String, Int64}()  # create all species dict: string --> int64
  for (st, index) in all_species2index_dict
    buffer *= "\n\tall_species_dict(\'$st\') = $index;"
  end
  # all species dictionary_reversed
  buffer *= "\n\n\t% all species dictionary_reversed"
  buffer *= "\n\tkey = keys(all_species_dict);" *
            "\n\tval = values(all_species_dict);" *
            "\n\tall_species_reversed_dict = containers.Map(val, key);  % Int --> String"

  # all rnx species array
  buffer *= "\n\n\t% all rnx species array"
  buffer *= "\n\trnx_species_array = {..." *
    join(["\n\t\t\'$st\'" for st in rnx_species_array], ", ...") *
    "  ...\n\t};"
  # all mRNA species array
  buffer *= "\n\n\t% all mRNA species array"
  buffer *= "\n\tmRNA_species_array = {..." *
    join(["\n\t\t\'$st\'" for st in mRNA_species_array], ", ...") *
    "  ...\n\t};"
  # all protein species array
  buffer *= "\n\n\t% all protein species array"
  buffer *= "\n\tprotein_species_array = {..." *
    join(["\n\t\t\'$st\'" for st in protein_species_array], ", ...") *
    "  ...\n\t};"
  # signaling reaction kinetic constants, reaction name: reactant(s)<rnx type:catalyst(s)>product(s)
  buffer *= "\n\n\t% signaling reaction kinetic constants, reaction name: reactant(s)<rnx type:catalyst(s)>product(s)"
  buffer *= "\n\tkcat_signaling = 1e-3*ones(1, $(length(all_rnx_list)));  % kcat[#reaction]: reaction name"
  for (id, rnx) in enumerate(all_rnx_list)
    buffer*= "\n\tkcat_signaling($id) = 1.0;  % kcat: $(rnx.rnxName)"
  end
  # Monod affinity constant in signaling reaction, nomenclature: MonodK_#reaction_Substrate
  buffer *= "\n\n\t% Monod affinity constant in signaling reaction, nomenclature: MonodK_#reaction_Substrate"
  buffer *= "\n\tMonodAffinityConstantDict = containers.Map();"
  for MK in Monod_affinity_constant_array
    # buffer *= "\n\t$MK = 1.0"
    # buffer *= "\n\tMonodAffinityConstant_dict(\'$MK\') = $MK"
    buffer *= "\n\tMonodAffinityConstantDict(\'$MK\') = 1.0;"
  end
  # W_value in transcription control term, nomenclature: W_targetmRNA_actor
  buffer *= "\n\n\n\t% W_value in transcription control term, nomenclature: W_targetmRNA_actor"
  buffer *= "\n\tW_value_dict = containers.Map();"
  for WS in W_string_array
    buffer *= "\n\tW_value_dict(\'$WS\') = 1.0;"
  end
  # disassociation constants in transfer function, nomenclature: KD_speciesName
  buffer *= "\n\n\t% disassociation constants in transfer function, nomenclature: KD_speciesName"
  buffer *= "\n\ttransferFunctionDisassociationConstantDict = containers.Map();"
  for DC in disassociation_const_array
    buffer *= "\n\ttransferFunctionDisassociationConstantDict(\'$DC\') = 1.0;"
  end
  # transcription specific correction factor
  buffer *= "\n\n\t% transcription specific correction factor"
  buffer *= "\n\ttranscriptionSpecificCorrectionFactorDict = containers.Map();"
  for mrna in mRNA_species_array
    buffer *= "\n\ttranscriptionSpecificCorrectionFactorDict(\'$mrna\') = 1.0;  % $mrna"
  end
  # background gene expression control term
  buffer *= "\n\n\t% background gene expression control term"
  buffer *= "\n\tbackgroundGeneExpressionControlTermDict = containers.Map();"
  for mrna in mRNA_species_array
    buffer *= "\n\tbackgroundGeneExpressionControlTermDict(\'$mrna\') = 0.1;  % $mrna"
  end
  # translation specific correction factor
  buffer *= "\n\n\t% translation specific correction factor"
  buffer *= "\n\ttranslationSpecificCorrectionFactor = containers.Map();"
  for pro in protein_species_array
    buffer *= "\n\ttranslationSpecificCorrectionFactor(\'$pro\') = 1.0;  % $pro"
  end
  # cooperativity number in transfer function
  buffer *= "\n\n\t% cooperativity number in transfer function"
  buffer *= "\n\tcooperativity = 1;"
  #buffer *= "\n\n\tbackground_mRNA_synthesis_rate_vector = 0.01*ones($length_TXTL)"
  # load txtl constants buffer
  buffer *= "\n"
  path_head = Base.@__DIR__
  if host_type == "bacteria"
    buffer *= replace(replace(m_include_constants_from_literature(joinpath(path_head,
        "txtl_constants_ecoli.jl"),"\n\t"), "#" => "%"), r"(  | \t|\t\t)% " => "  ; % ")
  else
    buffer *= replace(m_include_constants_from_literature(joinpath(path_head,
        "txtl_constants_hl60.jl"),"\n"), "#" => "%")
  end

  #---------------------------------
  # put all stuff in a Dictionary
  buffer *= "\n\n\t%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
  buffer *= "\n\t% put all stuff in a Dictionary"
  buffer *= "\n\tdataDictionary = containers.Map();"
  buffer *= "\n\tdataDictionary(\'stoichiometric_matrix\') = stoichiometric_matrix;"
  buffer *= "\n\tdataDictionary(\'all_species_dict\') = all_species_dict;"
  buffer *= "\n\tdataDictionary(\'all_species_reversed_dict\') = all_species_reversed_dict;"
  buffer *= "\n\tdataDictionary(\'rnx_species_array\') = rnx_species_array;"
  buffer *= "\n\tdataDictionary(\'mRNA_species_array\') = mRNA_species_array;"
  buffer *= "\n\tdataDictionary(\'protein_species_array\') = protein_species_array;"
  buffer *= "\n\tdataDictionary(\'initial_condition\') = initial_condition;"

  buffer *= "\n\tdataDictionary(\'kcat_signaling\') = kcat_signaling;"
  buffer *= "\n\tdataDictionary(\'MonodAffinityConstantDict\') = MonodAffinityConstantDict;"
  buffer *= "\n\tdataDictionary(\'W_value_dict\') = W_value_dict;"
  buffer *= "\n\tdataDictionary(\'transferFunctionDisassociationConstantDict\') = transferFunctionDisassociationConstantDict;"
  buffer *= "\n\tdataDictionary(\'transcriptionSpecificCorrectionFactorDict\') = transcriptionSpecificCorrectionFactorDict;"
  buffer *= "\n\tdataDictionary(\'backgroundGeneExpressionControlTermDict\') = backgroundGeneExpressionControlTermDict;"
  buffer *= "\n\tdataDictionary(\'translationSpecificCorrectionFactor\') = translationSpecificCorrectionFactor;"
  buffer *= "\n\tdataDictionary(\'cooperativity\') = cooperativity;"
  buffer *= "\n\tdataDictionary(\'RNAPConcentration\') = rnapII_concentration;"
  buffer *= "\n\tdataDictionary(\'RIBOConcentration\') = ribosome_concentration;"
  buffer *= "\n\tdataDictionary(\'degradationConstantmRNA\') = degradation_constant_mRNA;"
  buffer *= "\n\tdataDictionary(\'degradationConstantProtein\') = degradation_constant_protein;"
  buffer *= "\n\tdataDictionary(\'kcatTranscription\') = kcat_transcription;"
  buffer *= "\n\tdataDictionary(\'kcatTranslation\') = kcat_translation;"
  buffer *= "\n\tdataDictionary(\'avgGeneConcentration\') = avg_gene_concentration;"
  buffer *= "\n\tdataDictionary(\'transcriptionSaturationConstant\') = saturation_transcription;"
  buffer *= "\n\tdataDictionary(\'translationSaturationConstant\') = saturation_translation;"
  buffer *= "\n\tdataDictionary(\'specificGrowthRate\') = maximum_specific_growth_rate - death_rate_constant;"

  # buffer *= "\n\n\treturn dataDictionary"
  buffer *= "\n\nend"

  return buffer
end

function m_build_simulation_buffer(NoExtracellularSpecies::Int64)
  # buffer = build_copyright_header_buffer()
  buffer = "\n% set up ODE, get the derivatives
function dydt = Balances(t,y,dataDictionary)
  % correct for negatives
  y(find(y < 0)) = 0.0;

  % load data from data dictionary
  rnx_species = dataDictionary(\'rnx_species_array\');
  mRNA_species = dataDictionary(\'mRNA_species_array\');
  protein_species = dataDictionary(\'protein_species_array\');
  stoichiometry = dataDictionary(\'stoichiometric_matrix\');
  all_species_dict = dataDictionary(\'all_species_dict\');  % String --> number
  specific_growth_rate = dataDictionary(\'specificGrowthRate\');
  transcription_correction = dataDictionary(\'transcriptionSpecificCorrectionFactorDict\');
  translation_correction = dataDictionary(\'translationSpecificCorrectionFactor\');
  degradation_constant_mRNA = dataDictionary(\'degradationConstantmRNA\');
  degradation_constant_protein = dataDictionary(\'degradationConstantProtein\');

  % kinetics
  [rnx_rate, transcription_rate, translation_rate] = Kinetics(y, dataDictionary);

  % calculate the derivatives
  dydt = zeros(1, length(all_species_dict));

  % signaling reaction network
  X2 = zeros(1, (length(rnx_species)-$(NoExtracellularSpecies)));  % species inside the cell
  for i = $(NoExtracellularSpecies+1):length(rnx_species)  % get X2 value from y
    X2(i - $(NoExtracellularSpecies)) = y(all_species_dict(char(rnx_species(i))));
  end
  dX1 = stoichiometry(1:$(NoExtracellularSpecies), :) * rnx_rate;  % extracellular metabolites
  dX2 = stoichiometry($(NoExtracellularSpecies)+1:end, :)*rnx_rate - transpose(specific_growth_rate*X2);  % INTRA
  dX = [dX1; dX2];
  for i = 1:length(rnx_species)  % return derivatives back to dydt
    dydt(all_species_dict(char(rnx_species(i)))) = dX(i);
  end
  dydt(1) = specific_growth_rate * y(1);  % for Biomass

  % TXTL network
  for i = 1:length(mRNA_species)
  % seems easy to add additional terms to account for the signaling effect
    mRNA_id = all_species_dict(char(mRNA_species(i)));
    p_id = all_species_dict(char(protein_species(i)));
    dydt(mRNA_id) = transcription_rate(char(mRNA_species(i))) - (specific_growth_rate + transcription_correction(char(mRNA_species(i)))*degradation_constant_mRNA) * y(mRNA_id);
    dydt(p_id) = translation_rate(char(protein_species(i))) - (specific_growth_rate + translation_correction(char(protein_species(i)))*degradation_constant_protein) * y(p_id);
  end
  dydt = dydt';
end"
  return buffer
end

function m_build_solveODEBalances_buffer(all_species_array::Array,
  all_species_dict::Dict{String, Int}, mRNA_species::Array, protein_species::Array)

  buffer = "function [t, y] = SolveBalances(TStart, TStop, TStep, dataDictionary)" *
    "\n\tt = TStart:TStep:TStop;"
  buffer *= "\n\tall_species_reversed_dict = dataDictionary(\'all_species_reversed_dict\');"
  # initial conditions
  buffer *= "\n\t% initial species concentration"
  buffer *= "\n\ty0 = dataDictionary(\'initial_condition\');"
  # run simulation
  buffer *= "\n\n\t[t, y] = ode23s(@(t,y) Balances(t,y,dataDictionary), t, y0);"
  # plotting
  buffer *= "\n\n\t% plotting results" *
    "\n\tsubplt_col = ceil(length(y0)/4);" *
    "\n\tfor i = 1:length(y0)" *
    "\n\t\tsubplot(4, subplt_col, i)" *
    "\n\t\tplot(t, y(:, i))" *
    "\n\t\ttitle(char(all_species_reversed_dict(i)), 'Interpreter', 'none')" *
    "\n\tend"
  buffer *= "\nend"

  return buffer
end

# FBA data dictionary generation
function m_generate_FBA_data_dictionary(all_rnx_list::Array,
  rnx_species_array::Array, extra_species_num::Int)
  secrete_id = Set{Int64}()  # For initialize the coefficient array
  # data dictionary
  buffer = "function data_dictionary = DataDictionary(time_start,time_stop,time_step)"
  # stoichiometry
  buffer *= "\n\n\t% load stoichiometry"
  buffer *= "\n\tstoichiometric_matrix = load(\'stoichiometry.dat\');"
  # Setup default flux bounds array
  buffer *= "\n\t% Setup default flux bounds array"
  buffer *= "\n\tdefault_bounds_array = ["
  for (id, rnx) in enumerate(all_rnx_list)
    buffer *= "\n\t\t0 100.0; % $(id) $(rnx.rnxName)"
  end
  buffer *= "\n\t];"
  # Setup default species bounds array
  buffer *= "\n\n\t% Setup default species bounds array"
  buffer *= "\n\tspecies_bounds_array = ["
  for (id, sp) in enumerate(rnx_species_array)
    if id > extra_species_num
      buffer *= "\n\t\t0.0 0.0;  % $(id) $(sp)"
    else
      buffer *= "\n\t\t-1.0 1.0; % $(id) $(sp)"
    end
  end
  buffer *= "\n\t];"
  # Setup the objective coefficient array
  buffer *= "\n\n\t% Setup the objective coefficient array"
  buffer *= "\n\tobjective_coefficient_array = ["
  for (id, rnx) in enumerate(all_rnx_list)
    buffer *= "\n\t\t0.0, ... % $(id) $(rnx.rnxName)"
    if rnx.rnxType == "secrete"
        push!(secrete_id, id)
    end
  end
  buffer *= "\n\t];"
  # List of reation strings - used to write flux report
  buffer *= "\n\n\t% List of reation strings - used to write flux report"
  buffer *= "\n\tlist_of_reaction_strings = {..."
  for rnx in all_rnx_list
    buffer *= "\n\t\t\'$(rnx.rnxName)\', ..."
  end
  buffer *= "\n\t};"
  # List of metabolite strings - used to write flux report
  buffer *= "\n\n\t% List of metabolite strings - used to write flux report"
  buffer *= "\n\tlist_of_metabolite_symbols = [..."
  for sp in rnx_species_array
    buffer *= "\n\t\t\'$(sp)\'..."
  end
  buffer *= "\n\t];"
  # return block -
  buffer *= "\n\n"
  buffer *= "\t% =============================== DO NOT EDIT BELOW THIS LINE ============================== #\n"
  buffer *= "\tdata_dictionary = containers.Map();\n"
  buffer *= "\tdata_dictionary(\'stoichiometric_matrix\') = stoichiometric_matrix;\n"
  buffer *= "\tdata_dictionary(\'objective_coefficient_array\') = objective_coefficient_array;\n"
  buffer *= "\tdata_dictionary(\'default_flux_bounds_array\') = default_bounds_array;\n"
  buffer *= "\tdata_dictionary(\'species_bounds_array\') = species_bounds_array;\n"
  buffer *= "\tdata_dictionary(\'list_of_reaction_strings\') = list_of_reaction_strings;\n"
  buffer *= "\tdata_dictionary(\'list_of_metabolite_symbols\') = list_of_metabolite_symbols;\n"
  buffer *= "\tdata_dictionary(\'extra_species_num\') = $(extra_species_num);\n"
  buffer *= "\t% =============================== DO NOT EDIT ABOVE THIS LINE ============================== #\n"
  buffer *= "end\n"

  # maximize_product_dictionary
  buffer2 = "% maximize_product_dictionary"
  buffer2 *= "\nfunction dataDict = maximizeProductDictionary(time_start,time_stop,time_step)"
  # load the data dictionary -
  buffer2 *= "\n\t% load the data dictionary -"
  buffer2 *= "\n\tdataDict = DataDictionary(time_start,time_stop,time_step);"
  # Modify the data dictionary -
  buffer2 *= "\n\t% Modify the data dictionary -"
  buffer2 *= "\n\tobjective_coefficient_array = dataDict(\'objective_coefficient_array\');"
  secrete_id = sort(collect(secrete_id))
  for id in secrete_id
      buffer2 *= "\n\tobjective_coefficient_array($(id)) = -1;"
  end
  buffer2 *= "\n\tdataDict(\'objective_coefficient_array\') = objective_coefficient_array;"
  buffer2 *= "\nend"

  return buffer, buffer2
end


text:
"""
    _SOCtoPSDaff{T}(f::MOI.VectorAffineFunction{T}, g::MOI.ScalarAffineFunction{T})

Builds a VectorAffineFunction representing the upper (or lower) triangular part of the matrix
[ f[1]     f[2:end]' ]
[ f[2:end] g * I     ]
"""
function _SOCtoPSDaff(T::Type, F::Type{<:MOI.AbstractVectorFunction},
                      f::MOI.AbstractVectorFunction, g::MOI.AbstractScalarFunction)
    dim = MOI.output_dimension(f)
    n = div(dim * (dim+1), 2)
    h = MOIU.zero_with_output_dimension(F, n)
    f_scalars = MOIU.eachscalar(f)
    MOIU.operate_output_index!(+, T, trimap(1, 1), h, f_scalars[1])
    for i in 2:dim
        MOIU.operate_output_index!(+, T, trimap(1, i), h, f_scalars[i])
        MOIU.operate_output_index!(+, T, trimap(i, i), h, g)
    end
    return h
end

"""
The `SOCtoPSDBridge` transforms the second order cone constraint ``\\lVert x \\rVert \\le t`` into the semidefinite cone constraints
```math
\\begin{pmatrix}
  t & x^\\top\\\\
  x & tI
\\end{pmatrix} \\succeq 0
```
Indeed by the Schur Complement, it is positive definite iff
```math
\\begin{align*}
  tI & \\succ 0\\\\
  t - x^\\top (tI)^{-1} x & \\succ 0
\\end{align*}
```
which is equivalent to
```math
\\begin{align*}
  t & > 0\\\\
  t^2 & > x^\\top x
\\end{align*}
```
This bridge is not added by default by [`MOI.Bridges.full_bridge_optimizer`](@ref)
as bridging second order cone constraints to semidefinite constraints can be
achieved by the [`SOCRBridge`](@ref) followed by the [`RSOCtoPSDBridge`](@ref)
while creating a smaller semidefinite constraint.
"""
struct SOCtoPSDBridge{T, F, G} <: AbstractBridge
    dim::Int
    cr::MOI.ConstraintIndex{F, MOI.PositiveSemidefiniteConeTriangle}
end
function bridge_constraint(::Type{SOCtoPSDBridge{T, F, G}}, model::MOI.ModelLike, g::G,
                           s::MOI.SecondOrderCone) where {T, F, G}
    d = MOI.dimension(s)
    f = _SOCtoPSDaff(T, F, g, MOIU.eachscalar(g)[1])
    cr = MOI.add_constraint(model, f, MOI.PositiveSemidefiniteConeTriangle(d))
    return SOCtoPSDBridge{T, F, G}(d, cr)
end

function MOI.supports_constraint(
    ::Type{<:SOCtoPSDBridge}, ::Type{<:MOI.AbstractVectorFunction},
    ::Type{MOI.SecondOrderCone})
    return true
end
MOIB.added_constrained_variable_types(::Type{<:SOCtoPSDBridge}) = Tuple{DataType}[]
function MOIB.added_constraint_types(::Type{<:SOCtoPSDBridge{T, F}}) where {T, F}
    return [(F, MOI.PositiveSemidefiniteConeTriangle)]
end
function concrete_bridge_type(::Type{<:SOCtoPSDBridge{T}},
                              G::Type{<:MOI.AbstractVectorFunction},
                              ::Type{MOI.SecondOrderCone}) where T
    F = MOIU.promote_operation(vcat, T, MOIU.scalar_type(G), T)
    return SOCtoPSDBridge{T, F, G}
end

# Attributes, Bridge acting as a model
function MOI.get(
    ::SOCtoPSDBridge{T, F},
    ::MOI.NumberOfConstraints{F, MOI.PositiveSemidefiniteConeTriangle}) where {T, F}
    return 1
end
function MOI.get(
    bridge::SOCtoPSDBridge{T},
    ::MOI.ListOfConstraintIndices{F, MOI.PositiveSemidefiniteConeTriangle}) where {T, F}
    return [bridge.cr]
end

# References
function MOI.delete(model::MOI.ModelLike, c::SOCtoPSDBridge)
    MOI.delete(model, c.cr)
end

# Attributes, Bridge acting as a constraint
function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintFunction,
                 bridge::SOCtoPSDBridge{T, F, G}) where {T, F, G}
    f = MOI.get(model, attr, bridge.cr)
    g = MOIU.eachscalar(f)[trimap.(1, 1:bridge.dim)]
    return MOIU.convert_approx(G, g)
end
function MOI.get(::MOI.ModelLike, ::MOI.ConstraintSet, bridge::SOCtoPSDBridge)
    return MOI.SecondOrderCone(bridge.dim)
end
function MOI.get(model::MOI.ModelLike, a::MOI.ConstraintPrimal, c::SOCtoPSDBridge)
    return MOI.get(model, a, c.cr)[trimap.(1, 1:c.dim)]
end
function MOI.get(model::MOI.ModelLike, a::MOI.ConstraintDual, c::SOCtoPSDBridge)
    dual = MOI.get(model, a, c.cr)
    tdual = sum(i -> dual[trimap(i, i)], 1:c.dim)
    return [tdual; dual[trimap.(2:c.dim, 1)]*2]
end

"""
The `RSOCtoPSDBridge` transforms the second order cone constraint ``\\lVert x \\rVert \\le 2tu`` with ``u \\ge 0`` into the semidefinite cone constraints
```math
\\begin{pmatrix}
  t & x^\\top\\\\
  x & 2uI
\\end{pmatrix} \\succeq 0
```
Indeed by the Schur Complement, it is positive definite iff
```math
\\begin{align*}
  uI & \\succ 0\\\\
  t - x^\\top (2uI)^{-1} x & \\succ 0
\\end{align*}
```
which is equivalent to
```math
\\begin{align*}
  u & > 0\\\\
  2tu & > x^\\top x
\\end{align*}
```
"""
struct RSOCtoPSDBridge{T, F, G} <: AbstractBridge
    dim::Int
    cr::MOI.ConstraintIndex{F, MOI.PositiveSemidefiniteConeTriangle}
end
function bridge_constraint(::Type{RSOCtoPSDBridge{T, F, G}}, model::MOI.ModelLike, g::G,
                           set::MOI.RotatedSecondOrderCone) where {T, F, G}
    dim = MOI.dimension(set) - 1
    g_scalars = MOIU.eachscalar(g)
    h = MOIU.operate!(*, T, g_scalars[2], convert(T, 2))
    f = _SOCtoPSDaff(T, F, g_scalars[[1; 3:MOI.output_dimension(g)]], h)
    cr = MOI.add_constraint(model, f, MOI.PositiveSemidefiniteConeTriangle(dim))
    return RSOCtoPSDBridge{T, F, G}(dim, cr)
end

function MOI.supports_constraint(
    ::Type{<:RSOCtoPSDBridge}, ::Type{<:MOI.AbstractVectorFunction},
    ::Type{MOI.RotatedSecondOrderCone})
    return true
end
MOIB.added_constrained_variable_types(::Type{<:RSOCtoPSDBridge}) = Tuple{DataType}[]
function MOIB.added_constraint_types(::Type{<:RSOCtoPSDBridge{T, F}}) where {T, F}
    return [(F, MOI.PositiveSemidefiniteConeTriangle)]
end
function concrete_bridge_type(::Type{<:RSOCtoPSDBridge{T}},
                              G::Type{<:MOI.AbstractVectorFunction},
                              ::Type{MOI.RotatedSecondOrderCone}) where T
    S = MOIU.scalar_type(G)
    H = MOIU.promote_operation(*, T, T, S)
    F = MOIU.promote_operation(vcat, T, S, H, T)
    return RSOCtoPSDBridge{T, F, G}
end

# Attributes, Bridge acting as a model
function MOI.get(
    ::RSOCtoPSDBridge{T, F},
    ::MOI.NumberOfConstraints{F, MOI.PositiveSemidefiniteConeTriangle}) where {T, F}
    return 1
end
function MOI.get(
    bridge::RSOCtoPSDBridge{T},
    ::MOI.ListOfConstraintIndices{F, MOI.PositiveSemidefiniteConeTriangle}) where {T, F}
    return [bridge.cr]
end

# References
function MOI.delete(model::MOI.ModelLike, bridge::RSOCtoPSDBridge)
    MOI.delete(model, bridge.cr)
end

# Attributes, Bridge acting as a constraint
function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintFunction,
                 bridge::RSOCtoPSDBridge{T, F, G}) where {T, F, G}
    f_scalars = MOIU.eachscalar(MOI.get(model, attr, bridge.cr))
    t = f_scalars[1]
    u = MOIU.operate!(/, T, f_scalars[3], convert(T, 2))
    funcs = [t, u]
    for i in 2:bridge.dim
        push!(funcs, f_scalars[trimap(1, i)])
    end
    return MOIU.convert_approx(G, MOIU.vectorize(funcs))
end
function MOI.get(::MOI.ModelLike, ::MOI.ConstraintSet, bridge::RSOCtoPSDBridge)
    return MOI.RotatedSecondOrderCone(bridge.dim + 1)
end
function MOI.get(model::MOI.ModelLike, a::MOI.ConstraintPrimal, bridge::RSOCtoPSDBridge)
    x = MOI.get(model, MOI.ConstraintPrimal(), bridge.cr)[[trimap(1, 1); trimap(2, 2); trimap.(2:bridge.dim, 1)]]
    x[2] /= 2 # It is (2u*I)[1,1] so it needs to be divided by 2 to get u
    return x
end
function MOI.get(model::MOI.ModelLike, a::MOI.ConstraintDual, bridge::RSOCtoPSDBridge)
    dual = MOI.get(model, MOI.ConstraintDual(), bridge.cr)
    udual = sum(i -> dual[trimap(i, i)], 2:bridge.dim)
    return [dual[1]; 2udual; dual[trimap.(2:bridge.dim, 1)]*2]
end


text:
### A Pluto.jl notebook ###
# v0.15.1

using Markdown
using InteractiveUtils

# ╔═╡ 4cb24885-2f27-400e-a20d-b8152f4c59d3
using LinearAlgebra

# ╔═╡ 5e4375cd-a2a4-4027-88a0-a7ae49fbde1c
# we write a function that returns a root of func in the
# interval [low, high] with relative error tol, assuming such a root exists
function binarySearch(func, low, high, tol)

	# we repeat until the tolorance condition is met
	while abs(high - low)/2 > tol * abs(high + low)/2

		# we define mid as the midpoint of the interval [low, high]
		mid = (low + high)/2

		# we check if the value of func at high and mid have the same sign
		if func(high) * func(mid) > 0

			# and reassign [low, high] accordingly
		    high = mid
		else
		    low = mid
		end
	end

	# we return the avarage of low and high as the approximate root
	return (low + high)/2;
end

# ╔═╡ 4533a6ef-85b5-44e2-a381-f60182886903
# we write an improved version of binarySearch that evaluates func fewer times
# and checks if there is definitely a root in the initial interval
function binarySearchV2(func, low, high, tol)

  # we now also store the values of func at low and high
  f_low = func(low)
  f_high = func(high)

  # we check if the initial interval definitely contains a root
  if f_low * f_high > 0
      error("func(low) and func(high) must have different sign")
  end

  # the while loop is similar to before
  while abs(high - low)/2 > tol * abs(low + high)/2
    mid = (low + high)/2

    # we only have one evaluation of func per loop
    f_mid = func(mid)

    if f_high * f_mid > 0
        high = mid
        f_high = f_mid
        # in this case f_low stays the same
    else
        low = mid
        f_low = f_mid
        # in this case f_high stays the same
    end
  end

  return (low + high)/2
end

# ╔═╡ 977bcb22-8089-49f3-931f-be39cbcd950c
f(x) = exp(x)-4x;

# ╔═╡ 25a35c4d-76aa-406b-b6ba-6cbd960c08b9
zeta = 1e-7;

# ╔═╡ d6406dcf-fda4-4fc3-ba27-8f33ac760c80
md"Both implementations find the same approxiamte root in the interval $(0,1)$."

# ╔═╡ 067c0a97-da55-46b9-9476-396b7d60968f
binarySearch(f, 0, 1, zeta)

# ╔═╡ 8a15ea59-e2fc-4534-84dd-52b3c11abeac
binarySearchV2(f, 0, 1, zeta)

# ╔═╡ 031bddca-90a9-43cd-b5d4-e6ccb41f501e
md"Let's test the case where there is no root in the initial interval. Note that binarySearch does not test whether there is a root it gives a bad answer: it is an **unreliable** function."

# ╔═╡ e53a4200-0d48-4d0a-9b60-ac447eff62d4
binarySearch(f, 1, 2, zeta)

# ╔═╡ 57435047-508b-419a-a1b2-8e5cfa4cbfaf
md"BinarySearchV2 on the other hand includes the initial check."

# ╔═╡ a9fd5fbd-8c6e-45b5-9821-ba64175ec3e0
binarySearchV2(f, 1, 2, zeta)

# ╔═╡ aaf5575d-fb60-43b7-858c-46babc9401b1
md"Let's try it for a more complicated function that involves large matrices. We construct a $n\times n$ matrix with random numbers."

# ╔═╡ ebe48c7f-bae4-4262-b792-0581a0757e0d
n = 4000;

# ╔═╡ 5827d97a-5cf3-4578-90ba-11d07afb84d6
A = rand(n,n);

# ╔═╡ ba659a53-9947-4e2d-b615-c3dcd71b3c18
md"Make a symmetric matrix (for real matrices, M' is the transpose of M)."

# ╔═╡ f997e7f0-26d2-4f90-939b-4776a29d5f66
B = A + A';

# ╔═╡ 25142ed0-b0d2-440a-874b-e49fb7fd198c
g(x) = eigmax(exp.(x*B))-2n;

# ╔═╡ 8dfb5f2b-95ef-4701-b90e-2727aa2b2a32
md"""
This function computes a matrix with elements
$A_{ij} = exp(x B_{ij} )$, computes the largest eigenvalue $\lambda$ of $A$ and returns $\lambda - 2n$.
"""

# ╔═╡ f0fdc7d6-9aa3-4a6c-b2e0-f45706d4ae30
g(0), g(1)

# ╔═╡ fe08f690-93a7-46bb-a110-2cb917d12991
md"""
Since $g$ is continuous and $sign(g(0)) \neq sign(g(1))$, it must have a root in the interval $(0,1)$.
"""

# ╔═╡ 0ff24d6c-9032-4743-9400-f90a70b0bc11
binarySearchV2(g, 0, 1, zeta)

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
"""

# ╔═╡ Cell order:
# ╠═5e4375cd-a2a4-4027-88a0-a7ae49fbde1c
# ╠═4533a6ef-85b5-44e2-a381-f60182886903
# ╠═977bcb22-8089-49f3-931f-be39cbcd950c
# ╠═25a35c4d-76aa-406b-b6ba-6cbd960c08b9
# ╟─d6406dcf-fda4-4fc3-ba27-8f33ac760c80
# ╠═067c0a97-da55-46b9-9476-396b7d60968f
# ╠═8a15ea59-e2fc-4534-84dd-52b3c11abeac
# ╟─031bddca-90a9-43cd-b5d4-e6ccb41f501e
# ╠═e53a4200-0d48-4d0a-9b60-ac447eff62d4
# ╟─57435047-508b-419a-a1b2-8e5cfa4cbfaf
# ╠═a9fd5fbd-8c6e-45b5-9821-ba64175ec3e0
# ╟─aaf5575d-fb60-43b7-858c-46babc9401b1
# ╠═ebe48c7f-bae4-4262-b792-0581a0757e0d
# ╠═5827d97a-5cf3-4578-90ba-11d07afb84d6
# ╟─ba659a53-9947-4e2d-b615-c3dcd71b3c18
# ╠═f997e7f0-26d2-4f90-939b-4776a29d5f66
# ╠═4cb24885-2f27-400e-a20d-b8152f4c59d3
# ╠═25142ed0-b0d2-440a-874b-e49fb7fd198c
# ╟─8dfb5f2b-95ef-4701-b90e-2727aa2b2a32
# ╠═f0fdc7d6-9aa3-4a6c-b2e0-f45706d4ae30
# ╟─fe08f690-93a7-46bb-a110-2cb917d12991
# ╠═0ff24d6c-9032-4743-9400-f90a70b0bc11
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_path: 1a/Root Finding/binaryTest.jl

max_issues_repo_path: 1a/Root Finding/binaryTest.jl

max_forks_repo_path: 1a/Root Finding/binaryTest.jl


text:
# tests for the various forms of stParallel for the ScatteringTransform
# using Revise
using ScatteringTransform
using ContinuousWavelets
using AbstractFFTs, FFTW
using Test, LinearAlgebra, Statistics
using Flux, FourierFilterFlux, CUDA
using Zygote

include("pathTests.jl")
include("fluxtests.jl")


text:
# Define gradients for rand and ones to make the lipschitz loss autodiff'able
Zygote.@adjoint CUDA.ones(x...) = CUDA.ones(x...), _ -> map(_ -> nothing, x)
Zygote.@adjoint CUDA.rand(x...) = CUDA.rand(x...), _ -> map(_ -> nothing, x)

"""
    lipschitz1_gradient_loss(m, x_true, x_generated, batch_size)
Estimates 𝐄ₓ(‖∇ₓD(x)‖₂ - 1)², where x is sampled uniformly on lines between
points from the data distribution and the generators distribution
"""
function lipschitz1_gradient_loss(m, x_true, x_generated, batch_size)
    ξ = CUDA.rand(1, 1, 1, batch_size)
    x_interpolated = ξ .* x_true + (1.0f0 .- ξ) .* x_generated
    _, b = pullback(() -> m(x_interpolated), params(x_interpolated))
    grads = b(CUDA.ones(1, batch_size))
    sqddx = grads[x_interpolated] .^ 2
    mean((sqrt.(sum(sqddx, dims = (1, 2))) .- 1.0f0) .^ 2)
end

"""
    critic_loss(m, x_true, x_generated, batch_size, λ)
WGAN-GP relaxed critic loss with lagrange multiplier λ
"""
function critic_loss(m, x_true, x_generated, batch_size, λ)
    gp = lipschitz1_gradient_loss(m, x_true, x_generated, batch_size)
    return mean(m(x_generated)) - mean(m(x_true)) + λ * gp
end

function generator_loss(m, crit, z)
    -mean(crit(m(z)))
end

"""
    step_critic!(opt, m, x_true, x_generated, batch_size; λ = 10.0f0)
A single optimisation step for the critic, with λ gradient penalty factor.
"""
function step_critic!(opt, m, x_true, x_generated, batch_size; λ = 10.0f0)
    ps = params(m)
    loss, back = pullback(ps) do
        critic_loss(m, x_true, x_generated, batch_size, λ)
    end
    gs = back(1.0f0)
    Flux.update!(opt, ps, gs)
    return loss
end

"""
A single optimisation step for the generator
"""
function step_generator!(opt, m, crit, z)
    ps = params(m)
    loss, back = pullback(ps) do
        generator_loss(m, crit, z)
    end
    gs = back(1.0f0)
    Flux.update!(opt, ps, gs)
    return loss
end

# Alternatively, this could be done outside the pullback in step_critic!
# function interpolate_x(x_true, x_generated, batch_size)
#    ξ = CUDA.rand(1, 1, 1, batch_size)
#    x_interpolated = ξ .* x_true + (1.0f0 .- ξ) .* x_generated
# end
# Zygote.@nograd interpolate_x


text:
using SpecialFunctions, ForwardDiff, FiniteDifferences
include("shared.jl")

const BIG_FD    = central_fdm(10,1)
const BIG_FD_O2 = central_fdm(10,2)

xvbesselkdv(v,x)  = ForwardDiff.derivative(_v->BesselK.adbesselkxv(_v,x), v)
xvbesselkdv2(v,x) = ForwardDiff.derivative(_v->xvbesselkdv(_v,x), v)

beskxv(v,x) = besselk(v,x)*(x^v)
dxvbesselkdv(v, x)  = BIG_FD(_v->beskxv(_v,x), v)
dxvbesselkdv2(v, x) = BIG_FD_O2(_v->beskxv(_v,x), v)

fastdxvbesselkdv2(v, x) = (beskxv(v+2e-6, x) - 2*beskxv(v+1e-6, x) + beskxv(v, x))/1e-12

const VGRID = range(0.25, 10.0, length=101)        # to avoid integer v.
const XGRID = range(0.0,  50.0, length=201)[2:end] # to avoid zero x.

const BASELINE = [dxvbesselkdv2(z[1], z[2])     for z in Iterators.product(VGRID, XGRID)]
const OURSOL   = [xvbesselkdv2(z[1], z[2])      for z in Iterators.product(VGRID, XGRID)]
const FASTFD   = [fastdxvbesselkdv2(z[1], z[2]) for z in Iterators.product(VGRID, XGRID)]
const TOLS     = atolfun.(zip(BASELINE, OURSOL))
const FDTOLS   = atolfun.(zip(BASELINE, FASTFD))
const DIFTOLS  = log10.(TOLS) .- log10.(FDTOLS)

gnuplot_save_matrix!("../plotdata/atols_deriv2_xv_fd.csv",    FDTOLS, VGRID, XGRID)
gnuplot_save_matrix!("../plotdata/atols_deriv2_xv_ad.csv",      TOLS, VGRID, XGRID)
gnuplot_save_matrix!("../plotdata/atols_deriv2_xv_fdad.csv", DIFTOLS, VGRID, XGRID)


text:
# Autogenerated wrapper script for ABC_jll for x86_64-apple-darwin
export abc, abc_rc, libabc

using Readline_jll
JLLWrappers.@generate_wrapper_header("ABC")
JLLWrappers.@declare_executable_product(abc)
JLLWrappers.@declare_file_product(abc_rc)
JLLWrappers.@declare_library_product(libabc, "@rpath/libabc.so")
function __init__()
    JLLWrappers.@generate_init_header(Readline_jll)
    JLLWrappers.@init_executable_product(
        abc,
        "bin/abc",
    )

    JLLWrappers.@init_file_product(
        abc_rc,
        "abc.rc",
    )

    JLLWrappers.@init_library_product(
        libabc,
        "lib/libabc.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module LIBUSB

include("../deps/deps.jl")

const LibusbContext = Cvoid
const LibusbDevice = Cvoid
const LibusbDeviceHandle = Cvoid
const LibusbTransferCbFn = Cvoid

const LIBUSB_LOG_LEVEL_ERROR = Cint(1)
const LIBUSB_LOG_LEVEL_DEBUG = Cint(4)

const LIBUSB_SUCCESS = Cint(0)

const Cuint8  = UInt8
const Cuint16 = UInt16

struct LibusbDeviceDescriptor
    bLength::Cuint8            # Size of this descriptor (in bytes)
    bDescriptorType::Cuint8    # Descriptor type. Will have value \ref libusb_descriptor_type::LIBUSB_DT_DEVICE LIBUSB_DT_DEVICE in this context.
    bcdUSB::Cuint16            # USB specification release number in binary-coded decimal. A value of 0x0200 indicates USB 2.0, 0x0110 indicates USB 1.1, etc.
    bDeviceClass::Cuint8       # USB-IF class code for the device. See \ref libusb_class_code.
    bDeviceSubClass::Cuint8    # USB-IF subclass code for the device, qualified by the bDeviceClass value
    bDeviceProtocol::Cuint8    # USB-IF protocol code for the device, qualified by the bDeviceClass and bDeviceSubClass values
    bMaxPacketSize0::Cuint8    # Maximum packet size for endpoint 0
    idVendor::Cuint16          # USB-IF vendor ID
    idProduct::Cuint16         # USB-IF product ID
    bcdDevice::Cuint16         # Device release number in binary-coded decimal
    iManufacturer::Cuint8      # Index of string descriptor describing manufacturer
    iProduct::Cuint8           # Index of string descriptor describing product
    iSerialNumber::Cuint8      # Index of string descriptor containing device serial number
    bNumConfigurations::Cuint8 # Number of possible configurations
end

const LibusbTransferStatus = Cint # enum
const LIBUSB_TRANSFER_COMPLETED = 0
# 1 LIBUSB_TRANSFER_ERRO
# 2 LIBUSB_TRANSFER_TIMED_OUT
# 3 LIBUSB_TRANSFER_CANCELLED
# 4 LIBUSB_TRANSFER_STALL
# 5 LIBUSB_TRANSFER_NO_DEVICE
# 6 LIBUSB_TRANSFER_OVERFLOW

struct LibusbIsoPacketDescriptor
    length::Cuint # Length of data to request in this packet
    actual_length::Cuint # Amount of data that was actually transferred
    status::LibusbTransferStatus # Status code for this packet
end

mutable struct LibusbTransfer
    dev_handle::Ptr{LibusbDeviceHandle} # Handle of the device that this transfer will be submitted to
    flags::Cuint8 # A bitwise OR combination of \ref libusb_transfer_flags.
    endpoint::Cuchar # Address of the endpoint where this transfer will be sent.
    typ::Cuchar # Type of the endpoint from \ref libusb_transfer_type
    timeout::Cuint # Timeout for this transfer in milliseconds. A value of 0 indicates no timeout.
    status::LibusbTransferStatus
    #= The status of the transfer. Read-only, and only for use within
       transfer callback function.
       If this is an isochronous transfer, this field may read COMPLETED even
       if there were errors in the frames. Use the
       \ref libusb_iso_packet_descriptor::status "status" field in each packet
       to determine if errors occurred. =#
    length::Cint # Length of the data buffer
    actual_length::Cint
    #= Actual length of data that was transferred. Read-only, and only for
       use within transfer callback function. Not valid for isochronous
       endpoint transfers. =#
    callback::Ptr{LibusbTransferCbFn} # Callback function. This will be invoked when the transfer completes, fails, or is cancelled.
    user_data::Ptr{Cvoid} # User context data to pass to the callback function.
    buffer::Ptr{Cuchar} # Data buffer
    num_iso_packets::Cint # Number of isochronous packets. Only used for I/O with isochronous endpoints.
    iso_packet_desc::LibusbIsoPacketDescriptor # Isochronous packet descriptors, for isochronous transfers only.
end

struct Ctimeval
    tv_sec::Clong
    tv_usec::Clong
end

# int LIBUSB_CALL libusb_init(libusb_context **ctx)
function libusb_init(ctx::Ptr{Ptr{LibusbContext}})
    ccall((:libusb_init, libusb), Cint, (Ptr{Ptr{LibusbContext}},), ctx)
end

# libusb_device_handle * LIBUSB_CALL libusb_open_device_with_vid_pid(libusb_context *ctx, uint16_t vendor_id, uint16_t product_id)
function libusb_open_device_with_vid_pid(ctx::Ptr{LibusbContext}, vendor_id::Cuint16, product_id::Cuint16)::Ptr{LibusbDeviceHandle}
    ccall((:libusb_open_device_with_vid_pid, libusb), Ptr{LibusbDeviceHandle}, (Ptr{LibusbContext}, Cuint16, Cuint16), ctx, vendor_id, product_id)
end

# void LIBUSB_CALL libusb_set_debug(libusb_context *ctx, int level)
function libusb_set_debug(ctx::Ptr{LibusbContext}, level::Cint)
    ccall((:libusb_set_debug, libusb), Cvoid, (Ptr{LibusbContext}, Cint), ctx, level)
end

# const char * LIBUSB_CALL libusb_error_name(int errcode)
function libusb_error_name(errcode::Cint)
    error_name = ccall((:libusb_error_name, libusb), Cstring, (Cint,), errcode)
    unsafe_string(error_name)
end

# libusb_device * LIBUSB_CALL libusb_get_device(libusb_device_handle *dev_handle)
function libusb_get_device(dev_handle::Ptr{LibusbDeviceHandle})
    ccall((:libusb_get_device, libusb), Ptr{LibusbDevice}, (Ptr{LibusbDeviceHandle},), dev_handle)
end

# int LIBUSB_CALL libusb_open(libusb_device *dev, libusb_device_handle **dev_handle)
function libusb_open(dev::Ptr{LibusbDevice}, dev_handle::Ptr{Ptr{LibusbDeviceHandle}})
    ccall((:libusb_open, libusb), Cint, (Ptr{LibusbDevice}, Ptr{Ptr{LibusbDeviceHandle}}), dev, dev_handle)
end

# int LIBUSB_CALL libusb_kernel_driver_active(libusb_device_handle *dev_handle, int interface_number)
function libusb_kernel_driver_active(dev_handle::Ptr{LibusbDeviceHandle}, interface_number::Cint)
    ccall((:libusb_kernel_driver_active, libusb), Cint, (Ptr{LibusbDeviceHandle}, Cint), dev_handle, interface_number)
end

# int LIBUSB_CALL libusb_detach_kernel_driver(libusb_device_handle *dev_handle, int interface_number)
function libusb_detach_kernel_driver(dev_handle::Ptr{LibusbDeviceHandle}, interface_number::Cint)
    ccall((:libusb_detach_kernel_driver, libusb), Cint, (Ptr{LibusbDeviceHandle}, Cint), dev_handle, interface_number)
end

# int LIBUSB_CALL libusb_claim_interface(libusb_device_handle *dev_handle, int interface_number)
function libusb_claim_interface(dev_handle::Ptr{LibusbDeviceHandle}, interface_number::Cint)
    ccall((:libusb_claim_interface, libusb), Cint, (Ptr{LibusbDeviceHandle}, Cint), dev_handle, interface_number)
end

# int LIBUSB_CALL libusb_release_interface(libusb_device_handle *dev_handle, int interface_number)
function libusb_release_interface(dev_handle::Ptr{LibusbDeviceHandle}, interface_number::Cint)
    ccall((:libusb_release_interface, libusb), Cint, (Ptr{LibusbDeviceHandle}, Cint), dev_handle, interface_number)
end

# void LIBUSB_CALL libusb_close(libusb_device_handle *dev_handle)
function libusb_close(dev_handle::Ptr{LibusbDeviceHandle})
    ccall((:libusb_close, libusb), Cvoid, (Ptr{LibusbDeviceHandle},), dev_handle)
end

# void LIBUSB_CALL libusb_exit(libusb_context *ctx)
function libusb_exit(ctx::Ptr{LibusbContext})
    ccall((:libusb_exit, libusb), Cvoid, (Ptr{LibusbContext},), ctx)
end

# int LIBUSB_CALL libusb_reset_device(libusb_device_handle *dev_handle)
function libusb_reset_device(dev_handle::Ptr{LibusbDeviceHandle})
    ccall((:libusb_reset_device, libusb), Cint, (Ptr{LibusbDeviceHandle},), dev_handle)
end

# int LIBUSB_CALL libusb_get_device_descriptor(libusb_device *dev, struct libusb_device_descriptor *desc)
function libusb_get_device_descriptor(dev::Ptr{LibusbDevice}, desc::Ref{Ptr{LibusbDeviceDescriptor}})
    ccall((:libusb_get_device_descriptor, libusb), Cint, (Ptr{LibusbDevice}, Ptr{Ptr{LibusbDeviceDescriptor}}), dev, desc)
end

# struct libusb_transfer * LIBUSB_CALL libusb_alloc_transfer(int iso_packets)
function libusb_alloc_transfer(iso_packets::Cint)
    ccall((:libusb_alloc_transfer, libusb), Ptr{LibusbTransfer}, (Cint,), iso_packets)
end

# int LIBUSB_CALL libusb_submit_transfer(struct libusb_transfer *transfer)
function libusb_submit_transfer(transfer::Ptr{LibusbTransfer})
    ccall((:libusb_submit_transfer, libusb), Cint, (Ptr{LibusbTransfer},), transfer)
end

# int LIBUSB_CALL libusb_cancel_transfer(struct libusb_transfer *transfer)
function libusb_cancel_transfer(transfer::Ptr{LibusbTransfer})
    ccall((:libusb_cancel_transfer, libusb), Cint, (Ptr{LibusbTransfer},), transfer)
end

# void LIBUSB_CALL libusb_free_transfer(struct libusb_transfer *transfer)
function libusb_free_transfer(transfer::Ptr{LibusbTransfer})
    ccall((:libusb_free_transfer, libusb), Cvoid, (Ptr{LibusbTransfer},), transfer)
end

# int LIBUSB_CALL libusb_handle_events_completed(libusb_context *ctx, int *completed);
function libusb_handle_events_completed(ctx::Ptr{LibusbContext}, completed::Ptr{Cint})
    ccall((:libusb_handle_events_completed, libusb), Cint, (Ptr{LibusbContext}, Ptr{Cint}), ctx, completed)
end

# int API_EXPORTED libusb_handle_events_timeout_completed(libusb_context *ctx, struct timeval *tv, int *completed)
function libusb_handle_events_timeout_completed(ctx::Ptr{LibusbContext}, tv::Ref{Ctimeval}, completed::Ptr{Cint})
    ccall((:libusb_handle_events_timeout_completed, libusb), Cint, (Ptr{LibusbContext}, Ptr{Ctimeval}, Ptr{Cint}), ctx, tv, completed)
end

const LIBUSB_TRANSFER_TYPE_BULK = Cuchar(2)

# static inline void libusb_fill_bulk_transfer(struct libusb_transfer *transfer, libusb_device_handle *dev_handle, unsigned char endpoint, unsigned char *buffer, int length, libusb_transfer_cb_fn callback, void *user_data, unsigned int timeout)
function libusb_fill_bulk_transfer(transfer::Ptr{LibusbTransfer}, dev_handle::Ptr{LibusbDeviceHandle}, endpoint::Cuchar, buffer::Ptr{Cuchar}, length::Int, callback::Ptr{LibusbTransferCbFn}, user_data::Ptr{Cvoid}, timeout::Cuint)
    t = unsafe_load(transfer)
    t.dev_handle = dev_handle
    t.endpoint = endpoint
    t.typ = LIBUSB_TRANSFER_TYPE_BULK
    t.timeout = timeout
    t.buffer = buffer
    t.length = length
    t.user_data = user_data
    t.callback = callback
    unsafe_store!(transfer, t)
end

end # module LIBUSB


text:
using InteractiveUtils

struct DoubleNumWiltonSauterQStrat{R,S}
    outer_rule_far::R
    inner_rule_far::R
    outer_rule_near::R
    inner_rule_near::R
    sauter_schwab_common_tetr::S
    sauter_schwab_common_face::S
    sauter_schwab_common_edge::S
    sauter_schwab_common_vert::S
end

struct DoubleNumQStrat{R}
    outer_rule::R
    inner_rule::R
end

struct SauterSchwab3DQStrat{R,S}
    outer_rule::R
    inner_rule::R
    sauter_schwab_1D::S
    sauter_schwab_2D::S
    sauter_schwab_3D::S
    sauter_schwab_4D::S
end

defaultquadstrat(op, tfs, bfs) = defaultquadstrat(op, refspace(tfs), refspace(bfs))
macro defaultquadstrat(dop, body)
    @assert dop.head == :tuple
    @assert length(dop.args) == 3
    op = dop.args[1]
    tfs = dop.args[2]
    bfs = dop.args[3]
    ex = quote
        function BEAST.defaultquadstrat(::typeof($op), ::typeof($tfs), ::typeof($bfs))
            $body
        end
    end
    return esc(ex)
end

struct SingleNumQStrat
    quad_rule::Int
end

function quadinfo(op, tfs, bfs; quadstrat=defaultquadstrat(op, tfs, bfs))

    tels, tad = assemblydata(tfs)
    bels, bad = assemblydata(bfs)

    tref = refspace(tfs)
    bref = refspace(bfs)

    i, τ = 1, first(tels)
    j, σ = 1, first(bels)

    @show quadstrat
    println(@which BEAST.quaddata(op,tref,bref,tels,bels,quadstrat))

    qd = quaddata(op,tref,bref,tels,bels,quadstrat)
    println(@which quadrule(op,tref,bref,i,τ,j,σ,qd,quadstrat))

    nothing
end


text:
# All functions needed for the invasion experiments.
# Create the global species pool, carry out species introductions, and disturbance.

let speciespool = Individual[]
    #XXX should `speciespool` be part of the `world` object holding all grid cells?
    """
        initglobalspeciespool!()

    Initialise the foreign species pool (i.e. an array of random individuals that
    are used as the propagule source for invasion events.
    """
    function initglobalspeciespool!()
        for i in 1:setting("global-species-pool")
            push!(speciespool, createind())
        end
    end

    """
        invade!(patch, pressure)

    Select a given amount of individuals from the global species pool and add
    them to the patch.
    """
    function invade!(patch::Patch, pressure::Int)
        invaders = deepcopy(rand(speciespool, pressure))
        append!(patch.community, invaders)
    end

    """
        invade!(world)

    Introduce non-native species from the species pool to all patches marked
    as invasible.
    """
    global function invade!(world::Array{Patch,1})
        # This function has to be global to escape the let-block of the species pool
        (setting("propagule-pressure") == 0 || setting("global-species-pool") == 0) && return
        (length(speciespool) == 0) && initglobalspeciespool!()
        for patch in world
            patch.invasible && invade!(patch, setting("propagule-pressure"))
        end
    end
end

"""
    disturb!(patch, intensity)

Species-independent mortality due to disturbance on the given patch. The intensity
gives the mortality percentage.
"""
function disturb!(patch::Patch, intensity::Int)
    length(patch.community) <= 0 && return
    deaths = Integer(round(length(patch.community) * (intensity/100)))
    dead = unique(rand(1:length(patch.community), deaths))
    while length(dead) < deaths
        newdead = rand(1:length(patch.community), deaths-length(dead))
        dead = unique(vcat(dead, newdead))
    end
    sort!(dead)
    deleteat!(patch.community, dead)
end

"""
    disturb!(world)

Disturb all patches in the world.
"""
function disturb!(world::Array{Patch,1})
    (setting("disturbance") == 0) && return
    if setting("disturbance") > 100
        simlog("disturbance must be no more than 100%", 'w')
        updatesetting("disturbance", 100)
    end
    for patch in world
        (patch.isisland || !setting("static")) && disturb!(patch, setting("disturbance"))
    end
end

max_stars_repo_name: CCTB-Ecomods/gemm

max_issues_repo_name: CCTB-Ecomods/gemm

max_forks_repo_name: CCTB-Ecomods/gemm


text:
# TODO: use StatePrior here, too.

# a gibbs sampler to compare.

# general algorithm:
# * pull out data point x_i
# * compute p(z_i = k | x_i) \propto p(x_i | z_i=k) p(z_i=k)
# * sample z_i from this.

"""
    GibbsCRP

Non-parametric clustering via Gibbs sampling with a Chinese Restaurant Process.

# Constructors

    GibbsCRP(prior, α, x)

Initialize all data assigned to one component.

# Fields

* `prior::NormalInverseChisq`: Prior for each cluster
* `logα::Float64`: (log) pseudo observations reserved for a new cluster
* `components::Vector{Component}`: Clusters
* `empties::Vector{Int}`: Keeps track of which indices are empty and can be
   recycled
* `assignments::Vector{Int}`: Current cluster assignments for each data point
* `data::Vector{Float64}`: Data points to cluster

"""
mutable struct GibbsCRP
    prior::NormalInverseChisq
    logα::Float64
    components::Vector{Component}
    empties::Vector{Int}
    assignments::Vector{Int}
    data::Vector{Float64}
end

GibbsCRP(prior::NormalInverseChisq, α::Float64, x::Vector{Float64}) =
    GibbsCRP(prior,
             log(α),
             [Component(prior, suffstats(Normal, x)), Component(prior)],
             [2],
             ones(Int, length(x)),
             x)

function StatsBase.sample!(gc::GibbsCRP, i::Int)
    x = gc.data[i]
    old_comp = sub(gc.components[gc.assignments[i]], x)
    gc.components[gc.assignments[i]] = old_comp

    # recycle if empty
    isempty(old_comp) && push!(gc.empties, gc.assignments[i])

    log_probs =
        [isempty(comp) ?
           -Inf :
           logpdf(comp, x) + log(nobs(comp))
         for comp
         in gc.components]

    log_probs[gc.empties[1]] = logpdf(Component(gc.prior), x) + gc.logα

    new_k = sample(Weights(exp.(log_probs)))
    gc.components[new_k] = add(gc.components[new_k], x)
    gc.assignments[i] = new_k

    # clean up empties
    if new_k == gc.empties[1]
        popfirst!(gc.empties)
        if isempty(gc.empties)
            push!(gc.components, Component(gc.prior))
            push!(gc.empties, length(gc.components))
        end
    end

    return gc
end

function StatsBase.sample!(gc::GibbsCRP)
    isempty(gc.data) &&
        throw(ArgumentError("can't sample from Gibbs sampler with empty .data field"))
    for i in 1:length(gc.data)
        sample!(gc, i)
    end
    return gc
end

mutable struct GibbsCRPSamples
    g::GibbsCRP
    assignments::Matrix{Int}
    components::Vector{Vector{Component}}
end

function fit!(gc::GibbsCRP, progress=true;
              samples=5000, burnin=samples÷10)
    assignments = Matrix{Int}(length(gc.data), samples)
    components = Vector{Vector{Component}}()
    for _ in 1:burnin
        sample!(gc)
    end

    @showprogress (progress ? 1 : Inf) "Gibbs sampling" for i in 1:samples
        sample!(gc)
        assignments[:,i] .= gc.assignments
        push!(components, copy(gc.components))
    end
    GibbsCRPSamples(gc, assignments, components)
end

assignments(gf::GibbsCRPSamples) = gf.assignments
ncomponents_dist(gf::GibbsCRPSamples) =
    fit(Categorical, mapslices(x -> length(unique(x)), gf.assignments, 1))

marginal_log_posterior(g::GibbsCRP) = marginal_log_posterior(g.components, g.logα)