//#####################################################################
// Copyright 2013.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
// Class IDENTITY_MATRIX
//##################################################################### 
#ifndef __IDENTITY_MATRIX__
#define __IDENTITY_MATRIX__

#include <Core/Math_Tools/cube.h>
#include <Core/Matrices/GENERAL_MATRIX_BASE.h>
#include <Core/Matrices/MATRIX.h>
#include <Core/Matrices/SCALE_MATRIX.h>
#include <Core/Matrices/SYMMETRIC_MATRIX.h>
#include <Core/Vectors/VECTOR.h>
#include <Core/Vectors/ZERO_VECTOR.h>
#include <cmath>
namespace PhysBAM{
template<class T,int d>
class IDENTITY_MATRIX:public GENERAL_MATRIX_BASE<T,IDENTITY_MATRIX<T,d>>
{
public:
    static constexpr int m=d,n=d,fixed_m=m,fixed_n=n;
    static constexpr bool is_sym=true,is_diag=true;
    static constexpr int type=mat_id;
    typedef T SCALAR;
    IDENTITY_MATRIX()=default;

    explicit IDENTITY_MATRIX(INITIAL_SIZE mm,INITIAL_SIZE nn=INITIAL_SIZE(d)){}
    IDENTITY_MATRIX(const IDENTITY_MATRIX&)=default;

    int Rows() const
    {return m;}

    int Columns() const
    {return n;}

    T operator()(const int i,const int j) const
    {assert(i==j);return 1;}

    SCALE_MATRIX<T,d> operator-() const
    {return SCALE_MATRIX<T,d>(-1);}

    SCALE_MATRIX<T,d> operator+(T a) const
    {return SCALE_MATRIX<T,d>(1+a);}

    SCALE_MATRIX<T,d> operator-(T a) const
    {return SCALE_MATRIX<T,d>(1-a);}

    SCALE_MATRIX<T,d> operator*(T a) const
    {return SCALE_MATRIX<T,d>(a);}

    SCALE_MATRIX<T,d> operator/(T a) const
    {return SCALE_MATRIX<T,d>(1/a);}

    IDENTITY_MATRIX Transposed() const
    {return *this;}

    template<class A,class R=std::enable_if_t<(IS_MATRIX<A>::value||IS_VECTOR<A>::value)&&compat(m,A::fixed_m)>>
    A Transpose_Times(const A& a) const
    {return a;}

    template<class A,class R=std::enable_if_t<IS_MATRIX<A>::value&&compat(n,A::fixed_n)>>
    auto Times_Transpose(const A& a) const
    {return a.Transposed();}

    FIXED_NUMBER<T,d> Trace() const
    {return FIXED_NUMBER<T,d>();}

    IDENTITY_MATRIX Symmetric_Part() const
    {return *this;}

    SCALE_MATRIX<T,d> Twice_Symmetric_Part() const
    {return *this*2;}

    IDENTITY_MATRIX Diagonal_Part() const
    {return *this;}

    static IDENTITY_MATRIX Identity_Matrix()
    {return {};}

    T Determinant() const
    {
        if constexpr(m>0) return 1;
        else return T{};
    }

    T Max_Abs() const
    {
        if constexpr(m>0) return 1;
        else return T{};
    }

    IDENTITY_MATRIX Normal_Equations_Matrix() const
    {return *this;}

    IDENTITY_MATRIX Outer_Product_Matrix() const
    {return *this;}

    template<class A>
    auto Double_Contract(const A& a)
    {return ::PhysBAM::Double_Contract(*this,a);}
};
template<int i,int j,class T,int d> auto get(const IDENTITY_MATRIX<T,d>& a)
{
    if constexpr(i==j) return T{1};
    else return T{};
}
template<class T,int d> auto diag(const IDENTITY_MATRIX<T,d>& a,int i){return T{1};}
template<class T,int d> auto lower_(const IDENTITY_MATRIX<T,d>& a,int i,int j){return T{};}
template<class T,int d> auto upper_(const IDENTITY_MATRIX<T,d>& a,int i,int j){return T{};}

template<class T,int d>
inline std::ostream& operator<< (std::ostream& output_stream,const IDENTITY_MATRIX<T,d>& A)
{return output_stream<<SCALE_MATRIX<T,d>(1);}
}
#endif
