/* =========================STYLE-GUIDELINES===========================
 * All comments use MD-Syntax;
 * Formulas are written in KaTeX/TeX Syntax.
 * We adhere by the `kernel.org` style-guidelines.
 * See: https://www.kernel.org/doc/html/v4.10/process/coding-style.html
 * =============================AUTHORS================================
 * Created by Frederik Beimgraben, Minh Dan Cam and Lea Hornberger in
 * the winter term of 2024/25 at the faculty of computer science at
 * Reutlingen University.
 * After grading the full source code will be available at:
 * https://git.beimgraben.net/MMS-2024-25/MMS-Loesung-2024-25.git
 */

#include <stdio.h>
#include <stdarg.h>
#include <stdlib.h>

#include "MMS24-25.h"
#include "math.h"

// #region Error Handling and Debugging
void _error(char *fmt, ...) {
    fprintf(stderr, "< \x1B[31;4;1mERROR\x1B[0m : \x1B[34m'");

    va_list args;
    va_start(args, fmt);
    vfprintf(stderr, fmt, args);
    va_end(args);

    fprintf(stderr, "'\x1B[0m >\n");

#ifndef _TESTS_NONSTOP
    exit(EXIT_FAILURE);
#endif
}

void _debug(char *fmt, ...) {
    fprintf(stderr, "[ \x1B[37;4;1mDEBUG\x1B[0m ] \x1B[34m'");

    va_list args;
    va_start(args, fmt);
    vfprintf(stderr, fmt, args);
    va_end(args);

    fprintf(stderr, "'\x1B[0m\n");
}

void _warn(char *fmt, ...) {
    fprintf(stderr, "[ \x1B[33;4;1mWARN\x1B[0m ] \x1B[34m'");

    va_list args;
    va_start(args, fmt);
    vfprintf(stderr, fmt, args);
    va_end(args);

    fprintf(stderr, "'\x1B[0m\n");
}
// #endregion

// #region Aufgabe 1
double interpolateLine(
        double dX1,
        double dY1,
        double dX2,
        double dY2,
        double dXq
    ) {
    if (dX2 == dX1)
        _error("interpolateLine: dX1 and dX2 must differ!");

    double y = dY1 * (
        (dX2 - dXq) / (dX2 - dX1)
    ) + dY2 * (
        (dXq - dX1) / (dX2 - dX1)
    );

    #ifdef _DEBUG_
        // DEBUG: print the parameters and result
        _debug("interpolateLine(%lf, %lf, %lf, %lf, %lf) -> %lf\n",
           dX1,dY1,dX2,dY2,dXq,y);
    #endif

    return y;
}

void rescaleDoubleArray(
        double *pdIn,
        int iLen,
        double dMin,
        double dFactor,
        double *pdOut
    ) {
    /* Rescale an array of double values, so that the following rules apply:
     * Given $ |pdIn| = |pdOut| = iLen $ it should be, that:
     * $ min(pdIn) > dMin $, as well as
     * \[ \forall i, j \in (0,\dots,iLen-1):
     *        (pdOut_i - pdOut_j) = (pdIn_i - pdIn_j)  *dFactor \],
     * which naturally results from:
     * $ pdOut_i = pdIn_i  *dFactor + dMin $
     */
    // First find the old minimum
    double dOldMin = pdIn[0];
    for (int i = 1; i < iLen; i++) {
        if (pdIn[i] < dOldMin)
            dOldMin = pdIn[i];
    }

    // Calculate the offset
    double dOffset = dMin - dOldMin;

    // Rescale the array
    for (int i = 0; i < iLen; i++) {
        pdOut[i] = pdIn[i] * dFactor + dOffset;
    }
}

double *generateSineValues(
        int iNumValues,
        int iNumSamplesPerPeriod,
        double dAmp
    ) {
    /* Generate an array of `iNumValues` double value, so that the values
     * represent the Y-coordinates of a sine curve of which one full period
     * occupies `iNumSamplesPerPeriod` values.
     */

    if (iNumSamplesPerPeriod == 0)
        _error("*generateSineValues: iNumSamplesPerPeriod must be non-0!");

    if (iNumSamplesPerPeriod < 0)
        _warn("*generateSineValues: iNumSamplesPerPeriod is less than 0!");

    if (iNumValues < 0)
        _error("*generateSineValues: iNumValues must be equal to or greater than 0!");

    if (iNumValues == 0)
        _warn("*generateSineValues: iNumValues is 0!");


    // Initialize the output array
    double *sineArray = NULL;
    sineArray = (double*) malloc(sizeof(double) * iNumValues);

    if (sineArray == NULL)
        _error("*generateSineValues: Failed to allocate memory for Sine Array\n");

    for (int i=0; i<iNumValues;i++){
        sineArray[i] = dAmp * sin(((2*M_PI)/iNumSamplesPerPeriod) * i);
    }

    return sineArray;
}

void writeDoubleArray(double *pdOut, int iLen, char *pcOutName) {
    // Write a list of double values to a file
    FILE *file = fopen(pcOutName, "w");
    if (file == NULL)
        _error("writeDoubleArray: Error while opening file '%s'.\n", pcOutName);

    if (pdOut == NULL)
        _error("writeDoubleArray: pdOut is NULL\n");

    if (iLen == 0)
        _warn("writeDoubleArray: iLen is 0\n");

    for (int i = 0; i < iLen; i++) {
        fprintf(file, "%f\n", pdOut[i]);
    }

    fclose(file);
}

int readDoubleArray(char *pcInName, double **ppdIn) {
    char *line = NULL;
    size_t len = 0;
    ssize_t read;
    int i = 0;

    *ppdIn = NULL;

    FILE *file = fopen(pcInName, "r");
    if (file == NULL)
        _error("readDoubleArray: Error while opening file '%s'.\n", pcInName);

    while ((read = getline(&line, &len, file)) != -1) {
        if (*ppdIn == NULL) {
            *ppdIn = (double*) malloc(sizeof(double));
        } else {
            *ppdIn = (double*) realloc(*ppdIn, sizeof(double) *(i+1));
        }

        if (*ppdIn == NULL)
            _error("readDoubleArray: Failed to allocate or re-allocate memory for pdIn\n");

        (*ppdIn)[i] = atof(line);
        i++;
    }

    free(line);
    fclose(file);

    return i;
}

MMSignal *createSignal() {
    MMSignal *signal = NULL;
    signal = (MMSignal*) malloc(sizeof(MMSignal));

    if (signal == NULL)
        _error("createSignal: Failed to allocate memory for signal\n");

    signal->pdValues = NULL;
    signal->iNumValues = 0;
    signal->dArea = 0;
    signal->dMean = 0;
    signal->dStdDev = 0;
    signal->dMedian = 0;
    signal->pexExtrema = NULL;

    return signal;
}

LocalExtrema* copyLocalExtrema(LocalExtrema *pexIn) {
    if (pexIn == NULL)
        return NULL;

    LocalExtrema *localExtrema = (LocalExtrema*) malloc(sizeof(LocalExtrema));
    if (localExtrema == NULL)
        _error("copyLocalExtrema: Failed to allocate memory for localExtrema\n");

    localExtrema->iNumLocalExtrema = pexIn->iNumLocalExtrema;
    localExtrema->piLocalExtremaPos = (int*) malloc(sizeof(int) *pexIn->iNumLocalExtrema);
    if (localExtrema->piLocalExtremaPos == NULL)
        _error("copyLocalExtrema: Failed to allocate memory for localExtrema->piLocalExtremaPos\n");

    for (int i = 0; i < pexIn->iNumLocalExtrema; i++) {
        localExtrema->piLocalExtremaPos[i] = pexIn->piLocalExtremaPos[i];
    }

    return localExtrema;
}

Extrema* copyExtrema(Extrema *pexIn) {
    if (pexIn == NULL)
        return NULL;

    Extrema *extrema = (Extrema*) malloc(sizeof(Extrema));
    if (extrema == NULL)
        _error("copyExtrema: Failed to allocate memory for extrema\n");

    extrema->iGlobalMinPos = pexIn->iGlobalMinPos;
    extrema->iGlobalMaxPos = pexIn->iGlobalMaxPos;

    extrema->pexLocalMax = copyLocalExtrema(pexIn->pexLocalMax);
    extrema->pexLocalMin = copyLocalExtrema(pexIn->pexLocalMin);

    return extrema;
}

MMSignal *createSignalFromSignal(MMSignal *pmmsIn) {
    MMSignal *signal = createSignal();

    if (pmmsIn->pdValues != NULL) {
        // Initialize and copy the values (deep copy)
        signal->pdValues = (double*) malloc(sizeof(double) *pmmsIn->iNumValues);
        if (signal->pdValues == NULL)
            _error("createSignalFromSignal: Failed to allocate memory for signal->pdValues\n");

        for (int i = 0; i < pmmsIn->iNumValues; i++) {
            signal->pdValues[i] = pmmsIn->pdValues[i];
        }
    } else {
        _warn("createSignalFromSignal: pmmsIn->pdValues is NULL\n");
        signal->pdValues = NULL;
    }

    // Copy shallow items
    signal->iNumValues = pmmsIn->iNumValues;
    signal->dArea = pmmsIn->dArea;
    signal->dMean = pmmsIn->dMean;
    signal->dStdDev = pmmsIn->dStdDev;
    signal->dMedian = pmmsIn->dMedian;

    // Deep copy of Extrema
    if (pmmsIn->pexExtrema == NULL) {
        _warn("createSignalFromSignal: pmmsIn->pexExtrema is NULL\n");
        signal->pexExtrema = NULL;
        return signal;
    }

    signal->pexExtrema = copyExtrema(pmmsIn->pexExtrema);

    return signal;
}

MMSignal *createSignalWithDefault(int iArrayLen, double dDefaultValue) {
    MMSignal *signal = createSignal();

    signal->pdValues = (double*) malloc(sizeof(double) *iArrayLen);
    if (signal->pdValues == NULL)
        _error("createSignalWithDefault: Failed to allocate memory for signal->pdValues\n");

    for (int i = 0; i < iArrayLen; i++) {
        signal->pdValues[i] = dDefaultValue;
    }

    signal->iNumValues = iArrayLen;

    return signal;
}

MMSignal *createSignalFromDoubleArray(int iArrayLen, double *pdIn) {
    MMSignal *signal = createSignal();

    signal->pdValues = (double*) malloc(sizeof(double) *iArrayLen);
    if (signal->pdValues == NULL)
        _error("createSignalFromDoubleArray: Failed to allocate memory for signal->pdValues\n");

    for (int i = 0; i < iArrayLen; i++) {
        signal->pdValues[i] = pdIn[i];
    }

    signal->iNumValues = iArrayLen;

    return signal;
}

MMSignal *createSignalFromFile(char *pcInName) {
    double *pdIn = NULL;
    int iArrayLen = readDoubleArray(pcInName, &pdIn);

    MMSignal *signal = createSignal();

    signal->pdValues = (double*) malloc(sizeof(double) *iArrayLen);
    if (signal->pdValues == NULL)
        _error("createSignalFromFile: Failed to allocate memory for signal->pdValues\n");

    for (int i = 0; i < iArrayLen; i++) {
        signal->pdValues[i] = pdIn[i];
    }

    signal->iNumValues = iArrayLen;

    return signal;
}

void writeMMSignal(char *pcInName, MMSignal *pmmsIn) {
    if (pmmsIn->pdValues == NULL)
        _error("writeMMSignal: pmmsIn->pdValues is NULL\n");

    writeDoubleArray(
        pmmsIn->pdValues,
        pmmsIn->iNumValues,
        pcInName
    );
}

void deleteLocalExtrema(LocalExtrema *pexIn) {
    if (pexIn == NULL)
        return;

    free(pexIn->piLocalExtremaPos);
    free(pexIn);
}

void deleteExtrema(Extrema *pexIn) {
    if (pexIn == NULL)
        return;

    deleteLocalExtrema(pexIn->pexLocalMax);
    deleteLocalExtrema(pexIn->pexLocalMin);
    free(pexIn);
}

void deleteSignal(MMSignal *pmmsIn) {
    if (pmmsIn == NULL)
        return;

    if (pmmsIn->pdValues != NULL)
        free(pmmsIn->pdValues);

    deleteExtrema(pmmsIn->pexExtrema);

    free(pmmsIn);
}
// #endregion

// #region Aufgabe 2
double calculateArea(double *pdIn, int iLen) {
    double area = 0;

    // #FIXME: Which one should we use according to the task?
    for (int i = 0; i < iLen; i++) {
        /* Alternative would be:
         * area += (pdIn[i] + pdIn[i+1]) / 2;
         * with i < iLen-1
         * Which one to use is up to the implemention.
         */
        area += pdIn[i];
    }

    return area;
}

double calculateMean(double *pdIn, int iLen) {
    /* Here the above defined function `calculateArea` is used to calculate the mean, since:
     * \[ mean(pdIn) = \frac{1}{ |pdIn| } \sum_{i=1}^{ |pdIn| } x_i
     *               = \frac{area(x)}{ |pdIn| } \]
     */
    return calculateArea(pdIn, iLen) / iLen;
}

double calculateStddev(double *pdIn, int iLen) {
    /* The standard deviation is calculated as follows:
     * \[ stddev(pdIn) = \sqrt{ \frac{1}{ |pdIn| } \sum_{i=1}^{ |pdIn| } (x_i - mean(pdIn))^2 } \]
     */

    double mean = calculateMean(pdIn, iLen);
    double stdDev = 0;

    for (int i = 0; i < iLen; i++) {
        stdDev += pow(pdIn[i] - mean, 2);
    }

    return sqrt(stdDev / iLen);
}

double calculateMedian(double *pdIn, int iLen) {
    /* The median is calculated as follows:
     * \[ median(pdIn) = \begin{cases}
     *      pdIn_{\frac{|pdIn|}{2}} & \text{if } |pdIn| \text{ is even} \\
     *      \frac{pdIn_{\frac{|pdIn|}{2}} + pdIn_{\frac{|pdIn|}{2}+1}}{2} & \text{if } |pdIn| \text{ is odd}
     * \end{cases} \]
     */

    if (iLen % 2 == 0) {
        return pdIn[iLen/2];
    } else {
        return (pdIn[iLen/2] + pdIn[iLen/2+1]) / 2;
    }
}

LocalExtrema *initLocalExtrema(int iLen) {
    LocalExtrema *localExtrema = (LocalExtrema*) malloc(sizeof(LocalExtrema));
    if (localExtrema == NULL)
        _error("initLocalExtrema: Failed to allocate memory for localExtrema\n");

    localExtrema->iNumLocalExtrema = 0;
    localExtrema->piLocalExtremaPos = (int*) malloc(sizeof(int) *iLen);
    if (localExtrema->piLocalExtremaPos == NULL)
        _error("initLocalExtrema: Failed to allocate memory for localExtrema->piLocalExtremaPos\n");

    return localExtrema;
}

Extrema *initExtrema(double *pdIn, int iLen) {
    Extrema *extrema = (Extrema*) malloc(sizeof(Extrema));
    if (extrema == NULL)
        _error("initExtrema: Failed to allocate memory for extrema\n");

    extrema->pexLocalMax = initLocalExtrema(iLen);
    extrema->pexLocalMin = initLocalExtrema(iLen);

    extrema->iGlobalMaxPos = 0;
    extrema->iGlobalMinPos = 0;

    return extrema;
}

void initMMSignalFeatures(MMSignal *pmmsIn) {
    pmmsIn->dArea = 0;
    pmmsIn->dMean = 0;
    pmmsIn->dStdDev = 0;
    pmmsIn->dMedian = 0;
    pmmsIn->pexExtrema = initExtrema(pmmsIn->pdValues, pmmsIn->iNumValues);
}

Histogram *initHistogram(double *pdIn, int iLen, int iNumBins) {
    Histogram *histogram = (Histogram*) malloc(sizeof(Histogram));
    if (histogram == NULL)
        _error("initHistogram: Failed to allocate memory for histogram");

    if (iNumBins == 0)
        _error("initHistogram: iNumBins must be non-0");

    if (iLen == 0)
        _warn("initHistogram: iLen is 0\n");

    histogram->dIntervalMin = pdIn[0];
    histogram->dIntervalMax = pdIn[0];
    for (int i = 0; i < iLen; i++) {
        if (pdIn[i] < histogram->dIntervalMin)
            histogram->dIntervalMin = pdIn[i];
        if (pdIn[i] > histogram->dIntervalMax)
            histogram->dIntervalMax = pdIn[i];
    }

    #ifdef _DEBUG_
    _debug("initHistogram: histogram->dIntervalMin = %lf", histogram->dIntervalMin);
    _debug("initHistogram: histogram->dIntervalMax = %lf", histogram->dIntervalMax);
    #endif

    histogram->dBinWidth = (histogram->dIntervalMax - histogram->dIntervalMin) / iNumBins;
    histogram->iNumBins = iNumBins;

    #ifdef _DEBUG_
    _debug("initHistogram: histogram->dBinWidth = %lf", histogram->dBinWidth);
    #endif

    if (histogram->dBinWidth == 0)
        histogram->dBinWidth = 1;

    histogram->pdBinValues = (double*) malloc(sizeof(double) *iNumBins);
    if (histogram->pdBinValues == NULL)
        _error("initHistogram: Failed to allocate memory for histogram->pdBinValues");

    for (int i = 0; i < iNumBins; i++) {
        histogram->pdBinValues[i] = 0;
    }

    for (int i = 0; i < iLen; i++) {
        int bin = floor((pdIn[i] - histogram->dIntervalMin) / histogram->dBinWidth);
        if (bin == iNumBins)
            bin--;

        histogram->pdBinValues[bin]++;
    }

    for (int i = 0; i < iNumBins; i++) {
        histogram->pdBinValues[i] /= iLen;
    }

    #ifdef _DEBUG_
    // DEBUG: print the parameters and result
    _debug("initHistogram(%lf, %lf, %lf, %d)",
        histogram->dIntervalMin,histogram->dIntervalMax,histogram->dBinWidth,histogram->iNumBins);
    for (int i = 0; i < histogram->iNumBins; i++) {
        _debug("initHistogram: histogram->pdBinValues[%d] = %lf", i, histogram->pdBinValues[i]);
    }
    #endif

    return histogram;
}

void deleteHistogram(Histogram *phIn) {
    if (phIn == NULL)
        return;

    free(phIn->pdBinValues);
    free(phIn);
}

double calculateEntropy(Histogram *phisIn) {
    double entropy = 0;

    for (int i = 0; i < phisIn->iNumBins; i++) {
        if (phisIn->pdBinValues[i] != 0) {
            entropy += phisIn->pdBinValues[i] * log2(phisIn->pdBinValues[i]);
        }
    }

    return -entropy;
}

// #endregion

// #region Aufgabe 3

MMSignal *convoluteSignals(MMSignal *pmmsInA, MMSignal *pmmsInB) {
    // Calculate output length
    int iOutLen = pmmsInA->iNumValues + pmmsInB->iNumValues - 1;

    double *pdOut = (double*) malloc(sizeof(double) *iOutLen);
    if (pdOut == NULL)
        _error("convoluteSignals: Failed to allocate memory for pdOut\n");

    for (int i = 0; i < iOutLen; i++) {
        pdOut[i] = 0;
    }

    for (int i = 0; i < pmmsInA->iNumValues; i++) {
        for (int j = 0; j < pmmsInB->iNumValues; j++) {
            pdOut[i+j] += pmmsInA->pdValues[i] * pmmsInB->pdValues[j];
        }
    }

    MMSignal *signal = createSignalFromDoubleArray(iOutLen, pdOut);
    initMMSignalFeatures(signal);

    return signal;
}

// #FIXME: Check if this checks out
MMSignal *getPascalLine(int iLinenum) {
    if (iLinenum <= 0)
        _error("getPascalLine: iLinenum must be greater than 0\n");

    // Calculate a kernel signal using pascales triangle without using the factorial
    double *pdValues = NULL;
    pdValues = (double*) malloc(sizeof(double) * iLinenum);
    if (pdValues == NULL)
        _error("getPascalLine: Failed to allocate memory for pdValues\n");

    for (int i = 0; i < iLinenum; i++) {
        pdValues[i] = 1;
        for (int j = i-1; j > 0; j--) {
            pdValues[j] = pdValues[j] + pdValues[j-1];
        }
    }

    MMSignal *signal = createSignalFromDoubleArray(iLinenum, pdValues);
    initMMSignalFeatures(signal);

    return signal;
}

// #endregion

// #region Aufgabe 4
void computeDFT(
        int iLen,
    	double *pdRealIn, double *pdImgIn,
    	double *pdRealOut, double *pdImgOut,
    	int iDirection
    ) {
	// Check the value of iDirection
	if (iDirection != -1 && iDirection != 1)
        _error("computeDFT: iDirection must be either -1 or 1\n");

	for (int k = 0; k < iLen; k++) {
        pdRealOut[k] = 0;
        pdImgOut[k] = 0;

        // Compute the sum for the k-th output element
        for (int n = 0; n < iLen; n++) {
            double angle = iDirection * 2.0 * M_PI * k * n / iLen;
            double cosTerm = cos(angle);
            double sinTerm = sin(angle);

            pdRealOut[k] += pdRealIn[n] * cosTerm - pdImgIn[n] * sinTerm;
            pdImgOut[k] += pdRealIn[n] * sinTerm + pdImgIn[n] * cosTerm;
        }

        // Normalize the output
        if (iDirection == -1) {
            pdRealOut[k] /= iLen;
            pdImgOut[k] /= iLen;
        }
    }
}

void convertCart2Polar(double *pdRealIn, double *pdImgIn, double *pdRadiusOut, double *pdAngleOut, int iLen) {
    for (int i = 0; i < iLen; i++) {
        pdRadiusOut[i] = sqrt(pow(pdRealIn[i], 2) + pow(pdImgIn[i], 2));
        pdAngleOut[i] = atan2(pdImgIn[i], pdRealIn[i]);
    }
}

void convertPolar2Cart(double *pdRadiusIn, double *pdAngleIn, double *pdRealOut, double *pdImgOut, int iLen) {
    for (int i = 0; i < iLen; i++) {
        pdRealOut[i] = pdRadiusIn[i] * cos(pdAngleIn[i]);
        pdImgOut[i] = pdRadiusIn[i] * sin(pdAngleIn[i]);
    }
}
// #endregion