// <copyright file="ManagedFourierTransformProvider.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// https://numerics.mathdotnet.com
//
// Copyright (c) 2009-2018 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
// files (the "Software"), to deal in the Software without
// restriction, including without limitation the rights to use,
// copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following
// conditions:
//
// The above copyright notice and this permission notice shall be
// included in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
// OTHER DEALINGS IN THE SOFTWARE.
// </copyright>

using System;

using Complex = System.Numerics.Complex;

namespace IStation.Numerics.Providers.FourierTransform.Managed
{
    internal partial class ManagedFourierTransformProvider : IFourierTransformProvider
    {
        /// <summary>
        /// Try to find out whether the provider is available, at least in principle.
        /// Verification may still fail if available, but it will certainly fail if unavailable.
        /// </summary>
        public bool IsAvailable()
        {
            return true;
        }

        /// <summary>
        /// Initialize and verify that the provided is indeed available. If not, fall back to alternatives like the managed provider
        /// </summary>
        public void InitializeVerify()
        {
        }

        /// <summary>
        /// Frees memory buffers, caches and handles allocated in or to the provider.
        /// Does not unload the provider itself, it is still usable afterwards.
        /// </summary>
        public virtual void FreeResources()
        {
        }

        public override string ToString()
        {
            return "Managed";
        }

        public void Forward(Complex32[] samples, FourierTransformScaling scaling)
        {
            if (samples.Length.IsPowerOfTwo())
            {
                if (samples.Length >= 1024)
                {
                    Radix2ForwardParallel(samples);
                }
                else
                {
                    Radix2Forward(samples);
                }
            }
            else
            {
                BluesteinForward(samples);
            }

            switch (scaling)
            {
                case FourierTransformScaling.SymmetricScaling:
                {
                    HalfRescale(samples);
                    break;
                }
                case FourierTransformScaling.ForwardScaling:
                {
                    FullRescale(samples);
                    break;
                }
            }
        }

        public void Forward(Complex[] samples, FourierTransformScaling scaling)
        {
            if (samples.Length.IsPowerOfTwo())
            {
                if (samples.Length >= 1024)
                {
                    Radix2ForwardParallel(samples);
                }
                else
                {
                    Radix2Forward(samples);
                }
            }
            else
            {
                BluesteinForward(samples);
            }

            switch (scaling)
            {
                case FourierTransformScaling.SymmetricScaling:
                {
                    HalfRescale(samples);
                    break;
                }
                case FourierTransformScaling.ForwardScaling:
                {
                    FullRescale(samples);
                    break;
                }
            }
        }

        public void Backward(Complex32[] spectrum, FourierTransformScaling scaling)
        {
            if (spectrum.Length.IsPowerOfTwo())
            {
                if (spectrum.Length >= 1024)
                {
                    Radix2InverseParallel(spectrum);
                }
                else
                {
                    Radix2Inverse(spectrum);
                }
            }
            else
            {
                BluesteinInverse(spectrum);
            }

            switch (scaling)
            {
                case FourierTransformScaling.SymmetricScaling:
                {
                    HalfRescale(spectrum);
                    break;
                }
                case FourierTransformScaling.BackwardScaling:
                {
                    FullRescale(spectrum);
                    break;
                }
            }
        }

        public void Backward(Complex[] spectrum, FourierTransformScaling scaling)
        {
            if (spectrum.Length.IsPowerOfTwo())
            {
                if (spectrum.Length >= 1024)
                {
                    Radix2InverseParallel(spectrum);
                }
                else
                {
                    Radix2Inverse(spectrum);
                }
            }
            else
            {
                BluesteinInverse(spectrum);
            }

            switch (scaling)
            {
                case FourierTransformScaling.SymmetricScaling:
                {
                    HalfRescale(spectrum);
                    break;
                }
                case FourierTransformScaling.BackwardScaling:
                {
                    FullRescale(spectrum);
                    break;
                }
            }
        }

        public void ForwardReal(float[] samples, int n, FourierTransformScaling scaling)
        {
            // TODO: backport proper, optimized implementation from Iridium

            Complex32[] data = new Complex32[n];
            for (int i = 0; i < data.Length; i++)
            {
                data[i] = new Complex32(samples[i], 0.0f);
            }

            Forward(data, scaling);

            samples[0] = data[0].Real;
            samples[1] = 0f;
            for (int i = 1, j = 2; i < data.Length / 2; i++)
            {
                samples[j++] = data[i].Real;
                samples[j++] = data[i].Imaginary;
            }
            if (n.IsEven())
            {
                samples[n] = data[data.Length / 2].Real;
                samples[n + 1] = 0f;
            }
            else
            {
                samples[n - 1] = data[data.Length / 2].Real;
                samples[n] = data[data.Length / 2].Imaginary;
            }
        }

        public void ForwardReal(double[] samples, int n, FourierTransformScaling scaling)
        {
            // TODO: backport proper, optimized implementation from Iridium

            Complex[] data = new Complex[n];
            for (int i = 0; i < data.Length; i++)
            {
                data[i] = new Complex(samples[i], 0.0);
            }

            Forward(data, scaling);

            samples[0] = data[0].Real;
            samples[1] = 0d;
            for (int i = 1, j = 2; i < data.Length/2; i++)
            {
                samples[j++] = data[i].Real;
                samples[j++] = data[i].Imaginary;
            }
            if (n.IsEven())
            {
                samples[n] = data[data.Length/2].Real;
                samples[n+1] = 0d;
            }
            else
            {
                samples[n-1] = data[data.Length / 2].Real;
                samples[n] = data[data.Length / 2].Imaginary;
            }
        }

        public void BackwardReal(float[] spectrum, int n, FourierTransformScaling scaling)
        {
            // TODO: backport proper, optimized implementation from Iridium

            Complex32[] data = new Complex32[n];
            data[0] = new Complex32(spectrum[0], 0f);
            for (int i = 1, j = 2; i < data.Length / 2; i++)
            {
                data[i] = new Complex32(spectrum[j++], spectrum[j++]);
                data[data.Length - i] = data[i].Conjugate();
            }
            if (n.IsEven())
            {
                data[data.Length / 2] = new Complex32(spectrum[n], 0f);
            }
            else
            {
                data[data.Length / 2] = new Complex32(spectrum[n - 1], spectrum[n]);
                data[data.Length / 2 + 1] = data[data.Length / 2].Conjugate();
            }

            Backward(data, scaling);

            for (int i = 0; i < data.Length; i++)
            {
                spectrum[i] = data[i].Real;
            }
            spectrum[n] = 0f;
        }

        public void BackwardReal(double[] spectrum, int n, FourierTransformScaling scaling)
        {
            // TODO: backport proper, optimized implementation from Iridium

            Complex[] data = new Complex[n];
            data[0] = new Complex(spectrum[0], 0d);
            for (int i = 1, j = 2; i < data.Length/2; i++)
            {
                data[i] = new Complex(spectrum[j++], spectrum[j++]);
                data[data.Length - i] = data[i].Conjugate();
            }
            if (n.IsEven())
            {
                data[data.Length/2] = new Complex(spectrum[n], 0d);
            }
            else
            {
                data[data.Length/2] = new Complex(spectrum[n-1], spectrum[n]);
                data[data.Length/2 + 1] = data[data.Length/2].Conjugate();
            }

            Backward(data, scaling);

            for (int i = 0; i < data.Length; i++)
            {
                spectrum[i] = data[i].Real;
            }
            spectrum[n] = 0d;
        }

        public void ForwardMultidim(Complex32[] samples, int[] dimensions, FourierTransformScaling scaling)
        {
            throw new NotSupportedException();
        }

        public void ForwardMultidim(Complex[] samples, int[] dimensions, FourierTransformScaling scaling)
        {
            throw new NotSupportedException();
        }

        public void BackwardMultidim(Complex32[] spectrum, int[] dimensions, FourierTransformScaling scaling)
        {
            throw new NotSupportedException();
        }

        public void BackwardMultidim(Complex[] spectrum, int[] dimensions, FourierTransformScaling scaling)
        {
            throw new NotSupportedException();
        }

        /// <summary>
        /// Fully rescale the FFT result.
        /// </summary>
        /// <param name="samples">Sample Vector.</param>
        private static void FullRescale(Complex32[] samples)
        {
            var scalingFactor = (float)1.0 / samples.Length;
            for (int i = 0; i < samples.Length; i++)
            {
                samples[i] *= scalingFactor;
            }
        }

        /// <summary>
        /// Fully rescale the FFT result.
        /// </summary>
        /// <param name="samples">Sample Vector.</param>
        private static void FullRescale(Complex[] samples)
        {
            var scalingFactor = 1.0 / samples.Length;
            for (int i = 0; i < samples.Length; i++)
            {
                samples[i] *= scalingFactor;
            }
        }

        /// <summary>
        /// Half rescale the FFT result (e.g. for symmetric transforms).
        /// </summary>
        /// <param name="samples">Sample Vector.</param>
        private static void HalfRescale(Complex32[] samples)
        {
            var scalingFactor = (float)Math.Sqrt(1.0 / samples.Length);
            for (int i = 0; i < samples.Length; i++)
            {
                samples[i] *= scalingFactor;
            }
        }

        /// <summary>
        /// Fully rescale the FFT result (e.g. for symmetric transforms).
        /// </summary>
        /// <param name="samples">Sample Vector.</param>
        private static void HalfRescale(Complex[] samples)
        {
            var scalingFactor = Math.Sqrt(1.0 / samples.Length);
            for (int i = 0; i < samples.Length; i++)
            {
                samples[i] *= scalingFactor;
            }
        }
    }
}