%some examples for using the HOMC MATLAB package
%choose 1 through 5 to run the examples or choose 6 to quit
%for more detail and additional examples, see the following:
%J.Xu, HOMC: A MATLAB Package for Higher Order Markov Chain
clc
close all
while 1
    clear
    fprintf('Some Examples for Using HOMC MATLAB Package:\n\n')
    prompt1='Enter the item number from the list above or enter 6 to quit: ';
    fprintf('1. k-step Transition Probabilities\n')
    fprintf('2. Reduced First Order Chain\n')
    fprintf('3. Limiting Probability Distribution\n')
    fprintf('4. Ever-reaching Probabilities\n')
    fprintf('5. Mean First Passage Times\n\n')
    k=input(prompt1);
    fprintf('\n')
    if k<=3
        P(:,:,1)=[1/2 0 0 0; 1/2 0 1 0; 0 1 0 1; 0 0 0 0];
        P(:,:,2)=[0 0 1/2 1; 0 1/2 0 0; 1/2 1/2 0 0; 1/2 0 1/2 0];
        P(:,:,3)=[0 1 0 1; 1 0 1/2 0; 0 0 1/2 0; 0 0 0 0];
        P(:,:,4)=[0 0 0 0; 1 1 1 0; 0 0 0 1/2; 0 0 0 1/2];
    end
    if k==1
        fprintf('Let P below be the transition tensor of a second order chain:\n')
        P
        fprintf('To find the 10-step transition probabilities, we use bpow(P,10) to obtain:\n')
        fprintf('Press any key to continue...\n')
        pause
        bpow(P,10)
        fprintf('Press any key to continue...\n')
        pause
    elseif k==2
        fprintf('Let P below be the transition tensor of a second order chain:\n')
        P
        fprintf('To construct the transition matrix of the reduced first order chain, we use rcmat(P) to obtain:\n')
        fprintf('Press any key to continue...\n')
        pause
        rcmat(P)
        fprintf('Press any key to continue...\n')
        pause
    elseif k==3
        fprintf('Let P below be the transition tensor of a regular second order chain:\n')
        P
        fprintf('To find the limiting probability distribution, we first use Q=rcmat(P) to construct the reduced first order chain.\n')
        Q=rcmat(P);
        fprintf('Next, we use [v,d]=eigs(Q,1) to obtain an eigenvector of Q corresponding to the dominant eigenvalue 1.\n')
        [v,d]=eigs(Q,1);
        fprintf('Finally, we use t2mat(eyet(size(P)))*v/sum(v) to obtain the following limiting probability distribution:\n')
        fprintf('Press any key to continue...\n')
        pause
        t2mat(eyet(size(P)))*v/sum(v)
        fprintf('Press any key to continue...\n')
        pause
    elseif k==4
        P(:,:,1)=[1/2 1/3 1/2; 1/2 1/3 0; 0 1/3 1/2];
        P(:,:,2)=[1 0 1/2; 0 1 1/2; 0 0 0];
        P(:,:,3)=[0 0 1/2; 1 0 1/2; 0 1 0];
        fprintf('Let P below be the transition tensor of a second order chain:\n')
        P
        fprintf('To find the ever-reaching probabilities, we use F=erp(P,tol) with tol=1e-8 to obtain in long format:\n')
        fprintf('Press any key to continue...\n')
        pause
        F=erp(P,1e-8);
        fprintf('\n')
        for j=1:size(F,3)
            fprintf('F(:,:,%g)=\n\n',j)
            fprintf([repmat('%20.15f',1,size(F,3)) '\n'],F(:,:,j)')
            fprintf('\n')
        end
        fprintf('Press any key to continue...\n')
        pause
    elseif k==5
        P(:,:,1)=[0 1/2 0; 1 0 1; 0 1/2 0];
        P(:,:,2)=[0 1/2 0; 1 0 1; 0 1/2 0];
        P(:,:,3)=[0 1/2 0; 1 0 1; 0 1/2 0];
        prompt2='Enter 1 for the direct method, and 2 for the iterative method: ';
        fprintf('Let P below be the transition tensor of an ergodic second order chain:\n')
        P
        i=input(prompt2);
        if i==1
            fprintf('To find the mean first passage times with the direct method, we use mfptd(P) to obtain:\n')
            fprintf('Press any key to continue...\n')
            pause
            mfptd(P)
        elseif i==2
            fprintf('To find the mean first passage times with the iterative method, we use mu=mfpti(P,mu0,tol),\n')
            fprintf('where mu0=ones(size(P)) and tol=1e-6, to obtain in long format:\n')
            fprintf('Press any key to continue...\n')
            pause
            mu=mfpti(P,ones(size(P)),1e-6);
            fprintf('\n')
            for j=1:size(mu,3)
                fprintf('mu(:,:,%g)=\n\n',j)
                fprintf([repmat('%20.15f',1,size(mu,3)) '\n'],mu(:,:,j)')
                fprintf('\n')
            end
        end
        fprintf('Press any key to continue...\n')
        pause
    elseif k==6
        break
    end
    fprintf('\n')
end
    


