MNK

#include <windows.h>
#include <math.h>
#include <random>
#include <fstream>
#include <sstream>
#include <string>
#include <vector>
#include <cmath>
#include <ctime>

/////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////Не интерактивное меню////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////

#define isFunction 3 // выбор функции, 1 = cos(2*pi*x), 2 = pow(x,3)*5.0+pow(x,2)+5.0, 3 = sin(2*pi*x)*x

const double Pi = 3.14;//число пи

const int points = 6;//кол-во точек выборки
double mistake = 0.10; // модуль максимального значения ошибки

int typeOfError = 1;//тип ошибки: равномерная =1  или нормальная =2.

double coef_x, coef_y; // коэффициенты

int power_of_polynom = 5;//полином n-й степени
const double lambda = 0.0003;

std::mt19937 gen(time(0));
std::uniform_real_distribution<> uid(0, 1);

double currentError = 0.;
/////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////
inline double cosinus(double x) {
    return cos(2 * Pi*x / coef_x);
}
inline double cube(double x) {//не подогнал
    double cock = pow(x / coef_x, 3)*5.0 + pow(x / coef_x, 2) + 5.0 / coef_x;
    return cock;
}
inline double sinus(double x) {
    double cock = sin(2.0 * Pi*x / coef_x)*x / coef_x;
    return cock;
}

double WhatTheFunction(int f, double coef_y, int i)
{
    if (f == 1)
        return coef_y * cosinus(i);
    if (f == 2)
        return coef_y * cube(i);
    if (f == 3)
        return coef_y * sinus(i);
}


//теперь функция, которая рисует полином этот
inline double polynom(double vec_sol[], int sizeOfMatrix, double coef_x, double x) {
    double polynom = 0;
    for (int i = 0; i < sizeOfMatrix; i++) {
        polynom += vec_sol[i] * pow(x, i) ;
    }
    return polynom;
}

double* Gauss(double** matrix, int n) {
    //Метод Гаусса
    //Прямой ход, приведение к верхнетреугольному виду
    int i, j, k;
    double tmp;
    double* vec_sol = new double[n];

    for (i = 0; i<n; i++)
    {
        tmp = matrix[i][i];//Элемент на главной диагонали
        for (j = n + 1; j >= i; j--)//приведение элемента главной диагонали к 1, деление элементов строки
            matrix[i][j] /= tmp;
        for (j = i + 1; j<n; j++)
        {
            tmp = matrix[j][i];//зануление столбцов под главной диагональю
            for (k = n; k >= i; k--)
                matrix[j][k] -= tmp * matrix[i][k];
        }
    }
    //имеем верхнетреугольную матрицу с правым столбиком
    /*обратный ход*/
    vec_sol[n - 1] = matrix[n - 1][n];
    for (i = n - 2; i >= 0; i--)
    {
        vec_sol[i] = matrix[i][n];
        for (j = i + 1; j<n; j++)
            vec_sol[i] -= matrix[i][j] * vec_sol[j];
    }
    return vec_sol;
}

double NormalErr(double mistake) {
    double ksi = 0, randomNumber;
    for (int i = 0; i< 12; i++) {
        randomNumber = uid(gen);
        ksi += randomNumber;//сумма 12 случайно распределённых числел
    }
    double etha = ksi - 6;//Etha принадлежит отрезку от 0 до 1
    randomNumber = uid(gen);
    double randomValue = mistake * (2 * randomNumber - 1);
    return etha * randomValue;
}

void logging(int sizeOfMatrix, std::ofstream& loggi, double** matrix) {
    for (int i = 0; i < sizeOfMatrix; i++) {
        for (int j = 0; j < sizeOfMatrix + 1; j++) {
            loggi << matrix[i][j] << " ";
        }
        loggi << std::endl;
    }
    loggi << std::endl;
}

LRESULT CALLBACK WndProc(HWND hwnd, UINT iMsg, WPARAM wParam, LPARAM lParam)
{
    static int cxClient, cyClient;
    HDC hdc;
    PAINTSTRUCT ps;
    HPEN hPen, hPenOld;
    HBRUSH hBrush, hBrushOld;

    switch (iMsg)
    {

    case WM_CREATE:
        break;

    case WM_SIZE:
        cxClient = LOWORD(lParam);//ширина
        cyClient = HIWORD(lParam);//высота
        break;

    case WM_PAINT:
    {
        std::ofstream loggi;// для сохранения матриц в файл
        loggi.open("loggi.txt");

        double *xx = new double[points];
        double *yy = new double[points];
        double *tt = new double[points];

        hdc = BeginPaint(hwnd, &ps);//начало рисования

        SetMapMode(hdc, MM_ISOTROPIC);// режим рисования с единообразными осями
        SetWindowExtEx(hdc, cxClient, cyClient, 0);//установили размеры
        SetViewportExtEx(hdc, cxClient / 2, -cyClient / 2, 0);
        SetViewportOrgEx(hdc, cxClient / 2, cyClient / 2, 0);

        MoveToEx(hdc, -cxClient, 0, nullptr);//переместили курсор
        LineTo(hdc, cxClient, 0);//нарисовали одну ось
        MoveToEx(hdc, 0, -cyClient, nullptr);
        LineTo(hdc, 0, cyClient);//нарисовали другую


        hPen = CreatePen(PS_SOLID, 0, RGB(255, 0, 0));
        SelectObject(hdc, hPen);

        coef_x = cxClient;
        coef_y = cyClient;// / 2;

        MoveToEx(hdc, -cxClient, static_cast<int>(WhatTheFunction(isFunction, coef_y, -coef_x)), nullptr);//устанавливаем начальную позицию

        for (int i = -cxClient; i < cxClient; i++)//рисуем график функции
        {
            double y = static_cast<int>(WhatTheFunction(isFunction, coef_y, i));
            LineTo(hdc, i, (int)y);
        }

        for (int i = 0; i < points; i++)
        {
            // вычисление равномерного распределённого вектора X
            //генерируем точки для метода
            xx[i] = uid(gen)*coef_x * 2 - coef_x;// генерация элемента выборки, из отрезка [0,1] -> [-cxClient,cxClient]
            yy[i] = WhatTheFunction(isFunction, coef_y, xx[i]);

            //double mis = 0;// ошибка измерений

            //находим t[i] в зависимости от вида распределения ошибки
            if (typeOfError == 1)//имеет равномерное распределение
            {
                currentError = uid(gen) * 2 * mistake*coef_y - mistake*coef_y;
                //домножение на коэффициент (coef_y) оправдано
                tt[i] = yy[i] + currentError;
            }
            if (typeOfError == 2)//имеет нормальное распределение
            {
                currentError = NormalErr(mistake*coef_y);
                tt[i] = yy[i] + currentError;
            }

            hPen = CreatePen(PS_SOLID, 0, RGB(0, 0, 255));
            hPenOld = (HPEN)SelectObject(hdc, hPen);
            hBrush = CreateSolidBrush(RGB(0, 0, 255));
            hBrushOld = (HBRUSH)SelectObject(hdc, hBrush);
            Ellipse(hdc, xx[i] - 5, tt[i] - 5, xx[i] + 5, tt[i] + 5);

            //  xx[i] /= coef_x;
            //  yy[i] /= coef_y;

            SelectObject(hdc, hBrushOld);
            DeleteObject(hBrush);

            SelectObject(hdc, hPenOld);
            DeleteObject(hPen);
        }


        const int sizeOfMatrix = power_of_polynom + 1;
        double** matrix = new double*[sizeOfMatrix];
        for (int i = 0; i < sizeOfMatrix; i++) {
            matrix[i] = new double[sizeOfMatrix + 1];
        }

        //слау такого-то размера подготавливаем к работе
        for (int i = 0; i < sizeOfMatrix; i++) {
            for (int j = 0; j < sizeOfMatrix + 1; j++)
                matrix[i][j] = 0;
        }
        //вычисление всей матрицы, кроме последнего столбца
        for (int i = 0; i < sizeOfMatrix; i++)
        {
            for (int j = 0; j < sizeOfMatrix; j++)
            {
                for (int k = 0; k < points; k++)
                    matrix[i][j] += pow(xx[k], i + j);
            }
        }

        //вычисление последнего столбца (правых частей)
        for (int i = 0; i < sizeOfMatrix; i++) {
            for (int k = 0; k < points; k++)
                matrix[i][sizeOfMatrix] += tt[k] * pow(xx[k], i);
        }
        //логирование
        logging(sizeOfMatrix, loggi, matrix);
        double* vec_sol = Gauss(matrix, sizeOfMatrix);
        logging(sizeOfMatrix, loggi, matrix);

        //регуляризация
        for (int i = 0; i < sizeOfMatrix; i++) {
        matrix[i][i] += lambda;
        }

        double* vec_sol_reg = new double[sizeOfMatrix];
        vec_sol_reg = Gauss(matrix, sizeOfMatrix);

        //Выводим решения
        for (int i = 0; i < sizeOfMatrix; i++) {
        loggi << vec_sol[i] << " ";
        }
        loggi << std::endl;
        for (int i = 0; i < sizeOfMatrix; i++) {
        loggi << vec_sol_reg[i] << " ";
        }
        loggi.close();


        hPen = CreatePen(PS_SOLID, 0, RGB(0, 255, 0));
        hPenOld = (HPEN)SelectObject(hdc, hPen);
        hBrush = CreateSolidBrush(RGB(0, 255, 0));
        hBrushOld = (HBRUSH)SelectObject(hdc, hBrush);

        MoveToEx(hdc, -cxClient, (int)polynom(vec_sol, sizeOfMatrix, coef_x, -cxClient), nullptr);//устанавливаем начальную позицию
        for (int i = -cxClient; i < cxClient; i++)//рисуем график функции
        {
        double y = polynom(vec_sol, sizeOfMatrix, coef_x, i);
        LineTo(hdc, i, y);
        }



        hPen = CreatePen(PS_SOLID, 0, RGB(0, 0, 0));
        hPenOld = (HPEN)SelectObject(hdc, hPen);
        hBrush = CreateSolidBrush(RGB(0, 0, 0));
        hBrushOld = (HBRUSH)SelectObject(hdc, hBrush);
        MoveToEx(hdc, -cxClient, (int)polynom(vec_sol_reg, sizeOfMatrix, coef_x, -cxClient), nullptr);//устанавливаем начальную позицию
        for (int i = -cxClient; i < cxClient; i++)//рисуем график функции
        {
        double y = polynom(vec_sol_reg, sizeOfMatrix, coef_x, i);
        LineTo(hdc, i, y);
        }



        /*for (int i = 0; i < sizeOfMatrix; i++) {
        delete[] matrix[i];
        }
        delete[] matrix; */

        ///////////////////////////////////////////////////////////////////////////

        return 0;
    }
    case WM_DESTROY:
        PostQuitMessage(0);
        return 0;
    }
    return DefWindowProc(hwnd, iMsg, wParam, lParam);
}

int WINAPI WinMain(HINSTANCE hInstance, HINSTANCE hPrevInstance,
    PSTR szCmdLine, int iCmdShow)
{
    static TCHAR szAppName[] = TEXT("SineWave");
    HWND hwnd;
    MSG msg;

    WNDCLASSEX wndclass;
    wndclass.cbSize = sizeof(wndclass);
    wndclass.style = CS_HREDRAW | CS_VREDRAW;
    wndclass.lpfnWndProc = WndProc;
    wndclass.cbClsExtra = 0;
    wndclass.cbWndExtra = 0;
    wndclass.hInstance = hInstance;
    wndclass.hIcon = LoadIcon(nullptr, IDI_APPLICATION);
    wndclass.hCursor = LoadCursor(nullptr, IDC_ARROW);
    wndclass.hbrBackground = (HBRUSH)GetStockObject(WHITE_BRUSH);
    wndclass.lpszMenuName = TEXT("MYMENU");
    wndclass.lpszClassName = (szAppName);
    wndclass.hIconSm = LoadIcon(nullptr, IDI_APPLICATION);
    RegisterClassEx(&wndclass);
    hwnd = CreateWindow(szAppName, TEXT("Метод наименьших квадратов"),
        WS_OVERLAPPEDWINDOW,
        CW_USEDEFAULT, CW_USEDEFAULT,
        1000, 750,
        nullptr, nullptr, hInstance, nullptr);
    ShowWindow(hwnd, iCmdShow);
    UpdateWindow(hwnd);
    while (GetMessage(&msg, nullptr, 0, 0))
    {
        TranslateMessage(&msg);
        DispatchMessage(&msg);
    }
    return msg.wParam;
}