
function sm_run

clf;
close all 
clear all;
set(gcf, 'PaperPositionMode','auto','color', 'white');
set(gcf,'PaperPosition',[1.5 3 5 5])
hFig=axes('Position',[0.2 0.4 .6 .6],'Visible','off','Color',[.9 .9 .9],...
         'FontName','times',...
         'FontSize',8,...
         'XColor',[0 0 .0],...
         'YColor',[0 0 .0],...
         'ZColor',[0 0 .0]);



r=49.5;
shift=0;
flor = @(t,s) [ -10*(s(1)-s(2)); r*s(1)-s(2)-s(1).*s(3)+shift; -2.6666*s(3)+s(1).*s(2)];
options = odeset('RelTol',1e-4,'AbsTol',1e-5,'Events',@events);
%options = odeset('Events',@events);

[t,y,tau1,ye1,ie1]=ode45(flor,[0:0.01:1750],[0.001,0.,0.],options);

plot3(y(:,1), y(:,2), y(:,3),'Color', [0.8 0.8 0.8])
hold on
plot3(ye1(:,1), ye1(:,2), ye1(:,3),'.','Color', [0.8 0.0 0.0])
hold on
k1=length(y(:,2));
plot3(y(k1-100:k1,1), y(k1-100:k1,2), y(k1-100:k1,3),'Color', [1.0 0.0 0.0])
hold on

xlabel('x','FontSize',20);
ylabel('y','FontSize',20);
zlabel('z','FontSize',20);
axis on;
axis([-25 25 -30 30 0 55])
axis tight 
view(10,5)

set(gca,'box','off','FontSize',8);
 hFig=axes('Position',[0.2 0.1 .6 .2],'Visible','off','Color',[.9 .9 .9],...
         'FontName','times',...
         'FontSize',8,...
         'XColor',[0 0 .7],...
         'YColor',[0 0 .7],...
         'ZColor',[0 0 .7]);
box on

plot (t,y(:,1),'Color','b','LineWidth',1, 'Color', [0.8 0.8 0.8])
hold on
plot (tau1,ye1(:,1),'.', 'MarkerSize',15,'Color',[250./255  81./255  5./255])
hold on


b = logical(heaviside(ye1(:,1)'));
b2=b+1;
lookup_string1 = '01';
b1 = lookup_string1(b + 1)
 % find changes in vector,
% forcing a change at the beginning and at the end
%changes = find(diff([~b(1); b(:); ~b(end)]) ~= 0)


[TRANS_EST, EMIS_EST] = hmmestimate(b2,b2);
TRANS_EST'
[vec val] = eig(TRANS_EST');


%c=heaviside(ye1(:,1)');
%heaviside(ye1(:,1)');

xlabel('time','FontSize',20);
ylabel('x','FontSize',20);
axis tight
axis on

complex=calc_lz_complexity1(b,'primitive',0);
complex=calc_lz_complexity1(b,'primitive',0)/length(b);
[q1,q2]=calc_lz_complexity1(b,'primitive',0);
q1;
%celldisp(q2);
length(b);

T={'r', 'complex';r,complex};
disp(T)

end


function [value,isterminal,direction] = events(t,y);
if y(1)>0
    direction= [1];
else
    direction= [-1];
end
value= 10*(y(1)-y(2));
% this means x'=0 as x'=-10(x-y) in the Lorenz model 
isterminal=[0];
end 


