text:
# Copyright Contributors to the Pyro project.
# SPDX-License-Identifier: Apache-2.0

from collections import OrderedDict, namedtuple
from contextlib import ExitStack  # python 3
from enum import Enum

from jax import lax
import jax.numpy as jnp

import funsor
from numpyro.handlers import infer_config, trace as OrigTraceMessenger
from numpyro.primitives import Messenger, apply_stack, plate as OrigPlateMessenger

funsor.set_backend("jax")

__all__ = ["enum", "infer_config", "markov", "plate", "to_data", "to_funsor", "trace"]

##################################
# DimStack to store global state
##################################

# name_to_dim : dict, dim_to_name : dict, parents : tuple, iter_parents : tuple
class StackFrame(
    namedtuple(
        "StackFrame", ["name_to_dim", "dim_to_name", "parents", "iter_parents", "keep"]
    )
):
    def read(self, name, dim):
        found_name = self.dim_to_name.get(dim, name)
        found_dim = self.name_to_dim.get(name, dim)
        found = name in self.name_to_dim or dim in self.dim_to_name
        return found_name, found_dim, found

    def write(self, name, dim):
        assert name is not None and dim is not None
        self.dim_to_name[dim] = name
        self.name_to_dim[name] = dim

    def free(self, name, dim):
        self.dim_to_name.pop(dim, None)
        self.name_to_dim.pop(name, None)
        return name, dim

class DimType(Enum):
    """Enumerates the possible types of dimensions to allocate"""

    LOCAL = 0
    GLOBAL = 1
    VISIBLE = 2

DimRequest = namedtuple("DimRequest", ["dim", "dim_type"])
DimRequest.__new__.__defaults__ = (None, DimType.LOCAL)
NameRequest = namedtuple("NameRequest", ["name", "dim_type"])
NameRequest.__new__.__defaults__ = (None, DimType.LOCAL)

class DimStack:
    """
    Single piece of global state to keep track of the mapping between names and dimensions.

    Replaces the plate DimAllocator, the enum EnumAllocator, the stack in MarkovMessenger,
    _param_dims and _value_dims in EnumMessenger, and dim_to_symbol in msg['infer']
    """

    def __init__(self):
        self._stack = [
            StackFrame(
                name_to_dim=OrderedDict(),
                dim_to_name=OrderedDict(),
                parents=(),
                iter_parents=(),
                keep=False,
            )
        ]
        self._first_available_dim = -1
        self.outermost = None

    MAX_DIM = -25

    def set_first_available_dim(self, dim):
        assert dim is None or (self.MAX_DIM < dim < 0)
        old_dim, self._first_available_dim = self._first_available_dim, dim
        return old_dim

    def push(self, frame):
        self._stack.append(frame)

    def pop(self):
        assert len(self._stack) > 1, "cannot pop the global frame"
        return self._stack.pop()

    @property
    def current_frame(self):
        return self._stack[-1]

    @property
    def global_frame(self):
        return self._stack[0]

    def _gendim(self, name_request, dim_request):
        assert isinstance(name_request, NameRequest) and isinstance(
            dim_request, DimRequest
        )
        dim_type = dim_request.dim_type

        if name_request.name is None:
            fresh_name = f"_pyro_dim_{-dim_request.dim}"
        else:
            fresh_name = name_request.name

        conflict_frames = (
            (self.current_frame, self.global_frame)
            + self.current_frame.parents
            + self.current_frame.iter_parents
        )
        if dim_request.dim is None:
            fresh_dim = self._first_available_dim if dim_type != DimType.VISIBLE else -1
            fresh_dim = -1 if fresh_dim is None else fresh_dim
            while any(fresh_dim in p.dim_to_name for p in conflict_frames):
                fresh_dim -= 1
        else:
            fresh_dim = dim_request.dim

        if (
            fresh_dim < self.MAX_DIM
            or any(fresh_dim in p.dim_to_name for p in conflict_frames)
            or (dim_type == DimType.VISIBLE and fresh_dim <= self._first_available_dim)
        ):
            raise ValueError(f"Ran out of free dims during allocation for {fresh_name}")

        return fresh_name, fresh_dim

    def request(self, name, dim):
        assert isinstance(name, NameRequest) ^ isinstance(dim, DimRequest)
        if isinstance(dim, DimRequest):
            dim, dim_type = dim.dim, dim.dim_type
        elif isinstance(name, NameRequest):
            name, dim_type = name.name, name.dim_type

        read_frames = (
            (self.global_frame,)
            if dim_type != DimType.LOCAL
            else (self.current_frame,)
            + self.current_frame.parents
            + self.current_frame.iter_parents
            + (self.global_frame,)
        )

        # read dimension
        for frame in read_frames:
            name, dim, found = frame.read(name, dim)
            if found:
                break

        # generate fresh name or dimension
        if not found:
            name, dim = self._gendim(
                NameRequest(name, dim_type), DimRequest(dim, dim_type)
            )

            write_frames = (
                (self.global_frame,)
                if dim_type != DimType.LOCAL
                else (self.current_frame,)
                + (self.current_frame.parents if self.current_frame.keep else ())
            )

            # store the fresh dimension
            for frame in write_frames:
                frame.write(name, dim)

        return name, dim

_DIM_STACK = DimStack()  # only one global instance

#################################################
# Messengers that implement guts of enumeration
#################################################

class ReentrantMessenger(Messenger):
    def __init__(self, fn=None):
        self._ref_count = 0
        super().__init__(fn)

    # def __call__(self, fn):
    #     return functools.wraps(fn)(super().__call__(fn))

    def __enter__(self):
        self._ref_count += 1
        if self._ref_count == 1:
            super().__enter__()
        return self

    def __exit__(self, exc_type, exc_value, traceback):
        self._ref_count -= 1
        if self._ref_count == 0:
            super().__exit__(exc_type, exc_value, traceback)

class DimStackCleanupMessenger(ReentrantMessenger):
    def __init__(self, fn=None):
        self._saved_dims = ()
        return super().__init__(fn)

    def __enter__(self):
        if self._ref_count == 0 and _DIM_STACK.outermost is None:
            _DIM_STACK.outermost = self
            for name, dim in self._saved_dims:
                _DIM_STACK.global_frame.write(name, dim)
            self._saved_dims = ()
        return super().__enter__()

    def __exit__(self, *args, **kwargs):
        if self._ref_count == 1 and _DIM_STACK.outermost is self:
            _DIM_STACK.outermost = None
            for name, dim in reversed(
                tuple(_DIM_STACK.global_frame.name_to_dim.items())
            ):
                self._saved_dims += (_DIM_STACK.global_frame.free(name, dim),)
        return super().__exit__(*args, **kwargs)

class NamedMessenger(DimStackCleanupMessenger):
    def process_message(self, msg):
        if msg["type"] == "to_funsor":
            self._pyro_to_funsor(msg)
        elif msg["type"] == "to_data":
            self._pyro_to_data(msg)

    @staticmethod
    def _get_name_to_dim(batch_names, name_to_dim=None, dim_type=DimType.LOCAL):
        name_to_dim = OrderedDict() if name_to_dim is None else name_to_dim.copy()

        # interpret all names/dims as requests since we only run this function once
        for name in batch_names:
            dim = name_to_dim.get(name, None)
            name_to_dim[name] = (
                dim if isinstance(dim, DimRequest) else DimRequest(dim, dim_type)
            )

        # read dimensions and allocate fresh dimensions as necessary
        for name, dim_request in name_to_dim.items():
            name_to_dim[name] = _DIM_STACK.request(name, dim_request)[1]

        return name_to_dim

    @classmethod  # only depends on the global _DIM_STACK state, not self
    def _pyro_to_data(cls, msg):

        (funsor_value,) = msg["args"]
        name_to_dim = msg["kwargs"].setdefault("name_to_dim", OrderedDict())
        dim_type = msg["kwargs"].setdefault("dim_type", DimType.LOCAL)

        batch_names = tuple(funsor_value.inputs.keys())

        name_to_dim.update(
            cls._get_name_to_dim(
                batch_names, name_to_dim=name_to_dim, dim_type=dim_type
            )
        )
        msg["stop"] = True  # only need to run this once per to_data call

    @staticmethod
    def _get_dim_to_name(batch_shape, dim_to_name=None, dim_type=DimType.LOCAL):
        dim_to_name = OrderedDict() if dim_to_name is None else dim_to_name.copy()
        batch_dim = len(batch_shape)

        # interpret all names/dims as requests since we only run this function once
        for dim in range(-batch_dim, 0):
            name = dim_to_name.get(dim, None)
            # the time dimension on the left sometimes necessitates empty dimensions appearing
            # before they have been assigned a name
            if batch_shape[dim] == 1 and name is None:
                continue
            dim_to_name[dim] = (
                name if isinstance(name, NameRequest) else NameRequest(name, dim_type)
            )

        for dim, name_request in dim_to_name.items():
            dim_to_name[dim] = _DIM_STACK.request(name_request, dim)[0]

        return dim_to_name

    @classmethod  # only depends on the global _DIM_STACK state, not self
    def _pyro_to_funsor(cls, msg):

        if len(msg["args"]) == 2:
            raw_value, output = msg["args"]
        else:
            raw_value = msg["args"][0]
            output = msg["kwargs"].setdefault("output", None)
        dim_to_name = msg["kwargs"].setdefault("dim_to_name", OrderedDict())
        dim_type = msg["kwargs"].setdefault("dim_type", DimType.LOCAL)

        event_dim = len(output.shape) if output else 0
        try:
            batch_shape = raw_value.batch_shape  # TODO make make this more robust
        except AttributeError:
            batch_shape = raw_value.shape[: len(raw_value.shape) - event_dim]

        dim_to_name.update(
            cls._get_dim_to_name(
                batch_shape, dim_to_name=dim_to_name, dim_type=dim_type
            )
        )
        msg["stop"] = True  # only need to run this once per to_funsor call

class LocalNamedMessenger(NamedMessenger):
    """
    Handler for converting to/from funsors consistent with Pyro's positional batch dimensions.

    :param int history: The number of previous contexts visible from the
        current context. Defaults to 1. If zero, this is similar to
        :class:`pyro.plate`.
    :param bool keep: If true, frames are replayable. This is important
        when branching: if ``keep=True``, neighboring branches at the same
        level can depend on each other; if ``keep=False``, neighboring branches
        are independent (conditioned on their shared ancestors).
    """

    def __init__(self, fn=None, history=1, keep=False):
        self.history = history
        self.keep = keep
        self._iterable = None
        self._saved_frames = []
        self._iter_parents = ()
        super().__init__(fn)

    def generator(self, iterable):
        self._iterable = iterable
        return self

    def _get_iter_parents(self, frame):
        iter_parents = [frame]
        frontier = (frame,)
        while frontier:
            frontier = sum([p.iter_parents for p in frontier], ())
            iter_parents += frontier
        return tuple(iter_parents)

    def __iter__(self):
        assert self._iterable is not None
        self._iter_parents = self._get_iter_parents(_DIM_STACK.current_frame)
        with ExitStack() as stack:
            for value in self._iterable:
                stack.enter_context(self)
                yield value

    def __enter__(self):
        if self.keep and self._saved_frames:
            saved_frame = self._saved_frames.pop()
            name_to_dim, dim_to_name = saved_frame.name_to_dim, saved_frame.dim_to_name
        else:
            name_to_dim, dim_to_name = OrderedDict(), OrderedDict()

        frame = StackFrame(
            name_to_dim=name_to_dim,
            dim_to_name=dim_to_name,
            parents=tuple(
                reversed(_DIM_STACK._stack[len(_DIM_STACK._stack) - self.history :])
            ),
            iter_parents=tuple(self._iter_parents),
            keep=self.keep,
        )

        _DIM_STACK.push(frame)
        return super().__enter__()

    def __exit__(self, *args, **kwargs):
        if self.keep:
            # don't keep around references to other frames
            old_frame = _DIM_STACK.pop()
            saved_frame = StackFrame(
                name_to_dim=old_frame.name_to_dim,
                dim_to_name=old_frame.dim_to_name,
                parents=(),
                iter_parents=(),
                keep=self.keep,
            )
            self._saved_frames.append(saved_frame)
        else:
            _DIM_STACK.pop()
        return super().__exit__(*args, **kwargs)

class GlobalNamedMessenger(NamedMessenger):
    def __init__(self, fn=None):
        self._saved_globals = ()
        super().__init__(fn)

    def __enter__(self):
        if self._ref_count == 0:
            for name, dim in self._saved_globals:
                _DIM_STACK.global_frame.write(name, dim)
            self._saved_globals = ()
        return super().__enter__()

    def __exit__(self, *args, **kwargs):
        if self._ref_count == 1:
            for name, dim in self._saved_globals:
                _DIM_STACK.global_frame.free(name, dim)
        return super().__exit__(*args, **kwargs)

    def postprocess_message(self, msg):
        if msg["type"] == "to_funsor":
            self._pyro_post_to_funsor(msg)
        elif msg["type"] == "to_data":
            self._pyro_post_to_data(msg)

    def _pyro_post_to_funsor(self, msg):
        if msg["kwargs"]["dim_type"] in (DimType.GLOBAL, DimType.VISIBLE):
            for name in msg["value"].inputs:
                self._saved_globals += (
                    (name, _DIM_STACK.global_frame.name_to_dim[name]),
                )

    def _pyro_post_to_data(self, msg):
        if msg["kwargs"]["dim_type"] in (DimType.GLOBAL, DimType.VISIBLE):
            for name in msg["args"][0].inputs:
                self._saved_globals += (
                    (name, _DIM_STACK.global_frame.name_to_dim[name]),
                )

class BaseEnumMessenger(NamedMessenger):
    """
    Handles first_available_dim management, enum effects should inherit from this
    """

    def __init__(self, fn=None, first_available_dim=None):
        assert (
            first_available_dim is None or first_available_dim < 0
        ), first_available_dim
        self.first_available_dim = first_available_dim
        super().__init__(fn)

    def __enter__(self):
        if self._ref_count == 0 and self.first_available_dim is not None:
            self._prev_first_dim = _DIM_STACK.set_first_available_dim(
                self.first_available_dim
            )
        return super().__enter__()

    def __exit__(self, *args, **kwargs):
        if self._ref_count == 1 and self.first_available_dim is not None:
            _DIM_STACK.set_first_available_dim(self._prev_first_dim)
        return super().__exit__(*args, **kwargs)

##########################################
# User-facing handler implementations
##########################################

class plate(GlobalNamedMessenger):
    """
    An alternative implementation of :class:`numpyro.primitives.plate` primitive. Note
    that only this version is compatible with enumeration.

    There is also a context manager
    :func:`~numpyro.contrib.funsor.infer_util.plate_to_enum_plate`
    which converts `numpyro.plate` statements to this version.

    :param str name: Name of the plate.
    :param int size: Size of the plate.
    :param int subsample_size: Optional argument denoting the size of the mini-batch.
        This can be used to apply a scaling factor by inference algorithms. e.g.
        when computing ELBO using a mini-batch.
    :param int dim: Optional argument to specify which dimension in the tensor
        is used as the plate dim. If `None` (default), the rightmost available dim
        is allocated.
    """

    def __init__(
        self, name, size, subsample_size=None, dim=None, subsample_scale=1.0
    ):  # TODO: check correct subsample_scale!
        self.name = name
        self.size = size
        if dim is not None and dim >= 0:
            raise ValueError("dim arg must be negative.")
        self.dim, indices = OrigPlateMessenger._subsample(
            self.name, self.size, subsample_size, dim
        )
        self.subsample_size = indices.shape[0]
        self.subsample_scale = subsample_scale
        self._indices = funsor.Tensor(
            indices,
            OrderedDict([(self.name, funsor.Bint[self.subsample_size])]),
            self.subsample_size,
        )
        super(plate, self).__init__(None)

    def __enter__(self):
        super().__enter__()  # do this first to take care of globals recycling
        name_to_dim = (
            OrderedDict([(self.name, self.dim)])
            if self.dim is not None
            else OrderedDict()
        )
        indices = to_data(
            self._indices, name_to_dim=name_to_dim, dim_type=DimType.VISIBLE
        )
        # extract the dimension allocated by to_data to match plate's current behavior
        self.dim, self.indices = -len(indices.shape), indices.squeeze()
        return self.indices

    @staticmethod
    def _get_batch_shape(cond_indep_stack):
        n_dims = max(-f.dim for f in cond_indep_stack)
        batch_shape = [1] * n_dims
        for f in cond_indep_stack:
            batch_shape[f.dim] = f.size
        return tuple(batch_shape)

    def process_message(self, msg):
        if msg["type"] in ["to_funsor", "to_data"]:
            return super().process_message(msg)
        return OrigPlateMessenger.process_message(self, msg)

    def postprocess_message(self, msg):
        if msg["type"] in ["to_funsor", "to_data"]:
            return super().postprocess_message(msg)
        # NB: copied literally from original plate messenger, with self._indices is replaced
        # by self.indices
        if msg["type"] in ("subsample", "param") and self.dim is not None:
            event_dim = msg["kwargs"].get("event_dim")
            if event_dim is not None:
                assert event_dim >= 0
                dim = self.dim - event_dim
                shape = jnp.shape(msg["value"])
                if len(shape) >= -dim and shape[dim] != 1:
                    if shape[dim] != self.size:
                        if msg["type"] == "param":
                            statement = "numpyro.param({}, ..., event_dim={})".format(
                                msg["name"], event_dim
                            )
                        else:
                            statement = "numpyro.subsample(..., event_dim={})".format(
                                event_dim
                            )
                        raise ValueError(
                            "Inside numpyro.plate({}, {}, dim={}) invalid shape of {}: {}".format(
                                self.name, self.size, self.dim, statement, shape
                            )
                        )
                    if self.subsample_size < self.size:
                        value = msg["value"]
                        new_value = jnp.take(value, self.indices, dim)
                        msg["value"] = new_value

class enum(BaseEnumMessenger):
    """
    Enumerates in parallel over discrete sample sites marked
    ``infer={"enumerate": "parallel"}``.

    :param callable fn: Python callable with NumPyro primitives.
    :param int first_available_dim: The first tensor dimension (counting
        from the right) that is available for parallel enumeration. This
        dimension and all dimensions left may be used internally by Pyro.
        This should be a negative integer or None.
    """

    def process_message(self, msg):
        if (
            msg["type"] != "sample"
            or msg.get("done", False)
            or msg["is_observed"]
            or msg["infer"].get("expand", False)
            or msg["infer"].get("enumerate") != "parallel"
            or (not msg["fn"].has_enumerate_support)
        ):
            if msg["type"] == "control_flow":
                msg["kwargs"]["enum"] = True
                msg["kwargs"]["first_available_dim"] = self.first_available_dim
            return super().process_message(msg)

        if msg["infer"].get("num_samples", None) is not None:
            raise NotImplementedError("TODO implement multiple sampling")

        if msg["infer"].get("expand", False):
            raise NotImplementedError("expand=True not implemented")

        size = msg["fn"].enumerate_support(expand=False).shape[0]
        raw_value = jnp.arange(0, size)
        funsor_value = funsor.Tensor(
            raw_value, OrderedDict([(msg["name"], funsor.Bint[size])]), size
        )

        msg["value"] = to_data(funsor_value)
        msg["done"] = True

class trace(OrigTraceMessenger):
    """
    This version of :class:`~numpyro.handlers.trace` handler records
    information necessary to do packing after execution.

    Each sample site is annotated with a "dim_to_name" dictionary,
    which can be passed directly to :func:`to_funsor`.
    """

    def postprocess_message(self, msg):
        if msg["type"] == "sample":
            total_batch_shape = lax.broadcast_shapes(
                tuple(msg["fn"].batch_shape),
                jnp.shape(msg["value"])[: jnp.ndim(msg["value"]) - msg["fn"].event_dim],
            )
            msg["infer"]["dim_to_name"] = NamedMessenger._get_dim_to_name(
                total_batch_shape
            )
            msg["infer"]["name_to_dim"] = {
                name: dim for dim, name in msg["infer"]["dim_to_name"].items()
            }
        if msg["type"] in ("sample", "param"):
            super().postprocess_message(msg)

def markov(fn=None, history=1, keep=False):
    """
    Markov dependency declaration.

    This is a statistical equivalent of a memory management arena.

    :param callable fn: Python callable with NumPyro primitives.
    :param int history: The number of previous contexts visible from the
        current context. Defaults to 1. If zero, this is similar to
        :class:`numpyro.primitives.plate`.
    :param bool keep: If true, frames are replayable. This is important
        when branching: if ``keep=True``, neighboring branches at the same
        level can depend on each other; if ``keep=False``, neighboring branches
        are independent (conditioned on their shared ancestors).
    """
    if fn is not None and not callable(fn):  # Used as a generator
        return LocalNamedMessenger(fn=None, history=history, keep=keep).generator(
            iterable=fn
        )
    return LocalNamedMessenger(fn, history=history, keep=keep)

####################
# New primitives
####################

def to_funsor(x, output=None, dim_to_name=None, dim_type=DimType.LOCAL):
    """
    A primitive to convert a Python object to a :class:`~funsor.terms.Funsor`.

    :param x: An object.
    :param funsor.domains.Domain output: An optional output hint to uniquely
        convert a data to a Funsor (e.g. when `x` is a string).
    :param OrderedDict dim_to_name: An optional mapping from negative
        batch dimensions to name strings.
    :param int dim_type: Either 0, 1, or 2. This optional argument indicates
        a dimension should be treated as 'local', 'global', or 'visible',
        which can be used to interact with the global :class:`DimStack`.
    :return: A Funsor equivalent to `x`.
    :rtype: funsor.terms.Funsor
    """
    dim_to_name = OrderedDict() if dim_to_name is None else dim_to_name

    initial_msg = {
        "type": "to_funsor",
        "fn": lambda x, output, dim_to_name, dim_type: funsor.to_funsor(
            x, output=output, dim_to_name=dim_to_name
        ),
        "args": (x,),
        "kwargs": {"output": output, "dim_to_name": dim_to_name, "dim_type": dim_type},
        "value": None,
        "mask": None,
    }

    msg = apply_stack(initial_msg)
    return msg["value"]

def to_data(x, name_to_dim=None, dim_type=DimType.LOCAL):
    """
    A primitive to extract a python object from a :class:`~funsor.terms.Funsor`.

    :param ~funsor.terms.Funsor x: A funsor object
    :param OrderedDict name_to_dim: An optional inputs hint which maps
        dimension names from `x` to dimension positions of the returned value.
    :param int dim_type: Either 0, 1, or 2. This optional argument indicates
        a dimension should be treated as 'local', 'global', or 'visible',
        which can be used to interact with the global :class:`DimStack`.
    :return: A non-funsor equivalent to `x`.
    """
    name_to_dim = OrderedDict() if name_to_dim is None else name_to_dim

    initial_msg = {
        "type": "to_data",
        "fn": lambda x, name_to_dim, dim_type: funsor.to_data(
            x, name_to_dim=name_to_dim
        ),
        "args": (x,),
        "kwargs": {"name_to_dim": name_to_dim, "dim_type": dim_type},
        "value": None,
        "mask": None,
    }

    msg = apply_stack(initial_msg)
    return msg["value"]

max_stars_repo_name: aleatory-science/numpyro

max_issues_repo_name: aleatory-science/numpyro

max_forks_repo_name: aleatory-science/numpyro

reason: from numpy,import jax,from jax


text:
from keras import backend as K
import keras
import cv2
from Utils import *
from sklearn.model_selection import train_test_split
from keras.callbacks import ModelCheckpoint
from keras.optimizers import Adam
from keras.layers import *
from keras.applications.mobilenetv2 import MobileNetV2
import tensorflow as tf
import numpy as np
from keras.models import Model
from keras.models import model_from_json
import matplotlib.pyplot as plt
import os

#Variable Definition
img_w = 512
img_h = 512
channels = 3
classes = 1
info = 5
grid_w = 16
grid_h = 16

def load_model(strr):
    json_file = open(strr, 'r')
    loaded_model_json = json_file.read()
    json_file.close()
    loaded_model = model_from_json(loaded_model_json)
    return loaded_model

model = load_model('model/text_detect_model.json')
model.load_weights('model/text_detect.h5')

def predict_func(model , inp , iou , name):

    ans = model.predict(inp)

    #np.save('Results/ans.npy',ans)
    boxes = decode(ans[0] , img_w , img_h , iou)

    img = ((inp + 1)/2)
    img = img[0]
    #plt.imshow(img)
    #plt.show()

    for i in boxes:

        i = [int(x) for x in i]

        img = cv2.rectangle(img , (i[0] ,i[1]) , (i[2] , i[3]) , color = (0,255,0) , thickness = 2)

    plt.imshow(img)
    plt.show()

    cv2.imwrite(os.path.join('Results' , str(name) + '.jpg') , img*255.0)

# Run evaluation
for i in os.listdir('Test'):
    img = cv2.imread(os.path.join('Test',i))
    img = cv2.resize(img,(512,512))
    img = (img - 127.5)/127.5
    predict_func(model , np.expand_dims(img,axis= 0) , 0.5 , 'sample')

ext: py

max_stars_repo_name: dxlin/tf-yolo-textdetect

max_issues_repo_name: dxlin/tf-yolo-textdetect

max_forks_repo_name: dxlin/tf-yolo-textdetect

reason: import numpy


text:
## original code is from https://github.com/ShiyuLiang/odin-pytorch/blob/master/code/calMetric.py
## Modeified by Kimin Lee
from __future__ import print_function
import torch
from torch.autograd import Variable
import torch.nn as nn
import torch.nn.functional as F
import numpy as np
import torch.optim as optim
import torchvision
import torchvision.transforms as transforms
import numpy as np
import time
from scipy import misc

def tpr95(dir_name):
    #calculate the falsepositive error when tpr is 95%
    cifar = np.loadtxt('%s/confidence_Base_In.txt'%dir_name, delimiter=',')
    other = np.loadtxt('%s/confidence_Base_Out.txt'%dir_name, delimiter=',')
    Y1 = other
    X1 = cifar
    end = np.max([np.max(X1), np.max(Y1)])
    start = np.min([np.min(X1),np.min(Y1)])
    gap = (end- start)/200000 # precision:200000

    total = 0.0
    fpr = 0.0
    thresh = 0.0
    for delta in np.arange(start, end, gap):
        tpr = np.sum(np.sum(X1 >= delta)) / np.float(len(X1))
        error2 = np.sum(np.sum(Y1 > delta)) / np.float(len(Y1))
        if tpr <= 0.9505 and tpr >= 0.9495:
            thresh += delta
            fpr += error2
            total += 1
    if total == 0:
        print('corner case')
        fprBase = 1
    else:
        fprBase = fpr/total
        thresh = thresh/total

    return fprBase, thresh

def auroc(dir_name):
    #calculate the AUROC
    f1 = open('%s/Update_Base_ROC_tpr.txt'%dir_name, 'w')
    f2 = open('%s/Update_Base_ROC_fpr.txt'%dir_name, 'w')

    cifar = np.loadtxt('%s/confidence_Base_In.txt'%dir_name, delimiter=',')
    other = np.loadtxt('%s/confidence_Base_Out.txt'%dir_name, delimiter=',')
    Y1 = other
    X1 = cifar
    end = np.max([np.max(X1), np.max(Y1)])
    start = np.min([np.min(X1),np.min(Y1)])
    gap = (end- start)/200000

    aurocBase = 0.0
    fprTemp = 1.0
    for delta in np.arange(start, end, gap):
        tpr = np.sum(np.sum(X1 >= delta)) / np.float(len(X1))
        fpr = np.sum(np.sum(Y1 > delta)) / np.float(len(Y1))
        f1.write("{}\n".format(tpr))
        f2.write("{}\n".format(fpr))
        aurocBase += (-fpr+fprTemp)*tpr
        fprTemp = fpr

    return aurocBase

def auprIn(dir_name):
    #calculate the AUPR

    cifar = np.loadtxt('%s/confidence_Base_In.txt'%dir_name, delimiter=',')
    other = np.loadtxt('%s/confidence_Base_Out.txt'%dir_name, delimiter=',')
    precisionVec = []
    recallVec = []
    Y1 = other
    X1 = cifar
    end = np.max([np.max(X1), np.max(Y1)])
    start = np.min([np.min(X1),np.min(Y1)])
    gap = (end- start)/200000

    auprBase = 0.0
    recallTemp = 1.0
    for delta in np.arange(start, end, gap):
        tp = np.sum(np.sum(X1 >= delta)) / np.float(len(X1))
        fp = np.sum(np.sum(Y1 >= delta)) / np.float(len(Y1))
        if tp + fp == 0: continue
        precision = tp / (tp + fp)
        recall = tp
        precisionVec.append(precision)
        recallVec.append(recall)
        auprBase += (recallTemp-recall)*precision
        recallTemp = recall
    auprBase += recall * precision

    return auprBase

def auprOut(dir_name):
    #calculate the AUPR
    cifar = np.loadtxt('%s/confidence_Base_In.txt'%dir_name, delimiter=',')
    other = np.loadtxt('%s/confidence_Base_Out.txt'%dir_name, delimiter=',')
    Y1 = other
    X1 = cifar
    end = np.max([np.max(X1), np.max(Y1)])
    start = np.min([np.min(X1),np.min(Y1)])
    gap = (end- start)/200000

    auprBase = 0.0
    recallTemp = 1.0
    for delta in np.arange(end, start, -gap):
        fp = np.sum(np.sum(X1 < delta)) / np.float(len(X1))
        tp = np.sum(np.sum(Y1 < delta)) / np.float(len(Y1))
        if tp + fp == 0: break
        precision = tp / (tp + fp)
        recall = tp
        auprBase += (recallTemp-recall)*precision
        recallTemp = recall
    auprBase += recall * precision

    return auprBase

def detection(dir_name):
    #calculate the minimum detection error

    cifar = np.loadtxt('%s/confidence_Base_In.txt'%dir_name, delimiter=',')
    other = np.loadtxt('%s/confidence_Base_Out.txt'%dir_name, delimiter=',')
    Y1 = other
    X1 = cifar
    end = np.max([np.max(X1), np.max(Y1)])
    start = np.min([np.min(X1),np.min(Y1)])
    gap = (end- start)/200000

    errorBase = 1.0
    for delta in np.arange(start, end, gap):
        tpr = np.sum(np.sum(X1 < delta)) / np.float(len(X1))
        error2 = np.sum(np.sum(Y1 > delta)) / np.float(len(Y1))
        temp = errorBase
        errorBase = np.minimum(errorBase, (tpr+error2)/2.0)
        if temp != errorBase :
            threshold = delta

    return errorBase, threshold

def detection_at_thresh_50percent(dir_name):
    #calculate the minimum detection error

    cifar = np.loadtxt('%s/confidence_Base_In.txt'%dir_name, delimiter=',')
    other = np.loadtxt('%s/confidence_Base_Out.txt'%dir_name, delimiter=',')
    Y1 = other
    X1 = cifar

    tpr = np.sum(np.sum(X1 < 0.5)) / np.float(len(X1))
    error2 = np.sum(np.sum(Y1 > 0.5)) / np.float(len(Y1))
    errorBase = (tpr+error2)/2.0

    return errorBase

import matplotlib.pyplot as plt
def draw_graph(dir_name, acc):
    cifar = np.loadtxt('%s/confidence_Base_In.txt'%dir_name, delimiter=',')
    other = np.loadtxt('%s/confidence_Base_Out.txt'%dir_name, delimiter=',')
    Y1 = other
    X1 = cifar
    bins = np.arange(0,1.005, step=0.01)
    # bins = [0.0, 0.05, 0.1,0.15, 0.2, 0.25,0.3, 0.35,0.4,0.45, 0.5,0.55, 0.6, 0.65, 0.7, 0.75, 0.8,0.85, 0.9,0.95, 1.0]
    plt.hist(X1, bins=bins, alpha=0.5, color='black', label='in', rwidth= 1.4, weights=np.ones(len(X1)) / len(X1))
    plt.hist(Y1, bins=bins, alpha=0.5, color='orange', label='out', rwidth=1.4, weights=np.ones(len(Y1)) / len(Y1))
    plt.ylim([0,1])
    xticks = np.arange(0, 1.005, step=0.1)
    plt.xticks(xticks)
    plt.axvline(X1.mean(), color='black', linestyle='dashed', linewidth=1)
    plt.axvline(Y1.mean(), color='orange', linestyle='dashed', linewidth=1)
    plt.axvline(acc/100., color='red', linestyle='dashed', linewidth=2, label='acc:{:.2f}%'.format(acc))
    plt.legend(loc='upper right')
    plt.savefig('{}/graph_{}acc.png'.format(dir_name,acc))
    plt.close()

    cifar_sharing_node_correct = np.loadtxt('%s/confidence_Base_In_sharing_node_of_correct_case.txt'%dir_name, delimiter=',')
    cifar_sharing_node__wrong = np.loadtxt('%s/confidence_Base_In_sharing_node_of_wrong_case.txt'%dir_name, delimiter=',')
    other_sharing_node = np.loadtxt('%s/confidence_Base_Out_sharing_node.txt'%dir_name, delimiter=',')
    Y1_sharing = other_sharing_node
    X1_sharing = cifar_sharing_node_correct
    X2_sharing = cifar_sharing_node__wrong
    bins = np.arange(0,1.005, step=0.02)
    # bins = [0.0, 0.05, 0.1,0.15, 0.2, 0.25,0.3, 0.35,0.4,0.45, 0.5,0.55, 0.6, 0.65, 0.7, 0.75, 0.8,0.85, 0.9,0.95, 1.0]
    plt.hist(X1_sharing, bins=bins, alpha=0.3, color='black', label='in_right_sharing', rwidth= 0.7, weights=np.ones(len(X1_sharing)) / len(X1_sharing))
    plt.hist(X2_sharing, bins=bins, alpha=0.3, color='green', label='in_wrong_sharing', rwidth= 0.7, weights=np.ones(len(X2_sharing)) / len(X2_sharing))
    plt.hist(Y1_sharing, bins=bins, alpha=0.3, color='orange', label='out_sharing', rwidth=0.7, weights=np.ones(len(Y1_sharing)) / len(Y1_sharing))
    plt.ylim([0,1])
    plt.xticks(xticks)
    plt.legend(loc='upper right')
    plt.savefig('{}/graph_{}acc_OOD_node_split.png'.format(dir_name,acc))
    plt.close()

    cifar_sharing_node = np.loadtxt('%s/confidence_Base_In_sharing_node.txt'%dir_name, delimiter=',')
    other_sharing_node = np.loadtxt('%s/confidence_Base_Out_sharing_node.txt'%dir_name, delimiter=',')
    Y1_sharing = other_sharing_node
    X1_sharing = cifar_sharing_node
    # bins = [0.0, 0.05, 0.1,0.15, 0.2, 0.25,0.3, 0.35,0.4,0.45, 0.5,0.55, 0.6, 0.65, 0.7, 0.75, 0.8,0.85, 0.9,0.95, 1.0]
    # bins = [0.00, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.10, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.20]
    bins = np.arange(0,0.2, step=0.0033)
    # bins = [0.300, 0.325, 0.350, 0.375, 0.400, 0.425, 0.450, 0.475, 0.500, 0.525]
    plt.hist(X1_sharing, bins=bins, alpha=0.3, color='black', label='in_OOD_node', rwidth= 0.7, weights=np.ones(len(X1_sharing)) / len(X1_sharing))
    plt.hist(Y1_sharing, bins=bins, alpha=0.3, color='orange', label='out_OOD_node', rwidth=0.7, weights=np.ones(len(Y1_sharing)) / len(Y1_sharing))
    plt.ylim([0,1])
    plt.legend(loc='upper right')
    plt.savefig('{}/graph_{}acc_OOD_node.png'.format(dir_name,acc))
    plt.close()

    cifar_sharing_node = np.loadtxt('%s/confidence_Base_In_softmax_max_and_sharing_node.txt'%dir_name, delimiter=',')
    other_sharing_node = np.loadtxt('%s/confidence_Base_Out_softmax_max_and_sharing_node.txt'%dir_name, delimiter=',')
    X1_sharing_soft = cifar_sharing_node
    Y1_sharing_soft = other_sharing_node
    bins = np.arange(0,0.6, step=0.01)
    # bins = [0.0, 0.05, 0.1,0.15, 0.2, 0.25,0.3, 0.35,0.4,0.45, 0.5,0.55, 0.6, 0.65, 0.7, 0.75, 0.8,0.85, 0.9,0.95, 1.0]
    # bins = [0.00, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.10, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.20]
    # bins = [0.300, 0.325, 0.350, 0.375, 0.400, 0.425, 0.450, 0.475, 0.500, 0.525]
    plt.hist(X1_sharing_soft, bins=bins, alpha=0.3, color='black', label='in_soft(max,OOD_node)', rwidth= 0.7, weights=np.ones(len(X1_sharing)) / len(X1_sharing))
    plt.hist(Y1_sharing_soft, bins=bins, alpha=0.3, color='red', label='out_soft(max,OOD_node)', rwidth=0.7, weights=np.ones(len(Y1_sharing)) / len(Y1_sharing))
    plt.ylim([0,1])
    plt.legend(loc='upper right')
    plt.savefig('{}/graph_{}acc_soft_max_and_OOD_node.png'.format(dir_name,acc))
    plt.close()

    cifar_correct = np.loadtxt('%s/confidence_Base_In_correct.txt'%dir_name, delimiter=',')
    cifar_wrong = np.loadtxt('%s/confidence_Base_In_wrong.txt'%dir_name, delimiter=',')
    other = np.loadtxt('%s/confidence_Base_Out.txt'%dir_name, delimiter=',')
    Y1 = other
    X1_correcdt = cifar_correct
    X1_wrong = cifar_wrong
    # bins = [0.0, 0.05, 0.1,0.15, 0.2, 0.25,0.3, 0.35,0.4,0.45, 0.5,0.55, 0.6, 0.65, 0.7, 0.75, 0.8,0.85, 0.9,0.95, 1.0]
    bins = np.arange(0,1.005, step=0.02)
    plt.hist(X1_correcdt, bins=bins, alpha=0.3, color='black', label='in_right', rwidth= 0.7, weights=np.ones(len(X1_correcdt)) / len(X1_correcdt))
    plt.hist(X1_wrong, bins=bins, alpha=0.3, color='green', label='in_wrong', rwidth= 0.7, weights=np.ones(len(X1_wrong)) / len(X1_wrong))
    plt.hist(Y1, bins=bins, alpha=0.3, color='orange', label='out', rwidth=0.7, weights=np.ones(len(Y1)) / len(Y1))
    plt.ylim([0,1])
    plt.xticks(xticks)
    plt.axvline(X1_correcdt.mean(), color='black', linestyle='dashed', linewidth=1)
    plt.axvline(X1_wrong.mean(), color='green', linestyle='dashed', linewidth=1)
    plt.axvline(Y1.mean(), color='orange', linestyle='dashed', linewidth=1)
    plt.axvline(acc/100., color='red', linestyle='dashed', linewidth=2, label='acc:{:.1f}%'.format(acc))
    plt.legend(loc='upper right')
    plt.savefig('{}/graph_{}acc_right_wrong_split.png'.format(dir_name,acc))
    plt.close()

def metric(dir_name, acc=None, mode='all'):
    draw_graph(dir_name, acc)
    if mode == 'graph':
        print('graph draw_finish -> exit')
        exit()
    print("{:>34}".format("Performance of Baseline detector"))
    fprBase, thresh = tpr95(dir_name)
    print("{:20}{:13.3f}%".format("TNR at TPR 95%:",(1-fprBase)*100))
    print('threshold at TPR 95% : {}'.format(thresh))

    errorBase, threshold = detection(dir_name)
    print("{:20}{:13.3f}%".format("Detection acc:",(1-errorBase)*100))
    print('threshold at min(detection error) : {}'.format(threshold))
    errorBase_at_thresh_50percent = detection_at_thresh_50percent(dir_name)
    print('Detection acc at 0.5 threshold : {}'.format((1-errorBase_at_thresh_50percent) * 100.))

    aurocBase = auroc(dir_name)
    print("{:20}{:13.3f}%".format("AUROC:",aurocBase*100))

    auprinBase = auprIn(dir_name)
    print("{:20}{:13.3f}%".format("AUPR In:",auprinBase*100))
    auproutBase = auprOut(dir_name)
    print("{:20}{:13.3f}%".format("AUPR Out:",auproutBase*100))

    with open('{}/result.txt'.format(dir_name), 'w') as f :
        f.write("{:>34}\n".format("Performance of Baseline detector"))
        f.write("{:20}{:13.3f}%\n".format("TNR at TPR 95%:",(1-fprBase)*100))
        f.write("{:20}{:13.3f}\n".format("Threshold at TPR 95%:", thresh))
        f.write("{:20}{:13.3f}%\n".format("Detection acc:",(1-errorBase)*100))
        f.write("{:20}{:13.3f}\n".format("Threshold at Detectir acc.", threshold))
        f.write("{:20}{:13.3f}\n".format("Detectir acc at 0.5 thresh.", (1-errorBase_at_thresh_50percent)*100.))
        f.write("{:20}{:13.3f}%\n".format("AUROC:",aurocBase*100))
        f.write("{:20}{:13.3f}%\n".format("AUPR In:",auprinBase*100))
        f.write("{:20}{:13.3f}%\n".format("AUPR Out:",auproutBase*100))

reason: import numpy,from scipy


text:
import numpy as np

class SimLC(object):
    """Class for simulating light curves"""

    def __init__(self, seed=None):
        """Create a simlc object, setting some defaults

        Parameters:
            seed: a seed for the random number generator

        """

        # built-in models #
        self.builtin_models = ['powerlaw', 'broken_powerlaw', 'bending_powerlaw',
                                'step', 'constant', 'lorentz', 'user_array']

        # model containers #
        self.psd_models = []
        self.lag_models = []

        # seed random generator #
        self.rng = np.random.RandomState(seed)

    def add_model(self, model, params, clear=True, lag=False):
        """Add a model to the generating psd/lag models. Adds the model
            to self.psd_models or self.lag_models depending on lag

        Parameters:
            model: a callable model(freq, params), or a string
                from the builtin models
            params: model parameters
            clear: If true, clear all previously defined models
            lag: if true, this is a lag model, else it is a psd
        """

        # check the model is callable or in self.builtin_models #
        if isinstance(model, str):
            if not model in self.builtin_models:
                mtxt = ', '.join([x for x in self.builtin_models])
                raise ValueError('model is not in builtin_models\n' + mtxt)
            else:
                model = eval('SimLC.' + model)
        else:
            if not hasattr(model, '__call__'):
                raise ValueError('model is not callable')

        # add given model to the generating list #
        m_container = self.lag_models if lag else self.psd_models

        if clear:
            for k in m_container: m_container.remove(k)
        m_container.append([model, params])

    def calculate_model(self, freq, lag=False):
        """Calculate the psd/lag model using the models
            added with @add_model.

        Parameters:
            freq: frequency array.
            lag: if true, calculate the lag model, else do the psd

        Returns:
            model: array of the same length as freq, containing the model

        Note:
            the normalization of the returned psd/lag is taken from
            the input model parameters without any renormalization

        """

        # check we have models added #
        models = self.lag_models if lag else self.psd_models
        if len(models) == 0:
            raise ValueError('No models added. Use add_model to do so.')

        # calculate the model #
        freq = np.array(freq, np.float)
        mod = np.zeros_like(freq)
        for m,p in models:
            mod += m(freq, p)

        return mod

    def simulate(self, n, dt, mu, norm='var'):
        """Simulate a light curve using the psd model stored in
            self.psd_models, added with add_model

        Parameters:
            n: number of points in the light curve
            dt: time sampling
            mu: the light curve mean
            norm: string var|rms|leahy

        Returns:
            Nothing. The simulated light curve and time are stored
            in self.t and self.x

        Note:
            The normalized psd is stored in self.normalized_psd

        """

        # make sure the norm is known #
        if not (norm in ['var', 'leahy', 'rms']):
            raise ValueError('norm need to be one of var|rms|leahy')

        # calculate psd #
        freq = np.fft.rfftfreq(n, dt)
        psd = self.calculate_model(freq)
        self.normalized_psd = np.array([freq, psd])

        # get correct renoramlization #
        expon = {'var':0, 'leahy':1, 'rms':2}
        renorm = mu**expon[norm] * n/(2.*dt)
        psd *= renorm

        # inverse fft #
        ix = ( self.rng.randn(len(psd)) * np.sqrt(0.5*psd)
             + self.rng.randn(len(psd)) * np.sqrt(0.5*psd)*1j )
        ix[0] = n * mu
        ix[-1] = np.abs(ix[-1])
        self.x = np.fft.irfft(ix)
        self.t = np.arange(n) * dt * 1.0

    def simulate_pdf(self, n, dt, mu, norm='var', pdf='lognorm(s=0.5)'):
        """Simulate a light curve using the psd model stored in
            self.psd_models, added with add_model, and enforcing
            log-normal distribution.
            This uses the algorithm of Emmanoulopoulos+ (2013) MNRAS 433, 907–927

        Parameters:
            n: number of points in the light curve
            dt: time sampling
            mu: the light curve mean
            norm: string var|rms|leahy
            pdf: a string representing a probability distribution
                from scipy.stats. e.g. lognorm(s=0.5)

        Returns:
            Nothing. The simulated light curve and time are stored
            in self.t and self.x

        Note:
            The normalized psd is stored in self.normalized_psd

        """

        # make sure the norm is known #
        if not (norm in ['var', 'leahy', 'rms']):
            raise ValueError('norm need to be one of var|rms|leahy')

        import scipy.stats as st
        lognorm = eval('st.{}'.format(pdf))

        # calculate psd #
        freq = np.fft.rfftfreq(n, dt)
        psd = self.calculate_model(freq)
        self.normalized_psd = np.array([freq, psd])

        # get correct renoramlization #
        expon = {'var':0, 'leahy':1, 'rms':2}
        renorm = mu**expon[norm] * n/(2.*dt)
        psd *= renorm

        # do the algorithm #
        ix = ( self.rng.randn(len(psd)) * np.sqrt(0.5*psd)
             + self.rng.randn(len(psd)) * np.sqrt(0.5*psd)*1j )
        ix[0] = n * mu
        ix[-1] = np.abs(ix[-1])
        # step-1 #
        anorm  = np.abs(ix)

        # step-2; before loop starts #
        xsim   = lognorm.rvs(n)
        diff   = 1.0
        while diff > 1e-4:

            # step-2; main loop #
            psim   = np.angle(np.fft.rfft(xsim))

            # step-3
            xsim_j = np.fft.irfft(anorm * np.exp(1j*psim))

            # step-4
            xsim_j[np.argsort(xsim_j)] = xsim[np.argsort(xsim)]

            # check for convergenece & prepare for next loop #
            diff   = np.sum((xsim - xsim_j)**2)/n
            xsim[:]= xsim_j[:]

        # get the final light curve
        psim   = np.angle(np.fft.rfft(xsim))
        self.x = np.fft.irfft(anorm * np.exp(1j*psim))
        self.t = np.arange(n) * dt * 1.0

    def apply_lag(self, phase=False):
        """Apply lag in self.lag_models to the simulate self.c

        Parameters:
            phase: the lags in self.lag_models are in radians
            if true, else in seconds.

        Returns:
            Nothing. creates self.y for the shifted light curve
            and self.normalized_lag for the actual [freq, lag ] used

        Note:
            The lag vs freq is found by calling the functions in
            self.lag_models, filled by calling add_model(..., lag=True)
            The light curve in self.x is assumed to exist (i.e. self.simulate
            should've been called already)

        """

        # has simulate beed run? #
        try:
            self.x
        except:
            raise ValueError('self.x does not exists. Run simulate first')

        freq = self.normalized_psd[0]
        lag  = self.calculate_model(freq, lag=True)
        self.normalized_lag = np.array([freq, lag])
        self.y = SimLC.lag_array(self.x, lag, phase, freq)

    @staticmethod
    def lag_array(x, lag, phase=False, freq=None):
        """Shift the x by amount lag

        Args:
            x: light curve to be shifted
            lag: float or array of length len(n)/2+1
            phase: if True, lag is in radians. Otherwise in seconds
            freq: the frequency axis (used when phase is False)

        Returns:
            An array containing the shifted light curve.

        """
        nfq = len(x)/2 + 1
        if not isinstance(lag, np.ndarray): lag = np.repeat(lag, nfq)
        if nfq != len(lag):
            raise ValueError('Lag array does not match given array')

        if freq is not None and nfq != len(freq):
            raise ValueError('freq array does not match given array')

        phi = lag if phase else lag*2.*np.pi*freq
        xfft = np.fft.rfft(x)
        xfft[1:-1] *= np.exp(-1j*phi[1:-1])
        return np.fft.irfft(xfft)

    @staticmethod
    def add_noise(x, norm=None, seed=None, dt=1.0):
        """Add noise to lcurve x

        Parameters:
            norm: if None, add Poisson noise, else
                gaussian noise, with std=norm
            seed: random seed
            dt: used with norm is None. It gives the time samling
                of the light curve. Poisson noise is applied to
                the counts per bin. x in this case is the count rate.
                Counts/bin = Rate/sec * dt

        Returns:
            array similar to x with noise added

        """

        if seed is not None:
            np.random.seed(seed)

        if norm is None:
            xn = np.random.poisson(x*dt)/dt
        else:
            xn = np.random.randn(len(x)) * norm + x
        return xn

    @staticmethod
    def user_array(freq, params):
        """The model is given by the user directly as an array params

        Parameters:
            freq (np.ndarray): the frequency array
            params (np.ndarray): the model values

        Returns:
            array mod of same length as freq, containing the psd/lag model
        """

        if len(params) != len(freq):
            raise ValueError('params does not match freq: %d vs %d'%(len(params), len(freq)))
        return np.array(params)

    @staticmethod
    def powerlaw(freq, params):
        """A powerlaw model for the psd/lag

        Parameters:
            freq (np.ndarray): the frequency array
            params (list or array: [norm, indx]): parameters of the model

        Returns:
            array mod of same length as freq, containing the psd/lag model

        """

        if len(params) != 2:
            raise ValueError('powerlaw needs 2 params: norm, index')
        norm, indx = params
        freq = np.clip(freq, 1e-30, np.inf)
        mod = norm * freq**indx
        return mod

    @staticmethod
    def broken_powerlaw(freq, params):
        """A borken powerlaw model for the psd/lag.

        Parameters:
            freq (np.ndarray): the frequency array
            params (list or array: [norm, indx1, indx2, break]):
                parameters of the model

        Returns:
            array mod of same length as freq, containing the psd/lag model

        """

        if len(params) != 4:
            raise ValueError(('broek_powerlaw needs 4 params: norm, '
                      'index1, index2, break'))
        norm, a1, a2, brk = params
        freq = np.clip(freq, 1e-30, np.inf)
        fq_indx = freq < brk
        mod = np.zeros_like(freq)
        mod[fq_indx] = (freq[fq_indx]/brk) ** a1
        mod[~fq_indx] = (freq[~fq_indx]/brk) ** a2
        mod *= norm * brk**a2
        return mod

    @staticmethod
    def bending_powerlaw(freq, params):
        """A bending powerlaw model for the psd/lag.

        Parameters:
            freq (np.ndarray): the frequency array
            params (list or array: [norm, index, break]):
                parameters of the model. The index below the
                break is always 0

        Returns:
            array mod of same length as freq, containing the psd/lag model

        """

        if len(params) != 3:
            raise ValueError(('bending_powerlaw needs 3 params: norm, '
                      'index, break'))
        norm, a, brk = params
        freq = np.clip(freq, 1e-30, np.inf)
        mod = (norm/freq) * (1+(freq/brk)**(-a-1))**(-1)
        return mod

    @staticmethod
    def step(freq, params):
        """A step function model for the psd/lag.

        Parameters:
            freq (np.ndarray): the frequency array
            params: a list (or array) of 2 lists (arrays).
                The first holds the frequency bin boundaries.
                The second holds the values of psd|lag
                len(list_1) = len(list_2) + 1

        Returns:
            array mod of same length as freq, containing the psd/lag model

        """

        if len(params) != 2:
            raise ValueError('step needs a list of 2 arrays/lists')
        fqL, pars = params
        fqL, pars = np.array(fqL), np.array(pars)
        if len(fqL) != len(pars)+1:
            raise ValueError('fqL and pars do not match')

        ibins = np.digitize(freq, fqL, right=True) - 1
        ibins[ibins==-1] = 0
        ibins[ibins==len(fqL)-1] = len(fqL)-2
        mod = pars[ibins]
        return mod

    @staticmethod
    def constant(freq, params):
        """A constant model for the psd/lag

        Parameters:
            freq (np.ndarray): the frequency array
            params (float): the value of the constant

        Returns:
            array mod of same length as freq, containing the psd model

        """
        mod = freq*0 + params
        return mod

    @staticmethod
    def lorentz(freq, params):
        """A lorentzian model for the psd/lag

        Parameters:
            freq (np.ndarray): the frequency array
            params (list or array: [norm, fq_center, fq_sigma]):
                parameters of the model

        Returns:
            array mod of same length as freq, containing the psd/lag model

        """
        if len(params) != 3:
            raise ValueError('lorentz needs 3 pars: norm, fq_cent, fq_sig')
        norm, fq_center, fq_sigma = params
        mod = norm * (fq_sigma/(2*np.pi)) / (
                    (freq-fq_center)**2 + (fq_sigma/2)**2 )
        return mod

max_stars_repo_name: abduzoghbi/aztools

max_issues_repo_name: abduzoghbi/aztools

max_forks_repo_name: abduzoghbi/aztools

reason: import numpy,import scipy,from scipy


text:
import numpy as np
import pylab as plt
from matplotlib.colors import LinearSegmentedColormap, colorConverter
from src.BAGLE import model, model_fitter, munge
import dynesty
from dynesty import utils as dyutil
from dynesty import plotting as dyplot
from astropy.table import Table
from astropy.time import Time
from astropy import units
from astropy.coordinates import SkyCoord
from multiprocessing import Pool, cpu_count
import time
import pickle
import pdb
import os, shutil
import yaml
from scipy import spatial
from scipy.ndimage import gaussian_filter as norm_kde
from scipy.stats import gaussian_kde
import glob
import astropy.constants as const

def get_data_and_fitter(mnest_base):
    info_file = open(mnest_base + 'params.yaml', 'r')
    info = yaml.full_load(info_file)

    pho_dsets = info['phot_data']
    if 'astrom_data' in info.keys():
        ast_dsets = info['astrom_data']
    else:
        ast_dsets = []

    my_model = getattr(model, info['model'])
    my_data = munge.getdata2(info['target'].lower(),
                             phot_data=pho_dsets,
                             ast_data=ast_dsets)

    # Need this for backwards compatability.
    if 'use_phot_optional_params' in info:
        use_phot_opt = info['use_phot_optional_params']
    else:
        use_phot_opt = False

    # Load up the first fitter object to get the parameter names.
    fitter = model_fitter.PSPL_Solver(my_data, my_model,
                                      add_error_on_photometry = info['add_error_on_photometry'],
                                      multiply_error_on_photometry = info['multiply_error_on_photometry'],
                                      use_phot_optional_params=info['use_phot_optional_params'],
                                      outputfiles_basename = mnest_base)

    return fitter, my_data

def check_priors(fitter, target):
    # This is specifically for photometry + astrometry.

    names = fitter.fitter_param_names
    all_names = fitter.all_param_names

    #####
    # Priors
    #####
    # Number of draws from posterior.
    Nsamp = 5000

    # Store the prior distributions for each parameter.
    priors_dict = {}

    # First, insert all the independent priors.
    for i in range(len(names)):
        priors_dict[names[i]] = fitter.priors[names[i]].rvs(Nsamp)

    # Next, calculate all the derived priors from the independent priors.
    if 'log10_thetaE' in names:
        priors_dict['thetaE_amp'] = 10**priors_dict['log10_thetaE']
    else:
        priors_dict['thetaE_amp'] = priors_dict['thetaE']
    priors_dict['piE'] = np.hypot(priors_dict['piE_E'], priors_dict['piE_N'])
    priors_dict['piRel'] = priors_dict['piE'] * priors_dict['thetaE_amp']
    kappa_tmp = 4.0 * const.G / (const.c ** 2 * units.AU)
    kappa = kappa_tmp.to(units.mas / units.Msun,
                         equivalencies=units.dimensionless_angles()).value
    priors_dict['mL'] = priors_dict['thetaE_amp']**2 / (priors_dict['piRel'] * kappa)

    priors_dict['piL'] = priors_dict['piRel'] + priors_dict['piS']

    priors_dict['muRel'] = priors_dict['thetaE_amp']/(priors_dict['tE']/365.25)
    priors_dict['muRel_E'] = priors_dict['muRel'] * (priors_dict['piE_E']/priors_dict['piE'])
    priors_dict['muRel_N'] = priors_dict['muRel'] * (priors_dict['piE_N']/priors_dict['piE'])
    priors_dict['muL_E'] = priors_dict['muS_E'] - priors_dict['muRel_E']
    priors_dict['muL_N'] = priors_dict['muS_N'] - priors_dict['muRel_N']

    #####
    # Posteriors
    #####
    res = fitter.load_mnest_results(remake_fits=False)

    #####
    # Plot posteriors and priors together for comparison.
    #####
    fig, ax = plt.subplots(nrows=6, ncols=5, figsize=(16, 16))
    fax = ax.flatten()
    fig.subplots_adjust(top = 0.95, bottom = 0.05, left=0.05, right=0.95, hspace = 0.5, wspace = 0.3)
    plt.suptitle(target)
    for ii, name in enumerate(all_names):
        # Get the 3-sigma upper and lower bounds for the priors and the posteriors,
        # so we don't have crazy tails that make the x-range of the plot unreadable.
        sig3_hi = 0.9973
        sig3_lo = 1.0 - sig3_hi

        quantiles = [sig3_lo, sig3_hi]

        quants_prior = model_fitter.weighted_quantile(priors_dict[name], quantiles)
        quants_post = model_fitter.weighted_quantile(res[name], quantiles,
                                                     sample_weight=res['weights'])

        _min = min([quants_prior[0], quants_post[0]])
        _max = max([quants_prior[1], quants_post[1]])

        bins = np.linspace(_min, _max, 30)

        fax[ii].hist(priors_dict[name], bins=bins,
                    density=True, histtype='step', lw=2, color='blue',
                    label='Prior')
        fax[ii].hist(res[name], weights=res['weights'], bins=bins,
                     density=True, histtype='step', lw=2, color='red',
                     label='Posterior')
        fax[ii].yaxis.set_ticks([])
        fax[ii].yaxis.set_ticklabels([])
        fax[ii].set_ylabel(name)

    fax[0].legend()

    plt.savefig(target + '_posteriors_priors.png')

def check_priors_phot(fitter, target):
    # This is specifically for photometry + astrometry.

    names = fitter.fitter_param_names
    all_names = fitter.all_param_names

    #####
    # Priors
    #####
    # Number of draws from posterior.
    Nsamp = 5000

    # Store the prior distributions for each parameter.
    priors_dict = {}

    # First, insert all the independent priors.
    for i in range(len(names)):
        priors_dict[names[i]] = fitter.priors[names[i]].rvs(Nsamp)

    #####
    # Posteriors
    #####
    res = fitter.load_mnest_results(remake_fits=False)

    #####
    # Plot posteriors and priors together for comparison.
    #####
    fig, ax = plt.subplots(nrows=6, ncols=5, figsize=(16, 16))
    fax = ax.flatten()
    fig.subplots_adjust(top = 0.95, bottom = 0.05, left=0.05, right=0.95, hspace = 0.5, wspace = 0.3)
    plt.suptitle(target)
    for ii, name in enumerate(names):
        # Get the 3-sigma upper and lower bounds for the priors and the posteriors,
        # so we don't have crazy tails that make the x-range of the plot unreadable.
        sig3_hi = 0.9973
        sig3_lo = 1.0 - sig3_hi

        quantiles = [sig3_lo, sig3_hi]

        quants_prior = model_fitter.weighted_quantile(priors_dict[name], quantiles)
        quants_post = model_fitter.weighted_quantile(res[name], quantiles,
                                                     sample_weight=res['weights'])

        _min = min([quants_prior[0], quants_post[0]])
        _max = max([quants_prior[1], quants_post[1]])

        bins = np.linspace(_min, _max, 30)

        fax[ii].hist(priors_dict[name], bins=bins,
                    density=True, histtype='step', lw=2, color='blue',
                    label='Prior')
        fax[ii].hist(res[name], weights=res['weights'], bins=bins,
                     density=True, histtype='step', lw=2, color='red',
                     label='Posterior')
        fax[ii].yaxis.set_ticks([])
        fax[ii].yaxis.set_ticklabels([])
        fax[ii].set_ylabel(name)

    fax[0].legend()

    plt.savefig(target + '_posteriors_priors.png')

def check_priors_PSPL_Astrom_Par_Param4(fitter, params_from_post, target, posterior=False):
    """
    Param2, Param4
    """
    Nsamp = 5000
    param_dict = {}

    names = fitter.fitter_param_names
    all_names = fitter.all_param_names
    fig, ax = plt.subplots(nrows=6, ncols=3, figsize=(12, 12))
    fax = ax.flatten()
    fig.subplots_adjust(top = 0.95, bottom = 0.05, left=0.05, right=0.95, hspace = 0.5, wspace = 0.3)
    plt.suptitle(target)
    for i in range(len(names)):
        if names[i] in params_from_post:
            continue
        else:
            param_dict[names[i]] = fitter.priors[names[i]].rvs(Nsamp)
            fax[i].hist(fitter.priors[names[i]].rvs(Nsamp), bins=30,
                       density=True, histtype='step', lw=2, color='blue')
            fax[i].set_ylabel(names[i])

    # Posterior --> Prior
    binmids = []
    for bb in np.arange(len(fitter.post_param_bins)):
        binmids.append((fitter.post_param_bins[bb][:-1] + fitter.post_param_bins[bb][1:])/2)

    draws = np.zeros((len(params_from_post), Nsamp))
    for ii in np.arange(Nsamp):
        draws[:,ii] = fitter.sample_post(binmids, fitter.post_param_cdf, fitter.post_param_bininds)

    for i in np.arange(len(params_from_post)):
        param_dict[params_from_post[i]] = draws[i]
        idx = names.index(params_from_post[i])
        fax[idx].hist(draws[i], bins=30,
                     density=True, histtype='step', lw=2, color='blue',
                     label='Prior')
        fax[idx].set_ylabel(params_from_post[i])

    # Calculate piL, muS, muRel also
    piE = np.hypot(param_dict['piE_E'], param_dict['piE_N'])
    piRel = piE * param_dict['thetaE']
    kappa_tmp = 4.0 * const.G / (const.c ** 2 * units.AU)
    kappa = kappa_tmp.to(units.mas / units.Msun,
                         equivalencies=units.dimensionless_angles()).value
    mL = param_dict['thetaE'] ** 2 / (piRel * kappa)

    piL = piRel + param_dict['piS']

    muRel = param_dict['thetaE']/(param_dict['tE']/365.25)
    muRel_E = muRel * (param_dict['piE_E']/piE)
    muRel_N = muRel * (param_dict['piE_N']/piE)
    muL_E = param_dict['muS_E'] - muRel_E
    muL_N = param_dict['muS_N'] - muRel_N

    fax[len(names)].hist(mL, bins=30,
                        density=True, histtype='step', lw=2, color='blue')
    fax[len(names)].set_ylabel('mL')

    fax[len(names)+1].hist(piL, bins=30,
                          density=True, histtype='step', lw=2, color='blue')
    fax[len(names)+1].set_ylabel('piL')

    fax[len(names)+2].hist(piRel, bins=30,
                          density=True, histtype='step', lw=2, color='blue')
    fax[len(names)+2].set_ylabel('piRel')

    fax[len(names)+3].hist(muL_E, bins=30,
                          density=True, histtype='step', lw=2, color='blue')
    fax[len(names)+3].set_ylabel('muL_E')

    fax[len(names)+4].hist(muL_N, bins=30,
                          density=True, histtype='step', lw=2, color='blue')
    fax[len(names)+4].set_ylabel('muL_N')

    fax[len(names)+5].hist(muRel_E, bins=30,
                          density=True, histtype='step', lw=2, color='blue')
    fax[len(names)+5].set_ylabel('muRel_E')

    fax[len(names)+6].hist(muRel_N, bins=30,
                          density=True, histtype='step', lw=2, color='blue')
    fax[len(names)+6].set_ylabel('muRel_N')

    if posterior:
        res = fitter.load_mnest_results(remake_fits=False)
        for ii, name in enumerate(all_names):
            fax[ii].hist(res[name], weights=res['weights'], bins=30,
                        density=True, histtype='step', lw=2, color='red',
                        label='Posterior')
            fax[ii].yaxis.set_ticks([])
            fax[ii].yaxis.set_ticklabels([])
    fax[0].legend()

#    ax[len(names)].hist(piE, bins=30)
#    ax[len(names)].set_ylabel('piE')

    if posterior:
        plt.savefig(target + '_posteriors_priors.png')
    else:
        plt.savefig(target + '_priors.png')

def check_priors_PSPL_Astrom_Par_Param5(fitter, params_from_post, target, posterior=False):
    """
    Param2, Param4
    """
    Nsamp = 5000
    param_dict = {}

    names = fitter.fitter_param_names
    all_names = fitter.all_param_names
    fig, ax = plt.subplots(nrows=6, ncols=3, figsize=(12, 12))
    fax = ax.flatten()
    fig.subplots_adjust(top = 0.95, bottom = 0.05, left=0.05, right=0.95, hspace = 0.5, wspace = 0.3)
    plt.suptitle(target)
    for i in range(len(names)):
        if names[i] in params_from_post:
            continue
        else:
            param_dict[names[i]] = fitter.priors[names[i]].rvs(Nsamp)
            fax[i].hist(fitter.priors[names[i]].rvs(Nsamp), bins=30,
                       density=True, histtype='step', lw=2, color='blue')
            fax[i].set_ylabel(names[i])

    # Posterior --> Prior
    binmids = []
    for bb in np.arange(len(fitter.post_param_bins)):
        binmids.append((fitter.post_param_bins[bb][:-1] + fitter.post_param_bins[bb][1:])/2)

    draws = np.zeros((len(params_from_post), Nsamp))
    for ii in np.arange(Nsamp):
        draws[:,ii] = fitter.sample_post(binmids, fitter.post_param_cdf, fitter.post_param_bininds)

    for i in np.arange(len(params_from_post)):
        param_dict[params_from_post[i]] = draws[i]
        idx = names.index(params_from_post[i])
        fax[idx].hist(draws[i], bins=30,
                     density=True, histtype='step', lw=2, color='blue',
                     label='Prior')
        fax[idx].set_ylabel(params_from_post[i])

    # Calculate piL, muS, muRel also
    piE = np.hypot(param_dict['piE_E'], param_dict['piE_N'])
    piRel = piE * 10**param_dict['log10_thetaE']
    kappa_tmp = 4.0 * const.G / (const.c ** 2 * units.AU)
    kappa = kappa_tmp.to(units.mas / units.Msun,
                         equivalencies=units.dimensionless_angles()).value
    mL = (10**param_dict['log10_thetaE']) ** 2 / (piRel * kappa)

    piL = piRel + param_dict['piS']

    muRel = 10**param_dict['log10_thetaE']/(param_dict['tE']/365.25)
    muRel_E = muRel * (param_dict['piE_E']/piE)
    muRel_N = muRel * (param_dict['piE_N']/piE)
    muL_E = param_dict['muS_E'] - muRel_E
    muL_N = param_dict['muS_N'] - muRel_N

    fax[len(names)].hist(mL, bins=30,
                        density=True, histtype='step', lw=2, color='blue')
    fax[len(names)].set_ylabel('mL')

    fax[len(names)+1].hist(piL, bins=30,
                          density=True, histtype='step', lw=2, color='blue')
    fax[len(names)+1].set_ylabel('piL')

    fax[len(names)+2].hist(piRel, bins=30,
                          density=True, histtype='step', lw=2, color='blue')
    fax[len(names)+2].set_ylabel('piRel')

    fax[len(names)+3].hist(muL_E, bins=30,
                          density=True, histtype='step', lw=2, color='blue')
    fax[len(names)+3].set_ylabel('muL_E')

    fax[len(names)+4].hist(muL_N, bins=30,
                          density=True, histtype='step', lw=2, color='blue')
    fax[len(names)+4].set_ylabel('muL_N')

    fax[len(names)+5].hist(muRel_E, bins=30,
                          density=True, histtype='step', lw=2, color='blue')
    fax[len(names)+5].set_ylabel('muRel_E')

    fax[len(names)+6].hist(muRel_N, bins=30,
                          density=True, histtype='step', lw=2, color='blue')
    fax[len(names)+6].set_ylabel('muRel_N')

    if posterior:
        res = fitter.load_mnest_results(remake_fits=False)
        for ii, name in enumerate(all_names):
            fax[ii].hist(res[name], weights=res['weights'], bins=30,
                        density=True, histtype='step', lw=2, color='red',
                        label='Posterior')
            fax[ii].yaxis.set_ticks([])
            fax[ii].yaxis.set_ticklabels([])
    fax[0].legend()

#    ax[len(names)].hist(piE, bins=30)
#    ax[len(names)].set_ylabel('piE')

    if posterior:
        plt.savefig(target + '_posteriors_priors.png')
    else:
        plt.savefig(target + '_priors.png')

def MultiDimPrior_Gen(tab, pnames, Nbins, binlims=None, plot_1D=False):
    """
    Takes multinest output and histograms them.
    Returns the histogram bins, along with an array of the
    PDF, CDF, and indices of the non-zero histogram bins.

    Note: the number of bins Nbins needs to be balanced
    against the number of parameters, otherwise python
    will crash. Because the number of bins goes as
    (number of bins)^(number of parameters), small changes to
    either number can cause the program to fail where it didn't
    before. For example, 40 bins with 5 parameters works, but
    40 bins with 7 parameters will fail.

    Parameters
    ----------
    tab : multinest results
        output of model_fitter.load_mnest_results

    pnames : list of strings
        names of the parameters

    Nbins : int
        number of bins

    binlims : list of length pnames, each entry is list of length 2
        Limits of  bins. The nth entry corresponds to
        [lower limit, upper limit] for the nth parameter in pnames.
        If None, limits will be inferred from multinest.

    Return
    ------
    bins : list of lists
        bin edges

    savearr : array of shape (N nonzero histogram bins, N parameters + 2)
        columns are PDF, CDF, parameter1 index, ..., parameterN index
        where indices are multi-dimensional bin indices

    """
    Nrows = len(tab)
    Ncols = len(pnames)
    print(Nrows)

    # Organize the multinest output into a table.
    # Shape of intable is (Nsamples, Nparams)
    intable = np.zeros((Nrows, Ncols))
    for ii, pname in enumerate(pnames):
        intable[:,ii] = tab[pname].data

    # Set bin-edges for each column in the table.
    # Bins is a list of length Nparams, and each entry in the list is an array of length Nbins
    bins = []
    if binlims is None:
        for i in range(Ncols):
            bins.append(np.linspace(intable[:,i].min(), intable[:,i].max(), Nbins))
    else:
        for i in range(Ncols):
            lo = binlims[i][0]
            hi = binlims[i][1]
            bins.append(np.linspace(lo, hi, Nbins))

    # Setup the multidimensional histogram
    # H shape is an array (Nbins-1, Nbins-1, ... , Nbins-1) where the length is the number of dimensions
    H, edges = np.histogramdd(intable, bins=bins, weights=tab['weights'])
    Hnorm = H/np.sum(H)
    print(np.sum(H))
    # Indices where histogram entries != 0
    indpdf = np.transpose(np.array((Hnorm).nonzero()))

    # Checking binning weirdness:
    # 1. Rounding for .nonzero() not problem. Using density=False, i.e.
    # H, edges = np.histogramdd(intable, bins=bins, density=False)
    # gives the same results for number of nonzero bins.
    # 2. Turning KDE off in dynesty plotting and comparing. Not sure
    # why, plots looked exactly the same...
    # 3. Turns out, you have to be very precise in selecting bin edges.
    # If they are too wide, you get the weird random gaps. IDK why...
    # Use plot_1D=True to visualize everything first and make sure
    # it is all okay.
    Nrows = len(indpdf[:,0])
    cdf = np.zeros(Nrows, dtype='float')
    pdf = np.zeros(Nrows, dtype='float')

    # Calculate 1D cdf and pdf using only non-zero bins.
    for j in range(Nrows):
        cdf[j] = cdf[j-1] + Hnorm[tuple(indpdf[j,:])]
        pdf[j] = Hnorm[tuple(indpdf[j,:])]

    # Stack PDF, CDF, and non-zero histogram indices into a single table.
    savearr = indpdf[:,0].astype(int)
    for i in range(Ncols-1):savearr = np.column_stack((savearr.astype(int), indpdf[:,i+1]))
    savearr = np.column_stack((pdf, cdf, savearr))

    # Plot 1-D histograms compared to binned to make sure it all works out.
    if plot_1D:
        for ii in np.arange(Ncols):
            plt.figure(ii+1)
            hist, binedges = np.histogram(savearr[:,ii+2], weights=savearr[:,0], bins=Nbins-1, density=True)
            f = ((Nbins-1)/Nbins)/(bins[ii][1] - bins[ii][0])
            plt.step(bins[ii][:-1], hist * f, where='post', lw=2, ls='-', label='N-D hist', zorder=1)
            plt.hist(intable[:, ii], weights=tab['weights'], bins=bins[ii], lw=2, ls=':',
                     label='1-D hist', histtype='step', density=True, zorder=2)
            plt.axvline(binlims[ii][0])
            plt.axvline(binlims[ii][1])
            plt.title(pnames[ii])
            plt.legend()
            plt.savefig(pnames[ii] + '.png')

    return bins, savearr

def onedee_posterior(tab, colname, log = False, nbin = 50):
    # Make a histogram of the mass using the weights. This creates the
    # marginalized 1D posteriors.

    fontsize1 = 18
    fontsize2 = 14

    xmax = tab[colname].max()
    xmin = tab[colname].min()

    plt.figure(1)
    plt.clf()
    plt.subplots_adjust(left=0.15)
    if log is True:
        xliml = 0.9 * xmin
        xlimu = 1.1 * xmax
        bins = np.logspace(np.log10(xliml), np.log10(xlimu), nbin)

    if log is False:
        xliml = xmin - 0.1 * np.abs(xmin)
        xlimu = xmax + 0.1 * np.abs(xmax)
        bins = np.linspace(xliml, xlimu, nbin)

    n, foo, patch = plt.hist(tab[colname], normed=True,
                             histtype='step', weights=tab['weights'],
                             bins=bins, linewidth=1.5)
    bin_centers = (bins[:-1] + bins[1:]) / 2
    plt.xlabel(colname)
    plt.xlabel(xtitle, fontsize=fontsize1, labelpad=10)
    plt.xlim(xliml, xlimu)
    plt.ylabel('Relative probability', fontsize=fontsize1, labelpad=10)
    plt.xticks(fontsize=fontsize2)
    plt.yticks(fontsize=fontsize2)

    if log is True:
        plt.xscale('log')

    ##########
    # Calculate 3-sigma boundaries for mass limits.
    ##########
    sig1_hi = 0.682689
    sig1_lo = 1.0 - sig1_hi
    sig_med = 0.5
    sig2_hi = 0.9545
    sig2_lo = 1.0 - sig2_hi
    sig3_hi = 0.9973
    sig3_lo = 1.0 - sig3_hi

    quantiles = [sig3_lo, sig2_lo, sig1_lo, sig_med, sig1_hi, sig2_hi, sig3_hi]

    mass_quants = model_fitter.weighted_quantile(tab[colname], quantiles,
                                                 sample_weight=tab['weights'])

    for qq in range(len(quantiles)):
        print( 'Mass at {0:.1f}% quantiles:  M = {1:5.2f}'.format(quantiles[qq]*100, mass_quants[qq]))

    ax = plt.axis()
    plt.axvline(mass_quants[0], color='red', linestyle='--')
#     plt.text(mass_quants[0] + 0.05, 0.9*ax[3],
#              colname + '>{0:5.2f} with 99.7% confidence'.format(mass_quants[0]), fontsize=12)
    plt.axvline(mass_quants[3], color='black', linestyle='--')
#     plt.text(mass_quants[3] + 0.05, 0.8*ax[3],
#              colname + 'median = {0:5.2f}'.format(mass_quants[0]), fontsize=12)
    plt.axvline(mass_quants[-1], color='green', linestyle='--')
#     plt.text(mass_quants[-1] + 0.05, 0.7*ax[3],
#              colname + '<{0:5.2f} with 99.7% confidence'.format(mass_quants[-1]), fontsize=12)
#
    plt.show()

    ##########
    # Save figure
    ##########

#    outfile =  outdir + outfile
#    fileUtil.mkdir(outdir)
#    print( 'writing plot to file ' + outfile)
#
#    plt.savefig(outfile)

    return

def onedee_posterior_smooth(tab, colname, log = False, nbin = 500):
    # Make a histogram of the mass using the weights. This creates the
    # marginalized 1D posteriors.

    fontsize1 = 18
    fontsize2 = 14

    xmax = tab[colname].max()
    xmin = tab[colname].min()

    plt.figure(1)
    plt.clf()
    plt.subplots_adjust(left=0.23)
    if log is True:
        xliml = 0.9 * xmin
        xlimu = 1.1 * xmax
        bins = np.logspace(np.log10(xliml), np.log10(xlimu), nbin)

    if log is False:
        xliml = xmin - 0.1 * np.abs(xmin)
        xlimu = xmax + 0.1 * np.abs(xmax)
        bins = np.linspace(xliml, xlimu, nbin)

    n, b = np.histogram(tab[colname], bins = bins,
                        weights = tab['weights'], normed = True)

    n = norm_kde(n, 10.)
    b0 = 0.5 * (b[1:] + b[:-1])

    n, b, _ = plt.hist(b0, bins = b, weights = n)

    plt.xlabel(colname)
    xtitle = colname
    plt.xlabel(xtitle, fontsize=fontsize1, labelpad=10)
    plt.xlim(xliml, xlimu)
    plt.ylabel('Relative probability', fontsize=fontsize1, labelpad=10)
    plt.xticks(fontsize=fontsize2)
    plt.yticks(fontsize=fontsize2)

    if log is True:
        plt.xscale('log')

    ##########
    # Calculate 3-sigma boundaries for mass limits.
    ##########
    sig1_hi = 0.682689
    sig1_lo = 1.0 - sig1_hi
    sig_med = 0.5
    sig2_hi = 0.9545
    sig2_lo = 1.0 - sig2_hi
    sig3_hi = 0.9973
    sig3_lo = 1.0 - sig3_hi

    quantiles = [sig3_lo, sig2_lo, sig1_lo, sig_med, sig1_hi, sig2_hi, sig3_hi]

    mass_quants = model_fitter.weighted_quantile(tab[colname], quantiles,
                                                 sample_weight=tab['weights'])

    for qq in range(len(quantiles)):
        print( 'Mass at {0:.1f}% quantiles:  M = {1:5.2f}'.format(quantiles[qq]*100, mass_quants[qq]))

    ax = plt.axis()
    # plot median and +/- 3 sigma
    plt.axvline(mass_quants[0], color='k', linestyle='--')
    plt.axvline(mass_quants[3], color='k', linestyle='--')
    plt.axvline(mass_quants[-1], color='k', linestyle='--')
    plt.show()

    ##########
    # Save figure
    ##########

#    outfile =  outdir + outfile
#    fileUtil.mkdir(outdir)
#    print( 'writing plot to file ' + outfile)
#
#    plt.savefig(outfile)

    return

def twodee_posterior_points(tab, colname1, colname2, logx = False, logy = False):
    """
    Color of points = weight
    """
    fig = plt.figure(1, figsize=(8,6))
    plt.subplots_adjust(left = 0.15)
    plt.clf()
    plt.scatter(tab[colname1], tab[colname2], c = tab['weights'],
                alpha = 0.3, s = 1, cmap = 'viridis')
    plt.xlabel(colname1)
    plt.ylabel(colname2)
    plt.colorbar()

    if logx == True:
        plt.xscale('log')

    if logy == True:
        plt.yscale('log')

    plt.show()

def twodee_posterior_smooth(tab, colname1, colname2, logx = False, logy = False,
                            xnbin = 500, ynbin = 500,
                            xsval = 2, ysval = 2,
                            span=None, weights=None, levels=None,
                            color='gray', plot_datapoints=False, plot_density=True,
                            plot_contours=True, no_fill_contours=False, fill_contours=True,
                            contour_kwargs=None, contourf_kwargs=None, data_kwargs=None,
                            **kwargs):
    try:
        str_type = types.StringTypes
        float_type = types.FloatType
        int_type = types.IntType
    except:
        str_type = str
        float_type = float
        int_type = int

    """
    Basically the _hist2d function from dynesty, but with a few mods I made.
    https://github.com/joshspeagle/dynesty/blob/master/dynesty/plotting.py

    Parameters
    ----------
    tab : astropy table
        Read in output from load_mnest_results_to_tab

    colname1, 2 : str
        The stuff you wanna plot. Column names from the tab

    x, ysval : int (float?)
        smoothing value. Kinda inverse of xnbin and ynbin I think

    span : iterable with shape (ndim,), optional
        A list where each element is either a length-2 tuple containing
        lower and upper bounds or a float from `(0., 1.]` giving the
        fraction of (weighted) samples to include. If a fraction is provided,
        the bounds are chosen to be equal-tailed. An example would be::

            span = [(0., 10.), 0.95, (5., 6.)]

        Default is `0.999999426697` (5-sigma credible interval).

    weights : iterable with shape (nsamps,)
        Weights associated with the samples. Default is `None` (no weights).

    levels : iterable, optional
        The contour levels to draw. Default are `[0.5, 1, 1.5, 2]`-sigma.

    ax : `~matplotlib.axes.Axes`, optional
        An `~matplotlib.axes.axes` instance on which to add the 2-D histogram.
        If not provided, a figure will be generated.

    color : str, optional
        The `~matplotlib`-style color used to draw lines and color cells
        and contours. Default is `'gray'`.

    plot_datapoints : bool, optional
        Whether to plot the individual data points. Default is `False`.

    plot_density : bool, optional
        Whether to draw the density colormap. Default is `True`.

    plot_contours : bool, optional
        Whether to draw the contours. Default is `True`.

    no_fill_contours : bool, optional
        Whether to add absolutely no filling to the contours. This differs
        from `fill_contours=False`, which still adds a white fill at the
        densest points. Default is `False`.

    fill_contours : bool, optional
        Whether to fill the contours. Default is `True`.

    contour_kwargs : dict
        Any additional keyword arguments to pass to the `contour` method.

    contourf_kwargs : dict
        Any additional keyword arguments to pass to the `contourf` method.

    data_kwargs : dict
        Any additional keyword arguments to pass to the `plot` method when
        adding the individual data points.

    """

    x = tab[colname1]
    y = tab[colname2]
    weights = tab['weights']

    fig = plt.figure(1)
    plt.clf()
    ax = plt.gca()

    # Determine plotting bounds.
    data = [x, y]
    if span is None:
        span = [0.999999426697 for i in range(2)]
    span = list(span)
    if len(span) != 2:
        raise ValueError("Dimension mismatch between samples and span.")
    for i, _ in enumerate(span):
        try:
            xmin, xmax = span[i]
        except:
            q = [0.5 - 0.5 * span[i], 0.5 + 0.5 * span[i]]
            span[i] = dyutil.quantile(data[i], q, weights=weights)

    # The default "sigma" contour levels.
    if levels is None:
        levels = 1.0 - np.exp(-0.5 * np.arange(0.5, 2.1, 0.5) ** 2)

    # Color map for the density plot, over-plotted to indicate the
    # density of the points near the center.
    density_cmap = LinearSegmentedColormap.from_list(
        "density_cmap", [color, (1, 1, 1, 0)])

    # Color map used to hide the points at the high density areas.
    white_cmap = LinearSegmentedColormap.from_list(
        "white_cmap", [(1, 1, 1), (1, 1, 1)], N=2)

    # This "color map" is the list of colors for the contour levels if the
    # contours are filled.
    rgba_color = colorConverter.to_rgba(color)
    contour_cmap = [list(rgba_color) for l in levels] + [rgba_color]
    for i, l in enumerate(levels):
        contour_cmap[i][-1] *= float(i) / (len(levels)+1)

    xmax = tab[colname1].max()
    xmin = tab[colname1].min()

    ymax = tab[colname2].max()
    ymin = tab[colname2].min()

    if logx is True:
        xliml = 0.9 * xmin
        xlimu = 1.1 * xmax
        xbins = np.logspace(np.log10(xliml), np.log10(xlimu), xnbin)

    if logx is False:
        xliml = xmin - 0.1 * np.abs(xmin)
        xlimu = xmax + 0.1 * np.abs(xmax)
        xbins = np.linspace(xliml, xlimu, xnbin)

    if logy is True:
        yliml = 0.9 * ymin
        ylimu = 1.1 * ymax
        ybins = np.logspace(np.log10(yliml), np.log10(ylimu), ynbin)

    if logy is False:
        yliml = ymin - 0.1 * np.abs(ymin)
        ylimu = ymax + 0.1 * np.abs(ymax)
        ybins = np.linspace(yliml, ylimu, ynbin)

    # We'll make the 2D histogram to directly estimate the density.
    H, X, Y = np.histogram2d(x.flatten(), y.flatten(), bins=[xbins, ybins],
                             range=list(map(np.sort, span)),
                             weights=weights)
    # Smooth the results.
    H = norm_kde(H, [xsval, ysval])

    # Compute the density levels.
    Hflat = H.flatten()
    inds = np.argsort(Hflat)[::-1]
    Hflat = Hflat[inds]
    sm = np.cumsum(Hflat)
    sm /= sm[-1]
    V = np.empty(len(levels))
    for i, v0 in enumerate(levels):
        try:
            V[i] = Hflat[sm <= v0][-1]
        except:
            V[i] = Hflat[0]
    V.sort()
    m = (np.diff(V) == 0)
    if np.any(m) and plot_contours:
        logging.warning("Too few points to create valid contours.")
        logging.warning("Make xnbin or ynbin bigger!!!!!!!")
    while np.any(m):
        V[np.where(m)[0][0]] *= 1.0 - 1e-4
        m = (np.diff(V) == 0)
    V.sort()

    # Compute the bin centers.
    X1, Y1 = 0.5 * (X[1:] + X[:-1]), 0.5 * (Y[1:] + Y[:-1])

    # Extend the array for the sake of the contours at the plot edges.
    H2 = H.min() + np.zeros((H.shape[0] + 4, H.shape[1] + 4))
    H2[2:-2, 2:-2] = H
    H2[2:-2, 1] = H[:, 0]
    H2[2:-2, -2] = H[:, -1]
    H2[1, 2:-2] = H[0]
    H2[-2, 2:-2] = H[-1]
    H2[1, 1] = H[0, 0]
    H2[1, -2] = H[0, -1]
    H2[-2, 1] = H[-1, 0]
    H2[-2, -2] = H[-1, -1]
    X2 = np.concatenate([X1[0] + np.array([-2, -1]) * np.diff(X1[:2]), X1,
                         X1[-1] + np.array([1, 2]) * np.diff(X1[-2:])])
    Y2 = np.concatenate([Y1[0] + np.array([-2, -1]) * np.diff(Y1[:2]), Y1,
                         Y1[-1] + np.array([1, 2]) * np.diff(Y1[-2:])])

    # Plot the data points.
    if plot_datapoints:
        if data_kwargs is None:
            data_kwargs = dict()
        data_kwargs["color"] = data_kwargs.get("color", color)
        data_kwargs["ms"] = data_kwargs.get("ms", 2.0)
        data_kwargs["mec"] = data_kwargs.get("mec", "none")
        data_kwargs["alpha"] = data_kwargs.get("alpha", 0.1)
        ax.plot(x, y, "o", zorder=-1, rasterized=True, **data_kwargs)

    # Plot the base fill to hide the densest data points.
    if (plot_contours or plot_density) and not no_fill_contours:
        ax.contourf(X2, Y2, H2.T, [V.min(), H.max()],
                    cmap=white_cmap, antialiased=False)

    if plot_contours and fill_contours:
        if contourf_kwargs is None:
            contourf_kwargs = dict()
        contourf_kwargs["colors"] = contourf_kwargs.get("colors", contour_cmap)
        contourf_kwargs["antialiased"] = contourf_kwargs.get("antialiased",
                                                             False)
        ax.contourf(X2, Y2, H2.T, np.concatenate([[0], V, [H.max()*(1+1e-4)]]),
                    **contourf_kwargs)

    # Plot the density map. This can't be plotted at the same time as the
    # contour fills.
    elif plot_density:
        ax.pcolor(X, Y, H.max() - H.T, cmap=density_cmap)

    # Plot the contour edge colors.
    if plot_contours:
        if contour_kwargs is None:
            contour_kwargs = dict()
        contour_kwargs["colors"] = contour_kwargs.get("colors", color)
        ax.contour(X2, Y2, H2.T, V, **contour_kwargs)

    ax.set_xlabel(colname1)
    ax.set_ylabel(colname2)

    ax.set_xlim(span[0])
    ax.set_ylim(span[1])

    if logx is True:
        ax.set_xscale('log')

    if logy is True:
        ax.set_yscale('log')

    plt.show()

    return

def separate_mode_files(mnest_dir = './', mnest_root = 'aa_'):
    # FIXME: Is there a more intelligent way to deal with all the indices???
    # But it seems to work for now...
    mode_file = mnest_dir + mnest_root + 'post_separate.dat'

    empty_lines = []
    with open(mode_file, 'r') as orig_file:
        for num, line in enumerate(orig_file, start=0):
            if line == '\n':
                empty_lines.append(num)

    # Error checking
    if len(empty_lines) % 2 != 0:
        print('SOMETHING BAD HAPPENED!')

    idx_range = int(len(empty_lines)/2)

    orig_tab = np.loadtxt(mode_file)
    for idx in np.arange(idx_range):
        start_idx = empty_lines[idx * 2 + 1] + 1 - 2 * (idx + 1)
        if idx != np.arange(idx_range)[-1]:
            end_idx = empty_lines[idx * 2 + 2] - 2 * (idx + 1)
            np.savetxt(mnest_dir + mnest_root + 'mode' + str(idx) + '.dat', orig_tab[start_idx:end_idx])
        else:
            np.savetxt(mnest_dir + mnest_root + 'mode' + str(idx) + '.dat', orig_tab[start_idx:])

    return

def convert_gp_names(base):
    """
    Fix pre-April 1 2021 gp parameter names to match new names (commit de12879)
    Read in old fits file, and rename columns
    """
    # Gather all the .fits files that need to have columns renamed
    fits_files = glob.glob(base + '*.fits')

    fits_files = [x for x in fits_files if 'old' not in x]

    # Loop through and fix.
    for ff in fits_files:
        tab = Table.read(ff)
        for colname in tab.colnames:
            if 'omegaofour_So' in colname:
                tab.rename_column(colname, colname.replace('omegaofour_So', 'omega04_S0'))
            elif 'So' in colname:
                tab.rename_column(colname, colname.replace('So', 'S0'))
            elif 'omegaofour' in colname:
                tab.rename_column(colname, colname.replace('omegaofour', 'omega04'))
            elif 'omegao' in colname:
                tab.rename_column(colname, colname.replace('omegao', 'omega0'))

        # Move the original file by renaming it "old".
        os.system('mv ' + ff + ' ' + ff[:-5] + '_old.fits')

        try:
            tab.write(ff, overwrite=False)
        except:
            print('Exists already, skipping')

def convert_pre_2020apr_u0_sign(old_base, new_base):
    """
    Read in Multinest output from a pre 2020 April u0 sign
    convention and convert everything into a new convention
    (same as Gould+ 2004).
    """
    # Fetch the model class out of the YAML file.
    yaml_file = open(old_base + 'params.yaml', 'r')
    yaml_params = yaml.full_load(yaml_file)
    with open(new_base + 'params.yaml', 'w') as f:
        foo = yaml.dump(yaml_params, f)

    # Get the parameter names off the model class.
    model_class = getattr(model, yaml_params['model'])
    params = model_class.all_param_names

    # Find the u0_amp or beta parameters that need to be changed.
    # There could be multiples (some fit paras, others not).
    u0_cols = []
    for pp in range(len(params)):
        if ((params[pp] == 'u0_amp') or
            (params[pp] == 'beta')):
            u0_cols.append(pp)

    # # Fix the <base>.txt file
    # tab = Table.read(old_base + '.txt', format='ascii')
    # for cc in u0_cols:
    #     tab['col{0:d}'.format(cc+1+2)] *= -1.0
    # tab.write(new_base + '.txt', format='ascii.no_header')

    # # Fix the <base>summary.txt file
    # tab_a = Table.read(old_base + 'summary.txt', format='ascii')
    # for cc in u0_cols:
    #     mean_col = 0 * len(params) + 1 + cc
    #     std_col = 1 * len(params) + 1 + cc
    #     maxL_col = 2 * len(params) + 1 + cc
    #     maxAP_col = 3 * len(params) + 1 + cc

    #     tab_a['col{0:d}'.format(mean_col)] *= -1.0
    #     tab_a['col{0:d}'.format(maxL_col)] *= -1.0
    #     tab_a['col{0:d}'.format(maxAP_col)] *= -1.0
    # tab_a.write(new_base + 'summary.txt', format='ascii.no_header')

    # # Fix the *ev.dat file.
    # tab_e = Table.read(old_base + 'ev.dat', format='ascii')
    # for cc in u0_cols:
    #     tab_e['col{0:d}'.format(cc+1)] *= -1.0
    # tab_e.write(new_base + 'ev.dat', format='ascii.no_header')

    # # Fix the *phys_live.points file
    # tab_l = Table.read(old_base + 'phys_live.points', format='ascii')
    # for cc in u0_cols:
    #     tab_l['col{0:d}'.format(cc+1)] *= -1.0
    # tab_l.write(new_base + 'phys_live.points', format='ascii.no_header')

    # # Fix *post_equal_weights.dat file.
    # tab_pe = Table.read(old_base + 'post_equal_weights.dat', format='ascii')
    # for cc in u0_cols:
    #     tab_pe['col{0:d}'.format(cc+1)] *= -1.0
    # tab_pe.write(new_base + 'post_equal_weights.dat', format='ascii.no_header')

    # Fix *post_separate.dat
    _psep_old = open(old_base + 'post_separate.dat', 'r')
    _psep_new = open(new_base + 'post_separate.dat', 'w')
    for line in _psep_old:
        sline = line.strip()
        if sline == '':
            _psep_new.write(line)
        else:
            # Found a sample row.
            fields = sline.split()

            # Loop through u0 cols
            for cc in u0_cols:
                if fields[cc+2].startswith('-'):
                    fields[cc+2] = fields[cc+2][1:]
                else:
                    fields[cc+2] = '-' + fields[cc+2]

            _psep_new.write('  '.join(fields) + '\n')

    _psep_old.close()
    _psep_new.close()

    # # Fix *stats.dat
    # _psta_old = open(old_base + 'stats.dat', 'r')
    # _psta_new = open(new_base + 'stats.dat', 'w')
    # for line in _psta_old:
    #     sline = line.strip()

    #     for cc in u0_cols:
    #         # Looking for a match like '3 ' or '13 ':
    #         if sline.startswith(str(cc + 1) + ' '):
    #             fields = sline.split()
    #             u0_val_str = fields[1]

    #             # Determine new "sign" string
    #             if u0_val_str.startswith('-'):
    #                 new_sign_str = ' '
    #                 u0_abs_str = u0_val_str[1:]
    #             else:
    #                 new_sign_str = '-'
    #                 u0_abs_str = u0_val_str

    #             # Replace the sign character.
    #             ldx = line.find(u0_abs_str)
    #             line = line[:ldx-1] + new_sign_str + line[ldx:]

    #     _psta_new.write(line)

    # _psta_old.close()
    # _psta_new.close()

    # # Fix *resume.dat
    # shutil.copyfile(old_base + 'resume.dat', new_base + 'resume.dat')

    # Fix *live.points
    # PASS... I don't know what is in this file.

    return

reason: import numpy,from scipy,import astropy,from astropy


text:
# USAGE
# python /home/nmorales/cxgn/DroneImageScripts/ImageProcess/CalculatePhenotypeSurf.py --image_paths /folder/mypic1.png,/folder/mypic2.png --outfile_paths /export/mychoppedimages/outimage2.png,/export/mychoppedimages/outimage2.png

# import the necessary packages
import argparse
import imutils
import cv2
import numpy as np
import urllib.request
import math
from matplotlib import pyplot as plt

# construct the argument parse and parse the arguments
ap = argparse.ArgumentParser()
ap.add_argument("-i", "--image_paths", required=True, help="image path")
ap.add_argument("-o", "--outfile_paths", required=True, help="file path directory where the output will be saved")
ap.add_argument("-r", "--results_outfile_path", required=True, help="file path where results will be saved")
args = vars(ap.parse_args())

input_images = args["image_paths"]
outfile_paths = args["outfile_paths"]
results_ourfile = args["results_outfile_path"]
images = input_images.split(",")
outfiles = outfile_paths.split(",")

count = 0
for image in images:
    img = cv2.imread(image)
    img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)

    orb = cv2.xfeatures2d.SURF_create()
    kp, des = orb.detectAndCompute(img, None)
    print(kp)
    #print(des)

    kpsimage = cv2.drawKeypoints(img, kp, img, flags=cv2.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS)
    #cv2.imshow('image'+str(count),kpsimage)
    cv2.imwrite(outfiles[count], kpsimage)

    count += 1

#cv2.waitKey(0)

max_stars_repo_name: soybase/DroneImageScripts

max_issues_repo_name: soybase/DroneImageScripts

max_forks_repo_name: soybase/DroneImageScripts