// <copyright file="Generate.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://numerics.mathdotnet.com
// http://github.com/mathnet/mathnet-numerics
//
// Copyright (c) 2009-2016 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 System.Collections.Generic;
using System.Linq;
using IStation.Numerics.Distributions;
using IStation.Numerics.Random;
using IStation.Numerics.Threading;
using BigInteger = System.Numerics.BigInteger;
using Complex = System.Numerics.Complex;

namespace IStation.Numerics
{
    public static class Generate
    {
        /// <summary>
        /// Generate samples by sampling a function at the provided points.
        /// </summary>
        public static T[] Map<TA, T>(TA[] points, Func<TA, T> map)
        {
            var res = new T[points.Length];
            for (int i = 0; i < points.Length; i++)
            {
                res[i] = map(points[i]);
            }
            return res;
        }

        /// <summary>
        /// Generate a sample sequence by sampling a function at the provided point sequence.
        /// </summary>
        public static IEnumerable<T> MapSequence<TA, T>(IEnumerable<TA> points, Func<TA, T> map)
        {
            return points.Select(map);
        }

        /// <summary>
        /// Generate samples by sampling a function at the provided points.
        /// </summary>
        public static T[] Map2<TA, TB, T>(TA[] pointsA, TB[] pointsB, Func<TA, TB, T> map)
        {
            if (pointsA.Length != pointsB.Length)
            {
                throw new ArgumentException("The array arguments must have the same length.", nameof(pointsB));
            }

            var res = new T[pointsA.Length];
            for (int i = 0; i < res.Length; i++)
            {
                res[i] = map(pointsA[i], pointsB[i]);
            }
            return res;
        }

        /// <summary>
        /// Generate a sample sequence by sampling a function at the provided point sequence.
        /// </summary>
        public static IEnumerable<T> Map2Sequence<TA, TB, T>(IEnumerable<TA> pointsA, IEnumerable<TB> pointsB, Func<TA, TB, T> map)
        {
            return pointsA.Zip(pointsB, map);
        }

        /// <summary>
        /// Generate a linearly spaced sample vector of the given length between the specified values (inclusive).
        /// Equivalent to MATLAB linspace but with the length as first instead of last argument.
        /// </summary>
        public static double[] LinearSpaced(int length, double start, double stop)
        {
            if (length < 0)
            {
                throw new ArgumentOutOfRangeException(nameof(length));
            }

            if (length == 0) return new double[0];
            if (length == 1) return new[] { stop };

            double step = (stop - start)/(length - 1);

            var data = new double[length];
            for (int i = 0; i < data.Length; i++)
            {
                data[i] = start + i*step;
            }
            data[data.Length - 1] = stop;
            return data;
        }

        /// <summary>
        /// Generate samples by sampling a function at linearly spaced points between the specified values (inclusive).
        /// </summary>
        public static T[] LinearSpacedMap<T>(int length, double start, double stop, Func<double, T> map)
        {
            if (length < 0)
            {
                throw new ArgumentOutOfRangeException(nameof(length));
            }

            if (length == 0) return new T[0];
            if (length == 1) return new[] { map(stop) };

            double step = (stop - start)/(length - 1);

            var data = new T[length];
            for (int i = 0; i < data.Length; i++)
            {
                data[i] = map(start + i*step);
            }
            data[data.Length - 1] = map(stop);
            return data;
        }

        /// <summary>
        /// Generate a base 10 logarithmically spaced sample vector of the given length between the specified decade exponents (inclusive).
        /// Equivalent to MATLAB logspace but with the length as first instead of last argument.
        /// </summary>
        public static double[] LogSpaced(int length, double startExponent, double stopExponent)
        {
            if (length < 0)
            {
                throw new ArgumentOutOfRangeException(nameof(length));
            }

            if (length == 0) return new double[0];
            if (length == 1) return new[] { Math.Pow(10, stopExponent) };

            double step = (stopExponent - startExponent)/(length - 1);

            var data = new double[length];
            for (int i = 0; i < data.Length; i++)
            {
                data[i] = Math.Pow(10, startExponent + i*step);
            }
            data[data.Length - 1] = Math.Pow(10, stopExponent);
            return data;
        }

        /// <summary>
        /// Generate samples by sampling a function at base 10 logarithmically spaced points between the specified decade exponents (inclusive).
        /// </summary>
        public static T[] LogSpacedMap<T>(int length, double startExponent, double stopExponent, Func<double, T> map)
        {
            if (length < 0)
            {
                throw new ArgumentOutOfRangeException(nameof(length));
            }

            if (length == 0) return new T[0];
            if (length == 1) return new[] { map(Math.Pow(10, stopExponent)) };

            double step = (stopExponent - startExponent)/(length - 1);

            var data = new T[length];
            for (int i = 0; i < data.Length; i++)
            {
                data[i] = map(Math.Pow(10, startExponent + i*step));
            }
            data[data.Length - 1] = map(Math.Pow(10, stopExponent));
            return data;
        }

        /// <summary>
        /// Generate a linearly spaced sample vector within the inclusive interval (start, stop) and step 1.
        /// Equivalent to MATLAB colon operator (:).
        /// </summary>
        public static double[] LinearRange(int start, int stop)
        {
            if (start == stop) return new double[] { start };
            if (start < stop)
            {
                var data = new double[stop - start + 1];
                for (int i = 0; i < data.Length; i++)
                {
                    data[i] = start + i;
                }
                return data;
            }
            else
            {
                var data = new double[start - stop + 1];
                for (int i = 0; i < data.Length; i++)
                {
                    data[i] = start - i;
                }
                return data;
            }
        }

        /// <summary>
        /// Generate a linearly spaced sample vector within the inclusive interval (start, stop) and step 1.
        /// Equivalent to MATLAB colon operator (:).
        /// </summary>
        public static int[] LinearRangeInt32(int start, int stop)
        {
            if (start == stop) return new int[] { start };
            if (start < stop)
            {
                var data = new int[stop - start + 1];
                for (int i = 0; i < data.Length; i++)
                {
                    data[i] = start + i;
                }
                return data;
            }
            else
            {
                var data = new int[start - stop + 1];
                for (int i = 0; i < data.Length; i++)
                {
                    data[i] = start - i;
                }
                return data;
            }
        }

        /// <summary>
        /// Generate a linearly spaced sample vector within the inclusive interval (start, stop) and the provided step.
        /// The start value is aways included as first value, but stop is only included if it stop-start is a multiple of step.
        /// Equivalent to MATLAB double colon operator (::).
        /// </summary>
        public static double[] LinearRange(int start, int step, int stop)
        {
            if (start == stop) return new double[] { start };
            if (start < stop && step < 0 || start > stop && step > 0 || step == 0d)
            {
                return new double[0];
            }

            var data = new double[(stop - start)/step + 1];
            for (int i = 0; i < data.Length; i++)
            {
                data[i] = start + i*step;
            }
            return data;
        }

        /// <summary>
        /// Generate a linearly spaced sample vector within the inclusive interval (start, stop) and the provided step.
        /// The start value is aways included as first value, but stop is only included if it stop-start is a multiple of step.
        /// Equivalent to MATLAB double colon operator (::).
        /// </summary>
        public static int[] LinearRangeInt32(int start, int step, int stop)
        {
            if (start == stop) return new int[] { start };
            if (start < stop && step < 0 || start > stop && step > 0 || step == 0d)
            {
                return new int[0];
            }

            var data = new int[(stop - start) / step + 1];
            for (int i = 0; i < data.Length; i++)
            {
                data[i] = start + i * step;
            }
            return data;
        }

        /// <summary>
        /// Generate a linearly spaced sample vector within the inclusive interval (start, stop) and the provide step.
        /// The start value is aways included as first value, but stop is only included if it stop-start is a multiple of step.
        /// Equivalent to MATLAB double colon operator (::).
        /// </summary>
        public static double[] LinearRange(double start, double step, double stop)
        {
            if (start == stop) return new[] { start };
            if (start < stop && step < 0 || start > stop && step > 0 || step == 0d)
            {
                return new double[0];
            }

            var data = new double[(int)Math.Floor((stop - start)/step + 1d)];
            for (int i = 0; i < data.Length; i++)
            {
                data[i] = start + i*step;
            }
            return data;
        }

        /// <summary>
        /// Generate samples by sampling a function at linearly spaced points within the inclusive interval (start, stop) and the provide step.
        /// The start value is aways included as first value, but stop is only included if it stop-start is a multiple of step.
        /// </summary>
        public static T[] LinearRangeMap<T>(double start, double step, double stop, Func<double, T> map)
        {
            if (start == stop) return new[] { map(start) };
            if (start < stop && step < 0 || start > stop && step > 0 || step == 0d)
            {
                return new T[0];
            }

            var data = new T[(int)Math.Floor((stop - start)/step + 1d)];
            for (int i = 0; i < data.Length; i++)
            {
                data[i] = map(start + i*step);
            }
            return data;
        }

        /// <summary>
        /// Create a periodic wave.
        /// </summary>
        /// <param name="length">The number of samples to generate.</param>
        /// <param name="samplingRate">Samples per time unit (Hz). Must be larger than twice the frequency to satisfy the Nyquist criterion.</param>
        /// <param name="frequency">Frequency in periods per time unit (Hz).</param>
        /// <param name="amplitude">The length of the period when sampled at one sample per time unit. This is the interval of the periodic domain, a typical value is 1.0, or 2*Pi for angular functions.</param>
        /// <param name="phase">Optional phase offset.</param>
        /// <param name="delay">Optional delay, relative to the phase.</param>
        public static double[] Periodic(int length, double samplingRate, double frequency, double amplitude = 1.0, double phase = 0.0, int delay = 0)
        {
            if (length < 0)
            {
                throw new ArgumentOutOfRangeException(nameof(length));
            }

            double step = frequency/samplingRate*amplitude;
            phase = Euclid.Modulus(phase - delay*step, amplitude);

            var data = new double[length];
            for (int i = 0, k = 0; i < data.Length; i++, k++)
            {
                var x = phase + k*step;
                if (x >= amplitude)
                {
                    x %= amplitude;
                    phase = x;
                    k = 0;
                }

                data[i] = x;
            }
            return data;
        }

        /// <summary>
        /// Create a periodic wave.
        /// </summary>
        /// <param name="length">The number of samples to generate.</param>
        /// <param name="map">The function to apply to each of the values and evaluate the resulting sample.</param>
        /// <param name="samplingRate">Samples per time unit (Hz). Must be larger than twice the frequency to satisfy the Nyquist criterion.</param>
        /// <param name="frequency">Frequency in periods per time unit (Hz).</param>
        /// <param name="amplitude">The length of the period when sampled at one sample per time unit. This is the interval of the periodic domain, a typical value is 1.0, or 2*Pi for angular functions.</param>
        /// <param name="phase">Optional phase offset.</param>
        /// <param name="delay">Optional delay, relative to the phase.</param>
        public static T[] PeriodicMap<T>(int length, Func<double, T> map, double samplingRate, double frequency, double amplitude = 1.0, double phase = 0.0, int delay = 0)
        {
            if (length < 0)
            {
                throw new ArgumentOutOfRangeException(nameof(length));
            }

            double step = frequency/samplingRate*amplitude;
            phase = Euclid.Modulus(phase - delay*step, amplitude);

            var data = new T[length];
            for (int i = 0, k = 0; i < data.Length; i++, k++)
            {
                var x = phase + k*step;
                if (x >= amplitude)
                {
                    x %= amplitude;
                    phase = x;
                    k = 0;
                }

                data[i] = map(x);
            }
            return data;
        }

        /// <summary>
        /// Create an infinite periodic wave sequence.
        /// </summary>
        /// <param name="samplingRate">Samples per time unit (Hz). Must be larger than twice the frequency to satisfy the Nyquist criterion.</param>
        /// <param name="frequency">Frequency in periods per time unit (Hz).</param>
        /// <param name="amplitude">The length of the period when sampled at one sample per time unit. This is the interval of the periodic domain, a typical value is 1.0, or 2*Pi for angular functions.</param>
        /// <param name="phase">Optional phase offset.</param>
        /// <param name="delay">Optional delay, relative to the phase.</param>
        public static IEnumerable<double> PeriodicSequence(double samplingRate, double frequency, double amplitude = 1.0, double phase = 0.0, int delay = 0)
        {
            double step = frequency/samplingRate*amplitude;
            phase = Euclid.Modulus(phase - delay*step, amplitude);

            int k = 0;
            while (true)
            {
                var x = phase + (k++)*step;
                if (x >= amplitude)
                {
                    x %= amplitude;
                    phase = x;
                    k = 1;
                }

                yield return x;
            }
        }

        /// <summary>
        /// Create an infinite periodic wave sequence.
        /// </summary>
        /// <param name="map">The function to apply to each of the values and evaluate the resulting sample.</param>
        /// <param name="samplingRate">Samples per time unit (Hz). Must be larger than twice the frequency to satisfy the Nyquist criterion.</param>
        /// <param name="frequency">Frequency in periods per time unit (Hz).</param>
        /// <param name="amplitude">The length of the period when sampled at one sample per time unit. This is the interval of the periodic domain, a typical value is 1.0, or 2*Pi for angular functions.</param>
        /// <param name="phase">Optional phase offset.</param>
        /// <param name="delay">Optional delay, relative to the phase.</param>
        public static IEnumerable<T> PeriodicMapSequence<T>(Func<double, T> map, double samplingRate, double frequency, double amplitude = 1.0, double phase = 0.0, int delay = 0)
        {
            double step = frequency/samplingRate*amplitude;
            phase = Euclid.Modulus(phase - delay*step, amplitude);

            int k = 0;
            while (true)
            {
                var x = phase + (k++)*step;
                if (x >= amplitude)
                {
                    x %= amplitude;
                    phase = x;
                    k = 1;
                }

                yield return map(x);
            }
        }

        /// <summary>
        /// Create a Sine wave.
        /// </summary>
        /// <param name="length">The number of samples to generate.</param>
        /// <param name="samplingRate">Samples per time unit (Hz). Must be larger than twice the frequency to satisfy the Nyquist criterion.</param>
        /// <param name="frequency">Frequency in periods per time unit (Hz).</param>
        /// <param name="amplitude">The maximal reached peak.</param>
        /// <param name="mean">The mean, or DC part, of the signal.</param>
        /// <param name="phase">Optional phase offset.</param>
        /// <param name="delay">Optional delay, relative to the phase.</param>
        public static double[] Sinusoidal(int length, double samplingRate, double frequency, double amplitude, double mean = 0.0, double phase = 0.0, int delay = 0)
        {
            if (length < 0)
            {
                throw new ArgumentOutOfRangeException(nameof(length));
            }

            double step = frequency/samplingRate*Constants.Pi2;
            phase = (phase - delay*step)%Constants.Pi2;

            var data = new double[length];
            for (int i = 0; i < data.Length; i++)
            {
                data[i] = mean + amplitude*Math.Sin(phase + i*step);
            }
            return data;
        }

        /// <summary>
        /// Create an infinite Sine wave sequence.
        /// </summary>
        /// <param name="samplingRate">Samples per unit.</param>
        /// <param name="frequency">Frequency in samples per unit.</param>
        /// <param name="amplitude">The maximal reached peak.</param>
        /// <param name="mean">The mean, or DC part, of the signal.</param>
        /// <param name="phase">Optional phase offset.</param>
        /// <param name="delay">Optional delay, relative to the phase.</param>
        public static IEnumerable<double> SinusoidalSequence(double samplingRate, double frequency, double amplitude, double mean = 0.0, double phase = 0.0, int delay = 0)
        {
            double step = frequency/samplingRate*Constants.Pi2;
            phase = (phase - delay*step)%Constants.Pi2;

            while (true)
            {
                for (int i = 0; i < 1000; i++)
                {
                    yield return mean + amplitude*Math.Sin(phase + i*step);
                }
                phase = (phase + 1000*step)%Constants.Pi2;
            }
        }

        /// <summary>
        /// Create a periodic square wave, starting with the high phase.
        /// </summary>
        /// <param name="length">The number of samples to generate.</param>
        /// <param name="highDuration">Number of samples of the high phase.</param>
        /// <param name="lowDuration">Number of samples of the low phase.</param>
        /// <param name="lowValue">Sample value to be emitted during the low phase.</param>
        /// <param name="highValue">Sample value to be emitted during the high phase.</param>
        /// <param name="delay">Optional delay.</param>
        public static double[] Square(int length, int highDuration, int lowDuration, double lowValue, double highValue, int delay = 0)
        {
            var duration = highDuration + lowDuration;
            return PeriodicMap(length, x => x < highDuration ? highValue : lowValue, 1.0, 1.0/duration, duration, 0.0, delay);
        }

        /// <summary>
        /// Create an infinite periodic square wave sequence, starting with the high phase.
        /// </summary>
        /// <param name="highDuration">Number of samples of the high phase.</param>
        /// <param name="lowDuration">Number of samples of the low phase.</param>
        /// <param name="lowValue">Sample value to be emitted during the low phase.</param>
        /// <param name="highValue">Sample value to be emitted during the high phase.</param>
        /// <param name="delay">Optional delay.</param>
        public static IEnumerable<double> SquareSequence(int highDuration, int lowDuration, double lowValue, double highValue, int delay = 0)
        {
            var duration = highDuration + lowDuration;
            return PeriodicMapSequence(x => x < highDuration ? highValue : lowValue, 1.0, 1.0/duration, duration, 0.0, delay);
        }

        /// <summary>
        /// Create a periodic triangle wave, starting with the raise phase from the lowest sample.
        /// </summary>
        /// <param name="length">The number of samples to generate.</param>
        /// <param name="raiseDuration">Number of samples of the raise phase.</param>
        /// <param name="fallDuration">Number of samples of the fall phase.</param>
        /// <param name="lowValue">Lowest sample value.</param>
        /// <param name="highValue">Highest sample value.</param>
        /// <param name="delay">Optional delay.</param>
        public static double[] Triangle(int length, int raiseDuration, int fallDuration, double lowValue, double highValue, int delay = 0)
        {
            var duration = raiseDuration + fallDuration;
            var height = highValue - lowValue;
            var raise = height / raiseDuration;
            var fall = height / fallDuration;
            return PeriodicMap(length, x => x < raiseDuration ? lowValue + x*raise : highValue - (x-raiseDuration)*fall, 1.0, 1.0/duration, duration, 0.0, delay);
        }

        /// <summary>
        /// Create an infinite periodic triangle wave sequence, starting with the raise phase from the lowest sample.
        /// </summary>
        /// <param name="raiseDuration">Number of samples of the raise phase.</param>
        /// <param name="fallDuration">Number of samples of the fall phase.</param>
        /// <param name="lowValue">Lowest sample value.</param>
        /// <param name="highValue">Highest sample value.</param>
        /// <param name="delay">Optional delay.</param>
        public static IEnumerable<double> TriangleSequence(int raiseDuration, int fallDuration, double lowValue, double highValue, int delay = 0)
        {
            var duration = raiseDuration + fallDuration;
            var height = highValue - lowValue;
            var raise = height / raiseDuration;
            var fall = height / fallDuration;
            return PeriodicMapSequence(x => x < raiseDuration ? lowValue + x*raise : highValue - (x-raiseDuration)*fall, 1.0, 1.0/duration, duration, 0.0, delay);
        }

        /// <summary>
        /// Create a periodic sawtooth wave, starting with the lowest sample.
        /// </summary>
        /// <param name="length">The number of samples to generate.</param>
        /// <param name="period">Number of samples a full sawtooth period.</param>
        /// <param name="lowValue">Lowest sample value.</param>
        /// <param name="highValue">Highest sample value.</param>
        /// <param name="delay">Optional delay.</param>
        public static double[] Sawtooth(int length, int period, double lowValue, double highValue, int delay = 0)
        {
            var height = highValue - lowValue;
            return PeriodicMap(length, x => x + lowValue, 1.0, 1.0/period, height*period/(period-1), 0.0, delay);
        }

        /// <summary>
        /// Create an infinite periodic sawtooth wave sequence, starting with the lowest sample.
        /// </summary>
        /// <param name="period">Number of samples a full sawtooth period.</param>
        /// <param name="lowValue">Lowest sample value.</param>
        /// <param name="highValue">Highest sample value.</param>
        /// <param name="delay">Optional delay.</param>
        public static IEnumerable<double> SawtoothSequence(int period, double lowValue, double highValue, int delay = 0)
        {
            var height = highValue - lowValue;
            return PeriodicMapSequence(x => x + lowValue, 1.0, 1.0/period, height*period/(period-1), 0.0, delay);
        }

        /// <summary>
        /// Create an array with each field set to the same value.
        /// </summary>
        /// <param name="length">The number of samples to generate.</param>
        /// <param name="value">The value that each field should be set to.</param>
        public static T[] Repeat<T>(int length, T value)
        {
            if (length < 0)
            {
                throw new ArgumentOutOfRangeException(nameof(length));
            }

            var data = new T[length];
            CommonParallel.For(0, data.Length, 4096, (a, b) =>
            {
                for (int i = a; i < b; i++)
                {
                    data[i] = value;
                }
            });
            return data;
        }

        /// <summary>
        /// Create an infinite sequence where each element has the same value.
        /// </summary>
        /// <param name="value">The value that each element should be set to.</param>
        public static IEnumerable<T> RepeatSequence<T>(T value)
        {
            while (true)
            {
                yield return value;
            }
        }

        /// <summary>
        /// Create a Heaviside Step sample vector.
        /// </summary>
        /// <param name="length">The number of samples to generate.</param>
        /// <param name="amplitude">The maximal reached peak.</param>
        /// <param name="delay">Offset to the time axis.</param>
        public static double[] Step(int length, double amplitude, int delay)
        {
            if (length < 0)
            {
                throw new ArgumentOutOfRangeException(nameof(length));
            }

            var data = new double[length];
            for (int i = Math.Max(0, delay); i < data.Length; i++)
            {
                data[i] = amplitude;
            }
            return data;
        }

        /// <summary>
        /// Create an infinite Heaviside Step sample sequence.
        /// </summary>
        /// <param name="amplitude">The maximal reached peak.</param>
        /// <param name="delay">Offset to the time axis.</param>
        public static IEnumerable<double> StepSequence(double amplitude, int delay)
        {
            for (int i = 0; i < delay; i++)
            {
                yield return 0d;
            }

            while (true)
            {
                yield return amplitude;
            }
        }

        /// <summary>
        /// Create a Kronecker Delta impulse sample vector.
        /// </summary>
        /// <param name="length">The number of samples to generate.</param>
        /// <param name="amplitude">The maximal reached peak.</param>
        /// <param name="delay">Offset to the time axis. Zero or positive.</param>
        public static double[] Impulse(int length, double amplitude, int delay)
        {
            if (length < 0)
            {
                throw new ArgumentOutOfRangeException(nameof(length));
            }

            var data = new double[length];
            if (delay >= 0 && delay < length)
            {
                data[delay] = amplitude;
            }
            return data;
        }

        /// <summary>
        /// Create a Kronecker Delta impulse sample vector.
        /// </summary>
        /// <param name="amplitude">The maximal reached peak.</param>
        /// <param name="delay">Offset to the time axis, hence the sample index of the impulse.</param>
        public static IEnumerable<double> ImpulseSequence(double amplitude, int delay)
        {
            if (delay >= 0)
            {
                for (int i = 0; i < delay; i++)
                {
                    yield return 0d;
                }

                yield return amplitude;
            }

            while (true)
            {
                yield return 0d;
            }
        }

        /// <summary>
        /// Create a periodic Kronecker Delta impulse sample vector.
        /// </summary>
        /// <param name="length">The number of samples to generate.</param>
        /// <param name="period">impulse sequence period.</param>
        /// <param name="amplitude">The maximal reached peak.</param>
        /// <param name="delay">Offset to the time axis. Zero or positive.</param>
        public static double[] PeriodicImpulse(int length, int period, double amplitude, int delay)
        {
            if (length < 0)
            {
                throw new ArgumentOutOfRangeException(nameof(length));
            }

            var data = new double[length];
            delay = Euclid.Modulus(delay, period);
            while (delay < length)
            {
                data[delay] = amplitude;
                delay += period;
            }
            return data;
        }

        /// <summary>
        /// Create a Kronecker Delta impulse sample vector.
        /// </summary>
        /// <param name="period">impulse sequence period.</param>
        /// <param name="amplitude">The maximal reached peak.</param>
        /// <param name="delay">Offset to the time axis. Zero or positive.</param>
        public static IEnumerable<double> PeriodicImpulseSequence(int period, double amplitude, int delay)
        {
            delay = Euclid.Modulus(delay, period);

            for (int i = 0; i < delay; i++)
            {
                yield return 0d;
            }

            while (true)
            {
                yield return amplitude;

                for (int i = 1; i < period; i++)
                {
                    yield return 0d;
                }
            }
        }

        /// <summary>
        /// Generate samples generated by the given computation.
        /// </summary>
        public static T[] Unfold<T, TState>(int length, Func<TState, Tuple<T, TState>> f, TState state)
        {
            if (length < 0)
            {
                throw new ArgumentOutOfRangeException(nameof(length));
            }

            var data = new T[length];
            for (int i = 0; i < data.Length; i++)
            {
                Tuple<T, TState> next = f(state);
                data[i] = next.Item1;
                state = next.Item2;
            }
            return data;
        }

        /// <summary>
        /// Generate an infinite sequence generated by the given computation.
        /// </summary>
        public static IEnumerable<T> UnfoldSequence<T, TState>(Func<TState, Tuple<T, TState>> f, TState state)
        {
            while (true)
            {
                Tuple<T, TState> next = f(state);
                state = next.Item2;
                yield return next.Item1;
            }
        }

        /// <summary>
        /// Generate a Fibonacci sequence, including zero as first value.
        /// </summary>
        public static BigInteger[] Fibonacci(int length)
        {
            if (length < 0)
            {
                throw new ArgumentOutOfRangeException(nameof(length));
            }

            var data = new BigInteger[length];
            if (data.Length > 0)
            {
                data[0] = BigInteger.Zero;
            }
            if (data.Length > 1)
            {
                data[1] = BigInteger.One;
            }
            for (int i = 2; i < data.Length; i++)
            {
                data[i] = data[i - 1] + data[i - 2];
            }
            return data;
        }

        /// <summary>
        /// Generate an infinite Fibonacci sequence, including zero as first value.
        /// </summary>
        public static IEnumerable<BigInteger> FibonacciSequence()
        {
            BigInteger a = BigInteger.Zero;
            yield return a;

            BigInteger b = BigInteger.One;
            yield return b;

            while (true)
            {
                a = a + b;
                yield return a;

                b = a + b;
                yield return b;
            }
        }

        /// <summary>
        /// Create random samples, uniform between 0 and 1.
        /// Faster than other methods but with reduced guarantees on randomness.
        /// </summary>
        public static double[] Uniform(int length)
        {
            if (length < 0)
            {
                throw new ArgumentOutOfRangeException(nameof(length));
            }

            return SystemRandomSource.FastDoubles(length);
        }

        /// <summary>
        /// Create an infinite random sample sequence, uniform between 0 and 1.
        /// Faster than other methods but with reduced guarantees on randomness.
        /// </summary>
        public static IEnumerable<double> UniformSequence()
        {
            return SystemRandomSource.DoubleSequence();
        }


        /// <summary>
        /// Generate samples by sampling a function at samples from a probability distribution, uniform between 0 and 1.
        /// Faster than other methods but with reduced guarantees on randomness.
        /// </summary>
        public static T[] UniformMap<T>(int length, Func<double, T> map)
        {
            if (length < 0)
            {
                throw new ArgumentOutOfRangeException(nameof(length));
            }

            var samples = SystemRandomSource.FastDoubles(length);
            return Map(samples, map);
        }

        /// <summary>
        /// Generate a sample sequence by sampling a function at samples from a probability distribution, uniform between 0 and 1.
        /// Faster than other methods but with reduced guarantees on randomness.
        /// </summary>
        public static IEnumerable<T> UniformMapSequence<T>(Func<double, T> map)
        {
            return SystemRandomSource.DoubleSequence().Select(map);
        }

        /// <summary>
        /// Generate samples by sampling a function at sample pairs from a probability distribution, uniform between 0 and 1.
        /// Faster than other methods but with reduced guarantees on randomness.
        /// </summary>
        public static T[] UniformMap2<T>(int length, Func<double, double, T> map)
        {
            if (length < 0)
            {
                throw new ArgumentOutOfRangeException(nameof(length));
            }

            var samples1 = SystemRandomSource.FastDoubles(length);
            var samples2 = SystemRandomSource.FastDoubles(length);
            return Map2(samples1, samples2, map);
        }

        /// <summary>
        /// Generate a sample sequence by sampling a function at sample pairs from a probability distribution, uniform between 0 and 1.
        /// Faster than other methods but with reduced guarantees on randomness.
        /// </summary>
        public static IEnumerable<T> UniformMap2Sequence<T>(Func<double, double, T> map)
        {
            var rnd1 = SystemRandomSource.Default;
            for (int i = 0; i < 128; i++)
            {
                yield return map(rnd1.NextDouble(), rnd1.NextDouble());
            }

            var rnd2 = new System.Random(RandomSeed.Robust());
            while (true)
            {
                yield return map(rnd2.NextDouble(), rnd2.NextDouble());
            }
        }

        /// <summary>
        /// Create samples with independent amplitudes of standard distribution.
        /// </summary>
        public static double[] Standard(int length)
        {
            return Normal(length, 0.0, 1.0);
        }

        /// <summary>
        /// Create an infinite sample sequence with independent amplitudes of standard distribution.
        /// </summary>
        public static IEnumerable<double> StandardSequence()
        {
            return NormalSequence(0.0, 1.0);
        }

        /// <summary>
        /// Create samples with independent amplitudes of normal distribution and a flat spectral density.
        /// </summary>
        public static double[] Normal(int length, double mean, double standardDeviation)
        {
            if (length < 0)
            {
                throw new ArgumentOutOfRangeException(nameof(length));
            }

            var samples = new double[length];
            Distributions.Normal.Samples(SystemRandomSource.Default, samples, mean, standardDeviation);
            return samples;
        }

        /// <summary>
        /// Create an infinite sample sequence with independent amplitudes of normal distribution and a flat spectral density.
        /// </summary>
        public static IEnumerable<double> NormalSequence(double mean, double standardDeviation)
        {
            return Distributions.Normal.Samples(SystemRandomSource.Default, mean, standardDeviation);
        }

        /// <summary>
        /// Create random samples.
        /// </summary>
        public static double[] Random(int length, IContinuousDistribution distribution)
        {
            if (length < 0)
            {
                throw new ArgumentOutOfRangeException(nameof(length));
            }

            var samples = new double[length];
            distribution.Samples(samples);
            return samples;
        }

        /// <summary>
        /// Create an infinite random sample sequence.
        /// </summary>
        public static IEnumerable<double> Random(IContinuousDistribution distribution)
        {
            return distribution.Samples();
        }

        /// <summary>
        /// Create random samples.
        /// </summary>
        public static float[] RandomSingle(int length, IContinuousDistribution distribution)
        {
            if (length < 0)
            {
                throw new ArgumentOutOfRangeException(nameof(length));
            }

            var samples = new double[length];
            distribution.Samples(samples);
            return Map(samples, v => (float)v);
        }

        /// <summary>
        /// Create an infinite random sample sequence.
        /// </summary>
        public static IEnumerable<float> RandomSingle(IContinuousDistribution distribution)
        {
            return distribution.Samples().Select(v => (float)v);
        }

        /// <summary>
        /// Create random samples.
        /// </summary>
        public static Complex[] RandomComplex(int length, IContinuousDistribution distribution)
        {
            return RandomMap2(length, distribution, (r, i) => new Complex(r, i));
        }

        /// <summary>
        /// Create an infinite random sample sequence.
        /// </summary>
        public static IEnumerable<Complex> RandomComplex(IContinuousDistribution distribution)
        {
            return RandomMap2Sequence(distribution, (r, i) => new Complex(r, i));
        }

        /// <summary>
        /// Create random samples.
        /// </summary>
        public static Complex32[] RandomComplex32(int length, IContinuousDistribution distribution)
        {
            return RandomMap2(length, distribution, (r, i) => new Complex32((float)r, (float)i));
        }

        /// <summary>
        /// Create an infinite random sample sequence.
        /// </summary>
        public static IEnumerable<Complex32> RandomComplex32(IContinuousDistribution distribution)
        {
            return RandomMap2Sequence(distribution, (r, i) => new Complex32((float)r, (float)i));
        }

        /// <summary>
        /// Generate samples by sampling a function at samples from a probability distribution.
        /// </summary>
        public static T[] RandomMap<T>(int length, IContinuousDistribution distribution, Func<double, T> map)
        {
            if (length < 0)
            {
                throw new ArgumentOutOfRangeException(nameof(length));
            }

            var samples = new double[length];
            distribution.Samples(samples);
            return Map(samples, map);
        }

        /// <summary>
        /// Generate a sample sequence by sampling a function at samples from a probability distribution.
        /// </summary>
        public static IEnumerable<T> RandomMapSequence<T>(IContinuousDistribution distribution, Func<double, T> map)
        {
            return distribution.Samples().Select(map);
        }

        /// <summary>
        /// Generate samples by sampling a function at sample pairs from a probability distribution.
        /// </summary>
        public static T[] RandomMap2<T>(int length, IContinuousDistribution distribution, Func<double, double, T> map)
        {
            if (length < 0)
            {
                throw new ArgumentOutOfRangeException(nameof(length));
            }

            var samples1 = new double[length];
            var samples2 = new double[length];
            distribution.Samples(samples1);
            distribution.Samples(samples2);
            return Map2(samples1, samples2, map);
        }

        /// <summary>
        /// Generate a sample sequence by sampling a function at sample pairs from a probability distribution.
        /// </summary>
        public static IEnumerable<T> RandomMap2Sequence<T>(IContinuousDistribution distribution, Func<double, double, T> map)
        {
            return distribution.Samples().Zip(distribution.Samples(), map);
        }
    }
}