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

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 : Sélectionner tout - Visualiser dans une fenêtre à part
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 : Sélectionner tout - Visualiser dans une fenêtre à part
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