#include "MaxClass_Base.h"
#include "MaxBuffer.h"

#include <algorithm>
#include <cmath>
#include <vector>


#include <algorithm>
#include <cstdint>
#include <cmath>
#include <limits>
#include <thread>
#include <vector>

struct average_calculator
{
    enum class average_type { mean, rms, mean_hann, rms_hann };
    
    // Helper for making a window
    
    std::vector<double> make_hann_window(intptr_t window_size)
    {
        std::vector<double> window(window_size);
        
        double sum = 0.0;
        
        for (intptr_t i = 0; i < window_size; i++)
        {
            window[i] = 0.5 - (0.5 * cos(M_PI * 2. * ((double) i / (double) window_size)));
            sum += window[i];
        }
        
        sum = 1.0 / sum;
        
        for (intptr_t i = 0; i < window_size; i++)
            window[i] *= sum;
        
        return window;
    }
    
    // Retrieve samples from a temporary buffer
    
    static void get_samples(double *out, const double *buffer, intptr_t length, intptr_t offset, intptr_t nsamps)
    {
        intptr_t temp_offset = 0;
        intptr_t temp_nsamps = nsamps;
        intptr_t i;
        
        // Do not read before the buffer
        
        if (offset < 0)
        {
            temp_offset = -offset;
            temp_nsamps -= temp_offset;
            offset = 0;
        }
        
        if (temp_nsamps < 0)
        {
            temp_nsamps = 0;
            temp_offset = nsamps;
        }
        
        // Do not read beyond the buffer
        
        if (offset + temp_nsamps > length)
            temp_nsamps = length - offset;
        
        if (temp_nsamps < 0)
            temp_nsamps = 0;
        
        for (i = 0; i < temp_offset; i++)
            out[i] = buffer[0];
        
        for (i = 0; i < temp_nsamps; i++)
            out[i + temp_offset] = buffer[i + offset];
        
        for (i = temp_nsamps + temp_offset; i < nsamps; i++)
            out[i] = buffer[length - 1];
    }
    
    // Helpers for a window of samples
    
    template <bool Normalise, bool Root, bool MinusInfinity = false>
    struct calculation
    {
        double result(intptr_t window_size)
        {
            const double v = (Normalise ? m_value / window_size : m_value);
            return Root ? sqrt(v) : v;
        }
        
        double m_value = MinusInfinity ? -std::numeric_limits<double>::infinity() : 0.0;
    };
    
    template <bool Normalise, bool Root, bool MinusInfinity = false>
    struct window_calculation : calculation<Normalise, Root, MinusInfinity>
    {
        window_calculation(double *window) : m_window(window) {}
        double *m_window;
    };
    
    // Ops
    
    struct peak : calculation<false, false, true>
    {
        void operator()(double in, intptr_t)
        {
            const double mag = in;
            m_value = m_value < mag ? mag : m_value;
        }
    };
    
    struct average : calculation<true, false>
    {
        void operator()(double in, intptr_t) { m_value += in; }
    };
    
    struct rms : calculation<true, true>
    {
        void operator()(double in, intptr_t) { m_value += in * in; }
    };
    
    struct average_hann : window_calculation<false, false>
    {
        average_hann(double *window) : window_calculation(window) {}
        void operator()(double in, intptr_t idx) { m_value += in * m_window[idx]; }
    };
    
    struct rms_hann : window_calculation<false, true>
    {
        rms_hann(double *window) : window_calculation(window) {}
        void operator()(double in, intptr_t idx) { m_value += in * in * m_window[idx]; }
    };
    
    // Calculate one window
    
    template <class Calc, class...Args>
    static double get_window_calc(const double *input, double *samples, intptr_t in_length, intptr_t pos, intptr_t window_size, Args...args)
    {
        Calc op(args...);
        
        get_samples(samples, input, in_length, pos - (window_size >> 1), window_size);
        
        for (intptr_t i = 0; i < window_size; i++)
            op(samples[i], i);
        
        return op.result(window_size);
    }
    
    template <class Calc, class...Args>
    static void calc_op_loop(const double *input, double *output, intptr_t offset, intptr_t N, intptr_t length, intptr_t window_size, Args...args)
    {
        std::vector<double> samples(window_size);
        
        for (intptr_t i = offset; i < offset + N; i++)
            output[i] = get_window_calc<Calc>(input, samples.data(), length, i, window_size, args...);
    }
    
    template <class Calc, class...Args>
    void calc_op(double *input, double *output, intptr_t length, intptr_t window_size, Args...args)
    {
        int num_threads = std::thread::hardware_concurrency();
        intptr_t N = length / num_threads;
        
        std::vector<std::thread> threads;
        
        for (int i = 0; i < num_threads - 1; i++)
            threads.emplace_back(calc_op_loop<Calc, Args...>, input, output, i * N, N, length, window_size, args...);
        
        calc_op_loop<Calc>(input, output, (num_threads - 1) * N, length - (num_threads - 1) * N, length, window_size, args...);
        
        for (auto it = threads.begin(); it != threads.end() ; it++)
            it->join();
    }
    
    // Calculate ops for all samples
    
    void calc_peaks(double *input, double *output, intptr_t length, intptr_t window_size)
    {
        calc_op<peak>(input, output, length, window_size);
    }
    
    void calc_avg(double *input, double *output, intptr_t length, intptr_t window_size)
    {
        calc_op<average>(input, output, length, window_size);
    }
    
    void calc_rms(double *input, double *output, intptr_t length, intptr_t window_size)
    {
        calc_op<rms>(input, output, length, window_size);
    }
    
    void calc_rms_hann(double *input, double *output, intptr_t length, intptr_t window_size)
    {
        auto window = make_hann_window(window_size);
        calc_op<rms_hann>(input, output, length, window_size, window.data());
    }
    
    void calc_avg_hann(double *input, double *output, intptr_t length, intptr_t window_size)
    {
        auto window = make_hann_window(window_size);
        calc_op<average_hann>(input, output, length, window_size, window.data());
    }
    
    // Conversions
    
    static double dbtoa(double db)
    {
        return pow(10.0, db / 20.);
    }
    
    intptr_t ms_to_samples(double ms)
    {
        return static_cast<intptr_t>((m_sample_rate * ms) / 1000.0);
    }
    
    double samples_to_ms(intptr_t samples)
    {
        return (samples * 1000.0) / m_sample_rate;
    }
    
    average_calculator(double *input, intptr_t length, double sample_rate, double avg_time, average_type type, bool time_in_samps)
    : m_sample_rate(sample_rate)
    {
        intptr_t avg_window_size = time_in_samps ? avg_time : ms_to_samples(avg_time);
        
        m_buffer.resize(length);
        
        switch (type)
        {
            case average_type::mean:
                calc_avg(input, m_buffer.data(), length, avg_window_size);
                break;
                
            case average_type::mean_hann:
                calc_avg_hann(input, m_buffer.data(), length, avg_window_size);
                break;
                
            case average_type::rms:
                calc_rms(input, m_buffer.data(), length, avg_window_size);
                break;
                
            default:
                calc_rms_hann(input, m_buffer.data(), length, avg_window_size);
        }
    }
    
    size_t size() const { return m_buffer.size(); }
    const std::vector<double> get_buffer() const { return m_buffer; }
    const double& operator[](size_t idx) const { return m_buffer[idx]; }

    std::vector<double> m_buffer;
    double m_sample_rate = 0.0;
};


class FinesseSlices : public MaxClass_Base
{
public:
    
    static void classInit(t_class *c, t_symbol *nameSpace, const char *classname)
    {
        addMethod<FinesseSlices, &FinesseSlices::buffer>(c, "buffer");
        addMethod<FinesseSlices, &FinesseSlices::window>(c, "window");
        addMethod<FinesseSlices, &FinesseSlices::search>(c, "search");
        addMethod<FinesseSlices, &FinesseSlices::zerocross>(c, "zerocross");
        addMethod<FinesseSlices, &FinesseSlices::timeinsamps>(c, "timeinsamps");
        addMethod<FinesseSlices, &FinesseSlices::unordered>(c, "unordered");
        
        addMethod<FinesseSlices, &FinesseSlices::slices>(c, "list");

        addMethod<FinesseSlices, &FinesseSlices::assist>(c, "assist");
    }

    FinesseSlices(t_object *x, t_symbol *sym, long ac, t_atom *av)
    {
        mBuffer = ac > 0 ? atom_getsym(av + 0) : nullptr;
        mAverageTime = ac > 1 ? atom_getfloat(av + 1) : 10.0;
        mSearchTime = ac > 2 ? atom_getfloat(av + 2) : 20.0;
        mZeroCross = ac > 3 ? atom_getlong(av + 3) : 0;
        mTimeInSamples = ac > 4 ? atom_getlong(av + 4) : 0;
        mAllowUnordered = ac > 5 ? atom_getlong(av + 5) : 0;
        
        mOutlet = outlet_new(x, "list");
    }
    
    void assist(void *b, long m, long a, char *s)
    {
        if (m == ASSIST_OUTLET)
        {
            sprintf(s,"(list) New slices" );
        }
        else
        {
            sprintf(s,"(list) Slices / Set Buffer" );
        }
    }
    
    void buffer(t_symbol *s)
    {
        mBuffer = s;
    }
    
    void window(double time)
    {
        mAverageTime = time;
    }
    
    void search(double time)
    {
        mSearchTime = time;
    }
    
    void zerocross(t_atom_long flag)
    {
        mZeroCross = flag;
    }
    
    void timeinsamps(t_atom_long flag)
    {
        mTimeInSamples = flag;
    }
    
    void unordered(t_atom_long flag)
    {
        mAllowUnordered = flag;
    }
    
    void slices(t_symbol *sym, long ac, t_atom* av)
    {
        std::vector<double> input;
        std::vector<t_atom> outList(ac);
     
        MaxBufferAccess access(*this, mBuffer);
        
        if (!access.length())
            return;
        
        input.resize(access.length());
        access.read(input.data(), access.length(), 0, 0);
        
        // Get an RMS windowed summary of the buffer
        
        average_calculator average(input.data(), input.size(), access.sampleRate(), mAverageTime, average_calculator::average_type::rms_hann, mTimeInSamples);
        
        auto constrain = [&](long idx)
        {
            return std::max(0L, std::min(idx, static_cast<long>(access.length() - 1)));
        };
        
        auto search_samples = mTimeInSamples ? mSearchTime : average.ms_to_samples(mSearchTime);
        
        double last_value = -1.0;
        
        long valid_ac = ac;
        long next_arg = 0;

        for (long i = 0 ; i < ac; i++)
        {
            double value = atom_getfloat(av + i);
            
            // Reconsider position...
            
            long offset = constrain(std::round(mTimeInSamples ? value : average.ms_to_samples(value)));
            
            long lo = constrain(offset - (search_samples / 2));
            long hi = constrain(offset + (search_samples / 2));
            
            double min = average[offset];
            
            for (long j = lo + 1; j < hi - 1; j++)
            {
                if (average[j] < min && average[j] < average[j - 1] && average[j] < average[j + 1])
                {
                    min = average[j];
                    offset = j;
                    value = mTimeInSamples ? offset : average.samples_to_ms(offset);
                }
            }
            
            if (mZeroCross)
            {
                long cross1 = offset;
                long cross2 = offset;

                for (long j = offset; j > lo + 1; j--)
                {
                    if (!average[j] || (copysign(1.0, average[j]) == -copysign(1.0, average[j - 1])))
                    {
                        cross1 = j;
                        break;
                    }
                }
                
                for (long j = offset; j < hi - 1; j++)
                {
                    if (!average[j] || (copysign(1.0, average[j]) == -copysign(1.0, average[j + 1])))
                    {
                        cross2 = j;
                        break;
                    }
                }
                
                long cross = (offset - cross1) < (cross2 - offset) ? cross1: cross2;
                
                if (offset != cross)
                {
                    offset = cross;
                    value = mTimeInSamples ? offset : average.samples_to_ms(offset);
                }
            }
            
            if (mAllowUnordered || value > last_value)
            {
                atom_setfloat(outList.data() + next_arg, value);
                next_arg++;
            }
            else
                valid_ac--;
            
            last_value = value;
        }
        
        outlet_list(mOutlet, gensym("list"), ac, outList.data());
    }
                
    t_symbol *mBuffer;
    t_outlet *mOutlet;
    double mAverageTime;
    double mSearchTime;
    t_atom_long mZeroCross;
    t_atom_long mTimeInSamples;
    t_atom_long mAllowUnordered;
};

void ext_main(void *r)
{
    FinesseSlices::makeClass<FinesseSlices>(CLASS_BOX, "finesseslices~");
}