text:
import warnings

from scipy import sparse as sp
import numpy as np

from .utils import check_spmat, check_densemat

class Matrix:
    @property
    def shape(self):
        return self._data.shape

    @property
    def size(self):
        return self._data.size

class SparseMatrix(Matrix):
    """
    """
    def __init__(self, data, dtype):
        """
        """
        super().__init__()
        if data is None:
            self._data = sp.csr_matrix((0, 0), dtype=dtype)
        else:
            self._data = check_spmat(data, dtype=dtype)

        self.dtype = dtype

class DenseMatrix(Matrix):
    """
    """
    def __init__(self, data, dtype):
        """
        """
        super().__init__()
        if data is None:
            self._data = np.array([[]], dtype=dtype)
        else:
            self._data = check_densemat(data, dtype=dtype)

        self.dtype = dtype

class InteractionMatrix(SparseMatrix):
    """ matrix contains interaction between entity-entity
    """
    def __init__(self, data, is_implicit=True, transform_fn=None, dtype=np.float32):
        """
        Inputs:
            data (sparse matrix): input interaction
            is_implicit (bool): if interaction implicit. this case the each entry
                                of the matrix is the noisy proxy of the true interaciton
                                which assumes the binary {0, 1} coding.
            transform_fn (callable): data transform function when interaction is implicit
            dtype (np.dtype): numpy dtype for this interaction
        """
        SparseMatrix.__init__(self, data, dtype)

        # if true, observations are implicit proxy (confidence) for the binary data
        self.is_implicit = is_implicit
        # currently, only support sampled case for both implicit and implicit
        self.is_sampled = True
        self._transform_fn = transform_fn

    def transform(self):
        """ this is in-place method
        """
        if self.is_implicit:
            self._data = self._transform_fn(self._data)
        else:
            warnings.warn(
                '[Warning] sampled explicit interaction assume binary '
                'confidence. transformation is not applied.'
            )

    def transpose(self):
        """ transpose interaction matrix and return copied self
        """
        return InteractionMatrix(self._data.T.tocsr(),
                                 self.is_implicit,
                                 self._transform_fn,
                                 self.dtype)

    @property
    def is_sampled_explicit(self):
        return (not self.is_implicit) and self.is_sampled

class SparseFeatureMatrix(SparseMatrix):
    """
    """
    def __init__(self, data, dtype=np.float32):
        """
        """
        SparseMatrix.__init__(self, data, dtype)

class DenseFeatureMatrix(DenseMatrix):
    """
    """
    def __init__(self, data, dtype=np.float32):
        """
        """
        DenseMatrix.__init__(self, data, dtype)

max_stars_repo_name: eldrin/aarms

max_issues_repo_name: eldrin/aarms

max_forks_repo_name: eldrin/aarms

reason: import numpy,from scipy


text:
import curses

import numpy as np

from .widget import Widget, BORDER_STYLES

class ArrayWin(Widget):
    """
    A Widget whose state is stored an a numpy array and whose __getitem__ / __setitem__ use numpy indexing to update text.

    Other Parameters
    ----------------
    border: optional
        The border type; one of `nurses.widget.BORDER_STYLES`. (by default a widget has no border)
    border_color: optional
        A curses color_pair.  If a border is given, border_color will be the color of the border. (the default is `color`)

    Notes
    -----
    __getitem__ and __setitem__ call the respective buffer functions directly, so one can slice
    and write to a Widget as if it was a numpy array.
    """

    default_character = " "

    def __init__(self, *args, **kwargs):
        super().__init__(*args, **kwargs)
        self._buffer = None
        self._colors = None

    def update_geometry(self):
        if self.root is None:
            return

        super().update_geometry()

        if self._buffer is None:
            h, w = self.height, self.width
            self._buffer = np.full((h, w), self.default_character)
            self._colors = np.full((h, w), self.color)

            if self.has_border:
                self.border(self.border_style, self.border_color)

    @property
    def colors(self):
        return self._colors[1: -1, 1: -1] if self.has_border else self._colors

    @colors.setter
    def colors(self, array):
        if self.has_border:
            self._colors[1: -1, 1: -1] = array
        else:
            self._colors = array

    @property
    def buffer(self):
        return self._buffer[1: -1, 1: -1] if self.has_border else self._buffer

    @buffer.setter
    def buffer(self, array):
        if self.has_border:
            self._buffer[1: -1, 1: -1] = array
        else:
            self._buffer = array

    def _resize(self):
        if self.window is None:
            return

        if self.has_border:
            self._buffer[:, -1] = self._buffer[-1] = self.default_character  # Erase the right-most/bottom-most border in case widget expands

        height, width = self.height, self.width
        old_h, old_w = self._buffer.shape
        min_h, min_w = min(height, old_h), min(width, old_w)

        new_buffer = np.full((height, width), self.default_character)
        new_buffer[:min_h, :min_w] = self._buffer[:min_h, :min_w]

        new_colors = np.full((height, width), self.color)
        new_colors[:min_h, :min_w] = self._colors[:min_h, :min_w]

        self._buffer = new_buffer
        self._colors = new_colors

        super()._resize()

        if self.has_border:
            self.border(self.border_style, self.border_color)

    def push(self):
        """Write the buffers to the window.
        """
        it = np.nditer((self._buffer, self._colors), ["multi_index"])
        for char, color in it:

            # Newline character on the last corner of a window will advance the cursor out-of-bounds causing an error
            # TODO: Implement a more perfomant solution, either check the index, catch the error, or add an extra row to the windows
            if (ch := str(char)) == "\n":
                ch = " "

            y, x = it.multi_index
            self.window.addstr(y, x, ch, color)

    def refresh(self):
        self.push()
        super().refresh()

    def __getitem__(self, key):
        """
        `buffer.__getitem__` except offset if `self.has_border` is true
        (i.e., `buffer[1: -1, 1: -1].__getitem__` if `self.has_border`).
        """
        return self.buffer[key]

    def __setitem__(self, key, text):
        """
        Coerce `text` into a ndarray then call `buffer.__setitem__(key, text)`.

        Notes
        -----
        If `len(text) > 1`, `text` is coerced into an array or an array of arrays (depending on the presence of newlines).
        If the array's shape can't be cast to `self.buffer` it will be rotated and tried again (setting the text vertically).

        If `self.has_border` is true then indices will be offset automatically.
        (i.e., `buffer[1: -1, 1: -1].__setitem__` will be called instead)

        Examples
        --------
        >>> my_widget[2:4, :13] = "Hello, World!\\nI'm a widget!"
        """
        if "\n" in text:
            text = np.array(tuple(map(tuple, text.rstrip("\n").splitlines())))
        elif len(text) > 1:
            text = np.array(tuple(text))

        try:
            self.buffer[key] = text
        except ValueError:
            self.buffer[key] = np.rot90(text if len(text.shape) == 2 else text[None, ], -1)  # Try to fit the text vertically

    def border(self, style="light", color=None):
        """
        Draw a border on the edges of the widget.

        Parameters
        ----------
        style: optional
            The style of the border, can be one of `nurses.widget.BORDER_STYLES`. (the default is "light")

        color: optional
            The color of the border.  (the default is the widget's `color`)

        Notes
        -----
        Methods such as `__setitem__`, `roll`, `scroll`, `_resize` will take care to preserve the border
        as long as `has_border` is truth-y.  To disable this behavior set `has_border` to False or call
        this method with `read_only=False`.
        """
        if self._buffer is None:
            return

        self.border_style = style
        self.border_color = color

        ul, ur, v, h, ll, lr = BORDER_STYLES[style]

        b = self._buffer
        b[(0, -1), :] = h
        b[:, (0, -1)] = v
        b[ 0,  0] = ul
        b[ 0, -1] = ur
        b[-1,  0] = ll
        b[-1, -1] = lr

        c = self._colors
        c[0] = c[-1] = c[:, 0] = c[:, -1] = color or self.color

    def roll(self, shift=1, vertical=False):
        """
        Roll the contents of the widget. Items that roll beyond the last position are re-introduced at the first.

        Parameters
        ----------
        shift: optional
            Number of places to shift the contents. `shift` may be negative. (the default is 1)

        vertical: optional
            Whether to roll vertically.  (the default is `False`, i.e., rolls are horizontal by default)
        """
        axis = (-shift, 0) if vertical else (0, -shift)
        self.buffer = np.roll(self.buffer, axis, (0, 1))
        self.colors = np.roll(self.colors, axis, (0, 1))

    def scroll(self, lines=1):
        """
        Scroll the contents of the buffer upwards or downwards, erasing the last/first lines.

        Parameters
        ----------
        lines: optional
            Number of lines to scroll. To scroll down, lines should be negative. (the default is 1)
        """
        self.roll(lines, vertical=True)
        slice_ = slice(-lines, None) if lines > 0 else slice(None, -lines)
        self.buffer[slice_] = self.default_character
        self.colors[slice_] = self.color

ext: py

max_stars_repo_name: salt-die/nurses

max_issues_repo_name: salt-die/nurses

max_forks_repo_name: salt-die/nurses

reason: import numpy


text:
from typing import Dict, Tuple, Union

from collections import deque
import copy
from datetime import datetime
from enum import IntEnum
import io
import logging
import os

import numpy as np

from ..messages import *
from ..utils import trace

class MessageData(object):
    def __init__(self, message_type, params):
        self.message_type = message_type
        self.message_class = message_type_to_class.get(self.message_type, None)
        self.params = params
        self.messages = []

    def to_numpy(self, remove_nan_times: bool = True):
        """!
        @brief Convert the raw FusionEngine message data into numpy arrays that can be used for data analysis.

        On return, this function will add additional member variables to the class containing numpy arrays as returned
        by the `to_numpy()` function in the individual message data class (if supported).

        For example, if called on a @ref fusion_engine_client.messages.solution.PoseMessage "PoseMessage", this function
        will generate new members including `lla_deg`, a 3xN numpy array with the WGS-84 latitude, longitude, and
        altitude data for all available time epochs that can be used as follows:

        ```{.py}
        pose_data.to_numpy()
        mean_lla_deg = np.mean(pose_data.lla_deg, axis=1)
        ```

        @param remove_nan_times If `True`, remove entries whose P1 timestamps are `NaN` (if P1 time is available for
               this message type).
        """
        if hasattr(self.message_class, 'to_numpy'):
            have_cached_numpy_data = 'p1_time' in self.__dict__
            if have_cached_numpy_data:
                # If we don't have message data we can't do any conversion so the currently cached numpy data is as good
                # as it's gonna get. If it doesn't exist, so be it.
                if len(self.messages) == 0:
                    do_conversion = False
                else:
                    do_conversion = (len(self.messages) != len(self.p1_time) or
                                     float(self.messages[0].p1_time) != self.p1_time[0] or
                                     float(self.messages[-1].p1_time) != self.p1_time[-1])
            else:
                do_conversion = True

            if do_conversion:
                self.__dict__.update(self.message_class.to_numpy(self.messages))

                if remove_nan_times and 'p1_time' in self.__dict__:
                    is_nan = np.isnan(self.p1_time)
                    if np.any(is_nan):
                        keep_idx = ~is_nan
                        for key, value in self.__dict__.items():
                            if (key not in ('message_type', 'message_class', 'params', 'messages') and
                                    isinstance(value, np.ndarray)):
                                if len(value.shape) == 1:
                                    self.__dict__[key] = value[keep_idx]
                                elif len(value.shape) == 2:
                                    if value.shape[0] == len(is_nan):
                                        # Assuming first dimension is time.
                                        self.__dict__[key] = value[keep_idx, :]
                                    elif value.shape[1] == len(is_nan):
                                        # Assuming second dimension is time.
                                        self.__dict__[key] = value[:, keep_idx]
                                    else:
                                        # Unrecognized data shape.
                                        pass
                                else:
                                    # Unrecognized data shape.
                                    pass
        else:
            raise ValueError('Message type %s does not support numpy conversion.' %
                             MessageType.get_type_string(self.message_type))

class FileIndex(object):
    # Note: To reduce the index file size, we've made the following limitations:
    # - Fractional timestamp is floored so time 123.4 becomes 123. The data read should not assume that an entry's
    #   timestamp is its exact time
    RAW_DTYPE = np.dtype([('int', '<u4'), ('type', '<u2'), ('offset', '<u8')])

    DTYPE = np.dtype([('time', '<f8'), ('type', '<u2'), ('offset', '<u8')])

    @classmethod
    def load(cls, index_path):
        raw_data = np.fromfile(index_path, dtype=cls.RAW_DTYPE)
        return cls._from_raw(raw_data)

    @classmethod
    def save(cls, index_path, data: Union[np.ndarray, list]):
        if isinstance(data, np.ndarray) and data.dtype == cls.RAW_DTYPE:
            raw_data = data
            data = cls._from_raw(raw_data)
        else:
            if isinstance(data, list):
                data = np.array(data, dtype=cls.DTYPE)

            if data.dtype == cls.DTYPE:
                raw_data = cls._to_raw(data)
            else:
                raise ValueError('Unsupported array format.')

        raw_data.tofile(index_path)
        return data

    @classmethod
    def get_path(cls, data_path):
        return os.path.splitext(data_path)[0] + '.p1i'

    @classmethod
    def _from_raw(cls, raw_data):
        idx = raw_data['int'] == Timestamp._INVALID
        data = raw_data.astype(dtype=cls.DTYPE)
        data['time'][idx] = np.nan
        return data

    @classmethod
    def _to_raw(cls, data):
        time_sec = data['time']
        idx = np.isnan(time_sec)
        raw_data = data.astype(dtype=cls.RAW_DTYPE)
        raw_data['int'][idx] = Timestamp._INVALID
        return raw_data

class TimeAlignmentMode(IntEnum):
    NONE = 0
    DROP = 1
    INSERT = 2

class FileReader(object):
    logger = logging.getLogger('point_one.fusion_engine.analysis.file_reader')

    def __init__(self, path=None, regenerate_index=False):
        """!
        @brief Create a new reader instance.

        @param path The path to a FusionEngine binary file to open.
        """
        self.file = None
        self.file_size = 0
        self.data: Dict[MessageType, MessageData] = {}
        self.t0 = None
        self.system_t0 = None
        self.system_t0_ns = None

        self.index = None

        if path is not None:
            self.open(path, regenerate_index=regenerate_index)

    def open(self, path, regenerate_index=False):
        """!
        @brief Open a FusionEngine binary file.

        @param path The path to the file, or an existing Python file object.
        @param regenerate_index If `True`, delete the `.p1i` index file if it exists and regenerate it from the `.p1log`
               data file.
        """
        self.close()

        if isinstance(path, io.IOBase):
            self.file = path
        else:
            self.file = open(path, 'rb')

        self.file.seek(0, io.SEEK_END)
        self.file_size = self.file.tell()
        self.file.seek(0, 0)

        if self.file_size == 0:
            raise RuntimeError("File '%s' is empty." % path)

        # Load the data index file if present.
        index_path = FileIndex.get_path(self.file.name)
        have_index = os.path.exists(index_path)

        if have_index and regenerate_index:
            self.logger.debug("Deleting index file '%s'." % index_path)
            os.remove(index_path)
            have_index = False

        if have_index:
            self.logger.debug("Reading index file '%s'." % index_path)
            self.index = FileIndex.load(index_path)

            # If the index doesn't cover the full binary file, the user might have interrupted the read when it was
            # being generated. Delete it and create a new one.
            index_valid = True
            if len(self.index['time']) == 0:
                self.logger.warning("Index file '%s' is empty. Deleting." % index_path)
                index_valid = False
            else:
                last_offset = self.index['offset'][-1]
                if last_offset > self.file_size + MessageHeader.calcsize():
                    self.logger.warning("Last index entry past end of file. Deleting index.")
                    index_valid = False
                else:
                    # Read the header of the last entry to get its size.
                    self.file.seek(last_offset, io.SEEK_SET)
                    buffer = self.file.read(MessageHeader.calcsize())
                    self.file.seek(0, io.SEEK_SET)

                    header = MessageHeader()
                    header.unpack(buffer=buffer, warn_on_unrecognized=False)
                    message_size_bytes = MessageHeader.calcsize() + header.payload_size_bytes

                    index_size = last_offset + message_size_bytes
                    if index_size != self.file_size:
                        self.logger.warning("Size expected by index file does not match binary file. Deleting index.")
                        index_valid = False

            if not index_valid:
                os.remove(index_path)
                self.index = None
        else:
            self.index = None

        # Read the first message (with P1 time) in the file to set self.t0.
        #
        # Note that we explicitly set a start time since, if the time range is not specified, read() will include
        # messages that do not have P1 time. We want to make sure the 1 message is one with time.
        if self.index is None:
            self.t0 = None
            self.read(time_range=(0.0, None), require_p1time=True, max_messages=1, generate_index=False)
        else:
            idx = np.argmax(~np.isnan(self.index['time']))
            self.t0 = Timestamp(self.index['time'][idx])

        # Similarly, we also set the system t0 based on the first system-stamped (typically POSIX) message to appear in
        # the log, if any (profiling data, etc.). Unlike P1 time, since the index file does not contain system
        # timestamps, we have to do a read() even if an index exists. read() will use the index to at least speed up the
        # read operation.
        self.system_t0 = None
        self.system_t0_ns = None
        self.read(require_system_time=True, max_messages=1, ignore_index_max_messages=True, generate_index=False)

    def close(self):
        """!
        @brief Close the file.
        """
        if self.file is not None:
            self.file = None

    def generate_index(self):
        """!
        @brief Generate an index file for the current binary file if one does not already exist.
        """
        if self.index is None:
            # We'll read pose data (doesn't actually matter which). Store the currently cached data and restore it when
            # we're done. That way if the user already did a read (with generate_index == False), they don't have to
            # re-read the data if they try to use it again.
            prev_data = self.data.get(MessageType.POSE, None)

            self.read(message_types=[MessageType.POSE], max_messages=1, generate_index=True)

            if prev_data is not None:
                self.data[MessageType.POSE] = prev_data

    def read(self, message_types: Union[list, tuple, set] = None,
             time_range: Tuple[Union[float, Timestamp], Union[float, Timestamp]] = None, absolute_time: bool = False,
             max_messages: int = None, require_p1time: bool = False, require_system_time: bool = False,
             return_numpy: bool = False, keep_messages: bool = False, remove_nan_times: bool = True,
             time_align: TimeAlignmentMode = TimeAlignmentMode.NONE,
             aligned_message_types: Union[list, tuple, set] = None,
             generate_index: bool = True, show_progress: bool = False,
             ignore_index: bool = False, ignore_index_max_messages: bool = False) \
            -> Dict[MessageType, MessageData]:
        """!
        @brief Read data for one or more desired message types.

        The read data will be cached internally. Subsequent reads for the same data type will return the cached data.

        @note
        This function uses a data index file to speed up reads when available. If `generate_index == True` and no index
        file exists, one will be generated automatically. In order to do this, this function must read the entire data
        file, even if it could normally return early when `max_messages` or the end of `time_range` are reached.

        @param message_types A list of one or more @ref fusion_engine_client.messages.defs.MessageType "MessageTypes" to
               be returned. If `None` or an empty list, read all available messages.
        @param time_range An optional length-2 tuple specifying desired start and end bounds on the data timestamps.
               Both the start and end values may be set to `None` to read all data.
        @param absolute_time If `True`, interpret the timestamps in `time_range` as absolute P1 times. Otherwise, treat
               them as relative to the first message in the file.
        @param max_messages If set, read up to the specified maximum number of messages. Applies across all message
               types. If negative, read the last N messages.
        @param require_p1time If `True`, omit messages that do not contain valid P1 timestamps.
        @param require_system_time If `True`, omit messages that do not contain valid system timestamps.
        @param return_numpy If `True`, convert the results to numpy arrays for analysis.
        @param keep_messages If `return_numpy == True` and `keep_messages == False`, the raw data in the `messages`
               field will be cleared for each @ref MessageData object for which numpy conversion is supported.
        @param remove_nan_times If `True`, remove messages whose P1 timestamps are `NaN` when converting to numpy.
               Ignored if `return_numpy == False`.
        @param time_align The type of alignment to be performed (for data with P1 timestamps):
               - @ref TimeAlignmentMode.NONE - Do nothing
               - @ref TimeAlignmentMode.DROP - Drop messages at times when _all_ message types are not present
               - @ref TimeAlignmentMode.INSERT - Insert default-constructed messages for any message types not present
                 at a given time epoch
        @param aligned_message_types A list of message types for which time alignment will be performed. Any message
               types not present in the list will be left unmodified. If `None`, all message types will be aligned.
        @param show_progress If `True`, print the read progress every 10 MB (useful for large files).
        @param generate_index If `True` and an index file does not exist for this data file, read the entire data file
               and create an index file on the first call to this function. The file will be stored in the same
               directory as the input file.
        @param ignore_index If `True`, ignore the data index if loaded.
        @param ignore_index_max_messages If `True`, do not apply the `max_messages` limit when listing messages to be
               loaded from the data index file. `max_messages` _will_ still be applied after the data is decoded.

        @return A dictionary, keyed by @ref fusion_engine_client.messages.defs.MessageType "MessageType", containing
               @ref MessageData objects with the data read for each of the requested message types.
        """
        if message_types is None:
            message_types = []
        elif not isinstance(message_types, (list, tuple, set)):
            message_types = (message_types,)

        if time_range is None:
            time_range = (None, None)
            time_range_specified = False
        else:
            time_range_specified = (time_range[0] is not None or time_range[1] is not None)

        if max_messages is None:
            max_messages = 0

        params = {
            'time_range': time_range,
            'absolute_time': absolute_time,
            'max_messages': max_messages
        }

        # Allow the user to pass in a list of message classes for convenience and convert them to message types
        # automatically.
        message_types = [(t if isinstance(t, MessageType) else t.MESSAGE_TYPE) for t in message_types]

        # If the message type list is empty, read all messages.
        if len(message_types) == 0:
            message_types = list(message_type_to_class.keys())

        # If any of the requested types were already read from the file for the requested parameters, skip them.
        needed_message_types = [t for t in message_types
                                if (t not in self.data or self.data[t].params != params)]
        needed_message_types = set(needed_message_types)

        # Make cache entries for the messages to be read.
        message_class = {}
        for type in needed_message_types:
            self.data[type] = MessageData(message_type=type, params=params)

            message_class[type] = message_type_to_class.get(type, None)
            if message_class[type] is None:
                self.logger.warning('Decode not supported for message type %s. Omitting from output.' %
                                    MessageType.get_type_string(type))

        needed_message_types = [t for t in needed_message_types if message_class[t] is not None]

        system_time_messages_requested = any([t in messages_with_system_time for t in needed_message_types])

        # Create a dict with references to the requested types only to be returned below.
        result = {t: self.data[t] for t in message_types}

        num_needed = len(needed_message_types)
        if num_needed == 0:
            # Nothing to read. Return cached data.
            self.logger.debug('Requested data already cached. [# types=%d, start=%s, end=%s]' %
                              (len(message_types), str(time_range[0]), str(time_range[1])))
            return result
        elif self.file is None:
            raise IOError("File not open.")

        num_total = len(message_types)
        self.logger.debug('Reading data for %d message types. [cached=%d, start=%s, end=%s, max=%d]' %
                          (num_total, num_total - num_needed, str(time_range[0]), str(time_range[1]), max_messages))

        # Set the reference time used to compare timestamps against the user-specified time range. If we don't know t0
        # yet, it will be set later. If we already loaded an index file, t0 should have been set from that.
        if absolute_time:
            p1_reference_time_sec = 0.0
        elif self.t0 is not None:
            p1_reference_time_sec = float(self.t0)
        else:
            p1_reference_time_sec = None

        # For data which is timestamped in system time (CPU time, POSIX time, etc.) and not P1 time, we only support
        # relative time ranges (since absolute time is specified as P1 time, not system).
        if absolute_time:
            # Explicitly _not_ defining system_reference_time_sec. That way, if we attempt to use it below we'll hit
            # an exception.
            pass
        elif self.system_t0 is not None:
            system_reference_time_sec = self.system_t0
        else:
            system_reference_time_sec = None

        # If there's an index file, use it to determine the offsets to all the messages we're interested in.
        if self.index is not None and not ignore_index:
            type_idx = np.full_like(self.index['time'], False, dtype=bool)
            for type in needed_message_types:
                type_idx = np.logical_or(type_idx, self.index['type'] == type)

            # If t0 has never been set, this is probably the "first message" read done in open() to set t0. Ignore the
            # time range.
            if time_range_specified and self.t0 is not None:
                time_idx = np.full_like(self.index['time'], False, dtype=bool)
                limit_time = self.index['time'] - p1_reference_time_sec
                if time_range[0] is not None:
                    # Note: The index stores only the integer part of the timestamp.
                    time_idx = np.logical_and(time_idx, limit_time >= np.floor(time_range[0]))
                if time_range[1] is not None:
                    time_idx = np.logical_and(time_idx, limit_time <= time_range[1])

                # Messages with system time may have NAN timestamps in the index file since they can occur in a log
                # before P1 time is established. We'll allow any NAN times to pass this check, and we will decode and
                # filter them later based on their system timestamps.
                if system_time_messages_requested:
                    time_idx = np.logical_or(time_idx, np.isnan(self.index['time']))

                idx = np.logical_and(type_idx, time_idx)
            else:
                idx = type_idx

            data_index = self.index[idx]

            if not ignore_index_max_messages:
                if max_messages > 0:
                    data_index = data_index[:max_messages]
                elif max_messages < 0:
                    # If max is negative, take the last N entries.
                    data_index = data_index[max_messages:]

            generate_index = False
            data_offsets = data_index['offset']
        # Otherwise, seek to the start of the file and read all messages.
        else:
            data_offsets = None
            index_entries = []
            self.file.seek(0, 0)

            if generate_index:
                self.logger.debug('Reading all contents to generate index file.')

        # Read all messages meeting the criteria.
        HEADER_SIZE = MessageHeader.calcsize()

        total_bytes_read = 0
        bytes_used = 0
        message_count = 0
        index_count = 0

        if max_messages < 0:
            # Used for max_messages < 0 only.
            newest_messages = deque(maxlen=abs(max_messages))
        else:
            newest_messages = None

        last_print_bytes = 0
        start_time = datetime.now()
        while True:
            # Read the next message header. If we have an index, seek directly to the message.
            if data_offsets is not None:
                # All messages read.
                if index_count >= len(data_offsets):
                    break

                message_offset_bytes = data_offsets[index_count]
                self.file.seek(message_offset_bytes, io.SEEK_SET)
                index_count += 1
            else:
                message_offset_bytes = self.file.tell()

            buffer = self.file.read(HEADER_SIZE)
            if len(buffer) == 0:
                break

            # Deserialize the header.
            try:
                header = MessageHeader()
                header.unpack(buffer=buffer, warn_on_unrecognized=False)
                message_size_bytes = HEADER_SIZE + header.payload_size_bytes

                # Check if this is one of the message types we need. If not, continue to read the payload and then skip
                # the message below.
                try:
                    if data_offsets is not None:
                        message_needed = True
                    else:
                        message_needed = header.message_type.value in needed_message_types
                except AttributeError:
                    # If the message type was not recognized, header.message_type will be an int instead of MessageType
                    # enum. If it's not in MessageType, there's definitely no message class for it so skip it.
                    message_needed = False

                self.logger.debug('  Deserializing %s message @ %d. [length=%d B]%s' %
                                  (header.get_type_string(), message_offset_bytes, message_size_bytes,
                                   '' if message_needed else ' [skip]'))
            except Exception as e:
                self.logger.error('Error decoding header @ %d: %s' % (message_offset_bytes, repr(e)))
                break

            # Read the message payload and append it to the header.
            try:
                buffer += self.file.read(header.payload_size_bytes)
            except Exception as e:
                self.logger.error('Error reading %s payload @ %d: %s' %
                                  (header.get_type_string(), message_offset_bytes, repr(e)))
                break

            if len(buffer) != message_size_bytes:
                self.logger.warning('Unexpected EOF.')
                break

            # Validate the CRC.
            try:
                header.validate_crc(buffer)
            except ValueError as e:
                self.logger.error('Error reading %s message @ %d: %s' %
                                  (header.get_type_string(), message_offset_bytes, repr(e)))
                break

            total_bytes_read = self.file.tell()

            if total_bytes_read - last_print_bytes > 10e6:
                elapsed_sec = (datetime.now() - start_time).total_seconds()
                self.logger.log(logging.INFO if show_progress else logging.DEBUG,
                                '%sProcessed %d/%d bytes (%.1f%%). [elapsed=%.1f sec, rate=%.1f MB/s]' %
                                ('' if show_progress else '  ',
                                 total_bytes_read, self.file_size, 100.0 * float(total_bytes_read) / self.file_size,
                                 elapsed_sec, total_bytes_read / elapsed_sec / 1e6))
                last_print_bytes = total_bytes_read

            # If this isn't one of the requested messages, skip it. If we don't know t0 yet or we need the message's P1
            # time to generate an index, continue to unpack.
            if (not message_needed and
                self.t0 is not None and
                (not system_time_messages_requested or self.system_t0 is not None) and
                    not generate_index):
                continue

            # Now decode the payload.
            try:
                # Get the message class for this type and unpack the payload.
                cls = message_class.get(header.message_type, None)
                if cls is None:
                    is_reserved = int(header.message_type) >= int(MessageType.RESERVED)
                    self.logger.log(logging.DEBUG,
                                    '  Skipping message type %s @ %d. [needed=%s, reserved=%s]' %
                                    (header.get_type_string(), message_offset_bytes, repr(message_needed),
                                     repr(is_reserved)))
                    p1_time = None
                    contents = None
                else:
                    contents = cls()
                    contents.unpack(buffer=buffer, offset=HEADER_SIZE)

                    # Extract P1 and system times from this message, if applicable.
                    p1_time = contents.__dict__.get('p1_time', None)
                    system_time_ns = contents.__dict__.get('system_time_ns', None)
                    if system_time_ns is not None:
                        system_time_sec = system_time_ns * 1e-9
                    else:
                        system_time_sec = None

                # If we're building up an index file, add an entry for this message. If this is an unrecognized message
                # type, we won't have P1 time so we'll just insert a nan.
                if generate_index:
                    index_entries.append((float(p1_time) if p1_time is not None else np.nan, int(header.message_type),
                                          message_offset_bytes))

                if contents is None:
                    # Unrecognized type - skip.
                    continue

                # Store t0 if this is the first message with a (valid) timestamp.
                if p1_time is not None and p1_time:
                    if self.t0 is None:
                        self.logger.debug('Received first message. [type=%s, time=%s]' %
                                          (header.get_type_string(), str(p1_time)))
                        self.t0 = p1_time

                        if p1_reference_time_sec is None:
                            p1_reference_time_sec = float(p1_time)
                elif require_p1time:
                    continue

                if system_time_ns is not None:
                    if self.system_t0 is None:
                        self.logger.debug('Received first system-timestamped message. [type=%s, time=%.3f]' %
                                          (header.get_type_string(), system_time_sec))
                        self.system_t0 = system_time_sec
                        self.system_t0_ns = system_time_ns

                        if not absolute_time and system_reference_time_sec is None:
                            system_reference_time_sec = system_time_sec
                elif require_system_time:
                    continue

                # Now skip this message if we don't need it.
                if not message_needed:
                    continue

                # If this message has P1 time or system time and the user specified a time range, exclude it if it falls
                # outside the range. If it does not have time and a time range was specified, skip the message since we
                # can't be sure it's in the correct range.
                #
                # Note that system time does not support absolute time ranges.
                if time_range_specified:
                    # Select the appropriate timestamp and reference t0 value.
                    time_type = None
                    if (p1_time is None or not p1_time) and system_time_sec is None:
                        self.logger.debug('  Message does not contain time and time range specified. Discarding.')
                        continue
                    elif p1_time is not None and p1_time:
                        message_time_sec = p1_time
                        selected_ref_time_sec = p1_reference_time_sec
                        time_type = 'P1'
                    elif not absolute_time:
                        message_time_sec = system_time_sec
                        selected_ref_time_sec = system_reference_time_sec
                        time_type = 'system'
                    else:
                        self.logger.debug('  Message does not contain P1 time and absolute time range specified. '
                                          'Discarding.')
                        continue

                    # Now filter the message based on the selected time type.
                    time_offset_sec = float(message_time_sec) - selected_ref_time_sec
                    if time_range[0] is not None and time_offset_sec < float(time_range[0]):
                        self.logger.debug('  Message before requested time range. Discarding. [time=%s (%s)]' %
                                          (str(message_time_sec), time_type))
                        continue
                    elif time_range[1] is not None and time_offset_sec > float(time_range[1]):
                        # Assuming data is in order so if we pass the end of the time range we're done.
                        if generate_index:
                            self.logger.debug('  Message past requested time range. Discarding. [time=%s (%s)]' %
                                              (str(message_time_sec), time_type))
                            continue
                        else:
                            self.logger.debug('  Message past requested time range. Done reading. [time=%s (%s)]' %
                                              (str(message_time_sec), time_type))
                            break

                # Store the message.
                #
                # If we haven't reached the max message count yet, or it is 0 (disabled), store the message by type. If
                # max_messages is negative, the user requested the N latest messages. Add them to an N-length circular
                # buffer and pick off the newest ones when we're done with the file.
                message_count += 1
                if max_messages < 0:
                    newest_messages.append((header, contents))
                    bytes_used += message_size_bytes
                elif max_messages == 0 or message_count <= abs(max_messages):
                    self.data[header.message_type].messages.append(contents)
                    bytes_used += message_size_bytes

                if max_messages > 0:
                    if generate_index and message_count > max_messages:
                        self.logger.debug('  Max messages reached. Discarding. [# messages=%d]' % message_count)
                        continue
                    elif not generate_index and message_count == max_messages:
                        self.logger.debug('  Max messages reached. Done reading. [# messages=%d]' % message_count)
                        break
            except Exception as e:
                # Since the CRC passed we know we read the correct number of bytes, so if the decode fails because of
                # a format version mismatch or something, we can simply skip this message. No need to stop reading the
                # file.
                self.logger.warning('Error decoding %s payload @ %d: %s' %
                                    (header.get_type_string(), message_offset_bytes, repr(e)))

        # If the user only wanted the N newest messages, take them from the circular buffer now.
        if max_messages < 0:
            for entry in newest_messages:
                self.data[entry[0].message_type].messages.append(entry[1])

        # Print a summary.
        end_time = datetime.now()
        elapsed_sec = (end_time - start_time).total_seconds()
        self.logger.log(logging.INFO if show_progress else logging.DEBUG,
                        '%sRead %d bytes, used %d bytes. [elapsed=%s sec]' %
                        ('' if show_progress else '  ',
                         total_bytes_read, bytes_used,
                         ('%.1f' if elapsed_sec > 1.0 else '%e') % elapsed_sec))

        # Save the index file for faster reading in the future.
        if generate_index:
            index_path = FileIndex.get_path(self.file.name)
            self.logger.debug("Saving index file '%s' with %d entries." % (index_path, len(index_entries)))
            self.index = FileIndex.save(index_path, index_entries)

        # Time-align the data if requested.
        FileReader.time_align_data(result, mode=time_align, message_types=aligned_message_types)

        # Convert the resulting message data to numpy (if supported).
        if return_numpy:
            FileReader.to_numpy(result, keep_messages=keep_messages, remove_nan_times=remove_nan_times)

        # Done.
        return result

    def get_system_t0(self):
        return self.system_t0

    def get_system_t0_ns(self):
        return self.system_t0_ns

    @classmethod
    def time_align_data(cls, data: dict, mode: TimeAlignmentMode = TimeAlignmentMode.INSERT,
                        message_types: Union[list, tuple, set] = None):
        """!
        @brief Time-align messages of different types.

        @post
        `data` will be modified in-place. Message types that do not contain P1 time will be left unmodified.

        @param data A data `dict` as returned by @ref read().
        @param mode The type of alignment to be performed:
               - @ref TimeAlignmentMode.NONE - Do nothing
               - @ref TimeAlignmentMode.DROP - Drop messages at times when _all_ message types are not present
               - @ref TimeAlignmentMode.INSERT - Insert default-constructed messages for any message types not present
                 at a given time epoch
        @param message_types A list of message types for which alignment will be performed. Any message types not
               present in the list will be left unmodified. If `None`, all message types will be aligned.

        @return A modified `dict` with removed or inserted messages.
        """
        # Time alignment disabled - do nothing.
        if mode == TimeAlignmentMode.NONE:
            return data

        if message_types is not None:
            # Allow the user to pass in a list of message classes for convenience and convert them to message types
            # automatically.
            message_types = set([(t if isinstance(t, MessageType) else t.MESSAGE_TYPE) for t in message_types])

        # Pull out the P1 times for each message type. In drop mode, compute the intersection of all P1 timestamps. In
        # insert mode, make a list of all unique P1 timestamps.
        info_by_type = {}
        time_set = None
        for type, entry in data.items():
            default = entry.message_class()
            if 'p1_time' in default.__dict__ and (message_types is None or entry.message_type in message_types):
                p1_time = np.array([float(m.p1_time) for m in entry.messages])
                info_by_type[type] = {'p1_time': p1_time, 'messages': entry.messages, 'class': entry.message_class}

                if mode == TimeAlignmentMode.DROP:
                    if time_set is None:
                        time_set = p1_time
                    else:
                        time_set = np.intersect1d(time_set, p1_time)
                else:
                    if time_set is None:
                        time_set = p1_time
                    else:
                        time_set = np.hstack((time_set, p1_time))

        # In insertion mode, insert default-constructed objects for any missing timestamps.
        if mode == TimeAlignmentMode.INSERT:
            time_set = np.unique(time_set)

            for type, entry in info_by_type.items():
                # Locate the timestamps where we do/do not have data. For timestamps with data we store the index of the
                # corresponding message. For timestamps without we store -1.
                _, idx, all_idx = np.intersect1d(entry['p1_time'], time_set, return_indices=True)
                message_indices = np.full_like(time_set, -1, dtype=int)
                message_indices[all_idx] = idx

                # Now interlace messages with defaults as needed.
                messages = entry['messages']
                cls = entry['class']

                def _get_value(i):
                    message_idx = message_indices[i]
                    if message_idx >= 0:
                        return messages[message_idx]
                    else:
                        default = cls()
                        default.p1_time = time_set[i]
                        return default
                data[type].messages = [_get_value(i) for i in range(len(time_set))]
        # In drop mode, drop messages that aren't present across _all_ message types.
        elif mode == TimeAlignmentMode.DROP:
            for type, entry in info_by_type.items():
                _, idx, _ = np.intersect1d(entry['p1_time'], time_set, return_indices=True)
                data[type].messages = [entry['messages'][i] for i in idx]
        else:
            raise ValueError('Unrecognized alignment mode.')

        return data

    @classmethod
    def to_numpy(cls, data: dict, keep_messages: bool = True, remove_nan_times: bool = True):
        """!
        @brief Convert all (supported) messages in a data dictionary to numpy for analysis.

        See @ref MessageData.to_numpy().

        @param data A data `dict` as returned by @ref read().
        @param keep_messages If `False`, the raw data in the `messages` field will be cleared for each @ref MessageData
               object for which numpy conversion is supported.
        @param remove_nan_times If `True`, remove entries whose P1 timestamps are `NaN` (if P1 time is available for
               this message type).
        """
        for entry in data.values():
            try:
                entry.to_numpy(remove_nan_times=remove_nan_times)
                if not keep_messages:
                    entry.messages = []
            except ValueError:
                pass

max_stars_repo_name: PointOneNav/fusion-engine-client

max_issues_repo_name: PointOneNav/fusion-engine-client

max_forks_repo_name: PointOneNav/fusion-engine-client


text:
#! python3
# Runs a Dixon tissue segmentation for each case in the database.

from __future__ import print_function
from GetMetricFromElastixRegistration import GetFinalMetricFromElastixLogFile
from MultiAtlasSegmentation import MultiAtlasSegmentation
from ApplyBiasCorrection import ApplyBiasCorrection
import SimpleITK as sitk
import SitkImageManipulation as sitkIm
import DixonTissueSegmentation
import winshell
import numpy as np
import sys
import os

############################### CONFIGURATION #####################################
DEBUG = 1 # In debug mode, all the intermediate iamges are written.
USE_COSINES_AND_ORIGIN = 1
OVERWRITE_EXISTING_SEGMENTATIONS = 1

############################### TARGET FOLDER ###################################
# The target is the folder where the MRI images to be processed are. In the folder only
# folders with the case name should be found. Inside each case folder there must be a subfolder
# named "ForLibrary" with the dixon images called "case_I, case_O, case_W, case_F".
#targetPath = 'D:\\Martin\\Data\\MuscleSegmentation\\MarathonStudy\\PreMarathon\\AllWithLinks\\'
#targetPath = 'D:\\Martin\\Data\\MuscleSegmentation\\MarathonStudy\\PostMarathon\\AllWithLinks\\'
#targetPath = 'D:\\Martin\\Data\\MuscleSegmentation\\MarathonStudy\\PostMarathon\\NotSegmented\\'
#targetPath = 'D:\\Martin\\Data\\MuscleSegmentation\\DixonFovOK\\'
#targetPath = 'D:\\Martin\\Data\\MuscleSegmentation\\DixonFovOkTLCCases2020\\'
targetPath = 'D:\\Martin\\Data\\MuscleSegmentation\\MuscleStudyHipSpine\\CouchTo5kStudy\\'

subFolder = '\\ForLibraryNoCropping\\'
subFolder = '\\ForLibraryCropped\\'

# Cases to process, leave it empty to process all the cases in folder:
casesToSegment = ('C00007', 'C00019', 'C00020', 'C00057', 'C00077')
casesToSegment = ('C00024')
#casesToSegment = list()
# Look for the folders or shortcuts:
files = os.listdir(targetPath)
# It can be lnk with shortcuts or folders:
extensionShortcuts = 'lnk'
strForShortcut = '-> '
extensionImages = 'mhd'
inPhaseSuffix = '_I'#
outOfPhaseSuffix = '_O'#
waterSuffix = '_W'#
fatSuffix = '_F'#
suffixSegmentedImages = '_tissue_segmented'
suffixSkinFatImages = '_skin_fat'
suffixMuscleImages = '_muscle'
suffixFatFractionImages = '_fat_fraction'
dixonSuffixInOrder = (inPhaseSuffix, outOfPhaseSuffix, waterSuffix, fatSuffix)
for filenameInDir in files:
    dixonImages = []
    name, extension = os.path.splitext(filenameInDir)
    if (len(casesToSegment) == 0) or (name in casesToSegment):
        # if name is a lnk, get the path:
        if str(extension).endswith(extensionShortcuts):
            # This is a shortcut:
            shortcut = winshell.shortcut(targetPath + filenameInDir)
            indexStart = shortcut.as_string().find(strForShortcut)
            dataPath = shortcut.as_string()[indexStart+len(strForShortcut):] + '\\'
        else:
            dataPath = targetPath + filenameInDir + '\\'
        # Check if the images are available:
        filename = dataPath + subFolder + name + inPhaseSuffix + '.' + extensionImages
        outFilenameSegmented = dataPath + subFolder + name + suffixSegmentedImages + '.' + extensionImages
        outFilenameFatFraction = dataPath + subFolder + name + suffixFatFractionImages + '.' + extensionImages
        if (OVERWRITE_EXISTING_SEGMENTATIONS) or (not os.path.exists(outFilenameSegmented) and not os.path.exists(outFilenameFatFraction)):
            if os.path.exists(filename):
                # Process this image:
                print('Image to be processed: {0}\n'.format(name))
                # Add images in order:
                for suffix in dixonSuffixInOrder:
                    filename = dataPath + subFolder + name + suffix + '.' + extensionImages
                    dixonImages.append(sitk.Cast(sitk.ReadImage(filename), sitk.sitkFloat32))

                # Generate teh Dixon tissue image:
                segmentedImage = DixonTissueSegmentation.DixonTissueSegmentation(dixonImages)
                # Write image:
                sitk.WriteImage(segmentedImage, outFilenameSegmented, True)

                # Now create a skin fat mask:
                #skinFat = DixonTissueSegmentation.GetSkinFatFromTissueSegmentedImage(segmentedImage)
                #sitk.WriteImage(skinFat,
                #                dataPath + subFolder + name + suffixSkinFatImages + '.' + extensionImages, True)

                muscleMask = DixonTissueSegmentation.GetMuscleMaskFromTissueSegmentedImage(segmentedImage, vectorRadius = (4,4,3))
                sitk.WriteImage(muscleMask,
                                dataPath + subFolder + name + suffixMuscleImages + '.' + extensionImages, True)

max_stars_repo_name: mabelzunce/MuscleSegmentation

max_issues_repo_name: mabelzunce/MuscleSegmentation

max_forks_repo_name: mabelzunce/MuscleSegmentation


text:
from lenstronomy.LensModel.LightConeSim.light_cone import LightCone
from lenstronomy.LensModel.lens_model import LensModel
from astropy.cosmology import FlatLambdaCDM
from lenstronomy.Cosmo.lens_cosmo import LensCosmo
from lenstronomy.Util import util
from lenstronomy.LensModel import convergence_integrals
import numpy.testing as npt
import numpy as np

class TestLightCone(object):

    def setup(self):
        # define a cosmology
        cosmo = FlatLambdaCDM(H0=70, Om0=0.3, Ob0=0.05)
        self._cosmo = cosmo
        redshift_list = [0.1, 0.3, 0.8]  # list of redshift of the deflectors
        z_source = 2  # source redshift
        self._z_source = z_source
        # analytic profile class in multi plane
        self._lensmodel = LensModel(lens_model_list=['NFW', 'NFW', 'NFW'], lens_redshift_list=redshift_list,
                                    multi_plane=True, z_source_convention=z_source, cosmo=cosmo, z_source=z_source)
        # a single plane class from which the convergence/mass maps are computeded
        single_plane = LensModel(lens_model_list=['NFW'], multi_plane=False)
        # multi-plane class with three interpolation grids
        self._lens_model_interp = LensModel(lens_model_list=['INTERPOL', 'INTERPOL', 'INTERPOL'], lens_redshift_list=redshift_list,
                                    multi_plane=True, z_source_convention=z_source, cosmo=cosmo, z_source=z_source)
        # deflector parameterisation in units of reduced deflection angles to the source convention redshift
        logM_200_list = [8, 9, 10]  # log 10 halo masses of the three deflectors
        c_list = [20, 10, 8]  # concentrations of the three halos
        kwargs_lens = []
        kwargs_lens_interp = []
        grid_spacing = 0.01  # spacing of the convergence grid in units arc seconds
        x_grid, y_grid = util.make_grid(numPix=500, deltapix=grid_spacing)  # we create the grid coordinates centered at zero
        x_axes, y_axes = util.get_axes(x_grid, y_grid)  # we need the axes only for the interpolation
        mass_map_list = []
        grid_spacing_list_mpc = []
        for i, z in enumerate(redshift_list):  # loop through the three deflectors
            lens_cosmo = LensCosmo(z_lens=z, z_source=z_source, cosmo=cosmo)  # instance of LensCosmo, a class that manages cosmology relevant quantities of a lens
            alpha_Rs, Rs = lens_cosmo.nfw_physical2angle(M=10**(logM_200_list[i]), c=c_list[i])  # we turn the halo mass and concentration in reduced deflection angles and angles on the sky
            kwargs_nfw = {'Rs': Rs, 'alpha_Rs': alpha_Rs, 'center_x': 0, 'center_y': 0}  # lensing parameters of the NFW profile in lenstronomy conventions
            kwargs_lens.append(kwargs_nfw)
            kappa_map = single_plane.kappa(x_grid, y_grid, [kwargs_nfw])  # convergence map of a single NFW profile
            kappa_map = util.array2image(kappa_map)
            mass_map = lens_cosmo.sigma_crit_angle * kappa_map * grid_spacing ** 2  # projected mass per pixel on the gird
            mass_map_list.append(mass_map)
            npt.assert_almost_equal(np.log10(np.sum(mass_map)), logM_200_list[i], decimal=0)  # check whether the sum of mass roughtly correspoonds the mass definition
            grid_spacing_mpc = lens_cosmo.arcsec2phys_lens(grid_spacing)  # turn grid spacing from arcseconds into Mpc
            grid_spacing_list_mpc.append(grid_spacing_mpc)
            f_x, f_y = convergence_integrals.deflection_from_kappa_grid(kappa_map, grid_spacing)  # perform the deflection calculation from the convergence map
            f_ = convergence_integrals.potential_from_kappa_grid(kappa_map, grid_spacing)  # perform the lensing potential calculation from the convergence map (attention: arbitrary normalization)
            kwargs_interp = {'grid_interp_x': x_axes, 'grid_interp_y': y_axes, 'f_': f_, 'f_x': f_x, 'f_y': f_y}  # keyword arguments of the interpolation model
            kwargs_lens_interp.append(kwargs_interp)
        self.kwargs_lens = kwargs_lens
        self.kwargs_lens_interp = kwargs_lens_interp
        self.lightCone = LightCone(mass_map_list, grid_spacing_list_mpc, redshift_list)  # here we make the instance of the LightCone class based on the mass map, physical grid spacing and redshifts.

    def test_ray_shooting(self):
        beta_x, beta_y = self._lensmodel.ray_shooting(2., 1., self.kwargs_lens)
        beta_x_num, beta_y_num = self._lens_model_interp.ray_shooting(2., 1., self.kwargs_lens_interp)
        npt.assert_almost_equal(beta_x_num, beta_x, decimal=1)
        npt.assert_almost_equal(beta_y_num, beta_y, decimal=1)
        lens_model, kwargs_lens = self.lightCone.cone_instance(z_source=self._z_source, cosmo=self._cosmo, multi_plane=True)
        assert len(lens_model.lens_model_list) == 3
        beta_x_cone, beta_y_cone = lens_model.ray_shooting(2., 1., kwargs_lens)
        npt.assert_almost_equal(kwargs_lens[0]['grid_interp_x'], self.kwargs_lens_interp[0]['grid_interp_x'], decimal=5)
        npt.assert_almost_equal(kwargs_lens[0]['grid_interp_y'], self.kwargs_lens_interp[0]['grid_interp_y'], decimal=5)

        npt.assert_almost_equal(kwargs_lens[0]['f_x'], self.kwargs_lens_interp[0]['f_x'], decimal=5)
        npt.assert_almost_equal(kwargs_lens[0]['f_y'], self.kwargs_lens_interp[0]['f_y'], decimal=5)

        npt.assert_almost_equal(beta_x_cone, beta_x_num, decimal=5)
        npt.assert_almost_equal(beta_y_cone, beta_y_num, decimal=5)

    def test_deflection(self):
        alpha_x, alpha_y = self._lensmodel.alpha(2, 1, self.kwargs_lens)
        alpha_x_num, alpha_y_num = self._lens_model_interp.alpha(2, 1, self.kwargs_lens_interp)
        npt.assert_almost_equal(alpha_x_num, alpha_x, decimal=3)
        npt.assert_almost_equal(alpha_y_num, alpha_y, decimal=3)
        lens_model, kwargs_lens = self.lightCone.cone_instance(z_source=self._z_source, cosmo=self._cosmo, multi_plane=True)
        alpha_x_cone, alpha_y_cone = lens_model.alpha(2, 1, kwargs_lens)
        npt.assert_almost_equal(alpha_x_cone, alpha_x, decimal=3)
        npt.assert_almost_equal(alpha_y_cone, alpha_y, decimal=3)

    def test_arrival_time(self):
        x = np.array([1, 1])
        y = np.array([0, 1])
        f_ = self._lensmodel.arrival_time(x, y, self.kwargs_lens)
        f_num = self._lens_model_interp.arrival_time(x, y, self.kwargs_lens_interp)
        npt.assert_almost_equal(f_num[0] - f_num[1], f_[0] - f_[1], decimal=1)
        lens_model, kwargs_lens = self.lightCone.cone_instance(z_source=self._z_source, cosmo=self._cosmo, multi_plane=True)
        f_cone = lens_model.arrival_time(x, y, kwargs_lens)
        npt.assert_almost_equal(f_cone[0] - f_cone[1], f_[0] - f_[1], decimal=1)

reason: import numpy,from astropy


text:
import tensorflow as tf
import numpy as np
from dias.model.FPN_Model import Dias_FPN
from dias.model.Unet_Model import Dias_Unet
from dias.dataIO.data import IonoDataManager
import segmentation_models as sm
import matplotlib.pyplot as plt
import os

def train(cfgs):
    """
    Function for training DIAS.
    """
    print('Setting Model...')
    if cfgs['Model']['Type'] == 'Unet' or cfgs['Model']['Type'] == 'naiveUnet':
        model = Dias_Unet(cfgs)
    elif cfgs['Model']['Type'] == 'FPN':
        model = Dias_FPN(cfgs)
    dataManager = IonoDataManager(cfgs)
    # check data IO
    total_step = int(cfgs['Train']['TotalStep'])
    # print interval
    print_interval = int(cfgs['Train']['PrintInterval'])
    # plot interval
    plot_inverval = int(cfgs['Train']['PlotInterval'])
    # save interval
    save_inverval = int(cfgs['Train']['SaveInterval'])
    # plot image save directory
    img_save_dir = cfgs['Train']['ImgSaveDir']
    if os.path.exists(img_save_dir):
        pass
    else:
        os.makedirs(img_save_dir)
    # history log directory
    hist_log_dir = cfgs['Train']['HistLogDir']
    if os.path.exists(hist_log_dir):
        pass
    else:
        os.makedirs(hist_log_dir)
    # model save directory
    model_save_dir = cfgs['Train']['ModelSaveDir']
    if os.path.exists(model_save_dir):
        pass
    else:
        os.makedirs(model_save_dir)
    hist_list = list()
    print('Start Training')
    # start training
    for idx in range(total_step):
        x_train, y_train = dataManager.get_train_batch()
        hist = model.train_on_batch(x=x_train,y=y_train)
        hist_list.append(hist)

        # Print history information
        if idx % print_interval == 0:
            print('On Step {}'.format(idx))
            print(hist)

        # Save image result
        if idx % plot_inverval == 0:
            x_train, y_train = dataManager.get_train_batch()
            y_test = model.predict(x_train)
            plt.figure(figsize=(24,8))
            plt.subplot(1,3,1)
            plt.imshow(x_train[0,:,:,:])
            plt.subplot(1,3,2)
            plt.imshow(y_train[0,:,:,:])
            plt.subplot(1,3,3)
            plt.imshow(y_test[0,:,:,:])
            plt.savefig(img_save_dir+'STEP_{}.png'.format(idx),dpi=300)
            plt.close()

        # Save model
        if idx % save_inverval == 0:
            model.save(model_save_dir+'STEP_{}.model'.format(idx))
            np.save(model_save_dir+'hist.npy',hist_list)

    return

max_stars_repo_name: MrXiaoXiao/DIAS

max_issues_repo_name: MrXiaoXiao/DIAS

max_forks_repo_name: MrXiaoXiao/DIAS


text:
# -*- coding: utf-8 -*-
""" Example of calculation of g-functions using equal inlet temperatures.

    The g-functions of a field of 6x4 boreholes are calculated for boundary
    conditions of (a) uniform heat extraction rate, equal for all boreholes,
    (b) uniform borehole wall temperature along the boreholes, equal for all
    boreholes, and (c) equal inlet fluid temperature into all boreholes.

"""
from __future__ import division, print_function, absolute_import

import matplotlib.pyplot as plt
from matplotlib.ticker import AutoMinorLocator
import numpy as np
from scipy import pi

import pygfunction as gt

def main():
    # -------------------------------------------------------------------------
    # Simulation parameters
    # -------------------------------------------------------------------------

    # Borehole dimensions
    D = 4.0             # Borehole buried depth (m)
    H = 150.0           # Borehole length (m)
    r_b = 0.075         # Borehole radius (m)
    B = 7.5             # Borehole spacing (m)

    # Pipe dimensions
    rp_out = 0.0211     # Pipe outer radius (m)
    rp_in = 0.0147      # Pipe inner radius (m)
    D_s = 0.052         # Shank spacing (m)
    epsilon = 1.0e-6    # Pipe roughness (m)

    # Pipe positions
    # Single U-tube [(x_in, y_in), (x_out, y_out)]
    pos_pipes = [(-D_s, 0.), (D_s, 0.)]

    # Ground properties
    alpha = 1.0e-6      # Ground thermal diffusivity (m2/s)
    k_s = 2.0           # Ground thermal conductivity (W/m.K)

    # Grout properties
    k_g = 1.0           # Grout thermal conductivity (W/m.K)

    # Pipe properties
    k_p = 0.4           # Pipe thermal conductivity (W/m.K)

    # Fluid properties
    m_flow = 0.25       # Total fluid mass flow rate per borehole (kg/s)
    cp_f = 3977.        # Fluid specific isobaric heat capacity (J/kg.K)
    den_f = 1015.       # Fluid density (kg/m3)
    visc_f = 0.00203    # Fluid dynamic viscosity (kg/m.s)
    k_f = 0.492         # Fluid thermal conductivity (W/m.K)

    # Number of segments per borehole
    nSegments = 12

    # Geometrically expanding time vector.
    dt = 100*3600.                  # Time step
    tmax = 3000. * 8760. * 3600.    # Maximum time
    Nt = 50                         # Number of time steps
    ts = H**2/(9.*alpha)            # Bore field characteristic time
    time = gt.utilities.time_geometric(dt, tmax, Nt)

    # -------------------------------------------------------------------------
    # Borehole field
    # -------------------------------------------------------------------------

    # Field of 6x4 (n=24) boreholes
    N_1 = 6
    N_2 = 4
    boreField = gt.boreholes.rectangle_field(N_1, N_2, B, B, H, D, r_b)

    # -------------------------------------------------------------------------
    # Initialize pipe model
    # -------------------------------------------------------------------------

    # Pipe thermal resistance
    R_p = gt.pipes.conduction_thermal_resistance_circular_pipe(rp_in,
                                                               rp_out,
                                                               k_p)
    # Fluid to inner pipe wall thermal resistance (Single U-tube)
    h_f = gt.pipes.convective_heat_transfer_coefficient_circular_pipe(m_flow,
                                                                      rp_in,
                                                                      visc_f,
                                                                      den_f,
                                                                      k_f,
                                                                      cp_f,
                                                                      epsilon)
    R_f = 1.0/(h_f*2*pi*rp_in)

    # Single U-tube, same for all boreholes in the bore field
    UTubes = []
    for borehole in boreField:
        SingleUTube = gt.pipes.SingleUTube(pos_pipes, rp_in, rp_out,
                                           borehole, k_s, k_g, R_f + R_p)
        UTubes.append(SingleUTube)

    # -------------------------------------------------------------------------
    # Evaluate the g-functions for the borefield
    # -------------------------------------------------------------------------

    # Calculate the g-function for uniform heat extraction rate
    gfunc_uniform_Q = gt.gfunction.uniform_heat_extraction(
            boreField, time, alpha, disp=True)

    # Calculate the g-function for uniform borehole wall temperature
    gfunc_uniform_T = gt.gfunction.uniform_temperature(
            boreField, time, alpha, nSegments=nSegments, disp=True)

    # Calculate the g-function for equal inlet fluid temperature
    gfunc_equal_Tf_in = gt.gfunction.equal_inlet_temperature(
            boreField, UTubes, m_flow, cp_f, time, alpha,
            nSegments=nSegments, disp=True)

    # -------------------------------------------------------------------------
    # Plot g-functions
    # -------------------------------------------------------------------------

    plt.rc('figure')
    fig = plt.figure()
    ax1 = fig.add_subplot(111)
    # g-functions
    ax1.plot(np.log(time/ts), gfunc_uniform_Q,
             'k-', lw=1.5, label='Uniform heat extraction rate')
    ax1.plot(np.log(time/ts), gfunc_uniform_T,
             'b--', lw=1.5, label='Uniform borehole wall temperature')
    ax1.plot(np.log(time/ts), gfunc_equal_Tf_in,
             'r-.', lw=1.5, label='Equal inlet temperature')
    ax1.legend()
    # Axis labels
    ax1.set_xlabel(r'$ln(t/t_s)$')
    ax1.set_ylabel(r'$g(t/t_s)$')
    # Axis limits
    ax1.set_xlim([-10.0, 5.0])
    ax1.set_ylim([0., 50.])
    # Show minor ticks
    ax1.xaxis.set_minor_locator(AutoMinorLocator())
    ax1.yaxis.set_minor_locator(AutoMinorLocator())
    # Adjust to plot window
    plt.tight_layout()

    return

# Main function
if __name__ == '__main__':
    main()

max_stars_repo_name: mitchute/pygfunction

max_issues_repo_name: mitchute/pygfunction

max_forks_repo_name: mitchute/pygfunction


text:
import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.ticker as mtick
import matplotlib.patches as patches

from .config import config
from . import sv_html_formatters
from .sv_types import FeatureType, FeatureToProcess
from . import graph
from matplotlib.ticker import PercentFormatter

class GraphLegend(graph.Graph):
    def __init__(self, dataframe_report):
        styles = ["graph_base.mplstyle"]
        self.set_style(styles)

        fig = plt.figure(
            figsize=(config["Graphs"].getfloat("legend_width"), config["Graphs"].getfloat("legend_height")))
        axs = fig.add_axes([0, 0, 1, 1])
        axs.axis('off')
        scale = axs.transAxes.transform((1,1))
        scale = [ 1.0 / x for x in scale ]
        def to_fractions(pos_in_pix):
            return (pos_in_pix[0] * scale[0], pos_in_pix[1] * scale[1])

        def to_fractionsxy(x, y):
            return (x * scale[0], y * scale[1])

        def to_fraction_seq(seq, scalar):
            return [x * scalar for x in seq]

        cycle_colors = plt.rcParams['axes.prop_cycle'].by_key()['color']
        label_color = plt.rcParams['axes.labelcolor']

        gfx_x_source = 225
        gfx_x_compare = 350
        gfx_bar_y = 14
        gfx_line_y = 5
        bar_size = np.array([25,6])
        bar_text_offset = np.array([6,0])
        line_text_offset = np.array([0,-3])

        axs.add_patch(patches.Rectangle(to_fractionsxy(gfx_x_source, gfx_bar_y), bar_size[0] * scale[0], bar_size[1] * scale[1], facecolor=cycle_colors[0]))
        text1 = np.array([gfx_x_source,gfx_bar_y]) - bar_text_offset
        text1[1] += 1
        text1_elem = plt.text(text1[0] * scale[0], text1[1] * scale[1], "" + dataframe_report.source_name, fontsize=8, color=cycle_colors[0],
                              ha='right')

        # COMPARE value
        if dataframe_report.compare_name:
            axs.add_patch(
                patches.Rectangle(to_fractionsxy(gfx_x_compare, gfx_bar_y), bar_size[0] * scale[0], bar_size[1] * scale[1],
                                  facecolor=cycle_colors[1]))
            text2 = np.array([gfx_x_compare, gfx_bar_y]) + bar_text_offset
            text2[0] += bar_size[0]
            text2_elem = plt.text(text2[0] * scale[0], text2[1] * scale[1],
                                  "" + dataframe_report.compare_name, fontsize=8, color=cycle_colors[1])

        # TARGETS
        if dataframe_report.get_target_type() is not None:
            axs.add_line(
                matplotlib.lines.Line2D(to_fraction_seq([gfx_x_source, gfx_x_source + bar_size[0]], scale[0]),
                                        to_fraction_seq([gfx_line_y, gfx_line_y], scale[1]), lw=(1),
                                        color= graph.COLOR_TARGET_SOURCE, marker='o' ))
            text1[1] = gfx_line_y
            text1 += line_text_offset
            if dataframe_report.get_target_type() == FeatureType.TYPE_NUM:
                text_content = "Avg. " + dataframe_report._target["name"]
            else:
                text_content = "% " + dataframe_report._target["name"]
            text1_elem = plt.text(text1[0] * scale[0], text1[1] * scale[1], text_content, fontsize=8, color= graph.COLOR_TARGET_SOURCE,
                                  ha='right')

            if dataframe_report.compare_name and dataframe_report._target.get("compare"):
                axs.add_line(
                    matplotlib.lines.Line2D(
                        to_fraction_seq([gfx_x_compare, gfx_x_compare + bar_size[0]], scale[0]),
                        to_fraction_seq([gfx_line_y, gfx_line_y], scale[1]), lw=(1),
                        color= graph.COLOR_TARGET_COMPARE, marker='o'))

                text2[1] = gfx_line_y
                text2 += line_text_offset
                if dataframe_report.get_target_type() == FeatureType.TYPE_NUM:
                    text_content = "Avg. " + dataframe_report._target["name"]
                    #+ f" ({dataframe_report.compare_name})"
                else:
                    text_content = "% " + dataframe_report._target["name"]
                    #+ f" ({dataframe_report.compare_name})"
                text2_elem = plt.text(text2[0] * scale[0], text2[1] * scale[1], text_content, fontsize=8, color= graph.COLOR_TARGET_COMPARE)

        # transf = axs.transData.inverted()
        # bb = text1_elem.get_window_extent(renderer=fig.canvas.renderer)
        # bb_datacoords = bb.transformed(transf)

        self.graph_base64 = self.get_encoded_base64(fig)
        plt.close('all')
        return

max_stars_repo_name: mehrankamal/GML

max_issues_repo_name: mehrankamal/GML

max_forks_repo_name: mehrankamal/GML


text:
import numpy as np
import argparse, os
from utils import load_data, make_path
from pca import pca_train, pca_test
from svd import svd_train, svd_test
from isomap import isomap
from knn import knn

parser = argparse.ArgumentParser(description='Dimensionality Reduction')

parser.add_argument('--mode', type=str, default='train', help='train/test')
parser.add_argument('--train_data', type=str, default='sonar', help='train data source')
parser.add_argument('--test_data', type=str, default='sonar', help='test data source')
parser.add_argument('--alg', type=str, default='pca', help='pca/svd/isomap')
parser.add_argument('--output', type=str, default='output', help='result output dir')
parser.add_argument('--dim', type=int, default= 10, help='target dim')

args = parser.parse_args()
make_path(args.output)

train_data_path = os.path.join('data', args.train_data + '-train.txt')
test_data_path = os.path.join('data', args.test_data + '-test.txt')

def train():
    train_data, _ = load_data(train_data_path)
    if args.alg == 'pca':
        mean_vecs, eig_vecs = pca_train(np.array(train_data), args.dim)
        np.savez(os.path.join(args.output, 'pca' + str(args.dim) + '.npz'), mean_vecs=mean_vecs, eig_vecs=eig_vecs)
    elif args.alg == 'svd':
        v = svd_train(np.array(train_data), args.dim)
        np.savez(os.path.join(args.output, 'svd' + str(args.dim) + '.npz'), v=v)
    else:
        pass

def test():
    test_data, test_labels = load_data(test_data_path)
    train_data, train_labels = load_data(train_data_path)

    if args.alg == 'pca':
        saved_data = np.load(os.path.join(args.output, args.alg + str(args.dim) + '.npz'))
        mean_vecs = saved_data['mean_vecs']
        eig_vecs = saved_data['eig_vecs']

        reduction_test_data = pca_test(np.array(test_data), mean_vecs, eig_vecs)
        reduction_train_data = pca_test(np.array(train_data), mean_vecs, eig_vecs)
    elif args.alg == 'svd':
        saved_data = np.load(os.path.join(args.output, args.alg + str(args.dim) + '.npz'))
        v = saved_data['v']

        reduction_test_data = svd_test(np.array(test_data), v)
        reduction_train_data = svd_test(np.array(train_data), v)
    else:
        merge_data = train_data + test_data
        result = isomap(np.array(merge_data), n=args.dim)
        reduction_test_data = result[len(train_data):]
        reduction_train_data = result[:len(train_data)]

    acc = eval(reduction_test_data, test_labels, reduction_train_data, train_labels)
    print('ACC = ' + str(acc))

def eval(test_data, test_data_label, train_data, train_data_label):
    predict_label = []
    for x in test_data:
        label =  knn(x, train_data, train_data_label, 1)
        predict_label.append(label)

    right = 0
    size = len(test_data_label)
    for i in range(size):
        if predict_label[i] == test_data_label[i]:
            right = right + 1
    return right / size

if args.mode == 'train':
    train()
elif args.mode == 'test':
    test()

max_stars_repo_name: njucjc/machine-learning

max_issues_repo_name: njucjc/machine-learning

max_forks_repo_name: njucjc/machine-learning


text:
import cPickle
import numpy as np

def unpickle(file):
    fo = open(file, 'rb')
    dict = cPickle.load(fo)
    fo.close()
    return dict

def GetLabelsAndData(data_batch):
    y = data_batch['labels']
    X = data_batch['data']
    return y, X

def ZeroMean(X, N):
    b = np.mean(X,0)
    i=0
    while i < N:
        X[i] = X[i] - b
        i=i+1
    return X

max_stars_repo_path: SPN-DAG Project/Important/src/Data.py

max_stars_repo_name: atcbosselut/spn-discriminative

max_issues_repo_path: SPN-DAG Project/Important/src/Data.py

max_issues_repo_name: atcbosselut/spn-discriminative

max_forks_repo_path: SPN-DAG Project/Important/src/Data.py

max_forks_repo_name: atcbosselut/spn-discriminative


text:
import numpy
import matplotlib.pylab as plt

class LinearMaterialsParallel(object):

    """
    Solves the system:

        \div \frac{\rho}{\mu} k \grad (p + \rho g z) = 0  on the domain [x_0, x_1] \cross [z_0,z_1]

    Boundary conditions are given by:

       h(x_0,z,t) = h_0 [m]
       h(x_1,z,t) = h_L [m]

    Parameters are in units of:

       \rho : density, [kg/m^3]
       \mu  : viscosity, [kg / m s^2]
       K    : absolute permeability, [ m^2 ]
       g    : gravity, used in converting head to pressure, [ m / s^2 ]

    """

    def __init__(self, params=None):
        if params is None:
            params = dict()

        params.setdefault("x_0",0)
        params.setdefault("x_1",100)
        params.setdefault("z_0",0)
        params.setdefault("z_1",10)

        params.setdefault("k1",1.1847e-12)
        params.setdefault("k2",1.1847e-11)
        params.setdefault("rho",998.2)
        params.setdefault("mu",1.002e-3)

        params.setdefault("h_0",20.0)
        params.setdefault("h_L",19.0)

        params.setdefault("g",9.80665)
        params.setdefault("p_atm",101325.0)

        self.__dict__.update(params)

    def head(self, coords):

        """
        Compute the head at the x-values given by coords[:]

          h(x) = (h_L - h_0)*(x/L) + h0, i = 1,2

        """

        head = numpy.zeros(len(coords))
        # Compute the head at the interface.
        h0 = self.h_0
        hL = self.h_L
        L = self.x_1 - self.x_0
        K1 = self.rho * self.g * self.k1 / self.mu
        K2 = self.rho * self.g * self.k2 / self.mu
        for i in range(len(coords[:,0])):
            x = coords[i,0]
            head[i] = (hL - h0)*(x/L) + h0
        return head

def createFromXML(filename):

    # grab params from input file
    params = dict()

    import amanzi_xml.utils.io
    xml = amanzi_xml.utils.io.fromFile(filename)
    import amanzi_xml.utils.search as search

    #
    #  Domain Size
    #
    xyz = search.find_tag_path(xml, ["amanzi_input","mesh","generate","box",]).get("low_coordinates")
    params["x_0"] = float(xyz.split(',')[0])
    params["z_0"] = float(xyz.split(',')[2])
    xyz = search.find_tag_path(xml, ["amanzi_input","mesh","generate","box",]).get("high_coordinates")
    params["x_1"] = float(xyz.split(',')[0])
    params["z_1"] = float(xyz.split(',')[2])

    #
    #  Material Properties
    #
    strK = search.find_tag_path(xml, ["amanzi_input","materials","material,Front Material","permeability",]).get('x')
    params["k1"] = float(strK)
    strK = search.find_tag_path(xml, ["amanzi_input","materials","material,Back Material","permeability",]).get('x')
    params["k2"] = float(strK)
    strMu = search.find_tag_path(xml, ["amanzi_input","phases","liquid_phase","viscosity",]).text
    params["mu"] = float(strMu)
    strRho = search.find_tag_path(xml, ["amanzi_input","phases","liquid_phase","density",]).text
    params["rho"] = float(strRho)

    #
    #  Boundary Conditions
    #
    strh0 = search.find_tag_path(xml, ["amanzi_input","boundary_conditions","boundary_condition,LeftBC","liquid_phase","liquid_component","hydrostatic",]).get("value")
    params["h_0"] = float(strh0)
    strhL = search.find_tag_path(xml, ["amanzi_input","boundary_conditions","boundary_condition,RightBC","liquid_phase","liquid_component","hydrostatic",]).get("value")
    params["h_L"] = float(strhL)

    #
    #  Standard Gravity
    #
    params.setdefault("g",9.80665)

    # instantiate the class
    return LinearMaterialsParallel(params)

if __name__ == "__main__":

    # Instantiate the class
    lhh = LinearMaterialsParallel()

    # Get 11 equally spaced points: dx=(x_1-x_0)/10
    x = numpy.linspace(lhh.x_0,lhh.x_1,11)

    # Create space for a set of (x,z) points
    coords = numpy.zeros((11,2))
    # set x
    coords[:,0]=x
    # set z
    coords[:,1]=3

    # compute head.
    h1 = lhh.head(coords)

    # reset z
    coords[:,1]=7

    # compute head.
    h2 = lhh.head(coords)

    # plot
    plt.plot(x,h1)
    plt.plot(x,h2)
    plt.xlabel('x-coordinate [m]')
    plt.ylabel('Head [m]')

    # show the plot
    plt.show()

max_stars_repo_name: fmyuan/amanzi

max_issues_repo_name: fmyuan/amanzi

max_forks_repo_name: fmyuan/amanzi


text:
# -*- coding: utf-8 -*-
"""
@author: BELGADA Zainab
"""

import numpy as np
import pandas as pd
from flask import Flask, render_template, request
from datetime import datetime
from helper_functions import distance_de_Levenshtein, extraire_info

# Permet d'éviter un faux warning de pandas (false positive)
# source 1 : https://www.dataquest.io/blog/settingwithcopywarning/
# source 2 : https://stackoverflow.com/questions/20625582/how-to-deal-with-settingwithcopywarning-in-pandas
pd.set_option('mode.chained_assignment', None)

# Load dataset into a pandas dataframe
# data.csv est préparé en utilisant la fonction merge_datasets()
# définie dans le fichier preprocess.py
data = pd.read_csv('./data/data.csv', sep=',', header=0, low_memory=False)

app = Flask('laposte')

# Affiche la page d'accueil
@app.route('/')
def index():
    return render_template('index.html')

# Affiche les résultats de recherche
# Recherche par ville : On affiche soit la liste des bureaux de poste disponible
# dans la localité saisie par l'utilisateur OU bien une liste de suggestions de ville
# si l'utilisateur opte pour une recherche approchée
@app.route('/search')
def display_search_results():
    # extraire le nom de la ville recherché par l'utilisateur
    # et le transformé en majuscule (même format que le champ Localité dans le dataframe)
    nom_ville = request.args.get('nom_ville').upper()
    # Afficher un message personalisé si l'utilisateur saisi un nom de ville non valide
    # (chaîne de caractère vide. A améliorer pour tenir compte d'autre cas)
    if nom_ville.strip() == '':
        message = 'Veuillez saisir un nom de localité valide.'
        return render_template('results.html', message=message, bureaux_disponibles=None, nom_ville=nom_ville)
    # Si l'utilisateur souhaite faire une recherche par nom de ville incomplet
    # On affiche d'abord une liste des noms de ville qui contiennent la partie saisie
    # puis on classe les résultats obtenus par le score de Levenshtein
    elif request.args.get('resultats_approches') == 'true':
        # Liste des noms de ville (unique)
        villes_disponibles = pd.DataFrame({'Localité': data['Localité'].unique()})
        # Extraire les noms de ville qui contiennent le nom incomplet saisi par l'utilisateur
        suggestions_villes = villes_disponibles[villes_disponibles['Localité'].str.contains(nom_ville, case=False, regex=False)]
        suggestions_count = suggestions_villes.shape[0]
        # Message personalisé si on trouve pas de résultats similaires
        if suggestions_count == 0:
            message = 'Désolé ! Aucun résultat ne correspond à votre recherche.'
        # Sinon, classer les suggestion par score de Levenshtein
        else:
            # On définit une version "vectorisé" de la fonction distance_de_Levenshtein
            # afin de l'appliquer à chaque élément de l'array suggestions_villes['Localité']
            vect_distance_de_Levenshtein = np.vectorize(distance_de_Levenshtein)
            # Calculer Levenshtein Score
            suggestions_villes['Levenshtein_Score'] = vect_distance_de_Levenshtein(suggestions_villes['Localité'], nom_ville)
            suggestions_villes.sort_values('Levenshtein_Score', inplace=True)
            message = 'Nous avons trouvés {} {} à votre requête. '.format(suggestions_count, 'seule localité similaire' if suggestions_count == 1 else 'localités similaires')
            if suggestions_count == 1:
                message += 'La voici :'
            elif suggestions_count <= 10:
                message += 'Les voici classées par score de Levenshtein :'
            else:
                # Si y a plusieurs suggestions (> 10), afficher uniquement les 10 premières
                # (ayant une score faible, donc plus proches de la requête)
                suggestions_villes = suggestions_villes.head(n=10)
                message += 'Voici les top 10 suggestions classées par score de Levenshtein :'
        return render_template('levenshtein_results.html', message=message, nom_ville=nom_ville, suggestions_count=suggestions_count, suggestions_villes=suggestions_villes.to_dict('records'))
    else:
        # L'utilisateur ne coche pas l'option de recherche approcheé
        # donc souhaite de rechercher exactement le nom de ville saisi
        # Dans ce cas, on affiche directement la liste des bureaux de poste disponible dans cette ville
        ville_data = data[data['Localité'] == nom_ville]
        unique_bureaux = ville_data['Identifiant_du_site'].unique()
        count_results = unique_bureaux.shape[0]
        if count_results >= 1:
            message = 'Il y a ' + ('1 seul bureau' if count_results == 1 else '{} bureaux'.format(count_results))
            message += ' de poste à {}.'.format(nom_ville)
            bureaux_disponibles = extraire_info(data, unique_bureaux, ['Identifiant_du_site', 'Adresse'])
            return render_template('results.html', message=message, nom_ville=nom_ville, bureaux_disponibles=bureaux_disponibles)
        # Message personalisé si on trouve pas de résultats pour la localité saisie par l'utilisateur
        else:
            message = "Désolé ! Aucun bureau de poste n'a été trouvé dans cette localité."
            return render_template('results.html', message=message, nom_ville=nom_ville, bureaux_disponibles=None)

# Affiche des informations sur un bureau de poste spécifique
# ex. Localisation (avec lien Google Maps),
# Calendrier d'ouverture dans les prochains 15 jours...
@app.route('/details/<identifiant_du_site>')
def show_details(identifiant_du_site):
    # filtrer le dataset pour garder uniquement les lignes qui correspondent à
    # identifiant_du_site en question. NB: le dataset contient 3 mois de données par identifiant
    # (qui correspondent à l'horaire d'ouverture dans les prochains 3 mois)
    data_subset = data[data['Identifiant_du_site'] == identifiant_du_site]
    # Pour extraire des informations basiques (localisation, services disponibles...), on aura
    # besoin que de la première ligne.
    info = data_subset.iloc[0].to_dict()
    # Extraire les horaires d'ouvertures pour les prochains 15 jours
    # Pour ce faire, on aura besoin de la date d'aujourd'hui
    data_subset['Date_calendrier'] = pd.to_datetime(data_subset['Date_calendrier'], format='%Y-%m-%d')
    data_subset = data_subset.sort_values(['Date_calendrier'])
    current_date = datetime.today().strftime('%Y-%m-%d')
    horaires = data_subset[data_subset['Date_calendrier'] >= current_date].iloc[0:15]
    horaires['Date_calendrier'] = horaires['Date_calendrier'].dt.date
    horaires = horaires[['Date_calendrier', 'Plage_horaire_1', 'Plage_horaire_2',
                         'Heure_limite_dépôt_Courrier', 'Heure_limite_dépôt_Chrono', 'Heure_limite_dépôt_Colis']]
    # Afin d'afficher '-' au lieu de 'Nan' dans les horaires pour les jours de fermeture
    horaires.fillna('-', inplace=True)
    return render_template('details.html', info=info, days=horaires.to_dict('records'))

# Affiches les 100 premières villes du dataset classées par nombre de bureaux
# de poste disponibles
@app.route('/cities')
def show_cities():
    # nombre de villes disponibles dans le dataset
    number_of_cities = data['Localité'].unique().shape[0]
    # Calcule le nombre de bureaux de poste par ville
    # NB: on fait appel a nunique() au lieu de count() car chaque identifiant et représenté par
    # 3 mois de données (horaire d'ouverture dans le 3 prochains mois)
    top_n = data.groupby('Localité')['Identifiant_du_site'].nunique().sort_values(ascending=False)[:100].to_dict()
    # Classer les villes par le nombre de bureaux de poste
    top_n = sorted(list(top_n.items()), key=lambda tup: -tup[1])
    return render_template('cities.html', number_of_cities=number_of_cities, top_n=top_n)

if __name__ == '__main__':
    app.run(debug=True)

max_stars_repo_name: belzaina/postefrance

max_issues_repo_name: belzaina/postefrance

max_forks_repo_name: belzaina/postefrance


text:
# Copyright 2020 Adam Gleave
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#      http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

"""Test evaluating_rewards.tabular."""

from typing import Tuple

import hypothesis
from hypothesis import strategies as st
from hypothesis.extra import numpy as hp_numpy
import numpy as np
import pytest

from evaluating_rewards.distances import tabular

# pylint:disable=no-value-for-parameter
# pylint gets confused with hypothesis draw magic
@st.composite
def distribution(draw, shape) -> np.ndarray:
    """Search strategy for a probability distribution of given shape."""
    nonneg_elements = st.floats(min_value=0, max_value=1, allow_nan=False, allow_infinity=False)
    arr = draw(hp_numpy.arrays(np.float128, shape, elements=nonneg_elements, fill=st.nothing()))
    hypothesis.assume(np.any(arr > 0))
    return arr / np.sum(arr)

def numeric_float(max_abs: float = 1e3) -> st.SearchStrategy:
    """Search strategy for numeric (non-inf, non-NaN) floats with bounded absolute value."""
    return st.floats(min_value=-max_abs, max_value=max_abs, allow_nan=False, allow_infinity=False)

@st.composite
def arr_and_distribution(draw) -> Tuple[np.ndarray, np.ndarray]:
    """Search strategy for array and a probability distribution of same shape as array."""
    shape = draw(hp_numpy.array_shapes())
    arr = draw(
        hp_numpy.arrays(dtype=np.float128, shape=shape, elements=numeric_float(), fill=st.nothing())
    )
    dist = draw(distribution(shape))
    return arr, dist

@hypothesis.given(
    arr_dist=arr_and_distribution(), p=st.integers(min_value=1, max_value=10), scale=numeric_float()
)
def test_weighted_lp_norm(arr_dist: Tuple[np.ndarray, np.ndarray], p: int, scale: float) -> None:
    """Test for tabular.weighted_lp_norm."""
    arr, dist = arr_dist
    norm = tabular.lp_norm(arr, p, dist)
    # Non-negativity
    assert norm >= 0

    # Absolute homogeneity
    norm_scaled = tabular.lp_norm(scale * arr, p, dist)
    assert np.allclose(norm_scaled, abs(scale) * norm)

@st.composite
def reward(
    draw,
    n_states=st.integers(min_value=1, max_value=10),
    n_actions=st.integers(min_value=1, max_value=10),
) -> np.ndarray:
    """Search strategy for reward function with specified numbers of states and actions."""
    ns = draw(n_states)
    na = draw(n_actions)
    shape = (ns, na, ns)
    rew = draw(hp_numpy.arrays(np.float128, shape, elements=numeric_float(), fill=st.nothing()))
    return rew

_default_discount = st.floats(min_value=0.0, max_value=1.0)

@st.composite
def shaped_reward_pair(draw, discount=_default_discount) -> Tuple[np.ndarray, np.ndarray]:
    """Search strategy for a pair of rewards equivalent up to potential shaping."""
    rew = draw(reward())
    ns = rew.shape[0]
    potential = draw(
        hp_numpy.arrays(np.float128, (ns,), elements=numeric_float(), fill=st.nothing())
    )
    gamma = draw(discount)
    return rew, gamma, tabular.shape(rew, potential, gamma)

@st.composite
def equiv_reward_pair(draw, discount=_default_discount) -> Tuple[np.ndarray, np.ndarray]:
    """Search strategy for a pair of rewards equivalent up to shaping and positive rescaling."""
    base_rew, gamma, shaped_rew = draw(shaped_reward_pair(discount))
    scale = draw(st.floats(min_value=1e-2, max_value=1e2, exclude_min=True))
    equiv_rew = scale * shaped_rew
    return base_rew, gamma, equiv_rew

DESHAPE_FN = [
    tabular.singleton_shaping_canonical_reward,
    tabular.fully_connected_random_canonical_reward,
    tabular.fully_connected_greedy_canonical_reward,
]

@pytest.mark.parametrize("deshape_fn", DESHAPE_FN)
@hypothesis.given(base_shaped=shaped_reward_pair())
def test_deshape(deshape_fn: tabular.DeshapeFn, base_shaped) -> None:
    """Tests potential shaping invariance."""
    base_rew, discount, shaped_rew = base_shaped
    deshaped_base = deshape_fn(base_rew, discount)
    deshaped_shaped = deshape_fn(shaped_rew, discount)
    assert np.allclose(deshaped_base, deshaped_shaped, atol=1e-6)

def _check_scale(base_rew: np.ndarray, discount: float, deshape_fn: tabular.DeshapeFn) -> None:
    deshaped_base = deshape_fn(base_rew, discount)
    # Assume not near-zero rewards. This can cause problems with floating point errors
    # if things go slightly above/below round-to-zero thresholds.
    hypothesis.assume(np.linalg.norm(deshaped_base) > 1e-6)

@pytest.mark.parametrize("deshape_fn", DESHAPE_FN)
@hypothesis.given(base_equiv=equiv_reward_pair())
def test_canonical(deshape_fn: tabular.DeshapeFn, base_equiv) -> None:
    """Tests reward canonicalization."""
    base_rew, discount, equiv_rew = base_equiv
    _check_scale(base_rew, discount, deshape_fn)
    canon_base = tabular.canonical_reward(base_rew, discount, deshape_fn)
    canon_equiv = tabular.canonical_reward(equiv_rew, discount, deshape_fn)
    assert np.allclose(canon_base, canon_equiv, atol=1e-6)

@pytest.mark.parametrize("deshape_fn", DESHAPE_FN)
@hypothesis.given(base_equiv=equiv_reward_pair())
def test_canonical_dist(deshape_fn: tabular.DeshapeFn, base_equiv) -> None:
    """Test distance from canonicalized rewards."""
    base_rew, discount, equiv_rew = base_equiv
    _check_scale(base_rew, discount, deshape_fn)
    dist_equiv = tabular.canonical_reward_distance(base_rew, equiv_rew, discount, deshape_fn)
    assert np.allclose(dist_equiv, 0, atol=1e-6)

    dist_opposite = tabular.canonical_reward_distance(base_rew, -equiv_rew, discount, deshape_fn)
    # pylint:disable=comparison-with-callable
    if deshape_fn == tabular.fully_connected_random_canonical_reward:
        # Distance should be exactly equal to one for this function
        assert np.allclose(dist_opposite, 1.0, atol=1e-6)
    else:
        # For others it may be smaller since they take greedy policies, introducing an asymmetry
        # when rewards from a state differ depending on action. Usually it'll be >0.5, but there
        # are some pathological examples where it can be a lot lower.
        assert 0.01 < dist_opposite <= (1 + 1e-10)

@st.composite
def potential_only_reward(
    draw,
    n_states=st.integers(min_value=1, max_value=10),
    n_actions=st.integers(min_value=1, max_value=10),
    discount=_default_discount,
) -> np.ndarray:
    """Search strategy for a all-zero reward and a potential-shaping only reward."""
    ns = draw(n_states)
    na = draw(n_actions)
    rew = np.zeros((ns, na, ns))
    potential = draw(
        hp_numpy.arrays(np.float128, (ns,), elements=numeric_float(), fill=st.nothing())
    )
    gamma = draw(discount)
    shaped_rew = tabular.shape(rew, potential, gamma)
    return rew, gamma, shaped_rew

@pytest.mark.parametrize("deshape_fn", DESHAPE_FN)
@hypothesis.given(shaped_pair=potential_only_reward())
def test_canonical_dist_near_zero(deshape_fn: tabular.DeshapeFn, shaped_pair) -> None:
    """Test distance from zero-equivalent (i.e. potential-shaping only) rewards."""
    zero_rew, discount, shaped_rew = shaped_pair

    dist_equiv = tabular.canonical_reward_distance(zero_rew, shaped_rew, discount, deshape_fn)
    assert np.allclose(dist_equiv, 0, atol=1e-6)

    # Since it is equivalent to zero, it is also equivalent to the negative of itself!
    dist_opposite = tabular.canonical_reward_distance(zero_rew, -shaped_rew, discount, deshape_fn)
    assert np.allclose(dist_opposite, 0, atol=1e-6)

max_issues_repo_name: HumanCompatibleAI/evaluating-rewards