text:
from sklearn.compose._column_transformer import _get_transformer_list
from sklearn.utils import _get_column_indices, _safe_indexing
from sklearn.utils.validation import check_is_fitted, _check_feature_names_in
from sklearn.compose import ColumnTransformer as _ColumnTransformer

class ColumnTransformer(_ColumnTransformer):
    """Drop-in replacement for Scikit-learn's ColumnTransformer

    Scikit-learn's ColumnTransformer API provides for transformers to specify
    the names of the transformed features' columns using the
    get_feature_names_out API. But it assumes that all the original columns
    that were provided to the transformer was used by the transformer. These
    columns are no longer part of Scikit-learn's 'remainder' API that allows
    for unused columns to either be dropped or passed through. This makes it
    cumbersome to use columns to create new derived columns while also avoiding
    duplicate column names etc.

    This API extension to ColumnTransformer allows for transformers to
    implement a `get_feature_names_used` function that can return which of the
    input features were actually used by it and removed from the list of
    'remainder' columns. The columns not returned by the
    `get_feature_names_used` function will be part of the 'remainder' list and
    can utilize the 'remainder' API functionality even though they were
    provided to one or more transformers. By default, if a transformer does not
    implement this function, it assumes all columns provided to the transformer
    were used by it (which is consistent with the default Scikit-learn
    ColumnTransformer API).

    Parameters
    ----------
        See Scikit-learn's ColumnTransformer

    Attributes
    ----------
        See Scikit-learn's ColumnTransformer

    Examples
    --------

    In this example, ColumnSummer adds all columns in its input to create a
    single final column. The user can name the resulting column. If the summed
    column name matches one of the original input column names, then the column
    whose name matches the final column name is considered 'replaced' by the
    summed column. Else, the summed column is considered a 'new' column.

    In the first example, scikit-learn's FunctionTransformer is subclassed and
    get_feature_names_out and get_feature_names_used are implemented for
    derived class.

    >>> import numpy as np
    >>> import pandas as pd
    >>> from sklearn_transformer_extensions.compose import ColumnTransformer
    >>> from sklearn.preprocessing import FunctionTransformer
    >>> class ColumnSummer(FunctionTransformer):
    ...     def __init__(self, name):
    ...         super().__init__(lambda X: np.c_[np.sum(X, axis=1)])
    ...         self.name = name
    ...     def get_feature_names_out(self, _):
    ...         return [self.name]
    ...     def get_feature_names_used(self, input_features):
    ...         return [c for c in input_features if c == self.name]
    >>> ct = ColumnTransformer([("sum1", ColumnSummer(name='sum01'), [0, 1]),
    ...                         ("sum2", ColumnSummer(name='sum23'), [2, 3])],
    ...                        remainder='drop', verbose_feature_names_out=False)
    >>> X = pd.DataFrame(np.array([[0., 1., 2., 2.], [1., 1., 0., 1.]]),
    ...                  columns=['c1', 'c2', 'c3', 'c4'])
    >>> print(X)
        c1   c2   c3   c4
    0  0.0  1.0  2.0  2.0
    1  1.0  1.0  0.0  1.0
    >>> ct.fit_transform(X)
    array([[1., 4.],
           [2., 1.]])
    >>> ct.get_feature_names_out()
    array(['sum01', 'sum23'], dtype='<U5')
    >>> pd.DataFrame(ct.fit_transform(X), columns=ct.get_feature_names_out())
       sum01  sum23
    0    1.0    4.0
    1    2.0    1.0
    >>> ct = ColumnTransformer([("sum1", ColumnSummer(name='sum01'), [0, 1]),
    ...                         ("sum2", ColumnSummer(name='sum23'), [2, 3])],
    ...                        remainder='passthrough', verbose_feature_names_out=False)
    >>> ct.fit_transform(X)
    array([[1., 4., 0., 1., 2., 2.],
           [2., 1., 1., 1., 0., 1.]])
    >>> ct.get_feature_names_out()
    array(['sum01', 'sum23', 'c1', 'c2', 'c3', 'c4'], dtype=object)
    >>> pd.DataFrame(ct.fit_transform(X), columns=ct.get_feature_names_out())
       sum01  sum23   c1   c2   c3   c4
    0    1.0    4.0  0.0  1.0  2.0  2.0
    1    2.0    1.0  1.0  1.0  0.0  1.0
    >>> ct = ColumnTransformer([("sum1", ColumnSummer(name='c1'), [0, 1]),
    ...                         ("sum2", ColumnSummer(name='sum23'), [2, 3])],
    ...                        remainder='passthrough', verbose_feature_names_out=False)
    >>> pd.DataFrame(ct.fit_transform(X), columns=ct.get_feature_names_out())
        c1  sum23   c2   c3   c4
    0  1.0    4.0  1.0  2.0  2.0
    1  2.0    1.0  1.0  0.0  1.0

    In the second example, ColumnSummer is implemented by subclassing the
    extension to FunctionTransformer that is implemented in this package (see
    `FunctionTransformer`).

    >>> import numpy as np
    >>> import pandas as pd
    >>> from sklearn_transformer_extensions.compose import ColumnTransformer
    >>> from sklearn_transformer_extensions.preprocessing import FunctionTransformer
    >>> class ColumnSummer(FunctionTransformer):
    ...     def __init__(self, col_name):
    ...         super().__init__(lambda X: np.c_[np.sum(X, axis=1)], col_name=col_name)
    >>> ct = ColumnTransformer([("sum1", ColumnSummer(col_name='c1'), [0, 1]),
    ...                         ("sum2", ColumnSummer(col_name='sum23'), [2, 3])],
    ...                        remainder='passthrough', verbose_feature_names_out=False)
    >>> X = pd.DataFrame(np.array([[0., 1., 2., 2.], [1., 1., 0., 1.]]),
    ...                  columns=['c1', 'c2', 'c3', 'c4'])
    >>> pd.DataFrame(ct.fit_transform(X), columns=ct.get_feature_names_out())
        c1  sum23   c2   c3   c4
    0  1.0    4.0  1.0  2.0  2.0
    1  2.0    1.0  1.0  0.0  1.0
    """

    def _validate_remainder(self, X):
        """
        Validates ``remainder`` and defines ``_remainder`` targeting
        the remaining columns.
        """
        # yapf: disable
        is_transformer = (
            hasattr(self.remainder, "fit") or hasattr(self.remainder, "fit_transform")
        ) and hasattr(self.remainder, "transform")
        if self.remainder not in ("drop", "passthrough") and not is_transformer:
            raise ValueError(
                "The remainder keyword needs to be one of 'drop', "
                "'passthrough', or estimator. '%s' was passed instead"
                % self.remainder
            )
        # yapf: enable

        check_is_fitted(self)
        input_features = _check_feature_names_in(self)

        self._n_features = X.shape[1]
        used_cols = set()
        for name, trans, _ in self.transformers:
            columns = self._transformer_to_input_indices[name]
            columns = _safe_indexing(input_features, columns)

            used_cols_ = columns
            if hasattr(trans, "get_feature_names_used"):
                used_cols_ = trans.get_feature_names_used(columns)

            used_cols.update(used_cols_)

        used_cols_idx = set(
            [i for i, k in enumerate(input_features) if k in used_cols])
        remaining = sorted(set(range(self._n_features)) - used_cols_idx)
        self._remainder = ("remainder", self.remainder, remaining)
        self._transformer_to_input_indices["remainder"] = remaining

def make_column_transformer(
    *transformers,
    remainder="drop",
    sparse_threshold=0.3,
    n_jobs=None,
    verbose=False,
    verbose_feature_names_out=True,
):
    """Drop-in replacement for Scikit-learn's make_column_transformer"""

    transformer_list = _get_transformer_list(transformers)
    return ColumnTransformer(
        transformer_list,
        n_jobs=n_jobs,
        remainder=remainder,
        sparse_threshold=sparse_threshold,
        verbose=verbose,
        verbose_feature_names_out=verbose_feature_names_out,
    )


text:
import logging
import os
import numpy as np
import random
import zipfile
import zlib
from flarestack.shared import band_mask_cache_name
from flarestack.core.energy_pdf import EnergyPDF, read_e_pdf_dict
from flarestack.core.time_pdf import TimePDF, read_t_pdf_dict
from flarestack.core.spatial_pdf import SpatialPDF
from flarestack.utils.catalogue_loader import calculate_source_weight
from scipy import sparse, interpolate
from flarestack.shared import k_to_flux

logger = logging.getLogger(__name__)

def read_injector_dict(inj_dict):
    """Ensures that injection dictionaries remain backwards-compatible

    :param inj_dict: Injection Dictionary
    :return: Injection Dictionary compatible with new format
    """

    if inj_dict != {}:

        maps = [
            ("Injection Time PDF", "injection_sig_time_pdf"),
            ("Injection Energy PDF", "injection_energy_pdf"),
            ("injection_sig_energy_pdf", "injection_energy_pdf"),
            ("inj_energy_pdf", "injection_energy_pdf"),
            ("Poisson Smear?", "poisson_smear_bool"),
            ("injection_time_pdf", "injection_sig_time_pdf"),
            ("inj_sig_time_pdf", "injection_sig_time_pdf")
        ]

        for (old_key, new_key) in maps:

            if old_key in list(inj_dict.keys()):
                logger.warning("Deprecated inj_dict key '{0}' was used. Please use '{1}' in future.".format(
                    old_key, new_key))
                inj_dict[new_key] = inj_dict[old_key]

        pairs = [
            ("injection_energy_pdf", read_e_pdf_dict),
            ("injection_sig_time_pdf", read_t_pdf_dict)
        ]

        for (key, f) in pairs:
            if key in list(inj_dict.keys()):
                inj_dict[key] = f(inj_dict[key])

        if "injection_spatial_pdf" not in inj_dict.keys():
            inj_dict["injection_spatial_pdf"] = {}

        # if "injection_bkg_time_pdf" not in inj_dict.keys():
        #     logger.warning("No 'injection_bkg_time_pdf' was specified. A 'steady' pdf will be assumed.")
        #     inj_dict["injection_bkg_time_pdf"] = {"time_pdf_name": "steady"}

    return inj_dict

class BaseInjector:
    """Base Injector Class
    """

    subclasses = {}

    def __init__(self, season, sources, **kwargs):

        kwargs = read_injector_dict(kwargs)
        self.inj_kwargs = kwargs

        logger.info("Initialising Injector for {0}".format(season.season_name))
        self.injection_band_mask = dict()
        self.season = season
        self.season.load_background_model()

        self.sources = sources

        if len(sources) > 0:
            self.weight_scale = calculate_source_weight(self.sources)

        try:
            self.sig_time_pdf = TimePDF.create(kwargs["injection_sig_time_pdf"],
                                               season.get_time_pdf())
            # self.bkg_time_pdf = TimePDF.create(kwargs["injection_bkg_time_pdf"],
            #                                    season.get_time_pdf())
            self.energy_pdf = EnergyPDF.create(kwargs["injection_energy_pdf"])
            self.spatial_pdf = SpatialPDF(kwargs["injection_spatial_pdf"],
                                          season)
        except KeyError:
            raise Exception("Injection Arguments missing. \n "
                            "'injection_energy_pdf', 'injection_time_pdf',"
                            "and 'injection_spatial_pdf' are required. \n"
                            "Found: \n {0}".format(kwargs))

        if "poisson_smear_bool" in list(kwargs.keys()):
            self.poisson_smear = kwargs["poisson_smear_bool"]
        else:
            self.poisson_smear = True

        self.n_exp = np.nan

        try:
            self.fixed_n = kwargs["fixed_n"]
        except KeyError:
            self.fixed_n = np.nan

    def calculate_n_exp(self):

        all_n_exp = np.empty((len(self.sources), 1), dtype=np.dtype(
            [('source_name', 'a30'), ('n_exp', float)]))

        for i, source in enumerate(self.sources):
            all_n_exp[i]["source_name"] = source["source_name"]
            all_n_exp[i]["n_exp"] = self.calculate_n_exp_single(source)
        return all_n_exp

    def calculate_n_exp_single(self, source):
        raise NotImplementedError

    def get_n_exp_single(self, source):

        if not isinstance(source['source_name'], bytes):
            name = bytes(source['source_name'], encoding='utf8')
        else:
            name = source['source_name']

        return np.copy(self.n_exp[self.n_exp["source_name"] == name])

    def get_expectation(self, source, scale):
        return float(self.get_n_exp_single(source)["n_exp"]) * scale

    def update_sources(self, sources):
        """Reuses an injector with new sources

        :param sources: Sources to be added
        """
        self.sources = sources
        self.weight_scale = np.sum(
                self.sources["base_weight"] * self.sources["distance_mpc"]**-2)
        self.n_exp = self.calculate_n_exp()

    def create_dataset(self, scale, angular_error_modifier=None):
        """Create a dataset based on scrambled data for background, and Monte
        Carlo simulation for signal. Returns the composite dataset. The source
        flux can be scaled by the scale parameter.

        :param scale: Ratio of Injected Flux to source flux
        :param angular_error_modifier: AngularErrorModifier to change angular errors
        :return: Simulated dataset
        """
        bkg_events = self.season.simulate_background()

        if scale > 0.:
            sig_events = self.inject_signal(scale)
        else:
            sig_events = []

        if len(sig_events) > 0:
            simulated_data = np.concatenate((bkg_events, sig_events))
        else:
            simulated_data = bkg_events

        if angular_error_modifier is not None:

            simulated_data = angular_error_modifier.pull_correct_static(simulated_data)

        return simulated_data

    def inject_signal(self, scale):
        return

    @classmethod
    def register_subclass(cls, inj_name):
        """Adds a new subclass of EnergyPDF, with class name equal to
        "energy_pdf_name".
        """
        def decorator(subclass):
            BaseInjector.subclasses[inj_name] = subclass
            return subclass

        return decorator

    @classmethod
    def create(cls, season, sources, **kwargs):

        inj_dict = read_injector_dict(kwargs)

        if "injector_name" not in inj_dict.keys():
            return cls(season, sources, **inj_dict)

        inj_name = inj_dict["injector_name"]

        if inj_name not in BaseInjector.subclasses:
            raise ValueError(f'Bad Injector name {inj_name}. '
                             f'Available options are {BaseInjector.subclasses.keys()}')
        else:
            return BaseInjector.subclasses[inj_name](season, sources, **inj_dict)

    @staticmethod
    def get_dec_and_omega(source, bandwidth):
        # Sets half width of band
        dec_width = np.sin(np.deg2rad(bandwidth))

        sinDec = np.sin(source['dec_rad'])

        # Sets a declination band above and below the source
        min_dec = max(-1, sinDec - dec_width)
        max_dec = min(1., sinDec + dec_width)
        # Gives the solid angle coverage of the sky for the band
        omega = 2. * np.pi * (max_dec - min_dec)
        return np.arcsin(dec_width), np.arcsin(min_dec), np.arcsin(max_dec), omega

@BaseInjector.register_subclass('mc_injector')
class MCInjector(BaseInjector):
    """Core Injector Class, returns a dataset on which calculations can be
    performed. This base class is tailored for injection of MC into mock
    background. This can be either MC background, or scrambled real data.
    """

    subclasses = {}

    def __init__(self, season, sources, **kwargs):
        kwargs = read_injector_dict(kwargs)
        self._mc = season.get_mc()
        BaseInjector.__init__(self, season, sources, **kwargs)

        self.injection_declination_bandwidth = self.inj_kwargs.pop('injection_declination_bandwidth', 1.5)

        try:
            self.mc_weights = self.energy_pdf.weight_mc(self._mc)
            self.n_exp = self.calculate_n_exp()

        except KeyError:
            logger.warning("No Injection Arguments. Are you unblinding?")
            pass

    def select_mc_band(self, source):
        """For a given source, selects MC events within a declination band of
        width +/- 5 degrees that contains the source. Then returns the MC data
        subset containing only those MC events.

        :param source: Source to be simulated
        :return: mc (cut): Simulated events which lie within the band
        :return: omega: Solid Angle of the chosen band
        :return: band_mask: The mask which removes events outside band
        """

        dec_width, min_dec, max_dec, omega = self.get_dec_and_omega(source, self.injection_declination_bandwidth)

        band_mask = self.get_band_mask(source, min_dec, max_dec)

        return np.copy(self._mc[band_mask]), omega, band_mask

    def get_band_mask(self, source, min_dec, max_dec):
        # Checks if the mask has already been evaluated for the source
        # If not, creates the mask for this source, and saves it
        if source["source_name"] in list(self.injection_band_mask.keys()):
            band_mask = self.injection_band_mask[source["source_name"]]
        else:
            band_mask = np.logical_and(np.greater(self._mc["trueDec"], min_dec),
                                       np.less(self._mc["trueDec"], max_dec))
            self.injection_band_mask[source["source_name"]] = band_mask

        return band_mask

    def calculate_single_source(self, source, scale):
        """Calculate the weighted MC for a single source, given a flux scale
        and a distance scale.

        :param source:
        :param scale:
        :return:
        """
        # Selects MC events lying in a +/- 5 degree declination band
        source_mc, omega, band_mask = self.select_mc_band(source)

        source_mc = self.calculate_fluence(source, scale, source_mc,
                                           band_mask, omega)

        return source_mc

    def calculate_n_exp_single(self, source):
        return np.sum(self.calculate_single_source(source, 1.)["ow"])

    def calculate_fluence(self, source, scale, source_mc, band_mask, omega):
        """Function to calculate the fluence for a given source, and multiply
        the oneweights by this. After this step, the oneweight sum is equal
        to the expected neutrino number.

        :param source: Source to be calculated
        :param scale: Flux scale
        :param source_mc: MC that is close to source
        :param band_mask: Closeness mask for MC
        :param omega: Solid angle covered by MC mask
        :return: Modified source MC
        """
        # Calculate the effective injection time for simulation. Equal to
        # the overlap between the season and the injection time PDF for
        # the source, scaled if the injection PDF is not uniform in time.
        eff_inj_time = self.sig_time_pdf.effective_injection_time(source)

        # All injection fluxes are given in terms of k, equal to 1e-9
        inj_flux = k_to_flux(source["injection_weight_modifier"] * scale)

        # Fraction of total flux allocated to given source, assuming
        # standard candles with flux proportional to 1/d^2 multiplied by the
        # sources weight

        weight = calculate_source_weight(source) / self.weight_scale

        # Calculate the fluence, using the effective injection time.
        fluence = inj_flux * eff_inj_time * weight

        # Recalculates the oneweights to account for the declination
        # band, and the relative distance of the sources.
        # Multiplies by the fluence, to enable calculations of n_inj,
        # the expected number of injected events

        source_mc["ow"] = fluence * self.mc_weights[band_mask] / omega

        return source_mc

    def inject_signal(self, scale):
        """Randomly select simulated events from the Monte Carlo dataset to
        simulate a signal for each source. The source flux can be scaled by
        the scale parameter.

        :param scale: Ratio of Injected Flux to source flux.
        :return: Set of signal events for the given IC Season.
        """
        # Creates empty signal event array
        sig_events = np.empty((0, ), dtype=self.season.get_background_dtype())

        n_tot_exp = 0

        # Loop over each source to be simulated
        for i, source in enumerate(self.sources):

            # If a number of neutrinos to inject is specified, use that.
            # Otherwise, inject based on the flux scale as normal.

            if not np.isnan(self.fixed_n):
                n_inj = int(self.fixed_n)
            else:
                n_inj = self.get_expectation(source, scale)

            n_tot_exp += n_inj

            # Simulates poisson noise around the expectation value n_inj.
            if self.poisson_smear:
                n_s = int(np.random.poisson(n_inj))
            # If there is no poisson noise, rounds n_s down to nearest integer
            else:
                n_s = int(n_inj)

            try:
                f_n_inj = float(n_inj[0])
            except (TypeError, IndexError):
                f_n_inj = float(n_inj)

            logger.debug("Injected {0} events with an expectation of {1:.2f} events for {2}".format(
                n_s, f_n_inj, source["source_name"]
            ))

            #  If n_s = 0, skips simulation step.
            if n_s < 1:
                continue

            source_mc = self.calculate_single_source(source, scale)

            # Creates a normalised array of OneWeights
            p_select = source_mc['ow'] / np.sum(source_mc['ow'])

            # Creates an array with n_signal entries.
            # Each entry is a random integer between 0 and no. of sources.
            # The probability for each integer is equal to the OneWeight of
            # the corresponding source_path.
            ind = np.random.choice(len(source_mc['ow']), size=n_s, p=p_select)

            # Selects the sources corresponding to the random integer array
            sim_ev = source_mc[ind]

            # Rotates the Monte Carlo events onto the source_path
            sim_ev = self.spatial_pdf.rotate_to_position(
                sim_ev, source['ra_rad'], source['dec_rad']
            )

            # Generates times for each simulated event, drawing from the
            # Injector time PDF.
            sim_ev["time"] = self.sig_time_pdf.simulate_times(source, n_s)

            # Joins the new events to the signal events
            sig_events = np.concatenate(
                (sig_events,
                 sim_ev[list(self.season.get_background_dtype().names)])
            )

        return sig_events

@MCInjector.register_subclass("low_memory_injector")
class LowMemoryInjector(MCInjector):
    """For large numbers of sources O(~100), saving MC masks becomes
    increasingly burdensome. As a solution, the LowMemoryInjector should be
    used instead. It will be somewhat slower, but will have much more
    reasonable memory consumption.
    """

    def __init__(self, season, sources, **kwargs):
        self.split_cats = None
        self.injection_band_paths = None
        self.band_mask_cache = None
        self.band_mask_index = None

        MCInjector.__init__(self, season, sources, **kwargs)

    def calculate_n_exp(self):

        cats, paths, m_index, s_index = band_mask_cache_name(
            self.season, self.sources, self.injection_declination_bandwidth
        )
        self.split_cats = cats
        self.injection_band_paths = paths

        if np.sum([not os.path.isfile(x) for x in self.injection_band_paths]) > 0.:
            logger.info("No saved band masks found. These will have to be made first.")
            self.make_injection_band_mask()

        self.n_exp = np.zeros((len(self.sources), 1), dtype=np.dtype(
            [('source_name', 'a30'), ('n_exp', float),
             ('mask_index', int), ("source_index", int)]))

        self.n_exp["mask_index"] = np.array(m_index).reshape(len(m_index), 1)
        self.n_exp["source_index"] = np.array(s_index).reshape(len(s_index), 1)

        for i, source in enumerate(self.sources):
            self.n_exp[i]["source_name"] = source["source_name"]
            self.n_exp[i]["n_exp"] = self.calculate_n_exp_single(source)

        return self.n_exp

    def make_injection_band_mask(self):

        for j, cat in enumerate(self.split_cats):

            path = self.injection_band_paths[j]

            try:
                os.makedirs(os.path.dirname(path))
            except OSError:
                pass

            # Make mask
            injection_band_mask = sparse.lil_matrix((len(cat),
                                                     len(self._mc)), dtype=bool)
            for i, source in enumerate(cat):
                dec_width, min_dec, max_dec, omega = self.get_dec_and_omega(source,
                                                                            self.injection_declination_bandwidth)
                band_mask = np.logical_and(np.greater(self._mc["trueDec"], min_dec),
                                           np.less(self._mc["trueDec"], max_dec))
                injection_band_mask[i, :] = band_mask
            injection_band_mask = injection_band_mask.tocsr()
            sparse.save_npz(path, injection_band_mask)

            del injection_band_mask

            logger.info(f"Saving to {path}")

    def load_band_mask(self, index):
        path = self.injection_band_paths[index]
        # logger.debug(f'type(band_mask_cache) = {type(self.band_mask_cache)}')
        del self.band_mask_cache
        logger.debug(f'loading bandmask from {path}')
        self.band_mask_cache = sparse.load_npz(path)
        self.band_mask_index = index
        # return sparse.load_npz(path)

    def get_band_mask(self, source, min_dec, max_dec):

        entry = self.get_n_exp_single(source)
        if len(entry) != 1:
            raise ValueError(f"Length of found entries for {source['source_name']} "
                             f"is {len(entry)} but should be 1!")
        mask_index = entry["mask_index"]

        if not np.logical_and(not isinstance(self.band_mask_cache, type(None)),
                              self.band_mask_index == mask_index):
            try:
                self.load_band_mask(mask_index[0])
            except (zlib.error, zipfile.BadZipFile):
                self.make_injection_band_mask()
                self.load_band_mask(mask_index[0])

            # self.load_band_mask(mask_index[0])

        # band_mask = self.load_band_mask(mask_index[0])

        return self.band_mask_cache.getrow(entry["source_index"][0]).toarray(

        )[0]

@MCInjector.register_subclass("effective_area_injector")
class EffectiveAreaInjector(BaseInjector):
    """Class for injecting signal events by relying on effective areas rather
    than pre-existing Monte Carlo simulation. This Injector should be used
    for analysing public data, as no MC is provided.
    """

    def __init__(self, season, sources, **kwargs):
        self.effective_area_f = season.load_effective_area()
        self.energy_proxy_mapping = season.load_energy_proxy_mapping()
        self.angular_res_f = season.load_angular_resolution()
        BaseInjector.__init__(self, season, sources, **kwargs)
        self.n_exp = self.calculate_n_exp()
        self.conversion_cache = dict()

    def inject_signal(self, scale):

        # Creates empty signal event array
        sig_events = np.empty((0,),
                              dtype=self.season.get_background_dtype())

        n_tot_exp = 0

        # Loop over each source to be simulated
        for i, source in enumerate(self.sources):

            # If a number of neutrinos to inject is specified, use that.
            # Otherwise, inject based on the flux scale as normal.

            if not np.isnan(self.fixed_n):
                n_inj = int(self.fixed_n)
            else:
                n_inj = self.get_expectation(source, scale)

            n_tot_exp += n_inj

            # Simulates poisson noise around the expectation value n_inj.
            if self.poisson_smear:
                n_s = int(np.random.poisson(n_inj))
            # If there is no poisson noise, rounds n_s to nearest integer
            else:
                n_s = int(n_inj)

            logger.debug("Injected {0} events with an expectation of {1:.2f} events for {2}".format(
                n_s, n_inj if isinstance(n_inj, float) else float(n_inj[0]), source["source_name"]
            ))

            #  If n_s = 0, skips simulation step.
            if n_s < 1:
                continue

            sim_ev = np.empty(
                (n_s,), dtype=self.season.get_background_dtype())

            # Fills the energy proxy conversion cache

            if source["source_name"] not in self.conversion_cache.keys():
                self.calculate_energy_proxy(source)

            # Produces random seeds, converts this using a convolution
            # of energy pdf and approximated energy proxy mapping to
            # produce final energy proxy values

            convert_f = self.conversion_cache[source["source_name"]]

            random_fraction = [random.random() for _ in range(n_s)]

            sim_ev["logE"] = convert_f(random_fraction)

            # Simulates times according to Time PDF

            sim_ev["time"] = self.sig_time_pdf.simulate_times(source, n_s)
            sim_ev["sigma"] = self.angular_res_f(sim_ev["logE"]).copy()
            sim_ev["raw_sigma"] = sim_ev["sigma"].copy()

            sim_ev = self.spatial_pdf.simulate_distribution(source, sim_ev)

            sim_ev = sim_ev[list(
                self.season.get_background_dtype().names)].copy()
            #

            # Joins the new events to the signal events
            sig_events = np.concatenate((sig_events, sim_ev))

        sig_events = np.array(sig_events)

        return sig_events

    def calculate_single_source(self, source, scale):

        # Calculate the effective injection time for simulation. Equal to
        # the overlap between the season and the injection time PDF for
        # the source, scaled if the injection PDF is not uniform in time.
        eff_inj_time = self.sig_time_pdf.effective_injection_time(source)

        # All injection fluxes are given in terms of k, equal to 1e-9
        inj_flux = k_to_flux(source["injection_weight_modifier"] * scale)

        # Fraction of total flux allocated to given source, assuming
        # standard candles with flux proportional to 1/d^2 multiplied by the
        # sources weight

        weight = calculate_source_weight(source) / self.weight_scale

        # Calculate the fluence, using the effective injection time.
        fluence = inj_flux * eff_inj_time * weight

        def source_eff_area(e):
            return self.effective_area_f(
                np.log10(e), np.sin(source["dec_rad"])) * self.energy_pdf.f(e)

        int_eff_a = self.energy_pdf.integrate_over_E(source_eff_area)

        # Effective areas are given in m2, but flux is in per cm2

        int_eff_a *= 10**4

        n_inj = fluence*int_eff_a

        return n_inj

    def calculate_n_exp_single(self, source):
        return self.calculate_single_source(source, scale=1.)

    def calculate_energy_proxy(self, source):
        # Simulates energy proxy values

        def source_eff_area(log_e):
            return self.effective_area_f(log_e,
                                         np.sin(source["dec_rad"])) * \
                   self.energy_pdf.f(log_e) * \
                   self.energy_proxy_mapping(log_e)

        start_x = np.log10(self.energy_pdf.integral_e_min)

        x_vals = np.linspace(
            start_x + 1e-7,
            np.log10(self.energy_pdf.integral_e_max),
            100
        )[1:]

        y_vals = np.array([
            self.energy_pdf.integrate_over_E(source_eff_area, upper=np.exp(x))
            for x in x_vals
        ])
        y_vals /= max(y_vals)

        f = interpolate.interp1d([0.0] + list(y_vals), [start_x] + list(x_vals))
        self.conversion_cache[source["source_name"]] = f

class MockUnblindedInjector:
    """If the data is not really to be unblinded, then MockUnblindedInjector
    should be called. In this case, the create_dataset function simply returns
    one background scramble.
    """

    def __init__(self, season, sources=np.nan, **kwargs):
        self.season = season
        self._raw_data = season.get_exp_data()

    def create_dataset(self, scale, angular_error_modifier=None):
        """Returns a background scramble

        :return: Scrambled data
        """
        seed = int(123456)
        np.random.seed(seed)

        simulated_data = self.season.simulate_background()
        if angular_error_modifier is not None:
            simulated_data = angular_error_modifier.pull_correct_static(simulated_data)

        return simulated_data

class TrueUnblindedInjector:
    """If the data is unblinded, then UnblindedInjector should be called. In
    this case, the create_dataset function simply returns the unblinded dataset.
    """
    def __init__(self, season, sources, **kwargs):
        self.season = season

    def create_dataset(self, scale, angular_error_modifier=None):

        exp_data = self.season.get_exp_data()

        if angular_error_modifier is not None:
            exp_data = angular_error_modifier.pull_correct_static(exp_data)

        return exp_data

# if __name__ == "__main__":
#     from flarestack.data.icecube.ps_tracks.ps_v002_p01 import ps_v002_p01
#     data = ps_v002_p01[0]
#     from flarestack.analyses.agn_cores.shared_agncores import agncores_cat_dir
#     cat = np.load(
#         agncores_cat_dir +
#         "radioloud_2rxs_noBL_2000brightest_srcs_weight1.npy"
#     )
#
#     LowMemoryInjector(data, cat)

ext: py

max_stars_repo_name: grburgess/flarestack

max_issues_repo_name: grburgess/flarestack

max_forks_repo_name: grburgess/flarestack

reason: import numpy,from scipy


text:
# Run forward and backward pass on a two-spin `VectorSpinAttention` module and plot `phi`.

import matplotlib.pyplot as plt

import numpy as np
import torch

from afem.models import VectorSpinModel

num_spins, dim = 2, 64

# Setup vector-spin model.
model = VectorSpinModel(
    num_spins=num_spins,
    dim=dim,
    beta=1.0,
)

# Create random inputs and initial value for auxiliary variable (see `afem.models.VectorSpinModel`).
x = (torch.randn(1, num_spins, dim) / np.sqrt(dim)).requires_grad_()
t0 = 0.5*torch.ones(num_spins)

# Run forward pass and get output magnetizations, stationary t_star value, and approximate free energy.
out = model(x, t0=t0, return_afe=True, return_magnetizations=True)

print(
    f'✨ t_star {out.t_star.detach()}\n✨ magnetizations: {out.magnetizations}\n✨ approximate free energy: {out.afe.detach()}'
)

##############
# PLOT STUFF #
##############

def filter_array(a, threshold=50):
    idx = np.where(np.abs(a) > threshold)
    a[idx] = np.nan
    return a

fig = plt.figure()
ax = plt.axes()

t_star = out.t_star[0].detach().numpy()
t_min, t_max, t_step = 0.2, 3.0, 0.01
t_range = torch.arange(t_min, t_max, t_step)

grid_x, grid_y = torch.meshgrid(torch.arange(t_min, t_max, t_step), torch.arange(t_min, t_max, t_step))
grid = torch.stack((grid_x.reshape(-1), grid_y.reshape(-1)), dim=-1)

phis = np.array(model._phi(grid, x[:1, :, :].repeat(grid.shape[0], 1, 1)
                           ).detach().reshape(grid_x.shape[0], grid_x.shape[0]).numpy())

artist = ax.contourf(grid_x, grid_y, filter_array(phis), cmap='viridis', levels=20)
ax.scatter(t_star[0], t_star[1], color="red", s=50, marker="o")
ax.axvline(x=t_star[0], linestyle='--', color="red")
ax.axhline(y=t_star[1], linestyle='--', color="red")
ax.set_title(f'$\\varphi(t)$ for two spins (t = ({t_star[0]:.4f}, {t_star[1]:.4f}))')
ax.set_xlabel('$t_0$')
ax.set_ylabel('$t_1$')
plt.colorbar(artist)
# plt.savefig('plot.png', bbox_inches='tight')
plt.show()

max_stars_repo_name: mcbal/afem

max_issues_repo_name: mcbal/afem

max_forks_repo_name: mcbal/afem

reason: import numpy


text:
# coding: utf-8
# pylint: disable=invalid-name, no-member, too-many-locals
""" create multi-layered mesh on a unit circle """
# Copyright (c) Benyuan Liu. All rights reserved.
# Distributed under the (new) BSD License. See LICENSE.txt for more info.
from __future__ import division, absolute_import, print_function

import numpy as np
import matplotlib.pyplot as plt

class MeshCircle:
    """create meshes on uniform circle"""

    def __init__(self, n_fan=6, n_layer=8, n_el=16):
        """
        Parameters
        ----------
        n_fan : int
            number of fans (see the inner most layer)
        n_layer : int
            number of layers
        n_el : int
            number of boundary electrodes (default: 16)
        """
        self.n_fan = n_fan
        self.n_layer = n_layer
        self.n_el = n_el

        # number of points per-layer
        pts_per_layer = self.n_fan * np.arange(self.n_layer + 1)
        # The number of points of the initial layer (center point: 1 node)
        pts_per_layer[0] = 1
        self.pts_per_layer = pts_per_layer

        # starting point index of each layer
        # the initial layer (put at [-1]) must start with 0
        index = np.cumsum(pts_per_layer)
        index[-1] = 0
        self.index_per_layer = index

    def create(self):
        """create pts and tri"""
        pts = self._spawn_points()
        tri = self._spawn_elements()
        el_pos = self._get_electrodes()
        return pts, tri, el_pos

    def update(self, n_fan=8, n_layer=6, n_el=16):
        """update parameters"""
        self.n_fan = n_fan
        self.n_layer = n_layer
        self.n_el = n_el

    def _get_electrodes(self):
        """return the numbering of electrodes"""
        el_start = self.index_per_layer[self.n_layer - 1]
        el_len = self.pts_per_layer[self.n_layer]

        # place electrodes uniformly on the boundary
        n = np.linspace(
            el_start, el_start + el_len, num=self.n_el, endpoint=False, dtype=np.int
        )

        # for FMMU, electrodes should be placed clockwise
        # with 1 on the right (x+)
        n = n[::-1]  # reverse the order of electrodes
        n = np.roll(n, 1)  # rotate the index by 1
        return n

    def _spawn_points(self):
        """generate points"""
        # init points
        p = [0, 0]

        # divide r uniformly axial
        delta_r = 1.0 / self.n_layer

        for i in range(1, self.n_layer + 1):
            # increment points per-layer by fans
            n = i * self.n_fan
            r = i * delta_r
            # generate points on a layer
            pts = r * self._points_on_circle(n, offset=i)
            p = np.vstack([p, pts])

        return p

    @staticmethod
    def _points_on_circle(n, offset=0, offset_enabled=False):
        """generate points on unit circle"""
        fan_angle = 2 * np.pi / n
        a = np.array([i * fan_angle for i in range(n)])
        if offset_enabled:
            a += offset * (fan_angle / 2.0)
        pts = np.array([np.cos(a), np.sin(a)]).T

        return pts

    def _spawn_elements(self):
        """connect points fan-by-fan using a fixed pattern"""

        # element connections
        e = []
        for i in range(self.n_layer):
            e_layer = self._connect_layer(i)
            e.append(e_layer)

        return np.vstack(e)

    def _connect_layer(self, i):
        """
        generate connections on the i-th layer using points
        on the i-th and (i-1)-the layers.

        Notes
        -----
        make sure the triangles are counter-clock-wise (CCW)
        """

        # points per layer (ppl) in current and previous layer
        ppl_now = self.pts_per_layer[i + 1]
        ppl_pre = self.pts_per_layer[i]

        # starting index of current and previous layer
        index_now = self.index_per_layer[i]
        index_pre = self.index_per_layer[i - 1]

        e = []
        # A circle is divided into multiple fans, we assume that
        # points per layer per fan (plpf) is monotonically increasing
        point_plpf = i + 1
        # The connectivity pattern is the same for every fan,
        # but is different regarding whether the outer point (outer_now) is
        # interior point or boundary point on the current layer of a fan.
        # When the outer point is on the boundary,
        # the inner points (inner_now) is shared by adjacent fans,
        # thus k is not increasing and inner_now is reused.
        k = 0
        for j in range(ppl_now):

            # circular numbering in current layer
            outer_now = index_now + j
            outer_next = index_now + (j + 1) % ppl_now

            # circular numbering in previous layer
            inner_now = index_pre + k
            inner_next = index_pre + (k + 1) % ppl_pre

            # every (ppl_now/n_fan) points
            mode = j % point_plpf
            if mode == 0:
                # outer points is on the boundary
                ei = [outer_now, outer_next, inner_now]
                e.append(ei)
            else:
                # outer points is the interior points of a fan
                ei = [inner_now, outer_now, inner_next]
                e.append(ei)
                ei = [outer_now, outer_next, inner_next]
                e.append(ei)
                # advance k when outer points is a interior point
                k += 1

        return e

def demo():
    """demo using unit_circle_mesh"""
    model = MeshCircle()
    p, e, el_pos = model.create()

    # the order of inner-most triangles
    print(e[[0, 1, 2]])

    _, ax = plt.subplots(figsize=(6, 6))
    ax.plot(p[:, 0], p[:, 1], "ro", markersize=5)
    for i in range(p.shape[0]):
        ax.text(p[i, 0], p[i, 1], str(i))
    ax.set_xlim([-1.2, 1.2])
    ax.set_ylim([-1.2, 1.2])
    ax.grid("on")

    _, ax = plt.subplots(figsize=(6, 6))
    ax.triplot(p[:, 0], p[:, 1], e)
    ax.plot(p[el_pos, 0], p[el_pos, 1], "ro")
    for i, el in enumerate(el_pos):
        ax.text(p[el, 0], p[el, 1], str(i + 1))
    ax.set_xlim([-1.2, 1.2])
    ax.set_ylim([-1.2, 1.2])
    ax.grid("on")

    plt.show()

if __name__ == "__main__":
    demo()

max_stars_repo_name: ccfbeltran/pyEIT

max_issues_repo_name: ccfbeltran/pyEIT

max_forks_repo_name: ccfbeltran/pyEIT


text:
from QHO import QHO
import numpy as np

#Show that QHO can find numeric solution to QHO
for E in [0.95, 1.0, 1.05]:
    eqn = QHO(E, INITIAL = np.array([1,0], dtype= float), times = np.arange(0,5,0.05))
    eqn.solve()
    eqn.plot(0, title = "Solution for E = " + str(E),x_label = "x / m", y_label = "Psi")

#show how it compares to analytical solution
eqn = QHO(5.0, INITIAL = np.array([1,0], dtype= float), times = np.arange(0,6.4,0.05))
eqn.solve()
eqn.plot_analytical()

#shows how it can find eigen energies of QHO
eigen_energy = eqn.find_eigen_energy([33,31])
eigen_energy = eqn.brute_force_eigen_energy(eigen_energy, delta = 0.00001)
print("eigen energy found: ",eigen_energy)
eqn = QHO(eigen_energy)
eqn.solve()
eqn.plot_analytical()

#find eigen energies near 1,3,5
eqn = QHO(1., INITIAL = np.array([1,0], dtype= float), times = np.arange(0,6.4,0.05))

eigen_energy = eqn.brute_force_eigen_energy(1.5, delta = 0.01)
print("start from 1.5: ", eigen_energy)

eigen_energy = eqn.brute_force_eigen_energy(5.2, delta = 0.01)
print("start from 5.2: ", eigen_energy)

eqn = QHO(1., INITIAL = np.array([0,1], dtype= float), times = np.arange(0,6.4,0.05))

eigen_energy = eqn.brute_force_eigen_energy(3.1, delta = 0.01)
print("start from 3.1: ", eigen_energy)


text:
from typing import List, Optional

import cv2
import numpy as np

from .rasterizer import Rasterizer
from .render_context import RenderContext
from .satellite_image import get_sat_image_crop_scaled
from ..geometry import rotation33_as_yaw, transform_point

class SatelliteRasterizer(Rasterizer):
    """This rasterizer takes a satellite image in its constructor and a transform from world coordinates to this image.
    When you call rasterize, it will return a crop around the agent of interest with the agent's forward vector
    pointing right for the current timestep.
    """

    def __init__(
            self,
            render_context: RenderContext,
            map_im: np.ndarray,
            world_to_aerial: np.ndarray,
            interpolation: int = cv2.INTER_LINEAR,
    ):
        """

        Arguments:
            render_context (RenderContext): Render Context
            map_im (np.ndarray): Satellite image to crop from.
            world_to_aerial (np.ndarray): Transform to go from map coordinates to satellite image pixel coordinates.
        """
        self.render_context = render_context
        self.raster_size = render_context.raster_size_px
        self.pixel_size = render_context.pixel_size_m
        self.map_im = map_im
        self.world_to_aerial = world_to_aerial
        self.interpolation = interpolation
        self.map_pixel_scale = (
            1 / np.linalg.norm(world_to_aerial[0, 0:3]) + 1 / np.linalg.norm(world_to_aerial[1, 0:3])
        ) / 2

    def rasterize(
            self,
            history_frames: np.ndarray,
            history_agents: List[np.ndarray],
            history_tl_faces: List[np.ndarray],
            agent: Optional[np.ndarray] = None,
    ) -> np.ndarray:

        # Note 1: it looks like we are assuming that yaw in ecef == yaw in sat image
        if agent is None:
            ego_translation_m = history_frames[0]["ego_translation"]
            ego_yaw_rad = rotation33_as_yaw(history_frames[0]["ego_rotation"])

        else:
            ego_translation_m = np.append(agent["centroid"], history_frames[0]["ego_translation"][-1])
            ego_yaw_rad = agent["yaw"]

        raster_from_world = self.render_context.raster_from_world(ego_translation_m, ego_yaw_rad)
        world_from_raster = np.linalg.inv(raster_from_world)

        # Transform raster center to satellite coordinates (consider also z here)
        center_in_raster_px = np.asarray(self.raster_size) * (0.5, 0.5)
        center_in_world_m = transform_point(center_in_raster_px, world_from_raster)
        center_in_aerial_px = transform_point(np.append(center_in_world_m, ego_translation_m[2]), self.world_to_aerial)

        # Note 2: there is a negation here, unknown why this is necessary.
        # My best guess is because Y is flipped, maybe we can do this more elegantly.
        sat_im = get_sat_image_crop_scaled(
            self.map_im,
            self.raster_size,
            center_in_aerial_px,
            yaw=-ego_yaw_rad,
            pixel_size=self.pixel_size,
            sat_pixel_scale=self.map_pixel_scale,
            interpolation=self.interpolation,
        )
        if not self.render_context.set_origin_to_bottom:
            sat_im = sat_im[::-1]
        return sat_im.astype(np.float32) / 255

    def to_rgb(self, in_im: np.ndarray, **kwargs: dict) -> np.ndarray:
        return (in_im * 255).astype(np.uint8)

    def num_channels(self) -> int:
        return 3


text:
"""
pyAeroStruct_problem

Developers:
-----------
- Dr. Gaetan K. W. Kenway (GKWK)

History
-------
    v. 0.1    - Complete overall of AeroProblem (GKWK, 2014)
"""
# ======================================================================
# Imports
# ======================================================================
import numpy
from . pyAero_problem import AeroProblem
from . pyStruct_problem import StructProblem

class Error(Exception):
    """
    Format the error message in a box to make it clear this
    was a expliclty raised exception.
    """
    def __init__(self, message):
        msg = '\n+'+'-'*78+'+'+'\n' + '| AeroStructProblem Error: '
        i = 26
        for word in message.split():
            if len(word) + i + 1 > 78: # Finish line and start new one
                msg += ' '*(78-i)+'|\n| ' + word + ' '
                i = 1 + len(word)+1
            else:
                msg += word + ' '
                i += len(word)+1
        msg += ' '*(78-i) + '|\n' + '+'+'-'*78+'+'+'\n'
        print(msg)
        Exception.__init__(self)

class AeroStructProblem(object):
    """
    The main purpose of this class is to represent all relevant
    information for a coupled aero-structural analysis. To this end,
    it maintains a reference to an AaeroProblem and a StructProblem.

    Parameters
    ----------
    ap : AeroProblem class instance
        An instance of the AeroProblem class defining the aerodynamic
        part of the problem.

    sp : StructProblem class instance
        An instance of the StructProblem class defining the structural
        part of the problem
        """

    def __init__(self, ap, sp, **kwargs):

        if not isinstance(ap, AeroProblem):
            raise Error("The argument for 'ap' was not an AeroProblem!")
        if not isinstance(sp, StructProblem):
            raise Error("The argument for 'sp' was not a StructProblem!")

        self.AP = ap
        self.SP = sp

        # For consistency, we require that the name of the AP and SP
        # are the same:
        if self.AP.name != self.SP.name:
            raise Error('The name of the AeroProblem and the StructProblem \
            used to create this clsss *must* be the same')
        self.name = self.AP.name
        self.funcNames = {}

    def setDesignVars(self, x):
        """
        Set the variables in the x-dict for this object.

        Parameters
        ----------
        x : dict
            Dictionary of variables which may or may not contain the
            design variable names this object needs
            """
        self.AP.setDesignVars(x)
        self.SP.setDesignVars(x)

    def addVariablesPyOpt(self, optProb):
        """
        Add the current set of variables to the optProb object.

        Parameters
        ----------
        optProb : pyOpt_optimization class
            Optimization problem definition to which variables are added
            """
        self.AP.addVariablesPyOpt(optProb)
        self.SP.addVariablesPyOpt(optProb)

    def __getitem__(self, key):

        return self.funcNames[key]

    def evalFunctions(self, funcs, evalFuncs):
        """
        Evaluate functions of the AP and SP.
        """
        self.AP.evalFunctions(funcs, evalFuncs, ignoreMissing=True)
        self.SP.evalFunctions(funcs, evalFuncs, ignoreMissing=True)

    def evalFunctionsSens(self, funcsSens, evalFuncs):
        """
        Evaluate the sensitivity of the desired functions

        Parameters
        ----------
        funcsSens : dict
            Dictionary into which the function sensitivities are saved
        evalFuncs : iterable object containing strings
            The functions that the user wants evaluated
            """
        self.AP.evalFunctionsSens(funcsSens, evalFuncs, ignoreMissing=True)
        self.SP.evalFunctionsSens(funcsSens, evalFuncs, ignoreMissing=True)

max_stars_repo_name: camader/baseclasses

max_issues_repo_name: camader/baseclasses

max_forks_repo_name: camader/baseclasses


text:
import os

import anndata as ad
import numpy as np
import pandas as pd
import scipy.sparse

import tiledbsc.util as util

# ----------------------------------------------------------------
def describe_ann_file(
    input_path: str, show_summary=True, show_types=False, show_data=False
) -> None:
    """
    This is an anndata-describer that goes a bit beyond what `h5ls` does for us.
    In particular, it shows us that for one HDF5 file we have `anndata.X` being of type `numpy.ndarray`
    while for another HDF5 file we have `anndata.X` being of type `scipy.sparse.csr.csr_matrix`.  This is
    crucial information for building I/O logic that accepts a diversity of anndata HDF5 files.
    """
    anndata = ad.read_h5ad(input_path)
    anndata.var_names_make_unique()

    print()
    print(
        f"================================================================ {input_path}"
    )

    if show_summary:
        _describe_ann_file_show_summary(anndata, input_path)
    if show_types:
        _describe_ann_file_show_types(anndata, input_path)
    if show_data:
        _describe_ann_file_show_data(anndata, input_path)

# ----------------------------------------------------------------
def _describe_ann_file_show_summary(anndata: ad.AnnData, input_path: str) -> None:

    print()
    print("----------------------------------------------------------------")
    print("ANNDATA SUMMARY:")
    print(anndata)

    print("X SHAPE  ", anndata.X.shape)
    print("OBS  LEN ", len(anndata.obs))
    print("VAR  LEN ", len(anndata.var))

    print("OBS IS A", type(anndata.obs))
    for name in anndata.obs.keys():
        print("  ", name, anndata.obs[name].dtype)
    print("VAR IS A", type(anndata.var))
    for name in anndata.var.keys():
        print("  ", name, anndata.var[name].dtype)

    try:  # not all groups have raw X
        print("RAW X SHAPE  ", anndata.raw.X.shape)
        print(" RAW OBS  LEN ", len(anndata.raw.X.obs_names))
        print("RAW VAR  LEN ", len(anndata.raw.X.var_names))
    except AttributeError:
        pass

    print("OBS  KEYS", list(anndata.obs.keys()))
    print("VAR  KEYS", list(anndata.var.keys()))

    print("OBSM KEYS", list(anndata.obsm.keys()))
    print("VARM KEYS", list(anndata.varm.keys()))
    print("OBSP KEYS", list(anndata.obsp.keys()))
    print("VARP KEYS", list(anndata.varp.keys()))

    _describe_ann_file_show_uns_summary(anndata.uns, ["uns"])

# ----------------------------------------------------------------
def _describe_ann_file_show_types(anndata: ad.AnnData, input_path: str) -> None:

    print()
    print("----------------------------------------------------------------")
    print("ANNDATA FILE TYPES:")

    namewidth = 40

    X = anndata.X
    print("%-*s %s" % (namewidth, "X/data", type(X)))
    m, n = X.shape
    print("%-*s (%d, %d)" % (namewidth, "X/data shape", m, n))
    print("%-*s %s" % (namewidth, "X/data dtype", X.dtype))
    if isinstance(X, (scipy.sparse.csr_matrix, scipy.sparse.csc_matrix)):
        density = X.nnz / (m * n)
        print("%-*s %.4f" % (namewidth, "X/data density", density))

    has_raw = False
    try:  # not all groups have raw X
        X = anndata.raw.X
        has_raw = True
    except AttributeError:
        pass

    if has_raw:
        X = anndata.raw.X
        print("%-*s %s" % (namewidth, "X/raw", type(X)))
        m, n = X.shape
        print("%-*s (%d, %d)" % (namewidth, "X/raw shape", m, n))
        print("%-*s %s" % (namewidth, "X/data dtype", X.dtype))
        if isinstance(X, (scipy.sparse.csr_matrix, scipy.sparse.csc_matrix)):
            density = X.nnz / (m * n)
            print("%-*s %.4f" % (namewidth, "X/raw density", density))

    print("%-*s %s" % (namewidth, "obs", type(anndata.obs)))
    print("%-*s %s" % (namewidth, "var", type(anndata.var)))

    for k in anndata.obsm.keys():
        print("%-*s %s" % (namewidth, "obsm/" + k, type(anndata.obsm[k])))

    for k in anndata.varm.keys():
        print("%-*s %s" % (namewidth, "varm/" + k, type(anndata.varm[k])))

    for k in anndata.obsp.keys():
        print("%-*s %s" % (namewidth, "obsp/" + k, type(anndata.obsp[k])))

    for k in anndata.varp.keys():
        print("%-*s %s" % (namewidth, "varp/" + k, type(anndata.varp[k])))

    _describe_ann_file_show_uns_types(anndata.uns, ["uns"])

# ----------------------------------------------------------------
def _describe_ann_file_show_data(anndata: ad.AnnData, input_path: str) -> None:

    print()
    print("----------------------------------------------------------------")
    print("ANNDATA FILE DATA:")

    print("X DATA", type(anndata.X), anndata.X.shape)
    print(anndata.X)

    print()
    print("----------------------------------------------------------------")
    print("OBS DATA")
    print(anndata.obs)

    print()
    print("----------------------------------------------------------------")
    print("VAR DATA")
    print(anndata.var)

    print()
    print("----------------------------------------------------------------")
    for k in anndata.obsm.keys():
        print()
        d = anndata.obsm[k]
        print(f"OBSM/{k} DATA", type(d), d.shape)
        print(d)

    print()
    print("----------------------------------------------------------------")
    for k in anndata.varm.keys():
        print()
        d = anndata.varm[k]
        print(f"varm/{k} DATA", type(d), d.shape)
        print(d)

    print()
    print("----------------------------------------------------------------")
    for k in anndata.obsp.keys():
        print()
        d = anndata.obsp[k]
        print(f"obsp/{k} DATA", type(d), d.shape)
        print(d)

    print()
    print("----------------------------------------------------------------")
    for k in anndata.varp.keys():
        print()
        d = anndata.varp[k]
        print(f"varp/{k} DATA", type(d), d.shape)
        print(d)

    _describe_ann_file_show_uns_data(anndata.uns, ["uns"])

# ----------------------------------------------------------------
def _describe_ann_file_show_uns_summary(
    uns: ad.compat.OverloadedDict, parent_path_components
) -> None:
    """
    Recursively shows summary information about the anndata.uns structure.
    """
    for key in uns.keys():
        current_path_components = parent_path_components + [key]
        value = uns[key]
        display_name = os.path.sep.join(current_path_components)
        if isinstance(value, (dict, ad.compat.OverloadedDict)):
            _describe_ann_file_show_uns_summary(value, current_path_components)
        else:
            print(display_name)

# ----------------------------------------------------------------
def _describe_ann_file_show_uns_types(uns, parent_path_components) -> None:
    """
    Recursively shows data-type information about the anndata.uns structure.
    """
    namewidth = 40
    for key in uns.keys():
        current_path_components = parent_path_components + [key]
        value = uns[key]
        display_name = os.path.sep.join(current_path_components)
        if isinstance(value, (dict, ad.compat.OverloadedDict)):
            _describe_ann_file_show_uns_types(value, current_path_components)
        elif isinstance(value, np.ndarray):
            print(
                "%-*s" % (namewidth, display_name),
                value.shape,
                type(value),
                value.dtype,
            )
        elif isinstance(value, (scipy.sparse.csr_matrix, pd.DataFrame)):
            print("%-*s" % (namewidth, display_name), value.shape, type(value))
        else:
            print("%-*s" % (namewidth, display_name), type(value))

# ----------------------------------------------------------------
def _describe_ann_file_show_uns_data(uns, parent_path_components) -> None:
    """
    Recursively shows data contained within the anndata.uns structure.
    """
    namewidth = 40
    for key in uns.keys():
        current_path_components = parent_path_components + [key]
        value = uns[key]
        display_name = os.path.sep.join(current_path_components)
        if isinstance(value, (dict, ad.compat.OverloadedDict)):
            _describe_ann_file_show_uns_data(value, current_path_components)
        elif (
            isinstance(value, np.ndarray)
            or isinstance(value, scipy.sparse.csr_matrix)
            or isinstance(value, pd.DataFrame)
        ):
            print()
            print("%-*s" % (namewidth, display_name), type(value), value.shape)
            print(value)
        else:
            print()
            print("%-*s" % (namewidth, display_name), type(value))
            print(value)

# ----------------------------------------------------------------
def _decategoricalize(anndata: ad.AnnData) -> None:
    """
    Performs an in-place typecast into types that TileDB can persist.
    """

    # If the DataFrame contains only an index, just use it as is.
    if len(anndata.obs.columns) > 0:
        new_obs = pd.DataFrame.from_dict(
            {k: util._to_tiledb_supported_array_type(v) for k, v in anndata.obs.items()}
        )
    else:
        new_obs = anndata.obs
    if len(anndata.var.columns) > 0:
        new_var = pd.DataFrame.from_dict(
            {k: util._to_tiledb_supported_array_type(v) for k, v in anndata.var.items()}
        )
    else:
        new_var = anndata.var

    for key in anndata.obsm.keys():
        anndata.obsm[key] = util._to_tiledb_supported_array_type(anndata.obsm[key])
    for key in anndata.varm.keys():
        anndata.varm[key] = util._to_tiledb_supported_array_type(anndata.varm[key])
    for key in anndata.obsp.keys():
        anndata.obsp[key] = util._to_tiledb_supported_array_type(anndata.obsp[key])
    for key in anndata.varp.keys():
        anndata.varp[key] = util._to_tiledb_supported_array_type(anndata.varp[key])

    if anndata.raw is None:  # Some datasets have no raw.
        new_raw = None
    else:
        # Note there is some code-duplication here between cooked & raw.  However anndata.raw
        # has var not directly assignable ('AttributeError: can't set attribute'), and
        # anndata.AnnData and anndata.Raw have different constructor syntaxes, and raw doesn't
        # have obs or obsm or obsp -- so, it turns out to be simpler to just repeat ourselves a
        # little.

        new_raw_var = anndata.raw.var
        # If the DataFrame contains only an index, just use it as is.
        if len(anndata.raw.var.columns) > 0:
            new_raw_var = pd.DataFrame.from_dict(
                {
                    k: util._to_tiledb_supported_array_type(v)
                    for k, v in anndata.raw.var.items()
                }
            )

        for key in anndata.raw.varm.keys():
            anndata.raw.varm[key] = util._to_tiledb_supported_array_type(
                anndata.raw.varm[key]
            )

        new_raw = ad.Raw(
            anndata,
            X=anndata.raw.X,
            var=new_raw_var,
            varm=anndata.raw.varm,
        )

    anndata = ad.AnnData(
        X=anndata.X,
        dtype=None if anndata.X is None else anndata.X.dtype,  # some datasets have no X
        obs=new_obs,
        var=new_var,
        obsm=anndata.obsm,
        obsp=anndata.obsp,
        varm=anndata.varm,
        varp=anndata.varp,
        raw=new_raw,
        uns=anndata.uns,
    )

    return anndata

max_stars_repo_name: single-cell-data/TileDB-SingleCell

max_issues_repo_name: single-cell-data/TileDB-SingleCell

max_forks_repo_name: single-cell-data/TileDB-SingleCell

reason: import numpy,import scipy


text:
import matplotlib
# %matplotlib inline
import matplotlib.patches as patches
import matplotlib.pyplot as plt
import numpy as np
import os
import pandas as pd
import skimage.measure as measure
import sys
from step1 import *
from full_prep import lumTrans
from layers import nms,iou
datapath = '/home/zhaojie/zhaojie/Lung/DSB_Code/DSB2017-master/training/Data/ForTest/Preprocess/'
segpath = '/home/zhaojie/zhaojie/Lung/DSB_Code/DSB2017-master/training/Data/ForTest/TestSeg/'
pbbpath = '/home/zhaojie/zhaojie/Lung/code/detector_py3/results/dpn3d26/Test_Prediction_GGO/TestBbox-16-128-160/'
ct_lungsegpath = '/home/zhaojie/zhaojie/Lung/DSB_Code/DSB2017-master/training/Data/ForTest/TestLungSeg/'
def GetCenterOfNodule(seg_array):
    NoduleBox_list = []
    seg_array = measure.label(seg_array, 4)
    # print(np.max(pred_seg))
    for indexx in range(np.max(seg_array)):
        indexxx = indexx+1
        seg_array_copy = seg_array.copy()
        seg_array_copy[seg_array == indexxx] = 1
        seg_array_copy[seg_array != indexxx] = 0
        z1 = np.any(seg_array_copy, axis=(1, 2))
        ZDstart_slice, ZUend_slice = np.where(z1)[0][[0, -1]]
        z2 = np.any(seg_array_copy, axis=(0, 1))
        Lstart_slice, Rend_slice = np.where(z2)[0][[0, -1]]
        z3 = np.any(seg_array_copy, axis=(0, 2))
        Dstart_slice, Uend_slice = np.where(z3)[0][[0, -1]]
        NoduleBox = [(ZUend_slice + ZDstart_slice)//2,(Uend_slice + Dstart_slice)//2, (Rend_slice + Lstart_slice)//2, max((Uend_slice-Dstart_slice)//2, (Rend_slice-Lstart_slice)//2)]
        # print(ZDstart_slice,ZUend_slice,Dstart_slice,Uend_slice,Lstart_slice,Rend_slice + 1)
        # print('NoduleBox',NoduleBox)
        if NoduleBox[-1] !=0:
            NoduleBox_list.append(NoduleBox)
    return NoduleBox_list

imglist = sorted([f for f in os.listdir(pbbpath) if f.endswith('_pbb.npy')])
print('imglist',len(imglist))
for i in range(len(imglist)):
    img  = np.load(os.path.join(datapath,imglist[i].split('_')[0] + '_clean.npy'))
    imgLungmask  = np.load(os.path.join(ct_lungsegpath,imglist[i].split('_')[0] + '_Lungmask.npy'))
    pbb  = np.load(os.path.join(pbbpath,imglist[i]))
    save_dir = os.path.join('./GGONodulePrediction',imglist[i].split('_')[0])
    if not os.path.exists(save_dir):
        os.makedirs(save_dir)
    ##########将结节金标准保存
    # noduleMask  = np.load(os.path.join(ct_lungsegpath,imglist[i].split('_')[0] + '_Nodulemask.npy'))[0]
    # NoduleBox_list = GetCenterOfNodule(noduleMask)
    print('-----------------------------')
    # for seg_num in range(len(NoduleBox_list)):
        # boxS = NoduleBox_list[seg_num]
        # print('NoduleBox',boxS)#
        # ax = plt.subplot(1,1,1)
        # plt.imshow(img[0,boxS[0]],'gray')
        # plt.axis('off')
        # rect = patches.Rectangle((boxS[2]-boxS[3],boxS[1]-boxS[3]),boxS[3]*2,boxS[3]*2,linewidth=2,edgecolor='blue',facecolor='none')
        # ax.add_patch(rect)
        # plt.savefig(os.path.join(save_dir,imglist[i].split('_')[0] + '---' + str(boxS[0]) + '.png'))
        # plt.close()
    if pbb.shape[0] != 0 :
        max = np.max(pbb[:,0])
        thes = -10
        # thes = -1
        if max < -1:
            thes = max
        pbb = pbb[pbb[:,0]>=thes]
        pbb = nms(pbb,0.1)
        print('pbb',imglist[i].split('_')[0],pbb.shape, img.shape)

        for box_num in range(pbb.shape[0]):
            box = pbb[box_num].astype('int')[1:]#取第一个box的xyzd
            # if imgLungmask[0,box[0],box[1],box[2]] != 0:
            if box[0] < img.shape[1]:
                print('box',box)#array([169, 153,  69,  13])

                ax = plt.subplot(1,1,1)
                plt.imshow(img[0,box[0]],'gray')
                plt.axis('off')
                rect = patches.Rectangle((box[2]-box[3],box[1]-box[3]),box[3]*2,box[3]*2,linewidth=2,edgecolor='red',facecolor='none')
                ax.add_patch(rect)
                plt.savefig(os.path.join(save_dir,imglist[i].split('_')[0] + '---' + str(box[0]) + '.png'))
                if box_num > 0:
                    if box_num < pbb.shape[0]-1:#倒数第二层
                        if box[0] != pbb[box_num + 1].astype('int')[1]:#每个box和下一张不在同一层
                            plt.close()
                    elif box_num == pbb.shape[0]-1:#最后一层
                        plt.close()
                elif box[0] != pbb[box_num + 1].astype('int')[1]:#第一个box和下一张box不在同一层
                    plt.close()


text:
from mpi4py import MPI
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import random

comm = MPI.COMM_WORLD
rank = comm.Get_rank()
size = comm.Get_size()

# Function for means comparation
def are_means_equals(array1, array2):
    if len(array1) == len(array2):
        for i in range(len(array1)):
            if np.array_equal(array1[i], array2[i]) == False:
                return False
        return True
    else:
        return False

def kmeans(dataset, n):
    # ***** Algorithm start *****

    # Define some variables
    ds_dimensions = dataset.shape[1]
    ds_size = dataset.shape[0]

    # Select n random elements from the dataset
    means = []
    for i in range(n):
        means.append(dataset[random.randint(0,ds_size-1)])

    # Loop for defining clusters
    do = True
    new_means = means
    iterations = 0
    while are_means_equals(means, new_means) == False or do:
        means = new_means.copy()
        clusters = [[] for i in range(n)]

        for e in dataset:
            distances = []
            for m in means:
                distances.append(sum([(m[i]-e[i])**2 for i in range(ds_dimensions)]))
            clusters[distances.index(min(distances))].append(e)

        new_means = []
        for c in clusters:
            mean = []
            for i in range(ds_dimensions):
                mean.append(sum([point[i] for point in c])/len(c))
            new_means.append(np.array(mean))

        do = False
        iterations += 1

    for i in range(len(clusters)):
        clusters[i] = np.array(clusters[i])

    variances = [np.var(c) for c in clusters]
    average_variance = sum(variances)/len(variances)

    return [clusters, average_variance, iterations]

if __name__ == "__main__":
    if rank == 0:
        # Setting up the dataset and number of clusters (n) variable
        n = int(input('Introduce la cantidad de clusters: '))
        # n = 3
        filename = input('Introduce el nombre del archivo con el dataset a utilizar: ')
        csv = pd.read_csv(filename, header=None)
        csv = csv.loc[:, [i for i in range(csv.shape[1]) if ((csv[i].dtype == 'float64' or csv[i].dtype == 'int64'))]]
        # iris = pd.read_csv('/home/fernando-dlz/Documents/python/computo avanzado/multithreading/data/wine.data', header=None)
        print("Archivo cargado exitosamente.")

        cstm_colms = None
        while (cstm_colms != 'si' and cstm_colms != 'no'):
            cstm_colms = input('¿Deseas especificar las columnas del dataset a utilizar (de no ser asi se utilizarán todas las columnas)?: ')
            if cstm_colms == 'si':
                colms = input('Introduce las columnas a utilizar, separadas por coma (ej:1,3,4): ').split(',')
                colms = [int(e) for e in colms]
                print(f"Se utilizaran las columas {colms}")
                dataset = np.array(csv.iloc[:, colms])
            elif cstm_colms == 'no':
                print("Se utilizaran todas las columnas.")
                dataset = np.array(csv.iloc[:, :])
            else:
                print("Respuesta introducida invalida, opciones validas: [si/no]")

        # dataset = np.array(iris.loc[:, 2:3])
        ds_dimensions = dataset.shape[1]

    else:
        dataset = None
        n = 0

    dataset = comm.bcast(dataset, root=0)
    n = comm.bcast(n, root=0)

    result = kmeans(dataset, n)
    results = comm.gather(result, root=0)

    if rank == 0:
        # once all the workers finish its process we find the best result
        clusters_variances = [result[1] for result in results]
        index_best = clusters_variances.index(min(clusters_variances))

        clusters = results[index_best][0]

        print('Numero total de iteraciones antes de la convergencia: ' + str(results[index_best][2]))

        if ds_dimensions == 2 and n==3:
            # Graph designed to work specifically with 3 clusters and 2 dimensions
            plt.scatter(clusters[0][:, 0], clusters[0][:, 1], color='red', alpha=0.5)
            plt.scatter(clusters[1][:, 0], clusters[1][:, 1], color='blue', alpha=0.5)
            plt.scatter(clusters[2][:, 0], clusters[2][:, 1], color='green', alpha=0.5)
            plt.show()
        else:
            for c in clusters:
                print(c)

max_stars_repo_path: computo avanzado/mpi/kmeans-mpi.py

max_stars_repo_name: corahama/python

max_issues_repo_path: computo avanzado/mpi/kmeans-mpi.py

max_issues_repo_name: corahama/python

max_forks_repo_path: computo avanzado/mpi/kmeans-mpi.py

max_forks_repo_name: corahama/python


text:
# A wrapper to PredNet using better interface
import os
import numpy as np

from prednet import PredNet
from keras.models import Model, model_from_json
from keras import backend as K
from keras.layers import Input

class Agent():

    def read_from_json(self, json_file, weights_file):
        '''
        json_file (str): typically would be os.path.join(WEIGHTS_DIR, 'prednet_kitti_model.json')
        weights_file (str): typically would be os.path.join(WEIGHTS_DIR, 'tensorflow_weights/prednet_kitti_weights.hdf5')
        nt (int): number of images in a sequence
        '''
        f = open(json_file, 'r')
        json_string = f.read()
        f.close()
        self.train_model = model_from_json(json_string, custom_objects = {'PredNet': PredNet})
        self.train_model.load_weights(weights_file)

    def _build_test_prednet(self, output_mode='prediction'):
        # Create testing model (to output predictions)
        layer_config = self.train_model.layers[1].get_config()
        layer_config['output_mode'] = output_mode
        data_format = layer_config['data_format'] if 'data_format' in layer_config else layer_config['dim_ordering']
        self.test_prednet = PredNet(weights=self.train_model.layers[1].get_weights(), **layer_config)

        return self.test_prednet

    def output_multiple(self, seq_batch, cha_first=False, is_upscaled=True, output_mode=['prediction'], batch_size=10):
        '''
        output multiple modules' neural response
        '''
        feature = {}
        for key in output_mode:
            feature[key] = self.output(seq_batch, cha_first=cha_first, is_upscaled=is_upscaled, output_mode=key, batch_size=batch_size)

        return feature

    def output(self, seq_batch, cha_first=False, is_upscaled=True, output_mode='prediction', batch_size=10):
        '''
        input:
          seq_batch (np array, numeric, [number of sequences, number of images in each sequence, imshape[0], imshape[1], 3 channels]): if cha_first is false, the final three dimension should be imshape[0], imshape[1], 3 channels; if cha_first is ture, the final three dimensions should be 3 channels, imshape[0], imshape[1].
          is_upscaled (bool): True means the RGB value in the seq ranges from 0 to 255 and need to be normalized. The output seq_hat RGB values are in the same range as the input seq.
        '''
        self.test_prednet = self._build_test_prednet(output_mode)

        input_shape = list(self.train_model.layers[0].batch_input_shape[1:]) # find the input shape, (number of images, 3 channels, imshape[0], imshape[1]) if the channel_first = True
        input_shape[0] = seq_batch.shape[1]
        inputs = Input(shape=tuple(input_shape))
        predictions = self.test_prednet(inputs)
        test_model = Model(inputs=inputs, outputs=predictions)

        seq_wrapper = seq_batch

        if is_upscaled:
            seq_wrapper = seq_wrapper / 255

        if cha_first:
            seq_hat = test_model.predict(seq_wrapper, batch_size=batch_size)
        else:
            seq_tran = np.transpose(seq_wrapper, (0, 1, 4, 2, 3)) # make it channel first
            seq_hat = test_model.predict(seq_tran, batch_size=batch_size)
            seq_hat = np.transpose(seq_hat, (0, 1, 3, 4, 2)) # convert to original shape

        if (is_upscaled) and (output_mode == 'prediction'):
            seq_hat = seq_hat * 255

        return seq_hat

    def get_config(self):
        self._build_test_prednet()
        return self.test_prednet.get_config()

if __name__ == '__main__':
    import matplotlib.pyplot as plt
    import predusion.immaker as immaker
    from kitti_settings import *

    json_file = os.path.join(WEIGHTS_DIR, 'prednet_kitti_model.json')
    weights_file = os.path.join(WEIGHTS_DIR, 'tensorflow_weights/prednet_kitti_weights.hdf5')

    imshape = (128, 160)
    square = immaker.Square(imshape, background=150)
    seq_gener = immaker.Seq_gen()

    im = square.set_full_square(color=[255, 0, 0])
    seq_repeat = seq_gener.repeat(im, 4)[None, ...] # broacast one axis

    sub = Agent()
    sub.read_from_json(json_file, weights_file)

    ##### show the prediction
    seq_pred = sub.output(seq_repeat)
    f, ax = plt.subplots(2, 3, sharey=True, sharex=True)
    for i, sq_p, sq_r in zip(range(3), seq_pred[0], seq_repeat[0]):
        ax[0][i].imshow(sq_r.astype(int))
        ax[1][i].imshow(sq_p.astype(int))

    plt.tick_params(axis='both', which='both', bottom=False, top=False, left=False, right=False, labelbottom=False, labelleft=False)
    plt.show()

    ##### show the R2 neural activity
    r2 = sub.output(seq_repeat, output_mode='R2') # if output is not prediction, the output shape would be (number of images in a seq, a 3d tensor represent neural activation)
    print(r2.shape)

max_stars_repo_name: AgeYY/prednet

max_issues_repo_name: AgeYY/prednet

max_forks_repo_name: AgeYY/prednet


text:
from functools import reduce

import aiger_bv as BV
import networkx as nx
import numpy as np

def prob(circ, *, log=False):
    """Return probability that the output of circ is True given a valid run."""
    # Validate Input
    if circ.inputs:
        raise ValueError("All inputs must be randomized.")

    if len(circ.outputs) != 1:
        raise ValueError("Only support querying probability of single output.")
    output, *_ = circ.outputs
    if len(circ.omap[output]) != 1:
        raise ValueError("Only support single bit output queries.")

    # Note: circ represents Pr(query & valid)
    # Want: Pr(query | valid) = Pr(query & valid) / Pr(valid)
    #
    # Can answer Pr(valid) by creating another circuit such that:
    #
    #      query' = valid, and valid' = True.
    #
    valid_id = circ.circ.valid_id
    valid_test = circ.circ.aigbv.cone(valid_id)
    query_and_valid_test = (BV.uatom(1, output) & BV.uatom(1, valid_id)).aigbv
    query_and_valid_test <<= circ.circ.aigbv
    biases = circ.coin_biases

    result = _lprob(query_and_valid_test, biases)
    result -= _lprob(valid_test, biases)
    return result if log else np.exp(result)

def _lprob(circ, biases):
    # Convert circuit into labeled DAG based on BDD.
    from dfa import dfa2dict
    from bdd2dfa import to_dfa
    from aiger_bdd import to_bdd

    bdd, *_ = to_bdd(circ, renamer=lambda _, x: x)
    dag, root = dfa2dict(to_dfa(bdd, qdd=False))
    graph = nx.DiGraph()
    for node, (label, transitions) in dag.items():
        label = node.label()
        graph.add_node(node, label=label)

        if isinstance(label, bool):
            continue

        for token, node2 in transitions.items():
            graph.add_edge(node, node2, label=token)

    # View graph of MDD as circuit over LogSumExp
    lprobs = {'DUMMY': -float('inf')}
    for node in nx.topological_sort(graph.reverse()):
        val = graph.nodes[node]['label']
        if isinstance(val, bool):
            lprobs[node] = 0 if val else -float('inf')
        else:
            if graph.out_degree(node) == 1:
                left, *_ = graph.neighbors(node)
                right = 'DUMMY'
            else:
                left, right = graph.neighbors(node)

            # Swap if polarity switched.
            if not graph.edges[node, left]['label']:
                right, left = left, right

            _, idx = BV.aigbv.unpack_name(val)
            bias = float(biases[idx])
            log_biases = np.log(np.array([bias, 1 - bias]))
            kid_lprobs = np.array([lprobs[left], lprobs[right]])

            # Compute average likelihood in log scale.
            lprobs[node] = np.logaddexp(*(log_biases + kid_lprobs))

    return lprobs[node]

def coins_preimage(pcirc, *,
                   inputs=None,
                   outputs=None,
                   latchins=None,
                   latchouts=None) -> BV.UnsignedBVExpr:
    if not pcirc.has_coins:
        return BV.uatom(1, 1)

    # Default IO.
    if inputs is None:
        inputs = {}
    if outputs is None:
        outputs = {}
    if latchouts is None:
        latchouts = {}
    if latchins is None:
        latchins = pcirc.latch2init

    # Check interface.
    assert set(inputs) == pcirc.inputs
    assert set(outputs) == pcirc.outputs
    assert set(latchouts) == pcirc.latches
    assert set(latchins) == pcirc.latches

    # Encode IO.
    mappings = [
        (inputs, pcirc.circ.input_encodings),
        (outputs, pcirc.circ.output_encodings),
    ]
    for val_map, enc_map in mappings:
        val_map.update({k: e.encode(val_map[k]) for k, e in enc_map.items()})

    # Cutlatches and add to i/o.
    circ, relabels = pcirc.aigbv.cutlatches()
    relabels = {k: v for k, (v, _) in relabels.items()}  # Drop init info.
    inputs.update({relabels[k]: v for k, v in latchins.items()})
    outputs.update({relabels[k]: v for k, v in latchouts.items()})

    # Fix inputs.
    srcs = (
        BV.uatom(circ.imap[k].size, v).with_output(k).aigbv
        for k, v in inputs.items()
    )
    circ <<= reduce(lambda x, y: x | y, srcs)
    assert len(circ.inputs) == int(pcirc.has_coins)

    # Test for outputs.
    omap = circ.omap
    tests = (BV.uatom(omap[k].size, k) == v for k, v in outputs.items())
    circ >>= reduce(lambda x, y: x & y, tests).aigbv
    assert len(circ.outputs - {pcirc.circ.valid_id}) == bool(circ.outputs)

    # Test for valid.
    if len(circ.outputs) == 2:  # Force valid.
        left, right = (BV.uatom(1, o) for o in circ.outputs)
        circ >>= (left & right).aigbv

    assert len(circ.outputs) == 1
    return BV.UnsignedBVExpr(circ)

def find_coins(pcirc, *,
               inputs=None,
               outputs=None,
               latchins=None,
               latchouts=None):
    from aiger_sat.sat_bv import solve
    expr = coins_preimage(
        pcirc, inputs=inputs, outputs=outputs,
        latchins=latchins, latchouts=latchouts
    )
    model = solve(expr)
    return {k: BV.decode_int(v, signed=False) for k, v in model.items()}

__all__ = ['prob', 'coins_preimage', 'find_coins']

max_stars_repo_name: mvcisback/py-aiger-coins

max_issues_repo_name: mvcisback/py-aiger-coins

max_forks_repo_name: mvcisback/py-aiger-coins

reason: import numpy,import networkx


text:
"""
Submerged Breakwater
"""
from proteus import Domain, Context
from proteus.mprans import SpatialTools as st
from proteus import WaveTools as wt
import math
import numpy as np

opts=Context.Options([
    # predefined test cases
    ("water_level", 0.315, "Height of free surface above seabed"),
    # tank
    ("h", 0.75, "Height of domain in meters"),
    ("Ls", 2., "Distance of the front toe of the structure end from generation zone in wavelengths"),
    ("Lend", 3.5, "Distance of the back toe of the structure end from absorption zone in meters"),
    ("Lgen", 1., "Length of generation zone in wavelengths"),
    ("Labs", 2., "Length of absorption zone in wavelengths"),
    ("free_slip", True, "Should tank walls have free slip conditions "
                        "(otherwise, no slip conditions will be applied)."),
    # Geometry of a trapezoidal breakwater
    ("b", 0.25, "Width of breakwater at the top"),
    ("hs", 0.25, "Height of the breakwater"),
    ("slope1", 1./2., "Slope1 of the breakwater"),
    ("slope2", 2./3., "Slope2 of the breakwater"),
    ('porosity', 0.4, "Porosity of the medium"),
    ('d50', 0.058, "Mean diameter of the medium"),
    # waves
    ("wave_period", 1., "Period of the waves"),
    ("wave_height", 0.12, "Height of the waves"),
    ("wave_dir", np.array([1.,0.,0.]), "Direction of the waves (from left boundary)"),
    ("wave_wavelength",1.489, "Direction of the waves (from left boundary)"),
    ("wave_type", 'Fenton', "type of wave"),
    ('Ycoeff', np.array([0.23585651, 0.05332571, 0.01505139, 0.00508848,
                         0.00191607, 0.00078766, 0.00037795, 0.00027376]), 'Ycoeff only if Fenton is activated'),
    ('Bcoeff', np.array([0.24551497, 0.01803292, 0.00053150, -0.00000162,
                         0.00001059, 0.00000239, 0.00000023, 0.00000002]), 'Bcoeff only if Fenton is activated'),
    ("fast", True, "switch for fast cosh calculations in WaveTools"),
    # mesh refinement
    ("refinement", False, "Gradual refinement"),
    ("he", 0.04, "Set characteristic element size"),
    ("he_max", 10, "Set maximum characteristic element size"),
    ("he_max_water", 10, "Set maximum characteristic in water phase"),
    ("refinement_freesurface", 0.1,"Set area of constant refinement around free surface (+/- value)"),
    ("refinement_caisson", 0.,"Set area of constant refinement (Box) around caisson (+/- value)"),
    ("refinement_grading", np.sqrt(1.1*4./np.sqrt(3.))/np.sqrt(1.*4./np.sqrt(3)), "Grading of refinement/coarsening (default: 10% volume)"),
    # numerical options
    ("gen_mesh", True, "True: generate new mesh every time. False: do not generate mesh if file exists"),
    ("use_gmsh", True, "True: use Gmsh. False: use Triangle/Tetgen"),
    ("movingDomain", False, "True/False"),
    ("T", 30.0, "Simulation time"),
    ("dt_init", 0.001, "Initial time step"),
    ("dt_fixed", None, "Fixed (maximum) time step"),
    ("timeIntegration", "backwardEuler", "Time integration scheme (backwardEuler/VBDF)"),
    ("cfl", 0.5 , "Target cfl"),
    ("nsave",  5, "Number of time steps to save per second"),
    ("useRANS", 0, "RANS model"),
    ("useVF", 0., "Smoothing function at the water/air interface"),
    ("conservativeFlux", False, "Switches ON/OFF postprocessing on velocity field")])

# ----- Physical constants ----- #

rho_0 = 998.2
nu_0 = 1.004e-6
rho_1 = 1.205
nu_1 = 1.500e-5
sigma_01 = 0.0
g = np.array([0., -9.81, 0.])
gAbs = math.sqrt(sum(g**2))

# ----- CONTEXT ------ #

# waves
period = opts.wave_period
omega = 2*np.pi/opts.wave_period
height = opts.wave_height
mwl = depth = opts.water_level
direction = opts.wave_dir
wave = wt.MonochromaticWaves(period=period, waveHeight=height, mwl=mwl, depth=depth,
                             g=g, waveDir=direction,
                             wavelength=opts.wave_wavelength,
                             waveType=opts.wave_type,
                             Ycoeff=np.array(opts.Ycoeff),
                             Bcoeff=np.array(opts.Bcoeff),
                             Nf=len(opts.Bcoeff),
                             fast=opts.fast)
wavelength = wave.wavelength
waterLevel = opts.water_level

# ----- DOMAIN ----- #

domain = Domain.PlanarStraightLineGraphDomain()

# refinement
he = opts.he
smoothing = he*3.

# ----- TANK ------ #

hs = opts.hs
b = opts.b
slope1 = opts.slope1
slope2 = opts.slope2

x0 = opts.Ls*wavelength
x1 = x0 + hs/slope1
x2 = x1 + b
x3 = x2 + hs/slope2

obs = [[[x0, 0.],
        [x1, hs],
        [x2, hs],
        [x3, 0.]],]

obstacle_regions=[[x1+(x2-x1)/2., hs/2.]]

tank_dim = [x3+opts.Lend*wavelength, opts.h]

tank = st.TankWithObstacles2D(domain=domain, dim=tank_dim, obstacles=obs, hole=False,
                              obstacle_regions=obstacle_regions)

# ---- POROUS MEDIA ---- #

porosity = opts.porosity
voidFrac = 1.0-porosity
d50 = opts.d50
d15 = d50/1.2

term1 = 3.12*(10**-3.)
term2 = (gAbs/(nu_0**2.))**(2./3.)
term3 = (d15**2.)
Alpha1 = 1684+term1*term2*term3

term1 = -5.10*(10**-3.)
term2 = (gAbs/(nu_0**2.))**(1./3.)
term3 = (d15)
Beta1 = 1.72+1.57*math.exp(term1*term2*term3)
Alpha=Alpha1*nu_0*(voidFrac**2)/((porosity**3)*(d15**2))
Beta=Beta1*voidFrac/((porosity**3)*d15)

#Proteus scale in viscosity, so i need to divide alpha and beta by nu_0
dragAlpha_P=(porosity**2)*Alpha/nu_0
dragBeta=0.0 #(porosity**3)*Beta/nu_0

# ----- GENERATION / ABSORPTION LAYERS ----- #

tank.setSponge(x_n=opts.Lgen*wavelength, x_p=opts.Labs*wavelength)
dragAlpha = 10.*omega/nu_0
tank.setGenerationZones(x_n=True, waves=wave, dragAlpha=dragAlpha, smoothing = smoothing)
tank.setAbsorptionZones(x_p=True, dragAlpha = dragAlpha)
tank.setPorousZones(flags=tank.regionFlags[tank.regionIndice['obstacle1']], dragAlpha=dragAlpha_P, dragBeta=dragBeta, porosity=porosity)

# ----- BOUNDARY CONDITIONS ----- #

# Waves
tank.BC['x-'].setUnsteadyTwoPhaseVelocityInlet(wave, smoothing=smoothing, vert_axis=1)

# open top
tank.BC['y+'].setAtmosphere()

if opts.free_slip:
    tank.BC['y-'].setFreeSlip()
    tank.BC['x+'].setFreeSlip()

else:  # no slip
    tank.BC['y-'].setNoSlip()
    tank.BC['x+'].setNoSlip()

# sponge
tank.BC['sponge'].setNonMaterial()
tank.BC['obstacle1'].setNonMaterial()

for bc in tank.BC_list:
    bc.setFixedNodes()

# ----- GAUGES ----- #

# Layout 1
probes1 = []
probes1.append(((x0-0.47, 0., 0.), (x0-0.47, tank_dim[1], 0.)),)
probes1.append(((x0-0.35, 0., 0.), (x0-0.35, tank_dim[1], 0.)),)
for i in range(20):
    probes1.append(((x0+1.125+0.85+i*0.1, 0.0, 0.0), (x0+1.125+0.85+i*0.1, tank_dim[1], 0.0)),)

# Layout 2
probes2 = []
probes2.append(((x0+1.125+1.50, 0.0, 0.0),                (x0+1.125+1.50, tank_dim[1], 0.0)),)
probes2.append(((x0+1.125+1.50+0.56, 0.0, 0.0),           (x0+1.125+1.50+0.56, tank_dim[1], 0.0)),)
probes2.append(((x0+1.125+1.50+0.56+0.30, 0.0, 0.0),      (x0+1.125+1.50+0.56+0.30, tank_dim[1], 0.0)),)
probes2.append(((x0+1.125+1.50+0.56+0.30+0.16, 0.0, 0.0), (x0+1.125+1.50+0.56+0.30+0.16, tank_dim[1], 0.0)),)

# GenZone
genProbes = []
dx = opts.Lgen*wavelength/10.
for j in range(11):
    genProbes.append(((-opts.Lgen*wavelength+j*dx, 0., 0.), (-opts.Lgen*wavelength+j*dx, tank_dim[1], 0.)),)

columnLines1 = tuple(map(tuple,probes1))
columnLines2 = tuple(map(tuple,probes2))
columnLinesG = tuple(map(tuple,genProbes))

tank.attachLineIntegralGauges('vof', gauges=((('vof',),columnLines1),), fileName='line_integral_gauges_1.csv')
tank.attachLineIntegralGauges('vof', gauges=((('vof',),columnLines2),), fileName='line_integral_gauges_2.csv')
tank.attachLineIntegralGauges('vof', gauges=((('vof',),columnLinesG),), fileName='line_integral_gauges_G.csv')

# ----- ASSEMBLE DOMAIN ----- #

tank.facets = np.array([[[0, 1, 2, 3]]]+[[[i for i in range(8)]]]+[[[7, 8, 9, 6]]]+[[[4, 10, 11, 5]]])

domain.MeshOptions.use_gmsh = opts.use_gmsh
domain.MeshOptions.genMesh = opts.gen_mesh
domain.MeshOptions.he = opts.he
domain.use_gmsh = opts.use_gmsh
st.assembleDomain(domain)

# ----- REFINEMENT OPTIONS ----- #

import py2gmsh
from MeshRefinement import geometry_to_gmsh
mesh = geometry_to_gmsh(domain)

field_list = []
box_size = 0.07

box = py2gmsh.Fields.Box(mesh=mesh)
box.VIn = 0.02 #he-he/4.
box.VOut = he
box.XMin = -opts.Lgen*wavelength
box.XMax = tank_dim[0]+opts.Labs*wavelength
box.YMin = waterLevel-box_size
box.YMax = waterLevel+box_size
field_list += [box]

p0 = py2gmsh.Entity.Point([-opts.Lgen*wavelength, waterLevel+box_size, 0.], mesh=mesh)
p1 = py2gmsh.Entity.Point([tank_dim[0]+opts.Labs*wavelength, waterLevel+box_size, 0.], mesh=mesh)
p2 = py2gmsh.Entity.Point([-opts.Lgen*wavelength, waterLevel-box_size, 0.], mesh=mesh)
p3 = py2gmsh.Entity.Point([tank_dim[0]+opts.Labs*wavelength, waterLevel-box_size, 0.], mesh=mesh)
l1 = py2gmsh.Entity.Line([p0, p1], mesh=mesh)
l2 = py2gmsh.Entity.Line([p2, p3], mesh=mesh)

grading = 1.05
bl = py2gmsh.Fields.BoundaryLayer(mesh=mesh)
bl.hwall_n = 0.02 #he-he/4.
bl.ratio = grading
bl.EdgesList = [l1, l2]
field_list += [bl]

fmin = py2gmsh.Fields.Min(mesh=mesh)
fmin.FieldsList = field_list
mesh.setBackgroundField(fmin)

mesh.Options.Mesh.CharacteristicLengthMax = he

domain.MeshOptions.genMesh = opts.gen_mesh
domain.MeshOptions.use_gmsh = opts.use_gmsh
domain.use_gmsh = opts.use_gmsh

geofile = 'mesh'
mesh.writeGeo(geofile+'.geo')
domain.geofile = geofile

##########################################
# Numerical Options and other parameters #
##########################################

from math import *
from proteus import MeshTools, AuxiliaryVariables
import numpy
import proteus.MeshTools
from proteus import Domain
from proteus.Profiling import logEvent
from proteus.default_n import *
from proteus.ctransportCoefficients import smoothedHeaviside
from proteus.ctransportCoefficients import smoothedHeaviside_integral

#----------------------------------------------------
# Boundary conditions and other flags
#----------------------------------------------------
movingDomain=opts.movingDomain
checkMass=False
applyCorrection=True
applyRedistancing=True
freezeLevelSet=True

#----------------------------------------------------
# Time stepping and velocity
#----------------------------------------------------
weak_bc_penalty_constant = 10.0/nu_0#Re
dt_init = opts.dt_init
T = opts.T
nDTout = int(opts.T*opts.nsave)
timeIntegration = opts.timeIntegration
if nDTout > 0:
    dt_out= (T-dt_init)/nDTout
else:
    dt_out = 0
runCFL = opts.cfl
dt_fixed = opts.dt_fixed

#----------------------------------------------------

#  Discretization -- input options
useOldPETSc=False
useSuperlu = not True
spaceOrder = 1
useHex     = False
useRBLES   = 0.0
useMetrics = 1.0
useVF = opts.useVF #1.0
useOnlyVF = False
useRANS = opts.useRANS # 0 -- None
            # 1 -- K-Epsilon
            # 2 -- K-Omega, 1998
            # 3 -- K-Omega, 1988
# Input checks
if spaceOrder not in [1,2]:
    print "INVALID: spaceOrder" + spaceOrder
    sys.exit()

if useRBLES not in [0.0, 1.0]:
    print "INVALID: useRBLES" + useRBLES
    sys.exit()

if useMetrics not in [0.0, 1.0]:
    print "INVALID: useMetrics"
    sys.exit()

#  Discretization
nd = 2
if spaceOrder == 1:
    hFactor=1.0
    if useHex:
	 basis=C0_AffineLinearOnCubeWithNodalBasis
         elementQuadrature = CubeGaussQuadrature(nd,3)
         elementBoundaryQuadrature = CubeGaussQuadrature(nd-1,3)
    else:
    	 basis=C0_AffineLinearOnSimplexWithNodalBasis
         elementQuadrature = SimplexGaussQuadrature(nd,3)
         elementBoundaryQuadrature = SimplexGaussQuadrature(nd-1,3)
         #elementBoundaryQuadrature = SimplexLobattoQuadrature(nd-1,1)
elif spaceOrder == 2:
    hFactor=0.5
    if useHex:
	basis=C0_AffineLagrangeOnCubeWithNodalBasis
        elementQuadrature = CubeGaussQuadrature(nd,4)
        elementBoundaryQuadrature = CubeGaussQuadrature(nd-1,4)
    else:
	basis=C0_AffineQuadraticOnSimplexWithNodalBasis
        elementQuadrature = SimplexGaussQuadrature(nd,4)
        elementBoundaryQuadrature = SimplexGaussQuadrature(nd-1,4)

# Numerical parameters
sc = 0.5 # default: 0.5. Test: 0.25
sc_beta = 1.5 # default: 1.5. Test: 1.
epsFact_consrv_diffusion = 1. # default: 1.0. Test: 0.1
ns_forceStrongDirichlet = False
backgroundDiffusionFactor=0.01
if useMetrics:
    ns_shockCapturingFactor  = sc
    ns_lag_shockCapturing = True
    ns_lag_subgridError = True
    ls_shockCapturingFactor  = sc
    ls_lag_shockCapturing = True
    ls_sc_uref  = 1.0
    ls_sc_beta  = sc_beta
    vof_shockCapturingFactor = sc
    vof_lag_shockCapturing = True
    vof_sc_uref = 1.0
    vof_sc_beta = sc_beta
    rd_shockCapturingFactor  =sc
    rd_lag_shockCapturing = False
    epsFact_density    = 3.
    epsFact_viscosity  = epsFact_curvature  = epsFact_vof = epsFact_consrv_heaviside = epsFact_consrv_dirac = epsFact_density
    epsFact_redistance = 0.33
    epsFact_consrv_diffusion = epsFact_consrv_diffusion
    redist_Newton = True#False
    kappa_shockCapturingFactor = sc
    kappa_lag_shockCapturing = True
    kappa_sc_uref = 1.0
    kappa_sc_beta = sc_beta
    dissipation_shockCapturingFactor = sc
    dissipation_lag_shockCapturing = True
    dissipation_sc_uref = 1.0
    dissipation_sc_beta = sc_beta
else:
    ns_shockCapturingFactor  = 0.9
    ns_lag_shockCapturing = True
    ns_lag_subgridError = True
    ls_shockCapturingFactor  = 0.9
    ls_lag_shockCapturing = True
    ls_sc_uref  = 1.0
    ls_sc_beta  = 1.0
    vof_shockCapturingFactor = 0.9
    vof_lag_shockCapturing = True
    vof_sc_uref  = 1.0
    vof_sc_beta  = 1.0
    rd_shockCapturingFactor  = 0.9
    rd_lag_shockCapturing = False
    epsFact_density    = 1.5
    epsFact_viscosity  = epsFact_curvature  = epsFact_vof = epsFact_consrv_heaviside = epsFact_consrv_dirac = epsFact_density
    epsFact_redistance = 0.33
    epsFact_consrv_diffusion = 10.0
    redist_Newton = False#True
    kappa_shockCapturingFactor = 0.9
    kappa_lag_shockCapturing = True#False
    kappa_sc_uref  = 1.0
    kappa_sc_beta  = 1.0
    dissipation_shockCapturingFactor = 0.9
    dissipation_lag_shockCapturing = True#False
    dissipation_sc_uref  = 1.0
    dissipation_sc_beta  = 1.0

ns_nl_atol_res = 1e-6#max(1.0e-6,0.001*domain.MeshOptions.he**2)
vof_nl_atol_res = 1e-6#max(1.0e-6,0.001*domain.MeshOptions.he**2)
ls_nl_atol_res = 1e-6#max(1.0e-6,0.001*domain.MeshOptions.he**2)
mcorr_nl_atol_res = 1e-6#max(1.0e-6,0.0001*domain.MeshOptions.he**2)
rd_nl_atol_res = 1e-6#max(1.0e-6,0.01*domain.MeshOptions.he)
kappa_nl_atol_res = 1e-6#max(1.0e-6,0.001*domain.MeshOptions.he**2)
dissipation_nl_atol_res = 1e-6#max(1.0e-6,0.001*domain.MeshOptions.he**2)
mesh_nl_atol_res = 1e-6#max(1.0e-6,0.001*domain.MeshOptions.he**2)

#turbulence
ns_closure=0 #1-classic smagorinsky, 2-dynamic smagorinsky, 3 -- k-epsilon, 4 -- k-omega

if useRANS == 1:
    ns_closure = 3
elif useRANS >= 2:
    ns_closure == 4

def twpflowPressure_init(x, t):
    p_L = 0.0
    phi_L = tank_dim[nd-1] - waterLevel
    phi = x[nd-1] - waterLevel
    return p_L -g[nd-1]*(rho_0*(phi_L - phi)+(rho_1 -rho_0)*(smoothedHeaviside_integral(epsFact_consrv_heaviside*opts.he,phi_L)
                                                         -smoothedHeaviside_integral(epsFact_consrv_heaviside*opts.he,phi)))

max_stars_repo_name: AntoninoScala/air-water-vv

max_issues_repo_name: AntoninoScala/air-water-vv

max_forks_repo_name: AntoninoScala/air-water-vv


text:
#!/opt/conda/bin/python
from webtest import TestApp
import logging
import numpy as np
from io import StringIO

from bottle import Bottle, route, request, response, run, template, static_file, abort
from bikesim.simulation.suspension_simulation import simulate_damper_sweep
from bikesim.models.kinematics import BikeKinematics
import matplotlib.pyplot as plt

logging.basicConfig(level=logging.INFO)

app = Bottle()

def simulate_and_plot(damper_travel):
    bike_file = 'geometries/5010_bike.json'
    n_point = 20
    sag_array = np.linspace(0, 1, n_point)

    bike = BikeKinematics.from_json(bike_file)
    bike.damper_travel = damper_travel

    wheel_travel = simulate_damper_sweep(sag_array=sag_array, bike=bike)

    plt.plot(sag_array*100, wheel_travel, '-*')
    plt.xlabel('suspension sag [%]')
    plt.ylabel('wheel travel [m]')
    buf = StringIO()
    plt.savefig(buf, format='svg')
    return buf.getvalue()

@app.route('/')
def index():
    return static_file('html/index.html', root="./")

@app.route('/plot.svg')
def index():
    damper_travel = request.query.damper_travel
    try:
        damper_travel = float(damper_travel)
        response.content_type = 'image/svg+xml'
        return simulate_and_plot(damper_travel)
    except ValueError:
        abort(404, "Invalid damper travel: " + damper_travel)

def main():
    run(app, host='0.0.0.0', port=8080)

if __name__ == '__main__':
    main()