En espérant être sur le bon forum, j'aimerai obtenir de l'aide pour optimiser mon code matlab dans le but de raccourcir son temps d’exécution.
Le but du code, est d'examiner la sensibilité de la pression artérielle dans l'aorte aux 5 variables suivantes (selon un model d'un circuit analogique):
Cv : Tonus veineux
Rp: Résistance périphérique
Emax : Contractilité maximale du cœur
HR : Rythme cardiaque
Total BV : Volume total de sang
Pour examiner la sensibilité, on nous demande de modifier chaque fois une variable et de laisser les autres constantes; et chaque variable modifiée devient un vecteur composé de 50% à 250% de sa valeur de départ.
Voici mon code :
clc ; clear all ; close all
tic
% Flags:
Plot_flag = 1; % 0 = off , 1 = on
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Part 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Time parameters
HR = 80 ; %* [BPM]
T=(60/HR); %[sec]
dt = 0.01 ; % [sec] step period
Ts=(1/3)*T; % systolic time of the period
Heart_cycles = 20 ; % total heart cycles
N_cycle = round(T/dt); % number of steps per heart cycle
tot_time= Heart_cycles*T;
% Heart parameters
V0 = 15 ; % [ml] V0 for Plv calculation
Emax = 2.0 ; % contractility
N_Systole = round(Ts/dt) ; % number of steps per ventricle Systole
%Elasticity
Ed= 10/(120-V0); % diastolic elaticity- taken from ex4
En = zeros(1,N_cycle);
En(1:N_Systole)=0.5*(1+sin(2*pi*(1:N_Systole)/N_Systole-pi/2)); % given by instructions
E = max(Ed,Emax*En); % combine diastolic and systolic elasticity functions
% Vascular constants:
Ra = 0.1; % arterial resistance
Rp = 1.0; % peripheral resistance
Rv = 0.01; % venous filling resistance
Ca = 2.0; % arterial compliance
Cv = 300.0; % venous compliance
tot_Va= zeros(Heart_cycles,N_cycle);
tot_Vlv= zeros(Heart_cycles,N_cycle);
tot_Vv= zeros(Heart_cycles,N_cycle);
tot_Plv= zeros(Heart_cycles,N_cycle);
tot_Pv=zeros(Heart_cycles,N_cycle);
tot_Pa= zeros(Heart_cycles,N_cycle);
tot_Pao= zeros(Heart_cycles,N_cycle);
tot_Qp= zeros(Heart_cycles,N_cycle);
tot_Qv= zeros(Heart_cycles,N_cycle);
tot_Qlv= zeros(Heart_cycles,N_cycle);
% Initiate variables:
%Volume[ml]
Vlv(1) = 120; % left ventricle
Va(1) = 270; % arteries
Vv(1) = 2700; % veins
TotalBV=Vlv+Va+Vv;
%Pressure[mmHg]
Plv(1) = 0; % left ventricle
Pa(1) = 70; % arterial capacitor
Pv(1) = 9; % venous filling
Pao(1) = 70; % aorta
%Flow[ml\sec]
Qlv(1) = 0; % left ventricle (outflow)
Qp(1) = 0; % peripheral resistance
Qv(1) = 0; % ventricle filling (inflow)
% initialization of continuous variables:
tot_Pao=Pao(1);
%we changed the values of our paramters to be from 50% until 250% in
%function of the original values
template_percenteage=0.5:0.1:2.5;
samples_Cv=Cv*template_percenteage;
samples_Rp=Rp*template_percenteage;
samples_Emax=Emax*template_percenteage;
samples_HR=HR*template_percenteage;
samples_TotalBV=TotalBV*template_percenteage;
%we creatted matrices in goal to stock our mean of each parameters
samples_parameters = [samples_Cv;samples_Rp;samples_Emax;samples_HR;samples_TotalBV];
Average_Pao_Vec_Cv = [];
Average_Pao_Vec_Rp = [];
Average_Pao_Vec_Emax = [];
Average_Pao_Vec_HR = [];
Average_Pao_Vec_TotalBV = [];
for each_parameter = 1:length(samples_parameters)
for sample = 1:length(template_percenteage)
if each_parameter==1
Cv = samples_parameters(1,sample);
elseif each_parameter==2
Rp = samples_parameters(2,sample);
elseif each_parameter==3
Emax = samples_parameters(3,sample);
E = max(Ed,Emax*En);
elseif each_parameter==4
HR = samples_parameters(4,sample);
T=(60/HR); %*[sec]
ts=1/3*T;
N_cycle = floor(T/dt);
N_Systole =round(ts/dt);
En = zeros(1,N_cycle);
En(1:N_Systole)=0.5+0.5*sin(2*pi*(1:N_Systole)/N_Systole-pi/2);
E = max(Ed,Emax*En);
elseif each_parameter==5
Va= template_percenteage(sample)*270;
Vlv = template_percenteage(sample)*120;
Vv = template_percenteage(sample)*2700;
end
for cycleIdx = 1:Heart_cycles % main loop for each heart cycle
for stepInCycle = 2:N_cycle %calculating all variables for each cycle at N points:
Vv(stepInCycle)=Vv(stepInCycle-1)+dt*(Qp(stepInCycle-1)-Qv(stepInCycle-1));
Va(stepInCycle)=Va(stepInCycle-1)+dt*(Qlv(stepInCycle-1)-Qp(stepInCycle-1));
Vlv(stepInCycle)=Vlv(stepInCycle-1)+dt*(Qv(stepInCycle-1)-Qlv(stepInCycle-1));
Pv(stepInCycle) = Vv(stepInCycle)/Cv;
Pa(stepInCycle)=Va(stepInCycle)/Ca;
Plv(stepInCycle)=E(stepInCycle)*(Vlv(stepInCycle)-V0);
Qv(stepInCycle)=max(0,(Pv(stepInCycle)-Plv(stepInCycle))/Rv);
Qp(stepInCycle)=(Pa(stepInCycle)-Pv(stepInCycle))/Rp;
if Plv(stepInCycle)>=Pa(stepInCycle)
Pao(stepInCycle) = Plv(stepInCycle);
Qlv(stepInCycle) = (Pao(stepInCycle)-Pa(stepInCycle))/Ra;
else
Pao(stepInCycle) = Pa(stepInCycle);
Qlv(stepInCycle) = 0;
end
end
% Saving variables from cycle to continuous variables:
tot_Pao = [tot_Pao, Pao(2:N_cycle)];
% Update the initial variables before the next cycle:
Pv(1)=Pv(N_cycle);
Plv(1)=Plv(N_cycle);
Pa(1)=Pa(N_cycle);
Pao(1)=Pao(N_cycle);
Qp(1)=Qp(N_cycle);
Qlv(1)=Qlv(N_cycle);
Qv(1)=Qv(N_cycle);
Va(1)=Va(N_cycle);
Vv(1)=Vv(N_cycle);
Vlv(1)=Vlv(N_cycle);
if cycleIdx == Heart_cycles
if each_parameter == 1
Average_Pao_Vec_Cv = [Average_Pao_Vec_Cv, mean(Pao)];
elseif each_parameter == 2
Average_Pao_Vec_Rp = [Average_Pao_Vec_Rp, mean(Pao)];
elseif each_parameter == 3
Average_Pao_Vec_Emax = [Average_Pao_Vec_Emax, mean(Pao)];
elseif each_parameter == 4
Average_Pao_Vec_HR = [Average_Pao_Vec_HR, mean(Pao)];
elseif each_parameter == 5
Average_Pao_Vec_TotalBV = [Average_Pao_Vec_TotalBV, mean(Pao)];
end
end
end
end
end
if Plot_flag==1
figure(5);
plot(template_percenteage*100, Average_Pao_Vec_Cv/Average_Pao_Vec_Cv(6) ,'b');
hold on
plot(template_percenteage*100, Average_Pao_Vec_Rp/Average_Pao_Vec_Rp(6) ,'r');
hold on
plot(template_percenteage*100, Average_Pao_Vec_Emax/Average_Pao_Vec_Emax(6) ,'g');
hold on
plot(template_percenteage*100, Average_Pao_Vec_HR/Average_Pao_Vec_HR(6) ,'m');
hold on
plot(template_percenteage*100, Average_Pao_Vec_TotalBV/Average_Pao_Vec_TotalBV(6) ,'y');
title('The average Pao as function of the difference in percantage relative to the original value');
xlabel('The difference relative to the original value [%]');
ylabel('The average Pao [mmHg]');
legend('Cv','Rp','Emax','HR','total BV');
end
toc
Le code est assez "scolaire" et comme je débute sur matlab, je ne saisis pas encore toutes les astuces pour éviter toutes ces itérations.
le secret principal de Matlab est ... que tu dois écrire le moins de boucles possible, parce que Matlab déteste les boucles écrites explicitement (comme Python !) ! Au contraire, il faut tout vectoriser, quitte à effectuer des opérations supplémentaires inutiles. En fait, ce que Matlab sait faire le mieux, c'est les opérations matricielles (d'ailleurs Matlab veux dire matrix laboratory). Je m'explique. Si on prend un vecteur \(x\) de taille \(n\) pour lequel tu veux calculer (par exemple) \(y_i = x_{i+1} - x_i\), en C (par exemple), tu écrirais un truc du genre
for(int i=0; i<n-1; ++i){
y[i] = x[i+1] - x[i];
}
alors qu'en Matlab, tu vas simplement écrire
y = x(2:end) - x(1:end-1);
Si tu écris explicitement la boucle, tu vas mettre beaucoup plus de temps (fais-le !). Une première manière de faire les choses est donc de se représenter ton problème sous forme tensorielle. Par exemple, si tu as deux paramètres \(x, y\) à étudier, et que tu regardes une sortie \(f \equiv f(x,y) \), que tu as discrétisé tes intervalles de paramètres en \(x_i\) de taille \(n_x\) et \(y_j\) de taille \(n_y\), ton résultat est en fait une matrice \(F = (f_{ij} = f(x_i, y_j))_{ij}\) de taille \(n_i\times n_j\). Du coup, au lieu de boucler sur chacun de tes couples de paramètres, tu dois effectuer tes calculs directement sur ta matrice. Pour t'aider, tu peux regarder la documentation pour la fonction ndgrid et les exemples.
Ce que je te conseille pour l'instant, c'est de vectoriser progressivement en analysant l'ordre de tes boucles. Par exemple, peut-être que ta boucle en temps peut être mise à l'extérieur.
Avez-vous entendu parler de Julia ? Laissez-vous tenter ...
× Après avoir cliqué sur "Répondre" vous serez invité à vous connecter pour que votre message soit publié.
× Attention, ce sujet est très ancien. Le déterrer n'est pas forcément approprié. Nous te conseillons de créer un nouveau sujet pour poser ta question.
Avez-vous entendu parler de Julia ? Laissez-vous tenter ...