Modélisation sur MNT : optimisation de code

Bonjour à tous,

Je viens d'écrire un code qui permet de modéliser les transferts hydriques sur un bassin versant (oui je suis en hydrologie ;). Cela fonctionne mais c'est assez lent. J'ai vectorialisé au maximum les calculs, sauf une partie que je n'ai pas réussie. Je ne sais pas d'ailleurs si c'est possible... Je vais essayer d'expliquer clairement :

Pour chaque cellule de mon terrain j'ai calculé dans quelle direction s'écoule l'eau. Une matrice de la même dimension que mon terrain possède cette information. Si je mets cela en 1D pour simplifier :
valeurs = [1, 0, 1];
directions = [2, 0, 2];

Ce que j'attends c'est
valeurs = [0, 2, 0]
c'est à dire que la première valeur soit déplacée dans la seconde cellule, et la dernière valeur également. Et il faut que ce soit sommé !

Quelque chose de ce genre serait génial (mais ne fonctionne pas) :
valeurs2 = valeurs + valeurs(directions)
(ici valeurs2 est une matrice de transfert)
Mais impossible de mettre en argument un vecteur d'indices.
Ma solution actuelle est de traiter au cas par cas. Mais le for... ça prend du temps :s

Merci d'avance pour vos idées ou pistes...
Anthony.

Pour ma part, je ne trouve pas ton exemple très parlant... :?

Salut,

Je crois avoir compris, mais je reconnais qu'il faudrait que tu sois plus rigoureux dans ta description. Typiquement, tu as fait une erreur dans le vecteur direction que tu donnes pour que ce soit cohérent avec ce que tu dis.

Pour répondre à ta question, tu peux déjà découper ton problème en 3 parties, selon les 3 "directions" possibles. Je prendrai cette norme (différente de la tienne) pour être compréhensible par tous :
+1 pour aller à droite
0 pour rester sur place
-1 pour aller à gauche

Ce que tu peux faire, c'est traiter les directions les unes après les autres. Donc pour la direction "vers la droite" par exemple, tu peux repérer les indices des endroits où tu as +1 dans ta matrice :

Code:

---

indices = find(direction == +1)

---

Du coup, tu peux récupérer les éléments que tu dois déplacer vers la droite via

Code:

---

valeurs(indices)

---

Enfin, tu peux créer une nouvelle variable :

Code:

---

1 2	[l, c] = size(valeurs); valeurs_finales = zeros(l+2, c+2);

---

Tu remarqueras que j'ai augmenté la taille de la matrice pour pouvoir gérer les bords en créant une matrice plus grande.

Et, pour chaque direction (ici, vers la droite donc) :

Code:

---

valeurs_finales( 1+1:l+1 , 1+2:c+2 ) = valeurs_finales( 1+1:l+1 , 1+2:c+2 ) + valeurs(indices)

---

Une fois que tu as fait toutes tes directions (y compris les 0 pour dire qu'il n'y a pas de déplacements), alors tu supprimes les bords pour récupérer la matrice à la bonne taille :

Code:

---

valeurs_finales = valeurs_finales( 1+1:l+1 , 1+1:c+1 )

---

Cordialement,

ps: simple curiosité, c'est pour faire quoi ?

Bonjour et merci pour ta réponse.

Oui déjà une petite erreur dans mon post :
déplacement vaut 2 1 2 et non 2 0 2 sinon on a un indice 0.

Tachmou je comprends comment tu procèdes. Du coup cela se rapproche d'une solution que je n'avais pas encore testé mais qui me rebutait un peu : prendre le problème à l'envers... Au lieu de dire où l'on va c'est regarder d'où l'on vient. Mais comme je peux avoir à sommer différentes cases il faut que je fasse le calcul en plusieurs fois. (J'ai déjà essayé dans le sens normal en vectoriel comme la machine fait ça linéairement les derniers calculs fait sur une cellule supplantent ceux fait en premier).

Par exemple pour une matrice :
1 2 3
4 5 6
7 8 9
pour la case 5 je peux avoir des sources en 1, 2 et 3. Je ne peux pas le dessiner mais imaginez des flèches qui convergent de 1, 2, 3 à 5.
Je ne peut donc pas assigner à une case une seule origine.
En regardant le nombre d'origine max par cellule que je peux avoir je devrai pouvoir faire ce calcul moins de 8 fois (je n'ai pas de points de convergence totale).

Je vais explorer cette piste.
J'espérais qu'il y avait dans la syntaxe matlab une méthode que je ne connaisse pas encore qui permette une telle réalisation en une ligne ;) mais ce n'est pas le cas.

Il s'agit d'une modélisation de bassin versant.
J'ai un MNT. Grâce aux différentes altitudes j'en ai déduit les gradients de pente et les directions d'écoulement. A partir d'un point, définit comme exutoire (point de sortie d'un bassin, là ou passe toute l'eau qui ruisselle sur une surface donnée) je définit en "remontant" les écoulements le contour du bassin. Tout ça c'est fait. Ensuite je fais tomber de la pluie. Quelques lois physiques simples (sinon ce sont les équations de Barré de Saint Venant, très vite bien plus chaud) caractérisent mon écoulement. Tout ceci étant en fait pour tracer des particules arrachées aux versants lors d'épisodes de crue.

Voili voilou.
Pas de solution "miracle" donc ;) je vais essayer de creuser tout ça.

Salut,

Ok, je vois où tu veux en venir ^^

Je ne sais pas si on peut faire ça en une fois, je ne suis pas un pro de matlab, loin de là. Mais vu ce que tu donnes et ce que je sais faire, c'est comme ça que j'aurais procédé ;)

Sinon, pour le travail que tu es en train de faire, j'aurais peut-être une remarque idée assez vague en fait, mais qui peut peut-être te donner une autre piste : je n'ai aucune idée de ce qu'elle vaut. Ton projet me rappelle un peu un projet sur lequel je travaillais sur la détection de cachalot en mer, et je devais bosser sur la simulation de la propagation des ondes dans l'eau de mer... et j'avais justement développé un outil de simulation numérique sous matlab pour simuler la propagation d'onde à partir d'un point source dans l'eau de mer... et j'avais procédé de proche en proche en scindant mon espace... ce qui impliquait de lourds calculs --> un peu comme toi.

Après discussion avec des collègues, il s'est avéré qu'il aurait été plus intéressant de mettre mon système complet en équation... et ainsi de résoudre le système complet vectoriellement ! Après réflexion, le gain en calcul aurait été énorme, ça aurait "simplifié" le problème même si cela impliquait de travailler sur un modèle un peu moins "explicite"...

Du coup, mon idée/remarque allait dans le même sens : tu ne peux pas exprimer ton problème sous forme "d'équations différentielles généralisées" que tu pourrais balancer dans un solver directement ? Le gain en temps de calcul pourrait être... énorme et te permettre de simuler du coup des systèmes plus gros et plus complexes. Bref, c'est un boulot monstre que je te propose en remplacement, et difficile... mais qui peut être tellement bénéfique !

Mais si tu ne vois pas trop de quoi je veux parler : oublie ^^

Sinon, si tu continues sur la solution que je t'ai donnée plus haut, pose bien les choses et définies des bases à ton problème : Choisis tes directions, écris les à part, définie les correctement et applique les à ton problème. De même pour tes valeurs. Une fois que c'est fait, tu peux avancer sereinement. Parce que pour moi, ce n'est pas encore clair la manière dont tu procèdes.

Typiquement, une convention utilisée par le passé était :
8 1 2
7 0 3
6 5 4

Ce qui constitue ta "rose des vent". Maintenant, si la "direction" d'une case vaut 4, cela voudra dire que la valeur contenue dans cette case sera rajoutée dans la case en bat à droite dans une nouvelle matrice. Et ainsi de suite. Là, c'est un peu mieux défini par exemple. Mais il faut pousser le bouchon plus loin encore pour que tout soit bien clair. Ce travail de "définition" est primordial je trouve : ça simplifie tout.
Dans tous les cas, bon courage, c'est un projet super intéressant !

Cordialement,

Oui je vois ce que tu veux dire. Plutôt qu'une discrétisation temporelle un moyen d'arriver directement à l'état final. Dans mon cas je doute que cela soit possible, mais surtout je ne veux pas ! Ce sont toutes les étapes intermédiaires qui m'intéressent ;) L'intérêt pour moi est en fait le traçage, c'est à dire de voir où passent et où se concentrent les particules.

J'avais commencé au début à utiliser cette fameuse rose des vent. J'avais déduit tous mes écoulements à partir de ça. Sauf qu'à chaque fois il fallait que j'utilise une fonction de conversion. J'ai trouvé beaucoup plus rapide ma méthode actuelle qui est de mettre pour un MNT mxn les coordonnées dans une matrice mxnx2 avec i j les décallages de cellule genre -1;1.

Du coup j'ai trouvé une manière plus propre que ce que je disais au post précédent. Pour ne pas avoir de télescopage entre les différentes cellules, il faut faire le calcul en 4 fois. Cela reste du calcul matriciel donc très rapide.
Je note dans ton idée le problème des bords que j'aurai pu avoir.
Du coup pour une matrice mxn, je m'intéresse à i+1:m-1, j+1,n-1, et en 4 fois j'additione ce qui se déplace dans une matrice de transfert.
Je ne sais pas comment expliquer cela plus clairement que:
En prenant en 1 les cases que je considère, plutôt que :
111111
111111
111111
111111
111111
111111

je fais
000000
010100
000000
010100
000000
000000

Puis je décale à droite, en bas, puis à gauche. En 4 fois, les cellules de destination ne se recoupent pas ;)
Ca fonctionne dans ma tete, je vais voir si MatLab est d'accord, avec l'avis de Profiler sur les perfos ;)
Je posterai un message (peut-être plus clair :p) pour dire si cela fonctionne. Ce genre de problème peut-être rencontré par d'autres je pense.

Ta ta ta ta... je m'autocorrige encore.
Après avoir testé mon idée ci-dessus, je me suis rendu compte que... et bien non ce n'était pas assez écarté. Bref à ceux qui essayent de me lire (et comprendre) oubliez l'idée ci-dessus. La méthode la plus compacte (que je viens de réaliser) est de prendre le problème à l'envers.

Pour ceux qui rencontreraient le même probleme de transfert sur une grille de terrain. Pour quoi que ce soit. Un déplacement de troupeau de lapins ou autres : Lister pour chaque case les cases alentours qui pointent dessus. Prendre le nombre maximal 'max' de cases alentours qui peuvent pointer dans une seule (ce qui permet de limiter le temps de calcul, aussi vectoriel soit-il). Ainsi pour une matrice mxn, je crée une matrice chemins = ones(m x n x max). Pour chaque cellule je liste dans les différents niveaux les index des cellules qui me regardent. Bien sûr on ne remplit pas tout. Partout ailleurs il y aura des 1. La cellule 1x1 servira donc de poubelle. Pour moi le bord m'était inutile. Si ce n'est pas le cas, on pourra traiter ce cas individuellement à la fin.
Enfin avec matrice(chemins), niveau par niveau on a les éléments qui se déplacent. On somme les niveau. On soustraie à la matrice d'origine ce qui a bougé. On additionne à la matrice ce qui a bougé, et qui a été déplacé. ;)
(pour les indices utiliser ind2sub et sub2ind).

Mon code ne sera pas plus clair que ma tentative d'explication, il est noyé dans la partie "utile" du modèle. Ceux qui rencontreront le même problème devraient comprendre de quoi je parle. J'espère ;).

Merci tachmou pour tes conseils en tous cas ;)
Sur un run de test, ma fonction est passée de 60 secondes à 2 ou 3 :)
Encore un peu de Profiler pour soigner les fonctions graphiques trop gourmandes et ce sera top 8-)

Ciao !

cs809Jai lu ta réponse et je travaille sur un projet d'écoulement d'eau sur matlab et d'apres tes réponses, j'ai vu que tu t'y connaissais ( j'ai lancé des forums et j'ai pas eu de résultats).
Mon problème c'est que j'arrive pas à modéliser mon terrain avec déja une pente et une hauteur ( la hauteur je ne suis pas sure qu'elle sera utile, ça serait bien de l'avoir).

Merci d'avance pour ta réponse

Citation:

---

Envoyé par nouriy

cs809Jai lu ta réponse et je travaille sur un projet d'écoulement d'eau sur matlab et d'apres tes réponses, j'ai vu que tu t'y connaissais ( j'ai lancé des forums et j'ai pas eu de résultats).
Mon problème c'est que j'arrive pas à modéliser mon terrain avec déja une pente et une hauteur ( la hauteur je ne suis pas sure qu'elle sera utile, ça serait bien de l'avoir).

Merci d'avance pour ta réponse

---

La création d'un Modèle Numérique de Terrain est un problème très différent de ce qui a été évoqué plus haut (peut être que le titre de la discussion est un peu trompeur puisqu'il s'agit non pas d'un MNT directement mais d'un problème de routage de surface). Je pense que tu devrais créer un post en présentant tes entrées (sous quelle forme sont disponibles tes données de ton bassin versant, etc), et le format de sortie que tu souhaites (et proposer des solutions, c'est important ! )

je lance de pas une discussion à ce sujet.