/**++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 *
 *                                              PANIC
 *                              Portable Algorithms and Numerics In C++
 *
 *                       Scientific computing from scratch, with feeling.
 *
 * Copyright (c) 2026 Michelle Bausager
 *
 * This file is part of PANIC.
 *
 * PANIC is free software licensed under the GNU General Public License v3.0 or later.
 * You may redistribute and/or modify it under the terms of the GPL.
 *
 * PANIC is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
 * without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 * See the LICENSE file for the full license text.
 *
 * SPDX-License-Identifier: GPL-3.0-or-later
 *
 *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 *
 * Project Name:      PANIC
 * Module Name:       <module-name>
 * File Name:         <file-name>
 * Revision:          0.1.0
 * Date:              <yyyy-mm-dd>
 * Author:            Michelle Bausager
 *
 * Description:
 * <short description of what this file/module does>
 *
 *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
//---------------------------------------------------------------------------------------------------------------------------
// INCLUDE DESCRIPTION
//---------------------------------------------------------------------------------------------------------------------------
#include <config/omp.hpp>
#include <iostream> // std::cout, std::endl
#include <config/types.hpp> // include types to use
#include <config/constants.hpp> // Math constants
#include <tensor/vector.hpp>
#include <tensor/matrix.hpp>
#include <io/print_tensor.hpp>
#include <math/matmul.hpp>
#include <neural_network/layer/layer_dense.hpp>
#include <math/add.hpp>
#include <random/uniform.hpp>
#include <neural_network/layer/layer_dense.hpp>
#include <neural_network/model/model.hpp>

#include <math/mul.hpp>
#include <math/clip.hpp>
#include <neural_network/activation/activation_relu.hpp>
#include <tensor/generators/linspace.hpp>
#include <math/trigonometry/sin.hpp>
#include <neural_network/datasets/sine_data.hpp>
#include <neural_network/datasets/spiral_data.hpp>
#include <neural_network/datasets/vertical_data.hpp>
#include <math/exp.hpp>
#include <math/sub.hpp>
#include <neural_network/loss/loss_categorical_crossentropy.hpp>

#include <math.h>



// For omp tesing:
#include <chrono>
#include <cstddef>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <limits>

#if PANIC_HAS_OPENMP
    #include <omp.h>
#endif

//---------------------------------------------------------------------------------------------------------------------------
// DEFINE DESCRIPTION
//---------------------------------------------------------------------------------------------------------------------------
 // difinitions like:
// #define TEST_FALG              1











//---------------------------------------------------------------------------------------------------------------------------
// Function Name : benchmark_omp_min_work
//
// Description:
//      Measures an operation with one thread and with all available threads.
//
//      It then tests several possible omp_min_work values and estimates which
//      threshold best matches the measured results.
//---------------------------------------------------------------------------------------------------------------------------

void benchmark_omp_min_work(){

#if !PANIC_HAS_OPENMP

    std::cout << "OpenMP is not enabled.\n";

#else

    //-----------------------------------------------------------------------------------------------------------------------
    // OPENMP SETTINGS
    //-----------------------------------------------------------------------------------------------------------------------

    omp_set_dynamic(0);

    const int parallel_thread_count = omp_get_max_threads();


    //-----------------------------------------------------------------------------------------------------------------------
    // EDIT 1:
    // Sizes to test.
    //
    // This list is suitable for vector operations.
    //-----------------------------------------------------------------------------------------------------------------------

    const panic::types::uint_t sizes[] = {
        100,
        250,
        500,
        1000,
        2000,
        4000,
        8000,
        12000,
        16000,
        24000,
        32000,
        50000,
        75000,
        100000,
        150000,
        250000,
        500000,
        1000000,
        2000000,
        5000000
    };

    /*
    // Example sizes for square matrix multiplication:

    const panic::types::uint_t sizes[] = {
        4,
        6,
        8,
        10,
        12,
        16,
        20,
        24,
        32,
        48,
        64,
        96,
        128,
        192,
        256,
        384,
        512
    };
    */

    const std::size_t number_of_sizes =
        sizeof(sizes) / sizeof(sizes[0]);


    //-----------------------------------------------------------------------------------------------------------------------
    // EDIT 2:
    // Possible omp_min_work values to test.
    //
    // These are work values, not necessarily vector or matrix sizes.
    //-----------------------------------------------------------------------------------------------------------------------

	const std::uint64_t omp_min_work_values[] = {
	    0,
	    50,
	    100,
	    150,
	    200,
	    250,
	    300,
	    350,
	    400,
	    450,
	    500,
	    600,
	    750,
	    1000,
	    1500,
	    2000,
	    4000,
	    8000,
	    16000,
	    32000
	};

    /*
    // Example omp_min_work values for matrix multiplication:

    const std::uint64_t omp_min_work_values[] = {
        0,
        1000,
        5000,
        10000,
        25000,
        50000,
        100000,
        250000,
        500000,
        1000000,
        2000000,
        5000000,
        10000000,
        25000000,
        50000000,
        100000000
    };
    */

    const std::size_t number_of_omp_min_work_values =
        sizeof(omp_min_work_values) /
        sizeof(omp_min_work_values[0]);


    //-----------------------------------------------------------------------------------------------------------------------
    // BENCHMARK SETTINGS
    //-----------------------------------------------------------------------------------------------------------------------

    // The benchmark tries to perform approximately this much work for each size.
    const std::uint64_t target_total_work = 50000000;

    // Minimum and maximum number of repetitions for each size.
    const std::uint64_t minimum_repetitions = 3;
    const std::uint64_t maximum_repetitions = 1000;

    // Prevent the compiler from treating all calculated results as unused.
    volatile panic::types::real_t benchmark_sink = 0.0f;


    //-----------------------------------------------------------------------------------------------------------------------
    // ARRAYS FOR THE MEASURED RESULTS
    //-----------------------------------------------------------------------------------------------------------------------

    std::uint64_t measured_work[number_of_sizes];

    double one_thread_results[number_of_sizes];
    double parallel_results[number_of_sizes];


    //-----------------------------------------------------------------------------------------------------------------------
    // PRINT BENCHMARK INFORMATION
    //-----------------------------------------------------------------------------------------------------------------------

	std::cout
	    << std::right
	    << std::setw(12) << "Size"
	    << std::setw(16) << "Work"
	    << std::setw(14) << "Repetitions"
	    << std::setw(18) << "1 thread (us)"
	    << std::setw(18) << "Parallel (us)"
	    << std::setw(12) << "Speedup"
	    << std::setw(14) << "Fastest"
	    << "\n";

	std::cout
	    << std::string(104, '-')
	    << "\n";


    //-----------------------------------------------------------------------------------------------------------------------
    // TEST EVERY SIZE
    //-----------------------------------------------------------------------------------------------------------------------

    for (std::size_t size_index = 0;
         size_index < number_of_sizes;
         ++size_index){

        const panic::types::uint_t size =
            sizes[size_index];


        //-------------------------------------------------------------------------------------------------------------------
        // EDIT 3:
        // Initialize the vectors or matrices used by the operation.
        //-------------------------------------------------------------------------------------------------------------------

        // Vector-add example:

        panic::tensor::real_vector a(size, 1.0f);
        panic::tensor::real_vector b(size, 2.0f);
        panic::tensor::real_vector c(size);


        /*
        // Square-matrix multiplication example:

        panic::tensor::real_matrix a(
            size,
            size,
            0.01f
        );

        panic::tensor::real_matrix b(
            size,
            size,
            0.02f
        );

        panic::tensor::real_matrix c(
            size,
            size
        );
        */


        //-------------------------------------------------------------------------------------------------------------------
        // EDIT 4:
        // Calculate work in the same way as the function being tested.
        //-------------------------------------------------------------------------------------------------------------------

        // Vector operation:

        const std::uint64_t work =
            static_cast<std::uint64_t>(size);


        /*
        // Square matrix add:

        const std::uint64_t work =
            static_cast<std::uint64_t>(size) *
            static_cast<std::uint64_t>(size);
        */


        /*
        // Square matrix multiplication:

        const std::uint64_t work =
            static_cast<std::uint64_t>(size) *
            static_cast<std::uint64_t>(size) *
            static_cast<std::uint64_t>(size);
        */


        //-------------------------------------------------------------------------------------------------------------------
        // CALCULATE NUMBER OF REPETITIONS
        //-------------------------------------------------------------------------------------------------------------------

        std::uint64_t repetitions =
            target_total_work / work;

        if (repetitions < minimum_repetitions){
            repetitions = minimum_repetitions;
        }

        if (repetitions > maximum_repetitions){
            repetitions = maximum_repetitions;
        }


        //-------------------------------------------------------------------------------------------------------------------
        // MEASURE WITH ONE THREAD AND ALL THREADS
        //
        // test == 0: one OpenMP thread
        // test == 1: all available OpenMP threads
        //-------------------------------------------------------------------------------------------------------------------

        double measured_time_us[2] = {
            0.0,
            0.0
        };

        for (int test = 0; test < 2; ++test){

            if (test == 0){
                omp_set_num_threads(1);
            }
            else{
                omp_set_num_threads(
                    parallel_thread_count
                );
            }


            //----------------------------------------------------------------------------------------------------------------
            // WARM-UP
            //----------------------------------------------------------------------------------------------------------------

            for (std::size_t warmup = 0;
                 warmup < 2;
                 ++warmup){

                //------------------------------------------------------------------------------------------------------------
                // EDIT 5:
                // Put the operation being tested here.
                //------------------------------------------------------------------------------------------------------------

                panic::math::add(a, b, c);

                // Matrix multiplication:
                // panic::math::matmul(a, b, c);
            }


            //----------------------------------------------------------------------------------------------------------------
            // TIMED LOOP
            //----------------------------------------------------------------------------------------------------------------

            const std::chrono::steady_clock::time_point start =
                std::chrono::steady_clock::now();

            for (std::uint64_t repetition = 0;
                 repetition < repetitions;
                 ++repetition){

                //------------------------------------------------------------------------------------------------------------
                // EDIT 6:
                // Put the same operation here.
                //------------------------------------------------------------------------------------------------------------

                panic::math::add(a, b, c);

                // Matrix multiplication:
                // panic::math::matmul(a, b, c);
            }

            const std::chrono::steady_clock::time_point end =
                std::chrono::steady_clock::now();


            //----------------------------------------------------------------------------------------------------------------
            // CALCULATE AVERAGE TIME PER OPERATION
            //----------------------------------------------------------------------------------------------------------------

            const double total_time_us =
                std::chrono::duration<double, std::micro>(
                    end - start
                ).count();

            measured_time_us[test] =
                total_time_us /
                static_cast<double>(repetitions);


            //----------------------------------------------------------------------------------------------------------------
            // READ ONE RESULT
            //
            // Change this if the output cannot be accessed with c[size / 2].
            //----------------------------------------------------------------------------------------------------------------

            benchmark_sink += c[size / 2];

            /*
            // For a matrix:

            benchmark_sink += c(
                size / 2,
                size / 2
            );
            */
        }


        //-------------------------------------------------------------------------------------------------------------------
        // STORE RESULTS
        //-------------------------------------------------------------------------------------------------------------------

        const double one_thread_us =
            measured_time_us[0];

        const double parallel_us =
            measured_time_us[1];

        measured_work[size_index] = work;
        one_thread_results[size_index] = one_thread_us;
        parallel_results[size_index] = parallel_us;


        //-------------------------------------------------------------------------------------------------------------------
        // PRINT RESULTS FOR THIS SIZE
        //-------------------------------------------------------------------------------------------------------------------

        const double speedup =
            one_thread_us / parallel_us;

        const char* fastest;

        if (parallel_us < one_thread_us){
            fastest = "parallel";
        }
        else{
            fastest = "one_thread";
        }

		std::cout
		    << std::right
		    << std::setw(12) << size
		    << std::setw(16) << work
		    << std::setw(14) << repetitions
		    << std::setw(18) << std::fixed << std::setprecision(3)
		    << one_thread_us
		    << std::setw(18) << parallel_us
		    << std::setw(12) << speedup
		    << std::setw(14) << fastest
		    << "\n";
    }


    //-----------------------------------------------------------------------------------------------------------------------
    // TEST THE POSSIBLE OMP_MIN_WORK VALUES
    //-----------------------------------------------------------------------------------------------------------------------

    std::uint64_t best_omp_min_work = 0;

    double best_score =
        std::numeric_limits<double>::max();

    std::cout
        << "\n"
        << "omp_min_work"
        << ", average_slowdown"
        << "\n";

    for (std::size_t threshold_index = 0;
         threshold_index < number_of_omp_min_work_values;
         ++threshold_index){

        const std::uint64_t omp_min_work =
            omp_min_work_values[threshold_index];

        double score = 0.0;


        //-------------------------------------------------------------------------------------------------------------------
        // SEE WHICH VERSION THIS THRESHOLD WOULD SELECT
        //-------------------------------------------------------------------------------------------------------------------

        for (std::size_t size_index = 0;
             size_index < number_of_sizes;
             ++size_index){

            double selected_time_us;


            //----------------------------------------------------------------------------------------------------------------
            // This uses > because the PANIC functions currently use:
            //
            //     work > omp_min_work
            //----------------------------------------------------------------------------------------------------------------

            if (measured_work[size_index] > omp_min_work){

                // This threshold would select OpenMP.
                selected_time_us =
                    parallel_results[size_index];
            }
            else{

                // This threshold would select serial execution.
                selected_time_us =
                    one_thread_results[size_index];
            }


            //----------------------------------------------------------------------------------------------------------------
            // FIND THE FASTEST MEASURED VERSION FOR THIS SIZE
            //----------------------------------------------------------------------------------------------------------------

            double fastest_time_us =
                one_thread_results[size_index];

            if (
                parallel_results[size_index] <
                fastest_time_us
            ){
                fastest_time_us =
                    parallel_results[size_index];
            }


            //----------------------------------------------------------------------------------------------------------------
            // CALCULATE HOW MUCH SLOWER THE SELECTED VERSION IS
            //
            // 1.0 means the threshold selected the fastest version.
            // 1.1 means it was 10% slower than the fastest version.
            //----------------------------------------------------------------------------------------------------------------

            score +=
                selected_time_us /
                fastest_time_us;
        }


        //-------------------------------------------------------------------------------------------------------------------
        // CALCULATE THE AVERAGE SCORE
        //-------------------------------------------------------------------------------------------------------------------

        score /=
            static_cast<double>(number_of_sizes);


        //-------------------------------------------------------------------------------------------------------------------
        // PRINT THIS OMP_MIN_WORK RESULT
        //-------------------------------------------------------------------------------------------------------------------

        std::cout
            << omp_min_work
            << ", "
            << std::fixed
            << std::setprecision(4)
            << score
            << "\n";


        //-------------------------------------------------------------------------------------------------------------------
        // SAVE THE BEST OMP_MIN_WORK
        //-------------------------------------------------------------------------------------------------------------------

        if (score < best_score){

            best_score = score;
            best_omp_min_work = omp_min_work;
        }
    }


    //-----------------------------------------------------------------------------------------------------------------------
    // PRINT FINAL ESTIMATE
    //-----------------------------------------------------------------------------------------------------------------------

    std::cout
        << "\n"
        << "Best estimated omp_min_work: "
        << best_omp_min_work
        << "\n";

    std::cout
        << "Average slowdown score: "
        << std::fixed
        << std::setprecision(4)
        << best_score
        << "\n";

    std::cout
        << "\n"
        << "A score of 1.0000 means the threshold selected the\n"
        << "fastest measured version for every tested size.\n";


    //-----------------------------------------------------------------------------------------------------------------------
    // PRINT CHECKSUM
    //-----------------------------------------------------------------------------------------------------------------------

    std::cout
        << "Benchmark checksum: "
        << benchmark_sink
        << "\n";

#endif
}






//---------------------------------------------------------------------------------------------------------------------------
// VARIABLE DESCRIPTION
//---------------------------------------------------------------------------------------------------------------------------
// Variable difinition like:
panic::types::real_t x = 2.4;

panic::tensor::real_vector a(3);
panic::tensor::uint_vector b(3);
panic::tensor::int_vector c(3);

panic::tensor::real_matrix A(2,2, 1);
panic::tensor::uint_matrix B(2,2, 2);
panic::tensor::int_matrix C(2,2, 3);

//---------------------------------------------------------------------------------------------------------------------------
// FUNCTION PROTOTYPE
//---------------------------------------------------------------------------------------------------------------------------
// Function prototypes like:
// static void TestLeds(void);

int main(void) {

    //benchmark_omp_min_work();

    // Comment out benchmark_omp_min_work() when it is not needed.

    std::cout << "neural_network" << std::endl;

    panic::tensor::real_matrix X;
    panic::tensor::uint_vector y;

    panic::types::uint_t samples = 10;
    panic::types::uint_t classes = 3;


    // create spiral data
    panic::neural_network::spiral_data(samples, classes, X, y);

    // Initilise my model
    panic::neural_network::model mymodel;

    // Create Dense layer with 2 input features and 3 output values
    mymodel.add_layer_dense(2,3);
    
    // Create an activation ReLU layer 
    mymodel.add_activation_relu();

    // Create a second dense layer with 3 inputs and 3 outputs
    mymodel.add_layer_dense(3, 3);

    // Create activation softmax layer
    mymodel.add_activation_softmax();

    mymodel.forward(X); 

    panic::io::print_matrix(mymodel.outputs);




	return 0;
}