% -------------------------------------------------------------------------
% ANALISI SIMULATIVA DI ONDE VIAGGIANTI per il modello di viroterapia
%  - Reazione-diffusione 1D per (u,i,z): cellule tumorali non infette,
%    cellule tumorali infette, cellule immunitarie.
%  - Equilibri considerati: E1 (assenza tumore), E2 (solo tumore non infetto),
%    E3 (coesistenza, con formule corrette da errata corrige).
%  - Metodo: differenze finite secondo ordine (Neumann), method-of-lines,
%    integrazione temporale con ODE stiff solver (ode15s).
%  
%
% ISTRUZIONI:
%  1) Imposta i parametri nella sezione "PARAMETRI".
%  2) Esegui il file; i risultati grafici verranno visualizzati e stampati
%     in console. Per esportare figure aggiungi comandi 'exportgraphics'.
%  3) Per fare scansioni parametriche usare funzione wrapper (non inclusa
%     per brevita'), ma la struttura consente di integrarla facilmente.
%
% NOTA:
%   Il modello PDE e':
%      u_t = D_u u_{xx} + f1(u,i,z)
%      i_t = D_i i_{xx} + f2(u,i,z)
%      z_t = D_z z_{xx} + f3(u,i,z)
% -------------------------------------------------------------------------

clear; close all; clc;
fprintf('*** Avvio: analisi di onde viaggianti (viroterapia) ***\n\n');

%% -------------------- PARAMETRI (personalizzare per la tesi) ------------
par.p     = 1.87e-2;    % proliferazione tumore (1/h)
par.q     = 8.34e-3;    % mortalita' cellule infette (1/h)
par.q_z   = 7.50e-3;    % mortalita' cellule immunitarie (1/h)
par.S_z   = 5.00e-2;    % sorgente immunitaria basale (cells/(mm*h))
par.K     = 1.00e4;     % capacita' portante (cells/mm)
par.zeta  = 5.00e-1;    % forza citotossica immunitaria (1/h)
par.alpha = 5.00e-2;    % reclutamento immunitario indotto (1/h)
par.beta  = 3.00e-2;    % tasso di infezione virale (1/h)

% Coefficienti di diffusione (dimensioni spaziali arbitrarie; tarare per i casi)
D_u = 1.0e-3;   % diffusivita'  tumore non infetto (mm^2/h o unita'  adimensionale)
D_i = 1.0e-3;   % diffusivita'  tumore infetto
D_z = 1.0e-1;   % diffusivita'  linfociti (piu' mobile)

% Dominio spaziale
L = 200;        % lunghezza dominio (unita'  spaziali; es. mm)
N = 400;        % numero di nodi spaziali (consigliato: 200-1000 a seconda potenza)
x = linspace(0,L,N)'; dx = x(2)-x(1);

% Tempo di simulazione
tspan = [0, 2000];   % orizzonte temporale (h), adattare in base ai parametri

% Tolleranze numeriche e limiti
tol.zero = 1e-12;
MAXPOP = 1e6 * par.K;   % soglia di "blow-up" numerico (protezione)

% Scelta della condizione al bordo: 'Neumann' o 'Dirichlet'
BCtype = 'Neumann';  % raccomandato: Neumann (zero-flusso)

% Opzioni grafiche
doAnimation = false;   % true per animazione (puo' essere costosa)
plotSnapshots = true;
snapshotTimes = [0, 50, 150, 400, 800, tspan(2)]; % orari in cui salvare snapshot

% Opzioni solver
useJacobian = true;   % se true stimare e passare Jacobiano sparse a ode15s
useNonNegativity = true; % impone soluzioni non negative (puo' rallentare)

% -------------------- controlli preliminari -----------------------------
if N < 50
    warning('Numero di punti spaziali N troppo piccolo: risultati poco accurati.');
end
if dx <= 0
    error('Discretizzazione spaziale errata (dx<=0). Controllare L e N.');
end

%% -------------------- Calcolo equilibri (E1, E2, E3) --------------------
[E1, ok1, msg1] = equilibrium_E1(par, tol);
[E2, ok2, msg2] = equilibrium_E2(par, tol);
[E3, ok3, msg3] = equilibrium_E3(par, tol); 

fprintf('Equilibri calcolati:\n');
print_equilibrium('E1 (assenza tumore)', E1, ok1, msg1);
print_equilibrium('E2 (solo tumore non infetto)', E2, ok2, msg2);
print_equilibrium('E3 (infezione attiva)', E3, ok3, msg3);
fprintf('\n');

%% -------------------- Costruzione della Laplaciana 1D (Neumann) ----------
Lapl = laplacian_1d(N, dx, BCtype);   % matrix N x N approximating d^2/dx^2

% Precompute sparse diffusion operators (3-block diagonal)
I_N = speye(N);
Lu = D_u * Lapl;
Li = D_i * Lapl;
Lz = D_z * Lapl;

%% -------------------- Iniziali (condizioni iniziali spaziali) -----------
% Impostare diversi scenari: 'onda_iniziale' o 'perturbazione_locale'
scenario = 'front_seed'; % 'front_seed' / 'localized_seed' / 'random_small'

switch scenario
    case 'front_seed'
        % semina di tumore su sinistra: front iniziale che si propaga verso destra
        u0 = 0.9 * E2.u * (0.5 * (1 - tanh((x-20)/5))); % transizione su x~20
        i0 = 1e-6 * ones(N,1);                          % quasi zero infezione
        % opzionale: piccola inoculazione virale all'origine x=5
        i0 = i0 + 1e-3 * exp(-((x-5).^2)/(2*2^2));
        z0 = E1.z * ones(N,1);                          % immunita'  basale
    case 'localized_seed'
        % piccolo focolaio al centro
        u0 = 1e-6 * ones(N,1) + 0.5*E2.u * exp(-(x-L/4).^2/(2*(L/40)^2));
        i0 = 0.01 * exp(-(x-L/4).^2/(2*(L/80)^2));
        z0 = E1.z * ones(N,1);
    case 'random_small'
        rng(0);
        u0 = 1e-6 + 1e-4 * rand(N,1);
        i0 = 1e-6 + 1e-5 * rand(N,1);
        z0 = E1.z * ones(N,1) + 1e-3 * randn(N,1);
        z0 = max(z0, 0);
    otherwise
        error('Scenario di inizializzazione "%s" non riconosciuto.', scenario);
end

% concatenazione stato Y = [u; i; z] (lunghezza 3N)
Y0 = [u0; i0; z0];

% imposizione di non-negativita'  immediata
Y0(Y0<0) = 0;

%% -------------------- Rhs e Jacobiano per ODE integratore ----------------
rhs = @(t,Y) rhs_pde(t,Y,par,Lu,Li,Lz,N);

if useJacobian
    jac = @(t,Y) jacobian_pde(t,Y,par,Lapl,D_u,D_i,D_z,N);
else
    jac = [];
end

% Opzioni ODE
opts = odeset('RelTol',1e-7,'AbsTol',1e-9);
if ~isempty(jac), opts = odeset(opts,'Jacobian',jac); end
if useNonNegativity
    opts = odeset(opts,'NonNegative',1:3*N);
end
opts = odeset(opts,'Events',@(t,Y) event_stop_if_blowup(t,Y,MAXPOP));

% Scegli integratore (stiff preferito per RD)
solver = @ode15s;

%% -------------------- Integrazione (method-of-lines) --------------------
fprintf('Integrazione PDE (method-of-lines) su %d punti spaziali e t in [%g, %g]...\n', N, tspan(1), tspan(2));
tic;
try
    sol = solver(rhs, tspan, Y0, opts);
catch ME
    %warning('Integrazione fallita con ode15s: %s\nProvo con ode45 (non-stiff)...', ME.message);
    warning(sprintf('Integrazione fallita con ode15s: %s\nProvo con ode45 (non-stiff)...', ME.message));
    solver = @ode45;
    sol = solver(rhs, tspan, Y0, opts);
end
t_cpu = toc;
fprintf('Integrazione completata in %.2f s (solver %s). Numero di passi registrati: %d\n', t_cpu, func2str(solver), length(sol.x));

% Soluzione: sol.y is (3N) x Nt
T = sol.x;
Yt = sol.y;

% Riordina in matrici (Nt x N)
Nt = numel(T);
U = reshape(Yt(1:N,:), N, Nt)';   % Nt x N
I = reshape(Yt(N+1:2*N,:), N, Nt)';  % Nt x N
Zmat = reshape(Yt(2*N+1:3*N,:), N, Nt)'; % Nt x N

%% -------------------- Analisi onde: front tracking e velocita' ----------
% Scelta della specie da tracciare (tipicamente i o u)
trackSpecies = 'i';  % 'i' o 'u' (in genere i evidenzia bene l'onda di infezione)
switch trackSpecies
    case 'i', fieldMat = I; threshold_rel = 0.05;   % soglia 5% del massimo istantaneo
    case 'u', fieldMat = U; threshold_rel = 0.05;
    otherwise, error('Specie da tracciare non riconosciuta.');
end

% Calcola posizione frontalex(t) come x dove field crosses threshold
frontPositions = nan(1,Nt);
threshold_abs = nan(1,Nt);
for k=1:Nt
    arr = fieldMat(k,:);
    maxval = max(arr);
    threshold_abs(k) = threshold_rel * maxval;
    % Cerca coordinate dove arr crosses threshold in monotone front from left to right
    % prendo primo x where arr >= threshold_abs
    idx = find(arr >= threshold_abs(k), 1, 'last'); % 'last' for rightmost point above? We'll use centroid
    if isempty(idx)
        frontPositions(k) = NaN;
    else
        % linear interpolation to improve front position
        j = idx;
        if j==N
            frontPositions(k) = x(N);
        else
            % find first index where arr >= threshold starting from left
            jleft = find(arr >= threshold_abs(k), 1, 'first');
            if isempty(jleft)
                frontPositions(k) = NaN;
            else
                if jleft==1
                    frontPositions(k) = x(1);
                else
                    % interpolate between jleft-1 and jleft
                    a1 = arr(jleft-1); a2 = arr(jleft);
                    if abs(a2-a1) < eps
                        frontPositions(k) = x(jleft);
                    else
                        s = (threshold_abs(k)-a1)/(a2-a1);
                        frontPositions(k) = x(jleft-1) + s * (x(jleft)-x(jleft-1));
                    end
                end
            end
        end
    end
end

% Pulizia frontPositions: rimuovi NaN iniziali/transitori
valid = ~isnan(frontPositions);
if sum(valid) < 3
    warning('Poche posizioni frontali valide: impossibile stimare velocita''.');
    v_est = NaN;
else
    % considera solo parte non-transitoria: scarta i primi 10% del tempo
    kstart = find(valid,1,'first');
    kend = find(valid,1,'last');
    k0 = kstart + round(0.1*(kend-kstart));
    tik = k0:kend;
    tfit = T(tik)';
    xfit = frontPositions(tik)';
    % regressione lineare (x = v t + b)
    P = polyfit(tfit, xfit, 1);
    v_est = P(1);
    fprintf('Velocita  stimata dell''onda per %s: v circa %.5g (unita  spaziali / unita  temporali)\n', trackSpecies, v_est);
end

%% -------------------- GRAFICI -------------------------------------------
% 1) mappe spazio-temporali (imagesc)
figure('Name','Mappe spazio-temporali','Color','w','Units','normalized','Position',[0.05 0.05 0.9 0.8]);
subplot(3,1,1);
imagesc(x, T, U); axis xy; colorbar; ylabel('t'); xlabel('x'); title('u(x,t) (tumore non infetto)');
hold on;
plot(frontPositions, T, 'w-','LineWidth',1); % front trace
subplot(3,1,2);
imagesc(x, T, I); axis xy; colorbar; ylabel('t'); xlabel('x'); title('i(x,t) (tumore infetto)');
hold on; plot(frontPositions, T, 'w-','LineWidth',1);
subplot(3,1,3);
imagesc(x, T, Zmat); axis xy; colorbar; ylabel('t'); xlabel('x'); title('z(x,t) (immunita)');
hold on; plot(frontPositions, T, 'k-','LineWidth',1);

% 2) snapshot dei profili in alcuni tempi
if plotSnapshots
    figure('Name','Snapshot profili','Color','w','Units','normalized','Position',[0.05 0.05 0.9 0.6]);
    nSnap = numel(snapshotTimes);
    colors = lines(nSnap);
    subplot(3,1,1); hold on; box on; title('u(x,t) : snapshot'); xlabel('x'); ylabel('u');
    subplot(3,1,2); hold on; box on; title('i(x,t) : snapshot'); xlabel('x'); ylabel('i');
    subplot(3,1,3); hold on; box on; title('z(x,t) : snapshot'); xlabel('x'); ylabel('z');
    legendEntries = cell(nSnap,1);
    for k=1:nSnap
        % trova indice tempo piu' prossimo
        [~,idx] = min(abs(T - snapshotTimes(k)));
        subplot(3,1,1); plot(x, U(idx,:), 'Color', colors(k,:), 'LineWidth', 1.4);
        subplot(3,1,2); plot(x, I(idx,:), 'Color', colors(k,:), 'LineWidth', 1.4);
        subplot(3,1,3); plot(x, Zmat(idx,:), 'Color', colors(k,:), 'LineWidth', 1.4);
        legendEntries{k} = sprintf('t=%.1f', T(idx));
    end
    for s=1:3, subplot(3,1,s); legend(legendEntries,'Location','best'); end
end

% 3) evoluzione frontale e velocita' stimata
figure('Name','Posizione frontale e velocita''','Color','w');
plot(T, frontPositions, 'b-'); hold on;
xlabel('t'); ylabel('posizione frontale x_f(t)'); grid on;
if ~isnan(v_est)
    plot(T, v_est*T + frontPositions(find(~isnan(frontPositions),1,'first')), 'r--','LineWidth',1.2);
    legend('x_f(t)','retta di regressione');
    title(sprintf('Stima velocita  v circa %.4g', v_est));
else
    legend('x_f(t)');
    title('Stima velocita  non disponibile');
end

% 4) fase locale (u,i) in un sito rappresentativo (es: punto x = L/4)
[~, idx_site] = min(abs(x - L/4));
figure('Name','Fase locale u-i','Color','w');
plot(U(:,idx_site), I(:,idx_site), 'k-','LineWidth',1.2); grid on;
xlabel('u'); ylabel('i'); title(sprintf('Proiezione di fase (u,i) in x=%.2f', x(idx_site)));
hold on;
if ok1, plot(E1.u,E1.i,'ko','MarkerFaceColor',[.7 .7 .7]); end
if ok2, plot(E2.u,E2.i,'ks','MarkerFaceColor',[.5 .5 .5]); end
if ok3, plot(E3.u,E3.i,'k^','MarkerFaceColor',[.3 .3 .3]); end
legend('traiettoria locale','E1','E2','E3','Location','best');

% 5) se richiesto, animazione (opzionale)
if doAnimation
    hfig = figure('Name','Animazione spazio-temporale','Color','w');
    for k=1:round(Nt/4):Nt
        clf;
        subplot(3,1,1); plot(x, U(k,:),'LineWidth',1.6); ylim([min(U(:)) max(U(:))]); title(sprintf('u t=%.2f',T(k)));
        subplot(3,1,2); plot(x, I(k,:),'LineWidth',1.6); ylim([min(I(:)) max(I(:))]); title(sprintf('i t=%.2f',T(k)));
        subplot(3,1,3); plot(x, Zmat(k,:),'LineWidth',1.6); ylim([min(Zmat(:)) max(Zmat(:))]); title(sprintf('z t=%.2f',T(k)));
        drawnow;
    end
end

fprintf('\n*** Analisi completata. Verificare figure e console per diagnosi. ***\n');

%% ==================== FUNZIONI LOCALI ================================

function [E1, ok, msg] = equilibrium_E1(par, tol)
    if par.q_z <= tol.zero
        ok = false; E1 = struct('u',NaN,'i',NaN,'z',NaN); msg = 'q_z nullo: z* non definito'; return;
    end
    E1.u = 0; E1.i = 0; E1.z = par.S_z/par.q_z;
    ok = (E1.z > tol.zero); msg = ternary(ok, 'esiste (biologico)', 'z* <= 0');
end

function [E2, ok, msg] = equilibrium_E2(par, tol)
    if par.q_z <= tol.zero || par.p <= tol.zero
        ok=false; E2=struct('u',NaN,'i',NaN,'z',NaN); msg='p o q_z non validi'; return;
    end
    E2.u = par.K - (par.zeta*par.S_z)/(par.p*par.q_z);
    E2.i = 0;
    E2.z = par.S_z/par.q_z;
    ok = (E2.u > tol.zero) && (E2.z > tol.zero);
    msg = ternary(ok, 'esiste (biologico)', 'u* <= 0 o z* <= 0 (E2 non fisico)');
end

function [E3, ok, msg] = equilibrium_E3(par, tol)
    
    denom = (par.beta + par.p) * (par.alpha*par.zeta + par.beta*par.q_z);
    if denom <= tol.zero
        ok = false; E3 = struct('u',NaN,'i',NaN,'z',NaN); msg = 'Denominatore nullo o non positivo in i*'; return;
    end
    num = par.K * par.p * par.q_z * (par.beta - par.q) - par.zeta * par.S_z * (par.beta + par.p);
    i_star = num / denom;
    if par.q_z <= tol.zero
        ok=false; E3=struct('u',NaN,'i',NaN,'z',NaN); msg='q_z nullo: z* non definito'; return;
    end
    z_star = par.alpha*i_star / par.q_z + par.S_z / par.q_z;
    if par.beta <= tol.zero
        ok=false; E3=struct('u',NaN,'i',NaN,'z',NaN); msg='beta nullo o troppo piccolo: u* non definito'; return;
    end
    u_star = par.K * par.q / par.beta + (par.zeta/par.beta) * z_star;
    ok = (i_star > tol.zero) && (u_star > tol.zero) && (z_star > tol.zero);
    E3.u = u_star; E3.i = i_star; E3.z = z_star;
    msg = ternary(ok, 'esiste (biologico)', 'Almeno una componente <= 0 (E3 non fisico)');
end

function print_equilibrium(name, E, ok, msg)
    if ok
        fprintf('%s: u*=%.5g, i*=%.5g, z*=%.5g   [%s]\n', name, E.u, E.i, E.z, msg);
    else
        fprintf('%s: NON esiste - %s\n', name, msg);
    end
end

function L = laplacian_1d(N, dx, BCtype)
    % Matrice Laplaciana secondo ordine con condizioni di bordo
    e = ones(N,1);
    L = spdiags([e -2*e e], -1:1, N, N) / dx^2;
    switch lower(BCtype)
        case 'neumann'
            % ghost point reflection: second derivative at 1 = 2*(u2-u1)/dx^2
            L(1,1) = -2/dx^2; L(1,2) = 2/dx^2;
            L(N,N) = -2/dx^2; L(N,N-1) = 2/dx^2;
        case 'dirichlet'
            % Dirichlet zero: L as above but boundaries assume u(0)=u(N+1)=0
            % nothing to change for standard second-order formula if ghost zeros
        otherwise
            error('BCtype non supportato: usare ''Neumann'' o ''Dirichlet''.');
    end
end

function dY = rhs_pde(~, Y, par, Lu, Li, Lz, N)
    % RHS vettoriale per solver ODE: dY/dt = [D_u Lap u + f1; D_i Lap i + f2; D_z Lap z + f3]
    u = Y(1:N);
    i = Y(N+1:2*N);
    z = Y(2*N+1:3*N);

    % Reazioni punto a punto
    f1 = par.p .* u .* (1 - (u + i) ./ par.K) - (par.beta/par.K) .* u .* i - (par.zeta/par.K) .* u .* z;
    f2 = (par.beta/par.K) .* u .* i - par.q .* i - (par.zeta/par.K) .* i .* z;
    f3 = par.alpha .* i - par.q_z .* z + par.S_z;

    % Diffusione (lineare)
    du = Lu * u + f1;
    di = Li * i + f2;
    dz = Lz * z + f3;

    dY = [du; di; dz];
end

function J = jacobian_pde(~, Y, par, Lapl, D_u, D_i, D_z, N)
    % Jacobiano sparso completo (3N x 3N) per ode15s
    u = Y(1:N); i = Y(N+1:2*N); z = Y(2*N+1:3*N);

    % derivate puntuali della reazione
    % d f1 / d u
    df1_du = par.p - (2*par.p .* u)/par.K - (par.p .* i)/par.K - (par.beta .* i)/par.K - (par.zeta .* z)/par.K;
    df1_di = -(par.p + par.beta) .* u / par.K;
    df1_dz = -(par.zeta .* u) / par.K;

    df2_du = (par.beta/par.K) .* i;
    df2_di = (par.beta/par.K) .* u - par.q - (par.zeta/par.K) .* z;
    df2_dz = -(par.zeta/par.K) .* i;

    df3_du = zeros(N,1);
    df3_di = par.alpha * ones(N,1);
    df3_dz = -par.q_z * ones(N,1);

    % Matrici diagonali
    Df11 = spdiags(df1_du, 0, N, N);
    Df12 = spdiags(df1_di, 0, N, N);
    Df13 = spdiags(df1_dz, 0, N, N);

    Df21 = spdiags(df2_du, 0, N, N);
    Df22 = spdiags(df2_di, 0, N, N);
    Df23 = spdiags(df2_dz, 0, N, N);

    Df31 = spdiags(df3_du, 0, N, N);
    Df32 = spdiags(df3_di, 0, N, N);
    Df33 = spdiags(df3_dz, 0, N, N);

    % Sparse Jacobian blocks (3x3 blocks)
    J11 = D_u * Lapl + Df11;
    J12 = Df12;
    J13 = Df13;

    J21 = Df21;
    J22 = D_i * Lapl + Df22;
    J23 = Df23;

    J31 = Df31; % zero
    J32 = Df32;
    J33 = Df33;

    % Assemble sparse 3N x 3N
    J = [J11, J12, J13;
         J21, J22, J23;
         J31, J32, J33];
end

function [value, isterminal, direction] = event_stop_if_blowup(~, Y, MAXPOP)
    % Evento: ferma integrazione se qualche componente supera MAXPOP o diventa NaN
    if any(isnan(Y)) || any(isinf(Y))
        value = 0; isterminal = 1; direction = 0; return;
    end
    value = double(all(Y < MAXPOP));
    isterminal = 1; direction = -1;
end

function out = ternary(cond, a, b)
    if cond, out = a; else out = b; end
end