optimiser mon code (modèle de Potts)
Bonjour à tous,
Je souhaiterais optimiser une partie mon code, j'ai fait pas mal de changements qui ont apporté un gain de temps important (nottament j'ai enlevé tous les randi ou rand de mes boucles) mais là j'ai l'impression d'arriver un peu à mes limites de connaissance en programmation alors si vous avez des idées je prends :wink:
L'idée est de faire une simulation de potts avec deux types de cellules+un milieu extérieur. Chaque cellule est composée de plusieurs sites. L'appartenance de chaque site est designé dans la matrice sg. Et chaque cellule (groupement de sites) est associée à un type cellulaire déterminé dans la matrice type. L'idée à chaque indentation est de tirer un site au hasard dans la grille et de le remplacer par un site voisin (aussi tiré au hasard) si cela est favorable à l'énergie du système (si on diminue l'énergie). Si on diminue pas l'énergie alors la transition s'effectue selon une loi de probabilité.
ça c'est la première partie de mon code, je n'ai pas cherché à l'optimiser pour le moment car c'est pas celle qui bouffe le plus de temps (moins d'une minute), c'est la seconde partie qui est vraiment problématique
sg et vol sont des matrices que je génère dans un autre code avant
Code:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
| JMD=12;
JML=12; %type1-milieu
JDD=2; %type 1-1
JDL=12; %type1-2
JLL=2; %type 2-milieu
T=35; %temperature
%% INITIAL CONDITIONS
XMAX=100; %taille de la grille
YMAX=100;
type=zeros(length(sg),length(sg));
type1=zeros();type2=zeros();
type1(1)=0;type2(1)=0;
for i=1:max(max(sg)) %melange aleatoire de cellules
if rand(1) > 0.5
type1(end+1)=i;
else
type2(end+1)=i;
end
end
type1(1)=[];type2(1)=[];
for i=1:length(sg0)
for j=1:length(sg0)
if sg0(i,j)>0
res=find(type1(:)==sg0(i,j));
if isempty(res)==0
type(i,j)=2;
else
type(i,j)=1;
end
end
end
end
imagesc(sg)
MCS_MAX=1000;
rnd_xy=zeros(1,MCS_MAX*XMAX*YMAX);
rnd_xy=randi([2,length(sg)-1],2,MCS_MAX*XMAX*YMAX);
adjx=randi([-1,1],1,MCS_MAX*XMAX*YMAX);
adjy=randi([-1,1],1,MCS_MAX*XMAX*YMAX);
rnd_proba=rand(1,1,MCS_MAX*XMAX*YMAX);
k=0; |
Voici la seconde partie que j'aimerais optimiser. J'ai deux boucles for, une sur MCS et une autre sur t. J'aimerais à l'avenir lancer une simulation avec MCS_MAX=10^6 et XMAX*YMAX=10^4 et pour le moment sur ma machine je tourne à 0.09s/MCS ce qui donnerait environ 6jours de simulation '^^
J'ai lu qu'on pouvait vectoriser mais chaque indentation dépend du résultat de la précédente et aussi avec des if partout dans la boucle for ça aide pas '^^
Code:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
| for MCS=1:1:MCS_MAX
%MCS
%tic
if MCS>MCS_MAX-100 %annihilation step
T=3;
end
for t=1:1:XMAX*YMAX
%% ----------------- Random choosing of one site ------------------------ I
k=k+1;
rnd_x=rnd_xy(1,k);
rnd_y=rnd_xy(2,k);
c=sg(rnd_x,rnd_y);
%% ----------------- Random choosing of a second site adjacent to the first one ------------------------ II
nb_x=adjx(1,k);
if nb_x==0 %order I neighborhood
nb_y=adjy(1,k);
else
nb_y=0;
end
c_neighbor=sg(rnd_x+nb_x,rnd_y+nb_y);
%% --------------- Calculation of Interfacial Energy ------------------------- IV
if c ~= c_neighbor
e_adhesion=0;new_e_adhesion=0;
for i=-1:1
for j=-1:1
%tic
%toc
if ((logical(i==0)==0) || (logical(j==0)==0)) %si le site (ij) est different du site cible
if (c~= sg(rnd_x+i,rnd_y+j)) %si site 1 est une cellule differente du site (ij)
if type(rnd_x,rnd_y)==type(rnd_x+i,rnd_y+j);
e_adhesion = e_adhesion+JDD;
else
e_adhesion = e_adhesion+JDL;
end
end
if (c_neighbor ~= sg(rnd_x+i,rnd_y+j)) %si le site (ij) est different du site cible
if type(rnd_x+nb_x,rnd_y+nb_y)==type(rnd_x+i,rnd_y+j);
new_e_adhesion = new_e_adhesion+JDD;
else
new_e_adhesion = new_e_adhesion+JDL;
end
end
end
end
end
%% ---------------- Calculation of area contraints ------------------------------ V
e_vol=0;
if (c>0)
e_vol=e_vol+(2*AREA_T-2*vol(c)+1);
end
if (c_neighbor>0)
e_vol=e_vol+(2*vol(c_neighbor)-2*AREA_T+1);
end
%% ------------- Energy difference ----------------------- VI
e_all= new_e_adhesion - e_adhesion + e_vol;
if (e_all >= 0)
prob = exp(-e_all/T);
elseif (e_all <0)
prob = 1;
end
if (prob >= rnd_proba(k))
sg(rnd_x,rnd_y)=c_neighbor;
type(rnd_x,rnd_y)=type(rnd_x+nb_x,rnd_y+nb_y);
if (c>0)
vol(c)=vol(c)-1;
end
if (c_neighbor>0)
vol(c_neighbor)=vol(c_neighbor)+1;
end
end
end
end
%toc
end |
La seule chose que je vois pour le moment serait de paralléliser le calcul des deux if (ligne 44 et 52) qu'est ce que vous en pensez?
Merci d'avance
Bonne fin de journée
Fabien
