function SP_LR_2_SALNIKOVA
  % Variant 13

  p = [0, 0.353, 0.647];
  q = [0.59, 0.41];


  %______________________________________task 1:
  states = [ 'GGGRR', 'GGRGR', 'GGRRG', 'GRGGR', 'GRGRG', 'GRRGG', 'RGGGR', 'RGGRG', 'RGRGG', 'RRGGG' ];
  possibleStates = [];
  P = A = zeros(10, 10);

  % togda k n-mu sostoyaniyu mojno obraschat'sya kak  states(1+5*(n-1):5+5*(n-1))
  for N = 1:10
    numbers = [];
    mainState = states(1+5*(N-1):5+5*(N-1));

    for n = 1:10
      state = states(1+5*(n-1):5+5*(n-1));
      flag = g = r = 0;
      pq = 1;

      for m = 1:5
        mS = mainState(m);

        g += mS == 'G';
        r += mS == 'R';
        if mS ~= state(m);
          flag += 1;
  %______________________________________dlya perehodnoy matricy P:

          switch mS
            case 'G'
              pq *= p(g);
            case 'R'
              pq *= q(r);
          endswitch
  %______________________________________

        endif
      endfor

      if flag == 2 && pq ~= 0 % sostoyaniya doljny otlichat'sya na 2 simvola
        numbers = [numbers, n];
        P(N, n) = pq; % stroim perehodnuyu matricu P
        A(N, n) = 1;
      endif
    endfor

    possibleStates = [possibleStates; numbers];
  endfor

  printf('\n\n-------------\n  TASK NUMBER 1\n\n\n');
  possibleStates


  %______________________________________task 2:
  printf('\n\n-------------\n  TASK NUMBER 2\n\n\n');
  P


  %______________________________________task 3:
  xy = [0.9, 3.5; 1.4, 2; 2.5, 3; 2.2, 4; 2.1, 4.8; 1.5, 5; 0.5, 2.25; 0.6, 2.75; 0, 1.75; -0.5, 3.25];
    figure
  graph(triu(P), xy, num2cell(1:10));
  title('i to j');
    figure
  graph(tril(P), xy, num2cell(1:10));
  title('j to i');


  %______________________________________task 4:
  printf('\n\n-------------\n  TASK NUMBER 4\n\n\n');
  incs = incommun(A) % nesuschestvennye sostoyaniya
  if ~isempty(incs)
    fprintf("the chain isn't ergodic\n\n");
  endif


  %______________________________________task 5:
  printf('\n\n-------------\n  TASK NUMBER 5\n\n\n');
  p100 = P^100


  %______________________________________task 6:
  printf('\n\n-------------\n  TASK NUMBER 6\n\n\n');
  r = stationDist(P)


  %______________________________________task 7:
  deltas = [];
  for i = 1:10
    for j = 1:10
      deltas(i, j) = abs( p100(i, j) - r(i) );
    endfor
  endfor

  most_max_delta = max( max(deltas) );
  for j = 1:10
    MAX_DELTAS(j) = max( deltas(:, j) );
  endfor

  printf('\n\n-------------\n  TASK NUMBER 7\n\n\n');
  deltas
  MAX_DELTAS
  most_max_delta

  %______________________________________task 8:
  V = N = [];
  J = 1:10;
  for j = J
    NV = gen(j, 10000);
    N(j) = NV(1);
    V(j, 1:3) = NV(2:4);
  endfor

  printf('\n\n-------------\n  TASK NUMBER 8\n\n\n');
  TABLE = [J', N', r', V(:, 1:3)]
  %print_(TABLE, J);



  %______________________________________used functions:

  function incs = incommun(A) % task 4
    incs = [];

    for k = 1:10
      for i = 1:10
        for j = 1:10

          if i ~= j && A(i,k) + A(k,j) == 2
            A(i,j) = 1; % vnosim v nashu matricu smejnosti informatsiyu o tom, mojno li hot' kak-to iz odnoy vershiny proyti v druguyu
          endif

        endfor
      endfor
    endfor

    for j = 1:10
      indices_of_zeros = find(A(:, j) == 0);
      if length(indices_of_zeros) > 1
        incs(end+1) = j;
      endif
    endfor

  endfunction


  function q = stationDist(P) % task 6
      % 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

      A = P' - eye(10); % vychitaem edinichnuyu matricu
      A = [A, zeros(10, 1)]; % pripisyvaem 4-y stolbets s nulevymi svobodnymi chlenami
      A(10, :) = ones(1, 11); % 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(:, 11)' % 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

  endfunction

  function NV = gen(j, Nj) % task 8
      v = I = [];
      i = j; n = 1;

      while n <= Nj
           I(end+1) = i;

           amount_of_k = 0;
           for k = 1:n
             if I(k) == j
                amount_of_k += 1;
             endif
           endfor
           R = amount_of_k;
           v(n) = R/n;

           [value, min_n] = min( abs(v .- r(j)) );
           if value < 0.001
              Nj = min_n;
           endif

           E = rand(1);
           if  0 <= E && E < P(i,1)
             i = 1;
           elseif E < (P(i,1)+P(i,2))
             i = 2;
           elseif E < (P(i,1)+P(i,2)+P(i,3))
             i = 3;
           elseif E < (P(i,1)+P(i,2)+P(i,3)+P(i,4))
             i = 4;
           elseif E < (P(i,1)+P(i,2)+P(i,3)+P(i,4)+P(i,5))
             i = 5;
           elseif E < (P(i,1)+P(i,2)+P(i,3)+P(i,4)+P(i,5)+P(i,6))
             i = 6;
           elseif E < (P(i,1)+P(i,2)+P(i,3)+P(i,4)+P(i,5)+P(i,6)+P(i,7))
             i = 7;
           elseif E < (P(i,1)+P(i,2)+P(i,3)+P(i,4)+P(i,5)+P(i,6)+P(i,7)+P(i,8))
             i = 8;
           elseif E < (P(i,1)+P(i,2)+P(i,3)+P(i,4)+P(i,5)+P(i,6)+P(i,7)+P(i,8)+P(i,9))
             i = 9;
           else
             i = 10;
           endif

           n += 1;
       endwhile

       NV = [Nj, v(Nj), v(Nj-1), v(Nj-2)];
    endfunction

  %______________________________________konec osnovnoy funkcii
end



%______________________________________other functions:

function A = gauss(A) % task 6 (uf: stationDist)
  len = length(A)-1;

  for j = 1:len-1
    A = noDivOnZero(A, len, j);
    for i = (j+1):len
      d = -A(i, j) / A(j, j);
      A(i, :) += A(j, :) .* d;
    endfor
  endfor

  for i = len:-1:1
    a = A(i, i);
    A(i, :) = A(i, :) ./ a;
    xsum = 0;
    j = len-1;

    while j >= i
      xsum += A(i, j+1) * A(j+1, len+1);
      j -= 1;
    endwhile

    A(i, len+1) -= xsum;
  endfor

endfunction


function A = noDivOnZero(A, len, ij) % task 6 (of: gauss)

  C = [];
  flag = 0;
  k = 0;
  for i = ij:len

    if A(i, ij) == 0
      flag = 1;
    endif

    if flag == 0
      break
    elseif flag == 1
      flag = 2;
    else
      k = i;
      break
    endif

  endfor

  if flag ~= 0
    A([ij, k], :) = A([k, ij], :); % menyayem stroki mestami
  endif

endfunction


function graph(A, xy, labels) % task 3

  [i,j,w] = find(A);
  [ignore, p] = sort(max(i,j));
  i = i(p);
  j = j(p);
  w = w(p);

  X = [ xy(i,1) xy(j,1) repmat(NaN,size(i))]';
  Y = [ xy(i,2) xy(j,2) repmat(NaN,size(i))]';

  plot(X, Y, 'Color', 'c')
    hold on;
  plot(xy(:, 1), xy(:, 2), 'ro');

  if (max(w) - min(w) > 0) % vychislyayem serediny otrezkov:
      cxs = (xy(i, 1)  + xy(j, 1)) ./ 2;
      cys = (xy(i, 2)  + xy(j, 2)) ./ 2;
      weights = cell(size(cxs));

      for iw=1:length(w)
          weights{iw} = num2str(w(iw)); % prevrashayem tsifry v strochnye simvoly
      endfor
      text(cxs, cys, weights, 'FontSize', 14); % pishem vesa u seredin otrezkov
  endif

  if (nargin < 3)
      labels = [];
  endif
  if (~isempty(labels))
      text(xy(:, 1), xy(:, 2), labels, 'Color', 'red', 'FontSize', 24, 'FontWeight', 'bold'); % nanosim nomera vershin
  endif

endfunction

%function print_(TABLE, J) % task 8
%  printf('\nJ\n')
%  for j = J
%   printf('\n%i\n', myround(TABLE(j, 1), 5))
%  endfor
%  printf('\nN\n')
%  for j = J
%   printf('\n%i\n', myround(TABLE(j, 2), 5))
%  endfor
%  printf('\nr\n')
%  for j = J
%   printf('\n%i\n', myround(TABLE(j, 3), 5))
%  endfor
%  printf('\nV\n')
%  for j = J
%   printf('\n%i\n', myround(TABLE(j, 4), 5))
%  endfor
%  printf('\nV1\n')
%  for j = J
%   printf('\n%i\n', myround(TABLE(j, 5), 5))
%  endfor
%  printf('\nV2\n')
%  for j = J
%   printf('\n%i\n', myround(TABLE(j, 6), 5))
%  endfor
%endfunction

function R = myround(x, k) % okrugleniye do 5 znakov posle zapyatoy
  R = 0;

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

endfunction