text:
# This code is part of Qiskit.
#
# (C) Copyright IBM 2020.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.
# pylint: disable=invalid-name

"""
Generator models module.
"""

from abc import ABC, abstractmethod
from typing import Callable, Tuple, Union, List, Optional
from copy import copy
import numpy as np

from qiskit import QiskitError
from qiskit.quantum_info.operators import Operator
from qiskit_dynamics.models.operator_collections import (
    DenseOperatorCollection,
    SparseOperatorCollection,
)
from qiskit_dynamics.dispatch import Array
from qiskit_dynamics.signals import Signal, SignalList
from qiskit_dynamics.type_utils import to_array, to_csr
from .rotating_frame import RotatingFrame

class BaseGeneratorModel(ABC):
    r"""Defines an interface for a time-dependent
    linear differential equation of the form :math:`\dot{y}(t) = \Lambda(t, y)`,
    where :math:`\Lambda` is linear in :math:`y`. The core functionality is
    evaluation of :math:`\Lambda(t, y)`, as well as,
    if possible, evaluation of the map :math:`\Lambda(t, \cdot)`.

    Additionally, the class defines interfaces for:

        - Setting a "drift" term, representing the time-independent part of :math:`\Lambda`.
        - Setting a "rotating frame", specified either directly as a :class:`RotatingFrame`
        instance, or an operator from which a :class:`RotatingFrame` instance can be constructed.
        The exact meaning of this transformation is determined by the structure of
        :math:`\Lambda(t, y)`, and is therefore by handled by concrete subclasses.
        - Setting an internal "evaluation mode", to set the specific numerical methods to use
        when evaluating :math:`\Lambda(t, y)`.
    """

    @property
    @abstractmethod
    def dim(self) -> int:
        """Gets matrix dimension."""
        pass

    @property
    @abstractmethod
    def rotating_frame(self) -> RotatingFrame:
        """Get the rotating frame."""
        pass

    @rotating_frame.setter
    @abstractmethod
    def rotating_frame(self, rotating_frame: RotatingFrame):
        """Set the rotating frame; either an already instantiated :class:`RotatingFrame` object,
        or a valid argument for the constructor of :class:`RotatingFrame`.
        """
        pass

    @property
    def evaluation_mode(self) -> str:
        """Numerical evaluation mode of the model."""
        # pylint: disable=no-member
        return self._evaluation_mode

    @evaluation_mode.setter
    def evaluation_mode(self, new_mode: str):
        """Sets evaluation mode of model."""

        # default error message to be used by all subclasses.
        raise NotImplementedError(
            "Evaluation mode '"
            + str(new_mode)
            + "' is not supported. "
            + "Call help("
            + str(self.__class__.__name__)
            + ".evaluation_mode) "
            + "for available options."
        )

    @abstractmethod
    def get_operators(
        self, in_frame_basis: Optional[bool] = False
    ) -> Union[Array, Tuple[Array], Callable]:
        """Get the operators used in the model construction.

        Args:
            in_frame_basis: Flag indicating whether to return the operators
            in the basis in which the rotating frame operator is diagonal.
        Returns:
            The operators in the basis specified by `in_frame_basis`.
        """
        pass

    def get_drift(self, in_frame_basis: Optional[bool] = False) -> Array:
        """Get the drift term.

        Args:
            in_frame_basis: Flag for whether the returned drift should be
            in the basis in which the frame is diagonal.
        Returns:
            The drift term.
        """
        if not in_frame_basis and self.rotating_frame is not None:
            return self.rotating_frame.operator_out_of_frame_basis(self._drift)
        else:
            return self._drift

    def set_drift(
        self,
        new_drift: Array,
        operator_in_frame_basis: Optional[bool] = False,
    ):
        """Sets drift term. Note that if the model has a rotating frame this will override
        any contributions to the drift due to the frame transformation.

        Args:
            new_drift: The drift operator.
            operator_in_frame_basis: Whether `new_drift` is already in the rotating
            frame basis.
        """
        if new_drift is None:
            new_drift = np.zeros((self.dim, self.dim))

        new_drift = to_array(new_drift)

        if not operator_in_frame_basis and self.rotating_frame is not None:
            new_drift = self.rotating_frame.operator_into_frame_basis(new_drift)
        # pylint: disable = attribute-defined-outside-init
        self._drift = new_drift
        # pylint: disable=no-member
        if self._operator_collection is not None:
            # pylint: disable=no-member
            self._operator_collection.drift = new_drift

    @abstractmethod
    def evaluate(self, time: float, in_frame_basis: Optional[bool] = False) -> Array:
        r"""If possible, evaluate the map :math:`\Lambda(t, \cdot)`.

        Args:
            time: Time.
            in_frame_basis: Whether to return the result in the rotating frame basis.
        """
        pass

    @abstractmethod
    def evaluate_rhs(self, time: float, y: Array, in_frame_basis: Optional[bool] = False) -> Array:
        r"""Evaluate the right hand side :math:`\dot{y}(t) = \Lambda(t, y)`.

        Args:
            time: Time.
            y: State of the differential equation.
            in_frame_basis: Whether `y` is in the rotating frame basis and the results should be
                returned in the rotatign frame basis.
        """
        pass

    def copy(self):
        """Return a copy of self."""
        return copy(self)

    def __call__(
        self, time: float, y: Optional[Array] = None, in_frame_basis: Optional[bool] = False
    ) -> Array:
        r"""Evaluate generator RHS functions. If ``y is None``,
        attemps to evaluate :math:`\Lambda(t, \cdot)`, otherwise, calculates
        :math:`\Lambda(t, y)`

        Args:
            time: Time.
            y: Optional state.
            in_frame_basis: Whether or not to evaluate in the rotating frame basis.

        Returns:
            Array: Either the evaluated model or the RHS for the given y
        """

        if y is None:
            return self.evaluate(time, in_frame_basis=in_frame_basis)

        return self.evaluate_rhs(time, y, in_frame_basis=in_frame_basis)

class CallableGenerator(BaseGeneratorModel):
    r"""Specifies a linear matrix differential equation of the form
    :math:`\dot{y}=\Lambda(t, y)=G(t)y` with :math:`G(t)` passed as a callable function.
    """

    def __init__(
        self,
        generator: Callable,
        rotating_frame: Optional[Union[Operator, Array, RotatingFrame]] = None,
        drift: Optional[Union[Operator, Array]] = None,
        dim: Optional[int] = None,
    ):
        self._generator = generator
        self.rotating_frame = rotating_frame
        self._drift = drift
        self._evaluation_mode = "callable_generator"
        self._operator_collection = None
        self._dim = dim

    @property
    def dim(self) -> int:
        if self._dim is not None:
            return self._dim
        else:
            raise ValueError(
                "Dimension of CallableGenerator object should be specified at initialization."
            )

    def get_drift(self, in_frame_basis: Optional[bool] = False) -> Array:
        if in_frame_basis and self.rotating_frame is not None:
            return self.rotating_frame.operator_into_frame_basis(self._drift)
        else:
            return self._drift

    def get_operators(self, in_frame_basis: Optional[bool] = False) -> Callable:
        """`CallableGenerator` does not have decomposition in terms of operators."""
        raise QiskitError("CallableGenerator does not have decomposition in terms of operators.")

    def set_drift(
        self,
        new_drift: Array,
        operator_in_frame_basis: Optional[bool] = False,
    ):
        # subtracting the frame operator from the generator is handled at evaluation time.
        if operator_in_frame_basis and self.rotating_frame is not None:
            self._drift = self.rotating_frame.operator_out_of_frame_basis(new_drift)
        else:
            self._drift = new_drift

        if self.rotating_frame.frame_diag is not None:
            self._drift = self._drift + self.rotating_frame.operator_out_of_frame_basis(
                self.rotating_frame.frame_diag
            )

    @property
    def evaluation_mode(self) -> str:
        """CallableGenerator does not support evaluation_mode."""
        pass

    @evaluation_mode.setter
    def evaluation_mode(self, new_mode: str):
        """Setting the evaluation mode for CallableGenerator
        is not supported."""
        raise NotImplementedError(
            "Setting implementation mode for CallableGenerator is not supported."
        )

    @property
    def rotating_frame(self) -> RotatingFrame:
        """Return the frame."""
        return self._rotating_frame

    @rotating_frame.setter
    def rotating_frame(self, rotating_frame: Union[Operator, Array, RotatingFrame]):
        """Set the frame; either an already instantiated :class:`RotatingFrame` object
        a valid argument for the constructor of :class:`RotatingFrame`, or `None`.
        """
        self._rotating_frame = RotatingFrame(rotating_frame)

    def evaluate_rhs(self, time: float, y: Array, in_frame_basis: Optional[bool] = False) -> Array:
        return self.evaluate(time, in_frame_basis=in_frame_basis) @ y

    def evaluate(self, time: float, in_frame_basis: Optional[bool] = False) -> Array:
        """Evaluate the model in array format.

        Args:
            time: Time to evaluate the model
            in_frame_basis: Whether to evaluate in the basis in which the frame
                            operator is diagonal

        Returns:
            Array: The evaluated model.
        """

        # evaluate generator and map it into the rotating frame
        gen = self._generator(time)
        if self._drift is not None:
            gen = gen + self._drift
        return self.rotating_frame.generator_into_frame(
            time, gen, operator_in_frame_basis=False, return_in_frame_basis=in_frame_basis
        )

class GeneratorModel(BaseGeneratorModel):
    r"""A model for a a linear matrix differential equation in standard form.

    :class:`GeneratorModel` is a concrete instance of :class:`BaseGeneratorModel`, where the
    map :math:`\Lambda(t, y)` is explicitly constructed as:

    .. math::
        \Lambda(t, y) = G(t)y,

        G(t) = \sum_i s_i(t) G_i + G_d

    where the :math:`G_i` are matrices (represented by :class:`Operator`
    or :class:`Array` objects), the :math:`s_i(t)` are signals represented by
    a list of :class:`Signal` objects, and
    :math:`G_d` is the constant-in-time drift term of the generator.
    """

    def __init__(
        self,
        operators: Array,
        signals: Optional[Union[SignalList, List[Signal]]] = None,
        drift: Optional[Array] = None,
        rotating_frame: Optional[Union[Operator, Array, RotatingFrame]] = None,
        evaluation_mode: str = "dense",
    ):
        """Initialize.

        Args:
            operators: A list of operators :math:`G_i`.
            signals: Stores the terms :math:`s_i(t)`. While required for evaluation,
                     :class:`GeneratorModel` signals are not required at instantiation.
            rotating_frame: Rotating frame operator.
            evaluation_mode: Evaluation mode to use. See ``GeneratorModel.evaluation_mode``
                             for more details. Supported options are:

                                - 'dense' (DenseOperatorCollection)
                                - 'sparse' (SparseOperatorCollection)

        """

        # initialize internal operator representation
        self._operator_collection = None
        self._operators = to_array(operators)
        self._drift = None
        self._evaluation_mode = None
        self.set_drift(drift, operator_in_frame_basis=True)

        # set frame and transform operators into frame basis.
        self._rotating_frame = None
        self.rotating_frame = RotatingFrame(rotating_frame)

        # initialize signal-related attributes
        self._signals = None
        self.signals = signals

        self.evaluation_mode = evaluation_mode

    def get_operators(self, in_frame_basis: Optional[bool] = False) -> Array:
        if not in_frame_basis and self.rotating_frame is not None:
            return self.rotating_frame.operator_out_of_frame_basis(self._operators)
        else:
            return self._operators

    @property
    def dim(self) -> int:
        return self._operators.shape[-1]

    @property
    def evaluation_mode(self) -> str:
        """Numerical evaluation mode of the model.

        Available options:

            - 'dense': Stores/evaluates operators using dense Arrays.
            - 'sparse': stores/evaluates operators using scipy
            :class:`csr_matrix` types. Not compatible with JAX.
        """
        return super().evaluation_mode

    @evaluation_mode.setter
    def evaluation_mode(self, new_mode: str):
        """Set evaluation mode.

        Args:
            new_mode: String specifying new mode. Available options
                      are 'dense' and 'sparse'. See property doc string for details.

        Raises:
            NotImplementedError: if new_mode is not one of the above
            supported evaluation modes.
        """

        if new_mode == "dense":
            self._operators = to_array(self._operators)
            self._operator_collection = DenseOperatorCollection(
                self.get_operators(in_frame_basis=True), drift=self.get_drift(in_frame_basis=True)
            )
        elif new_mode == "sparse":
            self._operators = to_csr(self._operators)
            self._operator_collection = SparseOperatorCollection(
                self.get_operators(in_frame_basis=True), drift=self.get_drift(in_frame_basis=True)
            )
        else:
            # raise error with standard message
            # this is equivalent to calling the setter of the base class
            super(__class__, self.__class__).evaluation_mode.fset(self, new_mode)

        self._evaluation_mode = new_mode

    @property
    def signals(self) -> SignalList:
        """Return the signals in the model."""
        return self._signals

    @signals.setter
    def signals(self, signals: Union[SignalList, List[Signal]]):
        """Set the signals."""

        if signals is None:
            self._signals = None
        else:
            # if signals is a list, instantiate a SignalList
            if isinstance(signals, list):
                signals = SignalList(signals)

            # if it isn't a SignalList by now, raise an error
            if not isinstance(signals, SignalList):
                raise QiskitError("Signals specified in unaccepted format.")

            # verify signal length is same as operators
            if len(signals) != self.get_operators().shape[0]:
                raise QiskitError(
                    """Signals needs to have the same length as
                                    operators."""
                )

            self._signals = signals

    @property
    def rotating_frame(self) -> RotatingFrame:
        """Return the rotating frame."""
        return self._rotating_frame

    @rotating_frame.setter
    def rotating_frame(self, rotating_frame: Union[Operator, Array, RotatingFrame]):
        if self._rotating_frame is not None and self._rotating_frame.frame_diag is not None:
            self._drift = self._drift + Array(np.diag(self._rotating_frame.frame_diag))
            self._operators = self.rotating_frame.operator_out_of_frame_basis(self._operators)
            self._drift = self.rotating_frame.operator_out_of_frame_basis(self._drift)

        self._rotating_frame = RotatingFrame(rotating_frame)

        if self._rotating_frame.frame_diag is not None:
            self._operators = self.rotating_frame.operator_into_frame_basis(self._operators)
            self._drift = self.rotating_frame.operator_into_frame_basis(self._drift)
            self._drift = self._drift - Array(np.diag(self._rotating_frame.frame_diag))

        # Reset internal operator collection
        if self.evaluation_mode is not None:
            self.evaluation_mode = self.evaluation_mode

    def evaluate(self, time: float, in_frame_basis: Optional[bool] = False) -> Array:
        """Evaluate the model in array format as a matrix, independent of state.

        Args:
            time: Time to evaluate the model.
            in_frame_basis: Whether to evaluate in the basis in which the rotating frame
            operator is diagonal.

        Returns:
            Array: the evaluated model as a (n,n) matrix

        Raises:
            QiskitError: If model cannot be evaluated.

        """

        if self._signals is None:
            raise QiskitError("GeneratorModel cannot be evaluated without signals.")

        # pylint: disable=not-callable
        sig_vals = self._signals(time)

        # Evaluated in frame basis, but without rotations
        op_combo = self._operator_collection(sig_vals)

        # Apply rotations e^{-Ft}Ae^{Ft} in frame basis where F = D
        return self.rotating_frame.operator_into_frame(
            time, op_combo, operator_in_frame_basis=True, return_in_frame_basis=in_frame_basis
        )

    def evaluate_rhs(self, time: float, y: Array, in_frame_basis: Optional[bool] = False) -> Array:
        r"""Evaluate `G(t) @ y`.

        Args:
            time: Time to evaluate the model.
            y: Array specifying system state, in basis specified by
            in_frame_basis.
            in_frame_basis: Whether to evaluate in the basis in which the frame
                            operator is diagonal.

        Returns:
            Array defined by :math:`G(t) \times y`.

        Raises:
            QiskitError: If model cannot be evaluated.

        """

        if self._signals is None:
            raise QiskitError("""GeneratorModel cannot be evaluated without signals.""")

        sig_vals = self._signals.__call__(time)

        # Evaluated in frame basis, but without rotations e^{\pm Ft}
        op_combo = self._operator_collection(sig_vals)

        if self.rotating_frame is not None:
            # First, compute e^{tF}y as a pre-rotation in the frame basis
            out = self.rotating_frame.state_out_of_frame(
                time, y, y_in_frame_basis=in_frame_basis, return_in_frame_basis=True
            )
            # Then, compute the product Ae^{tF}y
            out = op_combo @ out
            # Finally, we have the full operator e^{-tF}Ae^{tF}y
            out = self.rotating_frame.state_into_frame(
                time, out, y_in_frame_basis=True, return_in_frame_basis=in_frame_basis
            )
        else:
            return op_combo @ y

        return out

max_stars_repo_name: haggaila/qiskit-dynamics

max_issues_repo_name: haggaila/qiskit-dynamics

max_forks_repo_name: haggaila/qiskit-dynamics


text:
import cv2
import numpy as np
import os

from os.path import isfile, join

frame_array = []
# inPath ='/home/jacobphillips/sandbox/ac-teach/videos/03-04-2020_11-09-09/12626/'
# inPath ='/home/jacobphillips/sandbox/ac-teach/videos/03-04-2020_11-09-09/12632/'
# inPath ='/home/jacobphillips/sandbox/ac-teach/videos/03-04-2020_11-09-09/56559/'
inPath ='/home/jacobphillips/sandbox/ac-teach/videos/03-04-2020_11-09-09/56562/'
files = [f for f in os.listdir(inPath) if f[-3:]=='png']
files.sort(key = lambda x: int(x[:-4]))

for filename in files:
    img = cv2.imread(os.path.join(inPath, filename))
    frame_array.append(img)
size = frame_array[0].shape[:3]
out = cv2.VideoWriter("video56k_2.avi", 0, 10, (500,500))

for frame in frame_array:
    out.write(frame)
out.release()
print("done")

ext: py

reason: import numpy


text:
#!/usr/bin/python

#default module
import numpy as np
import matplotlib.pyplot as plt
import scipy.fftpack as fftp
import scipy.optimize as opt
from scipy import signal
import sys
import os
import string

from bcolors import *
from funct import *
from histogram import *

def periodogram(t,c,name):
    F=[]
    A=[]

    window = np.kaiser(c.shape[-1],5)
    c_win  = c * window

    f,a=DO_FFT(c_win, 1)
    F.append(f)
    A.append(a)

    plt.plot(F[0],A[0],'k', markersize=0.1)
    plt.show()

    #try:
    #    plt.savefig("graph/"+name+"_FFT.png", figuresize=(8,6), dpi=320, format="png")
    #    print_ok("Graph saved in: "+"graph/"+name+"_FFT.png")
    #except IOError as IoE:
    #    print_fail("I/O Error! Erro number = {0}; {1}".format(IoE.errno,IoE.strerror))
    #    exit(IoE.errno)

if __name__ == "__main__":

    t,c=load_file(sys.argv[1])

    try:
        periodogram(t,c,"corr")
        print ""
        print_ok("Stage1 completed - FFT")
        print ""
    except:
        print_fail("Some eorros has been occured in the Stage1")
        exit(1)

    exit(0)

max_stars_repo_name: gruffoc/AFO

max_issues_repo_name: gruffoc/AFO

max_forks_repo_name: gruffoc/AFO

reason: import numpy,import scipy,from scipy


text:
from sstcam_sandbox import get_astri_2019, get_plot
from CHECLabPy.core.io import TIOReader
from CHECLabPy.plotting.waveform import CameraPixelWaveformPlotter
from CHECLabPy.calib import PixelMasks, get_calib_data
from CHECOnsky.calib import get_calib_data as get_onsky_calib_data
from CHECLabPy.utils.mapping import get_superpixel_mapping
import numpy as np
import pandas as pd
from tqdm import tqdm
import matplotlib as mpl
from IPython import embed

class CameraPlotter(CameraPixelWaveformPlotter):
    def __init__(self, mapping):
        super().__init__(mapping)

    def save(self, output_path, **kwargs):
        super().save(output_path, **kwargs)
        for ax in self.ax_dict.values():
            ax.lines = []

class GetPlot:
    def __init__(self, mapping):
        self.mapping = mapping
        self.plots = dict()

    def __getitem__(self, item):
        if item not in self.plots:
            self.plots[item] = CameraPlotter(self.mapping)
        return self.plots[item]

def main():
    paths = dict(
        d190502="/Volumes/ICYBOX/astri_onsky_archive/d2019-06-10_ledflashers/flasher_comparison_May-June/unit0pattern-low_r1.tio",
        d190508="/Volumes/ICYBOX/astri_onsky_archive/d2019-06-10_ledflashers/flasher_comparison_May-June/Run13270_r1.tio",
        d190610="/Volumes/ICYBOX/astri_onsky_archive/d2019-06-10_ledflashers/flasher_comparison_May-June/Run13401_r1.tio",
    )
    readers = {
        key: TIOReader(path, max_events=100)
        for key, path in paths.items()
    }
    waveforms = {
        key: reader[:].mean(0) for key, reader in readers.items()
    }
    mapping = get_superpixel_mapping(readers['d190502'].mapping)

    pixelmasks = dict(
        d190428T2149=PixelMasks(get_calib_data("pixel_masks/d190428T2149.dat")),
    )

    # Calculate coeff
    waveforms_coeff = waveforms['d190502']
    n_pixels, n_samples = waveforms_coeff.shape
    n_sp = n_pixels // 4
    waveforms_sp = waveforms_coeff.reshape(
        (n_sp, 4, n_samples)
    ).sum(1)
    coeff = (waveforms_sp.max(1).mean() / waveforms_sp.max(1))[:, None]

    sp_waveforms = {
        key: wfs.reshape((n_sp, 4, n_samples)).sum(1)
        for key, wfs in waveforms.items()
    }
    sp_waveforms_cal = {
        key: wfs * coeff
        for key, wfs in sp_waveforms.items()
    }

    trigger_off = dict()
    for key, wfs in sp_waveforms_cal.items():
        amplitude = wfs.max(1)
        median = np.median(amplitude)
        threshold = median * 1.10
        trigger_off[key] = amplitude > threshold

    trigger_off['d190610'][[
        355,356,357,358,359,382,375,460,459,457,458,465,289,489,
        502,254,247,154,144,56,46,39,24,25,76
    ]] = True

    for key, sp in trigger_off.items():
        print(f"{key}: f{np.where(sp)}")

    plots = GetPlot(mapping)
    plots['d190428T2149'].highlight_pixels(trigger_off['d190502'])
    plots['d190428T2149'].plot_waveforms(sp_waveforms['d190502'])
    plots['d190428T2149'].save(get_plot("d190611_led_flasher/led/d190502.pdf"))
    plots['d190428T2149'].plot_waveforms(sp_waveforms_cal['d190502'])
    plots['d190428T2149'].save(get_plot("d190611_led_flasher/coeff/d190502.pdf"))
    plots['d190428T2149'].reset_pixel_highlight()

    plots['d190428T2149'].highlight_pixels(trigger_off['d190508'])
    plots['d190428T2149'].plot_waveforms(sp_waveforms['d190508'])
    plots['d190428T2149'].save(get_plot("d190611_led_flasher/led/d190508.pdf"))
    plots['d190428T2149'].plot_waveforms(sp_waveforms_cal['d190508'])
    plots['d190428T2149'].save(get_plot("d190611_led_flasher/coeff/d190508.pdf"))
    plots['d190428T2149'].reset_pixel_highlight()

    plots['d190428T2149'].highlight_pixels(trigger_off['d190610'])
    plots['d190428T2149'].plot_waveforms(sp_waveforms['d190610'])
    plots['d190428T2149'].save(get_plot("d190611_led_flasher/led/d190610.pdf"))
    plots['d190428T2149'].plot_waveforms(sp_waveforms_cal['d190610'])
    plots['d190428T2149'].save(get_plot("d190611_led_flasher/coeff/d190610.pdf"))
    plots['d190428T2149'].reset_pixel_highlight()

    hv_off = [24, 453]
    hv_off_array = np.zeros(n_sp, dtype=bool)
    hv_off_array[hv_off] = True
    trigger_off_array = np.zeros(n_sp, dtype=bool)
    trigger_off_array[trigger_off['d190610']] = True
    df = pd.DataFrame(dict(
        superpixel=np.arange(n_sp),
        hv_on=(~hv_off_array).astype(int),
        trigger_off=trigger_off_array.astype(int),
    ))
    df.to_csv(
        get_plot("d190611_led_flasher/hv_list.txt"),
        sep='\t', index=False
    )

if __name__ == '__main__':
    main()

max_stars_repo_name: watsonjj/CHECLabPySB

max_issues_repo_name: watsonjj/CHECLabPySB

max_forks_repo_name: watsonjj/CHECLabPySB


text:
""" \package MOSS With Locations
Plagiarism detector for source code files.
Detects plagiarism using Karp Rabin Hashing and winnowing
Marks overlapping hashes to exactly show the parts which where copied
"""

import numpy as np
from pygments.lexers import guess_lexer_for_filename
import pygments
import sys
import os

q=1000000007

boiler_fingerprint = []

def tokenize(filename):
	"""!
	\details Tokenises a file according to it's extension. Replaces all variable identifiers by V, function names by F and strings by S
		Other words are taken as it is
		Stores starting index of each token in original and final file

	\param filename: Name of file with correct extension
	\return mapping: a list of tokens of the form [token , start index in original file, start index in tokenised file]
	"""
	with open(filename, 'r') as f:
		contents = f.read()
		tokens = list(guess_lexer_for_filename(filename, contents).get_tokens(contents))
		mapping = []
		n = len(tokens)

		a = 0
		b = 0

		for i in range(n):
			if tokens[i][0] == pygments.token.Name.Variable:
				mapping.append(['V', a, b])
				b = b + 1

			elif tokens[i][0] == pygments.token.Name.Function:
				mapping.append(['F', a, b])
				b = b + 1

			elif tokens[i][0] == pygments.token.Literal.String:
				mapping.append(['S',a,b])
				b = b + 1

			elif tokens[i][0] != pygments.token.Text and tokens[i][0] != pygments.token.Comment:
				mapping.append([tokens[i][1], a, b])
				b = b + len(tokens[i][1])
			a = a + len(tokens[i][1])
		return mapping

def intersection(lst1, lst2):
	"""!

	\details Finds hashes that are common to both lists and stores their location in both documents
	Finds similarity that is measured by

	sim(A,B) = number of hashes in intersection of both hash sets divided by minimum of the number of hashes in lst1 and lst2

	\param  lst1 : 1st list whose elements are of the form [hash, start location, end location]
	\param lst2: 2nd list whose elements are of the form [hash, start location, end location]
	\return l3: list of common hashes and their locations in both documents. This is a list whose elements are of the form
		[common hash, [start location in 1, end location in 1], [start location in 2, end location in 2]]
	\return sim: similarity measure evaluated
	"""
	l1h = [h[0] for h in lst1]
	l2h = [h[0] for h in lst2]
	l1loc = {h[0]:h[1:] for h in lst1}
	l2loc = {h[0]:h[1:] for h in lst2}

	l3h = list(set(l1h)&set(l2h))
	l3 = [[h, l1loc[h], l2loc[h]] for h in l3h]
	sim = len(l3)/min(len(set(l1h)), len(set(l2h)))
	return l3, sim

def intersection_minus_boilerplate(lst1, lst2):
	"""!
	\details Finds hashes that are common to both lists and stores their locations in both documents.
	Also considers the fingerprint of boilerplate code and those hashes are not considered as part of the documents while evaluating similarity

	Finds similarity that is measured by

	sim(A,B) = number of hashes in intersection of both hash sets divided by minimum of the number of hashes in lst1 and lst2
	\param  lst1 : 1st list whose elements are of the form [hash, start location, end location]
	\param lst2: 2nd list whose elements are of the form [hash, start location, end location]
	\return l3: list of common hashes and their locations in both documents. This is a list whose elements are of the form
		[common hash, [start location in 1, end location in 1], [start location in 2, end location in 2]]
	\return sim: similarity measure evaluated
	"""

	l1h = [h[0] for h in lst1]
	l2h = [h[0] for h in lst2]
	l1loc = {h[0]:h[1:] for h in lst1}
	l2loc = {h[0]:h[1:] for h in lst2}

	boiler = set([h[0] for h in boiler_fingerprint])

	l3h = list(set(l1h)&set(l2h))
	l3h_without_boiler = [h for h in l3h if h not in boiler]

	l1h_without_boiler = set([h for h in l1h if h not in boiler])
	l2h_without_boiler = set([h for h in l2h if h not in boiler])

	sim = len(l3h)/min(len(l1h_without_boiler), len(l2h_without_boiler))
	l3 = [[h, l1loc[h], l2loc[h]] for h in l3h]
	return l3, sim

def similarity(lst1, lst2):
	"""!
		\details Evaluates similarity as done in intersection function but doesn't return locations of common hashes
	"""
	l1h = [h[0] for h in lst1]
	l2h = [h[0] for h in lst2]
	l3h = list(set(l1h)&set(l2h))
	sim = len(l3h)/min(len(set(l1h)), len(set(l2h)))
	return sim

def GetHLoc(t,k):
	"""!

	\details Reads the file and tokenises it. Evaluate its k-grams.
	For each k-gram, Karp-Rabin hash value is evaluated and stored in a list H. Along with the hashes, the location information given by start and end location of the k-gram in both initial and tokenized files are stored

	\param t: file name
	\param k: k-gram parameter

	\return H: The list H of hashes with their locations
	"""
	H = []
	mapping = tokenize(t)
	current_length = 0
	kgram = ""
	kgrams = []
	start_locs = []
	end_locs = []
	i = 0

	while i < len(mapping):
		if kgram == "":
			start_locs += [(mapping[i][1], mapping[i][2])]

		if current_length + len(mapping[i][0]) > k:
			deficit = k - current_length
			kgram = kgram + mapping[i][0][:deficit]
			kgrams += [kgram]
			end_locs += [[mapping[i][1], mapping[i][2]]]
			mapping[i][0] = mapping[i][0][deficit:]
			current_length = 0
			kgram = ""

		elif current_length + len(mapping[i][0]) == k:
			kgram = kgram + mapping[i][0]
			kgrams += [kgram]
			end_locs += [[mapping[i][1], mapping[i][2]]]
			current_length = 0
			kgram = ""
			i = i + 1

		else:
			kgram = kgram + mapping[i][0]
			current_length += len(mapping[i][0])
			i = i + 1
			if i == len(mapping):
				kgrams += [kgram]
				end_locs += [[mapping[i-1][1], mapping[i-1][2]]]

	for j in range(len(kgrams)):
		h = 0
		for k in kgrams[j]: h = (256*h + ord(k))%q
		H.append([h,start_locs[j], end_locs[j]])
	return H

def Win(H,t,k):
	"""!
        \details Applies Winnowing algorithm on given list of hashes with a window size such that for every substring of length t we will pick atleast one hash.
        Also stores the locations of selected hashes

        \param H: List of (hash, document ID, location)
        \param t: Winnowing threshold parameter
        \param k: k-gram parameter used while calculating hashes

        \return HS: Selected (winnowed) hashes
	"""
	HS=[]
	w=t+1-k
	n=len(H)
	mI=-1
	pmI=-1

        if(len(H) < w):
                if (len(H) > 0):
                    HS.append(H[0])
                else:
                    HS.append(q)

	for i in range(0,len(H)-w+1):
		tm=9223372036854775807
		for j in range(i, i+w):
			if H[j][0]<=tm:
				mI=j
				tm=H[j][0]
		if mI != pmI:
			pmI=mI
			HS.append(H[mI])
	return HS

def moss_all_pairs(files, boilerplate, t, k):
	"""!
        \details Evaluates MOSS similarity and matching portions between each pair of files
        \param files: list of files
        \param t: Winnowing threshold parameter
        \param k: k-gram parameter
        \return C: similarity matrix such that C[i][j] denotes the similarity between the ith and jth file
        \return markings: markings matrix such that markings[i][j] is a list of charavter indices where matching k-grams begin for ith and jth file"""

	n = len(files)
	boiler_fingerprint = Win(GetHLoc(boilerplate,k),t,k)
	H = [GetHLoc(f,k) for f in files]
	HS = [Win(h,t,k) for h in H]

	C = np.identity(n)
	markings = []
	for i in range(n):
		marks = []
		for j in range(n):
			marks += [[]]
		markings += [marks]

	for i in range(n):
		for j in range(i+1, n):
			intersect,sim = intersection_minus_boilerplate(HS[i], HS[j])
			C[i][j] = sim
			C[j][i] = sim
			for h in intersect:
				markings[i][j] += [(h[1][0][0],h[1][1][0])]
				markings[j][i] += [(h[2][0][0],h[2][1][0])]

	return C, markings

folder = sys.argv[1]
print(os.getcwd())
bo = sys.argv[2]
files = os.listdir(folder)
paths = [folder + '/' + f for f in files]
t = 10
k = 5
print(moss_all_pairs(paths,bo,t,k))


text:
import numpy as np
from random import randint
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors

class MinesweeperGame:
    def __init__(self, length, width, num_mines):
        self.length = length
        self.width = width
        self.num_mines = num_mines
        self.grid = np.zeros((width, length))
        self.generate_new_grid()

    def get_adjacent_count(self, x, y):
        if not self.in_bounds(x, y):
            # Error out of bounds
            pass
        count = 0
        for i in range(x-1, x+2):
            for j in range(y-1, y+2):
                if (self.in_bounds(i, j) and
                        self.grid[i][j] == 9
                        and not (i == x and j == y)):
                    count += 1
        return count

    def in_bounds(self, x, y):
        return not ((x < 0 or x >= self.width)
                    or (y < 0 or y >= self.length))

    def generate_new_grid(self):
        self.grid = np.zeros((self.width, self.length))
        # Place the mines
        remaining_mines = self.num_mines
        while remaining_mines > 0:
            i = randint(0, self.width-1)
            j = randint(0, self.length-1)
            if self.grid[i][j] != 9:
                self.grid[i][j] = 9
                remaining_mines -= 1
        # Set the clues
        for i in range(0, self.width):
            for j in range(0, self.length):
                if self.grid[i][j] != 9:
                    self.grid[i][j] = self.get_adjacent_count(i, j)

    def query_tile(self, x, y):
        return self.grid[x][y]

    def render_grid(self):
        # cmap = colors.ListedColormap(['white', 'red', 'red', 'red', 'blue'])
        arraycopy = np.copy(self.grid)
        arraycopy = np.rot90(arraycopy,1)
        colors0 = plt.cm.binary(np.linspace(0, 1, 1))
        colors1 = plt.cm.Reds(np.linspace(.2, .8, 8))
        # colors1 = np.flip(colors1,0)
        colors2 = plt.cm.binary(np.linspace(1, 1, 1))

        # combine them and build a new colormap
        colors = np.vstack((colors0, colors1, colors2))
        mymap = mcolors.LinearSegmentedColormap.from_list('my_colormap', colors)
        '''
        bounds = [0, .5, 1.5, 2.5, 8, 10]
        norm = colors.BoundaryNorm(bounds, cmap.N)
        fig, ax = plt.subplots()
        ax.imshow(self.grid, cmap=cmap, norm=norm)
        '''
        plt.pcolor(arraycopy, cmap=mymap)
        plt.colorbar()
        plt.show()
        '''
        # draw gridlines
        ax.grid(which='major', axis='both',
                linestyle='-', color='k', linewidth=2)
        ax.set_xticks(np.arange(-.5, 10, 1))
        ax.set_yticks(np.arange(-.5, 10, 1))
        plt.show()
        '''

if __name__ == '__main__':
    ms_game = MinesweeperGame(10, 10, 25)
    ms_game.render_grid()

max_stars_repo_name: chris-gong/mine-sweeper

max_issues_repo_name: chris-gong/mine-sweeper

max_forks_repo_name: chris-gong/mine-sweeper


text:
from email.mime.text import MIMEText
from gate import Gate
import os
import os.path
from os import path
import time
import random
import signal
import numpy as np
import smtplib, ssl
import getpass
from socket import gaierror

class GateLab(Gate):

	exit_message = "Good bye !"

	def __init__(self, args):
		self.notdoneyet = True
		Gate.__init__(self, args)
		if self.go:
			self.handleExecutions()
		else:
			# Here we precise to execute the Gate version not to remove the lock of the already running process
			msg = "The file " + self.jobfile + " is already in process"
			self.exitMessage(msg)

	def getFakeList(self):
		fakeList = [{'identifier': 1, 'status': "Running"},
					{'identifier': "workflow-EPrOSZ", 'status': "Finished"},
					{'identifier': "workflow-DaNIel", 'status': "Held"},
					{'identifier': 4, 'status': "Killed"},
					{'identifier': 5, 'status': "Killed"},
					{'identifier': 6, 'status': "Killed"},
					]
		return fakeList

	def handleExecutions(self):
		n_free_slots = 0
		# Is it necessary to checkJobs here as we won't be able to launch several time the same jobfile ?
		# => workflowID won't be identical in two different jobfile
		self.tryNewSubmit()
		# if notdoneyet:
		while self.notdoneyet:
			result = input("Jobs are launched, do you want more to be launched ? yes / no\n")
			if result == "yes":
				# check if there is a free slot on VIP
				self.tryNewSubmit()
				# check if jobs are finished
				self.checkFinishedJobs()
				# check if jobs are held
				self.checkHeldJobs()
				# wait before new check
				if self.notdoneyet:
					time.sleep(5)
			else:
				break
		else:
			self.exitMessage("No more jobs to launch. Good bye !")

	def tryNewSubmit(self):
		n_free_slots = self.getFreeSlots()
		if n_free_slots:
			self.submitJobs(n_free_slots)

	# Retrieve running jobs and return available slots from vip list executions
	def getFreeSlots(self):
		n_jobs = 0

		# get list of jobs on vip
		if os.environ['DEBUG_VIP'] != "1":
			# retrieve workflowIDs
			execList = vip.list_executions()
		else:
			execList = self.getFakeList()

		for anExec in execList:
			if anExec['status'] == "Running":
				n_jobs += 1

		return self.maxExecsNb - n_jobs

	def submitJobs(self, n_free_slots):
		for i in range(n_free_slots):
			job = self.getNextJob()
			if job == False:
				print("No more job to launch, it's over")
				# save self.joblist to file before exiting
				self.saveJobList()
				# breaks jobs' execution's
				self.notdoneyet = false
			else:
				print("Starting job: ", job.workflowID)
				workflowID = self.launchExecution(job)
				#TODO : add a test if the launch has failed => do not change the job status
				# set the job status to submitted with a timestamp and set its workflowID
				self.setJobSubmitted(job, workflowID)
		# save self.joblist to file before exiting
		self.saveJobList()

	def checkFinishedJobs(self):
		# Just a short name
		l = self.joblist
		changed = False
		# get list of jobs on vip
		if os.environ['DEBUG_VIP'] != "1":
			# retrieve workflowIDs
			execList = vip.list_executions()
		else:
			execList = self.getFakeList()

		for anExec in execList:
			if anExec['status'] == "Finished" and (l.loc[l['workflowID']==anExec['identifier'],["finished"]]).values[0][0] == 0:
				print("Job id ", anExec['identifier'], " finished")
				changed = True
				self.setJobFinished(anExec['identifier'])
		if changed:
			self.saveJobList()

	def checkHeldJobs(self):
		# Just a short name
		l = self.joblist
		changed = False

		if os.environ['DEBUG_VIP'] != "1":
			# retrieve workflowIDs
			execList = vip.list_executions()
		else:
			execList = self.getFakeList()

		for anExec in execList:
			# Here we define that if the workflow of the exec_list is Held
			# we define the finished value of the joblist to 2
			# this imply that we have delt with this case (email sent)
			if anExec['status'] == "Held" and (l.loc[l['workflowID']==anExec['identifier'],["finished"]]).values[0][0] != 2:
				print("Job id ", anExec['identifier'], " held")
				changed = True
				self.setJobHeld(anExec['identifier'])
				self.sendMail(anExec['identifier'])
		if changed:
			self.saveJobList()

	def getNextJob(self):
		# return the first job in the list with status not submitted (0)
		for job in self.joblist.itertuples():
			if job.submitted == "0":
				return job
		return False

	def saveJobList(self):
		# save the self.joblist to a file
		self.joblist.to_csv(self.jobfile, index=None)

	def launchExecution(self, job):
		# job looks like this:
		# Pandas(Index=0, model='AF', source=61, particle='gamma', energy=0.2, primaries=100000000, seed=2614427, cpuParam=2, workflowID='workflow-mBt3pB', submitted=0, finished=0, downloaded=0)
		# you can access model with simply job.model
		gaterelease = self.config["gaterelease"]
		application = self.config["application"]
		CPUparam = self.config["CPUparam"]
		gateinput = self.config["gateinput"]
		macfile = self.config["macfile"]
		outputdir = self.config["outputdir"]
		# build the -a [...] string passed to Gate via Gatelab through parameter "phaseSpace"
		alias_string= "-a [Source_ID," + str(job.source) + "][particle," + str(job.particle) + "][energy," + str(job.energy) + "][nb," + str(job.primaries) + "][seed," + str(job.seed) + "]"

		executionName = "OpenDose_GateLab" + "_" + job.model + "_" + str(job.source) + "_" + str(job.particle) + "_" + str(job.energy) + "_" + str(job.primaries)
		datetime = time.strftime('%d-%m-%Y %H:%M')
		print (datetime + " - Launching execution " + executionName)

		if os.environ['DEBUG_VIP'] != "1":
			execID = vip.init_exec(application, executionName, {'CPUestimation':CPUparam, 'ParallelizationType':"stat", 'GateRelease':gaterelease, 'NumberOfParticles':job.particle, 'GateInput':gateinput, 'phaseSpace':alias_string})
		else:
			execID = "workflow-" + str(random.randrange(1000))
		#TODO : simulate a fake execution in a thread in a random time
		print ("execution id : " + execID)
		return execID

	def setJobSubmitted(self, job, workflowID):
		# set the job status to submitted with a timestamp and set its workflowID
		datetime = time.strftime('%d-%m-%Y %H:%M')
		self.joblist.loc[job.Index, ['workflowID']] = workflowID
		self.joblist.loc[job.Index, ['submitted']] = datetime

	def setJobFinished(self, workflowID):
		# set the job status to finished with a timestamp
		datetime = time.strftime('%d-%m-%Y %H:%M')
		self.joblist.loc[self.joblist['workflowID']==workflowID,["finished"]] = datetime

	def setJobHeld(self, workflowID):
		# set the job status to held
		self.joblist.loc[self.joblist['workflowID']==workflowID,["finished"]] = 2

	# Email sending via SSL (http://zetcode.com/python/smtplib/)
	def sendMail(self, workflowID):
		print("Email send to ", self.config['email'])
		smtp_server = "smtp.inserm.fr"
		port = 587

		password = getpass.getpass("Type your password and press enter : ")
		receiver = self.config['email']
		sender = self.config['email']
		user = self.config['email']

		msg = 	''.join(("Hi,\nYour workflow ",
						 workflowID,
						 " is in held status.\n",
						 "Regards"))

		msg = MIMEText(msg)

		msg['Subject'] = '[opendose] Held workflow'
		msg['From'] = 'opendose@gatelab.com'
		msg['To'] = receiver

		# Secure SSL context
		# Try to log in o server and send email
		try :
			with smtplib.SMTP("smtp.inserm.fr", port) as server:
			    server.starttls() # Secure the connection
			    server.login(user, password)
			    server.sendmail(sender, receiver, msg.as_string())
		except Exception as e:
			print(e)
		finally:
			print("Mail successfully sent to ", receiver)

	def handler(self, signum, frame):
		print(" Interruption signal catched")
		global exit_message
		self.notdoneyet = False
		exit_message = "Safe exit"

	def atExit(self):
		if path.exists(self.lock_file) and self.go:
			os.remove(self.lock_file)

max_stars_repo_name: CISI-INSERM/OpenDose-Computing

max_issues_repo_name: CISI-INSERM/OpenDose-Computing

max_forks_repo_name: CISI-INSERM/OpenDose-Computing


text:
import numpy as np
from scipy.constants import pi

from ... import njit

@njit
def steffen_3d(
    v_in,
    z_in,
    z_out,
    z_min_surface,
    z_max_surface,
    lower_extrapolation_with_gradient=False,
    upper_extrapolation_with_gradient=False,
):
    """
    Performs Steffen interpolation on each individual column.
    Steffen, M. (1990). A simple method for monotonic interpolation
    in one dimension. Astronomy and Astrophysics, 239, 443.

    `z_min_surface` and `z_max_surface` define the minimum and maximum height
    outside of which the interpolated value will be set to nan (useful for not
    interpolating pressure into solid ground for example).
    `lower_extrapolation_with_gradient` and `upper_extrapolation_with_gradient`
    enables the use of one-sided gradient when extrapolating outside the range
    of `z_in`, otherwise the limit value of `v_in` is used by default
    """
    assert v_in.shape[1:] == z_min_surface.shape
    k_max = v_in.shape[0]
    j_max = v_in.shape[1]
    i_max = v_in.shape[2]
    k_max_output = z_out.shape[0]
    k_max_minus = k_max - 1
    linear_slope = np.empty((k_max))
    v_out = np.empty((k_max_output, j_max, i_max))
    for i in range(i_max):
        for j in range(j_max):
            # first point
            delta_lower = z_in[1, j, i] - z_in[0, j, i]
            delta_upper = z_in[2, j, i] - z_in[1, j, i]
            if delta_lower < 0:
                raise Exception("Non-montonic increase in z_in")
            if delta_upper < 0:
                raise Exception("Non-montonic increase in z_in")
            slope_lower = (v_in[1, j, i] - v_in[0, j, i]) / delta_lower
            slope_upper = (v_in[2, j, i] - v_in[1, j, i]) / delta_upper

            weighted_slope = slope_lower * (
                1 + delta_lower / (delta_lower + delta_upper)
            ) - slope_upper * delta_lower / (delta_lower + delta_upper)
            if weighted_slope * slope_lower <= 0.0:
                linear_slope[0] = 0.0
            elif np.abs(weighted_slope) > 2 * np.abs(slope_lower):
                linear_slope[0] = 2.0 * slope_lower
            else:
                linear_slope[0] = weighted_slope

            # intermediate points
            for k in range(1, k_max_minus):
                delta_lower = z_in[k, j, i] - z_in[k - 1, j, i]
                delta_upper = z_in[k + 1, j, i] - z_in[k, j, i]
                slope_lower = (v_in[k, j, i] - v_in[k - 1, j, i]) / delta_lower
                slope_upper = (v_in[k + 1, j, i] - v_in[k, j, i]) / delta_upper
                weighted_slope = (
                    slope_lower * delta_upper + slope_upper * delta_lower
                ) / (delta_lower + delta_upper)

                if slope_lower * slope_upper <= 0.0:
                    linear_slope[k] = 0.0
                elif np.abs(weighted_slope) > 2.0 * np.abs(slope_lower):
                    linear_slope[k] = np.copysign(2.0, slope_lower) * min(
                        np.abs(slope_lower), np.abs(slope_upper)
                    )
                elif np.abs(weighted_slope) > 2.0 * np.abs(slope_upper):
                    linear_slope[k] = np.copysign(2.0, slope_lower) * min(
                        np.abs(slope_lower), np.abs(slope_upper)
                    )
                else:
                    linear_slope[k] = weighted_slope

            # last point
            delta_lower = z_in[k_max_minus - 1, j, i] - z_in[k_max_minus - 2, j, i]
            delta_upper = z_in[k_max_minus, j, i] - z_in[k_max_minus - 1, j, i]
            slope_lower = (
                v_in[k_max_minus - 1, j, i] - v_in[k_max_minus - 2, j, i]
            ) / delta_lower
            slope_upper = (
                v_in[k_max_minus, j, i] - v_in[k_max_minus - 1, j, i]
            ) / delta_upper
            weighted_slope = slope_upper * (
                1 + delta_upper / (delta_upper + delta_lower)
            ) - slope_lower * delta_upper / (delta_upper + delta_lower)
            if weighted_slope * slope_upper <= 0.0:
                linear_slope[k_max_minus] = 0.0
            elif np.abs(weighted_slope) > 2.0 * np.abs(slope_upper):
                linear_slope[k_max_minus] = 2.0 * slope_upper
            else:
                linear_slope[k_max_minus] = weighted_slope

            # loop over output points
            k_temp = 0
            for k_out in range(k_max_output):
                while (k_temp < k_max) and (z_in[k_temp, j, i] < z_out[k_out]):
                    k_temp = k_temp + 1
                if 0 < k_temp < k_max:
                    k_high = k_temp
                    k_low = k_high - 1
                    delta = z_in[k_high, j, i] - z_in[k_low, j, i]
                    slope = (v_in[k_high, j, i] - v_in[k_low, j, i]) / delta
                    a = (linear_slope[k_low] + linear_slope[k_high] - 2 * slope) / (
                        delta * delta
                    )
                    b = (
                        3 * slope - 2 * linear_slope[k_low] - linear_slope[k_high]
                    ) / delta
                    c = linear_slope[k_low]
                    d = v_in[k_low, j, i]
                    t_1 = z_out[k_out] - z_in[k_low, j, i]
                    t_2 = t_1 * t_1
                    t_3 = t_2 * t_1
                    v_out[k_out, j, i] = a * t_3 + b * t_2 + c * t_1 + d
                elif (k_temp == 0) and (z_out[k_out] >= z_min_surface[j, i]):
                    if lower_extrapolation_with_gradient:
                        v_out[k_out, j, i] = v_in[0, j, i] + linear_slope[0] * (
                            z_out[k_out] - z_in[0, j, i]
                        )
                    else:
                        v_out[k_out, j, i] = v_in[0, j, i]
                elif (k_temp == k_max) and (z_out[k_out] <= z_max_surface[j, i]):
                    if upper_extrapolation_with_gradient:
                        v_out[k_out, j, i] = v_in[k_max - 1, j, i] + linear_slope[
                            k_max - 1
                        ] * (z_out[k_out] - z_in[k_max - 1, j, i])
                    else:
                        v_out[k_out, j, i] = v_in[k_max - 1, j, i]
                else:
                    v_out[k_out, j, i] = np.nan
    return v_out

# This can probably be replaced by a generic Steffen interpolation function
# No extrapolation is performed
# Time axis present, but not lat, lon
@njit
def steffen_1d_no_ep_time(
    input_data,
    input_levels,
    output_level_array,
):
    """Performs Steffen interpolation on one individual column for
    each time step.
    Steffen, M. (1990). A simple method for monotonic interpolation in
    one dimension. Astronomy and Astrophysics, 239, 443."""
    t_max = input_data.shape[0]
    k_max = input_data.shape[1]
    k_max_output = output_level_array.shape[0]
    k_max_minus = k_max - 1
    linear_slope = np.empty((k_max))
    output_data = np.empty((t_max, k_max_output))
    # first point
    delta_lower = input_levels[1] - input_levels[0]
    delta_upper = input_levels[2] - input_levels[1]
    if delta_lower < 0:
        raise Exception("Non-montonic increase in input_levels")
    if delta_upper < 0:
        raise Exception("Non-montonic increase in input_levels")
    for time_index in range(t_max):
        slope_lower = (
            input_data[time_index, 1] - input_data[time_index, 0]
        ) / delta_lower
        slope_upper = (
            input_data[time_index, 2] - input_data[time_index, 1]
        ) / delta_upper
        weighted_slope = slope_lower * (
            1 + delta_lower / (delta_lower + delta_upper)
        ) - slope_upper * delta_lower / (delta_lower + delta_upper)
        if weighted_slope * slope_lower <= 0.0:
            linear_slope[0] = 0.0
        elif np.abs(weighted_slope) > 2 * np.abs(slope_lower):
            linear_slope[0] = 2.0 * slope_lower
        else:
            linear_slope[0] = weighted_slope

        # intermediate points
        for k in range(1, k_max_minus):
            delta_lower = input_levels[k] - input_levels[k - 1]
            delta_upper = input_levels[k + 1] - input_levels[k]
            slope_lower = (
                input_data[time_index, k] - input_data[time_index, k - 1]
            ) / delta_lower
            slope_upper = (
                input_data[time_index, k + 1] - input_data[time_index, k]
            ) / delta_upper
            weighted_slope = (slope_lower * delta_upper + slope_upper * delta_lower) / (
                delta_lower + delta_upper
            )

            if slope_lower * slope_upper <= 0.0:
                linear_slope[k] = 0.0
            elif np.abs(weighted_slope) > 2.0 * np.abs(slope_lower):
                linear_slope[k] = np.copysign(2.0, slope_lower) * min(
                    np.abs(slope_lower), np.abs(slope_upper)
                )
            elif np.abs(weighted_slope) > 2.0 * np.abs(slope_upper):
                linear_slope[k] = np.copysign(2.0, slope_lower) * min(
                    np.abs(slope_lower), np.abs(slope_upper)
                )
            else:
                linear_slope[k] = weighted_slope

        # last point
        delta_lower = input_levels[k_max_minus - 1] - input_levels[k_max_minus - 2]
        delta_upper = input_levels[k_max_minus] - input_levels[k_max_minus - 1]
        slope_lower = (
            input_data[time_index, k_max_minus - 1]
            - input_data[time_index, k_max_minus - 2]
        ) / delta_lower
        slope_upper = (
            input_data[time_index, k_max_minus]
            - input_data[time_index, k_max_minus - 1]
        ) / delta_upper
        weighted_slope = slope_upper * (
            1 + delta_upper / (delta_upper + delta_lower)
        ) - slope_lower * delta_upper / (delta_upper + delta_lower)
        if weighted_slope * slope_upper <= 0.0:
            linear_slope[k_max_minus] = 0.0
        elif np.abs(weighted_slope) > 2.0 * np.abs(slope_upper):
            linear_slope[k_max_minus] = 2.0 * slope_upper
        else:
            linear_slope[k_max_minus] = weighted_slope

        # loop over output points
        k_temp = 0
        for k_out in range(k_max_output):
            while (k_temp < k_max) and (
                input_levels[k_temp] < output_level_array[k_out]
            ):
                k_temp = k_temp + 1
            if 0 < k_temp < k_max:
                k_high = k_temp
                k_low = k_high - 1
                delta = input_levels[k_high] - input_levels[k_low]
                slope = (
                    input_data[time_index, k_high] - input_data[time_index, k_low]
                ) / delta
                a = (linear_slope[k_low] + linear_slope[k_high] - 2 * slope) / (
                    delta * delta
                )
                b = (3 * slope - 2 * linear_slope[k_low] - linear_slope[k_high]) / delta
                c = linear_slope[k_low]
                d = input_data[time_index, k_low]
                t_1 = output_level_array[k_out] - input_levels[k_low]
                t_2 = t_1 * t_1
                t_3 = t_2 * t_1
                output_data[time_index, k_out] = a * t_3 + b * t_2 + c * t_1 + d
            # Allow for small deviations in upper/lower levels
            elif (k_temp == 0) and (
                abs(output_level_array[k_out] - input_levels[k_temp]) < 1e-6
            ):
                output_data[time_index, k_out] = input_data[time_index, 0]
            elif (k_temp == k_max) and (
                abs(output_level_array[k_out] - input_levels[k_temp]) < 1e-6
            ):
                output_data[time_index, k_out] = input_data[time_index, k_max]
            else:
                output_data[time_index, k_out] = np.nan
    return output_data

def central_estimate(a_in):
    """take a one-sided difference at the edges, and a central difference elsewhere"""
    return np.concatenate(([a_in[0]], 0.5 * (a_in[1:-1] + a_in[2:]), [a_in[-1]]))

def cos_transition(absolute_input, transition_start, transition_end):
    """function that smoothly transitions from 1 to 0 using a
    cosine-shaped transition between start and end"""
    normalised_input = (absolute_input - transition_start) / (
        transition_end - transition_start
    )
    weight_factor = 1.0 * (normalised_input < 0.0) + (
        0.5 + 0.5 * np.cos(normalised_input * pi)
    ) * (1.0 - (normalised_input < 0.0) - (normalised_input > 1.0))
    return weight_factor

max_stars_repo_name: xychen-ocn/lagtraj

max_issues_repo_name: xychen-ocn/lagtraj

max_forks_repo_name: xychen-ocn/lagtraj

reason: import numpy,from scipy


text:
import argparse
import datetime
import pickle
from subprocess import call
import os
import requests
import io

from flymazerl.agents.classical import *
from flymazerl.agents.phenomenological import *
from flymazerl.gym.environment import *

import numpy as np
import pandas as pd

start_text = """
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
███████╗██╗     ██╗   ██╗███╗   ███╗ █████╗ ███████╗███████╗██████╗ ██╗
██╔════╝██║     ╚██╗ ██╔╝████╗ ████║██╔══██╗╚══███╔╝██╔════╝██╔══██╗██║
█████╗  ██║      ╚████╔╝ ██╔████╔██║███████║  ███╔╝ █████╗  ██████╔╝██║
██╔══╝  ██║       ╚██╔╝  ██║╚██╔╝██║██╔══██║ ███╔╝  ██╔══╝  ██╔══██╗██║
██║     ███████╗   ██║   ██║ ╚═╝ ██║██║  ██║███████╗███████╗██║  ██║███████╗
╚═╝     ╚══════╝   ╚═╝   ╚═╝     ╚═╝╚═╝  ╚═╝╚══════╝╚══════╝╚═╝  ╚═╝╚══════╝
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Developed by:
    Rishika Mohanta, Research Technician, Turner Lab, Janelia Research Campus
"""

print(start_text)

model_database = pd.read_csv("https://raw.githubusercontent.com/neurorishika/flymazerl/main/model_description.csv")

argument_parser = argparse.ArgumentParser(
    description=start_text
    + "Script to fit the parameters of the a variety of Classical RL agents to the a 2AFC maze task using Bayesian MCMC."
)
argument_parser.add_argument("--save_path", type=str, default="../fits/", help="Path to save the fits.")
argument_parser.add_argument(
    "--agent",
    type=str,
    default="CQES",
    help="The agent to fit the dynamics to (Ref: SHORTCODE column in model_description.csv)",
)
argument_parser.add_argument(
    "--action-set-data",
    type=str,
    default=".././data/action_set.csv",
    help="The file to load the action set data from.",
)
argument_parser.add_argument(
    "--reward-set-data",
    type=str,
    default=".././data/reward_set.csv",
    help="The file to load the reward set data from.",
)
argument_parser.add_argument("--n_chains", type=int, default=1, help="Number of MC-MC chains to sample")
argument_parser.add_argument("--n_burn", type=int, default=1000, help="Number of MC-MC burn-in steps per chain")
argument_parser.add_argument("--n_steps", type=int, default=2500, help="Number of MC-MC sampling steps per chain")
argument_parser.add_argument(
    "--target_acceptance_rate", type=float, default=0.8, help="Target acceptance rate for MC-MC"
)
argument_parser.add_argument(
    "--n_cores", type=int, default=1, help="Number of cores to use for parallelization",
)
argument_parser.add_argument("--n_splits", type=int, default=1, help="Number of subdatasets to split the data into")
args = argument_parser.parse_args()

assert args.agent in model_database.SHORTCODE.values, "Agent not found in model_description.csv"
assert args.save_path[-1] == "/", "Save path must end with a '/'"
assert os.path.exists(args.save_path), "Save path does not exist"

agentClass = eval(model_database.loc[model_database.SHORTCODE == args.agent, "AgentClass"].values[0])

timestring = datetime.datetime.now().strftime("%Y_%m_%d_%H_%M_%S")

print("Fitting agent: {}".format(agentClass.__name__))
print("Loading data from:\n{}\n{}".format(args.action_set_data, args.reward_set_data))

if args.n_splits == 1:
    if args.action_set_data.startswith("http"):
        # use requests to download the data
        response = requests.get(args.action_set_data)
        response.raise_for_status()
        action_set = np.loadtxt(io.BytesIO(response.content), delimiter=",", dtype=np.int32)
    else:
        action_set = np.loadtxt(args.action_set_data, delimiter=",", dtype=np.int32)
    if args.reward_set_data.startswith("http"):
        # use requests to download the data
        response = requests.get(args.reward_set_data)
        response.raise_for_status()
        reward_set = np.loadtxt(io.BytesIO(response.content), delimiter=",", dtype=np.int32)
    else:
        reward_set = np.loadtxt(args.reward_set_data, delimiter=",", dtype=np.int32)

    print("Loaded action and reward sets.")

    assert len(action_set) == len(reward_set), "Action and reward set must have the same length."
    n_sessions = len(action_set)
    n_trials_in_each_session = set([len(action_set[i]) for i in range(n_sessions)])

    assert len(n_trials_in_each_session) == 1, "All sessions must have the same number of trials."
    n_trials_per_session = n_trials_in_each_session.pop()

    print("Found {} sessions with {} trials each.".format(n_sessions, n_trials_per_session))

    fitAgent = agentClass(ymaze_static(n_trials_per_session))

    extra_code = model_database.loc[model_database.SHORTCODE == args.agent, "FittingExtras"].values[0]
    if extra_code != "None":
        print("Found extra code for agent initialization.")
        exec(extra_code)

    print("Agent successfully initialized. Starting fitting.")

    if model_database.loc[model_database.SHORTCODE == args.agent, "OptimizationAlgorithm"].values[0] == "NUTS":
        model, data = fitAgent.fit(
            action_set,
            reward_set,
            niters=args.n_steps,
            ntune=args.n_burn,
            nchains=args.n_chains,
            nparallel=args.n_cores,
            target_accept=args.target_acceptance_rate,
        )
    elif model_database.loc[model_database.SHORTCODE == args.agent, "OptimizationAlgorithm"].values[0] == "CG+NUTS":
        model, data = fitAgent.fit(
            action_set,
            reward_set,
            niters=args.n_steps,
            ntune=args.n_burn,
            nchains=args.n_chains,
            nparallel=args.n_cores,
        )
    else:
        raise ValueError("Optimization algorithm not recognized.")

    print("Fitting finished. Saving model.")

    data.to_netcdf(args.save_path + args.agent + "_" + timestring + ".nc")

    with open(args.save_path + "../fits/model_fitting_log.csv", "a") as buffer:
        buffer.write(
            args.agent + "," + timestring + "," + str(args.n_burn) + "," + str(args.n_steps) + "," + str(args.n_chains)
        )
        buffer.write("\n")

    print("Model saved.")

else:
    assert args.n_splits > 1, "Number of splits must be greater than 1."
    n_sessions = np.loadtxt(args.action_set_data, delimiter=",").shape[0]

    assert n_sessions % args.n_splits == 0, "Number of splits must divide the number of sessions."
    n_sessions_per_split = n_sessions // args.n_splits

    if args.action_set_data.startswith("http"):
        # use requests to download the data
        response = requests.get(args.action_set_data)
        response.raise_for_status()
        full_action_set = np.loadtxt(io.BytesIO(response.content), delimiter=",", dtype=np.int32)
    else:
        full_action_set = np.loadtxt(args.action_set_data, delimiter=",", dtype=np.int32)
    if args.reward_set_data.startswith("http"):
        # use requests to download the data
        response = requests.get(args.reward_set_data)
        response.raise_for_status()
        full_reward_set = np.loadtxt(io.BytesIO(response.content), delimiter=",", dtype=np.int32)
    else:
        full_reward_set = np.loadtxt(args.reward_set_data, delimiter=",", dtype=np.int32)

    print("Loaded action and reward sets.")

    assert len(full_action_set) == len(full_reward_set), "Action and reward set must have the same length."
    n_trials_in_each_session = set([len(full_action_set[i]) for i in range(n_sessions)])

    assert len(n_trials_in_each_session) == 1, "All sessions must have the same number of trials."
    n_trials_per_session = n_trials_in_each_session.pop()

    print("Found {} sessions with {} trials each.".format(n_sessions, n_trials_per_session))

    for i in range(args.n_splits):
        print("Starting split {} of {}.".format(i + 1, args.n_splits))

        action_set = full_action_set[i * n_sessions_per_split : (i + 1) * n_sessions_per_split]
        reward_set = full_reward_set[i * n_sessions_per_split : (i + 1) * n_sessions_per_split]

        n_sessions = len(action_set)

        fitAgent = agentClass(ymaze_static(n_trials_per_session))

        extra_code = model_database.loc[model_database.SHORTCODE == args.agent, "FittingExtras"].values[0]
        if extra_code != "None":
            print("Found extra code for agent initialization.")
            exec(extra_code)

        print("Agent successfully initialized. Starting fitting.")

        if model_database.loc[model_database.SHORTCODE == args.agent, "OptimizationAlgorithm"].values[0] == "NUTS":
            model, data = fitAgent.fit(
                action_set,
                reward_set,
                niters=args.n_steps,
                ntune=args.n_burn,
                nchains=args.n_chains,
                nparallel=args.n_cores,
                target_accept=args.target_acceptance_rate,
            )
        elif model_database.loc[model_database.SHORTCODE == args.agent, "OptimizationAlgorithm"].values[0] == "CG+NUTS":
            model, data = fitAgent.fit(
                action_set,
                reward_set,
                niters=args.n_steps,
                ntune=args.n_burn,
                nchains=args.n_chains,
                nparallel=args.n_cores,
            )
        else:
            raise ValueError("Optimization algorithm not recognized.")

        print("Fitting finished. Saving model.")

        data.to_netcdf(args.save_path + args.agent + "_" + str(i) + "_" + timestring + ".nc")

        print("Model saved.")

    with open(args.save_path + "model_fitting_log.csv", "a") as buffer:
        buffer.write(
            args.agent + "," + timestring + "," + str(args.n_burn) + "," + str(args.n_steps) + "," + str(args.n_chains)
        )
        buffer.write("\n")

    print("All splits finished.")

print("Fitting is complete. The model fitting log is available at {}.".format(args.save_path + "model_fitting_log.csv"))
print("The model is available at {}.".format(args.save_path + args.agent + "_" + timestring + ".pkl"))
print("Thank you for using flymazerl. Have a nice day :)")

max_stars_repo_name: neurorishika/FlYMazeRL

max_issues_repo_name: neurorishika/FlYMazeRL

max_forks_repo_name: neurorishika/FlYMazeRL


text:
"""Icestupa class function that returns additional field data_h to initialise ice radius and height and include drone validation
"""

import pandas as pd
import math
import numpy as np
from functools import lru_cache
from datetime import datetime
import logging
import coloredlogs

logger = logging.getLogger("__main__")

def get_calibration(site, input):
    df_c = pd.read_csv(
        input + "drone.csv",
        sep=",",
        header=0,
        parse_dates=["time"],
    )
    df_c = df_c.reset_index()
    return df_c

    # if site in ["guttannen21", "guttannen20"]:
    #     df_c = df_c.set_index("time").sort_index().reset_index()
    #     df_cam = pd.read_csv(
    #         input  + "cam_temp.csv",
    #         sep=",",
    #         header=0,
    #         parse_dates=["time"],
    #     )

    #     # Correct thermal cam temp.
    #     if site == "guttannen20":
    #         df_c = df_c.reset_index()
    #         mask = (df_cam["time"] >= datetime(2020, 1, 2)) & (
    #             df_cam["time"] <= datetime(2020, 2, 16)
    #         )  # No ice
    #         df_cam = df_cam.loc[mask]
    #         df_cam = df_cam.reset_index(drop=True)

    #     # Correct thermal cam temp.
    #     if site == "guttannen21":
    #         mask = (df_cam["time"] >= datetime(2020, 12, 5)) & (
    #             df_cam["time"] <= datetime(2021, 3, 25)
    #         )  # No ice
    #         mask2 = (df_cam["time"] >= datetime(2020, 12, 19, 12)) & (
    #             df_cam["time"] <= datetime(2020, 12, 26, 17)
    #         )  # No ice
    #         df_cam = df_cam.loc[mask & ~mask2]
    #         df_cam = df_cam.reset_index(drop=True)

    #     mask = df_cam["std"] >= 1  # Cloudy and sunny times
    #     mask &= df_cam["mean"] >= -7  # Bluish images
    #     df_cam = df_cam.loc[mask]
    #     df_cam = df_cam.reset_index(drop=True)
    #     df_cam.loc[
    #         (df_cam["cam_temp"] > 0), "cam_temp"
    #     ] = np.NaN  # Sunlight causes this
    #     df_cam = df_cam[["time", "cam_temp"]]
    #     df_cam = df_cam.dropna()
    #     df_cam = df_cam.set_index("time")

    #     return df_c, df_cam
    # else:


text:
# AUTOGENERATED! DO NOT EDIT! File to edit: 00_hagerstrand.ipynb (unless otherwise specified).

__all__ = ['Diffusion', 'SimpleDiffusion', 'AdvancedDiffusion']

# Cell
import sys
from random import randint
from random import uniform
import numpy as np
from scipy.spatial.distance import cdist
from skimage import data, io, filters
sys.setrecursionlimit(11500)

# Cell
class Diffusion(object):
    """General class for all types of diffusion"""
    #por lo pronto solo la creación del espacio se deriva a las clases hijas?
    def __init__(self,mif_size=5,pob=20,initial_diff=[(50,50)],
                p0=0.3, max_iter=15):
        self._pob = pob
        self._p0 = p0
        self.max_iter = max_iter
        self.mif_size = mif_size
        self.iteration = 0
        self._infected_pop = []
        self._tmp_adopted = []
        self._clean = False
        self._initial_diff = initial_diff
        self.time_series = []
        self.mif_size = mif_size

    def initialize_mif(self,mif_size):
        """Initialize the MIF"""
        x = np.linspace(0.5,mif_size - 0.5,mif_size)
        y = np.linspace(0.5,mif_size - 0.5,mif_size)
        xv,yv = np.meshgrid(x,y)
        points = np.array(list(zip(np.ravel(xv),np.ravel(yv))))
        center = np.array([[mif_size/2 + 0.5,mif_size/2 + 0.5]])
        #print(points)
        #print(center)
        dist = cdist(center,points)
        dist = dist/np.sum(dist)
        #Everything: has to be different to respect the user's p0
        # print(type(mif_size), type(mif_size/2), mif_size/2)
        dist.reshape(mif_size, mif_size)[int(mif_size/2 + 0.5), int(mif_size/2 + 0.5)] = self._p0
        dist = dist/np.sum(dist)
        return np.cumsum(dist)

    def _mif2delta(self,index):
        """Returns a tuple with the increments to get to the propagated frame."""

        return np.unravel_index(index,(self.mif_size,self.mif_size))

    def _select_from_mif(self):
        """Returns an address (pob_adress) from the MIF."""
        rnd = uniform(0,1)
        index = np.nonzero(self._mif>rnd)[0][0]
        return self._mif2delta(index)

    def _clean_adopters(self):
        """Clean and initialize before a new simulation."""

        self._infected_pop = []
        self._tmp_adopted = []
        self._pop_array = np.zeros((len(np.ravel(self.space)),self._pob),
                                    dtype=np.bool)
        self.time_series = []
        for c in self._initial_diff:
            self.space[c[0],c[1]] = 1
            #We also modify the original settlers:
            index = self._space2pop_index(c)
            self._pop_array[index][0] = True
            self._infected_pop.append((index,0))

        self._clean = False

# Cell
class SimpleDiffusion(Diffusion):
    """Simple model of spatial diffusion based on Hägerstrand.
    1.- Homogeneous and isotropic space
    2.- A single initial diffuser
    3.- ....other assumptions...
    :param N: int Number of rows in simulation space.
    :param M: int Number of columns in simulation space.
    :param mif_size: int MIF matrix (square) size (must be non).
    :param pob: int population in each cell.
    :param initial_diff: [(int,int)] Coordinate list of start diffusers.
    :param p0: float Probability of self-diffusion.
    :param max_iter: int Maximum number of iterations.
    :attribute space: np.array(M,N,dtype=np.int8) Available space.
    :attribute _pop_array: np.array(M*N,pob,dtype=np.bool) array of population in each cell
    :attribute _infected_pop: list (space_idx,int) List of the adopting cell indices.
                                                The first entry is the flattened index of the cell
                                                in the space array and the second is the number of
                                                the settler in pop_array. That is, the list of addresses
                                                of each infected resident.
    :attribute results: np.array((M,N,max_iter)) Save the results of each iteration.
    :attribute time_series: list int Propagations for each iteration.
    :attribute _clean: bool Indicates if we have saved results.
    """

    def __init__(self,N=100,M=100,mif_size=5,pob=20,initial_diff=[(50,50)],
                p0=0.3, max_iter=15):

        super().__init__(mif_size, pob, initial_diff, p0, max_iter)
        # super(SimpleDiffusion,self).__init__(mif_size,pob,initial_diff,
        #             p0, max_iter)
        self.M = M
        self.N = N
        self.space = np.zeros((self.N,self.M),dtype=np.int8)
        self._pop_array = np.zeros((len(np.ravel(self.space)),pob),
                                    dtype=np.bool)
        self.result = np.zeros((M,N,max_iter),dtype=np.int8)
        for c in initial_diff:
            if c[0] > M or c[1] > N:
                raise ValueError("The coordinates on the starting difusors do not belong to the space")
            #Modificamos también a los pobladores originales:
            index = self._space2pop_index(c)
            self._pop_array[index][0] = True
            self._infected_pop.append((index,0))
        if self.mif_size%2 == 0:
            raise ValueError("MIF size must be non")
        else:
            self._mif = self.initialize_mif(self.mif_size)

    def initialize_mif(self,mif_size):
        return super(SimpleDiffusion,self).initialize_mif(self.mif_size)

    def _propagate(self,pob_adress):
        """It propagates towards the inhabitant in pob_adress if it is non-adopter.
        :param pob_adress: (int,int) the address of the inhabitant to propagate.
                           The first entry is the index (flattened) in space and
                           the second is the number of the settler in the cell
        """

        #checo si es no-adoptante
        if self._pop_array[pob_adress[0]][pob_adress[1]] == False:
            self._pop_array[pob_adress[0]][pob_adress[1]] = True
            self._tmp_adopted.append(pob_adress)
            #print "infecté al "  + str(pob_adress)

        else:
            pass

    def _space2pop_index(self,index):
        """Transform the index of space into the index of the pop_array.
        :param index (int,int) the index to transform
        """
        # print(type(index), index)
        return np.ravel_multi_index(index,dims=(self.M,self.N))

    def _pop2space_index(self,index):
        """Return the tuple (i,j)  that corresponds to the flattened index."""
        return np.unravel_index(index, (self.M,self.N))

    def _mif2delta(self,index):
        """Returns a tuple with the increments to get to the propagated frame."""
        return super(SimpleDiffusion,self)._mif2delta(index)

    def _random_adress(self):
        """Returns a random address (pob_adress)."""
        return (randint(0,(self.M*self.N) - 1),randint(0,self._pob - 1))

    def _select_from_mif(self):
        """Returns an address (pob_adress) from the MIF."""
        return super(SimpleDiffusion,self)._select_from_mif()

    def _get_propagation_adress(self,adress):
        """Returns a pop_adress address propagated by the MIF"""

        #print "Propagó: " + str(adress)
        delta = self._select_from_mif()
        delta = (delta[0] - int(self.mif_size/2+0.5),delta[1] - int(self.mif_size/2+0.5))
        space_adress = self._pop2space_index(adress[0])
        prop_space_adress = (space_adress[0] + delta[0],
                              space_adress[1] + delta[1])
        try:
            habitant = randint(0,self._pob - 1)
            return (self._space2pop_index(prop_space_adress),habitant)
        except ValueError:
            return self._get_propagation_adress(adress)

    def _clean_adopters(self):
        """Clean and initialize before a new simulation."""
        return super(SimpleDiffusion,self)._clean_adopters()

    def spatial_diffusion(self):
        """Propagate the Hagerstrand way."""

        #If we already have results, we must clean and initialize
        if self._clean:
            self._clean_adopters()

        if self.iteration == (self.max_iter or
                              np.sum(self._pop_array) >= self.M*self.N*self._pob):
            print("finished")
            print("There are %i adopters out of a total of %i inhabitants" \
                    % (np.sum(self._pop_array),self.M*self.N*self._pob))
            print("The total number of iterations performed is: %i" % self.iteration)
            self.iteration = 0
            self._clean = True
            return None
        else:
            for adress in self._infected_pop:
                propagated_adress = self._get_propagation_adress(adress)
                self._propagate(propagated_adress)

            self._infected_pop.extend(self._tmp_adopted)
            #print "Hay %i adoptantes" % len(self._infected_pop)
            self.result[:,:,self.iteration] = np.sum(self._pop_array,
                                                axis=1).reshape(self.M,self.N)
            self.time_series.append(len(self._tmp_adopted))
            self.iteration += 1
            self._tmp_adopted = []
            return self.spatial_diffusion()

    def random_diffusion(self):
        """Randomly propagates in space."""

        #If we already have results, we must clean and initialize
        if self._clean:
            self._clean_adopters()

        if self.iteration == (self.max_iter or
                              np.sum(self._pop_array) >= self.M*self.N*self._pob):
            #self.space = np.sum(s._pop_array,axis=1).reshape(s.M,s.N)
            print("finished")
            print("There are %i adopters out of a total of %i inhabitants" \
                    % (np.sum(self._pop_array),self.M*self.N*self._pob))
            print("The total number of iterations performed is: %i" % self.iteration)
            self.iteration = 0
            self._clean = True
            return None
        else:
            for adress in self._infected_pop:
                rand_adress = self._random_adress()
                if adress == rand_adress:
                    #TODO: you have to change, it could happen to get twice the same
                    rand_adress = self._random_adress()

                self._propagate(rand_adress)

            self._infected_pop.extend(self._tmp_adopted)
            #print "Hay %i adoptantes" % len(self._infected_pop)
            self.result[:,:,self.iteration] = np.sum(self._pop_array,
                                                axis=1).reshape(self.M,self.N)
            self.time_series.append(len(self._tmp_adopted))
            self.iteration += 1
            self._tmp_adopted = []
            return self.random_diffusion()

    def mixed_diffusion(self,proportion=0.5):
        """ Mix the two types of diffusion.
            In each iteration he randomly chooses, according to proportion, the
            points that diffuse randomly and those that do so spatially.
             :param proportion: float Proportion of adopters who diffuse spatially.
        """

        if proportion < 0 or proportion > 1:
            raise ValueError("The proportion must be between 0 and 1.")

        #If we already have results, we must clean and initialize
        if self._clean:
            self._clean_adopters()

        if self.iteration == (self.max_iter or
                              np.sum(self._pop_array) >= self.M*self.N*self._pob):
            #self.space = np.sum(s._pop_array,axis=1).reshape(s.M,s.N)
            print("finished")
            print("There are %i adopters out of a total of %i inhabitants" \
                    % (np.sum(self._pop_array),self.M*self.N*self._pob))
            print("The total number of iterations performed is: %i" % self.iteration)
            self.iteration = 0
            self._clean = True
            return None
        else:
            for adress in self._infected_pop:
                rnd = uniform(0,1)
                if rnd <= proportion:
                    propagated_adress = self._get_propagation_adress(adress)
                    self._propagate(propagated_adress)
                else:
                    rand_adress = self._random_adress()
                    if adress == rand_adress:
                        #TODO: you have to change, it could happen to get twice the same
                        rand_adress = self._random_adress()

                    self._propagate(rand_adress)

            self._infected_pop.extend(self._tmp_adopted)
            #print "There are %i adopters %i len(self._infected_pop)
            self.result[:,:,self.iteration] = np.sum(self._pop_array,
                                                axis=1).reshape(self.M,self.N)
            self.time_series.append(len(self._tmp_adopted))
            self.iteration += 1
            self._tmp_adopted = []
            return self.mixed_diffusion(proportion)

# Cell
class AdvancedDiffusion(Diffusion):
    """Hägerstrand-based spatial diffusion model, with heterogeneous space.
    1.- Isotropic space
    2.- A single initial diffuser
    3.- .... Other assumptions ...
    :param N: int Number of rows and columns in the simulation space.
    :param mif_size: int MIF matrix size (square) (must be odd).
    :param pob: int maximum population at each cell.
    :param density: int Number of Number of initial population nuclei.
    :param amplitud: float Gaussian filter width to blur the population.
    :param initial_diff: [(int,int)] Coordinate list of start diffusers
    :param p0: float Auto-difussion probability
    :param max_iter: int Maximum number of iterations
    :attribute space: np.array(N,N,dtype=np.int8) Available space
    :attribute _pop_array: np.array(N*N,pob,dtype=np.bool) array of inhabitants in each cell
    :attribute _infected_pop: list (space_idx,int) List of adoptive cell indices.
                                The first entry is the flattened index of the cell in the space matrix
                                and the second is the number of the settler in pop_array. That is,
                                the list of addresses of each infected resident.
    :attribute results: np.array((N,N,max_iter)) Save results of each iteration.
    :attribute time_series: list int Propagation of each iteration.
    :attribute _clean: bool Indicates if Indicates if there are saved results.
    """

    def __init__(self,N=100,mif_size=5,pob=20,initial_diff=[(50,50)],
                p0=0.3, max_iter=25,density=20,amplitud=4.0):
        super(AdvancedDiffusion,self).__init__(mif_size,pob,initial_diff, p0,
                                                max_iter)
        self.N = N
        self.density = density
        self.amplitud = amplitud
        self.space = np.zeros((self.N,self.N),dtype=np.int8)
        points = self.N * np.random.random((2, self.density ** 2))
        self.space[(points[0]).astype(np.int), (points[1]).astype(np.int)] = 1
        self.space = filters.gaussian(self.space, sigma= self.N / (self.amplitud * self.density))
        # We rescale to the value of the maximum pop and convert to integer:
        self.space *= self._pob / self.space.max()
        self.space = self.space.astype(np.int8)
        self._pop_array = np.zeros((len(np.ravel(self.space)),self._pob),
                                    dtype=np.bool)
        self.result = np.zeros((self.N,self.N,max_iter),dtype=np.int8)
        for c in initial_diff:
            if c[0] > self.N or c[1] > self.N:
                raise ValueError("Coordinates of initial diffusers do not fall in space")
            # We also modify original settlers:
            index = self._space2pop_index(c)
            self._pop_array[index][0] = True
            self._infected_pop.append((index,0))

        if self.mif_size%2 == 0:
            raise ValueError("MIF size must be odd")
        else:
            self._mif = self.initialize_mif(self.mif_size)

    def _space2pop_index(self,index):
        """Transform the index of space into the index of the pop_array.
        :param index (int,int) index to transform
        """
        return np.ravel_multi_index(index,dims=(self.N,self.N))

    def _pop2space_index(self,index):
        """Returns the tuple (i, j) that corresponds to the flattened index."""
        return np.unravel_index(index,(self.N,self.N))

    def _mif2delta(self,index):
        """Returns the tuple with the increments to get to the propagated frame."""
        return super(AdvancedDiffusion,self)._mif2delta(index)

    def _select_from_mif(self):
        """Returns an address (pob_adress) from the MIF."""
        return super(AdvancedDiffusion,self)._select_from_mif()

    def _random_adress(self):
        """Returns a random address (pob_adress)."""
        i = randint(0,self.N - 1)
        j = randint(0,self.N - 1)
        pop_idx = self._space2pop_index((i,j))
        return (pop_idx,randint(0,self.space[i,j] - 1))

    def _get_propagation_adress(self,adress):
        """Returns an address propagated from the MIF (pop_adress)."""

        #print "Propagates: " + str(adress)
        delta = self._select_from_mif()
        delta = (delta[0] - int(self.mif_size/2+0.5),delta[1] - int(self.mif_size/2+0.5))
        space_adress = self._pop2space_index(adress[0])
        prop_space_adress = (space_adress[0] + delta[0],
                             space_adress[1] + delta[1])
        try:
            # print(prop_space_adress[0],prop_space_adress[1])
            # print(self.space[prop_space_adress[0],prop_space_adress[1]])
            habitant = randint(0,self.space[prop_space_adress[0],prop_space_adress[1]])
            return (self._space2pop_index(prop_space_adress),habitant)
        except ValueError as e:
            return self._get_propagation_adress(adress)

    def _propagate(self,pob_adress):
        """Propagates through inhabitant in pob_adress if it is not-adoptant.
        :param pob_adress: (int,int) The direction of inhabitant to propagate.
                            The first entry is the index (flattened) in space
                            and the second is the number of the settler in the cell
        """

        # Check if it is not-adoptant
        try:
            if self._pop_array[pob_adress[0]][pob_adress[1]] == False:
                self._pop_array[pob_adress[0]][pob_adress[1]] = True
                self._tmp_adopted.append(pob_adress)
            else:
                pass
        except IndexError:
            # This means we are infecting someone outside the space
            pass

    def _clean_adopters(self):
        """Clean and initialize before start a new simulation."""
        return super(AdvancedDiffusion,self)._clean_adopters()

    def spatial_diffusion(self):
        """Propagates Hagerstrand like."""

        # If we have results already, we must to clean and initialize
        if self._clean:
            self._clean_adopters()

        if self.iteration == (self.max_iter or
                              np.sum(self._pop_array) >= self.M*self.N*self._pob):
            print("Done")
            print("There are %i adoptants from a total of %i inhabitants" \
                    % (np.sum(self._pop_array),self.N * self.N * self._pob))
            print("The total number of iterations performed is %i" % self.iteration)
            self.iteration = 0
            self._clean = True
            return None
        else:
            for adress in self._infected_pop:
                propagated_adress = self._get_propagation_adress(adress)
                self._propagate(propagated_adress)

            self._infected_pop.extend(self._tmp_adopted)
            #print "Hay %i adoptantes" % len(self._infected_pop)
            self.result[:,:,self.iteration] = np.sum(self._pop_array,
                                                axis=1).reshape(self.N,self.N)
            self.time_series.append(len(self._tmp_adopted))
            self.iteration += 1
            self._tmp_adopted = []
            return self.spatial_diffusion()

    def random_diffusion(self):
        """Propagates randomly in space."""

        #Si ya tenemos resultados hay que limpiar e inicializar
        if self._clean:
            self._clean_adopters()

        if self.iteration == (self.max_iter or
                              np.sum(self._pop_array) >= self.N*self.N*self._pob):
            #self.space = np.sum(s._pop_array,axis=1).reshape(s.M,s.N)
            print("Done")
            print("There are %i adoptants from a total of %i inhabitantes" \
                    % (np.sum(self._pop_array),self.N*self.N*self._pob))
            print("The total number of iterations performed is %i" % self.iteration)
            self.iteration = 0
            self._clean = True
            return None
        else:
            for adress in self._infected_pop:
                rand_adress = self._random_adress()
                if adress == rand_adress:
                    #TODO: must change, it could obtain twice the same
                    rand_adress = self._random_adress()

                self._propagate(rand_adress)

            self._infected_pop.extend(self._tmp_adopted)
            self.result[:,:,self.iteration] = np.sum(self._pop_array,
                                                axis=1).reshape(self.N,self.N)
            self.time_series.append(len(self._tmp_adopted))
            self.iteration += 1
            self._tmp_adopted = []
            return self.random_diffusion()

max_forks_repo_name: JulioMorgado/hager-py


text:
import numpy as np

### CAST
# cast_[input measure]_[output measure]

def cast_pureDP_approxDP(epsilon):
    return epsilon, 0

def cast_pureDP_approxDP_fix_delta(epsilon, delta):
    """cast an `epsilon`-DP measurement to a (epsilon', `delta`)-DP measurement
    Subtlety 2: https://eprint.iacr.org/2018/277.pdf#page=3
    """
    from utils.logs import log_sub_exp
    # delta = exp(eps) - exp(eps_new)
    # eps_new = log(exp(eps) - delta)
    #         = log(exp(eps) - exp(log(delta)))
    #         = log_sub_exp(eps - log(delta))
    return log_sub_exp(epsilon, np.log(delta)), delta

def cast_pureDP_zCDP(epsilon):
    """cast a epsilon-DP measurement to a (xi, rho)-zCDP measurement
    Proposition 1.4: https://arxiv.org/pdf/1605.02065.pdf#subsubsection.1.2.1
    """
    return 0, epsilon ** 2 / 2
    # also valid under Lemma 8.3
    # return epsilon, 0

def cast_approxDP_approxzCDP(distance):
    """cast a (epsilon, delta)-DP measurement to a delta-approximate (xi, rho)-zCDP measurement
    Proposition 1.13: https://arxiv.org/pdf/1605.02065.pdf#subsubsection.1.2.6
    "for stability-based techniques"
    """
    epsilon, delta = distance
    return delta, 0, cast_pureDP_zCDP(epsilon)

def cast_approxDP_zCDP_fix_xi(distance, xi):
    """cast a (epsilon, delta)-DP measurement to a (xi, rho)-zCDP measurement where xi is fixed
    Lemma 3.5 https://arxiv.org/pdf/1605.02065.pdf#subsection.3.2
    """
    epsilon, delta = distance
    return xi, (np.sqrt(epsilon - xi + np.log(1 / delta)) - np.sqrt(np.log(1 / delta))) ** 2

def cast_approxDP_zCDP_fix_rho(distance, rho):
    """cast a (epsilon, delta)-DP measurement to a (rho, xi)-zCDP measurement where rho is fixed
    Lemma 3.5 https://arxiv.org/pdf/1605.02065.pdf#subsection.3.2
    """
    epsilon, delta = distance
    return epsilon - rho - np.sqrt(4 * rho * np.log(1 / delta)), rho

def cast_zCDP_approxDP_fix_delta(distance, delta):
    """cast a (xi, rho)-zCDP measurement to a (epsilon, delta)-DP measurement where delta is fixed
    Lemma 3.6: https://arxiv.org/pdf/1605.02065.pdf#subsection.3.2
    Looser version in Proposition 1.3
    """
    xi, rho = distance
    return xi + rho + np.sqrt(4 * rho * np.log(1 / delta)), delta

def cast_zCDP_approxDP_fix_epsilon(distance, epsilon):
    """cast a (xi, rho)-zCDP measurement to a (epsilon, delta)-DP measurement where epsilon is fixed
    Lemma 3.6: https://arxiv.org/pdf/1605.02065.pdf#subsection.3.2
    """
    xi, rho = distance
    a = (epsilon - xi - rho) / (2 * rho)

    # the third option for least dominates
    # least = np.sqrt(np.pi * rho)
    # least = 1 / (1 + a)
    least = 2 / (1 + a + np.sqrt((1 + a) ** 2 + 4 / (np.pi * rho)))

    return epsilon, np.exp(-a ** 2 * rho) * least

def cast_renyiDP_approxDP_original_fix_delta(alpha, epsilon, delta):
    """cast a (alpha, epsilon)-RDP measurement to a (epsilon, delta)-DP measurement
    Proposition 3: https://arxiv.org/pdf/1702.07476.pdf#section.4
    """
    return epsilon + np.log(1 / delta) / (alpha - 1), delta

def cast_renyiDP_approxDP_fix_delta(alpha, epsilon, delta):
    """cast a (alpha, epsilon)-RDP measurement to a (epsilon, delta)-DP measurement

    Theorem 21: https://arxiv.org/pdf/1905.09982.pdf

    `alpha` and `epsilon` may each be length-n vectors.
    The choice of alpha that minimizes the resulting epsilon is returned.
    """

    # Implementation from paper:
    # return epsilon + np.log((alpha - 1) / alpha) - (np.log(delta) + np.log(alpha)) / (alpha - 1), delta

    # Implementation from Opacus/Tensorflow:
    # with extra bounds checking and alpha selection
    from utils.renyi_sgm import get_privacy_spent
    return get_privacy_spent(alpha, epsilon, delta)

def cast_tCDP_approxDP_fix_delta(distance, delta):
    """cast a (rho, omega)-tCDP measurement to a (epsilon, delta)-DP measurement where delta is fixed
    Lemma 6: https://projects.iq.harvard.edu/files/privacytools/files/bun_mark_composable_.pdf#%5B%7B%22num%22%3A55%2C%22gen%22%3A0%7D%2C%7B%22name%22%3A%22Fit%22%7D%5D
    See Lemma 8 for more context (choice of alpha)
    """
    rho, omega = distance
    epsilon = None
    if np.log(1 / delta) <= (omega - 1) * 2 * rho:
        epsilon = rho + 2 * np.sqrt(rho + np.log(1/delta))
    else:
        epsilon = rho * omega + np.log(1 / delta) / (omega - 1)

    return epsilon, delta

max_stars_repo_name: opendp/prelim

max_issues_repo_name: opendp/prelim

max_forks_repo_name: opendp/prelim