% Laboration i FMM
%% Uppgift 1 & 2
clear all
% close all
load randwebs
load randinits

IVP = 225;        %vilka randweb & randinits som ska köras
A = randwebscell{IVP};
x0 = randinitscell{IVP};
NoNeg = odeset('NonNegative',1);

tmax = 1500;
tspan = [1 tmax];

[t, x] = (ode45(@(t,x) x.*A*x,tspan,x0,NoNeg)) ;

plot(t, x)


n = 0;
for i = 1:10
    if x(end, i) > 10^(-4)
        n= n+1;
    end
end
disp('Number of species alive')
n
%% Uppgift 3
alive = zeros(1,1000);
for i = 1:1000
    i
    A = randwebscell{i};
    X0 = randinitscell{i};
    [t,x] = ode45(@(t,x)x.*A*x,tspan,X0,odeset('NonNegative',1));
    for j = 1:10
        if x(end,j) > 10^-3
            alive(i) = alive(i) + 1;
        end
    end
end
%% Uppgift 4
histogram(alive)


%% uppgift 9 (att oscillera eller inte oscillera, det är frågan)
    % Den här testar bara en i taget, det var så vi jobbade när du var här
    % Se nästa för färdig version


grafNr = 10;
annat = A(grafNr, :)*x';
dx = (x(:,grafNr)).*annat';
figure
plot(t, dx)

if isempty(find(abs(dx(round((0.75*length(dx))):length(dx))) > 10^(-4))) == 1
    disp('Skiten oscillerar inte')
else
    disp('Oscillerar')
end


%% Riktiga uppgift 9, räkna hur många som oscillerar
    %verkar funka, men dubbelkolla gärna

n = 0;
for i = 1:10
    grafNr = i;
    annat = A(grafNr, :)*x';
    dx = (x(:,grafNr)).*annat';

   if isempty(find(abs(dx(round(0.75*(length(dx)))):length(dx))) > 10^(-4)) == 0
       n=n+1;
   end
end

disp('Antal oscillerande: ')
n


%% Uppgift 3 & 4
% Max tidsintervall?
tspan = [1 2500];


alive = zeros(1,1000);
sensitivity = 10^(-3)

for j= 1:1000
    A = randwebscell{j};
    x0 = randinitscell{j};

    [t, x] = (ode45(@(t,x) x.*A*x,tspan,x0,NoNeg)) ;

    for i = 1:10
        if x(end, i) > sensitivity
            alive(j) = alive(j) + 1;
        end
    end
end

figure
histogram(alive)

%Hur många har 10 levande?
tens = find(alive == 10)
nrOf10 = length(tens)

% Uppgift 5, plotta 10 överlevare
A = randwebscell{225};
x0 = randinitscell{225};
[t, x] = (ode45(@(t,x) x.*A*x,tspan,x0,NoNeg)) ;
figure
plot(t, x)
hold on
title('Nr 225, två som återuppstår vid t = 1500')