text:
struct TracerSector{StateType<:AbstractVector,StateTypeA<:AbstractVector,StateTypeB<:AbstractVector}
    # The indexed states in each subregion.  These define what is
    # meant by idx_A and idx_B below.
    indexer_A::UniqueVector{StateTypeA}
    indexer_B::UniqueVector{StateTypeB}

    # Quick lookup of what states are compatible with a given idx_A:
    # by_A[idx_A] contains (idx_B, idx) for each state that looks like
    # idx_A in region A.
    #
    # and vice versa for idx_B:
    # by_B[idx_B] contains a list of (idx_A, idx)
    #
    # NOTE: idx is the index of the full state in the original indexer
    by_A::Vector{Vector{Tuple{Int,Int}}}
    by_B::Vector{Vector{Tuple{Int,Int}}}

    # (idx, idx_A, idx_B) for each state in the sector
    backmap::Vector{Tuple{Int, Int, Int}}

    # To make sure any wavefunctions ψ we are given have the correct
    # number of elements
    original_basis_length::Int
end

struct Tracer{StateType<:AbstractVector,StateTypeA<:AbstractVector,StateTypeB<:AbstractVector}
    sectors::Vector{TracerSector{StateType,StateTypeA,StateTypeB}}

    function Tracer{StateType,StateTypeA,StateTypeB}(sites_A::AbstractVector{Int}, sites_B::AbstractVector{Int}, basis::AbstractVector{StateType}) where {StateType,StateTypeA,StateTypeB}
        nsites = length(sites_A) + length(sites_B)
        # See julia commit 1a28512e1b6cd8a for why array != range
        sort([sites_A; sites_B]) == [1:nsites;] || throw(ArgumentError("All sites must be accounted for"))

        # First figure out and organize all the basis states in A and
        # B, without regard to sector
        preliminary_indexer_A = UniqueVector{StateTypeA}()
        preliminary_indexer_B = UniqueVector{StateTypeB}()
        preliminary_backmap = Tuple{Int, Int}[]
        for state in basis
            length(state) == nsites || throw(ArgumentError("Invalid basis for the given number of sites"))
            idx_A = findfirst!(isequal(state[sites_A]), preliminary_indexer_A)
            idx_B = findfirst!(isequal(state[sites_B]), preliminary_indexer_B)
            push!(preliminary_backmap, (idx_A, idx_B))
        end
        preliminary_by_A = [Int[] for i in 1:length(preliminary_indexer_A)]
        preliminary_by_B = [Int[] for i in 1:length(preliminary_indexer_B)]
        for (idx, (idx_A, idx_B)) in enumerate(preliminary_backmap)
            push!(preliminary_by_A[idx_A], idx_B)
            push!(preliminary_by_B[idx_B], idx_A)
        end

        # Now figure out the independent sectors
        sectors = TracerSector{StateType,StateTypeA,StateTypeB}[]
        remaining_A = BitSet(1:length(preliminary_indexer_A))
        remaining_B = BitSet(1:length(preliminary_indexer_B))
        while !isempty(remaining_A)
            @assert !isempty(remaining_B)

            # Find all states in the current sector
            indexer_A_set = Set{StateTypeA}()
            indexer_B_set = Set{StateTypeB}()
            let
                idx_A_queue = Int[pop!(remaining_A)]
                idx_B_queue = Int[]
                while !(isempty(idx_A_queue) && isempty(idx_B_queue))
                    while !isempty(idx_A_queue)
                        idx_A = pop!(idx_A_queue)
                        @assert idx_A ∉ remaining_A
                        push!(indexer_A_set, preliminary_indexer_A[idx_A])
                        for idx_B in preliminary_by_A[idx_A]
                            if idx_B in remaining_B
                                delete!(remaining_B, idx_B)
                                push!(idx_B_queue, idx_B)
                            end
                        end
                    end
                    while !isempty(idx_B_queue)
                        idx_B = pop!(idx_B_queue)
                        @assert idx_B ∉ remaining_B
                        push!(indexer_B_set, preliminary_indexer_B[idx_B])
                        for idx_A in preliminary_by_B[idx_B]
                            if idx_A in remaining_A
                                delete!(remaining_A, idx_A)
                                push!(idx_A_queue, idx_A)
                            end
                        end
                    end
                end
            end

            # Sort things to be a bit more predictable
            indexer_A = UniqueVector{StateTypeA}(sort!(collect(indexer_A_set), by=(x -> (x...,))))
            indexer_B = UniqueVector{StateTypeB}(sort!(collect(indexer_B_set), by=(x -> (x...,))))

            # Construct by_A, by_B, and backmap
            backmap = Tuple{Int, Int, Int}[]
            for (idx, state) in enumerate(basis)
                if state[sites_A] ∉ indexer_A
                    continue
                end
                idx_A = findfirst(isequal(state[sites_A]), indexer_A)
                idx_B = findfirst(isequal(state[sites_B]), indexer_B)
                push!(backmap, (idx, idx_A, idx_B))
            end
            by_A = [Tuple{Int, Int}[] for i in 1:length(indexer_A)]
            by_B = [Tuple{Int, Int}[] for i in 1:length(indexer_B)]
            for (idx, idx_A, idx_B) in backmap
                push!(by_A[idx_A], (idx_B, idx))
                push!(by_B[idx_B], (idx_A, idx))
            end

            push!(sectors, TracerSector{StateType,StateTypeA,StateTypeB}(indexer_A, indexer_B, by_A, by_B, backmap, length(basis)))
        end
        @assert isempty(remaining_B)

        return new(sectors)
    end
end

subsystem_statetype(::Type{Vector{T}}, ::Integer) where {T} = Vector{T}
subsystem_statetype(::Type{SVector{D,T}}, n::Integer) where {D,T} = SVector{n,T}

# XXX NOTE: Not exactly type stable
function Tracer(hs::HilbertSpace, _sites_A::AbstractVector{Int})
    nsites = length(hs.lattice)

    all(0 .< _sites_A .<= nsites) || throw(ArgumentError("Subsystem contains sites that do not exist in the system."))

    # This first line fixes things up if we specify a site twice
    # in sites_A, or if it's not ordered.
    #
    # FIXME: Perhaps this outside this function
    sites_A = filter(i -> (i in _sites_A), 1:nsites)
    sites_B = setdiff(1:nsites, sites_A)
    @assert issorted(sites_B) # will be true according to setdiff docstring

    StateType = statetype(hs)
    StateTypeA = subsystem_statetype(StateType, length(sites_A))
    StateTypeB = subsystem_statetype(StateType, length(sites_B))

    Tracer{statetype(hs),StateTypeA,StateTypeB}(sites_A, sites_B, hs.indexer)
end

function diagsizes(tracer::Tracer)
    # Returns the number of matrices of each size that must be diagonalized to
    # calculate the QDL diagnostic of a single state
    rv = Dict{Int,Int}()
    for sector in tracer.sectors
        s = length(sector.indexer_A)
        rv[s] = get(rv, s, 0) + 1
    end
    return rv
end

function construct_psimat_block!(psimat::Matrix{T}, ts::TracerSector, ψ::AbstractVector{T}) where {T<:Number}
    length(ψ) == ts.original_basis_length || throw(DimensionMismatch("Wavefunction ψ has the wrong number of elements"))
    size(psimat) == (length(ts.indexer_A), length(ts.indexer_B)) || throw(DimensionMismatch("psimat has the wrong size"))
    fill!(psimat, zero(T))
    for b in 1:length(ts.indexer_B)
        for (a, i) in ts.by_B[b]
            psimat[a, b] = ψ[i]
        end
    end
    return psimat
end

function construct_ρ_A_block!(ρ_A::AbstractMatrix{T}, psimat::Matrix{T}, ts::TracerSector, ψ::AbstractVector{T}) where {T<:Number}
    size(ρ_A) == (length(ts.indexer_A), length(ts.indexer_A)) || throw(DimensionMismatch("Density matrix ρ_A has the wrong size"))
    construct_psimat_block!(psimat, ts, ψ)
    mul!(ρ_A, psimat, psimat')
    # It should be Hermitian by construction (FIXME: right?)
    #@assert sum(abs(ρ_A - ρ_A')) == 0
    return Hermitian(ρ_A)
end

construct_ρ_A_block!(ρ_A::AbstractMatrix{T}, ts::TracerSector, ψ::AbstractVector{T}) where {T<:Number} =
    construct_ρ_A_block!(ρ_A, Array{T}(undef, length(ts.indexer_A), length(ts.indexer_B)), ts, ψ)

function construct_ρ_A_block!(ρ_A::AbstractMatrix{T}, ts::TracerSector, ρ::AbstractMatrix{T}) where {T<:Number}
    size(ρ) == (ts.original_basis_length, ts.original_basis_length) || throw(DimensionMismatch("Density matrix ρ has the wrong size"))
    # FIXME: make sure ρ is Hermitian, or require Hermitian type
    size(ρ_A) == (length(ts.indexer_A), length(ts.indexer_A)) || throw(DimensionMismatch("Density matrix ρ_A has the wrong size"))
    fill!(ρ_A, zero(T))
    for idx_B in 1:length(ts.indexer_B)
        z = ts.by_B[idx_B]
        for (a2, i2) in z
            for (a1, i1) in z
                @inbounds ρ_A[a1, a2] += ρ[i1, i2]
            end
        end
    end
    #@assert sum(abs(ρ_A - ρ_A')) == 0
    return Hermitian(ρ_A)
end

function construct_ρ_A_block(ts::TracerSector, ψ_or_ρ::AbstractVecOrMat{T}) where {T<:Number}
    M = length(ts.indexer_A)
    return construct_ρ_A_block!(Array{T}(undef, M, M), ts, ψ_or_ρ)
end

# Apparently julia does not special-case eigvals() for small, Hermitian
# matrices.  When possible, we return the closed-form solution instead of
# invoking the full eigenvalue machinery.  This can lead to efficiency
# improvements since for some calculations many of our blocks of ρ_A will be
# 1x1 or 2x2.
function myeigvals(mat::Hermitian)
    if size(mat, 1) == 1
        # Trivial
        return [real(mat[1,1])]
    elseif size(mat, 2) == 2
        # Obtain the eigenvalues from the quadratic formula
        a = real(mat[1,1])
        d = real(mat[2,2])
        apd = a + d
        amd = a - d
        desc = sqrt(amd * amd + 4 * abs2(mat[1,2]))
        return [(apd - desc) / 2, (apd + desc) / 2]
    else
        # Call the full eigenvalue machinery
        return eigvals(mat)
    end
end

function entanglement_entropy(tracer::Tracer, ψ_or_ρ::AbstractVecOrMat{<:Number}, alpha::Real=1)
    if alpha == 1
        rv = 0.0
        for sector in tracer.sectors
            ρ_A = construct_ρ_A_block(sector, ψ_or_ρ)
            for v in myeigvals(ρ_A)
                v > 0 || continue
                rv -= v * log(v)
            end
        end
        return rv
    else
        s = 0.0
        for sector in tracer.sectors
            ρ_A = construct_ρ_A_block(sector, ψ_or_ρ)
            for v in myeigvals(ρ_A)
                v > 0 || continue
                s += v ^ alpha
            end
        end
        return log(s) / (1 - alpha)
    end
end

function entanglement_entropy(eigenvalues::AbstractVector{<:Real}, alpha::Real=1)
    if alpha == 1
        return sum(-v * log(v) for v in eigenvalues if v > 0)
    else
        return log(sum(v ^ alpha for v in eigenvalues if v > 0)) / (1 - alpha)
    end
end

ext: jl


text:
export AcquisitionData, kData, kdataSingleSlice, convertUndersampledData,
       changeEncodingSize2D, convert3dTo2d, samplingDensity,
       numContrasts, numChannels, numSlices, numRepititions, apodization!

"""
struct describing MRI acquisition data.

# Fields
* `sequenceInfo::Dict{Symbol,Any}`          - additional information on the pulse sequence
* `traj::Vector{Trajectory}`                - trajectories for each echo/contrast
* `kdata::Array{Matrix{ComplexF64},3}`      - each matrix contains data for one trajectory
                                              (1. dim k-space nodes, 2. dim coils)
                                              the outer dims describe:
                                              1. dim echoes, 2. dim slices, 3. dim repetitions
* `subsampleIndices::Vector{Array{Int64}}`  - indices sampled for each echo/contrast
* `encodingSize::Vector{Int64}`             - size of the underlying image matrix
* `fov::Vector{Float64}`                    - field of view in m
"""
mutable struct AcquisitionData
  sequenceInfo::Dict{Symbol,Any}
  traj::Vector{Trajectory}
  kdata::Array{Matrix{ComplexF64},3}
  subsampleIndices::Vector{Vector{Int64}}
  encodingSize::Vector{Int64}
  fov::Vector{Float64}
end

"""
    numChannels(acqData::AcquisitionData)

returns the number of channels/coils in acqData
"""
numChannels(acqData::AcquisitionData) = size(acqData.kdata[1],2)

"""
    numContrasts(acqData::AcquisitionData)

returns the number of contrasts/echoes in acqData
"""
numContrasts(acqData::AcquisitionData) = size(acqData.kdata,1)

"""
    numSlices(acqData::AcquisitionData)

returns the number of slices in acqData
"""
numSlices(acqData::AcquisitionData) = size(acqData.kdata,2)

"""
    numRepititions(acqData::AcquisitionData)

returns the number of repetitions in acqData
"""
numRepititions(acqData::AcquisitionData) = size(acqData.kdata,3)

"""
    AcquisitionData(tr::T,kdata::Array{Matrix{ComplexF64},3}; seqInfo=Dict{Symbol,Any}()
                        , idx=nothing, encodingSize=Int64[0,0,0], fov=Float64[0,0,0]
                        , kargs...) where T <: Union{Trajectory,Vector{Trajectory}}

constructor for `AcquisitionData`

# Arguments
* `tr <: Union{Trajectory,Vector{Trajectory}}` - trajectories
* `kdata::Array{Matrix{ComplexF64},3}`         - k-space data

the other fields of `AcquisitionData` can be passed as keyword arguments.
"""
function AcquisitionData(tr::T,kdata::Array{Matrix{ComplexF64},3}
                        ; seqInfo=Dict{Symbol,Any}()
                        , idx=nothing
                        , encodingSize=Int64[0,0,0]
                        , fov=Float64[0,0,0]
                        , kargs...) where T <: Union{Trajectory,Vector{Trajectory}}
  tr_vec = vec(tr)
  if idx != nothing
    subsampleIndices = idx
  else
    numContr = size(kdata,1)
    if length(tr_vec) == numContr
      subsampleIndices = [collect(1:size(kspaceNodes(tr_vec[echo]),2)) for echo=1:numContr]
    else
      numSamp = size(kspaceNodes(tr_vec[1]),2)
      subsampleIndices = [collect(1:numSamp) for echo=1:numContr]
    end
  end

  return AcquisitionData(seqInfo,tr_vec,kdata,subsampleIndices,encodingSize,fov)
end

function Images.pixelspacing(acqData::AcquisitionData)
  return [1.0,1.0,1.0]*Unitful.mm
  #return fov./encodingSize*Unitful.mm  #TODO: all needs to be properly initialized
end

"""
    trajectory(acqData::AcquisitionData,i::Int64=1)

returns the `i`-th trajectory contained in `acqData`.
"""
trajectory(acqData::AcquisitionData,i::Int64=1) = acqData.traj[i]

######################
# getting k-space data
######################
"""
    kData(acqData::AcquisitionData, echo::Int64=1, coil::Int64=1, slice::Int64=1;rep::Int64=1)

returns the k-space contained in `acqData` for given `echo`, `coil`, `slice` and `rep`(etition).
"""
function kData(acqData::AcquisitionData, echo::Int64=1, coil::Int64=1, slice::Int64=1;rep::Int64=1)
  return acqData.kdata[echo,slice,rep][:,coil]
end

"""
    multiEchoData(acqData::AcquisitionData, coil::Int64, slice::Int64;rep::Int64=1)

returns the k-space contained in `acqData` for all echoes and given `coil`, `slice` and `rep`(etition).
"""
function multiEchoData(acqData::AcquisitionData, coil::Int64, slice::Int64;rep::Int64=1)
  kdata = ComplexF64[]
  for echo=1:numContrasts(acqData)
    append!(kdata,acqData.kdata[echo,slice,rep][:,coil])
  end
  return kdata
end

"""
    multiCoilData(acqData::AcquisitionData, echo::Int64, slice::Int64;rep::Int64=1)

returns the k-space contained in `acqData` for all coils and given `echo`, `slice` and `rep`(etition).
"""
function multiCoilData(acqData::AcquisitionData, echo::Int64, slice::Int64;rep::Int64=1)
  return vec(acqData.kdata[echo,slice,rep])
end

"""
    multiCoilMultiEchoData(acqData::AcquisitionData, echo::Int64, slice::Int64;rep::Int64=1)

returns the k-space contained in `acqData` for all coils, echoes and given `slice` and `rep`(etition).
"""
function multiCoilMultiEchoData(acqData::AcquisitionData,slice::Int64;rep=1)
  kdata = ComplexF64[]
  for coil=1:numChannels(acqData)
    for echo=1:numContrasts(acqData)
      append!(kdata, acqData.kdata[echo,slice,rep][:,coil])
    end
  end
  return kdata
end

"""
    profileData(acqData::AcquisitionData, echo::Int64, slice::Int64, rep::Int, prof_tr::Int)

returns the profile-data `prof_tr` contained in `acqData` for given `echo`, `coil`, `slice` and `rep`(etition).
"""
function profileData(acqData::AcquisitionData, echo::Int64, slice::Int64, rep::Int, prof_tr::Int)
  tr = trajectory(acqData,echo)
  numSamp, numSl = numSamplingPerProfile(tr), numSlices(tr)
  numChan = numChannels(acqData)
  numProf = div(length(acqData.subsampleIndices[echo]),numSamp) #numProfiles(tr)
  if dims(tr)==2 || numSl==1
    kdata = reshape(multiCoilData(acqData,echo,slice;rep=rep),numSamp,numProf,numChan)
    prof_data = kdata[:,prof_tr,:]
  else
    kdata = reshape(multiCoilData(acqData,echo,1,rep=rep),numSamp,numProf,numSl,numChan)
    prof_data = kdata[:,prof_tr,slice,:]
  end
  return prof_data
end

######################################
# utilities to convert and edit acqData
######################################
"""
    convertUndersampledData(acqData::AcquisitionData)

converts undersampled AcquisitionData, where only profiles contained in
acqData.subsampleIndices are sampled,
into a format where trajectories only contain the sampled profiles.
"""
function convertUndersampledData(acqData::AcquisitionData)

  acqDataSub = deepcopy(acqData)
  numContr = numContrasts(acqData)

  # get number of nodes and reset idx
  numNodes = size(acqData.subsampleIndices,1)
  for echo=1:numContr
    acqDataSub.subsampleIndices[echo] = collect(1:length(acqData.subsampleIndices[echo]))
  end

  # replace trajectories by Undersampled Trajectories
  for i = 1:numContr
    tr = trajectory(acqDataSub,i)
    # assume that coils and slices experience the same trajectory
    tr.nodes = tr.nodes[:,acqData.subsampleIndices[i]]
    tr.cartesian = false
  end

  return acqDataSub
end

##################
# sampling weights
##################
"""
    samplingDensity(acqData::AcquisitionData,shape::Tuple)

returns the sampling density for all trajectories contained in `acqData`.
"""
function samplingDensity(acqData::AcquisitionData,shape::Tuple)
  numContr = numContrasts(acqData)
  weights = Array{Vector{ComplexF64}}(undef,numContr)
  for echo=1:numContr
    tr = trajectory(acqData,echo)
    if isCartesian(tr)
      nodes = kspaceNodes(tr)[:,acqData.subsampleIndices[echo]]
    else
      nodes = kspaceNodes(tr)
    end
    plan = NFFTPlan(nodes, shape,3, 1.25)
    weights[echo] = sqrt.(sdc(plan, iters=3))
  end
  return weights
end

#########################################################################
# convert acqData for a reconstruction with a encodingSize (resolution)
#########################################################################
"""
    changeEncodingSize2D(acqData::AcquisitionData,newEncodingSize::Vector{Int64})

changes the encoding size of 2d encoded `acqData` to `newEncodingSize`.
Returns a new `AcquisitionData` object.
"""
function changeEncodingSize2D(acqData::AcquisitionData,newEncodingSize::Vector{Int64})
  dest = deepcopy(acqData)
  changeEncodingSize2D!(dest,newEncodingSize)
end

"""
    changeEncodingSize2D!(acqData::AcquisitionData,newEncodingSize::Vector{Int64})

does the same thing as `changeEncodingSize2D` but acts in-place on `acqData`.
"""
function changeEncodingSize2D!(acqData::AcquisitionData,newEncodingSize::Vector{Int64})
  fac = acqData.encodingSize[1:2] ./ newEncodingSize[1:2]
  numContr = numContrasts(acqData)
  numSl = numSlices(acqData)
  numReps = numRepititions(acqData)
  idx = Vector{Vector{Int64}}(undef,numContr)

  for i=1:numContr
    tr = trajectory(acqData,i)
    nodes = fac .* kspaceNodes(tr)

    # filter out nodes with magnitude > 0.5
    idxX = findall(x->(x>=-0.5)&&(x<0.5), nodes[1,:])
    idxY = findall(x->(x>=-0.5)&&(x<0.5), nodes[2,:])
    idx[i] = intersect(idxX,idxY, acqData.subsampleIndices[i])

    tr.nodes = nodes[:,idx[i]]
    times = readoutTimes(tr)
    tr.times = times[idx[i]]
    acqData.subsampleIndices[i] = acqData.subsampleIndices[i][idx[i]]
  end

  # find relevant kspace data
  kdata2 = Array{Matrix{ComplexF64}}(undef,numContr,numSl,numReps)
  for rep=1:numReps
    for slice=1:numSl
      for echo=1:numContr
        kdata2[echo,slice,rep] = 1.0/prod(fac) * acqData.kdata[echo,slice,rep][idx[echo],:]
      end
    end
  end
  acqData.kdata = kdata2

  # adjust subsample Indices
  for echo=1:numContr
    acqData.subsampleIndices[echo] = collect(1:length(idx[echo]))
  end

  return acqData
end

"""
    convert3dTo2d(acqData::AcquisitionData)

convert the 3d encoded AcquisitionData `acqData` to the equivalent 2d AcquisitionData.
"""
function convert3dTo2d(acqData::AcquisitionData)
  numContr = numContrasts(acqData)
  numChan = numChannels(acqData)
  numSl = numSlices(trajectory(acqData,1))
  numReps = numRepititions(acqData)
  # check if all trajectories are cartesian
  for i=1:numContr
    if !isCartesian(trajectory(acqData,i))
      @error "conversion to 2d is not supported for non-cartesian data"
    end
  end

  # create 2d trajectories along phase encoding directions
  tr2d = Vector{Trajectory}(undef,numContr)
  for i=1:numContr
    tr3d = trajectory(acqData,i)
    # 1. arg (numProfiles=>y), 2. arg (numSamp=>x)
    tr2d[i] = CartesianTrajectory(numSlices(tr3d),numProfiles(tr3d),TE=echoTime(tr3d),AQ=acqTimePerProfile(tr3d))
  end

  # convert k-space data and place it in the appropriate array structure
  numSamp = numSamplingPerProfile(trajectory(acqData,1)) # assume the same number of samples for all contrasts
  kdata2d = Array{Matrix{ComplexF64}}(undef,numContr,numSamp, numReps)
  for i=1:numContr
    tr = trajectory(acqData,i)
    numProf = div( size(acqData.kdata[i,1,1],1), numSamp ) #numProfiles(tr)
    kdata_i = zeros(ComplexF64, numSamp, numProf, numChan, numReps)
    #convert
    F = 1/sqrt(numSamp)*FFTOp(ComplexF64, (numSamp,))
    for r=1:numReps
      for p=1:numProf # including slices
        for c=1:numChan
          # p_idx = (s-1)*numProf+p
          kdata_i[:,p,c,r] .= adjoint(F) * acqData.kdata[i,1,r][(p-1)*numSamp+1:p*numSamp,c]
        end
      end
    end
    # place kdata in transformed Array structure
    for r=1:numReps
      for j=1:numSamp
        kdata2d[i,j,r] = kdata_i[j,:,:,r] # numSl/numSamp*kdata_i[j,:,:,r]
      end
    end
  end

  # adapt subsampleIndices
  subsampleIndices2d = Vector{Vector{Int64}}(undef, numContr)
  for i=1:numContr
    idx = div.( acqData.subsampleIndices[i] .- 1, numSamp) .+ 1
    subsampleIndices2d[i] = sort(unique(idx))
  end

  return AcquisitionData(acqData.sequenceInfo, tr2d, kdata2d, subsampleIndices2d, acqData.encodingSize, acqData.fov)
end

hann(x) = 0.5*(1-cos(2*pi*(x-0.5)))

function apodization!(acqData::AcquisitionData)
    numContr = numContrasts(acqData)
    numSl = numSlices(acqData)
    numReps = numRepititions(acqData)

    for rep=1:numReps
      for slice=1:numSl
        for echo=1:numContr
          tr = trajectory(acqData,echo)
          nodes = kspaceNodes(tr)
          for k=1:size(acqData.kdata[echo,slice,rep],1)
            weight = hann(nodes[1,k])*hann(nodes[2,k])+0.5
            acqData.kdata[echo,slice,rep][k,:] .*= weight
          end
        end
      end
    end

    return acqData
end


text:
# Allow this to be set manually
function rasterpath()
    if haskey(ENV, "ECODATASOURCES_PATH") && isdir(ENV["ECODATASOURCES_PATH"])
        ENV["ECODATASOURCES_PATH"]
    else
        error("You must set `ENV[\"ECODATASOURCES_PATH\"]` to a path in your system")
    end
end

# function _raster_assets_folder()
#     project_path = dirname(something(Base.current_project(pwd()), Base.load_path_expand(LOAD_PATH[2])))
#     assets_folder = joinpath(project_path, "assets", "general")
#     ispath(assets_folder) || mkpath(assets_folder)
#     return assets_folder
# end

# The folder structure could be more or less deeply nested than this

# function _raster_assets_folder(::Type{TS}, ::Type{TD}) where {TS <: SDMDataSource, TD <: SDMDataSet}
#     project_path = dirname(something(Base.current_project(pwd()), Base.load_path_expand(LOAD_PATH[2])))
#     assets_folder = joinpath(project_path, "assets", string(TS), string(TD))
#     ispath(assets_folder) || mkpath(assets_folder)
#     return assets_folder
# end

function cleanup_assets()
    # May need an "are you sure"? - this could be a lot of GB of data to lose
    ispath(rasterpath()) && rm(rasterpath())
end

function cleanup_assets(T::Type)
    ispath(rasterpath(T)) && rm(rasterpath(T))
end

function cleanup_assets(::Type{TS}, ::Type{TD}) where {TS <: SDMDataSource, TD <: SDMDataSet}
    ispath(_raster_assets_folder(TS, TD)) && rm(_raster_assets_folder(TS, TD); recursive=false)
end


text:
"""
    get_total_costs(supply_chain::SupplyChain)

Gets the total costs of operating the supply chain.
"""
function get_total_costs(supply_chain::SupplyChain)
    check(supply_chain)
    return value(supply_chain.optimization_model[:total_costs])
end

"""
    get_total_fixed_costs(supply_chain::SupplyChain)

Gets the total fixed costs of operating the supply chain.
"""
function get_total_fixed_costs(supply_chain::SupplyChain)
    check(supply_chain)
    return value(supply_chain.optimization_model[:total_fixed_costs])
end

"""
    get_total_transportation_costs(supply_chain::SupplyChain)

Gets the total transportation costs of operating the supply chain.
"""
function get_total_transportation_costs(supply_chain::SupplyChain)
    check(supply_chain)
    return value(supply_chain.optimization_model[:total_transportation_costs])
end

"""
    get_production(supply_chain::SupplyChain, plant::Plant, product::Product, period=1)

Gets the amount of a given product produced at a given plant during a given period.
"""
function get_production(supply_chain::SupplyChain, plant::Plant, product::Product, period=1)
    check(supply_chain)
    return value(supply_chain.optimization_model[:produced][product, plant, period])
end

"""
    get_receipts(supply_chain::SupplyChain, storage::Storage, product::Product, period=1)

Gets the amount of a given product received at a given storage location at a given period.
"""
function get_receipts(supply_chain::SupplyChain, storage::Storage, product::Product, period=1)
    check(supply_chain)
    return value(supply_chain.optimization_model[:received][product, storage, period])
end

"""
    get_shipments(supply_chain::SupplyChain, storage::Storage, product::Product, period=1)

Gets the amount of a given product sent from a given storage location at a given period.
"""
function get_shipments(supply_chain::SupplyChain, storage::Storage, product::Product, period=1)
    check(supply_chain)
    return sum(value(supply_chain.optimization_model[:sent][product, l, period]) for l in get_lanes_out(supply_chain, storage))
end

"""
    get_shipments(supply_chain::SupplyChain, plant::Plant, product::Product, period=1)

Gets the amount of a given product sent from a given plant  at a given period.
"""
function get_shipments(supply_chain::SupplyChain, plant::Plant, product::Product, period=1)
    check(supply_chain)
    return sum(value(supply_chain.optimization_model[:sent][product, l, period]) for l in get_lanes_out(supply_chain, plant))
end

"""
    get_shipments(supply_chain::SupplyChain, supplier::Supplier, product::Product, period=1)

Gets the amount of a given product shipped from a given supplier at a given period.
"""
function get_shipments(supply_chain::SupplyChain, supplier::Supplier, product::Product, period=1)
    check(supply_chain)
    return value(supply_chain.optimization_model[:bought][product, supplier, period])
end

"""
    get_shipments(supply_chain::SupplyChain, lane::Lane, product::Product, period=1)

Gets the amount of a given product sent on a lane at a given period.
"""
function get_shipments(supply_chain::SupplyChain, lane::Lane, product::Product, period=1)
    check(supply_chain)
    return value(supply_chain.optimization_model[:sent][product, lane, period])
end

"""
    is_opened(supply_chain::SupplyChain, storage::Storage, period=1)

Gets whether a given storage location is opened during a given period.
"""
function is_opened(supply_chain::SupplyChain, storage::Storage, period=1)
    check(supply_chain)
    return value(supply_chain.optimization_model[:opened][storage, period]) ≈ 1.0
end

"""
    is_opened(supply_chain::SupplyChain, plant::Plant, period=1)

Gets whether a given plant is opened during a given period.
"""
function is_opened(supply_chain::SupplyChain, plant::Plant, period=1)
    check(supply_chain)
    return value(supply_chain.optimization_model[:opened][plant, period]) ≈ 1.0
end

"""
Gets whether a given storage location is opening during a given period.
"""
function is_opening(supply_chain::SupplyChain, storage::Storage, period=1)
    check(supply_chain)
    return value(supply_chain.optimization_model[:opening][storage, period]) ≈ 1.0
end

"""
Gets whether a given storage location is closing during a given period.
"""
function is_closing(supply_chain::SupplyChain, storage::Storage, period=1)
    check(supply_chain)
    return value(supply_chain.optimization_model[:closing][storage, period]) ≈ 1.0
end

"""
Gets the inventory of a product stored at the start of a period.
"""
function get_inventory_at_start(supply_chain::SupplyChain, storage::Storage, product::Product, period=1)
    check(supply_chain)
    return value(supply_chain.optimization_model[:stored_at_start][product, storage, period])
end

"""
Gets the inventory of a product stored at the end of a period.
"""
function get_inventory_at_end(supply_chain::SupplyChain, storage::Storage, product::Product, period=1)
    check(supply_chain)
    return value(supply_chain.optimization_model[:stored_at_end][product, storage, period])
end

function check(supply_chain)
    if isnothing(supply_chain.optimization_model) ||
        !((termination_status(supply_chain.optimization_model) == JuMP.MathOptInterface.OPTIMAL) || (primal_status(supply_chain.optimization_model) == JuMP.MathOptInterface.FEASIBLE_POINT))
        throw(ErrorException("The optimize_network! function must be called first."))
    end
end


text:
# a simple LP example
#
#   maximize x + y
#
#   s.t. 50 x + 24 y <= 2400
#        30 x + 33 y <= 2100
#        x >= 45, y >= 5
#
#   solution: x = 45, y = 6.25, objv = 51.25

using Xpress

#method = getparam(env, "Method")
#println("method = $method")

model = Xpress.Model("lp_01", :maximize)

# add variables
add_cvar!(model, 1.0, 45., Inf)  # x
add_cvar!(model, 1.0,  5., Inf)  # y
#update_model!(model)

# add constraints
add_constr!(model, [50., 24.], '<', 2400.)
add_constr!(model, [30., 33.], '<', 2100.)
#update_model!(model)

println(model)

# perform optimization
optimize(model)

# show results
info = get_optiminfo(model)
println()
println(info)

sol = get_solution(model)
println("soln = $(sol)")

objv = get_objval(model)
println("objv = $(objv)")

gc()  # test finalizers


text:
dir = "../OCTHaGOn_benchmarks/gams/global/"
filenames = readdir(dir)
GAMS_DIR = dir

include(OCTHaGOn.PROJECT_ROOT * "/test/tools/gams.jl")

pd = DataFrame()
n_max = 100
pd = CSV.read(dir * "problem_stats.csv", DataFrame)
pd = pd[pd.all_bounded .>= 0.5, :] # Use only benchmarks that are bounded...
pd = pd[pd.n_vars .<= n_max, :]    # and have fewer than n_max variables.

function gs_and_time!(m::Union{JuMP.Model, GlobalModel}, n_tries = 10)
    t1 = time()
    if m isa JuMP.Model
        optimize!(m);
    else
        set_param(m, :sample_coeff, 400)
        ct = 0
        while ct < n_tries
            try
                globalsolve!(m)
                ct += 1e5
            catch
                ct += 1
            end
        end
        if ct == n_tries
            throw(OCTHaGOnException("Globalsolve failed for $(m.name)."))
        end
    end
    println("Model solution time: " * string(time()-t1) * " seconds.")
    return
end

# Benchmarking loop
gms = Dict()
results = DataFrame(:name => [], :oct_feas => [], :oct_cost => [], :oct_time => [])
for filename in pd.name
    @info "Trying $(filename)"
    gm = GAMS_to_GlobalModel(dir, filename * ".gms")
    set_optimizer(gm, SOLVER_SILENT)
    gms[gm.name] = gm
    t1 = time()
    try
        gs_and_time!(gm)
        push!(results, [gm.name,  is_feasible(gm), gm.cost[end], time() - t1])
    catch
        if !isempty(gm.cost)
            push!(results, [gm.name,  is_feasible(gm), gm.cost[end], time() - t1])
        else
            push!(results, [gm.name,  false, Inf, time() - t1])
        end
    end
end
results.name = [na[1:end-4] for na in results.name] # Trimming the .gms string

# Updating and rewriting data
orig_data = CSV.read(dir * "problem_stats.csv", DataFrame)
idxs = [idx for idx in 1:size(orig_data, 1) if orig_data.name[idx] in results.name]
orig_data[idxs,[:oct_feas, :oct_cost, :oct_time]] = results[:,[:oct_feas, :oct_cost, :oct_time]]
CSV.write(dir * "problem_stats.csv", orig_data)


text:
using Base.Iterators: partition
using DelimitedFiles
using Flux
using Flux.Data.MNIST
using Flux.Optimise: update!
using Flux: logitbinarycrossentropy, glorot_normal
using Images
using Statistics
using Printf

animation_dir = "mnist-dcgan-animation"
result_dir = "mnist-dcgan-results"

mutable struct DCGAN
    noise_dim::Int64
    channels::Int64
    batch_size::Int64

    generator::Chain
    discriminator::Chain

    generator_optimizer
    discriminator_optimizer

    data::Vector{<: AbstractArray{Float32, 4}}

    animation_size::Pair{Int64, Int64}
    animation_noise::AbstractMatrix{Float32}

    train_steps::Int64
    verbose_freq::Int64
    generator_loss_hist::Vector{Float32}
    discriminator_loss_hist::Vector{Float32}
end

function DCGAN(; image_vector::Vector{<: AbstractMatrix}, noise_dim::Int64, channels::Int64, batch_size::Int64,
    generator::Chain, discriminator::Chain,
    animation_size::Pair{Int64, Int64}, verbose_freq::Int64)

    @assert (channels == 1) || (channels == 3)

    data = [reshape(reduce(hcat, channelview.(xs)), 28, 28, 1, :) for xs in partition(image_vector, batch_size)]
    data = [2f0 .* gpu(Float32.(xs)) .- 1f0 for xs in data]

    animation_noise = randn(Float32, noise_dim, prod(animation_size)) |> gpu

    DCGAN(noise_dim, channels, batch_size, generator, discriminator, ADAM(0.0001f0), ADAM(0.0001f0), data,
        animation_size, animation_noise, 0, verbose_freq, Vector{Float32}(), Vector{Float32}())
end

function generator_loss(fake_output)
    loss = mean(logitbinarycrossentropy.(fake_output, 1f0))
    return loss
end

function discriminator_loss(real_output, fake_output)
    real_loss = mean(logitbinarycrossentropy.(real_output, 1f0))
    fake_loss = mean(logitbinarycrossentropy.(fake_output, 0f0))
    loss = 0.5f0 * (real_loss +  fake_loss)
    return loss
end

function convert_to_image(image_array::Matrix{Float32}, channels::Int64)
    image_array = @. (image_array + 1f0) / 2f0
    if channels == 1
        return Gray.(image_array)
    else
        return colorview(RGB, image_array)
    end
end

function save_fake_image(dcgan::DCGAN)
    @eval Flux.istraining() = false
    fake_images = dcgan.generator(dcgan.animation_noise)
    @eval Flux.istraining() = true
    h, w, _, _ = size(fake_images)
    rows, cols = dcgan.animation_size.first, dcgan.animation_size.second
    tile_image = Matrix{Float32}(undef, h * rows, w * cols)
    for n in 0:prod(dcgan.animation_size) - 1
        j = n ÷ rows
        i = n % cols
        tile_image[j * h + 1:(j + 1) * h, i * w + 1:(i + 1) * w] = fake_images[:, :, :, n + 1] |> cpu
    end
    image = convert_to_image(tile_image, dcgan.channels)
    save(@sprintf("%s/steps_%06d.png", result_dir, dcgan.train_steps), image)
end

function train_discriminator!(dcgan::DCGAN, batch::AbstractArray{Float32, 4})
    noise = randn(Float32, dcgan.noise_dim, dcgan.batch_size) |> gpu
    fake_input = dcgan.generator(noise)
    loss(m) = discriminator_loss(m(batch), m(fake_input))
    disc_grad = gradient(()->loss(dcgan.discriminator), Flux.params(dcgan.discriminator))
    update!(dcgan.discriminator_optimizer, Flux.params(dcgan.discriminator), disc_grad)
    return loss(dcgan.discriminator)
end

function train_generator!(dcgan::DCGAN, batch::AbstractArray{Float32, 4})
    noise = randn(Float32, dcgan.noise_dim, dcgan.batch_size) |> gpu
    loss(m) = generator_loss(dcgan.discriminator(m(noise)))
    gen_grad = gradient(()->loss(dcgan.generator), Flux.params(dcgan.generator))
    update!(dcgan.generator_optimizer, Flux.params(dcgan.generator), gen_grad)
    return loss(dcgan.generator)
end

function train!(dcgan::DCGAN, epochs::Integer)
    for ep in 1:epochs
        @info "epoch $ep"
        for batch in dcgan.data
            disc_loss = train_discriminator!(dcgan, batch)
            gen_loss = train_generator!(dcgan, batch)

            if dcgan.train_steps % dcgan.verbose_freq == 0
                push!(dcgan.discriminator_loss_hist, disc_loss)
                push!(dcgan.generator_loss_hist, gen_loss)
                @info("Train step $(dcgan.train_steps), Discriminator loss: $(disc_loss), Generator loss: $(gen_loss)")
                # create fake images for animation
                save_fake_image(dcgan)
            end
            dcgan.train_steps += 1
        end
    end
end

function main()
    if !isdir(animation_dir)
        mkdir(animation_dir)
    end

    if !isdir(result_dir)
        mkdir(result_dir)
    end

    noise_dim = 100
    channels = 1

    generator = Chain(
        Dense(noise_dim, 7 * 7 * 256; initW = glorot_normal),
        BatchNorm(7 * 7 * 256, relu),
        x->reshape(x, 7, 7, 256, :),
        ConvTranspose((5, 5), 256 => 128; init = glorot_normal, stride = 1, pad = 2),
        BatchNorm(128, relu),
        ConvTranspose((4, 4), 128 => 64; init = glorot_normal, stride = 2, pad = 1),
        BatchNorm(64, relu),
        ConvTranspose((4, 4), 64 => channels, tanh; init = glorot_normal, stride = 2, pad = 1),
        ) |> gpu

    discriminator =  Chain(
        Conv((4, 4), channels => 64, leakyrelu; init = glorot_normal, stride = 2, pad = 1),
        Dropout(0.25),
        Conv((4, 4), 64 => 128, leakyrelu; init = glorot_normal, stride = 2, pad = 1),
        Dropout(0.25),
        x->reshape(x, 7 * 7 * 128, :),
        # drop sigmoid, and use logitbinarycrossentropy (it is more numerically stable)
        # https://github.com/FluxML/Flux.jl/issues/914
        Dense(7 * 7 * 128, 1; initW = glorot_normal)) |> gpu

    dcgan = DCGAN(; image_vector = MNIST.images(), noise_dim = noise_dim,
        channels = channels, batch_size = 128,
        generator = generator, discriminator = discriminator,
        animation_size = 6=>6, verbose_freq = 100)
    train!(dcgan, 30)

    open(@sprintf("%s/discriminator_loss.txt", result_dir), "w") do io
        writedlm(io, dcgan.discriminator_loss_hist)
    end

    open(@sprintf("%s/generator_loss.txt", result_dir), "w") do io
        writedlm(io, dcgan.generator_loss_hist)
    end

end

main()

max_stars_repo_name: matsueushi/fluxjl-dcgan

max_issues_repo_name: matsueushi/fluxjl-gan

max_forks_repo_name: matsueushi/fluxjl-gan


text:
import OpenCL.cl.CLArray

using LinearAlgebra

@testset "OpenCL.CLArray" begin

    @testset "OpenCL.CLArray constructors" begin
        for device in cl.devices()

            ctx = cl.Context(device)
            queue = cl.CmdQueue(ctx)
            hostarray = zeros(Float32, 128*64)
            A = CLArray(queue, hostarray)

            @test CLArray(queue, (:rw, :copy), hostarray) != nothing

            @test CLArray(queue, hostarray, flags=(:rw, :copy)) != nothing

            @test CLArray(queue, hostarray) != nothing

            @test CLArray(cl.Buffer(Float32, ctx, (:r, :copy), hostbuf=hostarray),
                          queue,
                          (128, 64)) != nothing

            @test copy(A) == A
        end
     end

    @testset "OpenCL.CLArray fill" begin
        for device in cl.devices()
            ctx = cl.Context(device)
            queue = cl.CmdQueue(ctx)

            @test cl.to_host(cl.fill(Float32, queue, Float32(0.5),
                                            32, 64)) == fill(Float32(0.5), 32, 64)
            @test cl.to_host(cl.zeros(Float32, queue, 64)) == zeros(Float32, 64)
            @test cl.to_host(cl.ones(Float32, queue, 64)) == ones(Float32, 64)
        end
     end

    @testset "OpenCL.CLArray core functions" begin
        for device in cl.devices()
            ctx = cl.Context(device)
            queue = cl.CmdQueue(ctx)
            A = CLArray(queue, rand(Float32, 128, 64))
            @test size(A) == (128, 64)
            @test ndims(A) == 2
            @test length(A) == 128*64
            # reshape
            B = reshape(A, 128*64)
            @test reshape(B, 128, 64) == A
            # transpose
            X = CLArray(queue, rand(Float32, 32, 32))
            B = cl.zeros(Float32, queue, 64, 128)
            ev = transpose!(B, A)
            cl.wait(ev)
            #@test cl.to_host(copy(A')) == cl.to_host(B)
        end
     end
end


text:
"""
CUR (columns-U-rows) 3-tensor decomposition.
"""
immutable CUR <: TensorDecomposition{Float64, 3}
    Cindex::Vector{Int64}
    Cweight::Vector{Int64}
    Rindex::Vector{NTuple{2, Int64}}
    Rweight::Vector{Int64}
    U::Matrix{Float64}
    error::Vector{Float64}

    function CUR(tnsr::StridedArray{Float64},
                 slab_axis::Integer,
                 fiber_axes::NTuple{2, Int64},
                 fiber_size::NTuple{2, Int64},
                 Cindex::Vector{Int64},
                 Cweight::Vector{Int64},
                 Rindex::Vector{Int64},
                 Rweight::Vector{Int64},
                 U::Matrix{Float64},
                 compute_u::Bool)

        res = zeros(0)
        if compute_u
            output_index = collect(1:3)
            output_index[slab_axis] = 4
            S = tensorcontract(slicedim(tnsr, slab_axis, Cindex), [1, 2, 3], U, [4, slab_axis], output_index)
            W = squeeze(permutedims(mapslices(slab -> slab[Rindex], tnsr, fiber_axes),
                                    [slab_axis, fiber_axes[1], fiber_axes[2]]), 3)
            S = tensorcontract(S, output_index, W, [slab_axis, 4], [1, 2, 3])
            res = mapslices(vecnorm, S - tnsr, fiber_axes)[:] ./ mapslices(vecnorm, tnsr, fiber_axes)[:]
        end
        new(Cindex, Cweight, [zip(ind2sub(fiber_size, Rindex)...)...], Rweight, U, res)
    end
end

"""
Calculates CUR decomposition for 3-mode tensors.
"""
function tensorcur3(tnsr::StridedArray,
                    c::Integer, r::Integer,
                    slab_axis::Integer=3;
                    compute_u::Bool=true)

    ndims(tnsr) == 3 || error("This method currently only supports 3-mode tensors.")
    slab_axis > 0 && slab_axis < 4 || error("Invalid slab_axis; slab_axis should be 1, 2, or 3")
    fiber_axes = tuple(circshift([1, 2, 3], 1 - slab_axis)[2:3]...)
    fiber_size = (size(tnsr, fiber_axes[1]), size(tnsr, fiber_axes[2]))
    T2 = tnsr .^ 2
    T2_sum = sum(T2)
    p = sum(T2, fiber_axes)[:] / T2_sum
    q = sum(T2, slab_axis)[:] / T2_sum

    Cindex = rand(Categorical(p), c)
    Rindex = rand(Categorical(q), r)
    Cweight = hist(Cindex, 0:size(tnsr, slab_axis))[2]
    Rweight = hist(Rindex, 0:prod(fiber_size))[2]
    Cindex = sort(unique(Cindex))
    Rindex = sort(unique(Rindex))
    Cweight = Cweight[Cindex]
    Rweight = Rweight[Rindex]

    U = compute_u ? Array(Float64, r, c) : zeros(0, 0)
    if compute_u
        P = (Cweight*Rweight') ./ (p[Cindex]*q[Rindex]')
        W = squeeze(permutedims(mapslices(slab -> slab[Rindex],
                                slicedim(tnsr, slab_axis, Cindex), fiber_axes),
                                [slab_axis, fiber_axes[1], fiber_axes[2]]), 3)
        U = pinv(W .* P) .* P'
    end

    return CUR(tnsr, slab_axis, fiber_axes, fiber_size, Cindex, Cweight, Rindex, Rweight, U, compute_u)
end


text:
export pdesolve

include("cont_lyap.jl")

# Converts a vector of constants/funs to a vector of Funs
function convert2fun{T<:Number,S<:Space}(f::Array{T},sp::S)
    ret=similar(f,Fun{S,T})
    for k=1:size(f,1),j=1:size(f,2)
        ret[k,j]=Fun(f[k,j],sp)
    end
    ret
end
convert2fun{T<:Fun}(f::Array{T},d::Space)=f
function convert2fun{S<:Space}(f::Array,sp::S)
    mytyp=Fun{S,promote_type(mapreduce(eltype,promote_type,f),eltype(sp))}

    ret=similar(f,mytyp)

    for k=1:size(f,1),j=1:size(f,2)
        #TODO: Check domains match for Funs
        ret[k,j]=Fun(f[k,j],sp)
    end
    ret
end

###
# pde_standardize_rhs manipulate "non-standard" rhs
# so it has a standardized form
###

function depiecereorient(f,d)
    vf=pieces(f)
    dd1=d[1];dd2=d[2]

    if length(vf)==4
        # assume we are in boudnary
        # boundary has positive orientation
        vf=Any[setdomain(reverseorientation(vf[4]),dd2),
               setdomain(vf[2],dd2),
               setdomain(vf[1],dd1),
               setdomain(reverseorientation(vf[3]),dd1)]
    else
        @assert length(vf)==2
        if isa(dd1,PeriodicDomain)
            vf=Any[setdomain(vf[1],dd1),setdomain(reverseorientation(vf[2]),dd1)]
        elseif isa(dd2,PeriodicDomain)
            vf=Any[setdomain(vf[1],dd2),setdomain(reverseorientation(vf[2]),dd2)]
        end
    end
    vf
end

function pde_standardize_rhs(A,f::Vector)
    indsBx=bcinds(A,1);indsBy=bcinds(A,2)

    rs=rangespace(A)

    # if the fun lives on the boundary, we need to get rid of
    # the boundary information
    # i.e., if it lives on (-1-im,-1+im) we want to convert it to
    # living on (-1,1)
    if isa(A,PDEOperatorSchur) &&
            !isempty(f) &&
            isa(f[1],Fun) &&
            domain(f[1])==∂(domain(A))
        @assert length(f)<=2
        @assert isa(f,Vector)

        #TODO: More elegent domain conversion
        vf=depiecereorient(f[1],domain(A))

        if length(f)==1
            f=vf
        else
            @assert length(f)==2
            f=Any[vf...;f[end]]
        end
    end
    f=pad(f,max(length(indsBx)+length(indsBy),length(f)))

    ds=domainspace(A)
    fx=isempty(indsBx)?Fun{spacetype(ds,2),eltype(A)}[]:convert2fun(f[indsBx],ds[2])
    fy=isempty(indsBy)?Fun{spacetype(ds,1),eltype(A)}[]:convert2fun(f[indsBy],ds[1])

    if length(f)==length(indsBx)+length(indsBy)+1
        F=Fun(f[end],rs)
    else
        F=zeros(eltype(A),rs)
    end

    fx,fy,F
end

function pde_standardize_rhs(A,f::Matrix)
    ∂A=∂(domain(A))
    indsBx=bcinds(A,1);indsBy=bcinds(A,2)

    if isa(A,PDEOperatorSchur) &&
                !isempty(f) &&
                all(g->isa(g,Fun)&&domain(g)==∂A,f[1,:])
        @assert size(f,1)==1
        d=domain(A)
        f2=Array(Any,length(indsBx)+length(indsBy),size(f,2))
        for k=1:size(f,2)
            f2[:,k]=depiecereorient(f[1,k],d)
        end
        f=f2
    end

    @assert length(indsBx)+length(indsBy)==size(f,1)

    fx=isempty(indsBx)?[]:convert2fun(f[indsBx,:],domainspace(A,2))
    fy=isempty(indsBy)?[]:convert2fun(f[indsBy,:],domainspace(A,1))
    F=fill(zeros(rangespace(A)),1,size(f,2))

    fx,fy,F
end

function pdesolve{T<:Operator}(A::Vector{T},f;tolerance::Float64=100000eps(),opts...)
    @assert tolerance>0
    maxit=11

    for k=5:maxit
        u=pdesolve(A,f,2^k;tolerance=tolerance/100,opts...)
        if norm(map(f->norm(f.coefficients),u.coefficients[end-2:end]))<tolerance
            return u
        end
    end
    error("Maximum number of iterations " * string(maxit) * "reached")
end

function pdesolve(A::AbstractPDEOperatorSchur,f::Vector;opts...)
    fx,fy,F=pde_standardize_rhs(A,f)
    ProductFun(cont_constrained_lyap(A,fx,fy,F;opts...),domainspace(A))
end

function pdesolve(A::AbstractPDEOperatorSchur,f::Matrix;opts...)
    fx,fy,F=pde_standardize_rhs(A,f)
    X=cont_constrained_lyap(A,fx,fy,F;opts...)
    ds=domainspace(A)
    #TODO: Change Float64
    ProductFun{typeof(ds[1]),
               typeof(ds[2]),
               typeof(ds),
               Float64}[ProductFun(X[:,k],ds) for k=1:size(X,2)]
end

pdesolve(A::AbstractPDEOperatorSchur,f::Union{Fun,MultivariateFun,Number},nx...;opts...)=pdesolve(A,[f],nx...;opts...)
pdesolve{T<:Operator}(A::Vector{T},f::Vector,n::Integer,n2...;opts...)=pdesolve(schurfact(A,n),f,n2...;opts...)
pdesolve{T<:Operator}(A::Vector{T},f::Union{Fun,MultivariateFun,Number},n...;opts...)=pdesolve(A,[f],n...;opts...)
pdesolve(A::Operator,f...;opts...)=pdesolve([A],f...;opts...)

\(A::AbstractPDEOperatorSchur,f::Union{MultivariateFun,Number,Fun,Array})=pdesolve(A,f)


text:
"""
    abstract type AbstractSplitter

An `AbstractSplitter` is used to split a [`ClusterTree`](@ref). The interface
requires the following methods:
- `should_split(clt,splitter)` : return a `Bool` determining if the
  `ClusterTree` should be further divided
- `split!(clt,splitter)` : perform the splitting of the `ClusterTree` handling
  the necessary data sorting.

See [`GeometricSplitter`](@ref) for an example of an implementation.
"""
abstract type AbstractSplitter end

"""
    should_split(clt::ClusterTree,splitter::AbstractSplitter)

Determine whether or not a `ClusterTree` should be further divided.
"""
function should_split(clt, splitter::AbstractSplitter)
    abstract_method(splitter)
end

"""
    split!(clt::ClusterTree,splitter::AbstractSplitter)

Divide `clt` using the strategy implemented by `splitter`.
"""
function split!(clt, splitter::AbstractSplitter)
    abstract_method(splitter)
end

"""
    struct DyadicSplitter

Used to split an `N` dimensional `ClusterTree` into `2^N` children until at most
`nmax` points are contained in node.
"""
Base.@kwdef struct DyadicSplitter <: AbstractSplitter
    nmax::Int = typemax(Int)
end

function should_split(node::ClusterTree, splitter::DyadicSplitter)
    return length(node) > splitter.nmax
end

function split!(parentcluster::ClusterTree, ::DyadicSplitter)
    d = ambient_dimension(parentcluster)
    clusters = [parentcluster]
    rec = container(parentcluster)
    rec_center = center(rec)
    for i = 1:d
        pos = rec_center[i]
        nel = length(clusters) #2^(i-1)
        for _ = 1:nel
            clt = popfirst!(clusters)
            append!(clusters, _binary_split!(clt, i, pos; parentcluster))
        end
    end
    return clusters
end

"""
    struct GeometricSplitter <: AbstractSplitter

Used to split a `ClusterTree` in half along the largest axis.
"""
Base.@kwdef struct GeometricSplitter <: AbstractSplitter
    nmax::Int = 50
end

should_split(node::ClusterTree, splitter::GeometricSplitter) = length(node) > splitter.nmax

function split!(cluster::ClusterTree, ::GeometricSplitter)
    rec = cluster.container
    wmax, imax = findmax(high_corner(rec) - low_corner(rec))
    left_node, right_node = _binary_split!(cluster, imax, low_corner(rec)[imax] + wmax / 2)
    return [left_node, right_node]
end

"""
    struct GeometricMinimalSplitter <: AbstractSplitter

Like [`GeometricSplitter`](@ref), but shrinks the children's containters.
"""
Base.@kwdef struct GeometricMinimalSplitter <: AbstractSplitter
    nmax::Int = 50
end

should_split(node::ClusterTree, splitter::GeometricMinimalSplitter) = length(node) > splitter.nmax

function split!(cluster::ClusterTree, ::GeometricMinimalSplitter)
    rec = cluster.container
    wmax, imax = findmax(high_corner(rec) - low_corner(rec))
    mid = low_corner(rec)[imax] + wmax / 2
    predicate = (x) -> x[imax] < mid
    left_node, right_node = _binary_split!(cluster, predicate)
    return [left_node, right_node]
end

"""
    struct PrincipalComponentSplitter <: AbstractSplitter
"""
Base.@kwdef struct PrincipalComponentSplitter <: AbstractSplitter
    nmax::Int = 50
end

should_split(node::ClusterTree, splitter::PrincipalComponentSplitter) = length(node) > splitter.nmax

function split!(cluster::ClusterTree, ::PrincipalComponentSplitter)
    pts = cluster._elements
    irange = cluster.index_range
    xc = center_of_mass(cluster)
    # compute covariance matrix for principal direction
    cov = sum(irange) do i
        x = coords(pts[i])
        (x - xc) * transpose(x - xc)
    end
    v = eigvecs(cov)[:, end]
    predicate = (x) -> dot(x - xc, v) < 0
    left_node, right_node = _binary_split!(cluster, predicate)
    return [left_node, right_node]
end

function center_of_mass(clt::ClusterTree)
    pts = clt._elements
    loc_idxs = clt.index_range
    # w    = clt.weights
    n = length(loc_idxs)
    # M    = isempty(w) ? n : sum(i->w[i],glob_idxs)
    # xc   = isempty(w) ? sum(i->pts[i]/M,glob_idxs) : sum(i->w[i]*pts[i]/M,glob_idxs)
    M = n
    xc = sum(i -> coords(pts[i]) / M, loc_idxs)
    return xc
end

"""
    struct CardinalitySplitter <: AbstractSplitter

Used to split a `ClusterTree` along the largest dimension if
`length(tree)>nmax`. The split is performed so the `data` is evenly distributed
amongst all children.
"""
Base.@kwdef struct CardinalitySplitter <: AbstractSplitter
    nmax::Int = 50
end

should_split(node::ClusterTree, splitter::CardinalitySplitter) = length(node) > splitter.nmax

function split!(cluster::ClusterTree, ::CardinalitySplitter)
    points = cluster._elements
    irange = cluster.index_range
    rec = container(cluster)
    _, imax = findmax(high_corner(rec) - low_corner(rec))
    med = median(coords(points[i])[imax] for i in irange) # the median along largest axis `imax`
    predicate = (x) -> x[imax] < med
    left_node, right_node = _binary_split!(cluster, predicate)
    return [left_node, right_node]
end

max_stars_repo_name: WaveProp/WavePropBase

max_issues_repo_name: WaveProp/WavePropBase

max_forks_repo_name: WaveProp/WavePropBase


text:
# Autogenerated wrapper script for Capstone_jll for x86_64-apple-darwin
export cstool, libcapstone

JLLWrappers.@generate_wrapper_header("Capstone")
JLLWrappers.@declare_executable_product(cstool)
JLLWrappers.@declare_library_product(libcapstone, "@rpath/libcapstone.4.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        cstool,
        "bin/cstool",
    )

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

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module nj

macro vers03x(ex)
   VERSION.minor == 3
end

macro vers03x_only(ex)
   @vers03x(ex)?esc(ex):nothing
end

macro vers04x(ex)
   VERSION.minor == 4
end

macro vers04x_only(ex)
   @vers04x(ex)?esc(ex):nothing
end

preserve = Array(Any,0);

function topExpr(mod::Module,paths::Array{ASCIIString,1})
   res = Expr(:toplevel,:(eval(x) = Core.eval(mod,x)),:(eval(m,x) = Core.eval(m,x)))
   for path in paths
      push!(res.args,:(include($path)))
   end
   return res;
end

include(mod::Module,path,args::Vector) = include(mod,path,UTF8String[args...])

function scriptify(mod::Module,filename::ASCIIString)
    ast = parse("function _(args...)\n" * readall(filename) * "end");
    args2 = Array(Any,0);
    paths = Array(ASCIIString,0);

    for astNode in ast.args[2].args
       if(typeof(astNode) == Expr && astNode.head == :call && astNode.args[1] == :include)
          push!(paths,astNode.args[2]);
       else
          push!(args2,astNode);
       end
    end
    res = topExpr(mod,paths);
    push!(res.args,Expr(ast.head,ast.args[1],Expr(ast.args[2].head,args2...)));
    return res;
end

newRegex(pattern) = Regex(pattern)
getPattern(re::Regex) = re.pattern
getRegexType() = Regex

@vers03x_only getDateTimeType() = typeof(nothing)
@vers04x_only getDateTimeType() = DateTime
@vers03x_only toDate(milliseconds) = nothing

@vers04x_only function toDate(milliseconds::Int64)
   return DateTime(1970) + Base.Dates.Millisecond(milliseconds)
end

@vers03x_only toMilliseconds(date) = nothing
@vers04x_only toMilliseconds(date) = Base.Dates.datetime2unix(date);

function getError(ex,bt)
   io = IOBuffer();
   Base.showerror(io,ex,bt);
   s = takebuf_string(io);
   return s;
end

#  The following copies the base definition of require and the supporting functions
#  but redefines pwd() to avoid the conflict caused by node's and julia's double
#  definition of uv_cwd which are not compatible.
#
#  This can be removed once the issue is resolved in base.

@vers04x_only function pwd()
    b = Array(UInt8,1024)
    len = Csize_t[length(b),]
    Base.uv_error(:getcwd, ccall((:uv_cwd,"libjulia"), Cint, (Ptr{UInt8}, Ptr{Csize_t}), b, len))
    bytestring(b[1:len[1]-1])
end

@vers03x_only function pwd()
    b = Array(Uint8,1024)
    len = Csize_t[length(b),]
    Base.uv_error(:getcwd, ccall((:uv_cwd,"libjulia"), Cint, (Ptr{Uint8}, Ptr{Csize_t}), b, len))
    bytestring(b[1:len[1]-1])
end

abspath(a::String) = normpath(isabspath(a) ? a : joinpath(pwd(),a))
abspath(a::String, b::String...) = abspath(joinpath(a,b...))

function find_in_path(name::String)
    isabspath(name) && return name
    isfile(name) && return abspath(name)
    base = name
    if endswith(name,".jl")
        base = name[1:end-3]
    else
        name = string(base,".jl")
        isfile(name) && return abspath(name)
    end
    for prefix in [Pkg.dir(); LOAD_PATH]
        path = joinpath(prefix, name)
        isfile(path) && return abspath(path)
        path = joinpath(prefix, base, "src", name)
        isfile(path) && return abspath(path)
        path = joinpath(prefix, name, "src", name)
        isfile(path) && return abspath(path)
    end
    return nothing
end

find_in_node1_path(name) = myid()==1 ?
    find_in_path(name) : remotecall_fetch(1, find_in_path, name)

function find_source_file(file)
    (isabspath(file) || isfile(file)) && return file
    file2 = find_in_path(file)
    file2 != nothing && return file2
    file2 = joinpath(JULIA_HOME, DATAROOTDIR, "julia", "base", file)
    isfile(file2) ? file2 : nothing
end

# Store list of files and their load time
package_list = Dict{ByteString,Float64}()
# to synchronize multiple tasks trying to require something
package_locks = Dict{ByteString,Any}()
require(fname::String) = require(bytestring(fname))
require(f::String, fs::String...) = (require(f); for x in fs require(x); end)

# only broadcast top-level (not nested) requires and reloads
toplevel_load = true

function require(name::String)
    path = find_in_node1_path(name)
    path == nothing && throw(ArgumentError("$name not found in path"))

    if myid() == 1 && toplevel_load
        refs = Any[ @spawnat p _require(path) for p in filter(x->x!=1, procs()) ]
        _require(path)
        for r in refs; wait(r); end
    else
        _require(path)
    end
    nothing
end

function _require(path)
    global toplevel_load
    if haskey(package_list,path)
        wait(package_locks[path])
    else
        last = toplevel_load
        toplevel_load = false
        try
            reload_path(path)
        finally
            toplevel_load = last
        end
    end
end

function reload(name::String)
    global toplevel_load
    path = find_in_node1_path(name)
    path == nothing && throw(ArgumentError("$name not found in path"))
    refs = nothing
    if myid() == 1 && toplevel_load
        refs = Any[ @spawnat p reload_path(path) for p in filter(x->x!=1, procs()) ]
    end
    last = toplevel_load
    toplevel_load = false
    try
        reload_path(path)
    finally
        toplevel_load = last
    end
    if refs !== nothing
        for r in refs; wait(r); end
    end
    nothing
end

# remote/parallel load

function source_path(default::Union(String,Void)="")
    t = current_task()
    while true
        s = t.storage
        if !is(s, nothing) && haskey(s, :SOURCE_PATH)
            return s[:SOURCE_PATH]
        end
        if is(t, t.parent)
            return default
        end
        t = t.parent
    end
end

macro __FILE__() source_path() end

function include_from_node1(path::String)
    prev = source_path(nothing)
    path = (prev == nothing) ? abspath(path) : joinpath(dirname(prev),path)
    tls = task_local_storage()
    tls[:SOURCE_PATH] = path
    local result
    try
        if myid()==1
            # sleep a bit to process file requests from other nodes
            nprocs()>1 && sleep(0.005)
            result = Core.include(path)
            nprocs()>1 && sleep(0.005)
        else
            result = include_string(remotecall_fetch(1, readall, path), path)
        end
    finally
        if prev == nothing
            delete!(tls, :SOURCE_PATH)
        else
            tls[:SOURCE_PATH] = prev
        end
    end
    result
end

function reload_path(path::String)
    had = haskey(package_list, path)
    if !had
        package_locks[path] = RemoteRef()
    end
    package_list[path] = time()
    tls = task_local_storage()
    prev = pop!(tls, :SOURCE_PATH, nothing)
    try
        eval(Main, :(Base.include_from_node1($path)))
    catch e
        had || delete!(package_list, path)
        rethrow(e)
    finally
        if prev != nothing
            tls[:SOURCE_PATH] = prev
        end
    end
    if !isready(package_locks[path])
        put!(package_locks[path],nothing)
    end
    nothing
end

#
# End loading.jl inclusion

function importModule(modulePath::ASCIIString)
   substitutedModulePath= reduce((x,y)->"$x$(Base.path_separator)$y",split(modulePath,"/"))
   require(substitutedModulePath);

   functionNames::Array{UTF8String,1} = Array(UTF8String,0);
   m::Module = Base.eval(Main,symbol(split(substitutedModulePath,Base.path_separator)[end]));

   for name in names(m)
      push!(functionNames,"$name");
   end
   return m,functionNames;
end

newTuple(args...)=tuple(args...)

end

max_stars_repo_name: Alfheim/node-julia

max_issues_repo_name: Alfheim/node-julia

max_forks_repo_name: Alfheim/node-julia


text:
module AMD

import GetC.@getCFun

typealias GLenum Cuint
typealias GLboolean Cuchar
typealias GLbitfield Cuint
typealias GLvoid Void
typealias GLbyte Cuchar
typealias GLshort Cshort
typealias GLint Cint
typealias GLubyte Cuchar
typealias GLushort Cushort
typealias GLuint Cuint
typealias GLsizei Csize_t
typealias GLfloat Cfloat
typealias GLclampf Cfloat
typealias GLdouble Cdouble
typealias GLclampd Cdouble
typealias GLchar Cchar
typealias GLint64 Clonglong
typealias GLuint64 Culonglong
typealias GLhalf Cushort
typealias GLhalfARB Cushort
typealias GLhalfNV Cushort
typealias GLsync Ptr{Void}
typealias Pointer Ptr{Void}
typealias GLsizeiptr Cint
typealias GLintptr Cptrdiff_t
const FACTOR_MAX_AMD = 0x901D
export FACTOR_MAX_AMD
const FACTOR_MIN_AMD = 0x901C
export FACTOR_MIN_AMD
const DEBUG_CATEGORY_API_ERROR_AMD = 0x9149
export DEBUG_CATEGORY_API_ERROR_AMD
const DEBUG_CATEGORY_APPLICATION_AMD = 0x914F
export DEBUG_CATEGORY_APPLICATION_AMD
const DEBUG_CATEGORY_DEPRECATION_AMD = 0x914B
export DEBUG_CATEGORY_DEPRECATION_AMD
const DEBUG_CATEGORY_OTHER_AMD = 0x9150
export DEBUG_CATEGORY_OTHER_AMD
const DEBUG_CATEGORY_PERFORMANCE_AMD = 0x914D
export DEBUG_CATEGORY_PERFORMANCE_AMD
const DEBUG_CATEGORY_SHADER_COMPILER_AMD = 0x914E
export DEBUG_CATEGORY_SHADER_COMPILER_AMD
const DEBUG_CATEGORY_UNDEFINED_BEHAVIOR_AMD = 0x914C
export DEBUG_CATEGORY_UNDEFINED_BEHAVIOR_AMD
const DEBUG_CATEGORY_WINDOW_SYSTEM_AMD = 0x914A
export DEBUG_CATEGORY_WINDOW_SYSTEM_AMD
const DEBUG_LOGGED_MESSAGES_AMD = 0x9145
export DEBUG_LOGGED_MESSAGES_AMD
const DEBUG_SEVERITY_HIGH_AMD = 0x9146
export DEBUG_SEVERITY_HIGH_AMD
const DEBUG_SEVERITY_LOW_AMD = 0x9148
export DEBUG_SEVERITY_LOW_AMD
const DEBUG_SEVERITY_MEDIUM_AMD = 0x9147
export DEBUG_SEVERITY_MEDIUM_AMD
const MAX_DEBUG_LOGGED_MESSAGES_AMD = 0x9144
export MAX_DEBUG_LOGGED_MESSAGES_AMD
const MAX_DEBUG_MESSAGE_LENGTH_AMD = 0x9143
export MAX_DEBUG_MESSAGE_LENGTH_AMD
const DEPTH_CLAMP_FAR_AMD = 0x901F
export DEPTH_CLAMP_FAR_AMD
const DEPTH_CLAMP_NEAR_AMD = 0x901E
export DEPTH_CLAMP_NEAR_AMD
const DATA_BUFFER_AMD = 0x9151
export DATA_BUFFER_AMD
const PERFORMANCE_MONITOR_AMD = 0x9152
export PERFORMANCE_MONITOR_AMD
const QUERY_OBJECT_AMD = 0x9153
export QUERY_OBJECT_AMD
const SAMPLER_OBJECT_AMD = 0x9155
export SAMPLER_OBJECT_AMD
const VERTEX_ARRAY_OBJECT_AMD = 0x9154
export VERTEX_ARRAY_OBJECT_AMD
const COUNTER_RANGE_AMD = 0x8BC1
export COUNTER_RANGE_AMD
const COUNTER_TYPE_AMD = 0x8BC0
export COUNTER_TYPE_AMD
const PERCENTAGE_AMD = 0x8BC3
export PERCENTAGE_AMD
const PERFMON_RESULT_AMD = 0x8BC6
export PERFMON_RESULT_AMD
const PERFMON_RESULT_AVAILABLE_AMD = 0x8BC4
export PERFMON_RESULT_AVAILABLE_AMD
const PERFMON_RESULT_SIZE_AMD = 0x8BC5
export PERFMON_RESULT_SIZE_AMD
const UNSIGNED_INT64_AMD = 0x8BC2
export UNSIGNED_INT64_AMD
const EXTERNAL_VIRTUAL_MEMORY_BUFFER_AMD = 0x9160
export EXTERNAL_VIRTUAL_MEMORY_BUFFER_AMD
const QUERY_BUFFER_AMD = 0x9192
export QUERY_BUFFER_AMD
const QUERY_BUFFER_BINDING_AMD = 0x9193
export QUERY_BUFFER_BINDING_AMD
const QUERY_RESULT_NO_WAIT_AMD = 0x9194
export QUERY_RESULT_NO_WAIT_AMD
const SUBSAMPLE_DISTANCE_AMD = 0x883F
export SUBSAMPLE_DISTANCE_AMD
const TEXTURE_CUBE_MAP_SEAMLESS = 0x884F
export TEXTURE_CUBE_MAP_SEAMLESS
const MAX_SPARSE_3D_TEXTURE_SIZE_AMD = 0x9199
export MAX_SPARSE_3D_TEXTURE_SIZE_AMD
const MAX_SPARSE_ARRAY_TEXTURE_LAYERS = 0x919A
export MAX_SPARSE_ARRAY_TEXTURE_LAYERS
const MAX_SPARSE_TEXTURE_SIZE_AMD = 0x9198
export MAX_SPARSE_TEXTURE_SIZE_AMD
const MIN_LOD_WARNING_AMD = 0x919C
export MIN_LOD_WARNING_AMD
const MIN_SPARSE_LEVEL_AMD = 0x919B
export MIN_SPARSE_LEVEL_AMD
const TEXTURE_STORAGE_SPARSE_BIT_AMD = 0x00000001
export TEXTURE_STORAGE_SPARSE_BIT_AMD
const VIRTUAL_PAGE_SIZE_X_AMD = 0x9195
export VIRTUAL_PAGE_SIZE_X_AMD
const VIRTUAL_PAGE_SIZE_Y_AMD = 0x9196
export VIRTUAL_PAGE_SIZE_Y_AMD
const VIRTUAL_PAGE_SIZE_Z_AMD = 0x9197
export VIRTUAL_PAGE_SIZE_Z_AMD
const REPLACE_VALUE_AMD = 0x874B
export REPLACE_VALUE_AMD
const SET_AMD = 0x874A
export SET_AMD
const STENCIL_BACK_OP_VALUE_AMD = 0x874D
export STENCIL_BACK_OP_VALUE_AMD
const STENCIL_OP_VALUE_AMD = 0x874C
export STENCIL_OP_VALUE_AMD
const CONTINUOUS_AMD = 0x9007
export CONTINUOUS_AMD
const DISCRETE_AMD = 0x9006
export DISCRETE_AMD
const INT_SAMPLER_BUFFER_AMD = 0x9002
export INT_SAMPLER_BUFFER_AMD
const SAMPLER_BUFFER_AMD = 0x9001
export SAMPLER_BUFFER_AMD
const TESSELLATION_FACTOR_AMD = 0x9005
export TESSELLATION_FACTOR_AMD
const TESSELLATION_MODE_AMD = 0x9004
export TESSELLATION_MODE_AMD
const UNSIGNED_INT_SAMPLER_BUFFER_AMD = 0x9003
export UNSIGNED_INT_SAMPLER_BUFFER_AMD
@getCFun "libGL" glDebugMessageEnableAMD glDebugMessageEnableAMD(category::GLenum, severity::GLenum, count::GLsizei, ids::Ptr{GLuint}, enabled::GLboolean)::Void
export glDebugMessageEnableAMD
@getCFun "libGL" glDebugMessageInsertAMD glDebugMessageInsertAMD(category::GLenum, severity::GLenum, id::GLuint, length::GLsizei, buf::Ptr{GLchar})::Void
export glDebugMessageInsertAMD
@getCFun "libGL" glDebugMessageCallbackAMD glDebugMessageCallbackAMD(callback::Ptr{Void}, userParam::Ptr{Void})::Void
export glDebugMessageCallbackAMD
@getCFun "libGL" glGetDebugMessageLogAMD glGetDebugMessageLogAMD(count::GLuint, bufsize::GLsizei, categories::Ptr{GLenum}, severities::Ptr{GLuint}, ids::Ptr{GLuint}, lengths::Ptr{GLsizei}, message::Ptr{GLchar})::Cuint
export glGetDebugMessageLogAMD
@getCFun "libGL" glBlendFuncIndexedAMD glBlendFuncIndexedAMD(buf::GLuint, src::GLenum, dst::GLenum)::Void
export glBlendFuncIndexedAMD
@getCFun "libGL" glBlendFuncSeparateIndexedAMD glBlendFuncSeparateIndexedAMD(buf::GLuint, srcRGB::GLenum, dstRGB::GLenum, srcAlpha::GLenum, dstAlpha::GLenum)::Void
export glBlendFuncSeparateIndexedAMD
@getCFun "libGL" glBlendEquationIndexedAMD glBlendEquationIndexedAMD(buf::GLuint, mode::GLenum)::Void
export glBlendEquationIndexedAMD
@getCFun "libGL" glBlendEquationSeparateIndexedAMD glBlendEquationSeparateIndexedAMD(buf::GLuint, modeRGB::GLenum, modeAlpha::GLenum)::Void
export glBlendEquationSeparateIndexedAMD
@getCFun "libGL" glMultiDrawArraysIndirectAMD glMultiDrawArraysIndirectAMD(mode::GLenum, indirect::Ptr{Void}, primcount::GLsizei, stride::GLsizei)::Void
export glMultiDrawArraysIndirectAMD
@getCFun "libGL" glMultiDrawElementsIndirectAMD glMultiDrawElementsIndirectAMD(mode::GLenum, type_::GLenum, indirect::Ptr{Void}, primcount::GLsizei, stride::GLsizei)::Void
export glMultiDrawElementsIndirectAMD
@getCFun "libGL" glGenNamesAMD glGenNamesAMD(identifier::GLenum, num::GLuint, names::Ptr{GLuint})::Void
export glGenNamesAMD
@getCFun "libGL" glDeleteNamesAMD glDeleteNamesAMD(identifier::GLenum, num::GLuint, names::Ptr{GLuint})::Void
export glDeleteNamesAMD
@getCFun "libGL" glIsNameAMD glIsNameAMD(identifier::GLenum, name::GLuint)::Bool
export glIsNameAMD
@getCFun "libGL" glGetPerfMonitorGroupsAMD glGetPerfMonitorGroupsAMD(numGroups::Ptr{GLint}, groupsSize::GLsizei, groups::Ptr{GLuint})::Void
export glGetPerfMonitorGroupsAMD
@getCFun "libGL" glGetPerfMonitorCountersAMD glGetPerfMonitorCountersAMD(group::GLuint, numCounters::Ptr{GLint}, maxActiveCounters::Ptr{GLint}, counterSize::GLsizei, counters::Ptr{GLuint})::Void
export glGetPerfMonitorCountersAMD
@getCFun "libGL" glGetPerfMonitorGroupStringAMD glGetPerfMonitorGroupStringAMD(group::GLuint, bufSize::GLsizei, length::Ptr{GLsizei}, groupString::Ptr{GLchar})::Void
export glGetPerfMonitorGroupStringAMD
@getCFun "libGL" glGetPerfMonitorCounterStringAMD glGetPerfMonitorCounterStringAMD(group::GLuint, counter::GLuint, bufSize::GLsizei, length::Ptr{GLsizei}, counterString::Ptr{GLchar})::Void
export glGetPerfMonitorCounterStringAMD
@getCFun "libGL" glGetPerfMonitorCounterInfoAMD glGetPerfMonitorCounterInfoAMD(group::GLuint, counter::GLuint, pname::GLenum, data::Ptr{Void})::Void
export glGetPerfMonitorCounterInfoAMD
@getCFun "libGL" glGenPerfMonitorsAMD glGenPerfMonitorsAMD(n::GLsizei, monitors::Ptr{GLuint})::Void
export glGenPerfMonitorsAMD
@getCFun "libGL" glDeletePerfMonitorsAMD glDeletePerfMonitorsAMD(n::GLsizei, monitors::Ptr{GLuint})::Void
export glDeletePerfMonitorsAMD
@getCFun "libGL" glSelectPerfMonitorCountersAMD glSelectPerfMonitorCountersAMD(monitor::GLuint, enable::GLboolean, group::GLuint, numCounters::GLint, counterList::Ptr{GLuint})::Void
export glSelectPerfMonitorCountersAMD
@getCFun "libGL" glBeginPerfMonitorAMD glBeginPerfMonitorAMD(monitor::GLuint)::Void
export glBeginPerfMonitorAMD
@getCFun "libGL" glEndPerfMonitorAMD glEndPerfMonitorAMD(monitor::GLuint)::Void
export glEndPerfMonitorAMD
@getCFun "libGL" glGetPerfMonitorCounterDataAMD glGetPerfMonitorCounterDataAMD(monitor::GLuint, pname::GLenum, dataSize::GLsizei, data::Ptr{GLuint}, bytesWritten::Ptr{GLint})::Void
export glGetPerfMonitorCounterDataAMD
@getCFun "libGL" glSetMultisamplefvAMD glSetMultisamplefvAMD(pname::GLenum, index::GLuint, val::Ptr{GLfloat})::Void
export glSetMultisamplefvAMD
@getCFun "libGL" glTexStorageSparseAMD glTexStorageSparseAMD(target::GLenum, internalFormat::GLenum, width::GLsizei, height::GLsizei, depth::GLsizei, layers::GLsizei, flags::GLbitfield)::Void
export glTexStorageSparseAMD
@getCFun "libGL" glTextureStorageSparseAMD glTextureStorageSparseAMD(texture::GLuint, target::GLenum, internalFormat::GLenum, width::GLsizei, height::GLsizei, depth::GLsizei, layers::GLsizei, flags::GLbitfield)::Void
export glTextureStorageSparseAMD
@getCFun "libGL" glStencilOpValueAMD glStencilOpValueAMD(face::GLenum, value::GLuint)::Void
export glStencilOpValueAMD
@getCFun "libGL" glTessellationFactorAMD glTessellationFactorAMD(factor::GLfloat)::Void
export glTessellationFactorAMD
@getCFun "libGL" glTessellationModeAMD glTessellationModeAMD(mode::GLenum)::Void
export glTessellationModeAMD

end


text:
"""
    transferentropy(source::AbstractArray{<:Real, 1}, response::AbstractArray{<:Real, 1};
        dim = 3, η = 1, τ = 1,
        estimator = VisitationFrequency(),
        n_subdivs = 3)

Compute transfer entropy from `source` to `target` using a `dim`-dimensional delay
reconstruction with prediction lag `η` and embedding lag `τ`,  inferring appropriate
bin sizes from the time series length. Extra dimensions are assigned to the ``T_{pp}``
component of the delay reconstruction (see [`TEVars`](@ref)), while 1D embeddings
are used for the ``T_f`` and ``S_pp`` component.

Returns one transfer entropy estimate per partition.
"""
function transferentropy(source::AbstractArray{<:Real, 1}, response::AbstractArray{<:Real, 1};
    dim = 3, η = 1, τ = 1,
    estimator = VisitationFrequency(),
    n_subdivs = 3)

    dim >= 3 || throw(ArgumentError("dim must be 3 or higher for regular TE"))

    # Distribute dimensions with as (Tf -> k = 1), (Tpp -> l = dim - 2), (Spp -> m = 1),
    # so that it is the states present/past of the target variable that gets a longer
    # history.
    k, l, m = 1, dim - 2, 1
    pts, vars = te_embed(source, response, k, l, m, η = η, τ = τ)

    # Determine appropriate binnings from time series length (according to
    # Krakovska et al. (2018)'s recommendations)
    L = length(source)
    n_subdivisions_coarsest = ceil(Int, L^(1/(dim + 1)))
    n_subdivisions_finest = n_subdivisions_coarsest + n_subdivs
    binning_scheme = map(n-> RectangularBinning(n), n_subdivisions_coarsest:n_subdivisions_finest)

    # Compute TE over different partitions
    # ====================================
    tes = map(b -> transferentropy(pts, vars, b, estimator), binning_scheme)

    return tes
end

"""
    transferentropy(source::AbstractArray{<:Real, 1}, response::AbstractArray{<:Real, 1},
        binning_scheme::Union{RectangularBinning, AbstractArray{RectangularBinning}};
        dim = 3, η = 1, τ = 1, estimator = VisitationFrequency())

Compute transfer entropy from `source` to `target` using a `dim`-dimensional delay
reconstruction with prediction lag `η` and embedding lag `τ`, inferring appropriate bin
sizes from time series length. Extra dimensions are assigned to the ``T_{pp}`` component of
the delay reconstruction (see [`TEVars`](@ref)), while 1D embeddings
are used for the ``T_f`` and ``S_pp`` component.

Returns one transfer entropy estimate per partition.
"""
function transferentropy(source::AbstractArray{<:Real, 1}, response::AbstractArray{<:Real, 1},
    binning_scheme::Union{RectangularBinning, AbstractArray{RectangularBinning}};
    dim = 3, η = 1, τ = 1, estimator = VisitationFrequency())

    dim >= 3 || throw(ArgumentError("dim must be 3 or higher for regular TE"))

    # Distribute dimensions with as (Tf -> k = 1), (Tpp -> l = dim - 2), (Spp -> m = 1),
    # so that it is the states present/past of the target variable that gets a longer
    # history.
    k, l, m = 1, dim - 2, 1
    pts, vars = te_embed(source, response, k, l, m, η = η, τ = τ)

    # Compute TE over different partitions
    # ====================================
    if binning_scheme isa RectangularBinning
        return transferentropy(pts, vars, binning_scheme, estimator)
    else
        return map(b -> transferentropy(pts, vars, b, estimator), binning_scheme)
    end
end

"""
    transferentropy(source::AbstractArray{<:Real, 1}, response::AbstractArray{<:Real, 1},
        k::Int, l::Int, m::Int;
        η = 1, τ = 1, estimator = VisitationFrequency(), n_subdivs = 3)

Compute transfer entropy from `source` to `target` using a `k + l + m`-dimensional delay
reconstruction with `k` being the dimension of the ``T_f`` component, `l` the
dimension of the ``T_{pp}`` component and `m` the dimension of the ``S_{pp}`` component
(see [`TEVars`](@ref)), using prediction lag `η` and embedding lag `τ`,
inferring appropriate bin sizes from time series length.

Returns one transfer entropy estimate per partition.
"""
function transferentropy(source::AbstractArray{<:Real, 1}, response::AbstractArray{<:Real, 1},
        k::Int, l::Int, m::Int;
        η = 1, τ = 1, estimator = VisitationFrequency(), n_subdivs = 3)

    k + l + m >= 3 || throw(ArgumentError("`dim = k + l + m` must be 3 or higher for regular TE"))

    pts, vars = te_embed(source, response, k, l, m, η = η, τ = τ)

    # Determine appropriate binnings from time series length (according to
    # Krakovska et al. (2018)'s recommendations)
    L = length(source)
    n_subdivisions_coarsest = ceil(Int, L^(1/(k + l + m + 1)))
    n_subdivisions_finest = n_subdivisions_coarsest + n_subdivs
    binning_scheme = map(n-> RectangularBinning(n), n_subdivisions_coarsest:n_subdivisions_finest)

    # Compute TE over different partitions
    # ====================================
    tes = map(b -> transferentropy(pts, vars, b, estimator), binning_scheme)

    return tes
end

"""
    transferentropy(source::AbstractArray{<:Real, 1}, response::AbstractArray{<:Real, 1},
        binning_scheme::Union{RectangularBinning, AbstractArray{RectangularBinning}},
        k::Int, l::Int, m::Int;
        dim = 3, η = 1, τ = 1, estimator = VisitationFrequency())

Compute transfer entropy from `source` to `target` using a `k + l + m`-dimensional delay
reconstruction with `k` being the dimension of the ``T_f`` component, `l` the dimension
of the ``T_{pp}`` component and `m` the dimension of the ``S_{pp}`` component
(see [`TEVars`](@ref)), using prediction lag `η` and embedding lag `τ`,
explicitly specifying the partition(s).

Returns one transfer entropy estimate per partition.
"""
function transferentropy(source::AbstractArray{<:Real, 1}, response::AbstractArray{<:Real, 1},
        binning_scheme::Union{RectangularBinning, AbstractArray{RectangularBinning}},
        k::Int, l::Int, m::Int;
        dim = 3, η = 1, τ = 1, estimator = VisitationFrequency())

    k + l + m >= 3 || throw(ArgumentError("`dim = k + l + m` must be 3 or higher for regular TE"))

    # Distribute dimensions with as (Tf -> k = 1), (Tpp -> l = dim - 2), (Spp -> m = 1),
    # so that it is the states present/past of the target variable that gets a longer
    # history.
    pts, vars = te_embed(source, response, k, l, m, η = η, τ = τ)

    # Compute TE over different partitions
    # ====================================
    if binning_scheme isa RectangularBinning
        return transferentropy(pts, vars, binning_scheme, estimator)
    else
        return map(b -> transferentropy(pts, vars, b, estimator), binning_scheme)
    end
end

export transferentropy


text:
# logbeta
using CairoMakie, LaTeXStrings, SpecialFunctions
let
    x = 0.1:0.1:15
    cycle = Cycle([:color, :linestyle, :marker], covary = true)
    set_theme!(Lines = (cycle = cycle,), Scatter = (cycle = cycle,))
    fig = Figure(resolution = (600, 400), font = "CMU Serif", backgroundcolor = :honeydew) # probably you need to install this font in your system
    ax = Axis(fig[1, 1], xlabel = L"x", ylabel = L"Y_{\nu}(x)", ylabelsize = 22,
        xlabelsize = 22, xgridstyle = :dash, ygridstyle = :dash, xtickalign = 1,
        xticksize = 10, ytickalign = 1, yticksize = 10, xlabelpadding = -10)
    for ν in 0:4
        CairoMakie.lines!(ax, x, logbeta.(ν, x), label = latexstring("Y_{$(ν)}(x)"), linewidth = 3)
    end
    axislegend(; position = :rb, nbanks = 2, framecolor = (:grey, 0.5))
    CairoMakie.ylims!(2.5, -10)
    display(fig)
end;

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
__precompile__(false)

module DSB2018

using Reexport

include("errors.jl")
@reexport using .Errors

include("data.jl")
@reexport using .Data

include("model.jl")
@reexport using .Model

end


text:
abstract Broadcast <: Op
type Add <: Broadcast; Add(;o...)=new(); end
type Mul <: Broadcast; Mul(;o...)=new(); end

ninputs(::Broadcast)=2
canoverwrite(::Add)=true        # add can overwrite but only if output same size as one of the inputs
back_reads_x(::Add)=false
back_reads_y(::Add)=false
canoverwrite(::Mul)=false       # mul needs x for back, so cannot overwrite
back_reads_x(::Mul)=true
back_reads_y(::Mul)=false

# TODO: test this with 2d, 4d, lstm, irnn, minibatch=1 etc.
# better yet prove that it will work.
function infersize(::Broadcast, D...)
    findfirst(d->!isa(d,Void), D) == 0 && return nothing
    maxd = 0; for d in D
        d != nothing && length(d) > maxd && (maxd = length(d))
    end
    A = cell(length(D)); for i=1:length(D)
        A[i] = (D[i] == nothing ? zeros(Int,maxd) : Int[D[i]...])
    end
    length(A[end]) == maxd || throw(DimensionMismatch("$A"))
    for i=1:maxd
        maxi = 0; for a in A
            length(a) >= i && a[i] >= maxi && (maxi = a[i])
        end
        for a in A
            length(a) < i ? continue :
            a[i] == maxi ? continue :
            a[i] == 1 ? continue :
            a[i] == 0 ? a[i] = maxi :
            throw(DimensionMismatch("$A"))
        end
        A[end][i] == maxi || throw(DimensionMismatch("$A"))
    end
    ntuple(i->ntuple(j->A[i][j], length(A[i])), length(A))
end

### LEVEL 1: forw/back Handle `nothing`

function forw(::Add,x1,x2,y; o...)
    x1==x2==nothing ? nothing :         # `nothing` represents the zero array to avoid unnecessary fill! in forw and copy! in stack
    x1==nothing ? broadcast!(+,y,x2) :	# use broadcast here instead of copy in case x,y have different sizes
    x2==nothing ? broadcast!(+,y,x1) :
    broadcast!(+,y,x1,x2)               # forw uses julia broadcast semantics
end

function forw(::Mul,x1,x2,y;o...)
    x1==nothing ? nothing :             # `nothing` represents zero for mul also
    x2==nothing ? nothing :             # TODO: the alternative is to have it represent the ones array?
    broadcast!(*,y,x1,x2)               # but I had trouble with that before for a reason I don't remember.
end

function back(::Add,dy,dx1,dx2; o...)
    dx1!=nothing && addback(dy,dx1)
    dx2!=nothing && addback(dy,dx2)
end

function back(::Mul,dy,dx1,dx2; x=nothing, o...)
    (x==nothing || length(x)!=2) && error("back(mul) needs inputs x1 and x2")
    (x1,x2) = x
    dx1!=nothing && (x2==nothing ? fillsync!(dx1,0) : mulback(dy,x2,dx1))
    dx2!=nothing && (x1==nothing ? fillsync!(dx2,0) : mulback(dy,x1,dx2))
end

### LEVEL 2: broadcast!, add/mulback: select best op based on size

import Base: broadcast!, sum!

@gpu function broadcast!{T}(f::Function, y::CudaArray{T}, x::CudaArray{T})
    if f===+
        if size(y)===size(x)
            # Base.warn_once(:ADD2FORW)
            y===x || copysync!(y,x)  # addirnn=11.95 addlstm=2.48 copyseq=11.51 rnnlm=22.84
        else
            Base.warn_once(:BADD2FORW_NOT_TESTED)
            fill!(y,0)
            cudnnAddTensor(x,y)
            # baddforw(x,y)
        end
    else
        error("$f not supported")
    end
    return y
end

@gpu function broadcast!{T}(f::Function, y::CudaArray{T}, x1::CudaArray{T}, x2::CudaArray{T})
    if f===+
        if size(y)===size(x1)===size(x2)
            # Base.warn_once(:ADD3FORW)
            # geam3(x1,x2,y)    # addirnn=12.18 addlstm=2.52 copyseq=11.59 rnnlm=22.97
            axpy3(x1,x2,y)      # addirnn=11.95 addlstm=2.48 copyseq=11.51 rnnlm=22.84
        else
            # Base.warn_once(:BADD3FORW)
            # baddforw(x1,x2,y)      # mnist2d=4.27 4d=15.62 addirnn=14.75 addlstm=2.99 copyseq=15.39 rnnlm=22.28
            cudnnAddTensor3(x1,x2,y) # mnist2d=3.60 4d=14.42 addirnn=11.95 addlstm=2.48 copyseq=11.51 rnnlm=22.84
        end
    elseif f===*
        if size(y)===size(x1)===size(x2)
            # Base.warn_once(:MUL3FORW)
            mul3(x1,x2,y)       # addlstm=2.48 copyseq=11.51 rnnlm=22.84
        else
            Base.warn_once(:BMUL3FORW_NOT_TESTED)
            bmulforw(x1,x2,y)
        end
    else
        error("$f not supported")
    end
    return y
end

# addback sets dx=dy except if different sizes appropriate dims of dy are summed up
function addback(dy,dx)
    if dy === dx
        # done
    elseif size(dy) == size(dx)
        # Base.warn_once(:ADDBACK)
        copysync!(dx,dy)                    # addirnn=11.95 addlstm=2.48 copyseq=11.51 rnnlm=22.84
    else
        # Base.warn_once(:BADDBACK)
        sum!(dx,dy)
    end
end

# warning: this works for sum(a,2) but not for sum(a,1)
@gpu function sum!{T}(dx::CudaArray{T},dy::CudaArray{T})
    # baddback(dy,dx)                   # mnist2d=3.70 4d=22.36 addirnn=14.35 addlstm=2.90 copyseq=11.32 rnnlm=20.94   nondeterministic.
    cudnnConvolutionBackwardBias(dy,dx) # mnist2d=3.60 4d=14.42 addirnn=11.95 addlstm=2.48 copyseq=11.51 rnnlm=22.84
end

# mulback has dx1=dy*x2 and dx2=dy*x1 with appropriate broadcasting sums
function mulback(dy,x1,dx2)
    if size(dy)==size(x1)==size(dx2)
        # Base.warn_once(:MULBACK)
        mul3(dy,x1,dx2)         # addlstm=2.48 copyseq=11.51 rnnlm=22.84
    else
        Base.warn_once(:BMULBACK_NOT_TESTED)
        bmulback(dy,x1,dx2)
    end
end

### LEVEL 3: Actual implementations

# 3a. same size add

function axpy3(a,b,c)
    size(a)==size(b)==size(c) || throw(DimensionMismatch())
    c===a ? axpy!(1,b,c) :
    c===b ? axpy!(1,a,c) :
    (copysync!(c,b); axpy!(1,a,c))
    gpusync(); return c
end

@gpu function geam3{T}(a::CudaArray{T},b::CudaArray{T},c::CudaArray{T})
    size(a)==size(b)==size(c) || throw(DimensionMismatch())
    CUBLAS.geam!('N','N',T(1),mat2d(a),T(1),mat2d(b),mat2d(c))
    gpusync(); return c
end

# 3b. same size mul

@gpu function mul3{T}(a::CudaArray{T},b::CudaArray{T},c::CudaArray{T})
    size(a)==size(b)==size(c) || throw(DimensionMismatch())
    T <: Float32 ? ccall((:mul32,libknet),Void,(Cint,Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}), length(c),a,b,c) :
    T <: Float64 ? ccall((:mul64,libknet),Void,(Cint,Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}),length(c),a,b,c) :
    error("$T not supported")
    gpusync(); return c
end

function mul3{T}(a::Array{T},b::Array{T},c::Array{T})
    size(a)==size(b)==size(c) || throw(DimensionMismatch())
    @inbounds for i=1:length(c); c[i]=a[i]*b[i]; end
    return c
end

# 3b. broadcasting add

@gpu function cudnnAddTensor3{T}(a::CudaArray{T},b::CudaArray{T},c::CudaArray{T})
    c===a ? cudnnAddTensor(b,c) :
    c===b ? cudnnAddTensor(a,c) :
    size(c)==size(b) ? (copysync!(c,b); cudnnAddTensor(a,c)) :
    size(c)==size(a) ? (copysync!(c,a); cudnnAddTensor(b,c)) :
    throw(DimensionMismatch())
    gpusync(); return c
end

# At least avoid alloc for cudadims (still paying copy cost)
_x1dims = CudaArray(Cint,8)
_x2dims = CudaArray(Cint,8)
_ydims = CudaArray(Cint,8)

@gpu function baddforw{T}(x::CudaArray{T}, y::CudaArray{T})
    ndims(x) <= ndims(y) <= 8 || error("Only xdims<=ydims<=8 supported")
    cudadims!(_x1dims, size(x), size(y))
    cudadims!(_ydims, size(y))
    T <: Float32 ? ccall((:badd2forw32,libknet),Void,(Cint,Cint,Ptr{Cint},Ptr{Cfloat}, Ptr{Cint},Ptr{Cfloat}), length(y),ndims(y),_x1dims,x,_ydims,y) :
    T <: Float64 ? ccall((:badd2forw64,libknet),Void,(Cint,Cint,Ptr{Cint},Ptr{Cdouble},Ptr{Cint},Ptr{Cdouble}),length(y),ndims(y),_x1dims,x,_ydims,y) :
    error("$T not supported")
    gpusync(); return y
end

@gpu function baddforw{T}(x1::CudaArray{T}, x2::CudaArray{T}, y::CudaArray{T})
    (ndims(x1) <= ndims(y) && ndims(x2) <= ndims(y) && ndims(y) <= 8) || error("Only xdims<=ydims<=8 supported")
    cudadims!(_x1dims, size(x1), size(y))
    cudadims!(_x2dims, size(x2), size(y))
    cudadims!(_ydims, size(y))
    T <: Float32 ? ccall((:badd3forw32,libknet),Void,(Cint,Cint,Ptr{Cint},Ptr{Cfloat}, Ptr{Cint},Ptr{Cfloat}, Ptr{Cint},Ptr{Cfloat}), length(y),ndims(y),_x1dims,x1,_x2dims,x2,_ydims,y) :
    T <: Float64 ? ccall((:badd3forw64,libknet),Void,(Cint,Cint,Ptr{Cint},Ptr{Cdouble},Ptr{Cint},Ptr{Cdouble},Ptr{Cint},Ptr{Cdouble}),length(y),ndims(y),_x1dims,x1,_x2dims,x2,_ydims,y) :
    error("$T not supported")
    gpusync(); return y
end

# 3c. broadcasting add back

@gpu function baddback{T}(dy::CudaArray{T}, dx::CudaArray{T})
    ndims(dx) <= ndims(dy) <= 8 || error("Only _x1dims<=ydims<=8 supported")
    cudadims!(_x1dims, size(dx), size(dy))
    cudadims!(_ydims, size(dy))
    fill!(dx,0)
    T <: Float32 ? ccall((:badd2back32,libknet),Void,(Cint,Cint,Ptr{Cint},Ptr{Cfloat}, Ptr{Cint},Ptr{Cfloat}), length(dy),ndims(dy),_ydims,dy,_x1dims,dx) :
    T <: Float64 ? ccall((:badd2back64,libknet),Void,(Cint,Cint,Ptr{Cint},Ptr{Cdouble},Ptr{Cint},Ptr{Cdouble}),length(dy),ndims(dy),_ydims,dy,_x1dims,dx) :
    error("$T not supported")
    gpusync();
end

@gpu function bmulback{T}(dy::CudaArray{T}, x1::CudaArray{T}, dx2::CudaArray{T})
    error(:BMULBACK_NOT_IMPLEMENTED)
    ndims(dx) <= ndims(dy) <= 8 || error("Only xdims<=ydims<=8 supported")
    cudadims!(_x1dims, size(dx), size(dy))
    cudadims!(_ydims, size(dy))
    fill!(dx,0)
    T <: Float32 ? ccall((:bmul3back32,libknet),Void,(Cint,Cint,Ptr{Cint},Ptr{Cfloat}, Ptr{Cint},Ptr{Cfloat}), length(dy),ndims(dy),_ydims,dy,_x1dims,dx) :
    T <: Float64 ? ccall((:bmul3back64,libknet),Void,(Cint,Cint,Ptr{Cint},Ptr{Cdouble},Ptr{Cint},Ptr{Cdouble}),length(dy),ndims(dy),_ydims,dy,_x1dims,dx) :
    error("$T not supported")
    gpusync();
end

@gpu function cudadims!(c::CudaArray{Cint}, xdims::Dims, ydims::Dims)
    d = ones(Cint, length(ydims))
    @inbounds for i=1:length(xdims)
        d[i] = xdims[i]
        d[i] == 1 || d[i] == ydims[i] || throw(DimensionMismatch("$xdims,$ydims"))
    end
    copysync!(c,1,d,1,length(d))
end

@gpu cudadims!(c::CudaArray{Cint},ydims::Dims)=copysync!(c, 1, Cint[ydims...], 1, length(ydims))


text:
export
    CachedElement,
    get_element,
    is_up_to_date,
    time_stamp,
    set_up_to_date!,
    set_element!,
    set_time_stamp!,
    update_element!

"""
    CachedElement{E}

A mutable container for caching things.
"""
mutable struct CachedElement{E}
    element::E
    is_up_to_date::Bool
    timestamp::Float64
end
CachedElement(element,flag) = CachedElement(element,flag,time())
CachedElement(element) = CachedElement(element,false)
"""
    get_element(n::CachedElement)

Retrieve the element stored in n.
"""
get_element(n::CachedElement) = n.element
"""
    is_up_to_date(n::CachedElement)

Check if the n is up to date.
"""
is_up_to_date(n::CachedElement) = n.is_up_to_date
"""
    time_stamp(n::CachedElement)

Return the time_stamp at which n was last modified.
"""
time_stamp(n::CachedElement) = n.timestamp
"""
    set_up_to_date!(n::CachedElement,val::Bool=true)

Set the time stamp of n.
"""
function set_up_to_date!(n::CachedElement,val::Bool=true)
    n.is_up_to_date = val
end
"""
    set_element!(n::CachedElement,g)

Set the element stored in n. Does NOT set the `is_up_to_date` flag.
"""
function set_element!(n::CachedElement,g)
    n.element = g
end
"""
    set_time_stamp!(n::CachedElement,t=time())

Set the time stamp of n.
"""
function set_time_stamp!(n::CachedElement,t=time())
    n.timestamp = t
end
"""
    update_element!(n::CachedElement,g)

Set the element of n to g, update the `is_up_to_date` flag and the time stamp.
"""
function update_element!(n::CachedElement,g)
    set_element!(n,g)
    set_up_to_date!(n,true)
    set_time_stamp!(n)
    return g
end

"""
    Base.copy(e::CachedElement)

Shares the e.element, since it doesn't need to be replaced until `set_element!`
is called. Copies `e.is_up_to_date` to preserve the cache state.
"""
Base.copy(e::CachedElement) = CachedElement(e.element,copy(is_up_to_date(e)),time())
Base.convert(::Type{CachedElement{T}},e::CachedElement) where {T} = CachedElement{T}(get_element(e),is_up_to_date(e),time())

const cached_element_accessor_interface = [:is_up_to_date,:time_stamp]
const cached_element_mutator_interface = [:update_element!,:set_time_stamp!,:set_element!,:set_up_to_date!]

export
    Toggle,
    set_toggle_status!,
    get_toggle_status

mutable struct Toggle
    status::Bool
end
function set_toggle_status!(t::Toggle,val=true)
    t.status = val
end
get_toggle_status(t::Toggle) = copy(t.status)

export
    Counter,
    set_counter_status!,
    get_counter_status

mutable struct Counter
    status::Int
end
function set_counter_status!(t::Counter,val::Int)
    t.status = val
end
get_counter_status(t::Counter) = copy(t.status)