#!/usr/bin/env julia

using Statistics

function t_val(group_fast, group_slow)
    if size(group_fast, 1) == 0 || size(group_slow, 1) == 0
        return 0
    end
    mean_slow = mean(group_slow)
    mean_fast = mean(group_fast)
    s_slow = var(group_slow, mean=mean_slow)
    s_fast = var(group_fast, mean=mean_fast)
    return (mean_slow - mean_fast) / sqrt(s_slow / size(group_slow, 1) + s_fast / size(group_fast, 1))
end

function filter_peaks(data::Matrix{T})::Matrix{T} where T
    m = mean(data, dims=2)
    blocks = data .> m
    blocks = hcat(.!blocks[:, 1], blocks)
    is_edge = blocks[:, 1:end-1] .!= blocks[:, 2:end]
    no_peaks = maximum([count(is_edge[row, :]) - 1 for row=1:size(data, 1)])
    peaks = zeros(T, size(data, 1), no_peaks)
    @Threads.threads for row_no=1:size(data, 1)
        row = data[row_no, :]
        separators = findall(is_edge[row_no, :])
        for block_no=1:size(separators, 1)-1
            lsep, rsep = separators[block_no:block_no+1]
            bm = mean(row[lsep:rsep])
            if bm > m[row_no]
                peaks[row_no, block_no] = maximum(row[lsep:rsep])
            else
                peaks[row_no, block_no] = minimum(row[lsep:rsep])
            end
        end
    end
    return peaks
end

function load_file(filename, type::Type{T})::Matrix{T} where T
    open(filename, "r") do file
        field_size = read(file, UInt8)
        if field_size != sizeof(T)
            throw(ArgumentError(format("Expected type of size {:d} but passed type '{}' has size {:d}", field_size, T, sizeof(T))))
        end
        rows = read(file, UInt64)
        cols = read(file, UInt64)
        data = Matrix{T}(undef, rows, cols)
        read!(file, data)
        return data
    end
end