//#####################################################################
// Copyright 2003-2007, Geoffrey Irving, Neil Molino, Craig Schroeder, Tamar Shinar, Jerry Talton.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
// Class DIAGONAL_MATRIX
//#####################################################################
#ifndef __DIAGONAL_MATRIX__
#define __DIAGONAL_MATRIX__
#include <Core/Matrices/GENERAL_MATRIX_BASE.h>
#include <Core/Matrices/MATRIX.h>
#include <Core/Matrices/MATRIX_OPERATIONS.h>
#include <Core/Vectors/SCALAR_POLICY.h>
#include <Core/Vectors/VECTOR.h>
#include <Core/Vectors/ZERO_VECTOR.h>

namespace PhysBAM{

template<class T,int d_input> class DIAGONAL_MATRIX;

template<class T,int d> struct IS_SCALAR_BLOCK<DIAGONAL_MATRIX<T,d> >:public IS_SCALAR_BLOCK<T>{};
template<class T,int d,class RW> struct IS_BINARY_IO_SAFE<DIAGONAL_MATRIX<T,d>,RW>:public IS_BINARY_IO_SAFE<T,RW>{};

template<class T,int d_input>
class DIAGONAL_MATRIX:public GENERAL_MATRIX_BASE<T,DIAGONAL_MATRIX<T,d_input>>
{
public:
    typedef T SCALAR;
    static constexpr int d=d_input,m=d_input,n=d_input;
    static constexpr int fixed_m=m,fixed_n=n;
    static constexpr bool is_sym=true,is_diag=true;
    static constexpr int type=mat_diag;

    VECTOR<T,d> x;

    DIAGONAL_MATRIX()=default;

    explicit DIAGONAL_MATRIX(INITIAL_SIZE mm,INITIAL_SIZE nn=INITIAL_SIZE(d)){}

    DIAGONAL_MATRIX(const DIAGONAL_MATRIX& A)=default;

    template<class T2> explicit
    DIAGONAL_MATRIX(const DIAGONAL_MATRIX<T2,d>& matrix_input)
        :x(matrix_input.x)
    {}

    explicit DIAGONAL_MATRIX(const VECTOR<T,d>& v)
        :x(v)
    {}

    template<class ...Args,class R=std::enable_if_t<(1+sizeof...(Args)==d)>>
    DIAGONAL_MATRIX(const T& a,Args&& ...args)
        :x{a,((T)args)...}
    {
    }

    DIAGONAL_MATRIX(const SCALE_MATRIX<T,d>& A)
    {
        x.Fill(A.x);
    }

    DIAGONAL_MATRIX(const IDENTITY_MATRIX<T,d>& A)
    {
        x.Fill(1);
    }

    DIAGONAL_MATRIX(const ZERO_MATRIX<T,d>& A)
    {
    }

    int Rows() const
    {return d;}

    int Columns() const
    {return d;}

    T& operator()(const int i)
    {return x(i);}

    const T& operator()(const int i) const
    {return x(i);}

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

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

    bool Valid_Index(const int i,const int j) const
    {return (unsigned)i<d && i==j;}

    T Determinant() const
    {
        if constexpr(d>0) return x.Product();
        else return T{};
    }

    DIAGONAL_MATRIX Inverse() const
    {return DIAGONAL_MATRIX(PhysBAM::Inverse(x));}

    VECTOR<T,d> Inverse_Times(const VECTOR<T,d>& v) const
    {return v/x;}

    DIAGONAL_MATRIX Cofactor_Matrix() const
    {
        if constexpr(d==0) return {};
        else if constexpr(d==1) return {1};
        else
            return []<int...k>(auto& u,std::integer_sequence<int,k...>)
                {
                    T a=u(0),b=u(d-1);
                    T l[]={a,(a*=u(k+1))...};
                    T r[]={b,(b*=u(d-2-k))...};
                    return DIAGONAL_MATRIX{r[d-2],(l[k]*r[d-3-k])...,l[d-2]};
                }(x,std::make_integer_sequence<int,d-2>{});
    }

    DIAGONAL_MATRIX Positive_Definite_Part() const
    {return DIAGONAL_MATRIX(Componentwise_Max(x,(T)0));}

    VECTOR<T,d> Vector() const
    {return x;}

    template<class RW> void Read(UNTYPED_ISTREAM<RW>& input)
    {Read_Binary(input,x);}

    template<class RW> void Write(UNTYPED_OSTREAM<RW>& output) const
    {Write_Binary(output,x);}

    void Fill_Uniform(RANDOM_NUMBERS<T>& rand,const T a,const T b)
    {rand.Fill_Uniform(x,a,b);}
//#####################################################################
};
// global functions

template<int i,int j,class T,int d> auto get(const DIAGONAL_MATRIX<T,d>& a)
{
    if constexpr(i==j) return a.x(i);
    else return T{};
}
template<class T,int d> auto diag(const DIAGONAL_MATRIX<T,d>& a,int i){return a.x(i);}
template<class T,int d> auto lower_(const DIAGONAL_MATRIX<T,d>& a,int i,int j){return T{};}
template<class T,int d> auto upper_(const DIAGONAL_MATRIX<T,d>& a,int i,int j){return T{};}

template<class T,int d> DIAGONAL_MATRIX<T,d>
Transpose_Times_Self(const DIAGONAL_MATRIX<T,d>& a)
{return a*a;}

template<class T,int d>
inline std::istream& operator>>(std::istream& input_stream,DIAGONAL_MATRIX<T,d>& A)
{return input_stream>>A.x;}

template<class T,int d>
inline std::ostream& operator<<(std::ostream& output_stream,const DIAGONAL_MATRIX<T,d>& A)
{return output_stream<<A.x;}

template<class T,int d>
inline DIAGONAL_MATRIX<T,d> log(const DIAGONAL_MATRIX<T,d>& A)
{return DIAGONAL_MATRIX<T,d>(log(A.x));}

template<class T,int d>
inline DIAGONAL_MATRIX<T,d> exp(const DIAGONAL_MATRIX<T,d>& A)
{return DIAGONAL_MATRIX<T,d>(exp(A.x));}

template<class T,int d>
inline DIAGONAL_MATRIX<T,d> abs(const DIAGONAL_MATRIX<T,d>& A)
{return DIAGONAL_MATRIX<T,d>(abs(A.x));}
}
#endif
