text:
using DocStringExtensions

"""
Compute ABR peak points on the basis of peak latencies in another condition.
The delay is recomputed on the basis of peak V delay if peak V is found.
"""
function calcDelayedPeakPoints(sig, pkTimes, delay, sampRate; epochStart=0)
    nPeaks = length(pkTimes)
    xpksPnts, xpksTimes = findPeaks(sig, sampRate, epochStart=epochStart)
    xtrsPnts, xtrsTimes = findTroughs(sig, sampRate, epochStart=epochStart)

    ## try to find peak V first
    minPeakVLat = pkTimes[5]+delay[1]
    maxPeakVLat = pkTimes[5]+delay[2]

    if ismissing(minPeakVLat) == false
        peakVTimes = xpksTimes[(xpksTimes .>= minPeakVLat) .& (xpksTimes .<= maxPeakVLat)]
        peakVPnts = xpksPnts[(xpksTimes .>= minPeakVLat) .& (xpksTimes .<= maxPeakVLat)]
        if length(peakVTimes) > 0
            ## find largest peak
            idx = findall(sig[peakVPnts] .== maximum(sig[peakVPnts]))[1]
            peakVTime = peakVTimes[idx]
            peakVDelay = peakVTime - pkTimes[5]
            ## recompute delay on the basis of peak V delay
            delay = (peakVDelay-0.00035, peakVDelay+0.00035)
        end
    end
    peakTimes = missings(Float64, nPeaks)
    peakPoints = missings(Int, nPeaks)
    for pk=1:nPeaks
        if ismissing(pkTimes[pk]) == false
            fooTimes = xpksTimes[(xpksTimes .>= pkTimes[pk]+delay[1]) .& (xpksTimes .<= pkTimes[pk]+delay[2])]
            fooPnts = xpksPnts[(xpksTimes.>=pkTimes[pk]+delay[1]) .& (xpksTimes.<=pkTimes[pk]+delay[2])]
        else
            fooTimes = (Float64)[]
            fooPnts = (Int)[]
        end

        if length(fooTimes) > 0
            #find largest peak
            idx = findall(sig[fooPnts] .== maximum(sig[fooPnts]))[1]
            peakTimes[pk] = fooTimes[idx]
            peakPoints[pk] = fooPnts[idx]
        else
            peakTimes[pk] = missing
            peakPoints[pk] = missing
        end
    end

    minPeakTroughLat = [0.25/1000, 0.12/1000, 0.25/1000, 0.12/1000, 0.25/1000]
    maxPeakTroughLatDiff = [0.75/1000, missing, 0.75/1000, missing, 1.5/1000]
    maxPeakTroughLat = minPeakTroughLat + maxPeakTroughLatDiff

    troughTimes = missings(Float64, nPeaks)
    troughPoints = missings(Int, nPeaks)
    minTroughLat = missings(Float64, nPeaks)
    maxTroughLat = missings(Float64, nPeaks)
    for pk=1:nPeaks
        thisMin = minPeakTroughLat[pk]
        if ismissing(maxPeakTroughLat[pk]) == true
            if pk+1 <= 5
                thisMax = peakTimes[pk+1] - peakTimes[pk]
            end
        else
            thisMax = maxPeakTroughLat[pk]
        end
        minTroughLat[pk] = peakTimes[pk]+thisMin
        maxTroughLat[pk] = peakTimes[pk]+thisMax

        if (ismissing(peakTimes[pk]+thisMin+thisMax) == false)
            fooTimes = xtrsTimes[(xtrsTimes .>= peakTimes[pk]+thisMin) .& (xtrsTimes .<= peakTimes[pk]+thisMax)]
            fooPnts = xtrsPnts[(xtrsTimes.>= peakTimes[pk]+thisMin) .& (xtrsTimes.<= peakTimes[pk]+thisMax)]
        else
            fooTimes = (Float64)[]
            fooPnts = (Int)[]
        end

        if length(fooTimes) > 0
            #find largest trough
            idx = findall(sig[fooPnts] .== minimum(sig[fooPnts]))[1]
            troughTimes[pk] = fooTimes[idx]
            troughPoints[pk] = fooPnts[idx]
        else
            troughTimes[pk] = missing
            troughPoints[pk] = missing
        end
    end

    peakLatencies = peakTimes
    troughLatencies = troughTimes

    peakAmps = missings(Float64, length(peakPoints))
    troughAmps = missings(Float64, length(peakPoints))
    for i=1:length(peakAmps)
        if ismissing(peakPoints[i]) == false
            peakAmps[i] = sig[1, peakPoints[i]]
        else
            peakAmps[i] = missing
        end

        if ismissing(troughPoints[i]) == false
            troughAmps[i] = sig[1, troughPoints[i]]
        else
            troughAmps[i] = missing
        end
    end

    peakTroughAmps = peakAmps-troughAmps

    waveLabels = ["I", "II", "III", "IV", "V"]
    df = DataFrame(wave=waveLabels, peakPoint=peakPoints, troughPoint=troughPoints, peakLatency=peakLatencies, troughLatency=troughLatencies, peakAmp=peakAmps, troughAmp=troughAmps, peakTroughAmp=peakTroughAmps,
                   minPeakLat=vcat(pkTimes[1:4].+delay[1], minPeakVLat),
                   maxPeakLat=vcat(pkTimes[1:4].+delay[2], maxPeakVLat),
                   minTroughLat=minTroughLat,
                   maxTroughLat=maxTroughLat)

    return df

end

"""
Given a vector of expected peak times (e.g. from a grand average),
find peak points within a time window of exp_pk_time +/- delay.

The function first finds peaks

$(SIGNATURES)

##### Arguments

* `sig::Union{AbstractMatrix{T}, AbstractVector{T}}`: the waveform for which the peaks and troughs are sought.
* `expPeakLatencies::AbstractVector{Real}`: vector of expected latencies of the peaks sought, in seconds.
* `winStart::AbstractVector{Real}`: the start of the time windows relative to the expected peak latencies in which to find peaks, in seconds.
* `winStop::AbstractVector{Real}`: the end of the time windows relative to the expected peak latencies in which to find peaks, in seconds.
* `sampRate::Real`: sampling rate of the waveform in Hz.
* `epochStart::Real`: the time, in seconds, at which the epoch starts relative to the starting time of the stimulus

##### Returns

* ``:

##### Examples

```julia

```

"""
function calcPeaksFromExpectedLatencies(sig::Union{AbstractMatrix{T}, AbstractVector{T}},
                                        expPeakLatencies::Union{AbstractVector{S}, AbstractVector{Union{S, Missing}}},
                                        winStart::AbstractVector{P}, winStop::AbstractVector{Q}, sampRate::Real;
                                        epochStart::Real=0,
                                        minPeakTroughLat::AbstractVector{G}=[0.25, 0.12, 0.25, 0.12, 0.25]./1000,
                                        maxPeakTroughLat::AbstractVector{H}=[1, 1, 1, 0.75, 2]./1000,
                                        minAbsLat::Union{AbstractVector{S}, AbstractVector{Union{R, Missing}}}=[missing for i=1:5],
                                        maxAbsLat::Union{AbstractVector{S}, AbstractVector{Union{W, Missing}}}=[missing for i=1:5],
                                        waveLabels::AbstractVector{String} = ["I", "II", "III", "IV", "V"]) where {T<:Real, S<:Real, P<:Real, Q<:Real, G<:Real, H<:Real, R<:Real, W<:Real}
    nPeaks = length(expPeakLatencies)
    xpksPnts, xpksTimes = findPeaks(sig, sampRate, epochStart=epochStart) #find all peaks in the waveform
    xtrsPnts, xtrsTimes = findTroughs(sig, sampRate, epochStart=epochStart) #find all troughs in the waveform

    #############################
    ## find peaks
    peakLatencies = missings(Float64, nPeaks)
    peakPoints = missings(Int, nPeaks)
    peakLatenciesEst = missings(Float64, nPeaks)
    peakPointsEst = missings(Int, nPeaks)

    searchWinStart = missings(Float64, nPeaks)
    searchWinStop = missings(Float64, nPeaks)
    for pk=1:nPeaks
        searchWinStart[pk] = max(expPeakLatencies[pk]+winStart[pk], ifelse(ismissing(minAbsLat[pk]), expPeakLatencies[pk]+winStart[pk], minAbsLat[pk]))
        searchWinStop[pk] = min(expPeakLatencies[pk]+winStop[pk], ifelse(ismissing(maxAbsLat[pk]), expPeakLatencies[pk]+winStop[pk], maxAbsLat[pk]))
    end
    for pk=1:nPeaks
        if ismissing(expPeakLatencies[pk]) == false
            candidateTimes = xpksTimes[(xpksTimes .>= searchWinStart[pk]) .& (xpksTimes .<= searchWinStop[pk])]
            candidatePnts = xpksPnts[(xpksTimes.>=searchWinStart[pk]) .& (xpksTimes.<=searchWinStop[pk])]
        else
            candidateTimes = (Float64)[]
            candidatePnts = (Int)[]
        end
        if length(candidateTimes) > 0
            #find largest peak
            idx = findall(sig[candidatePnts] .== maximum(sig[candidatePnts]))[1]
            peakLatencies[pk] = candidateTimes[idx]
            peakPoints[pk] = candidatePnts[idx]
            peakLatenciesEst[pk] = candidateTimes[idx]
            peakPointsEst[pk] = candidatePnts[idx]
        else
            peakLatencies[pk] = missing
            peakPoints[pk] = missing
            if ismissing(expPeakLatencies[pk]) == false
                idxStart = round(Int, (searchWinStart[pk]-epochStart)*sampRate)+1;
                idxStop = round(Int, (searchWinStop[pk]-epochStart)*sampRate)+1;
                idx = findfirst(sig[idxStart:idxStop] .== maximum(sig[idxStart:idxStop])) + idxStart-1
                peakPointsEst[pk] = idx
                peakLatenciesEst[pk] = ((idx-1)/sampRate)+epochStart
            else
                peakLatenciesEst[pk] = missing
                peakPointsEst[pk] = missing
            end
        end
    end

    ###############################
    ## determine trough windows from peak latencies
    troughLatencies = missings(Float64, nPeaks)
    troughPoints = missings(Int, nPeaks) #has to be float to possibly be missing
    troughLatenciesEst = missings(Float64, nPeaks)
    troughPointsEst = missings(Int, nPeaks) #has to be float to possibly be missing
    minTroughLat = missings(Float64, nPeaks)
    maxTroughLat = missings(Float64, nPeaks)
    minTroughLatEst = missings(Float64, nPeaks)
    maxTroughLatEst = missings(Float64, nPeaks)
    for pk=1:nPeaks
        thisMin = minPeakTroughLat[pk]
        thisMax = maxPeakTroughLat[pk]
        ## if ismissing(maxPeakTroughLat[pk]) == true
        ##     if pk+1 <= nPeaks
        ##         thisMax = peakLatencies[pk+1] - peakLatencies[pk]
        ##     end
        ## else
        ##     thisMax = maxPeakTroughLat[pk]
        ## end
        minTroughLat[pk] = peakLatencies[pk]+thisMin
        maxTroughLat[pk] = peakLatencies[pk]+thisMax
        minTroughLatEst[pk]= peakLatenciesEst[pk]+thisMin
        maxTroughLatEst[pk] =  peakLatenciesEst[pk]+thisMax

        if ismissing(peakLatencies[pk]) == false
            candidateTimes = xtrsTimes[(xtrsTimes .>= peakLatencies[pk]+thisMin) .& (xtrsTimes .<= peakLatencies[pk]+thisMax)]
            candidatePnts = xtrsPnts[(xtrsTimes.>= peakLatencies[pk]+thisMin) .& (xtrsTimes.<= peakLatencies[pk]+thisMax)]
        else
            candidateTimes = (Float64)[]
            candidatePnts = (Int)[]
        end

        if length(candidateTimes) > 0
            #find largest trough
            idx = findall(sig[candidatePnts] .== minimum(sig[candidatePnts]))[1]
            troughLatencies[pk] = candidateTimes[idx]
            troughPoints[pk] = candidatePnts[idx]
            troughLatenciesEst[pk] = candidateTimes[idx]
            troughPointsEst[pk] = candidatePnts[idx]
        else
            troughLatencies[pk] = missing
            troughPoints[pk] = missing
            if ismissing(expPeakLatencies[pk]) == false
                idxStart = round(Int, (minTroughLatEst[pk]-epochStart)*sampRate)+1;
                idxStop = round(Int, (maxTroughLatEst[pk]-epochStart)*sampRate)+1;
                idx = findfirst(sig[idxStart:idxStop] .== minimum(sig[idxStart:idxStop])) + idxStart-1
                troughPointsEst[pk] = idx
                troughLatenciesEst[pk] = ((idx-1)/sampRate)+epochStart
            else
                troughLatenciesEst[pk] = missing
                troughPointsEst[pk] = missing
            end
        end
    end

    #####################################
    ## measure peak and trough amplitudes
    peakAmps = missings(Float64, nPeaks)#zeros(nPeaks)
    troughAmps = missings(Float64, nPeaks)#zeros(nPeaks)
    peakAmpsEst = missings(Float64, nPeaks)#zeros(nPeaks)
    troughAmpsEst = missings(Float64, nPeaks)#zeros(nPeaks)
    for i=1:nPeaks
        if ismissing(peakPoints[i]) == false
            peakAmps[i] = sig[round(Int, peakPoints[i])]
            peakAmpsEst[i] = sig[round(Int, peakPoints[i])]
        else
            peakAmps[i] = missing
            if ismissing(expPeakLatencies[i]) == false
                #idxStart = round(Int, (expPeakLatencies[i]+winStart[i]-epochStart)*sampRate);
                #idxStop = round(Int, (expPeakLatencies[i]+winStop[i]-epochStart)*sampRate);
                peakAmpsEst[i] = sig[round(Int, peakPointsEst[i])] #maximum(sig[idxStart:idxStop])[1]
            else
                peakAmpsEst[i] = missing
            end
        end

        if ismissing(troughPoints[i]) == false
            troughAmps[i] = sig[round(Int, troughPoints[i])]
            troughAmpsEst[i] = sig[round(Int, troughPoints[i])]
        else
            troughAmps[i] = missing
            if ismissing(minTroughLatEst[i]) == false
                #idxStart = round(Int, (minTroughLat[i]-epochStart)*sampRate);
                #idxStop = round(Int, (maxTroughLat[i]-epochStart)*sampRate);
                troughAmpsEst[i] = sig[round(Int, troughPointsEst[i])] #minimum(sig[idxStart:idxStop])[1]
            else
                troughAmpsEst[i] = missing
            end
        end
    end

    ####################################
    ## measure peak-trough amplitudes
    peakTroughAmps = peakAmps-troughAmps
    peakTroughAmpsEst = peakAmpsEst-troughAmpsEst

    ## #############################################################
    ## Estimate peaks and troughs by finding point that maximizes
    ## peak-trough amplitude
    ## #############################################################
    peakPointsEstPTMax = missings(Int, nPeaks)
    troughPointsEstPTMax = missings(Int, nPeaks)
    peakAmpsEstPTMax = missings(Float64, nPeaks)
    troughAmpsEstPTMax = missings(Float64, nPeaks)
    peakLatenciesEstPTMax = missings(Float64, nPeaks)
    troughLatenciesEstPTMax = missings(Float64, nPeaks)

    for pk=1:nPeaks
        if ismissing(expPeakLatencies[pk]) == false
            minPeakIdx = round(Int, (searchWinStart[pk]-epochStart)*sampRate)
            maxPeakIdx = round(Int, (searchWinStop[pk]-epochStart)*sampRate)

            bestDiff = 0
            pkIdx = 1
            trIdx = 1
            for i=minPeakIdx:maxPeakIdx-1
                for j=i+round(Int, minPeakTroughLat[pk]*sampRate):i+round(Int, maxPeakTroughLat[pk]*sampRate)
                    currDiff = sig[i]-sig[j]
                    if currDiff > bestDiff
                        bestDiff = copy(currDiff)
                        pkIdx = i
                        trIdx = j
                    end
                end
            end

            peakPointsEstPTMax[pk] = pkIdx
            troughPointsEstPTMax[pk] = trIdx
            peakAmpsEstPTMax[pk] = sig[pkIdx]
            troughAmpsEstPTMax[pk] = sig[trIdx]
            peakLatenciesEstPTMax[pk] = (pkIdx-1)/sampRate+epochStart
            troughLatenciesEstPTMax[pk] = (trIdx-1)/sampRate+epochStart
        end
    end
    peakTroughAmpsEstPTMax = peakAmpsEstPTMax-troughAmpsEstPTMax

    ####################################
    ## store results in a dataframe
    df = DataFrame(wave=waveLabels, peakPoint=peakPoints, troughPoint=troughPoints, peakLatency=peakLatencies, troughLatency=troughLatencies, peakAmp=peakAmps, troughAmp=troughAmps, peakTroughAmp=peakTroughAmps, peakPointEst=peakPointsEst, peakLatencyEst=peakLatenciesEst, peakAmpEst=peakAmpsEst, troughPointEst=troughPointsEst, troughLatencyEst=troughLatenciesEst, troughAmpEst=troughAmpsEst, peakTroughAmpEst=peakTroughAmpsEst,
                   peakPointEstPTMax=peakPointsEstPTMax, peakLatencyEstPTMax=peakLatenciesEstPTMax, peakAmpEstPTMax=peakAmpsEstPTMax, troughPointEstPTMax=troughPointsEstPTMax, troughLatencyEstPTMax=troughLatenciesEstPTMax, troughAmpEstPTMax=troughAmpsEstPTMax, peakTroughAmpEstPTMax=peakTroughAmpsEstPTMax,
                   minPeakLat=expPeakLatencies+winStart,
                   maxPeakLat=expPeakLatencies+winStop,
                   minTroughLat=minTroughLat,
                   maxTroughLat=maxTroughLat)

    return df

end

"""
Given a vector of expected peak times (e.g. from a grand average),
find peak points within a time window of exp_pk_time +/- delay.

The difference from `calcPeaksFromExpectedLatencies` is that if troughs are not found inflection points are used to mark troughs.
The function first finds peaks.

$(SIGNATURES)

##### Arguments

* `sig::Union{AbstractMatrix{T}, AbstractVector{T}}`: the waveform for which the peaks and troughs are sought.
* `expPeakLatencies::AbstractVector{Real}`: vector of expected latencies of the peaks sought, in seconds.
* `winStart::AbstractVector{Real}`: the start of the time windows relative to the expected peak latencies in which to find peaks, in seconds.
* `winStop::AbstractVector{Real}`: the end of the time windows relative to the expected peak latencies in which to find peaks, in seconds.
* `sampRate::Real`: sampling rate of the waveform in Hz.
* `epochStart::Real`: the time, in seconds, at which the epoch starts relative to the starting time of the stimulus

##### Returns

* ``:

##### Examples

```julia

```

"""
function calcPeaksTroughsInflFromExpectedLatencies(sig::Union{AbstractMatrix{T}, AbstractVector{T}}, expPeakLatencies::Union{AbstractVector{S}, AbstractVector{Union{S, Missing}}}, winStart::AbstractVector{P}, winStop::AbstractVector{Q}, sampRate::Real; epochStart::Real=0, minPeakTroughLat::AbstractVector{G}=[0.25, 0.12, 0.25, 0.12, 0.25]./1000, maxPeakTroughLat::AbstractVector{H}=[1, 1, 1, 0.75, 2]./1000, waveLabels::AbstractVector{String} = ["I", "II", "III", "IV", "V"]) where {T<:Real, S<:Real, P<:Real, Q<:Real, G<:Real, H<:Real}
    nPeaks = length(expPeakLatencies)
    xpksPnts, xpksTimes = findPeaks(sig, sampRate, epochStart=epochStart) #find all peaks in the waveform
    xtrsPnts, xtrsTimes = findTroughs(sig, sampRate, epochStart=epochStart) #find all troughs in the waveform
    xinflPnts, xinflTimes = findInflections(sig, sampRate, epochStart=epochStart) #find all inflections in the waveform
    dy = diff(vec(sig))
    #############################
    ## find peaks
    peakLatencies = missings(Float64, nPeaks)
    peakPoints = missings(Int, nPeaks)
    peakLatenciesEst = missings(Float64, nPeaks)
    peakPointsEst = missings(Int, nPeaks)
    for pk=1:nPeaks
        if ismissing(expPeakLatencies[pk]) == false
            candidateTimes = xpksTimes[(xpksTimes .>= expPeakLatencies[pk]+winStart[pk]) .& (xpksTimes .<= expPeakLatencies[pk]+winStop[pk])]
            candidatePnts = xpksPnts[(xpksTimes.>=expPeakLatencies[pk]+winStart[pk]) .& (xpksTimes.<=expPeakLatencies[pk]+winStop[pk])]
        else
            candidateTimes = (Float64)[]
            candidatePnts = (Int)[]
        end
        if length(candidateTimes) > 0
            #find largest peak
            idx = findall(sig[candidatePnts] .== maximum(sig[candidatePnts]))[1]
            peakLatencies[pk] = candidateTimes[idx]
            peakPoints[pk] = candidatePnts[idx]
            peakLatenciesEst[pk] = candidateTimes[idx]
            peakPointsEst[pk] = candidatePnts[idx]
        else
            peakLatencies[pk] = missing
            peakPoints[pk] = missing
            if ismissing(expPeakLatencies[pk]) == false
                idxStart = round(Int, (expPeakLatencies[pk]+winStart[pk]-epochStart)*sampRate)+1;
                idxStop = round(Int, (expPeakLatencies[pk]+winStop[pk]-epochStart)*sampRate)+1;
                idx = findfirst(sig[idxStart:idxStop] .== maximum(sig[idxStart:idxStop])) + idxStart-1
                peakPointsEst[pk] = idx
                peakLatenciesEst[pk] = ((idx-1)/sampRate)+epochStart
            else
                peakLatenciesEst[pk] = missing
                peakPointsEst[pk] = missing
            end
        end
    end

    ###############################
    ## determine trough windows from peak latencies
    troughLatencies = missings(Float64, nPeaks)
    troughPoints = missings(Int, nPeaks) #has to be float to possibly be missing
    troughLatenciesEst = missings(Float64, nPeaks)
    troughPointsEst = missings(Int, nPeaks) #has to be float to possibly be missing
    minTroughLat = missings(Float64, nPeaks)
    maxTroughLat = missings(Float64, nPeaks)
    minTroughLatEst = missings(Float64, nPeaks)
    maxTroughLatEst = missings(Float64, nPeaks)
    for pk=1:nPeaks
        thisMin = minPeakTroughLat[pk]
        thisMax = maxPeakTroughLat[pk]
        ## if ismissing(maxPeakTroughLat[pk]) == true
        ##     if pk+1 <= nPeaks
        ##         thisMax = peakLatencies[pk+1] - peakLatencies[pk]
        ##     end
        ## else
        ##     thisMax = maxPeakTroughLat[pk]
        ## end
        minTroughLat[pk] = peakLatencies[pk]+thisMin
        maxTroughLat[pk] = peakLatencies[pk]+thisMax
        minTroughLatEst[pk]= peakLatenciesEst[pk]+thisMin
        maxTroughLatEst[pk] =  peakLatenciesEst[pk]+thisMax

        candidateTimes = (Float64)[]
        candidatePnts = (Int)[]
        if ismissing(peakLatencies[pk]) == false
            candidateTimes = xtrsTimes[(xtrsTimes .>= peakLatencies[pk]+thisMin) .& (xtrsTimes .<= peakLatencies[pk]+thisMax)]
            candidatePnts = xtrsPnts[(xtrsTimes.>= peakLatencies[pk]+thisMin) .& (xtrsTimes.<= peakLatencies[pk]+thisMax)]
            if length(candidateTimes) > 0
                #find largest trough
                idx = findall(sig[candidatePnts] .== minimum(sig[candidatePnts]))[1]
                troughLatencies[pk] = candidateTimes[idx]
                troughPoints[pk] = candidatePnts[idx]
                troughLatenciesEst[pk] = candidateTimes[idx]
                troughPointsEst[pk] = candidatePnts[idx]
            end
        end

        if (ismissing(troughPoints[pk]) == true) & (ismissing(peakLatencies[pk]) == false) #try inflection points
            candidateTimes = xinflTimes[(xinflTimes .>= peakLatencies[pk]+thisMin) .& (xinflTimes .<= peakLatencies[pk]+thisMax)]
            candidatePnts = xinflPnts[(xinflTimes.>=peakLatencies[pk]+thisMin) .& (xinflTimes.<=peakLatencies[pk]+thisMax)]
            if length(candidateTimes) > 0
                idx = findall(abs.(dy[candidatePnts .- 1]) .== minimum(abs.(dy[candidatePnts .- 1])))[1] #choose point where first derivative is closest to zero
                troughLatencies[pk] = candidateTimes[idx]
                troughPoints[pk] = candidatePnts[idx]
                troughLatenciesEst[pk] = candidateTimes[idx]
                troughPointsEst[pk] = candidatePnts[idx]
            end
        end

        if ismissing(troughPointsEst[pk]) == true
            if ismissing(expPeakLatencies[pk]) == false
                idxStart = round(Int, (minTroughLatEst[pk]-epochStart)*sampRate)+1;
                idxStop = round(Int, (maxTroughLatEst[pk]-epochStart)*sampRate)+1;
                idx = findfirst(sig[idxStart:idxStop] .== minimum(sig[idxStart:idxStop])) + idxStart-1
                troughPointsEst[pk] = idx
                troughLatenciesEst[pk] = ((idx-1)/sampRate)+epochStart
            else
                troughLatenciesEst[pk] = missing
                troughPointsEst[pk] = missing
            end
        end
    end

    #####################################
    ## measure peak and trough amplitudes
    peakAmps = missings(Float64, nPeaks)#zeros(nPeaks)
    troughAmps = missings(Float64, nPeaks)#zeros(nPeaks)
    peakAmpsEst = missings(Float64, nPeaks)#zeros(nPeaks)
    troughAmpsEst = missings(Float64, nPeaks)#zeros(nPeaks)
    for i=1:nPeaks
        if ismissing(peakPoints[i]) == false
            peakAmps[i] = sig[round(Int, peakPoints[i])]
            peakAmpsEst[i] = sig[round(Int, peakPoints[i])]
        else
            peakAmps[i] = missing
            if ismissing(expPeakLatencies[i]) == false
                #idxStart = round(Int, (expPeakLatencies[i]+winStart[i]-epochStart)*sampRate);
                #idxStop = round(Int, (expPeakLatencies[i]+winStop[i]-epochStart)*sampRate);
                peakAmpsEst[i] = sig[round(Int, peakPointsEst[i])] #maximum(sig[idxStart:idxStop])[1]
            else
                peakAmpsEst[i] = missing
            end
        end

        if ismissing(troughPoints[i]) == false
            troughAmps[i] = sig[round(Int, troughPoints[i])]
            troughAmpsEst[i] = sig[round(Int, troughPoints[i])]
        else
            troughAmps[i] = missing
            if ismissing(minTroughLatEst[i]) == false
                #idxStart = round(Int, (minTroughLat[i]-epochStart)*sampRate);
                #idxStop = round(Int, (maxTroughLat[i]-epochStart)*sampRate);
                troughAmpsEst[i] = sig[round(Int, troughPointsEst[i])] #minimum(sig[idxStart:idxStop])[1]
            else
                troughAmpsEst[i] = missing
            end
        end
    end

    ####################################
    ## measure peak-trough amplitudes
    peakTroughAmps = peakAmps-troughAmps
    peakTroughAmpsEst = peakAmpsEst-troughAmpsEst

    ## #############################################################
    ## Estimate peaks and troughs by finding point that maximizes
    ## peak-trough amplitude
    ## #############################################################
    peakPointsEstPTMax = missings(Int, nPeaks)
    troughPointsEstPTMax = missings(Int, nPeaks)
    peakAmpsEstPTMax = missings(Float64, nPeaks)
    troughAmpsEstPTMax = missings(Float64, nPeaks)
    peakLatenciesEstPTMax = missings(Float64, nPeaks)
    troughLatenciesEstPTMax = missings(Float64, nPeaks)

    for pk=1:nPeaks
        if ismissing(expPeakLatencies[pk]) == false
            minPeakIdx = round(Int, (expPeakLatencies[pk]+winStart[pk]-epochStart)*sampRate)
            maxPeakIdx = round(Int, (expPeakLatencies[pk]+winStop[pk]-epochStart)*sampRate)

            bestDiff = 0
            pkIdx = 1
            trIdx = 1
            for i=minPeakIdx:maxPeakIdx-1
                for j=i+round(Int, minPeakTroughLat[pk]*sampRate):i+round(Int, maxPeakTroughLat[pk]*sampRate)
                    currDiff = sig[i]-sig[j]
                    if currDiff > bestDiff
                        bestDiff = copy(currDiff)
                        pkIdx = i
                        trIdx = j
                    end
                end
            end

            peakPointsEstPTMax[pk] = pkIdx
            troughPointsEstPTMax[pk] = trIdx
            peakAmpsEstPTMax[pk] = sig[pkIdx]
            troughAmpsEstPTMax[pk] = sig[trIdx]
            peakLatenciesEstPTMax[pk] = (pkIdx-1)/sampRate+epochStart
            troughLatenciesEstPTMax[pk] = (trIdx-1)/sampRate+epochStart
        end
    end
    peakTroughAmpsEstPTMax = peakAmpsEstPTMax-troughAmpsEstPTMax

    ####################################
    ## store results in a dataframe
    df = DataFrame(wave=waveLabels, peakPoint=peakPoints, troughPoint=troughPoints, peakLatency=peakLatencies, troughLatency=troughLatencies, peakAmp=peakAmps, troughAmp=troughAmps, peakTroughAmp=peakTroughAmps, peakPointEst=peakPointsEst, peakLatencyEst=peakLatenciesEst, peakAmpEst=peakAmpsEst, troughPointEst=troughPointsEst, troughLatencyEst=troughLatenciesEst, troughAmpEst=troughAmpsEst, peakTroughAmpEst=peakTroughAmpsEst,
                   peakPointEstPTMax=peakPointsEstPTMax, peakLatencyEstPTMax=peakLatenciesEstPTMax, peakAmpEstPTMax=peakAmpsEstPTMax, troughPointEstPTMax=troughPointsEstPTMax, troughLatencyEstPTMax=troughLatenciesEstPTMax, troughAmpEstPTMax=troughAmpsEstPTMax, peakTroughAmpEstPTMax=peakTroughAmpsEstPTMax,
                   minPeakLat=expPeakLatencies+winStart,
                   maxPeakLat=expPeakLatencies+winStop,
                   minTroughLat=minTroughLat,
                   maxTroughLat=maxTroughLat)

    return df

end

function calcPeaksInflFromExpectedLatencies(sig::Union{AbstractMatrix{T}, AbstractVector{T}}, expPeakLatencies::Union{AbstractVector{S}, AbstractVector{Union{S, Missing}}}, winStart::AbstractVector{P}, winStop::AbstractVector{Q}, sampRate::Real; epochStart::Real=0, minPeakTroughLat::AbstractVector{G}=[0.25, 0.12, 0.25, 0.12, 0.25]./1000, maxPeakTroughLat::AbstractVector{H}=[1, 1, 1, 0.75, 2]./1000, waveLabels::AbstractVector{String} = ["I", "II", "III", "IV", "V"]) where {T<:Real, S<:Real, P<:Real, Q<:Real, G<:Real, H<:Real}
    nPeaks = length(expPeakLatencies)
    xpksPnts, xpksTimes = findPeaks(sig, sampRate, epochStart=epochStart) #find all peaks in the waveform
    xtrsPnts, xtrsTimes = findTroughs(sig, sampRate, epochStart=epochStart) #find all troughs in the waveform
    xinflPnts, xinflTimes = findInflections(sig, sampRate, epochStart=epochStart) #find all inflections in the waveform

    #############################
    ## find peaks
    peakLatencies = missings(Float64, nPeaks)
    peakPoints = missings(Int, nPeaks) #has to be float to possibly be missing
    peakLatenciesEst = missings(Float64, nPeaks) #in case a peak can't be found fall back to finding the maximum value
    peakPointsEst = missings(Int, nPeaks)
    for pk=1:nPeaks
        if ismissing(expPeakLatencies[pk]) == false
            candidateTimes = xpksTimes[(xpksTimes .>= expPeakLatencies[pk]+winStart[pk]) .& (xpksTimes .<= expPeakLatencies[pk]+winStop[pk])]
            candidatePnts = xpksPnts[(xpksTimes.>=expPeakLatencies[pk]+winStart[pk]) .& (xpksTimes.<=expPeakLatencies[pk]+winStop[pk])]
        else
            candidateTimes = (Float64)[]
            candidatePnts = (Int)[]
        end
        if length(candidateTimes) > 0
            #find largest peak
            idx = findall(sig[candidatePnts] .== maximum(sig[candidatePnts]))[1]
            peakLatencies[pk] = candidateTimes[idx]
            peakPoints[pk] = candidatePnts[idx]
            peakLatenciesEst[pk] = candidateTimes[idx]
            peakPointsEst[pk] = candidatePnts[idx]
        else
            ##try inflection points
            if ismissing(expPeakLatencies[pk]) == false
                candidateTimes = xinflTimes[(xinflTimes .>= expPeakLatencies[pk]+winStart[pk]) .& (xinflTimes .<= expPeakLatencies[pk]+winStop[pk])]
                candidatePnts = xinflPnts[(xinflTimes.>=expPeakLatencies[pk]+winStart[pk]) .& (xinflTimes.<=expPeakLatencies[pk]+winStop[pk])]
            else
                candidateTimes = (Float64)[]
                candidatePnts = (Int)[]
            end
            if length(candidateTimes) > 0
                #find largest peak
                idx = findall(sig[candidatePnts] .== maximum(sig[candidatePnts]))[1]
                peakLatencies[pk] = candidateTimes[idx]
                peakPoints[pk] = candidatePnts[idx]
                peakLatenciesEst[pk] = candidateTimes[idx]
                peakPointsEst[pk] = candidatePnts[idx]
            else
                peakLatencies[pk] = missing
                peakPoints[pk] = missing
                if ismissing(expPeakLatencies[pk]) == false
                    idxStart = round(Int, (expPeakLatencies[pk]+winStart[pk]-epochStart)*sampRate)+1;
                    idxStop = round(Int, (expPeakLatencies[pk]+winStop[pk]-epochStart)*sampRate)+1;
                    idx = findfirst(sig[idxStart:idxStop] .== maximum(sig[idxStart:idxStop])) + idxStart-1
                    peakPointsEst[pk] = idx
                    peakLatenciesEst[pk] = ((idx-1)/sampRate)+epochStart
                else
                    peakLatenciesEst[pk] = missing
                    peakPointsEst[pk] = missing
                end
            end
        end
    end

    ###############################
    ## determine trough windows from peak latencies
    troughLatencies = missings(Float64, nPeaks)
    troughPoints = missings(Int, nPeaks) #has to be float to possibly be missing
    troughLatenciesEst = missings(Float64, nPeaks)
    troughPointsEst = missings(Int, nPeaks) #has to be float to possibly be missing
    minTroughLat = missings(Float64, nPeaks)
    maxTroughLat = missings(Float64, nPeaks)
    minTroughLatEst = missings(Float64, nPeaks)
    maxTroughLatEst = missings(Float64, nPeaks)
    for pk=1:nPeaks
        thisMin = minPeakTroughLat[pk]
        thisMax = maxPeakTroughLat[pk]
        ## if ismissing(maxPeakTroughLat[pk]) == true
        ##     if pk+1 <= nPeaks
        ##         thisMax = peakLatencies[pk+1] - peakLatencies[pk]
        ##     end
        ## else
        ##     thisMax = maxPeakTroughLat[pk]
        ## end
        minTroughLat[pk] = peakLatencies[pk]+thisMin
        maxTroughLat[pk] =  peakLatencies[pk]+thisMax
        minTroughLatEst[pk] = peakLatenciesEst[pk]+thisMin
        maxTroughLatEst[pk] =  peakLatenciesEst[pk]+thisMax
        if ismissing(peakLatencies[pk]) == false
            candidateTimes = xtrsTimes[(xtrsTimes .>= peakLatencies[pk]+thisMin) .& (xtrsTimes .<= peakLatencies[pk]+thisMax)]
            candidatePnts = xtrsPnts[(xtrsTimes.>= peakLatencies[pk]+thisMin) .& (xtrsTimes.<= peakLatencies[pk]+thisMax)]
        else
            candidateTimes = (Float64)[]
            candidatePnts = (Int)[]
        end
        if length(candidateTimes) > 0
            #find largest trough
            idx = findall(sig[candidatePnts] .== minimum(sig[candidatePnts]))[1]
            troughLatencies[pk] = candidateTimes[idx]
            troughPoints[pk] = candidatePnts[idx]
            troughLatenciesEst[pk] = candidateTimes[idx]
            troughPointsEst[pk] = candidatePnts[idx]
        else
            troughLatencies[pk] = missing
            troughPoints[pk] = missing
            if ismissing(expPeakLatencies[pk]) == false
                idxStart = round(Int, (minTroughLatEst[pk]-epochStart)*sampRate)+1;
                idxStop = round(Int, (maxTroughLatEst[pk]-epochStart)*sampRate)+1;
                idx = findfirst(sig[idxStart:idxStop] .== minimum(sig[idxStart:idxStop])) + idxStart-1
                troughPointsEst[pk] = idx
                troughLatenciesEst[pk] = ((idx-1)/sampRate)+epochStart
            else
                troughLatenciesEst[pk] = missing
                troughPointsEst[pk] = missing
            end
        end
    end

    #####################################
    ## measure peak and trough amplitudes
    peakAmps = missings(Float64, nPeaks)
    troughAmps = missings(Float64, nPeaks)
    peakAmpsEst = missings(Float64, nPeaks)
    troughAmpsEst = missings(Float64, nPeaks)
    for i=1:nPeaks
        if ismissing(peakPoints[i]) == false
            peakAmps[i] = sig[round(Int, peakPoints[i])]
            peakAmpsEst[i] = sig[round(Int, peakPoints[i])]
        else
            peakAmps[i] = missing
            if ismissing(expPeakLatencies[i]) == false
                #idxStart = round(Int, (expPeakLatencies[i]+winStart[i]-epochStart)*sampRate);
                #idxStop = round(Int, (expPeakLatencies[i]+winStop[i]-epochStart)*sampRate);
                peakAmpsEst[i] = sig[round(Int, peakPointsEst[i])] #maximum(sig[idxStart:idxStop])[1]
            else
                peakAmpsEst[i] = missing
            end
        end

        if ismissing(troughPoints[i]) == false
            troughAmps[i] = sig[round(Int, troughPoints[i])]
            troughAmpsEst[i] = sig[round(Int, troughPoints[i])]
        else
            troughAmps[i] = missing
            if ismissing(minTroughLatEst[i]) == false
                #idxStart = round(Int, (minTroughLat[i]-epochStart)*sampRate);
                #idxStop = round(Int, (maxTroughLat[i]-epochStart)*sampRate);
                troughAmpsEst[i] = sig[round(Int, troughPointsEst[i])] #minimum(sig[idxStart:idxStop])[1]
            else
                troughAmpsEst[i] = missing
            end
        end
    end

    ####################################
    ## measure peak-trough amplitudes
    peakTroughAmps = peakAmps-troughAmps
    peakTroughAmpsEst = peakAmpsEst-troughAmpsEst

    ####################################
    ## store results in a dataframe
    df = DataFrame(wave=waveLabels, peakPoint=peakPoints, troughPoint=troughPoints, peakLatency=peakLatencies, troughLatency=troughLatencies, peakAmp=peakAmps, troughAmp=troughAmps, peakTroughAmp=peakTroughAmps, peakPointEst=peakPointsEst, peakLatencyEst=peakLatenciesEst, peakAmpEst=peakAmpsEst, troughPointEst=troughPointsEst, troughLatencyEst=troughLatenciesEst, troughAmpEst=troughAmpsEst, peakTroughAmpEst=peakTroughAmpsEst,
                   minPeakLat=expPeakLatencies+winStart,
                   maxPeakLat=expPeakLatencies+winStop,
                   minTroughLat=minTroughLat,
                   maxTroughLat=maxTroughLat)

    return df

end

##difference from calcPeaksInflFromExpectedLatencies is that this function only uses
## down-going inflections to find peaks. Also, the inflection point with the lowest absolute first derivative is chosen.
function calcPeaksInflFromExpectedLatencies2(sig::Union{AbstractMatrix{T}, AbstractVector{T}}, expPeakLatencies::Union{AbstractVector{S}, AbstractVector{Union{S, Missing}}}, winStart::AbstractVector{P}, winStop::AbstractVector{Q}, sampRate::Real; epochStart::Real=0, minPeakTroughLat::AbstractVector{G}=[0.25, 0.12, 0.25, 0.12, 0.25]./1000, maxPeakTroughLat::AbstractVector{H}=[1, 1, 1, 0.75, 2]./1000, waveLabels::AbstractVector{String} = ["I", "II", "III", "IV", "V"]) where {T<:Real, S<:Real, P<:Real, Q<:Real, G<:Real, H<:Real}
    nPeaks = length(expPeakLatencies)
    xpksPnts, xpksTimes = findPeaks(sig, sampRate, epochStart=epochStart) #find all peaks in the waveform
    xtrsPnts, xtrsTimes = findTroughs(sig, sampRate, epochStart=epochStart) #find all troughs in the waveform
    xinflPnts, xinflTimes = findInflections(sig, sampRate, epochStart=epochStart) #find all inflections in the waveform
    xinflPntsUp, xinflTimesUp = findInflectionsUp(sig, sampRate, epochStart=epochStart) #find all inflections in the waveform
    xinflPntsDown, xinflTimesDown = findInflectionsDown(sig, sampRate, epochStart=epochStart) #find all inflections in the waveform

    #############################
    ## find peaks
    peakLatencies = missings(Float64, nPeaks)
    peakPoints = missings(Int, nPeaks) #has to be float to possibly be missing
    peakLatenciesEst = missings(Float64, nPeaks) #in case a peak can't be found fall back to finding the maximum value
    peakPointsEst = missings(Int, nPeaks)
    for pk=1:nPeaks
        if ismissing(expPeakLatencies[pk]) == false
            candidateTimes = xpksTimes[(xpksTimes .>= expPeakLatencies[pk]+winStart[pk]) .& (xpksTimes .<= expPeakLatencies[pk]+winStop[pk])]
            candidatePnts = xpksPnts[(xpksTimes.>=expPeakLatencies[pk]+winStart[pk]) .& (xpksTimes.<=expPeakLatencies[pk]+winStop[pk])]
        else
            candidateTimes = (Float64)[]
            candidatePnts = (Int)[]
        end
        if length(candidateTimes) > 0
            #find largest peak
            idx = findall(sig[candidatePnts] .== maximum(sig[candidatePnts]))[1]
            peakLatencies[pk] = candidateTimes[idx]
            peakPoints[pk] = candidatePnts[idx]
            peakLatenciesEst[pk] = candidateTimes[idx]
            peakPointsEst[pk] = candidatePnts[idx]
        else
            ##try inflection points
            if ismissing(expPeakLatencies[pk]) == false
                candidateTimes = xinflTimesDown[(xinflTimesDown .>= expPeakLatencies[pk]+winStart[pk]) .& (xinflTimesDown .<= expPeakLatencies[pk]+winStop[pk])]
                candidatePnts = xinflPntsDown[(xinflTimesDown.>=expPeakLatencies[pk]+winStart[pk]) .& (xinflTimesDown.<=expPeakLatencies[pk]+winStop[pk])]
            else
                candidateTimes = (Float64)[]
                candidatePnts = (Int)[]
            end
            if length(candidateTimes) > 0
                ###find largest peak
                ##idx = findall(sig[candidatePnts] .== maximum(sig[candidatePnts]))[1]
                ## pick point where first derivative has minimum absolute value
                dy = diff(sig)
                idx = findall(abs.(dy[candidatePnts.-1]) .== minimum.(abs.(dy[candidatePnts.-1])))[1]#sig[candidatePnts] .== maximum(sig[candidatePnts]))[1]
                peakLatencies[pk] = candidateTimes[idx]
                peakPoints[pk] = candidatePnts[idx]
                peakLatenciesEst[pk] = candidateTimes[idx]
                peakPointsEst[pk] = candidatePnts[idx]
            else
                peakLatencies[pk] = missing
                peakPoints[pk] = missing
                if ismissing(expPeakLatencies[pk]) == false
                    idxStart = round(Int, (expPeakLatencies[pk]+winStart[pk]-epochStart)*sampRate)+1;
                    idxStop = round(Int, (expPeakLatencies[pk]+winStop[pk]-epochStart)*sampRate)+1;
                    idx = findfirst(sig[idxStart:idxStop] .== maximum(sig[idxStart:idxStop])) + idxStart-1
                    peakPointsEst[pk] = idx
                    peakLatenciesEst[pk] = ((idx-1)/sampRate)+epochStart
                else
                    peakLatenciesEst[pk] = missing
                    peakPointsEst[pk] = missing
                end
            end
        end
    end

    ###############################
    ## determine trough windows from peak latencies
    troughLatencies = missings(Float64, nPeaks)
    troughPoints = missings(Int, nPeaks) #has to be float to possibly be missing
    troughLatenciesEst = missings(Float64, nPeaks)
    troughPointsEst = missings(Int, nPeaks) #has to be float to possibly be missing
    minTroughLat = missings(Float64, nPeaks)
    maxTroughLat = missings(Float64, nPeaks)
    minTroughLatEst = missings(Float64, nPeaks)
    maxTroughLatEst = missings(Float64, nPeaks)
    for pk=1:nPeaks
        thisMin = minPeakTroughLat[pk]
        thisMax = maxPeakTroughLat[pk]
        ## if ismissing(maxPeakTroughLat[pk]) == true
        ##     if pk+1 <= nPeaks
        ##         thisMax = peakLatencies[pk+1] - peakLatencies[pk]
        ##     end
        ## else
        ##     thisMax = maxPeakTroughLat[pk]
        ## end
        minTroughLat[pk] = peakLatencies[pk]+thisMin
        maxTroughLat[pk] =  peakLatencies[pk]+thisMax
        minTroughLatEst[pk] = peakLatenciesEst[pk]+thisMin
        maxTroughLatEst[pk] =  peakLatenciesEst[pk]+thisMax
        if ismissing(peakLatencies[pk]) == false
            candidateTimes = xtrsTimes[(xtrsTimes .>= peakLatencies[pk]+thisMin) .& (xtrsTimes .<= peakLatencies[pk]+thisMax)]
            candidatePnts = xtrsPnts[(xtrsTimes.>= peakLatencies[pk]+thisMin) .& (xtrsTimes.<= peakLatencies[pk]+thisMax)]
        else
            candidateTimes = (Float64)[]
            candidatePnts = (Int)[]
        end
        if length(candidateTimes) > 0
            #find largest trough
            idx = findall(sig[candidatePnts] .== minimum(sig[candidatePnts]))[1]
            troughLatencies[pk] = candidateTimes[idx]
            troughPoints[pk] = candidatePnts[idx]
            troughLatenciesEst[pk] = candidateTimes[idx]
            troughPointsEst[pk] = candidatePnts[idx]
        else
            troughLatencies[pk] = missing
            troughPoints[pk] = missing
            if ismissing(expPeakLatencies[pk]) == false
                idxStart = round(Int, (minTroughLatEst[pk]-epochStart)*sampRate)+1;
                idxStop = round(Int, (maxTroughLatEst[pk]-epochStart)*sampRate)+1;
                idx = findfirst(sig[idxStart:idxStop] .== minimum(sig[idxStart:idxStop])) + idxStart-1
                troughPointsEst[pk] = idx
                troughLatenciesEst[pk] = ((idx-1)/sampRate)+epochStart
            else
                troughLatenciesEst[pk] = missing
                troughPointsEst[pk] = missing
            end
        end
    end

    #####################################
    ## measure peak and trough amplitudes
    peakAmps = missings(Float64, nPeaks)
    troughAmps = missings(Float64, nPeaks)
    peakAmpsEst = missings(Float64, nPeaks)
    troughAmpsEst = missings(Float64, nPeaks)
    for i=1:nPeaks
        if ismissing(peakPoints[i]) == false
            peakAmps[i] = sig[round(Int, peakPoints[i])]
            peakAmpsEst[i] = sig[round(Int, peakPoints[i])]
        else
            peakAmps[i] = missing
            if ismissing(expPeakLatencies[i]) == false
                #idxStart = round(Int, (expPeakLatencies[i]+winStart[i]-epochStart)*sampRate);
                #idxStop = round(Int, (expPeakLatencies[i]+winStop[i]-epochStart)*sampRate);
                peakAmpsEst[i] = sig[round(Int, peakPointsEst[i])] #maximum(sig[idxStart:idxStop])[1]
            else
                peakAmpsEst[i] = missing
            end
        end

        if ismissing(troughPoints[i]) == false
            troughAmps[i] = sig[round(Int, troughPoints[i])]
            troughAmpsEst[i] = sig[round(Int, troughPoints[i])]
        else
            troughAmps[i] = missing
            if ismissing(minTroughLatEst[i]) == false
                #idxStart = round(Int, (minTroughLat[i]-epochStart)*sampRate);
                #idxStop = round(Int, (maxTroughLat[i]-epochStart)*sampRate);
                troughAmpsEst[i] = sig[round(Int, troughPointsEst[i])] #minimum(sig[idxStart:idxStop])[1]
            else
                troughAmpsEst[i] = missing
            end
        end
    end

    ####################################
    ## measure peak-trough amplitudes
    peakTroughAmps = peakAmps-troughAmps
    peakTroughAmpsEst = peakAmpsEst-troughAmpsEst

    ####################################
    ## store results in a dataframe
    df = DataFrame(wave=waveLabels, peakPoint=peakPoints, troughPoint=troughPoints, peakLatency=peakLatencies, troughLatency=troughLatencies, peakAmp=peakAmps, troughAmp=troughAmps, peakTroughAmp=peakTroughAmps, peakPointEst=peakPointsEst, peakLatencyEst=peakLatenciesEst, peakAmpEst=peakAmpsEst, troughPointEst=troughPointsEst, troughLatencyEst=troughLatenciesEst, troughAmpEst=troughAmpsEst, peakTroughAmpEst=peakTroughAmpsEst,
                   minPeakLat=expPeakLatencies+winStart,
                   maxPeakLat=expPeakLatencies+winStop,
                   minTroughLat=minTroughLat,
                   maxTroughLat=maxTroughLat)

    return df

end

#difference from calcPeaksInflFromExpectedLatencies2 is that there are further constraint on search window (useful for example
# to make sure wave at a lower level has delayed latency compared to wave previously measure at a higher level
function calcPeaksInflFromExpectedLatencies3(sig::Union{AbstractMatrix{T}, AbstractVector{T}},
                                             expPeakLatencies::Union{AbstractVector{S}, AbstractVector{Union{S, Missing}}},
                                             winStart::AbstractVector{P}, winStop::AbstractVector{Q},
                                             sampRate::Real; epochStart::Real=0,
                                             minPeakTroughLat::AbstractVector{G}=[0.25, 0.12, 0.25, 0.12, 0.25]./1000,
                                             maxPeakTroughLat::AbstractVector{H}=[1, 1, 1, 0.75, 2]./1000,
                                             minAbsLat::Union{AbstractVector{S}, AbstractVector{Union{R, Missing}}}=[missing for i=1:5],
                                             maxAbsLat::Union{AbstractVector{S}, AbstractVector{Union{W, Missing}}}=[missing for i=1:5],
                                             waveLabels::AbstractVector{String} = ["I", "II", "III", "IV", "V"]) where {T<:Real, S<:Real, P<:Real, Q<:Real, G<:Real, H<:Real, R<:Real, W<:Real}
    nPeaks = length(expPeakLatencies)
    xpksPnts, xpksTimes = findPeaks(sig, sampRate, epochStart=epochStart) #find all peaks in the waveform
    xtrsPnts, xtrsTimes = findTroughs(sig, sampRate, epochStart=epochStart) #find all troughs in the waveform
    xinflPnts, xinflTimes = findInflections(sig, sampRate, epochStart=epochStart) #find all inflections in the waveform
    xinflPntsUp, xinflTimesUp = findInflectionsUp(sig, sampRate, epochStart=epochStart) #find all inflections in the waveform
    xinflPntsDown, xinflTimesDown = findInflectionsDown(sig, sampRate, epochStart=epochStart) #find all inflections in the waveform

    #############################
    ## find peaks
    peakLatencies = missings(Float64, nPeaks)
    peakPoints = missings(Int, nPeaks) #has to be float to possibly be missing
    peakLatenciesEst = missings(Float64, nPeaks) #in case a peak can't be found fall back to finding the maximum value
    peakPointsEst = missings(Int, nPeaks)

    searchWinStart = missings(Float64, nPeaks)
    searchWinStop = missings(Float64, nPeaks)
    for pk=1:nPeaks
        searchWinStart[pk] = max(expPeakLatencies[pk]+winStart[pk], ifelse(ismissing(minAbsLat[pk]), expPeakLatencies[pk]+winStart[pk], minAbsLat[pk]))
        searchWinStop[pk] = min(expPeakLatencies[pk]+winStop[pk], ifelse(ismissing(maxAbsLat[pk]), expPeakLatencies[pk]+winStop[pk], maxAbsLat[pk]))
    end
    for pk=1:nPeaks
        if ismissing(expPeakLatencies[pk]) == false
            candidateTimes = xpksTimes[(xpksTimes .>= searchWinStart[pk]) .& (xpksTimes .<= searchWinStop[pk])]
            candidatePnts = xpksPnts[(xpksTimes.>=searchWinStart[pk]) .& (xpksTimes.<=searchWinStop[pk])]
        else
            candidateTimes = (Float64)[]
            candidatePnts = (Int)[]
        end
        if length(candidateTimes) > 0
            #find largest peak
            idx = findall(sig[candidatePnts] .== maximum(sig[candidatePnts]))[1]
            peakLatencies[pk] = candidateTimes[idx]
            peakPoints[pk] = candidatePnts[idx]
            peakLatenciesEst[pk] = candidateTimes[idx]
            peakPointsEst[pk] = candidatePnts[idx]
        else
            ##try inflection points
            if ismissing(expPeakLatencies[pk]) == false
                candidateTimes = xinflTimesDown[(xinflTimesDown .>= searchWinStart[pk]) .& (xinflTimesDown .<= searchWinStop[pk])]
                candidatePnts = xinflPntsDown[(xinflTimesDown.>=searchWinStart[pk]) .& (xinflTimesDown.<=searchWinStop[pk])]
            else
                candidateTimes = (Float64)[]
                candidatePnts = (Int)[]
            end
            if length(candidateTimes) > 0
                ###find largest peak
                ##idx = findall(sig[candidatePnts] .== maximum(sig[candidatePnts]))[1]
                ## pick point where first derivative has minimum absolute value
                dy = diff(sig)
                idx = findall(abs.(dy[candidatePnts.-1]) .== minimum.(abs.(dy[candidatePnts.-1])))[1]#sig[candidatePnts] .== maximum(sig[candidatePnts]))[1]
                peakLatencies[pk] = candidateTimes[idx]
                peakPoints[pk] = candidatePnts[idx]
                peakLatenciesEst[pk] = candidateTimes[idx]
                peakPointsEst[pk] = candidatePnts[idx]
            else
                peakLatencies[pk] = missing
                peakPoints[pk] = missing
                if ismissing(expPeakLatencies[pk]) == false
                    idxStart = round(Int, (searchWinStart[pk]-epochStart)*sampRate)+1;
                    idxStop = round(Int, (searchWinStop[pk]-epochStart)*sampRate)+1;
                    idx = findfirst(sig[idxStart:idxStop] .== maximum(sig[idxStart:idxStop])) + idxStart-1
                    peakPointsEst[pk] = idx
                    peakLatenciesEst[pk] = ((idx-1)/sampRate)+epochStart
                else
                    peakLatenciesEst[pk] = missing
                    peakPointsEst[pk] = missing
                end
            end
        end
    end

    ###############################
    ## determine trough windows from peak latencies
    troughLatencies = missings(Float64, nPeaks)
    troughPoints = missings(Int, nPeaks) #has to be float to possibly be missing
    troughLatenciesEst = missings(Float64, nPeaks)
    troughPointsEst = missings(Int, nPeaks) #has to be float to possibly be missing
    minTroughLat = missings(Float64, nPeaks)
    maxTroughLat = missings(Float64, nPeaks)
    minTroughLatEst = missings(Float64, nPeaks)
    maxTroughLatEst = missings(Float64, nPeaks)
    for pk=1:nPeaks
        thisMin = minPeakTroughLat[pk]
        thisMax = maxPeakTroughLat[pk]
        ## if ismissing(maxPeakTroughLat[pk]) == true
        ##     if pk+1 <= nPeaks
        ##         thisMax = peakLatencies[pk+1] - peakLatencies[pk]
        ##     end
        ## else
        ##     thisMax = maxPeakTroughLat[pk]
        ## end
        minTroughLat[pk] = peakLatencies[pk]+thisMin
        maxTroughLat[pk] =  peakLatencies[pk]+thisMax
        minTroughLatEst[pk] = peakLatenciesEst[pk]+thisMin
        maxTroughLatEst[pk] =  peakLatenciesEst[pk]+thisMax
        if ismissing(peakLatencies[pk]) == false
            candidateTimes = xtrsTimes[(xtrsTimes .>= peakLatencies[pk]+thisMin) .& (xtrsTimes .<= peakLatencies[pk]+thisMax)]
            candidatePnts = xtrsPnts[(xtrsTimes.>= peakLatencies[pk]+thisMin) .& (xtrsTimes.<= peakLatencies[pk]+thisMax)]
        else
            candidateTimes = (Float64)[]
            candidatePnts = (Int)[]
        end
        if length(candidateTimes) > 0
            #find largest trough
            idx = findall(sig[candidatePnts] .== minimum(sig[candidatePnts]))[1]
            troughLatencies[pk] = candidateTimes[idx]
            troughPoints[pk] = candidatePnts[idx]
            troughLatenciesEst[pk] = candidateTimes[idx]
            troughPointsEst[pk] = candidatePnts[idx]
        else
            troughLatencies[pk] = missing
            troughPoints[pk] = missing
            if ismissing(expPeakLatencies[pk]) == false
                idxStart = round(Int, (minTroughLatEst[pk]-epochStart)*sampRate)+1;
                idxStop = round(Int, (maxTroughLatEst[pk]-epochStart)*sampRate)+1;
                idx = findfirst(sig[idxStart:idxStop] .== minimum(sig[idxStart:idxStop])) + idxStart-1
                troughPointsEst[pk] = idx
                troughLatenciesEst[pk] = ((idx-1)/sampRate)+epochStart
            else
                troughLatenciesEst[pk] = missing
                troughPointsEst[pk] = missing
            end
        end
    end

    #####################################
    ## measure peak and trough amplitudes
    peakAmps = missings(Float64, nPeaks)
    troughAmps = missings(Float64, nPeaks)
    peakAmpsEst = missings(Float64, nPeaks)
    troughAmpsEst = missings(Float64, nPeaks)
    for i=1:nPeaks
        if ismissing(peakPoints[i]) == false
            peakAmps[i] = sig[round(Int, peakPoints[i])]
            peakAmpsEst[i] = sig[round(Int, peakPoints[i])]
        else
            peakAmps[i] = missing
            if ismissing(expPeakLatencies[i]) == false
                #idxStart = round(Int, (expPeakLatencies[i]+winStart[i]-epochStart)*sampRate);
                #idxStop = round(Int, (expPeakLatencies[i]+winStop[i]-epochStart)*sampRate);
                peakAmpsEst[i] = sig[round(Int, peakPointsEst[i])] #maximum(sig[idxStart:idxStop])[1]
            else
                peakAmpsEst[i] = missing
            end
        end

        if ismissing(troughPoints[i]) == false
            troughAmps[i] = sig[round(Int, troughPoints[i])]
            troughAmpsEst[i] = sig[round(Int, troughPoints[i])]
        else
            troughAmps[i] = missing
            if ismissing(minTroughLatEst[i]) == false
                #idxStart = round(Int, (minTroughLat[i]-epochStart)*sampRate);
                #idxStop = round(Int, (maxTroughLat[i]-epochStart)*sampRate);
                troughAmpsEst[i] = sig[round(Int, troughPointsEst[i])] #minimum(sig[idxStart:idxStop])[1]
            else
                troughAmpsEst[i] = missing
            end
        end
    end

    ####################################
    ## measure peak-trough amplitudes
    peakTroughAmps = peakAmps-troughAmps
    peakTroughAmpsEst = peakAmpsEst-troughAmpsEst

    ####################################
    ## store results in a dataframe
    df = DataFrame(wave=waveLabels, peakPoint=peakPoints, troughPoint=troughPoints, peakLatency=peakLatencies, troughLatency=troughLatencies, peakAmp=peakAmps, troughAmp=troughAmps, peakTroughAmp=peakTroughAmps, peakPointEst=peakPointsEst, peakLatencyEst=peakLatenciesEst, peakAmpEst=peakAmpsEst, troughPointEst=troughPointsEst, troughLatencyEst=troughLatenciesEst, troughAmpEst=troughAmpsEst, peakTroughAmpEst=peakTroughAmpsEst,
                   minPeakLat=expPeakLatencies+winStart,
                   maxPeakLat=expPeakLatencies+winStop,
                   minTroughLat=minTroughLat,
                   maxTroughLat=maxTroughLat)

    return df

end

function findInflectionsUp(y::Union{AbstractMatrix{T}, AbstractVector{T}}, sampRate::Real; epochStart::Real=0) where {T<:Real}
    inflPntsUp = (Int)[]

    y = vec(y)
    dy = diff(y)
    ddy = diff(dy)
    sgnch = ddy[1:end-1].*ddy[2:end]
    for i=1:length(sgnch)
        if signbit(sgnch[i]) == true
            if dy[i+1] > 0
                push!(inflPntsUp, i+2)
            end
        end
    end

    inflTimesUp = zeros(length(inflPntsUp))

    for i=1:length(inflPntsUp)
        inflTimesUp[i] = (inflPntsUp[i]-1)/sampRate
    end

    inflTimesUp = inflTimesUp .+ epochStart

    return inflPntsUp, inflTimesUp

end

function findInflectionsDown(y::Union{AbstractMatrix{T}, AbstractVector{T}}, sampRate::Real; epochStart::Real=0) where {T<:Real}
    inflPntsDown = (Int)[]

    y = vec(y)
    dy = diff(y)
    ddy = diff(dy)
    sgnch = ddy[1:end-1].*ddy[2:end]
    for i=1:length(sgnch)
        if signbit(sgnch[i]) == true
            if dy[i+1] < 0
                push!(inflPntsDown, i+2)
            end
        end
    end

    inflTimesDown = zeros(length(inflPntsDown))

    for i=1:length(inflPntsDown)
        inflTimesDown[i] = (inflPntsDown[i]-1)/sampRate
    end

    inflTimesDown = inflTimesDown .+ epochStart

    return inflPntsDown, inflTimesDown

end

## function calcPeakPoints(sig, pkTimes, delay, sampRate; epochStart=0)

##     xpksPnts, xpksTimes = findPeaks(sig, sampRate, epochStart=epochStart)
##     xtrsPnts, xtrsTimes = findTroughs(sig, sampRate, epochStart=epochStart)

##     peakTimes = (Float64)[]
##     peakPoints = (Float64)[]
##     for pk=1:length(pkTimes)
##         fooTimes = xpksTimes[(xpksTimes .>= pkTimes[pk]+delay[1]) .& (xpksTimes .<= pkTimes[pk]+delay[2])]
##         fooPnts = xpksPnts[(xpksTimes.>=pkTimes[pk]+delay[1]) .& (xpksTimes.<=pkTimes[pk]+delay[2])]
##         if length(fooTimes) > 0
##             #find largest peak
##             idx = findall(sig[fooPnts] .== maximum(sig[fooPnts]))[1]
##             push!(peakTimes, fooTimes[idx])
##             push!(peakPoints, fooPnts[idx])
##         else
##             push!(peakTimes, missing)
##             push!(peakPoints, missing)
##         end
##     end

##     minPeakTroughLat = [0.25/1000, 0.12/1000, 0.25/1000, 0.12/1000, 0.25/1000]
##     maxPeakTroughLatDiff = [0.75/1000, missing, 0.75/1000, missing, 1.5/1000]
##     maxPeakTroughLat = minPeakTroughLat + maxPeakTroughLatDiff

##     troughTimes = (Float64)[]
##     troughPoints = (Float64)[]
##     for pk=1:length(pkTimes)
##         thisMin = minPeakTroughLat[pk]
##         if isnan(maxPeakTroughLat[pk]) == true
##             if pk+1 <= 5
##                 thisMax = peakTimes[pk+1] - peakTimes[pk]
##             end
##         else
##             thisMax = maxPeakTroughLat[pk]
##         end
##         fooTimes = xtrsTimes[(xtrsTimes .>= peakTimes[pk]+thisMin) .& (xtrsTimes .<= peakTimes[pk]+thisMax)]
##         fooPnts = xtrsPnts[(xtrsTimes.>= peakTimes[pk]+thisMin) .& (xtrsTimes.<= peakTimes[pk]+thisMax)]
##         if length(fooTimes) > 0
##             #find largest trough
##             idx = findall(sig[fooPnts] .== minimum(sig[fooPnts]))[1]
##             push!(troughTimes, fooTimes[idx])
##             push!(troughPoints, fooPnts[idx])
##         else
##             push!(troughTimes, missing)
##             push!(troughPoints, missing)
##         end
##     end

##     peakLatencies = peakTimes
##     troughLatencies = troughTimes

##     peakAmps = zeros(length(peakPoints))
##     troughAmps = zeros(length(peakPoints))
##     for i=1:length(peakAmps)
##         if isnan(peakPoints[i]) == false
##             peakAmps[i] = sig[1, round(Int, peakPoints[i])]
##         else
##             peakAmps[i] = missing
##         end

##         if isnan(troughPoints[i]) == false
##             troughAmps[i] = sig[1, round(Int, troughPoints[i])]
##         else
##             troughAmps[i] = missing
##         end
##     end

##     peakTroughAmps = peakAmps-troughAmps
##     prms = []## Dict("minPeakLat" => [minPeakILat, minPeakIILat, minPeakIIILat, minPeakIVLat, minPeakVLat],
##              ##    "maxPeakLat" => [maxPeakILat, maxPeakIILat, maxPeakIIILat, maxPeakIVLat, maxPeakVLat],
##              ##    "minTroughLat" => [minTroughILat, minTroughIILat, minTroughIIILat, minTroughIVLat, minTroughVLat],
##              ##    "maxTroughLat" => [maxTroughILat, maxTroughIILat, maxTroughIIILat, maxTroughIVLat, maxTroughVLat])

##     return peakPoints, troughPoints, peakLatencies, troughLatencies, peakAmps, troughAmps, peakTroughAmps, prms

## end


text:
"""
    mutable struct RkMatrix{T}

Representation of a rank `r` matrix `M` in outer product format `M =
A*adjoint(B)` where `A` has size `m × r` and `B` has size `n × r`.

The internal representation stores `A` and `B`, but `R.Bt` or `R.At` can be used
to get the respective adjoints.
"""
mutable struct RkMatrix{T} <: AbstractMatrix{T}
    A::Matrix{T}
    B::Matrix{T}
    function RkMatrix(A::Matrix{T},B::Matrix{T}) where {T}
        @assert size(A,2) == size(B,2) "second dimension of `A` and `B` must match"
        m,r = size(A)
        n  = size(B,1)
        if  r*(m+n) >= m*n
            @debug "Inefficient RkMatrix:" size(A) size(B)
            # error("Inefficient RkMatrix")
        end
        new{T}(A,B)
    end
end
RkMatrix(A,B) = RkMatrix(promote(A,B)...)

function Base.getindex(rmat::RkMatrix,i::Int,j::Int)
    msg = """
    method `getindex(::AbstractHMatrix,args...)` has been disabled to avoid
    performance pitfalls. Unless you made an explicit call to `getindex`, this
    error usually means that a linear algebra routine involving an
    `AbstractHMatrix` has fallen back to a generic implementation.
    """
    if ALLOW_GETINDEX[]
        r = rank(rmat)
        acc = zero(eltype(rmat))
        for k in 1:r
            acc += rmat.A[i,k] * conj(rmat.B[j,k])
        end
    else
        error(msg)
    end
    return acc
end

# some "fast" ways of computing a column of R and row of ajoint(R)
Base.getindex(R::RkMatrix, ::Colon, j::Int) = getcol(R,j)

"""
    getcol!(col,M::AbstractMatrix,j)

Fill the entries of `col` with column `j` of `M`.
"""
function getcol!(col,R::RkMatrix,j::Int)
    mul!(col,R.A,conj(view(R.B,j,:)))
end

"""
    getcol(M::AbstractMatrix,j)

Return a vector containing the `j`-th column of `M`.
"""
function getcol(R::RkMatrix,j::Int)
    m = size(R,1)
    T = eltype(R)
    col = zeros(T,m)
    getcol!(col,R,j)
end

function getcol!(col,Ra::Adjoint{<:Any,<:RkMatrix},j::Int)
    R = parent(Ra)
    mul!(col,R.B,conj(view(R.A,j,:)))
end

function getcol(Ra::Adjoint{<:Any,<:RkMatrix},j::Int)
    m = size(Ra,1)
    T = eltype(Ra)
    col = zeros(T,m)
    getcol!(col,Ra,j)
end

Base.getindex(Ra::Adjoint{<:Any,<:RkMatrix}, ::Colon, j::Int) = getcol(Ra,j)

"""
    RkMatrix(A::Vector{<:Vector},B::Vector{<:Vector})

Construct an `RkMatrix` from a vector of vectors. Assumes that `length(A) ==
length(B)`, which determines the rank, and that all vectors in `A` (resp. `B`) have the same length `m`
(resp. `n`).
"""
function RkMatrix(_A::Vector{V},_B::Vector{V}) where {V<:AbstractVector}
    T   = eltype(V)
    @assert length(_A) == length(_B)
    k   = length(_A)
    m   = length(first(_A))
    n   = length(first(_B))
    A   = Matrix{T}(undef,m,k)
    B  = Matrix{T}(undef,n,k)
    for i in 1:k
        copyto!(view(A,:,i),_A[i])
        copyto!(view(B,:,i),_B[i])
    end
    return RkMatrix(A,B)
end

Base.eltype(::RkMatrix{T}) where {T} = T
Base.size(rmat::RkMatrix)                                        = (size(rmat.A,1), size(rmat.B,1))
Base.length(rmat::RkMatrix)                                      = prod(size(rmat))
Base.isapprox(rmat::RkMatrix,B::AbstractArray,args...;kwargs...) = isapprox(Matrix(rmat),B,args...;kwargs...)

rank(M::RkMatrix) = size(M.A,2)

function Base.getproperty(R::RkMatrix,s::Symbol)
    if  s == :Bt
        return adjoint(R.B)
    elseif  s == :At
        return adjoint(R.A)
    else
        return getfield(R,s)
    end
end

"""
    hcat(M1::RkMatrix,M2::RkMatrix)

Concatenated `M1` and `M2` horizontally to produce a new `RkMatrix` of rank
`rank(M1)+rank(M2)`.
"""
function Base.hcat(M1::RkMatrix{T},M2::RkMatrix{T}) where {T}
    m,n  = size(M1)
    s,t  = size(M2)
    (m == s) || throw(ArgumentError("number of rows of each array must match: got  ($m,$s)"))
    r1   = size(M1.A,2)
    r2   = size(M2.A,2)
    A    = hcat(M1.A,M2.A)
    B1   = vcat(M1.B,zeros(T,t,r1))
    B2   = vcat(zeros(T,n,r2),M2.B)
    B    = hcat(B1,B2)
    return RkMatrix(A,B)
end

"""
    vcat(M1::RkMatrix,M2::RkMatrix)

Concatenated `M1` and `M2` vertically to produce a new `RkMatrix` of rank
`rank(M1)+rank(M2)`
"""
function Base.vcat(M1::RkMatrix{T},M2::RkMatrix{T}) where {T}
    m,n  = size(M1)
    s,t  = size(M2)
    n == t || throw(ArgumentError("number of columns of each array must match (got  ($n,$t))"))
    r1   = size(M1.A,2)
    r2   = size(M2.A,2)
    A1   = vcat(M1.A,zeros(T,s,r1))
    A2   = vcat(zeros(T,m,r2),M2.A)
    A    = hcat(A1,A2)
    B    = hcat(M1.B,M2.B)
    return RkMatrix(A,B)
end

Base.copy(R::RkMatrix) = RkMatrix(copy(R.A),copy(R.B))

function Base.Matrix(R::RkMatrix{<:Number})
    Matrix(R.A*R.Bt)
end
function Base.Matrix(R::RkMatrix{<:SMatrix})
    # collect must be used when we have a matrix of `SMatrix` because of this issue:
    # https://github.com/JuliaArrays/StaticArrays.jl/issues/966#issuecomment-943679214
    R.A*collect(R.Bt)
end

"""
    RkMatrix(F::SVD)

Construct an [`RkMatrix`](@ref) from an `SVD` factorization.
"""
function RkMatrix(F::SVD)
    A  = F.U*Diagonal(F.S)
    B  = copy(F.V)
    return RkMatrix(A,B)
end

"""
    RkMatrix(M::Matrix)

Construct an [`RkMatrix`](@ref) from a `Matrix` by passing through the full
`svd` of `M`.
"""
function RkMatrix(M::Matrix)
    F = svd(M)
    RkMatrix(F)
end

"""
    num_stored_elements(R::RkMatrix)

The number of entries stored in the representation. Note that this is *not*
`length(R)`.
"""
num_stored_elements(R::RkMatrix)        = size(R.A,2)*(sum(size(R)))

"""
    compression_ratio(R::RkMatrix)

The ratio of the uncompressed size of `R` to its compressed size in outer
product format.
"""
compression_ratio(R::RkMatrix) = prod(size(R)) / num_stored_elements(R)

# FIXME: to support scalar and tensorial problems, we currently allow for T
# to be something other than a plain number. If that is the case, we
# implement a (slow) multiplication algorithm by hand to circumvent a
# problem in LinearAlgebra for the generic mulplication mul!(C,A,B,a,b) when
# C and B are a vectors of static matrices, and A is a matrix of static
# matrices. Should eventually be removed.
function mul!(C::AbstractVector,Rk::RkMatrix{T},F::AbstractVector,a::Number,b::Number) where {T<:SMatrix}
    m,n = size(Rk)
    r   = rank(Rk)
    tmp = Rk.Bt*F
    rmul!(C,b)
    for k in 1:r
        tmp[k] *= a
        for i in 1:m
            C[i] += Rk.A[i,k]*tmp[k]
        end
    end
    return C
end

ext: jl

max_stars_repo_name: WaveProp/HMatrices

max_forks_repo_name: WaveProp/HMatrices


text:
module ShortestPaths
using SparseArrays: sparse
using LightGraphs
using LightGraphs: AbstractGraph, AbstractEdge
using LightGraphs.SimpleGraphs: AbstractSimpleGraph

using DataStructures:PriorityQueue, enqueue!, dequeue!

# TODO: figure out how we keep environmental params.
# struct LGEnvironment
#     threaded::Bool
#     parallel::Bool
#     LGEnvironment() = new(false, false)
# end

abstract type AbstractGraphResult end
abstract type AbstractGraphAlgorithm end

"""
    ShortestPathResult <: AbstractGraphResult

Concrete subtypes of `ShortestPathResult` contain the
results of a shortest-path calculation using a specific
[`ShortestPathAlgorithm`](@ref).

In general, the fields in these structs should not be
accessed directly; use the [`dists`](@ref) and
[`paths`](@ref) functions to obtain the results of the
calculation.
"""
abstract type ShortestPathResult <: AbstractGraphResult end

"""
    ShortestPathAlgorithm <: AbstractGraphAlgorithm

Concrete subtypes of `ShortestPathAlgorithm` are used to specify
the type of shortest path calculation used by [`shortest_paths`](@ref).
Some concrete subtypes (most notably [`Dijkstra`](@ref) have fields
that specify algorithm parameters.

See [`AStar`](@ref), [`BellmanFord`](@ref), [`BFS`](@ref),
[`DEspopoPape`](@ref), [`Dijkstra`](@ref), [`FloydWarshall`](@ref),
[`Johnson`](@ref), and [`SPFA`](@ref) for specific requirements and
usage details.
"""
abstract type ShortestPathAlgorithm <: AbstractGraphAlgorithm end

include("astar.jl")
include("bellman-ford.jl")
include("bfs.jl")
include("desopo-pape.jl")
include("dijkstra.jl")
include("floyd-warshall.jl")
include("johnson.jl")
include("spfa.jl")

################################
# Shortest Paths via algorithm #
################################
# if we don't pass in distances but specify an algorithm, use weights.
"""
    shortest_paths(g, s, distmx, alg)
    shortest_paths(g, s, t, alg)
    shortest_paths(g, s, alg)
    shortest_paths(g, s)
    shortest_paths(g)

Return a `ShortestPathResult` that allows construction of the shortest path
between sets of vertices in graph `g`. Depending on the algorithm specified,
other information may be required: (e.g., a distance matrix `distmx`, and/or
a target vertex `t`). Some algorithms will accept multiple source vertices
`s`; algorithms that do not accept any source vertex `s` produce all-pairs
shortest paths.

See `ShortestPathAlgorithm` for more details on the algorithm specifications.

### Implementation Notes
The elements of `distmx` may be of any type that has a [Total Ordering](https://en.m.wikipedia.org/wiki/Total_order)
and valid comparator, `zero` and `typemax` functions. Concretely, this means that
distance matrices containing complex numbers are invalid.

### Examples
```
g = path_graph(4)
w = zeros(4, 4)
for i in 1:3
    w[i, i+1] = 1.0
    w[i+1, i] = 1.0
end

s1 = shortest_paths(g)                   # `alg` defaults to `FloydWarshall`
s2 = shortest_paths(g, 1)                # `alg` defaults to `BFS`
s3 = shortest_paths(g, 1, w)             # `alg` defaults to `Dijkstra`
s4 = shortest_paths(g, 1, BellmanFord())
s5 = shortest_paths(g, 1, w, DEsopoPape())
```
"""
shortest_paths(g::AbstractGraph, s, alg::ShortestPathAlgorithm) =
    shortest_paths(g, s, weights(g), alg)

# If we don't specify an algorithm AND there are no dists, use BFS.
shortest_paths(g::AbstractGraph{T}, s::Integer) where {T<:Integer} = shortest_paths(g, s, BFS())
shortest_paths(g::AbstractGraph{T}, ss::AbstractVector) where {T<:Integer} = shortest_paths(g, ss, BFS())

# Full-formed methods.
"""
    paths(state[, v])
    paths(state[, vs])
    paths(state[, v, d]))

Given the output of a [`shortest_paths`](@ref) calculation of type
[`ShortestPathResult`](@ref), return a vector (indexed by vertex)
of the paths between the source vertex used to compute the shortest path
and a single destination vertex `v`, a vector of destination vertices `vs`,
or the entire graph.

For multiple destination vertices, each path is represented by a vector of
vertices on the path between the source and the destination. Nonexistent
paths will be indicated by an empty vector. For single destinations, the
path is represented by a single vector of vertices, and will be length 0
if the path does not exist.

For [`ShortestPathAlgorithm`](@ref)s that compute all shortest paths for all
pairs of vertices, `paths(state)` will return a vector (indexed by source
vertex) of vectors (indexed by destination vertex) of paths. `paths(state, v)`
will return a vector (indexed by destination vertex) of paths from source `v`
to all other vertices. In addition, `paths(state, v, d)` will return a vector
representing the path from vertex `v` to vertex `d`.
"""
function paths(state::ShortestPathResult, vs::AbstractVector{<:Integer})
    parents = state.parents
    T = eltype(parents)

    num_vs = length(vs)
    all_paths = Vector{Vector{T}}(undef, num_vs)
    for i = 1:num_vs
        all_paths[i] = Vector{T}()
        index = T(vs[i])
        if parents[index] != 0 || parents[index] == index
            while parents[index] != 0
                pushfirst!(all_paths[i], index)
                index = parents[index]
            end
            pushfirst!(all_paths[i], index)
        end
    end
    return all_paths
end

paths(state::ShortestPathResult, v::Integer) = paths(state, [v])[1]
paths(state::ShortestPathResult) = paths(state, [1:length(state.parents);])

"""
    dists(state[, v])

Given the output of a [`shortest_paths`](@ref) calculation of type
[`ShortestPathResult`](@ref), return a vector (indexed by vertex)
of the distances between the source vertex used to compute the
shortest path and a single destination vertex `v` or the entire graph.

For [`ShortestPathAlgorithm`](@ref)s that compute all-pairs shortest
paths, `dists(state)` will return a matrix (indexed by source and destination
vertices) of distances.
"""
dists(state::ShortestPathResult, v::Integer) = state.dists[v]
dists(state::ShortestPathResult) = state.dists

"""
    has_negative_weight_cycle(g[, distmx=weights(g), alg=BellmanFord()])

Given a graph `g`, an optional distance matrix `distmx`, and an optional
algorithm `alg` (one of [`BellmanFord`](@ref) or [`SPFA`](@ref)), return
`true` if any cycle detected in the graph has a negative weight.
# Examples

```jldoctest
julia> g = complete_graph(3);

julia> d = [1 -3 1; -3 1 1; 1 1 1];

julia> has_negative_weight_cycle(g, d)
true

julia> g = complete_graph(4);

julia> d = [1 1 -1 1; 1 1 -1 1; 1 1 1 1; 1 1 1 1];

julia> has_negative_weight_cycle(g, d, SPFA())
false
```
"""
has_negative_weight_cycle(g::AbstractGraph, distmx::AbstractMatrix=weights(g)) = has_negative_weight_cycle(g, distmx, BellmanFord())
has_negative_weight_cycle(g::AbstractSimpleGraph) = false

export ShortestPathAlgorithm
export paths, dists, shortest_paths, has_negative_weight_cycle
export Dijkstra, AStar, BellmanFord, FloydWarshall, DEsopoPape, Johnson, SPFA, BFS
export NegativeCycleError

end  # module


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Datamine

Rotation convention of the Datamine software (fixed to ZXZ axes).
"""
struct Datamine <: RotationConvention end
axesseq(::Type{Datamine}) = :ZXZ
orientation(::Type{Datamine}) = (:CW,:CW,:CW)
angleunits(::Type{Datamine}) = :DEG
mainaxis(::Type{Datamine}) = :X
isextrinsic(::Type{Datamine}) = false


text:
@testset "Spin spherical harmonic mode indices (s=$s)" for s in -2:2
    lmax = 10
    N = lmax + 1
    M = 2 * N - 1

    mode_seen = zeros(Int, N, M)
    for l in abs(s):lmax, m in (-l):l
        mode_seen[spinsph_mode(s, l, m)] += 1
    end
    @test mode_seen == spinsphones(Bool, N, M, s)

    # Test some concrete modes
    @test spinsph_mode(0, 0, 0) == CartesianIndex(1, 1)
    @test spinsph_mode(0, 1, -1) == CartesianIndex(1, 2)
    @test spinsph_mode(0, 1, 1) == CartesianIndex(1, 3)

    @test spinsph_mode(1, 1, 0) == CartesianIndex(1, 1)
    @test spinsph_mode(1, 1, -1) == CartesianIndex(1, 2)
    @test spinsph_mode(1, 1, 1) == CartesianIndex(1, 3)

    @test spinsph_mode(2, 2, 0) == CartesianIndex(1, 1)
    @test spinsph_mode(2, 2, -1) == CartesianIndex(1, 2)
    @test spinsph_mode(2, 2, 1) == CartesianIndex(1, 3)
end

@testset "Spin spherical harmonics: simple modes ($T)" for T in
                                                           [Complex{Float64}]
    lmax = 100

    N = lmax + 1
    Θ, Φ = sph_points(N)
    @test length(Θ) == N
    M = length(Φ)

    # Test s=1, l=1, m=0
    s, l, m = 1, 1, 0
    F = T[Y_s1_1_0(θ, ϕ) for θ in Θ, ϕ in Φ]
    @test size(F) == (N, M)
    C = spinsph_transform(F, s)
    @test C ≈ unit(spinsph_mode(s, l, m), size(C))
    F′ = spinsph_evaluate(C, s)
    @test F′ ≈ F
end

@testset "Spin spherical harmonics: simple modes (s=$s, $T)" for s in -2:2,
                                                                 T in [Float64,
                                                                  Complex{Float64}]

    T <: Real && s ≠ 0 && continue

    lmax = 10

    N = lmax + 1
    Θ, Φ = sph_points(N)
    @test length(Θ) == N
    M = length(Φ)

    lmax_test = 4
    for l in abs(s):lmax_test, m in (-l):l
        F = T[sYlm(T, s, l, m, θ, ϕ) for θ in Θ, ϕ in Φ]
        C = spinsph_transform(F, s)
        @test C ≈ unit(spinsph_mode(s, l, m), size(C))
        F′ = spinsph_evaluate(C, s)
        @test F′ ≈ F
    end
end

@testset "Spin spherical harmonics: linearity (s=$s, $T)" for s in -2:2,
                                                              T in [Float64,
                                                               Complex{Float64}]

    T <: Real && s ≠ 0 && continue

    lmax = 100

    N = lmax + 1
    Θ, Φ = sph_points(N)
    @test length(Θ) == N
    M = length(Φ)

    F = zeros(T, N, M)
    C = spinsph_transform(F, s)
    @test iszero(C)

    F1 = randn(T, N, M)
    F2 = randn(T, N, M)
    C1 = spinsph_transform(F1, s)
    C2 = spinsph_transform(F2, s)
    F = F1 + F2
    C = spinsph_transform(F, s)
    @test C ≈ C1 + C2

    F = randn(T, N, M)
    # Why doesn't this work with a complex scale factor? Is the
    # transform not linear?
    # α = randn(T)
    α = randn()
    C = spinsph_transform(F, s)
    Cα = spinsph_transform(α * F, s)
    @test Cα ≈ α * C
end

@testset "Spin spherical harmonics: duality (s=$s, $T)" for s in -2:2,
                                                            T in [Float64,
                                                             Complex{Float64}]

    T <: Real && s ≠ 0 && continue

    lmax = 100

    N = lmax + 1
    Θ, Φ = sph_points(N)
    @test length(Θ) == N
    M = length(Φ)

    F = randn(T, N, M)
    C = spinsph_transform(F, s)
    @test size(C) == size(F)
    @test eltype(C) == eltype(F)
    F′ = spinsph_evaluate(C, s)
    @test F′ ≈ F

    C = randn(T, N, M)
    F = spinsph_evaluate(C, s)
    @test size(F) == size(C)
    @test eltype(F) == eltype(C)
    C′ = spinsph_transform(F, s)
    @test C′ ≈ C
end

@testset "Spin spherical harmonics: orthonormality (s=$s, $T)" for s in -2:2,
                                                                   T in
                                                                   [Float64,
                                                                    Complex{Float64}]

    T <: Real && s ≠ 0 && continue

    lmax = 100
    atol = 4 / lmax^2

    N = lmax + 1
    Θ, Φ = sph_points(N)
    @test length(Θ) == N
    M = length(Φ)

    lmax_test = 4
    for l in abs(s):lmax_test, m in (-l):l
        C = zeros(T, N, M)
        C[spinsph_mode(s, l, m)] = 1
        F = spinsph_evaluate(C, s)
        for l′ in abs(s):lmax_test, m′ in (-l′):l′
            C′ = zeros(T, N, M)
            C′[spinsph_mode(s, l′, m′)] = 1
            F′ = spinsph_evaluate(C′, s)
            int = integrate(F, F′)
            δ = l == l′ && m == m′
            @test isapprox(int, δ; atol=atol)
        end
    end
end

@testset "Spin spherical harmonics: eth and ethbar (s=$s, $T)" for s in -2:2,
                                                                   T in
                                                                   [Complex{Float64}]

    RT = typeof(Real(one(T)))
    lmax = 10

    N = lmax + 1
    Θ, Φ = sph_points(N)
    @test length(Θ) == N
    M = length(Φ)

    lmax_test = 4

    for l in max(abs(s), abs(s + 1)):lmax_test, m in (-l):l
        F = T[ðsYlm(s, l, m, θ, ϕ) for θ in Θ, ϕ in Φ]
        C = spinsph_transform(F, s + 1)
        α = ifelse(m ≥ -s, 1, -1) * sqrt((l - s) * (l + s + 1))
        @test C ≈ unit(α, spinsph_mode(s + 1, l, m), size(C))
        F′ = spinsph_evaluate(C, s + 1)
        @test F′ ≈ F
    end

    for l in max(abs(s), abs(s - 1)):lmax_test, m in (-l):l
        F = T[ð̄sYlm(s, l, m, θ, ϕ) for θ in Θ, ϕ in Φ]
        C = spinsph_transform(F, s - 1)
        α = -ifelse(m > -s, 1, -1) * sqrt((l + s) * (l - s + 1))
        @test C ≈ unit(α, spinsph_mode(s - 1, l, m), size(C))
        F′ = spinsph_evaluate(C, s - 1)
        @test F′ ≈ F
    end

    for l in abs(s):lmax_test, m in (-l):l
        C = zeros(T, N, M)
        C[spinsph_mode(s, l, m)] = ifelse(m ≥ -s, 1, -1)

        ðC = spinsph_eth(C, s)
        ðF = spinsph_evaluate(ðC, s + 1)
        @test ðF .+ 1 ≈ T[ðsYlm(s, l, m, θ, ϕ) for θ in Θ, ϕ in Φ] .+ 1
    end

    for l in abs(s):lmax_test, m in (-l):l
        C = zeros(T, N, M)
        C[spinsph_mode(s, l, m)] = ifelse(m > -s, 1, -1)

        ð̄C = spinsph_ethbar(C, s)
        ð̄F = spinsph_evaluate(ð̄C, s - 1)
        @test ð̄F .+ 1 ≈ T[ð̄sYlm(s, l, m, θ, ϕ) for θ in Θ, ϕ in Φ] .+ 1
    end
end

@testset "Spin spherical harmonics: Laplacian ($T)" for T in [Complex{Float64}]
    lmax = 100

    N = lmax + 1
    Θ, Φ = sph_points(N)
    @test length(Θ) == N
    M = length(Φ)

    # Calculate Laplacian directly
    RT = typeof(Real(one(T)))
    F = randn(RT, N, M)
    C = sph_transform(F)
    ΔC = sph_laplace(C)
    ΔF = sph_evaluate(ΔC)

    # Calculate Laplacian via eth and ethbar
    F⁰ = T.(F)
    C⁰ = spinsph_transform(F⁰, 0)
    @test eltype(C⁰) === T
    ðC¹ = spinsph_eth(C⁰, 0)
    @test eltype(ðC¹) === T
    ð̄ðC⁰ = spinsph_ethbar(ðC¹, 1)
    @test eltype(ð̄ðC⁰) === T
    ð̄ðF⁰ = spinsph_evaluate(ð̄ðC⁰, 0)
    @test eltype(ð̄ðF⁰) === T
    @test ð̄ðF⁰ ≈ ΔF

    # Calculate Laplacian via ethbar and eth
    ð̄C⁻¹ = spinsph_ethbar(C⁰, 0)
    @test eltype(ð̄C⁻¹) === T
    ðð̄C⁰ = spinsph_eth(ð̄C⁻¹, -1)
    @test eltype(ðð̄C⁰) === T
    ðð̄F⁰ = spinsph_evaluate(ðð̄C⁰, 0)
    @test eltype(ð̄ðF⁰) === T
    @test ðð̄F⁰ ≈ ΔF
end

@testset "Spin spherical harmonics: grad ($T)" for T in [Float64]
    lmax = 10

    N = lmax + 1
    Θ, Φ = sph_points(N)
    @test length(Θ) == N
    M = length(Φ)

    lmax_test = 4

    for l in 0:lmax_test, m in (-l):l
        F = Complex{T}[sYlm(0, l, m, θ, ϕ) for θ in Θ, ϕ in Φ]
        ðF₀ = Complex{T}[ðsYlm(0, l, m, θ, ϕ) for θ in Θ, ϕ in Φ]
        ð̄F₀ = Complex{T}[ð̄sYlm(0, l, m, θ, ϕ) for θ in Θ, ϕ in Φ]

        C = spinsph_transform(F, 0)
        @test C ≈ unit(1, spinsph_mode(0, l, m), size(C))

        ðC = spinsph_eth(C, 0)
        ð̄C = spinsph_ethbar(C, 0)

        ðF = spinsph_evaluate(ðC, 1)
        ð̄F = spinsph_evaluate(ð̄C, -1)
        # Different sYlm modes need different signs:
        α = ifelse(m ≥ 0, 1, -1)
        ᾱ = ifelse(m > 0, 1, -1)

        if l == 0
            @test norm(ðF) ≤ sqrt(eps(T))
            @test norm(ð̄F) ≤ sqrt(eps(T))
        else
            @test α * ðF ≈ ðF₀
            @test ᾱ * ð̄F ≈ ð̄F₀
        end
    end

    for l in 0:lmax_test, m in (-l):l
        F = T[Ylm(l, m, θ, ϕ) for θ in Θ, ϕ in Φ]
        ðF₀ = SVector{2,T}[c2a(ðYlm(l, m, θ, ϕ)) for θ in Θ, ϕ in Φ]
        ð̄F₀ = SVector{2,T}[c2a(ð̄Ylm(l, m, θ, ϕ)) for θ in Θ, ϕ in Φ]
        ðC₀ = spinsph_transform(ðF₀, 1)
        ð̄C₀ = spinsph_transform(ð̄F₀, -1)
        ðC₀::Array{SVector{2,T},2}
        ð̄C₀::Array{SVector{2,T},2}

        C = spinsph_transform(F, 0)
        C::Array{T,2}
        for l′ in 0:lmax, m′ in (-l′):l′
            # Allow ±m modes to mix
            if (l′, abs(m′)) ≠ (l, abs(m))
                @test norm(C[spinsph_mode(0, l′, m′)]) ≤ sqrt(eps(T))
            end
        end

        ðC = spinsph_eth(C, 0)
        ð̄C = spinsph_ethbar(C, 0)
        ðC::Array{SVector{2,T},2}
        ð̄C::Array{SVector{2,T},2}

        if l == 0
            @test norm(ðC) ≤ sqrt(eps(T))
            @test norm(ð̄C) ≤ sqrt(eps(T))
        else
            @test ðC ≈ ðC₀
            @test ð̄C ≈ ð̄C₀
        end

        ðF = spinsph_evaluate(ðC, 1)
        ð̄F = spinsph_evaluate(ð̄C, -1)
        ðF::Array{SVector{2,T},2}
        ð̄F::Array{SVector{2,T},2}
        if l == 0
            @test norm(ðF) ≤ sqrt(eps(T))
            @test norm(ð̄F) ≤ sqrt(eps(T))
        else
            @test ðF ≈ ðF₀
            @test ð̄F ≈ ð̄F₀
        end
    end
end

@testset "Spin spherical harmonics: div ($T)" for T in [Float64]
    lmax = 10

    N = lmax + 1
    Θ, Φ = sph_points(N)
    @test length(Θ) == N
    M = length(Φ)

    lmax_test = 4

    for l in 0:lmax_test, m in (-l):l
        ðF = Complex{T}[ðsYlm(0, l, m, θ, ϕ) for θ in Θ, ϕ in Φ]
        ð̄F = Complex{T}[ð̄sYlm(0, l, m, θ, ϕ) for θ in Θ, ϕ in Φ]
        ΔF₀ = Complex{T}[-l * (l + 1) * sYlm(0, l, m, θ, ϕ) for θ in Θ, ϕ in Φ]

        ðC = spinsph_transform(ðF, 1)
        ð̄C = spinsph_transform(ð̄F, -1)
        ð̄ðC = spinsph_ethbar(ðC, 1)
        ðð̄C = spinsph_eth(ð̄C, -1)

        ð̄ðF = spinsph_evaluate(ð̄ðC, 0)
        ðð̄F = spinsph_evaluate(ðð̄C, 0)
        # Different sYlm modes need different signs:
        α = ifelse(m ≥ 0, 1, -1)
        ᾱ = ifelse(m > 0, 1, -1)

        if l == 0
            @test norm(ð̄ðF) ≤ sqrt(eps(T))
            @test norm(ðð̄F) ≤ sqrt(eps(T))
        else
            @test α * ð̄ðF ≈ ΔF₀
            @test ᾱ * ðð̄F ≈ ΔF₀
        end
    end

    for l in 0:lmax_test, m in (-l):l
        ðF = SVector{2,T}[c2a(ðYlm(l, m, θ, ϕ)) for θ in Θ, ϕ in Φ]
        ð̄F = SVector{2,T}[c2a(ð̄Ylm(l, m, θ, ϕ)) for θ in Θ, ϕ in Φ]
        ΔF₀ = T[-l * (l + 1) * Ylm(l, m, θ, ϕ) for θ in Θ, ϕ in Φ]
        ΔC₀ = spinsph_transform(ΔF₀, 0)
        ΔC₀::Array{T,2}

        ðC = spinsph_transform(ðF, 1)
        ð̄C = spinsph_transform(ð̄F, -1)
        ðC::Array{SVector{2,T},2}
        ð̄C::Array{SVector{2,T},2}
        ð̄ðC = spinsph_ethbar(ðC, 1)
        ðð̄C = spinsph_eth(ð̄C, -1)
        ð̄ðC::Array{T,2}
        ðð̄C::Array{T,2}

        if l == 0
            @test norm(ð̄ðC) ≤ sqrt(eps(T))
            @test norm(ðð̄C) ≤ sqrt(eps(T))
        else
            @test ð̄ðC ≈ ΔC₀
            @test ðð̄C ≈ ΔC₀
        end

        ð̄ðF = spinsph_evaluate(ð̄ðC, 0)
        ðð̄F = spinsph_evaluate(ðð̄C, 0)
        ð̄ðF::Array{T,2}
        ðð̄F::Array{T,2}

        if l == 0
            @test norm(ð̄ðF) ≤ sqrt(eps(T))
            @test norm(ðð̄F) ≤ sqrt(eps(T))
        else
            @test ð̄ðF ≈ ΔF₀
            @test ðð̄F ≈ ΔF₀
        end
    end
end

@testset "Spin spherical harmonics: Laplacian ($T)" for T in [Float64]
    lmax = 10

    N = lmax + 1
    Θ, Φ = sph_points(N)
    @test length(Θ) == N
    M = length(Φ)

    for iter in 1:10
        F = randn(T, N, M)

        C = spinsph_transform(F, 0)
        ðC = spinsph_eth(C, 0)
        ð̄C = spinsph_ethbar(C, 0)
        ð̄ðC = spinsph_ethbar(ðC, 1)
        ðð̄C = spinsph_eth(ð̄C, -1)

        for l in 0:lmax, m in (-l):l
            if l == 0
                @test norm(ð̄ðC[spinsph_mode(0, l, m)]) ≤ sqrt(eps(T))
                @test norm(ðð̄C[spinsph_mode(0, l, m)]) ≤ sqrt(eps(T))
            else
                @test ð̄ðC[spinsph_mode(0, l, m)] ≈
                      -l * (l + 1) * C[spinsph_mode(0, l, m)]
                @test ðð̄C[spinsph_mode(0, l, m)] ≈
                      -l * (l + 1) * C[spinsph_mode(0, l, m)]
            end
        end

        ð̄ðF = spinsph_evaluate(ð̄ðC, 0)
        ðð̄F = spinsph_evaluate(ðð̄C, 0)

        C₀ = sph_transform(F)
        ΔC₀ = sph_laplace(C₀)
        ΔF₀ = sph_evaluate(ΔC₀)

        @test ð̄ðF ≈ ΔF₀
        @test ðð̄F ≈ ΔF₀
    end
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule CUDNN_CUDA112_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("CUDNN_CUDA112")
JLLWrappers.@generate_main_file("CUDNN_CUDA112", UUID("d433cf69-735b-5d73-bb80-2d42f39ec054"))
end  # module CUDNN_CUDA112_jll


text:
using ConvOptDL
using Test
using StatsBase

@testset "data_loader.jl" begin
    dloader = FewShotDataLoader("./test_data.jls")
    meta_sample = sample(dloader, 2)
    @test size(meta_sample) == 2
end


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

"""
    ConstValueShape{T} <: AbstractValueShape

A `ConstValueShape` describes the shape of constant values of type `T`.

Constructor:

    ConstValueShape(value)

`value` may be of arbitrary type, e.g. a constant scalar value or array:

    ConstValueShape(4.2)
    ConstValueShape([11 21; 12 22])

Shapes of constant values have zero degrees of freedom (see
[`totalndof`](@ref)).

See also the documentation of [`AbstractValueShape`](@ref).
"""
struct ConstValueShape{T} <: AbstractValueShape
    value::T
end

export ConstValueShape

@inline Base.size(shape::ConstValueShape) = size(shape.value)
@inline Base.length(shape::ConstValueShape) = length(shape.value)

@inline Base.eltype(shape::ConstValueShape) = Int

@inline totalndof(::ConstValueShape) = 0

# ToDo/Decision: Return a copy instead?
(shape::ConstValueShape)(::UndefInitializer) = shape.value

@static if VERSION < v"1.3"
    # Workaround for Julia issue #14919
    @inline (shape::ConstValueShape)(data::AbstractVector{<:Real}) =
        _apply_shape_to_data(shape, data)
end

const ConstAccessor = ValueAccessor{ConstValueShape{T}} where {T}

@inline Base.getindex(data::AbstractVector{<:Real}, va::ConstAccessor) = view(data, va)

@inline Base.view(::AbstractVector, va::ConstAccessor) = va.shape.value

@inline _bcasted_getindex(data::AbstractVectorOfSimilarVectors{<:Real}, va::ConstAccessor) =
    Fill(va.shape.value, size(data,1))

Base.copy(instance::VSBroadcasted2{typeof(getindex),AbstractVectorOfSimilarVectors{<:Real},Ref{<:ConstAccessor}}) =
    _bcasted_getindex(instance.args[1], instance.args[2][])


text:
importall Bukdu

type WelcomeController <: ApplicationController
end

index(::WelcomeController) = "hello world"

Router() do
    get("/", WelcomeController, index)
end

loaded = []
Endpoint() do
    push!(loaded, 1)
    plug(Plug.Static, at= "/", from= "../examples/public")
    plug(Plug.Logger, level=:error)
    plug(Router)
end

using Base.Test
conn = (Endpoint)("/index.html")
@test 200 == conn.status
@test "text/html" == conn.resp_header["Content-Type"]
@test startswith(String(conn.resp_body), "<!DOCTYPE html>")

conn = (Endpoint)("/")
@test 200 == conn.status
@test startswith(String(conn.resp_body), "<!DOCTYPE html>")

conn = (Endpoint)("/js/vue.min.js")
@test 200 == conn.status
@test "application/javascript" == conn.resp_header["Content-Type"]
@test 76807 == sizeof(conn.resp_body)

@test_throws NoRouteError (Endpoint)("/js/not_found")

@test [1] == loaded

reload(Endpoint)

@test [1, 1] == loaded