MMEEP 2

// MMEEP2.cpp: определяет точку входа для консольного приложения.
//

#include "stdafx.h"
#include <iostream>
#include <armadillo>

using namespace std;
using namespace arma;

struct node
{
    int ind;//виробничих потреб промисловості;
    int agr;//виробничих потреб сільського господарства;
    int con;//споживання
    int gross_output;
    int Calc_Gross_Output() { return gross_output = ind + agr + con; }
};

inline void GrossOutput(node & nd , string str) {
    cout << "Gross output of " + str + ": " << nd.Calc_Gross_Output() << endl;
}

inline mat BuildTechMat(const node & a, const node & b) {
    mat A = zeros(2, 2);
    A(0, 0) = 1.0 * a.ind / a.gross_output;
    A(0, 1) = 1.0 * a.agr / a.gross_output;

    A(1, 0) = 1.0 * b.ind / b.gross_output;
    A(1, 1) = 1.0 * b.agr / b.gross_output;

    return A;
}

void Step_1() {
    node industry, agriculture;

    industry.ind = 1150;
    industry.agr = 920;
    industry.con = 1575;

    agriculture.ind = 820;
    agriculture.agr = 650;
    agriculture.con = 1230;

    GrossOutput(industry, "industry");
    GrossOutput(agriculture, "agriculture");
    cout << endl;

    mat A = BuildTechMat(industry, agriculture);
    cout << "A = \n" << A << endl;

    // Модель міжгалузевої залежності цін (модель врівноважених цін) p = pA + s
    // (p1 p2) = (a11 * p1 + a21 * p2     a12 * p1 + a22 * p2) + (s1 s2)
    //
    // p1 = a11 * p1 + a21 * p2 + s1
    // p2 = a12 * p1 + a22 * p2 + s2
    //
    // (1 - a11) * p1 - a21 * p2 = s1
    // - a12 * p1 + (1 - a22) * p2 = s2
    //
    // Bp = s
    //

    mat B = zeros(2, 2);
    B(0, 0) = 1 - A(0, 0), B(0, 1) = -A(1, 0);
    B(1, 0) = -A(0, 1), B(1, 1) = 1 - A(1, 1);

    vec s = zeros(2);
    s(0) = 0.15;
    s(1) = 0.4;

    vec ans = solve(B, s);

    cout << "Price of industry product: " << ans(0) << endl;
    cout << "Price of agriculture product: " << ans(1) << endl;
    cout << endl;
}

vec Coefficients_of_characteristic_poly(mat33 A) {
    vec res = zeros(4);
    res(0) = -1; // л^3
    res(1) = A(0, 0) + A(1, 1) + A(2, 2); // л^2
    res(2) = A(2, 0) * A(0, 2) + A(2, 1) * A(1, 2) + A(1, 0) * A(0, 1) - A(0, 0) * A(1, 1) - A(0, 0) * A(2, 2) - A(1, 1) * A(2, 2); // л^1
    res(3) = A(0, 0) * A(1, 1) * A(2, 2) + A(1, 0) * A(2, 1) * A(0, 2) + A(0, 1) * A(1, 2) * A(2, 0)
            - A(2, 0) * A(0, 2) * A(1, 1) - A(2, 1) * A(1, 2) * A(0, 0) - A(1, 0) * A(0, 1) * A(2, 2); // л^0
    return res *= -1;
}


void Step_2() {
    mat A = zeros(3, 3);

    A(0, 0) = 0.65, A(0, 1) = 0.1, A(0, 2) = 0.25;
    A(1, 0) = 0.35, A(1, 1) = 0.35, A(1, 2) = 0.2;
    A(2, 0) = 0.1, A(2, 1) = 0.45, A(2, 2) = 0.1;

    vec coef = Coefficients_of_characteristic_poly(A);
    cout << "Coefficients of the characteristic polynomial:\n" << coef << endl;

    cx_vec eigval;
    cx_mat eigvec;
    eig_gen(eigval, eigvec, A);
    cout << "Eigenvalues:\n" << eigval << endl;

    double frob = eigval(0).real();
    int frob_i = 0;
    for (int i = 1; i < eigval.n_rows ; ++i) {
        if (eigval(i).real() > frob) {
            frob = eigval(i).real();
            frob_i = i;
        }
    }
    cout << "Frobenius number: " << frob << endl << endl;

    cout << "Right Frobenius vector:\n" << eigvec.col(frob_i) << endl;

    mat At = A.t();
    eig_gen(eigval, eigvec, At);
    frob = eigval(0).real();
    frob_i = 0;
    for (int i = 1; i < eigval.n_rows; ++i) {
        if (eigval(i).real() > frob) {
            frob = eigval(i).real();
            frob_i = i;
        }
    }
    cout << "Left Frobenius vector:\n" << eigvec.col(frob_i) << endl;

    mat B = inv(eye(A.n_rows, A.n_cols) - A);
    cout << "B:\n" << B << endl;



    mat C = eye(A.n_rows, A.n_cols) - B;
    mat P = eye(A.n_rows, A.n_cols);
    for (int i = 1; ; ++i, P *= A, C += P) {
        if (i > 100000){
            cout << "The series is divergent (N <= 100000)\n";
            break;
        }
        if (-C.min() < 0.01) {
            cout << "The series is convergence for N = " << i << endl;
            cout << "B - E - A - ... = " << endl << -C << endl;
            break;
        }
    }

    vec y = zeros(3);
    y(0) = 900, y(1) = 1500, y(2) = 1000;
    vec x = solve(eye(A.n_rows, A.n_cols) - A, y);
    cout << "(E - A)x = y" << endl << "x = " << endl;
    cout << x << endl;
    cout << (x.min() > 0 ? "Matrix A is productive" : "Matrix A is non-productive") << endl;
}

int main(int argc, const char **argv) {
    Step_1();
    Step_2();

    system("pause");
    return 0;
}