using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.IO;


namespace FractalDimension
{
    public class Point
    {
        public double x;
        public double y;
        public Point(double _x, double _y)
        {
            x = _x;
            y = _y;
        }
    }
    public class FractalDimension
    {
        // Вычисление среднего
        private static double calcMean(double[] v, int beg, int N)
        {
            double sum = 0.0;
            for (int i = 0; i < N; i++)
            {
                sum = sum + v[i + beg];
            }
            double mean = sum / N;
            return mean;
        }

        /*
            Вычисление R/S для блока данных
            \arg v массив с блоком данных;
            \arg boxSize размер блока;
            \arg beg начало блока в массиве.
        */
        private static double calcRS(double[] v, int beg, int boxSize)
        {
            double RS = 0.0;
            if (v.Length != 0 && boxSize > 0)
            {
                double min;
                double max;
                double runningSum;
                double runningSumSqr;

                double mean = calcMean(v, beg, boxSize);

                min = 0.0;
                max = 0.0;
                runningSum = 0.0;
                runningSumSqr = 0.0;
                for (int i = 0; i < boxSize; i++)
                {
                    double devFromMean = v[i + beg] - mean;
                    runningSum = runningSum + devFromMean;
                    runningSumSqr = runningSumSqr + (devFromMean * devFromMean);
                    if (runningSum < min)
                        min = runningSum;
                    if (runningSum > max)
                        max = runningSum;
                }
                double variance = runningSumSqr / Convert.ToDouble(boxSize);
                double stdDev = Math.Sqrt(variance);

                double range = max - min;
                RS = range / stdDev;
            }
            return RS;
        }
       
        /*
            Вычисление среднего R/S по всем блокам размера boxSize.
            В этой реализации применяется допущение, что N делится нацело на
            boxSize.
            \arg v массив данных для которого надо вычислить R/S;
            \arg N размер массива;
            \arg boxSize размер боксов, по которым будет считаться R/S.
        */
        private static double calcRSAve(double[] v, int N, int boxSize)
        {
            int i;
            double RS;
            double RSSum;
            double RSAve =0.0;

            int numBoxes = N/boxSize;
            if (numBoxes > 0)
            {
                RSSum = 0;
                for (i = 0; i+boxSize <= N; i = i + boxSize)
                {
                    RS = calcRS(v, i, boxSize);
                    RSSum = RSSum + RS;
                }
                RSAve = RSSum / Convert.ToDouble(numBoxes);
            }
            return RSAve;
        }

        /*
           Вычисление показателя Херста. Ряд делится на блоки размера равного степени
           2-ки (напр., 8, 16, 32, 64, ..). Для каждого блока вычисляется R/S и для
           всей совокупности блоков вычисляется среднее R/S. Изначально n = N, на каждой
           итерации алгоритма n уменьшается в 2 раза, пока не станет равно 8. 
           Набор полученных точек наносим на логарифмические оси и проводим 
           линейную регрессию по методу наименьших квадратов.
           \arg v массив данных для которого надо вычислить R/S;
           \arg N размер массива.
        */
        public static double calcHurstEst(double[] v, int N)
        {
            if (v.Length != 0 && N > 0)
            {
                List<Point> list = new List<Point>();
                int minBox = 8;
                int boxSize;
                for (boxSize = N; boxSize >= minBox; boxSize = (boxSize / 2))
                {
                    double RSAve = calcRSAve(v, N, boxSize);
                    double logRSAve = Math.Log(RSAve, 2);
                    double logBoxSize = Math.Log(boxSize, 2);
                    Point p = new Point(logBoxSize, logRSAve);
                    list.Add(p);
                }
                return lregress(list);
            }
            return -1;
        }

        /*
           Вычисление среднего по компоненте x для списка точек
        */
        private static double meanX(List<Point> list)
        {
            double mean = 0.0;
            double sum = 0.0;

            foreach (Point p in list)
            {
                sum = sum + p.x;
            }
            mean = sum / list.Count;
            return mean;
        }
        /*
           Вычисление среднего по компоненте y для списка точек
        */
        private static double meanY(List<Point> list)
        {
            double mean = 0.0;
            double sum = 0.0;

            foreach (Point p in list)
            {
                sum = sum + p.y;
            }
            mean = sum / list.Count;
            return mean;
        }
        /*
            Линейная регрессия
        */
        private static double lregress(List<Point> list)
        {
            int N = list.Count;
            if (N > 0)
            {
                double muX = meanX(list);
                double muY = meanY(list);

                //     N-1
                //     ---
                //     \   (Xi - meanX)(Yi - meanY)
                //     /__ 
                //     i=0
                // b = -----------------------------
                //     N-1
                //     ---             2
                //     \   (Xi - meanX)
                //     /__ 
                //     i=0
                //

                double SSxy = 0;
                double SSxx = 0;
                for (int i = 0; i < N; i++)
                {
                    double subX = (list[i].x - muX);
                    double subY = (list[i].y - muY);
                    SSxy = SSxy + subX * subY;
                    SSxx = SSxx + subX * subX;
                }

                double b = SSxy / SSxx;

                double a = muY - b * muX;
                return b;
            }
            return -1;
        }
        /*
            Рандомизированная линейная регрессия
        */
        private static double lregressRandomized(List<Point> list, int randMax)
        {
            int N = list.Count;
            if (N > 0)
            {
                double muX = meanX(list);
                double muY = meanY(list);


                //     N-1
                //     ---
                //     \   n(Xi - meanX)(Yi - meanY)
                //     /__ 
                //     i=0
                // b = -----------------------------
                //     N-1
                //     ---             2
                //     \   n(Xi - meanX)
                //     /__ 
                //     i=0
                // n - для каждого i генерируется случайно из (0,5)

                double SSxy = 0;
                double SSxx = 0;
                Random rand = new Random();
                for (int i = 0; i < N; i++)
                {
                    int n = rand.Next(0, randMax);
                    double subX = (list[i].x - muX);
                    double subY = (list[i].y - muY);
                    SSxy = SSxy + n * subX * subY;
                    SSxx = SSxx + n * subX * subX;
                }

                double b = SSxy / SSxx;

                double a = muY - b * muX;
                return b;
            }
            return -1;
        }
        /*
            Вычисление показателя Херста по рандомизированному алгоритму R/S-анализа.
            Используется рандомизированная линейная регрессия.
        */
        public static double calcHurstEstRandomized(double[] v, int N, int randMax)
        {
            if (v.Length != 0 && N > 0)
            {
                List<Point> list = new List<Point>();
                int minBox = 2;
                int boxSize;
                for (boxSize = N; boxSize >= minBox; boxSize = (boxSize / 2))
                {
                    double RSAve = calcRSAve(v, N, boxSize);
                    double logRSAve = Math.Log(RSAve, 2);
                    double logBoxSize = Math.Log(boxSize, 2);
                    Point p = new Point(logBoxSize, logRSAve);
                    list.Add(p);
                }
                return lregressRandomized(list, randMax);
            }
            return -1;
        }
    }
    public class Generator
    {
        //помеха в виде ступеньки
        public static double[] addNoise1(double[] a, double height)
        {
            for (int i = 300; i < 800; i++)
            {
                a[i] = a[i] + height;
            }
            return a;
        }
        //помеха в виде ряда ступенек
        public static double[] addNoise2(double[] a, double height)
        {
            for (int j = 1; j < 9; j++)
            {
                for (int i = 0; i < 50; i++)
                {
                    a[i] = a[100 * j + i] + height;
                }
            }
            return a;
        }
        //помеха в виде ступеньки с фрактальной структурой
        public static double[] addNoise3(double[] a, double[] b)
        {

            for (int i = 700; i < 800; i++)
            {
                a[i] = a[i] + b[i];
            }
            return a;
        }

        /*
            генерация фрактального ряда с H=0,5. На нем R/S должен выдавать около 5,5.
            \arg rand1 - равномерно распределенная на [0..1] последовательность чисел.
            \arg noise - вид помехи
            \arg height - высота ступени в случае ступенчатых помех
            \arg frIns - фрактальная вставка в случае фрактальной помехи
        */
        public static double[] gen(Random rand1, int noise, double height, double[] frIns)
        {
            double[] a = new double[1025];
            a[0] = 50.0;
            double t = 1.0 / 1024.0;
            for (int i = 1; i <= 1024; i++)
            {
                double r = rand1.NextDouble();
                double f = rand1.NextDouble();
                a[i] = a[i - 1] + Math.Sqrt(t) * (Math.Cos(2 * Math.PI * f) * Math.Sqrt((-1) * 2 * Math.Log(r, Math.E)));
            }
            switch (noise)
            {
                // нет помехи, ничего не делаем
                case 0:
                    break;
                // ступенчатая помеха высоты height
                case 1:
                    a = addNoise1(a, height);
                    break;
                // последовательность ступенек высоты height
                case 2:
                    a = addNoise2(a, height);
                    break;
                // фрактальная вставка
                case 3:
                    a = addNoise3(a, frIns);
                    break;
                default:
                    break;
            }
            double[] b = new double[1024];
            for (int i = 0; i < 1023; i++)
            {
                b[i] = Math.Log(a[i + 1], 2) - Math.Log(a[i], 2);
            }
            return b;
        }

        /*
            генерация фрактального ряда с произвольным H. Желательно брать больше 0,4, иначе метод
            дает плохой результат.
            \arg H - фрактальная размерность ряда. 
        */
        public static double[] genFractDim(double H)
        {
            Random rand1 = new Random();
            double[] e = new double[80000];
            for (int j = 0; j < 80000; j++)
            {
                double r = rand1.NextDouble();
                double f = rand1.NextDouble();
                e[j] = Math.Cos(2 * Math.PI * f) * Math.Sqrt((-1) * 2 * Math.Log(r, Math.E));
            }
            double[] a = new double[1024];
            for (int j = 0; j < 1024; j++)
            {
                a[j] = 0;
                for (int n = 1; n < 25001; n++)
                {
                    a[j] = a[j] + KHx(H, n / 100.0) * e[j * 50 - n + 25000];
                }
            }
            return a;
        }
        //вспомогательная функция.
        private static double KHx(double H, double x)
        {
            if (x > 1) return (Math.Exp((H - 0.5) * Math.Log(x, Math.E)) - Math.Exp((H - 0.5) * Math.Log(x - 1, Math.E)));
            else return (Math.Exp((H - 0.5) * Math.Log(x, Math.E)));
        }
    }

    class Program
    {
        //пример применения (для чистого ряда с фрактальной размерностью 0,5)
        // Если надо брать ряд с произвольной фрактальной размерностью, то
        // пользуемся функцией Generator.genFractDim(H). Но, в этом случае
        // генерация займет около полу минуты. Для вычисления среднего по 
        // по 1000 испытаний, времени требуется очень много. Рекомендуется снизить
        // количество испытаний до 50.
        static void Main(string[] args)
        {
            Console.WriteLine("Выберите действие:");
            Console.WriteLine("1) Тестирование R/S и Rand R/S на ряде с произвольной фрактальной размерностью");
            Console.WriteLine("2) Тестирование ряда с h=0.5 на ступенчатой помехе");
            Console.WriteLine("3) Тестирование ряда с h=0.5 на помехе из ряда ступенек");
            Console.WriteLine("4) Тестирование ряда с h=0.5 на помехе с фрактальной вставкой");
            int testNum = Convert.ToInt32(Console.ReadLine());
            switch (testNum)
            {
                case 1:
                    Console.Write("Введите показатель Херста (дробная часть через зпт). H = ");
                    double testH = Convert.ToDouble(Console.ReadLine());
                    double[] a = Generator.genFractDim(testH);
                    Console.Write("R/S = ");
                    Console.WriteLine(FractalDimension.calcHurstEst(a, 1024));
                    Console.Write("Rand R/S = ");
                    Console.WriteLine(FractalDimension.calcHurstEstRandomized(a, 1024, 5));
                    break;
                case 2:
                    double[] x = new double[1000];
                    double[] y = new double[1000];
                    Console.Write("Высота ступени = ");
                    double height = Convert.ToDouble(Console.ReadLine());
                    for (int i = 0; i < 1000; i++)
                    {
                        Random rand = new Random();
                        a = Generator.gen(rand, 1, height, null);
                        double H = FractalDimension.calcHurstEst(a, 1024);
                        double H2 = FractalDimension.calcHurstEstRandomized(a, 1024, 5);
                        x[i] = H;
                        y[i] = H2;
                    }
                    // считаем среднее
                    double sumx = 0.0;
                    double sumy = 0.0;
                    for (int i = 0; i < 1000; i++)
                    {
                        sumx = sumx + x[i];
                        sumy = sumy + y[i];
                    }
                    double meanx = sumx / 1000.0;
                    double meany = sumy / 1000.0;
                    // выводим среднее на экран
                    Console.WriteLine("");
                    Console.WriteLine("Без помехи результат:");
                    Console.WriteLine("R/S = 0,55");
                    Console.WriteLine("Rand R/S = 0,57");
                    Console.WriteLine("");
                    Console.WriteLine("Cо ступенчатой помехой:");
                    Console.Write("R/S = ");
                    Console.WriteLine(meanx);
                    Console.Write("Rand R/S = ");
                    Console.WriteLine(meany);
                    break;
                case 3:
                    x = new double[1000];
                    y = new double[1000];
                    Console.Write("Высота ступени = ");
                    height = Convert.ToDouble(Console.ReadLine());
                    for (int i = 0; i < 1000; i++)
                    {
                        Random rand = new Random();
                        a = Generator.gen(rand, 2, height, null);
                        double H = FractalDimension.calcHurstEst(a, 1024);
                        double H2 = FractalDimension.calcHurstEstRandomized(a, 1024, 5);
                        x[i] = H;
                        y[i] = H2;
                    }
                    // считаем среднее
                    sumx = 0.0;
                    sumy = 0.0;
                    for (int i = 0; i < 1000; i++)
                    {
                        sumx = sumx + x[i];
                        sumy = sumy + y[i];
                    }
                    meanx = sumx / 1000.0;
                    meany = sumy / 1000.0;
                    // выводим среднее на экран
                    Console.WriteLine("");
                    Console.WriteLine("Без помехи результат:");
                    Console.WriteLine("R/S = 0,55");
                    Console.WriteLine("Rand R/S = 0,57");
                    Console.WriteLine("");
                    Console.WriteLine("C помехой из последовательности ступенек:");
                    Console.Write("R/S = ");
                    Console.WriteLine(meanx);
                    Console.Write("Rand R/S = ");
                    Console.WriteLine(meany);
                    break;
                case 4:
                    x = new double[1000];
                    y = new double[1000];
                    Console.Write("Фрактальная размерность вставки (дробное значение вводится через зпт) = ");
                    testH = Convert.ToDouble(Console.ReadLine());
                    double[] b = Generator.genFractDim(testH);
                    for (int i = 0; i < 1000; i++)
                    {
                        Random rand = new Random();
                        a = Generator.gen(rand, 3, 0, b);
                        double H = FractalDimension.calcHurstEst(a, 1024);
                        double H2 = FractalDimension.calcHurstEstRandomized(a, 1024, 5);
                        x[i] = H;
                        y[i] = H2;
                    }
                    // считаем среднее
                    sumx = 0.0;
                    sumy = 0.0;
                    for (int i = 0; i < 1000; i++)
                    {
                        sumx = sumx + x[i];
                        sumy = sumy + y[i];
                    }
                    meanx = sumx / 1000.0;
                    meany = sumy / 1000.0;
                    // выводим среднее на экран
                    Console.WriteLine("");
                    Console.WriteLine("Без помехи результат:");
                    Console.WriteLine("R/S = 0,55");
                    Console.WriteLine("Rand R/S = 0,57");
                    Console.WriteLine("");
                    Console.WriteLine("C помехой с фрактальной вставкой:");
                    Console.Write("R/S = ");
                    Console.WriteLine(meanx);
                    Console.Write("Rand R/S = ");
                    Console.WriteLine(meany);
                    break;
                default:
                    break;
            }
            Console.ReadLine();
        }
    }
}