SP_LR_3_SALNIKOVA

function SP_LR_3_SALNIKOVA

    pn1 = 0.325, pn2 = 0.325, pm1 = 0.504, y1 = 0.25, y2 = 0.10 %13 var
    iN=100;
    tH = zeros(iN,1);

    num_of_Sn1 = num_of_Sn2 = num_of_Sn3 = num_of_Sm0 = num_of_Sm1 = 0;
    S = {'S_N(1)','S_N(2)','S_N(3)','S_M(0)','S_M(1)'};
    Sch = cell(iN,1);
    Si = zeros(iN,1);
    tL1 = zeros(iN,1); tL2 = zeros(iN,1);
    C = zeros(iN,2);
    tW = zeros(iN,1);
    num_t = zeros(iN,1);
    type_t = cell(iN,1);
    type = cell(2*iN,1);
    tn = zeros(2*iN,1); tu = zeros(2*iN,1);
    Des1 = zeros(2*iN,1); Des2 = zeros(2*iN,1);
    Sch{1} = S{1};
    Si(1) = 1;
    lambda=0.1;
    tL1(1) = exprnd(1/lambda); tL2(1) = -1;
    C(1,1) = 1;
    tW(1) = tL1(1);
    num_t(1) = 1;
    type_t{1} = 'N';
    type{1} = 'N';
    tn(1) = 0; tu(1) = tL1(1);
    TN = tL1(1);
    jN = 1; TM = 0; jM = 0;
    j = 1;
    for i = 2:iN
        r = rand;
        switch type_t{i-1}
            case 'N'
                if r < pn1
                    tH(i) = tW(i-1) + tH(i-1);
                    Sch{i} = S{1};
                    Si(i) = 1;
                    lambda=y1*jN+2*y1*jM+0.1;
                    tL1(i) = exprnd(1/lambda);
                    tL2(i) = -1;
                    Des1(num_t(i-1)) = j+1; Des2(num_t(i-1)) = -1;
                    j = j + 1;
                    type{j} = 'N';
                    tn(j) = tH(i);
                    tu(j) = tH(i) + tL1(i);
                    TN = TN + tL1(i);
                    jN = jN + 1;
                    C(i,1) = C(i-1,1);
                    C(i,2) = C(i-1,2);
                elseif r < pn1 + pn2 %&& r > pn1
                    tH(i) = tW(i-1) + tH(i-1);
                    Sch{i} = S{2};
                    Si(i) = 2;
                    lambda=y1*jN+2*y1*jM+0.1;
                    tL1(i) = exprnd(1/lambda);
                    tL2(i) = exprnd(1/lambda);
                    Des1(num_t(i-1)) = j+1; Des2(num_t(i-1)) = j+2;
                    j = j + 1;
                    type{j} = 'N';
                    tn(j) = tH(i);
                    type{j+1} = 'N';
                    tn(j+1) = tH(i);
                    if tL1(i) < tL2(i)
                        tu(j) = tH(i) + tL1(i);
                        j = j + 1;
                        tu(j) = tH(i) + tL2(i);
                    else
                        tu(j) = tH(i) + tL2(i);
                        j = j + 1;
                        tu(j) = tH(i) + tL1(i);
                    end
                    TN = TN + tL1(i) + tL2(i);
                    jN = jN + 2;
                    C(i,1) = C(i-1,1) + 1;
                    C(i,2) = C(i-1,2);
                else
                    tH(i) = tW(i-1) + tH(i-1);
                    Sch{i} = S{3};
                    Si(i) = 3;
                    lambda=y1*jN+2*y1*jM+0.1;
                    mu=3*y2*jN+y2*jM;
                    tL1(i) = exprnd(1/lambda);
                    tL2(i) = exprnd(1/mu);
                    Des1(num_t(i-1)) = j+1; Des2(num_t(i-1)) = j+2;
                    j = j + 1;
                    if tL1(i) < tL2(i)
                        type{j} = 'N';
                        tn(j) = tH(i);
                        tu(j) = tH(i) + tL1(i);
                        j = j + 1;
                        type{j} = 'M';
                        tn(j) = tH(i);
                        tu(j) = tH(i) + tL2(i);
                    else
                        type{j} = 'M';
                        tn(j) = tH(i);
                        tu(j) = tH(i) + tL2(i);
                        j = j + 1;
                        type{j} = 'N';
                        tn(j) = tH(i);
                        tu(j) = tH(i) + tL1(i);
                    end
                    TN = TN + tL1(i);
                    jN = jN + 1;
                    TM = TM + tL2(i);
                    jM = jM + 1;
                    C(i,1) = C(i-1,1);
                    C(i,2) = C(i-1,2) + 1;
                end
            case 'M'
                if r > pm1
                    tH(i) = tW(i-1) + tH(i-1);
                    Sch{i} = S{4};
                    Si(i) = 4;
                    tL1(i) = -1;
                    tL2(i) = -1;
                    C(i,1) = C(i-1,1);
                    C(i,2) = C(i-1,2) - 1;
                    Des1(num_t(i-1)) = -1; Des2(num_t(i-1)) = -1;

                else
                    tH(i) = tW(i-1) + tH(i-1);
                    Sch{i} = S{5};
                    Si(i) = 5;
                    mu=3*y2*jN+y2*jM;
                    tL1(i) = exprnd(1/mu);
                    tL2(i) = -1;
                    Des1(num_t(i-1)) = j+1; Des2(num_t(i-1)) = -1;
                    j = j + 1;
                    type{j} = 'M';
                    tn(j) = tH(i);
                    tu(j) = tH(i) + tL1(i);
                    TM = TM + tL1(i);
                    jM = jM + 1;
                    C(i,1) = C(i-1,1);
                    C(i,2) = C(i-1,2);
                end
        end

        tw = tu - tH(i);
        q = find(tw == min( tw( tw > 0 ) ), 1);
        tW(i) = tw(q);
        num_t(i) = q;
        type_t{i} = type{q};
    end

    Sn = zeros(1,5);
    for i = 1:5
        Sn(i) = sum( Si == i );
    end
    Sfreq = Sn ./ iN;
    TN = TN / jN;
    if jM > 0
        TM = TM / jM;
    end

    printf('\n\n---TABLE1---\nNomer sobytia:\n')
    print_it(1:iN)
    printf('\nMoment nastuplenia sobytia:\n')
    print_it(tH)
    printf('\nTip sobytia:\n')
    print_is(Sch)
    printf('\nVremya jizni 1:\n')
    print_it(tL1)
    printf('\nVremya jizni 2:\n')
    print_it(tL2)
    printf('\nSostoyanie sistemy posle sobutia:\n')
    for i = 1:length(C)
      printf('\n(%i, %i)', round_it(C(i, 1), 5), round_it(C(i, 2), 5));
    endfor
    printf('\nVremya ojidania do sled. sobytia:\n')
    print_it(tW)
    printf("\nNomer ob'ekta, u kotorogo ranshe konec jizni:\n")
    print_it(num_t)
    printf('\nEgo vid:\n')
    print_is(type_t)

    printf("\n\n---TABLE2---\nNomer ob'ekta:\n")
    print_it(1:2*iN)
    printf("\nVid ob'ekta:\n")
    print_is(type)
    printf("\nMoment poyavlenia ob'ekta:\n")
    print_it(tn)
    printf("\nVremya jizni ob'ekta:\n")
    print_it(tu-tn)
    printf("\nMoment ischeznovenia ob'ekta:\n")
    print_it(tu)
    printf('\nPotomok 1:\n')
    print_it(Des1)
    printf('\nPotomok 2:\n')
    print_it(Des2)

    printf("\nSredneye vremya jizni N-ob'ekta:\n")
    disp(TN);
    printf("\nSredneye vremya jizni M-ob'ekta:\n")
    disp(TM);
    printf("\nChislo proizoshedshih sobytiy kajdogo vida:\n")
    printf("\n Sn1: %i\n Sn2: %i\n Sn3: %i\n Sm0: %i\n Sm1: %i\n", Sn(1), Sn(2), Sn(3), Sn(4), Sn(5));
    printf("\nChislo poyavivshihsya ob'ektov tipa N:\n")
    disp(jN);
    printf("\nChislo poyavivshihsya ob'ektov tipa M:\n")
    disp(jM);
    printf("\nOtnositelnye chastoty:\n")
    disp(Sfreq);
end


function print_it(arr)

  for i = 1:length(arr)
   printf('\n%i', round_it(arr(i), 5));
  endfor

end

% okrugleniye do 5 znakov posle zapyatoy:
function R = round_it(x, k)
  R = 0; R = round(x*10^k)/10^k;
endfunction

function print_is(arr)

  for i = 1:length(arr)
   printf('\n%s', arr{i});
  endfor

end