Pamhyr2/src/Solver/RubarBE.py

# RubarBE.py -- Pamhyr
# Copyright (C) 2024-2025  INRAE
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <https://www.gnu.org/licenses/>.

# -*- coding: utf-8 -*-

import os
import logging
import math
from decimal import Decimal, InvalidOperation, localcontext
import numpy as np

from tools import timer, trace, old_pamhyr_date_to_timestamp, logger_exception

from Solver.CommandLine import CommandLineSolver
from View.RunSolver.Translate import RunSolverTranslate

from Model.Results.Results import Results
from Model.Results.River.River import River, Reach, Profile

logger = logging.getLogger()


class Rubar3(CommandLineSolver):
    _type = "rubar3"

    def __init__(self, name):
        super(Rubar3, self).__init__(name)

        self._type = "rubar3"

        self._cmd_input = ""
        self._cmd_solver = "@path @args @input"
        self._cmd_output = ""

    @classmethod
    def default_parameters(cls):
        # Keep historical rubarbe_* keys for existing Rubar3 studies.

        lst = [
            ("rubarbe_cfl", "0.50000E+00"),
            ("rubarbe_condam", "1"),
            ("rubarbe_condav", "3"),
            ("rubarbe_regime", "0"),
            ("rubarbe_iodev", "n"),
            ("rubarbe_iodebord", "o"),
            ("rubarbe_iostockage", "o"),
            ("rubarbe_iopdt", "y"),
            ("rubarbe_iovis", "n"),
            ("rubarbe_rep", "n"),
            ("rubarbe_tinit", "000:00:00:00"),
            ("rubarbe_tmax", "001:00:00:00"),  # 1 day
            ("rubarbe_tiopdt", "000:00:00:00"),
            ("rubarbe_dt", "5.0"),
            ("rubarbe_ts", "999:99:99:00"),
            ("rubarbe_dtsauv", "00:00:05:00"),
            ("rubarbe_psave", "00:00:05:00"),
            ("rubarbe_fdeb1", "1"),
            ("rubarbe_fdeb2", "10"),
            ("rubarbe_fdeb3", "100"),
            ("rubarbe_tf_1", "y"),
            ("rubarbe_tf_2", "y"),
            ("rubarbe_tf_3", "y"),
            ("rubarbe_tf_4", "y"),
            ("rubarbe_tf_5", "y"),
            ("rubarbe_tf_6", "n"),
            ("rubarbe_trased", "n"),
            ("rubarbe_optfpc", "0"),
            # trased parameters
            ("rubarbe_ros", "2650.0"),
            ("rubarbe_dm", "0.1"),
            ("rubarbe_segma", "1.0"),
            # Sediment parameters
            ("rubarbe_sediment_ros", "2650.0"),
            ("rubarbe_sediment_por", "0.4"),
            ("rubarbe_sediment_dcharg", "0.0"),
            ("rubarbe_sediment_halfa", "1.0"),
            ("rubarbe_sediment_mult_1", "1.0"),
            ("rubarbe_sediment_mult_2", ""),
            ("rubarbe_sediment_mult_3", ""),
            ("rubarbe_sediment_mult_4", ""),
            ("rubarbe_sediment_mult_5", ""),
            ("rubarbe_sediment_visc", "1.3E-06"),
            ("rubarbe_sediment_opts", "3"),
            ("rubarbe_sediment_odchar", "0"),
            ("rubarbe_sediment_unisol", "1"),
            ("rubarbe_sediment_typdef", "3"),
            ("rubarbe_sediment_depot", "2"),
            ("rubarbe_sediment_choixc", "2"),
            ("rubarbe_sediment_option", "2"),
            ("rubarbe_sediment_capsol", "1"),
            ("rubarbe_sediment_bmiu", "0.85"),
            ("rubarbe_sediment_demix", "0"),
            ("rubarbe_sediment_defond", "1"),
            ("rubarbe_sediment_varcons", "1"),
            ("rubarbe_sediment_dchard", "0.0"),
            ("rubarbe_sediment_dchars", "0.0"),
        ]

        return lst

    @classmethod
    def checkers(cls):
        lst = [
        ]

        return lst

    ##########
    # Export #
    ##########

    def cmd_args(self, study):
        lst = super(Rubar3, self).cmd_args(study)

        return lst

    def input_param(self):
        name = self._study.name
        return f"{name}"

    def output_param(self):
        name = self._study.name
        return f"profil.{name}"

    def log_file(self):
        name = self._study.name
        return f"geomac-i.{name}"

    def export(self, study, repertory, qlog=None):
        self._study = study
        name = study.name.replace(" ", "_")

        self._export_donnee(study, repertory, qlog, name=name)
        if self._parameter_flag(study, "trased") == "O":
            if not self._export_ts(study, repertory, qlog, name=name):
                return False
        self._export_geomac_i(study, repertory, qlog, name=name)
        self._export_mail(study, repertory, qlog, name=name)
        self._export_devers(study, repertory, qlog, name=name)
        self._export_condin(study, repertory, qlog, name=name)
        self._export_stricklers(study, repertory, qlog, name=name)
        self._export_hydro(study, repertory, qlog, name=name)
        self._export_condav(study, repertory, qlog, name=name)
        # self._export_abshyd(study, repertory, qlog, name=name)

        return True

    # Notice Rubar3/RubarBE, annex 1: records are positional, not keyed.
    # Each record starts with an A50 label, including records of integers.
    _donnee_fields = (
        ("cfl", "real"),
        ("condam", "integer"), ("condav", "integer"),
        ("regime", "integer"), ("iodev", "flags"),
        ("iopdt", "flag"), ("iovis", "flag"), ("rep", "flag"),
        ("tinit", "time"), ("tmax", "time"), ("tiopdt", "time"),
        ("dt", "real"), ("ts", "time"), ("dtsauv", "time"),
        ("psave", "time"), ("fdeb1", "real"), ("fdeb2", "real"),
        ("fdeb3", "real"),
        ("tf_1", "flag"), ("tf_2", "flag"), ("tf_3", "flag"),
        ("tf_4", "flag"), ("tf_5", "flag"), ("tf_6", "flag"),
        ("trased", "flag"),
    )
    _ts_fields = (
        ("ros", "real"), ("por", "real"), ("dcharg", "real"),
        ("halfa", "real"), ("mult_1", "multipliers"),
        ("visc", "real"), ("opts", "law"), ("odchar", "integer"),
        ("unisol", "integer"), ("typdef", "integer"),
        ("depot", "depot"), ("choixc", "integer"),
        ("option", "options"), ("capsol", "integer"),
        ("bmiu", "real"), ("demix", "integer"),
        ("defond", "integer"), ("varcons", "integer"),
    )

    def _parameter_values(self, study):
        return {p.name: str(p.value).strip()
                for p in study.river.get_params(self.type).parameters}

    @staticmethod
    def _required(values, key):
        value = values.get(key, "")
        if not value:
            raise ValueError(f"Missing Rubar parameter: {key}")
        return value

    @staticmethod
    def _flag(value, key, mixed=False):
        result = value.upper()
        if result == "Y":
            result = "O"
        if result not in (("O", "N", "M") if mixed else ("O", "N")):
            raise ValueError(f"{key}: expected O/Y or N" +
                             (" or M" if mixed else ""))
        return result

    def _parameter_flag(self, study, suffix):
        key = "rubarbe_" + suffix
        return self._flag(self._required(self._parameter_values(study), key),
                          key)

    @staticmethod
    def _integer(value, key, width=1):
        if not value.isascii() or not value.isdigit() or len(value) > width:
            raise ValueError(f"{key}: expected an integer in {width} columns")
        return value.rjust(width)

    @staticmethod
    def _real(value, key, scientific=False):
        try:
            number = Decimal(value.replace("d", "E").replace("D", "E"))
        except InvalidOperation:
            raise ValueError(f"{key}: invalid real value {value!r}") from None
        if not number.is_finite():
            raise ValueError(f"{key}: expected a finite real value")
        if scientific:
            # DONNEE: E12.5. Reject overflow instead of shifting columns.
            floating = float(number)
            if not math.isfinite(floating) or (number and floating == 0):
                raise ValueError(f"{key}: real value outside E12.5 range")
            result = f"{floating:12.5E}"
            width = 12
        else:
            # TS: F10.0 input accepts an explicit decimal point/exponent.
            # Preserve the value exactly; do not round porosity etc. to int.
            with localcontext() as context:
                context.prec = max(28, len(number.as_tuple().digits))
                number = number.normalize()
            result = str(number).upper()
            if len(result) > 10:
                result = format(number, "E")
            width = 10
        if len(result) > width:
            raise ValueError(f"{key}: value {value!r} exceeds {width} columns")
        return result.rjust(width)

    @staticmethod
    def _write_parameters(repertory, filename, lines):
        # Validate all records before opening/truncating the destination.
        path = os.path.join(repertory, filename)
        with open(path, "w", encoding="ascii") as f:
            f.writelines(lines)

    def _export_donnee(self, study, repertory, qlog=None, name="0"):
        values = self._parameter_values(study)
        fields = list(self._donnee_fields)
        # The diagnostic transport-capacity extension is optional as a block.
        if values.get("rubarbe_optfpc", ""):
            fields += [("optfpc", "integer"), ("ros", "real"),
                       ("dm", "real"), ("segma", "real")]
        elif any(values.get("rubarbe_" + key, "")
                 for key in ("ros", "dm", "segma")):
            raise ValueError("rubarbe_optfpc: required before ROS/DM/SEGMA")
        lines = []
        for suffix, kind in fields:
            key = "rubarbe_" + suffix
            value = self._required(values, key)
            if kind == "integer":
                result = self._integer(value, key)
            elif kind in ("flag", "flags"):
                result = self._flag(value, key, mixed=suffix == "iopdt")
                if kind == "flags":
                    # Missing optional flags mean O (notice section 3.1).
                    # Write both explicitly so STOCKAGE cannot shift left.
                    for extra in ("iodebord", "iostockage"):
                        extra_key = "rubarbe_" + extra
                        result += self._flag(values.get(extra_key) or "O",
                                             extra_key)
            else:
                if kind == "time" and ":" in value:
                    parts = value.split(":")
                    if len(parts) != 4 or not all(p.isdigit() for p in parts):
                        raise ValueError(
                            f"{key}: expected days:hours:minutes:seconds")
                    value = str(old_pamhyr_date_to_timestamp(value))
                result = self._real(value, key, scientific=True)
            lines.append(f"{key:<50}{result}\n")
        if qlog is not None:
            qlog.put("Export DONNEE file")
        self._write_parameters(repertory, f"donnee.{name}", lines)

    def _export_ts(self, study, repertory, qlog=None, name="0"):
        trad = RunSolverTranslate()
        values = self._parameter_values(study)
        prefix = "rubarbe_sediment_"
        fields = list(self._ts_fields)
        demix = self._required(values, prefix + "demix")
        if demix not in ("0", "1"):
            raise ValueError(f"{prefix}demix: expected 0 or 1")
        if demix == "1":
            fields += [("dchard", "real"), ("dchars", "real")]
        lines = []
        for suffix, kind in fields:
            key = prefix + suffix
            value = self._required(values, key)
            if kind == "real":
                result = self._real(value, key)
            elif kind == "multipliers":
                result = self._real(value, key)
                extras = [values.get(prefix + f"mult_{i}", "")
                          for i in range(2, 6)]
                # Preserve empty interior F10.0 fields; omit trailing fields
                # to retain the solver's defaults for optional coefficients.
                while extras and not extras[-1]:
                    extras.pop()
                for i, extra in enumerate(extras, 2):
                    result += (self._real(extra, prefix + f"mult_{i}")
                               if extra else " " * 10)
            elif kind == "options":
                if len(value) != 3 or any(c not in "12" for c in value):
                    message = trad["invalid_options"].format(key=key)
                    logger.warning(message)
                    if qlog is not None:
                        qlog.put(trad["warning"].format(message=message))
                    return False
                result = value
            elif kind == "law":
                # Notice section 3.10 / 8.7: 4--7 are not listed.
                if value not in {str(i) for i in [1, 2, 3, *range(8, 31)]}:
                    message = trad["unsupported_transport_law"].format(
                        key=key, value=value)
                    logger.warning(message)
                    if qlog is not None:
                        qlog.put(trad["warning"].format(message=message))
                    return False
                result = self._integer(value, key, width=2)
            else:
                result = self._integer(value, key,
                                       width=10 if kind == "depot" else 1)
            lines.append(f"{key:<50}{result}\n")
        if qlog is not None:
            qlog.put("Export TS file")
        self._write_parameters(repertory, f"ts.{name}", lines)
        return True

    def _export_geomac_i(self, study, repertory, qlog, name="0"):
        if qlog is not None:
            qlog.put("Export GEOMAC-i file")

        with open(
                os.path.join(
                    repertory, f"geomac-i.{name}"
                ), "w+"
        ) as f:
            for edge in study.river.enable_edges():
                reach = edge.reach
                n_profiles = len(reach)
                time = 0.0

                f.write(f"{n_profiles:>5} {time:>11.3f}\n")

                ind = 1
                for profile in reach.profiles:
                    rk = profile.rk
                    n_points = len(profile)

                    f.write(f"{ind:>4} {rk:>11.3f} {n_points:>4}\n")

                    station = profile.get_station()
                    for i, point in enumerate(profile.points):
                        label = point.name.lower()
                        if label != "":
                            if label[0] == "r":
                                label = label[1].upper()
                            else:
                                label = " "
                        else:
                            label = " "

                        y = station[i]
                        z = point.z
                        dcs = 1.0
                        scs = 1.0
                        tmcs = 0.0

                        f.write(
                            f"{label} {y:>11.5f}" +
                            f"{z:>13.5f}{dcs:>15.10f}" +
                            f"{scs:>15.10f}{tmcs:>15.5f}" +
                            "\n"
                        )

                    ind += 1

    def _export_mail(self, study, repertory, qlog, name="0"):
        if qlog is not None:
            qlog.put("Export MAIL file")

        with open(
                os.path.join(
                    repertory, f"mail.{name}"
                ), "w+"
        ) as f:
            for edge in study.river.enable_edges():
                reach = edge.reach
                lm = len(reach) + 1
                f.write(f"{lm:>13}\n")

                for mails in [reach.inter_profiles_rk(),
                              reach.get_rk()]:
                    ind = 0
                    for mail in mails:
                        f.write(f"{mail:15.3f}")

                        ind += 1
                        if ind % 3 == 0:
                            f.write("\n")

                    if ind % 3 != 0:
                        f.write("\n")

    def _export_abshyd(self, study, repertory, qlog, name="0"):
        if qlog is not None:
            qlog.put("Export ABSHYD file")

        with open(
                os.path.join(
                    repertory, f"abshyd.{name}"
                ), "w+"
        ) as f:
            reach_ind = 1
            for edge in study.river.enable_edges():
                reach = edge.reach
                lm = len(reach) + 1
                f.write(f"{lm:>13}\n")

                ind = 1
                for mail in reach.inter_profiles_rk():
                    f.write(f"{ind:>4} {mail:15.3f} {reach_ind:>4}\n")
                    ind += 1

                reach_ind += 1

    def _export_devers(self, study, repertory, qlog, name="0"):
        if qlog is not None:
            qlog.put("Export DEVERS file")

        with open(
                os.path.join(
                    repertory, f"devers.{name}"
                ), "w+"
        ) as f:
            reach_ind = 1
            for edge in study.river.enable_edges():
                reach = edge.reach
                lm = len(reach) + 1
                f.write(f"{lm:>13}\n")

                ind = 1
                for mail in reach.inter_profiles_rk():
                    f.write(f"{ind:>4} {0.0:15.3f} {0.0:15.3f} {1.0:>15.3f}\n")
                    ind += 1

                reach_ind += 1

    def _export_stricklers(self, study, repertory, qlog, name="0"):
        self._export_frot(study, repertory, qlog, name=name, version="")
        self._export_frot(study, repertory, qlog, name=name, version="2")

    def _export_frot(self, study, repertory, qlog, name="0", version=""):
        if qlog is not None:
            qlog.put(f"Export FROT{version} file")

        with open(
                os.path.join(
                    repertory, f"frot{version}.{name}"
                ), "w+"
        ) as f:
            for edge in study.river.enable_edges():
                reach = edge.reach
                lm = len(reach) + 1
                f.write(f"{lm:>6}\n")

                lst = list(filter(
                    lambda f: f.is_full_defined(),
                    edge.frictions.frictions
                ))

                rk_min = 9999999.9
                rk_max = -9999999.9
                coeff_min = -1.0
                coeff_max = -1.0
                for s in lst:  # TODO optimise ?
                    if s.begin_rk > rk_max:
                        rk_max = s.begin_rk
                        coeff_max = s.begin_strickler
                    if s.begin_rk < rk_min:
                        rk_min = s.begin_rk
                        coeff_min = s.begin_strickler
                    if s.end_rk > rk_max:
                        rk_max = s.end_rk
                        coeff_max = s.end_strickler
                    if s.end_rk < rk_min:
                        rk_min = s.end_rk
                        coeff_min = s.end_strickler

                def get_stricklers_from_rk(rk, lst):

                    coeff = None
                    if rk > rk_max:
                        coeff = coeff_max
                    elif rk < rk_min:
                        coeff = coeff_min
                    else:
                        for s in lst:
                            if (rk >= s.begin_rk and rk <= s.end_rk or
                                    rk <= s.begin_rk and rk >= s.end_rk):
                                coeff = s.begin_strickler  # TODO: inerpolate
                                break

                    # TODO interpolation if rk is not in frictons

                    if coeff is None:
                        logger.error(
                            "Study frictions are not fully defined"
                        )
                        return None

                    return coeff.medium if version == "2" else coeff.minor

                ind = 1
                for mail in edge.reach.inter_profiles_rk():
                    coef = get_stricklers_from_rk(mail, lst)
                    if coef is not None:
                        f.write(f"{ind:>6} {coef:>12.5f}")

                    ind += 1
                    f.write("\n")

    def _export_condin(self, study, repertory, qlog, name="0"):
        if qlog is not None:
            qlog.put("Export CONDIN file")

        with open(
                os.path.join(
                    repertory, f"condin.{name}"
                ), "w+"
        ) as f:
            for edge in study.river.enable_edges():
                reach = edge.reach

                f.write(f"{0.0}\n")

                ics = study.river.initial_conditions.get(edge)
                data = self._export_condin_init_data(ics)

                profiles = reach.profiles
                first = profiles[0]
                last = profiles[-1]

                if first not in data or last not in data:
                    logger.error(
                        "Study initial condition is not fully defined"
                    )
                    return

                f_h_s = self._export_condin_profile_height_speed(first, data)
                l_h_s = self._export_condin_profile_height_speed(last, data)

                # First mail
                f.write(f"{1:>5} {f_h_s[0]} {f_h_s[1]}\n")

                ind = 2
                it = iter(profiles)
                prev = next(it)

                prev_h, prev_s = f_h_s
                for profile in it:
                    if profile not in data:
                        ind += 1
                        continue

                    cur_h, cur_s = self._export_condin_profile_height_speed(
                        profile, data
                    )

                    # Mean of height and speed
                    h = (prev_h + cur_h) / 2
                    s = (prev_s + cur_s) / 2

                    f.write(f"{ind:>5} {h} {s}\n")

                    prev_h, prev_s = cur_h, cur_s
                    ind += 1

                # Last mail
                f.write(f"{ind:>5} {f_h_s[0]} {f_h_s[1]}\n")

    def _export_condin_init_data(self, ics):
        data = {}

        for d in ics.data:
            data[d['rk']] = (
                d['elevation'],
                d['discharge'],
            )

        return data

    def _export_condin_profile_height_speed(self, profile, data):
        z = data[profile][0]
        q = data[profile][1]

        height = z - profile.z_min()
        speed = profile.speed(q, z)

        return height, speed

    def _export_hydro(self, study, repertory, qlog, name="0"):
        if qlog is not None:
            qlog.put("Export HYDRO file")

        with open(
                os.path.join(
                    repertory, f"hydro.{name}"
                ), "w+"
        ) as f:
            bcs = []
            for edge in study.river.enable_edges():
                for bound in study.river.boundary_condition.get_tab("liquid"):
                    # BC is an hydrogramme
                    if bound.bctype == "TD" or bound.bctype == "PC":
                        # BC is on input node of this reach
                        if bound.node == edge.node1:
                            bcs.append(bound)

            for bc in bcs:
                f.write(f"{len(bc)}\n")
                for d0, d1 in bc.data:
                    f.write(f"{d1} {d0}\n")

    def _export_condav(self, study, repertory, qlog, name="0"):
        if qlog is not None:
            qlog.put("Export CONDAV file")

        with open(
                os.path.join(
                    repertory, f"condav.{name}"
                ), "w+"
        ) as f:
            bcs = []
            for edge in study.river.enable_edges():
                for bound in study.river.boundary_condition.get_tab("liquid"):
                    # BC is an hydrogramme
                    if bound.bctype == "ZD" or bound.bctype == "TZ":
                        # BC is on input node of this reach
                        if bound.node == edge.node2:
                            bcs.append(bound)

            for bc in bcs:
                f.write(f"{len(bc)}\n")
                if bc.bctype == "ZD":
                    for d0, d1 in bc.data:
                        f.write(f"{d1} {d0}\n")
                else:
                    for d0, d1 in bc.data:
                        f.write(f"{d0} {d1}\n")

    def read_profil(self, study, fname, results, qlog=None, name="0"):
        logger.info(f"read: profil.{name}")
        with open(fname, "r") as f:
            reachs = []
            for i, edge in enumerate(study.river.enable_edges()):
                reach = edge.reach
                # Add results reach to reach list
                res_reach = results.river.add(i)
                reachs.append((res_reach, len(reach) + 1))
            results.set("nb_reach", f"{len(reachs)}")
            results.set("nb_profile", f"{sum([r[1] - 1 for r in reachs])}")

            def read_data_line(f):
                line = f.readline().split()
                if len(line) > 3:
                    return tuple(map(float, line))
                else:
                    return (0.0, 0.0, 0.0, 0.0)

            def set_and_compute_limites(reach, ind, z, q, s):
                reach.set(ind, timestamp, "Z", z)
                reach.set(ind, timestamp, "Q", q)
                reach.set(ind, timestamp, "V", s)

                profile = reach.profile(ind)
                limits = profile.geometry.get_water_limits(z)
                reach.set(
                    ind, timestamp, "water_limits", limits
                )

            ts = set()
            timestamp = next(filter(
                lambda p: p.name == 'rubarbe_tinit',
                study.river.get_params(self._type).parameters
            )).value
            if timestamp.count(':') == 3:
                timestamp = old_pamhyr_date_to_timestamp(timestamp)
            timestamp = float(timestamp)
            ts.add(timestamp)

            # add initial condition
            for r, edge in enumerate(study.river.enable_edges()):
                reach = edge.reach
                ics = study.river.initial_conditions.get(edge)
                q = ics.get_discharge()
                z = ics.get_elevation()
                k = 0
                for i, j in zip(z, q):
                    v = reach.profiles[k].speed(i, j)
                    set_and_compute_limites(reachs[r][0], k, i, j, v)
                    k += 1

            # start read

            end = False
            while True:
                line = f.readline()
                if line == "":
                    break

                timestamp = float(line)
                ts.add(timestamp)

                logger.info(f"read: profil.{name}: timestamp = {timestamp}")

                for reach, lm in reachs:
                    # First profile
                    h, s, q, z = read_data_line(f)
                    set_and_compute_limites(reach, 0, z+h, q, s)

                    # For each profile
                    ind = 1
                    ph, ps, pq, pz = read_data_line(f)
                    while ind < lm - 2:
                        h, s, q, z = read_data_line(f)
                        set_and_compute_limites(
                            reach, ind,
                            (pz + z + ph + h) / 2,
                            (pq + q) / 2,
                            (ps + s) / 2
                        )

                        ph = h
                        ps = s
                        pq = q
                        pz = z
                        ind += 1

                    # Last profile
                    h, s, q, z = read_data_line(f)
                    set_and_compute_limites(reach, ind, z+h, q, s)

            # Tables and plots use timestamp rows and global profile columns.
            timestamps = sorted(ts)
            profiles = []
            for reach, _ in reachs:
                offset = len(profiles)
                profiles.extend(reach.profiles)
                reach.set_global_index(range(offset, len(profiles)))

            results.set("timestamps", timestamps)
            results.compute_timestamp_index()
            results.set("table", {
                key: np.array([
                    [profile.get_ts_key(t, key) for profile in profiles]
                    for t in timestamps
                ], dtype=float)
                for key in ("Z", "Q", "V")
            })

    @timer
    def write_bin(self, study, fname, results, qlog=None, name="0"):
        logger.info(f"write_bin: Start writing '{fname}' ...")

        with open(fname, "wb") as f:
            def newline(j): return np.asarray([j], dtype=np.int32).tofile(f)
            def endline(j): return np.asarray([j], dtype=np.int32).tofile(f)

            def write_int(i):
                newline(len(i)*4)
                np.array(i, dtype=np.int32).tofile(f)
                endline(len(i)*4)

            def write_float(i):
                newline(len(i)*4)
                np.array(i, dtype=np.float32).tofile(f)
                endline(len(i)*4)

            def write_float64(i):
                newline(len(i)*8)
                np.array(i, dtype=np.float64).tofile(f)
                endline(len(i)*8)

            def write_float64(i):
                newline(len(i)*8)
                np.array(i, dtype=np.float64).tofile(f)
                endline(len(i)*8)

            def write_data(npts, t, a, val):
                newline(npts*4 + 13)
                np.array(npts, dtype=np.int32).tofile(f)
                np.array(t, dtype=np.float64).tofile(f)
                np.array(bytearray(a.encode()), dtype=np.byte).tofile(f)
                np.array(val, dtype=np.float32).tofile(f)
                endline(npts*4 + 13)

            ts_list = sorted(results.get("timestamps"))
            profiles = []
            for r in results.river.reachs:
                profiles += r.profiles

            # Meta data (1st line)
            nb_reach = len(results.river)
            nb_profile = len(profiles)
            write_int([nb_reach, nb_profile, "82"])

            # Reach information (2nd line)
            is1 = []
            is2 = []
            ltmp = []
            i = 1
            for r in results.river.reachs:
                j = i+len(r)-1
                ltmp += [i, j]
                is1.append(i)
                is2.append(j)
                i += i+len(r)
            write_int(ltmp)

            # X (3rd line)
            rk = []
            for r in results.river.reachs:
                rk += r.geometry.get_rk()
            write_float(rk)

            # Z and Y (4th line)
            # ltmp = []
            # for r in results.river.reachs:
            #     for p in r.prifiles:
            #         ltpm.append(p.rk)
            write_float(3*nb_profile*[0.0])

            # Data

            for timestamp in ts_list:
                q = list(
                    map(
                        lambda p: p.get_ts_key(timestamp, "Q"),
                        profiles
                    )
                )
                write_data(nb_profile, timestamp, "Q", q)
                z = list(
                    map(
                        lambda p: p.get_ts_key(timestamp, "Z"),
                        profiles
                    )
                )
                write_data(nb_profile, timestamp, "Z", z)
            logger.info(f"write_bin: ... end with {len(ts_list)} timestamps")

    @timer
    def results(self, study, repertory, qlog=None, name=None):
        results = Results(
            study=study,
            solver=self,
            repertory=repertory,
            name=name,
        )
        if name is None:
            name = study.name.replace(" ", "_")

        results_file = f"profil.{name}"

        fname = os.path.join(repertory, results_file)
        if not os.path.isfile(fname):
            logger.warning(f"Result file {results_file} does not exist")
            return None
        try:
            self.read_profil(study, fname, results, qlog, name=name)
        except Exception as e:
            logger.error(f"Failed to read results")
            logger_exception(e)
            return None

        return results


class RubarBE(Rubar3):
    _type = "rubarbe"

    def __init__(self, name):
        super(RubarBE, self).__init__(name)

        self._type = "rubarbe"

        self._cmd_input = ""
        self._cmd_solver = "@path @args @input"
        self._cmd_output = ""

    @classmethod
    def default_parameters(cls):
        # New BE studies evolve the bed. Keep Rubar3 defaults unchanged.
        overrides = {
            "rubarbe_trased": "y",
            "rubarbe_sediment_opts": "1",  # Meyer-Peter and Mueller
            "rubarbe_sediment_option": "212",  # no Ikeda, given Shields, total
        }
        return [(key, overrides.get(key, value))
                for key, value in super().default_parameters()]