//#####################################################################
// Copyright 2022, Craig Schroeder.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
#include <Core/Symbolics/SYMBOLIC_SPARSE_MATRIX_HELPER.h>
#include <Core/Data_Structures/HASHTABLE.h>
#include <Core/Data_Structures/TUPLE.h>
#include <Core/Krylov_Solvers/KRYLOV_SYSTEM_BASE.h>
#include <Core/Krylov_Solvers/KRYLOV_VECTOR_BASE.h>
#include <Core/Math_Tools/cube.h>
#include <Core/Math_Tools/sqr.h>
#include <Core/Read_Write/OCTAVE_OUTPUT.h>
#include <Core/Symbolics/SYMBOLIC.h>
#include <Core/Utilities/DEBUG_CAST.h>
#include <Core/Vectors/SCALAR_POLICY.h>
#include <cmath>
#include <complex>
#include <iosfwd>
#include <string>
namespace PhysBAM{

void SYMBOLIC_SPARSE_MATRIX_HELPER::Set(double x)
{
    coeffs.Clean_Memory();
}

void SYMBOLIC_SPARSE_MATRIX_HELPER::Print(double x,std::ostream& o) const
{
    o<<coeffs;
}

SYMBOLIC_BASE* SYMBOLIC_SPARSE_MATRIX_HELPER::Clone() const
{
    return new SYMBOLIC_SPARSE_MATRIX_HELPER;
}

void SYMBOLIC_SPARSE_MATRIX_HELPER::Binary(double x,const SYMBOLIC& a,const SYMBOLIC& b,type_2 type)
{
    ARRAY<std::pair<int,double>> empty;

    const auto* A=dynamic_cast<const SYMBOLIC_SPARSE_MATRIX_HELPER*>(a.z);
    const auto* B=dynamic_cast<const SYMBOLIC_SPARSE_MATRIX_HELPER*>(b.z);
    if(!A && !B) return;
    const auto& Ac=A?A->coeffs:empty;
    const auto& Bc=B?B->coeffs:empty;
    switch(type)
    {
        case type_2::add:
            if(Ac.m || Bc.m)
            {
                int i=0,j=0;
                while(i<Ac.m && j<Bc.m)
                {
                    if(Ac(i).first<Bc(j).first) coeffs.Append(Ac(i++));
                    else if(Ac(i).first>Bc(j).first) coeffs.Append(Bc(j++));
                    else
                    {
                        coeffs.Append({Ac(i).first,Ac(i).second+Bc(j++).second});
                        i++;
                    }
                }
                while(i<Ac.m) coeffs.Append(Ac(i++));
                while(j<Bc.m) coeffs.Append(Bc(j++));
            }
            break;

        case type_2::sub:
            if(Ac.m || Bc.m)
            {
                int i=0,j=0;
                while(i<Ac.m && j<Bc.m)
                {
                    if(Ac(i).first<Bc(j).first) coeffs.Append(Ac(i++));
                    else if(Ac(i).first>Bc(j).first)
                    {
                        coeffs.Append(Bc(j++));
                        coeffs.Last().second*=-1;
                    }
                    else
                    {
                        coeffs.Append({Ac(i).first,Ac(i).second-Bc(j++).second});
                        i++;
                    }
                }
                while(i<Ac.m) coeffs.Append(Ac(i++));
                while(j<Bc.m)
                {
                    coeffs.Append(Bc(j++));
                    coeffs.Last().second*=-1;
                }
            }
            break;

        case type_2::mul:
            if(Ac.m)
            {
                assert(!Bc.m);
                coeffs=Ac;
                for(auto&i:coeffs) i.second*=b.x;
            }
            else if(Bc.m)
            {
                coeffs=Bc;
                for(auto&i:coeffs) i.second*=a.x;
            }
            break;

        case type_2::div:
            assert(!Bc.m);
            if(Ac.m)
            {
                coeffs=Ac;
                for(auto&i:coeffs) i.second/=b.x;
            }
            break;

        default:
            assert(!Ac.m && !Bc.m);
            break;
    }

    int k=0;
    for(int i=0;i<coeffs.m;i++)
        if(coeffs(i).second)
            coeffs(k++)=coeffs(i);
    coeffs.Resize(k);
}

void SYMBOLIC_SPARSE_MATRIX_HELPER::Unary(double x,const SYMBOLIC& a,type_1 type)
{
    const auto* A=dynamic_cast<const SYMBOLIC_SPARSE_MATRIX_HELPER*>(a.z);
    if(!A) return;
    const auto& Ac=A->coeffs;
    switch(type)
    {
        case type_1::neg:
            if(Ac.m)
            {
                coeffs=Ac;
                for(auto&a:coeffs) a.second=-a.second;
            }
            break;

        case type_1::copy:
        case type_1::pos:
            coeffs=Ac;
            break;

        default:
            assert(!Ac.m);
            break;
    }
}
}
