Skeleton/src/Algebra/Matrices/NormalMatrices/FUnitaryMatrix.cs

using System.Numerics;
using Skeleton.Algebra.AdditionalProperties;
using Skeleton.Algebra.AdditionalProperties.Extensions;
using Skeleton.Algebra.Extensions;
using Skeleton.Algebra.FixableObjects;

namespace Skeleton.Algebra;

public static partial class CategoryOf<TDimension>
{
    /// <summary>
    /// Matrices whose inverse is its hermitian
    /// AA* = I = A*A 
    /// </summary>
    public class FUnitaryMatrix :
        OnField<Complex>.FNormalMatrix,
        IUnitaryInvariantMatrix<FUnitaryMatrix>
    {
        /// <inheritdoc />
        public FUnitaryMatrix() => ResetOne();

        /// <inheritdoc />
        public FUnitaryMatrix(IEnumerable<OnField<Complex>.FVector> es, bool row = true) : base(es, row) =>
            this.PostConstructionFix();

        /// <inheritdoc />
        public FUnitaryMatrix(params Complex[] es) : base(es) => this.PostConstructionFix();

        /// <inheritdoc />
        public FUnitaryMatrix(OnField<Complex>.FMatrix u, OnField<Complex>.FMatrix d) : base(u, d) =>
            this.PostConstructionFix();

        /// <inheritdoc />
        public FUnitaryMatrix(OnField<Complex>.FNormalMatrix n) : base(n) => 
            this.PostConstructionFix();
        /// <summary>
        /// id matrix as unitary matrix
        /// </summary>
        public new static readonly FUnitaryMatrix One = new();

        /// <summary>
        /// inverse of the matrix
        /// </summary>
        /// <returns></returns>
        public new FUnitaryMatrix Inv() => Hermitian();

        /// <summary>
        /// multiplication between unitary matrices
        /// </summary>
        /// <param name="self"></param>
        /// <param name="other"></param>
        /// <returns></returns>
        public static FUnitaryMatrix operator *(FUnitaryMatrix self, FUnitaryMatrix other) =>
            new(self.AsMatrix * other);

        /// <summary>
        /// transpose + conjugate
        /// </summary>
        /// <returns></returns>
        public new FUnitaryMatrix Hermitian() => new(U, D.Hermitian());

        internal override void InternalFix()
        {
            for (int i = 1; i <= FDim; i++)
                D[i, i] /= D[i, i].Magnitude;
        }

        private bool UnitaryNeedFix => 
            !D.Diagonals.Aggregate((a, b) => a * b).Magnitude.IsEqualApprox(1);
        internal override bool NeedFix => base.NeedFix || UnitaryNeedFix;

        /// <summary>
        ///     a conjugate on b => a b inv(a)
        /// </summary>
        /// <param name="a"></param>
        /// <typeparam name="TUnitaryInvariant"></typeparam>
        /// <returns></returns>
        public TUnitaryInvariant ConjugateOn<TUnitaryInvariant>(TUnitaryInvariant a)
            where TUnitaryInvariant : OnField<Complex>.FMatrix, IUnitaryInvariantMatrix<TUnitaryInvariant>, new()
        {
            TUnitaryInvariant res = new();
            res.Elements = (this * a * Hermitian()).ToArray();
            res.Fix();
            return res;
        }

        /// <summary>
        /// conjugate of real is itself<br/>
        /// so if homomorphism h: R -> K kept the conj<br/>
        /// h(x) = h(conj(x)) = conj(h(x))<br/>
        /// h will map real diag matrix to hermitian matrix in K
        /// </summary>
        /// <param name="a"></param>
        /// <returns></returns>
        public FHermitianMatrix ConjugateOn(OnField<double>.FDiagonalMatrix a)
        {
            FHermitianMatrix g = a.AsComplex().MatrixCast<FHermitianMatrix>();
            return ConjugateOn(g);
        }

        /// <summary>
        /// log of unitary matrix
        /// </summary>
        /// <returns></returns>
        public FLieUnitaryMatrix Log()
        {
            OnField<Complex>.FMatrix logd = new();
            for (int i = 1; i <= FDim; i++)
                logd[i, i] = Complex.Log(D[i, i]);
            return new(U, logd);
        }

        /*
         * /// <inheritdoc />
         * public override JDPair<OnField<Complex>.FMatrix> EigenDecomposition() => this.UnitaryDecomposition();
         * */

        internal override void UnitaryFix(OnField<Complex>.FMatrix a) => a.UnitaryFix();
        protected override bool NeedPreNormalFix() => !Property.IsUnitary;
        protected override void PreNormalConstructionFix()
        {
            if (Property.IsDiagonal)
            {
                for (int i = 1; i <= Dim; i++)
                    this[i, i] /= (this[i, i] * Complex.Conjugate(this[i, i]));
            }
            else
            {
                this[1] = this[1].Normalize;
                for (int i = 2; i <= Dim; i++)
                {
                    for (int j = 1; j < i; j++)
                        this[i] -= this[i].ProjectionOn(this[j]);
                    this[i] = this[i].Normalize;
                }
            }
        }
    }
}