text:
AnyIntType = Union{Int8, Int32, Int64, Int128, Int}
AnyDoubleType = Union{Int8, Int32, Int64, Int128, Int, Float16, Float32, Float64}

@enum AggOperations begin

    MEAN = 0
    MEDIAN = 1
    EACH = 2

end

@enum FilterTypes begin

    BUTTERWORTH = 0
    CHEBYSHEV_TYPE_1 = 1
    BESSEL = 2

end

@enum WindowFunctions begin

    NO_WINDOW = 0
    HANNING = 1
    HAMMING = 2
    BLACKMAN_HARRIS = 3

end

@enum DetrendOperations begin

    # NONE = 0 # will lead to invalid enum redefinitions since enums entries are global constants in fact
    # its a temp hack, lets wait for fixed enums in julia
    CONSTANT = 1
    LINEAR = 2

end

function perform_lowpass(data, sampling_rate::AnyIntType, cutoff::AnyDoubleType, order::AnyIntType,
    filter_type::AnyIntType, ripple::AnyDoubleType)
    ec = STATUS_OK
    # due to this bug https://github.com/JuliaLang/julia/issues/29602 libname should be hardcoded
    if Sys.iswindows()
        ec = ccall((:perform_lowpass, "DataHandler.dll"), Cint, (Ptr{Float64}, Cint, Cint, Float64, Cint, Cint, Float64),
            data, length(data), Int32(sampling_rate), Float64(cutoff), Int32(order), Int32(filter_type), Float64(ripple))
    elseif Sys.isapple()
        ec = ccall((:perform_lowpass, "libDataHandler.dylib"), Cint, (Ptr{Float64}, Cint, Cint, Float64, Cint, Cint, Float64),
            data, length(data), Int32(sampling_rate), Float64(cutoff), Int32(order), Int32(filter_type), Float64(ripple))
    else
        ec = ccall((:perform_lowpass, "libDataHandler.so"), Cint, (Ptr{Float64}, Cint, Cint, Float64, Cint, Cint, Float64),
            data, length(data), Int32(sampling_rate), Float64(cutoff), Int32(order), Int32(filter_type), Float64(ripple))
    end
    if ec != Integer(STATUS_OK)
        throw(BrainFlowError(string("Error in perform_lowpass ", ec), ec))
    end
    return
end

function perform_highpass(data, sampling_rate::AnyIntType, cutoff::AnyDoubleType, order::AnyIntType,
    filter_type::AnyIntType, ripple::AnyDoubleType)
    ec = STATUS_OK
    # due to this bug https://github.com/JuliaLang/julia/issues/29602 libname should be hardcoded
    if Sys.iswindows()
        ec = ccall((:perform_highpass, "DataHandler.dll"), Cint, (Ptr{Float64}, Cint, Cint, Float64, Cint, Cint, Float64),
            data, length(data), Int32(sampling_rate), Float64(cutoff), Int32(order), Int32(filter_type), Float64(ripple))
    elseif Sys.isapple()
        ec = ccall((:perform_highpass, "libDataHandler.dylib"), Cint, (Ptr{Float64}, Cint, Cint, Float64, Cint, Cint, Float64),
            data, length(data), Int32(sampling_rate), Float64(cutoff), Int32(order), Int32(filter_type), Float64(ripple))
    else
        ec = ccall((:perform_highpass, "libDataHandler.so"), Cint, (Ptr{Float64}, Cint, Cint, Float64, Cint, Cint, Float64),
            data, length(data), Int32(sampling_rate), Float64(cutoff), Int32(order), Int32(filter_type), Float64(ripple))
    end
    if ec != Integer(STATUS_OK)
        throw(BrainFlowError(string("Error in perform_highpass ", ec), ec))
    end
    return
end

function perform_bandpass(data, sampling_rate::AnyIntType, center_freq::AnyDoubleType,
    band_width::AnyDoubleType, order::AnyIntType, filter_type::AnyIntType, ripple::AnyDoubleType)
    ec = STATUS_OK
    # due to this bug https://github.com/JuliaLang/julia/issues/29602 libname should be hardcoded
    if Sys.iswindows()
        ec = ccall((:perform_bandpass, "DataHandler.dll"), Cint, (Ptr{Float64}, Cint, Cint, Float64, Float64, Cint, Cint, Float64),
            data, length(data), Int32(sampling_rate), Float64(center_freq), Float64(band_width), Int32(order), Int32(filter_type), Float64(ripple))
    elseif Sys.isapple()
        ec = ccall((:perform_bandpass, "libDataHandler.dylib"), Cint, (Ptr{Float64}, Cint, Cint, Float64, Float64, Cint, Cint, Float64),
            data, length(data), Int32(sampling_rate), Float64(center_freq), Float64(band_width), Int32(order), Int32(filter_type), Float64(ripple))
    else
        ec = ccall((:perform_bandpass, "libDataHandler.so"), Cint, (Ptr{Float64}, Cint, Cint, Float64, Float64, Cint, Cint, Float64),
            data, length(data), Int32(sampling_rate), Float64(center_freq), Float64(band_width), Int32(order), Int32(filter_type), Float64(ripple))
    end
    if ec != Integer(STATUS_OK)
        throw(BrainFlowError(string("Error in perform_bandpass ", ec), ec))
    end
    return
end

function perform_bandstop(data, sampling_rate::AnyIntType, center_freq::AnyDoubleType,
    band_width::AnyDoubleType, order::AnyIntType, filter_type::AnyIntType, ripple::AnyDoubleType)
    ec = STATUS_OK
    # due to this bug https://github.com/JuliaLang/julia/issues/29602 libname should be hardcoded
    if Sys.iswindows()
        ec = ccall((:perform_bandstop, "DataHandler.dll"), Cint, (Ptr{Float64}, Cint, Cint, Float64, Float64, Cint, Cint, Float64),
            data, length(data), Int32(sampling_rate), Float64(center_freq), Float64(band_width), Int32(order), Int32(filter_type), Float64(ripple))
    elseif Sys.isapple()
        ec = ccall((:perform_bandstop, "libDataHandler.dylib"), Cint, (Ptr{Float64}, Cint, Cint, Float64, Float64, Cint, Cint, Float64),
            data, length(data), Int32(sampling_rate), Float64(center_freq), Float64(band_width), Int32(order), Int32(filter_type), Float64(ripple))
    else
        ec = ccall((:perform_bandstop, "libDataHandler.so"), Cint, (Ptr{Float64}, Cint, Cint, Float64, Float64, Cint, Cint, Float64),
            data, length(data), Int32(sampling_rate), Float64(center_freq), Float64(band_width), Int32(order), Int32(filter_type), Float64(ripple))
    end
    if ec != Integer(STATUS_OK)
        throw(BrainFlowError(string("Error in perform_bandstop ", ec), ec))
    end
    return
end

function perform_rolling_filter(data, period::AnyIntType, operation::AnyIntType)
    ec = STATUS_OK
    # due to this bug https://github.com/JuliaLang/julia/issues/29602 libname should be hardcoded
    if Sys.iswindows()
        ec = ccall((:perform_rolling_filter, "DataHandler.dll"), Cint, (Ptr{Float64}, Cint, Cint, Cint),
            data, length(data), Int32(period), Int32(operation))
    elseif Sys.isapple()
        ec = ccall((:perform_rolling_filter, "libDataHandler.dylib"), Cint, (Ptr{Float64}, Cint, Cint, Cint),
            data, length(data), Int32(period), Int32(operation))
    else
        ec = ccall((:perform_rolling_filter, "libDataHandler.so"), Cint, (Ptr{Float64}, Cint, Cint, Cint),
            data, length(data), Int32(period), Int32(operation))
    end
    if ec != Integer(STATUS_OK)
        throw(BrainFlowError(string("Error in perform_rolling_filter ", ec), ec))
    end
    return
end

function detrend(data, operation::AnyIntType)
    ec = STATUS_OK
    # due to this bug https://github.com/JuliaLang/julia/issues/29602 libname should be hardcoded
    if Sys.iswindows()
        ec = ccall((:detrend, "DataHandler.dll"), Cint, (Ptr{Float64}, Cint, Cint),
            data, length(data), Int32(operation))
    elseif Sys.isapple()
        ec = ccall((:detrend, "libDataHandler.dylib"), Cint, (Ptr{Float64}, Cint, Cint),
            data, length(data), Int32(operation))
    else
        ec = ccall((:detrend, "libDataHandler.so"), Cint, (Ptr{Float64}, Cint, Cint),
            data, length(data), Int32(operation))
    end
    if ec != Integer(STATUS_OK)
        throw(BrainFlowError(string("Error in detrend ", ec), ec))
    end
    return
end

function perform_wavelet_denoising(data, wavelet::String, decomposition_level::AnyIntType)
    ec = STATUS_OK
    # due to this bug https://github.com/JuliaLang/julia/issues/29602 libname should be hardcoded
    if Sys.iswindows()
        ec = ccall((:perform_wavelet_denoising, "DataHandler.dll"), Cint, (Ptr{Float64}, Cint, Ptr{UInt8}, Cint),
            data, length(data), wavelet, Int32(decomposition_level))
    elseif Sys.isapple()
        ec = ccall((:perform_wavelet_denoising, "libDataHandler.dylib"), Cint, (Ptr{Float64}, Cint, Ptr{UInt8}, Cint),
            data, length(data), wavelet, Int32(decomposition_level))
    else
        ec = ccall((:perform_wavelet_denoising, "libDataHandler.so"), Cint, (Ptr{Float64}, Cint, Ptr{UInt8}, Cint),
            data, length(data), wavelet, Int32(decomposition_level))
    end
    if ec != Integer(STATUS_OK)
        throw(BrainFlowError(string("Error in perform_wavelet_denoising ", ec), ec))
    end
    return
end

function perform_downsampling(data, period::AnyIntType, operation::AnyIntType)
    len = Integer(floor(length(data) / period))
    downsampled_data = Vector{Float64}(undef, len)
    ec = STATUS_OK
    # due to this bug https://github.com/JuliaLang/julia/issues/29602 libname should be hardcoded
    if Sys.iswindows()
        ec = ccall((:perform_downsampling, "DataHandler.dll"), Cint, (Ptr{Float64}, Cint, Cint, Cint, Ptr{Float64}),
            data, length(data), Int32(period), Int32(operation), downsampled_data)
    elseif Sys.isapple()
        ec = ccall((:perform_downsampling, "libDataHandler.dylib"), Cint, (Ptr{Float64}, Cint, Cint, Cint, Ptr{Float64}),
            data, length(data), Int32(period), Int32(operation), downsampled_data)
    else
        ec = ccall((:perform_downsampling, "libDataHandler.so"), Cint, (Ptr{Float64}, Cint, Cint, Cint, Ptr{Float64}),
            data, length(data), Int32(period), Int32(operation), downsampled_data)
    end
    if ec != Integer(STATUS_OK)
        throw(BrainFlowError(string("Error in perform_downsampling ", ec), ec))
    end
    downsampled_data
end

function write_file(data, file_name::String, file_mode::String)
    shape = size(data)
    flatten = reshape(transpose(data), (1, shape[1] * shape[2]))
    flatten = copy(flatten)
    ec = STATUS_OK
    # due to this bug https://github.com/JuliaLang/julia/issues/29602 libname should be hardcoded
    if Sys.iswindows()
        ec = ccall((:write_file, "DataHandler.dll"), Cint, (Ptr{Float64}, Cint, Cint, Ptr{UInt8}, Ptr{UInt8}),
            flatten, shape[1], shape[2], file_name, file_mode)
    elseif Sys.isapple()
        ec = ccall((:write_file, "libDataHandler.dylib"), Cint, (Ptr{Float64}, Cint, Cint, Ptr{UInt8}, Ptr{UInt8}),
            flatten, shape[1], shape[2], file_name, file_mode)
    else
        ec = ccall((:write_file, "libDataHandler.so"), Cint, (Ptr{Float64}, Cint, Cint, Ptr{UInt8}, Ptr{UInt8}),
            flatten, shape[1], shape[2], file_name, file_mode)
    end
    if ec != Integer(STATUS_OK)
        throw(BrainFlowError(string("Error in write_file ", ec), ec))
    end
    return
end

function get_nearest_power_of_two(value::AnyIntType)
    power_of_two = Vector{Cint}(undef, 1)
    ec = STATUS_OK
    # due to this bug https://github.com/JuliaLang/julia/issues/29602 libname should be hardcoded
    if Sys.iswindows()
        ec = ccall((:get_nearest_power_of_two, "DataHandler.dll"), Cint, (Cint, Ptr{Cint}),
                Int32(value), power_of_two)
    elseif Sys.isapple()
        ec = ccall((:get_nearest_power_of_two, "libDataHandler.dylib"), Cint, (Cint, Ptr{Cint}),
                Int32(value), power_of_two)
    else
        ec = ccall((:get_nearest_power_of_two, "libDataHandler.so"), Cint, (Cint, Ptr{Cint}),
                Int32(value), power_of_two)
    end
    if ec != Integer(STATUS_OK)
        throw(BrainFlowError(string("Error in get_nearest_power_of_two ", ec), ec))
    end

    power_of_two[1]
end

function read_file(file_name::String)
    num_elements = Vector{Cint}(undef, 1)
    ec = STATUS_OK
    # due to this bug https://github.com/JuliaLang/julia/issues/29602 libname should be hardcoded
    if Sys.iswindows()
        ec = ccall((:get_num_elements_in_file, "DataHandler.dll"), Cint, (Ptr{UInt8}, Ptr{Cint}),
                file_name, num_elements)
    elseif Sys.isapple()
        ec = ccall((:get_num_elements_in_file, "libDataHandler.dylib"), Cint, (Ptr{UInt8}, Ptr{Cint}),
                file_name, num_elements)
    else
        ec = ccall((:get_num_elements_in_file, "libDataHandler.so"), Cint, (Ptr{UInt8}, Ptr{Cint}),
                file_name, num_elements)
    end
    if ec != Integer(STATUS_OK)
        throw(BrainFlowError(string("Error in get_num_elements_in_file ", ec), ec))
    end

    data = Vector{Float64}(undef, num_elements[1])
    num_cols = Vector{Cint}(undef, 1)
    num_rows = Vector{Cint}(undef, 1)

    if Sys.iswindows()
        ec = ccall((:read_file, "DataHandler.dll"), Cint, (Ptr{Float64}, Ptr{Cint}, Ptr{Cint}, Ptr{UInt8}, Cint),
                data, num_rows, num_cols, file_name, num_elements[1])
    elseif Sys.isapple()
        ec = ccall((:read_file, "libDataHandler.dylib"), Cint, (Ptr{Float64}, Ptr{Cint}, Ptr{Cint}, Ptr{UInt8}, Cint),
                data, num_rows, num_cols, file_name, num_elements[1])
    else
        ec = ccall((:read_file, "libDataHandler.so"), Cint, (Ptr{Float64}, Ptr{Cint}, Ptr{Cint}, Ptr{UInt8}, Cint),
                data, num_rows, num_cols, file_name, num_elements[1])
    end
    if ec != Integer(STATUS_OK)
        throw(BrainFlowError(string("Error in read_file ", ec), ec))
    end
    value = transpose(reshape(data, (num_cols[1], num_rows[1])))
    value
end

function perform_wavelet_transform(data, wavelet::String, decomposition_level::AnyIntType)
    wavelet_coeffs = Vector{Float64}(undef, length(data) + 2 * (40 + 1))
    lengths = Vector{Cint}(undef, decomposition_level + 1)
    ec = STATUS_OK
    # due to this bug https://github.com/JuliaLang/julia/issues/29602 libname should be hardcoded
    if Sys.iswindows()
        ec = ccall((:perform_wavelet_transform, "DataHandler.dll"), Cint, (Ptr{Float64}, Cint, Ptr{UInt8}, Cint, Ptr{Float64}, Ptr{Cint}),
            data, length(data), wavelet, Int32(decomposition_level), wavelet_coeffs, lengths)
    elseif Sys.isapple()
        ec = ccall((:perform_wavelet_transform, "libDataHandler.dylib"), Cint, (Ptr{Float64}, Cint, Ptr{UInt8}, Cint, Ptr{Float64}, Ptr{Cint}),
            data, length(data), wavelet, Int32(decomposition_level), wavelet_coeffs, lengths)
    else
        ec = ccall((:perform_wavelet_transform, "libDataHandler.so"), Cint, (Ptr{Float64}, Cint, Ptr{UInt8}, Cint, Ptr{Float64}, Ptr{Cint}),
            data, length(data), wavelet, Int32(decomposition_level), wavelet_coeffs, lengths)
    end
    if ec != Integer(STATUS_OK)
        throw(BrainFlowError(string("Error in perform_wavelet_transform ", ec), ec))
    end
    wavelet_coeffs[1:sum(lengths)], lengths
end

function perform_inverse_wavelet_transform(wavelet_output, original_data_len::AnyIntType, wavelet::String, decomposition_level::AnyIntType)
    original_data = Vector{Float64}(undef, original_data_len)
    ec = STATUS_OK
    # due to this bug https://github.com/JuliaLang/julia/issues/29602 libname should be hardcoded
    if Sys.iswindows()
        ec = ccall((:perform_inverse_wavelet_transform, "DataHandler.dll"), Cint, (Ptr{Float64}, Cint, Ptr{UInt8}, Cint, Ptr{Float64}, Ptr{Float64}),
            wavelet_output[1], Int32(original_data_len), wavelet, Int32(decomposition_level), wavelet_output[2], original_data)
    elseif Sys.isapple()
        ec = ccall((:perform_inverse_wavelet_transform, "libDataHandler.dylib"), Cint, (Ptr{Float64}, Cint, Ptr{UInt8}, Cint, Ptr{Float64}, Ptr{Float64}),
            wavelet_output[1], Int32(original_data_len), wavelet, Int32(decomposition_level), wavelet_output[2], original_data)
    else
        ec = ccall((:perform_inverse_wavelet_transform, "libDataHandler.so"), Cint, (Ptr{Float64}, Cint, Ptr{UInt8}, Cint, Ptr{Float64}, Ptr{Float64}),
            wavelet_output[1], Int32(original_data_len), wavelet, Int32(decomposition_level), wavelet_output[2], original_data)
    end
    if ec != Integer(STATUS_OK)
        throw(BrainFlowError(string("Error in perform_inverse_wavelet_transform ", ec), ec))
    end
    original_data
end

function perform_fft(data, window::AnyIntType)

    function is_power_of_two(value)
        (value != 0) && (value & (value - 1) == 0)
    end

    if !is_power_of_two(length(data))
        throw(BrainFlowError(string("Data Len must be power of two ", INVALID_ARGUMENTS_ERROR), INVALID_ARGUMENTS_ERROR))
    end

    temp_re = Vector{Float64}(undef, Integer(length(data) / 2) + 1)
    temp_im = Vector{Float64}(undef, Integer(length(data) / 2) + 1)
    res = Vector{Complex}(undef, Integer(length(data) / 2) + 1)

    ec = STATUS_OK
    # due to this bug https://github.com/JuliaLang/julia/issues/29602 libname should be hardcoded
    if Sys.iswindows()
        ec = ccall((:perform_fft, "DataHandler.dll"), Cint, (Ptr{Float64}, Cint, Cint, Ptr{Float64}, Ptr{Float64}),
            data, length(data), Int32(window), temp_re, temp_im)
    elseif Sys.isapple()
        ec = ccall((:perform_fft, "libDataHandler.dylib"), Cint, (Ptr{Float64}, Cint, Cint, Ptr{Float64}, Ptr{Float64}),
            data, length(data), Int32(window), temp_re, temp_im)
    else
        ec = ccall((:perform_fft, "libDataHandler.so"), Cint, (Ptr{Float64}, Cint, Cint, Ptr{Float64}, Ptr{Float64}),
            data, length(data), Int32(window), temp_re, temp_im)
    end
    if ec != Integer(STATUS_OK)
        throw(BrainFlowError(string("Error in perform_fft ", ec), ec))
    end
    for i in 1:Integer(length(data) / 2) + 1
        res[i] = Complex(temp_re[i], temp_im[i])
    end
    res
end

function perform_ifft(data)

    temp_re = Vector{Float64}(undef, length(data))
    temp_im = Vector{Float64}(undef, length(data))
    res = Vector{Float64}(undef, 2 * (length(data) - 1))

    for i in 1:length(data)
        temp_re[i] = data[i].re
        temp_im[i] = data[i].im
    end
    ec = STATUS_OK
    # due to this bug https://github.com/JuliaLang/julia/issues/29602 libname should be hardcoded
    if Sys.iswindows()
        ec = ccall((:perform_ifft, "DataHandler.dll"), Cint, (Ptr{Float64}, Ptr{Float64}, Cint, Ptr{Float64}),
            temp_re, temp_im, length(res), res)
    elseif Sys.isapple()
        ec = ccall((:perform_ifft, "libDataHandler.dylib"), Cint, (Ptr{Float64}, Ptr{Float64}, Cint, Ptr{Float64}),
            temp_re, temp_im, length(res), res)
    else
        ec = ccall((:perform_ifft, "libDataHandler.so"), Cint, (Ptr{Float64}, Ptr{Float64}, Cint, Ptr{Float64}),
            temp_re, temp_im, length(res), res)
    end
    if ec != Integer(STATUS_OK)
        throw(BrainFlowError(string("Error in perform_ifft ", ec), ec))
    end
    res
end

function get_psd(data, sampling_rate::AnyIntType, window::AnyIntType)

    function is_power_of_two(value)
        (value != 0) && (value & (value - 1) == 0)
    end

    if !is_power_of_two(length(data))
        throw(BrainFlowError(string("Data Len must be power of two ", INVALID_ARGUMENTS_ERROR), INVALID_ARGUMENTS_ERROR))
    end

    temp_ampls = Vector{Float64}(undef, Integer(length(data) / 2) + 1)
    temp_freqs = Vector{Float64}(undef, Integer(length(data) / 2) + 1)

    ec = STATUS_OK
    # due to this bug https://github.com/JuliaLang/julia/issues/29602 libname should be hardcoded
    if Sys.iswindows()
        ec = ccall((:get_psd, "DataHandler.dll"), Cint, (Ptr{Float64}, Cint, Cint, Cint, Ptr{Float64}, Ptr{Float64}),
            data, length(data), Int32(sampling_rate), Int32(window), temp_ampls, temp_freqs)
    elseif Sys.isapple()
        ec = ccall((:get_psd, "libDataHandler.dylib"), Cint, (Ptr{Float64}, Cint, Cint, Cint, Ptr{Float64}, Ptr{Float64}),
            data, length(data), Int32(sampling_rate), Int32(window), temp_ampls, temp_freqs)
    else
        ec = ccall((:get_psd, "libDataHandler.so"), Cint, (Ptr{Float64}, Cint, Cint, Cint, Ptr{Float64}, Ptr{Float64}),
            data, length(data), Int32(sampling_rate), Int32(window), temp_ampls, temp_freqs)
    end
    if ec != Integer(STATUS_OK)
        throw(BrainFlowError(string("Error in get_psd ", ec), ec))
    end
    temp_ampls, temp_freqs
end

function get_psd_welch(data, nfft::AnyIntType, overlap::AnyIntType, sampling_rate::AnyIntType, window::AnyIntType)

    function is_power_of_two(value)
        (value != 0) && (value & (value - 1) == 0)
    end

    if !is_power_of_two(nfft)
        throw(BrainFlowError(string("nfft must be power of two ", INVALID_ARGUMENTS_ERROR), INVALID_ARGUMENTS_ERROR))
    end

    temp_ampls = Vector{Float64}(undef, Integer(nfft / 2) + 1)
    temp_freqs = Vector{Float64}(undef, Integer(nfft / 2) + 1)

    ec = STATUS_OK
    # due to this bug https://github.com/JuliaLang/julia/issues/29602 libname should be hardcoded
    if Sys.iswindows()
        ec = ccall((:get_psd_welch, "DataHandler.dll"), Cint, (Ptr{Float64}, Cint, Cint, Cint, Cint, Cint, Ptr{Float64}, Ptr{Float64}),
            data, length(data), Int32(nfft), Int32(overlap), Int32(sampling_rate), Int32(window), temp_ampls, temp_freqs)
    elseif Sys.isapple()
        ec = ccall((:get_psd_welch, "libDataHandler.dylib"), Cint, (Ptr{Float64}, Cint, Cint, Cint, Cint, Cint, Ptr{Float64}, Ptr{Float64}),
            data, length(data), Int32(nfft), Int32(overlap), Int32(sampling_rate), Int32(window), temp_ampls, temp_freqs)
    else
        ec = ccall((:get_psd_welch, "libDataHandler.so"), Cint, (Ptr{Float64}, Cint, Cint, Cint, Cint, Cint, Ptr{Float64}, Ptr{Float64}),
            data, length(data), Int32(nfft), Int32(overlap), Int32(sampling_rate), Int32(window), temp_ampls, temp_freqs)
    end
    if ec != Integer(STATUS_OK)
        throw(BrainFlowError(string("Error in get_psd_welch ", ec), ec))
    end
    temp_ampls, temp_freqs
end

function get_log_psd_welch(data, nfft::AnyIntType, overlap::AnyIntType, sampling_rate::AnyIntType, window::AnyIntType)

    function is_power_of_two(value)
        (value != 0) && (value & (value - 1) == 0)
    end

    if !is_power_of_two(nfft)
        throw(BrainFlowError(string("nfft must be power of two ", INVALID_ARGUMENTS_ERROR), INVALID_ARGUMENTS_ERROR))
    end

    temp_ampls = Vector{Float64}(undef, Integer(nfft / 2) + 1)
    temp_freqs = Vector{Float64}(undef, Integer(nfft / 2) + 1)

    ec = STATUS_OK
    # due to this bug https://github.com/JuliaLang/julia/issues/29602 libname should be hardcoded
    if Sys.iswindows()
        ec = ccall((:get_log_psd_welch, "DataHandler.dll"), Cint, (Ptr{Float64}, Cint, Cint, Cint, Cint, Cint, Ptr{Float64}, Ptr{Float64}),
            data, length(data), Int32(nfft), Int32(overlap), Int32(sampling_rate), Int32(window), temp_ampls, temp_freqs)
    elseif Sys.isapple()
        ec = ccall((:get_log_psd_welch, "libDataHandler.dylib"), Cint, (Ptr{Float64}, Cint, Cint, Cint, Cint, Cint, Ptr{Float64}, Ptr{Float64}),
            data, length(data), Int32(nfft), Int32(overlap), Int32(sampling_rate), Int32(window), temp_ampls, temp_freqs)
    else
        ec = ccall((:get_log_psd_welch, "libDataHandler.so"), Cint, (Ptr{Float64}, Cint, Cint, Cint, Cint, Cint, Ptr{Float64}, Ptr{Float64}),
            data, length(data), Int32(nfft), Int32(overlap), Int32(sampling_rate), Int32(window), temp_ampls, temp_freqs)
    end
    if ec != Integer(STATUS_OK)
        throw(BrainFlowError(string("Error in get_log_psd_welch ", ec), ec))
    end
    temp_ampls, temp_freqs
end

function get_log_psd(data, sampling_rate::AnyIntType, window::AnyIntType)

    function is_power_of_two(value)
        (value != 0) && (value & (value - 1) == 0)
    end

    if !is_power_of_two(length(data))
        throw(BrainFlowError(string("Data Len must be power of two ", INVALID_ARGUMENTS_ERROR), INVALID_ARGUMENTS_ERROR))
    end

    temp_ampls = Vector{Float64}(undef, Integer(length(data) / 2) + 1)
    temp_freqs = Vector{Float64}(undef, Integer(length(data) / 2) + 1)

    ec = STATUS_OK
    # due to this bug https://github.com/JuliaLang/julia/issues/29602 libname should be hardcoded
    if Sys.iswindows()
        ec = ccall((:get_log_psd, "DataHandler.dll"), Cint, (Ptr{Float64}, Cint, Cint, Cint, Ptr{Float64}, Ptr{Float64}),
            data, length(data), Int32(sampling_rate), Int32(window), temp_re, temp_im)
    elseif Sys.isapple()
        ec = ccall((:get_log_psd, "libDataHandler.dylib"), Cint, (Ptr{Float64}, Cint, Cint, Cint, Ptr{Float64}, Ptr{Float64}),
            data, length(data), Int32(sampling_rate), Int32(window), temp_re, temp_im)
    else
        ec = ccall((:get_log_psd, "libDataHandler.so"), Cint, (Ptr{Float64}, Cint, Cint, Cint, Ptr{Float64}, Ptr{Float64}),
            data, length(data), Int32(sampling_rate), Int32(window), temp_re, temp_im)
    end
    if ec != Integer(STATUS_OK)
        throw(BrainFlowError(string("Error in get_log_psd ", ec), ec))
    end
    temp_ampls, temp_freqs
end

function get_band_power(psd, freq_start::AnyDoubleType, freq_end::AnyDoubleType)
    band_power = Vector{Float64}(undef, 1)

    ec = STATUS_OK
    # due to this bug https://github.com/JuliaLang/julia/issues/29602 libname should be hardcoded
    if Sys.iswindows()
        ec = ccall((:get_band_power, "DataHandler.dll"), Cint, (Ptr{Float64}, Ptr{Float64}, Cint, Float64, Float64, Ptr{Float64}),
            psd[1], psd[2], length(psd[1]), Float64(freq_start), Float64(freq_end), band_power)
    elseif Sys.isapple()
        ec = ccall((:get_band_power, "libDataHandler.dylib"), Cint, (Ptr{Float64}, Ptr{Float64}, Cint, Float64, Float64, Ptr{Float64}),
            psd[1], psd[2], length(psd[1]), Float64(freq_start), Float64(freq_end), band_power)
    else
        ec = ccall((:get_band_power, "libDataHandler.so"), Cint, (Ptr{Float64}, Ptr{Float64}, Cint, Float64, Float64, Ptr{Float64}),
            psd[1], psd[2], length(psd[1]), Float64(freq_start), Float64(freq_end), band_power)
    end
    if ec != Integer(STATUS_OK)
        throw(BrainFlowError(string("Error in get_band_power ", ec), ec))
    end
    band_power[1]
end

max_stars_repo_name: pareeknikhil/brainflow

max_issues_repo_name: pareeknikhil/brainflow

max_forks_repo_name: pareeknikhil/brainflow


text:
OPT_PROBS = Dict()

#### Himmelblau

OPT_PROBS["himmelblau"] = Dict()
OPT_PROBS["himmelblau"]["array"] = Dict()

function himmelblau!(x)
    fx = (x[1]^2 + x[2] - 11)^2 + (x[1] + x[2]^2 - 7)^2
    return fx
end
function himmelblau_batched_f!(X)
    F = map(himmelblau!, X)
    return F
end
function himmelblau_batched_f!(F, X)
    map!(himmelblau!, F, X)
    return F
end
function himmelblau_g!(∇f, x)
    ∇f[1] = 4.0 * x[1]^3 + 4.0 * x[1] * x[2] -
        44.0 * x[1] + 2.0 * x[1] + 2.0 * x[2]^2 - 14.0
    ∇f[2] = 2.0 * x[1]^2 + 2.0 * x[2] - 22.0 +
        4.0 * x[1] * x[2] + 4.0 * x[2]^3 - 28.0 * x[2]
    ∇f
end

function himmelblau_fg!(∇f, x)
    ∇f = himmelblau_g!(∇f, x)
    fx = himmelblau!(x)
    return fx, ∇f
end
function himmelblau_fgh!(∇f, ∇²f, x)
    ∇²f = himmelblau_h!(∇²f, x)
    fx, ∇f = himmelblau_fg!(∇f, x)
    return fx, ∇f, ∇²f
end
function himmelblau_h!(∇²f, x)
    ∇²f[1, 1] = 12.0 * x[1]^2 + 4.0 * x[2] - 44.0 + 2.0
    ∇²f[1, 2] = 4.0 * x[1] + 4.0 * x[2]
    ∇²f[2, 1] = ∇²f[1, 2]
    ∇²f[2, 2] = 2.0 + 4.0 * x[1] + 12.0 * x[2]^2 - 28.0
    return ∇²f
end
function himmelblau_hv!(hv, x, v)
    hv[1] = (12.0 * x[1]^2 + 4.0 * x[2] - 44.0 + 2.0)*v[1] + (4.0 * x[1] + 4.0 * x[2])*v[2]
    hv[2] =  (4.0 * x[1] + 4.0 * x[2])*v[1] + (2.0 + 4.0 * x[1] + 12.0 * x[2]^2 - 28.0)*v[2]
    return hv
end

OPT_PROBS["himmelblau"]["array"]["x0"] = [3.0, 1.0]
OPT_PROBS["himmelblau"]["array"]["mutating"] = ScalarObjective(himmelblau!, himmelblau_g!, himmelblau_fg!, himmelblau_fgh!, himmelblau_h!, himmelblau_hv!, himmelblau_batched_f!, nothing)

### Exponential

exponential!(x) = exp((2.0 - x[1])^2) + exp((3.0 - x[2])^2)
function exponential_g!(g, x)
    g[1] = -2.0 * (2.0 - x[1]) * exp((2.0 - x[1])^2)
    g[2] = -2.0 * (3.0 - x[2]) * exp((3.0 - x[2])^2)
    return g
end
function exponential_h!(H, x)
    H[1, 1] = 2.0 * exp((2.0 - x[1])^2) * (2.0 * x[1]^2 - 8.0 * x[1] + 9)
    H[1, 2] = 0.0
    H[2, 1] = 0.0
    H[2, 2] = 2.0 * exp((3.0 - x[2])^2) * (2.0 * x[2]^2 - 12.0 * x[2] + 19)
    return H
end
function exponential_hv!(Hv, x, v)
    Hv[1, 1] = (2.0 * exp((2.0 - x[1])^2) * (2.0 * x[1]^2 - 8.0 * x[1] + 9))*v[1]
    Hv[2, 2] = (2.0 * exp((3.0 - x[2])^2) * (2.0 * x[2]^2 - 12.0 * x[2] + 19))*v[2]
    return Hv
end
function exponential_fg!(g, x)
    fx = exponential!(x)
    g = exponential_g!(g, x)
    return fx, g
end
function exponential_fgh!(g, H, x)
    fx, g = exponential_fg!(g, x)
    H = exponential_h!(H, x)
    return fx, g, H
end

OPT_PROBS["exponential"] = Dict()
OPT_PROBS["exponential"]["array"] = Dict()
# Byttet om på x og H
OPT_PROBS["exponential"]["array"]["x0"] = [0.0, 0.0]
OPT_PROBS["exponential"]["array"]["mutating"] = ScalarObjective(exponential!, exponential_g!, exponential_fg!, exponential_fgh!, exponential_h!, exponential_hv!, nothing, nothing)

exponential_prb = Dict()
exponential_prb["twicediffed!"] = TwiceDiffed(exponential!)
exponential_prb["oncediffed!"] = OnceDiffed(exponential!)
exponential_prb["initial_x"] = [0.0, 0.0]
exponential_prb["minimizer"] = [2.0, 3.0]
exponential_prb["minimum"] = exponential!(nothing, nothing, [2.0, 3.0])
problems["unconstrained"]["exponential"] = exponential_prb

ext: jl


text:
const OPENSSL_LIB = "libcrypto"

md5init() = ccall((:OpenSSL_add_all_digests, OPENSSL_LIB), Void, ())
md5cleanup() = ccall((:EVP_cleanup, OPENSSL_LIB), Void, ())

function md5(v...)
    iob = IOBuffer()
    write(iob, v...)
    ctx = ccall((:EVP_MD_CTX_create, OPENSSL_LIB), Ptr{Void}, ())
    try
        md = ccall((:EVP_get_digestbyname, OPENSSL_LIB), Ptr{Void}, (Ptr{UInt8},), "MD5")
        (md == C_NULL) && error("Unknown message digest $name")

        ccall((:EVP_DigestInit_ex, OPENSSL_LIB), Void, (Ptr{Void}, Ptr{Void}, Ptr{Void}), ctx, md, C_NULL)

        data = takebuf_array(iob)
        ccall((:EVP_DigestUpdate, OPENSSL_LIB), Void, (Ptr{Void}, Ptr{UInt8}, UInt), ctx, data, length(data))

        size = ccall((:EVP_MD_size, OPENSSL_LIB), UInt, (Ptr{Void},), md)
        uval = Array(UInt8, size)

        ccall((:EVP_DigestFinal_ex, OPENSSL_LIB), Void, (Ptr{Void}, Ptr{UInt8}, Ptr{UInt}), ctx, uval, C_NULL)

        return uval
    finally
        ccall((:EVP_MD_CTX_destroy, OPENSSL_LIB), Void, (Ptr{Void},), ctx)
    end
end


text:
# DGA(q::Float64, m2::Float64, hatZA::Float64, tildeZA::Float64, DtildeZA::Float64, t::Float64) =
#     -1 *
#     hatZA *(hatZA^2+δ)^-1 *
#     (m2 + (q^2) * (Rb(exp(-2t) * Λ^-2 * (q^2)) + 1))^-2 *
#     (q^2) *
#     (
#         -2 * exp(-2t) * Λ^-2 * (q^2) * Rbp(exp(-2t) * Λ^-2 * (q^2)) +
#         DtildeZA * tildeZA *(tildeZA^2+δ)^-1 * Rb(exp(-2 * t) * Λ^-2 * (q^2))
#     )

invhatZa(
    ZA::Interpolations.Extrapolation,
    x::Float64,
    m2::Float64,
    t::Float64,
) =
    ZA(max(1, sqrt(x + δ * exp(2t) * Λ^2)))^-1 *
    (1 + m2 / (x + δ * exp(2t) * Λ^2))

invhatZa(
    ZA::Interpolations.Extrapolation,
    x::Vector{Float64},
    m2::Float64,
    t::Float64,
) =
    ZA(max.(1, sqrt.(x .+ δ * exp(2t) * Λ^2))) .^ -1 .*
    (1 .+ m2 ./ (x .+ δ * exp(2t) * Λ^2))

DGA(
    q::Float64,
    m2::Float64,
    ZA::Interpolations.Extrapolation,
    tildeZA::Float64,
    DtildeZA::Float64,
    t::Float64,
) =
    -1 *
    ZA(sqrt(q^2 + δ * exp(2t) * Λ^2))^-1 *
    (1 + m2 / (q^2 + δ * exp(2t) * Λ^2)) *
    (m2 + (q^2) * (Rb(exp(-2t) * Λ^-2 * (q^2)) + 1))^-2 *
    (q^2) *
    (
        -2 * exp(-2t) * Λ^-2 * (q^2) * Rbp(exp(-2t) * Λ^-2 * (q^2)) +
        DtildeZA * tildeZA^-1 * Rb(exp(-2 * t) * Λ^-2 * (q^2))
    )

DGc(
    q::Float64,
    Zc::Interpolations.Extrapolation,
    tildeZc::Float64,
    DtildeZc::Float64,
    t::Float64,
) =
    -1 *
    Zc(sqrt(q^2+ δ * Λ^2 * exp(2 * t)))^-1 *
    (q^2)^-1 *
    (1 + Rb(exp(-2t) * Λ^-2 * (q^2)))^-2 *
    (
        -2 * exp(-2t) * Λ^-2 * (q^2) * Rbp(exp(-2t) * Λ^-2 * (q^2)) +
        DtildeZc * tildeZc^-1 * Rb(exp(-2 * t) * Λ^-2 * (q^2))
    )

# GAqmp(q::Float64, p::Float64, costh::Float64, m2::Float64, hatZAqmp::Float64, t::Float64) =
#     hatZAqmp *(hatZAqmp^2+δ)^-1  *
#     (m2 + (1 + Rb(exp(-2t) * Λ^-2 * Qmp(q, p, costh))) * Qmp(q, p, costh))^-1

GAqmp(
    q::Float64,
    p::Vector{Float64},
    costh::Float64,
    m2::Float64,
    ZA::Interpolations.Extrapolation,
    t::Float64,
) =
    ZA(max.(sqrt.(Qmp(q, p, costh) .+ δ * exp(2t) * Λ^2), 1)) .^ -1 .*
    (1 .+ m2 ./ max.(Qmp(q, p, costh) .+ δ * exp(2t) * Λ^2, 1)) .*
    (
        m2 .+
        (1 .+ Rb(exp(-2t) * Λ^-2 .* Qmp(q, p, costh))) .* Qmp(q, p, costh)
    ) .^ -1

GAqmp(
    q::Float64,
    p::Float64,
    costh::Float64,
    m2::Float64,
    ZA::Interpolations.Extrapolation,
    t::Float64,
) =
    ZA(max(sqrt(Qmp(q, p, costh) + δ * exp(2t) * Λ^2), 1))^-1 *
    (1 + m2 / max(Qmp(q, p, costh) + δ * exp(2t) * Λ^2, 1)) *
    (m2 + (1 + Rb(exp(-2t) * Λ^-2 * Qmp(q, p, costh))) * Qmp(q, p, costh))^-1

Gcqmp(
    q::Float64,
    p::Vector{Float64},
    costh::Float64,
    Zc::Interpolations.Extrapolation,
    t::Float64,
) =
    Zc(max.(sqrtQmp(q, p, costh), 1)) .^ -1 .*
    (1 .+ Rb(exp(-2t) .* Λ^-2 .* Qmp(q, p, costh))) .^ -1 .*
    Qmp(q, p, costh) .^ -1

Gcqmp(
    q::Float64,
    p::Float64,
    costh::Float64,
    Zc::Interpolations.Extrapolation,
    t::Float64,
) =
    Zc(max(sqrt(Qmp(q, p, costh)+ δ * Λ^2 * exp(2 * t)), 1))^-1 *
    (1 + Rb(exp(-2t) * Λ^-2 * Qmp(q, p, costh)))^-1 *
    Qmp(q, p, costh)^-1

Gc(
    q::Float64,
    p::Vector{Float64},
    x1::Float64,
    x2::Float64,
    Zc::Interpolations.Extrapolation,
    Qpf::Function,
    t::Float64,
) =
    Zc(max.(Qpf(q, p, x1, x2), 1)) .^ -1 .*
    (1 .+ Rb(exp(-2t) .* Λ^-2 .* Qpf(q, p, x1, x2))) .^ -1 .*
    Qpf(q, p, x1, x2) .^ -1

Gc(
    q::Float64,
    p::Float64,
    x1::Float64,
    x2::Float64,
    Zc::Interpolations.Extrapolation,
    Qpf::Function,
    t::Float64,
) =
    Zc(max(Qpf(q, p, x1, x2), 1))^-1 *
    (1 + Rb(exp(-2t) * Λ^-2 * Qpf(q, p, x1, x2)))^-1 *
    Qpf(q, p, x1, x2)^-1

max_stars_repo_name: Yangyang-Tan/FRG-Code-Collection

max_issues_repo_name: Yangyang-Tan/FRG-Code-Collection

max_forks_repo_name: Yangyang-Tan/FRG-Code-Collection


text:
struct Bezier <: ReferenceFEName end

const bezier = Bezier()

struct BezierRefFE{D} <: ReferenceFE{D}
  reffe::ReferenceFE{D}
  node_to_own_node::Vector{Int}
  shapefuns::AbstractVector{<:Field}
  metadata
end

function BezierRefFE(::Type{T},p::Polytope{D},orders) where {D,T}
  reffe = LagrangianRefFE(T,p,orders)
  nodes = get_node_coordinates(reffe)
  prebasis = get_prebasis(reffe)
  node_to_own_node = compute_node_to_bezier_node(prebasis,nodes)
  shapefuns = berstein_basis(prebasis,p)
  shapefuns = collect(lazy_map(Reindex(shapefuns),node_to_own_node))
  lagreffaces = reffe.reffe.metadata
  reffaces = _convert_reffaces(T,lagreffaces)
  tuple( reffaces ... )
  BezierRefFE(reffe,node_to_own_node,shapefuns,reffaces)
end

function ReferenceFE(
  polytope::Polytope,
  ::Bezier,
  ::Type{T},
  orders::Union{Integer,Tuple{Vararg{Integer}}}) where T

  BezierRefFE(T,polytope,orders)
end

get_node_coordinates(reffe::BezierRefFE) = get_node_coordinates(reffe.reffe)

get_node_to_own_node(reffe::BezierRefFE) = reffe.node_to_own_node

get_metadata(reffe::BezierRefFE) = reffe.metadata

get_shapefuns(reffe::BezierRefFE) = reffe.shapefuns

get_prebasis(reffe::BezierRefFE) = reffe.prebasis

function (==)(a::BezierRefFE{D},b::BezierRefFE{D}) where D
  t = true
  a.reffe == b.reffe
  na = get_node_to_own_node(a)
  nb = get_node_to_own_node(a)
  t = t && (na == nb)
  t
end

function (==)(a::BezierRefFE,b::BezierRefFE)
  false
end

function _convert_reffaces(::Type{T},reffaces) where T
  _reffaces = []
  for d_reffaces in reffaces
    _d_reffaces = [ _lagrange_to_bezier_fe(T,reffe) for reffe in d_reffaces ]
    push!(_reffaces,_d_reffaces)
  end
  tuple(_reffaces...)
end

function _lagrange_to_bezier_fe(::Type{T},reffe::LagrangianRefFE) where T
  BezierRefFE(T,get_polytope(reffe),get_orders(reffe))
end

function compute_node_to_bezier_node(prebasis::MonomialBasis{D,T},nodes) where {D,T}
  orders = get_orders(prebasis)
  terms = _coords_to_terms(nodes,orders)
  _prebasis = MonomialBasis{D}(T,orders,terms)
  _exps = get_exponents(_prebasis)
  exps = get_exponents(prebasis)
  [ findfirst( isequal(i), exps) for i in _exps ]
end

## Bernstein Basis

function _bernstein_term(p,a,i)
  @assert i ≤ p
  @assert a ≤ p
  if ( i ≤ a ≤ p )
    f = factorial(p) ÷ ( factorial(i)*factorial(p-i) )
    b1 = binomial( p-i, a-i )
    s = (-1)^(a-i)
    f*b1*s
  else
    0
  end
end

function _bernstein_term(p,a,b,i,j)
  @assert i+j ≤ p
  @assert a+b ≤ p
  if ( i ≤ a ≤ p-j ) && ( j ≤ b ≤ p-a )
    f = factorial(p) ÷ ( factorial(i)*factorial(j)*factorial(p-i-j) )
    b1 = binomial( p-i-j, a-i )
    b2 = binomial( p-a-j, b-j )
    s = (-1)^(a+b-i-j)
    f*b1*b2*s
  else
    0
  end
end

function _bernstein_term(p,a,b,c,i,j,k)
  @assert i+j+k ≤ p
  @assert a+b+c ≤ p
  if ( i ≤ a ≤ p-j-k ) && ( j ≤ b ≤ p-a-k ) && ( k ≤ c ≤ p-a-b )
    p!,i!,j!,k! = factorial(p),factorial(i),factorial(j),factorial(k)
    f = p! ÷ ( i!*j!*k!*factorial(p-i-j-k) )
    b1 = binomial( p-i-j-k, a-i )
    b2 = binomial( p-a-j-k, b-j )
    b3 = binomial( p-a-b-k, c-k )
    s = (-1)^(a+b+c-i-j-k)
    f*b1*b2*b3*s
  else
    0
  end
end

function _bernstein_term_n_cube(p,a,b,i,j)
  _bernstein_term(p,a,i)*_bernstein_term(p,b,j)
end

function _bernstein_term_n_cube(p,a,b,c,i,j,k)
  _bernstein_term(p,a,i)*_bernstein_term(p,b,j)*_bernstein_term(p,c,k)
end

function _bernstein_term(
  p::Polytope,
  order::Integer,
  a::NTuple{N},
  i::NTuple{N}) where N

  args = (order,a...,i...)
  if is_simplex(p)
    _bernstein_term(args...)
  else
    _bernstein_term_n_cube(args...)
  end
end

function _berstein_matrix(prebasis,p::Polytope)
  o = get_order(prebasis)
  e = get_exponents(prebasis)
  M = CartesianIndices( (length(e),length(e)) )
  collect(lazy_map( i-> _bernstein_term(p,o,e[i.I[1]],e[i.I[2]]), M))
end

function berstein_basis(prebasis,polytope::Polytope)
  C = _berstein_matrix(prebasis,polytope)
  linear_combination( C, prebasis )
end

berstein_basis(::MonomialBasis{0},::Polytope) = [ConstantField(1)]

function rationalize_bernstein_basis(basis,weights)
  ( basis .* weights ) ./ ( basis ⋅ weights )
end


text:
using .ImageIO

function log_image(lg::MLFLogger, key::AbstractString, obj::AbstractArray{<:Colorant}; step=nothing)
    mktempdir(;prefix="jlmlf_") do dirname
        fname = joinpath(dirname, key)
        fname = joinpath(fname, "image.png")
        FileIO.save(fname, obj)
        log_image(lg, key, fname; step=step)
    end
end

process(lg::MLFLogger, name::AbstractString, obj::AbstractArray{<:Colorant}, step::Int) = log_image(lg, name, obj; step=step)


text:
function default_in_partition(sites::Tuple{Int}, p::Integer, nparts::Integer)
  return p == mod1(sites[1], nparts)
end

# i ≥ j
_f(i, j) = (i - 1) * i ÷ 2 + j
function default_in_partition(sites::NTuple{2,Int}, p, nparts)
  i, j = sites
  if i ≤ j
    return p == mod1(_f(j, i), nparts)
  end
  return p == mod1(_f(i, j), nparts)
end

function default_in_partition_sorted(sites::NTuple{4,Int}, p, nparts)
  i, j, k, l = sites
  if j == k
    return default_in_partition((j,), p, nparts)
  end
  return default_in_partition((i, j), p, nparts)
end

function default_in_partition(sites::NTuple{4,Int}, p, nparts)
  return default_in_partition_sorted(sort(sites), p, nparts)
end

max_stars_repo_name: mtfishman/ITensorParallel

max_issues_repo_name: mtfishman/ITensorParallel

max_forks_repo_name: mtfishman/ITensorParallel


text:
using Pkg
Pkg.activate("GraphFusedElasticNet.jl")
# Pkg.instantiate()

using DataFrames
using Printf
using DelimitedFiles
using JSON
using NPZ
using GraphFusedElasticNet

#
debug = false
if debug
    split_to_fit = 0
else
    split_to_fit =  parse(Int, ENV["SLURM_PROCID"])
end

##
splits = readtable("processed_data/splits_qua.csv")
split = splits[split_to_fit + 1, :]

##
fname = @sprintf("productivity_splits/%02d.csv", split_to_fit)
splitdata = readdlm(fname, ',')
println("Fitting to data $(fname)...")

## find be[s]t lambda
fname = @sprintf("modelfit_metrics/cvloss_%02d.csv", split_to_fit)
best_lams = readtable(fname)
best_row = argmax(best_lams.final_pred)
λsl1 = best_lams.λsl1[best_row]
λsl2 = best_lams.λsl2[best_row]
λtl1 = best_lams.λtl1[best_row]
λtl2 = best_lams.λtl2[best_row]

##
edges_df = readtable("../gfen-reproduce/processed_data/spatiotemporal_graph.csv")

##
num_nodes = size(splitdata, 1)
num_edges = size(edges_df, 1)

## load data and map estimate for fast mixing
ϵ = 1e-3
s = splitdata[:, 1] .+ ϵ
a = splitdata[:, 2] .+ 2ϵ

## fit MAP model
trails = JSON.parsefile("processed_data/trails.json")
istemp = Bool.(trails["istemp"])
wts = Float64.(trails["wts"])
ptr = Int.(trails["ptr"])
brks = Int.(trails["brks"])
lambdasl1 = [w * (t ? λtl1 : λsl1) for (t, w) in zip(istemp, wts)]
lambdasl2 = [w * (t ? λtl2 : λsl2) for (t, w) in zip(istemp, wts)]

##
modelopts = Dict{Symbol, Any}(
    :admm_balance_every => 20,
    :admm_init_penalty => 1.0,
    :admm_residual_balancing => true,
    :admm_adaptive_inflation => true,
    :reltol => 1e-4,
    :admm_min_penalty => 0.25,
    :admm_max_penalty => 32.0,
    :abstol => 1e-7,
    :save_norms => false,
    :save_loss => false
)
fitopts = Dict{Symbol, Any}(
    :walltime => 3600.0,
    :parallel => true
)

## define and train model
fname = @sprintf("../gfen-reproduce/best_betas/betas_%02d.csv", split_to_fit)

println("Fitting MAP model...")
map_mod = BinomialGFEN(ptr, brks, lambdasl1, lambdasl2; modelopts...)
@time converged = fit!(map_mod, s, a ; fitopts...)
println("MAP estimate converged: $(converged)")
println("Saving MAP estimate to $fname... in Float16")
open(fname, "w") do io
    writedlm(io, Float16.(map_mod.beta), ',')
end

##

println("Loading MAP model...")
mcmc_init = vec(readdlm(fname, ','))

##

edges = [(r.vertex1 + 1, r.vertex2 + 1) for r in eachrow(edges_df)]
tv1 = [(r.temporal == 1) ? λtl1 : λsl1 for r in eachrow(edges_df)]
tv2 = [(r.temporal == 1) ? λtl2 : λsl2 for r in eachrow(edges_df)]

println("Initializing MCMC chain...")
bmod = BayesianBinomialGFEN(edges, tv1=tv1, tv2=tv2)

## fit model
# runs took anytime between 3 and 7 hours
n = 5500
thinning = 5
burnin = 4500

# n = 200
# thinning = 1
# burnin = 0

#
init = mcmc_init  # zeros(size(mcmc_init))
@time chain = sample_chain(bmod, s, a, n, burnin=burnin, thinning=thinning, init=init, verbose=false, async=true)

##
fname = @sprintf("best_betas_bayesian/%02d.npy", split_to_fit)
println("Saving chain $(fname) in Float16...")

# open(fname, "w") do io
#     writedlm(io, Float16.(chain), ',')
# end
npzwrite(fname, Float16.(chain))

max_stars_repo_name: mauriciogtec/gfen-reproduce

max_issues_repo_name: mauriciogtec/gfen-reproduce

max_forks_repo_name: mauriciogtec/gfen-reproduce


text:
####
#### Implementation for Hamiltonian dynamics trajectories
####
#### Developers' Notes
####
#### Not all functions that use `rng` require a fallback function with `GLOBAL_RNG`
#### as default. In short, only those exported to other libries need such a fallback
#### function. Internal uses shall always use the explict `rng` version. (Kai Xu 6/Jul/19)

"""
$(TYPEDEF)
A transition that contains the phase point and
other statistics of the transition.

# Fields
$(TYPEDFIELDS)
"""
struct Transition{P<:PhasePoint, NT<:NamedTuple}
    "Phase-point for the transition."
    z       ::  P
    "Statistics related to the transition, e.g. energy."
    stat    ::  NT
end

"Returns the statistics for transition `t`."
stat(t::Transition) = t.stat

abstract type AbstractMCMCKernel end

abstract type AbstractTerminationCriterion end

abstract type StaticTerminationCriterion <: AbstractTerminationCriterion end

abstract type DynamicTerminationCriterion <: AbstractTerminationCriterion end

"""
$(TYPEDEF)
Static HMC with a fixed number of leapfrog steps.

# Fields
$(TYPEDFIELDS)

# References
1. Neal, R. M. (2011). MCMC using Hamiltonian dynamics. Handbook of Markov chain Monte Carlo, 2(11), 2. ([arXiv](https://arxiv.org/pdf/1206.1901))
"""
struct FixedNSteps <: StaticTerminationCriterion
    "Number of steps to simulate, i.e. length of trajectory will be `L + 1`."
    L::Int
end

"""
$(TYPEDEF)
Standard HMC implementation with a fixed integration time.

# Fields
$(TYPEDFIELDS)

# References
1. Neal, R. M. (2011). MCMC using Hamiltonian dynamics. Handbook of Markov chain Monte Carlo, 2(11), 2. ([arXiv](https://arxiv.org/pdf/1206.1901))
"""
struct FixedIntegrationTime{F<:AbstractFloat} <: StaticTerminationCriterion
    "Total length of the trajectory, i.e. take `floor(λ / integrator_step_size)` number of leapfrog steps."
    λ::F
end

##
## Sampling methods for trajectories.
##

"How to sample a phase-point from the simulated trajectory."
abstract type AbstractTrajectorySampler end

"Samples the end-point of the trajectory."
struct EndPointTS <: AbstractTrajectorySampler end

"""
$(TYPEDEF)

Trajectory slice sampler carried during the building of the tree.
It contains the slice variable and the number of acceptable condidates in the tree.

# Fields

$(TYPEDFIELDS)
"""
struct SliceTS{F<:AbstractFloat} <: AbstractTrajectorySampler
    "Sampled candidate `PhasePoint`."
    zcand   ::  PhasePoint
    "Slice variable in log-space."
    ℓu      ::  F
    "Number of acceptable candidates, i.e. those with probability larger than slice variable `u`."
    n       ::  Int
end

Base.show(io::IO, s::SliceTS) = print(io, "SliceTS(ℓu=$(s.ℓu), n=$(s.n))")

"""
$(TYPEDEF)

Multinomial trajectory sampler carried during the building of the tree.
It contains the weight of the tree, defined as the total probabilities of the leaves.

# Fields

$(TYPEDFIELDS)
"""
struct MultinomialTS{F<:AbstractFloat} <: AbstractTrajectorySampler
    "Sampled candidate `PhasePoint`."
    zcand   ::  PhasePoint
    "Total energy for the given tree, i.e. the sum of energies of all leaves."
    ℓw      ::  F
end

"""
$(TYPEDEF)

Slice sampler for the starting single leaf tree.
Slice variable is initialized.
"""
SliceTS(rng::AbstractRNG, z0::PhasePoint) = SliceTS(z0, log(rand(rng)) - energy(z0), 1)

"""
$(TYPEDEF)

Multinomial sampler for the starting single leaf tree.
(Log) weights for leaf nodes are their (unnormalised) Hamiltonian energies.

Ref: https://github.com/stan-dev/stan/blob/develop/src/stan/mcmc/hmc/nuts/base_nuts.hpp#L226
"""
MultinomialTS(rng::AbstractRNG, z0::PhasePoint) = MultinomialTS(z0, zero(energy(z0)))

"""
$(TYPEDEF)

Create a slice sampler for a single leaf tree:
- the slice variable is copied from the passed-in sampler `s` and
- the number of acceptable candicates is computed by comparing the slice variable against the current energy.
"""
function SliceTS(s::SliceTS, H0::AbstractFloat, zcand::PhasePoint)
    return SliceTS(zcand, s.ℓu, (s.ℓu <= -energy(zcand)) ? 1 : 0)
end

"""
$(TYPEDEF)

Multinomial sampler for a trajectory consisting only a leaf node.
- tree weight is the (unnormalised) energy of the leaf.
"""
function MultinomialTS(s::MultinomialTS, H0::AbstractFloat, zcand::PhasePoint)
    return MultinomialTS(zcand, H0 - energy(zcand))
end

function combine(rng::AbstractRNG, s1::SliceTS, s2::SliceTS)
    @assert s1.ℓu == s2.ℓu "Cannot combine two slice sampler with different slice variable"
    n = s1.n + s2.n
    zcand = rand(rng) < s1.n / n ? s1.zcand : s2.zcand
    return SliceTS(zcand, s1.ℓu, n)
end

function combine(zcand::PhasePoint, s1::SliceTS, s2::SliceTS)
    @assert s1.ℓu == s2.ℓu "Cannot combine two slice sampler with different slice variable"
    n = s1.n + s2.n
    return SliceTS(zcand, s1.ℓu, n)
end

function combine(rng::AbstractRNG, s1::MultinomialTS, s2::MultinomialTS)
    ℓw = logaddexp(s1.ℓw, s2.ℓw)
    zcand = rand(rng) < exp(s1.ℓw - ℓw) ? s1.zcand : s2.zcand
    return MultinomialTS(zcand, ℓw)
end

function combine(zcand::PhasePoint, s1::MultinomialTS, s2::MultinomialTS)
    ℓw = logaddexp(s1.ℓw, s2.ℓw)
    return MultinomialTS(zcand, ℓw)
end

mh_accept(rng::AbstractRNG, s::SliceTS, s′::SliceTS) = rand(rng) < min(1, s′.n / s.n)

function mh_accept(rng::AbstractRNG, s::MultinomialTS, s′::MultinomialTS)
    return rand(rng) < min(1, exp(s′.ℓw - s.ℓw))
end

"""
$(TYPEDEF)

Numerically simulated Hamiltonian trajectories.
"""
struct Trajectory{
    TS<:AbstractTrajectorySampler,
    I<:AbstractIntegrator,
    TC<:AbstractTerminationCriterion,
}
    "Integrator used to simulate trajectory."
    integrator::I
    "Criterion to terminate the simulation."
    termination_criterion::TC
end

Trajectory{TS}(integrator::I, termination_criterion::TC) where {TS, I, TC} =
    Trajectory{TS, I, TC}(integrator, termination_criterion)

ConstructionBase.constructorof(::Type{<:Trajectory{TS}}) where {TS} = Trajectory{TS}

function Base.show(io::IO, τ::Trajectory{TS}) where {TS}
    print(io, "Trajectory{$TS}(integrator=$(τ.integrator), tc=$(τ.termination_criterion))")
end

nsteps(τ::Trajectory{TS, I, TC}) where {TS, I, TC<:FixedNSteps} = τ.termination_criterion.L
nsteps(τ::Trajectory{TS, I, TC}) where {TS, I, TC<:FixedIntegrationTime} =
    max(1, floor(Int, τ.termination_criterion.λ / nom_step_size(τ.integrator)))

##
## Kernel interface
##

struct HMCKernel{R, T<:Trajectory} <: AbstractMCMCKernel
    refreshment::R
    τ::T
end

HMCKernel(τ::Trajectory) = HMCKernel(FullMomentumRefreshment(), τ)

"""
$(SIGNATURES)

Make a MCMC transition from phase point `z` using the trajectory `τ` under Hamiltonian `h`.

NOTE: This is a RNG-implicit fallback function for `transition(GLOBAL_RNG, τ, h, z)`
"""
function transition(τ::Trajectory, h::Hamiltonian, z::PhasePoint)
    return transition(GLOBAL_RNG, τ, h, z)
end

###
### Actual trajectory implementations
###

function transition(
    rng::Union{AbstractRNG, AbstractVector{<:AbstractRNG}},
    τ::Trajectory{TS, I, TC},
    h::Hamiltonian,
    z::PhasePoint,
) where {TS<:AbstractTrajectorySampler, I, TC<:StaticTerminationCriterion}
    H0 = energy(z)

    z′, is_accept, α = sample_phasepoint(rng, τ, h, z)
    # Do the actual accept / reject
    z = accept_phasepoint!(z, z′, is_accept)    # NOTE: this function changes `z′` in place in matrix-parallel mode
    # Reverse momentum variable to preserve reversibility
    z = PhasePoint(z.θ, -z.r, z.ℓπ, z.ℓκ)
    H = energy(z)
    tstat = merge(
        (
            n_steps=nsteps(τ),
            is_accept=is_accept,
            acceptance_rate=α,
            log_density=z.ℓπ.value,
            hamiltonian_energy=H,
            hamiltonian_energy_error=H - H0,
        ),
        stat(τ.integrator),
    )
    return Transition(z, tstat)
end

# Return the accepted phase point
function accept_phasepoint!(z::T, z′::T, is_accept::Bool) where {T<:PhasePoint{<:AbstractVector}}
    if is_accept
        return z′
    else
        return z
    end
end
function accept_phasepoint!(z::T, z′::T, is_accept) where {T<:PhasePoint{<:AbstractMatrix}}
    # Revert unaccepted proposals in `z′`
    is_reject = (!).(is_accept)
    if any(is_reject)
        # Convert logical indexing to number indexing to support CUDA.jl
        # NOTE: for x::CuArray, x[:,Vector{Bool}]  is NOT supported
        #                       x[:,CuVector{Int}] is NOT supported
        #                       x[:,Vector{Int}]   is     supported
        is_reject = findall(is_reject) |> Array
        z′.θ[:,is_reject] = z.θ[:,is_reject]
        z′.r[:,is_reject] = z.r[:,is_reject]
        z′.ℓπ.value[is_reject] = z.ℓπ.value[is_reject]
        z′.ℓπ.gradient[:,is_reject] = z.ℓπ.gradient[:,is_reject]
        z′.ℓκ.value[is_reject] = z.ℓκ.value[is_reject]
        z′.ℓκ.gradient[:,is_reject] = z.ℓκ.gradient[:,is_reject]
    end
    # Always return `z′` as any unaccepted proposal is already reverted
    # NOTE: This in place treatment of `z′` is for memory efficient consideration.
    #       We can also copy `z′ and avoid mutating the original `z′`. But this is
    #       not efficient and immutability of `z′` is not important in this local scope.
    return z′
end

### Use end-point from the trajectory as a proposal and apply MH correction

function sample_phasepoint(rng, τ::Trajectory{EndPointTS}, h, z)
    z′ = step(τ.integrator, h, z, nsteps(τ))
    is_accept, α = mh_accept_ratio(rng, energy(z), energy(z′))
    return z′, is_accept, α
end

### Multinomial sampling from trajectory

function randcat(rng::AbstractRNG, zs::AbstractVector{<:PhasePoint}, unnorm_ℓp::AbstractVector)
    p = exp.(unnorm_ℓp .- logsumexp(unnorm_ℓp))
    i = randcat(rng, p)
    return zs[i]
end

# zs is in the form of Vector{PhasePoint{Matrix}} and has shape [n_steps][dim, n_chains]
function randcat(rng, zs::AbstractVector{<:PhasePoint}, unnorm_ℓP::AbstractMatrix)
    z = similar(first(zs))
    P = exp.(unnorm_ℓP .- logsumexp(unnorm_ℓP; dims=2)) # (n_chains, n_steps)
    is = randcat(rng, P')
    foreach(enumerate(is)) do (i_chain, i_step)
        zi = zs[i_step]
        z.θ[:,i_chain] = zi.θ[:,i_chain]
        z.r[:,i_chain] = zi.r[:,i_chain]
        z.ℓπ.value[i_chain] = zi.ℓπ.value[i_chain]
        z.ℓπ.gradient[:,i_chain] = zi.ℓπ.gradient[:,i_chain]
        z.ℓκ.value[i_chain] = zi.ℓκ.value[i_chain]
        z.ℓκ.gradient[:,i_chain] = zi.ℓκ.gradient[:,i_chain]
    end
    return z
end

function sample_phasepoint(rng, τ::Trajectory{MultinomialTS}, h, z)
    n_steps = abs(nsteps(τ))
    # TODO: Deal with vectorized-mode generically.
    #       Currently the direction of multiple chains are always coupled
    n_steps_fwd = rand_coupled(rng, 0:n_steps)
    zs_fwd = step(τ.integrator, h, z, n_steps_fwd; fwd=true, full_trajectory=Val(true))
    n_steps_bwd = n_steps - n_steps_fwd
    zs_bwd = step(τ.integrator, h, z, n_steps_bwd; fwd=false, full_trajectory=Val(true))
    zs = vcat(reverse(zs_bwd)..., z, zs_fwd...)
    ℓweights = -energy.(zs)
    if eltype(ℓweights) <: AbstractVector
        ℓweights = cat(ℓweights...; dims=2)
    end
    unnorm_ℓprob = ℓweights
    z′ = randcat(rng, zs, unnorm_ℓprob)
    # Computing adaptation statistics for dual averaging as done in NUTS
    Hs = -ℓweights
    ΔH = Hs .- energy(z)
    α = exp.(min.(0, -ΔH))  # this is a matrix for vectorized mode and a vector otherwise
    α = typeof(α) <: AbstractVector ? mean(α) : vec(mean(α; dims=2))
    return z′, true, α
end

###
### Advanced HMC implementation with (adaptive) dynamic trajectory length.
###

##
## Variants of no-U-turn criteria
##

"""
$(TYPEDEF)
Classic No-U-Turn criterion as described in Eq. (9) in [1].

Informally, this will terminate the trajectory expansion if continuing
the simulation either forwards or backwards in time will decrease the
distance between the left-most and right-most positions.

# Fields
$(TYPEDFIELDS)

# References
1. Hoffman, M. D., & Gelman, A. (2014). The No-U-Turn Sampler: adaptively setting path lengths in Hamiltonian Monte Carlo. Journal of Machine Learning Research, 15(1), 1593-1623. ([arXiv](http://arxiv.org/abs/1111.4246))
"""
Base.@kwdef struct ClassicNoUTurn{F<:AbstractFloat} <: DynamicTerminationCriterion
    max_depth::Int=10
    Δ_max::F=1000.0
end

"""
$(TYPEDEF)
Generalised No-U-Turn criterion as described in Section A.4.2 in [1].

# Fields
$(TYPEDFIELDS)

# References
1. Betancourt, M. (2017). A Conceptual Introduction to Hamiltonian Monte Carlo. [arXiv preprint arXiv:1701.02434](https://arxiv.org/abs/1701.02434).
"""
Base.@kwdef struct GeneralisedNoUTurn{F<:AbstractFloat} <: DynamicTerminationCriterion
    max_depth::Int=10
    Δ_max::F=1000.0
end

"""
$(TYPEDEF)
Generalised No-U-Turn criterion as described in Section A.4.2 in [1] with
added U-turn check as described in [2].

# Fields
$(TYPEDFIELDS)

# References
1. Betancourt, M. (2017). A Conceptual Introduction to Hamiltonian Monte Carlo. [arXiv preprint arXiv:1701.02434](https://arxiv.org/abs/1701.02434).
2. [https://github.com/stan-dev/stan/pull/2800](https://github.com/stan-dev/stan/pull/2800)
"""
Base.@kwdef struct StrictGeneralisedNoUTurn{F<:AbstractFloat} <: DynamicTerminationCriterion
    max_depth::Int=10
    Δ_max::F=1000.0
end
GeneralisedNoUTurn(tc::StrictGeneralisedNoUTurn) = GeneralisedNoUTurn(tc.max_depth, tc.Δ_max)

struct TurnStatistic{T}
    "Integral or sum of momenta along the integration path."
    rho::T
end
TurnStatistic() = TurnStatistic(undef)

TurnStatistic(::ClassicNoUTurn, ::PhasePoint) = TurnStatistic()
TurnStatistic(::Union{GeneralisedNoUTurn, StrictGeneralisedNoUTurn}, z::PhasePoint) =
    TurnStatistic(z.r)

combine(ts::TurnStatistic{T}, ::TurnStatistic{T}) where {T<:UndefInitializer} = ts
combine(tsl::T, tsr::T) where {T<:TurnStatistic} = TurnStatistic(tsl.rho + tsr.rho)

###
### The doubling tree algorithm for expanding trajectory.
###

"""
    Termination

Termination reasons
- `dynamic`: due to stoping criteria
- `numerical`: due to large energy deviation from starting (possibly numerical errors)
"""
struct Termination
    dynamic::Bool
    numerical::Bool
end

Base.show(io::IO, d::Termination) = print(io, "Termination(dynamic=$(d.dynamic), numerical=$(d.numerical))")
Base.:*(d1::Termination, d2::Termination) = Termination(d1.dynamic || d2.dynamic, d1.numerical || d2.numerical)
isterminated(d::Termination) = d.dynamic || d.numerical

"""
$(SIGNATURES)

Check termination of a Hamiltonian trajectory.
"""
function Termination(s::SliceTS, nt::Trajectory, H0::F, H′::F) where {F<:AbstractFloat}
    return Termination(false, !(s.ℓu < nt.termination_criterion.Δ_max + -H′))
end

"""
$(SIGNATURES)

Check termination of a Hamiltonian trajectory.
"""
function Termination(s::MultinomialTS, nt::Trajectory, H0::F, H′::F) where {F<:AbstractFloat}
    return Termination(false, !(-H0 < nt.termination_criterion.Δ_max + -H′))
end

"""
A full binary tree trajectory with only necessary leaves and information stored.
"""
struct BinaryTree
    zleft   # left most leaf node
    zright  # right most leaf node
    ts      # turn statistics
    sum_α   # MH stats, i.e. sum of MH accept prob for all leapfrog steps
    nα      # total # of leap frog steps, i.e. phase points in a trajectory
    ΔH_max  # energy in tree with largest absolute different from initial energy
end

"""
    maxabs(a, b)

Return the value with the largest absolute value.
"""
maxabs(a, b) = abs(a) > abs(b) ? a : b

"""
$(SIGNATURES)
Merge a left tree `treeleft` and a right tree `treeright` under given Hamiltonian `h`,
then draw a new candidate sample and update related statistics for the resulting tree.
"""
function combine(treeleft::BinaryTree, treeright::BinaryTree)
    return BinaryTree(
        treeleft.zleft,
        treeright.zright,
        combine(treeleft.ts, treeright.ts),
        treeleft.sum_α + treeright.sum_α,
        treeleft.nα + treeright.nα,
        maxabs(treeleft.ΔH_max, treeright.ΔH_max),
    )
end

"""
$(SIGNATURES)
Detect U turn for two phase points (`zleft` and `zright`) under given Hamiltonian `h`
using the (original) no-U-turn cirterion.

Ref: https://arxiv.org/abs/1111.4246, https://arxiv.org/abs/1701.02434
"""
function isterminated(::ClassicNoUTurn, h::Hamiltonian, t::BinaryTree)
    # z0 is starting point and z1 is ending point
    z0, z1 = t.zleft, t.zright
    Δθ = z1.θ - z0.θ
    s = (dot(Δθ, ∂H∂r(h, -z0.r)) >= 0) || (dot(-Δθ, ∂H∂r(h, z1.r)) >= 0)
    return Termination(s, false)
end

"""
$(SIGNATURES)
Detect U turn for two phase points (`zleft` and `zright`) under given Hamiltonian `h`
using the generalised no-U-turn criterion.

Ref: https://arxiv.org/abs/1701.02434
"""
function isterminated(::GeneralisedNoUTurn, h::Hamiltonian, t::BinaryTree)
    rho = t.ts.rho
    s = generalised_uturn_criterion(rho, ∂H∂r(h, t.zleft.r), ∂H∂r(h, t.zright.r))
    return Termination(s, false)
end

"""
$(SIGNATURES)
Detect U turn for two phase points (`zleft` and `zright`) under given Hamiltonian `h`
using the generalised no-U-turn criterion with additional U-turn checks.

Ref: https://arxiv.org/abs/1701.02434 https://github.com/stan-dev/stan/pull/2800
"""
function isterminated(
    tc::StrictGeneralisedNoUTurn, h::Hamiltonian, t, tleft, tright
)
    # (Non-strict) generalised U-turn check
    s1 = isterminated(GeneralisedNoUTurn(tc), h, t)

    # U-turn checks for left and right subtree
    # See https://discourse.mc-stan.org/t/nuts-misses-u-turns-runs-in-circles-until-max-treedepth/9727/33
    # for a visualisation.
    s2 = check_left_subtree(h, t, tleft, tright)
    s3 = check_right_subtree(h, t, tleft, tright)
    return s1 * s2 * s3
end

"""
$(SIGNATURES)
Do a U-turn check between the leftmost phase point of `t` and the leftmost
phase point of `tright`, the right subtree.
"""
function check_left_subtree(
    h::Hamiltonian, t::T, tleft::T, tright::T
) where {T<:BinaryTree}
    rho = tleft.ts.rho + tright.zleft.r
    s = generalised_uturn_criterion(rho, ∂H∂r(h, t.zleft.r), ∂H∂r(h, tright.zleft.r))
    return Termination(s, false)
end

"""
$(SIGNATURES)
Do a U-turn check between the rightmost phase point of `t` and the rightmost
phase point of `tleft`, the left subtree.
"""
function check_right_subtree(
    h::Hamiltonian, t::T, tleft::T, tright::T
) where {T<:BinaryTree}
    rho = tleft.zright.r + tright.ts.rho
    s = generalised_uturn_criterion(rho, ∂H∂r(h, tleft.zright.r), ∂H∂r(h, t.zright.r))
    return Termination(s, false)
end

function generalised_uturn_criterion(rho, p_sharp_minus, p_sharp_plus)
    return (dot(rho, p_sharp_minus) <= 0) || (dot(rho, p_sharp_plus) <= 0)
end

function isterminated(tc::TC, h::Hamiltonian, t::BinaryTree, _tleft, _tright) where {TC<:Union{ClassicNoUTurn, GeneralisedNoUTurn}}
    return isterminated(tc, h, t)
end

"Recursivly build a tree for a given depth `j`."
function build_tree(
    rng::AbstractRNG,
    nt::Trajectory{TS, I, TC},
    h::Hamiltonian,
    z::PhasePoint,
    sampler::AbstractTrajectorySampler,
    v::Int,
    j::Int,
    H0::AbstractFloat,
) where {TS<:AbstractTrajectorySampler, I<:AbstractIntegrator, TC<:DynamicTerminationCriterion}
    if j == 0
        # Base case - take one leapfrog step in the direction v.
        z′ = step(nt.integrator, h, z, v)
        H′ = energy(z′)
        ΔH = H′ - H0
        α′ = exp(min(0, -ΔH))
        sampler′ = TS(sampler, H0, z′)
        return BinaryTree(z′, z′, TurnStatistic(nt.termination_criterion, z′), α′, 1, ΔH), sampler′, Termination(sampler′, nt, H0, H′)
    else
        # Recursion - build the left and right subtrees.
        tree′, sampler′, termination′ = build_tree(rng, nt, h, z, sampler, v, j - 1, H0)
        # Expand tree if not terminated
        if !isterminated(termination′)
            # Expand left
            if v == -1
                tree′′, sampler′′, termination′′ = build_tree(rng, nt, h, tree′.zleft, sampler, v, j - 1, H0) # left tree
                treeleft, treeright = tree′′, tree′
            # Expand right
            else
                tree′′, sampler′′, termination′′ = build_tree(rng, nt, h, tree′.zright, sampler, v, j - 1, H0) # right tree
                treeleft, treeright = tree′, tree′′
            end
            tree′ = combine(treeleft, treeright)
            sampler′ = combine(rng, sampler′, sampler′′)
            termination′ = termination′ * termination′′ * isterminated(nt.termination_criterion, h, tree′, treeleft, treeright)
        end
        return tree′, sampler′, termination′
    end
end

function transition(
    rng::AbstractRNG,
    τ::Trajectory{TS, I, TC},
    h::Hamiltonian,
    z0::PhasePoint,
) where {TS<:AbstractTrajectorySampler, I<:AbstractIntegrator, TC<:DynamicTerminationCriterion}
    H0 = energy(z0)
    tree = BinaryTree(z0, z0, TurnStatistic(τ.termination_criterion, z0), zero(H0), zero(Int), zero(H0))
    sampler = TS(rng, z0)
    termination = Termination(false, false)
    zcand = z0

    j = 0
    while !isterminated(termination) && j < τ.termination_criterion.max_depth
        # Sample a direction; `-1` means left and `1` means right
        v = rand(rng, [-1, 1])
        if v == -1
            # Create a tree with depth `j` on the left
            tree′, sampler′, termination′ = build_tree(rng, τ, h, tree.zleft, sampler, v, j, H0)
            treeleft, treeright = tree′, tree
        else
            # Create a tree with depth `j` on the right
            tree′, sampler′, termination′ = build_tree(rng, τ, h, tree.zright, sampler, v, j, H0)
            treeleft, treeright = tree, tree′
        end
        # Perform a MH step and increse depth if not terminated
        if !isterminated(termination′)
            j = j + 1   # increment tree depth
            if mh_accept(rng, sampler, sampler′)
                zcand = sampler′.zcand
            end
        end
        # Combine the proposed tree and the current tree (no matter terminated or not)
        tree = combine(treeleft, treeright)
        # Update sampler
        sampler = combine(zcand, sampler, sampler′)
        # update termination
        termination = termination * termination′ * isterminated(τ.termination_criterion, h, tree, treeleft, treeright)
    end

    H = energy(zcand)
    tstat = merge(
        (
            n_steps=tree.nα,
            is_accept=true,
            acceptance_rate=tree.sum_α / tree.nα,
            log_density=zcand.ℓπ.value,
            hamiltonian_energy=H,
            hamiltonian_energy_error=H - H0,
            max_hamiltonian_energy_error=tree.ΔH_max,
            tree_depth=j,
            numerical_error=termination.numerical,
        ),
        stat(τ.integrator),
    )

    return Transition(zcand, tstat)
end

###
### Initialisation of step size
###

"""
A single Hamiltonian integration step.

NOTE: this function is intended to be used in `find_good_stepsize` only.
"""
function A(h, z, ϵ)
    z′ = step(Leapfrog(ϵ), h, z)
    H′ = energy(z′)
    return z′, H′
end

"Find a good initial leap-frog step-size via heuristic search."
function find_good_stepsize(
    rng::AbstractRNG,
    h::Hamiltonian,
    θ::AbstractVector{T};
    max_n_iters::Int=100,
) where {T<:Real}
    # Initialize searching parameters
    ϵ′ = ϵ = T(0.1)
    a_min, a_cross, a_max = T(0.25), T(0.5), T(0.75) # minimal, crossing, maximal accept ratio
    d = T(2.0)
    # Create starting phase point
    r = rand(rng, h.metric) # sample momentum variable
    z = phasepoint(h, θ, r)
    H = energy(z)

    # Make a proposal phase point to decide direction
    z′, H′ = A(h, z, ϵ)
    ΔH = H - H′ # compute the energy difference; `exp(ΔH)` is the MH accept ratio
    direction = ΔH > log(a_cross) ? 1 : -1

    # Crossing step: increase/decrease ϵ until accept ratio cross a_cross.
    for _ = 1:max_n_iters
        # `direction` being  `1` means MH ratio too high
        #     - this means our step size is too small, thus we increase
        # `direction` being `-1` means MH ratio too small
        #     - this means our step szie is too large, thus we decrease
        ϵ′ = direction == 1 ? d * ϵ : 1 / d * ϵ
        z′, H′ = A(h, z, ϵ)
        ΔH = H - H′
        DEBUG && @debug "Crossing step" direction H′ ϵ "α = $(min(1, exp(ΔH)))"
        if (direction == 1) && !(ΔH > log(a_cross))
            break
        elseif (direction == -1) && !(ΔH < log(a_cross))
            break
        else
            ϵ = ϵ′
        end
    end
    # Note after the for loop,
    # `ϵ` and `ϵ′` are the two neighbour step sizes across `a_cross`.

    # Bisection step: ensure final accept ratio: a_min < a < a_max.
    # See https://en.wikipedia.org/wiki/Bisection_method

    ϵ, ϵ′ = ϵ < ϵ′ ? (ϵ, ϵ′) : (ϵ′, ϵ)  # ensure ϵ < ϵ′;
    # Here we want to use a value between these two given the
    # criteria that this value also gives us a MH ratio between `a_min` and `a_max`.
    # This condition is quite mild and only intended to avoid cases where
    # the middle value of `ϵ` and `ϵ′` is too extreme.
    for _ = 1:max_n_iters
        ϵ_mid = middle(ϵ, ϵ′)
        z′, H′ = A(h, z, ϵ_mid)
        ΔH = H - H′
        DEBUG && @debug "Bisection step" H′ ϵ_mid "α = $(min(1, exp(ΔH)))"
        if (exp(ΔH) > a_max)
            ϵ = ϵ_mid
        elseif (exp(ΔH) < a_min)
            ϵ′ = ϵ_mid
        else
            ϵ = ϵ_mid
            break
        end
    end

    return ϵ
end

function find_good_stepsize(
    h::Hamiltonian,
    θ::AbstractVector{<:AbstractFloat};
    max_n_iters::Int=100,
)
    return find_good_stepsize(GLOBAL_RNG, h, θ; max_n_iters=max_n_iters)
end

"Perform MH acceptance based on energy, i.e. negative log probability."
function mh_accept_ratio(
    rng::AbstractRNG,
    Horiginal::T,
    Hproposal::T,
) where {T<:AbstractFloat}
    α = min(one(T), exp(Horiginal - Hproposal))
    accept = rand(rng, T) < α
    return accept, α
end

function mh_accept_ratio(
    rng::Union{AbstractRNG, AbstractVector{<:AbstractRNG}},
    Horiginal::AbstractVector{<:T},
    Hproposal::AbstractVector{<:T},
) where {T<:AbstractFloat}
    α = min.(one(T), exp.(Horiginal .- Hproposal))
    # NOTE: There is a chance that sharing the RNG over multiple
    #       chains for accepting / rejecting might couple
    #       the chains. We need to revisit this more rigirously
    #       in the future. See discussions at
    #       https://github.com/TuringLang/AdvancedHMC.jl/pull/166#pullrequestreview-367216534
    accept = rand(rng, T, length(Horiginal)) .< α
    return accept, α
end


text:
using NES

env = SMBEnv("smb", :RIGHT_ONLY)
current_state = reset!(env)
done = false
actions = [env.action_map |> keys |> rand for _=1:400]

for action ∈ actions
    global current_state, done
    current_state, reward, done, info = step!(env, action)
    sleep(0.001)
    render(env) |> display
    done && break
end


text:
# @everywhere module_dir = "/Users/brandonbehring/Desktop/Leap_Frog_2019"
# @everywhere push!(LOAD_PATH, "/Users/brandonbehring/Desktop/Leap_Frog_2019")
# # @everywhere module_dir ="/Users/brandonbehring/Desktop/Leap_Frog_2019"
# # @everywhere push!(LOAD_PATH, $module_dir)
# # @everywhere thisDir = dirname(@__FILE__())
# # @everywhere any(path -> path == thisDir, LOAD_PATH) || push!(LOAD_PATH, thisDir)
# include("escape.jl")
# @everywhere using .escape

using ProgressMeter
using Printf
using MATLAB
using BSON: @save, @load
include("escape_num_dimer.jl")

@everywhere begin
    include("escape_num_dimer.jl")
    # push!(LOAD_PATH, "/Users/brandonbehring/Desktop/Leap_Frog_2019")
    push!(LOAD_PATH, "/home/brandon_behring/Desktop/Leap_Frog_2019")
    using ProgressMeter
    using Printf
    using MATLAB
    using .escape_num_dimer
    t_end = exp(5)
    # width = sqrt(3); height = 1
    width = 2.5; height=5
    n_iter_Q = 100;n_iter_P = 2*n_iter_Q; N = n_iter_P * n_iter_Q;
    ArrP = range(-width, stop = width, length = n_iter_P)
    ArrQ = range(-height, stop = height, length = n_iter_Q)
    mesh_P=[P for Q in ArrQ, P in ArrP]
    mesh_Q=[Q for Q in ArrQ, P in ArrP]
    mesh = [(P, Q) for Q in ArrQ, P in ArrP]
    mesh_list = reshape(mesh, 1, :)
    t_end = t_end * ones(n_iter_Q,n_iter_P);
    H=.2
    location = "/mnt/bdd38f66-9ece-451a-b915-952523c139d2/Escape/"
    # location = "/Users/brandonbehring/Desktop/"
end

    @everywhere Energy = H* ones(n_iter_Q,n_iter_P);
    # @everywhere Energy = H* ones(N);

    @everywhere h = replace(@sprintf("%.8f",H), "." => "_")
    @everywhere nQ = @sprintf("%d",n_iter_Q)
    @everywhere file_name = location * "Dimer_time_before_Escape_" * h *nQ* ".fig"
    println(file_name)
    h_BSON="escape_data"*h*nQ*".bson"
    num_until_exit = @showprogress pmap(escape_num_dimer.escape_exit_num_dimer, mesh, t_end, Energy)

    logz=log1p.(num_until_exit)

    SAVE_DATA=Dict("escape_times"=>num_until_exit, "mesh"=>mesh, "ArrP"=>ArrP, "ArrQ"=>ArrQ, "Energy"=>H)
    @save h_BSON SAVE_DATA

    mat"figure();set(gcf, 'Position',  [0, 0, 1500, 1500]); hold on;"
    # # # # mat"plot($P_0,$Q_0 ,'b.','MarkerSize',3)"
    # # # # mat"plot($P_1,$Q_1 ,'r.','MarkerSize',3)"
    # # # # mat"plot($P_2,$Q_2 ,'g.','MarkerSize',3)"
    # # # # mat"plot($P_4,$Q_4 ,'c.','MarkerSize',3)"
    # # # # mat"plot($P_5,$Q_5 ,'y.','MarkerSize',3)"
    mat"imagesc([-$width,$width],[-$height,$height ],$logz)"
    mat"colorbar"

    mat"axis([ -$width,$width,-$height,$height ])"
    # # # mat"axis([ -1,1,-1,1])"
    mat"savefig($file_name)"


text:
"""
A Julia interface to [nmap](https://nmap.org/) that allows scan and output manipulation

NMAP.jl aims to provide developers with idiomatic nmap bindings to aid the creation of security audit tools

## Scanner

## Scan

## run!

## Examples

**Create scanners**

Example command: `nmap localhost 127.0.0.1 scanme.nmap.org -p 21-25,80,443,8080 -sV -O -sC --packet-trace -T0`
```julia
scanner = NMAP.Scanner(
    NMAP.targets("localhost", "127.0.0.1", "scanme.nmap.org"),
    NMAP.ports("21-25", "80", "443", "8080"),
    NMAP.service_info(),
    NMAP.os_detection(),
    NMAP.default_script(),
    NMAP.packet_trace(),
    NMAP.timingtemplate(NMAP.paranoid))
```

**Execute a scan**
```
scan = NMAP.run!(scanner)
```
"""
module NMAP

    export Scan
    export Scanner

    using Dates
    using EzXML
    using JSON
    using Lazy: Lazy, @forward
    using Printf
    using Revise
    using StructTypes
    using XMLDict

    include("unmarshalling.jl")
    using .Marshalling: Marshalling, Marsh, Leaf, unmarshall
    include("scan.jl")
    include("scanner.jl")
    include("options.jl")
    #include("methods.jl")

end


text:
# Copyright 2014 Mitchell Kember. Subject to the MIT License.
# Google Code Jam Qualification Round 2014
# Problem D: Deceitful War

# Returns the points Naomi can score if she plays War optimally.
function points_war(msn, msk)
	points = 0
	mskc = copy(msk)
	for mn in msn
		i = findfirst(mk -> mk > mn, mskc)
		if i == 0
			shift!(mskc)
			points += 1
		else
			deleteat!(mskc, i)
		end
	end
	points
end

# Returns the points Naomi can score if she plays Deceitful War optimally.
# Assumes that msn and msk are sorted.
# Empties both msn and msk.
function points_deceitfulwar!(msn, msk)
	points = 0
	while !isempty(msn)
		i = findfirst(mn -> mn > msk[end], msn)
		if i == 0
			shift!(msn)
			pop!(msk)
		else
			deleteat!(msn, i)
			pop!(msk)
			points += 1
		end
	end
	points
end

function main()
	n = int(readline())
	for i = 1:n
		int(readline()) # number of masses
		msn = sort(float(split(readline())))
		msk = sort(float(split(readline())))
		pw = points_war(msn, msk)
		pdw = points_deceitfulwar!(msn, msk)
		println("Case #$i: $pdw $pw")
	end
end

main()


text:
export log_best, log_all
using Logging
using Printf
using Dates

function log_best(id::String, seed::Int64, eval_count::Int64, max_fit::Float64,
                  best::Chromosome, ea::Function, ctype::DataType, log_gen::Bool)
    if log_gen
        logstr = @sprintf("E: %s %d %d %0.5f %d %d",
                          id, seed, eval_count, max_fit,
                          sum([n.active for n in best.nodes]),
                          length(best.nodes))
        if Config.log_config
            logstr = string(logstr,
                            @sprintf("%s %s %s",
                                     string(ea), string(ctype), Config.to_string()))
        end
        @info logstr timestamp=Dates.now(UTC)
    end
    if Config.save_best
        @info(@sprintf("C: %s", string(best.genes)))
    end
    flush(current_logger().stream)
    nothing
end

function log_all(id::String, seed::Int64, eval_count::Int64, max_fit::Float64,
                  best::Chromosome, ea::Function, ctype::DataType, log_gen::Bool)
    log_best(id, seed, eval_count, max_fit, best, ea, ctype, true)
end


text:
module AST
using MLStyle
using MLStyle.AbstractPatterns
using MLStyle.Err
export @matchast, @capture, Capture

@nospecialize
struct Capture end

function MLStyle.pattern_uncall(
    ::Type{Capture},
    self::Function,
    type_params::AbstractArray,
    type_args::AbstractArray,
    args::AbstractArray,
)
    isempty(type_params) || error("A Capture requires no type params.")
    isempty(type_args) || error("A Capture pattern requires no type arguments.")
    length(args) === 1 || error("A Capture pattern accepts one argument.")
    function extract(::Any, ::Int, scope::ChainDict{Symbol, Symbol}, ::Any)
        ret = Expr(:call, Dict)
        for_chaindict(scope) do k, v
            k = QuoteNode(k)
            push!(ret.args, :($k => $v))
        end
        ret
    end
    decons(extract, [self(args[1])])
end

function matchast(target, actions, source::LineNumberNode, mod::Module)
    @switch actions begin
        @case Expr(:quote, Expr(:block, stmts...))
        last_lnode = source
        cbl = Expr(:block)
        for stmt in stmts
            @switch stmt begin
                @case ::LineNumberNode
                last_lnode = stmt
                continue
                @case :($a => $b)
                push!(cbl.args, last_lnode, :($(Expr(:quote, a)) => $b))
                continue
                @case _
                throw(
                    SyntaxError(
                        "Malformed ast template, should be formed as `a => b`, at $(string(last_lnode)).",
                    ),
                )
            end
        end
        return gen_match(target, cbl, source, mod)
        @case _
        msg = "Malformed ast template, the second arg should be a block with a series of pairs(`a => b`), at $(string(source))."
        throw(SyntaxError(msg))
    end
end

"""
An eye candy of `@match` for AST matching.

e.g.,

```julia
    @matchast :(1 + 1) quote
        \$a + 1 => a
    end # 1

    @matchast :(f(a, b)) quote
        \$(Expr(:call, :f, :a, :b)) =>
         dosomething
    end
```
"""
macro matchast(template, actions)
    res = matchast(template, actions, __source__, __module__)
    res = init_cfg(res)
    return esc(res)
end

"""
@capture template
@capture template expr

Template matching for expressions.

```julia
julia> @capture f(\$x)  :(f(2))
Dict{Symbol,Int64} with 1 entry:
  :x => 2
```

If the template doesn't match input AST, return `nothing`.
"""
:(@capture)

macro capture(template)
    farg = gensym("expression")
    fbody = gen_capture(template, farg, __source__, __module__) |> init_cfg
    fhead = Expr(:call, farg, farg)
    esc(Expr(:function, fhead, fbody))
end

macro capture(template, ex)
    res = gen_capture(template, ex, __source__, __module__)
    res = init_cfg(res)
    return esc(res)
end

function gen_capture(template::Any, ex::Any, source::LineNumberNode, mod::Module)
    template = Expr(:quote, template)
    sym = :__SCOPE_CAPTURE__
    p_capture_scope = Expr(:call, Capture, sym)
    p_whole = Expr(:&&, template, p_capture_scope)
    tbl = Expr(:block, :($p_whole => $sym))

    gen_match(ex, tbl, source, mod)
end
@specialize
end


text:
import RCall

try
    RCall.R"""
        dir.create(path = Sys.getenv("R_LIBS_USER"), showWarnings = FALSE, recursive = TRUE)
        .libPaths(Sys.getenv("R_LIBS_USER"))
        if ( !require("DMwR") ) {
            install.packages(
                "DMwR",
                lib = Sys.getenv("R_LIBS_USER"),
                repos = "http://cran.r-project.org/",
                type="binary",
                )
        }
    """
catch e
    @warn(string("ignoring error: "), e)
    RCall.R"""
        dir.create(path = Sys.getenv("R_LIBS_USER"), showWarnings = FALSE, recursive = TRUE)
        .libPaths(Sys.getenv("R_LIBS_USER"))
        if ( !require("DMwR") ) {
            install.packages(
                "DMwR",
                lib = Sys.getenv("R_LIBS_USER"),
                repos = "http://cran.r-project.org/",
                type="source",
                )
        }
    """
end

result_rcall_library_DMwR = RCall.R"""
    library("DMwR")
"""

@info(
    string("Result of calling library(\"DMwR\"): "),
    result_rcall_library_DMwR,
    )


text:
using MRC, TranscodingStreams
using Test

@testset "MRC.jl" begin
    include("utils.jl")
    include("header.jl")
    include("io.jl")
end


text:
using DrWatson
@quickactivate :AlbedoProperties

include(srcdir("plotting", "recipes.jl"))

# width correspoding to half-column width
figx = 12
figy = 8
PyPlot.rc("figure", figsize = (figx, figy))
PyPlot.rc("legend", fontsize = 24, framealpha = 0.8, handlelength = 1)
PyPlot.rc("grid", color = "0.75", alpha = 0.5, lw = 1.5)

max_stars_repo_name: Datseris/EarthAlbedoSymmetry

max_issues_repo_name: Datseris/EarthAlbedoSymmetry

max_forks_repo_name: Datseris/EarthAlbedoSymmetry


text:
using MultilaneDriving
using Base.Test
using POMDPToolbox

using StaticArrays
using POMDPs

using Gallium
using Traceur

rng = MersenneTwister(7)

dt = 0.1

driver = NoCrashMLDriver(dt)

car = HighwayVehicle()
carstate = SVector(0.0, 0.0, 20.0)
carstate2 = SVector(100.0, 0.0, 20.0)

road = Highway(3, 3.0, 100.0)

context = HighwayDrivingContext(car, [carstate, carstate2], road)

# @trace action_distribution(driver, MLDriverState(context, nothing, context), 1, AccelSpeed)
# @code_warntype action(driver.long, context, 1, LongAccel)
ad = @inferred action_distribution(driver, MLDriverState(context, nothing, context), 1, AccelSpeed)

a = rand(MersenneTwister(4), ad)

propagate(car, carstate, a, context, dt)

s = MLDState([carstate, SVector(50.0, 0.0, 20.0)], [driver], [MLDriverState(context, nothing, context)])
p = MLDPOMDP()
sp, r = @inferred generate_sr(p, s, AccelSpeed(1.0, 0.1), MersenneTwister(5))

as = @inferred actions(p, sp)
a = @inferred rand(rng, as)

# @show [rand(rng, as) for i in 1:10]
# hr = HistoryRecorder(max_steps=10, rng=MersenneTwister(4))
# simulate(hr, p, RandomPolicy(p, rng=MersenneTwister(5)))


text:
#!/usr/bin/env julia

using SpinQubits, ArgParse

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table! s begin
        "--length", "-L"
            arg_type = Int
	        default = 4
        "--singlet", "-s"
            arg_type = Bool
            default = false
        "--sigmas", "-σ"
		    arg_type = String
	        default = "0.0,0.0,0.0"
        "--dirspec"
	        arg_type = String
	        default = ""
	    "--nReals"
	        arg_type = Int
	        default = 10000
    end

    return parse_args(s)
end

args = parse_commandline()
sigmas = args["sigmas"]
sigmas = replace(sigmas, "[" => "")
sigmas = replace(sigmas, "]" => "")
sigArray = parse.(Float64,split(sigmas,","))
disGam = sigArray[2]
spacing = maximum(sigArray) == 0.0 ? 0.01 : 0.03
realSigArray = [sigArray[1], 0.0, sigArray[3]]

calculateFidelities(4, 0.01, 0.0, zeros(3), 10, 0.01; verbose=false) # Compilation run

#saveFidelities( 4, 0.01, 0.1, zeros(3), 10, 0.01; singlet=false)

saveFidelities(args["length"], 0.01, disGam, realSigArray, args["nReals"], spacing ;index=args["dirspec"], singlet=args["singlet"])


text:
"""
	$(SIGNATURES)

Tranform parameter into bounds specified by `LinearTransformation`. Calibrated values only.
Out of bounds values are pushed inside bounds with a warning.
"""
function transform_param(tr :: LinearTransformation{F1}, p :: AbstractParam) where {F1, F2}
    pValue = clamp.(calibrated_value(p), calibrated_lb(p), calibrated_ub(p));
    # if any(value(p) .> p.ub)  ||  any(value(p) .< p.lb)
    #     @warn "Values out of bounds for $(p.name)."
    #     value = min.(p.ub, max.(p.lb, value(p)));
    # else
    #     value = value(p);
    # end
    return lb(tr) .+ (ub(tr) .- lb(tr)) .* (pValue .- calibrated_lb(p)) ./ (calibrated_ub(p) .- calibrated_lb(p));
end

"""
	$(SIGNATURES)

Undo parameter transformation. Calibrated values only.

`value` must be the same size as `value(p)`.
"""
function untransform_param(tr :: LinearTransformation{F1},
    p :: AbstractParam, pValue) where F1

    @assert size(pValue) == size(calibrated_value(p))  "Size mismatch: $(size(pValue)) vs $(size(calibrated_value(p))) for $p";

    if any(pValue .> ub(tr))  ||  any(pValue .< lb(tr))
        @warn """
            Values for $p out of bounds:
            $pValue
            """;
    end
    # @assert all(pValue .<= ub(tr))  "Values to high: $pValue  vs  $(ub(tr))"
    # @assert all(pValue .>= lb(tr))

    outV = clamp.(
        calibrated_lb(p) .+ (calibrated_ub(p) .- calibrated_lb(p)) .* (pValue .- lb(tr)) ./ (ub(tr) .- lb(tr)),
        calibrated_lb(p), calibrated_ub(p));

    @assert size(outV) == size(pValue);
    return outV
end

# function untransform(tr :: LinearTransformation{F1}) where F1

"""
	$(SIGNATURES)

Bounds on transformed parameters.
"""
function transform_bounds(tr :: LinearTransformation{F1}) where F1
    return lb(tr), ub(tr)
end

# -------


text:
# should escape column-names:
# http://stackoverflow.com/questions/2901453/sql-standard-to-escape-column-names
# https://www.sqlite.org/lang_keywords.html
ident(identifier::Symbol) = string(identifier)

function ident(identifier::Expr)
    @assert identifier.head == :.
    @assert length(identifier.args) == 2
    table,colname = identifier.args
    @assert isa(table, Symbol)
    @assert isa(colname, Expr)
    @assert colname.head == :quote
    @assert length(colname.args) == 1
    @assert isa(colname.args[1], Symbol)
    string(identifier)
end

# should process this properly in the future
_sqlexpr(ex::QueryArg) = "$ex"

translatesql(q::Symbol, offset::Int) = ident(q) # table-name
_translatesubquery(q::Symbol, offset::Int) = ident(q)
_translatesubquery(q::QueryNode, offset::Int) = "($(translatesql(q, offset)))"

_selectarg(a::Symbol) = string(a) # assume it corresponds to a column-name

function _selectarg(a::Expr)
    if a.head == :kw # newcol=col (SELECT col AS newcol)
        @assert length(a.args) == 2
        newcol,expr = a.args
        @assert isa(newcol, Symbol)
        return "$(_sqlexpr(expr)) AS $newcol"
    else
        @assert a.head == :. # table.columnname
        return ident(a)
    end
end

"Returns `true` is the last GROUP BY argument a `having(...)` expression"
_groupbyhaving(arg::QueryArg) =
    isa(arg, Expr) && arg.head == :call && arg.args[1] == :having

_orderbyterm(term::Symbol) = ident(term)
function _orderbyterm(term::QueryArg)
    if isa(term, Expr) && term.head == :call
        @assert length(term.args) == 2 "invalid orderby term: $term"
        order, identifier = term.args
        order == :asc && return "$(ident(identifier)) ASC"
        order == :desc && return "$(ident(identifier)) DESC"
    else
        return ident(term)
    end
    error("Unable to parse orderby term: $term")
end

function _parsejoinargs(q::JoinNode, offset::Int)
    nargs = length(q.args)
    @assert nargs > 0
    @assert isa(q.args[1], Symbol) "only support joining to a table-alias"
    indent = " " ^ offset
    on = by = ""
    source = string(q.args[1])
    if length(q.args) > 1 # has join constraints (e.g. join(t2, on(expr)))
        @assert length(q.args) == 2 "only allow for one constraint"
        constraint = q.args[2]
        @assert isa(constraint, Expr)
        @assert constraint.head == :call
        cons = exf(constraint); args = exfargs(constraint)
        if cons == :on
            # TODO: parse the arguments here
            on = "\n$(indent)ON $(join(args, "\n $indent  AND "))"
        else
            @assert cons == :by
            columnnames = join(map(ident,args),",\n       $indent")
            by = "\n$(indent)USING ($columnnames)"
        end
    end
    source, on, by
end

function _parsejoin(q::LeftJoinNode, offset::Int=0)
    source, on, by = _parsejoinargs(q, offset+17)
    "LEFT JOIN", source, on, by
end

function _parsejoin(q::OuterJoinNode, offset::Int=0)
    source, on, by = _parsejoinargs(q, offset+23)
    "LEFT OUTER JOIN", source, on, by
end

function _parsejoin(q::InnerJoinNode, offset::Int=0)
    source, on, by = _parsejoinargs(q, offset+18)
    "INNER JOIN", source, on, by
end

function _parsejoin(q::CrossJoinNode, offset::Int=0)
    source, on, by = _parsejoinargs(q, offset+18)
    "CROSS JOIN", source, on, by
end

function translatesql(q::FilterNode, offset::Int=0)
    @assert length(q.args) > 0 "you shouldn't filter by nothing"
    indent = " " ^ offset
    source = _translatesubquery(q.input, offset+8)
    # TODO: should properly parse q.args
    conditions = join(q.args, "\n$(indent)   AND ")
    "SELECT *\n $indent FROM $source\n$indent WHERE $conditions"
end

function translatesql(q::SelectNode, offset::Int=0)
    @assert length(q.args) > 0 "you shouldn't select by nothing"
    indent = " " ^ offset
    source = _translatesubquery(q.input, offset+8)
    resultcolumns = join(map(_selectarg, q.args), ",\n $indent      ")
    "SELECT $resultcolumns\n $indent FROM $source"
end

function translatesql(q::DistinctNode, offset::Int=0)
    @assert length(q.args) > 0 "you shouldn't select by nothing"
    indent = " " ^ offset; pad = " " ^ 9
    source = _translatesubquery(q.input, offset+8)
    resultcolumns = join(map(_selectarg, q.args), ",\n $indent$pad      ")
    "SELECT DISTINCT $resultcolumns\n $indent FROM $source"
end

function translatesql(q::GroupbyNode, offset::Int=0)
    @assert length(q.args) > 0 "you shouldn't groupby nothing"
    indent = " " ^ offset
    source = _translatesubquery(q.input, offset+10)
    lastarg = q.args[end]
    groupby = if _groupbyhaving(lastarg)
        groupbyargs = q.args[1:end-1]
        @assert length(groupbyargs) > 0
        havingargs = exfargs(lastarg)
        @assert length(havingargs) > 0
        # TODO: should properly parse havingargs
        conditions = "HAVING $(join(havingargs, "\n $indent    AND "))"
        "GROUP BY $(join(map(ident, groupbyargs), ", "))\n $conditions"
    else
        "GROUP BY $(join(map(ident, q.args), ", "))"
    end
    "  SELECT *\n   $indent FROM $source\n$(indent)$groupby"
end

function translatesql(q::OrderbyNode, offset::Int=0)
    @assert length(q.args) > 0 "you shouldn't order by nothing"
    indent = " " ^ offset
    source = _translatesubquery(q.input, offset+10)
    orderby = join(map(_orderbyterm, q.args), ",\n $indent        ")
    "  SELECT *\n   $indent FROM $source\n$(indent)ORDER BY $orderby"
end

function translatesql(q::LimitNode, offset::Int=0)
    indent = " " ^ offset
    source = _translatesubquery(q.input, offset+8)
    @assert length(q.limit) == 1
    @assert isa(q.limit[1], Integer)
    "SELECT *\n $indent FROM $source\n$(indent) LIMIT $(q.limit[1])"
end

function translatesql(q::OffsetNode, offset::Int=0)
    indent = " " ^ offset
    source = _translatesubquery(q.input, offset+8)
    @assert length(q.offset) == 1
    @assert isa(q.offset[1], Integer)
    limit = "LIMIT -1 OFFSET $(q.offset[1])"
    "SELECT *\n $indent FROM $source\n$(indent) $limit"
end

function translatesql(q::JoinNode, offset::Int=0)
    indent = " " ^ offset
    join, table, on, by = _parsejoin(q, offset)
    source = _translatesubquery(q.input, offset+8)
    "SELECT *\n $indent FROM $source\n       $indent$join $table$on$by"
end


text:
### Dependence

immutable Dependence{S<:Variable, T<:Variable}
  index::Integer
  source::S
  target::T
end

Dependence{S<:Variable, T<:Variable}(source::S, target::T) = Dependence(0, source, target)

typealias DependenceVector{D<:Dependence} Vector{D}

edge_index(d::Dependence) = d.index
source(d::Dependence) = e.source
target(d::Dependence) = e.target

is_indexed(d::Dependence) = d.index > 0 ? true : false

revedge{S<:Variable, T<:Variable}(d::Dependence{S, T}) = Dependence(d.index, d.target, d.source)

Base.show(io::IO, d::Dependence) = print(io, "Dependence [$(d.index)]: $(d.target.key) | $(d.source.key)")


text:
# This file was generated, do not modify it. # hide
using MLJ, Statistics, PrettyPrinting
X, y = @load_iris;


text:
using CatBoost
using Documenter

makedocs(; modules=[CatBoost], sitename="CatBoost.jl", authors="Beacon Biosignals, Inc.",
         pages=["API Documentation" => "index.md"])

deploydocs(; repo="github.com/beacon-biosignals/CatBoost.jl.git", push_preview=true,
           devbranch="main")


text:
"""
`AbstractGeoArray` wraps an array (or location of an array) and metadata
about its contents. It may be memory ([`GeoArray`](@ref)) or disk-backed
([`NCDarray`](@ref), [`GDALarray`](@ref), [`GRDarray`](@ref)).

`AbstractGeoArray`s inherit from [`AbstractDimArray`]($DDarraydocs)
from DimensionalData.jl. They can be indexed as regular Julia arrays or with
DimensionalData.jl [`Dimension`]($DDdimdocs)s. They will plot as a heatmap in
Plots.jl with correct coordinates and labels, even after slicing with
`getindex` or `view`. `getindex` on a `AbstractGeoArray` will always return
a memory-backed `GeoArray`.
"""
abstract type AbstractGeoArray{T,N,D,A} <: AbstractDimensionalArray{T,N,D,A} end

# Marker singleton for lazy loaded arrays, only used for broadcasting.
# Can be removed when DiskArrays.jl is used everywhere
struct LazyArray{T,N} <: AbstractArray{T,N} end

# Interface methods ###########################################################

"""
    missingval(x)

Returns the value representing missing data in the dataset """ function missingval end
missingval(x) = missing
missingval(A::AbstractGeoArray) = A.missingval

"""
    crs(x)

Get the crs projection of a dim or for the `Lat`/`Lon` dims of an array.
"""
function crs end
crs(A::AbstractGeoArray) =
    if hasdim(A, Lat)
        crs(dims(A, Lat))
    elseif hasdim(A, Lon)
        crs(dims(A, Lon))
    else
        error("No Lat or Lon dimension, crs not available")
    end
crs(dim::Dimension) = crs(mode(dim), dim)

"""
    usercrs(x)

Get the user facing crs projection of a [`Projected`](@ref) mode
dim or for the `Lat`/`Lon` dims of an array.

This is used to convert [`Selector`]($DDselectordocs) values form the user defined
projection to the underlying projection, and to show plot axes in the user projection.
"""
function usercrs end
usercrs(A::AbstractGeoArray) =
    if hasdim(A, Lat)
        usercrs(dims(A, Lat))
    elseif hasdim(A, Lon)
        usercrs(dims(A, Lon))
    else
        error("No Lat or Lon dimension, usercrs not available")
    end
usercrs(dim::Dimension) = usercrs(mode(dim), dim)

"""
    dimcrs(x)

Get the index crs projection of a [`Converted`](@ref) mode dim or
for the `Lat`/`Lon` dims of an array.

This is used for NetCDF and other formats where the underlying projection
of the data may  be different to what is contained in the vector index.
"""
function dimcrs end
dimcrs(A::AbstractGeoArray) =
    if hasdim(A, Lat)
        dimcrs(dims(A, Lat))
    elseif hasdim(A, Lon)
        dimcrs(dims(A, Lon))
    else
        error("No Lat or Lon dimension, dimcrs not available")
    end
dimcrs(dim::Dimension) = dimcrs(mode(dim), dim)

# DimensionalData methods

units(A::AbstractGeoArray) = getmeta(A, :units, nothing)

# Rebuild all types of AbstractGeoArray as GeoArray
rebuild(A::AbstractGeoArray, data, dims::Tuple, refdims, name, metadata, missingval=missingval(A)) =
    GeoArray(data, dims, refdims, name, metadata, missingval)
rebuild(A::AbstractGeoArray; data=data(A), dims=dims(A), refdims=refdims(A), name=name(A), metadata=metadata(A), missingval=missingval(A)) =
    GeoArray(data, dims, refdims, name, metadata, missingval)

Base.parent(A::AbstractGeoArray) = data(A)

"""
Abstract supertype for all memory-backed GeoArrays where the data is an array.
"""
abstract type MemGeoArray{T,N,D,A} <: AbstractGeoArray{T,N,D,A} end

"""
Abstract supertype for all disk-backed GeoArrays.
For these the data is lazyily loaded from disk.
"""
abstract type DiskGeoArray{T,N,D,A} <: AbstractGeoArray{T,N,D,A} end

filename(A::DiskGeoArray) = A.filename

"""
    data(f, A::DiskGeoArray)

Run method `f` on the data source object for A, as passed by the
`withdata` method for the array. The only requirement of the
object is that it has an `Array` method that returns the data as an array.
"""
data(A::DiskGeoArray) = withsourcedata(Array, A)

# Base methods

Base.size(A::DiskGeoArray) = A.size

Base.getindex(A::DiskGeoArray, i1::StandardIndices, i2::StandardIndices, I::StandardIndices...) =
    rebuildgetindex(A, i1, i2, I...)
Base.getindex(A::DiskGeoArray, i1::Integer, i2::Integer, I::Vararg{<:Integer}) =
    rawgetindex(A, i1, i2, I...)
# Linear indexing returns Array
Base.@propagate_inbounds Base.getindex(A::DiskGeoArray, i::Union{Colon,AbstractVector{<:Integer}}) =
    rawgetindex(A, i)
# Exempt 1D DimArrays
Base.@propagate_inbounds Base.getindex(A::DiskGeoArray{<:Any,1}, i::Union{Colon,AbstractVector{<:Integer}}) =
    rebuildgetindex(A, i)

rawgetindex(A, I...) =
    withsourcedata(A) do data
        readwindowed(data, I...)
    end

rebuildgetindex(A, I...) =
    withsourcedata(A) do data
        dims_, refdims_ = slicedims(dims(A), refdims(A), I)
        data = readwindowed(data, I...)
        rebuild(A, data, dims_, refdims_)
    end

Base.write(A::T) where T <: DiskGeoArray = write(filename(A), A)
Base.write(filename::AbstractString, A::T) where T <: DiskGeoArray =
    write(filename, T, A)

# Concrete implementation ######################################################

"""
    GeoArray(A::AbstractArray{T,N}, dims::Tuple;
             refdims=(), name=Symbol(""), metadata=nothing, missingval=missing)
    GeoArray(A::AbstractArray{T,N};
             dims, refdims=(), name=Symbol(""), metadata=nothing, missingval=missing)
    GeoArray(A::AbstractGeoArray; [data=data(A), dims=dims(A), refdims=refdims(A),
             name=name(A), metadata=metadata(A), missingval=missingval(A)]) =

A generic, memory-backed spatial array type. All [`AbstractGeoArray`](@ref) are
converted to `GeoArray` when indexed or otherwise transformed.

# Keyword Arguments

- `name`: `Symbol` name for the array.
- `dims`: `Tuple` of `Dimension`s for the array.
- `refdims`: `Tuple of` position `Dimension`s the array was sliced from,
  defaulting to `()`.
- `missingval`: Value reprsenting missing values, defaulting to `missing`.
  can be passed it.
- `metadata`: [`Metadata`](@ref) object for the array, or `nothing`.

## Example

```julia
A = GRDarray(gdalarray; name="surfacetemp")
# Select Australia using lat/lon coords, whatever the crs is underneath.
A[Lat(Between(-10, -43), Lon(Between(113, 153)))
```
"""
struct GeoArray{T,N,D<:Tuple,R<:Tuple,A<:AbstractArray{T,N},Na<:Symbol,Me,Mi} <: MemGeoArray{T,N,D,A}
    data::A
    dims::D
    refdims::R
    name::Na
    metadata::Me
    missingval::Mi
end
GeoArray(A::AbstractArray, dims, refdims, name::String, metadata, missingval=missing) = begin
    @warn "The GeoArray `name` field is now a Symbol"
    GeoArray(A, dims, refdims, Symbol(name), metadata, missingval)
end
@inline GeoArray(A::AbstractArray, dims::Tuple;
                 refdims=(), name=Symbol(""), metadata=nothing, missingval=missing) =
    GeoArray(A, formatdims(A, dims), refdims, name, metadata, missingval)
@inline GeoArray(A::AbstractArray; dims, refdims=(), name=Symbol(""), metadata=nothing, missingval=missing) =
    GeoArray(A, formatdims(A, dims), refdims, name, metadata, missingval)
@inline GeoArray(A::AbstractGeoArray; data=data(A), dims=dims(A), refdims=refdims(A),
                 name=name(A), metadata=metadata(A), missingval=missingval(A)) =
    GeoArray(data, dims, refdims, name, metadata, missingval)

Base.@propagate_inbounds Base.setindex!(A::GeoArray, x, I::Vararg{DimensionalData.StandardIndices}) =
    setindex!(data(A), x, I...)

Base.convert(::Type{GeoArray}, array::AbstractGeoArray) = GeoArray(array)