class BlueKenue:
    def __init__(self, filename, mode):
        """!
        @brief BlueKenue file object
        @param filename <str>: path to BlueKenue file
        @param mode <str>: `r` for read, `w` or `x` for write mode
        """
        self.filename = filename
        self.mode = mode
        self.file = None
        self.header = None

    def __enter__(self):
        self.file = open(self.filename, self.mode, encoding='ISO-8859-1')
        return self

    def __exit__(self, exc_type, exc_val, exc_tb):
        self.file.close() def get_shape_type(input_filename):
def get_lines(input_filename, shape_type):
def get_open_polylines(input_filename):
def get_polygons(input_filename):
def get_all_fields(input_filename):
def get_attribute_names(input_filename):
def get_numeric_attribute_names(input_filename):
def get_points(input_filename, indices=None, with_z=False):
def write_shp_points_z(output_filename, z_name, points):
def write_shp_lines(output_filename, shape_type, lines, attribute_name): class SerafinData:
    def __init__(self, job_id, filename, language):
        self.job_id = job_id
        self.language = language
        self.filename = filename
        self.index = None
        self.triangles = {}
        self.header = None
        self.time = []  # <[float]>
        self.time_second = []  # <[datetime.timedelta]>  FIXME: should be renamed differently!
        self.start_time = None

        self.selected_vars = []
        self.selected_vars_names = {}
        self.selected_time_indices = []
        self.equations = []
        self.us_equation = None
        self.to_single = False

        self.operator = None
        self.metadata = {}
        # * select single layer: 'layer_selection': an integer
        # * vertical aggregation: 'vertical_operator': 'Min', 'Max' or 'Mean'

    def read(self):
        try:
            with Serafin.Read(self.filename, self.language) as input_stream:
                input_stream.read_header()
                input_stream.get_time()

                self.header = input_stream.header.copy()
                self.time = input_stream.time[:]
        except PermissionError:
            raise Serafin.SerafinRequestError('Permission denied (Is the file opened by another application?).')

        if self.header.date is not None:
            try:
                year, month, day, hour, minute, second = self.header.date
                self.start_time = datetime.datetime(year, month, day, hour, minute, second)
            except ValueError:
                logger.warning('Date seems invalid, replaced by default date.')
        if self.start_time is None:
            self.start_time = datetime.datetime(1900, 1, 1, 0, 0, 0)
        self.time_second = list(map(lambda x: datetime.timedelta(seconds=x), self.time))
        self.selected_vars = self.header.var_IDs[:]
        self.selected_vars_names = {var_id: (var_name, var_unit) for (var_id, var_name, var_unit)
                                    in zip(self.header.var_IDs, self.header.var_names, self.header.var_units)}
        self.selected_time_indices = list(range(len(self.time)))
        return self.header.is_2d

    def copy(self):
        copy_data = SerafinData(self.job_id, self.filename, self.language)
        copy_data.index = self.index
        copy_data.triangles = self.triangles
        copy_data.header = self.header
        copy_data.time = self.time
        copy_data.start_time = self.start_time
        copy_data.time_second = self.time_second
        copy_data.metadata = self.metadata

        copy_data.selected_vars = self.selected_vars[:]
        copy_data.selected_vars_names = deepcopy(self.selected_vars_names)
        copy_data.selected_time_indices = self.selected_time_indices[:]
        copy_data.equations = self.equations[:]
        copy_data.us_equation = self.us_equation
        copy_data.to_single = self.to_single
        copy_data.operator = self.operator
        return copy_data

    def default_output_header(self):
        output_header = self.header.copy()
        output_header.empty_variables()
        for var_ID in self.selected_vars:
            var_name, var_unit = self.selected_vars_names[var_ID]
            output_header.add_variable(var_ID, var_name, var_unit)
        if self.to_single:
            output_header.to_single_precision()
        return output_header

    def build_2d_output_header(self):
        output_header = self.header.copy_as_2d()
        output_header.empty_variables()
        for var_ID in self.selected_vars:
            var_name, var_unit = self.selected_vars_names[var_ID]
            output_header.add_variable(var_ID, var_name, var_unit)
        if self.to_single:
            output_header.to_single_precision()
        return output_header

    def transform_mesh(self, transformations):
        """!
        @brief Apply transformations on mesh nodes (only in 2D)
        @param transformations <[geom.transformation.Transformation]>: list of successive transformations
        """
        self.header.transform_mesh(transformations)
        self.index = None
        self.triangles = {} class PolylineData:
    def __init__(self):
        self.lines = []
        self.fields = []

    def __len__(self):
        return len(self.lines)

    def add_line(self, line):
        self.lines.append(line)

    def set_fields(self, fields):
        self.fields = fields[:]

    def is_empty(self):
        return len(self.lines) == 0

    def id_are_unique(self):
        ids = [line.id for line in self.lines]
        return len(ids) == len(set(ids)) class PointData:
    def __init__(self):
        self.points = []
        self.attributes = []
        self.fields = []
        self.fields_name = []
        self.attributes_decoded = []

    def __len__(self):
        return len(self.points)

    def add_point(self, point):
        self.points.append(point)

    def add_attribute(self, attribute):
        self.attributes.append(attribute)
        decoded = []
        for a in attribute:
            if type(a) == bytes:
                decoded.append(a.decode('latin-1'))
            decoded.append(str(a))
        self.attributes_decoded.append(decoded)

    def set_fields(self, fields):
        self.fields = fields[:]
        for f in fields:
            name = f[0]
            if type(name) == bytes:
                self.fields_name.append(name.decode('latin-1'))
            else:
                self.fields_name.append(name)

    def is_empty(self):
        return len(self.points) == 0 class CSVData:
    def __init__(self, filename, header=None, out_name='', separator=''):
        self.filename = filename
        self.out_name = ''
        self.metadata = {}
        self.separator = ''

        if header is None:  # read existing file
            self.separator = separator
            self.table = []
            self.out_name = out_name
            with open(out_name, 'r') as f:
                for line in f.readlines():
                    self.table.append(line.rstrip().split(self.separator))
        else:
            self.table = [header]

    def add_row(self, row):
        self.table.append(row)

    def write(self, filename, separator):
        with open(filename, 'w') as output_stream:
            for line in self.table:
                output_stream.write(separator.join(line))
                output_stream.write('\n')
        self.out_name = filename
        self.separator = separator class FluxCalculator:
    """!
    Compute flux across sections (integral along lines) from a Serafin input stream
    """

    LINE_INTEGRAL, DOUBLE_LINE_INTEGRAL, LINE_FLUX, AREA_FLUX, MASS_FLUX = 0, 1, 2, 3, 4

    def __init__(self, flux_type, var_IDs, input_stream, section_names, sections, time_sampling_frequency):
        self.flux_type = flux_type
        self.input_stream = input_stream
        self.section_names = section_names
        self.sections = sections

        self.var_IDs = var_IDs

        self.time_indices = range(0, len(input_stream.time), time_sampling_frequency)

        self.mesh = None
        self.intersections = []

    def construct_triangles(self, iter_pbar=lambda x, unit: x):
        """!
        Construct triangular elements (index construction)
        @param iter_pbar: iterable progress bar
        """
        self.mesh = TriangularVectorField(self.input_stream.header, True, iter_pbar)

    def construct_intersections(self):
        """!
        Construct the intersections between the mesh and all input sections
        """
        for section in self.sections:
            self.intersections.append(self.mesh.section_intersection(section))

    def flux_in_frame(self, intersections, values):
        """!
        @brief Do the flux computation in a single frame, depending on the flux type
        @param values <[numpy.1D-array]>: list of the values of the scalar/vector fields
        @return <float>: The value of the flux
        """
        if self.flux_type == FluxCalculator.LINE_INTEGRAL:
            return TriangularVectorField.line_integral(intersections, values[0])
        elif self.flux_type == FluxCalculator.DOUBLE_LINE_INTEGRAL:
            return TriangularVectorField.line_double_integral(intersections, values[0], values[1])
        elif self.flux_type == FluxCalculator.LINE_FLUX:
            return TriangularVectorField.line_flux(intersections, values[0], values[1])
        elif self.flux_type == FluxCalculator.AREA_FLUX:
            return TriangularVectorField.area_flux(intersections, values[0], values[1], values[2])
        else:
            return TriangularVectorField.mass_flux(intersections, values[0], values[1], values[2], values[3])

    def run(self, iter_pbar=lambda x, unit: x, fmt_float=settings.FMT_FLOAT):
        """!
        Separate the major part of the computation, allowing a GUI override
        @param iter_pbar: iterable progress bar
        """
        result = []
        for time_index in iter_pbar(self.time_indices, unit='frames'):
            i_result = [str(self.input_stream.time[time_index])]
            values = []
            for var_ID in self.var_IDs:
                values.append(self.input_stream.read_var_in_frame(time_index, var_ID))

            for j in range(len(self.sections)):
                intersections = self.intersections[j]
                flux = self.flux_in_frame(intersections, values)
                i_result.append(fmt_float.format(flux))
            result.append(i_result)
        return result

    def write_csv(self, result, output_stream, separator):
        output_stream.write('time')
        for name in self.section_names:
            output_stream.write(separator)
            output_stream.write(name)
        output_stream.write('\n')

        for line in result:
            output_stream.write(separator.join(line))
            output_stream.write('\n') class PossibleFluxComputation:
    """!
    @brief Determine which flux computations are possible from a set of variables and which flux type it is
    """
    COMMON_FLUXES = {
        'Cross sectional area (m²)': [('H',)],
        'Liquid flux (m³/s)': [('U', 'V', 'H'), ('I', 'J'), ('M', 'H'), ('Q',)],
        'Solid flux TOTAL (m³/s)': [('QSX', 'QSY'), ('QS',)],
        'Solid flux BEDLOAD (m³/s)': [('QSBLX', 'QSBLY'), ('QSBL',)],
        'Solid flux SUSPENSION (m³/s)': [('QSSUSPX', 'QSSUSPY'), ('QSSUSP',)],
    }

    def __init__(self, var_IDs, var_names):
        """!
        Determine all possible flux computations
        @param var_IDs <[str]>: variables identifiers
        @param var_names <[bytes]>: variable names
        """
        self.var_IDs = var_IDs
        self.var_names = var_names

    def __iter__(self):
        """!Iterate over all the possible flux computations"""
        # Add common liquid and solid fluxes
        for title, list_var_IDs in PossibleFluxComputation.COMMON_FLUXES.items():
            for var_IDs in list_var_IDs:
                if set(var_IDs).issubset(self.var_IDs):
                    yield '%s: (%s)' % (title, ', '.join(var_IDs))

        # Add specific solid fluxes for each sediment classes
        for name in self.var_names:
            str_name = name.decode(SLF_EIT).strip()
            if set(['U', 'V', 'H']).issubset(self.var_IDs):
                if 'TRACEUR' in str_name or 'TRACER' in str_name:
                    yield 'Solid flux (kg/s): (U, V, H, %s)' % str_name
            if 'QS CLASS' in str_name:
                yield 'Solid flux TOTAL (m³/s): (%s)' % str_name
            if 'QS BEDLOAD CL' in str_name:
                yield 'Solid flux BEDLOAD (m³/s): (%s)' % str_name
            if 'QS SUSP. CL' in str_name:
                yield 'Solid flux SUSPENSION (m³/s): (%s)' % str_name

    @staticmethod
    def get_variables(possible_flux):
        """!
        Get variables from a possible flux
        @param possible_flux <str>: text describing a possible flux
        @return <[str]>: list of var_IDs
        """
        return list(possible_flux.split(':')[1].split('(')[1][:-1].split(', '))

    @staticmethod
    def get_flux_type(var_IDs):
        """!
        Get flux type
        @param <[str]>: list of varIDs (between 1 and 4)
        @return <int>: flux type identifier
        """
        nb_vars = len(var_IDs)
        if nb_vars == 1:
            flux_type = FluxCalculator.LINE_INTEGRAL
        elif nb_vars == 2:
            if var_IDs[0] == 'M':
                flux_type = FluxCalculator.DOUBLE_LINE_INTEGRAL
            else:
                flux_type = FluxCalculator.LINE_FLUX
        elif nb_vars == 3:
            flux_type = FluxCalculator.AREA_FLUX
        elif nb_vars == 4:
            flux_type = FluxCalculator.MASS_FLUX
        else:
            NotImplementedError('Flux with more than 4 variables is not possible')
        logger.debug('Flux type %s' % flux_type)
        return flux_type

    @staticmethod
    def common_fluxes():
        """!
        @brief Get list of var_IDs corresponding to common fluxes
        """
        list_common_fluxes = []
        for _, list_var_IDs in PossibleFluxComputation.COMMON_FLUXES.items():
            for var_IDs in list_var_IDs:
                list_common_fluxes.append(list(var_IDs))
        return list_common_fluxes class TriangularVectorField(Mesh2D):
    """!
    @brief The representation of Mesh2D in Serafin file when one computes the flux across sections of some vectors

    The integral of scalars along a line can also be computed.

    The flux across a section (or integral along a line) is the sum of the flux (or integral) on
    all intersected segments between the mesh and the section.
    """
    def __init__(self, *args, **kwargs):
        super().__init__(*args, **kwargs)

    def section_intersection(self, section):
        """!
        @brief Return the intersections (normal vectors and interpolators) of the mesh with a open polyline
        @param section <geom.geometry.Polyline>: An open polyline
        @return <dict>: The list of tuples (normal vector, interpolator) of every intersected segments in triangles
        """
        intersections = {}
        potential_elements = self.get_intersecting_elements(section.bounds())
        for i, j, k in potential_elements:
            t = self.triangles[i, j, k]
            is_intersected, t_intersections = section.linestring_intersection(t)
            if is_intersected:
                interpolator = Interpolator(t)
                intersections[i, j, k] = []
                for intersection in t_intersections:
                    line = []   # the list of tuple (normal_vector, interpolator) for all start/end/turning points
                    prev_x, prev_y = None, None
                    for x, y in intersection.coords:
                        if prev_x is None:  # the first point doesn't have a normal vector
                            line.append(([0, 0], interpolator.get_interpolator_at(x, y)))
                        else:
                            line.append(([prev_y-y, x-prev_x],
                                         interpolator.get_interpolator_at(x, y)))
                        prev_x, prev_y = x, y
                    intersections[i, j, k].append(line)
        return intersections

    @staticmethod
    def line_integral(intersections, scalar_flux_values):
        """!
        @brief The line integral of a scalar field along a line (not really a flux)
        """
        flux = 0
        for i, j, k in intersections:  # iterating through triangles in the intersection
            f_vals = scalar_flux_values[[i, j, k]]

            for endpoints in intersections[i, j, k]:
                endpoints_f = []

                for normal, interpolator in endpoints:  # pre-compute the values on each start/turning/end points
                    endpoints_f.append(interpolator.dot(f_vals))

                for p in range(len(endpoints)-1):  # iterating through intersecting segments inside the triangle
                    flux += (endpoints_f[p] + endpoints_f[p+1]) * np.linalg.norm(endpoints[p+1][0])
        return flux / 2

    @staticmethod
    def line_double_integral(intersections, height_values, scalar_flux_values):
        """!
        @brief The line integral of a scalar field (product of two scalars) across a section (not really a flux)
        """
        flux = 0
        for i, j, k in intersections:  # iterating through triangles in the intersection
            h_vals = height_values[[i, j, k]]
            f_vals = scalar_flux_values[[i, j, k]]

            for endpoints in intersections[i, j, k]:
                endpoints_h, endpoints_f = [], []

                for normal, interpolator in endpoints:  # pre-compute the values on each start/turning/end points
                    endpoints_h.append(interpolator.dot(h_vals))
                    endpoints_f.append(interpolator.dot(f_vals))

                for p in range(len(endpoints)-1):  # iterating through intersecting segments inside the triangle
                    first_h, second_h = endpoints_h[p], endpoints_h[p+1]
                    first_f, second_f = endpoints_f[p], endpoints_f[p+1]

                    flux += (2 * (first_f * first_h + second_f * second_h)
                             + (first_f * second_h + second_f * first_h)) * np.linalg.norm(endpoints[p+1][0])

        return flux / 6

    @staticmethod
    def line_flux(intersections, x_vector_values, y_vector_values):
        """!
        @brief The flux of a vector field across a line
        """
        flux = 0
        for i, j, k in intersections:  # iterating through triangles in the intersection
            x_vals = x_vector_values[[i, j, k]]
            y_vals = y_vector_values[[i, j, k]]

            for endpoints in intersections[i, j, k]:
                endpoints_x, endpoints_y = [], []

                for normal, interpolator in endpoints:  # pre-compute the values on each start/turning/end points
                    endpoints_x.append(interpolator.dot(x_vals))
                    endpoints_y.append(interpolator.dot(y_vals))

                for p in range(len(endpoints)-1):  # iterating through intersecting segments inside the triangle
                    normal = endpoints[p+1][0]

                    first_normal = np.dot([endpoints_x[p], endpoints_y[p]], normal)
                    second_normal = np.dot([endpoints_x[p+1], endpoints_y[p+1]], normal)

                    flux += first_normal + second_normal

        return flux / 2

    @staticmethod
    def area_flux(intersections, x_vector_values, y_vector_values, height_values):
        """!
        @brief The flux of a vector field across a section (the surface height is a scalar field)
        """
        flux = 0
        for i, j, k in intersections:  # iterating through triangles in the intersection
            x_vals = x_vector_values[[i, j, k]]
            y_vals = y_vector_values[[i, j, k]]
            h_vals = height_values[[i, j, k]]

            for endpoints in intersections[i, j, k]:
                endpoints_x, endpoints_y, endpoints_h = [], [], []

                for normal, interpolator in endpoints:  # pre-compute the values on each start/turning/end points
                    endpoints_x.append(interpolator.dot(x_vals))
                    endpoints_y.append(interpolator.dot(y_vals))
                    endpoints_h.append(interpolator.dot(h_vals))

                for p in range(len(endpoints)-1):  # iterating through intersecting segments inside the triangle
                    normal = endpoints[p+1][0]

                    first_normal = np.dot([endpoints_x[p], endpoints_y[p]], normal)
                    second_normal = np.dot([endpoints_x[p+1], endpoints_y[p+1]], normal)

                    first_h, second_h = endpoints_h[p], endpoints_h[p+1]

                    flux += 2 * (first_normal * first_h + second_normal * second_h) \
                              + (first_normal * second_h + second_normal * first_h)

        return flux / 6

    @staticmethod
    def mass_flux(intersections, x_vector_values, y_vector_values, height_values, density_values):
        """!
        @brief The mass-flux of a vector field across a section (with surface height and density, two scalars fields)
        """
        flux = 0
        for i, j, k in intersections:  # iterating through triangles in the intersection
            x_vals = x_vector_values[[i, j, k]]
            y_vals = y_vector_values[[i, j, k]]
            h_vals = height_values[[i, j, k]]
            d_vals = density_values[[i, j, k]]

            for endpoints in intersections[i, j, k]:
                endpoints_x, endpoints_y, endpoints_h, endpoints_d = [], [], [], []

                for normal, interpolator in endpoints:  # pre-compute the values on each start/turning/end points
                    endpoints_x.append(interpolator.dot(x_vals))
                    endpoints_y.append(interpolator.dot(y_vals))
                    endpoints_h.append(interpolator.dot(h_vals))
                    endpoints_d.append(interpolator.dot(d_vals))

                for p in range(len(endpoints)-1):  # iterating through intersecting segments inside the triangle
                    normal = endpoints[p+1][0]

                    first_normal = np.dot([endpoints_x[p], endpoints_y[p]], normal)
                    second_normal = np.dot([endpoints_x[p+1], endpoints_y[p+1]], normal)

                    first_h, second_h = endpoints_h[p], endpoints_h[p+1]
                    first_d, second_d = endpoints_d[p], endpoints_d[p+1]

                    flux += 9 * (first_normal * first_h * first_d + second_normal * second_h * second_d) + \
                            (2*first_normal+second_normal) * (2*first_h+second_h) * (2*first_d+second_d) + \
                            (first_normal+2*second_normal) * (first_h+2*second_h) * (first_d+2*second_d)

        return flux / 72 class MeshInterpolator(Mesh2D):
    def __init__(self, *args, **kwargs):
        super().__init__(*args, **kwargs)

    def get_point_interpolators(self, points):
        nb_points = len(points)
        is_inside = [False] * nb_points
        point_interpolators = [None] * nb_points

        for index, (x, y) in enumerate(points):
            potential_elements = self.get_intersecting_elements((x, y, x, y))
            if not potential_elements:
                continue
            for i, j, k in potential_elements:
                t = self.triangles[i, j, k]
                is_in, point_interpolator = Interpolator(t).is_in_triangle(x, y)
                if is_in:
                    is_inside[index] = True
                    point_interpolators[index] = ((i, j, k), point_interpolator)
                    break

        return is_inside, point_interpolators

    def _get_line_interpolators(self, line):
        intersections = []
        internal_points = []  # line interpolators without intersections

        # record the distance offset before the first intersection point
        offset = 0
        found_intersection = False

        for right, up, segment in line.segments():  # for every segment, sort intersection points
            segment_intersections = []
            potential_elements = self.get_intersecting_elements(segment.bounds())
            for i, j, k in potential_elements:
                t = self.triangles[i, j, k]
                is_intersected, t_intersections = segment.linestring_intersection(t)
                if is_intersected:
                    interpolator = Interpolator(t)
                    for intersection in t_intersections:
                        for x, y in intersection.coords:
                            segment_intersections.append((x, y, (i, j, k), interpolator.get_interpolator_at(x, y)))

            # first sort by y, then sort by x
            if up:
                segment_intersections.sort(key=lambda x: x[1])
            else:
                segment_intersections.sort(key=lambda x: x[1], reverse=True)
            if right:
                segment_intersections.sort()
            else:
                segment_intersections.sort(reverse=True)

            intersections.extend(segment_intersections)
            if not segment_intersections:
                continue

            internal_points.append(segment_intersections[0])
            internal_points.append(segment_intersections[-1])

            if not found_intersection:
                first_point, second_point = list(segment.coords())
                if not segment_intersections:
                    offset += np.linalg.norm(np.array(second_point) - np.array(first_point))
                else:
                    found_intersection = True
                    first_intersection = np.array(segment_intersections[0][:2])
                    offset += np.linalg.norm(first_intersection - np.array(first_point))

        # if the intersection is continuous, every internal point or turning point has at least two duplicates
        prev_x, prev_y = None, None
        duplicates = 0
        to_remove = [False] * len(intersections)
        for i, (x, y, _, __) in enumerate(intersections):
            if i == 0:  # the start and end points are not duplicated
                continue
            if prev_x is None:
                prev_x, prev_y = x, y
                continue
            if x == prev_x and y == prev_y:
                to_remove[i] = True
                duplicates += 1
            else:
                if duplicates == 0:  # no duplicate found, the intersection is discontinuous
                    return [], [], [], []
                duplicates = 0
                prev_x, prev_y = x, y

        intersections = [intersections[i] for i in range(len(intersections)) if not to_remove[i]]

        # trim internal points from 2n+2 to n+1
        if internal_points:
            internal_points = internal_points[0:-1:2] + [internal_points[-1]]

        # compute cumulative distance
        distance = offset
        distances = [offset]
        for i in range(len(intersections)-1):
            first, second = intersections[i+1], intersections[i]
            distance += np.linalg.norm([second[0] - first[0], second[1] - first[1]])
            distances.append(distance)

        distance = offset
        distances_internal = [offset]
        for i in range(len(internal_points)-1):
            first, second = internal_points[i+1], internal_points[i]
            distance += np.linalg.norm([second[0] - first[0], second[1] - first[1]])
            distances_internal.append(distance)

        return intersections, distances, internal_points, distances_internal

    def get_line_interpolators(self, lines):
        nb_nonempty = 0
        indices_nonempty = []
        line_interpolators = []
        line_interpolators_internal = []

        for i, line in enumerate(lines):
            line_interpolator, distance, line_interpolator_internal, distance_internal = self._get_line_interpolators(line)

            if line_interpolator:
                nb_nonempty += 1
                indices_nonempty.append(i)

            line_interpolators.append((line_interpolator, distance))
            line_interpolators_internal.append((line_interpolator_internal, distance_internal))

        return nb_nonempty, indices_nonempty, line_interpolators, line_interpolators_internal

    @staticmethod
    def interpolate_along_lines(input_stream, selected_vars, selected_time_indices, indices_nonempty,
                                line_interpolators, fmt_float):
        for u, id_line in enumerate(indices_nonempty):
            line_interpolator, distances = line_interpolators[id_line]

            for v, time_index in enumerate(selected_time_indices):
                time_value = input_stream.time[time_index]
                var_values = []
                for var in selected_vars:
                    var_values.append(input_stream.read_var_in_frame(time_index, var))

                for (x, y, (i, j, k), interpolator), distance in zip(line_interpolator, distances):
                    row = [str(id_line+1), str(time_value), fmt_float.format(x), fmt_float.format(y),
                           fmt_float.format(distance)]

                    for i_var, var in enumerate(selected_vars):
                        values = var_values[i_var]
                        row.append(fmt_float.format(interpolator.dot(values[[i, j, k]])))
                    yield u, v, row

    @staticmethod
    def project_lines(input_stream, selected_vars, time_index, indices_nonempty, max_distance,
                      reference, line_interpolators, fmt_float):
        var_values = []
        for var in selected_vars:
            var_values.append(input_stream.read_var_in_frame(time_index, var))

        for u, id_line in enumerate(indices_nonempty):
            line_interpolator, _ = line_interpolators[id_line]

            for x, y, (i, j, k), interpolator in line_interpolator:
                d = reference.project(x, y)
                if d <= 0 or d >= max_distance:
                    continue
                row = [str(id_line+1), fmt_float.format(x), fmt_float.format(y), fmt_float.format(d)]
                for i_var, var in enumerate(selected_vars):
                    values = var_values[i_var]
                    row.append(fmt_float.format(interpolator.dot(values[[i, j, k]])))
                yield u, row class Serafin:
    """!
    @brief A Serafin object corresponds to a single Serafin in file IO stream

    # Attributes:
    - language <str>: language for variables ('fr' or 'en')
    - mode <str>: file open mode
    - file <_io.BufferedReader>: input Serafin stream
    - filename <str>: path to file
    - file_size <int>: size of file
    """
    def __init__(self, filename, mode, language):
        self.language = language
        self.mode = mode

        self.file = None  # will be opened when called using `with` statement
        self.filename = filename
        self.file_size = None

    def __enter__(self):
        self.file = open(self.filename, self.mode)
        return self

    def __exit__(self, exc_type, exc_val, exc_tb):
        self.file.close()
        return False def get_available_variables(input_variables, is_2d):
def get_necessary_equations(known_var_IDs, needed_var_IDs, is_2d, us_equation=None):
def do_calculations_in_frame(equations, input_serafin, time_index, selected_output_IDs,
                             output_float_type, is_2d, us_equation, ori_values={}): class TruncatedTriangularPrisms(Mesh2D):
    """!
    @brief The representation of Mesh2D in Serafin file when one computes the volume inside polygons of some variables

    A truncated triangular prism is a triangular prism with non-parallel bases.
    Its geometrical properties are entirely determined by the base triangles.

    A prism is constructed, on the top of the base triangle,
    by giving the height of the nodes in the upper triangle, called the values of the triangle.

    A prism can be intersected with a 2D polygon. In this case, only the base triangle is considered.

    The volume of the prism is always equal to the surface area of the base triangle times
    the interpolated value of the centroid of the upper triangle.
    When intersected with a 2D polygon, the volume of the intersection is equal to
    the surface area of the intersection (polygon or multipolygon) times
    the interpolated value of the centroid of the intersection.
    """
    def __init__(self, *args, **kwargs):
        super().__init__(*args, **kwargs)

    def polygon_intersection_strict(self, polygon):
        """!
        @brief Return the weight carried by the triangle nodes entirely contained in polygon
        @param polygon <geom.geometry.Polyline>: A polygon
        @return <numpy.1D-array>: The weight carried by the triangle nodes
        """
        potential_elements = self.get_intersecting_elements(polygon.bounds())
        weight = np.zeros((self.nb_points,), dtype=np.float64)
        for i, j, k in potential_elements:
            t = self.triangles[i, j, k]
            if polygon.contains(t):
                weight[[i, j, k]] += t.area
        return weight / 3.0

    def polygon_intersection(self, polygon):
        """!
        @brief Return the weight carried by the triangle nodes entirely contained in polygon and info about boundary triangles
        @param polygon <geom.geometry.Polyline>: A polygon
        @return <numpy.1D-array, dict>: The weight carried by the triangle nodes, and the dictionary of tuple (area, centroid value) for boundary triangles
        """
        potential_elements = self.get_intersecting_elements(polygon.bounds())
        weight = np.zeros((self.nb_points,), dtype=np.float64)
        triangle_polygon_intersection = {}
        for i, j, k in potential_elements:
            t = self.triangles[i, j, k]
            if polygon.contains(t):
                weight[[i, j, k]] += t.area
            else:
                is_intersected, intersection = polygon.polygon_intersection(t)
                if is_intersected:
                    centroid = intersection.centroid
                    interpolator = Interpolator(t).get_interpolator_at(centroid.x, centroid.y)
                    triangle_polygon_intersection[i, j, k] = (intersection.area, interpolator)
        return weight / 3.0, triangle_polygon_intersection

    def polygon_intersection_all(self, polygon):
        """!
        @brief Return all triangles entirely contained in polygon and all boundary triangle-polygon intersections
        @param polygon <geom.geometry.Polyline>: A polygon
        @return <dict, dict>: The dictionaries of all triangles contained in polygon, and of tuples (base triangle, intersection) for boundary triangles
        """
        potential_elements = self.get_intersecting_elements(polygon.bounds())
        weight = np.zeros((self.nb_points,), dtype=np.float64)
        triangles = {}
        triangle_polygon_net_intersection = {}
        triangle_polygon_intersection = {}
        for i, j, k in potential_elements:
            t = self.triangles[i, j, k]
            if polygon.contains(t):
                area = t.area
                vertices = tuple(map(np.array, list(t.exterior.coords)[:-1]))
                triangles[i, j, k] = (vertices, area)
                weight[[i, j, k]] += area
            else:
                is_intersected, intersection = polygon.polygon_intersection(t)
                if is_intersected:
                    vertices = tuple(map(np.array, list(t.exterior.coords)[:-1]))
                    area = t.area
                    centroid = intersection.centroid
                    interpolator = Interpolator(t).get_interpolator_at(centroid.x, centroid.y)
                    triangle_polygon_net_intersection[i, j, k] = (intersection.area, interpolator)
                    triangle_polygon_intersection[i, j, k] = (vertices, area, intersection)
        return weight / 3.0, triangle_polygon_net_intersection, triangles, triangle_polygon_intersection

    @staticmethod
    def boundary_volume_in_polygon(triangle_polygon_intersection, variable):
        """!
        @brief Return the total volume in all triangle-polygon intersections
        @param triangle_polygon_intersection <dict>: All triangle-polygon intersections defined by boundary triangles
        @param variable <numpy.1D-array>: The values of the variable on all nodes
        @return <numpy.float>: The total volume in all triangle-polygon intersections
        """
        volume = 0
        for i, j, k in triangle_polygon_intersection:
            area, interpolator = triangle_polygon_intersection[i, j, k]
            volume += area * interpolator.dot(variable[[i, j, k]])
        return volume

    @staticmethod
    def superior_prism_volume(vertices, area, values):
        """!
        @brief Return the volume in the half-space z > 0 of the prism with the given base triangle and values
        @param vertices <list>: The coordinates of the three vertices defining the base triangle
        @param area <float>: The area of the base triangle
        @param values <numpy.1D-array>: The values of the variable on the three nodes of the triangle
        @return <float>: The volume of the prism in the half-space z > 0
        """
        # eliminate special cases
        if min(values) >= 0:
            return area * sum(values) / 3.0
        elif max(values) <= 0:
            return 0

        p1, p2, p3 = vertices
        # remaining cases: triangle crosses the plane z = 0
        nb_points_superior = sum(values > 0)
        if nb_points_superior == 1:  # positive tetrahedron
            (z_bottom, _, _), (z_middle, _, _), (z_top, _, _) = sorted(zip(values, [1, 2, 3], [p1, p2, p3]))
            return area * z_top ** 3 / 3 / (z_top - z_middle) / (z_top - z_bottom)
        else:  # negative tetrahedron
            (z_bottom, _, bottom), (z_middle, _, middle), (z_top, _, top) = sorted(zip(values, [1, 2, 3], [p1, p2, p3]))
            volume_total = area * sum(values) / 3.0
            volume_negative = area * z_bottom ** 3 / 3 / (z_top - z_bottom) / (z_middle - z_bottom)
            return volume_total - volume_negative

    @staticmethod
    def superior_prism_volume_in_intersection(polygon, vertices, area, intersection, values):
        """!
        @brief Return the volume of the prism in the half-space z > 0 and inside the polygon
        @param polygon <shapely.geometry.Polygon>: The volume-defining polygon
        @param vertices <list>: The coordinates of the three vertices defining the base triangle
        @param area <float>: The area of the base triangle
        @param values <numpy.1D-array>: The values of the variable on the three nodes of the triangle
        @return <float>: The volume of the prism in the half-space z > 0 and inside the polygon
        """
        # eliminate special cases
        if min(values) >= 0:
            return area * sum(values) / 3.0
        elif max(values) <= 0:
            return 0
        p1, p2, p3 = vertices
        # remaining cases: triangle crosses the plane z = 0
        nb_points_superior = sum(values > 0)
        if nb_points_superior == 1:  # positive tetrahedron
            (z_bottom, _, bottom), (z_middle, _, middle), (z_top, _, top) = sorted(zip(values, [1, 2, 3], [p1, p2, p3]))
            top_middle_intersect = (z_top * middle - z_middle * top) / (z_top - z_middle)
            top_bottom_intersect = (z_top * bottom - z_bottom * top) / (z_top - z_bottom)
            new_triangle = geom.Polygon([top, top_middle_intersect, top_bottom_intersect])
            if polygon.contains(new_triangle):
                return new_triangle.area * z_top / 3.0
            else:
                is_intersected, intersection = polygon.polygon_intersection(new_triangle)
                if not is_intersected:
                    return 0
                centroid = intersection.centroid
                height = Interpolator(new_triangle).get_interpolator_at(centroid.x, centroid.y).dot(np.array([z_top, 0, 0]))
                return intersection.area * height
        else:  # negative tetrahedron
            (z_bottom, _, bottom), (z_middle, _, middle), (z_top, _, top) = sorted(zip(values, [1, 2, 3], [p1, p2, p3]))
            bottom_middle_intersect = (z_middle * bottom - z_bottom * middle) / (z_middle - z_bottom)
            bottom_top_intersect = (z_top * bottom - z_bottom * top) / (z_top - z_bottom)
            new_triangle = geom.Polygon([bottom, bottom_middle_intersect, bottom_top_intersect])
            volume_negative = new_triangle.area * z_bottom / 3.0
            centroid = intersection.centroid
            interpolator = Interpolator(new_triangle)
            new_values = np.array([z_bottom, 0, 0])
            height = interpolator.get_interpolator_at(centroid.x, centroid.y).dot(new_values)
            intersected_volume = intersection.area * height  # volume net intersected

            if polygon.contains(new_triangle):
                # volume sup = (volume net intersected) - (volume negative tetrahedron)
                return intersected_volume - volume_negative
            else:
                # volume sup = (volume net intersected)
                #              - (volume negative tetrahedron - volume in tetrahedron but not in polygon)
                is_intersected, new_difference = geometry.Polyline.triangle_difference(new_triangle, polygon)
                if not is_intersected:
                    return 0
                centroid = new_difference.centroid
                height = interpolator.get_interpolator_at(centroid.x, centroid.y).dot(new_values)
                volume_tetrahedron_not_in_polygon = new_difference.area * height
                return intersected_volume - (volume_negative - volume_tetrahedron_not_in_polygon) class VolumeCalculator:
    """!
    Compute volumes inside polygons from a Serafin input stream
    """

    NET_STRICT, NET, POSITIVE = 0, 1, 2
    INIT_VALUE = '+'  # special variable ID for the option 'Initial values of the first variable'

    def __init__(self, volume_type, var_ID, second_var_ID, input_stream, polynames, polygons,
                 time_sampling_frequency):
        self.volume_type = volume_type
        self.input_stream = input_stream
        self.polynames = polynames
        self.polygons = polygons

        self.var_ID = var_ID
        self.second_var_ID = second_var_ID

        self.time_indices = range(0, len(input_stream.time), time_sampling_frequency)

        self.mesh = None
        self.weights = []

        self.init_values = None
        if self.second_var_ID == VolumeCalculator.INIT_VALUE:
            self.init_values = input_stream.read_var_in_frame(0, self.var_ID)

    def construct_triangles(self, iter_pbar=lambda x, unit: x):
        self.mesh = TruncatedTriangularPrisms(self.input_stream.header, True, iter_pbar)

    def construct_weights(self, iter_pbar=lambda x, unit: x):
        """!
        Construct the point weights/intersections etc. depending on the volume type, for every polygons
        """
        if self.volume_type == VolumeCalculator.NET_STRICT:
            for poly in iter_pbar(self.polygons, unit='polygons'):
                self.weights.append(self.mesh.polygon_intersection_strict(poly))
        elif self.volume_type == VolumeCalculator.NET:
            for poly in iter_pbar(self.polygons, unit='polygons'):
                weight, triangle_polygon_intersection = self.mesh.polygon_intersection(poly)
                self.weights.append((weight, triangle_polygon_intersection))
        elif self.volume_type == VolumeCalculator.POSITIVE:
            for poly in iter_pbar(self.polygons, unit='polygons'):
                self.weights.append(self.mesh.polygon_intersection_all(poly))

    def volume_in_frame_in_polygon(self, weight, values, polygon):
        """!
        @brief Do the volume computation in a single frame, depending on the volume type
        @param weight <tuple>: the point weights/intersections etc. depending on the polygon and on the volume type
        @param values <numpy.1D-array>: the values of the variable for which the volume will be computed
        @param polygon <geom.geometry.Polygon>: the polygon in which the volume will be computed
        @return <float>: The value of the volume
        """
        if self.volume_type == VolumeCalculator.NET_STRICT:
            return weight.dot(values)
        elif self.volume_type == VolumeCalculator.NET:
            strict_weight, triangle_polygon_intersection = weight
            volume_inside = strict_weight.dot(values)
            volume_boundary = TruncatedTriangularPrisms.boundary_volume_in_polygon(triangle_polygon_intersection, values)
            return volume_inside + volume_boundary
        else:
            strict_weight, triangle_polygon_net_intersection, triangles, triangle_polygon_intersection = weight
            volume_net = strict_weight.dot(values)
            volume_net += TruncatedTriangularPrisms.boundary_volume_in_polygon(triangle_polygon_net_intersection, values)

            volume_positive = 0
            for a, b, c in triangles:
                vertices, area = triangles[a, b, c]
                volume_positive += TruncatedTriangularPrisms.superior_prism_volume(vertices, area, values[[a, b, c]])
            for a, b, c in triangle_polygon_intersection:
                vertices, area, intersection = triangle_polygon_intersection[a, b, c]
                volume_positive += TruncatedTriangularPrisms.superior_prism_volume_in_intersection(polygon, vertices,
                                                                                                   area, intersection,
                                                                                                   values[[a, b, c]])
            return volume_net, volume_positive, volume_net - volume_positive

    def read_values_in_frame(self, time_index):
        """!
        Read variable values in a single frame, depending on the first/second variable choice
        """
        values = self.input_stream.read_var_in_frame(time_index, self.var_ID)
        if self.second_var_ID is not None:
            if self.second_var_ID == VolumeCalculator.INIT_VALUE:
                values -= self.init_values
            else:
                second_values = self.input_stream.read_var_in_frame(time_index, self.second_var_ID)
                values -= second_values
        return values

    def run(self, fmt_float=settings.FMT_FLOAT):
        """!
        Separate the major part of the computation, allowing a GUI override
        """
        result = []
        for time_index in self.time_indices:
            i_result = [str(self.input_stream.time[time_index])]
            values = self.read_values_in_frame(time_index)

            for j in range(len(self.polygons)):
                weight = self.weights[j]
                volume = self.volume_in_frame_in_polygon(weight, values, self.polygons[j])
                if self.volume_type == VolumeCalculator.POSITIVE:
                    for v in volume:
                        i_result.append(fmt_float.format(v))
                else:
                    i_result.append(fmt_float.format(volume))
            result.append(i_result)
        return result

    def get_csv_header(self):
        header = ['time']
        if self.volume_type == VolumeCalculator.POSITIVE:
            for name in self.polynames:
                header.append(name)
                header.append(name + ' POSITIVE')
                header.append(name + ' NEGATIVE')
        else:
            for name in self.polynames:
                header.append(name)
        return header

    def write_csv(self, result, output_stream, separator):
        output_stream.write(separator.join(self.get_csv_header()))
        output_stream.write('\n')

        for line in result:
            output_stream.write(separator.join(line))
            output_stream.write('\n') def process_output_options(input_file, job_id, extension, suffix, in_source_folder, dir_path, double_name):
    input_path, input_name = os.path.split(input_file)
    input_rootname = os.path.splitext(input_name)[0]
    if double_name:
        output_name = input_rootname + '_' + job_id + suffix + extension
    else:
        output_name = input_rootname + suffix + extension
    if in_source_folder:
        filename = os.path.join(input_path, output_name)
    else:
        filename = os.path.join(dir_path, output_name)
    return filename def process_geom_output_options(input_file, job_id, extension, suffix, in_source_folder, dir_path, double_name):
    input_path, input_name = os.path.split(input_file)
    input_rootname = os.path.splitext(input_name)[0]
    if double_name:
        output_name = input_rootname + '_' + job_id + suffix + extension
    else:
        output_name = input_rootname + suffix + extension
    if in_source_folder:
        path = os.path.join(input_path, 'gis')
        if not os.path.exists(path):
            os.mkdir(path)
        filename = os.path.join(path, output_name)
    else:
        filename = os.path.join(dir_path, output_name)
    return filename def process_vtk_output_options(input_file, job_id, time_index, suffix, in_source_folder, dir_path, double_name):
    input_path, input_name = os.path.split(input_file)
    input_rootname = os.path.splitext(input_name)[0]
    if double_name:
        output_name = input_rootname + '_' + job_id + suffix + '_' + str(time_index) + '.vtk'
    else:
        output_name = input_rootname + suffix + '_' + str(time_index) + '.vtk'
    if in_source_folder:
        path = os.path.join(input_path, 'vtk')
        if not os.path.exists(path):
            os.mkdir(path)
        filename = os.path.join(path, output_name)
    else:
        filename = os.path.join(dir_path, output_name)
    return filename
from datetime import datetime
from multiprocessing import Process, Queue
from shapefile import ShapefileException
from shapely.geometry import Polygon
from pyteltools.conf import settings
from pyteltools.geom import BlueKenue, Shapefile
from pyteltools.slf.datatypes import SerafinData, PolylineData, PointData, CSVData
from pyteltools.slf.flux import FluxCalculator, PossibleFluxComputation, TriangularVectorField
from pyteltools.slf.interpolation import MeshInterpolator
from pyteltools.slf import Serafin
from pyteltools.slf.variables import do_calculations_in_frame, get_available_variables, \
    get_necessary_equations, new_variables_from_US
from pyteltools.slf.volume import TruncatedTriangularPrisms, VolumeCalculator
from .nodes_op import VerticalAggregationNode  # use only class constants VERTICAL_OPERATIONS
from .util import process_output_options, process_geom_output_options, process_vtk_output_options
import numpy as np
import os
import pyteltools.slf.misc as operations
        return self.done_queue.get()


def worker(input_queue, output_queue):
    for func, args in iter(input_queue.get, 'STOP'):
        result = func(*args)
        output_queue.put(result)


def success_message(node_name, job_id, info='', second_job_id=''):
    return '== %s - SUCCESS == %s (%s%s)%s' % (datetime.now().strftime('%Y-%m-%d %H:%M:%S'), node_name, job_id,
                                               ' AND %s' % second_job_id if second_job_id else '',
                                               ': %s.' % info if info else '')


def fail_message(reason, node_name, job_id, second_job_id=''):
    return '== %s - FAIL == %s (%s%s): %s.' % \
          (datetime.now().strftime('%Y-%m-%d %H:%M:%S'), node_name, job_id,
           ' AND %s' % second_job_id if second_job_id else '',
           reason)


def read_slf_2d(node_id, fid, filename, language, job_id):
    data = SerafinData(job_id, filename, language)
    try:
        is_2d = data.read()
    except (Serafin.SerafinRequestError, Serafin.SerafinValidationError) as e:
        return False, node_id, fid, data, fail_message(e.message, 'Load Serafin 2D', job_id)
    if not is_2d:
        success, message = False, fail_message('file is not 2D', 'Load Serafin 2D', job_id)
    else:
        success, message = True, success_message('Load Serafin 2D', job_id)
    return success, node_id, fid, data, message


def read_slf_3d(node_id, fid, filename, language, job_id):
    data = SerafinData(job_id, filename, language)
    try:
        is_2d = data.read()
    except (Serafin.SerafinRequestError, Serafin.SerafinValidationError) as e:
        return False, node_id, fid, data, fail_message(e.message, 'Load Serafin 3D', job_id)
    if is_2d:
        success, message = False, fail_message('file is not 3D', 'Load Serafin 3D', job_id)
    else:
        success, message = True, success_message('Load Serafin 3D', job_id)
    return success, node_id, fid, data, message


def read_slf_reference(node_id, filename, language):
    data = SerafinData('Ref', filename, language)
    try:
        is_2d = data.read()
    except (Serafin.SerafinRequestError, Serafin.SerafinValidationError) as e:
        return False, node_id, data, fail_message(e.message, 'Load Reference Serafin', 'all')
    if not is_2d:
        success, message = False, fail_message('file is not 2D', 'Load Reference Serafin', 'all')
    elif len(data.time) != 1:
        success, message = False, fail_message('the file has more than one frame', 'Load Reference Serafin', 'all')
    else:
        success, message = True, success_message('Load Reference Serafin', 'all')
    return success, node_id, data, message


def read_polygons(node_id, filename):
    try:
        with open(filename):
            pass
    except FileNotFoundError:
        message = fail_message('file not found', 'Load 2D Polygons', 'all')
        return False, node_id, None, message
    except PermissionError:
        message = fail_message('access denied', 'Load 2D Polygons', 'all')
        return False, node_id, None, message

    data = PolylineData()
    if filename.endswith('.i2s'):
        with BlueKenue.Read(filename) as f:
            f.read_header()
            for i, poly in enumerate(f.get_polygons()):
                poly.set_id('Polygon %i' % (i + 1))
                data.add_line(poly)
        data.set_fields(['Value'])
    else:
        try:
            for i, poly in enumerate(Shapefile.get_polygons(filename)):
                poly.set_id('Polygon %i' % (i + 1))
                data.add_line(poly)
            data.set_fields(Shapefile.get_all_fields(filename))
        except ShapefileException as e:
            message = fail_message(e, 'Load 2D Polygons', 'all')
            return False, node_id, None, message

    if data.is_empty():
        message = fail_message('the file does not contain any polygon', 'Load 2D Polygons', 'all')
        return False, node_id, None, message

    return True, node_id, data, success_message('Load 2D Polygons', 'all')


def read_polylines(node_id, filename):
    try:
        with open(filename):
            pass
    except FileNotFoundError:
        message = fail_message('file not found', 'Load 2D Open Polylines', 'all')
        return False, node_id, None, message
    except PermissionError:
        message = fail_message('access denied', 'Load 2D Open Polylines', 'all')
        return False, node_id, None, message

    data = PolylineData()
    if filename.endswith('.i2s'):
        with BlueKenue.Read(filename) as f:
            f.read_header()
            for i, poly in enumerate(f.get_open_polylines()):
                poly.set_id('Line %i' % (i + 1))
                data.add_line(poly)
        data.set_fields(['Value'])
    else:
        try:

            for i, poly in enumerate(Shapefile.get_open_polylines(filename)):