#ifdef _MSC_VER  //If the compiler is the Visual Studio Compiler
#define NOMINMAX //solves a problem with std::max on Windows
#pragma warning(disable: 4297) //remove throw warnings
#endif

#include "common/definitions.h"
#include <iostream>
#include <vector>

#include "common/SpinSystem.h"

const bool debug = false;

const char* get_version()
{
    return SPINTRUM_VERSION;
}

void print_software_info()
{
    std::cout<<"Spintrum version "<<get_version()<<". ";
    std::cout<<"This software is provided under a license that should have been delivered with its package in the file LICENSE.txt."
               "For more information, contact Samer Afach: samer [[at]] afach.de."<<std::endl<<std::endl;
}

/**
 * A quantum mechanics operation (populate spins, tip spins, etc...) is a struct that contains a specific, predefined number of parameters.
 * Each operation has an Identifier, which is simply an integer. From this identifier, the code knows what parameters to seek
 * and how many parameters should be expected.
 * There are integer parameters, real parameters, and array of real parameters.
 */
struct SpinOperation
{
    long     OperationIdentifier;
    long     Params_int_count;     //number of integer parameters
    long*    Params_int;           //integer parameters
    long     Params_real_count;    //number of real parameters
    double*  Params_real;          //real parameters
    long     Params_array_count;   //number of arrays
    long*    Params_array_lengths; //lengths of suppled arrays
    char**   Params_array_names;   //names of the parameters of the arrays
    double** Params_array;         //array data
};
const static int OPERATION__THERMAL_POPULATE = 0;
const static int OPERATION__TIP_SPINS = 1;
const static int OPERATION__SET_HAMILTONIAN = 2;
const static int OPERATION__EVOLVE_TIME_INDEPENDENT = 3;
const static int OPERATION__INIT_TIME_INDEPENDENT_EVOLUTION = 4;
const static int OPERATION__EVOLVE_TIME_DEPENDENT = 5;

typedef double Real;

typedef SpinSystem<Real,std::complex<Real>,VEC,CX_VEC,MAT,CX_MAT> SpinSystemType;

struct SpinInfo
{
    CX_MAT jCouplings;
    CX_VEC gammas;
    G_MAT<int> multiplicities;
};

void devectorizeJCouplings(const long numOfSpins, const double* const jCouplingsInput, CX_MAT& jCouplingsOutput)
{
    if(jCouplingsInput == 0)
    {
        jCouplingsOutput.resize(0,0);
        return;
    }
    int i = 0; //element of given array
    int j = 0; //element of current row
    int row = 0;
    while (i < numOfSpins*(numOfSpins-1)/2)
    {
        if(j <= row)
        {
            j++;
            continue;
        }
        if(j >= numOfSpins)
        {
            j = 0;
            row++;
            continue;
        }
        jCouplingsOutput.at(row,j) = jCouplingsInput[i];
        j++;
        i++;
    }
}
extern "C" int InitializeSimulator(double* gammas, double* jCouplings, int* spinMultiplicities, int numOfSpins, bool doPrintOutput, void** SpinSystemObject)
{
    try
    {
        print_software_info();
        //    std::cout<<"Number of spins: "<<numOfSpins<<std::endl;
        SpinInfo spinInfo;
        spinInfo.gammas.resize(numOfSpins,1);
        spinInfo.multiplicities.resize(numOfSpins,1);
        spinInfo.jCouplings.resize(numOfSpins,numOfSpins);
        for(int i = 0; i < numOfSpins; i++)
        {
            //        std::cout<<gammas[i]<<"\t";
            spinInfo.gammas[i] = gammas[i];
        }
        //    std::cout<<spinInfo.gammas<<std::endl;
        if(spinMultiplicities == 0)
        {
            for(int i = 0; i < numOfSpins; i++)
                spinInfo.multiplicities[i] = 2;
        }
        else
        {
            for(int i = 0; i < numOfSpins; i++)
                spinInfo.multiplicities[i] = spinMultiplicities[i];
        }

        devectorizeJCouplings(numOfSpins, jCouplings, spinInfo.jCouplings);
        //    std::cout<<spinInfo.jCouplings<<std::endl;

        SpinSystemType *Spins = new SpinSystem<Real,std::complex<Real>,VEC,CX_VEC,MAT,CX_MAT>();
        Spins->print = doPrintOutput;
        Spins->init(spinInfo.gammas,spinInfo.jCouplings,spinInfo.multiplicities);
        *SpinSystemObject = static_cast<void*>(Spins);
        return 0;
    }
    catch(std::exception &ex)
    {
        std::cout<<"An exception was thrown while trying to initialize the spin system. Exception message is: " << ex.what() << std::endl;
        return -1;
    }
    catch(...)
    {
        std::cout<<"An unhandled exception was thrown while trying to initialize the spin system."<< std::endl;
        return -2;
    }
}

extern "C" int DenitializeSimulator(void* SpinSystemObject)
{
    try
    {
        delete (static_cast<SpinSystemType*>(SpinSystemObject));
        return 0;
    }
    catch(std::exception &ex)
    {
        std::cout<<"An exception was thrown while trying to delete the spin system from memory. Exception message is: " << ex.what() << std::endl;
        return -1;
    }
    catch(...)
    {
        std::cout<<"An unhandled exception was thrown while trying to delete the spin system from memory."<< std::endl;
        return -2;
    }
}


extern "C" int SetGammasInSimulator(double* gammas, void* SpinSystemObject)
{
    try
    {
        CX_VEC gammasVec((static_cast<SpinSystemType*>(SpinSystemObject))->getNumOfSpins(),1);
        for(int i = 0; i < gammasVec.size(); i++)
        {
            gammasVec[i] = gammas[i];
        }
        (static_cast<SpinSystemType*>(SpinSystemObject))->setGammas(gammasVec);
        return 0;
    }
    catch(std::exception &ex)
    {
        std::cout<<"An exception was thrown while trying to delete the spin system from memory. Exception message is: " << ex.what() << std::endl;
        return -1;
    }
    catch(...)
    {
        std::cout<<"An unhandled exception was thrown while trying to delete the spin system from memory."<< std::endl;
        return -2;
    }
}

extern "C" int SetJCouplingsInSimulator(double* jCouplings, void* SpinSystemObject)
{
    try
    {
        CX_MAT jCouplingsMat((static_cast<SpinSystemType*>(SpinSystemObject))->getNumOfSpins(),(static_cast<SpinSystemType*>(SpinSystemObject))->getNumOfSpins());
        devectorizeJCouplings((static_cast<SpinSystemType*>(SpinSystemObject))->getNumOfSpins(), jCouplings, jCouplingsMat);
        (static_cast<SpinSystemType*>(SpinSystemObject))->setJCouplings(jCouplingsMat);
        return 0;
    }
    catch(std::exception &ex)
    {
        std::cout<<"An exception was thrown while trying to delete the spin system from memory. Exception message is: " << ex.what() << std::endl;
        return -1;
    }
    catch(...)
    {
        std::cout<<"An unhandled exception was thrown while trying to delete the spin system from memory."<< std::endl;
        return -2;
    }
}
void _do_operation_evolve_time_dependent(SpinOperation* op, SpinSystemType* SpinsPtr, double* outputArray, long &currentDataPoint)
{
    if(debug) std::cout<<"Printing received data:"<<std::endl;
    if(op->Params_array_count == 0)
    {
        throw std::runtime_error("There's no data provided for the time-dependent evolution.");
    }
    else
    {
        //make sure that all lengths are equal
        long initialLength = op->Params_array_lengths[0];
        for(long i = 1; i < op->Params_array_count; i++)
        {
            if(op->Params_array_lengths[i] != initialLength)
            {
                throw std::length_error("Data arrays in time-dependent evolution are expected all to be of the same length.");
            }
        }
    }
    std::vector<std::string> varsNames(static_cast<std::vector<std::string>::size_type>(op->Params_array_count));
    std::copy(op->Params_array_names,op->Params_array_names+op->Params_array_count,varsNames.begin()); //copy names to vector<string>

    CX_MAT evData(op->Params_array_lengths[0],op->Params_array_count);
    for(long i = 0; i < op->Params_array_count; i++)
    {
        std::copy(&op->Params_array[i][0],&op->Params_array[i][op->Params_array_lengths[i]],evData.begin()+evData.rows()*i);
    }
    std::set<char> dirString;
    if(op->Params_int[1]) dirString.insert('x');
    if(op->Params_int[2]) dirString.insert('y');
    if(op->Params_int[3]) dirString.insert('z');
    std::vector<Real> tr = SpinsPtr->evolveDensityMatrix(op->Params_real[0],varsNames,evData,static_cast<int>(op->Params_int[0]),dirString);
    std::copy(std::make_move_iterator(tr.begin()), std::make_move_iterator(tr.end()), outputArray + currentDataPoint);
    currentDataPoint += op->Params_array_lengths[0]*static_cast<long>(dirString.size());
}

long ExecuteOperations(void* SpinSystemObject, void** ops, int numOfOperations, double* output, bool debug)
{
    if(debug)
    {
        std::cout<<"Number of operations: "<<numOfOperations<<std::endl;
    }

    SpinSystemType* SpinsPtr = static_cast<SpinSystemType*>(SpinSystemObject);
    long currentDataPoint = 0; //for a simulation with many operations that consecutively fill the output, this keeps track of the current point
    for(int i = 0; i < numOfOperations; i++) {
        void** currOperationBytePosition = (static_cast<void**>(ops));
        long OperationIdentifier = *(static_cast<long*>(currOperationBytePosition[i]));
        if(debug) std::cout<<"Identifier: "<<OperationIdentifier<<std::endl;
        if(OperationIdentifier == OPERATION__THERMAL_POPULATE) {
            SpinOperation* op;
            op = static_cast<SpinOperation*>(currOperationBytePosition[i]);
            SpinsPtr->populate_thermal(op->Params_real[0],op->Params_real[1],op->Params_real[2],op->Params_real[3]);
        }
        else if (OperationIdentifier == OPERATION__TIP_SPINS) {
            SpinOperation* op;
            op = static_cast<SpinOperation*>(currOperationBytePosition[i]);
            SpinsPtr->tipSpins(static_cast<int>(op->Params_int[0]),op->Params_real[0]);
        }
        else if (OperationIdentifier == OPERATION__SET_HAMILTONIAN) {
            SpinOperation* op;
            op = static_cast<SpinOperation*>(currOperationBytePosition[i]);
            SpinsPtr->calculateHamiltonians(op->Params_real[0],op->Params_real[1],op->Params_real[2]);
        }
        else if (OperationIdentifier == OPERATION__EVOLVE_TIME_INDEPENDENT) {
            SpinOperation* op;
            op = static_cast<SpinOperation*>(currOperationBytePosition[i]);
            std::vector<Real> trVals = SpinsPtr->solveMagnetization_parallel(op->Params_int[0],static_cast<int>(op->Params_int[1]));
            std::copy(trVals.begin(), trVals.end(), output + currentDataPoint);
            currentDataPoint += op->Params_int[0];
        }
        else if (OperationIdentifier == OPERATION__INIT_TIME_INDEPENDENT_EVOLUTION) {
            SpinOperation* op;
            op = static_cast<SpinOperation*>(currOperationBytePosition[i]);
            SpinsPtr->initTimeIndependentEvolution(op->Params_real[0],static_cast<int>(op->Params_int[0]));
        }
        else if (OperationIdentifier == OPERATION__EVOLVE_TIME_DEPENDENT) {
            _do_operation_evolve_time_dependent(static_cast<SpinOperation*>(currOperationBytePosition[i]), SpinsPtr, output, currentDataPoint);
        }
        else {
            std::cout << "Error: An unknown operation was requested... Exiting." << std::endl;
            exit(1);
        }
    }
    return currentDataPoint;
}

extern "C" long CreateSignalFromInitializedSimulator(void* SpinSystemObject, void** ops, int numOfOperations, long numOfPoints, double* output)
{
    try {
        if(debug)
        {
            std::cout<<"Number of operations: "<<numOfOperations<<std::endl;
            std::cout<<"Number of points: "<<numOfPoints<<std::endl;
        }

        long processedNumOfPoints = ExecuteOperations(SpinSystemObject,ops,numOfOperations,output,debug);
        if (debug) std::cout<<processedNumOfPoints<<std::endl;
        return processedNumOfPoints;
    }
    catch(std::exception& ex) {
        std::cerr << "An exception was thrown by Spintrum. It says: " << ex.what() << std::endl;
        return 0;
    }
    catch(...) {
        std::cerr << "An unknown exception was thrown by Spintrum. " << std::endl;
        return 0;
    }
}
void _fill_spin_info(int numOfSpins, double* gammas, double* jCouplings, int* spinMultiplicities, SpinInfo& spinInfo)
{
    spinInfo.gammas.resize(numOfSpins,1);
    spinInfo.multiplicities.resize(numOfSpins,1);
    spinInfo.jCouplings.resize(numOfSpins,numOfSpins);
    for(int i = 0; i < numOfSpins; i++) {
        spinInfo.gammas[i] = gammas[i];
    }
    if(spinMultiplicities == 0) {
        for(int i = 0; i < numOfSpins; i++)
            spinInfo.multiplicities[i] = 2;
    }
    else {
        for(int i = 0; i < numOfSpins; i++)
            spinInfo.multiplicities[i] = spinMultiplicities[i];
    }
    if(debug) std::cout<<"Gammas:"<<std::endl;
    if(debug) std::cout<<spinInfo.gammas<<std::endl;
    devectorizeJCouplings(numOfSpins, jCouplings, spinInfo.jCouplings);
    if(debug) std::cout<<"J-couplings:"<<std::endl;
    if(debug) std::cout<<spinInfo.jCouplings<<std::endl;
}

extern "C" long SimulateSignal(void** ops, double* gammas, double* jCouplings, int* spinMultiplicities, int numOfSpins, int numOfOperations, long numOfPoints, double* output)
{
    try {
        print_software_info();
        if(debug) {
            std::cout<<"Number of operations: "<<numOfOperations<<std::endl;
            std::cout<<"Number of points: "<<numOfPoints<<std::endl;
            std::cout<<"Number of spins: "<<numOfSpins<<std::endl;
        }
        SpinInfo spinInfo;
        _fill_spin_info(numOfSpins, gammas, jCouplings, spinMultiplicities, spinInfo);
        SpinSystem<Real,std::complex<Real>,VEC,CX_VEC,MAT,CX_MAT > Spins(spinInfo.gammas,spinInfo.jCouplings,spinInfo.multiplicities);
        long processedNumOfPoints = ExecuteOperations(static_cast<void*>(&Spins),ops,numOfOperations,output,debug);
        if (debug) std::cout<<processedNumOfPoints<<std::endl;
        return processedNumOfPoints;
    }
    catch(std::exception& ex) {
        std::cerr << "An exception was thrown by Spintrum. It says: " << ex.what() << std::endl;
        return 0;
    }
    catch(...) {
        std::cerr << "An unknown exception was thrown by Spintrum. " << std::endl;
        return 0;
    }
}