Skip to content

Interface for an external routing model - #500

Closed
maximejay wants to merge 23 commits into
DassHydro:mainfrom
maximejay:interface_gamma_routing_model
Closed

maximejay wants to merge 23 commits into
DassHydro:mainfrom
maximejay:interface_gamma_routing_model

Conversation

@maximejay

Copy link
Copy Markdown
Collaborator

Hello,
Here is a PR to build an interface with external routing model.

One can now define a routing_module as "zero" like the snow_module.
When this routing_module is set to zero, a new variable in response%qt is allocated to (mesh%nac, setup%ntime_step) (otherwise it is allocated to (1,1) and set to -99.). Then Smash will return the elementary discharge at every active cells in response%qt.
response%qt. may be used as input of any routing model.

Moreover, a new differentiation rule has been added: base_forward_run(parameters.control.x)(output.response.qt)
This produce the gradient of qt compare to the control vector. This gradient can be computed using forward_run_b0 (mw_forward.f90) and bound to an other gradients from the external routing model (with differentiation rule like routing_forward_run(cost)(output.response.qt)) <=> routing_forward_run(cost)(input_qt)) ).

The baseline has been regenerated.

Everything are ok except 2 tests, but the diff are very small.
For unknown reason 2 variables has some slight modification:
- custom_bayesian_optimize.zero-gr4-lr.custom_set_1.sim_q : mean_diff = 0.00424
- optimize.zero-gr6-lr.distributed.sim_q : mean_diff = 0.00404

May be we need to document that also.

Thanks

maximejay and others added 10 commits December 17, 2025 18:01
Extract qt in output%response%q if zeros routing module
…xternal routing model such as Gamma.

Add a second rule for the differentiation of Smash: output%response%qt / parametersDT
Simplify _constant.py to handle these names correctly
For unknown reason 2 variable has some slight modification:
- custom_bayesian_optimize.zero-gr4-lr.custom_set_1.sim_q : mean_diff = 0.00424
- optimize.zero-gr6-lr.distributed.sim_q : mean_diff = 0.00404
@maximejay
maximejay requested a review from inoelloc January 13, 2026 17:16
@maximejay maximejay added the enhancement New feature or request label Jan 13, 2026

@inoelloc inoelloc left a comment

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Voila une première passe. Je pense qu'on peut rediscuter de l'implementation pour récupérer les gradients par rapport aux débits pour etre plus flexible. De cette question de options_b et options_d a nettoyer et le truc de créer une nouvelle subroutine à différencier pour eviter le b0 et d0.

Il doit y avoir un truc avec ta version de ruff ou je ne sais quoi mais j'ai l'impression qu'on a pas la meme version pour faire le formattage et le check. Essaye peut etre un pip install --upgrade ruff avant de refaire un make format et check

Comment thread smash/_constant.py Outdated
Comment on lines +128 to +132
# ~ MODULE_RR_PARAMETERS = dict(
# ~ **SNOW_MODULE_RR_PARAMETERS,
# ~ **HYDROLOGICAL_MODULE_RR_PARAMETERS,
# ~ **ROUTING_MODULE_RR_PARAMETERS,
# ~ )

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Supprimer

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Super merci,
j'ai fais un un upgrade de ruff, mais ca n'a rien changé ...

Comment thread smash/_constant.py Outdated
Comment thread smash/_constant.py Outdated
Comment on lines +49 to +58
!~ When conditionning this allocatation, tapenade force
!~ its value to zeros before calling SIMULATION_B...
if (setup%routing_module == "zero") then
allocate (this%qt(mesh%nac, setup%ntime_step))
this%qt = -99._sp
else
!save memory
allocate (this%qt(1, 1))
this%qt = -99._sp
end if

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Est ce que finalement ca serait pas plus simple d'avoir une variable qui s'appelle q_domain. Comme ca, on pourrait elargir le code en disant qu'on peut récupérer les gradients du débit routé de smash et avec le routage zero, on retombe sur le debit elementaire ?

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Et est ce que ca serait pas possible d'allouer par défaut à (1, 1) quand on initialise le modèle. Puis quand on lance un run si on demande le gradient par rapport aux debits, on realloue plutot que routing zero ? Désolé, j'essaye de voir jusqu'ou on peut aller pour etendre ce couplage smash par rapport à la dernière reu a Aix.

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Alors oui, on peut renommer la variable en qdomain, mais je trouve que ca porte à confusion avec qdomain du type derivé mwd_return.f90. En plus cette variable sera visibles par les utilisateurs. Sinon on l'appelle qac (Q Active Cells) ?
A la place d'allouer en fonction du module de routage zeros, on ajoute une ou deux option dans setup:

return_grad_qe = Bool
return_grad_q = Bool

ou une seule variable:

return_opt_grad = "zero" | qe" |  "q" (default="zero") Ca donne quoi None quand ca passe en fortran ?

On alloue qac en fonction des variables précédente:
allocate (this%qac(mesh%nac, setup%ntime_step))

ensuite dans mw_forward on test (ajuster les tests en fonction du choix de setup):

 if (setup%return_grad_qe == .true.) then
            do i = 1, mesh%nac
                output%response%qac(i, time_step) = checkpoint_variable%ac_qtz(i, setup%nqz)
            end do
  end if
 if (setup%return_grad_q == .true.) then
            do i = 1, mesh%nac
                output%response%qac(i, time_step) = checkpoint_variable%ac_qz(i, setup%nqz)
            end do
  end if

ll faut aussi vérifier dans Python dans model/_standardize.py que soit return_grad_q ou return_grad_qe soit True et pas les deux. Donc peut être plus simple de mettre une seule variable dans setup: return_opt_grad

Qu'est ce que tu penses de cette option ?

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

J'ai opté pour la solution avec dans setup la nouvelle variable: return_opt_grad =["none" | "qe" | "q"] !
C'est pas mal, tu me dis...

du coup on peut avoir les grads des param par rapport aux débits élémentaires ou aux débits sur le domaine.

Comment on lines +56 to +57
type(OptionsDT) :: option_d
option_d = options

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Je viens de me rendre compte que ca remonte à un commit de 2023, l'ajout de options_d et options_b aux routines générées par Tapenade .. Je ne sais pas du tout pourquoi ca ne plante pas .. Mais mieux vaut tard que jamais. Par contre, l'idée de la fonction dans mw_forward.f90 et de juste wrapper à l'identique la subroutine tapenade. Je prefererais que options_d et options_b soient passés en argument de cette fonction. On modifiera dans Python pour ajouter une copy d'options dans wrap_forward_run

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oui ok, mais j'ai eu un pb avec la copy de ces types dérivés en Python je crois... C'est pour cela que j'avais choisi cette solution. Je vais retenter.

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ok

Comment thread smash/fcore/meson.build
c_sources, f90_sources, f90wrap_f90_sources, f2py_f90wrap_sources,
c_args: ['-O3', c_ignore_warnings],
fortran_args: ['-O3', '-march=native', '-cpp', fortran_ignore_warnings],
#fortran_args: ['-Wall', '-Wextra', '-fmax-errors=1', '-cpp', '-g', '-fcheck=all', '-fbacktrace', fortran_ignore_warnings],

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Supprimer

Suggested change
#fortran_args: ['-Wall', '-Wextra', '-fmax-errors=1', '-cpp', '-g', '-fcheck=all', '-fbacktrace', fortran_ignore_warnings],

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ok

Comment thread smash/tests/generate_baseline.py Outdated
Comment thread smash/tests/generate_baseline.py Outdated
Comment thread smash/tests/generate_baseline.py Outdated
Comment thread tapenade/generate_tapenade.py Outdated
Comment on lines +81 to +82
"-head",
r"base_forward_run(parameters.control.x)\(output.response.qt)",

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Je me demande si ca serait pas plus simple de créer une nouvelle routine, je sais pas : base_forward_run_q et on met la regle de diff sur cette subroutine. Ca permettra d'tre plus clair, de pas avoir un base_forward_run_b0 et ca sera exactement le base_forward_run mais sans le calcul de la fonction cout.

@maximejay maximejay Feb 3, 2026

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oui, bonne idée !
J'ai fait comme tu as proposé. Ca fonctionne bien.

- New diff rule for q (qdomain) and qe (elemntary discharge)
- rename response.qt -> response.qac
- New forward routine to handle new diff rules
@maximejay

Copy link
Copy Markdown
Collaborator Author

Salut François,
Voici une nouvelle proposition, avec une option dans le setup pour sortir les gradients par rapport à qt ou q sur toute la grille.
Après quelques réflexion supplémentaires: Les sorties sont copiées dans la variable response%qac. En fonction du mode choisi, c'est un peu redondant avec la sortie qdomain (output option) et qt (internal flux). Mais on a pas vraiment le choix sur ce point. A se demander si on garde qdomain du coup ?

Toujours un problème de make check/format... j'ai pourtant fais la mise )à jour des paquets.

ou souhaites tu que l'on documente cette option pour renvoyer les nouveaux gradients dans la doc ?

A++

@inoelloc
inoelloc self-requested a review February 4, 2026 10:01
Comment thread smash/fcore/forward/md_simulation.f90
Comment thread smash/fcore/forward/md_simulation.f90 Outdated

end do

if (setup%return_opt_grad == "qe") then

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Pourquoi nommer qe et pas qt ?

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Pour Q élémentaire (1 pixel). On peut changer de nom si il faut.
Peut-être on peu trouver un meilleur terme d'ailleurs ?

@inoelloc

Copy link
Copy Markdown
Collaborator

Je suis pas sur de bien comprendre comment on récupère les gradients au final. Quelle fonction de l'API renvoie les gradients avec cette option ?
Je vais m'occuper du formatage du code, des conflits et des détails restants

@maximejay

Copy link
Copy Markdown
Collaborator Author

Salut,

Les gradients renvoyés sont aux choix:

  • qe/parametters -> débits élémentaires (sortie du module hydrologique) par rapport aux paramètres
  • q/parametters -> débits dans le réseaux hydro par rapport aux paramètres

L'option qui permet de contrôler cela s'appelle : setup%return_opt_grad

L'API Fortran renvoie les gradients via les fonctions suivantes:

  • forward_run_q_d dans mw_forward.f90 -> TLM
  • forward_run_q_b dans mw_forward.f90 -> Adjoint

L'API Python est la fonction :

  • _get_parameters_q_b dans core/simulation/optimize/_tools.py
  • on appelle le fortran via import: from smash.fcore._mw_forward import forward_run_q_b as wrap_forward_run_q_b

Pour avoir les gradient scoté Python j'utilise la fonction

parameters_q_b = smash.core.simulation.optimize._tools._get_parameters_q_b(
        model, parameters, wrap_options, wrap_returns
    )

C'est là ou un utilisateur de Smash ne peut rien faire car l'appel de cette dernière est complexe. Je me suis inspiré des fonctions dans core/simulation/optimize/optimize.py pour en déduire la fonction suivante utilisé dans Gamma. Je pense que ce serait bien de l'intégrer dans l'API Python pour que les utilisateur puisse utiliser ce gradient. De même pour le gradient de base Cost/parameters ? Voici la fonction python:

def _get_smash_gradient(model, control_smash):

    from smash.fcore._mwd_options import OptionsDT
    from smash.fcore._mwd_returns import ReturnsDT
    from smash.core.model._build_model import _map_dict_to_fortran_derived_type
    from smash.core.simulation.optimize._standardize import _standardize_optimize_args
    from smash.fcore._mwd_parameters_manipulation import (
        parameters_to_control as wrap_parameters_to_control,
    )
    from smash.core.simulation.optimize._tools import (
        _handle_bayesian_optimize_control_prior,
    )

    (
        mapping,
        optimizer,
        optimize_options,
        cost_options,
        common_options,
        return_options,
        callback,
    ) = _standardize_optimize_args(
        model,
        mapping="distributed",
        optimizer="lbfgsb",
        optimize_options={
            "parameters": control_smash["ParamList"],
            "bounds": control_smash["bounds"],
        },
        cost_options=None,
        common_options=None,
        return_options=None,
        callback=None,
    )

    wrap_options = OptionsDT(
        model.setup,
        model.mesh,
        cost_options["njoc"],
        cost_options["njrc"],
    )

    wrap_returns = ReturnsDT(
        model.setup,
        model.mesh,
        return_options["nmts"],
        return_options["fkeys"],
    )

    # % Map optimize_options dict to derived type
    _map_dict_to_fortran_derived_type(optimize_options, wrap_options.optimize)

    # % Map cost_options dict to derived type
    _map_dict_to_fortran_derived_type(
        cost_options, wrap_options.cost, skip=["control_prior"]
    )

    # % Map common_options dict to derived type
    _map_dict_to_fortran_derived_type(common_options, wrap_options.comm)

    # % Map return_options dict to derived type
    _map_dict_to_fortran_derived_type(return_options, wrap_returns)

    parameters = model._parameters.copy()

    wrap_parameters_to_control(
        model.setup,
        model.mesh,
        model._input_data,
        parameters,
        wrap_options,
    )

    parameters_q_b = smash.core.simulation.optimize._tools._get_parameters_q_b(
        model, parameters, wrap_options, wrap_returns
    )

    grad = parameters_q_b.control.x.copy()

    return grad

@inoelloc

Copy link
Copy Markdown
Collaborator

Ok ok je vois. Ca me semble bizarre d'ajouter un argument au setup qu'un utilisateur ne peut pas utiliser. Il nous faut la fonction API publique pour recuperer les gradients. Sinon, on reste au niveau du dev. On pourrait avoir quelque chose comme ca :

model.backward_run(mapping='uniform')
cost = model.output.cost
grad = model.output.grad

On garde tjs les tableaux dans model.output.response pour q ou qt. Je garderais qt qui est le nom de la variable dans le code, dans la documentation et dans les flux internes et faudrait allouer en sortie la variable output.grad à la bonne dimension.
Et pour moi, vu comment le code a été pensé, je mettrais dans cette fonction l'option de récupérer les gradients par rapport à la fonction cout ou par rapport à une autre variable d’intérêt. Peut etre :

# Grad Cost (par defaut)
model.backward_run(mapping='uniform')

# Grad Q
model.backward_run(mapping='uniform', grad_kind="q")

Si grad_kind est pas clair, on peut changer, peu importe ca

@nghi-truyen nghi-truyen added this to the Release v1.3.0 milestone Mar 18, 2026
@maximejay

maximejay commented Mar 31, 2026

Copy link
Copy Markdown
Collaborator Author

Salut,

J'ai créé une nouvelle fonction dans la class model, qui s'appelle backward_run. Cette fonction prend 2 arguments:

  • le mapping

  • grad_mode: une valeur parmis ["j","q","qt"]

  • si grad_mode == j, on renvoie le gradient dJ/dX

  • si grad_mode == q, on renvoie le gradient dQ/dX

  • si grad_mode == qt, on renvoie le gradient dQt/dX

Lorsque que l'on appelle backward_run(), l'objet returns_options est complété en fonction de la valeur de grad_mode.
Coté Fortran, l'objet returnsDT a deux nouvelles propriétés:

  • q_domain_kind : str, ["j","q","qt"]
  • q_domain_kind_flag : le flag associé qui est seulement utilisé pour tester si on demande un gradient q ou qt

J'ai nettoyé le code (setup et constant.py) des anciennes variables (cf commit plus haut).

j'ai ajouté de la doc dans core/simulation/_doc.py.

Il faudrait faire un script pour tester et valider ces gradients...

Voici un script pour tester:

import smash

mapping = "uniform"

setup_cance, mesh_cance = smash.factory.load_dataset("Cance")

Avec routage

mapping = "distributed"

setup_cance, mesh_cance = smash.factory.load_dataset("Cance")

smash_model = smash.Model(setup_cance, mesh_cance)
smash_model.forward_run()
grad1, name1 = smash_model.backward_run(mapping=mapping, optimizer="lbfgsb", grad_mode="j")

grad2, name2 = smash_model.backward_run(mapping=mapping, optimizer="lbfgsb", grad_mode="q")
grad3, name3 = smash_model.backward_run(mapping=mapping, optimizer="lbfgsb", grad_mode="qt")

La normalisation des param dépend du type d'optimiseur, cela change les gradients

grad3, name3 = smash_model.backward_run(mapping="uniform", optimizer="sbs", grad_mode="qt")
grad3, name3 = smash_model.backward_run(mapping="uniform", optimizer="lbfgsb", grad_mode="qt")

sans routage

setup_cance["routing_module"] = "zero"
smash_model_zero_routing = smash.Model(setup_cance, mesh_cance)
smash_model_zero_routing.forward_run()
grad4, name4 = smash_model_zero_routing.backward_run(mapping=mapping, grad_mode="qt")
grad5, name5 = smash_model_zero_routing.backward_run(mapping=mapping, grad_mode="q")

# grad4 et grad5 sont identiques (pas de routage, q=qt)
# grad3 et grad4 sont identiques, seulement grad3 est un tableau plus grand, contenant les les gradient pour llr (0.0)

…peu pas récupérer les débits qt sur la grille ...

Pour l'instant c'est fait dans la class model diretement dans forward_run...

Suppression de option_b/_d dans base_forward_run_d et _b
@maximejay

Copy link
Copy Markdown
Collaborator Author

Un dernier commit, car je me suis aperçu qu'il fallait aussi allouer response%qac dans forward_run afin de pouvoir récupérer les débits qt sur la grille. Pour l'instant j'ai mis ca direct dans la class model -> forward_run. Il faudra peut-être le déplacer dans les standardize ?

Merci
A++

…ith previous version of smash. Get qt in the same way than q, i.e with q_domain flag.

Add flag to get the gradient for q and qt vs X, and allocate output%response%qac
- New input_derivatives parameters to pass the gradients of an external model to the adjoint of Smash: the full gradient of the model chain in therefore obtained : dj/dX=dj/dq*dq/dx. This fix issue for calibrating smash when coupled with the gamma routing model
@maximejay

maximejay commented May 6, 2026

Copy link
Copy Markdown
Collaborator Author

Finalisation du travail:

  • Nouvelle variable input_derivatives permettant de faire passer le résultat de dérivés adjointe d'un autre model vers l'adjoint de smash. Cela permet de coupler les 2 gradients ensemble et d'accomplir le calcul de dJ/dX=dJ/dq*dq/dX.
  • Ajout d'une condition dans io.model.read_model() afin de réallouer le tableau qac si besoin + test son existence pour un compatibilité descendante

@inoelloc inoelloc removed this from the Release v1.3.0 milestone Sep 4, 2026
@inoelloc

inoelloc commented Sep 4, 2026

Copy link
Copy Markdown
Collaborator

Done by PR #541

@inoelloc inoelloc closed this Sep 4, 2026
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

enhancement New feature or request

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants