%function SP_LR_2_KHMELEV
%% Вариант 16
%
%p1 = 0.662;
%p2 = 0.208;
%p3 = 0.13;
%q1 = 0;
%q2 = 1;
%end


function SP_LR_1_KHMELEV

  % task 1
  printf('\n\n-------------\n  TASK NUMBER 1\n\n\n');

  P = [0.173, 0.287, 0.540; 0, 0.587, 0.413; 0, 0.159, 0.841];
  n_min = 100;
  delta = [];
  n = 2;

  while n <= n_min

    P(:, :, n) = P(:, :, n-1) * P(:, :, 1); % 3-merny massiv, kotory soderjit vse matricy stepeney ot 1 do n_min
    delta(n) = max(max(abs( P(:, :, n) - P(:, :, n-1) ))); % nahodim maksimalnoe otklonenie sredi vseh veroyatnostey mejdu nastoyacshey matricey i predyduschey

    if delta(n) < 0.00001
      n_min = n;
    endif

    n += 1;
  endwhile


  printf('\n\n %i Matricies:\n\n', n_min);
  for n = 1:n_min

    disp(P(:, :, n));
    printf('\n');

  endfor

  printf('\n\n And deltas from 2 to %i:\n\n', n_min);
  for n = 2:n_min

    printf('   %i\n', myround(delta(n), 5));

  endfor



  % task 2

  % q - stolbets
  % qP = q => qP - q = 0 => q(P - E) = 0, gde E - edinichnaya matrica, prichem v nashem kode eto eye(n)
  % q = (r1, r2, r3)
  % eye(n) - edinichnaya matrica n x n
  printf('\n\n-------------\n  TASK NUMBER 2\n\n\n');

  A = P(:, :, 1)' - eye(3); % vychitaem edinichnuyu matricu
  A = [A, zeros(3, 1)]; % pripisyvaem 4-y stolbets s nulevymi svobodnymi chlenami
  A(3,:) = ones(1, 4); % t.k. vectory lineyno zavisimye, to vycherkivaem poslednyuyu stroku
  % pri etom u nas est' uslovie normirovki: r1 + r2 + r3 = 1 -- stavim ego na mesto vycherknutoy stroki

  A = gauss(A); % reshaem SLAU

  printf('\n\n');
  q = A(:, 4)' % q - naydennoe stacionarnoe raspredelenie
  qP = q*P(:, :, 1) % Prover 04ka:

  if myround(q, 5) == myround(qP, 5)
    disp(' => q = qP')
  else
    disp(' => q =/= qP')
  endif



  % task 3

  % p(0) = ( p1(0), p2(0), p3(0) )
  % p(n) = p(0)*P^n
  % prichem v nashem kode P^n - eto "P(:, :, n)"
  printf('\n\n-------------\n  TASK NUMBER 3\n\n\n');

  states = eye(3); % potomu chto sostoyaniya (1, 0, 0), (0, 1, 0), (0, 0, 1) - eto stolbtsy edinichnoy matricy
  for p0 = states
    printf('\n\n');

    p0 = p0';
    ProbDistr(p0);

  endfor



  % task 4

  % orientiruyus' po kartinke na doske
  printf('\n\n-------------\n  TASK NUMBER 4\n\n\n');

  Seq = [];
  for i = 1:3

    Seq = [Seq; gen(i)];
    printf('\n\n %i. State numbers for i_0 = %i through 100 steps:\n\n', i, i);
    disp(Seq(i, 1:100));

  endfor


  % task 5
  printf('\n\n-------------\n  TASK NUMBER 5\n\n\n');

  R = [];
  for state = 1:3
    for i = 1:3
      R(state, i) = sum(i == Seq(state,:));
    endfor
  endfor

  V = R/1000;
  for i = 1:3

    printf('\n\n Frequencies for i_0 = %i:  (%i, %i, %i)\n', i, R(i, 1), R(i, 2), R(i, 3));
    printf('\n\n Relative frequencies for i_0 = %i:  (%i, %i, %i)\n\n', i, V(i, 1), V(i, 2), V(i, 3));

  endfor


  figure('Name','Relative frequencies');
  x = [1, 2, 3];
  plot(x, V(1, :), '-or', x, V(2, :), '-og', x, V(3, :), '-ob');
    title('  Relative frequencies      red - 1; green - 2, blue - 3'), grid on








  %_____________________________________________________________________________


  function ProbDistr(p0) % task 3
      n = 1;
      m_min = 100;

      while n <= m_min

        p(n, :) = p0 * P(:, :, n);
        delta_(n) = max(abs( p(n, :) - q ));

        if delta_(n) < 0.00001
          m_min = n;
        endif

        n += 1;
      endwhile


      printf('\n\n   Probability distributions through %i steps for (%i, %i, %i):\n\n', m_min, p0(1), p0(2), p0(3));
      for n = 1:m_min

        disp(myround(p(n, :), 5));
        printf('\n');

      endfor


      printf('\n\n And deltas from 0 to %i:\n\n', m_min);
      delta0 = max(abs( p0 - q) );
      delta_ = [delta0, delta_];
      for n = 1:m_min

        printf('   %i\n', myround(delta_(n), 5));

      endfor

      printf('   %i\n', delta_(end)); % pishu otdelno, chtoby  ne okruglyat' do nulya (nastolko malen'kaya delta)

   endfunction



   function I = gen(i) % task 4
      I = [];
      E = 0;

      for n = 1:1000
        E = rand(1);

        if E < P(i, 1, 1)
           i = 1;
         elseif E < P(i, 1, 1) + P(i, 2, 1)
           i = 2;
         else
           i = 3;
         endif

         I(n) = i;
      endfor

   endfunction


end



function R = myround(x, k) % округление до 5 знаков после запятой
  R = 0;

  R = round(x*10^k)/10^k;

endfunction



% for task 4
%        switch i
%
%          case 1
%            if E < P(1, 1, 1)
%              i = 1;
%            elseif E < P(1, 1, 1) + P(1, 2, 1)
%              i = 2;
%            else
%              i = 3;
%            endif
%            I(n) = i;
%
%
%          case 2
%            if E < P(2, 1, 1)
%              i = 1;
%            elseif E < P(2, 1, 1) + P(2, 2, 1)
%              i = 2;
%            else
%              i = 3;
%            endif
%            I(n) = i;
%
%
%          case 3
%            if E < P(3, 1, 1)
%              i = 1;
%            elseif E < P(3, 1, 1) + P(3, 2, 1)
%              i = 2;
%            else
%              i = 3;
%            endif
%            I(n) = i;